// 29/06/2012 by Graeme Watt . // ROOT script to reproduce Figures 13 and 14 of arXiv:1205.4024v2. // Plot PDF uncertainties versus number of random predictions. // Use Z, W+/W-, t-tbar and H(120GeV) total cross-section // predictions for LHC at 7 TeV centre-of-mass energy. // Symmetrised (as in paper) or asymmetric options. { gROOT->Reset(); gROOT->SetStyle("Plain"); gStyle->SetPadTickX(1); gStyle->SetPadTickY(1); random = new TRandom3(); // random number generator with default seed // Options to symmetrise (or not) uncertainties and random predictions. bool lSymmetrise=true; // symmetric (as used in Figures 13 and 14 of paper) //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) // Loop over different/same random numbers. // iSame = 1 => Figure 13 of paper. // iSame = 0 => Figure 14 of paper. for (int iSame=0;iSame<=1;iSame++) { bool lSameRand; if (iSame) lSameRand = true; // same random numbers for different Npdf values (Figure 13) else lSameRand = false; // different random numbers for different Npdf values (Figure 14) // Specify predictions for best-fit and eigenvector PDF sets. // 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 neigen=20; double xs[2*neigen+1]; if (ixs==1) { // sigma(Z) * BR(Z->ll) in nb double xs[] = {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 } if (ixs==2) { // sigma(W+)/sigma(W-) cross-section ratio double xs[] = {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 } if (ixs==3) { // sigma(t-tbar) in pb for m_t = 171.3 GeV double xs[] = {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 } if (ixs==4) { // sigma(gg->H) in pb for M_H = 120 GeV double xs[] = {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 } // Calculate the symmetric and asymmetric Hessian uncertainties // using Eqs. (2.7), (2.8) and (2.9) of arXiv:1205.4024v2. double errsym = 0.0, errplus = 0.0, errminus = 0.0; for (int ieigen=1;ieigen<=neigen;ieigen++) { errsym += (xs[2*ieigen-1]-xs[2*ieigen])**2; 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 = 0.5*sqrt(errsym); errplus = sqrt(errplus); errminus = sqrt(errminus); cout << "ixs = " << ixs << ", Hessian: " << xs[0] << "+/-" << errsym << " (+" << errplus << ", -" << errminus << ")" << endl; // Specify minimum and maximum values of Npdf for random predictions. const int NpdfMin = 2; const int NpdfMax = 1000; grHess = new TGraphAsymmErrors(NpdfMax-NpdfMin+1); grHessP = new TGraph(NpdfMax-NpdfMin+1); grHessM = new TGraph(NpdfMax-NpdfMin+1); grRand = new TGraphErrors(NpdfMax-NpdfMin+1); grRandP = new TGraph(NpdfMax-NpdfMin+1); grRandM = new TGraph(NpdfMax-NpdfMin+1); // Choose random numbers from a Gaussian distribution // with mean zero and variance one. double r; // single random number double *rand = new double[neigen*NpdfMax+1]; // array of random numbers for (int ir=0;ir<=neigen*NpdfMax;ir++) { rand[ir] = random->Gaus(0.0,1.0); } // Generate the Npdf random predictions where Npdf = NpdfMin, ..., NpdfMax. for (int Npdf=NpdfMin;Npdf<=NpdfMax;Npdf++) { // Construct Npdf random predictions using either // Eq. (6.4) or Eq. (6.5) of arXiv:1205.4024v2. double *xsrand = new double[Npdf+1]; for (int ipdf=1;ipdf<=Npdf;ipdf++) { xsrand[ipdf] = xs[0]; // Loop over number of eigenvectors. for (int ieigen=1;ieigen<=neigen;ieigen++) { // Same/different random numbers for different Npdf values. if (lSameRand) r = rand[(ipdf-1)*neigen+ieigen]; // same for different Npdf else r = random->Gaus(0.0,1.0); // different for different Npdf if (lSymmetrise) xsrand[ipdf] += 0.5*r*TMath::Abs(xs[2*ieigen-1]-xs[2*ieigen]); else { // not symmetrised if (r<0.0) xsrand[ipdf] -= r*(xs[2*ieigen]-xs[0]); // negative direction else xsrand[ipdf] += r*(xs[2*ieigen-1]-xs[0]); // positive direction } } //cout << ipdf << " " << xsrand[ipdf] << endl; } // Calculate the average and standard deviation // over the Npdf random predictions. // Use Eqs. (2.13) or (2.14) of arXiv:1205.4024v2 with Nrep->Npdf. double xsav = 0.0, xssd = 0.0; for (int ipdf=1;ipdf<=Npdf;ipdf++) { xsav += xsrand[ipdf]; xssd += xsrand[ipdf]**2; } xsav /= Npdf; xssd /= Npdf; xssd = Npdf/(Npdf-1.0)*(xssd-xsav**2); if (xssd>0.0) xssd = sqrt(xssd); else xssd = 0.0; // Fill the TGraphs with everything normalised to the best-fit prediction. grHess->SetPoint(Npdf-NpdfMin,Npdf,xs[0]/xs[0]); if (lSymmetrise) { grHess->SetPointError(Npdf-NpdfMin,0.0,0.0,errsym/xs[0],errsym/xs[0]); grHessP->SetPoint(Npdf-NpdfMin,Npdf,(xs[0]+errsym)/xs[0]); grHessM->SetPoint(Npdf-NpdfMin,Npdf,(xs[0]-errsym)/xs[0]); } else { grHess->SetPointError(Npdf-NpdfMin,0.0,0.0,errminus/xs[0],errplus/xs[0]); grHessP->SetPoint(Npdf-NpdfMin,Npdf,(xs[0]+errplus)/xs[0]); grHessM->SetPoint(Npdf-NpdfMin,Npdf,(xs[0]-errminus)/xs[0]); } grRand->SetPoint(Npdf-NpdfMin,Npdf,xsav/xs[0]); grRand->SetPointError(Npdf-NpdfMin,Npdf,xssd/xs[0]); grRandP->SetPoint(Npdf-NpdfMin,Npdf,(xsav+xssd)/xs[0]); grRandM->SetPoint(Npdf-NpdfMin,Npdf,(xsav-xssd)/xs[0]); } // end of loop over Npdf // Print out average and standard deviation for Npdf = NpdfMax. cout << "ixs = " << ixs << ", Random: " << xsav << "+/-" << xssd << endl; // Specify maximum percentage uncertainty used for range of y-axis. double percentMax; if (ixs==1) percentMax=6; // Z^0 if (ixs==2) percentMax=3; // W^+/W^- if (ixs==3) percentMax=10; // t-tbar if (ixs==4) percentMax=5; // gg->H(120GeV) // Make a plot versus Npdf normalised to best-fit of average // and standard deviation compared to Hessian uncertainty. canvas = new TCanvas(); canvas->SetLogx(); canvas->SetLeftMargin(0.15); canvas->SetBottomMargin(0.15); h = new TH2D("h","h",50,NpdfMin,NpdfMax,50,1.0-percentMax/100.0,1.0+percentMax/100.0); h->SetStats(kFALSE); h->SetTitle(";Number of random predictions, N_{pdf};Ratio to best-fit prediction"); h->GetXaxis()->SetTitleOffset(1.3); h->GetYaxis()->SetTitleOffset(1.3); h->GetXaxis()->SetTitleSize(h->GetXaxis()->GetTitleSize()*1.3); h->GetYaxis()->SetTitleSize(h->GetYaxis()->GetTitleSize()*1.3); h->GetXaxis()->SetLabelSize(h->GetXaxis()->GetLabelSize()*1.3); h->GetYaxis()->SetLabelSize(h->GetYaxis()->GetLabelSize()*1.3); h->Draw("axis"); // Draw the average and s.d. of the random predictions. grRand->Draw("3"); grRand->Draw("LX"); grRandP->Draw("L"); grRandM->Draw("L"); grRand->SetFillStyle(1001); grRand->SetLineStyle(1); grRand->SetLineColor(kBlack); grRand->SetFillColor(kYellow); grRand->SetLineWidth(2); grRandP->SetLineColor(kBlack); grRandM->SetLineColor(kBlack); grRandP->SetLineStyle(1); grRandM->SetLineStyle(1); // Draw the Hessian uncertainty band. grHess->Draw("3"); grHess->Draw("LX"); grHessP->Draw("L"); grHessM->Draw("L"); grHess->SetFillColor(kBlue); grHess->SetFillStyle(3395); grHess->SetLineStyle(2); grHess->SetLineWidth(2); grHess->SetLineColor(kBlue); grHessP->SetLineColor(kBlue); grHessM->SetLineColor(kBlue); grHessP->SetLineStyle(2); grHessM->SetLineStyle(2); // Draw the legend. leg = new TLegend(0.25,0.7,0.85,0.85,""); leg->SetFillColor(0); leg->AddEntry(grHess,"MSTW 2008 NLO PDFs (68% C.L.)","LF"); leg->AddEntry(grRand,"Average and s.d. over N_{pdf} predictions","LF"); leg->Draw(); // Draw text labels. text = new TLatex(); text->SetNDC(); text->SetTextAlign(22); if (ixs==1) text->DrawLatex(0.5,0.95,"NLO Z^{0} cross section at the LHC (#sqrt{s} = 7 TeV)"); if (ixs==2) text->DrawLatex(0.5,0.95,"NLO W^{+}/W^{-} cross-section ratio at the LHC (#sqrt{s} = 7 TeV)"); if (ixs==3) text->DrawLatex(0.5,0.96,"NLO t#bar{t} cross section 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->SetTextSize(text->GetTextSize()/1.15); if (lSameRand) text->DrawLatex(0.52,0.23,"Same random numbers for each value of N_{pdf}"); else text->DrawLatex(0.52,0.22,"Different random numbers for each value of N_{pdf}"); text->SetTextSize(text->GetTextSize()*1.15); canvas->RedrawAxis(); // Print plot to file. string outfile="conv_"; if (lSameRand) outfile+="samerand_"; else outfile+="diffrand_"; 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 iSame } // end loop over ixs }