diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 6d7490c7a1353e9f3229e17fb68ae497fc30c011..2ced7dee3575fc8a3d15e0b0384bcba048b5fe9d 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -119,6 +119,7 @@ include:
- local: 'benchmarks/dis/config.yml'
#- local: 'benchmarks/dvmp/config.yml'
- local: 'benchmarks/dvcs/config.yml'
+ - local: 'benchmarks/lambda/config.yml'
- local: 'benchmarks/tcs/config.yml'
- local: 'benchmarks/u_omega/config.yml'
- local: 'benchmarks/single/config.yml'
@@ -131,6 +132,7 @@ summary:
- "demp:results"
- "dis:results"
- "dvcs:results"
+ - "lambda:results"
- "tcs:results"
- "u_omega:results"
- "single:results"
diff --git a/Snakefile b/Snakefile
index 73b03c2b57ca94f25bb0d80d54cf8748a32d3b8f..1825ce1a250623bf12a318232661217040448544 100644
--- a/Snakefile
+++ b/Snakefile
@@ -42,4 +42,5 @@ ddsim \
include: "benchmarks/diffractive_vm/Snakefile"
include: "benchmarks/dis/Snakefile"
+include: "benchmarks/lambda/Snakefile"
include: "benchmarks/demp/Snakefile"
diff --git a/benchmarks/lambda/Snakefile b/benchmarks/lambda/Snakefile
new file mode 100644
index 0000000000000000000000000000000000000000..e565e8ea16ee54351cfafb6002f240bc8b570551
--- /dev/null
+++ b/benchmarks/lambda/Snakefile
@@ -0,0 +1,64 @@
+rule lambda_generate:
+ input:
+ script="benchmarks/lambda/analysis/gen_lambda_decay.cxx",
+ params:
+ NEVENTS_GEN=100000,
+ output:
+ GEN_FILE="results/lambda/lambda_decay_{P}GeV.hepmc"
+ shell:
+ """
+mkdir -p results/lambda
+root -l -b -q '{input.script}({params.NEVENTS_GEN},0,"{output.GEN_FILE}",{wildcards.P},{wildcards.P})'
+"""
+
+rule lambda_simulate:
+ input:
+ GEN_FILE="results/lambda/lambda_decay_{P}GeV.hepmc"
+ params:
+ PHYSICS_LIST="FTFP_BERT"
+ output:
+ SIM_FILE="results/lambda/{DETECTOR_CONFIG}_sim_lambda_dec_{P}GeV.edm4hep.root"
+ shell:
+ """
+if [[ {wildcards.P} -gt 225 ]]; then
+ NEVENTS_SIM=1000
+else
+ NEVENTS_SIM=1000
+fi
+# 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 lambda_recon:
+ input:
+ SIM_FILE="results/lambda/{DETECTOR_CONFIG}_sim_lambda_dec_{P}GeV.edm4hep.root"
+ output:
+ REC_FILE="results/lambda/{DETECTOR_CONFIG}_rec_lambda_dec_{P}GeV.edm4hep.root"
+ shell:
+ """
+if [[ {wildcards.P} -gt 225 ]]; then
+ NEVENTS_REC=1000
+else
+ NEVENTS_REC=1000
+fi
+eicrecon {input.SIM_FILE} -Ppodio:output_file={output.REC_FILE} -Pdd4hep:xml_files=$DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml -Ppodio:output_include_collections=MCParticles,HcalFarForwardZDCClusters,HcalFarForwardZDCRecHits,HcalFarForwardZDCSubcellHits -Pjana:nevents=$NEVENTS_REC
+"""
+
+rule lambda_analysis:
+ input:
+ expand("results/lambda/{DETECTOR_CONFIG}_rec_lambda_dec_{P}GeV.edm4hep.root",
+ P=[100, 125, 150,175, 200, 225, 250, 275],
+ DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]),
+ script="benchmarks/lambda/analysis/lambda_plots.py",
+ output:
+ results_dir=directory("results/lambda/results_{DETECTOR_CONFIG}_lambda_dec"),
+ shell:
+ """
+mkdir -p {output.results_dir}
+python {input.script} {output.results_dir}
+"""
diff --git a/benchmarks/lambda/analysis/gen_lambda_decay.cxx b/benchmarks/lambda/analysis/gen_lambda_decay.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..567eda5bcf5497f1b64596745ef5820e7a54eff9
--- /dev/null
+++ b/benchmarks/lambda/analysis/gen_lambda_decay.cxx
@@ -0,0 +1,239 @@
+#include "HepMC3/GenEvent.h"
+#include "HepMC3/ReaderAscii.h"
+#include "HepMC3/WriterAscii.h"
+#include "HepMC3/Print.h"
+
+#include "TRandom3.h"
+#include "TVector3.h"
+
+#include <TDatabasePDG.h>
+#include <TParticlePDG.h>
+
+#include <iostream>
+#include <random>
+#include <TMath.h>
+
+using namespace HepMC3;
+
+std::tuple<double, int, double> GetParticleInfo(TDatabasePDG* pdg, TString particle_name)
+{
+ TParticlePDG *particle = pdg->GetParticle(particle_name);
+ const double mass = particle->Mass();
+ const int pdgID = particle->PdgCode();
+ const double lifetime = particle->Lifetime();
+ return std::make_tuple(mass, pdgID, lifetime);
+}
+// Calculates the decay length of a particle. Samples from an exponential decay.
+double GetDecayLength(TRandom3* r1, double lifetime, double mass, double momentum_magnitude)
+{
+ double c_speed = TMath::C() * 1000.; // speed of light im mm/sec
+ double average_decay_length = (momentum_magnitude/mass) * lifetime * c_speed;
+ return r1->Exp(average_decay_length);
+}
+
+// Generate single lambda mesons and decay them to a neutron + 2 photons
+void gen_lambda_decay(int n_events = 100000, UInt_t seed = 0, char* out_fname = "lambda_decay.hepmc",
+ double p_min = 100., // in GeV/c
+ double p_max = 275.) // in GeV/c
+{
+
+ const double theta_min = 0.0; // in mRad
+ const double theta_max = 3.0; // in mRad
+ //const double p_min = 100.; // in GeV/c
+ //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
+ TRandom3 *r1 = new TRandom3(seed); //Default = 0, which uses clock to set seed
+ cout<<"Random number seed is "<<r1->GetSeed()<<"!"<<endl;
+
+ // Getting generated particle information
+ TDatabasePDG *pdg = new TDatabasePDG();
+
+ auto lambda_info = GetParticleInfo(pdg, "Lambda0");
+ double lambda_mass = std::get<0>(lambda_info);
+ int lambda_pdgID = std::get<1>(lambda_info);
+ double lambda_lifetime = std::get<2>(lambda_info);
+
+ auto neutron_info = GetParticleInfo(pdg, "neutron");
+ double neutron_mass = std::get<0>(neutron_info);
+ int neutron_pdgID = std::get<1>(neutron_info);
+
+ auto pi0_info = GetParticleInfo(pdg, "pi0");
+ double pi0_mass = std::get<0>(pi0_info);
+ int pi0_pdgID = std::get<1>(pi0_info);
+ double pi0_lifetime = std::get<2>(pi0_info);
+
+ auto photon_info = GetParticleInfo(pdg, "gamma");
+ 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++) {
+
+ //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 2212 - proton
+ GenParticlePtr p1 =
+ std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
+ GenParticlePtr p2 = std::make_shared<GenParticle>(
+ FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
+
+ // Define momentum with respect to EIC proton beam direction
+ Double_t lambda_p = r1->Uniform(p_min, p_max);
+ Double_t lambda_phi = r1->Uniform(0.0, 2.0 * M_PI);
+ Double_t lambda_th = r1->Uniform(theta_min/1000., theta_max/1000.); // Divide by 1000 for radians
+ Double_t lambda_px = lambda_p * TMath::Cos(lambda_phi) * TMath::Sin(lambda_th);
+ Double_t lambda_py = lambda_p * TMath::Sin(lambda_phi) * TMath::Sin(lambda_th);
+ Double_t lambda_pz = lambda_p * TMath::Cos(lambda_th);
+ Double_t lambda_E = TMath::Sqrt(lambda_p*lambda_p + lambda_mass*lambda_mass);
+
+ // Rotate to lab coordinate system
+ TVector3 lambda_pvec(lambda_px, lambda_py, lambda_pz);
+ double cross_angle = -25./1000.; // in Rad
+ TVector3 pbeam_dir(TMath::Sin(cross_angle), 0, TMath::Cos(cross_angle)); //proton beam direction
+ lambda_pvec.RotateY(cross_angle); // Theta is returned positive, beam in negative X
+
+ // type 2 is state that will decay
+ GenParticlePtr p_lambda = std::make_shared<GenParticle>(
+ FourVector(lambda_pvec.X(), lambda_pvec.Y(), lambda_pvec.Z(), lambda_E), lambda_pdgID, 2 );
+
+ // Generating lambda particle, will be generated at origin
+ // Must have input electron + proton for vertex
+ GenVertexPtr lambda_initial_vertex = std::make_shared<GenVertex>();
+ lambda_initial_vertex->add_particle_in(p1);
+ lambda_initial_vertex->add_particle_in(p2);
+ lambda_initial_vertex->add_particle_out(p_lambda);
+ evt.add_vertex(lambda_initial_vertex);
+
+ // Generate neutron + pi0 in lambda rest frame
+ TLorentzVector neutron_rest, pi0_rest;
+
+ // Generating uniformly along a sphere
+ double cost_neutron_rest = r1->Uniform(-1,1);
+ double th_neutron_rest = TMath::ACos(cost_neutron_rest);
+ double sint_neutron_rest = TMath::Sin(th_neutron_rest);
+
+ double phi_neutron_rest = r1->Uniform(-1.*TMath::Pi(),1.*TMath::Pi());
+ double cosp_neutron_rest = TMath::Cos(phi_neutron_rest);
+ double sinp_neutron_rest = TMath::Sin(phi_neutron_rest);
+
+ // Calculate energy of each particle in the lambda rest frame
+ // See problem 3.19 in Introduction to Elementary Particles, 2nd edition by D. Griffiths
+ double E_neutron_rest = (-TMath::Power(pi0_mass, 2.) + TMath::Power(lambda_mass, 2.) + TMath::Power(neutron_mass, 2.) ) / (2. * lambda_mass) ;
+ double E_pi0_rest = (-TMath::Power(neutron_mass, 2.) + TMath::Power(lambda_mass, 2.) + TMath::Power(pi0_mass, 2.) ) / (2. * lambda_mass) ;
+
+ // Both particles will have the same momentum, so just use neutron variables
+ double momentum_rest = TMath::Sqrt( E_neutron_rest*E_neutron_rest - neutron_mass*neutron_mass );
+
+ neutron_rest.SetE(E_neutron_rest);
+ neutron_rest.SetPx( momentum_rest * sint_neutron_rest * cosp_neutron_rest );
+ neutron_rest.SetPy( momentum_rest * sint_neutron_rest * sinp_neutron_rest );
+ neutron_rest.SetPz( momentum_rest * cost_neutron_rest );
+
+ pi0_rest.SetE(E_pi0_rest);
+ pi0_rest.SetPx( -neutron_rest.Px() );
+ pi0_rest.SetPy( -neutron_rest.Py() );
+ pi0_rest.SetPz( -neutron_rest.Pz() );
+
+ // Boost neutron & pion to lab frame
+ TLorentzVector lambda_lab(lambda_pvec.X(), lambda_pvec.Y(), lambda_pvec.Z(), lambda_E);
+ TVector3 lambda_boost = lambda_lab.BoostVector();
+ TLorentzVector neutron_lab, pi0_lab;
+ neutron_lab = neutron_rest;
+ neutron_lab.Boost(lambda_boost);
+ pi0_lab = pi0_rest;
+ pi0_lab.Boost(lambda_boost);
+
+ // Calculating position for lambda decay
+ TVector3 lambda_unit = lambda_lab.Vect().Unit();
+ double lambda_decay_length = GetDecayLength(r1, lambda_lifetime, lambda_mass, lambda_lab.P());
+ TVector3 lambda_decay_position = lambda_unit * lambda_decay_length;
+ double lambda_decay_time = lambda_decay_length / lambda_lab.Beta() ; // Decay time in lab frame in length units (mm)
+
+ // Generating vertex for lambda decay
+ GenParticlePtr p_neutron = std::make_shared<GenParticle>(
+ FourVector(neutron_lab.Px(), neutron_lab.Py(), neutron_lab.Pz(), neutron_lab.E()), neutron_pdgID, 1 );
+
+ GenParticlePtr p_pi0 = std::make_shared<GenParticle>(
+ FourVector(pi0_lab.Px(), pi0_lab.Py(), pi0_lab.Pz(), pi0_lab.E()), pi0_pdgID, 2 );
+
+ GenVertexPtr v_lambda_decay = std::make_shared<GenVertex>(FourVector(lambda_decay_position.X(), lambda_decay_position.Y(), lambda_decay_position.Z(), lambda_decay_time));
+ v_lambda_decay->add_particle_in(p_lambda);
+ v_lambda_decay->add_particle_out(p_neutron);
+ v_lambda_decay->add_particle_out(p_pi0);
+
+ evt.add_vertex(v_lambda_decay);
+
+ // Generate two photons from pi0 decay
+ TLorentzVector gamma1_rest, gamma2_rest;
+
+ // Generating uniformly along a sphere
+ double cost_gamma1_rest = r1->Uniform(-1,1);
+ double th_gamma1_rest = TMath::ACos(cost_gamma1_rest);
+ double sint_gamma1_rest = TMath::Sin(th_gamma1_rest);
+
+ double phi_gamma1_rest = r1->Uniform(-1.*TMath::Pi(),1.*TMath::Pi());
+ double cosp_gamma1_rest = TMath::Cos(phi_gamma1_rest);
+ double sinp_gamma1_rest = TMath::Sin(phi_gamma1_rest);
+
+ // Photons are massless so they each get equal energies
+ gamma1_rest.SetE(pi0_mass/2.);
+ gamma1_rest.SetPx( (pi0_mass/2.)*sint_gamma1_rest*cosp_gamma1_rest );
+ gamma1_rest.SetPy( (pi0_mass/2.)*sint_gamma1_rest*sinp_gamma1_rest );
+ gamma1_rest.SetPz( (pi0_mass/2.)*cost_gamma1_rest );
+
+ gamma2_rest.SetE(pi0_mass/2.);
+ gamma2_rest.SetPx( -gamma1_rest.Px() );
+ gamma2_rest.SetPy( -gamma1_rest.Py() );
+ gamma2_rest.SetPz( -gamma1_rest.Pz() );
+
+ // Boost neutron & pion to lab frame
+ TVector3 pi0_boost = pi0_lab.BoostVector();
+ TLorentzVector gamma1_lab, gamma2_lab;
+ gamma1_lab = gamma1_rest;
+ gamma1_lab.Boost(pi0_boost);
+ gamma2_lab = gamma2_rest;
+ gamma2_lab.Boost(pi0_boost);
+
+ GenParticlePtr p_gamma1 = std::make_shared<GenParticle>(
+ FourVector(gamma1_lab.Px(), gamma1_lab.Py(), gamma1_lab.Pz(), gamma1_lab.E()), photon_pdgID, 1 );
+
+ GenParticlePtr p_gamma2 = std::make_shared<GenParticle>(
+ FourVector(gamma2_lab.Px(), gamma2_lab.Py(), gamma2_lab.Pz(), gamma2_lab.E()), photon_pdgID, 1 );
+
+ // Generate pi0 at same position as the lambda. Approximating pi0 decay as instantaneous
+ GenVertexPtr v_pi0_decay = std::make_shared<GenVertex>(FourVector(lambda_decay_position.X(), lambda_decay_position.Y(), lambda_decay_position.Z(), lambda_decay_time));
+ v_pi0_decay->add_particle_in(p_pi0);
+ v_pi0_decay->add_particle_out(p_gamma1);
+ v_pi0_decay->add_particle_out(p_gamma2);
+
+ //std::cout<< lambda_pvec.Angle(pbeam_dir) << " " << neutron_lab.Angle(pbeam_dir) << " " << gamma1_lab.Angle(pbeam_dir) << " " << gamma2_lab.Angle(pbeam_dir) << std::endl;
+
+ evt.add_vertex(v_pi0_decay);
+
+ if (events_parsed == 0) {
+ std::cout << "First event: " << std::endl;
+ Print::listing(evt);
+ }
+ double zdc_z=35800;
+ TVector3 extrap_neutron=lambda_decay_position+neutron_lab.Vect()*((zdc_z-pbeam_dir.Dot(lambda_decay_position))/(pbeam_dir.Dot(neutron_lab.Vect())));
+ 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)<zdc_z)
+ hepmc_output.write_event(evt);
+ if (events_parsed % 1000 == 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/lambda/analysis/lambda_plots.py b/benchmarks/lambda/analysis/lambda_plots.py
new file mode 100644
index 0000000000000000000000000000000000000000..10260be6482aaf73c1ed8b0ee20b8b0690ae9011
--- /dev/null
+++ b/benchmarks/lambda/analysis/lambda_plots.py
@@ -0,0 +1,368 @@
+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)
+
+outdir=sys.argv[1]+"/"
+try:
+ import os
+ os.mkdir(outdir[:-1])
+except:
+ pass
+import uproot as ur
+arrays_sim={}
+momenta=100, 125, 150, 175,200,225,250,275
+for p in momenta:
+ filename=f'results/lambda/epic_zdc_sipm_on_tile_only_rec_lambda_dec_{p}GeV.edm4hep.root'
+ print("opening file", filename)
+ events = ur.open(filename+':events')
+ arrays_sim[p] = events.arrays()[:-1] #remove last event, which for some reason is blank
+ import gc
+ gc.collect()
+ print("read", filename)
+
+def gauss(x, A,mu, sigma):
+ return A * np.exp(-(x-mu)**2/(2*sigma**2))
+
+#keep track of the number of clusters per event
+nclusters={}
+
+for p in momenta:
+ plt.figure()
+ nclusters[p]=[]
+ for i in range(len(arrays_sim[p])):
+ nclusters[p].append(len(arrays_sim[p]["HcalFarForwardZDCClusters.position.x"][i]))
+ nclusters[p]=np.array(nclusters[p])
+ plt.hist(nclusters[p],bins=20, range=(0,20))
+ plt.xlabel("number of clusters")
+ plt.yscale('log')
+ plt.title(f"$p_\Lambda={p}$ GeV")
+ plt.ylim(1)
+ plt.savefig(outdir+f"nclust_{p}GeV_recon.pdf")
+ print("saved file ", outdir+f"nclust_{p}GeV_recon.pdf")
+
+
+
+pt_truth={}
+theta_truth={}
+
+for p in momenta:
+ #get the truth value of theta* and pt*
+ px=arrays_sim[p]["MCParticles.momentum.x"][:,2]
+ py=arrays_sim[p]["MCParticles.momentum.y"][:,2]
+ pz=arrays_sim[p]["MCParticles.momentum.z"][:,2]
+ tilt=-0.025
+ pt_truth[p]=np.hypot(px*np.cos(tilt)-pz*np.sin(tilt), py)
+ theta_truth[p]=np.arctan2(pt_truth[p],pz*np.cos(tilt)+px*np.sin(tilt))
+
+
+#create an array with the same shape as the cluster-level arrays
+is_neutron_cand={}
+for p in momenta:
+ is_neutron_cand[p]=(0*arrays_sim[p][f"HcalFarForwardZDCClusters.energy"]).to_list()
+
+ #largest_eigenvalues
+ for i in range(len(arrays_sim[p])):
+ pars=arrays_sim[p]["_HcalFarForwardZDCClusters_shapeParameters"][i]
+ index_of_max=-1
+ max_val=0
+ eigs=[]
+ shape_params_per_cluster=7
+ for j in range(len(pars)//shape_params_per_cluster):
+ largest_eigenvalue=max(pars[shape_params_per_cluster*j+4:shape_params_per_cluster*j+7])
+ eigs.append(largest_eigenvalue)
+ if(largest_eigenvalue>max_val):
+ max_val=largest_eigenvalue
+ index_of_max=j
+ if index_of_max >=0:
+ is_neutron_cand[p][i][index_of_max]=1
+ eigs.sort()
+
+ is_neutron_cand[p]=ak.Array(is_neutron_cand[p])
+
+
+#with the position of the vertex determined by assuming the mass of the pi0
+#corrected pt* and theta* recon
+pt_recon_corr={}
+theta_recon_corr={}
+mass_recon_corr={}
+mass_pi0_recon_corr={}
+pi0_converged={}
+zvtx_recon={}
+
+maxZ=35800
+for p in momenta:
+ xvtx=0
+ yvtx=0
+ zvtx=0
+
+ for iteration in range(20):
+
+ #compute the value of theta* using the clusters in the ZDC
+ xc=arrays_sim[p][f"HcalFarForwardZDCClusters.position.x"]
+ yc=arrays_sim[p][f"HcalFarForwardZDCClusters.position.y"]
+ zc=arrays_sim[p][f"HcalFarForwardZDCClusters.position.z"]
+ E=arrays_sim[p][f"HcalFarForwardZDCClusters.energy"]
+ #apply correction to the neutron candidates only
+ A,B,C=-0.0756, -1.91, 2.30
+ neutron_corr=(1-is_neutron_cand[p])+is_neutron_cand[p]/(1+A+B/np.sqrt(E)+C/E)
+ E=E*neutron_corr
+
+ E_recon=np.sum(E, axis=-1)
+ pabs=np.sqrt(E**2-is_neutron_cand[p]*.9406**2)
+ tilt=-0.025
+ xcp=xc*np.cos(tilt)-zc*np.sin(tilt)
+ ycp=yc
+ zcp=zc*np.cos(tilt)+xc*np.sin(tilt)
+ rcp=np.sqrt(xcp**2+ycp**2+zcp**2)
+
+ ux=(xcp-xvtx)
+ uy=(ycp-yvtx)
+ uz=(zcp-zvtx)
+
+ norm=np.sqrt(ux**2+uy**2+uz**2)
+ ux=ux/norm
+ uy=uy/norm
+ uz=uz/norm
+
+ px_recon,py_recon,pz_recon=np.sum(pabs*ux, axis=-1),np.sum(pabs*uy, axis=-1),np.sum(pabs*uz, axis=-1)
+
+ pt_recon_corr[p]=np.hypot(px_recon,py_recon)
+ theta_recon_corr[p]=np.arctan2(pt_recon_corr[p], pz_recon)
+
+ mass_recon_corr[p]=np.sqrt((E_recon)**2\
+ -(px_recon)**2\
+ -(py_recon)**2\
+ -(pz_recon)**2)
+ mass_pi0_recon_corr[p]=np.sqrt(np.sum(pabs*(1-is_neutron_cand[p]), axis=-1)**2\
+ -np.sum(pabs*ux*(1-is_neutron_cand[p]), axis=-1)**2\
+ -np.sum(pabs*uy*(1-is_neutron_cand[p]), axis=-1)**2\
+ -np.sum(pabs*uz*(1-is_neutron_cand[p]), axis=-1)**2)
+ alpha=1
+ if iteration==0:
+ u=np.sqrt(px_recon**2+py_recon**2+pz_recon**2)
+ ux=px_recon/u
+ uy=py_recon/u
+ uz=pz_recon/u
+ zeta=1/2
+ zvtx=maxZ*np.power(zeta,alpha)
+ xvtx=ux/uz*zvtx
+ yvtx=uy/uz*zvtx
+ else :
+ u=np.sqrt(px_recon**2+py_recon**2+pz_recon**2)
+ ux=px_recon/u
+ uy=py_recon/u
+ uz=pz_recon/u
+ s=2*(mass_pi0_recon_corr[p]<0.135)-1
+ zeta=np.power(zvtx/maxZ, 1/alpha)
+ zeta=zeta+s*1/2**(1+iteration)
+ zvtx=maxZ*np.power(zeta,alpha)
+ xvtx=ux/uz*zvtx
+ yvtx=uy/uz*zvtx
+ #print(zvtx)
+ pi0_converged[p]=np.abs(mass_pi0_recon_corr[p]-0.135)<0.01
+ zvtx_recon[p]=zvtx
+
+fig,axs=plt.subplots(1,3, figsize=(24, 8))
+plt.sca(axs[0])
+plt.title(f"$E_{{\\Lambda}}=100-275$ GeV")
+x=[]
+y=[]
+for p in momenta:
+ accept=(nclusters[p]==3) &(pi0_converged[p])
+ x+=list(theta_truth[p][accept]*1000)
+ y+=list(theta_recon_corr[p][accept]*1000)
+plt.scatter(x,y)
+plt.xlabel("$\\theta^{*\\rm truth}_{\\Lambda}$ [mrad]")
+plt.ylabel("$\\theta^{*\\rm recon}_{\\Lambda}$ [mrad]")
+plt.xlim(0,3.2)
+plt.ylim(0,3.2)
+
+plt.sca(axs[1])
+plt.title(f"$E_{{\\Lambda}}=100-275$ GeV")
+y,x,_=plt.hist(y-np.array(x), bins=50, range=(-1,1))
+bc=(x[1:]+x[:-1])/2
+
+from scipy.optimize import curve_fit
+slc=abs(bc)<0.3
+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))
+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]")
+plt.ylabel("events")
+
+plt.sca(axs[2])
+sigmas=[]
+dsigmas=[]
+for p in momenta:
+
+ accept=(nclusters[p]==3) &(pi0_converged[p])
+ y,x=np.histogram((theta_recon_corr[p]-theta_truth[p])[accept]*1000, bins=100, range=(-1,1))
+ bc=(x[1:]+x[:-1])/2
+
+ from scipy.optimize import curve_fit
+ slc=abs(bc)<0.3
+ fnc=gauss
+ p0=(100, 0, 0.06)
+ #print(bc[slc],y[slc])
+ sigma=np.sqrt(y[slc])+(y[slc]==0)
+ coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+
+ x=np.linspace(-1, 1)
+ sigmas.append(coeff[2])
+ dsigmas.append(np.sqrt(var_matrix[2][2]))
+plt.ylim(0, 0.3)
+
+plt.errorbar(momenta, sigmas, dsigmas, ls='', marker='o', color='k')
+x=np.linspace(100, 275, 100)
+plt.plot(x, 3/np.sqrt(x), color='tab:orange')
+plt.text(170, .23, "YR requirement:\n 3 mrad/$\\sqrt{E}$")
+plt.xlabel("$E_{\\Lambda}$ [GeV]")
+plt.ylabel("$\\sigma[\\theta^*_{\\Lambda}]$ [mrad]")
+plt.tight_layout()
+plt.savefig(outdir+"thetastar_recon.pdf")
+#plt.show()
+
+
+fig,axs=plt.subplots(1,3, figsize=(24, 8))
+plt.sca(axs[0])
+plt.title(f"$E_{{\\Lambda}}=100-275$ GeV")
+x=[]
+y=[]
+for p in momenta:
+ accept=(nclusters[p]==3) &(pi0_converged[p])
+ x+=list(arrays_sim[p]['MCParticles.vertex.z'][:,3][accept]/1000)
+ y+=list(zvtx_recon[p][accept]/1000)
+plt.scatter(x,y)
+#print(x,y)
+plt.xlabel("$z^{\\rm truth}_{\\rm vtx}$ [m]")
+plt.ylabel("$z^{\\rm recon}_{\\rm vtx}$ [m]")
+plt.xlim(0,40)
+plt.ylim(0,40)
+
+plt.sca(axs[1])
+plt.title(f"$E_{{\\Lambda}}=100-275$ GeV")
+y,x,_=plt.hist(y-np.array(x), bins=50, range=(-10,10))
+bc=(x[1:]+x[:-1])/2
+
+from scipy.optimize import curve_fit
+slc=abs(bc)<5
+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))
+x=np.linspace(-5, 5)
+plt.plot(x, gauss(x, *coeff), color='tab:orange')
+print(coeff[2], np.sqrt(var_matrix[2][2]))
+plt.xlabel("$z^{*\\rm recon}_{\\rm vtx}-z^{*\\rm truth}_{\\rm vtx}$ [m]")
+plt.ylabel("events")
+
+plt.sca(axs[2])
+sigmas=[]
+dsigmas=[]
+xvals=[]
+for p in momenta:
+
+ accept=(nclusters[p]==3) &(pi0_converged[p])
+ a=list((zvtx_recon[p]-arrays_sim[p]['MCParticles.vertex.z'][:,3])[accept]/1000)
+ y,x=np.histogram(a, bins=100, range=(-10,10))
+ bc=(x[1:]+x[:-1])/2
+
+ from scipy.optimize import curve_fit
+ slc=abs(bc)<5
+ fnc=gauss
+ p0=(100, 0, 1)
+ #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))
+ sigmas.append(coeff[2])
+ dsigmas.append(np.sqrt(var_matrix[2][2]))
+ xvals.append(p)
+ except:
+ print("fit failed")
+plt.ylim(0, 2)
+
+plt.errorbar(xvals, sigmas, dsigmas, ls='', marker='o', color='k')
+x=np.linspace(100, 275, 100)
+
+avg=np.sum(sigmas/np.array(dsigmas)**2)/np.sum(1/np.array(dsigmas)**2)
+plt.axhline(avg, color='tab:orange')
+plt.text(150, 1.25,f"$\\sigma\\approx${avg:.1f} m")
+
+plt.xlabel("$E_{\\Lambda}$ [GeV]")
+plt.ylabel("$\\sigma[z_{\\rm vtx}]$ [m]")
+plt.tight_layout()
+plt.savefig(outdir+"zvtx_recon.pdf")
+#plt.show()
+
+p=100
+fig,axs=plt.subplots(1,2, figsize=(16, 8))
+plt.sca(axs[0])
+lambda_mass=1.115683
+vals=[]
+for p in momenta:
+ accept=(nclusters[p]==3) &(pi0_converged[p])
+ vals+=list(mass_recon_corr[p][accept])
+
+y,x,_= plt.hist(vals, bins=100, range=(1.0, 1.25))
+bc=(x[1:]+x[:-1])/2
+plt.axvline(lambda_mass, ls='--', color='tab:green', lw=3)
+plt.text(lambda_mass+.01, np.max(y)*1.05, "PDG mass", color='tab:green')
+plt.xlabel("$m_{\\Lambda}^{\\rm recon}$ [GeV]")
+plt.ylim(0, np.max(y)*1.2)
+plt.xlim(1.0, 1.25)
+
+from scipy.optimize import curve_fit
+slc=abs(bc-lambda_mass)<0.07
+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))
+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]))
+plt.xlabel("$m^{\\rm recon}_{\\Lambda}$ [GeV]")
+plt.ylabel("events")
+plt.title(f"$E_{{\\Lambda}}=100-275$ GeV")
+
+plt.sca(axs[1])
+xvals=[]
+sigmas=[]
+dsigmas=[]
+for p in momenta:
+ accept=(nclusters[p]==3) &(pi0_converged[p])
+ y,x= np.histogram(mass_recon_corr[p][accept], bins=100, range=(0.6,1.4))
+ bc=(x[1:]+x[:-1])/2
+
+ from scipy.optimize import curve_fit
+ slc=abs(bc-lambda_mass)<0.07
+ fnc=gauss
+ 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)))
+ x=np.linspace(0.8, 1.3, 200)
+ sigmas.append(coeff[2])
+ dsigmas.append(np.sqrt(var_matrix[2][2]))
+ xvals.append(p)
+ except:
+ print("fit failed")
+
+plt.errorbar(xvals, sigmas, dsigmas, ls='', marker='o', color='k')
+avg=np.sum(sigmas/np.array(dsigmas)**2)/np.sum(1/np.array(dsigmas)**2)
+plt.axhline(avg, color='tab:orange')
+plt.text(150, 0.01,f"$\\sigma\\approx${avg*1000:.0f} MeV")
+plt.xlabel("$E_{\\Lambda}$ [GeV]")
+plt.ylabel("$\\sigma[m_{\\Lambda}]$ [GeV]")
+plt.ylim(0, 0.02)
+plt.tight_layout()
+plt.savefig(outdir+"lambda_mass_rec.pdf")
diff --git a/benchmarks/lambda/config.yml b/benchmarks/lambda/config.yml
new file mode 100644
index 0000000000000000000000000000000000000000..1813e8d15b58ecdc6fe206ab856d549a2798796a
--- /dev/null
+++ b/benchmarks/lambda/config.yml
@@ -0,0 +1,28 @@
+
+lambda:simulate:
+ stage: simulate
+ extends: .phy_benchmark
+ needs: ["common:detector"]
+ parallel:
+ matrix:
+ - P: 100
+ - P: 125
+ - P: 150
+ - P: 175
+ - P: 200
+ - P: 225
+ - P: 250
+ - P: 275
+ timeout: 6 hours
+ script:
+ - snakemake --cores 1 results/lambda/epic_zdc_sipm_on_tile_only_rec_lambda_dec_${P}GeV.edm4hep.root
+ retry:
+ max: 2
+ when:
+ - runner_system_failure
+
+lambda:results:
+ stage: collect
+ needs: ["lambda:simulate"]
+ script:
+ - python benchmarks/lambda/analysis/lambda_plots.py results/lambda/results_epic_zdc_sipm_on_tile_only_lambda_dec
diff --git a/benchmarks/lambda/test_lambda_recon.py b/benchmarks/lambda/test_lambda_recon.py
new file mode 100644
index 0000000000000000000000000000000000000000..4b92d36e4aef3b0194e6418fcddd30f086cb198d
--- /dev/null
+++ b/benchmarks/lambda/test_lambda_recon.py
@@ -0,0 +1,32 @@
+import uproot as ur
+filename=f"lambda_recon.root"
+events = ur.open(f'{filename}:events')
+arrays = events.arrays()
+
+#get the truth value of theta*
+px=arrays_sim["MCParticles.momentum.x"][:,2]
+py=arrays_sim["MCParticles.momentum.y"][:,2]
+pz=arrays_sim["MCParticles.momentum.z"][:,2]
+tilt=-0.025
+pt=np.hypot(px*np.cos(tilt)-pz*np.sin(tilt), py)
+theta=np.arctan2(pt,pz*np.cos(tilt)+px*np.sin(tilt))
+
+#compute the value of theta* using the clusters in the ZDC
+xc=arrays["HcalFarForwardZDCClusters.position.x"]
+yc=arrays["HcalFarForwardZDCClusters.position.y"]
+zc=arrays["HcalFarForwardZDCClusters.position.z"]
+E=arrays_sim["HcalFarForwardZDCClusters.energy"]
+
+rc=np.sqrt(xc**2+yc**2+zc**2)
+xcp=xc*np.cos(tilt)-zc*np.sin(tilt)
+ycp=yc
+zcp=zc*np.cos(tilt)+xc*np.sin(tilt)
+
+E=arrays_sim["HcalFarForwardZDCClusters.energy"][i]
+
+px_recon,py_recon,pz_recon=np.sum(E*xcp/rc, axis=-1),np.sum(E*ycp/rc, axis=-1),np.sum(E*zcp/rc, axis=-1)
+pt_recon=np.hypot(px_recon,py_recon)
+theta_recon=np.arctan2(pt_recon, np.sum(E, axis=-1))
+
+
+