Insert tau

.gitlab-ci.yml

@@ -143,6 +143,7 @@ include:
   - local: 'benchmarks/timing/config.yml'
   - local: 'benchmarks/b0_tracker/config.yml'
   - local: 'benchmarks/insert_muon/config.yml'
+  - local: 'benchmarks/insert_tau/config.yml'
   - local: 'benchmarks/zdc_sigma/config.yml'
   - local: 'benchmarks/zdc_lambda/config.yml'
   - local: 'benchmarks/insert_neutron/config.yml'
@@ -169,6 +170,7 @@ deploy_results:
     - "collect_results:zdc"
     - "collect_results:zdc_lyso"
     - "collect_results:insert_muon"
+    - "collect_results:insert_tau"
     - "collect_results:zdc_photon"
     - "collect_results:zdc_pi0"
     - "collect_results:femc_electron"

Snakefile

@@ -8,13 +8,14 @@ include: "benchmarks/material_scan/Snakefile"
 include: "benchmarks/tracking_performances/Snakefile"
 include: "benchmarks/tracking_performances_dis/Snakefile"
 include: "benchmarks/lfhcal/Snakefile"
+include: "benchmarks/zdc_lyso/Snakefile"
 include: "benchmarks/insert_muon/Snakefile"
 include: "benchmarks/zdc_lambda/Snakefile"
-include: "benchmarks/zdc_lyso/Snakefile"
 include: "benchmarks/zdc_photon/Snakefile"
 include: "benchmarks/zdc_pi0/Snakefile"
 include: "benchmarks/zdc_sigma/Snakefile"
 include: "benchmarks/insert_neutron/Snakefile"
+include: "benchmarks/insert_tau/Snakefile"
 include: "benchmarks/femc_electron/Snakefile"
 include: "benchmarks/femc_photon/Snakefile"
 include: "benchmarks/femc_pi0/Snakefile"

benchmarks/insert_muon/analysis/muon_plots.py

@@ -1,83 +1,83 @@
-import numpy as np, pandas as pd, matplotlib.pyplot as plt, matplotlib as mpl, awkward as ak, sys
-import mplhep as hep
-hep.style.use("CMS")
-
-plt.rcParams['figure.facecolor']='white'
-plt.rcParams['savefig.facecolor']='white'
-plt.rcParams['savefig.bbox']='tight'
-
-plt.rcParams["figure.figsize"] = (7, 7)
-
-config=sys.argv[1].split("/")[1]  #results/{config}/{benchmark_name}
-outdir=sys.argv[1]+"/"
-try:
-    import os
-    os.mkdir(outdir[:-1])
-except:
-    pass
-
-def Landau(x, normalization,location,stdev):
-    #print(type(x))
-    u=(x-location)*3.591/stdev/2.355
-    renormalization = 1.64872*normalization
-    return renormalization * np.exp(-u/2 - np.exp(-u)/2)
-
-import uproot as ur
-arrays_sim={}
-momenta=50,
-for p in momenta:
-    arrays_sim[p] = ur.concatenate({
-        f'sim_output/insert_muon/{config}_sim_mu-_{p}GeV_{index}.edm4hep.root': 'events'
-        for index in range(5)
-    })
-
-for array in arrays_sim.values():
-    tilt=-0.025
-    px=array['MCParticles.momentum.x'][:,2]
-    py=array['MCParticles.momentum.y'][:,2]
-    pz=array['MCParticles.momentum.z'][:,2]
-    p=np.sqrt(px**2+py**2+pz**2)
-
-    pxp=px*np.cos(tilt)-pz*np.sin(tilt)
-    pyp=py
-    pzp=pz*np.cos(tilt)+px*np.sin(tilt)
-
-    array['eta_truth']=1/2*np.log((p+pzp)/(p-pzp))
-    array['phi_truth']=np.arctan2(pyp,pxp)
-
-for p in 50,:
-    E=arrays_sim[p]["HcalEndcapPInsertHits.energy"]
-    y, x,_=plt.hist(1e3*ak.flatten(E),bins=100, range=(0, 1.2), histtype='step')
-    bc=(x[1:]+x[:-1])/2
-    from scipy.optimize import curve_fit
-    slc=abs(bc-.48)<.15
-    fnc=Landau
-    p0=[100, .5, .05]
-    #print(list(y), list(x))
-    coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,
-                                 sigma=list(np.sqrt(y[slc])+(y[slc]==0)), maxfev=10000)
-    print(coeff)
-    xx=np.linspace(0,.7, 100)
-    MIP=coeff[1]/1000
-    plt.plot(xx, fnc(xx,*coeff), label=f'Landau fit:\nMIP={coeff[1]*1e3:.0f}$\\pm${1e3*np.sqrt(var_matrix[1][1]):.0f} keV')
-    plt.xlabel("hit energy [MeV]")
-    plt.ylabel("hits")
-    plt.title(f"$E_{{\\mu^-}}=${p} GeV")
-    plt.legend(fontsize=20)
-    plt.savefig(outdir+"/MIP.pdf")
-
-    plt.figure(figsize=(10,7))
-    array=arrays_sim[p]
-    bins=30; r=((-np.pi, np.pi),(2.8, 4.2))
-    selection=np.sum(array["HcalEndcapPInsertHits.energy"]>0.5*MIP,axis=-1)>0
-    h1, xedges, yedges = np.histogram2d(list(array[selection]['phi_truth']),list(array[selection]['eta_truth']), bins=bins, range=r)
-    h2, xedges, yedges = np.histogram2d(list(array['phi_truth']),list(array['eta_truth']), bins=bins, range=r)
-
-    h = h1 / h2
-    pc=plt.pcolor(xedges, yedges, h.T,linewidth=0)
-    plt.xlabel("$\\phi^*$ [rad]")
-    plt.ylabel("$\\eta^*$")
-    cb = plt.colorbar(pc)
-    cb.set_label("acceptance")
-    plt.title(f"$E_{{\\mu^-}}=${p} GeV")
-    plt.savefig(outdir+"/acceptance.pdf")
+import numpy as np, pandas as pd, matplotlib.pyplot as plt, matplotlib as mpl, awkward as ak, sys
+import mplhep as hep
+hep.style.use("CMS")
+
+plt.rcParams['figure.facecolor']='white'
+plt.rcParams['savefig.facecolor']='white'
+plt.rcParams['savefig.bbox']='tight'
+
+plt.rcParams["figure.figsize"] = (7, 7)
+
+config=sys.argv[1].split("/")[1]  #results/{config}/{benchmark_name}
+outdir=sys.argv[1]+"/"
+try:
+    import os
+    os.mkdir(outdir[:-1])
+except:
+    pass
+
+def Landau(x, normalization,location,stdev):
+    #print(type(x))
+    u=(x-location)*3.591/stdev/2.355
+    renormalization = 1.64872*normalization
+    return renormalization * np.exp(-u/2 - np.exp(-u)/2)
+
+import uproot as ur
+arrays_sim={}
+momenta=50,
+for p in momenta:
+    arrays_sim[p] = ur.concatenate({
+        f'sim_output/insert_muon/{config}_sim_mu-_{p}GeV_{index}.edm4hep.root': 'events'
+        for index in range(5)
+    })
+
+for array in arrays_sim.values():
+    tilt=-0.025
+    px=array['MCParticles.momentum.x'][:,2]
+    py=array['MCParticles.momentum.y'][:,2]
+    pz=array['MCParticles.momentum.z'][:,2]
+    p=np.sqrt(px**2+py**2+pz**2)
+    
+    pxp=px*np.cos(tilt)-pz*np.sin(tilt)
+    pyp=py
+    pzp=pz*np.cos(tilt)+px*np.sin(tilt)
+    
+    array['eta_truth']=1/2*np.log((p+pzp)/(p-pzp))
+    array['phi_truth']=np.arctan2(pyp,pxp)
+    
+for p in 50,:
+    E=arrays_sim[p]["HcalEndcapPInsertHits.energy"]
+    y, x,_=plt.hist(1e3*ak.flatten(E),bins=100, range=(0, 1.2), histtype='step')
+    bc=(x[1:]+x[:-1])/2
+    from scipy.optimize import curve_fit
+    slc=abs(bc-.48)<.15
+    fnc=Landau
+    p0=[100, .5, .05]
+    #print(list(y), list(x))
+    coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,
+                                 sigma=list(np.sqrt(y[slc])+(y[slc]==0)), maxfev=10000)
+    print(coeff)
+    xx=np.linspace(0,.7, 100)
+    MIP=coeff[1]/1000
+    plt.plot(xx, fnc(xx,*coeff), label=f'Landau fit:\nMIP={coeff[1]*1e3:.0f}$\\pm${1e3*np.sqrt(var_matrix[1][1]):.0f} keV')
+    plt.xlabel("hit energy [MeV]")
+    plt.ylabel("hits")
+    plt.title(f"$E_{{\\mu^-}}=${p} GeV")
+    plt.legend(fontsize=20)
+    plt.savefig(outdir+"/MIP.pdf")
+
+    plt.figure(figsize=(10,7))
+    array=arrays_sim[p]
+    bins=30; r=((-np.pi, np.pi),(2.8, 4.2))
+    selection=np.sum(array["HcalEndcapPInsertHits.energy"]>0.5*MIP,axis=-1)>0
+    h1, xedges, yedges = np.histogram2d(list(array[selection]['phi_truth']),list(array[selection]['eta_truth']), bins=bins, range=r)
+    h2, xedges, yedges = np.histogram2d(list(array['phi_truth']),list(array['eta_truth']), bins=bins, range=r)
+
+    h = h1 / h2
+    pc=plt.pcolor(xedges, yedges, h.T,linewidth=0)
+    plt.xlabel("$\\phi^*$ [rad]")
+    plt.ylabel("$\\eta^*$")
+    cb = plt.colorbar(pc)
+    cb.set_label("acceptance")
+    plt.title(f"$E_{{\\mu^-}}=${p} GeV")
+    plt.savefig(outdir+"/acceptance.pdf")

benchmarks/insert_tau/Snakefile

@@ -0,0 +1,61 @@
+rule insert_tau_generate:
+	input:
+                script="benchmarks/insert_tau/analysis/gen_particles.cxx",
+	params:
+		NEVENTS_GEN=1000,
+		th_max=7.0,
+		th_min=1.7
+	output:
+		GEN_FILE="sim_output/insert_tau/tau-_{P}GeV.hepmc",
+	shell:
+		"""
+root -l -b -q '{input.script}({params.NEVENTS_GEN},"{output.GEN_FILE}", "tau-", {params.th_min}, {params.th_max}, 0., 360., {wildcards.P})'
+"""
+
+rule insert_tau_simulate:
+	input:
+		GEN_FILE="sim_output/insert_tau/tau-_{P}GeV.hepmc"
+		warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root"
+	params:
+		PHYSICS_LIST="FTFP_BERT"
+	output:
+		SIM_FILE="sim_output/insert_tau/{DETECTOR_CONFIG}_sim_tau-_{P}GeV_{INDEX}.edm4hep.root"
+	shell:
+		"""
+NEVENTS_SIM=1000
+# Running simulation
+npsim \
+   --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \
+   --numberOfEvents $NEVENTS_SIM \
+   --skipNEvents $(( $NEVENTS_SIM * {wildcards.INDEX} )) \
+   --physicsList {params.PHYSICS_LIST} \
+   --inputFiles {input.GEN_FILE} \
+   --outputFile {output.SIM_FILE}
+"""
+
+rule insert_tau_recon:
+        input:
+                SIM_FILE="sim_output/insert_tau/{DETECTOR_CONFIG}_sim_tau-_{P}GeV_{INDEX}.edm4hep.root"
+        output:
+                REC_FILE="sim_output/insert_tau/{DETECTOR_CONFIG}_rec_tau-_{P}GeV_{INDEX}.edm4hep.root"
+        shell:
+                """
+NEVENTS_REC=1000
+eicrecon {input.SIM_FILE} -Ppodio:output_file={output.REC_FILE} -Pdd4hep:xml_files=$DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml -Ppodio:output_collections=MCParticles,HcalEndcapPInsertRecHits,HcalEndcapPInsertClusters,HcalEndcapPInsertSubcellHits,EcalEndcapPInsertRecHits,EcalEndcapPInsertClusters,EcalEndcapPClusters,LFHCALClusters  -Pjana:nevents=$NEVENTS_REC
+"""
+
+rule insert_tau_analysis:
+	input:
+                expand("sim_output/insert_tau/{DETECTOR_CONFIG}_rec_tau-_{P}GeV.edm4hep_{INDEX}.root",
+		    P=[20,30, 40,50,60, 80, 100],
+		    DETECTOR_CONFIG=["{DETECTOR_CONFIG}"],
+		    INDEX=range(5),
+                 ),
+                script="benchmarks/insert_tau/analysis/tau_plots.py",
+	output:
+		results_dir=directory("results/{DETECTOR_CONFIG}/insert_tau"),
+	shell:
+		"""
+mkdir -p {output.results_dir}
+python {input.script} {output.results_dir}
+"""

benchmarks/insert_tau/Snakefile~

@@ -0,0 +1,57 @@
+rule insert_tau_generate:
+	input:
+                script="benchmarks/insert_tau/analysis/gen_particles.cxx",
+	params:
+		NEVENTS_GEN=1000,
+		th_max=7.0,
+		th_min=1.7
+	output:
+		GEN_FILE="sim_output/insert_tau/tau-_{P}GeV.hepmc"
+	shell:
+		"""
+root -l -b -q '{input.script}({params.NEVENTS_GEN},"{output.GEN_FILE}", "tau-", {params.th_min}, {params.th_max}, 0., 360., {wildcards.P})'
+"""
+
+rule insert_tau_simulate:
+	input:
+		GEN_FILE="sim_output/insert_tau/tau-_{P}GeV.hepmc"
+	params:
+		PHYSICS_LIST="FTFP_BERT"
+	output:
+		SIM_FILE="sim_output/insert_tau/{DETECTOR_CONFIG}_sim_tau-_{P}GeV.edm4hep.root"
+	shell:
+		"""
+NEVENTS_SIM=1000
+# Running simulation
+npsim \
+   --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \
+   --numberOfEvents $NEVENTS_SIM \
+   --physicsList {params.PHYSICS_LIST} \
+   --inputFiles {input.GEN_FILE} \
+   --outputFile {output.SIM_FILE}
+"""
+
+rule insert_tau_recon:
+        input:
+                SIM_FILE="sim_output/insert_tau/{DETECTOR_CONFIG}_sim_tau-_{P}GeV.edm4hep.root"
+        output:
+                REC_FILE="sim_output/insert_tau/{DETECTOR_CONFIG}_rec_tau-_{P}GeV.edm4hep.root"
+        shell:
+                """
+NEVENTS_REC=1000
+eicrecon {input.SIM_FILE} -Ppodio:output_file={output.REC_FILE} -Pdd4hep:xml_files=$DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml -Ppodio:output_collections=MCParticles,HcalEndcapPInsertRecHits,HcalEndcapPInsertClusters,HcalEndcapPInsertSubcellHits,EcalEndcapPInsertRecHits,EcalEndcapPInsertClusters  -Pjana:nevents=$NEVENTS_REC
+"""
+
+rule insert_tau_analysis:
+	input:
+                expand("sim_output/insert_tau/{DETECTOR_CONFIG}_rec_tau-_{P}GeV.edm4hep.root",
+		    P=[20,50, 100, 200],
+		    DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]),
+                script="benchmarks/insert_tau/analysis/tau_plots.py",
+	output:
+		results_dir=directory("results/{DETECTOR_CONFIG}/insert_tau"),
+	shell:
+		"""
+mkdir -p {output.results_dir}
+python {input.script} {output.results_dir}
+"""