Skip to content
Snippets Groups Projects
Commit 8bba950e authored by Sebouh Paul's avatar Sebouh Paul Committed by Whitney Armstrong
Browse files

added ZDC SiPM-on-tile reconstruction benchmark

parent 173a92de
Branches
No related tags found
1 merge request!311added ZDC SiPM-on-tile reconstruction benchmark
Hide whitespace changes
Inline Side-by-side
Showing
with 1259 additions and 0 deletions
benchmarks/far_forward/analysis/analysis_zdc_sipmontile_particles.cxx 0 → 100644
+ 307
− 0
View file @ 8bba950e
////////////////////////////////////////
// Read ROOT output file
// Plot variables
// Jihee Kim 09/2021
////////////////////////////////////////
#include "ROOT/RDataFrame.hxx"
#include <iostream>
#include <vector>
#include "edm4hep/MCParticleCollection.h"
#include "edm4hep/SimCalorimeterHitCollection.h"
#include "edm4eic/ClusterCollection.h"
#include "edm4eic/ClusterData.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TGaxis.h"
#include "TPaveStats.h"
#include "TMath.h"
#include "TH1.h"
#include "TF1.h"
#include "TH1D.h"
#include "TVector3.h"
R__LOAD_LIBRARY(libedm4eic.so)
#include "fmt/core.h"
#include "DD4hep/Detector.h"
#include "DDG4/Geant4Data.h"
#include "DDRec/CellIDPositionConverter.h"
using ROOT::RDataFrame;
void analysis_zdc_sipmontile_particles(
const std::string& input_fname = "sim_zdc_sipmontile_uniform.edm4hep.root",
const std::string& results_path = "results/far_forward/zdc_sipmontile/"
) {
// Setting for graphs
TStyle* kStyle = new TStyle("kStyle","Kim's Style");
kStyle->SetOptStat("emr");
kStyle->SetOptTitle(0);
kStyle->SetOptFit(1101);
kStyle->SetStatColor(0);
kStyle->SetStatW(0.15);
kStyle->SetStatH(0.10);
kStyle->SetStatX(0.9);
kStyle->SetStatY(0.9);
kStyle->SetStatBorderSize(1);
kStyle->SetLabelSize(0.045,"xyz");
kStyle->SetTitleSize(0.045,"xyz");
kStyle->SetTitleOffset(1.2,"y");
kStyle->SetLineWidth(2);
kStyle->SetTitleFont(42,"xyz");
kStyle->SetLabelFont(42,"xyz");
kStyle->SetCanvasDefW(700);
kStyle->SetCanvasDefH(500);
kStyle->SetCanvasColor(0);
kStyle->SetPadTickX(1);
kStyle->SetPadTickY(1);
kStyle->SetPadGridX(1);
kStyle->SetPadGridY(1);
kStyle->SetPadLeftMargin(0.1);
kStyle->SetPadRightMargin(0.1);
kStyle->SetPadTopMargin(0.1);
kStyle->SetPadBottomMargin(0.1);
TGaxis::SetMaxDigits(3);
gROOT->SetStyle("kStyle");
ROOT::EnableImplicitMT();
ROOT::RDataFrame d0("events", input_fname);
// Thrown Energy [GeV]
auto E_gen = [](const std::vector<edm4hep::MCParticleData> & input) {
auto p = input[2];
auto energy = TMath::Sqrt(p.momentum.x * p.momentum.x + p.momentum.y * p.momentum.y + p.momentum.z * p.momentum.z + p.mass * p.mass);
return energy;
};
// Theta [rad]
auto Theta_gen = [](const std::vector<edm4hep::MCParticleData> & input) {
auto p = input[2];
auto theta = std::atan2(std::hypot(p.momentum.x, p.momentum.y), p.momentum.z);
return theta;
};
auto Thetarec = [](const std::vector<edm4eic::ClusterData> & input) {
std::vector<float> output;
for (const auto &c: input) {
auto r = c.position;
auto theta = std::atan2(std::hypot(r.x, r.y), r.z);
output.push_back(theta);
}
return output;
};
// Phi [rad]
auto Phi_gen = [](const std::vector<edm4hep::MCParticleData> & input) {
auto p = input[2];
auto phi = std::atan2(p.momentum.y, p.momentum.x);
if (phi < 0) phi += 2.*M_PI;
return phi;
};
auto Phirec = [](const std::vector<edm4eic::ClusterData> & input) {
std::vector<float> output;
for (const auto &c: input) {
auto r = c.position;
auto phi = std::atan2(r.y, r.x);
if (phi < 0) phi += 2.*M_PI;
output.push_back(phi);
}
return output;
};
// Eta
auto Eta_gen = [](const std::vector<edm4hep::MCParticleData> & input) {
auto p = input[2];
auto theta = std::atan2(std::hypot(p.momentum.x, p.momentum.y), p.momentum.z);
auto eta = -std::log(std::tan(theta/2));
return eta;
};
auto Etarec = [](const std::vector<edm4eic::ClusterData> & input) {
std::vector<float> output;
for (const auto &c: input) {
auto r = c.position;
auto theta = std::atan2(std::hypot(r.x, r.y), r.z);
auto eta = -std::log(std::tan(theta/2));
output.push_back(eta);
}
return output;
};
// Define variables
auto d1 = d0
// Truth
.Define("E_gen", E_gen, {"MCParticles"})
.Define("Theta_gen", Theta_gen, {"MCParticles"})
.Define("Phi_gen", Phi_gen, {"MCParticles"})
.Define("Eta_gen", Eta_gen, {"MCParticles"})
// Ecal
/*.Define("E_Ecal", {"ZDCEcalClusters.energy"})
.Define("Theta_Ecal", Thetarec, {"ZDCEcalClusters"})
.Define("Phi_Ecal", Phirec, {"ZDCEcalClusters"})
.Define("Eta_Ecal", Etarec, {"ZDCEcalClusters"})*/
// HCal
/*.Define("E_Hcal", {"HcalFarForwardZDCClusters.energy"})
.Define("Theta_Hcal", Thetarec, {"HcalFarForwardZDCClusters"})
.Define("Phi_Hcal", Phirec, {"HcalFarForwardZDCClusters"})
.Define("Eta_Hcal", Etarec, {"HcalFarForwardZDCClusters"})*/
;
// Histogram boundaries
const auto E_nbins = 100;
const auto E_max = 50.0;
const auto theta_nbins = 50;
const auto theta_max = 0.050;
const auto phi_nbins = 90;
const auto phi_min = 0.0;
const auto phi_max = 2.*M_PI;
const auto eta_nbins = 40;
const auto eta_min = 6.0;
const auto eta_max = 10.0;
// Define Histograms
auto hE_gen = d1.Histo1D({"hE_gen", "Thrown Energy; E_{gen} [GeV]; Events", E_nbins, 0.0, E_max}, "E_gen");
auto hTheta_gen = d1.Histo1D({"hTheta_gen", "Thrown Theta; #theta_{gen} [rad]; Events", theta_nbins, 0.0, theta_max}, "Theta_gen");
auto hPhi_gen = d1.Histo1D({"hPhi_gen", "Thrown Phi; #phi_{gen} [rad]; Events", phi_nbins, phi_min, phi_max}, "Phi_gen");
auto hEta_gen = d1.Histo1D({"hEta_gen", "Thrown Eta; #eta_{gen}; Events", eta_nbins, eta_min, eta_max}, "Eta_gen");
auto hPhiTheta_gen = d1.Histo2D({"hPhiTheta_gen", "Thrown Phi vs Theta; #theta_{gen} [rad]; #phi_{gen}", theta_nbins, 0.0, theta_max, phi_nbins, phi_min, phi_max},
"Theta_gen", "Phi_gen");
/* auto hE_Ecal = d1.Histo1D({"hE_Ecal", "Reconstructed Ecal Energy; E_{rec} [GeV]; Events", E_nbins, 0.0, E_max}, "E_Ecal");
auto hTheta_Ecal = d1.Histo1D({"hTheta_Ecal", "Reconstructed Ecal Theta; #theta_{rec} [rad]; Events", theta_nbins, 0.0, theta_max}, "Theta_Ecal");
auto hPhi_Ecal = d1.Histo1D({"hPhi_Ecal", "Reconstructed Ecal Phi; #phi_{rec} [rad]; Events", phi_nbins, phi_min, phi_max}, "Phi_Ecal");
auto hEta_Ecal = d1.Histo1D({"hEta_Ecal", "Reconstructed Ecal Eta; #eta_{rec}; Events", eta_nbins, eta_min, eta_max}, "Eta_Ecal");
auto hPhiTheta_Ecal = d1.Histo2D({"hPhiTheta_Ecal", "Reconstructed Ecal Phi vs Theta; #theta_{rec} [rad]; #phi_{rec}", theta_nbins, 0.0, theta_max, phi_nbins, phi_min, phi_max},
"Theta_Ecal", "Phi_Ecal");*/
/*
auto hE_Hcal = d1.Histo1D({"hE_Hcal", "Reconstructed Hcal Energy; E_{rec} [GeV]; Events", E_nbins, 0.0, E_max}, "E_Hcal");
auto hTheta_Hcal = d1.Histo1D({"hTheta_Hcal", "Reconstructed Hcal Theta; #theta_{rec} [rad]; Events", theta_nbins, 0.0, theta_max}, "Theta_Hcal");
auto hPhi_Hcal = d1.Histo1D({"hPhi_Hcal", "Reconstructed Hcal Phi; #phi_{rec} [rad]; Events", phi_nbins, phi_min, phi_max}, "Phi_Hcal");
auto hEta_Hcal = d1.Histo1D({"hEta_Hcal", "Reconstructed Hcal Eta; #eta_{rec}; Events", eta_nbins, eta_min, eta_max}, "Eta_Hcal");
auto hPhiTheta_Hcal = d1.Histo2D({"hPhiTheta_Hcal", "Reconstructed Hcal Phi vs Theta; #theta_{rec} [rad]; #phi_{rec}", theta_nbins, 0.0, theta_max, phi_nbins, phi_min, phi_max},
"Theta_Hcal", "Phi_Hcal");
*/
// Event Counts
auto nevents_thrown = d1.Count();
std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n";
// Draw Histograms
{
TCanvas* c1 = new TCanvas("c1", "c1");
auto h = hE_gen->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c1->SaveAs(TString(results_path + "/zdc_E_gen.png"));
c1->SaveAs(TString(results_path + "/zdc_E_gen.pdf"));
}
{
TCanvas* c2 = new TCanvas("c2", "c2");
auto h = hTheta_gen->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c2->SaveAs(TString(results_path + "/zdc_Theta_gen.png"));
c2->SaveAs(TString(results_path + "/zdc_Theta_gen.pdf"));
}
{
TCanvas* c3 = new TCanvas("c3", "c3");
auto h = hPhi_gen->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c3->SaveAs(TString(results_path + "/zdc_Phi_gen.png"));
c3->SaveAs(TString(results_path + "/zdc_Phi_gen.pdf"));
}
{
TCanvas* c4 = new TCanvas("c4", "c4");
auto h = hEta_gen->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c4->SaveAs(TString(results_path + "/zdc_Eta_gen.png"));
c4->SaveAs(TString(results_path + "/zdc_Eta_gen.pdf"));
}
{
TCanvas* c5 = new TCanvas("c5", "c5");
auto h = hPhiTheta_gen->DrawCopy("COLZ");
c5->SaveAs(TString(results_path + "/zdc_PhiTheta_gen.png"));
c5->SaveAs(TString(results_path + "/zdc_PhiTheta_gen.pdf"));
}
/*
// Draw Histograms Ecal
{
TCanvas* c1 = new TCanvas("c1", "c1");
auto h = hE_Ecal->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c1->SaveAs(TString(results_path + "/zdc_E_Ecal.png"));
c1->SaveAs(TString(results_path + "/zdc_E_Ecal.pdf"));
}
{
TCanvas* c2 = new TCanvas("c2", "c2");
auto h = hTheta_Ecal->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c2->SaveAs(TString(results_path + "/zdc_Theta_Ecal.png"));
c2->SaveAs(TString(results_path + "/zdc_Theta_Ecal.pdf"));
}
{
TCanvas* c3 = new TCanvas("c3", "c3");
auto h = hPhi_Ecal->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c3->SaveAs(TString(results_path + "/zdc_Phi_Ecal.png"));
c3->SaveAs(TString(results_path + "/zdc_Phi_Ecal.pdf"));
}
{
TCanvas* c4 = new TCanvas("c4", "c4");
auto h = hEta_Ecal->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c4->SaveAs(TString(results_path + "/zdc_Eta_Ecal.png"));
c4->SaveAs(TString(results_path + "/zdc_Eta_Ecal.pdf"));
}
{
TCanvas* c5 = new TCanvas("c5", "c5");
auto h = hPhiTheta_Ecal->DrawCopy("COLZ");
c5->SaveAs(TString(results_path + "/zdc_PhiTheta_Ecal.png"));
c5->SaveAs(TString(results_path + "/zdc_PhiTheta_Ecal.pdf"));
}
*/
/*
// Draw Histograms Hcal
{
TCanvas* c1 = new TCanvas("c1", "c1");
auto h = hE_Hcal->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c1->SaveAs(TString(results_path + "/zdc_E_Hcal.png"));
c1->SaveAs(TString(results_path + "/zdc_E_Hcal.pdf"));
}
{
TCanvas* c2 = new TCanvas("c2", "c2");
auto h = hTheta_Hcal->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c2->SaveAs(TString(results_path + "/zdc_Theta_Hcal.png"));
c2->SaveAs(TString(results_path + "/zdc_Theta_Hcal.pdf"));
}
{
TCanvas* c3 = new TCanvas("c3", "c3");
auto h = hPhi_Hcal->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c3->SaveAs(TString(results_path + "/zdc_Phi_Hcal.png"));
c3->SaveAs(TString(results_path + "/zdc_Phi_Hcal.pdf"));
}
{
TCanvas* c4 = new TCanvas("c4", "c4");
auto h = hEta_Hcal->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c4->SaveAs(TString(results_path + "/zdc_Eta_Hcal.png"));
c4->SaveAs(TString(results_path + "/zdc_Eta_Hcal.pdf"));
}
{
TCanvas* c5 = new TCanvas("c5", "c5");
auto h = hPhiTheta_Hcal->DrawCopy("COLZ");
c5->SaveAs(TString(results_path + "/zdc_PhiTheta_Hcal.png"));
c5->SaveAs(TString(results_path + "/zdc_PhiTheta_Hcal.pdf"));
}*/
}
benchmarks/far_forward/analysis/gen_zdc_sipmontile_particles.cxx 0 → 100644
+ 97
− 0
View file @ 8bba950e
//////////////////////////////////////////////////////////////
// ZDC detector
// Single Particle dataset
// J. Kim 09/2021
//////////////////////////////////////////////////////////////
#include "HepMC3/GenEvent.h"
#include "HepMC3/Print.h"
#include "HepMC3/ReaderAscii.h"
#include "HepMC3/WriterAscii.h"
#include <TMath.h>
#include <cmath>
#include <iostream>
#include <math.h>
#include <random>
#include <TRandom.h>
#include <Math/Vector3D.h>
#include <Math/RotationY.h>
using namespace HepMC3;
void gen_zdc_sipmontile_particles(int n_events = 1e6, std::string particle = "neutron", double p_start = 20.0, double p_end = 140.0, const std::string& out_fname = "zdc_neutron.hepmc") {
WriterAscii hepmc_output(out_fname);
int events_parsed = 0;
GenEvent evt(Units::GEV, Units::MM);
// Random number generator
TRandom* r1 = new TRandom();
// Set crossing angle [rad]
double crossing_angle = -0.025;
// Set polar angle range [rad]
double theta_min = 0.0;
double theta_max = 0.005;
for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
// FourVector(px,py,pz,e,pdgid,status)
// 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
Double_t p_mag = r1->Uniform(p_start, p_end);
Double_t phi = r1->Uniform(0.0, 2.0 * M_PI);
Double_t theta = acos(r1->Uniform(cos(theta_max), cos(theta_min)));
Double_t p0_x = p_mag * std::cos(phi) * std::sin(theta);
Double_t p0_y = p_mag * std::sin(phi) * std::sin(theta);
Double_t p0_z = p_mag * std::cos(theta);
// Rotate the vector in Y by crossing angle when particles are being generated
ROOT::Math::XYZVector p0{p0_x, p0_y, p0_z};
ROOT::Math::RotationY r(crossing_angle);
auto p = r * p0;
auto px = p.X();
auto py = p.Y();
auto pz = p.Z();
// FourVector(px,py,pz,e,pdgid,status)
// status - type 1 is final state
// pdgid 22 - photon massless
// pdgid 2112 - neutron 939.565 MeV/c^2
GenParticlePtr p3;
if (particle == "neutron") {
p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, std::hypot(p_mag, 0.939565)), 2112, 1);
} else if (particle == "photon") {
p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, p_mag), 22, 1);
} else if (particle == "muon") {
p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, std::hypot(p_mag, 0.1057)), 13, 1);
} else {
std::cout << "Particle type " << particle << " not recognized!" << std::endl;
exit(-1);
}
// Create vertex
GenVertexPtr v1 = std::make_shared<GenVertex>();
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;
}
benchmarks/far_forward/analysis/make_plots_zdc_sipmontile_hexplit.py 0 → 100644
+ 29
− 0
View file @ 8bba950e
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import sys
infile=sys.argv[1]
outdir=sys.argv[2]
df=pd.read_csv(infile)
print(df.columns)
#first the px over pz
plt.hist(df.eval("(x_reco_rw/z_reco_rw-x_truth/z_truth)*1000"), bins=np.linspace(-2,2, 20), histtype='step', label='HEXPLIT')
plt.hist(df.eval("(x_reco/z_reco-x_truth/z_truth)*1000"), bins=np.linspace(-2,2, 20), histtype='step', label='baseline')
plt.xlabel("$x_{\\rm reco}/z_{\\rm reco}-x_{\\rm truth}/z_{\\rm reco}$ [mrad]")
plt.ylabel("events")
plt.legend()
print("saving to file,",outdir+"/residuals_x.png")
plt.savefig(outdir+"/residuals_x.png")
plt.figure()
plt.hist(df.eval("(y_reco_rw/z_reco_rw-y_truth/z_truth)*1000"), bins=np.linspace(-2, 2, 20), histtype='step', label='HEXPLIT')
plt.hist(df.eval("(y_reco/z_reco-y_truth/z_truth)*1000"), bins=np.linspace(-2, 2, 20), histtype='step', label='baseline')
plt.xlabel("$y_{\\rm reco}/z_{\\rm reco}-y_{\\rm truth}/z_{\\rm reco}$ [mrad]")
plt.ylabel("events")
plt.legend()
print("saving to file,",outdir+"/residuals_y.png")
plt.savefig(outdir+"/residuals_y.png")
print([np.std(df[a][abs(df[a])<70]) for a in "dx dx_rw dy dy_rw".split()])
benchmarks/far_forward/config.yml
+ 15
− 0
View file @ 8bba950e
......@@ -17,12 +17,27 @@ ZDC_far_forward:
stage: run
needs: ["far_forward:compile"]
script:
- print_env.sh
- bash benchmarks/far_forward/run_zdc.sh --particle $PARTICLE
parallel:
matrix:
- PARTICLE: ["neutron", "photon"]
allow_failure: true
ZDC_sipmontile:
extends: .rec_benchmark
variables:
DETECTOR_CONFIG: epic_zdc_sipm_on_tile_only
stage: run
needs: ["far_forward:compile"]
script:
- print_env.sh
- bash benchmarks/far_forward/run_zdc_sipmontile.sh --particle $PARTICLE
parallel:
matrix:
- PARTICLE: ["neutron", "photon"]
allow_failure: true
far_forward:collect:
stage: collect
needs:
......
benchmarks/far_forward/hexplit/hexplit.py 0 → 100644
+ 566
− 0
View file @ 8bba950e
import uproot as ur, numpy as np, pandas as pd
c_in_mm_per_ns = 299.792458
#default values. Can be configured
Emin=0.000472*0.1 # 1/10 a MIP
tmax=150 # 150 ns
def get_xyzr_reco_no_reweighting(arrays, event, w0=4, weight_by_granularity=True, prefix="ZDC", MIP=0.00047):
x=arrays[f'{prefix}HitsReco.position.x'][event]
y=arrays[f'{prefix}HitsReco.position.y'][event]
z=arrays[f'{prefix}HitsReco.position.z'][event]
E=arrays[f'{prefix}HitsReco.energy'][event]
t=arrays[f'{prefix}HitsReco.time'][event] - z/c_in_mm_per_ns #correct for time of flight
sl=arrays[f'{prefix}HitsReco.dimension.x'][event]/2
#make selection of hits to be used
slc=(E>Emin) & (t<tmax)
x=x[slc]
y=y[slc]
z=z[slc]
E=E[slc]
sl=sl[slc]
if type(w0)!=float or w0 !=0 :
w=w0+np.log((E+.0000001)/sum(E))
w=w*(w>0)
else :
w=E
if weight_by_granularity:
w=w/sl**2
sumw=np.sum(w+.0000001, axis=-1)
x_reco=np.sum(x*w, axis=-1)/sumw
y_reco=np.sum(y*w, axis=-1)/sumw
z_reco=np.sum(z*w, axis=-1)/sumw
r_reco=np.hypot(x_reco,y_reco)
return [x_reco,y_reco,z_reco,r_reco]
def get_xyzr_truth(arrays, event, w0=4, weight_by_granularity=True, prefix="ZDC"):
#first determine the z value used in the recon
z=arrays[f'{prefix}HitsReco.position.z'][event]
E=arrays[f'{prefix}HitsReco.energy'][event]
t=arrays[f'{prefix}HitsReco.time'][event] - z/c_in_mm_per_ns #correct for time of flight
sl=arrays[f'{prefix}HitsReco.dimension.x'][event]
#make selection
slc=(E>Emin) & (t<tmax)
z=z[slc]
E=E[slc]
t=t[slc]
sl=sl[slc]
if type(w0)!= float or w0 !=0:
w=w0+np.log((E+.0000001)/sum(E))
w=w*(w>0)
else :
w=E
if weight_by_granularity:
w=w/sl**2
z_reco=np.sum(z*w, axis=-1)/np.sum(w+.000000001, axis=-1)
px=arrays["MCParticles.momentum.x"][event,2]
py=arrays["MCParticles.momentum.y"][event,2]
pz=arrays["MCParticles.momentum.z"][event,2]
x_truth=px/pz*z_reco
y_truth=py/pz*z_reco
z_truth=z_reco
r_truth=np.hypot(x_truth,y_truth)
return [x_truth,y_truth,z_truth,r_truth]
sqrt3=np.sqrt(3)
sqrt5=np.sqrt(5)
def mx(a,b):
return a*(a>b)+b*(b>=a)
#with H3
def get_xyzr_reco_reweighted_H3(arrays, event, w0=5, weight_by_granularity=True, prefix="ZDC", MIP=0.000472):
x=arrays[f'{prefix}HitsReco.position.x'][event]
y=arrays[f'{prefix}HitsReco.position.y'][event]
z=arrays[f'{prefix}HitsReco.position.z'][event]
E=arrays[f'{prefix}HitsReco.energy'][event]
#lay=arrays['HcalEndcapPInsertHitsReco.layer'][event]
sl=arrays[f'{prefix}HitsReco.dimension.x'][event]/2
t=arrays[f'{prefix}HitsReco.time'][event] - z/c_in_mm_per_ns #correct for time of flight
slc=(E>Emin) & (t<tmax)
x=x[slc]
y=y[slc]
z=z[slc]
E=E[slc]
sl=sl[slc]
minz=min(z)
dz=min(z[z!=minz])-minz
Etot=sum(E)
if type(w0) !=float or w0 != 0.0:
thresh=Etot*np.exp(-np.max(w0))
else :
thresh = 0
#print("thresh", thresh)
xnew=[]
ynew=[]
znew=[]
Enew=[]
slnew=[]
phi=np.linspace(0, np.pi*5/3, 6)
cph=np.cos(phi)
sph=np.sin(phi)
roll_cph=np.roll(cph,1)
roll_sph=np.roll(sph,1)
#centers of the triangles minus that of the hexagon
dx_tri=(cph+roll_cph)/3
dy_tri=(sph+roll_sph)/3
for i in range(len(x)):
if E[i]<thresh:
continue
neighbors_found=0
Eneighbors=[0,0,0,0,0,0]
for j in range(len(x)):
if abs(z[i]-z[j])>dz*1.1 or E[j]<thresh or z[i]==z[j]:
continue
dx=(x[j]-x[i])/sl[i]
dy=(y[j]-y[i])/sl[i]
if abs(dx)>1.1 or abs(dy)>1.1:
continue
tol=0.01
if neighbors_found==6:
break
for k in range(6):
if Eneighbors[k]:
continue
if abs(dx-cph[k])<tol and abs(dy-sph[k])<tol:
#print("found neighbor")
Eneighbors[k]=E[j]
#print(Eneighbors, E[j])
neighbors_found+=1
break
#a=thresh*b
Eneighbors=np.array(Eneighbors)
reweight_energy=mx(Eneighbors,MIP)*mx(np.roll(Eneighbors,1),MIP)
reweight_energy/=sum(reweight_energy)
for k in range(6):
if E[i]*reweight_energy[k]<thresh:
continue
xnew.append(x[i]+sl[i]*dx_tri[k])
ynew.append(y[i]+sl[i]*dy_tri[k])
znew.append(z[i])
Enew.append(E[i]*reweight_energy[k])
slnew.append(sl[i]/sqrt5)
xnew=np.array(xnew)
ynew=np.array(ynew)
znew=np.array(znew)
Enew=np.array(Enew)
slnew=np.array(slnew)
if type(w0)!=float or w0 !=0:
w=w0+np.log((Enew+.0000001)/Etot)
w=w*(w>0)
else :
w=Enew
if weight_by_granularity:
w=w/slnew**2
sumw=np.sum(w+.0000001, axis=-1)
x_reco=np.sum(xnew*w, axis=-1)/sumw
y_reco=np.sum(ynew*w, axis=-1)/sumw
z_reco=np.sum(znew*w, axis=-1)/sumw
r_reco=np.hypot(x_reco,y_reco)
return [x_reco,y_reco,z_reco,r_reco]
#with H4
def get_xyzr_reco_reweighted_H4(arrays, event, w0=6, weight_by_granularity=True, prefix="ZDC", MIP=0.000472):
x=arrays[f'{prefix}HitsReco.position.x'][event]
y=arrays[f'{prefix}HitsReco.position.y'][event]
z=arrays[f'{prefix}HitsReco.position.z'][event]
E=arrays[f'{prefix}HitsReco.energy'][event]
t=arrays[f'{prefix}HitsReco.time'][event] - z/c_in_mm_per_ns #correct for time of flight
sl=arrays[f'{prefix}HitsReco.dimension.x'][event]/2
lx=arrays[f'{prefix}HitsReco.local.x'][event]
ly=arrays[f'{prefix}HitsReco.local.y'][event]
lz=arrays[f'{prefix}HitsReco.local.z'][event]
# time and energy cuts
slc=(E>Emin) & (t<tmax)
x=x[slc]
y=y[slc]
z=z[slc]
E=E[slc]
sl=sl[slc]
lx=lx[slc]
ly=ly[slc]
lz=lz[slc]
tol=0.01
#determine the spacing between layers
if len(lz)==0:
return None
minz=min(lz)
dz=min(lz[lz>minz+tol])-minz
Etot=sum(E)
if type(w0) !=float or w0 != 0.0:
thresh=Etot*np.exp(-np.max(w0))
else :
thresh = 0
xnew=[]
ynew=[]
znew=[]
Enew=[]
slnew=[]
phi=np.linspace(0, np.pi*5/3, 6)
cph=np.cos(phi)
sph=np.sin(phi)
for i in range(len(x)):
if E[i]<thresh:
continue
# there are twelve neighboring cell positions where we need to determine the energy
Eneighbors=[0,0,0,0,0,0,0,0,0,0,0,0]
for j in range(len(x)):
if abs(lz[i]-lz[j])>dz*2.1 or E[j]<thresh or j == i :
continue
dx=(lx[j]-lx[i])/sl[i]
dy=(ly[j]-ly[i])/sl[i]
if abs(dx)>1.6 or abs(dy)>1.6:
continue
tol=0.01
for k in range(6):
if abs(dx-1.5*cph[k])<tol and abs(dy-1.5*sph[k])<tol:
Eneighbors[k]+=E[j]
break
for k in range(6):
if abs(dx+sqrt3/2*sph[k])<tol and abs(dy-sqrt3/2*cph[k])<tol:
Eneighbors[k+6]+=E[j]
break
#print(Eneighbors)
Eneighbors=mx(np.array(Eneighbors),MIP)
reweight_energy_1=Eneighbors[:6]*np.roll(Eneighbors[6:],1)*np.roll(Eneighbors[6:],2)
reweight_energy_2=Eneighbors[6:]*np.roll(Eneighbors[6:],-1)*np.roll(Eneighbors[6:],1)
reweight_energy=np.concatenate([reweight_energy_1, reweight_energy_2])
reweight_energy/=sum(reweight_energy)
for k in range(6):
if E[i]*reweight_energy[k]< thresh:
continue
xnew.append(x[i]+sl[i]*0.75*cph[k])
ynew.append(y[i]+sl[i]*0.75*sph[k])
znew.append(z[i])
Enew.append(E[i]*reweight_energy[k])
slnew.append(sl[i]/sqrt5)
for k in range(6):
if E[i]*reweight_energy[k+6]<thresh:
continue
xnew.append(x[i]+sl[i]*sqrt3/4*-sph[k])
ynew.append(y[i]+sl[i]*sqrt3/4*cph[k])
znew.append(z[i])
Enew.append(E[i]*reweight_energy[k+6])
slnew.append(sl[i]/sqrt5)
xnew=np.array(xnew)
ynew=np.array(ynew)
znew=np.array(znew)
Enew=np.array(Enew)
slnew=np.array(slnew)
if type(w0)!=float or w0 !=0:
w=w0+np.log((Enew+.0000001)/Etot)
w=w*(w>0)
else :
w=Enew
if weight_by_granularity:
w=w/slnew**2
sumw=np.sum(w+.0000001, axis=-1)
x_reco=np.sum(xnew*w, axis=-1)/sumw
y_reco=np.sum(ynew*w, axis=-1)/sumw
z_reco=np.sum(znew*w, axis=-1)/sumw
r_reco=np.hypot(x_reco,y_reco)
return [x_reco,y_reco,z_reco,r_reco]
sqrt2=np.sqrt(2)
#with S2
def get_xyzr_reco_reweighted_S2(arrays, event, w0=6, weight_by_granularity=True, prefix="ZDC", MIP=0.000472):
x=arrays[f'{prefix}HitsReco.position.x'][event]
y=arrays[f'{prefix}HitsReco.position.y'][event]
z=arrays[f'{prefix}HitsReco.position.z'][event]
E=arrays[f'{prefix}HitsReco.energy'][event]
#side length.
sl=arrays[f'{prefix}HitsReco.dimension.x'][event]
t=arrays[f'{prefix}HitsReco.time'][event] - z/c_in_mm_per_ns #correct for time of flight
slc=(E>Emin) & (t<tmax)
x=x[slc]
y=y[slc]
z=z[slc]
E=E[slc]
sl=sl[slc]
minz=min(z)
dz=min(z[z!=minz])-minz
Etot=sum(E)
if type(w0) !=float or w0 != 0.0:
thresh=Etot*np.exp(-np.max(w0))
else :
thresh = 0
xnew=[]
ynew=[]
znew=[]
Enew=[]
slnew=[]
phi=np.linspace(np.pi/4, np.pi*7/4, 4)
cph=np.cos(phi)
sph=np.sin(phi)
#print("bbb")
for i in range(len(x)):
if E[i]<thresh:
continue
neighbors_found = 0
Eneighbors=[0,0,0,0]
for j in range(len(x)):
if abs(z[i]-z[j])>dz*1.1 or E[j]<thresh or j == i or z[i]==z[j] :
continue
dx=(x[j]-x[i])/sl[i]
dy=(y[j]-y[i])/sl[i]
if abs(dx)>0.6 or abs(dy)>0.6:
continue
tol=0.01
for k in range(4):
if abs(dx-sqrt2*cph[k]/2)<tol and abs(dy-sqrt2*sph[k]/2)<tol:
Eneighbors[k]+=E[j]
neighbors_found+=1
break
if neighbors_found == 8:
break
#print("??")
Eneighbors=mx(np.array(Eneighbors),MIP)
#print("Eneighbors:", Eneighbors)
reweight_energy=Eneighbors/sum(Eneighbors)
#print("reweight_energy:", reweight_energy)
#print(reweight_energy)
#print(Eneighbors)
#reweight_energy/=sum(reweight_energy)
#print("aaa")
for k in range(4):
if E[i]*reweight_energy[k]< thresh:
continue
xnew.append(x[i]+sl[i]*sqrt2/4*cph[k])
ynew.append(y[i]+sl[i]*sqrt2/4*sph[k])
znew.append(z[i])
Enew.append(E[i]*reweight_energy[k])
slnew.append(sl[i]/2)
xnew=np.array(xnew)
ynew=np.array(ynew)
znew=np.array(znew)
Enew=np.array(Enew)
slnew=np.array(slnew)
#print("fff")
if type(w0)!=float or w0 !=0:
w=w0+np.log((Enew+.0000001)/Etot)
w=w*(w>0)
else :
w=Enew
if weight_by_granularity:
w=w/slnew**2
#print("ggg")
sumw=np.sum(w+.0000001, axis=-1)
x_reco=np.sum(xnew*w, axis=-1)/sumw
y_reco=np.sum(ynew*w, axis=-1)/sumw
z_reco=np.sum(znew*w, axis=-1)/sumw
r_reco=np.hypot(x_reco,y_reco)
#print("hh")
return [x_reco,y_reco,z_reco,r_reco]
if __name__ == "__main__":
import argparse
parser = argparse.ArgumentParser(
prog='process_showers',
description='determines the position of the showers',
epilog='Text at the bottom of help')
parser.add_argument('infile')
parser.add_argument('outfile')
parser.add_argument('--H4', help="flag for using H4 HEX-SPLIT",
action='store_true') # on/off flag
parser.add_argument('--H3', help="flag for using H3 HEX-SPLIT",
action='store_true') # on/off flag
parser.add_argument('--S2', help="flag for using S2 HEX-SPLIT",
action='store_true') # on/off flag
parser.add_argument('-n', '--nevents', help="number of events to run", default=-1, type=int)
parser.add_argument("-p", '--prefix', help="prefix for the detector type (hit type)", default="ZDC")
parser.add_argument("-s", '--skip', help="number of events to skip", default=0, type=int)
parser.add_argument("--w0_rw", help="w0 used in reweighted", default=5, type=float)
parser.add_argument("--w0_nrw", help="w0 used in reweighted", default=4, type=float)
parser.add_argument("--w0_use_range", help="use a range for w0", action="store_true")
parser.add_argument("--MIP", help="MIP value in GeV", default=0.000472, type=float)
args = parser.parse_args()
import sys
infile=args.infile
outfile=args.outfile
useH3reweighting=args.H3
useH4reweighting=args.H4
useS2reweighting=args.S2
nevents=args.nevents
prefix=args.prefix
first_event=args.skip
arrays=ur.open(f'{infile}:events').arrays()
w0=args.w0_rw
MIP=args.MIP
w0_nrw=args.w0_nrw
w0_use_range=args.w0_use_range
if w0_use_range:
w0_nrw=np.array([[a] for a in np.linspace(3.0, 8.0, 21)])
w0=np.array([[a] for a in np.linspace(3.0, 8.0, 21)])
r_truths=[]
x_truths=[]
y_truths=[]
z_truths=[]
r_recos=[]
x_recos=[]
y_recos=[]
z_recos=[]
r_recos_rw=[]
x_recos_rw=[]
y_recos_rw=[]
z_recos_rw=[]
drs=[]
drs_rw=[]
dxs=[]
dys=[]
dxs_rw=[]
dys_rw=[]
Es=[]
mc_pzs=[]
w0s=[]
if nevents==-1:
nevents=len(arrays)
for event in range(first_event,first_event+nevents):
if len(arrays[f'{prefix}HitsReco.position.x'][event])== 0:
print("warning, empty event. skipping")
continue
#print(arrays['ZDCHitsReco.energy'][event])
#print(len(arrays['ZDCHitsReco.position.x'][event]))
x_reco, y_reco, z_reco, r_reco=get_xyzr_reco_no_reweighting(arrays, event, w0=w0_nrw, weight_by_granularity=True, prefix=prefix)
if useH3reweighting:
x_reco_rw, y_reco_rw, z_reco_rw, r_reco_rw=get_xyzr_reco_reweighted_H3(arrays, event, w0=w0, MIP=MIP, weight_by_granularity=True, prefix=prefix)
elif useH4reweighting:
ret=get_xyzr_reco_reweighted_H4(arrays, event, w0=w0, MIP=MIP, weight_by_granularity=True, prefix=prefix)
if ret is not None:
x_reco_rw, y_reco_rw, z_reco_rw, r_reco_rw=ret
else : continue
elif useS2reweighting:
x_reco_rw, y_reco_rw, z_reco_rw, r_reco_rw=get_xyzr_reco_reweighted_S2(arrays, event, w0=w0, MIP=MIP, weight_by_granularity=True, prefix=prefix)
x_truth, y_truth, z_truth, r_truth=get_xyzr_truth(arrays, event, w0=w0, weight_by_granularity=True, prefix=prefix)
#print(r_truth, r_reco)
drs.append(r_reco-r_truth)
dxs.append(x_reco-x_truth)
dys.append(y_reco-y_truth)
r_recos.append(r_reco)
x_recos.append(x_reco)
y_recos.append(y_reco)
z_recos.append(z_reco)
r_truths.append(r_truth)
x_truths.append(x_truth)
y_truths.append(y_truth)
z_truths.append(z_truth)
if useH3reweighting or useH4reweighting or useS2reweighting:
drs_rw.append(r_reco_rw-r_truth)
dxs_rw.append(x_reco_rw-x_truth)
dys_rw.append(y_reco_rw-y_truth)
r_recos_rw.append(r_reco_rw)
x_recos_rw.append(x_reco_rw)
y_recos_rw.append(y_reco_rw)
z_recos_rw.append(z_reco_rw)
Es.append(sum(arrays[f'{prefix}HitsReco.energy'][event]))
mc_pzs.append(arrays[f'MCParticles.momentum.z'][event,2])
if w0_use_range:
w0s.append(w0)
#except Exception as e:
# print(e)
# e#pass
if event%10==0:
print(f"{infile}: done with event {event}/{nevents}")
if not w0_use_range:
d=dict(E=Es, dr=drs, dy=dys, dx=dxs, mc_pz=mc_pzs, x_truth=x_truths, y_truth=y_truths, z_truth=z_truths,
x_reco=x_recos, y_reco=y_recos,z_reco=z_recos,r_reco=r_recos)
if useH3reweighting or useH4reweighting or useS2reweighting:
d["dr_rw"]=drs_rw
d["dx_rw"]=dxs_rw
d["dy_rw"]=dys_rw
d['r_reco_rw']=r_recos_rw
d['x_reco_rw']=x_recos_rw
d['y_reco_rw']=y_recos_rw
d['z_reco_rw']=z_recos_rw
else :
w0s = [a[0] for a in w0] # flatten array
d=dict(E=Es, dr=drs, dy=dys, dx=dxs, mc_pz=mc_pzs, x_truth=x_truths, y_truth=y_truths, z_truth=z_truths, r_truth=r_truths, x_reco=x_recos, y_reco=y_recos, z_reco=z_recos, r_reco=r_recos)
if useH3reweighting or useH4reweighting or useS2reweighting:
d["dr_rw"]=drs_rw
d["dx_rw"]=dxs_rw
d["dy_rw"]=dys_rw
d['r_reco_rw']=r_recos_rw
d['x_reco_rw']=x_recos_rw
d['y_reco_rw']=y_recos_rw
d['z_reco_rw']=z_recos_rw
for key in list(d.keys()):
#print(key)
#print(type(d[key]))
#print(d[key][0])
if key not in "E mc_pz w0s".split():
print("making new columns")
for i in range(len(w0s)):
d[f"{key}_w0_{w0s[i]}".replace(".", "pt")]=[d[key][j][i] for j in range(len(d[key]))]
print({a:len(d[a]) for a in d})
if ".csv" in outfile:
pd.DataFrame(d).to_csv(outfile)
if ".pkl" in outfile:
pd.DataFrame(d).to_pickle(outfile)
benchmarks/far_forward/options/zdc_sipmontile_reconstruction.py 0 → 100644
+ 91
− 0
View file @ 8bba950e
import os
import ROOT
from Gaudi.Configuration import *
from GaudiKernel.SystemOfUnits import mm, MeV, rad, ns
from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
from Configurables import PodioInput
from Configurables import Jug__Digi__CalorimeterHitDigi as CalHitDigi
from Configurables import Jug__Reco__CalorimeterHitReco as CalHitReco
from Configurables import Jug__Reco__CalorimeterHitsMerger as CalHitsMerge
from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
kwargs = dict()
# input arguments through environment variables
kwargs['input'] = os.environ.get('JUGGLER_SIM_FILE', '../epic/zdc_sipmontile_neutrons.edm4hep.root')
kwargs['output'] = os.environ.get('JUGGLER_REC_FILE', 'rec_zdc_sipmontile_neutrons.root')
kwargs['compact'] = os.environ.get('DETECTOR_COMPACT_PATH', '../epic/epic.xml')
kwargs['nev'] = int(os.environ.get('JUGGLER_N_EVENTS', 100))
kwargs['PbSci_sf'] = float(os.environ.get('ZDC_PbSCI_SAMP_FRAC', 1.0))
if kwargs['nev'] < 1:
f = ROOT.TFile(kwargs['input'])
kwargs['nev'] = f.events.GetEntries()
print(kwargs)
geo_service = GeoSvc('GeoSvc', detectors=kwargs['compact'].split(','), OutputLevel=INFO)
podioevent = EICDataSvc('EventDataSvc', inputs=kwargs['input'].split(','), OutputLevel=DEBUG)
sim_colls = [
'MCParticles',
'HcalFarForwardZDCHits',
'HcalFarForwardZDCHitsContributions',
]
podin = PodioInput('PodioReader', collections=sim_colls, OutputLevel=DEBUG)
podout = PodioOutput('out', filename=kwargs['output'])
# Digitization
ci_zdc_daq = dict(\
dynamicRangeADC=200.*MeV,
capacityADC=32768,
pedestalMean=400,
pedestalSigma=10)
SiPMonTile_digi = CalHitDigi('SiPMonTile_digi',
inputHitCollection='HcalFarForwardZDCHits',
outputHitCollection='HcalFarForwardZDCRawHits',
**ci_zdc_daq)
# Reconstruction
SiPMonTile_reco = CalHitReco('SiPMonTile_reco',
inputHitCollection=SiPMonTile_digi.outputHitCollection,
outputHitCollection='HcalFarForwardZDCHitsReco',
readoutClass='HcalFarForwardZDCHits',**ci_zdc_daq,
thresholdFactor=0.0,
samplingFraction="1.0")
# localDetFields=['system'])
# Clustering
'''SiPMonTile_cl = IslandCluster('SiPMonTile_cl',
inputHitCollection=SiPMonTile_reco.outputHitCollection,
outputProtoClusterCollection='HcalFarForwardZDCProtoClusters',
minClusterCenterEdep=3.*MeV,
minClusterHitEdep=0.1*MeV,
localDistXY=[50*mm, 50*mm],
splitCluster=True)
SiPMonTile_clreco = RecoCoG('SiPMonTile_clreco',
inputProtoClusterCollection = SiPMonTile_cl.outputProtoClusterCollection,
outputClusterCollection='HcalFarForwardZDCClusters',
logWeightBase=3.6)
'''
podout.outputCommands = ['keep *']
ApplicationMgr(
TopAlg=[podin,
SiPMonTile_digi,
SiPMonTile_reco,
# SiPMonTile_cl, SiPMonTile_clreco,
podout],
EvtSel='NONE',
EvtMax=kwargs['nev'],
ExtSvc=[podioevent],
OutputLevel=DEBUG
)
benchmarks/far_forward/run_zdc_sipmontile.sh 0 → 100644
+ 154
− 0
View file @ 8bba950e
#!/bin/bash
function print_the_help {
echo "USAGE: ${0} [--rec-only] [--ana-only] [--particle=neutron] "
echo " OPTIONS: "
echo " --rec-only only run gaudi reconstruction part"
echo " --ana-only only run analysis macros part"
exit
}
PARTICLE="neutron"
REC_ONLY=
ANALYSIS_ONLY=
POSITIONAL=()
while [[ $# -gt 0 ]]
do
key="$1"
case $key in
-h|--help)
shift # past argument
print_the_help
;;
--rec-only)
REC_ONLY=1
shift # past value
;;
--ana-only)
ANALYSIS_ONLY=1
shift # past value
;;
--particle)
PARTICLE="$2"
shift # past argument
shift # past value
;;
*) # unknown option
#POSITIONAL+=("$1") # save it in an array for later
echo "unknown option $1"
print_the_help
shift # past argument
;;
esac
done
set -- "${POSITIONAL[@]}" # restore positional parameters
print_env.sh
## To run the reconstruction, we need the following global variables:
## - DETECTOR: the detector package we want to use for this benchmark
## - DETECTOR_VERSION: the detector package we want to use for this benchmark
## - DETECTOR_PATH: full path to the detector definitions
##
## You can ready options/env.sh for more in-depth explanations of the variables
## and how they can be controlled.
export DETECTOR_PATH=${DETECTOR_PATH}
if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then
export JUGGLER_N_EVENTS=100
fi
if [[ ! -n "${E_start}" ]] ; then
export E_start=10.0
fi
if [[ ! -n "${E_end}" ]] ; then
export E_end=10.0
fi
#if [[ ! -n "${ZDC_SCI_SAMP_FRAC}" ]] ; then
# export ZDC_PbSCI_SAMP_FRAC=1.0
#fi
export DETECTOR_CONFIG="epic_zdc_sipm_on_tile_only"
echo "DETECTOR_CONFIG is set to "$DETECTOR_CONFIG
export DETECTOR_COMPACT_PATH=${DETECTOR_PATH}/${DETECTOR_CONFIG}.xml
export JUGGLER_FILE_NAME_TAG="zdc_sipmontile_${PARTICLE}"
export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.edm4hep.root"
export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
export RESULTS_PATH="results/far_forward/zdc_sipmontile/${PARTICLE}"
export PHYSICS_LIST="FTFP_BERT"
mkdir -p ${RESULTS_PATH}
echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
echo "DETECTOR = ${DETECTOR}"
if [[ -z "${REC_ONLY}" && -z "${ANALYSIS_ONLY}" ]] ;
then
echo "Generating input events"
root -b -q "benchmarks/far_forward/analysis/gen_zdc_sipmontile_particles.cxx(${JUGGLER_N_EVENTS}, \"${PARTICLE}\", ${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running script: generating input events"
exit 1
fi
echo "Running Geant4 simulation"
npsim \
--numberOfEvents ${JUGGLER_N_EVENTS} \
--physicsList ${PHYSICS_LIST} \
--compactFile ${DETECTOR_PATH}/${DETECTOR_CONFIG}.xml \
--inputFiles ${JUGGLER_FILE_NAME_TAG}.hepmc \
--outputFile ${JUGGLER_SIM_FILE}
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running npdet: Geant4 simulation"
exit 1
fi
fi
rootls -t ${JUGGLER_SIM_FILE}
if [[ -z "${ANALYSIS_ONLY}" ]] ;
then
gaudirun.py benchmarks/far_forward/options/zdc_sipmontile_reconstruction.py
status=$?
if [[ "$status" -ne "0" && "$status" -ne "4" ]] ; then
echo "ERROR running juggler, got $status"
exit 1
fi
fi
rootls -t ${JUGGLER_REC_FILE}
HEXPLIT_OUT_FILE=${JUGGLER_REC_FILE}_hexplit.csv
python benchmarks/far_forward/hexplit/hexplit.py ${JUGGLER_REC_FILE} ${HEXPLIT_OUT_FILE} --prefix HcalFarForwardZDC --w0_rw 5 --w0_nrw 5 --H4
ls $HEXPLIT_OUT_FILE
echo "Running analysis root scripts"
root -b -q "benchmarks/far_forward/analysis/analysis_zdc_sipmontile_particles.cxx+(\"${JUGGLER_REC_FILE}\", \"${RESULTS_PATH}\")"
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running root analysis script"
exit 1
fi
echo "making plots of hexplit residuals"
python benchmarks/far_forward/analysis/make_plots_zdc_sipmontile_hexplit.py $HEXPLIT_OUT_FILE ${RESULTS_PATH}
root_filesize=$(stat --format=%s "${JUGGLER_REC_FILE}")
if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then
# file must be less than 10 MB to upload
if [[ "${root_filesize}" -lt "10000000" ]] ; then
cp ${JUGGLER_SIM_FILE} results/far_forward/zdc_sipmontile/.
cp ${JUGGLER_REC_FILE} results/far_forward/zdc_sipmontile/.
fi
fi
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment