CMS 3D CMS Logo

MuonMCClassifier.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: MuonMCClassifier
4 // Class: MuonMCClassifier
5 //
28 //
29 // Original Author: G.Petrucciani and G.Abbiendi
30 // Created: Sun Nov 16 16:14:09 CET 2008
31 // revised: 3/Aug/2017
32 //
33 
34 // system include files
35 #include <memory>
36 #include <set>
37 
38 // user include files
41 
46 
48 
52 
57 
60 
62 
63 //
64 // class decleration
66 public:
67  explicit MuonMCClassifier(const edm::ParameterSet &);
68  ~MuonMCClassifier() override;
69 
70 private:
71  void produce(edm::Event &, const edm::EventSetup &) override;
74 
79 
82 
85 
89 
92 
97 
99  int flavour(int pdgId) const;
100 
102  template <typename T>
105  const std::vector<T> &values,
106  const std::string &label) const;
107 
109  if (tp.isNonnull() && tp->parentVertex().isNonnull() && !tp->parentVertex()->sourceTracks().empty()) {
110  return tp->parentVertex()->sourceTracks()[0];
111  } else {
112  return TrackingParticleRef();
113  }
114  }
115 
121  const TrackingParticleRef &momRef,
123 };
124 
126  : muonsToken_(consumes<edm::View<reco::Muon> >(iConfig.getParameter<edm::InputTag>("muons"))),
127  hasMuonCut_(iConfig.existsAs<std::string>("muonPreselection")),
128  muonCut_(hasMuonCut_ ? iConfig.getParameter<std::string>("muonPreselection") : ""),
129  trackingParticlesToken_(
130  consumes<TrackingParticleCollection>(iConfig.getParameter<edm::InputTag>("trackingParticles"))),
131  muAssocToken_(
132  consumes<reco::MuonToTrackingParticleAssociator>(iConfig.getParameter<edm::InputTag>("associatorLabel"))),
133  decayRho_(iConfig.getParameter<double>("decayRho")),
134  decayAbsZ_(iConfig.getParameter<double>("decayAbsZ")),
135  linkToGenParticles_(iConfig.getParameter<bool>("linkToGenParticles")),
136  genParticles_(linkToGenParticles_ ? iConfig.getParameter<edm::InputTag>("genParticles") : edm::InputTag("NONE"))
137 
138 {
139  std::string trackType = iConfig.getParameter<std::string>("trackType");
140  if (trackType == "inner")
142  else if (trackType == "outer")
144  else if (trackType == "global")
146  else if (trackType == "segments")
148  else if (trackType == "glb_or_trk")
150  else
151  throw cms::Exception("Configuration") << "Track type '" << trackType << "' not supported.\n";
152  if (linkToGenParticles_) {
153  genParticlesToken_ = consumes<reco::GenParticleCollection>(genParticles_);
154  }
155 
156  produces<edm::ValueMap<int> >();
157  produces<edm::ValueMap<int> >("ext");
158  produces<edm::ValueMap<int> >("flav");
159  produces<edm::ValueMap<int> >("hitsPdgId");
160  produces<edm::ValueMap<int> >("G4processType"); // Geant process producing the particle
161  produces<edm::ValueMap<int> >("momPdgId");
162  produces<edm::ValueMap<int> >("momFlav");
163  produces<edm::ValueMap<int> >("momStatus");
164  produces<edm::ValueMap<int> >("gmomPdgId");
165  produces<edm::ValueMap<int> >("gmomFlav");
166  produces<edm::ValueMap<int> >("hmomFlav"); // heaviest mother flavour
167  produces<edm::ValueMap<int> >("tpId");
168  produces<edm::ValueMap<int> >("tpEv");
169  produces<edm::ValueMap<int> >("tpBx");
170  produces<edm::ValueMap<float> >("signp");
171  produces<edm::ValueMap<float> >("pt");
172  produces<edm::ValueMap<float> >("eta");
173  produces<edm::ValueMap<float> >("phi");
174  produces<edm::ValueMap<float> >("prodRho");
175  produces<edm::ValueMap<float> >("prodZ");
176  produces<edm::ValueMap<float> >("momRho");
177  produces<edm::ValueMap<float> >("momZ");
178  produces<edm::ValueMap<float> >("tpAssoQuality");
179  if (linkToGenParticles_) {
180  produces<reco::GenParticleCollection>("secondaries");
181  produces<edm::Association<reco::GenParticleCollection> >("toPrimaries");
182  produces<edm::Association<reco::GenParticleCollection> >("toSecondaries");
183  }
184 }
185 
187 
190  iEvent.getByToken(muonsToken_, muons);
191 
194 
196  if (linkToGenParticles_) {
198  }
199 
201  iEvent.getByToken(muAssocToken_, associatorBase);
202  const reco::MuonToTrackingParticleAssociator *assoByHits = associatorBase.product();
203 
204  reco::MuonToSimCollection recSimColl;
205  reco::SimToMuonCollection simRecColl;
206  edm::LogVerbatim("MuonMCClassifier") << "\n ***************************************************************** ";
207  edm::LogVerbatim("MuonMCClassifier") << " RECO MUON association, type: " << trackType_;
208  edm::LogVerbatim("MuonMCClassifier") << " ***************************************************************** \n";
209 
211  if (!hasMuonCut_) {
212  // all muons
213  for (size_t i = 0, n = muons->size(); i < n; ++i) {
214  edm::RefToBase<reco::Muon> rmu = muons->refAt(i);
215  selMuons.push_back(rmu);
216  }
217  } else {
218  // filter, fill refvectors, associate
219  // I pass through pat::Muon so that I can access muon id selectors
220  for (size_t i = 0, n = muons->size(); i < n; ++i) {
221  edm::RefToBase<reco::Muon> rmu = muons->refAt(i);
222  if (muonCut_(pat::Muon(rmu)))
223  selMuons.push_back(rmu);
224  }
225  }
226 
228  for (size_t i = 0, n = trackingParticles->size(); i < n; ++i) {
230  }
231 
232  assoByHits->associateMuons(recSimColl, simRecColl, selMuons, trackType_, allTPs);
233 
234  // for global muons without hits on muon detectors, look at the linked standalone muon
235  reco::MuonToSimCollection UpdSTA_recSimColl;
236  reco::SimToMuonCollection UpdSTA_simRecColl;
237  if (trackType_ == reco::GlobalTk) {
238  edm::LogVerbatim("MuonMCClassifier") << "\n ***************************************************************** ";
239  edm::LogVerbatim("MuonMCClassifier") << " STANDALONE (UpdAtVtx) MUON association ";
240  edm::LogVerbatim("MuonMCClassifier") << " ***************************************************************** \n";
241  assoByHits->associateMuons(UpdSTA_recSimColl, UpdSTA_simRecColl, selMuons, reco::OuterTk, allTPs);
242  }
243 
244  typedef reco::MuonToSimCollection::const_iterator r2s_it;
245  typedef reco::SimToMuonCollection::const_iterator s2r_it;
246 
247  size_t nmu = muons->size();
248  edm::LogVerbatim("MuonMCClassifier") << "\n There are " << nmu << " reco::Muons.";
249 
250  std::vector<int> classif(nmu, 0), ext(nmu, 0);
251  std::vector<int> hitsPdgId(nmu, 0), G4processType(nmu, 0), momPdgId(nmu, 0), gmomPdgId(nmu, 0), momStatus(nmu, 0);
252  std::vector<int> flav(nmu, 0), momFlav(nmu, 0), gmomFlav(nmu, 0), hmomFlav(nmu, 0);
253  std::vector<int> tpId(nmu, -1), tpBx(nmu, 999), tpEv(nmu, 999);
254  std::vector<float> signp(nmu, 0.0), pt(nmu, 0.0), eta(nmu, -99.), phi(nmu, -99.);
255  std::vector<float> prodRho(nmu, 0.0), prodZ(nmu, 0.0), momRho(nmu, 0.0), momZ(nmu, 0.0);
256  std::vector<float> tpAssoQuality(nmu, -1);
257 
258  std::unique_ptr<reco::GenParticleCollection> secondaries; // output collection of secondary muons
259  std::map<TrackingParticleRef, int> tpToSecondaries; // map from tp to (index+1) in output collection
260  std::vector<int> muToPrimary(nmu, -1), muToSecondary(nmu, -1); // map from input into (index) in output, -1 for null
262  secondaries = std::make_unique<reco::GenParticleCollection>();
263 
264  // loop on reco muons
265  for (size_t i = 0; i < nmu; ++i) {
267  if (hasMuonCut_ && (std::find(selMuons.begin(), selMuons.end(), mu) == selMuons.end())) {
268  LogTrace("MuonMCClassifier") << "\n reco::Muon # " << i
269  << "didn't pass the selection. classified as -99 and skipped";
270  classif[i] = -99;
271  continue;
272  } else
273  edm::LogVerbatim("MuonMCClassifier") << "\n reco::Muon # " << i;
274 
276  edm::RefToBase<reco::Muon> muMatchBack;
277  r2s_it match = recSimColl.find(mu);
278  s2r_it matchback;
279  if (match != recSimColl.end()) {
280  // match->second is vector, front is first element, first is the ref (second would be the quality)
281  tp = match->second.front().first;
282  tpId[i] = tp.isNonnull() ? tp.key() : -1; // we check, even if null refs should not appear here at all
283  tpAssoQuality[i] = match->second.front().second;
284  s2r_it matchback = simRecColl.find(tp);
285  if (matchback != simRecColl.end()) {
286  muMatchBack = matchback->second.front().first;
287  } else {
288  edm::LogWarning("MuonMCClassifier") << "\n***WARNING: This I do NOT understand: why no match back? *** \n";
289  }
290  } else if ((trackType_ == reco::GlobalTk) && mu->isGlobalMuon()) {
291  // perform a second attempt, matching with the standalone muon
292  r2s_it matchSta = UpdSTA_recSimColl.find(mu);
293  if (matchSta != UpdSTA_recSimColl.end()) {
294  tp = matchSta->second.front().first;
295  tpId[i] = tp.isNonnull() ? tp.key() : -1; // we check, even if null refs should not appear here at all
296  tpAssoQuality[i] = matchSta->second.front().second;
297  s2r_it matchback = UpdSTA_simRecColl.find(tp);
298  if (matchback != UpdSTA_simRecColl.end()) {
299  muMatchBack = matchback->second.front().first;
300  } else {
301  edm::LogWarning("MuonMCClassifier")
302  << "\n***WARNING: This I do NOT understand: why no match back in UpdSTA? *** \n";
303  }
304  }
305  } else {
306  edm::LogVerbatim("MuonMCClassifier") << "\t No matching TrackingParticle is found ";
307  }
308 
309  if (tp.isNonnull()) {
310  bool isGhost = muMatchBack != mu;
311  if (isGhost)
312  edm::LogVerbatim("MuonMCClassifier") << "\t *** This seems a Duplicate muon ! classif[i] will be < 0 ***";
313 
314  // identify signal and pileup TP
315  tpBx[i] = tp->eventId().bunchCrossing();
316  tpEv[i] = tp->eventId().event();
317 
318  hitsPdgId[i] = tp->pdgId();
319  prodRho[i] = tp->vertex().Rho();
320  prodZ[i] = tp->vertex().Z();
321 
322  // added info on GEANT process producing the TrackingParticle
323  const std::vector<SimVertex> &G4Vs = tp->parentVertex()->g4Vertices();
324  G4processType[i] = G4Vs[0].processType();
325 
326  signp[i] = tp->charge() * tp->p();
327  pt[i] = tp->pt();
328  eta[i] = tp->eta();
329  phi[i] = tp->phi();
330 
331  // Try to extract mother and grand mother of this muon.
332  // Unfortunately, SIM and GEN histories require diffent code :-(
333  if (!tp->genParticles().empty()) { // Muon is in GEN
334  reco::GenParticleRef genp = tp->genParticles()[0];
335  reco::GenParticleRef genMom = genp->numberOfMothers() > 0 ? genp->motherRef() : reco::GenParticleRef();
336  reco::GenParticleRef mMom = genMom;
337 
338  if (genMom.isNonnull()) {
339  if (genMom->pdgId() != tp->pdgId()) {
340  momPdgId[i] = genMom->pdgId();
341  momStatus[i] = genMom->status();
342  momRho[i] = genMom->vertex().Rho();
343  momZ[i] = genMom->vz();
344  } else {
345  // if mother has the same identity look backwards for the real mother (it may happen in radiative decays)
346  int jm = 0;
347  while (mMom->pdgId() == tp->pdgId()) {
348  jm++;
349 
350  if (mMom->numberOfMothers() > 0) {
351  mMom = mMom->motherRef();
352  }
353  LogTrace("MuonMCClassifier") << "\t\t backtracking mother " << jm << ", pdgId = " << mMom->pdgId()
354  << ", status= " << mMom->status();
355  }
356  genMom = mMom; // redefine genMom
357  momPdgId[i] = genMom->pdgId();
358  momStatus[i] = genMom->status();
359  momRho[i] = genMom->vertex().Rho();
360  momZ[i] = genMom->vz();
361  }
362  edm::LogVerbatim("MuonMCClassifier")
363  << "\t Particle pdgId = " << hitsPdgId[i] << ", (Event,Bx) = "
364  << "(" << tpEv[i] << "," << tpBx[i] << ")"
365  << "\n\t q*p = " << signp[i] << ", pT = " << pt[i] << ", eta = " << eta[i] << ", phi = " << phi[i]
366  << "\n\t produced at vertex rho = " << prodRho[i] << ", z = " << prodZ[i]
367  << ", (GEANT4 process = " << G4processType[i] << ")"
368  << "\n\t has GEN mother pdgId = " << momPdgId[i] << " (status = " << momStatus[i] << ")";
369 
370  reco::GenParticleRef genGMom = genMom->numberOfMothers() > 0 ? genMom->motherRef() : reco::GenParticleRef();
371 
372  if (genGMom.isNonnull()) {
373  gmomPdgId[i] = genGMom->pdgId();
374  edm::LogVerbatim("MuonMCClassifier") << "\t\t mother prod. vertex rho = " << momRho[i]
375  << ", z = " << momZ[i] << ", grand-mom pdgId = " << gmomPdgId[i];
376  }
377  // in this case, we might want to know the heaviest mom flavour
378  for (reco::GenParticleRef nMom = genMom;
379  nMom.isNonnull() &&
380  std::abs(nMom->pdgId()) >= 100; // stop when we're no longer looking at hadrons or mesons
381  nMom = nMom->numberOfMothers() > 0 ? nMom->motherRef() : reco::GenParticleRef()) {
382  int flav = flavour(nMom->pdgId());
383  if (hmomFlav[i] < flav)
384  hmomFlav[i] = flav;
385  edm::LogVerbatim("MuonMCClassifier") << "\t\t backtracking flavour: mom pdgId = " << nMom->pdgId()
386  << ", flavour = " << flav << ", heaviest so far = " << hmomFlav[i];
387  }
388  } // if (genMom.isNonnull())
389 
390  else { // mother is null ??
391  edm::LogWarning("MuonMCClassifier")
392  << "\t GenParticle with pdgId = " << hitsPdgId[i] << ", (Event,Bx) = "
393  << "(" << tpEv[i] << "," << tpBx[i] << ")"
394  << "\n\t q*p = " << signp[i] << ", pT = " << pt[i] << ", eta = " << eta[i] << ", phi = " << phi[i]
395  << "\n\t produced at vertex rho = " << prodRho[i] << ", z = " << prodZ[i]
396  << ", (GEANT4 process = " << G4processType[i] << ")"
397  << "\n\t has NO mother!";
398  }
399 
400  } else { // Muon is in SIM Only
402  if (simMom.isNonnull()) {
403  momPdgId[i] = simMom->pdgId();
404  momRho[i] = simMom->vertex().Rho();
405  momZ[i] = simMom->vertex().Z();
406  edm::LogVerbatim("MuonMCClassifier")
407  << "\t Particle pdgId = " << hitsPdgId[i] << ", (Event,Bx) = "
408  << "(" << tpEv[i] << "," << tpBx[i] << ")"
409  << "\n\t q*p = " << signp[i] << ", pT = " << pt[i] << ", eta = " << eta[i] << ", phi = " << phi[i]
410  << "\n\t produced at vertex rho = " << prodRho[i] << ", z = " << prodZ[i]
411  << ", (GEANT4 process = " << G4processType[i] << ")"
412  << "\n\t has SIM mother pdgId = " << momPdgId[i] << " produced at rho = " << simMom->vertex().Rho()
413  << ", z = " << simMom->vertex().Z();
414 
415  if (!simMom->genParticles().empty()) {
416  momStatus[i] = simMom->genParticles()[0]->status();
417  reco::GenParticleRef genGMom =
418  (simMom->genParticles()[0]->numberOfMothers() > 0 ? simMom->genParticles()[0]->motherRef()
420  if (genGMom.isNonnull())
421  gmomPdgId[i] = genGMom->pdgId();
422  edm::LogVerbatim("MuonMCClassifier")
423  << "\t\t SIM mother is in GEN (status " << momStatus[i] << "), grand-mom id = " << gmomPdgId[i];
424  } else {
425  momStatus[i] = -1;
426  TrackingParticleRef simGMom = getTpMother(simMom);
427  if (simGMom.isNonnull())
428  gmomPdgId[i] = simGMom->pdgId();
429  edm::LogVerbatim("MuonMCClassifier") << "\t\t SIM mother is in SIM only, grand-mom id = " << gmomPdgId[i];
430  }
431  } else {
432  edm::LogVerbatim("MuonMCClassifier")
433  << "\t Particle pdgId = " << hitsPdgId[i] << ", (Event,Bx) = "
434  << "(" << tpEv[i] << "," << tpBx[i] << ")"
435  << "\n\t q*p = " << signp[i] << ", pT = " << pt[i] << ", eta = " << eta[i] << ", phi = " << phi[i]
436  << "\n\t produced at vertex rho = " << prodRho[i] << ", z = " << prodZ[i]
437  << ", (GEANT4 process = " << G4processType[i] << ")"
438  << "\n\t has NO mother!";
439  }
440  }
441  momFlav[i] = flavour(momPdgId[i]);
442  gmomFlav[i] = flavour(gmomPdgId[i]);
443 
444  // Check first IF this is a muon at all
445  if (abs(tp->pdgId()) != 13) {
446  if (abs(tp->pdgId()) == 11) {
447  classif[i] = isGhost ? -11 : 11;
448  ext[i] = isGhost ? -11 : 11;
449  edm::LogVerbatim("MuonMCClassifier") << "\t This is electron/positron. classif[i] = " << classif[i];
450  } else {
451  classif[i] = isGhost ? -1 : 1;
452  ext[i] = isGhost ? -1 : 1;
453  edm::LogVerbatim("MuonMCClassifier") << "\t This is not a muon. Sorry. classif[i] = " << classif[i];
454  }
455 
456  continue;
457  }
458 
459  // Is this SIM muon also a GEN muon, with a mother?
460  if (!tp->genParticles().empty() && (momPdgId[i] != 0)) {
461  if (abs(momPdgId[i]) < 100 && (abs(momPdgId[i]) != 15)) {
462  classif[i] = isGhost ? -4 : 4;
463  flav[i] = 13;
464  ext[i] = 10;
465  edm::LogVerbatim("MuonMCClassifier") << "\t This seems PRIMARY MUON ! classif[i] = " << classif[i];
466  } else if (momFlav[i] == 4 || momFlav[i] == 5 || momFlav[i] == 15) {
467  classif[i] = isGhost ? -3 : 3;
468  flav[i] = momFlav[i];
469  if (momFlav[i] == 15)
470  ext[i] = 9; // tau->mu
471  else if (momFlav[i] == 5)
472  ext[i] = 8; // b->mu
473  else if (hmomFlav[i] == 5)
474  ext[i] = 7; // b->c->mu or b->tau->mu
475  else
476  ext[i] = 6; // c->mu
477  edm::LogVerbatim("MuonMCClassifier") << "\t This seems HEAVY FLAVOUR ! classif[i] = " << classif[i];
478  } else {
479  classif[i] = isGhost ? -2 : 2;
480  flav[i] = momFlav[i];
481  edm::LogVerbatim("MuonMCClassifier") << "\t This seems LIGHT FLAVOUR ! classif[i] = " << classif[i];
482  }
483  } else {
484  classif[i] = isGhost ? -2 : 2;
485  flav[i] = momFlav[i];
486  edm::LogVerbatim("MuonMCClassifier") << "\t This seems LIGHT FLAVOUR ! classif[i] = " << classif[i];
487  }
488  // extended classification
489  if (momPdgId[i] == 0)
490  ext[i] = 2; // if it has no mom, it's not a primary particle so it won't be in ppMuX
491  else if (abs(momPdgId[i]) < 100)
492  ext[i] = (momFlav[i] == 15 ? 9 : 10); // primary mu, or tau->mu
493  else if (momFlav[i] == 5)
494  ext[i] = 8; // b->mu
495  else if (momFlav[i] == 4)
496  ext[i] = (hmomFlav[i] == 5 ? 7 : 6); // b->c->mu and c->mu
497  else if (momStatus[i] != -1) { // primary light particle
498  int id = std::abs(momPdgId[i]);
499  if (id != /*pi+*/ 211 && id != /*K+*/ 321 && id != 130 /*K0L*/)
500  ext[i] = 5; // other light particle, possibly short-lived
501  else if (prodRho[i] < decayRho_ && std::abs(prodZ[i]) < decayAbsZ_)
502  ext[i] = 4; // decay a la ppMuX (primary pi/K within a cylinder)
503  else
504  ext[i] = 3; // late decay that wouldn't be in ppMuX
505  } else
506  ext[i] = 2; // decay of non-primary particle, would not be in ppMuX
507  if (isGhost)
508  ext[i] = -ext[i];
509 
510  if (linkToGenParticles_ && std::abs(ext[i]) >= 2) {
511  // Link to the genParticle if possible, but not decays in flight (in ppMuX they're in GEN block, but they have wrong parameters)
512  if (!tp->genParticles().empty() && std::abs(ext[i]) >= 5) {
513  if (genParticles.id() != tp->genParticles().id()) {
514  throw cms::Exception("Configuration")
515  << "Product ID mismatch between the genParticle collection (" << genParticles_ << ", id "
516  << genParticles.id() << ") and the references in the TrackingParticles (id " << tp->genParticles().id()
517  << ")\n";
518  }
519  muToPrimary[i] = tp->genParticles()[0].key();
520  } else {
521  // Don't put the same trackingParticle twice!
522  int &indexPlus1 = tpToSecondaries[tp]; // will create a 0 if the tp is not in the list already
523  if (indexPlus1 == 0)
524  indexPlus1 = convertAndPush(*tp, *secondaries, getTpMother(tp), genParticles) + 1;
525  muToSecondary[i] = indexPlus1 - 1;
526  }
527  }
528  edm::LogVerbatim("MuonMCClassifier") << "\t Extended classification code = " << ext[i];
529  } // if (tp.isNonnull())
530  } // end loop on reco muons
531 
532  writeValueMap(iEvent, muons, classif, "");
533  writeValueMap(iEvent, muons, ext, "ext");
534  writeValueMap(iEvent, muons, flav, "flav");
535  writeValueMap(iEvent, muons, tpId, "tpId");
536  writeValueMap(iEvent, muons, tpBx, "tpBx");
537  writeValueMap(iEvent, muons, tpEv, "tpEv");
538  writeValueMap(iEvent, muons, hitsPdgId, "hitsPdgId");
539  writeValueMap(iEvent, muons, G4processType, "G4processType");
540  writeValueMap(iEvent, muons, momPdgId, "momPdgId");
541  writeValueMap(iEvent, muons, momStatus, "momStatus");
542  writeValueMap(iEvent, muons, momFlav, "momFlav");
543  writeValueMap(iEvent, muons, gmomPdgId, "gmomPdgId");
544  writeValueMap(iEvent, muons, gmomFlav, "gmomFlav");
545  writeValueMap(iEvent, muons, hmomFlav, "hmomFlav");
546  writeValueMap(iEvent, muons, signp, "signp");
547  writeValueMap(iEvent, muons, pt, "pt");
548  writeValueMap(iEvent, muons, eta, "eta");
549  writeValueMap(iEvent, muons, phi, "phi");
550  writeValueMap(iEvent, muons, prodRho, "prodRho");
551  writeValueMap(iEvent, muons, prodZ, "prodZ");
552  writeValueMap(iEvent, muons, momRho, "momRho");
553  writeValueMap(iEvent, muons, momZ, "momZ");
554  writeValueMap(iEvent, muons, tpAssoQuality, "tpAssoQuality");
555 
556  if (linkToGenParticles_) {
557  edm::OrphanHandle<reco::GenParticleCollection> secHandle = iEvent.put(std::move(secondaries), "secondaries");
560  std::unique_ptr<edm::Association<reco::GenParticleCollection> > outPri(
562  std::unique_ptr<edm::Association<reco::GenParticleCollection> > outSec(
564  edm::Association<reco::GenParticleCollection>::Filler fillPri(*outPri), fillSec(*outSec);
565  fillPri.insert(muons, muToPrimary.begin(), muToPrimary.end());
566  fillSec.insert(muons, muToSecondary.begin(), muToSecondary.end());
567  fillPri.fill();
568  fillSec.fill();
569  iEvent.put(std::move(outPri), "toPrimaries");
570  iEvent.put(std::move(outSec), "toSecondaries");
571  }
572 }
573 
574 template <typename T>
577  const std::vector<T> &values,
578  const std::string &label) const {
579  using namespace edm;
580  using namespace std;
581  unique_ptr<ValueMap<T> > valMap(new ValueMap<T>());
582  typename edm::ValueMap<T>::Filler filler(*valMap);
583  filler.insert(handle, values.begin(), values.end());
584  filler.fill();
585  iEvent.put(std::move(valMap), label);
586 }
587 
589  int flav = std::abs(pdgId);
590  // for quarks, leptons and bosons except gluons, take their pdgId
591  // muons and taus have themselves as flavour
592  if (flav <= 37 && flav != 21)
593  return flav;
594  // look for barions
595  int bflav = ((flav / 1000) % 10);
596  if (bflav != 0)
597  return bflav;
598  // look for mesons
599  int mflav = ((flav / 100) % 10);
600  if (mflav != 0)
601  return mflav;
602  return 0;
603 }
604 
605 // push secondary in collection.
606 // if it has a primary mother link to it.
609  const TrackingParticleRef &simMom,
611  out.push_back(reco::GenParticle(tp.charge(), tp.p4(), tp.vertex(), tp.pdgId(), tp.status(), true));
612  if (simMom.isNonnull() && !simMom->genParticles().empty()) {
613  if (genParticles.id() != simMom->genParticles().id()) {
614  throw cms::Exception("Configuration")
615  << "Product ID mismatch between the genParticle collection (" << genParticles_ << ", id " << genParticles.id()
616  << ") and the references in the TrackingParticles (id " << simMom->genParticles().id() << ")\n";
617  }
618  out.back().addMother(simMom->genParticles()[0]);
619  }
620  return out.size() - 1;
621 }
622 
623 //define this as a plug-in
edm::InputTag genParticles_
bool linkToGenParticles_
Create a link to the generator level particles.
Log< level::Info, true > LogVerbatim
std::vector< GenParticle > GenParticleCollection
collection of GenParticles
std::map< edm::RefToBase< reco::Muon >, std::vector< std::pair< TrackingParticleRef, double > >, RefToBaseSort > MuonToSimCollection
Definition: MuonTrackType.h:37
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
void associateMuons(MuonToSimCollection &recoToSim, SimToMuonCollection &simToReco, const edm::RefToBaseVector< reco::Muon > &muons, MuonTrackType type, const edm::RefVector< TrackingParticleCollection > &tpColl) const
genp
produce generated paricles in acceptance #
ProductID id() const
Accessor for product ID.
Definition: Ref.h:244
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
edm::Ref< GenParticleCollection > GenParticleRef
persistent reference to a GenParticle
T const * product() const
Definition: Handle.h:70
void insert(const H &h, I begin, I end)
Definition: ValueMap.h:53
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:238
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
#define LogTrace(id)
reco::MuonTrackType trackType_
Track to use.
std::map< TrackingParticleRef, std::vector< std::pair< edm::RefToBase< reco::Muon >, double > > > SimToMuonCollection
Definition: MuonTrackType.h:38
~MuonMCClassifier() override
char const * label
const_iterator begin() const
int iEvent
Definition: GenABIO.cc:224
Definition: Muon.py:1
int convertAndPush(const TrackingParticle &tp, reco::GenParticleCollection &out, const TrackingParticleRef &momRef, const edm::Handle< reco::GenParticleCollection > &genParticles) const
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
TrackingParticleRef getTpMother(TrackingParticleRef tp)
const_iterator end() const
MuonMCClassifier(const edm::ParameterSet &)
edm::InputTag associatorLabel_
The Associations.
edm::EDGetTokenT< reco::GenParticleCollection > genParticlesToken_
void writeValueMap(edm::Event &iEvent, const edm::Handle< edm::View< reco::Muon > > &handle, const std::vector< T > &values, const std::string &label) const
Write a ValueMap<int> in the event.
edm::EDGetTokenT< reco::MuonToTrackingParticleAssociator > muAssocToken_
int flavour(int pdgId) const
Returns the flavour given a pdg id code.
double decayRho_
Cylinder to use to decide if a decay is early or late.
fixed size matrix
edm::EDGetTokenT< edm::View< reco::Muon > > muonsToken_
The RECO objects.
HLT enums.
void push_back(const RefToBase< T > &)
void push_back(value_type const &ref)
Add a Ref<C, T> to the RefVector.
Definition: RefVector.h:67
Monte Carlo truth information used for tracking validation.
std::pair< typename Association::data_type::first_type, double > match(Reference key, Association association, bool bestMatchByMaxValue)
Generic matching function.
Definition: Utils.h:10
std::vector< TrackingParticle > TrackingParticleCollection
Log< level::Warning, false > LogWarning
MuonTrackType
Definition: MuonTrackType.h:28
Definition: memstream.h:15
void produce(edm::Event &, const edm::EventSetup &) override
Analysis-level muon class.
Definition: Muon.h:51
edm::Ref< TrackingParticleCollection > TrackingParticleRef
def move(src, dest)
Definition: eostools.py:511
edm::EDGetTokenT< TrackingParticleCollection > trackingParticlesToken_
The TrackingParticle objects.
StringCutObjectSelector< pat::Muon > muonCut_