CMSSW Reference Manual

00001 // DijetMass.cc
00002 // Description:  Some Basic validation plots for jets.
00003 // Author: Robert M. Harris
00004 // Date:  30 - August - 2006
00005 // 
00006 #include "RecoJets/JetAnalyzers/interface/DijetMass.h"
00007 #include "RecoJets/JetAlgorithms/interface/JetAlgoHelper.h"
00008 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
00009 #include "DataFormats/JetReco/interface/CaloJet.h"
00010 #include "DataFormats/JetReco/interface/GenJet.h"
00011 #include "PhysicsTools/Utilities/interface/deltaR.h"
00012 #include "FWCore/Framework/interface/Event.h"
00013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00014 #include <TROOT.h>
00015 #include <TSystem.h>
00016 #include <TFile.h>
00017 #include <TCanvas.h>
00018 #include <cmath>
00019 using namespace edm;
00020 using namespace reco;
00021 using namespace std;
00022 
00023 // Get the algorithm of the jet collections we will read from the .cfg file 
00024 // which defines the value of the strings CaloJetAlgorithm and GenJetAlgorithm.
00025 DijetMass::DijetMass( const ParameterSet & cfg ) :
00026   PtHistMax( cfg.getParameter<double>( "PtHistMax" ) ),
00027   GenType( cfg.getParameter<string>( "GenType" ) )
00028   {
00029 }
00030 
00031 void DijetMass::beginJob( const EventSetup & ) {
00032   cout << "DijetMass: Maximum bin edge for Pt Hists = " << PtHistMax << endl;
00033   numJets=2;
00034 
00035   //Initialize some stuff
00036   evtCount = 0;
00037 
00038   // Open the histogram file and book some associated histograms
00039   m_file=new TFile("DijetMassHistos.root","RECREATE"); 
00040   
00041   //Simple histos
00042   
00043   //MC5 cal
00044   ptMC5cal  = TH1F("ptMC5cal","p_{T} of leading CaloJets (MC5)",50,0.0,PtHistMax);
00045   etaMC5cal = TH1F("etaMC5cal","#eta of leading CaloJets (MC5)",23,-1.0,1.0);
00046   phiMC5cal = TH1F("phiMC5cal","#phi of leading CaloJets (MC5)",72,-M_PI, M_PI);
00047   m2jMC5cal = TH1F("m2jMC5cal","Dijet Mass of leading CaloJets (MC5)",100,0.0,2*PtHistMax);
00048 
00049   //MC5 gen
00050   ptMC5gen = TH1F("ptMC5gen","p_{T} of leading genJets (MC5)",50,0.0,PtHistMax);
00051   etaMC5gen = TH1F("etaMC5gen","#eta of leading genJets (MC5)",23,-1.0,1.0);
00052   phiMC5gen = TH1F("phiMC5gen","#phi of leading genJets (MC5)",72,-M_PI, M_PI);
00053   m2jMC5gen = TH1F("m2jMC5gen","Dijet Mass of leading genJets (MC5)",100,0.0,2*PtHistMax);
00054 
00055   //MC5 cor
00056   ptMC5cor  = TH1F("ptMC5cor","p_{T} of leading Corrected CaloJets (MC5)",50,0.0,PtHistMax);
00057   etaMC5cor = TH1F("etaMC5cor","#eta of leading Corrected CaloJets (MC5)",23,-1.0,1.0);
00058   phiMC5cor = TH1F("phiMC5cor","#phi of leading Corrected CaloJets (MC5)",72,-M_PI, M_PI);
00059   m2jMC5cor = TH1F("m2jMC5cor","Dijet Mass of leading Corrected CaloJets (MC5)",100,0.0,2*PtHistMax);
00060 
00061   //IC5 cal
00062   ptIC5cal  = TH1F("ptIC5cal","p_{T} of leading CaloJets (IC5)",50,0.0,PtHistMax);
00063   etaIC5cal = TH1F("etaIC5cal","#eta of leading CaloJets (IC5)",23,-1.0,1.0);
00064   phiIC5cal = TH1F("phiIC5cal","#phi of leading CaloJets (IC5)",72,-M_PI, M_PI);
00065   m2jIC5cal = TH1F("m2jIC5cal","Dijet Mass of leading CaloJets (IC5)",100,0.0,2*PtHistMax);
00066 
00067   //IC5 gen
00068   ptIC5gen = TH1F("ptIC5gen","p_{T} of leading genJets (IC5)",50,0.0,PtHistMax);
00069   etaIC5gen = TH1F("etaIC5gen","#eta of leading genJets (IC5)",23,-1.0,1.0);
00070   phiIC5gen = TH1F("phiIC5gen","#phi of leading genJets (IC5)",72,-M_PI, M_PI);
00071   m2jIC5gen = TH1F("m2jIC5gen","Dijet Mass of leading genJets (IC5)",100,0.0,2*PtHistMax);
00072 
00073   //IC5 cor
00074   ptIC5cor  = TH1F("ptIC5cor","p_{T} of leading Corrected CaloJets (IC5)",50,0.0,PtHistMax);
00075   etaIC5cor = TH1F("etaIC5cor","#eta of leading Corrected CaloJets (IC5)",23,-1.0,1.0);
00076   phiIC5cor = TH1F("phiIC5cor","#phi of leading Corrected CaloJets (IC5)",72,-M_PI, M_PI);
00077   m2jIC5cor = TH1F("m2jIC5cor","Dijet Mass of leading Corrected CaloJets (IC5)",100,0.0,2*PtHistMax);
00078 
00079   //KT10 cal
00080   ptKT10cal  = TH1F("ptKT10cal","p_{T} of leading CaloJets (KT10)",50,0.0,PtHistMax);
00081   etaKT10cal = TH1F("etaKT10cal","#eta of leading CaloJets (KT10)",23,-1.0,1.0);
00082   phiKT10cal = TH1F("phiKT10cal","#phi of leading CaloJets (KT10)",72,-M_PI, M_PI);
00083   m2jKT10cal = TH1F("m2jKT10cal","Dijet Mass of leading CaloJets (KT10)",100,0.0,2*PtHistMax);
00084 
00085   //KT10 gen
00086   ptKT10gen = TH1F("ptKT10gen","p_{T} of leading genJets (KT10)",50,0.0,PtHistMax);
00087   etaKT10gen = TH1F("etaKT10gen","#eta of leading genJets (KT10)",23,-1.0,1.0);
00088   phiKT10gen = TH1F("phiKT10gen","#phi of leading genJets (KT10)",72,-M_PI, M_PI);
00089   m2jKT10gen = TH1F("m2jKT10gen","Dijet Mass of leading genJets (KT10)",100,0.0,2*PtHistMax);
00090 
00091 
00092   //Matched jets Analysis Histograms for MC5 CaloJets only
00093   dR = TH1F("dR","Leading genJets dR with matched CaloJet",100,0,0.5);
00094   respVsPt = TProfile("respVsPt","CaloJet Response of Leading genJets in Barrel",100,0.0,PtHistMax/2); 
00095   dRcor = TH1F("dRcor","CorJets dR with matched CaloJet",100,0.0,0.01);
00096   corRespVsPt = TProfile("corRespVsPt","Corrected CaloJet Response of Leading genJets in Barrel",100,0.0,PtHistMax/2); 
00097 }
00098 
00099 void DijetMass::analyze( const Event& evt, const EventSetup& es ) {
00100 
00101   evtCount++;
00102   math::XYZTLorentzVector p4jet[2], p4gen[2], p4cal[2], p4cor[2];
00103   int jetInd;
00104   Handle<CaloJetCollection> caloJets;
00105   Handle<GenJetCollection> genJets;
00106 
00107   //Fill Simple Histos
00108   
00109   //MC5 cal
00110   evt.getByLabel( "midPointCone5CaloJets", caloJets );
00111   jetInd = 0;
00112   for( CaloJetCollection::const_iterator cal = caloJets->begin(); cal != caloJets->end() && jetInd<2; ++ cal ) {
00113     p4jet[jetInd] = cal->p4();
00114     jetInd++;
00115   }
00116   if(jetInd==2&&abs(p4jet[0].eta())<1.0&&abs(p4jet[1].eta())<1.0){
00117      m2jMC5cal.Fill( (p4jet[0]+p4jet[1]).mass() ); 
00118      ptMC5cal.Fill( p4jet[0].Pt() ); ptMC5cal.Fill( p4jet[1].Pt() );  
00119      etaMC5cal.Fill( p4jet[0].eta() ); etaMC5cal.Fill( p4jet[1].eta() );  
00120      phiMC5cal.Fill( p4jet[0].phi() ); phiMC5cal.Fill( p4jet[1].phi() );  
00121    }
00122 
00123   //MC5 gen
00124   evt.getByLabel( "midPointCone5GenJets", genJets );
00125   jetInd = 0;
00126   for( GenJetCollection::const_iterator gen = genJets->begin(); gen != genJets->end() && jetInd<2; ++ gen ) {
00127     p4jet[jetInd] = gen->p4();
00128     jetInd++;
00129   }
00130   if(jetInd==2&&abs(p4jet[0].eta())<1.0&&abs(p4jet[1].eta())<1.0){
00131      m2jMC5gen.Fill( (p4jet[0]+p4jet[1]).mass() ); 
00132      ptMC5gen.Fill( p4jet[0].Pt() ); ptMC5gen.Fill( p4jet[1].Pt() );  
00133      etaMC5gen.Fill( p4jet[0].eta() ); etaMC5gen.Fill( p4jet[1].eta() );  
00134      phiMC5gen.Fill( p4jet[0].phi() ); phiMC5gen.Fill( p4jet[1].phi() );  
00135    }
00136 
00137   //MC5 cal
00138   evt.getByLabel( "corJetMcone5", caloJets );
00139   jetInd = 0;
00140   for( CaloJetCollection::const_iterator cal = caloJets->begin(); cal != caloJets->end() && jetInd<2; ++ cal ) {
00141     p4jet[jetInd] = cal->p4();
00142     jetInd++;
00143   }
00144   if(jetInd==2&&abs(p4jet[0].eta())<1.0&&abs(p4jet[1].eta())<1.0){
00145      m2jMC5cor.Fill( (p4jet[0]+p4jet[1]).mass() ); 
00146      ptMC5cor.Fill( p4jet[0].Pt() ); ptMC5cor.Fill( p4jet[1].Pt() );  
00147      etaMC5cor.Fill( p4jet[0].eta() ); etaMC5cor.Fill( p4jet[1].eta() );  
00148      phiMC5cor.Fill( p4jet[0].phi() ); phiMC5cor.Fill( p4jet[1].phi() );  
00149    }
00150 
00151 
00152   //IC5 cal
00153   evt.getByLabel( "iterativeCone5CaloJets", caloJets );
00154   jetInd = 0;
00155   for( CaloJetCollection::const_iterator cal = caloJets->begin(); cal != caloJets->end() && jetInd<2; ++ cal ) {
00156     p4jet[jetInd] = cal->p4();
00157     jetInd++;
00158   }
00159   if(jetInd==2&&abs(p4jet[0].eta())<1.0&&abs(p4jet[1].eta())<1.0){
00160      m2jIC5cal.Fill( (p4jet[0]+p4jet[1]).mass() ); 
00161      ptIC5cal.Fill( p4jet[0].Pt() ); ptIC5cal.Fill( p4jet[1].Pt() );  
00162      etaIC5cal.Fill( p4jet[0].eta() ); etaIC5cal.Fill( p4jet[1].eta() );  
00163      phiIC5cal.Fill( p4jet[0].phi() ); phiIC5cal.Fill( p4jet[1].phi() );  
00164    }
00165 
00166 
00167   //IC5 gen
00168   evt.getByLabel( "iterativeCone5GenJets", genJets );
00169   jetInd = 0;
00170   for( GenJetCollection::const_iterator gen = genJets->begin(); gen != genJets->end() && jetInd<2; ++ gen ) {
00171     p4jet[jetInd] = gen->p4();
00172     jetInd++;
00173   }
00174   if(jetInd==2&&abs(p4jet[0].eta())<1.0&&abs(p4jet[1].eta())<1.0){
00175      m2jIC5gen.Fill( (p4jet[0]+p4jet[1]).mass() ); 
00176      ptIC5gen.Fill( p4jet[0].Pt() ); ptIC5gen.Fill( p4jet[1].Pt() );  
00177      etaIC5gen.Fill( p4jet[0].eta() ); etaIC5gen.Fill( p4jet[1].eta() );  
00178      phiIC5gen.Fill( p4jet[0].phi() ); phiIC5gen.Fill( p4jet[1].phi() );  
00179    }
00180 
00181   //IC5 cor
00182   evt.getByLabel( "corJetIcone5", caloJets );
00183   jetInd = 0;
00184   for( CaloJetCollection::const_iterator cal = caloJets->begin(); cal != caloJets->end() && jetInd<2; ++ cal ) {
00185     p4jet[jetInd] = cal->p4();
00186     jetInd++;
00187   }
00188   if(jetInd==2&&abs(p4jet[0].eta())<1.0&&abs(p4jet[1].eta())<1.0){
00189      m2jIC5cor.Fill( (p4jet[0]+p4jet[1]).mass() ); 
00190      ptIC5cor.Fill( p4jet[0].Pt() ); ptIC5cor.Fill( p4jet[1].Pt() );  
00191      etaIC5cor.Fill( p4jet[0].eta() ); etaIC5cor.Fill( p4jet[1].eta() );  
00192      phiIC5cor.Fill( p4jet[0].phi() ); phiIC5cor.Fill( p4jet[1].phi() );  
00193    }
00194 
00195   //KT10 cal
00196   evt.getByLabel( "ktCaloJets", caloJets );
00197   jetInd = 0;
00198   for( CaloJetCollection::const_iterator cal = caloJets->begin(); cal != caloJets->end() && jetInd<2; ++ cal ) {
00199     p4jet[jetInd] = cal->p4();
00200     jetInd++;
00201   }
00202   if(jetInd==2&&abs(p4jet[0].eta())<1.0&&abs(p4jet[1].eta())<1.0){
00203      m2jKT10cal.Fill( (p4jet[0]+p4jet[1]).mass() ); 
00204      ptKT10cal.Fill( p4jet[0].Pt() ); ptKT10cal.Fill( p4jet[1].Pt() );  
00205      etaKT10cal.Fill( p4jet[0].eta() ); etaKT10cal.Fill( p4jet[1].eta() );  
00206      phiKT10cal.Fill( p4jet[0].phi() ); phiKT10cal.Fill( p4jet[1].phi() );  
00207    }
00208 
00209   //KT10 gen
00210   evt.getByLabel( "ktGenJets", genJets );
00211   jetInd = 0;
00212   for( GenJetCollection::const_iterator gen = genJets->begin(); gen != genJets->end() && jetInd<2; ++ gen ) {
00213     p4jet[jetInd] = gen->p4();
00214     jetInd++;
00215   }
00216   if(jetInd==2&&abs(p4jet[0].eta())<1.0&&abs(p4jet[1].eta())<1.0){
00217      m2jKT10gen.Fill( (p4jet[0]+p4jet[1]).mass() ); 
00218      ptKT10gen.Fill( p4jet[0].Pt() ); ptKT10gen.Fill( p4jet[1].Pt() );  
00219      etaKT10gen.Fill( p4jet[0].eta() ); etaKT10gen.Fill( p4jet[1].eta() );  
00220      phiKT10gen.Fill( p4jet[0].phi() ); phiKT10gen.Fill( p4jet[1].phi() );  
00221    }
00222 
00223   //Matching for MC5 Jets: leading genJets matched to any CaloJet
00224   evt.getByLabel( "midPointCone5GenJets", genJets );
00225   evt.getByLabel( "midPointCone5CaloJets", caloJets );
00226   jetInd = 0;
00227   double dRmin[2];
00228   for( GenJetCollection::const_iterator gen = genJets->begin(); gen != genJets->end() && jetInd<numJets; ++ gen ) { //leading genJets
00229     p4gen[jetInd] = gen->p4();    //Gen 4-vector
00230     dRmin[jetInd]=1000.0;
00231     for( CaloJetCollection::const_iterator cal = caloJets->begin(); cal != caloJets->end(); ++ cal ) { //all CaloJets
00232        double delR = deltaR( cal->eta(), cal->phi(), gen->eta(), gen->phi() ); 
00233        if(delR<dRmin[jetInd]){
00234          dRmin[jetInd]=delR;           //delta R of match
00235          p4cal[jetInd] = cal->p4();  //Matched Cal 4-vector
00236        }
00237     }
00238     dR.Fill(dRmin[jetInd]);
00239     if(dRmin[jetInd]>0.5)cout << "MC5 Match Warning: dR=" <<dRmin[jetInd]<<", GenPt="<<p4gen[jetInd].Pt()<<", CalPt="<<p4cal[jetInd].Pt()<<endl;
00240     jetInd++;    
00241   }
00242   //Fill Resp vs Pt profile histogram with response of two leading gen jets
00243   for( jetInd=0; jetInd<numJets; ++jetInd ){
00244     if(fabs(p4gen[jetInd].eta())<1.){
00245       respVsPt.Fill(p4gen[jetInd].Pt(), p4cal[jetInd].Pt()/p4gen[jetInd].Pt() );
00246     }
00247   }
00248 
00249   //Find the Corrected CaloJets that match the two uncorrected CaloJets
00250   evt.getByLabel( "corJetMcone5", caloJets );
00251   for( jetInd=0; jetInd<numJets; ++jetInd ){
00252     bool found=kFALSE;
00253     for( CaloJetCollection::const_iterator cor = caloJets->begin(); cor != caloJets->end() && !found; ++ cor ) { //all corrected CaloJets
00254        double delR = deltaR( cor->eta(), cor->phi(), p4cal[jetInd].eta(),  p4cal[jetInd].phi()); 
00255        if(delR<0.01){
00256          dRmin[jetInd]=delR;           //delta R of match
00257          p4cor[jetInd] = cor->p4();  //Matched Cal 4-vector
00258          found=kTRUE;
00259          dRcor.Fill(dRmin[jetInd]);
00260        }
00261     }
00262     if(!found)cout << "Warning: corrected jet not found. jetInd=" << jetInd << endl;
00263   }
00264   //Fill Resp vs Pt profile histogram with corrected response of two leading gen jets
00265   for( jetInd=0; jetInd<numJets; ++jetInd ){
00266     if(fabs(p4gen[jetInd].eta())<1.){
00267      corRespVsPt.Fill(p4gen[jetInd].Pt(), p4cor[jetInd].Pt()/p4gen[jetInd].Pt() ); 
00268     }
00269   }
00270    
00271 }
00272 
00273 void DijetMass::endJob() {
00274 
00275   //Write out the histogram file.
00276   m_file->Write(); 
00277 
00278 }
00279 #include "FWCore/Framework/interface/MakerMacros.h"
00280 DEFINE_FWK_MODULE(DijetMass);