Skip to content

Fix orphan sim clusters #47883

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
cmsbuild merged 3 commits into from
Apr 21, 2025
Merged

Fix orphan sim clusters #47883

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
510bbfb
Fix bug causing orphan SimClusters.
rovere Apr 4, 2025
2271dfb
Improve CaloParticleDebugger printouts.
rovere Apr 4, 2025
a912731
Improve SimTauCPLinkProducer debug printouts.
rovere Apr 16, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 4 additions & 3 deletions SimCalorimetry/HGCalAssociatorProducers/plugins/SimTauCPLinkProducer.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,8 +69,9 @@ void SimTauProducer::buildSimTau(SimTauCPLink& t,
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];
LogDebug("SimTauProducer")
.format(" TO BE SAVED {}, key {}, barcode {}", resonance_idx, gen_particle_key, gen_particle_barcode);
auto const& found_in_caloparticles = std::find_if(caloPartVec.begin(), caloPartVec.end(), [&](const auto& p) {
return p.g4Tracks()[0].genpartIndex() == gen_particle_barcode;
});
Expand All @@ -87,13 +88,13 @@ void SimTauProducer::buildSimTau(SimTauCPLink& t,
} 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);
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());
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);
}
}
Expand Down
49 changes: 28 additions & 21 deletions SimDataFormats/CaloAnalysis/interface/SimTauCPLink.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,31 +31,38 @@ class SimTauCPLink {

enum decayModes {
kNull = -1,
kOneProng0PiZero,
kOneProng1PiZero,
kOneProng2PiZero,
kOneProng3PiZero,
kOneProngNPiZero,
kTwoProng0PiZero,
kTwoProng1PiZero,
kTwoProng2PiZero,
kTwoProng3PiZero,
kTwoProngNPiZero,
kThreeProng0PiZero,
kThreeProng1PiZero,
kThreeProng2PiZero,
kThreeProng3PiZero,
kThreeProngNPiZero,
kRareDecayMode,
kElectron,
kMuon
kOneProng0PiZero, // 0
kOneProng1PiZero, // 1
kOneProng2PiZero, // 2
kOneProng3PiZero, // 3
kOneProngNPiZero, // 4
kTwoProng0PiZero, // 5
kTwoProng1PiZero, // 6
kTwoProng2PiZero, // 7
kTwoProng3PiZero, // 8
kTwoProngNPiZero, // 9
kThreeProng0PiZero, // 10
kThreeProng1PiZero, // 11
kThreeProng2PiZero, // 12
kThreeProng3PiZero, // 13
kThreeProngNPiZero, // 14
kRareDecayMode, // 15
kElectron, // 16
kMuon // 17
};

void dump(void) const {
LogDebug("SimTauProducer")
.format("Decay mode: {} ", buildDecayModes())
.format("Leaves: {} ", leaves.size())
.format("Resonances: {}", resonances.size());
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());
.format("L {} {} CP: {} GenP idx: {}",
l.pdgId(),
l.resonance_idx(),
(int)((l.calo_particle_idx() == -1) ? -1 : calo_particle_leaves[l.calo_particle_idx()].key()),
l.gen_particle_idx());
}
for (auto const &r : resonances) {
LogDebug("SimTauProducer").format("R {} {}", r.first, r.second);
Expand All @@ -78,7 +85,7 @@ class SimTauCPLink {
}
}

int buildDecayModes() {
int buildDecayModes() const {
int numElectrons = 0;
int numMuons = 0;
int numHadrons = 0;
Expand Down
10 changes: 8 additions & 2 deletions SimGeneral/CaloAnalysis/plugins/CaloTruthAccumulator.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -160,7 +160,8 @@ namespace {
simHitBarcodeToIndex_(simHitBarcodeToIndex),
simTrackDetIdEnergyMap_(simTrackDetIdEnergyMap),
vertex_time_map_(vertex_time_map),
selector_(selector) {}
selector_(selector),
insideCP_(false) {}
template <typename Vertex, typename Graph>
void discover_vertex(Vertex u, const Graph &g) {
// If we reach the vertex 0, it means that we are backtracking with respect
Expand All @@ -173,7 +174,7 @@ namespace {
auto trackIdx = vertex_property.simTrack->trackId();
IfLogDebug(DEBUG, messageCategoryGraph_)
<< " Found " << simHitBarcodeToIndex_.count(trackIdx) << " associated simHits" << std::endl;
if (simHitBarcodeToIndex_.count(trackIdx)) {
if (insideCP_ && simHitBarcodeToIndex_.count(trackIdx)) {
output_.pSimClusters->emplace_back(*vertex_property.simTrack);
auto &simcluster = output_.pSimClusters->back();
std::unordered_map<uint32_t, float> acc_energy;
Expand All @@ -193,10 +194,13 @@ namespace {
auto edge_property = get(edge_weight, g, e);
IfLogDebug(DEBUG, messageCategoryGraph_) << "Considering CaloParticle: " << edge_property.simTrack->trackId();
if (selector_(edge_property)) {
insideCP_ = true;
IfLogDebug(DEBUG, messageCategoryGraph_) << "Adding CaloParticle: " << edge_property.simTrack->trackId();
output_.pCaloParticles->emplace_back(*(edge_property.simTrack));
output_.pCaloParticles->back().setSimTime(vertex_time_map_[(edge_property.simTrack)->vertIndex()]);
caloParticles_.sc_start_.push_back(output_.pSimClusters->size());
} else {
insideCP_ = false;
}
}
}
Expand All @@ -209,6 +213,7 @@ namespace {
auto edge_property = get(edge_weight, g, e);
if (selector_(edge_property)) {
caloParticles_.sc_stop_.push_back(output_.pSimClusters->size());
insideCP_ = false;
}
}
}
Expand All @@ -220,6 +225,7 @@ namespace {
std::unordered_map<int, std::map<int, float>> &simTrackDetIdEnergyMap_;
std::unordered_map<uint32_t, float> &vertex_time_map_;
Selector selector_;
bool insideCP_;
};
} // namespace

Expand Down
32 changes: 22 additions & 10 deletions SimGeneral/Debugging/plugins/CaloParticleDebugger.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -163,17 +163,19 @@ void CaloParticleDebugger::analyze(const edm::Event& iEvent, const edm::EventSet
auto const& t = tracks[i];
LogVerbatim("CaloParticleDebuggerSimTracks")
<< i << "\t" << t.trackId() << "\t" << t << " Crossed Boundary: " << t.crossedBoundary()
<< " Vtx: " << t.vertIndex() << " isFromBackScattering: " << t.isFromBackScattering()
<< " isPrimary: " << t.isPrimary() << " ParentID: " << t.getPrimaryID()
<< " Position Boundary: " << t.getPositionAtBoundary() << " Momentum Boundary: " << t.getMomentumAtBoundary()
<< " Vtx: " << t.vertIndex() << " Momemtum Origin: " << t.momentum();
<< " Momemtum Origin: " << t.momentum();
trackid_to_track_index[t.trackId()] = i;
}

LogVerbatim("CaloParticleDebuggerGenParticles") << "\n\n**Printing GenParticles information **";
LogVerbatim("CaloParticleDebuggerGenParticles") << "BARCODE\tPDGID\tMOMENTUM(x,y,z)\tVertex(x,y,z)";
LogVerbatim("CaloParticleDebuggerGenParticles") << "BARCODE\tID\tPDGID\tMOMENTUM(x,y,z)\tVertex(x,y,z)";
for (size_t i = 0; i < genParticles.size(); ++i) {
auto const& gp = genParticles[i];
LogVerbatim("CaloParticleDebuggerGenParticles")
<< genBarcodes[i] << "\t" << gp.pdgId() << "\t" << gp.momentum() << "\t" << gp.vertex();
<< genBarcodes[i] << "\t" << i << "\t" << gp.pdgId() << "\t" << gp.momentum() << "\t" << gp.vertex();
}

LogVerbatim("CaloParticleDebuggerSimVertices") << "\n\n**Printing SimVertex information **";
Expand All @@ -190,32 +192,38 @@ void CaloParticleDebugger::analyze(const edm::Event& iEvent, const edm::EventSet
}

LogVerbatim("CaloParticleDebuggerCaloParticles") << "\n\n**Printing CaloParticles information **";
int totalSimCl_in_event = 0;
std::set<int> simClusters_in_CaloParticles;
for (size_t i = 0; i < calopart.size(); ++i) {
auto const& cp = calopart[i];
LogVerbatim("CaloParticleDebuggerCaloParticles") << "\n\n"
<< i << "\tType: " << cp.pdgId() << "\tEnergy: " << cp.energy()
<< "\tBarcode: " << cp.g4Tracks()[0].genpartIndex()
<< " G4_trackID: " << cp.g4Tracks()[0].trackId(); // << cp ;
double total_sim_energy = 0.;
double total_cp_energy = 0.;
LogVerbatim("CaloParticleDebuggerCaloParticles") << "--> Overall simclusters in CP: " << cp.simClusters().size();
totalSimCl_in_event += cp.simClusters().size();
// All the next mess just to print the simClusters ordered
auto const& simcs = cp.simClusters();
for (size_t j = 0; j < simcs.size(); ++j) {
LogVerbatim("CaloParticleDebuggerCaloParticles") << *(simcs[j]);
simClusters_in_CaloParticles.insert(simcs[j]->g4Tracks()[0].trackId());
}

for (auto const& sc : cp.simClusters()) {
for (auto const& cl : sc->hits_and_fractions()) {
total_sim_energy += detIdToTotalSimEnergy[cl.first] * cl.second;
total_cp_energy += cp.energy() * cl.second;
}
}
LogVerbatim("CaloParticleDebuggerCaloParticles")
<< "--> Overall SC energy (sum using sim energies): " << total_sim_energy;
LogVerbatim("CaloParticleDebuggerCaloParticles")
<< "--> Overall SC energy (sum using CaloP energies): " << total_cp_energy;
}
LogVerbatim("CaloParticleDebuggerCaloParticles")
<< "--> Overall SCs attached to CPs: " << totalSimCl_in_event;
LogVerbatim("CaloParticleDebuggerCaloParticles")
<< "--> Overall SCs in the event: " << simclusters.size();
LogVerbatim("CaloParticleDebuggerCaloParticles")
<< "--> Orphan SCs in the event: " << simclusters.size() - totalSimCl_in_event;

LogVerbatim("CaloParticleDebuggerSimClusters") << "\n\n**Printing SimClusters information **";
for (size_t i = 0; i < simclusters.size(); ++i) {
Expand All @@ -224,9 +232,13 @@ void CaloParticleDebugger::analyze(const edm::Event& iEvent, const edm::EventSet
<< "\n\n"
<< i << "\tType: " << simcl.pdgId() << "\tEnergy: " << simcl.energy() << "\tKey: " << i; // << simcl ;
auto const& simtrack = simcl.g4Tracks()[0];
LogVerbatim("CaloParticleDebuggerSimClusters") << " Crossed Boundary: " << simtrack.crossedBoundary()
<< " Position Boundary: " << simtrack.getPositionAtBoundary()
<< " Momentum Boundary: " << simtrack.getMomentumAtBoundary();
LogVerbatim("CaloParticleDebuggerSimClusters") << "\n Primary GenParticleID: " << simtrack.getPrimaryID()
<< "\n Crossed Boundary: " << simtrack.crossedBoundary()
<< "\n Position Boundary: " << simtrack.getPositionAtBoundary()
<< "\n Momentum Boundary: " << simtrack.getMomentumAtBoundary();
if (simClusters_in_CaloParticles.find(simcl.g4Tracks()[0].trackId()) == simClusters_in_CaloParticles.end()) {
LogVerbatim("CaloParticleDebuggerSimClusters") << " Orphan SimCluster: " << simtrack.trackId();
}
double total_sim_energy = 0.;
LogVerbatim("CaloParticleDebuggerSimClusters") << "--> Overall simclusters's size: " << simcl.numberOfRecHits();
for (auto const& cl : simcl.hits_and_fractions()) {
Expand Down
81 changes: 48 additions & 33 deletions SimGeneral/Debugging/python/caloParticleDebugger_cfg.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,60 +2,75 @@

# The line below always has to be included to make VarParsing work
from FWCore.ParameterSet.VarParsing import VarParsing
options = VarParsing ('analysis')

options = VarParsing("analysis")
options.parseArguments()

process = cms.Process("Demo")

process.load("FWCore.MessageService.MessageLogger_cfi")
process.load('Configuration.Geometry.GeometryExtended2026D96_cff')
process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff')
process.load("Configuration.Geometry.GeometryExtendedRun4D110Reco_cff")
process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff")
from Configuration.AlCa.GlobalTag import GlobalTag
process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:phase2_realistic_T21', '')

process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(-1) )
process.GlobalTag = GlobalTag(process.GlobalTag, "auto:phase2_realistic_T21", "")

process.maxEvents = cms.untracked.PSet(input=cms.untracked.int32(-1))

input_filename = 'default.root' if len(options.inputFiles) == 0 else options.inputFiles[0]
#input_filename='step2SingleElectronPt15Eta1p7_2p7_SimTracksters.root'
#input_filename='step2FineCaloSingleElectronPt15Eta1p7_2p7_SimTracksters.root'
#input_filename='step2SingleElectronPt15Eta1p7_2p7_CBWEAndSimTracksters.root'
#input_filename='step2FineCaloSingleElectronPt15Eta1p7_2p7_CBWEAndSimTracksters.root'

process.source = cms.Source("PoolSource",
inputCommands = cms.untracked.vstring(['keep *',
'drop l1tEMTFHit2016Extras_simEmtfDigis_CSC_HLT',
'drop l1tEMTFHit2016Extras_simEmtfDigis_RPC_HLT',
'drop l1tEMTFHit2016s_simEmtfDigis__HLT',
'drop l1tEMTFTrack2016Extras_simEmtfDigis__HLT',
'drop l1tEMTFTrack2016s_simEmtfDigis__HLT']),
# replace 'myfile.root' with the source file you want to use
fileNames = cms.untracked.vstring(
# 'file:/data/rovere/HGCAL/study/CMSSW_9_4_0/src/SimGeneral/Debugging/test/20800.0_FourMuPt1_200+FourMuPt_1_200_pythia8_2023D20_GenSimHLBeamSpotFull+DigiFull_2023D20+RecoFullGlobal_2023D20+HARVESTFullGlobal_2023D20/step2.root'
# 'file:/data/rovere/HGCAL/study/CMSSW_9_4_0/src/SimGeneral/Debugging/test/20824.0_TTbar_13+TTbar_13TeV_TuneCUETP8M1_2023D20_GenSimHLBeamSpotFull+DigiFull_2023D20+RecoFullGlobal_2023D20+HARVESTFullGlobal_2023D20/step2.root'
# 'file:/data/rovere/HGCAL/study/CMSSW_9_4_0/src/SimGeneral/Debugging/test/20002.0_SingleElectronPt35+SingleElectronPt35_pythia8_2023D17_GenSimHLBeamSpotFull+DigiFullTrigger_2023D17+RecoFullGlobal_2023D17+HARVESTFullGlobal_2023D17/step2.root'
# 'file:/data/rovere/HGCAL/study/CMSSW_9_4_0/src/SimGeneral/Debugging/test/20016.0_SingleGammaPt35Extended+DoubleGammaPt35Extended_pythia8_2023D17_GenSimHLBeamSpotFull+DigiFullTrigger_2023D17+RecoFullGlobal_2023D17+HARVESTFullGlobal_2023D17/step2.root'
# 'file:/data/rovere/HGCAL/study/CMSSW_9_4_0/src/SimGeneral/Debugging/test/20088.0_SinglePiPt25Eta1p7_2p7+SinglePiPt25Eta1p7_2p7_2023D17_GenSimHLBeamSpotFull+DigiFullTrigger_2023D17+RecoFullGlobal_2023D17+HARVESTFullGlobal_2023D17/step2.root'
'file:%s'%input_filename
input_filename = (
"default.root" if len(options.inputFiles) == 0 else options.inputFiles[0]
)
# input_filename='step2SingleElectronPt15Eta1p7_2p7_SimTracksters.root'
# input_filename='step2FineCaloSingleElectronPt15Eta1p7_2p7_SimTracksters.root'
# input_filename='step2SingleElectronPt15Eta1p7_2p7_CBWEAndSimTracksters.root'
# input_filename='step2FineCaloSingleElectronPt15Eta1p7_2p7_CBWEAndSimTracksters.root'

)
process.source = cms.Source(
"PoolSource",
inputCommands=cms.untracked.vstring(
[
"keep *",
"drop l1tEMTFHit2016Extras_simEmtfDigis_CSC_HLT",
"drop l1tEMTFHit2016Extras_simEmtfDigis_RPC_HLT",
"drop l1tEMTFHit2016s_simEmtfDigis__HLT",
"drop l1tEMTFTrack2016Extras_simEmtfDigis__HLT",
"drop l1tEMTFTrack2016s_simEmtfDigis__HLT",
]
),
# replace 'myfile.root' with the source file you want to use
fileNames=cms.untracked.vstring(
# 'file:/data/rovere/HGCAL/study/CMSSW_9_4_0/src/SimGeneral/Debugging/test/20800.0_FourMuPt1_200+FourMuPt_1_200_pythia8_2023D20_GenSimHLBeamSpotFull+DigiFull_2023D20+RecoFullGlobal_2023D20+HARVESTFullGlobal_2023D20/step2.root'
# 'file:/data/rovere/HGCAL/study/CMSSW_9_4_0/src/SimGeneral/Debugging/test/20824.0_TTbar_13+TTbar_13TeV_TuneCUETP8M1_2023D20_GenSimHLBeamSpotFull+DigiFull_2023D20+RecoFullGlobal_2023D20+HARVESTFullGlobal_2023D20/step2.root'
# 'file:/data/rovere/HGCAL/study/CMSSW_9_4_0/src/SimGeneral/Debugging/test/20002.0_SingleElectronPt35+SingleElectronPt35_pythia8_2023D17_GenSimHLBeamSpotFull+DigiFullTrigger_2023D17+RecoFullGlobal_2023D17+HARVESTFullGlobal_2023D17/step2.root'
# 'file:/data/rovere/HGCAL/study/CMSSW_9_4_0/src/SimGeneral/Debugging/test/20016.0_SingleGammaPt35Extended+DoubleGammaPt35Extended_pythia8_2023D17_GenSimHLBeamSpotFull+DigiFullTrigger_2023D17+RecoFullGlobal_2023D17+HARVESTFullGlobal_2023D17/step2.root'
# 'file:/data/rovere/HGCAL/study/CMSSW_9_4_0/src/SimGeneral/Debugging/test/20088.0_SinglePiPt25Eta1p7_2p7+SinglePiPt25Eta1p7_2p7_2023D17_GenSimHLBeamSpotFull+DigiFullTrigger_2023D17+RecoFullGlobal_2023D17+HARVESTFullGlobal_2023D17/step2.root'
"file:%s" % input_filename
),
)

process.load("SimGeneral.Debugging.caloParticleDebugger_cfi")

# MessageLogger customizations
process.MessageLogger.cerr.enable = False
process.MessageLogger.cout.enable = False
labels = ['SimTracks', 'SimVertices', 'GenParticles', 'TrackingParticles', 'CaloParticles', 'SimClusters']
labels = [
"SimTracks",
"SimVertices",
"GenParticles",
"TrackingParticles",
"CaloParticles",
"SimClusters",
]
messageLogger = dict()
for category in labels:
main_key = '%sMessageLogger'%(category)
category_key = 'CaloParticleDebugger%s'%(category)
main_key = "%sMessageLogger" % (category)
category_key = "CaloParticleDebugger%s" % (category)
messageLogger[main_key] = dict(
filename = '%s_%s.log' % (input_filename.replace('.root',''), category),
threshold = 'INFO',
default = dict(limit=0)
)
filename="%s_%s.log" % (input_filename.replace(".root", ""), category),
threshold="INFO",
default=dict(limit=0),
)
messageLogger[main_key][category_key] = dict(limit=-1)
# First create defaults
setattr(process.MessageLogger.files, category, dict())
Expand Down
5 changes: 4 additions & 1 deletion SimGeneral/MixingModule/interface/DecayGraph.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -103,7 +103,10 @@ namespace {
std::string graphviz_vertex(const VertexProperty &v) {
std::ostringstream oss;
oss << "{id: " << (v.simTrack ? v.simTrack->trackId() : 0) << ",\\ntype: " << (v.simTrack ? v.simTrack->type() : 0)
<< ",\\nchits: " << v.cumulative_simHits << "}";
<< ",\\nPVtxIdx?: " << (v.simTrack ? v.simTrack->vertIndex() : -1)
<< ",\\nPrim?: " << (v.simTrack ? v.simTrack->isPrimary() : -1)
<< ",\\nXB?: " << (v.simTrack ? v.simTrack->crossedBoundary() : -1) << ",\\nchits: " << v.cumulative_simHits
<< "}";
return oss.str();
}

Expand Down