diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index e38469abdfac48d654d8770793167866df3159a3..49ecec331ea4c095132634b4f229289f5e33429e 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -120,9 +120,6 @@ include:
- local: 'benchmarks/dis/config.yml'
#- local: 'benchmarks/dvmp/config.yml'
- local: 'benchmarks/dvcs/config.yml'
- - local: 'benchmarks/lambda/config.yml'
- - local: 'benchmarks/neutron/config.yml'
- - local: 'benchmarks/sigma/config.yml'
- local: 'benchmarks/tcs/config.yml'
- local: 'benchmarks/u_omega/config.yml'
- local: 'benchmarks/single/config.yml'
@@ -135,9 +132,6 @@ summary:
- "demp:results"
- "dis:results"
- "dvcs:results"
- - "lambda:results"
- - "neutron:results"
- - "sigma:results"
- "tcs:results"
- "u_omega:results"
- "single:results"
diff --git a/Snakefile b/Snakefile
index 10998fe929face1020c51e955e0cd5ea5970e1a8..73b03c2b57ca94f25bb0d80d54cf8748a32d3b8f 100644
--- a/Snakefile
+++ b/Snakefile
@@ -42,7 +42,4 @@ ddsim \
include: "benchmarks/diffractive_vm/Snakefile"
include: "benchmarks/dis/Snakefile"
-include: "benchmarks/lambda/Snakefile"
-include: "benchmarks/neutron/Snakefile"
include: "benchmarks/demp/Snakefile"
-include: "benchmarks/sigma/Snakefile"
\ No newline at end of file
diff --git a/benchmarks/lambda/Snakefile b/benchmarks/lambda/Snakefile
deleted file mode 100644
index 0793621db7d0bcbbab541ec8f372e514be5658e4..0000000000000000000000000000000000000000
--- a/benchmarks/lambda/Snakefile
+++ /dev/null
@@ -1,63 +0,0 @@
-rule lambda_generate:
- input:
- script="benchmarks/lambda/analysis/gen_lambda_decay.cxx",
- params:
- NEVENTS_GEN=100000,
- output:
- GEN_FILE="sim_output/lambda/lambda_decay_{P}GeV.hepmc"
- shell:
- """
-root -l -b -q '{input.script}({params.NEVENTS_GEN},0,"{output.GEN_FILE}",{wildcards.P},{wildcards.P})'
-"""
-
-rule lambda_simulate:
- input:
- GEN_FILE="sim_output/lambda/lambda_decay_{P}GeV.hepmc"
- params:
- PHYSICS_LIST="FTFP_BERT"
- output:
- SIM_FILE="sim_output/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="sim_output/lambda/{DETECTOR_CONFIG}_sim_lambda_dec_{P}GeV.edm4hep.root"
- output:
- REC_FILE="sim_output/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_collections=MCParticles,HcalFarForwardZDCClusters,HcalFarForwardZDCRecHits,HcalFarForwardZDCSubcellHits -Pjana:nevents=$NEVENTS_REC
-"""
-
-rule lambda_analysis:
- input:
- expand("sim_output/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/{DETECTOR_CONFIG}/lambda"),
- 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
deleted file mode 100644
index 567eda5bcf5497f1b64596745ef5820e7a54eff9..0000000000000000000000000000000000000000
--- a/benchmarks/lambda/analysis/gen_lambda_decay.cxx
+++ /dev/null
@@ -1,239 +0,0 @@
-#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
deleted file mode 100644
index dffccb1c2bb7efecf266960ef997c7b812f685c9..0000000000000000000000000000000000000000
--- a/benchmarks/lambda/analysis/lambda_plots.py
+++ /dev/null
@@ -1,371 +0,0 @@
-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/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=[]
-xvals=[]
-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)
- 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, 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_{\\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
deleted file mode 100644
index 877d7016db41b2c9b0b99b7af617ae2bdbef4918..0000000000000000000000000000000000000000
--- a/benchmarks/lambda/config.yml
+++ /dev/null
@@ -1,34 +0,0 @@
-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 sim_output/lambda/epic_zdc_sipm_on_tile_only_rec_lambda_dec_${P}GeV.edm4hep.root
- retry:
- max: 2
- when:
- - runner_system_failure
-
-lambda:analyze:
- stage: analyze
- extends: .phy_benchmark
- needs: ["lambda:simulate"]
- script:
- - snakemake --cores 1 results/epic_zdc_sipm_on_tile_only/lambda
-
-lambda:results:
- stage: collect
- needs: ["lambda:analyze"]
- script:
- - ls -al
diff --git a/benchmarks/lambda/test_lambda_recon.py b/benchmarks/lambda/test_lambda_recon.py
deleted file mode 100644
index 4b92d36e4aef3b0194e6418fcddd30f086cb198d..0000000000000000000000000000000000000000
--- a/benchmarks/lambda/test_lambda_recon.py
+++ /dev/null
@@ -1,32 +0,0 @@
-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))
-
-
-
diff --git a/benchmarks/neutron/Snakefile b/benchmarks/neutron/Snakefile
deleted file mode 100644
index 5b5658f955e5f64cd84c895b271c035259a73c79..0000000000000000000000000000000000000000
--- a/benchmarks/neutron/Snakefile
+++ /dev/null
@@ -1,57 +0,0 @@
-rule neutron_generate:
- input:
- script="benchmarks/neutron/analysis/gen_particles.cxx",
- params:
- NEVENTS_GEN=1000,
- th_max=5.7,
- th_min=2.0
- output:
- GEN_FILE="sim_output/neutron/neutron_{P}GeV.hepmc"
- shell:
- """
-root -l -b -q '{input.script}({params.NEVENTS_GEN},"{output.GEN_FILE}", "neutron", {params.th_min}, {params.th_max}, 0., 360., {wildcards.P})'
-"""
-
-rule neutron_simulate:
- input:
- GEN_FILE="sim_output/neutron/neutron_{P}GeV.hepmc"
- params:
- PHYSICS_LIST="FTFP_BERT"
- output:
- SIM_FILE="sim_output/neutron/{DETECTOR_CONFIG}_sim_neutron_{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 neutron_recon:
- input:
- SIM_FILE="sim_output/neutron/{DETECTOR_CONFIG}_sim_neutron_{P}GeV.edm4hep.root"
- output:
- REC_FILE="sim_output/neutron/{DETECTOR_CONFIG}_rec_neutron_{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 -Pjana:nevents=$NEVENTS_REC
-"""
-
-rule neutron_analysis:
- input:
- expand("sim_output/neutron/{DETECTOR_CONFIG}_rec_neutron_{P}GeV.edm4hep.root",
- P=[20, 30, 40, 50, 60, 70, 80],
- DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]),
- script="benchmarks/neutron/analysis/neutron_plots.py",
- output:
- results_dir=directory("results/{DETECTOR_CONFIG}/neutron"),
- shell:
- """
-mkdir -p {output.results_dir}
-python {input.script} {output.results_dir}
-"""
diff --git a/benchmarks/neutron/analysis/gen_particles.cxx b/benchmarks/neutron/analysis/gen_particles.cxx
deleted file mode 100644
index 0877a85223bbc1af27af91e7dedfb212de7cb20b..0000000000000000000000000000000000000000
--- a/benchmarks/neutron/analysis/gen_particles.cxx
+++ /dev/null
@@ -1,127 +0,0 @@
-#include "HepMC3/GenEvent.h"
-#include "HepMC3/ReaderAscii.h"
-#include "HepMC3/WriterAscii.h"
-#include "HepMC3/Print.h"
-
-#include "TRandom3.h"
-#include "TVector3.h"
-
-#include <iostream>
-#include <random>
-#include <cmath>
-#include <math.h>
-#include <TMath.h>
-#include <TDatabasePDG.h>
-#include <TParticlePDG.h>
-
-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<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 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<GenParticle>(
- 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<GenVertex>();
- 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/neutron/analysis/neutron_plots.py b/benchmarks/neutron/analysis/neutron_plots.py
deleted file mode 100644
index 0f174c8d7a7c1d93dab439fe00c29b2ce665a343..0000000000000000000000000000000000000000
--- a/benchmarks/neutron/analysis/neutron_plots.py
+++ /dev/null
@@ -1,295 +0,0 @@
-import numpy as np, pandas as pd, matplotlib.pyplot as plt, matplotlib as mpl, awkward as ak, sys, uproot as ur
-import mplhep as hep
-hep.style.use("CMS")
-
-plt.rcParams['figure.facecolor']='white'
-plt.rcParams['savefig.facecolor']='white'
-plt.rcParams['savefig.bbox']='tight'
-
-plt.rcParams["figure.figsize"] = (7, 7)
-config=sys.argv[1].split("/")[1] #results/{config}/neutron
-outdir=sys.argv[1]+"/"
-try:
- import os
- os.mkdir(outdir[:-1])
-except:
- pass
-
-#read files
-arrays_sim={}
-for p in 20,30,40,50,60,70,80:
- arrays_sim[p] = ur.open(f'sim_output/neutron/{config}_rec_neutron_{p}GeV.edm4hep.root:events')\
- .arrays()
-
-def gauss(x, A,mu, sigma):
- return A * np.exp(-(x-mu)**2/(2*sigma**2))
-
-#get the truth pseudorapidity and energy
-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['E_truth']=np.hypot(p, 0.9406)
- array['eta_truth']=1/2*np.log((p+pzp)/(p-pzp))
- array['theta_truth']=np.arccos(pzp/p)
-
-#
-# get reconstructed theta as avg of theta of cluster centers, weighted by energy
-for array in arrays_sim.values():
- tilt=-0.025
- x=array['HcalEndcapPInsertClusters.position.x']
- y=array['HcalEndcapPInsertClusters.position.y']
- z=array['HcalEndcapPInsertClusters.position.z']
- E=array['HcalEndcapPInsertClusters.energy']
- r=np.sqrt(x**2+y**2+z**2)
-
- xp=x*np.cos(tilt)-z*np.sin(tilt)
- yp=y
- zp=z*np.cos(tilt)+x*np.sin(tilt)
-
- w=E
-
- array['theta_recon']=np.sum(np.arccos(zp/r)*w, axis=-1)/np.sum(w, axis=-1)
- array['eta_recon']=-np.log(np.tan(array['theta_recon']/2))
-
-
-#plot theta residuals:
-print("making theta recon plot")
-from scipy.optimize import curve_fit
-
-fig, axs=plt.subplots(1,2, figsize=(16,8))
-plt.sca(axs[0])
-p=40
-eta_min=3.4; eta_max=3.6
-y,x,_=plt.hist(1000*(arrays_sim[p]['theta_recon']-arrays_sim[p]['theta_truth'])\
- [(arrays_sim[p]['eta_truth']>eta_min)&(arrays_sim[p]['eta_truth']<eta_max)], bins=50,
- range=(-10,10), histtype='step')
-bc=(x[1:]+x[:-1])/2
-slc=abs(bc)<3
-# try:
-fnc=gauss
-sigma=np.sqrt(y[slc])+(y[slc]==0)
-p0=(100, 0, 5)
-coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
-xx=np.linspace(-5,5,100)
-plt.plot(xx,fnc(xx,*coeff))
-# except:
-# pass
-plt.xlabel("$\\theta_{rec}-\\theta_{truth}$ [mrad]")
-plt.ylabel("events")
-plt.title(f"$p={p}$ GeV, ${eta_min}<\\eta<{eta_max}$")
-
-r=[3.2,3.4,3.6,3.8,4.0]
-for eta_min, eta_max in zip(r[:-1],r[1:]):
- xvals=[]
- sigmas=[]
- dsigmas=[]
- for p in 20,30,40, 50, 60, 70, 80:
- y,x=np.histogram(1000*(arrays_sim[p]['theta_recon']-arrays_sim[p]['theta_truth'])\
- [(arrays_sim[p]['eta_truth']>eta_min)&(arrays_sim[p]['eta_truth']<eta_max)],
- bins=50, range=(-10,10))
- bc=(x[1:]+x[:-1])/2
- slc=abs(bc)<3
- fnc=gauss
- p0=(100, 0, 5)
- #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:
- pass
- plt.sca(axs[1])
- plt.errorbar(xvals, sigmas, dsigmas, ls='', marker='o', label=f"${eta_min}<\\eta<{eta_max}$")
- if eta_min==3.4:
- fnc=lambda E, a, b: np.hypot(a,b/np.sqrt(E))
- p0=[.002,.05]
- coeff, var_matrix = curve_fit(fnc, xvals, sigmas, p0=p0,sigma=dsigmas)
- xx=np.linspace(15, 85, 100)
- axs[1].plot(xx, fnc(xx,*coeff), color='tab:purple',ls='--',
- label=f'fit ${eta_min:.1f}<\\eta<{eta_max:.1f}$:\n'+\
- f'({coeff[0]:.2f}$\\oplus\\frac{{{coeff[1]:.1f}}}{{\\sqrt{{E}}}}$) mrad')
-plt.xlabel("$p_{n}$ [GeV]")
-plt.ylabel("$\\sigma[\\theta]$ [mrad]")
-plt.ylim(0, 5)
-plt.legend()
-plt.tight_layout()
-plt.savefig(outdir+"neutron_theta_recon.pdf")
-
-#now determine the energy recon parameters
-pvals=[]
-resvals=[]
-reserrs=[]
-wvals=[]
-svals=[]
-Euncorr=[]
-
-print("determining the energy recon correction parameters")
-fig,axs=plt.subplots(1,2, figsize=(20,10))
-eta_min=3.4;eta_max=3.6
-for p in 20, 30,40,50,60,70, 80:
- best_res=1000
- res_err=1000
- best_s=1000
- wrange=np.linspace(0.8, 1.2, 41)
- coeff_best=None
-
- wbest=0
- a=arrays_sim[p]
- h=np.sum(a[f'HcalEndcapPInsertClusters.energy'], axis=-1)
- e=np.sum(a[f'EcalEndcapPInsertClusters.energy'], axis=-1)
- for w in wrange:
-
- r=(e+h*w)[(h>0)&(a['eta_truth']>eta_min)&(a['eta_truth']<eta_max)]
- y,x=np.histogram(r,bins=50)
- bcs=(x[1:]+x[:-1])/2
- fnc=gauss
- slc=abs(bcs-np.mean(r)*1.25)<2*np.std(r)
- sigma=np.sqrt(y[slc])+0.5*(y[slc]==0)
- p0=(100, np.mean(r), np.std(r))
- try:
- coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
- res=np.abs(coeff[2]/coeff[1])
-
- if res<best_res:
- best_res=res
- coeff_best=coeff
- res_err=res*np.hypot(np.sqrt(var_matrix[2][2])/coeff[2], np.sqrt(var_matrix[1][1])/coeff[1])
- wbest=w
- best_s=np.hypot(p,0.9406)/coeff[1]
- Euncorr_best=coeff[1]
- except :
- print("fit failed")
-
- if p==50:
- r=(e+h*wbest)[(h>0)&(a['eta_truth']>3.4)&(a['eta_truth']<3.6)]
- plt.sca(axs[0])
- y, x, _= plt.hist(r, histtype='step', bins=50)
- xx=np.linspace(20, 55, 100)
- plt.plot(xx,fnc(xx, *coeff_best), ls='-')
- plt.xlabel("$E_{uncorr}=w\\times E_{Hcal}+E_{Ecal}$ [GeV]")
- plt.title(f"p=50 GeV, ${eta_min}<\\eta<{eta_max}$, w={wbest:.2f}")
- plt.axvline(np.sqrt(50**2+.9406**2), color='g', ls=':')
- plt.text(40, max(y)*0.9, "generated\nenergy", color='g', fontsize=20)
-
- Euncorr.append(Euncorr_best)
- resvals.append(best_res)
- reserrs.append(res_err)
- pvals.append(p)
- svals.append(best_s)
- wvals.append(wbest)
-
-pvals=np.array(pvals)
-svals=np.array(svals)
-Euncorr=np.array(Euncorr)
-plt.sca(axs[1])
-plt.plot(Euncorr,wvals, label="w")
-w_avg=np.mean(wvals)
-plt.axhline(w_avg, label=f'w avg: {w_avg:.2f}', ls=':')
-plt.plot(Euncorr,svals, label="s")
-m=(np.sum(svals*Euncorr)*len(Euncorr)-np.sum(Euncorr)*np.sum(svals))/(np.sum(Euncorr**2)*len(Euncorr)-np.sum(Euncorr)**2)
-b=np.mean(svals)-np.mean(Euncorr)*m
-plt.plot(Euncorr,Euncorr*m+b, label=f"s fit: ${m:.4f}E_{{uncorr}}+{b:.2f}$", ls=':')
-
-plt.xlabel("$E_{uncorr}=w\\times E_{Hcal}+E_{Ecal}$ [GeV]")
-plt.title("$E_{n,recon}=s\\times(w\\times E_{Hcal}+E_{Ecal})$")
-plt.ylabel('parameter values')
-plt.legend()
-plt.ylim(0)
-plt.savefig(outdir+"neutron_energy_params.pdf")
-
-#now make the energy plot
-print("making energy recon plot")
-fig, axs=plt.subplots(1,3, figsize=(24,8))
-partitions=[3.2,3.4, 3.6, 3.8, 4.0]
-
-for eta_min, eta_max in zip(partitions[:-1],partitions[1:]):
- pvals=[]
- resvals=[]
- reserrs=[]
- scalevals=[]
- scaleerrs=[]
- for p in 20, 30,40,50,60,70, 80:
- best_res=1000
- res_err=1000
-
- wrange=np.linspace(30, 70, 30)*0.0257
-
- w=w_avg
- a=arrays_sim[p]
- h=np.sum(a[f'HcalEndcapPInsertClusters.energy'], axis=-1)
- e=np.sum(a[f'EcalEndcapPInsertClusters.energy'], axis=-1)
- #phi=a['phi_truth']
- uncorr=(e+h*w)
- s=-0.0047*uncorr+1.64
- r=uncorr*s #reconstructed energy with correction
- r=r[(h>0)&(a['eta_truth']>eta_min)&(a['eta_truth']<eta_max)]#&(abs(phi)>np.pi/2)]
- y,x=np.histogram(r,bins=50)
- bcs=(x[1:]+x[:-1])/2
- fnc=gauss
- slc=abs(bcs-np.mean(r)*1.25)<2*np.std(r)
- sigma=np.sqrt(y[slc])+0.5*(y[slc]==0)
- p0=(100, np.mean(r), np.std(r))
- try:
- coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
- res=np.abs(coeff[2]/coeff[1])
-
- if res<best_res:
- best_res=res
- res_err=res*np.hypot(np.sqrt(var_matrix[2][2])/coeff[2], np.sqrt(var_matrix[1][1])/coeff[1])
- wbest=w
- scale=coeff[1]/np.sqrt(p**2+0.9406**2)
- dscale=np.sqrt(var_matrix[1][1]/np.sqrt(p**2+0.9406**2))
- except :
- print("fit failed")
- if p==50 and eta_min==3.4:
- plt.sca(axs[0])
- plt.errorbar(bcs, y, np.sqrt(y)+(y==0),marker='o', ls='', )
- plt.title(f'p={p} GeV, ${eta_min}<\\eta<{eta_max}$')
- plt.xlabel("$E_{recon}$ [GeV]")
- plt.ylabel("events")
- #plt.ylim(0)
- xx=np.linspace(40, 70, 50)
- plt.plot(xx, fnc(xx, *coeff))
- resvals.append(best_res)
- reserrs.append(res_err)
- scalevals.append(scale)
- scaleerrs.append(dscale)
- pvals.append(p)
- plt.sca(axs[1])
- plt.errorbar(pvals, resvals, reserrs, marker='o', ls='', label=f"${eta_min}<\\eta<{eta_max}$")
- #plt.ylim(0)
- plt.ylabel("$\\sigma[E]/\\mu[E]$")
- plt.xlabel("$p_{n}$ [GeV]")
-
- plt.sca(axs[2])
- plt.errorbar(pvals, scalevals, scaleerrs, marker='o', ls='', label=f"${eta_min}<\\eta<{eta_max}$")
-
-
- plt.ylabel("$\\mu[E]/E$")
- if eta_min==3.4:
- fnc=lambda E, b: b/np.sqrt(E)
- p0=[.5]
- coeff, var_matrix = curve_fit(fnc, pvals, resvals, p0=p0,sigma=np.array(reserrs))
- xx=np.linspace(15, 85, 100)
- axs[1].plot(xx, fnc(xx,*coeff), color='tab:purple',ls='--',
- label=f'fit ${eta_min:.1f}<\\eta<{eta_max:.1f}$: '+\
- f'$\\frac{{{coeff[0]*100:.0f}\\%}}{{\\sqrt{{E}}}}$')
-axs[2].set_xlabel("$p_n$ [GeV]")
-axs[2].axhline(1, ls='--', color='0.5', alpha=0.7)
-axs[0].set_ylim(0)
-axs[1].set_ylim(0, 0.35)
-axs[2].set_ylim(0)
-axs[1].legend(fontsize=20)
-axs[2].legend(fontsize=20)
-plt.tight_layout()
-plt.savefig(outdir+"neutron_energy_recon.pdf")
diff --git a/benchmarks/neutron/config.yml b/benchmarks/neutron/config.yml
deleted file mode 100644
index 304b8311a47eac6ed00782815936edb5b047a9d9..0000000000000000000000000000000000000000
--- a/benchmarks/neutron/config.yml
+++ /dev/null
@@ -1,34 +0,0 @@
-neutron:simulate:
- stage: simulate
- extends: .phy_benchmark
- needs: ["common:detector"]
- 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/neutron/epic_craterlake_rec_neutron_${P}GeV.edm4hep.root
- retry:
- max: 2
- when:
- - runner_system_failure
-
-neutron:analyze:
- stage: analyze
- extends: .phy_benchmark
- needs: ["neutron:simulate"]
- script:
- - mkdir -p results/epic_craterlake
- - python benchmarks/neutron/analysis/neutron_plots.py results/epic_craterlake/neutron
-
-neutron:results:
- stage: collect
- needs: ["neutron:analyze"]
- script:
- - ls -al
diff --git a/benchmarks/sigma/Snakefile b/benchmarks/sigma/Snakefile
deleted file mode 100644
index 77b64df8eb4c7c5b35c34be38b7184e58b25c99e..0000000000000000000000000000000000000000
--- a/benchmarks/sigma/Snakefile
+++ /dev/null
@@ -1,63 +0,0 @@
-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:
- """
-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}
-"""
diff --git a/benchmarks/sigma/analysis/gen_sigma_decay.cxx b/benchmarks/sigma/analysis/gen_sigma_decay.cxx
deleted file mode 100644
index 4da92801257e861bca80cc6638722fff360019c7..0000000000000000000000000000000000000000
--- a/benchmarks/sigma/analysis/gen_sigma_decay.cxx
+++ /dev/null
@@ -1,305 +0,0 @@
-#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;
-}
diff --git a/benchmarks/sigma/analysis/sigma_plots.py b/benchmarks/sigma/analysis/sigma_plots.py
deleted file mode 100644
index 1dbf47e2d50fb848f13b619e83051f23dee9fa67..0000000000000000000000000000000000000000
--- a/benchmarks/sigma/analysis/sigma_plots.py
+++ /dev/null
@@ -1,480 +0,0 @@
-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")
diff --git a/benchmarks/sigma/benchmark.json b/benchmarks/sigma/benchmark.json
deleted file mode 100644
index 7c475ddbaf393ea27bec5fbf4003953746b546af..0000000000000000000000000000000000000000
--- a/benchmarks/sigma/benchmark.json
+++ /dev/null
@@ -1,6 +0,0 @@
-{
- "name": "Sigma0 in ZDC",
- "title": "Sigma0 Benchmark",
- "description": "Benchmark for measuring sigma0 reconstruction using its decay products in the ZDC",
- "target": "0.8"
-}
diff --git a/benchmarks/sigma/config.yml b/benchmarks/sigma/config.yml
deleted file mode 100644
index 46845c37ce945ce5b0b85fd76212dada83bec2c7..0000000000000000000000000000000000000000
--- a/benchmarks/sigma/config.yml
+++ /dev/null
@@ -1,34 +0,0 @@
-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