Skip to content

Commit

Permalink
feat: adding some functionality to material mapping checks (acts-proj…
Browse files Browse the repository at this point in the history
…ect#3266)

This PR adds some (very minor) updates to the Material mapping and
associated tests:

- it adds "mat_r" such r:z plots can be made relatively quickly
- it adds intersection position to allow projection plots
- it adds 2 convenience python files to invest the output

ROOT hashes will change for the material mapping, will update.
  • Loading branch information
asalzburger authored and paulgessinger committed Jun 14, 2024
1 parent 6f5eb85 commit 84e7d1b
Show file tree
Hide file tree
Showing 7 changed files with 474 additions and 11 deletions.
2 changes: 1 addition & 1 deletion Core/src/Material/BinnedSurfaceMaterialAccumulater.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -113,7 +113,7 @@ void Acts::BinnedSurfaceMaterialAccumulater::accumulate(
"Surface material is not found, inconsistent configuration.");
}
// Accumulate the material - remember the touched bin
auto tBin = accMaterial->second.accumulate(mi.position, mi.materialSlab,
auto tBin = accMaterial->second.accumulate(mi.intersection, mi.materialSlab,
mi.pathCorrection);
touchedMapBins.insert(MapBin(&(accMaterial->second), tBin));
}
Expand Down
5 changes: 3 additions & 2 deletions Core/src/Material/MaterialInteractionAssignment.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -81,12 +81,13 @@ Acts::MaterialInteractionAssignment::assign(
}
}

// Assign the material interaction
// Assign the material interaction - position stays, but intersection is
// updated
MaterialInteraction assignedMaterialInteraction = materialInteraction;
assignedMaterialInteraction.pathCorrection = pathCorrection;
assignedMaterialInteraction.surface = surface;
assignedMaterialInteraction.position = position;
assignedMaterialInteraction.direction = direction;
assignedMaterialInteraction.intersection = position;
assignedMaterialInteraction.intersectionID = intersectionID;
// Check for possible reassignment
if (is + 1u < intersectedSurfaces.size() &&
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -146,6 +146,8 @@ class RootMaterialTrackWriter
std::vector<float> m_step_y;
/// step z position
std::vector<float> m_step_z;
/// step r position
std::vector<float> m_step_r;
/// step x (end) position (optional)
std::vector<float> m_step_ex;
/// step y (end) position (optional)
Expand Down Expand Up @@ -175,12 +177,16 @@ class RootMaterialTrackWriter
std::vector<std::uint64_t> m_sur_id;
/// Type of the surface associated with the step
std::vector<int32_t> m_sur_type;
/// x position of the center of the surface associated with the step
/// x position of the surface intersection associated with the step
std::vector<float> m_sur_x;
/// y position of the center of the surface associated with the step
/// y position of the surface intersection associated with the step
std::vector<float> m_sur_y;
/// z position of the center of the surface associated with the step
/// z position of the surface intersection associated with the step
std::vector<float> m_sur_z;
/// r of the position of the surface intersection associated with the step
std::vector<float> m_sur_r;
/// the distance to the surface associated with the step
std::vector<float> m_sur_distance;
/// path correction when associating material to the given surface
std::vector<float> m_sur_pathCorrection;
/// Min range of the surface associated with the step
Expand Down
13 changes: 13 additions & 0 deletions Examples/Io/Root/src/RootMaterialTrackWriter.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@
#include "Acts/Surfaces/SurfaceBounds.hpp"
#include "Acts/Utilities/Intersection.hpp"
#include "Acts/Utilities/Logger.hpp"
#include "Acts/Utilities/VectorHelpers.hpp"
#include "ActsExamples/Framework/AlgorithmContext.hpp"

#include <algorithm>
Expand Down Expand Up @@ -76,6 +77,7 @@ RootMaterialTrackWriter::RootMaterialTrackWriter(
m_outputTree->Branch("mat_x", &m_step_x);
m_outputTree->Branch("mat_y", &m_step_y);
m_outputTree->Branch("mat_z", &m_step_z);
m_outputTree->Branch("mat_r", &m_step_r);
m_outputTree->Branch("mat_dx", &m_step_dx);
m_outputTree->Branch("mat_dy", &m_step_dy);
m_outputTree->Branch("mat_dz", &m_step_dz);
Expand All @@ -100,6 +102,8 @@ RootMaterialTrackWriter::RootMaterialTrackWriter(
m_outputTree->Branch("sur_x", &m_sur_x);
m_outputTree->Branch("sur_y", &m_sur_y);
m_outputTree->Branch("sur_z", &m_sur_z);
m_outputTree->Branch("sur_r", &m_sur_r);
m_outputTree->Branch("sur_distance", &m_sur_distance);
m_outputTree->Branch("sur_pathCorrection", &m_sur_pathCorrection);
m_outputTree->Branch("sur_range_min", &m_sur_range_min);
m_outputTree->Branch("sur_range_max", &m_sur_range_max);
Expand Down Expand Up @@ -143,6 +147,7 @@ ProcessCode RootMaterialTrackWriter::writeT(
m_step_x.clear();
m_step_y.clear();
m_step_z.clear();
m_step_r.clear();
m_step_ex.clear();
m_step_ey.clear();
m_step_ez.clear();
Expand All @@ -161,6 +166,8 @@ ProcessCode RootMaterialTrackWriter::writeT(
m_sur_x.clear();
m_sur_y.clear();
m_sur_z.clear();
m_sur_r.clear();
m_sur_distance.clear();
m_sur_pathCorrection.clear();
m_sur_range_min.clear();
m_sur_range_max.clear();
Expand Down Expand Up @@ -210,6 +217,7 @@ ProcessCode RootMaterialTrackWriter::writeT(
m_step_x.reserve(mints);
m_step_y.reserve(mints);
m_step_z.reserve(mints);
m_step_r.reserve(mints);
m_step_ex.reserve(mints);
m_step_ey.reserve(mints);
m_step_ez.reserve(mints);
Expand All @@ -228,6 +236,8 @@ ProcessCode RootMaterialTrackWriter::writeT(
m_sur_x.reserve(mints);
m_sur_y.reserve(mints);
m_sur_z.reserve(mints);
m_sur_r.reserve(mints);
m_sur_distance.reserve(mints);
m_sur_pathCorrection.reserve(mints);
m_sur_range_min.reserve(mints);
m_sur_range_max.reserve(mints);
Expand Down Expand Up @@ -261,6 +271,7 @@ ProcessCode RootMaterialTrackWriter::writeT(
m_step_x.push_back(mint.position.x());
m_step_y.push_back(mint.position.y());
m_step_z.push_back(mint.position.z());
m_step_r.push_back(perp(mint.position));
m_step_dx.push_back(direction.x());
m_step_dy.push_back(direction.y());
m_step_dz.push_back(direction.z());
Expand Down Expand Up @@ -288,6 +299,8 @@ ProcessCode RootMaterialTrackWriter::writeT(
m_sur_x.push_back(mint.intersection.x());
m_sur_y.push_back(mint.intersection.y());
m_sur_z.push_back(mint.intersection.z());
m_sur_r.push_back(perp(mint.intersection));
m_sur_distance.push_back((mint.position - mint.intersection).norm());
} else if (surface != nullptr) {
auto sfIntersection =
surface
Expand Down
10 changes: 5 additions & 5 deletions Examples/Python/tests/root_file_hashes.txt
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ test_itk_seeding__performance_seeding.root: 78ebda54cd0f026ba4b7f316724ffd946de5
test_itk_seeding__particles.root: 0b6f4ad438010ac48803d48ed98e80b5e87d310bae6c2c02b16cd94d7a4d7d07
test_itk_seeding__particles_simulation.root: ef0246069aa697019f28a8b270a68de95312cae5f2f2c74848566c3ce4f70363
test_propagation__propagation_steps.root: ba2e20d66f9f58850a9248bfaaafecbf0e8257020bb879b4572b783917f4d19b
test_material_recording__geant4_material_tracks.root: 95dc92ab12a4389c3839f03878f366c45a0c4c0035a9f35478980da98ae0f6dd
test_material_recording__geant4_material_tracks.root: c022b9362249b29f57a07926b20644e3ab4ab8ebcf03f773fbf46c446fc1a0a1
test_truth_tracking_kalman[generic-0.0]__trackstates_fitter.root: 26640d70f7fab870e59666b0915569a6f8f82af68c63e5505548ffa9d24a3212
test_truth_tracking_kalman[generic-0.0]__tracksummary_fitter.root: 11b2e2a50343c636fa977175a30220805412d3200e164ae4c3f439fe2087fb88
test_truth_tracking_kalman[generic-0.0]__performance_track_finder.root: ada9cda2ec3c648b144bdaa081d6eff923c80f3d484c4196006e847428cf67a8
Expand All @@ -34,10 +34,10 @@ test_truth_tracking_gsf[generic]__tracksummary_gsf.root: 618a47def602f4a9b387ed2
test_truth_tracking_gsf[odd]__trackstates_gsf.root: 1b70e6e24af0c18f80d42c39e05dd20189bf4d1fcb48851f807cf723f461d7c0
test_truth_tracking_gsf[odd]__tracksummary_gsf.root: ee806a7cddbc4bd41a74dceead569762954e8161af221ad60d637c6264e41752
test_particle_gun__particles.root: 5fe7dda2933ee6b9615b064d192322fe07831133cd998e5ed99a3b992b713a10
test_material_mapping__material-map_tracks.root: b1138566d8d51579dce07f80166f05986942cf78c76b36875aea9b4b8b9bb957
test_material_mapping__propagation-material.root: 2ec28ec8c10860ce0ca03b0b8abbfc9b52cef779d0113fddfe05949a8c362ec0
test_volume_material_mapping__material-map-volume_tracks.root: 5092bbb52f15e4db7466e5a58c515794073e460af2411be2aa4ff001367928e4
test_volume_material_mapping__propagation-volume-material.root: 89894f11a436038f6511d8b126f11ca2899d18b2b048a97c1b6c296ed85faddf
test_material_mapping__material-map_tracks.root: 938b1a855369e9304401cb10d2751df3fd7acf32a2573f2057eb1691cd94edf3
test_material_mapping__propagation-material.root: d081c2c54403941bc7876fa552425eef29b496883dcdaab01113462895d0905e
test_volume_material_mapping__material-map-volume_tracks.root: 98e212d32ca054fa3d01af4167c1f49755a139d43b82c57908197f5985e0a4ff
test_volume_material_mapping__propagation-volume-material.root: cc3a2a6df914526be7e4a0c33db615a9589441b6bfb763271ab7ce7ea7f8789e
test_digitization_example[smeared]__measurements.root: 03184f20b7b8e69f40424f76e2fd0620fb9f319d23ad69aa4bdd25c1b5cd482e
test_digitization_example[geometric]__measurements.root: bc5d158478d7fdc4e09728cc9b4609afaa62be9dfcc1162685d03160f4d89495
test_digitization_example_input[smeared]__particles.root: 5fe7dda2933ee6b9615b064d192322fe07831133cd998e5ed99a3b992b713a10
Expand Down
239 changes: 239 additions & 0 deletions Examples/Scripts/MaterialMapping/material_comparison.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,239 @@
import ROOT
import argparse
import math

from ROOT import (
TCanvas,
TFile,
TTree,
gDirectory,
gStyle,
TH1D,
TH2F,
TProfile,
TRatioPlot,
)


def TH1D_from_TProf(tprof):
h = TH1D(
tprof.GetName() + "_th1",
tprof.GetTitle(),
tprof.GetNbinsX(),
tprof.GetXaxis().GetXmin(),
tprof.GetXaxis().GetXmax(),
)
for i in range(tprof.GetNbinsX()):
if tprof.GetBinContent(i + 1) == 0.0:
h.SetBinContent(i + 1, 0.0)
h.SetBinError(i + 1, 1000.0)
continue
h.SetBinContent(i + 1, tprof.GetBinContent(i + 1))
h.SetBinError(i + 1, tprof.GetBinError(i + 1))
return h


if "__main__" == __name__:
p = argparse.ArgumentParser()

p.add_argument(
"-e", "--entries", type=int, default=100000, help="Number of events to process"
)
p.add_argument(
"-i",
"--input",
type=str,
nargs="+",
default="",
help="Input files with material tracks",
)
p.add_argument(
"-o",
"--output",
type=str,
default="",
help="Output file with produced material overview plots",
)
p.add_argument(
"-l",
"--labels",
type=str,
nargs="+",
default="",
help="The labels for the input files",
)
p.add_argument(
"-v",
"--variables",
type=str,
nargs="+",
default=["t_X0", "t_L0"],
help="The variables to be plotted",
)
p.add_argument(
"-m",
"--max",
type=float,
nargs="+",
default=[4.5, 1.5],
help="The variables to be plotted",
)
p.add_argument(
"-a",
"--axes",
type=str,
nargs="+",
default=["v_eta", "v_phi"],
help="The axes versus which the variables to be plotted",
)
p.add_argument(
"--eta",
type=float,
nargs=2,
default=[-3.0, 3.0],
help="Eta range for the plotting",
)
p.add_argument("--eta-bins", type=int, default=60, help="Eta bins for the plotting")
p.add_argument(
"--phi",
type=float,
nargs=2,
default=[-math.pi, math.pi],
help="Phi range for the plotting",
)
p.add_argument("--phi-bins", type=int, default=72, help="Phi bins for the plotting")
p.add_argument(
"-d",
"--detray-input",
type=str,
default="",
help="Optional detray input csv file",
)

args = p.parse_args()

ttrees = []

if len(args.input) != len(args.labels):
print("** ERROR ** The number of input files and labels must match")
exit(1)

# The histogram map
h_dict = {}
tfiles = []

ofile = ROOT.TFile.Open(args.output, "RECREATE") if args.output else None

input_files = args.input

# Detray csv file
if args.detray_input:
input_files.append(args.detray_input)

# Loop over the files and create the comparison histograms
for fi, ifile in enumerate(input_files):
# Special treament for detray input
if ifile == args.detray_input:
print("Reading detray input from file: ", args.detray_input)
ttree = TTree("t", "material read from detray")
ttree.ReadFile(args.detray_input, "v_eta/D:v_phi:t_X0:t_L0:t_X0:t_L0")
lbl = "detray"
else:
# Get the file
tfile = ROOT.TFile.Open(ifile)
tfiles.append(tfile)
ttree = tfile.Get("material-tracks")
ttrees.append(ttree)
# The label
lbl = args.labels[fi]

# Loop over the variables and axes
for iv, var in enumerate(args.variables):
for ax in args.axes:
if ax == "v_eta":
bins = args.eta_bins
low = args.eta[0]
high = args.eta[1]
cut = f"v_phi > {args.phi[0]} && v_phi < {args.phi[1]}"
elif ax == "v_phi":
bins = args.phi_bins
low = args.phi[0]
high = args.phi[1]
cut = f"v_eta > {args.eta[0]} && v_eta < {args.eta[1]}"
# Some naming magic
hcmmd = f"{var}:{ax}>>"
hname = f"{var}_vs_{ax}"
hfname = f"{lbl}_{var}_vs_{ax}"
hrange = f"({bins},{low},{high},100,0,{args.max[iv]})"
htitle = f"{var} vs {ax}"
ttree.Draw(hcmmd + hfname + hrange, cut, "")
h = ROOT.gDirectory.Get(hfname)
# Fill into comparison histogram
if h_dict.get(hname):
h_dict[hname].append(h)
else:
h_dict[hname] = [h]

# Write to file
if ofile:
ofile.cd()
h.Write()

colors = [
ROOT.kBlack,
ROOT.kRed + 2,
ROOT.kBlue - 4,
ROOT.kGreen + 1,
ROOT.kYellow - 2,
]
markers = [
ROOT.kFullCircle,
ROOT.kFullCircle,
ROOT.kFullTriangleUp,
ROOT.kFullTriangleDown,
ROOT.kFullDiamond,
]

# Now create the comparison histograms
c = ROOT.TCanvas("Comparison", "Comparison", 1200, 800)
c.Divide(2, 2)

# Remove the stat box
gStyle.SetOptStat(0)

# Memory garbage collection, thanks ROOT
hist_memory_pool = []

ic = 0
for hname in h_dict:
ic += 1
c.cd(ic)
h_list = h_dict[hname]
h_ref = None
h_prof_ref = None
for ih, h in enumerate(h_list):
h_prof = TH1D_from_TProf(h.ProfileX())
hist_memory_pool.append(h_prof)
h_prof.SetObjectStat(0)
h_prof.SetLineColor(colors[ih])
h_prof.SetMarkerColor(colors[ih])
h_prof.SetMarkerSize(0.5)
h_prof.SetMarkerStyle(markers[ih])
h_prof.GetYaxis().SetRangeUser(0.0, 1.3 * h_prof.GetMaximum())
if ih == 0:
h_ref = h
h_prof_ref = h_prof
else:
h_ratio = ROOT.TRatioPlot(h_prof_ref, h_prof)
h_ratio.SetGraphDrawOpt("pe")
h_ratio.SetSeparationMargin(0.005)
drawOption = "e,same" if ih > 1 else "e"
h_ratio.Draw(drawOption)
h_ratio.GetLowerRefGraph().SetLineColor(colors[ih])
h_ratio.GetLowerRefGraph().SetMarkerColor(colors[ih])
h_ratio.GetLowerRefGraph().SetMarkerSize(0.5)
h_ratio.GetLowerRefGraph().SetMarkerStyle(markers[ih])
c.Update()
hist_memory_pool.append(h_ratio)
h_prof.Draw(("same" if ih > 0 else ""))
c.SaveAs(f"{hname}.png")
Loading

0 comments on commit 84e7d1b

Please sign in to comment.