CMSSW Reference Manual

00001 #include "Validation/RecoMET/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 "PluginManager/ModuleDef.h"
00014 #include "FWCore/Framework/interface/Event.h"
00015 #include "FWCore/Framework/interface/ESHandle.h"
00016 #include "DataFormats/Common/interface/Handle.h"
00017 //#include "FWCore/Framework/interface/Handle.h"
00018 #include "FWCore/Framework/interface/EventSetup.h"
00019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00020 
00021 #include "DataFormats/JetReco/interface/CaloJet.h"
00022 #include "DataFormats/METReco/interface/CaloMET.h"
00023 #include "DataFormats/METReco/interface/GenMET.h"
00024 #include "DataFormats/METReco/interface/CaloMETCollection.h"
00025 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
00026 #include "DataFormats/METReco/interface/GenMETCollection.h"
00027 #include "DataFormats/METReco/interface/CaloMETCollection.h"
00028 
00029 #include "DataFormats/Math/interface/LorentzVector.h"
00030 #include "RecoMET/METAlgorithms/interface/CaloSpecificAlgo.h"
00031 
00032 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
00033 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
00034 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
00035 
00036 #include "DataFormats/DetId/interface/DetId.h"
00037 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
00038 
00039 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00040 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00041 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00042 
00043 //#include "Geometry/Vector/interface/GlobalPoint.h"
00044 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00045 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00046 
00047 #include <vector>
00048 #include <utility>
00049 #include <ostream>
00050 #include <fstream>
00051 #include <string>
00052 #include <algorithm>
00053 #include <cmath>
00054 #include <memory>
00055 #include <TLorentzVector.h>
00056 #include "DQMServices/Core/interface/DQMStore.h"
00057 
00058 using namespace reco;
00059 using namespace std;
00060 
00061 CaloTowerAnalyzer::CaloTowerAnalyzer(const edm::ParameterSet & iConfig)
00062 {
00063 
00064   // outputFile_          = iConfig.getUntrackedParameter<std::string>("OutputFile");
00065   geometryFile_        = iConfig.getUntrackedParameter<std::string>("GeometryFile");
00066   caloTowersLabel_     = iConfig.getParameter<edm::InputTag>("CaloTowersLabel");
00067   debug_               = iConfig.getParameter<bool>("Debug");
00068   dumpGeometry_        = iConfig.getParameter<bool>("DumpGeometry");
00069 
00070   //  if (outputFile_.size() > 0)
00071   //  edm::LogInfo("OutputInfo") << " MET/CaloTower Task histograms will be saved to '" << outputFile_.c_str() << "'";
00072   //else edm::LogInfo("OutputInfo") << " MET/CaloTower Task histograms will NOT be saved";
00073 
00074 }
00075 
00076 void CaloTowerAnalyzer::beginJob(const edm::EventSetup& iSetup)
00077 {
00078   Nevents = 0;
00079   // get ahold of back-end interface
00080   dbe_ = edm::Service<DQMStore>().operator->();
00081 
00082   if (dbe_) {
00083     dbe_->setCurrentFolder("RecoMETV/METTask/CaloTowers/geometry");
00084 
00085     me["hCT_ieta_iphi_etaMap"]      = dbe_->book2D("METTask_CT_ieta_iphi_etaMap","",83,-41,42, 72,1,73);
00086     me["hCT_ieta_iphi_phiMap"]      = dbe_->book2D("METTask_CT_ieta_iphi_phiMap","",83,-41,42, 72,1,73);
00087     me["hCT_ieta_detaMap"]          = dbe_->book1D("METTask_CT_ieta_detaMap","", 83, -41, 42);
00088     me["hCT_ieta_dphiMap"]          = dbe_->book1D("METTask_CT_ieta_dphiMap","", 83, -41, 42);
00089 
00090     // Initialize bins for geometry to -999 because z = 0 is a valid entry
00091     for (int i=1; i<=83; i++) {
00092       me["hCT_ieta_detaMap"]->setBinContent(i,-999);
00093       me["hCT_ieta_dphiMap"]->setBinContent(i,-999);
00094       for (int j=1; j<=73; j++) {
00095         me["hCT_ieta_iphi_etaMap"]->setBinContent(i,j,-999);
00096         me["hCT_ieta_iphi_phiMap"]->setBinContent(i,j,-999);
00097       }
00098     }
00099     TString dirName = "RecoMETV/METTask/CaloTowers/";
00100     TString label(caloTowersLabel_.label()); 
00101     if(label=="towerMaker") dirName += "SchemeB";
00102     else if(label=="calotoweroptmaker") dirName += "Optimized";
00103     else dirName += label;
00104     dbe_->setCurrentFolder((string)dirName); 
00105     
00106     //--Store number of events used
00107     me["hCT_Nevents"]          = dbe_->book1D("METTask_CT_Nevents","",1,0,1);  
00108     //--Data integrated over all events and stored by CaloTower(ieta,iphi) 
00109     me["hCT_et_ieta_iphi"]          = dbe_->book2D("METTask_CT_et_ieta_iphi","",83,-41,42, 72,1,73);  
00110     me["hCT_Minet_ieta_iphi"]          = dbe_->book2D("METTask_CT_Minet_ieta_iphi","",83,-41,42, 72,1,73);  
00111     me["hCT_Maxet_ieta_iphi"]          = dbe_->book2D("METTask_CT_Maxet_ieta_iphi","",83,-41,42, 72,1,73);  
00112     for (int i = 1; i<=83; i++)
00113       for (int j = 1; j<=73; j++)
00114         {
00115           me["hCT_Minet_ieta_iphi"]->setBinContent(i,j,14E3);
00116           me["hCT_Maxet_ieta_iphi"]->setBinContent(i,j,-999);
00117         }
00118 
00119     me["hCT_emEt_ieta_iphi"]        = dbe_->book2D("METTask_CT_emEt_ieta_iphi","",83,-41,42, 72,1,73);  
00120     me["hCT_hadEt_ieta_iphi"]       = dbe_->book2D("METTask_CT_hadEt_ieta_iphi","",83,-41,42, 72,1,73);  
00121     me["hCT_energy_ieta_iphi"]      = dbe_->book2D("METTask_CT_energy_ieta_iphi","",83,-41,42, 72,1,73);  
00122     me["hCT_outerEnergy_ieta_iphi"] = dbe_->book2D("METTask_CT_outerEnergy_ieta_iphi","",83,-41,42, 72,1,73);  
00123     me["hCT_hadEnergy_ieta_iphi"]   = dbe_->book2D("METTask_CT_hadEnergy_ieta_iphi","",83,-41,42, 72,1,73);  
00124     me["hCT_emEnergy_ieta_iphi"]    = dbe_->book2D("METTask_CT_emEnergy_ieta_iphi","",83,-41,42, 72,1,73);  
00125     me["hCT_Occ_ieta_iphi"]         = dbe_->book2D("METTask_CT_Occ_ieta_iphi","",83,-41,42, 72,1,73);  
00126     //--Data over eta-rings
00127     // CaloTower values
00128     me["hCT_etvsieta"]          = dbe_->book2D("METTask_CT_etvsieta","", 83,-41,42, 10001,0,1001);  
00129     me["hCT_Minetvsieta"]          = dbe_->book2D("METTask_CT_Minetvsieta","", 83,-41,42, 10001,0,1001);  
00130     me["hCT_Maxetvsieta"]          = dbe_->book2D("METTask_CT_Maxetvsieta","", 83,-41,42, 10001,0,1001);  
00131     me["hCT_emEtvsieta"]        = dbe_->book2D("METTask_CT_emEtvsieta","",83,-41,42, 10001,0,1001);  
00132     me["hCT_hadEtvsieta"]       = dbe_->book2D("METTask_CT_hadEtvsieta","",83,-41,42, 10001,0,1001);  
00133     me["hCT_energyvsieta"]      = dbe_->book2D("METTask_CT_energyvsieta","",83,-41,42, 10001,0,1001);  
00134     me["hCT_outerEnergyvsieta"] = dbe_->book2D("METTask_CT_outerEnergyvsieta","",83,-41,42, 10001,0,1001);  
00135     me["hCT_hadEnergyvsieta"]   = dbe_->book2D("METTask_CT_hadEnergyvsieta","",83,-41,42, 10001,0,1001);  
00136     me["hCT_emEnergyvsieta"]    = dbe_->book2D("METTask_CT_emEnergyvsieta","",83,-41,42, 10001,0,1001);  
00137     // Integrated over phi
00138     me["hCT_Occvsieta"]         = dbe_->book2D("METTask_CT_Occvsieta","",83,-41,42, 84,0,84);  
00139     me["hCT_SETvsieta"]         = dbe_->book2D("METTask_CT_SETvsieta","",83,-41,42, 20001,0,2001);  
00140     me["hCT_METvsieta"]         = dbe_->book2D("METTask_CT_METvsieta","",83,-41,42, 20001,0,2001);  
00141     me["hCT_METPhivsieta"]      = dbe_->book2D("METTask_CT_METPhivsieta","",83,-41,42, 80,-4,4);  
00142     me["hCT_MExvsieta"]         = dbe_->book2D("METTask_CT_MExvsieta","",83,-41,42, 10001,-500,501);  
00143     me["hCT_MEyvsieta"]         = dbe_->book2D("METTask_CT_MEyvsieta","",83,-41,42, 10001,-500,501);  
00144   }
00145 
00146   // Inspect Setup for CaloTower Geometry
00147   FillGeometry(iSetup);
00148 
00149 }
00150 
00151 void CaloTowerAnalyzer::FillGeometry(const edm::EventSetup& iSetup)
00152 {
00153 
00154   // ==========================================================
00155   // Retrieve!
00156   // ==========================================================
00157 
00158   const CaloSubdetectorGeometry* geom;
00159 
00160   try {
00161 
00162     edm::ESHandle<CaloGeometry> pG;
00163     iSetup.get<CaloGeometryRecord>().get(pG);
00164     const CaloGeometry cG = *pG;
00165     geom = cG.getSubdetectorGeometry(DetId::Calo,1);
00166 
00167   } catch (...) {
00168 
00169     edm::LogInfo("OutputInfo") << "Failed to retrieve an Event Setup Handle, Aborting METTask/"
00170                           << "CaloTowerAnalyzer::FillGeometry!"; return;
00171 
00172   }
00173     
00174   // ==========================================================
00175   // Fill Histograms!
00176   // ==========================================================
00177 
00178   vector<DetId> ids = geom->getValidDetIds(DetId::Calo,1);
00179   vector<DetId>::iterator i;
00180 
00181   // Loop Over all CaloTower DetId's
00182   int ndetid = 0;
00183   for (i = ids.begin(); i != ids.end(); i++) {
00184 
00185     ndetid++;
00186 
00187     const CaloCellGeometry* cell = geom->getGeometry(*i);
00188     CaloTowerDetId ctId(i->rawId());
00189     //GlobalPoint p = cell->getPosition();
00190       
00191     int Tower_ieta = ctId.ieta();
00192     int Tower_iphi = ctId.iphi();
00193     double Tower_eta = cell->getPosition().eta();
00194     double Tower_phi = cell->getPosition().phi();
00195       
00196     me["hCT_ieta_iphi_etaMap"]->setBinContent(Tower_ieta+42, Tower_iphi, Tower_eta);
00197     me["hCT_ieta_iphi_phiMap"]->setBinContent(Tower_ieta+42, Tower_iphi, (Tower_phi*180.0/M_PI) );
00198       
00199   } // end loop over DetId's
00200 
00201   // Set the Cell Size for each (ieta, iphi) Bin
00202   double currentLowEdge_eta = 0; //double currentHighEdge_eta = 0;
00203   
00204   for (int ieta=1; ieta<=41 ; ieta++) {
00205 
00206     int ieta_ = 42 + ieta;
00207     double eta = me["hCT_ieta_iphi_etaMap"]->getBinContent(ieta_,3);
00208     double phi = me["hCT_ieta_iphi_phiMap"]->getBinContent(ieta_,3);
00209     double deta = 2.0*(eta-currentLowEdge_eta);
00210     deta = ((float)((int)(1.0E3*deta + 0.5)))/1.0E3;
00211     double dphi = 2.0*phi;
00212     if (ieta==40 || ieta==41) dphi = 20;
00213     if (ieta<=39 && ieta>=21) dphi = 10;
00214     if (ieta<=20) dphi = 5;
00215     // BS: This is WRONG...need to correct overlap 
00216     if (ieta==28) deta = 0.218;
00217     if (ieta==29) deta= 0.096;      
00218     currentLowEdge_eta += deta;
00219 
00220     // BS: This is WRONG...need to correct overlap 
00221     if (ieta==29) currentLowEdge_eta = 2.964;
00222     me["hCT_ieta_detaMap"]->setBinContent(ieta_,deta); // positive rings
00223     me["hCT_ieta_dphiMap"]->setBinContent(ieta_,dphi); // positive rings
00224     me["hCT_ieta_detaMap"]->setBinContent(42-ieta,deta); // negative rings
00225     me["hCT_ieta_dphiMap"]->setBinContent(42-ieta,dphi); // negative rings
00226 
00227   } // end loop over ieta
00228 
00229 }
00230 
00231 void CaloTowerAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00232 {
00233   //----------GREG & CHRIS' idea---///
00234    float ETTowerMin = -1; //GeV
00235    float METRingMin = -2; // GeV
00236 
00237   Nevents++;
00238   me["hCT_Nevents"]->Fill(0);
00239 
00240   // ==========================================================
00241   // Retrieve!
00242   // ==========================================================
00243 
00244   
00245   //rcr// const CaloTowerCollection *towerCollection;
00246 
00247 
00248   /* rcr
00249   edm::Handle<reco::CandidateCollection> to;
00250   iEvent.getByLabel(caloTowersLabel_,to);
00251   if (!to.isValid()) {
00252     edm::LogInfo("OutputInfo") << "Failed to retrieve an Event Handle, Aborting METTask/"
00253                                << "CaloTowerAnalyzer::analyze!"; return;
00254   } else {
00255     const CandidateCollection *towers = (CandidateCollection *)to.product();
00256     reco::CandidateCollection::const_iterator tower = towers->begin();
00257     edm::Ref<CaloTowerCollection> towerRef = tower->get<CaloTowerRef>();
00258     towerCollection = towerRef.product();
00259   }
00260 
00261   rcr */
00262 
00263   edm::Handle<edm::View<Candidate> > towers;
00264   iEvent.getByLabel(caloTowersLabel_, towers);
00265   edm::View<Candidate>::const_iterator towerCand = towers->begin();
00266   
00267   // ==========================================================
00268   // Fill Histograms!
00269   // ==========================================================
00270 
00271   //rcr edm::LogInfo("OutputInfo") << "There are " << towerCollection->size() << " CaloTowers";
00272   //rcr CaloTowerCollection::const_iterator calotower;
00273   
00274   int CTmin_iphi = 99, CTmax_iphi = -99;
00275   int CTmin_ieta = 99, CTmax_ieta = -99;
00276 
00277   TLorentzVector vMET_EtaRing[83];
00278   int ActiveRing[83];
00279   int NActiveTowers[83];
00280   double SET_EtaRing[83];
00281   double MinEt_EtaRing[83];
00282   double MaxEt_EtaRing[83];
00283   for (int i=0;i<83; i++) 
00284     {
00285       ActiveRing[i] = 0;
00286       NActiveTowers[i] = 0;
00287       SET_EtaRing[i] = 0;
00288       MinEt_EtaRing[i] = 0;
00289       MaxEt_EtaRing[i] = 0;
00290     }
00291 
00292   //rcr for (calotower = towerCollection->begin(); calotower != towerCollection->end(); calotower++) {
00293     
00294   for ( ; towerCand != towers->end(); towerCand++)
00295     {
00296       const Candidate* candidate = &(*towerCand);
00297       if (candidate) 
00298         {
00299           const CaloTower* calotower = dynamic_cast<const CaloTower*> (candidate);
00300           if (calotower){
00301           //math::RhoEtaPhiVector Momentum = calotower->momentum();
00302           double Tower_ET = calotower->et();
00303           double Tower_Energy  = calotower->energy();
00304           double Tower_Eta = calotower->eta();
00305           double Tower_Phi = calotower->phi();
00306           double Tower_EMEnergy = calotower->emEnergy();
00307           double Tower_HadEnergy = calotower->hadEnergy();
00308           double Tower_OuterEnergy = calotower->outerEnergy();
00309           double Tower_EMEt = calotower->emEt();
00310           double Tower_HadEt = calotower->hadEt();
00311           //int Tower_EMLV1 = calotower->emLvl1();
00312           //int Tower_HadLV1 = calotower->hadLv11();
00313           int Tower_ieta = calotower->id().ieta();
00314           int Tower_iphi = calotower->id().iphi();
00315           int EtaRing = 41+Tower_ieta;
00316           ActiveRing[EtaRing] = 1;
00317           NActiveTowers[EtaRing]++;
00318           SET_EtaRing[EtaRing]+=Tower_ET;
00319           TLorentzVector v_;
00320           v_.SetPtEtaPhiE(Tower_ET, 0, Tower_Phi, Tower_ET);
00321           if (Tower_ET>ETTowerMin)
00322             vMET_EtaRing[EtaRing]-=v_;
00323           
00324           // Fill Histograms
00325           me["hCT_Occ_ieta_iphi"]->Fill(Tower_ieta,Tower_iphi);
00326           me["hCT_et_ieta_iphi"]->Fill(Tower_ieta,Tower_iphi,Tower_ET);
00327           me["hCT_emEt_ieta_iphi"]->Fill(Tower_ieta,Tower_iphi,Tower_EMEt);
00328           me["hCT_hadEt_ieta_iphi"]->Fill(Tower_ieta,Tower_iphi,Tower_HadEt);
00329           me["hCT_energy_ieta_iphi"]->Fill(Tower_ieta,Tower_iphi,Tower_Energy);
00330           me["hCT_outerEnergy_ieta_iphi"]->Fill(Tower_ieta,Tower_iphi,Tower_OuterEnergy);
00331           me["hCT_hadEnergy_ieta_iphi"]->Fill(Tower_ieta,Tower_iphi,Tower_HadEnergy);
00332           me["hCT_emEnergy_ieta_iphi"]->Fill(Tower_ieta,Tower_iphi,Tower_EMEnergy);
00333           
00334           me["hCT_etvsieta"]->Fill(Tower_ieta, Tower_ET);
00335           me["hCT_emEtvsieta"]->Fill(Tower_ieta, Tower_EMEt);
00336           me["hCT_hadEtvsieta"]->Fill(Tower_ieta,Tower_HadEt);
00337           me["hCT_energyvsieta"]->Fill(Tower_ieta,Tower_Energy);
00338           me["hCT_outerEnergyvsieta"]->Fill(Tower_ieta,Tower_OuterEnergy);
00339           me["hCT_hadEnergyvsieta"]->Fill(Tower_ieta ,Tower_HadEnergy);
00340           me["hCT_emEnergyvsieta"]->Fill(Tower_ieta,Tower_EMEnergy);
00341 
00342           if (Tower_ET > MaxEt_EtaRing[EtaRing])
00343             MaxEt_EtaRing[EtaRing] = Tower_ET;
00344           if (Tower_ET < MinEt_EtaRing[EtaRing] && Tower_ET>0)
00345             MinEt_EtaRing[EtaRing] = Tower_ET;
00346           
00347           
00348           if (Tower_ieta < CTmin_ieta) CTmin_ieta = Tower_ieta;
00349           if (Tower_ieta > CTmax_ieta) CTmax_ieta = Tower_ieta;
00350           if (Tower_iphi < CTmin_iphi) CTmin_iphi = Tower_iphi;
00351           if (Tower_iphi > CTmax_iphi) CTmax_iphi = Tower_iphi;
00352           } //end if (calotower) ..
00353         } // end if(candidate) ...
00354       
00355     } // end loop over towers
00356   
00357       // Fill eta-ring MET quantities
00358   for (int iEtaRing=0; iEtaRing<83; iEtaRing++)
00359     { 
00360       me["hCT_Minetvsieta"]->Fill(iEtaRing-41, MinEt_EtaRing[iEtaRing]);  
00361       me["hCT_Maxetvsieta"]->Fill(iEtaRing-41, MaxEt_EtaRing[iEtaRing]);  
00362       
00363       if (ActiveRing[iEtaRing])
00364         {
00365           if (vMET_EtaRing[iEtaRing].Pt()>METRingMin)
00366             {
00367               me["hCT_METPhivsieta"]->Fill(iEtaRing-41, vMET_EtaRing[iEtaRing].Phi());
00368               me["hCT_MExvsieta"]->Fill(iEtaRing-41, vMET_EtaRing[iEtaRing].Px());
00369               me["hCT_MEyvsieta"]->Fill(iEtaRing-41, vMET_EtaRing[iEtaRing].Py());
00370               me["hCT_METvsieta"]->Fill(iEtaRing-41, vMET_EtaRing[iEtaRing].Pt());
00371             }
00372           me["hCT_SETvsieta"]->Fill(iEtaRing-41, SET_EtaRing[iEtaRing]);
00373           me["hCT_Occvsieta"]->Fill(iEtaRing-41, NActiveTowers[iEtaRing]);
00374         }
00375     }
00376   
00377   edm::LogInfo("OutputInfo") << "CT ieta range: " << CTmin_ieta << " " << CTmax_ieta;
00378   edm::LogInfo("OutputInfo") << "CT iphi range: " << CTmin_iphi << " " << CTmax_iphi;
00379   
00380 }
00381 
00382 
00383 void CaloTowerAnalyzer::DumpGeometry()
00384 {
00385   
00386   ofstream dump(geometryFile_.c_str());
00387   
00388   dump << "Tower Definitions: " << endl << endl;
00389   
00390   dump.width(15); dump << left << "ieta bin";
00391   dump.width(15); dump << left << "Eta";
00392   //dump.width(15); dump << left << "Phi";
00393   dump.width(15); dump << left << "dEta";
00394   dump.width(15); dump << left << "dPhi" << endl;
00395 
00396   int max_ieta_bin = me["hCT_ieta_iphi_etaMap"]->getNbinsX();
00397   for (int i = 1; i <= max_ieta_bin; i++) {
00398 
00399     dump.width(15); dump << left << i;
00400     dump.width(15); dump << left << me["hCT_ieta_iphi_etaMap"]->getBinContent(i,1);
00401     //dump.width(15); dump << left << me["hCT_ieta_iphi_phiMap"]->getBinContent(i,1);
00402     dump.width(15); dump << left << me["hCT_ieta_detaMap"]->getBinContent(i);
00403     dump.width(15); dump << left << me["hCT_ieta_dphiMap"]->getBinContent(i) << endl;
00404     
00405   }
00406   
00407   dump.close();
00408   
00409 }
00410 
00411 void CaloTowerAnalyzer::endJob()
00412 {
00413   // Store the DAQ Histograms
00414   //if (outputFile_.size() > 0 && dbe_)
00415   //  dbe_->save(outputFile_);
00416 
00417   // Dump Geometry Info to a File
00418   if (dumpGeometry_); DumpGeometry();
00419 
00420 }
CaloTowerAnalyzer.cc
