Skip to content

Commit 5994b27

Browse files
author
Michał Sieczczyński
committed
Adding benchmark sampF
1 parent 865731a commit 5994b27

File tree

3 files changed

+274
-0
lines changed

3 files changed

+274
-0
lines changed

‎benchmarks/nhcal_sampling_fraction/Snakefile‎

Lines changed: 56 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,56 @@
1+
import os
2+
3+
DETECTOR_PATH = os.environ["DETECTOR_PATH"]
4+
DETECTOR_CONFIG = "epic_backward_hcal_only.xml"
5+
6+
rule nhcal_sampling_fraction_simulate:
7+
output:
8+
"sim_output/nhcal_sampling_fraction/{PARTICLE}/E{ENERGY}GeV/sim.edm4hep.root",
9+
params:
10+
N_EVENTS=1000,
11+
shell:
12+
"""
13+
set -m
14+
mkdir -p sim_output/nhcal_sampling_fraction/{wildcards.PARTICLE}/E{wildcards.ENERGY}GeV
15+
exec ddsim \
16+
--compactFile {DETECTOR_PATH}/{DETECTOR_CONFIG} \
17+
--numberOfEvents {params.N_EVENTS} \
18+
--random.seed $RANDOM \
19+
--enableGun \
20+
-v WARNING \
21+
--gun.particle {wildcards.PARTICLE} \
22+
--gun.thetaMin 120*degree \
23+
--gun.thetaMax 180*degree \
24+
--gun.distribution uniform \
25+
--gun.energy "{wildcards.ENERGY}*GeV" \
26+
--outputFile {output}
27+
"""
28+
29+
30+
rule nhcal_sampling_fraction_combine:
31+
input:
32+
lambda wildcards: expand(
33+
"sim_output/nhcal_sampling_fraction/{PARTICLE}/E{ENERGY}GeV/sim.edm4hep.root",
34+
ENERGY=["1", "2", "5", "10"],
35+
PARTICLE=["pi-", "neutron", "e-"],
36+
),
37+
wildcard_constraints:
38+
ENERGY=r"\d+"
39+
output:
40+
f"sim_output/nhcal_sampling_fraction/sim_combined.edm4hep.root",
41+
shell:
42+
"""
43+
hadd -f {output} {input}
44+
"""
45+
46+
rule nhcal_sampling_fraction_analysis:
47+
input:
48+
combined="sim_output/nhcal_sampling_fraction/sim_combined.edm4hep.root",
49+
script="benchmarks/nhcal_sampling_fraction/scripts/sampling_fraction_analysis.cxx",
50+
output:
51+
pdf=f"results/nhcal_sampling_fraction/analysis.pdf",
52+
png=f"results/nhcal_sampling_fraction/analysis.png",
53+
shell:
54+
"""
55+
root -l -b -q '{input.script}("{input.combined}","{output.pdf}","{output.png}","{DETECTOR_PATH}/{DETECTOR_CONFIG}")'
56+
"""

‎benchmarks/nhcal_sampling_fraction/config.yml‎

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
sim:nhcal_sampling_fraction:
2+
extends: .det_benchmark
3+
stage: simulate
4+
parallel:
5+
matrix:
6+
- ENERGY: ["1GeV", "2GeV", "5GeV", "10GeV"]
7+
PARTICLE: ["pi-", "neutron", "e-"]
8+
script:
9+
- snakemake --cores 5 sim_output/nhcal_sampling_fraction/${PARTICLE}/E${ENERGY}/sim.edm4hep.root
10+
11+
bench:nhcal_sampling_fraction:
12+
extends: .det_benchmark
13+
stage: benchmarks
14+
needs:
15+
- "sim:nhcal_sampling_fraction"
16+
script:
17+
- snakemake --cores 1 results/nhcal_sampling_fraction/analysis.pdf
18+
19+
collect_results:nhcal_sampling_fraction:
20+
extends: .det_benchmark
21+
stage: collect
22+
needs:
23+
- "bench:nhcal_sampling_fraction"
24+
script:
25+
- ls -lrht
26+
- mv results{,_save}/
27+
- snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output results/nhcal_sampling_fraction/analysis.pdf
28+
- mv results{_save,}/
29+

‎benchmarks/nhcal_sampling_fraction/scripts/sampling_fraction_analysis.cxx‎

Lines changed: 189 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,189 @@
1+
#include <cmath>
2+
#include <fstream>
3+
#include <iostream>
4+
#include <string>
5+
#include <vector>
6+
7+
#include "ROOT/RDataFrame.hxx"
8+
#include <TH1D.h>
9+
#include <TH2D.h>
10+
#include <TRandom3.h>
11+
#include <TFile.h>
12+
#include <TChain.h>
13+
#include <TTree.h>
14+
#include <TMath.h>
15+
#include <TVector3.h>
16+
#include <TCanvas.h>
17+
18+
#include "TROOT.h"
19+
#include "TRandom.h"
20+
#include "TH3.h"
21+
22+
23+
#include "DD4hep/Detector.h"
24+
#include "DDRec/CellIDPositionConverter.h"
25+
26+
#include <podio/Frame.h>
27+
#include <podio/CollectionBase.h>
28+
#include "podio/ROOTReader.h"
29+
//#include <podio/ROOTFrameReader.h>
30+
#include "podio/CollectionIDTable.h"
31+
#include "podio/ObjectID.h"
32+
33+
#include "edm4hep/MCParticleCollection.h"
34+
#include "edm4hep/MCParticleCollectionData.h"
35+
#include "edm4hep/MCParticle.h"
36+
#include "edm4hep/MCParticleData.h"
37+
38+
#include "edm4hep/SimCalorimeterHitCollectionData.h"
39+
#include "edm4hep/SimCalorimeterHitCollection.h"
40+
#include "edm4hep/SimCalorimeterHitData.h"
41+
#include "edm4hep/SimCalorimeterHit.h"
42+
43+
#include "edm4hep/CalorimeterHit.h"
44+
#include "edm4hep/CalorimeterHitCollectionData.h"
45+
#include "edm4hep/CalorimeterHitCollection.h"
46+
#include "edm4hep/CalorimeterHitData.h"
47+
#include "edm4hep/CalorimeterHit.h"
48+
#include "edm4hep/CalorimeterHitObj.h"
49+
50+
#include "edm4eic/ClusterCollection.h"
51+
#include "edm4eic/Cluster.h"
52+
#include "edm4eic/ClusterData.h"
53+
54+
#include "edm4eic/CalorimeterHit.h"
55+
#include "edm4eic/CalorimeterHitCollectionData.h"
56+
#include "edm4eic/CalorimeterHitCollection.h"
57+
#include "edm4eic/CalorimeterHitData.h"
58+
#include "edm4eic/CalorimeterHit.h"
59+
#include "edm4eic/CalorimeterHitObj.h"
60+
61+
#include <edm4eic/vector_utils_legacy.h>
62+
#include <edm4hep/Vector3f.h>
63+
64+
using namespace std;
65+
using namespace ROOT;
66+
using namespace TMath;
67+
using namespace edm4hep;
68+
69+
dd4hep::Detector* det = NULL;
70+
71+
int sampling_fraction_analysis(const string &filename, string outname_pdf, string outname_png, TString compact_file)
72+
{
73+
74+
podio::ROOTReader *reader = new podio::ROOTReader();
75+
reader->openFile(filename);
76+
unsigned nEvents = reader->getEntries("events");
77+
cout << "Number of events: " << nEvents << endl;
78+
79+
det = &(dd4hep::Detector::getInstance());
80+
det->fromCompact(compact_file.Data());
81+
det->volumeManager();
82+
det->apply("DD4hepVolumeManager", 0, 0);
83+
84+
dd4hep::rec::CellIDPositionConverter cellid_converter(*det);
85+
auto idSpec = det->readout("HcalEndcapNHits").idSpec();
86+
auto decoder = idSpec.decoder();
87+
const int slice_index = decoder->index("slice");
88+
if (slice_index < 0) {
89+
cerr << "ERROR: 'slice' field not found in cell ID spec!" << endl;
90+
return 1;
91+
}
92+
93+
TH2D *h_sampF_e = new TH2D("h_sampF_e", "nHCal sampling fraction vs. energy (e-); E [GeV]; sampling fraction; counts",
94+
50, 0.0, 12.0, 50, 0.0, 1.0);
95+
TProfile *p_sampF_e = new TProfile("p_sampF_e", "nHCal sampling fraction vs. energy (e-); E [GeV]; sampling fraction",
96+
50, 0.0, 12.0, 0.0, 1.0);
97+
98+
TH2D *h_sampF_n = new TH2D("h_sampF_n", "nHCal sampling fraction vs. energy (neutron); E [GeV]; sampling fraction; counts",
99+
50, 0.0, 12.0, 50, 0.0, 1.0);
100+
TProfile *p_sampF_n = new TProfile("p_sampF_n", "nHCal sampling fraction vs. energy (neutron); E [GeV]; sampling fraction",
101+
50, 0.0, 12.0, 0.0, 1.0);
102+
103+
TH2D *h_sampF_pi = new TH2D("h_sampF_pi", "nHCal sampling fraction vs. energy (#pi-); E [GeV]; sampling fraction; counts",
104+
50, 0.0, 12.0, 50, 0.0, 1.0);
105+
TProfile *p_sampF_pi = new TProfile("p_sampF_pi", "nHCal sampling fraction vs. energy (#pi-); E [GeV]; sampling fraction",
106+
50, 0.0, 12.0, 0.0, 1.0);
107+
108+
109+
110+
for (unsigned ev = 0; ev < nEvents; ev++)
111+
{
112+
double hit_Esum = 0;
113+
double hit_scint_Esum = 0;
114+
double singlePart_Ekin = 0;
115+
116+
auto frameData = reader->readNextEntry(podio::Category::Event);
117+
if (!frameData)
118+
{
119+
cerr << "Invalid FrameData at event " << ev << endl;
120+
continue;
121+
}
122+
123+
podio::Frame frame(std::move(frameData));
124+
125+
auto& MCParticles_coll = frame.get<edm4hep::MCParticleCollection>("MCParticles");
126+
auto& SimCalorimeterHit_coll = frame.get<edm4hep::SimCalorimeterHitCollection>("HcalEndcapNHits");
127+
128+
if (!SimCalorimeterHit_coll.isValid())
129+
{
130+
cerr << "HcalEndcapNHits collection is not valid!" << endl;
131+
}
132+
if (!MCParticles_coll.isValid())
133+
{
134+
cerr << "MCParticles collection is not valid!" << endl;
135+
}
136+
137+
edm4hep::MCParticle mcpart = MCParticles_coll.at(0);
138+
singlePart_Ekin = mcpart.getEnergy(); //-mcpart.getMass();
139+
int pdg = mcpart.getPDG();
140+
141+
for (const auto& hit : SimCalorimeterHit_coll)
142+
{
143+
const int hit_slice = decoder->get(hit.getCellID(), slice_index);
144+
hit_Esum += hit.getEnergy();
145+
if(hit_slice == 3) hit_scint_Esum += hit.getEnergy();
146+
}
147+
148+
if (pdg == 11) // e-
149+
{
150+
h_sampF_e->Fill(singlePart_Ekin, hit_scint_Esum/hit_Esum);
151+
p_sampF_e->Fill(singlePart_Ekin, hit_scint_Esum/hit_Esum);
152+
}
153+
else if (pdg == -211) // pi-
154+
{
155+
h_sampF_pi->Fill(singlePart_Ekin, hit_scint_Esum/hit_Esum);
156+
p_sampF_pi->Fill(singlePart_Ekin, hit_scint_Esum/hit_Esum);
157+
}
158+
else if (pdg == 2112) // neutron
159+
{
160+
h_sampF_n->Fill(singlePart_Ekin, hit_scint_Esum/hit_Esum);
161+
p_sampF_n->Fill(singlePart_Ekin, hit_scint_Esum/hit_Esum);
162+
}
163+
}
164+
165+
delete reader;
166+
167+
TProfile *p_e_over_pi = (TProfile*) p_sampF_e->Clone("p_e_over_pi");
168+
p_e_over_pi->SetTitle("e/h ratio;E [GeV];e/h");
169+
p_e_over_pi->Divide(p_sampF_pi);
170+
171+
TCanvas *canvas = new TCanvas("canvas", "canvas", 1600, 800);
172+
canvas->Divide(2,2);
173+
canvas->cd(1);
174+
h_sampF_e->Draw("COLZ");
175+
p_sampF_e->Draw("SAME");
176+
canvas->cd(2);
177+
h_sampF_pi->Draw("COLZ");
178+
p_sampF_pi->Draw("SAME");
179+
canvas->cd(3);
180+
h_sampF_n->Draw("COLZ");
181+
p_sampF_n->Draw("SAME");
182+
canvas->cd(4);
183+
p_e_over_pi->Draw("HIST");
184+
185+
canvas->SaveAs(outname_pdf.c_str());
186+
canvas->SaveAs(outname_png.c_str());
187+
188+
return 0;
189+
}

0 commit comments

Comments
 (0)