analysis/HERWIG_UE/makeSpec.C

0001 // makeSpec.C
0002 //
0003 // Draws the HERWIG UE cross-checks against STAR and against the PPG10
0004 // sPHENIX result. PYTHIA has been dropped; instead the two particle sources
0005 // written by HerwigSTARUECheck are overlaid:
0006 //   "_hepmc" -- generator-level HepMC record
0007 //   "_g4"    -- PHG4TruthInfoContainer sPHENIX primaries
0008 // Both live in the same per-sample output file.
0009 //
0010 // The sumET panel additionally shows the PPG10 unfolded data (points + band)
0011 // and the PPG10 HERWIG truth curve.
0012
0013 const int NSAMP = 7;
0014 const std::string samp[NSAMP] = {"MB", "Jet5", "Jet12", "Jet20", "Jet30", "Jet40", "Jet50"};
0015 const double HerwigCS[NSAMP] = {3.1909e+10, 1.8437e+08, 1.132355e+06, 5.2613e+04, 2.0694e+03, 1.0510e+02, 5.2089};
0016
0017 bool csWeightsAlreadyInFiles = true;
0018
0019 // ---- the two particle sources written by HerwigSTARUECheck -------------
0020 const int NVERS = 2;
0021 const std::string versTag[NVERS] = {"hepmc", "g4"};
0022 const std::string versLabel[NVERS] = {"HERWIG (HepMC)", "HERWIG (sPHENIX primary)"};
0023 const int versColor[NVERS] = {kBlue, kOrange + 7};
0024 const int versMarker[NVERS] = {33, 34};
0025 // h_nevents bins in the new module: 1 = processed, 2 = HepMC node found,
0026 // 3 = G4TruthInfo node found. Normalize each source by its own bin.
0027 const int versNevtBin[NVERS] = {2, 3};
0028
0029 const std::string inDir = "/sphenix/tg/tg01/jets/bkimelman/Herwig_UE";
0030 const std::string prodTag = "Aug5";  // production tag in the input file names
0031
0032 namespace
0033 {
0034   void styleObj(TH1 *h, int color, int marker)
0035   {
0036     if (!h) return;
0037     h->SetMarkerStyle(marker);
0038     h->SetMarkerColor(color);
0039     h->SetLineColor(color);
0040     h->SetMarkerSize(1);
0041   }
0042 }  // namespace
0043
0044 void makeSpec()
0045 {
0046   gStyle->SetOptStat(0);
0047
0048   // ------------------------------------------------------------------
0049   // PPG10 reference: unfolded sPHENIX data + HERWIG truth
0050   // ------------------------------------------------------------------
0051   TFile *fPPG = new TFile("/sphenix/u/bkimelman/PPG10_Prelim/plot_result_hist_output_efrac_bkg_cut_run28_iter_3_1000toys.root", "READ");
0052
0053   TH1D *hPPG = (TH1D *) fPPG->Get("unfold_hist_calib_dijet_reweight_trim_10_2_etEffCorrected_2");
0054   hPPG->SetMarkerStyle(20);
0055   hPPG->SetMarkerColor(kBlack);
0056   hPPG->SetLineColor(kBlack);
0057   hPPG->SetLineWidth(2);
0058   hPPG->SetMarkerSize(1);
0059
0060   TGraphAsymmErrors *grPPG = (TGraphAsymmErrors *) fPPG->Get("Graph");
0061   grPPG->SetFillColorAlpha(kBlack, 0.60);
0062   grPPG->SetFillStyle(1001);
0063   grPPG->SetLineWidth(2);
0064   grPPG->SetMarkerStyle(20);
0065   grPPG->SetMarkerSize(1);
0066   grPPG->SetMarkerColor(kBlack);
0067   grPPG->SetLineColor(kBlack);
0068   grPPG->RemovePoint(0);
0069
0070   TH1D *hPPGH = (TH1D *) fPPG->Get("truth_hist_herwig");
0071   styleObj(hPPGH, kCyan + 1, 41);
0072   hPPGH->SetLineWidth(2);
0073
0074   // ------------------------------------------------------------------
0075   // STAR references
0076   // ------------------------------------------------------------------
0077   TFile *fS = new TFile((inDir + "/HEPData-ins709170-v1-Table_2.root").c_str(), "READ");
0078   TGraphAsymmErrors *grS = (TGraphAsymmErrors *) fS->Get("Table 2/Graph1D_y1");
0079   grS->SetMarkerStyle(20);
0080   grS->SetMarkerColor(kBlack);
0081   grS->SetLineColor(kBlack);
0082
0083   TFile *fN = new TFile((inDir + "/HEPData-ins1771348-v1-Figure_3.root").c_str(), "READ");
0084   TGraphAsymmErrors *grN = (TGraphAsymmErrors *) fN->Get("Figure 3/Graph1D_y1");
0085   grN->SetMarkerStyle(20);
0086   grN->SetMarkerColor(kBlack);
0087   grN->SetLineColor(kBlack);
0088
0089   // ------------------------------------------------------------------
0090   // accumulate the HERWIG samples, once per particle source
0091   // ------------------------------------------------------------------
0092   TH1D *spec[NVERS] = {nullptr, nullptr};  // pi+ spectrum (MB only)
0093   TH2D *N[NVERS] = {nullptr, nullptr};     // transverse dNch/detadphi
0094   TH2D *E[NVERS] = {nullptr, nullptr};     // transverse sumET, R = 0.4
0095
0096   for (int i = 0; i < NSAMP; i++)
0097   {
0098     std::cout << "working on sample " << i << ": " << samp[i] << std::endl;
0099
0100     TFile *fSamp = new TFile(std::format("{}/Herwig_UE_{}_{}.root", inDir, prodTag, samp[i]).c_str(), "READ");
0101     if (!fSamp || fSamp->IsZombie())
0102     {
0103       std::cerr << "  could not open input for " << samp[i] << ", skipping" << std::endl;
0104       continue;
0105     }
0106
0107     TH1D *ev = (TH1D *) fSamp->Get("h_nevents");
0108
0109     for (int v = 0; v < NVERS; v++)
0110     {
0111       const std::string t = "_" + versTag[v];
0112
0113       const double nev = ev->GetBinContent(versNevtBin[v]);
0114       if (nev <= 0)
0115       {
0116         std::cerr << "  " << samp[i] << ": no events for source " << versTag[v]
0117                   << ", skipping" << std::endl;
0118         continue;
0119       }
0120
0121       // per-sample stitching weight still to be applied here
0122       const double sampWeight = csWeightsAlreadyInFiles ? 1.0 : HerwigCS[i];
0123       // extra factor needed to undo a weight already baked into the file
0124       const double undoWeight = csWeightsAlreadyInFiles ? HerwigCS[i] : 1.0;
0125
0126       if (i == 0)
0127       {
0128         spec[v] = (TH1D *) fSamp->Get(("pip_spec" + t).c_str());
0129         if (spec[v])
0130         {
0131           spec[v] = (TH1D *) spec[v]->Clone(("spec_" + versTag[v]).c_str());
0132           spec[v]->SetDirectory(nullptr);
0133           spec[v]->Scale(1.0 / (nev * undoWeight));
0134         }
0135       }
0136
0137       TH2D *NTmp = (TH2D *) fSamp->Get(("p_dens_trans_pt02" + t).c_str());
0138       if (NTmp)
0139       {
0140         NTmp->Scale(sampWeight / nev);
0141         if (!N[v])
0142         {
0143           N[v] = (TH2D *) NTmp->Clone(("N_" + versTag[v]).c_str());
0144           N[v]->SetDirectory(nullptr);
0145         }
0146         else
0147         {
0148           N[v]->Add(NTmp);
0149         }
0150       }
0151
0152       TH2D *ETmp = (TH2D *) fSamp->Get(("p_sumET_trans_R04" + t).c_str());
0153       if (ETmp)
0154       {
0155         ETmp->Scale(sampWeight / nev);
0156         if (!E[v])
0157         {
0158           E[v] = (TH2D *) ETmp->Clone(("E_" + versTag[v]).c_str());
0159           E[v]->SetDirectory(nullptr);
0160         }
0161         else
0162         {
0163           E[v]->Add(ETmp);
0164         }
0165       }
0166     }
0167
0168     fSamp->Close();
0169   }
0170
0171   // ------------------------------------------------------------------
0172   // profiles and styling
0173   // ------------------------------------------------------------------
0174   TProfile *pN[NVERS] = {nullptr, nullptr};
0175   TProfile *pE[NVERS] = {nullptr, nullptr};
0176
0177   for (int v = 0; v < NVERS; v++)
0178   {
0179     if (spec[v])
0180     {
0181       spec[v]->Scale(1.0, "width");
0182       styleObj(spec[v], versColor[v], versMarker[v]);
0183     }
0184
0185     if (N[v])
0186     {
0187       pN[v] = N[v]->ProfileX(("pN_" + versTag[v]).c_str());
0188       pN[v]->Scale(1.0 / (2.0 * 2.0 * TMath::Pi() / 3.0));
0189       styleObj(pN[v], versColor[v], versMarker[v]);
0190     }
0191
0192     if (E[v])
0193     {
0194       pE[v] = E[v]->ProfileX(("pE_" + versTag[v]).c_str());
0195       styleObj(pE[v], versColor[v], versMarker[v]);
0196     }
0197   }
0198
0199   TCanvas *c1 = new TCanvas();
0200
0201   // ------------------------------------------------------------------
0202   // transverse charged-particle density vs STAR PRD 101
0203   // ------------------------------------------------------------------
0204   c1->Clear();
0205   c1->SetLogy(0);
0206
0207   TProfile *pNFrame = pN[0] ? pN[0] : pN[1];
0208   pNFrame->GetYaxis()->SetRangeUser(0.0, 1.5);
0209   pNFrame->Draw("P");
0210   grN->Draw("PSAME");
0211   for (int v = 0; v < NVERS; v++)
0212   {
0213     if (pN[v]) pN[v]->Draw("PSAME");
0214   }
0215
0216   TLegend *leg = new TLegend(0.5, 0.5, 0.85, 0.85);
0217   leg->AddEntry(grN, "STAR PRD 101", "P");
0218   for (int v = 0; v < NVERS; v++)
0219   {
0220     if (pN[v]) leg->AddEntry(pN[v], versLabel[v].c_str(), "P");
0221   }
0222   leg->Draw("same");
0223
0224   c1->SaveAs((inDir + "/multDens.pdf").c_str());
0225
0226   // ------------------------------------------------------------------
0227   // transverse sumET vs PPG10 data and PPG10 HERWIG
0228   // ------------------------------------------------------------------
0229   c1->Clear();
0230   c1->SetLogy(0);
0231
0232   TProfile *pEFrame = pE[0] ? pE[0] : pE[1];
0233   pEFrame->GetYaxis()->SetRangeUser(0.0, 0.85);
0234   pEFrame->GetYaxis()->SetTitle("#LT#Sigma E_{T}/#delta#eta#delta#phi#GT [GeV]");
0235   pEFrame->Draw("P");
0236   grPPG->Draw("E2 SAME");
0237   hPPG->Draw("E1 SAME");
0238   hPPGH->Draw("PSAME");
0239   for (int v = 0; v < NVERS; v++)
0240   {
0241     if (pE[v]) pE[v]->Draw("PSAME");
0242   }
0243
0244   leg = new TLegend(0.15, 0.15, 0.45, 0.4);
0245   leg->AddEntry(grPPG, "sPHENIX PPG10 data", "PE");
0246   leg->AddEntry(hPPGH, "HERWIG PPG10", "P");
0247   for (int v = 0; v < NVERS; v++)
0248   {
0249     if (pE[v]) leg->AddEntry(pE[v], versLabel[v].c_str(), "P");
0250   }
0251   leg->Draw("same");
0252
0253   c1->SaveAs((inDir + "/sumET.pdf").c_str());
0254
0255   // ------------------------------------------------------------------
0256   // pi+ spectrum vs STAR PLB 637, with MC/data ratio panel
0257   // ------------------------------------------------------------------
0258   c1->Clear();
0259
0260   TPad *topPad = new TPad("topPad", "", 0.0, 0.4, 1.0, 1.0);
0261   topPad->SetTopMargin(0.05);
0262   topPad->SetRightMargin(0.05);
0263   topPad->SetBottomMargin(0.0);
0264   topPad->SetLogy();
0265   topPad->Draw();
0266
0267   TPad *bottomPad = new TPad("bottomPad", "", 0.0, 0.0, 1.0, 0.4);
0268   bottomPad->SetTopMargin(0.0);
0269   bottomPad->SetRightMargin(0.05);
0270   bottomPad->SetBottomMargin(0.16);
0271   bottomPad->Draw();
0272
0273   topPad->cd();
0274
0275   // set the y range from data AND MC together, so a mis-normalized MC shows
0276   // up as an off-scale curve instead of silently pushing the data off the pad
0277   double ylo = 1e300, yhi = -1e300;
0278   for (int i = 0; i < grS->GetN(); i++)
0279   {
0280     const double y = grS->GetY()[i];
0281     if (y > 0)
0282     {
0283       ylo = std::min(ylo, y);
0284       yhi = std::max(yhi, y);
0285     }
0286   }
0287   for (int v = 0; v < NVERS; v++)
0288   {
0289     if (!spec[v]) continue;
0290     for (int i = 1; i <= spec[v]->GetNbinsX(); i++)
0291     {
0292       const double y = spec[v]->GetBinContent(i);
0293       if (y > 0)
0294       {
0295         ylo = std::min(ylo, y);
0296         yhi = std::max(yhi, y);
0297       }
0298     }
0299   }
0300
0301   TH1D *specFrame = spec[0] ? spec[0] : spec[1];
0302   specFrame->GetYaxis()->SetRangeUser(0.2 * ylo, 5.0 * yhi);
0303   specFrame->Draw("P");
0304   grS->Draw("P SAME");
0305   for (int v = 0; v < NVERS; v++)
0306   {
0307     if (spec[v]) spec[v]->Draw("P SAME");
0308   }
0309
0310   leg = new TLegend(0.5, 0.5, 0.85, 0.85);
0311   leg->AddEntry(grS, "STAR PLB 637", "P");
0312   for (int v = 0; v < NVERS; v++)
0313   {
0314     if (spec[v]) leg->AddEntry(spec[v], versLabel[v].c_str(), "P");
0315   }
0316   leg->Draw("same");
0317
0318   bottomPad->cd();
0319
0320   TGraphAsymmErrors *grSRat = new TGraphAsymmErrors();
0321   grSRat->SetFillColorAlpha(kYellow, 1.0);
0322   grSRat->SetFillStyle(1001);
0323   grSRat->SetLineWidth(2);
0324   grSRat->SetMarkerStyle(20);
0325   grSRat->SetMarkerSize(0);
0326   grSRat->SetMarkerColor(kYellow);
0327   grSRat->SetLineColor(kYellow);
0328   for (int i = 0; i < grS->GetN(); i++)
0329   {
0330     double x = grS->GetX()[i];
0331     double y = grS->GetY()[i];
0332
0333     grSRat->AddPoint(x, 1.0);
0334     grSRat->SetPointError(i, grS->GetErrorXlow(i), grS->GetErrorXhigh(i),
0335                           grS->GetErrorYlow(i) / y, grS->GetErrorYhigh(i) / y);
0336   }
0337   grSRat->Print();
0338
0339   TH1D *specRat[NVERS] = {nullptr, nullptr};
0340   for (int v = 0; v < NVERS; v++)
0341   {
0342     if (!spec[v]) continue;
0343     specRat[v] = (TH1D *) spec[v]->Clone(("specRat_" + versTag[v]).c_str());
0344     for (int i = 1; i <= spec[v]->GetNbinsX(); i++)
0345     {
0346       const double num = spec[v]->GetBinContent(i);
0347       const double den = grS->GetY()[i - 1];
0348       if (den == 0.0 || num == 0.0)
0349       {
0350         specRat[v]->SetBinContent(i, 0.0);
0351         specRat[v]->SetBinError(i, 0.0);
0352         continue;
0353       }
0354       specRat[v]->SetBinContent(i, num / den);
0355       specRat[v]->SetBinError(i, (num / den) * spec[v]->GetBinError(i) / num);
0356     }
0357   }
0358
0359   TLine *l = new TLine(0.3, 1.0, 10.0, 1.0);
0360   l->SetLineColor(kBlack);
0361
0362   TH1D *ratFrame = specRat[0] ? specRat[0] : specRat[1];
0363   ratFrame->GetYaxis()->SetRangeUser(0.0, 2.0);
0364   ratFrame->GetYaxis()->SetTitle("MC / Data");
0365   ratFrame->SetTitle("");
0366   ratFrame->Draw("P");
0367   grSRat->Draw("E2 SAME");
0368   l->Draw("same");
0369   for (int v = 0; v < NVERS; v++)
0370   {
0371     if (specRat[v]) specRat[v]->Draw("P SAME");
0372   }
0373
0374   c1->SaveAs((inDir + "/pip_spec.pdf").c_str());
0375 }