diff --git a/Configuration/ProcessModifiers/python/ticl_v3_cff.py b/Configuration/ProcessModifiers/python/ticl_v3_cff.py new file mode 100644 index 0000000000000..77f0f70c70404 --- /dev/null +++ b/Configuration/ProcessModifiers/python/ticl_v3_cff.py @@ -0,0 +1,5 @@ +import FWCore.ParameterSet.Config as cms + +# This modifier is for running TICL v3. + +ticl_v3 = cms.Modifier() diff --git a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py index 5b879ca1f4310..584230afffa2f 100644 --- a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py +++ b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py @@ -507,6 +507,30 @@ def condition(self, fragment, stepList, key, hasHarvest): upgradeWFs['ticl_FastJet'].step3 = {'--procModifiers': 'fastJetTICL'} upgradeWFs['ticl_FastJet'].step4 = {'--procModifiers': 'fastJetTICL'} +class UpgradeWorkflow_ticl_v3(UpgradeWorkflow): + def setup_(self, step, stepName, stepDict, k, properties): + if 'RecoGlobal' in step: + stepDict[stepName][k] = merge([self.step3, stepDict[step][k]]) + if 'HARVESTGlobal' in step: + stepDict[stepName][k] = merge([self.step4, stepDict[step][k]]) + def condition(self, fragment, stepList, key, hasHarvest): + return (fragment=="TTbar_14TeV" or 'CloseByP' in fragment or 'Eta1p7_2p7' in fragment) and '2026' in key +upgradeWFs['ticl_v3'] = UpgradeWorkflow_ticl_v3( + steps = [ + 'RecoGlobal', + 'HARVESTGlobal' + ], + PU = [ + 'RecoGlobal', + 'HARVESTGlobal' + ], + suffix = '_ticl_v3', + offset = 0.203, +) +upgradeWFs['ticl_v3'].step3 = {'--procModifiers': 'ticl_v3'} +upgradeWFs['ticl_v3'].step4 = {'--procModifiers': 'ticl_v3'} + + # Track DNN workflows class UpgradeWorkflow_trackdnn(UpgradeWorkflow): def setup_(self, step, stepName, stepDict, k, properties): diff --git a/DataFormats/HGCalReco/interface/Trackster.h b/DataFormats/HGCalReco/interface/Trackster.h index c5e7a067c4196..242b1080272dc 100644 --- a/DataFormats/HGCalReco/interface/Trackster.h +++ b/DataFormats/HGCalReco/interface/Trackster.h @@ -38,7 +38,7 @@ namespace ticl { Trackster() : iterationIndex_(0), - seedIndex_(0), + seedIndex_(-1), time_(0.f), timeError_(-1.f), regressed_energy_(0.f), diff --git a/HLTrigger/Configuration/python/HLT_75e33/modules/ticlTrackstersMerge_cfi.py b/HLTrigger/Configuration/python/HLT_75e33/modules/ticlTrackstersMerge_cfi.py index 00ac8ca9ade75..2bab985afc1cc 100644 --- a/HLTrigger/Configuration/python/HLT_75e33/modules/ticlTrackstersMerge_cfi.py +++ b/HLTrigger/Configuration/python/HLT_75e33/modules/ticlTrackstersMerge_cfi.py @@ -1,8 +1,7 @@ import FWCore.ParameterSet.Config as cms -ticlTrackstersMerge = cms.EDProducer("TrackstersMergeProducer", +ticlTrackstersMerge = cms.EDProducer("TrackstersMergeProducerV3", cosangle_align = cms.double(0.9945), - debug = cms.bool(True), e_over_h_threshold = cms.double(1), eid_input_name = cms.string('input'), eid_min_cluster_energy = cms.double(1), diff --git a/RecoHGCal/TICL/interface/PatternRecognitionAlgoBase.h b/RecoHGCal/TICL/interface/PatternRecognitionAlgoBase.h index c47aca80a7d92..507777de13d29 100644 --- a/RecoHGCal/TICL/interface/PatternRecognitionAlgoBase.h +++ b/RecoHGCal/TICL/interface/PatternRecognitionAlgoBase.h @@ -13,6 +13,7 @@ #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "DataFormats/Common/interface/ValueMap.h" #include "RecoHGCal/TICL/interface/GlobalCache.h" +#include "RecoHGCal/TICL/interface/commons.h" #include "FWCore/Framework/interface/ConsumesCollector.h" #include "PhysicsTools/TensorFlow/interface/TensorFlow.h" @@ -54,8 +55,6 @@ namespace ticl { std::vector& result, std::unordered_map>& seedToTracksterAssociation) = 0; - enum VerbosityLevel { None = 0, Basic, Advanced, Expert, Guru }; - protected: int algo_verbosity_; }; diff --git a/RecoHGCal/TICL/interface/commons.h b/RecoHGCal/TICL/interface/commons.h index ec52182bd49d6..1ec1956cee38a 100644 --- a/RecoHGCal/TICL/interface/commons.h +++ b/RecoHGCal/TICL/interface/commons.h @@ -7,6 +7,10 @@ namespace ticl { + //constants + constexpr double mpion = 0.13957; + constexpr float mpion2 = mpion * mpion; + inline Trackster::ParticleType tracksterParticleTypeFromPdgId(int pdgId, int charge) { if (pdgId == 111) { return Trackster::ParticleType::neutral_pion; @@ -33,42 +37,8 @@ namespace ticl { } } - static void addTrackster( - const int& index, - const std::vector, std::pair>>& lcVec, - const std::vector& inputClusterMask, - const float& fractionCut_, - const float& energy, - const int& pdgId, - const int& charge, - const edm::ProductID& seed, - const Trackster::IterationIndex iter, - std::vector& output_mask, - std::vector& result) { - if (lcVec.empty()) - return; - - Trackster tmpTrackster; - tmpTrackster.zeroProbabilities(); - tmpTrackster.vertices().reserve(lcVec.size()); - tmpTrackster.vertex_multiplicity().reserve(lcVec.size()); - for (auto const& [lc, energyScorePair] : lcVec) { - if (inputClusterMask[lc.index()] > 0) { - double fraction = energyScorePair.first / lc->energy(); - if (fraction < fractionCut_) - continue; - tmpTrackster.vertices().push_back(lc.index()); - output_mask[lc.index()] -= fraction; - tmpTrackster.vertex_multiplicity().push_back(1. / fraction); - } - } - - tmpTrackster.setIdProbability(tracksterParticleTypeFromPdgId(pdgId, charge), 1.f); - tmpTrackster.setRegressedEnergy(energy); - tmpTrackster.setIteration(iter); - tmpTrackster.setSeed(seed, index); - result.emplace_back(tmpTrackster); - } + // verbosity levels for ticl algorithms + enum VerbosityLevel { None = 0, Basic, Advanced, Expert, Guru }; } // namespace ticl diff --git a/RecoHGCal/TICL/plugins/HGCGraph.cc b/RecoHGCal/TICL/plugins/HGCGraph.cc index 964f8b530928e..948e7c2acae1e 100644 --- a/RecoHGCal/TICL/plugins/HGCGraph.cc +++ b/RecoHGCal/TICL/plugins/HGCGraph.cc @@ -71,7 +71,7 @@ void HGCGraphT::makeAndConnectDoublets(const TILES &histo, endEtaBin = std::min(entryEtaBin + etaWindow + 1, nEtaBins); startPhiBin = entryPhiBin - phiWindow; endPhiBin = entryPhiBin + phiWindow + 1; - if (verbosity_ > Guru) { + if (verbosity_ > ticl::VerbosityLevel::Guru) { LogDebug("HGCGraph") << " Entrance eta, phi: " << origin_eta << ", " << origin_phi << " entryEtaBin: " << entryEtaBin << " entryPhiBin: " << entryPhiBin << " globalBin: " << firstLayerHisto.globalBin(origin_eta, origin_phi) @@ -96,7 +96,7 @@ void HGCGraphT::makeAndConnectDoublets(const TILES &histo, : (il <= lastLayerFH) ? siblings_maxRSquared[1] : siblings_maxRSquared[2]; const int etaLimitIncreaseWindowBin = innerLayerHisto.etaBin(etaLimitIncreaseWindow); - if (verbosity_ > Advanced) { + if (verbosity_ > ticl::VerbosityLevel::Advanced) { LogDebug("HGCGraph") << "Limit of Eta for increase: " << etaLimitIncreaseWindow << " at etaBin: " << etaLimitIncreaseWindowBin << std::endl; } @@ -105,7 +105,7 @@ void HGCGraphT::makeAndConnectDoublets(const TILES &histo, auto offset = ieta * nPhiBins; for (int iphi_it = startPhiBin; iphi_it < endPhiBin; ++iphi_it) { int iphi = ((iphi_it % nPhiBins + nPhiBins) % nPhiBins); - if (verbosity_ > Guru) { + if (verbosity_ > ticl::VerbosityLevel::Advanced) { LogDebug("HGCGraph") << "Inner Global Bin: " << (offset + iphi) << " on layers I/O: " << currentInnerLayerId << "/" << currentOuterLayerId << " with clusters: " << innerLayerHisto[offset + iphi].size() << std::endl; @@ -113,7 +113,7 @@ void HGCGraphT::makeAndConnectDoublets(const TILES &histo, for (auto innerClusterId : innerLayerHisto[offset + iphi]) { // Skip masked clusters if (mask[innerClusterId] == 0.) { - if (verbosity_ > Advanced) + if (verbosity_ > ticl::VerbosityLevel::Advanced) LogDebug("HGCGraph") << "Skipping inner masked cluster " << innerClusterId << std::endl; continue; } @@ -128,7 +128,7 @@ void HGCGraphT::makeAndConnectDoublets(const TILES &histo, if (isGlobal && ieta > etaLimitIncreaseWindowBin) { etaWindow++; phiWindow++; - if (verbosity_ > Advanced) { + if (verbosity_ > ticl::VerbosityLevel::Advanced) { LogDebug("HGCGraph") << "Eta and Phi window increased by one" << std::endl; } } @@ -144,7 +144,7 @@ void HGCGraphT::makeAndConnectDoublets(const TILES &histo, // account for all other cases, since we add in // between a full nPhiBins slot. auto ophi = ((iphi + phiRange - phiWindow) % nPhiBins + nPhiBins) % nPhiBins; - if (verbosity_ > Guru) { + if (verbosity_ > ticl::VerbosityLevel::Guru) { LogDebug("HGCGraph") << "Outer Global Bin: " << (oeta * nPhiBins + ophi) << " on layers I/O: " << currentInnerLayerId << "/" << currentOuterLayerId << " with clusters: " << innerLayerHisto[oeta * nPhiBins + ophi].size() @@ -153,14 +153,14 @@ void HGCGraphT::makeAndConnectDoublets(const TILES &histo, for (auto outerClusterId : outerLayerHisto[oeta * nPhiBins + ophi]) { // Skip masked clusters if (mask[outerClusterId] == 0.) { - if (verbosity_ > Advanced) + if (verbosity_ > ticl::VerbosityLevel::Advanced) LogDebug("HGCGraph") << "Skipping outer masked cluster " << outerClusterId << std::endl; continue; } auto doubletId = allDoublets_.size(); if (maxDeltaTime != -1 && !areTimeCompatible(innerClusterId, outerClusterId, layerClustersTime, maxDeltaTime)) { - if (verbosity_ > Advanced) + if (verbosity_ > ticl::VerbosityLevel::Advanced) LogDebug("HGCGraph") << "Rejecting doublets due to timing!" << std::endl; continue; } @@ -175,7 +175,7 @@ void HGCGraphT::makeAndConnectDoublets(const TILES &histo, allDoublets_.emplace_back( innerClusterId, outerClusterId, doubletId, &layerClusters, r.index, false); } - if (verbosity_ > Advanced) { + if (verbosity_ > ticl::VerbosityLevel::Advanced) { LogDebug("HGCGraph") << "Creating doubletsId: " << doubletId << " layerLink in-out: [" << currentInnerLayerId << ", " << currentOuterLayerId << "] clusterLink in-out: [" << innerClusterId << ", " @@ -184,18 +184,19 @@ void HGCGraphT::makeAndConnectDoublets(const TILES &histo, isOuterClusterOfDoublets_[outerClusterId].push_back(doubletId); auto &neigDoublets = isOuterClusterOfDoublets_[innerClusterId]; auto &thisDoublet = allDoublets_[doubletId]; - if (verbosity_ > Expert) { + if (verbosity_ > ticl::VerbosityLevel::Expert) { LogDebug("HGCGraph") << "Checking compatibility of doubletId: " << doubletId << " with all possible inners doublets link by the innerClusterId: " << innerClusterId << std::endl; } - bool isRootDoublet = thisDoublet.checkCompatibilityAndTag(allDoublets_, - neigDoublets, - r.directionAtOrigin, - minCosTheta, - minCosPointing, - verbosity_ > Advanced); + bool isRootDoublet = + thisDoublet.checkCompatibilityAndTag(allDoublets_, + neigDoublets, + r.directionAtOrigin, + minCosTheta, + minCosPointing, + verbosity_ > ticl::VerbosityLevel::Advanced); if (isRootDoublet and checkDistanceRootDoubletVsSeed) { if (reco::deltaR2(layerClusters[innerClusterId].eta(), layerClusters[innerClusterId].phi(), @@ -217,7 +218,7 @@ void HGCGraphT::makeAndConnectDoublets(const TILES &histo, } } // #ifdef FP_DEBUG - if (verbosity_ > None) { + if (verbosity_ > ticl::VerbosityLevel::None) { LogDebug("HGCGraph") << "number of Root doublets " << theRootDoublets_.size() << " over a total number of doublets " << allDoublets_.size() << std::endl; } diff --git a/RecoHGCal/TICL/plugins/HGCGraph.h b/RecoHGCal/TICL/plugins/HGCGraph.h index bc127ba6c73d1..2b597db404fe4 100644 --- a/RecoHGCal/TICL/plugins/HGCGraph.h +++ b/RecoHGCal/TICL/plugins/HGCGraph.h @@ -59,7 +59,6 @@ class HGCGraphT { isOuterClusterOfDoublets_.shrink_to_fit(); } void setVerbosity(int level) { verbosity_ = level; } - enum VerbosityLevel { None = 0, Basic, Advanced, Expert, Guru }; private: std::vector allDoublets_; diff --git a/RecoHGCal/TICL/plugins/LinkingAlgoBase.h b/RecoHGCal/TICL/plugins/LinkingAlgoBase.h new file mode 100644 index 0000000000000..46340a42e567f --- /dev/null +++ b/RecoHGCal/TICL/plugins/LinkingAlgoBase.h @@ -0,0 +1,50 @@ +#ifndef RecoHGCal_TICL_LinkingAlgoBase_H__ +#define RecoHGCal_TICL_LinkingAlgoBase_H__ + +#include +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "DataFormats/HGCalReco/interface/Trackster.h" +#include "DataFormats/HGCalReco/interface/TICLCandidate.h" +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/TrackReco/interface/TrackFwd.h" +#include "DataFormats/MuonReco/interface/Muon.h" +#include "MagneticField/Engine/interface/MagneticField.h" +#include "TrackingTools/GeomPropagators/interface/Propagator.h" +#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h" +#include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h" + +namespace edm { + class Event; + class EventSetup; +} // namespace edm + +namespace ticl { + class LinkingAlgoBase { + public: + LinkingAlgoBase(const edm::ParameterSet& conf) : algo_verbosity_(conf.getParameter("algo_verbosity")) {} + + virtual ~LinkingAlgoBase(){}; + + virtual void initialize(const HGCalDDDConstants* hgcons, + const hgcal::RecHitTools rhtools, + const edm::ESHandle bfieldH, + const edm::ESHandle propH) = 0; + + virtual void linkTracksters(const edm::Handle> tkH, + const edm::ValueMap& tkTime, + const edm::ValueMap& tkTimeErr, + const edm::ValueMap& tkTimeQual, + const std::vector& muons, + const edm::Handle> tsH, + std::vector& resultTracksters) = 0; + + static void fillPSetDescription(edm::ParameterSetDescription& desc) { desc.add("algo_verbosity", 0); }; + + protected: + int algo_verbosity_; + }; +} // namespace ticl + +#endif diff --git a/RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.cc b/RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.cc new file mode 100644 index 0000000000000..b139e102cdfcd --- /dev/null +++ b/RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.cc @@ -0,0 +1,557 @@ +#include +#include +#include "RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.h" + +#include "DataFormats/GeometrySurface/interface/BoundDisk.h" +#include "DataFormats/HGCalReco/interface/Common.h" + +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h" + +#include "RecoParticleFlow/PFProducer/interface/PFMuonAlgo.h" + +using namespace ticl; + +LinkingAlgoByDirectionGeometric::LinkingAlgoByDirectionGeometric(const edm::ParameterSet &conf) + : LinkingAlgoBase(conf), + del_tk_ts_layer1_(conf.getParameter("delta_tk_ts_layer1")), + del_tk_ts_int_(conf.getParameter("delta_tk_ts_interface")), + del_ts_em_had_(conf.getParameter("delta_ts_em_had")), + del_ts_had_had_(conf.getParameter("delta_ts_had_had")), + timing_quality_threshold_(conf.getParameter("track_time_quality_threshold")), + pid_threshold_(conf.getParameter("pid_threshold")), + energy_em_over_total_threshold_(conf.getParameter("energy_em_over_total_threshold")), + filter_on_categories_(conf.getParameter>("filter_hadronic_on_categories")), + cutTk_(conf.getParameter("cutTk")) {} + +LinkingAlgoByDirectionGeometric::~LinkingAlgoByDirectionGeometric() {} + +void LinkingAlgoByDirectionGeometric::initialize(const HGCalDDDConstants *hgcons, + const hgcal::RecHitTools rhtools, + const edm::ESHandle bfieldH, + const edm::ESHandle propH) { + hgcons_ = hgcons; + rhtools_ = rhtools; + buildLayers(); + + bfield_ = bfieldH; + propagator_ = propH; +} + +math::XYZVector LinkingAlgoByDirectionGeometric::propagateTrackster(const Trackster &t, + const unsigned idx, + float zVal, + std::array &tracksterTiles) { + // needs only the positive Z co-ordinate of the surface to propagate to + // the correct sign is calculated inside according to the barycenter of trackster + Vector const &baryc = t.barycenter(); + Vector directnv = t.eigenvectors(0); + + // barycenter as direction for tracksters w/ poor PCA + // propagation still done to get the cartesian coords + // which are anyway converted to eta, phi in linking + // -> can be simplified later + + //FP: disable PCA propagation for the moment and fallback to barycenter position + // if (t.eigenvalues()[0] / t.eigenvalues()[1] < 20) + directnv = baryc.unit(); + + zVal *= (baryc.Z() > 0) ? 1 : -1; + + double par = (zVal - baryc.Z()) / directnv.Z(); + double xOnSurface = par * directnv.X() + baryc.X(); + double yOnSurface = par * directnv.Y() + baryc.Y(); + Vector tPoint(xOnSurface, yOnSurface, zVal); + if (tPoint.Eta() > 0) + tracksterTiles[1].fill(tPoint.Eta(), tPoint.Phi(), idx); + + else if (tPoint.Eta() < 0) + tracksterTiles[0].fill(tPoint.Eta(), tPoint.Phi(), idx); + + return tPoint; +} + +void LinkingAlgoByDirectionGeometric::findTrackstersInWindow( + const std::vector> &seedingCollection, + const std::array &tracksterTiles, + double delta, + unsigned trackstersSize, + std::vector> &resultCollection, + bool useMask = false) { + // Finds tracksters in tracksterTiles within an eta-phi window + // (given by delta) of the objects (track/trackster) in the seedingCollection. + // Element i in resultCollection is the vector of trackster + // indices found close to the i-th object in the seedingCollection. + // If specified, Tracksters are masked once found as close to an object. + std::vector mask(trackstersSize, 0); + + for (auto &i : seedingCollection) { + auto seed_eta = i.first.Eta(); + auto seed_phi = i.first.Phi(); + unsigned seedId = i.second; + auto sideZ = seed_eta > 0; //forward or backward region + const TICLLayerTile &tile = tracksterTiles[sideZ]; + double eta_min = std::max(abs(seed_eta) - delta, (double)TileConstants::minEta); + double eta_max = std::min(abs(seed_eta) + delta, (double)TileConstants::maxEta); + + std::array search_box = tile.searchBoxEtaPhi(eta_min, eta_max, seed_phi - delta, seed_phi + delta); + if (search_box[2] > search_box[3]) + search_box[3] += TileConstants::nPhiBins; + + for (int eta_i = search_box[0]; eta_i <= search_box[1]; ++eta_i) { + for (int phi_i = search_box[2]; phi_i <= search_box[3]; ++phi_i) { + const auto &in_box = tile[tile.globalBin(eta_i, (phi_i % TileConstants::nPhiBins))]; + for (const unsigned t_i : in_box) { + if (!mask[t_i]) { + resultCollection[seedId].push_back(t_i); + if (useMask) + mask[t_i] = 1; + } + } + } + } + } // seeding collection loop +} + +bool LinkingAlgoByDirectionGeometric::timeAndEnergyCompatible(double &total_raw_energy, + const reco::Track &track, + const Trackster &trackster, + const float &tkT, + const float &tkTErr, + const float &tkTimeQual) { + double threshold = std::min(0.2 * trackster.raw_energy(), 10.0); + + bool energyCompatible = (total_raw_energy + trackster.raw_energy() < track.p() + threshold); + // 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 tsT = trackster.time(); + double tsTErr = trackster.timeError(); + + bool timeCompatible = false; + + if (tsT == -99. or tkTimeQual < timing_quality_threshold_) + timeCompatible = true; + else { + timeCompatible = (std::abs(tsT - tkT) < maxDeltaT_ * sqrt(tsTErr * tsTErr + tkTErr * tkTErr)); + } + + if (LinkingAlgoBase::algo_verbosity_ > VerbosityLevel::Advanced) { + if (!(energyCompatible)) + LogDebug("LinkingAlgoByDirectionGeometric") + << "energy incompatible : track p " << track.p() << " trackster energy " << trackster.raw_energy() << "\n"; + if (!(timeCompatible)) + LogDebug("LinkingAlgoByDirectionGeometric") << "time incompatible : track time " << tkT << " +/- " << tkTErr + << " trackster time " << tsT << " +/- " << tsTErr << "\n"; + } + return energyCompatible && timeCompatible; +} + +void LinkingAlgoByDirectionGeometric::recordTrackster(const unsigned ts, //trackster index + const std::vector &tracksters, + const edm::Handle> tsH, + std::vector &ts_mask, + double &energy_in_candidate, + TICLCandidate &candidate) { + if (ts_mask[ts]) + return; + candidate.addTrackster(edm::Ptr(tsH, ts)); + ts_mask[ts] = 1; + energy_in_candidate += tracksters[ts].raw_energy(); +} + +void LinkingAlgoByDirectionGeometric::dumpLinksFound(std::vector> &resultCollection, + const char *label) const { +#ifdef EDM_ML_DEBUG + if (!(LinkingAlgoBase::algo_verbosity_ > VerbosityLevel::Advanced)) + return; + + LogDebug("LinkingAlgoByDirectionGeometric") << "All links found - " << label << "\n"; + LogDebug("LinkingAlgoByDirectionGeometric") << "(seed can either be a track or trackster depending on the step)\n"; + for (unsigned i = 0; i < resultCollection.size(); ++i) { + LogDebug("LinkingAlgoByDirectionGeometric") << "seed " << i << " - tracksters : "; + const auto &links = resultCollection[i]; + for (unsigned j = 0; j < links.size(); ++j) { + LogDebug("LinkingAlgoByDirectionGeometric") << j; + } + LogDebug("LinkingAlgoByDirectionGeometric") << "\n"; + } +#endif // EDM_ML_DEBUG +} + +void LinkingAlgoByDirectionGeometric::buildLayers() { + // build disks at HGCal front & EM-Had interface for track propagation + + float zVal = hgcons_->waferZ(1, true); + std::pair rMinMax = hgcons_->rangeR(zVal, true); + + float zVal_interface = rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z(); + std::pair rMinMax_interface = hgcons_->rangeR(zVal_interface, true); + + for (int iSide = 0; iSide < 2; ++iSide) { + float zSide = (iSide == 0) ? (-1. * zVal) : zVal; + firstDisk_[iSide] = + std::make_unique(Disk::build(Disk::PositionType(0, 0, zSide), + Disk::RotationType(), + SimpleDiskBounds(rMinMax.first, rMinMax.second, zSide - 0.5, zSide + 0.5)) + .get()); + + zSide = (iSide == 0) ? (-1. * zVal_interface) : zVal_interface; + interfaceDisk_[iSide] = std::make_unique( + Disk::build(Disk::PositionType(0, 0, zSide), + Disk::RotationType(), + SimpleDiskBounds(rMinMax_interface.first, rMinMax_interface.second, zSide - 0.5, zSide + 0.5)) + .get()); + } +} + +void LinkingAlgoByDirectionGeometric::linkTracksters(const edm::Handle> tkH, + const edm::ValueMap &tkTime, + const edm::ValueMap &tkTimeErr, + const edm::ValueMap &tkTimeQual, + const std::vector &muons, + const edm::Handle> tsH, + std::vector &resultLinked) { + const auto &tracks = *tkH; + const auto &tracksters = *tsH; + + auto bFieldProd = bfield_.product(); + const Propagator &prop = (*propagator_); + + // propagated point collections + // elements in the propagated points collecions are used + // to look for potential linkages in the appropriate tiles + std::vector> trackPColl; // propagated track points and index of track in collection + std::vector> tkPropIntColl; // tracks propagated to lastLayerEE + std::vector> tsPropIntColl; // Tracksters in CE-E, propagated to lastLayerEE + std::vector> tsHadPropIntColl; // Tracksters in CE-H, propagated to lastLayerEE + trackPColl.reserve(tracks.size()); + tkPropIntColl.reserve(tracks.size()); + tsPropIntColl.reserve(tracksters.size()); + tsHadPropIntColl.reserve(tracksters.size()); + // tiles, element 0 is bw, 1 is fw + std::array tracksterPropTiles = {}; // all Tracksters, propagated to layer 1 + std::array tsPropIntTiles = {}; // all Tracksters, propagated to lastLayerEE + std::array tsHadPropIntTiles = {}; // Tracksters in CE-H, propagated to lastLayerEE + + // filters (true for) anything but EM + auto isHadron = [&](const Trackster &t) -> bool { + auto cumulative_prob = 0.; + for (const auto index : filter_on_categories_) { + cumulative_prob += t.id_probabilities(index); + } + return ((cumulative_prob <= pid_threshold_) and (t.raw_em_energy() != t.raw_energy())) or + (t.raw_em_energy() < energy_em_over_total_threshold_ * t.raw_energy()); + }; + if (LinkingAlgoBase::algo_verbosity_ > VerbosityLevel::Advanced) + LogDebug("LinkingAlgoByDirectionGeometric") << "------- Geometric Linking ------- \n"; + + // Propagate tracks + std::vector candidateTrackIds; + candidateTrackIds.reserve(tracks.size()); + for (unsigned i = 0; i < tracks.size(); ++i) { + const auto &tk = tracks[i]; + + reco::TrackRef trackref = reco::TrackRef(tkH, i); + // also veto tracks associated to muons + int muId = PFMuonAlgo::muAssocToTrack(trackref, muons); + if (LinkingAlgoBase::algo_verbosity_ > VerbosityLevel::Advanced) + LogDebug("LinkingAlgoByDirectionGeometric") + << "track " << i << " - eta " << tk.eta() << " phi " << tk.phi() << " time " << tkTime[reco::TrackRef(tkH, i)] + << " time qual " << tkTimeQual[reco::TrackRef(tkH, i)] << " muid " << muId << "\n"; + if (!cutTk_((tk)) or muId != -1) + continue; + + // record tracks that can used to make a ticlcandidate + candidateTrackIds.push_back(i); + + // don't consider tracks below 2 GeV for linking + if (std::sqrt(tk.p() * tk.p() + ticl::mpion2) < tkEnergyCut_) + continue; + + int iSide = int(tk.eta() > 0); + const auto &fts = trajectoryStateTransform::outerFreeState((tk), bFieldProd); + // to the HGCal front + const auto &tsos = prop.propagate(fts, firstDisk_[iSide]->surface()); + if (tsos.isValid()) { + Vector trackP(tsos.globalPosition().x(), tsos.globalPosition().y(), tsos.globalPosition().z()); + trackPColl.emplace_back(trackP, i); + } + // to lastLayerEE + const auto &tsos_int = prop.propagate(fts, interfaceDisk_[iSide]->surface()); + if (tsos_int.isValid()) { + Vector trackP(tsos_int.globalPosition().x(), tsos_int.globalPosition().y(), tsos_int.globalPosition().z()); + tkPropIntColl.emplace_back(trackP, i); + } + } // Tracks + tkPropIntColl.shrink_to_fit(); + trackPColl.shrink_to_fit(); + candidateTrackIds.shrink_to_fit(); + // Propagate tracksters + for (unsigned i = 0; i < tracksters.size(); ++i) { + const auto &t = tracksters[i]; + if (LinkingAlgoBase::algo_verbosity_ > VerbosityLevel::Advanced) + LogDebug("LinkingAlgoByDirectionGeometric") + << "trackster " << i << " - eta " << t.barycenter().eta() << " phi " << t.barycenter().phi() << " time " + << t.time() << " energy " << t.raw_energy() << "\n"; + + // to HGCal front + float zVal = hgcons_->waferZ(1, true); + auto tsP = propagateTrackster(t, i, zVal, tracksterPropTiles); + + // to lastLayerEE + zVal = rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z(); + tsP = propagateTrackster(t, i, zVal, tsPropIntTiles); + + if (!isHadron(t)) // EM tracksters + tsPropIntColl.emplace_back(tsP, i); + else { // HAD + tsHadPropIntTiles[(t.barycenter().Z() > 0) ? 1 : 0].fill(tsP.Eta(), tsP.Phi(), i); + tsHadPropIntColl.emplace_back(tsP, i); + } + } // TS + tsPropIntColl.shrink_to_fit(); + tsHadPropIntColl.shrink_to_fit(); + // Track - Trackster link finding + // step 3: tracks -> all tracksters, at layer 1 + + std::vector> tsNearTk(tracks.size()); + findTrackstersInWindow(trackPColl, tracksterPropTiles, del_tk_ts_layer1_, tracksters.size(), tsNearTk); + + // step 4: tracks -> all tracksters, at lastLayerEE + + std::vector> tsNearTkAtInt(tracks.size()); + findTrackstersInWindow(tkPropIntColl, tsPropIntTiles, del_tk_ts_int_, tracksters.size(), tsNearTkAtInt); + + // Trackster - Trackster link finding + // step 2: tracksters EM -> HAD, at lastLayerEE + + std::vector> tsNearAtInt(tracksters.size()); + findTrackstersInWindow(tsPropIntColl, tsHadPropIntTiles, del_ts_em_had_, tracksters.size(), tsNearAtInt, true); + + // step 1: tracksters HAD -> HAD, at lastLayerEE + + std::vector> tsHadNearAtInt(tracksters.size()); + findTrackstersInWindow(tsHadPropIntColl, tsHadPropIntTiles, del_ts_had_had_, tracksters.size(), tsHadNearAtInt, true); +#ifdef EDM_ML_DEBUG + + dumpLinksFound(tsNearTk, "track -> tracksters at layer 1"); + dumpLinksFound(tsNearTkAtInt, "track -> tracksters at lastLayerEE"); + dumpLinksFound(tsNearAtInt, "EM -> HAD tracksters at lastLayerEE"); + dumpLinksFound(tsHadNearAtInt, "HAD -> HAD tracksters at lastLayerEE"); +#endif //EDM_ML_DEBUG + // make final collections + + std::vector chargedCandidates; + std::vector chargedHadronsFromTk; + std::vector chargedMask(tracksters.size(), 0); + for (unsigned &i : candidateTrackIds) { + if (tsNearTk[i].empty() && tsNearTkAtInt[i].empty()) { // nothing linked to track, make charged hadrons + TICLCandidate chargedHad; + const auto &tk = tracks[i]; + chargedHad.setCharge(tk.charge()); + chargedHad.setPdgId(211 * tk.charge()); + chargedHad.setTrackPtr(edm::Ptr(tkH, i)); + float energy = std::sqrt(tk.p() * tk.p() + ticl::mpion2); + chargedHad.setRawEnergy(energy); + math::PtEtaPhiMLorentzVector p4Polar(tk.pt(), tk.eta(), tk.phi(), ticl::mpion); + chargedHad.setP4(p4Polar); + chargedHadronsFromTk.push_back(chargedHad); + continue; + } + + TICLCandidate chargedCandidate; + double total_raw_energy = 0.; + + auto tkRef = reco::TrackRef(tkH, i); + auto track_time = tkTime[tkRef]; + auto track_timeErr = tkTimeErr[tkRef]; + auto track_timeQual = tkTimeQual[tkRef]; + + for (const unsigned ts3_idx : tsNearTk[i]) { // tk -> ts + if (timeAndEnergyCompatible( + total_raw_energy, tracks[i], tracksters[ts3_idx], track_time, track_timeErr, track_timeQual)) { + recordTrackster(ts3_idx, tracksters, tsH, chargedMask, total_raw_energy, chargedCandidate); + } + for (const unsigned ts2_idx : tsNearAtInt[ts3_idx]) { // ts_EM -> ts_HAD + if (timeAndEnergyCompatible( + total_raw_energy, tracks[i], tracksters[ts2_idx], track_time, track_timeErr, track_timeQual)) { + recordTrackster(ts2_idx, tracksters, tsH, chargedMask, total_raw_energy, chargedCandidate); + } + for (const unsigned ts1_idx : tsHadNearAtInt[ts2_idx]) { // ts_HAD -> ts_HAD + if (timeAndEnergyCompatible( + total_raw_energy, tracks[i], tracksters[ts1_idx], track_time, track_timeErr, track_timeQual)) { + recordTrackster(ts1_idx, tracksters, tsH, chargedMask, total_raw_energy, chargedCandidate); + } + } + } + for (const unsigned ts1_idx : tsHadNearAtInt[ts3_idx]) { // ts_HAD -> ts_HAD + if (timeAndEnergyCompatible( + total_raw_energy, tracks[i], tracksters[ts1_idx], track_time, track_timeErr, track_timeQual)) { + recordTrackster(ts1_idx, tracksters, tsH, chargedMask, total_raw_energy, chargedCandidate); + } + } + } + for (const unsigned ts4_idx : tsNearTkAtInt[i]) { // do the same for tk -> ts links at the interface + if (timeAndEnergyCompatible( + total_raw_energy, tracks[i], tracksters[ts4_idx], track_time, track_timeErr, track_timeQual)) { + recordTrackster(ts4_idx, tracksters, tsH, chargedMask, total_raw_energy, chargedCandidate); + } + for (const unsigned ts2_idx : tsNearAtInt[ts4_idx]) { + if (timeAndEnergyCompatible( + total_raw_energy, tracks[i], tracksters[ts2_idx], track_time, track_timeErr, track_timeQual)) { + recordTrackster(ts2_idx, tracksters, tsH, chargedMask, total_raw_energy, chargedCandidate); + } + for (const unsigned ts1_idx : tsHadNearAtInt[ts2_idx]) { + if (timeAndEnergyCompatible( + total_raw_energy, tracks[i], tracksters[ts1_idx], track_time, track_timeErr, track_timeQual)) { + recordTrackster(ts1_idx, tracksters, tsH, chargedMask, total_raw_energy, chargedCandidate); + } + } + } + for (const unsigned ts1_idx : tsHadNearAtInt[ts4_idx]) { + if (timeAndEnergyCompatible( + total_raw_energy, tracks[i], tracksters[ts1_idx], track_time, track_timeErr, track_timeQual)) { + recordTrackster(ts1_idx, tracksters, tsH, chargedMask, total_raw_energy, chargedCandidate); + } + } + } + + // do not create a candidate if no tracksters were added to candidate + // can happen if all the tracksters linked to that track were already masked + if (!chargedCandidate.tracksters().empty()) { + chargedCandidate.setTrackPtr(edm::Ptr(tkH, i)); + chargedCandidates.push_back(chargedCandidate); + } else { // create charged hadron + TICLCandidate chargedHad; + const auto &tk = tracks[i]; + chargedHad.setCharge(tk.charge()); + chargedHad.setPdgId(211 * tk.charge()); + chargedHad.setTrackPtr(edm::Ptr(tkH, i)); + float energy = std::sqrt(tk.p() * tk.p() + ticl::mpion2); + chargedHad.setRawEnergy(energy); + math::PtEtaPhiMLorentzVector p4Polar(tk.pt(), tk.eta(), tk.phi(), mpion); + chargedHad.setP4(p4Polar); + chargedHadronsFromTk.push_back(chargedHad); + } + } + + std::vector neutralCandidates; + std::vector neutralMask(tracksters.size(), 0); + for (unsigned i = 0; i < tracksters.size(); ++i) { + if (chargedMask[i]) + continue; + + TICLCandidate neutralCandidate; + if (tsNearAtInt[i].empty() && tsHadNearAtInt[i].empty() && !neutralMask[i]) { // nothing linked to this ts + neutralCandidate.addTrackster(edm::Ptr(tsH, i)); + neutralMask[i] = 1; + neutralCandidates.push_back(neutralCandidate); + continue; + } + if (!neutralMask[i]) { + neutralCandidate.addTrackster(edm::Ptr(tsH, i)); + neutralMask[i] = 1; + } + for (const unsigned ts2_idx : tsNearAtInt[i]) { + if (chargedMask[ts2_idx]) + continue; + if (!neutralMask[ts2_idx]) { + neutralCandidate.addTrackster(edm::Ptr(tsH, ts2_idx)); + neutralMask[ts2_idx] = 1; + } + for (const unsigned ts1_idx : tsHadNearAtInt[ts2_idx]) { + if (chargedMask[ts1_idx]) + continue; + if (!neutralMask[ts1_idx]) { + neutralCandidate.addTrackster(edm::Ptr(tsH, ts1_idx)); + neutralMask[ts1_idx] = 1; + } + } + } + for (const unsigned ts1_idx : tsHadNearAtInt[i]) { + if (chargedMask[ts1_idx]) + continue; + if (!neutralMask[ts1_idx]) { + neutralCandidate.addTrackster(edm::Ptr(tsH, ts1_idx)); + neutralMask[ts1_idx] = 1; + } + } + // filter empty candidates + if (!neutralCandidate.tracksters().empty()) { + neutralCandidates.push_back(neutralCandidate); + } + } + + // set other attributes of created candidates + for (auto &cand : chargedCandidates) { + bool isHAD = false; + double rawE = 0.; + const auto track = cand.trackPtr(); + for (const auto &ts : cand.tracksters()) { + // isHAD if atleast one trackster is not EM + if (isHadron(*ts)) + isHAD = true; + rawE += ts->raw_energy(); + } + auto pdgID = isHAD ? 211 : 11; + + cand.setCharge(track->charge()); + cand.setPdgId(pdgID * track->charge()); + cand.setRawEnergy(rawE); + math::XYZTLorentzVector p4(rawE * track->momentum().unit().x(), + rawE * track->momentum().unit().y(), + rawE * track->momentum().unit().z(), + rawE); + cand.setP4(p4); + } + + for (auto &cand : neutralCandidates) { + bool isHAD = false; + double rawE = 0.; + const auto track = cand.trackPtr(); + double wtSum_baryc[3] = {0}; + for (const auto &ts : cand.tracksters()) { + if (isHadron(*ts)) + isHAD = true; + rawE += ts->raw_energy(); + wtSum_baryc[0] += (ts->raw_energy()) * (ts->barycenter().x()); + wtSum_baryc[1] += (ts->raw_energy()) * (ts->barycenter().y()); + wtSum_baryc[2] += (ts->raw_energy()) * (ts->barycenter().z()); + } + Vector combined_baryc(wtSum_baryc[0] / rawE, wtSum_baryc[1] / rawE, wtSum_baryc[2] / rawE); + auto pdgID = isHAD ? 130 : 22; + auto const &magnitude = isHAD ? std::sqrt(rawE * rawE - ticl::mpion2) : rawE; + + cand.setCharge(0); + cand.setPdgId(pdgID); + cand.setRawEnergy(rawE); + math::XYZTLorentzVector p4(magnitude * combined_baryc.unit().x(), + magnitude * combined_baryc.unit().y(), + magnitude * combined_baryc.unit().z(), + rawE); + + cand.setP4(p4); + } + + resultLinked.insert(std::end(resultLinked), std::begin(neutralCandidates), std::end(neutralCandidates)); + resultLinked.insert(std::end(resultLinked), std::begin(chargedCandidates), std::end(chargedCandidates)); + resultLinked.insert(std::end(resultLinked), std::begin(chargedHadronsFromTk), std::end(chargedHadronsFromTk)); + +} // linkTracksters + +void LinkingAlgoByDirectionGeometric::fillPSetDescription(edm::ParameterSetDescription &desc) { + desc.add("cutTk", + "1.48 < abs(eta) < 3.0 && pt > 1. && quality(\"highPurity\") && " + "hitPattern().numberOfLostHits(\"MISSING_OUTER_HITS\") < 5"); + desc.add("delta_tk_ts_layer1", 0.02); + desc.add("delta_tk_ts_interface", 0.03); + desc.add("delta_ts_em_had", 0.03); + desc.add("delta_ts_had_had", 0.03); + desc.add("track_time_quality_threshold", 0.5); + desc.add("pid_threshold", 0.5); + desc.add("energy_em_over_total_threshold", 0.9); + desc.add>("filter_hadronic_on_categories", {0, 1}); + LinkingAlgoBase::fillPSetDescription(desc); +} diff --git a/RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.h b/RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.h new file mode 100644 index 0000000000000..e1f8d3ac9aaa2 --- /dev/null +++ b/RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.h @@ -0,0 +1,117 @@ +#ifndef RecoHGCal_TICL_LinkingAlgoByDirectionGeometric_H__ +#define RecoHGCal_TICL_LinkingAlgoByDirectionGeometric_H__ + +#include +#include +#include "RecoHGCal/TICL/plugins/LinkingAlgoBase.h" +#include "RecoHGCal/TICL/interface/commons.h" + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" + +#include "DataFormats/Math/interface/Vector3D.h" +#include "DataFormats/GeometrySurface/interface/BoundDisk.h" +#include "DataFormats/HGCalReco/interface/TICLLayerTile.h" + +#include "TrackingTools/Records/interface/TrackingComponentsRecord.h" + +#include "CommonTools/Utils/interface/StringCutObjectSelector.h" + +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" + +namespace ticl { + class LinkingAlgoByDirectionGeometric final : public LinkingAlgoBase { + public: + LinkingAlgoByDirectionGeometric(const edm::ParameterSet &conf); + ~LinkingAlgoByDirectionGeometric() override; + + void initialize(const HGCalDDDConstants *hgcons, + const hgcal::RecHitTools rhtools, + const edm::ESHandle bfieldH, + const edm::ESHandle propH) override; + + void linkTracksters(const edm::Handle>, + const edm::ValueMap &, + const edm::ValueMap &, + const edm::ValueMap &, + const std::vector &, + const edm::Handle>, + std::vector &) override; + + static void fillPSetDescription(edm::ParameterSetDescription &desc); + + private: + typedef math::XYZVector Vector; + + void buildLayers(); + + math::XYZVector propagateTrackster(const Trackster &t, + const unsigned idx, + float zVal, + std::array &tracksterTiles); + + void findTrackstersInWindow(const std::vector> &seedingCollection, + const std::array &tracksterTiles, + double delta, + unsigned trackstersSize, + std::vector> &resultCollection, + bool useMask); + + bool timeAndEnergyCompatible(double &total_raw_energy, + const reco::Track &track, + const Trackster &trackster, + const float &tkTime, + const float &tkTimeErr, + const float &tkTimeQual); + + void recordTrackster(const unsigned ts, // trackster index + const std::vector &tracksters, + const edm::Handle> tsH, + std::vector &ts_mask, + double &energy_in_candidate, + TICLCandidate &candidate); + + void addTracksterIfCompatible(const unsigned ts, + const edm::Handle> tsH, + const unsigned tk, + const reco::TrackRef &tkRef, + const std::vector &tracksters, + const std::vector &tracks, + std::vector &ts_mask, + double &energy_in_candidate, + TICLCandidate &candidate, + const edm::ValueMap &tkTime, + const edm::ValueMap &tkTimeErr, + const edm::ValueMap &tkTimeQual); + + void dumpLinksFound(std::vector> &resultCollection, const char *label) const; + + const double tkEnergyCut_ = 2.0; + const double maxDeltaT_ = 3.; + const double del_tk_ts_layer1_; + const double del_tk_ts_int_; + const double del_ts_em_had_; + const double del_ts_had_had_; + + const double timing_quality_threshold_; + const double pid_threshold_; + const double energy_em_over_total_threshold_; + const std::vector filter_on_categories_; + + const StringCutObjectSelector cutTk_; + std::once_flag initializeGeometry_; + + const HGCalDDDConstants *hgcons_; + + std::unique_ptr firstDisk_[2]; + std::unique_ptr interfaceDisk_[2]; + + hgcal::RecHitTools rhtools_; + + edm::ESHandle bfield_; + edm::ESHandle propagator_; + }; +} // namespace ticl +#endif diff --git a/RecoHGCal/TICL/plugins/LinkingAlgoFactory.cc b/RecoHGCal/TICL/plugins/LinkingAlgoFactory.cc new file mode 100644 index 0000000000000..ddea59017633f --- /dev/null +++ b/RecoHGCal/TICL/plugins/LinkingAlgoFactory.cc @@ -0,0 +1,10 @@ +#include "FWCore/ParameterSet/interface/ValidatedPluginFactoryMacros.h" +#include "FWCore/ParameterSet/interface/ValidatedPluginMacros.h" +#include "LinkingAlgoFactory.h" +#include "LinkingAlgoByDirectionGeometric.h" + +EDM_REGISTER_VALIDATED_PLUGINFACTORY(LinkingAlgoFactory, "LinkingAlgoFactory"); + +DEFINE_EDM_VALIDATED_PLUGIN(LinkingAlgoFactory, + ticl::LinkingAlgoByDirectionGeometric, + "LinkingAlgoByDirectionGeometric"); diff --git a/RecoHGCal/TICL/plugins/LinkingAlgoFactory.h b/RecoHGCal/TICL/plugins/LinkingAlgoFactory.h new file mode 100644 index 0000000000000..5fa59121cd685 --- /dev/null +++ b/RecoHGCal/TICL/plugins/LinkingAlgoFactory.h @@ -0,0 +1,10 @@ +#ifndef RecoHGCAL_TICL_LinkingAlgoFactory_h +#define RecoHGCAL_TICL_LinkingAlgoFactory_h + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/PluginManager/interface/PluginFactory.h" +#include "LinkingAlgoBase.h" + +using LinkingAlgoFactory = edmplugin::PluginFactory; + +#endif diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc index 0fda8fcd9fc5b..c90c4d0ce8fee 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc @@ -66,7 +66,7 @@ void PatternRecognitionbyCA::makeTracksters( theGraph_->setVerbosity(PatternRecognitionAlgoBaseT::algo_verbosity_); theGraph_->clear(); - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::None) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::None) { LogDebug("HGCPatternRecoByCA") << "Making Tracksters with CA" << std::endl; } @@ -119,7 +119,7 @@ void PatternRecognitionbyCA::makeTracksters( effective_cluster_idx.insert(innerCluster); effective_cluster_idx.insert(outerCluster); - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { LogDebug("HGCPatternRecoByCA") << " New doublet " << doublet << " for trackster: " << result.size() << " InnerCl " << innerCluster << " " << input.layerClusters[innerCluster].x() << " " << input.layerClusters[innerCluster].y() << " " @@ -214,7 +214,7 @@ void PatternRecognitionbyCA::makeTracksters( const auto &t = tmpTracksters[selectedTrackstersIds[i]]; for (auto const lcId : t.vertices()) { layer_cluster_usage[lcId]++; - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Basic) + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Basic) LogDebug("HGCPatternRecoByCA") << "LayerID: " << lcId << " count: " << (int)layer_cluster_usage[lcId] << std::endl; } @@ -228,7 +228,7 @@ void PatternRecognitionbyCA::makeTracksters( assert(trackster.vertices().size() <= trackster.vertex_multiplicity().size()); for (size_t i = 0; i < trackster.vertices().size(); ++i) { trackster.vertex_multiplicity()[i] = layer_cluster_usage[trackster.vertices(i)]; - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Basic) + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Basic) LogDebug("HGCPatternRecoByCA") << "LayerID: " << trackster.vertices(i) << " count: " << (int)trackster.vertex_multiplicity(i) << std::endl; } @@ -254,7 +254,7 @@ void PatternRecognitionbyCA::makeTracksters( emptyTrackstersFromSeedsTRK(result, seedToTracksterAssociation, input.regions[0].collectionID); } - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { for (auto &trackster : result) { LogDebug("HGCPatternRecoByCA") << "Trackster characteristics: " << std::endl; LogDebug("HGCPatternRecoByCA") << "Size: " << trackster.vertices().size() << std::endl; @@ -285,7 +285,7 @@ void PatternRecognitionbyCA::mergeTrackstersTRK( for (unsigned int j = 1; j < numberOfTrackstersInSeed; ++j) { auto &thisTrackster = input[tracksters[j]]; updated_size += thisTrackster.vertices().size(); - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Basic) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Basic) { LogDebug("HGCPatternRecoByCA") << "Updated size: " << updated_size << std::endl; } outTrackster.vertices().reserve(updated_size); diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc index 16b0ab927f483..b15a60274e76b 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc @@ -56,7 +56,7 @@ void PatternRecognitionbyCLUE3D::dumpTiles(const TILES &tiles) const { for (int phi = 0; phi < nPhiBin; phi++) { int iphi = ((phi % nPhiBin + nPhiBin) % nPhiBin); if (!tiles[layer][offset + iphi].empty()) { - if (this->algo_verbosity_ > this->Advanced) { + if (this->algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Layer: " << layer << " ieta: " << ieta << " phi: " << phi << " " << tiles[layer][offset + iphi].size(); } @@ -70,7 +70,7 @@ template void PatternRecognitionbyCLUE3D::dumpTracksters(const std::vector> &layerIdx2layerandSoa, const int eventNumber, const std::vector &tracksters) const { - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "[evt, tracksterId, cells, prob_photon, prob_ele, prob_chad, prob_nhad, layer_i, x_i, y_i, eta_i, phi_i, " "energy_i, radius_i, rho_i, z_extension, delta_tr, delta_lyr, isSeed_i"; @@ -82,7 +82,7 @@ void PatternRecognitionbyCLUE3D::dumpTracksters(const std::vector::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D_NTP") << std::setw(4) << eventNumber << sep << std::setw(4) << num << sep << std::setw(4) << t.vertices().size() << sep << std::setw(8) << t.id_probability(ticl::Trackster::ParticleType::photon) << sep << std::setw(8) @@ -105,7 +105,7 @@ template void PatternRecognitionbyCLUE3D::dumpClusters(const TILES &tiles, const std::vector> &layerIdx2layerandSoa, const int eventNumber) const { - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "[evt, lyr, Seed, x, y, z, r/|z|, eta, phi, " "etab, phib, cells, enrgy, e/rho, rho, z_ext, " " dlt_tr, dlt_lyr, " @@ -116,7 +116,7 @@ void PatternRecognitionbyCLUE3D::dumpClusters(const TILES &tiles, auto const &thisLayer = clusters_[layer]; int num = 0; for (auto v : thisLayer.x) { - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << std::setw(4) << eventNumber << ", " << std::setw(3) << layer << ", " << std::setw(4) << thisLayer.isSeed[num] << ", " << std::setprecision(3) << std::fixed << v << ", " << thisLayer.y[num] @@ -140,7 +140,7 @@ void PatternRecognitionbyCLUE3D::dumpClusters(const TILES &tiles, // Skip masked layer clusters if ((layerandSoa.first == -1) && (layerandSoa.second == -1)) continue; - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "lcIdx: " << lcIdx << " on Layer: " << layerandSoa.first << " SOA: " << layerandSoa.second; } @@ -157,7 +157,7 @@ void PatternRecognitionbyCLUE3D::makeTracksters( return; const int eventNumber = input.ev.eventAuxiliary().event(); - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "New Event"; } @@ -169,7 +169,7 @@ void PatternRecognitionbyCLUE3D::makeTracksters( // Also, layers inside the HGCAL geometry start from 1. for (unsigned int i = 0; i < rhtools_.lastLayer(); ++i) { layersPosZ_.push_back(rhtools_.getPositionLayer(i + 1).z()); - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Layer " << i << " located at Z: " << layersPosZ_.back(); } } @@ -182,7 +182,7 @@ void PatternRecognitionbyCLUE3D::makeTracksters( unsigned int layerIdx = 0; for (auto const &lc : input.layerClusters) { if (input.mask[layerIdx] == 0.) { - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Skipping masked cluster: " << layerIdx; } layerIdx2layerandSoa.emplace_back(-1, -1); @@ -217,7 +217,7 @@ void PatternRecognitionbyCLUE3D::makeTracksters( // below, too. float radius_x = sqrt((sum_sqr_x - (sum_x * sum_x) * invClsize) * invClsize); float radius_y = sqrt((sum_sqr_y - (sum_y * sum_y) * invClsize) * invClsize); - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "cluster rx: " << std::setw(5) << radius_x << ", ry: " << std::setw(5) << radius_y << ", r: " << std::setw(5) << (radius_x + radius_y) << ", cells: " << std::setw(4) @@ -230,7 +230,7 @@ void PatternRecognitionbyCLUE3D::makeTracksters( // Silicon case if (rhtools_.isSilicon(detId)) { radius_x = radius_y = rhtools_.getRadiusToSide(detId); - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Single cell cluster in silicon, rx: " << std::setw(5) << radius_x << ", ry: " << std::setw(5) << radius_y; } @@ -238,7 +238,7 @@ void PatternRecognitionbyCLUE3D::makeTracksters( auto const &point = rhtools_.getPosition(detId); auto const &eta_phi_window = rhtools_.getScintDEtaDPhi(detId); radius_x = radius_y = point.perp() * eta_phi_window.second; - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Single cell cluster in scintillator. rx: " << std::setw(5) << radius_x << ", ry: " << std::setw(5) << radius_y << ", eta-span: " << std::setw(5) << eta_phi_window.first << ", phi-span: " << std::setw(5) @@ -280,7 +280,7 @@ void PatternRecognitionbyCLUE3D::makeTracksters( } auto nTracksters = findAndAssignTracksters(input.tiles, layerIdx2layerandSoa); - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Reconstructed " << nTracksters << " tracksters" << std::endl; dumpClusters(input.tiles, layerIdx2layerandSoa, eventNumber); } @@ -290,15 +290,15 @@ void PatternRecognitionbyCLUE3D::makeTracksters( for (unsigned int layer = 0; layer < clusters_.size(); ++layer) { const auto &thisLayer = clusters_[layer]; - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Examining Layer: " << layer; } for (unsigned int lc = 0; lc < thisLayer.x.size(); ++lc) { - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Trackster " << thisLayer.clusterIndex[lc]; } if (thisLayer.clusterIndex[lc] >= 0) { - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << " adding lcIdx: " << thisLayer.layerClusterOriginalIdx[lc]; } result[thisLayer.clusterIndex[lc]].vertices().push_back(thisLayer.layerClusterOriginalIdx[lc]); @@ -328,7 +328,7 @@ void PatternRecognitionbyCLUE3D::makeTracksters( // run energy regression and ID energyRegressionAndID(input.layerClusters, input.tfSession, result); - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { for (auto const &t : result) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Barycenter: " << t.barycenter(); edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "LCs: " << t.vertices().size(); @@ -338,13 +338,13 @@ void PatternRecognitionbyCLUE3D::makeTracksters( } // Dump Tracksters information - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { dumpTracksters(layerIdx2layerandSoa, eventNumber, result); } // Reset internal clusters_ structure of array for next event reset(); - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << std::endl; } } @@ -520,13 +520,13 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( float deltaLayersZ = std::abs(layersPosZ_[maxLayer % lastLayerPerSide] - layersPosZ_[minLayer % lastLayerPerSide]); for (int currentLayer = minLayer; currentLayer <= maxLayer; currentLayer++) { - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "RefLayer: " << layerId << " SoaIDX: " << i; edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "NextLayer: " << currentLayer; } const auto &tileOnLayer = tiles[currentLayer]; bool onSameLayer = (currentLayer == layerId); - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "onSameLayer: " << onSameLayer; } const int etaWindow = 2; @@ -535,7 +535,7 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( int etaBinMax = std::min(tileOnLayer.etaBin(clustersOnLayer.eta[i]) + etaWindow, nEtaBin); int phiBinMin = tileOnLayer.phiBin(clustersOnLayer.phi[i]) - phiWindow; int phiBinMax = tileOnLayer.phiBin(clustersOnLayer.phi[i]) + phiWindow; - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "eta: " << clustersOnLayer.eta[i]; edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "phi: " << clustersOnLayer.phi[i]; edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "etaBinMin: " << etaBinMin << ", etaBinMax: " << etaBinMax; @@ -543,12 +543,12 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( } for (int ieta = etaBinMin; ieta <= etaBinMax; ++ieta) { auto offset = ieta * nPhiBin; - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "offset: " << offset; } for (int iphi_it = phiBinMin; iphi_it <= phiBinMax; ++iphi_it) { int iphi = ((iphi_it % nPhiBin + nPhiBin) % nPhiBin); - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "iphi: " << iphi; edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Entries in tileBin: " << tileOnLayer[offset + iphi].size(); @@ -557,13 +557,13 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( auto const &layerandSoa = layerIdx2layerandSoa[otherClusterIdx]; // Skip masked layer clusters if ((layerandSoa.first == -1) && (layerandSoa.second == -1)) { - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Skipping masked layerIdx " << otherClusterIdx; } continue; } auto const &clustersLayer = clusters_[layerandSoa.first]; - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "OtherLayer: " << layerandSoa.first << " SoaIDX: " << layerandSoa.second; edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "OtherEta: " << clustersLayer.eta[layerandSoa.second]; @@ -600,7 +600,7 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( clustersLayer.eta[layerandSoa.second], clustersLayer.phi[layerandSoa.second]) < densityEtaPhiDistanceSqr_); } - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Distance[eta,phi]: " << reco::deltaR2(clustersOnLayer.eta[i], clustersOnLayer.phi[i], @@ -624,7 +624,7 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( clustersLayer.energy[layerandSoa.second]; clustersOnLayer.rho[i] += energyToAdd; clustersOnLayer.z_extension[i] = deltaLayersZ; - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Adding " << energyToAdd << " partial " << clustersOnLayer.rho[i]; } @@ -634,7 +634,7 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( } // end of loop over eta-bin region } // end of loop on the sibling layers if (rescaleDensityByZ_) { - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Rescaling original density: " << clustersOnLayer.rho[i] << " by Z: " << deltaLayersZ << " to final density/cm: " << clustersOnLayer.rho[i] / deltaLayersZ; @@ -658,7 +658,7 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( }; for (unsigned int i = 0; i < numberOfClusters; i++) { - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Starting searching nearestHigher on " << layerId << " with rho: " << clustersOnLayer.rho[i] << " at eta, phi: " << tiles[layerId].etaBin(clustersOnLayer.eta[i]) << ", " @@ -693,7 +693,7 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( auto offset = ieta * nPhiBin; for (int iphi_it = phiBinMin; iphi_it <= phiBinMax; ++iphi_it) { int iphi = ((iphi_it % nPhiBin + nPhiBin) % nPhiBin); - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Searching nearestHigher on " << currentLayer << " eta, phi: " << ieta << ", " << iphi_it << " " << iphi << " " << offset << " " << (offset + iphi); @@ -725,7 +725,7 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( (clustersOnOtherLayer.rho[layerandSoa.second] == clustersOnLayer.rho[i] && clustersOnOtherLayer.layerClusterOriginalIdx[layerandSoa.second] > clustersOnLayer.layerClusterOriginalIdx[i]); - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Searching nearestHigher on " << currentLayer << " with rho: " << clustersOnOtherLayer.rho[layerandSoa.second] @@ -745,7 +745,7 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( } // End of loop on layers bool foundNearestInFiducialVolume = (i_delta != maxDelta); - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "i_delta: " << i_delta << " passed: " << foundNearestInFiducialVolume << " " << i_nearestHigher.first << " " << i_nearestHigher.second << " distances: " << nearest_distances.first << ", " @@ -790,7 +790,7 @@ int PatternRecognitionbyCLUE3D::findAndAssignTracksters( bool isOutlier = (clustersOnLayer.delta[i].first > outlierMultiplier_ * critical_transverse_distance) && (clustersOnLayer.rho[i] < criticalDensity_); if (isSeed) { - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Found seed on Layer " << layer << " SOAidx: " << i << " assigned ClusterIdx: " << nTracksters; } @@ -799,7 +799,7 @@ int PatternRecognitionbyCLUE3D::findAndAssignTracksters( localStack.emplace_back(layer, i); } else if (!isOutlier) { auto [lyrIdx, soaIdx] = clustersOnLayer.nearestHigher[i]; - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Found follower on Layer " << layer << " SOAidx: " << i << " attached to cluster on layer: " << lyrIdx << " SOAidx: " << soaIdx; @@ -807,7 +807,7 @@ int PatternRecognitionbyCLUE3D::findAndAssignTracksters( if (lyrIdx >= 0) clusters_[lyrIdx].followers[soaIdx].emplace_back(layer, i); } else { - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Found Outlier on Layer " << layer << " SOAidx: " << i << " with rho: " << clustersOnLayer.rho[i] << " and delta: " << clustersOnLayer.delta[i].first << ", " << clustersOnLayer.delta[i].second; diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyFastJet.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyFastJet.cc index 21eba683e5ab2..145728c881e16 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyFastJet.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyFastJet.cc @@ -42,13 +42,13 @@ PatternRecognitionbyFastJet::PatternRecognitionbyFastJet(const edm::Param template void PatternRecognitionbyFastJet::buildJetAndTracksters(std::vector &fjInputs, std::vector &result) { - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Basic) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Basic) { edm::LogVerbatim("PatternRecogntionbyFastJet") << "Creating FastJet with " << fjInputs.size() << " LayerClusters in input"; } fastjet::ClusterSequence sequence(fjInputs, JetDefinition(antikt_algorithm, antikt_radius_)); auto jets = fastjet::sorted_by_pt(sequence.inclusive_jets(0)); - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Basic) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Basic) { edm::LogVerbatim("PatternRecogntionbyFastJet") << "FastJet produced " << jets.size() << " jets/trackster"; } @@ -59,14 +59,14 @@ void PatternRecognitionbyFastJet::buildJetAndTracksters(std::vector::algo_verbosity_ > PatternRecognitionAlgoBaseT::Basic) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Basic) { edm::LogVerbatim("PatternRecogntionbyFastJet") << "Jet has " << pj.constituents().size() << " components that are stored in trackster " << trackster_idx; } } trackster_idx++; } else { - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecogntionbyFastJet") << "Jet with " << pj.constituents().size() << " constituents discarded since too small wrt " << minNumLayerCluster_; @@ -105,18 +105,18 @@ void PatternRecognitionbyFastJet::makeTracksters( const auto &tileOnLayer = input.tiles[currentLayer]; for (int ieta = 0; ieta <= nEtaBin; ++ieta) { auto offset = ieta * nPhiBin; - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecogntionbyFastJet") << "offset: " << offset; } for (int iphi = 0; iphi <= nPhiBin; ++iphi) { - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecogntionbyFastJet") << "iphi: " << iphi; edm::LogVerbatim("PatternRecogntionbyFastJet") << "Entries in tileBin: " << tileOnLayer[offset + iphi].size(); } for (auto clusterIdx : tileOnLayer[offset + iphi]) { // Skip masked layer clusters if (input.mask[clusterIdx] == 0.) { - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Advanced) { edm::LogVerbatim("PatternRecogntionbyFastJet") << "Skipping masked layerIdx " << clusterIdx; } continue; @@ -144,7 +144,7 @@ void PatternRecognitionbyFastJet::makeTracksters( // run energy regression and ID energyRegressionAndID(input.layerClusters, input.tfSession, result); - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Basic) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > VerbosityLevel::Basic) { for (auto const &t : result) { edm::LogVerbatim("PatternRecogntionbyFastJet") << "Barycenter: " << t.barycenter(); edm::LogVerbatim("PatternRecogntionbyFastJet") << "LCs: " << t.vertices().size(); diff --git a/RecoHGCal/TICL/plugins/SeedingRegionAlgoBase.h b/RecoHGCal/TICL/plugins/SeedingRegionAlgoBase.h index 5721e2df4227b..c71a1aad7e748 100644 --- a/RecoHGCal/TICL/plugins/SeedingRegionAlgoBase.h +++ b/RecoHGCal/TICL/plugins/SeedingRegionAlgoBase.h @@ -33,8 +33,6 @@ namespace ticl { static void fillPSetDescription(edm::ParameterSetDescription& desc) { desc.add("algo_verbosity", 0); } - enum VerbosityLevel { None = 0, Basic, Advanced, Expert, Guru }; - protected: int algo_verbosity_; int algoId_; diff --git a/RecoHGCal/TICL/plugins/SimTrackstersProducer.cc b/RecoHGCal/TICL/plugins/SimTrackstersProducer.cc index 2abe1a0f7be40..b0510642313d8 100644 --- a/RecoHGCal/TICL/plugins/SimTrackstersProducer.cc +++ b/RecoHGCal/TICL/plugins/SimTrackstersProducer.cc @@ -41,6 +41,17 @@ class SimTrackstersProducer : public edm::stream::EDProducer<> { static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); void produce(edm::Event&, const edm::EventSetup&) override; + void addTrackster(const int& index, + const std::vector, std::pair>>& lcVec, + const std::vector& inputClusterMask, + const float& fractionCut_, + const float& energy, + const int& pdgId, + const int& charge, + const edm::ProductID& seed, + const Trackster::IterationIndex iter, + std::vector& output_mask, + std::vector& result); private: std::string detector_; @@ -99,6 +110,43 @@ void SimTrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& des descriptions.addWithDefaultLabel(desc); } +void SimTrackstersProducer::addTrackster( + const int& index, + const std::vector, std::pair>>& lcVec, + const std::vector& inputClusterMask, + const float& fractionCut_, + const float& energy, + const int& pdgId, + const int& charge, + const edm::ProductID& seed, + const Trackster::IterationIndex iter, + std::vector& output_mask, + std::vector& result) { + if (lcVec.empty()) + return; + + Trackster tmpTrackster; + tmpTrackster.zeroProbabilities(); + tmpTrackster.vertices().reserve(lcVec.size()); + tmpTrackster.vertex_multiplicity().reserve(lcVec.size()); + for (auto const& [lc, energyScorePair] : lcVec) { + if (inputClusterMask[lc.index()] > 0) { + double fraction = energyScorePair.first / lc->energy(); + if (fraction < fractionCut_) + continue; + tmpTrackster.vertices().push_back(lc.index()); + output_mask[lc.index()] -= fraction; + tmpTrackster.vertex_multiplicity().push_back(1. / fraction); + } + } + + tmpTrackster.setIdProbability(tracksterParticleTypeFromPdgId(pdgId, charge), 1.f); + tmpTrackster.setRegressedEnergy(energy); + tmpTrackster.setIteration(iter); + tmpTrackster.setSeed(seed, index); + result.emplace_back(tmpTrackster); +} + void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) { auto result = std::make_unique(); auto output_mask = std::make_unique>(); diff --git a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc index b76a1584f9686..38663b72f4d8a 100644 --- a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc @@ -2,23 +2,49 @@ #include "FWCore/Framework/interface/stream/EDProducer.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/PluginDescription.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" #include "FWCore/Utilities/interface/ESGetToken.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/Framework/interface/ConsumesCollector.h" #include "DataFormats/CaloRecHit/interface/CaloCluster.h" #include "DataFormats/HGCalReco/interface/Common.h" #include "DataFormats/HGCalReco/interface/TICLLayerTile.h" #include "DataFormats/HGCalReco/interface/Trackster.h" -#include "DataFormats/HGCalReco/interface/TICLSeedingRegion.h" #include "DataFormats/TrackReco/interface/Track.h" -#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h" +#include "DataFormats/MuonReco/interface/Muon.h" +#include "DataFormats/GeometrySurface/interface/BoundDisk.h" +#include "DataFormats/HGCalReco/interface/TICLCandidate.h" +#include "DataFormats/TrackReco/interface/TrackFwd.h" +#include "DataFormats/Math/interface/Vector3D.h" + #include "RecoHGCal/TICL/interface/GlobalCache.h" #include "PhysicsTools/TensorFlow/interface/TfGraphRecord.h" #include "PhysicsTools/TensorFlow/interface/TensorFlow.h" #include "PhysicsTools/TensorFlow/interface/TfGraphDefWrapper.h" -#include "DataFormats/HGCalReco/interface/TICLCandidate.h" +#include "RecoHGCal/TICL/plugins/LinkingAlgoBase.h" +#include "RecoHGCal/TICL/plugins/LinkingAlgoFactory.h" +#include "RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.h" + +#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h" +#include "PhysicsTools/TensorFlow/interface/TensorFlow.h" + +#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h" +#include "TrackingTools/GeomPropagators/interface/Propagator.h" +#include "TrackingTools/Records/interface/TrackingComponentsRecord.h" + +#include "MagneticField/Engine/interface/MagneticField.h" +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" + +#include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h" +#include "Geometry/Records/interface/IdealGeometryRecord.h" +#include "Geometry/CommonDetUnit/interface/GeomDet.h" #include "TrackstersPCA.h" @@ -31,8 +57,18 @@ class TrackstersMergeProducer : public edm::stream::EDProducer<> { void produce(edm::Event &, const edm::EventSetup &) override; static void fillDescriptions(edm::ConfigurationDescriptions &descriptions); + // static methods for handling the global cache + static std::unique_ptr initializeGlobalCache(const edm::ParameterSet &); + static void globalEndJob(TrackstersCache *); + + void beginJob(); + void endJob(); + + void beginRun(edm::Run const &iEvent, edm::EventSetup const &es) override; + private: typedef ticl::Trackster::IterationIndex TracksterIterIndex; + typedef math::XYZVector Vector; void fillTile(TICLTracksterTiles &, const std::vector &, TracksterIterIndex); @@ -43,19 +79,26 @@ class TrackstersMergeProducer : public edm::stream::EDProducer<> { void assignTimeToCandidates(std::vector &resultCandidates) const; void dumpTrackster(const Trackster &) const; - const edm::EDGetTokenT> tracksterstrkem_token_; - const edm::EDGetTokenT> trackstersem_token_; - const edm::EDGetTokenT> tracksterstrk_token_; - const edm::EDGetTokenT> trackstershad_token_; - const edm::EDGetTokenT> seedingTrk_token_; + std::unique_ptr linkingAlgo_; + + const edm::EDGetTokenT> tracksters_clue3d_token_; const edm::EDGetTokenT> clusters_token_; const edm::EDGetTokenT>> clustersTime_token_; const edm::EDGetTokenT> tracks_token_; - const edm::ESGetToken geometry_token_; + const edm::EDGetTokenT> tracks_time_token_; + const edm::EDGetTokenT> tracks_time_quality_token_; + const edm::EDGetTokenT> tracks_time_err_token_; + const edm::EDGetTokenT> muons_token_; const std::string tfDnnLabel_; const edm::ESGetToken tfDnnToken_; const tensorflow::Session *tfSession_; + const edm::ESGetToken geometry_token_; + const std::string detector_; + const std::string propName_; + + const edm::ESGetToken bfield_token_; + const edm::ESGetToken propagator_token_; const bool optimiseAcrossTracksters_; const int eta_bin_window_; const int phi_bin_window_; @@ -73,34 +116,45 @@ class TrackstersMergeProducer : public edm::stream::EDProducer<> { const double resol_calo_scale_had_; const double resol_calo_offset_em_; const double resol_calo_scale_em_; - const bool debug_; const std::string eidInputName_; const std::string eidOutputNameEnergy_; const std::string eidOutputNameId_; const float eidMinClusterEnergy_; const int eidNLayers_; const int eidNClusters_; + std::once_flag initializeGeometry_; + + const HGCalDDDConstants *hgcons_; + + std::unique_ptr firstDisk_[2]; tensorflow::Session *eidSession_; hgcal::RecHitTools rhtools_; static constexpr int eidNFeatures_ = 3; + + edm::ESGetToken hdc_token_; }; TrackstersMergeProducer::TrackstersMergeProducer(const edm::ParameterSet &ps) - : tracksterstrkem_token_(consumes>(ps.getParameter("tracksterstrkem"))), - trackstersem_token_(consumes>(ps.getParameter("trackstersem"))), - tracksterstrk_token_(consumes>(ps.getParameter("tracksterstrk"))), - trackstershad_token_(consumes>(ps.getParameter("trackstershad"))), - seedingTrk_token_(consumes>(ps.getParameter("seedingTrk"))), + : tracksters_clue3d_token_(consumes>(ps.getParameter("trackstersclue3d"))), clusters_token_(consumes>(ps.getParameter("layer_clusters"))), clustersTime_token_( consumes>>(ps.getParameter("layer_clustersTime"))), tracks_token_(consumes>(ps.getParameter("tracks"))), - geometry_token_(esConsumes()), + tracks_time_token_(consumes>(ps.getParameter("tracksTime"))), + tracks_time_quality_token_(consumes>(ps.getParameter("tracksTimeQual"))), + tracks_time_err_token_(consumes>(ps.getParameter("tracksTimeErr"))), + muons_token_(consumes>(ps.getParameter("muons"))), tfDnnLabel_(ps.getParameter("tfDnnLabel")), tfDnnToken_(esConsumes(edm::ESInputTag("", tfDnnLabel_))), tfSession_(nullptr), + geometry_token_(esConsumes()), + detector_(ps.getParameter("detector")), + propName_(ps.getParameter("propagator")), + bfield_token_(esConsumes()), + propagator_token_( + esConsumes(edm::ESInputTag("", propName_))), optimiseAcrossTracksters_(ps.getParameter("optimiseAcrossTracksters")), eta_bin_window_(ps.getParameter("eta_bin_window")), phi_bin_window_(ps.getParameter("phi_bin_window")), @@ -118,7 +172,6 @@ TrackstersMergeProducer::TrackstersMergeProducer(const edm::ParameterSet &ps) resol_calo_scale_had_(ps.getParameter("resol_calo_scale_had")), resol_calo_offset_em_(ps.getParameter("resol_calo_offset_em")), resol_calo_scale_em_(ps.getParameter("resol_calo_scale_em")), - debug_(ps.getParameter("debug")), eidInputName_(ps.getParameter("eid_input_name")), eidOutputNameEnergy_(ps.getParameter("eid_output_name_energy")), eidOutputNameId_(ps.getParameter("eid_output_name_id")), @@ -128,8 +181,33 @@ TrackstersMergeProducer::TrackstersMergeProducer(const edm::ParameterSet &ps) eidSession_(nullptr) { produces>(); produces>(); + + std::string detectorName_ = (detector_ == "HFNose") ? "HGCalHFNoseSensitive" : "HGCalEESensitive"; + hdc_token_ = + esConsumes(edm::ESInputTag("", detectorName_)); + + auto linkingPSet = ps.getParameter("linkingPSet"); + auto algoType = linkingPSet.getParameter("type"); + linkingAlgo_ = LinkingAlgoFactory::get()->create(algoType, linkingPSet); } +void TrackstersMergeProducer::beginJob() {} + +void TrackstersMergeProducer::endJob(){}; + +void TrackstersMergeProducer::beginRun(edm::Run const &iEvent, edm::EventSetup const &es) { + edm::ESHandle hdc = es.getHandle(hdc_token_); + hgcons_ = hdc.product(); + + edm::ESHandle geom = es.getHandle(geometry_token_); + rhtools_.setGeometry(*geom); + + edm::ESHandle bfield = es.getHandle(bfield_token_); + edm::ESHandle propagator = es.getHandle(propagator_token_); + + linkingAlgo_->initialize(hgcons_, rhtools_, bfield, propagator); +}; + void TrackstersMergeProducer::fillTile(TICLTracksterTiles &tracksterTile, const std::vector &tracksters, TracksterIterIndex tracksterIteration) { @@ -168,377 +246,120 @@ void TrackstersMergeProducer::dumpTrackster(const Trackster &t) const { } void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es) { - edm::ESHandle geom = es.getHandle(geometry_token_); - rhtools_.setGeometry(*geom); auto resultTrackstersMerged = std::make_unique>(); auto resultCandidates = std::make_unique>(); tfSession_ = es.getData(tfDnnToken_).getSession(); - TICLTracksterTiles tracksterTile; - std::vector usedTrackstersMerged; - std::vector indexInMergedCollTRKEM; - std::vector indexInMergedCollEM; - std::vector indexInMergedCollTRK; - std::vector indexInMergedCollHAD; - std::vector usedSeeds; + edm::Handle> trackstersclue3d_h; + evt.getByToken(tracksters_clue3d_token_, trackstersclue3d_h); - // associating seed to the index of the trackster in the merged collection and the iteration that found it - std::map>> seedToTracksterAssociator; edm::Handle> track_h; evt.getByToken(tracks_token_, track_h); const auto &tracks = *track_h; - edm::Handle> cluster_h; - evt.getByToken(clusters_token_, cluster_h); - const auto &layerClusters = *cluster_h; - - edm::Handle>> clustersTime_h; - evt.getByToken(clustersTime_token_, clustersTime_h); - const auto &layerClustersTimes = *clustersTime_h; - - edm::Handle> trackstersem_h; - evt.getByToken(trackstersem_token_, trackstersem_h); - const auto &trackstersEM = *trackstersem_h; - - edm::Handle> tracksterstrkem_h; - evt.getByToken(tracksterstrkem_token_, tracksterstrkem_h); - const auto &trackstersTRKEM = *tracksterstrkem_h; - - edm::Handle> tracksterstrk_h; - evt.getByToken(tracksterstrk_token_, tracksterstrk_h); - const auto &trackstersTRK = *tracksterstrk_h; - - edm::Handle> trackstershad_h; - evt.getByToken(trackstershad_token_, trackstershad_h); - const auto &trackstersHAD = *trackstershad_h; - - edm::Handle> seedingTrk_h; - evt.getByToken(seedingTrk_token_, seedingTrk_h); - const auto &seedingTrk = *seedingTrk_h; - usedSeeds.resize(tracks.size(), false); - - fillTile(tracksterTile, trackstersTRKEM, TracksterIterIndex::TRKEM); - fillTile(tracksterTile, trackstersEM, TracksterIterIndex::EM); - fillTile(tracksterTile, trackstersTRK, TracksterIterIndex::TRKHAD); - fillTile(tracksterTile, trackstersHAD, TracksterIterIndex::HAD); - - auto totalNumberOfTracksters = - trackstersTRKEM.size() + trackstersTRK.size() + trackstersEM.size() + trackstersHAD.size(); - resultTrackstersMerged->reserve(totalNumberOfTracksters); - usedTrackstersMerged.resize(totalNumberOfTracksters, false); - indexInMergedCollTRKEM.reserve(trackstersTRKEM.size()); - indexInMergedCollEM.reserve(trackstersEM.size()); - indexInMergedCollTRK.reserve(trackstersTRK.size()); - indexInMergedCollHAD.reserve(trackstersHAD.size()); - - if (debug_) { - printTrackstersDebug(trackstersTRKEM, "tracksterTRKEM"); - printTrackstersDebug(trackstersEM, "tracksterEM"); - printTrackstersDebug(trackstersTRK, "tracksterTRK"); - printTrackstersDebug(trackstersHAD, "tracksterHAD"); - } - - for (auto const &t : trackstersTRKEM) { - indexInMergedCollTRKEM.push_back(resultTrackstersMerged->size()); - seedToTracksterAssociator[t.seedIndex()].emplace_back(resultTrackstersMerged->size(), TracksterIterIndex::TRKEM); - resultTrackstersMerged->push_back(t); - } + const auto &layerClusters = evt.get(clusters_token_); + const auto &layerClustersTimes = evt.get(clustersTime_token_); + const auto &muons = evt.get(muons_token_); + const auto &trackTime = evt.get(tracks_time_token_); + const auto &trackTimeErr = evt.get(tracks_time_err_token_); + const auto &trackTimeQual = evt.get(tracks_time_quality_token_); + + // Linking + linkingAlgo_->linkTracksters( + track_h, trackTime, trackTimeErr, trackTimeQual, muons, trackstersclue3d_h, *resultCandidates); + + // Print debug info + LogDebug("TrackstersMergeProducer") << "Results from the linking step : " << std::endl + << "No. of Tracks : " << tracks.size() + << " No. of Tracksters : " << (*trackstersclue3d_h).size() << std::endl + << "(neutral candidates have track id -1)" << std::endl; + + std::vector &candidates = *resultCandidates; + for (const auto &cand : candidates) { + auto track_ptr = cand.trackPtr(); + auto trackster_ptrs = cand.tracksters(); + +#ifdef EDM_ML_DEBUG + auto track_idx = track_ptr.get() - (edm::Ptr(track_h, 0)).get(); + track_idx = (track_ptr.isNull()) ? -1 : track_idx; + LogDebug("TrackstersMergeProducer") << "PDG ID " << cand.pdgId() << " charge " << cand.charge() << " p " << cand.p() + << std::endl; + LogDebug("TrackstersMergeProducer") << "track id (p) : " << track_idx << " (" + << (track_ptr.isNull() ? -1 : track_ptr->p()) << ") " + << " trackster ids (E) : "; +#endif + + // Merge included tracksters + ticl::Trackster outTrackster; + auto updated_size = 0; + for (const auto &ts_ptr : trackster_ptrs) { +#ifdef EDM_ML_DEBUG + auto ts_idx = ts_ptr.get() - (edm::Ptr(trackstersclue3d_h, 0)).get(); + LogDebug("TrackstersMergeProducer") << ts_idx << " (" << ts_ptr->raw_energy() << ") "; +#endif + + auto &thisTrackster = *ts_ptr; + updated_size += thisTrackster.vertices().size(); + outTrackster.vertices().reserve(updated_size); + outTrackster.vertex_multiplicity().reserve(updated_size); + std::copy(std::begin(thisTrackster.vertices()), + std::end(thisTrackster.vertices()), + std::back_inserter(outTrackster.vertices())); + std::copy(std::begin(thisTrackster.vertex_multiplicity()), + std::end(thisTrackster.vertex_multiplicity()), + std::back_inserter(outTrackster.vertex_multiplicity())); + } - for (auto const &t : trackstersEM) { - indexInMergedCollEM.push_back(resultTrackstersMerged->size()); - resultTrackstersMerged->push_back(t); - } + LogDebug("TrackstersMergeProducer") << std::endl; - for (auto const &t : trackstersTRK) { - indexInMergedCollTRK.push_back(resultTrackstersMerged->size()); - seedToTracksterAssociator[t.seedIndex()].emplace_back(resultTrackstersMerged->size(), TracksterIterIndex::TRKHAD); - resultTrackstersMerged->push_back(t); - } + // Find duplicate LCs + auto &orig_vtx = outTrackster.vertices(); + auto vtx_sorted{orig_vtx}; + std::sort(std::begin(vtx_sorted), std::end(vtx_sorted)); + for (unsigned int iLC = 1; iLC < vtx_sorted.size(); ++iLC) { + if (vtx_sorted[iLC] == vtx_sorted[iLC - 1]) { + // Clean up duplicate LCs + const auto lcIdx = vtx_sorted[iLC]; + const auto firstEl = std::find(orig_vtx.begin(), orig_vtx.end(), lcIdx); + const auto firstPos = std::distance(std::begin(orig_vtx), firstEl); + auto iDup = std::find(std::next(firstEl), orig_vtx.end(), lcIdx); + while (iDup != orig_vtx.end()) { + orig_vtx.erase(iDup); + outTrackster.vertex_multiplicity().erase(outTrackster.vertex_multiplicity().begin() + + std::distance(std::begin(orig_vtx), iDup)); + outTrackster.vertex_multiplicity()[firstPos] -= 1; + iDup = std::find(std::next(firstEl), orig_vtx.end(), lcIdx); + }; + } + } - for (auto const &t : trackstersHAD) { - indexInMergedCollHAD.push_back(resultTrackstersMerged->size()); - resultTrackstersMerged->push_back(t); + outTrackster.zeroProbabilities(); + if (!track_ptr.isNull()) { + outTrackster.setSeed(track_h.id(), track_ptr.get() - (edm::Ptr(track_h, 0)).get()); + if (std::abs(cand.pdgId()) == 11) + outTrackster.setIdProbability(ticl::Trackster::ParticleType::electron, 1.f); + else + outTrackster.setIdProbability(ticl::Trackster::ParticleType::charged_hadron, 1.f); + } else { + if (cand.pdgId() == 22) + outTrackster.setIdProbability(ticl::Trackster::ParticleType::photon, 1.f); + else + outTrackster.setIdProbability(ticl::Trackster::ParticleType::neutral_hadron, 1.f); + } + if (!outTrackster.vertices().empty()) + resultTrackstersMerged->push_back(outTrackster); } assignPCAtoTracksters(*resultTrackstersMerged, layerClusters, layerClustersTimes, rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z()); - energyRegressionAndID(layerClusters, tfSession_, *resultTrackstersMerged); - - printTrackstersDebug(*resultTrackstersMerged, "TrackstersMergeProducer"); - - auto trackstersMergedHandle = evt.put(std::move(resultTrackstersMerged)); - - // TICL Candidate creation - // We start from neutrals first - - // Photons - for (unsigned i = 0; i < trackstersEM.size(); ++i) { - auto mergedIdx = indexInMergedCollEM[i]; - usedTrackstersMerged[mergedIdx] = true; - const auto &t = trackstersEM[i]; //trackster - TICLCandidate tmpCandidate; - tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, mergedIdx)); - tmpCandidate.setCharge(0); - tmpCandidate.setPdgId(22); - tmpCandidate.setRawEnergy(t.raw_energy()); - math::XYZTLorentzVector p4(t.raw_energy() * t.barycenter().unit().x(), - t.raw_energy() * t.barycenter().unit().y(), - t.raw_energy() * t.barycenter().unit().z(), - t.raw_energy()); - tmpCandidate.setP4(p4); - resultCandidates->push_back(tmpCandidate); - } - - // Neutral Hadrons - constexpr double mpion = 0.13957; - constexpr float mpion2 = mpion * mpion; - for (unsigned i = 0; i < trackstersHAD.size(); ++i) { - auto mergedIdx = indexInMergedCollHAD[i]; - usedTrackstersMerged[mergedIdx] = true; - const auto &t = trackstersHAD[i]; //trackster - TICLCandidate tmpCandidate; - tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, mergedIdx)); - tmpCandidate.setCharge(0); - tmpCandidate.setPdgId(130); - tmpCandidate.setRawEnergy(t.raw_energy()); - float momentum = std::sqrt(t.raw_energy() * t.raw_energy() - mpion2); - math::XYZTLorentzVector p4(momentum * t.barycenter().unit().x(), - momentum * t.barycenter().unit().y(), - momentum * t.barycenter().unit().z(), - t.raw_energy()); - tmpCandidate.setP4(p4); - resultCandidates->push_back(tmpCandidate); - } - - // Charged Particles - for (unsigned i = 0; i < trackstersTRKEM.size(); ++i) { - auto mergedIdx = indexInMergedCollTRKEM[i]; - if (!usedTrackstersMerged[mergedIdx]) { - const auto &t = trackstersTRKEM[i]; //trackster - auto trackIdx = t.seedIndex(); - auto const &track = tracks[trackIdx]; - if (!usedSeeds[trackIdx] and t.raw_energy() > 0) { - usedSeeds[trackIdx] = true; - usedTrackstersMerged[mergedIdx] = true; - - std::vector trackstersTRKwithSameSeed; - std::vector trackstersTRKEMwithSameSeed; - - for (const auto &tracksterIterationPair : seedToTracksterAssociator[trackIdx]) { - if (tracksterIterationPair.first != mergedIdx and !usedTrackstersMerged[tracksterIterationPair.first] and - trackstersMergedHandle->at(tracksterIterationPair.first).raw_energy() > 0.) { - if (tracksterIterationPair.second == TracksterIterIndex::TRKEM) { - trackstersTRKEMwithSameSeed.push_back(tracksterIterationPair.first); - } else if (tracksterIterationPair.second == TracksterIterIndex::TRKHAD) { - trackstersTRKwithSameSeed.push_back(tracksterIterationPair.first); - } - } - } - - float tracksterTotalRawPt = t.raw_pt(); - std::vector haloTrackstersTRKIdx; - bool foundCompatibleTRK = false; - - for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { - usedTrackstersMerged[otherTracksterIdx] = true; - tracksterTotalRawPt += trackstersMergedHandle->at(otherTracksterIdx).raw_pt(); - - // Check the X,Y,Z barycenter and merge if they are very close (halo) - if ((t.barycenter() - trackstersMergedHandle->at(otherTracksterIdx).barycenter()).mag2() < - halo_max_distance2_) { - haloTrackstersTRKIdx.push_back(otherTracksterIdx); - - } else { - foundCompatibleTRK = true; - } - } - - //check if there is 1-to-1 relationship - if (trackstersTRKEMwithSameSeed.empty()) { - if (foundCompatibleTRK) { - TICLCandidate tmpCandidate; - tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, mergedIdx)); - double raw_energy = t.raw_energy(); - - tmpCandidate.setCharge(track.charge()); - tmpCandidate.setTrackPtr(edm::Ptr(track_h, trackIdx)); - tmpCandidate.setPdgId(211 * track.charge()); - for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { - tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); - raw_energy += trackstersMergedHandle->at(otherTracksterIdx).raw_energy(); - } - tmpCandidate.setRawEnergy(raw_energy); - math::XYZTLorentzVector p4(raw_energy * track.momentum().unit().x(), - raw_energy * track.momentum().unit().y(), - raw_energy * track.momentum().unit().z(), - raw_energy); - tmpCandidate.setP4(p4); - resultCandidates->push_back(tmpCandidate); - - } else { - TICLCandidate tmpCandidate; - tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, mergedIdx)); - double raw_energy = t.raw_energy(); - tmpCandidate.setCharge(track.charge()); - tmpCandidate.setTrackPtr(edm::Ptr(track_h, trackIdx)); - for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { - tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); - raw_energy += trackstersMergedHandle->at(otherTracksterIdx).raw_energy(); - } - tmpCandidate.setPdgId(11 * track.charge()); - - tmpCandidate.setRawEnergy(raw_energy); - math::XYZTLorentzVector p4(raw_energy * track.momentum().unit().x(), - raw_energy * track.momentum().unit().y(), - raw_energy * track.momentum().unit().z(), - raw_energy); - tmpCandidate.setP4(p4); - resultCandidates->push_back(tmpCandidate); - } - - } else { - // if 1-to-many find closest trackster in momentum - int closestTrackster = mergedIdx; - float minPtDiff = std::abs(t.raw_pt() - track.pt()); - for (auto otherTracksterIdx : trackstersTRKEMwithSameSeed) { - auto thisPt = tracksterTotalRawPt + trackstersMergedHandle->at(otherTracksterIdx).raw_pt() - t.raw_pt(); - closestTrackster = std::abs(thisPt - track.pt()) < minPtDiff ? otherTracksterIdx : closestTrackster; - } - usedTrackstersMerged[closestTrackster] = true; - - if (foundCompatibleTRK) { - TICLCandidate tmpCandidate; - tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, closestTrackster)); - double raw_energy = trackstersMergedHandle->at(closestTrackster).raw_energy(); - - tmpCandidate.setCharge(track.charge()); - tmpCandidate.setTrackPtr(edm::Ptr(track_h, trackIdx)); - tmpCandidate.setPdgId(211 * track.charge()); - for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { - tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); - raw_energy += trackstersMergedHandle->at(otherTracksterIdx).raw_energy(); - } - tmpCandidate.setRawEnergy(raw_energy); - float momentum = std::sqrt(raw_energy * raw_energy - mpion2); - math::XYZTLorentzVector p4(momentum * track.momentum().unit().x(), - momentum * track.momentum().unit().y(), - momentum * track.momentum().unit().z(), - raw_energy); - tmpCandidate.setP4(p4); - resultCandidates->push_back(tmpCandidate); - - } else { - TICLCandidate tmpCandidate; - tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, closestTrackster)); - double raw_energy = trackstersMergedHandle->at(closestTrackster).raw_energy(); - - tmpCandidate.setCharge(track.charge()); - tmpCandidate.setTrackPtr(edm::Ptr(track_h, trackIdx)); - for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { - tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); - raw_energy += trackstersMergedHandle->at(otherTracksterIdx).raw_energy(); - } - tmpCandidate.setPdgId(11 * track.charge()); - tmpCandidate.setRawEnergy(raw_energy); - math::XYZTLorentzVector p4(raw_energy * track.momentum().unit().x(), - raw_energy * track.momentum().unit().y(), - raw_energy * track.momentum().unit().z(), - raw_energy); - tmpCandidate.setP4(p4); - resultCandidates->push_back(tmpCandidate); - } - // Promote all other TRKEM tracksters as photons with their energy. - for (auto otherTracksterIdx : trackstersTRKEMwithSameSeed) { - auto tmpIndex = (otherTracksterIdx != closestTrackster) ? otherTracksterIdx : mergedIdx; - TICLCandidate photonCandidate; - const auto &otherTrackster = trackstersMergedHandle->at(tmpIndex); - auto gammaEnergy = otherTrackster.raw_energy(); - photonCandidate.setCharge(0); - photonCandidate.setPdgId(22); - photonCandidate.setRawEnergy(gammaEnergy); - math::XYZTLorentzVector gammaP4(gammaEnergy * otherTrackster.barycenter().unit().x(), - gammaEnergy * otherTrackster.barycenter().unit().y(), - gammaEnergy * otherTrackster.barycenter().unit().z(), - gammaEnergy); - photonCandidate.setP4(gammaP4); - photonCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, tmpIndex)); - resultCandidates->push_back(photonCandidate); - } - } - } - } - } //end of loop over trackstersTRKEM - - for (unsigned i = 0; i < trackstersTRK.size(); ++i) { - auto mergedIdx = indexInMergedCollTRK[i]; - const auto &t = trackstersTRK[i]; //trackster - - if (!usedTrackstersMerged[mergedIdx] and t.raw_energy() > 0) { - auto trackIdx = t.seedIndex(); - auto const &track = tracks[trackIdx]; - if (!usedSeeds[trackIdx]) { - usedSeeds[trackIdx] = true; - usedTrackstersMerged[mergedIdx] = true; - TICLCandidate tmpCandidate; - tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, mergedIdx)); - tmpCandidate.setCharge(track.charge()); - tmpCandidate.setTrackPtr(edm::Ptr(track_h, trackIdx)); - tmpCandidate.setPdgId(211 * track.charge()); - tmpCandidate.setRawEnergy(t.raw_energy()); - float momentum = std::sqrt(t.raw_energy() * t.raw_energy() - mpion2); - math::XYZTLorentzVector p4(momentum * track.momentum().unit().x(), - momentum * track.momentum().unit().y(), - momentum * track.momentum().unit().z(), - t.raw_energy()); - tmpCandidate.setP4(p4); - resultCandidates->push_back(tmpCandidate); - } - } - } - // For all seeds that have 0-energy tracksters whose track is not marked as used, create a charged hadron with the track information. - for (auto const &s : seedingTrk) { - if (usedSeeds[s.index] == false) { - auto const &track = tracks[s.index]; - // emit a charged hadron - TICLCandidate tmpCandidate; - tmpCandidate.setCharge(track.charge()); - tmpCandidate.setTrackPtr(edm::Ptr(track_h, s.index)); - tmpCandidate.setPdgId(211 * track.charge()); - float energy = std::sqrt(track.p() * track.p() + mpion2); - tmpCandidate.setRawEnergy(energy); - math::PtEtaPhiMLorentzVector p4Polar(track.pt(), track.eta(), track.phi(), mpion); - tmpCandidate.setP4(p4Polar); - resultCandidates->push_back(tmpCandidate); - usedSeeds[s.index] = true; - } - } - - // for all general tracks (high purity, pt > 1), check if they have been used: if not, promote them as charged hadrons - for (unsigned i = 0; i < tracks.size(); ++i) { - auto const &track = tracks[i]; - if (track.pt() > track_min_pt_ and track.quality(reco::TrackBase::highPurity) and - track.missingOuterHits() < track_max_missing_outerhits_ and std::abs(track.outerEta()) > track_min_eta_ and - std::abs(track.outerEta()) < track_max_eta_ and usedSeeds[i] == false) { - // emit a charged hadron - TICLCandidate tmpCandidate; - tmpCandidate.setCharge(track.charge()); - tmpCandidate.setTrackPtr(edm::Ptr(track_h, i)); - tmpCandidate.setPdgId(211 * track.charge()); - float energy = std::sqrt(track.p() * track.p() + mpion2); - tmpCandidate.setRawEnergy(energy); - math::PtEtaPhiMLorentzVector p4Polar(track.pt(), track.eta(), track.phi(), mpion); - tmpCandidate.setP4(p4Polar); - resultCandidates->push_back(tmpCandidate); - usedSeeds[i] = true; - } - } // Compute timing assignTimeToCandidates(*resultCandidates); + evt.put(std::move(resultTrackstersMerged)); evt.put(std::move(resultCandidates)); } @@ -707,9 +528,7 @@ void TrackstersMergeProducer::assignTimeToCandidates(std::vector } void TrackstersMergeProducer::printTrackstersDebug(const std::vector &tracksters, const char *label) const { - if (!debug_) - return; - +#ifdef EDM_ML_DEBUG int counter = 0; for (auto const &t : tracksters) { LogDebug("TrackstersMergeProducer") @@ -731,18 +550,26 @@ void TrackstersMergeProducer::printTrackstersDebug(const std::vector } LogDebug("TrackstersMergeProducer") << std::endl; } +#endif } void TrackstersMergeProducer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) { edm::ParameterSetDescription desc; - desc.add("tracksterstrkem", edm::InputTag("ticlTrackstersTrkEM")); - desc.add("trackstersem", edm::InputTag("ticlTrackstersEM")); - desc.add("tracksterstrk", edm::InputTag("ticlTrackstersTrk")); - desc.add("trackstershad", edm::InputTag("ticlTrackstersHAD")); - desc.add("seedingTrk", edm::InputTag("ticlSeedingTrk")); + + edm::ParameterSetDescription linkingDesc; + linkingDesc.addNode(edm::PluginDescription("type", "LinkingAlgoByDirectionGeometric", true)); + desc.add("linkingPSet", linkingDesc); + + desc.add("trackstersclue3d", edm::InputTag("ticlTrackstersCLUE3DHigh")); desc.add("layer_clusters", edm::InputTag("hgcalLayerClusters")); desc.add("layer_clustersTime", edm::InputTag("hgcalLayerClusters", "timeLayerCluster")); desc.add("tracks", edm::InputTag("generalTracks")); + desc.add("tracksTime", edm::InputTag("tofPID:t0")); + desc.add("tracksTimeQual", edm::InputTag("mtdTrackQualityMVA:mtdQualMVA")); + desc.add("tracksTimeErr", edm::InputTag("tofPID:sigmat0")); + desc.add("muons", edm::InputTag("muons1stStep")); + desc.add("detector", "HGCAL"); + desc.add("propagator", "PropagatorWithMaterial"); desc.add("optimiseAcrossTracksters", true); desc.add("eta_bin_window", 1); desc.add("phi_bin_window", 1); @@ -760,7 +587,6 @@ void TrackstersMergeProducer::fillDescriptions(edm::ConfigurationDescriptions &d desc.add("resol_calo_scale_had", 0.15); desc.add("resol_calo_offset_em", 1.5); desc.add("resol_calo_scale_em", 0.15); - desc.add("debug", true); desc.add("tfDnnLabel", "tracksterSelectionTf"); desc.add("eid_input_name", "input"); desc.add("eid_output_name_energy", "output/regressed_energy"); @@ -771,5 +597,4 @@ void TrackstersMergeProducer::fillDescriptions(edm::ConfigurationDescriptions &d descriptions.add("trackstersMergeProducer", desc); } -#include "FWCore/Framework/interface/MakerMacros.h" DEFINE_FWK_MODULE(TrackstersMergeProducer); diff --git a/RecoHGCal/TICL/plugins/TrackstersMergeProducerV3.cc b/RecoHGCal/TICL/plugins/TrackstersMergeProducerV3.cc new file mode 100644 index 0000000000000..c5c7162ee2e8c --- /dev/null +++ b/RecoHGCal/TICL/plugins/TrackstersMergeProducerV3.cc @@ -0,0 +1,752 @@ +#include // unique_ptr +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/Utilities/interface/ESGetToken.h" + +#include "DataFormats/CaloRecHit/interface/CaloCluster.h" +#include "DataFormats/HGCalReco/interface/Common.h" +#include "DataFormats/HGCalReco/interface/TICLLayerTile.h" +#include "DataFormats/HGCalReco/interface/Trackster.h" +#include "DataFormats/HGCalReco/interface/TICLSeedingRegion.h" +#include "DataFormats/TrackReco/interface/Track.h" +#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h" +#include "RecoHGCal/TICL/interface/GlobalCache.h" + +#include "PhysicsTools/TensorFlow/interface/TfGraphRecord.h" +#include "PhysicsTools/TensorFlow/interface/TensorFlow.h" +#include "PhysicsTools/TensorFlow/interface/TfGraphDefWrapper.h" + +#include "DataFormats/HGCalReco/interface/TICLCandidate.h" + +#include "TrackstersPCA.h" + +using namespace ticl; + +class TrackstersMergeProducerV3 : public edm::stream::EDProducer<> { +public: + explicit TrackstersMergeProducerV3(const edm::ParameterSet &ps); + ~TrackstersMergeProducerV3() override{}; + void produce(edm::Event &, const edm::EventSetup &) override; + static void fillDescriptions(edm::ConfigurationDescriptions &descriptions); + +private: + typedef ticl::Trackster::IterationIndex TracksterIterIndex; + + void fillTile(TICLTracksterTiles &, const std::vector &, TracksterIterIndex); + + void energyRegressionAndID(const std::vector &layerClusters, + const tensorflow::Session *, + std::vector &result) const; + void printTrackstersDebug(const std::vector &, const char *label) const; + void assignTimeToCandidates(std::vector &resultCandidates) const; + void dumpTrackster(const Trackster &) const; + + const edm::EDGetTokenT> tracksterstrkem_token_; + const edm::EDGetTokenT> trackstersem_token_; + const edm::EDGetTokenT> tracksterstrk_token_; + const edm::EDGetTokenT> trackstershad_token_; + const edm::EDGetTokenT> seedingTrk_token_; + const edm::EDGetTokenT> clusters_token_; + const edm::EDGetTokenT>> clustersTime_token_; + const edm::EDGetTokenT> tracks_token_; + const edm::ESGetToken geometry_token_; + const std::string tfDnnLabel_; + const edm::ESGetToken tfDnnToken_; + const tensorflow::Session *tfSession_; + + const bool optimiseAcrossTracksters_; + const int eta_bin_window_; + const int phi_bin_window_; + const double pt_sigma_high_; + const double pt_sigma_low_; + const double halo_max_distance2_; + const double track_min_pt_; + const double track_min_eta_; + const double track_max_eta_; + const int track_max_missing_outerhits_; + const double cosangle_align_; + const double e_over_h_threshold_; + const double pt_neutral_threshold_; + const double resol_calo_offset_had_; + const double resol_calo_scale_had_; + const double resol_calo_offset_em_; + const double resol_calo_scale_em_; + const std::string eidInputName_; + const std::string eidOutputNameEnergy_; + const std::string eidOutputNameId_; + const float eidMinClusterEnergy_; + const int eidNLayers_; + const int eidNClusters_; + + tensorflow::Session *eidSession_; + hgcal::RecHitTools rhtools_; + + static constexpr int eidNFeatures_ = 3; +}; + +TrackstersMergeProducerV3::TrackstersMergeProducerV3(const edm::ParameterSet &ps) + : tracksterstrkem_token_(consumes>(ps.getParameter("tracksterstrkem"))), + trackstersem_token_(consumes>(ps.getParameter("trackstersem"))), + tracksterstrk_token_(consumes>(ps.getParameter("tracksterstrk"))), + trackstershad_token_(consumes>(ps.getParameter("trackstershad"))), + seedingTrk_token_(consumes>(ps.getParameter("seedingTrk"))), + clusters_token_(consumes>(ps.getParameter("layer_clusters"))), + clustersTime_token_( + consumes>>(ps.getParameter("layer_clustersTime"))), + tracks_token_(consumes>(ps.getParameter("tracks"))), + geometry_token_(esConsumes()), + tfDnnLabel_(ps.getParameter("tfDnnLabel")), + tfDnnToken_(esConsumes(edm::ESInputTag("", tfDnnLabel_))), + tfSession_(nullptr), + optimiseAcrossTracksters_(ps.getParameter("optimiseAcrossTracksters")), + eta_bin_window_(ps.getParameter("eta_bin_window")), + phi_bin_window_(ps.getParameter("phi_bin_window")), + pt_sigma_high_(ps.getParameter("pt_sigma_high")), + pt_sigma_low_(ps.getParameter("pt_sigma_low")), + halo_max_distance2_(ps.getParameter("halo_max_distance2")), + track_min_pt_(ps.getParameter("track_min_pt")), + track_min_eta_(ps.getParameter("track_min_eta")), + track_max_eta_(ps.getParameter("track_max_eta")), + track_max_missing_outerhits_(ps.getParameter("track_max_missing_outerhits")), + cosangle_align_(ps.getParameter("cosangle_align")), + e_over_h_threshold_(ps.getParameter("e_over_h_threshold")), + pt_neutral_threshold_(ps.getParameter("pt_neutral_threshold")), + resol_calo_offset_had_(ps.getParameter("resol_calo_offset_had")), + resol_calo_scale_had_(ps.getParameter("resol_calo_scale_had")), + resol_calo_offset_em_(ps.getParameter("resol_calo_offset_em")), + resol_calo_scale_em_(ps.getParameter("resol_calo_scale_em")), + eidInputName_(ps.getParameter("eid_input_name")), + eidOutputNameEnergy_(ps.getParameter("eid_output_name_energy")), + eidOutputNameId_(ps.getParameter("eid_output_name_id")), + eidMinClusterEnergy_(ps.getParameter("eid_min_cluster_energy")), + eidNLayers_(ps.getParameter("eid_n_layers")), + eidNClusters_(ps.getParameter("eid_n_clusters")), + eidSession_(nullptr) { + produces>(); + produces>(); +} + +void TrackstersMergeProducerV3::fillTile(TICLTracksterTiles &tracksterTile, + const std::vector &tracksters, + TracksterIterIndex tracksterIteration) { + int tracksterId = 0; + for (auto const &t : tracksters) { + tracksterTile.fill(tracksterIteration, t.barycenter().eta(), t.barycenter().phi(), tracksterId); + LogDebug("TrackstersMergeProducerV3") << "Adding tracksterId: " << tracksterId << " into bin [eta,phi]: [ " + << tracksterTile[tracksterIteration].etaBin(t.barycenter().eta()) << ", " + << tracksterTile[tracksterIteration].phiBin(t.barycenter().phi()) + << "] for iteration: " << tracksterIteration << std::endl; + + tracksterId++; + } +} + +void TrackstersMergeProducerV3::dumpTrackster(const Trackster &t) const { + auto e_over_h = (t.raw_em_pt() / ((t.raw_pt() - t.raw_em_pt()) != 0. ? (t.raw_pt() - t.raw_em_pt()) : 1.)); + LogDebug("TrackstersMergeProducerV3") + << "\nTrackster raw_pt: " << t.raw_pt() << " raw_em_pt: " << t.raw_em_pt() << " eoh: " << e_over_h + << " barycenter: " << t.barycenter() << " eta,phi (baricenter): " << t.barycenter().eta() << ", " + << t.barycenter().phi() << " eta,phi (eigen): " << t.eigenvectors(0).eta() << ", " << t.eigenvectors(0).phi() + << " pt(eigen): " << std::sqrt(t.eigenvectors(0).Unit().perp2()) * t.raw_energy() << " seedID: " << t.seedID() + << " seedIndex: " << t.seedIndex() << " size: " << t.vertices().size() << " average usage: " + << (std::accumulate(std::begin(t.vertex_multiplicity()), std::end(t.vertex_multiplicity()), 0.) / + (float)t.vertex_multiplicity().size()) + << " raw_energy: " << t.raw_energy() << " regressed energy: " << t.regressed_energy() + << " probs(ga/e/mu/np/cp/nh/am/unk): "; + for (auto const &p : t.id_probabilities()) { + LogDebug("TrackstersMergeProducerV3") << std::fixed << p << " "; + } + LogDebug("TrackstersMergeProducerV3") << " sigmas: "; + for (auto const &s : t.sigmas()) { + LogDebug("TrackstersMergeProducerV3") << s << " "; + } + LogDebug("TrackstersMergeProducerV3") << std::endl; +} + +void TrackstersMergeProducerV3::produce(edm::Event &evt, const edm::EventSetup &es) { + edm::ESHandle geom = es.getHandle(geometry_token_); + rhtools_.setGeometry(*geom); + auto resultTrackstersMerged = std::make_unique>(); + auto resultCandidates = std::make_unique>(); + + tfSession_ = es.getData(tfDnnToken_).getSession(); + + TICLTracksterTiles tracksterTile; + std::vector usedTrackstersMerged; + std::vector indexInMergedCollTRKEM; + std::vector indexInMergedCollEM; + std::vector indexInMergedCollTRK; + std::vector indexInMergedCollHAD; + std::vector usedSeeds; + + // associating seed to the index of the trackster in the merged collection and the iteration that found it + std::map>> seedToTracksterAssociator; + + edm::Handle> track_h; + evt.getByToken(tracks_token_, track_h); + const auto &tracks = *track_h; + + const auto &layerClusters = evt.get(clusters_token_); + const auto &layerClustersTimes = evt.get(clustersTime_token_); + const auto &trackstersEM = evt.get(trackstersem_token_); + const auto &trackstersTRKEM = evt.get(tracksterstrkem_token_); + const auto &trackstersTRK = evt.get(tracksterstrk_token_); + const auto &trackstersHAD = evt.get(trackstershad_token_); + const auto &seedingTrk = evt.get(seedingTrk_token_); + + usedSeeds.resize(tracks.size(), false); + + fillTile(tracksterTile, trackstersTRKEM, TracksterIterIndex::TRKEM); + fillTile(tracksterTile, trackstersEM, TracksterIterIndex::EM); + fillTile(tracksterTile, trackstersTRK, TracksterIterIndex::TRKHAD); + fillTile(tracksterTile, trackstersHAD, TracksterIterIndex::HAD); + + auto totalNumberOfTracksters = + trackstersTRKEM.size() + trackstersTRK.size() + trackstersEM.size() + trackstersHAD.size(); + resultTrackstersMerged->reserve(totalNumberOfTracksters); + usedTrackstersMerged.resize(totalNumberOfTracksters, false); + indexInMergedCollTRKEM.reserve(trackstersTRKEM.size()); + indexInMergedCollEM.reserve(trackstersEM.size()); + indexInMergedCollTRK.reserve(trackstersTRK.size()); + indexInMergedCollHAD.reserve(trackstersHAD.size()); + + printTrackstersDebug(trackstersTRKEM, "tracksterTRKEM"); + printTrackstersDebug(trackstersEM, "tracksterEM"); + printTrackstersDebug(trackstersTRK, "tracksterTRK"); + printTrackstersDebug(trackstersHAD, "tracksterHAD"); + + for (auto const &t : trackstersTRKEM) { + indexInMergedCollTRKEM.push_back(resultTrackstersMerged->size()); + seedToTracksterAssociator[t.seedIndex()].emplace_back(resultTrackstersMerged->size(), TracksterIterIndex::TRKEM); + resultTrackstersMerged->push_back(t); + } + + for (auto const &t : trackstersEM) { + indexInMergedCollEM.push_back(resultTrackstersMerged->size()); + resultTrackstersMerged->push_back(t); + } + + for (auto const &t : trackstersTRK) { + indexInMergedCollTRK.push_back(resultTrackstersMerged->size()); + seedToTracksterAssociator[t.seedIndex()].emplace_back(resultTrackstersMerged->size(), TracksterIterIndex::TRKHAD); + resultTrackstersMerged->push_back(t); + } + + for (auto const &t : trackstersHAD) { + indexInMergedCollHAD.push_back(resultTrackstersMerged->size()); + resultTrackstersMerged->push_back(t); + } + + assignPCAtoTracksters(*resultTrackstersMerged, + layerClusters, + layerClustersTimes, + rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z()); + energyRegressionAndID(layerClusters, tfSession_, *resultTrackstersMerged); + + printTrackstersDebug(*resultTrackstersMerged, "TrackstersMergeProducerV3"); + + auto trackstersMergedHandle = evt.put(std::move(resultTrackstersMerged)); + + // TICL Candidate creation + // We start from neutrals first + + // Photons + for (unsigned i = 0; i < trackstersEM.size(); ++i) { + auto mergedIdx = indexInMergedCollEM[i]; + usedTrackstersMerged[mergedIdx] = true; + const auto &t = trackstersEM[i]; //trackster + TICLCandidate tmpCandidate; + tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, mergedIdx)); + tmpCandidate.setCharge(0); + tmpCandidate.setPdgId(22); + tmpCandidate.setRawEnergy(t.raw_energy()); + math::XYZTLorentzVector p4(t.raw_energy() * t.barycenter().unit().x(), + t.raw_energy() * t.barycenter().unit().y(), + t.raw_energy() * t.barycenter().unit().z(), + t.raw_energy()); + tmpCandidate.setP4(p4); + resultCandidates->push_back(tmpCandidate); + } + + // Neutral Hadrons + constexpr double mpion = 0.13957; + constexpr float mpion2 = mpion * mpion; + for (unsigned i = 0; i < trackstersHAD.size(); ++i) { + auto mergedIdx = indexInMergedCollHAD[i]; + usedTrackstersMerged[mergedIdx] = true; + const auto &t = trackstersHAD[i]; //trackster + TICLCandidate tmpCandidate; + tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, mergedIdx)); + tmpCandidate.setCharge(0); + tmpCandidate.setPdgId(130); + tmpCandidate.setRawEnergy(t.raw_energy()); + float momentum = std::sqrt(t.raw_energy() * t.raw_energy() - mpion2); + math::XYZTLorentzVector p4(momentum * t.barycenter().unit().x(), + momentum * t.barycenter().unit().y(), + momentum * t.barycenter().unit().z(), + t.raw_energy()); + tmpCandidate.setP4(p4); + resultCandidates->push_back(tmpCandidate); + } + + // Charged Particles + for (unsigned i = 0; i < trackstersTRKEM.size(); ++i) { + auto mergedIdx = indexInMergedCollTRKEM[i]; + if (!usedTrackstersMerged[mergedIdx]) { + const auto &t = trackstersTRKEM[i]; //trackster + auto trackIdx = t.seedIndex(); + auto const &track = tracks[trackIdx]; + if (!usedSeeds[trackIdx] and t.raw_energy() > 0) { + usedSeeds[trackIdx] = true; + usedTrackstersMerged[mergedIdx] = true; + + std::vector trackstersTRKwithSameSeed; + std::vector trackstersTRKEMwithSameSeed; + + for (const auto &tracksterIterationPair : seedToTracksterAssociator[trackIdx]) { + if (tracksterIterationPair.first != mergedIdx and !usedTrackstersMerged[tracksterIterationPair.first] and + trackstersMergedHandle->at(tracksterIterationPair.first).raw_energy() > 0.) { + if (tracksterIterationPair.second == TracksterIterIndex::TRKEM) { + trackstersTRKEMwithSameSeed.push_back(tracksterIterationPair.first); + } else if (tracksterIterationPair.second == TracksterIterIndex::TRKHAD) { + trackstersTRKwithSameSeed.push_back(tracksterIterationPair.first); + } + } + } + + float tracksterTotalRawPt = t.raw_pt(); + std::vector haloTrackstersTRKIdx; + bool foundCompatibleTRK = false; + + for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { + usedTrackstersMerged[otherTracksterIdx] = true; + tracksterTotalRawPt += trackstersMergedHandle->at(otherTracksterIdx).raw_pt(); + + // Check the X,Y,Z barycenter and merge if they are very close (halo) + if ((t.barycenter() - trackstersMergedHandle->at(otherTracksterIdx).barycenter()).mag2() < + halo_max_distance2_) { + haloTrackstersTRKIdx.push_back(otherTracksterIdx); + + } else { + foundCompatibleTRK = true; + } + } + + //check if there is 1-to-1 relationship + if (trackstersTRKEMwithSameSeed.empty()) { + if (foundCompatibleTRK) { + TICLCandidate tmpCandidate; + tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, mergedIdx)); + double raw_energy = t.raw_energy(); + + tmpCandidate.setCharge(track.charge()); + tmpCandidate.setTrackPtr(edm::Ptr(track_h, trackIdx)); + tmpCandidate.setPdgId(211 * track.charge()); + for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { + tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); + raw_energy += trackstersMergedHandle->at(otherTracksterIdx).raw_energy(); + } + tmpCandidate.setRawEnergy(raw_energy); + math::XYZTLorentzVector p4(raw_energy * track.momentum().unit().x(), + raw_energy * track.momentum().unit().y(), + raw_energy * track.momentum().unit().z(), + raw_energy); + tmpCandidate.setP4(p4); + resultCandidates->push_back(tmpCandidate); + + } else { + TICLCandidate tmpCandidate; + tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, mergedIdx)); + double raw_energy = t.raw_energy(); + tmpCandidate.setCharge(track.charge()); + tmpCandidate.setTrackPtr(edm::Ptr(track_h, trackIdx)); + for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { + tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); + raw_energy += trackstersMergedHandle->at(otherTracksterIdx).raw_energy(); + } + tmpCandidate.setPdgId(11 * track.charge()); + + tmpCandidate.setRawEnergy(raw_energy); + math::XYZTLorentzVector p4(raw_energy * track.momentum().unit().x(), + raw_energy * track.momentum().unit().y(), + raw_energy * track.momentum().unit().z(), + raw_energy); + tmpCandidate.setP4(p4); + resultCandidates->push_back(tmpCandidate); + } + + } else { + // if 1-to-many find closest trackster in momentum + int closestTrackster = mergedIdx; + float minPtDiff = std::abs(t.raw_pt() - track.pt()); + for (auto otherTracksterIdx : trackstersTRKEMwithSameSeed) { + auto thisPt = tracksterTotalRawPt + trackstersMergedHandle->at(otherTracksterIdx).raw_pt() - t.raw_pt(); + closestTrackster = std::abs(thisPt - track.pt()) < minPtDiff ? otherTracksterIdx : closestTrackster; + } + usedTrackstersMerged[closestTrackster] = true; + + if (foundCompatibleTRK) { + TICLCandidate tmpCandidate; + tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, closestTrackster)); + double raw_energy = trackstersMergedHandle->at(closestTrackster).raw_energy(); + + tmpCandidate.setCharge(track.charge()); + tmpCandidate.setTrackPtr(edm::Ptr(track_h, trackIdx)); + tmpCandidate.setPdgId(211 * track.charge()); + for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { + tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); + raw_energy += trackstersMergedHandle->at(otherTracksterIdx).raw_energy(); + } + tmpCandidate.setRawEnergy(raw_energy); + float momentum = std::sqrt(raw_energy * raw_energy - mpion2); + math::XYZTLorentzVector p4(momentum * track.momentum().unit().x(), + momentum * track.momentum().unit().y(), + momentum * track.momentum().unit().z(), + raw_energy); + tmpCandidate.setP4(p4); + resultCandidates->push_back(tmpCandidate); + + } else { + TICLCandidate tmpCandidate; + tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, closestTrackster)); + double raw_energy = trackstersMergedHandle->at(closestTrackster).raw_energy(); + + tmpCandidate.setCharge(track.charge()); + tmpCandidate.setTrackPtr(edm::Ptr(track_h, trackIdx)); + for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { + tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); + raw_energy += trackstersMergedHandle->at(otherTracksterIdx).raw_energy(); + } + tmpCandidate.setPdgId(11 * track.charge()); + tmpCandidate.setRawEnergy(raw_energy); + math::XYZTLorentzVector p4(raw_energy * track.momentum().unit().x(), + raw_energy * track.momentum().unit().y(), + raw_energy * track.momentum().unit().z(), + raw_energy); + tmpCandidate.setP4(p4); + resultCandidates->push_back(tmpCandidate); + } + // Promote all other TRKEM tracksters as photons with their energy. + for (auto otherTracksterIdx : trackstersTRKEMwithSameSeed) { + auto tmpIndex = (otherTracksterIdx != closestTrackster) ? otherTracksterIdx : mergedIdx; + TICLCandidate photonCandidate; + const auto &otherTrackster = trackstersMergedHandle->at(tmpIndex); + auto gammaEnergy = otherTrackster.raw_energy(); + photonCandidate.setCharge(0); + photonCandidate.setPdgId(22); + photonCandidate.setRawEnergy(gammaEnergy); + math::XYZTLorentzVector gammaP4(gammaEnergy * otherTrackster.barycenter().unit().x(), + gammaEnergy * otherTrackster.barycenter().unit().y(), + gammaEnergy * otherTrackster.barycenter().unit().z(), + gammaEnergy); + photonCandidate.setP4(gammaP4); + photonCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, tmpIndex)); + resultCandidates->push_back(photonCandidate); + } + } + } + } + } //end of loop over trackstersTRKEM + + for (unsigned i = 0; i < trackstersTRK.size(); ++i) { + auto mergedIdx = indexInMergedCollTRK[i]; + const auto &t = trackstersTRK[i]; //trackster + + if (!usedTrackstersMerged[mergedIdx] and t.raw_energy() > 0) { + auto trackIdx = t.seedIndex(); + auto const &track = tracks[trackIdx]; + if (!usedSeeds[trackIdx]) { + usedSeeds[trackIdx] = true; + usedTrackstersMerged[mergedIdx] = true; + TICLCandidate tmpCandidate; + tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, mergedIdx)); + tmpCandidate.setCharge(track.charge()); + tmpCandidate.setTrackPtr(edm::Ptr(track_h, trackIdx)); + tmpCandidate.setPdgId(211 * track.charge()); + tmpCandidate.setRawEnergy(t.raw_energy()); + float momentum = std::sqrt(t.raw_energy() * t.raw_energy() - mpion2); + math::XYZTLorentzVector p4(momentum * track.momentum().unit().x(), + momentum * track.momentum().unit().y(), + momentum * track.momentum().unit().z(), + t.raw_energy()); + tmpCandidate.setP4(p4); + resultCandidates->push_back(tmpCandidate); + } + } + } + // For all seeds that have 0-energy tracksters whose track is not marked as used, create a charged hadron with the track information. + for (auto const &s : seedingTrk) { + if (usedSeeds[s.index] == false) { + auto const &track = tracks[s.index]; + // emit a charged hadron + TICLCandidate tmpCandidate; + tmpCandidate.setCharge(track.charge()); + tmpCandidate.setTrackPtr(edm::Ptr(track_h, s.index)); + tmpCandidate.setPdgId(211 * track.charge()); + float energy = std::sqrt(track.p() * track.p() + mpion2); + tmpCandidate.setRawEnergy(energy); + math::PtEtaPhiMLorentzVector p4Polar(track.pt(), track.eta(), track.phi(), mpion); + tmpCandidate.setP4(p4Polar); + resultCandidates->push_back(tmpCandidate); + usedSeeds[s.index] = true; + } + } + + // for all general tracks (high purity, pt > 1), check if they have been used: if not, promote them as charged hadrons + for (unsigned i = 0; i < tracks.size(); ++i) { + auto const &track = tracks[i]; + if (track.pt() > track_min_pt_ and track.quality(reco::TrackBase::highPurity) and + track.missingOuterHits() < track_max_missing_outerhits_ and std::abs(track.outerEta()) > track_min_eta_ and + std::abs(track.outerEta()) < track_max_eta_ and usedSeeds[i] == false) { + // emit a charged hadron + TICLCandidate tmpCandidate; + tmpCandidate.setCharge(track.charge()); + tmpCandidate.setTrackPtr(edm::Ptr(track_h, i)); + tmpCandidate.setPdgId(211 * track.charge()); + float energy = std::sqrt(track.p() * track.p() + mpion2); + tmpCandidate.setRawEnergy(energy); + math::PtEtaPhiMLorentzVector p4Polar(track.pt(), track.eta(), track.phi(), mpion); + tmpCandidate.setP4(p4Polar); + resultCandidates->push_back(tmpCandidate); + usedSeeds[i] = true; + } + } + + // Compute timing + assignTimeToCandidates(*resultCandidates); + + evt.put(std::move(resultCandidates)); +} + +void TrackstersMergeProducerV3::energyRegressionAndID(const std::vector &layerClusters, + const tensorflow::Session *eidSession, + std::vector &tracksters) const { + // Energy regression and particle identification strategy: + // + // 1. Set default values for regressed energy and particle id for each trackster. + // 2. Store indices of tracksters whose total sum of cluster energies is above the + // eidMinClusterEnergy_ (GeV) treshold. Inference is not applied for soft tracksters. + // 3. When no trackster passes the selection, return. + // 4. Create input and output tensors. The batch dimension is determined by the number of + // selected tracksters. + // 5. Fill input tensors with layer cluster features. Per layer, clusters are ordered descending + // by energy. Given that tensor data is contiguous in memory, we can use pointer arithmetic to + // fill values, even with batching. + // 6. Zero-fill features for empty clusters in each layer. + // 7. Batched inference. + // 8. Assign the regressed energy and id probabilities to each trackster. + // + // Indices used throughout this method: + // i -> batch element / trackster + // j -> layer + // k -> cluster + // l -> feature + + // set default values per trackster, determine if the cluster energy threshold is passed, + // and store indices of hard tracksters + std::vector tracksterIndices; + for (int i = 0; i < (int)tracksters.size(); i++) { + // calculate the cluster energy sum (2) + // note: after the loop, sumClusterEnergy might be just above the threshold + // which is enough to decide whether to run inference for the trackster or + // not + float sumClusterEnergy = 0.; + for (const unsigned int &vertex : tracksters[i].vertices()) { + sumClusterEnergy += (float)layerClusters[vertex].energy(); + // there might be many clusters, so try to stop early + if (sumClusterEnergy >= eidMinClusterEnergy_) { + // set default values (1) + tracksters[i].setRegressedEnergy(0.f); + tracksters[i].zeroProbabilities(); + tracksterIndices.push_back(i); + break; + } + } + } + + // do nothing when no trackster passes the selection (3) + int batchSize = (int)tracksterIndices.size(); + if (batchSize == 0) { + return; + } + + // create input and output tensors (4) + tensorflow::TensorShape shape({batchSize, eidNLayers_, eidNClusters_, eidNFeatures_}); + tensorflow::Tensor input(tensorflow::DT_FLOAT, shape); + tensorflow::NamedTensorList inputList = {{eidInputName_, input}}; + static constexpr int inputDimension = 4; + + std::vector outputs; + std::vector outputNames; + if (!eidOutputNameEnergy_.empty()) { + outputNames.push_back(eidOutputNameEnergy_); + } + if (!eidOutputNameId_.empty()) { + outputNames.push_back(eidOutputNameId_); + } + + // fill input tensor (5) + for (int i = 0; i < batchSize; i++) { + const Trackster &trackster = tracksters[tracksterIndices[i]]; + + // per layer, we only consider the first eidNClusters_ clusters in terms of + // energy, so in order to avoid creating large / nested structures to do + // the sorting for an unknown number of total clusters, create a sorted + // list of layer cluster indices to keep track of the filled clusters + std::vector clusterIndices(trackster.vertices().size()); + for (int k = 0; k < (int)trackster.vertices().size(); k++) { + clusterIndices[k] = k; + } + sort(clusterIndices.begin(), clusterIndices.end(), [&layerClusters, &trackster](const int &a, const int &b) { + return layerClusters[trackster.vertices(a)].energy() > layerClusters[trackster.vertices(b)].energy(); + }); + + // keep track of the number of seen clusters per layer + std::vector seenClusters(eidNLayers_); + + // loop through clusters by descending energy + for (const int &k : clusterIndices) { + // get features per layer and cluster and store the values directly in the input tensor + const reco::CaloCluster &cluster = layerClusters[trackster.vertices(k)]; + int j = rhtools_.getLayerWithOffset(cluster.hitsAndFractions()[0].first) - 1; + if (j < eidNLayers_ && seenClusters[j] < eidNClusters_) { + // get the pointer to the first feature value for the current batch, layer and cluster + float *features = &input.tensor()(i, j, seenClusters[j], 0); + + // fill features + *(features++) = float(cluster.energy() / float(trackster.vertex_multiplicity(k))); + *(features++) = float(std::abs(cluster.eta())); + *(features) = float(cluster.phi()); + + // increment seen clusters + seenClusters[j]++; + } + } + + // zero-fill features of empty clusters in each layer (6) + for (int j = 0; j < eidNLayers_; j++) { + for (int k = seenClusters[j]; k < eidNClusters_; k++) { + float *features = &input.tensor()(i, j, k, 0); + for (int l = 0; l < eidNFeatures_; l++) { + *(features++) = 0.f; + } + } + } + } + + // run the inference (7) + tensorflow::run(const_cast(eidSession), inputList, outputNames, &outputs); + + // store regressed energy per trackster (8) + if (!eidOutputNameEnergy_.empty()) { + // get the pointer to the energy tensor, dimension is batch x 1 + float *energy = outputs[0].flat().data(); + + for (const int &i : tracksterIndices) { + tracksters[i].setRegressedEnergy(*(energy++)); + } + } + + // store id probabilities per trackster (8) + if (!eidOutputNameId_.empty()) { + // get the pointer to the id probability tensor, dimension is batch x id_probabilities.size() + int probsIdx = !eidOutputNameEnergy_.empty(); + float *probs = outputs[probsIdx].flat().data(); + + for (const int &i : tracksterIndices) { + tracksters[i].setProbabilities(probs); + probs += tracksters[i].id_probabilities().size(); + } + } +} + +void TrackstersMergeProducerV3::assignTimeToCandidates(std::vector &resultCandidates) const { + for (auto &cand : resultCandidates) { + if (cand.tracksters().size() > 1) { // For single-trackster candidates the timing is already set + float time = 0.f; + float timeErr = 0.f; + for (const auto &tr : cand.tracksters()) { + if (tr->timeError() > 0) { + auto invTimeESq = pow(tr->timeError(), -2); + time += tr->time() * invTimeESq; + timeErr += invTimeESq; + } + } + if (timeErr > 0) { + timeErr = 1. / timeErr; + + cand.setTime(time * timeErr); + cand.setTimeError(sqrt(timeErr)); + } + } + } +} + +void TrackstersMergeProducerV3::printTrackstersDebug(const std::vector &tracksters, + const char *label) const { +#ifdef EDM_ML_DEBUG + int counter = 0; + for (auto const &t : tracksters) { + LogDebug("TrackstersMergeProducerV3") + << counter++ << " TrackstersMergeProducerV3 (" << label << ") obj barycenter: " << t.barycenter() + << " eta,phi (baricenter): " << t.barycenter().eta() << ", " << t.barycenter().phi() + << " eta,phi (eigen): " << t.eigenvectors(0).eta() << ", " << t.eigenvectors(0).phi() + << " pt(eigen): " << std::sqrt(t.eigenvectors(0).Unit().perp2()) * t.raw_energy() << " seedID: " << t.seedID() + << " seedIndex: " << t.seedIndex() << " size: " << t.vertices().size() << " average usage: " + << (std::accumulate(std::begin(t.vertex_multiplicity()), std::end(t.vertex_multiplicity()), 0.) / + (float)t.vertex_multiplicity().size()) + << " raw_energy: " << t.raw_energy() << " regressed energy: " << t.regressed_energy() + << " probs(ga/e/mu/np/cp/nh/am/unk): "; + for (auto const &p : t.id_probabilities()) { + LogDebug("TrackstersMergeProducerV3") << std::fixed << p << " "; + } + LogDebug("TrackstersMergeProducerV3") << " sigmas: "; + for (auto const &s : t.sigmas()) { + LogDebug("TrackstersMergeProducerV3") << s << " "; + } + LogDebug("TrackstersMergeProducerV3") << std::endl; + } +#endif +} + +void TrackstersMergeProducerV3::fillDescriptions(edm::ConfigurationDescriptions &descriptions) { + edm::ParameterSetDescription desc; + desc.add("tracksterstrkem", edm::InputTag("ticlTrackstersTrkEM")); + desc.add("trackstersem", edm::InputTag("ticlTrackstersEM")); + desc.add("tracksterstrk", edm::InputTag("ticlTrackstersTrk")); + desc.add("trackstershad", edm::InputTag("ticlTrackstersHAD")); + desc.add("seedingTrk", edm::InputTag("ticlSeedingTrk")); + desc.add("layer_clusters", edm::InputTag("hgcalLayerClusters")); + desc.add("layer_clustersTime", edm::InputTag("hgcalLayerClusters", "timeLayerCluster")); + desc.add("tracks", edm::InputTag("generalTracks")); + desc.add("optimiseAcrossTracksters", true); + desc.add("eta_bin_window", 1); + desc.add("phi_bin_window", 1); + desc.add("pt_sigma_high", 2.); + desc.add("pt_sigma_low", 2.); + desc.add("halo_max_distance2", 4.); + desc.add("track_min_pt", 1.); + desc.add("track_min_eta", 1.48); + desc.add("track_max_eta", 3.); + desc.add("track_max_missing_outerhits", 5); + desc.add("cosangle_align", 0.9945); + desc.add("e_over_h_threshold", 1.); + desc.add("pt_neutral_threshold", 2.); + desc.add("resol_calo_offset_had", 1.5); + desc.add("resol_calo_scale_had", 0.15); + desc.add("resol_calo_offset_em", 1.5); + desc.add("resol_calo_scale_em", 0.15); + desc.add("tfDnnLabel", "tracksterSelectionTf"); + desc.add("eid_input_name", "input"); + desc.add("eid_output_name_energy", "output/regressed_energy"); + desc.add("eid_output_name_id", "output/id_probabilities"); + desc.add("eid_min_cluster_energy", 1.); + desc.add("eid_n_layers", 50); + desc.add("eid_n_clusters", 10); + descriptions.add("trackstersMergeProducerV3", desc); +} + +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(TrackstersMergeProducerV3); diff --git a/RecoHGCal/TICL/python/iterativeTICL_cff.py b/RecoHGCal/TICL/python/iterativeTICL_cff.py index 66c6bd3a1959e..69c9989ca8955 100644 --- a/RecoHGCal/TICL/python/iterativeTICL_cff.py +++ b/RecoHGCal/TICL/python/iterativeTICL_cff.py @@ -12,22 +12,19 @@ from RecoHGCal.TICL.ticlLayerTileProducer_cfi import ticlLayerTileProducer from RecoHGCal.TICL.pfTICLProducer_cfi import pfTICLProducer as _pfTICLProducer from RecoHGCal.TICL.trackstersMergeProducer_cfi import trackstersMergeProducer as _trackstersMergeProducer +from RecoHGCal.TICL.trackstersMergeProducerV3_cfi import trackstersMergeProducerV3 as _trackstersMergeProducerV3 from RecoHGCal.TICL.tracksterSelectionTf_cfi import * ticlLayerTileTask = cms.Task(ticlLayerTileProducer) ticlTrackstersMerge = _trackstersMergeProducer.clone() -ticlTracksterMergeTask = cms.Task(ticlTrackstersMerge) - +ticlTrackstersMergeV3 = _trackstersMergeProducerV3.clone() pfTICL = _pfTICLProducer.clone() ticlPFTask = cms.Task(pfTICL) ticlIterationsTask = cms.Task( - ticlTrkEMStepTask - ,ticlEMStepTask - ,ticlTrkStepTask - ,ticlHADStepTask + ticlCLUE3DHighStepTask ) from Configuration.ProcessModifiers.clue3D_cff import clue3D @@ -36,15 +33,32 @@ from Configuration.ProcessModifiers.fastJetTICL_cff import fastJetTICL fastJetTICL.toModify(ticlIterationsTask, func=lambda x : x.add(ticlFastJetStepTask)) +from Configuration.ProcessModifiers.ticl_v3_cff import ticl_v3 +ticl_v3.toModify(ticlIterationsTask, func=lambda x : x.add( ticlTrkEMStepTask + ,ticlEMStepTask + ,ticlTrkStepTask + ,ticlHADStepTask) ) ticlIterLabels = [_step.itername.value() for _iteration in ticlIterationsTask for _step in _iteration if (_step._TypedParameterizable__type == "TrackstersProducer")] -iterTICLTask = cms.Task(ticlLayerTileTask +ticlTracksterMergeTask = cms.Task(ticlTrackstersMerge) +ticlTracksterMergeTaskV3 = cms.Task(ticlTrackstersMergeV3) + +ticl_v3.toModify(pfTICL, ticlCandidateSrc = "ticlTrackstersMergeV3") + +mergeTICLTask = cms.Task(ticlLayerTileTask ,ticlIterationsTask ,ticlTracksterMergeTask - ,ticlPFTask ) + +ticl_v3.toModify(mergeTICLTask, func=lambda x : x.add(ticlTracksterMergeTaskV3)) ticlIterLabelsMerge = ticlIterLabels + ["Merge"] +ticlIterLabelsMergeV3 = ticlIterLabels + ["MergeV3"] +ticl_v3.toModify(ticlIterLabelsMerge, func=lambda x : x.extend(ticlIterLabelsMergeV3)) + +iterTICLTask = cms.Task(mergeTICLTask + ,ticlPFTask) + ticlLayerTileHFNose = ticlLayerTileProducer.clone( detector = 'HFNose' ) diff --git a/Validation/HGCalValidation/python/hgcalPlots.py b/Validation/HGCalValidation/python/hgcalPlots.py index 76e3194541726..50dc4bc37cb0e 100644 --- a/Validation/HGCalValidation/python/hgcalPlots.py +++ b/Validation/HGCalValidation/python/hgcalPlots.py @@ -2495,7 +2495,7 @@ def append_hgcalTrackstersPlots(collection = 'ticlTrackstersMerge', name_collect loopSubFolders=False, purpose=PlotPurpose.Timing #,page=tsToCP_linking.replace('TSToCP_','TICL-') - ,page=tsToCP_linking.replace('TSToCP_','Test-TICL').replace('linking','') + ,page=tsToCP_linking.replace('TSToCP_','TICL-linking').replace('linking','') ,section=name_collection) )