Added benchmark

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/
 
 rule fetch_epic:
     output:

benchmarks/zdc_lyso/README.md

@@ -0,0 +1,6 @@
+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.
+

benchmarks/zdc_lyso/Snakefile

@@ -0,0 +1,88 @@
+rule all:
+    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."],
+               theta_min=["0"],
+               theta_max=["0.3"],
+	       phi_min=["0"],
+	       phi_max=["360"])
+
+
+rule ecal_lyso_sim_hepmc:
+    input:
+        gen_par = "gen_particles.cxx",
+    output:
+        hepmcfile="data/{particle}_{beam_energy}GeV_theta_{theta_min}deg_thru_{theta_max}deg.hepmc",
+    log:
+        "log/{particle}_{beam_energy}GeV_theta_{theta_min}deg_thru_{theta_max}deg.log",
+    wildcard_constraints:
+
+    params:
+        physics_list="FTFP_BERT",
+        num_events=1000,
+    shell:
+        """
+root -l -b -q '{input.gen_par}({params.num_events}, "{output.hepmcfile}", "{wildcards.particle}", {wildcards.theta_min}, {wildcards.theta_max}, 0, 360, {wildcards.beam_energy})'
+"""
+
+
+rule ecal_lyso_sim:
+    input:
+        "data/{particle}_{beam_energy}GeV_theta_{theta_min}deg_thru_{theta_max}deg.hepmc",
+    output:
+        "data/{DETECTOR}_sim_{particle}_{beam_energy}GeV_theta_{theta_min}deg_thru_{theta_max}deg.edm4hep.root",
+    log:
+        "log/{DETECTOR}_sim_{particle}_{beam_energy}GeV_theta_{theta_min}deg_thru_{theta_max}deg.log",
+    wildcard_constraints:
+
+    params:
+        physics_list="FTFP_BERT",
+        num_events=1000,
+
+    shell:
+        """
+npsim \
+  --compactFile $DETECTOR_PATH/{wildcards.DETECTOR}.xml \
+  --numberOfEvents {params.num_events} \
+  --physicsList {params.physics_list} \
+  --inputFiles {input} \
+  --outputFile {output}  || exit
+"""
+
+
+rule ecal_lyso_reco:
+    input:
+        "data/{DETECTOR}_sim_{particle}_{beam_energy}GeV_theta_{theta_min}deg_thru_{theta_max}deg.edm4hep.root",
+    output:
+        "data/{DETECTOR}_reco_{particle}_{beam_energy}GeV_theta_{theta_min}deg_thru_{theta_max}deg.edm4hep.root",
+    log:
+        "log/{DETECTOR}_reco_{particle}_{beam_energy}GeV_theta_{theta_min}deg_thru_{theta_max}deg.log",
+    wildcard_constraints:
+
+    params:
+        DETECTOR="epic_craterlake",
+    shell:
+        """
+eicrecon -Ppodio:output_include_collections=HcalFarForwardZDCRawHits,HcalFarForwardZDCRecHits,HcalFarForwardZDCClusters,EcalFarForwardZDCRawHits,EcalFarForwardZDCRecHits,EcalFarForwardZDCClusters,MCParticles {input}
+mv podio_output.root {output}
+"""
+
+
+rule zdc_analysis:
+	input:
+		"data/epic_craterlake_reco_gamma_0.005GeV_theta_0deg_thru_0.3deg.edm4hep.root",
+		"data/epic_craterlake_reco_gamma_0.01GeV_theta_0deg_thru_0.3deg.edm4hep.root",
+		"data/epic_craterlake_reco_gamma_0.05GeV_theta_0deg_thru_0.3deg.edm4hep.root",
+		"data/epic_craterlake_reco_gamma_0.1GeV_theta_0deg_thru_0.3deg.edm4hep.root",
+		"data/epic_craterlake_reco_gamma_0.5GeV_theta_0deg_thru_0.3deg.edm4hep.root",
+		"data/epic_craterlake_reco_gamma_1GeV_theta_0deg_thru_0.3deg.edm4hep.root",
+
+		script="analysis.py",
+	output:
+		"results/plots.pdf",
+
+	shell:
+		"""
+python {input.script}
+"""

benchmarks/zdc_lyso/analysis/analysis.py

@@ -0,0 +1,154 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import mplhep as hep
+import uproot
+import pandas as pd
+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) ) 
+
+Energy = [0.005, 0.01, 0.05, 0.1, 0.5, 1]
+
+df = pd.DataFrame({})
+for eng in Energy:
+    tree = uproot.open(f'data/epic_craterlake_reco_gamma_{eng}GeV_theta_0deg_thru_0.3deg.edm4hep.root')['events']
+    ecal_reco_energy = tree['EcalFarForwardZDCClusters/EcalFarForwardZDCClusters.energy'].array()
+    #hcal_reco_energy = tree['HcalFarForwardZDCClusters/HcalFarForwardZDCClusters.energy'].array()
+
+    tree = uproot.open(f'data/epic_craterlake_sim_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()
+
+    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)})
+
+    df = pd.concat([df,ecal_reco_energy,ecal_sim_energy,hcal_sim_energy],axis=1)
+
+
+mu = []
+sigma = []
+resolution = []
+fig1, ax = plt.subplots(3,2,figsize=(20,10))
+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 = [item[0] for item in df[f'ecal_reco_energy_{eng}'] if item]
+    hist, x = np.histogram(temp,bins=15)
+    x = x[1:]/2 + x[:-1]/2
+    plt.scatter(x,hist,marker='o')
+    coeff, covar = curve_fit(gaussian,x,hist,p0=(200,eng,eng/2))
+    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(Energy, mu, c='b')
+ax1.plot([0.0045,1],[0.0045,1],c='black',label='x=y')
+ax1.set_ylabel('Reconstructed Energy (GeV)')
+ax1.legend()
+ax1.set_title('ECal Craterlake Energy Reconstruction')
+
+ax2.plot(Energy, resolution, c='r')
+ax2.scatter(Energy, resolution, c='r')
+ax2.set_ylabel('Resolution (%)')
+ax2.set_xlabel('Gamma Energy (GeV)')
+ax2.set_xscale('log')
+plt.savefig('results/Energy_resolution.pdf')
+
+plt.show()
+
+
+htower = []
+herr = []
+hmean = []
+hhits = []
+emean = []
+ehits = []
+etower = []
+eerr = []
+fig3, ax = plt.subplots(2,3,figsize=(20,10))
+for i in range(6):
+    plt.sca(ax[i//3,i%3])
+
+    eng = int(Energy[i]*1000)
+
+    energy = np.array([sum(item) for item in df[f'hcal_sim_energy_{eng}'] if item])#df.eval(f'hcal_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="HCal")
+    hmean.append(sum(energy)*1000/len(energy))
+    hhits.append(len(energy[energy!=0]))
+    energy = np.array([sum(item) for item in df[f'ecal_sim_energy_{eng}'] if item])#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)*1000/len(energy))
+    ehits.append(len(energy[energy!=0]))
+    plt.legend()
+    plt.xscale('log')
+    plt.xlim(1e0,1e3)
+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,hmean,label='HCal Active Layers',fmt='-o')
+plt.errorbar(np.array(Energy)*1000,np.array(hmean)*47.619,label='HCal',fmt='-o')
+plt.errorbar(np.array(Energy)*1000,np.array(hmean)*47.619+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,np.array(hmean)/np.array(emean),label='HCal Active Layers/ECal',fmt='-o')
+plt.errorbar(np.array(Energy)*1000,np.array(hmean)/np.array(emean)*47.619,label='HCal/ECal',fmt='-o')
+plt.legend()
+plt.ylabel('Leakage Ratio')
+plt.xlabel('Truth Energy (MeV)')
+plt.savefig('results/Energy_ratio_and_Leakage.pdf')
+plt.tight_layout()
+plt.show()
+
+fig5 = plt.figure()
+plt.scatter(Energy,np.array(hhits)/1000*100,label='HCal Hits')
+plt.scatter(Energy,np.array(ehits)/1000*100,label='ECal Hits')
+plt.plot(Energy,np.array(hhits)/1000*100)
+plt.plot(Energy,np.array(ehits)/1000*100)
+plt.plot(Energy,np.array(hhits)/np.array(ehits)*100,label='HCal / ECal')
+plt.scatter(Energy,np.array(hhits)/np.array(ehits)*100)
+
+plt.legend()
+plt.xlabel('Simulation Truth Gamma Energy (GeV)')
+plt.ylabel('Simulation Hits at ZDC (%)')
+plt.savefig('results/Hits.pdf')
+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("results/plots.pdf") as pdf:
+    pdf.savefig(fig1)
+    pdf.savefig(fig2)
+    pdf.savefig(fig3)
+    pdf.savefig(fig4)
+    pdf.savefig(fig5)

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;
+}