// 27/11/2012 by Graeme Watt . // Updated at NNLO for plots to appear in Annual Review [arXiv:1301.6754]. // ROOT script to demonstrate statistical combination of MSTW, CTEQ // and NNPDF predictions as an alternative to the PDF4LHC "envelope". // Use Z, W+/W-, t-tbar total cross-section predictions for LHC // at 7 TeV centre-of-mass energy, and H(126GeV) for LHC at 8 TeV. // Default values of alpha_s with no alpha_s uncertainties included. void combinationAR(void) { gROOT->Reset(); gROOT->SetStyle("Plain"); gStyle->SetPadTickY(1); gStyle->SetEndErrorSize(5); random = new TRandom3(); // random number generator with default seed // Options to symmetrise (or not) uncertainties and random predictions. bool lSymmetrise=true; // symmetric (then average equals best-fit) //bool lSymmetrise=false; // asymmetric (then average differs from best-fit) // Loop over all four cross sections: Z, W+/W-, t-tbar, H(126GeV). for (int ixs=1;ixs<=4;ixs++) { //int ixs=1; // Z //int ixs=2; // W+/W- //int ixs=3; // t-tbar //int ixs=4; // H(126GeV) // Specify predictions for best-fit and eigenvector PDF sets (MSTW/CTEQ). // Use LHC total cross sections at NNLO for Annual Review. const int neigenM=20, neigenC=25, Nrep=100; double xsM[2*neigenM+1], xsC[2*neigenC+1], xsN[Nrep+1]; if (ixs==1) { // sigma(Z) * BR(Z->ll) in nb double xsM[] = {0.9586,0.956163,0.961024,0.956452,0.960093,0.958733,0.958508,0.958655,0.958413,0.960343,0.957342,0.957201,0.959695,0.962894,0.955995,0.958598,0.958444,0.973046,0.946165,0.960344,0.956398,0.960049,0.957846,0.959131,0.958198,0.95316,0.962591,0.959257,0.95791,0.957686,0.955528,0.959098,0.958171,0.959202,0.957153,0.959837,0.956499,0.958706,0.958178,0.958594,0.958662}; // MSTW08 NNLO (68% C.L.) best-fit PDF set and 2*20 eigenvector PDF sets double xsC[] = {0.965772,0.95202,0.977193,0.958822,0.973713,0.978612,0.952095,0.972743,0.960638,0.965932,0.965545,0.944523,0.982221,0.971782,0.959831,0.969963,0.963359,0.961813,0.971673,0.960614,0.970711,0.965549,0.964668,0.969732,0.951739,0.970662,0.959644,0.964061,0.967957,0.967073,0.963767,0.966914,0.96491,0.967547,0.965095,0.962639,0.960675,0.973624,0.961488,0.967709,0.958263,0.97379,0.976469,0.963737,0.967786,0.961289,0.966009,0.959498,0.965872,0.966662,0.962505}; // CT10 NNLO (90% C.L.) best-fit PDF set and 2*25 eigenvector PDF sets double xsN[] = {0.958297,0.957505,0.9407,0.958234,0.949411,0.953493,0.97465,0.95362,0.946743,0.953314,0.954094,0.970197,0.949628,0.93952,0.953774,0.940197,0.937347,0.9578,0.960256,0.948021,0.979274,0.976545,0.955994,0.955746,0.987039,0.958066,0.961704,0.965811,0.950259,0.955397,0.940136,0.957104,0.954646,0.956744,0.975113,0.96163,0.962208,0.967413,0.955562,0.971205,0.954094,0.957304,0.949812,0.948602,0.958745,0.945081,0.961967,0.9489,0.983306,0.94165,0.971562,0.946307,0.956059,0.958731,0.97471,0.929019,0.950305,0.946006,0.964624,0.973059,0.952807,0.95198,0.966037,0.957391,0.966249,0.963228,0.963539,0.951701,0.965922,0.961833,0.934362,0.975526,0.970751,0.96671,0.990717,0.960102,0.964582,0.963189,0.949974,0.969126,0.963086,0.956987,0.95165,0.975204,0.965551,0.950032,0.971696,0.948336,0.952926,0.956911,0.976996,0.953294,0.963967,0.956386,0.948783,0.971335,0.947516,0.945066,0.95356,0.962484,0.957102}; // NNPDF2.3 NNLO (68% C.L.) average PDF set and 100 replica PDF sets } if (ixs==2) { // sigma(W+)/sigma(W-) cross-section ratio double xsM[] = {1.42879,1.42839,1.42951,1.4345,1.42506,1.42794,1.42991,1.43143,1.42562,1.4303,1.42822,1.42469,1.43246,1.42602,1.43085,1.42665,1.43524,1.4292,1.42851,1.42889,1.4298,1.42892,1.42921,1.42679,1.43025,1.42891,1.43013,1.43351,1.42703,1.4285,1.43034,1.43246,1.42823,1.42834,1.42997,1.4358,1.42321,1.43148,1.42939,1.43108,1.42479}; // MSTW08 NNLO (68% C.L.) best-fit PDF set and 2*20 eigenvector PDF sets double xsC[] = {1.46155,1.46674,1.45765,1.45971,1.46433,1.45798,1.46588,1.46374,1.4606,1.46211,1.46167,1.46409,1.46079,1.45329,1.46967,1.45757,1.4655,1.46574,1.45677,1.45517,1.4673,1.46113,1.4631,1.46136,1.46349,1.45703,1.46513,1.46176,1.46203,1.46465,1.45825,1.45773,1.46479,1.46644,1.45795,1.45681,1.46585,1.468,1.45749,1.46702,1.45262,1.44557,1.47078,1.4666,1.45698,1.46557,1.46023,1.46423,1.46082,1.46244,1.46521}; // CT10 NNLO (90% C.L.) best-fit PDF set and 2*25 eigenvector PDF sets double xsN[] = {1.45101,1.45884,1.45606,1.47196,1.43315,1.45529,1.47144,1.45493,1.44171,1.46634,1.43796,1.4194,1.45065,1.45185,1.42044,1.45843,1.45774,1.45114,1.47274,1.44222,1.43906,1.44588,1.43201,1.46264,1.46005,1.4558,1.43693,1.44625,1.44388,1.44788,1.46417,1.44214,1.44025,1.44927,1.46135,1.47908,1.44926,1.45022,1.45176,1.44942,1.44494,1.44975,1.42755,1.45911,1.46549,1.45161,1.45549,1.46296,1.45536,1.48455,1.44471,1.44444,1.4538,1.4501,1.46177,1.45985,1.43231,1.47255,1.44245,1.4604,1.43845,1.44989,1.45004,1.45591,1.46688,1.43674,1.44308,1.45875,1.44576,1.435,1.4489,1.45118,1.45411,1.44171,1.4488,1.42533,1.48041,1.459,1.46667,1.44664,1.44652,1.4553,1.44967,1.44012,1.46558,1.44697,1.44891,1.45358,1.45963,1.44836,1.43005,1.44161,1.46518,1.42992,1.45164,1.484,1.45577,1.46802,1.44452,1.44391,1.44352}; // NNPDF2.3 NNLO (68% C.L.) average PDF set and 100 replica PDF sets } if (ixs==3) { // sigma(t-tbar) in pb for m_t = 173.18 GeV double xsM[] = {155.005,155.67,154.343,154.747,155.196,155.029,154.983,154.994,155.018,155.102,154.938,155.072,154.951,154.482,155.319,155.087,154.787,151.705,157.998,154.38,155.784,154.433,155.256,154.207,155.386,157.099,153.214,154.592,155.613,154.513,155.767,154.526,155.179,154.746,155.513,154.301,155.748,154.49,155.578,154.893,155.572}; // MSTW08 NNLO (68% C.L.) best-fit PDF set and 2*20 eigenvector PDF sets double xsC[] = {155.154,155.75,154.664,155.587,154.659,154.305,156.062,154.246,155.801,155.068,155.215,162.328,149.647,153.647,156.602,157.339,153.149,155.886,153.979,155.817,154.496,154.875,155.683,152.315,160.705,151.081,157.974,156.082,153.726,151.729,159.49,154.342,155.695,155.778,154.501,156.992,156.115,153.143,156.32,154.334,157.84,155.793,153.682,155.733,154.639,156.548,154.864,156.392,154.605,154.663,153.666}; // CT10 NNLO (90% C.L.) best-fit PDF set and 2*25 eigenvector PDF sets double xsN[] = {157.421,154.65,159.991,154.318,153.803,154.36,153.95,156.961,160.725,154.861,156.7,154.651,153.744,155.728,154.107,161.764,169.722,158.422,159.066,164.248,151.253,154.648,159.602,162.677,144.704,161.164,158.104,155.659,162.607,162.361,159.566,156.555,159.037,159.37,149.218,161.486,155.296,158.477,158.209,157.96,164.681,157.086,163.346,163.216,152.204,159.916,153.98,157.15,150.745,162.996,154.753,156.094,151.334,149.521,156.795,155.899,159.051,157.213,157.855,153.09,159.206,162.169,156.73,155.081,157.088,155.409,151.203,158.299,155.101,159.521,168.741,154.064,153.027,153.279,148.277,154.677,155.972,157.557,163.037,160.442,155.995,157.982,160.746,157.757,152.175,159.1,156.032,159.681,165.702,160.679,153.881,160.572,155.66,155.32,165.038,153.237,159.437,162.117,159.096,158.76,157.889}; // NNPDF2.3 NNLO (68% C.L.) average PDF set and 100 replica PDF sets } if (ixs==4) { // sigma(gg->H) in pb for M_H = 126 GeV double xsM[] = {19.9534,19.9485,19.9315,19.862,20.0123,19.9497,19.9172,19.9365,19.9459,19.8992,19.9736,19.9071,19.9446,19.73,20.0647,20.0054,19.7947,19.8809,19.9814,19.9514,19.942,19.9472,19.9443,19.9542,19.9233,20.0479,19.8555,19.9195,20.0178,19.8885,19.914,19.9444,19.9401,19.949,19.9448,19.9098,19.9434,19.8513,19.9289,19.9421,19.912}; // MSTW08 NNLO (68% C.L.) best-fit PDF set and 2*20 eigenvector PDF sets double xsC[] = {19.8592,19.8422,19.8498,19.8763,19.8254,19.4771,20.2206,19.747,19.9155,19.8046,19.8762,19.8618,19.7798,19.7544,19.9303,19.9088,19.7887,19.9783,19.6654,19.8721,19.8294,19.8038,19.932,20.0436,19.431,20.0163,19.7194,19.8161,19.9029,20.0021,19.64,19.8505,19.8461,19.9201,19.783,19.8993,19.6327,19.785,19.8723,19.8483,19.7453,19.7737,19.9367,19.8389,19.8557,19.7803,19.8484,19.9439,19.754,19.9555,19.7122}; // CT10 NNLO (90% C.L.) best-fit PDF set and 2*25 eigenvector PDF sets double xsN[] = {21.1866,21.263,20.6108,20.5776,21.0038,20.8755,20.9238,21.5741,21.2487,20.8147,21.2531,21.5387,21.0741,21.0117,21.3194,20.4308,21.1834,20.8977,21.2197,21.4249,21.0964,20.7703,21.441,21.2594,20.4547,21.4152,21.6142,21.365,21.4786,21.1767,20.8769,21.3449,21.288,21.0598,21.2587,21.3279,21.0074,21.8011,21.1899,21.1354,21.3771,21.6492,21.3852,21.2945,20.9858,21.1507,21.1612,21.1875,21.1546,21.6365,20.9622,21.1719,20.984,21.3474,21.2064,20.6459,21.4007,21.5164,21.3695,21.2517,21.1852,21.0628,21.3177,20.8672,21.2707,21.0252,21.0498,21.435,21.0149,21.3909,21.2588,21.3729,20.9784,20.5647,21.0142,20.9029,20.8349,20.9327,21.1774,21.0006,21.1695,21.424,21.064,21.1269,21.0848,21.5118,21.1517,21.2378,21.3925,21.2915,21.3414,21.199,21.0812,21.4985,21.5595,20.5622,20.9718,21.4239,21.1835,21.6834,21.1473}; // NNPDF2.3 NNLO (68% C.L.) average PDF set and 100 replica PDF sets } // Calculate usual PDF uncertainties for the three groups. double *uncM = getPDFuncertainty(2*neigenM,xsM,false,false); // MSTW08 double *uncC = getPDFuncertainty(2*neigenC,xsC,false,true); // CT10 double *uncN = getPDFuncertainty(Nrep,xsN,true,false); // NNPDF2.3 // Generate random predictions for MSTW08 and CT10 separately. const int Npdf = Nrep; // choose equal to number of NNPDF replicas double *xsMrand = getRandomPredictions(neigenM,xsM,random,Npdf,lSymmetrise,false); double *xsCrand = getRandomPredictions(neigenC,xsC,random,Npdf,lSymmetrise,true); // Calculate average and s.d. for MSTW08 and CT10 separately. double *uncMrand = getPDFuncertainty(Npdf,xsMrand,true,false); // MSTW08 double *uncCrand = getPDFuncertainty(Npdf,xsCrand,true,false); // CT10 // Compare best-fit and Hessian uncertainty with average // and s.d. over random predictions for MSTW08 and CT10. // If lSymmetrise=true, then should be equal for large Npdf. // Differences for small Npdf are due to statistical fluctuation. cout << "ixs = " << ixs << ", MSTWhess: " << uncM[0] << "+/-" << uncM[3] << " (+" << uncM[1] << ", -" << uncM[2] << ")" << endl; cout << "ixs = " << ixs << ", MSTWrand: " << uncMrand[0] << "+/-" << uncMrand[3] << endl; cout << "ixs = " << ixs << ", CT10hess: " << uncC[0] << "+/-" << uncC[3] << " (+" << uncC[1] << ", -" << uncC[2] << ")" << endl; cout << "ixs = " << ixs << ", CT10rand: " << uncCrand[0] << "+/-" << uncCrand[3] << endl; cout << "ixs = " << ixs << ", NNPDF2.3: " << uncN[0] << "+/-" << uncN[3] << endl; // Calculate PDF4LHC envelope and midpoint for the three groups. if (lSymmetrise) { double upper = TMath::Max(TMath::Max(uncM[0]+uncM[3],uncC[0]+uncC[3]),uncN[0]+uncN[3]); double lower = TMath::Min(TMath::Min(uncM[0]-uncM[3],uncC[0]-uncC[3]),uncN[0]-uncN[3]); } else { double upper = TMath::Max(TMath::Max(uncM[0]+uncM[1],uncC[0]+uncC[1]),uncN[0]+uncN[1]); double lower = TMath::Min(TMath::Min(uncM[0]-uncM[2],uncC[0]-uncC[2]),uncN[0]-uncN[2]); } double midpoint = 0.5*(upper+lower); double envelope = 0.5*(upper-lower); cout << "ixs = " << ixs << ", Envelope: " << midpoint << "+/-" << envelope << endl; // Fill array with equal number of random predictions from three groups. double *xsCombined = new double[Npdf+Npdf+Nrep+1]; for (int ipdf=1;ipdf<=Npdf;ipdf++) { xsCombined[ipdf] = xsMrand[ipdf]; // MSTW08 } for (int ipdf=1;ipdf<=Npdf;ipdf++) { xsCombined[Npdf+ipdf] = xsCrand[ipdf]; // CT10 } for (int irep=1;irep<=Nrep;irep++) { xsCombined[Npdf+Npdf+irep] = xsN[irep]; // NNPDF2.3 } // Calculate average and s.d. for combination of MSTW08, CT10 and NNPDF2.3. double *uncCombined = getPDFuncertainty(Npdf+Npdf+Nrep,xsCombined,true,false); cout << "ixs = " << ixs << ", Combined: " << uncCombined[0] << "+/-" << uncCombined[3] << endl; // Make a plot using the numbers printed out above. // Show MSTW08, CT10, NNPDF2.3 as markers with error bars. // Open markers = Hessian, closed markers = random. // Horizontal lines for envelope and statistical combination. // Solid lines = envelope, dashed lines = combination. grM = new TGraphAsymmErrors(1); grMrand = new TGraphErrors(1); grC = new TGraphAsymmErrors(1); grCrand = new TGraphErrors(1); grN = new TGraphErrors(1); // MSTW08: grM->SetPoint(0,0.9,uncM[0]); grMrand->SetPoint(0,1.1,uncMrand[0]); if (lSymmetrise) grM->SetPointError(0,0,0,uncM[3],uncM[3]); else grM->SetPointError(0,0,0,uncM[2],uncM[1]); grMrand->SetPointError(0,0,uncMrand[3]); grM->SetMarkerStyle(24); grMrand->SetMarkerStyle(20); grM->SetMarkerSize(2); grMrand->SetMarkerSize(2); grM->SetLineWidth(2); grMrand->SetLineWidth(2); // CT10: grC->SetPoint(0,1.9,uncC[0]); grCrand->SetPoint(0,2.1,uncCrand[0]); if (lSymmetrise) grC->SetPointError(0,0,0,uncC[3],uncC[3]); else grC->SetPointError(0,0,0,uncC[2],uncC[1]); grCrand->SetPointError(0,0,uncCrand[3]); grC->SetMarkerStyle(25); grCrand->SetMarkerStyle(21); grC->SetMarkerColor(kBlue); grCrand->SetMarkerColor(kBlue); grC->SetLineColor(kBlue); grCrand->SetLineColor(kBlue); grC->SetMarkerSize(2); grCrand->SetMarkerSize(2); grC->SetLineWidth(2); grCrand->SetLineWidth(2); // NNPDF2.3: grN->SetPoint(0,3,uncN[0]); grN->SetPointError(0,0,uncN[3]); grN->SetMarkerStyle(22); grN->SetMarkerColor(kRed); grN->SetLineColor(kRed); grN->SetMarkerSize(2); grN->SetLineWidth(2); // Set range of vertical axis. double ymin, ymax; if (ixs==1) { ymin = 0.93; ymax = 1.00; } if (ixs==2) { ymin = 1.40; ymax = 1.50; } if (ixs==3) { ymin = 142; ymax = 166; } if (ixs==4) { ymin = 18.5; ymax = 22.; } // Draw canvas and axes. canvas = new TCanvas(); canvas->SetLeftMargin(0.15); canvas->SetBottomMargin(0.15); h = new TH2D("h","h",50,0.0,4.0,50,ymin,ymax); h->SetStats(kFALSE); h->Draw("axis"); h->SetTitle(""); if (ixs==1) h->SetYTitle("#sigma_{Z^{0}} #upoint B(Z^{0} #rightarrow l^{+}l^{-}) (nb)"); if (ixs==2) h->SetYTitle("R_{#pm} #equiv #sigma_{W^{+}} / #sigma_{W^{-}}"); if (ixs==3) h->SetYTitle("#sigma_{t#bar{t}} (pb)"); if (ixs==4) h->SetYTitle("#sigma_{H} (pb)"); h->GetYaxis()->SetTitleSize(1.5*h->GetYaxis()->GetTitleSize()); h->GetXaxis()->SetNdivisions(0); // Draw markers for the three groups. grM->Draw("P"); grMrand->Draw("P"); grC->Draw("P"); grCrand->Draw("P"); grN->Draw("P"); // Draw text labels. text = new TLatex(); text->SetTextAlign(22); text->SetNDC(1); if (ixs==1) text->DrawLatex(0.5,0.95,"NNLO Z^{0} #rightarrow l^{+}l^{-} at the LHC (#sqrt{s} = 7 TeV)"); if (ixs==2) text->DrawLatex(0.5,0.95,"NNLO W^{+}/W^{-} ratio at the LHC (#sqrt{s} = 7 TeV)"); if (ixs==3) text->DrawLatex(0.5,0.96,"NNLO+NNLL t#bar{t} cross sections at the LHC (#sqrt{s} = 7 TeV)"); if (ixs==4) text->DrawLatex(0.5,0.95,"NNLO gg#rightarrowH at the LHC (#sqrt{s} = 8 TeV) for M_{H} = 126 GeV"); text->SetNDC(0); text->DrawLatex(1,ymin+0.1*(ymax-ymin),"MSTW08"); text->SetTextSize(text->GetTextSize()/1.5); text->DrawLatex(1,ymin+0.04*(ymax-ymin),"#alpha_{S}(M_{Z}^{2}) = 0.1171"); text->SetTextSize(text->GetTextSize()*1.5); text->SetTextColor(kBlue); text->DrawLatex(2,ymin+0.1*(ymax-ymin),"CT10"); text->SetTextSize(text->GetTextSize()/1.5); text->DrawLatex(2,ymin+0.04*(ymax-ymin),"#alpha_{S}(M_{Z}^{2}) = 0.1180"); text->SetTextSize(text->GetTextSize()*1.5); text->SetTextColor(kRed); text->DrawLatex(3,ymin+0.1*(ymax-ymin),"NNPDF2.3"); text->SetTextSize(text->GetTextSize()/1.5); text->DrawLatex(3,ymin+0.04*(ymax-ymin),"#alpha_{S}(M_{Z}^{2}) = 0.1190"); text->SetTextSize(text->GetTextSize()*1.5); text->SetTextColor(kBlack); text->SetTextSize(text->GetTextSize()/1.2); text->SetTextAlign(22); text->DrawLatex(2,ymin-0.1*(ymax-ymin),"#splitline{Open markers: usual best-fit and 68% C.L. Hessian uncertainty.}{Closed markers: average and s.d. over random predictions.}"); text->SetTextColor(kGreen+2); text->DrawLatex(2,ymin+0.94*(ymax-ymin),"Solid lines: envelope and midpoint."); text->SetTextColor(kMagenta); text->DrawLatex(2,ymin+0.88*(ymax-ymin),"Dashed lines: statistical combination."); text->SetTextColor(kBlack); text->SetTextSize(text->GetTextSize()*1.2); text->SetNDC(1); text->SetTextAlign(22); text->SetTextSize(text->GetTextSize()/1.5); if (ixs==3) text->DrawLatex(0.5,0.3,"Top++ (v1.4), #mu_{R} = #mu_{F} = m_{t}^{pole} = 173.18 GeV"); if (ixs==4) text->DrawLatex(0.5,0.3,"ggh@nnlo (v1.4.1), #mu_{R} = #mu_{F} = M_{H} / 2"); text->SetTextSize(text->GetTextSize()*1.5); text->SetTextAlign(); text->SetTextAlign(22); text->SetTextSize(text->GetTextSize()/2); text->SetTextAngle(90); text->DrawLatex(0.97,0.5,"G. Watt (November 2012)"); text->SetTextAngle(0); text->SetTextAlign(); text->SetTextSize(text->GetTextSize()*2); // Draw horizontal lines for envelope and midpoint. lMidpoint = new TLine(0,midpoint,4,midpoint); lMidpoint->SetLineColor(kGreen+2); lMidpoint->SetLineWidth(1); lMidpoint->Draw(); lUpperEnvelope = new TLine(0,upper,4,upper); lUpperEnvelope->SetLineColor(kGreen+2); lUpperEnvelope->SetLineWidth(1); lUpperEnvelope->Draw(); lLowerEnvelope = new TLine(0,lower,4,lower); lLowerEnvelope->SetLineColor(kGreen+2); lLowerEnvelope->SetLineWidth(1); lLowerEnvelope->Draw(); // Draw horizontal lines for statistical combination. lAverage = new TLine(0,uncCombined[0],4,uncCombined[0]); lAverage->SetLineStyle(2); lAverage->SetLineWidth(2); lAverage->SetLineColor(kMagenta); lAverage->Draw(); lUpperCombined = new TLine(0,uncCombined[0]+uncCombined[3],4,uncCombined[0]+uncCombined[3]); lUpperCombined->SetLineStyle(2); lUpperCombined->SetLineWidth(2); lUpperCombined->SetLineColor(kMagenta); lUpperCombined->Draw(); lLowerCombined = new TLine(0,uncCombined[0]-uncCombined[3],4,uncCombined[0]-uncCombined[3]); lLowerCombined->SetLineStyle(2); lLowerCombined->SetLineWidth(2); lLowerCombined->SetLineColor(kMagenta); lLowerCombined->Draw(); // Print plot to file. string outfile="combination_"; if (ixs==1) outfile+="z0"; if (ixs==2) outfile+="wpwmratio"; if (ixs==3) outfile+="ttbar"; if (ixs==4) outfile+="ggh126"; outfile+=".pdf"; // or .eps or .png etc. canvas->Print(outfile.c_str()); } // end loop over ixs } // Function to calculate either Monte Carlo or Hessian PDF uncertainties. // Pass an array xs[imem] (imem=0,...,Nmem). l90cl to correct for 90% C.L. double *getPDFuncertainty(int nmem, double *xs, bool lMonteCarlo, bool l90cl) { double *ret = new double[4]; double central=0.0,errsym=0.0,errplus=0.0,errminus=0.0; if (lMonteCarlo) { // Calculate the average and standard deviation // using Eqs. (2.13) and (2.14) of arXiv:1205.4024v2. double xsav = 0.0, xssd = 0.0; for (int imem=1;imem<=nmem;imem++) { xsav += xs[imem]; xssd += xs[imem]**2; } xsav /= nmem; xssd /= nmem; xssd = nmem/(nmem-1.0)*(xssd-xsav**2); if (xssd>0.0) xssd = sqrt(xssd); else xssd = 0.0; central = xsav; errplus = errminus = errsym = xssd; } else { // Calculate the symmetric and asymmetric Hessian uncertainties // using Eqs. (2.7), (2.8) and (2.9) of arXiv:1205.4024v2. central = xs[0]; for (int ieigen=1;ieigen<=nmem/2;ieigen++) { errplus += (TMath::Max(TMath::Max(xs[2*ieigen-1]-xs[0],xs[2*ieigen]-xs[0]),0.0))**2; errminus += (TMath::Max(TMath::Max(xs[0]-xs[2*ieigen-1],xs[0]-xs[2*ieigen]),0.0))**2; errsym += (xs[2*ieigen-1]-xs[2*ieigen])**2; } errsym = 0.5*sqrt(errsym); errplus = sqrt(errplus); errminus = sqrt(errminus); if (l90cl) { // correct from 90% C.L. to 68% C.L. (for CT10) errsym /= 1.64485; errplus /= 1.64485; errminus /= 1.64485; } } // Return central value and asymmetric/symmetric PDF uncertainties. ret[0] = central; ret[1] = errplus; ret[2] = errminus; ret[3] = errsym; return ret; } // Function to generate Npdf random predictions using // either Eq. (6.4) or Eq. (6.5) of arXiv:1205.4024v2. // Pass an array xs[imem] (imem=0,...,2*neigen). double *getRandomPredictions(int neigen, double *xs, TRandom3 *random, int Npdf, bool lSymmetrise, bool l90cl) { double *xsrand = new double[Npdf+1]; double scale=1.0; if (l90cl) scale=1.64485; // correct from 90% C.L. to 68% C.L. (for CT10) for (int ipdf=1;ipdf<=Npdf;ipdf++) { xsrand[ipdf] = xs[0]; // Loop over number of eigenvectors. for (int ieigen=1;ieigen<=neigen;ieigen++) { double r = random->Gaus(0.0,1.0); // Gaussian random number if (lSymmetrise) xsrand[ipdf] += 0.5*r*TMath::Abs(xs[2*ieigen-1]-xs[2*ieigen])/scale; else { // not symmetrised if (r<0.0) xsrand[ipdf] -= r*(xs[2*ieigen]-xs[0])/scale; // negative direction else xsrand[ipdf] += r*(xs[2*ieigen-1]-xs[0])/scale; // positive direction } } } // // Simpler calculation of random predictions if only considering // // one observable in isolation (as in the above example). Then can // // simply choose random predictions from a Gaussian distribution // // with mean xs[0] and standard deviation equal to Hessian uncertainty. // // But if considering more than one observable, it is necessary to use // // the eigenvectors as done above to preserve the correlations. // double *unc = getPDFuncertainty(2*neigen,xs,false,l90cl); // for (int ipdf=1;ipdf<=Npdf;ipdf++) { // xsrand[ipdf] = xs[0]; // double r = random->Gaus(0.0,1.0); // Gaussian random number // if (lSymmetrise) xsrand[ipdf] += r*unc[3]; // else { // not symmetrised // if (r<0.0) xsrand[ipdf] += r*unc[2]; // negative direction // else xsrand[ipdf] += r*unc[1]; // positive direction // } // } return xsrand; }