Added benchmark for single neutron in the ZDC (#99)

.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"

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"

benchmarks/zdc_neutron/README.md

+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)
+
+

benchmarks/zdc_neutron/Snakefile

+# 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}")'
+"""

benchmarks/zdc_neutron/analysis/fwd_neutrons_geant.C

+//------------------
+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()));
+
+}

benchmarks/zdc_neutron/analysis/fwd_neutrons_recon.C

+//------------------
+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())); 
+
+}

benchmarks/zdc_neutron/config.yml

+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
+

benchmarks/zdc_neutron/gen_forward_neutrons.cxx

+#include "HepMC3/GenEvent.h"
+#include "HepMC3/ReaderAscii.h"
+#include "HepMC3/WriterAscii.h"
+#include "HepMC3/Print.h"
+
+#include "TRandom3.h"
+#include "TVector3.h"
+
+#include <TDatabasePDG.h>
+#include <TParticlePDG.h>
+
+#include <iostream>
+#include <random>
+#include <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;
+}