doc/WorkBook/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 using namespace reco;
0008
0009   // Create histo file and book histograms
0010   TFile histofile("jet_hists.root","RECREATE");
0011   TH1F* h_pt = new TH1F("pt","Leading Jets pT (GeV)",100,0.0,500.0);
0012   TH1F* h_eta = new TH1F("eta","Leading Jets eta",100,-5.0,5.0);
0013   TH1F* h_phi = new TH1F("phi","Leading Jets phi",72,-3.141527,3.141527);
0014   TH1F* h_m2j = new TH1F("m2j","Dijet Mass",100,0.0,1000.0);
0015
0016   //Get event tree and jet collection branch
0017   TTree *tree = (TTree*)file->Get("Events");
0018   std::vector<CaloJet> CaloJetCollection;
0019   TBranch *branch = tree->GetBranch(tree->GetAlias("MC5CaloJet"));
0020   branch->SetAddress(&CaloJetCollection);
0021
0022   // Loop over events
0023   for ( unsigned int index = 0; index < tree->GetEntries(); ++index ) {
0024     double px[2], py[2], pz[2], E[2];
0025     std::cout << "Entry index: " << index << std::endl;
0026     branch->GetEntry(index);
0027     int numJets = CaloJetCollection.size();
0028     std::cout << "Num Jets: " << numJets << std::endl;
0029
0030     //Loop over jets
0031     for ( unsigned int jetIndex = 0; jetIndex < CaloJetCollection.size(); ++jetIndex ) {
0032       std::cout << "jet" << jetIndex  ;
0033       CaloJet* Jet = (CaloJet*)CaloJetCollection[jetIndex];
0034
0035       //Get and printout jet pt, eta, phi for all jets
0036       double pt = Jet->ptUncached();    std::cout << ": pt=" << pt;
0037       double eta = Jet->etaUncached();  std::cout << ", eta=" << eta;
0038       double phi = Jet->phiUncached();  std::cout << ", phi=" << phi << std::endl;
0039
0040       if(jetIndex<2)
0041       {
0042
0043         //Fill Histograms for two highest pt jets
0044         h_pt->Fill(pt);
0045     h_eta->Fill(eta);
0046     h_phi->Fill(phi);
0047
0048        //Get Lorentz Vector components of two highest pt jets
0049        px[jetIndex] = Jet->px();
0050        py[jetIndex] = Jet->py();
0051        pz[jetIndex] = Jet->pz();
0052        E[jetIndex]  = Jet->energy();
0053       }
0054     }
0055     //Printout Dijet Mass and Fill Dijet Mass histogram
0056     if( numJets >= 2 ){
0057       double DijetMass = sqrt( pow(E[0]+E[1],2) - pow(px[0]+px[1],2)
0058                                                 - pow(py[0]+py[1],2)
0059                                                 - pow(pz[0]+pz[1],2) );
0060       std::cout << "Dijet Mass = " << DijetMass  << std::endl;
0061       h_m2j->Fill(DijetMass);
0062
0063     }
0064   }
0065   // save histograms
0066   histofile.Write();
0067   histofile.Close();
0068 }