CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
EmDQM.cc
Go to the documentation of this file.
1 // Header file for this //
5 
7 // Collaborating Class Header //
15 //#include "DataFormats/EgammaCandidates/interface/Photon.h"
18 //#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
22 
29 // Root include files //
31 #include "TFile.h"
32 #include "TDirectory.h"
33 #include "TH1F.h"
34 #include <iostream>
35 #include <string>
36 #include <Math/VectorUtil.h>
37 using namespace ROOT::Math::VectorUtil ;
38 
39 
41 // Constructor //
44 {
45 
47  dbe->setVerbose(0);
48 
49 
51  // Read from configuration file //
53  dirname_="HLT/HLTEgammaValidation/"+pset.getParameter<std::string>("@module_label");
54  dbe->setCurrentFolder(dirname_);
55 
56  triggerobjwithrefs = pset.getParameter<edm::InputTag>("triggerobject");
57  // paramters for generator study
58  reqNum = pset.getParameter<unsigned int>("reqNum");
59  pdgGen = pset.getParameter<int>("pdgGen");
60  genEtaAcc = pset.getParameter<double>("genEtaAcc");
61  genEtAcc = pset.getParameter<double>("genEtAcc");
62  // plotting paramters (untracked because they don't affect the physics)
63  plotPtMin = pset.getUntrackedParameter<double>("PtMin",0.);
64  plotPtMax = pset.getUntrackedParameter<double>("PtMax",1000.);
65  plotEtaMax = pset.getUntrackedParameter<double>("EtaMax", 2.7);
66  plotPhiMax = pset.getUntrackedParameter<double>("EtaMax", 3.15);
67  plotBins = pset.getUntrackedParameter<unsigned int>("Nbins",40);
68  useHumanReadableHistTitles = pset.getUntrackedParameter<bool>("useHumanReadableHistTitles", false);
69 
70  //preselction cuts
71  gencutCollection_= pset.getParameter<edm::InputTag>("cutcollection");
72  gencut_ = pset.getParameter<int>("cutnum");
73 
75  // Read in the Vector of Parameter Sets. //
76  // Information for each filter-step //
78  std::vector<edm::ParameterSet> filters =
79  pset.getParameter<std::vector<edm::ParameterSet> >("filters");
80 
81  int i = 0;
82  for(std::vector<edm::ParameterSet>::iterator filterconf = filters.begin() ; filterconf != filters.end() ; filterconf++)
83  {
84 
85  theHLTCollectionLabels.push_back(filterconf->getParameter<edm::InputTag>("HLTCollectionLabels"));
86  theHLTOutputTypes.push_back(filterconf->getParameter<int>("theHLTOutputTypes"));
87  // Grab the human-readable name, if it is not specified, use the Collection Label
88  theHLTCollectionHumanNames.push_back(filterconf->getUntrackedParameter<std::string>("HLTCollectionHumanName",theHLTCollectionLabels[i].label()));
89 
90  std::vector<double> bounds = filterconf->getParameter<std::vector<double> >("PlotBounds");
91  // If the size of plot "bounds" vector != 2, abort
92  assert(bounds.size() == 2);
93  plotBounds.push_back(std::pair<double,double>(bounds[0],bounds[1]));
94  isoNames.push_back(filterconf->getParameter<std::vector<edm::InputTag> >("IsoCollections"));
95  // If the size of the isoNames vector is not greater than zero, abort
96  assert(isoNames.back().size()>0);
97  if (isoNames.back().at(0).label()=="none") {
98  plotiso.push_back(false);
99  } else {
100  plotiso.push_back(true);
101  }
102  i++;
103  } // END of loop over parameter sets
104 
105  // Record number of HLTCollectionLabels
106  numOfHLTCollectionLabels = theHLTCollectionLabels.size();
107 
108 }
109 
110 
112 // method called once each job just before starting event loop //
114 void
116 {
117  //edm::Service<TFileService> fs;
118  dbe->setCurrentFolder(dirname_);
119 
121  // Set up Histogram of Effiency vs Step. //
122  // theHLTCollectionLabels is a vector of InputTags //
123  // from the configuration file. //
125 
126  std::string histName="total_eff";
127  std::string histTitle = "total events passing";
128  // This plot will have bins equal to 2+(number of
129  // HLTCollectionLabels in the config file)
130  total = dbe->book1D(histName.c_str(),histTitle.c_str(),numOfHLTCollectionLabels+2,0,numOfHLTCollectionLabels+2);
131  total->setBinLabel(numOfHLTCollectionLabels+1,"Total");
132  total->setBinLabel(numOfHLTCollectionLabels+2,"Gen");
133  for (unsigned int u=0; u<numOfHLTCollectionLabels; u++){total->setBinLabel(u+1,theHLTCollectionLabels[u].label().c_str());}
134 
135  histName="total_eff_MC_matched";
136  histTitle="total events passing (mc matched)";
137  totalmatch = dbe->book1D(histName.c_str(),histTitle.c_str(),numOfHLTCollectionLabels+2,0,numOfHLTCollectionLabels+2);
138  totalmatch->setBinLabel(numOfHLTCollectionLabels+1,"Total");
139  totalmatch->setBinLabel(numOfHLTCollectionLabels+2,"Gen");
140  for (unsigned int u=0; u<numOfHLTCollectionLabels; u++){totalmatch->setBinLabel(u+1,theHLTCollectionLabels[u].label().c_str());}
141 
142  MonitorElement* tmphisto;
143  MonitorElement* tmpiso;
144 
146  // Set up generator-level histograms //
148  std::string pdgIdString;
149  switch(pdgGen) {
150  case 11:
151  pdgIdString="Electron";break;
152  case 22:
153  pdgIdString="Photon";break;
154  default:
155  pdgIdString="Particle";
156  }
157 
158  histName = "gen_et";
159  histTitle= "E_{T} of " + pdgIdString + "s" ;
160  etgen = dbe->book1D(histName.c_str(),histTitle.c_str(),plotBins,plotPtMin,plotPtMax);
161  histName = "gen_eta";
162  histTitle= "#eta of "+ pdgIdString +"s " ;
163  etagen = dbe->book1D(histName.c_str(),histTitle.c_str(),plotBins,-plotEtaMax,plotEtaMax);
164  histName = "gen_phi";
165  histTitle= "#phi of "+ pdgIdString +"s " ;
166  phigen = dbe->book1D(histName.c_str(),histTitle.c_str(),plotBins,-plotPhiMax,plotPhiMax);
167 
168 
169 
171  // Set up histograms of HLT objects //
173 
174  // Determine what strings to use for histogram titles
175  std::vector<std::string> HltHistTitle;
176  if ( theHLTCollectionHumanNames.size() == numOfHLTCollectionLabels && useHumanReadableHistTitles ) {
177  HltHistTitle = theHLTCollectionHumanNames;
178  } else {
179  for (unsigned int i =0; i < numOfHLTCollectionLabels; i++) {
180  HltHistTitle.push_back(theHLTCollectionLabels[i].label());
181  }
182  }
183 
184  for(unsigned int i = 0; i< numOfHLTCollectionLabels ; i++){
185  // Et distribution of HLT objects passing filter i
186  histName = theHLTCollectionLabels[i].label()+"et_all";
187  histTitle = HltHistTitle[i]+" Et (ALL)";
188  tmphisto = dbe->book1D(histName.c_str(),histTitle.c_str(),plotBins,plotPtMin,plotPtMax);
189  ethist.push_back(tmphisto);
190 
191  // Eta distribution of HLT objects passing filter i
192  histName = theHLTCollectionLabels[i].label()+"eta_all";
193  histTitle = HltHistTitle[i]+" #eta (ALL)";
194  tmphisto = dbe->book1D(histName.c_str(),histTitle.c_str(),plotBins,-plotEtaMax,plotEtaMax);
195  etahist.push_back(tmphisto);
196 
197  // Phi distribution of HLT objects passing filter i
198  histName = theHLTCollectionLabels[i].label()+"phi_all";
199  histTitle = HltHistTitle[i]+" #phi (ALL)";
200  tmphisto = dbe->book1D(histName.c_str(),histTitle.c_str(),plotBins,-plotPhiMax,plotPhiMax);
201  phihist.push_back(tmphisto);
202 
203 
204  // Et distribution of gen object matching HLT object passing filter i
205  histName = theHLTCollectionLabels[i].label()+"et_MC_matched";
206  histTitle = HltHistTitle[i]+" Et (MC matched)";
207  tmphisto = dbe->book1D(histName.c_str(),histTitle.c_str(),plotBins,plotPtMin,plotPtMax);
208  ethistmatch.push_back(tmphisto);
209 
210  // Eta distribution of gen object matching HLT object passing filter i
211  histName = theHLTCollectionLabels[i].label()+"eta_MC_matched";
212  histTitle = HltHistTitle[i]+" #eta (MC matched)";
213  tmphisto = dbe->book1D(histName.c_str(),histTitle.c_str(),plotBins,-plotEtaMax,plotEtaMax);
214  etahistmatch.push_back(tmphisto);
215 
216  // Phi distribution of gen object matching HLT object passing filter i
217  histName = theHLTCollectionLabels[i].label()+"phi_MC_matched";
218  histTitle = HltHistTitle[i]+" #phi (MC matched)";
219  tmphisto = dbe->book1D(histName.c_str(),histTitle.c_str(),plotBins,-plotPhiMax,plotPhiMax);
220  phihistmatch.push_back(tmphisto);
221 
222 
223 
224  // Et distribution of HLT object that is closest delta-R match to sorted gen particle(s)
225  histName = theHLTCollectionLabels[i].label()+"et";
226  histTitle = HltHistTitle[i]+" Et";
227  tmphisto = dbe->book1D(histName.c_str(),histTitle.c_str(),plotBins,plotPtMin,plotPtMax);
228  histEtOfHltObjMatchToGen.push_back(tmphisto);
229 
230  // eta distribution of HLT object that is closest delta-R match to sorted gen particle(s)
231  histName = theHLTCollectionLabels[i].label()+"eta";
232  histTitle = HltHistTitle[i]+" eta";
233  tmphisto = dbe->book1D(histName.c_str(),histTitle.c_str(),plotBins,-plotEtaMax,plotEtaMax);
234  histEtaOfHltObjMatchToGen.push_back(tmphisto);
235 
236  // phi distribution of HLT object that is closest delta-R match to sorted gen particle(s)
237  histName = theHLTCollectionLabels[i].label()+"phi";
238  histTitle = HltHistTitle[i]+" phi";
239  tmphisto = dbe->book1D(histName.c_str(),histTitle.c_str(),plotBins,-plotPhiMax,plotPhiMax);
240  histPhiOfHltObjMatchToGen.push_back(tmphisto);
241 
242 
243 
244  if (!plotiso[i]) {
245  tmpiso = NULL;
246  etahistiso.push_back(tmpiso);
247  phihistiso.push_back(tmpiso);
248  ethistiso.push_back(tmpiso);
249  etahistisomatch.push_back(tmpiso);
250  phihistisomatch.push_back(tmpiso);
251  ethistisomatch.push_back(tmpiso);
252  histEtaIsoOfHltObjMatchToGen.push_back(tmpiso);
253  histPhiIsoOfHltObjMatchToGen.push_back(tmpiso);
254  histEtIsoOfHltObjMatchToGen.push_back(tmpiso);
255  } else {
256  // 2D plot: Isolation values vs eta for all objects
257  histName = theHLTCollectionLabels[i].label()+"eta_isolation_all";
258  histTitle = HltHistTitle[i]+" isolation vs #eta (all)";
259  tmpiso = dbe->book2D(histName.c_str(),histTitle.c_str(),plotBins,-plotEtaMax,plotEtaMax,plotBins,plotBounds[i].first,plotBounds[i].second);
260  etahistiso.push_back(tmpiso);
261 
262  // 2D plot: Isolation values vs phi for all objects
263  histName = theHLTCollectionLabels[i].label()+"phi_isolation_all";
264  histTitle = HltHistTitle[i]+" isolation vs #phi (all)";
265  tmpiso = dbe->book2D(histName.c_str(),histTitle.c_str(),plotBins,-plotPhiMax,plotPhiMax,plotBins,plotBounds[i].first,plotBounds[i].second);
266  phihistiso.push_back(tmpiso);
267 
268 
269  // 2D plot: Isolation values vs et for all objects
270  histName = theHLTCollectionLabels[i].label()+"et_isolation_all";
271  histTitle = HltHistTitle[i]+" isolation vs Et (all)";
272  tmpiso = dbe->book2D(histName.c_str(),histTitle.c_str(),plotBins,plotPtMin,plotPtMax,plotBins,plotBounds[i].first,plotBounds[i].second);
273  ethistiso.push_back(tmpiso);
274 
275  // 2D plot: Isolation values vs eta for matched objects
276  histName = theHLTCollectionLabels[i].label()+"eta_isolation_MC_matched";
277  histTitle = HltHistTitle[i]+" isolation vs #eta (mc matched)";
278  tmpiso = dbe->book2D(histName.c_str(),histTitle.c_str(),plotBins,-plotEtaMax,plotEtaMax,plotBins,plotBounds[i].first,plotBounds[i].second);
279  etahistisomatch.push_back(tmpiso);
280 
281  // 2D plot: Isolation values vs phi for matched objects
282  histName = theHLTCollectionLabels[i].label()+"phi_isolation_MC_matched";
283  histTitle = HltHistTitle[i]+" isolation vs #phi (mc matched)";
284  tmpiso = dbe->book2D(histName.c_str(),histTitle.c_str(),plotBins,-plotPhiMax,plotPhiMax,plotBins,plotBounds[i].first,plotBounds[i].second);
285  phihistisomatch.push_back(tmpiso);
286 
287 
288  // 2D plot: Isolation values vs et for matched objects
289  histName = theHLTCollectionLabels[i].label()+"et_isolation_MC_matched";
290  histTitle = HltHistTitle[i]+" isolation vs Et (mc matched)";
291  tmpiso = dbe->book2D(histName.c_str(),histTitle.c_str(),plotBins,plotPtMin,plotPtMax,plotBins,plotBounds[i].first,plotBounds[i].second);
292  ethistisomatch.push_back(tmpiso);
293 
294 
295  // 2D plot: Isolation values vs eta for HLT object that
296  // is closest delta-R match to sorted gen particle(s)
297  histName = theHLTCollectionLabels[i].label()+"eta_isolation";
298  histTitle = HltHistTitle[i]+" isolation vs #eta";
299  tmpiso = dbe->book2D(histName.c_str(),histTitle.c_str(),plotBins,-plotEtaMax,plotEtaMax,plotBins,plotBounds[i].first,plotBounds[i].second);
300  histEtaIsoOfHltObjMatchToGen.push_back(tmpiso);
301 
302  // 2D plot: Isolation values vs phi for HLT object that
303  // is closest delta-R match to sorted gen particle(s)
304  histName = theHLTCollectionLabels[i].label()+"phi_isolation";
305  histTitle = HltHistTitle[i]+" isolation vs #phi";
306  tmpiso = dbe->book2D(histName.c_str(),histTitle.c_str(),plotBins,-plotPhiMax,plotPhiMax,plotBins,plotBounds[i].first,plotBounds[i].second);
307  histPhiIsoOfHltObjMatchToGen.push_back(tmpiso);
308 
309  // 2D plot: Isolation values vs et for HLT object that
310  // is closest delta-R match to sorted gen particle(s)
311  histName = theHLTCollectionLabels[i].label()+"et_isolation";
312  histTitle = HltHistTitle[i]+" isolation vs Et";
313  tmpiso = dbe->book2D(histName.c_str(),histTitle.c_str(),plotBins,plotPtMin,plotPtMax,plotBins,plotBounds[i].first,plotBounds[i].second);
314  histEtIsoOfHltObjMatchToGen.push_back(tmpiso);
315  } // END of HLT histograms
316 
317  }
318 }
319 
320 
322 // Destructor //
325 }
326 
327 
329 // method called to for each event //
331 void
333 {
334 
336  // Check if there's enough gen particles //
337  // of interest //
340  event.getByLabel(gencutCollection_,cutCounter);
341  if (cutCounter->size() < (unsigned int)gencut_) {
342  //edm::LogWarning("EmDQM") << "Less than "<< gencut_ <<" gen particles with pdgId=" << pdgGen;
343  return;
344  }
345 
346 
347  // fill L1 and HLT info
348  // get objects possed by each filter
350  event.getByLabel(triggerobjwithrefs,triggerObj);
351  if(!triggerObj.isValid()) {
352  edm::LogWarning("EmDQM") << "triggerobjwithrefs=" << triggerobjwithrefs;
353  return;
354  }
355 
356 
358  // Decide if this was an event of interest. //
359  // Did the highest energy particles happen //
360  // to have |eta| < 2.5 ? Then continue. //
363  event.getByLabel("genParticles", genParticles);
364  if(!genParticles.isValid()) {
365  edm::LogWarning("EmDQM") << "genParticles invalid";
366  return;
367  }
368 
369  std::vector<reco::GenParticle> allSortedGenParticles;
370 
371  for(edm::View<reco::GenParticle>::const_iterator currentGenParticle = genParticles->begin(); currentGenParticle != genParticles->end(); currentGenParticle++){
372 
373  if ( !( abs((*currentGenParticle).pdgId())==pdgGen && (*currentGenParticle).status()==1 && (*currentGenParticle).et() > 2.0) ) continue;
374 
375  reco::GenParticle tmpcand( *(currentGenParticle) );
376  allSortedGenParticles.push_back(tmpcand);
377  }
378 
379  std::sort(allSortedGenParticles.begin(), allSortedGenParticles.end(),pTGenComparator_);
380 
381  // Were enough high energy gen particles found?
382  if (allSortedGenParticles.size() < gencut_) {
383  // if no, throw event away
384  return;
385  }
386 
387  // We now have a sorted collection of all generated particles
388  // with pdgId = pdgGen.
389  // Loop over them to see if the top gen particles have eta within acceptance
390  bool keepEvent = true;
391  for (unsigned int i = 0 ; i < gencut_ ; i++ ) {
392  bool inECALgap = fabs(allSortedGenParticles[i].eta()) > 1.4442 && fabs(allSortedGenParticles[i].eta()) < 1.556;
393  if ( (fabs(allSortedGenParticles[i].eta()) > genEtaAcc) || inECALgap ) {
394  //edm::LogWarning("EmDQM") << "Throwing event away. Gen particle with pdgId="<< allSortedGenParticles[i].pdgId() <<"; et="<< allSortedGenParticles[i].et() <<"; and eta="<< allSortedGenParticles[i].eta() <<" beyond acceptance.";
395  keepEvent=false;
396  break;
397  }
398  }
399  if (!keepEvent) {
400  return;
401  }
402 
403 
404  // It was an event worth keeping. Continue.
405 
407  // Fill the bin labeled "Total" //
408  // This will be the number of events looked at. //
410  total->Fill(numOfHLTCollectionLabels+0.5);
411  totalmatch->Fill(numOfHLTCollectionLabels+0.5);
412 
413 
415  // Fill generator info //
417  // the gencut_ highest Et generator objects of the preselected type are our matches
418 
419  std::vector<reco::Particle> sortedGen;
420  for(edm::View<reco::Candidate>::const_iterator genpart = cutCounter->begin(); genpart != cutCounter->end();genpart++){
421  reco::Particle tmpcand( genpart->charge(), genpart->p4(), genpart->vertex(),genpart->pdgId(),genpart->status() );
422  sortedGen.push_back(tmpcand);
423  }
424  std::sort(sortedGen.begin(),sortedGen.end(),pTComparator_ );
425 
426  // Now the collection of gen particles is sorted by pt.
427  // So, remove all particles from the collection so that we
428  // only have the top "1 thru gencut_" particles in it
429  sortedGen.erase(sortedGen.begin()+gencut_,sortedGen.end());
430  if (gencut_ != sortedGen.size() ){
431  return;
432  }
433 
434 
435  for (unsigned int i = 0 ; i < gencut_ ; i++ ) {
436  etgen ->Fill( sortedGen[i].et() ); //validity has been implicitily checked by the cut on gencut_ above
437  etagen->Fill( sortedGen[i].eta() );
438  phigen->Fill( sortedGen[i].phi() );
439  } // END of loop over Generated particles
440  if (gencut_ >= reqNum) total->Fill(numOfHLTCollectionLabels+1.5); // this isn't really needed anymore keep for backward comp.
441  if (gencut_ >= reqNum) totalmatch->Fill(numOfHLTCollectionLabels+1.5); // this isn't really needed anymore keep for backward comp.
442 
443 
445  // Loop over filter modules //
447  for(unsigned int n=0; n < numOfHLTCollectionLabels ; n++) {
448  // These numbers are from the Parameter Set, such as:
449  // theHLTOutputTypes = cms.uint32(100)
450  switch(theHLTOutputTypes[n])
451  {
452  case trigger::TriggerL1NoIsoEG: // Non-isolated Level 1
453  fillHistos<l1extra::L1EmParticleCollection>(triggerObj,event,n,sortedGen);break;
454  case trigger::TriggerL1IsoEG: // Isolated Level 1
455  fillHistos<l1extra::L1EmParticleCollection>(triggerObj,event,n,sortedGen);break;
456  case trigger::TriggerPhoton: // Photon
457  fillHistos<reco::RecoEcalCandidateCollection>(triggerObj,event,n,sortedGen);break;
458  case trigger::TriggerElectron: // Electron
459  fillHistos<reco::ElectronCollection>(triggerObj,event,n,sortedGen);break;
460  case trigger::TriggerCluster: // TriggerCluster
461  fillHistos<reco::RecoEcalCandidateCollection>(triggerObj,event,n,sortedGen);break;
462  default:
463  throw(cms::Exception("Release Validation Error") << "HLT output type not implemented: theHLTOutputTypes[n]" );
464  }
465  } // END of loop over filter modules
466 }
467 
468 
470 // fillHistos //
471 // Called by analyze method. //
473 template <class T> void EmDQM::fillHistos(edm::Handle<trigger::TriggerEventWithRefs>& triggerObj,const edm::Event& iEvent ,unsigned int n,std::vector<reco::Particle>& sortedGen)
474 {
475  std::vector<edm::Ref<T> > recoecalcands;
476  if ( ( triggerObj->filterIndex(theHLTCollectionLabels[n])>=triggerObj->size() )){ // only process if available
477  return;
478  }
479 
481  // Retrieve saved filter objects //
483  triggerObj->getObjects(triggerObj->filterIndex(theHLTCollectionLabels[n]),theHLTOutputTypes[n],recoecalcands);
484  //Danger: special case, L1 non-isolated
485  // needs to be merged with L1 iso
486  if (theHLTOutputTypes[n] == trigger::TriggerL1NoIsoEG){
487  std::vector<edm::Ref<T> > isocands;
488  triggerObj->getObjects(triggerObj->filterIndex(theHLTCollectionLabels[n]),trigger::TriggerL1IsoEG,isocands);
489  if (isocands.size()>0)
490  {
491  for (unsigned int i=0; i < isocands.size(); i++)
492  recoecalcands.push_back(isocands[i]);
493  }
494  } // END of if theHLTOutputTypes == 82
495 
496 
497  if (recoecalcands.size() < 1){ // stop if no object passed the previous filter
498  return;
499  }
500 
501 
502  if (recoecalcands.size() >= reqNum )
503  total->Fill(n+0.5);
504 
505 
507  // check for validity //
508  // prevents crash in CMSSW_3_1_0_pre6 //
510  for (unsigned int j=0; j<recoecalcands.size(); j++){
511  if(!( recoecalcands.at(j).isAvailable())){
512  edm::LogError("EmDQM") << "Event content inconsistent: TriggerEventWithRefs contains invalid Refs" << std::endl << "invalid refs for: " << theHLTCollectionLabels[n].label();
513  return;
514  }
515  }
516 
518  // Loop over the Generated Particles, and find the //
519  // closest HLT object match. //
521  for (unsigned int i=0; i < gencut_; i++) {
522  math::XYZVector currentGenParticleMomentum = sortedGen[i].momentum();
523 
524  float closestDeltaR = 0.5;
525  int closestEcalCandIndex = -1;
526  for (unsigned int j=0; j<recoecalcands.size(); j++) {
527  float deltaR = DeltaR(recoecalcands[j]->momentum(),currentGenParticleMomentum);
528 
529  if (deltaR < closestDeltaR) {
530  closestDeltaR = deltaR;
531  closestEcalCandIndex = j;
532  }
533  }
534 
535  // If an HLT object was found within some delta-R
536  // of this gen particle, store it in a histogram
537  if ( closestEcalCandIndex >= 0 ) {
538  histEtOfHltObjMatchToGen[n] ->Fill( recoecalcands[closestEcalCandIndex]->et() );
539  histEtaOfHltObjMatchToGen[n]->Fill( recoecalcands[closestEcalCandIndex]->eta() );
540  histPhiOfHltObjMatchToGen[n]->Fill( recoecalcands[closestEcalCandIndex]->phi() );
541 
542  // Also store isolation info
543  if (n+1 < numOfHLTCollectionLabels){ // can't plot beyond last
544  if (plotiso[n+1] ){ // only plot if requested in config
545  for (unsigned int j = 0 ; j < isoNames[n+1].size() ;j++ ){
547  iEvent.getByLabel(isoNames[n+1].at(j),depMap);
548  if (depMap.isValid()){ //Map may not exist if only one candidate passes a double filter
549  typename edm::AssociationMap<edm::OneToValue< T , float > >::const_iterator mapi = depMap->find(recoecalcands[closestEcalCandIndex]);
550  if (mapi!=depMap->end()) { // found candidate in isolation map!
551  histEtaIsoOfHltObjMatchToGen[n+1]->Fill( recoecalcands[closestEcalCandIndex]->eta(),mapi->val);
552  histPhiIsoOfHltObjMatchToGen[n+1]->Fill( recoecalcands[closestEcalCandIndex]->phi(),mapi->val);
553  histEtIsoOfHltObjMatchToGen[n+1] ->Fill( recoecalcands[closestEcalCandIndex]->et(), mapi->val);
554  }
555  }
556  }
557  }
558  }
559  } // END of if closestEcalCandIndex >= 0
560  }
561 
562 
563 
565  // Loop over all HLT objects in this filter step, and //
566  // fill histograms. //
568  // bool foundAllMatches = false;
569  // unsigned int numOfHLTobjectsMatched = 0;
570  for (unsigned int i=0; i<recoecalcands.size(); i++) {
571  /*
572  // See if this HLT object has a gen-level match
573  float closestGenParticleDr = 99.0;
574  for(unsigned int j =0; j < gencut_; j++) {
575  math::XYZVector currentGenParticle = sortedGen[j].momentum();
576 
577  double currentDeltaR = DeltaR(recoecalcands[i]->momentum(),currentGenParticle);
578  if ( currentDeltaR < closestGenParticleDr ) {
579  closestGenParticleDr = currentDeltaR;
580  }
581  }
582  // If this HLT object did not have a gen particle match, go to next HLT object
583  if ( !(fabs(closestGenParticleDr < 0.3)) ) continue;
584 
585  numOfHLTobjectsMatched++;
586  if (numOfHLTobjectsMatched >= gencut_) foundAllMatches=true;
587  */
588  // Fill HLT object histograms
589  ethist[n] ->Fill(recoecalcands[i]->et() );
590  etahist[n]->Fill(recoecalcands[i]->eta() );
591  phihist[n]->Fill(recoecalcands[i]->phi() );
592 
594  // Plot isolation variables (show the not-yet-cut //
595  // isolation, i.e. associated to next filter) //
597  if ( n+1 < numOfHLTCollectionLabels ) { // can't plot beyond last
598  if (plotiso[n+1]) {
599  for (unsigned int j = 0 ; j < isoNames[n+1].size() ;j++ ){
601  iEvent.getByLabel(isoNames[n+1].at(j),depMap);
602  if (depMap.isValid()){ //Map may not exist if only one candidate passes a double filter
603  typename edm::AssociationMap<edm::OneToValue< T , float > >::const_iterator mapi = depMap->find(recoecalcands[i]);
604  if (mapi!=depMap->end()){ // found candidate in isolation map!
605  etahistiso[n+1]->Fill(recoecalcands[i]->eta(),mapi->val);
606  phihistiso[n+1]->Fill(recoecalcands[i]->phi(),mapi->val);
607  ethistiso[n+1]->Fill(recoecalcands[i]->et(),mapi->val);
608  }
609  }
610  }
611  }
612  } // END of if n+1 < then the number of hlt collections
613  }
614 
615 
617  // Fill mc matched objects into histograms //
619  unsigned int mtachedMcParts = 0;
620  float mindist=0.3;
621  if(n==0) mindist=0.5; //low L1-resolution => allow wider matching
622  for(unsigned int i =0; i < gencut_; i++){
623  //match generator candidate
624  bool matchThis= false;
625  math::XYZVector candDir=sortedGen[i].momentum();
626  unsigned int closest = 0;
627  double closestDr = 1000.;
628  for(unsigned int trigOb = 0 ; trigOb < recoecalcands.size(); trigOb++){
629  double dr = DeltaR(recoecalcands[trigOb]->momentum(),candDir);
630  if (dr < closestDr) {
631  closestDr = dr;
632  closest = trigOb;
633  }
634  if (closestDr > mindist) { // it's not really a "match" if it's that far away
635  closest = -1;
636  } else {
637  mtachedMcParts++;
638  matchThis = true;
639  }
640  }
641  if ( !matchThis ) continue; // only plot matched candidates
642  // fill coordinates of mc particle matching trigger object
643  ethistmatch[n] ->Fill( sortedGen[i].et() );
644  etahistmatch[n]->Fill( sortedGen[i].eta() );
645  phihistmatch[n]->Fill( sortedGen[i].phi() );
647  // Plot isolation variables (show the not-yet-cut //
648  // isolation, i.e. associated to next filter) //
650  if (n+1 < numOfHLTCollectionLabels){ // can't plot beyond last
651  if (plotiso[n+1] ){ // only plot if requested in config
652  for (unsigned int j = 0 ; j < isoNames[n+1].size() ;j++ ){
654  iEvent.getByLabel(isoNames[n+1].at(j),depMap);
655  if (depMap.isValid()){ //Map may not exist if only one candidate passes a double filter
656  typename edm::AssociationMap<edm::OneToValue< T , float > >::const_iterator mapi = depMap->find(recoecalcands[closest]);
657  if (mapi!=depMap->end()){ // found candidate in isolation map!
658  etahistisomatch[n+1]->Fill(sortedGen[i].eta(),mapi->val);
659  phihistisomatch[n+1]->Fill(sortedGen[i].phi(),mapi->val);
660  ethistisomatch[n+1]->Fill(sortedGen[i].et(),mapi->val);
661  }
662  }
663  }
664  }
665  } // END of if n+1 < then the number of hlt collections
666  }
667  // fill total mc matched efficiency
668  if (mtachedMcParts >= reqNum )
669  totalmatch->Fill(n+0.5);
670 
671 
672 }
673 
674 
676 // method called once each job just after ending the event loop //
679 
680 }
681 
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:86
void beginJob()
Definition: EmDQM.cc:115
void fillHistos(edm::Handle< trigger::TriggerEventWithRefs > &, const edm::Event &, unsigned int, std::vector< reco::Particle > &)
Definition: EmDQM.cc:473
const std::string & label
Definition: MVAComputer.cc:186
~EmDQM()
Destructor.
Definition: EmDQM.cc:324
Definition: deltaR.h:51
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
#define abs(x)
Definition: mlp_lapack.h:159
#define NULL
Definition: scimark2.h:8
Definition: EmDQM.h:22
T eta() const
std::vector< TPRegexp > filters
Definition: eve_filter.cc:25
int iEvent
Definition: GenABIO.cc:243
int j
Definition: DBlmapReader.cc:9
tuple pset
Definition: CrabTask.py:85
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
bool isValid() const
Definition: HandleBase.h:76
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:359
EmDQM(const edm::ParameterSet &pset)
Constructor.
Definition: EmDQM.cc:43
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:31
void analyze(const edm::Event &event, const edm::EventSetup &)
Definition: EmDQM.cc:332
void endJob()
Definition: EmDQM.cc:678
list at
Definition: asciidump.py:428
Definition: DDAxes.h:10