Skip to content

Pfpart cmssw 11 3 0 pre3 #47

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
jkiesele wants to merge 5 commits into
base: pepr_CMSSW_11_3_0_pre3
Choose a base branch
from
Open

Pfpart cmssw 11 3 0 pre3 #47

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
5883aa6
Merged pepr_CMSSW_11_3_0_pre3 from repository kdlong with cms-merge-t…
jkiesele Jun 18, 2021
0efca38
Merge branch 'pepr_CMSSW_11_3_0_pre3' of https://github.com/kdlong/cm…
jkiesele Jun 18, 2021
985e422
added some members to pftruthparticle
jkiesele Jun 18, 2021
793b3db
a step towards pf particle
jkiesele Jun 21, 2021
83ce949
fixes
jkiesele Jun 21, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 3 additions & 1 deletion DPGAnalysis/PFNanoAOD/python/pfTruth_cff.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down
237 changes: 181 additions & 56 deletions RecoParticleFlow/PFTruthProducer/plugins/PFTruthParticleProducer.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 <set>

//
// 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 <softkillcone>
* and tracks that contribute with less than a relative <softkillrelenergy> to
* the total particle energy are removed.
*
* This splitting of neutral particles depends on the granularity of the SimClusters
* in the respective subdetector.
* <add other possible options as configuration parameters>
*
* In case of charged particles, there are multiple options, e.g. when to split off bremsstrahlung
* for an electron.
*
*
*
*
*/

typedef edm::AssociationMap<edm::OneToMany<
TrackingParticleCollection, SimClusterCollection>> TrackingParticleToSimCluster;
Expand All @@ -38,20 +78,27 @@ class PFTruthParticleProducer : public edm::global::EDProducer<> {

private:
void produce(edm::StreamID, edm::Event &, const edm::EventSetup &) const override;
std::set<TrackingParticleRef> findTrackingParticleMatch(
std::unordered_map<unsigned int, TrackingParticleRef>& trackIdToTPRef, SimClusterRef scRef) const;

const SimTrack* getRoot(const SimTrack * st,
const std::vector<SimTrack> & simtracks,
const std::vector<SimVertex> & simvertices,
const std::map<int,int>& trackIdToTrackIdxAsso ) const;

edm::EDGetTokenT<TrackingParticleCollection> tpCollectionToken_;
edm::EDGetTokenT<SimClusterCollection> scCollectionToken_;
edm::EDGetTokenT<CaloParticleCollection> cpCollectionToken_;
edm::EDGetTokenT<std::vector<SimVertex> > svCollectionToken_;
edm::EDGetTokenT<std::vector<SimTrack> > stCollectionToken_;


};

PFTruthParticleProducer::PFTruthParticleProducer(const edm::ParameterSet &pset)
: tpCollectionToken_(consumes<TrackingParticleCollection>(pset.getParameter<edm::InputTag>("trackingParticles"))),
scCollectionToken_(consumes<SimClusterCollection>(pset.getParameter<edm::InputTag>("simClusters")))
cpCollectionToken_(consumes<CaloParticleCollection>(pset.getParameter<edm::InputTag>("caloParticles"))),
svCollectionToken_(consumes<std::vector<SimVertex> >(pset.getParameter<edm::InputTag>("simVertices"))),
stCollectionToken_(consumes<std::vector<SimTrack> >(pset.getParameter<edm::InputTag>("simTracks")))
{
produces<PFTruthParticleCollection>();
produces<edm::Association<PFTruthParticleCollection>>("trackingPartToPFTruth");
produces<edm::Association<PFTruthParticleCollection>>("simClusToPFTruth");
}

PFTruthParticleProducer::~PFTruthParticleProducer() {}
Expand All @@ -60,72 +107,150 @@ PFTruthParticleProducer::~PFTruthParticleProducer() {}
// member functions
//

std::set<TrackingParticleRef> PFTruthParticleProducer::findTrackingParticleMatch(
std::unordered_map<unsigned int, TrackingParticleRef>& trackIdToTPRef, SimClusterRef scRef) const {
std::set<TrackingParticleRef> 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<SimTrack> & simtracks,
const std::vector<SimVertex> & simvertices,
const std::map<int,int>& 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<TrackingParticleCollection> tpCollection;
iEvent.getByToken(tpCollectionToken_, tpCollection);

edm::Handle<SimClusterCollection> scCollection;
iEvent.getByToken(scCollectionToken_, scCollection);

auto out = std::make_unique<PFTruthParticleCollection>();
// Not used right now, but useful for the trackID lookup later
std::unordered_map<unsigned int, TrackingParticleRef> trackIdToTPRef;

std::vector<int> 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<CaloParticleCollection> cpCollection;
iEvent.getByToken(cpCollectionToken_, cpCollection);

edm::Handle<std::vector<SimVertex> > svCollection;
iEvent.getByToken(svCollectionToken_, svCollection);

edm::Handle<std::vector<SimTrack> > stCollection;
iEvent.getByToken(stCollectionToken_, stCollection);

std::map<int,int> trackIdToTrackIdxAsso;
for(size_t i=0;i<stCollection->size();i++){
trackIdToTrackIdxAsso[stCollection->at(i).trackId()] = i;
}

//unrealistic / ideal PF Truth: create one PFParticle form Geant handover aka CaloParticle

std::vector<size_t> 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<int> 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<PFTruthParticleCollection> 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<edm::Association<PFTruthParticleCollection>>(pfTruthHand);
auto scAssoc = std::make_unique<edm::Association<PFTruthParticleCollection>>(pfTruthHand);

edm::Association<PFTruthParticleCollection>::Filler tpFiller(*tpAssoc);
tpFiller.insert(tpCollection, trackingIndices.begin(), trackingIndices.end());
tpFiller.fill();

edm::Association<PFTruthParticleCollection>::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
Expand Down
Loading