From 5b432632b08b5160efa0e5e99c715d6a692e1c81 Mon Sep 17 00:00:00 2001
From: Barak Schmookler <bschmookler1@gmail.com>
Date: Sat, 9 Nov 2024 20:51:24 -0600
Subject: [PATCH] Added benchmark for single neutron in the ZDC (#99)
---
.gitlab-ci.yml | 2 +
Snakefile | 1 +
benchmarks/zdc_neutron/README.md | 10 +
benchmarks/zdc_neutron/Snakefile | 65 ++++
.../zdc_neutron/analysis/fwd_neutrons_geant.C | 331 ++++++++++++++++++
.../zdc_neutron/analysis/fwd_neutrons_recon.C | 306 ++++++++++++++++
benchmarks/zdc_neutron/config.yml | 38 ++
.../zdc_neutron/gen_forward_neutrons.cxx | 102 ++++++
8 files changed, 855 insertions(+)
create mode 100644 benchmarks/zdc_neutron/README.md
create mode 100644 benchmarks/zdc_neutron/Snakefile
create mode 100644 benchmarks/zdc_neutron/analysis/fwd_neutrons_geant.C
create mode 100644 benchmarks/zdc_neutron/analysis/fwd_neutrons_recon.C
create mode 100644 benchmarks/zdc_neutron/config.yml
create mode 100644 benchmarks/zdc_neutron/gen_forward_neutrons.cxx
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 4d77de9b..e313531b 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -136,6 +136,7 @@ include:
- local: 'benchmarks/lfhcal/config.yml'
- local: 'benchmarks/zdc/config.yml'
- local: 'benchmarks/zdc_lyso/config.yml'
+ - local: 'benchmarks/zdc_neutron/config.yml'
- local: 'benchmarks/zdc_photon/config.yml'
- local: 'benchmarks/zdc_pi0/config.yml'
- local: 'benchmarks/material_scan/config.yml'
@@ -170,6 +171,7 @@ deploy_results:
- "collect_results:tracking_performances_dis"
- "collect_results:zdc"
- "collect_results:zdc_lyso"
+ - "collect_results:zdc_neutron"
- "collect_results:insert_muon"
- "collect_results:insert_tau"
- "collect_results:zdc_photon"
diff --git a/Snakefile b/Snakefile
index b913cc80..23ae5169 100644
--- a/Snakefile
+++ b/Snakefile
@@ -9,6 +9,7 @@ include: "benchmarks/tracking_performances/Snakefile"
include: "benchmarks/tracking_performances_dis/Snakefile"
include: "benchmarks/lfhcal/Snakefile"
include: "benchmarks/zdc_lyso/Snakefile"
+include: "benchmarks/zdc_neutron/Snakefile"
include: "benchmarks/insert_muon/Snakefile"
include: "benchmarks/zdc_lambda/Snakefile"
include: "benchmarks/zdc_photon/Snakefile"
diff --git a/benchmarks/zdc_neutron/README.md b/benchmarks/zdc_neutron/README.md
new file mode 100644
index 00000000..623d1465
--- /dev/null
+++ b/benchmarks/zdc_neutron/README.md
@@ -0,0 +1,10 @@
+Detector Benchmark for single neutron in the ZDC
+=================================================
+
+## Overview
+This benchmark generates single-neutron events. The neutrons are generated with 100 GeV total momentum; polar angles of 0-6 mRad with respect to the proton/ion beam direction, uniform over cosine of the polar angle; and uniform azimuthal angles with respect to the proton/ion beam direction. The benchmark creates acceptance and reconstruction performance plots.
+
+## Contacts
+[Barak Schmookler](baraks@ucr.edu)
+
+
diff --git a/benchmarks/zdc_neutron/Snakefile b/benchmarks/zdc_neutron/Snakefile
new file mode 100644
index 00000000..4b4c0d36
--- /dev/null
+++ b/benchmarks/zdc_neutron/Snakefile
@@ -0,0 +1,65 @@
+# Generate the single neutrons and put them into a HepMC file
+rule zdc_neutron_hepmc:
+ input:
+ script = "benchmarks/zdc_neutron/gen_forward_neutrons.cxx",
+ output:
+ hepmcfile="sim_output/zdc_neutron/{DETECTOR_CONFIG}/fwd_neutrons.hepmc",
+ params:
+ num_events=1000,
+ shell:
+ """
+root -l -b -q '{input.script}({params.num_events}, 0, "{output.hepmcfile}")'
+"""
+
+# Run the generated events through the Geant simulation
+rule zdc_neutron_sim:
+ input:
+ hepmcfile="sim_output/zdc_neutron/{DETECTOR_CONFIG}/fwd_neutrons.hepmc",
+ warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root",
+ output:
+ "sim_output/zdc_neutron/{DETECTOR_CONFIG}/fwd_neutrons.edm4hep.root",
+ params:
+ num_events=100,
+ shell:
+ """
+set -m # monitor mode to prevent lingering processes
+exec npsim \
+ --runType batch \
+ -v WARNING \
+ --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \
+ --numberOfEvents {params.num_events} \
+ --inputFiles {input.hepmcfile} \
+ --outputFile {output}
+"""
+
+
+# Process the file produced in the previous step through EICRecon
+rule zdc_neutron_reco:
+ input:
+ "sim_output/zdc_neutron/{DETECTOR_CONFIG}/fwd_neutrons.edm4hep.root",
+ output:
+ "sim_output/zdc_neutron/{DETECTOR_CONFIG}/fwd_neutrons.edm4eic.root",
+ shell:
+ """
+set -m # monitor mode to prevent lingering processes
+exec env DETECTOR_CONFIG={wildcards.DETECTOR_CONFIG} \
+ eicrecon {input} -Ppodio:output_file={output} \
+ -Ppodio:output_collections=MCParticles,EcalFarForwardZDCRawHits,EcalFarForwardZDCRecHits,EcalFarForwardZDCClusters,HcalFarForwardZDCRawHits,HcalFarForwardZDCRecHits,HcalFarForwardZDCClusters,ReconstructedFarForwardZDCNeutrons
+"""
+
+
+# Run the analysis scripts
+rule zdc_neutron_analyses:
+ input:
+ geant_script = "benchmarks/zdc_neutron/analysis/fwd_neutrons_geant.C",
+ data_geant = "sim_output/zdc_neutron/{DETECTOR_CONFIG}/fwd_neutrons.edm4hep.root",
+ recon_script = "benchmarks/zdc_neutron/analysis/fwd_neutrons_recon.C",
+ data_recon = "sim_output/zdc_neutron/{DETECTOR_CONFIG}/fwd_neutrons.edm4eic.root",
+ output:
+ geant_analysis_out = "results/zdc_neutron/{DETECTOR_CONFIG}/fwd_neutrons_geant.pdf",
+ recon_analysis_out = "results/zdc_neutron/{DETECTOR_CONFIG}/fwd_neutrons_recon.pdf",
+ shell:
+ """
+root -l -b -q '{input.geant_script}("{input.data_geant}","{output.geant_analysis_out}")'
+root -l -b -q '{input.recon_script}("{input.data_recon}","{output.recon_analysis_out}")'
+"""
diff --git a/benchmarks/zdc_neutron/analysis/fwd_neutrons_geant.C b/benchmarks/zdc_neutron/analysis/fwd_neutrons_geant.C
new file mode 100644
index 00000000..0e4c7e92
--- /dev/null
+++ b/benchmarks/zdc_neutron/analysis/fwd_neutrons_geant.C
@@ -0,0 +1,331 @@
+//------------------
+int get_layer_number(TVector3 pos){
+
+ // Get Layer position wrt proton beam direction
+ pos.RotateY(0.025);
+
+ auto local_z = pos.Z() - 35.8*1000.; //in mm
+
+ // Convert to layer number
+ // local_z = 22.25 + (layer_number - 1)*24.9
+ int layer_number = (int)std::round( (local_z - 22.25)/24.9 ) + 1;
+
+ return layer_number;
+}
+
+//------------------
+void fwd_neutrons_geant(std::string inputfile, std::string outputfile){
+
+ //Define Style
+ gStyle->SetOptStat(0);
+ gStyle->SetPadBorderMode(0);
+ gStyle->SetFrameBorderMode(0);
+ gStyle->SetFrameLineWidth(2);
+ gStyle->SetLabelSize(0.035,"X");
+ gStyle->SetLabelSize(0.035,"Y");
+ //gStyle->SetLabelOffset(0.01,"X");
+ //gStyle->SetLabelOffset(0.01,"Y");
+ gStyle->SetTitleXSize(0.04);
+ gStyle->SetTitleXOffset(0.9);
+ gStyle->SetTitleYSize(0.04);
+ gStyle->SetTitleYOffset(0.9);
+
+ //Define histograms
+ TH2 *h1_neut = new TH2D("h1_neut","Neutron true energy vs. polar angle",100,0.6,2.2,100,0,200);
+ h1_neut->GetXaxis()->SetTitle("#theta [deg]"); h1_neut->GetXaxis()->CenterTitle();
+ h1_neut->GetYaxis()->SetTitle("E [GeV]"); h1_neut->GetYaxis()->CenterTitle();
+
+ TH2 *h2_neut = new TH2D("h2_neut","Neutron true azimuthal angle vs. polar angle around p axis",100,-0.1,12,100,-200,200);
+ h2_neut->GetXaxis()->SetTitle("#theta^{*} [mRad]"); h2_neut->GetXaxis()->CenterTitle();
+ h2_neut->GetYaxis()->SetTitle("#phi^{*} [deg]"); h2_neut->GetYaxis()->CenterTitle();
+
+ TH2 *h2a_neut = new TH2D("h2a_neut","Neutron true azimuthal angle vs. cosine of polar angle around p axis",100,0.99996,1,100,-200,200);
+ h2a_neut->GetXaxis()->SetTitle("cos(#theta^{*})"); h2a_neut->GetXaxis()->CenterTitle();
+ h2a_neut->GetYaxis()->SetTitle("#phi^{*} [deg]"); h2a_neut->GetYaxis()->CenterTitle();
+
+ //Require HCal hit energy sum > 1 GeV
+ TH2 *h3_neut = new TH2D("h3_neut","Neutron true azimuthal angle vs. polar angle around p axis",100,-0.1,12,100,-200,200);
+ h3_neut->GetXaxis()->SetTitle("#theta^{*} [mRad]"); h3_neut->GetXaxis()->CenterTitle();
+ h3_neut->GetYaxis()->SetTitle("#phi^{*} [deg]"); h3_neut->GetYaxis()->CenterTitle();
+
+ TH2 *h4_neut = new TH2D("h4_neut","Neutron local hit position at ZDC HCal front face",100,-300,300,100,-300,300);
+ h4_neut->GetXaxis()->SetTitle("x [mm]"); h4_neut->GetXaxis()->CenterTitle();
+ h4_neut->GetYaxis()->SetTitle("y [mm]"); h4_neut->GetYaxis()->CenterTitle();
+
+ TH1 *h1_ecal = new TH1D("h1_ecal","Total true hit energy sum",100,-0.1,2);
+ h1_ecal->GetXaxis()->SetTitle("E [GeV]"); h1_ecal->GetXaxis()->CenterTitle();
+ h1_ecal->SetLineColor(kRed);h1_ecal->SetLineWidth(2);
+
+ TH1 *h1_hcal = new TH1D("h1_hcal","Total true hit energy sum",100,-0.1,2);
+ h1_hcal->GetXaxis()->SetTitle("E [GeV]"); h1_hcal->GetXaxis()->CenterTitle();
+ h1_hcal->SetLineColor(kBlue);h1_hcal->SetLineWidth(2);
+
+ //Cut on theta*
+ TH1 *h2_ecal = new TH1D("h2_ecal","Total true hit energy sum",100,-0.1,2);
+ h2_ecal->GetXaxis()->SetTitle("E [GeV]"); h2_ecal->GetXaxis()->CenterTitle();
+ h2_ecal->SetLineColor(kRed);h2_ecal->SetLineWidth(2);
+
+ //Cut on theta*
+ TH1 *h2_hcal = new TH1D("h2_hcal","Total true hit energy sum",100,-0.1,2);
+ h2_hcal->GetXaxis()->SetTitle("E [GeV]"); h2_hcal->GetXaxis()->CenterTitle();
+ h2_hcal->SetLineColor(kBlue);h2_hcal->SetLineWidth(2);
+
+ //Cut on theta*
+ TH1 *h2a_ecal = new TH1D("h2a_ecal","Total true hit energy sum",200,-0.1,30);
+ h2a_ecal->GetXaxis()->SetTitle("E [GeV]"); h2a_ecal->GetXaxis()->CenterTitle();
+ h2a_ecal->SetLineColor(kRed);h2a_ecal->SetLineWidth(2);
+
+ //Cut on theta*
+ TH1 *h2a_hcal = new TH1D("h2a_hcal","Total true hit energy sum",200,-0.1,30);
+ h2a_hcal->GetXaxis()->SetTitle("E [GeV]"); h2a_hcal->GetXaxis()->CenterTitle();
+ h2a_hcal->SetLineColor(kBlue);h2a_hcal->SetLineWidth(2);
+
+ //Cut on theta*
+ TH2 *h2_cal_both = new TH2D("h2_cal_both","Total true hit energy sum: ECal vs. HCal",100,-0.1,2.5,100,-0.1,30);
+ h2_cal_both->GetXaxis()->SetTitle("E_{HCal.} [GeV]"); h2_cal_both->GetXaxis()->CenterTitle();
+ h2_cal_both->GetYaxis()->SetTitle("E_{ECal.} [GeV]"); h2_cal_both->GetYaxis()->CenterTitle();
+
+ TH1 *h3a_hcal = new TH1D("h3a_hcal","Cells w/ non-zero energy deposit ",100,0,64);
+ h3a_hcal->GetXaxis()->SetTitle("Layer Number in ZDC HCal"); h3a_hcal->GetXaxis()->CenterTitle();
+ h3a_hcal->GetYaxis()->SetTitle("Number of cells hit per event");h3a_hcal->GetYaxis()->CenterTitle();
+ h3a_hcal->SetLineColor(kBlue);h3a_hcal->SetLineWidth(3);
+
+ //Cut on theta* < 3.5 mRad
+ TH1 *h3b_hcal = new TH1D("h3b_hcal","Cells w/ non-zero energy deposit ",100,0,64);
+ h3b_hcal->GetXaxis()->SetTitle("Layer Number in ZDC HCal"); h3b_hcal->GetXaxis()->CenterTitle();
+ h3b_hcal->GetYaxis()->SetTitle("Number of cells hit per event");h3b_hcal->GetYaxis()->CenterTitle();
+ h3b_hcal->SetLineColor(kBlue);h3b_hcal->SetLineWidth(3);
+
+ //Cut on 3.5 mRad < theta* < 6 mRad
+ TH1 *h3c_hcal = new TH1D("h3c_hcal","Cells w/ non-zero energy deposit ",100,0,64);
+ h3c_hcal->GetXaxis()->SetTitle("Layer Number in ZDC HCal"); h3c_hcal->GetXaxis()->CenterTitle();
+ h3c_hcal->GetYaxis()->SetTitle("Number of cells hit per event");h3c_hcal->GetYaxis()->CenterTitle();
+ h3c_hcal->SetLineColor(kBlue);h3c_hcal->SetLineWidth(3);
+
+ //Read ROOT file
+ TFile* file = new TFile(inputfile.c_str());
+ TTree *tree = (TTree*) file->Get("events");
+
+ //Set cut
+ float edep_cut = 1.; //In GeV
+
+ cout<<"Total number of events to analyze is "<<tree->GetEntries()<<endl;
+
+ //Create Array Reader
+ TTreeReader tr(tree);
+
+ TTreeReaderArray<int> gen_status(tr, "MCParticles.generatorStatus");
+ TTreeReaderArray<int> gen_pid(tr, "MCParticles.PDG");
+ TTreeReaderArray<float> gen_px(tr, "MCParticles.momentum.x");
+ TTreeReaderArray<float> gen_py(tr, "MCParticles.momentum.y");
+ TTreeReaderArray<float> gen_pz(tr, "MCParticles.momentum.z");
+ TTreeReaderArray<double> gen_mass(tr, "MCParticles.mass");
+ TTreeReaderArray<float> gen_charge(tr, "MCParticles.charge");
+ TTreeReaderArray<double> gen_vx(tr, "MCParticles.vertex.x");
+ TTreeReaderArray<double> gen_vy(tr, "MCParticles.vertex.y");
+ TTreeReaderArray<double> gen_vz(tr, "MCParticles.vertex.z");
+
+ //ZDC LYSO ECal
+ TTreeReaderArray<float> ecal_hit_e(tr,"EcalFarForwardZDCHits.energy");
+
+ //ZDC SiPM-on-tile HCal
+ TTreeReaderArray<float> hcal_hit_e(tr, "HcalFarForwardZDCHits.energy");
+ TTreeReaderArray<float> hcal_hit_x(tr, "HcalFarForwardZDCHits.position.x");
+ TTreeReaderArray<float> hcal_hit_y(tr, "HcalFarForwardZDCHits.position.y");
+ TTreeReaderArray<float> hcal_hit_z(tr, "HcalFarForwardZDCHits.position.z");
+
+ //Other variables
+ int counter(0),counter_3p5(0),counter_3p5to6(0);
+
+ TLorentzVector neut_true; //True neutron in lab coordinates
+ TLorentzVector neut_true_rot; //True neutron wrt proton beam direction
+
+ float ecal_e_tot(0);
+ float hcal_e_tot(0);
+
+ //Loop over events
+ while (tr.Next()) {
+
+ if(counter%1000==0) cout<<"Analyzing event "<<counter<<endl;
+ counter++;
+
+ //Reset variables
+ ecal_e_tot = 0;
+ hcal_e_tot = 0;
+
+ //Loop over generated particles, select primary neutron
+ for(int igen=0;igen<gen_status.GetSize();igen++){
+ if(gen_status[igen]==1 && gen_pid[igen]==2112){
+
+ neut_true.SetXYZM(gen_px[igen],gen_py[igen],gen_pz[igen],gen_mass[igen]);
+ h1_neut->Fill(neut_true.Theta()*TMath::RadToDeg(),neut_true.E());
+
+ //Wrt proton beam direction
+ neut_true_rot = neut_true;
+ neut_true_rot.RotateY(0.025);
+ h2_neut->Fill(neut_true_rot.Theta()*1000.,neut_true_rot.Phi()*TMath::RadToDeg()); //Theta* in mRad
+ h2a_neut->Fill(std::cos(neut_true_rot.Theta()),neut_true_rot.Phi()*TMath::RadToDeg());
+
+ }
+ } //End loop over generated particles
+
+ //Loop over ECal hits
+ for(int ihit=0;ihit<ecal_hit_e.GetSize();ihit++){
+ ecal_e_tot += ecal_hit_e[ihit];
+ }
+
+ //Loop over HCal hits
+ for(int ihit=0;ihit<hcal_hit_e.GetSize();ihit++){
+ hcal_e_tot += hcal_hit_e[ihit];
+
+ //Cell hit position in global coordinates
+ double hit_x = (double) hcal_hit_x[ihit];
+ double hit_y = (double) hcal_hit_y[ihit];
+ double hit_z = (double) hcal_hit_z[ihit];
+
+ TVector3 hit_pos_det(hit_x,hit_y,hit_z);
+ auto layer_number = get_layer_number(hit_pos_det);
+
+ h3a_hcal->Fill(layer_number);
+ if( neut_true_rot.Theta()*1000. < 3.5 ) h3b_hcal->Fill(layer_number);
+ if( neut_true_rot.Theta()*1000. > 3.5 && neut_true_rot.Theta()*1000 < 6.) h3c_hcal->Fill(layer_number);
+
+ } //End loop over HCal hits
+
+ //Fill histograms
+ h1_ecal->Fill(ecal_e_tot);
+ h1_hcal->Fill(hcal_e_tot);
+
+ if(hcal_e_tot > edep_cut){
+ h3_neut->Fill(neut_true_rot.Theta()*1000.,neut_true_rot.Phi()*TMath::RadToDeg()); //Theta* in mRad
+
+ //Projection to front face of HCal
+ auto proj_x = 35.8 * 1000. * tan(neut_true_rot.Theta()) * cos(neut_true_rot.Phi()) ;
+ auto proj_y = 35.8 * 1000. * tan(neut_true_rot.Theta()) * sin(neut_true_rot.Phi()) ;
+
+ h4_neut->Fill(proj_x,proj_y);
+ }
+
+ if( neut_true_rot.Theta()*1000. < 3.5 ){
+ h2_ecal->Fill(ecal_e_tot);
+ h2_hcal->Fill(hcal_e_tot);
+ h2a_ecal->Fill(ecal_e_tot);
+ h2a_hcal->Fill(hcal_e_tot);
+
+ h2_cal_both->Fill(hcal_e_tot,ecal_e_tot);
+
+ counter_3p5++;
+ }
+
+ if( neut_true_rot.Theta()*1000. > 3.5 && neut_true_rot.Theta()*1000 < 6.){
+ counter_3p5to6++;
+ }
+
+ } //End loop over events
+
+ //Scale histograms
+ h3a_hcal->Scale(1./counter);
+ h3b_hcal->Scale(1./counter_3p5);
+ h3c_hcal->Scale(1./counter_3p5to6);
+
+ //Make plots
+ TCanvas *c1 = new TCanvas("c1");
+ h1_neut->Draw("colz");
+
+ TCanvas *c2 = new TCanvas("c2");
+ h2_neut->Draw("colz");
+
+ TCanvas *c2a = new TCanvas("c2a");
+ h2a_neut->Draw("colz");
+
+ TCanvas *c3 = new TCanvas("c3");
+ h1_ecal->Draw("");
+ h1_hcal->Draw("same");
+
+ TLegend *leg3 = new TLegend(0.6,0.6,0.85,0.8);
+ leg3->SetBorderSize(0);leg3->SetFillStyle(0);
+ leg3->AddEntry(h1_ecal,"Sum of ZDC Ecal hit energies","l");
+ leg3->AddEntry(h1_hcal,"Sum of ZDC Hcal hit energies","l");
+ leg3->Draw();
+
+ TCanvas *c4 = new TCanvas("c4");
+ h3_neut->Draw("colz");
+
+ TLatex *tex4 = new TLatex(7,150,Form("Hcal hit energy sum > %.1f GeV",edep_cut));
+ tex4->SetTextSize(0.035);
+ tex4->Draw();
+
+ TCanvas *c4a = new TCanvas("c4a");
+ h4_neut->Draw("colz");
+
+ TLatex *tex4a = new TLatex(50,250,Form("Hcal hit energy sum > %.1f GeV",edep_cut));
+ tex4a->SetTextSize(0.035);
+ tex4a->Draw();
+
+ TCanvas *c5 = new TCanvas("c5");
+ h2_ecal->Draw("");
+ h2_hcal->Draw("same");
+
+ TLegend *leg5 = new TLegend(0.5,0.6,0.9,0.8);
+ leg5->SetBorderSize(0);leg5->SetFillStyle(0);
+ leg5->SetHeader("Require neutron #theta^{*} (wrt proton beam) < 3.5 mRad");
+ leg5->AddEntry(h1_ecal,"Sum of ZDC Ecal hit energies","l");
+ leg5->AddEntry(h1_hcal,"Sum of ZDC Hcal hit energies","l");
+ leg5->Draw();
+
+ TCanvas *c5a = new TCanvas("c5a");
+ c5a->SetLogy();
+ h2a_ecal->Draw("");
+ h2a_hcal->Draw("same");
+
+ leg5->Draw();
+
+ TCanvas *c5b = new TCanvas("c5b");
+ c5b->SetLogz();
+ h2_cal_both->Draw("colz");
+
+ TLatex *tex5b = new TLatex(0,28,"Require neutron #theta^{*} (wrt proton beam) < 3.5 mRad");
+ tex5b->SetTextSize(0.03);
+ tex5b->Draw();
+
+ TCanvas *c6a = new TCanvas("c6a");
+ h3a_hcal->Draw();
+
+ TLegend *leg6a = new TLegend(0.5,0.6,0.9,0.8);
+ leg6a->SetBorderSize(0);leg6a->SetFillStyle(0);
+ leg6a->SetHeader("Average over all events");
+ leg6a->Draw();
+
+ TCanvas *c6b = new TCanvas("c6b");
+ h3b_hcal->Draw();
+
+ TLegend *leg6b = new TLegend(0.5,0.6,0.9,0.8);
+ leg6b->SetBorderSize(0);leg6b->SetFillStyle(0);
+ leg6b->SetHeader("Average over events w/ neutron #theta^{*} < 3.5 mRad");
+ leg6b->Draw();
+
+ TCanvas *c6c = new TCanvas("c6c");
+ h3c_hcal->Draw();
+
+ TLegend *leg6c = new TLegend(0.5,0.6,0.9,0.8);
+ leg6c->SetBorderSize(0);leg6c->SetFillStyle(0);
+ leg6c->SetHeader("Average over events w/ neutron 3.5 < #theta^{*} < 6 mRad");
+ leg6c->Draw();
+
+ //Print plots to file
+ c1->Print(Form("%s[",outputfile.c_str()));
+ c1->Print(Form("%s",outputfile.c_str()));
+ c2->Print(Form("%s",outputfile.c_str()));
+ c2a->Print(Form("%s",outputfile.c_str()));
+ c3->Print(Form("%s",outputfile.c_str()));
+ c4->Print(Form("%s",outputfile.c_str()));
+ c4a->Print(Form("%s",outputfile.c_str()));
+ c5->Print(Form("%s",outputfile.c_str()));
+ c5a->Print(Form("%s",outputfile.c_str()));
+ c5b->Print(Form("%s",outputfile.c_str()));
+ c6a->Print(Form("%s",outputfile.c_str()));
+ c6b->Print(Form("%s",outputfile.c_str()));
+ c6c->Print(Form("%s",outputfile.c_str()));
+ c6c->Print(Form("%s]",outputfile.c_str()));
+
+}
diff --git a/benchmarks/zdc_neutron/analysis/fwd_neutrons_recon.C b/benchmarks/zdc_neutron/analysis/fwd_neutrons_recon.C
new file mode 100644
index 00000000..8a8a7b90
--- /dev/null
+++ b/benchmarks/zdc_neutron/analysis/fwd_neutrons_recon.C
@@ -0,0 +1,306 @@
+//------------------
+void fwd_neutrons_recon(std::string inputfile, std::string outputfile){
+
+ //Define Style
+ gStyle->SetOptStat(0);
+ gStyle->SetPadBorderMode(0);
+ gStyle->SetFrameBorderMode(0);
+ gStyle->SetFrameLineWidth(2);
+ gStyle->SetLabelSize(0.035,"X");
+ gStyle->SetLabelSize(0.035,"Y");
+ //gStyle->SetLabelOffset(0.01,"X");
+ //gStyle->SetLabelOffset(0.01,"Y");
+ gStyle->SetTitleXSize(0.04);
+ gStyle->SetTitleXOffset(0.9);
+ gStyle->SetTitleYSize(0.04);
+ gStyle->SetTitleYOffset(0.9);
+
+ //Define histograms
+ TH2 *h1_neut = new TH2D("h1_neut","Neutron true energy vs. polar angle",100,0.6,2.2,100,0,200);
+ h1_neut->GetXaxis()->SetTitle("#theta [deg]"); h1_neut->GetXaxis()->CenterTitle();
+ h1_neut->GetYaxis()->SetTitle("E [GeV]"); h1_neut->GetYaxis()->CenterTitle();
+
+ TH2 *h2_neut = new TH2D("h2_neut","Neutron true azimuthal angle vs. polar angle around p axis",100,-0.1,12,100,-200,200);
+ h2_neut->GetXaxis()->SetTitle("#theta^{*} [mRad]"); h2_neut->GetXaxis()->CenterTitle();
+ h2_neut->GetYaxis()->SetTitle("#phi^{*} [deg]"); h2_neut->GetYaxis()->CenterTitle();
+
+ TH2 *h2a_neut = new TH2D("h2a_neut","Neutron true azimuthal angle vs. cosine of polar angle around p axis",100,0.99996,1,100,-200,200);
+ h2a_neut->GetXaxis()->SetTitle("cos(#theta^{*})"); h2a_neut->GetXaxis()->CenterTitle();
+ h2a_neut->GetYaxis()->SetTitle("#phi^{*} [deg]"); h2a_neut->GetYaxis()->CenterTitle();
+
+ TH1 *h1_ecal_adc = new TH1D("h1_ecal_adc","ECal ADC amplitude spectrum",1000,-0.1,35000);
+ h1_ecal_adc->GetXaxis()->SetTitle("ADC Channel"); h1_ecal_adc->GetXaxis()->CenterTitle();
+ h1_ecal_adc->SetLineColor(kRed);h1_ecal_adc->SetLineWidth(2);
+
+ TH1 *h1_hcal_adc = new TH1D("h1_hcal_adc","HCal ADC amplitude spectrum",1000,-0.1,35000);
+ h1_hcal_adc->GetXaxis()->SetTitle("ADC Channel"); h1_hcal_adc->GetXaxis()->CenterTitle();
+ h1_hcal_adc ->SetLineColor(kBlue);h1_hcal_adc->SetLineWidth(2);
+
+ TH1 *h1_ecal = new TH1D("h1_ecal","Total reconstructed hit energy sum",100,-0.1,2);
+ h1_ecal->GetXaxis()->SetTitle("E [GeV]"); h1_ecal->GetXaxis()->CenterTitle();
+ h1_ecal->SetLineColor(kRed);h1_ecal->SetLineWidth(2);
+
+ TH1 *h1_hcal = new TH1D("h1_hcal","Total reconstructed hit energy sum",100,-0.1,2);
+ h1_hcal->GetXaxis()->SetTitle("E [GeV]"); h1_hcal->GetXaxis()->CenterTitle();
+ h1_hcal->SetLineColor(kBlue);h1_hcal->SetLineWidth(2);
+
+ //Cut on theta*
+ TH1 *h2_ecal = new TH1D("h2_ecal","Total reconstructed hit energy sum",100,-0.1,2);
+ h2_ecal->GetXaxis()->SetTitle("E [GeV]"); h2_ecal->GetXaxis()->CenterTitle();
+ h2_ecal->SetLineColor(kRed);h2_ecal->SetLineWidth(2);
+
+ //Cut on theta*
+ TH1 *h2_hcal = new TH1D("h2_hcal","Total reconstructed hit energy sum",100,-0.1,2);
+ h2_hcal->GetXaxis()->SetTitle("E [GeV]"); h2_hcal->GetXaxis()->CenterTitle();
+ h2_hcal->SetLineColor(kBlue);h2_hcal->SetLineWidth(2);
+
+ //Cut on theta*
+ TH1 *h3_ecal = new TH1D("h3_ecal","Number of reconstructed clusters in ECal",5,0,5);
+ h3_ecal->GetXaxis()->SetTitle("Number of clusters"); h3_ecal->GetXaxis()->CenterTitle();
+ h3_ecal->GetXaxis()->SetNdivisions(405);h3_ecal->GetXaxis()->CenterLabels();
+ h3_ecal->SetLineColor(kRed);h3_ecal->SetLineWidth(2);
+
+ //Cut on theta*
+ TH1 *h3_hcal = new TH1D("h3_hcal","Number of reconstructed clusters in HCal",5,0,5);
+ h3_hcal->GetXaxis()->SetTitle("Number of clusters"); h3_hcal->GetXaxis()->CenterTitle();
+ h3_hcal->GetXaxis()->SetNdivisions(405);h3_hcal->GetXaxis()->CenterLabels();
+ h3_hcal->SetLineColor(kBlue);h3_hcal->SetLineWidth(2);
+
+ //Cut on theta*
+ TH1 *h4_hcal = new TH1D("h4_hcal","HCal cluster energy: Events w/ 1 HCal Cluster",100,-0.1,120);
+ h4_hcal->GetXaxis()->SetTitle("E [GeV]"); h4_hcal->GetXaxis()->CenterTitle();
+ h4_hcal->SetLineColor(kBlue);h4_hcal->SetLineWidth(2);
+
+ //Cut on generated neutron theta* + 1 HCal Cluster
+ TH1 *h1_neut_rec = new TH1D("h1_neut_rec","Reconstructed Neutron Energy",100,-0.1,120);
+ h1_neut_rec->GetXaxis()->SetTitle("E [GeV]"); h1_neut_rec->GetXaxis()->CenterTitle();
+ h1_neut_rec->SetLineColor(kBlue);h1_neut_rec->SetLineWidth(2);
+
+ //Cut on generated neutron theta* + 1 HCal Cluster
+ TH1 *h2_neut_rec = new TH2D("h2_neut_rec","Reconstructed Neutron local hit position at ZDC HCal front face",100,-200,200,100,-200,200);
+ h2_neut_rec->GetXaxis()->SetTitle("x [mm]"); h2_neut_rec->GetXaxis()->CenterTitle();
+ h2_neut_rec->GetYaxis()->SetTitle("y [mm]"); h2_neut_rec->GetYaxis()->CenterTitle();
+
+ //Read ROOT file
+ TFile* file = new TFile(inputfile.c_str());
+ TTree *tree = (TTree*) file->Get("events");
+
+ cout<<"Total number of events to analyze is "<<tree->GetEntries()<<endl;
+
+ //Create Array Reader
+ TTreeReader tr(tree);
+
+ TTreeReaderArray<int> gen_status(tr, "MCParticles.generatorStatus");
+ TTreeReaderArray<int> gen_pid(tr, "MCParticles.PDG");
+ TTreeReaderArray<float> gen_px(tr, "MCParticles.momentum.x");
+ TTreeReaderArray<float> gen_py(tr, "MCParticles.momentum.y");
+ TTreeReaderArray<float> gen_pz(tr, "MCParticles.momentum.z");
+ TTreeReaderArray<double> gen_mass(tr, "MCParticles.mass");
+ TTreeReaderArray<float> gen_charge(tr, "MCParticles.charge");
+ TTreeReaderArray<double> gen_vx(tr, "MCParticles.vertex.x");
+ TTreeReaderArray<double> gen_vy(tr, "MCParticles.vertex.y");
+ TTreeReaderArray<double> gen_vz(tr, "MCParticles.vertex.z");
+
+ //ZDC LYSO ECal
+ TTreeReaderArray<int> ecal_adc(tr,"EcalFarForwardZDCRawHits.amplitude");
+ TTreeReaderArray<float> ecal_hit_e(tr,"EcalFarForwardZDCRecHits.energy");
+ TTreeReaderArray<float> ecal_cluster_energy(tr, "EcalFarForwardZDCClusters.energy");
+
+ //ZDC SiPM-on-tile HCal
+ TTreeReaderArray<int> hcal_adc(tr,"HcalFarForwardZDCRawHits.amplitude");
+ TTreeReaderArray<float> hcal_hit_e(tr, "HcalFarForwardZDCRecHits.energy");
+ TTreeReaderArray<float> hcal_cluster_energy(tr, "HcalFarForwardZDCClusters.energy");
+
+ //Reconstructed neutron quantity
+ TTreeReaderArray<float> rec_neutron_energy(tr,"ReconstructedFarForwardZDCNeutrons.energy");
+ TTreeReaderArray<float> rec_neutron_px(tr,"ReconstructedFarForwardZDCNeutrons.momentum.x");
+ TTreeReaderArray<float> rec_neutron_py(tr,"ReconstructedFarForwardZDCNeutrons.momentum.y");
+ TTreeReaderArray<float> rec_neutron_pz(tr,"ReconstructedFarForwardZDCNeutrons.momentum.z");
+
+ //Other variables
+ int counter(0);
+
+ TLorentzVector neut_true; //True neutron in lab coordinates
+ TLorentzVector neut_true_rot; //True neutron wrt proton beam direction
+
+ float ecal_e_tot(0);
+ float hcal_e_tot(0);
+
+ //Loop over events
+ while (tr.Next()) {
+
+ if(counter%100==0) cout<<"Analyzing event "<<counter<<endl;
+ counter++;
+
+ //Reset variables
+ ecal_e_tot = 0;
+ hcal_e_tot = 0;
+
+ //Loop over generated particles, select primary neutron
+ for(int igen=0;igen<gen_status.GetSize();igen++){
+ if(gen_status[igen]==1 && gen_pid[igen]==2112){
+
+ neut_true.SetXYZM(gen_px[igen],gen_py[igen],gen_pz[igen],gen_mass[igen]);
+ h1_neut->Fill(neut_true.Theta()*TMath::RadToDeg(),neut_true.E());
+
+ //Wrt proton beam direction
+ neut_true_rot = neut_true;
+ neut_true_rot.RotateY(0.025);
+ h2_neut->Fill(neut_true_rot.Theta()*1000.,neut_true_rot.Phi()*TMath::RadToDeg()); //Theta* in mRad
+ h2a_neut->Fill(std::cos(neut_true_rot.Theta()),neut_true_rot.Phi()*TMath::RadToDeg());
+
+ }
+ } //End loop over generated particles
+
+ //Loop over ECal ADC hits (Raw hits may have different length than Rec hits due to zero supression)
+ for(int iadc=0;iadc<ecal_adc.GetSize();iadc++){
+ h1_ecal_adc->Fill(ecal_adc[iadc]);
+ }
+
+ //Loop over ECal hits
+ for(int ihit=0;ihit<ecal_hit_e.GetSize();ihit++){
+ ecal_e_tot += ecal_hit_e[ihit];
+ }
+
+ //Loop over ECal ADC hits (Raw hits may have different length than Rec hits due to zero supression)
+ for(int iadc=0;iadc<hcal_adc.GetSize();iadc++){
+ h1_hcal_adc->Fill(hcal_adc[iadc]);
+ }
+
+ //Loop over HCal hits
+ for(int ihit=0;ihit<hcal_hit_e.GetSize();ihit++){
+ hcal_e_tot += hcal_hit_e[ihit];
+ }
+
+ //ECal cluster size
+ int ecal_clus_size = ecal_cluster_energy.GetSize();
+
+ //HCal cluster size
+ int hcal_clus_size = hcal_cluster_energy.GetSize();
+
+ //Fill histograms for total energy and clusters
+ h1_ecal->Fill(ecal_e_tot);
+ h1_hcal->Fill(hcal_e_tot);
+
+ if( neut_true_rot.Theta()*1000. < 3.5 ){
+ h2_ecal->Fill(ecal_e_tot);
+ h2_hcal->Fill(hcal_e_tot);
+
+ h3_ecal->Fill(ecal_clus_size);
+ h3_hcal->Fill(hcal_clus_size);
+
+ //HCal cluster energy -- 1 cluster events
+ if(hcal_clus_size==1) h4_hcal->Fill(hcal_cluster_energy[0]);
+ }
+
+ //Reconstructed neutron(s)
+ for(int irec=0;irec<rec_neutron_energy.GetSize();irec++){
+ if(neut_true_rot.Theta()*1000. < 3.5 && hcal_clus_size==1){
+ h1_neut_rec->Fill(rec_neutron_energy[irec]);
+
+ TVector3 neut_rec(rec_neutron_px[irec],rec_neutron_py[irec],rec_neutron_pz[irec]);
+
+ //Wrt proton beam direction
+ TVector3 neut_rec_rot = neut_rec;
+ neut_rec_rot.RotateY(0.025);
+
+ //Projection of reconstructed neutron to front face of HCal
+ auto proj_x = 35.8 * 1000. * tan(neut_rec_rot.Theta()) * cos(neut_rec_rot.Phi()) ;
+ auto proj_y = 35.8 * 1000. * tan(neut_rec_rot.Theta()) * sin(neut_rec_rot.Phi()) ;
+ h2_neut_rec->Fill(proj_x,proj_y);
+ }
+ } //End loop over reconstructed neutrons
+
+ } //End loop over events
+
+ //Make plots
+ TCanvas *c1 = new TCanvas("c1");
+ h1_neut->Draw("colz");
+
+ TCanvas *c2 = new TCanvas("c2");
+ h2_neut->Draw("colz");
+
+ TCanvas *c2a = new TCanvas("c2a");
+ h2a_neut->Draw("colz");
+
+ TCanvas *c3a = new TCanvas("c3a");
+ c3a->SetLogy();
+ h1_ecal_adc->Draw();
+
+ TCanvas *c3b = new TCanvas("c3b");
+ c3b->SetLogy();
+ h1_hcal_adc->Draw();
+
+ TCanvas *c4 = new TCanvas("c4");
+ h1_ecal->Draw("");
+ h1_hcal->Draw("same");
+
+ TLegend *leg4 = new TLegend(0.6,0.6,0.85,0.8);
+ leg4->SetBorderSize(0);leg4->SetFillStyle(0);
+ leg4->AddEntry(h1_ecal,"Sum of digitized ZDC Ecal hit energies","l");
+ leg4->AddEntry(h1_hcal,"Sum of digitized ZDC Hcal hit energies","l");
+ leg4->Draw();
+
+ TCanvas *c5 = new TCanvas("c5");
+ h2_ecal->Draw("");
+ h2_hcal->Draw("same");
+
+ TLegend *leg5 = new TLegend(0.5,0.6,0.9,0.8);
+ leg5->SetBorderSize(0);leg5->SetFillStyle(0);
+ leg5->SetHeader("Require neutron #theta^{*} (wrt proton beam) < 3.5 mRad");
+ leg5->AddEntry(h1_ecal,"Sum of digitized ZDC Ecal hit energies","l");
+ leg5->AddEntry(h1_hcal,"Sum of digitized ZDC Hcal hit energies","l");
+ leg5->Draw();
+
+ TCanvas *c6a = new TCanvas("c6a");
+ h3_ecal->Draw();
+
+ TLegend *leg6 = new TLegend(0.5,0.7,0.9,0.9);
+ leg6->SetBorderSize(0);leg6->SetFillStyle(0);
+ leg6->SetHeader("Require neutron #theta^{*} (wrt proton beam) < 3.5 mRad");
+ leg6->Draw();
+
+ TCanvas *c6b = new TCanvas("c6b");
+ h3_hcal->Draw();
+ leg6->Draw();
+
+ TCanvas *c7 = new TCanvas("c7");
+ h4_hcal->Draw();
+
+ TLegend *leg7 = new TLegend(0.15,0.7,0.5,0.9);
+ leg7->SetBorderSize(0);leg7->SetFillStyle(0);
+ leg7->SetHeader("Require neutron #theta^{*} (wrt proton beam) < 3.5 mRad");
+ leg7->Draw();
+
+ TCanvas *c8 = new TCanvas("c8");
+ h1_neut_rec->Draw();
+
+ TLegend *leg8 = new TLegend(0.15,0.7,0.7,0.9);
+ leg8->SetBorderSize(0);leg8->SetFillStyle(0);
+ leg8->SetHeader("Generated neutron #theta^{*} < 3.5 mRad + 1 HCal cluster");
+ leg8->Draw();
+
+ TCanvas *c9 = new TCanvas("c9");
+ h2_neut_rec->Draw("colz");
+
+ TLatex *tex9 = new TLatex(-150,150,"Generated neutron #theta^{*} < 3.5 mRad + 1 HCal cluster");
+ tex9->SetTextSize(0.03);
+ tex9->Draw();
+
+ //Print plots to file
+ c1->Print(Form("%s[",outputfile.c_str()));
+ c1->Print(Form("%s",outputfile.c_str()));
+ c2->Print(Form("%s",outputfile.c_str()));
+ c2a->Print(Form("%s",outputfile.c_str()));
+ c3a->Print(Form("%s",outputfile.c_str()));
+ c3b->Print(Form("%s",outputfile.c_str()));
+ c4->Print(Form("%s",outputfile.c_str()));
+ c5->Print(Form("%s",outputfile.c_str()));
+ c6a->Print(Form("%s",outputfile.c_str()));
+ c6b->Print(Form("%s",outputfile.c_str()));
+ c7->Print(Form("%s",outputfile.c_str()));
+ c8->Print(Form("%s",outputfile.c_str()));
+ c9->Print(Form("%s",outputfile.c_str()));
+ c9->Print(Form("%s]",outputfile.c_str()));
+
+}
diff --git a/benchmarks/zdc_neutron/config.yml b/benchmarks/zdc_neutron/config.yml
new file mode 100644
index 00000000..715a9b20
--- /dev/null
+++ b/benchmarks/zdc_neutron/config.yml
@@ -0,0 +1,38 @@
+sim:zdc_neutron:
+ extends: .det_benchmark
+ stage: simulate
+ script:
+ - snakemake --cache --cores 1 sim_output/zdc_neutron/epic_craterlake/fwd_neutrons.edm4eic.root
+ retry:
+ max: 2
+ when:
+ - runner_system_failure
+
+bench:zdc_neutron:
+ extends: .det_benchmark
+ stage: benchmarks
+ needs:
+ - ["sim:zdc_neutron"]
+ script:
+ - snakemake --cores 1 results/zdc_neutron/epic_craterlake/fwd_neutrons_geant.pdf
+
+collect_results:zdc_neutron:
+ extends: .det_benchmark
+ stage: collect
+ needs:
+ - "bench:zdc_neutron"
+ script:
+ - ls -lrht
+ - mv results{,_save}/ # move results directory out of the way to preserve it
+ - snakemake --cores 1 --delete-all-output results/zdc_neutron/epic_craterlake/fwd_neutrons_geant.pdf
+ - mv results{_save,}/
+ # convert to png
+ - |
+ gs -sDEVICE=pngalpha -dUseCropBox -r144 \
+ -o 'results/zdc_neutron/epic_craterlake/geant_plots_%03d.png' \
+ results/zdc_neutron/epic_craterlake/fwd_neutrons_geant.pdf
+ - |
+ gs -sDEVICE=pngalpha -dUseCropBox -r144 \
+ -o 'results/zdc_neutron/epic_craterlake/recon_plots_%03d.png' \
+ results/zdc_neutron/epic_craterlake/fwd_neutrons_recon.pdf
+
diff --git a/benchmarks/zdc_neutron/gen_forward_neutrons.cxx b/benchmarks/zdc_neutron/gen_forward_neutrons.cxx
new file mode 100644
index 00000000..90a3eafe
--- /dev/null
+++ b/benchmarks/zdc_neutron/gen_forward_neutrons.cxx
@@ -0,0 +1,102 @@
+#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 <cmath>
+#include <math.h>
+#include <TMath.h>
+
+using namespace HepMC3;
+
+//Generate single neutrons in the forward detector region.
+void gen_forward_neutrons(int n_events = 10000, UInt_t seed = 0, const char* out_fname = "fwd_neutrons.hepmc")
+{
+
+ double theta_min = 0.0; //in mRad
+ double theta_max = 6.0; //in mRad
+ double cost_min = std::cos(theta_max/1000.) ; //Minimum value of cos(theta)
+ double cost_max = std::cos(theta_min/1000.) ; //Maximum value of cos(theta)
+
+ double p_min = 100.; //in GeV/c
+ double p_max = 100.; //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();
+ TParticlePDG *particle = pdg->GetParticle("neutron");
+ 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 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 p = r1->Uniform(p_min,p_max); //Uniform momentum generation
+ Double_t phi = r1->Uniform(0.0, 2.0 * M_PI); //Uniform phi generation
+
+ Double_t cost = r1->Uniform(cost_min,cost_max); //Uniform cos(theta) generation
+ Double_t th = std::acos(cost);
+
+ Double_t px = p * std::cos(phi) * std::sin(th);
+ Double_t py = p * std::sin(phi) * std::sin(th);
+ Double_t pz = p * std::cos(th);
+ Double_t E = sqrt(p*p + mass*mass);
+
+ //Rotate to lab coordinate system
+ TVector3 pvec(px,py,pz);
+ double cross_angle = -25./1000.; //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
+ GenParticlePtr p3 = std::make_shared<GenParticle>(
+ FourVector(pvec.X(), pvec.Y(), pvec.Z(), E), pdgID, 1 );
+
+ 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 % 10000 == 0) {
+ std::cout << "Event: " << events_parsed << std::endl;
+ }
+ evt.clear();
+ }
+ hepmc_output.close();
+ std::cout << "Events parsed and written: " << events_parsed << std::endl;
+}
--
GitLab