File indexing completed on 2024-04-06 12:25:36
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015 #ifndef __CINT__
0016 #include <vector>
0017 #include "TChain.h"
0018 #include "TH1.h"
0019 #include "TFile.h"
0020 #include <iostream>
0021 #include "DataFormats/JetReco/interface/CaloJet.h"
0022 #endif
0023
0024 void analyzeJetsChain()
0025 {
0026
0027
0028 TFile histofile("jet_hists.root","RECREATE");
0029 TH1F* h_pt = new TH1F("pt","Leading Jets pT (GeV)",100,0.0,1000.0);
0030 TH1F* h_eta = new TH1F("eta","Leading Jets eta",100,-5.0,5.0);
0031 TH1F* h_phi = new TH1F("phi","Leading Jets phi",72,-3.141527,3.141527);
0032 TH1F* h_m2j = new TH1F("m2j","Dijet Mass",100,0.0,1000.0);
0033
0034
0035 std::vector<reco::CaloJet> CaloJetCollection;
0036
0037 #ifndef __CINT__
0038
0039 TChain chain("Events");
0040 chain.Add("evtgen_jets.root");
0041 chain.Add("evtgen_jets2.root");
0042 #endif
0043
0044
0045 Int_t nevent = chain.GetEntries();
0046
0047
0048 chain.SetBranchStatus("*",0);
0049 chain.SetBranchStatus("recoCaloJets*",1);
0050
0051
0052 for ( int index = 0; index < nevent; ++index ) {
0053
0054
0055 int current = chain.LoadTree(index);
0056
0057
0058 if (current==0) {
0059
0060
0061 chain.GetEvent(index);
0062
0063
0064 chain.SetBranchAddress(chain.GetAlias("MC5CaloJet"),&CaloJetCollection);
0065 }
0066
0067
0068 chain.GetEvent(index);
0069
0070 double px[2], py[2], pz[2], E[2];
0071 std::cout << "Entry index: " << index << std::endl;
0072
0073 int numJets = CaloJetCollection.size();
0074 std::cout << "Num Jets: " << numJets << std::endl;
0075
0076
0077 for ( unsigned int jetIndex = 0; jetIndex < CaloJetCollection.size(); ++jetIndex ) {
0078 std::cout << "jet" << jetIndex ;
0079 #ifndef __CINT__
0080 reco::CaloJet* Jet= &(CaloJetCollection[jetIndex]);
0081 #else
0082 reco::CaloJet* Jet = (reco::CaloJet*)CaloJetCollection[jetIndex];
0083 #endif
0084
0085
0086 double pt = Jet->pt(); std::cout << ": pt=" << pt;
0087 double eta = Jet->eta(); std::cout << ", eta=" << eta;
0088 double phi = Jet->phi(); std::cout << ", phi=" << phi << std::endl;
0089
0090 if(jetIndex<2)
0091 {
0092
0093
0094 h_pt->Fill(pt);
0095 h_eta->Fill(eta);
0096 h_phi->Fill(phi);
0097
0098
0099 px[jetIndex] = Jet->px();
0100 py[jetIndex] = Jet->py();
0101 pz[jetIndex] = Jet->pz();
0102 E[jetIndex] = Jet->energy();
0103 }
0104 }
0105
0106 if( numJets >= 2 ){
0107 double DijetMass = sqrt( pow(E[0]+E[1],2) - pow(px[0]+px[1],2)
0108 - pow(py[0]+py[1],2)
0109 - pow(pz[0]+pz[1],2) );
0110 std::cout << "Dijet Mass = " << DijetMass << std::endl;
0111 h_m2j->Fill(DijetMass);
0112
0113 }
0114 }
0115
0116 histofile.Write();
0117 histofile.Close();
0118 }