From 70e4fa24e087b517243cc03e98a715fcc9c49cb1 Mon Sep 17 00:00:00 2001 From: Jacopo Fanini Date: Thu, 27 Aug 2026 08:33:02 +0200 Subject: [PATCH 1/3] [Gaudi] New EDM4hep reader and plotter --- .gitignore | 1 + .../options/readPlotEDM4hep.py | 180 ++++++++++++++++++ 2 files changed, 181 insertions(+) create mode 100644 GaudiTutorial/RandomNoiseDigitizer/options/readPlotEDM4hep.py diff --git a/.gitignore b/.gitignore index 6725208..769c1eb 100644 --- a/.gitignore +++ b/.gitignore @@ -251,3 +251,4 @@ test/gaudi_opts/testConverterConstants.py # Files produced during running examples *root *png +*pdf diff --git a/GaudiTutorial/RandomNoiseDigitizer/options/readPlotEDM4hep.py b/GaudiTutorial/RandomNoiseDigitizer/options/readPlotEDM4hep.py new file mode 100644 index 0000000..be23de8 --- /dev/null +++ b/GaudiTutorial/RandomNoiseDigitizer/options/readPlotEDM4hep.py @@ -0,0 +1,180 @@ +import os + +# Disable ROOT web display on remote machines +os.environ["ROOT_WEBDISPLAY"] = "off" + +import podio +import ROOT + + +input_file = "../../data/simpleCalo_noiseDigitizer.root" + +# detector geometry parameters +n_layers = 20 +calo_x = 1000.0 # mm +calo_y = 1000.0 # mm +calo_z = 2000.0 # mm +layer_thickness = 100.0 # mm +n_cells_x = 10 +n_cells_y = 10 + +# Histos for transverse energy profile histograms for each layer +h_transverse = [] +for layer in range(n_layers): + h = ROOT.TH2F( + f"h_transverse_{layer}", + f"Layer {layer}", + n_cells_x, -calo_x / 2, calo_x / 2, + n_cells_y, -calo_y / 2, calo_y / 2, + ) + h.SetStats(0) + h_transverse.append(h) + +# Histo for longitudinal energy profile +h_longitudinal = ROOT.TH1F( + "h_longitudinal", + "Longitudinal energy profile;Layer;Energy [MeV]", + n_layers, -0.5, n_layers - 0.5 +) + + +# Read EDM4hep file using podio and python bindings +reader = podio.root_io.Reader(input_file) + +n_events = 0 +cell_diffs = [] + +# Event loop +for event in reader.get("events"): + # Get relevant collections + simhits = event.get("simplecaloRO") + digihits = event.get("CaloDigiHits") + # Dictionaries for CellID-based sim-digi comparison + sim_by_cell = {} + digi_by_cell = {} + + for hit in digihits: + # Store hit by CellID for sim-digi comparison + digi_by_cell[hit.getCellID()] = hit + + # Get hit energy and position + energy = hit.getEnergy() * 1000.0 # GeV -> MeV + pos = hit.getPosition() + + # Determine calorimeter layer from z + layer = int((pos.z + calo_z / 2.0) / layer_thickness) + if 0 <= layer < n_layers: + h_transverse[layer].Fill(pos.x, pos.y, energy) + h_longitudinal.Fill(layer, energy) + + for hit in simhits: + sim_by_cell[hit.getCellID()] = hit + + common_cells = sim_by_cell.keys() & digi_by_cell.keys() + + for cellid in common_cells: + sim_energy = sim_by_cell[cellid].getEnergy() * 1000.0 + digi_energy = digi_by_cell[cellid].getEnergy() * 1000.0 + cell_diffs.append(digi_energy - sim_energy) + + n_events += 1 + + +# Average profiles over events +for h in h_transverse: + h.Scale(1.0 / n_events) +h_longitudinal.Scale(1.0 / n_events) + +# --------------------------------- +# Longitudinal profile +canvas_longitudinal = ROOT.TCanvas( + "canvas_longitudinal", # ROOT object name + "Longitudinal energy profile", # Canvas title + 800, # Width in pixels + 600 # Height in pxels +) +canvas_longitudinal.SetLogy() # Set log scale +h_longitudinal.SetLineColor(ROOT.kRed) +h_longitudinal.Draw("HIST") +canvas_longitudinal.SaveAs("longitudinal_profile.pdf") + +# --------------------------------- +# Transverse per-layer profiles + +# Use same energy scale for each histogram +max_energy = max(h.GetMaximum() for h in h_transverse) +for h in h_transverse: + h.SetMinimum(0.0) + h.SetMaximum(max_energy) + +# Draw all layers in one canvas +canvas_transverse = ROOT.TCanvas( + "canvas_transverse", + "Transverse energy profile by layer", + 1500, + 1100 +) +canvas_transverse.Divide(5, 4) # Split canvas in 5 columns x 4 rows +for layer, h in enumerate(h_transverse): + + canvas_transverse.cd(layer + 1) # Select corresponding canvas pad + + ROOT.gPad.SetRightMargin(0.05) + ROOT.gPad.SetLeftMargin(0.12) + ROOT.gPad.SetBottomMargin(0.12) + + h.Draw("COL") # 2D histogram with color map +canvas_transverse.SaveAs("transverse_profiles_layers.pdf") + +# --------------------------------- +# Sim-digi cell energy difference + +diff_min = min(cell_diffs) +diff_max = max(cell_diffs) +# Histogram adapts to min, max +h_cell_diff = ROOT.TH1F( + "h_cell_diff", + "Cell-by-cell digitization effect;" + "E_{digi} - E_{sim} [MeV];Cells", + 100, + diff_min, + diff_max +) +h_cell_diff.SetStats(0) + +for diff in cell_diffs: + h_cell_diff.Fill(diff) + +canvas_energy_diff = ROOT.TCanvas( + "canvas_energy_diff", + "Energy difference", + 800, + 600 +) +h_cell_diff.SetLineColor(ROOT.kGreen) + +h_cell_diff.Fit("gaus") # Gaussian fit +fit = h_cell_diff.GetFunction("gaus") +mean = fit.GetParameter(1) +sigma = fit.GetParameter(2) + +h_cell_diff.Draw("HIST") +fit.Draw("SAME") + +legend = ROOT.TLegend(0.60, 0.70, 0.88, 0.88) +legend.AddEntry(h_cell_diff, "Cell energy difference", "l") +legend.AddEntry(fit, "Gaussian fit", "l") +legend.AddEntry(0, f"#mu = {mean:.4g} MeV", "") +legend.AddEntry(0, f"#sigma = {sigma:.4g} MeV", "") +legend.Draw() + +canvas_energy_diff.SaveAs("energy_difference.pdf") + +# Save canvases to a ROOT file +output_file = ROOT.TFile("plots.root", "RECREATE") + +canvas_longitudinal.Write() +canvas_transverse.Write() +canvas_energy_diff.Write() + +output_file.Close() \ No newline at end of file From 71c13c48306e6773db053c150159f0f02bfc5ecc Mon Sep 17 00:00:00 2001 From: Jacopo Fanini Date: Sat, 29 Aug 2026 17:26:23 +0200 Subject: [PATCH 2/3] Implement ALC's suggestions --- .../options/readPlotEDM4hep.py | 35 ++++++++++++------- 1 file changed, 23 insertions(+), 12 deletions(-) diff --git a/GaudiTutorial/RandomNoiseDigitizer/options/readPlotEDM4hep.py b/GaudiTutorial/RandomNoiseDigitizer/options/readPlotEDM4hep.py index be23de8..f0c6b07 100644 --- a/GaudiTutorial/RandomNoiseDigitizer/options/readPlotEDM4hep.py +++ b/GaudiTutorial/RandomNoiseDigitizer/options/readPlotEDM4hep.py @@ -3,11 +3,22 @@ # Disable ROOT web display on remote machines os.environ["ROOT_WEBDISPLAY"] = "off" +import argparse import podio import ROOT +from dd4hep import dd4hep -input_file = "../../data/simpleCalo_noiseDigitizer.root" +parser = argparse.ArgumentParser() +parser.add_argument("-i", "--input-file", required=True, help="Input ROOT EDM4hep file") +parser.add_argument("-o", "--output-folder", default="../../data/digitizer-plots", help="Output folder for ROOT and pdf files") +args = parser.parse_args() + +#input = "../../data/simpleCalo_noiseDigitizer.root" +input_file = args.input_file +output_folder = args.output_folder + +os.makedirs(output_folder, exist_ok=True) # detector geometry parameters n_layers = 20 @@ -20,7 +31,7 @@ # Histos for transverse energy profile histograms for each layer h_transverse = [] -for layer in range(n_layers): +for layer in range(1, n_layers + 1): h = ROOT.TH2F( f"h_transverse_{layer}", f"Layer {layer}", @@ -34,10 +45,9 @@ h_longitudinal = ROOT.TH1F( "h_longitudinal", "Longitudinal energy profile;Layer;Energy [MeV]", - n_layers, -0.5, n_layers - 0.5 + n_layers, 0.5, n_layers + 0.5 ) - # Read EDM4hep file using podio and python bindings reader = podio.root_io.Reader(input_file) @@ -61,10 +71,11 @@ energy = hit.getEnergy() * 1000.0 # GeV -> MeV pos = hit.getPosition() - # Determine calorimeter layer from z - layer = int((pos.z + calo_z / 2.0) / layer_thickness) - if 0 <= layer < n_layers: - h_transverse[layer].Fill(pos.x, pos.y, energy) + # Determine calorimeter layer and fill histograms + decoder = dd4hep.BitFieldCoder("calolayer:5,abslayer:1,x:-10,y:-10") + layer = decoder.get(hit.getCellID(), "calolayer") + if 1 <= layer <= n_layers: + h_transverse[layer - 1].Fill(pos.x, pos.y, energy) h_longitudinal.Fill(layer, energy) for hit in simhits: @@ -96,7 +107,7 @@ canvas_longitudinal.SetLogy() # Set log scale h_longitudinal.SetLineColor(ROOT.kRed) h_longitudinal.Draw("HIST") -canvas_longitudinal.SaveAs("longitudinal_profile.pdf") +canvas_longitudinal.SaveAs(os.path.join(output_folder, "longitudinal_profile.pdf")) # --------------------------------- # Transverse per-layer profiles @@ -124,7 +135,7 @@ ROOT.gPad.SetBottomMargin(0.12) h.Draw("COL") # 2D histogram with color map -canvas_transverse.SaveAs("transverse_profiles_layers.pdf") +canvas_transverse.SaveAs(os.path.join(output_folder, "transverse_profiles_layers.pdf")) # --------------------------------- # Sim-digi cell energy difference @@ -168,10 +179,10 @@ legend.AddEntry(0, f"#sigma = {sigma:.4g} MeV", "") legend.Draw() -canvas_energy_diff.SaveAs("energy_difference.pdf") +canvas_energy_diff.SaveAs(os.path.join(output_folder, "energy_difference.pdf")) # Save canvases to a ROOT file -output_file = ROOT.TFile("plots.root", "RECREATE") +output_file = ROOT.TFile(os.path.join(output_folder, "plots.root"), "RECREATE") canvas_longitudinal.Write() canvas_transverse.Write() From 7e11d079fc50ab09e7abee2ba977947c29fbc6dd Mon Sep 17 00:00:00 2001 From: Jacopo Fanini Date: Tue, 1 Sep 2026 08:14:27 +0200 Subject: [PATCH 3/3] Formatting --- .../options/readPlotEDM4hep.py | 46 +++++++++++++------ 1 file changed, 32 insertions(+), 14 deletions(-) diff --git a/GaudiTutorial/RandomNoiseDigitizer/options/readPlotEDM4hep.py b/GaudiTutorial/RandomNoiseDigitizer/options/readPlotEDM4hep.py index f0c6b07..58e23f4 100644 --- a/GaudiTutorial/RandomNoiseDigitizer/options/readPlotEDM4hep.py +++ b/GaudiTutorial/RandomNoiseDigitizer/options/readPlotEDM4hep.py @@ -1,3 +1,21 @@ +# +# Copyright (c) 2020-2024 Key4hep-Project. +# +# This file is part of Key4hep. +# See https://key4hep.github.io/key4hep-doc/ for further info. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# import os # Disable ROOT web display on remote machines @@ -6,7 +24,7 @@ import argparse import podio import ROOT -from dd4hep import dd4hep +from dd4hep import dd4hep parser = argparse.ArgumentParser() @@ -15,7 +33,7 @@ args = parser.parse_args() #input = "../../data/simpleCalo_noiseDigitizer.root" -input_file = args.input_file +input_file = args.input_file output_folder = args.output_folder os.makedirs(output_folder, exist_ok=True) @@ -40,7 +58,7 @@ ) h.SetStats(0) h_transverse.append(h) - + # Histo for longitudinal energy profile h_longitudinal = ROOT.TH1F( "h_longitudinal", @@ -56,38 +74,38 @@ # Event loop for event in reader.get("events"): - # Get relevant collections + # Get relevant collections simhits = event.get("simplecaloRO") digihits = event.get("CaloDigiHits") - # Dictionaries for CellID-based sim-digi comparison + # Dictionaries for CellID-based sim-digi comparison sim_by_cell = {} digi_by_cell = {} - + for hit in digihits: # Store hit by CellID for sim-digi comparison digi_by_cell[hit.getCellID()] = hit - + # Get hit energy and position energy = hit.getEnergy() * 1000.0 # GeV -> MeV pos = hit.getPosition() - + # Determine calorimeter layer and fill histograms decoder = dd4hep.BitFieldCoder("calolayer:5,abslayer:1,x:-10,y:-10") layer = decoder.get(hit.getCellID(), "calolayer") if 1 <= layer <= n_layers: - h_transverse[layer - 1].Fill(pos.x, pos.y, energy) + h_transverse[layer - 1].Fill(pos.x, pos.y, energy) h_longitudinal.Fill(layer, energy) - + for hit in simhits: sim_by_cell[hit.getCellID()] = hit - + common_cells = sim_by_cell.keys() & digi_by_cell.keys() - + for cellid in common_cells: sim_energy = sim_by_cell[cellid].getEnergy() * 1000.0 digi_energy = digi_by_cell[cellid].getEnergy() * 1000.0 cell_diffs.append(digi_energy - sim_energy) - + n_events += 1 @@ -138,7 +156,7 @@ canvas_transverse.SaveAs(os.path.join(output_folder, "transverse_profiles_layers.pdf")) # --------------------------------- -# Sim-digi cell energy difference +# Sim-digi cell energy difference diff_min = min(cell_diffs) diff_max = max(cell_diffs)