JetsAndHF/jetReaderExamples/jetReader_TTreeReader.C

0001
0002 // Example script to read jet branch, find associated jet constituents, and confirm that constituents return the jet kinematics
0003 // Author: B. Page (bpage@bnl.gov)
0004 //
0005 // Usage: root -l -q jetReader_TTreeReader.C'("/path/to/eicrecon/output/file")'
0006
0007 void jetReader_TTreeReader(TString infile="root://dtn-eic.jlab.org//volatile/eic/EPIC/RECO/25.05.0/epic_craterlake/DIS/NC/18x275/minQ2=10/pythia8NCDIS_18x275_minQ2=10_beamEffects_xAngle=-0.025_hiDiv_5.18*.eicrecon.edm4eic.root") // 1950
0008 {
0009   // Input
0010   TChain *mychain = new TChain("events");
0011   mychain->Add(infile);
0012
0013   // Output
0014   TFile *ofile = TFile::Open("test.hist.root","RECREATE");
0015
0016   // TTreeReader
0017   TTreeReader tree_reader(mychain);
0018
0019   // Reco Jets
0020   TTreeReaderArray<int> recoType = {tree_reader, "ReconstructedChargedJets.type"};
0021   TTreeReaderArray<float> recoNRG = {tree_reader, "ReconstructedChargedJets.energy"};
0022   TTreeReaderArray<int> recoPDG = {tree_reader, "ReconstructedChargedJets.PDG"};
0023   TTreeReaderArray<float> recoMomX = {tree_reader, "ReconstructedChargedJets.momentum.x"};
0024   TTreeReaderArray<float> recoMomY = {tree_reader, "ReconstructedChargedJets.momentum.y"};
0025   TTreeReaderArray<float> recoMomZ = {tree_reader, "ReconstructedChargedJets.momentum.z"};
0026   TTreeReaderArray<float> recoM = {tree_reader, "ReconstructedChargedJets.mass"};
0027   TTreeReaderArray<unsigned int> partsBegin = {tree_reader, "ReconstructedChargedJets.particles_begin"};
0028   TTreeReaderArray<unsigned int> partsEnd = {tree_reader, "ReconstructedChargedJets.particles_end"};
0029
0030   TTreeReaderArray<int> recoPartIndex = {tree_reader, "_ReconstructedChargedJets_particles.index"};
0031
0032   // Reconstructed Particles
0033   TTreeReaderArray<float> recoPartMomX = {tree_reader, "ReconstructedChargedParticles.momentum.x"};
0034   TTreeReaderArray<float> recoPartMomY = {tree_reader, "ReconstructedChargedParticles.momentum.y"};
0035   TTreeReaderArray<float> recoPartMomZ = {tree_reader, "ReconstructedChargedParticles.momentum.z"};
0036   TTreeReaderArray<float> recoPartM = {tree_reader, "ReconstructedChargedParticles.mass"};
0037   TTreeReaderArray<int> recoPartPDG = {tree_reader, "ReconstructedChargedParticles.PDG"};
0038   TTreeReaderArray<float> recoPartNRG = {tree_reader, "ReconstructedChargedParticles.energy"};
0039
0040   // Uncomment the following two lines if using an eic-shell version older than 25.12.0-stable
0041   // Refer https://chat.epic-eic.org/main/pl/n9mbqtf4fiyptjz9q13f7m8xne
0042
0043   //TTreeReaderArray<unsigned int> recoPartAssocRec = {tree_reader, "ReconstructedChargedParticleAssociations.recID"}; // Reco <-> MCParticle
0044   //TTreeReaderArray<unsigned int> recoPartAssocSim = {tree_reader, "ReconstructedChargedParticleAssociations.simID"};
0045
0046   // updated code after eic-shell --version 25.12.0-stable
0047   TTreeReaderArray<int> recoPartAssocRec = {tree_reader, "_ReconstructedChargedParticleAssociations_rec.index"}; // Reco <-> MCParticle
0048   TTreeReaderArray<int> recoPartAssocSim = {tree_reader, "_ReconstructedChargedParticleAssociations_sim.index"};
0049
0050   TTreeReaderArray<float> recoPartAssocWeight = {tree_reader, "ReconstructedChargedParticleAssociations.weight"};
0051
0052   // Generated Jets
0053   TTreeReaderArray<int> genType = {tree_reader, "GeneratedChargedJets.type"};
0054   TTreeReaderArray<float> genNRG = {tree_reader, "GeneratedChargedJets.energy"};
0055   TTreeReaderArray<int> genPDG = {tree_reader, "GeneratedChargedJets.PDG"};
0056   TTreeReaderArray<float> genMomX = {tree_reader, "GeneratedChargedJets.momentum.x"};
0057   TTreeReaderArray<float> genMomY = {tree_reader, "GeneratedChargedJets.momentum.y"};
0058   TTreeReaderArray<float> genMomZ = {tree_reader, "GeneratedChargedJets.momentum.z"};
0059   TTreeReaderArray<float> genM = {tree_reader, "GeneratedChargedJets.mass"};
0060   TTreeReaderArray<unsigned int> genPartsBegin = {tree_reader, "GeneratedChargedJets.particles_begin"};
0061   TTreeReaderArray<unsigned int> genPartsEnd = {tree_reader, "GeneratedChargedJets.particles_end"};
0062
0063   TTreeReaderArray<int> genPartIndex = {tree_reader, "_GeneratedChargedJets_particles.index"};
0064
0065   // MC
0066   TTreeReaderArray<float> mcMomX = {tree_reader, "GeneratedParticles.momentum.x"};
0067   TTreeReaderArray<float> mcMomY = {tree_reader, "GeneratedParticles.momentum.y"};
0068   TTreeReaderArray<float> mcMomZ = {tree_reader, "GeneratedParticles.momentum.z"};
0069   TTreeReaderArray<float> mcM = {tree_reader, "GeneratedParticles.mass"};
0070   TTreeReaderArray<float> mcE = {tree_reader, "GeneratedParticles.energy"};
0071   TTreeReaderArray<int> pdg = {tree_reader, "GeneratedParticles.PDG"};
0072
0073   TTreeReaderArray<int> mcGenStat = {tree_reader, "MCParticles.generatorStatus"};
0074   TTreeReaderArray<double> mcMomXPart = {tree_reader, "MCParticles.momentum.x"};
0075   TTreeReaderArray<double> mcMomYPart = {tree_reader, "MCParticles.momentum.y"};
0076   TTreeReaderArray<double> mcMomZPart = {tree_reader, "MCParticles.momentum.z"};
0077   TTreeReaderArray<double> mcMPart = {tree_reader, "MCParticles.mass"};
0078   TTreeReaderArray<int> pdgMCPart = {tree_reader, "MCParticles.PDG"};
0079
0080
0081   // Define Histograms
0082   // Reco
0083   TH1D *numRecoJetsEventHist = new TH1D("numRecoJetsEvent","",20,0.,20.);
0084   TH1D *numRecoJetsNoElecEventHist = new TH1D("numRecoJetsNoElecEvent","",20,0.,20.);
0085
0086   TH1D *recoJetEHist = new TH1D("recoJetE","",300,0.,300.);
0087   TH2D *recoJetEvsEtaHist = new TH2D("recoJetEvsEta","",100,-5.,5.,300,0.,300.);
0088   TH2D *recoJetPhiVsEtaHist = new TH2D("recoJetPhiVsEta","",100,-5.,5.,100,-TMath::Pi(),TMath::Pi());
0089
0090   TH1D *recoJetENoElecHist = new TH1D("recoJetENoElec","",300,0.,300.);
0091   TH2D *recoJetEvsEtaNoElecHist = new TH2D("recoJetEvsEtaNoElec","",100,-5.,5.,300,0.,300.);
0092   TH2D *recoJetPhiVsEtaNoElecHist = new TH2D("recoJetPhiVsEtaNoElec","",100,-5.,5.,100,-TMath::Pi(),TMath::Pi());
0093
0094   TH1D *recoJetEElecHist = new TH1D("recoJetEElec","",300,0.,300.);
0095   TH2D *recoJetEvsEtaElecHist = new TH2D("recoJetEvsEtaElec","",100,-5.,5.,300,0.,300.);
0096   TH2D *recoJetPhiVsEtaElecHist = new TH2D("recoJetPhiVsEtaElec","",100,-5.,5.,100,-TMath::Pi(),TMath::Pi());
0097
0098   TH1D *constituentRecoPHist = new TH1D("constituentRecoP","",500,0.,50.);
0099   TH2D *constituentRecoPVsEtaHist = new TH2D("constituentRecoPVsEta","",100,-5.,5.,500,0.,50.);
0100   TH2D *constituentRecoPhiVsEtaHist = new TH2D("constituentRecoPhiVsEta","",100,-5.,5.,100,-TMath::Pi(),TMath::Pi());
0101
0102   TH1D *constituentRecoTrueDeltaRHist = new TH1D("constituentRecoTrueDeltaR","",5000,0.,5.);
0103
0104   TH2D *constituentRecoVsTruePHist = new TH2D("constituentRecoVsTrueP","",500,0.,50.,500,0.,50.);
0105   TH2D *constituentRecoVsTrueEtaHist = new TH2D("constituentRecoVsTrueEta","",100,-5.,5.,100,-5.,5.);
0106   TH2D *constituentRecoVsTruePhiHist = new TH2D("constituentRecoVsTruePhi","",100,-TMath::Pi(),TMath::Pi(),100,-TMath::Pi(),TMath::Pi());
0107
0108   TH1D *constituentPResHist = new TH1D("constituentPRes","",2000,-10.,10.);
0109
0110   TH2D *constituentRecoPVsEtaNoPIDHist = new TH2D("constituentRecoPVsEtaNoPID","",100,-5.,5.,500,0.,50.);
0111   TH2D *constituentAssocPartPVsEtaNoPIDHist = new TH2D("constituentAssocPartPVsEtaNoPID","",100,-5.,5.,500,0.,50.);
0112
0113   TH1D *constituentPDGElecHist = new TH1D("constituentPDGElec","",3000,0.,3000.);
0114
0115   TH1D *constituentPResElecHist = new TH1D("constituentPResElec","",2000,-10.,10.);
0116   TH2D *constituentRecoPVsEtaElecHist = new TH2D("constituentRecoPVsEtaElec","",100,-5.,5.,500,0.,50.);
0117   TH2D *constituentAssocPartPVsEtaElecHist = new TH2D("constituentAssocPartPVsEtaElec","",100,-5.,5.,500,0.,50.);
0118
0119   TH1D *constituentPResNoPIDElecHist = new TH1D("constituentPResNoPIDElec","",2000,-10.,10.);
0120   TH2D *constituentRecoPVsEtaNoPIDElecHist = new TH2D("constituentRecoPVsEtaNoPIDElec","",100,-5.,5.,500,0.,50.);
0121   TH2D *constituentAssocPartPVsEtaNoPIDElecHist = new TH2D("constituentAssocPartPVsEtaNoPIDElec","",100,-5.,5.,500,0.,50.);
0122
0123   TH1D *constituentPResBadPIDElecHist = new TH1D("constituentPResBadPIDElec","",2000,-10.,10.);
0124   TH2D *constituentRecoPVsEtaBadPIDElecHist = new TH2D("constituentRecoPVsEtaBadPIDElec","",100,-5.,5.,500,0.,50.);
0125   TH2D *constituentAssocPartPVsEtaBadPIDElecHist = new TH2D("constituentAssocPartPVsEtaBadPIDElec","",100,-5.,5.,500,0.,50.);
0126
0127   TH1D *constituentPDGPionHist = new TH1D("constituentPDGPion","",3000,0.,3000.);
0128
0129   TH1D *constituentPResPionHist = new TH1D("constituentPResPion","",2000,-10.,10.);
0130   TH2D *constituentRecoPVsEtaPionHist = new TH2D("constituentRecoPVsEtaPion","",100,-5.,5.,500,0.,50.);
0131   TH2D *constituentAssocPartPVsEtaPionHist = new TH2D("constituentAssocPartPVsEtaPion","",100,-5.,5.,500,0.,50.);
0132
0133   TH1D *constituentPResNoPIDPionHist = new TH1D("constituentPResNoPIDPion","",2000,-10.,10.);
0134   TH2D *constituentRecoPVsEtaNoPIDPionHist = new TH2D("constituentRecoPVsEtaNoPIDPion","",100,-5.,5.,500,0.,50.);
0135   TH2D *constituentAssocPartPVsEtaNoPIDPionHist = new TH2D("constituentAssocPartPVsEtaNoPIDPion","",100,-5.,5.,500,0.,50.);
0136
0137   TH1D *constituentPResBadPIDPionHist = new TH1D("constituentPResBadPIDPion","",2000,-10.,10.);
0138   TH2D *constituentRecoPVsEtaBadPIDPionHist = new TH2D("constituentRecoPVsEtaBadPIDPion","",100,-5.,5.,500,0.,50.);
0139   TH2D *constituentAssocPartPVsEtaBadPIDPionHist = new TH2D("constituentAssocPartPVsEtaBadPIDPion","",100,-5.,5.,500,0.,50.);
0140
0141   TH1D *constituentPResKaonHist = new TH1D("constituentPResKaon","",2000,-10.,10.);
0142   TH2D *constituentRecoPVsEtaKaonHist = new TH2D("constituentRecoPVsEtaKaon","",100,-5.,5.,500,0.,50.);
0143   TH2D *constituentAssocPartPVsEtaKaonHist = new TH2D("constituentAssocPartPVsEtaKaon","",100,-5.,5.,500,0.,50.);
0144
0145   TH1D *constituentPResNoPIDKaonHist = new TH1D("constituentPResNoPIDKaon","",2000,-10.,10.);
0146   TH2D *constituentRecoPVsEtaNoPIDKaonHist = new TH2D("constituentRecoPVsEtaNoPIDKaon","",100,-5.,5.,500,0.,50.);
0147   TH2D *constituentAssocPartPVsEtaNoPIDKaonHist = new TH2D("constituentAssocPartPVsEtaNoPIDKaon","",100,-5.,5.,500,0.,50.);
0148
0149   TH1D *constituentPResBadPIDKaonHist = new TH1D("constituentPResBadPIDKaon","",2000,-10.,10.);
0150   TH2D *constituentRecoPVsEtaBadPIDKaonHist = new TH2D("constituentRecoPVsEtaBadPIDKaon","",100,-5.,5.,500,0.,50.);
0151   TH2D *constituentAssocPartPVsEtaBadPIDKaonHist = new TH2D("constituentAssocPartPVsEtaBadPIDKaon","",100,-5.,5.,500,0.,50.);
0152
0153   TH1D *constituentPResProtonHist = new TH1D("constituentPResProton","",2000,-10.,10.);
0154   TH2D *constituentRecoPVsEtaProtonHist = new TH2D("constituentRecoPVsEtaProton","",100,-5.,5.,500,0.,50.);
0155   TH2D *constituentAssocPartPVsEtaProtonHist = new TH2D("constituentAssocPartPVsEtaProton","",100,-5.,5.,500,0.,50.);
0156
0157   TH1D *constituentPResNoPIDProtonHist = new TH1D("constituentPResNoPIDProton","",2000,-10.,10.);
0158   TH2D *constituentRecoPVsEtaNoPIDProtonHist = new TH2D("constituentRecoPVsEtaNoPIDProton","",100,-5.,5.,500,0.,50.);
0159   TH2D *constituentAssocPartPVsEtaNoPIDProtonHist = new TH2D("constituentAssocPartPVsEtaNoPIDProton","",100,-5.,5.,500,0.,50.);
0160
0161   TH1D *constituentPResBadPIDProtonHist = new TH1D("constituentPResBadPIDProton","",2000,-10.,10.);
0162   TH2D *constituentRecoPVsEtaBadPIDProtonHist = new TH2D("constituentRecoPVsEtaBadPIDProton","",100,-5.,5.,500,0.,50.);
0163   TH2D *constituentAssocPartPVsEtaBadPIDProtonHist = new TH2D("constituentAssocPartPVsEtaBadPIDProton","",100,-5.,5.,500,0.,50.);
0164
0165   TH1D *numRecoJetPartsHist = new TH1D("numRecoJetParts","",20,0.,20.);
0166   TH2D *recoJetEvsPartESumHist = new TH2D("recoJetEvsPartESum","",3000,0.,300.,3000,0.,300.);
0167   TH1D *recoJetEDiffHist = new TH1D("recoJetEDiff","",20000,-0.001,0.001);
0168
0169   TH2D *recoJetEvsEtaBadHist = new TH2D("recoJetEvsEtaBad","",100,-5.,5.,300,0.,300.);
0170   TH2D *recoJetPhiVsEtaBadHist = new TH2D("recoJetPhiVsEtaBad","",100,-5.,5.,100,-TMath::Pi(),TMath::Pi());
0171
0172   // Gen
0173   TH1D *numGenJetsEventHist = new TH1D("numGenJetsEvent","",20,0.,20.);
0174   TH1D *numGenJetsNoElecEventHist = new TH1D("numGenJetsNoElecEvent","",20,0.,20.);
0175
0176   TH1D *genJetEHist = new TH1D("genJetE","",300,0.,300.);
0177   TH2D *genJetEvsEtaHist = new TH2D("genJetEvsEta","",100,-5.,5.,300,0.,300.);
0178   TH2D *genJetPhiVsEtaHist = new TH2D("genJetPhiVsEta","",100,-5.,5.,100,-TMath::Pi(),TMath::Pi());
0179
0180   TH1D *genJetENoElecHist = new TH1D("genJetENoElec","",300,0.,300.);
0181   TH2D *genJetEvsEtaNoElecHist = new TH2D("genJetEvsEtaNoElec","",100,-5.,5.,300,0.,300.);
0182   TH2D *genJetPhiVsEtaNoElecHist = new TH2D("genJetPhiVsEtaNoElec","",100,-5.,5.,100,-TMath::Pi(),TMath::Pi());
0183
0184   TH1D *genJetEElecHist = new TH1D("genJetEElec","",300,0.,300.);
0185   TH2D *genJetEvsEtaElecHist = new TH2D("genJetEvsEtaElec","",100,-5.,5.,300,0.,300.);
0186   TH2D *genJetPhiVsEtaElecHist = new TH2D("genJetPhiVsEtaElec","",100,-5.,5.,100,-TMath::Pi(),TMath::Pi());
0187
0188   TH1D *numGenJetPartsHist = new TH1D("numGenJetParts","",20,0.,20.);
0189   TH2D *genJetPartEvsEtaHist = new TH2D("genJetPartEvsEta","",100,-5.,5.,300,0.,300.);
0190   TH2D *genJetPartPhiVsEtaHist = new TH2D("genJetPartPhiVsEta","",100,-5.,5.,100,-TMath::Pi(),TMath::Pi());
0191
0192   TH2D *genJetEvsPartESumHist = new TH2D("genJetEvsPartESum","",3000,0.,300.,3000,0.,300.);
0193   TH1D *genJetEDiffHist = new TH1D("genJetEDiff","",20000,-0.001,0.001);
0194
0195   TH2D *genJetEvsEtaBadHist = new TH2D("genJetEvsEtaBad","",100,-5.,5.,300,0.,300.);
0196   TH2D *genJetPhiVsEtaBadHist = new TH2D("genJetPhiVsEtaBad","",100,-5.,5.,100,-TMath::Pi(),TMath::Pi());
0197
0198
0199   // Loop Through Events
0200   int NEVENTS = 0;
0201   while(tree_reader.Next()) {
0202
0203     if(NEVENTS%10000 == 0) cout << "Events Processed: " << NEVENTS << endl;
0204
0205     // Analyze Reonstructed Jets
0206     int numRecoJetsNoElec = 0;
0207     numRecoJetsEventHist->Fill(recoType.GetSize());
0208     for(unsigned int i=0; i<recoType.GetSize(); i++)
0209       {
0210     TVector3 jetMom(recoMomX[i],recoMomY[i],recoMomZ[i]);
0211
0212     if(TMath::Abs(jetMom.PseudoRapidity()) > 2.5)
0213       continue;
0214
0215     recoJetEHist->Fill(recoNRG[i]);
0216     recoJetEvsEtaHist->Fill(jetMom.PseudoRapidity(),recoNRG[i]);
0217     recoJetPhiVsEtaHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
0218
0219     // Check if Jet Contains an Electron - Use Particle Matching to Find True PID
0220     bool noElectron = true;
0221     for(unsigned int m=partsBegin[i]; m<partsEnd[i]; m++) // Loop over jet constituents
0222       {
0223         int elecIndex = -1;
0224         double elecIndexWeight = -1.0;
0225         int chargePartIndex = recoPartIndex[m]; // ReconstructedChargedParticle Index for m'th Jet Component
0226         for(unsigned int n=0; n<recoPartAssocRec.GetSize(); n++) // Loop Over All ReconstructedChargedParticleAssociations
0227           {
0228         if(recoPartAssocRec[n] == chargePartIndex) // Select Entry Matching the ReconstructedChargedParticle Index
0229           {
0230             if(recoPartAssocWeight[n] > elecIndexWeight) // Find Particle with Greatest Weight = Contributed Most Hits to Track
0231               {
0232             elecIndex = recoPartAssocSim[n]; // Get Index of MCParticle Associated with ReconstructedChargedParticle
0233             elecIndexWeight = recoPartAssocWeight[n];
0234               }
0235           }
0236           }
0237
0238         if(pdgMCPart[elecIndex] == 11) // Test if Matched Particle is an Electron
0239           noElectron = false;
0240       }
0241
0242     if(noElectron)
0243       {
0244         recoJetENoElecHist->Fill(recoNRG[i]);
0245         recoJetEvsEtaNoElecHist->Fill(jetMom.PseudoRapidity(),recoNRG[i]);
0246         recoJetPhiVsEtaNoElecHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
0247
0248         numRecoJetsNoElec++;
0249       }
0250
0251     if(!noElectron)
0252       {
0253         recoJetEElecHist->Fill(recoNRG[i]);
0254         recoJetEvsEtaElecHist->Fill(jetMom.PseudoRapidity(),recoNRG[i]);
0255         recoJetPhiVsEtaElecHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
0256       }
0257
0258
0259     // Look at Constituents
0260     double esum = 0.0;
0261     double eRes = -999.;
0262     for(unsigned int j=partsBegin[i]; j<partsEnd[i]; j++)
0263       {
0264         // partsbegin and partsEnd specify the entries from _ReconstructedChargedJets_particles.index that make up the jet
0265         // _ReconstructedChargedJets_particles.index stores the ReconstructedChargedParticles index of the jet constituent
0266
0267         TVector3 recoPartMom(recoPartMomX[recoPartIndex[j]],recoPartMomY[recoPartIndex[j]],recoPartMomZ[recoPartIndex[j]]);
0268         double mM = recoPartM[recoPartIndex[j]];
0269         double mE = recoPartNRG[recoPartIndex[j]];
0270         int mPDG = recoPartPDG[recoPartIndex[j]];
0271
0272         esum += mE;
0273
0274         constituentRecoPHist->Fill(recoPartMom.Mag());
0275         constituentRecoPVsEtaHist->Fill(recoPartMom.PseudoRapidity(),recoPartMom.Mag());
0276         constituentRecoPhiVsEtaHist->Fill(recoPartMom.PseudoRapidity(),recoPartMom.Phi());
0277
0278         // Find Associated MC Particle (Association with Largest Weight)
0279         int simIndex = -1;
0280         double simIndexWeight = -1.0;
0281         int chargePartIndex = recoPartIndex[j];
0282         for(unsigned int n=0; n<recoPartAssocRec.GetSize(); n++)
0283           {
0284         if(recoPartAssocRec[n] == chargePartIndex)
0285           {
0286             if(recoPartAssocWeight[n] > simIndexWeight)
0287               {
0288             simIndex = recoPartAssocSim[n];
0289             simIndexWeight = recoPartAssocWeight[n];
0290               }
0291           }
0292           }
0293
0294         // Define Matching Truth Particle
0295         TVector3 genPartMom(mcMomXPart[simIndex],mcMomYPart[simIndex],mcMomZPart[simIndex]);
0296         double mPartM = mcMPart[simIndex];
0297         int mPartPDG = pdgMCPart[simIndex];
0298         int mPartGenStat = mcGenStat[simIndex];
0299
0300         double momRes = (recoPartMom.Mag() - genPartMom.Mag())/genPartMom.Mag();
0301         constituentPResHist->Fill(momRes);
0302
0303         double dEta = recoPartMom.PseudoRapidity() - genPartMom.PseudoRapidity();
0304         double dPhi = TVector2::Phi_mpi_pi(recoPartMom.Phi() - genPartMom.Phi());
0305         double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
0306
0307         constituentRecoTrueDeltaRHist->Fill(dR);
0308
0309         constituentRecoVsTruePHist->Fill(genPartMom.Mag(),recoPartMom.Mag());
0310         constituentRecoVsTrueEtaHist->Fill(genPartMom.PseudoRapidity(),recoPartMom.PseudoRapidity());
0311         constituentRecoVsTruePhiHist->Fill(genPartMom.Phi(),recoPartMom.Phi());
0312
0313         // Look at PID
0314         if(mPDG == 0)
0315           {
0316         constituentRecoPVsEtaNoPIDHist->Fill(recoPartMom.PseudoRapidity(),recoPartMom.Mag());
0317         constituentAssocPartPVsEtaNoPIDHist->Fill(genPartMom.PseudoRapidity(),genPartMom.Mag());
0318           }
0319
0320         if(TMath::Abs(mPartPDG) == 11)
0321           {
0322         constituentPDGElecHist->Fill(TMath::Abs(mPDG));
0323
0324         constituentPResElecHist->Fill(momRes);
0325         constituentRecoPVsEtaElecHist->Fill(recoPartMom.PseudoRapidity(),recoPartMom.Mag());
0326         constituentAssocPartPVsEtaElecHist->Fill(genPartMom.PseudoRapidity(),genPartMom.Mag());
0327
0328         if(mPDG == 0)
0329           {
0330             constituentPResNoPIDElecHist->Fill(momRes);
0331             constituentRecoPVsEtaNoPIDElecHist->Fill(recoPartMom.PseudoRapidity(),recoPartMom.Mag());
0332             constituentAssocPartPVsEtaNoPIDElecHist->Fill(genPartMom.PseudoRapidity(),genPartMom.Mag());
0333           }
0334
0335         if(mPDG != mPartPDG)
0336           {
0337             constituentPResBadPIDElecHist->Fill(momRes);
0338             constituentRecoPVsEtaBadPIDElecHist->Fill(recoPartMom.PseudoRapidity(),recoPartMom.Mag());
0339             constituentAssocPartPVsEtaBadPIDElecHist->Fill(genPartMom.PseudoRapidity(),genPartMom.Mag());
0340           }
0341           }
0342         if(TMath::Abs(mPartPDG) == 211)
0343           {
0344         constituentPDGPionHist->Fill(TMath::Abs(mPDG));
0345
0346         constituentPResPionHist->Fill(momRes);
0347         constituentRecoPVsEtaPionHist->Fill(recoPartMom.PseudoRapidity(),recoPartMom.Mag());
0348         constituentAssocPartPVsEtaPionHist->Fill(genPartMom.PseudoRapidity(),genPartMom.Mag());
0349
0350         if(mPDG == 0)
0351           {
0352             constituentPResNoPIDPionHist->Fill(momRes);
0353             constituentRecoPVsEtaNoPIDPionHist->Fill(recoPartMom.PseudoRapidity(),recoPartMom.Mag());
0354             constituentAssocPartPVsEtaNoPIDPionHist->Fill(genPartMom.PseudoRapidity(),genPartMom.Mag());
0355           }
0356
0357         if(mPDG != mPartPDG)
0358           {
0359             constituentPResBadPIDPionHist->Fill(momRes);
0360             constituentRecoPVsEtaBadPIDPionHist->Fill(recoPartMom.PseudoRapidity(),recoPartMom.Mag());
0361             constituentAssocPartPVsEtaBadPIDPionHist->Fill(genPartMom.PseudoRapidity(),genPartMom.Mag());
0362           }
0363           }
0364         if(TMath::Abs(mPartPDG) == 321)
0365           {
0366         constituentPResKaonHist->Fill(momRes);
0367         constituentRecoPVsEtaKaonHist->Fill(recoPartMom.PseudoRapidity(),recoPartMom.Mag());
0368         constituentAssocPartPVsEtaKaonHist->Fill(genPartMom.PseudoRapidity(),genPartMom.Mag());
0369
0370         if(mPDG == 0)
0371           {
0372             constituentPResNoPIDKaonHist->Fill(momRes);
0373             constituentRecoPVsEtaNoPIDKaonHist->Fill(recoPartMom.PseudoRapidity(),recoPartMom.Mag());
0374             constituentAssocPartPVsEtaNoPIDKaonHist->Fill(genPartMom.PseudoRapidity(),genPartMom.Mag());
0375           }
0376
0377         if(mPDG != mPartPDG)
0378           {
0379             constituentPResBadPIDKaonHist->Fill(momRes);
0380             constituentRecoPVsEtaBadPIDKaonHist->Fill(recoPartMom.PseudoRapidity(),recoPartMom.Mag());
0381             constituentAssocPartPVsEtaBadPIDKaonHist->Fill(genPartMom.PseudoRapidity(),genPartMom.Mag());
0382           }
0383           }
0384         if(TMath::Abs(mPartPDG) == 2212)
0385           {
0386         constituentPResProtonHist->Fill(momRes);
0387         constituentRecoPVsEtaProtonHist->Fill(recoPartMom.PseudoRapidity(),recoPartMom.Mag());
0388         constituentAssocPartPVsEtaProtonHist->Fill(genPartMom.PseudoRapidity(),genPartMom.Mag());
0389
0390         if(mPDG == 0)
0391           {
0392             constituentPResNoPIDProtonHist->Fill(momRes);
0393             constituentRecoPVsEtaNoPIDProtonHist->Fill(recoPartMom.PseudoRapidity(),recoPartMom.Mag());
0394             constituentAssocPartPVsEtaNoPIDProtonHist->Fill(genPartMom.PseudoRapidity(),genPartMom.Mag());
0395           }
0396
0397         if(mPDG != mPartPDG)
0398           {
0399             constituentPResBadPIDProtonHist->Fill(momRes);
0400             constituentRecoPVsEtaBadPIDProtonHist->Fill(recoPartMom.PseudoRapidity(),recoPartMom.Mag());
0401             constituentAssocPartPVsEtaBadPIDProtonHist->Fill(genPartMom.PseudoRapidity(),genPartMom.Mag());
0402           }
0403           }
0404       }
0405
0406     numRecoJetsNoElecEventHist->Fill(numRecoJetsNoElec);
0407     numRecoJetPartsHist->Fill(partsEnd[i] - partsBegin[i]);
0408     recoJetEvsPartESumHist->Fill(recoNRG[i],esum);
0409     recoJetEDiffHist->Fill(recoNRG[i]-esum);
0410
0411     if(TMath::Abs(esum - recoNRG[i]) > 0.000001)
0412       {
0413         recoJetEvsEtaBadHist->Fill(jetMom.PseudoRapidity(),recoNRG[i]);
0414         recoJetPhiVsEtaBadHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
0415       }
0416       }
0417
0418
0419     // Analyze Generated Jets
0420     int numGenJetsNoElec = 0;
0421     numGenJetsEventHist->Fill(genType.GetSize());
0422     for(unsigned int i=0; i<genType.GetSize(); i++)
0423       {
0424     TVector3 jetMom(genMomX[i],genMomY[i],genMomZ[i]);
0425
0426     if(TMath::Abs(jetMom.PseudoRapidity()) > 2.5)
0427       continue;
0428
0429     genJetEHist->Fill(genNRG[i]);
0430     genJetEvsEtaHist->Fill(jetMom.PseudoRapidity(),genNRG[i]);
0431     genJetPhiVsEtaHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
0432
0433     double esumG = 0.0;
0434     bool noGenElectron = true;
0435     for(unsigned int j=genPartsBegin[i]; j<genPartsEnd[i]; j++)
0436       {
0437         double mX = mcMomX[genPartIndex[j]];
0438         double mY = mcMomY[genPartIndex[j]];
0439         double mZ = mcMomZ[genPartIndex[j]];
0440         double mM = mcM[genPartIndex[j]];
0441         int mPDG = pdg[genPartIndex[j]];
0442
0443         double tmpE = TMath::Sqrt(mX*mX + mY*mY + mZ*mZ + mM*mM);
0444
0445         esumG += tmpE;
0446
0447         if(mPDG == 11)
0448           noGenElectron = false;
0449
0450         TVector3 partMom(mX,mY,mZ);
0451
0452         genJetPartEvsEtaHist->Fill(partMom.PseudoRapidity(),tmpE);
0453         genJetPartPhiVsEtaHist->Fill(partMom.PseudoRapidity(),partMom.Phi());
0454       }
0455
0456     if(noGenElectron)
0457       {
0458         genJetENoElecHist->Fill(genNRG[i]);
0459         genJetEvsEtaNoElecHist->Fill(jetMom.PseudoRapidity(),genNRG[i]);
0460         genJetPhiVsEtaNoElecHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
0461
0462         numGenJetsNoElec++;
0463       }
0464
0465     if(!noGenElectron)
0466       {
0467         genJetEElecHist->Fill(genNRG[i]);
0468         genJetEvsEtaElecHist->Fill(jetMom.PseudoRapidity(),genNRG[i]);
0469         genJetPhiVsEtaElecHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
0470       }
0471
0472     numGenJetsNoElecEventHist->Fill(numGenJetsNoElec);
0473     numGenJetPartsHist->Fill(genPartsEnd[i] - genPartsBegin[i]);
0474     genJetEvsPartESumHist->Fill(genNRG[i],esumG);
0475     genJetEDiffHist->Fill(genNRG[i]-esumG);
0476
0477     if(TMath::Abs(esumG - genNRG[i]) > 0.000001)
0478       {
0479         genJetEvsEtaBadHist->Fill(jetMom.PseudoRapidity(),genNRG[i]);
0480         genJetPhiVsEtaBadHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
0481       }
0482       }
0483
0484     NEVENTS++;
0485   }
0486
0487   ofile->Write();
0488   ofile->Close();
0489
0490 }