JetProducers/test/analyzeJets.C

0001 { //File: analyzeJets.C
0002   //Author: R. Harris from example of Oliver Gutsche
0003   //Date: June 16, 2006.
0004   //Description: Framework lite analysis of jets in root
0005   //Pre-requisites: Loading and enabling framework lite libraries and event file with analyzeJets_head.C
0006
0007   // Create histo file and book histograms
0008   TFile histofile("jet_hists.root","RECREATE");
0009   TH1F* h_pt = new TH1F("pt","Leading Jets pT (GeV)",100,0.0,1000.0);
0010   TH1F* h_eta = new TH1F("eta","Leading Jets eta",100,-5.0,5.0);
0011   TH1F* h_phi = new TH1F("phi","Leading Jets phi",72,-3.141527,3.141527);
0012   TH1F* h_m2j = new TH1F("m2j","Dijet Mass",100,0.0,1000.0);
0013
0014   //Get event tree and jet collection branch
0015   TTree *tree = (TTree*)file.Get("Events");
0016   std::vector<reco::CaloJet> CaloJetCollection;
0017   TBranch *branch = tree->GetBranch(tree->GetAlias("MC5CaloJet"));
0018
0019   branch->SetAddress(&CaloJetCollection);
0020
0021   // Loop over events
0022   for ( unsigned int index = 0; index < tree->GetEntries(); ++index ) {
0023     double px[2], py[2], pz[2], E[2];
0024     std::cout << "Entry index: " << index << std::endl;
0025     branch->GetEntry(index);
0026     int numJets = CaloJetCollection.size();
0027     std::cout << "Num Jets: " << numJets << std::endl;
0028
0029     //Loop over jets
0030     for ( unsigned int jetIndex = 0; jetIndex < CaloJetCollection.size(); ++jetIndex ) {
0031       std::cout << "jet" << jetIndex  ;
0032       reco::CaloJet* Jet = (reco::CaloJet*)CaloJetCollection[jetIndex];
0033
0034       //Get and printout jet pt, eta, phi for all jets
0035       double pt = Jet->pt();    std::cout << ": pt=" << pt;
0036       double eta = Jet->eta();  std::cout << ", eta=" << eta;
0037       double phi = Jet->phi();  std::cout << ", phi=" << phi << std::endl;
0038
0039       if(jetIndex<2)
0040       {
0041
0042         //Fill Histograms for two highest pt jets
0043         h_pt->Fill(pt);
0044     h_eta->Fill(eta);
0045     h_phi->Fill(phi);
0046
0047        //Get Lorentz Vector components of two highest pt jets
0048        px[jetIndex] = Jet->px();
0049        py[jetIndex] = Jet->py();
0050        pz[jetIndex] = Jet->pz();
0051        E[jetIndex]  = Jet->energy();
0052       }
0053     }
0054     //Printout Dijet Mass and Fill Dijet Mass histogram
0055     if( numJets >= 2 ){
0056       double DijetMass = sqrt( pow(E[0]+E[1],2) - pow(px[0]+px[1],2)
0057                                                 - pow(py[0]+py[1],2)
0058                                                 - pow(pz[0]+pz[1],2) );
0059       std::cout << "Dijet Mass = " << DijetMass  << std::endl;
0060       h_m2j->Fill(DijetMass);
0061
0062     }
0063   }
0064   // save histograms
0065   histofile.Write();
0066   histofile.Close();
0067 }