Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
D0EffScan.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file D0EffScan.C
1 #include <cstdlib>
2 #include <iostream>
3 #include <map>
4 #include <string>
5 #include <vector>
6 
7 #include "TChain.h"
8 #include "TFile.h"
9 #include "TTree.h"
10 #include "TString.h"
11 #include "TObjString.h"
12 #include "TSystem.h"
13 #include "TROOT.h"
14 
15 /*******************************/
16 /* A*E Scan for D0/D0bar */
17 /* Thomas Marshall, UCLA, 2022 */
18 /* rosstom@ucla.edu */
19 /*******************************/
20 
21 void D0EffScan()
22 {
23  // input directory where .root files are from KFParticle and QAG4SimulationTruthDecay output
24  string inDir = "/sphenix/user/rosstom/analysis/HF-Particle/KFParticle_sPHENIX/Run40Acceptance082922/";
25  TString truthFilePath;
26  TString KFPFilePath;
27 
28  Int_t nFiles = 237;
29  string fileModuleName = "Run40Acceptance082922";
30 
31  Float_t etaAccept = 0.5; // cut on each track abs(eta) to be within acceptance
32  Float_t pTAccept = 0.2; // cut on each track pT to be within acceptance
33 
34  Float_t totalTruthEvents = 0;
35  Float_t insideKinematic = 0;
36  Float_t insideGeometric = 0;
37  Float_t insideBoth = 0;
38  Float_t totalReconstructed = 0;
39 
40  // pT binning for reconstruction efficiency measurement
41  // Total number of truth events
42  Float_t truthTotal_01 = 0;
43  Float_t truthTotal_12 = 0;
44  Float_t truthTotal_23 = 0;
45  Float_t truthTotal_34 = 0;
46  Float_t truthTotal_45 = 0;
47  Float_t truthTotal_5 = 0;
48 
49  // Number of truth events within geometric/kinematic acceptance
50  Float_t truthAccept_01 = 0;
51  Float_t truthAccept_12 = 0;
52  Float_t truthAccept_23 = 0;
53  Float_t truthAccept_34 = 0;
54  Float_t truthAccept_45 = 0;
55  Float_t truthAccept_5 = 0;
56 
57  // Number of events that were reconstructed by KFParticle
58  Float_t reconstructed_01 = 0;
59  Float_t reconstructed_12 = 0;
60  Float_t reconstructed_23 = 0;
61  Float_t reconstructed_34 = 0;
62  Float_t reconstructed_45 = 0;
63  Float_t reconstructed_5 = 0;
64 
65  // eta binning for geometric acceptance measurement
66  // Uses same naming conventions as above
67  Float_t truthTotal_n11n1 = 0;
68  Float_t truthTotal_n1n08 = 0;
69  Float_t truthTotal_n08n06 = 0;
70  Float_t truthTotal_n06n04 = 0;
71  Float_t truthTotal_n04n02 = 0;
72  Float_t truthTotal_n02n00 = 0;
73  Float_t truthTotal_p11p1 = 0;
74  Float_t truthTotal_p1p08 = 0;
75  Float_t truthTotal_p08p06 = 0;
76  Float_t truthTotal_p06p04 = 0;
77  Float_t truthTotal_p04p02 = 0;
78  Float_t truthTotal_p02p00 = 0;
79  Float_t truthTotal_p11 = 0;
80  Float_t truthTotal_n11 = 0;
81 
82  Float_t truthAccept_n11n1 = 0;
83  Float_t truthAccept_n1n08 = 0;
84  Float_t truthAccept_n08n06 = 0;
85  Float_t truthAccept_n06n04 = 0;
86  Float_t truthAccept_n04n02 = 0;
87  Float_t truthAccept_n02n00 = 0;
88  Float_t truthAccept_p11p1 = 0;
89  Float_t truthAccept_p1p08 = 0;
90  Float_t truthAccept_p08p06 = 0;
91  Float_t truthAccept_p06p04 = 0;
92  Float_t truthAccept_p04p02 = 0;
93  Float_t truthAccept_p02p00 = 0;
94  Float_t truthAccept_p11 = 0;
95  Float_t truthAccept_n11 = 0;
96 
97  Float_t reconstructed_n11n1 = 0;
98  Float_t reconstructed_n1n08 = 0;
99  Float_t reconstructed_n08n06 = 0;
100  Float_t reconstructed_n06n04 = 0;
101  Float_t reconstructed_n04n02 = 0;
102  Float_t reconstructed_n02n00 = 0;
103  Float_t reconstructed_p11p1 = 0;
104  Float_t reconstructed_p1p08 = 0;
105  Float_t reconstructed_p08p06 = 0;
106  Float_t reconstructed_p06p04 = 0;
107  Float_t reconstructed_p04p02 = 0;
108  Float_t reconstructed_p02p00 = 0;
109  Float_t reconstructed_p11 = 0;
110  Float_t reconstructed_n11 = 0;
111 
112 
113  // Loop over output files from each subjob to keep reconstructed events with their corresponding truth events
114  for (int i = 235; i < nFiles; ++i)
115  {
116  std::cout << "Starting File Number: " << i << std::endl;
117  string fNumb;
118  if (i < 10)
119  {
120  fNumb = "00" + std::to_string(i);
121  }
122  else if (i < 100)
123  {
124  fNumb = "0" + std::to_string(i);
125  }
126  else
127  {
128  fNumb = std::to_string(i);
129  }
130 
131  string truthFile = "outputData_" + fileModuleName + "_00" + fNumb + ".root";
132  string KFPFile = "outputData_" + fileModuleName + "_KFP_00" + fNumb + ".root";
133  truthFilePath = inDir + truthFile;
134  KFPFilePath = inDir + KFPFile;
135 
136  TFile *truthInput(0);
137  truthInput = TFile::Open(truthFilePath);
138  TTree *truthTree = (TTree*)truthInput->Get("QAG4SimulationTruthDecay");
139  TFile *KFPInput(0);
140  KFPInput = TFile::Open(KFPFilePath);
141  TTree *recoTree = (TTree*)KFPInput->Get("DecayTree");
142 
143  Float_t truthMotherPx, truthMotherPy, truthMotherPz, truthMotherPe;
144  Float_t truthT1Px, truthT1Py, truthT1Pz, truthT1Pe;
145  Float_t truthT2Px, truthT2Py, truthT2Pz, truthT2Pe;
146  Float_t truthT1Pt, truthT2Pt;
147  Float_t truthT1Eta, truthT2Eta;
148  Float_t recoMotherPx, recoMotherPy, recoMotherPz, recoMotherPe;
149  Float_t KFPTruthT1Px, KFPTruthT1Py, KFPTruthT1Pz, KFPTruthT1Pe;
150  Float_t KFPTruthT2Px, KFPTruthT2Py, KFPTruthT2Pz, KFPTruthT2Pe;
151  Float_t truthMotherPt, truthMotherEta;
152  Int_t truthMotherBarcode;
153  Float_t recoMotherMass;
154 
155  truthTree->SetBranchAddress("mother_px", &truthMotherPx);
156  truthTree->SetBranchAddress("mother_py", &truthMotherPy);
157  truthTree->SetBranchAddress("mother_pz", &truthMotherPz);
158  truthTree->SetBranchAddress("mother_pE", &truthMotherPe);
159  truthTree->SetBranchAddress("mother_pT", &truthMotherPt);
160  truthTree->SetBranchAddress("mother_eta", &truthMotherEta);
161  truthTree->SetBranchAddress("mother_barcode", &truthMotherBarcode);
162  truthTree->SetBranchAddress("track_1_px", &truthT1Px);
163  truthTree->SetBranchAddress("track_1_py", &truthT1Py);
164  truthTree->SetBranchAddress("track_1_pz", &truthT1Pz);
165  truthTree->SetBranchAddress("track_1_pT", &truthT1Pt);
166  truthTree->SetBranchAddress("track_2_px", &truthT2Px);
167  truthTree->SetBranchAddress("track_2_py", &truthT2Py);
168  truthTree->SetBranchAddress("track_2_pz", &truthT2Pz);
169  truthTree->SetBranchAddress("track_2_pT", &truthT2Pt);
170  truthTree->SetBranchAddress("track_1_eta", &truthT1Eta);
171  truthTree->SetBranchAddress("track_2_eta", &truthT2Eta);
172 
173  recoTree->SetBranchAddress("D0_px", &recoMotherPx);
174  recoTree->SetBranchAddress("D0_py", &recoMotherPy);
175  recoTree->SetBranchAddress("D0_pz", &recoMotherPz);
176  recoTree->SetBranchAddress("D0_pE", &recoMotherPe);
177  recoTree->SetBranchAddress("D0_mass", &recoMotherMass);
178  recoTree->SetBranchAddress("track_1_true_px", &KFPTruthT1Px);
179  recoTree->SetBranchAddress("track_1_true_py", &KFPTruthT1Py);
180  recoTree->SetBranchAddress("track_1_true_pz", &KFPTruthT1Pz);
181  recoTree->SetBranchAddress("track_2_true_px", &KFPTruthT2Px);
182  recoTree->SetBranchAddress("track_2_true_py", &KFPTruthT2Py);
183  recoTree->SetBranchAddress("track_2_true_pz", &KFPTruthT2Pz);
184 
185  vector<Int_t> reconstructedBarcodes;
186  vector<Int_t> missedBarcodes;
187  vector<int> usedRecoEntries;
188 
189  Float_t min_deltaPx, min_deltaPy, min_deltaPz, min_deltaPe;
190  int minEntry;
191 
192  for (int j = 0; j < truthTree->GetEntries(); ++j)
193  {
194  truthTree->GetEntry(j);
195  totalTruthEvents += 1;
196 
197  if (truthT1Pt > pTAccept && truthT2Pt > pTAccept) insideKinematic += 1;
198  if (abs(truthT1Eta) < etaAccept && abs(truthT2Eta) < etaAccept) insideGeometric += 1;
199  if (abs(truthT1Eta) < etaAccept && abs(truthT2Eta) < etaAccept && truthT1Pt > pTAccept && truthT2Pt > pTAccept) insideBoth += 1;
200 
201  if (truthMotherPt < 1) truthTotal_01 += 1;
202  if (truthMotherPt < 2 && truthMotherPt >= 1) truthTotal_12 += 1;
203  if (truthMotherPt < 3 && truthMotherPt >= 2) truthTotal_23 += 1;
204  if (truthMotherPt < 4 && truthMotherPt >= 3) truthTotal_34 += 1;
205  if (truthMotherPt < 5 && truthMotherPt >= 4) truthTotal_45 += 1;
206  if (truthMotherPt >= 5) truthTotal_5 += 1;
207 
208  if (truthMotherEta < -1) truthTotal_n11 += 1;
209  if (truthMotherEta < -1 && truthMotherEta >= -1.1) truthTotal_n11n1 += 1;
210  if (truthMotherEta < -0.8 && truthMotherEta >= -1) truthTotal_n1n08 += 1;
211  if (truthMotherEta < -0.6 && truthMotherEta >= -0.8) truthTotal_n08n06 += 1;
212  if (truthMotherEta < -0.4 && truthMotherEta >= -0.6) truthTotal_n06n04 += 1;
213  if (truthMotherEta < -0.2 && truthMotherEta >= -0.4) truthTotal_n04n02 += 1;
214  if (truthMotherEta < 0.0 && truthMotherEta >= -0.2) truthTotal_n02n00 += 1;
215  if (truthMotherEta < 0.2 && truthMotherEta >= 0.0) truthTotal_p02p00 += 1;
216  if (truthMotherEta < 0.4 && truthMotherEta >= 0.2) truthTotal_p04p02 += 1;
217  if (truthMotherEta < 0.6 && truthMotherEta >= 0.4) truthTotal_p06p04 += 1;
218  if (truthMotherEta < 0.8 && truthMotherEta >= 0.6) truthTotal_p08p06 += 1;
219  if (truthMotherEta < 1.0 && truthMotherEta >= 0.8) truthTotal_p1p08 += 1;
220  if (truthMotherEta < 1.1 && truthMotherEta >= 1) truthTotal_p11p1 += 1;
221  if (truthMotherEta >= 1) truthTotal_p11 += 1;
222 
223  if (abs(truthT1Eta) < etaAccept && abs(truthT2Eta) < etaAccept && truthT1Pt > pTAccept && truthT2Pt > pTAccept)
224  {
225  if (truthMotherPt < 1) truthAccept_01 += 1;
226  if (truthMotherPt < 2 && truthMotherPt >= 1) truthAccept_12 += 1;
227  if (truthMotherPt < 3 && truthMotherPt >= 2) truthAccept_23 += 1;
228  if (truthMotherPt < 4 && truthMotherPt >= 3) truthAccept_34 += 1;
229  if (truthMotherPt < 5 && truthMotherPt >= 4) truthAccept_45 += 1;
230  if (truthMotherPt >= 5) truthAccept_5 += 1;
231 
232  if (truthMotherEta < -1) truthAccept_n11 += 1;
233  if (truthMotherEta < -1 && truthMotherEta >= -1.1) truthAccept_n11n1 += 1;
234  if (truthMotherEta < -0.8 && truthMotherEta >= -1) truthAccept_n1n08 += 1;
235  if (truthMotherEta < -0.6 && truthMotherEta >= -0.8) truthAccept_n08n06 += 1;
236  if (truthMotherEta < -0.4 && truthMotherEta >= -0.6) truthAccept_n06n04 += 1;
237  if (truthMotherEta < -0.2 && truthMotherEta >= -0.4) truthAccept_n04n02 += 1;
238  if (truthMotherEta < 0.0 && truthMotherEta >= -0.2) truthAccept_n02n00 += 1;
239  if (truthMotherEta < 0.2 && truthMotherEta >= 0.0) truthAccept_p02p00 += 1;
240  if (truthMotherEta < 0.4 && truthMotherEta >= 0.2) truthAccept_p04p02 += 1;
241  if (truthMotherEta < 0.6 && truthMotherEta >= 0.4) truthAccept_p06p04 += 1;
242  if (truthMotherEta < 0.8 && truthMotherEta >= 0.6) truthAccept_p08p06 += 1;
243  if (truthMotherEta < 1.0 && truthMotherEta >= 0.8) truthAccept_p1p08 += 1;
244  if (truthMotherEta < 1.1 && truthMotherEta >= 1) truthAccept_p11p1 += 1;
245  if (truthMotherEta >= 1) truthAccept_p11 += 1;
246  }
247 
248  bool matchPx_t1r1;
249  bool matchPy_t1r1;
250  bool matchPz_t1r1;
251  bool matchPx_t1r2;
252  bool matchPy_t1r2;
253  bool matchPz_t1r2;
254  bool matchPx_t2r1;
255  bool matchPy_t2r1;
256  bool matchPz_t2r1;
257  bool matchPx_t2r2;
258  bool matchPy_t2r2;
259  bool matchPz_t2r2;
260 
261  bool matches = false;
262 
263  for (int k = 0; k < recoTree->GetEntries(); ++k)
264  {
265  recoTree->GetEntry(k);
266  matchPx_t1r1 = (KFPTruthT1Px == truthT1Px);
267  matchPy_t1r1 = (KFPTruthT1Py == truthT1Py);
268  matchPz_t1r1 = (KFPTruthT1Pz == truthT1Pz);
269  matchPx_t1r2 = (KFPTruthT2Px == truthT1Px);
270  matchPy_t1r2 = (KFPTruthT2Py == truthT1Py);
271  matchPz_t1r2 = (KFPTruthT2Pz == truthT1Pz);
272  matchPx_t2r1 = (KFPTruthT1Px == truthT2Px);
273  matchPy_t2r1 = (KFPTruthT1Py == truthT2Py);
274  matchPz_t2r1 = (KFPTruthT1Pz == truthT2Pz);
275  matchPx_t2r2 = (KFPTruthT2Px == truthT2Px);
276  matchPy_t2r2 = (KFPTruthT2Py == truthT2Py);
277  matchPz_t2r2 = (KFPTruthT2Pz == truthT2Pz);
278 
279  // check the truth information from KFParticle and QAG4SimulationTruthDecay to match a reconstructed candidate with its truth information
280  // if match is found, that means a true D0->PiK event was reconstructed in KFParticle
281  if (matchPx_t1r1 && matchPy_t1r1 && matchPz_t1r1)
282  {
283  if (matchPx_t2r2 && matchPy_t2r2 && matchPz_t2r2)
284  {
285  std::cout << "Got a match" << std::endl;
286  matches = true;
287  minEntry = k;
288  break;
289  }
290  }
291  else if (matchPx_t1r2 && matchPy_t1r2 && matchPz_t1r2)
292  {
293  if (matchPx_t2r1 && matchPy_t2r1 && matchPz_t2r1)
294  {
295  std::cout << "Got a match" << std::endl;
296  matches = true;
297  minEntry = k;
298  break;
299  }
300  }
301  }
302  if (matches)
303  {
304  if (usedRecoEntries.size() > 0)
305  {
306  bool recoCandidateUsed = false;
307  for (int ent : usedRecoEntries)
308  {
309  if (ent == minEntry)
310  {
311  std::cout << "Candidate has already been used, skipping this candidate" << std::endl;
312  recoCandidateUsed = true;
313  }
314  }
315  if (recoCandidateUsed) continue;
316  else
317  {
318  usedRecoEntries.push_back(minEntry);
319 
320  if (abs(truthT1Eta) < etaAccept && abs(truthT2Eta) < etaAccept && truthT1Pt > pTAccept && truthT2Pt > pTAccept)
321  {
322  totalReconstructed += 1;
323 
324  if (truthMotherPt < 1) reconstructed_01 += 1;
325  if (truthMotherPt < 2 && truthMotherPt >= 1) reconstructed_12 += 1;
326  if (truthMotherPt < 3 && truthMotherPt >= 2) reconstructed_23 += 1;
327  if (truthMotherPt < 4 && truthMotherPt >= 3) reconstructed_34 += 1;
328  if (truthMotherPt < 5 && truthMotherPt >= 4) reconstructed_45 += 1;
329  if (truthMotherPt >= 5) reconstructed_5 += 1;
330 
331  if (truthMotherEta < -1.1) reconstructed_n11 += 1;
332  if (truthMotherEta < -1 && truthMotherEta >= -1.1) reconstructed_n11n1 += 1;
333  if (truthMotherEta < -0.8 && truthMotherEta >= -1) reconstructed_n1n08 += 1;
334  if (truthMotherEta < -0.6 && truthMotherEta >= -0.8) reconstructed_n08n06 += 1;
335  if (truthMotherEta < -0.4 && truthMotherEta >= -0.6) reconstructed_n06n04 += 1;
336  if (truthMotherEta < -0.2 && truthMotherEta >= -0.4) reconstructed_n04n02 += 1;
337  if (truthMotherEta < 0.0 && truthMotherEta >= -0.2) reconstructed_n02n00 += 1;
338  if (truthMotherEta < 0.2 && truthMotherEta >= 0.0) reconstructed_p02p00 += 1;
339  if (truthMotherEta < 0.4 && truthMotherEta >= 0.2) reconstructed_p04p02 += 1;
340  if (truthMotherEta < 0.6 && truthMotherEta >= 0.4) reconstructed_p06p04 += 1;
341  if (truthMotherEta < 0.8 && truthMotherEta >= 0.6) reconstructed_p08p06 += 1;
342  if (truthMotherEta < 1.0 && truthMotherEta >= 0.8) reconstructed_p1p08 += 1;
343  if (truthMotherEta < 1.1 && truthMotherEta >= 1) reconstructed_p11p1 += 1;
344  if (truthMotherEta >= 1) reconstructed_p11 += 1;
345  }
346  }
347  }
348  else
349  {
350  usedRecoEntries.push_back(minEntry);
351 
352  if (abs(truthT1Eta) < etaAccept && abs(truthT2Eta) < etaAccept && truthT1Pt > pTAccept && truthT2Pt > pTAccept)
353  {
354  totalReconstructed += 1;
355 
356  if (truthMotherPt < 1) reconstructed_01 += 1;
357  if (truthMotherPt < 2 && truthMotherPt >= 1) reconstructed_12 += 1;
358  if (truthMotherPt < 3 && truthMotherPt >= 2) reconstructed_23 += 1;
359  if (truthMotherPt < 4 && truthMotherPt >= 3) reconstructed_34 += 1;
360  if (truthMotherPt < 5 && truthMotherPt >= 4) reconstructed_45 += 1;
361  if (truthMotherPt >= 5) reconstructed_5 += 1;
362 
363  if (truthMotherEta < -1.1) reconstructed_n11 += 1;
364  if (truthMotherEta < -1 && truthMotherEta >= -1.1) reconstructed_n11n1 += 1;
365  if (truthMotherEta < -0.8 && truthMotherEta >= -1) reconstructed_n1n08 += 1;
366  if (truthMotherEta < -0.6 && truthMotherEta >= -0.8) reconstructed_n08n06 += 1;
367  if (truthMotherEta < -0.4 && truthMotherEta >= -0.6) reconstructed_n06n04 += 1;
368  if (truthMotherEta < -0.2 && truthMotherEta >= -0.4) reconstructed_n04n02 += 1;
369  if (truthMotherEta < 0.0 && truthMotherEta >= -0.2) reconstructed_n02n00 += 1;
370  if (truthMotherEta < 0.2 && truthMotherEta >= 0.0) reconstructed_p02p00 += 1;
371  if (truthMotherEta < 0.4 && truthMotherEta >= 0.2) reconstructed_p04p02 += 1;
372  if (truthMotherEta < 0.6 && truthMotherEta >= 0.4) reconstructed_p06p04 += 1;
373  if (truthMotherEta < 0.8 && truthMotherEta >= 0.6) reconstructed_p08p06 += 1;
374  if (truthMotherEta < 1.0 && truthMotherEta >= 0.8) reconstructed_p1p08 += 1;
375  if (truthMotherEta < 1.1 && truthMotherEta >= 1) reconstructed_p11p1 += 1;
376  if (truthMotherEta >= 1) reconstructed_p11 += 1;
377  }
378  }
379  }
380  }
381  }
382  std::cout << "truth mother pT efficiencies" << std::endl;
383  std::cout << "0-1 GeV: " << reconstructed_01/truthAccept_01 << std::endl;
384  std::cout << "1-2 GeV: " << reconstructed_12/truthAccept_12 << std::endl;
385  std::cout << "2-3 GeV: " << reconstructed_23/truthAccept_23 << std::endl;
386  std::cout << "3-4 GeV: " << reconstructed_34/truthAccept_34 << std::endl;
387  std::cout << "4-5 GeV: " << reconstructed_45/truthAccept_45 << std::endl;
388  std::cout << "5+ GeV: " << reconstructed_5/truthAccept_5 << std::endl;
389 
390  std::cout << "" << std::endl;
391  std::cout << "Number of truth events w/in truth track acceptance" << std::endl;
392 
393  std::cout << "0-1 GeV: " << truthAccept_01 << std::endl;
394  std::cout << "1-2 GeV: " << truthAccept_12 << std::endl;
395  std::cout << "2-3 GeV: " << truthAccept_23 << std::endl;
396  std::cout << "3-4 GeV: " << truthAccept_34 << std::endl;
397  std::cout << "4-5 GeV: " << truthAccept_45 << std::endl;
398  std::cout << "5+ GeV: " << truthAccept_5 << std::endl;
399 
400  std::cout << "" << std::endl;
401 
402  std::cout << "truth mother eta acceptance" << std::endl;
403  std::cout << "eta < -1.1 : " << reconstructed_n11/truthAccept_n11 << std::endl;
404  std::cout << "-1.1 <= eta < -1.0 : " << reconstructed_n11n1/truthAccept_n11n1 << std::endl;
405  std::cout << "-1 <= eta < -0.8 : " << reconstructed_n1n08/truthAccept_n1n08 << std::endl;
406  std::cout << "-0.8 <= eta < -0.6 : " << reconstructed_n08n06/truthAccept_n08n06 << std::endl;
407  std::cout << "-0.6 <= eta < -0.4 : " << reconstructed_n06n04/truthAccept_n06n04 << std::endl;
408  std::cout << "-0.4 <= eta < -0.2 : " << reconstructed_n04n02/truthAccept_n04n02 << std::endl;
409  std::cout << "-0.2 <= eta < 0.0 : " << reconstructed_n02n00/truthAccept_n02n00 << std::endl;
410  std::cout << "0.0 <= eta < 0.2 : " << reconstructed_p02p00/truthAccept_p02p00 << std::endl;
411  std::cout << "0.2 <= eta < 0.4 : " << reconstructed_p04p02/truthAccept_p04p02 << std::endl;
412  std::cout << "0.4 <= eta < 0.6 : " << reconstructed_p06p04/truthAccept_p06p04 << std::endl;
413  std::cout << "0.6 <= eta < 0.8 : " << reconstructed_p08p06/truthAccept_p08p06 << std::endl;
414  std::cout << "0.8 <= eta < 1.0 : " << reconstructed_p1p08/truthAccept_p1p08 << std::endl;
415  std::cout << "1.0 <= eta < 1.1 : " << reconstructed_p11p1/truthAccept_p11p1 << std::endl;
416  std::cout << "eta >= 1.1 : " << reconstructed_p11/truthAccept_p11 << std::endl;
417 
418  std::cout << "" << std::endl;
419  std::cout << "Number of truth events w/in truth track acceptance" << std::endl;
420 
421  std::cout << "eta < -1.1 : " << truthAccept_n11 << std::endl;
422  std::cout << "-1.1 <= eta < -1.0 : " << truthAccept_n11n1 << std::endl;
423  std::cout << "-1 <= eta < -0.8 : " << truthAccept_n1n08 << std::endl;
424  std::cout << "-0.8 <= eta < -0.6 : " << truthAccept_n08n06 << std::endl;
425  std::cout << "-0.6 <= eta < -0.4 : " << truthAccept_n06n04 << std::endl;
426  std::cout << "-0.4 <= eta < -0.2 : " << truthAccept_n04n02 << std::endl;
427  std::cout << "-0.2 <= eta < 0.0 : " << truthAccept_n02n00 << std::endl;
428  std::cout << "0.0 <= eta < 0.2 : " << truthAccept_p02p00 << std::endl;
429  std::cout << "0.2 <= eta < 0.4 : " << truthAccept_p04p02 << std::endl;
430  std::cout << "0.4 <= eta < 0.6 : " << truthAccept_p06p04 << std::endl;
431  std::cout << "0.6 <= eta < 0.8 : " << truthAccept_p08p06 << std::endl;
432  std::cout << "0.8 <= eta < 1.0 : " << truthAccept_p1p08 << std::endl;
433  std::cout << "1.0 <= eta < 1.1 : " << truthAccept_p11p1 << std::endl;
434  std::cout << "eta >= 1.1 : " << truthAccept_p11 << std::endl;
435 
436  std::cout << "Total Number of Truth Events: " << totalTruthEvents << std::endl;
437  std::cout << "Fraction of events w/ truth track pT > " << pTAccept << " GeV: " << insideKinematic/totalTruthEvents << std::endl;
438  std::cout << "Fraction of events w/ truth track |eta| < " << etaAccept << ": " << insideGeometric/totalTruthEvents << std::endl;
439  std::cout << "Fraction of events w/ both: " << insideBoth/totalTruthEvents << std::endl;
440  std::cout << "Total Integrated Efficiency: " << totalReconstructed/insideBoth << std::endl;
441 
442  const Int_t n_pT = 6;
443  Double_t pT_bin[n_pT] = {0.5,1.5,2.5,3.5,4.5,5.5};
444  Double_t pT_efficiency[n_pT] = {reconstructed_01/truthAccept_01,reconstructed_12/truthAccept_12,reconstructed_23/truthAccept_23,reconstructed_34/truthAccept_34,reconstructed_45/truthAccept_45,reconstructed_5/truthAccept_5};
445  Double_t pT_bin_widths[n_pT] = {0.5,0.5,0.5,0.5,0.5,0.5};
446  Double_t pT_efficiency_err[n_pT] = {sqrt(truthAccept_01)/truthAccept_01,sqrt(truthAccept_12)/truthAccept_12,sqrt(truthAccept_23)/truthAccept_23,sqrt(truthAccept_34)/truthAccept_34,sqrt(truthAccept_45)/truthAccept_45,sqrt(truthAccept_5)/truthAccept_5};
447 
448  TCanvas *c1 = new TCanvas("D0_Efficiency_pT", "D0_Efficiency_pT", 1800, 900);
449  TGraphErrors* D0_Efficiency_pT = new TGraphErrors(n_pT, pT_bin, pT_efficiency, pT_bin_widths, pT_efficiency_err);
450  D0_Efficiency_pT->SetTitle("D^{0} Reconstruction Efficiency vs Truth Mother p_{T};Truth p_{T} [GeV/c];Reconstruction Efficiency");
451  D0_Efficiency_pT->SetLineColor(kBlue);
452  D0_Efficiency_pT->GetXaxis()->SetLimits(0.,6.0);
453  D0_Efficiency_pT->SetMinimum(0.);
454  D0_Efficiency_pT->SetMaximum(1.0);
455  D0_Efficiency_pT->SetMarkerStyle(22);
456  D0_Efficiency_pT->SetMarkerColor(kBlue);
457  D0_Efficiency_pT->SetMarkerSize(2);
458  D0_Efficiency_pT->SetLineWidth(2);
459  D0_Efficiency_pT->Draw("ACP");
460  auto legend = new TLegend(0.1,0.7,0.48,0.9);
461  legend->AddEntry("","#it{#bf{sPHENIX}} Simulation","");
462  legend->AddEntry("", "p+p #rightarrow D^{0} (#pi^{+}K^{-}) or #bar{D}^{0} (#pi^{-}K^{+})", "");
463  legend->AddEntry("", "#sqrt{s}=200 GeV", "");
464  legend->AddEntry("", "True Track pT > 0.2 GeV, |#eta| < 0.5", "");
465  legend->SetMargin(0);
466  legend->Draw();
467  c1->SaveAs("D0_Efficiency_pT.png");
468  c1->Close();
469 
470  const Int_t n_eta = 14;
471  Double_t eta_bin[n_eta] = {-1.15,-1.05,-0.9,-0.7,-0.5,-0.3,-0.1,0.1,0.3,0.5,0.7,0.9,1.05,1.15};
472  Double_t eta_efficiency[n_eta] = {reconstructed_n11/truthAccept_n11, reconstructed_n11n1/truthAccept_n11n1, reconstructed_n1n08/truthAccept_n1n08, reconstructed_n08n06/truthAccept_n08n06,
473  reconstructed_n06n04/truthAccept_n06n04, reconstructed_n04n02/truthAccept_n04n02, reconstructed_n02n00/truthAccept_n02n00, reconstructed_p02p00/truthAccept_p02p00,
474  reconstructed_p04p02/truthAccept_p04p02, reconstructed_p06p04/truthAccept_p06p04, reconstructed_p08p06/truthAccept_p08p06, reconstructed_p1p08/truthAccept_p1p08,
475  reconstructed_p11p1/truthAccept_p11p1, reconstructed_p11/truthAccept_p11};
476  Double_t eta_bin_widths[n_eta] = {0.05,0.05,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.05,0.05};
477  Double_t eta_efficiency_err[n_eta] = {sqrt(truthAccept_n11)/truthAccept_n11,sqrt(truthAccept_n11n1)/truthAccept_n11n1,sqrt(truthAccept_n1n08)/truthAccept_n1n08,sqrt(truthAccept_n08n06)/truthAccept_n08n06,
478  sqrt(truthAccept_n06n04)/truthAccept_n06n04,sqrt(truthAccept_n04n02)/truthAccept_n04n02,sqrt(truthAccept_n02n00)/truthAccept_n02n00,sqrt(truthAccept_p02p00)/truthAccept_p02p00,
479  sqrt(truthAccept_p04p02)/truthAccept_p04p02,sqrt(truthAccept_p06p04)/truthAccept_p06p04,sqrt(truthAccept_p08p06)/truthAccept_p08p06,sqrt(truthAccept_p1p08)/truthAccept_p1p08,
480  sqrt(truthAccept_p11p1)/truthAccept_p11p1,sqrt(truthAccept_p11)/truthAccept_p11};
481 
482  TCanvas *c2 = new TCanvas("D0_Efficiency_eta", "D0_Efficiency_eta", 1800, 900);
483  TGraphErrors* D0_Efficiency_eta = new TGraphErrors(n_eta, eta_bin, eta_efficiency, eta_bin_widths, eta_efficiency_err);
484  D0_Efficiency_eta->SetTitle("D^{0} Reconstruction Efficiency vs Truth Mother #eta;Truth #eta;Reconstruction Efficiency");
485  D0_Efficiency_eta->SetLineColor(kBlue);
486  D0_Efficiency_eta->GetXaxis()->SetLimits(-1.25,1.25);
487  D0_Efficiency_eta->SetMinimum(0.);
488  D0_Efficiency_eta->SetMaximum(1.0);
489  D0_Efficiency_eta->SetMarkerStyle(22);
490  D0_Efficiency_eta->SetMarkerColor(kBlue);
491  D0_Efficiency_eta->SetMarkerSize(2);
492  D0_Efficiency_eta->SetLineWidth(2);
493  D0_Efficiency_eta->Draw("ACP");
494  auto legend2 = new TLegend(0.1,0.7,0.48,0.9);
495  legend2->AddEntry("","#it{#bf{sPHENIX}} Simulation","");
496  legend2->AddEntry("", "p+p #rightarrow D^{0} (#pi^{+}K^{-}) or #bar{D}^{0} (#pi^{-}K^{+})", "");
497  legend2->AddEntry("", "#sqrt{s}=200 GeV", "");
498  legend2->AddEntry("", "True Track pT > 0.2 GeV, |#eta| < 0.5", "");
499  legend2->SetMargin(0);
500  legend2->Draw();
501  c2->SaveAs("D0_Efficiency_eta.png");
502  c2->Close();
503 }