Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Draw_BDiJetImbalance_ncoll.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Draw_BDiJetImbalance_ncoll.C
1 
2 //
3 // Draw the Di b-jet <xj> vs Ncoll
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 // 29 Jun 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 #include <TGraphErrors.h>
28 
29 #include <iostream>
30 
31 using namespace std;
32 
33 void calc_mxj(TH1D* htmp, double &mxj, double &mxj_e)
34 {
35 
36  double sum = 0;
37  double w = 0;
38  double err = 0;
39  for (int i = 1; i <= htmp->GetNbinsX(); i++)
40  {
41  double c = htmp->GetBinContent(i);
42  if (c > 0)
43  {
44  sum += c * htmp->GetBinCenter(i);
45  w += c;
46  err += TMath::Power(c * htmp->GetBinError(i), 2);
47  }
48  }
49  double mean = sum / w;
50  err = TMath::Sqrt(err) / w;
51 
52  mxj = mean;
53  mxj_e = err;
54 }
55 
56 void Draw_BDiJetImbalance_ncoll()
57 {
58  //==========================================================================//
59  // SET RUNNING CONDITIONS
60  //==========================================================================//
61  bool print_plots = true;
62 
63  const char* ppFile = "dibjet_imbalance_histos_pp.root";
64  const char* auauFile = "dibjet_imbalance_histos_AuAu0-10.root";
65 
66  const int NCENT = 4;
67  const char* centlim[] = {"0-10%", "10-20%", "20-40%", "40-60%", "60-92%"};
68  double ncoll[] = {962, 603, 296, 94, 15};
69  double lcent[] = {549, 344, 338, 107, 27}; // nn luminosity (for scaling)
70 
71 
72  //==========================================================================//
73  // DECLARE VARIABLES
74  //==========================================================================//
75 
76  TH1D* hxj_pp;
77  TH1D* hxj_auau[NCENT];
78 
79  TGraphErrors *gmxj = new TGraphErrors();
80  gmxj->SetMarkerStyle(kFullCircle);
81  gmxj->SetMarkerColor(kBlack);
82  gmxj->SetLineColor(kBlack);
83 
84 
85  double mxj_pp;
86  double mxj_e_pp;
87 
88  double mxj_auau[NCENT];
89  double mxj_e_auau[NCENT];
90 
91  //-- other
92  TFile *fin;
93  char name[500];
94 
95  //==========================================================================//
96  // READ RESULTS FROM FILE
97  //==========================================================================//
98  cout << endl;
99  cout << "--> Reading results from file" << endl;
100 
101  cout << "----> Reading pp from " << ppFile << endl;
102  fin = TFile::Open(ppFile);
103  if (!fin)
104  {
105  cout << "ERROR!! Unable to read " << ppFile << endl;
106  return;
107  }
108 
109  sprintf(name, "hdijet_xj_fix");
110  hxj_pp = (TH1D*) fin->Get(name);
111  if ( !hxj_pp )
112  {
113  cout << "ERRRO! Unable to find " << name << " in " << ppFile << endl;
114  return;
115  }
116  hxj_pp->SetDirectory(0);
117  hxj_pp->SetName("hxj_pp");
118 
119  fin->Close();
120  delete fin;
121 
122 
123 
124  cout << "----> Reading auau from " << auauFile << endl;
125  fin = TFile::Open(auauFile);
126  if (!fin)
127  {
128  cout << "ERROR!! Unable to read " << auauFile << endl;
129  return;
130  }
131 
132  sprintf(name, "hdijet_xj_fix");
133  hxj_auau[0] = (TH1D*) fin->Get(name);
134  if ( !hxj_auau[0] )
135  {
136  cout << "ERRRO! Unable to find " << name << " in " << auauFile << endl;
137  return;
138  }
139  hxj_auau[0]->SetDirectory(0);
140  hxj_auau[0]->SetName("hxj_auau_0");
141 
142  fin->Close();
143  delete fin;
144 
145  //==========================================================================//
146  // SCALE UNCERTAINTIES FOR OTHER CENTRALITY BINS
147  //==========================================================================//
148  cout << endl;
149  cout << "--> Scaling AuAu centrality bins" << endl;
150 
151  for (int icent = 1; icent < NCENT; icent++)
152  {
153  sprintf(name, "hxj_auau_%i", icent);
154  hxj_auau[icent] = (TH1D*) hxj_auau[0]->Clone(name);
155 
156  double sf = 1./TMath::Sqrt(lcent[icent] / lcent[0]);
157  for (int ix = 1; ix <= hxj_auau[icent]->GetNbinsX(); ix++)
158  {
159  double be = hxj_auau[icent]->GetBinError(ix);
160 
161  hxj_auau[icent]->SetBinError(ix, be * sf);
162  } // ix
163  }
164 
165 
166  //==========================================================================//
167  // CALCULATE <XJ>
168  //==========================================================================//
169  cout << endl;
170  cout << "--> Calculating <x_j>" << endl;
171 
172  //-- pp
173  calc_mxj(hxj_pp, mxj_pp, mxj_e_pp);
174 
175  gmxj->SetPoint(0, 1, mxj_pp);
176  gmxj->SetPointError(0, 0, mxj_e_pp);
177 
178  //-- AuAu
179  for (int i = 0; i < NCENT; i++)
180  {
181  calc_mxj(hxj_auau[i], mxj_auau[i], mxj_e_auau[i]);
182 
183  gmxj->SetPoint(i + 1, ncoll[i], mxj_auau[i]);
184  gmxj->SetPointError(i + 1, 0, mxj_e_auau[i]);
185 
186  } // i
187 
188 
189  //==========================================================================//
190  // PLOT OBJECTS
191  //==========================================================================//
192 
193  TH1D* haxis = new TH1D("haxis", ";N_{coll};#LTx_{j}#GT", 1200, -100, 1100);
194  // haxis->SetMinimum(0.551);
195  // haxis->SetMaximum(0.749);
196  haxis->SetMinimum(0.601);
197  haxis->SetMaximum(0.779);
198 
199  TLatex lbl;
200  lbl.SetTextAlign(22);
201 
202  //==========================================================================//
203  // PLOT
204  //==========================================================================//
205 
206  TCanvas *cxj = new TCanvas("cxj", "xj", 800, 600);
207 
208  cxj->cd();
209  // gPad->SetLogx();
210  haxis->Draw();
211 
212  gmxj->Draw("P");
213 
214  // cent labels
215  lbl.DrawLatex(1, mxj_pp + 0.02, "p+p");
216 
217  for (int i = 0; i < NCENT; i++)
218  lbl.DrawLatex(ncoll[i], mxj_auau[i] + 0.02 + 0.01*(i%2), centlim[i]);
219 
220  TLegend *leg = new TLegend(0.10, 0.20, 0.70, 0.55);
221  leg->SetFillStyle(0);
222  leg->SetBorderSize(0);
223  leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation", "");
224  leg->AddEntry("", "#sqrt{s_{_{NN}}} = 200 GeV", "");
225  leg->AddEntry("", "Di b-jets (Pythia8, CTEQ6L)", "");
226  leg->AddEntry("", "anti-k_{T}, R = 0.4, |#eta|<0.7", "");
227  leg->AddEntry("", "p+p: 200 pb^{-1}, 60% Eff., 40% Pur.", "");
228  leg->AddEntry("", "Au+Au: 240B Events, 40% Eff., 40% Pur.", "");
229  // leg->AddEntry("", "200 pb^{-1} p+p, 240B Au+Au", "");
230  leg->Draw("same");
231 
232  //==========================================================================//
233  // PRINT PLOTS
234  //==========================================================================//
235  if ( print_plots )
236  {
237  cout << endl;
238  cout << "--> Printing results" << endl;
239 
240  cout << "p+p <xj> = " << mxj_pp << " +/- " << mxj_e_pp << endl;
241  for (int i = 0; i < NCENT; i++)
242  {
243  cout << "Au+Au " << centlim[i] << " <xj> = " << mxj_auau[i] << " +/- " << mxj_e_auau[i] << endl;
244  }
245 
246  cxj->Print("dibjet_imbalance_ncoll.pdf");
247  cxj->Print("dibjet_imbalance_ncoll.C");
248  cxj->Print("dibjet_imbalance_ncoll.root");
249  }
250 
251 
252 
253 
254 
255 }