CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MuIsoValidation.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // MuIsoValidation.cc
3 // Package: Muon Isolation Validation
4 // Class: MuIsoValidation
5 //
6 /*
7 
8 
9  Description: Muon Isolation Validation 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  requireCombinedMuon = iConfig.getUntrackedParameter<bool>("requireCombinedMuon");
68  dirName = iConfig.getParameter<std::string>("directory");
69  // subDirName = iConfig.getParameter<std::string>("@module_label");
70 
71  // dirName += subDirName;
72 
73  //--------Initialize tags-------
74  Muon_Tag = iConfig.getUntrackedParameter<edm::InputTag>("Global_Muon_Label");
75 
76  //-------Initialize counters----------------
77  nEvents = 0;
78  nIncMuons = 0;
79  nCombinedMuons = 0;
80 
81  InitStatics();
82 
83  //Set up DAQ
84  dbe = 0;
86 
87  //------"allocate" space for the data vectors-------
88 
89  /*
90  h_1D is a 2D vector with indices [var][muon#]
91  cd_plots is a 2D vector with indices [var][muon#]
92  h_2D is a 3D vector with indices [var][var][muon#]
93  p_2D is a 3D vector with indices [var][var][muon#]
94  */
95  //NOTE:the total number of muons and events is initially unknown,
96  // so that dimension is not initialized. Hence, theMuonData
97  // needs no resizing.
98 
99  h_1D.resize (NUM_VARS);
100  cd_plots.resize(NUM_VARS);
101  // h_2D.resize(NUM_VARS, vector<MonitorElement*> (NUM_VARS));
102  p_2D.resize(NUM_VARS, vector<MonitorElement*>(NUM_VARS));
103 
104  dbe->cd();
105 }
106 
107 //
108 //----------Destructor-----------------
109 //
111 
112  //Deallocate memory
113 
114 }
115 
116 //
117 //------------Methods-----------------------
118 //
120 
121  //-----------Initialize primatives-----------
122  S_BIN_WIDTH = 1.0;//in GeV
123  L_BIN_WIDTH = 2.0;//in GeV
125  NUM_LOG_BINS = 15;
126  LOG_BINNING_RATIO = 1.1;//ratio by which each bin is wider than the last for log binning
127  //i.e. bin widths are (x), (r*x), (r^2*x), ..., (r^(nbins)*x)
128 
129 
130  //-------Initialize Titles---------
131  title_sam = "";//"[Sample b-jet events] ";
132  title_cone = "";//" [in R=0.3 IsoDeposit Cone]";
133  //The above two pieces of info will be printed on the title of the whole page,
134  //not for each individual histogram
135  title_cd = "C.D. of ";
136 
137  //-------"Allocate" memory for vectors
138  main_titles.resize(NUM_VARS);
139  axis_titles.resize(NUM_VARS);
140  names.resize(NUM_VARS);
141  param.resize(NUM_VARS, vector<double>(3) );
142  isContinuous.resize(NUM_VARS);
143 
144  //-----Titles of the plots-----------
145  main_titles[0 ] = "Total Tracker Momentum";
146  main_titles[1 ] = "Total EM Cal Energy";
147  main_titles[2 ] = "Total Had Cal Energy";
148  main_titles[3 ] = "Total HO Cal Energy";
149  main_titles[4 ] = "Number of Tracker Tracks";
150  main_titles[5 ] = "Number of Jets around Muon";
151  main_titles[6 ] = "Tracker p_{T} within veto cone";
152  main_titles[7 ] = "EM E_{T} within veto cone";
153  main_titles[8 ] = "Had E_{T} within veto cone";
154  main_titles[9 ] = "HO E_{T} within veto cone";
155  main_titles[10] = "Muon p_{T}";
156  main_titles[11] = "Muon #eta";
157  main_titles[12] = "Muon #phi";
158  main_titles[13] = "Average Momentum per Track ";
159  main_titles[14] = "Weighted Energy";
160 
161  //------Titles on the X or Y axis------------
162  axis_titles[0 ] = "#Sigma p_{T} (GeV)";
163  axis_titles[1 ] = "#Sigma E_{T}^{EM} (GeV)";
164  axis_titles[2 ] = "#Sigma E_{T}^{Had} (GeV)";
165  axis_titles[3 ] = "#Sigma E_{T}^{HO} (GeV)";
166  axis_titles[4 ] = "N_{Tracks}";
167  axis_titles[5 ] = "N_{Jets}";
168  axis_titles[6 ] = "#Sigma p_{T,veto} (GeV)";
169  axis_titles[7 ] = "#Sigma E_{T,veto}^{EM} (GeV)";
170  axis_titles[8 ] = "#Sigma E_{T,veto}^{Had} (GeV)";
171  axis_titles[9 ] = "#Sigma E_{T,veto}^{HO} (GeV)";
172  axis_titles[10] = "p_{T,#mu} (GeV)";
173  axis_titles[11] = "#eta_{#mu}";
174  axis_titles[12] = "#phi_{#mu}";
175  axis_titles[13] = "#Sigma p_{T} / N_{Tracks} (GeV)";
176  axis_titles[14] = "(1.5) X #Sigma E_{T}^{EM} + #Sigma E_{T}^{Had}";
177 
178  //-----------Names given for the root file----------
179  names[0 ] = "sumPt";
180  names[1 ] = "emEt";
181  names[2 ] = "hadEt";
182  names[3 ] = "hoEt";
183  names[4 ] = "nTracks";
184  names[5 ] = "nJets";
185  names[6 ] = "trackerVetoPt";
186  names[7 ] = "emVetoEt";
187  names[8 ] = "hadVetoEt";
188  names[9 ] = "hoVetoEt";
189  names[10] = "muonPt";
190  names[11] = "muonEta";
191  names[12] = "muonPhi";
192  names[13] = "avgPt";
193  names[14] = "weightedEt";
194 
195  //----------Parameters for binning of histograms---------
196  //param[var][0] is the number of bins
197  //param[var][1] is the low edge of the low bin
198  //param[var][2] is the high edge of the high bin
199  //
200  // maximum value------,
201  // |
202  // V
203  param[0 ][0]= (int)( 20.0/S_BIN_WIDTH); param[0 ][1]= 0.0; param[0 ][2]= param[0 ][0]*S_BIN_WIDTH;
204  param[1 ][0]= (int)( 20.0/S_BIN_WIDTH); param[1 ][1]= 0.0; param[1 ][2]= param[1 ][0]*S_BIN_WIDTH;
205  param[2 ][0]= (int)( 20.0/S_BIN_WIDTH); param[2 ][1]= 0.0; param[2 ][2]= param[2 ][0]*S_BIN_WIDTH;
206  param[3 ][0]= 20; param[3 ][1]= 0.0; param[3 ][2]= 2.0;
207  param[4 ][0]= 16; param[4 ][1]= -0.5; param[4 ][2]= param[4 ][0]-0.5;
208  param[5 ][0]= 4; param[5 ][1]= -0.5; param[5 ][2]= param[5 ][0]-0.5;
209  param[6 ][0]= (int)( 40.0/S_BIN_WIDTH); param[6 ][1]= 0.0; param[6 ][2]= param[6 ][0]*S_BIN_WIDTH;
210  param[7 ][0]= 20; param[7 ][1]= 0.0; param[7 ][2]= 10.0;
211  param[8 ][0]= (int)( 20.0/S_BIN_WIDTH); param[8 ][1]= 0.0; param[8 ][2]= param[8 ][0]*S_BIN_WIDTH;
212  param[9 ][0]= 20; param[9 ][1]= 0.0; param[9 ][2]= 5.0;
213  param[10][0]= (int)( 40.0/S_BIN_WIDTH); param[10][1]= 0.0; param[10][2]= param[10][0]*S_BIN_WIDTH;
214  param[11][0]= 24; param[11][1]= -2.4; param[11][2]= 2.4;
215  param[12][0]= 32; param[12][1]= -3.2; param[12][2]= 3.2;
216  param[13][0]= (int)( 15.0/S_BIN_WIDTH); param[13][1]= 0.0; param[13][2]= param[13][0]*S_BIN_WIDTH;
217  param[14][0]= (int)( 20.0/S_BIN_WIDTH); param[14][1]= 0.0; param[14][2]= param[14][0]*S_BIN_WIDTH;
218 
219  //--------------Is the variable continuous (i.e. non-integer)?-------------
220  //---------(Log binning will only be used for continuous variables)--------
221  isContinuous[0 ] = 1;
222  isContinuous[1 ] = 1;
223  isContinuous[2 ] = 1;
224  isContinuous[3 ] = 1;
225  isContinuous[4 ] = 0;
226  isContinuous[5 ] = 0;
227  isContinuous[6 ] = 1;
228  isContinuous[7 ] = 1;
229  isContinuous[8 ] = 1;
230  isContinuous[9 ] = 1;
231  isContinuous[10] = 1;
232  isContinuous[11] = 1;
233  isContinuous[12] = 1;
234  isContinuous[13] = 1;
235  isContinuous[14] = 1;
236 
237 }
238 
239 
240 // ------------ method called for each event ------------
242 
243  ++nEvents;
244  edm::LogInfo("Tutorial") << "\nInvestigating event #" << nEvents<<"\n";
245 
246  // Get Muon Collection
247  edm::Handle<edm::View<reco::Muon> > muonsHandle; //
248  iEvent.getByLabel(Muon_Tag, muonsHandle);
249 
250  //Fill event entry in histogram of number of muons
251  edm::LogInfo("Tutorial") << "Number of Muons: " << muonsHandle->size();
252  theMuonData = muonsHandle->size();
254 
255  //Fill historgams concerning muon isolation
256  uint iMuon=0;
257  dbe->setCurrentFolder(dirName.c_str());
258  for (MuonIterator muon = muonsHandle->begin(); muon != muonsHandle->end(); ++muon, ++iMuon ) {
259  ++nIncMuons;
260  if (requireCombinedMuon) {
261  if (muon->combinedMuon().isNull()) continue;
262  }
263  ++nCombinedMuons;
264  RecordData(muon);
265  FillHistos();
266  }
267  dbe->cd();
268 
269 }
270 
271 //---------------Record data for a signle muon's data---------------------
273 
274 
275  theData[0] = muon->isolationR03().sumPt;
276  theData[1] = muon->isolationR03().emEt;
277  theData[2] = muon->isolationR03().hadEt;
278  theData[3] = muon->isolationR03().hoEt;
279 
280  theData[4] = muon->isolationR03().nTracks;
281  theData[5] = muon->isolationR03().nJets;
282  theData[6] = muon->isolationR03().trackerVetoPt;
283  theData[7] = muon->isolationR03().emVetoEt;
284  theData[8] = muon->isolationR03().hadVetoEt;
285  theData[9] = muon->isolationR03().hoVetoEt;
286 
287  theData[10] = muon->pt();
288  theData[11] = muon->eta();
289  theData[12] = muon->phi();
290 
291  // make sure nTracks != 0 before filling this one
292  if (theData[4] != 0) theData[13] = (double)theData[0] / (double)theData[4];
293  else theData[13] = -99;
294 
295  theData[14] = 1.5 * theData[1] + theData[2];
296 
297 }
298 
299 // ------------ method called once each job just before starting event loop ------------
300 void
302 {
303 
304  edm::LogInfo("Tutorial") << "\n#########################################\n\n"
305  << "Lets get started! "
306  << "\n\n#########################################\n";
307  dbe->setCurrentFolder(dirName.c_str());
308  InitHistos();
309  dbe->cd();
310 
311 }
312 
313 // ------------ method called once each job just after ending the event loop ------------
314 void
316 
317  // check if ME still there (and not killed by MEtoEDM for memory saving)
318  if( dbe )
319  {
320  // check existence of first histo in the list
321  if (! dbe->get(dirName+"/nMuons")) return;
322  }
323  else
324  return;
325 
326  edm::LogInfo("Tutorial") << "\n#########################################\n\n"
327  << "Total Number of Events: " << nEvents
328  << "\nTotal Number of Muons: " << nIncMuons
329  << "\n\n#########################################\n"
330  << "\nInitializing Histograms...\n";
331 
332  edm::LogInfo("Tutorial") << "\nIntializing Finished. Filling...\n";
333  NormalizeHistos();
334  edm::LogInfo("Tutorial") << "\nFilled. Saving...\n";
335  // dbe->save(rootfilename); // comment out for incorporation
336  edm::LogInfo("Tutorial") << "\nSaved. Peace, homie, I'm out.\n";
337 
338 }
339 
341 
342  //---initialize number of muons histogram---
343  h_nMuons = dbe->book1D("nMuons", title_sam + "Number of Muons", 20, 0., 20.);
344  h_nMuons->setAxisTitle("Number of Muons",XAXIS);
345  h_nMuons->setAxisTitle("Fraction of Events",YAXIS);
346 
347 
348  //---Initialize 1D Histograms---
349  for(int var = 0; var < NUM_VARS; var++){
350  h_1D[var] = dbe->book1D(
351  names[var],
352  title_sam + main_titles[var] + title_cone,
353  (int)param[var][0],
354  param[var][1],
355  param[var][2]
356  );
357  cd_plots[var] = dbe->book1D(
358  names[var] + "_cd",
360  (int)param[var][0],
361  param[var][1],
362  param[var][2]
363  );
364 
365  h_1D[var]->setAxisTitle(axis_titles[var],XAXIS);
366  h_1D[var]->setAxisTitle("Fraction of Muons",YAXIS);
367  GetTH1FromMonitorElement(h_1D[var])->Sumw2();
368 
369  cd_plots[var]->setAxisTitle(axis_titles[var],XAXIS);
370  cd_plots[var]->setAxisTitle("Fraction of Muons",YAXIS);
371  GetTH1FromMonitorElement(cd_plots[var])->Sumw2();
372 
373  }//Finish 1D
374 
375  //---Initialize 2D Histograms---
376  for(int var1 = 0; var1 < NUM_VARS; var1++){
377  for(int var2 = 0; var2 < NUM_VARS; var2++){
378  if(var1 == var2) continue;
379 
380  /* h_2D[var1][var2] = dbe->book2D(
381  names[var1] + "_" + names[var2] + "_s",
382  //title is in "y-var vs. x-var" format
383  title_sam + main_titles[var2] + " <vs> " + main_titles[var1] + title_cone,
384  (int)param[var1][0],
385  param[var1][1],
386  param[var1][2],
387  (int)param[var2][0],
388  param[var2][1],
389  param[var2][2]
390  );
391  */
392  //Monitor elements is weird and takes y axis parameters as well
393  //as x axis parameters for a 1D profile plot
394  p_2D[var1][var2] = dbe->bookProfile(
395  names[var1] + "_" + names[var2],
396  title_sam + main_titles[var2] + " <vs> " + main_titles[var1] + title_cone,
397  (int)param[var1][0],
398  param[var1][1],
399  param[var1][2],
400  (int)param[var2][0], //documentation says this is disregarded
401  param[var2][1], //does this do anything?
402  param[var2][2], //does this do anything?
403  " " //profile errors = spread/sqrt(num_datums)
404  );
405 
406  if(LOG_BINNING_ENABLED && isContinuous[var1]){
407  Double_t * bin_edges = new Double_t[NUM_LOG_BINS+1];
408  // nbins+1 because there is one more edge than there are bins
409  MakeLogBinsForProfile(bin_edges, param[var1][1], param[var1][2]);
410  GetTProfileFromMonitorElement(p_2D[var1][var2])->
411  SetBins(NUM_LOG_BINS, bin_edges);
412  delete[] bin_edges;
413  }
414  /* h_2D[var1][var2]->setAxisTitle(axis_titles[var1],XAXIS);
415  h_2D[var1][var2]->setAxisTitle(axis_titles[var2],YAXIS);
416  GetTH2FromMonitorElement(h_2D[var1][var2])->Sumw2();
417  */
418 
419  p_2D[var1][var2]->setAxisTitle(axis_titles[var1],XAXIS);
420  p_2D[var1][var2]->setAxisTitle(axis_titles[var2],YAXIS);
421  // GetTProfileFromMonitorElement(p_2D[var1][var2])->Sumw2();
422  }
423  }//Finish 2D
424 
425 
426 
427  //avg pT not defined for zero tracks.
428  //MonitorElement is inflxible and won't let me change the
429  //number of bins! I guess all I'm doing here is changing
430  //range of the x axis when it is printed, not the actual
431  //bins that are filled
432  p_2D[4][9]->setAxisRange(0.5,15.5,XAXIS);
433 
434 }
435 
436 void MuIsoValidation::MakeLogBinsForProfile(Double_t* bin_edges, const double min,
437  const double max){
438 
439  const double &r = LOG_BINNING_RATIO;
440  const int &nbins = NUM_LOG_BINS;
441 
442  const double first_bin_width = (r > 1.0) ? //so we don't divide by zero
443  (max - min)*(1-r)/(1-pow(r,nbins)) :
444  (max - min)/nbins;
445 
446  bin_edges[0] = min;
447  bin_edges[1] = min + first_bin_width;
448  for(int n = 2; n<nbins; ++n){
449  bin_edges[n] = bin_edges[n-1] + (bin_edges[n-1] - bin_edges[n-2])*r;
450  }
451  bin_edges[nbins] = max;
452 }
453 
455  for(int var=0; var<NUM_VARS; var++){
456  //turn cd_plots into CDF's
457  //underflow -> bin #0. overflow -> bin #(nbins+1)
458  //0th bin doesn't need changed
459 
460  double entries = GetTH1FromMonitorElement(h_1D[var])->GetEntries();
461  if (entries==0)continue;
462  int n_max = int(param[var][0])+1;
463  for(int n=1; n<=n_max; ++n){
464  cd_plots[var]->setBinContent(n, cd_plots[var]->getBinContent(n) + cd_plots[var]->getBinContent(n-1)); //Integrate.
465  }
466  //----normalize------
467  if (requireCombinedMuon) {
468  GetTH1FromMonitorElement(h_1D[var])->Scale(1./entries);
469  GetTH1FromMonitorElement(cd_plots[var])->Scale(1./entries);
470  }
471  else {
472  GetTH1FromMonitorElement(h_1D[var])->Scale(1./entries);
473  GetTH1FromMonitorElement(cd_plots[var])->Scale(1./entries);
474  }
475  }
476 }
477 
479 
480  int overFlowBin;
481  double overFlow = 0;
482 
483  //----------Fill 1D histograms---------------
484  for(int var=0; var<NUM_VARS; var++){
485  h_1D[var]->Fill(theData[var]);
486  cd_plots[var]->Fill(theData[var]);//right now, this is a regular PDF (just like h_1D)
487  if (theData[var] > param[var][2]) {
488  // fill the overflow bin
489  overFlowBin = (int) param[var][0] + 1;
490  overFlow = GetTH1FromMonitorElement(h_1D[var])->GetBinContent(overFlowBin);
491  GetTH1FromMonitorElement(h_1D[var])->SetBinContent(overFlowBin, overFlow + 1);
492  }
493  }//Finish 1D
494 
495  //----------Fill 2D histograms---------------
496  for(int var1=0; var1<NUM_VARS; ++var1){
497  for(int var2=0; var2<NUM_VARS; ++var2){
498  if(var1 == var2) continue;
499  //change below to regular int interating!
500  // h_2D[var1][var2]->Fill(theData[var1], theData[var2]);
501  p_2D[var1][var2]->Fill(theData[var1], theData[var2]);
502  }
503  }//Finish 2D
504 }
505 
507  return me->getTH1();
508 }
509 
511  return me->getTH2F();
512 }
513 
515  return me->getTProfile();
516 }
517 
518 
519 //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::vector< int > isContinuous
TH1 * GetTH1FromMonitorElement(MonitorElement *me)
std::vector< std::vector< MonitorElement * > > p_2D
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
Definition: DQMStore.cc:519
std::vector< std::string > names
std::string title_cd
void cd(void)
go to top directory (ie. root)
Definition: DQMStore.cc:214
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
edm::View< reco::Muon >::const_iterator MuonIterator
std::vector< std::string > main_titles
TH2 * GetTH2FromMonitorElement(MonitorElement *me)
double theData[NUM_VARS]
virtual void beginJob()
MuIsoValidation(const edm::ParameterSet &)
#define min(a, b)
Definition: mlp_lapack.h:161
std::vector< MonitorElement * > h_1D
TProfile * GetTProfileFromMonitorElement(MonitorElement *me)
void Fill(long long x)
MonitorElement * h_nMuons
virtual void analyze(const edm::Event &, const edm::EventSetup &)
int iEvent
Definition: GenABIO.cc:243
std::string title_sam
const T & max(const T &a, const T &b)
std::vector< MonitorElement * > cd_plots
std::string dirName
std::string title_cone
TH1 * getTH1(void) const
MonitorElement * bookProfile(const char *name, const char *title, int nchX, double lowX, double highX, int nchY, double lowY, double highY, const char *option="s")
Definition: DQMStore.cc:833
MonitorElement * get(const std::string &path) const
get ME from full pathname (e.g. &quot;my/long/dir/my_histo&quot;)
Definition: DQMStore.cc:1270
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:359
static const int NUM_VARS
std::vector< std::vector< double > > param
TProfile * getTProfile(void) const
void MakeLogBinsForProfile(Double_t *bin_edges, const double min, const double max)
std::vector< std::string > axis_titles
TH2F * getTH2F(void) const
void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)
edm::InputTag Muon_Tag
virtual void endJob()
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:237