CMSSW Reference Manual

00001 // Producer for validation histograms for CaloJet objects
00002 // F. Ratnikov, Sept. 7, 2006
00003 // $Id: CaloJetTester.cc,v 1.8 2008/03/04 19:52:42 fedor Exp $
00004 
00005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00006 #include "FWCore/Framework/interface/Event.h"
00007 #include "DataFormats/Common/interface/Handle.h"
00008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00009 #include "FWCore/ServiceRegistry/interface/Service.h"
00010 
00011 #include "DataFormats/Math/interface/deltaR.h"
00012 
00013 #include "DQMServices/Core/interface/DQMStore.h"
00014 #include "DQMServices/Core/interface/MonitorElement.h"
00015 
00016 #include "DataFormats/JetReco/interface/CaloJet.h"
00017 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
00018 #include "DataFormats/JetReco/interface/GenJet.h"
00019 #include "DataFormats/JetReco/interface/GenJetCollection.h"
00020 
00021 #include "RecoJets/JetAlgorithms/interface/JetMatchingTools.h"
00022 
00023 #include "CaloJetTester.h"
00024 
00025 using namespace edm;
00026 using namespace reco;
00027 using namespace std;
00028 
00029 namespace {
00030   bool is_B (const reco::Jet& fJet) {return fabs (fJet.eta()) < 1.4;}
00031   bool is_E (const reco::Jet& fJet) {return fabs (fJet.eta()) >= 1.4 && fabs (fJet.eta()) < 3.;}
00032   bool is_F (const reco::Jet& fJet) {return fabs (fJet.eta()) >= 3.;}
00033 }
00034 
00035 CaloJetTester::CaloJetTester(const edm::ParameterSet& iConfig)
00036   : mInputCollection (iConfig.getParameter<edm::InputTag>( "src" )),
00037     mInputGenCollection (iConfig.getParameter<edm::InputTag>( "srcGen" )),
00038     mOutputFile (iConfig.getUntrackedParameter<string>("outputFile", "")),
00039     mMatchGenPtThreshold (iConfig.getParameter<double>("genPtThreshold")),
00040     mGenEnergyFractionThreshold (iConfig.getParameter<double>("genEnergyFractionThreshold")),
00041     mReverseEnergyFractionThreshold (iConfig.getParameter<double>("reverseEnergyFractionThreshold")),
00042     mRThreshold (iConfig.getParameter<double>("RThreshold"))
00043 {
00044   mEta = mPhi = mE = mP = mPt = mMass = mConstituents
00045     = mEtaFirst = mPhiFirst = mEFirst = mPtFirst 
00046     = mMaxEInEmTowers = mMaxEInHadTowers 
00047     = mHadEnergyInHO = mHadEnergyInHB = mHadEnergyInHF = mHadEnergyInHE 
00048     = mEmEnergyInEB = mEmEnergyInEE = mEmEnergyInHF 
00049     = mEnergyFractionHadronic = mEnergyFractionEm 
00050     = mN90
00051     = mAllGenJetsPt = mMatchedGenJetsPt = mAllGenJetsEta = mMatchedGenJetsEta 
00052     = mGenJetMatchEnergyFraction = mReverseMatchEnergyFraction = mRMatch
00053     = mDeltaEta = mDeltaPhi = mEScale = mDeltaE
00054     = 0;
00055   
00056   DQMStore* dbe = &*edm::Service<DQMStore>();
00057   if (dbe) {
00058     dbe->setCurrentFolder("RecoJetsV/CaloJetTask_" + mInputCollection.label());
00059     mEta = dbe->book1D("Eta", "Eta", 100, -5, 5); 
00060     mPhi = dbe->book1D("Phi", "Phi", 70, -3.5, 3.5); 
00061     mE = dbe->book1D("E", "E", 100, 0, 500); 
00062     mP = dbe->book1D("P", "P", 100, 0, 500); 
00063     mPt = dbe->book1D("Pt", "Pt", 100, 0, 50); 
00064     mMass = dbe->book1D("Mass", "Mass", 100, 0, 25); 
00065     mConstituents = dbe->book1D("Constituents", "# of Constituents", 100, 0, 100); 
00066     //
00067     mEtaFirst = dbe->book1D("EtaFirst", "EtaFirst", 100, -5, 5); 
00068     mPhiFirst = dbe->book1D("PhiFirst", "PhiFirst", 70, -3.5, 3.5); 
00069     mEFirst = dbe->book1D("EFirst", "EFirst", 100, 0, 1000); 
00070     mPtFirst = dbe->book1D("PtFirst", "PtFirst", 100, 0, 500); 
00071     //
00072     mMaxEInEmTowers = dbe->book1D("MaxEInEmTowers", "MaxEInEmTowers", 100, 0, 100); 
00073     mMaxEInHadTowers = dbe->book1D("MaxEInHadTowers", "MaxEInHadTowers", 100, 0, 100); 
00074     mHadEnergyInHO = dbe->book1D("HadEnergyInHO", "HadEnergyInHO", 100, 0, 10); 
00075     mHadEnergyInHB = dbe->book1D("HadEnergyInHB", "HadEnergyInHB", 100, 0, 50); 
00076     mHadEnergyInHF = dbe->book1D("HadEnergyInHF", "HadEnergyInHF", 100, 0, 50); 
00077     mHadEnergyInHE = dbe->book1D("HadEnergyInHE", "HadEnergyInHE", 100, 0, 100); 
00078     mEmEnergyInEB = dbe->book1D("EmEnergyInEB", "EmEnergyInEB", 100, 0, 50); 
00079     mEmEnergyInEE = dbe->book1D("EmEnergyInEE", "EmEnergyInEE", 100, 0, 50); 
00080     mEmEnergyInHF = dbe->book1D("EmEnergyInHF", "EmEnergyInHF", 120, -20, 100); 
00081     mEnergyFractionHadronic = dbe->book1D("EnergyFractionHadronic", "EnergyFractionHadronic", 120, -0.1, 1.1); 
00082     mEnergyFractionEm = dbe->book1D("EnergyFractionEm", "EnergyFractionEm", 120, -0.1, 1.1); 
00083     mN90 = dbe->book1D("N90", "N90", 50, 0, 50); 
00084     //
00085     double log10PtMin = 0.5;
00086     double log10PtMax = 4.;
00087     int log10PtBins = 14;
00088     double etaMin = -5.;
00089     double etaMax = 5.;
00090     int etaBins = 50;
00091     mAllGenJetsPt = dbe->book1D("GenJetLOGpT", "GenJet LOG(pT_gen)", 
00092                                 log10PtBins, log10PtMin, log10PtMax);
00093     mMatchedGenJetsPt = dbe->book1D("MatchedGenJetLOGpT", "MatchedGenJet LOG(pT_gen)", 
00094                                     log10PtBins, log10PtMin, log10PtMax);
00095     mAllGenJetsEta = dbe->book2D("GenJetEta", "GenJet Eta vs LOG(pT_gen)", 
00096                                  log10PtBins, log10PtMin, log10PtMax, etaBins, etaMin, etaMax);
00097     mMatchedGenJetsEta = dbe->book2D("MatchedGenJetEta", "MatchedGenJet Eta vs LOG(pT_gen)", 
00098                                      log10PtBins, log10PtMin, log10PtMax, etaBins, etaMin, etaMax);
00099     mGenJetMatchEnergyFraction  = dbe->book3D("GenJetMatchEnergyFraction", "GenJetMatchEnergyFraction vs LOG(pT_gen) vs eta", 
00100                                               log10PtBins, log10PtMin, log10PtMax, etaBins, etaMin, etaMax, 101, 0, 1.01);
00101     mReverseMatchEnergyFraction  = dbe->book3D("ReverseMatchEnergyFraction", "ReverseMatchEnergyFraction vs LOG(pT_gen) vs eta", 
00102                                                log10PtBins, log10PtMin, log10PtMax, etaBins, etaMin, etaMax, 101, 0, 1.01);
00103     mRMatch  = dbe->book3D("RMatch", "delta(R)(Gen-Calo) vs LOG(pT_gen) vs eta", 
00104                            log10PtBins, log10PtMin, log10PtMax, etaBins, etaMin, etaMax, 60, 0, 3);
00105     mDeltaEta = dbe->book3D("DeltaEta", "DeltaEta vs LOG(pT_gen) vs eta", 
00106                               log10PtBins, log10PtMin, log10PtMax, etaBins, etaMin, etaMax, 100, -1, 1);
00107     mDeltaPhi = dbe->book3D("DeltaPhi", "DeltaPhi vs LOG(pT_gen) vs eta", 
00108                               log10PtBins, log10PtMin, log10PtMax, etaBins, etaMin, etaMax, 100, -1, 1);
00109     mEScale = dbe->book3D("EScale", "EnergyScale vs LOG(pT_gen) vs eta", 
00110                             log10PtBins, log10PtMin, log10PtMax, etaBins, etaMin, etaMax, 100, 0, 2);
00111     mDeltaE = dbe->book3D("DeltaE", "DeltaE vs LOG(pT_gen) vs eta", 
00112                             log10PtBins, log10PtMin, log10PtMax, etaBins, etaMin, etaMax, 2000, -200, 200);
00113   }
00114 
00115   if (mOutputFile.empty ()) {
00116     LogInfo("OutputInfo") << " CaloJet histograms will NOT be saved";
00117   } 
00118   else {
00119     LogInfo("OutputInfo") << " CaloJethistograms will be saved to file:" << mOutputFile;
00120   }
00121 }
00122    
00123 CaloJetTester::~CaloJetTester()
00124 {
00125 }
00126 
00127 void CaloJetTester::beginJob(const edm::EventSetup& c){
00128 }
00129 
00130 void CaloJetTester::endJob() {
00131  if (!mOutputFile.empty() && &*edm::Service<DQMStore>()) edm::Service<DQMStore>()->save (mOutputFile);
00132 }
00133 
00134 
00135 void CaloJetTester::analyze(const edm::Event& mEvent, const edm::EventSetup& mSetup)
00136 {
00137   Handle<CaloJetCollection> caloJets;
00138   mEvent.getByLabel(mInputCollection, caloJets);
00139   CaloJetCollection::const_iterator jet = caloJets->begin ();
00140   int jetIndex = 0;
00141   for (; jet != caloJets->end (); jet++, jetIndex++) {
00142     if (mEta) mEta->Fill (jet->eta());
00143     if (mPhi) mPhi->Fill (jet->phi());
00144     if (mE) mE->Fill (jet->energy());
00145     if (mP) mP->Fill (jet->p());
00146     if (mPt) mPt->Fill (jet->pt());
00147     if (mMass) mMass->Fill (jet->mass());
00148     if (mConstituents) mConstituents->Fill (jet->nConstituents());
00149     if (jet == caloJets->begin ()) { // first jet
00150       if (mEtaFirst) mEtaFirst->Fill (jet->eta());
00151       if (mPhiFirst) mPhiFirst->Fill (jet->phi());
00152       if (mEFirst) mEFirst->Fill (jet->energy());
00153       if (mPtFirst) mPtFirst->Fill (jet->pt());
00154     }
00155     if (mMaxEInEmTowers) mMaxEInEmTowers->Fill (jet->maxEInEmTowers());
00156     if (mMaxEInHadTowers) mMaxEInHadTowers->Fill (jet->maxEInHadTowers());
00157     if (mHadEnergyInHO) mHadEnergyInHO->Fill (jet->hadEnergyInHO());
00158     if (mHadEnergyInHB) mHadEnergyInHB->Fill (jet->hadEnergyInHB());
00159     if (mHadEnergyInHF) mHadEnergyInHF->Fill (jet->hadEnergyInHF());
00160     if (mHadEnergyInHE) mHadEnergyInHE->Fill (jet->hadEnergyInHE());
00161     if (mEmEnergyInEB) mEmEnergyInEB->Fill (jet->emEnergyInEB());
00162     if (mEmEnergyInEE) mEmEnergyInEE->Fill (jet->emEnergyInEE());
00163     if (mEmEnergyInHF) mEmEnergyInHF->Fill (jet->emEnergyInHF());
00164     if (mEnergyFractionHadronic) mEnergyFractionHadronic->Fill (jet->energyFractionHadronic());
00165     if (mEnergyFractionEm) mEnergyFractionEm->Fill (jet->emEnergyFraction());
00166     if (mN90) mN90->Fill (jet->n90());
00167   }
00168 
00169   // now match CaloJets to GenJets
00170   JetMatchingTools jetMatching (mEvent);
00171   if (!(mInputGenCollection.label().empty())) {
00172     Handle<GenJetCollection> genJets;
00173     mEvent.getByLabel(mInputGenCollection, genJets);
00174 
00175     std::vector <std::vector <const reco::GenParticle*> > genJetConstituents (genJets->size());
00176     std::vector <std::vector <const reco::GenParticle*> > caloJetConstituents (caloJets->size());
00177     if (mRThreshold > 0) { 
00178     }
00179     else {
00180       for (unsigned iGenJet = 0; iGenJet < genJets->size(); ++iGenJet) {
00181         genJetConstituents [iGenJet] = jetMatching.getGenParticles ((*genJets) [iGenJet]);
00182       }
00183       
00184       for (unsigned iCaloJet = 0; iCaloJet < caloJets->size(); ++iCaloJet) {
00185         caloJetConstituents [iCaloJet] = jetMatching.getGenParticles ((*caloJets) [iCaloJet], false);
00186       }
00187     }
00188 
00189     for (unsigned iGenJet = 0; iGenJet < genJets->size(); ++iGenJet) {
00190       const GenJet& genJet = (*genJets) [iGenJet];
00191       double genJetPt = genJet.pt();
00192       if (fabs(genJet.eta()) > 5.) continue; // out of detector 
00193       if (genJetPt < mMatchGenPtThreshold) continue; // no low momentum 
00194       double logPtGen = log10 (genJetPt);
00195       mAllGenJetsPt->Fill (logPtGen);
00196       mAllGenJetsEta->Fill (logPtGen, genJet.eta());
00197       if (caloJets->size() <= 0) continue; // no CaloJets - nothing to match
00198       if (mRThreshold > 0) {
00199         unsigned iCaloJetBest = 0;
00200         double deltaRBest = 999.;
00201         for (unsigned iCaloJet = 0; iCaloJet < caloJets->size(); ++iCaloJet) {
00202           double dR = deltaR (genJet.eta(), genJet.phi(), (*caloJets) [iCaloJet].eta(), (*caloJets) [iCaloJet].phi());
00203           if (deltaRBest < mRThreshold && dR < mRThreshold && genJet.pt() > 5.) {
00204             std::cout << "Yet another matched jet for GenJet pt=" << genJet.pt()
00205                       << " previous CaloJet pt/dr: " << (*caloJets) [iCaloJetBest].pt() << '/' << deltaRBest
00206                       << " new CaloJet pt/dr: " << (*caloJets) [iCaloJet].pt() << '/' << dR
00207                       << std::endl;
00208           }
00209           if (dR < deltaRBest) {
00210             iCaloJetBest = iCaloJet;
00211             deltaRBest = dR;
00212           }
00213         }
00214         mRMatch->Fill (logPtGen, genJet.eta(), deltaRBest);
00215         if (deltaRBest < mRThreshold) { // Matched
00216           fillMatchHists (genJet, (*caloJets) [iCaloJetBest]);
00217         }
00218       }
00219       else {
00220         unsigned iCaloJetBest = 0;
00221         double energyFractionBest = 0.;
00222         for (unsigned iCaloJet = 0; iCaloJet < caloJets->size(); ++iCaloJet) {
00223           double energyFraction = jetMatching.overlapEnergyFraction (genJetConstituents [iGenJet], 
00224                                                                      caloJetConstituents [iCaloJet]);
00225           if (energyFraction > energyFractionBest) {
00226             iCaloJetBest = iCaloJet;
00227             energyFractionBest = energyFraction;
00228           }
00229         }
00230         mGenJetMatchEnergyFraction->Fill (logPtGen, genJet.eta(), energyFractionBest);
00231         if (energyFractionBest > mGenEnergyFractionThreshold) { // good enough
00232           double reverseEnergyFraction = jetMatching.overlapEnergyFraction (caloJetConstituents [iCaloJetBest], 
00233                                                                             genJetConstituents [iGenJet]);
00234           mReverseMatchEnergyFraction->Fill (logPtGen, genJet.eta(), reverseEnergyFraction);
00235           if (reverseEnergyFraction > mReverseEnergyFractionThreshold) { // Matched
00236             fillMatchHists (genJet, (*caloJets) [iCaloJetBest]);
00237           }
00238         }
00239       }
00240     }
00241   }
00242 }
00243 
00244 void CaloJetTester::fillMatchHists (const reco::GenJet& fGenJet, const reco::CaloJet& fCaloJet) {
00245   double logPtGen = log10 (fGenJet.pt());
00246   mMatchedGenJetsPt->Fill (logPtGen);
00247   mMatchedGenJetsEta->Fill (logPtGen, fGenJet.eta());
00248   mDeltaEta->Fill (logPtGen, fGenJet.eta(), fCaloJet.eta()-fGenJet.eta());
00249   mDeltaPhi->Fill (logPtGen, fGenJet.eta(), fCaloJet.phi()-fGenJet.phi());
00250   mEScale->Fill (logPtGen, fGenJet.eta(), fCaloJet.energy()/fGenJet.energy());
00251   mDeltaE->Fill (logPtGen, fGenJet.eta(), fCaloJet.energy()-fGenJet.energy());
00252 }