HF-Particle/KFParticle_sPHENIX/D0EffScan.C

0001 #include <cstdlib>
0002 #include <iostream>
0003 #include <map>
0004 #include <string>
0005 #include <vector>
0006
0007 #include "TChain.h"
0008 #include "TFile.h"
0009 #include "TTree.h"
0010 #include "TString.h"
0011 #include "TObjString.h"
0012 #include "TSystem.h"
0013 #include "TROOT.h"
0014
0015 /*******************************/
0016 /*    A*E Scan for D0/D0bar    */
0017 /* Thomas Marshall, UCLA, 2022 */
0018 /*       rosstom@ucla.edu      */
0019 /*******************************/
0020
0021 void D0EffScan()
0022 {
0023   // input directory where .root files are from KFParticle and QAG4SimulationTruthDecay output
0024   string inDir = "/sphenix/user/rosstom/analysis/HF-Particle/KFParticle_sPHENIX/Run40Acceptance082922/";
0025   TString truthFilePath;
0026   TString KFPFilePath;
0027
0028   Int_t nFiles = 237;
0029   string fileModuleName = "Run40Acceptance082922";
0030
0031   Float_t etaAccept = 0.5; // cut on each track abs(eta) to be within acceptance
0032   Float_t pTAccept = 0.2; // cut on each track pT to be within acceptance
0033
0034   Float_t totalTruthEvents = 0;
0035   Float_t insideKinematic = 0;
0036   Float_t insideGeometric = 0;
0037   Float_t insideBoth = 0;
0038   Float_t totalReconstructed = 0;
0039
0040   // pT binning for reconstruction efficiency measurement
0041   // Total number of truth events
0042   Float_t truthTotal_01 = 0;
0043   Float_t truthTotal_12 = 0;
0044   Float_t truthTotal_23 = 0;
0045   Float_t truthTotal_34 = 0;
0046   Float_t truthTotal_45 = 0;
0047   Float_t truthTotal_5 = 0;
0048
0049   // Number of truth events within geometric/kinematic acceptance
0050   Float_t truthAccept_01 = 0;
0051   Float_t truthAccept_12 = 0;
0052   Float_t truthAccept_23 = 0;
0053   Float_t truthAccept_34 = 0;
0054   Float_t truthAccept_45 = 0;
0055   Float_t truthAccept_5 = 0;
0056
0057   // Number of events that were reconstructed by KFParticle
0058   Float_t reconstructed_01 = 0;
0059   Float_t reconstructed_12 = 0;
0060   Float_t reconstructed_23 = 0;
0061   Float_t reconstructed_34 = 0;
0062   Float_t reconstructed_45 = 0;
0063   Float_t reconstructed_5 = 0;
0064
0065   // eta binning for geometric acceptance measurement
0066   // Uses same naming conventions as above
0067   Float_t truthTotal_n11n1 = 0;
0068   Float_t truthTotal_n1n08 = 0;
0069   Float_t truthTotal_n08n06 = 0;
0070   Float_t truthTotal_n06n04 = 0;
0071   Float_t truthTotal_n04n02 = 0;
0072   Float_t truthTotal_n02n00 = 0;
0073   Float_t truthTotal_p11p1 = 0;
0074   Float_t truthTotal_p1p08 = 0;
0075   Float_t truthTotal_p08p06 = 0;
0076   Float_t truthTotal_p06p04 = 0;
0077   Float_t truthTotal_p04p02 = 0;
0078   Float_t truthTotal_p02p00 = 0;
0079   Float_t truthTotal_p11 = 0;
0080   Float_t truthTotal_n11 = 0;
0081
0082   Float_t truthAccept_n11n1 = 0;
0083   Float_t truthAccept_n1n08 = 0;
0084   Float_t truthAccept_n08n06 = 0;
0085   Float_t truthAccept_n06n04 = 0;
0086   Float_t truthAccept_n04n02 = 0;
0087   Float_t truthAccept_n02n00 = 0;
0088   Float_t truthAccept_p11p1 = 0;
0089   Float_t truthAccept_p1p08 = 0;
0090   Float_t truthAccept_p08p06 = 0;
0091   Float_t truthAccept_p06p04 = 0;
0092   Float_t truthAccept_p04p02 = 0;
0093   Float_t truthAccept_p02p00 = 0;
0094   Float_t truthAccept_p11 = 0;
0095   Float_t truthAccept_n11 = 0;
0096
0097   Float_t reconstructed_n11n1 = 0;
0098   Float_t reconstructed_n1n08 = 0;
0099   Float_t reconstructed_n08n06 = 0;
0100   Float_t reconstructed_n06n04 = 0;
0101   Float_t reconstructed_n04n02 = 0;
0102   Float_t reconstructed_n02n00 = 0;
0103   Float_t reconstructed_p11p1 = 0;
0104   Float_t reconstructed_p1p08 = 0;
0105   Float_t reconstructed_p08p06 = 0;
0106   Float_t reconstructed_p06p04 = 0;
0107   Float_t reconstructed_p04p02 = 0;
0108   Float_t reconstructed_p02p00 = 0;
0109   Float_t reconstructed_p11 = 0;
0110   Float_t reconstructed_n11 = 0;
0111
0112
0113   // Loop over output files from each subjob to keep reconstructed events with their corresponding truth events
0114   for (int i = 235; i < nFiles; ++i)
0115   {
0116     std::cout << "Starting File Number: " << i << std::endl;
0117     string fNumb;
0118     if (i < 10)
0119     {
0120       fNumb = "00" + std::to_string(i);
0121     }
0122     else if (i < 100)
0123     {
0124       fNumb = "0" + std::to_string(i);
0125     }
0126     else
0127     {
0128       fNumb = std::to_string(i);
0129     }
0130
0131     string truthFile = "outputData_" + fileModuleName + "_00" + fNumb + ".root";
0132     string KFPFile = "outputData_" + fileModuleName + "_KFP_00" + fNumb + ".root";
0133     truthFilePath = inDir + truthFile;
0134     KFPFilePath = inDir + KFPFile;
0135
0136     TFile *truthInput(0);
0137     truthInput = TFile::Open(truthFilePath);
0138     TTree *truthTree = (TTree*)truthInput->Get("QAG4SimulationTruthDecay");
0139     TFile *KFPInput(0);
0140     KFPInput = TFile::Open(KFPFilePath);
0141     TTree *recoTree = (TTree*)KFPInput->Get("DecayTree");
0142
0143     Float_t truthMotherPx, truthMotherPy, truthMotherPz, truthMotherPe;
0144     Float_t truthT1Px, truthT1Py, truthT1Pz, truthT1Pe;
0145     Float_t truthT2Px, truthT2Py, truthT2Pz, truthT2Pe;
0146     Float_t truthT1Pt, truthT2Pt;
0147     Float_t truthT1Eta, truthT2Eta;
0148     Float_t recoMotherPx, recoMotherPy, recoMotherPz, recoMotherPe;
0149     Float_t KFPTruthT1Px, KFPTruthT1Py, KFPTruthT1Pz, KFPTruthT1Pe;
0150     Float_t KFPTruthT2Px, KFPTruthT2Py, KFPTruthT2Pz, KFPTruthT2Pe;
0151     Float_t truthMotherPt, truthMotherEta;
0152     Int_t truthMotherBarcode;
0153     Float_t recoMotherMass;
0154
0155     truthTree->SetBranchAddress("mother_px", &truthMotherPx);
0156     truthTree->SetBranchAddress("mother_py", &truthMotherPy);
0157     truthTree->SetBranchAddress("mother_pz", &truthMotherPz);
0158     truthTree->SetBranchAddress("mother_pE", &truthMotherPe);
0159     truthTree->SetBranchAddress("mother_pT", &truthMotherPt);
0160     truthTree->SetBranchAddress("mother_eta", &truthMotherEta);
0161     truthTree->SetBranchAddress("mother_barcode", &truthMotherBarcode);
0162     truthTree->SetBranchAddress("track_1_px", &truthT1Px);
0163     truthTree->SetBranchAddress("track_1_py", &truthT1Py);
0164     truthTree->SetBranchAddress("track_1_pz", &truthT1Pz);
0165     truthTree->SetBranchAddress("track_1_pT", &truthT1Pt);
0166     truthTree->SetBranchAddress("track_2_px", &truthT2Px);
0167     truthTree->SetBranchAddress("track_2_py", &truthT2Py);
0168     truthTree->SetBranchAddress("track_2_pz", &truthT2Pz);
0169     truthTree->SetBranchAddress("track_2_pT", &truthT2Pt);
0170     truthTree->SetBranchAddress("track_1_eta", &truthT1Eta);
0171     truthTree->SetBranchAddress("track_2_eta", &truthT2Eta);
0172
0173     recoTree->SetBranchAddress("D0_px", &recoMotherPx);
0174     recoTree->SetBranchAddress("D0_py", &recoMotherPy);
0175     recoTree->SetBranchAddress("D0_pz", &recoMotherPz);
0176     recoTree->SetBranchAddress("D0_pE", &recoMotherPe);
0177     recoTree->SetBranchAddress("D0_mass", &recoMotherMass);
0178     recoTree->SetBranchAddress("track_1_true_px", &KFPTruthT1Px);
0179     recoTree->SetBranchAddress("track_1_true_py", &KFPTruthT1Py);
0180     recoTree->SetBranchAddress("track_1_true_pz", &KFPTruthT1Pz);
0181     recoTree->SetBranchAddress("track_2_true_px", &KFPTruthT2Px);
0182     recoTree->SetBranchAddress("track_2_true_py", &KFPTruthT2Py);
0183     recoTree->SetBranchAddress("track_2_true_pz", &KFPTruthT2Pz);
0184
0185     vector<Int_t> reconstructedBarcodes;
0186     vector<Int_t> missedBarcodes;
0187     vector<int> usedRecoEntries;
0188
0189     Float_t min_deltaPx, min_deltaPy, min_deltaPz, min_deltaPe;
0190     int minEntry;
0191
0192     for (int j = 0; j < truthTree->GetEntries(); ++j)
0193     {
0194       truthTree->GetEntry(j);
0195       totalTruthEvents += 1;
0196
0197       if (truthT1Pt > pTAccept && truthT2Pt > pTAccept) insideKinematic += 1;
0198       if (abs(truthT1Eta) < etaAccept && abs(truthT2Eta) < etaAccept) insideGeometric += 1;
0199       if (abs(truthT1Eta) < etaAccept && abs(truthT2Eta) < etaAccept && truthT1Pt > pTAccept && truthT2Pt > pTAccept) insideBoth += 1;
0200
0201       if (truthMotherPt < 1) truthTotal_01 += 1;
0202       if (truthMotherPt < 2 && truthMotherPt >= 1) truthTotal_12 += 1;
0203       if (truthMotherPt < 3 && truthMotherPt >= 2) truthTotal_23 += 1;
0204       if (truthMotherPt < 4 && truthMotherPt >= 3) truthTotal_34 += 1;
0205       if (truthMotherPt < 5 && truthMotherPt >= 4) truthTotal_45 += 1;
0206       if (truthMotherPt >= 5) truthTotal_5 += 1;
0207
0208       if (truthMotherEta < -1) truthTotal_n11 += 1;
0209       if (truthMotherEta < -1 && truthMotherEta >= -1.1) truthTotal_n11n1 += 1;
0210       if (truthMotherEta < -0.8 && truthMotherEta >= -1) truthTotal_n1n08 += 1;
0211       if (truthMotherEta < -0.6 && truthMotherEta >= -0.8) truthTotal_n08n06 += 1;
0212       if (truthMotherEta < -0.4 && truthMotherEta >= -0.6) truthTotal_n06n04 += 1;
0213       if (truthMotherEta < -0.2 && truthMotherEta >= -0.4) truthTotal_n04n02 += 1;
0214       if (truthMotherEta < 0.0 && truthMotherEta >= -0.2) truthTotal_n02n00 += 1;
0215       if (truthMotherEta < 0.2 && truthMotherEta >= 0.0) truthTotal_p02p00 += 1;
0216       if (truthMotherEta < 0.4 && truthMotherEta >= 0.2) truthTotal_p04p02 += 1;
0217       if (truthMotherEta < 0.6 && truthMotherEta >= 0.4) truthTotal_p06p04 += 1;
0218       if (truthMotherEta < 0.8 && truthMotherEta >= 0.6) truthTotal_p08p06 += 1;
0219       if (truthMotherEta < 1.0 && truthMotherEta >= 0.8) truthTotal_p1p08 += 1;
0220       if (truthMotherEta < 1.1 && truthMotherEta >= 1) truthTotal_p11p1 += 1;
0221       if (truthMotherEta >= 1) truthTotal_p11 += 1;
0222
0223       if (abs(truthT1Eta) < etaAccept && abs(truthT2Eta) < etaAccept && truthT1Pt > pTAccept && truthT2Pt > pTAccept)
0224       {
0225         if (truthMotherPt < 1) truthAccept_01 += 1;
0226         if (truthMotherPt < 2 && truthMotherPt >= 1) truthAccept_12 += 1;
0227         if (truthMotherPt < 3 && truthMotherPt >= 2) truthAccept_23 += 1;
0228         if (truthMotherPt < 4 && truthMotherPt >= 3) truthAccept_34 += 1;
0229         if (truthMotherPt < 5 && truthMotherPt >= 4) truthAccept_45 += 1;
0230         if (truthMotherPt >= 5) truthAccept_5 += 1;
0231
0232         if (truthMotherEta < -1) truthAccept_n11 += 1;
0233         if (truthMotherEta < -1 && truthMotherEta >= -1.1) truthAccept_n11n1 += 1;
0234         if (truthMotherEta < -0.8 && truthMotherEta >= -1) truthAccept_n1n08 += 1;
0235         if (truthMotherEta < -0.6 && truthMotherEta >= -0.8) truthAccept_n08n06 += 1;
0236         if (truthMotherEta < -0.4 && truthMotherEta >= -0.6) truthAccept_n06n04 += 1;
0237         if (truthMotherEta < -0.2 && truthMotherEta >= -0.4) truthAccept_n04n02 += 1;
0238         if (truthMotherEta < 0.0 && truthMotherEta >= -0.2) truthAccept_n02n00 += 1;
0239         if (truthMotherEta < 0.2 && truthMotherEta >= 0.0) truthAccept_p02p00 += 1;
0240         if (truthMotherEta < 0.4 && truthMotherEta >= 0.2) truthAccept_p04p02 += 1;
0241         if (truthMotherEta < 0.6 && truthMotherEta >= 0.4) truthAccept_p06p04 += 1;
0242         if (truthMotherEta < 0.8 && truthMotherEta >= 0.6) truthAccept_p08p06 += 1;
0243         if (truthMotherEta < 1.0 && truthMotherEta >= 0.8) truthAccept_p1p08 += 1;
0244         if (truthMotherEta < 1.1 && truthMotherEta >= 1) truthAccept_p11p1 += 1;
0245         if (truthMotherEta >= 1) truthAccept_p11 += 1;
0246       }
0247
0248       bool matchPx_t1r1;
0249       bool matchPy_t1r1;
0250       bool matchPz_t1r1;
0251       bool matchPx_t1r2;
0252       bool matchPy_t1r2;
0253       bool matchPz_t1r2;
0254       bool matchPx_t2r1;
0255       bool matchPy_t2r1;
0256       bool matchPz_t2r1;
0257       bool matchPx_t2r2;
0258       bool matchPy_t2r2;
0259       bool matchPz_t2r2;
0260
0261       bool matches = false;
0262
0263       for (int k = 0; k < recoTree->GetEntries(); ++k)
0264       {
0265         recoTree->GetEntry(k);
0266         matchPx_t1r1 = (KFPTruthT1Px == truthT1Px);
0267         matchPy_t1r1 = (KFPTruthT1Py == truthT1Py);
0268         matchPz_t1r1 = (KFPTruthT1Pz == truthT1Pz);
0269         matchPx_t1r2 = (KFPTruthT2Px == truthT1Px);
0270         matchPy_t1r2 = (KFPTruthT2Py == truthT1Py);
0271         matchPz_t1r2 = (KFPTruthT2Pz == truthT1Pz);
0272         matchPx_t2r1 = (KFPTruthT1Px == truthT2Px);
0273         matchPy_t2r1 = (KFPTruthT1Py == truthT2Py);
0274         matchPz_t2r1 = (KFPTruthT1Pz == truthT2Pz);
0275         matchPx_t2r2 = (KFPTruthT2Px == truthT2Px);
0276         matchPy_t2r2 = (KFPTruthT2Py == truthT2Py);
0277         matchPz_t2r2 = (KFPTruthT2Pz == truthT2Pz);
0278
0279         // check the truth information from KFParticle and QAG4SimulationTruthDecay to match a reconstructed candidate with its truth information
0280         // if match is found, that means a true D0->PiK event was reconstructed in KFParticle
0281         if (matchPx_t1r1 && matchPy_t1r1 && matchPz_t1r1)
0282         {
0283           if (matchPx_t2r2 && matchPy_t2r2 && matchPz_t2r2)
0284           {
0285             std::cout << "Got a match" << std::endl;
0286             matches = true;
0287             minEntry = k;
0288             break;
0289           }
0290         }
0291         else if (matchPx_t1r2 && matchPy_t1r2 && matchPz_t1r2)
0292         {
0293           if (matchPx_t2r1 && matchPy_t2r1 && matchPz_t2r1)
0294           {
0295             std::cout << "Got a match" << std::endl;
0296             matches = true;
0297             minEntry = k;
0298             break;
0299           }
0300         }
0301       }
0302       if (matches)
0303       {
0304         if (usedRecoEntries.size() > 0)
0305         {
0306           bool recoCandidateUsed = false;
0307           for (int ent : usedRecoEntries)
0308           {
0309             if (ent == minEntry)
0310         {
0311               std::cout << "Candidate has already been used, skipping this candidate" << std::endl;
0312               recoCandidateUsed = true;
0313         }
0314           }
0315           if (recoCandidateUsed) continue;
0316           else
0317           {
0318             usedRecoEntries.push_back(minEntry);
0319
0320             if (abs(truthT1Eta) < etaAccept && abs(truthT2Eta) < etaAccept && truthT1Pt > pTAccept && truthT2Pt > pTAccept)
0321             {
0322               totalReconstructed += 1;
0323
0324               if (truthMotherPt < 1) reconstructed_01 += 1;
0325               if (truthMotherPt < 2 && truthMotherPt >= 1) reconstructed_12 += 1;
0326               if (truthMotherPt < 3 && truthMotherPt >= 2) reconstructed_23 += 1;
0327               if (truthMotherPt < 4 && truthMotherPt >= 3) reconstructed_34 += 1;
0328               if (truthMotherPt < 5 && truthMotherPt >= 4) reconstructed_45 += 1;
0329               if (truthMotherPt >= 5) reconstructed_5 += 1;
0330
0331               if (truthMotherEta < -1.1) reconstructed_n11 += 1;
0332               if (truthMotherEta < -1 && truthMotherEta >= -1.1) reconstructed_n11n1 += 1;
0333               if (truthMotherEta < -0.8 && truthMotherEta >= -1) reconstructed_n1n08 += 1;
0334               if (truthMotherEta < -0.6 && truthMotherEta >= -0.8) reconstructed_n08n06 += 1;
0335               if (truthMotherEta < -0.4 && truthMotherEta >= -0.6) reconstructed_n06n04 += 1;
0336               if (truthMotherEta < -0.2 && truthMotherEta >= -0.4) reconstructed_n04n02 += 1;
0337               if (truthMotherEta < 0.0 && truthMotherEta >= -0.2) reconstructed_n02n00 += 1;
0338               if (truthMotherEta < 0.2 && truthMotherEta >= 0.0) reconstructed_p02p00 += 1;
0339               if (truthMotherEta < 0.4 && truthMotherEta >= 0.2) reconstructed_p04p02 += 1;
0340               if (truthMotherEta < 0.6 && truthMotherEta >= 0.4) reconstructed_p06p04 += 1;
0341               if (truthMotherEta < 0.8 && truthMotherEta >= 0.6) reconstructed_p08p06 += 1;
0342               if (truthMotherEta < 1.0 && truthMotherEta >= 0.8) reconstructed_p1p08 += 1;
0343               if (truthMotherEta < 1.1 && truthMotherEta >= 1) reconstructed_p11p1 += 1;
0344               if (truthMotherEta >= 1) reconstructed_p11 += 1;
0345             }
0346           }
0347         }
0348         else
0349         {
0350           usedRecoEntries.push_back(minEntry);
0351
0352           if (abs(truthT1Eta) < etaAccept && abs(truthT2Eta) < etaAccept && truthT1Pt > pTAccept && truthT2Pt > pTAccept)
0353           {
0354             totalReconstructed += 1;
0355
0356             if (truthMotherPt < 1) reconstructed_01 += 1;
0357             if (truthMotherPt < 2 && truthMotherPt >= 1) reconstructed_12 += 1;
0358             if (truthMotherPt < 3 && truthMotherPt >= 2) reconstructed_23 += 1;
0359             if (truthMotherPt < 4 && truthMotherPt >= 3) reconstructed_34 += 1;
0360             if (truthMotherPt < 5 && truthMotherPt >= 4) reconstructed_45 += 1;
0361             if (truthMotherPt >= 5) reconstructed_5 += 1;
0362
0363             if (truthMotherEta < -1.1) reconstructed_n11 += 1;
0364             if (truthMotherEta < -1 && truthMotherEta >= -1.1) reconstructed_n11n1 += 1;
0365             if (truthMotherEta < -0.8 && truthMotherEta >= -1) reconstructed_n1n08 += 1;
0366             if (truthMotherEta < -0.6 && truthMotherEta >= -0.8) reconstructed_n08n06 += 1;
0367             if (truthMotherEta < -0.4 && truthMotherEta >= -0.6) reconstructed_n06n04 += 1;
0368             if (truthMotherEta < -0.2 && truthMotherEta >= -0.4) reconstructed_n04n02 += 1;
0369             if (truthMotherEta < 0.0 && truthMotherEta >= -0.2) reconstructed_n02n00 += 1;
0370             if (truthMotherEta < 0.2 && truthMotherEta >= 0.0) reconstructed_p02p00 += 1;
0371             if (truthMotherEta < 0.4 && truthMotherEta >= 0.2) reconstructed_p04p02 += 1;
0372             if (truthMotherEta < 0.6 && truthMotherEta >= 0.4) reconstructed_p06p04 += 1;
0373             if (truthMotherEta < 0.8 && truthMotherEta >= 0.6) reconstructed_p08p06 += 1;
0374             if (truthMotherEta < 1.0 && truthMotherEta >= 0.8) reconstructed_p1p08 += 1;
0375             if (truthMotherEta < 1.1 && truthMotherEta >= 1) reconstructed_p11p1 += 1;
0376             if (truthMotherEta >= 1) reconstructed_p11 += 1;
0377           }
0378         }
0379       }
0380     }
0381   }
0382   std::cout << "truth mother pT efficiencies" << std::endl;
0383   std::cout << "0-1 GeV: " << reconstructed_01/truthAccept_01 << std::endl;
0384   std::cout << "1-2 GeV: " << reconstructed_12/truthAccept_12 << std::endl;
0385   std::cout << "2-3 GeV: " << reconstructed_23/truthAccept_23 << std::endl;
0386   std::cout << "3-4 GeV: " << reconstructed_34/truthAccept_34 << std::endl;
0387   std::cout << "4-5 GeV: " << reconstructed_45/truthAccept_45 << std::endl;
0388   std::cout << "5+ GeV: " << reconstructed_5/truthAccept_5 << std::endl;
0389
0390   std::cout << "" << std::endl;
0391   std::cout << "Number of truth events w/in truth track acceptance" << std::endl;
0392
0393   std::cout << "0-1 GeV: " << truthAccept_01 << std::endl;
0394   std::cout << "1-2 GeV: " << truthAccept_12 << std::endl;
0395   std::cout << "2-3 GeV: " << truthAccept_23 << std::endl;
0396   std::cout << "3-4 GeV: " << truthAccept_34 << std::endl;
0397   std::cout << "4-5 GeV: " << truthAccept_45 << std::endl;
0398   std::cout << "5+ GeV: " << truthAccept_5 << std::endl;
0399
0400   std::cout << "" << std::endl;
0401
0402   std::cout << "truth mother eta acceptance" << std::endl;
0403   std::cout << "eta < -1.1 : " << reconstructed_n11/truthAccept_n11 << std::endl;
0404   std::cout << "-1.1 <= eta < -1.0 : " << reconstructed_n11n1/truthAccept_n11n1 << std::endl;
0405   std::cout << "-1 <= eta < -0.8 : " << reconstructed_n1n08/truthAccept_n1n08 << std::endl;
0406   std::cout << "-0.8 <= eta < -0.6 : " << reconstructed_n08n06/truthAccept_n08n06 << std::endl;
0407   std::cout << "-0.6 <= eta < -0.4 : " << reconstructed_n06n04/truthAccept_n06n04 << std::endl;
0408   std::cout << "-0.4 <= eta < -0.2 : " << reconstructed_n04n02/truthAccept_n04n02 << std::endl;
0409   std::cout << "-0.2 <= eta < 0.0 : " << reconstructed_n02n00/truthAccept_n02n00 << std::endl;
0410   std::cout << "0.0 <= eta < 0.2 : " << reconstructed_p02p00/truthAccept_p02p00 << std::endl;
0411   std::cout << "0.2 <= eta < 0.4 : " << reconstructed_p04p02/truthAccept_p04p02 << std::endl;
0412   std::cout << "0.4 <= eta < 0.6 : " << reconstructed_p06p04/truthAccept_p06p04 << std::endl;
0413   std::cout << "0.6 <= eta < 0.8 : " << reconstructed_p08p06/truthAccept_p08p06 << std::endl;
0414   std::cout << "0.8 <= eta < 1.0 : " << reconstructed_p1p08/truthAccept_p1p08 << std::endl;
0415   std::cout << "1.0 <= eta < 1.1 : " << reconstructed_p11p1/truthAccept_p11p1 << std::endl;
0416   std::cout << "eta >= 1.1 : " << reconstructed_p11/truthAccept_p11 << std::endl;
0417
0418   std::cout << "" << std::endl;
0419   std::cout << "Number of truth events w/in truth track acceptance" << std::endl;
0420
0421   std::cout << "eta < -1.1 : " << truthAccept_n11 << std::endl;
0422   std::cout << "-1.1 <= eta < -1.0 : " << truthAccept_n11n1 << std::endl;
0423   std::cout << "-1 <= eta < -0.8 : " << truthAccept_n1n08 << std::endl;
0424   std::cout << "-0.8 <= eta < -0.6 : " << truthAccept_n08n06 << std::endl;
0425   std::cout << "-0.6 <= eta < -0.4 : " << truthAccept_n06n04 << std::endl;
0426   std::cout << "-0.4 <= eta < -0.2 : " << truthAccept_n04n02 << std::endl;
0427   std::cout << "-0.2 <= eta < 0.0 : " << truthAccept_n02n00 << std::endl;
0428   std::cout << "0.0 <= eta < 0.2 : " << truthAccept_p02p00 << std::endl;
0429   std::cout << "0.2 <= eta < 0.4 : " << truthAccept_p04p02 << std::endl;
0430   std::cout << "0.4 <= eta < 0.6 : " << truthAccept_p06p04 << std::endl;
0431   std::cout << "0.6 <= eta < 0.8 : " << truthAccept_p08p06 << std::endl;
0432   std::cout << "0.8 <= eta < 1.0 : " << truthAccept_p1p08 << std::endl;
0433   std::cout << "1.0 <= eta < 1.1 : " << truthAccept_p11p1 << std::endl;
0434   std::cout << "eta >= 1.1 : " << truthAccept_p11 << std::endl;
0435
0436   std::cout << "Total Number of Truth Events: " << totalTruthEvents << std::endl;
0437   std::cout << "Fraction of events w/ truth track pT > " << pTAccept << " GeV: " << insideKinematic/totalTruthEvents << std::endl;
0438   std::cout << "Fraction of events w/ truth track |eta| < " << etaAccept << ": " << insideGeometric/totalTruthEvents << std::endl;
0439   std::cout << "Fraction of events w/ both: " << insideBoth/totalTruthEvents << std::endl;
0440   std::cout << "Total Integrated Efficiency: " << totalReconstructed/insideBoth << std::endl;
0441
0442   const Int_t n_pT = 6;
0443   Double_t pT_bin[n_pT] = {0.5,1.5,2.5,3.5,4.5,5.5};
0444   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};
0445   Double_t pT_bin_widths[n_pT] = {0.5,0.5,0.5,0.5,0.5,0.5};
0446   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};
0447
0448   TCanvas *c1 = new TCanvas("D0_Efficiency_pT", "D0_Efficiency_pT", 1800, 900);
0449   TGraphErrors* D0_Efficiency_pT = new TGraphErrors(n_pT, pT_bin, pT_efficiency, pT_bin_widths, pT_efficiency_err);
0450   D0_Efficiency_pT->SetTitle("D^{0} Reconstruction Efficiency vs Truth Mother p_{T};Truth p_{T} [GeV/c];Reconstruction Efficiency");
0451   D0_Efficiency_pT->SetLineColor(kBlue);
0452   D0_Efficiency_pT->GetXaxis()->SetLimits(0.,6.0);
0453   D0_Efficiency_pT->SetMinimum(0.);
0454   D0_Efficiency_pT->SetMaximum(1.0);
0455   D0_Efficiency_pT->SetMarkerStyle(22);
0456   D0_Efficiency_pT->SetMarkerColor(kBlue);
0457   D0_Efficiency_pT->SetMarkerSize(2);
0458   D0_Efficiency_pT->SetLineWidth(2);
0459   D0_Efficiency_pT->Draw("ACP");
0460   auto legend = new TLegend(0.1,0.7,0.48,0.9);
0461   legend->AddEntry("","#it{#bf{sPHENIX}} Simulation","");
0462   legend->AddEntry("", "p+p #rightarrow D^{0} (#pi^{+}K^{-}) or #bar{D}^{0} (#pi^{-}K^{+})", "");
0463   legend->AddEntry("", "#sqrt{s}=200 GeV", "");
0464   legend->AddEntry("", "True Track pT > 0.2 GeV, |#eta| < 0.5", "");
0465   legend->SetMargin(0);
0466   legend->Draw();
0467   c1->SaveAs("D0_Efficiency_pT.png");
0468   c1->Close();
0469
0470   const Int_t n_eta = 14;
0471   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};
0472   Double_t eta_efficiency[n_eta] = {reconstructed_n11/truthAccept_n11, reconstructed_n11n1/truthAccept_n11n1, reconstructed_n1n08/truthAccept_n1n08, reconstructed_n08n06/truthAccept_n08n06,
0473                                    reconstructed_n06n04/truthAccept_n06n04, reconstructed_n04n02/truthAccept_n04n02, reconstructed_n02n00/truthAccept_n02n00, reconstructed_p02p00/truthAccept_p02p00,
0474                                    reconstructed_p04p02/truthAccept_p04p02, reconstructed_p06p04/truthAccept_p06p04, reconstructed_p08p06/truthAccept_p08p06, reconstructed_p1p08/truthAccept_p1p08,
0475                                    reconstructed_p11p1/truthAccept_p11p1, reconstructed_p11/truthAccept_p11};
0476   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};
0477   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,
0478                                        sqrt(truthAccept_n06n04)/truthAccept_n06n04,sqrt(truthAccept_n04n02)/truthAccept_n04n02,sqrt(truthAccept_n02n00)/truthAccept_n02n00,sqrt(truthAccept_p02p00)/truthAccept_p02p00,
0479                                        sqrt(truthAccept_p04p02)/truthAccept_p04p02,sqrt(truthAccept_p06p04)/truthAccept_p06p04,sqrt(truthAccept_p08p06)/truthAccept_p08p06,sqrt(truthAccept_p1p08)/truthAccept_p1p08,
0480                                        sqrt(truthAccept_p11p1)/truthAccept_p11p1,sqrt(truthAccept_p11)/truthAccept_p11};
0481
0482   TCanvas *c2 = new TCanvas("D0_Efficiency_eta", "D0_Efficiency_eta", 1800, 900);
0483   TGraphErrors* D0_Efficiency_eta = new TGraphErrors(n_eta, eta_bin, eta_efficiency, eta_bin_widths, eta_efficiency_err);
0484   D0_Efficiency_eta->SetTitle("D^{0} Reconstruction Efficiency vs Truth Mother #eta;Truth #eta;Reconstruction Efficiency");
0485   D0_Efficiency_eta->SetLineColor(kBlue);
0486   D0_Efficiency_eta->GetXaxis()->SetLimits(-1.25,1.25);
0487   D0_Efficiency_eta->SetMinimum(0.);
0488   D0_Efficiency_eta->SetMaximum(1.0);
0489   D0_Efficiency_eta->SetMarkerStyle(22);
0490   D0_Efficiency_eta->SetMarkerColor(kBlue);
0491   D0_Efficiency_eta->SetMarkerSize(2);
0492   D0_Efficiency_eta->SetLineWidth(2);
0493   D0_Efficiency_eta->Draw("ACP");
0494   auto legend2 = new TLegend(0.1,0.7,0.48,0.9);
0495   legend2->AddEntry("","#it{#bf{sPHENIX}} Simulation","");
0496   legend2->AddEntry("", "p+p #rightarrow D^{0} (#pi^{+}K^{-}) or #bar{D}^{0} (#pi^{-}K^{+})", "");
0497   legend2->AddEntry("", "#sqrt{s}=200 GeV", "");
0498   legend2->AddEntry("", "True Track pT > 0.2 GeV, |#eta| < 0.5", "");
0499   legend2->SetMargin(0);
0500   legend2->Draw();
0501   c2->SaveAs("D0_Efficiency_eta.png");
0502   c2->Close();
0503 }