CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CaloTowerAnalyzer.cc
Go to the documentation of this file.
2 // author: Bobby Scurlock, University of Florida
3 // first version 12/18/2006
4 // modified: Mike Schmitt
5 // date: 02.28.2007
6 // note: code rewrite
7 
11 
18 
26 
29 
34 
36 
39 
43 
44 //#include "Geometry/Vector/interface/GlobalPoint.h"
47 
48 #include <vector>
49 #include <utility>
50 #include <ostream>
51 #include <fstream>
52 #include <string>
53 #include <algorithm>
54 #include <cmath>
55 #include <memory>
56 #include <TLorentzVector.h>
58 
59 using namespace reco;
60 
62 {
63 
64  caloTowersLabel_ = iConfig.getParameter<edm::InputTag>("CaloTowersLabel");
65  HLTResultsLabel_ = iConfig.getParameter<edm::InputTag>("HLTResultsLabel");
66  HBHENoiseFilterResultLabel_ = iConfig.getParameter<edm::InputTag>("HBHENoiseFilterResultLabel");
67 
68  if(iConfig.exists("HLTBitLabels"))
69  HLTBitLabel_ = iConfig.getParameter<std::vector<edm::InputTag> >("HLTBitLabels");
70 
71  debug_ = iConfig.getParameter<bool>("Debug");
72  finebinning_ = iConfig.getUntrackedParameter<bool>("FineBinning");
73  allhist_ = iConfig.getUntrackedParameter<bool>("AllHist");
74  FolderName_ = iConfig.getUntrackedParameter<std::string>("FolderName");
75 
76  hltselection_ = iConfig.getUntrackedParameter<bool>("HLTSelection");
77 
78 }
79 
80 
81 void CaloTowerAnalyzer::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup)
82 {
83  Nevents = 0;
84  // get ahold of back-end interface
86 
87  if (dbe_) {
88 
89  dbe_->setCurrentFolder(FolderName_);
90 
91  //Store number of events which pass each HLT bit
92  for(unsigned int i = 0 ; i < HLTBitLabel_.size() ; i++ )
93  {
94  if( HLTBitLabel_[i].label().size() )
95  {
96  hCT_NEvents_HLT.push_back( dbe_->book1D("METTask_CT_"+HLTBitLabel_[i].label(),HLTBitLabel_[i].label(),2,-0.5,1.5) );
97  }
98  }
99 
100  //--Store number of events used
101  hCT_Nevents = dbe_->book1D("METTask_CT_Nevents","",1,0,1);
102  //--Data integrated over all events and stored by CaloTower(ieta,iphi)
103  hCT_et_ieta_iphi = dbe_->book2D("METTask_CT_et_ieta_iphi","",83,-41,42, 72,1,73);
104  hCT_et_ieta_iphi->getTH2F()->SetOption("colz");
105  hCT_et_ieta_iphi->setAxisTitle("ieta",1);
106  hCT_et_ieta_iphi->setAxisTitle("ephi",2);
107 
108  hCT_emEt_ieta_iphi = dbe_->book2D("METTask_CT_emEt_ieta_iphi","",83,-41,42, 72,1,73);
109  hCT_emEt_ieta_iphi->getTH2F()->SetOption("colz");
110  hCT_emEt_ieta_iphi->setAxisTitle("ieta",1);
111  hCT_emEt_ieta_iphi->setAxisTitle("ephi",2);
112  hCT_hadEt_ieta_iphi = dbe_->book2D("METTask_CT_hadEt_ieta_iphi","",83,-41,42, 72,1,73);
113  hCT_hadEt_ieta_iphi->getTH2F()->SetOption("colz");
114  hCT_hadEt_ieta_iphi->setAxisTitle("ieta",1);
115  hCT_hadEt_ieta_iphi->setAxisTitle("ephi",2);
116  hCT_outerEt_ieta_iphi = dbe_->book2D("METTask_CT_outerEt_ieta_iphi","",83,-41,42, 72,1,73);
117  hCT_outerEt_ieta_iphi->getTH2F()->SetOption("colz");
118  hCT_outerEt_ieta_iphi->setAxisTitle("ieta",1);
119  hCT_outerEt_ieta_iphi->setAxisTitle("ephi",2);
120  hCT_Occ_ieta_iphi = dbe_->book2D("METTask_CT_Occ_ieta_iphi","",83,-41,42, 72,1,73);
121  hCT_Occ_ieta_iphi->getTH2F()->SetOption("colz");
122  hCT_Occ_ieta_iphi->setAxisTitle("ieta",1);
123  hCT_Occ_ieta_iphi->setAxisTitle("ephi",2);
124 
125  hCT_Occ_EM_Et_ieta_iphi = dbe_->book2D("METTask_CT_Occ_EM_Et_ieta_iphi","",83,-41,42, 72,1,73);
126  hCT_Occ_EM_Et_ieta_iphi->getTH2F()->SetOption("colz");
127  hCT_Occ_EM_Et_ieta_iphi->setAxisTitle("ieta",1);
128  hCT_Occ_EM_Et_ieta_iphi->setAxisTitle("ephi",2);
129 
130  hCT_Occ_HAD_Et_ieta_iphi = dbe_->book2D("METTask_CT_Occ_HAD_Et_ieta_iphi","",83,-41,42, 72,1,73);
131  hCT_Occ_HAD_Et_ieta_iphi->getTH2F()->SetOption("colz");
132  hCT_Occ_HAD_Et_ieta_iphi->setAxisTitle("ieta",1);
133  hCT_Occ_HAD_Et_ieta_iphi->setAxisTitle("ephi",2);
134 
135  hCT_Occ_Outer_Et_ieta_iphi = dbe_->book2D("METTask_CT_Occ_Outer_Et_ieta_iphi","",83,-41,42, 72,1,73);
136  hCT_Occ_Outer_Et_ieta_iphi->getTH2F()->SetOption("colz");
137  hCT_Occ_Outer_Et_ieta_iphi->setAxisTitle("ieta",1);
138  hCT_Occ_Outer_Et_ieta_iphi->setAxisTitle("ephi",2);
139 
140  //--Data over eta-rings
141 
142  // CaloTower values
143  if(allhist_){
144  if(finebinning_)
145  {
146 
147  hCT_etvsieta = dbe_->book2D("METTask_CT_etvsieta","", 83,-41,42, 10001,0,1001);
148  hCT_Minetvsieta = dbe_->book2D("METTask_CT_Minetvsieta","", 83,-41,42, 10001,0,1001);
149  hCT_Maxetvsieta = dbe_->book2D("METTask_CT_Maxetvsieta","", 83,-41,42, 10001,0,1001);
150  hCT_emEtvsieta = dbe_->book2D("METTask_CT_emEtvsieta","",83,-41,42, 10001,0,1001);
151  hCT_hadEtvsieta = dbe_->book2D("METTask_CT_hadEtvsieta","",83,-41,42, 10001,0,1001);
152  hCT_outerEtvsieta = dbe_->book2D("METTask_CT_outerEtvsieta","",83,-41,42, 10001,0,1001);
153  // Integrated over phi
154 
155  hCT_Occvsieta = dbe_->book2D("METTask_CT_Occvsieta","",83,-41,42, 84,0,84);
156  hCT_SETvsieta = dbe_->book2D("METTask_CT_SETvsieta","",83,-41,42, 20001,0,2001);
157  hCT_METvsieta = dbe_->book2D("METTask_CT_METvsieta","",83,-41,42, 20001,0,2001);
158  hCT_METPhivsieta = dbe_->book2D("METTask_CT_METPhivsieta","",83,-41,42, 80,-4,4);
159  hCT_MExvsieta = dbe_->book2D("METTask_CT_MExvsieta","",83,-41,42, 10001,-500,501);
160  hCT_MEyvsieta = dbe_->book2D("METTask_CT_MEyvsieta","",83,-41,42, 10001,-500,501);
161  }
162  else
163  {
164 
165  if(allhist_){
166  hCT_etvsieta = dbe_->book2D("METTask_CT_etvsieta","", 83,-41,42, 200,-0.5,999.5);
167  hCT_Minetvsieta = dbe_->book2D("METTask_CT_Minetvsieta","", 83,-41,42, 200,-0.5,999.5);
168  hCT_Maxetvsieta = dbe_->book2D("METTask_CT_Maxetvsieta","", 83,-41,42, 200,-0.5,999.5);
169  hCT_emEtvsieta = dbe_->book2D("METTask_CT_emEtvsieta","",83,-41,42, 200,-0.5,999.5);
170  hCT_hadEtvsieta = dbe_->book2D("METTask_CT_hadEtvsieta","",83,-41,42, 200,-0.5,999.5);
171  hCT_outerEtvsieta = dbe_->book2D("METTask_CT_outerEtvsieta","",83,-41,42, 80,-0.5,399.5);
172  // Integrated over phi
173  }
174 
175  hCT_Occvsieta = dbe_->book2D("METTask_CT_Occvsieta","",83,-41,42, 73,-0.5,72.5);
176  hCT_SETvsieta = dbe_->book2D("METTask_CT_SETvsieta","",83,-41,42, 200,-0.5,1999.5);
177  hCT_METvsieta = dbe_->book2D("METTask_CT_METvsieta","",83,-41,42, 200,-0.5,1999.5);
178  hCT_METPhivsieta = dbe_->book2D("METTask_CT_METPhivsieta","",83,-41,42, 80,-4,4);
179  hCT_MExvsieta = dbe_->book2D("METTask_CT_MExvsieta","",83,-41,42, 100,-499.5,499.5);
180  hCT_MEyvsieta = dbe_->book2D("METTask_CT_MEyvsieta","",83,-41,42, 100,-499.5,499.5);
181 
182  }
183  } // allhist
184  }
185 }
186 
188 {
189  // Get HLT Results
190  edm::Handle<edm::TriggerResults> TheHLTResults;
191  iEvent.getByLabel( HLTResultsLabel_ , TheHLTResults);
192 
193  bool EventPasses = true;
194 
195  // Make sure handle is valid
196  if( TheHLTResults.isValid() && hltselection_ )
197  {
198  //Get HLT Names
199  const edm::TriggerNames & TheTriggerNames = iEvent.triggerNames(*TheHLTResults);
200 
201  for( unsigned int index = 0 ; index < HLTBitLabel_.size(); index++)
202  {
203  if( HLTBitLabel_[index].label().size() )
204  {
205  //Change the default value since HLT requirement has been issued by the user
206  if( index == 0 ) EventPasses = false;
207  //Get the HLT bit and check to make sure it is valid
208  unsigned int bit = TheTriggerNames.triggerIndex( HLTBitLabel_[index].label().c_str());
209  if( bit < TheHLTResults->size() )
210  {
211  //If any of the HLT names given by the user accept, then the event passes
212  if( TheHLTResults->accept( bit ) && !TheHLTResults->error( bit ) )
213  {
214  EventPasses = true;
215  hCT_NEvents_HLT[index]->Fill(1);
216  }
217  else
218  hCT_NEvents_HLT[index]->Fill(0);
219  }
220  else
221  {
222  edm::LogInfo("OutputInfo")
223  << "The HLT Trigger Name : " << HLTBitLabel_[index].label() << " is not valid for this trigger table " << std::endl;
224  }
225  }
226  }
227  }
228 
229  if( !EventPasses && hltselection_ )
230  return;
231 
232  //----------GREG & CHRIS' idea---///
233  float ETTowerMin = -1; //GeV
234  float METRingMin = -2; // GeV
235 
236  Nevents++;
237  hCT_Nevents->Fill(0);
238 
239  // ==========================================================
240  // Retrieve!
241  // ==========================================================
242 
244  iEvent.getByLabel(caloTowersLabel_, towers);
245 
246  if( (!towers.isValid())) {
247  edm::LogInfo("")<<"CaloTowers "<< caloTowersLabel_<<" not found!"<<std::endl;
248  return;
249  }
250 
251 
252  edm::Handle<bool> HBHENoiseFilterResultHandle;
253  iEvent.getByLabel(HBHENoiseFilterResultLabel_, HBHENoiseFilterResultHandle);
254  bool HBHENoiseFilterResult = *HBHENoiseFilterResultHandle;
255  if (!HBHENoiseFilterResultHandle.isValid()) {
256  LogDebug("") << "CaloTowerAnalyzer: Could not find HBHENoiseFilterResult" << std::endl;
257  }
258 
259 
260  bool bHcalNoiseFilter = HBHENoiseFilterResult;
261 
262  if(!bHcalNoiseFilter) return;
263 
264  edm::View<Candidate>::const_iterator towerCand = towers->begin();
265 
266  // ==========================================================
267  // Fill Histograms!
268  // ==========================================================
269 
270  int CTmin_iphi = 99, CTmax_iphi = -99;
271  int CTmin_ieta = 99, CTmax_ieta = -99;
272 
273  TLorentzVector vMET_EtaRing[83];
274  int ActiveRing[83];
275  int NActiveTowers[83];
276  double SET_EtaRing[83];
277  double MinEt_EtaRing[83];
278  double MaxEt_EtaRing[83];
279  for (int i=0;i<83; i++)
280  {
281  ActiveRing[i] = 0;
282  NActiveTowers[i] = 0;
283  SET_EtaRing[i] = 0;
284  MinEt_EtaRing[i] = 0;
285  MaxEt_EtaRing[i] = 0;
286  }
287 
288  //rcr for (calotower = towerCollection->begin(); calotower != towerCollection->end(); calotower++) {
289 
290  for ( ; towerCand != towers->end(); towerCand++)
291  {
292  const Candidate* candidate = &(*towerCand);
293  if (candidate)
294  {
295  const CaloTower* calotower = dynamic_cast<const CaloTower*> (candidate);
296  if (calotower){
297  //math::RhoEtaPhiVector Momentum = calotower->momentum();
298  double Tower_ET = calotower->et();
299  //double Tower_Energy = calotower->energy();
300  // double Tower_Eta = calotower->eta();
301  double Tower_Phi = calotower->phi();
302  //double Tower_EMEnergy = calotower->emEnergy();
303  //double Tower_HadEnergy = calotower->hadEnergy();
304  double Tower_OuterEt = calotower->outerEt();
305  double Tower_EMEt = calotower->emEt();
306  double Tower_HadEt = calotower->hadEt();
307  //int Tower_EMLV1 = calotower->emLvl1();
308  //int Tower_HadLV1 = calotower->hadLv11();
309  int Tower_ieta = calotower->id().ieta();
310  int Tower_iphi = calotower->id().iphi();
311  int EtaRing = 41+Tower_ieta;
312  ActiveRing[EtaRing] = 1;
313  NActiveTowers[EtaRing]++;
314  SET_EtaRing[EtaRing]+=Tower_ET;
315  TLorentzVector v_;
316  v_.SetPtEtaPhiE(Tower_ET, 0, Tower_Phi, Tower_ET);
317  if (Tower_ET>ETTowerMin)
318  vMET_EtaRing[EtaRing]-=v_;
319 
320  // Fill Histograms
321  hCT_Occ_ieta_iphi->Fill(Tower_ieta,Tower_iphi);
322  if (calotower->emEt() > 0 && calotower->emEt() + calotower->hadEt() > 0.3)
323  hCT_Occ_EM_Et_ieta_iphi->Fill(Tower_ieta,Tower_iphi);
324  if (calotower->hadEt() > 0 && calotower->emEt() + calotower->hadEt() > 0.3)
325  hCT_Occ_HAD_Et_ieta_iphi->Fill(Tower_ieta,Tower_iphi);
326  if (calotower->outerEt() > 0 && calotower->emEt() + calotower->hadEt() > 0.3)
327  hCT_Occ_Outer_Et_ieta_iphi->Fill(Tower_ieta,Tower_iphi);
328 
329  hCT_et_ieta_iphi->Fill(Tower_ieta,Tower_iphi,Tower_ET);
330  hCT_emEt_ieta_iphi->Fill(Tower_ieta,Tower_iphi,Tower_EMEt);
331  hCT_hadEt_ieta_iphi->Fill(Tower_ieta,Tower_iphi,Tower_HadEt);
332  hCT_outerEt_ieta_iphi->Fill(Tower_ieta,Tower_iphi,Tower_OuterEt);
333 
334  if (allhist_){
335  hCT_etvsieta->Fill(Tower_ieta, Tower_ET);
336  hCT_emEtvsieta->Fill(Tower_ieta, Tower_EMEt);
337  hCT_hadEtvsieta->Fill(Tower_ieta,Tower_HadEt);
338  hCT_outerEtvsieta->Fill(Tower_ieta,Tower_OuterEt);
339  }
340 
341  if (Tower_ET > MaxEt_EtaRing[EtaRing])
342  MaxEt_EtaRing[EtaRing] = Tower_ET;
343  if (Tower_ET < MinEt_EtaRing[EtaRing] && Tower_ET>0)
344  MinEt_EtaRing[EtaRing] = Tower_ET;
345 
346 
347  if (Tower_ieta < CTmin_ieta) CTmin_ieta = Tower_ieta;
348  if (Tower_ieta > CTmax_ieta) CTmax_ieta = Tower_ieta;
349  if (Tower_iphi < CTmin_iphi) CTmin_iphi = Tower_iphi;
350  if (Tower_iphi > CTmax_iphi) CTmax_iphi = Tower_iphi;
351  } //end if (calotower) ..
352  } // end if(candidate) ...
353 
354  } // end loop over towers
355 
356  // Fill eta-ring MET quantities
357  if (allhist_){
358  for (int iEtaRing=0; iEtaRing<83; iEtaRing++)
359  {
360  hCT_Minetvsieta->Fill(iEtaRing-41, MinEt_EtaRing[iEtaRing]);
361  hCT_Maxetvsieta->Fill(iEtaRing-41, MaxEt_EtaRing[iEtaRing]);
362 
363  if (ActiveRing[iEtaRing])
364  {
365  if (vMET_EtaRing[iEtaRing].Pt()>METRingMin)
366  {
367  hCT_METPhivsieta->Fill(iEtaRing-41, vMET_EtaRing[iEtaRing].Phi());
368  hCT_MExvsieta->Fill(iEtaRing-41, vMET_EtaRing[iEtaRing].Px());
369  hCT_MEyvsieta->Fill(iEtaRing-41, vMET_EtaRing[iEtaRing].Py());
370  hCT_METvsieta->Fill(iEtaRing-41, vMET_EtaRing[iEtaRing].Pt());
371  }
372  hCT_SETvsieta->Fill(iEtaRing-41, SET_EtaRing[iEtaRing]);
373  hCT_Occvsieta->Fill(iEtaRing-41, NActiveTowers[iEtaRing]);
374  }
375  } // ietaring
376  } // allhist_
377 
378  edm::LogInfo("OutputInfo") << "CT ieta range: " << CTmin_ieta << " " << CTmax_ieta;
379  edm::LogInfo("OutputInfo") << "CT iphi range: " << CTmin_iphi << " " << CTmax_iphi;
380 
381 }
382 
384 {
385 }
#define LogDebug(id)
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
Definition: View.h:81
virtual edm::TriggerNames const & triggerNames(edm::TriggerResults const &triggerResults) const
Definition: Event.cc:204
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
Definition: DQMStore.cc:872
double hadEt() const
Definition: CaloTower.h:83
double outerEt() const
Definition: CaloTower.h:84
bool exists(std::string const &parameterName) const
checks if a parameter exists
virtual float phi() const GCC11_FINAL
momentum azimuthal angle
CaloTowerAnalyzer(const edm::ParameterSet &)
int iEvent
Definition: GenABIO.cc:243
unsigned int triggerIndex(std::string const &name) const
Definition: TriggerNames.cc:32
int iphi() const
get the tower iphi
bool isValid() const
Definition: HandleBase.h:76
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:390
DQMStore * dbe_
CaloTowerDetId id() const
Definition: CaloTower.h:70
virtual void analyze(const edm::Event &, const edm::EventSetup &)
virtual void beginRun(const edm::Run &, const edm::EventSetup &)
double et(double vtxZ) const
Definition: CaloTower.h:99
int ieta() const
get the tower ieta
TH2F * getTH2F(void) const
MonitorElement * book2D(const char *name, const char *title, int nchX, double lowX, double highX, int nchY, double lowY, double highY)
Book 2D histogram.
Definition: DQMStore.cc:1000
virtual void endJob()
double emEt() const
Definition: CaloTower.h:82
tuple size
Write out results.
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:584
Definition: Run.h:41