diff --git a/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py b/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py index 867c8a52e9b09..cd378f8f543b4 100644 --- a/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py +++ b/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py @@ -35,6 +35,7 @@ 'keep *_ticlSimTracksters_*_*', 'keep *_ticlSimTICLCandidates_*_*', 'keep *_ticlSimTrackstersFromCP_*_*', + 'keep *_SimTau*_*_*' ) ) TICL_FEVT.outputCommands.extend(TICL_RECO.outputCommands) diff --git a/SimCalorimetry/HGCalAssociatorProducers/plugins/SimTauCPLinkProducer.cc b/SimCalorimetry/HGCalAssociatorProducers/plugins/SimTauCPLinkProducer.cc new file mode 100644 index 0000000000000..356a6b49d456a --- /dev/null +++ b/SimCalorimetry/HGCalAssociatorProducers/plugins/SimTauCPLinkProducer.cc @@ -0,0 +1,133 @@ +// -*- C++ -*- +// +// +// Authors: Marco Rovere, Andreas Gruber +// Created: Mon, 16 Oct 2023 14:24:35 GMT +// +// + +// system include files +#include + +// 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/StreamID.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include "DataFormats/HGCalReco/interface/TICLCandidate.h" +#include "SimDataFormats/CaloAnalysis/interface/CaloParticleFwd.h" +#include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h" +#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h" +#include "SimDataFormats/CaloAnalysis/interface/SimTauCPLink.h" + +class SimTauProducer : public edm::global::EDProducer<> { +public: + explicit SimTauProducer(const edm::ParameterSet&); + ~SimTauProducer() override = default; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + +private: + void produce(edm::StreamID, edm::Event&, edm::EventSetup const&) const override; + + void buildSimTau(SimTauCPLink&, + uint8_t, + int, + const reco::GenParticle&, + int, + edm::Handle>, + const std::vector&) const; + // ----------member data --------------------------- + const edm::EDGetTokenT> caloParticles_token_; + const edm::EDGetTokenT> genParticles_token_; + const edm::EDGetTokenT> genBarcodes_token_; +}; + +SimTauProducer::SimTauProducer(const edm::ParameterSet& iConfig) + : caloParticles_token_(consumes>(iConfig.getParameter("caloParticles"))), + genParticles_token_(consumes(iConfig.getParameter("genParticles"))), + genBarcodes_token_(consumes>(iConfig.getParameter("genBarcodes"))) { + produces>(); +} + +void SimTauProducer::buildSimTau(SimTauCPLink& t, + uint8_t generation, + int resonance_idx, + const reco::GenParticle& gen_particle, + int gen_particle_key, + edm::Handle> calo_particle_h, + const std::vector& gen_particle_barcodes) const { + const auto& caloPartVec = *calo_particle_h; + auto& daughters = gen_particle.daughterRefVector(); + bool is_leaf = (daughters.empty()); + + if (is_leaf) { + LogDebug("SimTauProducer").format(" TO BE SAVED {} ", resonance_idx); + auto const& gen_particle_barcode = gen_particle_barcodes[gen_particle_key]; + auto const& found_in_caloparticles = std::find_if(caloPartVec.begin(), caloPartVec.end(), [&](const auto& p) { + return p.g4Tracks()[0].genpartIndex() == gen_particle_barcode; + }); + if (found_in_caloparticles != caloPartVec.end()) { + auto calo_particle_idx = (found_in_caloparticles - caloPartVec.begin()); + t.calo_particle_leaves.push_back(CaloParticleRef(calo_particle_h, calo_particle_idx)); + t.leaves.push_back( + {gen_particle.pdgId(), resonance_idx, (int)t.calo_particle_leaves.size() - 1, gen_particle_key}); + LogDebug("SimTauProducer").format(" CP {} {}", calo_particle_idx, caloPartVec[calo_particle_idx]); + } else { + t.leaves.push_back({gen_particle.pdgId(), resonance_idx, -1, gen_particle_key}); + } + return; + } else if (generation != 0) { + t.resonances.push_back({gen_particle.pdgId(), resonance_idx}); + resonance_idx = t.resonances.size() - 1; + LogDebug("SimTauProducer").format(" RESONANCE/INTERMEDIATE {} ", resonance_idx); + } + + ++generation; + for (auto daughter = daughters.begin(); daughter != daughters.end(); ++daughter) { + int gen_particle_key = (*daughter).key(); + LogDebug("SimTauProducer").format(" gen {} {} {} ", generation, gen_particle_key, (*daughter)->pdgId()); + buildSimTau(t, generation, resonance_idx, *(*daughter), gen_particle_key, calo_particle_h, gen_particle_barcodes); + } +} + +void SimTauProducer::produce(edm::StreamID, edm::Event& iEvent, edm::EventSetup const& iSetup) const { + using namespace edm; + + auto caloParticles_h = iEvent.getHandle(caloParticles_token_); + const auto& genParticles = iEvent.get(genParticles_token_); + const auto& genBarcodes = iEvent.get(genBarcodes_token_); + + auto tauDecayVec = std::make_unique>(); + for (auto const& g : genParticles) { + auto const& flags = g.statusFlags(); + if (std::abs(g.pdgId()) == 15 and flags.isPrompt() and flags.isDecayedLeptonHadron()) { + SimTauCPLink t; + buildSimTau(t, 0, -1, g, -1, caloParticles_h, genBarcodes); + t.decayMode = t.buildDecayModes(); +#ifdef EDM_ML_DEBUG + t.dumpFullDecay(); + t.dump(); +#endif + (*tauDecayVec).push_back(t); + } + } + iEvent.put(std::move(tauDecayVec)); +} + +void SimTauProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("caloParticles", edm::InputTag("mix", "MergedCaloTruth")); + desc.add("genParticles", edm::InputTag("genParticles")); + desc.add("genBarcodes", edm::InputTag("genParticles")); + descriptions.add("SimTauProducer", desc); +} + +DEFINE_FWK_MODULE(SimTauProducer); diff --git a/SimDataFormats/CaloAnalysis/interface/SimTauCPLink.h b/SimDataFormats/CaloAnalysis/interface/SimTauCPLink.h new file mode 100644 index 0000000000000..81237ad745a7d --- /dev/null +++ b/SimDataFormats/CaloAnalysis/interface/SimTauCPLink.h @@ -0,0 +1,157 @@ +#ifndef SimDataFormats_SimTauCPLink_h +#define SimDataFormats_SimTauCPLink_h + +#include "DataFormats/HGCalReco/interface/TICLCandidate.h" +#include "SimDataFormats/CaloAnalysis/interface/CaloParticleFwd.h" +#include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h" +#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +class SimTauCPLink { +public: + SimTauCPLink(){}; + ~SimTauCPLink(){}; + struct DecayNav { + int pdgId_; + int resonance_idx_; + int calo_particle_idx_; + int gen_particle_idx_; + int pdgId() const { return pdgId_; } + int resonance_idx() const { return resonance_idx_; } + int calo_particle_idx() const { return calo_particle_idx_; } + int gen_particle_idx() const { return gen_particle_idx_; } + + private: + }; + + std::vector> resonances; + std::vector leaves; + CaloParticleRefVector calo_particle_leaves; + int decayMode; + + enum decayModes { + kNull = -1, + kOneProng0PiZero, + kOneProng1PiZero, + kOneProng2PiZero, + kOneProng3PiZero, + kOneProngNPiZero, + kTwoProng0PiZero, + kTwoProng1PiZero, + kTwoProng2PiZero, + kTwoProng3PiZero, + kTwoProngNPiZero, + kThreeProng0PiZero, + kThreeProng1PiZero, + kThreeProng2PiZero, + kThreeProng3PiZero, + kThreeProngNPiZero, + kRareDecayMode, + kElectron, + kMuon + }; + + void dump(void) const { + for (auto const &l : leaves) { + LogDebug("SimTauProducer") + .format( + "L {} {} CP: {} GenP idx: {}", l.pdgId(), l.resonance_idx(), l.calo_particle_idx(), l.gen_particle_idx()); + } + for (auto const &r : resonances) { + LogDebug("SimTauProducer").format("R {} {}", r.first, r.second); + } + } + + void dumpDecay(const std::pair &entry) const { + if (entry.second == -1) { // No intermediate mother. + LogDebug("SimTauProducer").format("{} {}", entry.first, entry.second); + } else { + LogDebug("SimTauProducer").format("{} {} coming from: ", entry.first, entry.second); + auto const &mother = resonances[entry.second]; + dumpDecay(mother); + } + } + + void dumpFullDecay(void) const { + for (auto const &leaf : leaves) { + dumpDecay({leaf.pdgId(), leaf.resonance_idx()}); + } + } + + int buildDecayModes() { + int numElectrons = 0; + int numMuons = 0; + int numHadrons = 0; + int numPhotons = 0; + for (auto leaf : leaves) { + int pdg_id = abs(leaf.pdgId()); + switch (pdg_id) { + case 22: + numPhotons++; + break; + case 11: + numElectrons++; + break; + case 13: + numMuons++; + break; + case 16: + break; + default: + numHadrons++; + } + } + + if (numElectrons == 1) + return kElectron; + else if (numMuons == 1) + return kMuon; + switch (numHadrons) { + case 1: + switch (numPhotons) { + case 0: + return kOneProng0PiZero; + case 2: + return kOneProng1PiZero; + case 4: + return kOneProng2PiZero; + case 6: + return kOneProng3PiZero; + default: + return kOneProngNPiZero; + } + case 2: + switch (numPhotons) { + case 0: + return kTwoProng0PiZero; + case 2: + return kTwoProng1PiZero; + case 4: + return kTwoProng2PiZero; + case 6: + return kTwoProng3PiZero; + default: + return kTwoProngNPiZero; + } + case 3: + switch (numPhotons) { + case 0: + return kThreeProng0PiZero; + case 2: + return kThreeProng1PiZero; + case 4: + return kThreeProng2PiZero; + case 6: + return kThreeProng3PiZero; + default: + return kThreeProngNPiZero; + } + default: + return kRareDecayMode; + } + } + +private: +}; + +#endif //SimTauCPLink diff --git a/SimDataFormats/CaloAnalysis/src/classes.h b/SimDataFormats/CaloAnalysis/src/classes.h index 4ae821b748489..e176ec45d4b66 100644 --- a/SimDataFormats/CaloAnalysis/src/classes.h +++ b/SimDataFormats/CaloAnalysis/src/classes.h @@ -10,3 +10,5 @@ #include "SimDataFormats/CaloAnalysis/interface/MtdSimLayerClusterFwd.h" #include "SimDataFormats/CaloAnalysis/interface/MtdSimTrackster.h" #include "SimDataFormats/CaloAnalysis/interface/MtdSimTracksterFwd.h" +#include "SimDataFormats/CaloAnalysis/interface/SimTauCPLink.h" +#include "DataFormats/Common/interface/Wrapper.h" diff --git a/SimDataFormats/CaloAnalysis/src/classes_def.xml b/SimDataFormats/CaloAnalysis/src/classes_def.xml index 35ee606d29a93..35eca588e8cb9 100644 --- a/SimDataFormats/CaloAnalysis/src/classes_def.xml +++ b/SimDataFormats/CaloAnalysis/src/classes_def.xml @@ -52,6 +52,12 @@ + + + + + + diff --git a/SimGeneral/Debugging/plugins/BuildFile.xml b/SimGeneral/Debugging/plugins/BuildFile.xml index 6660b275280f6..bb457e1eaf2b6 100644 --- a/SimGeneral/Debugging/plugins/BuildFile.xml +++ b/SimGeneral/Debugging/plugins/BuildFile.xml @@ -9,4 +9,6 @@ + + diff --git a/SimGeneral/Debugging/plugins/SimTauAnalyzer.cc b/SimGeneral/Debugging/plugins/SimTauAnalyzer.cc new file mode 100644 index 0000000000000..96fbfea2c3c2a --- /dev/null +++ b/SimGeneral/Debugging/plugins/SimTauAnalyzer.cc @@ -0,0 +1,65 @@ +// -*- C++ -*- +// +// +// Original Author: Andreas Gruber +// Created: Mon, 16 Oct 2023 14:24:35 GMT +// +// + +// system include files +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/one/EDAnalyzer.h" +#include "FWCore/Framework/interface/ConsumesCollector.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "TH1.h" + +#include "SimDataFormats/CaloAnalysis/interface/CaloParticleFwd.h" +#include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h" +#include "SimDataFormats/CaloAnalysis/interface/SimTauCPLink.h" + +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "CommonTools/UtilAlgos/interface/TFileService.h" + +class SimTauAnalyzer : public edm::one::EDAnalyzer { +public: + explicit SimTauAnalyzer(const edm::ParameterSet&); + ~SimTauAnalyzer() override = default; + +private: + void analyze(const edm::Event&, const edm::EventSetup&) override; + + const edm::EDGetTokenT> simTau_token_; + TH1D* DM_histo; +}; + +SimTauAnalyzer::SimTauAnalyzer(const edm::ParameterSet& iConfig) + : simTau_token_(consumes>(iConfig.getParameter("simTau"))) { + usesResource(TFileService::kSharedResource); + edm::Service fs; + DM_histo = fs->make("DM_histo", "DM_histo", 20, -1, 19); +} + +// ------------ method called for each event ------------ +void SimTauAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) { + edm::Handle> simTau_h; + iEvent.getByToken(simTau_token_, simTau_h); + + const auto& simTaus = *simTau_h; + + for (auto const& simTau : simTaus) { +#ifdef EDM_ML_DEBUG + simTau.dumpFullDecay(); + simTau.dump(); +#endif + DM_histo->Fill(simTau.decayMode); + } +} + +//define this as a plug-in +DEFINE_FWK_MODULE(SimTauAnalyzer); diff --git a/SimGeneral/Debugging/python/caloParticleDebugger_cfg.py b/SimGeneral/Debugging/python/caloParticleDebugger_cfg.py index f60c86ca5bea0..5b90c445f23a3 100644 --- a/SimGeneral/Debugging/python/caloParticleDebugger_cfg.py +++ b/SimGeneral/Debugging/python/caloParticleDebugger_cfg.py @@ -8,7 +8,7 @@ process = cms.Process("Demo") process.load("FWCore.MessageService.MessageLogger_cfi") -process.load('Configuration.Geometry.GeometryExtended2026D96Reco_cff') +process.load('Configuration.Geometry.GeometryExtended2026D96_cff') process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') from Configuration.AlCa.GlobalTag import GlobalTag process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:phase2_realistic_T21', '') diff --git a/SimGeneral/Debugging/test/simTauAnalyzer_cfg.py b/SimGeneral/Debugging/test/simTauAnalyzer_cfg.py new file mode 100644 index 0000000000000..5625bf1bccb57 --- /dev/null +++ b/SimGeneral/Debugging/test/simTauAnalyzer_cfg.py @@ -0,0 +1,27 @@ +import FWCore.ParameterSet.Config as cms +import sys + +process = cms.Process("SimTauAnalyzer") + +process.load("FWCore.MessageService.MessageLogger_cfi") +process.MessageLogger.cerr.threshold = "DEBUG" +process.MessageLogger.debugModules = ["SimTauCPLink", "SimTauProducer", "SimTauAnalyzer"] + +process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(-1) + ) +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring( + 'file:/data/agruber/patatrack/CMSSW_14_0_0_pre0/src/24036.0_ZTT_14TeV+2026D96/step3.root' + #'file:SimTauProducer_test.root' + ) + ) + +process.TFileService = cms.Service("TFileService", + fileName = cms.string('test.root') + ) + +process.SimTauAnalyzer = cms.EDAnalyzer('SimTauAnalyzer', +simTau = cms.InputTag('SimTauCPLink') + ) + +process.p = cms.Path(process.SimTauAnalyzer) diff --git a/Validation/Configuration/python/hgcalSimValid_cff.py b/Validation/Configuration/python/hgcalSimValid_cff.py index 142fda660ed81..89ff13ad37371 100644 --- a/Validation/Configuration/python/hgcalSimValid_cff.py +++ b/Validation/Configuration/python/hgcalSimValid_cff.py @@ -10,6 +10,7 @@ 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 SimCalorimetry.HGCalAssociatorProducers.SimTauProducer_cfi import * from Validation.HGCalValidation.simhitValidation_cff import * from Validation.HGCalValidation.digiValidation_cff import * @@ -38,7 +39,8 @@ simTsAssocByEnergyScoreProducer, simTracksterHitLCAssociatorByEnergyScoreProducer, tracksterSimTracksterAssociationLinking, tracksterSimTracksterAssociationPR, tracksterSimTracksterAssociationLinkingbyCLUE3D, tracksterSimTracksterAssociationPRbyCLUE3D, - tracksterSimTracksterAssociationLinkingPU, tracksterSimTracksterAssociationPRPU + tracksterSimTracksterAssociationLinkingPU, tracksterSimTracksterAssociationPRPU, + SimTauProducer ) hgcalValidation = cms.Sequence(hgcalSimHitValidationEE