analysis/HERWIG_UE/PlotUETransverse.C

0001 // PlotUETransverse.C
0002 //
0003 // Reproduce Fig. 4 (top right) of the HERWIG7 RHIC tune paper
0004 // [arXiv:2411.16897]: Transverse-region <dNch/(deta dphi)> (pT > 0.2 GeV/c)
0005 // vs leading-jet pT, HERWIG (Nashville) vs STAR, with an MC/data ratio panel.
0006 //
0007 // STAR data points are PRD 101, 052004 (2020), Figure 2 Transverse region
0008 // (HEPData / Rivet ref STAR_2019_I1771348, d01-x01-y03).
0009 //
0010 //   root -l 'PlotUETransverse.C("herwig_ue_check.root")'
0011
0012 #include <TCanvas.h>
0013 #include <TFile.h>
0014 #include <TGraphAsymmErrors.h>
0015 #include <TH1D.h>
0016 #include <TLatex.h>
0017 #include <TLegend.h>
0018 #include <TLine.h>
0019 #include <TPad.h>
0020 #include <TProfile.h>
0021 #include <TStyle.h>
0022
0023 #include <cmath>
0024 #include <iostream>
0025
0026 void PlotUETransverse(const char *mcfile = "herwig_ue_check.root",
0027                       const char *outname = "ue_transverse_nashville_vs_star")
0028 {
0029   gStyle->SetOptStat(0);
0030   gStyle->SetOptTitle(0);
0031
0032   // ------------------------------------------------------------------
0033   // STAR PRD 101, 052004: Transverse <dNch/detadphi>, pT > 0.2 GeV/c
0034   // ------------------------------------------------------------------
0035   const int nb = 8;
0036   const double edges[nb + 1] = {5., 7., 9., 11., 15., 20., 25., 35., 45.};
0037   const double val[nb] = {0.7585, 0.7454, 0.7195, 0.6641,
0038                           0.5665, 0.5415, 0.5205, 0.4944};
0039   const double stat[nb] = {2.0e-4, 1.0e-4, 2.0e-4, 2.0e-4,
0040                            3.0e-4, 8.0e-4, 1.5e-3, 6.0e-3};
0041   const double sysDn[nb] = {0.060, 0.070, 0.070, 0.050,
0042                             0.033, 0.034, 0.027, 0.033};
0043   const double sysUp[nb] = {0.060, 0.040, 0.050, 0.050,
0044                             0.050, 0.040, 0.032, 0.032};
0045
0046   TGraphAsymmErrors *gData = new TGraphAsymmErrors(nb);
0047   TGraphAsymmErrors *gDataSys = new TGraphAsymmErrors(nb);  // sys band
0048   for (int i = 0; i < nb; ++i)
0049   {
0050     const double xc = 0.5 * (edges[i] + edges[i + 1]);
0051     const double xw = 0.5 * (edges[i + 1] - edges[i]);
0052     gData->SetPoint(i, xc, val[i]);
0053     gData->SetPointError(i, 0., 0., stat[i], stat[i]);
0054     gDataSys->SetPoint(i, xc, val[i]);
0055     gDataSys->SetPointError(i, xw, xw, sysDn[i], sysUp[i]);
0056   }
0057
0058   // ------------------------------------------------------------------
0059   // HERWIG Nashville prediction from the Fun4All module
0060   // ------------------------------------------------------------------
0061   TFile *fin = TFile::Open(mcfile, "READ");
0062   if (!fin || fin->IsZombie())
0063   {
0064     std::cerr << "cannot open " << mcfile << std::endl;
0065     return;
0066   }
0067   TProfile *prof = dynamic_cast<TProfile *>(fin->Get("p_dens_trans_pt02"));
0068   if (!prof)
0069   {
0070     std::cerr << "p_dens_trans_pt02 not found in " << mcfile << std::endl;
0071     return;
0072   }
0073   TH1D *hMC = prof->ProjectionX("h_mc_trans");
0074   hMC->SetDirectory(nullptr);
0075
0076   // ratio MC/data with the data uncertainty band around unity
0077   TH1D *hRatio = (TH1D *) hMC->Clone("h_ratio");
0078   hRatio->SetDirectory(nullptr);
0079   TGraphAsymmErrors *gBand = new TGraphAsymmErrors(nb);
0080   for (int i = 0; i < nb; ++i)
0081   {
0082     const double d = val[i];
0083     hRatio->SetBinContent(i + 1, hMC->GetBinContent(i + 1) / d);
0084     hRatio->SetBinError(i + 1, hMC->GetBinError(i + 1) / d);
0085     const double xc = 0.5 * (edges[i] + edges[i + 1]);
0086     const double xw = 0.5 * (edges[i + 1] - edges[i]);
0087     const double eDn = std::sqrt(stat[i] * stat[i] + sysDn[i] * sysDn[i]) / d;
0088     const double eUp = std::sqrt(stat[i] * stat[i] + sysUp[i] * sysUp[i]) / d;
0089     gBand->SetPoint(i, xc, 1.);
0090     gBand->SetPointError(i, xw, xw, eDn, eUp);
0091   }
0092
0093   // ------------------------------------------------------------------
0094   // draw
0095   // ------------------------------------------------------------------
0096   TCanvas *c = new TCanvas("c", "", 700, 750);
0097   TPad *pTop = new TPad("pTop", "", 0., 0.33, 1., 1.);
0098   TPad *pBot = new TPad("pBot", "", 0., 0., 1., 0.33);
0099   pTop->SetBottomMargin(0.02);
0100   pBot->SetTopMargin(0.03);
0101   pBot->SetBottomMargin(0.30);
0102   pTop->Draw();
0103   pBot->Draw();
0104
0105   pTop->cd();
0106   TH1D *frame = new TH1D("frame", ";;#LTdN_{ch}/d#etad#phi#GT", nb, edges);
0107   frame->GetYaxis()->SetRangeUser(0., 1.1);
0108   frame->GetYaxis()->SetTitleSize(0.06);
0109   frame->GetYaxis()->SetLabelSize(0.05);
0110   frame->GetXaxis()->SetLabelSize(0.);
0111   frame->Draw("AXIS");
0112
0113   gDataSys->SetFillColorAlpha(kOrange - 2, 0.5);
0114   gDataSys->Draw("2 SAME");
0115   gData->SetMarkerStyle(20);
0116   gData->SetMarkerSize(1.1);
0117   gData->SetMarkerColor(kBlack);
0118   gData->SetLineColor(kBlack);
0119   gData->Draw("P SAME");
0120
0121   hMC->SetLineColor(kRed + 1);
0122   hMC->SetMarkerColor(kRed + 1);
0123   hMC->SetMarkerStyle(24);
0124   hMC->SetLineWidth(2);
0125   hMC->Draw("E1 SAME");
0126
0127   TLegend *leg = new TLegend(0.45, 0.62, 0.88, 0.87);
0128   leg->SetBorderSize(0);
0129   leg->SetFillStyle(0);
0130   leg->AddEntry(gData, "STAR, PRD 101, 052004", "p");
0131   leg->AddEntry(gDataSys, "syst. uncertainty", "f");
0132   leg->AddEntry(hMC, "HERWIG 7.3 Nashville (sPHENIX prod.)", "lp");
0133   leg->Draw();
0134
0135   TLatex tx;
0136   tx.SetNDC();
0137   tx.SetTextSize(0.05);
0138   tx.DrawLatex(0.14, 0.85, "p+p #sqrt{s} = 200 GeV");
0139   tx.DrawLatex(0.14, 0.78, "Transverse region, p_{T}^{ch} > 0.2 GeV/c");
0140   tx.SetTextSize(0.042);
0141   tx.DrawLatex(0.14, 0.71, "anti-k_{T} R = 0.6, |#eta_{jet}| < 0.4, |#eta_{ch}| < 1");
0142
0143   pBot->cd();
0144   TH1D *rframe = new TH1D("rframe",
0145                           ";leading jet p_{T} [GeV/c];MC / data", nb, edges);
0146   rframe->GetYaxis()->SetRangeUser(0.5, 1.5);
0147   rframe->GetYaxis()->SetNdivisions(505);
0148   rframe->GetYaxis()->SetTitleSize(0.11);
0149   rframe->GetYaxis()->SetTitleOffset(0.5);
0150   rframe->GetYaxis()->SetLabelSize(0.10);
0151   rframe->GetXaxis()->SetTitleSize(0.12);
0152   rframe->GetXaxis()->SetLabelSize(0.10);
0153   rframe->Draw("AXIS");
0154
0155   gBand->SetFillColorAlpha(kOrange - 2, 0.5);
0156   gBand->Draw("2 SAME");
0157   TLine unity(edges[0], 1., edges[nb], 1.);
0158   unity.SetLineStyle(2);
0159   unity.DrawClone("SAME");
0160
0161   hRatio->SetLineColor(kRed + 1);
0162   hRatio->SetMarkerColor(kRed + 1);
0163   hRatio->SetMarkerStyle(24);
0164   hRatio->SetLineWidth(2);
0165   hRatio->Draw("E1 SAME");
0166
0167   c->SaveAs(Form("%s.pdf", outname));
0168   c->SaveAs(Form("%s.png", outname));
0169
0170   // quick numerical summary
0171   std::cout << "\n  bin [GeV/c]      STAR        HERWIG      MC/data\n";
0172   for (int i = 0; i < nb; ++i)
0173   {
0174     printf("  %4.0f - %4.0f    %.4f      %.4f      %.3f\n",
0175            edges[i], edges[i + 1], val[i], hMC->GetBinContent(i + 1),
0176            hMC->GetBinContent(i + 1) / val[i]);
0177   }
0178 }