Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Draw_BDiJetImbalance.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Draw_BDiJetImbalance.C
1 
2 //
3 // Draw the Di b-jet momentum imbalance
4 //
5 // Useful references:
6 // CMS di b-jet: http://cds.cern.ch/record/2202805/files/HIN-16-005-pas.pdf
7 //
9 //
10 // Darren McGlinchey
11 // 6 Jan 2017
12 //
14 
15 #include <TFile.h>
16 #include <TChain.h>
17 #include <TTree.h>
18 #include <TH1.h>
19 #include <TCanvas.h>
20 #include <TLatex.h>
21 #include <TStyle.h>
22 #include <TMath.h>
23 #include <TLegend.h>
24 #include <TRandom3.h>
25 #include <TString.h>
26 #include <TSystem.h>
27 
28 #include <iostream>
29 
30 using namespace std;
31 
32 void Draw_BDiJetImbalance()
33 {
34  gStyle->SetOptStat(0);
35  gStyle->SetOptTitle(0);
36 
37  //==========================================================================//
38  // SET RUNNING CONDITIONS
39  //==========================================================================//
40  bool print_plots = true;
41 
42  bool is_pp = false;
43 
44  // const char* inFile = "HFtag_bjet.root";
45  // const char* inFile = "HFtag_bjet_pp.root";
46  // const char* inFile = "/phenix/plhf/dcm07e/sPHENIX/bjetsims/test.root";
47  // if (!is_pp)
48  // const char* inFile = "HFtag_bjet_AuAu0-10.root";
49 
50  const char* inDir = "/phenix/plhf/dcm07e/sPHENIX/bjetsims/ana";
51 
52  // // desired nn integrated luminosity [pb^{-1}]
53  // double lumiDesired = 175;
54  // if ( !is_pp )
55  // lumiDesired = 229; // AuAu 0-10% central (10B events)
56  // desired nn integrated luminosity [pb^{-1}]
57  double lumiDesired = 200;
58  TString scol("p + p #sqrt{s_{_{NN}}} = 200 GeV");
59  TString slumi(Form("#int L dt = %.0f pb^{-1} |z| < 10 cm", lumiDesired));
60  if ( !is_pp )
61  {
62  // equivelant nn luminosity:
63  // N_{evt}^{NN} / sigma_NN = N_{evt}^{AA}*Ncoll / sigma_NN
64  // = 10e9 * 962 / 42e9 pb = 229 pb^-1
65  double evntDesired = 24e9; // AuAu 0-10% central
66  lumiDesired = evntDesired * 962 / 42e9; // AuAu 0-10% central
67  scol = "0-10% Au + Au #sqrt{s_{_{NN}}}=200 GeV";
68  slumi = Form("%.0fB events |z| < 10 cm", evntDesired / 1.e9);
69  }
70 
71  // luminosity per file generated [pb^{-1}]
72  double lumiPerFile = 99405 / (4.641e-06 * 1e12);
73  double lumiRead = 0;
74 
75  double pT1min = 20;
76  double pT2min = 10;
77  // double etamax = 1.0;
78  double etamax = 0.7;
79 
80  // double bJetEff = 0.5; // b-jet tagging efficiency
81  double bJetEff = 0.60; // p+p b-jet tagging efficiency
82  if ( !is_pp )
83  bJetEff = 0.40;
84 
85  double bJetPur = 0.40; // b-jet tagging purity
86 
87  double RAA = 1.0;
88  if (!is_pp)
89  RAA = 0.6;
90 
91  const int NETACUT = 3;
92  double etacut[] = {1.3, 1.0, 0.7};
93  double etaColor[] = {kBlue, kRed, kBlack};
94 
95  const int PTCUT = 3;
96  double ptcut[] = {10, 20, 40};
97  double ptColor[] = {kBlue, kRed, kBlack};
98 
99  //==========================================================================//
100  // DECLARE VARIABLES
101  //==========================================================================//
102 
103  //-- tree variables
104  TChain *ttree;
105  Int_t event;
106  Int_t truthjet_n;
107  Int_t truthjet_parton_flavor[100];
108  Int_t truthjet_hadron_flavor[100];
109  Float_t truthjet_pt[100];
110  Float_t truthjet_eta[100];
111  Float_t truthjet_phi[100];
112 
113  //-- histograms
114 
115  TH1D* hjet_pT = new TH1D("hjet_pT", ";p_{T}^{jet} [GeV/c]", 20, 0, 100);
116  TH1D* hjet_eta = new TH1D("hjet_eta", ";#eta^{jet}", 40, -2., 2.);
117  TH1D* hjet_phi = new TH1D("hjet_phi", ";#phi^{jet}", 40, -4, 4);
118 
119  TH1D* hdijet_pT1 = new TH1D("hdijet_pT1", ";p_{T,1}^{jet} [GeV/c]", 20, 0, 100);
120  TH1D* hdijet_pT2 = new TH1D("hdijet_pT2", ";p_{T,2}^{jet} [GeV/c]", 20, 0, 100);
121  TH1D* hdijet_xj = new TH1D("hdijet_xj", ";x_{j} = p_{T,2} / p_{T,1}; event fraction", 10, 0, 1);
122  TH1D* hdijet_dphi = new TH1D("hdijet_dphi", ";|#Delta#phi_{12}|; event fraction", 10, 0, TMath::Pi());
123 
124  TH1D* hdijet_sing_pT1 = new TH1D("hdijet_sing_pT1",
125  ";p_{T,1} [GeV/c];Fraction of Dijet / Single Accepted",
126  20, 00, 100);
127  TH1D* hdijet_sing_eta2 = new TH1D("hdijet_sing_eta2",
128  ";#eta_{2};Fraction of Dijet / Single Accepted",
129  15, -1.5, 1.5);
130  TH1D* hdijet_eff_pT1[3];
131  TH1D* hdijet_eff_eta2[3];
132  for (int i = 0; i < NETACUT; i++)
133  {
134  hdijet_eff_pT1[i] = new TH1D(Form("hdijet_eff_pT1_%i", i),
135  ";p_{T,1} [GeV/c];Fraction of Dijet / Single Accepted",
136  20, 00, 100);
137  hdijet_eff_pT1[i]->SetLineColor(etaColor[i]);
138  hdijet_eff_pT1[i]->SetLineWidth(2);
139  } // i
140  for (int i = 0; i < PTCUT; i++)
141  {
142  hdijet_eff_eta2[i] = new TH1D(Form("hdijet_eff_eta2_%i", i),
143  ";#eta;Fraction of Dijet / Single Accepted",
144  15, -1.5, 1.5);
145  hdijet_eff_eta2[i]->SetLineColor(etaColor[i]);
146  hdijet_eff_eta2[i]->SetLineWidth(2);
147  } // i
148 
149  hjet_pT->Sumw2();
150  hjet_eta->Sumw2();
151  hjet_phi->Sumw2();
152 
153  hdijet_pT1->Sumw2();
154  hdijet_pT2->Sumw2();
155  hdijet_dphi->Sumw2();
156  hdijet_xj->Sumw2();
157 
158  // for fixing the uncertainties
159  TH1D* hdijet_xj_fix;
160 
161  //-- other
162  TRandom3 *rand = new TRandom3();
163 
164  bool acc[100];
165 
166  //==========================================================================//
167  // LOAD TTREE
168  //==========================================================================//
169  // cout << endl;
170  // cout << "--> Reading data from " << inFile << endl;
171 
172  // TFile *fin = TFile::Open(inFile);
173  // if (!fin)
174  // {
175  // cout << "ERROR!! Unable to open " << inFile << endl;
176  // return;
177  // }
178 
179  // ttree = (TTree*) fin->Get("ttree");
180  // if (!ttree)
181  // {
182  // cout << "ERROR!! Unable to find ttree in " << inFile << endl;
183  // return;
184  // }
185 
186  cout << endl;
187  cout << "--> Reading data from dir: " << inDir << endl;
188 
189  // calculate the number of files necessary for the desired luminosity
190  int nfiles = (int)(lumiDesired / lumiPerFile);
191  cout << " desired luminosity : " << lumiDesired << endl;
192  cout << " luminosity per file: " << lumiPerFile << endl;
193  cout << " N files needed : " << nfiles << endl;
194 
195 
196  // check the number of files found in the directory
197  int nfound = (gSystem->GetFromPipe(Form("ls %s/*.root | wc -l", inDir))).Atoi();
198  cout << " N files found : " << nfound << endl;
199 
200  if (nfound < nfiles)
201  {
202  cout << "WARNING!! There are not enough files for the desired luminosity." << endl;
203  cout << " Will scale the statistical uncertainties accordingly." << endl;
204  // return;
205  }
206 
207  // chain the desired files
208  ttree = new TChain("ttree");
209  int nread = 0;
210 
211  TString fileList = gSystem->GetFromPipe(Form("ls %s/*.root", inDir));
212 
213  int spos = fileList.First("\n");
214  while (spos > 0 && nread < nfiles)
215  {
216  TString tfile = fileList(0, spos);
217  // cout << tsub << endl;
218 
219  ttree->Add(tfile.Data());
220 
221 
222  // chop off the file
223  fileList = fileList(spos + 1, fileList.Length());
224  spos = fileList.First("\n");
225  nread++;
226  }
227 
228  cout << " successfully read in " << nread << " files" << endl;
229  lumiRead = lumiPerFile * nread;
230 
231 
232  ttree->SetBranchAddress("event", &event);
233  ttree->SetBranchAddress("truthjet_n", &truthjet_n);
234  ttree->SetBranchAddress("truthjet_parton_flavor", truthjet_parton_flavor);
235  ttree->SetBranchAddress("truthjet_hadron_flavor", truthjet_hadron_flavor);
236  ttree->SetBranchAddress("truthjet_pt", truthjet_pt);
237  ttree->SetBranchAddress("truthjet_eta", truthjet_eta);
238  ttree->SetBranchAddress("truthjet_phi", truthjet_phi);
239 
240  Long64_t nentries = ttree->GetEntries();
241 
242 
243  //==========================================================================//
244  // GET MOMENTUM IMBALANCE
245  //==========================================================================//
246  cout << endl;
247  cout << "--> Running over " << nentries << endl;
248 
249  int iprint = 0;
250  for (Long64_t ientry = 0; ientry < nentries; ientry++)
251  {
252  ttree->GetEntry(ientry);
253 
254  if (ientry % 100000 == 0) cout << "----> Processing event " << ientry << endl;
255 
256  //-- apply acceptance to each jet
257  for (int i = 0; i < truthjet_n; i++)
258  acc[i] = (rand->Uniform() < bJetEff);
259  // acc[i] = (rand->Uniform() < bJetEff) && (rand->Uniform() < RAA);
260 
261  //-- get kinematics for all b-jets
262  for (int i = 0; i < truthjet_n; i++)
263  {
264  if (TMath::Abs(truthjet_parton_flavor[i]) != 5)
265  continue;
266 
267  // single b-jet kinematics
268  hjet_pT->Fill(truthjet_pt[i]);
269  hjet_eta->Fill(truthjet_eta[i]);
270  hjet_phi->Fill(truthjet_phi[i]);
271 
272  for (int j = i + 1; j < truthjet_n; j++)
273  {
274  if (TMath::Abs(truthjet_parton_flavor[j]) != 5)
275  continue;
276 
277 
278  // find the leading and subleading jets
279  int i1, i2;
280  double pT1, pT2;
281  double phi1, phi2;
282  double eta1, eta2;
283  bool acc1, acc2;
284 
285  if (truthjet_pt[i] > truthjet_pt[j])
286  {
287  i1 = i;
288  i2 = j;
289  }
290  else
291  {
292  i1 = j;
293  i2 = i;
294  }
295 
296  pT1 = truthjet_pt[i1];
297  phi1 = truthjet_phi[i1];
298  eta1 = truthjet_eta[i1];
299  acc1 = acc[i1];
300 
301  pT2 = truthjet_pt[i2];
302  phi2 = truthjet_phi[i2];
303  eta2 = truthjet_eta[i2];
304  acc2 = acc[i2];
305 
306 
307 
308  if (iprint < 20)
309  {
310  cout << " ientry: " << ientry << endl
311  << " i1:" << i1 << " pT1:" << pT1 << " eta1:" << eta1 << " acc1:" << acc1 << endl
312  << " i2:" << i2 << " pT2:" << pT2 << " eta2:" << eta2 << " acc2:" << acc2 << endl;
313  iprint++;
314  }
315 
316 
317 
318  // check if we accept the leading jet
319  // only take RAA hit once
320  if ( pT1 < pT1min || TMath::Abs(eta1) > etamax || !acc1 || rand->Uniform() > RAA)
321  continue;
322 
323  // calculate the delta phi
324  double dphi = TMath::Abs(TMath::ATan2(TMath::Sin(phi1 - phi2), TMath::Cos(phi1 - phi2)));
325  if (dphi > TMath::Pi())
326  cout << " " << truthjet_phi[i] << " " << truthjet_phi[j] << " " << dphi << endl;
327 
328  // calculate efficiency for finding the dijet
329  hdijet_sing_pT1->Fill(pT1);
330 
331  // cut on the subleading jet
332  if ( pT2 < pT2min || TMath::Abs(eta2) > etamax || !acc2 )
333  continue;
334 
335  hdijet_sing_eta2->Fill(eta2);
336 
337 
338  // fill efficiency for finding the dijet
339  for (int k = 0; k < NETACUT; k++)
340  {
341  if ( TMath::Abs(eta2) < etacut[k])
342  hdijet_eff_pT1[k]->Fill(pT1);
343  }
344  for (int k = 0; k < PTCUT; k++)
345  {
346  if (pT2 > ptcut[k])
347  hdijet_eff_eta2[k]->Fill(eta2);
348  }
349 
350  hdijet_dphi->Fill(dphi);
351 
352  if ( dphi < 2.*TMath::Pi() / 3)
353  continue;
354 
355 
356  // fill diagnostic histograms
357 
358  hdijet_pT1->Fill(pT1);
359  hdijet_pT2->Fill(pT2);
360  hdijet_xj->Fill(pT2 / pT1);
361 
362  } // j
363  } // i
364 
365  } // ientry
366 
367  // first get some statistics
368  cout << " N b dijets: " << hdijet_dphi->Integral() << endl;
369  cout << " N b dijets w / dphi cut: " << hdijet_xj->Integral() << endl;
370 
371  // calculate efficiencies
372  for (int i = 0; i < NETACUT; i++)
373  {
374  hdijet_eff_pT1[i]->Divide(hdijet_sing_pT1);
375  cout << " i: " << hdijet_eff_pT1[i]->GetMaximum() << endl;
376  }
377  for (int i = 0; i < PTCUT; i++)
378  {
379  hdijet_eff_eta2[i]->Divide(hdijet_sing_eta2);
380  }
381 
382  //-- normalize to integral 1
383  hdijet_dphi->Scale(1. / hdijet_dphi->Integral());
384  hdijet_xj->Scale(1. / hdijet_xj->Integral());
385 
386  //-- scale statistical uncertainties if not enough simulated events
387  if ( lumiRead < lumiDesired )
388  {
389  cout << endl;
390  cout << "-- Scaling statistical uncertainties by:" << endl;
391  cout << " sqrt(" << lumiRead << "/" << lumiDesired << ") = "
392  << TMath::Sqrt(lumiRead / lumiDesired) << endl;
393  for (int ix = 1; ix <= hdijet_xj->GetNbinsX(); ix++)
394  {
395  double be = hdijet_xj->GetBinError(ix);
396  be = be * TMath::Sqrt(lumiRead / lumiDesired);
397  hdijet_xj->SetBinError(ix, be);
398  } // ix
399  }
400 
401  //-- fix statistical uncertainties for purity
402  hdijet_xj_fix = (TH1D*) hdijet_xj->Clone("hdijet_xj_fix");
403  hdijet_xj_fix->SetLineColor(kRed);
404  hdijet_xj_fix->SetMarkerColor(kRed);
405  for (int ix = 1; ix <= hdijet_xj->GetNbinsX(); ix++)
406  {
407  double bc = hdijet_xj->GetBinContent(ix);
408  double be = hdijet_xj->GetBinError(ix);
409 
410  // increase the bin error due to the purity
411  // need to square it, once for each bjet
412  be = be / TMath::Sqrt(bJetPur * bJetPur);
413  hdijet_xj_fix->SetBinError(ix, be);
414  } // ix
415 
416 
417 
418  //==========================================================================//
419  // <x_j>
420  //==========================================================================//
421  cout << endl;
422  cout << "--> Calculating <x_j>" << endl;
423 
424  double sum = 0;
425  double w = 0;
426  double err = 0;
427  for (int i = 1; i <= hdijet_xj->GetNbinsX(); i++)
428  {
429  double c = hdijet_xj->GetBinContent(i);
430  if (c > 0)
431  {
432  sum += c * hdijet_xj->GetBinCenter(i);
433  w += c;
434  err += TMath::Power(c * hdijet_xj->GetBinError(i), 2);
435  }
436  }
437  double mean = sum / w;
438  err = TMath::Sqrt(err) / w;
439 
440  cout << " <x_j>=" << mean << " +/- " << err << endl;
441 
442  //==========================================================================//
443  // PLOT OBJECTS
444  //==========================================================================//
445  cout << endl;
446  cout << "-- > Plotting" << endl;
447 
448  hdijet_xj->SetMarkerStyle(kFullCircle);
449  hdijet_xj->SetMarkerColor(kBlack);
450  hdijet_xj->SetLineColor(kBlack);
451  // hdijet_xj->GetYaxis()->SetTitleFont(63);
452  // hdijet_xj->GetYaxis()->SetTitleSize(24);
453  // hdijet_xj->GetYaxis()->SetTitleOffset(1.4);
454  // hdijet_xj->GetYaxis()->SetLabelFont(63);
455  // hdijet_xj->GetYaxis()->SetLabelSize(20);
456  // hdijet_xj->GetXaxis()->SetTitleFont(63);
457  // hdijet_xj->GetXaxis()->SetTitleSize(24);
458  // hdijet_xj->GetXaxis()->SetTitleOffset(1.0);
459  // hdijet_xj->GetXaxis()->SetLabelFont(63);
460  // hdijet_xj->GetXaxis()->SetLabelSize(20);
461 
462  hdijet_xj_fix->SetMarkerStyle(kFullCircle);
463  hdijet_xj_fix->SetMarkerColor(kBlack);
464  hdijet_xj_fix->SetLineColor(kBlack);
465 
466  hdijet_dphi->SetMarkerStyle(kFullCircle);
467  hdijet_dphi->SetMarkerColor(kBlack);
468  hdijet_dphi->SetLineColor(kBlack);
469  hdijet_dphi->GetYaxis()->SetTitleFont(63);
470  hdijet_dphi->GetYaxis()->SetTitleSize(24);
471  hdijet_dphi->GetYaxis()->SetTitleOffset(1.4);
472  hdijet_dphi->GetYaxis()->SetLabelFont(63);
473  hdijet_dphi->GetYaxis()->SetLabelSize(20);
474  hdijet_dphi->GetXaxis()->SetTitleFont(63);
475  hdijet_dphi->GetXaxis()->SetTitleSize(24);
476  hdijet_dphi->GetXaxis()->SetTitleOffset(1.0);
477  hdijet_dphi->GetXaxis()->SetLabelFont(63);
478  hdijet_dphi->GetXaxis()->SetLabelSize(20);
479 
480  hdijet_pT1->SetLineColor(kBlue);
481  hdijet_pT2->SetLineColor(kRed);
482 
483  TH1D* heff_axis_pt = new TH1D("heff_axis_pt",
484  "; p_{T, 1} [GeV / c]; Fraction of Dijet / Single Accepted",
485  20, 00, 100);
486  heff_axis_pt->SetMinimum(0);
487  heff_axis_pt->SetMaximum(1.);
488  heff_axis_pt->GetYaxis()->SetTitleFont(63);
489  heff_axis_pt->GetYaxis()->SetTitleSize(24);
490  heff_axis_pt->GetYaxis()->SetTitleOffset(1.4);
491  heff_axis_pt->GetYaxis()->SetLabelFont(63);
492  heff_axis_pt->GetYaxis()->SetLabelSize(20);
493  heff_axis_pt->GetXaxis()->SetTitleFont(63);
494  heff_axis_pt->GetXaxis()->SetTitleSize(24);
495  heff_axis_pt->GetXaxis()->SetTitleOffset(1.0);
496  heff_axis_pt->GetXaxis()->SetLabelFont(63);
497  heff_axis_pt->GetXaxis()->SetLabelSize(20);
498 
499  TH1D* heff_axis_eta = new TH1D("heff_axis_eta",
500  "; #eta_{2}; Fraction of Dijet / Single Accepted",
501  150, -1.5, 1.5);
502  heff_axis_eta->SetMinimum(0);
503  heff_axis_eta->SetMaximum(1.);
504  heff_axis_eta->GetYaxis()->SetTitleFont(63);
505  heff_axis_eta->GetYaxis()->SetTitleSize(24);
506  heff_axis_eta->GetYaxis()->SetTitleOffset(1.4);
507  heff_axis_eta->GetYaxis()->SetLabelFont(63);
508  heff_axis_eta->GetYaxis()->SetLabelSize(20);
509  heff_axis_eta->GetXaxis()->SetTitleFont(63);
510  heff_axis_eta->GetXaxis()->SetTitleSize(24);
511  heff_axis_eta->GetXaxis()->SetTitleOffset(1.0);
512  heff_axis_eta->GetXaxis()->SetLabelFont(63);
513  heff_axis_eta->GetXaxis()->SetLabelSize(20);
514 
515  for (int i = 0; i < NETACUT; i++)
516  hdijet_eff_pT1[i]->SetLineColor(etaColor[i]);
517 
518 //-- legends
519  TLegend *leg_jetpt = new TLegend(0.6, 0.5, 0.95, 0.95);
520  leg_jetpt->SetFillStyle(0);
521  leg_jetpt->SetBorderSize(0);
522  leg_jetpt->SetTextFont(63);
523  leg_jetpt->SetTextSize(12);
524  leg_jetpt->AddEntry(hjet_pT, "Inclusive b - jet", "L");
525  leg_jetpt->AddEntry(hdijet_pT1, "leading b - jet", "L");
526  leg_jetpt->AddEntry(hdijet_pT2, "subleading b - jet", "L");
527 
528 
529 //-- other
530  TLatex ltxt;
531  ltxt.SetNDC();
532  ltxt.SetTextFont(63);
533  ltxt.SetTextSize(20);
534 
535 //==========================================================================//
536 // PLOT
537 //==========================================================================//
538 
539 // plot b-jet kinematics
540  TCanvas *cjet = new TCanvas("cjet", "b - jet", 1200, 400);
541  cjet->SetMargin(0, 0, 0, 0);
542  cjet->Divide(3, 1);
543 
544  cjet->GetPad(1)->SetMargin(0.12, 0.02, 0.12, 0.02);
545  cjet->GetPad(2)->SetMargin(0.12, 0.02, 0.12, 0.02);
546  cjet->GetPad(3)->SetMargin(0.12, 0.02, 0.12, 0.02);
547 
548  cjet->cd(1);
549  gPad->SetLogy();
550  hjet_pT->GetYaxis()->SetRangeUser(0.5, 1.5 * hjet_pT->GetMaximum());
551  hjet_pT->GetXaxis()->SetRangeUser(0, 50);
552  hjet_pT->Draw();
553  // hdijet_pT1->Draw("same");
554  // hdijet_pT2->Draw("same");
555 
556  // leg_jetpt->Draw("same");
557 
558  cjet->cd(2);
559  hjet_eta->Draw("HIST");
560 
561  cjet->cd(3);
562  hjet_phi->Draw("HIST");
563 
564 
565 // plot dijet eff
566  TCanvas *ceff = new TCanvas("ceff", "eff", 1200, 600);
567  ceff->Divide(2, 1);
568  ceff->GetPad(1)->SetMargin(0.12, 0.02, 0.12, 0.02);
569  ceff->GetPad(1)->SetTicks(1, 1);
570  ceff->GetPad(2)->SetMargin(0.12, 0.02, 0.12, 0.02);
571  ceff->GetPad(2)->SetTicks(1, 1);
572 
573  ceff->cd(1);
574 
575  heff_axis_pt->Draw();
576 
577  for (int i = 0; i < NETACUT; i++)
578  hdijet_eff_pT1[i]->Draw("L same");
579 
580  ceff->cd(2);
581 
582  heff_axis_eta->Draw();
583 
584  for (int i = 0; i < NETACUT; i++)
585  hdijet_eff_eta2[i]->Draw("L same");
586 
587 // plot dijet checks
588  TCanvas *cdijet2 = new TCanvas("cdijet2", "dijet", 1200, 600);
589  cdijet2->SetMargin(0, 0, 0, 0);
590  cdijet2->Divide(2, 1);
591 
592  cdijet2->GetPad(1)->SetMargin(0.12, 0.02, 0.12, 0.02);
593  cdijet2->GetPad(2)->SetMargin(0.12, 0.02, 0.12, 0.02);
594 
595  cdijet2->cd(1);
596  hdijet_xj->Draw();
597 
598  ltxt.DrawLatex(0.2, 0.92, "Di b-jets (Pythia8)");
599  ltxt.DrawLatex(0.2, 0.86, "p + p #sqrt{s_{_{NN}}}=200 GeV");
600  ltxt.DrawLatex(0.2, 0.80, "anti - k_ {T}, R = 0.4");
601  ltxt.DrawLatex(0.2, 0.74, Form("p_{T, 1} > % .0f GeV", pT1min));
602  ltxt.DrawLatex(0.2, 0.68, Form("p_{T, 2} > % .0f GeV", pT2min));
603  ltxt.DrawLatex(0.2, 0.62, Form(" | #Delta#phi_{12}| > 2#pi/3"));
604 
605  cdijet2->cd(2);
606  hdijet_dphi->Draw();
607 
608 
609 
610 
611 
612 
613  // make this one consistent with sPHENIX style
614  TCanvas *cdijet = new TCanvas("cdijet", "dijet", 600, 600);
615  // cdijet->SetMargin(0.12, 0.02, 0.12, 0.02);
616  // cdijet->SetTicks(1, 1);
617  cdijet->cd();
618  hdijet_xj_fix->GetYaxis()->SetRangeUser(0, 0.3);
619  hdijet_xj_fix->Draw();
620  // hdijet_xj->Draw("same");
621 
622 
623  TLegend *legsphenix = new TLegend(0.15, 0.5, 0.5, 0.9);
624  legsphenix->SetFillStyle(0);
625  legsphenix->SetBorderSize(0);
626  legsphenix->AddEntry("", "#it{#bf{sPHENIX}} Simulation", "");
627  legsphenix->AddEntry("", "Di b-jets (Pythia8, CTEQ6L)", "");
628  // if (is_pp)
629  // {
630  // legsphenix->AddEntry("", "p + p #sqrt{s_{_{NN}}} = 200 GeV", "");
631  // legsphenix->AddEntry("", Form("#int L dt = %.0f pb^{-1} |z| < 10 cm", lumiDesired), "");
632  // }
633  // else
634  // {
635  // legsphenix->AddEntry("", "0-10% Au + Au #sqrt{s_{_{NN}}}=200 GeV", "");
636  // legsphenix->AddEntry("", "10B events |z| < 10 cm", "");
637  // }
638  legsphenix->AddEntry("", scol.Data(), "");
639  legsphenix->AddEntry("", slumi.Data(), "");
640  legsphenix->AddEntry("", "anti-k_{T}, R = 0.4", "");
641  legsphenix->AddEntry("", Form(" |#eta| < %.1f", etamax), "");
642  legsphenix->AddEntry("", Form(" |#Delta#phi_{12}| > 2#pi/3"), "");
643  legsphenix->AddEntry("", Form("p_{T, 1} > % .0f GeV", pT1min), "");
644  legsphenix->AddEntry("", Form("p_{T, 2} > % .0f GeV", pT2min), "");
645  legsphenix->AddEntry("", Form("%.0f%% b-jet Eff, %.0f%% b-jet Pur.", bJetEff * 100., bJetPur * 100.), "");
646  legsphenix->Draw("same");
647 
648  // ltxt.DrawLatex(0.2, 0.92, "Di b-jets (Pythia8)");
649  // if (is_pp)
650  // {
651  // ltxt.DrawLatex(0.2, 0.86, "p + p #sqrt{s_{_{NN}}} = 200 GeV");
652  // ltxt.DrawLatex(0.2, 0.80, "#int L dt = 175 pb^{-1} |z| < 10 cm");
653  // }
654  // else
655  // {
656  // ltxt.DrawLatex(0.2, 0.86, "0-10% Au + Au #sqrt{s_{_{NN}}}=200 GeV");
657  // ltxt.DrawLatex(0.2, 0.80, "10B events |z| < 10 cm");
658  // }
659  // ltxt.DrawLatex(0.2, 0.74, "anti-k_{T}, R = 0.4");
660  // ltxt.DrawLatex(0.2, 0.68, Form("p_{T, 1} > % .0f GeV", pT1min));
661  // ltxt.DrawLatex(0.2, 0.62, Form("p_{T, 2} > % .0f GeV", pT2min));
662  // ltxt.DrawLatex(0.2, 0.56, Form(" |#eta| < %.0f", etamax));
663  // ltxt.DrawLatex(0.2, 0.50, Form(" |#Delta#phi_{12}| > 2#pi/3"));
664 
665 
666  //==========================================================================//
667  // PRINT PLOTS
668  //==========================================================================//
669  if (print_plots)
670  {
671  if (is_pp)
672  {
673  // cjet->Print("bjet_kinematics_pp.pdf");
674  // cdijet2->Print("dibjet_2panel_pp.pdf");
675  cdijet->Print("dibjet_imbalance_pp.pdf");
676  cdijet->Print("dibjet_imbalance_pp.C");
677  cdijet->Print("dibjet_imbalance_pp.root");
678 
679  TFile *fout = new TFile("dibjet_imbalance_histos_pp.root","RECREATE");
680  fout->cd();
681  hdijet_xj->Write();
682  hdijet_xj_fix->Write();
683  fout->Close();
684  delete fout;
685 
686  }
687  else
688  {
689  cdijet->Print("dibjet_imbalance_AuAu0-10.pdf");
690  cdijet->Print("dibjet_imbalance_AuAu0-10.C");
691  cdijet->Print("dibjet_imbalance_AuAu0-10.root");
692 
693  TFile *fout = new TFile("dibjet_imbalance_histos_AuAu0-10.root","RECREATE");
694  fout->cd();
695  hdijet_xj->Write();
696  hdijet_xj_fix->Write();
697  fout->Close();
698  delete fout;
699  }
700  }
701 
702 }