// 07/08/2012 by Graeme Watt . // 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 and H(120GeV) total cross-section // predictions for LHC at 7 TeV centre-of-mass energy. // Default values of alpha_s with no alpha_s uncertainties included. void combination(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(120GeV). 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(120GeV) // Specify predictions for best-fit and eigenvector PDF sets (MSTW/CTEQ). // Use 7 TeV LHC total cross sections at NLO with settings given in: // G. Watt, JHEP 1109 (2011) 069 [arXiv:1106.5788 [hep-ph]]. const int neigenM=20, neigenC=26, 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.930756,0.929768,0.931549,0.933702,0.927479,0.930614,0.931,0.931299,0.930364,0.931465,0.929849,0.934523,0.928419,0.931329,0.929959,0.92964,0.931202,0.931188,0.926651,0.932011,0.9295,0.920232,0.948637,0.935878,0.927351,0.932776,0.929763,0.933286,0.928693,0.933065,0.925094,0.931464,0.930514,0.929681,0.931358,0.929209,0.931297,0.929037,0.929783,0.930843,0.930678}; // MSTW08 NLO (68% C.L.) best-fit PDF set and 2*20 eigenvector PDF sets double xsC[] = {0.928505,0.916708,0.939455,0.923139,0.933992,0.91964,0.936276,0.909667,0.945215,0.924481,0.9336,0.925838,0.929904,0.916821,0.94398,0.929996,0.926529,0.932349,0.926864,0.934785,0.926091,0.927511,0.928529,0.932065,0.917169,0.929681,0.924882,0.92851,0.926958,0.928047,0.928378,0.928188,0.929731,0.923866,0.926958,0.922374,0.931885,0.929916,0.928337,0.934955,0.924643,0.929927,0.928966,0.927379,0.930253,0.930438,0.934398,0.928594,0.916449,0.925317,0.9222,0.926754,0.922719}; // CT10 NLO (90% C.L.) best-fit PDF set and 2*26 eigenvector PDF sets double xsN[] = {0.919927,0.946335,0.92309,0.936325,0.917945,0.922983,0.906422,0.925994,0.923316,0.91854,0.916318,0.87893,0.91341,0.937759,0.917391,0.919634,0.857901,0.892099,0.911519,0.904914,0.908949,0.94506,0.915149,0.925105,0.927887,0.907808,0.92323,0.940284,0.929494,0.922581,0.919157,0.910698,0.930242,0.923061,0.930356,0.911946,0.92524,0.849002,0.932164,0.94948,0.872119,0.90896,0.926945,0.937368,0.920167,0.91952,0.900806,0.931916,0.910141,0.939411,0.936523,0.947663,0.897334,0.914154,0.9304,0.940281,0.940333,0.924865,0.926969,0.911667,0.926796,0.919871,0.922628,0.927066,0.918152,0.928042,1.01088,0.912117,0.909852,0.914898,0.931579,0.919797,0.909069,0.920774,0.9285,0.920135,0.909436,0.931235,0.911116,0.914808,0.931542,0.928558,0.919938,0.915897,0.908006,0.906398,0.902192,0.932876,0.905347,0.925007,0.905557,0.923553,0.924193,0.93734,0.916363,0.906656,0.943977,0.912012,0.914994,0.903633,0.924454}; // NNPDF2.1 NLO (68% C.L.) average PDF set and 100 replica PDF sets } if (ixs==2) { // sigma(W+)/sigma(W-) cross-section ratio double xsM[] = {1.42181,1.42031,1.42315,1.41888,1.42525,1.42072,1.42363,1.421,1.42242,1.42436,1.41815,1.41783,1.42435,1.42537,1.41712,1.42845,1.41938,1.42098,1.42299,1.42285,1.42119,1.42145,1.42227,1.42161,1.42207,1.42011,1.42259,1.42048,1.42833,1.42489,1.41888,1.42513,1.42052,1.42446,1.42046,1.41714,1.42501,1.42019,1.42596,1.41819,1.42387}; // MSTW08 NLO (68% C.L.) best-fit PDF set and 2*20 eigenvector PDF sets double xsC[] = {1.46389,1.46825,1.46002,1.46069,1.4674,1.46858,1.46021,1.46432,1.46447,1.4613,1.46726,1.46114,1.4654,1.46718,1.45866,1.44793,1.47859,1.45356,1.46946,1.46578,1.46286,1.45925,1.46754,1.45934,1.46597,1.46947,1.45297,1.46213,1.4656,1.46236,1.46495,1.46149,1.46519,1.46761,1.46185,1.45328,1.46792,1.47081,1.45877,1.47256,1.44992,1.45928,1.46734,1.47092,1.45586,1.45841,1.47167,1.46014,1.46494,1.46434,1.46611,1.46575,1.46872}; // CT10 NLO (90% C.L.) best-fit PDF set and 2*26 eigenvector PDF sets double xsN[] = {1.46467,1.48551,1.45113,1.47065,1.44906,1.47228,1.45313,1.41277,1.45912,1.40199,1.46021,1.48301,1.47499,1.49261,1.4667,1.46892,1.47172,1.4741,1.4603,1.48205,1.44872,1.4792,1.46394,1.45994,1.47463,1.4911,1.4554,1.46579,1.47118,1.47808,1.4654,1.4716,1.46528,1.47226,1.48864,1.46611,1.43971,1.44907,1.49131,1.46162,1.4585,1.45649,1.43899,1.4785,1.48978,1.49986,1.45036,1.46977,1.45711,1.46257,1.43594,1.50782,1.46493,1.48181,1.45826,1.39701,1.42294,1.45091,1.44807,1.45359,1.48381,1.45414,1.4457,1.4697,1.43001,1.47082,1.48535,1.47103,1.4678,1.48313,1.45286,1.46057,1.47031,1.46427,1.45031,1.49005,1.48065,1.47711,1.48929,1.47649,1.44233,1.47442,1.44191,1.47011,1.48711,1.44855,1.45728,1.47973,1.46847,1.47766,1.44908,1.48013,1.47482,1.45605,1.46519,1.46082,1.45865,1.47367,1.46891,1.47999,1.46351}; // NNPDF2.1 NLO (68% C.L.) average PDF set and 100 replica PDF sets } if (ixs==3) { // sigma(t-tbar) in pb for m_t = 171.3 GeV double xsM[] = {167.923056,168.551562,167.406749,168.033103,167.80091,167.845431,168.04314,168.021827,167.858295,167.818941,168.071355,167.305616,168.309152,167.99163,167.84255,167.783292,167.975365,166.141893,170.098048,168.121997,167.719809,170.711342,163.401085,165.997897,169.18785,166.837479,168.455445,167.22972,168.597499,166.794823,169.745934,167.511321,168.088551,168.172107,167.77762,168.762133,167.470761,168.233623,168.3793,168.212157,167.847682}; // MSTW08 NLO (68% C.L.) best-fit PDF set and 2*20 eigenvector PDF sets double xsC[] = {157.805147,158.38331,157.27522,158.451481,157.14953,157.635498,157.906654,163.903496,152.341654,158.295274,157.188016,157.919071,157.750843,160.521159,153.908261,157.123006,158.480565,157.853368,157.619993,153.879257,159.623011,155.083547,160.328445,153.568409,163.822901,157.080241,159.339565,157.067794,158.488159,158.834921,157.060628,155.699465,159.964312,161.709838,156.767237,158.782153,157.058952,157.685869,157.929456,155.698981,159.789017,159.540199,156.522153,158.670536,156.969648,156.938415,158.931058,156.708219,161.704013,157.613114,157.842466,158.224794,155.208717}; // CT10 NLO (90% C.L.) best-fit PDF set and 2*26 eigenvector PDF sets double xsN[] = {170.138295,161.64164,167.385076,165.947941,173.343179,173.289344,175.001462,168.946216,165.44022,166.712552,173.3325,159.271666,167.700912,162.109561,174.487195,169.168534,170.900326,182.382001,174.273009,174.201859,170.752398,170.244037,173.290546,176.068216,168.465762,181.069794,168.327752,156.859844,164.647954,165.87623,167.636497,171.422921,166.191682,164.28855,168.450492,173.033775,170.274312,167.175368,166.606023,159.737856,175.246068,174.411359,171.938776,164.543621,169.762302,177.083757,180.469708,170.403454,171.429707,166.176921,165.631746,173.161979,168.047875,169.301462,166.312281,155.257302,159.760887,169.290965,162.691854,174.812711,170.002867,166.668565,165.893517,169.487995,171.924214,170.987962,175.318333,177.044793,174.484275,180.245011,169.704083,172.479221,173.21057,171.651001,169.993846,169.840021,175.152344,175.023789,181.078822,168.275713,159.792776,177.212653,168.79597,178.41011,174.485761,173.829007,171.431706,160.022807,181.995618,168.632564,179.261869,170.194041,169.422545,164.554784,170.048485,175.100337,168.355411,155.282379,172.160614,171.790989,172.292281}; // NNPDF2.1 NLO (68% C.L.) average PDF set and 100 replica PDF sets } if (ixs==4) { // sigma(gg->H) in pb for M_H = 120 GeV double xsM[] = {12.413,12.4427,12.3883,12.4275,12.3979,12.4266,12.3953,12.4242,12.4058,12.3916,12.4435,12.2898,12.4898,12.4461,12.3684,12.3703,12.429,12.3935,12.4068,12.4024,12.4231,12.4587,12.3313,12.39,12.4279,12.4058,12.4166,12.4041,12.4348,12.3445,12.5372,12.405,12.4167,12.4319,12.4023,12.4082,12.4147,12.4659,12.3046,12.4011,12.4167}; // MSTW08 NLO (68% C.L.) best-fit PDF set and 2*20 eigenvector PDF sets double xsC[] = {11.6705,11.6778,11.6639,11.6941,11.6453,11.8671,11.498,11.8073,11.5312,11.7143,11.6125,11.6711,11.6696,11.6142,11.7665,11.6265,11.7099,11.4996,11.7685,11.5968,11.7065,11.723,11.6377,11.807,11.5335,11.66,11.6911,11.6868,11.6574,11.6323,11.6997,11.7077,11.6383,11.7202,11.6103,11.7034,11.7125,11.675,11.6638,11.6692,11.7124,11.7156,11.6289,11.692,11.6473,11.6216,11.7395,11.7338,11.5556,11.6688,11.6338,11.7423,11.5096}; // CT10 NLO (90% C.L.) best-fit PDF set and 2*26 eigenvector PDF sets double xsN[] = {12.3539,12.0462,12.4863,12.4906,12.4829,12.5601,12.4711,12.6781,12.2124,12.4622,12.3309,11.6069,12.4823,12.0067,12.4828,12.4951,12.5903,12.0661,12.254,12.1867,12.1854,12.2323,12.3774,12.1636,12.3597,12.4908,12.4333,11.9027,12.353,12.4479,12.5804,12.5277,12.4354,12.4118,12.1013,12.4656,12.2614,12.3781,12.472,12.0925,12.612,12.4706,12.4592,12.4691,12.2744,11.8719,12.5018,12.2253,12.2037,11.9433,12.5576,12.5745,12.5311,12.5003,12.3384,12.3075,12.2471,12.432,12.2149,12.6705,12.4419,12.5751,12.455,11.867,12.3541,12.257,12.5383,12.5499,12.5487,12.2214,12.4559,12.4744,12.1947,12.1733,12.3296,12.4697,12.3214,12.3521,12.322,12.4444,12.33,12.3998,12.3121,12.501,12.3105,12.4338,12.4131,12.3075,12.2528,12.3014,12.5679,12.4549,12.4173,12.1788,12.3625,12.4711,12.3649,12.0919,12.352,12.587,12.1928}; // NNPDF2.1 NLO (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.1 // 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.1: " << 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.1 } // Calculate average and s.d. for combination of MSTW08, CT10 and NNPDF2.1. 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.1 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.1: 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.89; ymax = 0.96; } if (ixs==2) { ymin = 1.38; ymax = 1.52; } if (ixs==3) { ymin = 145; ymax = 185; } if (ixs==4) { ymin = 11.; ymax = 13.; } // 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,"NLO Z^{0} #rightarrow l^{+}l^{-} at the LHC (#sqrt{s} = 7 TeV)"); if (ixs==2) text->DrawLatex(0.5,0.95,"NLO W^{+}/W^{-} ratio at the LHC (#sqrt{s} = 7 TeV)"); if (ixs==3) text->DrawLatex(0.5,0.96,"NLO t#bar{t} cross sections at the LHC (#sqrt{s} = 7 TeV)"); if (ixs==4) text->DrawLatex(0.5,0.95,"NLO gg#rightarrowH at the LHC (#sqrt{s} = 7 TeV) for M_{H} = 120 GeV"); text->SetNDC(0); text->DrawLatex(1,ymin+0.1*(ymax-ymin),"MSTW08"); text->SetTextColor(kBlue); text->DrawLatex(2,ymin+0.1*(ymax-ymin),"CT10"); text->SetTextColor(kRed); text->DrawLatex(3,ymin+0.1*(ymax-ymin),"NNPDF2.1"); 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); // 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+="ggh120"; outfile+=".png"; // or .eps or .pdf 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; }