From c72c399b3bba78d55114f6cafd18a7ebe79cd47b Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Wed, 2 Oct 2024 08:39:47 -0400 Subject: [PATCH] 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 << " ("<