From fda7e5bef9098b2430dce5b927d526f772022eb5 Mon Sep 17 00:00:00 2001 From: Sebouh Paul Date: Mon, 23 Sep 2024 18:23:47 -0400 Subject: [PATCH 01/24] Muon insert (#71) Co-authored-by: Dmitry Kalinkin --- .gitlab-ci.yml | 3 +- Snakefile | 2 + benchmarks/insert_muon/Snakefile | 57 ++++++++ .../insert_muon/analysis/gen_particles.cxx | 127 ++++++++++++++++++ benchmarks/insert_muon/analysis/muon_plots.py | 86 ++++++++++++ benchmarks/insert_muon/config.yml | 30 +++++ 6 files changed, 304 insertions(+), 1 deletion(-) create mode 100644 benchmarks/insert_muon/Snakefile create mode 100644 benchmarks/insert_muon/analysis/gen_particles.cxx create mode 100644 benchmarks/insert_muon/analysis/muon_plots.py create mode 100644 benchmarks/insert_muon/config.yml diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 4b3c9c8e..fb5b57e3 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -147,11 +147,11 @@ include: - local: 'benchmarks/pid/config.yml' - local: 'benchmarks/timing/config.yml' - local: 'benchmarks/b0_tracker/config.yml' + - local: 'benchmarks/insert_muon/config.yml' - local: 'benchmarks/zdc_sigma/config.yml' - local: 'benchmarks/zdc_lambda/config.yml' - local: 'benchmarks/insert_neutron/config.yml' - deploy_results: allow_failure: true stage: deploy @@ -171,6 +171,7 @@ deploy_results: - "collect_results:tracking_performances_dis" - "collect_results:zdc" - "collect_results:zdc_lyso" + - "collect_results:insert_muon" - "collect_results:zdc_photon" - "collect_results:zdc_pi0" script: diff --git a/Snakefile b/Snakefile index 58357199..bce37f8d 100644 --- a/Snakefile +++ b/Snakefile @@ -7,6 +7,8 @@ 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" include: "benchmarks/zdc_photon/Snakefile" diff --git a/benchmarks/insert_muon/Snakefile b/benchmarks/insert_muon/Snakefile new file mode 100644 index 00000000..b36026b0 --- /dev/null +++ b/benchmarks/insert_muon/Snakefile @@ -0,0 +1,57 @@ +rule insert_muon_generate: + input: + script="benchmarks/insert_muon/analysis/gen_particles.cxx", + params: + NEVENTS_GEN=5000, + th_max=7.0, + th_min=1.7 + output: + GEN_FILE="sim_output/insert_muon/mu-_{P}GeV.hepmc" + shell: + """ +root -l -b -q '{input.script}({params.NEVENTS_GEN},"{output.GEN_FILE}", "mu-", {params.th_min}, {params.th_max}, 0., 360., {wildcards.P})' +""" + +rule insert_muon_simulate: + input: + GEN_FILE="sim_output/insert_muon/mu-_{P}GeV.hepmc" + params: + PHYSICS_LIST="FTFP_BERT" + output: + SIM_FILE="sim_output/insert_muon/{DETECTOR_CONFIG}_sim_mu-_{P}GeV.edm4hep.root" + shell: + """ +NEVENTS_SIM=5000 +# 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_muon_recon: + input: + SIM_FILE="sim_output/insert_muon/{DETECTOR_CONFIG}_sim_mu-_{P}GeV.edm4hep.root" + output: + REC_FILE="sim_output/insert_muon/{DETECTOR_CONFIG}_rec_mu-_{P}GeV.edm4hep.root" + shell: + """ +NEVENTS_REC=5000 +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", + P=[50], + DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]), + script="benchmarks/insert_muon/analysis/muon_plots.py", + output: + results_dir=directory("results/{DETECTOR_CONFIG}/insert_muon"), + shell: + """ +mkdir -p {output.results_dir} +python {input.script} {output.results_dir} +""" diff --git a/benchmarks/insert_muon/analysis/gen_particles.cxx b/benchmarks/insert_muon/analysis/gen_particles.cxx new file mode 100644 index 00000000..0877a852 --- /dev/null +++ b/benchmarks/insert_muon/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/insert_muon/analysis/muon_plots.py b/benchmarks/insert_muon/analysis/muon_plots.py new file mode 100644 index 00000000..4c45aac3 --- /dev/null +++ b/benchmarks/insert_muon/analysis/muon_plots.py @@ -0,0 +1,86 @@ +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") diff --git a/benchmarks/insert_muon/config.yml b/benchmarks/insert_muon/config.yml new file mode 100644 index 00000000..bf69f555 --- /dev/null +++ b/benchmarks/insert_muon/config.yml @@ -0,0 +1,30 @@ +sim:insert_muon: + stage: simulate + extends: .det_benchmark + needs: ["common:detector"] + parallel: + matrix: + - P: 50 + timeout: 1 hours + script: + - snakemake --cores 1 sim_output/insert_muon/epic_craterlake_sim_mu-_${P}GeV.edm4hep.root + retry: + max: 2 + when: + - runner_system_failure + +bench:insert_muon: + stage: benchmarks + extends: .det_benchmark + needs: ["sim:insert_muon"] + script: + - snakemake --cores 1 results/epic_craterlake/insert_muon + +collect_results:insert_muon: + stage: collect + needs: ["bench:insert_muon"] + 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/insert_muon + - mv results{_save,}/ From 58f4b0a0b2e078581ba296969c2523b3b7f7a396 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Tue, 24 Sep 2024 06:17:28 -0400 Subject: [PATCH 02/24] .gitlab-ci.yml: don't build detector (#69) --- .gitlab-ci.yml | 13 +++---------- benchmarks/insert_muon/config.yml | 1 - benchmarks/zdc_photon/config.yml | 1 - benchmarks/zdc_pi0/config.yml | 1 - 4 files changed, 3 insertions(+), 13 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index fb5b57e3..a1dc5d8c 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -3,7 +3,6 @@ image: ${BENCHMARKS_REGISTRY}/${BENCHMARKS_CONTAINER}:${BENCHMARKS_TAG} variables: DETECTOR: epic DETECTOR_CONFIG: epic_craterlake - DETECTOR_REPOSITORYURL: 'https://github.com/eic/epic.git' GITHUB_SHA: '' GITHUB_REPOSITORY: '' @@ -89,13 +88,6 @@ common:setup: echo "BENCHMARKS_CONTAINER: ${BENCHMARKS_CONTAINER}" echo "BENCHMARKS_REGISTRY: ${BENCHMARKS_REGISTRY}" - source setup/bin/env.sh && ./setup/bin/install_common.sh - - -common:detector: - stage: initialize - needs: ["common:setup"] - script: - - source .local/bin/env.sh && build_detector.sh - mkdir_local_data_link sim_output - mkdir -p results - mkdir -p config @@ -103,7 +95,7 @@ common:detector: get_data: stage: data_init - needs: ["common:detector"] + needs: ["common:setup"] script: - source .local/bin/env.sh - ln -s "${LOCAL_DATA_PATH}/sim_output" sim_output @@ -112,10 +104,11 @@ get_data: .det_benchmark: needs: - - ["get_data","common:detector"] + - ["get_data","common:setup"] before_script: - mc config host add S3 https://eics3.sdcc.bnl.gov:9000 ${S3_ACCESS_KEY} ${S3_SECRET_KEY} - source .local/bin/env.sh + - source /opt/detector/epic-main/setup.sh - ls -lrtha - ln -s "${LOCAL_DATA_PATH}/sim_output" sim_output - ln -s "${LOCAL_DATA_PATH}/datasets/data" data diff --git a/benchmarks/insert_muon/config.yml b/benchmarks/insert_muon/config.yml index bf69f555..6f7681f1 100644 --- a/benchmarks/insert_muon/config.yml +++ b/benchmarks/insert_muon/config.yml @@ -1,7 +1,6 @@ sim:insert_muon: stage: simulate extends: .det_benchmark - needs: ["common:detector"] parallel: matrix: - P: 50 diff --git a/benchmarks/zdc_photon/config.yml b/benchmarks/zdc_photon/config.yml index c8e8bdb5..946650a3 100644 --- a/benchmarks/zdc_photon/config.yml +++ b/benchmarks/zdc_photon/config.yml @@ -1,7 +1,6 @@ sim:zdc_photon: stage: simulate extends: .det_benchmark - needs: ["common:detector"] parallel: matrix: - P: 20 diff --git a/benchmarks/zdc_pi0/config.yml b/benchmarks/zdc_pi0/config.yml index 4059d236..a2fcb5de 100644 --- a/benchmarks/zdc_pi0/config.yml +++ b/benchmarks/zdc_pi0/config.yml @@ -1,7 +1,6 @@ sim:zdc_pi0: stage: simulate extends: .det_benchmark - needs: ["common:detector"] parallel: matrix: - P: 60 From f04a3693fe86c04b4314ecc25df4778431da8ecc Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Fri, 27 Sep 2024 09:01:38 -0400 Subject: [PATCH 03/24] 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 04/24] 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 05/24] 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,}/ From 8116047447bcfa6aa493ea3dd4f6ab3e724b5c7f Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Mon, 7 Oct 2024 10:50:28 -0400 Subject: [PATCH 06/24] zdc_lyso: set maxfev=10000 for curve_fit (#77) --- benchmarks/zdc_lyso/analysis/analysis.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/benchmarks/zdc_lyso/analysis/analysis.py b/benchmarks/zdc_lyso/analysis/analysis.py index eace6112..8994fa24 100644 --- a/benchmarks/zdc_lyso/analysis/analysis.py +++ b/benchmarks/zdc_lyso/analysis/analysis.py @@ -82,7 +82,7 @@ def rotateY(xdata, zdata, angle): hist, x = np.histogram(temp,bins=np.linspace(min(temp),max(temp)+np.std(abs(temp)),2*int(np.sqrt(len(temp))))) x = x[1:]/2 + x[:-1]/2 plt.errorbar(x,hist,yerr=np.sqrt(hist),fmt='-o',label='Cluster') - coeff, covar = curve_fit(gaussian,x[1:],hist[1:],p0=(max(hist[x>=np.std(abs(temp))]),np.mean(temp[temp!=0]),np.std(temp[temp!=0]))) + coeff, covar = curve_fit(gaussian,x[1:],hist[1:],p0=(max(hist[x>=np.std(abs(temp))]),np.mean(temp[temp!=0]),np.std(temp[temp!=0])),maxfev=10000) #plt.plot(np.linspace(coeff[1]-3*coeff[2],coeff[1]+3*coeff[2],50),gaussian(np.linspace(coeff[1]-3*coeff[2],coeff[1]+3*coeff[2],50),*coeff)) mu.append(coeff[1]) sigma.append(coeff[2]) @@ -91,7 +91,7 @@ def rotateY(xdata, zdata, angle): hist, x = np.histogram(temp,bins=np.linspace(min(temp),max(temp)+np.std(abs(temp)),2*int(np.sqrt(len(temp))))) x = x[1:]/2 + x[:-1]/2 plt.errorbar(x,hist,yerr=np.sqrt(hist),fmt='-o',label='Digitization') - coeff, covar = curve_fit(gaussian,x[1:],hist[1:],p0=(max(hist[x>=np.std(abs(temp))]),np.mean(temp[temp!=0]),np.std(temp[temp!=0]))) + coeff, covar = curve_fit(gaussian,x[1:],hist[1:],p0=(max(hist[x>=np.std(abs(temp))]),np.mean(temp[temp!=0]),np.std(temp[temp!=0])),maxfev=10000) #plt.plot(np.linspace(coeff[1]-3*coeff[2],coeff[1]+3*coeff[2],50),gaussian(np.linspace(coeff[1]-3*coeff[2],coeff[1]+3*coeff[2],50),*coeff)) mu.append(coeff[1]) sigma.append(coeff[2]) @@ -100,7 +100,7 @@ def rotateY(xdata, zdata, angle): hist, x = np.histogram(temp,bins=np.linspace(min(temp),max(temp)+np.std(abs(temp)),2*int(np.sqrt(len(temp))))) x = x[1:]/2 + x[:-1]/2 plt.errorbar(x,hist,yerr=np.sqrt(hist),fmt='-o',label='Simulation') - coeff, covar = curve_fit(gaussian,x[1:],hist[1:],p0=(max(hist[x>=np.std(abs(temp))]),np.mean(temp[temp!=0]),np.std(temp[temp!=0]))) + coeff, covar = curve_fit(gaussian,x[1:],hist[1:],p0=(max(hist[x>=np.std(abs(temp))]),np.mean(temp[temp!=0]),np.std(temp[temp!=0])),maxfev=10000) #plt.plot(np.linspace(coeff[1]-3*coeff[2],coeff[1]+3*coeff[2],50),gaussian(np.linspace(coeff[1]-3*coeff[2],coeff[1]+3*coeff[2],50),*coeff)) mu.append(coeff[1]) sigma.append(coeff[2]) From 8de512aaa6642aedcdd0d0d29eb42353cd14568d Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Mon, 7 Oct 2024 18:39:00 -0400 Subject: [PATCH 07/24] Snakefile: fix caching in fetch_epic (#78) --- Snakefile | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Snakefile b/Snakefile index 54323284..9fbb633a 100644 --- a/Snakefile +++ b/Snakefile @@ -34,6 +34,9 @@ def get_remote_path(path): rule fetch_epic: output: filepath="EPIC/{PATH}" + params: + # wildcards are not included in hash for caching, we need to add them as params + PATH=lambda wildcards: wildcards.PATH cache: True retries: 3 shell: """ From 8c7d0fda39839c76d52ac683de0be4eecbc7cb89 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Tue, 8 Oct 2024 22:28:52 -0400 Subject: [PATCH 08/24] **/config.yml: collect_results shoud extend from .det_benchmark (#80) --- benchmarks/femc_electron/config.yml | 5 +++-- benchmarks/femc_photon/config.yml | 5 +++-- benchmarks/femc_pi0/config.yml | 6 +++--- benchmarks/insert_muon/config.yml | 5 +++-- benchmarks/insert_neutron/config.yml | 5 +++-- benchmarks/zdc_lambda/config.yml | 5 +++-- benchmarks/zdc_photon/config.yml | 5 +++-- benchmarks/zdc_pi0/config.yml | 7 +++---- benchmarks/zdc_sigma/config.yml | 5 +++-- 9 files changed, 27 insertions(+), 21 deletions(-) diff --git a/benchmarks/femc_electron/config.yml b/benchmarks/femc_electron/config.yml index a600454d..7b6f1043 100644 --- a/benchmarks/femc_electron/config.yml +++ b/benchmarks/femc_electron/config.yml @@ -1,6 +1,6 @@ sim:femc_electron: - stage: simulate extends: .det_benchmark + stage: simulate parallel: matrix: - P: 10 @@ -20,13 +20,14 @@ sim:femc_electron: - runner_system_failure bench:femc_electron: - stage: benchmarks extends: .det_benchmark + stage: benchmarks needs: ["sim:femc_electron"] script: - snakemake --cores 1 results/epic_craterlake/femc_electron collect_results:femc_electron: + extends: .det_benchmark stage: collect needs: ["bench:femc_electron"] script: diff --git a/benchmarks/femc_photon/config.yml b/benchmarks/femc_photon/config.yml index a367e509..af710897 100644 --- a/benchmarks/femc_photon/config.yml +++ b/benchmarks/femc_photon/config.yml @@ -1,6 +1,6 @@ sim:femc_photon: - stage: simulate extends: .det_benchmark + stage: simulate parallel: matrix: - P: 10 @@ -20,13 +20,14 @@ sim:femc_photon: - runner_system_failure bench:femc_photon: - stage: benchmarks extends: .det_benchmark + stage: benchmarks needs: ["sim:femc_photon"] script: - snakemake --cores 1 results/epic_craterlake/femc_photon collect_results:femc_photon: + extends: .det_benchmark stage: collect needs: ["bench:femc_photon"] script: diff --git a/benchmarks/femc_pi0/config.yml b/benchmarks/femc_pi0/config.yml index bcbd90af..46e6ed9e 100644 --- a/benchmarks/femc_pi0/config.yml +++ b/benchmarks/femc_pi0/config.yml @@ -1,6 +1,6 @@ sim:femc_pi0: - stage: simulate extends: .det_benchmark + stage: simulate parallel: matrix: - P: 10 @@ -11,7 +11,6 @@ sim:femc_pi0: - 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: @@ -20,13 +19,14 @@ sim:femc_pi0: - runner_system_failure bench:femc_pi0: - stage: benchmarks extends: .det_benchmark + stage: benchmarks needs: ["sim:femc_pi0"] script: - snakemake --cores 1 results/epic_craterlake/femc_pi0 collect_results:femc_pi0: + extends: .det_benchmark stage: collect needs: ["bench:femc_pi0"] script: diff --git a/benchmarks/insert_muon/config.yml b/benchmarks/insert_muon/config.yml index 252cd6e8..76899c1d 100644 --- a/benchmarks/insert_muon/config.yml +++ b/benchmarks/insert_muon/config.yml @@ -1,6 +1,6 @@ sim:insert_muon: - stage: simulate extends: .det_benchmark + stage: simulate parallel: matrix: - P: 50 @@ -12,13 +12,14 @@ sim:insert_muon: - runner_system_failure bench:insert_muon: - stage: benchmarks extends: .det_benchmark + stage: benchmarks needs: ["sim:insert_muon"] script: - snakemake --cores 1 results/epic_craterlake/insert_muon collect_results:insert_muon: + extends: .det_benchmark stage: collect needs: ["bench:insert_muon"] script: diff --git a/benchmarks/insert_neutron/config.yml b/benchmarks/insert_neutron/config.yml index 4c723e68..eedac937 100644 --- a/benchmarks/insert_neutron/config.yml +++ b/benchmarks/insert_neutron/config.yml @@ -1,6 +1,6 @@ sim:insert_neutron: - stage: simulate extends: .det_benchmark + stage: simulate parallel: matrix: - P: 20 @@ -18,13 +18,14 @@ sim:insert_neutron: - runner_system_failure bench:insert_neutron: - stage: benchmarks extends: .det_benchmark + stage: benchmarks needs: ["sim:insert_neutron"] script: - snakemake --cores 1 results/epic_craterlake/insert_neutron collect_results:insert_neutron: + extends: .det_benchmark stage: collect needs: ["bench:insert_neutron"] script: diff --git a/benchmarks/zdc_lambda/config.yml b/benchmarks/zdc_lambda/config.yml index e62b0d9c..e6b16031 100644 --- a/benchmarks/zdc_lambda/config.yml +++ b/benchmarks/zdc_lambda/config.yml @@ -1,6 +1,6 @@ sim:zdc_lambda: - stage: simulate extends: .det_benchmark + stage: simulate parallel: matrix: - P: 100 @@ -19,13 +19,14 @@ sim:zdc_lambda: - runner_system_failure bench:zdc_lambda: - stage: benchmarks extends: .det_benchmark + stage: benchmarks needs: ["sim:zdc_lambda"] script: - snakemake --cores 1 results/epic_zdc_sipm_on_tile_only/zdc_lambda collect_results:zdc_lambda: + extends: .det_benchmark stage: collect needs: ["bench:zdc_lambda"] script: diff --git a/benchmarks/zdc_photon/config.yml b/benchmarks/zdc_photon/config.yml index 4d205eb2..aee4f6c6 100644 --- a/benchmarks/zdc_photon/config.yml +++ b/benchmarks/zdc_photon/config.yml @@ -1,6 +1,6 @@ sim:zdc_photon: - stage: simulate extends: .det_benchmark + stage: simulate parallel: matrix: - P: 20 @@ -19,14 +19,15 @@ sim:zdc_photon: - runner_system_failure bench:zdc_photon: - stage: benchmarks extends: .det_benchmark + stage: benchmarks needs: ["sim:zdc_photon"] script: - mkdir -p results/epic_zdc_sipm_on_tile_only - python benchmarks/zdc_photon/analysis/zdc_photon_plots.py results/epic_zdc_sipm_on_tile_only/zdc_photon collect_results:zdc_photon: + extends: .det_benchmark stage: collect extends: .det_benchmark needs: ["bench:zdc_photon"] diff --git a/benchmarks/zdc_pi0/config.yml b/benchmarks/zdc_pi0/config.yml index ae2c1fd9..f0b8993e 100644 --- a/benchmarks/zdc_pi0/config.yml +++ b/benchmarks/zdc_pi0/config.yml @@ -1,6 +1,6 @@ sim:zdc_pi0: - stage: simulate extends: .det_benchmark + stage: simulate parallel: matrix: - P: 60 @@ -16,17 +16,16 @@ sim:zdc_pi0: - runner_system_failure bench:zdc_pi0: - stage: benchmarks - extends: .det_benchmark + stage: benchmarks needs: ["sim:zdc_pi0"] script: - mkdir -p results/epic_zdc_sipm_on_tile_only - python benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py results/epic_zdc_sipm_on_tile_only/zdc_pi0 collect_results:zdc_pi0: - stage: collect extends: .det_benchmark + stage: collect needs: ["bench:zdc_pi0"] script: - ls -al diff --git a/benchmarks/zdc_sigma/config.yml b/benchmarks/zdc_sigma/config.yml index 5a490510..af204d42 100644 --- a/benchmarks/zdc_sigma/config.yml +++ b/benchmarks/zdc_sigma/config.yml @@ -1,6 +1,6 @@ sim:zdc_sigma: - stage: simulate extends: .det_benchmark + stage: simulate parallel: matrix: - P: 100 @@ -19,13 +19,14 @@ sim:zdc_sigma: - runner_system_failure bench:zdc_sigma: - stage: benchmarks extends: .det_benchmark + stage: benchmarks needs: ["sim:zdc_sigma"] script: - snakemake --cores 1 results/epic_zdc_sipm_on_tile_only/zdc_sigma collect_results:zdc_sigma: + extends: .det_benchmark stage: collect needs: ["bench:zdc_sigma"] script: From a8c6d875240d46ab38944f12388119d0d7a6b029 Mon Sep 17 00:00:00 2001 From: Peter Steinberg Date: Thu, 10 Oct 2024 13:20:31 -0400 Subject: [PATCH 09/24] first commit of LFHCAL benchmarks (#48) Co-authored-by: Dmitry Kalinkin --- .gitlab-ci.yml | 1 + Snakefile | 1 + benchmarks/lfhcal/LFHCAL_Performance.C | 156 ++++++++++++++++++ benchmarks/lfhcal/README.md | 11 ++ benchmarks/lfhcal/Snakefile | 145 +++++++++++++++++ benchmarks/lfhcal/config.yml | 48 ++++++ benchmarks/lfhcal/doCompare_widebins_mom.C | 176 +++++++++++++++++++++ 7 files changed, 538 insertions(+) create mode 100644 benchmarks/lfhcal/LFHCAL_Performance.C create mode 100644 benchmarks/lfhcal/README.md create mode 100644 benchmarks/lfhcal/Snakefile create mode 100644 benchmarks/lfhcal/config.yml create mode 100644 benchmarks/lfhcal/doCompare_widebins_mom.C diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index a56c79d2..f25c74e9 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -132,6 +132,7 @@ include: - local: 'benchmarks/tracking_performances_dis/config.yml' - local: 'benchmarks/barrel_ecal/config.yml' - local: 'benchmarks/barrel_hcal/config.yml' + - local: 'benchmarks/lfhcal/config.yml' - local: 'benchmarks/zdc/config.yml' - local: 'benchmarks/zdc_lyso/config.yml' - local: 'benchmarks/zdc_photon/config.yml' diff --git a/Snakefile b/Snakefile index 9fbb633a..1390e867 100644 --- a/Snakefile +++ b/Snakefile @@ -7,6 +7,7 @@ include: "benchmarks/ecal_gaps/Snakefile" include: "benchmarks/material_scan/Snakefile" include: "benchmarks/tracking_performances/Snakefile" include: "benchmarks/tracking_performances_dis/Snakefile" +include: "benchmarks/lfhcal/Snakefile" include: "benchmarks/insert_muon/Snakefile" include: "benchmarks/zdc_lambda/Snakefile" include: "benchmarks/zdc_lyso/Snakefile" diff --git a/benchmarks/lfhcal/LFHCAL_Performance.C b/benchmarks/lfhcal/LFHCAL_Performance.C new file mode 100644 index 00000000..dcfb7595 --- /dev/null +++ b/benchmarks/lfhcal/LFHCAL_Performance.C @@ -0,0 +1,156 @@ +// Code to extract the Tracking Performances +// Shyam Kumar; INFN Bari, Italy +// shyam.kumar@ba.infn.it; shyam.kumar@cern.ch + +#include "TGraphErrors.h" +#include "TF1.h" +#include "TRandom.h" +#include "TCanvas.h" +#include "TLegend.h" +#include "TMath.h" +#include "TVector3.h" + +#define mpi 0.139 // 1.864 GeV/c^2 + +void LFHCAL_Performance(TString filename="tracking_output",TString particle="pi-", double mom=0.1, Double_t pTcut = 0.0, TString name = "") +{ + + // style of the plot + gStyle->SetPalette(1); + gStyle->SetOptTitle(1); + gStyle->SetTitleOffset(.85,"X");gStyle->SetTitleOffset(.85,"Y"); + gStyle->SetTitleSize(.05,"X");gStyle->SetTitleSize(.05,"Y"); + gStyle->SetLabelSize(.04,"X");gStyle->SetLabelSize(.04,"Y"); + gStyle->SetHistLineWidth(2); + gStyle->SetOptFit(1); + gStyle->SetOptStat(1); + + TString dir = ""; + TString dist_dir_mom = "mom_resol"; + + bool debug=true; + // Tree with reconstructed tracks + const int nbins_eta = 5; + int nfiles = 100; + double eta[nbins_eta+1]={1.2,1.5,2,2.5,3,3.5}; + double pt[nbins_eta+1]={0.5,1.0,2.0,5.0,10.0,20.1}; + TH1D *histp[nbins_eta]; + + + for (int i=0; iSetTitle(Form("%1.1f < #eta < %1.1f && p = %1.1f ",eta[i],eta[i+1],mom)); + histp[i]->SetName(Form("hist_mom_%1.1f_%1.1f_pmax_%1.1f",mom,eta[i],eta[i+1])); + } + + TFile* file = TFile::Open(filename.Data()); + if (!file) {printf("file not found !!!"); return;} + TTreeReader myReader("events", file); // name of tree and file + if (debug) cout<<"Filename: "<GetName()<<"\t NEvents: "< charge(myReader, "MCParticles.charge"); + TTreeReaderArray vx_mc(myReader, "MCParticles.vertex.x"); + TTreeReaderArray vy_mc(myReader, "MCParticles.vertex.y"); + TTreeReaderArray vz_mc(myReader, "MCParticles.vertex.z"); + TTreeReaderArray px_mc(myReader, "MCParticles.momentum.x"); + TTreeReaderArray py_mc(myReader, "MCParticles.momentum.y"); + TTreeReaderArray pz_mc(myReader, "MCParticles.momentum.z"); + TTreeReaderArray status(myReader, "MCParticles.generatorStatus"); + TTreeReaderArray pdg(myReader, "MCParticles.PDG"); + + TTreeReaderArray pe_lc(myReader, "LFHCALClusters.energy"); + TTreeReaderArray px_lc(myReader, "LFHCALClusters.position.x"); + TTreeReaderArray py_lc(myReader, "LFHCALClusters.position.y"); + TTreeReaderArray pz_lc(myReader, "LFHCALClusters.position.z"); + TTreeReaderArray pe_ec(myReader, "EcalEndcapPClusters.energy"); + TTreeReaderArray px_ec(myReader, "EcalEndcapPClusters.position.x"); + TTreeReaderArray py_ec(myReader, "EcalEndcapPClusters.position.y"); + TTreeReaderArray pz_ec(myReader, "EcalEndcapPClusters.position.z"); + + int count =0; + int matchId = 1; // Always matched track assigned the index 0 + while (myReader.Next()) + { + std::cout << "events = " << count++ << std::endl; + for (int j = 0; j < pdg.GetSize(); ++j) + { + if (status[j] !=1 && pdg.GetSize()!=1) continue; + Double_t pzmc = pz_mc[j]; + Double_t ptmc = sqrt(px_mc[j]*px_mc[j]+py_mc[j]*py_mc[j]); + Double_t pmc = sqrt(px_mc[j]*px_mc[j]+py_mc[j]*py_mc[j]+pz_mc[j]*pz_mc[j]); // 1./(q/p); similar to prec + Double_t etamc = -1.0*TMath::Log(TMath::Tan((TMath::ACos(pzmc/fabs(pmc)))/2)); + Double_t phimc = TMath::ATan2(py_mc[j],px_mc[j]); + std::cout << "neutron p=" << pmc << " pt=" << ptmc << std::endl; + + if (fabs(ptmc) < pTcut) continue; + + float l_px_tot=0; + float l_py_tot=0; + float l_pz_tot=0; + float l_e_tot=0; + + std::cout << "LFHCAL nclus=" << px_lc.GetSize() << " ECAL nclus=" << pe_ec.GetSize() << std::endl; + for (int jl = 0;jlcd(); + for (int ibin=0; ibinWrite(); + fout_mom->Close(); + +} + + + + diff --git a/benchmarks/lfhcal/README.md b/benchmarks/lfhcal/README.md new file mode 100644 index 00000000..5f56298b --- /dev/null +++ b/benchmarks/lfhcal/README.md @@ -0,0 +1,11 @@ +Detector Benchmark for LFHCAL +=============================== + +## Overview +This benchmark generates events with a range of particle species at fixed energies, all aimed at, or forward of, the LFHCAL, and produces basic plots of clustering performance, both in isolation and combined with the FEMCAL. + + +## Contacts +[Peter Steinberg](steinberg@bnl.gov) + + diff --git a/benchmarks/lfhcal/Snakefile b/benchmarks/lfhcal/Snakefile new file mode 100644 index 00000000..459392c0 --- /dev/null +++ b/benchmarks/lfhcal/Snakefile @@ -0,0 +1,145 @@ +rule lfhcal_sim: + input: + steering_file="EPIC/EVGEN/SINGLE/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.steer", + warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root", + output: + "sim_output/lfhcal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.edm4hep.root", + log: + "sim_output/lfhcal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.edm4hep.root.log", + wildcard_constraints: + ENERGY="[0-9]+[kMG]eV", + PARTICLE="(neutron|pi-|gamma)", + PHASE_SPACE="3to50deg", + INDEX="\d{4}", + params: + N_EVENTS=1000 + shell: + """ +ddsim \ + --runType batch \ + --enableGun \ + --steeringFile "{input.steering_file}" \ + --random.seed 1{wildcards.INDEX} \ + --filter.tracker edep0 \ + -v INFO \ + --numberOfEvents {params.N_EVENTS} \ + --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \ + --outputFile {output} +""" + + +rule lfhcal_recon: + input: + "sim_output/lfhcal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.edm4hep.root", + output: + "sim_output/lfhcal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.eicrecon.tree.edm4eic.root", + log: + "sim_output/lfhcal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.eicrecon.tree.edm4eic.root.log", + wildcard_constraints: + INDEX="\d{4}", + shell: """ +env DETECTOR_CONFIG={wildcards.DETECTOR_CONFIG} \ + eicrecon {input} -Ppodio:output_file={output} \ + -Ppodio:output_collections=MCParticles,ReconstructedParticles,LFHCALTruthClusters,LFHCALClusters,LFHCALHits,EcalEndcapPTruthClusters,EcalEndcapPClusters,EcalEndcapPHits +""" + +rule lfhcal_at_momentum: + input: + script="benchmarks/lfhcal/LFHCAL_Performance.C", + # TODO pass as a file list? + sim=lambda wildcards: + expand( + "sim_output/lfhcal/{DETECTOR_CONFIG}/{{PARTICLE}}/{ENERGY}/{PHASE_SPACE}/{{PARTICLE}}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.eicrecon.tree.edm4eic.root", + DETECTOR_CONFIG="epic_craterlake", + ENERGY=f"{float(wildcards.MOMENTUM):.0f}GeV" if float(wildcards.MOMENTUM) >= 1 else f"{float(wildcards.MOMENTUM) * 1000:.0f}MeV", + PHASE_SPACE=["3to50deg"], + INDEX=range(1), + ) + if wildcards.CAMPAIGN == "local" else + ancient(expand( + "EPIC/RECO/{{CAMPAIGN}}/epic_craterlake/SINGLE/{{PARTICLE}}/{ENERGY}/{PHASE_SPACE}/{{PARTICLE}}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.eicrecon.tree.edm4eic.root", + DETECTOR_CONFIG="epic_craterlake", + ENERGY=f"{float(wildcards.MOMENTUM):.0f}GeV" if float(wildcards.MOMENTUM) >= 1 else f"{float(wildcards.MOMENTUM) * 1000:.0f}MeV", + PHASE_SPACE=["3to50deg"], + INDEX=range(1), + )), + output: + "{CAMPAIGN}/{PARTICLE}/mom/lfhcal_mom_{MOMENTUM}_mom_resol_{PARTICLE}.root", + combined_root=temp("{CAMPAIGN}/lfhcal_sim_{MOMENTUM}_{PARTICLE}.root"), + shell: + """ +hadd {output.combined_root} {input.sim} +cd {wildcards.CAMPAIGN} +root -l -b -q ../{input.script}'("../{output.combined_root}", "{wildcards.PARTICLE}", {wildcards.MOMENTUM}, 0.15)' +""" + +rule lfhcal_summary_at_eta: + input: + expand( + [ + "{{CAMPAIGN}}/{{PARTICLE}}/mom/lfhcal_mom_{MOMENTUM:.1f}_mom_resol_{{PARTICLE}}.root", + ], + MOMENTUM=[0.5, 1.0, 2.0, 5.0, 10.0, 20.0], + ), + script="benchmarks/lfhcal/doCompare_widebins_mom.C", + output: + "{CAMPAIGN}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{PARTICLE}_{ETA_MIN}_eta_{ETA_MAX}.png", + "{CAMPAIGN}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{PARTICLE}_{ETA_MIN}_eta_{ETA_MAX}.root", + shell: + """ +if [[ "{wildcards.CAMPAIGN}" == "local" ]]; then + set +e + EPIC_VERSION="${{DETECTOR_VERSION:-}}" + EICRECON_VERSION="$(eicrecon -v | sed -n -e 's/.*\(v[0-9\.]\+\).*/\\1/p')" + # Legacy detection + : ${{EPIC_VERSION:="$(echo $DETECTOR_PATH | sed -n -e 's/.*epic-\([^-/]\+\).*/\\1/p')"}} + set -e + + echo "ePIC version: $EPIC_VERSION" + echo "EICrecon version: $EICRECON_VERSION" + EXTRA_LEGEND="ePIC $EPIC_VERSION / EICrecon $EICRECON_VERSION" +else + EXTRA_LEGEND="ePIC Simulation {wildcards.CAMPAIGN}" +fi +cd {wildcards.CAMPAIGN} +root -l -b -q ../{input.script}'("{wildcards.PARTICLE}", {wildcards.ETA_MIN}, {wildcards.ETA_MAX}, 1., true, "'"$EXTRA_LEGEND"'")' +""" + +LFHCAL_ETA_BINS = [1.2,1.5,2,2.5,3,3.5] + +rule lfhcal: + input: + lambda wildcards: expand( + [ + "{{CAMPAIGN}}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{PARTICLE}_{ETA_BIN}.png", + "{{CAMPAIGN}}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{PARTICLE}_{ETA_BIN}.root", + ], + ETA_BIN=[f"{eta_min:.1f}_eta_{eta_max:.1f}" for eta_min, eta_max in zip(LFHCAL_ETA_BINS[:-1], LFHCAL_ETA_BINS[1:])], + PARTICLE=["neutron", "pi-", "gamma"] if wildcards.CAMPAIGN == "local" else ["pi-"], + ) + output: + directory("results/lfhcal/{CAMPAIGN}/") + shell: + """ +mkdir {output} +cp {input} {output} +""" + + +rule lfhcal_local: + input: + "results/lfhcal/local", + + +rule lfhcal_campaigns: + input: + expand( + "results/lfhcal/{CAMPAIGN}", + CAMPAIGN=[ + "23.10.0", + "23.11.0", + "23.12.0", + "24.03.1", + "24.04.0", + ], + ) diff --git a/benchmarks/lfhcal/config.yml b/benchmarks/lfhcal/config.yml new file mode 100644 index 00000000..f7b9ad3a --- /dev/null +++ b/benchmarks/lfhcal/config.yml @@ -0,0 +1,48 @@ +sim:lfhcal: + extends: .det_benchmark + stage: simulate + parallel: + matrix: + - PARTICLE: ["neutron", "pi-"] + MOMENTUM: ["500MeV", "1GeV", "2GeV", "5GeV", "10GeV", "20GeV"] + script: + - | + snakemake --cache --cores 5 \ + sim_output/lfhcal/epic_craterlake/${PARTICLE}/${MOMENTUM}/3to50deg/${PARTICLE}_${MOMENTUM}_3to50deg.0000.eicrecon.tree.edm4eic.root + +bench:lfhcal: + extends: .det_benchmark + stage: benchmarks + needs: + - ["sim:lfhcal"] + script: + - snakemake --cores 3 lfhcal_local + +collect_results:lfhcal: + extends: .det_benchmark + stage: collect + needs: + - "bench:lfhcal" + script: + - ls -lrht + - mv results{,_save}/ # move results directory out of the way to preserve it + - snakemake --cores 1 --delete-all-output lfhcal_local + - mv results{_save,}/ + +bench:lfhcal_campaigns: + extends: .det_benchmark + stage: benchmarks + #when: manual + script: + - snakemake --cores 1 lfhcal_campaigns + +collect_results:lfhcal_campaigns: + extends: .det_benchmark + stage: collect + needs: + - "bench:lfhcal_campaigns" + script: + - ls -lrht + - mv results{,_save}/ # move results directory out of the way to preserve it + - snakemake --cores 1 --delete-all-output lfhcal_campaigns + - mv results{_save,}/ diff --git a/benchmarks/lfhcal/doCompare_widebins_mom.C b/benchmarks/lfhcal/doCompare_widebins_mom.C new file mode 100644 index 00000000..5258e5ef --- /dev/null +++ b/benchmarks/lfhcal/doCompare_widebins_mom.C @@ -0,0 +1,176 @@ +// Code to compare the tracking performances: Truth seeding vs real seeding +// Shyam Kumar; shyam.kumar@ba.infn.it; shyam055119@gmail.com + +#include "TGraphErrors.h" +#include "TF1.h" +#include "TRandom.h" +#include "TCanvas.h" +#include "TLegend.h" +#include "TMath.h" +#define mpi 0.139 // 1.864 GeV/c^2 + +//void draw_req_Mom(double etamin, double etamax, double xmin=0., double xmax=0.); +void doCompare_widebins_mom(TString particle = "pi-",double etamin=-1.0, double etamax=1.0, double range =0.3, Bool_t drawreq=1, TString extra_legend = "") // name = p, pt for getting p or pt dependence fitted results +{ + + //=== style of the plot========= + gStyle->SetPalette(1); + gStyle->SetOptTitle(1); + gStyle->SetTitleOffset(1.0,"XY"); + gStyle->SetTitleSize(.04,"XY"); + gStyle->SetLabelSize(.04,"XY"); + gStyle->SetHistLineWidth(2); + gStyle->SetOptFit(1); + gStyle->SetOptStat(1); + + const Int_t nfiles = 6; + double mom[nfiles] ={0.5,1.0,2.0,5.0,10.0,20.0}; + std::vector momV, momresolV, err_momresolV; + momV.clear(); momresolV.clear(); err_momresolV.clear(); + TString symbolname = ""; + if (particle == "pi-") symbolname = "#pi^{-}"; + else symbolname = particle; + ofstream outfile; + outfile.open ("Mom_resol.txt",ios_base::app); + + TF1 *f1=new TF1("f1","FitMomentumResolution",0.,30.0,2); + f1->SetParLimits(0,0.,0.1); + f1->SetParLimits(1,0.,5.0); + + TCanvas *c_mom = new TCanvas("cmom","cmom",1400,1000); + c_mom->SetMargin(0.10, 0.05 ,0.1,0.05); + c_mom->SetGridy(); + + //Reading the root file + TFile *fmom[nfiles]; + TGraphErrors *gr_mom; + TMultiGraph *mgMom; + TLegend *lmom; + mgMom = new TMultiGraph("mgMom",";p (GeV/c); #sigmap/p %"); + + lmom = new TLegend(0.65,0.80,0.90,0.93); + lmom->SetTextSize(0.03); + lmom->SetBorderSize(0); + lmom->SetHeader(extra_legend.Data(), "C"); + lmom->AddEntry((TObject*)0, Form("%s, %1.1f < #eta < %1.1f", symbolname.Data(), etamin, etamax), "C"); + + TF1 *func = new TF1("func","gaus",-0.5,0.5); + + for (int i =0; iSetMargin(0.10, 0.05 ,0.1,0.07); + + //pi-/mom/lfhcal_mom_20.0_mom_resol_pi-.root + fmom[i] = TFile::Open(Form("./%s/mom/lfhcal_mom_%1.1f_mom_resol_%s.root",particle.Data(),mom[i],particle.Data())); + + TH1D *hist = (TH1D*) fmom[i]->Get(Form("hist_mom_%1.1f_%1.1f_pmax_%1.1f",mom[i],etamin,etamax)); + hist->Rebin(2); + hist->SetName(Form("hist_mom_%1.1f_%1.1f_eta_%1.1f_%s",mom[i],etamin,etamax,particle.Data())); + hist->SetTitle(Form("Momentum = %1.1f && %1.1f<#eta<%1.1f;#Delta p/p; Entries(a.u.)",particle.Data(),mom[i],etamin,etamax)); + + double mu = hist->GetMean(); + double sigma = hist->GetRMS(); + double sigma_err = hist->GetRMSError(); + /* + hist->GetXaxis()->SetRangeUser(-1.0*range,1.0*range); + func->SetRange(mu-2.0*sigma,mu+2.0*sigma); // fit with in 2 sigma range + hist->Fit(func,"NR+"); + mu = func->GetParameter(1); + sigma = func->GetParameter(2); + func->SetRange(mu-2.0*sigma,mu+2.0*sigma); + hist->Fit(func,"R+"); + float par2 = func->GetParameter(2)*100; + float par2_err = func->GetParError(2)*100; + */ + momV.push_back(mom[i]); + momresolV.push_back(sigma); + err_momresolV.push_back(sigma_err); + + cp->cd(); + hist->Draw(); + //cp->SaveAs(Form("Debug_Plots/%s/mom/mom_resol_mom%1.1f_%1.1f_eta_%1.1f.png",particle.Data(),mom[i],etamin,etamax)); + } // all files + + const int size = momV.size(); + double p[size], err_p[size], sigma_p[size], err_sigma_p[size]; + + for (int i=0; iSetName("grseed"); + gr1->SetMarkerStyle(25); + gr1->SetMarkerColor(kBlue); + gr1->SetMarkerSize(2.0); + gr1->SetTitle(";p (GeV/c);#sigmap/p"); + gr1->GetXaxis()->CenterTitle(); + gr1->GetYaxis()->CenterTitle(); + + + mgMom->Add(gr1); + c_mom->cd(); + mgMom->GetXaxis()->SetRangeUser(0.40,20.2); + mgMom->GetYaxis()->SetRangeUser(0.0,1.50*TMath::MaxElement(gr1->GetN(),gr1->GetY())); // 50% more of the maximum value on yaxis + mgMom->Draw("AP"); + lmom->AddEntry(gr1,"Nominal"); + lmom->Draw("same"); + //draw_req_Mom(etamin,etamax,0.,mgMom->GetXaxis()->GetXmax()); + c_mom->SaveAs(Form("Final_Results/%s/mom/lfhcal_mom_resol_%s_%1.1f_eta_%1.1f.png",particle.Data(),particle.Data(),etamin,etamax)); + + // Write the numbers in output file for comparisons + outfile << extra_legend << endl; + outfile<<"Etamin"<GetN(); ++i){ + double x,y; + gr1->GetPoint(i,x,y); + outfile<cd(); + mgMom->SetName(Form("mom_resol_%1.1f_eta_%1.1f",etamin,etamax)); + mgMom->Write(); + fout->Close(); +} + +//===Fit Momentum Resolution +float FitMomentumResolution(Double_t *x, Double_t *par) +{ + float func = sqrt(par[0]*par[0]*x[0]*x[0]+par[1]*par[1]); + return func; +} + +//From Yellow report from section 11.2.2 + +// 1.2,1.5,2,2.5,3,3.5 + +/* +void draw_req_Mom(double etamin, double etamax, double xmin=0., double xmax=0.) +{ + + TF1 *dd4hep_p; + if (etamin >= -3.5 && etamax <= -2.5) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.1*x)^2+2.0^2)",xmin,xmax); + else if (etamin >= -2.5 && etamax <= -1.0) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.05*x)^2+1.0^2)",xmin,xmax); + else if (etamin >= -1.0 && etamax <= 1.0) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.05*x)^2+0.5^2)",xmin,xmax); + else if (etamin >= 1.0 && etamax <= 2.5) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.05*x)^2+1.0^2)",xmin,xmax); + else if (etamin >= 2.5 && etamax <= 3.5) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.1*x)^2+2.0^2)",xmin,xmax); + else return; + dd4hep_p->SetLineStyle(7); + dd4hep_p->SetLineColor(kMagenta); + dd4hep_p->SetLineWidth(3.0); + dd4hep_p->Draw("same"); + + TLegend *l= new TLegend(0.70,0.75,0.90,0.80); + l->SetTextSize(0.03); + l->SetBorderSize(0); + l->AddEntry(dd4hep_p,"PWGReq","l"); + l->Draw("same"); + } +*/ From fee8b16aa55210e4b73ca4da81f69306362633c8 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Thu, 10 Oct 2024 22:04:13 -0400 Subject: [PATCH 10/24] zdc_sigma: use maxfev=10000 (#79) --- benchmarks/zdc_sigma/analysis/sigma_plots.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/benchmarks/zdc_sigma/analysis/sigma_plots.py b/benchmarks/zdc_sigma/analysis/sigma_plots.py index fc422ed6..7ded93da 100644 --- a/benchmarks/zdc_sigma/analysis/sigma_plots.py +++ b/benchmarks/zdc_sigma/analysis/sigma_plots.py @@ -236,7 +236,7 @@ def gauss(x, A,mu, sigma): fnc=gauss p0=[100, 0, 0.5] coeff, var_matrix = curve_fit(fnc, bc[slc], y[slc], p0=p0, - sigma=np.sqrt(y[slc])+(y[slc]==0)) + sigma=np.sqrt(y[slc])+(y[slc]==0), maxfev=10000) x=np.linspace(-1, 1) plt.plot(x, gauss(x, *coeff), color='tab:orange') plt.xlabel("$\\theta^{*\\rm recon}_{\\Sigma}-\\theta^{*\\rm truth}_{\\Sigma}$ [mrad]") @@ -259,7 +259,7 @@ def gauss(x, A,mu, sigma): #print(bc[slc],y[slc]) sigma=np.sqrt(y[slc])+(y[slc]==0) try: - coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,sigma=list(sigma)) + coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000) sigmas.append(np.abs(coeff[2])) dsigmas.append(np.sqrt(var_matrix[2][2])) xvals.append(p) @@ -307,7 +307,7 @@ def gauss(x, A,mu, sigma): fnc=gauss p0=[100, 0, 1] coeff, var_matrix = curve_fit(fnc, bc[slc], y[slc], p0=p0, - sigma=np.sqrt(y[slc])+(y[slc]==0)) + sigma=np.sqrt(y[slc])+(y[slc]==0), maxfev=10000) x=np.linspace(-5, 5) plt.plot(x, gauss(x, *coeff), color='tab:orange') plt.xlabel("$z^{*\\rm recon}_{\\rm vtx}-z^{*\\rm truth}_{\\rm vtx}$ [m]") @@ -331,7 +331,7 @@ def gauss(x, A,mu, sigma): #print(bc[slc],y[slc]) sigma=np.sqrt(y[slc])+(y[slc]==0) try: - coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,sigma=list(sigma)) + coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000) sigmas.append(abs(coeff[2])) dsigmas.append(np.sqrt(var_matrix[2][2])) xvals.append(p) @@ -373,7 +373,7 @@ def gauss(x, A,mu, sigma): fnc=gauss p0=[100, lambda_mass, 0.03] coeff, var_matrix = curve_fit(fnc, bc[slc], y[slc], p0=p0, - sigma=np.sqrt(y[slc])+(y[slc]==0)) + sigma=np.sqrt(y[slc])+(y[slc]==0), maxfev=10000) x=np.linspace(0.8, 1.3, 200) plt.plot(x, gauss(x, *coeff), color='tab:orange') print(coeff[2], np.sqrt(var_matrix[2][2])) @@ -396,7 +396,7 @@ def gauss(x, A,mu, sigma): p0=[100, lambda_mass, 0.03] try: coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0, - sigma=list(np.sqrt(y[slc])+(y[slc]==0))) + sigma=list(np.sqrt(y[slc])+(y[slc]==0)), maxfev=10000) x=np.linspace(0.8, 1.3, 200) sigmas.append(coeff[2]) dsigmas.append(np.sqrt(var_matrix[2][2])) @@ -437,7 +437,7 @@ def gauss(x, A,mu, sigma): fnc=gauss p0=[100, sigma_mass, 0.03] coeff, var_matrix = curve_fit(fnc, bc[slc], y[slc], p0=p0, - sigma=np.sqrt(y[slc])+(y[slc]==0)) + sigma=np.sqrt(y[slc])+(y[slc]==0), maxfev=10000) x=np.linspace(0.8, 1.3, 200) plt.plot(x, gauss(x, *coeff), color='tab:orange') print(coeff[2], np.sqrt(var_matrix[2][2])) @@ -460,7 +460,7 @@ def gauss(x, A,mu, sigma): p0=[100, sigma_mass, 0.03] try: coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0, - sigma=list(np.sqrt(y[slc])+(y[slc]==0))) + sigma=list(np.sqrt(y[slc])+(y[slc]==0)), maxfev=10000) sigmas.append(abs(coeff[2])) dsigmas.append(np.sqrt(var_matrix[2][2])) xvals.append(p) From 3675fd60f511ef67659476f77155d1d1bcbac02b Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Fri, 11 Oct 2024 18:51:42 -0400 Subject: [PATCH 11/24] lfhcal: simulate photons in sim step (#85) --- benchmarks/lfhcal/config.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/lfhcal/config.yml b/benchmarks/lfhcal/config.yml index f7b9ad3a..a6dfcbe4 100644 --- a/benchmarks/lfhcal/config.yml +++ b/benchmarks/lfhcal/config.yml @@ -3,7 +3,7 @@ sim:lfhcal: stage: simulate parallel: matrix: - - PARTICLE: ["neutron", "pi-"] + - PARTICLE: ["gamma", "neutron", "pi-"] MOMENTUM: ["500MeV", "1GeV", "2GeV", "5GeV", "10GeV", "20GeV"] script: - | From 56ad9ae3af7f7d13ccafb39ad35de346e35c1123 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sun, 13 Oct 2024 02:22:41 -0400 Subject: [PATCH 12/24] Delete benchmarks/femc_photon/analysis/femc_photon_plots.py~ --- .../analysis/femc_photon_plots.py~ | 171 ------------------ 1 file changed, 171 deletions(-) delete mode 100644 benchmarks/femc_photon/analysis/femc_photon_plots.py~ diff --git a/benchmarks/femc_photon/analysis/femc_photon_plots.py~ b/benchmarks/femc_photon/analysis/femc_photon_plots.py~ deleted file mode 100644 index 4227b4e3..00000000 --- a/benchmarks/femc_photon/analysis/femc_photon_plots.py~ +++ /dev/null @@ -1,171 +0,0 @@ -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'] Date: Thu, 10 Oct 2024 21:47:15 -0400 Subject: [PATCH 13/24] benchmarks/femc_*: scale event count with momentum --- benchmarks/femc_electron/Snakefile | 20 ++++++++++++++------ benchmarks/femc_photon/Snakefile | 20 ++++++++++++++------ benchmarks/femc_pi0/Snakefile | 20 ++++++++++++++------ 3 files changed, 42 insertions(+), 18 deletions(-) diff --git a/benchmarks/femc_electron/Snakefile b/benchmarks/femc_electron/Snakefile index 0e8dd084..c2a83698 100644 --- a/benchmarks/femc_electron/Snakefile +++ b/benchmarks/femc_electron/Snakefile @@ -1,8 +1,15 @@ +def get_n_events(wildcards): + energy = float(wildcards.P) + n_events = 1000 + n_events = int(n_events // ((energy / 20) ** 0.5)) + return n_events + + rule femc_electron_generate: input: script="benchmarks/femc_electron/analysis/gen_particles.cxx", params: - NEVENTS_GEN=1000, + N_EVENTS=get_n_events, th_max=28, th_min=2.0 output: @@ -10,23 +17,23 @@ rule femc_electron_generate: 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})' +root -l -b -q '{input.script}({params.N_EVENTS},"{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: + N_EVENTS=get_n_events, 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 \ + --numberOfEvents {params.N_EVENTS} \ --physicsList {params.PHYSICS_LIST} \ --inputFiles {input.GEN_FILE} \ --outputFile {output.SIM_FILE} @@ -37,10 +44,11 @@ rule femc_electron_recon: 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" + params: + N_EVENTS=get_n_events, 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 +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={params.N_EVENTS} """ rule femc_electron_analysis: diff --git a/benchmarks/femc_photon/Snakefile b/benchmarks/femc_photon/Snakefile index 40524365..cf976044 100644 --- a/benchmarks/femc_photon/Snakefile +++ b/benchmarks/femc_photon/Snakefile @@ -1,8 +1,15 @@ +def get_n_events(wildcards): + energy = float(wildcards.P) + n_events = 1000 + n_events = int(n_events // ((energy / 20) ** 0.5)) + return n_events + + rule femc_photon_generate: input: script="benchmarks/femc_photon/analysis/gen_particles.cxx", params: - NEVENTS_GEN=1000, + N_EVENTS=get_n_events, th_max=28, th_min=2.0 output: @@ -10,23 +17,23 @@ rule femc_photon_generate: 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})' +root -l -b -q '{input.script}({params.N_EVENTS},"{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: + N_EVENTS=get_n_events, 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 \ + --numberOfEvents {params.N_EVENTS} \ --physicsList {params.PHYSICS_LIST} \ --inputFiles {input.GEN_FILE} \ --outputFile {output.SIM_FILE} @@ -37,10 +44,11 @@ rule femc_photon_recon: 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" + params: + N_EVENTS=get_n_events, 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 +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={params.N_EVENTS} """ rule femc_photon_analysis: diff --git a/benchmarks/femc_pi0/Snakefile b/benchmarks/femc_pi0/Snakefile index 4936040c..49497f06 100644 --- a/benchmarks/femc_pi0/Snakefile +++ b/benchmarks/femc_pi0/Snakefile @@ -1,8 +1,15 @@ +def get_n_events(wildcards): + energy = float(wildcards.P) + n_events = 1000 + n_events = int(n_events // ((energy / 20) ** 0.5)) + return n_events + + rule femc_pi0_generate: input: script="benchmarks/femc_pi0/analysis/gen_particles.cxx", params: - NEVENTS_GEN=1000, + N_EVENTS=get_n_events, th_max=28, th_min=2.0 output: @@ -10,23 +17,23 @@ rule femc_pi0_generate: 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})' +root -l -b -q '{input.script}({params.N_EVENTS},"{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: + N_EVENTS=get_n_events, 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 \ + --numberOfEvents {params.N_EVENTS} \ --physicsList {params.PHYSICS_LIST} \ --inputFiles {input.GEN_FILE} \ --outputFile {output.SIM_FILE} @@ -37,10 +44,11 @@ rule femc_pi0_recon: 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" + params: + N_EVENTS=get_n_events, 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 +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={params.N_EVENTS} """ rule femc_pi0_analysis: From 8f14f3171ebe59f935943c59404d6f840aae550f Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Fri, 11 Oct 2024 03:37:06 -0400 Subject: [PATCH 14/24] benchmarks/femc_*: use spaces consistently --- benchmarks/femc_electron/Snakefile | 72 +++++++++++++++--------------- benchmarks/femc_photon/Snakefile | 72 +++++++++++++++--------------- benchmarks/femc_pi0/Snakefile | 72 +++++++++++++++--------------- 3 files changed, 108 insertions(+), 108 deletions(-) diff --git a/benchmarks/femc_electron/Snakefile b/benchmarks/femc_electron/Snakefile index c2a83698..13e8b737 100644 --- a/benchmarks/femc_electron/Snakefile +++ b/benchmarks/femc_electron/Snakefile @@ -6,30 +6,30 @@ def get_n_events(wildcards): rule femc_electron_generate: - input: - script="benchmarks/femc_electron/analysis/gen_particles.cxx", - params: - N_EVENTS=get_n_events, - th_max=28, - th_min=2.0 - output: - GEN_FILE="sim_output/femc_electron/e-_{P}GeV.hepmc" - shell: - """ + input: + script="benchmarks/femc_electron/analysis/gen_particles.cxx", + params: + N_EVENTS=get_n_events, + 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.N_EVENTS},"{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: - N_EVENTS=get_n_events, - PHYSICS_LIST="FTFP_BERT" - output: - SIM_FILE="sim_output/femc_electron/{DETECTOR_CONFIG}_sim_e-_{P}GeV.edm4hep.root" - shell: - """ + input: + GEN_FILE="sim_output/femc_electron/e-_{P}GeV.hepmc" + params: + N_EVENTS=get_n_events, + PHYSICS_LIST="FTFP_BERT" + output: + SIM_FILE="sim_output/femc_electron/{DETECTOR_CONFIG}_sim_e-_{P}GeV.edm4hep.root" + shell: + """ # Running simulation npsim \ --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \ @@ -40,27 +40,27 @@ npsim \ """ 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" - params: - N_EVENTS=get_n_events, - shell: - """ + 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" + params: + N_EVENTS=get_n_events, + shell: + """ 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={params.N_EVENTS} """ 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: - """ + 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_photon/Snakefile b/benchmarks/femc_photon/Snakefile index cf976044..d92a3cd1 100644 --- a/benchmarks/femc_photon/Snakefile +++ b/benchmarks/femc_photon/Snakefile @@ -6,30 +6,30 @@ def get_n_events(wildcards): rule femc_photon_generate: - input: - script="benchmarks/femc_photon/analysis/gen_particles.cxx", - params: - N_EVENTS=get_n_events, - th_max=28, - th_min=2.0 - output: - GEN_FILE="sim_output/femc_photon/photon_{P}GeV.hepmc" - shell: - """ + input: + script="benchmarks/femc_photon/analysis/gen_particles.cxx", + params: + N_EVENTS=get_n_events, + 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.N_EVENTS},"{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: - N_EVENTS=get_n_events, - PHYSICS_LIST="FTFP_BERT" - output: - SIM_FILE="sim_output/femc_photon/{DETECTOR_CONFIG}_sim_photon_{P}GeV.edm4hep.root" - shell: - """ + input: + GEN_FILE="sim_output/femc_photon/photon_{P}GeV.hepmc" + params: + N_EVENTS=get_n_events, + PHYSICS_LIST="FTFP_BERT" + output: + SIM_FILE="sim_output/femc_photon/{DETECTOR_CONFIG}_sim_photon_{P}GeV.edm4hep.root" + shell: + """ # Running simulation npsim \ --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \ @@ -40,27 +40,27 @@ npsim \ """ 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" - params: - N_EVENTS=get_n_events, - shell: - """ + 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" + params: + N_EVENTS=get_n_events, + shell: + """ 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={params.N_EVENTS} """ 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: - """ + 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_pi0/Snakefile b/benchmarks/femc_pi0/Snakefile index 49497f06..9b8713ab 100644 --- a/benchmarks/femc_pi0/Snakefile +++ b/benchmarks/femc_pi0/Snakefile @@ -6,30 +6,30 @@ def get_n_events(wildcards): rule femc_pi0_generate: - input: - script="benchmarks/femc_pi0/analysis/gen_particles.cxx", - params: - N_EVENTS=get_n_events, - th_max=28, - th_min=2.0 - output: - GEN_FILE="sim_output/femc_pi0/pi0_{P}GeV.hepmc" - shell: - """ + input: + script="benchmarks/femc_pi0/analysis/gen_particles.cxx", + params: + N_EVENTS=get_n_events, + 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.N_EVENTS},"{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: - N_EVENTS=get_n_events, - PHYSICS_LIST="FTFP_BERT" - output: - SIM_FILE="sim_output/femc_pi0/{DETECTOR_CONFIG}_sim_pi0_{P}GeV.edm4hep.root" - shell: - """ + input: + GEN_FILE="sim_output/femc_pi0/pi0_{P}GeV.hepmc" + params: + N_EVENTS=get_n_events, + PHYSICS_LIST="FTFP_BERT" + output: + SIM_FILE="sim_output/femc_pi0/{DETECTOR_CONFIG}_sim_pi0_{P}GeV.edm4hep.root" + shell: + """ # Running simulation npsim \ --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \ @@ -40,27 +40,27 @@ npsim \ """ 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" - params: - N_EVENTS=get_n_events, - shell: - """ + 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" + params: + N_EVENTS=get_n_events, + shell: + """ 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={params.N_EVENTS} """ 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: - """ + 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} """ From 0e836a1cab9906c3fd1753ea75c2d69bdbb2e8f2 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Tue, 8 Oct 2024 19:01:50 -0400 Subject: [PATCH 15/24] backwards_ecal: add energy resolution benchmark --- benchmarks/backwards_ecal/backwards_ecal.org | 116 +++++++++++++++++++ 1 file changed, 116 insertions(+) diff --git a/benchmarks/backwards_ecal/backwards_ecal.org b/benchmarks/backwards_ecal/backwards_ecal.org index d3ffca5e..e2acf583 100644 --- a/benchmarks/backwards_ecal/backwards_ecal.org +++ b/benchmarks/backwards_ecal/backwards_ecal.org @@ -113,6 +113,122 @@ for energy in energies: )) #+end_src +** Energy resolution + +#+begin_src jupyter-python +fig, axs = plt.subplots(2, 4, sharex=True, sharey=True, figsize=(15, 6)) + +fig.suptitle(PLOT_TITLE) + +axs = np.ravel(np.array(axs)) + +sigmas_rel_FWHM_cb = {} +fractions_below = {} + +for ix, energy in enumerate(energies): + for use_clusters in [False, True]: + energy_value = float(energy.replace("GeV", "").replace("MeV", "e-3")) + if use_clusters: + clf_label = "leading cluster" + else: + clf_label = "sum all hits" + def clf(events): + if use_clusters: + return ak.drop_none(ak.max(events["EcalEndcapNClusters.energy"], axis=-1)) / energy_value + else: + return ak.sum(events["EcalEndcapNRecHits.energy"], axis=-1) / energy_value + e_pred = clf(e_eval[energy]) + + plt.sca(axs[ix]) + counts, bins, patches = plt.hist(e_pred, weights=np.full_like(e_pred, 1.0 / ak.num(e_pred, axis=0)), bins=np.linspace(0.01, 1.01, 101), label=rf"$e^-$ {clf_label}", hatch=None if use_clusters else r"xxx", alpha=0.8 if use_clusters else 1.) + plt.title(f"{energy}") + + e_over_p = (bins[1:] + bins[:-1]) / 2 + import scipy.stats + def f(x, n, beta, m, loc, scale): + return n * scipy.stats.crystalball.pdf(x, beta, m, loc, scale) + p0 = (np.sum(counts[10:]), 2., 3., 0.95, 0.05) + + try: + import scipy.optimize + par, pcov = scipy.optimize.curve_fit(f, e_over_p[5:], counts[5:], p0=p0, maxfev=10000) + except RuntimeError: + par = None + plt.plot(e_over_p, f(e_over_p, *par), label=rf"Crystal Ball fit", color="tab:green" if use_clusters else "green", lw=0.8) + + def summarize_fit(par): + _, _, _, loc_cb, scale_cb = par + # Calculate FWHM + y_max = np.max(f(np.linspace(0., 1., 100), *par)) + f_prime = lambda x: f(x, *par) - y_max / 2 + x_plus, = scipy.optimize.root(f_prime, loc_cb + scale_cb).x + x_minus, = scipy.optimize.root(f_prime, loc_cb - scale_cb).x + color = "cyan" if use_clusters else "orange" + plt.axvline(x_minus, ls="--", lw=0.75, color=patches[0].get_facecolor(), label=r"$\mu - $FWHM") + plt.axvline(x_plus, ls=":", lw=0.75, color=patches[0].get_facecolor(), label=r"$\mu + $FWHM") + fwhm = (x_plus - x_minus) / loc_cb + sigma_rel_FWHM_cb = fwhm / 2 / np.sqrt(2 * np.log(2)) + + cutoff_x = loc_cb - fwhm + fraction_below = np.sum(counts[e_over_p < cutoff_x]) / ak.num(e_pred, axis=0) + + return sigma_rel_FWHM_cb, fraction_below + + sigma_rel_FWHM_cb, fraction_below = summarize_fit(par) + sigmas_rel_FWHM_cb.setdefault(clf_label, {})[energy] = sigma_rel_FWHM_cb + fractions_below.setdefault(clf_label, {})[energy] = fraction_below + + plt.legend() + plt.xlabel("$E/p$", loc="right") + plt.ylabel("Event yield", loc="top") + +fig.savefig(output_dir / f"resolution_plots.pdf", bbox_inches="tight") +plt.show() +plt.close(fig) + +plt.figure() +energy_values = np.array([float(energy.replace("GeV", "").replace("MeV", "e-3")) for energy in energies]) + +for clf_label, sigma_rel_FWHM_cb in sigmas_rel_FWHM_cb.items(): + sigma_over_e = np.array([sigma_rel_FWHM_cb[energy] for energy in energies]) * 100 # convert to % + + def f(energy, stochastic, constant): + return np.sqrt((stochastic / np.sqrt(energy)) ** 2 + constant ** 2) + cond = energy_values >= 0.5 + try: + import scipy.optimize + par, pcov = scipy.optimize.curve_fit(f, energy_values[cond], sigma_over_e[cond], maxfev=10000) + except RuntimeError: + par = None + stochastic, constant = par + + plt.plot( + energy_values, + sigma_over_e, + marker=".", + label=f"{clf_label}" + ) + plt.plot( + energy_values[cond], + f(energy_values[cond], *par), + color="black", + ls="--", + lw=0.5, + label=f"{clf_label}, ${np.ceil(stochastic * 10) / 10:.1f}\% / \sqrt{{E}} \oplus {np.ceil(constant * 10) / 10:.1f}\%$", + ) +plt.plot( + energy_values, + np.sqrt((1 / energy_values) ** 2 + (1 / np.sqrt(energy_values)) ** 2 + 1 ** 2), + color="black", label=r"YR requirement $1\% / E \oplus 2.5\% / \sqrt{E} \oplus 1\%$", +) +plt.title(INPUT_PATH_FORMAT) +plt.legend() +plt.xlabel("Energy, GeV", loc="right") +plt.ylabel(r"$\sigma_{E} / E$ derived from FWHM, %", loc="top") +plt.savefig(output_dir / f"resolution.pdf", bbox_inches="tight") +plt.show() +#+end_src + ** Pion rejection #+begin_src jupyter-python From dfd7510f6afdbc44b1383031aee50640c16636c5 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Tue, 8 Oct 2024 20:15:35 -0400 Subject: [PATCH 16/24] backwards_ecal: output png --- benchmarks/backwards_ecal/backwards_ecal.org | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/benchmarks/backwards_ecal/backwards_ecal.org b/benchmarks/backwards_ecal/backwards_ecal.org index e2acf583..27f1094c 100644 --- a/benchmarks/backwards_ecal/backwards_ecal.org +++ b/benchmarks/backwards_ecal/backwards_ecal.org @@ -183,6 +183,7 @@ for ix, energy in enumerate(energies): plt.ylabel("Event yield", loc="top") fig.savefig(output_dir / f"resolution_plots.pdf", bbox_inches="tight") +fig.savefig(output_dir / f"resolution_plots.png", bbox_inches="tight") plt.show() plt.close(fig) @@ -226,6 +227,7 @@ plt.legend() plt.xlabel("Energy, GeV", loc="right") plt.ylabel(r"$\sigma_{E} / E$ derived from FWHM, %", loc="top") plt.savefig(output_dir / f"resolution.pdf", bbox_inches="tight") +plt.savefig(output_dir / f"resolution.png", bbox_inches="tight") plt.show() #+end_src @@ -292,10 +294,13 @@ for ix, energy in enumerate(energies): plt.ylabel("Pion rejection factor") fig.savefig(output_dir / f"pred.pdf", bbox_inches="tight") +fig.savefig(output_dir / f"pred.png", bbox_inches="tight") plt.close(fig) fig_log.savefig(output_dir / f"pred_log.pdf", bbox_inches="tight") +fig_log.savefig(output_dir / f"pred_log.png", bbox_inches="tight") fig_log.show() fig_roc.savefig(output_dir / f"roc.pdf", bbox_inches="tight") +fig_roc.savefig(output_dir / f"roc.png", bbox_inches="tight") fig_roc.show() plt.figure() @@ -312,5 +317,6 @@ plt.legend() plt.xlabel("Energy, GeV") plt.ylabel("Pion rejection at 95%") plt.savefig(output_dir / f"pion_rej.pdf", bbox_inches="tight") +plt.savefig(output_dir / f"pion_rej.png", bbox_inches="tight") plt.show() #+end_src From 939acc4559f55e8d3dfd2db18271bc4684e976e9 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sun, 13 Oct 2024 02:18:27 -0400 Subject: [PATCH 17/24] Add a basic .pre-commit-config.yaml --- .pre-commit-config.yaml | 7 +++++++ benchmarks/zdc_photon/config.yml | 1 - 2 files changed, 7 insertions(+), 1 deletion(-) create mode 100644 .pre-commit-config.yaml diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 00000000..15f050fc --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,7 @@ +repos: +- repo: https://github.com/pre-commit/pre-commit-hooks + rev: v5.0.0 + hooks: + - id: check-added-large-files + args: ['--maxkb=128'] + - id: check-yaml diff --git a/benchmarks/zdc_photon/config.yml b/benchmarks/zdc_photon/config.yml index aee4f6c6..d62cef61 100644 --- a/benchmarks/zdc_photon/config.yml +++ b/benchmarks/zdc_photon/config.yml @@ -29,7 +29,6 @@ bench:zdc_photon: collect_results:zdc_photon: extends: .det_benchmark stage: collect - extends: .det_benchmark needs: ["bench:zdc_photon"] script: - ls -al From 2a39b59cea82b8884cbd8dabd9530e6f2ad7b0cf Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Tue, 15 Oct 2024 12:38:31 -0400 Subject: [PATCH 18/24] backwards_ecal: add simulation caching (#83) --- benchmarks/backwards_ecal/Snakefile | 30 ++++++++++++++++++++++++++--- 1 file changed, 27 insertions(+), 3 deletions(-) diff --git a/benchmarks/backwards_ecal/Snakefile b/benchmarks/backwards_ecal/Snakefile index 003d0ab7..da90f2e2 100644 --- a/benchmarks/backwards_ecal/Snakefile +++ b/benchmarks/backwards_ecal/Snakefile @@ -4,11 +4,29 @@ def get_n_events(wildcards): n_events = int(n_events // (energy ** 0.5)) return n_events +import functools +import json +import ctypes.util +import warnings +from snakemake.logging import logger + +@functools.cache +def get_spack_package_hash(package_name): + try: + ver_info = json.loads(subprocess.check_output(["spack", "find", "--json", package_name])) + return ver_info[0]["package_hash"] + except FileNotFoundError as e: + logger.warning("Spack is not installed") + return "" + except subprocess.CalledProcessError as e: + print(e) + return "" rule backwards_ecal_sim: input: steering_file=ancient("EPIC/EVGEN/SINGLE/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.steer"), warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root", + geometry_lib=os.environ["DETECTOR_PATH"] + "/../../lib/" + ctypes.util.find_library("epic"), output: "sim_output/backwards_ecal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.edm4hep.root", log: @@ -19,7 +37,13 @@ rule backwards_ecal_sim: PHASE_SPACE="(3to50|45to135|130to177)deg", INDEX="\d{4}", params: - N_EVENTS=get_n_events + N_EVENTS=get_n_events, + SEED=lambda wildcards: "1" + wildcards.INDEX, + DETECTOR_PATH=os.environ["DETECTOR_PATH"], + DETECTOR_CONFIG=lambda wildcards: wildcards.DETECTOR_CONFIG, + DD4HEP_HASH=get_spack_package_hash("dd4hep"), + NPSIM_HASH=get_spack_package_hash("npsim"), + cache: True shell: """ set -m # monitor mode to prevent lingering processes @@ -27,11 +51,11 @@ exec ddsim \ --runType batch \ --enableGun \ --steeringFile "{input.steering_file}" \ - --random.seed 1{wildcards.INDEX} \ + --random.seed {params.SEED} \ --filter.tracker edep0 \ -v WARNING \ --numberOfEvents {params.N_EVENTS} \ - --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \ + --compactFile {params.DETECTOR_PATH}/{params.DETECTOR_CONFIG}.xml \ --outputFile {output} """ From 6e99ae8b291211ef519a21eb414cf833fc7e4c31 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 19 Oct 2024 09:24:52 -0400 Subject: [PATCH 19/24] .gitlab-ci.yml: introduce global SNAKEMAKE_FLAGS, don't re-run based on mtime (#92) --- .gitlab-ci.yml | 1 + benchmarks/backgrounds/config.yml | 6 +++--- benchmarks/backwards_ecal/config.yml | 6 +++--- benchmarks/barrel_ecal/config.yml | 20 +++++++++---------- benchmarks/ecal_gaps/config.yml | 6 +++--- benchmarks/femc_electron/config.yml | 6 +++--- benchmarks/femc_photon/config.yml | 6 +++--- benchmarks/femc_pi0/config.yml | 6 +++--- benchmarks/insert_muon/config.yml | 6 +++--- benchmarks/insert_neutron/config.yml | 6 +++--- benchmarks/lfhcal/config.yml | 8 ++++---- benchmarks/material_scan/config.yml | 2 +- benchmarks/tracking_performances/config.yml | 8 ++++---- .../tracking_performances_dis/config.yml | 6 +++--- benchmarks/zdc_lambda/config.yml | 6 +++--- benchmarks/zdc_lyso/config.yml | 4 ++-- benchmarks/zdc_photon/config.yml | 2 +- benchmarks/zdc_pi0/config.yml | 2 +- benchmarks/zdc_sigma/config.yml | 6 +++--- 19 files changed, 57 insertions(+), 56 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index f25c74e9..7a19fd1a 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -5,6 +5,7 @@ variables: DETECTOR_CONFIG: epic_craterlake GITHUB_SHA: '' GITHUB_REPOSITORY: '' + SNAKEMAKE_FLAGS: '--cache --rerun-triggers code input params software-env' workflow: name: '$PIPELINE_NAME' diff --git a/benchmarks/backgrounds/config.yml b/benchmarks/backgrounds/config.yml index c4d084d7..79b864bb 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 $SNAKEMAKE_FLAGS --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 @@ -19,7 +19,7 @@ bench:backgrounds_emcal_backwards: - ln -s $LOCAL_DATA_PATH/input input - export PYTHONUSERBASE=$LOCAL_DATA_PATH/deps - pip install -r benchmarks/backgrounds/requirements.txt - - snakemake --cores 8 backgrounds_ecal_backwards + - snakemake $SNAKEMAKE_FLAGS --cores 8 backgrounds_ecal_backwards collect_results:backgrounds: extends: .det_benchmark @@ -29,5 +29,5 @@ collect_results:backgrounds: script: - ls -lrht - mv results{,_save}/ # move results directory out of the way to preserve it - - snakemake --cores 1 --delete-all-output backgrounds_ecal_backwards + - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output backgrounds_ecal_backwards - mv results{_save,}/ diff --git a/benchmarks/backwards_ecal/config.yml b/benchmarks/backwards_ecal/config.yml index 9eff9b58..7048f8e7 100644 --- a/benchmarks/backwards_ecal/config.yml +++ b/benchmarks/backwards_ecal/config.yml @@ -16,7 +16,7 @@ sim:backwards_ecal: ] script: - | - snakemake --cache --cores 5 sim_output/backwards_ecal/${DETECTOR_CONFIG}/${PARTICLE}/${MOMENTUM}/130to177deg/flag + snakemake $SNAKEMAKE_FLAGS --cores 5 sim_output/backwards_ecal/${DETECTOR_CONFIG}/${PARTICLE}/${MOMENTUM}/130to177deg/flag bench:backwards_ecal: extends: .det_benchmark @@ -26,7 +26,7 @@ bench:backwards_ecal: script: - export PYTHONUSERBASE=$LOCAL_DATA_PATH/deps - pip install -r benchmarks/backwards_ecal/requirements.txt - - snakemake --cores 1 backwards_ecal + - snakemake $SNAKEMAKE_FLAGS --cores 1 backwards_ecal collect_results:backwards_ecal: extends: .det_benchmark @@ -36,5 +36,5 @@ collect_results:backwards_ecal: script: - ls -lrht - mv results{,_save}/ # move results directory out of the way to preserve it - - snakemake --cores 1 --delete-all-output backwards_ecal + - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output backwards_ecal - mv results{_save,}/ diff --git a/benchmarks/barrel_ecal/config.yml b/benchmarks/barrel_ecal/config.yml index a236f68b..45fc3e1c 100644 --- a/benchmarks/barrel_ecal/config.yml +++ b/benchmarks/barrel_ecal/config.yml @@ -2,13 +2,13 @@ sim:emcal_barrel_pions: extends: .det_benchmark stage: simulate script: - - snakemake --cores 1 $DETECTOR_CONFIG/sim_output/sim_emcal_barrel_{piplus,piminus}_energies5.0_5.0.edm4hep.root + - snakemake $SNAKEMAKE_FLAGS --cores 1 $DETECTOR_CONFIG/sim_output/sim_emcal_barrel_{piplus,piminus}_energies5.0_5.0.edm4hep.root sim:emcal_barrel_pi0: extends: .det_benchmark stage: simulate script: - - snakemake --cores 1 $DETECTOR_CONFIG/sim_output/sim_emcal_barrel_pi0_energies5.0_5.0.edm4hep.root + - snakemake $SNAKEMAKE_FLAGS --cores 1 $DETECTOR_CONFIG/sim_output/sim_emcal_barrel_pi0_energies5.0_5.0.edm4hep.root sim:emcal_barrel_electrons: extends: .det_benchmark @@ -16,20 +16,20 @@ sim:emcal_barrel_electrons: script: - if [[ "$RUN_EXTENDED_BENCHMARK" == "true" ]] ; then snakemake --cores $DETECTOR_CONFIG/results/energy_scan/emcal_barrel_electron_fsam_scan.png; fi - export JUGGLER_N_EVENTS=400 - - snakemake --cores 1 $DETECTOR_CONFIG/sim_output/sim_emcal_barrel_electron_energies5.0_5.0.edm4hep.root + - snakemake $SNAKEMAKE_FLAGS --cores 1 $DETECTOR_CONFIG/sim_output/sim_emcal_barrel_electron_energies5.0_5.0.edm4hep.root sim:emcal_barrel_photons: extends: .det_benchmark stage: simulate script: - if [[ "$RUN_EXTENDED_BENCHMARK" == "true" ]] ; then snakemake --cores $DETECTOR_CONFIG/results/energy_scan/emcal_barrel_proton_fsam_scan.png; fi - - snakemake --cores 1 $DETECTOR_CONFIG/sim_output/sim_emcal_barrel_photon_energies5.0_5.0.edm4hep.root + - snakemake $SNAKEMAKE_FLAGS --cores 1 $DETECTOR_CONFIG/sim_output/sim_emcal_barrel_photon_energies5.0_5.0.edm4hep.root sim:emcal_barrel_pion_rejection: extends: .det_benchmark stage: simulate script: - - snakemake --cores 1 $DETECTOR_CONFIG/sim_output/sim_emcal_barrel_{piminus,electron}_energies1.0_18.0.edm4hep.root + - snakemake $SNAKEMAKE_FLAGS --cores 1 $DETECTOR_CONFIG/sim_output/sim_emcal_barrel_{piminus,electron}_energies1.0_18.0.edm4hep.root calib:emcal_barrel_electrons: extends: .det_benchmark @@ -39,7 +39,7 @@ calib:emcal_barrel_electrons: script: - ls -lhtR sim_output/ - rootls -t sim_output/sim_emcal_barrel_electron_energies5.0_5.0.edm4hep.root - - snakemake --cores 1 $DETECTOR_CONFIG/results/emcal_barrel_electron_calibration.json + - snakemake $SNAKEMAKE_FLAGS --cores 1 $DETECTOR_CONFIG/results/emcal_barrel_electron_calibration.json - mv sim_output/sim_emcal_barrel_electron.edm4hep.root results/. - echo "JSON file(s) from analysis:" ; cat results/*.json @@ -49,7 +49,7 @@ bench:emcal_barrel_pions: needs: - ["sim:emcal_barrel_pions"] script: - - snakemake --cores 1 $DETECTOR_CONFIG/results/emcal_barrel_pions_Ethr.png + - snakemake $SNAKEMAKE_FLAGS --cores 1 $DETECTOR_CONFIG/results/emcal_barrel_pions_Ethr.png bench:emcal_barrel_electrons_scan: extends: .det_benchmark @@ -66,7 +66,7 @@ bench:emcal_barrel_pi0: - ["sim:emcal_barrel_pi0", "calib:emcal_barrel_electrons"] script: - echo "JSON file(s) from analysis:" ; cat results/*.json - - snakemake --cores 1 epic_craterlake/results/Barrel_emcal_pi0.json + - snakemake $SNAKEMAKE_FLAGS --cores 1 epic_craterlake/results/Barrel_emcal_pi0.json bench:emcal_barrel_photons: extends: .det_benchmark @@ -76,7 +76,7 @@ bench:emcal_barrel_photons: script: - ls -lhtR sim_output/ - rootls -t sim_output/sim_emcal_barrel_photon_energies5.0_5.0.edm4hep.root - - snakemake --cores 1 $DETECTOR_CONFIG/results/emcal_barrel_photon_calibration.json + - snakemake $SNAKEMAKE_FLAGS --cores 1 $DETECTOR_CONFIG/results/emcal_barrel_photon_calibration.json - mv sim_output/sim_emcal_barrel_photon.edm4hep.root results/. - if [[ "$RUN_EXTENDED_BENCHMARK" == "true" ]] ; then snakemake --cores $DETECTOR_CONFIG/results/energy_scan/emcal_barrel_proton_fsam_scan.png; fi @@ -89,7 +89,7 @@ bench:emcal_barrel_pion_rejection: - ls -lhtR sim_output/ - rootls -t $DETECTOR_CONFIG/sim_output/sim_emcal_barrel_piminus_energies1.0_18.0.edm4hep.root - rootls -t $DETECTOR_CONFIG/sim_output/sim_emcal_barrel_electron_energies1.0_18.0.edm4hep.root - - snakemake --cores 1 $DETECTOR_CONFIG/results/Barrel_emcal_pion_rej.json + - snakemake $SNAKEMAKE_FLAGS --cores 1 $DETECTOR_CONFIG/results/Barrel_emcal_pion_rej.json collect_results:barrel_ecal: extends: .det_benchmark diff --git a/benchmarks/ecal_gaps/config.yml b/benchmarks/ecal_gaps/config.yml index 9d967c87..9c0a8a09 100644 --- a/benchmarks/ecal_gaps/config.yml +++ b/benchmarks/ecal_gaps/config.yml @@ -4,7 +4,7 @@ sim:ecal_gaps: script: - mkdir -p $LOCAL_DATA_PATH/input - ln -s $LOCAL_DATA_PATH/input input - - snakemake --cache --cores 10 ecal_gaps --omit-from ecal_gaps + - snakemake $SNAKEMAKE_FLAGS --cache --cores 10 ecal_gaps --omit-from ecal_gaps bench:ecal_gaps: extends: .det_benchmark @@ -15,7 +15,7 @@ bench:ecal_gaps: - ln -s $LOCAL_DATA_PATH/input input - export PYTHONUSERBASE=$LOCAL_DATA_PATH/deps - pip install -r benchmarks/ecal_gaps/requirements.txt - - snakemake --cores 8 ecal_gaps + - snakemake $SNAKEMAKE_FLAGS --cores 8 ecal_gaps collect_results:ecal_gaps: extends: .det_benchmark @@ -25,5 +25,5 @@ collect_results:ecal_gaps: script: - ls -lrht - mv results{,_save}/ # move results directory out of the way to preserve it - - snakemake --cores 1 --delete-all-output ecal_gaps + - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output ecal_gaps - mv results{_save,}/ diff --git a/benchmarks/femc_electron/config.yml b/benchmarks/femc_electron/config.yml index 7b6f1043..11c8deea 100644 --- a/benchmarks/femc_electron/config.yml +++ b/benchmarks/femc_electron/config.yml @@ -13,7 +13,7 @@ sim:femc_electron: - P: 80 timeout: 1 hours script: - - snakemake --cores 1 sim_output/femc_electron/epic_craterlake_rec_e-_${P}GeV.edm4eic.root + - snakemake $SNAKEMAKE_FLAGS --cores 1 sim_output/femc_electron/epic_craterlake_rec_e-_${P}GeV.edm4eic.root retry: max: 2 when: @@ -24,7 +24,7 @@ bench:femc_electron: stage: benchmarks needs: ["sim:femc_electron"] script: - - snakemake --cores 1 results/epic_craterlake/femc_electron + - snakemake $SNAKEMAKE_FLAGS --cores 1 results/epic_craterlake/femc_electron collect_results:femc_electron: extends: .det_benchmark @@ -33,5 +33,5 @@ collect_results: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 + - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output results/epic_craterlake/femc_electron - mv results{_save,}/ diff --git a/benchmarks/femc_photon/config.yml b/benchmarks/femc_photon/config.yml index af710897..dca2694e 100644 --- a/benchmarks/femc_photon/config.yml +++ b/benchmarks/femc_photon/config.yml @@ -13,7 +13,7 @@ sim:femc_photon: - P: 80 timeout: 1 hours script: - - snakemake --cores 1 sim_output/femc_photon/epic_craterlake_rec_photon_${P}GeV.edm4eic.root + - snakemake $SNAKEMAKE_FLAGS --cores 1 sim_output/femc_photon/epic_craterlake_rec_photon_${P}GeV.edm4eic.root retry: max: 2 when: @@ -24,7 +24,7 @@ bench:femc_photon: stage: benchmarks needs: ["sim:femc_photon"] script: - - snakemake --cores 1 results/epic_craterlake/femc_photon + - snakemake $SNAKEMAKE_FLAGS --cores 1 results/epic_craterlake/femc_photon collect_results:femc_photon: extends: .det_benchmark @@ -33,5 +33,5 @@ collect_results: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 + - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output results/epic_craterlake/femc_photon - mv results{_save,}/ diff --git a/benchmarks/femc_pi0/config.yml b/benchmarks/femc_pi0/config.yml index 46e6ed9e..787cb583 100644 --- a/benchmarks/femc_pi0/config.yml +++ b/benchmarks/femc_pi0/config.yml @@ -12,7 +12,7 @@ sim:femc_pi0: - P: 70 - P: 80 script: - - snakemake --cores 1 sim_output/femc_pi0/epic_craterlake_rec_pi0_${P}GeV.edm4eic.root + - snakemake $SNAKEMAKE_FLAGS --cores 1 sim_output/femc_pi0/epic_craterlake_rec_pi0_${P}GeV.edm4eic.root retry: max: 2 when: @@ -23,7 +23,7 @@ bench:femc_pi0: stage: benchmarks needs: ["sim:femc_pi0"] script: - - snakemake --cores 1 results/epic_craterlake/femc_pi0 + - snakemake $SNAKEMAKE_FLAGS --cores 1 results/epic_craterlake/femc_pi0 collect_results:femc_pi0: extends: .det_benchmark @@ -32,5 +32,5 @@ collect_results: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 + - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output results/epic_craterlake/femc_pi0 - mv results{_save,}/ diff --git a/benchmarks/insert_muon/config.yml b/benchmarks/insert_muon/config.yml index 76899c1d..4e3272af 100644 --- a/benchmarks/insert_muon/config.yml +++ b/benchmarks/insert_muon/config.yml @@ -5,7 +5,7 @@ sim:insert_muon: matrix: - P: 50 script: - - snakemake --cores 5 sim_output/insert_muon/epic_craterlake_sim_mu-_${P}GeV_{0,1,2,3,4}.edm4hep.root + - snakemake $SNAKEMAKE_FLAGS --cores 5 sim_output/insert_muon/epic_craterlake_sim_mu-_${P}GeV_{0,1,2,3,4}.edm4hep.root retry: max: 2 when: @@ -16,7 +16,7 @@ bench:insert_muon: stage: benchmarks needs: ["sim:insert_muon"] script: - - snakemake --cores 1 results/epic_craterlake/insert_muon + - snakemake $SNAKEMAKE_FLAGS --cores 1 results/epic_craterlake/insert_muon collect_results:insert_muon: extends: .det_benchmark @@ -25,5 +25,5 @@ collect_results:insert_muon: 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/insert_muon + - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output results/epic_craterlake/insert_muon - mv results{_save,}/ diff --git a/benchmarks/insert_neutron/config.yml b/benchmarks/insert_neutron/config.yml index eedac937..90435261 100644 --- a/benchmarks/insert_neutron/config.yml +++ b/benchmarks/insert_neutron/config.yml @@ -11,7 +11,7 @@ sim:insert_neutron: - P: 70 - P: 80 script: - - snakemake --cores 5 sim_output/insert_neutron/epic_craterlake_rec_neutron_${P}GeV_{0,1,2,3,4}.edm4eic.root + - snakemake $SNAKEMAKE_FLAGS --cores 5 sim_output/insert_neutron/epic_craterlake_rec_neutron_${P}GeV_{0,1,2,3,4}.edm4eic.root retry: max: 2 when: @@ -22,7 +22,7 @@ bench:insert_neutron: stage: benchmarks needs: ["sim:insert_neutron"] script: - - snakemake --cores 1 results/epic_craterlake/insert_neutron + - snakemake $SNAKEMAKE_FLAGS --cores 1 results/epic_craterlake/insert_neutron collect_results:insert_neutron: extends: .det_benchmark @@ -31,5 +31,5 @@ collect_results:insert_neutron: 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/insert_neutron + - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output results/epic_craterlake/insert_neutron - mv results{_save,}/ diff --git a/benchmarks/lfhcal/config.yml b/benchmarks/lfhcal/config.yml index a6dfcbe4..8dfc7130 100644 --- a/benchmarks/lfhcal/config.yml +++ b/benchmarks/lfhcal/config.yml @@ -16,7 +16,7 @@ bench:lfhcal: needs: - ["sim:lfhcal"] script: - - snakemake --cores 3 lfhcal_local + - snakemake $SNAKEMAKE_FLAGS --cores 3 lfhcal_local collect_results:lfhcal: extends: .det_benchmark @@ -26,7 +26,7 @@ collect_results:lfhcal: script: - ls -lrht - mv results{,_save}/ # move results directory out of the way to preserve it - - snakemake --cores 1 --delete-all-output lfhcal_local + - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output lfhcal_local - mv results{_save,}/ bench:lfhcal_campaigns: @@ -34,7 +34,7 @@ bench:lfhcal_campaigns: stage: benchmarks #when: manual script: - - snakemake --cores 1 lfhcal_campaigns + - snakemake $SNAKEMAKE_FLAGS --cores 1 lfhcal_campaigns collect_results:lfhcal_campaigns: extends: .det_benchmark @@ -44,5 +44,5 @@ collect_results:lfhcal_campaigns: script: - ls -lrht - mv results{,_save}/ # move results directory out of the way to preserve it - - snakemake --cores 1 --delete-all-output lfhcal_campaigns + - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output lfhcal_campaigns - mv results{_save,}/ diff --git a/benchmarks/material_scan/config.yml b/benchmarks/material_scan/config.yml index f8bc0822..3c87305e 100644 --- a/benchmarks/material_scan/config.yml +++ b/benchmarks/material_scan/config.yml @@ -2,7 +2,7 @@ bench:material_scan: extends: .det_benchmark stage: benchmarks script: - - snakemake --cores 1 epic_craterlake/results/material_scan_details.pdf + - snakemake $SNAKEMAKE_FLAGS --cores 1 epic_craterlake/results/material_scan_details.pdf collect_results:material_scan: extends: .det_benchmark diff --git a/benchmarks/tracking_performances/config.yml b/benchmarks/tracking_performances/config.yml index 90a0d0f6..2419000b 100644 --- a/benchmarks/tracking_performances/config.yml +++ b/benchmarks/tracking_performances/config.yml @@ -18,7 +18,7 @@ bench:tracking_performance: needs: - ["sim:tracking_performance"] script: - - snakemake --cores 3 tracking_performance_local + - snakemake $SNAKEMAKE_FLAGS --cores 3 tracking_performance_local collect_results:tracking_performance: extends: .det_benchmark @@ -27,14 +27,14 @@ collect_results:tracking_performance: - "bench:tracking_performance" script: - ls -lrht - - snakemake --cores 1 --delete-all-output tracking_performance_local + - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output tracking_performance_local bench:tracking_performance_campaigns: extends: .det_benchmark stage: benchmarks #when: manual script: - - snakemake --cores 1 tracking_performance_campaigns + - snakemake $SNAKEMAKE_FLAGS --cores 1 tracking_performance_campaigns collect_results:tracking_performance_campaigns: extends: .det_benchmark @@ -44,5 +44,5 @@ collect_results:tracking_performance_campaigns: script: - ls -lrht - mv results{,_save}/ # move results directory out of the way to preserve it - - snakemake --cores 1 --delete-all-output tracking_performance_campaigns + - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output tracking_performance_campaigns - mv results{_save,}/ diff --git a/benchmarks/tracking_performances_dis/config.yml b/benchmarks/tracking_performances_dis/config.yml index 961df549..066e8df7 100644 --- a/benchmarks/tracking_performances_dis/config.yml +++ b/benchmarks/tracking_performances_dis/config.yml @@ -2,7 +2,7 @@ sim:tracking_performances_dis: extends: .det_benchmark stage: simulate script: - - snakemake --cache --cores 5 results/tracking_performances_dis/epic_craterlake_tracking_only/pythia8NCDIS_18x275_minQ2=1_combined_5/hists.root + - snakemake $SNAKEMAKE_FLAGS --cores 5 results/tracking_performances_dis/epic_craterlake_tracking_only/pythia8NCDIS_18x275_minQ2=1_combined_5/hists.root retry: max: 2 when: @@ -14,7 +14,7 @@ bench:tracking_performances_dis: needs: - ["sim:tracking_performances_dis"] script: - - snakemake --cores 1 trk_dis_run_locally_trk_only + - snakemake $SNAKEMAKE_FLAGS --cores 1 trk_dis_run_locally_trk_only # avoid uploading intermediate results - find results/tracking_performances_dis/ -mindepth 2 -maxdepth 2 -type d ! -name "*combined*" | xargs rm -rfv @@ -26,7 +26,7 @@ collect_results:tracking_performances_dis: script: - ls -lrht - mv results{,_save}/ # move results directory out of the way to preserve it - - snakemake --cores 1 --delete-all-output trk_dis_run_locally_trk_only + - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output trk_dis_run_locally_trk_only - mv results{_save,}/ # convert to png - | diff --git a/benchmarks/zdc_lambda/config.yml b/benchmarks/zdc_lambda/config.yml index e6b16031..609942f5 100644 --- a/benchmarks/zdc_lambda/config.yml +++ b/benchmarks/zdc_lambda/config.yml @@ -12,7 +12,7 @@ sim:zdc_lambda: - P: 250 - P: 275 script: - - snakemake --cores 5 sim_output/zdc_lambda/epic_zdc_sipm_on_tile_only_rec_lambda_dec_${P}GeV_{0,1,2,3,4}.edm4eic.root + - snakemake $SNAKEMAKE_FLAGS --cores 5 sim_output/zdc_lambda/epic_zdc_sipm_on_tile_only_rec_lambda_dec_${P}GeV_{0,1,2,3,4}.edm4eic.root retry: max: 2 when: @@ -23,7 +23,7 @@ bench:zdc_lambda: stage: benchmarks needs: ["sim:zdc_lambda"] script: - - snakemake --cores 1 results/epic_zdc_sipm_on_tile_only/zdc_lambda + - snakemake $SNAKEMAKE_FLAGS --cores 1 results/epic_zdc_sipm_on_tile_only/zdc_lambda collect_results:zdc_lambda: extends: .det_benchmark @@ -32,5 +32,5 @@ collect_results:zdc_lambda: script: - ls -al - mv results{,_save}/ # move results directory out of the way to preserve it - - snakemake --cores 1 --delete-all-output results/epic_zdc_sipm_on_tile_only/zdc_lambda + - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output results/epic_zdc_sipm_on_tile_only/zdc_lambda - mv results{_save,}/ diff --git a/benchmarks/zdc_lyso/config.yml b/benchmarks/zdc_lyso/config.yml index eb0a763b..990aaef2 100644 --- a/benchmarks/zdc_lyso/config.yml +++ b/benchmarks/zdc_lyso/config.yml @@ -2,7 +2,7 @@ sim:zdc_lyso: extends: .det_benchmark stage: simulate script: - - snakemake --cache --cores 1 zdc_lyso_local + - snakemake $SNAKEMAKE_FLAGS --cores 1 zdc_lyso_local retry: max: 2 when: @@ -16,5 +16,5 @@ collect_results:zdc_lyso: script: - ls -lrht - mv results{,_save}/ # move results directory out of the way to preserve it - - snakemake --cores 1 --delete-all-output zdc_lyso_local + - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output zdc_lyso_local - mv results{_save,}/ diff --git a/benchmarks/zdc_photon/config.yml b/benchmarks/zdc_photon/config.yml index d62cef61..9814792a 100644 --- a/benchmarks/zdc_photon/config.yml +++ b/benchmarks/zdc_photon/config.yml @@ -12,7 +12,7 @@ sim:zdc_photon: - P: 200 - P: 275 script: - - snakemake --cores 5 sim_output/zdc_photon/epic_zdc_sipm_on_tile_only_rec_zdc_photon_${P}GeV_{0,1,2,3,4}.edm4eic.root + - snakemake $SNAKEMAKE_FLAGS --cores 5 sim_output/zdc_photon/epic_zdc_sipm_on_tile_only_rec_zdc_photon_${P}GeV_{0,1,2,3,4}.edm4eic.root retry: max: 2 when: diff --git a/benchmarks/zdc_pi0/config.yml b/benchmarks/zdc_pi0/config.yml index f0b8993e..0f9a68c4 100644 --- a/benchmarks/zdc_pi0/config.yml +++ b/benchmarks/zdc_pi0/config.yml @@ -9,7 +9,7 @@ sim:zdc_pi0: - P: 130 - P: 160 script: - - snakemake --cores 5 sim_output/zdc_pi0/epic_zdc_sipm_on_tile_only_rec_zdc_pi0_${P}GeV_{0,1,2,3,4}.edm4eic.root + - snakemake $SNAKEMAKE_FLAGS --cores 5 sim_output/zdc_pi0/epic_zdc_sipm_on_tile_only_rec_zdc_pi0_${P}GeV_{0,1,2,3,4}.edm4eic.root retry: max: 2 when: diff --git a/benchmarks/zdc_sigma/config.yml b/benchmarks/zdc_sigma/config.yml index af204d42..17a84f39 100644 --- a/benchmarks/zdc_sigma/config.yml +++ b/benchmarks/zdc_sigma/config.yml @@ -12,7 +12,7 @@ sim:zdc_sigma: - P: 250 - P: 275 script: - - snakemake --cores 5 sim_output/zdc_sigma/epic_zdc_sipm_on_tile_only_rec_sigma_dec_${P}GeV_{0,1,2,3,4}.edm4eic.root + - snakemake $SNAKEMAKE_FLAGS --cores 5 sim_output/zdc_sigma/epic_zdc_sipm_on_tile_only_rec_sigma_dec_${P}GeV_{0,1,2,3,4}.edm4eic.root retry: max: 2 when: @@ -23,7 +23,7 @@ bench:zdc_sigma: stage: benchmarks needs: ["sim:zdc_sigma"] script: - - snakemake --cores 1 results/epic_zdc_sipm_on_tile_only/zdc_sigma + - snakemake $SNAKEMAKE_FLAGS --cores 1 results/epic_zdc_sipm_on_tile_only/zdc_sigma collect_results:zdc_sigma: extends: .det_benchmark @@ -32,5 +32,5 @@ collect_results:zdc_sigma: script: - ls -al - mv results{,_save}/ # move results directory out of the way to preserve it - - snakemake --cores 1 --delete-all-output results/epic_zdc_sipm_on_tile_only/zdc_sigma + - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output results/epic_zdc_sipm_on_tile_only/zdc_sigma - mv results{_save,}/ From df6a0b3f9709a63dace41e3ef24d0f3e056ad489 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 19 Oct 2024 22:26:14 -0400 Subject: [PATCH 20/24] .gitlab-ci.yml: remove --rerun-triggers It appears that presence of mtime in the list is simply ignored by Snakemake. Probably, it's not possible to run with ignoring mtimes. --- .gitlab-ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 7a19fd1a..36abfc88 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -5,7 +5,7 @@ variables: DETECTOR_CONFIG: epic_craterlake GITHUB_SHA: '' GITHUB_REPOSITORY: '' - SNAKEMAKE_FLAGS: '--cache --rerun-triggers code input params software-env' + SNAKEMAKE_FLAGS: '--cache' workflow: name: '$PIPELINE_NAME' From 6e29bea1f8f93d12231c1d583956a65414b4a652 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 19 Oct 2024 22:27:27 -0400 Subject: [PATCH 21/24] backwards_ecal: enable caching for the reco step --- benchmarks/backwards_ecal/Snakefile | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/benchmarks/backwards_ecal/Snakefile b/benchmarks/backwards_ecal/Snakefile index da90f2e2..0edae3a7 100644 --- a/benchmarks/backwards_ecal/Snakefile +++ b/benchmarks/backwards_ecal/Snakefile @@ -69,9 +69,13 @@ rule backwards_ecal_recon: "sim_output/backwards_ecal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.eicrecon.tree.edm4eic.root.log", wildcard_constraints: INDEX="\d{4}", + params: + DETECTOR_CONFIG=lambda wildcards: wildcards.DETECTOR_CONFIG, + EICRECON_HASH=get_spack_package_hash("eicrecon"), + cache: True shell: """ set -m # monitor mode to prevent lingering processes -exec env DETECTOR_CONFIG={wildcards.DETECTOR_CONFIG} \ +exec env DETECTOR_CONFIG={params.DETECTOR_CONFIG} \ eicrecon {input} -Ppodio:output_file={output} \ -Ppodio:output_collections=MCParticles,EcalEndcapNRecHits,EcalEndcapNClusters """ From 89c8682bb736143aef08b2daea6c560177def569 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Wed, 23 Oct 2024 00:42:32 -0400 Subject: [PATCH 22/24] femc_*: use maxfev=10000 --- benchmarks/femc_electron/analysis/femc_electron_plots.py | 6 +++--- benchmarks/femc_photon/analysis/femc_photon_plots.py | 6 +++--- benchmarks/femc_pi0/analysis/femc_pi0_plots.py | 6 +++--- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/benchmarks/femc_electron/analysis/femc_electron_plots.py b/benchmarks/femc_electron/analysis/femc_electron_plots.py index cc34633d..8979e0c7 100644 --- a/benchmarks/femc_electron/analysis/femc_electron_plots.py +++ b/benchmarks/femc_electron/analysis/femc_electron_plots.py @@ -132,7 +132,7 @@ def gauss(x, A,mu, sigma): 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)) + coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000) #res=np.abs(coeff[2]/coeff[1]) if p==50: xx=np.linspace(15*p/20,22*p/20, 100) @@ -151,7 +151,7 @@ def gauss(x, A,mu, sigma): 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) +coeff, var_matrix = curve_fit(fnc, pvals, res, p0=p0, sigma=dres, maxfev=10000) xx=np.linspace(7, 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() @@ -190,7 +190,7 @@ def gauss(x, A,mu, sigma): p0=(100, p, 3) try: - coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma)) + coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000) if abs(coeff[1])>100 or np.sqrt(var_matrix[1][1])>100: continue pvals.append(p) diff --git a/benchmarks/femc_photon/analysis/femc_photon_plots.py b/benchmarks/femc_photon/analysis/femc_photon_plots.py index 8170f4a1..176eb663 100644 --- a/benchmarks/femc_photon/analysis/femc_photon_plots.py +++ b/benchmarks/femc_photon/analysis/femc_photon_plots.py @@ -131,7 +131,7 @@ def gauss(x, A,mu, sigma): 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)) + coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000) #res=np.abs(coeff[2]/coeff[1]) if p==50: xx=np.linspace(15*p/20,22*p/20, 100) @@ -150,7 +150,7 @@ def gauss(x, A,mu, sigma): 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) +coeff, var_matrix = curve_fit(fnc, pvals, res, p0=p0, sigma=dres, maxfev=10000) xx=np.linspace(7, 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() @@ -188,7 +188,7 @@ def gauss(x, A,mu, sigma): 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)) + coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000) if abs(coeff[1])>100 or np.sqrt(var_matrix[1][1])>100: continue pvals.append(p) diff --git a/benchmarks/femc_pi0/analysis/femc_pi0_plots.py b/benchmarks/femc_pi0/analysis/femc_pi0_plots.py index cdf9d231..0bdf26e3 100644 --- a/benchmarks/femc_pi0/analysis/femc_pi0_plots.py +++ b/benchmarks/femc_pi0/analysis/femc_pi0_plots.py @@ -167,7 +167,7 @@ def gauss(x, A,mu, sigma): 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)) + coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000) #res=np.abs(coeff[2]/coeff[1]) if p==50: xx=np.linspace(15*p/20,22*p/20, 100) @@ -186,7 +186,7 @@ def gauss(x, A,mu, sigma): 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) +coeff, var_matrix = curve_fit(fnc, pvals, res, p0=p0, sigma=dres, maxfev=10000) 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() @@ -225,7 +225,7 @@ def gauss(x, A,mu, sigma): p0=(100, p, 3) try: - coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma)) + coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000) if abs(coeff[1])>100 or np.sqrt(var_matrix[1][1])>100: continue pvals.append(p) From 958f6cb777739330ff27b0a26bdf734c88b08a7b Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Wed, 23 Oct 2024 00:46:22 -0400 Subject: [PATCH 23/24] insert_*: use maxfev=10000 --- benchmarks/insert_muon/analysis/muon_plots.py | 2 +- benchmarks/insert_neutron/analysis/neutron_plots.py | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/benchmarks/insert_muon/analysis/muon_plots.py b/benchmarks/insert_muon/analysis/muon_plots.py index 5c81b52d..ac2aa100 100644 --- a/benchmarks/insert_muon/analysis/muon_plots.py +++ b/benchmarks/insert_muon/analysis/muon_plots.py @@ -55,7 +55,7 @@ def Landau(x, normalization,location,stdev): 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))) + sigma=list(np.sqrt(y[slc])+(y[slc]==0)), maxfev=10000) print(coeff) xx=np.linspace(0,.7, 100) MIP=coeff[1]/1000 diff --git a/benchmarks/insert_neutron/analysis/neutron_plots.py b/benchmarks/insert_neutron/analysis/neutron_plots.py index b9d2002a..424295cf 100644 --- a/benchmarks/insert_neutron/analysis/neutron_plots.py +++ b/benchmarks/insert_neutron/analysis/neutron_plots.py @@ -79,7 +79,7 @@ def gauss(x, A,mu, sigma): fnc=gauss sigma=np.sqrt(y[slc])+(y[slc]==0) p0=(100, 0, 5) -coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,sigma=list(sigma)) +coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000) xx=np.linspace(-5,5,100) plt.plot(xx,fnc(xx,*coeff)) # except: @@ -104,7 +104,7 @@ def gauss(x, A,mu, sigma): #print(bc[slc],y[slc]) sigma=np.sqrt(y[slc])+(y[slc]==0) try: - coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,sigma=list(sigma)) + coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000) sigmas.append(np.abs(coeff[2])) dsigmas.append(np.sqrt(var_matrix[2][2])) xvals.append(p) @@ -151,7 +151,7 @@ def gauss(x, A,mu, sigma): sigma=np.sqrt(y[slc])+0.5*(y[slc]==0) p0=(100, np.mean(r), np.std(r)) try: - coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma)) + coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000) res=np.abs(coeff[2]/coeff[1]) if res Date: Wed, 23 Oct 2024 00:51:22 -0400 Subject: [PATCH 24/24] zdc_*: use maxfev=10000 --- benchmarks/zdc_lambda/analysis/lambda_plots.py | 12 ++++++------ benchmarks/zdc_photon/analysis/zdc_photon_plots.py | 8 ++++---- benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py | 12 ++++++------ 3 files changed, 16 insertions(+), 16 deletions(-) diff --git a/benchmarks/zdc_lambda/analysis/lambda_plots.py b/benchmarks/zdc_lambda/analysis/lambda_plots.py index bf855d9d..931a44b3 100644 --- a/benchmarks/zdc_lambda/analysis/lambda_plots.py +++ b/benchmarks/zdc_lambda/analysis/lambda_plots.py @@ -191,7 +191,7 @@ def gauss(x, A,mu, sigma): fnc=gauss p0=[100, 0, 0.05] coeff, var_matrix = curve_fit(fnc, bc[slc], y[slc], p0=p0, - sigma=np.sqrt(y[slc])+(y[slc]==0)) + sigma=np.sqrt(y[slc])+(y[slc]==0), maxfev=10000) x=np.linspace(-1, 1) plt.plot(x, gauss(x, *coeff), color='tab:orange') plt.xlabel("$\\theta^{*\\rm recon}_{\\Lambda}-\\theta^{*\\rm truth}_{\\Lambda}$ [mrad]") @@ -214,7 +214,7 @@ def gauss(x, A,mu, sigma): #print(bc[slc],y[slc]) sigma=np.sqrt(y[slc])+(y[slc]==0) try: - coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,sigma=list(sigma)) + coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000) sigmas.append(coeff[2]) dsigmas.append(np.sqrt(var_matrix[2][2])) xvals.append(p) @@ -259,7 +259,7 @@ def gauss(x, A,mu, sigma): fnc=gauss p0=[100, 0, 1] coeff, var_matrix = curve_fit(fnc, bc[slc], y[slc], p0=p0, - sigma=np.sqrt(y[slc])+(y[slc]==0)) + sigma=np.sqrt(y[slc])+(y[slc]==0), maxfev=10000) x=np.linspace(-5, 5) plt.plot(x, gauss(x, *coeff), color='tab:orange') print(coeff[2], np.sqrt(var_matrix[2][2])) @@ -284,7 +284,7 @@ def gauss(x, A,mu, sigma): #print(bc[slc],y[slc]) sigma=np.sqrt(y[slc])+(y[slc]==0) try: - coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,sigma=list(sigma)) + coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000) sigmas.append(coeff[2]) dsigmas.append(np.sqrt(var_matrix[2][2])) xvals.append(p) @@ -327,7 +327,7 @@ def gauss(x, A,mu, sigma): fnc=gauss p0=[100, lambda_mass, 0.04] coeff, var_matrix = curve_fit(fnc, bc[slc], y[slc], p0=p0, - sigma=np.sqrt(y[slc])+(y[slc]==0)) + sigma=np.sqrt(y[slc])+(y[slc]==0), maxfev=10000) x=np.linspace(0.8, 1.3, 200) plt.plot(x, gauss(x, *coeff), color='tab:orange') print(coeff[2], np.sqrt(var_matrix[2][2])) @@ -350,7 +350,7 @@ def gauss(x, A,mu, sigma): p0=[100, lambda_mass, 0.05] try: coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0, - sigma=list(np.sqrt(y[slc])+(y[slc]==0))) + sigma=list(np.sqrt(y[slc])+(y[slc]==0)), maxfev=10000) x=np.linspace(0.8, 1.3, 200) sigmas.append(coeff[2]) dsigmas.append(np.sqrt(var_matrix[2][2])) diff --git a/benchmarks/zdc_photon/analysis/zdc_photon_plots.py b/benchmarks/zdc_photon/analysis/zdc_photon_plots.py index 94a93101..a2cf5339 100644 --- a/benchmarks/zdc_photon/analysis/zdc_photon_plots.py +++ b/benchmarks/zdc_photon/analysis/zdc_photon_plots.py @@ -57,7 +57,7 @@ def gauss(x, A,mu, sigma): p0=[100, p, 10] #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))) + sigma=list(np.sqrt(y[slc])+(y[slc]==0)), maxfev=10000) if p==100: xx=np.linspace(p*0.75,p*1.25, 100) plt.plot(xx, fnc(xx,*coeff)) @@ -78,7 +78,7 @@ def gauss(x, A,mu, sigma): fnc=lambda E,a: a/np.sqrt(E) #pvals, resvals, dresvals coeff, var_matrix = curve_fit(fnc, pvals, resvals, p0=(1,), - sigma=dresvals) + sigma=dresvals, maxfev=10000) xx=np.linspace(15, 275, 100) plt.plot(xx, fnc(xx, *coeff), label=f'fit: $\\frac{{{coeff[0]*100:.0f}\\%}}{{\\sqrt{{E}}}}$') @@ -129,7 +129,7 @@ def gauss(x, A,mu, sigma): p0=[100, 0, 0.1] #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))) + sigma=list(np.sqrt(y[slc])+(y[slc]==0)), maxfev=10000) if p==100: xx=np.linspace(-0.5,0.5, 100) plt.plot(xx, fnc(xx,*coeff)) @@ -143,7 +143,7 @@ def gauss(x, A,mu, sigma): fnc=lambda E,a, b: np.hypot(a/np.sqrt(E), b) #pvals, resvals, dresvals coeff, var_matrix = curve_fit(fnc, pvals, resvals, p0=(1,.1), - sigma=dresvals) + sigma=dresvals, maxfev=10000) xx=np.linspace(15, 275, 100) diff --git a/benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py b/benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py index 1b54430e..6d0cdddd 100644 --- a/benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py +++ b/benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py @@ -57,7 +57,7 @@ def gauss(x, A,mu, sigma): p0=[100, p, 10] #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))) + sigma=list(np.sqrt(y[slc])+(y[slc]==0)), maxfev=10000) if p==100: xx=np.linspace(p*0.5,p*1.5, 100) plt.plot(xx, fnc(xx,*coeff)) @@ -76,7 +76,7 @@ def gauss(x, A,mu, sigma): fnc=lambda E,a: a/np.sqrt(E) #pvals, resvals, dresvals coeff, var_matrix = curve_fit(fnc, pvals, resvals, p0=(1,), - sigma=dresvals) + sigma=dresvals, maxfev=10000) xx=np.linspace(55, 200, 100) plt.plot(xx, fnc(xx, *coeff), label=f'fit: $\\frac{{{coeff[0]:.2f}\\%}}{{\\sqrt{{E}}}}$') plt.legend() @@ -133,7 +133,7 @@ def gauss(x, A,mu, sigma): p0=[100, 0, 0.1] #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))) + sigma=list(np.sqrt(y[slc])+(y[slc]==0)), maxfev=10000) if p==100: xx=np.linspace(-0.5,0.5, 100) plt.plot(xx, fnc(xx,*coeff)) @@ -148,7 +148,7 @@ def gauss(x, A,mu, sigma): fnc=lambda E,a: a/np.sqrt(E) #pvals, resvals, dresvals coeff, var_matrix = curve_fit(fnc, pvals, resvals, p0=(1,), - sigma=dresvals) + sigma=dresvals, maxfev=10000) xx=np.linspace(55, 200, 100) @@ -201,7 +201,7 @@ def gauss(x, A,mu, sigma): p0=[100, .135, 0.2] #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))) + sigma=list(np.sqrt(y[slc])+(y[slc]==0)), maxfev=10000) if p==100: xx=np.linspace(0,0.2) plt.plot(xx, fnc(xx,*coeff)) @@ -218,7 +218,7 @@ def gauss(x, A,mu, sigma): fnc=lambda E,a,b: a+b*E #pvals, resvals, dresvals coeff, var_matrix = curve_fit(fnc, pvals, resvals, p0=(1,1), - sigma=dresvals) + sigma=dresvals, maxfev=10000) xx=np.linspace(55, 200, 100) #plt.plot(xx, fnc(xx, *coeff), label=f'fit: ${coeff[0]*1000:.1f}+{coeff[1]*1000:.4f}\\times E$ MeV') plt.plot(xx, fnc(xx, *coeff), label=f'fit: $({coeff[0]*1000:.1f}+{coeff[1]*1000:.4f}\\times [E\,in\,GeV])$ MeV')