diff --git a/DPGAnalysis/HGCalNanoAOD/plugins/BuildFile.xml b/DPGAnalysis/HGCalNanoAOD/plugins/BuildFile.xml index 314b2fcab95ae..a947945653565 100644 --- a/DPGAnalysis/HGCalNanoAOD/plugins/BuildFile.xml +++ b/DPGAnalysis/HGCalNanoAOD/plugins/BuildFile.xml @@ -5,6 +5,9 @@ + + + diff --git a/DPGAnalysis/HGCalNanoAOD/plugins/SimTracksterTableProducer.cc b/DPGAnalysis/HGCalNanoAOD/plugins/SimTracksterTableProducer.cc index 03134e1c0d065..85941f9541516 100644 --- a/DPGAnalysis/HGCalNanoAOD/plugins/SimTracksterTableProducer.cc +++ b/DPGAnalysis/HGCalNanoAOD/plugins/SimTracksterTableProducer.cc @@ -48,6 +48,7 @@ class SimTracksterTableProducer : public edm::global::EDProducer<> { const size_t nSimTracksters = simTrackstersHandle.isValid() ? simTrackstersHandle->size() : 0; static constexpr float default_value = std::numeric_limits::quiet_NaN(); + static constexpr int default_int_value = -1; std::vector boundaryX(nSimTracksters, default_value); std::vector boundaryY(nSimTracksters, default_value); @@ -61,6 +62,8 @@ class SimTracksterTableProducer : public edm::global::EDProducer<> { std::vector simTime(nSimTracksters, default_value); std::vector genPt(nSimTracksters, default_value); std::vector mass(nSimTracksters, default_value); + std::vector caloParticleIdx(nSimTracksters, default_int_value); + std::vector isPU(nSimTracksters, default_int_value); if ((simTrackstersHandle.isValid() && caloParticlesHandle.isValid() && simClustersHandle.isValid() && cpToSCMapHandle.isValid()) || @@ -71,7 +74,7 @@ class SimTracksterTableProducer : public edm::global::EDProducer<> { const auto& cpToSCMap = *cpToSCMapHandle; //utility lambda for filling vectors - auto fillVectors = [&](const auto& obj, size_t iSim, float time) { + auto fillVectors = [&](const auto& obj, size_t iSim, float time, int cpIdx) { const auto& simTrack = obj.g4Tracks()[0]; const auto caloPt = obj.pt(); const auto simHitSumEnergy = obj.simEnergy(); @@ -90,6 +93,12 @@ class SimTracksterTableProducer : public edm::global::EDProducer<> { simEnergy[iSim] = simHitSumEnergy; genPt[iSim] = caloPt; mass[iSim] = caloMass; + + // PU flag: a non-zero event/bunch-crossing on the + // seed g4Track means the CaloParticle/SimCluster came from a pileup interaction. + isPU[iSim] = (simTrack.eventId().event() != 0 || simTrack.eventId().bunchCrossing() != 0) ? 1 : 0; + + caloParticleIdx[iSim] = cpIdx; }; for (size_t iSim = 0; iSim < simTracksters.size(); ++iSim) { @@ -97,19 +106,22 @@ class SimTracksterTableProducer : public edm::global::EDProducer<> { float time = default_value; if (simT.seedID() == caloParticlesHandle.id()) { + const int cpIdx = static_cast(simT.seedIndex()); const auto& cp = caloParticles[simT.seedIndex()]; time = cp.simTime(); - fillVectors(cp, iSim, time); + fillVectors(cp, iSim, time, cpIdx); } else { const auto& sc = simClusters[simT.seedIndex()]; + int cpIdx = default_int_value; //SCtoCP map not availalbe, use CPtoSC map instead - for (const auto& [cpIdx, scVec] : cpToSCMap) { + for (const auto& [cpIdxCandidate, scVec] : cpToSCMap) { if (std::ranges::find(scVec, simT.seedIndex()) != scVec.end()) { - time = caloParticles[cpIdx].simTime(); + cpIdx = static_cast(cpIdxCandidate); + time = caloParticles[cpIdxCandidate].simTime(); break; //dont need to check further } } - fillVectors(sc, iSim, time); + fillVectors(sc, iSim, time, cpIdx); } } } @@ -125,7 +137,7 @@ class SimTracksterTableProducer : public edm::global::EDProducer<> { simTrackstersTable->addColumn( "boundaryEta", boundaryEta, "CaloVolume boundary pseudorapidity of associated Simobject", precision_); simTrackstersTable->addColumn( - "boundaryPhi", boundaryEta, "CaloVolume boundary phi of associated Simobject", precision_); + "boundaryPhi", boundaryPhi, "CaloVolume boundary phi of associated Simobject", precision_); simTrackstersTable->addColumn( "boundaryPx", boundaryPx, "X component of momentum at CaloVolume boundary of associated Simobject", precision_); simTrackstersTable->addColumn( @@ -135,6 +147,10 @@ class SimTracksterTableProducer : public edm::global::EDProducer<> { simTrackstersTable->addColumn("simTime", simTime, "Sim-Time of simulated object [ns]", precision_); simTrackstersTable->addColumn("genPt", genPt, "Gen-pT associated with SimObject", precision_); simTrackstersTable->addColumn("mass", mass, "mass associated with SimObject", precision_); + simTrackstersTable->addColumn( + "caloParticleIdx", caloParticleIdx, "Index of the parent CaloParticle in the CaloParticle collection"); + simTrackstersTable->addColumn( + "isPU", isPU, "PU flag of the associated Simobject: 1 = pileup (non-zero event/bx), 0 = otherwise"); event.put(std::move(simTrackstersTable), tableName_); } diff --git a/DPGAnalysis/HGCalNanoAOD/plugins/TICLCandidateExtraTableProducer.cc b/DPGAnalysis/HGCalNanoAOD/plugins/TICLCandidateExtraTableProducer.cc index e1bd73f5d0f23..acddb1891c018 100644 --- a/DPGAnalysis/HGCalNanoAOD/plugins/TICLCandidateExtraTableProducer.cc +++ b/DPGAnalysis/HGCalNanoAOD/plugins/TICLCandidateExtraTableProducer.cc @@ -1,3 +1,5 @@ +#include +#include #include "PhysicsTools/NanoAOD/interface/SimpleFlatTableProducer.h" #include "DataFormats/HGCalReco/interface/TICLCandidate.h" #include "DataFormats/HGCalReco/interface/Trackster.h" @@ -16,11 +18,7 @@ #include "TrackingTools/GeomPropagators/interface/Propagator.h" #include "TrackingTools/Records/interface/TrackingComponentsRecord.h" #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h" - -// -// One-to-many: TICLCandidate -> linked Tracksters -// Or SimTICLCandidate --> SimTracksters -// +#include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h" class TICLCandidateExtraTableProducer : public SimpleFlatTableProducerBase> { public: @@ -31,6 +29,20 @@ class TICLCandidateExtraTableProducer : public SimpleFlatTableProducerBase>(params), tracksters_token_(consumes>(params.getParameter("tracksters"))), tracks_token_(consumes>(params.getParameter("tracks"))), + hasLinkedTracksters_(params.existsAs("linkedTracksters")), + linkedTracksters_token_(hasLinkedTracksters_ ? consumes>>( + params.getParameter("linkedTracksters")) + : edm::EDGetTokenT>>()), + hasPUInfo_(params.existsAs("caloParticles") && + params.existsAs("caloParticleToSimClustersMap")), + caloParticles_token_( + hasPUInfo_ ? consumes>(params.getParameter("caloParticles")) + : edm::EDGetTokenT>()), + caloParticleToSimClustersMap_token_( + hasPUInfo_ ? consumes>>( + params.getParameter("caloParticleToSimClustersMap")) + : edm::EDGetTokenT>>()), + produceGeneralTrackBoundary_(params.getParameter("produceGeneralTrackBoundary")), detector_(params.getParameter("detector")), propName_(params.getParameter("propagator")), geometry_token_(esConsumes()), @@ -38,6 +50,27 @@ class TICLCandidateExtraTableProducer : public SimpleFlatTableProducerBase(edm::ESInputTag("", propName_))), hdc_token_(esConsumes( edm::ESInputTag("", (detector_ == "HFNose") ? "HGCalHFNoseSensitive" : "HGCalEESensitive"))) { + // Only declare/emit the linkedTracksters product for instances that actually have it + // configured (reco TICLCandidates). Otherwise, two module instances of this same class + // (reco + sim) would both try to produce a main FlatTable literally named + // "linkedTracksters", which NanoAOD's output module rejects as "multiple main tables". + if (hasLinkedTracksters_) { + produces("linkedTracksters"); + } + + // Full track-index / GSF-track-index lists per candidate (all tracks, not just the first + // via trackPtr()/gsftrackPtr()). Always produced by every instance; the product name is + // prefixed with this->name_ ("TICLCandidates"/"SimTICLCandidates") so the two module + // instances never collide on the same product name. + produces(this->name_ + "TrackIdxs"); + produces(this->name_ + "GsfTrackIdxs"); + + // GeneralTrack HGCal-boundary extension: independent of TICLCandidate linking, one row per + // track in the full tracks collection. Optional enable on exactly one module instance. + if (produceGeneralTrackBoundary_) { + produces("trackBoundary"); + } + if (params.existsAs("collectionVariables")) { edm::ParameterSet const& collectionVarsPSet = params.getParameter("collectionVariables"); for (const auto& coltablename : collectionVarsPSet.getParameterNamesForType()) { @@ -89,6 +122,48 @@ class TICLCandidateExtraTableProducer : public SimpleFlatTableProducerBaseaddColumn("nLinkedTracksters", + std::vector(table_size, 0), + "Number of Tracksters linked to candidate before final linking/merging"); + + auto linkedTable = std::make_unique(0, "linkedTracksters", false, false); + std::vector emptyLinkedIdx; + linkedTable->addColumn("tracksterIndex", emptyLinkedIdx, "Index of linked Trackster"); + linkedTable->setDoc("Tracksters linked to candidate before final linking/merging"); + iEvent.put(std::move(linkedTable), "linkedTracksters"); + } + + // Full track-index / GSF-track-index lists per candidate + out->addColumn( + "nTrackIdxs", std::vector(table_size, 0), "Number of generalTracks associated with candidate"); + { + auto trackIdxsTable = std::make_unique(0, this->name_ + "TrackIdxs", false, false); + std::vector emptyIdx; + trackIdxsTable->addColumn("trackIndex", emptyIdx, "Index of associated generalTrack"); + trackIdxsTable->setDoc("Full list of generalTrack indices associated with the candidate"); + iEvent.put(std::move(trackIdxsTable), this->name_ + "TrackIdxs"); + } + + out->addColumn( + "nGsfTrackIdxs", std::vector(table_size, 0), "Number of GSFTracks associated with candidate"); + { + auto gsfTrackIdxsTable = std::make_unique(0, this->name_ + "GsfTrackIdxs", false, false); + std::vector emptyIdx; + gsfTrackIdxsTable->addColumn("trackIndex", emptyIdx, "Index of associated GSFTrack"); + gsfTrackIdxsTable->setDoc("Full list of GSFTrack indices associated with the candidate"); + iEvent.put(std::move(gsfTrackIdxsTable), this->name_ + "GsfTrackIdxs"); + } + + // isPU (sim candidates only); + if (hasPUInfo_) { + out->addColumn("isPU", + std::vector(table_size, -1), + "PU flag of the candidate's parent CaloParticle: 1 = pileup (non-zero event/bx), 0 = " + "otherwise, -1 = unresolved"); + } + if (out->nColumns() > 0) { out->setDoc(this->doc_); iEvent.put(std::move(out)); @@ -105,6 +180,56 @@ class TICLCandidateExtraTableProducer : public SimpleFlatTableProducerBasesize() : 0; + auto trackBoundaryTable = + std::make_unique(nAllTracks, "GeneralTrack", /*singleton*/ false, /*extension*/ true); + std::vector hgcal_x(nAllTracks), hgcal_y(nAllTracks), hgcal_z(nAllTracks); + std::vector hgcal_eta(nAllTracks), hgcal_phi(nAllTracks); + std::vector hgcal_px(nAllTracks), hgcal_py(nAllTracks), hgcal_pz(nAllTracks); + + if (allTracks_h.isValid()) { + const auto& allTracks = *allTracks_h; + for (size_t i = 0; i < allTracks.size(); ++i) { + const auto& track = allTracks[i]; + int iSide = int(track.eta() > 0); + const auto& fts = trajectoryStateTransform::outerFreeState(track, &bfield); + const auto& tsos = propagator.propagate(fts, firstDisk[iSide]->surface()); + if (tsos.isValid()) { + const auto& globalPos = tsos.globalPosition(); + const auto& globalMom = tsos.globalMomentum(); + hgcal_x[i] = globalPos.x(); + hgcal_y[i] = globalPos.y(); + hgcal_z[i] = globalPos.z(); + hgcal_eta[i] = globalPos.eta(); + hgcal_phi[i] = globalPos.phi(); + hgcal_px[i] = globalMom.x(); + hgcal_py[i] = globalMom.y(); + hgcal_pz[i] = globalMom.z(); + } else { + hgcal_x[i] = hgcal_y[i] = hgcal_z[i] = kInvalidBoundaryValue; + hgcal_eta[i] = hgcal_phi[i] = kInvalidBoundaryValue; + hgcal_px[i] = hgcal_py[i] = hgcal_pz[i] = kInvalidBoundaryValue; + } + } + } + + trackBoundaryTable->addColumn("hgcal_x", hgcal_x, "Track X position at HGCal boundary"); + trackBoundaryTable->addColumn("hgcal_y", hgcal_y, "Track Y position at HGCal boundary"); + trackBoundaryTable->addColumn("hgcal_z", hgcal_z, "Track Z position at HGCal boundary"); + trackBoundaryTable->addColumn("hgcal_eta", hgcal_eta, "Track eta at HGCal boundary"); + trackBoundaryTable->addColumn("hgcal_phi", hgcal_phi, "Track phi at HGCal boundary"); + trackBoundaryTable->addColumn("hgcal_px", hgcal_px, "Track Px at HGCal boundary"); + trackBoundaryTable->addColumn("hgcal_py", hgcal_py, "Track Py at HGCal boundary"); + trackBoundaryTable->addColumn("hgcal_pz", hgcal_pz, "Track Pz at HGCal boundary"); + iEvent.put(std::move(trackBoundaryTable), "trackBoundary"); + } + if (!prod.isValid() && this->skipNonExistingSrc_) { writeEmptyTables(iEvent, 0); return; @@ -124,6 +249,30 @@ class TICLCandidateExtraTableProducer : public SimpleFlatTableProducerBase>> linkedTracksters_h; + if (hasLinkedTracksters_) { + linkedTracksters_h = iEvent.getHandle(linkedTracksters_token_); + if (!linkedTracksters_h.isValid() && this->skipNonExistingSrc_) { + writeEmptyTables(iEvent, prod->size()); + return; + } + } + + // isPU (sim candidates only): resolved via the seedID/seedIndex of the candidate's first + // constituent Trackster, same CaloParticle-resolution logic as SimTracksterTableProducer, + edm::Handle> caloParticles_h; + edm::Handle>> cpToSCMap_h; + if (hasPUInfo_) { + caloParticles_h = iEvent.getHandle(caloParticles_token_); + cpToSCMap_h = iEvent.getHandle(caloParticleToSimClustersMap_token_); + if ((!caloParticles_h.isValid() || !cpToSCMap_h.isValid()) && this->skipNonExistingSrc_) { + writeEmptyTables(iEvent, prod->size()); + return; + } + } + const auto& candidates = *prod; const size_t table_size = candidates.size(); @@ -138,7 +287,25 @@ class TICLCandidateExtraTableProducer : public SimpleFlatTableProducerBase track_boundaryEta, track_boundaryPhi; std::vector track_boundaryPx, track_boundaryPy, track_boundaryPz; - for (const auto& cand : candidates) { + std::vector linkedCounts; + linkedCounts.reserve(table_size); + std::vector linkedTracksterIndices; + + std::vector trackIdxCounts; + trackIdxCounts.reserve(table_size); + std::vector allTrackIndices; + + std::vector gsfTrackIdxCounts; + gsfTrackIdxCounts.reserve(table_size); + std::vector allGsfTrackIndices; + + std::vector isPU; + if (hasPUInfo_) { + isPU.reserve(table_size); + } + + for (size_t i = 0; i < candidates.size(); ++i) { + const auto& cand = candidates[i]; const auto& children = cand.tracksters(); counts.push_back(children.size()); coltablesize += children.size(); @@ -147,13 +314,13 @@ class TICLCandidateExtraTableProducer : public SimpleFlatTableProducerBase 0); const auto& fts = trajectoryStateTransform::outerFreeState(track, &bfield); @@ -192,6 +359,70 @@ class TICLCandidateExtraTableProducer : public SimpleFlatTableProducerBasesize()) { + const auto& linked = (*linkedTracksters_h)[i]; + linkedCounts.push_back(linked.size()); + for (auto idx : linked) { + linkedTracksterIndices.push_back(idx); + } + } else { + linkedCounts.push_back(0); + } + + // Full track-index / GSF-track-index lists (all tracks, not just the first via + // trackPtr()/gsftrackPtr()). cand.trackPtrs()/gsfTrackPtrs() are Ptrs into whatever + // collection the candidate was built against -- key() gives the row index directly. + const auto& candTrackPtrs = cand.trackPtrs(); + trackIdxCounts.push_back(candTrackPtrs.size()); + for (const auto& tp : candTrackPtrs) { + allTrackIndices.push_back(tp.key()); + } + + const auto& candGsfTrackPtrs = cand.gsfTrackPtrs(); + gsfTrackIdxCounts.push_back(candGsfTrackPtrs.size()); + for (const auto& gp : candGsfTrackPtrs) { + allGsfTrackIndices.push_back(gp.key()); + } + + // isPU: resolve the parent CaloParticle from the first constituent Trackster's + // seedID/seedIndex, then read its g4Track eventId/bunchCrossing + if (hasPUInfo_) { + int puFlag = -1; + if (caloParticles_h.isValid() && cpToSCMap_h.isValid() && !children.empty()) { + const auto& firstTrackster = tracksters[children[0].key()]; + const auto& caloParticles = *caloParticles_h; + const auto& cpToSCMap = *cpToSCMap_h; + + const CaloParticle* cp = nullptr; + if (firstTrackster.seedID() == caloParticles_h.id()) { + const auto seedIdx = firstTrackster.seedIndex(); + if (seedIdx >= 0 && static_cast(seedIdx) < caloParticles.size()) { + cp = &caloParticles[seedIdx]; + } + } else { + const auto seedIdx = firstTrackster.seedIndex(); + for (const auto& [cpIdx, scVec] : cpToSCMap) { + if (seedIdx >= 0 && + std::find(scVec.begin(), scVec.end(), static_cast(seedIdx)) != scVec.end()) { + if (static_cast(cpIdx) < caloParticles.size()) { + cp = &caloParticles[cpIdx]; + } + break; + } + } + } + + if (cp && !cp->g4Tracks().empty()) { + const auto& simTrack = cp->g4Tracks()[0]; + puFlag = (simTrack.eventId().event() != 0 || simTrack.eventId().bunchCrossing() != 0) ? 1 : 0; + } + } + isPU.push_back(puFlag); + } } for (const auto& coltable : coltables_) { @@ -225,6 +456,49 @@ class TICLCandidateExtraTableProducer : public SimpleFlatTableProducerBaseaddColumn( + "nLinkedTracksters", linkedCounts, "Number of Tracksters linked to candidate before final linking/merging"); + + auto linkedTable = + std::make_unique(linkedTracksterIndices.size(), "linkedTracksters", false, false); + linkedTable->addColumn("tracksterIndex", linkedTracksterIndices, "Index of linked Trackster"); + linkedTable->setDoc("Tracksters linked to candidate before final linking/merging"); + iEvent.put(std::move(linkedTable), "linkedTracksters"); + } + + // Full track-index / GSF-track-index lists: count columns on the main table + flattened + // sub-tables, always emitted (see constructor comment). Product names are prefixed with + // this->name_ so the reco and sim module instances never collide. + out->addColumn("nTrackIdxs", trackIdxCounts, "Number of generalTracks associated with candidate"); + { + auto trackIdxsTable = + std::make_unique(allTrackIndices.size(), this->name_ + "TrackIdxs", false, false); + trackIdxsTable->addColumn("trackIndex", allTrackIndices, "Index of associated generalTrack"); + trackIdxsTable->setDoc("Full list of generalTrack indices associated with the candidate"); + iEvent.put(std::move(trackIdxsTable), this->name_ + "TrackIdxs"); + } + + out->addColumn("nGsfTrackIdxs", gsfTrackIdxCounts, "Number of GSFTracks associated with candidate"); + { + auto gsfTrackIdxsTable = + std::make_unique(allGsfTrackIndices.size(), this->name_ + "GsfTrackIdxs", false, false); + gsfTrackIdxsTable->addColumn("trackIndex", allGsfTrackIndices, "Index of associated GSFTrack"); + gsfTrackIdxsTable->setDoc("Full list of GSFTrack indices associated with the candidate"); + iEvent.put(std::move(gsfTrackIdxsTable), this->name_ + "GsfTrackIdxs"); + } + + // isPU output: one scalar per candidate, directly on the extension table. + // Only emitted for instances configured with caloParticles + caloParticleToSimClustersMap (sim). + if (hasPUInfo_) { + out->addColumn("isPU", + isPU, + "PU flag of the candidate's parent CaloParticle: 1 = pileup (non-zero event/bx), 0 = " + "otherwise, -1 = unresolved"); + } + if (out->nColumns() > 0) { out->setDoc(this->doc_); iEvent.put(std::move(out)); @@ -242,8 +516,25 @@ class TICLCandidateExtraTableProducer : public SimpleFlatTableProducerBase("tracksters", edm::InputTag("ticlTrackstersCLUE3DHigh")); desc.add("tracks", edm::InputTag("generalTracks")); + desc.addOptional("linkedTracksters") + ->setComment("Pre-merge/pre-linking Trackster links (reco TICLCandidates only, e.g. "); + desc.addOptional("caloParticles") + ->setComment( + "CaloParticle collection (sim candidates only), used together with " + "caloParticleToSimClustersMap to compute the isPU column. Omit for reco candidates."); + desc.addOptional("caloParticleToSimClustersMap") + ->setComment( + "CaloParticle-to-SimCluster index map (sim candidates only), e.g. the unlabeled product from " + "SimTrackstersProducer (same module as the simTracksters src). Required together with " + "caloParticles to compute the isPU column. Omit for reco candidates."); desc.add("detector", "HGCAL"); desc.add("propagator", "PropagatorWithMaterial"); + desc.add("produceGeneralTrackBoundary", false) + ->setComment( + "If true, also emit an extension of the GeneralTrack table with HGCal-boundary " + "propagated position/momentum for every track in the tracks collection (unconditional, " + "not restricted to candidate-linked tracks), mirroring TICLDumper's track_hgcal_* branches. " + "Enable on exactly one module instance to avoid duplicate GeneralTrack extensions."); edm::ParameterSetDescription coltable; coltable.add("name", "hltTiclCandidate"); @@ -272,6 +563,12 @@ class TICLCandidateExtraTableProducer : public SimpleFlatTableProducerBase> tracksters_token_; const edm::EDGetTokenT> tracks_token_; + const bool hasLinkedTracksters_; + const edm::EDGetTokenT>> linkedTracksters_token_; + const bool hasPUInfo_; + const edm::EDGetTokenT> caloParticles_token_; + const edm::EDGetTokenT>> caloParticleToSimClustersMap_token_; + const bool produceGeneralTrackBoundary_; const std::string detector_; const std::string propName_; const edm::ESGetToken geometry_token_; diff --git a/DPGAnalysis/HGCalNanoAOD/plugins/TrackMuonInfoProducer.cc b/DPGAnalysis/HGCalNanoAOD/plugins/TrackMuonInfoProducer.cc new file mode 100644 index 0000000000000..a9c5754b30200 --- /dev/null +++ b/DPGAnalysis/HGCalNanoAOD/plugins/TrackMuonInfoProducer.cc @@ -0,0 +1,117 @@ +// TrackMuonInfoProducer.cc +#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 "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/TrackReco/interface/TrackFwd.h" +#include "DataFormats/MuonReco/interface/Muon.h" +#include "DataFormats/MuonReco/interface/MuonFwd.h" + +#include "DataFormats/Common/interface/ValueMap.h" +#include "FWCore/Utilities/interface/EDGetToken.h" + +#include "RecoParticleFlow/PFProducer/interface/PFMuonAlgo.h" + +class TrackMuonInfoProducer : public edm::global::EDProducer<> { +public: + explicit TrackMuonInfoProducer(const edm::ParameterSet& cfg) + : tracksToken_(consumes>(cfg.getParameter("tracks"))), + muonsToken_(consumes(cfg.getParameter("muons"))) { + produces>("isMuon"); + produces>("isTrackerMuon"); + produces>("muonDtHits"); + produces>("muonCscHits"); + produces>("muonType"); + } + + void produce(edm::StreamID, edm::Event& evt, const edm::EventSetup&) const override { + edm::Handle> tracksH; + evt.getByToken(tracksToken_, tracksH); + edm::Handle muonsH; + evt.getByToken(muonsToken_, muonsH); + + const auto& tracks = *tracksH; + const auto& muons = *muonsH; + + std::vector isMuonV; + isMuonV.reserve(tracks.size()); + std::vector isTrackerMuonV; + isTrackerMuonV.reserve(tracks.size()); + std::vector muonDtV; + muonDtV.reserve(tracks.size()); + std::vector muonCscV; + muonCscV.reserve(tracks.size()); + std::vector muonTypeV; + muonTypeV.reserve(tracks.size()); + + for (size_t i = 0; i < tracks.size(); ++i) { + reco::TrackRef tref(tracksH, i); + int muId = -1; + for (size_t im = 0; im < muons.size(); ++im) { + if (muons[im].track().isNonnull() && muons[im].track() == tref) { + muId = (int)im; + break; + } + } + if (muId == -1) { + isMuonV.push_back(-1); + isTrackerMuonV.push_back(-1); + muonDtV.push_back(-1); + muonCscV.push_back(-1); + muonTypeV.push_back(-1); + } else { + reco::MuonRef mref(muonsH, muId); + bool ismu = PFMuonAlgo::isMuon(mref); + isMuonV.push_back(ismu ? 1 : 0); + isTrackerMuonV.push_back((*muonsH)[muId].isTrackerMuon() ? 1 : 0); + int dtHits = 0, cscHits = 0; + if (mref->standAloneMuon().isNonnull()) { + auto st = mref->standAloneMuon(); + dtHits = st->hitPattern().numberOfValidMuonDTHits(); + cscHits = st->hitPattern().numberOfValidMuonCSCHits(); + } + muonDtV.push_back(dtHits); + muonCscV.push_back(cscHits); + muonTypeV.push_back(mref->type()); + } + } + + auto vm_isMuon = std::make_unique>(); + edm::ValueMap::Filler f_isMuon(*vm_isMuon); + f_isMuon.insert(tracksH, isMuonV.begin(), isMuonV.end()); + f_isMuon.fill(); + evt.put(std::move(vm_isMuon), "isMuon"); + + auto vm_isTracker = std::make_unique>(); + edm::ValueMap::Filler f_isTracker(*vm_isTracker); + f_isTracker.insert(tracksH, isTrackerMuonV.begin(), isTrackerMuonV.end()); + f_isTracker.fill(); + evt.put(std::move(vm_isTracker), "isTrackerMuon"); + + auto vm_dt = std::make_unique>(); + edm::ValueMap::Filler f_dt(*vm_dt); + f_dt.insert(tracksH, muonDtV.begin(), muonDtV.end()); + f_dt.fill(); + evt.put(std::move(vm_dt), "muonDtHits"); + + auto vm_csc = std::make_unique>(); + edm::ValueMap::Filler f_csc(*vm_csc); + f_csc.insert(tracksH, muonCscV.begin(), muonCscV.end()); + f_csc.fill(); + evt.put(std::move(vm_csc), "muonCscHits"); + + auto vm_type = std::make_unique>(); + edm::ValueMap::Filler f_type(*vm_type); + f_type.insert(tracksH, muonTypeV.begin(), muonTypeV.end()); + f_type.fill(); + evt.put(std::move(vm_type), "muonType"); + } + +private: + const edm::EDGetTokenT> tracksToken_; + const edm::EDGetTokenT muonsToken_; +}; + +DEFINE_FWK_MODULE(TrackMuonInfoProducer); diff --git a/DPGAnalysis/HGCalNanoAOD/python/HGCalNanoAOD_cff.py b/DPGAnalysis/HGCalNanoAOD/python/HGCalNanoAOD_cff.py index e39f799a7caac..9435a73fabe7c 100644 --- a/DPGAnalysis/HGCalNanoAOD/python/HGCalNanoAOD_cff.py +++ b/DPGAnalysis/HGCalNanoAOD/python/HGCalNanoAOD_cff.py @@ -5,6 +5,9 @@ from DPGAnalysis.HGCalNanoAOD.hgcalTICLCandidates_cfi import * from DPGAnalysis.HGCalNanoAOD.hgcalTICLSuperClusters_cfi import * from DPGAnalysis.HGCalNanoAOD.hgcalLayerClusters_cfi import * +from DPGAnalysis.HGCalNanoAOD.hgcalGeneralTracks_cfi import * +from DPGAnalysis.HGCalNanoAOD.hgcalGSFTracks_cfi import * +from DPGAnalysis.HGCalNanoAOD.hgcalGen_cfi import * ###################################### # Offline HGCAL NanoAOD Tables @@ -14,6 +17,8 @@ hgcalTrackstersTableSequence + ticlCandidateTable + ticlCandidateExtraTable + + hgcalGeneralTracksTableSequence + + hgcalGSFTracksTableSequence ) # Store additional validation objects @@ -23,6 +28,7 @@ + ticlSimCandidateTable + ticlSimCandidateExtraTable + hgcalLayerClustersTableSequence + + hgcalGenSequence ) ###################################### diff --git a/DPGAnalysis/HGCalNanoAOD/python/hgcalGSFTracks_cfi.py b/DPGAnalysis/HGCalNanoAOD/python/hgcalGSFTracks_cfi.py new file mode 100644 index 0000000000000..e2c42b23bd2d2 --- /dev/null +++ b/DPGAnalysis/HGCalNanoAOD/python/hgcalGSFTracks_cfi.py @@ -0,0 +1,49 @@ +import FWCore.ParameterSet.Config as cms +from PhysicsTools.NanoAOD.common_cff import * + +# GSF Tracks table for NanoAOD +gsfTracksTable = cms.EDProducer( + "SimpleGsfTrackFlatTableProducer", + skipNonExistingSrc=cms.bool(True), + src=cms.InputTag("electronGsfTracks"), + cut=cms.string(""), + name=cms.string("GSFTrack"), + doc=cms.string("GSF tracks"), + extension=cms.bool(False), + variables=cms.PSet( + # Basic kinematics + pt=Var("pt()", "float", doc="GSF track p_T [GeV]"), + p=Var("p()", "float", doc="GSF track momentum magnitude [GeV]"), + eta=Var("eta()", "float", doc="GSF track pseudorapidity"), + phi=Var("phi()", "float", doc="GSF track phi angle [rad]"), + charge=Var("charge()", "int", doc="GSF track charge"), + + # Position + vx=Var("vx()", "float", doc="GSF track vertex x [cm]"), + vy=Var("vy()", "float", doc="GSF track vertex y [cm]"), + vz=Var("vz()", "float", doc="GSF track vertex z [cm]"), + + nhits=Var("recHitsSize()", "int", doc="GSF track nHits"), + ptMode=Var("ptMode()", "float", doc="GSF track p_T Mode [GeV]"), + ptModeError=Var("ptModeError()", "float", doc="GSF track p_T Mode Error"), + pxMode=Var("pxMode()", "float", doc="GSF track px Mode [GeV]"), + pyMode=Var("pyMode()", "float", doc="GSF track py Mode [GeV]"), + pzMode=Var("pzMode()", "float", doc="GSF track pz Mode [GeV]"), + pMode=Var("pMode()", "float", doc="GSF track momentum Mode [GeV]"), + etaMode=Var("etaMode()", "float", doc="GSF track Mode pseudorapidity"), + etaModeError=Var("etaModeError()", "float", doc="GSF track Mode pseudorapidity Error"), + phiMode=Var("phiMode()", "float", doc="GSF track phi angle Mode [rad]"), + phiModeError=Var("phiModeError()", "float", doc="GSF track phi angle Mode Error"), + chargeMode=Var("chargeMode()", "int", doc="GSF track charge Mode"), + lambdaMode=Var("lambdaMode()", "float", doc="GSF track lambda Mode"), + lambdaModeError=Var("lambdaModeError()", "float", doc="GSF track lambda Mode Error"), + thetaMode=Var("thetaMode()", "float", doc="GSF track theta Mode"), + thetaModeError=Var("thetaModeError()", "float", doc="GSF track theta Mode Error"), + qoverpMode=Var("qoverpMode()", "float", doc="GSF track qoverp Mode"), + qoverpModeError=Var("qoverpModeError()", "float", doc="GSF track qoverp Mode Error"), + + ), +) + +# Sequence for gsf tracks +hgcalGSFTracksTableSequence = cms.Sequence(gsfTracksTable) diff --git a/DPGAnalysis/HGCalNanoAOD/python/hgcalGen_cfi.py b/DPGAnalysis/HGCalNanoAOD/python/hgcalGen_cfi.py new file mode 100644 index 0000000000000..11ddb2c34b1a7 --- /dev/null +++ b/DPGAnalysis/HGCalNanoAOD/python/hgcalGen_cfi.py @@ -0,0 +1,28 @@ +import FWCore.ParameterSet.Config as cms +from PhysicsTools.NanoAOD.common_cff import * +from PhysicsTools.NanoAOD.nano_cff import nanoMetadata + +# GenParticles table +hgcalGenPartTable = cms.EDProducer( + "SimpleGenParticleFlatTableProducer", + skipNonExistingSrc=cms.bool(True), + src=cms.InputTag("genParticles"), + cut=cms.string("status==1 && pt > 0.5"), # stable gen particles with pT > 0.5 + name=cms.string("HGCalGenPart"), + doc=cms.string("Gen particles information"), + singleton=cms.bool(False), + variables=cms.PSet( + pt=Var("pt", "float", doc="gen particle pT"), + eta=Var("eta", "float", doc="gen particle eta"), + phi=Var("phi", "float", doc="gen particle phi"), + mass=Var("mass", "float", doc="gen particle mass"), + pdgId=Var("pdgId", "int", doc="gen particle PDG ID"), + status=Var("status", "int", doc="gen particle status"), + energy=Var("energy", "float", doc="gen particle energy"), + vx=Var("vx", "float", doc="gen particle vx"), + vy=Var("vy", "float", doc="gen particle vy"), + vz=Var("vz", "float", doc="gen particle vz"), + ), +) + +hgcalGenSequence = cms.Sequence(hgcalGenPartTable) diff --git a/DPGAnalysis/HGCalNanoAOD/python/hgcalGeneralTracks_cfi.py b/DPGAnalysis/HGCalNanoAOD/python/hgcalGeneralTracks_cfi.py new file mode 100644 index 0000000000000..59dfe4cb49874 --- /dev/null +++ b/DPGAnalysis/HGCalNanoAOD/python/hgcalGeneralTracks_cfi.py @@ -0,0 +1,53 @@ +import FWCore.ParameterSet.Config as cms +from PhysicsTools.NanoAOD.common_cff import * + +trackMuonInfo = cms.EDProducer("TrackMuonInfoProducer", + tracks = cms.InputTag("generalTracks"), + muons = cms.InputTag("muons") +) + +# General Tracks table for NanoAOD +generalTracksTable = cms.EDProducer( + "SimpleTrackFlatTableProducer", + skipNonExistingSrc=cms.bool(True), + src=cms.InputTag("generalTracks"), + cut=cms.string(""), + name=cms.string("GeneralTrack"), + doc=cms.string("General reconstructed tracks"), + extension=cms.bool(False), + variables=cms.PSet( + # Basic kinematics + pt=Var("pt()", "float", doc="track p_T [GeV]"), + p=Var("p()", "float", doc="track momentum magnitude [GeV]"), + eta=Var("eta()", "float", doc="track pseudorapidity"), + phi=Var("phi()", "float", doc="track phi angle [rad]"), + charge=Var("charge()", "int", doc="track charge"), + trackLambda=Var("lambda()", "float", doc="track lambda"), + + # Position + vx=Var("vx()", "float", doc="track vertex x [cm]"), + vy=Var("vy()", "float", doc="track vertex y [cm]"), + vz=Var("vz()", "float", doc="track vertex z [cm]"), + + # Quality metrics + nhits=Var("numberOfValidHits()", "uint16", doc="number of valid hits"), + missingOuterHits=Var("missingOuterHits()", "uint8", doc="number of missing outer hits"), + + # Error parameters + ptErr=Var("ptError()", "float", doc="track p_T error [GeV]"), + etaErr=Var("etaError()", "float", doc="track eta error"), + phiErr=Var("phiError()", "float", doc="track phi error"), + lambdaErr=Var("lambdaError()", "float", doc="track lambda error"), + qoverpErr=Var("qoverpError()", "float", doc="q/p error"), + ), + externalVariables = cms.PSet( + isMuon = ExtVar(cms.InputTag("trackMuonInfo", "isMuon"),"int", doc="1 if PFMuonAlgo::isMuon, 0 otherwise, -1 if no muon"), + isTrackerMuon = ExtVar(cms.InputTag("trackMuonInfo", "isTrackerMuon"),"int", doc="1 if tracker muon, -1 if no muon"), + muon_dt_hits = ExtVar(cms.InputTag("trackMuonInfo", "muonDtHits"),"int", doc="ST muon DT hits (-1 if no muon)"), + muon_csc_hits = ExtVar(cms.InputTag("trackMuonInfo", "muonCscHits"), "int", doc="ST muon CSC hits (-1 if no muon)"), + muon_type = ExtVar(cms.InputTag("trackMuonInfo", "muonType"), "int", doc="muon type (-1 if no muon)"), + ) +) + +# Sequence for general tracks +hgcalGeneralTracksTableSequence = cms.Sequence(trackMuonInfo + generalTracksTable) diff --git a/DPGAnalysis/HGCalNanoAOD/python/hgcalTICLCandidates_cfi.py b/DPGAnalysis/HGCalNanoAOD/python/hgcalTICLCandidates_cfi.py index 9aee5bdcd5335..16df4c7ac13d1 100644 --- a/DPGAnalysis/HGCalNanoAOD/python/hgcalTICLCandidates_cfi.py +++ b/DPGAnalysis/HGCalNanoAOD/python/hgcalTICLCandidates_cfi.py @@ -49,6 +49,7 @@ ), ) + # SimTICLCandidates for validation ticlSimCandidateTable = cms.EDProducer( "TICLCandidateTableProducer", @@ -91,18 +92,24 @@ MTDtimeError=Var("MTDtimeError", "float", doc="Trackster associated MTD time error, meaningful only for offline reconstruction"), trackIdx=Var("trackPtr().key", "int", - doc="Index of generalTrack associated with TICLCandidate") + doc="Index of generalTrack associated with SimTICLCandidate"), + gsftrackIdx=Var("gsftrackPtr().key", "int", + doc="Index of GSFTrack associated with SimTICLCandidate"), + ), ) # Extra tables for linking tracksters to candidates ticlCandidateExtraTable = cms.EDProducer( "TICLCandidateExtraTableProducer", + skipNonExistingSrc=cms.bool(True), src=cms.InputTag("ticlCandidate"), name=cms.string("TICLCandidates"), doc=cms.string("TICLCandidates extra table with linked Tracksters"), tracksters=cms.InputTag("ticlCandidate"), tracks=cms.InputTag("generalTracks"), + linkedTracksters=cms.InputTag("ticlCandidate", "linkedTracksters"), + produceGeneralTrackBoundary=cms.bool(True), detector=cms.string("HGCAL"), propagator=cms.string("PropagatorWithMaterial"), collectionVariables=cms.PSet( @@ -118,11 +125,14 @@ ticlSimCandidateExtraTable = cms.EDProducer( "TICLCandidateExtraTableProducer", + skipNonExistingSrc=cms.bool(True), src=cms.InputTag("ticlSimTracksters"), name=cms.string("SimTICLCandidates"), doc=cms.string("TICLCandidates extra table with linked Tracksters"), tracksters=cms.InputTag("ticlSimTracksters"), tracks=cms.InputTag("generalTracks"), + caloParticles=cms.InputTag("mix", "MergedCaloTruth"), + caloParticleToSimClustersMap=cms.InputTag("ticlSimTracksters"), detector=cms.string("HGCAL"), propagator=cms.string("PropagatorWithMaterial"), collectionVariables=cms.PSet( @@ -135,4 +145,3 @@ ), ), ) - diff --git a/RecoHGCal/TICL/interface/TICLInterpretationAlgoBase.h b/RecoHGCal/TICL/interface/TICLInterpretationAlgoBase.h index 7dfbec58adad5..1d91b0d5f3bf4 100644 --- a/RecoHGCal/TICL/interface/TICLInterpretationAlgoBase.h +++ b/RecoHGCal/TICL/interface/TICLInterpretationAlgoBase.h @@ -91,7 +91,8 @@ namespace ticl { edm::Handle inputTiming_h, std::vector& resultTracksters, std::vector& resultCandidate, - std::vector& maskedTracksters) = 0; + std::vector& maskedTracksters, + std::vector>& linkedResultTracksters) = 0; virtual void initialize(const HGCalDDDConstants* hgcons, const ticlgeom::Tools rhtools, diff --git a/RecoHGCal/TICL/plugins/GNNInterpretationAlgo.cc b/RecoHGCal/TICL/plugins/GNNInterpretationAlgo.cc index bae60ee643ea3..1e752ce40f5b8 100644 --- a/RecoHGCal/TICL/plugins/GNNInterpretationAlgo.cc +++ b/RecoHGCal/TICL/plugins/GNNInterpretationAlgo.cc @@ -332,7 +332,8 @@ void GNNInterpretationAlgo::makeCandidates(const Inputs& input, edm::Handle inputTiming_h, std::vector& resultTracksters, std::vector& resultCandidate, - std::vector& maskedTracksters) { + std::vector& maskedTracksters, + std::vector>& linkedResultTracksters) { const auto& tracks = *input.tracksHandle; const auto& maskTracks = input.maskedTracks; const auto& tracksters = input.tracksters; @@ -525,6 +526,7 @@ void GNNInterpretationAlgo::makeCandidates(const Inputs& input, } } // Build output tracksters + linkedResultTracksters.reserve(input.tracksters.size()); for (unsigned trkId = 0; trkId < trackToTracksters.size(); ++trkId) { if (trackToTracksters[trkId].empty()) @@ -545,13 +547,15 @@ void GNNInterpretationAlgo::makeCandidates(const Inputs& input, 1.f); resultTracksters.push_back(std::move(merged)); + linkedResultTracksters.push_back(trackToTracksters[trkId]); } } // Add unlinked tracksters - for (unsigned i = 0; i < tracksters.size(); ++i) { - if (tracksterAvailable[i]) - resultTracksters.push_back(tracksters[i]); + for (auto iTrackster = 0u; iTrackster < input.tracksters.size(); iTrackster++) { + if (tracksterAvailable[iTrackster]) + resultTracksters.push_back(tracksters[iTrackster]); + linkedResultTracksters.push_back({iTrackster}); } } diff --git a/RecoHGCal/TICL/plugins/GNNInterpretationAlgo.h b/RecoHGCal/TICL/plugins/GNNInterpretationAlgo.h index 4fa73fb2c5b64..33f5e769967da 100644 --- a/RecoHGCal/TICL/plugins/GNNInterpretationAlgo.h +++ b/RecoHGCal/TICL/plugins/GNNInterpretationAlgo.h @@ -44,7 +44,8 @@ namespace ticl { edm::Handle inputTiming_h, std::vector &resultTracksters, std::vector &resultCandidate, - std::vector &maskedTracksters) override; + std::vector &maskedTracksters, + std::vector> &linkedResultTracksters) override; void initialize(const HGCalDDDConstants *hgcons, const ticlgeom::Tools rhtools, diff --git a/RecoHGCal/TICL/plugins/GeneralInterpretationAlgo.cc b/RecoHGCal/TICL/plugins/GeneralInterpretationAlgo.cc index 13a0f79b0b91a..3f3d2894df91e 100644 --- a/RecoHGCal/TICL/plugins/GeneralInterpretationAlgo.cc +++ b/RecoHGCal/TICL/plugins/GeneralInterpretationAlgo.cc @@ -197,7 +197,8 @@ void GeneralInterpretationAlgo::makeCandidates(const Inputs &input, edm::Handle inputTiming_h, std::vector &resultTracksters, std::vector &resultCandidate, - std::vector &maskedTracksters) { + std::vector &maskedTracksters, + std::vector> &linkedResultTracksters) { bool useMTDTiming = inputTiming_h.isValid(); const auto tkH = input.tracksHandle; const auto maskTracks = input.maskedTracks; @@ -372,13 +373,14 @@ void GeneralInterpretationAlgo::makeCandidates(const Inputs &input, } trackstersInTrackIndices[i] = chargedCandidate; } - + linkedResultTracksters.reserve(input.tracksters.size()); for (size_t iTrack = 0; iTrack < trackstersInTrackIndices.size(); iTrack++) { if (!trackstersInTrackIndices[iTrack].empty()) { if (trackstersInTrackIndices[iTrack].size() == 1) { auto tracksterId = trackstersInTrackIndices[iTrack][0]; resultCandidate[iTrack] = resultTracksters.size(); resultTracksters.push_back(input.tracksters[tracksterId]); + linkedResultTracksters.push_back(trackstersInTrackIndices[iTrack]); } else { // in this case mergeTracksters() clears the pid probabilities and the regressed energy is not set // TODO: fix probabilities when CNN will be splitted @@ -398,12 +400,14 @@ void GeneralInterpretationAlgo::makeCandidates(const Inputs &input, else resultTracksters.back().setIdProbability(ticl::Trackster::ParticleType::electron, 1.f); } + linkedResultTracksters.push_back(trackstersInTrackIndices[iTrack]); } } - for (size_t iTrackster = 0; iTrackster < input.tracksters.size(); iTrackster++) { + for (auto iTrackster = 0u; iTrackster < input.tracksters.size(); iTrackster++) { if (chargedMask[iTrackster]) { resultTracksters.push_back(input.tracksters[iTrackster]); + linkedResultTracksters.push_back({iTrackster}); } } }; diff --git a/RecoHGCal/TICL/plugins/GeneralInterpretationAlgo.h b/RecoHGCal/TICL/plugins/GeneralInterpretationAlgo.h index ba805139131d1..20c853daecddb 100644 --- a/RecoHGCal/TICL/plugins/GeneralInterpretationAlgo.h +++ b/RecoHGCal/TICL/plugins/GeneralInterpretationAlgo.h @@ -21,7 +21,8 @@ namespace ticl { edm::Handle inputTiming_h, std::vector &resultTracksters, std::vector &resultCandidate, - std::vector &maskedTracksters) override; + std::vector &maskedTracksters, + std::vector> &linkedResultTracksters) override; void initialize(const HGCalDDDConstants *hgcons, const ticlgeom::Tools rhtools, diff --git a/RecoHGCal/TICL/plugins/MuonInterpretationAlgo.cc b/RecoHGCal/TICL/plugins/MuonInterpretationAlgo.cc index 173d7dc27a771..6bd2577be6ad1 100644 --- a/RecoHGCal/TICL/plugins/MuonInterpretationAlgo.cc +++ b/RecoHGCal/TICL/plugins/MuonInterpretationAlgo.cc @@ -44,7 +44,8 @@ void MuonInterpretationAlgo::makeCandidates(const Inputs &input, edm::Handle /*inputTiming_h*/, std::vector &resultTracksters, std::vector &resultCandidate, - std::vector &maskedTracksters) { + std::vector &maskedTracksters, + std::vector> &linkedResultTracksters) { const auto &tracks = *input.tracksHandle; const auto &maskTracks = input.maskedTracks; const auto &tracksters = input.tracksters; diff --git a/RecoHGCal/TICL/plugins/MuonInterpretationAlgo.h b/RecoHGCal/TICL/plugins/MuonInterpretationAlgo.h index 5ab3542060c33..5281a33753d55 100644 --- a/RecoHGCal/TICL/plugins/MuonInterpretationAlgo.h +++ b/RecoHGCal/TICL/plugins/MuonInterpretationAlgo.h @@ -43,7 +43,8 @@ namespace ticl { edm::Handle inputTiming_h, std::vector &resultTracksters, std::vector &resultCandidate, - std::vector &maskedTracksters) override; + std::vector &maskedTracksters, + std::vector> &linkedResultTracksters) override; void initialize(const HGCalDDDConstants *hgcons, const ticlgeom::Tools rhtools, diff --git a/RecoHGCal/TICL/plugins/TICLCandidateProducer.cc b/RecoHGCal/TICL/plugins/TICLCandidateProducer.cc index a27ea74c17aa9..0dcfcbb106327 100644 --- a/RecoHGCal/TICL/plugins/TICLCandidateProducer.cc +++ b/RecoHGCal/TICL/plugins/TICLCandidateProducer.cc @@ -209,6 +209,7 @@ TICLCandidateProducer::TICLCandidateProducer(const edm::ParameterSet &ps, const } produces>(); + produces>>("linkedTracksters"); // New trackster collection after linking produces>(); @@ -376,7 +377,7 @@ void TICLCandidateProducer::produce(edm::Event &evt, const edm::EventSetup &es) // muon), or a rejection (kMuonRejected: the trajectory points to a shower). std::vector maskedInputTracksters(generalTrackstersSpan.size(), false); muonInterpretationAlgo_->makeCandidates( - muonInput, inputTiming_h, *resultTracksters, muonInTrackIndices, maskedInputTracksters); + muonInput, inputTiming_h, *resultTracksters, muonInTrackIndices, maskedInputTracksters, *linkedResultTracksters); // A track the muon pass rejected is not a muon: route it back to the general pass so // it is reconstructed there (and no muon candidate is built for it below). @@ -396,7 +397,7 @@ void TICLCandidateProducer::produce(edm::Event &evt, const edm::EventSetup &es) tracks_h, generalTrackMask); generalInterpretationAlgo_->makeCandidates( - input, inputTiming_h, *resultTracksters, trackstersInTrackIndices, maskedInputTracksters); + input, inputTiming_h, *resultTracksters, trackstersInTrackIndices, maskedInputTracksters, *linkedResultTracksters); assignPCAtoTracksters(*resultTracksters, layerClusters, @@ -411,6 +412,7 @@ void TICLCandidateProducer::produce(edm::Event &evt, const edm::EventSetup &es) std::vector maskTracksters(resultTracksters->size(), true); edm::OrphanHandle> resultTracksters_h = evt.put(std::move(resultTracksters)); + auto linkedTracksters = std::make_unique>>(); // Muon candidates: energy from the track momentum (pdgId 13), attaching the MIP // trackster the muon pass associated (if any) and masking it so it is not re-emitted. @@ -440,6 +442,7 @@ void TICLCandidateProducer::produce(edm::Event &evt, const edm::EventSetup &es) if (tracksterId != -1 and !maskTracksters.empty()) { auto tracksterPtr = edm::Ptr(resultTracksters_h, tracksterId); TICLCandidate chargedCandidate(trackPtr, tracksterPtr); + linkedTracksters->push_back((*linkedResultTracksters)[tracksterId]); resultCandidates->push_back(chargedCandidate); maskTracksters[tracksterId] = false; } @@ -452,6 +455,7 @@ void TICLCandidateProducer::produce(edm::Event &evt, const edm::EventSetup &es) edm::Ptr tracksterPtr(resultTracksters_h, iTrackster); edm::Ptr trackPtr; TICLCandidate neutralCandidate(trackPtr, tracksterPtr); + linkedTracksters->push_back((*linkedResultTracksters)[iTrackster]); resultCandidates->push_back(neutralCandidate); } } @@ -511,6 +515,7 @@ void TICLCandidateProducer::produce(edm::Event &evt, const edm::EventSetup &es) assignTimeToCandidates(*resultCandidates, tracks_h, inputTimingView, getPathLength); evt.put(std::move(resultCandidates)); + evt.put(std::move(linkedTracksters), "linkedTracksters"); } template