Skip to content
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

Jump to bottom

Calo particle validation with graph #22596

Merged
merged 28 commits into from Apr 3, 2018
Merged

Calo particle validation with graph #22596

Show file tree
Hide file tree
Changes from 27 commits
Commits
Show all changes
28 commits
Select commit Hold shift + click to select a range
8a1ec7e
Add CaloParticle Debug package
rovere Nov 23, 2017
d263af2
Minor stylistic changes
rovere Nov 24, 2017
6c5f05e
First working implementation of CaloP with Graphs.
rovere Dec 4, 2017
d2e1c22
Cleanup unused structures
rovere Dec 4, 2017
66b57e8
More code cleanup
rovere Dec 4, 2017
fbdd6bf
Code reformatting with clang
rovere Dec 5, 2017
b97765b
Remove energy from the edge properties
rovere Dec 5, 2017
73b8869
Final version
rovere Dec 7, 2017
f8db957
Remove gap in vertex numbering in the final graph
rovere Dec 11, 2017
d2950d3
Merge header into cc file
rovere Dec 13, 2017
77836a5
Simplify the Debugging plugin
rovere Dec 13, 2017
6a80062
Actually use the CaloTruthAccumulator with Graph
rovere Dec 13, 2017
b573803
Remove unnecessary cfi file
rovere Dec 15, 2017
13a77cf
Use meaningful filenames
rovere Dec 15, 2017
b91785a
Address some review comments
rovere Dec 15, 2017
59e3314
Better format for Graphviz output
rovere Dec 15, 2017
0e15c13
Fix GCC flags and dot output
rovere Dec 19, 2017
89857cb
Add cumulative hits to edges in dot-graphs
rovere Dec 19, 2017
5b394a8
Monitor the difference between {,self}Energy
rovere Feb 9, 2018
5c9cf9d
Order collection while debugging
rovere Feb 9, 2018
c853c2a
Debug also SimClusters
rovere Feb 19, 2018
2e157bc
Improve CaloParticle Validation
rovere Mar 13, 2018
fd504b2
Add more CaloParticles
rovere Mar 13, 2018
d759937
Recover CaloParticle into SimPFProducer
rovere Mar 12, 2018
506b3a4
Actually implement changes of 4f139f44cb515c3cad7137fa6ca22c38ce96f772
rovere Mar 12, 2018
59c4c11
Apply code-checks
rovere Mar 13, 2018
1bd7e49
Remove unused dependency
rovere Mar 13, 2018
8b0ab5c
Implement Review comments 1
rovere Mar 14, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
82 changes: 48 additions & 34 deletions RecoParticleFlow/PFSimProducer/plugins/SimPFProducer.cc
View file
Open in desktop
Expand Up @@ -56,22 +56,22 @@
#include "FWCore/Utilities/interface/transform.h"

class SimPFProducer : public edm::global::EDProducer<> {
public:
public:
SimPFProducer(const edm::ParameterSet&);
~SimPFProducer() override { }

void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
private:

private:
// parameters
const double superClusterThreshold_, neutralEMThreshold_, neutralHADThreshold_;
const bool useTiming_;

// inputs
const edm::EDGetTokenT<edm::View<reco::PFRecTrack> > pfRecTracks_;
const edm::EDGetTokenT<edm::View<reco::Track> > tracks_;
const edm::EDGetTokenT<edm::View<reco::Track> > gsfTracks_;
const edm::EDGetTokenT<reco::MuonCollection> muons_;
const edm::EDGetTokenT<reco::MuonCollection> muons_;
const edm::EDGetTokenT<edm::ValueMap<float>> srcTrackTime_, srcTrackTimeError_;
const edm::EDGetTokenT<edm::ValueMap<float>> srcGsfTrackTime_, srcGsfTrackTimeError_;
const edm::EDGetTokenT<TrackingParticleCollection> trackingParticles_;
Expand All @@ -86,7 +86,7 @@ class SimPFProducer : public edm::global::EDProducer<> {
DEFINE_FWK_MODULE(SimPFProducer);

namespace {

template<typename T>
uint64_t hashSimInfo(const T& simTruth,size_t i = 0) {
uint64_t evtid = simTruth.eventId().rawId();
Expand Down Expand Up @@ -119,15 +119,15 @@ SimPFProducer::SimPFProducer(const edm::ParameterSet& conf) :
produces<reco::PFCandidateCollection>();
}

void SimPFProducer::produce(edm::StreamID, edm::Event& evt, const edm::EventSetup& es) const {
void SimPFProducer::produce(edm::StreamID, edm::Event& evt, const edm::EventSetup& es) const {
//get associators
std::vector<edm::Handle<reco::TrackToTrackingParticleAssociator> > associators;
for( const auto& token : associators_ ) {
associators.emplace_back();
auto& back = associators.back();
evt.getByToken(token,back);
}

//get PFRecTrack
edm::Handle<edm::View<reco::PFRecTrack> > PFTrackCollectionH;
evt.getByToken(pfRecTracks_,PFTrackCollectionH);
Expand All @@ -137,7 +137,7 @@ void SimPFProducer::produce(edm::StreamID, edm::Event& evt, const edm::EventSetu
const auto ptr = PFTrackCollection.ptrAt(i);
PFTrackToGeneralTrack.insert(ptr->trackRef().key());
}

//get track collections
edm::Handle<edm::View<reco::Track> > TrackCollectionH;
evt.getByToken(tracks_, TrackCollectionH);
Expand All @@ -160,7 +160,7 @@ void SimPFProducer::produce(edm::StreamID, edm::Event& evt, const edm::EventSetu
evt.getByToken(srcGsfTrackTime_, gsfTrackTimeH);
evt.getByToken(srcGsfTrackTimeError_, gsfTrackTimeErrH);
}

//get tracking particle collections
edm::Handle<TrackingParticleCollection> TPCollectionH;
evt.getByToken(trackingParticles_, TPCollectionH);
Expand All @@ -179,15 +179,15 @@ void SimPFProducer::produce(edm::StreamID, edm::Event& evt, const edm::EventSetu
evt.getByToken(simClusters_,SimClustersH);
const std::vector<reco::PFCluster>& SimClusters = *SimClustersH;

std::unordered_map<uint64_t,size_t> hashToSimCluster;
std::unordered_map<uint64_t,size_t> hashToSimCluster;

for( unsigned i = 0; i < SimClustersTruth.size(); ++i ) {
const auto& simTruth = SimClustersTruth[i];
hashToSimCluster[hashSimInfo(simTruth)] = i;
}

// associate the reco tracks / gsf Tracks
std::vector<reco::RecoToSimCollection> associatedTracks, associatedTracksGsf;
std::vector<reco::RecoToSimCollection> associatedTracks, associatedTracksGsf;
for( auto associator : associators ) {
associatedTracks.emplace_back(associator->associateRecoToSim(TrackCollectionH, TPCollectionH));
//associatedTracksGsf.emplace_back(associator->associateRecoToSim(GsfTrackCollectionH, TPCollectionH));
Expand All @@ -213,18 +213,18 @@ void SimPFProducer::produce(edm::StreamID, edm::Event& evt, const edm::EventSetu
auto pdgId = std::abs(simc->pdgId());
edm::Ref<std::vector<reco::PFCluster> > clusterRef(SimClustersH,simc.key());
if( ( (pdgId == 22 || pdgId == 11) && clusterRef->energy() > neutralEMThreshold_) ||
clusterRef->energy() > neutralHADThreshold_ ) {
clusterRef->energy() > neutralHADThreshold_ ) {
good_simclusters.push_back(isc);
etot += clusterRef->energy();
pttot += clusterRef->pt();
pttot += clusterRef->pt();
auto bec = std::make_unique<reco::PFBlockElementCluster>(clusterRef,reco::PFBlockElement::HGCAL);
block.addElement(bec.get());
simCluster2Block[simc.key()] = icp;
simCluster2BlockIndex[simc.key()] = bec->index();
caloParticle2SimCluster.emplace(icp,simc.key());
caloParticle2SimCluster.emplace(CaloParticles[icp].g4Tracks()[0].trackId(), simc.key());
}
}

auto pdgId = std::abs(CaloParticles[icp].pdgId());

caloParticle2SuperCluster.push_back(-1);
Expand All @@ -245,12 +245,13 @@ void SimPFProducer::produce(edm::StreamID, edm::Event& evt, const edm::EventSetu
superclusters->emplace_back(etot,seedpos,seed,clusters);
}
}

auto blocksHandle = evt.put(std::move(blocks));
auto superClustersHandle = evt.put(std::move(superclusters),"perfect");

// list tracks so we can mark them as used and/or fight over them
std::vector<bool> usedTrack(TrackCollection.size(),false),
//usedGsfTrack(GsfTrackCollection.size(),false),
std::vector<bool> usedTrack(TrackCollection.size(),false),
//usedGsfTrack(GsfTrackCollection.size(),false),
usedSimCluster(SimClusters.size(),false);

auto candidates = std::make_unique<reco::PFCandidateCollection>();
Expand All @@ -267,14 +268,14 @@ void SimPFProducer::produce(edm::StreamID, edm::Event& evt, const edm::EventSetu
if( assoc_tps == associatedTracks.back().end() ) continue;
// assured now that we are matched to a set of tracks
const auto& matches = assoc_tps->val;

const auto absPdgId = std::abs(matches[0].first->pdgId());
const auto charge = tkRef->charge();
const auto three_mom = tkRef->momentum();
constexpr double mpion2 = 0.13957*0.13957;
double energy = std::sqrt(three_mom.mag2() + mpion2);
math::XYZTLorentzVector trk_p4(three_mom.x(),three_mom.y(),three_mom.z(),energy);

reco::PFCandidate::ParticleType part_type;

switch( absPdgId ) {
Expand All @@ -287,22 +288,22 @@ void SimPFProducer::produce(edm::StreamID, edm::Event& evt, const edm::EventSetu
default:
part_type = reco::PFCandidate::h;
}

candidates->emplace_back(charge, trk_p4, part_type);
auto& candidate = candidates->back();

candidate.setTrackRef(tkRef.castTo<reco::TrackRef>());

if (useTiming_) candidate.setTime( (*trackTimeH)[tkRef], (*trackTimeErrH)[tkRef] );

// bind to cluster if there is one and try to gather conversions, etc
for( const auto& match : matches ) {
for( const auto& match : matches ) {
uint64_t hash = hashSimInfo(*(match.first));
if( hashToSimCluster.count(hash) ) {
if( hashToSimCluster.count(hash) ) {
auto simcHash = hashToSimCluster[hash];
if( !usedSimCluster[simcHash] ) {
if( simCluster2Block.count(simcHash) &&

if( !usedSimCluster[simcHash] ) {
if( simCluster2Block.count(simcHash) &&
simCluster2BlockIndex.count(simcHash) ) {
size_t block = simCluster2Block.find(simcHash)->second;
size_t blockIdx = simCluster2BlockIndex.find(simcHash)->second;
Expand All @@ -324,21 +325,34 @@ void SimPFProducer::produce(edm::StreamID, edm::Event& evt, const edm::EventSetu
usedSimCluster[ref.key()] = true;
}
}

//*TODO* cluster time is not reliable at the moment, so just keep time from the track if available

}
}
}
}
// Now try to include also electrons that have been reconstructed using the GraphCaloParticles
else if (caloParticle2SimCluster.count(match.first->g4Tracks()[0].trackId())) {
auto range = caloParticle2SimCluster.equal_range(match.first->g4Tracks()[0].trackId());
for (auto it = range.first; it != range.second; ++it) {
if (!usedSimCluster[it->second]) {
usedSimCluster[it->second] = true;
size_t block = simCluster2Block.find(it->second)->second;
size_t blockIdx = simCluster2BlockIndex.find(it->second)->second;
edm::Ref<reco::PFBlockCollection> blockRef(blocksHandle,block);
candidate.addElementInBlock(blockRef,blockIdx);
}
}
}
}
usedTrack[tkRef.key()] = true;
usedTrack[tkRef.key()] = true;
// remove tracks already used by muons
if( MuonTrackToGeneralTrack.count(itk) || absPdgId == 13)
candidates->pop_back();
}

// now loop over the non-collected clusters in blocks
// now loop over the non-collected clusters in blocks
// and turn them into neutral hadrons or photons
const auto& theblocks = *blocksHandle;
for( unsigned ibl = 0; ibl < theblocks.size(); ++ibl ) {
Expand Down Expand Up @@ -367,6 +381,6 @@ void SimPFProducer::produce(edm::StreamID, edm::Event& evt, const edm::EventSetu
}
}
}

evt.put(std::move(candidates));
}