TruthGeneration/macros/Draw_HFJetTruth_InvMass.C

0001 // $Id: $
0002
0003 /*!
0004  * \file Draw_HFJetTruth_InvMass.C
0005  * \brief
0006  * \author Jin Huang <jhuang@bnl.gov>
0007  * \version $Revision:   $
0008  * \date $Date: $
0009  */
0010
0011 #include <TFile.h>
0012 #include <TGraphAsymmErrors.h>
0013 #include <TGraphErrors.h>
0014 #include <TLatex.h>
0015 #include <TLine.h>
0016 #include <TString.h>
0017 #include <TTree.h>
0018 #include <cassert>
0019 #include <cmath>
0020 #include "SaveCanvas.C"
0021 #include "sPhenixStyle.C"
0022
0023 TFile *_file0 = NULL;
0024 TTree *T = NULL;
0025
0026 void Draw_HFJetTruth_InvMass(const TString infile =
0027                                  //                         "/sphenix/user/jinhuang/HF-jet/event_gen/200pp_pythia8_CTEQ6L_20GeV/200pp_pythia8_CTEQ6L_20GeV_ALL.cfg_eneg_DSTReader.root",
0028                              //                     double int_lumi = 210715 / 5.533e-05 / 1e9, const double dy = 0.6 * 2)
0029
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
0033 //"/sphenix/user/jinhuang/HF-jet/event_gen/200pp_pythia8_CTEQ6L/200pp_pythia8_CTEQ6L_111ALL.cfg_eneg_DSTReader.root",
0034 //double int_lumi = 789908/4.631e-03 / 1e9, const double dy = 0.6*2)
0035
0036 //Draw_HFJetTruth_InvMass(const TString infile =
0037 //    "../macros3/G4sPHENIXCells.root_DSTReader.root",
0038 //    double int_lumi = 789908/4.631e-03 / 1e9, const double dy = 0.6*2)
0039 //Draw_HFJetTruth_InvMass(const TString infile =
0040 //    "../macros2/G4sPHENIXCells.root_DSTReader.root",
0041 //    const double int_lumi = 1842152 / 5.801e-03 / 1e9,
0042 //    const double dy = 0.6 * 2)
0043 {
0044   SetsPhenixStyle();
0045   gStyle->SetOptStat(0);
0046   gStyle->SetOptFit(1111);
0047   TVirtualFitter::SetDefaultFitter("Minuit2");
0048
0049   gSystem->Load("libg4eval.so");
0050   gSystem->Load("libg4dst.so");
0051
0052   if (!_file0)
0053   {
0054     TString chian_str = infile;
0055     chian_str.ReplaceAll("ALL", "*");
0056
0057     TChain *t = new TChain("T");
0058     const int n = t->Add(chian_str);
0059
0060     int_lumi *= n;
0061     //      cout << "Loaded " << n << " root files with " << chian_str << endl;
0062     cout << "int_lumi = " << int_lumi << " pb^-1 from " << n << " root files with " << chian_str << endl;
0063     assert(n > 0);
0064
0065     T = t;
0066
0067     _file0 = new TFile;
0068     _file0->SetName(infile);
0069   }
0070
0071   // step1: get the cross section
0072   //    DrawCrossSection(int_lumi, dy);
0073
0074   // step2: ploting
0075
0076   // Draw_HFJetTruth_InvMass_DrawCrossSection_PR(infile);
0077   //  CrossSection2RAA_Proposal(infile);
0078   //      CrossSection2RAA(infile);
0079   const double acceptance1 = 2. * (1.1 - .4);
0080   // CrossSection2RAA(infile, false, acceptance1);
0081   CrossSection2RAA(infile, false, acceptance1, true);
0082
0083   //  const double acceptance2 = 2.* (0.85 - .4);
0084   //  CrossSection2RAA(infile, false, acceptance2);
0085
0086   //  CrossSection2v2(infile, false, .7, acceptance);
0087   //  CrossSection2v2(infile, false, .4, acceptance);
0088 }
0089
0090 void DrawCrossSection(double int_lumi, const double dy)
0091 {
0092   TCut basicDijetEventQA("(AntiKt_Truth_r04[n_AntiKt_Truth_r04-1].get_pt())>15 && (AntiKt_Truth_r04[n_AntiKt_Truth_r04-2].get_pt())>10 && abs(AntiKt_Truth_r04[n_AntiKt_Truth_r04-1].get_eta())<.6 && abs(AntiKt_Truth_r04[n_AntiKt_Truth_r04-2].get_eta())<.6 && cos(AntiKt_Truth_r04[n_AntiKt_Truth_r04-1].get_phi() - AntiKt_Truth_r04[n_AntiKt_Truth_r04-2].get_phi())<-.5");
0093
0094   cout << "Build event selection of " << (const char *) basicDijetEventQA << endl;
0095
0096   T->Draw(">>EventList", basicDijetEventQA);
0097   TEventList *elist = gDirectory->GetObjectChecked("EventList", "TEventList");
0098   assert(elist);
0099   cout << elist->GetN() << " / " << T->GetEntriesFast() << " events selected" << endl;
0100   T->SetEventList(elist);
0101
0102   T->SetAlias("DiJetInvMass", "sqrt( (AntiKt_Truth_r04[n_AntiKt_Truth_r04-1].get_e() + AntiKt_Truth_r04[n_AntiKt_Truth_r04-2].get_e())**2 - (AntiKt_Truth_r04[n_AntiKt_Truth_r04-1].get_px() + AntiKt_Truth_r04[n_AntiKt_Truth_r04-2].get_px())**2 - (AntiKt_Truth_r04[n_AntiKt_Truth_r04-1].get_py() + AntiKt_Truth_r04[n_AntiKt_Truth_r04-2].get_py())**2 - (AntiKt_Truth_r04[n_AntiKt_Truth_r04-1].get_pz() + AntiKt_Truth_r04[n_AntiKt_Truth_r04-2].get_pz())**2  )");
0103
0104   TH1 *hall = new TH1F("hall", ";m_{12} [GeV]", 200, 0, 200);
0105   TH1 *h_b = new TH1F("h_b", ";m_{12} [GeV]", 200, 0, 200);
0106   TH1 *h_bq = new TH1F("h_bq", ";m_{12} [GeV]", 200, 0, 200);
0107   TH1 *h_bh = new TH1F("h_bh", ";m_{12} [GeV]", 200, 0, 200);
0108   TH1 *h_bh5 = new TH1F("h_bh5", ";m_{12} [GeV]", 200, 0, 200);
0109   TH1 *h_ch5 = new TH1F("h_ch5", ";m_{12} [GeV]", 200, 0, 200);
0110   TH1 *h_c = new TH1F("h_c", ";m_{12} [GeV]", 200, 0, 200);
0111
0112   cout << "DiJetInvMass>>hall" << endl;
0113   T->Draw("DiJetInvMass>>hall",
0114           "",
0115           "goff");
0116
0117   cout << "DiJetInvMass>>h_b" << endl;
0118   T->Draw("DiJetInvMass>>h_b",
0119           "abs(AntiKt_Truth_r04[n_AntiKt_Truth_r04-1].get_property(1000))==5 && abs(AntiKt_Truth_r04[n_AntiKt_Truth_r04-2].get_property(1000))==5",
0120           "goff");
0121   //  T->Draw(
0122   //      "AntiKt_Truth_r04.get_pt()* AntiKt_Truth_r04.get_property(1001) >>h_bq",
0123   //      "AntiKt_Truth_r04.get_pt()>15 && abs(AntiKt_Truth_r04.get_eta())<0.6 && abs(AntiKt_Truth_r04.get_property(1000))==5",
0124   //      "goff");
0125
0126   //    T->Draw("AntiKt_Truth_r04.get_pt()>>h_bh",
0127   //        "AntiKt_Truth_r04.get_pt()>15 && abs(AntiKt_Truth_r04.get_eta())<0.6 && AntiKt_Truth_r04.get_property(1010)==5",
0128   //        "goff");
0129
0130   //  T->Draw("AntiKt_Truth_r04.get_pt()>>h_bh5",
0131   //          "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",
0132   //          "goff");
0133   //
0134   //  T->Draw("AntiKt_Truth_r04.get_pt()>>h_c",
0135   //          "AntiKt_Truth_r04.get_pt()>15 && abs(AntiKt_Truth_r04.get_eta())<0.6 && abs(AntiKt_Truth_r04.get_property(1000))==4",
0136   //          "goff");
0137   //  T->Draw("AntiKt_Truth_r04.get_pt()>>h_ch5",
0138   //          "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",
0139   //          "goff");
0140
0141   Convert2CrossSection(hall, int_lumi, dy);
0142   Convert2CrossSection(h_b, int_lumi, dy);
0143   //  Convert2CrossSection(h_bq, int_lumi, dy);
0144   //  Convert2CrossSection(h_bh5, int_lumi, dy);
0145   //  Convert2CrossSection(h_ch5, int_lumi, dy);
0146   //  Convert2CrossSection(h_c, int_lumi, dy);
0147
0148   h_b->SetLineColor(kBlue);
0149   h_b->SetMarkerColor(kBlue);
0150
0151   h_bq->SetLineColor(kBlue + 3);
0152   h_bq->SetMarkerColor(kBlue + 3);
0153
0154   h_bh5->SetLineColor(kBlue + 3);
0155   h_bh5->SetMarkerColor(kBlue + 3);
0156   h_bh5->SetMarkerStyle(kStar);
0157
0158   h_c->SetLineColor(kRed);
0159   h_c->SetMarkerColor(kRed);
0160
0161   h_ch5->SetLineColor(kRed + 3);
0162   h_ch5->SetMarkerColor(kRed + 3);
0163   h_ch5->SetMarkerStyle(kStar);
0164
0165   //  TGraphAsymmErrors *gr_fonll_b = GetFONLL_B();
0166   //  gr_fonll_b->SetFillColor(kBlue - 7);
0167   //  gr_fonll_b->SetFillStyle(3002);
0168   //  TGraphAsymmErrors *gr_fonll_c = GetFONLL_C();
0169   //  gr_fonll_c->SetFillColor(kRed - 7);
0170   //  gr_fonll_c->SetFillStyle(3003);
0171   //  TGraph *gr_phenix = GetPHENIX_jet();
0172   //  gr_phenix->SetLineColor(kGray + 2);
0173   //  gr_phenix->SetMarkerColor(kBlack);
0174   //  gr_phenix->SetMarkerStyle(kOpenCross);
0175
0176   TCanvas *c1 = new TCanvas("Draw_HFJetTruth_InvMass_DrawCrossSection", "Draw_HFJetTruth_InvMass_DrawCrossSection", 1000, 860);
0177   //  c1->Divide(2, 2);
0178   int idx = 1;
0179   TPad *p;
0180
0181   p = (TPad *) c1->cd(idx++);
0182   c1->Update();
0183   p->SetGridx(0);
0184   p->SetGridy(0);
0185   p->SetLogy();
0186
0187   hall->Draw();
0188
0189   h_b->Draw("same");
0190   //  h_bh5->Draw("same");
0191   //  h_bq->Draw("same");
0192   //  h_c->Draw("same");
0193   //  h_ch5->Draw("same");
0194
0195   //  hall->GetXaxis()->SetRangeUser(30, 160);
0196   hall->GetYaxis()->SetTitle(
0197       "d^{3}#sigma/(d#eta_{1}d#eta_{2}dm_{12}) [pb/(GeV)]");
0198
0199   TLegend *leg = new TLegend(0.45, 0.6, 0.95, 0.93);
0200   leg->SetFillColor(kWhite);
0201   leg->SetFillStyle(1001);
0202   leg->SetHeader("#splitline{#it{#bf{sPHENIX}} fast simulation, #it{p}+#it{p} #sqrt{s} = 200 GeV}{|#eta_{1,2}|<0.6, |#Delta#phi_{1,2}|>2#pi/3, p_{T,1}>15 GeV/c, p_{T,2}>10 GeV/c}");
0203   leg->AddEntry(hall, "Inclusive jet, PYTHIA8, Truth, anti-k_{t}, R=0.4",
0204                 "lpe");
0205   //  leg->AddEntry(h_c, "c-quark jet, Pythia8, Truth, anti-k_{t}, R=0.4", "lpe");
0206   leg->AddEntry(h_b, "b-quark jet, PYTHIA8, Truth, anti-k_{t}, R=0.4", "lpe");
0207   //  leg->AddEntry(h_bq, "Leading b-quark in jet, Pythia8, anti-k_{t}, R=0.4", "lpe");
0208   //  leg->AddEntry(h_bh5,
0209   //                "b-hadron jet, Pythia8, Truth, anti-k_{t}, R=0.4, p_{T, b-hadron}>5 GeV/c",
0210   //                "lpe");
0211   //  leg->AddEntry(gr_phenix, "PHENIX inclusive jet, PRL 116, 122301 (2016)", "ple");
0212   //  leg->AddEntry(gr_fonll_c, "c-quark, FONLL v1.3.2, CTEQ6.6, |y|<0.6", "f");
0213   //  leg->AddEntry(gr_fonll_b, "b-quark, FONLL v1.3.2, CTEQ6.6, |y|<0.6", "f");
0214   leg->Draw();
0215
0216   SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
0217 }
0218
0219 void Draw_HFJetTruth_InvMass_DrawCrossSection_PR(const TString infile)
0220 {
0221   TFile *f = TFile::Open(infile + "Draw_HFJetTruth_InvMass_DrawCrossSection.root");
0222   assert(f);
0223
0224   TH1F *hall = (TH1F *) f->GetObjectChecked("hall", "TH1F");
0225   assert(hall);
0226   TH1F *h_b = (TH1F *) f->GetObjectChecked("h_b", "TH1F");
0227   assert(h_b);
0228
0229   hall->SetMarkerColor(kBlack);
0230   hall->SetLineColor(kBlack);
0231   hall->SetMarkerStyle(kFullCircle);
0232
0233   h_b->SetMarkerColor(kBlue + 2);
0234   h_b->SetLineColor(kBlue + 2);
0235   h_b->SetMarkerStyle(kFullCircle);
0236
0237   TGraph *gr_star = GetSTAR2019();
0238   gr_star->SetLineColor(kGray + 2);
0239   gr_star->SetLineWidth(3);
0240   gr_star->SetMarkerColor(kGray + 2);
0241   gr_star->SetMarkerSize(2);
0242   gr_star->SetMarkerStyle(kFullSquare);
0243
0244   TCanvas *c1 = new TCanvas("Draw_HFJetTruth_InvMass_DrawCrossSection_PR", "Draw_HFJetTruth_InvMass_DrawCrossSection_PR", 1000, 860);
0245   //  c1->Divide(2, 2);
0246   int idx = 1;
0247   TPad *p;
0248
0249   p = (TPad *) c1->cd(idx++);
0250   c1->Update();
0251   p->SetLogy();
0252
0253   TH1 *hframe = p->DrawFrame(35, 1e-2, 140, 1e6);
0254   hframe->SetTitle(";m_{12} [GeV/c^{2}];d^{3}#sigma/(d#eta_{1}d#eta_{2}dm_{12}) [pb/(GeV/c^{2})]");
0255
0256
0257   gr_star->Draw("pe");
0258   hall->Draw("same");
0259   h_b->Draw("same");
0260
0261   TLegend *leg = new TLegend(0.2, 0.7, 0.95, 0.92);
0262   leg->SetFillColor(kWhite);
0263   leg->SetFillStyle(1001);
0264   //  leg->SetHeader("#splitline{#it{#bf{sPHENIX }} Simulation}{#it{p}+#it{p}, #sqrt{s} = 200 GeV, |#eta|<0.6}");
0265   leg->SetHeader("#splitline{#it{#bf{sPHENIX}} fast simulation, #it{p}+#it{p} #sqrt{s} = 200 GeV}{|#eta_{1,2}|<0.6, |#Delta#phi_{1,2}|>2#pi/3, p_{T,1}>15 GeV/c, p_{T,2}>10 GeV/c}");
0266   leg->AddEntry("", "",
0267                 "");
0268   leg->AddEntry(hall, "Inclusive jet, PYTHIA8 + CTEQ6L, anti-k_{T} R=0.4",
0269                 "lpe");
0270   leg->AddEntry(gr_star, "STAR, PRD95 (2017) no.7, 071103, R=0.6", "ple");
0271   leg->AddEntry(h_b, "#it{b}-quark jet, PYTHIA8 + CTEQ6L, anti-k_{T} R=0.4", "lpe");
0272
0273   leg->Draw();
0274
0275   SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
0276 }
0277
0278 void CrossSection2RAA(const TString infile, const bool use_AA_jet_trigger = true, const double dy = .7 * 2, const bool b3yr = false)
0279 {
0280   TFile *f = TFile::Open(infile + "Draw_HFJetTruth_InvMass_DrawCrossSection.root");
0281   assert(f);
0282
0283   TH1F *hall = (TH1F *) f->GetObjectChecked("hall", "TH1F");
0284   assert(hall);
0285   TH1F *h_b = (TH1F *) f->GetObjectChecked("h_b", "TH1F");
0286   assert(h_b);
0287
0288   TString s_suffix(use_AA_jet_trigger ? "_AAJetTriggered" : "");
0289   s_suffix += b3yr ? "_3yr" : "";
0290   s_suffix += Form("_deta%.2f", dy / 2);
0291
0292   const double b_jet_RAA = 0.4;
0293   const double target_RAA_ratio = 4;;
0294
0295   const double pp_eff = 0.6;
0296   const double pp_purity = 0.4;
0297   const double AuAu_eff = 0.4;
0298   const double AuAu_purity = 0.4;
0299
0300   ////////////////////////////
0301   // 5-year lumi in [sPH-TRG-000]
0302   ////////////////////////////
0303
0304   const double pp_lumi = b3yr ? 62 : 200;                          // pb^-1 [sPH-TRG-000], rounded up from 197 pb^-1
0305   const double pp_inelastic_crosssec = 42e-3 / 1e-12;  // 42 mb in pb [sPH-TRG-000]
0306
0307   const double AuAu_MB_Evt = use_AA_jet_trigger ? (b3yr ? 320e9 : 550e9) : (b3yr ? 142e9 : 240e9);  // [sPH-TRG-000], depending on whether jet trigger applied in AA collisions
0308   const double pAu_MB_Evt = 600e9;                                // [sPH-TRG-000]
0309
0310   const double AuAu_Ncoll_C0_10 = 960.2;  // [DOI:?10.1103/PhysRevC.87.034911?]
0311   const double AuAu_Ncoll_C0_20 = 770.6;  // [DOI:?10.1103/PhysRevC.91.064904?]
0312   const double AuAu_Ncoll_C0_100 = 250;   // pb^-1 [sPH-TRG-000]
0313   const double pAu_Ncoll_C0_100 = 4.7;    // pb^-1 [sPH-TRG-000]
0314
0315   const double AuAu_eq_lumi_C0_10 = AuAu_MB_Evt * 0.1 * AuAu_Ncoll_C0_10 / pp_inelastic_crosssec;  //
0316   const double AuAu_eq_lumi_C0_20 = AuAu_MB_Evt * 0.2 * AuAu_Ncoll_C0_20 / pp_inelastic_crosssec;  //
0317   const double AuAu_eq_lumi_C0_100 = AuAu_MB_Evt * 1 * AuAu_Ncoll_C0_100 / pp_inelastic_crosssec;  //
0318
0319   const double pAu_eq_lumi_C0_100 = pAu_MB_Evt * 1 * pAu_Ncoll_C0_100 / pp_inelastic_crosssec;  //
0320
0321   cout << "CrossSection2RAA integrated luminosity assumptions in pb^-1: " << endl;
0322   cout << "\t"
0323        << "pp_lumi = " << pp_lumi << endl;
0324   cout << "\t"
0325        << "AuAu_eq_lumi_C0_10 = " << AuAu_eq_lumi_C0_10 << endl;
0326   cout << "\t"
0327        << "AuAu_eq_lumi_C0_20 = " << AuAu_eq_lumi_C0_20 << endl;
0328   cout << "\t"
0329        << "AuAu_eq_lumi_C0_100 = " << AuAu_eq_lumi_C0_100 << endl;
0330   cout << "\t"
0331        << "pAu_eq_lumi_C0_100 = " << pAu_eq_lumi_C0_100 << endl;
0332
0333   TCanvas *c1 = new TCanvas("Draw_HFJetTruth_InvMass_CrossSection2RAA_Ratio" + s_suffix, "Draw_HFJetTruth_InvMass_CrossSection2RAA_Ratio" + s_suffix, 1100, 800);
0334   c1->Divide(1, 1);
0335   int idx = 1;
0336   TPad *p;
0337
0338   p = (TPad *) c1->cd(idx++);
0339   c1->Update();
0340   //  p->SetGridx(0);
0341   //  p->SetGridy(0);
0342   //  p->SetLogy();
0343
0344   // in sPH-HF-2017-001, dijet purity efficiency was conservatively assumed to be the product of that of two single b-jets,
0345   // i.e. without benifit of the enhanace of 2nd b-jet abundance after the 1st jet in the event is tagged as a b-jet
0346   TH1 *g_pp = CrossSection2RelUncert(h_b, 1., dy, pp_lumi * pp_eff * pp_purity * pp_eff * pp_purity);
0347   TH1 *g_AA = CrossSection2RelUncert(h_b, b_jet_RAA, dy, AuAu_eq_lumi_C0_10 * AuAu_eff * AuAu_purity * AuAu_eff * AuAu_purity);
0348
0349   g_pp->SetLineColor(kRed);
0350   g_AA->SetLineColor(kBlue);
0351
0352   g_pp->Draw();
0353   g_AA->Draw("same");
0354   SaveCanvas(c1, infile + "_" + TString(c1->GetName()), kTRUE);
0355
0356   TCanvas *c1 = new TCanvas("Draw_HFJetTruth_InvMass_CrossSection2RAA" + s_suffix, "Draw_HFJetTruth_InvMass_CrossSection2RAA" + s_suffix, 1100, 800);
0357   c1->Divide(1, 1);
0358   int idx = 1;
0359   TPad *p;
0360
0361   p = (TPad *) c1->cd(idx++);
0362   c1->Update();
0363
0364   p->DrawFrame(30, 0, 75, 1.3)->SetTitle(";Di-jet invariant mass [GeV/#it{c}^{2}];#it{R}^{bb}_{AA}");
0365
0366   TGraphErrors *ge_RAA = GetRAA(g_pp, g_AA);
0367
0368   ge_RAA->SetLineWidth(4);
0369   ge_RAA->SetMarkerStyle(kFullCircle);
0370   ge_RAA->SetMarkerSize(3);
0371
0372   ge_RAA->Draw("pe");
0373   ge_RAA->Print();
0374
0375   TLegend *leg = new TLegend(.0, .7, .85, .93);
0376   leg->SetFillStyle(0);
0377   leg->AddEntry("", "#it{#bf{sPHENIX}} Projection, #it{b}-jet, 0-10% Au+Au, Year 1-3", "");
0378   leg->AddEntry("", Form("|#eta_{1,2}|<%.1f, |#Delta#phi_{1,2}|>2#pi/3, p_{T,1}>15 GeV/c, p_{T,2}>10 GeV/c", dy / 2), "");
0379   leg->AddEntry("", Form("#it{p}+#it{p}: %.0f pb^{-1} samp., %.0f%% Eff., %.0f%% Pur.", pp_lumi, pp_eff * 100, pp_purity * 100), "");
0380   leg->AddEntry("", Form("Au+Au: %.0fnb^{-1} rec., %.0f%% Eff., %.0f%% Pur.", '%', AuAu_MB_Evt/6.8252 / 1e9, AuAu_eff * 100, AuAu_purity * 100), "");
0381   leg->Draw();
0382
0383 //  TLegend *leg2 = new TLegend(.17, .70, .88, .77);
0384 //  leg2->AddEntry(ge_RAA, "#it{b}-jet #it{R}_{AA}, Au+Au 0-10%C, #sqrt{s_{NN}}=200 GeV", "pl");
0385 //  leg2->Draw();
0386
0387   SaveCanvas(c1, infile + "_" + TString(c1->GetName()), kTRUE);
0388
0389   TCanvas *c1 = new TCanvas("Draw_HFJetTruth_InvMass_CrossSection2RAA_Theory" + s_suffix, "Draw_HFJetTruth_InvMass_CrossSection2RAA_Theory" + s_suffix, 1100, 800);
0390   c1->Divide(1, 1);
0391   int idx = 1;
0392   TPad *p;
0393
0394   p = (TPad *) c1->cd(idx++);
0395   c1->Update();
0396
0397   p->DrawFrame(35, 0, 75, 1.4)->SetTitle(";Di-jet invariant mass [GeV];#it{R}^{bb}_{AA}");
0398
0399   //
0400   TGraph *g18 = GetRAAKang2019(1, 1.8);
0401   TGraph *g20 = GetRAAKang2019(1, 2.0);
0402   TGraph *g22 = GetRAAKang2019(1, 2.2);
0403   //
0404   g18->SetLineColor(kBlue + 3);
0405   g20->SetLineColor(kGreen + 3);
0406   g22->SetLineColor(kRed + 3);
0407   g18->SetLineWidth(3);
0408   g20->SetLineWidth(3);
0409   g22->SetLineWidth(3);
0410   //
0411   g18->Draw("l");
0412   g20->Draw("l");
0413   g22->Draw("l");
0414   //
0415   ge_RAA->DrawClone("pe");
0416   leg->DrawClone();
0417 //  leg2->DrawClone();
0418   //
0419   TLegend *leg1 = new TLegend(.19, .6, .92, .7);
0420 //  leg1->AddEntry("", "Kang et al, PRD #bf{99}, 034006 (2019), m=m_{b}, rad.+col.", "");
0421   leg1->SetHeader("Kang et al, PRD #bf{99}, 034006 (2019), m=m_{b}, rad.+col.");
0422   leg1->Draw();
0423   TLegend *leg11 = new TLegend(.35, .5, .55, .62);
0424   leg11->AddEntry(g18, "g_{med} = 1.8", "l");
0425   leg11->AddEntry(g20, "g_{med} = 2.0", "l");
0426   leg11->AddEntry(g22, "g_{med} = 2.2", "l");
0427   leg11->Draw();
0428
0429   SaveCanvas(c1, infile + "_" + TString(c1->GetName()), kTRUE);
0430
0431
0432   TCanvas *c1 = new TCanvas("Draw_HFJetTruth_InvMass_CrossSection2RAARatio_Theory" + s_suffix, "Draw_HFJetTruth_InvMass_CrossSection2RAARatio_Theory" + s_suffix, 1100, 800);
0433   c1->Divide(1, 1);
0434   int idx = 1;
0435   TPad *p;
0436
0437   p = (TPad *) c1->cd(idx++);
0438   c1->Update();
0439
0440   p->DrawFrame(35, 0, 75, 15)->SetTitle(";Di-jet invariant mass [GeV];#it{R}^{bb}_{AA}/#it{R}^{jj}_{AA}");
0441   TLine * l = new TLine(35, 1, 75, 1);
0442 //  l->SetLineStyle(kDashed);
0443   l->Draw();
0444
0445   //
0446   TGraph *g18 = GetRAARatioKang2019(1, 1.8);
0447   TGraph *g20 = GetRAARatioKang2019(1, 2.0);
0448   TGraph *g22 = GetRAARatioKang2019(1, 2.2);
0449   //
0450   g18->SetLineColor(kBlue + 3);
0451   g20->SetLineColor(kGreen + 3);
0452   g22->SetLineColor(kRed + 3);
0453   g18->SetLineWidth(3);
0454   g20->SetLineWidth(3);
0455   g22->SetLineWidth(3);
0456   //
0457   g18->Draw("l");
0458   g20->Draw("l");
0459   g22->Draw("l");
0460   //
0461   TGraphErrors * ge_RAARatio = (TGraphErrors *) ge_RAA->DrawClone("pe");
0462   assert(ge_RAARatio);
0463   for (int bin = 0; bin<ge_RAARatio->GetN();++bin)
0464   {
0465     ge_RAARatio->GetY()[bin] *= target_RAA_ratio/b_jet_RAA;
0466     ge_RAARatio->GetEY()[bin] *= target_RAA_ratio/b_jet_RAA;
0467   }
0468
0469   leg->DrawClone();
0470 //  leg2->DrawClone();
0471   //
0472   TLegend *leg1 = new TLegend(.2, .6, .92, .70);
0473 //  leg1->AddEntry("", "Kang et al, PRD #bf{99}, 034006 (2019), m=m_{b}, rad.+col.", "");
0474   leg1->SetHeader("Kang et al, PRD #bf{99}, 034006 (2019), m=m_{b}, rad.+col.");
0475   leg1->Draw();
0476   TLegend *leg11 = new TLegend(.35, .5, .55, .62);
0477   leg11->AddEntry(g18, "g_{med} = 1.8", "l");
0478   leg11->AddEntry(g20, "g_{med} = 2.0", "l");
0479   leg11->AddEntry(g22, "g_{med} = 2.2", "l");
0480   leg11->Draw();
0481
0482   SaveCanvas(c1, infile + "_" + TString(c1->GetName()), kTRUE);
0483
0484
0485 }
0486
0487 TGraphErrors *GetRAA(TH1 *h_pp, TH1 *h_AA)
0488 {
0489   int n_bin = 0;
0490
0491   double xs[1000] = {0};
0492   double ys[1000] = {0};
0493   double eys[1000] = {0};
0494
0495   assert(h_pp);
0496   assert(h_AA);
0497   assert(h_pp->GetNbinsX() == h_AA->GetNbinsX());
0498
0499   for (int i = 1; i <= h_pp->GetNbinsX(); ++i)
0500   {
0501     if (h_pp->GetBinError(i) > 0 && h_pp->GetBinError(i) < 1 && h_AA->GetBinError(i) > 0 && h_AA->GetBinError(i) < 1)
0502     {
0503       xs[n_bin] = h_pp->GetXaxis()->GetBinCenter(i);
0504
0505       ys[n_bin] = h_AA->GetBinContent(i) / h_pp->GetBinContent(i);
0506
0507       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));
0508
0509       // final uncertainty cut
0510       if (eys[n_bin] < .4)
0511         n_bin += 1;
0512     }
0513   }
0514
0515   TGraphErrors *ge = new TGraphErrors(n_bin, xs, ys, NULL, eys);
0516   ge->SetName(TString("RAA_") + h_AA->GetName());
0517
0518   return ge;
0519 }
0520
0521 TH1 *CrossSection2RelUncert(const TH1F *h_cross,
0522                             const double suppression,
0523                             const double deta,
0524                             const double pp_quiv_int_lum)
0525 {
0526   assert(h_cross);
0527   TH1 *
0528       h_ratio = (TH1 *)
0529                     h_cross->Clone(TString(h_cross->GetName()) + Form("_copy%d", rand()));
0530
0531   //convert to count per bin
0532   h_ratio->Scale(deta * deta * h_ratio->GetXaxis()->GetBinWidth(0) * pp_quiv_int_lum * suppression);
0533
0534   // define the x-binning
0535   h_ratio->Rebin(5);
0536
0537   for (int i = 1; i <= h_ratio->GetNbinsX(); ++i)
0538   {
0539     const double yield = h_ratio->GetBinContent(i);
0540
0541     if (yield > 0)
0542     {
0543       h_ratio->SetBinContent(i, suppression);
0544
0545       h_ratio->SetBinError(i, suppression / sqrt(yield));
0546     }
0547     else
0548     {
0549       h_ratio->SetBinContent(i, 0);
0550
0551       h_ratio->SetBinError(i, 0);
0552     }
0553   }
0554
0555   h_ratio->GetYaxis()->SetTitle("Relative Cross Section and Uncertainty");
0556
0557   return h_ratio;
0558 }
0559
0560 void Convert2CrossSection(TH1 *h, const double int_lumi, const double dy)
0561 {
0562   h->Sumw2();
0563
0564   for (int i = 1; i <= h->GetXaxis()->GetNbins(); ++i)
0565   {
0566     const double m1 = h->GetXaxis()->GetBinLowEdge(i);
0567     const double m2 = h->GetXaxis()->GetBinUpEdge(i);
0568
0569     //      const double dpT2 =  pT2*pT2 - pT1*pT1;
0570     const double dm = m2 - m1;
0571
0572     //      const double count2sigma = 1./dpT2/dy/int_lumi;
0573     const double count2sigma = 1. / dm / int_lumi / dy / dy;
0574
0575     h->SetBinContent(i, h->GetBinContent(i) * count2sigma);
0576     h->SetBinError(i, h->GetBinError(i) * count2sigma);
0577   }
0578
0579   //  h->GetYaxis()->SetTitle("#sigma/(dp_{T}^{2} dy) (pb/(GeV/c)^{2})");
0580
0581   h->SetMarkerStyle(kFullCircle);
0582 }
0583
0584 TGraphAsymmErrors *GetSTAR2019(const TString hepdata = "/sphenix/user/jinhuang/HF-jet/event_gen/HEPData-ins1493842-v1-Table_4.root")
0585 {
0586   ////  https://doi.org/10.17182/hepdata.77208.v1/t4
0587   //  Measurement of the cross section and longitudinal double-spin asymmetry for di-jet production in polarized pp collisions at s = 200 GeV
0588   //
0589   //  The STAR collaboration
0590   //  Adamczyk, L. , Adkins, J.K. , Agakishiev, G. , Aggarwal, M.M. , Ahammed, Z. , Alekseev, I. , Anderson, D.M. , Aoyama, R. , Aparin, A. , Arkhipkin, D.
0591   //  No Journal Information, 2016
0592   //  https://doi.org/10.17182/hepdata.77208
0593
0594   TFile *_file0 = TFile::Open(hepdata);
0595   _file0->cd("Table 4");
0596   TGraphAsymmErrors *gae = (TGraphAsymmErrors *)
0597                                gDirectory->GetObjectChecked("Graph1D_y1", "TGraphAsymmErrors");
0598   assert(gae);
0599
0600   for (int bin = 0; bin < gae->GetN(); ++bin)
0601   {
0602     (gae->GetY())[bin] *= 1e6;                                // ub -> pb
0603     gae->SetPointEYhigh(bin, gae->GetErrorYhigh(bin) * 1e6);  // ub -> pb
0604     gae->SetPointEYlow(bin, gae->GetErrorYlow(bin) * 1e6);    // ub -> pb
0605   }
0606   gae->SetName("gaeSTARDiJetMass");
0607
0608   return gae;
0609 }
0610
0611 TGraph *GetRAAKang2019(const int mass = 1, const double g = 1.8)
0612 {
0613   //  m12    2mb,g=1.8    1mb,g=1.8    2mb,g=2.0    1mb,g=2.0    2mb,g=2.2    1mb,g=2.2
0614   //  17     0.959318     0.813631     0.937975     0.718446     0.904169     0.583889
0615   //  19     0.924265     0.734466     0.887045     0.618448     0.831249     0.469604
0616   //  21     0.889933     0.676922     0.838959     0.550413     0.765768     0.397484
0617   //  23     0.852022     0.624044      0.78657     0.489034     0.696008     0.335986
0618   //  25     0.813713     0.572249     0.736696     0.435365     0.634859     0.287682
0619   //  27     0.829553     0.588906     0.752964     0.444663     0.647904     0.287249
0620   //  29     0.787303     0.554276     0.702613     0.412568     0.592341      0.26188
0621   //  31     0.767304     0.519823     0.674582     0.374897     0.555677      0.22774
0622   //  33     0.774209     0.531997      0.68111     0.383097     0.559728       0.2299
0623   //  35     0.733919     0.503672     0.633725     0.359385     0.509993     0.213934
0624   //  37     0.763963     0.528822     0.665819     0.373186     0.537606     0.214043
0625   //  39     0.741506     0.564112     0.645253     0.416396      0.52631     0.253036
0626   //  41      0.74574     0.497485     0.638774     0.336391     0.500389     0.179804
0627   //  43     0.725435     0.496459     0.617035     0.336753     0.479825     0.179778
0628   //  45     0.714991     0.472252     0.599281     0.309914     0.453767      0.15839
0629   //  47     0.680055     0.443568     0.559106     0.287112     0.413458     0.145234
0630   //  49     0.679489     0.450226     0.557737      0.28911     0.410027     0.142904
0631   //  51     0.652643     0.420073     0.525378     0.262171     0.375742     0.124363
0632   //  53     0.648852     0.421612     0.520426     0.260902     0.368284     0.121749
0633   //  55     0.659723     0.424185     0.522925     0.255923      0.36165     0.114938
0634   //  57     0.608096     0.394726     0.474039      0.23987     0.324512     0.109253
0635   //  59     0.599733     0.385379     0.462397     0.228792     0.309607     0.100507
0636   //  61     0.602495     0.398467     0.466165     0.238347     0.312724     0.103957
0637   //  63      0.58698     0.369703     0.441268      0.20886     0.281604    0.0844847
0638   //  65      0.58804      0.39896     0.449849     0.238734     0.296312     0.102839
0639   //  67     0.568857     0.369444     0.423411     0.208794     0.264683    0.0833302
0640   //  69     0.569308     0.374558     0.421737     0.210743     0.259612    0.0834087
0641   //  71     0.570937     0.379757     0.418615     0.213217     0.253412    0.0841172
0642   //  73     0.571426     0.368965     0.404248     0.201111     0.234995    0.0773808
0643   //  75     0.444747     0.285571     0.308362     0.153737      0.17623    0.0575338
0644   //  77     0.518565     0.353508     0.368773     0.198289     0.215024    0.0772062
0645   //  79      0.52314     0.354647      0.36598     0.194516     0.207028    0.0732502
0646   //  81     0.462486     0.306947     0.313368     0.163098     0.169892    0.0593907
0647   //  83     0.494343     0.338822     0.337503     0.182332     0.182904    0.0665247
0648   //  85     0.477053     0.321068     0.315372     0.167256     0.163585    0.0589566
0649   //  87     0.425877     0.287827     0.276778      0.14819     0.140167    0.0513555
0650   //  89     0.449382     0.310267     0.291343      0.16052     0.146033    0.0555645
0651   //  91     0.426898     0.294991     0.270933     0.150014     0.131644    0.0508301
0652   //  93      0.58802     0.451299     0.406374     0.248587     0.210967    0.0874347
0653   //  95     0.367013     0.245387     0.216773      0.11698    0.0963793    0.0374309
0654   //  97     0.363587     0.251356     0.213642     0.120293    0.0932188    0.0381195
0655   //  99     0.349846     0.240864     0.197647     0.111828    0.0824488    0.0345796
0656   //  101     0.294855      0.20189     0.158814    0.0907104    0.0634189    0.0274697
0657   //  103     0.264038      0.18211     0.137505    0.0801353    0.0531663    0.0239006
0658
0659   const int n_Data = 44;
0660
0661   const double m12[] = {
0662       17,
0663       19,
0664       21,
0665       23,
0666       25,
0667       27,
0668       29,
0669       31,
0670       33,
0671       35,
0672       37,
0673       39,
0674       41,
0675       43,
0676       45,
0677       47,
0678       49,
0679       51,
0680       53,
0681       55,
0682       57,
0683       59,
0684       61,
0685       63,
0686       65,
0687       67,
0688       69,
0689       71,
0690       73,
0691       75,
0692       77,
0693       79,
0694       81,
0695       83,
0696       85,
0697       87,
0698       89,
0699       91,
0700       93,
0701       95,
0702       97,
0703       99,
0704       101,
0705       103};
0706
0707   const double RAA_2mb_g18[] = {
0708       0.959318,
0709       0.924265,
0710       0.889933,
0711       0.852022,
0712       0.813713,
0713       0.829553,
0714       0.787303,
0715       0.767304,
0716       0.774209,
0717       0.733919,
0718       0.763963,
0719       0.741506,
0720       0.74574,
0721       0.725435,
0722       0.714991,
0723       0.680055,
0724       0.679489,
0725       0.652643,
0726       0.648852,
0727       0.659723,
0728       0.608096,
0729       0.599733,
0730       0.602495,
0731       0.58698,
0732       0.58804,
0733       0.568857,
0734       0.569308,
0735       0.570937,
0736       0.571426,
0737       0.444747,
0738       0.518565,
0739       0.52314,
0740       0.462486,
0741       0.494343,
0742       0.477053,
0743       0.425877,
0744       0.449382,
0745       0.426898,
0746       0.58802,
0747       0.367013,
0748       0.363587,
0749       0.349846,
0750       0.294855,
0751       0.264038};
0752
0753   const double RAA_1mb_g_18[] = {
0754       0.813631,
0755       0.734466,
0756       0.676922,
0757       0.624044,
0758       0.572249,
0759       0.588906,
0760       0.554276,
0761       0.519823,
0762       0.531997,
0763       0.503672,
0764       0.528822,
0765       0.564112,
0766       0.497485,
0767       0.496459,
0768       0.472252,
0769       0.443568,
0770       0.450226,
0771       0.420073,
0772       0.421612,
0773       0.424185,
0774       0.394726,
0775       0.385379,
0776       0.398467,
0777       0.369703,
0778       0.39896,
0779       0.369444,
0780       0.374558,
0781       0.379757,
0782       0.368965,
0783       0.285571,
0784       0.353508,
0785       0.354647,
0786       0.306947,
0787       0.338822,
0788       0.321068,
0789       0.287827,
0790       0.310267,
0791       0.294991,
0792       0.451299,
0793       0.245387,
0794       0.251356,
0795       0.240864,
0796       0.20189,
0797       0.18211};
0798
0799   const double RAA_2mb_g_20[] = {
0800       0.937975,
0801       0.887045,
0802       0.838959,
0803       0.78657,
0804       0.736696,
0805       0.752964,
0806       0.702613,
0807       0.674582,
0808       0.68111,
0809       0.633725,
0810       0.665819,
0811       0.645253,
0812       0.638774,
0813       0.617035,
0814       0.599281,
0815       0.559106,
0816       0.557737,
0817       0.525378,
0818       0.520426,
0819       0.522925,
0820       0.474039,
0821       0.462397,
0822       0.466165,
0823       0.441268,
0824       0.449849,
0825       0.423411,
0826       0.421737,
0827       0.418615,
0828       0.404248,
0829       0.308362,
0830       0.368773,
0831       0.36598,
0832       0.313368,
0833       0.337503,
0834       0.315372,
0835       0.276778,
0836       0.291343,
0837       0.270933,
0838       0.406374,
0839       0.216773,
0840       0.213642,
0841       0.197647,
0842       0.158814,
0843       0.137505};
0844
0845   const double RAA_1mb_g_20[] = {
0846       0.718446,
0847       0.618448,
0848       0.550413,
0849       0.489034,
0850       0.435365,
0851       0.444663,
0852       0.412568,
0853       0.374897,
0854       0.383097,
0855       0.359385,
0856       0.373186,
0857       0.416396,
0858       0.336391,
0859       0.336753,
0860       0.309914,
0861       0.287112,
0862       0.28911,
0863       0.262171,
0864       0.260902,
0865       0.255923,
0866       0.23987,
0867       0.228792,
0868       0.238347,
0869       0.20886,
0870       0.238734,
0871       0.208794,
0872       0.210743,
0873       0.213217,
0874       0.201111,
0875       0.153737,
0876       0.198289,
0877       0.194516,
0878       0.163098,
0879       0.182332,
0880       0.167256,
0881       0.14819,
0882       0.16052,
0883       0.150014,
0884       0.248587,
0885       0.11698,
0886       0.120293,
0887       0.111828,
0888       0.0907104,
0889       0.0801353};
0890
0891   const double RAA_2mb_g_22[] = {
0892       0.904169,
0893       0.831249,
0894       0.765768,
0895       0.696008,
0896       0.634859,
0897       0.647904,
0898       0.592341,
0899       0.555677,
0900       0.559728,
0901       0.509993,
0902       0.537606,
0903       0.52631,
0904       0.500389,
0905       0.479825,
0906       0.453767,
0907       0.413458,
0908       0.410027,
0909       0.375742,
0910       0.368284,
0911       0.36165,
0912       0.324512,
0913       0.309607,
0914       0.312724,
0915       0.281604,
0916       0.296312,
0917       0.264683,
0918       0.259612,
0919       0.253412,
0920       0.234995,
0921       0.17623,
0922       0.215024,
0923       0.207028,
0924       0.169892,
0925       0.182904,
0926       0.163585,
0927       0.140167,
0928       0.146033,
0929       0.131644,
0930       0.210967,
0931       0.0963793,
0932       0.0932188,
0933       0.0824488,
0934       0.0634189,
0935       0.0531663};
0936
0937   const double RAA_1mb_g_22[] = {
0938       0.583889,
0939       0.469604,
0940       0.397484,
0941       0.335986,
0942       0.287682,
0943       0.287249,
0944       0.26188,
0945       0.22774,
0946       0.2299,
0947       0.213934,
0948       0.214043,
0949       0.253036,
0950       0.179804,
0951       0.179778,
0952       0.15839,
0953       0.145234,
0954       0.142904,
0955       0.124363,
0956       0.121749,
0957       0.114938,
0958       0.109253,
0959       0.100507,
0960       0.103957,
0961       0.0844847,
0962       0.102839,
0963       0.0833302,
0964       0.0834087,
0965       0.0841172,
0966       0.0773808,
0967       0.0575338,
0968       0.0772062,
0969       0.0732502,
0970       0.0593907,
0971       0.0665247,
0972       0.0589566,
0973       0.0513555,
0974       0.0555645,
0975       0.0508301,
0976       0.0874347,
0977       0.0374309,
0978       0.0381195,
0979       0.0345796,
0980       0.0274697,
0981       0.0239006};
0982
0983   TGraph *gRAAKang2019(NULL);
0984
0985   if (mass == 1)
0986   {
0987     if (g == 1.8)
0988     {
0989       gRAAKang2019 = new TGraph(n_Data, m12, RAA_1mb_g_18);
0990     }
0991     else if (g == 2)
0992     {
0993       gRAAKang2019 = new TGraph(n_Data, m12, RAA_1mb_g_20);
0994     }
0995     else if (g == 2.2)
0996     {
0997       gRAAKang2019 = new TGraph(n_Data, m12, RAA_1mb_g_22);
0998     }
0999   }
1000   else if (mass == 2)
1001   {
1002     if (g == 1.8)
1003     {
1004       gRAAKang2019 = new TGraph(n_Data, m12, RAA_2mb_g_18);
1005     }
1006     else if (g == 2)
1007     {
1008       gRAAKang2019 = new TGraph(n_Data, m12, RAA_2mb_g_20);
1009     }
1010     else if (g == 2.2)
1011     {
1012       gRAAKang2019 = new TGraph(n_Data, m12, RAA_2mb_g_22);
1013     }
1014   }
1015
1016   cout <<"mass = "<<mass<<", g="<<g<<endl;
1017   assert(gRAAKang2019);
1018   gRAAKang2019->Print();
1019
1020   return gRAAKang2019;
1021 }
1022
1023 TGraph *GetRAARatioKang2019(const int mass = 1, const double g = 1.8)
1024 {
1025   //  m12     g=1.8        g=2.0        g=2.2
1026   //  17      4.10724      7.40351      16.6273
1027   //  19      4.44519      7.83383      18.0237
1028   //  21      4.57977       7.8564      17.5948
1029   //  23      4.29797      7.07926      15.2511
1030   //  25      3.58437      5.58706       11.445
1031   //  27      3.81735      5.93123      12.1819
1032   //  29      2.88061      4.20894      8.22601
1033   //  31       3.0618       4.4508      8.61545
1034   //  33      2.94732      4.17463      7.85741
1035   //  35      2.92894      4.20026       8.0066
1036   //  37      2.70991      3.76532      6.85238
1037   //  39      2.91857      4.25505      8.23559
1038   //  41      2.39766      3.14728      5.31672
1039   //  43      2.41189      3.19676      5.40647
1040   //  45       2.0471      2.56789      4.11909
1041   //  47      2.04743      2.57154      4.09681
1042   //  49      1.97205      2.43492       3.7728
1043   //  51      1.67187      1.96796      2.89214
1044   //  53      1.75599       2.0696      2.98919
1045   //  55      1.69526      1.94817      2.71746
1046   //  57      1.02735      1.09074      1.46593
1047   //  59      1.56319      1.76551      2.39532
1048   //  61      1.60436      1.82217      2.45074
1049   //  63      1.39696      1.47696      1.82356
1050   //  65      1.45008      1.62795      2.14264
1051   //  67      1.36533      1.43573      1.73491
1052   //  69      1.40584       1.4939      1.80942
1053   //  71      1.38037      1.45466      1.75179
1054   //  73      1.27568      1.28553      1.49002
1055   //  75      1.03127      1.04286      1.19033
1056   //  77      1.20987      1.25797      1.47642
1057   //  79      1.27784      1.31767      1.51129
1058   //  81      1.06928      1.06112       1.1719
1059   //  83      1.22018       1.2414        1.386
1060   //  85      1.15586      1.13797      1.23096
1061   //  87      1.00112     0.966856      1.02192
1062   //  89      1.07953      1.05446       1.1226
1063   //  91      1.08107      1.06196       1.1313
1064   //  93       1.6204      1.70796      1.87276
1065   //  95     0.905063     0.833268     0.839565
1066   //  97      0.94488     0.884202     0.892716
1067   //  99     0.956238     0.882098     0.880081
1068   //  101     0.855383     0.777236     0.768626
1069   //  103     0.832522     0.754264     0.742384
1070
1071   const int n_Data = 44;
1072
1073   const double m12[] = {
1074       17,
1075       19,
1076       21,
1077       23,
1078       25,
1079       27,
1080       29,
1081       31,
1082       33,
1083       35,
1084       37,
1085       39,
1086       41,
1087       43,
1088       45,
1089       47,
1090       49,
1091       51,
1092       53,
1093       55,
1094       57,
1095       59,
1096       61,
1097       63,
1098       65,
1099       67,
1100       69,
1101       71,
1102       73,
1103       75,
1104       77,
1105       79,
1106       81,
1107       83,
1108       85,
1109       87,
1110       89,
1111       91,
1112       93,
1113       95,
1114       97,
1115       99,
1116       101,
1117       103};
1118
1119   const double RAA_1mb_g_18[] = {
1120       4.10724,
1121       4.44519,
1122       4.57977,
1123       4.29797,
1124       3.58437,
1125       3.81735,
1126       2.88061,
1127       3.0618,
1128       2.94732,
1129       2.92894,
1130       2.70991,
1131       2.91857,
1132       2.39766,
1133       2.41189,
1134       2.0471,
1135       2.04743,
1136       1.97205,
1137       1.67187,
1138       1.75599,
1139       1.69526,
1140       1.02735,
1141       1.56319,
1142       1.60436,
1143       1.39696,
1144       1.45008,
1145       1.36533,
1146       1.40584,
1147       1.38037,
1148       1.27568,
1149       1.03127,
1150       1.20987,
1151       1.27784,
1152       1.06928,
1153       1.22018,
1154       1.15586,
1155       1.00112,
1156       1.07953,
1157       1.08107,
1158       1.6204,
1159       0.905063,
1160       0.94488,
1161       0.956238,
1162       0.855383,
1163       0.832522};
1164
1165   const double RAA_1mb_g_20[] = {
1166       7.40351,
1167       7.83383,
1168       7.8564,
1169       7.07926,
1170       5.58706,
1171       5.93123,
1172       4.20894,
1173       4.4508,
1174       4.17463,
1175       4.20026,
1176       3.76532,
1177       4.25505,
1178       3.14728,
1179       3.19676,
1180       2.56789,
1181       2.57154,
1182       2.43492,
1183       1.96796,
1184       2.0696,
1185       1.94817,
1186       1.09074,
1187       1.76551,
1188       1.82217,
1189       1.47696,
1190       1.62795,
1191       1.43573,
1192       1.4939,
1193       1.45466,
1194       1.28553,
1195       1.04286,
1196       1.25797,
1197       1.31767,
1198       1.06112,
1199       1.2414,
1200       1.13797,
1201       .966856,
1202       1.05446,
1203       1.06196,
1204       1.70796,
1205       .833268,
1206       .884202,
1207       .882098,
1208       0.777236,
1209       0.754264};
1210
1211   const double RAA_1mb_g_22[] = {
1212       16.6273,
1213       18.0237,
1214       17.5948,
1215       15.2511,
1216       11.445,
1217       12.1819,
1218       8.22601,
1219       8.61545,
1220       7.85741,
1221       8.0066,
1222       6.85238,
1223       8.23559,
1224       5.31672,
1225       5.40647,
1226       4.11909,
1227       4.09681,
1228       3.7728,
1229       2.89214,
1230       2.98919,
1231       2.71746,
1232       1.46593,
1233       2.39532,
1234       2.45074,
1235       1.82356,
1236       2.14264,
1237       1.73491,
1238       1.80942,
1239       1.75179,
1240       1.49002,
1241       1.19033,
1242       1.47642,
1243       1.51129,
1244       1.1719,
1245       1.386,
1246       1.23096,
1247       1.02192,
1248       1.1226,
1249       1.1313,
1250       1.87276,
1251       0.839565,
1252       0.892716,
1253       0.880081,
1254       0.768626,
1255       0.742384};
1256
1257   TGraph *gRAARatioKang2019(NULL);
1258
1259   if (mass == 1)
1260   {
1261     if (g == 1.8)
1262     {
1263       gRAARatioKang2019 = new TGraph(n_Data, m12, RAA_1mb_g_18);
1264     }
1265     else if (g == 2)
1266     {
1267       gRAARatioKang2019 = new TGraph(n_Data, m12, RAA_1mb_g_20);
1268     }
1269     else if (g == 2.2)
1270     {
1271       gRAARatioKang2019 = new TGraph(n_Data, m12, RAA_1mb_g_22);
1272     }
1273   }
1274
1275   assert(gRAARatioKang2019);
1276
1277   return gRAARatioKang2019;
1278 }