From f04a3693fe86c04b4314ecc25df4778431da8ecc Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Fri, 27 Sep 2024 09:01:38 -0400 Subject: [PATCH 1/3] minor typo in benchmarks/backgrounds/config.yml (#73) --- benchmarks/backgrounds/config.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/backgrounds/config.yml b/benchmarks/backgrounds/config.yml index 1f9775d8..c4d084d7 100644 --- a/benchmarks/backgrounds/config.yml +++ b/benchmarks/backgrounds/config.yml @@ -5,7 +5,7 @@ sim:backgrounds: - mkdir -p $LOCAL_DATA_PATH/input - ln -s $LOCAL_DATA_PATH/input input - | - snakemake -cache --cores 2 \ + snakemake --cache --cores 2 \ sim_output/$DETECTOR_CONFIG/backgrounds/EPIC/EVGEN/BACKGROUNDS/BEAMGAS/electron/GETaLM1.0.0-1.0/10GeV/GETaLM1.0.0-1.0_ElectronBeamGas_10GeV_foam_emin10keV_run001.edm4hep.root \ sim_output/$DETECTOR_CONFIG/backgrounds/EPIC/EVGEN/DIS/NC/10x100/minQ2=1/pythia8NCDIS_10x100_minQ2=1_beamEffects_xAngle=-0.025_hiDiv_1.edm4hep.root \ sim_output/$DETECTOR_CONFIG/backgrounds/EPIC/EVGEN/BACKGROUNDS/BEAMGAS/proton/pythia8.306-1.0/100GeV/pythia8.306-1.0_ProtonBeamGas_100GeV_run001.edm4hep.root From c72c399b3bba78d55114f6cafd18a7ebe79cd47b Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Wed, 2 Oct 2024 08:39:47 -0400 Subject: [PATCH 2/3] parallelize zdc and insert simulations (#74) --- benchmarks/insert_muon/Snakefile | 20 ++- benchmarks/insert_muon/analysis/muon_plots.py | 169 +++++++++--------- benchmarks/insert_muon/config.yml | 3 +- benchmarks/insert_neutron/Snakefile | 20 ++- .../insert_neutron/analysis/neutron_plots.py | 6 +- benchmarks/insert_neutron/config.yml | 3 +- benchmarks/zdc_lambda/Snakefile | 30 ++-- .../zdc_lambda/analysis/gen_lambda_decay.cxx | 15 +- .../zdc_lambda/analysis/lambda_plots.py | 11 +- benchmarks/zdc_lambda/config.yml | 3 +- benchmarks/zdc_photon/Snakefile | 34 ++-- .../zdc_photon/analysis/zdc_photon_plots.py | 13 +- benchmarks/zdc_photon/config.yml | 3 +- benchmarks/zdc_pi0/Snakefile | 20 ++- benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py | 11 +- benchmarks/zdc_pi0/config.yml | 3 +- benchmarks/zdc_sigma/Snakefile | 28 ++- .../zdc_sigma/analysis/gen_sigma_decay.cxx | 15 +- benchmarks/zdc_sigma/analysis/sigma_plots.py | 11 +- benchmarks/zdc_sigma/config.yml | 3 +- 20 files changed, 199 insertions(+), 222 deletions(-) diff --git a/benchmarks/insert_muon/Snakefile b/benchmarks/insert_muon/Snakefile index b36026b0..a9a32168 100644 --- a/benchmarks/insert_muon/Snakefile +++ b/benchmarks/insert_muon/Snakefile @@ -14,17 +14,19 @@ root -l -b -q '{input.script}({params.NEVENTS_GEN},"{output.GEN_FILE}", "mu-", { rule insert_muon_simulate: input: - GEN_FILE="sim_output/insert_muon/mu-_{P}GeV.hepmc" + GEN_FILE="sim_output/insert_muon/mu-_{P}GeV.hepmc", + warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root", params: PHYSICS_LIST="FTFP_BERT" output: - SIM_FILE="sim_output/insert_muon/{DETECTOR_CONFIG}_sim_mu-_{P}GeV.edm4hep.root" + SIM_FILE="sim_output/insert_muon/{DETECTOR_CONFIG}_sim_mu-_{P}GeV_{INDEX}.edm4hep.root" shell: """ -NEVENTS_SIM=5000 +NEVENTS_SIM=1000 # Running simulation npsim \ --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \ + --skipNEvents $(( $NEVENTS_SIM * {wildcards.INDEX} )) \ --numberOfEvents $NEVENTS_SIM \ --physicsList {params.PHYSICS_LIST} \ --inputFiles {input.GEN_FILE} \ @@ -33,20 +35,22 @@ npsim \ rule insert_muon_recon: input: - SIM_FILE="sim_output/insert_muon/{DETECTOR_CONFIG}_sim_mu-_{P}GeV.edm4hep.root" + SIM_FILE="sim_output/insert_muon/{DETECTOR_CONFIG}_sim_mu-_{P}GeV_{INDEX}.edm4hep.root" output: - REC_FILE="sim_output/insert_muon/{DETECTOR_CONFIG}_rec_mu-_{P}GeV.edm4hep.root" + REC_FILE="sim_output/insert_muon/{DETECTOR_CONFIG}_rec_mu-_{P}GeV_{INDEX}.edm4hep.root" shell: """ -NEVENTS_REC=5000 +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_muon_analysis: input: - expand("sim_output/insert_muon/{DETECTOR_CONFIG}_sim_mu-_{P}GeV.edm4hep.root", + expand("sim_output/insert_muon/{DETECTOR_CONFIG}_sim_mu-_{P}GeV_{INDEX}.edm4hep.root", P=[50], - DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]), + DETECTOR_CONFIG=["{DETECTOR_CONFIG}"], + INDEX=range(5), + ), script="benchmarks/insert_muon/analysis/muon_plots.py", output: results_dir=directory("results/{DETECTOR_CONFIG}/insert_muon"), diff --git a/benchmarks/insert_muon/analysis/muon_plots.py b/benchmarks/insert_muon/analysis/muon_plots.py index 4c45aac3..5c81b52d 100644 --- a/benchmarks/insert_muon/analysis/muon_plots.py +++ b/benchmarks/insert_muon/analysis/muon_plots.py @@ -1,86 +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: - filename=f'sim_output/insert_muon/{config}_sim_mu-_{p}GeV.edm4hep.root' - print("opening file", filename) - events = ur.open(filename+':events') - arrays_sim[p] = events.arrays() - import gc - gc.collect() - print("read", filename) - -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))) - 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))) + 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") diff --git a/benchmarks/insert_muon/config.yml b/benchmarks/insert_muon/config.yml index 6f7681f1..252cd6e8 100644 --- a/benchmarks/insert_muon/config.yml +++ b/benchmarks/insert_muon/config.yml @@ -4,9 +4,8 @@ sim:insert_muon: parallel: matrix: - P: 50 - timeout: 1 hours script: - - snakemake --cores 1 sim_output/insert_muon/epic_craterlake_sim_mu-_${P}GeV.edm4hep.root + - snakemake --cores 5 sim_output/insert_muon/epic_craterlake_sim_mu-_${P}GeV_{0,1,2,3,4}.edm4hep.root retry: max: 2 when: diff --git a/benchmarks/insert_neutron/Snakefile b/benchmarks/insert_neutron/Snakefile index 55e3ff9e..136c56e2 100644 --- a/benchmarks/insert_neutron/Snakefile +++ b/benchmarks/insert_neutron/Snakefile @@ -15,17 +15,19 @@ root -l -b -q '{input.script}({params.NEVENTS_GEN},"{output.GEN_FILE}", "neutron rule insert_neutron_simulate: input: - GEN_FILE="sim_output/insert_neutron/neutron_{P}GeV.hepmc" + GEN_FILE="sim_output/insert_neutron/neutron_{P}GeV.hepmc", + warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root", params: PHYSICS_LIST="FTFP_BERT" output: - SIM_FILE="sim_output/insert_neutron/{DETECTOR_CONFIG}_sim_neutron_{P}GeV.edm4hep.root" + SIM_FILE="sim_output/insert_neutron/{DETECTOR_CONFIG}_sim_neutron_{P}GeV_{INDEX}.edm4hep.root" shell: """ -NEVENTS_SIM=1000 +NEVENTS_SIM=200 # Running simulation npsim \ --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \ + --skipNEvents $(( $NEVENTS_SIM * {wildcards.INDEX} )) \ --numberOfEvents $NEVENTS_SIM \ --physicsList {params.PHYSICS_LIST} \ --inputFiles {input.GEN_FILE} \ @@ -34,20 +36,22 @@ npsim \ rule insert_neutron_recon: input: - SIM_FILE="sim_output/insert_neutron/{DETECTOR_CONFIG}_sim_neutron_{P}GeV.edm4hep.root" + SIM_FILE="sim_output/insert_neutron/{DETECTOR_CONFIG}_sim_neutron_{P}GeV_{INDEX}.edm4hep.root" output: - REC_FILE="sim_output/insert_neutron/{DETECTOR_CONFIG}_rec_neutron_{P}GeV.edm4eic.root" + REC_FILE="sim_output/insert_neutron/{DETECTOR_CONFIG}_rec_neutron_{P}GeV_{INDEX}.edm4eic.root" shell: """ -NEVENTS_REC=1000 +NEVENTS_REC=200 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_neutron_analysis: input: - expand("sim_output/insert_neutron/{DETECTOR_CONFIG}_rec_neutron_{P}GeV.edm4eic.root", + expand("sim_output/insert_neutron/{DETECTOR_CONFIG}_rec_neutron_{P}GeV_{INDEX}.edm4eic.root", P=[20, 30, 40, 50, 60, 70, 80], - DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]), + DETECTOR_CONFIG=["{DETECTOR_CONFIG}"], + INDEX=range(5), + ), script="benchmarks/insert_neutron/analysis/neutron_plots.py", output: results_dir=directory("results/{DETECTOR_CONFIG}/insert_neutron"), diff --git a/benchmarks/insert_neutron/analysis/neutron_plots.py b/benchmarks/insert_neutron/analysis/neutron_plots.py index 9314a974..b9d2002a 100644 --- a/benchmarks/insert_neutron/analysis/neutron_plots.py +++ b/benchmarks/insert_neutron/analysis/neutron_plots.py @@ -18,8 +18,10 @@ #read files arrays_sim={} for p in 20,30,40,50,60,70,80: - arrays_sim[p] = ur.open(f'sim_output/insert_neutron/{config}_rec_neutron_{p}GeV.edm4eic.root:events')\ - .arrays() + arrays_sim[p] = ur.concatenate({ + f'sim_output/insert_neutron/{config}_rec_neutron_{p}GeV_{index}.edm4eic.root': 'events' + for index in range(5) + }) def gauss(x, A,mu, sigma): return A * np.exp(-(x-mu)**2/(2*sigma**2)) diff --git a/benchmarks/insert_neutron/config.yml b/benchmarks/insert_neutron/config.yml index aac0a732..4c723e68 100644 --- a/benchmarks/insert_neutron/config.yml +++ b/benchmarks/insert_neutron/config.yml @@ -10,9 +10,8 @@ sim:insert_neutron: - P: 60 - P: 70 - P: 80 - timeout: 1 hours script: - - snakemake --cores 1 sim_output/insert_neutron/epic_craterlake_rec_neutron_${P}GeV.edm4eic.root + - snakemake --cores 5 sim_output/insert_neutron/epic_craterlake_rec_neutron_${P}GeV_{0,1,2,3,4}.edm4eic.root retry: max: 2 when: diff --git a/benchmarks/zdc_lambda/Snakefile b/benchmarks/zdc_lambda/Snakefile index 21563330..9ffd1701 100644 --- a/benchmarks/zdc_lambda/Snakefile +++ b/benchmarks/zdc_lambda/Snakefile @@ -2,7 +2,7 @@ rule zdc_lambda_generate: input: script="benchmarks/zdc_lambda/analysis/gen_lambda_decay.cxx", params: - NEVENTS_GEN=100000, + NEVENTS_GEN=1000, output: GEN_FILE="sim_output/zdc_lambda/lambda_decay_{P}GeV.hepmc" shell: @@ -12,21 +12,19 @@ root -l -b -q '{input.script}({params.NEVENTS_GEN},0,"{output.GEN_FILE}",{wildca rule zdc_lambda_simulate: input: - GEN_FILE="sim_output/zdc_lambda/lambda_decay_{P}GeV.hepmc" + GEN_FILE="sim_output/zdc_lambda/lambda_decay_{P}GeV.hepmc", + warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root", params: PHYSICS_LIST="FTFP_BERT" output: - SIM_FILE="sim_output/zdc_lambda/{DETECTOR_CONFIG}_sim_lambda_dec_{P}GeV.edm4hep.root" + SIM_FILE="sim_output/zdc_lambda/{DETECTOR_CONFIG}_sim_lambda_dec_{P}GeV_{INDEX}.edm4hep.root" shell: """ -if [[ {wildcards.P} -gt 220 ]]; then - NEVENTS_SIM=500 -else - NEVENTS_SIM=1000 -fi +NEVENTS_SIM=200 # Running simulation npsim \ --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \ + --skipNEvents $(( $NEVENTS_SIM * {wildcards.INDEX} )) \ --numberOfEvents $NEVENTS_SIM \ --physicsList {params.PHYSICS_LIST} \ --inputFiles {input.GEN_FILE} \ @@ -35,24 +33,22 @@ npsim \ rule zdc_lambda_recon: input: - SIM_FILE="sim_output/zdc_lambda/{DETECTOR_CONFIG}_sim_lambda_dec_{P}GeV.edm4hep.root" + SIM_FILE="sim_output/zdc_lambda/{DETECTOR_CONFIG}_sim_lambda_dec_{P}GeV_{INDEX}.edm4hep.root" output: - REC_FILE="sim_output/zdc_lambda/{DETECTOR_CONFIG}_rec_lambda_dec_{P}GeV.edm4eic.root" + REC_FILE="sim_output/zdc_lambda/{DETECTOR_CONFIG}_rec_lambda_dec_{P}GeV_{INDEX}.edm4eic.root" shell: """ -if [[ {wildcards.P} -gt 220 ]]; then - NEVENTS_REC=500 -else - NEVENTS_REC=1000 -fi +NEVENTS_REC=200 eicrecon {input.SIM_FILE} -Ppodio:output_file={output.REC_FILE} -Pdd4hep:xml_files=$DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml -Ppodio:output_collections=MCParticles,HcalFarForwardZDCClusters,HcalFarForwardZDCRecHits,HcalFarForwardZDCSubcellHits -Pjana:nevents=$NEVENTS_REC """ rule zdc_lambda_analysis: input: - expand("sim_output/zdc_lambda/{DETECTOR_CONFIG}_rec_lambda_dec_{P}GeV.edm4eic.root", + expand("sim_output/zdc_lambda/{DETECTOR_CONFIG}_rec_lambda_dec_{P}GeV_{INDEX}.edm4eic.root", P=[100, 125, 150,175, 200, 225, 250, 275], - DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]), + DETECTOR_CONFIG=["{DETECTOR_CONFIG}"], + INDEX=range(5), + ), script="benchmarks/zdc_lambda/analysis/lambda_plots.py", output: results_dir=directory("results/{DETECTOR_CONFIG}/zdc_lambda"), diff --git a/benchmarks/zdc_lambda/analysis/gen_lambda_decay.cxx b/benchmarks/zdc_lambda/analysis/gen_lambda_decay.cxx index 567eda5b..1b399dda 100644 --- a/benchmarks/zdc_lambda/analysis/gen_lambda_decay.cxx +++ b/benchmarks/zdc_lambda/analysis/gen_lambda_decay.cxx @@ -43,7 +43,6 @@ void gen_lambda_decay(int n_events = 100000, UInt_t seed = 0, char* out_fname = //const double p_max = 275.; // in GeV/c WriterAscii hepmc_output(out_fname); - int events_parsed = 0; GenEvent evt(Units::GEV, Units::MM); // Random number generator @@ -71,10 +70,11 @@ void gen_lambda_decay(int n_events = 100000, UInt_t seed = 0, char* out_fname = double photon_mass = std::get<0>(photon_info); int photon_pdgID = std::get<1>(photon_info); - for (events_parsed = 0; events_parsed < n_events; events_parsed++) { + int accepted_events = 0; + while (accepted_events < n_events) { //Set the event number - evt.set_event_number(events_parsed); + evt.set_event_number(accepted_events); // FourVector(px,py,pz,e,pdgid,status) // type 4 is beam @@ -218,7 +218,7 @@ void gen_lambda_decay(int n_events = 100000, UInt_t seed = 0, char* out_fname = evt.add_vertex(v_pi0_decay); - if (events_parsed == 0) { + if (accepted_events == 0) { std::cout << "First event: " << std::endl; Print::listing(evt); } @@ -227,13 +227,14 @@ void gen_lambda_decay(int n_events = 100000, UInt_t seed = 0, char* out_fname = TVector3 extrap_gamma1=lambda_decay_position+gamma1_lab.Vect()*((zdc_z-pbeam_dir.Dot(lambda_decay_position))/(pbeam_dir.Dot(gamma1_lab.Vect()))); TVector3 extrap_gamma2=lambda_decay_position+gamma2_lab.Vect()*((zdc_z-pbeam_dir.Dot(lambda_decay_position))/(pbeam_dir.Dot(gamma2_lab.Vect()))); if (extrap_neutron.Angle(pbeam_dir)<0.004 && extrap_gamma1.Angle(pbeam_dir)<0.004 && extrap_gamma2.Angle(pbeam_dir)<0.004 && lambda_decay_position.Dot(pbeam_dir)(photon_info); int photon_pdgID = std::get<1>(photon_info); - for (events_parsed = 0; events_parsed < n_events; events_parsed++) { + int accepted_events = 0; + while (accepted_events < n_events) { //Set the event number - evt.set_event_number(events_parsed); + evt.set_event_number(accepted_events); // FourVector(px,py,pz,e,pdgid,status) // type 4 is beam @@ -281,7 +280,7 @@ void gen_sigma_decay(int n_events = 100000, UInt_t seed = 0, char* out_fname = " evt.add_vertex(v_pi0_decay); - if (events_parsed == 0) { + if (accepted_events == 0) { std::cout << "First event: " << std::endl; Print::listing(evt); } @@ -294,12 +293,12 @@ void gen_sigma_decay(int n_events = 100000, UInt_t seed = 0, char* out_fname = " hepmc_output.write_event(evt); accepted_events ++; } - if (events_parsed % 1000 == 0) { - std::cout << "Event: " << events_parsed << " ("< Date: Wed, 2 Oct 2024 11:46:27 -0400 Subject: [PATCH 3/3] Femc electron. photon, and pi0 benchmarks (#75) --- .gitlab-ci.yml | 7 +- Snakefile | 5 +- benchmarks/femc_electron/Snakefile | 58 ++++ .../analysis/femc_electron_plots.py | 222 +++++++++++++++ .../femc_electron/analysis/gen_particles.cxx | 127 +++++++++ benchmarks/femc_electron/config.yml | 36 +++ benchmarks/femc_photon/Snakefile | 58 ++++ .../femc_photon/analysis/femc_photon_plots.py | 220 +++++++++++++++ .../analysis/femc_photon_plots.py~ | 171 ++++++++++++ .../femc_photon/analysis/gen_particles.cxx | 127 +++++++++ benchmarks/femc_photon/config.yml | 36 +++ benchmarks/femc_pi0/Snakefile | 58 ++++ .../femc_pi0/analysis/femc_pi0_plots.py | 257 ++++++++++++++++++ .../femc_pi0/analysis/gen_particles.cxx | 127 +++++++++ benchmarks/femc_pi0/config.yml | 36 +++ 15 files changed, 1542 insertions(+), 3 deletions(-) create mode 100644 benchmarks/femc_electron/Snakefile create mode 100644 benchmarks/femc_electron/analysis/femc_electron_plots.py create mode 100644 benchmarks/femc_electron/analysis/gen_particles.cxx create mode 100644 benchmarks/femc_electron/config.yml create mode 100644 benchmarks/femc_photon/Snakefile create mode 100644 benchmarks/femc_photon/analysis/femc_photon_plots.py create mode 100644 benchmarks/femc_photon/analysis/femc_photon_plots.py~ create mode 100644 benchmarks/femc_photon/analysis/gen_particles.cxx create mode 100644 benchmarks/femc_photon/config.yml create mode 100644 benchmarks/femc_pi0/Snakefile create mode 100644 benchmarks/femc_pi0/analysis/femc_pi0_plots.py create mode 100644 benchmarks/femc_pi0/analysis/gen_particles.cxx create mode 100644 benchmarks/femc_pi0/config.yml diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index a1dc5d8c..a56c79d2 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -144,7 +144,9 @@ include: - local: 'benchmarks/zdc_sigma/config.yml' - local: 'benchmarks/zdc_lambda/config.yml' - local: 'benchmarks/insert_neutron/config.yml' - + - local: 'benchmarks/femc_electron/config.yml' + - local: 'benchmarks/femc_photon/config.yml' + - local: 'benchmarks/femc_pi0/config.yml' deploy_results: allow_failure: true stage: deploy @@ -167,6 +169,9 @@ deploy_results: - "collect_results:insert_muon" - "collect_results:zdc_photon" - "collect_results:zdc_pi0" + - "collect_results:femc_electron" + - "collect_results:femc_photon" + - "collect_results:femc_pi0" script: - echo "deploy results!" - find results -print | sort | tee summary.txt diff --git a/Snakefile b/Snakefile index bce37f8d..54323284 100644 --- a/Snakefile +++ b/Snakefile @@ -7,7 +7,6 @@ include: "benchmarks/ecal_gaps/Snakefile" include: "benchmarks/material_scan/Snakefile" include: "benchmarks/tracking_performances/Snakefile" include: "benchmarks/tracking_performances_dis/Snakefile" -include: "benchmarks/zdc_lyso/Snakefile" include: "benchmarks/insert_muon/Snakefile" include: "benchmarks/zdc_lambda/Snakefile" include: "benchmarks/zdc_lyso/Snakefile" @@ -15,7 +14,9 @@ include: "benchmarks/zdc_photon/Snakefile" include: "benchmarks/zdc_pi0/Snakefile" include: "benchmarks/zdc_sigma/Snakefile" include: "benchmarks/insert_neutron/Snakefile" - +include: "benchmarks/femc_electron/Snakefile" +include: "benchmarks/femc_photon/Snakefile" +include: "benchmarks/femc_pi0/Snakefile" use_s3 = config["remote_provider"].lower() == "s3" use_xrootd = config["remote_provider"].lower() == "xrootd" diff --git a/benchmarks/femc_electron/Snakefile b/benchmarks/femc_electron/Snakefile new file mode 100644 index 00000000..0e8dd084 --- /dev/null +++ b/benchmarks/femc_electron/Snakefile @@ -0,0 +1,58 @@ +rule femc_electron_generate: + input: + script="benchmarks/femc_electron/analysis/gen_particles.cxx", + params: + NEVENTS_GEN=1000, + th_max=28, + th_min=2.0 + output: + GEN_FILE="sim_output/femc_electron/e-_{P}GeV.hepmc" + shell: + """ +mkdir -p sim_output/femc_electron +root -l -b -q '{input.script}({params.NEVENTS_GEN},"{output.GEN_FILE}", "e-", {params.th_min}, {params.th_max}, 0., 360., {wildcards.P})' +""" + +rule femc_electron_simulate: + input: + GEN_FILE="sim_output/femc_electron/e-_{P}GeV.hepmc" + params: + PHYSICS_LIST="FTFP_BERT" + output: + SIM_FILE="sim_output/femc_electron/{DETECTOR_CONFIG}_sim_e-_{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 femc_electron_recon: + input: + SIM_FILE="sim_output/femc_electron/{DETECTOR_CONFIG}_sim_e-_{P}GeV.edm4hep.root" + output: + REC_FILE="sim_output/femc_electron/{DETECTOR_CONFIG}_rec_e-_{P}GeV.edm4eic.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,EcalEndcapPInsertRecHits,EcalEndcapPClusters -Pjana:nevents=$NEVENTS_REC +""" + +rule femc_electron_analysis: + input: + expand("sim_output/femc_electron/{DETECTOR_CONFIG}_rec_e-_{P}GeV.edm4eic.root", + P=[10, 20, 30, 40, 50, 60, 70, 80], + DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]), + script="benchmarks/femc_electron/analysis/femc_electron_plots.py", + output: + results_dir=directory("results/{DETECTOR_CONFIG}/femc_electron"), + shell: + """ +mkdir -p {output.results_dir} +python {input.script} {output.results_dir} +""" diff --git a/benchmarks/femc_electron/analysis/femc_electron_plots.py b/benchmarks/femc_electron/analysis/femc_electron_plots.py new file mode 100644 index 00000000..cc34633d --- /dev/null +++ b/benchmarks/femc_electron/analysis/femc_electron_plots.py @@ -0,0 +1,222 @@ +import numpy as np, pandas as pd, matplotlib.pyplot as plt, matplotlib as mpl, awkward as ak, sys, uproot as ur +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 + +import uproot as ur +arrays_sim={p:ur.open(f'sim_output/femc_electron/{config}_rec_e-_{p}GeV.edm4eic.root:events').arrays() for p in (10, 20, 30, 40, 50, 60,70,80)} + +for p in arrays_sim: + array=arrays_sim[p] + tilt=-.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['nclust_endcap']=[len(array['EcalEndcapPClusters.energy'][i]) for i in range(len(array))] + +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) + +#number of clusters +plt.figure() +for eta_min, eta_max, field in (1.5, 2.8, 'nclust_endcap'),: + for p in arrays_sim: + array=arrays_sim[p] + plt.hist(array[field][(array['eta_truth']>eta_min)&(array['eta_truth']eta_min)]['EcalEndcapPClusters.energy']), bins=bins) + bcs=(x[1:]+x[:-1])/2 + + fnc=gauss + slc=abs(bcs-p)<3 + sigma=np.sqrt(y[slc])+0.5*(y[slc]==0) + p0=(100, p, 3) + + try: + coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma)) + if abs(coeff[1])>100 or np.sqrt(var_matrix[1][1])>100: + continue + pvals.append(p) + res.append(abs(coeff[2])/coeff[1]) + dres.append(np.sqrt(var_matrix[2][2])/coeff[1]) + scale.append(abs(coeff[1])/p) + dscale.append(np.sqrt(var_matrix[1][1])/p) + except: + pass + plt.sca(axs[0]) + plt.errorbar(pvals, 100*np.array(res), 100*np.array(dres), ls='', marker='o', label=f'${eta_min:.1f}<\\eta<{eta_max:.1f}$') + + plt.sca(axs[1]) + plt.errorbar(pvals, 100*np.array(scale), 100*np.array(dscale), ls='', marker='o', label=f'${eta_min:.1f}<\\eta<{eta_max:.1f}$') + +plt.sca(axs[0]) +plt.legend() +plt.ylim(0) +plt.ylabel("E resolution [%]") +plt.xlabel("E truth [GeV]") + +plt.sca(axs[1]) +plt.ylabel("energy scale [%]") +plt.xlabel("E truth [GeV]") +plt.axhline(100, color='0.5', alpha=0.5, ls='--') +plt.ylim(0, 110) + +plt.tight_layout() +plt.savefig(outdir+"/energy_res_eta_partitions.pdf") diff --git a/benchmarks/femc_electron/analysis/gen_particles.cxx b/benchmarks/femc_electron/analysis/gen_particles.cxx new file mode 100644 index 00000000..0877a852 --- /dev/null +++ b/benchmarks/femc_electron/analysis/gen_particles.cxx @@ -0,0 +1,127 @@ +#include "HepMC3/GenEvent.h" +#include "HepMC3/ReaderAscii.h" +#include "HepMC3/WriterAscii.h" +#include "HepMC3/Print.h" + +#include "TRandom3.h" +#include "TVector3.h" + +#include +#include +#include +#include +#include +#include +#include + +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 = 1000, + 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 + int useCrossingAngle=1 // 0= no rotation, 1 = -25 mrad + ) +{ + 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(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4); + GenParticlePtr p2 = std::make_shared( + 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.*useCrossingAngle; //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( + 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(); + 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; +} diff --git a/benchmarks/femc_electron/config.yml b/benchmarks/femc_electron/config.yml new file mode 100644 index 00000000..a600454d --- /dev/null +++ b/benchmarks/femc_electron/config.yml @@ -0,0 +1,36 @@ +sim:femc_electron: + stage: simulate + extends: .det_benchmark + parallel: + matrix: + - P: 10 + - P: 20 + - P: 30 + - P: 40 + - P: 50 + - P: 60 + - P: 70 + - P: 80 + timeout: 1 hours + script: + - snakemake --cores 1 sim_output/femc_electron/epic_craterlake_rec_e-_${P}GeV.edm4eic.root + retry: + max: 2 + when: + - runner_system_failure + +bench:femc_electron: + stage: benchmarks + extends: .det_benchmark + needs: ["sim:femc_electron"] + script: + - snakemake --cores 1 results/epic_craterlake/femc_electron + +collect_results:femc_electron: + stage: collect + needs: ["bench:femc_electron"] + script: + - ls -al + - mv results{,_save}/ # move results directory out of the way to preserve it + - snakemake --cores 1 --delete-all-output results/epic_craterlake/femc_electron + - mv results{_save,}/ diff --git a/benchmarks/femc_photon/Snakefile b/benchmarks/femc_photon/Snakefile new file mode 100644 index 00000000..40524365 --- /dev/null +++ b/benchmarks/femc_photon/Snakefile @@ -0,0 +1,58 @@ +rule femc_photon_generate: + input: + script="benchmarks/femc_photon/analysis/gen_particles.cxx", + params: + NEVENTS_GEN=1000, + th_max=28, + th_min=2.0 + output: + GEN_FILE="sim_output/femc_photon/photon_{P}GeV.hepmc" + shell: + """ +mkdir -p sim_output/femc_photon +root -l -b -q '{input.script}({params.NEVENTS_GEN},"{output.GEN_FILE}", "gamma", {params.th_min}, {params.th_max}, 0., 360., {wildcards.P})' +""" + +rule femc_photon_simulate: + input: + GEN_FILE="sim_output/femc_photon/photon_{P}GeV.hepmc" + params: + PHYSICS_LIST="FTFP_BERT" + output: + SIM_FILE="sim_output/femc_photon/{DETECTOR_CONFIG}_sim_photon_{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 femc_photon_recon: + input: + SIM_FILE="sim_output/femc_photon/{DETECTOR_CONFIG}_sim_photon_{P}GeV.edm4hep.root" + output: + REC_FILE="sim_output/femc_photon/{DETECTOR_CONFIG}_rec_photon_{P}GeV.edm4eic.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,EcalEndcapPInsertRecHits,EcalEndcapPClusters -Pjana:nevents=$NEVENTS_REC +""" + +rule femc_photon_analysis: + input: + expand("sim_output/femc_photon/{DETECTOR_CONFIG}_rec_photon_{P}GeV.edm4eic.root", + P=[10, 20, 30, 40, 50, 60, 70, 80], + DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]), + script="benchmarks/femc_photon/analysis/femc_photon_plots.py", + output: + results_dir=directory("results/{DETECTOR_CONFIG}/femc_photon"), + shell: + """ +mkdir -p {output.results_dir} +python {input.script} {output.results_dir} +""" diff --git a/benchmarks/femc_photon/analysis/femc_photon_plots.py b/benchmarks/femc_photon/analysis/femc_photon_plots.py new file mode 100644 index 00000000..8170f4a1 --- /dev/null +++ b/benchmarks/femc_photon/analysis/femc_photon_plots.py @@ -0,0 +1,220 @@ +import numpy as np, pandas as pd, matplotlib.pyplot as plt, matplotlib as mpl, awkward as ak, sys, uproot as ur +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 + +import uproot as ur +arrays_sim={p:ur.open(f'sim_output/femc_photon/{config}_rec_photon_{p}GeV.edm4eic.root:events').arrays() for p in (10, 20, 30, 40, 50, 60,70,80)} + +for p in arrays_sim: + array=arrays_sim[p] + tilt=-.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['nclust_endcap']=[len(array['EcalEndcapPClusters.energy'][i]) for i in range(len(array))] + +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) + +#number of clusters +plt.figure() +for eta_min, eta_max, field in (1.5, 2.8, 'nclust_endcap'),: + for p in arrays_sim: + array=arrays_sim[p] + plt.hist(array[field][(array['eta_truth']>eta_min)&(array['eta_truth']eta_min)]['EcalEndcapPClusters.energy']), bins=bins) + bcs=(x[1:]+x[:-1])/2 + + fnc=gauss + slc=abs(bcs-p)<3 + sigma=np.sqrt(y[slc])+0.5*(y[slc]==0) + p0=(100, p, 3) + try: + coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma)) + if abs(coeff[1])>100 or np.sqrt(var_matrix[1][1])>100: + continue + pvals.append(p) + res.append(abs(coeff[2])/coeff[1]) + dres.append(np.sqrt(var_matrix[2][2])/coeff[1]) + scale.append(abs(coeff[1])/p) + dscale.append(np.sqrt(var_matrix[1][1])/p) + except: + pass + plt.sca(axs[0]) + plt.errorbar(pvals, 100*np.array(res), 100*np.array(dres), ls='', marker='o', label=f'${eta_min:.1f}<\\eta<{eta_max:.1f}$') + + plt.sca(axs[1]) + plt.errorbar(pvals, 100*np.array(scale), 100*np.array(dscale), ls='', marker='o', label=f'${eta_min:.1f}<\\eta<{eta_max:.1f}$') + +plt.sca(axs[0]) +plt.legend() +plt.ylim(0) +plt.ylabel("E resolution [%]") +plt.xlabel("E truth [GeV]") + +plt.sca(axs[1]) +plt.ylabel("energy scale [%]") +plt.xlabel("E truth [GeV]") +plt.axhline(100, color='0.5', alpha=0.5, ls='--') +plt.ylim(0, 110) + +plt.tight_layout() +plt.savefig(outdir+"/energy_res_eta_partitions.pdf") diff --git a/benchmarks/femc_photon/analysis/femc_photon_plots.py~ b/benchmarks/femc_photon/analysis/femc_photon_plots.py~ new file mode 100644 index 00000000..4227b4e3 --- /dev/null +++ b/benchmarks/femc_photon/analysis/femc_photon_plots.py~ @@ -0,0 +1,171 @@ +import numpy as np, pandas as pd, matplotlib.pyplot as plt, matplotlib as mpl, awkward as ak, sys, uproot as ur +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 + +import uproot as ur +arrays_sim={p:ur.open(f'sim_output/femc_photon/epic_craterlake_rec_e-_{p}GeV.edm4hep.root:events').arrays() for p in (20, 30, 40, 50, 60,70,80)} + +for p in arrays_sim: + array=arrays_sim[p] + tilt=-.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['nclust_endcap']=[len(array['EcalEndcapPClusters.energy'][i]) for i in range(len(array))] + +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) + +#number of clusters +plt.figure() +for eta_min, eta_max, field in (1.5, 2.8, 'nclust_endcap'),: + for p in arrays_sim: + array=arrays_sim[p] + plt.hist(array[field][(array['eta_truth']>eta_min)&(array['eta_truth'] +#include +#include +#include +#include +#include +#include + +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 = 1000, + 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 + int useCrossingAngle=1 // 0= no rotation, 1 = -25 mrad + ) +{ + 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(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4); + GenParticlePtr p2 = std::make_shared( + 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.*useCrossingAngle; //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( + 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(); + 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; +} diff --git a/benchmarks/femc_photon/config.yml b/benchmarks/femc_photon/config.yml new file mode 100644 index 00000000..a367e509 --- /dev/null +++ b/benchmarks/femc_photon/config.yml @@ -0,0 +1,36 @@ +sim:femc_photon: + stage: simulate + extends: .det_benchmark + parallel: + matrix: + - P: 10 + - P: 20 + - P: 30 + - P: 40 + - P: 50 + - P: 60 + - P: 70 + - P: 80 + timeout: 1 hours + script: + - snakemake --cores 1 sim_output/femc_photon/epic_craterlake_rec_photon_${P}GeV.edm4eic.root + retry: + max: 2 + when: + - runner_system_failure + +bench:femc_photon: + stage: benchmarks + extends: .det_benchmark + needs: ["sim:femc_photon"] + script: + - snakemake --cores 1 results/epic_craterlake/femc_photon + +collect_results:femc_photon: + stage: collect + needs: ["bench:femc_photon"] + script: + - ls -al + - mv results{,_save}/ # move results directory out of the way to preserve it + - snakemake --cores 1 --delete-all-output results/epic_craterlake/femc_photon + - mv results{_save,}/ diff --git a/benchmarks/femc_pi0/Snakefile b/benchmarks/femc_pi0/Snakefile new file mode 100644 index 00000000..4936040c --- /dev/null +++ b/benchmarks/femc_pi0/Snakefile @@ -0,0 +1,58 @@ +rule femc_pi0_generate: + input: + script="benchmarks/femc_pi0/analysis/gen_particles.cxx", + params: + NEVENTS_GEN=1000, + th_max=28, + th_min=2.0 + output: + GEN_FILE="sim_output/femc_pi0/pi0_{P}GeV.hepmc" + shell: + """ +mkdir -p sim_output/femc_pi0 +root -l -b -q '{input.script}({params.NEVENTS_GEN},"{output.GEN_FILE}", "pi0", {params.th_min}, {params.th_max}, 0., 360., {wildcards.P})' +""" + +rule femc_pi0_simulate: + input: + GEN_FILE="sim_output/femc_pi0/pi0_{P}GeV.hepmc" + params: + PHYSICS_LIST="FTFP_BERT" + output: + SIM_FILE="sim_output/femc_pi0/{DETECTOR_CONFIG}_sim_pi0_{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 femc_pi0_recon: + input: + SIM_FILE="sim_output/femc_pi0/{DETECTOR_CONFIG}_sim_pi0_{P}GeV.edm4hep.root" + output: + REC_FILE="sim_output/femc_pi0/{DETECTOR_CONFIG}_rec_pi0_{P}GeV.edm4eic.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,EcalEndcapPInsertRecHits,EcalEndcapPClusters -Pjana:nevents=$NEVENTS_REC +""" + +rule femc_pi0_analysis: + input: + expand("sim_output/femc_pi0/{DETECTOR_CONFIG}_rec_pi0_{P}GeV.edm4eic.root", + P=[10, 20, 30, 40, 50, 60, 70, 80], + DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]), + script="benchmarks/femc_pi0/analysis/femc_pi0_plots.py", + output: + results_dir=directory("results/{DETECTOR_CONFIG}/femc_pi0"), + shell: + """ +mkdir -p {output.results_dir} +python {input.script} {output.results_dir} +""" diff --git a/benchmarks/femc_pi0/analysis/femc_pi0_plots.py b/benchmarks/femc_pi0/analysis/femc_pi0_plots.py new file mode 100644 index 00000000..cdf9d231 --- /dev/null +++ b/benchmarks/femc_pi0/analysis/femc_pi0_plots.py @@ -0,0 +1,257 @@ +import numpy as np, pandas as pd, matplotlib.pyplot as plt, matplotlib as mpl, awkward as ak, sys, uproot as ur +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 + +import uproot as ur +arrays_sim={p:ur.open(f'sim_output/femc_pi0/{config}_rec_pi0_{p}GeV.edm4eic.root:events').arrays() for p in (10, 20, 30, 40, 50, 60,70,80)} + +for p in arrays_sim: + array=arrays_sim[p] + tilt=-.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['nclust_endcap']=[len(array['EcalEndcapPClusters.energy'][i]) for i in range(len(array))] + +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) + +#number of clusters +plt.figure() +for eta_min, eta_max, field in (1.5, 2.8, 'nclust_endcap'),: + for p in arrays_sim: + array=arrays_sim[p] + plt.hist(array[field][(array['eta_truth']>eta_min)&(array['eta_truth']eta_min)&(array['eta_truth']=5))/tot) + plt.errorbar(pvals, f0, label="$n_{cl}$=0") + plt.errorbar(pvals, f1, label="$n_{cl}$=1") + plt.errorbar(pvals, f2, label="$n_{cl}$=2") + plt.errorbar(pvals, f3, label="$n_{cl}$=3") + plt.errorbar(pvals, f4, label="$n_{cl}$=4") + plt.errorbar(pvals, f5p, label="$n_{cl}\\geq$5") + plt.legend() + plt.ylabel("fraction of events") + plt.xlabel("$E_{\\pi^0}$ [GeV]") + plt.ylim(0, 1) + plt.legend(fontsize=20, ncol=2) + plt.savefig(outdir+f"/nclust.pdf") + + + +#number of hits per cluster +fig, axs=plt.subplots(1,2, figsize=(16,8)) +avgs=[] +stds=[] +pvals=[] + +for p in arrays_sim: + + a=arrays_sim[p] + n=[] + nn=-a['EcalEndcapPClusters.hits_begin']+a['EcalEndcapPClusters.hits_end'] + E=a['EcalEndcapPClusters.energy'] + for evt in range(len(array)): + if len(E[evt])==0: + continue + maxE=np.max(E[evt]) + found=False + for i in range(len(E[evt])): + if E[evt][i]==maxE: + n.append(nn[evt][i]) + found=True + break + #if not found: + # n.append(0) + + if p ==50: + plt.sca(axs[0]) + y,x,_=plt.hist(n, range=(0,100), bins=100, histtype='step', label=f"E={p} GeV") + plt.ylabel("events") + plt.xlabel("# hits in cluster") + plt.title(f"$\\pi^0$, E={p} GeV") + pvals.append(p) + avgs.append(np.mean(n)) + stds.append(np.std(n)) + +plt.sca(axs[1]) +plt.errorbar(pvals, avgs, stds, marker='o',ls='') +plt.xlabel("E [GeV]") +plt.ylabel("# hits in cluster [mean$\\pm$std]") +plt.ylim(0) +plt.savefig(outdir+"/nhits_per_cluster.pdf") + + +#energy resolution +def gauss(x, A,mu, sigma): + return A * np.exp(-(x-mu)**2/(2*sigma**2)) +from scipy.optimize import curve_fit + +fig, axs=plt.subplots(1,3, figsize=(24,8)) +pvals=[] +res=[] +dres=[] +scale=[] +dscale=[] +for p in arrays_sim: + bins=np.linspace(15*p/20,22*p/20, 50) + if p==50: + plt.sca(axs[0]) + plt.title(f"E={p} GeV") + y,x,_=plt.hist(np.sum(arrays_sim[p]['EcalEndcapPClusters.energy'], axis=-1), bins=bins, histtype='step') + plt.ylabel("events") + plt.xlabel("$E^{rec}_{\\pi^0}$ [GeV]") + else: + y,x=np.histogram(np.sum(arrays_sim[p]['EcalEndcapPClusters.energy'], axis=-1), bins=bins) + bcs=(x[1:]+x[:-1])/2 + + fnc=gauss + slc=abs(bcs-p)<3 + sigma=np.sqrt(y[slc])+0.5*(y[slc]==0) + p0=(100, p, 3) + + coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma)) + #res=np.abs(coeff[2]/coeff[1]) + if p==50: + xx=np.linspace(15*p/20,22*p/20, 100) + + plt.plot(xx, fnc(xx,*coeff), label=f"$\\sigma_E/E={abs(coeff[2])/coeff[1]*100:.1f}\%$") + plt.axvline(p, color='g', ls='--', alpha=0.7) + plt.legend() + #plt.xlim(0,60) + #plt.show() + pvals.append(p) + res.append(abs(coeff[2])/coeff[1]) + dres.append(np.sqrt(var_matrix[2][2])/coeff[1]) + scale.append(abs(coeff[1])/p) + dscale.append(np.sqrt(var_matrix[1][1])/p) +plt.sca(axs[1]) +plt.errorbar(pvals, 100*np.array(res), 100*np.array(dres), ls='', marker='o') +fnc = lambda E, a, b: np.hypot(a,b/np.sqrt(E)) +p0=(.05, .12) +coeff, var_matrix = curve_fit(fnc, pvals, res, p0=p0,sigma=dres) +xx=np.linspace(15, 85, 100) +plt.plot(xx, 100*fnc(xx,*coeff), label=f'fit:{100*coeff[0]:.1f}%$\\oplus\\frac{{{100*coeff[1]:.0f}\\%}}{{\\sqrt{{E}}}}$') +plt.legend() +plt.ylim(0) +plt.ylabel("E resolution [%]") +plt.xlabel("E truth [GeV]") +plt.sca(axs[2]) + +plt.errorbar(pvals, 100*np.array(scale), 100*np.array(dscale), ls='', marker='o') +plt.ylabel("energy scale [%]") +plt.xlabel("E truth [GeV]") +plt.axhline(100, color='0.5', alpha=0.5, ls='--') +plt.ylim(0, 110) +plt.tight_layout() +plt.savefig(outdir+"/energy_res.pdf") + +#energy res in eta bins +fig, axs=plt.subplots(1,2, figsize=(16,8)) + +partitions=[1.5, 2.0, 2.5, 3.0] +for eta_min, eta_max in zip(partitions[:-1], partitions[1:]): + pvals=[] + res=[] + dres=[] + scale=[] + dscale=[] + for p in arrays_sim: + bins=np.linspace(15*p/20,22*p/20, 50) + eta_truth=arrays_sim[p]['eta_truth'] + y,x=np.histogram(np.sum(arrays_sim[p][(eta_trutheta_min)]['EcalEndcapPClusters.energy'], axis=-1), bins=bins) + bcs=(x[1:]+x[:-1])/2 + + fnc=gauss + slc=abs(bcs-p)<3 + sigma=np.sqrt(y[slc])+0.5*(y[slc]==0) + p0=(100, p, 3) + + try: + coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma)) + if abs(coeff[1])>100 or np.sqrt(var_matrix[1][1])>100: + continue + pvals.append(p) + res.append(abs(coeff[2])/coeff[1]) + dres.append(np.sqrt(var_matrix[2][2])/coeff[1]) + scale.append(abs(coeff[1])/p) + dscale.append(np.sqrt(var_matrix[1][1])/p) + except: + pass + plt.sca(axs[0]) + plt.errorbar(pvals, 100*np.array(res), 100*np.array(dres), ls='', marker='o', label=f'${eta_min:.1f}<\\eta<{eta_max:.1f}$') + + plt.sca(axs[1]) + plt.errorbar(pvals, 100*np.array(scale), 100*np.array(dscale), ls='', marker='o', label=f'${eta_min:.1f}<\\eta<{eta_max:.1f}$') + +plt.sca(axs[0]) +plt.legend() +plt.ylim(0) +plt.ylabel("E resolution [%]") +plt.xlabel("E truth [GeV]") + +plt.sca(axs[1]) +plt.ylabel("energy scale [%]") +plt.xlabel("E truth [GeV]") +plt.axhline(100, color='0.5', alpha=0.5, ls='--') +plt.ylim(0, 110) + +plt.tight_layout() +plt.savefig(outdir+"/energy_res_eta_partitions.pdf") diff --git a/benchmarks/femc_pi0/analysis/gen_particles.cxx b/benchmarks/femc_pi0/analysis/gen_particles.cxx new file mode 100644 index 00000000..0877a852 --- /dev/null +++ b/benchmarks/femc_pi0/analysis/gen_particles.cxx @@ -0,0 +1,127 @@ +#include "HepMC3/GenEvent.h" +#include "HepMC3/ReaderAscii.h" +#include "HepMC3/WriterAscii.h" +#include "HepMC3/Print.h" + +#include "TRandom3.h" +#include "TVector3.h" + +#include +#include +#include +#include +#include +#include +#include + +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 = 1000, + 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 + int useCrossingAngle=1 // 0= no rotation, 1 = -25 mrad + ) +{ + 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(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4); + GenParticlePtr p2 = std::make_shared( + 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.*useCrossingAngle; //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( + 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(); + 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; +} diff --git a/benchmarks/femc_pi0/config.yml b/benchmarks/femc_pi0/config.yml new file mode 100644 index 00000000..bcbd90af --- /dev/null +++ b/benchmarks/femc_pi0/config.yml @@ -0,0 +1,36 @@ +sim:femc_pi0: + stage: simulate + extends: .det_benchmark + parallel: + matrix: + - P: 10 + - P: 20 + - P: 30 + - P: 40 + - P: 50 + - P: 60 + - P: 70 + - P: 80 + timeout: 1 hours + script: + - snakemake --cores 1 sim_output/femc_pi0/epic_craterlake_rec_pi0_${P}GeV.edm4eic.root + retry: + max: 2 + when: + - runner_system_failure + +bench:femc_pi0: + stage: benchmarks + extends: .det_benchmark + needs: ["sim:femc_pi0"] + script: + - snakemake --cores 1 results/epic_craterlake/femc_pi0 + +collect_results:femc_pi0: + stage: collect + needs: ["bench:femc_pi0"] + script: + - ls -al + - mv results{,_save}/ # move results directory out of the way to preserve it + - snakemake --cores 1 --delete-all-output results/epic_craterlake/femc_pi0 + - mv results{_save,}/