HFDiJetMomImbalance/macros/Draw_BDiJetImbalance.C

0001 ////////////////////////////////////////////////////////////////////////////////
0002 //
0003 // Draw the Di b-jet momentum imbalance
0004 //
0005 // Useful references:
0006 // CMS di b-jet: http://cds.cern.ch/record/2202805/files/HIN-16-005-pas.pdf
0007 //
0008 ////////////////////////////////////////////////////////////////////////////////
0009 //
0010 // Darren McGlinchey
0011 // 6 Jan 2017
0012 //
0013 ////////////////////////////////////////////////////////////////////////////////
0014
0015 #include <TFile.h>
0016 #include <TChain.h>
0017 #include <TTree.h>
0018 #include <TH1.h>
0019 #include <TCanvas.h>
0020 #include <TLatex.h>
0021 #include <TStyle.h>
0022 #include <TMath.h>
0023 #include <TLegend.h>
0024 #include <TRandom3.h>
0025 #include <TString.h>
0026 #include <TSystem.h>
0027
0028 #include <iostream>
0029
0030 using namespace std;
0031
0032 void Draw_BDiJetImbalance()
0033 {
0034   gStyle->SetOptStat(0);
0035   gStyle->SetOptTitle(0);
0036
0037   //==========================================================================//
0038   // SET RUNNING CONDITIONS
0039   //==========================================================================//
0040   bool print_plots = true;
0041
0042   bool is_pp = false;
0043
0044   // const char* inFile = "HFtag_bjet.root";
0045   // const char* inFile = "HFtag_bjet_pp.root";
0046   // const char* inFile = "/phenix/plhf/dcm07e/sPHENIX/bjetsims/test.root";
0047   // if (!is_pp)
0048   //   const char* inFile = "HFtag_bjet_AuAu0-10.root";
0049
0050   const char* inDir = "/phenix/plhf/dcm07e/sPHENIX/bjetsims/ana";
0051
0052   // // desired nn integrated luminosity [pb^{-1}]
0053   // double lumiDesired = 175;
0054   // if ( !is_pp )
0055   //   lumiDesired = 229; // AuAu 0-10% central (10B events)
0056   // desired nn integrated luminosity [pb^{-1}]
0057   double lumiDesired = 200;
0058   TString scol("p + p #sqrt{s_{_{NN}}} = 200 GeV");
0059   TString slumi(Form("#int L dt = %.0f pb^{-1} |z| < 10 cm", lumiDesired));
0060   if ( !is_pp )
0061   {
0062     // equivelant nn luminosity:
0063     // N_{evt}^{NN} / sigma_NN = N_{evt}^{AA}*Ncoll / sigma_NN
0064     // = 10e9 * 962 / 42e9 pb = 229 pb^-1
0065     double evntDesired = 24e9; // AuAu 0-10% central
0066     lumiDesired = evntDesired * 962 / 42e9; // AuAu 0-10% central
0067     scol = "0-10% Au + Au #sqrt{s_{_{NN}}}=200 GeV";
0068     slumi = Form("%.0fB events |z| < 10 cm", evntDesired / 1.e9);
0069   }
0070
0071   // luminosity per file generated [pb^{-1}]
0072   double lumiPerFile = 99405 / (4.641e-06 * 1e12);
0073   double lumiRead = 0;
0074
0075   double pT1min = 20;
0076   double pT2min = 10;
0077   // double etamax = 1.0;
0078   double etamax = 0.7;
0079
0080   // double bJetEff = 0.5; // b-jet tagging efficiency
0081   double bJetEff = 0.60; // p+p b-jet tagging efficiency
0082   if ( !is_pp )
0083     bJetEff = 0.40;
0084
0085   double bJetPur = 0.40; // b-jet tagging purity
0086
0087   double RAA = 1.0;
0088   if (!is_pp)
0089     RAA = 0.6;
0090
0091   const int NETACUT = 3;
0092   double etacut[] = {1.3, 1.0, 0.7};
0093   double etaColor[] = {kBlue, kRed, kBlack};
0094
0095   const int PTCUT = 3;
0096   double ptcut[] = {10, 20, 40};
0097   double ptColor[] = {kBlue, kRed, kBlack};
0098
0099   //==========================================================================//
0100   // DECLARE VARIABLES
0101   //==========================================================================//
0102
0103   //-- tree variables
0104   TChain *ttree;
0105   Int_t           event;
0106   Int_t           truthjet_n;
0107   Int_t           truthjet_parton_flavor[100];
0108   Int_t           truthjet_hadron_flavor[100];
0109   Float_t         truthjet_pt[100];
0110   Float_t         truthjet_eta[100];
0111   Float_t         truthjet_phi[100];
0112
0113   //-- histograms
0114
0115   TH1D* hjet_pT = new TH1D("hjet_pT", ";p_{T}^{jet} [GeV/c]", 20, 0, 100);
0116   TH1D* hjet_eta = new TH1D("hjet_eta", ";#eta^{jet}", 40, -2., 2.);
0117   TH1D* hjet_phi = new TH1D("hjet_phi", ";#phi^{jet}", 40, -4, 4);
0118
0119   TH1D* hdijet_pT1 = new TH1D("hdijet_pT1", ";p_{T,1}^{jet} [GeV/c]", 20, 0, 100);
0120   TH1D* hdijet_pT2 = new TH1D("hdijet_pT2", ";p_{T,2}^{jet} [GeV/c]", 20, 0, 100);
0121   TH1D* hdijet_xj = new TH1D("hdijet_xj", ";x_{j} = p_{T,2} / p_{T,1}; event fraction", 10, 0, 1);
0122   TH1D* hdijet_dphi = new TH1D("hdijet_dphi", ";|#Delta#phi_{12}|; event fraction", 10, 0, TMath::Pi());
0123
0124   TH1D* hdijet_sing_pT1 = new TH1D("hdijet_sing_pT1",
0125                                    ";p_{T,1} [GeV/c];Fraction of Dijet / Single Accepted",
0126                                    20, 00, 100);
0127   TH1D* hdijet_sing_eta2 = new TH1D("hdijet_sing_eta2",
0128                                     ";#eta_{2};Fraction of Dijet / Single Accepted",
0129                                     15, -1.5, 1.5);
0130   TH1D* hdijet_eff_pT1[3];
0131   TH1D* hdijet_eff_eta2[3];
0132   for (int i = 0; i < NETACUT; i++)
0133   {
0134     hdijet_eff_pT1[i] = new TH1D(Form("hdijet_eff_pT1_%i", i),
0135                                  ";p_{T,1} [GeV/c];Fraction of Dijet / Single Accepted",
0136                                  20, 00, 100);
0137     hdijet_eff_pT1[i]->SetLineColor(etaColor[i]);
0138     hdijet_eff_pT1[i]->SetLineWidth(2);
0139   } // i
0140   for (int i = 0; i < PTCUT; i++)
0141   {
0142     hdijet_eff_eta2[i] = new TH1D(Form("hdijet_eff_eta2_%i", i),
0143                                   ";#eta;Fraction of Dijet / Single Accepted",
0144                                   15, -1.5, 1.5);
0145     hdijet_eff_eta2[i]->SetLineColor(etaColor[i]);
0146     hdijet_eff_eta2[i]->SetLineWidth(2);
0147   } // i
0148
0149   hjet_pT->Sumw2();
0150   hjet_eta->Sumw2();
0151   hjet_phi->Sumw2();
0152
0153   hdijet_pT1->Sumw2();
0154   hdijet_pT2->Sumw2();
0155   hdijet_dphi->Sumw2();
0156   hdijet_xj->Sumw2();
0157
0158   // for fixing the uncertainties
0159   TH1D* hdijet_xj_fix;
0160
0161   //-- other
0162   TRandom3 *rand = new TRandom3();
0163
0164   bool acc[100];
0165
0166   //==========================================================================//
0167   // LOAD TTREE
0168   //==========================================================================//
0169   // cout << endl;
0170   // cout << "--> Reading data from " << inFile << endl;
0171
0172   // TFile *fin = TFile::Open(inFile);
0173   // if (!fin)
0174   // {
0175   //   cout << "ERROR!! Unable to open " << inFile << endl;
0176   //   return;
0177   // }
0178
0179   // ttree = (TTree*) fin->Get("ttree");
0180   // if (!ttree)
0181   // {
0182   //   cout << "ERROR!! Unable to find ttree in " << inFile << endl;
0183   //   return;
0184   // }
0185
0186   cout << endl;
0187   cout << "--> Reading data from dir: " << inDir << endl;
0188
0189   // calculate the number of files necessary for the desired luminosity
0190   int nfiles = (int)(lumiDesired / lumiPerFile);
0191   cout << "  desired luminosity : " << lumiDesired << endl;
0192   cout << "  luminosity per file: " << lumiPerFile << endl;
0193   cout << "  N files needed     : " << nfiles << endl;
0194
0195
0196   // check the number of files found in the directory
0197   int nfound = (gSystem->GetFromPipe(Form("ls %s/*.root | wc -l", inDir))).Atoi();
0198   cout << "  N files found      : " << nfound << endl;
0199
0200   if (nfound < nfiles)
0201   {
0202     cout << "WARNING!! There are not enough files for the desired luminosity." << endl;
0203     cout << "          Will scale the statistical uncertainties accordingly." << endl;
0204     // return;
0205   }
0206
0207   // chain the desired files
0208   ttree = new TChain("ttree");
0209   int nread = 0;
0210
0211   TString fileList = gSystem->GetFromPipe(Form("ls %s/*.root", inDir));
0212
0213   int spos = fileList.First("\n");
0214   while (spos > 0 && nread < nfiles)
0215   {
0216     TString tfile = fileList(0, spos);
0217     // cout << tsub << endl;
0218
0219     ttree->Add(tfile.Data());
0220
0221
0222     // chop off the file
0223     fileList = fileList(spos + 1, fileList.Length());
0224     spos = fileList.First("\n");
0225     nread++;
0226   }
0227
0228   cout << " successfully read in " << nread << " files" << endl;
0229   lumiRead = lumiPerFile * nread;
0230
0231
0232   ttree->SetBranchAddress("event", &event);
0233   ttree->SetBranchAddress("truthjet_n", &truthjet_n);
0234   ttree->SetBranchAddress("truthjet_parton_flavor", truthjet_parton_flavor);
0235   ttree->SetBranchAddress("truthjet_hadron_flavor", truthjet_hadron_flavor);
0236   ttree->SetBranchAddress("truthjet_pt", truthjet_pt);
0237   ttree->SetBranchAddress("truthjet_eta", truthjet_eta);
0238   ttree->SetBranchAddress("truthjet_phi", truthjet_phi);
0239
0240   Long64_t nentries = ttree->GetEntries();
0241
0242
0243   //==========================================================================//
0244   // GET MOMENTUM IMBALANCE
0245   //==========================================================================//
0246   cout << endl;
0247   cout << "--> Running over " << nentries << endl;
0248
0249   int iprint = 0;
0250   for (Long64_t ientry = 0; ientry < nentries; ientry++)
0251   {
0252     ttree->GetEntry(ientry);
0253
0254     if (ientry % 100000 == 0) cout << "----> Processing event " << ientry << endl;
0255
0256     //-- apply acceptance to each jet
0257     for (int i = 0; i < truthjet_n; i++)
0258       acc[i] = (rand->Uniform() < bJetEff);
0259       // acc[i] = (rand->Uniform() < bJetEff) && (rand->Uniform() < RAA);
0260
0261     //-- get kinematics for all b-jets
0262     for (int i = 0; i < truthjet_n; i++)
0263     {
0264       if (TMath::Abs(truthjet_parton_flavor[i]) != 5)
0265         continue;
0266
0267       // single b-jet kinematics
0268       hjet_pT->Fill(truthjet_pt[i]);
0269       hjet_eta->Fill(truthjet_eta[i]);
0270       hjet_phi->Fill(truthjet_phi[i]);
0271
0272       for (int j = i + 1; j < truthjet_n; j++)
0273       {
0274         if (TMath::Abs(truthjet_parton_flavor[j]) != 5)
0275           continue;
0276
0277
0278         // find the leading and subleading jets
0279         int i1, i2;
0280         double pT1, pT2;
0281         double phi1, phi2;
0282         double eta1, eta2;
0283         bool acc1, acc2;
0284
0285         if (truthjet_pt[i] > truthjet_pt[j])
0286         {
0287           i1 = i;
0288           i2 = j;
0289         }
0290         else
0291         {
0292           i1 = j;
0293           i2 = i;
0294         }
0295
0296         pT1 = truthjet_pt[i1];
0297         phi1 = truthjet_phi[i1];
0298         eta1 = truthjet_eta[i1];
0299         acc1 = acc[i1];
0300
0301         pT2 = truthjet_pt[i2];
0302         phi2 = truthjet_phi[i2];
0303         eta2 = truthjet_eta[i2];
0304         acc2 = acc[i2];
0305
0306
0307
0308         if (iprint < 20)
0309         {
0310           cout << " ientry: " << ientry << endl
0311                << "        i1:" << i1 << " pT1:" << pT1 << " eta1:" << eta1 << " acc1:" << acc1 << endl
0312                << "        i2:" << i2 << " pT2:" << pT2 << " eta2:" << eta2 << " acc2:" << acc2 << endl;
0313           iprint++;
0314         }
0315
0316
0317
0318         // check if we accept the leading jet
0319         // only take RAA hit once
0320         if ( pT1 < pT1min || TMath::Abs(eta1) > etamax || !acc1  || rand->Uniform() > RAA)
0321           continue;
0322
0323         // calculate the delta phi
0324         double dphi = TMath::Abs(TMath::ATan2(TMath::Sin(phi1 - phi2), TMath::Cos(phi1 - phi2)));
0325         if (dphi > TMath::Pi())
0326           cout << " " << truthjet_phi[i] << " " << truthjet_phi[j] << " " << dphi << endl;
0327
0328         // calculate efficiency for finding the dijet
0329         hdijet_sing_pT1->Fill(pT1);
0330
0331         // cut on the subleading jet
0332         if ( pT2 < pT2min || TMath::Abs(eta2) > etamax || !acc2 )
0333           continue;
0334
0335         hdijet_sing_eta2->Fill(eta2);
0336
0337
0338         // fill efficiency for finding the dijet
0339         for (int k = 0; k < NETACUT; k++)
0340         {
0341           if ( TMath::Abs(eta2) < etacut[k])
0342             hdijet_eff_pT1[k]->Fill(pT1);
0343         }
0344         for (int k = 0; k < PTCUT; k++)
0345         {
0346           if (pT2 > ptcut[k])
0347             hdijet_eff_eta2[k]->Fill(eta2);
0348         }
0349
0350         hdijet_dphi->Fill(dphi);
0351
0352         if ( dphi < 2.*TMath::Pi() / 3)
0353           continue;
0354
0355
0356         // fill diagnostic histograms
0357
0358         hdijet_pT1->Fill(pT1);
0359         hdijet_pT2->Fill(pT2);
0360         hdijet_xj->Fill(pT2 / pT1);
0361
0362       } // j
0363     } // i
0364
0365   } // ientry
0366
0367   // first get some statistics
0368   cout << " N b dijets: " << hdijet_dphi->Integral() << endl;
0369   cout << " N b dijets w / dphi cut: " << hdijet_xj->Integral() << endl;
0370
0371   // calculate efficiencies
0372   for (int i = 0; i < NETACUT; i++)
0373   {
0374     hdijet_eff_pT1[i]->Divide(hdijet_sing_pT1);
0375     cout << " i: " << hdijet_eff_pT1[i]->GetMaximum() << endl;
0376   }
0377   for (int i = 0; i < PTCUT; i++)
0378   {
0379     hdijet_eff_eta2[i]->Divide(hdijet_sing_eta2);
0380   }
0381
0382   //-- normalize to integral 1
0383   hdijet_dphi->Scale(1. / hdijet_dphi->Integral());
0384   hdijet_xj->Scale(1. / hdijet_xj->Integral());
0385
0386   //-- scale statistical uncertainties if not enough simulated events
0387   if ( lumiRead < lumiDesired )
0388   {
0389     cout << endl;
0390     cout << "-- Scaling statistical uncertainties by:" << endl;
0391     cout << "  sqrt(" << lumiRead << "/" << lumiDesired << ") = "
0392          << TMath::Sqrt(lumiRead / lumiDesired) << endl;
0393     for (int ix = 1; ix <= hdijet_xj->GetNbinsX(); ix++)
0394     {
0395       double be = hdijet_xj->GetBinError(ix);
0396       be = be * TMath::Sqrt(lumiRead / lumiDesired);
0397       hdijet_xj->SetBinError(ix, be);
0398     } // ix
0399   }
0400
0401   //-- fix statistical uncertainties for purity
0402   hdijet_xj_fix = (TH1D*) hdijet_xj->Clone("hdijet_xj_fix");
0403   hdijet_xj_fix->SetLineColor(kRed);
0404   hdijet_xj_fix->SetMarkerColor(kRed);
0405   for (int ix = 1; ix <= hdijet_xj->GetNbinsX(); ix++)
0406   {
0407     double bc = hdijet_xj->GetBinContent(ix);
0408     double be = hdijet_xj->GetBinError(ix);
0409
0410     // increase the bin error due to the purity
0411     // need to square it, once for each bjet
0412     be = be / TMath::Sqrt(bJetPur * bJetPur);
0413     hdijet_xj_fix->SetBinError(ix, be);
0414   } // ix
0415
0416
0417
0418   //==========================================================================//
0419   // <x_j>
0420   //==========================================================================//
0421   cout << endl;
0422   cout << "--> Calculating <x_j>" << endl;
0423
0424   double sum = 0;
0425   double w = 0;
0426   double err = 0;
0427   for (int i = 1; i <= hdijet_xj->GetNbinsX(); i++)
0428   {
0429     double c = hdijet_xj->GetBinContent(i);
0430     if (c > 0)
0431     {
0432       sum += c * hdijet_xj->GetBinCenter(i);
0433       w += c;
0434       err += TMath::Power(c * hdijet_xj->GetBinError(i), 2);
0435     }
0436   }
0437   double mean = sum / w;
0438   err = TMath::Sqrt(err) / w;
0439
0440   cout << "  <x_j>=" << mean << " +/- " << err << endl;
0441
0442   //==========================================================================//
0443   // PLOT OBJECTS
0444   //==========================================================================//
0445   cout << endl;
0446   cout << "-- > Plotting" << endl;
0447
0448   hdijet_xj->SetMarkerStyle(kFullCircle);
0449   hdijet_xj->SetMarkerColor(kBlack);
0450   hdijet_xj->SetLineColor(kBlack);
0451   // hdijet_xj->GetYaxis()->SetTitleFont(63);
0452   // hdijet_xj->GetYaxis()->SetTitleSize(24);
0453   // hdijet_xj->GetYaxis()->SetTitleOffset(1.4);
0454   // hdijet_xj->GetYaxis()->SetLabelFont(63);
0455   // hdijet_xj->GetYaxis()->SetLabelSize(20);
0456   // hdijet_xj->GetXaxis()->SetTitleFont(63);
0457   // hdijet_xj->GetXaxis()->SetTitleSize(24);
0458   // hdijet_xj->GetXaxis()->SetTitleOffset(1.0);
0459   // hdijet_xj->GetXaxis()->SetLabelFont(63);
0460   // hdijet_xj->GetXaxis()->SetLabelSize(20);
0461
0462   hdijet_xj_fix->SetMarkerStyle(kFullCircle);
0463   hdijet_xj_fix->SetMarkerColor(kBlack);
0464   hdijet_xj_fix->SetLineColor(kBlack);
0465
0466   hdijet_dphi->SetMarkerStyle(kFullCircle);
0467   hdijet_dphi->SetMarkerColor(kBlack);
0468   hdijet_dphi->SetLineColor(kBlack);
0469   hdijet_dphi->GetYaxis()->SetTitleFont(63);
0470   hdijet_dphi->GetYaxis()->SetTitleSize(24);
0471   hdijet_dphi->GetYaxis()->SetTitleOffset(1.4);
0472   hdijet_dphi->GetYaxis()->SetLabelFont(63);
0473   hdijet_dphi->GetYaxis()->SetLabelSize(20);
0474   hdijet_dphi->GetXaxis()->SetTitleFont(63);
0475   hdijet_dphi->GetXaxis()->SetTitleSize(24);
0476   hdijet_dphi->GetXaxis()->SetTitleOffset(1.0);
0477   hdijet_dphi->GetXaxis()->SetLabelFont(63);
0478   hdijet_dphi->GetXaxis()->SetLabelSize(20);
0479
0480   hdijet_pT1->SetLineColor(kBlue);
0481   hdijet_pT2->SetLineColor(kRed);
0482
0483   TH1D* heff_axis_pt = new TH1D("heff_axis_pt",
0484                                 "; p_{T, 1} [GeV / c]; Fraction of Dijet / Single Accepted",
0485                                 20, 00, 100);
0486   heff_axis_pt->SetMinimum(0);
0487   heff_axis_pt->SetMaximum(1.);
0488   heff_axis_pt->GetYaxis()->SetTitleFont(63);
0489   heff_axis_pt->GetYaxis()->SetTitleSize(24);
0490   heff_axis_pt->GetYaxis()->SetTitleOffset(1.4);
0491   heff_axis_pt->GetYaxis()->SetLabelFont(63);
0492   heff_axis_pt->GetYaxis()->SetLabelSize(20);
0493   heff_axis_pt->GetXaxis()->SetTitleFont(63);
0494   heff_axis_pt->GetXaxis()->SetTitleSize(24);
0495   heff_axis_pt->GetXaxis()->SetTitleOffset(1.0);
0496   heff_axis_pt->GetXaxis()->SetLabelFont(63);
0497   heff_axis_pt->GetXaxis()->SetLabelSize(20);
0498
0499   TH1D* heff_axis_eta = new TH1D("heff_axis_eta",
0500                                  "; #eta_{2}; Fraction of Dijet / Single Accepted",
0501                                  150, -1.5, 1.5);
0502   heff_axis_eta->SetMinimum(0);
0503   heff_axis_eta->SetMaximum(1.);
0504   heff_axis_eta->GetYaxis()->SetTitleFont(63);
0505   heff_axis_eta->GetYaxis()->SetTitleSize(24);
0506   heff_axis_eta->GetYaxis()->SetTitleOffset(1.4);
0507   heff_axis_eta->GetYaxis()->SetLabelFont(63);
0508   heff_axis_eta->GetYaxis()->SetLabelSize(20);
0509   heff_axis_eta->GetXaxis()->SetTitleFont(63);
0510   heff_axis_eta->GetXaxis()->SetTitleSize(24);
0511   heff_axis_eta->GetXaxis()->SetTitleOffset(1.0);
0512   heff_axis_eta->GetXaxis()->SetLabelFont(63);
0513   heff_axis_eta->GetXaxis()->SetLabelSize(20);
0514
0515   for (int i = 0; i < NETACUT; i++)
0516     hdijet_eff_pT1[i]->SetLineColor(etaColor[i]);
0517
0518 //-- legends
0519   TLegend *leg_jetpt = new TLegend(0.6, 0.5, 0.95, 0.95);
0520   leg_jetpt->SetFillStyle(0);
0521   leg_jetpt->SetBorderSize(0);
0522   leg_jetpt->SetTextFont(63);
0523   leg_jetpt->SetTextSize(12);
0524   leg_jetpt->AddEntry(hjet_pT, "Inclusive b - jet", "L");
0525   leg_jetpt->AddEntry(hdijet_pT1, "leading b - jet", "L");
0526   leg_jetpt->AddEntry(hdijet_pT2, "subleading b - jet", "L");
0527
0528
0529 //-- other
0530   TLatex ltxt;
0531   ltxt.SetNDC();
0532   ltxt.SetTextFont(63);
0533   ltxt.SetTextSize(20);
0534
0535 //==========================================================================//
0536 // PLOT
0537 //==========================================================================//
0538
0539 // plot b-jet kinematics
0540   TCanvas *cjet = new TCanvas("cjet", "b - jet", 1200, 400);
0541   cjet->SetMargin(0, 0, 0, 0);
0542   cjet->Divide(3, 1);
0543
0544   cjet->GetPad(1)->SetMargin(0.12, 0.02, 0.12, 0.02);
0545   cjet->GetPad(2)->SetMargin(0.12, 0.02, 0.12, 0.02);
0546   cjet->GetPad(3)->SetMargin(0.12, 0.02, 0.12, 0.02);
0547
0548   cjet->cd(1);
0549   gPad->SetLogy();
0550   hjet_pT->GetYaxis()->SetRangeUser(0.5, 1.5 * hjet_pT->GetMaximum());
0551   hjet_pT->GetXaxis()->SetRangeUser(0, 50);
0552   hjet_pT->Draw();
0553   // hdijet_pT1->Draw("same");
0554   // hdijet_pT2->Draw("same");
0555
0556   // leg_jetpt->Draw("same");
0557
0558   cjet->cd(2);
0559   hjet_eta->Draw("HIST");
0560
0561   cjet->cd(3);
0562   hjet_phi->Draw("HIST");
0563
0564
0565 // plot dijet eff
0566   TCanvas *ceff = new TCanvas("ceff", "eff", 1200, 600);
0567   ceff->Divide(2, 1);
0568   ceff->GetPad(1)->SetMargin(0.12, 0.02, 0.12, 0.02);
0569   ceff->GetPad(1)->SetTicks(1, 1);
0570   ceff->GetPad(2)->SetMargin(0.12, 0.02, 0.12, 0.02);
0571   ceff->GetPad(2)->SetTicks(1, 1);
0572
0573   ceff->cd(1);
0574
0575   heff_axis_pt->Draw();
0576
0577   for (int i = 0; i < NETACUT; i++)
0578     hdijet_eff_pT1[i]->Draw("L same");
0579
0580   ceff->cd(2);
0581
0582   heff_axis_eta->Draw();
0583
0584   for (int i = 0; i < NETACUT; i++)
0585     hdijet_eff_eta2[i]->Draw("L same");
0586
0587 // plot dijet checks
0588   TCanvas *cdijet2 = new TCanvas("cdijet2", "dijet", 1200, 600);
0589   cdijet2->SetMargin(0, 0, 0, 0);
0590   cdijet2->Divide(2, 1);
0591
0592   cdijet2->GetPad(1)->SetMargin(0.12, 0.02, 0.12, 0.02);
0593   cdijet2->GetPad(2)->SetMargin(0.12, 0.02, 0.12, 0.02);
0594
0595   cdijet2->cd(1);
0596   hdijet_xj->Draw();
0597
0598   ltxt.DrawLatex(0.2, 0.92, "Di b-jets (Pythia8)");
0599   ltxt.DrawLatex(0.2, 0.86, "p + p #sqrt{s_{_{NN}}}=200 GeV");
0600   ltxt.DrawLatex(0.2, 0.80, "anti - k_ {T}, R = 0.4");
0601   ltxt.DrawLatex(0.2, 0.74, Form("p_{T, 1} > % .0f GeV", pT1min));
0602   ltxt.DrawLatex(0.2, 0.68, Form("p_{T, 2} > % .0f GeV", pT2min));
0603   ltxt.DrawLatex(0.2, 0.62, Form(" | #Delta#phi_{12}| > 2#pi/3"));
0604
0605   cdijet2->cd(2);
0606   hdijet_dphi->Draw();
0607
0608
0609
0610
0611
0612
0613   // make this one consistent with sPHENIX style
0614   TCanvas *cdijet = new TCanvas("cdijet", "dijet", 600, 600);
0615   // cdijet->SetMargin(0.12, 0.02, 0.12, 0.02);
0616   // cdijet->SetTicks(1, 1);
0617   cdijet->cd();
0618   hdijet_xj_fix->GetYaxis()->SetRangeUser(0, 0.3);
0619   hdijet_xj_fix->Draw();
0620   // hdijet_xj->Draw("same");
0621
0622
0623   TLegend *legsphenix = new TLegend(0.15, 0.5, 0.5, 0.9);
0624   legsphenix->SetFillStyle(0);
0625   legsphenix->SetBorderSize(0);
0626   legsphenix->AddEntry("", "#it{#bf{sPHENIX}} Simulation", "");
0627   legsphenix->AddEntry("", "Di b-jets (Pythia8, CTEQ6L)", "");
0628   // if (is_pp)
0629   // {
0630   //   legsphenix->AddEntry("", "p + p #sqrt{s_{_{NN}}} = 200 GeV", "");
0631   //   legsphenix->AddEntry("", Form("#int L dt = %.0f pb^{-1} |z| < 10 cm", lumiDesired), "");
0632   // }
0633   // else
0634   // {
0635   //   legsphenix->AddEntry("", "0-10% Au + Au #sqrt{s_{_{NN}}}=200 GeV", "");
0636   //   legsphenix->AddEntry("", "10B events |z| < 10 cm", "");
0637   // }
0638   legsphenix->AddEntry("", scol.Data(), "");
0639   legsphenix->AddEntry("", slumi.Data(), "");
0640   legsphenix->AddEntry("", "anti-k_{T}, R = 0.4", "");
0641   legsphenix->AddEntry("", Form(" |#eta| < %.1f", etamax), "");
0642   legsphenix->AddEntry("", Form(" |#Delta#phi_{12}| > 2#pi/3"), "");
0643   legsphenix->AddEntry("", Form("p_{T, 1} > % .0f GeV", pT1min), "");
0644   legsphenix->AddEntry("", Form("p_{T, 2} > % .0f GeV", pT2min), "");
0645   legsphenix->AddEntry("", Form("%.0f%% b-jet Eff, %.0f%% b-jet Pur.", bJetEff * 100., bJetPur * 100.), "");
0646   legsphenix->Draw("same");
0647
0648   // ltxt.DrawLatex(0.2, 0.92, "Di b-jets (Pythia8)");
0649   // if (is_pp)
0650   // {
0651   //   ltxt.DrawLatex(0.2, 0.86, "p + p #sqrt{s_{_{NN}}} = 200 GeV");
0652   //   ltxt.DrawLatex(0.2, 0.80, "#int L dt = 175 pb^{-1} |z| < 10 cm");
0653   // }
0654   // else
0655   // {
0656   //   ltxt.DrawLatex(0.2, 0.86, "0-10% Au + Au #sqrt{s_{_{NN}}}=200 GeV");
0657   //   ltxt.DrawLatex(0.2, 0.80, "10B events |z| < 10 cm");
0658   // }
0659   // ltxt.DrawLatex(0.2, 0.74, "anti-k_{T}, R = 0.4");
0660   // ltxt.DrawLatex(0.2, 0.68, Form("p_{T, 1} > % .0f GeV", pT1min));
0661   // ltxt.DrawLatex(0.2, 0.62, Form("p_{T, 2} > % .0f GeV", pT2min));
0662   // ltxt.DrawLatex(0.2, 0.56, Form(" |#eta| < %.0f", etamax));
0663   // ltxt.DrawLatex(0.2, 0.50, Form(" |#Delta#phi_{12}| > 2#pi/3"));
0664
0665
0666   //==========================================================================//
0667   // PRINT PLOTS
0668   //==========================================================================//
0669   if (print_plots)
0670   {
0671     if (is_pp)
0672     {
0673       // cjet->Print("bjet_kinematics_pp.pdf");
0674       // cdijet2->Print("dibjet_2panel_pp.pdf");
0675       cdijet->Print("dibjet_imbalance_pp.pdf");
0676       cdijet->Print("dibjet_imbalance_pp.C");
0677       cdijet->Print("dibjet_imbalance_pp.root");
0678
0679       TFile *fout = new TFile("dibjet_imbalance_histos_pp.root","RECREATE");
0680       fout->cd();
0681       hdijet_xj->Write();
0682       hdijet_xj_fix->Write();
0683       fout->Close();
0684       delete fout;
0685
0686     }
0687     else
0688     {
0689       cdijet->Print("dibjet_imbalance_AuAu0-10.pdf");
0690       cdijet->Print("dibjet_imbalance_AuAu0-10.C");
0691       cdijet->Print("dibjet_imbalance_AuAu0-10.root");
0692
0693       TFile *fout = new TFile("dibjet_imbalance_histos_AuAu0-10.root","RECREATE");
0694       fout->cd();
0695       hdijet_xj->Write();
0696       hdijet_xj_fix->Write();
0697       fout->Close();
0698       delete fout;
0699     }
0700   }
0701
0702 }