Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:22:47 

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