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   HcalNoiseSummaryTag_ = iConfig.getParameter<edm::InputTag>("HcalNoiseSummary");
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     //--Data over eta-rings
00125 
00126     // CaloTower values
00127     if(allhist_){
00128     if(finebinning_)
00129       {
00130 
00131         hCT_etvsieta          = dbe_->book2D("METTask_CT_etvsieta","", 83,-41,42, 10001,0,1001);  
00132         hCT_Minetvsieta       = dbe_->book2D("METTask_CT_Minetvsieta","", 83,-41,42, 10001,0,1001);  
00133         hCT_Maxetvsieta       = dbe_->book2D("METTask_CT_Maxetvsieta","", 83,-41,42, 10001,0,1001);  
00134         hCT_emEtvsieta        = dbe_->book2D("METTask_CT_emEtvsieta","",83,-41,42, 10001,0,1001);  
00135         hCT_hadEtvsieta       = dbe_->book2D("METTask_CT_hadEtvsieta","",83,-41,42, 10001,0,1001);  
00136         hCT_outerEtvsieta = dbe_->book2D("METTask_CT_outerEtvsieta","",83,-41,42, 10001,0,1001);  
00137         // Integrated over phi
00138 
00139         hCT_Occvsieta         = dbe_->book2D("METTask_CT_Occvsieta","",83,-41,42, 84,0,84);  
00140         hCT_SETvsieta         = dbe_->book2D("METTask_CT_SETvsieta","",83,-41,42, 20001,0,2001);  
00141         hCT_METvsieta         = dbe_->book2D("METTask_CT_METvsieta","",83,-41,42, 20001,0,2001);  
00142         hCT_METPhivsieta      = dbe_->book2D("METTask_CT_METPhivsieta","",83,-41,42, 80,-4,4);  
00143         hCT_MExvsieta         = dbe_->book2D("METTask_CT_MExvsieta","",83,-41,42, 10001,-500,501);  
00144         hCT_MEyvsieta         = dbe_->book2D("METTask_CT_MEyvsieta","",83,-41,42, 10001,-500,501);  
00145       }
00146     else 
00147       {
00148         
00149         if(allhist_){
00150         hCT_etvsieta          = dbe_->book2D("METTask_CT_etvsieta","", 83,-41,42, 200,-0.5,999.5);
00151         hCT_Minetvsieta       = dbe_->book2D("METTask_CT_Minetvsieta","", 83,-41,42, 200,-0.5,999.5);
00152         hCT_Maxetvsieta       = dbe_->book2D("METTask_CT_Maxetvsieta","", 83,-41,42, 200,-0.5,999.5);
00153         hCT_emEtvsieta        = dbe_->book2D("METTask_CT_emEtvsieta","",83,-41,42, 200,-0.5,999.5);
00154         hCT_hadEtvsieta       = dbe_->book2D("METTask_CT_hadEtvsieta","",83,-41,42, 200,-0.5,999.5);
00155         hCT_outerEtvsieta = dbe_->book2D("METTask_CT_outerEtvsieta","",83,-41,42, 80,-0.5,399.5);
00156         // Integrated over phi
00157         }
00158 
00159         hCT_Occvsieta         = dbe_->book2D("METTask_CT_Occvsieta","",83,-41,42, 73,-0.5,72.5);
00160         hCT_SETvsieta         = dbe_->book2D("METTask_CT_SETvsieta","",83,-41,42, 200,-0.5,1999.5);
00161         hCT_METvsieta         = dbe_->book2D("METTask_CT_METvsieta","",83,-41,42, 200,-0.5,1999.5);
00162         hCT_METPhivsieta      = dbe_->book2D("METTask_CT_METPhivsieta","",83,-41,42, 80,-4,4);
00163         hCT_MExvsieta         = dbe_->book2D("METTask_CT_MExvsieta","",83,-41,42, 100,-499.5,499.5);
00164         hCT_MEyvsieta         = dbe_->book2D("METTask_CT_MEyvsieta","",83,-41,42, 100,-499.5,499.5);
00165         
00166       }
00167     } // allhist
00168   }
00169 }
00170 
00171 void CaloTowerAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00172 {
00173   // Get HLT Results
00174   edm::Handle<edm::TriggerResults> TheHLTResults;
00175   iEvent.getByLabel( HLTResultsLabel_ , TheHLTResults);
00176 
00177   bool EventPasses = true;
00178 
00179   // Make sure handle is valid
00180   if( TheHLTResults.isValid() && hltselection_ )
00181     {
00182       //Get HLT Names
00183       const edm::TriggerNames & TheTriggerNames = iEvent.triggerNames(*TheHLTResults);
00184       
00185       for( unsigned int index = 0 ; index < HLTBitLabel_.size(); index++)
00186         {
00187           if( HLTBitLabel_[index].label().size() )
00188             {
00189               //Change the default value since HLT requirement has been issued by the user
00190               if( index == 0 ) EventPasses = false; 
00191               //Get the HLT bit and check to make sure it is valid
00192               unsigned int bit = TheTriggerNames.triggerIndex( HLTBitLabel_[index].label().c_str());
00193               if( bit < TheHLTResults->size() )
00194                 {
00195                   //If any of the HLT names given by the user accept, then the event passes
00196                   if( TheHLTResults->accept( bit ) && !TheHLTResults->error( bit ) )
00197                     {
00198                       EventPasses = true;
00199                       hCT_NEvents_HLT[index]->Fill(1);
00200                     }  
00201                   else 
00202                     hCT_NEvents_HLT[index]->Fill(0);
00203                 }
00204               else
00205                 {
00206                   edm::LogInfo("OutputInfo") 
00207                     << "The HLT Trigger Name : " << HLTBitLabel_[index].label() << " is not valid for this trigger table " << std::endl;
00208                 }
00209             }
00210         }
00211     }
00212 
00213   if( !EventPasses && hltselection_ ) 
00214     return;
00215   
00216   //----------GREG & CHRIS' idea---///
00217   float ETTowerMin = -1; //GeV
00218   float METRingMin = -2; // GeV
00219   
00220   Nevents++;
00221   hCT_Nevents->Fill(0);
00222 
00223   // ==========================================================
00224   // Retrieve!
00225   // ==========================================================
00226 
00227   edm::Handle<edm::View<Candidate> > towers;
00228   iEvent.getByLabel(caloTowersLabel_, towers);
00229 
00230   if( (!towers.isValid())) {
00231     edm::LogInfo("")<<"CaloTowers "<< caloTowersLabel_<<" not found!"<<std::endl;
00232     return;
00233   }
00234 
00235 
00236   edm::Handle<HcalNoiseSummary> HNoiseSummary;
00237   iEvent.getByLabel(HcalNoiseSummaryTag_,HNoiseSummary);
00238   if (!HNoiseSummary.isValid()) {
00239     LogDebug("") << "CaloTowerAnalyzer: Could not find Hcal NoiseSummary product" << std::endl;
00240   }
00241 
00242   bool bHcalNoiseFilter      = HNoiseSummary->passLooseNoiseFilter();
00243   if(!bHcalNoiseFilter) return;
00244 
00245   edm::View<Candidate>::const_iterator towerCand = towers->begin();
00246   
00247   // ==========================================================
00248   // Fill Histograms!
00249   // ==========================================================
00250 
00251   int CTmin_iphi = 99, CTmax_iphi = -99;
00252   int CTmin_ieta = 99, CTmax_ieta = -99;
00253 
00254   TLorentzVector vMET_EtaRing[83];
00255   int ActiveRing[83];
00256   int NActiveTowers[83];
00257   double SET_EtaRing[83];
00258   double MinEt_EtaRing[83];
00259   double MaxEt_EtaRing[83];
00260   for (int i=0;i<83; i++) 
00261     {
00262       ActiveRing[i] = 0;
00263       NActiveTowers[i] = 0;
00264       SET_EtaRing[i] = 0;
00265       MinEt_EtaRing[i] = 0;
00266       MaxEt_EtaRing[i] = 0;
00267     }
00268 
00269   //rcr for (calotower = towerCollection->begin(); calotower != towerCollection->end(); calotower++) {
00270     
00271   for ( ; towerCand != towers->end(); towerCand++)
00272     {
00273       const Candidate* candidate = &(*towerCand);
00274       if (candidate) 
00275         {
00276           const CaloTower* calotower = dynamic_cast<const CaloTower*> (candidate);
00277           if (calotower){
00278           //math::RhoEtaPhiVector Momentum = calotower->momentum();
00279           double Tower_ET = calotower->et();
00280           //double Tower_Energy  = calotower->energy();
00281           //      double Tower_Eta = calotower->eta();
00282           double Tower_Phi = calotower->phi();
00283           //double Tower_EMEnergy = calotower->emEnergy();
00284           //double Tower_HadEnergy = calotower->hadEnergy();
00285           double Tower_OuterEt = calotower->outerEt();
00286           double Tower_EMEt = calotower->emEt();
00287           double Tower_HadEt = calotower->hadEt();
00288           //int Tower_EMLV1 = calotower->emLvl1();
00289           //int Tower_HadLV1 = calotower->hadLv11();
00290           int Tower_ieta = calotower->id().ieta();
00291           int Tower_iphi = calotower->id().iphi();
00292           int EtaRing = 41+Tower_ieta;
00293           ActiveRing[EtaRing] = 1;
00294           NActiveTowers[EtaRing]++;
00295           SET_EtaRing[EtaRing]+=Tower_ET;
00296           TLorentzVector v_;
00297           v_.SetPtEtaPhiE(Tower_ET, 0, Tower_Phi, Tower_ET);
00298           if (Tower_ET>ETTowerMin)
00299             vMET_EtaRing[EtaRing]-=v_;
00300           
00301           // Fill Histograms
00302           hCT_Occ_ieta_iphi->Fill(Tower_ieta,Tower_iphi);
00303           hCT_et_ieta_iphi->Fill(Tower_ieta,Tower_iphi,Tower_ET);
00304           hCT_emEt_ieta_iphi->Fill(Tower_ieta,Tower_iphi,Tower_EMEt);
00305           hCT_hadEt_ieta_iphi->Fill(Tower_ieta,Tower_iphi,Tower_HadEt);
00306           hCT_outerEt_ieta_iphi->Fill(Tower_ieta,Tower_iphi,Tower_OuterEt);
00307 
00308           if (allhist_){
00309           hCT_etvsieta->Fill(Tower_ieta, Tower_ET);
00310           hCT_emEtvsieta->Fill(Tower_ieta, Tower_EMEt);
00311           hCT_hadEtvsieta->Fill(Tower_ieta,Tower_HadEt);
00312           hCT_outerEtvsieta->Fill(Tower_ieta,Tower_OuterEt);
00313           }
00314 
00315           if (Tower_ET > MaxEt_EtaRing[EtaRing])
00316             MaxEt_EtaRing[EtaRing] = Tower_ET;
00317           if (Tower_ET < MinEt_EtaRing[EtaRing] && Tower_ET>0)
00318             MinEt_EtaRing[EtaRing] = Tower_ET;
00319           
00320           
00321           if (Tower_ieta < CTmin_ieta) CTmin_ieta = Tower_ieta;
00322           if (Tower_ieta > CTmax_ieta) CTmax_ieta = Tower_ieta;
00323           if (Tower_iphi < CTmin_iphi) CTmin_iphi = Tower_iphi;
00324           if (Tower_iphi > CTmax_iphi) CTmax_iphi = Tower_iphi;
00325           } //end if (calotower) ..
00326         } // end if(candidate) ...
00327       
00328     } // end loop over towers
00329   
00330       // Fill eta-ring MET quantities
00331   if (allhist_){
00332   for (int iEtaRing=0; iEtaRing<83; iEtaRing++)
00333     { 
00334       hCT_Minetvsieta->Fill(iEtaRing-41, MinEt_EtaRing[iEtaRing]);  
00335       hCT_Maxetvsieta->Fill(iEtaRing-41, MaxEt_EtaRing[iEtaRing]);  
00336       
00337       if (ActiveRing[iEtaRing])
00338         {
00339           if (vMET_EtaRing[iEtaRing].Pt()>METRingMin)
00340             {
00341               hCT_METPhivsieta->Fill(iEtaRing-41, vMET_EtaRing[iEtaRing].Phi());
00342               hCT_MExvsieta->Fill(iEtaRing-41, vMET_EtaRing[iEtaRing].Px());
00343               hCT_MEyvsieta->Fill(iEtaRing-41, vMET_EtaRing[iEtaRing].Py());
00344               hCT_METvsieta->Fill(iEtaRing-41, vMET_EtaRing[iEtaRing].Pt());
00345             }
00346           hCT_SETvsieta->Fill(iEtaRing-41, SET_EtaRing[iEtaRing]);
00347           hCT_Occvsieta->Fill(iEtaRing-41, NActiveTowers[iEtaRing]);
00348         }
00349     } // ietaring
00350   }   // allhist_
00351   
00352   edm::LogInfo("OutputInfo") << "CT ieta range: " << CTmin_ieta << " " << CTmax_ieta;
00353   edm::LogInfo("OutputInfo") << "CT iphi range: " << CTmin_iphi << " " << CTmax_iphi;
00354   
00355 }
00356 
00357 void CaloTowerAnalyzer::endJob()
00358 {
00359 }
/data/refman/pasoursint/CMSSW_4_1_8_patch9/src/DQMOffline/JetMET/src/CaloTowerAnalyzer.cc
