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="")
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> recoClustIndex = {tree_reader, "_ReconstructedChargedJets_clusters.index"};
0031   TTreeReaderArray<int> recoTrackIndex = {tree_reader, "_ReconstructedChargedJets_tracks.index"};
0032   TTreeReaderArray<int> recoPartIndex = {tree_reader, "_ReconstructedChargedJets_particles.index"};
0033
0034   // Reconstructed Particles
0035   TTreeReaderArray<float> recoPartMomX = {tree_reader, "ReconstructedChargedParticles.momentum.x"};
0036   TTreeReaderArray<float> recoPartMomY = {tree_reader, "ReconstructedChargedParticles.momentum.y"};
0037   TTreeReaderArray<float> recoPartMomZ = {tree_reader, "ReconstructedChargedParticles.momentum.z"};
0038   TTreeReaderArray<float> recoPartM = {tree_reader, "ReconstructedChargedParticles.mass"};
0039
0040   // Generated Jets
0041   TTreeReaderArray<int> genType = {tree_reader, "GeneratedChargedJets.type"};
0042   TTreeReaderArray<float> genNRG = {tree_reader, "GeneratedChargedJets.energy"};
0043   TTreeReaderArray<int> genPDG = {tree_reader, "GeneratedChargedJets.PDG"};
0044   TTreeReaderArray<float> genMomX = {tree_reader, "GeneratedChargedJets.momentum.x"};
0045   TTreeReaderArray<float> genMomY = {tree_reader, "GeneratedChargedJets.momentum.y"};
0046   TTreeReaderArray<float> genMomZ = {tree_reader, "GeneratedChargedJets.momentum.z"};
0047   TTreeReaderArray<float> genM = {tree_reader, "GeneratedChargedJets.mass"};
0048   TTreeReaderArray<unsigned int> genPartsBegin = {tree_reader, "GeneratedChargedJets.particles_begin"};
0049   TTreeReaderArray<unsigned int> genPartsEnd = {tree_reader, "GeneratedChargedJets.particles_end"};
0050
0051   TTreeReaderArray<int> genPartIndex = {tree_reader, "_GeneratedChargedJets_particles.index"};
0052
0053   // MC
0054   TTreeReaderArray<float> mcMomX = {tree_reader, "GeneratedParticles.momentum.x"};
0055   TTreeReaderArray<float> mcMomY = {tree_reader, "GeneratedParticles.momentum.y"};
0056   TTreeReaderArray<float> mcMomZ = {tree_reader, "GeneratedParticles.momentum.z"};
0057   TTreeReaderArray<float> mcM = {tree_reader, "GeneratedParticles.mass"};
0058   TTreeReaderArray<int> pdg = {tree_reader, "GeneratedParticles.PDG"};
0059
0060   // Define Histograms
0061   // Reco
0062   TH1D *numRecoJetsEventHist = new TH1D("numRecoJetsEvent","",20,0.,20.);
0063   TH2D *recoJetEvsEtaHist = new TH2D("recoJetEvsEta","",100,-5.,5.,300,0.,300.);
0064   TH2D *recoJetPhiVsEtaHist = new TH2D("recoJetPhiVsEta","",100,-5.,5.,100,-TMath::Pi(),TMath::Pi());
0065
0066   TH1D *numRecoJetPartsHist = new TH1D("numRecoJetParts","",20,0.,20.);
0067   TH2D *recoJetPartEvsEtaHist = new TH2D("recoJetPartEvsEta","",100,-5.,5.,300,0.,300.);
0068   TH2D *recoJetPartPhiVsEtaHist = new TH2D("recoJetPartPhiVsEta","",100,-5.,5.,100,-TMath::Pi(),TMath::Pi());
0069
0070   TH1D *recoJetPartDeltaRHist = new TH1D("recoJetPartDeltaR","",5000,0.,5.);
0071
0072   TH2D *recoJetEvsPartESumHist = new TH2D("recoJetEvsPartESum","",3000,0.,300.,3000,0.,300.);
0073   TH1D *recoJetEDiffHist = new TH1D("recoJetEDiff","",500,-10.,10.);
0074
0075   TH2D *recoJetEvsEtaBadHist = new TH2D("recoJetEvsEtaBad","",100,-5.,5.,300,0.,300.);
0076   TH2D *recoJetPhiVsEtaBadHist = new TH2D("recoJetPhiVsEtaBad","",100,-5.,5.,100,-TMath::Pi(),TMath::Pi());
0077
0078   // Gen
0079   TH1D *numGenJetsEventHist = new TH1D("numGenJetsEvent","",20,0.,20.);
0080   TH2D *genJetEvsEtaHist = new TH2D("genJetEvsEta","",100,-5.,5.,300,0.,300.);
0081   TH2D *genJetPhiVsEtaHist = new TH2D("genJetPhiVsEta","",100,-5.,5.,100,-TMath::Pi(),TMath::Pi());
0082
0083   TH1D *numGenJetPartsHist = new TH1D("numGenJetParts","",20,0.,20.);
0084   TH2D *genJetPartEvsEtaHist = new TH2D("genJetPartEvsEta","",100,-5.,5.,300,0.,300.);
0085   TH2D *genJetPartPhiVsEtaHist = new TH2D("genJetPartPhiVsEta","",100,-5.,5.,100,-TMath::Pi(),TMath::Pi());
0086
0087   TH1D *genJetPartDeltaRHist = new TH1D("genJetPartDeltaR","",5000,0.,5.);
0088
0089   TH2D *genJetEvsPartESumHist = new TH2D("genJetEvsPartESum","",3000,0.,300.,3000,0.,300.);
0090   TH1D *genJetEDiffHist = new TH1D("genJetEDiff","",500,-10.,10.);
0091
0092   TH2D *genJetEvsEtaBadHist = new TH2D("genJetEvsEtaBad","",100,-5.,5.,300,0.,300.);
0093   TH2D *genJetPhiVsEtaBadHist = new TH2D("genJetPhiVsEtaBad","",100,-5.,5.,100,-TMath::Pi(),TMath::Pi());
0094
0095
0096   // Loop Through Events
0097   int NEVENTS = 0;
0098   while(tree_reader.Next()) {
0099
0100     if(NEVENTS%10000 == 0) cout << "Events Processed: " << NEVENTS << endl;
0101
0102     // Analyze Reonstructed Jets
0103     numRecoJetsEventHist->Fill(recoType.GetSize());
0104     for(unsigned int i=0; i<recoType.GetSize(); i++)
0105       {
0106     TVector3 jetMom(recoMomX[i],recoMomY[i],recoMomZ[i]);
0107
0108     recoJetEvsEtaHist->Fill(jetMom.PseudoRapidity(),recoNRG[i]);
0109     recoJetPhiVsEtaHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
0110
0111     double esum = 0.0;
0112     for(unsigned int j=partsBegin[i]; j<partsEnd[i]; j++)
0113       {
0114         // partsbegin and partsEnd specify the entries from _ReconstructedChargedJets_particles.index that make up the jet
0115         // _ReconstructedChargedJets_particles.index stores the ReconstructedChargedParticles index of the jet constituent
0116         double mX = recoPartMomX[recoPartIndex[j]];
0117         double mY = recoPartMomY[recoPartIndex[j]];
0118         double mZ = recoPartMomZ[recoPartIndex[j]];
0119         double mM = recoPartM[recoPartIndex[j]];
0120
0121         double tmpE = TMath::Sqrt(mX*mX + mY*mY + mZ*mZ + mM*mM);
0122
0123         esum += tmpE;
0124
0125         TVector3 partMom(mX,mY,mZ);
0126
0127         recoJetPartEvsEtaHist->Fill(partMom.PseudoRapidity(),tmpE);
0128         recoJetPartPhiVsEtaHist->Fill(partMom.PseudoRapidity(),partMom.Phi());
0129
0130         // Distance between jet axis and particle
0131         double dEta = jetMom.PseudoRapidity() - partMom.PseudoRapidity();
0132         double dPhi = TVector2::Phi_mpi_pi(jetMom.Phi() - partMom.Phi());
0133         double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
0134
0135         recoJetPartDeltaRHist->Fill(dR);
0136       }
0137
0138     numRecoJetPartsHist->Fill(partsEnd[i] - partsBegin[i]);
0139     recoJetEvsPartESumHist->Fill(recoNRG[i],esum);
0140     recoJetEDiffHist->Fill(recoNRG[i]-esum);
0141
0142     if(TMath::Abs(esum - recoNRG[i]) > 0.00001)
0143       {
0144         recoJetEvsEtaBadHist->Fill(jetMom.PseudoRapidity(),recoNRG[i]);
0145         recoJetPhiVsEtaBadHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
0146       }
0147       }
0148
0149
0150     // Analyze Generated Jets
0151     numRecoJetsEventHist->Fill(genType.GetSize());
0152     for(unsigned int i=0; i<genType.GetSize(); i++)
0153       {
0154     TVector3 jetMom(genMomX[i],genMomY[i],genMomZ[i]);
0155
0156     genJetEvsEtaHist->Fill(jetMom.PseudoRapidity(),genNRG[i]);
0157     genJetPhiVsEtaHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
0158
0159     double esumG = 0.0;
0160     for(unsigned int j=genPartsBegin[i]; j<genPartsEnd[i]; j++)
0161       {
0162         double mX = mcMomX[genPartIndex[j]];
0163         double mY = mcMomY[genPartIndex[j]];
0164         double mZ = mcMomZ[genPartIndex[j]];
0165         double mM = mcM[genPartIndex[j]];
0166
0167         double tmpE = TMath::Sqrt(mX*mX + mY*mY + mZ*mZ + mM*mM);
0168
0169         esumG += tmpE;
0170
0171         TVector3 partMom(mX,mY,mZ);
0172
0173         genJetPartEvsEtaHist->Fill(partMom.PseudoRapidity(),tmpE);
0174         genJetPartPhiVsEtaHist->Fill(partMom.PseudoRapidity(),partMom.Phi());
0175
0176         // Distance between jet axis and particle
0177         double dEta = jetMom.PseudoRapidity() - partMom.PseudoRapidity();
0178         double dPhi = TVector2::Phi_mpi_pi(jetMom.Phi() - partMom.Phi());
0179         double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
0180
0181         genJetPartDeltaRHist->Fill(dR);
0182       }
0183
0184     numGenJetPartsHist->Fill(genPartsEnd[i] - genPartsBegin[i]);
0185     genJetEvsPartESumHist->Fill(genNRG[i],esumG);
0186     genJetEDiffHist->Fill(genNRG[i]-esumG);
0187
0188     if(TMath::Abs(esumG - genNRG[i]) > 0.00001)
0189       {
0190         genJetEvsEtaBadHist->Fill(jetMom.PseudoRapidity(),genNRG[i]);
0191         genJetPhiVsEtaBadHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
0192       }
0193       }
0194
0195     NEVENTS++;
0196   }
0197
0198   ofile->Write();
0199   ofile->Close();
0200
0201 }