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