diff --git a/benchmarks/femc_electron/Snakefile b/benchmarks/femc_electron/Snakefile index 49b2995..0a11626 100644 --- a/benchmarks/femc_electron/Snakefile +++ b/benchmarks/femc_electron/Snakefile @@ -1,6 +1,6 @@ -rule generate: +rule femc_electron_generate: input: - script="src/gen_particles.cxx", + script="benchmarks/femc_electron/analysis/gen_particles.cxx", params: NEVENTS_GEN=1000, th_max=28, @@ -13,7 +13,7 @@ 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 simulate: +rule femc_electron_simulate: input: GEN_FILE="sim_output/femc_electron/e-_{P}GeV.hepmc" params: @@ -32,7 +32,7 @@ npsim \ --outputFile {output.SIM_FILE} """ -rule recon: +rule femc_electron_recon: input: SIM_FILE="sim_output/femc_electron/{DETECTOR_CONFIG}_sim_e-_{P}GeV.edm4hep.root" output: @@ -43,7 +43,7 @@ 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 analysis: +rule femc_electron_analysis: input: expand("sim_output/femc_electron/{DETECTOR_CONFIG}_rec_e-_{P}GeV.edm4hep.root", P=[20, 30, 40, 50, 60, 70, 80], diff --git a/benchmarks/femc_electron/analysis/femc_electron_plots.py b/benchmarks/femc_electron/analysis/femc_electron_plots.py index 55bf7da..4b7800f 100644 --- a/benchmarks/femc_electron/analysis/femc_electron_plots.py +++ b/benchmarks/femc_electron/analysis/femc_electron_plots.py @@ -17,7 +17,7 @@ pass import uproot as ur -arrays_sim={p:ur.open(f'/Users/spaul/EIC/detector_benchmarks/sim_output/femc_electron/epic_craterlake_rec_e-_{p}GeV.edm4hep.root:events').arrays() for p in (20, 30, 40, 50, 60,70,80)} +arrays_sim={p:ur.open(f'/sim_output/femc_electron/{config}_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] @@ -66,24 +66,42 @@ pvals=[] #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)): + 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(ak.flatten(-a['EcalEndcapPClusters.hits_begin']+a['EcalEndcapPClusters.hits_end']), - range=(0,100), bins=21, histtype='step', label=f"E={p} GeV") + y,x,_=plt.hist(n, range=(0,100), bins=100, histtype='step', label=f"E={p} GeV") plt.ylabel("events") - plt.xlabel("# hits per Ecal cluster") + plt.xlabel("# hits in cluster") plt.title(f"e-, E={p} GeV") - - n=ak.flatten(-a['EcalEndcapPClusters.hits_begin']+a['EcalEndcapPClusters.hits_end']) 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 per Ecal cluster [mean$\\pm$std]") +plt.ylabel("# hits in cluster [mean$\\pm$std]") plt.ylim(0) plt.savefig(outdir+"/nhits_per_cluster.pdf") diff --git a/benchmarks/femc_photon/Snakefile b/benchmarks/femc_photon/Snakefile new file mode 100644 index 0000000..c5259f1 --- /dev/null +++ b/benchmarks/femc_photon/Snakefile @@ -0,0 +1,59 @@ +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.edm4hep.root" + shell: + """ +NEVENTS_REC=1000 +eicrecon {input.SIM_FILE} -Ppodio:output_file={output.REC_FILE} -Pdd4hep:xml_files=$DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml -Ppodio:output_collections=MCParticles,HcalEndcapPInsertRecHits,HcalEndcapPInsertClusters,HcalEndcapPInsertSubcellHits,EcalEndcapPInsertRecHits,EcalEndcapPInsertClusters,EcalEndcapPInsertRecHits,EcalEndcapPClusters -Pjana:nevents=$NEVENTS_REC +""" + +rule femc_photon_analysis: + input: + expand("sim_output/femc_photon/{DETECTOR_CONFIG}_rec_photon_{P}GeV.edm4hep.root", + P=[20, 30, 40, 50, 60, 70, 80], + DETECTOR_CONFIG=["{DETECTOR_CONFIG}"], + particle=["{particle}"]), + 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 0000000..7eb5e41 --- /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/{config}_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']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 0000000..9ad6f6a --- /dev/null +++ b/benchmarks/femc_photon/config.yml @@ -0,0 +1,35 @@ +sim:femc_photon: + stage: simulate + extends: .det_benchmark + parallel: + matrix: + - 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 0000000..de2ef6d --- /dev/null +++ b/benchmarks/femc_pi0/Snakefile @@ -0,0 +1,59 @@ +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.edm4hep.root" + shell: + """ +NEVENTS_REC=1000 +eicrecon {input.SIM_FILE} -Ppodio:output_file={output.REC_FILE} -Pdd4hep:xml_files=$DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml -Ppodio:output_collections=MCParticles,HcalEndcapPInsertRecHits,HcalEndcapPInsertClusters,HcalEndcapPInsertSubcellHits,EcalEndcapPInsertRecHits,EcalEndcapPInsertClusters,EcalEndcapPInsertRecHits,EcalEndcapPClusters -Pjana:nevents=$NEVENTS_REC +""" + +rule femc_pi0_analysis: + input: + expand("sim_output/femc_pi0/{DETECTOR_CONFIG}_rec_pi0_{P}GeV.edm4hep.root", + P=[20, 30, 40, 50, 60, 70, 80], + DETECTOR_CONFIG=["{DETECTOR_CONFIG}"], + particle=["{particle}"]), + 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 0000000..7e24b2a --- /dev/null +++ b/benchmarks/femc_pi0/analysis/femc_pi0_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'/Users/spaul/EIC/detector_benchmarks/sim_output/femc_electron/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_pi0/config.yml b/benchmarks/femc_pi0/config.yml new file mode 100644 index 0000000..f6215da --- /dev/null +++ b/benchmarks/femc_pi0/config.yml @@ -0,0 +1,35 @@ +sim:femc_pi0: + stage: simulate + extends: .det_benchmark + parallel: + matrix: + - 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,}/