diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 036912933ce205529468c311fa8ab0a8983768ab..3b7997b011f59e5abc0f49f2b9ec48f498daaa8f 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -135,6 +135,7 @@ include:
- local: 'benchmarks/barrel_ecal/config.yml'
- local: 'benchmarks/barrel_hcal/config.yml'
- local: 'benchmarks/zdc/config.yml'
+ - local: 'benchmarks/zdc_lyso/config.yml'
- local: 'benchmarks/material_maps/config.yml'
- local: 'benchmarks/material_scan/config.yml'
- local: 'benchmarks/pid/config.yml'
diff --git a/Snakefile b/Snakefile
index e7a5435022346663a04feda247d73313795d8e99..26cfc5d0427ea6f40f10cedbe5425a30822dc1fd 100644
--- a/Snakefile
+++ b/Snakefile
@@ -3,7 +3,7 @@ include: "benchmarks/barrel_ecal/Snakefile"
include: "benchmarks/ecal_gaps/Snakefile"
include: "benchmarks/material_scan/Snakefile"
include: "benchmarks/tracking_performances/Snakefile"
-
+include: "benchmarks/zdc_lyso/Snakefile"
rule fetch_epic:
output:
diff --git a/benchmarks/zdc_lyso/README.md b/benchmarks/zdc_lyso/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..3d53843ee1834f406ec7e3947f2bcac222109ff6
--- /dev/null
+++ b/benchmarks/zdc_lyso/README.md
@@ -0,0 +1,11 @@
+Detector Benchmark for LYSO ZDC
+===============================
+
+## Overview
+This benchmark generates events with single low-energy (5 MeV - 1 GeV) photons. The photons are generated with angles of 0-5.24 mRad with respect to the proton/ion beam direction. The benchmark creates performance plots of the LYSO ZDC acceptance and energy reconstruction resolution.
+
+## Contacts
+[JiaJun Huang](jhuan328@ucr.edu)
+[Barak Schmookler](baraks@ucr.edu)
+
+
diff --git a/benchmarks/zdc_lyso/Snakefile b/benchmarks/zdc_lyso/Snakefile
new file mode 100644
index 0000000000000000000000000000000000000000..5df026dffb99a669f9d8858d93473f4ca1d6b3cc
--- /dev/null
+++ b/benchmarks/zdc_lyso/Snakefile
@@ -0,0 +1,86 @@
+import os
+
+
+rule ecal_lyso_sim_hepmc:
+ input:
+ script = "benchmarks/zdc_lyso/gen_particles.cxx",
+ output:
+ hepmcfile="data/{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc",
+ log:
+ "data/{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc.log",
+ params:
+ num_events=1000,
+ shell:
+ """
+root -l -b -q '{input.script}({params.num_events}, "{output.hepmcfile}", "{wildcards.PARTICLE}", {wildcards.THETA_MIN}, {wildcards.THETA_MAX}, 0, 360, {wildcards.BEAM_ENERGY})'
+"""
+
+
+rule ecal_lyso_sim:
+ input:
+ hepmcfile="data/{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc",
+ warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root",
+ output:
+ "sim_output/zdc_lyso/{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.edm4hep.root",
+ log:
+ "sim_output/zdc_lyso/{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.edm4hep.root.log",
+ params:
+ num_events=1000,
+ shell:
+ """
+npsim \
+ --runType batch \
+ -v WARNING \
+ --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \
+ --numberOfEvents {params.num_events} \
+ --inputFiles {input.hepmcfile} \
+ --outputFile {output}
+"""
+
+
+rule ecal_lyso_reco:
+ input:
+ "sim_output/zdc_lyso/{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.edm4hep.root",
+ output:
+ "sim_output/zdc_lyso/{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.eicrecon.tree.edm4eic.root",
+ log:
+ "sim_output/zdc_lyso/{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.eicrecon.tree.edm4eic.root.log",
+ shell:
+ """
+eicrecon -Ppodio:output_collections=HcalFarForwardZDCRawHits,HcalFarForwardZDCRecHits,HcalFarForwardZDCClusters,EcalFarForwardZDCRawHits,EcalFarForwardZDCRecHits,EcalFarForwardZDCClusters,MCParticles {input}
+mv podio_output.root {output}
+"""
+
+
+rule zdc_analysis:
+ input:
+ expand("sim_output/zdc_lyso/{{DETECTOR_CONFIG}}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.eicrecon.tree.edm4eic.root",
+ PARTICLE=["gamma"],
+ BEAM_ENERGY=["0.005", "0.01", "0.05", "0.1", "0.5", "1.0"],
+ THETA_MIN=["0"],
+ THETA_MAX=["0.3"]),
+ script="benchmarks/zdc_lyso/analysis/analysis.py",
+ output:
+ "results/{DETECTOR_CONFIG}/zdc_lyso/plots.pdf",
+ shell:
+ """
+python {input.script}
+"""
+
+
+# Examples of invocation
+rule create_all_hepmc:
+ input:
+ expand("data/{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc",
+ PARTICLE=["gamma"],
+ BEAM_ENERGY=["0.005", "0.01", "0.05", "0.1", "0.5", "1.0"],
+ THETA_MIN=["0"],
+ THETA_MAX=["0.3"])
+
+
+rule run_all_locally:
+ input:
+ "results/" + os.environ["DETECTOR_CONFIG"] + "/zdc_lyso/plots.pdf",
+ message:
+ "See output in {input[0]}"
+
diff --git a/benchmarks/zdc_lyso/analysis/analysis.py b/benchmarks/zdc_lyso/analysis/analysis.py
new file mode 100644
index 0000000000000000000000000000000000000000..04bd233fd43f540fc09405a0ed544aad173cbb67
--- /dev/null
+++ b/benchmarks/zdc_lyso/analysis/analysis.py
@@ -0,0 +1,195 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import mplhep as hep
+import uproot
+import pandas as pd
+import os
+from scipy.optimize import curve_fit
+from matplotlib.backends.backend_pdf import PdfPages
+
+plt.figure()
+hep.set_style(hep.style.CMS)
+hep.set_style("CMS")
+
+def gaussian(x, amp, mean, sigma):
+ return amp * np.exp( -(x - mean)**2 / (2*sigma**2) )
+
+def rotateY(xdata, zdata, angle):
+ s = np.sin(angle)
+ c = np.cos(angle)
+ rotatedz = c*zdata - s*xdata
+ rotatedx = s*zdata + c*xdata
+ return rotatedx, rotatedz
+
+Energy = [0.005, 0.01, 0.05, 0.1, 0.5, 1.0]
+q0 = [5, 10, 40, 90, 400, 700]
+q1 = [0.5, 0.5, 0.9, 5, 10, 20]
+
+df = pd.DataFrame({})
+for eng in Energy:
+ tree = uproot.open(f'sim_output/zdc_lyso/{os.environ["DETECTOR_CONFIG"]}_gamma_{eng}GeV_theta_0deg_thru_0.3deg.eicrecon.tree.edm4eic.root')['events']
+ ecal_reco_energy = tree['EcalFarForwardZDCClusters/EcalFarForwardZDCClusters.energy'].array()
+ #hcal_reco_energy = tree['HcalFarForwardZDCClusters/HcalFarForwardZDCClusters.energy'].array()
+
+ tree = uproot.open(f'sim_output/zdc_lyso/{os.environ["DETECTOR_CONFIG"]}_gamma_{eng}GeV_theta_0deg_thru_0.3deg.edm4hep.root')['events']
+ ecal_sim_energy = tree['EcalFarForwardZDCHits/EcalFarForwardZDCHits.energy'].array()
+ hcal_sim_energy = tree['HcalFarForwardZDCHits/HcalFarForwardZDCHits.energy'].array()
+
+ par_x = tree['MCParticles/MCParticles.momentum.x'].array()[:,2]
+ par_y = tree['MCParticles/MCParticles.momentum.y'].array()[:,2]
+ par_z = tree['MCParticles/MCParticles.momentum.z'].array()[:,2]
+
+ eng = int(eng*1000)
+
+ ecal_reco_energy = pd.DataFrame({f'ecal_reco_energy_{eng}': np.array(ecal_reco_energy.tolist(), dtype=object)})
+ #hcal_reco_energy = pd.DataFrame({f'hcal_reco_energy_{eng}': np.array(hcal_reco_energy.tolist(), dtype=object)})
+ ecal_sim_energy = pd.DataFrame({f'ecal_sim_energy_{eng}': np.array(ecal_sim_energy.tolist(), dtype=object)})
+ hcal_sim_energy = pd.DataFrame({f'hcal_sim_energy_{eng}': np.array(hcal_sim_energy.tolist(), dtype=object)})
+
+ par_x = pd.DataFrame({f'par_x_{eng}': np.array(par_x.tolist(), dtype=object)})
+ par_y = pd.DataFrame({f'par_y_{eng}': np.array(par_y.tolist(), dtype=object)})
+ par_z = pd.DataFrame({f'par_z_{eng}': np.array(par_z.tolist(), dtype=object)})
+
+ df = pd.concat([df,ecal_reco_energy,ecal_sim_energy,hcal_sim_energy,par_x,par_y,par_z],axis=1)
+
+mu = []
+sigma = []
+resolution = []
+
+fig1, ax = plt.subplots(3,2,figsize=(20,10))
+#fig1.suptitle('ZDC ECal Cluster Energy Reconstruction')
+
+plt.tight_layout()
+for i in range(6):
+ plt.sca(ax[i%3,i//3])
+ eng = int(Energy[i]*1000)
+ plt.title(f'Gamma Energy: {eng} MeV')
+ temp = np.array([sum(item) if (item != 0) else 0 for item in df[f'ecal_reco_energy_{eng}']])
+ hist, x = np.histogram(np.array(temp)*1000,bins=30)
+ x = x[1:]/2 + x[:-1]/2
+ plt.errorbar(x,hist,yerr=np.sqrt(hist),fmt='-o')
+ temp = np.array([item[0] for item in df[f'ecal_reco_energy_{eng}'] if item])
+ hist, x = np.histogram(np.array(temp)*1000,bins=30)
+ x = x[1:]/2 + x[:-1]/2
+ coeff, covar = curve_fit(gaussian,x,hist,p0=(200,q0[i],q1[i]),maxfev = 80000)
+ plt.plot(np.linspace(coeff[1]-5*coeff[2],coeff[1]+5*coeff[2],50),gaussian(np.linspace(coeff[1]-5*coeff[2],coeff[1]+5*coeff[2],50),*coeff)
+ ,label = f'$\mu$ = {coeff[1]:.3f} $\pm$ {covar[1][1]:.3f}\n$\sigma$ = {np.abs(coeff[2]):.3f} $\pm$ {covar[2][2]:.3f}\nResolution = {np.abs(coeff[2])*100/coeff[1]:.2f}%')
+
+ plt.xlabel('Energy (MeV)')
+ plt.legend()
+ mu.append(coeff[1])
+ sigma.append(coeff[2])
+ resolution.append(np.abs(coeff[2])*100/coeff[1])
+#plt.savefig('results/Energy_reconstruction_cluster.pdf')
+#plt.show()
+
+fig2, (ax1,ax2) = plt.subplots(2,1,figsize=(15,10),sharex=True)
+
+plt.tight_layout()
+# Plot data on primary axis
+ax1.scatter(np.array(Energy)*1000, mu, c='b')
+ax1.plot([4.5,1000],[4.5,1000],c='black',label='x=y')
+ax1.set_ylabel('Reconstructed Energy (MeV)')
+ax1.set_yscale('log')
+ax1.legend()
+ax1.set_title('ECal Craterlake Cluster Energy Reconstruction')
+
+ax2.plot(np.array(Energy)*1000, resolution, c='r')
+ax2.scatter(np.array(Energy)*1000, resolution, c='r')
+ax2.set_ylabel('Resolution (%)')
+ax2.set_xlabel('Gamma Energy (MeV)')
+ax2.set_xscale('log')
+#plt.savefig('results/Energy_resolution.pdf')
+#plt.show()
+
+htower = []
+herr = []
+hmean = []
+hhits = []
+hhits_cut = []
+emean = []
+ehits = []
+etower = []
+eerr = []
+ehits_cut = []
+
+fig3, ax = plt.subplots(2,3,figsize=(20,10))
+fig3.suptitle('ZDC Simulation Energy Reconstruction')
+for i in range(6):
+ plt.sca(ax[i//3,i%3])
+ eng = int(Energy[i]*1000)
+
+ x = df[f'par_x_{eng}'].astype(float).to_numpy()
+ y = df[f'par_y_{eng}'].astype(float).to_numpy()
+ z = df[f'par_z_{eng}'].astype(float).to_numpy()
+ x, z = rotateY(x,z, 0.025)
+ theta = np.arccos(z/np.sqrt((x**2+y**2+z**2)))*1000
+ condition = theta <= 3.5
+
+ plt.title(f'Gamma Energy: {eng} MeV')
+ energy1 = np.array([sum(item) if (item != 0) else 0 for item in df[f'hcal_sim_energy_{eng}']])#df.eval(f'hcal_sim_energy_{eng}').apply(lambda row: sum(row))
+ hist, x = np.histogram(energy1*1000,bins=np.logspace(0,3,200))
+ x = x[1:]/2 + x[:-1]/2
+ plt.plot(x,hist,marker='o',label="HCal")
+ hhits.append(len(energy1[energy1!=0]))
+ condition1 = energy1!=0
+ hhits_cut.append(len(energy1[condition & condition1])/len(condition[condition==True]))
+ energy = np.array([sum(item) if (item != 0) else 0 for item in df[f'ecal_sim_energy_{eng}']])#df.eval(f'ecal_sim_energy_{eng}').apply(lambda row: sum(row))
+ hist, x = np.histogram(energy*1000,bins=np.logspace(0,3,200))
+ x = x[1:]/2 + x[:-1]/2
+ plt.plot(x,hist,marker='o',label="ECal")
+ emean.append(sum(energy[energy!=0])*1000/len(energy[energy!=0]))
+ hmean.append(sum(energy1[energy!=0])*1000/len(energy[energy!=0]))
+ condition1 = energy!=0
+ ehits_cut.append(len(energy[condition & condition1])/len(condition[condition==True]))
+ ehits.append(len(energy[energy!=0]))
+ plt.legend()
+ plt.xscale('log')
+ plt.xlim(1e0,1e3)
+
+plt.xlabel('Energy (MeV)')
+#plt.savefig('results/Energy_deposition.pdf')
+#plt.show()
+
+fig4, ax = plt.subplots(2,1,sharex=True,gridspec_kw={'height_ratios': [2,1]})
+plt.sca(ax[0])
+plt.errorbar(np.array(Energy)*1000,np.array(hmean)*47.619+np.array(emean),label='HCal+ECal',fmt='-o')
+plt.errorbar(np.array(Energy)*1000,emean,label='ECal',fmt='-o')
+plt.legend()
+plt.yscale('log')
+plt.xscale('log')
+plt.ylabel('Simulation Energy (MeV)')
+plt.sca(ax[1])
+plt.errorbar(np.array(Energy)*1000,(1 - np.array(emean)/(np.array(hmean)*47.619+np.array(emean)))*100,label='Total/ECal',fmt='-o')
+plt.legend()
+plt.ylabel('Fraction of energy\n deposited in Hcal (%)')
+plt.xlabel('Truth Energy (MeV)')
+#plt.savefig('results/Energy_ratio_and_Leakage.pdf')
+plt.tight_layout()
+#plt.show()
+
+fig5 = plt.figure()
+plt.errorbar(np.array(Energy)*1000,np.array(hhits)/1000*100,label='HCal Hits',fmt='-o')
+plt.errorbar(np.array(Energy)*1000,np.array(ehits)/1000*100,label='ECal Hits',fmt='-o')
+#plt.errorbar(np.array(Energy)*1000,np.array(hhits)/np.array(ehits)*100,label='HCal / ECal',fmt='-o',c='b')
+
+plt.errorbar(np.array(Energy)*1000,np.array(hhits_cut)*100,label='HCal Hits with 3.5 mRad cut',fmt='-^')
+plt.errorbar(np.array(Energy)*1000,np.array(ehits_cut)*100,label='ECal Hits with 3.5 mRad cut',fmt='-^')
+#plt.errorbar(np.array(Energy)*1000,np.array(hhits_cut)/np.array(ehits_cut)*100,label='HCal / ECal with 3.5 mRad cut',fmt='-^',c='b')
+### 3mrad cuts
+
+plt.legend()
+plt.xlabel('Simulation Truth Gamma Energy (MeV)')
+plt.ylabel('Fraction of Events with non-zero energy (%)')
+#plt.savefig('results/Hits.pdf')
+plt.xscale('log')
+#plt.show()
+
+#pdfs = ['results/Energy_reconstruction_cluster.pdf','results/Energy_resolution.pdf','results/Energy_deposition.pdf','results/Energy_ratio_and_Leakage.pdf','results/Hits.pdf']
+with PdfPages(f'results/{os.environ["DETECTOR_CONFIG"]}/zdc_lyso/plots.pdf') as pdf:
+ pdf.savefig(fig1)
+ pdf.savefig(fig2)
+ pdf.savefig(fig3)
+ pdf.savefig(fig4)
+ pdf.savefig(fig5)
+
diff --git a/benchmarks/zdc_lyso/config.yml b/benchmarks/zdc_lyso/config.yml
new file mode 100644
index 0000000000000000000000000000000000000000..7810e0cf5e50e82196685c0e3aeafe913840bcb5
--- /dev/null
+++ b/benchmarks/zdc_lyso/config.yml
@@ -0,0 +1,17 @@
+sim:zdc_lyso:
+ extends: .det_benchmark
+ stage: simulate
+ script:
+ - snakemake --cores 1 run_all_locally
+ retry:
+ max: 2
+ when:
+ - runner_system_failure
+
+collect_results:zdc_lyso:
+ extends: .det_benchmark
+ stage: collect
+ needs:
+ - "sim:zdc_lyso"
+ script:
+ - ls -lrht
diff --git a/benchmarks/zdc_lyso/gen_particles.cxx b/benchmarks/zdc_lyso/gen_particles.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..b84b117096acdd36caa1ea669b0b69d683ecf9bc
--- /dev/null
+++ b/benchmarks/zdc_lyso/gen_particles.cxx
@@ -0,0 +1,126 @@
+#include "HepMC3/GenEvent.h"
+#include "HepMC3/ReaderAscii.h"
+#include "HepMC3/WriterAscii.h"
+#include "HepMC3/Print.h"
+
+#include "TRandom3.h"
+#include "TVector3.h"
+
+#include <iostream>
+#include <random>
+#include <cmath>
+#include <math.h>
+#include <TMath.h>
+#include <TDatabasePDG.h>
+#include <TParticlePDG.h>
+
+using namespace HepMC3;
+
+// Generate single electron as input to the Insert simulation.
+// --
+// We generate events with a constant polar angle with respect to
+// the proton direction and then rotate so that the events are given
+// in normal lab coordinate system
+// --
+void gen_particles(
+ int n_events = 10,
+ const char* out_fname = "gen_particles.hepmc",
+ TString particle_name = "e-",
+ double th_min = 3., // Minimum polar angle, in degrees
+ double th_max = 3., // Maximum polar angle, in degrees
+ double phi_min = 0., // Minimum azimuthal angle, in degrees
+ double phi_max = 360., // Maximum azimuthal angle, in degrees
+ double p = 10., // Momentum in GeV/c
+ int dist = 0 //Momentum distribution: 0=fixed, 1=uniform, 2=Gaussian
+ )
+{
+ WriterAscii hepmc_output(out_fname);
+ int events_parsed = 0;
+ GenEvent evt(Units::GEV, Units::MM);
+
+ // Random number generator
+ TRandom3 *r1 = new TRandom3(0); //Use time as random seed
+
+ // Getting generated particle information
+ TDatabasePDG *pdg = new TDatabasePDG();
+ TParticlePDG *particle = pdg->GetParticle(particle_name);
+ const double mass = particle->Mass();
+ const int pdgID = particle->PdgCode();
+
+ for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
+
+ //Set the event number
+ evt.set_event_number(events_parsed);
+
+ // FourVector(px,py,pz,e,pdgid,status)
+ // type 4 is beam
+ // pdgid 11 - electron
+ // pdgid 111 - pi0
+ // pdgid 2212 - proton
+ GenParticlePtr p1 =
+ std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
+ GenParticlePtr p2 = std::make_shared<GenParticle>(
+ FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
+
+ // Define momentum with respect to proton direction
+ double phi = r1->Uniform(phi_min*TMath::DegToRad(),phi_max*TMath::DegToRad());
+ double th = r1->Uniform(th_min*TMath::DegToRad(),th_max*TMath::DegToRad());
+
+ //Total momentum distribution
+ double pevent = -1;
+ if(dist==0){ //fixed
+ pevent = p;
+ }
+ else if(dist==1){ //Uniform: +-50% variation
+ pevent = p*(1. + r1->Uniform(-0.5,0.5) );
+ }
+ else if(dist==2){ //Gaussian: Sigma = 0.1*mean
+ while(pevent<0) //Avoid negative values
+ pevent = r1->Gaus(p,0.1*p);
+ }
+
+ double px = pevent * std::cos(phi) * std::sin(th);
+ double py = pevent * std::sin(phi) * std::sin(th);
+ double pz = pevent * std::cos(th);
+ TVector3 pvec(px,py,pz);
+
+ //Rotate to lab coordinate system
+ double cross_angle = -25./1000.; //in Rad
+ TVector3 pbeam_dir(sin(cross_angle),0,cos(cross_angle)); //proton beam direction
+ pvec.RotateY(-pbeam_dir.Theta()); // Theta is returned positive, beam in negative X
+ // type 1 is final state
+ // pdgid 11 - electron 0.510 MeV/c^2
+ GenParticlePtr p3 = std::make_shared<GenParticle>(
+ FourVector(
+ pvec.X(), pvec.Y(), pvec.Z(),
+ sqrt(pevent*pevent + (mass * mass))),
+ pdgID, 1);
+
+ //If wanted, set non-zero vertex
+ double vx = 0.;
+ double vy = 0.;
+ double vz = 0.;
+ double vt = 0.;
+
+ GenVertexPtr v1 = std::make_shared<GenVertex>();
+ evt.shift_position_by(FourVector(vx, vy, vz, vt));
+ v1->add_particle_in(p1);
+ v1->add_particle_in(p2);
+
+ v1->add_particle_out(p3);
+ evt.add_vertex(v1);
+
+ if (events_parsed == 0) {
+ std::cout << "First event: " << std::endl;
+ Print::listing(evt);
+ }
+
+ hepmc_output.write_event(evt);
+ if (events_parsed % 100 == 0) {
+ std::cout << "Event: " << events_parsed << std::endl;
+ }
+ evt.clear();
+ }
+ hepmc_output.close();
+ std::cout << "Events parsed and written: " << events_parsed << std::endl;
+}