Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:12:47 

0001 // $Id: $
0002 
0003 /*!
0004  * \file Draw_HFJetTruth.C
0005  * \brief 
0006  * \author Jin Huang <jhuang@bnl.gov>
0007  * \version $Revision:   $
0008  * \date $Date: $
0009  */
0010 
0011 #include <TChain.h>
0012 #include <TFile.h>
0013 #include <TGraphAsymmErrors.h>
0014 #include <TGraphErrors.h>
0015 #include <TLatex.h>
0016 #include <TLegend.h>
0017 #include <TLine.h>
0018 #include <TString.h>
0019 #include <TTree.h>
0020 #include <TVirtualFitter.h>
0021 #include <cassert>
0022 #include <cmath>
0023 #include <vector>
0024 #include "SaveCanvas.C"
0025 #include "sPhenixStyle.C"
0026 
0027 using namespace std;
0028 
0029 TGraphAsymmErrors *
0030 GetFONLL_C();
0031 TGraphAsymmErrors *
0032 GetFONLL_B();
0033 TGraph *pQCDModel_HuangKangVitev(const double g);
0034 TGraph *LIDO_Ke();
0035 TGraph *
0036 GetPHENIX_jet();
0037 
0038 TFile *_file0 = NULL;
0039 TTree *T = NULL;
0040 
0041 TH1 *CrossSection2RelUncert(const TH1F *h_cross,
0042                             const double suppression,
0043                             const double deta,
0044                             const double pp_quiv_int_lum)
0045 {
0046   assert(h_cross);
0047   TH1 *
0048       h_ratio = (TH1 *)
0049                     h_cross->Clone(TString(h_cross->GetName()) + Form("_copy%d", rand()));
0050 
0051   //convert to count per bin
0052   h_ratio->Scale(deta * h_ratio->GetXaxis()->GetBinWidth(0) * pp_quiv_int_lum * suppression);
0053   h_ratio->Rebin(5);
0054 
0055   for (int i = 1; i <= h_ratio->GetNbinsX(); ++i)
0056   {
0057     const double yield = h_ratio->GetBinContent(i);
0058 
0059     if (yield > 0)
0060     {
0061       h_ratio->SetBinContent(i, suppression);
0062 
0063       h_ratio->SetBinError(i, suppression / sqrt(yield));
0064     }
0065     else
0066     {
0067       h_ratio->SetBinContent(i, 0);
0068 
0069       h_ratio->SetBinError(i, 0);
0070     }
0071   }
0072 
0073   h_ratio->GetYaxis()->SetTitle("Relative Cross Section and Uncertainty");
0074 
0075   return h_ratio;
0076 }
0077 
0078 TGraphErrors *GetRAA(TH1 *h_pp, TH1 *h_AA)
0079 {
0080   int n_bin = 0;
0081 
0082   double xs[1000] = {0};
0083   double ys[1000] = {0};
0084   double eys[1000] = {0};
0085 
0086   assert(h_pp);
0087   assert(h_AA);
0088   assert(h_pp->GetNbinsX() == h_AA->GetNbinsX());
0089 
0090   for (int i = 1; i <= h_pp->GetNbinsX(); ++i)
0091   {
0092     if (h_pp->GetBinError(i) > 0 && h_pp->GetBinError(i) < 1 && h_AA->GetBinError(i) > 0 && h_AA->GetBinError(i) < 1)
0093     {
0094       xs[n_bin] = h_pp->GetXaxis()->GetBinCenter(i);
0095 
0096       ys[n_bin] = h_AA->GetBinContent(i) / h_pp->GetBinContent(i);
0097 
0098       eys[n_bin] = ys[n_bin] * sqrt(pow(h_AA->GetBinError(i) / h_AA->GetBinContent(i), 2) + pow(h_pp->GetBinError(i) / h_pp->GetBinContent(i), 2));
0099 
0100       n_bin += 1;
0101     }
0102   }
0103 
0104   TGraphErrors *ge = new TGraphErrors(n_bin, xs, ys, NULL, eys);
0105   ge->SetName(TString("RAA_") + h_AA->GetName());
0106 
0107   return ge;
0108 }
0109 
0110 void CrossSection2RAA(const TString infile, const bool use_AA_jet_trigger = true, const double dy = .7 * 2)
0111 {
0112   TFile *f = TFile::Open(infile + "Draw_HFJetTruth_DrawCrossSection.root");
0113   assert(f);
0114 
0115   TH1F *hall = (TH1F *) f->GetObjectChecked("hall", "TH1F");
0116   assert(hall);
0117   TH1F *h_b = (TH1F *) f->GetObjectChecked("h_b", "TH1F");
0118   assert(h_b);
0119 
0120   TString s_suffix(use_AA_jet_trigger ? "_AAJetTriggered" : "_3yr");
0121   s_suffix += TString(Form("_deta%.2f", dy / 2));
0122 
0123   const double b_jet_RAA = 0.6;
0124 
0125   const double pp_eff = 0.6;
0126   const double pp_purity = 0.4;
0127   const double AuAu_eff = 0.4;
0128   const double AuAu_purity = 0.4;
0129 
0130   ////////////////////////////
0131   // 5-year lumi in [sPH-TRG-000]
0132   ////////////////////////////
0133 
0134   const double pp_lumi = 62;                           // pb^-1 [sPH-TRG-000], rounded up from 197 pb^-1
0135   const double pp_inelastic_crosssec = 42e-3 / 1e-12;  // 42 mb in pb [sPH-TRG-000]
0136 
0137   const double AuAu_MB_Evt = use_AA_jet_trigger ? 218e9 : 142e9;  // [sPH-TRG-000], depending on whether jet trigger applied in AA collisions
0138   const double pAu_MB_Evt = 600e9;                                // [sPH-TRG-000]
0139 
0140   const double AuAu_Ncoll_C0_10 = 960.2;  // [DOI:?10.1103/PhysRevC.87.034911?]
0141   const double AuAu_Ncoll_C0_20 = 770.6;  // [DOI:?10.1103/PhysRevC.91.064904?]
0142   const double AuAu_Ncoll_C0_100 = 250;   // pb^-1 [sPH-TRG-000]
0143   const double pAu_Ncoll_C0_100 = 4.7;    // pb^-1 [sPH-TRG-000]
0144 
0145   const double AuAu_eq_lumi_C0_10 = AuAu_MB_Evt * 0.1 * AuAu_Ncoll_C0_10 / pp_inelastic_crosssec;  //
0146   const double AuAu_eq_lumi_C0_20 = AuAu_MB_Evt * 0.2 * AuAu_Ncoll_C0_20 / pp_inelastic_crosssec;  //
0147   const double AuAu_eq_lumi_C0_100 = AuAu_MB_Evt * 1 * AuAu_Ncoll_C0_100 / pp_inelastic_crosssec;  //
0148 
0149   const double pAu_eq_lumi_C0_100 = pAu_MB_Evt * 1 * pAu_Ncoll_C0_100 / pp_inelastic_crosssec;  //
0150 
0151   cout << "CrossSection2RAA integrated luminosity assumptions in pb^-1: " << endl;
0152   cout << "\t"
0153        << "pp_lumi = " << pp_lumi << endl;
0154   cout << "\t"
0155        << "AuAu_eq_lumi_C0_10 = " << AuAu_eq_lumi_C0_10 << endl;
0156   cout << "\t"
0157        << "AuAu_eq_lumi_C0_20 = " << AuAu_eq_lumi_C0_20 << endl;
0158   cout << "\t"
0159        << "AuAu_eq_lumi_C0_100 = " << AuAu_eq_lumi_C0_100 << endl;
0160   cout << "\t"
0161        << "pAu_eq_lumi_C0_100 = " << pAu_eq_lumi_C0_100 << endl;
0162 
0163   TCanvas *c1 = new TCanvas("Draw_HFJetTruth_CrossSection2RAA_Ratio" + s_suffix, "Draw_HFJetTruth_CrossSection2RAA_Ratio" + s_suffix, 700, 600);
0164   c1->Divide(1, 1);
0165   int idx = 1;
0166   TPad *p;
0167 
0168   p = (TPad *) c1->cd(idx++);
0169   c1->Update();
0170   //  p->SetGridx(0);
0171   //  p->SetGridy(0);
0172   //  p->SetLogy();
0173 
0174   TH1 *g_pp = CrossSection2RelUncert(h_b, 1., dy, pp_lumi * pp_eff * pp_purity);
0175   TH1 *g_AA = CrossSection2RelUncert(h_b, b_jet_RAA, dy, AuAu_eq_lumi_C0_10 * AuAu_eff * AuAu_purity);
0176 
0177   g_pp->SetLineColor(kRed);
0178   g_AA->SetLineColor(kBlue);
0179 
0180   g_pp->Draw();
0181   g_AA->Draw("same");
0182   SaveCanvas(c1, infile + "_" + TString(c1->GetName()), kTRUE);
0183 
0184   c1 = new TCanvas("Draw_HFJetTruth_CrossSection2RAA" + s_suffix, "Draw_HFJetTruth_CrossSection2RAA" + s_suffix, 1100, 800);
0185   c1->Divide(1, 1);
0186   idx = 1;
0187 
0188   p = (TPad *) c1->cd(idx++);
0189   c1->Update();
0190 
0191   p->DrawFrame(11, 0, 45, 1.2)
0192       ->SetTitle(";Transverse Momentum [GeV/#it{c}];#it{R}_{AA}");
0193 
0194   TGraphErrors *ge_RAA = GetRAA(g_pp, g_AA);
0195 
0196   ge_RAA->SetLineWidth(4);
0197   ge_RAA->SetMarkerStyle(kFullCircle);
0198   ge_RAA->SetMarkerSize(3);
0199 
0200   ge_RAA->Draw("pe");
0201   ge_RAA->Print();
0202 
0203   TLegend *leg = new TLegend(.0, .67, .85, .91);
0204   leg->SetFillStyle(0);
0205   leg->AddEntry("", "#it{#bf{sPHENIX}} Projection", "");
0206   //  leg->AddEntry("", Form("PYTHIA-8 #it{b}-jet, Anti-k_{T} R=0.4, |#eta|<%.2f, CTEQ6L", dy / 2), "");
0207   leg->AddEntry("", Form("#it{b}-jet Anti-k_{T} R=0.4, 0-10%% Au+Au, Year 1-3"), "");
0208   leg->AddEntry("", Form("#it{p}+#it{p}: %.0fpb^{-1} samp., %.0f%% Eff., %.0f%% Pur.", pp_lumi, pp_eff * 100, pp_purity * 100), "");
0209   leg->AddEntry("", Form("Au+Au: %.0fnb^{-1} rec., %.0f%% Eff., %.0f%% Pur.", AuAu_MB_Evt / 6.8252 / 1e9, AuAu_eff * 100, AuAu_purity * 100), "");
0210   leg->Draw();
0211 
0212   TLegend *leg2 = new TLegend(.17, .70, .88, .77);
0213   leg2->AddEntry(ge_RAA, "#it{b}-jet #it{R}_{AA}, Au+Au 0-10%C, #sqrt{s_{NN}}=200 GeV", "pl");
0214   //  leg2->Draw();
0215 
0216   SaveCanvas(c1, infile + "_" + TString(c1->GetName()), kTRUE);
0217 
0218   c1 = new TCanvas("Draw_HFJetTruth_CrossSection2RAA_Theory" + s_suffix, "Draw_HFJetTruth_CrossSection2RAA_Theory" + s_suffix, 1100, 800);
0219   c1->Divide(1, 1);
0220   idx = 1;
0221   *p;
0222 
0223   p = (TPad *) c1->cd(idx++);
0224   c1->Update();
0225 
0226   p->DrawFrame(11, 0, 45, 1.2)
0227       ->SetTitle(";p_{T} [GeV];#it{R}_{AA}");
0228 
0229   TGraph *gLIDO = LIDO_Ke();
0230   gLIDO->SetFillColor(kAzure - 9);
0231   gLIDO->Draw("3");
0232 
0233   TGraph *g20 = pQCDModel_HuangKangVitev(2.0);
0234   TGraph *g22 = pQCDModel_HuangKangVitev(2.2);
0235 
0236   g20->SetLineColor(kGreen + 3);
0237   g22->SetLineColor(kRed + 3);
0238   g20->SetLineWidth(3);
0239   g22->SetLineWidth(3);
0240 
0241   g20->Draw("l");
0242   g22->Draw("l");
0243 
0244   ge_RAA->DrawClone("pe");
0245   leg->DrawClone();
0246   //  leg2->DrawClone();
0247 
0248   TLegend *leg1 = new TLegend(.195, .17, .7, .3);
0249   //  leg1->SetHeader("#splitline{pQCD, Phys.Lett. B726 (2013) 251-256}{#sqrt{s_{NN}}=200 GeV, #it{b}-jet R=0.4}");
0250   leg1->SetHeader("pQCD, Phys.Lett. B726 (2013) 251-256");
0251   leg1->AddEntry(g20, "g^{med} = 2.0", "l");
0252   leg1->AddEntry(g22, "g^{med} = 2.2", "l");
0253   leg1->Draw();
0254 
0255   TLegend *leg_lido = new TLegend(.2, .32, .55, .37);
0256   leg_lido->AddEntry(gLIDO, "LIDO, arXiv:2008.07622 [nucl-th]", "f");
0257   leg_lido->Draw();
0258 
0259   SaveCanvas(c1, infile + "_" + TString(c1->GetName()), kTRUE);
0260 }
0261 
0262 TGraphErrors *GetZgPrediction(const TString data_config = "me_bbg")
0263 {
0264   // Inverting the mass hierarchy of jet quenching effects with prompt bb-jet substructure, Hai Tao Li(Los Alamos Natl. Lab., Theor. Div.), Ivan Vitev(Los Alamos Natl. Lab., Theor. Div.)
0265   // DOI: 10.1016/j.physletb.2019.04.052 (publication)
0266   // Li, Haitao <haitaoli@lanl.gov>
0267   // Here are the data for the predicted zg distributions.  For all the files the first column is zg value. The second and third columns are the higher and lower limits of the uncertainty bands.
0268 
0269   std::vector<double> zgs;
0270   std::vector<double> ymaxs;
0271   std::vector<double> ymins;
0272 
0273   std::vector<double> ys;
0274   std::vector<double> eys;
0275 
0276   if (data_config == "me_bbg")
0277   {
0278     zgs = {0.1,
0279            0.11,
0280            0.12,
0281            0.13,
0282            0.14,
0283            0.15,
0284            0.16,
0285            0.17,
0286            0.18,
0287            0.19,
0288            0.2,
0289            0.21,
0290            0.22,
0291            0.23,
0292            0.24,
0293            0.25,
0294            0.26,
0295            0.27,
0296            0.28,
0297            0.29,
0298            0.3,
0299            0.31,
0300            0.32,
0301            0.33,
0302            0.34,
0303            0.35,
0304            0.36,
0305            0.37,
0306            0.38,
0307            0.39,
0308            0.4,
0309            0.41,
0310            0.42,
0311            0.43,
0312            0.44,
0313            0.45,
0314            0.46,
0315            0.47,
0316            0.48,
0317            0.49,
0318            0.5};
0319 
0320     ymaxs = {11.564966612100944,
0321              10.228774883485626,
0322              9.040211253937503,
0323              8.00426321286891,
0324              7.085316588030289,
0325              6.304057794111209,
0326              5.620901023346734,
0327              5.020252217399414,
0328              4.471418358587862,
0329              4.006557331982638,
0330              3.6029776961294937,
0331              3.2746530582558067,
0332              2.979082995195153,
0333              2.717548121923332,
0334              2.484039405367637,
0335              2.2743520629543212,
0336              2.0857993148942358,
0337              1.916461379592078,
0338              1.7644938665497385,
0339              1.6282027691535934,
0340              1.5063707131889736,
0341              1.3971555437989853,
0342              1.2987288248009365,
0343              1.2107212053236827,
0344              1.131621650833046,
0345              1.0605838573254032,
0346              0.9969458228908805,
0347              0.941899699579719,
0348              0.8934983992713553,
0349              0.8508943534347049,
0350              0.8130937637364388,
0351              0.7797214265194943,
0352              0.7506054526040817,
0353              0.7255735799970017,
0354              0.7043026708717507,
0355              0.6866658867399881,
0356              0.6724855329238311,
0357              0.6615441971815835,
0358              0.6538330585599965,
0359              0.6493379751029307,
0360              0.6477595200387519};
0361     ymins = {8.872885875265245,
0362              8.111748205961844,
0363              7.432473620102464,
0364              6.797408173274767,
0365              6.186769713061456,
0366              5.608088106815033,
0367              5.084027720317393,
0368              4.610111938188973,
0369              4.187760399167046,
0370              3.7940340683954212,
0371              3.43117457064948,
0372              3.093119987479635,
0373              2.7861576826055696,
0374              2.5128579366528596,
0375              2.270447062142767,
0376              2.055432334380884,
0377              1.8652346826114092,
0378              1.6966233161676936,
0379              1.5473576814988068,
0380              1.4149243923530526,
0381              1.2976881530156943,
0382              1.1933534033076159,
0383              1.1006617243222991,
0384              1.018247797082035,
0385              0.9450088178795155,
0386              0.8800255813710026,
0387              0.8223509153369435,
0388              0.7710124745428798,
0389              0.7244026618659434,
0390              0.6749851216252735,
0391              0.6388960527749336,
0392              0.6074597301268359,
0393              0.5801575320223376,
0394              0.5567239908569284,
0395              0.5368940589741338,
0396              0.5204904428643295,
0397              0.5073303543085547,
0398              0.4972605201973995,
0399              0.4901298441121407,
0400              0.4859016247798608,
0401              0.4844717994196453};
0402   }
0403   else if (data_config == "bbg_zg")
0404   {
0405     zgs = {0.1,
0406            0.11,
0407            0.12,
0408            0.13,
0409            0.14,
0410            0.15,
0411            0.16,
0412            0.17,
0413            0.18,
0414            0.19,
0415            0.2,
0416            0.21,
0417            0.22,
0418            0.23,
0419            0.24,
0420            0.25,
0421            0.26,
0422            0.27,
0423            0.28,
0424            0.29,
0425            0.3,
0426            0.31,
0427            0.32,
0428            0.33,
0429            0.34,
0430            0.35,
0431            0.36,
0432            0.37,
0433            0.38,
0434            0.39,
0435            0.4,
0436            0.41,
0437            0.42,
0438            0.43,
0439            0.44,
0440            0.45,
0441            0.46,
0442            0.47,
0443            0.48,
0444            0.49,
0445            0.5};
0446 
0447     ymaxs = {9.762434858721084,
0448              8.47616359817982,
0449              7.492676243945642,
0450              6.654010010935872,
0451              5.913925919575108,
0452              5.27028638644749,
0453              4.713862130785145,
0454              4.2338140867044,
0455              3.8359176524435012,
0456              3.5411710604457696,
0457              3.2780402997684046,
0458              3.0428481843445376,
0459              2.8313528228811093,
0460              2.64115838360015,
0461              2.4819574262280666,
0462              2.338567310713753,
0463              2.2090001831145045,
0464              2.091727244023134,
0465              1.9854179724262009,
0466              1.888790011179277,
0467              1.8008939425055928,
0468              1.7209392536426298,
0469              1.6482358047754673,
0470              1.582167173395633,
0471              1.5221883540849073,
0472              1.467818836564728,
0473              1.4186285882228657,
0474              1.3742395815301207,
0475              1.3343146951863025,
0476              1.2985549222721922,
0477              1.266698029272377,
0478              1.2385132145874334,
0479              1.2137978220476016,
0480              1.1923789991503813,
0481              1.1741068974344762,
0482              1.1588532760513377,
0483              1.1465131698787123,
0484              1.1370016290107832,
0485              1.1302536619705346,
0486              1.1262226895802339,
0487              1.124882446490406};
0488     ymins = {7.5042886917903475,
0489              6.698904277963574,
0490              6.13325229167674,
0491              5.642204967225616,
0492              5.184676746368783,
0493              4.767489013960836,
0494              4.391074470600173,
0495              4.0530921099573165,
0496              3.747687167866155,
0497              3.4368852671586723,
0498              3.1391359733756863,
0499              2.872328619918876,
0500              2.637262134930488,
0501              2.430875459603178,
0502              2.24954656566354,
0503              2.090040255967163,
0504              1.9492781360305134,
0505              1.8247013582708238,
0506              1.7141703670450044,
0507              1.61589253533197,
0508              1.528349666705045,
0509              1.4502575139246554,
0510              1.3805229398832255,
0511              1.318209032398018,
0512              1.2625160741157493,
0513              1.2127547858323617,
0514              1.1683292385783044,
0515              1.1287292083339666,
0516              1.0935097910299434,
0517              1.0622884717784657,
0518              1.034734064269136,
0519              1.0105633112172412,
0520              0.9895286076610499,
0521              0.9714187985395711,
0522              0.9560602609384745,
0523              0.9433012733218156,
0524              0.9330225001525349,
0525              0.9251255015159976,
0526              0.9195356933679992,
0527              0.9162028180740642,
0528              0.9150954615243209};
0529   }
0530   else if (data_config == "bbg_zg_ratio")
0531   {
0532     zgs = {0.1,
0533            0.11,
0534            0.12,
0535            0.13,
0536            0.14,
0537            0.15,
0538            0.16,
0539            0.17,
0540            0.18,
0541            0.19,
0542            0.2,
0543            0.21,
0544            0.22,
0545            0.23,
0546            0.24,
0547            0.25,
0548            0.26,
0549            0.27,
0550            0.28,
0551            0.29,
0552            0.3,
0553            0.31,
0554            0.32,
0555            0.33,
0556            0.34,
0557            0.35,
0558            0.36,
0559            0.37,
0560            0.38,
0561            0.39,
0562            0.4,
0563            0.41,
0564            0.42,
0565            0.43,
0566            0.44,
0567            0.45,
0568            0.46,
0569            0.47,
0570            0.48,
0571            0.49,
0572            0.5};
0573 
0574     ymaxs = {1.151982890151877,
0575              1.1664826852504495,
0576              1.172251389758727,
0577              1.1741255109207875,
0578              1.1714608638124098,
0579              1.1681087817987694,
0580              1.1632495306605235,
0581              1.155513671471374,
0582              1.1397719364148045,
0583              1.127051000823432,
0584              1.1119535300939027,
0585              1.0948424254600155,
0586              1.0760601346878573,
0587              1.055878624941244,
0588              1.0345055306535231,
0589              1.0120360659743335,
0590              0.9890362761525937,
0591              0.9655962994186251,
0592              0.9421308194393136,
0593              0.9186226313722144,
0594              0.895503201582456,
0595              0.8725753749800488,
0596              0.8501317891243333,
0597              0.8282275629034408,
0598              0.8069766690984213,
0599              0.7865577997621885,
0600              0.7669973631487963,
0601              0.7481674017140337,
0602              0.7304488459867333,
0603              0.7015931864543251,
0604              0.688687461757641,
0605              0.6766716415600889,
0606              0.6657247620008147,
0607              0.6559757617778346,
0608              0.6474142123454508,
0609              0.6401392451833807,
0610              0.6341568146936112,
0611              0.6294320803863227,
0612              0.6260810257194698,
0613              0.6241574732399557,
0614              0.6234333158108032};
0615     ymins = {1.096464697907084,
0616              1.1148057864579233,
0617              1.1277713781468222,
0618              1.136191604396884,
0619              1.1385826379515291,
0620              1.135722028529819,
0621              1.1284916532508826,
0622              1.1169591429988097,
0623              1.103666045140106,
0624              1.0854505401169579,
0625              1.0654328782085087,
0626              1.0442210459999413,
0627              1.0223533047937827,
0628              1.0001575721016631,
0629              0.9777309842954044,
0630              0.9550386037145693,
0631              0.9325222701437524,
0632              0.9101431388163338,
0633              0.8882340053200084,
0634              0.8667312777949217,
0635              0.8459583146041589,
0636              0.8256251575966925,
0637              0.8059969846816892,
0638              0.7869325271199851,
0639              0.7684747631437073,
0640              0.750777564978232,
0641              0.733824152526011,
0642              0.7174732270599817,
0643              0.7020666530820847,
0644              0.6931802179916633,
0645              0.6795859695799639,
0646              0.6666218950521552,
0647              0.6546221226148139,
0648              0.6438152297250573,
0649              0.6342654764486364,
0650              0.6261010548445681,
0651              0.6193616959253557,
0652              0.6140461556406446,
0653              0.6102701796326937,
0654              0.6080839713862733,
0655              0.6072750596686409};
0656   }
0657   else
0658   {
0659     cout << __PRETTY_FUNCTION__ << "Invalid configuration : " << data_config << endl;
0660     exit(1);
0661   }
0662 
0663   double integral = 0;
0664   for (int i = 0; i < zgs.size(); ++i)
0665   {
0666     assert(i < ymaxs.size());
0667     assert(i < ymins.size());
0668 
0669     ys.push_back(0.5 * (ymaxs[i] + ymins[i]));
0670     eys.push_back(0.5 * fabs(ymaxs[i] - ymins[i]));
0671 
0672     integral += (zgs[1] - zgs[0]) * ys[i];
0673   }
0674   cout << __PRETTY_FUNCTION__ << " configuration : " << data_config << " integral = " << integral << endl;
0675 
0676   TGraphErrors *ge = new TGraphErrors(zgs.size(), zgs.data(), ys.data(), NULL, eys.data());
0677   ge->SetName(TString("zg_") + data_config);
0678 
0679   return ge;
0680 }
0681 
0682 TH1 *CrossSection2zg_HistogramMaker(const TH1F *h_cross,
0683                                     const double suppression,
0684                                     const double deta, const double pt_max,
0685                                     const double pp_quiv_int_lum, const TString data_config)
0686 {
0687   assert(h_cross);
0688   TH1 *
0689       h_yield = (TH1 *)
0690                     h_cross->Clone(TString(h_cross->GetName()) + Form("_copy%d", rand()));
0691 
0692   //convert to count per bin
0693   h_yield->Scale(deta * h_yield->GetXaxis()->GetBinWidth(0) * pp_quiv_int_lum * suppression);
0694   //  h_ratio->Rebin(5);
0695 
0696   double sum_yield = 0;
0697   for (int i = 1; i <= h_yield->GetNbinsX(); ++i)
0698   {
0699     const double yield = h_yield->GetBinContent(i);
0700     if (h_yield->GetXaxis()->GetBinCenter(i) < pt_max) sum_yield += yield;
0701   }
0702 
0703   TGraphErrors *ge_prediction = GetZgPrediction(data_config);
0704   const double bin_width = 0.05;
0705   TH1 *h_zg = new TH1F("hzg_" + data_config, "z_{g}", 0.4 / bin_width, 0.1, 0.5);
0706 
0707   for (int i = 1; i <= h_zg->GetNbinsX(); ++i)
0708   {
0709     const double x_center = h_zg->GetBinCenter(i);
0710 
0711     const double prediction = ge_prediction->Eval(x_center);
0712 
0713     h_zg->SetBinContent(i, prediction);
0714 
0715     const double bin_yield = bin_width * prediction * sum_yield;
0716 
0717     h_zg->SetBinError(i, prediction / sqrt(bin_yield));
0718   }
0719 
0720   return h_zg;
0721 }
0722 
0723 void CrossSection2zg(const TString infile, const bool use_AA_jet_trigger = true, const double dy = .7 * 2)
0724 {
0725   TFile *f = TFile::Open(infile + "Draw_HFJetTruth_DrawCrossSection.root");
0726   assert(f);
0727 
0728   TH1F *hall = (TH1F *) f->GetObjectChecked("hall", "TH1F");
0729   assert(hall);
0730   TH1F *h_b = (TH1F *) f->GetObjectChecked("h_b", "TH1F");
0731   assert(h_b);
0732 
0733   TString s_suffix(use_AA_jet_trigger ? "_AAJetTriggered" : "_3yr");
0734   s_suffix += TString(Form("_deta%.2f", dy / 2));
0735 
0736   const double pt_max = 30;  // max pT to be consistent with DOI: 10.1016/j.physletb.2019.04.052 (publication)
0737 
0738   const double b_jet_RAA = 0.6;
0739 
0740   const double pp_eff = 0.6;
0741   const double pp_purity = 0.4;
0742   const double AuAu_eff = 0.4;
0743   const double AuAu_purity = 0.4;
0744 
0745   ////////////////////////////
0746   // 5-year lumi in [sPH-TRG-000]
0747   ////////////////////////////
0748 
0749   const double pp_lumi = 62;                           // pb^-1 [sPH-TRG-000], rounded up from 197 pb^-1
0750   const double pp_inelastic_crosssec = 42e-3 / 1e-12;  // 42 mb in pb [sPH-TRG-000]
0751 
0752   const double AuAu_MB_Evt = use_AA_jet_trigger ? 218e9 : 142e9;  // [sPH-TRG-000], depending on whether jet trigger applied in AA collisions
0753   const double pAu_MB_Evt = 600e9;                                // [sPH-TRG-000]
0754 
0755   const double AuAu_Ncoll_C0_10 = 960.2;  // [DOI:?10.1103/PhysRevC.87.034911?]
0756   const double AuAu_Ncoll_C0_20 = 770.6;  // [DOI:?10.1103/PhysRevC.91.064904?]
0757   const double AuAu_Ncoll_C0_100 = 250;   // pb^-1 [sPH-TRG-000]
0758   const double pAu_Ncoll_C0_100 = 4.7;    // pb^-1 [sPH-TRG-000]
0759 
0760   const double AuAu_eq_lumi_C0_10 = AuAu_MB_Evt * 0.1 * AuAu_Ncoll_C0_10 / pp_inelastic_crosssec;  //
0761   const double AuAu_eq_lumi_C0_20 = AuAu_MB_Evt * 0.2 * AuAu_Ncoll_C0_20 / pp_inelastic_crosssec;  //
0762   const double AuAu_eq_lumi_C0_100 = AuAu_MB_Evt * 1 * AuAu_Ncoll_C0_100 / pp_inelastic_crosssec;  //
0763 
0764   const double pAu_eq_lumi_C0_100 = pAu_MB_Evt * 1 * pAu_Ncoll_C0_100 / pp_inelastic_crosssec;  //
0765 
0766   cout << "CrossSection2zg integrated luminosity assumptions in pb^-1: " << endl;
0767   cout << "\t"
0768        << "pp_lumi = " << pp_lumi << endl;
0769   cout << "\t"
0770        << "AuAu_eq_lumi_C0_10 = " << AuAu_eq_lumi_C0_10 << endl;
0771   cout << "\t"
0772        << "AuAu_eq_lumi_C0_20 = " << AuAu_eq_lumi_C0_20 << endl;
0773   cout << "\t"
0774        << "AuAu_eq_lumi_C0_100 = " << AuAu_eq_lumi_C0_100 << endl;
0775   cout << "\t"
0776        << "pAu_eq_lumi_C0_100 = " << pAu_eq_lumi_C0_100 << endl;
0777 
0778   TGraphErrors *bbg_zg_prediction = GetZgPrediction("bbg_zg");
0779   TGraphErrors *me_bbg_prediction = GetZgPrediction("me_bbg");
0780   TGraphErrors *bbg_zg_ratio_prediction = GetZgPrediction("bbg_zg_ratio");
0781 
0782   TCanvas *c1 = new TCanvas("Draw_HFJetTruth_CrossSection2zg_check" + s_suffix, "Draw_HFJetTruth_CrossSection2zg_check" + s_suffix, 700, 600);
0783   c1->Divide(1, 1);
0784   int idx = 1;
0785   TPad *p;
0786 
0787   p = (TPad *) c1->cd(idx++);
0788   c1->Update();
0789   //  p->SetGridx(0);
0790   //  p->SetGridy(0);
0791   //  p->SetLogy();
0792 
0793   TH1 *g_pp = CrossSection2zg_HistogramMaker(h_b, 1., dy, pt_max, pp_lumi * pp_eff * pp_purity, "bbg_zg");
0794   TH1 *g_AA = CrossSection2zg_HistogramMaker(h_b, b_jet_RAA, dy, pt_max, AuAu_eq_lumi_C0_10 * AuAu_eff * AuAu_purity, "me_bbg");
0795 
0796   g_pp->SetLineColor(kRed);
0797   g_AA->SetLineColor(kBlue);
0798 
0799   g_pp->Draw();
0800   g_AA->Draw("same");
0801   SaveCanvas(c1, infile + "_" + TString(c1->GetName()), kTRUE);
0802   //
0803   //
0804   c1 = new TCanvas("Draw_HFJetTruth_CrossSection2zg" + s_suffix, "Draw_HFJetTruth_CrossSection2RAA" + s_suffix, 1100, 800);
0805   c1->Divide(1, 1);
0806   idx = 1;
0807 
0808   p = (TPad *) c1->cd(idx++);
0809   c1->Update();
0810 
0811   p->DrawFrame(0.1, 0, .5, 13)
0812       ->SetTitle(";Jet z_{g};1/N_{jet} dN/dz_{g}");
0813 
0814   //  me_bbg_prediction->SetFillColor(kBlack - 4);
0815   //  me_bbg_prediction->Draw("pe");
0816 
0817   g_AA->SetLineColor(kBlack);
0818   g_AA->SetMarkerColor(kBlack);
0819   g_AA->SetLineWidth(4);
0820   g_AA->SetMarkerStyle(kFullCircle);
0821   g_AA->SetMarkerSize(3);
0822 
0823   g_AA->Draw("pe same");
0824   g_AA->Print();
0825 
0826   g_pp->SetLineColor(kRed + 2);
0827   g_pp->SetMarkerColor(kRed + 2);
0828   g_pp->SetLineWidth(4);
0829   g_pp->SetMarkerStyle(kFullCircle);
0830   g_pp->SetMarkerSize(3);
0831 
0832   g_pp->Draw("pe same");
0833   g_pp->Print();
0834 
0835   TLegend *leg = new TLegend(.0, .77, .85, .91);
0836   leg->SetFillStyle(0);
0837   leg->AddEntry("", "#it{#bf{sPHENIX}} Projection, Year 1-3", "");
0838   //  leg->AddEntry("", Form("PYTHIA-8 #it{b}-jet, Anti-k_{T} R=0.4, |#eta|<%.2f, CTEQ6L", dy / 2), "");
0839   leg->AddEntry("", Form("#it{b}-jet Anti-k_{T} R=0.4, 15<p_{T}<%.0f GeV, z_{cut}>0.1, #beta=0", pt_max), "");
0840   leg->Draw();
0841 
0842   TLegend *leg2 = new TLegend(.2, .6, .85, .75);
0843   leg2->AddEntry(g_pp, Form("#it{p}+#it{p}: %.0fpb^{-1} samp., %.0f%% Eff., %.0f%% Pur.", pp_lumi, pp_eff * 100, pp_purity * 100), "pe");
0844   leg2->AddEntry(g_AA, Form("Au+Au: %.0fnb^{-1} rec., %.0f%% Eff., %.0f%% Pur.", AuAu_MB_Evt / 6.8252 / 1e9, AuAu_eff * 100, AuAu_purity * 100), "pe");
0845   leg2->SetFillStyle(0);
0846   leg2->Draw();
0847 
0848   SaveCanvas(c1, infile + "_" + TString(c1->GetName()), kTRUE);
0849 
0850   c1 = new TCanvas("Draw_HFJetTruth_CrossSection2zg_ratio" + s_suffix, "Draw_HFJetTruth_CrossSection2zg_ratio" + s_suffix, 1100, 800);
0851   c1->Divide(1, 1);
0852   idx = 1;
0853 
0854   p = (TPad *) c1->cd(idx++);
0855   c1->Update();
0856 
0857   p->DrawFrame(0.1, .55, .5, 1.5)
0858       ->SetTitle(";Jet z_{g};P_{AA}(z_{g})/P_{pp}(z_{g})");
0859 
0860   bbg_zg_ratio_prediction->SetFillColor(kBlue - 5);
0861   bbg_zg_ratio_prediction->Draw("3");
0862 
0863   TH1 *g_AA_ratio = (TH1 *) g_AA->Clone("g_AA_ratio");
0864   g_AA_ratio->Divide(g_pp);
0865 
0866   for (int i = 0; i <= g_AA_ratio->GetNbinsX(); ++i)
0867   {
0868     g_AA_ratio->SetBinContent(i, 1.);  // set projection to center onto 1.0
0869   }
0870 
0871   g_AA_ratio->Draw("pe same X0");
0872   g_AA_ratio->Print();
0873 
0874   leg = new TLegend(.0, .67, .85, .91);
0875   leg->SetFillStyle(0);
0876   leg->AddEntry("", "#it{#bf{sPHENIX}} Projection, Year 1-3", "");
0877   //  leg->AddEntry("", Form("PYTHIA-8 #it{b}-jet, Anti-k_{T} R=0.4, |#eta|<%.2f, CTEQ6L", dy / 2), "");
0878   leg->AddEntry("", Form("#it{b}-jet Anti-k_{T} R=0.4, 15<p_{T}<%.0f GeV, z_{cut}>0.1, #beta=0", pt_max), "");
0879   leg->AddEntry("", Form("#it{p}+#it{p}: %.0fpb^{-1} samp., %.0f%% Eff., %.0f%% Pur.", pp_lumi, pp_eff * 100, pp_purity * 100), "");
0880   leg->AddEntry("", Form("Au+Au: %.0fnb^{-1} rec., %.0f%% Eff., %.0f%% Pur.", AuAu_MB_Evt / 6.8252 / 1e9, AuAu_eff * 100, AuAu_purity * 100), "");
0881   leg->Draw();
0882 
0883   leg2 = new TLegend(.2, .2, .75, .3);
0884   leg2->AddEntry(bbg_zg_ratio_prediction, Form("Li, Vitev, PLB 793 (2019)"), "f");
0885   leg2->AddEntry("", Form("b#rightarrowbg, g=2.0#pm0.1 MLL"), "");
0886   leg2->SetFillStyle(0);
0887   leg2->Draw();
0888 
0889   SaveCanvas(c1, infile + "_" + TString(c1->GetName()), kTRUE);
0890 }
0891 
0892 TGraph *CrossSection2v2Uncert(const TH1F *h_cross,
0893                               const double suppression,
0894                               const double deta,
0895                               const double pp_quiv_int_lum,
0896                               const double ep_resolution = 1,
0897                               const double pt_shift = 0)
0898 {
0899   assert(h_cross);
0900   TH1 *
0901       h_ratio = (TH1 *)
0902                     h_cross->Clone(TString(h_cross->GetName()) + Form("_copyv2%d", rand()));
0903 
0904   //convert to count per bin
0905   h_ratio->Scale(deta * h_ratio->GetXaxis()->GetBinWidth(0) * pp_quiv_int_lum * suppression);
0906   h_ratio->Rebin(5);
0907 
0908   vector<double> pts;
0909   vector<double> v2s;
0910   vector<double> v2es;
0911 
0912   for (int i = 1; i <= h_ratio->GetNbinsX(); ++i)
0913   {
0914     const double yield = h_ratio->GetBinContent(i);
0915 
0916     if (yield > 100)
0917     {
0918       h_ratio->SetBinContent(i, 0);
0919 
0920       h_ratio->SetBinError(i, 1. / sqrt(2 * yield) / ep_resolution);  // err(v2) = 1/ (sqrt(2) *Significance * Resolution)
0921 
0922       pts.push_back(h_ratio->GetBinCenter(i) + pt_shift);
0923       v2s.push_back(h_ratio->GetBinContent(i));
0924       v2es.push_back(h_ratio->GetBinError(i));
0925     }
0926     else
0927     {
0928       h_ratio->SetBinContent(i, 0);
0929 
0930       h_ratio->SetBinError(i, 0);
0931     }
0932   }
0933 
0934   h_ratio->GetYaxis()->SetTitle("v2 and uncertainty");
0935 
0936   TGraph *gr = new TGraphErrors(pts.size(), &pts[0], &v2s[0], 0, &v2es[0]);
0937   gr->SetName(TString("ge_") + h_ratio->GetName());
0938 
0939   gr->SetLineWidth(3);
0940   gr->SetMarkerStyle(kFullCircle);
0941   gr->SetMarkerSize(2);
0942 
0943   return gr;
0944 }
0945 
0946 void CrossSection2v2(const TString infile, const bool use_AA_jet_trigger = true, const double ep_resolution = 0.5, const double dy = .7 * 2)
0947 {
0948   TFile *f = TFile::Open(infile + "Draw_HFJetTruth_DrawCrossSection.root");
0949   assert(f);
0950 
0951   TH1F *hall = (TH1F *) f->GetObjectChecked("hall", "TH1F");
0952   assert(hall);
0953   TH1F *h_b = (TH1F *) f->GetObjectChecked("h_b", "TH1F");
0954   assert(h_b);
0955 
0956   const double b_jet_RAA = 0.6;
0957 
0958   const double pp_eff = 0.6;
0959   const double pp_purity = 0.4;
0960   const double AuAu_eff = 0.4;
0961   const double AuAu_purity = 0.4;
0962 
0963   ////////////////////////////
0964   // 5-year lumi in [sPH-TRG-000]
0965   ////////////////////////////
0966 
0967   const double pp_lumi = 62;                           // pb^-1 [sPH-TRG-000], rounded up from 197 pb^-1
0968   const double pp_inelastic_crosssec = 42e-3 / 1e-12;  // 42 mb in pb [sPH-TRG-000]
0969 
0970   const double AuAu_MB_Evt = use_AA_jet_trigger ? 550e9 : 142e9;  // [sPH-TRG-000], depending on whether jet trigger applied in AA collisions
0971   const double pAu_MB_Evt = 600e9;                                // [sPH-TRG-000]
0972 
0973   const double AuAu_Ncoll_C0_10 = 960.2;  // [DOI:?10.1103/PhysRevC.87.034911?]
0974   const double AuAu_Ncoll_C0_20 = 770.6;  // [DOI:?10.1103/PhysRevC.91.064904?]
0975   const double AuAu_Ncoll_C10_20 = 603;   // [sPH-HF-2017-001-v1]
0976   const double AuAu_Ncoll_C20_40 = 296;   // [sPH-HF-2017-001-v1]
0977   const double AuAu_Ncoll_C10_40 = (AuAu_Ncoll_C10_20 * 10 + AuAu_Ncoll_C20_40 * 20) / 30;
0978   const double AuAu_Ncoll_C40_60 = 94;   //  [sPH-HF-2017-001-v1]
0979   const double AuAu_Ncoll_C60_92 = 15;   //  [sPH-HF-2017-001-v1]
0980   const double AuAu_Ncoll_C0_100 = 250;  // pb^-1 [sPH-TRG-000]
0981   const double pAu_Ncoll_C0_100 = 4.7;   // pb^-1 [sPH-TRG-000]
0982 
0983   const double AuAu_eq_lumi_C0_10 = AuAu_MB_Evt * 0.1 * AuAu_Ncoll_C0_10 / pp_inelastic_crosssec;  //
0984   const double AuAu_eq_lumi_C0_20 = AuAu_MB_Evt * 0.2 * AuAu_Ncoll_C0_20 / pp_inelastic_crosssec;  //
0985   const double AuAu_eq_lumi_C0_100 = AuAu_MB_Evt * 1 * AuAu_Ncoll_C0_100 / pp_inelastic_crosssec;  //
0986 
0987   const double AuAu_eq_lumi_C10_20 = AuAu_MB_Evt * .1 * AuAu_Ncoll_C10_20 / pp_inelastic_crosssec;          //
0988   const double AuAu_eq_lumi_C20_40 = AuAu_MB_Evt * .2 * AuAu_Ncoll_C20_40 / pp_inelastic_crosssec;          //
0989   const double AuAu_eq_lumi_C10_40 = AuAu_MB_Evt * .3 * AuAu_Ncoll_C10_40 / pp_inelastic_crosssec;          //
0990   const double AuAu_eq_lumi_C40_60 = AuAu_MB_Evt * .2 * AuAu_Ncoll_C40_60 / pp_inelastic_crosssec;          //
0991   const double AuAu_eq_lumi_C60_92 = AuAu_MB_Evt * (.92 - .6) * AuAu_Ncoll_C60_92 / pp_inelastic_crosssec;  //
0992 
0993   const double pAu_eq_lumi_C0_100 = pAu_MB_Evt * 1 * pAu_Ncoll_C0_100 / pp_inelastic_crosssec;  //
0994 
0995   ;
0996 
0997   TString s_suffix(use_AA_jet_trigger ? "_AAJetTriggered" : "_3yr");
0998   s_suffix += Form("_EPR%.1f", ep_resolution);
0999   s_suffix += Form("_deta%.2f", dy / 2);
1000 
1001   cout << "CrossSection2v2 integrated luminosity assumptions in pb^-1: " << endl;
1002   cout << "\t"
1003        << "pp_lumi = " << pp_lumi << endl;
1004   cout << "\t"
1005        << "AuAu_eq_lumi_C0_10 = " << AuAu_eq_lumi_C0_10 << endl;
1006   cout << "\t"
1007        << "AuAu_eq_lumi_C0_20 = " << AuAu_eq_lumi_C0_20 << endl;
1008   cout << "\t"
1009        << "AuAu_eq_lumi_C0_100 = " << AuAu_eq_lumi_C0_100 << endl;
1010   cout << "\t"
1011        << "AuAu_eq_lumi_C10_20 = " << AuAu_eq_lumi_C10_20 << endl;
1012   cout << "\t"
1013        << "AuAu_eq_lumi_C20_40 = " << AuAu_eq_lumi_C20_40 << endl;
1014   cout << "\t"
1015        << "AuAu_eq_lumi_C10_40 = " << AuAu_eq_lumi_C10_40 << endl;
1016   cout << "\t"
1017        << "AuAu_eq_lumi_C40_60 = " << AuAu_eq_lumi_C40_60 << endl;
1018   cout << "\t"
1019        << "AuAu_eq_lumi_C60_92 = " << AuAu_eq_lumi_C60_92 << endl;
1020   cout << "\t"
1021        << "pAu_eq_lumi_C0_100 = " << pAu_eq_lumi_C0_100 << endl;
1022 
1023   TGraph *g_AA_C0_10 = CrossSection2v2Uncert(h_b, b_jet_RAA, dy, AuAu_eq_lumi_C0_10 * AuAu_eff * AuAu_purity, ep_resolution, -1 * .7);
1024   TGraph *g_AA_C0_20 = CrossSection2v2Uncert(h_b, b_jet_RAA, dy, AuAu_eq_lumi_C0_20 * AuAu_eff * AuAu_purity, ep_resolution, 0);
1025   TGraph *g_AA_C10_20 = CrossSection2v2Uncert(h_b, b_jet_RAA, dy, AuAu_eq_lumi_C10_20 * AuAu_eff * AuAu_purity, ep_resolution, 0);
1026   TGraph *g_AA_C20_40 = CrossSection2v2Uncert(h_b, b_jet_RAA, dy, AuAu_eq_lumi_C20_40 * AuAu_eff * AuAu_purity, ep_resolution, 1 * .7);
1027   TGraph *g_AA_C10_40 = CrossSection2v2Uncert(h_b, b_jet_RAA, dy, AuAu_eq_lumi_C10_40 * AuAu_eff * AuAu_purity, ep_resolution, 1 * .7);
1028   TGraph *g_AA_C40_60 = CrossSection2v2Uncert(h_b, b_jet_RAA, dy, AuAu_eq_lumi_C40_60 * AuAu_eff * AuAu_purity, ep_resolution, 2 * .7);
1029   TGraph *g_AA_C60_92 = CrossSection2v2Uncert(h_b, b_jet_RAA, dy, AuAu_eq_lumi_C60_92 * AuAu_eff * AuAu_purity, ep_resolution, 3 * .7);
1030   //
1031   g_AA_C0_10->SetLineColor(kBlue + 3);
1032   g_AA_C10_20->SetLineColor(kAzure + 3);
1033   g_AA_C20_40->SetLineColor(kTeal + 3);
1034   g_AA_C10_40->SetLineColor(kTeal + 3);
1035   g_AA_C40_60->SetLineColor(kSpring + 3);
1036 
1037   g_AA_C0_10->SetMarkerColor(kBlue + 3);
1038   g_AA_C10_20->SetMarkerColor(kAzure + 3);
1039   g_AA_C20_40->SetMarkerColor(kTeal + 3);
1040   g_AA_C10_40->SetMarkerColor(kTeal + 3);
1041   g_AA_C40_60->SetMarkerColor(kSpring + 3);
1042 
1043   g_AA_C0_10->SetMarkerStyle(kFullCircle);
1044   g_AA_C10_20->SetMarkerStyle(kFullSquare);
1045   g_AA_C20_40->SetMarkerStyle(kFullDiamond);
1046   g_AA_C10_40->SetMarkerStyle(kFullDiamond);
1047   g_AA_C40_60->SetMarkerStyle(kFullCross);
1048 
1049   g_AA_C0_10->SetLineWidth(4);
1050   g_AA_C0_10->SetMarkerSize(3);
1051   g_AA_C10_40->SetLineWidth(4);
1052   g_AA_C10_40->SetMarkerSize(3);
1053   g_AA_C40_60->SetLineWidth(4);
1054   g_AA_C40_60->SetMarkerSize(3);
1055   //
1056   TCanvas *c1 = new TCanvas("Draw_HFJetTruth_CrossSection2v2" + s_suffix, "Draw_HFJetTruth_CrossSection2v2" + s_suffix, 1100, 800);
1057   c1->Divide(1, 1);
1058   int idx = 1;
1059   TPad *p;
1060 
1061   p = (TPad *) c1->cd(idx++);
1062   c1->Update();
1063 
1064   p->DrawFrame(15, -.1, 40, .3)
1065       ->SetTitle(";p_{T} [GeV];v_{2}");
1066   //
1067   g_AA_C0_10->Draw("pe");
1068   //  g_AA_C10_20->Draw("pe");
1069   //  g_AA_C20_40->Draw("pe");
1070   g_AA_C10_40->Draw("pe");
1071   //  g_AA_C40_60->Draw("pe");
1072 
1073   //  g_AA_C20_40->Draw("same");
1074   //
1075   //  ge_RAA->SetLineWidth(3);
1076   //  ge_RAA->SetMarkerStyle(kFullCircle);
1077   //  ge_RAA->SetMarkerSize(2);
1078   //
1079   //  ge_RAA->Draw("pe");
1080   //  ge_RAA->Print();
1081   //
1082   TLegend *leg = new TLegend(.0, .67, .85, .91);
1083   leg->SetFillStyle(0);
1084   leg->AddEntry("", "#it{#bf{sPHENIX}} Projection, Au+Au, Year 1-3", "");
1085   leg->AddEntry("", Form("#it{b}-jet, Anti-k_{T} R=0.4, |#eta|<%.2f,Res(#Psi_{2})=%.1f", dy / 2, ep_resolution), "");
1086   leg->AddEntry("", Form("%.0fnb^{-1} rec. Au+Au, %.0f%% Eff., %.0f%% Pur.", AuAu_MB_Evt / 6.8252 / 1e9, AuAu_eff * 100, AuAu_purity * 100), "");
1087   leg->Draw();
1088   //
1089   TLegend *leg2 = new TLegend(.19, .5, 1, .67);
1090   //  leg2->SetHeader(Form("#it{b}-jet v_{2} Projection, #it{R}_{AA, #it{b}-jet}=%.1f, Res(#Psi_{2})=%.1f", b_jet_RAA, ep_resolution));
1091   leg2->AddEntry(g_AA_C0_10, "0-10% Au+Au", "pl");
1092   //    leg2->AddEntry(g_AA_C10_20, "Au+Au 10-20%C", "pl");
1093   //    leg2->AddEntry(g_AA_C20_40, "Au+Au 20-40%C", "pl");
1094   leg2->AddEntry(g_AA_C10_40, "10-40% Au+Au ", "pl");
1095   //    leg2->AddEntry(g_AA_C40_60, "Au+Au 40-60%C", "pl");
1096   leg2->SetFillStyle(0);
1097   leg2->Draw();
1098 
1099   SaveCanvas(c1, infile + "_" + TString(c1->GetName()), kTRUE);
1100 }
1101 
1102 void Convert2CrossSection(TH1 *h, const double int_lumi, const double dy)
1103 {
1104   h->Sumw2();
1105 
1106   for (int i = 1; i <= h->GetXaxis()->GetNbins(); ++i)
1107   {
1108     const double pT1 = h->GetXaxis()->GetBinLowEdge(i);
1109     const double pT2 = h->GetXaxis()->GetBinUpEdge(i);
1110 
1111     //      const double dpT2 =  pT2*pT2 - pT1*pT1;
1112     const double dpT = pT2 - pT1;
1113 
1114     //      const double count2sigma = 1./dpT2/dy/int_lumi;
1115     const double count2sigma = 1. / dpT / dy / int_lumi;
1116 
1117     h->SetBinContent(i, h->GetBinContent(i) * count2sigma);
1118     h->SetBinError(i, h->GetBinError(i) * count2sigma);
1119   }
1120 
1121   //  h->GetYaxis()->SetTitle("#sigma/(dp_{T}^{2} dy) (pb/(GeV/c)^{2})");
1122 
1123   h->SetMarkerStyle(kFullCircle);
1124 }
1125 
1126 TGraphAsymmErrors *
1127 GetFONLL_B()
1128 {
1129   //# Job started on: Wed Aug 24 04:28:34 CEST 2016 .
1130   //# FONLL heavy quark hadroproduction cross section, calculated on Wed Aug 24 04:28:42 CEST 2016
1131   //# FONLL version and perturbative order: ## FONLL v1.3.2 fonll [ds/dpt^2dy (pb/GeV^2)]
1132   //# quark = bottom
1133   //# final state = quark
1134   //# ebeam1 = 100, ebeam2 = 100
1135   //# PDF set = CTEQ6.6
1136   //# ptmin = 10
1137   //# ptmax = 30
1138   //# etamin = -.6
1139   //# etamax = .6
1140   //# Uncertainties from scales, masses combined quadratically
1141   //# cross section is ds/dpt (pb/GeV)
1142   //# pt      central      min       max       min_sc     max_sc     min_mass   max_mass
1143   //  10.0000 4.1645e+03 2.9134e+03 6.0916e+03 2.9634e+03 6.0553e+03 3.8146e+03 4.5365e+03
1144   //  11.0526 2.4830e+03 1.7438e+03 3.6234e+03 1.7667e+03 3.6072e+03 2.3005e+03 2.6743e+03
1145   //  12.1053 1.5046e+03 1.0602e+03 2.1900e+03 1.0710e+03 2.1826e+03 1.4071e+03 1.6056e+03
1146   //  13.1579 9.2707e+02 6.5506e+02 1.3454e+03 6.6029e+02 1.3418e+03 8.7394e+02 9.8137e+02
1147   //  14.2105 5.8027e+02 4.1098e+02 8.3990e+02 4.1358e+02 8.3817e+02 5.5067e+02 6.1015e+02
1148   //  15.2632 3.6865e+02 2.6159e+02 5.3224e+02 2.6290e+02 5.3139e+02 3.5193e+02 3.8532e+02
1149   //  16.3158 2.3744e+02 1.6874e+02 3.4198e+02 1.6941e+02 3.4156e+02 2.2785e+02 2.4686e+02
1150   //  17.3684 1.5484e+02 1.1016e+02 2.2257e+02 1.1051e+02 2.2236e+02 1.4930e+02 1.6021e+02
1151   //  18.4211 1.0215e+02 7.2725e+01 1.4657e+02 7.2904e+01 1.4646e+02 9.8911e+01 1.0522e+02
1152   //  19.4737 6.8101e+01 4.8502e+01 9.7587e+01 4.8594e+01 9.7535e+01 6.6205e+01 6.9866e+01
1153   //  20.5263 4.5818e+01 3.2631e+01 6.5588e+01 3.2678e+01 6.5562e+01 4.4707e+01 4.6827e+01
1154   //  21.5789 3.1111e+01 2.2148e+01 4.4506e+01 2.2171e+01 4.4494e+01 3.0461e+01 3.1684e+01
1155   //  22.6316 2.1287e+01 1.5142e+01 3.0444e+01 1.5153e+01 3.0439e+01 2.0909e+01 2.1607e+01
1156   //  23.6842 1.4663e+01 1.0419e+01 2.0974e+01 1.0424e+01 2.0972e+01 1.4446e+01 1.4838e+01
1157   //  24.7368 1.0175e+01 7.2188e+00 1.4562e+01 7.2214e+00 1.4562e+01 1.0053e+01 1.0266e+01
1158   //  25.7895 7.1026e+00 5.0296e+00 1.0176e+01 5.0307e+00 1.0175e+01 7.0360e+00 7.1470e+00
1159   //  26.8421 4.9796e+00 3.5183e+00 7.1439e+00 3.5187e+00 7.1438e+00 4.9454e+00 4.9976e+00
1160   //  27.8947 3.5103e+00 2.4738e+00 5.0451e+00 2.4739e+00 5.0451e+00 3.4945e+00 3.5143e+00
1161   //  28.9474 2.4892e+00 1.7492e+00 3.5856e+00 1.7493e+00 3.5856e+00 2.4835e+00 2.4892e+00
1162   //  30.0000 1.7693e+00 1.2394e+00 2.5556e+00 1.2394e+00 2.5556e+00 1.7633e+00 1.7693e+00
1163 
1164   const double deta = 0.6 * 2;
1165 
1166   const double pts[] =
1167       {
1168 
1169           10.0000, 11.0526, 12.1053, 13.1579, 14.2105, 15.2632, 16.3158, 17.3684,
1170           18.4211, 19.4737, 20.5263, 21.5789, 22.6316, 23.6842, 24.7368, 25.7895,
1171           26.8421, 27.8947, 28.9474, 30.0000
1172 
1173       };
1174 
1175   const double centrals[] =
1176       {4.1645e+03, 2.4830e+03, 1.5046e+03, 9.2707e+02, 5.8027e+02, 3.6865e+02,
1177        2.3744e+02, 1.5484e+02, 1.0215e+02, 6.8101e+01, 4.5818e+01, 3.1111e+01,
1178        2.1287e+01, 1.4663e+01, 1.0175e+01, 7.1026e+00, 4.9796e+00, 3.5103e+00,
1179        2.4892e+00, 1.7693e+00};
1180   const double min[] =
1181       {4.1645e+03 - 2.9134e+03, 2.4830e+03 - 1.7438e+03, 1.5046e+03 - 1.0602e+03,
1182        9.2707e+02 - 6.5506e+02, 5.8027e+02 - 4.1098e+02, 3.6865e+02 - 2.6159e+02, 2.3744e+02 - 1.6874e+02, 1.5484e+02 - 1.1016e+02,
1183        1.0215e+02 - 7.2725e+01, 6.8101e+01 - 4.8502e+01, 4.5818e+01 - 3.2631e+01, 3.1111e+01 - 2.2148e+01, 2.1287e+01 - 1.5142e+01,
1184        1.4663e+01 - 1.0419e+01, 1.0175e+01 - 7.2188e+00, 7.1026e+00 - 5.0296e+00, 4.9796e+00 - 3.5183e+00, 3.5103e+00 - 2.4738e+00,
1185        2.4892e+00 - 1.7492e+00, 1.7693e+00 - 1.2394e+00};
1186 
1187   const double max[] =
1188       {6.0916e+03 - 4.1645e+03, 3.6234e+03 - 2.4830e+03, 2.1900e+03 - 1.5046e+03,
1189        1.3454e+03 - 9.2707e+02, 8.3990e+02 - 5.8027e+02, 5.3224e+02 - 3.6865e+02, 3.4198e+02 - 2.3744e+02, 2.2257e+02 - 1.5484e+02,
1190        1.4657e+02 - 1.0215e+02, 9.7587e+01 - 6.8101e+01, 6.5588e+01 - 4.5818e+01, 4.4506e+01 - 3.1111e+01, 3.0444e+01 - 2.1287e+01,
1191        2.0974e+01 - 1.4663e+01, 1.4562e+01 - 1.0175e+01, 1.0176e+01 - 7.1026e+00, 7.1439e+00 - 4.9796e+00, 5.0451e+00 - 3.5103e+00,
1192        3.5856e+00 - 2.4892e+00, 2.5556e+00 - 1.7693e+00};
1193 
1194   TGraphAsymmErrors *gr = new TGraphAsymmErrors(20, pts, centrals, 0, 0, min,
1195                                                 max);
1196 
1197   for (int i = 0; i < gr->GetN(); ++i)
1198   {
1199     (gr->GetY())[i] /= deta;
1200     (gr->GetEYhigh())[i] /= deta;
1201     (gr->GetEYlow())[i] /= deta;
1202   }
1203 
1204   return gr;
1205 }
1206 
1207 TGraphAsymmErrors *
1208 GetFONLL_C()
1209 {
1210   //# Job started on: Thu Aug 25 05:08:57 CEST 2016 .
1211   //# FONLL heavy quark hadroproduction cross section, calculated on Thu Aug 25 05:09:06 CEST 2016
1212   //# FONLL version and perturbative order: ## FONLL v1.3.2 fonll [ds/dpt^2dy (pb/GeV^2)]
1213   //# quark = charm
1214   //# final state = quark
1215   //# ebeam1 = 100, ebeam2 = 100
1216   //# PDF set = CTEQ6.6
1217   //# ptmin = 10
1218   //# ptmax = 30
1219   //# etamin = -.6
1220   //# etamax = .6
1221   //# Uncertainties from scales, masses combined quadratically
1222   //# cross section is ds/dpt (pb/GeV)
1223   //# pt      central      min       max       min_sc     max_sc     min_mass   max_mass
1224   //  10.0000 9.7368e+03 6.5051e+03 1.5528e+04 6.5641e+03 1.5522e+04 9.1225e+03 1.0006e+04
1225   //  11.0526 4.9069e+03 3.2787e+03 7.7308e+03 3.3256e+03 7.7226e+03 4.5190e+03 5.1221e+03
1226   //  12.1053 2.5743e+03 1.7187e+03 4.0147e+03 1.7525e+03 4.0069e+03 2.3363e+03 2.7237e+03
1227   //  13.1579 1.3990e+03 9.3274e+02 2.1641e+03 9.5597e+02 2.1578e+03 1.2537e+03 1.4976e+03
1228   //  14.2105 7.8344e+02 5.2138e+02 1.2037e+03 5.3701e+02 1.1989e+03 6.9430e+02 8.4730e+02
1229   //  15.2632 4.5087e+02 2.9938e+02 6.8902e+02 3.0978e+02 6.8543e+02 3.9570e+02 4.9204e+02
1230   //  16.3158 2.6550e+02 1.7584e+02 4.0398e+02 1.8275e+02 4.0140e+02 2.3098e+02 2.9210e+02
1231   //  17.3684 1.5978e+02 1.0551e+02 2.4232e+02 1.1011e+02 2.4049e+02 1.3793e+02 1.7705e+02
1232   //  18.4211 9.7903e+01 6.4459e+01 1.4811e+02 6.7524e+01 1.4682e+02 8.3915e+01 1.0919e+02
1233   //  19.4737 6.0993e+01 4.0030e+01 9.2103e+01 4.2085e+01 9.1200e+01 5.1941e+01 6.8435e+01
1234   //  20.5263 3.8605e+01 2.5248e+01 5.8239e+01 2.6631e+01 5.7606e+01 3.2685e+01 4.3550e+01
1235   //  21.5789 2.4727e+01 1.6108e+01 3.7286e+01 1.7044e+01 3.6842e+01 2.0821e+01 2.8038e+01
1236   //  22.6316 1.6036e+01 1.0402e+01 2.4185e+01 1.1040e+01 2.3872e+01 1.3433e+01 1.8271e+01
1237   //  23.6842 1.0532e+01 6.8006e+00 1.5898e+01 7.2389e+00 1.5678e+01 8.7774e+00 1.2053e+01
1238   //  24.7368 6.9729e+00 4.4803e+00 1.0540e+01 4.7838e+00 1.0385e+01 5.7809e+00 8.0136e+00
1239   //  25.7895 4.6560e+00 2.9760e+00 7.0510e+00 3.1876e+00 6.9411e+00 3.8397e+00 5.3729e+00
1240   //  26.8421 3.1409e+00 1.9970e+00 4.7681e+00 2.1454e+00 4.6902e+00 2.5775e+00 3.6384e+00
1241   //  27.8947 2.1346e+00 1.3500e+00 3.2500e+00 1.4543e+00 3.1946e+00 1.7438e+00 2.4816e+00
1242   //  28.9474 1.4557e+00 9.1535e-01 2.2239e+00 9.8868e-01 2.1845e+00 1.1839e+00 1.6985e+00
1243   //  30.0000 9.9826e-01 6.2400e-01 1.5310e+00 6.7571e-01 1.5029e+00 8.0844e-01 1.1690e+00
1244 
1245   const double deta = 0.6 * 2;
1246 
1247   const double pts[] =
1248       {
1249           //         pt
1250           10.0000, 11.0526, 12.1053, 13.1579, 14.2105, 15.2632, 16.3158, 17.3684,
1251           18.4211, 19.4737, 20.5263, 21.5789, 22.6316, 23.6842, 24.7368, 25.7895,
1252           26.8421, 27.8947, 28.9474, 30.0000
1253 
1254       };
1255 
1256   const double centrals[] =
1257       {
1258 
1259           //            central
1260           9.7368e+03, 4.9069e+03, 2.5743e+03, 1.3990e+03, 7.8344e+02, 4.5087e+02,
1261           2.6550e+02, 1.5978e+02, 9.7903e+01, 6.0993e+01, 3.8605e+01, 2.4727e+01,
1262           1.6036e+01, 1.0532e+01, 6.9729e+00, 4.6560e+00, 3.1409e+00, 2.1346e+00,
1263           1.4557e+00, 9.9826e-01
1264 
1265       };
1266   const double min[] =
1267       {
1268 
1269           //       central           min
1270           9.7368e+03 - 6.5051e+03, 4.9069e+03 - 3.2787e+03, 2.5743e+03 - 1.7187e+03, 1.3990e+03 - 9.3274e+02, 7.8344e+02 - 5.2138e+02,
1271           4.5087e+02 - 2.9938e+02, 2.6550e+02 - 1.7584e+02, 1.5978e+02 - 1.0551e+02, 9.7903e+01 - 6.4459e+01, 6.0993e+01 - 4.0030e+01,
1272           3.8605e+01 - 2.5248e+01, 2.4727e+01 - 1.6108e+01, 1.6036e+01 - 1.0402e+01, 1.0532e+01 - 6.8006e+00, 6.9729e+00 - 4.4803e+00,
1273           4.6560e+00 - 2.9760e+00, 3.1409e+00 - 1.9970e+00, 2.1346e+00 - 1.3500e+00, 1.4557e+00 - 9.1535e-01, 9.9826e-01 - 6.2400e-01
1274 
1275       };
1276 
1277   const double max[] =
1278       {
1279 
1280           //            max             central
1281           1.5528e+04 - 9.7368e+03, 7.7308e+03 - 4.9069e+03, 4.0147e+03 - 2.5743e+03, 2.1641e+03 - 1.3990e+03, 1.2037e+03 - 7.8344e+02,
1282           6.8902e+02 - 4.5087e+02, 4.0398e+02 - 2.6550e+02, 2.4232e+02 - 1.5978e+02, 1.4811e+02 - 9.7903e+01, 9.2103e+01 - 6.0993e+01,
1283           5.8239e+01 - 3.8605e+01, 3.7286e+01 - 2.4727e+01, 2.4185e+01 - 1.6036e+01, 1.5898e+01 - 1.0532e+01, 1.0540e+01 - 6.9729e+00,
1284           7.0510e+00 - 4.6560e+00, 4.7681e+00 - 3.1409e+00, 3.2500e+00 - 2.1346e+00, 2.2239e+00 - 1.4557e+00, 1.5310e+00 - 9.9826e-01
1285 
1286       };
1287 
1288   TGraphAsymmErrors *gr = new TGraphAsymmErrors(20, pts, centrals, 0, 0, min,
1289                                                 max);
1290 
1291   for (int i = 0; i < gr->GetN(); ++i)
1292   {
1293     (gr->GetY())[i] /= deta;
1294     (gr->GetEYhigh())[i] /= deta;
1295     (gr->GetEYlow())[i] /= deta;
1296   }
1297 
1298   return gr;
1299 }
1300 
1301 TGraph *pQCDModel_HuangKangVitev(const double g)
1302 {
1303   // arXiv:1306.0909v2 [hep-ph]
1304   // Preliminary for sPHENIX energy
1305 
1306   //  R = 0.4;
1307   if (g == 2.0)
1308   {
1309     //  g=2.0
1310     //  10.130739 0.8087402
1311     //  13.8755   0.7102165
1312     //  17.915709 0.64113617
1313     //  19.950817 0.62270004
1314     //  21.838287 0.59506375
1315     //  24.758118 0.5618901
1316     //  28.149998 0.53331053
1317     //  30.303171 0.5194738
1318     //  35.701473 0.52399224
1319     //  40.21508  0.54508865
1320     //  49.182514 0.53758913
1321     //  54.492188 0.5338267
1322     //  59.890347 0.52914274
1323 
1324     const double pT[] =
1325         {
1326             10.130739,
1327             13.8755,
1328             17.915709,
1329             19.950817,
1330             21.838287,
1331             24.758118,
1332             28.149998,
1333             30.303171,
1334             35.701473,
1335             40.21508,
1336             49.182514,
1337             54.492188,
1338             59.890347};
1339 
1340     const double RAA[] =
1341         {
1342             0.8087402,
1343             0.7102165,
1344             0.64113617,
1345             0.62270004,
1346             0.59506375,
1347             0.5618901,
1348             0.53331053,
1349             0.5194738,
1350             0.52399224,
1351             0.54508865,
1352             0.53758913,
1353             0.5338267,
1354             0.52914274};
1355 
1356     return new TGraph(13, pT, RAA);
1357   }
1358   else if (g == 2.2)
1359   {
1360     //  g=2.2
1361     //  10.040924 0.72499925
1362     //  11.750852 0.6623964
1363     //  15.820038 0.56018674
1364     //  18.946156 0.51412654
1365     //  19.860464 0.50491005
1366     //  21.010588 0.484647
1367     //  24.343283 0.44410512
1368     //  29.888279 0.39800778
1369     //  35.876522 0.4006767
1370     //  40.154125 0.42085648
1371     //  47.941647 0.41521555
1372     //  55.64066  0.40865576
1373     //  59.91793  0.4076699
1374 
1375     const double pT[] =
1376         {
1377             10.040924,
1378             11.750852,
1379             15.820038,
1380             18.946156,
1381             19.860464,
1382             21.010588,
1383             24.343283,
1384             29.888279,
1385             35.876522,
1386             40.154125,
1387             47.941647,
1388             55.64066,
1389             59.91793};
1390 
1391     const double RAA[] =
1392         {
1393             0.72499925,
1394             0.6623964,
1395             0.56018674,
1396             0.51412654,
1397             0.50491005,
1398             0.484647,
1399             0.44410512,
1400             0.39800778,
1401             0.4006767,
1402             0.42085648,
1403             0.41521555,
1404             0.40865576,
1405             0.4076699};
1406 
1407     return new TGraph(13, pT, RAA);
1408   }
1409   else
1410   {
1411     assert(0);
1412   }
1413 
1414   return nullptr;
1415 }
1416 
1417 TGraph *LIDO_Ke()
1418 {
1419   //# pT [GeV]  pT-lower  pT-upper  Raa[mu=2*pi*T]  Raa[mu=1.5*pi*T]
1420   //9.0   8.0   10.0  0.704 0.586
1421   //12.0  10.0  14.0  0.632 0.552
1422   //16.0  14.0  18.0  0.544 0.469
1423   //21.5  18.0  25.0  0.519 0.368
1424   //30.0  25.0  35.0  0.417 0.341
1425   //40.0  35.0  45.0  0.355 0.282
1426   //52.5  45.0  60.0  0.317 0.186
1427   //70.0  60.0  80.0  0.286 0.143
1428 
1429   //  Once you unzip the LIDO-heavy-jets.tar.gz file, you will find data files for inclusive, D-jet and B-jet Raa for R=0.2, 0.3, and 0.4. Each file's format has been commented at the first line:
1430   //  # pT    pT-lower    pT-upper    Raa[mu=2*pi*T]    Raa[mu=1.5*pi*T]
1431   //  These five columns are: the jet pT, lower and upper bounds of the jet pT bins. And the two Raa values obtained with two choices of the coupling strength parameters in our model.
1432   //
1433   //  The Raa range enclosed by these two choices of coupling parameter only shows the sensitivity of the calculation to the variation of this single parameter in our model. We are working on a more systematic understanding of our theoretical uncertainty recently.
1434   // 2008.07622 [nucl-th]
1435 
1436   const vector<double> pT{
1437       9.0,
1438       12.0,
1439       16.0,
1440       21.5,
1441       30.0,
1442       40.0,
1443       52.5,
1444       70.0};
1445 
1446   const vector<double> Raa_2piT{
1447       0.704,
1448       0.632,
1449       0.544,
1450       0.519,
1451       0.417,
1452       0.355,
1453       0.317,
1454       0.286};
1455 
1456   const vector<double> Raa_1p5piT{
1457       0.586,
1458       0.552,
1459       0.469,
1460       0.368,
1461       0.341,
1462       0.282,
1463       0.186,
1464       0.143};
1465 
1466   vector<double> Raa(pT.size());
1467   vector<double> Raa_Error(pT.size());
1468 
1469   for (int i = 0; i < pT.size(); ++i)
1470   {
1471     Raa[i] = 0.5 * (Raa_2piT[i] + Raa_1p5piT[i]);
1472     Raa_Error[i] = 0.5 * fabs(Raa_2piT[i] - Raa_1p5piT[i]);
1473   }
1474 
1475   return new TGraphErrors(pT.size(), pT.data(), Raa.data(), nullptr, Raa_Error.data());
1476 }
1477 
1478 TGraph *
1479 GetPHENIX_jet()
1480 {
1481   //  PPG184
1482   //  p+p d^2sigma/dpT deta [mb / GeV]
1483   //
1484   //  pT low     pT high     yval        ystat[%] ysyst[%]
1485   //
1486   //  12.2       14.5        0.0001145   1.0      15.6
1487   //  14.5       17.3        3.774e-05   1.3      16.3
1488   //  17.3       20.7        1.263e-05   1.7      15.6
1489   //  20.7       24.7        4.230e-06   2.3      17.3
1490   //  24.7       29.4        1.224e-06   3.6      19.3
1491   //  29.4       35.1        2.962e-07   6.1      21.0
1492   //  35.1       41.9        5.837e-08   8.9      24.3
1493   //  41.9       50.0        8.882e-09   11       32.3
1494 
1495   const double pt_l[] =
1496       {
1497 
1498           //          pT low
1499 
1500           12.2, 14.5, 17.3, 20.7, 24.7, 29.4, 35.1, 41.9
1501 
1502       };
1503 
1504   const double pt_h[] =
1505       {
1506 
1507           //              pT high
1508 
1509           14.5, 17.3, 20.7, 24.7, 29.4, 35.1, 41.9, 50.0
1510 
1511       };
1512   double yval[] =
1513       {
1514 
1515           //          yval
1516 
1517           0.0001145, 3.774e-05, 1.263e-05, 4.230e-06, 1.224e-06, 2.962e-07,
1518           5.837e-08, 8.882e-09
1519 
1520       };
1521   const double ystat[] =
1522       {
1523 
1524           //          ystat[%]
1525 
1526           1.0, 1.3, 1.7, 2.3, 3.6, 6.1, 8.9, 11
1527 
1528       };
1529 
1530   const double ysyst[] =
1531       {
1532 
1533           //            ysyst[%]
1534 
1535           15.6, 16.3, 15.6, 17.3, 19.3, 21.0, 24.3, 32.3
1536 
1537       };
1538 
1539   double pT_c[100] =
1540       {0};
1541   double pT_e[100] =
1542       {0};
1543   double y_e[100] =
1544       {0};
1545 
1546   for (int i = 0; i < 8; ++i)
1547   {
1548     pT_c[i] = 0.5 * (pt_l[i] + pt_h[i]);
1549     pT_e[i] = 0.5 * (pt_h[i] - pt_l[i]);
1550     yval[i] *= 1e9;  // mb -> pb
1551     y_e[i] = yval[i] * sqrt(ystat[i] * ystat[i] + ysyst[i] * ysyst[i]) / 100;
1552   }
1553 
1554   TGraphErrors *gr = new TGraphErrors(8, pT_c, yval, pT_e, y_e);
1555 
1556   return gr;
1557 }
1558 
1559 void DrawCrossSection(double int_lumi, const double dy)
1560 {
1561   TH1 *hall = new TH1F("hall", ";p_{T} [GeV/c]", 100, 0, 100);
1562   TH1 *h_b = new TH1F("h_b", ";p_{T} [GeV/c]", 100, 0, 100);
1563   TH1 *h_bq = new TH1F("h_bq", ";p_{T} [GeV/c]", 100, 0, 100);
1564   TH1 *h_bh = new TH1F("h_bh", ";p_{T} [GeV/c]", 100, 0, 100);
1565   TH1 *h_bh5 = new TH1F("h_bh5", ";p_{T} [GeV/c]", 100, 0, 100);
1566   TH1 *h_ch5 = new TH1F("h_ch5", ";p_{T} [GeV/c]", 100, 0, 100);
1567   TH1 *h_c = new TH1F("h_c", ";p_{T} [GeV/c]", 100, 0, 100);
1568 
1569   T->Draw("AntiKt_Truth_r04.get_pt()>>hall",
1570           "AntiKt_Truth_r04.get_pt()>15 && abs(AntiKt_Truth_r04.get_eta())<0.6",
1571           "goff");
1572 
1573   T->Draw("AntiKt_Truth_r04.get_pt()>>h_b",
1574           "AntiKt_Truth_r04.get_pt()>15 && abs(AntiKt_Truth_r04.get_eta())<0.6 && abs(AntiKt_Truth_r04.get_property(1000))==5",
1575           "goff");
1576   T->Draw(
1577       "AntiKt_Truth_r04.get_pt()* AntiKt_Truth_r04.get_property(1001) >>h_bq",
1578       "AntiKt_Truth_r04.get_pt()>15 && abs(AntiKt_Truth_r04.get_eta())<0.6 && abs(AntiKt_Truth_r04.get_property(1000))==5",
1579       "goff");
1580 
1581   //    T->Draw("AntiKt_Truth_r04.get_pt()>>h_bh",
1582   //        "AntiKt_Truth_r04.get_pt()>15 && abs(AntiKt_Truth_r04.get_eta())<0.6 && AntiKt_Truth_r04.get_property(1010)==5",
1583   //        "goff");
1584 
1585   T->Draw("AntiKt_Truth_r04.get_pt()>>h_bh5",
1586           "AntiKt_Truth_r04.get_pt()>15 && abs(AntiKt_Truth_r04.get_eta())<0.6 && abs(AntiKt_Truth_r04.get_property(1010))==5 &&  AntiKt_Truth_r04.get_property(1011) * AntiKt_Truth_r04.get_pt() > 5",
1587           "goff");
1588   //
1589   T->Draw("AntiKt_Truth_r04.get_pt()>>h_c",
1590           "AntiKt_Truth_r04.get_pt()>15 && abs(AntiKt_Truth_r04.get_eta())<0.6 && abs(AntiKt_Truth_r04.get_property(1000))==4",
1591           "goff");
1592   T->Draw("AntiKt_Truth_r04.get_pt()>>h_ch5",
1593           "AntiKt_Truth_r04.get_pt()>15 && abs(AntiKt_Truth_r04.get_eta())<0.6 && abs(AntiKt_Truth_r04.get_property(1010))==4 &&  AntiKt_Truth_r04.get_property(1011) * AntiKt_Truth_r04.get_pt() > 5",
1594           "goff");
1595 
1596   Convert2CrossSection(hall, int_lumi, dy);
1597   Convert2CrossSection(h_b, int_lumi, dy);
1598   Convert2CrossSection(h_bq, int_lumi, dy);
1599   Convert2CrossSection(h_bh5, int_lumi, dy);
1600   Convert2CrossSection(h_ch5, int_lumi, dy);
1601   Convert2CrossSection(h_c, int_lumi, dy);
1602 
1603   h_b->SetLineColor(kBlue);
1604   h_b->SetMarkerColor(kBlue);
1605 
1606   h_bq->SetLineColor(kBlue + 3);
1607   h_bq->SetMarkerColor(kBlue + 3);
1608 
1609   h_bh5->SetLineColor(kBlue + 3);
1610   h_bh5->SetMarkerColor(kBlue + 3);
1611   h_bh5->SetMarkerStyle(kStar);
1612 
1613   h_c->SetLineColor(kRed);
1614   h_c->SetMarkerColor(kRed);
1615 
1616   h_ch5->SetLineColor(kRed + 3);
1617   h_ch5->SetMarkerColor(kRed + 3);
1618   h_ch5->SetMarkerStyle(kStar);
1619 
1620   TGraphAsymmErrors *gr_fonll_b = GetFONLL_B();
1621   gr_fonll_b->SetFillColor(kBlue - 7);
1622   gr_fonll_b->SetFillStyle(3002);
1623   TGraphAsymmErrors *gr_fonll_c = GetFONLL_C();
1624   gr_fonll_c->SetFillColor(kRed - 7);
1625   gr_fonll_c->SetFillStyle(3003);
1626   TGraph *gr_phenix = GetPHENIX_jet();
1627   gr_phenix->SetLineColor(kGray + 2);
1628   gr_phenix->SetMarkerColor(kBlack);
1629   gr_phenix->SetMarkerStyle(kOpenCross);
1630 
1631   TCanvas *c1 = new TCanvas("Draw_HFJetTruth_DrawCrossSection", "Draw_HFJetTruth_DrawCrossSection", 1000, 860);
1632   //  c1->Divide(2, 2);
1633   int idx = 1;
1634   TPad *p;
1635 
1636   p = (TPad *) c1->cd(idx++);
1637   c1->Update();
1638   p->SetGridx(0);
1639   p->SetGridy(0);
1640   p->SetLogy();
1641 
1642   hall->Draw();
1643 
1644   gr_fonll_b->Draw("3");
1645   gr_fonll_c->Draw("3");
1646   gr_phenix->Draw("pe");
1647 
1648   h_b->Draw("same");
1649   h_bh5->Draw("same");
1650   h_bq->Draw("same");
1651   h_c->Draw("same");
1652   h_ch5->Draw("same");
1653 
1654   hall->GetXaxis()->SetRangeUser(12, 60);
1655   hall->GetYaxis()->SetTitle(
1656       "d^{2}#sigma/(dp_{T}dy), d^{2}#sigma/(dp_{T}d#eta) [pb/(GeV/c)]");
1657 
1658   TLegend *leg = new TLegend(0.45, 0.6, 0.99, 0.95);
1659   leg->SetFillColor(kWhite);
1660   leg->SetFillStyle(1001);
1661   leg->SetHeader("#it{p}+#it{p} collisions @ sPHENIX, #sqrt{s} = 200 GeV, |#eta|<0.6");
1662   leg->AddEntry(hall, "Inclusive jet, Pythia8, Truth, anti-k_{t}, R=0.4",
1663                 "lpe");
1664   leg->AddEntry(h_c, "c-quark jet, Pythia8, Truth, anti-k_{t}, R=0.4", "lpe");
1665   leg->AddEntry(h_b, "b-quark jet, Pythia8, Truth, anti-k_{t}, R=0.4", "lpe");
1666   leg->AddEntry(h_bq, "Leading b-quark in jet, Pythia8, anti-k_{t}, R=0.4", "lpe");
1667   leg->AddEntry(h_bh5,
1668                 "b-hadron jet, Pythia8, Truth, anti-k_{t}, R=0.4, p_{T, b-hadron}>5 GeV/c",
1669                 "lpe");
1670   leg->AddEntry(gr_phenix, "PHENIX inclusive jet, PRL 116, 122301 (2016)", "ple");
1671   leg->AddEntry(gr_fonll_c, "c-quark, FONLL v1.3.2, CTEQ6.6, |y|<0.6", "f");
1672   leg->AddEntry(gr_fonll_b, "b-quark, FONLL v1.3.2, CTEQ6.6, |y|<0.6", "f");
1673   leg->Draw();
1674 
1675   SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
1676 }
1677 
1678 void Draw_HFJetTruth_DrawCrossSection_PR(const TString infile)
1679 {
1680   TFile *f = TFile::Open(infile + "Draw_HFJetTruth_DrawCrossSection.root");
1681   assert(f);
1682 
1683   TH1F *hall = (TH1F *) f->GetObjectChecked("hall", "TH1F");
1684   assert(hall);
1685   TH1F *h_b = (TH1F *) f->GetObjectChecked("h_b", "TH1F");
1686   assert(h_b);
1687 
1688   hall->SetMarkerColor(kBlack);
1689   hall->SetLineColor(kBlack);
1690   hall->SetMarkerStyle(kFullCircle);
1691 
1692   h_b->SetMarkerColor(kBlue + 2);
1693   h_b->SetLineColor(kBlue + 2);
1694   h_b->SetMarkerStyle(kFullCircle);
1695 
1696   TGraphAsymmErrors *gr_fonll_b = GetFONLL_B();
1697   gr_fonll_b->SetFillColor(kBlue - 7);
1698   gr_fonll_b->SetFillStyle(1001);
1699   TGraphAsymmErrors *gr_fonll_c = GetFONLL_C();
1700   gr_fonll_c->SetFillColor(kRed - 7);
1701   gr_fonll_c->SetFillStyle(3003);
1702   TGraph *gr_phenix = GetPHENIX_jet();
1703   gr_phenix->SetLineColor(kGray + 2);
1704   gr_phenix->SetLineWidth(3);
1705   gr_phenix->SetMarkerColor(kGray + 2);
1706   gr_phenix->SetMarkerSize(2);
1707   gr_phenix->SetMarkerStyle(kFullSquare);
1708 
1709   TCanvas *c1 = new TCanvas("Draw_HFJetTruth_DrawCrossSection_PR", "Draw_HFJetTruth_DrawCrossSection_PR", 1000, 860);
1710   //  c1->Divide(2, 2);
1711   int idx = 1;
1712   TPad *p;
1713 
1714   p = (TPad *) c1->cd(idx++);
1715   c1->Update();
1716   p->SetLogy();
1717 
1718   TH1 *hframe = p->DrawFrame(12, 0.1, 70, 1e8);
1719   hframe->SetTitle(";p_{T} [GeV/c];d^{2}#sigma/(dp_{T}d#eta) [pb/(GeV/c)]");
1720 
1721   gr_fonll_b->Draw("3");
1722   gr_phenix->Draw("pe");
1723 
1724   hall->Draw("same");
1725   h_b->Draw("same");
1726 
1727   TLegend *leg = new TLegend(0.2, 0.7, 0.95, 0.92);
1728   leg->SetFillColor(kWhite);
1729   leg->SetFillStyle(1001);
1730   //  leg->SetHeader("#splitline{#it{#bf{sPHENIX }} Simulation}{#it{p}+#it{p}, #sqrt{s} = 200 GeV, |#eta|<0.6}");
1731   leg->SetHeader("#it{#bf{sPHENIX }} Simulation, #it{p}+#it{p} #sqrt{s} = 200 GeV, |#eta|<0.6");
1732   leg->AddEntry(hall, "Inclusive jet, PYTHIA8 + CTEQ6L, anti-k_{T} R=0.4",
1733                 "lpe");
1734   leg->AddEntry(h_b, "#it{b}-quark jet, PYTHIA8 + CTEQ6L, anti-k_{T} R=0.4", "lpe");
1735   leg->AddEntry(gr_phenix, "PHENIX inclusive jet, PRL 116, 122301 (2016)", "ple");
1736   leg->AddEntry(gr_fonll_b, "#it{b}-quark, FONLL v1.3.2, CTEQ6.6", "f");
1737   leg->Draw();
1738 
1739   SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
1740 }
1741 
1742 void CrossSection2RAA_Proposal(const TString infile)
1743 {
1744   const double b_jet_RAA = 0.6;
1745 
1746   ////////////////////////////
1747   // 5-year lumi in [sPH-TRG-000]
1748   ////////////////////////////
1749 
1750   const double pp_inelastic_crosssec = 42e-3 / 1e-12;  // 42 mb in pb [sPH-TRG-000]
1751   const double AuAu_Ncoll_C0_10 = 960.2;               // [DOI:?10.1103/PhysRevC.87.034911?]
1752   const double AuAu_Ncoll_C0_20 = 770.6;               // [DOI:?10.1103/PhysRevC.91.064904?]
1753   const double AuAu_Ncoll_C0_100 = 250;                // pb^-1 [sPH-TRG-000]
1754   const double pAu_Ncoll_C0_100 = 4.7;                 // pb^-1 [sPH-TRG-000]
1755 
1756   ////////////////////////////
1757   // 2-year lumi in sPHENIX proposal
1758   ////////////////////////////
1759 
1760   //  const double pp_lumi_proposal = 630;                                                                 // Figure 4.36:
1761   //  const double AuAu_eq_lumi_C0_20_proposal = 0.6e12 * 0.2 * AuAu_Ncoll_C0_20 / pp_inelastic_crosssec;  //
1762   const double pp_lumi_proposal = 175;                                                                 // Table 4.1
1763   const double AuAu_eq_lumi_C0_20_proposal = 0.1e12 * 0.2 * AuAu_Ncoll_C0_20 / pp_inelastic_crosssec;  //
1764   const double dy_proposal = 2.;
1765   const double eff_proposal = 0.5;
1766 
1767   cout << "CrossSection2RAA_Proposal integrated luminosity assumptions in pb^-1: " << endl;
1768   cout << "\t"
1769        << "pp_lumi_proposal = " << pp_lumi_proposal << endl;
1770   cout << "\t"
1771        << "AuAu_eq_lumi_C0_20_proposal = " << AuAu_eq_lumi_C0_20_proposal << endl;
1772 
1773   TFile *f = TFile::Open(infile + "Draw_HFJetTruth_DrawCrossSection.root");
1774   assert(f);
1775 
1776   TH1F *h_b = (TH1F *) f->GetObjectChecked("h_b", "TH1F");
1777 
1778   assert(h_b);
1779 
1780   TCanvas *c1 = new TCanvas("Draw_HFJetTruth_CrossSection2RAA_Proposal", "Draw_HFJetTruth_CrossSection2RAA_Proposal", 1000, 860);
1781   c1->Divide(2, 2);
1782   int idx = 1;
1783   TPad *p;
1784 
1785   p = (TPad *) c1->cd(idx++);
1786   c1->Update();
1787   //  p->SetGridx(0);
1788   //  p->SetGridy(0);
1789   p->SetLogy();
1790 
1791   h_b->GetYaxis()->SetTitle(
1792       "d^{2}#sigma/(dp_{T}d#eta) [pb/(GeV/c)]");
1793   h_b->Draw();
1794 
1795   p = (TPad *) c1->cd(idx++);
1796   c1->Update();
1797   //  p->SetGridx(0);
1798   //  p->SetGridy(0);
1799   p->SetLogy();
1800 
1801   TH1 *h_b_int = (TH1 *) h_b->Clone(TString(h_b->GetName()) + "_IntrgratedCount");
1802   double integral = 0;
1803   for (int i = h_b_int->GetNbinsX(); i >= 1; --i)
1804   {
1805     const double cs = h_b_int->GetBinContent(i);
1806 
1807     integral += cs * h_b_int->GetXaxis()->GetBinWidth(i) * dy_proposal * 0.2 * AuAu_Ncoll_C0_20 / pp_inelastic_crosssec;
1808 
1809     h_b_int->SetBinContent(i, integral);
1810     h_b_int->SetBinError(i, 0);
1811   }
1812 
1813   h_b_int->GetYaxis()->SetTitle(
1814       "(Count > p_{T} cut)/Event 0-20% AuAu");
1815   h_b_int->Draw();
1816 
1817   p = (TPad *) c1->cd(idx++);
1818   c1->Update();
1819   //  p->SetGridx(0);
1820   //  p->SetGridy(0);
1821   //  p->SetLogy();
1822 
1823   TH1 *g_pp = CrossSection2RelUncert(h_b, 1., dy_proposal, pp_lumi_proposal * eff_proposal);
1824   TH1 *g_AA = CrossSection2RelUncert(h_b, b_jet_RAA, dy_proposal, AuAu_eq_lumi_C0_20_proposal * eff_proposal);
1825 
1826   g_pp->Draw();
1827   g_AA->Draw("same");
1828 
1829   p = (TPad *) c1->cd(idx++);
1830   c1->Update();
1831   //  p->SetGridx(0);
1832   //  p->SetGridy(0);
1833   //  p->SetLogy();
1834 
1835   p->DrawFrame(0, 0, 100, 1.2)
1836       ->SetTitle(";Transverse Momentum [GeV/#it{c}];#it{R}_{AA}");
1837   TLatex t;
1838   t.DrawLatex(10, 1, Form("#splitline{pp lumi %.0f pb^{-1}}{#splitline{AuAu 0-20 in 100B MB}{eff %.1f  }}", pp_lumi_proposal, eff_proposal));
1839 
1840   TGraphErrors *ge_RAA = GetRAA(g_pp, g_AA);
1841   ge_RAA->Draw("pe*");
1842 
1843   SaveCanvas(c1, infile + "_" + TString(c1->GetName()), kTRUE);
1844 }
1845 
1846 void Draw_HFJetTruth3yr(const TString infile =
1847                             //                         "/sphenix/user/jinhuang/HF-jet/event_gen/200pp_pythia8_CTEQ6L_20GeV/200pp_pythia8_CTEQ6L_20GeV_ALL.cfg_eneg_DSTReader.root",
1848                         //                     double int_lumi = 210715 / 5.533e-05 / 1e9, const double dy = 0.6 * 2)
1849                         "/sphenix/user/jinhuang/HF-jet/event_gen/200pp_pythia8_CTEQ6L_7GeV/200pp_pythia8_CTEQ6L_7GeV_ALL.cfg_eneg_DSTReader.root",
1850                         double int_lumi = 891093 / 3.332e-02 / 1e9, const double dy = 0.6 * 2)
1851 //"/sphenix/user/jinhuang/HF-jet/event_gen/200pp_pythia8_CTEQ6L/200pp_pythia8_CTEQ6L_111ALL.cfg_eneg_DSTReader.root",
1852 //double int_lumi = 789908/4.631e-03 / 1e9, const double dy = 0.6*2)
1853 //Draw_HFJetTruth(const TString infile =
1854 //    "../macros3/G4sPHENIXCells.root_DSTReader.root",
1855 //    double int_lumi = 789908/4.631e-03 / 1e9, const double dy = 0.6*2)
1856 //Draw_HFJetTruth(const TString infile =
1857 //    "../macros2/G4sPHENIXCells.root_DSTReader.root",
1858 //    const double int_lumi = 1842152 / 5.801e-03 / 1e9,
1859 //    const double dy = 0.6 * 2)
1860 {
1861   SetsPhenixStyle();
1862   gStyle->SetOptStat(0);
1863   gStyle->SetOptFit(1111);
1864   TVirtualFitter::SetDefaultFitter("Minuit2");
1865 
1866   gSystem->Load("libg4eval.so");
1867 
1868   if (!_file0)
1869   {
1870     TString chian_str = infile;
1871     chian_str.ReplaceAll("ALL", "*");
1872 
1873     TChain *t = new TChain("T");
1874     const int n = t->Add(chian_str);
1875 
1876     int_lumi *= n;
1877     //      cout << "Loaded " << n << " root files with " << chian_str << endl;
1878     cout << "int_lumi = " << int_lumi << " pb^-1 from " << n << " root files with " << chian_str << endl;
1879     assert(n > 0);
1880 
1881     T = t;
1882 
1883     _file0 = new TFile;
1884     _file0->SetName(infile);
1885   }
1886 
1887   //      DrawCrossSection(int_lumi, dy);
1888   //  Draw_HFJetTruth_DrawCrossSection_PR(infile);
1889   //  CrossSection2RAA_Proposal(infile);
1890   //    CrossSection2RAA(infile);
1891   const double acceptance1 = 2. * (1.1 - .4);
1892 
1893   CrossSection2zg(infile, false, acceptance1);
1894 
1895   //  CrossSection2RAA(infile, false, acceptance1);
1896 
1897   //  const double acceptance2 = 2.* (0.85 - .4);
1898   //  CrossSection2RAA(infile, false, acceptance2);
1899 
1900   //  CrossSection2v2(infile, false, .5, acceptance1);
1901   //  CrossSection2v2(infile, false, .4, acceptance);
1902 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!