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