diff --git a/DPGAnalysis/PFNanoAOD/python/pfTruth_cff.py b/DPGAnalysis/PFNanoAOD/python/pfTruth_cff.py index 80107c1b9ed68..e0bf7a206dc9b 100644 --- a/DPGAnalysis/PFNanoAOD/python/pfTruth_cff.py +++ b/DPGAnalysis/PFNanoAOD/python/pfTruth_cff.py @@ -4,8 +4,10 @@ from DPGAnalysis.TrackNanoAOD.trackingParticles_cff import trackingParticleTable pfTruthParticles = cms.EDProducer("PFTruthParticleProducer", - simClusters = cms.InputTag("mix:MergedCaloTruth"), trackingParticles= cms.InputTag("mix:MergedTrackTruth"), + caloParticles= cms.InputTag("mix:MergedCaloTruth"), + simVertices= cms.InputTag("SimVertices"), + simTracks= cms.InputTag("SimTracks"), ) pfTruthTable = cms.EDProducer("SimplePFTruthParticleFlatTableProducer", diff --git a/RecoParticleFlow/PFTruthProducer/plugins/PFTruthParticleProducer.cc b/RecoParticleFlow/PFTruthProducer/plugins/PFTruthParticleProducer.cc index 317af99038ed8..0b5177a0604a4 100644 --- a/RecoParticleFlow/PFTruthProducer/plugins/PFTruthParticleProducer.cc +++ b/RecoParticleFlow/PFTruthProducer/plugins/PFTruthParticleProducer.cc @@ -16,17 +16,57 @@ #include "FWCore/MessageLogger/interface/MessageLogger.h" #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h" +#include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h" +#include "SimDataFormats/CaloAnalysis/interface/CaloParticleFwd.h" #include "SimDataFormats/CaloAnalysis/interface/SimClusterFwd.h" #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h" #include "SimDataFormats/PFAnalysis/interface/PFTruthParticle.h" #include "SimDataFormats/PFAnalysis/interface/PFTruthParticleFwd.h" + #include "FWCore/Utilities/interface/EDGetToken.h" #include // -// class decleration +// class declaration // +/* + * + * produces in two steps PFTruthParticles. + * This producer considers truth particles one-by-one and does not + * account for possible overlaps. These will need to be resolved in a second producer + * (PFTruthParticleMerger). + * The advantage is that the plugin here can run on every single event, independent of + * pileup mixing and event density, while the second producer can be run after mixing in pileup + * + * Step A) + * Produces an ideal PFTruthParticle, without considering detector limitations. + * This is done by connecting every trackingParticle and SimCluster that originated + * from a particle passed to Geant back to it. + * No configuration options needed. + * The tracking particles should be selected by one of the standard tracking selectors, e.g.: + * https://github.com/cms-sw/cmssw/blob/master/Validation/RecoTrack/python/TrackingParticleSelectionForEfficiency_cfi.py + * + * Step B) + * The ideal PFTruthParticle gets split according to detector limitations, e.g. if + * a neutral pion splits to 2 photons and leaves to distinct photon showers, + * this step would produce 2 photons. + * + * For neutral and charged particles, stray clusters with deltaR larger than + * and tracks that contribute with less than a relative to + * the total particle energy are removed. + * + * This splitting of neutral particles depends on the granularity of the SimClusters + * in the respective subdetector. + * + * + * In case of charged particles, there are multiple options, e.g. when to split off bremsstrahlung + * for an electron. + * + * + * + * + */ typedef edm::AssociationMap> TrackingParticleToSimCluster; @@ -38,20 +78,27 @@ class PFTruthParticleProducer : public edm::global::EDProducer<> { private: void produce(edm::StreamID, edm::Event &, const edm::EventSetup &) const override; - std::set findTrackingParticleMatch( - std::unordered_map& trackIdToTPRef, SimClusterRef scRef) const; + + const SimTrack* getRoot(const SimTrack * st, + const std::vector & simtracks, + const std::vector & simvertices, + const std::map& trackIdToTrackIdxAsso ) const; edm::EDGetTokenT tpCollectionToken_; - edm::EDGetTokenT scCollectionToken_; + edm::EDGetTokenT cpCollectionToken_; + edm::EDGetTokenT > svCollectionToken_; + edm::EDGetTokenT > stCollectionToken_; + + }; PFTruthParticleProducer::PFTruthParticleProducer(const edm::ParameterSet &pset) : tpCollectionToken_(consumes(pset.getParameter("trackingParticles"))), - scCollectionToken_(consumes(pset.getParameter("simClusters"))) + cpCollectionToken_(consumes(pset.getParameter("caloParticles"))), + svCollectionToken_(consumes >(pset.getParameter("simVertices"))), + stCollectionToken_(consumes >(pset.getParameter("simTracks"))) { produces(); - produces>("trackingPartToPFTruth"); - produces>("simClusToPFTruth"); } PFTruthParticleProducer::~PFTruthParticleProducer() {} @@ -60,72 +107,150 @@ PFTruthParticleProducer::~PFTruthParticleProducer() {} // member functions // -std::set PFTruthParticleProducer::findTrackingParticleMatch( - std::unordered_map& trackIdToTPRef, SimClusterRef scRef) const { - std::set trackingParticles; - for (auto& track : scRef->g4Tracks()) { - unsigned int trackId = track.trackId(); - if (trackIdToTPRef.find(trackId) != trackIdToTPRef.end()) - trackingParticles.insert(trackIdToTPRef[trackId]); + +const SimTrack* PFTruthParticleProducer::getRoot( + const SimTrack * st, + const std::vector & simtracks, + const std::vector & simvertices, + const std::map& trackIdToTrackIdxAsso ) const{ + + const SimTrack* out=st; + int vidx = out->vertIndex();//starting point + while(true){ + if(vidx<0) + break; //no vertex + const SimVertex& vertex = simvertices.at(vidx); + int stid = vertex.parentIndex();//this is Geant track ID, not vector index + int stidx = trackIdToTrackIdxAsso.at(stid); //get vector index + if(stidx<0) + break; + const SimTrack & simtrack = simtracks.at(stidx); + vidx = simtrack.vertIndex(); + out=&simtrack; } - return trackingParticles; + return out; } + // ------------ method called to produce the data ------------ void PFTruthParticleProducer::produce(edm::StreamID, edm::Event &iEvent, const edm::EventSetup &iSetup) const { edm::Handle tpCollection; iEvent.getByToken(tpCollectionToken_, tpCollection); - edm::Handle scCollection; - iEvent.getByToken(scCollectionToken_, scCollection); - - auto out = std::make_unique(); - // Not used right now, but useful for the trackID lookup later - std::unordered_map trackIdToTPRef; - - std::vector trackingIndices; - for (size_t i = 0; i < tpCollection->size(); i++) { - TrackingParticleRef trackingParticle(tpCollection, i); - PFTruthParticle particle; - particle.addTrackingParticle(trackingParticle); - particle.setCharge(trackingParticle->charge()); - particle.setPdgId(trackingParticle->pdgId()); - particle.setP4(trackingParticle->p4()); - out->push_back(particle); - trackingIndices.push_back(i); - for (auto& track : trackingParticle->g4Tracks()) { - unsigned int trackId = track.trackId(); - trackIdToTPRef[trackId] = trackingParticle; + edm::Handle cpCollection; + iEvent.getByToken(cpCollectionToken_, cpCollection); + + edm::Handle > svCollection; + iEvent.getByToken(svCollectionToken_, svCollection); + + edm::Handle > stCollection; + iEvent.getByToken(stCollectionToken_, stCollection); + + std::map trackIdToTrackIdxAsso; + for(size_t i=0;isize();i++){ + trackIdToTrackIdxAsso[stCollection->at(i).trackId()] = i; + } + + //unrealistic / ideal PF Truth: create one PFParticle form Geant handover aka CaloParticle + + std::vector matchedtptocp; + PFTruthParticleCollection idealPFTruth; + + for(const auto& cp: *cpCollection){ + + TrackingParticleRefVector tprefs; + const SimClusterRefVector& screfs = cp.simClusters(); + + //match the tracking particle(s) + for(size_t i=0; i< tpCollection->size(); i++){ + if(std::find(matchedtptocp.begin(),matchedtptocp.end(), i) != matchedtptocp.end()) + continue; //already matched + + TrackingParticleRef tpref(tpCollection, i); + + //go through whole history + auto stp = getRoot(& tpref->g4Tracks().at(0), + *stCollection,*svCollection, trackIdToTrackIdxAsso); + + //match + if(stp->trackId() == cp.g4Tracks().at(0).trackId() + && stp->eventId() == cp.g4Tracks().at(0).eventId()){ + matchedtptocp.push_back(i); + tprefs.push_back(tpref); + } } + PFTruthParticle pftp(tprefs,screfs); + pftp.setP4(cp.p4()); + pftp.setPdgId(cp.pdgId()); + pftp.setCharge(cp.charge()); + int vidx = cp.g4Tracks().at(0).vertIndex(); + auto vertpos = svCollection->at(vidx).position(); + + pftp.setVertex(PFTruthParticle::LorentzVectorF( + vertpos.x(),vertpos.y(),vertpos.z(),vertpos.t())); + + idealPFTruth.push_back(pftp); } - std::vector scIndices; - for (size_t i = 0; i < scCollection->size(); i++) { - SimClusterRef simclus(scCollection, i); - // Doesn't really do much anyway - //findTrackingParticleMatch(trackIdToTPRef, simclus)) - // For now just randomly assign the SCs as a dummy - int tpIdx = i % out->size(); - auto& pf = out->at(tpIdx); - pf.addSimCluster(simclus); - scIndices.push_back(tpIdx); + + // --- now we have the ideal PF truth, far from realistic. + // --- this is where the actual "meat" of the implementation starts + + double softkill_relthreshold = 0.001;//just to avoid ambiguities. + double hardsplit_deltaR = 0.3; + + + std::unique_ptr PFtruth; + + //apply splitting particle by particle + //this is where the physics happens + for(const auto & pftp: idealPFTruth){ + + + //consider the ones without tracking particles // this is not the same as the ones without charge, + // because e.g. a very early converting photon would have charge 0, but two tracking particles associated + if(! pftp.trackingParticles().size()) { + for(const auto& sc: pftp.simClusters()){ + //create PF particle per SimCluster + //vertex position will remain the initial CP one (timing) + auto npftp = pftp; + npftp.trackingParticles().clear(); + npftp.simClusters().clear(); + npftp.addSimCluster(sc); + //this is now defined to be a neutral PF Particle. Even if the SC happens to be charged + npftp.setCharge(0); + npftp.setP4(sc->impactMomentum());//that's all we measure + npftp.setPdgId(sc->pdgId()); //to be discussed + PFtruth->push_back(npftp); + } + continue; + } + //remaining: the ones that have tracking particles associated + + + //some more cases, like far away brem... + + + //trim + auto trimmed_pfp = pftp; + SimClusterRefVector trimmersc; + for(auto & sc: trimmed_pfp.simClusters()){ + if(sc->impactMomentum().E() > trimmed_pfp.p4().E() * softkill_relthreshold){ + trimmersc.push_back(sc); + } + } + trimmed_pfp.setSimClusters(trimmersc); + //no split + PFtruth->push_back(trimmed_pfp); } - auto pfTruthHand = iEvent.put(std::move(out)); - auto tpAssoc = std::make_unique>(pfTruthHand); - auto scAssoc = std::make_unique>(pfTruthHand); - edm::Association::Filler tpFiller(*tpAssoc); - tpFiller.insert(tpCollection, trackingIndices.begin(), trackingIndices.end()); - tpFiller.fill(); - edm::Association::Filler scFiller(*scAssoc); - scFiller.insert(scCollection, scIndices.begin(), scIndices.end()); - scFiller.fill(); - iEvent.put(std::move(tpAssoc), "trackingPartToPFTruth"); - iEvent.put(std::move(scAssoc), "simClusToPFTruth"); + iEvent.put(std::move(PFtruth), "PFTruthParticles"); + + } // define this as a plug-in diff --git a/SimDataFormats/PFAnalysis/interface/PFTruthParticle.h b/SimDataFormats/PFAnalysis/interface/PFTruthParticle.h index 50ea4f28c3902..01198ca496664 100644 --- a/SimDataFormats/PFAnalysis/interface/PFTruthParticle.h +++ b/SimDataFormats/PFAnalysis/interface/PFTruthParticle.h @@ -22,9 +22,8 @@ class PFTruthParticle { friend std::ostream& operator<<(std::ostream& s, PFTruthParticle const& tp); public: - typedef math::XYZTLorentzVectorD LorentzVector; ///< Lorentz vector - typedef math::XYZPointD Point; ///< point in the space - typedef math::XYZVectorD Vector; ///< point in the space + typedef math::XYZTLorentzVectorF LorentzVectorF; ///< Lorentz vector + typedef math::XYZVectorF VectorF; ///< point in the space /** @brief Default constructor. Note that the object will be useless until it is provided * with a SimTrack and parent TrackingVertex. @@ -43,12 +42,15 @@ class PFTruthParticle { void setSimClusters(const SimClusterRefVector& refs); void setPdgId(int pdgId); void setCharge(int charge); - void setP4(LorentzVector p4); + void setP4(LorentzVectorF p4); + void setVertex(LorentzVectorF vertex); void addSimCluster(const SimClusterRef sc); void addTrackingParticle(const TrackingParticleRef tp); SimClusterRefVector& simClusters() { return simClusters_; } TrackingParticleRefVector& trackingParticles() { return trackingParticles_; } + const SimClusterRefVector& simClusters() const { return simClusters_; } + const TrackingParticleRefVector& trackingParticles() const { return trackingParticles_; } size_t nSimCluster() const { return simClusters_.size(); } size_t nTrackingParticle() const { return trackingParticles_.size(); } @@ -64,58 +66,61 @@ class PFTruthParticle { const std::vector& g4Tracks() const { return g4Tracks_; } - /// @brief Electric charge. Note this is taken from the first SimTrack only. + /// @brief Electric charge. float charge() const { return charge_; } - /// @brief Four-momentum Lorentz vector. Note this is taken from the first SimTrack only. - const LorentzVector& p4() const { return p4_; } + /// @brief Four-momentum Lorentz vector. + const LorentzVectorF& p4() const { return p4_; } + + /// @brief Vertex XYZT. + const LorentzVectorF& vertex() const { return vertex_; } /// @brief spatial momentum vector - Vector momentum() const { return p4().Vect(); } + VectorF momentum() const { return p4().Vect(); } - /// @brief Magnitude of momentum vector. Note this is taken from the first SimTrack only. + /// @brief Magnitude of momentum vector. double p() const { return p4().P(); } - /// @brief Energy. Note this is taken from the first SimTrack only. + /// @brief Energy. double energy() const { return p4().E(); } - /// @brief Transverse energy. Note this is taken from the first SimTrack only. + /// @brief Transverse energy. double et() const { return p4().Et(); } - /// @brief Mass. Note this is taken from the first SimTrack only. + /// @brief Mass. double mass() const { return p4().M(); } - /// @brief Mass squared. Note this is taken from the first SimTrack only. + /// @brief Mass squared. double massSqr() const { return pow(mass(), 2); } - /// @brief Transverse mass. Note this is taken from the first SimTrack only. + /// @brief Transverse mass. double mt() const { return p4().Mt(); } - /// @brief Transverse mass squared. Note this is taken from the first SimTrack only. + /// @brief Transverse mass squared. double mtSqr() const { return p4().Mt2(); } - /// @brief x coordinate of momentum vector. Note this is taken from the first SimTrack only. + /// @brief x coordinate of momentum vector. double px() const { return p4().Px(); } - /// @brief y coordinate of momentum vector. Note this is taken from the first SimTrack only. + /// @brief y coordinate of momentum vector. double py() const { return p4().Py(); } - /// @brief z coordinate of momentum vector. Note this is taken from the first SimTrack only. + /// @brief z coordinate of momentum vector. double pz() const { return p4().Pz(); } - /// @brief Transverse momentum. Note this is taken from the first SimTrack only. + /// @brief Transverse momentum. double pt() const { return p4().Pt(); } - /// @brief Momentum azimuthal angle. Note this is taken from the first SimTrack only. + /// @brief Momentum azimuthal angle. double phi() const { return p4().Phi(); } - /// @brief Momentum polar angle. Note this is taken from the first SimTrack only. + /// @brief Momentum polar angle. double theta() const { return p4().Theta(); } - /// @brief Momentum pseudorapidity. Note this is taken from the first SimTrack only. + /// @brief Momentum pseudorapidity. double eta() const { return p4().Eta(); } - /// @brief Rapidity. Note this is taken from the first SimTrack only. + /// @brief Rapidity. double rapidity() const { return p4().Rapidity(); } /// @brief Same as rapidity(). @@ -125,7 +130,7 @@ class PFTruthParticle { /// references to G4 and reco::GenParticle tracks int charge_; int pdgId_; - LorentzVector p4_; + LorentzVectorF p4_, vertex_; std::vector g4Tracks_; reco::GenParticleRefVector genParticles_; SimClusterRefVector simClusters_; diff --git a/SimDataFormats/PFAnalysis/src/PFTruthParticle.cc b/SimDataFormats/PFAnalysis/src/PFTruthParticle.cc index 9ca899f639d03..b533b09c250e4 100644 --- a/SimDataFormats/PFAnalysis/src/PFTruthParticle.cc +++ b/SimDataFormats/PFAnalysis/src/PFTruthParticle.cc @@ -43,6 +43,8 @@ void PFTruthParticle::setPdgId(int pdgId) { pdgId_ = pdgId; } void PFTruthParticle::setCharge(int charge) { charge_ = charge; } -void PFTruthParticle::setP4(LorentzVector p4) { p4_ = p4; } +void PFTruthParticle::setP4(LorentzVectorF p4) { p4_ = p4; } + +void PFTruthParticle::setVertex(LorentzVectorF vertex) { vertex_ = vertex; } void PFTruthParticle::addG4Track(const SimTrack& t) { g4Tracks_.push_back(t); }