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.

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

Speedup Track Validation #33462

Merged
merged 12 commits into from Apr 23, 2021
16 changes: 16 additions & 0 deletions SimDataFormats/Track/interface/UniqueSimTrackId.h
@@ -0,0 +1,16 @@
#ifndef SimDataFormatsTrackUniqueSimTrackId_H
#define SimDataFormatsTrackUniqueSimTrackId_H

#include "FWCore/Utilities/interface/hash_combine.h"

#include "SimDataFormats/EncodedEventId/interface/EncodedEventId.h"
#include <tuple>

using UniqueSimTrackId = std::pair<uint32_t, EncodedEventId>;
struct UniqueSimTrackIdHash {
std::size_t operator()(UniqueSimTrackId const &s) const noexcept {
return edm::hash_value(s.first, s.second.rawId());
}
};

#endif
Expand Up @@ -12,10 +12,10 @@
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"

class SimHitTPAssociationProducer : public edm::global::EDProducer<> {
class SimHitTPAssociationProducer final : public edm::global::EDProducer<> {
public:
typedef std::pair<TrackingParticleRef, TrackPSimHitRef> SimHitTPPair;
typedef std::vector<SimHitTPPair> SimHitTPAssociationList;
using SimHitTPPair = std::pair<TrackingParticleRef, TrackPSimHitRef>;
using SimHitTPAssociationList = std::vector<SimHitTPPair>;

explicit SimHitTPAssociationProducer(const edm::ParameterSet &);
~SimHitTPAssociationProducer() override;
Expand Down
44 changes: 28 additions & 16 deletions SimGeneral/TrackingAnalysis/plugins/SimHitTPAssociationProducer.cc
@@ -1,16 +1,24 @@
#include <fstream>
#include <iostream>
#include <memory>
#include <utility>
#include <algorithm>
#include <map>
#include <vector>

#include "FWCore/Utilities/interface/InputTag.h"

#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/global/EDProducer.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ParameterSetfwd.h"
#include "FWCore/Utilities/interface/EDGetToken.h"

#include "DataFormats/Common/interface/Handle.h"
#include "SimDataFormats/Track/interface/UniqueSimTrackId.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"

#include "SimGeneral/TrackingAnalysis/interface/SimHitTPAssociationProducer.h"

Expand All @@ -36,27 +44,31 @@ void SimHitTPAssociationProducer::produce(edm::StreamID, edm::Event &iEvent, con
iEvent.getByToken(_trackingParticleSrc, TPCollectionH);

// prepare temporary map between SimTrackId and TrackingParticle index
std::map<std::pair<size_t, EncodedEventId>, TrackingParticleRef> mapping;
for (TrackingParticleCollection::size_type itp = 0; itp < TPCollectionH.product()->size(); ++itp) {
TrackingParticleRef trackingParticle(TPCollectionH, itp);
std::unordered_map<UniqueSimTrackId, TrackingParticleRef, UniqueSimTrackIdHash> mapping;
auto const &tpColl = *TPCollectionH.product();
for (TrackingParticleCollection::size_type itp = 0, size = tpColl.size(); itp < size; ++itp) {
auto const &trackingParticle = tpColl[itp];
TrackingParticleRef trackingParticleRef(TPCollectionH, itp);
// SimTracks inside TrackingParticle
EncodedEventId eid(trackingParticle->eventId());
for (auto itrk = trackingParticle->g4Track_begin(); itrk != trackingParticle->g4Track_end(); ++itrk) {
std::pair<uint32_t, EncodedEventId> trkid(itrk->trackId(), eid);
mapping.insert(std::make_pair(trkid, trackingParticle));
EncodedEventId eid(trackingParticle.eventId());
for (auto const &trk : trackingParticle.g4Tracks()) {
UniqueSimTrackId trkid(trk.trackId(), eid);
mapping.insert(std::make_pair(trkid, trackingParticleRef));
}
}

// PSimHits
for (auto const &psit : _simHitSrc) {
edm::Handle<edm::PSimHitContainer> PSimHitCollectionH;
iEvent.getByToken(psit, PSimHitCollectionH);
for (unsigned int psimHit = 0; psimHit != PSimHitCollectionH->size(); ++psimHit) {
TrackPSimHitRef pSimHitRef(PSimHitCollectionH, psimHit);
std::pair<uint32_t, EncodedEventId> simTkIds(pSimHitRef->trackId(), pSimHitRef->eventId());
auto const &pSimHitCollection = *PSimHitCollectionH;
for (unsigned int psimHitI = 0, size = pSimHitCollection.size(); psimHitI < size; ++psimHitI) {
TrackPSimHitRef pSimHitRef(PSimHitCollectionH, psimHitI);
auto const &pSimHit = pSimHitCollection[psimHitI];
UniqueSimTrackId simTkIds(pSimHit.trackId(), pSimHit.eventId());
auto ipos = mapping.find(simTkIds);
if (ipos != mapping.end()) {
simHitTPList->push_back(std::make_pair(ipos->second, pSimHitRef));
simHitTPList->emplace_back(ipos->second, pSimHitRef);
}
}
}
Expand Down
Expand Up @@ -10,6 +10,8 @@
#include "SimTracker/TrackAssociation/interface/ParametersDefinerForTP.h"
#include <SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h>

#include <memory>

class CosmicParametersDefinerForTP : public ParametersDefinerForTP {
public:
CosmicParametersDefinerForTP(){};
Expand All @@ -22,6 +24,11 @@ class CosmicParametersDefinerForTP : public ParametersDefinerForTP {
const edm::EventSetup &iSetup,
const TrackingParticleRef &tpr) const override;

std::tuple<TrackingParticle::Vector, TrackingParticle::Point> momentumAndVertex(
const edm::Event &iEvent, const edm::EventSetup &iSetup, const TrackingParticleRef &tpr) const override {
return std::make_tuple(momentum(iEvent, iSetup, tpr), vertex(iEvent, iSetup, tpr));
}

TrackingParticle::Vector momentum(const edm::Event &iEvent,
const edm::EventSetup &iSetup,
const Charge ch,
Expand All @@ -43,7 +50,7 @@ class CosmicParametersDefinerForTP : public ParametersDefinerForTP {
}

std::unique_ptr<ParametersDefinerForTP> clone() const override {
return std::unique_ptr<CosmicParametersDefinerForTP>(new CosmicParametersDefinerForTP(*this));
return std::make_unique<CosmicParametersDefinerForTP>(*this);
}

private:
Expand Down
13 changes: 12 additions & 1 deletion SimTracker/TrackAssociation/interface/ParametersDefinerForTP.h
Expand Up @@ -62,10 +62,21 @@ class ParametersDefinerForTP {
return vertex(iEvent, iSetup, tp.charge(), tp.vertex(), tp.p4());
}

virtual std::tuple<TrackingParticle::Vector, TrackingParticle::Point> momentumAndVertex(
const edm::Event &iEvent, const edm::EventSetup &iSetup, const TrackingParticleRef &tpr) const {
return momentumAndVertex(iEvent, iSetup, tpr->charge(), tpr->vertex(), tpr->p4());
}

std::tuple<TrackingParticle::Vector, TrackingParticle::Point> momentumAndVertex(const edm::Event &iEvent,
const edm::EventSetup &iSetup,
const Charge ch,
const Point &vtx,
const LorentzVector &lv) const;

virtual void initEvent(edm::Handle<SimHitTPAssociationProducer::SimHitTPAssociationList> simHitsTPAssocToSet) {}

virtual std::unique_ptr<ParametersDefinerForTP> clone() const {
return std::unique_ptr<ParametersDefinerForTP>(new ParametersDefinerForTP(*this));
return std::make_unique<ParametersDefinerForTP>(*this);
}

edm::InputTag beamSpotInputTag_;
Expand Down
36 changes: 36 additions & 0 deletions SimTracker/TrackAssociation/src/ParametersDefinerForTP.cc
Expand Up @@ -83,4 +83,40 @@ TrackingParticle::Point ParametersDefinerForTP::vertex(const edm::Event &iEvent,
return vertex;
}

std::tuple<TrackingParticle::Vector, TrackingParticle::Point> ParametersDefinerForTP::momentumAndVertex(
const edm::Event &iEvent,
const edm::EventSetup &iSetup,
const Charge charge,
const Point &vtx,
const LorentzVector &lv) const {
using namespace edm;

edm::ESHandle<MagneticField> theMF;
iSetup.get<IdealMagneticFieldRecord>().get(theMF);

edm::Handle<reco::BeamSpot> bs;
iEvent.getByLabel(beamSpotInputTag_, bs);

TrackingParticle::Point vertex(bs->x0(), bs->y0(), bs->z0());
TrackingParticle::Vector momentum(0, 0, 0);

FreeTrajectoryState ftsAtProduction(GlobalPoint(vtx.x(), vtx.y(), vtx.z()),
GlobalVector(lv.x(), lv.y(), lv.z()),
TrackCharge(charge),
theMF.product());

TSCBLBuilderNoMaterial tscblBuilder;
TrajectoryStateClosestToBeamLine tsAtClosestApproach =
tscblBuilder(ftsAtProduction, *bs); // as in TrackProducerAlgorithm
if (tsAtClosestApproach.isValid()) {
GlobalPoint v = tsAtClosestApproach.trackStateAtPCA().position();
vertex = TrackingParticle::Point(v.x(), v.y(), v.z());
GlobalVector p = tsAtClosestApproach.trackStateAtPCA().momentum();
momentum = TrackingParticle::Vector(p.x(), p.y(), p.z());
;
}

return std::make_tuple(momentum, vertex);
}

TYPELOOKUP_DATA_REG(ParametersDefinerForTP);
Expand Up @@ -22,30 +22,47 @@
#include "DataFormats/Phase2TrackerCluster/interface/Phase2TrackerCluster1D.h"

#include "SimDataFormats/Track/interface/SimTrackContainer.h"
#include "SimDataFormats/Track/interface/UniqueSimTrackId.h"
#include "SimDataFormats/TrackerDigiSimLink/interface/StripDigiSimLink.h"
#include "SimDataFormats/TrackerDigiSimLink/interface/PixelDigiSimLink.h"
#include "DataFormats/Phase2TrackerDigi/interface/Phase2TrackerDigi.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
#include "SimTracker/TrackerHitAssociation/interface/ClusterTPAssociation.h"

class ClusterTPAssociationProducer : public edm::global::EDProducer<> {
namespace {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@VinInn , sorry I may be miss something but the same definition is in the header and here in the source.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

which definition?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I did not modify any header

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am sorry, the correct question: the same code included both inside SimHitTPAssociationProducer.cc and ClusterTPAssociationProducer.cc

namespace {
using TkId = std::pair<uint32_t, EncodedEventId>;
struct TkIdHash {
std::size_t operator()(TkId const& s) const noexcept {
std::size_t h1 = std::hash<uint32_t>{}(s.first);
std::size_t h2 = std::hash<uint32_t>{}(s.second.rawId());
return h1 ^ (h2 << 1);
}
};
} // namespace

I understand , that it is a correct and the most performant variant but Is it unavoidable to repeat these lines in each producer? For curiosity, how much you win?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes, can be factorized somewhere. (Where?)
What I win is that I will not need to recompile the whole universe (if I include it in EncodedEventId for instance) and I do not pollute the global namespace.
Actually I do not know where to locate it .... Any suggestion? One has to avoid additional dependencies (and of course the name shall be changed to be more unique etc)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Given the dependence chain (SimTrack::trackId(), EncodedEventId), a new header in SimDataFormats/Track/interface could be somewhat natural place.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Any suggestion for the name of the alias of the pair and for the hash struct (and the file itself)?
I hope that alias and struct can be located in the same file....

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just thinking out loud , maybe something along SimTrackAndEventId, SimTrackAndEncodedEventId, SimTrackEventIdPair (all with ...Hash)?

I'm not super happy with any of them though (and more often than not I find good names difficult to come up).

I hope that alias and struct can be located in the same file....

Given the usage I don't see strong motivations to use different files.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What about: UniqueSimTrackId (or SimTrackUniqueId) ?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What about: UniqueSimTrackId (or SimTrackUniqueId) ?

I would be fine with either, UniqueSimTrackId sounds maybe a bit better.


template <typename T>
void getSimTrackId(std::vector<UniqueSimTrackId>& simTrkId,
const edm::Handle<edm::DetSetVector<T> >& simLinks,
const DetId& detId,
uint32_t channel) {
auto isearch = simLinks->find(detId);
if (isearch != simLinks->end()) {
// Loop over DigiSimLink in this det unit
edm::DetSet<T> link_detset = (*isearch);
for (typename edm::DetSet<T>::const_iterator it = link_detset.data.begin(); it != link_detset.data.end(); ++it) {
if (channel == it->channel()) {
simTrkId.emplace_back(it->SimTrackId(), it->eventId());
}
}
}
}

} // namespace

class ClusterTPAssociationProducer final : public edm::global::EDProducer<> {
public:
typedef std::vector<OmniClusterRef> OmniClusterCollection;
using OmniClusterCollection = std::vector<OmniClusterRef>;

explicit ClusterTPAssociationProducer(const edm::ParameterSet&);
~ClusterTPAssociationProducer() override;
~ClusterTPAssociationProducer() override = default;

static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

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

template <typename T>
std::vector<std::pair<uint32_t, EncodedEventId> > getSimTrackId(const edm::Handle<edm::DetSetVector<T> >& simLinks,
const DetId& detId,
uint32_t channel) const;

edm::EDGetTokenT<edm::DetSetVector<PixelDigiSimLink> > sipixelSimLinksToken_;
edm::EDGetTokenT<edm::DetSetVector<StripDigiSimLink> > sistripSimLinksToken_;
edm::EDGetTokenT<edm::DetSetVector<PixelDigiSimLink> > siphase2OTSimLinksToken_;
Expand Down Expand Up @@ -73,8 +90,6 @@ ClusterTPAssociationProducer::ClusterTPAssociationProducer(const edm::ParameterS
produces<ClusterTPAssociation>();
}

ClusterTPAssociationProducer::~ClusterTPAssociationProducer() {}

void ClusterTPAssociationProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
edm::ParameterSetDescription desc;
desc.add<edm::InputTag>("simTrackSrc", edm::InputTag("g4SimHits"));
Expand Down Expand Up @@ -121,22 +136,23 @@ void ClusterTPAssociationProducer::produce(edm::StreamID, edm::Event& iEvent, co
auto clusterTPList = std::make_unique<ClusterTPAssociation>(TPCollectionH);

// prepare temporary map between SimTrackId and TrackingParticle index
std::map<std::pair<size_t, EncodedEventId>, TrackingParticleRef> mapping;
for (TrackingParticleCollection::size_type itp = 0; itp < TPCollectionH.product()->size(); ++itp) {
TrackingParticleRef trackingParticle(TPCollectionH, itp);

std::unordered_map<UniqueSimTrackId, TrackingParticleRef, UniqueSimTrackIdHash> mapping;
auto const& tpColl = *TPCollectionH.product();
for (TrackingParticleCollection::size_type itp = 0; itp < tpColl.size(); ++itp) {
TrackingParticleRef trackingParticleRef(TPCollectionH, itp);
auto const& trackingParticle = tpColl[itp];
// SimTracks inside TrackingParticle
EncodedEventId eid(trackingParticle->eventId());
EncodedEventId eid(trackingParticle.eventId());
//size_t index = 0;
for (std::vector<SimTrack>::const_iterator itrk = trackingParticle->g4Track_begin();
itrk != trackingParticle->g4Track_end();
++itrk) {
std::pair<uint32_t, EncodedEventId> trkid(itrk->trackId(), eid);
for (auto const& trk : trackingParticle.g4Tracks()) {
UniqueSimTrackId trkid(trk.trackId(), eid);
//std::cout << "creating map for id: " << trkid.first << " with tp: " << trackingParticle.key() << std::endl;
mapping.insert(std::make_pair(trkid, trackingParticle));
mapping.insert(std::make_pair(trkid, trackingParticleRef));
}
}

std::unordered_set<UniqueSimTrackId, UniqueSimTrackIdHash> simTkIds;
std::vector<UniqueSimTrackId> trkid;
if (foundPixelClusters) {
// Pixel Clusters
clusterTPList->addKeyID(pixelClusters.id());
Expand All @@ -150,20 +166,16 @@ void ClusterTPAssociationProducer::produce(edm::StreamID, edm::Event& iEvent, co
const SiPixelCluster& cluster = (*di);
edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> c_ref = edmNew::makeRefTo(pixelClusters, di);

std::set<std::pair<uint32_t, EncodedEventId> > simTkIds;
simTkIds.clear();
for (int irow = cluster.minPixelRow(); irow <= cluster.maxPixelRow(); ++irow) {
for (int icol = cluster.minPixelCol(); icol <= cluster.maxPixelCol(); ++icol) {
uint32_t channel = PixelChannelIdentifier::pixelToChannel(irow, icol);
std::vector<std::pair<uint32_t, EncodedEventId> > trkid(
getSimTrackId<PixelDigiSimLink>(sipixelSimLinks, detId, channel));
if (trkid.empty())
continue;
trkid.clear();
getSimTrackId<PixelDigiSimLink>(trkid, sipixelSimLinks, detId, channel);
simTkIds.insert(trkid.begin(), trkid.end());
}
}
for (std::set<std::pair<uint32_t, EncodedEventId> >::const_iterator iset = simTkIds.begin();
iset != simTkIds.end();
iset++) {
for (auto iset = simTkIds.begin(); iset != simTkIds.end(); iset++) {
auto ipos = mapping.find(*iset);
if (ipos != mapping.end()) {
//std::cout << "cluster in detid: " << detid << " from tp: " << ipos->second.key() << " " << iset->first << std::endl;
Expand All @@ -190,20 +202,16 @@ void ClusterTPAssociationProducer::produce(edm::StreamID, edm::Event& iEvent, co
const SiStripCluster& cluster = (*di);
edm::Ref<edmNew::DetSetVector<SiStripCluster>, SiStripCluster> c_ref = edmNew::makeRefTo(stripClusters, di);

std::set<std::pair<uint32_t, EncodedEventId> > simTkIds;
simTkIds.clear();
int first = cluster.firstStrip();
int last = first + cluster.amplitudes().size();

for (int istr = first; istr < last; ++istr) {
std::vector<std::pair<uint32_t, EncodedEventId> > trkid(
getSimTrackId<StripDigiSimLink>(sistripSimLinks, detId, istr));
if (trkid.empty())
continue;
trkid.clear();
getSimTrackId<StripDigiSimLink>(trkid, sistripSimLinks, detId, istr);
simTkIds.insert(trkid.begin(), trkid.end());
}
for (std::set<std::pair<uint32_t, EncodedEventId> >::const_iterator iset = simTkIds.begin();
iset != simTkIds.end();
iset++) {
for (auto iset = simTkIds.begin(); iset != simTkIds.end(); iset++) {
auto ipos = mapping.find(*iset);
if (ipos != mapping.end()) {
//std::cout << "cluster in detid: " << detid << " from tp: " << ipos->second.key() << " " << iset->first << std::endl;
Expand Down Expand Up @@ -233,20 +241,16 @@ void ClusterTPAssociationProducer::produce(edm::StreamID, edm::Event& iEvent, co
edm::Ref<edmNew::DetSetVector<Phase2TrackerCluster1D>, Phase2TrackerCluster1D> c_ref =
edmNew::makeRefTo(phase2OTClusters, di);

std::set<std::pair<uint32_t, EncodedEventId> > simTkIds;
simTkIds.clear();

for (unsigned int istr(0); istr < cluster.size(); ++istr) {
uint32_t channel = Phase2TrackerDigi::pixelToChannel(cluster.firstRow() + istr, cluster.column());
std::vector<std::pair<uint32_t, EncodedEventId> > trkid(
getSimTrackId<PixelDigiSimLink>(siphase2OTSimLinks, detId, channel));
if (trkid.empty())
continue;
trkid.clear();
getSimTrackId<PixelDigiSimLink>(trkid, siphase2OTSimLinks, detId, channel);
simTkIds.insert(trkid.begin(), trkid.end());
}

for (std::set<std::pair<uint32_t, EncodedEventId> >::const_iterator iset = simTkIds.begin();
iset != simTkIds.end();
iset++) {
for (auto iset = simTkIds.begin(); iset != simTkIds.end(); iset++) {
auto ipos = mapping.find(*iset);
if (ipos != mapping.end()) {
clusterTPList->emplace_back(OmniClusterRef(c_ref), ipos->second);
Expand All @@ -260,26 +264,6 @@ void ClusterTPAssociationProducer::produce(edm::StreamID, edm::Event& iEvent, co
iEvent.put(std::move(clusterTPList));
}

template <typename T>
std::vector<std::pair<uint32_t, EncodedEventId> >
//std::pair<uint32_t, EncodedEventId>
ClusterTPAssociationProducer::getSimTrackId(const edm::Handle<edm::DetSetVector<T> >& simLinks,
const DetId& detId,
uint32_t channel) const {
//std::pair<uint32_t, EncodedEventId> simTrkId;
std::vector<std::pair<uint32_t, EncodedEventId> > simTrkId;
auto isearch = simLinks->find(detId);
if (isearch != simLinks->end()) {
// Loop over DigiSimLink in this det unit
edm::DetSet<T> link_detset = (*isearch);
for (typename edm::DetSet<T>::const_iterator it = link_detset.data.begin(); it != link_detset.data.end(); ++it) {
if (channel == it->channel()) {
simTrkId.push_back(std::make_pair(it->SimTrackId(), it->eventId()));
}
}
}
return simTrkId;
}
#include "FWCore/PluginManager/interface/ModuleDef.h"
#include "FWCore/Framework/interface/MakerMacros.h"

Expand Down