CMS 3D CMS Logo

List of all members | Public Member Functions | Private Member Functions | Private Attributes
MuonMCClassifier Class Reference

#include <MuonAnalysis/MuonAssociators/src/MuonMCClassifier.cc>

Inheritance diagram for MuonMCClassifier:
edm::stream::EDProducer<>

Public Member Functions

 MuonMCClassifier (const edm::ParameterSet &)
 
 ~MuonMCClassifier () override
 
- Public Member Functions inherited from edm::stream::EDProducer<>
 EDProducer ()=default
 
 EDProducer (const EDProducer &)=delete
 
bool hasAbilityToProduceInBeginLumis () const final
 
bool hasAbilityToProduceInBeginProcessBlocks () const final
 
bool hasAbilityToProduceInBeginRuns () const final
 
bool hasAbilityToProduceInEndLumis () const final
 
bool hasAbilityToProduceInEndProcessBlocks () const final
 
bool hasAbilityToProduceInEndRuns () const final
 
const EDProduceroperator= (const EDProducer &)=delete
 

Private Member Functions

int convertAndPush (const TrackingParticle &tp, reco::GenParticleCollection &out, const TrackingParticleRef &momRef, const edm::Handle< reco::GenParticleCollection > &genParticles) const
 
int flavour (int pdgId) const
 Returns the flavour given a pdg id code. More...
 
TrackingParticleRef getTpMother (TrackingParticleRef tp)
 
void produce (edm::Event &, const edm::EventSetup &) override
 
template<typename T >
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. More...
 

Private Attributes

edm::InputTag associatorLabel_
 The Associations. More...
 
double decayAbsZ_
 
double decayRho_
 Cylinder to use to decide if a decay is early or late. More...
 
edm::InputTag genParticles_
 
edm::EDGetTokenT< reco::GenParticleCollectiongenParticlesToken_
 
bool hasMuonCut_
 
bool linkToGenParticles_
 Create a link to the generator level particles. More...
 
edm::EDGetTokenT< reco::MuonToTrackingParticleAssociatormuAssocToken_
 
StringCutObjectSelector< pat::MuonmuonCut_
 
edm::EDGetTokenT< edm::View< reco::Muon > > muonsToken_
 The RECO objects. More...
 
edm::EDGetTokenT< TrackingParticleCollectiontrackingParticlesToken_
 The TrackingParticle objects. More...
 
reco::MuonTrackType trackType_
 Track to use. More...
 

Additional Inherited Members

- Public Types inherited from edm::stream::EDProducer<>
using CacheTypes = CacheContexts< T... >
 
using GlobalCache = typename CacheTypes::GlobalCache
 
using HasAbility = AbilityChecker< T... >
 
using InputProcessBlockCache = typename CacheTypes::InputProcessBlockCache
 
using LuminosityBlockCache = typename CacheTypes::LuminosityBlockCache
 
using LuminosityBlockContext = LuminosityBlockContextT< LuminosityBlockCache, RunCache, GlobalCache >
 
using LuminosityBlockSummaryCache = typename CacheTypes::LuminosityBlockSummaryCache
 
using RunCache = typename CacheTypes::RunCache
 
using RunContext = RunContextT< RunCache, GlobalCache >
 
using RunSummaryCache = typename CacheTypes::RunSummaryCache
 

Detailed Description

CLASSIFICATION: For each RECO Muon, match to SIM particle, and then:

In any case, if the TP is not preferentially matched back to the same RECO muon, label as Ghost (flip the classification)

FLAVOUR:

Definition at line 65 of file MuonMCClassifier.cc.

Constructor & Destructor Documentation

◆ MuonMCClassifier()

MuonMCClassifier::MuonMCClassifier ( const edm::ParameterSet iConfig)
explicit

Definition at line 125 of file MuonMCClassifier.cc.

References Exception, genParticles_, genParticlesToken_, edm::ParameterSet::getParameter(), reco::GlbOrTrk, reco::GlobalTk, reco::InnerTk, linkToGenParticles_, reco::OuterTk, reco::Segments, AlCaHLTBitMon_QueryRunRegistry::string, PbPb_ZMuSkimMuonDPG_cff::trackType, and trackType_.

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") : ""),
130  consumes<TrackingParticleCollection>(iConfig.getParameter<edm::InputTag>("trackingParticles"))),
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 }
edm::InputTag genParticles_
bool linkToGenParticles_
Create a link to the generator level particles.
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
bool existsAs(std::string const &parameterName, bool trackiness=true) const
checks if a parameter exists as a given type
Definition: ParameterSet.h:171
reco::MuonTrackType trackType_
Track to use.
edm::EDGetTokenT< reco::GenParticleCollection > genParticlesToken_
edm::EDGetTokenT< reco::MuonToTrackingParticleAssociator > muAssocToken_
double decayRho_
Cylinder to use to decide if a decay is early or late.
edm::EDGetTokenT< edm::View< reco::Muon > > muonsToken_
The RECO objects.
edm::EDGetTokenT< TrackingParticleCollection > trackingParticlesToken_
The TrackingParticle objects.
StringCutObjectSelector< pat::Muon > muonCut_

◆ ~MuonMCClassifier()

MuonMCClassifier::~MuonMCClassifier ( )
override

Definition at line 186 of file MuonMCClassifier.cc.

186 {}

Member Function Documentation

◆ convertAndPush()

int MuonMCClassifier::convertAndPush ( const TrackingParticle tp,
reco::GenParticleCollection out,
const TrackingParticleRef momRef,
const edm::Handle< reco::GenParticleCollection > &  genParticles 
) const
private

Convert TrackingParticle into GenParticle, save into output collection, if mother is primary set reference to it, return index in output collection

Definition at line 607 of file MuonMCClassifier.cc.

References Exception, AJJGenJetFilter_cfi::genParticles, genParticles_, edm::Ref< C, T, F >::id(), edm::Ref< C, T, F >::isNonnull(), MillePedeFileConverter_cfg::out, and cmsswSequenceInfo::tp.

Referenced by produce().

610  {
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 }
edm::InputTag genParticles_

◆ flavour()

int MuonMCClassifier::flavour ( int  pdgId) const
private

Returns the flavour given a pdg id code.

Definition at line 588 of file MuonMCClassifier.cc.

References funct::abs(), and EgammaValidation_cff::pdgId.

Referenced by produce().

588  {
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 }
Abs< T >::type abs(const T &t)
Definition: Abs.h:22

◆ getTpMother()

TrackingParticleRef MuonMCClassifier::getTpMother ( TrackingParticleRef  tp)
inlineprivate

Definition at line 108 of file MuonMCClassifier.cc.

References cmsswSequenceInfo::tp.

Referenced by produce().

108  {
109  if (tp.isNonnull() && tp->parentVertex().isNonnull() && !tp->parentVertex()->sourceTracks().empty()) {
110  return tp->parentVertex()->sourceTracks()[0];
111  } else {
112  return TrackingParticleRef();
113  }
114  }
edm::Ref< TrackingParticleCollection > TrackingParticleRef

◆ produce()

void MuonMCClassifier::produce ( edm::Event iEvent,
const edm::EventSetup iSetup 
)
overrideprivate

Definition at line 188 of file MuonMCClassifier.cc.

References funct::abs(), reco::MuonToTrackingParticleAssociator::associateMuons(), edm::RefToBaseVector< T >::begin(), convertAndPush(), decayAbsZ_, decayRho_, edm::RefToBaseVector< T >::end(), PVValHelper::eta, Exception, edm::helper::Filler< Map >::fill(), spr::find(), flavour(), EgammaValidation_cff::genp, AJJGenJetFilter_cfi::genParticles, genParticles_, genParticlesToken_, getTpMother(), reco::GlobalTk, hasMuonCut_, mps_fire::i, iEvent, edm::helper::Filler< Map >::insert(), edm::Ref< C, T, F >::isNonnull(), linkToGenParticles_, LogTrace, match(), eostools::move(), amptDefaultParameters_cff::mu, muAssocToken_, muonCut_, PDWG_BPHSkim_cff::muons, muonsToken_, dqmiodumpmetadata::n, reco::OuterTk, phi, edm::Handle< T >::product(), DiDispStaMuonMonitor_cfi::pt, edm::RefToBaseVector< T >::push_back(), edm::RefVector< C, T, F >::push_back(), cmsswSequenceInfo::tp, muonClassificationByHits_cfi::trackingParticles, trackingParticlesToken_, trackType_, and writeValueMap().

188  {
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 }
edm::InputTag genParticles_
bool linkToGenParticles_
Create a link to the generator level particles.
Log< level::Info, true > LogVerbatim
std::map< edm::RefToBase< reco::Muon >, std::vector< std::pair< TrackingParticleRef, double > >, RefToBaseSort > MuonToSimCollection
Definition: MuonTrackType.h:37
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 #
edm::Ref< GenParticleCollection > GenParticleRef
persistent reference to a GenParticle
T const * product() const
Definition: Handle.h:70
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
const_iterator begin() const
int iEvent
Definition: GenABIO.cc:224
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
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.
edm::EDGetTokenT< edm::View< reco::Muon > > muonsToken_
The RECO objects.
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
std::pair< typename Association::data_type::first_type, double > match(Reference key, Association association, bool bestMatchByMaxValue)
Generic matching function.
Definition: Utils.h:10
Log< level::Warning, false > LogWarning
Definition: memstream.h:15
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_

◆ writeValueMap()

template<typename T >
void MuonMCClassifier::writeValueMap ( edm::Event iEvent,
const edm::Handle< edm::View< reco::Muon > > &  handle,
const std::vector< T > &  values,
const std::string &  label 
) const
private

Write a ValueMap<int> in the event.

Definition at line 575 of file MuonMCClassifier.cc.

References trigObjTnPSource_cfi::filler, patZpeak::handle, iEvent, label, eostools::move(), and contentValuesCheck::values.

Referenced by produce().

578  {
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 }
char const * label
int iEvent
Definition: GenABIO.cc:224
HLT enums.
def move(src, dest)
Definition: eostools.py:511

Member Data Documentation

◆ associatorLabel_

edm::InputTag MuonMCClassifier::associatorLabel_
private

The Associations.

Definition at line 87 of file MuonMCClassifier.cc.

◆ decayAbsZ_

double MuonMCClassifier::decayAbsZ_
private

Definition at line 91 of file MuonMCClassifier.cc.

Referenced by produce().

◆ decayRho_

double MuonMCClassifier::decayRho_
private

Cylinder to use to decide if a decay is early or late.

Definition at line 91 of file MuonMCClassifier.cc.

Referenced by produce().

◆ genParticles_

edm::InputTag MuonMCClassifier::genParticles_
private

Definition at line 95 of file MuonMCClassifier.cc.

Referenced by convertAndPush(), MuonMCClassifier(), and produce().

◆ genParticlesToken_

edm::EDGetTokenT<reco::GenParticleCollection> MuonMCClassifier::genParticlesToken_
private

Definition at line 96 of file MuonMCClassifier.cc.

Referenced by MuonMCClassifier(), and produce().

◆ hasMuonCut_

bool MuonMCClassifier::hasMuonCut_
private

A preselection cut for the muon. I pass through pat::Muon so that I can access muon id selectors

Definition at line 77 of file MuonMCClassifier.cc.

Referenced by produce().

◆ linkToGenParticles_

bool MuonMCClassifier::linkToGenParticles_
private

Create a link to the generator level particles.

Definition at line 94 of file MuonMCClassifier.cc.

Referenced by MuonMCClassifier(), and produce().

◆ muAssocToken_

edm::EDGetTokenT<reco::MuonToTrackingParticleAssociator> MuonMCClassifier::muAssocToken_
private

Definition at line 88 of file MuonMCClassifier.cc.

Referenced by produce().

◆ muonCut_

StringCutObjectSelector<pat::Muon> MuonMCClassifier::muonCut_
private

Definition at line 78 of file MuonMCClassifier.cc.

Referenced by produce().

◆ muonsToken_

edm::EDGetTokenT<edm::View<reco::Muon> > MuonMCClassifier::muonsToken_
private

The RECO objects.

Definition at line 73 of file MuonMCClassifier.cc.

Referenced by produce().

◆ trackingParticlesToken_

edm::EDGetTokenT<TrackingParticleCollection> MuonMCClassifier::trackingParticlesToken_
private

The TrackingParticle objects.

Definition at line 84 of file MuonMCClassifier.cc.

Referenced by produce().

◆ trackType_

reco::MuonTrackType MuonMCClassifier::trackType_
private

Track to use.

Definition at line 81 of file MuonMCClassifier.cc.

Referenced by MuonMCClassifier(), and produce().