Skip to content

Commit

Permalink
Merge pull request cms-sw#10 from Abhirikshma/geometric_linking_CMSSW…
Browse files Browse the repository at this point in the history 
…_12_4_0_pre2

Geometric linking CMSSW_12_4_0_pre2
  • Loading branch information
@felicepantaleo
felicepantaleo authored Apr 27, 2022
2 parents f4269d7 + dc6bcc6 commit 82bf232
Show file tree
Hide file tree
Showing 4 changed files with 86 additions and 24 deletions.
6 changes: 5 additions & 1 deletion RecoHGCal/TICL/plugins/LinkingAlgoBase.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,11 @@ namespace ticl {
const edm::ESHandle<Propagator> propH) = 0;

virtual void linkTracksters(const edm::Handle<std::vector<reco::Track>> tkH,
const edm::Handle<std::vector<reco::Muon>> muH,
const edm::ValueMap<float> &tkTime,
const edm::ValueMap<float> &tkTimeErr,
const edm::ValueMap<float> &tkTimeQual,
const double trackTimeQualThreshold,
const std::vector<reco::Muon> &muons,
const StringCutObjectSelector<reco::Track> cutTk,
const edm::Handle<std::vector<Trackster>> tsH,
std::vector<TICLCandidate>& resultTracksters) = 0;
Expand Down
71 changes: 50 additions & 21 deletions RecoHGCal/TICL/plugins/LinkingAlgoByPCAGeometric.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,11 @@ void LinkingAlgoByPCAGeometric::buildLayers() {
}

void LinkingAlgoByPCAGeometric::linkTracksters(const edm::Handle<std::vector<reco::Track>> tkH,
const edm::Handle<std::vector<reco::Muon>> muH,
const edm::ValueMap<float> &tkTime,
const edm::ValueMap<float> &tkTimeErr,
const edm::ValueMap<float> &tkTimeQual,
const double tkTimeQualThreshold,
const std::vector<reco::Muon> &muons,
const StringCutObjectSelector<reco::Track> cutTk,
const edm::Handle<std::vector<Trackster>> tsH,
std::vector<TICLCandidate> &resultLinked) {
Expand All @@ -108,8 +112,7 @@ void LinkingAlgoByPCAGeometric::linkTracksters(const edm::Handle<std::vector<rec

const auto &tracks = *tkH;
const auto &tracksters = *tsH;
const auto &muons = *muH;
std::cout << "N mu: " << muons.size() << std::endl;

auto bFieldProd = bfield_.product();
const Propagator &prop = (*propagator_);

Expand Down Expand Up @@ -146,7 +149,7 @@ void LinkingAlgoByPCAGeometric::linkTracksters(const edm::Handle<std::vector<rec
reco::TrackRef trackref = reco::TrackRef(tkH, i);
// also veto tracks associated to muons
int muId = PFMuonAlgo::muAssocToTrack(trackref, muons);
std::cout << "track (eta)" << i << " (" << tk.eta() <<")" << " muid " << muId << std::endl;
//std::cout << "track (eta)" << i << " (" << tk.eta() <<") time " << tkTime[reco::TrackRef(tkH, i)] << " time qual " << tkTimeQual[reco::TrackRef(tkH, i)] << " muid " << muId << std::endl;
if (!cutTk((tk)) or muId != -1) {
continue;
}
Expand Down Expand Up @@ -176,6 +179,7 @@ void LinkingAlgoByPCAGeometric::linkTracksters(const edm::Handle<std::vector<rec
// Propagate tracksters
for (unsigned i = 0; i < tracksters.size(); ++i) {
const auto &t = tracksters[i];
//std::cout << "trackster " << i << " time " << t.time() << " energy " << t.raw_energy() << std::endl;
Vector directnv = t.eigenvectors(0);

if (abs(directnv.Z()) < 0.00001)
Expand Down Expand Up @@ -422,69 +426,94 @@ void LinkingAlgoByPCAGeometric::linkTracksters(const edm::Handle<std::vector<rec
}

TICLCandidate chargedCandidate;
double total_raw_pt = 0.;
double total_raw_energy = 0.;
auto energyCompatible = [&](const Trackster & ts, const reco::Track tk) -> bool {
// compatible if accumulated energy does not
// exceed track momentum by more than threshold
double threshold = std::min(0.2*ts.raw_energy(), 10.0);
if (!(total_raw_energy + ts.raw_energy() < tk.p() + threshold))
std::cout << "track p : " << tk.p() << " trackster energy : " << ts.raw_energy() << std::endl;
return (total_raw_energy + ts.raw_energy() < tk.p() + threshold);
};
auto timeCompatible = [&](const Trackster & ts, const reco::TrackRef tk) -> bool {
// compatible if trackster time is within 3sigma of
// track time; compatible if either: no time assigned
// to trackster or track time quality is below threshold
double maxDeltaT = 3.;
double tsT = ts.time();
double tsTErr = ts.timeError();
double tkT = tkTime[tk];
double tkTErr = tkTimeErr[tk];

if (tsT == -99. or tkTimeQual[tk] < tkTimeQualThreshold) return true;
if (!(std::abs(tsT - tkT) < maxDeltaT * sqrt(tsTErr * tsTErr + tkTErr * tkTErr)))
std::cout << "track time : " << tkT << " trackster time : " << tsT << std::endl;
return (std::abs(tsT - tkT) < maxDeltaT * sqrt(tsTErr * tsTErr + tkTErr * tkTErr));
};
auto tkRef = reco::TrackRef(tkH, i);

for (const unsigned ts3_idx : tracksters_near[i]) { // tk -> ts
if (!chargedMask[ts3_idx]) {
if ((total_raw_pt += tsH->at(ts3_idx).raw_pt()) > tracks[i].pt()) continue;
if (!energyCompatible(tracksters[ts3_idx], tracks[i]) or !timeCompatible(tracksters[ts3_idx], tkRef)) continue;
chargedCandidate.addTrackster(edm::Ptr<Trackster>(tsH, ts3_idx));
chargedMask[ts3_idx] = 1;
total_raw_pt += tsH->at(ts3_idx).raw_pt();
total_raw_energy += tracksters[ts3_idx].raw_energy();
}
for (const unsigned ts2_idx : tsNearAtInt[ts3_idx]) { // ts_EM -> ts_HAD
if (!chargedMask[ts2_idx]) {
if ((total_raw_pt += tsH->at(ts2_idx).raw_pt()) > tracks[i].pt()) continue;
if (!energyCompatible(tracksters[ts2_idx], tracks[i]) or !timeCompatible(tracksters[ts2_idx], tkRef)) continue;
chargedCandidate.addTrackster(edm::Ptr<Trackster>(tsH, ts2_idx));
chargedMask[ts2_idx] = 1;
total_raw_pt += tsH->at(ts2_idx).raw_pt();
total_raw_energy += tracksters[ts2_idx].raw_energy();
}
for (const unsigned ts1_idx : tsHadNearAtInt[ts2_idx]) { // ts_HAD -> ts_HAD
if (!chargedMask[ts1_idx]) {
if ((total_raw_pt += tsH->at(ts1_idx).raw_pt()) > tracks[i].pt()) continue;
if (!energyCompatible(tracksters[ts1_idx], tracks[i]) or !timeCompatible(tracksters[ts1_idx], tkRef)) continue;
chargedCandidate.addTrackster(edm::Ptr<Trackster>(tsH, ts1_idx));
chargedMask[ts1_idx] = 1;
total_raw_pt += tsH->at(ts1_idx).raw_pt();
total_raw_energy += tracksters[ts1_idx].raw_energy();
}
}
}
for (const unsigned ts1_idx : tsHadNearAtInt[ts3_idx]) { // ts_HAD -> ts_HAD
if (!chargedMask[ts1_idx]) {
if ((total_raw_pt += tsH->at(ts1_idx).raw_pt()) > tracks[i].pt()) continue;
if (!energyCompatible(tracksters[ts1_idx], tracks[i]) or !timeCompatible(tracksters[ts1_idx], tkRef)) continue;
chargedCandidate.addTrackster(edm::Ptr<Trackster>(tsH, ts1_idx));
chargedMask[ts1_idx] = 1;
total_raw_pt += tsH->at(ts1_idx).raw_pt();
total_raw_energy += tracksters[ts1_idx].raw_energy();
}
}
}

for (const unsigned ts4_idx : tsNearTkAtInt[i]) { // do the same for tk -> ts links at the interface
if (!chargedMask[ts4_idx]) {
if ((total_raw_pt += tsH->at(ts4_idx).raw_pt()) > tracks[i].pt()) continue;
if (!energyCompatible(tracksters[ts4_idx], tracks[i]) or !timeCompatible(tracksters[ts4_idx], tkRef)) continue;
chargedCandidate.addTrackster(edm::Ptr<Trackster>(tsH, ts4_idx));
chargedMask[ts4_idx] = 1;
total_raw_pt += tsH->at(ts4_idx).raw_pt();
total_raw_energy += tracksters[ts4_idx].raw_energy();
}
for (const unsigned ts2_idx : tsNearAtInt[ts4_idx]) {
if (!chargedMask[ts2_idx]) {
if ((total_raw_pt += tsH->at(ts2_idx).raw_pt()) > tracks[i].pt()) continue;
if (!energyCompatible(tracksters[ts2_idx], tracks[i]) or !timeCompatible(tracksters[ts2_idx], tkRef)) continue;
chargedCandidate.addTrackster(edm::Ptr<Trackster>(tsH, ts2_idx));
chargedMask[ts2_idx] = 1;
total_raw_pt += tsH->at(ts2_idx).raw_pt();
total_raw_energy += tracksters[ts2_idx].raw_energy();
}
for (const unsigned ts1_idx : tsHadNearAtInt[ts2_idx]) {
if (!chargedMask[ts1_idx]) {
if ((total_raw_pt += tsH->at(ts1_idx).raw_pt()) > tracks[i].pt()) continue;
if (!energyCompatible(tracksters[ts1_idx], tracks[i]) or !timeCompatible(tracksters[ts1_idx], tkRef)) continue;
chargedCandidate.addTrackster(edm::Ptr<Trackster>(tsH, ts1_idx));
chargedMask[ts1_idx] = 1;
total_raw_pt += tsH->at(ts1_idx).raw_pt();
total_raw_energy += tracksters[ts1_idx].raw_energy();
}
}
}
for (const unsigned ts1_idx : tsHadNearAtInt[ts4_idx]) {
if (!chargedMask[ts1_idx]) {
if ((total_raw_pt += tsH->at(ts1_idx).raw_pt()) > tracks[i].pt()) continue;
if (!energyCompatible(tracksters[ts1_idx], tracks[i]) or !timeCompatible(tracksters[ts1_idx], tkRef)) continue;
chargedCandidate.addTrackster(edm::Ptr<Trackster>(tsH, ts1_idx));
chargedMask[ts1_idx] = 1;
total_raw_pt += tsH->at(ts1_idx).raw_pt();
total_raw_energy += tracksters[ts1_idx].raw_energy();
}
}
}
Expand Down
6 changes: 5 additions & 1 deletion RecoHGCal/TICL/plugins/LinkingAlgoByPCAGeometric.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,11 @@ namespace ticl {
const edm::ESHandle<Propagator> propH) override;

void linkTracksters(const edm::Handle<std::vector<reco::Track>>,
const edm::Handle<std::vector<reco::Muon>>,
const edm::ValueMap<float>&,
const edm::ValueMap<float>&,
const edm::ValueMap<float>&,
const double,
const std::vector<reco::Muon>&,
const StringCutObjectSelector<reco::Track>,
const edm::Handle<std::vector<Trackster>>,
std::vector<TICLCandidate> &) override;
Expand Down
27 changes: 26 additions & 1 deletion RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,6 +101,9 @@ class TrackstersMergeProducer : public edm::stream::EDProducer<> {
const edm::EDGetTokenT<std::vector<reco::CaloCluster>> clusters_token_;
const edm::EDGetTokenT<edm::ValueMap<std::pair<float, float>>> clustersTime_token_;
const edm::EDGetTokenT<std::vector<reco::Track>> tracks_token_;
const edm::EDGetTokenT<edm::ValueMap<float>> tracks_time_token_;
const edm::EDGetTokenT<edm::ValueMap<float>> tracks_time_quality_token_;
const edm::EDGetTokenT<edm::ValueMap<float>> tracks_time_err_token_;
const edm::EDGetTokenT<std::vector<reco::Muon>> muons_token_;
const std::string tfDnnLabel_;
const edm::ESGetToken<TfGraphDefWrapper, TfGraphRecord> tfDnnToken_;
Expand All @@ -123,6 +126,7 @@ class TrackstersMergeProducer : public edm::stream::EDProducer<> {
const double track_min_eta_;
const double track_max_eta_;
const int track_max_missing_outerhits_;
const double trackTimeQualThreshold_;
const double cosangle_align_;
const double e_over_h_threshold_;
const double pt_neutral_threshold_;
Expand Down Expand Up @@ -166,6 +170,9 @@ TrackstersMergeProducer::TrackstersMergeProducer(const edm::ParameterSet &ps)
clustersTime_token_(
consumes<edm::ValueMap<std::pair<float, float>>>(ps.getParameter<edm::InputTag>("layer_clustersTime"))),
tracks_token_(consumes<std::vector<reco::Track>>(ps.getParameter<edm::InputTag>("tracks"))),
tracks_time_token_(consumes<edm::ValueMap<float>>(ps.getParameter<edm::InputTag>("tracksTime"))),
tracks_time_quality_token_(consumes<edm::ValueMap<float>>(ps.getParameter<edm::InputTag>("tracksTimeQual"))),
tracks_time_err_token_(consumes<edm::ValueMap<float>>(ps.getParameter<edm::InputTag>("tracksTimeErr"))),
muons_token_(consumes<std::vector<reco::Muon>>(ps.getParameter<edm::InputTag>("muons"))),
tfDnnLabel_(ps.getParameter<std::string>("tfDnnLabel")),
tfDnnToken_(esConsumes(edm::ESInputTag("", tfDnnLabel_))),
Expand All @@ -187,6 +194,7 @@ TrackstersMergeProducer::TrackstersMergeProducer(const edm::ParameterSet &ps)
track_min_eta_(ps.getParameter<double>("track_min_eta")),
track_max_eta_(ps.getParameter<double>("track_max_eta")),
track_max_missing_outerhits_(ps.getParameter<int>("track_max_missing_outerhits")),
trackTimeQualThreshold_(ps.getParameter<double>("timingQualityThreshold")),
cosangle_align_(ps.getParameter<double>("cosangle_align")),
e_over_h_threshold_(ps.getParameter<double>("e_over_h_threshold")),
pt_neutral_threshold_(ps.getParameter<double>("pt_neutral_threshold")),
Expand Down Expand Up @@ -333,9 +341,22 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es

edm::Handle<std::vector<reco::Muon>> muons_h;
evt.getByToken(muons_token_, muons_h);
const auto &muons = *muons_h;

edm::Handle<edm::ValueMap<float>> trackTime_h;
evt.getByToken(tracks_time_token_, trackTime_h);
const auto &trackTime = *trackTime_h;

edm::Handle<edm::ValueMap<float>> trackTimeErr_h;
evt.getByToken(tracks_time_err_token_, trackTimeErr_h);
const auto &trackTimeErr = *trackTimeErr_h;

edm::Handle<edm::ValueMap<float>> trackTimeQual_h;
evt.getByToken(tracks_time_quality_token_, trackTimeQual_h);
const auto &trackTimeQual = *trackTimeQual_h;
// Linking
auto resultTrackstersLinked = std::make_unique<std::vector<TICLCandidate>>();
linkingAlgo_->linkTracksters(track_h, muons_h, cutTk_, trackstersclue3d_h, *resultTrackstersLinked);
linkingAlgo_->linkTracksters(track_h, trackTime, trackTimeErr, trackTimeQual, trackTimeQualThreshold_, muons, cutTk_, trackstersclue3d_h, *resultTrackstersLinked);

// Print debug info
if (debug_) {
Expand Down Expand Up @@ -954,6 +975,9 @@ void TrackstersMergeProducer::fillDescriptions(edm::ConfigurationDescriptions &d
desc.add<edm::InputTag>("layer_clusters", edm::InputTag("hgcalLayerClusters"));
desc.add<edm::InputTag>("layer_clustersTime", edm::InputTag("hgcalLayerClusters", "timeLayerCluster"));
desc.add<edm::InputTag>("tracks", edm::InputTag("generalTracks"));
desc.add<edm::InputTag>("tracksTime", edm::InputTag("tofPID:t0"));
desc.add<edm::InputTag>("tracksTimeQual", edm::InputTag("mtdTrackQualityMVA:mtdQualMVA"));
desc.add<edm::InputTag>("tracksTimeErr", edm::InputTag("tofPID:sigmat0"));
desc.add<edm::InputTag>("muons", edm::InputTag("muons1stStep"));
desc.add<std::string>("detector", "HGCAL");
desc.add<std::string>("propagator", "PropagatorWithMaterial");
Expand All @@ -971,6 +995,7 @@ void TrackstersMergeProducer::fillDescriptions(edm::ConfigurationDescriptions &d
desc.add<double>("track_min_eta", 1.48);
desc.add<double>("track_max_eta", 3.);
desc.add<int>("track_max_missing_outerhits", 5);
desc.add<double>("timingQualityThreshold", 0.5);
desc.add<double>("cosangle_align", 0.9945);
desc.add<double>("e_over_h_threshold", 1.);
desc.add<double>("pt_neutral_threshold", 2.);
Expand Down

0 comments on commit 82bf232

Please sign in to comment.