Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Draw_HFJetTruth.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Draw_HFJetTruth.C
1 // $Id: $
2 
11 #include <TFile.h>
12 #include <TGraph.h>
13 #include <TGraphAsymmErrors.h>
14 #include <TGraphErrors.h>
15 #include <TLatex.h>
16 #include <TLine.h>
17 #include <TString.h>
18 #include <TTree.h>
19 #include <cassert>
20 #include <cmath>
21 #include "SaveCanvas.C"
22 #include "sPhenixStyle.C"
23 
24 TFile *_file0 = NULL;
25 TTree *T = NULL;
26 
27 void Draw_HFJetTruth(const TString infile =
28  // "/sphenix/user/jinhuang/HF-jet/event_gen/200pp_pythia8_CTEQ6L_20GeV/200pp_pythia8_CTEQ6L_20GeV_ALL.cfg_eneg_DSTReader.root",
29  // double int_lumi = 210715 / 5.533e-05 / 1e9, const double dy = 0.6 * 2)
30  "/sphenix/user/jinhuang/HF-jet/event_gen/200pp_pythia8_CTEQ6L_7GeV/200pp_pythia8_CTEQ6L_7GeV_ALL.cfg_eneg_DSTReader.root",
31  double int_lumi = 891093 / 3.332e-02 / 1e9, const double dy = 0.6 * 2)
32 //"/sphenix/user/jinhuang/HF-jet/event_gen/200pp_pythia8_CTEQ6L/200pp_pythia8_CTEQ6L_111ALL.cfg_eneg_DSTReader.root",
33 //double int_lumi = 789908/4.631e-03 / 1e9, const double dy = 0.6*2)
34 //Draw_HFJetTruth(const TString infile =
35 // "../macros3/G4sPHENIXCells.root_DSTReader.root",
36 // double int_lumi = 789908/4.631e-03 / 1e9, const double dy = 0.6*2)
37 //Draw_HFJetTruth(const TString infile =
38 // "../macros2/G4sPHENIXCells.root_DSTReader.root",
39 // const double int_lumi = 1842152 / 5.801e-03 / 1e9,
40 // const double dy = 0.6 * 2)
41 {
42  SetsPhenixStyle();
43  gStyle->SetOptStat(0);
44  gStyle->SetOptFit(1111);
45  TVirtualFitter::SetDefaultFitter("Minuit2");
46 
47  gSystem->Load("libg4eval.so");
48 
49  if (!_file0)
50  {
51  TString chian_str = infile;
52  chian_str.ReplaceAll("ALL", "*");
53 
54  TChain *t = new TChain("T");
55  const int n = t->Add(chian_str);
56 
57  int_lumi *= n;
58  // cout << "Loaded " << n << " root files with " << chian_str << endl;
59  cout << "int_lumi = " << int_lumi << " pb^-1 from " << n << " root files with " << chian_str << endl;
60  assert(n > 0);
61 
62  T = t;
63 
64  _file0 = new TFile;
65  _file0->SetName(infile);
66  }
67 
68 // DrawCrossSection(int_lumi, dy);
69  // Draw_HFJetTruth_DrawCrossSection_PR(infile);
70  // CrossSection2RAA_Proposal(infile);
71  // CrossSection2RAA(infile);
72  const double acceptance1 = 2.* (1.1 - .4);
73  CrossSection2RAA(infile, false, acceptance1);
74 
75 // const double acceptance2 = 2.* (0.85 - .4);
76 // CrossSection2RAA(infile, false, acceptance2);
77 
78 // CrossSection2v2(infile, false, .7, acceptance);
79 // CrossSection2v2(infile, false, .4, acceptance);
80 }
81 
82 void DrawCrossSection(double int_lumi, const double dy)
83 {
84  TH1 *hall = new TH1F("hall", ";p_{T} [GeV/c]", 100, 0, 100);
85  TH1 *h_b = new TH1F("h_b", ";p_{T} [GeV/c]", 100, 0, 100);
86  TH1 *h_bq = new TH1F("h_bq", ";p_{T} [GeV/c]", 100, 0, 100);
87  TH1 *h_bh = new TH1F("h_bh", ";p_{T} [GeV/c]", 100, 0, 100);
88  TH1 *h_bh5 = new TH1F("h_bh5", ";p_{T} [GeV/c]", 100, 0, 100);
89  TH1 *h_ch5 = new TH1F("h_ch5", ";p_{T} [GeV/c]", 100, 0, 100);
90  TH1 *h_c = new TH1F("h_c", ";p_{T} [GeV/c]", 100, 0, 100);
91 
92  T->Draw("AntiKt_Truth_r04.get_pt()>>hall",
93  "AntiKt_Truth_r04.get_pt()>15 && abs(AntiKt_Truth_r04.get_eta())<0.6",
94  "goff");
95 
96  T->Draw("AntiKt_Truth_r04.get_pt()>>h_b",
97  "AntiKt_Truth_r04.get_pt()>15 && abs(AntiKt_Truth_r04.get_eta())<0.6 && abs(AntiKt_Truth_r04.get_property(1000))==5",
98  "goff");
99  T->Draw(
100  "AntiKt_Truth_r04.get_pt()* AntiKt_Truth_r04.get_property(1001) >>h_bq",
101  "AntiKt_Truth_r04.get_pt()>15 && abs(AntiKt_Truth_r04.get_eta())<0.6 && abs(AntiKt_Truth_r04.get_property(1000))==5",
102  "goff");
103 
104  // T->Draw("AntiKt_Truth_r04.get_pt()>>h_bh",
105  // "AntiKt_Truth_r04.get_pt()>15 && abs(AntiKt_Truth_r04.get_eta())<0.6 && AntiKt_Truth_r04.get_property(1010)==5",
106  // "goff");
107 
108  T->Draw("AntiKt_Truth_r04.get_pt()>>h_bh5",
109  "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",
110  "goff");
111  //
112  T->Draw("AntiKt_Truth_r04.get_pt()>>h_c",
113  "AntiKt_Truth_r04.get_pt()>15 && abs(AntiKt_Truth_r04.get_eta())<0.6 && abs(AntiKt_Truth_r04.get_property(1000))==4",
114  "goff");
115  T->Draw("AntiKt_Truth_r04.get_pt()>>h_ch5",
116  "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",
117  "goff");
118 
119  Convert2CrossSection(hall, int_lumi, dy);
120  Convert2CrossSection(h_b, int_lumi, dy);
121  Convert2CrossSection(h_bq, int_lumi, dy);
122  Convert2CrossSection(h_bh5, int_lumi, dy);
123  Convert2CrossSection(h_ch5, int_lumi, dy);
124  Convert2CrossSection(h_c, int_lumi, dy);
125 
126  h_b->SetLineColor(kBlue);
127  h_b->SetMarkerColor(kBlue);
128 
129  h_bq->SetLineColor(kBlue + 3);
130  h_bq->SetMarkerColor(kBlue + 3);
131 
132  h_bh5->SetLineColor(kBlue + 3);
133  h_bh5->SetMarkerColor(kBlue + 3);
134  h_bh5->SetMarkerStyle(kStar);
135 
136  h_c->SetLineColor(kRed);
137  h_c->SetMarkerColor(kRed);
138 
139  h_ch5->SetLineColor(kRed + 3);
140  h_ch5->SetMarkerColor(kRed + 3);
141  h_ch5->SetMarkerStyle(kStar);
142 
143  TGraphAsymmErrors *gr_fonll_b = GetFONLL_B();
144  gr_fonll_b->SetFillColor(kBlue - 7);
145  gr_fonll_b->SetFillStyle(3002);
146  TGraphAsymmErrors *gr_fonll_c = GetFONLL_C();
147  gr_fonll_c->SetFillColor(kRed - 7);
148  gr_fonll_c->SetFillStyle(3003);
149  TGraph *gr_phenix = GetPHENIX_jet();
150  gr_phenix->SetLineColor(kGray + 2);
151  gr_phenix->SetMarkerColor(kBlack);
152  gr_phenix->SetMarkerStyle(kOpenCross);
153 
154  TCanvas *c1 = new TCanvas("Draw_HFJetTruth_DrawCrossSection", "Draw_HFJetTruth_DrawCrossSection", 1000, 860);
155  // c1->Divide(2, 2);
156  int idx = 1;
157  TPad *p;
158 
159  p = (TPad *) c1->cd(idx++);
160  c1->Update();
161  p->SetGridx(0);
162  p->SetGridy(0);
163  p->SetLogy();
164 
165  hall->Draw();
166 
167  gr_fonll_b->Draw("3");
168  gr_fonll_c->Draw("3");
169  gr_phenix->Draw("pe");
170 
171  h_b->Draw("same");
172  h_bh5->Draw("same");
173  h_bq->Draw("same");
174  h_c->Draw("same");
175  h_ch5->Draw("same");
176 
177  hall->GetXaxis()->SetRangeUser(12, 60);
178  hall->GetYaxis()->SetTitle(
179  "d^{2}#sigma/(dp_{T}dy), d^{2}#sigma/(dp_{T}d#eta) [pb/(GeV/c)]");
180 
181  TLegend *leg = new TLegend(0.45, 0.6, 0.99, 0.95);
182  leg->SetFillColor(kWhite);
183  leg->SetFillStyle(1001);
184  leg->SetHeader("p+p collisions @ sPHENIX, #sqrt{s} = 200 GeV, |#eta|<0.6");
185  leg->AddEntry(hall, "Inclusive jet, Pythia8, Truth, anti-k_{t}, R=0.4",
186  "lpe");
187  leg->AddEntry(h_c, "c-quark jet, Pythia8, Truth, anti-k_{t}, R=0.4", "lpe");
188  leg->AddEntry(h_b, "b-quark jet, Pythia8, Truth, anti-k_{t}, R=0.4", "lpe");
189  leg->AddEntry(h_bq, "Leading b-quark in jet, Pythia8, anti-k_{t}, R=0.4", "lpe");
190  leg->AddEntry(h_bh5,
191  "b-hadron jet, Pythia8, Truth, anti-k_{t}, R=0.4, p_{T, b-hadron}>5 GeV/c",
192  "lpe");
193  leg->AddEntry(gr_phenix, "PHENIX inclusive jet, PRL 116, 122301 (2016)", "ple");
194  leg->AddEntry(gr_fonll_c, "c-quark, FONLL v1.3.2, CTEQ6.6, |y|<0.6", "f");
195  leg->AddEntry(gr_fonll_b, "b-quark, FONLL v1.3.2, CTEQ6.6, |y|<0.6", "f");
196  leg->Draw();
197 
198  SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
199 }
200 
201 void Draw_HFJetTruth_DrawCrossSection_PR(const TString infile)
202 {
203  TFile *f = TFile::Open(infile + "Draw_HFJetTruth_DrawCrossSection.root");
204  assert(f);
205 
206  TH1F *hall = (TH1F *) f->GetObjectChecked("hall", "TH1F");
207  assert(hall);
208  TH1F *h_b = (TH1F *) f->GetObjectChecked("h_b", "TH1F");
209  assert(h_b);
210 
211  hall->SetMarkerColor(kBlack);
212  hall->SetLineColor(kBlack);
213  hall->SetMarkerStyle(kFullCircle);
214 
215  h_b->SetMarkerColor(kBlue + 2);
216  h_b->SetLineColor(kBlue + 2);
217  h_b->SetMarkerStyle(kFullCircle);
218 
219  TGraphAsymmErrors *gr_fonll_b = GetFONLL_B();
220  gr_fonll_b->SetFillColor(kBlue - 7);
221  gr_fonll_b->SetFillStyle(1001);
222  TGraphAsymmErrors *gr_fonll_c = GetFONLL_C();
223  gr_fonll_c->SetFillColor(kRed - 7);
224  gr_fonll_c->SetFillStyle(3003);
225  TGraph *gr_phenix = GetPHENIX_jet();
226  gr_phenix->SetLineColor(kGray + 2);
227  gr_phenix->SetLineWidth(3);
228  gr_phenix->SetMarkerColor(kGray + 2);
229  gr_phenix->SetMarkerSize(2);
230  gr_phenix->SetMarkerStyle(kFullSquare);
231 
232  TCanvas *c1 = new TCanvas("Draw_HFJetTruth_DrawCrossSection_PR", "Draw_HFJetTruth_DrawCrossSection_PR", 1000, 860);
233  // c1->Divide(2, 2);
234  int idx = 1;
235  TPad *p;
236 
237  p = (TPad *) c1->cd(idx++);
238  c1->Update();
239  p->SetLogy();
240 
241  TH1 *hframe = p->DrawFrame(12, 0.1, 70, 1e8);
242  hframe->SetTitle(";p_{T} [GeV/c];d^{2}#sigma/(dp_{T}d#eta) [pb/(GeV/c)]");
243 
244  gr_fonll_b->Draw("3");
245  gr_phenix->Draw("pe");
246 
247  hall->Draw("same");
248  h_b->Draw("same");
249 
250  TLegend *leg = new TLegend(0.2, 0.7, 0.95, 0.92);
251  leg->SetFillColor(kWhite);
252  leg->SetFillStyle(1001);
253  // leg->SetHeader("#splitline{#it{#bf{sPHENIX }} Simulation}{p+p, #sqrt{s} = 200 GeV, |#eta|<0.6}");
254  leg->SetHeader("#it{#bf{sPHENIX }} Simulation, #it{p}+#it{p} #sqrt{s} = 200 GeV, |#eta|<0.6");
255  leg->AddEntry(hall, "Inclusive jet, PYTHIA8 + CTEQ6L, anti-k_{T} R=0.4",
256  "lpe");
257  leg->AddEntry(h_b, "#it{b}-quark jet, PYTHIA8 + CTEQ6L, anti-k_{T} R=0.4", "lpe");
258  leg->AddEntry(gr_phenix, "PHENIX inclusive jet, PRL 116, 122301 (2016)", "ple");
259  leg->AddEntry(gr_fonll_b, "#it{b}-quark, FONLL v1.3.2, CTEQ6.6", "f");
260  leg->Draw();
261 
262  SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
263 }
264 
265 void CrossSection2RAA_Proposal(const TString infile)
266 {
267  const double b_jet_RAA = 0.6;
268 
270  // 5-year lumi in [sPH-TRG-000]
272 
273  const double pp_inelastic_crosssec = 42e-3 / 1e-12; // 42 mb in pb [sPH-TRG-000]
274  const double AuAu_Ncoll_C0_10 = 960.2; // [DOI:?10.1103/PhysRevC.87.034911?]
275  const double AuAu_Ncoll_C0_20 = 770.6; // [DOI:?10.1103/PhysRevC.91.064904?]
276  const double AuAu_Ncoll_C0_100 = 250; // pb^-1 [sPH-TRG-000]
277  const double pAu_Ncoll_C0_100 = 4.7; // pb^-1 [sPH-TRG-000]
278 
280  // 2-year lumi in sPHENIX proposal
282 
283  // const double pp_lumi_proposal = 630; // Figure 4.36:
284  // const double AuAu_eq_lumi_C0_20_proposal = 0.6e12 * 0.2 * AuAu_Ncoll_C0_20 / pp_inelastic_crosssec; //
285  const double pp_lumi_proposal = 175; // Table 4.1
286  const double AuAu_eq_lumi_C0_20_proposal = 0.1e12 * 0.2 * AuAu_Ncoll_C0_20 / pp_inelastic_crosssec; //
287  const double dy_proposal = 2.;
288  const double eff_proposal = 0.5;
289  const double dy_proposal = 2;
290 
291  cout << "CrossSection2RAA_Proposal integrated luminosity assumptions in pb^-1: " << endl;
292  cout << "\t"
293  << "pp_lumi_proposal = " << pp_lumi_proposal << endl;
294  cout << "\t"
295  << "AuAu_eq_lumi_C0_20_proposal = " << AuAu_eq_lumi_C0_20_proposal << endl;
296 
297  TFile *f = TFile::Open(infile + "Draw_HFJetTruth_DrawCrossSection.root");
298  assert(f);
299 
300  TH1F *h_b = (TH1F *) f->GetObjectChecked("h_b", "TH1F");
301 
302  assert(h_b);
303 
304  TCanvas *c1 = new TCanvas("Draw_HFJetTruth_CrossSection2RAA_Proposal", "Draw_HFJetTruth_CrossSection2RAA_Proposal", 1000, 860);
305  c1->Divide(2, 2);
306  int idx = 1;
307  TPad *p;
308 
309  p = (TPad *) c1->cd(idx++);
310  c1->Update();
311  // p->SetGridx(0);
312  // p->SetGridy(0);
313  p->SetLogy();
314 
315  h_b->GetYaxis()->SetTitle(
316  "d^{2}#sigma/(dp_{T}d#eta) [pb/(GeV/c)]");
317  h_b->Draw();
318 
319  p = (TPad *) c1->cd(idx++);
320  c1->Update();
321  // p->SetGridx(0);
322  // p->SetGridy(0);
323  p->SetLogy();
324 
325  TH1 *h_b_int = (TH1 *) h_b->Clone(TString(h_b->GetName()) + "_IntrgratedCount");
326  double integral = 0;
327  for (int i = h_b_int->GetNbinsX(); i >= 1; --i)
328  {
329  const double cs = h_b_int->GetBinContent(i);
330 
331  integral += cs * h_b_int->GetXaxis()->GetBinWidth(i) * dy_proposal * 0.2 * AuAu_Ncoll_C0_20 / pp_inelastic_crosssec;
332 
333  const double cs = h_b_int->SetBinContent(i, integral);
334  const double cs = h_b_int->SetBinError(i, 0);
335  }
336 
337  h_b_int->GetYaxis()->SetTitle(
338  "(Count > p_{T} cut)/Event 0-20% AuAu");
339  h_b_int->Draw();
340 
341  p = (TPad *) c1->cd(idx++);
342  c1->Update();
343  // p->SetGridx(0);
344  // p->SetGridy(0);
345  // p->SetLogy();
346 
347  TH1 *g_pp = CrossSection2RelUncert(h_b, 1., dy_proposal, pp_lumi_proposal * eff_proposal);
348  TH1 *g_AA = CrossSection2RelUncert(h_b, b_jet_RAA, dy_proposal, AuAu_eq_lumi_C0_20_proposal * eff_proposal);
349 
350  g_pp->Draw();
351  g_AA->Draw("same");
352 
353  p = (TPad *) c1->cd(idx++);
354  c1->Update();
355  // p->SetGridx(0);
356  // p->SetGridy(0);
357  // p->SetLogy();
358 
359  p->DrawFrame(0, 0, 100, 1.2)
360  ->SetTitle(";Transverse Momentum [GeV/#it{c}];#it{R}_{AA}");
361  TLatex t;
362  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));
363 
364  TGraphErrors *ge_RAA = GetRAA(g_pp, g_AA);
365  ge_RAA->Draw("pe*");
366 
367  SaveCanvas(c1, infile + "_" + TString(c1->GetName()), kTRUE);
368 }
369 
370 void CrossSection2RAA(const TString infile, const bool use_AA_jet_trigger = true, const double dy = .7 * 2)
371 {
372  TFile *f = TFile::Open(infile + "Draw_HFJetTruth_DrawCrossSection.root");
373  assert(f);
374 
375  TH1F *hall = (TH1F *) f->GetObjectChecked("hall", "TH1F");
376  assert(hall);
377  TH1F *h_b = (TH1F *) f->GetObjectChecked("h_b", "TH1F");
378  assert(h_b);
379 
380  const TString s_suffix(use_AA_jet_trigger ? "_AAJetTriggered" : "");
381  s_suffix += Form("_deta%.2f",dy/2);
382 
383  const double b_jet_RAA = 0.6;
384 
385  const double pp_eff = 0.6;
386  const double pp_purity = 0.4;
387  const double AuAu_eff = 0.4;
388  const double AuAu_purity = 0.4;
389 
391  // 5-year lumi in [sPH-TRG-000]
393 
394  const double pp_lumi = 200; // pb^-1 [sPH-TRG-000], rounded up from 197 pb^-1
395  const double pp_inelastic_crosssec = 42e-3 / 1e-12; // 42 mb in pb [sPH-TRG-000]
396 
397  const double AuAu_MB_Evt = use_AA_jet_trigger ? 550e9 : 240e9; // [sPH-TRG-000], depending on whether jet trigger applied in AA collisions
398  const double pAu_MB_Evt = 600e9; // [sPH-TRG-000]
399 
400  const double AuAu_Ncoll_C0_10 = 960.2; // [DOI:?10.1103/PhysRevC.87.034911?]
401  const double AuAu_Ncoll_C0_20 = 770.6; // [DOI:?10.1103/PhysRevC.91.064904?]
402  const double AuAu_Ncoll_C0_100 = 250; // pb^-1 [sPH-TRG-000]
403  const double pAu_Ncoll_C0_100 = 4.7; // pb^-1 [sPH-TRG-000]
404 
405  const double AuAu_eq_lumi_C0_10 = AuAu_MB_Evt * 0.1 * AuAu_Ncoll_C0_10 / pp_inelastic_crosssec; //
406  const double AuAu_eq_lumi_C0_20 = AuAu_MB_Evt * 0.2 * AuAu_Ncoll_C0_20 / pp_inelastic_crosssec; //
407  const double AuAu_eq_lumi_C0_100 = AuAu_MB_Evt * 1 * AuAu_Ncoll_C0_100 / pp_inelastic_crosssec; //
408 
409  const double pAu_eq_lumi_C0_100 = pAu_MB_Evt * 1 * pAu_Ncoll_C0_100 / pp_inelastic_crosssec; //
410 
411  cout << "CrossSection2RAA integrated luminosity assumptions in pb^-1: " << endl;
412  cout << "\t"
413  << "pp_lumi = " << pp_lumi << endl;
414  cout << "\t"
415  << "AuAu_eq_lumi_C0_10 = " << AuAu_eq_lumi_C0_10 << endl;
416  cout << "\t"
417  << "AuAu_eq_lumi_C0_20 = " << AuAu_eq_lumi_C0_20 << endl;
418  cout << "\t"
419  << "AuAu_eq_lumi_C0_100 = " << AuAu_eq_lumi_C0_100 << endl;
420  cout << "\t"
421  << "pAu_eq_lumi_C0_100 = " << pAu_eq_lumi_C0_100 << endl;
422 
423  TCanvas *c1 = new TCanvas("Draw_HFJetTruth_CrossSection2RAA_Ratio" + s_suffix, "Draw_HFJetTruth_CrossSection2RAA_Ratio" + s_suffix, 700, 600);
424  c1->Divide(1, 1);
425  int idx = 1;
426  TPad *p;
427 
428  p = (TPad *) c1->cd(idx++);
429  c1->Update();
430  // p->SetGridx(0);
431  // p->SetGridy(0);
432  // p->SetLogy();
433 
434  TH1 *g_pp = CrossSection2RelUncert(h_b, 1., dy, pp_lumi * pp_eff * pp_purity);
435  TH1 *g_AA = CrossSection2RelUncert(h_b, b_jet_RAA, dy, AuAu_eq_lumi_C0_10 * AuAu_eff * AuAu_purity);
436 
437  g_pp->SetLineColor(kRed);
438  g_AA->SetLineColor(kBlue);
439 
440  g_pp->Draw();
441  g_AA->Draw("same");
442  SaveCanvas(c1, infile + "_" + TString(c1->GetName()), kTRUE);
443 
444  TCanvas *c1 = new TCanvas("Draw_HFJetTruth_CrossSection2RAA" + s_suffix, "Draw_HFJetTruth_CrossSection2RAA" + s_suffix, 700, 600);
445  c1->Divide(1, 1);
446  int idx = 1;
447  TPad *p;
448 
449  p = (TPad *) c1->cd(idx++);
450  c1->Update();
451 
452  p->DrawFrame(11, 0, 50, 1.2)
453  ->SetTitle(";Transverse Momentum [GeV/#it{c}];#it{R}_{AA}");
454 
455  TGraphErrors *ge_RAA = GetRAA(g_pp, g_AA);
456 
457  ge_RAA->SetLineWidth(3);
458  ge_RAA->SetMarkerStyle(kFullCircle);
459  ge_RAA->SetMarkerSize(2);
460 
461  ge_RAA->Draw("pe");
462  ge_RAA->Print();
463 
464  TLegend *leg = new TLegend(.0, .76, .85, .93);
465  leg->SetFillStyle(0);
466  leg->AddEntry("", "#it{#bf{sPHENIX }} Simulation", "");
467  leg->AddEntry("", Form("PYTHIA-8 #it{b}-jet, Anti-k_{T} R=0.4, |#eta|<%.2f, CTEQ6L", dy / 2), "");
468  leg->AddEntry("", Form("#it{p}+#it{p}: %.0f pb^{-1}, %.0f%% Eff., %.0f%% Pur.", pp_lumi, pp_eff * 100, pp_purity * 100), "");
469  leg->AddEntry("", Form("Au+Au: %.0fB col., %.0f%% Eff., %.0f%% Pur.", '%', AuAu_MB_Evt / 1e9, AuAu_eff * 100, AuAu_purity * 100), "");
470  leg->Draw();
471 
472  TLegend *leg2 = new TLegend(.17, .70, .88, .77);
473  leg2->AddEntry(ge_RAA, "#it{b}-jet #it{R}_{AA}, Au+Au 0-10%C, #sqrt{s_{NN}}=200 GeV", "pl");
474  leg2->Draw();
475 
476  SaveCanvas(c1, infile + "_" + TString(c1->GetName()), kTRUE);
477 
478  TCanvas *c1 = new TCanvas("Draw_HFJetTruth_CrossSection2RAA_Theory" + s_suffix, "Draw_HFJetTruth_CrossSection2RAA_Theory" + s_suffix, 700, 600);
479  c1->Divide(1, 1);
480  int idx = 1;
481  TPad *p;
482 
483  p = (TPad *) c1->cd(idx++);
484  c1->Update();
485 
486  p->DrawFrame(11, 0, 50, 1.2)
487  ->SetTitle(";Transverse Momentum [GeV/#it{c}];#it{R}_{AA}");
488 
489  TGraph *g20 = pQCDModel_HuangKangVitev(2.0);
490  TGraph *g22 = pQCDModel_HuangKangVitev(2.2);
491 
492  g20->SetLineColor(kGreen + 3);
493  g22->SetLineColor(kRed + 3);
494  g20->SetLineWidth(3);
495  g22->SetLineWidth(3);
496 
497  g20->Draw("l");
498  g22->Draw("l");
499 
500  ge_RAA->DrawClone("pe");
501  leg->DrawClone();
502  leg2->DrawClone();
503 
504  TLegend *leg1 = new TLegend(.2, .20, .85, .37);
505  leg1->SetHeader("#splitline{pQCD, Phys.Lett. B726 (2013) 251-256}{#sqrt{s_{NN}}=200 GeV, #it{b}-jet R=0.4}");
506  leg1->AddEntry("", "", "");
507  leg1->AddEntry(g20, "g^{med} = 2.0", "l");
508  leg1->AddEntry(g22, "g^{med} = 2.2", "l");
509  leg1->Draw();
510 
511  SaveCanvas(c1, infile + "_" + TString(c1->GetName()), kTRUE);
512 }
513 
514 void CrossSection2v2(const TString infile, const bool use_AA_jet_trigger = true,const double ep_resolution = 0.7, const double dy = .7 * 2)
515 {
516  TFile *f = TFile::Open(infile + "Draw_HFJetTruth_DrawCrossSection.root");
517  assert(f);
518 
519  TH1F *hall = (TH1F *) f->GetObjectChecked("hall", "TH1F");
520  assert(hall);
521  TH1F *h_b = (TH1F *) f->GetObjectChecked("h_b", "TH1F");
522  assert(h_b);
523 
524  const double b_jet_RAA = 0.6;
525 
526  const double pp_eff = 0.6;
527  const double pp_purity = 0.4;
528  const double AuAu_eff = 0.4;
529  const double AuAu_purity = 0.4;
530 
532  // 5-year lumi in [sPH-TRG-000]
534 
535  const double pp_lumi = 200; // pb^-1 [sPH-TRG-000], rounded up from 197 pb^-1
536  const double pp_inelastic_crosssec = 42e-3 / 1e-12; // 42 mb in pb [sPH-TRG-000]
537 
538  const double AuAu_MB_Evt = use_AA_jet_trigger ? 550e9 : 240e9; // [sPH-TRG-000], depending on whether jet trigger applied in AA collisions
539  const double pAu_MB_Evt = 600e9; // [sPH-TRG-000]
540 
541  const double AuAu_Ncoll_C0_10 = 960.2; // [DOI:?10.1103/PhysRevC.87.034911?]
542  const double AuAu_Ncoll_C0_20 = 770.6; // [DOI:?10.1103/PhysRevC.91.064904?]
543  const double AuAu_Ncoll_C10_20 = 603; // [sPH-HF-2017-001-v1]
544  const double AuAu_Ncoll_C20_40 = 296; // [sPH-HF-2017-001-v1]
545  const double AuAu_Ncoll_C10_40 = (AuAu_Ncoll_C10_20*10 + AuAu_Ncoll_C20_40*20)/30;
546  const double AuAu_Ncoll_C40_60 = 94; // [sPH-HF-2017-001-v1]
547  const double AuAu_Ncoll_C60_92 = 15; // [sPH-HF-2017-001-v1]
548  const double AuAu_Ncoll_C0_100 = 250; // pb^-1 [sPH-TRG-000]
549  const double pAu_Ncoll_C0_100 = 4.7; // pb^-1 [sPH-TRG-000]
550 
551  const double AuAu_eq_lumi_C0_10 = AuAu_MB_Evt * 0.1 * AuAu_Ncoll_C0_10 / pp_inelastic_crosssec; //
552  const double AuAu_eq_lumi_C0_20 = AuAu_MB_Evt * 0.2 * AuAu_Ncoll_C0_20 / pp_inelastic_crosssec; //
553  const double AuAu_eq_lumi_C0_100 = AuAu_MB_Evt * 1 * AuAu_Ncoll_C0_100 / pp_inelastic_crosssec; //
554 
555  const double AuAu_eq_lumi_C10_20 = AuAu_MB_Evt * .1 * AuAu_Ncoll_C10_20 / pp_inelastic_crosssec; //
556  const double AuAu_eq_lumi_C20_40 = AuAu_MB_Evt * .2 * AuAu_Ncoll_C20_40 / pp_inelastic_crosssec; //
557  const double AuAu_eq_lumi_C10_40 = AuAu_MB_Evt * .3 * AuAu_Ncoll_C10_40 / pp_inelastic_crosssec; //
558  const double AuAu_eq_lumi_C40_60 = AuAu_MB_Evt * .2 * AuAu_Ncoll_C40_60 / pp_inelastic_crosssec; //
559  const double AuAu_eq_lumi_C60_92 = AuAu_MB_Evt * (.92 - .6) * AuAu_Ncoll_C60_92 / pp_inelastic_crosssec; //
560 
561  const double pAu_eq_lumi_C0_100 = pAu_MB_Evt * 1 * pAu_Ncoll_C0_100 / pp_inelastic_crosssec; //
562 
563  ;
564 
565 
566  const TString s_suffix(use_AA_jet_trigger ? "_AAJetTriggered" : "");
567  s_suffix += Form("_EPR%.1f",ep_resolution);
568  s_suffix += Form("_deta%.2f",dy/2);
569 
570  cout << "CrossSection2v2 integrated luminosity assumptions in pb^-1: " << endl;
571  cout << "\t"
572  << "pp_lumi = " << pp_lumi << endl;
573  cout << "\t"
574  << "AuAu_eq_lumi_C0_10 = " << AuAu_eq_lumi_C0_10 << endl;
575  cout << "\t"
576  << "AuAu_eq_lumi_C0_20 = " << AuAu_eq_lumi_C0_20 << endl;
577  cout << "\t"
578  << "AuAu_eq_lumi_C0_100 = " << AuAu_eq_lumi_C0_100 << endl;
579  cout << "\t"
580  << "AuAu_eq_lumi_C10_20 = " << AuAu_eq_lumi_C10_20 << endl;
581  cout << "\t"
582  << "AuAu_eq_lumi_C20_40 = " << AuAu_eq_lumi_C20_40 << endl;
583  cout << "\t"
584  << "AuAu_eq_lumi_C10_40 = " << AuAu_eq_lumi_C10_40 << endl;
585  cout << "\t"
586  << "AuAu_eq_lumi_C40_60 = " << AuAu_eq_lumi_C40_60 << endl;
587  cout << "\t"
588  << "AuAu_eq_lumi_C60_92 = " << AuAu_eq_lumi_C60_92 << endl;
589  cout << "\t"
590  << "pAu_eq_lumi_C0_100 = " << pAu_eq_lumi_C0_100 << endl;
591 
592  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);
593  TGraph *g_AA_C0_20 = CrossSection2v2Uncert(h_b, b_jet_RAA, dy, AuAu_eq_lumi_C0_20 * AuAu_eff * AuAu_purity, ep_resolution, 0);
594  TGraph *g_AA_C10_20 = CrossSection2v2Uncert(h_b, b_jet_RAA, dy, AuAu_eq_lumi_C10_20 * AuAu_eff * AuAu_purity, ep_resolution, 0);
595  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);
596  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);
597  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);
598  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);
599  //
600  g_AA_C0_10->SetLineColor(kBlue+3);
601  g_AA_C10_20->SetLineColor(kAzure+3);
602  g_AA_C20_40->SetLineColor(kTeal+3);
603  g_AA_C10_40->SetLineColor(kTeal+3);
604  g_AA_C40_60->SetLineColor(kSpring+3);
605 
606  g_AA_C0_10->SetMarkerColor(kBlue+3);
607  g_AA_C10_20->SetMarkerColor(kAzure+3);
608  g_AA_C20_40->SetMarkerColor(kTeal+3);
609  g_AA_C10_40->SetMarkerColor(kTeal+3);
610  g_AA_C40_60->SetMarkerColor(kSpring+3);
611 
612  g_AA_C0_10->SetMarkerStyle(kFullCircle);
613  g_AA_C10_20->SetMarkerStyle(kFullSquare);
614  g_AA_C20_40->SetMarkerStyle(kFullDiamond);
615  g_AA_C10_40->SetMarkerStyle(kFullDiamond);
616  g_AA_C40_60->SetMarkerStyle(kFullCross);
617 
618 
619  //
620  TCanvas *c1 = new TCanvas("Draw_HFJetTruth_CrossSection2v2" + s_suffix, "Draw_HFJetTruth_CrossSection2v2" + s_suffix, 700, 600);
621  c1->Divide(1, 1);
622  int idx = 1;
623  TPad *p;
624 
625  p = (TPad *) c1->cd(idx++);
626  c1->Update();
627 
628  p->DrawFrame(15, -.1, 40, .3)
629  ->SetTitle(";Transverse Momentum [GeV/#it{c}];v_{2}");
630  //
631  g_AA_C0_10->Draw("pe");
632 // g_AA_C10_20->Draw("pe");
633 // g_AA_C20_40->Draw("pe");
634  g_AA_C10_40->Draw("pe");
635 // g_AA_C40_60->Draw("pe");
636 
637  // g_AA_C20_40->Draw("same");
638  //
639  // ge_RAA->SetLineWidth(3);
640  // ge_RAA->SetMarkerStyle(kFullCircle);
641  // ge_RAA->SetMarkerSize(2);
642  //
643  // ge_RAA->Draw("pe");
644  // ge_RAA->Print();
645  //
646  TLegend *leg = new TLegend(.0, .78, .85, .93);
647  leg->SetFillStyle(0);
648  leg->AddEntry("", "#it{#bf{sPHENIX }} Simulation", "");
649  leg->AddEntry("", Form("PYTHIA-8 #it{b}-jet, Anti-k_{T} R=0.4, |#eta|<%.2f, CTEQ6L", dy / 2), "");
650  leg->AddEntry("", Form("Au+Au: %.0fB col., %.0f%% Eff., %.0f%% Pur.", '%', AuAu_MB_Evt / 1e9, AuAu_eff * 100, AuAu_purity * 100), "");
651  leg->Draw();
652  //
653  TLegend *leg2 = new TLegend(.19, .55, 1, .78);
654  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));
655  leg2->AddEntry(g_AA_C0_10, "Au+Au 0-10%C", "pl");
656 // leg2->AddEntry(g_AA_C10_20, "Au+Au 10-20%C", "pl");
657 // leg2->AddEntry(g_AA_C20_40, "Au+Au 20-40%C", "pl");
658  leg2->AddEntry(g_AA_C10_40, "Au+Au 10-40%C", "pl");
659 // leg2->AddEntry(g_AA_C40_60, "Au+Au 40-60%C", "pl");
660  leg2->SetFillStyle(0);
661  leg2->Draw();
662 
663  SaveCanvas(c1, infile + "_" + TString(c1->GetName()), kTRUE);
664 }
665 
666 TGraphErrors *GetRAA(TH1 *h_pp, TH1 *h_AA)
667 {
668  int n_bin = 0;
669 
670  double xs[1000] = {0};
671  double ys[1000] = {0};
672  double eys[1000] = {0};
673 
674  assert(h_pp);
675  assert(h_AA);
676  assert(h_pp->GetNbinsX() == h_AA->GetNbinsX());
677 
678  for (int i = 1; i <= h_pp->GetNbinsX(); ++i)
679  {
680  if (h_pp->GetBinError(i) > 0 && h_pp->GetBinError(i) < 1 && h_AA->GetBinError(i) > 0 && h_AA->GetBinError(i) < 1)
681  {
682  xs[n_bin] = h_pp->GetXaxis()->GetBinCenter(i);
683 
684  ys[n_bin] = h_AA->GetBinContent(i) / h_pp->GetBinContent(i);
685 
686  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));
687 
688  n_bin += 1;
689  }
690  }
691 
692  TGraphErrors *ge = new TGraphErrors(n_bin, xs, ys, NULL, eys);
693  ge->SetName(TString("RAA_") + h_AA->GetName());
694 
695  return ge;
696 }
697 
698 TH1 *CrossSection2RelUncert(const TH1F *h_cross,
699  const double suppression,
700  const double deta,
701  const double pp_quiv_int_lum)
702 {
703  assert(h_cross);
704  TH1 *
705  h_ratio = (TH1 *)
706  h_cross->Clone(TString(h_cross->GetName()) + Form("_copy%d", rand()));
707 
708  //convert to count per bin
709  h_ratio->Scale(deta * h_ratio->GetXaxis()->GetBinWidth(0) * pp_quiv_int_lum * suppression);
710  h_ratio->Rebin(5);
711 
712  for (int i = 1; i <= h_ratio->GetNbinsX(); ++i)
713  {
714  const double yield = h_ratio->GetBinContent(i);
715 
716  if (yield > 0)
717  {
718  h_ratio->SetBinContent(i, suppression);
719 
720  h_ratio->SetBinError(i, suppression / sqrt(yield));
721  }
722  else
723  {
724  h_ratio->SetBinContent(i, 0);
725 
726  h_ratio->SetBinError(i, 0);
727  }
728  }
729 
730  h_ratio->GetYaxis()->SetTitle("Relative Cross Section and Uncertainty");
731 
732  return h_ratio;
733 }
734 
735 TGraph *CrossSection2v2Uncert(const TH1F *h_cross,
736  const double suppression,
737  const double deta,
738  const double pp_quiv_int_lum,
739  const double ep_resolution = 1,
740  const double pt_shift = 0
741 )
742 {
743  assert(h_cross);
744  TH1 *
745  h_ratio = (TH1 *)
746  h_cross->Clone(TString(h_cross->GetName()) + Form("_copyv2%d", rand()));
747 
748  //convert to count per bin
749  h_ratio->Scale(deta * h_ratio->GetXaxis()->GetBinWidth(0) * pp_quiv_int_lum * suppression);
750  h_ratio->Rebin(5);
751 
752  vector<double> pts;
753  vector<double> v2s;
754  vector<double> v2es;
755 
756  for (int i = 1; i <= h_ratio->GetNbinsX(); ++i)
757  {
758  const double yield = h_ratio->GetBinContent(i);
759 
760  if (yield > 100)
761  {
762  h_ratio->SetBinContent(i, 0);
763 
764  h_ratio->SetBinError(i, 1. / sqrt(2 * yield) / ep_resolution); // err(v2) = 1/ (sqrt(2) *Significance * Resolution)
765 
766  pts.push_back(h_ratio->GetBinCenter(i) + pt_shift);
767  v2s.push_back(h_ratio->GetBinContent(i));
768  v2es.push_back(h_ratio->GetBinError(i));
769  }
770  else
771  {
772  h_ratio->SetBinContent(i, 0);
773 
774  h_ratio->SetBinError(i, 0);
775  }
776  }
777 
778  h_ratio->GetYaxis()->SetTitle("v2 and uncertainty");
779 
780  TGraph *gr = new TGraphErrors(pts.size(), &pts[0], &v2s[0], 0, &v2es[0]);
781  gr->SetName(TString("ge_") + h_ratio->GetName());
782 
783  gr->SetLineWidth(3);
784  gr->SetMarkerStyle(kFullCircle);
785  gr->SetMarkerSize(2);
786 
787  return gr;
788 }
789 
790 void Convert2CrossSection(TH1 *h, const double int_lumi, const double dy)
791 {
792  h->Sumw2();
793 
794  for (int i = 1; i <= h->GetXaxis()->GetNbins(); ++i)
795  {
796  const double pT1 = h->GetXaxis()->GetBinLowEdge(i);
797  const double pT2 = h->GetXaxis()->GetBinUpEdge(i);
798 
799  // const double dpT2 = pT2*pT2 - pT1*pT1;
800  const double dpT = pT2 - pT1;
801 
802  // const double count2sigma = 1./dpT2/dy/int_lumi;
803  const double count2sigma = 1. / dpT / dy / int_lumi;
804 
805  h->SetBinContent(i, h->GetBinContent(i) * count2sigma);
806  h->SetBinError(i, h->GetBinError(i) * count2sigma);
807  }
808 
809  // h->GetYaxis()->SetTitle("#sigma/(dp_{T}^{2} dy) (pb/(GeV/c)^{2})");
810 
811  h->SetMarkerStyle(kFullCircle);
812 }
813 
814 TGraphAsymmErrors *
815 GetFONLL_B()
816 {
817  //# Job started on: Wed Aug 24 04:28:34 CEST 2016 .
818  //# FONLL heavy quark hadroproduction cross section, calculated on Wed Aug 24 04:28:42 CEST 2016
819  //# FONLL version and perturbative order: ## FONLL v1.3.2 fonll [ds/dpt^2dy (pb/GeV^2)]
820  //# quark = bottom
821  //# final state = quark
822  //# ebeam1 = 100, ebeam2 = 100
823  //# PDF set = CTEQ6.6
824  //# ptmin = 10
825  //# ptmax = 30
826  //# etamin = -.6
827  //# etamax = .6
828  //# Uncertainties from scales, masses combined quadratically
829  //# cross section is ds/dpt (pb/GeV)
830  //# pt central min max min_sc max_sc min_mass max_mass
831  // 10.0000 4.1645e+03 2.9134e+03 6.0916e+03 2.9634e+03 6.0553e+03 3.8146e+03 4.5365e+03
832  // 11.0526 2.4830e+03 1.7438e+03 3.6234e+03 1.7667e+03 3.6072e+03 2.3005e+03 2.6743e+03
833  // 12.1053 1.5046e+03 1.0602e+03 2.1900e+03 1.0710e+03 2.1826e+03 1.4071e+03 1.6056e+03
834  // 13.1579 9.2707e+02 6.5506e+02 1.3454e+03 6.6029e+02 1.3418e+03 8.7394e+02 9.8137e+02
835  // 14.2105 5.8027e+02 4.1098e+02 8.3990e+02 4.1358e+02 8.3817e+02 5.5067e+02 6.1015e+02
836  // 15.2632 3.6865e+02 2.6159e+02 5.3224e+02 2.6290e+02 5.3139e+02 3.5193e+02 3.8532e+02
837  // 16.3158 2.3744e+02 1.6874e+02 3.4198e+02 1.6941e+02 3.4156e+02 2.2785e+02 2.4686e+02
838  // 17.3684 1.5484e+02 1.1016e+02 2.2257e+02 1.1051e+02 2.2236e+02 1.4930e+02 1.6021e+02
839  // 18.4211 1.0215e+02 7.2725e+01 1.4657e+02 7.2904e+01 1.4646e+02 9.8911e+01 1.0522e+02
840  // 19.4737 6.8101e+01 4.8502e+01 9.7587e+01 4.8594e+01 9.7535e+01 6.6205e+01 6.9866e+01
841  // 20.5263 4.5818e+01 3.2631e+01 6.5588e+01 3.2678e+01 6.5562e+01 4.4707e+01 4.6827e+01
842  // 21.5789 3.1111e+01 2.2148e+01 4.4506e+01 2.2171e+01 4.4494e+01 3.0461e+01 3.1684e+01
843  // 22.6316 2.1287e+01 1.5142e+01 3.0444e+01 1.5153e+01 3.0439e+01 2.0909e+01 2.1607e+01
844  // 23.6842 1.4663e+01 1.0419e+01 2.0974e+01 1.0424e+01 2.0972e+01 1.4446e+01 1.4838e+01
845  // 24.7368 1.0175e+01 7.2188e+00 1.4562e+01 7.2214e+00 1.4562e+01 1.0053e+01 1.0266e+01
846  // 25.7895 7.1026e+00 5.0296e+00 1.0176e+01 5.0307e+00 1.0175e+01 7.0360e+00 7.1470e+00
847  // 26.8421 4.9796e+00 3.5183e+00 7.1439e+00 3.5187e+00 7.1438e+00 4.9454e+00 4.9976e+00
848  // 27.8947 3.5103e+00 2.4738e+00 5.0451e+00 2.4739e+00 5.0451e+00 3.4945e+00 3.5143e+00
849  // 28.9474 2.4892e+00 1.7492e+00 3.5856e+00 1.7493e+00 3.5856e+00 2.4835e+00 2.4892e+00
850  // 30.0000 1.7693e+00 1.2394e+00 2.5556e+00 1.2394e+00 2.5556e+00 1.7633e+00 1.7693e+00
851 
852  const double deta = 0.6 * 2;
853 
854  const double pts[] =
855  {
856 
857  10.0000, 11.0526, 12.1053, 13.1579, 14.2105, 15.2632, 16.3158, 17.3684,
858  18.4211, 19.4737, 20.5263, 21.5789, 22.6316, 23.6842, 24.7368, 25.7895,
859  26.8421, 27.8947, 28.9474, 30.0000
860 
861  };
862 
863  const double centrals[] =
864  {4.1645e+03, 2.4830e+03, 1.5046e+03, 9.2707e+02, 5.8027e+02, 3.6865e+02,
865  2.3744e+02, 1.5484e+02, 1.0215e+02, 6.8101e+01, 4.5818e+01, 3.1111e+01,
866  2.1287e+01, 1.4663e+01, 1.0175e+01, 7.1026e+00, 4.9796e+00, 3.5103e+00,
867  2.4892e+00, 1.7693e+00};
868  const double min[] =
869  {4.1645e+03 - 2.9134e+03, 2.4830e+03 - 1.7438e+03, 1.5046e+03 - 1.0602e+03,
870  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,
871  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,
872  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,
873  2.4892e+00 - 1.7492e+00, 1.7693e+00 - 1.2394e+00};
874 
875  const double max[] =
876  {6.0916e+03 - 4.1645e+03, 3.6234e+03 - 2.4830e+03, 2.1900e+03 - 1.5046e+03,
877  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,
878  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,
879  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,
880  3.5856e+00 - 2.4892e+00, 2.5556e+00 - 1.7693e+00};
881 
882  TGraphAsymmErrors *gr = new TGraphAsymmErrors(20, pts, centrals, 0, 0, min,
883  max);
884 
885  for (int i = 0; i < gr->GetN(); ++i)
886  {
887  (gr->GetY())[i] /= deta;
888  (gr->GetEYhigh())[i] /= deta;
889  (gr->GetEYlow())[i] /= deta;
890  }
891 
892  return gr;
893 }
894 
895 TGraphAsymmErrors *
896 GetFONLL_C()
897 {
898  //# Job started on: Thu Aug 25 05:08:57 CEST 2016 .
899  //# FONLL heavy quark hadroproduction cross section, calculated on Thu Aug 25 05:09:06 CEST 2016
900  //# FONLL version and perturbative order: ## FONLL v1.3.2 fonll [ds/dpt^2dy (pb/GeV^2)]
901  //# quark = charm
902  //# final state = quark
903  //# ebeam1 = 100, ebeam2 = 100
904  //# PDF set = CTEQ6.6
905  //# ptmin = 10
906  //# ptmax = 30
907  //# etamin = -.6
908  //# etamax = .6
909  //# Uncertainties from scales, masses combined quadratically
910  //# cross section is ds/dpt (pb/GeV)
911  //# pt central min max min_sc max_sc min_mass max_mass
912  // 10.0000 9.7368e+03 6.5051e+03 1.5528e+04 6.5641e+03 1.5522e+04 9.1225e+03 1.0006e+04
913  // 11.0526 4.9069e+03 3.2787e+03 7.7308e+03 3.3256e+03 7.7226e+03 4.5190e+03 5.1221e+03
914  // 12.1053 2.5743e+03 1.7187e+03 4.0147e+03 1.7525e+03 4.0069e+03 2.3363e+03 2.7237e+03
915  // 13.1579 1.3990e+03 9.3274e+02 2.1641e+03 9.5597e+02 2.1578e+03 1.2537e+03 1.4976e+03
916  // 14.2105 7.8344e+02 5.2138e+02 1.2037e+03 5.3701e+02 1.1989e+03 6.9430e+02 8.4730e+02
917  // 15.2632 4.5087e+02 2.9938e+02 6.8902e+02 3.0978e+02 6.8543e+02 3.9570e+02 4.9204e+02
918  // 16.3158 2.6550e+02 1.7584e+02 4.0398e+02 1.8275e+02 4.0140e+02 2.3098e+02 2.9210e+02
919  // 17.3684 1.5978e+02 1.0551e+02 2.4232e+02 1.1011e+02 2.4049e+02 1.3793e+02 1.7705e+02
920  // 18.4211 9.7903e+01 6.4459e+01 1.4811e+02 6.7524e+01 1.4682e+02 8.3915e+01 1.0919e+02
921  // 19.4737 6.0993e+01 4.0030e+01 9.2103e+01 4.2085e+01 9.1200e+01 5.1941e+01 6.8435e+01
922  // 20.5263 3.8605e+01 2.5248e+01 5.8239e+01 2.6631e+01 5.7606e+01 3.2685e+01 4.3550e+01
923  // 21.5789 2.4727e+01 1.6108e+01 3.7286e+01 1.7044e+01 3.6842e+01 2.0821e+01 2.8038e+01
924  // 22.6316 1.6036e+01 1.0402e+01 2.4185e+01 1.1040e+01 2.3872e+01 1.3433e+01 1.8271e+01
925  // 23.6842 1.0532e+01 6.8006e+00 1.5898e+01 7.2389e+00 1.5678e+01 8.7774e+00 1.2053e+01
926  // 24.7368 6.9729e+00 4.4803e+00 1.0540e+01 4.7838e+00 1.0385e+01 5.7809e+00 8.0136e+00
927  // 25.7895 4.6560e+00 2.9760e+00 7.0510e+00 3.1876e+00 6.9411e+00 3.8397e+00 5.3729e+00
928  // 26.8421 3.1409e+00 1.9970e+00 4.7681e+00 2.1454e+00 4.6902e+00 2.5775e+00 3.6384e+00
929  // 27.8947 2.1346e+00 1.3500e+00 3.2500e+00 1.4543e+00 3.1946e+00 1.7438e+00 2.4816e+00
930  // 28.9474 1.4557e+00 9.1535e-01 2.2239e+00 9.8868e-01 2.1845e+00 1.1839e+00 1.6985e+00
931  // 30.0000 9.9826e-01 6.2400e-01 1.5310e+00 6.7571e-01 1.5029e+00 8.0844e-01 1.1690e+00
932 
933  const double deta = 0.6 * 2;
934 
935  const double pts[] =
936  {
937  // pt
938  10.0000, 11.0526, 12.1053, 13.1579, 14.2105, 15.2632, 16.3158, 17.3684,
939  18.4211, 19.4737, 20.5263, 21.5789, 22.6316, 23.6842, 24.7368, 25.7895,
940  26.8421, 27.8947, 28.9474, 30.0000
941 
942  };
943 
944  const double centrals[] =
945  {
946 
947  // central
948  9.7368e+03, 4.9069e+03, 2.5743e+03, 1.3990e+03, 7.8344e+02, 4.5087e+02,
949  2.6550e+02, 1.5978e+02, 9.7903e+01, 6.0993e+01, 3.8605e+01, 2.4727e+01,
950  1.6036e+01, 1.0532e+01, 6.9729e+00, 4.6560e+00, 3.1409e+00, 2.1346e+00,
951  1.4557e+00, 9.9826e-01
952 
953  };
954  const double min[] =
955  {
956 
957  // central min
958  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,
959  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,
960  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,
961  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
962 
963  };
964 
965  const double max[] =
966  {
967 
968  // max central
969  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,
970  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,
971  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,
972  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
973 
974  };
975 
976  TGraphAsymmErrors *gr = new TGraphAsymmErrors(20, pts, centrals, 0, 0, min,
977  max);
978 
979  for (int i = 0; i < gr->GetN(); ++i)
980  {
981  (gr->GetY())[i] /= deta;
982  (gr->GetEYhigh())[i] /= deta;
983  (gr->GetEYlow())[i] /= deta;
984  }
985 
986  return gr;
987 }
988 
989 TGraph *pQCDModel_HuangKangVitev(const double g)
990 {
991  // arXiv:1306.0909v2 [hep-ph]
992  // Preliminary for sPHENIX energy
993 
994  // R = 0.4;
995  if (g == 2.0)
996  {
997  // g=2.0
998  // 10.130739 0.8087402
999  // 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 }