CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MuonIsolationDQM.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // MuonIsolationDQM.cc
3 // Package: Muon Isolation DQM
4 // Class: MuonIsolationDQM
5 //
6 /*
7 
8 
9  Description: Muon Isolation DQM class
10 
11  Implementation: This code will accept a data set and generate plots of
12  various quantities relevent to the Muon Isolation module. We will
13  be using the IsoDeposit class, *not* the MuonIsolation struct.
14 
15  The sequence of events is...
16  * initalize statics (which variables to plot, axis limtis, etc.)
17  * run over events
18  * for each event, run over the muons in that event
19  *record IsoDeposit data
20  * transfer data to histograms, profile plots
21  * save histograms to a root file
22 
23  Easy-peasy.
24 
25 */
26 //
27 // Original Author: "C. Jess Riedel", UC Santa Barbara
28 // Created: Tue Jul 17 15:58:24 CDT 2007
29 //
30 
31 //Class header file
33 
34 //System included files
35 #include <memory>
36 #include <string>
37 #include <typeinfo>
38 #include <utility>
39 
40 //Root included files
41 #include "TH1.h"
42 #include "TH2.h"
43 #include "TProfile.h"
44 
45 //Event framework included files
49 
50 //Other included files
52 
53 //Using declarations
54 using std::vector;
55 using std::pair;
56 using std::string;
57 
58 
59 
60 //
61 //-----------------Constructors---------------------
62 //
64 {
65 
66  // rootfilename = iConfig.getUntrackedParameter<string>("rootfilename"); // comment out for inclusion
67  requireSTAMuon = iConfig.getUntrackedParameter<bool>("requireSTAMuon");
68  requireTRKMuon = iConfig.getUntrackedParameter<bool>("requireTRKMuon");
69  requireGLBMuon = iConfig.getUntrackedParameter<bool>("requireGLBMuon");
70  dirName = iConfig.getParameter<std::string>("directory");
71  // subDirName = iConfig.getParameter<std::string>("@module_label");
72 
73  // dirName += subDirName;
74 
75  //--------Initialize tags-------
76  Muon_Tag = iConfig.getUntrackedParameter<edm::InputTag>("Global_Muon_Label");
77 
78  //-------Initialize counters----------------
79  nEvents = 0;
80  nSTAMuons = 0;
81  nTRKMuons = 0;
82  nGLBMuons = 0;
83 
84  InitStatics();
85 
86  //Set up DAQ
87  dbe = 0;
89 
90  //------"allocate" space for the data vectors-------
91 
92  /*
93  h_1D is a 2D vector with indices [var][muon#]
94  cd_plots is a 2D vector with indices [var][muon#]
95  */
96  //NOTE:the total number of muons and events is initially unknown,
97  // so that dimension is not initialized. Hence, theMuonData
98  // needs no resizing.
99 
100  h_1D.resize (NUM_VARS);
101  /* cd_plots.resize(NUM_VARS); */
102 
103  dbe->cd();
104 }
105 
106 //
107 //----------Destructor-----------------
108 //
110 
111  //Deallocate memory
112 
113 }
114 
115 //
116 //------------Methods-----------------------
117 //
119 
120  //-----------Initialize primatives-----------
121  S_BIN_WIDTH = 1.0;//in GeV
122  L_BIN_WIDTH = 2.0;//in GeV
124  NUM_LOG_BINS = 15;
125  LOG_BINNING_RATIO = 1.1;//ratio by which each bin is wider than the last for log binning
126  //i.e. bin widths are (x), (r*x), (r^2*x), ..., (r^(nbins)*x)
127 
128 
129  //-------Initialize Titles---------
130  title_sam = "";//"[Sample b-jet events] ";
131  title_cone = "";//" [in R=0.3 IsoDeposit Cone]";
132  //The above two pieces of info will be printed on the title of the whole page,
133  //not for each individual histogram
134  // title_cd = "C.D. of ";
135 
136  //-------"Allocate" memory for vectors
137  main_titles.resize(NUM_VARS);
138  axis_titles.resize(NUM_VARS);
139  names.resize(NUM_VARS);
140  param.resize(NUM_VARS, vector<double>(3) );
141  isContinuous.resize(NUM_VARS);
142 
143  //-----Titles of the plots-----------
144  main_titles[0 ] = "Total Tracker Momentum, #Delta R = 0.3";
145  main_titles[1 ] = "Total EM Cal Energy, #Delta R = 0.3";
146  main_titles[2 ] = "Total Had Cal Energy, #Delta R = 0.3";
147  main_titles[3 ] = "Total HO Cal Energy, #Delta R = 0.3";
148  main_titles[4 ] = "Number of Tracker Tracks, #Delta R = 0.3";
149  main_titles[5 ] = "Number of Jets around Muon, #Delta R = 0.3";
150  main_titles[6 ] = "Tracker p_{T} within veto cone, #Delta R = 0.3";
151  main_titles[7 ] = "EM E_{T} within veto cone, #Delta R = 0.3";
152  main_titles[8 ] = "Had E_{T} within veto cone, #Delta R = 0.3";
153  main_titles[9 ] = "HO E_{T} within veto cone, #Delta R = 0.3";
154  main_titles[10] = "Average Momentum per Track, #Delta R = 0.3";
155  main_titles[11] = "Weighted Energy, #Delta R = 0.3";
156 
157  main_titles[12] = "Total Tracker Momentum, #Delta R = 0.5";
158  main_titles[13] = "Total EM Cal Energy, #Delta R = 0.5";
159  main_titles[14] = "Total Had Cal Energy, #Delta R = 0.5";
160  main_titles[15] = "Total HO Cal Energy, #Delta R = 0.5";
161  main_titles[16] = "Number of Tracker Tracks, #Delta R = 0.5";
162  main_titles[17] = "Number of Jets around Muon, #Delta R = 0.5";
163  main_titles[18] = "Tracker p_{T} within veto cone, #Delta R = 0.5";
164  main_titles[19] = "EM E_{T} within veto cone, #Delta R = 0.5";
165  main_titles[20] = "Had E_{T} within veto cone, #Delta R = 0.5";
166  main_titles[21] = "HO E_{T} within veto cone, #Delta R = 0.5";
167  main_titles[22] = "Average Momentum per Track, #Delta R = 0.5";
168  main_titles[23] = "Weighted Energy, #Delta R = 0.5";
169 
170  //------Titles on the X or Y axis------------
171  axis_titles[0 ] = "#Sigma p_{T} (GeV)";
172  axis_titles[1 ] = "#Sigma E_{T}^{EM} (GeV)";
173  axis_titles[2 ] = "#Sigma E_{T}^{Had} (GeV)";
174  axis_titles[3 ] = "#Sigma E_{T}^{HO} (GeV)";
175  axis_titles[4 ] = "N_{Tracks}";
176  axis_titles[5 ] = "N_{Jets}";
177  axis_titles[6 ] = "#Sigma p_{T,veto} (GeV)";
178  axis_titles[7 ] = "#Sigma E_{T,veto}^{EM} (GeV)";
179  axis_titles[8 ] = "#Sigma E_{T,veto}^{Had} (GeV)";
180  axis_titles[9 ] = "#Sigma E_{T,veto}^{HO} (GeV)";
181  axis_titles[10] = "#Sigma p_{T} / N_{Tracks} (GeV)";
182  axis_titles[11] = "(1.5) X #Sigma E_{T}^{EM} + #Sigma E_{T}^{Had}";
183 
184  axis_titles[12] = "#Sigma p_{T} (GeV)";
185  axis_titles[13] = "#Sigma E_{T}^{EM} (GeV)";
186  axis_titles[14] = "#Sigma E_{T}^{Had} (GeV)";
187  axis_titles[15] = "#Sigma E_{T}^{HO} (GeV)";
188  axis_titles[16] = "N_{Tracks}";
189  axis_titles[17] = "N_{Jets}";
190  axis_titles[18] = "#Sigma p_{T,veto} (GeV)";
191  axis_titles[19] = "#Sigma E_{T,veto}^{EM} (GeV)";
192  axis_titles[20] = "#Sigma E_{T,veto}^{Had} (GeV)";
193  axis_titles[21] = "#Sigma E_{T,veto}^{HO} (GeV)";
194  axis_titles[22] = "#Sigma p_{T} / N_{Tracks} (GeV)";
195  axis_titles[23] = "(1.5) X #Sigma E_{T}^{EM} + #Sigma E_{T}^{Had}";
196 
197  //-----------Names given for the root file----------
198  names[0 ] = "sumPt_R03";
199  names[1 ] = "emEt_R03";
200  names[2 ] = "hadEt_R03";
201  names[3 ] = "hoEt_R03";
202  names[4 ] = "nTracks_R03";
203  names[5 ] = "nJets_R03";
204  names[6 ] = "trackerVetoPt_R03";
205  names[7 ] = "emVetoEt_R03";
206  names[8 ] = "hadVetoEt_R03";
207  names[9 ] = "hoVetoEt_R03";
208  names[10] = "avgPt_R03";
209  names[11] = "weightedEt_R03";
210 
211  names[12] = "sumPt_R05";
212  names[13] = "emEt_R05";
213  names[14] = "hadEt_R05";
214  names[15] = "hoEt_R05";
215  names[16] = "nTracks_R05";
216  names[17] = "nJets_R05";
217  names[18] = "trackerVetoPt_R05";
218  names[19] = "emVetoEt_R05";
219  names[20] = "hadVetoEt_R05";
220  names[21] = "hoVetoEt_R05";
221  names[22] = "avgPt_R05";
222  names[23] = "weightedEt_R05";
223 
224  //----------Parameters for binning of histograms---------
225  //param[var][0] is the number of bins
226  //param[var][1] is the low edge of the low bin
227  //param[var][2] is the high edge of the high bin
228  //
229  // maximum value------,
230  // |
231  // V
232  param[0 ][0]= (int)( 20.0/S_BIN_WIDTH); param[0 ][1]= 0.0; param[0 ][2]= param[0 ][0]*S_BIN_WIDTH;
233  param[1 ][0]= (int)( 20.0/S_BIN_WIDTH); param[1 ][1]= 0.0; param[1 ][2]= param[1 ][0]*S_BIN_WIDTH;
234  param[2 ][0]= (int)( 20.0/S_BIN_WIDTH); param[2 ][1]= 0.0; param[2 ][2]= param[2 ][0]*S_BIN_WIDTH;
235  param[3 ][0]= 20; param[3 ][1]= 0.0; param[3 ][2]= 2.0;
236  param[4 ][0]= 16; param[4 ][1]= -0.5; param[4 ][2]= param[4 ][0]-0.5;
237  param[5 ][0]= 4; param[5 ][1]= -0.5; param[5 ][2]= param[5 ][0]-0.5;
238  param[6 ][0]= (int)( 40.0/S_BIN_WIDTH); param[6 ][1]= 0.0; param[6 ][2]= param[6 ][0]*S_BIN_WIDTH;
239  param[7 ][0]= 20; param[7 ][1]= 0.0; param[7 ][2]= 10.0;
240  param[8 ][0]= (int)( 20.0/S_BIN_WIDTH); param[8 ][1]= 0.0; param[8 ][2]= param[8 ][0]*S_BIN_WIDTH;
241  param[9 ][0]= 20; param[9 ][1]= 0.0; param[9 ][2]= 5.0;
242  param[10][0]= (int)( 15.0/S_BIN_WIDTH); param[10][1]= 0.0; param[10][2]= param[10][0]*S_BIN_WIDTH;
243  param[11][0]= (int)( 20.0/S_BIN_WIDTH); param[11][1]= 0.0; param[11][2]= param[11][0]*S_BIN_WIDTH;
244 
245  param[12][0]= (int)( 20.0/S_BIN_WIDTH); param[12][1]= 0.0; param[12][2]= param[12][0]*S_BIN_WIDTH;
246  param[13][0]= (int)( 20.0/S_BIN_WIDTH); param[13][1]= 0.0; param[13][2]= param[13][0]*S_BIN_WIDTH;
247  param[14][0]= (int)( 20.0/S_BIN_WIDTH); param[14][1]= 0.0; param[14][2]= param[14][0]*S_BIN_WIDTH;
248  param[15][0]= 20; param[15][1]= 0.0; param[15][2]= 2.0;
249  param[16][0]= 16; param[16][1]= -0.5; param[16][2]= param[16][0]-0.5;
250  param[17][0]= 4; param[17][1]= -0.5; param[17][2]= param[17][0]-0.5;
251  param[18][0]= (int)( 40.0/S_BIN_WIDTH); param[18][1]= 0.0; param[18][2]= param[18][0]*S_BIN_WIDTH;
252  param[19][0]= 20; param[19][1]= 0.0; param[19][2]= 10.0;
253  param[20][0]= (int)( 20.0/S_BIN_WIDTH); param[20][1]= 0.0; param[20][2]= param[20][0]*S_BIN_WIDTH;
254  param[21][0]= 20; param[21][1]= 0.0; param[21][2]= 5.0;
255  param[22][0]= (int)( 15.0/S_BIN_WIDTH); param[22][1]= 0.0; param[22][2]= param[22][0]*S_BIN_WIDTH;
256  param[23][0]= (int)( 20.0/S_BIN_WIDTH); param[23][1]= 0.0; param[23][2]= param[23][0]*S_BIN_WIDTH;
257 
258  //--------------Is the variable continuous (i.e. non-integer)?-------------
259  //---------(Log binning will only be used for continuous variables)--------
260  isContinuous[0 ] = 1;
261  isContinuous[1 ] = 1;
262  isContinuous[2 ] = 1;
263  isContinuous[3 ] = 1;
264  isContinuous[4 ] = 0;
265  isContinuous[5 ] = 0;
266  isContinuous[6 ] = 1;
267  isContinuous[7 ] = 1;
268  isContinuous[8 ] = 1;
269  isContinuous[9 ] = 1;
270  isContinuous[10] = 1;
271  isContinuous[11] = 1;
272 
273  isContinuous[12] = 1;
274  isContinuous[13] = 1;
275  isContinuous[14] = 1;
276  isContinuous[15] = 1;
277  isContinuous[16] = 0;
278  isContinuous[17] = 0;
279  isContinuous[18] = 1;
280  isContinuous[19] = 1;
281  isContinuous[20] = 1;
282  isContinuous[21] = 1;
283  isContinuous[22] = 1;
284  isContinuous[23] = 1;
285 
286 }
287 
288 
289 // ------------ method called for each event ------------
291 
292  ++nEvents;
293  edm::LogInfo("Tutorial") << "\nInvestigating event #" << nEvents<<"\n";
294 
295  // Get Muon Collection
296  edm::Handle<edm::View<reco::Muon> > muonsHandle; //
297  iEvent.getByLabel(Muon_Tag, muonsHandle);
298 
299  //Fill event entry in histogram of number of muons
300  edm::LogInfo("Tutorial") << "Number of Muons: " << muonsHandle->size();
301  theMuonData = muonsHandle->size();
303 
304  //Fill historgams concerning muon isolation
305  uint iMuon=0;
306  dbe->setCurrentFolder(dirName.c_str());
307  for (MuonIterator muon = muonsHandle->begin(); muon != muonsHandle->end(); ++muon, ++iMuon ) {
308  // ++nMuons;
309  if (requireSTAMuon && muon->isStandAloneMuon()) {
310  ++nSTAMuons;
311  RecordData(muon);
312  FillHistos();
313  }
314  else if (requireTRKMuon && muon->isTrackerMuon()) {
315  ++nTRKMuons;
316  RecordData(muon);
317  FillHistos();
318  }
319  else if (requireGLBMuon && muon->isGlobalMuon()) {
320  ++nGLBMuons;
321  RecordData(muon);
322  FillHistos();
323  }
324  }
325  dbe->cd();
326 
327 }
328 
329 //---------------Record data for a signle muon's data---------------------
331 
332 
333  theData[0] = muon->isolationR03().sumPt;
334  theData[1] = muon->isolationR03().emEt;
335  theData[2] = muon->isolationR03().hadEt;
336  theData[3] = muon->isolationR03().hoEt;
337 
338  theData[4] = muon->isolationR03().nTracks;
339  theData[5] = muon->isolationR03().nJets;
340  theData[6] = muon->isolationR03().trackerVetoPt;
341  theData[7] = muon->isolationR03().emVetoEt;
342  theData[8] = muon->isolationR03().hadVetoEt;
343  theData[9] = muon->isolationR03().hoVetoEt;
344 
345  // make sure nTracks != 0 before filling this one
346  if (theData[4] != 0) theData[10] = (double)theData[0] / (double)theData[4];
347  else theData[10] = -99;
348 
349  theData[11] = 1.5 * theData[1] + theData[2];
350 
351  theData[12] = muon->isolationR05().sumPt;
352  theData[13] = muon->isolationR05().emEt;
353  theData[14] = muon->isolationR05().hadEt;
354  theData[15] = muon->isolationR05().hoEt;
355 
356  theData[16] = muon->isolationR05().nTracks;
357  theData[17] = muon->isolationR05().nJets;
358  theData[18] = muon->isolationR05().trackerVetoPt;
359  theData[19] = muon->isolationR05().emVetoEt;
360  theData[20] = muon->isolationR05().hadVetoEt;
361  theData[21] = muon->isolationR05().hoVetoEt;
362 
363  // make sure nTracks != 0 before filling this one
364  if (theData[16] != 0) theData[22] = (double)theData[12] / (double)theData[16];
365  else theData[22] = -99;
366 
367  theData[23] = 1.5 * theData[13] + theData[14];
368 
369 }
370 
371 // ------------ method called once each job just before starting event loop ------------
372 void
374 {
375 
376  edm::LogInfo("Tutorial") << "\n#########################################\n\n"
377  << "Lets get started! "
378  << "\n\n#########################################\n";
379  dbe->setCurrentFolder(dirName.c_str());
380  InitHistos();
381  dbe->cd();
382 
383 }
384 
385 // ------------ method called once each job just after ending the event loop ------------
386 void
388 
389  // check if ME still there (and not killed by MEtoEDM for memory saving)
390  if( dbe )
391  {
392  // check existence of first histo in the list
393  if (! dbe->get(dirName+"/nMuons")) return;
394  }
395  else
396  return;
397 
398  edm::LogInfo("Tutorial") << "\n#########################################\n\n"
399  << "Total Number of Events: " << nEvents
400  << "\n\n#########################################\n"
401  << "\nInitializing Histograms...\n";
402 
403  edm::LogInfo("Tutorial") << "\nIntializing Finished. Filling...\n";
404  NormalizeHistos();
405  edm::LogInfo("Tutorial") << "\nFilled. Saving...\n";
406  // dbe->save(rootfilename); // comment out for incorporation
407  edm::LogInfo("Tutorial") << "\nSaved. Peace, homie, I'm out.\n";
408 
409 }
410 
412 
413  //---initialize number of muons histogram---
414  h_nMuons = dbe->book1D("nMuons", title_sam + "Number of Muons", 20, 0., 20.);
415  h_nMuons->setAxisTitle("Number of Muons",XAXIS);
416  h_nMuons->setAxisTitle("Fraction of Events",YAXIS);
417 
418 
419  //---Initialize 1D Histograms---
420  for(int var = 0; var < NUM_VARS; var++){
421  h_1D[var] = dbe->book1D(
422  names[var],
423  title_sam + main_titles[var] + title_cone,
424  (int)param[var][0],
425  param[var][1],
426  param[var][2]
427  );
428  /* cd_plots[var] = dbe->book1D(
429  names[var] + "_cd",
430  title_sam + title_cd + main_titles[var] + title_cone,
431  (int)param[var][0],
432  param[var][1],
433  param[var][2]
434  );
435  */
436 
437  h_1D[var]->setAxisTitle(axis_titles[var],XAXIS);
438  // h_1D[var]->setAxisTitle("Fraction of Muons",YAXIS);
439  GetTH1FromMonitorElement(h_1D[var])->Sumw2();
440 
441  /* cd_plots[var]->setAxisTitle(axis_titles[var],XAXIS);
442  cd_plots[var]->setAxisTitle("Fraction of Muons",YAXIS);
443  GetTH1FromMonitorElement(cd_plots[var])->Sumw2();
444  */
445 
446  }//Finish 1D
447 
448  //avg pT not defined for zero tracks.
449  //MonitorElement is inflxible and won't let me change the
450  //number of bins! I guess all I'm doing here is changing
451  //range of the x axis when it is printed, not the actual
452  //bins that are filled
453  // p_2D[4][9]->setAxisRange(0.5,15.5,XAXIS);
454 
455 }
456 
458  for(int var=0; var<NUM_VARS; var++){
459  //turn cd_plots into CDF's
460  //underflow -> bin #0. overflow -> bin #(nbins+1)
461  //0th bin doesn't need changed
462 
463  double entries = GetTH1FromMonitorElement(h_1D[var])->GetEntries();
464 
465  /* int n_max = int(param[var][0])+1;
466  for(int n=1; n<=n_max; ++n){
467  cd_plots[var]->setBinContent(n, cd_plots[var]->getBinContent(n) + cd_plots[var]->getBinContent(n-1)); //Integrate.
468  }
469  */
470  //----normalize------
471  /* if (requireCombinedMuon) {
472  GetTH1FromMonitorElement(h_1D[var])->Scale(1./entries);
473  GetTH1FromMonitorElement(cd_plots[var])->Scale(1./entries);
474  }
475  else {
476  */
477  GetTH1FromMonitorElement(h_1D[var])->Scale(1./entries);
478  // GetTH1FromMonitorElement(cd_plots[var])->Scale(1./entries);
479  }
480 }
481 
483 
484  int overFlowBin;
485  double overFlow = 0;
486 
487  //----------Fill 1D histograms---------------
488  for(int var=0; var<NUM_VARS; var++){
489  h_1D[var]->Fill(theData[var]);
490  // cd_plots[var]->Fill(theData[var]);//right now, this is a regular PDF (just like h_1D)
491  if (theData[var] > param[var][2]) {
492  // fill the overflow bin
493  overFlowBin = (int) param[var][0] + 1;
494  overFlow = GetTH1FromMonitorElement(h_1D[var])->GetBinContent(overFlowBin);
495  GetTH1FromMonitorElement(h_1D[var])->SetBinContent(overFlowBin, overFlow + 1);
496  }
497  }//Finish 1D
498 
499 }
500 
502  return me->getTH1();
503 }
504 
505 //define this as a plug-in
void RecordData(MuonIterator muon)
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
std::string dirName
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
Definition: DQMStore.cc:514
MonitorElement * h_nMuons
void cd(void)
go to top directory (ie. root)
Definition: DQMStore.cc:209
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
MuonIsolationDQM(const edm::ParameterSet &)
std::vector< MonitorElement * > h_1D
std::vector< std::string > main_titles
edm::View< reco::Muon >::const_iterator MuonIterator
std::vector< std::string > names
void Fill(long long x)
std::vector< std::vector< double > > param
int iEvent
Definition: GenABIO.cc:243
virtual void beginJob(void)
std::string title_sam
TH1 * getTH1(void) const
std::string title_cone
std::vector< int > isContinuous
MonitorElement * get(const std::string &path) const
get ME from full pathname (e.g. &quot;my/long/dir/my_histo&quot;)
Definition: DQMStore.cc:1265
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:355
virtual void analyze(const edm::Event &, const edm::EventSetup &)
std::vector< std::string > axis_titles
edm::InputTag Muon_Tag
static const int NUM_VARS
virtual void endJob()
double theData[NUM_VARS]
void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:232
TH1 * GetTH1FromMonitorElement(MonitorElement *me)