CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MuonPFAnalyzer.cc
Go to the documentation of this file.
1 
7 //Base class
10 
12 
15 
16 //System included files
17 #include <memory>
18 #include <string>
19 #include <typeinfo>
20 #include <utility>
21 
22 //Root included files
23 #include "TH1.h"
24 #include "TH2.h"
25 #include "TProfile.h"
26 
27 
28 //Event framework included files
32 
37 
38 using namespace edm;
39 using namespace std;
40 using namespace reco;
41 
42 
44 
45  LogTrace("MuonPFAnalyzer") << "[MuonPFAnalyzer] Initializing configuration from parameterset.\n";
46 
47  theGenLabel_ = consumes<GenParticleCollection>(pSet.getParameter<InputTag>("inputTagGenParticles"));
48  theRecoLabel_ = consumes<MuonCollection> (pSet.getParameter<InputTag>("inputTagMuonReco"));
49  theVertexLabel_ = consumes<VertexCollection> (pSet.getParameter<InputTag>("inputTagVertex"));
50  theBeamSpotLabel_ = consumes<BeamSpot> (pSet.getParameter<InputTag>("inputTagBeamSpot"));
51 
52  theHighPtTh = pSet.getParameter<double>("highPtThreshold");
53  theRecoGenR = pSet.getParameter<double>("recoGenDeltaR");
54  theIsoCut = pSet.getParameter<double>("relCombIsoCut");
55  theRunOnMC = pSet.getParameter<bool>("runOnMC");
56 
57  theFolder = pSet.getParameter<string>("folder");
58 
59  theMuonKinds.push_back(""); // all TUNEP/PF muons
60  theMuonKinds.push_back("Tight"); // tight TUNEP/PF muons
61  theMuonKinds.push_back("TightIso"); // tight/iso TUNEP/PF muons
62 
63 
64 }
65 
66 
68 {
69 
70  LogTrace("MuonPFAnalyzer") <<
71  "[MuonPFAnalyzer] Destructor called.\n";
72 
73 }
74 
75 
76 
77 // ------------ method called when starting to processes a run ------------
78 
80 
81  LogTrace("MuonPFAnalyzer") <<
82  "[MuonPFAnalyzer] Booking histograms.\n";
83 
84  //Set up DAQ
85  theDbe = 0;
86  theDbe = edm::Service<DQMStore>().operator->();
87  theDbe->cd();
88 
89  if(theRunOnMC)
90  {
91  bookHistos("PF");
92  bookHistos("PFTight");
93  bookHistos("PFTightIso");
94  bookHistos("TUNEP");
95  bookHistos("TUNEPTight");
96  bookHistos("TUNEPTightIso");
97  }
98 
99  bookHistos("PFvsTUNEP");
100  bookHistos("PFvsTUNEPTight");
101  bookHistos("PFvsTUNEPTightIso");
102 
103 
104 }
105 
107  const EventSetup& context)
108 {
109 
111  event.getByToken(theRecoLabel_, muons);
112 
114  event.getByToken(theGenLabel_, genMuons);
115 
117  event.getByToken(theBeamSpotLabel_, beamSpot);
118 
120  event.getByToken(theVertexLabel_, vertex);
121 
122  const Vertex primaryVertex = getPrimaryVertex(vertex, beamSpot);
123 
124  recoToGenMatch(muons, genMuons);
125 
126  RecoGenCollection::const_iterator recoGenIt = theRecoGen.begin();
127  RecoGenCollection::const_iterator recoGenEnd = theRecoGen.end();
128 
129  for (;recoGenIt!=recoGenEnd;++recoGenIt)
130  {
131 
132  const Muon *muon = recoGenIt->first;
133  TrackRef tunePTrack = muon->tunePMuonBestTrack();
134 
135  const GenParticle *genMuon = recoGenIt->second;
136 
137  vector<string>::const_iterator kindIt = theMuonKinds.begin();
138  vector<string>::const_iterator kindEnd = theMuonKinds.end();
139 
140  for (;kindIt!=kindEnd;++kindIt)
141  {
142 
143  const string & kind = (*kindIt);
144 
145  if (kind.find("Tight") != string::npos &&
146  !muon::isTightMuon((*muon), primaryVertex)) continue;
147 
148  if (kind.find("Iso") != string::npos &&
149  combRelIso(muon) > theIsoCut) continue;
150 
151  if (theRunOnMC && genMuon && !muon->innerTrack().isNull() ) // has matched gen muon
152  {
153 
154  if (!tunePTrack.isNull())
155  {
156 
157  string group = "TUNEP" + kind;
158 
159  float pt = tunePTrack->pt();
160  float phi = tunePTrack->phi();
161  float eta = tunePTrack->eta();
162 
163  float genPt = genMuon->pt();
164  float genPhi = genMuon->p4().phi();
165  float genEta = genMuon->p4().eta();
166 
167  float dPtOverPt = (pt / genPt) - 1;
168 
169  if (pt < theHighPtTh)
170  {
171 
172  fillInRange(getPlot(group,"code"),1,muonTrackType(muon, false));
173  fillInRange(getPlot(group,"deltaPtOverPt"),1,dPtOverPt);
174  }
175  else
176  {
177  fillInRange(getPlot(group,"codeHighPt"),1,muonTrackType(muon, false));
178  fillInRange(getPlot(group,"deltaPtOverPtHighPt"),1,dPtOverPt);
179  }
180 
181  fillInRange(getPlot(group,"deltaPt"),1,(pt - genPt));
182  fillInRange(getPlot(group,"deltaPhi"),1,fDeltaPhi(genPhi,phi));
183  fillInRange(getPlot(group,"deltaEta"),1,genEta - eta);
184 
185  }
186 
187  if (muon->isPFMuon())
188  {
189 
190  string group = "PF" + kind;
191 
192  // Assumes that default in muon is PF
193  float pt = muon->pt();
194  float phi = muon->p4().phi();
195  float eta = muon->p4().eta();
196 
197  float genPt = genMuon->pt();
198  float genPhi = genMuon->p4().phi();
199  float genEta = genMuon->p4().eta();
200 
201  float dPtOverPt = (pt / genPt) - 1;
202 
203  if (pt < theHighPtTh)
204  {
205  fillInRange(getPlot(group,"code"),1,muonTrackType(muon, true));
206  fillInRange(getPlot(group,"deltaPtOverPt"),1,dPtOverPt);
207  }
208  else
209  {
210  fillInRange(getPlot(group,"codeHighPt"),1,muonTrackType(muon, true));
211  fillInRange(getPlot(group,"deltaPtOverPtHighPt"),1,dPtOverPt);
212  }
213 
214 
215  fillInRange(getPlot(group,"deltaPt"),1,pt - genPt);
216  fillInRange(getPlot(group,"deltaPhi"),1,fDeltaPhi(genPhi,phi));
217  fillInRange(getPlot(group,"deltaEta"),1,genEta - eta);
218 
219  }
220 
221  }
222 
223 
224 
225  if (muon->isPFMuon() && !tunePTrack.isNull() &&
226  !muon->innerTrack().isNull()) // Compare PF with TuneP + Tracker
227  { // No gen matching needed
228 
229  string group = "PFvsTUNEP" + kind;
230 
231  float pt = tunePTrack->pt();
232  float phi = tunePTrack->phi();
233  float eta = tunePTrack->eta();
234 
235  // Assumes that default in muon is PF
236  float pfPt = muon->pt();
237  float pfPhi = muon->p4().phi();
238  float pfEta = muon->p4().eta();
239  float dPtOverPt = (pfPt / pt) - 1; // TUNEP vs PF pt used as denum.
240 
241 
242  if (pt < theHighPtTh)
243  {
244  fillInRange(getPlot(group,"code"),2,
245  muonTrackType(muon, false),muonTrackType(muon, true));
246  fillInRange(getPlot(group,"deltaPtOverPt"),1,dPtOverPt);
247  }
248  else
249  {
250  fillInRange(getPlot(group,"codeHighPt"),
251  2,muonTrackType(muon, false),muonTrackType(muon, true));
252  fillInRange(getPlot(group,"deltaPtOverPtHighPt"),1,dPtOverPt);
253  }
254 
255  fillInRange(getPlot(group,"deltaPt"),1,pfPt - pt);
256  fillInRange(getPlot(group,"deltaPhi"),1,fDeltaPhi(pfPhi,phi));
257  fillInRange(getPlot(group,"deltaEta"),1,pfEta - eta);
258 
259 
260  if (theRunOnMC && genMuon) // has a matched gen muon
261 
262  {
263 
264  float genPt = genMuon->pt();
265  float dPtOverPtGen = (pt / genPt) - 1;
266  float dPtOverPtGenPF = (pfPt / genPt) - 1;
267 
268  if (pt < theHighPtTh)
269  {
270  fillInRange(getPlot(group,"deltaPtOverPtPFvsTUNEP"),
271  2,dPtOverPtGen,dPtOverPtGenPF);
272  }
273  else
274  {
275  fillInRange(getPlot(group,"deltaPtOverPtHighPtPFvsTUNEP"),
276  2,dPtOverPtGen,dPtOverPtGenPF);
277  }
278  }
279 
280  }
281 
282  }
283 
284  }
285 
286 }
287 
288 
289 
290 
291 void MuonPFAnalyzer::bookHistos(const string & group) {
292 
293 
294 
295  LogTrace("MuonPFAnalyzer") << "[MuonPFAnalyzer] Booking histos for group :"
296  << group << "\n";
297 
298  theDbe->setCurrentFolder(string(theFolder) + group);
299 
300 
301  bool isPFvsTUNEP = group.find("PFvsTUNEP") != string::npos;
302 
303  string hName;
304 
305 
306  hName = "deltaPtOverPt" + group;
307  thePlots[group]["deltaPtOverPt"] = theDbe->book1D(hName.c_str(),hName.c_str(),101,-1.01,1.01);
308 
309  hName = "deltaPtOverPtHighPt" + group;
310  thePlots[group]["deltaPtOverPtHighPt"] = theDbe->book1D(hName.c_str(),hName.c_str(),101,-1.01,1.01);
311 
312  hName = "deltaPt" + group;
313  thePlots[group]["deltaPt"] = theDbe->book1D(hName.c_str(),hName.c_str(),201.,-10.25,10.25);
314 
315  hName = "deltaPhi"+group;
316  thePlots[group]["deltaPhi"] = theDbe->book1D(hName.c_str(),hName.c_str(),51.,0,.0102);
317 
318  hName = "deltaEta"+group;
319  thePlots[group]["deltaEta"] = theDbe->book1D(hName.c_str(),hName.c_str(),101.,-.00505,.00505);
320 
321 
322 
323  if (isPFvsTUNEP) {
324 
325 
326  hName = "code"+group;
327  MonitorElement * plot = theDbe->book2D(hName.c_str(),hName.c_str(),7,-.5,6.5,7,-.5,6.5);
328  thePlots[group]["code"] = plot;
329  setCodeLabels(plot,1);
330  setCodeLabels(plot,2);
331 
332  hName = "codeHighPt"+group;
333  plot = theDbe->book2D(hName.c_str(),hName.c_str(),7,-.5,6.5,7,-.5,6.5);
334  thePlots[group]["codeHighPt"] = plot;
335  setCodeLabels(plot,1);
336  setCodeLabels(plot,2);
337 
338 
339  if (theRunOnMC)
340  {
341  hName = "deltaPtOverPtPFvsTUNEP" + group;
342  thePlots[group]["deltaPtOverPtPFvsTUNEP"] =
343  theDbe->book2D(hName.c_str(),hName.c_str(),
344  101,-1.01,1.01,101,-1.01,1.01);
345 
346  hName = "deltaPtOverPtHighPtPFvsTUNEP" + group;
347  thePlots[group]["deltaPtOverPtHighPtPFvsTUNEP"] =
348  theDbe->book2D(hName.c_str(),hName.c_str(),
349  101,-1.01,1.01,101,-1.01,1.01);
350  }
351  } else {
352  hName = "code"+group;
353  MonitorElement * plot = theDbe->book1D(hName.c_str(),hName.c_str(),7,-.5,6.5);
354  thePlots[group]["code"] = plot;
355  setCodeLabels(plot,1);
356 
357  hName = "codeHighPt"+group;
358  plot = theDbe->book1D(hName.c_str(),hName.c_str(),7,-.5,6.5);
359  thePlots[group]["codeHighPt"] = plot;
360  setCodeLabels(plot,1);
361  }
362 
363 }
364 
365 
367  const string & type) {
368 
369  map<string,map<string,MonitorElement *> >::iterator groupIt = thePlots.find(group);
370  if (groupIt == thePlots.end()) {
371  LogTrace("MuonPFAnalyzer") << "[MuonPFAnalyzer] GROUP : " << group
372  << " is not a valid plot group. Returning 0.\n";
373  return 0;
374  }
375 
376  map<string,MonitorElement *>::iterator typeIt = groupIt->second.find(type);
377  if (typeIt == groupIt->second.end()) {
378  LogTrace("MuonPFAnalyzer") << "[MuonPFAnalyzer] TYPE : " << type
379  << " is not a valid type for GROUP : " << group
380  << ". Returning 0.\n";
381  return 0;
382  }
383 
384  return typeIt->second;
385 
386 }
387 
388 
390 {
391 
392  MuonIsolation iso = muon->isolationR03();
393  float combRelIso = (iso.emEt + iso.hadEt + iso.sumPt) / muon->pt();
394 
395  return combRelIso;
396 
397 }
398 
399 
400 inline float MuonPFAnalyzer::fDeltaPhi(float phi1, float phi2) {
401 
402  float fPhiDiff = fabs(acos(cos(phi1-phi2)));
403  return fPhiDiff;
404 
405 }
406 
407 
409 {
410 
411  TAxis *axis = 0;
412 
413  TH1 * histo = plot->getTH1();
414  if(!histo) return;
415 
416  if (nAxis==1)
417  axis =histo->GetXaxis();
418  else if (nAxis == 2)
419  axis =histo->GetYaxis();
420 
421  if(!axis) return;
422 
423  axis->SetBinLabel(1,"Inner Track");
424  axis->SetBinLabel(2,"Outer Track");
425  axis->SetBinLabel(3,"Combined");
426  axis->SetBinLabel(4,"TPFMS");
427  axis->SetBinLabel(5,"Picky");
428  axis->SetBinLabel(6,"DYT");
429  axis->SetBinLabel(7,"None");
430 
431 }
432 
433 
434 void MuonPFAnalyzer::fillInRange(MonitorElement *plot, int nAxis, double x, double y)
435 {
436 
437  TH1 * histo = plot->getTH1();
438 
439  TAxis *axis[2] = {0, 0};
440  axis[0] = histo->GetXaxis();
441  if (nAxis == 2)
442  axis[1] = histo->GetYaxis();
443 
444  double value[2] = {0, 0};
445  value[0] = x;
446  value[1] = y;
447 
448  for (int i = 0;i<nAxis;++i)
449  {
450  double min = axis[i]->GetXmin();
451  double max = axis[i]->GetXmax();
452 
453  if (value[i] <= min)
454  value[i] = axis[i]->GetBinCenter(1);
455 
456  if (value[i] >= max)
457  value[i] = axis[i]->GetBinCenter(axis[i]->GetNbins());
458  }
459 
460  if (nAxis == 2)
461  plot->Fill(value[0],value[1]);
462  else
463  plot->Fill(value[0]);
464 
465 }
466 
467 
468 int MuonPFAnalyzer::muonTrackType(const Muon * muon, bool usePF) {
469 
470  switch ( usePF ? muon->muonBestTrackType() : muon->tunePMuonBestTrackType() ) {
471  case Muon::InnerTrack :
472  return 0;
473  case Muon::OuterTrack :
474  return 1;
475  case Muon::CombinedTrack :
476  return 2;
477  case Muon::TPFMS :
478  return 3;
479  case Muon::Picky :
480  return 4;
481  case Muon::DYT :
482  return 5;
483  case Muon::None :
484  return 6;
485  }
486 
487  return 6;
488 
489 }
490 
491 
494 {
495 
496  theRecoGen.clear();
497 
498  if (muons.isValid())
499  {
500 
501  MuonCollection::const_iterator muonIt = muons->begin();
502  MuonCollection::const_iterator muonEnd = muons->end();
503 
504  for(; muonIt!=muonEnd; ++muonIt)
505  {
506 
507  float bestDR = 999.;
508  const GenParticle *bestGen = 0;
509 
510  if (theRunOnMC && gens.isValid())
511  {
512 
513  GenParticleCollection::const_iterator genIt = gens->begin();
514  GenParticleCollection::const_iterator genEnd = gens->end();
515 
516  for(; genIt!=genEnd; ++genIt)
517  {
518 
519  if (abs(genIt->pdgId()) == 13 )
520  {
521 
522  float muonPhi = muonIt->phi();
523  float muonEta = muonIt->eta();
524 
525  float genPhi = genIt->phi();
526  float genEta = genIt->eta();
527 
528  float dR = deltaR(muonEta,muonPhi,
529  genEta,genPhi);
530 
531  if (dR < theRecoGenR && dR < bestDR)
532  {
533  bestDR = dR;
534  bestGen = &(*genIt);
535  }
536 
537  }
538 
539  }
540  }
541 
542  theRecoGen.push_back(RecoGenPair(&(*muonIt), bestGen));
543 
544  }
545  }
546 
547 }
548 
551 {
552 
553  Vertex::Point posVtx;
554  Vertex::Error errVtx;
555 
556  bool hasPrimaryVertex = false;
557 
558  if (vertex.isValid())
559  {
560 
561  vector<Vertex>::const_iterator vertexIt = vertex->begin();
562  vector<Vertex>::const_iterator vertexEnd = vertex->end();
563 
564  for (;vertexIt!=vertexEnd;++vertexIt)
565  {
566  if (vertexIt->isValid() &&
567  !vertexIt->isFake())
568  {
569  posVtx = vertexIt->position();
570  errVtx = vertexIt->error();
571  hasPrimaryVertex = true;
572  break;
573  }
574  }
575  }
576 
577  if ( !hasPrimaryVertex ) {
578 
579  LogInfo("MuonPFAnalyzer") <<
580  "[MuonPFAnalyzer] PrimaryVertex not found, use BeamSpot position instead.\n";
581 
582  posVtx = beamSpot->position();
583  errVtx(0,0) = beamSpot->BeamWidthX();
584  errVtx(1,1) = beamSpot->BeamWidthY();
585  errVtx(2,2) = beamSpot->sigmaZ();
586 
587  }
588 
589  const Vertex primaryVertex(posVtx,errVtx);
590 
591  return primaryVertex;
592 
593 }
594 
595 
596 //define this as a plug-in
int muonTrackType(const reco::Muon *muon, bool usePF)
type
Definition: HCALResponse.h:21
float hadEt
hcal sum-Et
Definition: MuonIsolation.h:9
T getParameter(std::string const &) const
int i
Definition: DBlmapReader.cc:9
const reco::Vertex getPrimaryVertex(edm::Handle< reco::VertexCollection > &vertex, edm::Handle< reco::BeamSpot > &beamSpot)
~MuonPFAnalyzer()
Destructor.
float sumPt
sum-pt of tracks
Definition: MuonIsolation.h:7
virtual MuonTrackType muonBestTrackType() const
Definition: Muon.h:64
virtual const LorentzVector & p4() const GCC11_FINAL
four-momentum Lorentz vector
virtual TrackRef innerTrack() const
Definition: Muon.h:48
math::Error< dimension >::type Error
covariance error matrix (3x3)
Definition: Vertex.h:43
DEFINE_FWK_MODULE(HiMixingModule)
T eta() const
void setCodeLabels(MonitorElement *plot, int nAxis)
void Fill(long long x)
void bookHistos()
Definition: Histogram.h:33
MonitorElement * getPlot(const std::string &group, const std::string &type)
bool isNull() const
Checks for null.
Definition: Ref.h:247
virtual void analyze(const edm::Event &, const edm::EventSetup &)
Perform the PF - TUNEP muon analysis.
void fillInRange(MonitorElement *plot, int nAxis, double x, double y=0)
virtual const_iterator begin() const
first daughter const_iterator
float combRelIso(const reco::Muon *muon)
float emEt
ecal sum-Et
Definition: MuonIsolation.h:8
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
TH1 * getTH1(void) const
math::XYZPoint Point
point in the space
Definition: Vertex.h:39
EventID const & min(EventID const &lh, EventID const &rh)
Definition: EventID.h:132
void bookHistos(const std::string &group)
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
#define LogTrace(id)
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
float fDeltaPhi(float phi1, float phi2)
virtual void beginRun(edm::Run const &, edm::EventSetup const &)
Initialize an book plots.
virtual MuonTrackType tunePMuonBestTrackType() const
Definition: Muon.h:66
tuple muons
Definition: patZpeak.py:38
void recoToGenMatch(edm::Handle< reco::MuonCollection > &reco, edm::Handle< reco::GenParticleCollection > &gen)
std::pair< const reco::Muon *, const reco::GenParticle * > RecoGenPair
bool isPFMuon() const
Definition: Muon.h:222
virtual TrackRef tunePMuonBestTrack() const
Definition: Muon.h:65
Definition: DDAxes.h:10
bool isTightMuon(const reco::Muon &, const reco::Vertex &)
virtual float pt() const GCC11_FINAL
transverse momentum
MuonPFAnalyzer(const edm::ParameterSet &)
Constructor.
def fillInRange
Definition: dataDML.py:89
const MuonIsolation & isolationR03() const
Definition: Muon.h:158
Definition: Run.h:41
EventID const & max(EventID const &lh, EventID const &rh)
Definition: EventID.h:137
Definition: DDAxes.h:10