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