JetProducers/test/analyzeJetsChain.C

0001 //File: analyzeJets.C
0002 //Author: R. Harris
0003 //Date: June 22, 2006.
0004 //Description: Framework lite analysis of jets in root using TChain
0005 //Pre-requisites: To use with CINT:
0006 //                   .x analyzeJets_headChain.C
0007 //                   .x analyzeJetsChain.C
0008 //                To use with ACliC in compiled mode
0009 //                   .x prepareCompile.C
0010 //                   .x analyzeJets_headChain.C
0011 //                   .x analyzeJetsChain.C++
0012 //     You will need to replace evtgen_jets*.root with
0013 //     your chain of root files in analyzeJets*.C
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   // Create histo file and book histograms
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   // Declare CaloJetCollection.
0035   std::vector<reco::CaloJet> CaloJetCollection;
0036
0037   #ifndef __CINT__
0038     // For the compiled version we need to define the chain here
0039     TChain chain("Events");
0040     chain.Add("evtgen_jets.root");
0041     chain.Add("evtgen_jets2.root");
0042    #endif
0043
0044    // Number of entries in chain
0045   Int_t   nevent = chain.GetEntries();
0046
0047   // Tell root we only want the CaloJets branches.
0048   chain.SetBranchStatus("*",0);
0049   chain.SetBranchStatus("recoCaloJets*",1);
0050
0051   // Loop over events
0052   for ( int index = 0; index < nevent; ++index ) {
0053
0054     // Load the TTree corresponding to event index and return the entry with respect to that tree.
0055     int current = chain.LoadTree(index);
0056
0057     // Check if we are on the first entry in the tree
0058     if (current==0) {
0059
0060        //Read the first entry in this new tree, needed to set branch address.
0061        chain.GetEvent(index);
0062
0063        //Set the branch address for this new tree
0064        chain.SetBranchAddress(chain.GetAlias("MC5CaloJet"),&CaloJetCollection);
0065     }
0066
0067     // Read the event.
0068     chain.GetEvent(index);
0069
0070     double px[2], py[2], pz[2], E[2];
0071     std::cout << "Entry index: " << index << std::endl;
0072     //chain.SetBranchAddress("CaloJets_midPointCone5CaloJets.obj",&CaloJetCollection);
0073     int numJets = CaloJetCollection.size();
0074     std::cout << "Num Jets: " << numJets << std::endl;
0075
0076     //Loop over jets
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       //Get and printout jet pt, eta, phi for all jets
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         //Fill Histograms for two highest pt jets
0094         h_pt->Fill(pt);
0095     h_eta->Fill(eta);
0096     h_phi->Fill(phi);
0097
0098        //Get Lorentz Vector components of two highest pt jets
0099        px[jetIndex] = Jet->px();
0100        py[jetIndex] = Jet->py();
0101        pz[jetIndex] = Jet->pz();
0102        E[jetIndex]  = Jet->energy();
0103       }
0104     }
0105     //Printout Dijet Mass and Fill Dijet Mass histogram
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   // save histograms
0116   histofile.Write();
0117   histofile.Close();
0118 }