Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    MuonMCClassifier
00004 // Class:      MuonMCClassifier
00005 // 
00028 //
00029 // Original Author:  Nov 16 16:12 (lxplus231.cern.ch)
00030 //         Created:  Sun Nov 16 16:14:09 CET 2008
00031 // $Id: MuonMCClassifier.cc,v 1.7 2011/01/24 11:24:32 gpetrucc Exp $
00032 //
00033 //
00034 
00035 
00036 // system include files
00037 #include <memory>
00038 #include <set>
00039 
00040 // user include files
00041 #include "FWCore/Framework/interface/Frameworkfwd.h"
00042 #include "FWCore/Framework/interface/EDProducer.h"
00043 
00044 #include "FWCore/Framework/interface/Event.h"
00045 #include "FWCore/Framework/interface/MakerMacros.h"
00046 #include "FWCore/Framework/interface/EventSetup.h"
00047 #include "FWCore/Framework/interface/ESHandle.h"
00048 
00049 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00050 
00051 #include "DataFormats/Common/interface/View.h"
00052 #include "DataFormats/Common/interface/ValueMap.h"
00053 #include "DataFormats/Common/interface/Association.h"
00054 
00055 #include "DataFormats/MuonReco/interface/Muon.h"
00056 #include "DataFormats/PatCandidates/interface/Muon.h"
00057 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
00058 
00059 #include "SimMuon/MCTruth/interface/MuonAssociatorByHits.h"
00060 #include "SimTracker/Records/interface/TrackAssociatorRecord.h"
00061 #include <SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h>
00062 
00063 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
00064 
00065 #include <boost/foreach.hpp>
00066 #define foreach BOOST_FOREACH
00067 
00068 //
00069 // class decleration
00070 class MuonMCClassifier : public edm::EDProducer {
00071     public:
00072         explicit MuonMCClassifier(const edm::ParameterSet&);
00073         ~MuonMCClassifier();
00074 
00075     private:
00076         virtual void produce(edm::Event&, const edm::EventSetup&);
00078         edm::InputTag muons_;
00079 
00082         bool hasMuonCut_;
00083         StringCutObjectSelector<pat::Muon> muonCut_;
00084  
00086         MuonAssociatorByHits::MuonTrackType trackType_;
00087 
00089         edm::InputTag trackingParticles_;
00090 
00092         std::string associatorLabel_;
00093 
00095         double decayRho_, decayAbsZ_;
00096 
00098         bool linkToGenParticles_; 
00099         edm::InputTag genParticles_; 
00100 
00102         int flavour(int pdgId) const ;
00103 
00105         template<typename T>
00106         void writeValueMap(edm::Event &iEvent,
00107                 const edm::Handle<edm::View<reco::Muon> > & handle,
00108                 const std::vector<T> & values,
00109                 const std::string    & label) const ;
00110 
00111         TrackingParticleRef getTpMother(TrackingParticleRef tp) {
00112             if (tp.isNonnull() && tp->parentVertex().isNonnull() && !tp->parentVertex()->sourceTracks().empty()) {
00113                 return tp->parentVertex()->sourceTracks()[0];
00114             } else {
00115                 return TrackingParticleRef();
00116             }
00117         }
00118 
00119         const HepMC::GenParticle * getGpMother(const HepMC::GenParticle *gp) {
00120             if (gp != 0) {
00121                 const HepMC::GenVertex *vtx = gp->production_vertex();
00122                 if (vtx != 0 && vtx->particles_in_size() > 0) {
00123                     return *vtx->particles_in_const_begin();
00124                 }
00125             }
00126             return 0;
00127         }
00128 
00130         int fetch(const edm::Handle<std::vector<int> > & genBarcodes, int barcode) const;
00131 
00135         int convertAndPush(const TrackingParticle &tp, 
00136                            reco::GenParticleCollection &out, 
00137                            const TrackingParticleRef &momRef, 
00138                            const edm::Handle<reco::GenParticleCollection> & genParticles, 
00139                            const edm::Handle<std::vector<int> > & genBarcodes) const ;
00140 };
00141 
00142 MuonMCClassifier::MuonMCClassifier(const edm::ParameterSet &iConfig) :
00143     muons_(iConfig.getParameter<edm::InputTag>("muons")),
00144     hasMuonCut_(iConfig.existsAs<std::string>("muonPreselection")),
00145     muonCut_(hasMuonCut_ ? iConfig.getParameter<std::string>("muonPreselection") : ""),
00146     trackingParticles_(iConfig.getParameter<edm::InputTag>("trackingParticles")),
00147     associatorLabel_(iConfig.getParameter< std::string >("associatorLabel")),
00148     decayRho_(iConfig.getParameter<double>("decayRho")),
00149     decayAbsZ_(iConfig.getParameter<double>("decayAbsZ")),
00150     linkToGenParticles_(iConfig.getParameter<bool>("linkToGenParticles")),
00151     genParticles_(linkToGenParticles_ ? iConfig.getParameter<edm::InputTag>("genParticles") : edm::InputTag("NONE"))
00152 {
00153     std::string trackType = iConfig.getParameter< std::string >("trackType");
00154     if (trackType == "inner") trackType_ = MuonAssociatorByHits::InnerTk;
00155     else if (trackType == "outer") trackType_ = MuonAssociatorByHits::OuterTk;
00156     else if (trackType == "global") trackType_ = MuonAssociatorByHits::GlobalTk;
00157     else if (trackType == "segments") trackType_ = MuonAssociatorByHits::Segments;
00158     else throw cms::Exception("Configuration") << "Track type '" << trackType << "' not supported.\n";
00159 
00160     produces<edm::ValueMap<int> >(); 
00161     produces<edm::ValueMap<int> >("ext"); 
00162     produces<edm::ValueMap<int> >("flav"); 
00163     produces<edm::ValueMap<int> >("hitsPdgId"); 
00164     produces<edm::ValueMap<int> >("momPdgId"); 
00165     produces<edm::ValueMap<int> >("momFlav"); 
00166     produces<edm::ValueMap<int> >("momStatus"); 
00167     produces<edm::ValueMap<int> >("gmomPdgId"); 
00168     produces<edm::ValueMap<int> >("gmomFlav"); 
00169     produces<edm::ValueMap<int> >("hmomFlav"); // heaviest mother flavour
00170     produces<edm::ValueMap<int> >("tpId");
00171     produces<edm::ValueMap<float> >("prodRho"); 
00172     produces<edm::ValueMap<float> >("prodZ"); 
00173     produces<edm::ValueMap<float> >("momRho"); 
00174     produces<edm::ValueMap<float> >("momZ"); 
00175     produces<edm::ValueMap<float> >("tpAssoQuality");
00176     if (linkToGenParticles_) {
00177         produces<reco::GenParticleCollection>("secondaries");
00178         produces<edm::Association<reco::GenParticleCollection> >("toPrimaries");
00179         produces<edm::Association<reco::GenParticleCollection> >("toSecondaries");
00180     }
00181 }
00182 
00183 MuonMCClassifier::~MuonMCClassifier() 
00184 {
00185 }
00186 
00187 void
00188 MuonMCClassifier::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00189 {
00190     edm::LogVerbatim("MuonMCClassifier") <<"\n sono in MuonMCClassifier !";
00191   
00192     edm::Handle<edm::View<reco::Muon> > muons; 
00193     iEvent.getByLabel(muons_, muons);
00194 
00195     edm::Handle<TrackingParticleCollection> trackingParticles;
00196     iEvent.getByLabel(trackingParticles_, trackingParticles);
00197 
00198     edm::Handle<reco::GenParticleCollection> genParticles;
00199     edm::Handle<std::vector<int> > genBarcodes;
00200     if (linkToGenParticles_) {
00201         iEvent.getByLabel(genParticles_, genParticles);
00202         iEvent.getByLabel(genParticles_, genBarcodes);
00203     }
00204 
00205     edm::ESHandle<TrackAssociatorBase> associatorBase;
00206     iSetup.get<TrackAssociatorRecord>().get(associatorLabel_, associatorBase);
00207     const MuonAssociatorByHits * assoByHits = dynamic_cast<const MuonAssociatorByHits *>(associatorBase.product());
00208     if (assoByHits == 0) throw cms::Exception("Configuration") << "The Track Associator with label '" << associatorLabel_ << "' is not a MuonAssociatorByHits.\n";
00209 
00210     MuonAssociatorByHits::MuonToSimCollection recSimColl;
00211     MuonAssociatorByHits::SimToMuonCollection simRecColl;
00212     edm::LogVerbatim("MuonMCClassifier") <<"\n ***************************************************************** ";
00213     edm::LogVerbatim("MuonMCClassifier") <<  " RECO MUON association, type:  "<< trackType_;
00214     edm::LogVerbatim("MuonMCClassifier") <<  " ***************************************************************** \n";
00215 
00216     edm::RefToBaseVector<reco::Muon> selMuons;
00217     if (!hasMuonCut_) {
00218         // all muons
00219         selMuons = muons->refVector();
00220     } else {
00221         // filter, fill refvectors, associate
00222         // I pass through pat::Muon so that I can access muon id selectors
00223         for (size_t i = 0, n = muons->size(); i < n; ++i) {
00224             edm::RefToBase<reco::Muon> rmu = muons->refAt(i);
00225             if (muonCut_(pat::Muon(rmu))) selMuons.push_back(rmu);
00226         }
00227     }
00228 
00229     edm::RefVector<TrackingParticleCollection> allTPs;
00230     for (size_t i = 0, n = trackingParticles->size(); i < n; ++i) {
00231         allTPs.push_back(TrackingParticleRef(trackingParticles,i));
00232     }
00233 
00234     assoByHits->associateMuons(recSimColl, simRecColl, selMuons, trackType_, allTPs, &iEvent, &iSetup);
00235 
00236     // for global muons without hits on muon detectors, look at the linked standalone muon
00237     MuonAssociatorByHits::MuonToSimCollection UpdSTA_recSimColl;
00238     MuonAssociatorByHits::SimToMuonCollection UpdSTA_simRecColl;
00239     if (trackType_ == MuonAssociatorByHits::GlobalTk) {
00240       edm::LogVerbatim("MuonMCClassifier") <<"\n ***************************************************************** ";
00241       edm::LogVerbatim("MuonMCClassifier") <<  " STANDALONE (UpdAtVtx) MUON association ";
00242       edm::LogVerbatim("MuonMCClassifier") <<  " ***************************************************************** \n";
00243       assoByHits->associateMuons(UpdSTA_recSimColl, UpdSTA_simRecColl, selMuons, MuonAssociatorByHits::OuterTk, 
00244                                  allTPs, &iEvent, &iSetup);
00245     }
00246 
00247     typedef MuonAssociatorByHits::MuonToSimCollection::const_iterator r2s_it;
00248     typedef MuonAssociatorByHits::SimToMuonCollection::const_iterator s2r_it;
00249 
00250     size_t nmu = muons->size();
00251     edm::LogVerbatim("MuonMCClassifier") <<"\n There are "<<nmu<<" reco::Muons.";
00252 
00253     std::vector<int> classif(nmu, 0), ext(nmu, 0);
00254     std::vector<int> hitsPdgId(nmu, 0), momPdgId(nmu, 0), gmomPdgId(nmu, 0), momStatus(nmu, 0);
00255     std::vector<int> flav(nmu, 0),      momFlav(nmu, 0),  gmomFlav(nmu, 0), hmomFlav(nmu, 0);
00256     std::vector<int> tpId(nmu, -1);
00257     std::vector<float> prodRho(nmu, 0.0), prodZ(nmu, 0.0), momRho(nmu, 0.0), momZ(nmu, 0.0);
00258     std::vector<float> tpAssoQuality(nmu, -1);
00259 
00260     std::auto_ptr<reco::GenParticleCollection> secondaries;     // output collection of secondary muons
00261     std::map<TrackingParticleRef, int>         tpToSecondaries; // map from tp to (index+1) in output collection
00262     std::vector<int> muToPrimary(nmu, -1), muToSecondary(nmu, -1); // map from input into (index) in output, -1 for null
00263     if (linkToGenParticles_) secondaries.reset(new reco::GenParticleCollection());
00264 
00265     for(size_t i = 0; i < nmu; ++i) {
00266         edm::LogVerbatim("MuonMCClassifier") <<"\n reco::Muons # "<<i;
00267         edm::RefToBase<reco::Muon> mu = muons->refAt(i);
00268         if (hasMuonCut_ && (std::find(selMuons.begin(), selMuons.end(), mu) == selMuons.end()) ) {
00269             edm::LogVerbatim("MuonMCClassifier") <<"\t muon didn't pass the selection. classified as -99 and skipped";
00270             classif[i] = -99; continue;
00271         }
00272 
00273         TrackingParticleRef        tp;
00274         edm::RefToBase<reco::Muon> muMatchBack;
00275         r2s_it match = recSimColl.find(mu);
00276         s2r_it matchback;
00277         if (match != recSimColl.end()) {
00278             edm::LogVerbatim("MuonMCClassifier") <<"\t RtS matched Ok...";
00279             // match->second is vector, front is first element, first is the ref (second would be the quality)
00280             tp = match->second.front().first;
00281             tpId[i]          = tp.isNonnull() ? tp.key() : -1; // we check, even if null refs should not appear here at all
00282             tpAssoQuality[i] = match->second.front().second;
00283             s2r_it matchback = simRecColl.find(tp);
00284             if (matchback != simRecColl.end()) {
00285                 muMatchBack = matchback->second.front().first;
00286             } else {
00287                 edm::LogWarning("MuonMCClassifier") << "\n***WARNING:  This I do NOT understand: why no match back? *** \n";
00288             }
00289         } else if ((trackType_ == MuonAssociatorByHits::GlobalTk) &&
00290                     mu->isGlobalMuon()) {
00291             // perform a second attempt, matching with the standalone muon
00292             r2s_it matchSta = UpdSTA_recSimColl.find(mu);
00293             if (matchSta != UpdSTA_recSimColl.end()) {
00294                 edm::LogVerbatim("MuonMCClassifier") <<"\t RtS matched Ok... from the UpdSTA_recSimColl ";
00295                 tp    = matchSta->second.front().first;
00296                 tpId[i]          = tp.isNonnull() ? tp.key() : -1; // we check, even if null refs should not appear here at all
00297                 tpAssoQuality[i] = matchSta->second.front().second;
00298                 s2r_it matchback = UpdSTA_simRecColl.find(tp);
00299                 if (matchback != UpdSTA_simRecColl.end()) {
00300                     muMatchBack = matchback->second.front().first;
00301                 } else {
00302                     edm::LogWarning("MuonMCClassifier") << "\n***WARNING:  This I do NOT understand: why no match back in UpdSTA? *** \n";
00303                 }
00304             }
00305         } 
00306         if (tp.isNonnull()) {
00307             bool isGhost = muMatchBack != mu;
00308             if (isGhost) edm::LogVerbatim("MuonMCClassifier") <<"\t This seems a GHOST ! classif[i] will be < 0";
00309 
00310             hitsPdgId[i] = tp->pdgId();
00311             prodRho[i]   = tp->vertex().Rho(); 
00312             prodZ[i]     = tp->vertex().Z();
00313             edm::LogVerbatim("MuonMCClassifier") <<"\t TP pdgId = "<<hitsPdgId[i] << ", vertex rho = " << prodRho[i] << ", z = " << prodZ[i];
00314 
00315             // Try to extract mother and grand mother of this muon.
00316             // Unfortunately, SIM and GEN histories require diffent code :-(
00317             if (!tp->genParticle().empty()) { // Muon is in GEN
00318                 const HepMC::GenParticle * genMom = getGpMother(tp->genParticle()[0].get());
00319                 if (genMom) {
00320                     momPdgId[i]  = genMom->pdg_id();
00321                     momStatus[i] = genMom->status();
00322                     if (genMom->production_vertex()) {
00323                         const HepMC::ThreeVector & momVtx = genMom->production_vertex()->point3d();
00324                         momRho[i] = momVtx.perp() * 0.1; momZ[i] = momVtx.z() * 0.1; // HepMC is in mm!
00325                     }
00326                     edm::LogVerbatim("MuonMCClassifier") << "\t Particle pdgId = "<<hitsPdgId[i] << " produced at rho = " << prodRho[i] << ", z = " << prodZ[i] << ", has GEN mother pdgId = " << momPdgId[i];
00327                     const HepMC::GenParticle * genGMom = getGpMother(genMom);
00328                     if (genGMom) {
00329                         gmomPdgId[i] = genGMom->pdg_id();
00330                         edm::LogVerbatim("MuonMCClassifier") << "\t\t mother prod. vertex rho = " << momRho[i] << ", z = " << momZ[i] << ", grand-mom pdgId = " << gmomPdgId[i];
00331                     }
00332                     // in this case, we might want to know the heaviest mom flavour
00333                     for (const HepMC::GenParticle *nMom = genMom; 
00334                             nMom != 0 && abs(nMom->pdg_id()) >= 100; // stop when we're no longer looking at hadrons or mesons
00335                             nMom = getGpMother(nMom)) {
00336                         int flav = flavour(nMom->pdg_id());
00337                         if (hmomFlav[i] < flav) hmomFlav[i] = flav; 
00338                         edm::LogVerbatim("MuonMCClassifier") << "\t\t backtracking flavour: mom pdgId = "<<nMom->pdg_id()<< ", flavour = " << flav << ", heaviest so far = " << hmomFlav[i];
00339                     }
00340                 }
00341             } else { // Muon is in SIM Only
00342                 TrackingParticleRef simMom = getTpMother(tp);
00343                 if (simMom.isNonnull()) {
00344                     momPdgId[i] = simMom->pdgId();
00345                     momRho[i] = simMom->vertex().Rho();
00346                     momZ[i]   = simMom->vertex().Z();
00347                     edm::LogVerbatim("MuonMCClassifier") << "\t Particle pdgId = "<<hitsPdgId[i] << " produced at rho = " << prodRho[i] << ", z = " << prodZ[i] << 
00348                                                             ", has SIM mother pdgId = " << momPdgId[i] << " produced at rho = " << simMom->vertex().Rho() << ", z = " << simMom->vertex().Z();
00349                     if (!simMom->genParticle().empty()) {
00350                         momStatus[i] = simMom->genParticle()[0]->status();
00351                         const HepMC::GenParticle * genGMom = getGpMother(simMom->genParticle()[0].get());
00352                         if (genGMom) gmomPdgId[i] = genGMom->pdg_id();
00353                         edm::LogVerbatim("MuonMCClassifier") << "\t\t SIM mother is in GEN (status " << momStatus[i] << "), grand-mom id = " << gmomPdgId[i];
00354                     } else {
00355                         momStatus[i] = -1;
00356                         TrackingParticleRef simGMom = getTpMother(simMom);
00357                         if (simGMom.isNonnull()) gmomPdgId[i] = simGMom->pdgId();
00358                         edm::LogVerbatim("MuonMCClassifier") << "\t\t SIM mother is in SIM only, grand-mom id = " << gmomPdgId[i];
00359                     }
00360                 } else {
00361                   edm::LogVerbatim("MuonMCClassifier") << "\t Particle pdgId = "<<hitsPdgId[i] << " produced at rho = " << prodRho[i] << ", z = " << prodZ[i] << ", has NO mother!";
00362                 }
00363             }
00364             momFlav[i]  = flavour(momPdgId[i]);
00365             gmomFlav[i] = flavour(gmomPdgId[i]);
00366 
00367             // Check first IF this is a muon at all
00368             if (abs(tp->pdgId()) != 13) {
00369                 classif[i] = isGhost ? -1 : 1;
00370                 ext[i]     = isGhost ? -1 : 1;
00371                 edm::LogVerbatim("MuonMCClassifier") <<"\t This is not a muon. Sorry. classif[i] = " << classif[i];
00372                 continue;
00373             }
00374 
00375             // Is this SIM muon also a GEN muon, with a mother?
00376             if (!tp->genParticle().empty() && (momPdgId[i] != 0)) {
00377                 if (abs(momPdgId[i]) < 100 && (abs(momPdgId[i]) != 15)) {
00378                     classif[i] = isGhost ? -4 : 4;
00379                     flav[i] = (abs(momPdgId[i]) == 15 ? 15 : 13);
00380                     edm::LogVerbatim("MuonMCClassifier") <<"\t This seems PRIMARY MUON ! classif[i] = " << classif[i];
00381                     ext[i] = 10;
00382                 } else if (momFlav[i] == 4 || momFlav[i] == 5 || momFlav[i] == 15) {
00383                     classif[i] = isGhost ? -3 : 3;
00384                     flav[i]    = momFlav[i];
00385                     if (momFlav[i] == 15)      ext[i] = 9; // tau->mu
00386                     else if (momFlav[i] == 5)  ext[i] = 8; // b->mu
00387                     else if (hmomFlav[i] == 5) ext[i] = 7; // b->c->mu
00388                     else                       ext[i] = 6; // c->mu
00389                     edm::LogVerbatim("MuonMCClassifier") <<"\t This seems HEAVY FLAVOUR ! classif[i] = " << classif[i];
00390                 } else {
00391                     classif[i] = isGhost ? -2 : 2;
00392                     flav[i]    = momFlav[i];
00393                     edm::LogVerbatim("MuonMCClassifier") <<"\t This seems LIGHT FLAVOUR ! classif[i] = " << classif[i];
00394                 }
00395             } else {
00396                 classif[i] = isGhost ? -2 : 2;
00397                 flav[i]    = momFlav[i];
00398                 edm::LogVerbatim("MuonMCClassifier") <<"\t This seems LIGHT FLAVOUR ! classif[i] = " << classif[i];
00399             }
00400             // extended classification
00401             if (momPdgId[i] == 0) ext[i] = 2; // if it has no mom, it's not a primary particle so it won't be in ppMuX
00402             else if (abs(momPdgId[i]) < 100) ext[i] = (momFlav[i] == 15 ? 9 : 10); // primary mu, or tau->mu
00403             else if (momFlav[i] == 5) ext[i] = 8; // b->mu
00404             else if (momFlav[i] == 4) ext[i] = (hmomFlav[i] == 5 ? 7 : 6); // b->c->mu and c->mu
00405             else if (momStatus[i] != -1) { // primary light particle
00406                 int id = abs(momPdgId[i]);
00407                 if (id != /*pi+*/211 && id != /*K+*/321 && id != 130 /*K0L*/)  ext[i] = 5; // other light particle, possibly short-lived
00408                 else if (prodRho[i] < decayRho_ && abs(prodZ[i]) < decayAbsZ_) ext[i] = 4; // decay a la ppMuX (primary pi/K within a cylinder)
00409                 else                                                           ext[i] = 3; // late decay that wouldn't be in ppMuX
00410             } else ext[i] = 2; // decay of non-primary particle, would not be in ppMuX
00411             if (isGhost) ext[i] = -ext[i];
00412 
00413             if (linkToGenParticles_ && abs(ext[i]) >= 2) {
00414                 // Link to the genParticle if possible, but not decays in flight (in ppMuX they're in GEN block, but they have wrong parameters)
00415                 if (!tp->genParticle().empty() && abs(ext[i]) >= 5) {
00416                     muToPrimary[i] = fetch(genBarcodes, tp->genParticle()[0]->barcode());
00417                 } else {
00418                     // Don't put the same trackingParticle twice!
00419                     int &indexPlus1 = tpToSecondaries[tp]; // will create a 0 if the tp is not in the list already
00420                     if (indexPlus1 == 0) indexPlus1 = convertAndPush(*tp, *secondaries, getTpMother(tp), genParticles, genBarcodes) + 1;
00421                     muToSecondary[i] = indexPlus1 - 1; 
00422                 }
00423             }
00424             edm::LogVerbatim("MuonMCClassifier") <<"\t Extended classification code = " << ext[i];
00425         }
00426     }
00427     
00428     writeValueMap(iEvent, muons, classif,   "");
00429     writeValueMap(iEvent, muons, ext,       "ext");
00430     writeValueMap(iEvent, muons, flav,      "flav");
00431     writeValueMap(iEvent, muons, tpId,      "tpId");
00432     writeValueMap(iEvent, muons, hitsPdgId, "hitsPdgId");
00433     writeValueMap(iEvent, muons, momPdgId,  "momPdgId");
00434     writeValueMap(iEvent, muons, momStatus, "momStatus");
00435     writeValueMap(iEvent, muons, momFlav,   "momFlav");
00436     writeValueMap(iEvent, muons, gmomPdgId, "gmomPdgId");
00437     writeValueMap(iEvent, muons, gmomFlav,  "gmomFlav");
00438     writeValueMap(iEvent, muons, hmomFlav,  "hmomFlav");
00439     writeValueMap(iEvent, muons, prodRho,   "prodRho");
00440     writeValueMap(iEvent, muons, prodZ,     "prodZ");
00441     writeValueMap(iEvent, muons, momRho,    "momRho");
00442     writeValueMap(iEvent, muons, momZ,      "momZ");
00443     writeValueMap(iEvent, muons, tpAssoQuality, "tpAssoQuality");
00444 
00445     if (linkToGenParticles_) {
00446         edm::OrphanHandle<reco::GenParticleCollection> secHandle = iEvent.put(secondaries, "secondaries");
00447         edm::RefProd<reco::GenParticleCollection> priRP(genParticles); 
00448         edm::RefProd<reco::GenParticleCollection> secRP(secHandle);
00449         std::auto_ptr<edm::Association<reco::GenParticleCollection> > outPri(new edm::Association<reco::GenParticleCollection>(priRP));
00450         std::auto_ptr<edm::Association<reco::GenParticleCollection> > outSec(new edm::Association<reco::GenParticleCollection>(secRP));
00451         edm::Association<reco::GenParticleCollection>::Filler fillPri(*outPri), fillSec(*outSec);
00452         fillPri.insert(muons, muToPrimary.begin(),   muToPrimary.end());
00453         fillSec.insert(muons, muToSecondary.begin(), muToSecondary.end());
00454         fillPri.fill(); fillSec.fill();
00455         iEvent.put(outPri, "toPrimaries");
00456         iEvent.put(outSec, "toSecondaries");
00457     }
00458 }    
00459 
00460 template<typename T>
00461 void
00462 MuonMCClassifier::writeValueMap(edm::Event &iEvent,
00463         const edm::Handle<edm::View<reco::Muon> > & handle,
00464         const std::vector<T> & values,
00465         const std::string    & label) const 
00466 {
00467     using namespace edm; 
00468     using namespace std;
00469     auto_ptr<ValueMap<T> > valMap(new ValueMap<T>());
00470     typename edm::ValueMap<T>::Filler filler(*valMap);
00471     filler.insert(handle, values.begin(), values.end());
00472     filler.fill();
00473     iEvent.put(valMap, label);
00474 }
00475 
00476 int
00477 MuonMCClassifier::flavour(int pdgId) const {
00478     int flav = abs(pdgId);
00479     // for quarks, leptons and bosons except gluons, take their pdgId
00480     // muons and taus have themselves as flavour
00481     if (flav <= 37 && flav != 21) return flav;
00482     // look for barions
00483     int bflav = ((flav / 1000) % 10);
00484     if (bflav != 0) return bflav;
00485     // look for mesons
00486     int mflav = ((flav / 100) % 10);
00487     if (mflav != 0) return mflav;
00488     return 0;
00489 }
00490 
00491 int MuonMCClassifier::fetch(const edm::Handle<std::vector<int> > & genBarcodes, int barcode) const
00492 {
00493     std::vector<int>::const_iterator it = std::find(genBarcodes->begin(), genBarcodes->end(), barcode);
00494     return (it == genBarcodes->end()) ? -1 : (it - genBarcodes->begin());
00495 }
00496 
00497 // push secondary in collection.
00498 // if it has a primary mother link to it. 
00499 int MuonMCClassifier::convertAndPush(const TrackingParticle &tp, 
00500                                      reco::GenParticleCollection &out,
00501                                      const TrackingParticleRef & simMom, 
00502                                      const edm::Handle<reco::GenParticleCollection> & genParticles,
00503                                      const edm::Handle<std::vector<int> > & genBarcodes) const {
00504     out.push_back(reco::GenParticle(tp.charge(), tp.p4(), tp.vertex(), tp.pdgId(), tp.status(), true));
00505     if (simMom.isNonnull() && !simMom->genParticle().empty()) {
00506         int momIdx = fetch(genBarcodes, simMom->genParticle()[0]->barcode());
00507         if (momIdx != -1)  out.back().addMother(reco::GenParticleRef(genParticles, momIdx));
00508     }
00509     return out.size()-1;
00510 }
00511 
00512 //define this as a plug-in
00513 DEFINE_FWK_MODULE(MuonMCClassifier);
/data/refman/pasoursint/CMSSW_5_3_10_patch2/src/MuonAnalysis/MuonAssociators/plugins/MuonMCClassifier.cc
