From e72f7828ea1e206c23bac028c5af8dacc5eaa387 Mon Sep 17 00:00:00 2001 From: Felice Pantaleo Date: Thu, 27 Jul 2023 18:27:24 +0200 Subject: [PATCH 01/12] Improve MC truth information in TICL + Associators + TICLDumper - Add best associated RecoTrack to SimTracksters, by using SimTrack-TrackingParticle associator and TrackingParticle-RecoTrack associator. Add SimTrackster from PU by considering as PU all the remaining contribution from the LCs after the SimTracksterCP/SC building - Creates SimTICLCandidates, object represeting the best PF reconstruction possible with the available recoObjects. It contains all the SimTracksters beloning to the same CP and the recoTrack associated to the CP - Associators: compute score between SimTracksters and RecoTrackster - TICLDumper: EDAnalyzer to save TICL information in root file Co-authored-by: Wahid Redjeb Co-authored-by: Aurora Perego --- .../HGCalReco/interface/TICLCandidate.h | 9 + DataFormats/HGCalReco/interface/Trackster.h | 76 +- DataFormats/HGCalReco/src/classes_def.xml | 3 +- .../python/RecoHGCal_EventContent_cff.py | 15 +- RecoHGCal/TICL/interface/commons.h | 54 + RecoHGCal/TICL/plugins/BuildFile.xml | 2 + .../LinkingAlgoByDirectionGeometric.cc | 9 +- .../plugins/LinkingAlgoByDirectionGeometric.h | 11 +- .../TICL/plugins/SimTrackstersProducer.cc | 313 ++- RecoHGCal/TICL/plugins/TICLDumper.cc | 2102 +++++++++++++++++ .../TICL/plugins/TrackstersMergeProducer.cc | 17 +- .../TICL/python/customiseTICLFromReco.py | 41 +- .../TSToSimTSAssociatorByEnergyScoreImpl.cc | 18 +- .../TSToSimTSAssociatorByEnergyScoreImpl.h | 2 +- .../plugins/TSToSimTSAssociatorEDProducer.cc | 8 +- ...ToSimTSHitLCAssociatorByEnergyScoreImpl.cc | 259 ++ ...SToSimTSHitLCAssociatorByEnergyScoreImpl.h | 78 + ...mTSHitLCAssociatorByEnergyScoreProducer.cc | 67 + .../TSToSimTSHitLCAssociatorEDProducer.cc | 103 + .../python/TSToSimTSAssociation_cfi.py | 56 +- ...racksterToSimTracksterAssociatorBaseImpl.h | 2 +- .../TracksterToSimTracksterHitLCAssociator.h | 57 + ...terToSimTracksterHitLCAssociatorBaseImpl.h | 58 + .../TracksterToSimTracksterHitLCAssociator.cc | 5 + ...erToSimTracksterHitLCAssociatorBaseImpl.cc | 34 + SimDataFormats/Associations/src/classes.h | 5 +- .../Associations/src/classes_def.xml | 3 + SimDataFormats/Track/BuildFile.xml | 1 + .../Track/interface/UniqueSimTrackId.h | 6 +- SimDataFormats/Track/src/classes.h | 1 + SimDataFormats/Track/src/classes_def.xml | 7 + .../plugins/SimHitTPAssociationProducer.cc | 10 +- .../python/SimTracker_EventContent_cff.py | 3 +- .../Configuration/python/hgcalSimValid_cff.py | 11 + 34 files changed, 3324 insertions(+), 122 deletions(-) create mode 100644 RecoHGCal/TICL/plugins/TICLDumper.cc create mode 100644 SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorByEnergyScoreImpl.cc create mode 100644 SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorByEnergyScoreImpl.h create mode 100644 SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorByEnergyScoreProducer.cc create mode 100644 SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorEDProducer.cc create mode 100644 SimDataFormats/Associations/interface/TracksterToSimTracksterHitLCAssociator.h create mode 100644 SimDataFormats/Associations/interface/TracksterToSimTracksterHitLCAssociatorBaseImpl.h create mode 100644 SimDataFormats/Associations/src/TracksterToSimTracksterHitLCAssociator.cc create mode 100644 SimDataFormats/Associations/src/TracksterToSimTracksterHitLCAssociatorBaseImpl.cc diff --git a/DataFormats/HGCalReco/interface/TICLCandidate.h b/DataFormats/HGCalReco/interface/TICLCandidate.h index ab9b59fdab6ee..25cedc49df025 100644 --- a/DataFormats/HGCalReco/interface/TICLCandidate.h +++ b/DataFormats/HGCalReco/interface/TICLCandidate.h @@ -53,7 +53,16 @@ class TICLCandidate : public reco::LeafCandidate { return idProbabilities_[(int)type]; } + inline const std::array& idProbabilities() const { return idProbabilities_; } + + void zeroProbabilities() { + for (auto& p : idProbabilities_) { + p = 0.f; + } + } + void setIdProbabilities(const std::array& idProbs) { idProbabilities_ = idProbs; } + inline void setIdProbability(ParticleType type, float value) { idProbabilities_[int(type)] = value; } private: float time_; diff --git a/DataFormats/HGCalReco/interface/Trackster.h b/DataFormats/HGCalReco/interface/Trackster.h index 242b1080272dc..648442584f516 100644 --- a/DataFormats/HGCalReco/interface/Trackster.h +++ b/DataFormats/HGCalReco/interface/Trackster.h @@ -18,7 +18,7 @@ namespace ticl { class Trackster { public: - typedef math::XYZVector Vector; + typedef math::XYZVectorF Vector; enum IterationIndex { TRKEM = 0, EM, TRKHAD, HAD, MIP, SIM, SIM_CP }; @@ -37,21 +37,20 @@ namespace ticl { enum class PCAOrdering { ascending = 0, descending }; Trackster() - : iterationIndex_(0), - seedIndex_(-1), - time_(0.f), - timeError_(-1.f), + : barycenter_({0.f, 0.f, 0.f}), regressed_energy_(0.f), raw_energy_(0.f), + time_(0.f), + timeError_(-1.f), raw_em_energy_(0.f), + id_probabilities_{}, raw_pt_(0.f), raw_em_pt_(0.f), - barycenter_({0., 0., 0.}), - eigenvalues_{{0.f, 0.f, 0.f}}, - sigmas_{{0.f, 0.f, 0.f}}, - sigmasPCA_{{0.f, 0.f, 0.f}} { - zeroProbabilities(); - } + seedIndex_(-1), + eigenvalues_{}, + sigmas_{}, + sigmasPCA_{}, + iterationIndex_(0) {} inline void setIteration(const Trackster::IterationIndex i) { iterationIndex_ = i; } std::vector &vertices() { return vertices_; } @@ -73,6 +72,9 @@ namespace ticl { inline void setRawPt(float value) { raw_pt_ = value; } inline void setRawEmPt(float value) { raw_em_pt_ = value; } inline void setBarycenter(Vector value) { barycenter_ = value; } + inline void setTrackIdx(int index) { track_idx_ = index; } + int trackIdx() const { return track_idx_; } + inline void fillPCAVariables(Eigen::Vector3d &eigenvalues, Eigen::Matrix3d &eigenvectors, Eigen::Vector3d &sigmas, @@ -147,20 +149,23 @@ namespace ticl { } private: - // TICL iteration producing the trackster - uint8_t iterationIndex_; + Vector barycenter_; + float regressed_energy_; + float raw_energy_; + // -99, -1 if not available. ns units otherwise + float time_; + float timeError_; + float raw_em_energy_; + + // trackster ID probabilities + std::array id_probabilities_; // The vertices of the DAG are the indices of the // 2d objects in the global collection std::vector vertices_; std::vector vertex_multiplicity_; - - // The edges connect two vertices together in a directed doublet - // ATTENTION: order matters! - // A doublet generator should create edges in which: - // the first element is on the inner layer and - // the outer element is on the outer layer. - std::vector > edges_; + float raw_pt_; + float raw_em_pt_; // Product ID of the seeding collection used to create the Trackster. // For GlobalSeeding the ProductID is set to 0. For track-based seeding @@ -173,31 +178,22 @@ namespace ticl { // created the trackster. For track-based seeding the pointer to the track // can be cooked using the previous ProductID and this index. int seedIndex_; + int track_idx_ = -1; - // We also need the pointer to the original seeding region ?? - // something like: - // int seedingRegionIdx; - - // -99, -1 if not available. ns units otherwise - float time_; - float timeError_; - - // regressed energy - float regressed_energy_; - float raw_energy_; - float raw_em_energy_; - float raw_pt_; - float raw_em_pt_; - - // PCA Variables - Vector barycenter_; - std::array eigenvalues_; std::array eigenvectors_; + std::array eigenvalues_; std::array sigmas_; std::array sigmasPCA_; - // trackster ID probabilities - std::array id_probabilities_; + // The edges connect two vertices together in a directed doublet + // ATTENTION: order matters! + // A doublet generator should create edges in which: + // the first element is on the inner layer and + // the outer element is on the outer layer. + std::vector > edges_; + + // TICL iteration producing the trackster + uint8_t iterationIndex_; }; typedef std::vector TracksterCollection; diff --git a/DataFormats/HGCalReco/src/classes_def.xml b/DataFormats/HGCalReco/src/classes_def.xml index 9ac5b29b4f2e2..18c2cf05b934d 100644 --- a/DataFormats/HGCalReco/src/classes_def.xml +++ b/DataFormats/HGCalReco/src/classes_def.xml @@ -1,6 +1,7 @@ - + + diff --git a/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py b/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py index bcf826da90773..867c8a52e9b09 100644 --- a/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py +++ b/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py @@ -19,7 +19,12 @@ 'keep *_ticlTrackstersHFNoseMIP_*_*', 'keep *_ticlTrackstersHFNoseHAD_*_*', 'keep *_ticlTrackstersHFNoseMerge_*_*',] + - ['keep *_pfTICL_*_*'] + ['keep *_pfTICL_*_*'] + + ['keep CaloParticles_mix_*_*', 'keep SimClusters_mix_*_*'] + + ['keep *_layerClusterSimClusterAssociationProducer_*_*','keep *_layerClusterCaloParticleAssociationProducer_*_*', 'keep *_layerClusterSimTracksterAssociationProducer_*_*'] + + ['keep *_tracksterSimTracksterAssociationLinking_*_*' ,'keep *_tracksterSimTracksterAssociationPR_*_*'] + + ['keep *_tracksterSimTracksterAssociationLinkingPU_*_*' ,'keep *_tracksterSimTracksterAssociationPRPU_*_*'] + + ['keep *_tracksterSimTracksterAssociationLinkingbyCLUE3D_*_*', 'keep *_tracksterSimTracksterAssociationPRbyCLUE3D_*_*'] ) ) TICL_RECO.outputCommands.extend(TICL_AOD.outputCommands) @@ -28,6 +33,7 @@ TICL_FEVT = cms.PSet( outputCommands = cms.untracked.vstring( 'keep *_ticlSimTracksters_*_*', + 'keep *_ticlSimTICLCandidates_*_*', 'keep *_ticlSimTrackstersFromCP_*_*', ) ) @@ -48,6 +54,13 @@ def cleanOutputAndSet(outputModule, ticl_outputCommads): 'keep *_layerClusterSimClusterAssociationProducer_*_*', 'keep *_layerClusterCaloParticleAssociationProducer_*_*', 'keep *_randomEngineStateProducer_*_*', + 'keep *_layerClusterSimTracksterAssociationProducer_*_*', + 'keep *_tracksterSimTracksterAssociationLinking_*_*', + 'keep *_tracksterSimTracksterAssociationPR_*_*', + 'keep *_tracksterSimTracksterAssociationLinkingPU_*_*', + 'keep *_tracksterSimTracksterAssociationPRPU_*_*', + 'keep *_tracksterSimTracksterAssociationLinkingbyCLUE3D_*_*', + 'keep *_tracksterSimTracksterAssociationPRbyCLUE3D_*_*', ]) if hasattr(process, 'FEVTDEBUGEventContent'): diff --git a/RecoHGCal/TICL/interface/commons.h b/RecoHGCal/TICL/interface/commons.h index 1ec1956cee38a..86b92f4c925fe 100644 --- a/RecoHGCal/TICL/interface/commons.h +++ b/RecoHGCal/TICL/interface/commons.h @@ -4,12 +4,28 @@ #include "DataFormats/Common/interface/Ref.h" #include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h" #include "DataFormats/HGCalReco/interface/Trackster.h" +#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h" namespace ticl { //constants constexpr double mpion = 0.13957; constexpr float mpion2 = mpion * mpion; + typedef math::XYZVectorF Vector; + enum LayerType { + + CE_E_120 = 0, + CE_E_200 = 1, + CE_E_300 = 2, + CE_H_120_F = 3, + CE_H_200_F = 4, + CE_H_300_F = 5, + CE_H_200_C = 6, + CE_H_300_C = 7, + CE_H_SCINT_C = 8, + EnumSize = 9 + + }; inline Trackster::ParticleType tracksterParticleTypeFromPdgId(int pdgId, int charge) { if (pdgId == 111) { @@ -40,6 +56,44 @@ namespace ticl { // verbosity levels for ticl algorithms enum VerbosityLevel { None = 0, Basic, Advanced, Expert, Guru }; + inline int returnClusterType(DetId& lc_seed, const hgcal::RecHitTools& rhtools_) { + auto layer_number = rhtools_.getLayerWithOffset(lc_seed); + auto thickness = rhtools_.getSiThickIndex(lc_seed); + auto isEELayer = (layer_number <= rhtools_.lastLayerEE(false)); + auto isScintillator = rhtools_.isScintillator(lc_seed); + auto isFine = (layer_number <= rhtools_.lastLayerEE(false) + 7); + + if (isScintillator) { + return CE_H_SCINT_C; + } + if (isEELayer) { + if (thickness == 0) { + return CE_E_120; + } else if (thickness == 1) { + return CE_E_200; + } else if (thickness == 2) { + return CE_E_300; + } + } else { + if (isFine) { + if (thickness == 0) { + return CE_H_120_F; + } else if (thickness == 1) { + return CE_H_200_F; + } else if (thickness == 2) { + return CE_H_300_F; + } + } else { + if (thickness == 1) { + return CE_H_200_C; + } else if (thickness == 2) { + return CE_H_300_C; + } + } + } + return -1; + }; + } // namespace ticl #endif diff --git a/RecoHGCal/TICL/plugins/BuildFile.xml b/RecoHGCal/TICL/plugins/BuildFile.xml index 8c631381e21b4..b15bebb41e80b 100644 --- a/RecoHGCal/TICL/plugins/BuildFile.xml +++ b/RecoHGCal/TICL/plugins/BuildFile.xml @@ -29,6 +29,8 @@ + + diff --git a/RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.cc b/RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.cc index 1b23bf8154a93..bc90431b5264d 100644 --- a/RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.cc +++ b/RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.cc @@ -10,6 +10,7 @@ #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h" #include "RecoParticleFlow/PFProducer/interface/PFMuonAlgo.h" +#include "DataFormats/HGCalReco/interface/Trackster.h" using namespace ticl; @@ -36,10 +37,10 @@ void LinkingAlgoByDirectionGeometric::initialize(const HGCalDDDConstants *hgcons propagator_ = propH; } -math::XYZVector LinkingAlgoByDirectionGeometric::propagateTrackster(const Trackster &t, - const unsigned idx, - float zVal, - std::array &tracksterTiles) { +Vector 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(); diff --git a/RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.h b/RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.h index 82139a8f1aa5c..6598a852a1f1c 100644 --- a/RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.h +++ b/RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.h @@ -20,6 +20,7 @@ #include "CommonTools/Utils/interface/StringCutObjectSelector.h" #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" +#include "DataFormats/HGCalReco/interface/Trackster.h" namespace ticl { class LinkingAlgoByDirectionGeometric final : public LinkingAlgoBase { @@ -44,14 +45,14 @@ namespace ticl { static void fillPSetDescription(edm::ParameterSetDescription &desc); private: - typedef math::XYZVector Vector; + using Vector = ticl::Trackster::Vector; void buildLayers(); - math::XYZVector propagateTrackster(const Trackster &t, - const unsigned idx, - float zVal, - std::array &tracksterTiles); + Vector propagateTrackster(const Trackster &t, + const unsigned idx, + float zVal, + std::array &tracksterTiles); void findTrackstersInWindow(const std::vector> &seedingCollection, const std::array &tracksterTiles, diff --git a/RecoHGCal/TICL/plugins/SimTrackstersProducer.cc b/RecoHGCal/TICL/plugins/SimTrackstersProducer.cc index 4e479a908ec96..34eae24cb0c07 100644 --- a/RecoHGCal/TICL/plugins/SimTrackstersProducer.cc +++ b/RecoHGCal/TICL/plugins/SimTrackstersProducer.cc @@ -11,11 +11,16 @@ #include "FWCore/MessageLogger/interface/MessageLogger.h" #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" #include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "DataFormats/Common/interface/OrphanHandle.h" #include "DataFormats/CaloRecHit/interface/CaloCluster.h" #include "DataFormats/ParticleFlowReco/interface/PFCluster.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 "CommonTools/Utils/interface/StringCutObjectSelector.h" #include "DataFormats/Common/interface/ValueMap.h" #include "SimDataFormats/Associations/interface/LayerClusterToSimClusterAssociator.h" @@ -25,6 +30,12 @@ #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h" #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h" +#include "SimDataFormats/Track/interface/UniqueSimTrackId.h" +#include "SimDataFormats/Vertex/interface/SimVertex.h" + +#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h" +#include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h" + #include "RecoHGCal/TICL/interface/commons.h" #include "TrackstersPCA.h" @@ -32,6 +43,7 @@ #include #include #include +#include using namespace ticl; @@ -41,17 +53,26 @@ 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, + void makePUTrackster(const std::vector& inputClusterMask, + std::vector& output_mask, + std::vector& result, + const edm::ProductID seed, + int loop_index); + + 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 float fractionCut_, + const float energy, + const int pdgId, + const int charge, + float time, + const edm::ProductID seed, const Trackster::IterationIndex iter, std::vector& output_mask, - std::vector& result); + std::vector& result, + int& loop_index, + const bool add = false); private: std::string detector_; @@ -67,8 +88,19 @@ class SimTrackstersProducer : public edm::stream::EDProducer<> { const edm::EDGetTokenT associatorMapCaloParticleToReco_token_; const edm::ESGetToken geom_token_; hgcal::RecHitTools rhtools_; - const double fractionCut_; + const float fractionCut_; + const float qualityCutTrack_; + const edm::EDGetTokenT> trackingParticleToken_; + const edm::EDGetTokenT> simVerticesToken_; + + const edm::EDGetTokenT> recoTracksToken_; + const StringCutObjectSelector cutTk_; + + const edm::EDGetTokenT associatormapStRsToken_; + const edm::EDGetTokenT associatormapRtSsToken_; + const edm::EDGetTokenT associationSimTrackToTPToken_; }; + DEFINE_FWK_MODULE(SimTrackstersProducer); SimTrackstersProducer::SimTrackstersProducer(const edm::ParameterSet& ps) @@ -84,16 +116,25 @@ SimTrackstersProducer::SimTrackstersProducer(const edm::ParameterSet& ps) associatorMapCaloParticleToReco_token_( consumes(ps.getParameter("layerClusterCaloParticleAssociator"))), geom_token_(esConsumes()), - fractionCut_(ps.getParameter("fractionCut")) { + fractionCut_(ps.getParameter("fractionCut")), + qualityCutTrack_(ps.getParameter("qualityCutTrack")), + trackingParticleToken_( + consumes>(ps.getParameter("trackingParticles"))), + simVerticesToken_(consumes>(ps.getParameter("simVertices"))), + recoTracksToken_(consumes>(ps.getParameter("recoTracks"))), + cutTk_(ps.getParameter("cutTk")), + associatormapStRsToken_(consumes(ps.getParameter("tpToTrack"))), + associationSimTrackToTPToken_(consumes(ps.getParameter("simTrackToTPMap"))) { produces(); produces>(); produces("fromCPs"); + produces("PU"); produces>("fromCPs"); produces>>(); + produces>(); } void SimTrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { - // hgcalMultiClusters edm::ParameterSetDescription desc; desc.add("detector", "HGCAL"); desc.add("layer_clusters", edm::InputTag("hgcalMergeLayerClusters")); @@ -105,28 +146,59 @@ void SimTrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& des edm::InputTag("layerClusterSimClusterAssociationProducer")); desc.add("layerClusterCaloParticleAssociator", edm::InputTag("layerClusterCaloParticleAssociationProducer")); + desc.add("recoTracks", edm::InputTag("generalTracks")); + desc.add("cutTk", + "1.48 < abs(eta) < 3.0 && pt > 1. && quality(\"highPurity\") && " + "hitPattern().numberOfLostHits(\"MISSING_OUTER_HITS\") < 5"); + desc.add("tpToTrack", edm::InputTag("trackingParticleRecoTrackAsssociation")); + + desc.add("trackingParticles", edm::InputTag("mix", "MergedTrackTruth")); + desc.add("simVertices", edm::InputTag("g4SimHits")); + + desc.add("simTrackToTPMap", edm::InputTag("simHitTPAssocProducer", "simTrackToTP")); desc.add("fractionCut", 0.); + desc.add("qualityCutTrack", 0.75); descriptions.addWithDefaultLabel(desc); } +void SimTrackstersProducer::makePUTrackster(const std::vector& inputClusterMask, + std::vector& output_mask, + std::vector& result, + const edm::ProductID seed, + int loop_index) { + Trackster tmpTrackster; + for (size_t i = 0; i < output_mask.size(); i++) { + const float remaining_fraction = output_mask[i]; + if (remaining_fraction > 0.) { + tmpTrackster.vertices().push_back(i); + tmpTrackster.vertex_multiplicity().push_back(1. / remaining_fraction); + } + } + tmpTrackster.setSeed(seed, 0); + result.push_back(tmpTrackster); +} void SimTrackstersProducer::addTrackster( - const int& index, + 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 float fractionCut_, + const float energy, + const int pdgId, + const int charge, + const float time, + const edm::ProductID seed, const Trackster::IterationIndex iter, std::vector& output_mask, - std::vector& result) { - if (lcVec.empty()) + std::vector& result, + int& loop_index, + const bool add) { + Trackster tmpTrackster; + if (lcVec.empty()) { + result[index] = tmpTrackster; return; + } - Trackster tmpTrackster; - tmpTrackster.zeroProbabilities(); tmpTrackster.vertices().reserve(lcVec.size()); tmpTrackster.vertex_multiplicity().reserve(lcVec.size()); for (auto const& [lc, energyScorePair] : lcVec) { @@ -144,15 +216,23 @@ void SimTrackstersProducer::addTrackster( tmpTrackster.setRegressedEnergy(energy); tmpTrackster.setIteration(iter); tmpTrackster.setSeed(seed, index); - result.emplace_back(tmpTrackster); + if (add) { + result[index] = tmpTrackster; + loop_index += 1; + } else { + result.push_back(tmpTrackster); + } } void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) { auto result = std::make_unique(); auto output_mask = std::make_unique>(); auto result_fromCP = std::make_unique(); + auto resultPU = std::make_unique(); auto output_mask_fromCP = std::make_unique>(); auto cpToSc_SimTrackstersMap = std::make_unique>>(); + auto result_ticlCandidates = std::make_unique>(); + const auto& layerClusters = evt.get(clusters_token_); const auto& layerClustersTimes = evt.get(clustersTime_token_); const auto& inputClusterMask = evt.get(filtered_layerclusters_mask_token_); @@ -160,21 +240,38 @@ void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) output_mask_fromCP->resize(layerClusters.size(), 1.f); const auto& simclusters = evt.get(simclusters_token_); - const auto& caloparticles = evt.get(caloparticles_token_); + edm::Handle> caloParticles_h; + evt.getByToken(caloparticles_token_, caloParticles_h); + const auto& caloparticles = *caloParticles_h; const auto& simClustersToRecoColl = evt.get(associatorMapSimClusterToReco_token_); const auto& caloParticlesToRecoColl = evt.get(associatorMapCaloParticleToReco_token_); + const auto& simVertices = evt.get(simVerticesToken_); + + edm::Handle> trackingParticles_h; + evt.getByToken(trackingParticleToken_, trackingParticles_h); + edm::Handle> recoTracks_h; + evt.getByToken(recoTracksToken_, recoTracks_h); + const auto& TPtoRecoTrackMap = evt.get(associatormapStRsToken_); + const auto& simTrackToTPMap = evt.get(associationSimTrackToTPToken_); + const auto& recoTracks = *recoTracks_h; const auto& geom = es.getData(geom_token_); rhtools_.setGeometry(geom); const auto num_simclusters = simclusters.size(); result->reserve(num_simclusters); // Conservative size, will call shrink_to_fit later const auto num_caloparticles = caloparticles.size(); - result_fromCP->reserve(num_caloparticles); - + result_fromCP->resize(num_caloparticles); + std::map SimClusterToCaloParticleMap; + int loop_index = 0; for (const auto& [key, lcVec] : caloParticlesToRecoColl) { auto const& cp = *(key); auto cpIndex = &cp - &caloparticles[0]; + for (const auto& scRef : cp.simClusters()) { + auto const& sc = *(scRef); + auto const scIndex = &sc - &simclusters[0]; + SimClusterToCaloParticleMap[scIndex] = cpIndex; + } auto regr_energy = cp.energy(); std::vector scSimTracksterIdx; @@ -183,7 +280,7 @@ void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) // Create a Trackster from the object entering HGCal if (cp.g4Tracks()[0].crossedBoundary()) { regr_energy = cp.g4Tracks()[0].getMomentumAtBoundary().energy(); - + float time = cp.g4Tracks()[0].getPositionAtBoundary().t(); addTrackster(cpIndex, lcVec, inputClusterMask, @@ -191,10 +288,12 @@ void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) regr_energy, cp.pdgId(), cp.charge(), + time, key.id(), ticl::Trackster::SIM, *output_mask, - *result); + *result, + loop_index); } else { for (const auto& scRef : cp.simClusters()) { const auto& it = simClustersToRecoColl.find(scRef); @@ -211,10 +310,12 @@ void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) sc.g4Tracks()[0].getMomentumAtBoundary().energy(), sc.pdgId(), sc.charge(), + sc.g4Tracks()[0].getPositionAtBoundary().t(), scRef.id(), ticl::Trackster::SIM, *output_mask, - *result); + *result, + loop_index); if (result->empty()) continue; @@ -225,7 +326,7 @@ void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) } scSimTracksterIdx.shrink_to_fit(); } - + float time = simVertices[cp.g4Tracks()[0].vertIndex()].position().t(); // Create a Trackster from any CP addTrackster(cpIndex, lcVec, @@ -234,29 +335,175 @@ void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) regr_energy, cp.pdgId(), cp.charge(), + time, key.id(), ticl::Trackster::SIM_CP, *output_mask_fromCP, - *result_fromCP); + *result_fromCP, + loop_index, + true); if (result_fromCP->empty()) continue; - const auto index = result_fromCP->size() - 1; + const auto index = loop_index - 1; if (cpToSc_SimTrackstersMap->find(index) == cpToSc_SimTrackstersMap->end()) { (*cpToSc_SimTrackstersMap)[index] = scSimTracksterIdx; } } - + // TODO: remove time computation from PCA calculation and + // store time from boundary position in simTracksters ticl::assignPCAtoTracksters( *result, layerClusters, layerClustersTimes, rhtools_.getPositionLayer(rhtools_.lastLayerEE(doNose_)).z()); result->shrink_to_fit(); ticl::assignPCAtoTracksters( *result_fromCP, layerClusters, layerClustersTimes, rhtools_.getPositionLayer(rhtools_.lastLayerEE(doNose_)).z()); - result_fromCP->shrink_to_fit(); - evt.put(std::move(result)); + makePUTrackster(inputClusterMask, *output_mask, *resultPU, caloParticles_h.id(), 0); + + auto simTrackToRecoTrack = [&](UniqueSimTrackId simTkId) -> std::pair { + int trackIdx = -1; + float quality = 0.f; + auto ipos = simTrackToTPMap.mapping.find(simTkId); + if (ipos != simTrackToTPMap.mapping.end()) { + auto jpos = TPtoRecoTrackMap.find((ipos->second)); + if (jpos != TPtoRecoTrackMap.end()) { + auto& associatedRecoTracks = jpos->val; + if (!associatedRecoTracks.empty()) { + // associated reco tracks are sorted by decreasing quality + if (associatedRecoTracks[0].second > qualityCutTrack_) { + trackIdx = &(*associatedRecoTracks[0].first) - &recoTracks[0]; + quality = associatedRecoTracks[0].second; + } + } + } + } + return {trackIdx, quality}; + }; + + // Creating the map from TrackingParticle to SimTrackstersFromCP + auto& simTrackstersFromCP = *result_fromCP; + for (unsigned int i = 0; i < simTrackstersFromCP.size(); ++i) { + if (simTrackstersFromCP[i].vertices().empty()) + continue; + const auto& simTrack = caloparticles[simTrackstersFromCP[i].seedIndex()].g4Tracks()[0]; + UniqueSimTrackId simTkIds(simTrack.trackId(), simTrack.eventId()); + auto bestAssociatedRecoTrack = simTrackToRecoTrack(simTkIds); + if (bestAssociatedRecoTrack.first != -1 and bestAssociatedRecoTrack.second > qualityCutTrack_) { + auto trackIndex = bestAssociatedRecoTrack.first; + simTrackstersFromCP[i].setTrackIdx(trackIndex); + } + } + + auto& simTracksters = *result; + // Creating the map from TrackingParticle to SimTrackster + std::unordered_map> TPtoSimTracksterMap; + for (unsigned int i = 0; i < simTracksters.size(); ++i) { + const auto& simTrack = (simTracksters[i].seedID() == caloParticles_h.id()) + ? caloparticles[simTracksters[i].seedIndex()].g4Tracks()[0] + : simclusters[simTracksters[i].seedIndex()].g4Tracks()[0]; + UniqueSimTrackId simTkIds(simTrack.trackId(), simTrack.eventId()); + auto bestAssociatedRecoTrack = simTrackToRecoTrack(simTkIds); + if (bestAssociatedRecoTrack.first != -1 and bestAssociatedRecoTrack.second > qualityCutTrack_) { + auto trackIndex = bestAssociatedRecoTrack.first; + simTracksters[i].setTrackIdx(trackIndex); + } + } + + edm::OrphanHandle> simTracksters_h = evt.put(std::move(result)); + + result_ticlCandidates->resize(result_fromCP->size()); + std::vector toKeep; + for (size_t i = 0; i < simTracksters_h->size(); ++i) { + const auto& simTrackster = (*simTracksters_h)[i]; + int cp_index = (simTrackster.seedID() == caloParticles_h.id()) + ? simTrackster.seedIndex() + : SimClusterToCaloParticleMap[simTrackster.seedIndex()]; + auto const& tCP = (*result_fromCP)[cp_index]; + if (!tCP.vertices().empty()) { + auto trackIndex = tCP.trackIdx(); + + auto& cand = (*result_ticlCandidates)[cp_index]; + cand.addTrackster(edm::Ptr(simTracksters_h, i)); + if (trackIndex != -1 and (trackIndex < 0 or trackIndex >= (long int)recoTracks.size())) { + } + cand.setTime((*result_fromCP)[cp_index].time()); + cand.setTimeError(0); + if (trackIndex != -1 && caloparticles[cp_index].charge() != 0) + cand.setTrackPtr(edm::Ptr(recoTracks_h, trackIndex)); + toKeep.push_back(cp_index); + } + } + + auto isHad = [](int pdgId) { + pdgId = std::abs(pdgId); + if (pdgId == 111) + return false; + return (pdgId > 100 and pdgId < 900) or (pdgId > 1000 and pdgId < 9000); + }; + + for (size_t i = 0; i < result_ticlCandidates->size(); ++i) { + auto cp_index = (*result_fromCP)[i].seedIndex(); + if (cp_index < 0) + continue; + auto& cand = (*result_ticlCandidates)[i]; + const auto& cp = caloparticles[cp_index]; + float rawEnergy = 0.f; + float regressedEnergy = 0.f; + + for (const auto& trackster : cand.tracksters()) { + rawEnergy += trackster->raw_energy(); + regressedEnergy += trackster->regressed_energy(); + } + cand.setRawEnergy(rawEnergy); + + auto pdgId = cp.pdgId(); + auto charge = cp.charge(); + if (cand.trackPtr().isNonnull() and charge == 0) { + } + if (cand.trackPtr().isNonnull() and charge != 0) { + auto const& track = cand.trackPtr().get(); + if (std::abs(pdgId) == 13) { + cand.setPdgId(pdgId); + } else { + cand.setPdgId((isHad(pdgId) ? 211 : 11) * charge); + } + cand.setCharge(charge); + math::XYZTLorentzVector p4(regressedEnergy * track->momentum().unit().x(), + regressedEnergy * track->momentum().unit().y(), + regressedEnergy * track->momentum().unit().z(), + regressedEnergy); + cand.setP4(p4); + } else { // neutral candidates + cand.setPdgId(isHad(pdgId) ? 130 : 22); + cand.setCharge(0); + + auto particleType = tracksterParticleTypeFromPdgId(cand.pdgId(), 1); + cand.setIdProbability(particleType, 1.f); + + const auto& simTracksterFromCP = (*result_fromCP)[i]; + float regressedEnergy = simTracksterFromCP.regressed_energy(); + math::XYZTLorentzVector p4(regressedEnergy * simTracksterFromCP.barycenter().unit().x(), + regressedEnergy * simTracksterFromCP.barycenter().unit().y(), + regressedEnergy * simTracksterFromCP.barycenter().unit().z(), + regressedEnergy); + cand.setP4(p4); + } + } + + std::vector all_nums(result_fromCP->size()); // vector containing all caloparticles indexes + std::iota(all_nums.begin(), all_nums.end(), 0); // fill the vector with consecutive numbers starting from 0 + + std::vector toRemove; + std::set_difference(all_nums.begin(), all_nums.end(), toKeep.begin(), toKeep.end(), std::back_inserter(toRemove)); + std::sort(toRemove.begin(), toRemove.end(), [](int x, int y) { return x > y; }); + for (auto const& r : toRemove) { + result_fromCP->erase(result_fromCP->begin() + r); + result_ticlCandidates->erase(result_ticlCandidates->begin() + r); + } + evt.put(std::move(result_ticlCandidates)); evt.put(std::move(output_mask)); evt.put(std::move(result_fromCP), "fromCPs"); + evt.put(std::move(resultPU), "PU"); evt.put(std::move(output_mask_fromCP), "fromCPs"); evt.put(std::move(cpToSc_SimTrackstersMap)); } diff --git a/RecoHGCal/TICL/plugins/TICLDumper.cc b/RecoHGCal/TICL/plugins/TICLDumper.cc new file mode 100644 index 0000000000000..7897679e35e45 --- /dev/null +++ b/RecoHGCal/TICL/plugins/TICLDumper.cc @@ -0,0 +1,2102 @@ +// Original Authors: Philipp Zehetner, Wahid Redjeb + +#include "TTree.h" +#include "TFile.h" + +#include +#include +#include +#include + +#include // unique_ptr +#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/one/EDAnalyzer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ServiceRegistry/interface/Service.h" + +#include "DataFormats/CaloRecHit/interface/CaloCluster.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/DetId/interface/DetId.h" +#include "DataFormats/Math/interface/Point3D.h" +#include "DataFormats/GeometrySurface/interface/BoundDisk.h" +#include "DataFormats/HGCalReco/interface/Common.h" +#include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h" +#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h" + +#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.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 "Geometry/CaloGeometry/interface/CaloGeometry.h" +#include "Geometry/Records/interface/CaloGeometryRecord.h" +#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h" + +#include "SimDataFormats/Associations/interface/TracksterToSimTracksterHitLCAssociator.h" +#include "RecoHGCal/TICL/interface/commons.h" + +// TFileService +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "CommonTools/UtilAlgos/interface/TFileService.h" + +class TICLDumper : public edm::one::EDAnalyzer { +public: + explicit TICLDumper(const edm::ParameterSet&); + ~TICLDumper() override; + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + typedef math::XYZVector Vector; + typedef std::vector Vec; + +private: + void beginJob() override; + void beginRun(const edm::Run&, const edm::EventSetup&) override; + + void initialize(const HGCalDDDConstants* hgcons, + const hgcal::RecHitTools rhtools, + const edm::ESHandle bfieldH, + const edm::ESHandle propH); + void buildLayers(); + + void analyze(const edm::Event&, const edm::EventSetup&) override; + void endRun(edm::Run const& iEvent, edm::EventSetup const&) override{}; + void endJob() override; + + // Define Tokens + const edm::EDGetTokenT> tracksters_token_; + const edm::EDGetTokenT> layer_clusters_token_; + const edm::EDGetTokenT> ticl_candidates_token_; + const edm::EDGetTokenT> tracks_token_; + const edm::EDGetTokenT> tracks_mask_token_; + const edm::EDGetTokenT> tracks_time_token_; + const edm::EDGetTokenT> tracks_time_quality_token_; + const edm::EDGetTokenT> tracks_time_err_token_; + const edm::EDGetTokenT> hgcaltracks_x_token_; + const edm::EDGetTokenT> hgcaltracks_y_token_; + const edm::EDGetTokenT> hgcaltracks_z_token_; + const edm::EDGetTokenT> hgcaltracks_eta_token_; + const edm::EDGetTokenT> hgcaltracks_phi_token_; + const edm::EDGetTokenT> hgcaltracks_px_token_; + const edm::EDGetTokenT> hgcaltracks_py_token_; + const edm::EDGetTokenT> hgcaltracks_pz_token_; + const edm::EDGetTokenT> tracksters_merged_token_; + const edm::EDGetTokenT>> clustersTime_token_; + const edm::EDGetTokenT> tracksterSeeds_token_; + edm::ESGetToken caloGeometry_token_; + const edm::EDGetTokenT> simTracksters_SC_token_; + const edm::EDGetTokenT> simTracksters_CP_token_; + const edm::EDGetTokenT> simTracksters_PU_token_; + const edm::EDGetTokenT> simTICLCandidate_token_; + const edm::EDGetTokenT tsRecoToSimSC_token_; + const edm::EDGetTokenT tsSimToRecoSC_token_; + const edm::EDGetTokenT tsRecoToSimCP_token_; + const edm::EDGetTokenT tsSimToRecoCP_token_; + const edm::EDGetTokenT MergeRecoToSimSC_token_; + const edm::EDGetTokenT MergeSimToRecoSC_token_; + const edm::EDGetTokenT MergeRecoToSimCP_token_; + const edm::EDGetTokenT MergeSimToRecoCP_token_; + const edm::EDGetTokenT MergeRecoToSimPU_token_; + const edm::EDGetTokenT MergeSimToRecoPU_token_; + const edm::EDGetTokenT> simclusters_token_; + const edm::EDGetTokenT> caloparticles_token_; + + const edm::ESGetToken geometry_token_; + const std::string detector_; + const std::string propName_; + const edm::ESGetToken bfield_token_; + const edm::ESGetToken propagator_token_; + hgcal::RecHitTools rhtools_; + edm::ESGetToken hdc_token_; + const HGCalDDDConstants* hgcons_; + std::unique_ptr firstDisk_[2]; + std::unique_ptr interfaceDisk_[2]; + edm::ESHandle bfield_; + edm::ESHandle propagator_; + bool saveLCs_; + bool saveCLUE3DTracksters_; + bool saveTrackstersMerged_; + bool saveSimTrackstersSC_; + bool saveSimTrackstersCP_; + bool saveTICLCandidate_; + bool saveSimTICLCandidate_; + bool saveTracks_; + bool saveAssociations_; + + // Output tree + TTree* tree_; + + void clearVariables(); + + // Variables for branches + unsigned int ev_event_; + unsigned int ntracksters_; + unsigned int nclusters_; + unsigned int stsSC_ntracksters_; + unsigned int stsCP_ntracksters_; + size_t nsimTrackstersSC; + size_t nsimTrackstersCP; + + std::vector trackster_time; + std::vector trackster_timeError; + std::vector trackster_regressed_energy; + std::vector trackster_raw_energy; + std::vector trackster_raw_em_energy; + std::vector trackster_raw_pt; + std::vector trackster_raw_em_pt; + std::vector trackster_barycenter_x; + std::vector trackster_barycenter_y; + std::vector trackster_barycenter_z; + std::vector trackster_EV1; + std::vector trackster_EV2; + std::vector trackster_EV3; + std::vector trackster_eVector0_x; + std::vector trackster_eVector0_y; + std::vector trackster_eVector0_z; + std::vector trackster_sigmaPCA1; + std::vector trackster_sigmaPCA2; + std::vector trackster_sigmaPCA3; + std::vector trackster_barycenter_eta; + std::vector trackster_barycenter_phi; + std::vector> trackster_id_probabilities; + std::vector> trackster_vertices_indexes; + std::vector> trackster_vertices_x; + std::vector> trackster_vertices_y; + std::vector> trackster_vertices_z; + std::vector> trackster_vertices_time; + std::vector> trackster_vertices_timeErr; + std::vector> trackster_vertices_energy; + std::vector> trackster_vertices_correctedEnergy; + std::vector> trackster_vertices_correctedEnergyUncertainty; + std::vector> trackster_vertices_multiplicity; + + std::vector stsSC_trackster_time; + std::vector stsSC_trackster_timeError; + std::vector stsSC_trackster_regressed_energy; + std::vector stsSC_trackster_regressed_pt; + std::vector stsSC_trackster_raw_energy; + std::vector stsSC_trackster_raw_em_energy; + std::vector stsSC_trackster_raw_pt; + std::vector stsSC_trackster_raw_em_pt; + std::vector stsSC_trackster_barycenter_x; + std::vector stsSC_trackster_barycenter_y; + std::vector stsSC_trackster_barycenter_z; + std::vector stsSC_trackster_barycenter_eta; + std::vector stsSC_trackster_barycenter_phi; + std::vector stsSC_trackster_EV1; + std::vector stsSC_trackster_EV2; + std::vector stsSC_trackster_EV3; + std::vector stsSC_trackster_eVector0_x; + std::vector stsSC_trackster_eVector0_y; + std::vector stsSC_trackster_eVector0_z; + std::vector stsSC_trackster_sigmaPCA1; + std::vector stsSC_trackster_sigmaPCA2; + std::vector stsSC_trackster_sigmaPCA3; + std::vector stsSC_pdgID; + std::vector stsSC_trackIdx; + std::vector stsSC_trackTime; + std::vector stsSC_boundaryX; + std::vector stsSC_boundaryY; + std::vector stsSC_boundaryZ; + std::vector stsSC_boundaryEta; + std::vector stsSC_boundaryPhi; + std::vector stsSC_boundaryPx; + std::vector stsSC_boundaryPy; + std::vector stsSC_boundaryPz; + std::vector stsSC_track_boundaryX; + std::vector stsSC_track_boundaryY; + std::vector stsSC_track_boundaryZ; + std::vector stsSC_track_boundaryEta; + std::vector stsSC_track_boundaryPhi; + std::vector stsSC_track_boundaryPx; + std::vector stsSC_track_boundaryPy; + std::vector stsSC_track_boundaryPz; + std::vector> stsSC_trackster_id_probabilities; + std::vector> stsSC_trackster_vertices_indexes; + std::vector> stsSC_trackster_vertices_x; + std::vector> stsSC_trackster_vertices_y; + std::vector> stsSC_trackster_vertices_z; + std::vector> stsSC_trackster_vertices_time; + std::vector> stsSC_trackster_vertices_timeErr; + std::vector> stsSC_trackster_vertices_energy; + std::vector> stsSC_trackster_vertices_correctedEnergy; + std::vector> stsSC_trackster_vertices_correctedEnergyUncertainty; + std::vector> stsSC_trackster_vertices_multiplicity; + std::vector stsCP_trackster_time; + std::vector stsCP_trackster_timeError; + std::vector stsCP_trackster_regressed_energy; + std::vector stsCP_trackster_regressed_pt; + std::vector stsCP_trackster_raw_energy; + std::vector stsCP_trackster_raw_em_energy; + std::vector stsCP_trackster_raw_pt; + std::vector stsCP_trackster_raw_em_pt; + std::vector stsCP_trackster_barycenter_x; + std::vector stsCP_trackster_barycenter_y; + std::vector stsCP_trackster_barycenter_z; + std::vector stsCP_trackster_barycenter_eta; + std::vector stsCP_trackster_barycenter_phi; + std::vector stsCP_trackster_EV1; + std::vector stsCP_trackster_EV2; + std::vector stsCP_trackster_EV3; + std::vector stsCP_trackster_eVector0_x; + std::vector stsCP_trackster_eVector0_y; + std::vector stsCP_trackster_eVector0_z; + std::vector stsCP_trackster_sigmaPCA1; + std::vector stsCP_trackster_sigmaPCA2; + std::vector stsCP_trackster_sigmaPCA3; + std::vector stsCP_pdgID; + std::vector stsCP_trackIdx; + std::vector stsCP_trackTime; + std::vector stsCP_boundaryX; + std::vector stsCP_boundaryY; + std::vector stsCP_boundaryZ; + std::vector stsCP_boundaryEta; + std::vector stsCP_boundaryPhi; + std::vector stsCP_boundaryPx; + std::vector stsCP_boundaryPy; + std::vector stsCP_boundaryPz; + std::vector stsCP_track_boundaryX; + std::vector stsCP_track_boundaryY; + std::vector stsCP_track_boundaryZ; + std::vector stsCP_track_boundaryEta; + std::vector stsCP_track_boundaryPhi; + std::vector stsCP_track_boundaryPx; + std::vector stsCP_track_boundaryPy; + std::vector stsCP_track_boundaryPz; + std::vector> stsCP_trackster_id_probabilities; + std::vector> stsCP_trackster_vertices_indexes; + std::vector> stsCP_trackster_vertices_x; + std::vector> stsCP_trackster_vertices_y; + std::vector> stsCP_trackster_vertices_z; + std::vector> stsCP_trackster_vertices_time; + std::vector> stsCP_trackster_vertices_timeErr; + std::vector> stsCP_trackster_vertices_energy; + std::vector> stsCP_trackster_vertices_correctedEnergy; + std::vector> stsCP_trackster_vertices_correctedEnergyUncertainty; + std::vector> stsCP_trackster_vertices_multiplicity; + + std::vector simTICLCandidate_raw_energy; + std::vector simTICLCandidate_regressed_energy; + std::vector> simTICLCandidate_simTracksterCPIndex; + std::vector simTICLCandidate_boundaryX; + std::vector simTICLCandidate_boundaryY; + std::vector simTICLCandidate_boundaryZ; + std::vector simTICLCandidate_boundaryPx; + std::vector simTICLCandidate_boundaryPy; + std::vector simTICLCandidate_boundaryPz; + std::vector simTICLCandidate_trackTime; + std::vector simTICLCandidate_trackBeta; + std::vector simTICLCandidate_caloParticleMass; + std::vector simTICLCandidate_pdgId; + std::vector simTICLCandidate_charge; + std::vector simTICLCandidate_track_in_candidate; + + // from TICLCandidate, product of linking + size_t nCandidates; + std::vector candidate_charge; + std::vector candidate_pdgId; + std::vector candidate_energy; + std::vector candidate_px; + std::vector candidate_py; + std::vector candidate_pz; + std::vector candidate_time; + std::vector candidate_time_err; + std::vector> candidate_id_probabilities; + std::vector> tracksters_in_candidate; + std::vector track_in_candidate; + + // merged tracksters + size_t nTrackstersMerged; + std::vector tracksters_merged_time; + std::vector tracksters_merged_timeError; + std::vector tracksters_merged_regressed_energy; + std::vector tracksters_merged_raw_energy; + std::vector tracksters_merged_raw_em_energy; + std::vector tracksters_merged_raw_pt; + std::vector tracksters_merged_raw_em_pt; + std::vector tracksters_merged_barycenter_x; + std::vector tracksters_merged_barycenter_y; + std::vector tracksters_merged_barycenter_z; + std::vector tracksters_merged_barycenter_eta; + std::vector tracksters_merged_barycenter_phi; + std::vector tracksters_merged_EV1; + std::vector tracksters_merged_EV2; + std::vector tracksters_merged_EV3; + std::vector tracksters_merged_eVector0_x; + std::vector tracksters_merged_eVector0_y; + std::vector tracksters_merged_eVector0_z; + std::vector tracksters_merged_sigmaPCA1; + std::vector tracksters_merged_sigmaPCA2; + std::vector tracksters_merged_sigmaPCA3; + std::vector> tracksters_merged_vertices_indexes; + std::vector> tracksters_merged_vertices_x; + std::vector> tracksters_merged_vertices_y; + std::vector> tracksters_merged_vertices_z; + std::vector> tracksters_merged_vertices_time; + std::vector> tracksters_merged_vertices_timeErr; + std::vector> tracksters_merged_vertices_energy; + std::vector> tracksters_merged_vertices_correctedEnergy; + std::vector> tracksters_merged_vertices_correctedEnergyUncertainty; + std::vector> tracksters_merged_vertices_multiplicity; + std::vector> tracksters_merged_id_probabilities; + + // associations + std::vector> trackstersCLUE3D_recoToSim_SC; + std::vector> trackstersCLUE3D_recoToSim_SC_score; + std::vector> trackstersCLUE3D_recoToSim_SC_sharedE; + std::vector> trackstersCLUE3D_simToReco_SC; + std::vector> trackstersCLUE3D_simToReco_SC_score; + std::vector> trackstersCLUE3D_simToReco_SC_sharedE; + + std::vector> trackstersCLUE3D_recoToSim_CP; + std::vector> trackstersCLUE3D_recoToSim_CP_score; + std::vector> trackstersCLUE3D_recoToSim_CP_sharedE; + std::vector> trackstersCLUE3D_simToReco_CP; + std::vector> trackstersCLUE3D_simToReco_CP_score; + std::vector> trackstersCLUE3D_simToReco_CP_sharedE; + + std::vector> MergeTracksters_recoToSim_SC; + std::vector> MergeTracksters_recoToSim_SC_score; + std::vector> MergeTracksters_recoToSim_SC_sharedE; + std::vector> MergeTracksters_simToReco_SC; + std::vector> MergeTracksters_simToReco_SC_score; + std::vector> MergeTracksters_simToReco_SC_sharedE; + + std::vector> MergeTracksters_recoToSim_CP; + std::vector> MergeTracksters_recoToSim_CP_score; + std::vector> MergeTracksters_recoToSim_CP_sharedE; + std::vector> MergeTracksters_simToReco_CP; + std::vector> MergeTracksters_simToReco_CP_score; + std::vector> MergeTracksters_simToReco_CP_sharedE; + + std::vector> MergeTracksters_recoToSim_PU; + std::vector> MergeTracksters_recoToSim_PU_score; + std::vector> MergeTracksters_recoToSim_PU_sharedE; + std::vector> MergeTracksters_simToReco_PU; + std::vector> MergeTracksters_simToReco_PU_score; + std::vector> MergeTracksters_simToReco_PU_sharedE; + + std::vector cluster_seedID; + std::vector cluster_energy; + std::vector cluster_correctedEnergy; + std::vector cluster_correctedEnergyUncertainty; + std::vector cluster_position_x; + std::vector cluster_position_y; + std::vector cluster_position_z; + std::vector cluster_position_eta; + std::vector cluster_position_phi; + std::vector cluster_layer_id; + std::vector cluster_type; + std::vector cluster_time; + std::vector cluster_timeErr; + std::vector cluster_number_of_hits; + + std::vector track_id; + std::vector track_hgcal_x; + std::vector track_hgcal_y; + std::vector track_hgcal_z; + std::vector track_hgcal_px; + std::vector track_hgcal_py; + std::vector track_hgcal_pz; + std::vector track_hgcal_eta; + std::vector track_hgcal_phi; + std::vector track_hgcal_pt; + std::vector track_pt; + std::vector track_quality; + std::vector track_missing_outer_hits; + std::vector track_charge; + std::vector track_time; + std::vector track_time_quality; + std::vector track_time_err; + std::vector track_nhits; + + TTree* trackster_tree_; + TTree* cluster_tree_; + TTree* candidate_tree_; + TTree* tracksters_merged_tree_; + TTree* associations_tree_; + TTree* simtrackstersSC_tree_; + TTree* simtrackstersCP_tree_; + TTree* tracks_tree_; + TTree* simTICLCandidate_tree; +}; + +void TICLDumper::clearVariables() { + // event info + ntracksters_ = 0; + nclusters_ = 0; + + trackster_time.clear(); + trackster_timeError.clear(); + trackster_regressed_energy.clear(); + trackster_raw_energy.clear(); + trackster_raw_em_energy.clear(); + trackster_raw_pt.clear(); + trackster_raw_em_pt.clear(); + trackster_barycenter_x.clear(); + trackster_barycenter_y.clear(); + trackster_barycenter_z.clear(); + trackster_EV1.clear(); + trackster_EV2.clear(); + trackster_EV3.clear(); + trackster_eVector0_x.clear(); + trackster_eVector0_y.clear(); + trackster_eVector0_z.clear(); + trackster_sigmaPCA1.clear(); + trackster_sigmaPCA2.clear(); + trackster_sigmaPCA3.clear(); + trackster_barycenter_eta.clear(); + trackster_barycenter_phi.clear(); + trackster_id_probabilities.clear(); + trackster_vertices_indexes.clear(); + trackster_vertices_x.clear(); + trackster_vertices_y.clear(); + trackster_vertices_z.clear(); + trackster_vertices_time.clear(); + trackster_vertices_timeErr.clear(); + trackster_vertices_energy.clear(); + trackster_vertices_correctedEnergy.clear(); + trackster_vertices_correctedEnergyUncertainty.clear(); + trackster_vertices_multiplicity.clear(); + + stsSC_trackster_time.clear(); + stsSC_trackster_timeError.clear(); + stsSC_trackster_regressed_energy.clear(); + stsSC_trackster_regressed_pt.clear(); + stsSC_trackster_raw_energy.clear(); + stsSC_trackster_raw_em_energy.clear(); + stsSC_trackster_raw_pt.clear(); + stsSC_trackster_raw_em_pt.clear(); + stsSC_trackster_barycenter_x.clear(); + stsSC_trackster_barycenter_y.clear(); + stsSC_trackster_barycenter_z.clear(); + stsSC_trackster_EV1.clear(); + stsSC_trackster_EV2.clear(); + stsSC_trackster_EV3.clear(); + stsSC_trackster_eVector0_x.clear(); + stsSC_trackster_eVector0_y.clear(); + stsSC_trackster_eVector0_z.clear(); + stsSC_trackster_sigmaPCA1.clear(); + stsSC_trackster_sigmaPCA2.clear(); + stsSC_trackster_sigmaPCA3.clear(); + stsSC_trackster_barycenter_eta.clear(); + stsSC_trackster_barycenter_phi.clear(); + stsSC_pdgID.clear(); + stsSC_trackIdx.clear(); + stsSC_trackTime.clear(); + stsSC_boundaryX.clear(); + stsSC_boundaryY.clear(); + stsSC_boundaryZ.clear(); + stsSC_boundaryEta.clear(); + stsSC_boundaryPhi.clear(); + stsSC_boundaryPx.clear(); + stsSC_boundaryPy.clear(); + stsSC_boundaryPz.clear(); + stsSC_track_boundaryX.clear(); + stsSC_track_boundaryY.clear(); + stsSC_track_boundaryZ.clear(); + stsSC_track_boundaryEta.clear(); + stsSC_track_boundaryPhi.clear(); + stsSC_track_boundaryPx.clear(); + stsSC_track_boundaryPy.clear(); + stsSC_track_boundaryPz.clear(); + stsSC_trackster_id_probabilities.clear(); + stsSC_trackster_vertices_indexes.clear(); + stsSC_trackster_vertices_x.clear(); + stsSC_trackster_vertices_y.clear(); + stsSC_trackster_vertices_z.clear(); + stsSC_trackster_vertices_time.clear(); + stsSC_trackster_vertices_timeErr.clear(); + stsSC_trackster_vertices_energy.clear(); + stsSC_trackster_vertices_correctedEnergy.clear(); + stsSC_trackster_vertices_correctedEnergyUncertainty.clear(); + stsSC_trackster_vertices_multiplicity.clear(); + + stsCP_trackster_time.clear(); + stsCP_trackster_timeError.clear(); + stsCP_trackster_regressed_energy.clear(); + stsCP_trackster_regressed_pt.clear(); + stsCP_trackster_raw_energy.clear(); + stsCP_trackster_raw_em_energy.clear(); + stsCP_trackster_raw_pt.clear(); + stsCP_trackster_raw_em_pt.clear(); + stsCP_trackster_barycenter_x.clear(); + stsCP_trackster_barycenter_y.clear(); + stsCP_trackster_barycenter_z.clear(); + stsCP_trackster_sigmaPCA1.clear(); + stsCP_trackster_sigmaPCA2.clear(); + stsCP_trackster_sigmaPCA3.clear(); + stsCP_trackster_barycenter_eta.clear(); + stsCP_trackster_barycenter_phi.clear(); + stsCP_pdgID.clear(); + stsCP_trackIdx.clear(); + stsCP_trackTime.clear(); + stsCP_boundaryX.clear(); + stsCP_boundaryY.clear(); + stsCP_boundaryZ.clear(); + stsCP_boundaryEta.clear(); + stsCP_boundaryPhi.clear(); + stsCP_boundaryPx.clear(); + stsCP_boundaryPy.clear(); + stsCP_boundaryPz.clear(); + stsCP_track_boundaryX.clear(); + stsCP_track_boundaryY.clear(); + stsCP_track_boundaryZ.clear(); + stsCP_track_boundaryEta.clear(); + stsCP_track_boundaryPhi.clear(); + stsCP_track_boundaryPx.clear(); + stsCP_track_boundaryPy.clear(); + stsCP_track_boundaryPz.clear(); + stsCP_trackster_id_probabilities.clear(); + stsCP_trackster_vertices_indexes.clear(); + stsCP_trackster_vertices_x.clear(); + stsCP_trackster_vertices_y.clear(); + stsCP_trackster_vertices_z.clear(); + stsCP_trackster_vertices_time.clear(); + stsCP_trackster_vertices_timeErr.clear(); + stsCP_trackster_vertices_energy.clear(); + stsCP_trackster_vertices_correctedEnergy.clear(); + stsCP_trackster_vertices_correctedEnergyUncertainty.clear(); + stsCP_trackster_vertices_multiplicity.clear(); + + simTICLCandidate_raw_energy.clear(); + simTICLCandidate_regressed_energy.clear(); + simTICLCandidate_simTracksterCPIndex.clear(); + simTICLCandidate_boundaryX.clear(); + simTICLCandidate_boundaryY.clear(); + simTICLCandidate_boundaryZ.clear(); + simTICLCandidate_boundaryPx.clear(); + simTICLCandidate_boundaryPy.clear(); + simTICLCandidate_boundaryPz.clear(); + simTICLCandidate_trackTime.clear(); + simTICLCandidate_trackBeta.clear(); + simTICLCandidate_caloParticleMass.clear(); + simTICLCandidate_pdgId.clear(); + simTICLCandidate_charge.clear(); + simTICLCandidate_track_in_candidate.clear(); + + nCandidates = 0; + candidate_charge.clear(); + candidate_pdgId.clear(); + candidate_energy.clear(); + candidate_px.clear(); + candidate_py.clear(); + candidate_pz.clear(); + candidate_time.clear(); + candidate_time_err.clear(); + candidate_id_probabilities.clear(); + tracksters_in_candidate.clear(); + track_in_candidate.clear(); + + nTrackstersMerged = 0; + tracksters_merged_time.clear(); + tracksters_merged_timeError.clear(); + tracksters_merged_regressed_energy.clear(); + tracksters_merged_raw_energy.clear(); + tracksters_merged_raw_em_energy.clear(); + tracksters_merged_raw_pt.clear(); + tracksters_merged_raw_em_pt.clear(); + tracksters_merged_barycenter_x.clear(); + tracksters_merged_barycenter_y.clear(); + tracksters_merged_barycenter_z.clear(); + tracksters_merged_barycenter_eta.clear(); + tracksters_merged_barycenter_phi.clear(); + tracksters_merged_EV1.clear(); + tracksters_merged_EV2.clear(); + tracksters_merged_EV3.clear(); + tracksters_merged_eVector0_x.clear(); + tracksters_merged_eVector0_y.clear(); + tracksters_merged_eVector0_z.clear(); + tracksters_merged_sigmaPCA1.clear(); + tracksters_merged_sigmaPCA2.clear(); + tracksters_merged_sigmaPCA3.clear(); + tracksters_merged_id_probabilities.clear(); + tracksters_merged_time.clear(); + tracksters_merged_timeError.clear(); + tracksters_merged_regressed_energy.clear(); + tracksters_merged_raw_energy.clear(); + tracksters_merged_raw_em_energy.clear(); + tracksters_merged_raw_pt.clear(); + tracksters_merged_raw_em_pt.clear(); + + tracksters_merged_vertices_indexes.clear(); + tracksters_merged_vertices_x.clear(); + tracksters_merged_vertices_y.clear(); + tracksters_merged_vertices_z.clear(); + tracksters_merged_vertices_time.clear(); + tracksters_merged_vertices_timeErr.clear(); + tracksters_merged_vertices_energy.clear(); + tracksters_merged_vertices_correctedEnergy.clear(); + tracksters_merged_vertices_correctedEnergyUncertainty.clear(); + tracksters_merged_vertices_multiplicity.clear(); + + trackstersCLUE3D_recoToSim_SC.clear(); + trackstersCLUE3D_recoToSim_SC_score.clear(); + trackstersCLUE3D_recoToSim_SC_sharedE.clear(); + trackstersCLUE3D_simToReco_SC.clear(); + trackstersCLUE3D_simToReco_SC_score.clear(); + trackstersCLUE3D_simToReco_SC_sharedE.clear(); + + trackstersCLUE3D_recoToSim_CP.clear(); + trackstersCLUE3D_recoToSim_CP_score.clear(); + trackstersCLUE3D_recoToSim_CP_sharedE.clear(); + trackstersCLUE3D_simToReco_CP.clear(); + trackstersCLUE3D_simToReco_CP_score.clear(); + trackstersCLUE3D_simToReco_CP_sharedE.clear(); + + MergeTracksters_recoToSim_SC.clear(); + MergeTracksters_recoToSim_SC_score.clear(); + MergeTracksters_recoToSim_SC_sharedE.clear(); + MergeTracksters_simToReco_SC.clear(); + MergeTracksters_simToReco_SC_score.clear(); + MergeTracksters_simToReco_SC_sharedE.clear(); + + MergeTracksters_recoToSim_CP.clear(); + MergeTracksters_recoToSim_CP_score.clear(); + MergeTracksters_recoToSim_CP_sharedE.clear(); + MergeTracksters_simToReco_CP.clear(); + MergeTracksters_simToReco_CP_score.clear(); + MergeTracksters_simToReco_CP_sharedE.clear(); + + MergeTracksters_recoToSim_PU.clear(); + MergeTracksters_recoToSim_PU_score.clear(); + MergeTracksters_recoToSim_PU_sharedE.clear(); + MergeTracksters_simToReco_PU.clear(); + MergeTracksters_simToReco_PU_score.clear(); + MergeTracksters_simToReco_PU_sharedE.clear(); + + nsimTrackstersSC = 0; + + cluster_seedID.clear(); + cluster_energy.clear(); + cluster_correctedEnergy.clear(); + cluster_correctedEnergyUncertainty.clear(); + cluster_position_x.clear(); + cluster_position_y.clear(); + cluster_position_z.clear(); + cluster_position_eta.clear(); + cluster_position_phi.clear(); + cluster_layer_id.clear(); + cluster_type.clear(); + cluster_time.clear(); + cluster_timeErr.clear(); + cluster_number_of_hits.clear(); + + track_id.clear(); + track_hgcal_x.clear(); + track_hgcal_y.clear(); + track_hgcal_z.clear(); + track_hgcal_eta.clear(); + track_hgcal_phi.clear(); + track_hgcal_px.clear(); + track_hgcal_py.clear(); + track_hgcal_pz.clear(); + track_hgcal_pt.clear(); + track_quality.clear(); + track_pt.clear(); + track_missing_outer_hits.clear(); + track_charge.clear(); + track_time.clear(); + track_time_quality.clear(); + track_time_err.clear(); + track_nhits.clear(); +}; + +TICLDumper::TICLDumper(const edm::ParameterSet& ps) + : tracksters_token_(consumes>(ps.getParameter("trackstersclue3d"))), + layer_clusters_token_(consumes>(ps.getParameter("layerClusters"))), + ticl_candidates_token_(consumes>(ps.getParameter("ticlcandidates"))), + tracks_token_(consumes>(ps.getParameter("tracks"))), + tracks_time_token_(consumes>(ps.getParameter("tracksTime"))), + tracks_time_quality_token_(consumes>(ps.getParameter("tracksTimeQual"))), + tracks_time_err_token_(consumes>(ps.getParameter("tracksTimeErr"))), + tracksters_merged_token_( + consumes>(ps.getParameter("trackstersmerged"))), + clustersTime_token_( + consumes>>(ps.getParameter("layer_clustersTime"))), + caloGeometry_token_(esConsumes()), + simTracksters_SC_token_( + consumes>(ps.getParameter("simtrackstersSC"))), + simTracksters_CP_token_( + consumes>(ps.getParameter("simtrackstersCP"))), + simTracksters_PU_token_( + consumes>(ps.getParameter("simtrackstersPU"))), + simTICLCandidate_token_( + consumes>(ps.getParameter("simTICLCandidates"))), + tsRecoToSimSC_token_( + consumes(ps.getParameter("recoToSimAssociatorSC"))), + tsSimToRecoSC_token_( + consumes(ps.getParameter("simToRecoAssociatorSC"))), + tsRecoToSimCP_token_( + consumes(ps.getParameter("recoToSimAssociatorCP"))), + tsSimToRecoCP_token_( + consumes(ps.getParameter("simToRecoAssociatorCP"))), + MergeRecoToSimSC_token_(consumes( + ps.getParameter("MergerecoToSimAssociatorSC"))), + MergeSimToRecoSC_token_(consumes( + ps.getParameter("MergesimToRecoAssociatorSC"))), + MergeRecoToSimCP_token_(consumes( + ps.getParameter("MergerecoToSimAssociatorCP"))), + MergeSimToRecoCP_token_(consumes( + ps.getParameter("MergesimToRecoAssociatorCP"))), + MergeRecoToSimPU_token_(consumes( + ps.getParameter("MergerecoToSimAssociatorPU"))), + MergeSimToRecoPU_token_(consumes( + ps.getParameter("MergesimToRecoAssociatorPU"))), + simclusters_token_(consumes(ps.getParameter("simclusters"))), + caloparticles_token_(consumes(ps.getParameter("caloparticles"))), + geometry_token_(esConsumes()), + detector_(ps.getParameter("detector")), + propName_(ps.getParameter("propagator")), + bfield_token_(esConsumes()), + propagator_token_( + esConsumes(edm::ESInputTag("", propName_))), + saveLCs_(ps.getParameter("saveLCs")), + saveCLUE3DTracksters_(ps.getParameter("saveCLUE3DTracksters")), + saveTrackstersMerged_(ps.getParameter("saveTrackstersMerged")), + saveSimTrackstersSC_(ps.getParameter("saveSimTrackstersSC")), + saveSimTrackstersCP_(ps.getParameter("saveSimTrackstersCP")), + saveTICLCandidate_(ps.getParameter("saveSimTICLCandidate")), + saveSimTICLCandidate_(ps.getParameter("saveSimTICLCandidate")), + saveTracks_(ps.getParameter("saveTracks")), + saveAssociations_(ps.getParameter("saveAssociations")) { + std::string detectorName_ = (detector_ == "HFNose") ? "HGCalHFNoseSensitive" : "HGCalEESensitive"; + hdc_token_ = + esConsumes(edm::ESInputTag("", detectorName_)); +}; + +TICLDumper::~TICLDumper() { clearVariables(); }; + +void TICLDumper::beginRun(edm::Run const&, edm::EventSetup const& es) { + const CaloGeometry& geom = es.getData(caloGeometry_token_); + rhtools_.setGeometry(geom); + + edm::ESHandle hdc = es.getHandle(hdc_token_); + hgcons_ = hdc.product(); + edm::ESHandle bfield_ = es.getHandle(bfield_token_); + edm::ESHandle propagator = es.getHandle(propagator_token_); + initialize(hgcons_, rhtools_, bfield_, propagator); +} + +// Define tree and branches +void TICLDumper::beginJob() { + edm::Service fs; + if (saveCLUE3DTracksters_) { + trackster_tree_ = fs->make("tracksters", "TICL tracksters"); + trackster_tree_->Branch("event", &ev_event_); + trackster_tree_->Branch("NClusters", &nclusters_); + trackster_tree_->Branch("NTracksters", &ntracksters_); + trackster_tree_->Branch("time", &trackster_time); + trackster_tree_->Branch("timeError", &trackster_timeError); + trackster_tree_->Branch("regressed_energy", &trackster_regressed_energy); + trackster_tree_->Branch("raw_energy", &trackster_raw_energy); + trackster_tree_->Branch("raw_em_energy", &trackster_raw_em_energy); + trackster_tree_->Branch("raw_pt", &trackster_raw_pt); + trackster_tree_->Branch("raw_em_pt", &trackster_raw_em_pt); + trackster_tree_->Branch("barycenter_x", &trackster_barycenter_x); + trackster_tree_->Branch("barycenter_y", &trackster_barycenter_y); + trackster_tree_->Branch("barycenter_z", &trackster_barycenter_z); + trackster_tree_->Branch("barycenter_eta", &trackster_barycenter_eta); + trackster_tree_->Branch("barycenter_phi", &trackster_barycenter_phi); + trackster_tree_->Branch("EV1", &trackster_EV1); + trackster_tree_->Branch("EV2", &trackster_EV2); + trackster_tree_->Branch("EV3", &trackster_EV3); + trackster_tree_->Branch("eVector0_x", &trackster_eVector0_x); + trackster_tree_->Branch("eVector0_y", &trackster_eVector0_y); + trackster_tree_->Branch("eVector0_z", &trackster_eVector0_z); + trackster_tree_->Branch("sigmaPCA1", &trackster_sigmaPCA1); + trackster_tree_->Branch("sigmaPCA2", &trackster_sigmaPCA2); + trackster_tree_->Branch("sigmaPCA3", &trackster_sigmaPCA3); + trackster_tree_->Branch("id_probabilities", &trackster_id_probabilities); + trackster_tree_->Branch("vertices_indexes", &trackster_vertices_indexes); + trackster_tree_->Branch("vertices_x", &trackster_vertices_x); + trackster_tree_->Branch("vertices_y", &trackster_vertices_y); + trackster_tree_->Branch("vertices_z", &trackster_vertices_z); + trackster_tree_->Branch("vertices_time", &trackster_vertices_time); + trackster_tree_->Branch("vertices_timeErr", &trackster_vertices_timeErr); + trackster_tree_->Branch("vertices_energy", &trackster_vertices_energy); + trackster_tree_->Branch("vertices_correctedEnergy", &trackster_vertices_correctedEnergy); + trackster_tree_->Branch("vertices_correctedEnergyUncertainty", &trackster_vertices_correctedEnergyUncertainty); + trackster_tree_->Branch("vertices_multiplicity", &trackster_vertices_multiplicity); + } + if (saveLCs_) { + cluster_tree_ = fs->make("clusters", "TICL tracksters"); + cluster_tree_->Branch("seedID", &cluster_seedID); + cluster_tree_->Branch("energy", &cluster_energy); + cluster_tree_->Branch("correctedEnergy", &cluster_correctedEnergy); + cluster_tree_->Branch("correctedEnergyUncertainty", &cluster_correctedEnergyUncertainty); + cluster_tree_->Branch("position_x", &cluster_position_x); + cluster_tree_->Branch("position_y", &cluster_position_y); + cluster_tree_->Branch("position_z", &cluster_position_z); + cluster_tree_->Branch("position_eta", &cluster_position_eta); + cluster_tree_->Branch("position_phi", &cluster_position_phi); + cluster_tree_->Branch("cluster_layer_id", &cluster_layer_id); + cluster_tree_->Branch("cluster_type", &cluster_type); + cluster_tree_->Branch("cluster_time", &cluster_time); + cluster_tree_->Branch("cluster_timeErr", &cluster_timeErr); + cluster_tree_->Branch("cluster_number_of_hits", &cluster_number_of_hits); + } + if (saveTICLCandidate_) { + candidate_tree_ = fs->make("candidates", "TICL candidates"); + candidate_tree_->Branch("NCandidates", &nCandidates); + candidate_tree_->Branch("candidate_charge", &candidate_charge); + candidate_tree_->Branch("candidate_pdgId", &candidate_pdgId); + candidate_tree_->Branch("candidate_id_probabilities", &candidate_id_probabilities); + candidate_tree_->Branch("candidate_time", &candidate_time); + candidate_tree_->Branch("candidate_timeErr", &candidate_time_err); + candidate_tree_->Branch("candidate_energy", &candidate_energy); + candidate_tree_->Branch("candidate_px", &candidate_px); + candidate_tree_->Branch("candidate_py", &candidate_py); + candidate_tree_->Branch("candidate_pz", &candidate_pz); + candidate_tree_->Branch("track_in_candidate", &track_in_candidate); + candidate_tree_->Branch("tracksters_in_candidate", &tracksters_in_candidate); + } + if (saveTrackstersMerged_) { + tracksters_merged_tree_ = fs->make("trackstersMerged", "TICL tracksters merged"); + tracksters_merged_tree_->Branch("event", &ev_event_); + tracksters_merged_tree_->Branch("time", &tracksters_merged_time); + tracksters_merged_tree_->Branch("timeError", &tracksters_merged_timeError); + tracksters_merged_tree_->Branch("regressed_energy", &tracksters_merged_regressed_energy); + tracksters_merged_tree_->Branch("raw_energy", &tracksters_merged_raw_energy); + tracksters_merged_tree_->Branch("raw_em_energy", &tracksters_merged_raw_em_energy); + tracksters_merged_tree_->Branch("raw_pt", &tracksters_merged_raw_pt); + tracksters_merged_tree_->Branch("raw_em_pt", &tracksters_merged_raw_em_pt); + tracksters_merged_tree_->Branch("NTrackstersMerged", &nTrackstersMerged); + tracksters_merged_tree_->Branch("barycenter_x", &tracksters_merged_barycenter_x); + tracksters_merged_tree_->Branch("barycenter_y", &tracksters_merged_barycenter_y); + tracksters_merged_tree_->Branch("barycenter_z", &tracksters_merged_barycenter_z); + tracksters_merged_tree_->Branch("barycenter_eta", &tracksters_merged_barycenter_eta); + tracksters_merged_tree_->Branch("barycenter_phi", &tracksters_merged_barycenter_phi); + tracksters_merged_tree_->Branch("EV1", &tracksters_merged_EV1); + tracksters_merged_tree_->Branch("EV2", &tracksters_merged_EV2); + tracksters_merged_tree_->Branch("EV3", &tracksters_merged_EV3); + tracksters_merged_tree_->Branch("eVector0_x", &tracksters_merged_eVector0_x); + tracksters_merged_tree_->Branch("eVector0_y", &tracksters_merged_eVector0_y); + tracksters_merged_tree_->Branch("eVector0_z", &tracksters_merged_eVector0_z); + tracksters_merged_tree_->Branch("sigmaPCA1", &tracksters_merged_sigmaPCA1); + tracksters_merged_tree_->Branch("sigmaPCA2", &tracksters_merged_sigmaPCA2); + tracksters_merged_tree_->Branch("sigmaPCA3", &tracksters_merged_sigmaPCA3); + tracksters_merged_tree_->Branch("id_probabilities", &tracksters_merged_id_probabilities); + tracksters_merged_tree_->Branch("vertices_indexes", &tracksters_merged_vertices_indexes); + tracksters_merged_tree_->Branch("vertices_x", &tracksters_merged_vertices_x); + tracksters_merged_tree_->Branch("vertices_y", &tracksters_merged_vertices_y); + tracksters_merged_tree_->Branch("vertices_z", &tracksters_merged_vertices_z); + tracksters_merged_tree_->Branch("vertices_time", &tracksters_merged_vertices_time); + tracksters_merged_tree_->Branch("vertices_timeErr", &tracksters_merged_vertices_timeErr); + tracksters_merged_tree_->Branch("vertices_energy", &tracksters_merged_vertices_energy); + tracksters_merged_tree_->Branch("vertices_correctedEnergy", &tracksters_merged_vertices_correctedEnergy); + tracksters_merged_tree_->Branch("vertices_correctedEnergyUncertainty", + &tracksters_merged_vertices_correctedEnergyUncertainty); + tracksters_merged_tree_->Branch("vertices_multiplicity", &tracksters_merged_vertices_multiplicity); + } + if (saveAssociations_) { + associations_tree_ = fs->make("associations", "Associations"); + associations_tree_->Branch("tsCLUE3D_recoToSim_SC", &trackstersCLUE3D_recoToSim_SC); + associations_tree_->Branch("tsCLUE3D_recoToSim_SC_score", &trackstersCLUE3D_recoToSim_SC_score); + associations_tree_->Branch("tsCLUE3D_recoToSim_SC_sharedE", &trackstersCLUE3D_recoToSim_SC_sharedE); + associations_tree_->Branch("tsCLUE3D_simToReco_SC", &trackstersCLUE3D_simToReco_SC); + associations_tree_->Branch("tsCLUE3D_simToReco_SC_score", &trackstersCLUE3D_simToReco_SC_score); + associations_tree_->Branch("tsCLUE3D_simToReco_SC_sharedE", &trackstersCLUE3D_simToReco_SC_sharedE); + + associations_tree_->Branch("tsCLUE3D_recoToSim_CP", &trackstersCLUE3D_recoToSim_CP); + associations_tree_->Branch("tsCLUE3D_recoToSim_CP_score", &trackstersCLUE3D_recoToSim_CP_score); + associations_tree_->Branch("tsCLUE3D_recoToSim_CP_sharedE", &trackstersCLUE3D_recoToSim_CP_sharedE); + associations_tree_->Branch("tsCLUE3D_simToReco_CP", &trackstersCLUE3D_simToReco_CP); + associations_tree_->Branch("tsCLUE3D_simToReco_CP_score", &trackstersCLUE3D_simToReco_CP_score); + associations_tree_->Branch("tsCLUE3D_simToReco_CP_sharedE", &trackstersCLUE3D_simToReco_CP_sharedE); + + associations_tree_->Branch("Mergetstracksters_recoToSim_SC", &MergeTracksters_recoToSim_SC); + associations_tree_->Branch("Mergetstracksters_recoToSim_SC_score", &MergeTracksters_recoToSim_SC_score); + associations_tree_->Branch("Mergetstracksters_recoToSim_SC_sharedE", &MergeTracksters_recoToSim_SC_sharedE); + associations_tree_->Branch("Mergetstracksters_simToReco_SC", &MergeTracksters_simToReco_SC); + associations_tree_->Branch("Mergetstracksters_simToReco_SC_score", &MergeTracksters_simToReco_SC_score); + associations_tree_->Branch("Mergetstracksters_simToReco_SC_sharedE", &MergeTracksters_simToReco_SC_sharedE); + + associations_tree_->Branch("Mergetracksters_recoToSim_CP", &MergeTracksters_recoToSim_CP); + associations_tree_->Branch("Mergetracksters_recoToSim_CP_score", &MergeTracksters_recoToSim_CP_score); + associations_tree_->Branch("Mergetracksters_recoToSim_CP_sharedE", &MergeTracksters_recoToSim_CP_sharedE); + associations_tree_->Branch("Mergetracksters_simToReco_CP", &MergeTracksters_simToReco_CP); + associations_tree_->Branch("Mergetracksters_simToReco_CP_score", &MergeTracksters_simToReco_CP_score); + associations_tree_->Branch("Mergetracksters_simToReco_CP_sharedE", &MergeTracksters_simToReco_CP_sharedE); + + associations_tree_->Branch("Mergetracksters_recoToSim_PU", &MergeTracksters_recoToSim_PU); + associations_tree_->Branch("Mergetracksters_recoToSim_PU_score", &MergeTracksters_recoToSim_PU_score); + associations_tree_->Branch("Mergetracksters_recoToSim_PU_sharedE", &MergeTracksters_recoToSim_PU_sharedE); + associations_tree_->Branch("Mergetracksters_simToReco_PU", &MergeTracksters_simToReco_PU); + associations_tree_->Branch("Mergetracksters_simToReco_PU_score", &MergeTracksters_simToReco_PU_score); + associations_tree_->Branch("Mergetracksters_simToReco_PU_sharedE", &MergeTracksters_simToReco_PU_sharedE); + } + + if (saveSimTrackstersSC_) { + simtrackstersSC_tree_ = fs->make("simtrackstersSC", "TICL simTracksters SC"); + simtrackstersSC_tree_->Branch("event", &ev_event_); + simtrackstersSC_tree_->Branch("NTracksters", &stsSC_ntracksters_); + simtrackstersSC_tree_->Branch("time", &stsSC_trackster_time); + simtrackstersSC_tree_->Branch("timeError", &stsSC_trackster_timeError); + simtrackstersSC_tree_->Branch("regressed_energy", &stsSC_trackster_regressed_energy); + simtrackstersSC_tree_->Branch("regressed_pt", &stsSC_trackster_regressed_pt); + simtrackstersSC_tree_->Branch("raw_energy", &stsSC_trackster_raw_energy); + simtrackstersSC_tree_->Branch("raw_em_energy", &stsSC_trackster_raw_em_energy); + simtrackstersSC_tree_->Branch("raw_pt", &stsSC_trackster_raw_pt); + simtrackstersSC_tree_->Branch("raw_em_pt", &stsSC_trackster_raw_em_pt); + simtrackstersSC_tree_->Branch("barycenter_x", &stsSC_trackster_barycenter_x); + simtrackstersSC_tree_->Branch("barycenter_y", &stsSC_trackster_barycenter_y); + simtrackstersSC_tree_->Branch("barycenter_z", &stsSC_trackster_barycenter_z); + simtrackstersSC_tree_->Branch("barycenter_eta", &stsSC_trackster_barycenter_eta); + simtrackstersSC_tree_->Branch("barycenter_phi", &stsSC_trackster_barycenter_phi); + simtrackstersSC_tree_->Branch("EV1", &stsSC_trackster_EV1); + simtrackstersSC_tree_->Branch("EV2", &stsSC_trackster_EV2); + simtrackstersSC_tree_->Branch("EV3", &stsSC_trackster_EV3); + simtrackstersSC_tree_->Branch("eVector0_x", &stsSC_trackster_eVector0_x); + simtrackstersSC_tree_->Branch("eVector0_y", &stsSC_trackster_eVector0_y); + simtrackstersSC_tree_->Branch("eVector0_z", &stsSC_trackster_eVector0_z); + simtrackstersSC_tree_->Branch("sigmaPCA1", &stsSC_trackster_sigmaPCA1); + simtrackstersSC_tree_->Branch("sigmaPCA2", &stsSC_trackster_sigmaPCA2); + simtrackstersSC_tree_->Branch("sigmaPCA3", &stsSC_trackster_sigmaPCA3); + simtrackstersSC_tree_->Branch("pdgID", &stsSC_pdgID); + simtrackstersSC_tree_->Branch("trackIdx", &stsSC_trackIdx); + simtrackstersSC_tree_->Branch("trackTime", &stsSC_trackTime); + simtrackstersSC_tree_->Branch("boundaryX", &stsSC_boundaryX); + simtrackstersSC_tree_->Branch("boundaryY", &stsSC_boundaryY); + simtrackstersSC_tree_->Branch("boundaryZ", &stsSC_boundaryZ); + simtrackstersSC_tree_->Branch("boundaryEta", &stsSC_boundaryEta); + simtrackstersSC_tree_->Branch("boundaryPhi", &stsSC_boundaryPhi); + simtrackstersSC_tree_->Branch("boundaryPx", &stsSC_boundaryPx); + simtrackstersSC_tree_->Branch("boundaryPy", &stsSC_boundaryPy); + simtrackstersSC_tree_->Branch("boundaryPz", &stsSC_boundaryPz); + simtrackstersSC_tree_->Branch("track_boundaryX", &stsSC_track_boundaryX); + simtrackstersSC_tree_->Branch("track_boundaryY", &stsSC_track_boundaryY); + simtrackstersSC_tree_->Branch("track_boundaryZ", &stsSC_track_boundaryZ); + simtrackstersSC_tree_->Branch("track_boundaryEta", &stsSC_track_boundaryEta); + simtrackstersSC_tree_->Branch("track_boundaryPhi", &stsSC_track_boundaryPhi); + simtrackstersSC_tree_->Branch("track_boundaryPx", &stsSC_track_boundaryPx); + simtrackstersSC_tree_->Branch("track_boundaryPy", &stsSC_track_boundaryPy); + simtrackstersSC_tree_->Branch("track_boundaryPz", &stsSC_track_boundaryPz); + simtrackstersSC_tree_->Branch("id_probabilities", &stsSC_trackster_id_probabilities); + simtrackstersSC_tree_->Branch("vertices_indexes", &stsSC_trackster_vertices_indexes); + simtrackstersSC_tree_->Branch("vertices_x", &stsSC_trackster_vertices_x); + simtrackstersSC_tree_->Branch("vertices_y", &stsSC_trackster_vertices_y); + simtrackstersSC_tree_->Branch("vertices_z", &stsSC_trackster_vertices_z); + simtrackstersSC_tree_->Branch("vertices_time", &stsSC_trackster_vertices_time); + simtrackstersSC_tree_->Branch("vertices_timeErr", &stsSC_trackster_vertices_timeErr); + simtrackstersSC_tree_->Branch("vertices_energy", &stsSC_trackster_vertices_energy); + simtrackstersSC_tree_->Branch("vertices_correctedEnergy", &stsSC_trackster_vertices_correctedEnergy); + simtrackstersSC_tree_->Branch("vertices_correctedEnergyUncertainty", + &stsSC_trackster_vertices_correctedEnergyUncertainty); + simtrackstersSC_tree_->Branch("vertices_multiplicity", &stsSC_trackster_vertices_multiplicity); + simtrackstersSC_tree_->Branch("NsimTrackstersSC", &nsimTrackstersSC); + } + if (saveSimTrackstersCP_) { + simtrackstersCP_tree_ = fs->make("simtrackstersCP", "TICL simTracksters CP"); + simtrackstersCP_tree_->Branch("event", &ev_event_); + simtrackstersCP_tree_->Branch("NTracksters", &stsCP_ntracksters_); + simtrackstersCP_tree_->Branch("time", &stsCP_trackster_time); + simtrackstersCP_tree_->Branch("timeError", &stsCP_trackster_timeError); + simtrackstersCP_tree_->Branch("regressed_energy", &stsCP_trackster_regressed_energy); + simtrackstersCP_tree_->Branch("regressed_pt", &stsCP_trackster_regressed_pt); + simtrackstersCP_tree_->Branch("raw_energy", &stsCP_trackster_raw_energy); + simtrackstersCP_tree_->Branch("raw_em_energy", &stsCP_trackster_raw_em_energy); + simtrackstersCP_tree_->Branch("raw_pt", &stsCP_trackster_raw_pt); + simtrackstersCP_tree_->Branch("raw_em_pt", &stsCP_trackster_raw_em_pt); + simtrackstersCP_tree_->Branch("barycenter_x", &stsCP_trackster_barycenter_x); + simtrackstersCP_tree_->Branch("barycenter_y", &stsCP_trackster_barycenter_y); + simtrackstersCP_tree_->Branch("barycenter_z", &stsCP_trackster_barycenter_z); + simtrackstersCP_tree_->Branch("barycenter_eta", &stsCP_trackster_barycenter_eta); + simtrackstersCP_tree_->Branch("barycenter_phi", &stsCP_trackster_barycenter_phi); + simtrackstersCP_tree_->Branch("pdgID", &stsCP_pdgID); + simtrackstersCP_tree_->Branch("trackIdx", &stsCP_trackIdx); + simtrackstersCP_tree_->Branch("trackTime", &stsCP_trackTime); + simtrackstersCP_tree_->Branch("boundaryX", &stsCP_boundaryX); + simtrackstersCP_tree_->Branch("boundaryY", &stsCP_boundaryY); + simtrackstersCP_tree_->Branch("boundaryZ", &stsCP_boundaryZ); + simtrackstersCP_tree_->Branch("boundaryEta", &stsCP_boundaryEta); + simtrackstersCP_tree_->Branch("boundaryPhi", &stsCP_boundaryPhi); + simtrackstersCP_tree_->Branch("boundaryPx", &stsCP_boundaryPx); + simtrackstersCP_tree_->Branch("boundaryPy", &stsCP_boundaryPy); + simtrackstersCP_tree_->Branch("boundaryPz", &stsCP_boundaryPz); + simtrackstersCP_tree_->Branch("track_boundaryX", &stsCP_track_boundaryX); + simtrackstersCP_tree_->Branch("track_boundaryY", &stsCP_track_boundaryY); + simtrackstersCP_tree_->Branch("track_boundaryZ", &stsCP_track_boundaryZ); + simtrackstersCP_tree_->Branch("track_boundaryEta", &stsCP_track_boundaryEta); + simtrackstersCP_tree_->Branch("track_boundaryPhi", &stsCP_track_boundaryPhi); + simtrackstersCP_tree_->Branch("track_boundaryPx", &stsCP_track_boundaryPx); + simtrackstersCP_tree_->Branch("track_boundaryPy", &stsCP_track_boundaryPy); + simtrackstersCP_tree_->Branch("track_boundaryPz", &stsCP_track_boundaryPz); + simtrackstersCP_tree_->Branch("EV1", &stsCP_trackster_EV1); + simtrackstersCP_tree_->Branch("EV2", &stsCP_trackster_EV2); + simtrackstersCP_tree_->Branch("EV3", &stsCP_trackster_EV3); + simtrackstersCP_tree_->Branch("eVector0_x", &stsCP_trackster_eVector0_x); + simtrackstersCP_tree_->Branch("eVector0_y", &stsCP_trackster_eVector0_y); + simtrackstersCP_tree_->Branch("eVector0_z", &stsCP_trackster_eVector0_z); + simtrackstersCP_tree_->Branch("sigmaPCA1", &stsCP_trackster_sigmaPCA1); + simtrackstersCP_tree_->Branch("sigmaPCA2", &stsCP_trackster_sigmaPCA2); + simtrackstersCP_tree_->Branch("sigmaPCA3", &stsCP_trackster_sigmaPCA3); + simtrackstersCP_tree_->Branch("id_probabilities", &stsCP_trackster_id_probabilities); + simtrackstersCP_tree_->Branch("vertices_indexes", &stsCP_trackster_vertices_indexes); + simtrackstersCP_tree_->Branch("vertices_x", &stsCP_trackster_vertices_x); + simtrackstersCP_tree_->Branch("vertices_y", &stsCP_trackster_vertices_y); + simtrackstersCP_tree_->Branch("vertices_z", &stsCP_trackster_vertices_z); + simtrackstersCP_tree_->Branch("vertices_time", &stsCP_trackster_vertices_time); + simtrackstersCP_tree_->Branch("vertices_timeErr", &stsCP_trackster_vertices_timeErr); + simtrackstersCP_tree_->Branch("vertices_energy", &stsCP_trackster_vertices_energy); + simtrackstersCP_tree_->Branch("vertices_correctedEnergy", &stsCP_trackster_vertices_correctedEnergy); + simtrackstersCP_tree_->Branch("vertices_correctedEnergyUncertainty", + &stsCP_trackster_vertices_correctedEnergyUncertainty); + simtrackstersCP_tree_->Branch("vertices_multiplicity", &stsCP_trackster_vertices_multiplicity); + } + + if (saveTracks_) { + tracks_tree_ = fs->make("tracks", "Tracks"); + tracks_tree_->Branch("event", &ev_event_); + tracks_tree_->Branch("track_id", &track_id); + tracks_tree_->Branch("track_hgcal_pt", &track_hgcal_pt); + tracks_tree_->Branch("track_pt", &track_pt); + tracks_tree_->Branch("track_missing_outer_hits", &track_missing_outer_hits); + tracks_tree_->Branch("track_quality", &track_quality); + tracks_tree_->Branch("track_charge", &track_charge); + tracks_tree_->Branch("track_time", &track_time); + tracks_tree_->Branch("track_time_quality", &track_time_quality); + tracks_tree_->Branch("track_time_err", &track_time_err); + tracks_tree_->Branch("track_nhits", &track_nhits); + } + + if (saveSimTICLCandidate_) { + simTICLCandidate_tree = fs->make("simTICLCandidate", "Sim TICL Candidate"); + simTICLCandidate_tree->Branch("simTICLCandidate_raw_energy", &simTICLCandidate_raw_energy); + simTICLCandidate_tree->Branch("simTICLCandidate_regressed_energy", &simTICLCandidate_regressed_energy); + simTICLCandidate_tree->Branch("simTICLCandidate_simTracksterCPIndex", &simTICLCandidate_simTracksterCPIndex); + simTICLCandidate_tree->Branch("simTICLCandidate_boundaryX", &simTICLCandidate_boundaryX); + simTICLCandidate_tree->Branch("simTICLCandidate_boundaryY", &simTICLCandidate_boundaryY); + simTICLCandidate_tree->Branch("simTICLCandidate_boundaryZ", &simTICLCandidate_boundaryZ); + simTICLCandidate_tree->Branch("simTICLCandidate_boundaryPx", &simTICLCandidate_boundaryPx); + simTICLCandidate_tree->Branch("simTICLCandidate_boundaryPy", &simTICLCandidate_boundaryPy); + simTICLCandidate_tree->Branch("simTICLCandidate_boundaryPz", &simTICLCandidate_boundaryPz); + simTICLCandidate_tree->Branch("simTICLCandidate_trackTime", &simTICLCandidate_trackTime); + simTICLCandidate_tree->Branch("simTICLCandidate_trackBeta", &simTICLCandidate_trackBeta); + simTICLCandidate_tree->Branch("simTICLCandidate_caloParticleMass", &simTICLCandidate_caloParticleMass); + simTICLCandidate_tree->Branch("simTICLCandidate_pdgId", &simTICLCandidate_pdgId); + simTICLCandidate_tree->Branch("simTICLCandidate_charge", &simTICLCandidate_charge); + simTICLCandidate_tree->Branch("simTICLCandidate_track_in_candidate", &simTICLCandidate_track_in_candidate); + } +} + +void TICLDumper::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 TICLDumper::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; +} +void TICLDumper::analyze(const edm::Event& event, const edm::EventSetup& setup) { + ev_event_ += 1; + clearVariables(); + auto bFieldProd = bfield_.product(); + const Propagator& prop = (*propagator_); + //get all the tracksters + edm::Handle> tracksters_handle; + event.getByToken(tracksters_token_, tracksters_handle); + const auto& tracksters = *tracksters_handle; + + //get all the layer clusters + edm::Handle> layer_clusters_h; + event.getByToken(layer_clusters_token_, layer_clusters_h); + const auto& clusters = *layer_clusters_h; + + edm::Handle>> clustersTime_h; + event.getByToken(clustersTime_token_, clustersTime_h); + const auto& layerClustersTimes = *clustersTime_h; + + //TICL Candidate + edm::Handle> candidates_h; + event.getByToken(ticl_candidates_token_, candidates_h); + const auto& ticlcandidates = *candidates_h; + + //Track + edm::Handle> tracks_h; + event.getByToken(tracks_token_, tracks_h); + const auto& tracks = *tracks_h; + + edm::Handle> trackTime_h; + event.getByToken(tracks_time_token_, trackTime_h); + const auto& trackTime = *trackTime_h; + + edm::Handle> trackTimeErr_h; + event.getByToken(tracks_time_err_token_, trackTimeErr_h); + const auto& trackTimeErr = *trackTimeErr_h; + + edm::Handle> trackTimeQual_h; + event.getByToken(tracks_time_quality_token_, trackTimeQual_h); + const auto& trackTimeQual = *trackTimeQual_h; + + //Tracksters merged + edm::Handle> tracksters_merged_h; + event.getByToken(tracksters_merged_token_, tracksters_merged_h); + const auto& trackstersmerged = *tracksters_merged_h; + + // simTracksters from SC + edm::Handle> simTrackstersSC_h; + event.getByToken(simTracksters_SC_token_, simTrackstersSC_h); + const auto& simTrackstersSC = *simTrackstersSC_h; + + // simTracksters from CP + edm::Handle> simTrackstersCP_h; + event.getByToken(simTracksters_CP_token_, simTrackstersCP_h); + const auto& simTrackstersCP = *simTrackstersCP_h; + + // simTracksters from PU + edm::Handle> simTrackstersPU_h; + event.getByToken(simTracksters_PU_token_, simTrackstersPU_h); + const auto& simTrackstersPU = *simTrackstersPU_h; + + edm::Handle> simTICLCandidates_h; + event.getByToken(simTICLCandidate_token_, simTICLCandidates_h); + const auto& simTICLCandidates = *simTICLCandidates_h; + + // trackster reco to sim SC + edm::Handle tsRecoToSimSC_h; + event.getByToken(tsRecoToSimSC_token_, tsRecoToSimSC_h); + auto const& tsRecoSimSCMap = *tsRecoToSimSC_h; + + // sim simTrackster SC to reco trackster + edm::Handle tsSimToRecoSC_h; + event.getByToken(tsSimToRecoSC_token_, tsSimToRecoSC_h); + auto const& tsSimToRecoSCMap = *tsSimToRecoSC_h; + + // trackster reco to sim CP + edm::Handle tsRecoToSimCP_h; + event.getByToken(tsRecoToSimCP_token_, tsRecoToSimCP_h); + auto const& tsRecoSimCPMap = *tsRecoToSimCP_h; + + // sim simTrackster CP to reco trackster + edm::Handle tsSimToRecoCP_h; + event.getByToken(tsSimToRecoCP_token_, tsSimToRecoCP_h); + auto const& tsSimToRecoCPMap = *tsSimToRecoCP_h; + + edm::Handle mergetsRecoToSimSC_h; + event.getByToken(MergeRecoToSimSC_token_, mergetsRecoToSimSC_h); + auto const& MergetsRecoSimSCMap = *mergetsRecoToSimSC_h; + + // sim simTrackster SC to reco trackster + edm::Handle mergetsSimToRecoSC_h; + event.getByToken(MergeSimToRecoSC_token_, mergetsSimToRecoSC_h); + auto const& MergetsSimToRecoSCMap = *mergetsSimToRecoSC_h; + + // trackster reco to sim CP + edm::Handle mergetsRecoToSimCP_h; + event.getByToken(MergeRecoToSimCP_token_, mergetsRecoToSimCP_h); + auto const& MergetsRecoSimCPMap = *mergetsRecoToSimCP_h; + + // sim simTrackster CP to reco trackster + edm::Handle mergetsSimToRecoCP_h; + event.getByToken(MergeSimToRecoCP_token_, mergetsSimToRecoCP_h); + auto const& MergetsSimToRecoCPMap = *mergetsSimToRecoCP_h; + + // trackster reco to sim PU + edm::Handle mergetsRecoToSimPU_h; + event.getByToken(MergeRecoToSimPU_token_, mergetsRecoToSimPU_h); + auto const& MergetsRecoSimPUMap = *mergetsRecoToSimPU_h; + + // sim simTrackster PU to reco trackster + edm::Handle mergetsSimToRecoPU_h; + event.getByToken(MergeSimToRecoPU_token_, mergetsSimToRecoPU_h); + auto const& MergetsSimToRecoPUMap = *mergetsSimToRecoPU_h; + + edm::Handle> caloparticles_h; + event.getByToken(caloparticles_token_, caloparticles_h); + const auto& caloparticles = *caloparticles_h; + + const auto& simclusters = event.get(simclusters_token_); + + ntracksters_ = tracksters.size(); + nclusters_ = clusters.size(); + + for (auto trackster_iterator = tracksters.begin(); trackster_iterator != tracksters.end(); ++trackster_iterator) { + //per-trackster analysis + trackster_time.push_back(trackster_iterator->time()); + trackster_timeError.push_back(trackster_iterator->timeError()); + trackster_regressed_energy.push_back(trackster_iterator->regressed_energy()); + trackster_raw_energy.push_back(trackster_iterator->raw_energy()); + trackster_raw_em_energy.push_back(trackster_iterator->raw_em_energy()); + trackster_raw_pt.push_back(trackster_iterator->raw_pt()); + trackster_raw_em_pt.push_back(trackster_iterator->raw_em_pt()); + trackster_barycenter_x.push_back(trackster_iterator->barycenter().x()); + trackster_barycenter_y.push_back(trackster_iterator->barycenter().y()); + trackster_barycenter_z.push_back(trackster_iterator->barycenter().z()); + trackster_barycenter_eta.push_back(trackster_iterator->barycenter().eta()); + trackster_barycenter_phi.push_back(trackster_iterator->barycenter().phi()); + trackster_EV1.push_back(trackster_iterator->eigenvalues()[0]); + trackster_EV2.push_back(trackster_iterator->eigenvalues()[1]); + trackster_EV3.push_back(trackster_iterator->eigenvalues()[2]); + trackster_eVector0_x.push_back((trackster_iterator->eigenvectors()[0]).x()); + trackster_eVector0_y.push_back((trackster_iterator->eigenvectors()[0]).y()); + trackster_eVector0_z.push_back((trackster_iterator->eigenvectors()[0]).z()); + trackster_sigmaPCA1.push_back(trackster_iterator->sigmasPCA()[0]); + trackster_sigmaPCA2.push_back(trackster_iterator->sigmasPCA()[1]); + trackster_sigmaPCA3.push_back(trackster_iterator->sigmasPCA()[2]); + std::vector id_probs; + for (size_t i = 0; i < 8; i++) + id_probs.push_back(trackster_iterator->id_probabilities(i)); + trackster_id_probabilities.push_back(id_probs); + + // Clusters + std::vector vertices_indexes; + std::vector vertices_x; + std::vector vertices_y; + std::vector vertices_z; + std::vector vertices_time; + std::vector vertices_timeErr; + std::vector vertices_energy; + std::vector vertices_correctedEnergy; + std::vector vertices_correctedEnergyUncertainty; + for (auto idx : trackster_iterator->vertices()) { + vertices_indexes.push_back(idx); + auto associated_cluster = (*layer_clusters_h)[idx]; + vertices_x.push_back(associated_cluster.x()); + vertices_y.push_back(associated_cluster.y()); + vertices_z.push_back(associated_cluster.z()); + vertices_energy.push_back(associated_cluster.energy()); + vertices_correctedEnergy.push_back(associated_cluster.correctedEnergy()); + vertices_correctedEnergyUncertainty.push_back(associated_cluster.correctedEnergyUncertainty()); + vertices_time.push_back(layerClustersTimes.get(idx).first); + vertices_timeErr.push_back(layerClustersTimes.get(idx).second); + } + trackster_vertices_indexes.push_back(vertices_indexes); + trackster_vertices_x.push_back(vertices_x); + trackster_vertices_y.push_back(vertices_y); + trackster_vertices_z.push_back(vertices_z); + trackster_vertices_time.push_back(vertices_time); + trackster_vertices_timeErr.push_back(vertices_timeErr); + trackster_vertices_energy.push_back(vertices_energy); + trackster_vertices_correctedEnergy.push_back(vertices_correctedEnergy); + trackster_vertices_correctedEnergyUncertainty.push_back(vertices_correctedEnergyUncertainty); + + // Multiplicity + std::vector vertices_multiplicity; + for (auto multiplicity : trackster_iterator->vertex_multiplicity()) { + vertices_multiplicity.push_back(multiplicity); + } + trackster_vertices_multiplicity.push_back(vertices_multiplicity); + } + + stsSC_ntracksters_ = simTrackstersSC.size(); + using CaloObjectVariant = std::variant; + for (auto trackster_iterator = simTrackstersSC.begin(); trackster_iterator != simTrackstersSC.end(); + ++trackster_iterator) { + //per-trackster analysis + stsSC_trackster_time.push_back(trackster_iterator->time()); + stsSC_trackster_timeError.push_back(trackster_iterator->timeError()); + stsSC_trackster_regressed_energy.push_back(trackster_iterator->regressed_energy()); + stsSC_trackster_raw_energy.push_back(trackster_iterator->raw_energy()); + stsSC_trackster_raw_em_energy.push_back(trackster_iterator->raw_em_energy()); + stsSC_trackster_raw_pt.push_back(trackster_iterator->raw_pt()); + stsSC_trackster_raw_em_pt.push_back(trackster_iterator->raw_em_pt()); + stsSC_trackster_barycenter_x.push_back(trackster_iterator->barycenter().x()); + stsSC_trackster_barycenter_y.push_back(trackster_iterator->barycenter().y()); + stsSC_trackster_barycenter_z.push_back(trackster_iterator->barycenter().z()); + stsSC_trackster_barycenter_eta.push_back(trackster_iterator->barycenter().eta()); + stsSC_trackster_barycenter_phi.push_back(trackster_iterator->barycenter().phi()); + stsSC_trackster_EV1.push_back(trackster_iterator->eigenvalues()[0]); + stsSC_trackster_EV2.push_back(trackster_iterator->eigenvalues()[1]); + stsSC_trackster_EV3.push_back(trackster_iterator->eigenvalues()[2]); + stsSC_trackster_eVector0_x.push_back((trackster_iterator->eigenvectors()[0]).x()); + stsSC_trackster_eVector0_y.push_back((trackster_iterator->eigenvectors()[0]).y()); + stsSC_trackster_eVector0_z.push_back((trackster_iterator->eigenvectors()[0]).z()); + stsSC_trackster_sigmaPCA1.push_back(trackster_iterator->sigmasPCA()[0]); + stsSC_trackster_sigmaPCA2.push_back(trackster_iterator->sigmasPCA()[1]); + stsSC_trackster_sigmaPCA3.push_back(trackster_iterator->sigmasPCA()[2]); + stsSC_pdgID.push_back(simclusters[trackster_iterator->seedIndex()].pdgId()); + + CaloObjectVariant caloObj; + if (trackster_iterator->seedID() == caloparticles_h.id()) { + caloObj = caloparticles[trackster_iterator->seedIndex()]; + } else { + caloObj = simclusters[trackster_iterator->seedIndex()]; + } + + auto const& simTrack = std::visit([](auto&& obj) { return obj.g4Tracks()[0]; }, caloObj); + auto const& caloPt = std::visit([](auto&& obj) { return obj.pt(); }, caloObj); + stsSC_trackster_regressed_pt.push_back(caloPt); + if (simTrack.crossedBoundary()) { + stsSC_boundaryX.push_back(simTrack.getPositionAtBoundary().x()); + stsSC_boundaryY.push_back(simTrack.getPositionAtBoundary().y()); + stsSC_boundaryZ.push_back(simTrack.getPositionAtBoundary().z()); + stsSC_boundaryEta.push_back(simTrack.getPositionAtBoundary().eta()); + stsSC_boundaryPhi.push_back(simTrack.getPositionAtBoundary().phi()); + stsSC_boundaryPx.push_back(simTrack.getMomentumAtBoundary().x()); + stsSC_boundaryPy.push_back(simTrack.getMomentumAtBoundary().y()); + stsSC_boundaryPz.push_back(simTrack.getMomentumAtBoundary().z()); + } else { + stsSC_boundaryX.push_back(-999); + stsSC_boundaryY.push_back(-999); + stsSC_boundaryZ.push_back(-999); + stsSC_boundaryEta.push_back(-999); + stsSC_boundaryPhi.push_back(-999); + stsSC_boundaryPx.push_back(-999); + stsSC_boundaryPy.push_back(-999); + stsSC_boundaryPz.push_back(-999); + } + auto const trackIdx = trackster_iterator->trackIdx(); + stsSC_trackIdx.push_back(trackIdx); + if (trackIdx != -1) { + const auto& track = tracks[trackIdx]; + + int iSide = int(track.eta() > 0); + + const auto& fts = trajectoryStateTransform::outerFreeState((track), bFieldProd); + // to the HGCal front + const auto& tsos = prop.propagate(fts, firstDisk_[iSide]->surface()); + if (tsos.isValid()) { + const auto& globalPos = tsos.globalPosition(); + const auto& globalMom = tsos.globalMomentum(); + stsSC_track_boundaryX.push_back(globalPos.x()); + stsSC_track_boundaryY.push_back(globalPos.y()); + stsSC_track_boundaryZ.push_back(globalPos.z()); + stsSC_track_boundaryEta.push_back(globalPos.eta()); + stsSC_track_boundaryPhi.push_back(globalPos.phi()); + stsSC_track_boundaryPx.push_back(globalMom.x()); + stsSC_track_boundaryPy.push_back(globalMom.y()); + stsSC_track_boundaryPz.push_back(globalMom.z()); + stsSC_trackTime.push_back(track.t0()); + } else { + stsSC_track_boundaryX.push_back(-999); + stsSC_track_boundaryY.push_back(-999); + stsSC_track_boundaryZ.push_back(-999); + stsSC_track_boundaryEta.push_back(-999); + stsSC_track_boundaryPhi.push_back(-999); + stsSC_track_boundaryPx.push_back(-999); + stsSC_track_boundaryPy.push_back(-999); + stsSC_track_boundaryPz.push_back(-999); + } + } else { + stsSC_track_boundaryX.push_back(-999); + stsSC_track_boundaryY.push_back(-999); + stsSC_track_boundaryZ.push_back(-999); + stsSC_track_boundaryEta.push_back(-999); + stsSC_track_boundaryPhi.push_back(-999); + stsSC_track_boundaryPx.push_back(-999); + stsSC_track_boundaryPy.push_back(-999); + stsSC_track_boundaryPz.push_back(-999); + } + + std::vector id_probs; + for (size_t i = 0; i < 8; i++) + id_probs.push_back(trackster_iterator->id_probabilities(i)); + stsSC_trackster_id_probabilities.push_back(id_probs); + + // Clusters + std::vector vertices_indexes; + std::vector vertices_x; + std::vector vertices_y; + std::vector vertices_z; + std::vector vertices_time; + std::vector vertices_timeErr; + std::vector vertices_energy; + std::vector vertices_correctedEnergy; + std::vector vertices_correctedEnergyUncertainty; + for (auto idx : trackster_iterator->vertices()) { + vertices_indexes.push_back(idx); + auto associated_cluster = (*layer_clusters_h)[idx]; + vertices_x.push_back(associated_cluster.x()); + vertices_y.push_back(associated_cluster.y()); + vertices_z.push_back(associated_cluster.z()); + vertices_energy.push_back(associated_cluster.energy()); + vertices_correctedEnergy.push_back(associated_cluster.correctedEnergy()); + vertices_correctedEnergyUncertainty.push_back(associated_cluster.correctedEnergyUncertainty()); + vertices_time.push_back(layerClustersTimes.get(idx).first); + vertices_timeErr.push_back(layerClustersTimes.get(idx).second); + } + stsSC_trackster_vertices_indexes.push_back(vertices_indexes); + stsSC_trackster_vertices_x.push_back(vertices_x); + stsSC_trackster_vertices_y.push_back(vertices_y); + stsSC_trackster_vertices_z.push_back(vertices_z); + stsSC_trackster_vertices_time.push_back(vertices_time); + stsSC_trackster_vertices_timeErr.push_back(vertices_timeErr); + stsSC_trackster_vertices_energy.push_back(vertices_energy); + stsSC_trackster_vertices_correctedEnergy.push_back(vertices_correctedEnergy); + stsSC_trackster_vertices_correctedEnergyUncertainty.push_back(vertices_correctedEnergyUncertainty); + + // Multiplicity + std::vector vertices_multiplicity; + for (auto multiplicity : trackster_iterator->vertex_multiplicity()) { + vertices_multiplicity.push_back(multiplicity); + } + stsSC_trackster_vertices_multiplicity.push_back(vertices_multiplicity); + } + + stsCP_ntracksters_ = simTrackstersCP.size(); + + for (auto trackster_iterator = simTrackstersCP.begin(); trackster_iterator != simTrackstersCP.end(); + ++trackster_iterator) { + //per-trackster analysis + stsCP_trackster_time.push_back(trackster_iterator->time()); + stsCP_trackster_timeError.push_back(trackster_iterator->timeError()); + stsCP_trackster_regressed_energy.push_back(trackster_iterator->regressed_energy()); + stsCP_trackster_raw_energy.push_back(trackster_iterator->raw_energy()); + stsCP_trackster_raw_em_energy.push_back(trackster_iterator->raw_em_energy()); + stsCP_trackster_raw_pt.push_back(trackster_iterator->raw_pt()); + stsCP_trackster_raw_em_pt.push_back(trackster_iterator->raw_em_pt()); + stsCP_trackster_barycenter_x.push_back(trackster_iterator->barycenter().x()); + stsCP_trackster_barycenter_y.push_back(trackster_iterator->barycenter().y()); + stsCP_trackster_barycenter_z.push_back(trackster_iterator->barycenter().z()); + stsCP_trackster_barycenter_eta.push_back(trackster_iterator->barycenter().eta()); + stsCP_trackster_barycenter_phi.push_back(trackster_iterator->barycenter().phi()); + stsCP_trackster_EV1.push_back(trackster_iterator->eigenvalues()[0]); + stsCP_trackster_EV2.push_back(trackster_iterator->eigenvalues()[1]); + stsCP_trackster_EV3.push_back(trackster_iterator->eigenvalues()[2]); + stsCP_trackster_eVector0_x.push_back((trackster_iterator->eigenvectors()[0]).x()); + stsCP_trackster_eVector0_y.push_back((trackster_iterator->eigenvectors()[0]).y()); + stsCP_trackster_eVector0_z.push_back((trackster_iterator->eigenvectors()[0]).z()); + stsCP_trackster_sigmaPCA1.push_back(trackster_iterator->sigmasPCA()[0]); + stsCP_trackster_sigmaPCA2.push_back(trackster_iterator->sigmasPCA()[1]); + stsCP_trackster_sigmaPCA3.push_back(trackster_iterator->sigmasPCA()[2]); + stsCP_pdgID.push_back(caloparticles[trackster_iterator->seedIndex()].pdgId()); + CaloObjectVariant caloObj; + if (trackster_iterator->seedID() == caloparticles_h.id()) { + caloObj = caloparticles[trackster_iterator->seedIndex()]; + } else { + caloObj = simclusters[trackster_iterator->seedIndex()]; + } + + auto const& simTrack = std::visit([](auto&& obj) { return obj.g4Tracks()[0]; }, caloObj); + auto const& caloPt = std::visit([](auto&& obj) { return obj.pt(); }, caloObj); + stsCP_trackster_regressed_pt.push_back(caloPt); + + if (simTrack.crossedBoundary()) { + stsCP_boundaryX.push_back(simTrack.getPositionAtBoundary().x()); + stsCP_boundaryY.push_back(simTrack.getPositionAtBoundary().y()); + stsCP_boundaryZ.push_back(simTrack.getPositionAtBoundary().z()); + stsCP_boundaryEta.push_back(simTrack.getPositionAtBoundary().eta()); + stsCP_boundaryPhi.push_back(simTrack.getPositionAtBoundary().phi()); + stsCP_boundaryPx.push_back(simTrack.getMomentumAtBoundary().x()); + stsCP_boundaryPy.push_back(simTrack.getMomentumAtBoundary().y()); + stsCP_boundaryPz.push_back(simTrack.getMomentumAtBoundary().z()); + } else { + stsCP_boundaryX.push_back(-999); + stsCP_boundaryY.push_back(-999); + stsCP_boundaryZ.push_back(-999); + stsCP_boundaryEta.push_back(-999); + stsCP_boundaryPhi.push_back(-999); + stsCP_boundaryPx.push_back(-999); + stsCP_boundaryPy.push_back(-999); + stsCP_boundaryPz.push_back(-999); + } + auto const trackIdx = trackster_iterator->trackIdx(); + stsCP_trackIdx.push_back(trackIdx); + if (trackIdx != -1) { + const auto& track = tracks[trackIdx]; + + int iSide = int(track.eta() > 0); + + const auto& fts = trajectoryStateTransform::outerFreeState((track), bFieldProd); + // to the HGCal front + const auto& tsos = prop.propagate(fts, firstDisk_[iSide]->surface()); + if (tsos.isValid()) { + const auto& globalPos = tsos.globalPosition(); + const auto& globalMom = tsos.globalMomentum(); + stsCP_track_boundaryX.push_back(globalPos.x()); + stsCP_track_boundaryY.push_back(globalPos.y()); + stsCP_track_boundaryZ.push_back(globalPos.z()); + stsCP_track_boundaryEta.push_back(globalPos.eta()); + stsCP_track_boundaryPhi.push_back(globalPos.phi()); + stsCP_track_boundaryPx.push_back(globalMom.x()); + stsCP_track_boundaryPy.push_back(globalMom.y()); + stsCP_track_boundaryPz.push_back(globalMom.z()); + stsCP_trackTime.push_back(track.t0()); + } else { + stsCP_track_boundaryX.push_back(-999); + stsCP_track_boundaryY.push_back(-999); + stsCP_track_boundaryZ.push_back(-999); + stsCP_track_boundaryEta.push_back(-999); + stsCP_track_boundaryPhi.push_back(-999); + stsCP_track_boundaryPx.push_back(-999); + stsCP_track_boundaryPy.push_back(-999); + stsCP_track_boundaryPz.push_back(-999); + } + } else { + stsCP_track_boundaryX.push_back(-999); + stsCP_track_boundaryY.push_back(-999); + stsCP_track_boundaryZ.push_back(-999); + stsCP_track_boundaryEta.push_back(-999); + stsCP_track_boundaryPhi.push_back(-999); + stsCP_track_boundaryPx.push_back(-999); + stsCP_track_boundaryPy.push_back(-999); + stsCP_track_boundaryPz.push_back(-999); + } + std::vector id_probs; + for (size_t i = 0; i < 8; i++) + id_probs.push_back(trackster_iterator->id_probabilities(i)); + stsCP_trackster_id_probabilities.push_back(id_probs); + + // Clusters + std::vector vertices_indexes; + std::vector vertices_x; + std::vector vertices_y; + std::vector vertices_z; + std::vector vertices_time; + std::vector vertices_timeErr; + std::vector vertices_energy; + std::vector vertices_correctedEnergy; + std::vector vertices_correctedEnergyUncertainty; + for (auto idx : trackster_iterator->vertices()) { + vertices_indexes.push_back(idx); + auto associated_cluster = (*layer_clusters_h)[idx]; + vertices_x.push_back(associated_cluster.x()); + vertices_y.push_back(associated_cluster.y()); + vertices_z.push_back(associated_cluster.z()); + vertices_energy.push_back(associated_cluster.energy()); + vertices_correctedEnergy.push_back(associated_cluster.correctedEnergy()); + vertices_correctedEnergyUncertainty.push_back(associated_cluster.correctedEnergyUncertainty()); + vertices_time.push_back(layerClustersTimes.get(idx).first); + vertices_timeErr.push_back(layerClustersTimes.get(idx).second); + } + stsCP_trackster_vertices_indexes.push_back(vertices_indexes); + stsCP_trackster_vertices_x.push_back(vertices_x); + stsCP_trackster_vertices_y.push_back(vertices_y); + stsCP_trackster_vertices_z.push_back(vertices_z); + stsCP_trackster_vertices_time.push_back(vertices_time); + stsCP_trackster_vertices_timeErr.push_back(vertices_timeErr); + stsCP_trackster_vertices_energy.push_back(vertices_energy); + stsCP_trackster_vertices_correctedEnergy.push_back(vertices_correctedEnergy); + stsCP_trackster_vertices_correctedEnergyUncertainty.push_back(vertices_correctedEnergyUncertainty); + + // Multiplicity + std::vector vertices_multiplicity; + for (auto multiplicity : trackster_iterator->vertex_multiplicity()) { + vertices_multiplicity.push_back(multiplicity); + } + stsCP_trackster_vertices_multiplicity.push_back(vertices_multiplicity); + } + + simTICLCandidate_track_in_candidate.resize(simTICLCandidates.size(), -1); + for (size_t i = 0; i < simTICLCandidates.size(); ++i) { + auto const& cand = simTICLCandidates[i]; + + simTICLCandidate_raw_energy.push_back(cand.rawEnergy()); + simTICLCandidate_regressed_energy.push_back(cand.p4().energy()); + simTICLCandidate_pdgId.push_back(cand.pdgId()); + simTICLCandidate_charge.push_back(cand.charge()); + std::vector tmpIdxVec; + for (auto const& simTS : cand.tracksters()) { + auto trackster_idx = simTS.get() - (edm::Ptr(simTrackstersSC_h, 0)).get(); + tmpIdxVec.push_back(trackster_idx); + } + simTICLCandidate_simTracksterCPIndex.push_back(tmpIdxVec); + tmpIdxVec.clear(); + auto const& trackPtr = cand.trackPtr(); + if (!trackPtr.isNull()) { + auto const& track = *trackPtr; + int iSide = int(track.eta() > 0); + int tk_idx = trackPtr.get() - (edm::Ptr(tracks_h, 0)).get(); + simTICLCandidate_track_in_candidate[i] = tk_idx; + + const auto& fts = trajectoryStateTransform::outerFreeState((track), bFieldProd); + // to the HGCal front + const auto& tsos = prop.propagate(fts, firstDisk_[iSide]->surface()); + if (tsos.isValid()) { + const auto& globalPos = tsos.globalPosition(); + const auto& globalMom = tsos.globalMomentum(); + simTICLCandidate_boundaryX.push_back(globalPos.x()); + simTICLCandidate_boundaryY.push_back(globalPos.y()); + simTICLCandidate_boundaryZ.push_back(globalPos.z()); + simTICLCandidate_boundaryPx.push_back(globalMom.x()); + simTICLCandidate_boundaryPy.push_back(globalMom.y()); + simTICLCandidate_boundaryPz.push_back(globalMom.z()); + simTICLCandidate_trackTime.push_back(track.t0()); + simTICLCandidate_trackBeta.push_back(track.beta()); + } else { + simTICLCandidate_boundaryX.push_back(-999); + simTICLCandidate_boundaryY.push_back(-999); + simTICLCandidate_boundaryZ.push_back(-999); + simTICLCandidate_boundaryPx.push_back(-999); + simTICLCandidate_boundaryPy.push_back(-999); + simTICLCandidate_boundaryPz.push_back(-999); + simTICLCandidate_trackTime.push_back(-999); + simTICLCandidate_trackBeta.push_back(-999); + } + } else { + simTICLCandidate_boundaryX.push_back(-999); + simTICLCandidate_boundaryY.push_back(-999); + simTICLCandidate_boundaryZ.push_back(-999); + simTICLCandidate_boundaryPx.push_back(-999); + simTICLCandidate_boundaryPy.push_back(-999); + simTICLCandidate_boundaryPz.push_back(-999); + simTICLCandidate_trackTime.push_back(-999); + simTICLCandidate_trackBeta.push_back(-999); + } + } + + int c_id = 0; + + for (auto cluster_iterator = clusters.begin(); cluster_iterator != clusters.end(); ++cluster_iterator) { + auto lc_seed = cluster_iterator->seed(); + cluster_seedID.push_back(lc_seed); + cluster_energy.push_back(cluster_iterator->energy()); + cluster_correctedEnergy.push_back(cluster_iterator->correctedEnergy()); + cluster_correctedEnergyUncertainty.push_back(cluster_iterator->correctedEnergyUncertainty()); + cluster_position_x.push_back(cluster_iterator->x()); + cluster_position_y.push_back(cluster_iterator->y()); + cluster_position_z.push_back(cluster_iterator->z()); + cluster_position_eta.push_back(cluster_iterator->eta()); + cluster_position_phi.push_back(cluster_iterator->phi()); + auto haf = cluster_iterator->hitsAndFractions(); + auto layerId = rhtools_.getLayerWithOffset(haf[0].first); + cluster_layer_id.push_back(layerId); + uint32_t number_of_hits = cluster_iterator->hitsAndFractions().size(); + cluster_number_of_hits.push_back(number_of_hits); + cluster_type.push_back(ticl::returnClusterType(lc_seed, rhtools_)); + cluster_timeErr.push_back(layerClustersTimes.get(c_id).second); + cluster_time.push_back(layerClustersTimes.get(c_id).first); + c_id += 1; + } + + tracksters_in_candidate.resize(ticlcandidates.size()); + track_in_candidate.resize(ticlcandidates.size(), -1); + nCandidates = ticlcandidates.size(); + for (int i = 0; i < static_cast(ticlcandidates.size()); ++i) { + const auto& candidate = ticlcandidates[i]; + candidate_charge.push_back(candidate.charge()); + candidate_pdgId.push_back(candidate.pdgId()); + candidate_energy.push_back(candidate.energy()); + candidate_px.push_back(candidate.px()); + candidate_py.push_back(candidate.py()); + candidate_pz.push_back(candidate.pz()); + candidate_time.push_back(candidate.time()); + candidate_time_err.push_back(candidate.timeError()); + std::vector id_probs; + for (int j = 0; j < 8; j++) { + ticl::Trackster::ParticleType type = static_cast(j); + id_probs.push_back(candidate.id_probability(type)); + } + candidate_id_probabilities.push_back(id_probs); + + auto trackster_ptrs = candidate.tracksters(); + auto track_ptr = candidate.trackPtr(); + for (const auto& ts_ptr : trackster_ptrs) { + auto ts_idx = ts_ptr.get() - (edm::Ptr(tracksters_handle, 0)).get(); + tracksters_in_candidate[i].push_back(ts_idx); + } + if (track_ptr.isNull()) + continue; + int tk_idx = track_ptr.get() - (edm::Ptr(tracks_h, 0)).get(); + track_in_candidate[i] = tk_idx; + } + + nTrackstersMerged = trackstersmerged.size(); + for (auto trackster_iterator = trackstersmerged.begin(); trackster_iterator != trackstersmerged.end(); + ++trackster_iterator) { + tracksters_merged_time.push_back(trackster_iterator->time()); + tracksters_merged_timeError.push_back(trackster_iterator->timeError()); + tracksters_merged_regressed_energy.push_back(trackster_iterator->regressed_energy()); + tracksters_merged_raw_energy.push_back(trackster_iterator->raw_energy()); + tracksters_merged_raw_em_energy.push_back(trackster_iterator->raw_em_energy()); + tracksters_merged_raw_pt.push_back(trackster_iterator->raw_pt()); + tracksters_merged_raw_em_pt.push_back(trackster_iterator->raw_em_pt()); + tracksters_merged_barycenter_x.push_back(trackster_iterator->barycenter().x()); + tracksters_merged_barycenter_y.push_back(trackster_iterator->barycenter().y()); + tracksters_merged_barycenter_z.push_back(trackster_iterator->barycenter().z()); + tracksters_merged_barycenter_eta.push_back(trackster_iterator->barycenter().eta()); + tracksters_merged_barycenter_phi.push_back(trackster_iterator->barycenter().phi()); + tracksters_merged_EV1.push_back(trackster_iterator->eigenvalues()[0]); + tracksters_merged_EV2.push_back(trackster_iterator->eigenvalues()[1]); + tracksters_merged_EV3.push_back(trackster_iterator->eigenvalues()[2]); + tracksters_merged_eVector0_x.push_back((trackster_iterator->eigenvectors()[0]).x()); + tracksters_merged_eVector0_y.push_back((trackster_iterator->eigenvectors()[0]).y()); + tracksters_merged_eVector0_z.push_back((trackster_iterator->eigenvectors()[0]).z()); + tracksters_merged_sigmaPCA1.push_back(trackster_iterator->sigmasPCA()[0]); + tracksters_merged_sigmaPCA2.push_back(trackster_iterator->sigmasPCA()[1]); + tracksters_merged_sigmaPCA3.push_back(trackster_iterator->sigmasPCA()[2]); + + std::vector id_probs; + for (size_t i = 0; i < 8; i++) + id_probs.push_back(trackster_iterator->id_probabilities(i)); + tracksters_merged_id_probabilities.push_back(id_probs); + + std::vector vertices_indexes; + std::vector vertices_x; + std::vector vertices_y; + std::vector vertices_z; + std::vector vertices_time; + std::vector vertices_timeErr; + std::vector vertices_energy; + std::vector vertices_correctedEnergy; + std::vector vertices_correctedEnergyUncertainty; + for (auto idx : trackster_iterator->vertices()) { + vertices_indexes.push_back(idx); + auto associated_cluster = (*layer_clusters_h)[idx]; + vertices_x.push_back(associated_cluster.x()); + vertices_y.push_back(associated_cluster.y()); + vertices_z.push_back(associated_cluster.z()); + vertices_energy.push_back(associated_cluster.energy()); + vertices_correctedEnergy.push_back(associated_cluster.correctedEnergy()); + vertices_correctedEnergyUncertainty.push_back(associated_cluster.correctedEnergyUncertainty()); + vertices_time.push_back(layerClustersTimes.get(idx).first); + vertices_timeErr.push_back(layerClustersTimes.get(idx).second); + } + tracksters_merged_vertices_indexes.push_back(vertices_indexes); + tracksters_merged_vertices_x.push_back(vertices_x); + tracksters_merged_vertices_y.push_back(vertices_y); + tracksters_merged_vertices_z.push_back(vertices_z); + tracksters_merged_vertices_time.push_back(vertices_time); + tracksters_merged_vertices_timeErr.push_back(vertices_timeErr); + tracksters_merged_vertices_energy.push_back(vertices_energy); + tracksters_merged_vertices_correctedEnergy.push_back(vertices_correctedEnergy); + tracksters_merged_vertices_correctedEnergyUncertainty.push_back(vertices_correctedEnergyUncertainty); + } + + // Tackster reco->sim associations + trackstersCLUE3D_recoToSim_SC.resize(tracksters.size()); + trackstersCLUE3D_recoToSim_SC_score.resize(tracksters.size()); + trackstersCLUE3D_recoToSim_SC_sharedE.resize(tracksters.size()); + for (size_t i = 0; i < tracksters.size(); ++i) { + const edm::Ref tsRef(tracksters_handle, i); + + // CLUE3D -> STS-SC + const auto stsSC_iter = tsRecoSimSCMap.find(tsRef); + if (stsSC_iter != tsRecoSimSCMap.end()) { + const auto& stsSCassociated = stsSC_iter->val; + for (auto& sts : stsSCassociated) { + auto sts_id = (sts.first).get() - (edm::Ref(simTrackstersSC_h, 0)).get(); + trackstersCLUE3D_recoToSim_SC[i].push_back(sts_id); + trackstersCLUE3D_recoToSim_SC_score[i].push_back(sts.second.second); + trackstersCLUE3D_recoToSim_SC_sharedE[i].push_back(sts.second.first); + } + } + } + + // SimTracksters + nsimTrackstersSC = simTrackstersSC.size(); + trackstersCLUE3D_simToReco_SC.resize(nsimTrackstersSC); + trackstersCLUE3D_simToReco_SC_score.resize(nsimTrackstersSC); + trackstersCLUE3D_simToReco_SC_sharedE.resize(nsimTrackstersSC); + for (size_t i = 0; i < nsimTrackstersSC; ++i) { + const edm::Ref stsSCRef(simTrackstersSC_h, i); + + // STS-SC -> CLUE3D + const auto ts_iter = tsSimToRecoSCMap.find(stsSCRef); + if (ts_iter != tsSimToRecoSCMap.end()) { + const auto& tsAssociated = ts_iter->val; + for (auto& ts : tsAssociated) { + auto ts_idx = (ts.first).get() - (edm::Ref(tracksters_handle, 0)).get(); + trackstersCLUE3D_simToReco_SC[i].push_back(ts_idx); + trackstersCLUE3D_simToReco_SC_score[i].push_back(ts.second.second); + trackstersCLUE3D_simToReco_SC_sharedE[i].push_back(ts.second.first); + } + } + } + + // Tackster reco->sim associations + trackstersCLUE3D_recoToSim_CP.resize(tracksters.size()); + trackstersCLUE3D_recoToSim_CP_score.resize(tracksters.size()); + trackstersCLUE3D_recoToSim_CP_sharedE.resize(tracksters.size()); + for (size_t i = 0; i < tracksters.size(); ++i) { + const edm::Ref tsRef(tracksters_handle, i); + + // CLUE3D -> STS-CP + const auto stsCP_iter = tsRecoSimCPMap.find(tsRef); + if (stsCP_iter != tsRecoSimCPMap.end()) { + const auto& stsCPassociated = stsCP_iter->val; + for (auto& sts : stsCPassociated) { + auto sts_id = (sts.first).get() - (edm::Ref(simTrackstersCP_h, 0)).get(); + trackstersCLUE3D_recoToSim_CP[i].push_back(sts_id); + trackstersCLUE3D_recoToSim_CP_score[i].push_back(sts.second.second); + trackstersCLUE3D_recoToSim_CP_sharedE[i].push_back(sts.second.first); + } + } + } + + // SimTracksters + nsimTrackstersCP = simTrackstersCP.size(); + trackstersCLUE3D_simToReco_CP.resize(nsimTrackstersCP); + trackstersCLUE3D_simToReco_CP_score.resize(nsimTrackstersCP); + trackstersCLUE3D_simToReco_CP_sharedE.resize(nsimTrackstersCP); + for (size_t i = 0; i < nsimTrackstersCP; ++i) { + const edm::Ref stsCPRef(simTrackstersCP_h, i); + + // STS-CP -> CLUE3D + const auto ts_iter = tsSimToRecoCPMap.find(stsCPRef); + if (ts_iter != tsSimToRecoCPMap.end()) { + const auto& tsAssociated = ts_iter->val; + for (auto& ts : tsAssociated) { + auto ts_idx = (ts.first).get() - (edm::Ref(tracksters_handle, 0)).get(); + trackstersCLUE3D_simToReco_CP[i].push_back(ts_idx); + trackstersCLUE3D_simToReco_CP_score[i].push_back(ts.second.second); + trackstersCLUE3D_simToReco_CP_sharedE[i].push_back(ts.second.first); + } + } + } + + // Tackster reco->sim associations + MergeTracksters_recoToSim_SC.resize(trackstersmerged.size()); + MergeTracksters_recoToSim_SC_score.resize(trackstersmerged.size()); + MergeTracksters_recoToSim_SC_sharedE.resize(trackstersmerged.size()); + for (size_t i = 0; i < trackstersmerged.size(); ++i) { + const edm::Ref tsRef(tracksters_merged_h, i); + + // CLUE3D -> STS-SC + const auto stsSC_iter = MergetsRecoSimSCMap.find(tsRef); + if (stsSC_iter != MergetsRecoSimSCMap.end()) { + const auto& stsSCassociated = stsSC_iter->val; + for (auto& sts : stsSCassociated) { + auto sts_id = (sts.first).get() - (edm::Ref(simTrackstersSC_h, 0)).get(); + MergeTracksters_recoToSim_SC[i].push_back(sts_id); + MergeTracksters_recoToSim_SC_score[i].push_back(sts.second.second); + MergeTracksters_recoToSim_SC_sharedE[i].push_back(sts.second.first); + } + } + } + + // SimTracksters + nsimTrackstersSC = simTrackstersSC.size(); + MergeTracksters_simToReco_SC.resize(nsimTrackstersSC); + MergeTracksters_simToReco_SC_score.resize(nsimTrackstersSC); + MergeTracksters_simToReco_SC_sharedE.resize(nsimTrackstersSC); + for (size_t i = 0; i < nsimTrackstersSC; ++i) { + const edm::Ref stsSCRef(simTrackstersSC_h, i); + + // STS-SC -> CLUE3D + const auto ts_iter = MergetsSimToRecoSCMap.find(stsSCRef); + if (ts_iter != MergetsSimToRecoSCMap.end()) { + const auto& tsAssociated = ts_iter->val; + for (auto& ts : tsAssociated) { + auto ts_idx = (ts.first).get() - (edm::Ref(tracksters_merged_h, 0)).get(); + MergeTracksters_simToReco_SC[i].push_back(ts_idx); + MergeTracksters_simToReco_SC_score[i].push_back(ts.second.second); + MergeTracksters_simToReco_SC_sharedE[i].push_back(ts.second.first); + } + } + } + + // Tackster reco->sim associations + MergeTracksters_recoToSim_CP.resize(trackstersmerged.size()); + MergeTracksters_recoToSim_CP_score.resize(trackstersmerged.size()); + MergeTracksters_recoToSim_CP_sharedE.resize(trackstersmerged.size()); + for (size_t i = 0; i < trackstersmerged.size(); ++i) { + const edm::Ref tsRef(tracksters_merged_h, i); + + // CLUE3D -> STS-CP + const auto stsCP_iter = MergetsRecoSimCPMap.find(tsRef); + if (stsCP_iter != MergetsRecoSimCPMap.end()) { + const auto& stsCPassociated = stsCP_iter->val; + for (auto& sts : stsCPassociated) { + auto sts_id = (sts.first).get() - (edm::Ref(simTrackstersCP_h, 0)).get(); + MergeTracksters_recoToSim_CP[i].push_back(sts_id); + MergeTracksters_recoToSim_CP_score[i].push_back(sts.second.second); + MergeTracksters_recoToSim_CP_sharedE[i].push_back(sts.second.first); + } + } + } + + // Tackster reco->sim associations + MergeTracksters_recoToSim_PU.resize(trackstersmerged.size()); + MergeTracksters_recoToSim_PU_score.resize(trackstersmerged.size()); + MergeTracksters_recoToSim_PU_sharedE.resize(trackstersmerged.size()); + for (size_t i = 0; i < trackstersmerged.size(); ++i) { + const edm::Ref tsRef(tracksters_merged_h, i); + + // CLUE3D -> STS-PU + const auto stsPU_iter = MergetsRecoSimPUMap.find(tsRef); + if (stsPU_iter != MergetsRecoSimPUMap.end()) { + const auto& stsPUassociated = stsPU_iter->val; + for (auto& sts : stsPUassociated) { + auto sts_id = (sts.first).get() - (edm::Ref(simTrackstersPU_h, 0)).get(); + MergeTracksters_recoToSim_PU[i].push_back(sts_id); + MergeTracksters_recoToSim_PU_score[i].push_back(sts.second.second); + MergeTracksters_recoToSim_PU_sharedE[i].push_back(sts.second.first); + } + } + } + + // SimTracksters + nsimTrackstersCP = simTrackstersCP.size(); + MergeTracksters_simToReco_CP.resize(nsimTrackstersCP); + MergeTracksters_simToReco_CP_score.resize(nsimTrackstersCP); + MergeTracksters_simToReco_CP_sharedE.resize(nsimTrackstersCP); + for (size_t i = 0; i < nsimTrackstersCP; ++i) { + const edm::Ref stsCPRef(simTrackstersCP_h, i); + + // STS-CP -> TrackstersMerge + const auto ts_iter = MergetsSimToRecoCPMap.find(stsCPRef); + if (ts_iter != MergetsSimToRecoCPMap.end()) { + const auto& tsAssociated = ts_iter->val; + for (auto& ts : tsAssociated) { + auto ts_idx = (ts.first).get() - (edm::Ref(tracksters_merged_h, 0)).get(); + MergeTracksters_simToReco_CP[i].push_back(ts_idx); + MergeTracksters_simToReco_CP_score[i].push_back(ts.second.second); + MergeTracksters_simToReco_CP_sharedE[i].push_back(ts.second.first); + } + } + } + + // SimTracksters + auto nsimTrackstersPU = simTrackstersPU.size(); + MergeTracksters_simToReco_PU.resize(nsimTrackstersPU); + MergeTracksters_simToReco_PU_score.resize(nsimTrackstersPU); + MergeTracksters_simToReco_PU_sharedE.resize(nsimTrackstersPU); + for (size_t i = 0; i < nsimTrackstersPU; ++i) { + const edm::Ref stsPURef(simTrackstersPU_h, i); + + // STS-PU -> Tracksters Merge + const auto ts_iter = MergetsSimToRecoPUMap.find(stsPURef); + if (ts_iter != MergetsSimToRecoPUMap.end()) { + const auto& tsAssociated = ts_iter->val; + for (auto& ts : tsAssociated) { + auto ts_idx = (ts.first).get() - (edm::Ref(tracksters_merged_h, 0)).get(); + MergeTracksters_simToReco_PU[i].push_back(ts_idx); + MergeTracksters_simToReco_PU_score[i].push_back(ts.second.second); + MergeTracksters_simToReco_PU_sharedE[i].push_back(ts.second.first); + } + } + } + + //Tracks + for (size_t i = 0; i < tracks.size(); i++) { + const auto& track = tracks[i]; + reco::TrackRef trackref = reco::TrackRef(tracks_h, i); + int iSide = int(track.eta() > 0); + const auto& fts = trajectoryStateTransform::outerFreeState((track), bFieldProd); + // to the HGCal front + const auto& tsos = prop.propagate(fts, firstDisk_[iSide]->surface()); + if (tsos.isValid()) { + const auto& globalPos = tsos.globalPosition(); + const auto& globalMom = tsos.globalMomentum(); + track_id.push_back(i); + track_hgcal_x.push_back(globalPos.x()); + track_hgcal_y.push_back(globalPos.y()); + track_hgcal_z.push_back(globalPos.z()); + track_hgcal_eta.push_back(globalPos.eta()); + track_hgcal_phi.push_back(globalPos.phi()); + track_hgcal_px.push_back(globalMom.x()); + track_hgcal_py.push_back(globalMom.y()); + track_hgcal_pz.push_back(globalMom.z()); + track_hgcal_pt.push_back(globalMom.perp()); + track_pt.push_back(track.pt()); + track_quality.push_back(track.quality(reco::TrackBase::highPurity)); + track_missing_outer_hits.push_back(track.missingOuterHits()); + track_charge.push_back(track.charge()); + track_time.push_back(trackTime[trackref]); + track_time_quality.push_back(trackTimeQual[trackref]); + track_time_err.push_back(trackTimeErr[trackref]); + track_nhits.push_back(tracks[i].recHitsSize()); + } + } + + if (saveCLUE3DTracksters_) + trackster_tree_->Fill(); + if (saveLCs_) + cluster_tree_->Fill(); + if (saveTICLCandidate_) + candidate_tree_->Fill(); + if (saveTrackstersMerged_) + tracksters_merged_tree_->Fill(); + if (saveAssociations_) + associations_tree_->Fill(); + if (saveSimTrackstersSC_) + simtrackstersSC_tree_->Fill(); + if (saveSimTrackstersCP_) + simtrackstersCP_tree_->Fill(); + if (saveTracks_) + tracks_tree_->Fill(); + if (saveSimTICLCandidate_) + simTICLCandidate_tree->Fill(); +} + +void TICLDumper::endJob() {} + +void TICLDumper::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("trackstersclue3d", edm::InputTag("ticlTrackstersCLUE3DHigh")); + desc.add("layerClusters", edm::InputTag("hgcalMergeLayerClusters")); + desc.add("layer_clustersTime", edm::InputTag("hgcalMergeLayerClusters", "timeLayerCluster")); + desc.add("ticlcandidates", edm::InputTag("ticlTrackstersMerge")); + 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("trackstersmerged", edm::InputTag("ticlTrackstersMerge")); + desc.add("simtrackstersSC", edm::InputTag("ticlSimTracksters")); + desc.add("simtrackstersCP", edm::InputTag("ticlSimTracksters", "fromCPs")); + desc.add("simtrackstersPU", edm::InputTag("ticlSimTracksters", "PU")); + desc.add("simTICLCandidates", edm::InputTag("ticlSimTracksters")); + desc.add("recoToSimAssociatorSC", + edm::InputTag("tracksterSimTracksterAssociationPRbyCLUE3D", "recoToSim")); + desc.add("simToRecoAssociatorSC", + edm::InputTag("tracksterSimTracksterAssociationPRbyCLUE3D", "simToReco")); + desc.add("recoToSimAssociatorCP", + edm::InputTag("tracksterSimTracksterAssociationLinkingbyCLUE3D", "recoToSim")); + desc.add("simToRecoAssociatorCP", + edm::InputTag("tracksterSimTracksterAssociationLinkingbyCLUE3D", "simToReco")); + desc.add("MergerecoToSimAssociatorSC", + edm::InputTag("tracksterSimTracksterAssociationPR", "recoToSim")); + desc.add("MergesimToRecoAssociatorSC", + edm::InputTag("tracksterSimTracksterAssociationPR", "simToReco")); + desc.add("MergerecoToSimAssociatorCP", + edm::InputTag("tracksterSimTracksterAssociationLinking", "recoToSim")); + desc.add("MergesimToRecoAssociatorCP", + edm::InputTag("tracksterSimTracksterAssociationLinking", "simToReco")); + desc.add("MergerecoToSimAssociatorPU", + edm::InputTag("tracksterSimTracksterAssociationLinkingPU", "recoToSim")); + desc.add("MergesimToRecoAssociatorPU", + edm::InputTag("tracksterSimTracksterAssociationLinkingPU", "simToReco")); + desc.add("simclusters", edm::InputTag("mix", "MergedCaloTruth")); + desc.add("caloparticles", edm::InputTag("mix", "MergedCaloTruth")); + desc.add("detector", "HGCAL"); + desc.add("propagator", "PropagatorWithMaterial"); + + desc.add("saveLCs", true); + desc.add("saveCLUE3DTracksters", true); + desc.add("saveTrackstersMerged", true); + desc.add("saveSimTrackstersSC", true); + desc.add("saveSimTrackstersCP", true); + desc.add("saveTICLCandidate", true); + desc.add("saveSimTICLCandidate", true); + desc.add("saveTracks", true); + desc.add("saveAssociations", true); + descriptions.add("ticlDumper", desc); +} + +DEFINE_FWK_MODULE(TICLDumper); diff --git a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc index ddc6484ec3eed..932a829beece5 100644 --- a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc @@ -280,9 +280,9 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es 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; +#ifdef EDM_ML_DEBUG LogDebug("TrackstersMergeProducer") << "PDG ID " << cand.pdgId() << " charge " << cand.charge() << " p " << cand.p() << std::endl; LogDebug("TrackstersMergeProducer") << "track id (p) : " << track_idx << " (" @@ -292,6 +292,7 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es // Merge included tracksters ticl::Trackster outTrackster; + outTrackster.setTrackIdx(track_idx); auto updated_size = 0; for (const auto &ts_ptr : trackster_ptrs) { #ifdef EDM_ML_DEBUG @@ -335,8 +336,6 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es } outTrackster.zeroProbabilities(); - if (!track_ptr.isNull()) - outTrackster.setSeed(track_h.id(), track_ptr.get() - (edm::Ptr(track_h, 0)).get()); if (!outTrackster.vertices().empty()) { resultTrackstersMerged->push_back(outTrackster); } @@ -534,19 +533,17 @@ void TrackstersMergeProducer::assignTimeToCandidates(std::vector 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; + float invTimeErr = 0.f; for (const auto &tr : cand.tracksters()) { if (tr->timeError() > 0) { auto invTimeESq = pow(tr->timeError(), -2); time += tr->time() * invTimeESq; - timeErr += invTimeESq; + invTimeErr += invTimeESq; } } - if (timeErr > 0) { - timeErr = 1. / timeErr; - - cand.setTime(time * timeErr); - cand.setTimeError(sqrt(timeErr)); + if (invTimeErr > 0) { + cand.setTime(time / invTimeErr); + cand.setTimeError(sqrt(1.f / invTimeErr)); } } } diff --git a/RecoHGCal/TICL/python/customiseTICLFromReco.py b/RecoHGCal/TICL/python/customiseTICLFromReco.py index 89d34942429c6..60f2fa260cf6e 100644 --- a/RecoHGCal/TICL/python/customiseTICLFromReco.py +++ b/RecoHGCal/TICL/python/customiseTICLFromReco.py @@ -2,6 +2,7 @@ from RecoHGCal.TICL.iterativeTICL_cff import * from RecoLocalCalo.HGCalRecProducers.hgcalLayerClusters_cff import hgcalLayerClustersEE, hgcalLayerClustersHSi, hgcalLayerClustersHSci from RecoLocalCalo.HGCalRecProducers.hgcalMergeLayerClusters_cfi import hgcalMergeLayerClusters +from RecoHGCal.TICL.ticlDumper_cfi import ticlDumper # Validation from Validation.HGCalValidation.HGCalValidator_cfi import * from RecoLocalCalo.HGCalRecProducers.hgcalRecHitMapProducer_cfi import hgcalRecHitMapProducer @@ -11,10 +12,12 @@ # Automatic addition of the customisation function from RecoHGCal.Configuration.RecoHGCal_EventContent_cff from RecoHGCal.Configuration.RecoHGCal_EventContent_cff import customiseHGCalOnlyEventContent +from SimCalorimetry.HGCalAssociatorProducers.simTracksterAssociatorByEnergyScore_cfi import simTracksterAssociatorByEnergyScore as simTsAssocByEnergyScoreProducer +from SimCalorimetry.HGCalAssociatorProducers.TSToSimTSAssociation_cfi import tracksterSimTracksterAssociationLinking, tracksterSimTracksterAssociationPR, tracksterSimTracksterAssociationLinkingbyCLUE3D, tracksterSimTracksterAssociationPRbyCLUE3D, tracksterSimTracksterAssociationLinkingPU, tracksterSimTracksterAssociationPRPU def customiseTICLFromReco(process): -# TensorFlow ESSource + # TensorFlow ESSource process.TFESSource = cms.Task(process.trackdnn_source) process.hgcalLayerClustersTask = cms.Task(process.hgcalLayerClustersEE, @@ -34,11 +37,20 @@ def customiseTICLFromReco(process): process.layerClusterCaloParticleAssociationProducer, process.scAssocByEnergyScoreProducer, process.layerClusterSimClusterAssociationProducer, - ) + process.simTsAssocByEnergyScoreProducer, + process.simTracksterHitLCAssociatorByEnergyScoreProducer, + process.tracksterSimTracksterAssociationLinking, + process.tracksterSimTracksterAssociationPR, + process.tracksterSimTracksterAssociationLinkingbyCLUE3D, + process.tracksterSimTracksterAssociationPRbyCLUE3D, + process.tracksterSimTracksterAssociationLinkingPU, + process.tracksterSimTracksterAssociationPRPU + ) + process.TICL_Validator = cms.Task(process.hgcalValidator) process.TICL_Validation = cms.Path(process.TICL_ValidationProducers, process.TICL_Validator - ) + ) # Path and EndPath definitions process.FEVTDEBUGHLToutput_step = cms.EndPath(process.FEVTDEBUGHLToutput) process.DQMoutput_step = cms.EndPath(process.DQMoutput) @@ -48,7 +60,28 @@ def customiseTICLFromReco(process): process.TICL_Validation, process.FEVTDEBUGHLToutput_step, process.DQMoutput_step) -#call to customisation function customiseHGCalOnlyEventContent imported from RecoHGCal.Configuration.RecoHGCal_EventContent_cff +# call to customisation function customiseHGCalOnlyEventContent imported from RecoHGCal.Configuration.RecoHGCal_EventContent_cff process = customiseHGCalOnlyEventContent(process) return process + + +def customiseTICLForDumper(process): + + process.ticlDumper = ticlDumper.clone( + saveLCs=True, + saveCLUE3DTracksters=True, + saveTrackstersMerged=True, + saveSimTrackstersSC=True, + saveSimTrackstersCP=True, + saveTICLCandidate=True, + saveSimTICLCandidate=True, + saveTracks=True, + saveAssociations=True, + ) + process.TFileService = cms.Service("TFileService", + fileName=cms.string("histo.root") + ) + process.FEVTDEBUGHLToutput_step = cms.EndPath( + process.FEVTDEBUGHLToutput + process.ticlDumper) + return process diff --git a/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSAssociatorByEnergyScoreImpl.cc b/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSAssociatorByEnergyScoreImpl.cc index bdb3ec90dc582..6494153d27438 100644 --- a/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSAssociatorByEnergyScoreImpl.cc +++ b/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSAssociatorByEnergyScoreImpl.cc @@ -128,7 +128,7 @@ hgcal::association TSToSimTSAssociatorByEnergyScoreImpl::makeConnections( tssInSimTrackster[st.clusterId].tracksterIdToEnergyAndScore[tsId].first += lcFractionInTs * st.fraction * layerClusters[lcId].energy(); //TS_i -> ST_j, ST_k, ... - stsInTrackster[tsId].emplace_back(st.clusterId, 0.f); + stsInTrackster[tsId].emplace_back(st.clusterId, std::make_pair(0.f, 0.f)); } } } // End loop over LayerClusters in Trackster @@ -294,9 +294,10 @@ hgcal::association TSToSimTSAssociatorByEnergyScoreImpl::makeConnections( // SimTrackster, assigned score 1 if (tracksters[tsId].raw_energy() == 0. && !stsInTrackster[tsId].empty()) { for (auto& stPair : stsInTrackster[tsId]) { - stPair.second = 1.; + stPair.second.second = 1.; LogDebug("TSToSimTSAssociatorByEnergyScoreImpl") - << "TracksterId:\t " << tsId << "\tST id:\t" << stPair.first << "\tscore\t " << stPair.second << "\n"; + << "TracksterId:\t " << tsId << "\tST id:\t" << stPair.first << "\tenergy" << stPair.second.first + << "\tscore\t " << stPair.second.second << "\n"; } continue; } @@ -332,13 +333,13 @@ hgcal::association TSToSimTSAssociatorByEnergyScoreImpl::makeConnections( stFraction = findLCIt->fraction; } } - stPair.second += + stPair.second.second += (lcFractionInTs - stFraction) * (lcFractionInTs - stFraction) * lcEnergyWeight * invTracksterEnergyWeight; #ifdef EDM_ML_DEBUG LogDebug("TSToSimTSAssociatorByEnergyScoreImpl") << "lcId:\t" << (uint32_t)lcId << "\ttracksterId:\t" << tsId << "\ttsFraction,stFraction:\t" << lcFractionInTs << ", " << stFraction << "\tlcEnergyWeight:\t" << lcEnergyWeight << "\tcurrent score:\t" - << stPair.second << "\tinvTracksterEnergyWeight:\t" << invTracksterEnergyWeight << "\n"; + << stPair.second.second << "\tinvTracksterEnergyWeight:\t" << invTracksterEnergyWeight << "\n"; #endif } } // End of loop over LayerClusters in Trackster @@ -446,9 +447,10 @@ hgcal::RecoToSimCollectionSimTracksters TSToSimTSAssociatorByEnergyScoreImpl::as LogDebug("TSToSimTSAssociatorByEnergyScoreImpl") << "Trackster Id:\t" << tsId << "\tSimTrackster id:\t" << stPair.first << "\tscore:\t" << stPair.second << "\n"; // Fill AssociationMap - returnValue.insert(edm::Ref(tCH, tsId), // Ref to TS - std::make_pair(edm::Ref(sTCH, stPair.first), - stPair.second) // Pair + returnValue.insert( + edm::Ref(tCH, tsId), // Ref to TS + std::make_pair(edm::Ref(sTCH, stPair.first), //Pair + std::make_pair(stPair.second.first, stPair.second.second)) // Pair ); } } diff --git a/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSAssociatorByEnergyScoreImpl.h b/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSAssociatorByEnergyScoreImpl.h index 0de2bb23bd1da..067a5c544c16b 100644 --- a/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSAssociatorByEnergyScoreImpl.h +++ b/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSAssociatorByEnergyScoreImpl.h @@ -56,7 +56,7 @@ namespace hgcal { // the SimTracksters (via their ids (stIds)) that share at least one LayerCluster. In that pair // it stores the score (tsId->(stId,score)). Keep in mind that the association is not unique, since there could be // several instances of the same SimTrackster from several related SimClusters that each contributed to the same Trackster. - typedef std::vector>> tracksterToSimTrackster; + typedef std::vector>>> tracksterToSimTrackster; // This is used to save the simTracksterOnLayer structure for all simTracksters. // It is not exactly what is returned outside, but out of its entries, the output object is build. typedef std::vector simTracksterToTrackster; diff --git a/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSAssociatorEDProducer.cc b/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSAssociatorEDProducer.cc index 26bdb30f95e4f..53a9ed1e5998f 100644 --- a/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSAssociatorEDProducer.cc +++ b/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSAssociatorEDProducer.cc @@ -44,8 +44,8 @@ class TSToSimTSAssociatorEDProducer : public edm::global::EDProducer<> { }; TSToSimTSAssociatorEDProducer::TSToSimTSAssociatorEDProducer(const edm::ParameterSet &pset) { - produces(); - produces(); + produces("simToReco"); + produces("recoToSim"); TSCollectionToken_ = consumes(pset.getParameter("label_tst")); SimTSCollectionToken_ = consumes(pset.getParameter("label_simTst")); @@ -87,8 +87,8 @@ void TSToSimTSAssociatorEDProducer::produce(edm::StreamID, edm::Event &iEvent, c auto rts = std::make_unique(recSimColl); auto str = std::make_unique(simRecColl); - iEvent.put(std::move(rts)); - iEvent.put(std::move(str)); + iEvent.put(std::move(rts), "recoToSim"); + iEvent.put(std::move(str), "simToReco"); } // define this as a plug-in diff --git a/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorByEnergyScoreImpl.cc b/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorByEnergyScoreImpl.cc new file mode 100644 index 0000000000000..d3bac4802178d --- /dev/null +++ b/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorByEnergyScoreImpl.cc @@ -0,0 +1,259 @@ +#include + +#include "TSToSimTSHitLCAssociatorByEnergyScoreImpl.h" +#include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +TSToSimTSHitLCAssociatorByEnergyScoreImpl::TSToSimTSHitLCAssociatorByEnergyScoreImpl( + edm::EDProductGetter const& productGetter, + bool hardScatterOnly, + std::shared_ptr recHitTools, + const std::unordered_map* hitMap) + : hardScatterOnly_(hardScatterOnly), recHitTools_(recHitTools), hitMap_(hitMap), productGetter_(&productGetter) { + layers_ = recHitTools_->lastLayerBH(); +} + +hgcal::association_t TSToSimTSHitLCAssociatorByEnergyScoreImpl::makeConnections( + const edm::Handle& tCH, + const edm::Handle& lCCH, + const edm::Handle& sCCH, + const edm::Handle& cPCH, + const edm::Handle& sTCH) const { + // Get collections + const auto& tracksters = *tCH.product(); + const auto& layerClusters = *lCCH.product(); + const auto& sC = *sCCH.product(); + const auto& cP = *cPCH.product(); + const auto cPHandle_id = cPCH.id(); + const auto& simTSs = *sTCH.product(); + const auto nTracksters = tracksters.size(); + const auto nSimTracksters = simTSs.size(); + + std::unordered_map>> detIdSimTSId_Map; + std::unordered_map>> detIdSimClusterId_Map; + std::unordered_map>> detIdCaloParticleId_Map; + std::unordered_map>> detIdToRecoTSId_Map; + + hgcal::sharedEnergyAndScore_t recoToSim_sharedEnergyAndScore; + hgcal::sharedEnergyAndScore_t simToReco_sharedEnergyAndScore; + + recoToSim_sharedEnergyAndScore.resize(nTracksters); + simToReco_sharedEnergyAndScore.resize(nSimTracksters); + + for (size_t i = 0; i < nTracksters; ++i) + recoToSim_sharedEnergyAndScore[i].resize(nSimTracksters); + for (size_t i = 0; i < nSimTracksters; ++i) + simToReco_sharedEnergyAndScore[i].resize(nTracksters); + + // fill sim maps + + for (size_t i = 0; i < sC.size(); ++i) { + for (const auto& haf : sC[i].hits_and_fractions()) { + detIdSimClusterId_Map[haf.first].emplace_back(i, haf.second); + } + } + + for (size_t i = 0; i < cP.size(); ++i) { + for (const auto& sc : cP[i].simClusters()) { + for (const auto& haf : sc->hits_and_fractions()) { + auto hitId = haf.first; + auto found = std::find_if(detIdCaloParticleId_Map[hitId].begin(), + detIdCaloParticleId_Map[hitId].end(), + [=](const std::pair& v) { return v.first == static_cast(i); }); + if (found == detIdCaloParticleId_Map[hitId].end()) + detIdCaloParticleId_Map[haf.first].emplace_back(i, haf.second); + else + found->second += haf.second; + } + } + } + + for (size_t i = 0; i < nSimTracksters; ++i) { + const auto& lcsInSimTrackster = simTSs[i].vertices(); + const auto& multiplicities = simTSs[i].vertex_multiplicity(); + for (size_t j = 0; j < lcsInSimTrackster.size(); ++j) { + const auto& v = lcsInSimTrackster[j]; + float fraction = 1.f / multiplicities[j]; + for (const auto& haf : layerClusters[v].hitsAndFractions()) { + detIdSimTSId_Map[haf.first].emplace_back(i, haf.second * fraction); + } + } + } + + // fill reco map + + for (size_t i = 0; i < nTracksters; ++i) { + const auto& lcsInSimTrackster = tracksters[i].vertices(); + const auto& multiplicities = tracksters[i].vertex_multiplicity(); + for (size_t j = 0; j < lcsInSimTrackster.size(); ++j) { + const auto& v = lcsInSimTrackster[j]; + float fraction = 1.f / multiplicities[j]; + for (const auto& haf : layerClusters[v].hitsAndFractions()) { + detIdToRecoTSId_Map[haf.first].emplace_back(i, haf.second * fraction); + } + } + } + + std::vector denominator_simToReco(nSimTracksters, 0.f); + std::vector> numerator_simToReco(nSimTracksters); + std::vector> sharedEnergy(nSimTracksters); + + for (size_t i = 0; i < nSimTracksters; ++i) { + numerator_simToReco[i].resize(nTracksters, 0.f); + sharedEnergy[i].resize(nTracksters, 0.f); + + const auto seedIndex = simTSs[i].seedIndex(); + const auto& lcsInSimTrackster = simTSs[i].vertices(); + + for (const auto& v : lcsInSimTrackster) { + for (const auto& haf : layerClusters[v].hitsAndFractions()) { + const auto hitId = haf.first; + float simFraction = 0.f; + + std::vector>::iterator found; + if (simTSs[i].seedID() == cPHandle_id) { + found = std::find_if(detIdSimTSId_Map[hitId].begin(), + detIdSimTSId_Map[hitId].end(), + [=](const std::pair& v) { return v.first == seedIndex; }); + if (found != detIdSimTSId_Map[hitId].end()) { + const auto iLC = std::find(simTSs[i].vertices().begin(), simTSs[i].vertices().end(), v); + const auto lcFraction = + 1.f / simTSs[i].vertex_multiplicity(std::distance(std::begin(simTSs[i].vertices()), iLC)); + simFraction = found->second * lcFraction; + } + } else { + found = std::find_if(detIdSimClusterId_Map[hitId].begin(), + detIdSimClusterId_Map[hitId].end(), + [=](const std::pair& v) { return v.first == seedIndex; }); + if (found != detIdSimClusterId_Map[hitId].end()) { + simFraction = found->second; + } + } + + float hitEnergy = hitMap_->find(hitId)->second->energy(); + float hitEnergySquared = hitEnergy * hitEnergy; + float simFractionSquared = simFraction * simFraction; + denominator_simToReco[i] += simFractionSquared * hitEnergySquared; + for (size_t j = 0; j < nTracksters; ++j) { + float recoFraction = 0.f; + + auto found_reco = + std::find_if(detIdToRecoTSId_Map[hitId].begin(), + detIdToRecoTSId_Map[hitId].end(), + [=](const std::pair& v) { return v.first == static_cast(j); }); + if (found_reco != detIdToRecoTSId_Map[hitId].end()) + recoFraction = found_reco->second; + numerator_simToReco[i][j] += + std::min(simFractionSquared, (simFraction - recoFraction) * (simFraction - recoFraction)) * + hitEnergySquared; + sharedEnergy[i][j] += std::min(simFraction, recoFraction) * hitEnergy; + } + } + } + } + + std::vector denominator_recoToSim(nTracksters, 0.f); + std::vector> numerator_recoToSim(nTracksters); + + for (unsigned int i = 0; i < nTracksters; ++i) { + numerator_recoToSim[i].resize(nSimTracksters, 0.f); + const auto& lcsInTrackster = tracksters[i].vertices(); + for (const auto& v : lcsInTrackster) { + for (const auto& haf : layerClusters[v].hitsAndFractions()) { + const auto hitId = haf.first; + float recoFraction = 0.f; + + auto found = std::find_if(detIdToRecoTSId_Map[hitId].begin(), + detIdToRecoTSId_Map[hitId].end(), + [=](const std::pair& v) { return v.first == static_cast(i); }); + if (found != detIdToRecoTSId_Map[hitId].end()) + recoFraction = found->second; + + float hitEnergy = hitMap_->find(hitId)->second->energy(); + float hitEnergySquared = hitEnergy * hitEnergy; + float recoFractionSquared = recoFraction * recoFraction; + denominator_recoToSim[i] += recoFractionSquared * hitEnergySquared; + + for (size_t j = 0; j < nSimTracksters; ++j) { + float simFraction = 0.f; + + auto found_sim = std::find_if(detIdSimTSId_Map[hitId].begin(), + detIdSimTSId_Map[hitId].end(), + [=](const std::pair& v) { return v.first == static_cast(j); }); + if (found_sim != detIdSimTSId_Map[hitId].end()) + simFraction = found_sim->second; + numerator_recoToSim[i][j] += + std::min(recoFractionSquared, (simFraction - recoFraction) * (simFraction - recoFraction)) * + hitEnergySquared; + } + } + } + } + + // compute score + + for (size_t i = 0; i < nSimTracksters; ++i) { + for (size_t j = 0; j < nTracksters; ++j) { + simToReco_sharedEnergyAndScore[i][j].first = sharedEnergy[i][j]; + simToReco_sharedEnergyAndScore[i][j].second = numerator_simToReco[i][j] / denominator_simToReco[i]; + recoToSim_sharedEnergyAndScore[j][i].first = sharedEnergy[i][j]; + recoToSim_sharedEnergyAndScore[j][i].second = numerator_recoToSim[j][i] / denominator_recoToSim[j]; + } + } + + return {recoToSim_sharedEnergyAndScore, simToReco_sharedEnergyAndScore}; +} + +hgcal::RecoToSimCollectionSimTracksters TSToSimTSHitLCAssociatorByEnergyScoreImpl::associateRecoToSim( + const edm::Handle& tCH, + const edm::Handle& lCCH, + const edm::Handle& sCCH, + const edm::Handle& cPCH, + const edm::Handle& sTCH) const { + hgcal::RecoToSimCollectionSimTracksters returnValue(productGetter_); + const auto& links = makeConnections(tCH, lCCH, sCCH, cPCH, sTCH); + const auto& recoToSim_sharedEnergyAndScore = std::get<0>(links); + for (std::size_t tsId = 0; tsId < recoToSim_sharedEnergyAndScore.size(); ++tsId) { + std::size_t numSimTracksters = recoToSim_sharedEnergyAndScore[tsId].size(); + for (std::size_t simTsId = 0; simTsId < numSimTracksters; ++simTsId) { + LogDebug("TSToSimTSHitLCAssociatorByEnergyScoreImpl") + << " Trackster Id:\t" << tsId << "\tSimTrackster id:\t" << recoToSim_sharedEnergyAndScore[tsId][simTsId].first + << "\tscore:\t" << recoToSim_sharedEnergyAndScore[tsId][simTsId].second << "\n"; + // Fill AssociationMap + returnValue.insert( + edm::Ref(tCH, tsId), // Ref to TS + std::make_pair( + edm::Ref(sTCH, simTsId), + std::make_pair(recoToSim_sharedEnergyAndScore[tsId][simTsId].first, + recoToSim_sharedEnergyAndScore[tsId][simTsId].second)) // Pair + ); + } + } + return returnValue; +} + +hgcal::SimToRecoCollectionSimTracksters TSToSimTSHitLCAssociatorByEnergyScoreImpl::associateSimToReco( + const edm::Handle& tCH, + const edm::Handle& lCCH, + const edm::Handle& sCCH, + const edm::Handle& cPCH, + const edm::Handle& sTCH) const { + hgcal::SimToRecoCollectionSimTracksters returnValue(productGetter_); + const auto& links = makeConnections(tCH, lCCH, sCCH, cPCH, sTCH); + const auto& simToReco_sharedEnergyAndScore = std::get<1>(links); + for (std::size_t simTsId = 0; simTsId < simToReco_sharedEnergyAndScore.size(); ++simTsId) { + std::size_t numTracksters = simToReco_sharedEnergyAndScore[simTsId].size(); + for (std::size_t tsId = 0; tsId < numTracksters; ++tsId) { + LogDebug("TSToSimTSHitLCAssociatorByEnergyScoreImpl") + << "Trackster Id:\t" << tsId << "\tSimTrackster id:\t" << simTsId << " Shared energy " + << simToReco_sharedEnergyAndScore[simTsId][tsId].first << "\tscore:\t" + << simToReco_sharedEnergyAndScore[simTsId][tsId].second << "\n"; + // Fill AssociationMap + returnValue.insert(edm::Ref(sTCH, simTsId), + std::make_pair(edm::Ref(tCH, tsId), + std::make_pair(simToReco_sharedEnergyAndScore[simTsId][tsId].first, + simToReco_sharedEnergyAndScore[simTsId][tsId].second))); + } + } + return returnValue; +} diff --git a/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorByEnergyScoreImpl.h b/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorByEnergyScoreImpl.h new file mode 100644 index 0000000000000..ff6d95d7f2b8d --- /dev/null +++ b/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorByEnergyScoreImpl.h @@ -0,0 +1,78 @@ +// Original Author: Leonardo Cristella + +#include +#include +#include +#include // shared_ptr + +#include "DataFormats/ForwardDetId/interface/HGCalDetId.h" +#include "DataFormats/HGCRecHit/interface/HGCRecHit.h" +#include "SimDataFormats/Associations/interface/TracksterToSimTracksterHitLCAssociator.h" +#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h" + +namespace edm { + class EDProductGetter; +} + +namespace hgcal { + + struct detIdInfoInCluster { + bool operator==(const detIdInfoInCluster &o) const { return clusterId == o.clusterId; }; + long unsigned int clusterId; + float fraction; + }; + + struct detIdInfoInTrackster { + bool operator==(const detIdInfoInTrackster &o) const { return tracksterId == o.tracksterId; }; + unsigned int tracksterId; + long unsigned int clusterId; + float fraction; + }; + + struct caloParticleOnLayer { + unsigned int caloParticleId; + float energy = 0; + std::vector> hits_and_fractions; + std::unordered_map> layerClusterIdToEnergyAndScore; + }; + + // This object connects a Trackster, identified through its id (tsId), with a vector of pairs containing all + // the SimTracksters (via their ids (stIds)) that share at least one LayerCluster. In that pair + // it stores the score (tsId->(stId,score)). Keep in mind that the association is not unique, since there could be + // several instances of the same SimTrackster from several related SimClusters that each contributed to the same Trackster. +} // namespace hgcal + +class TSToSimTSHitLCAssociatorByEnergyScoreImpl : public hgcal::TracksterToSimTracksterHitLCAssociatorBaseImpl { +public: + explicit TSToSimTSHitLCAssociatorByEnergyScoreImpl(edm::EDProductGetter const &, + bool, + std::shared_ptr, + const std::unordered_map *); + + hgcal::association_t makeConnections(const edm::Handle &tCH, + const edm::Handle &lCCH, + const edm::Handle &sCCH, + const edm::Handle &cPCH, + const edm::Handle &sTCH) const; + + hgcal::RecoToSimCollectionSimTracksters associateRecoToSim( + const edm::Handle &tCH, + const edm::Handle &lCCH, + const edm::Handle &sCCH, + const edm::Handle &cPCH, + const edm::Handle &sTCH) const override; + + hgcal::SimToRecoCollectionSimTracksters associateSimToReco( + const edm::Handle &tCH, + const edm::Handle &lCCH, + const edm::Handle &sCCH, + const edm::Handle &cPCH, + const edm::Handle &sTCH) const override; + +private: + const bool hardScatterOnly_; + std::shared_ptr recHitTools_; + const std::unordered_map *hitMap_; + unsigned layers_; + edm::EDProductGetter const *productGetter_; +}; diff --git a/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorByEnergyScoreProducer.cc b/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorByEnergyScoreProducer.cc new file mode 100644 index 0000000000000..5c2f66d9c1d0e --- /dev/null +++ b/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorByEnergyScoreProducer.cc @@ -0,0 +1,67 @@ +// Original author: Leonardo Cristella + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/global/EDProducer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/EDGetToken.h" +#include "FWCore/Utilities/interface/ESGetToken.h" + +#include "SimDataFormats/Associations/interface/TracksterToSimTracksterHitLCAssociator.h" +#include "TSToSimTSHitLCAssociatorByEnergyScoreImpl.h" + +class TSToSimTSHitLCAssociatorByEnergyScoreProducer : public edm::global::EDProducer<> { +public: + explicit TSToSimTSHitLCAssociatorByEnergyScoreProducer(const edm::ParameterSet &); + ~TSToSimTSHitLCAssociatorByEnergyScoreProducer() override; + + static void fillDescriptions(edm::ConfigurationDescriptions &descriptions); + +private: + void produce(edm::StreamID, edm::Event &, const edm::EventSetup &) const override; + edm::EDGetTokenT> hitMap_; + edm::ESGetToken caloGeometry_; + const bool hardScatterOnly_; + std::shared_ptr rhtools_; +}; + +TSToSimTSHitLCAssociatorByEnergyScoreProducer::TSToSimTSHitLCAssociatorByEnergyScoreProducer(const edm::ParameterSet &ps) + : hitMap_(consumes>(ps.getParameter("hitMapTag"))), + caloGeometry_(esConsumes()), + hardScatterOnly_(ps.getParameter("hardScatterOnly")) { + rhtools_.reset(new hgcal::RecHitTools()); + + // Register the product + produces(); +} + +TSToSimTSHitLCAssociatorByEnergyScoreProducer::~TSToSimTSHitLCAssociatorByEnergyScoreProducer() {} + +void TSToSimTSHitLCAssociatorByEnergyScoreProducer::produce(edm::StreamID, + edm::Event &iEvent, + const edm::EventSetup &es) const { + edm::ESHandle geom = es.getHandle(caloGeometry_); + rhtools_->setGeometry(*geom); + + const auto hitMap = &iEvent.get(hitMap_); + + auto impl = std::make_unique( + iEvent.productGetter(), hardScatterOnly_, rhtools_, hitMap); + auto toPut = std::make_unique(std::move(impl)); + iEvent.put(std::move(toPut)); +} + +void TSToSimTSHitLCAssociatorByEnergyScoreProducer::fillDescriptions(edm::ConfigurationDescriptions &cfg) { + edm::ParameterSetDescription desc; + desc.add("hitMapTag", edm::InputTag("hgcalRecHitMapProducer")); + desc.add("hardScatterOnly", true); + + cfg.add("simTracksterHitLCAssociatorByEnergyScore", desc); +} + +//define this as a plug-in +DEFINE_FWK_MODULE(TSToSimTSHitLCAssociatorByEnergyScoreProducer); diff --git a/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorEDProducer.cc b/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorEDProducer.cc new file mode 100644 index 0000000000000..e66c29a89a91b --- /dev/null +++ b/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorEDProducer.cc @@ -0,0 +1,103 @@ +// +// Original Author: Leonardo Cristella +// Created: Wed Mar 30 10:52:11 CET 2022 +// +// + +// system include files +#include +#include + +// user include files +#include "FWCore/Framework/interface/global/EDProducer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/Framework/interface/ESHandle.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "SimDataFormats/Associations/interface/TracksterToSimTracksterHitLCAssociator.h" + +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "DataFormats/HGCalReco/interface/Trackster.h" +#include "SimDataFormats/CaloAnalysis/interface/CaloParticleFwd.h" + +#include "FWCore/Utilities/interface/EDGetToken.h" + +class TSToSimTSHitLCAssociatorEDProducer : public edm::global::EDProducer<> { +public: + explicit TSToSimTSHitLCAssociatorEDProducer(const edm::ParameterSet &); + ~TSToSimTSHitLCAssociatorEDProducer() override; + +private: + void produce(edm::StreamID, edm::Event &, const edm::EventSetup &) const override; + + edm::EDGetTokenT TSCollectionToken_; + edm::EDGetTokenT SimTSCollectionToken_; + edm::EDGetTokenT SimTSFromCPCollectionToken_; + edm::EDGetTokenT LCCollectionToken_; + edm::EDGetTokenT SCCollectionToken_; + edm::EDGetTokenT CPCollectionToken_; + edm::EDGetTokenT>> simTrackstersMap_; + hgcal::validationType valType_; + edm::EDGetTokenT associatorToken_; +}; + +TSToSimTSHitLCAssociatorEDProducer::TSToSimTSHitLCAssociatorEDProducer(const edm::ParameterSet &pset) { + produces("simToReco"); + produces("recoToSim"); + + TSCollectionToken_ = consumes(pset.getParameter("label_tst")); + SimTSCollectionToken_ = consumes(pset.getParameter("label_simTst")); + LCCollectionToken_ = consumes(pset.getParameter("label_lcl")); + SCCollectionToken_ = consumes(pset.getParameter("label_scl")); + CPCollectionToken_ = consumes(pset.getParameter("label_cp")); + associatorToken_ = + consumes(pset.getParameter("associator")); +} + +TSToSimTSHitLCAssociatorEDProducer::~TSToSimTSHitLCAssociatorEDProducer() {} + +void TSToSimTSHitLCAssociatorEDProducer::produce(edm::StreamID, + edm::Event &iEvent, + const edm::EventSetup &iSetup) const { + using namespace edm; + + edm::Handle theAssociator; + iEvent.getByToken(associatorToken_, theAssociator); + + Handle TSCollection; + iEvent.getByToken(TSCollectionToken_, TSCollection); + + Handle SimTSCollection; + iEvent.getByToken(SimTSCollectionToken_, SimTSCollection); + + Handle LCCollection; + iEvent.getByToken(LCCollectionToken_, LCCollection); + + Handle SCCollection; + iEvent.getByToken(SCCollectionToken_, SCCollection); + + Handle CPCollection; + iEvent.getByToken(CPCollectionToken_, CPCollection); + + // associate TS and SimTS + LogTrace("AssociatorValidator") << "Calling associateRecoToSim method\n"; + + hgcal::RecoToSimCollectionSimTracksters recSimColl = + theAssociator->associateRecoToSim(TSCollection, LCCollection, SCCollection, CPCollection, SimTSCollection); + + LogTrace("AssociatorValidator") << "Calling associateSimToReco method\n"; + hgcal::SimToRecoCollectionSimTracksters simRecColl = + theAssociator->associateSimToReco(TSCollection, LCCollection, SCCollection, CPCollection, SimTSCollection); + + auto rts = std::make_unique(recSimColl); + auto str = std::make_unique(simRecColl); + + iEvent.put(std::move(rts), "recoToSim"); + iEvent.put(std::move(str), "simToReco"); +} + +DEFINE_FWK_MODULE(TSToSimTSHitLCAssociatorEDProducer); diff --git a/SimCalorimetry/HGCalAssociatorProducers/python/TSToSimTSAssociation_cfi.py b/SimCalorimetry/HGCalAssociatorProducers/python/TSToSimTSAssociation_cfi.py index 6acf5bf94a386..1d012306ab86d 100644 --- a/SimCalorimetry/HGCalAssociatorProducers/python/TSToSimTSAssociation_cfi.py +++ b/SimCalorimetry/HGCalAssociatorProducers/python/TSToSimTSAssociation_cfi.py @@ -1,8 +1,58 @@ import FWCore.ParameterSet.Config as cms -tracksterSimTracksterAssociation = cms.EDProducer("TSToSimTSAssociatorEDProducer", - associator = cms.InputTag('simTsAssocByEnergyScoreProducer'), +tracksterSimTracksterAssociationLinking = cms.EDProducer("TSToSimTSHitLCAssociatorEDProducer", + associator = cms.InputTag('simTracksterHitLCAssociatorByEnergyScoreProducer'), label_tst = cms.InputTag("ticlTrackstersMerge"), + label_simTst = cms.InputTag("ticlSimTracksters", "fromCPs"), + label_lcl = cms.InputTag("hgcalMergeLayerClusters"), + label_scl = cms.InputTag("mix", "MergedCaloTruth"), + label_cp = cms.InputTag("mix","MergedCaloTruth"), +) + +tracksterSimTracksterAssociationPR = cms.EDProducer("TSToSimTSHitLCAssociatorEDProducer", + associator = cms.InputTag('simTracksterHitLCAssociatorByEnergyScoreProducer'), + label_tst = cms.InputTag("ticlTrackstersMerge"), + label_simTst = cms.InputTag("ticlSimTracksters"), + label_lcl = cms.InputTag("hgcalMergeLayerClusters"), + label_scl = cms.InputTag("mix", "MergedCaloTruth"), + label_cp = cms.InputTag("mix","MergedCaloTruth"), +) + + +tracksterSimTracksterAssociationLinkingbyCLUE3D = cms.EDProducer("TSToSimTSHitLCAssociatorEDProducer", + associator = cms.InputTag('simTracksterHitLCAssociatorByEnergyScoreProducer'), + label_tst = cms.InputTag("ticlTrackstersCLUE3DHigh"), + label_simTst = cms.InputTag("ticlSimTracksters", "fromCPs"), + label_lcl = cms.InputTag("hgcalMergeLayerClusters"), + label_scl = cms.InputTag("mix", "MergedCaloTruth"), + label_cp = cms.InputTag("mix","MergedCaloTruth"), +) + +tracksterSimTracksterAssociationPRbyCLUE3D = cms.EDProducer("TSToSimTSHitLCAssociatorEDProducer", + associator = cms.InputTag('simTracksterHitLCAssociatorByEnergyScoreProducer'), + label_tst = cms.InputTag("ticlTrackstersCLUE3DHigh"), label_simTst = cms.InputTag("ticlSimTracksters"), - label_lcl = cms.InputTag("hgcalMergeLayerClusters") + label_lcl = cms.InputTag("hgcalMergeLayerClusters"), + label_scl = cms.InputTag("mix", "MergedCaloTruth"), + label_cp = cms.InputTag("mix","MergedCaloTruth"), ) + +tracksterSimTracksterAssociationLinkingPU = cms.EDProducer("TSToSimTSHitLCAssociatorEDProducer", + associator = cms.InputTag('simTracksterHitLCAssociatorByEnergyScoreProducer'), + label_tst = cms.InputTag("ticlTrackstersMerge"), + label_simTst = cms.InputTag("ticlSimTracksters", "PU"), + label_lcl = cms.InputTag("hgcalMergeLayerClusters"), + label_scl = cms.InputTag("mix", "MergedCaloTruth"), + label_cp = cms.InputTag("mix","MergedCaloTruth"), +) + +tracksterSimTracksterAssociationPRPU = cms.EDProducer("TSToSimTSHitLCAssociatorEDProducer", + associator = cms.InputTag('simTracksterHitLCAssociatorByEnergyScoreProducer'), + label_tst = cms.InputTag("ticlTrackstersMerge"), + label_simTst = cms.InputTag("ticlSimTracksters", "PU"), + label_lcl = cms.InputTag("hgcalMergeLayerClusters"), + label_scl = cms.InputTag("mix", "MergedCaloTruth"), + label_cp = cms.InputTag("mix","MergedCaloTruth"), +) + + diff --git a/SimDataFormats/Associations/interface/TracksterToSimTracksterAssociatorBaseImpl.h b/SimDataFormats/Associations/interface/TracksterToSimTracksterAssociatorBaseImpl.h index c1db97827ef8e..08839cc729cc6 100644 --- a/SimDataFormats/Associations/interface/TracksterToSimTracksterAssociatorBaseImpl.h +++ b/SimDataFormats/Associations/interface/TracksterToSimTracksterAssociatorBaseImpl.h @@ -23,7 +23,7 @@ namespace hgcal { edm::OneToManyWithQualityGeneric>> SimToRecoCollectionSimTracksters; typedef edm::AssociationMap< - edm::OneToManyWithQualityGeneric> + edm::OneToManyWithQualityGeneric>> RecoToSimCollectionSimTracksters; class TracksterToSimTracksterAssociatorBaseImpl { diff --git a/SimDataFormats/Associations/interface/TracksterToSimTracksterHitLCAssociator.h b/SimDataFormats/Associations/interface/TracksterToSimTracksterHitLCAssociator.h new file mode 100644 index 0000000000000..5162e2bd08ec2 --- /dev/null +++ b/SimDataFormats/Associations/interface/TracksterToSimTracksterHitLCAssociator.h @@ -0,0 +1,57 @@ +#ifndef SimDataFormats_Associations_TracksterToSimTracksterHitLCAssociator_h +#define SimDataFormats_Associations_TracksterToSimTracksterHitLCAssociator_h +// Original Author: Leonardo Cristella + +// system include files +#include + +// user include files + +#include "SimDataFormats/Associations/interface/TracksterToSimTracksterHitLCAssociatorBaseImpl.h" + +namespace hgcal { + + class TracksterToSimTracksterHitLCAssociator { + public: + TracksterToSimTracksterHitLCAssociator(std::unique_ptr); + TracksterToSimTracksterHitLCAssociator() = default; + TracksterToSimTracksterHitLCAssociator(TracksterToSimTracksterHitLCAssociator &&) = default; + TracksterToSimTracksterHitLCAssociator &operator=(TracksterToSimTracksterHitLCAssociator &&) = default; + TracksterToSimTracksterHitLCAssociator(const TracksterToSimTracksterHitLCAssociator &) = delete; + const TracksterToSimTracksterHitLCAssociator &operator=(const TracksterToSimTracksterHitLCAssociator &) = delete; + + ~TracksterToSimTracksterHitLCAssociator() = default; + + hgcal::association_t makeConnections(const edm::Handle &tCH, + const edm::Handle &lCCH, + const edm::Handle &sCCH, + const edm::Handle &cPCH, + const edm::Handle &sTCH) const { + return m_impl->makeConnections(tCH, lCCH, sCCH, cPCH, sTCH); + } + /// Associate a Trackster to SimClusters + hgcal::RecoToSimCollectionSimTracksters associateRecoToSim( + const edm::Handle &tCH, + const edm::Handle &lCCH, + const edm::Handle &sCCH, + const edm::Handle &cPCH, + const edm::Handle &sTCH) const { + return m_impl->associateRecoToSim(tCH, lCCH, sCCH, cPCH, sTCH); + }; + + /// Associate a SimCluster to Tracksters + hgcal::SimToRecoCollectionSimTracksters associateSimToReco( + const edm::Handle &tCH, + const edm::Handle &lCCH, + const edm::Handle &sCCH, + const edm::Handle &cPCH, + const edm::Handle &sTCH) const { + return m_impl->associateSimToReco(tCH, lCCH, sCCH, cPCH, sTCH); + } + + private: + std::unique_ptr m_impl; + }; +} // namespace hgcal + +#endif diff --git a/SimDataFormats/Associations/interface/TracksterToSimTracksterHitLCAssociatorBaseImpl.h b/SimDataFormats/Associations/interface/TracksterToSimTracksterHitLCAssociatorBaseImpl.h new file mode 100644 index 0000000000000..405657555f032 --- /dev/null +++ b/SimDataFormats/Associations/interface/TracksterToSimTracksterHitLCAssociatorBaseImpl.h @@ -0,0 +1,58 @@ +#ifndef SimDataFormats_Associations_TracksterToSimTracksterHitLCAssociatorBaseImpl_h +#define SimDataFormats_Associations_TracksterToSimTracksterHitLCAssociatorBaseImpl_h + +#include "DataFormats/Common/interface/Handle.h" +#include "DataFormats/Common/interface/AssociationMap.h" +#include "DataFormats/HGCalReco/interface/Trackster.h" +#include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h" + +#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h" +typedef std::vector SimClusterCollection; +#include "SimDataFormats/CaloAnalysis/interface/CaloParticleFwd.h" + +namespace hgcal { + + enum validationType { Linking = 0, PatternRecognition, PatternRecognition_CP }; + + typedef std::vector>> sharedEnergyAndScore_t; + // This is used to save the simTracksterOnLayer structure for all simTracksters. + // It is not exactly what is returned outside, but out of its entries, the output object is build. + typedef std::tuple association_t; + + typedef edm::AssociationMap< + edm::OneToManyWithQualityGeneric>> + SimToRecoCollectionSimTracksters; + typedef SimToRecoCollectionSimTracksters RecoToSimCollectionSimTracksters; + + class TracksterToSimTracksterHitLCAssociatorBaseImpl { + public: + /// Constructor + TracksterToSimTracksterHitLCAssociatorBaseImpl(); + /// Destructor + virtual ~TracksterToSimTracksterHitLCAssociatorBaseImpl(); + + hgcal::association_t makeConnections(const edm::Handle &tCH, + const edm::Handle &lCCH, + const edm::Handle &sCCH, + const edm::Handle &cPCH, + const edm::Handle &sTCH) const; + + /// Associate a Trackster to SimClusters + virtual hgcal::RecoToSimCollectionSimTracksters associateRecoToSim( + const edm::Handle &tCH, + const edm::Handle &lCCH, + const edm::Handle &sCCH, + const edm::Handle &cPCH, + const edm::Handle &sTCH) const; + + /// Associate a SimCluster to Tracksters + virtual hgcal::SimToRecoCollectionSimTracksters associateSimToReco( + const edm::Handle &tCH, + const edm::Handle &lCCH, + const edm::Handle &sCCH, + const edm::Handle &cPCH, + const edm::Handle &sTCH) const; + }; +} // namespace hgcal + +#endif diff --git a/SimDataFormats/Associations/src/TracksterToSimTracksterHitLCAssociator.cc b/SimDataFormats/Associations/src/TracksterToSimTracksterHitLCAssociator.cc new file mode 100644 index 0000000000000..8ff90ed3d7123 --- /dev/null +++ b/SimDataFormats/Associations/src/TracksterToSimTracksterHitLCAssociator.cc @@ -0,0 +1,5 @@ +#include "SimDataFormats/Associations/interface/TracksterToSimTracksterHitLCAssociator.h" + +hgcal::TracksterToSimTracksterHitLCAssociator::TracksterToSimTracksterHitLCAssociator( + std::unique_ptr ptr) + : m_impl(std::move(ptr)) {} diff --git a/SimDataFormats/Associations/src/TracksterToSimTracksterHitLCAssociatorBaseImpl.cc b/SimDataFormats/Associations/src/TracksterToSimTracksterHitLCAssociatorBaseImpl.cc new file mode 100644 index 0000000000000..6b169b9724e41 --- /dev/null +++ b/SimDataFormats/Associations/src/TracksterToSimTracksterHitLCAssociatorBaseImpl.cc @@ -0,0 +1,34 @@ +#include "SimDataFormats/Associations/interface/TracksterToSimTracksterHitLCAssociatorBaseImpl.h" + +namespace hgcal { + TracksterToSimTracksterHitLCAssociatorBaseImpl::TracksterToSimTracksterHitLCAssociatorBaseImpl(){}; + TracksterToSimTracksterHitLCAssociatorBaseImpl::~TracksterToSimTracksterHitLCAssociatorBaseImpl(){}; + + hgcal::association_t TracksterToSimTracksterHitLCAssociatorBaseImpl::makeConnections( + const edm::Handle &tCH, + const edm::Handle &lCCH, + const edm::Handle &sCCH, + const edm::Handle &cPCH, + const edm::Handle &sTCH) const { + return hgcal::association_t(); + } + + hgcal::RecoToSimCollectionSimTracksters TracksterToSimTracksterHitLCAssociatorBaseImpl::associateRecoToSim( + const edm::Handle &tCH, + const edm::Handle &lCCH, + const edm::Handle &sCCH, + const edm::Handle &cPCH, + const edm::Handle &sTCH) const { + return hgcal::RecoToSimCollectionSimTracksters(); + } + + hgcal::SimToRecoCollectionSimTracksters TracksterToSimTracksterHitLCAssociatorBaseImpl::associateSimToReco( + const edm::Handle &tCH, + const edm::Handle &lCCH, + const edm::Handle &sCCH, + const edm::Handle &cPCH, + const edm::Handle &sTCH) const { + return hgcal::SimToRecoCollectionSimTracksters(); + } + +} // namespace hgcal diff --git a/SimDataFormats/Associations/src/classes.h b/SimDataFormats/Associations/src/classes.h index b3e3f80f3e930..8fa4f488c47a4 100644 --- a/SimDataFormats/Associations/src/classes.h +++ b/SimDataFormats/Associations/src/classes.h @@ -12,6 +12,7 @@ #include "SimDataFormats/Associations/interface/TracksterToSimClusterAssociator.h" #include "SimDataFormats/Associations/interface/MultiClusterToCaloParticleAssociator.h" #include "SimDataFormats/Associations/interface/TracksterToSimTracksterAssociator.h" +#include "SimDataFormats/Associations/interface/TracksterToSimTracksterHitLCAssociator.h" #include "SimDataFormats/Associations/interface/TTTrackTruthPair.h" #include "SimDataFormats/Associations/interface/LayerClusterToSimTracksterAssociator.h" @@ -34,7 +35,9 @@ namespace SimDataFormats_Associations { edm::Wrapper dummy9; - edm::Wrapper dummy10; + edm::Wrapper dummy10; + + edm::Wrapper dummy11; reco::VertexSimToRecoCollection vstrc; reco::VertexSimToRecoCollection::const_iterator vstrci; diff --git a/SimDataFormats/Associations/src/classes_def.xml b/SimDataFormats/Associations/src/classes_def.xml index 83e5364cac2a8..7e50c17559255 100644 --- a/SimDataFormats/Associations/src/classes_def.xml +++ b/SimDataFormats/Associations/src/classes_def.xml @@ -24,6 +24,9 @@ + + + diff --git a/SimDataFormats/Track/BuildFile.xml b/SimDataFormats/Track/BuildFile.xml index a9f8fc715dac6..06527492a4075 100644 --- a/SimDataFormats/Track/BuildFile.xml +++ b/SimDataFormats/Track/BuildFile.xml @@ -1,3 +1,4 @@ + diff --git a/SimDataFormats/Track/interface/UniqueSimTrackId.h b/SimDataFormats/Track/interface/UniqueSimTrackId.h index 57edf4735f572..cb918e44c5d0e 100644 --- a/SimDataFormats/Track/interface/UniqueSimTrackId.h +++ b/SimDataFormats/Track/interface/UniqueSimTrackId.h @@ -2,7 +2,8 @@ #define SimDataFormatsTrackUniqueSimTrackId_H #include "FWCore/Utilities/interface/hash_combine.h" - +#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h" +#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h" #include "SimDataFormats/EncodedEventId/interface/EncodedEventId.h" #include @@ -13,4 +14,7 @@ struct UniqueSimTrackIdHash { } }; +struct SimTrackToTPMap { + std::unordered_map mapping; +}; #endif diff --git a/SimDataFormats/Track/src/classes.h b/SimDataFormats/Track/src/classes.h index 3c4798d038ae2..f7773a626adaf 100644 --- a/SimDataFormats/Track/src/classes.h +++ b/SimDataFormats/Track/src/classes.h @@ -1,6 +1,7 @@ #include "SimDataFormats/Track/interface/CoreSimTrack.h" #include "SimDataFormats/Track/interface/SimTrack.h" #include "SimDataFormats/Track/interface/SimTrackContainer.h" +#include "SimDataFormats/Track/interface/UniqueSimTrackId.h" #include "DataFormats/Common/interface/Wrapper.h" #include diff --git a/SimDataFormats/Track/src/classes_def.xml b/SimDataFormats/Track/src/classes_def.xml index d76801ad762c7..c53804b2299e5 100644 --- a/SimDataFormats/Track/src/classes_def.xml +++ b/SimDataFormats/Track/src/classes_def.xml @@ -13,4 +13,11 @@ + + + + + + + diff --git a/SimGeneral/TrackingAnalysis/plugins/SimHitTPAssociationProducer.cc b/SimGeneral/TrackingAnalysis/plugins/SimHitTPAssociationProducer.cc index e52bb255be8a4..c8658041c088a 100644 --- a/SimGeneral/TrackingAnalysis/plugins/SimHitTPAssociationProducer.cc +++ b/SimGeneral/TrackingAnalysis/plugins/SimHitTPAssociationProducer.cc @@ -27,6 +27,7 @@ SimHitTPAssociationProducer::SimHitTPAssociationProducer(const edm::ParameterSet _trackingParticleSrc( consumes(cfg.getParameter("trackingParticleSrc"))) { produces(); + produces("simTrackToTP"); std::vector tags = cfg.getParameter>("simHitSrc"); _simHitSrc.reserve(tags.size()); for (auto const &tag : tags) { @@ -44,7 +45,7 @@ void SimHitTPAssociationProducer::produce(edm::StreamID, edm::Event &iEvent, con iEvent.getByToken(_trackingParticleSrc, TPCollectionH); // prepare temporary map between SimTrackId and TrackingParticle index - std::unordered_map mapping; + auto simTrackToTPMap = std::make_unique(); auto const &tpColl = *TPCollectionH.product(); for (TrackingParticleCollection::size_type itp = 0, size = tpColl.size(); itp < size; ++itp) { auto const &trackingParticle = tpColl[itp]; @@ -53,7 +54,7 @@ void SimHitTPAssociationProducer::produce(edm::StreamID, edm::Event &iEvent, con EncodedEventId eid(trackingParticle.eventId()); for (auto const &trk : trackingParticle.g4Tracks()) { UniqueSimTrackId trkid(trk.trackId(), eid); - mapping.insert(std::make_pair(trkid, trackingParticleRef)); + simTrackToTPMap->mapping.insert(std::make_pair(trkid, trackingParticleRef)); } } @@ -66,8 +67,8 @@ void SimHitTPAssociationProducer::produce(edm::StreamID, edm::Event &iEvent, con TrackPSimHitRef pSimHitRef(PSimHitCollectionH, psimHitI); auto const &pSimHit = pSimHitCollection[psimHitI]; UniqueSimTrackId simTkIds(pSimHit.trackId(), pSimHit.eventId()); - auto ipos = mapping.find(simTkIds); - if (ipos != mapping.end()) { + auto ipos = simTrackToTPMap->mapping.find(simTkIds); + if (ipos != simTrackToTPMap->mapping.end()) { simHitTPList->emplace_back(ipos->second, pSimHitRef); } } @@ -75,6 +76,7 @@ void SimHitTPAssociationProducer::produce(edm::StreamID, edm::Event &iEvent, con std::sort(simHitTPList->begin(), simHitTPList->end(), simHitTPAssociationListGreater); iEvent.put(std::move(simHitTPList)); + iEvent.put(std::move(simTrackToTPMap), "simTrackToTP"); } #include "FWCore/Framework/interface/MakerMacros.h" diff --git a/SimTracker/Configuration/python/SimTracker_EventContent_cff.py b/SimTracker/Configuration/python/SimTracker_EventContent_cff.py index 7b1b86089b7a0..8b687fbb43944 100644 --- a/SimTracker/Configuration/python/SimTracker_EventContent_cff.py +++ b/SimTracker/Configuration/python/SimTracker_EventContent_cff.py @@ -15,7 +15,8 @@ 'keep *_assocOutInConversionTracks_*_*', 'keep *_assocInOutConversionTracks_*_*', 'keep *_TTClusterAssociatorFromPixelDigis_*_*', - 'keep *_TTStubAssociatorFromPixelDigis_*_*') + 'keep *_TTStubAssociatorFromPixelDigis_*_*', + 'keep *_simHitTPAssocProducer_*_*') ) # For phase2 premixing switch the sim digi collections to the ones including pileup diff --git a/Validation/Configuration/python/hgcalSimValid_cff.py b/Validation/Configuration/python/hgcalSimValid_cff.py index 680caacb61ff6..142fda660ed81 100644 --- a/Validation/Configuration/python/hgcalSimValid_cff.py +++ b/Validation/Configuration/python/hgcalSimValid_cff.py @@ -1,10 +1,15 @@ import FWCore.ParameterSet.Config as cms from SimCalorimetry.HGCalSimProducers.hgcHitAssociation_cfi import lcAssocByEnergyScoreProducer, scAssocByEnergyScoreProducer +from SimCalorimetry.HGCalAssociatorProducers.simTracksterAssociatorByEnergyScore_cfi import simTracksterAssociatorByEnergyScore as simTsAssocByEnergyScoreProducer +from SimCalorimetry.HGCalAssociatorProducers.layerClusterSimTracksterAssociatorByEnergyScore_cfi import layerClusterSimTracksterAssociatorByEnergyScore as lcSimTSAssocByEnergyScoreProducer from SimCalorimetry.HGCalAssociatorProducers.LCToCPAssociation_cfi import layerClusterCaloParticleAssociation as layerClusterCaloParticleAssociationProducer +from SimCalorimetry.HGCalAssociatorProducers.simTracksterHitLCAssociatorByEnergyScore_cfi import simTracksterHitLCAssociatorByEnergyScore as simTracksterHitLCAssociatorByEnergyScoreProducer from SimCalorimetry.HGCalAssociatorProducers.LCToSCAssociation_cfi import layerClusterSimClusterAssociation as layerClusterSimClusterAssociationProducer +from SimCalorimetry.HGCalAssociatorProducers.LCToSimTSAssociation_cfi import layerClusterSimTracksterAssociation as layerClusterSimTracksterAssociationProducer from SimCalorimetry.HGCalAssociatorProducers.LCToCPAssociation_cfi import layerClusterCaloParticleAssociationHFNose as layerClusterCaloParticleAssociationProducerHFNose from SimCalorimetry.HGCalAssociatorProducers.LCToSCAssociation_cfi import layerClusterSimClusterAssociationHFNose as layerClusterSimClusterAssociationProducerHFNose +from SimCalorimetry.HGCalAssociatorProducers.TSToSimTSAssociation_cfi import tracksterSimTracksterAssociationLinking, tracksterSimTracksterAssociationPR,tracksterSimTracksterAssociationLinkingbyCLUE3D, tracksterSimTracksterAssociationPRbyCLUE3D, tracksterSimTracksterAssociationLinkingPU, tracksterSimTracksterAssociationPRPU from Validation.HGCalValidation.simhitValidation_cff import * from Validation.HGCalValidation.digiValidation_cff import * @@ -12,6 +17,7 @@ from Validation.HGCalValidation.hgcalHitValidation_cff import * from RecoHGCal.TICL.SimTracksters_cff import * + from Validation.HGCalValidation.HGCalValidator_cfi import hgcalValidator from Validation.RecoParticleFlow.PFJetValidation_cff import pfJetValidation1 as _hgcalPFJetValidation @@ -28,6 +34,11 @@ hgcalAssociators = cms.Task(lcAssocByEnergyScoreProducer, layerClusterCaloParticleAssociationProducer, scAssocByEnergyScoreProducer, layerClusterSimClusterAssociationProducer, + lcSimTSAssocByEnergyScoreProducer, layerClusterSimTracksterAssociationProducer, + simTsAssocByEnergyScoreProducer, simTracksterHitLCAssociatorByEnergyScoreProducer, + tracksterSimTracksterAssociationLinking, tracksterSimTracksterAssociationPR, + tracksterSimTracksterAssociationLinkingbyCLUE3D, tracksterSimTracksterAssociationPRbyCLUE3D, + tracksterSimTracksterAssociationLinkingPU, tracksterSimTracksterAssociationPRPU ) hgcalValidation = cms.Sequence(hgcalSimHitValidationEE From 5a4c78e9baf873a7841e298091cb1d61e9f2066b Mon Sep 17 00:00:00 2001 From: Wahid Redjeb Date: Wed, 20 Sep 2023 10:37:21 +0200 Subject: [PATCH 02/12] add getLayerType to recHitTools --- RecoHGCal/TICL/plugins/TICLDumper.cc | 2 +- .../HGCalRecAlgos/interface/RecHitTools.h | 1 + .../HGCalRecAlgos/src/RecHitTools.cc | 45 +++++++++++++++++++ 3 files changed, 47 insertions(+), 1 deletion(-) diff --git a/RecoHGCal/TICL/plugins/TICLDumper.cc b/RecoHGCal/TICL/plugins/TICLDumper.cc index 7897679e35e45..f8ce64f56f66c 100644 --- a/RecoHGCal/TICL/plugins/TICLDumper.cc +++ b/RecoHGCal/TICL/plugins/TICLDumper.cc @@ -1688,7 +1688,7 @@ void TICLDumper::analyze(const edm::Event& event, const edm::EventSetup& setup) cluster_layer_id.push_back(layerId); uint32_t number_of_hits = cluster_iterator->hitsAndFractions().size(); cluster_number_of_hits.push_back(number_of_hits); - cluster_type.push_back(ticl::returnClusterType(lc_seed, rhtools_)); + cluster_type.push_back(rhtools_.getLayerType(lc_seed)); cluster_timeErr.push_back(layerClustersTimes.get(c_id).second); cluster_time.push_back(layerClustersTimes.get(c_id).first); c_id += 1; diff --git a/RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h b/RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h index 987d4c6db89a9..891beb4de7f81 100644 --- a/RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h +++ b/RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h @@ -51,6 +51,7 @@ namespace hgcal { unsigned int getLayer(ForwardSubdetector type) const; unsigned int getLayer(const DetId&) const; unsigned int getLayerWithOffset(const DetId&) const; + int getLayerType(const DetId& id) const; std::pair getWafer(const DetId&) const; std::pair getCell(const DetId&) const; diff --git a/RecoLocalCalo/HGCalRecAlgos/src/RecHitTools.cc b/RecoLocalCalo/HGCalRecAlgos/src/RecHitTools.cc index e55845499bb28..34787bbebfa9a 100644 --- a/RecoLocalCalo/HGCalRecAlgos/src/RecHitTools.cc +++ b/RecoLocalCalo/HGCalRecAlgos/src/RecHitTools.cc @@ -63,6 +63,17 @@ namespace { return ddd; } + enum LayerType { + CE_E_120 = 0, + CE_E_200 = 1, + CE_E_300 = 2, + CE_H_120 = 3, + CE_H_200 = 4, + CE_H_300 = 5, + CE_H_SCINT = 6, + EnumSize = 7 + }; + } // namespace void RecHitTools::setGeometry(const CaloGeometry& geom) { @@ -425,6 +436,40 @@ bool RecHitTools::isHalfCell(const DetId& id) const { return ishalf; } +int RecHitTools::getLayerType(const DetId& id) const { + auto layer_number = getLayerWithOffset(id); + auto thickness = getSiThickIndex(id); + auto geomNose = + static_cast(geom_->getSubdetectorGeometry(DetId::Forward, ForwardSubdetector::HFNose)); + auto isNose = geomNose ? true : false; + auto isEELayer = (layer_number <= lastLayerEE(isNose)); + auto isScint = isScintillator(id); + int layerType = -1; + + if (isScint) { + layerType = CE_H_SCINT; + } + if (isEELayer) { + if (thickness == 0) { + layerType = CE_E_120; + } else if (thickness == 1) { + layerType = CE_E_200; + } else if (thickness == 2) { + layerType = CE_E_300; + } + } else { + if (thickness == 0) { + layerType = CE_H_120; + } else if (thickness == 1) { + layerType = CE_H_200; + } else if (thickness == 2) { + layerType = CE_H_300; + } + } + assert(layerType != -1); + return layerType; +} + bool RecHitTools::isSilicon(const DetId& id) const { return (id.det() == DetId::HGCalEE || id.det() == DetId::HGCalHSi || (id.det() == DetId::Forward && id.subdetId() == static_cast(HFNose))); From 74754892edc81186b40411a2a5dab167ba721948 Mon Sep 17 00:00:00 2001 From: Wahid Redjeb Date: Mon, 25 Sep 2023 13:19:46 +0200 Subject: [PATCH 03/12] fix indentation --- SimDataFormats/Track/src/classes_def.xml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/SimDataFormats/Track/src/classes_def.xml b/SimDataFormats/Track/src/classes_def.xml index c53804b2299e5..f3548ad1215ff 100644 --- a/SimDataFormats/Track/src/classes_def.xml +++ b/SimDataFormats/Track/src/classes_def.xml @@ -13,9 +13,9 @@ - - - + + + From bab1995574ab156cffbbcd7478546aaebafdeccb Mon Sep 17 00:00:00 2001 From: Wahid Redjeb Date: Tue, 26 Sep 2023 19:51:34 +0200 Subject: [PATCH 04/12] fix dicionary for UniqueSimTrackId --- RecoHGCal/TICL/plugins/SimTrackstersProducer.cc | 3 ++- SimDataFormats/Track/src/classes.h | 1 - SimDataFormats/Track/src/classes_def.xml | 7 ------- .../interface/UniqueSimTrackId.h | 5 +++-- SimDataFormats/TrackingAnalysis/src/classes.h | 1 + SimDataFormats/TrackingAnalysis/src/classes_def.xml | 6 ++++++ .../plugins/SimHitTPAssociationProducer.cc | 2 +- 7 files changed, 13 insertions(+), 12 deletions(-) rename SimDataFormats/{Track => TrackingAnalysis}/interface/UniqueSimTrackId.h (78%) diff --git a/RecoHGCal/TICL/plugins/SimTrackstersProducer.cc b/RecoHGCal/TICL/plugins/SimTrackstersProducer.cc index 34eae24cb0c07..c062f36a48fa8 100644 --- a/RecoHGCal/TICL/plugins/SimTrackstersProducer.cc +++ b/RecoHGCal/TICL/plugins/SimTrackstersProducer.cc @@ -30,10 +30,11 @@ #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h" #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h" -#include "SimDataFormats/Track/interface/UniqueSimTrackId.h" #include "SimDataFormats/Vertex/interface/SimVertex.h" #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h" +#include "SimDataFormats/TrackingAnalysis/interface/UniqueSimTrackId.h" + #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h" #include "RecoHGCal/TICL/interface/commons.h" diff --git a/SimDataFormats/Track/src/classes.h b/SimDataFormats/Track/src/classes.h index f7773a626adaf..3c4798d038ae2 100644 --- a/SimDataFormats/Track/src/classes.h +++ b/SimDataFormats/Track/src/classes.h @@ -1,7 +1,6 @@ #include "SimDataFormats/Track/interface/CoreSimTrack.h" #include "SimDataFormats/Track/interface/SimTrack.h" #include "SimDataFormats/Track/interface/SimTrackContainer.h" -#include "SimDataFormats/Track/interface/UniqueSimTrackId.h" #include "DataFormats/Common/interface/Wrapper.h" #include diff --git a/SimDataFormats/Track/src/classes_def.xml b/SimDataFormats/Track/src/classes_def.xml index f3548ad1215ff..d76801ad762c7 100644 --- a/SimDataFormats/Track/src/classes_def.xml +++ b/SimDataFormats/Track/src/classes_def.xml @@ -13,11 +13,4 @@ - - - - - - - diff --git a/SimDataFormats/Track/interface/UniqueSimTrackId.h b/SimDataFormats/TrackingAnalysis/interface/UniqueSimTrackId.h similarity index 78% rename from SimDataFormats/Track/interface/UniqueSimTrackId.h rename to SimDataFormats/TrackingAnalysis/interface/UniqueSimTrackId.h index cb918e44c5d0e..0ee529ae2383b 100644 --- a/SimDataFormats/Track/interface/UniqueSimTrackId.h +++ b/SimDataFormats/TrackingAnalysis/interface/UniqueSimTrackId.h @@ -4,10 +4,10 @@ #include "FWCore/Utilities/interface/hash_combine.h" #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h" #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h" -#include "SimDataFormats/EncodedEventId/interface/EncodedEventId.h" #include using UniqueSimTrackId = std::pair; + struct UniqueSimTrackIdHash { std::size_t operator()(UniqueSimTrackId const &s) const noexcept { return edm::hash_value(s.first, s.second.rawId()); @@ -15,6 +15,7 @@ struct UniqueSimTrackIdHash { }; struct SimTrackToTPMap { - std::unordered_map mapping; + std::unordered_map mapping; }; + #endif diff --git a/SimDataFormats/TrackingAnalysis/src/classes.h b/SimDataFormats/TrackingAnalysis/src/classes.h index 8f99fdd31214d..ff56da99e359b 100644 --- a/SimDataFormats/TrackingAnalysis/src/classes.h +++ b/SimDataFormats/TrackingAnalysis/src/classes.h @@ -2,6 +2,7 @@ #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h" #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h" #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h" +#include "SimDataFormats/TrackingAnalysis/interface/UniqueSimTrackId.h" #include "DataFormats/TrackReco/interface/Track.h" #include "DataFormats/Common/interface/AssociationMapHelpers.h" #include "DataFormats/Common/interface/Wrapper.h" diff --git a/SimDataFormats/TrackingAnalysis/src/classes_def.xml b/SimDataFormats/TrackingAnalysis/src/classes_def.xml index 2d8f00696a58f..713960a5c64e2 100755 --- a/SimDataFormats/TrackingAnalysis/src/classes_def.xml +++ b/SimDataFormats/TrackingAnalysis/src/classes_def.xml @@ -34,4 +34,10 @@ + + + + + + diff --git a/SimGeneral/TrackingAnalysis/plugins/SimHitTPAssociationProducer.cc b/SimGeneral/TrackingAnalysis/plugins/SimHitTPAssociationProducer.cc index c8658041c088a..d6cc8a4aaba68 100644 --- a/SimGeneral/TrackingAnalysis/plugins/SimHitTPAssociationProducer.cc +++ b/SimGeneral/TrackingAnalysis/plugins/SimHitTPAssociationProducer.cc @@ -15,7 +15,7 @@ #include "FWCore/Utilities/interface/EDGetToken.h" #include "DataFormats/Common/interface/Handle.h" -#include "SimDataFormats/Track/interface/UniqueSimTrackId.h" +#include "SimDataFormats/TrackingAnalysis/interface/UniqueSimTrackId.h" #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h" #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h" #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h" From 432e161f0b5301135a3f7e79da8b4756ab5aff0e Mon Sep 17 00:00:00 2001 From: Wahid Redjeb Date: Tue, 26 Sep 2023 19:55:06 +0200 Subject: [PATCH 05/12] code format --- SimDataFormats/TrackingAnalysis/interface/UniqueSimTrackId.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimDataFormats/TrackingAnalysis/interface/UniqueSimTrackId.h b/SimDataFormats/TrackingAnalysis/interface/UniqueSimTrackId.h index 0ee529ae2383b..63b080e751d41 100644 --- a/SimDataFormats/TrackingAnalysis/interface/UniqueSimTrackId.h +++ b/SimDataFormats/TrackingAnalysis/interface/UniqueSimTrackId.h @@ -15,7 +15,7 @@ struct UniqueSimTrackIdHash { }; struct SimTrackToTPMap { - std::unordered_map mapping; + std::unordered_map mapping; }; #endif From 10f22afd12147531360207d38c6c07497ebe9d88 Mon Sep 17 00:00:00 2001 From: Wahid Redjeb Date: Tue, 26 Sep 2023 20:53:50 +0200 Subject: [PATCH 06/12] remove hepmc dependency from SimDataFormats/Track --- SimDataFormats/Track/BuildFile.xml | 1 - 1 file changed, 1 deletion(-) diff --git a/SimDataFormats/Track/BuildFile.xml b/SimDataFormats/Track/BuildFile.xml index 06527492a4075..a9f8fc715dac6 100644 --- a/SimDataFormats/Track/BuildFile.xml +++ b/SimDataFormats/Track/BuildFile.xml @@ -1,4 +1,3 @@ - From 7ab5acc9eca6849f7cd7e2e8c746d5fcb9173c80 Mon Sep 17 00:00:00 2001 From: Wahid Redjeb Date: Tue, 26 Sep 2023 22:04:59 +0200 Subject: [PATCH 07/12] fix UniqueSimTrackid includes --- SimTracker/TrackAssociation/plugins/DigiSimLinkPruner.cc | 2 +- .../plugins/ClusterTPAssociationProducer.cc | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/SimTracker/TrackAssociation/plugins/DigiSimLinkPruner.cc b/SimTracker/TrackAssociation/plugins/DigiSimLinkPruner.cc index c72c4e3373c12..56df949182417 100644 --- a/SimTracker/TrackAssociation/plugins/DigiSimLinkPruner.cc +++ b/SimTracker/TrackAssociation/plugins/DigiSimLinkPruner.cc @@ -30,10 +30,10 @@ #include "FWCore/Utilities/interface/StreamID.h" #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h" +#include "SimDataFormats/TrackingAnalysis/interface/UniqueSimTrackId.h" #include "DataFormats/Common/interface/DetSetVector.h" #include "SimDataFormats/TrackerDigiSimLink/interface/PixelDigiSimLink.h" #include "SimDataFormats/TrackerDigiSimLink/interface/StripDigiSimLink.h" -#include "SimDataFormats/Track/interface/UniqueSimTrackId.h" #include diff --git a/SimTracker/TrackerHitAssociation/plugins/ClusterTPAssociationProducer.cc b/SimTracker/TrackerHitAssociation/plugins/ClusterTPAssociationProducer.cc index a4c7e49c539d0..33c8965f4d029 100644 --- a/SimTracker/TrackerHitAssociation/plugins/ClusterTPAssociationProducer.cc +++ b/SimTracker/TrackerHitAssociation/plugins/ClusterTPAssociationProducer.cc @@ -22,12 +22,12 @@ #include "DataFormats/Phase2TrackerCluster/interface/Phase2TrackerCluster1D.h" #include "SimDataFormats/Track/interface/SimTrackContainer.h" -#include "SimDataFormats/Track/interface/UniqueSimTrackId.h" #include "SimDataFormats/TrackerDigiSimLink/interface/StripDigiSimLink.h" #include "SimDataFormats/TrackerDigiSimLink/interface/PixelDigiSimLink.h" #include "DataFormats/Phase2TrackerDigi/interface/Phase2TrackerDigi.h" #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h" #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h" +#include "SimDataFormats/TrackingAnalysis/interface/UniqueSimTrackId.h" #include "SimTracker/TrackerHitAssociation/interface/ClusterTPAssociation.h" namespace { From 3b29edbdee130120c8265a815feee1e072a57327 Mon Sep 17 00:00:00 2001 From: Wahid Redjeb Date: Wed, 27 Sep 2023 10:47:00 +0200 Subject: [PATCH 08/12] fix LogDebug in TSToSimTSAssociatorByEnergyScoreImpl --- .../plugins/TSToSimTSAssociatorByEnergyScoreImpl.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSAssociatorByEnergyScoreImpl.cc b/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSAssociatorByEnergyScoreImpl.cc index 6494153d27438..d5307c9f94cab 100644 --- a/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSAssociatorByEnergyScoreImpl.cc +++ b/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSAssociatorByEnergyScoreImpl.cc @@ -445,7 +445,7 @@ hgcal::RecoToSimCollectionSimTracksters TSToSimTSAssociatorByEnergyScoreImpl::as for (size_t tsId = 0; tsId < stsInTrackster.size(); ++tsId) { for (auto& stPair : stsInTrackster[tsId]) { LogDebug("TSToSimTSAssociatorByEnergyScoreImpl") << "Trackster Id:\t" << tsId << "\tSimTrackster id:\t" - << stPair.first << "\tscore:\t" << stPair.second << "\n"; + << stPair.first << "\tscore:\t" << stPair.second.second << "\n"; // Fill AssociationMap returnValue.insert( edm::Ref(tCH, tsId), // Ref to TS From b1d3b46cb00810dd5103ec18f67b400878fcd6f3 Mon Sep 17 00:00:00 2001 From: Wahid Redjeb Date: Mon, 2 Oct 2023 17:05:35 +0200 Subject: [PATCH 09/12] remove returnClusterType function --- RecoHGCal/TICL/interface/commons.h | 52 ------------------------------ 1 file changed, 52 deletions(-) diff --git a/RecoHGCal/TICL/interface/commons.h b/RecoHGCal/TICL/interface/commons.h index 86b92f4c925fe..528687d8e063c 100644 --- a/RecoHGCal/TICL/interface/commons.h +++ b/RecoHGCal/TICL/interface/commons.h @@ -12,20 +12,6 @@ namespace ticl { constexpr double mpion = 0.13957; constexpr float mpion2 = mpion * mpion; typedef math::XYZVectorF Vector; - enum LayerType { - - CE_E_120 = 0, - CE_E_200 = 1, - CE_E_300 = 2, - CE_H_120_F = 3, - CE_H_200_F = 4, - CE_H_300_F = 5, - CE_H_200_C = 6, - CE_H_300_C = 7, - CE_H_SCINT_C = 8, - EnumSize = 9 - - }; inline Trackster::ParticleType tracksterParticleTypeFromPdgId(int pdgId, int charge) { if (pdgId == 111) { @@ -56,44 +42,6 @@ namespace ticl { // verbosity levels for ticl algorithms enum VerbosityLevel { None = 0, Basic, Advanced, Expert, Guru }; - inline int returnClusterType(DetId& lc_seed, const hgcal::RecHitTools& rhtools_) { - auto layer_number = rhtools_.getLayerWithOffset(lc_seed); - auto thickness = rhtools_.getSiThickIndex(lc_seed); - auto isEELayer = (layer_number <= rhtools_.lastLayerEE(false)); - auto isScintillator = rhtools_.isScintillator(lc_seed); - auto isFine = (layer_number <= rhtools_.lastLayerEE(false) + 7); - - if (isScintillator) { - return CE_H_SCINT_C; - } - if (isEELayer) { - if (thickness == 0) { - return CE_E_120; - } else if (thickness == 1) { - return CE_E_200; - } else if (thickness == 2) { - return CE_E_300; - } - } else { - if (isFine) { - if (thickness == 0) { - return CE_H_120_F; - } else if (thickness == 1) { - return CE_H_200_F; - } else if (thickness == 2) { - return CE_H_300_F; - } - } else { - if (thickness == 1) { - return CE_H_200_C; - } else if (thickness == 2) { - return CE_H_300_C; - } - } - } - return -1; - }; - } // namespace ticl #endif From f1f4d7ca40e6d6d14ae3985e8844abf5bb20bd76 Mon Sep 17 00:00:00 2001 From: Wahid Redjeb Date: Mon, 2 Oct 2023 17:06:34 +0200 Subject: [PATCH 10/12] add assert to check on multiplicity --- .../plugins/TSToSimTSHitLCAssociatorByEnergyScoreImpl.cc | 3 +++ 1 file changed, 3 insertions(+) diff --git a/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorByEnergyScoreImpl.cc b/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorByEnergyScoreImpl.cc index d3bac4802178d..04673620ad29d 100644 --- a/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorByEnergyScoreImpl.cc +++ b/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorByEnergyScoreImpl.cc @@ -72,8 +72,10 @@ hgcal::association_t TSToSimTSHitLCAssociatorByEnergyScoreImpl::makeConnections( const auto& lcsInSimTrackster = simTSs[i].vertices(); const auto& multiplicities = simTSs[i].vertex_multiplicity(); for (size_t j = 0; j < lcsInSimTrackster.size(); ++j) { + assert(multiplicities[j] > 0.f); const auto& v = lcsInSimTrackster[j]; float fraction = 1.f / multiplicities[j]; + for (const auto& haf : layerClusters[v].hitsAndFractions()) { detIdSimTSId_Map[haf.first].emplace_back(i, haf.second * fraction); } @@ -86,6 +88,7 @@ hgcal::association_t TSToSimTSHitLCAssociatorByEnergyScoreImpl::makeConnections( const auto& lcsInSimTrackster = tracksters[i].vertices(); const auto& multiplicities = tracksters[i].vertex_multiplicity(); for (size_t j = 0; j < lcsInSimTrackster.size(); ++j) { + assert(multiplicities[j] > 0.f); const auto& v = lcsInSimTrackster[j]; float fraction = 1.f / multiplicities[j]; for (const auto& haf : layerClusters[v].hitsAndFractions()) { From b2a002b3e2f4ffce80cf41d4bff3f722a5981906 Mon Sep 17 00:00:00 2001 From: Wahid Redjeb Date: Mon, 2 Oct 2023 17:19:34 +0200 Subject: [PATCH 11/12] code-format --- .../plugins/TSToSimTSHitLCAssociatorByEnergyScoreImpl.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorByEnergyScoreImpl.cc b/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorByEnergyScoreImpl.cc index 04673620ad29d..0b68d6772b667 100644 --- a/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorByEnergyScoreImpl.cc +++ b/SimCalorimetry/HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorByEnergyScoreImpl.cc @@ -75,7 +75,6 @@ hgcal::association_t TSToSimTSHitLCAssociatorByEnergyScoreImpl::makeConnections( assert(multiplicities[j] > 0.f); const auto& v = lcsInSimTrackster[j]; float fraction = 1.f / multiplicities[j]; - for (const auto& haf : layerClusters[v].hitsAndFractions()) { detIdSimTSId_Map[haf.first].emplace_back(i, haf.second * fraction); } From a98834cadd6dce32df3f130f3d93d85c179c771f Mon Sep 17 00:00:00 2001 From: Wahid Redjeb Date: Mon, 9 Oct 2023 17:15:13 +0200 Subject: [PATCH 12/12] rename function getLayerType to getCellType --- RecoHGCal/TICL/plugins/TICLDumper.cc | 2 +- RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h | 2 +- RecoLocalCalo/HGCalRecAlgos/src/RecHitTools.cc | 4 ++-- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/RecoHGCal/TICL/plugins/TICLDumper.cc b/RecoHGCal/TICL/plugins/TICLDumper.cc index f8ce64f56f66c..dc0c5741d42fb 100644 --- a/RecoHGCal/TICL/plugins/TICLDumper.cc +++ b/RecoHGCal/TICL/plugins/TICLDumper.cc @@ -1688,7 +1688,7 @@ void TICLDumper::analyze(const edm::Event& event, const edm::EventSetup& setup) cluster_layer_id.push_back(layerId); uint32_t number_of_hits = cluster_iterator->hitsAndFractions().size(); cluster_number_of_hits.push_back(number_of_hits); - cluster_type.push_back(rhtools_.getLayerType(lc_seed)); + cluster_type.push_back(rhtools_.getCellType(lc_seed)); cluster_timeErr.push_back(layerClustersTimes.get(c_id).second); cluster_time.push_back(layerClustersTimes.get(c_id).first); c_id += 1; diff --git a/RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h b/RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h index 891beb4de7f81..49924b72d054b 100644 --- a/RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h +++ b/RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h @@ -51,7 +51,7 @@ namespace hgcal { unsigned int getLayer(ForwardSubdetector type) const; unsigned int getLayer(const DetId&) const; unsigned int getLayerWithOffset(const DetId&) const; - int getLayerType(const DetId& id) const; + int getCellType(const DetId& id) const; std::pair getWafer(const DetId&) const; std::pair getCell(const DetId&) const; diff --git a/RecoLocalCalo/HGCalRecAlgos/src/RecHitTools.cc b/RecoLocalCalo/HGCalRecAlgos/src/RecHitTools.cc index 34787bbebfa9a..2477e728795f5 100644 --- a/RecoLocalCalo/HGCalRecAlgos/src/RecHitTools.cc +++ b/RecoLocalCalo/HGCalRecAlgos/src/RecHitTools.cc @@ -63,7 +63,7 @@ namespace { return ddd; } - enum LayerType { + enum CellType { CE_E_120 = 0, CE_E_200 = 1, CE_E_300 = 2, @@ -436,7 +436,7 @@ bool RecHitTools::isHalfCell(const DetId& id) const { return ishalf; } -int RecHitTools::getLayerType(const DetId& id) const { +int RecHitTools::getCellType(const DetId& id) const { auto layer_number = getLayerWithOffset(id); auto thickness = getSiThickIndex(id); auto geomNose =