Resolve "Check ZDC particle gun acceptance"

benchmarks/zdc/run_zdc_particles.sh

@@ -38,16 +38,18 @@ if [[ ! -n  "${JUGGLER_DETECTOR}" ]] ; then
 fi
 
 if [[ ! -n  "${JUGGLER_N_EVENTS}" ]] ; then 
-  export JUGGLER_N_EVENTS=1000
+  export JUGGLER_N_EVENTS=5000
 fi
 
+#a place where I can run my events when I want to. - AMJ
+#export JUGGLER_N_EVENTS=1000
 
 if [[ ! -n  "${E_start}" ]] ; then
-  export E_start=5.0
+  export E_start=125.0
 fi
 
 if [[ ! -n  "${E_end}" ]] ; then
-  export E_end=5.0
+  export E_end=145.0
 fi
 
 export JUGGLER_FILE_NAME_TAG="zdc_${PARTICLE}"

benchmarks/zdc/scripts/analysis_zdc_particles.cxx

@@ -8,6 +8,8 @@
 
 #include "edm4hep/MCParticleCollection.h"
 #include "edm4hep/SimCalorimeterHitCollection.h"
+#include "eicd/CalorimeterHitCollection.h"
+#include "eicd/CalorimeterHitData.h"
 
 #include "TCanvas.h"
 #include "TStyle.h"
@@ -60,6 +62,30 @@ void analysis_zdc_particles(
     return sampled / thrown;
   };
 
+  // Hit position X
+  auto hit_x_position = [&] (const std::vector<edm4hep::SimCalorimeterHitData>& hits) {
+	std::vector<double> result;
+	for(const auto& h: hits)
+		result.push_back(h.position.x); //mm
+  return result; 
+  };
+
+  // Hit position Y
+  auto hit_y_position = [&] (const std::vector<edm4hep::SimCalorimeterHitData>& hits) {
+        std::vector<double> result;
+        for(const auto& h: hits)
+		result.push_back(h.position.y); //mm
+  return result;
+  };
+
+  // Hit position Z
+  auto hit_z_position = [&] (const std::vector<edm4hep::SimCalorimeterHitData>& hits) {
+        std::vector<double> result;
+        for(const auto& h: hits)
+		result.push_back(h.position.z); //mm
+  return result;
+  };
+
   // Define variables
   auto d1 = d0
     .Define("Ethr", Ethr, {"MCParticles"})
@@ -69,6 +95,8 @@ void analysis_zdc_particles(
     .Define("nhits_Hcal", nhits, {"ZDCHcalHits"})
     .Define("Esim_Hcal", Esim, {"ZDCHcalHits"})
     .Define("fsam_Hcal", fsam, {"Esim_Hcal", "Ethr"})
+    .Define("hit_x_position", hit_x_position, {"ZDCHcalHits"})
+	  .Define("hit_y_position", hit_y_position, {"ZDCHcalHits"})
   ;
 
   // Define Histograms
@@ -96,6 +124,9 @@ void analysis_zdc_particles(
       {"hfsam_Hcal", "Hcal Sampling Fraction; Hcal Sampling Fraction; Events", 100, 0.0, 0.1},
       "fsam_Hcal");
 
+  auto hXY_HitMap_Hcal = d1.Histo2D({"hXY_HitMap_Hcal", "hit position Y vs. X histogram; hit position X [mm]; hit position Y [mm]", 8,-1300,-500, 8, -400, 400}, "hit_x_position", "hit_y_position");
+
+
   // Event Counts
   auto nevents_thrown = d1.Count();
   std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n";
@@ -187,4 +218,10 @@ void analysis_zdc_particles(
     c4->SaveAs(TString(results_path + "/zdc_fsam_Hcal.png"));
     c4->SaveAs(TString(results_path + "/zdc_fsam_Hcal.pdf"));
   }
+  {
+    TCanvas* c5 = new TCanvas("c5", "c5", 600, 500);
+    hXY_HitMap_Hcal->Draw("COLZ");
+    c5->SaveAs(TString(results_path + "/zdc_xy_map_Hcal.png"));
+    c5->SaveAs(TString(results_path + "/zdc_xy_map_Hcal.pdf"));
+  }
 }

benchmarks/zdc/scripts/gen_zdc_particles.cxx

@@ -21,7 +21,7 @@
 
 using namespace HepMC3;
 
-void gen_zdc_particles(int n_events = 1e6, const std::string& particle = "neutron", double p_start = 0.0, double p_end = 30.0, const std::string& out_fname = "./data/zdc_neutrons.hepmc") {
+void gen_zdc_particles(int n_events = 1e6, const std::string& particle = "neutron", double p_start = 125.0, double p_end = 145.0, const std::string& out_fname = "./data/zdc_neutrons.hepmc") {
   WriterAscii hepmc_output(out_fname);
   int events_parsed = 0;
   GenEvent evt(Units::GEV, Units::MM);
@@ -36,8 +36,8 @@ void gen_zdc_particles(int n_events = 1e6, const std::string& particle = "neutro
   // detector
   // https://indico.bnl.gov/event/7449/contributions/35966/attachments/27177/41430/EIC-DWG-Calo-03192020.pdf
   // See a figure on slide 26
-  double cos_theta_min = std::cos(M_PI * (0.0 / 180.0));
-  double cos_theta_max = std::cos(M_PI * (5.0 / 180.0));
+  double cos_theta_min = std::cos(0.0);
+  double cos_theta_max = std::cos(0.005);
 
   for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
     // FourVector(px,py,pz,e,pdgid,status)

benchmarks/zdc/simple_info_plot_histograms.cxx

@@ -150,6 +150,8 @@ void simple_info_plot_histograms(const char* fname = "sim_output/output_zdc_phot
   auto h5 = d1.Histo1D({"h5", "detector ID; detector ID; Events", 3,-0.5,2.5}, "detID");
   auto h6 = d1.Histo1D({"h6", "volume ID; volume ID; Events", 100,0,50000000}, "volID");
 
+  auto h7 = d1.Histo2D({"h7", "hit position Y vs. X histogram; hit position X [mm]; hit position Y [mm]", 100,-30,30, 100, -30, 30}, "hit_x_position", "hit_y_position");
+
   auto n0 = d1.Filter([](int n){ return (n>0); },{"nhits"}).Count();
 
   d1.Snapshot("info_EVENT","sim_output/info_zdc_photons.edm4hep.root");
@@ -209,6 +211,10 @@ void simple_info_plot_histograms(const char* fname = "sim_output/output_zdc_phot
   h6->DrawClone();
   c4->SaveAs("sim_output/detID_volID_histo_zdc_photons.png");
 
+  TCanvas *c5 = new TCanvas("c5","c5", 500, 500);
+  h7->Draw("COLZ");
+  c5->SaveAs("results/far_forward/zdc/zdc_x_y_hit_map.png");
+
   if(*n0<1) {
     std::quick_exit(1);
   }