Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-07-11 08:52:20 

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

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!