Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:30:46 

0001 // -*- C++ -*-
0002 //
0003 // Package:    SimMuon/MCTruth
0004 // Class:      MuonSimClassifier
0005 //
0006 /*
0007 
0008 
0009  CLASSIFICATION: For each RECO Muon, match to SIM particle, and then:
0010   - If the SIM is not a Muon, label as Punchthrough (1) except if it is an
0011  electron or positron (11)
0012   - If the SIM is a Muon, then look at it's provenance.
0013      A) the SIM muon is also a GEN muon, whose parent is NOT A HADRON AND NOT A
0014  TAU
0015         -> classify as "primary" (4).
0016      B) the SIM muon is also a GEN muon, whose parent is HEAVY FLAVOURED HADRON
0017  OR A TAU
0018         -> classify as "heavy flavour" (3)
0019      C) classify as "light flavour/decay" (2)
0020 
0021   In any case, if the TP is not preferentially matched back to the same RECO
0022  muon, label as Ghost (flip the classification)
0023 
0024 
0025  FLAVOUR:
0026   - for non-muons: 0
0027   - for primary muons: 13
0028   - for non primary muons: flavour of the mother: std::abs(pdgId) of heaviest
0029  quark, or 15 for tau
0030 
0031 */
0032 //
0033 // Original Author:  G.Petrucciani and G.Abbiendi
0034 //         Created:  Sun Nov 16 16:14:09 CET 2008
0035 //         revised:  3/Aug/2017
0036 //
0037 
0038 // system include files
0039 #include <memory>
0040 #include <set>
0041 
0042 // user include files
0043 #include "FWCore/Framework/interface/Frameworkfwd.h"
0044 #include "FWCore/Framework/interface/stream/EDProducer.h"
0045 
0046 #include "FWCore/Framework/interface/ESHandle.h"
0047 #include "FWCore/Framework/interface/Event.h"
0048 #include "FWCore/Framework/interface/EventSetup.h"
0049 #include "FWCore/Framework/interface/MakerMacros.h"
0050 
0051 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0052 
0053 #include "DataFormats/Common/interface/Association.h"
0054 #include "DataFormats/Common/interface/ValueMap.h"
0055 #include "DataFormats/Common/interface/View.h"
0056 
0057 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0058 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
0059 #include "DataFormats/MuonReco/interface/Muon.h"
0060 #include "DataFormats/MuonReco/interface/MuonSimInfo.h"
0061 
0062 #include "SimDataFormats/Associations/interface/MuonToTrackingParticleAssociator.h"
0063 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
0064 
0065 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0066 
0067 //
0068 // class decleration
0069 class MuonSimClassifier : public edm::stream::EDProducer<> {
0070 public:
0071   explicit MuonSimClassifier(const edm::ParameterSet &);
0072   ~MuonSimClassifier() override;
0073 
0074 private:
0075   void produce(edm::Event &, const edm::EventSetup &) override;
0076   /// The RECO objects
0077   edm::EDGetTokenT<edm::View<reco::Muon>> muonsToken_;
0078 
0079   /// Track to use
0080   reco::MuonTrackType trackType_;
0081 
0082   /// The TrackingParticle objects
0083   edm::EDGetTokenT<TrackingParticleCollection> trackingParticlesToken_;
0084 
0085   /// The Associations
0086   edm::InputTag associatorLabel_;
0087   edm::EDGetTokenT<reco::MuonToTrackingParticleAssociator> muAssocToken_;
0088 
0089   /// Cylinder to use to decide if a decay is early or late
0090   double decayRho_, decayAbsZ_;
0091 
0092   /// Create a link to the generator level particles
0093   bool linkToGenParticles_;
0094   edm::InputTag genParticles_;
0095   edm::EDGetTokenT<reco::GenParticleCollection> genParticlesToken_;
0096 
0097   /// Returns the flavour given a pdg id code
0098   int flavour(int pdgId) const;
0099 
0100   /// Write a ValueMap<int> in the event
0101   template <typename T>
0102   void writeValueMap(edm::Event &iEvent,
0103                      const edm::Handle<edm::View<reco::Muon>> &handle,
0104                      const std::vector<T> &values,
0105                      const std::string &label) const;
0106 
0107   TrackingParticleRef getTpMother(TrackingParticleRef tp) {
0108     if (tp.isNonnull() && tp->parentVertex().isNonnull() && !tp->parentVertex()->sourceTracks().empty()) {
0109       return tp->parentVertex()->sourceTracks()[0];
0110     } else {
0111       return TrackingParticleRef();
0112     }
0113   }
0114 
0115   /// Convert TrackingParticle into GenParticle, save into output collection,
0116   /// if mother is primary set reference to it,
0117   /// return index in output collection
0118   int convertAndPush(const TrackingParticle &tp,
0119                      reco::GenParticleCollection &out,
0120                      const TrackingParticleRef &momRef,
0121                      const edm::Handle<reco::GenParticleCollection> &genParticles) const;
0122 };
0123 
0124 MuonSimClassifier::MuonSimClassifier(const edm::ParameterSet &iConfig)
0125     : muonsToken_(consumes<edm::View<reco::Muon>>(iConfig.getParameter<edm::InputTag>("muons"))),
0126       trackingParticlesToken_(
0127           consumes<TrackingParticleCollection>(iConfig.getParameter<edm::InputTag>("trackingParticles"))),
0128       muAssocToken_(
0129           consumes<reco::MuonToTrackingParticleAssociator>(iConfig.getParameter<edm::InputTag>("associatorLabel"))),
0130       decayRho_(iConfig.getParameter<double>("decayRho")),
0131       decayAbsZ_(iConfig.getParameter<double>("decayAbsZ")),
0132       linkToGenParticles_(iConfig.getParameter<bool>("linkToGenParticles")),
0133       genParticles_(linkToGenParticles_ ? iConfig.getParameter<edm::InputTag>("genParticles") : edm::InputTag())
0134 
0135 {
0136   std::string trackType = iConfig.getParameter<std::string>("trackType");
0137   if (trackType == "inner")
0138     trackType_ = reco::InnerTk;
0139   else if (trackType == "outer")
0140     trackType_ = reco::OuterTk;
0141   else if (trackType == "global")
0142     trackType_ = reco::GlobalTk;
0143   else if (trackType == "segments")
0144     trackType_ = reco::Segments;
0145   else if (trackType == "glb_or_trk")
0146     trackType_ = reco::GlbOrTrk;
0147   else
0148     throw cms::Exception("Configuration") << "Track type '" << trackType << "' not supported.\n";
0149   if (linkToGenParticles_) {
0150     genParticlesToken_ = consumes<reco::GenParticleCollection>(genParticles_);
0151   }
0152 
0153   produces<edm::ValueMap<reco::MuonSimInfo>>();
0154   if (linkToGenParticles_) {
0155     produces<reco::GenParticleCollection>("secondaries");
0156     produces<edm::Association<reco::GenParticleCollection>>("toPrimaries");
0157     produces<edm::Association<reco::GenParticleCollection>>("toSecondaries");
0158   }
0159 }
0160 
0161 MuonSimClassifier::~MuonSimClassifier() {}
0162 
0163 void dumpFormatedInfo(const reco::MuonSimInfo &simInfo) {
0164   return;
0165   LogTrace("MuonSimClassifier") << "\t Particle pdgId = " << simInfo.pdgId << ", (Event,Bx) = "
0166                                 << "(" << simInfo.tpEvent << "," << simInfo.tpBX << ")"
0167                                 << "\n\t   q*p = " << simInfo.charge * simInfo.p4.P() << ", pT = " << simInfo.p4.pt()
0168                                 << ", eta = " << simInfo.p4.eta() << ", phi = " << simInfo.p4.phi()
0169                                 << "\n\t   produced at vertex rho = " << simInfo.vertex.Rho()
0170                                 << ", z = " << simInfo.vertex.Z() << ", (GEANT4 process = " << simInfo.g4processType
0171                                 << ")\n";
0172 }
0173 
0174 void MuonSimClassifier::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0175   edm::Handle<edm::View<reco::Muon>> muons;
0176   iEvent.getByToken(muonsToken_, muons);
0177 
0178   edm::Handle<TrackingParticleCollection> trackingParticles;
0179   iEvent.getByToken(trackingParticlesToken_, trackingParticles);
0180 
0181   edm::Handle<reco::GenParticleCollection> genParticles;
0182   if (linkToGenParticles_) {
0183     iEvent.getByToken(genParticlesToken_, genParticles);
0184   }
0185 
0186   edm::Handle<reco::MuonToTrackingParticleAssociator> associatorBase;
0187   iEvent.getByToken(muAssocToken_, associatorBase);
0188   const reco::MuonToTrackingParticleAssociator *assoByHits = associatorBase.product();
0189 
0190   reco::MuonToSimCollection recSimColl;
0191   reco::SimToMuonCollection simRecColl;
0192   LogTrace("MuonSimClassifier") << "\n "
0193                                    "***************************************************************** ";
0194   LogTrace("MuonSimClassifier") << " RECO MUON association, type:  " << trackType_;
0195   LogTrace("MuonSimClassifier") << " ******************************************"
0196                                    "*********************** \n";
0197 
0198   edm::RefToBaseVector<reco::Muon> allMuons;
0199   for (size_t i = 0, n = muons->size(); i < n; ++i) {
0200     allMuons.push_back(muons->refAt(i));
0201   }
0202 
0203   edm::RefVector<TrackingParticleCollection> allTPs;
0204   for (size_t i = 0, n = trackingParticles->size(); i < n; ++i) {
0205     allTPs.push_back(TrackingParticleRef(trackingParticles, i));
0206   }
0207 
0208   assoByHits->associateMuons(recSimColl, simRecColl, allMuons, trackType_, allTPs);
0209 
0210   // for global muons without hits on muon detectors, look at the linked
0211   // standalone muon
0212   reco::MuonToSimCollection updSTA_recSimColl;
0213   reco::SimToMuonCollection updSTA_simRecColl;
0214   if (trackType_ == reco::GlobalTk) {
0215     LogTrace("MuonSimClassifier") << "\n "
0216                                      "***************************************************************** ";
0217     LogTrace("MuonSimClassifier") << " STANDALONE (UpdAtVtx) MUON association ";
0218     LogTrace("MuonSimClassifier") << " ****************************************"
0219                                      "************************* \n";
0220     assoByHits->associateMuons(updSTA_recSimColl, updSTA_simRecColl, allMuons, reco::OuterTk, allTPs);
0221   }
0222 
0223   typedef reco::MuonToSimCollection::const_iterator r2s_it;
0224   typedef reco::SimToMuonCollection::const_iterator s2r_it;
0225 
0226   size_t nmu = muons->size();
0227   LogTrace("MuonSimClassifier") << "\n There are " << nmu << " reco::Muons.";
0228 
0229   std::vector<reco::MuonSimInfo> simInfo;
0230 
0231   std::unique_ptr<reco::GenParticleCollection> secondaries;  // output collection of secondary muons
0232   std::map<TrackingParticleRef, int> tpToSecondaries;        // map from tp to (index+1) in output collection
0233   std::vector<int> muToPrimary(nmu, -1), muToSecondary(nmu,
0234                                                        -1);  // map from input into (index) in output, -1 for null
0235   if (linkToGenParticles_)
0236     secondaries = std::make_unique<reco::GenParticleCollection>();
0237 
0238   // loop on reco muons
0239   for (size_t i = 0; i < nmu; ++i) {
0240     simInfo.push_back(reco::MuonSimInfo());
0241     LogTrace("MuonSimClassifier") << "\n reco::Muon # " << i;
0242 
0243     TrackingParticleRef tp;
0244     edm::RefToBase<reco::Muon> muMatchBack;
0245     r2s_it match = recSimColl.find(allMuons.at(i));
0246     s2r_it matchback;
0247     if (match != recSimColl.end()) {
0248       // match->second is vector, front is first element, first is the ref
0249       // (second would be the quality)
0250       tp = match->second.front().first;
0251       simInfo[i].tpId = tp.isNonnull() ? tp.key() : -1;  // we check, even if null refs should not appear here at all
0252       simInfo[i].tpAssoQuality = match->second.front().second;
0253       s2r_it matchback = simRecColl.find(tp);
0254       if (matchback != simRecColl.end()) {
0255         muMatchBack = matchback->second.front().first;
0256       } else {
0257         LogTrace("MuonSimClassifier") << "\n***WARNING:  This I do NOT understand: why no match back? "
0258                                          "*** \n";
0259       }
0260     } else {
0261       if ((trackType_ == reco::GlobalTk) && allMuons.at(i)->isGlobalMuon()) {
0262         // perform a second attempt, matching with the standalone muon
0263         r2s_it matchSta = updSTA_recSimColl.find(allMuons.at(i));
0264         if (matchSta != updSTA_recSimColl.end()) {
0265           tp = matchSta->second.front().first;
0266           simInfo[i].tpId = tp.isNonnull() ? tp.key() : -1;  // we check, even if null refs
0267                                                              // should not appear here at all
0268           simInfo[i].tpAssoQuality = matchSta->second.front().second;
0269           s2r_it matchback = updSTA_simRecColl.find(tp);
0270           if (matchback != updSTA_simRecColl.end()) {
0271             muMatchBack = matchback->second.front().first;
0272           } else {
0273             LogTrace("MuonSimClassifier") << "\n***WARNING:  This I do NOT understand: why no match back "
0274                                              "in updSTA? *** \n";
0275           }
0276         }
0277       } else {
0278         LogTrace("MuonSimClassifier") << "\t No matching TrackingParticle is found ";
0279       }
0280     }
0281 
0282     if (tp.isNonnull()) {
0283       bool isGhost = muMatchBack != allMuons.at(i);
0284       if (isGhost)
0285         LogTrace("MuonSimClassifier") << "\t *** This seems a Duplicate muon ! "
0286                                          "classif[i] will be < 0 ***";
0287 
0288       // identify signal and pileup TP
0289       simInfo[i].tpBX = tp->eventId().bunchCrossing();
0290       simInfo[i].tpEvent = tp->eventId().event();
0291 
0292       simInfo[i].pdgId = tp->pdgId();
0293       simInfo[i].vertex = tp->vertex();
0294 
0295       // added info on GEANT process producing the TrackingParticle
0296       const std::vector<SimVertex> &g4Vs = tp->parentVertex()->g4Vertices();
0297       simInfo[i].g4processType = g4Vs[0].processType();
0298 
0299       simInfo[i].charge = tp->charge();
0300       simInfo[i].p4 = tp->p4();
0301 
0302       // Try to extract mother and grand mother of this muon.
0303       // Unfortunately, SIM and GEN histories require diffent code :-(
0304       if (!tp->genParticles().empty()) {  // Muon is in GEN
0305         reco::GenParticleRef genp = tp->genParticles()[0];
0306         reco::GenParticleRef genMom = genp->numberOfMothers() > 0 ? genp->motherRef() : reco::GenParticleRef();
0307         reco::GenParticleRef mMom = genMom;
0308 
0309         if (genMom.isNonnull()) {
0310           if (genMom->pdgId() != tp->pdgId()) {
0311             simInfo[i].motherPdgId = genMom->pdgId();
0312             simInfo[i].motherStatus = genMom->status();
0313             simInfo[i].motherVertex = genMom->vertex();
0314           } else {
0315             // if mother has the same identity look backwards for the real
0316             // mother (it may happen in radiative decays)
0317             int jm = 0;
0318             while (mMom->pdgId() == tp->pdgId()) {
0319               jm++;
0320               if (mMom->numberOfMothers() > 0) {
0321                 mMom = mMom->motherRef();
0322               } else {
0323                 LogTrace("MuonSimClassifier") << "\t No Mother is found ";
0324                 break;
0325               }
0326 
0327               LogTrace("MuonSimClassifier") << "\t\t backtracking mother " << jm << ", pdgId = " << mMom->pdgId()
0328                                             << ", status= " << mMom->status();
0329             }
0330             genMom = mMom;  // redefine genMom
0331             simInfo[i].motherPdgId = genMom->pdgId();
0332             simInfo[i].motherStatus = genMom->status();
0333             simInfo[i].motherVertex = genMom->vertex();
0334           }
0335           dumpFormatedInfo(simInfo[i]);
0336           LogTrace("MuonSimClassifier") << "\t   has GEN mother pdgId = " << simInfo[i].motherPdgId
0337                                         << " (status = " << simInfo[i].motherStatus << ")";
0338 
0339           reco::GenParticleRef genGMom = genMom->numberOfMothers() > 0 ? genMom->motherRef() : reco::GenParticleRef();
0340 
0341           if (genGMom.isNonnull()) {
0342             simInfo[i].grandMotherPdgId = genGMom->pdgId();
0343             LogTrace("MuonSimClassifier")
0344                 << "\t\t mother prod. vertex rho = " << simInfo[i].motherVertex.Rho()
0345                 << ", z = " << simInfo[i].motherVertex.Z() << ", grand-mom pdgId = " << simInfo[i].grandMotherPdgId;
0346           }
0347           // in this case, we might want to know the heaviest mom flavour
0348           for (reco::GenParticleRef nMom = genMom;
0349                nMom.isNonnull() && std::abs(nMom->pdgId()) >= 100;  // stop when we're no longer
0350                                                                     // looking at hadrons or mesons
0351                nMom = nMom->numberOfMothers() > 0 ? nMom->motherRef() : reco::GenParticleRef()) {
0352             int flav = flavour(nMom->pdgId());
0353             if (simInfo[i].heaviestMotherFlavour < flav)
0354               simInfo[i].heaviestMotherFlavour = flav;
0355             LogTrace("MuonSimClassifier")
0356                 << "\t\t backtracking flavour: mom pdgId = " << nMom->pdgId() << ", flavour = " << flav
0357                 << ", heaviest so far = " << simInfo[i].heaviestMotherFlavour;
0358           }
0359         } else {  // mother is null ??
0360           dumpFormatedInfo(simInfo[i]);
0361           LogTrace("MuonSimClassifier") << "\t   has NO mother!";
0362         }
0363       } else {  // Muon is in SIM Only
0364         TrackingParticleRef simMom = getTpMother(tp);
0365         if (simMom.isNonnull()) {
0366           simInfo[i].motherPdgId = simMom->pdgId();
0367           simInfo[i].motherVertex = simMom->vertex();
0368           dumpFormatedInfo(simInfo[i]);
0369           LogTrace("MuonSimClassifier") << "\t   has SIM mother pdgId = " << simInfo[i].motherPdgId
0370                                         << " produced at rho = " << simMom->vertex().Rho()
0371                                         << ", z = " << simMom->vertex().Z();
0372 
0373           if (!simMom->genParticles().empty()) {
0374             simInfo[i].motherStatus = simMom->genParticles()[0]->status();
0375             reco::GenParticleRef genGMom =
0376                 (simMom->genParticles()[0]->numberOfMothers() > 0 ? simMom->genParticles()[0]->motherRef()
0377                                                                   : reco::GenParticleRef());
0378             if (genGMom.isNonnull())
0379               simInfo[i].grandMotherPdgId = genGMom->pdgId();
0380             LogTrace("MuonSimClassifier") << "\t\t SIM mother is in GEN (status " << simInfo[i].motherStatus
0381                                           << "), grand-mom id = " << simInfo[i].grandMotherPdgId;
0382           } else {
0383             simInfo[i].motherStatus = -1;
0384             TrackingParticleRef simGMom = getTpMother(simMom);
0385             if (simGMom.isNonnull())
0386               simInfo[i].grandMotherPdgId = simGMom->pdgId();
0387             LogTrace("MuonSimClassifier")
0388                 << "\t\t SIM mother is in SIM only, grand-mom id = " << simInfo[i].grandMotherPdgId;
0389           }
0390         } else {
0391           dumpFormatedInfo(simInfo[i]);
0392           LogTrace("MuonSimClassifier") << "\t   has NO mother!";
0393         }
0394       }
0395       simInfo[i].motherFlavour = flavour(simInfo[i].motherPdgId);
0396       simInfo[i].grandMotherFlavour = flavour(simInfo[i].grandMotherPdgId);
0397 
0398       // Check first IF this is a muon at all
0399       if (std::abs(tp->pdgId()) != 13) {
0400         if (std::abs(tp->pdgId()) == 11) {
0401           simInfo[i].primaryClass = isGhost ? reco::MuonSimType::GhostElectron : reco::MuonSimType::MatchedElectron;
0402           simInfo[i].extendedClass =
0403               isGhost ? reco::ExtendedMuonSimType::ExtGhostElectron : reco::ExtendedMuonSimType::ExtMatchedElectron;
0404           LogTrace("MuonSimClassifier") << "\t This is electron/positron. classif[i] = " << simInfo[i].primaryClass;
0405         } else {
0406           simInfo[i].primaryClass =
0407               isGhost ? reco::MuonSimType::GhostPunchthrough : reco::MuonSimType::MatchedPunchthrough;
0408           simInfo[i].extendedClass = isGhost ? reco::ExtendedMuonSimType::ExtGhostPunchthrough
0409                                              : reco::ExtendedMuonSimType::ExtMatchedPunchthrough;
0410           LogTrace("MuonSimClassifier") << "\t This is not a muon. Sorry. classif[i] = " << simInfo[i].primaryClass;
0411         }
0412         continue;
0413       }
0414 
0415       // Is this SIM muon also a GEN muon, with a mother?
0416       if (!tp->genParticles().empty() && (simInfo[i].motherPdgId != 0)) {
0417         if (std::abs(simInfo[i].motherPdgId) < 100 && (std::abs(simInfo[i].motherPdgId) != 15)) {
0418           simInfo[i].primaryClass =
0419               isGhost ? reco::MuonSimType::GhostPrimaryMuon : reco::MuonSimType::MatchedPrimaryMuon;
0420           simInfo[i].flavour = 13;
0421           simInfo[i].extendedClass = isGhost ? reco::ExtendedMuonSimType::GhostMuonFromGaugeOrHiggsBoson
0422                                              : reco::ExtendedMuonSimType::MatchedMuonFromGaugeOrHiggsBoson;
0423           LogTrace("MuonSimClassifier") << "\t This seems PRIMARY MUON ! classif[i] = " << simInfo[i].primaryClass;
0424         } else if (simInfo[i].motherFlavour == 4 || simInfo[i].motherFlavour == 5 || simInfo[i].motherFlavour == 15) {
0425           simInfo[i].primaryClass =
0426               isGhost ? reco::MuonSimType::GhostMuonFromHeavyFlavour : reco::MuonSimType::MatchedMuonFromHeavyFlavour;
0427           simInfo[i].flavour = simInfo[i].motherFlavour;
0428           if (simInfo[i].motherFlavour == 15)
0429             simInfo[i].extendedClass =
0430                 isGhost ? reco::ExtendedMuonSimType::GhostMuonFromTau : reco::ExtendedMuonSimType::MatchedMuonFromTau;
0431           else if (simInfo[i].motherFlavour == 5)
0432             simInfo[i].extendedClass =
0433                 isGhost ? reco::ExtendedMuonSimType::GhostMuonFromB : reco::ExtendedMuonSimType::MatchedMuonFromB;
0434           else if (simInfo[i].heaviestMotherFlavour == 5)
0435             simInfo[i].extendedClass =
0436                 isGhost ? reco::ExtendedMuonSimType::GhostMuonFromBtoC : reco::ExtendedMuonSimType::MatchedMuonFromBtoC;
0437           else
0438             simInfo[i].extendedClass =
0439                 isGhost ? reco::ExtendedMuonSimType::GhostMuonFromC : reco::ExtendedMuonSimType::MatchedMuonFromC;
0440           LogTrace("MuonSimClassifier") << "\t This seems HEAVY FLAVOUR ! classif[i] = " << simInfo[i].primaryClass;
0441         } else {
0442           simInfo[i].primaryClass =
0443               isGhost ? reco::MuonSimType::GhostMuonFromLightFlavour : reco::MuonSimType::MatchedMuonFromLightFlavour;
0444           simInfo[i].flavour = simInfo[i].motherFlavour;
0445           LogTrace("MuonSimClassifier") << "\t This seems LIGHT FLAVOUR ! classif[i] = " << simInfo[i].primaryClass;
0446         }
0447       } else {
0448         simInfo[i].primaryClass =
0449             isGhost ? reco::MuonSimType::GhostMuonFromLightFlavour : reco::MuonSimType::MatchedMuonFromLightFlavour;
0450         simInfo[i].flavour = simInfo[i].motherFlavour;
0451         LogTrace("MuonSimClassifier") << "\t This seems LIGHT FLAVOUR ! classif[i] = " << simInfo[i].primaryClass;
0452       }
0453 
0454       // extended classification
0455       // don't we override previous decisions?
0456       if (simInfo[i].motherPdgId == 0)
0457         // if it has no mom, it's not a primary particle so it won't be in ppMuX
0458         simInfo[i].extendedClass = isGhost ? reco::ExtendedMuonSimType::GhostMuonFromNonPrimaryParticle
0459                                            : reco::ExtendedMuonSimType::MatchedMuonFromNonPrimaryParticle;
0460       else if (std::abs(simInfo[i].motherPdgId) < 100) {
0461         if (simInfo[i].motherFlavour == 15)
0462           simInfo[i].extendedClass =
0463               isGhost ? reco::ExtendedMuonSimType::GhostMuonFromTau : reco::ExtendedMuonSimType::MatchedMuonFromTau;
0464         else
0465           simInfo[i].extendedClass = isGhost ? reco::ExtendedMuonSimType::GhostMuonFromGaugeOrHiggsBoson
0466                                              : reco::ExtendedMuonSimType::MatchedMuonFromGaugeOrHiggsBoson;
0467       } else if (simInfo[i].motherFlavour == 5)
0468         simInfo[i].extendedClass =
0469             isGhost ? reco::ExtendedMuonSimType::GhostMuonFromB : reco::ExtendedMuonSimType::MatchedMuonFromB;
0470       else if (simInfo[i].motherFlavour == 4) {
0471         if (simInfo[i].heaviestMotherFlavour == 5)
0472           simInfo[i].extendedClass =
0473               isGhost ? reco::ExtendedMuonSimType::GhostMuonFromBtoC : reco::ExtendedMuonSimType::MatchedMuonFromBtoC;
0474         else
0475           simInfo[i].extendedClass =
0476               isGhost ? reco::ExtendedMuonSimType::GhostMuonFromC : reco::ExtendedMuonSimType::MatchedMuonFromC;
0477       } else if (simInfo[i].motherStatus != -1) {  // primary light particle
0478         int id = std::abs(simInfo[i].motherPdgId);
0479         if (id != /*pi+*/ 211 && id != /*K+*/ 321 && id != 130 /*K0L*/)
0480           // other light particle, possibly short-lived
0481           simInfo[i].extendedClass = isGhost ? reco::ExtendedMuonSimType::GhostMuonFromOtherLight
0482                                              : reco::ExtendedMuonSimType::MatchedMuonFromOtherLight;
0483         else if (simInfo[i].vertex.Rho() < decayRho_ && std::abs(simInfo[i].vertex.Z()) < decayAbsZ_)
0484           // decay a la ppMuX (primary pi/K within a cylinder)
0485           simInfo[i].extendedClass = isGhost ? reco::ExtendedMuonSimType::GhostMuonFromPiKppMuX
0486                                              : reco::ExtendedMuonSimType::MatchedMuonFromPiKppMuX;
0487         else
0488           // late decay that wouldn't be in ppMuX
0489           simInfo[i].extendedClass = isGhost ? reco::ExtendedMuonSimType::GhostMuonFromPiKNotppMuX
0490                                              : reco::ExtendedMuonSimType::MatchedMuonFromPiKNotppMuX;
0491       } else
0492         // decay of non-primary particle, would not be in ppMuX
0493         simInfo[i].extendedClass = isGhost ? reco::ExtendedMuonSimType::GhostMuonFromNonPrimaryParticle
0494                                            : reco::ExtendedMuonSimType::MatchedMuonFromNonPrimaryParticle;
0495 
0496       if (linkToGenParticles_ && std::abs(simInfo[i].extendedClass) >= 2) {
0497         // Link to the genParticle if possible, but not decays in flight (in
0498         // ppMuX they're in GEN block, but they have wrong parameters)
0499         if (!tp->genParticles().empty() && std::abs(simInfo[i].extendedClass) >= 5) {
0500           if (genParticles.id() != tp->genParticles().id()) {
0501             throw cms::Exception("Configuration")
0502                 << "Product ID mismatch between the genParticle collection (" << genParticles_ << ", id "
0503                 << genParticles.id() << ") and the references in the TrackingParticles (id " << tp->genParticles().id()
0504                 << ")\n";
0505           }
0506           muToPrimary[i] = tp->genParticles()[0].key();
0507         } else {
0508           // Don't put the same trackingParticle twice!
0509           int &indexPlus1 = tpToSecondaries[tp];  // will create a 0 if the tp is
0510                                                   // not in the list already
0511           if (indexPlus1 == 0)
0512             indexPlus1 = convertAndPush(*tp, *secondaries, getTpMother(tp), genParticles) + 1;
0513           muToSecondary[i] = indexPlus1 - 1;
0514         }
0515       }
0516       LogTrace("MuonSimClassifier") << "\t Extended classification code = " << simInfo[i].extendedClass;
0517     } else {  // if (tp.isNonnull())
0518       simInfo[i].primaryClass = reco::MuonSimType::NotMatched;
0519       simInfo[i].extendedClass = reco::ExtendedMuonSimType::ExtNotMatched;
0520     }
0521   }  // end loop on reco muons
0522 
0523   writeValueMap(iEvent, muons, simInfo, "");
0524 
0525   if (linkToGenParticles_) {
0526     edm::OrphanHandle<reco::GenParticleCollection> secHandle = iEvent.put(std::move(secondaries), "secondaries");
0527     edm::RefProd<reco::GenParticleCollection> priRP(genParticles);
0528     edm::RefProd<reco::GenParticleCollection> secRP(secHandle);
0529     std::unique_ptr<edm::Association<reco::GenParticleCollection>> outPri(
0530         new edm::Association<reco::GenParticleCollection>(priRP));
0531     std::unique_ptr<edm::Association<reco::GenParticleCollection>> outSec(
0532         new edm::Association<reco::GenParticleCollection>(secRP));
0533     edm::Association<reco::GenParticleCollection>::Filler fillPri(*outPri), fillSec(*outSec);
0534     fillPri.insert(muons, muToPrimary.begin(), muToPrimary.end());
0535     fillSec.insert(muons, muToSecondary.begin(), muToSecondary.end());
0536     fillPri.fill();
0537     fillSec.fill();
0538     iEvent.put(std::move(outPri), "toPrimaries");
0539     iEvent.put(std::move(outSec), "toSecondaries");
0540   }
0541 }
0542 
0543 template <typename T>
0544 void MuonSimClassifier::writeValueMap(edm::Event &iEvent,
0545                                       const edm::Handle<edm::View<reco::Muon>> &handle,
0546                                       const std::vector<T> &values,
0547                                       const std::string &label) const {
0548   using namespace edm;
0549   using namespace std;
0550   unique_ptr<ValueMap<T>> valMap(new ValueMap<T>());
0551   typename edm::ValueMap<T>::Filler filler(*valMap);
0552   filler.insert(handle, values.begin(), values.end());
0553   filler.fill();
0554   iEvent.put(std::move(valMap), label);
0555 }
0556 
0557 int MuonSimClassifier::flavour(int pdgId) const {
0558   int flav = std::abs(pdgId);
0559   // for quarks, leptons and bosons except gluons, take their pdgId
0560   // muons and taus have themselves as flavour
0561   if (flav <= 37 && flav != 21)
0562     return flav;
0563   // look for barions
0564   int bflav = ((flav / 1000) % 10);
0565   if (bflav != 0)
0566     return bflav;
0567   // look for mesons
0568   int mflav = ((flav / 100) % 10);
0569   if (mflav != 0)
0570     return mflav;
0571   return 0;
0572 }
0573 
0574 // push secondary in collection.
0575 // if it has a primary mother link to it.
0576 int MuonSimClassifier::convertAndPush(const TrackingParticle &tp,
0577                                       reco::GenParticleCollection &out,
0578                                       const TrackingParticleRef &simMom,
0579                                       const edm::Handle<reco::GenParticleCollection> &genParticles) const {
0580   out.push_back(reco::GenParticle(tp.charge(), tp.p4(), tp.vertex(), tp.pdgId(), tp.status(), true));
0581   if (simMom.isNonnull() && !simMom->genParticles().empty()) {
0582     if (genParticles.id() != simMom->genParticles().id()) {
0583       throw cms::Exception("Configuration")
0584           << "Product ID mismatch between the genParticle collection (" << genParticles_ << ", id " << genParticles.id()
0585           << ") and the references in the TrackingParticles (id " << simMom->genParticles().id() << ")\n";
0586     }
0587     out.back().addMother(simMom->genParticles()[0]);
0588   }
0589   return out.size() - 1;
0590 }
0591 
0592 // define this as a plug-in
0593 DEFINE_FWK_MODULE(MuonSimClassifier);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!