Back to home page

sPhenix code displayed by LXR
 
 

    

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

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 <TFile.h>
0012 #include <TGraph.h>
0013 #include <TGraphAsymmErrors.h>
0014 #include <TGraphErrors.h>
0015 #include <TLatex.h>
0016 #include <TLine.h>
0017 #include <TString.h>
0018 #include <TTree.h>
0019 #include <cassert>
0020 #include <cmath>
0021 #include "SaveCanvas.C"
0022 #include "sPhenixStyle.C"
0023 
0024 TFile *_file0 = NULL;
0025 TTree *T = NULL;
0026 
0027 void Draw_HFJetTruth(const TString infile =
0028                          //                         "/sphenix/user/jinhuang/HF-jet/event_gen/200pp_pythia8_CTEQ6L_20GeV/200pp_pythia8_CTEQ6L_20GeV_ALL.cfg_eneg_DSTReader.root",
0029                      //                     double int_lumi = 210715 / 5.533e-05 / 1e9, const double dy = 0.6 * 2)
0030                      "/sphenix/user/jinhuang/HF-jet/event_gen/200pp_pythia8_CTEQ6L_7GeV/200pp_pythia8_CTEQ6L_7GeV_ALL.cfg_eneg_DSTReader.root",
0031                      double int_lumi = 891093 / 3.332e-02 / 1e9, const double dy = 0.6 * 2)
0032 //"/sphenix/user/jinhuang/HF-jet/event_gen/200pp_pythia8_CTEQ6L/200pp_pythia8_CTEQ6L_111ALL.cfg_eneg_DSTReader.root",
0033 //double int_lumi = 789908/4.631e-03 / 1e9, const double dy = 0.6*2)
0034 //Draw_HFJetTruth(const TString infile =
0035 //    "../macros3/G4sPHENIXCells.root_DSTReader.root",
0036 //    double int_lumi = 789908/4.631e-03 / 1e9, const double dy = 0.6*2)
0037 //Draw_HFJetTruth(const TString infile =
0038 //    "../macros2/G4sPHENIXCells.root_DSTReader.root",
0039 //    const double int_lumi = 1842152 / 5.801e-03 / 1e9,
0040 //    const double dy = 0.6 * 2)
0041 {
0042   SetsPhenixStyle();
0043   gStyle->SetOptStat(0);
0044   gStyle->SetOptFit(1111);
0045   TVirtualFitter::SetDefaultFitter("Minuit2");
0046 
0047   gSystem->Load("libg4eval.so");
0048 
0049   if (!_file0)
0050   {
0051     TString chian_str = infile;
0052     chian_str.ReplaceAll("ALL", "*");
0053 
0054     TChain *t = new TChain("T");
0055     const int n = t->Add(chian_str);
0056 
0057     int_lumi *= n;
0058     //      cout << "Loaded " << n << " root files with " << chian_str << endl;
0059     cout << "int_lumi = " << int_lumi << " pb^-1 from " << n << " root files with " << chian_str << endl;
0060     assert(n > 0);
0061 
0062     T = t;
0063 
0064     _file0 = new TFile;
0065     _file0->SetName(infile);
0066   }
0067 
0068 //      DrawCrossSection(int_lumi, dy);
0069   //  Draw_HFJetTruth_DrawCrossSection_PR(infile);
0070   //  CrossSection2RAA_Proposal(infile);
0071   //    CrossSection2RAA(infile);
0072   const double acceptance1 = 2.* (1.1 - .4);
0073   CrossSection2RAA(infile, false, acceptance1);
0074 
0075 //  const double acceptance2 = 2.* (0.85 - .4);
0076 //  CrossSection2RAA(infile, false, acceptance2);
0077 
0078 //  CrossSection2v2(infile, false, .7, acceptance);
0079 //  CrossSection2v2(infile, false, .4, acceptance);
0080 }
0081 
0082 void DrawCrossSection(double int_lumi, const double dy)
0083 {
0084   TH1 *hall = new TH1F("hall", ";p_{T} [GeV/c]", 100, 0, 100);
0085   TH1 *h_b = new TH1F("h_b", ";p_{T} [GeV/c]", 100, 0, 100);
0086   TH1 *h_bq = new TH1F("h_bq", ";p_{T} [GeV/c]", 100, 0, 100);
0087   TH1 *h_bh = new TH1F("h_bh", ";p_{T} [GeV/c]", 100, 0, 100);
0088   TH1 *h_bh5 = new TH1F("h_bh5", ";p_{T} [GeV/c]", 100, 0, 100);
0089   TH1 *h_ch5 = new TH1F("h_ch5", ";p_{T} [GeV/c]", 100, 0, 100);
0090   TH1 *h_c = new TH1F("h_c", ";p_{T} [GeV/c]", 100, 0, 100);
0091 
0092   T->Draw("AntiKt_Truth_r04.get_pt()>>hall",
0093           "AntiKt_Truth_r04.get_pt()>15 && abs(AntiKt_Truth_r04.get_eta())<0.6",
0094           "goff");
0095 
0096   T->Draw("AntiKt_Truth_r04.get_pt()>>h_b",
0097           "AntiKt_Truth_r04.get_pt()>15 && abs(AntiKt_Truth_r04.get_eta())<0.6 && abs(AntiKt_Truth_r04.get_property(1000))==5",
0098           "goff");
0099   T->Draw(
0100       "AntiKt_Truth_r04.get_pt()* AntiKt_Truth_r04.get_property(1001) >>h_bq",
0101       "AntiKt_Truth_r04.get_pt()>15 && abs(AntiKt_Truth_r04.get_eta())<0.6 && abs(AntiKt_Truth_r04.get_property(1000))==5",
0102       "goff");
0103 
0104   //    T->Draw("AntiKt_Truth_r04.get_pt()>>h_bh",
0105   //        "AntiKt_Truth_r04.get_pt()>15 && abs(AntiKt_Truth_r04.get_eta())<0.6 && AntiKt_Truth_r04.get_property(1010)==5",
0106   //        "goff");
0107 
0108   T->Draw("AntiKt_Truth_r04.get_pt()>>h_bh5",
0109           "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",
0110           "goff");
0111   //
0112   T->Draw("AntiKt_Truth_r04.get_pt()>>h_c",
0113           "AntiKt_Truth_r04.get_pt()>15 && abs(AntiKt_Truth_r04.get_eta())<0.6 && abs(AntiKt_Truth_r04.get_property(1000))==4",
0114           "goff");
0115   T->Draw("AntiKt_Truth_r04.get_pt()>>h_ch5",
0116           "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",
0117           "goff");
0118 
0119   Convert2CrossSection(hall, int_lumi, dy);
0120   Convert2CrossSection(h_b, int_lumi, dy);
0121   Convert2CrossSection(h_bq, int_lumi, dy);
0122   Convert2CrossSection(h_bh5, int_lumi, dy);
0123   Convert2CrossSection(h_ch5, int_lumi, dy);
0124   Convert2CrossSection(h_c, int_lumi, dy);
0125 
0126   h_b->SetLineColor(kBlue);
0127   h_b->SetMarkerColor(kBlue);
0128 
0129   h_bq->SetLineColor(kBlue + 3);
0130   h_bq->SetMarkerColor(kBlue + 3);
0131 
0132   h_bh5->SetLineColor(kBlue + 3);
0133   h_bh5->SetMarkerColor(kBlue + 3);
0134   h_bh5->SetMarkerStyle(kStar);
0135 
0136   h_c->SetLineColor(kRed);
0137   h_c->SetMarkerColor(kRed);
0138 
0139   h_ch5->SetLineColor(kRed + 3);
0140   h_ch5->SetMarkerColor(kRed + 3);
0141   h_ch5->SetMarkerStyle(kStar);
0142 
0143   TGraphAsymmErrors *gr_fonll_b = GetFONLL_B();
0144   gr_fonll_b->SetFillColor(kBlue - 7);
0145   gr_fonll_b->SetFillStyle(3002);
0146   TGraphAsymmErrors *gr_fonll_c = GetFONLL_C();
0147   gr_fonll_c->SetFillColor(kRed - 7);
0148   gr_fonll_c->SetFillStyle(3003);
0149   TGraph *gr_phenix = GetPHENIX_jet();
0150   gr_phenix->SetLineColor(kGray + 2);
0151   gr_phenix->SetMarkerColor(kBlack);
0152   gr_phenix->SetMarkerStyle(kOpenCross);
0153 
0154   TCanvas *c1 = new TCanvas("Draw_HFJetTruth_DrawCrossSection", "Draw_HFJetTruth_DrawCrossSection", 1000, 860);
0155   //  c1->Divide(2, 2);
0156   int idx = 1;
0157   TPad *p;
0158 
0159   p = (TPad *) c1->cd(idx++);
0160   c1->Update();
0161   p->SetGridx(0);
0162   p->SetGridy(0);
0163   p->SetLogy();
0164 
0165   hall->Draw();
0166 
0167   gr_fonll_b->Draw("3");
0168   gr_fonll_c->Draw("3");
0169   gr_phenix->Draw("pe");
0170 
0171   h_b->Draw("same");
0172   h_bh5->Draw("same");
0173   h_bq->Draw("same");
0174   h_c->Draw("same");
0175   h_ch5->Draw("same");
0176 
0177   hall->GetXaxis()->SetRangeUser(12, 60);
0178   hall->GetYaxis()->SetTitle(
0179       "d^{2}#sigma/(dp_{T}dy), d^{2}#sigma/(dp_{T}d#eta) [pb/(GeV/c)]");
0180 
0181   TLegend *leg = new TLegend(0.45, 0.6, 0.99, 0.95);
0182   leg->SetFillColor(kWhite);
0183   leg->SetFillStyle(1001);
0184   leg->SetHeader("p+p collisions @ sPHENIX, #sqrt{s} = 200 GeV, |#eta|<0.6");
0185   leg->AddEntry(hall, "Inclusive jet, Pythia8, Truth, anti-k_{t}, R=0.4",
0186                 "lpe");
0187   leg->AddEntry(h_c, "c-quark jet, Pythia8, Truth, anti-k_{t}, R=0.4", "lpe");
0188   leg->AddEntry(h_b, "b-quark jet, Pythia8, Truth, anti-k_{t}, R=0.4", "lpe");
0189   leg->AddEntry(h_bq, "Leading b-quark in jet, Pythia8, anti-k_{t}, R=0.4", "lpe");
0190   leg->AddEntry(h_bh5,
0191                 "b-hadron jet, Pythia8, Truth, anti-k_{t}, R=0.4, p_{T, b-hadron}>5 GeV/c",
0192                 "lpe");
0193   leg->AddEntry(gr_phenix, "PHENIX inclusive jet, PRL 116, 122301 (2016)", "ple");
0194   leg->AddEntry(gr_fonll_c, "c-quark, FONLL v1.3.2, CTEQ6.6, |y|<0.6", "f");
0195   leg->AddEntry(gr_fonll_b, "b-quark, FONLL v1.3.2, CTEQ6.6, |y|<0.6", "f");
0196   leg->Draw();
0197 
0198   SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
0199 }
0200 
0201 void Draw_HFJetTruth_DrawCrossSection_PR(const TString infile)
0202 {
0203   TFile *f = TFile::Open(infile + "Draw_HFJetTruth_DrawCrossSection.root");
0204   assert(f);
0205 
0206   TH1F *hall = (TH1F *) f->GetObjectChecked("hall", "TH1F");
0207   assert(hall);
0208   TH1F *h_b = (TH1F *) f->GetObjectChecked("h_b", "TH1F");
0209   assert(h_b);
0210 
0211   hall->SetMarkerColor(kBlack);
0212   hall->SetLineColor(kBlack);
0213   hall->SetMarkerStyle(kFullCircle);
0214 
0215   h_b->SetMarkerColor(kBlue + 2);
0216   h_b->SetLineColor(kBlue + 2);
0217   h_b->SetMarkerStyle(kFullCircle);
0218 
0219   TGraphAsymmErrors *gr_fonll_b = GetFONLL_B();
0220   gr_fonll_b->SetFillColor(kBlue - 7);
0221   gr_fonll_b->SetFillStyle(1001);
0222   TGraphAsymmErrors *gr_fonll_c = GetFONLL_C();
0223   gr_fonll_c->SetFillColor(kRed - 7);
0224   gr_fonll_c->SetFillStyle(3003);
0225   TGraph *gr_phenix = GetPHENIX_jet();
0226   gr_phenix->SetLineColor(kGray + 2);
0227   gr_phenix->SetLineWidth(3);
0228   gr_phenix->SetMarkerColor(kGray + 2);
0229   gr_phenix->SetMarkerSize(2);
0230   gr_phenix->SetMarkerStyle(kFullSquare);
0231 
0232   TCanvas *c1 = new TCanvas("Draw_HFJetTruth_DrawCrossSection_PR", "Draw_HFJetTruth_DrawCrossSection_PR", 1000, 860);
0233   //  c1->Divide(2, 2);
0234   int idx = 1;
0235   TPad *p;
0236 
0237   p = (TPad *) c1->cd(idx++);
0238   c1->Update();
0239   p->SetLogy();
0240 
0241   TH1 *hframe = p->DrawFrame(12, 0.1, 70, 1e8);
0242   hframe->SetTitle(";p_{T} [GeV/c];d^{2}#sigma/(dp_{T}d#eta) [pb/(GeV/c)]");
0243 
0244   gr_fonll_b->Draw("3");
0245   gr_phenix->Draw("pe");
0246 
0247   hall->Draw("same");
0248   h_b->Draw("same");
0249 
0250   TLegend *leg = new TLegend(0.2, 0.7, 0.95, 0.92);
0251   leg->SetFillColor(kWhite);
0252   leg->SetFillStyle(1001);
0253   //  leg->SetHeader("#splitline{#it{#bf{sPHENIX }} Simulation}{p+p, #sqrt{s} = 200 GeV, |#eta|<0.6}");
0254   leg->SetHeader("#it{#bf{sPHENIX }} Simulation, #it{p}+#it{p} #sqrt{s} = 200 GeV, |#eta|<0.6");
0255   leg->AddEntry(hall, "Inclusive jet, PYTHIA8 + CTEQ6L, anti-k_{T} R=0.4",
0256                 "lpe");
0257   leg->AddEntry(h_b, "#it{b}-quark jet, PYTHIA8 + CTEQ6L, anti-k_{T} R=0.4", "lpe");
0258   leg->AddEntry(gr_phenix, "PHENIX inclusive jet, PRL 116, 122301 (2016)", "ple");
0259   leg->AddEntry(gr_fonll_b, "#it{b}-quark, FONLL v1.3.2, CTEQ6.6", "f");
0260   leg->Draw();
0261 
0262   SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
0263 }
0264 
0265 void CrossSection2RAA_Proposal(const TString infile)
0266 {
0267   const double b_jet_RAA = 0.6;
0268 
0269   ////////////////////////////
0270   // 5-year lumi in [sPH-TRG-000]
0271   ////////////////////////////
0272 
0273   const double pp_inelastic_crosssec = 42e-3 / 1e-12;  // 42 mb in pb [sPH-TRG-000]
0274   const double AuAu_Ncoll_C0_10 = 960.2;               // [DOI:?10.1103/PhysRevC.87.034911?]
0275   const double AuAu_Ncoll_C0_20 = 770.6;               // [DOI:?10.1103/PhysRevC.91.064904?]
0276   const double AuAu_Ncoll_C0_100 = 250;                // pb^-1 [sPH-TRG-000]
0277   const double pAu_Ncoll_C0_100 = 4.7;                 // pb^-1 [sPH-TRG-000]
0278 
0279   ////////////////////////////
0280   // 2-year lumi in sPHENIX proposal
0281   ////////////////////////////
0282 
0283   //  const double pp_lumi_proposal = 630;                                                                 // Figure 4.36:
0284   //  const double AuAu_eq_lumi_C0_20_proposal = 0.6e12 * 0.2 * AuAu_Ncoll_C0_20 / pp_inelastic_crosssec;  //
0285   const double pp_lumi_proposal = 175;                                                                 // Table 4.1
0286   const double AuAu_eq_lumi_C0_20_proposal = 0.1e12 * 0.2 * AuAu_Ncoll_C0_20 / pp_inelastic_crosssec;  //
0287   const double dy_proposal = 2.;
0288   const double eff_proposal = 0.5;
0289   const double dy_proposal = 2;
0290 
0291   cout << "CrossSection2RAA_Proposal integrated luminosity assumptions in pb^-1: " << endl;
0292   cout << "\t"
0293        << "pp_lumi_proposal = " << pp_lumi_proposal << endl;
0294   cout << "\t"
0295        << "AuAu_eq_lumi_C0_20_proposal = " << AuAu_eq_lumi_C0_20_proposal << endl;
0296 
0297   TFile *f = TFile::Open(infile + "Draw_HFJetTruth_DrawCrossSection.root");
0298   assert(f);
0299 
0300   TH1F *h_b = (TH1F *) f->GetObjectChecked("h_b", "TH1F");
0301 
0302   assert(h_b);
0303 
0304   TCanvas *c1 = new TCanvas("Draw_HFJetTruth_CrossSection2RAA_Proposal", "Draw_HFJetTruth_CrossSection2RAA_Proposal", 1000, 860);
0305   c1->Divide(2, 2);
0306   int idx = 1;
0307   TPad *p;
0308 
0309   p = (TPad *) c1->cd(idx++);
0310   c1->Update();
0311   //  p->SetGridx(0);
0312   //  p->SetGridy(0);
0313   p->SetLogy();
0314 
0315   h_b->GetYaxis()->SetTitle(
0316       "d^{2}#sigma/(dp_{T}d#eta) [pb/(GeV/c)]");
0317   h_b->Draw();
0318 
0319   p = (TPad *) c1->cd(idx++);
0320   c1->Update();
0321   //  p->SetGridx(0);
0322   //  p->SetGridy(0);
0323   p->SetLogy();
0324 
0325   TH1 *h_b_int = (TH1 *) h_b->Clone(TString(h_b->GetName()) + "_IntrgratedCount");
0326   double integral = 0;
0327   for (int i = h_b_int->GetNbinsX(); i >= 1; --i)
0328   {
0329     const double cs = h_b_int->GetBinContent(i);
0330 
0331     integral += cs * h_b_int->GetXaxis()->GetBinWidth(i) * dy_proposal * 0.2 * AuAu_Ncoll_C0_20 / pp_inelastic_crosssec;
0332 
0333     const double cs = h_b_int->SetBinContent(i, integral);
0334     const double cs = h_b_int->SetBinError(i, 0);
0335   }
0336 
0337   h_b_int->GetYaxis()->SetTitle(
0338       "(Count > p_{T} cut)/Event 0-20% AuAu");
0339   h_b_int->Draw();
0340 
0341   p = (TPad *) c1->cd(idx++);
0342   c1->Update();
0343   //  p->SetGridx(0);
0344   //  p->SetGridy(0);
0345   //  p->SetLogy();
0346 
0347   TH1 *g_pp = CrossSection2RelUncert(h_b, 1., dy_proposal, pp_lumi_proposal * eff_proposal);
0348   TH1 *g_AA = CrossSection2RelUncert(h_b, b_jet_RAA, dy_proposal, AuAu_eq_lumi_C0_20_proposal * eff_proposal);
0349 
0350   g_pp->Draw();
0351   g_AA->Draw("same");
0352 
0353   p = (TPad *) c1->cd(idx++);
0354   c1->Update();
0355   //  p->SetGridx(0);
0356   //  p->SetGridy(0);
0357   //  p->SetLogy();
0358 
0359   p->DrawFrame(0, 0, 100, 1.2)
0360       ->SetTitle(";Transverse Momentum [GeV/#it{c}];#it{R}_{AA}");
0361   TLatex t;
0362   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));
0363 
0364   TGraphErrors *ge_RAA = GetRAA(g_pp, g_AA);
0365   ge_RAA->Draw("pe*");
0366 
0367   SaveCanvas(c1, infile + "_" + TString(c1->GetName()), kTRUE);
0368 }
0369 
0370 void CrossSection2RAA(const TString infile, const bool use_AA_jet_trigger = true, const double dy = .7 * 2)
0371 {
0372   TFile *f = TFile::Open(infile + "Draw_HFJetTruth_DrawCrossSection.root");
0373   assert(f);
0374 
0375   TH1F *hall = (TH1F *) f->GetObjectChecked("hall", "TH1F");
0376   assert(hall);
0377   TH1F *h_b = (TH1F *) f->GetObjectChecked("h_b", "TH1F");
0378   assert(h_b);
0379 
0380   const TString s_suffix(use_AA_jet_trigger ? "_AAJetTriggered" : "");
0381   s_suffix += Form("_deta%.2f",dy/2);
0382 
0383   const double b_jet_RAA = 0.6;
0384 
0385   const double pp_eff = 0.6;
0386   const double pp_purity = 0.4;
0387   const double AuAu_eff = 0.4;
0388   const double AuAu_purity = 0.4;
0389 
0390   ////////////////////////////
0391   // 5-year lumi in [sPH-TRG-000]
0392   ////////////////////////////
0393 
0394   const double pp_lumi = 200;                          // pb^-1 [sPH-TRG-000], rounded up from 197 pb^-1
0395   const double pp_inelastic_crosssec = 42e-3 / 1e-12;  // 42 mb in pb [sPH-TRG-000]
0396 
0397   const double AuAu_MB_Evt = use_AA_jet_trigger ? 550e9 : 240e9;  // [sPH-TRG-000], depending on whether jet trigger applied in AA collisions
0398   const double pAu_MB_Evt = 600e9;                                // [sPH-TRG-000]
0399 
0400   const double AuAu_Ncoll_C0_10 = 960.2;  // [DOI:?10.1103/PhysRevC.87.034911?]
0401   const double AuAu_Ncoll_C0_20 = 770.6;  // [DOI:?10.1103/PhysRevC.91.064904?]
0402   const double AuAu_Ncoll_C0_100 = 250;   // pb^-1 [sPH-TRG-000]
0403   const double pAu_Ncoll_C0_100 = 4.7;    // pb^-1 [sPH-TRG-000]
0404 
0405   const double AuAu_eq_lumi_C0_10 = AuAu_MB_Evt * 0.1 * AuAu_Ncoll_C0_10 / pp_inelastic_crosssec;  //
0406   const double AuAu_eq_lumi_C0_20 = AuAu_MB_Evt * 0.2 * AuAu_Ncoll_C0_20 / pp_inelastic_crosssec;  //
0407   const double AuAu_eq_lumi_C0_100 = AuAu_MB_Evt * 1 * AuAu_Ncoll_C0_100 / pp_inelastic_crosssec;  //
0408 
0409   const double pAu_eq_lumi_C0_100 = pAu_MB_Evt * 1 * pAu_Ncoll_C0_100 / pp_inelastic_crosssec;  //
0410 
0411   cout << "CrossSection2RAA integrated luminosity assumptions in pb^-1: " << endl;
0412   cout << "\t"
0413        << "pp_lumi = " << pp_lumi << endl;
0414   cout << "\t"
0415        << "AuAu_eq_lumi_C0_10 = " << AuAu_eq_lumi_C0_10 << endl;
0416   cout << "\t"
0417        << "AuAu_eq_lumi_C0_20 = " << AuAu_eq_lumi_C0_20 << endl;
0418   cout << "\t"
0419        << "AuAu_eq_lumi_C0_100 = " << AuAu_eq_lumi_C0_100 << endl;
0420   cout << "\t"
0421        << "pAu_eq_lumi_C0_100 = " << pAu_eq_lumi_C0_100 << endl;
0422 
0423   TCanvas *c1 = new TCanvas("Draw_HFJetTruth_CrossSection2RAA_Ratio" + s_suffix, "Draw_HFJetTruth_CrossSection2RAA_Ratio" + s_suffix, 700, 600);
0424   c1->Divide(1, 1);
0425   int idx = 1;
0426   TPad *p;
0427 
0428   p = (TPad *) c1->cd(idx++);
0429   c1->Update();
0430   //  p->SetGridx(0);
0431   //  p->SetGridy(0);
0432   //  p->SetLogy();
0433 
0434   TH1 *g_pp = CrossSection2RelUncert(h_b, 1., dy, pp_lumi * pp_eff * pp_purity);
0435   TH1 *g_AA = CrossSection2RelUncert(h_b, b_jet_RAA, dy, AuAu_eq_lumi_C0_10 * AuAu_eff * AuAu_purity);
0436 
0437   g_pp->SetLineColor(kRed);
0438   g_AA->SetLineColor(kBlue);
0439 
0440   g_pp->Draw();
0441   g_AA->Draw("same");
0442   SaveCanvas(c1, infile + "_" + TString(c1->GetName()), kTRUE);
0443 
0444   TCanvas *c1 = new TCanvas("Draw_HFJetTruth_CrossSection2RAA" + s_suffix, "Draw_HFJetTruth_CrossSection2RAA" + s_suffix, 700, 600);
0445   c1->Divide(1, 1);
0446   int idx = 1;
0447   TPad *p;
0448 
0449   p = (TPad *) c1->cd(idx++);
0450   c1->Update();
0451 
0452   p->DrawFrame(11, 0, 50, 1.2)
0453       ->SetTitle(";Transverse Momentum [GeV/#it{c}];#it{R}_{AA}");
0454 
0455   TGraphErrors *ge_RAA = GetRAA(g_pp, g_AA);
0456 
0457   ge_RAA->SetLineWidth(3);
0458   ge_RAA->SetMarkerStyle(kFullCircle);
0459   ge_RAA->SetMarkerSize(2);
0460 
0461   ge_RAA->Draw("pe");
0462   ge_RAA->Print();
0463 
0464   TLegend *leg = new TLegend(.0, .76, .85, .93);
0465   leg->SetFillStyle(0);
0466   leg->AddEntry("", "#it{#bf{sPHENIX }} Simulation", "");
0467   leg->AddEntry("", Form("PYTHIA-8 #it{b}-jet, Anti-k_{T} R=0.4, |#eta|<%.2f, CTEQ6L", dy / 2), "");
0468   leg->AddEntry("", Form("#it{p}+#it{p}: %.0f pb^{-1}, %.0f%% Eff., %.0f%% Pur.", pp_lumi, pp_eff * 100, pp_purity * 100), "");
0469   leg->AddEntry("", Form("Au+Au: %.0fB col., %.0f%% Eff., %.0f%% Pur.", '%', AuAu_MB_Evt / 1e9, AuAu_eff * 100, AuAu_purity * 100), "");
0470   leg->Draw();
0471 
0472   TLegend *leg2 = new TLegend(.17, .70, .88, .77);
0473   leg2->AddEntry(ge_RAA, "#it{b}-jet #it{R}_{AA}, Au+Au 0-10%C, #sqrt{s_{NN}}=200 GeV", "pl");
0474   leg2->Draw();
0475 
0476   SaveCanvas(c1, infile + "_" + TString(c1->GetName()), kTRUE);
0477 
0478   TCanvas *c1 = new TCanvas("Draw_HFJetTruth_CrossSection2RAA_Theory" + s_suffix, "Draw_HFJetTruth_CrossSection2RAA_Theory" + s_suffix, 700, 600);
0479   c1->Divide(1, 1);
0480   int idx = 1;
0481   TPad *p;
0482 
0483   p = (TPad *) c1->cd(idx++);
0484   c1->Update();
0485 
0486   p->DrawFrame(11, 0, 50, 1.2)
0487       ->SetTitle(";Transverse Momentum [GeV/#it{c}];#it{R}_{AA}");
0488 
0489   TGraph *g20 = pQCDModel_HuangKangVitev(2.0);
0490   TGraph *g22 = pQCDModel_HuangKangVitev(2.2);
0491 
0492   g20->SetLineColor(kGreen + 3);
0493   g22->SetLineColor(kRed + 3);
0494   g20->SetLineWidth(3);
0495   g22->SetLineWidth(3);
0496 
0497   g20->Draw("l");
0498   g22->Draw("l");
0499 
0500   ge_RAA->DrawClone("pe");
0501   leg->DrawClone();
0502   leg2->DrawClone();
0503 
0504   TLegend *leg1 = new TLegend(.2, .20, .85, .37);
0505   leg1->SetHeader("#splitline{pQCD, Phys.Lett. B726 (2013) 251-256}{#sqrt{s_{NN}}=200 GeV, #it{b}-jet R=0.4}");
0506   leg1->AddEntry("", "", "");
0507   leg1->AddEntry(g20, "g^{med} = 2.0", "l");
0508   leg1->AddEntry(g22, "g^{med} = 2.2", "l");
0509   leg1->Draw();
0510 
0511   SaveCanvas(c1, infile + "_" + TString(c1->GetName()), kTRUE);
0512 }
0513 
0514 void CrossSection2v2(const TString infile, const bool use_AA_jet_trigger = true,const double ep_resolution = 0.7, const double dy = .7 * 2)
0515 {
0516   TFile *f = TFile::Open(infile + "Draw_HFJetTruth_DrawCrossSection.root");
0517   assert(f);
0518 
0519   TH1F *hall = (TH1F *) f->GetObjectChecked("hall", "TH1F");
0520   assert(hall);
0521   TH1F *h_b = (TH1F *) f->GetObjectChecked("h_b", "TH1F");
0522   assert(h_b);
0523 
0524   const double b_jet_RAA = 0.6;
0525 
0526   const double pp_eff = 0.6;
0527   const double pp_purity = 0.4;
0528   const double AuAu_eff = 0.4;
0529   const double AuAu_purity = 0.4;
0530 
0531   ////////////////////////////
0532   // 5-year lumi in [sPH-TRG-000]
0533   ////////////////////////////
0534 
0535   const double pp_lumi = 200;                          // pb^-1 [sPH-TRG-000], rounded up from 197 pb^-1
0536   const double pp_inelastic_crosssec = 42e-3 / 1e-12;  // 42 mb in pb [sPH-TRG-000]
0537 
0538   const double AuAu_MB_Evt = use_AA_jet_trigger ? 550e9 : 240e9;  // [sPH-TRG-000], depending on whether jet trigger applied in AA collisions
0539   const double pAu_MB_Evt = 600e9;                                // [sPH-TRG-000]
0540 
0541   const double AuAu_Ncoll_C0_10 = 960.2;  // [DOI:?10.1103/PhysRevC.87.034911?]
0542   const double AuAu_Ncoll_C0_20 = 770.6;  // [DOI:?10.1103/PhysRevC.91.064904?]
0543   const double AuAu_Ncoll_C10_20 = 603;   // [sPH-HF-2017-001-v1]
0544   const double AuAu_Ncoll_C20_40 = 296;   // [sPH-HF-2017-001-v1]
0545   const double AuAu_Ncoll_C10_40 =  (AuAu_Ncoll_C10_20*10 + AuAu_Ncoll_C20_40*20)/30;
0546   const double AuAu_Ncoll_C40_60 = 94;    //  [sPH-HF-2017-001-v1]
0547   const double AuAu_Ncoll_C60_92 = 15;    //  [sPH-HF-2017-001-v1]
0548   const double AuAu_Ncoll_C0_100 = 250;   // pb^-1 [sPH-TRG-000]
0549   const double pAu_Ncoll_C0_100 = 4.7;    // pb^-1 [sPH-TRG-000]
0550 
0551   const double AuAu_eq_lumi_C0_10 = AuAu_MB_Evt * 0.1 * AuAu_Ncoll_C0_10 / pp_inelastic_crosssec;  //
0552   const double AuAu_eq_lumi_C0_20 = AuAu_MB_Evt * 0.2 * AuAu_Ncoll_C0_20 / pp_inelastic_crosssec;  //
0553   const double AuAu_eq_lumi_C0_100 = AuAu_MB_Evt * 1 * AuAu_Ncoll_C0_100 / pp_inelastic_crosssec;  //
0554 
0555   const double AuAu_eq_lumi_C10_20 = AuAu_MB_Evt * .1 * AuAu_Ncoll_C10_20 / pp_inelastic_crosssec;          //
0556   const double AuAu_eq_lumi_C20_40 = AuAu_MB_Evt * .2 * AuAu_Ncoll_C20_40 / pp_inelastic_crosssec;          //
0557   const double AuAu_eq_lumi_C10_40 = AuAu_MB_Evt * .3 * AuAu_Ncoll_C10_40 / pp_inelastic_crosssec;          //
0558   const double AuAu_eq_lumi_C40_60 = AuAu_MB_Evt * .2 * AuAu_Ncoll_C40_60 / pp_inelastic_crosssec;          //
0559   const double AuAu_eq_lumi_C60_92 = AuAu_MB_Evt * (.92 - .6) * AuAu_Ncoll_C60_92 / pp_inelastic_crosssec;  //
0560 
0561   const double pAu_eq_lumi_C0_100 = pAu_MB_Evt * 1 * pAu_Ncoll_C0_100 / pp_inelastic_crosssec;  //
0562 
0563   ;
0564 
0565 
0566   const TString s_suffix(use_AA_jet_trigger ? "_AAJetTriggered" : "");
0567   s_suffix += Form("_EPR%.1f",ep_resolution);
0568   s_suffix += Form("_deta%.2f",dy/2);
0569 
0570   cout << "CrossSection2v2 integrated luminosity assumptions in pb^-1: " << endl;
0571   cout << "\t"
0572        << "pp_lumi = " << pp_lumi << endl;
0573   cout << "\t"
0574        << "AuAu_eq_lumi_C0_10 = " << AuAu_eq_lumi_C0_10 << endl;
0575   cout << "\t"
0576        << "AuAu_eq_lumi_C0_20 = " << AuAu_eq_lumi_C0_20 << endl;
0577   cout << "\t"
0578        << "AuAu_eq_lumi_C0_100 = " << AuAu_eq_lumi_C0_100 << endl;
0579   cout << "\t"
0580        << "AuAu_eq_lumi_C10_20 = " << AuAu_eq_lumi_C10_20 << endl;
0581   cout << "\t"
0582        << "AuAu_eq_lumi_C20_40 = " << AuAu_eq_lumi_C20_40 << endl;
0583   cout << "\t"
0584        << "AuAu_eq_lumi_C10_40 = " << AuAu_eq_lumi_C10_40 << endl;
0585   cout << "\t"
0586        << "AuAu_eq_lumi_C40_60 = " << AuAu_eq_lumi_C40_60 << endl;
0587   cout << "\t"
0588        << "AuAu_eq_lumi_C60_92 = " << AuAu_eq_lumi_C60_92 << endl;
0589   cout << "\t"
0590        << "pAu_eq_lumi_C0_100 = " << pAu_eq_lumi_C0_100 << endl;
0591 
0592   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);
0593   TGraph *g_AA_C0_20 = CrossSection2v2Uncert(h_b, b_jet_RAA, dy, AuAu_eq_lumi_C0_20 * AuAu_eff * AuAu_purity, ep_resolution, 0);
0594   TGraph *g_AA_C10_20 = CrossSection2v2Uncert(h_b, b_jet_RAA, dy, AuAu_eq_lumi_C10_20 * AuAu_eff * AuAu_purity, ep_resolution, 0);
0595   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);
0596   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);
0597   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);
0598   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);
0599   //
0600   g_AA_C0_10->SetLineColor(kBlue+3);
0601   g_AA_C10_20->SetLineColor(kAzure+3);
0602   g_AA_C20_40->SetLineColor(kTeal+3);
0603   g_AA_C10_40->SetLineColor(kTeal+3);
0604   g_AA_C40_60->SetLineColor(kSpring+3);
0605 
0606   g_AA_C0_10->SetMarkerColor(kBlue+3);
0607   g_AA_C10_20->SetMarkerColor(kAzure+3);
0608   g_AA_C20_40->SetMarkerColor(kTeal+3);
0609   g_AA_C10_40->SetMarkerColor(kTeal+3);
0610   g_AA_C40_60->SetMarkerColor(kSpring+3);
0611 
0612   g_AA_C0_10->SetMarkerStyle(kFullCircle);
0613   g_AA_C10_20->SetMarkerStyle(kFullSquare);
0614   g_AA_C20_40->SetMarkerStyle(kFullDiamond);
0615   g_AA_C10_40->SetMarkerStyle(kFullDiamond);
0616   g_AA_C40_60->SetMarkerStyle(kFullCross);
0617 
0618 
0619   //
0620   TCanvas *c1 = new TCanvas("Draw_HFJetTruth_CrossSection2v2" + s_suffix, "Draw_HFJetTruth_CrossSection2v2" + s_suffix, 700, 600);
0621   c1->Divide(1, 1);
0622   int idx = 1;
0623   TPad *p;
0624 
0625   p = (TPad *) c1->cd(idx++);
0626   c1->Update();
0627 
0628   p->DrawFrame(15, -.1, 40, .3)
0629       ->SetTitle(";Transverse Momentum [GeV/#it{c}];v_{2}");
0630   //
0631   g_AA_C0_10->Draw("pe");
0632 //  g_AA_C10_20->Draw("pe");
0633 //  g_AA_C20_40->Draw("pe");
0634   g_AA_C10_40->Draw("pe");
0635 //  g_AA_C40_60->Draw("pe");
0636 
0637   //  g_AA_C20_40->Draw("same");
0638   //
0639   //  ge_RAA->SetLineWidth(3);
0640   //  ge_RAA->SetMarkerStyle(kFullCircle);
0641   //  ge_RAA->SetMarkerSize(2);
0642   //
0643   //  ge_RAA->Draw("pe");
0644   //  ge_RAA->Print();
0645   //
0646     TLegend *leg = new TLegend(.0, .78, .85, .93);
0647     leg->SetFillStyle(0);
0648     leg->AddEntry("", "#it{#bf{sPHENIX }} Simulation", "");
0649     leg->AddEntry("", Form("PYTHIA-8 #it{b}-jet, Anti-k_{T} R=0.4, |#eta|<%.2f, CTEQ6L", dy / 2), "");
0650     leg->AddEntry("", Form("Au+Au: %.0fB col., %.0f%% Eff., %.0f%% Pur.", '%', AuAu_MB_Evt / 1e9, AuAu_eff * 100, AuAu_purity * 100), "");
0651     leg->Draw();
0652   //
0653     TLegend *leg2 = new TLegend(.19, .55, 1, .78);
0654     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));
0655     leg2->AddEntry(g_AA_C0_10, "Au+Au 0-10%C", "pl");
0656 //    leg2->AddEntry(g_AA_C10_20, "Au+Au 10-20%C", "pl");
0657 //    leg2->AddEntry(g_AA_C20_40, "Au+Au 20-40%C", "pl");
0658     leg2->AddEntry(g_AA_C10_40, "Au+Au 10-40%C", "pl");
0659 //    leg2->AddEntry(g_AA_C40_60, "Au+Au 40-60%C", "pl");
0660     leg2->SetFillStyle(0);
0661     leg2->Draw();
0662 
0663   SaveCanvas(c1, infile + "_" + TString(c1->GetName()), kTRUE);
0664 }
0665 
0666 TGraphErrors *GetRAA(TH1 *h_pp, TH1 *h_AA)
0667 {
0668   int n_bin = 0;
0669 
0670   double xs[1000] = {0};
0671   double ys[1000] = {0};
0672   double eys[1000] = {0};
0673 
0674   assert(h_pp);
0675   assert(h_AA);
0676   assert(h_pp->GetNbinsX() == h_AA->GetNbinsX());
0677 
0678   for (int i = 1; i <= h_pp->GetNbinsX(); ++i)
0679   {
0680     if (h_pp->GetBinError(i) > 0 && h_pp->GetBinError(i) < 1 && h_AA->GetBinError(i) > 0 && h_AA->GetBinError(i) < 1)
0681     {
0682       xs[n_bin] = h_pp->GetXaxis()->GetBinCenter(i);
0683 
0684       ys[n_bin] = h_AA->GetBinContent(i) / h_pp->GetBinContent(i);
0685 
0686       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));
0687 
0688       n_bin += 1;
0689     }
0690   }
0691 
0692   TGraphErrors *ge = new TGraphErrors(n_bin, xs, ys, NULL, eys);
0693   ge->SetName(TString("RAA_") + h_AA->GetName());
0694 
0695   return ge;
0696 }
0697 
0698 TH1 *CrossSection2RelUncert(const TH1F *h_cross,
0699                             const double suppression,
0700                             const double deta,
0701                             const double pp_quiv_int_lum)
0702 {
0703   assert(h_cross);
0704   TH1 *
0705       h_ratio = (TH1 *)
0706                     h_cross->Clone(TString(h_cross->GetName()) + Form("_copy%d", rand()));
0707 
0708   //convert to count per bin
0709   h_ratio->Scale(deta * h_ratio->GetXaxis()->GetBinWidth(0) * pp_quiv_int_lum * suppression);
0710   h_ratio->Rebin(5);
0711 
0712   for (int i = 1; i <= h_ratio->GetNbinsX(); ++i)
0713   {
0714     const double yield = h_ratio->GetBinContent(i);
0715 
0716     if (yield > 0)
0717     {
0718       h_ratio->SetBinContent(i, suppression);
0719 
0720       h_ratio->SetBinError(i, suppression / sqrt(yield));
0721     }
0722     else
0723     {
0724       h_ratio->SetBinContent(i, 0);
0725 
0726       h_ratio->SetBinError(i, 0);
0727     }
0728   }
0729 
0730   h_ratio->GetYaxis()->SetTitle("Relative Cross Section and Uncertainty");
0731 
0732   return h_ratio;
0733 }
0734 
0735 TGraph *CrossSection2v2Uncert(const TH1F *h_cross,
0736                            const double suppression,
0737                            const double deta,
0738                            const double pp_quiv_int_lum,
0739                            const double ep_resolution = 1,
0740                            const double pt_shift = 0
0741 )
0742 {
0743   assert(h_cross);
0744   TH1 *
0745       h_ratio = (TH1 *)
0746                     h_cross->Clone(TString(h_cross->GetName()) + Form("_copyv2%d", rand()));
0747 
0748   //convert to count per bin
0749   h_ratio->Scale(deta * h_ratio->GetXaxis()->GetBinWidth(0) * pp_quiv_int_lum * suppression);
0750   h_ratio->Rebin(5);
0751 
0752   vector<double> pts;
0753   vector<double> v2s;
0754   vector<double> v2es;
0755 
0756   for (int i = 1; i <= h_ratio->GetNbinsX(); ++i)
0757   {
0758     const double yield = h_ratio->GetBinContent(i);
0759 
0760     if (yield > 100)
0761     {
0762       h_ratio->SetBinContent(i, 0);
0763 
0764       h_ratio->SetBinError(i, 1. / sqrt(2 * yield) / ep_resolution);  // err(v2) = 1/ (sqrt(2) *Significance * Resolution)
0765 
0766       pts.push_back(h_ratio->GetBinCenter(i) + pt_shift);
0767       v2s.push_back(h_ratio->GetBinContent(i));
0768       v2es.push_back(h_ratio->GetBinError(i));
0769     }
0770     else
0771     {
0772       h_ratio->SetBinContent(i, 0);
0773 
0774       h_ratio->SetBinError(i, 0);
0775     }
0776   }
0777 
0778   h_ratio->GetYaxis()->SetTitle("v2 and uncertainty");
0779 
0780   TGraph *gr = new TGraphErrors(pts.size(), &pts[0], &v2s[0], 0, &v2es[0]);
0781   gr->SetName(TString("ge_") + h_ratio->GetName());
0782 
0783   gr->SetLineWidth(3);
0784   gr->SetMarkerStyle(kFullCircle);
0785   gr->SetMarkerSize(2);
0786 
0787   return gr;
0788 }
0789 
0790 void Convert2CrossSection(TH1 *h, const double int_lumi, const double dy)
0791 {
0792   h->Sumw2();
0793 
0794   for (int i = 1; i <= h->GetXaxis()->GetNbins(); ++i)
0795   {
0796     const double pT1 = h->GetXaxis()->GetBinLowEdge(i);
0797     const double pT2 = h->GetXaxis()->GetBinUpEdge(i);
0798 
0799     //      const double dpT2 =  pT2*pT2 - pT1*pT1;
0800     const double dpT = pT2 - pT1;
0801 
0802     //      const double count2sigma = 1./dpT2/dy/int_lumi;
0803     const double count2sigma = 1. / dpT / dy / int_lumi;
0804 
0805     h->SetBinContent(i, h->GetBinContent(i) * count2sigma);
0806     h->SetBinError(i, h->GetBinError(i) * count2sigma);
0807   }
0808 
0809   //  h->GetYaxis()->SetTitle("#sigma/(dp_{T}^{2} dy) (pb/(GeV/c)^{2})");
0810 
0811   h->SetMarkerStyle(kFullCircle);
0812 }
0813 
0814 TGraphAsymmErrors *
0815 GetFONLL_B()
0816 {
0817   //# Job started on: Wed Aug 24 04:28:34 CEST 2016 .
0818   //# FONLL heavy quark hadroproduction cross section, calculated on Wed Aug 24 04:28:42 CEST 2016
0819   //# FONLL version and perturbative order: ## FONLL v1.3.2 fonll [ds/dpt^2dy (pb/GeV^2)]
0820   //# quark = bottom
0821   //# final state = quark
0822   //# ebeam1 = 100, ebeam2 = 100
0823   //# PDF set = CTEQ6.6
0824   //# ptmin = 10
0825   //# ptmax = 30
0826   //# etamin = -.6
0827   //# etamax = .6
0828   //# Uncertainties from scales, masses combined quadratically
0829   //# cross section is ds/dpt (pb/GeV)
0830   //# pt      central      min       max       min_sc     max_sc     min_mass   max_mass
0831   //  10.0000 4.1645e+03 2.9134e+03 6.0916e+03 2.9634e+03 6.0553e+03 3.8146e+03 4.5365e+03
0832   //  11.0526 2.4830e+03 1.7438e+03 3.6234e+03 1.7667e+03 3.6072e+03 2.3005e+03 2.6743e+03
0833   //  12.1053 1.5046e+03 1.0602e+03 2.1900e+03 1.0710e+03 2.1826e+03 1.4071e+03 1.6056e+03
0834   //  13.1579 9.2707e+02 6.5506e+02 1.3454e+03 6.6029e+02 1.3418e+03 8.7394e+02 9.8137e+02
0835   //  14.2105 5.8027e+02 4.1098e+02 8.3990e+02 4.1358e+02 8.3817e+02 5.5067e+02 6.1015e+02
0836   //  15.2632 3.6865e+02 2.6159e+02 5.3224e+02 2.6290e+02 5.3139e+02 3.5193e+02 3.8532e+02
0837   //  16.3158 2.3744e+02 1.6874e+02 3.4198e+02 1.6941e+02 3.4156e+02 2.2785e+02 2.4686e+02
0838   //  17.3684 1.5484e+02 1.1016e+02 2.2257e+02 1.1051e+02 2.2236e+02 1.4930e+02 1.6021e+02
0839   //  18.4211 1.0215e+02 7.2725e+01 1.4657e+02 7.2904e+01 1.4646e+02 9.8911e+01 1.0522e+02
0840   //  19.4737 6.8101e+01 4.8502e+01 9.7587e+01 4.8594e+01 9.7535e+01 6.6205e+01 6.9866e+01
0841   //  20.5263 4.5818e+01 3.2631e+01 6.5588e+01 3.2678e+01 6.5562e+01 4.4707e+01 4.6827e+01
0842   //  21.5789 3.1111e+01 2.2148e+01 4.4506e+01 2.2171e+01 4.4494e+01 3.0461e+01 3.1684e+01
0843   //  22.6316 2.1287e+01 1.5142e+01 3.0444e+01 1.5153e+01 3.0439e+01 2.0909e+01 2.1607e+01
0844   //  23.6842 1.4663e+01 1.0419e+01 2.0974e+01 1.0424e+01 2.0972e+01 1.4446e+01 1.4838e+01
0845   //  24.7368 1.0175e+01 7.2188e+00 1.4562e+01 7.2214e+00 1.4562e+01 1.0053e+01 1.0266e+01
0846   //  25.7895 7.1026e+00 5.0296e+00 1.0176e+01 5.0307e+00 1.0175e+01 7.0360e+00 7.1470e+00
0847   //  26.8421 4.9796e+00 3.5183e+00 7.1439e+00 3.5187e+00 7.1438e+00 4.9454e+00 4.9976e+00
0848   //  27.8947 3.5103e+00 2.4738e+00 5.0451e+00 2.4739e+00 5.0451e+00 3.4945e+00 3.5143e+00
0849   //  28.9474 2.4892e+00 1.7492e+00 3.5856e+00 1.7493e+00 3.5856e+00 2.4835e+00 2.4892e+00
0850   //  30.0000 1.7693e+00 1.2394e+00 2.5556e+00 1.2394e+00 2.5556e+00 1.7633e+00 1.7693e+00
0851 
0852   const double deta = 0.6 * 2;
0853 
0854   const double pts[] =
0855       {
0856 
0857           10.0000, 11.0526, 12.1053, 13.1579, 14.2105, 15.2632, 16.3158, 17.3684,
0858           18.4211, 19.4737, 20.5263, 21.5789, 22.6316, 23.6842, 24.7368, 25.7895,
0859           26.8421, 27.8947, 28.9474, 30.0000
0860 
0861       };
0862 
0863   const double centrals[] =
0864       {4.1645e+03, 2.4830e+03, 1.5046e+03, 9.2707e+02, 5.8027e+02, 3.6865e+02,
0865        2.3744e+02, 1.5484e+02, 1.0215e+02, 6.8101e+01, 4.5818e+01, 3.1111e+01,
0866        2.1287e+01, 1.4663e+01, 1.0175e+01, 7.1026e+00, 4.9796e+00, 3.5103e+00,
0867        2.4892e+00, 1.7693e+00};
0868   const double min[] =
0869       {4.1645e+03 - 2.9134e+03, 2.4830e+03 - 1.7438e+03, 1.5046e+03 - 1.0602e+03,
0870        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,
0871        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,
0872        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,
0873        2.4892e+00 - 1.7492e+00, 1.7693e+00 - 1.2394e+00};
0874 
0875   const double max[] =
0876       {6.0916e+03 - 4.1645e+03, 3.6234e+03 - 2.4830e+03, 2.1900e+03 - 1.5046e+03,
0877        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,
0878        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,
0879        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,
0880        3.5856e+00 - 2.4892e+00, 2.5556e+00 - 1.7693e+00};
0881 
0882   TGraphAsymmErrors *gr = new TGraphAsymmErrors(20, pts, centrals, 0, 0, min,
0883                                                 max);
0884 
0885   for (int i = 0; i < gr->GetN(); ++i)
0886   {
0887     (gr->GetY())[i] /= deta;
0888     (gr->GetEYhigh())[i] /= deta;
0889     (gr->GetEYlow())[i] /= deta;
0890   }
0891 
0892   return gr;
0893 }
0894 
0895 TGraphAsymmErrors *
0896 GetFONLL_C()
0897 {
0898   //# Job started on: Thu Aug 25 05:08:57 CEST 2016 .
0899   //# FONLL heavy quark hadroproduction cross section, calculated on Thu Aug 25 05:09:06 CEST 2016
0900   //# FONLL version and perturbative order: ## FONLL v1.3.2 fonll [ds/dpt^2dy (pb/GeV^2)]
0901   //# quark = charm
0902   //# final state = quark
0903   //# ebeam1 = 100, ebeam2 = 100
0904   //# PDF set = CTEQ6.6
0905   //# ptmin = 10
0906   //# ptmax = 30
0907   //# etamin = -.6
0908   //# etamax = .6
0909   //# Uncertainties from scales, masses combined quadratically
0910   //# cross section is ds/dpt (pb/GeV)
0911   //# pt      central      min       max       min_sc     max_sc     min_mass   max_mass
0912   //  10.0000 9.7368e+03 6.5051e+03 1.5528e+04 6.5641e+03 1.5522e+04 9.1225e+03 1.0006e+04
0913   //  11.0526 4.9069e+03 3.2787e+03 7.7308e+03 3.3256e+03 7.7226e+03 4.5190e+03 5.1221e+03
0914   //  12.1053 2.5743e+03 1.7187e+03 4.0147e+03 1.7525e+03 4.0069e+03 2.3363e+03 2.7237e+03
0915   //  13.1579 1.3990e+03 9.3274e+02 2.1641e+03 9.5597e+02 2.1578e+03 1.2537e+03 1.4976e+03
0916   //  14.2105 7.8344e+02 5.2138e+02 1.2037e+03 5.3701e+02 1.1989e+03 6.9430e+02 8.4730e+02
0917   //  15.2632 4.5087e+02 2.9938e+02 6.8902e+02 3.0978e+02 6.8543e+02 3.9570e+02 4.9204e+02
0918   //  16.3158 2.6550e+02 1.7584e+02 4.0398e+02 1.8275e+02 4.0140e+02 2.3098e+02 2.9210e+02
0919   //  17.3684 1.5978e+02 1.0551e+02 2.4232e+02 1.1011e+02 2.4049e+02 1.3793e+02 1.7705e+02
0920   //  18.4211 9.7903e+01 6.4459e+01 1.4811e+02 6.7524e+01 1.4682e+02 8.3915e+01 1.0919e+02
0921   //  19.4737 6.0993e+01 4.0030e+01 9.2103e+01 4.2085e+01 9.1200e+01 5.1941e+01 6.8435e+01
0922   //  20.5263 3.8605e+01 2.5248e+01 5.8239e+01 2.6631e+01 5.7606e+01 3.2685e+01 4.3550e+01
0923   //  21.5789 2.4727e+01 1.6108e+01 3.7286e+01 1.7044e+01 3.6842e+01 2.0821e+01 2.8038e+01
0924   //  22.6316 1.6036e+01 1.0402e+01 2.4185e+01 1.1040e+01 2.3872e+01 1.3433e+01 1.8271e+01
0925   //  23.6842 1.0532e+01 6.8006e+00 1.5898e+01 7.2389e+00 1.5678e+01 8.7774e+00 1.2053e+01
0926   //  24.7368 6.9729e+00 4.4803e+00 1.0540e+01 4.7838e+00 1.0385e+01 5.7809e+00 8.0136e+00
0927   //  25.7895 4.6560e+00 2.9760e+00 7.0510e+00 3.1876e+00 6.9411e+00 3.8397e+00 5.3729e+00
0928   //  26.8421 3.1409e+00 1.9970e+00 4.7681e+00 2.1454e+00 4.6902e+00 2.5775e+00 3.6384e+00
0929   //  27.8947 2.1346e+00 1.3500e+00 3.2500e+00 1.4543e+00 3.1946e+00 1.7438e+00 2.4816e+00
0930   //  28.9474 1.4557e+00 9.1535e-01 2.2239e+00 9.8868e-01 2.1845e+00 1.1839e+00 1.6985e+00
0931   //  30.0000 9.9826e-01 6.2400e-01 1.5310e+00 6.7571e-01 1.5029e+00 8.0844e-01 1.1690e+00
0932 
0933   const double deta = 0.6 * 2;
0934 
0935   const double pts[] =
0936       {
0937           //         pt
0938           10.0000, 11.0526, 12.1053, 13.1579, 14.2105, 15.2632, 16.3158, 17.3684,
0939           18.4211, 19.4737, 20.5263, 21.5789, 22.6316, 23.6842, 24.7368, 25.7895,
0940           26.8421, 27.8947, 28.9474, 30.0000
0941 
0942       };
0943 
0944   const double centrals[] =
0945       {
0946 
0947           //            central
0948           9.7368e+03, 4.9069e+03, 2.5743e+03, 1.3990e+03, 7.8344e+02, 4.5087e+02,
0949           2.6550e+02, 1.5978e+02, 9.7903e+01, 6.0993e+01, 3.8605e+01, 2.4727e+01,
0950           1.6036e+01, 1.0532e+01, 6.9729e+00, 4.6560e+00, 3.1409e+00, 2.1346e+00,
0951           1.4557e+00, 9.9826e-01
0952 
0953       };
0954   const double min[] =
0955       {
0956 
0957           //       central           min
0958           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,
0959           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,
0960           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,
0961           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
0962 
0963       };
0964 
0965   const double max[] =
0966       {
0967 
0968           //            max             central
0969           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,
0970           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,
0971           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,
0972           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
0973 
0974       };
0975 
0976   TGraphAsymmErrors *gr = new TGraphAsymmErrors(20, pts, centrals, 0, 0, min,
0977                                                 max);
0978 
0979   for (int i = 0; i < gr->GetN(); ++i)
0980   {
0981     (gr->GetY())[i] /= deta;
0982     (gr->GetEYhigh())[i] /= deta;
0983     (gr->GetEYlow())[i] /= deta;
0984   }
0985 
0986   return gr;
0987 }
0988 
0989 TGraph *pQCDModel_HuangKangVitev(const double g)
0990 {
0991   // arXiv:1306.0909v2 [hep-ph]
0992   // Preliminary for sPHENIX energy
0993 
0994   //  R = 0.4;
0995   if (g == 2.0)
0996   {
0997     //  g=2.0
0998     //  10.130739 0.8087402
0999     //  13.8755   0.7102165
1000     //  17.915709 0.64113617
1001     //  19.950817 0.62270004
1002     //  21.838287 0.59506375
1003     //  24.758118 0.5618901
1004     //  28.149998 0.53331053
1005     //  30.303171 0.5194738
1006     //  35.701473 0.52399224
1007     //  40.21508  0.54508865
1008     //  49.182514 0.53758913
1009     //  54.492188 0.5338267
1010     //  59.890347 0.52914274
1011 
1012     const double pT[] =
1013         {
1014             10.130739,
1015             13.8755,
1016             17.915709,
1017             19.950817,
1018             21.838287,
1019             24.758118,
1020             28.149998,
1021             30.303171,
1022             35.701473,
1023             40.21508,
1024             49.182514,
1025             54.492188,
1026             59.890347};
1027 
1028     const double RAA[] =
1029         {
1030             0.8087402,
1031             0.7102165,
1032             0.64113617,
1033             0.62270004,
1034             0.59506375,
1035             0.5618901,
1036             0.53331053,
1037             0.5194738,
1038             0.52399224,
1039             0.54508865,
1040             0.53758913,
1041             0.5338267,
1042             0.52914274};
1043 
1044     return new TGraph(13, pT, RAA);
1045   }
1046   else if (g == 2.2)
1047   {
1048     //  g=2.2
1049     //  10.040924 0.72499925
1050     //  11.750852 0.6623964
1051     //  15.820038 0.56018674
1052     //  18.946156 0.51412654
1053     //  19.860464 0.50491005
1054     //  21.010588 0.484647
1055     //  24.343283 0.44410512
1056     //  29.888279 0.39800778
1057     //  35.876522 0.4006767
1058     //  40.154125 0.42085648
1059     //  47.941647 0.41521555
1060     //  55.64066  0.40865576
1061     //  59.91793  0.4076699
1062 
1063     const double pT[] =
1064         {
1065             10.040924,
1066             11.750852,
1067             15.820038,
1068             18.946156,
1069             19.860464,
1070             21.010588,
1071             24.343283,
1072             29.888279,
1073             35.876522,
1074             40.154125,
1075             47.941647,
1076             55.64066,
1077             59.91793};
1078 
1079     const double RAA[] =
1080         {
1081             0.72499925,
1082             0.6623964,
1083             0.56018674,
1084             0.51412654,
1085             0.50491005,
1086             0.484647,
1087             0.44410512,
1088             0.39800778,
1089             0.4006767,
1090             0.42085648,
1091             0.41521555,
1092             0.40865576,
1093             0.4076699};
1094 
1095     return new TGraph(13, pT, RAA);
1096   }
1097   else
1098   {
1099     assert(0);
1100   }
1101 }
1102 
1103 TGraph *
1104 GetPHENIX_jet()
1105 {
1106   //  PPG184
1107   //  p+p d^2sigma/dpT deta [mb / GeV]
1108   //
1109   //  pT low     pT high     yval        ystat[%] ysyst[%]
1110   //
1111   //  12.2       14.5        0.0001145   1.0      15.6
1112   //  14.5       17.3        3.774e-05   1.3      16.3
1113   //  17.3       20.7        1.263e-05   1.7      15.6
1114   //  20.7       24.7        4.230e-06   2.3      17.3
1115   //  24.7       29.4        1.224e-06   3.6      19.3
1116   //  29.4       35.1        2.962e-07   6.1      21.0
1117   //  35.1       41.9        5.837e-08   8.9      24.3
1118   //  41.9       50.0        8.882e-09   11       32.3
1119 
1120   const double pt_l[] =
1121       {
1122 
1123           //          pT low
1124 
1125           12.2, 14.5, 17.3, 20.7, 24.7, 29.4, 35.1, 41.9
1126 
1127       };
1128 
1129   const double pt_h[] =
1130       {
1131 
1132           //              pT high
1133 
1134           14.5, 17.3, 20.7, 24.7, 29.4, 35.1, 41.9, 50.0
1135 
1136       };
1137   const double yval[] =
1138       {
1139 
1140           //          yval
1141 
1142           0.0001145, 3.774e-05, 1.263e-05, 4.230e-06, 1.224e-06, 2.962e-07,
1143           5.837e-08, 8.882e-09
1144 
1145       };
1146   const double ystat[] =
1147       {
1148 
1149           //          ystat[%]
1150 
1151           1.0, 1.3, 1.7, 2.3, 3.6, 6.1, 8.9, 11
1152 
1153       };
1154 
1155   const double ysyst[] =
1156       {
1157 
1158           //            ysyst[%]
1159 
1160           15.6, 16.3, 15.6, 17.3, 19.3, 21.0, 24.3, 32.3
1161 
1162       };
1163 
1164   double pT_c[100] =
1165       {0};
1166   double pT_e[100] =
1167       {0};
1168   double y_e[100] =
1169       {0};
1170 
1171   for (int i = 0; i < 8; ++i)
1172   {
1173     pT_c[i] = 0.5 * (pt_l[i] + pt_h[i]);
1174     pT_e[i] = 0.5 * (pt_h[i] - pt_l[i]);
1175     yval[i] *= 1e9;  // mb -> pb
1176     y_e[i] = yval[i] * sqrt(ystat[i] * ystat[i] + ysyst[i] * ysyst[i]) / 100;
1177   }
1178 
1179   TGraphErrors *gr = new TGraphErrors(8, pT_c, yval, pT_e, y_e);
1180 
1181   return gr;
1182 }
[ 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!