Reference Manual

Go to the documentation of this file.
00001 #include "DQMOffline/JetMET/interface/CaloTowerAnalyzer.h"
00002 // author: Bobby Scurlock, University of Florida
00003 // first version 12/18/2006
00004 // modified: Mike Schmitt
00005 // date: 02.28.2007
00006 // note: code rewrite
00007 
00008 #include "FWCore/Framework/interface/MakerMacros.h"
00009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00010 #include "FWCore/ServiceRegistry/interface/Service.h"
00011 
00012 #include "FWCore/PluginManager/interface/ModuleDef.h"
00013 #include "FWCore/Framework/interface/Event.h"
00014 #include "FWCore/Framework/interface/ESHandle.h"
00015 #include "DataFormats/Common/interface/Handle.h"
00016 #include "FWCore/Framework/interface/EventSetup.h"
00017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00018 
00019 #include "DataFormats/JetReco/interface/CaloJet.h"
00020 #include "DataFormats/METReco/interface/CaloMET.h"
00021 #include "DataFormats/METReco/interface/GenMET.h"
00022 #include "DataFormats/METReco/interface/CaloMETCollection.h"
00023 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
00024 #include "DataFormats/METReco/interface/GenMETCollection.h"
00025 #include "DataFormats/METReco/interface/CaloMETCollection.h"
00026 
00027 #include "DataFormats/Math/interface/LorentzVector.h"
00028 #include "RecoMET/METAlgorithms/interface/CaloSpecificAlgo.h"
00029 
00030 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
00031 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
00032 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
00033 #include "DataFormats/Common/interface/TriggerResults.h"
00034 
00035 #include "FWCore/Common/interface/TriggerNames.h"
00036 
00037 #include "DataFormats/DetId/interface/DetId.h"
00038 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
00039 
00040 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00041 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00042 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00043 
00044 //#include "Geometry/Vector/interface/GlobalPoint.h"
00045 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00046 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00047 
00048 #include <vector>
00049 #include <utility>
00050 #include <ostream>
00051 #include <fstream>
00052 #include <string>
00053 #include <algorithm>
00054 #include <cmath>
00055 #include <memory>
00056 #include <TLorentzVector.h>
00057 #include "DQMServices/Core/interface/DQMStore.h"
00058 
00059 using namespace reco;
00060 
00061 CaloTowerAnalyzer::CaloTowerAnalyzer(const edm::ParameterSet & iConfig)
00062 {
00063 
00064   caloTowersLabel_            = iConfig.getParameter<edm::InputTag>("CaloTowersLabel");
00065   HLTResultsLabel_            = iConfig.getParameter<edm::InputTag>("HLTResultsLabel");  
00066   HBHENoiseFilterResultLabel_ = iConfig.getParameter<edm::InputTag>("HBHENoiseFilterResultLabel");
00067 
00068   if(iConfig.exists("HLTBitLabels"))
00069     HLTBitLabel_         = iConfig.getParameter<std::vector<edm::InputTag> >("HLTBitLabels");
00070   
00071   debug_               = iConfig.getParameter<bool>("Debug");
00072   finebinning_         = iConfig.getUntrackedParameter<bool>("FineBinning"); 
00073   allhist_             = iConfig.getUntrackedParameter<bool>("AllHist"); 
00074   FolderName_          = iConfig.getUntrackedParameter<std::string>("FolderName");
00075 
00076   hltselection_        = iConfig.getUntrackedParameter<bool>("HLTSelection"); 
00077   
00078 }
00079 
00080 
00081 void CaloTowerAnalyzer::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup)
00082 {
00083   Nevents = 0;
00084   // get ahold of back-end interface
00085   dbe_ = edm::Service<DQMStore>().operator->();
00086 
00087   if (dbe_) {
00088  
00089     dbe_->setCurrentFolder(FolderName_); 
00090 
00091     //Store number of events which pass each HLT bit 
00092     for(unsigned int i = 0 ; i < HLTBitLabel_.size() ; i++ )
00093       {
00094         if( HLTBitLabel_[i].label().size() )
00095           {
00096             hCT_NEvents_HLT.push_back( dbe_->book1D("METTask_CT_"+HLTBitLabel_[i].label(),HLTBitLabel_[i].label(),2,-0.5,1.5) );
00097           }
00098       }
00099     
00100     //--Store number of events used
00101     hCT_Nevents          = dbe_->book1D("METTask_CT_Nevents","",1,0,1);  
00102     //--Data integrated over all events and stored by CaloTower(ieta,iphi) 
00103     hCT_et_ieta_iphi          = dbe_->book2D("METTask_CT_et_ieta_iphi","",83,-41,42, 72,1,73);  
00104     hCT_et_ieta_iphi->getTH2F()->SetOption("colz");
00105     hCT_et_ieta_iphi->setAxisTitle("ieta",1);
00106     hCT_et_ieta_iphi->setAxisTitle("ephi",2);
00107 
00108     hCT_emEt_ieta_iphi        = dbe_->book2D("METTask_CT_emEt_ieta_iphi","",83,-41,42, 72,1,73);  
00109     hCT_emEt_ieta_iphi->getTH2F()->SetOption("colz");
00110     hCT_emEt_ieta_iphi->setAxisTitle("ieta",1);
00111     hCT_emEt_ieta_iphi->setAxisTitle("ephi",2);
00112     hCT_hadEt_ieta_iphi       = dbe_->book2D("METTask_CT_hadEt_ieta_iphi","",83,-41,42, 72,1,73);  
00113     hCT_hadEt_ieta_iphi->getTH2F()->SetOption("colz");
00114     hCT_hadEt_ieta_iphi->setAxisTitle("ieta",1);
00115     hCT_hadEt_ieta_iphi->setAxisTitle("ephi",2);
00116     hCT_outerEt_ieta_iphi = dbe_->book2D("METTask_CT_outerEt_ieta_iphi","",83,-41,42, 72,1,73);  
00117     hCT_outerEt_ieta_iphi->getTH2F()->SetOption("colz");
00118     hCT_outerEt_ieta_iphi->setAxisTitle("ieta",1);
00119     hCT_outerEt_ieta_iphi->setAxisTitle("ephi",2);
00120     hCT_Occ_ieta_iphi         = dbe_->book2D("METTask_CT_Occ_ieta_iphi","",83,-41,42, 72,1,73);  
00121     hCT_Occ_ieta_iphi->getTH2F()->SetOption("colz");
00122     hCT_Occ_ieta_iphi->setAxisTitle("ieta",1);
00123     hCT_Occ_ieta_iphi->setAxisTitle("ephi",2);
00124 
00125     hCT_Occ_EM_Et_ieta_iphi         = dbe_->book2D("METTask_CT_Occ_EM_Et_ieta_iphi","",83,-41,42, 72,1,73);
00126     hCT_Occ_EM_Et_ieta_iphi->getTH2F()->SetOption("colz");
00127     hCT_Occ_EM_Et_ieta_iphi->setAxisTitle("ieta",1);
00128     hCT_Occ_EM_Et_ieta_iphi->setAxisTitle("ephi",2);
00129 
00130     hCT_Occ_HAD_Et_ieta_iphi         = dbe_->book2D("METTask_CT_Occ_HAD_Et_ieta_iphi","",83,-41,42, 72,1,73);
00131     hCT_Occ_HAD_Et_ieta_iphi->getTH2F()->SetOption("colz");
00132     hCT_Occ_HAD_Et_ieta_iphi->setAxisTitle("ieta",1);
00133     hCT_Occ_HAD_Et_ieta_iphi->setAxisTitle("ephi",2);
00134 
00135     hCT_Occ_Outer_Et_ieta_iphi         = dbe_->book2D("METTask_CT_Occ_Outer_Et_ieta_iphi","",83,-41,42, 72,1,73);
00136     hCT_Occ_Outer_Et_ieta_iphi->getTH2F()->SetOption("colz");
00137     hCT_Occ_Outer_Et_ieta_iphi->setAxisTitle("ieta",1);
00138     hCT_Occ_Outer_Et_ieta_iphi->setAxisTitle("ephi",2);
00139 
00140     //--Data over eta-rings
00141 
00142     // CaloTower values
00143     if(allhist_){
00144     if(finebinning_)
00145       {
00146 
00147         hCT_etvsieta          = dbe_->book2D("METTask_CT_etvsieta","", 83,-41,42, 10001,0,1001);  
00148         hCT_Minetvsieta       = dbe_->book2D("METTask_CT_Minetvsieta","", 83,-41,42, 10001,0,1001);  
00149         hCT_Maxetvsieta       = dbe_->book2D("METTask_CT_Maxetvsieta","", 83,-41,42, 10001,0,1001);  
00150         hCT_emEtvsieta        = dbe_->book2D("METTask_CT_emEtvsieta","",83,-41,42, 10001,0,1001);  
00151         hCT_hadEtvsieta       = dbe_->book2D("METTask_CT_hadEtvsieta","",83,-41,42, 10001,0,1001);  
00152         hCT_outerEtvsieta = dbe_->book2D("METTask_CT_outerEtvsieta","",83,-41,42, 10001,0,1001);  
00153         // Integrated over phi
00154 
00155         hCT_Occvsieta         = dbe_->book2D("METTask_CT_Occvsieta","",83,-41,42, 84,0,84);  
00156         hCT_SETvsieta         = dbe_->book2D("METTask_CT_SETvsieta","",83,-41,42, 20001,0,2001);  
00157         hCT_METvsieta         = dbe_->book2D("METTask_CT_METvsieta","",83,-41,42, 20001,0,2001);  
00158         hCT_METPhivsieta      = dbe_->book2D("METTask_CT_METPhivsieta","",83,-41,42, 80,-4,4);  
00159         hCT_MExvsieta         = dbe_->book2D("METTask_CT_MExvsieta","",83,-41,42, 10001,-500,501);  
00160         hCT_MEyvsieta         = dbe_->book2D("METTask_CT_MEyvsieta","",83,-41,42, 10001,-500,501);  
00161       }
00162     else 
00163       {
00164         
00165         if(allhist_){
00166         hCT_etvsieta          = dbe_->book2D("METTask_CT_etvsieta","", 83,-41,42, 200,-0.5,999.5);
00167         hCT_Minetvsieta       = dbe_->book2D("METTask_CT_Minetvsieta","", 83,-41,42, 200,-0.5,999.5);
00168         hCT_Maxetvsieta       = dbe_->book2D("METTask_CT_Maxetvsieta","", 83,-41,42, 200,-0.5,999.5);
00169         hCT_emEtvsieta        = dbe_->book2D("METTask_CT_emEtvsieta","",83,-41,42, 200,-0.5,999.5);
00170         hCT_hadEtvsieta       = dbe_->book2D("METTask_CT_hadEtvsieta","",83,-41,42, 200,-0.5,999.5);
00171         hCT_outerEtvsieta = dbe_->book2D("METTask_CT_outerEtvsieta","",83,-41,42, 80,-0.5,399.5);
00172         // Integrated over phi
00173         }
00174 
00175         hCT_Occvsieta         = dbe_->book2D("METTask_CT_Occvsieta","",83,-41,42, 73,-0.5,72.5);
00176         hCT_SETvsieta         = dbe_->book2D("METTask_CT_SETvsieta","",83,-41,42, 200,-0.5,1999.5);
00177         hCT_METvsieta         = dbe_->book2D("METTask_CT_METvsieta","",83,-41,42, 200,-0.5,1999.5);
00178         hCT_METPhivsieta      = dbe_->book2D("METTask_CT_METPhivsieta","",83,-41,42, 80,-4,4);
00179         hCT_MExvsieta         = dbe_->book2D("METTask_CT_MExvsieta","",83,-41,42, 100,-499.5,499.5);
00180         hCT_MEyvsieta         = dbe_->book2D("METTask_CT_MEyvsieta","",83,-41,42, 100,-499.5,499.5);
00181         
00182       }
00183     } // allhist
00184   }
00185 }
00186 
00187 void CaloTowerAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00188 {
00189   // Get HLT Results
00190   edm::Handle<edm::TriggerResults> TheHLTResults;
00191   iEvent.getByLabel( HLTResultsLabel_ , TheHLTResults);
00192 
00193   bool EventPasses = true;
00194 
00195   // Make sure handle is valid
00196   if( TheHLTResults.isValid() && hltselection_ )
00197     {
00198       //Get HLT Names
00199       const edm::TriggerNames & TheTriggerNames = iEvent.triggerNames(*TheHLTResults);
00200       
00201       for( unsigned int index = 0 ; index < HLTBitLabel_.size(); index++)
00202         {
00203           if( HLTBitLabel_[index].label().size() )
00204             {
00205               //Change the default value since HLT requirement has been issued by the user
00206               if( index == 0 ) EventPasses = false; 
00207               //Get the HLT bit and check to make sure it is valid
00208               unsigned int bit = TheTriggerNames.triggerIndex( HLTBitLabel_[index].label().c_str());
00209               if( bit < TheHLTResults->size() )
00210                 {
00211                   //If any of the HLT names given by the user accept, then the event passes
00212                   if( TheHLTResults->accept( bit ) && !TheHLTResults->error( bit ) )
00213                     {
00214                       EventPasses = true;
00215                       hCT_NEvents_HLT[index]->Fill(1);
00216                     }  
00217                   else 
00218                     hCT_NEvents_HLT[index]->Fill(0);
00219                 }
00220               else
00221                 {
00222                   edm::LogInfo("OutputInfo") 
00223                     << "The HLT Trigger Name : " << HLTBitLabel_[index].label() << " is not valid for this trigger table " << std::endl;
00224                 }
00225             }
00226         }
00227     }
00228 
00229   if( !EventPasses && hltselection_ ) 
00230     return;
00231   
00232   //----------GREG & CHRIS' idea---///
00233   float ETTowerMin = -1; //GeV
00234   float METRingMin = -2; // GeV
00235   
00236   Nevents++;
00237   hCT_Nevents->Fill(0);
00238 
00239   // ==========================================================
00240   // Retrieve!
00241   // ==========================================================
00242 
00243   edm::Handle<edm::View<Candidate> > towers;
00244   iEvent.getByLabel(caloTowersLabel_, towers);
00245 
00246   if( (!towers.isValid())) {
00247     edm::LogInfo("")<<"CaloTowers "<< caloTowersLabel_<<" not found!"<<std::endl;
00248     return;
00249   }
00250 
00251 
00252   edm::Handle<bool> HBHENoiseFilterResultHandle;
00253   iEvent.getByLabel(HBHENoiseFilterResultLabel_, HBHENoiseFilterResultHandle);
00254   bool HBHENoiseFilterResult = *HBHENoiseFilterResultHandle;
00255   if (!HBHENoiseFilterResultHandle.isValid()) {
00256     LogDebug("") << "CaloTowerAnalyzer: Could not find HBHENoiseFilterResult" << std::endl;
00257   }
00258 
00259 
00260   bool bHcalNoiseFilter = HBHENoiseFilterResult;
00261 
00262   if(!bHcalNoiseFilter) return;
00263 
00264   edm::View<Candidate>::const_iterator towerCand = towers->begin();
00265   
00266   // ==========================================================
00267   // Fill Histograms!
00268   // ==========================================================
00269 
00270   int CTmin_iphi = 99, CTmax_iphi = -99;
00271   int CTmin_ieta = 99, CTmax_ieta = -99;
00272 
00273   TLorentzVector vMET_EtaRing[83];
00274   int ActiveRing[83];
00275   int NActiveTowers[83];
00276   double SET_EtaRing[83];
00277   double MinEt_EtaRing[83];
00278   double MaxEt_EtaRing[83];
00279   for (int i=0;i<83; i++) 
00280     {
00281       ActiveRing[i] = 0;
00282       NActiveTowers[i] = 0;
00283       SET_EtaRing[i] = 0;
00284       MinEt_EtaRing[i] = 0;
00285       MaxEt_EtaRing[i] = 0;
00286     }
00287 
00288   //rcr for (calotower = towerCollection->begin(); calotower != towerCollection->end(); calotower++) {
00289     
00290   for ( ; towerCand != towers->end(); towerCand++)
00291     {
00292       const Candidate* candidate = &(*towerCand);
00293       if (candidate) 
00294         {
00295           const CaloTower* calotower = dynamic_cast<const CaloTower*> (candidate);
00296           if (calotower){
00297           //math::RhoEtaPhiVector Momentum = calotower->momentum();
00298           double Tower_ET = calotower->et();
00299           //double Tower_Energy  = calotower->energy();
00300           //      double Tower_Eta = calotower->eta();
00301           double Tower_Phi = calotower->phi();
00302           //double Tower_EMEnergy = calotower->emEnergy();
00303           //double Tower_HadEnergy = calotower->hadEnergy();
00304           double Tower_OuterEt = calotower->outerEt();
00305           double Tower_EMEt = calotower->emEt();
00306           double Tower_HadEt = calotower->hadEt();
00307           //int Tower_EMLV1 = calotower->emLvl1();
00308           //int Tower_HadLV1 = calotower->hadLv11();
00309           int Tower_ieta = calotower->id().ieta();
00310           int Tower_iphi = calotower->id().iphi();
00311           int EtaRing = 41+Tower_ieta;
00312           ActiveRing[EtaRing] = 1;
00313           NActiveTowers[EtaRing]++;
00314           SET_EtaRing[EtaRing]+=Tower_ET;
00315           TLorentzVector v_;
00316           v_.SetPtEtaPhiE(Tower_ET, 0, Tower_Phi, Tower_ET);
00317           if (Tower_ET>ETTowerMin)
00318             vMET_EtaRing[EtaRing]-=v_;
00319           
00320           // Fill Histograms                                                                                                                                                                                                   
00321           hCT_Occ_ieta_iphi->Fill(Tower_ieta,Tower_iphi);
00322           if (calotower->emEt() > 0 && calotower->emEt() + calotower->hadEt() > 0.3)
00323             hCT_Occ_EM_Et_ieta_iphi->Fill(Tower_ieta,Tower_iphi);
00324           if (calotower->hadEt() > 0 && calotower->emEt() + calotower->hadEt() > 0.3)
00325             hCT_Occ_HAD_Et_ieta_iphi->Fill(Tower_ieta,Tower_iphi);
00326           if (calotower->outerEt() > 0 && calotower->emEt() + calotower->hadEt() > 0.3)
00327             hCT_Occ_Outer_Et_ieta_iphi->Fill(Tower_ieta,Tower_iphi);
00328 
00329           hCT_et_ieta_iphi->Fill(Tower_ieta,Tower_iphi,Tower_ET);
00330           hCT_emEt_ieta_iphi->Fill(Tower_ieta,Tower_iphi,Tower_EMEt);
00331           hCT_hadEt_ieta_iphi->Fill(Tower_ieta,Tower_iphi,Tower_HadEt);
00332           hCT_outerEt_ieta_iphi->Fill(Tower_ieta,Tower_iphi,Tower_OuterEt);
00333 
00334           if (allhist_){
00335           hCT_etvsieta->Fill(Tower_ieta, Tower_ET);
00336           hCT_emEtvsieta->Fill(Tower_ieta, Tower_EMEt);
00337           hCT_hadEtvsieta->Fill(Tower_ieta,Tower_HadEt);
00338           hCT_outerEtvsieta->Fill(Tower_ieta,Tower_OuterEt);
00339           }
00340 
00341           if (Tower_ET > MaxEt_EtaRing[EtaRing])
00342             MaxEt_EtaRing[EtaRing] = Tower_ET;
00343           if (Tower_ET < MinEt_EtaRing[EtaRing] && Tower_ET>0)
00344             MinEt_EtaRing[EtaRing] = Tower_ET;
00345           
00346           
00347           if (Tower_ieta < CTmin_ieta) CTmin_ieta = Tower_ieta;
00348           if (Tower_ieta > CTmax_ieta) CTmax_ieta = Tower_ieta;
00349           if (Tower_iphi < CTmin_iphi) CTmin_iphi = Tower_iphi;
00350           if (Tower_iphi > CTmax_iphi) CTmax_iphi = Tower_iphi;
00351           } //end if (calotower) ..
00352         } // end if(candidate) ...
00353       
00354     } // end loop over towers
00355   
00356       // Fill eta-ring MET quantities
00357   if (allhist_){
00358   for (int iEtaRing=0; iEtaRing<83; iEtaRing++)
00359     { 
00360       hCT_Minetvsieta->Fill(iEtaRing-41, MinEt_EtaRing[iEtaRing]);  
00361       hCT_Maxetvsieta->Fill(iEtaRing-41, MaxEt_EtaRing[iEtaRing]);  
00362       
00363       if (ActiveRing[iEtaRing])
00364         {
00365           if (vMET_EtaRing[iEtaRing].Pt()>METRingMin)
00366             {
00367               hCT_METPhivsieta->Fill(iEtaRing-41, vMET_EtaRing[iEtaRing].Phi());
00368               hCT_MExvsieta->Fill(iEtaRing-41, vMET_EtaRing[iEtaRing].Px());
00369               hCT_MEyvsieta->Fill(iEtaRing-41, vMET_EtaRing[iEtaRing].Py());
00370               hCT_METvsieta->Fill(iEtaRing-41, vMET_EtaRing[iEtaRing].Pt());
00371             }
00372           hCT_SETvsieta->Fill(iEtaRing-41, SET_EtaRing[iEtaRing]);
00373           hCT_Occvsieta->Fill(iEtaRing-41, NActiveTowers[iEtaRing]);
00374         }
00375     } // ietaring
00376   }   // allhist_
00377   
00378   edm::LogInfo("OutputInfo") << "CT ieta range: " << CTmin_ieta << " " << CTmax_ieta;
00379   edm::LogInfo("OutputInfo") << "CT iphi range: " << CTmin_iphi << " " << CTmax_iphi;
00380   
00381 }
00382 
00383 void CaloTowerAnalyzer::endJob()
00384 {
00385 }
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/DQMOffline/JetMET/src/CaloTowerAnalyzer.cc
