Sigma zdc (#20)

* added sigma generator script

* require lambda decay to occur before ZDC when generating sigma0 events

* Update gen_sigma_decay.cxx

* added Snakefile; benchmakr.json

* added sigma to global snakefile

* updated snake file to use the results directory

* updated sigma_plots.py

* use sim_output instead of results for large files

* Update benchmarks/sigma/config.yml

Co-authored-by: Wouter Deconinck <wdconinc@gmail.com>

* Update sigma_plots.py

* Update .gitlab-ci.yml

* Update .gitlab-ci.yml

* Update sigma_plots.py

* Update config.yml

increased time limit

* Update Snakefile

removed /   (not sure if this will change anything

* removed hardcoded config in sigma_plots.py

* added a space at the end of one line.  idk if this fixes anything

---------

Co-authored-by: Wouter Deconinck <wdconinc@gmail.com>

.gitlab-ci.yml

@@ -120,6 +120,7 @@ include:
     #- local: 'benchmarks/dvmp/config.yml'
   - local: 'benchmarks/dvcs/config.yml'
   - local: 'benchmarks/lambda/config.yml'
+  - local: 'benchmarks/sigma/config.yml'
   - local: 'benchmarks/tcs/config.yml'
   - local: 'benchmarks/u_omega/config.yml'
   - local: 'benchmarks/single/config.yml'
@@ -133,6 +134,7 @@ summary:
     - "dis:results"
     - "dvcs:results"
     - "lambda:results"
+    - "sigma:results"
     - "tcs:results"
     - "u_omega:results"
     - "single:results"

Snakefile

@@ -44,3 +44,4 @@ include: "benchmarks/diffractive_vm/Snakefile"
 include: "benchmarks/dis/Snakefile"
 include: "benchmarks/lambda/Snakefile"
 include: "benchmarks/demp/Snakefile"
+include: "benchmarks/sigma/Snakefile"
\ No newline at end of file

benchmarks/sigma/Snakefile

+rule sigma_generate:
+	input:
+                script="benchmarks/sigma/analysis/gen_sigma_decay.cxx",
+	params:
+		NEVENTS_GEN=100000,
+	output:
+		GEN_FILE="sim_output/sigma/sigma_decay_{P}GeV.hepmc"
+	shell:
+		"""
+mkdir -p results/sigma
+root -l -b -q '{input.script}({params.NEVENTS_GEN},0,"{output.GEN_FILE}",{wildcards.P},{wildcards.P})'
+"""
+
+rule sigma_simulate:
+	input:
+		GEN_FILE="sim_output/sigma/sigma_decay_{P}GeV.hepmc"
+	params:
+		PHYSICS_LIST="FTFP_BERT"
+	output:
+		SIM_FILE="sim_output/sigma/{DETECTOR_CONFIG}_sim_sigma_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 sigma_recon:
+        input:
+                SIM_FILE="sim_output/sigma/{DETECTOR_CONFIG}_sim_sigma_dec_{P}GeV.edm4hep.root"
+        output:
+                REC_FILE="sim_output/sigma/{DETECTOR_CONFIG}_rec_sigma_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_collections=MCParticles,HcalFarForwardZDCClusters,HcalFarForwardZDCRecHits,HcalFarForwardZDCSubcellHits  -Pjana:nevents=$NEVENTS_REC
+"""
+
+rule sigma_analysis:
+	input:
+                expand("sim_output/sigma/{DETECTOR_CONFIG}_rec_sigma_dec_{P}GeV.edm4hep.root",
+		    P=[100, 125, 150,175, 200, 225, 250, 275],
+		    DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]),
+                script="benchmarks/sigma/analysis/sigma_plots.py",
+	output:
+		results_dir=directory("results/{DETECTOR_CONFIG}/sigma"),
+	shell:
+		"""
+mkdir -p {output.results_dir}
+python {input.script} {output.results_dir}
+"""

benchmarks/sigma/analysis/gen_sigma_decay.cxx

+#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 sigma baryons and decay them to a neutron + 2 photons
+void gen_sigma_decay(int n_events = 100000, UInt_t seed = 0, char* out_fname = "sigma_decay.hepmc",
+		      double p_min = 100., // in GeV/c
+		      double p_max = 275.) // in GeV/c
+{
+  int accepted_events=0;
+  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 sigma_info = GetParticleInfo(pdg, "Sigma0");
+  double sigma_mass = std::get<0>(sigma_info);
+  int sigma_pdgID = std::get<1>(sigma_info);
+  double sigma_lifetime = std::get<2>(sigma_info);
+  
+  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 sigma_p     = r1->Uniform(p_min, p_max);
+    Double_t sigma_phi   = r1->Uniform(0.0, 2.0 * M_PI);
+    Double_t sigma_th    = r1->Uniform(theta_min/1000., theta_max/1000.); // Divide by 1000 for radians
+    Double_t sigma_px    = sigma_p * TMath::Cos(sigma_phi) * TMath::Sin(sigma_th);
+    Double_t sigma_py    = sigma_p * TMath::Sin(sigma_phi) * TMath::Sin(sigma_th);
+    Double_t sigma_pz    = sigma_p * TMath::Cos(sigma_th);
+    Double_t sigma_E     = TMath::Sqrt(sigma_p*sigma_p + sigma_mass*sigma_mass);
+    
+    
+    TVector3 sigma_pvec(sigma_px, sigma_py, sigma_pz);
+    
+    double cross_angle = -25./1000.; // in Rad
+    TVector3 pbeam_dir(TMath::Sin(cross_angle), 0, TMath::Cos(cross_angle)); //proton beam direction
+    sigma_pvec.RotateY(cross_angle); // Theta is returned positive, beam in negative X
+
+    // type 2 is state that will decay
+    GenParticlePtr p_sigma = std::make_shared<GenParticle>(
+        FourVector(sigma_pvec.X(), sigma_pvec.Y(), sigma_pvec.Z(), sigma_E), sigma_pdgID, 2 );
+    // Generating sigma particle, will be generated at origin
+        // Must have input electron + proton for vertex
+        GenVertexPtr sigma_initial_vertex = std::make_shared<GenVertex>();
+        sigma_initial_vertex->add_particle_in(p1);
+        sigma_initial_vertex->add_particle_in(p2);
+        sigma_initial_vertex->add_particle_out(p_sigma);
+        evt.add_vertex(sigma_initial_vertex);
+
+        // Generate lambda + gamma in sigma rest frame
+        TLorentzVector lambda_rest, gamma_rest;
+
+        // Generating uniformly along a sphere
+        double cost_lambda_rest = r1->Uniform(-1,1);
+        double th_lambda_rest = TMath::ACos(cost_lambda_rest);
+        double sint_lambda_rest = TMath::Sin(th_lambda_rest);
+
+        double phi_lambda_rest = r1->Uniform(-1.*TMath::Pi(),1.*TMath::Pi());
+        double cosp_lambda_rest = TMath::Cos(phi_lambda_rest);
+        double sinp_lambda_rest = TMath::Sin(phi_lambda_rest);
+
+        // Calculate energy of each particle in the sigma rest frame
+        // See problem 3.19 in Introduction to Elementary Particles, 2nd edition by D. Griffiths
+        double E_lambda_rest = (-TMath::Power(photon_mass, 2.) + TMath::Power(sigma_mass, 2.) + TMath::Power(lambda_mass, 2.) ) / (2. * sigma_mass) ;
+        double E_gamma_rest = (-TMath::Power(lambda_mass, 2.) + TMath::Power(sigma_mass, 2.) + TMath::Power(photon_mass, 2.) ) / (2. * sigma_mass) ;
+
+        // Both particles will have the same momentum, so just use lambda variables
+        double momentum_rest = TMath::Sqrt( E_lambda_rest*E_lambda_rest - lambda_mass*lambda_mass );
+
+        lambda_rest.SetE(E_lambda_rest);
+        lambda_rest.SetPx( momentum_rest * sint_lambda_rest * cosp_lambda_rest );
+        lambda_rest.SetPy( momentum_rest * sint_lambda_rest * sinp_lambda_rest );
+        lambda_rest.SetPz( momentum_rest * cost_lambda_rest );
+
+        gamma_rest.SetE(E_gamma_rest);
+        gamma_rest.SetPx( -lambda_rest.Px() );
+        gamma_rest.SetPy( -lambda_rest.Py() );
+        gamma_rest.SetPz( -lambda_rest.Pz() );
+
+        // Boost lambda & pion to lab frame
+        TLorentzVector sigma_lab(sigma_pvec.X(), sigma_pvec.Y(), sigma_pvec.Z(), sigma_E);
+        TVector3 sigma_boost = sigma_lab.BoostVector();
+        TLorentzVector lambda_lab, gamma_lab;
+        lambda_lab = lambda_rest;
+        lambda_lab.Boost(sigma_boost);
+        gamma_lab = gamma_rest;
+        gamma_lab.Boost(sigma_boost);
+        
+    // Calculating position for sigma decay
+    TVector3 sigma_unit = sigma_lab.Vect().Unit();
+    double sigma_decay_length = GetDecayLength(r1, sigma_lifetime, sigma_mass, sigma_lab.P());
+    TVector3 sigma_decay_position = sigma_unit * sigma_decay_length;
+    double sigma_decay_time = sigma_decay_length / sigma_lab.Beta() ; // Decay time in lab frame in length units (mm)
+  
+    // Generating vertex for sigma decay
+    GenParticlePtr p_lambda = std::make_shared<GenParticle>(
+      FourVector(lambda_lab.Px(), lambda_lab.Py(), lambda_lab.Pz(), lambda_lab.E()), lambda_pdgID, 2 );
+
+    GenParticlePtr p_gamma = std::make_shared<GenParticle>(
+      FourVector(gamma_lab.Px(), gamma_lab.Py(), gamma_lab.Pz(), gamma_lab.E()), photon_pdgID, 1 );
+
+    GenVertexPtr v_sigma_decay = std::make_shared<GenVertex>(FourVector(sigma_decay_position.X(), sigma_decay_position.Y(), sigma_decay_position.Z(), sigma_decay_time));
+    v_sigma_decay->add_particle_in(p_sigma);
+    v_sigma_decay->add_particle_out(p_lambda);
+    v_sigma_decay->add_particle_out(p_gamma);
+
+    evt.add_vertex(v_sigma_decay);
+
+    // 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
+    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
+    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_gamma=sigma_decay_position+gamma_lab.Vect()*((zdc_z-pbeam_dir.Dot(sigma_decay_position))/(pbeam_dir.Dot(gamma_lab.Vect())));
+    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 && extrap_gamma.Angle(pbeam_dir)<0.004 && lambda_decay_position.Dot(pbeam_dir)<zdc_z){
+      hepmc_output.write_event(evt);
+      accepted_events ++;
+    }
+    if (events_parsed % 1000 == 0) {
+      std::cout << "Event: " << events_parsed << " ("<<accepted_events<<" accepted)"<< std::endl;
+    }
+    evt.clear();
+  }
+  hepmc_output.close();
+
+  std::cout << "Events generated: " << events_parsed << " ("<<accepted_events<<" accepted)"<< std::endl;
+}

benchmarks/sigma/analysis/sigma_plots.py

+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'sim_output/sigma/epic_zdc_sipm_on_tile_only_rec_sigma_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_\Sigma={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=[]
+        #Must make sure this doesn't get messed up if someone changes the number of shape parameters in EICrecon.
+        nClust=nclusters[p][i]
+        nShapePars=len(pars)//nClust
+        for j in range(nClust):
+            largest_eigenvalue=max(pars[nShapePars*j+4:nShapePars*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])
+
+import ROOT
+
+lambda_mass=1.115683
+pi0_mass=0.1349768
+pt_recon_corr={}
+theta_recon_corr={}
+mass_recon_corr={}
+mass_lambda_recon_corr={}
+mass_pi0_recon_corr={}
+pi0_converged={}
+zvtx_recon={}
+
+#run this event-by-event:
+maxZ=35800
+for p in momenta:
+    pt_recon_corr[p]=[]
+    theta_recon_corr[p]=[]
+    mass_recon_corr[p]=[]
+    mass_lambda_recon_corr[p]=[]
+    mass_pi0_recon_corr[p]=[]
+    zvtx_recon[p]=[]
+    for evt in range(len(arrays_sim[p])):
+        if nclusters[p][evt]!=4:
+            nan=-1
+            pt_recon_corr[p].append(nan)
+            theta_recon_corr[p].append(nan)
+            mass_recon_corr[p].append(nan)
+            mass_lambda_recon_corr[p].append(nan)
+            mass_pi0_recon_corr[p].append(nan)
+            zvtx_recon[p].append(nan)
+            continue
+        xc=arrays_sim[p][f"HcalFarForwardZDCClusters.position.x"][evt]
+        yc=arrays_sim[p][f"HcalFarForwardZDCClusters.position.y"][evt]
+        zc=arrays_sim[p][f"HcalFarForwardZDCClusters.position.z"][evt]
+        E=arrays_sim[p][f"HcalFarForwardZDCClusters.energy"][evt]
+        
+        #apply correction to the neutron candidates only
+        A,B,C=-0.0756, -1.91,  2.30
+        neutron_corr=(1-is_neutron_cand[p][evt])+is_neutron_cand[p][evt]/(1+A+B/np.sqrt(E)+C/E)
+        E=E*neutron_corr
+
+        pabs=np.sqrt(E**2-is_neutron_cand[p][evt]*.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)
+        
+        #search for the combination of photons that would give the best lambda mass
+        pt_best=-999
+        theta_best=-999
+        mass_lambda_best=-999
+        mass_sigma_best=-999
+        mass_pi0_best=-999
+        zvtx_best=-999
+        for hypothesis in range(4):
+            if is_neutron_cand[p][evt][hypothesis]:
+                continue
+            
+            xvtx=0
+            yvtx=0
+            zvtx=0
+            #find the vertex position that reconstructs the pi0 mass
+            for iteration in range(20):
+                tot=ROOT.TLorentzVector(0,0,0,0)
+                Lambda=ROOT.TLorentzVector(0,0,0,0)
+                pi0=ROOT.TLorentzVector(0,0,0,0)
+
+                for i in range(4):
+
+                    if i!=hypothesis:
+                        ux=xcp[i]-xvtx
+                        uy=ycp[i]-yvtx
+                        uz=zcp[i]-zvtx
+                    else:
+                        ux=xcp[i]
+                        uy=ycp[i]
+                        uz=zcp[i]
+                    u=np.sqrt(ux**2+uy**2+uz**2)
+                    ux/=u
+                    uy/=u
+                    uz/=u
+
+                    P=ROOT.TLorentzVector(pabs[i]*ux, pabs[i]*uy, pabs[i]*uz, E[i])
+                    tot+=P
+                    if not is_neutron_cand[p][evt][i] and i!=hypothesis:
+                        pi0+=P
+                    if i!=hypothesis:
+                        Lambda+=P
+                alpha=1
+                if iteration==0:
+                    zeta=1/2
+                    zvtx=maxZ*np.power(zeta,alpha)
+                    xvtx=Lambda.X()/Lambda.Z()*zvtx
+                    yvtx=Lambda.Y()/Lambda.Z()*zvtx
+                else :
+                    s=2*(pi0.M()<pi0_mass)-1
+                    zeta=np.power(zvtx/maxZ, 1/alpha)
+                    zeta=zeta+s*1/2**(1+iteration)
+                    zvtx=maxZ*np.power(zeta,alpha)
+                    xvtx=Lambda.X()/Lambda.Z()*zvtx
+                    yvtx=Lambda.Y()/Lambda.Z()*zvtx
+
+            if abs(Lambda.M()-lambda_mass)< abs(mass_lambda_best-lambda_mass):
+                pt_best=tot.Pt()
+                theta_best=tot.Theta()
+                mass_lambda_best=Lambda.M()
+                mass_sigma_best=tot.M()
+                mass_pi0_best=pi0.M()
+                zvtx_best=zvtx
+                
+        pt_recon_corr[p].append(pt_best)
+        theta_recon_corr[p].append(theta_best)
+        mass_recon_corr[p].append(mass_sigma_best)
+        mass_lambda_recon_corr[p].append(mass_lambda_best)
+        mass_pi0_recon_corr[p].append(mass_pi0_best)
+        zvtx_recon[p].append(zvtx_best)
+    pt_recon_corr[p]=ak.Array(pt_recon_corr[p])
+    theta_recon_corr[p]=ak.Array(theta_recon_corr[p])
+    mass_recon_corr[p]=ak.Array(mass_recon_corr[p])
+    mass_lambda_recon_corr[p]=ak.Array(mass_lambda_recon_corr[p])
+    mass_pi0_recon_corr[p]=ak.Array(mass_pi0_recon_corr[p])
+    zvtx_recon[p]=ak.Array(zvtx_recon[p])
+        
+#now make plots
+
+#reconstructed vertex position plot
+fig,axs=plt.subplots(1,3, figsize=(24, 8))
+plt.sca(axs[0])
+plt.title(f"$E_{{\\Sigma}}=100-275$ GeV")
+x=[]
+y=[]
+for p in momenta:
+    accept=(nclusters[p]==4)# &(pi0_converged[p])
+    x+=list(theta_truth[p][accept]*1000)
+    y+=list(theta_recon_corr[p][accept]*1000)
+#print(x)
+plt.scatter(x,y)
+plt.xlabel("$\\theta^{*\\rm truth}_{\\Sigma}$ [mrad]")
+plt.ylabel("$\\theta^{*\\rm recon}_{\\Sigma}$ [mrad]")
+plt.xlim(0,3.2)
+plt.ylim(0,3.2)
+
+plt.sca(axs[1])
+plt.title(f"$E_{{\\Sigma}}=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.6
+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))
+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]")
+plt.ylabel("events")
+
+plt.sca(axs[2])
+sigmas=[]
+dsigmas=[]
+xvals=[]
+for p in momenta:
+    
+    accept=(nclusters[p]==4)
+    y,x=np.histogram((theta_recon_corr[p]-theta_truth[p])[accept]*1000, bins=100, range=(-0.5,0.5))
+    bc=(x[1:]+x[:-1])/2
+
+    from scipy.optimize import curve_fit
+    slc=abs(bc)<0.3
+    fnc=gauss
+    p0=(100, 0, 0.15)
+    #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(np.abs(coeff[2]))
+        dsigmas.append(np.sqrt(var_matrix[2][2]))
+        xvals.append(p)
+    except:
+        print(f"fit failed for p={p}")
+print(xvals)
+print(sigmas)
+plt.ylim(0, 0.3)
+
+plt.errorbar(xvals, 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_{\\Sigma}$ [GeV]")
+plt.ylabel("$\\sigma[\\theta^*_{\\Sigma}]$ [mrad]")
+plt.tight_layout()
+plt.savefig(outdir+"thetastar_recon.pdf")
+
+#reconstructed vertex position plot
+fig,axs=plt.subplots(1,3, figsize=(24, 8))
+plt.sca(axs[0])
+plt.title(f"$E_{{\\Sigma}}=100-275$ GeV")
+x=[]
+y=[]
+for p in momenta:
+    accept=(nclusters[p]==4)&(abs(mass_pi0_recon_corr[p]-pi0_mass)<.01)
+    tilt=-0.025
+    x+=list(arrays_sim[p]['MCParticles.vertex.z'][:,5][accept]*np.cos(tilt)/1000
+            +np.sin(tilt)*arrays_sim[p]['MCParticles.vertex.z'][:,5][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_{{\\Sigma}}=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')
+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]==4)&(abs(mass_pi0_recon_corr[p]-pi0_mass)<.01)
+    a=list((zvtx_recon[p]-arrays_sim[p]['MCParticles.vertex.z'][:,5])[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(abs(coeff[2]))
+        dsigmas.append(np.sqrt(var_matrix[2][2]))
+        xvals.append(p)
+    except:
+        print(f"fit failed for p={p}")
+plt.ylim(0, 3)
+
+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_{\\Sigma}$ [GeV]")
+plt.ylabel("$\\sigma[z_{\\rm vtx}]$ [m]")
+plt.tight_layout()
+plt.savefig(outdir+"zvtx_recon.pdf")
+        
+#lambda mass reconstruction
+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]==4)&(abs(mass_pi0_recon_corr[p]-pi0_mass)<.01)
+    vals+=list(mass_lambda_recon_corr[p][accept])
+
+y,x,_= plt.hist(vals, bins=100, range=(0.9, 1.3))
+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(0.9, 1.3)
+
+from scipy.optimize import curve_fit
+slc=abs(bc-lambda_mass)<0.05
+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))
+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_{{\\Sigma}}=100-275$ GeV")
+
+plt.sca(axs[1])
+xvals=[]
+sigmas=[]
+dsigmas=[]
+for p in momenta:
+    accept=(nclusters[p]==4)&(abs(mass_pi0_recon_corr[p]-pi0_mass)<.01)
+    y,x= np.histogram(mass_lambda_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.05
+    fnc=gauss
+    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)))
+        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_{\\Sigma}$ [GeV]")
+plt.ylabel("$\\sigma[m_{\\Lambda}]$ [GeV]")
+plt.ylim(0, 0.02)
+plt.tight_layout()
+plt.savefig(outdir+"lambda_mass_rec_from_sigma_decay.pdf")
+
+#sigma mass reconstruction
+p=100
+fig,axs=plt.subplots(1,2, figsize=(16, 8))
+plt.sca(axs[0])
+sigma_mass=1.192
+vals=[]
+for p in momenta:
+    accept=(nclusters[p]==4)&(abs(mass_pi0_recon_corr[p]-pi0_mass)<.01)
+    vals+=list(mass_recon_corr[p][accept])
+
+y,x,_= plt.hist(vals, bins=100, range=(1.0, 1.4))
+bc=(x[1:]+x[:-1])/2
+plt.axvline(sigma_mass, ls='--', color='tab:green', lw=3)
+plt.text(sigma_mass+.01, np.max(y)*1.05, "PDG mass", color='tab:green')
+plt.xlabel("$m_{\\Sigma}^{\\rm recon}$ [GeV]")
+plt.ylim(0, np.max(y)*1.2)
+plt.xlim(1.0, 1.45)
+
+from scipy.optimize import curve_fit
+slc=abs(bc-sigma_mass)<0.02
+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))
+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}_{\\Sigma}$ [GeV]")
+plt.ylabel("events")
+plt.title(f"$E_{{\\Sigma}}=100-275$ GeV")
+
+plt.sca(axs[1])
+xvals=[]
+sigmas=[]
+dsigmas=[]
+for p in momenta:
+    accept=(nclusters[p]==4)&(abs(mass_pi0_recon_corr[p]-pi0_mass)<.01)
+    y,x= np.histogram(mass_recon_corr[p][accept], bins=100, range=(1.0,1.4))
+    bc=(x[1:]+x[:-1])/2
+
+    from scipy.optimize import curve_fit
+    slc=abs(bc-sigma_mass)<0.02
+    fnc=gauss
+    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)))
+        sigmas.append(abs(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_{\\Sigma}$ [GeV]")
+plt.ylabel("$\\sigma[m_{\\Sigma}]$ [GeV]")
+plt.ylim(0, 0.1)
+plt.tight_layout()
+plt.savefig(outdir+"sigma_mass_rec.pdf")

benchmarks/sigma/benchmark.json

+{
+  "name": "Sigma0 in ZDC",
+  "title": "Sigma0 Benchmark",
+    "description": "Benchmark for measuring sigma0 reconstruction using its decay products in the ZDC",
+  "target": "0.8"
+}

benchmarks/sigma/config.yml

+sigma: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 sim_output/sigma/epic_zdc_sipm_on_tile_only_rec_sigma_dec_${P}GeV.edm4hep.root 
+  retry:
+    max: 2
+    when:
+      - runner_system_failure
+
+sigma:analyze:
+  stage: analyze
+  extends: .phy_benchmark
+  needs: ["sigma:simulate"]
+  script:
+    - snakemake --cores 1 results/epic_zdc_sipm_on_tile_only/sigma 
+
+sigma:results:
+  stage: collect
+  needs: ["sigma:analyze"]
+  script:
+    - ls -al