From af5f7c132916f39f515a2884d5e6b8178c538a78 Mon Sep 17 00:00:00 2001
From: Whitney Armstrong <warmstrong@anl.gov>
Date: Thu, 26 Aug 2021 21:22:31 -0500
Subject: [PATCH] Squashed commit of the following:

commit cbee19dbe412292d4f2209a567f1743fb3a8cd68
Author: Alexander Jentsch <ajentsch@bnl.gov>
Date:   Mon Aug 16 15:23:55 2021 -0400

    This is really not the place to do this test. Should be done with final output file after reco.

commit 582f72b6f0aeda542a1ab3672fcdc2b460207aab
Author: Alexander Jentsch <ajentsch@bnl.gov>
Date:   Mon Aug 16 14:52:58 2021 -0400

    Something got messed up in the rebase. Fixed and now testing original functionality.

commit 510cb6cb5290823c7bb518e9b06d30ca0c0b7d50
Merge: 0d80f47 548e97b
Author: Alexander Jentsch <ajentsch@bnl.gov>
Date:   Mon Aug 16 14:40:25 2021 -0400

    Merge branch '50-b0tracker-acceptance-tests' of https://eicweb.phy.anl.gov/EIC/benchmarks/detector_benchmarks into 50-b0tracker-acceptance-tests
    Still trying to add MC vs. accepted information - rebase issue.

commit 0d80f4731fed7d490a2ec502dd6a4bb97880ecbd
Author: Alexander Jentsch <ajentsch@bnl.gov>
Date:   Mon Aug 16 14:37:54 2021 -0400

    Trying to grab MC data so I can make some very simple acceptance plots.

commit 548e97b5b5705204a320a076b0c6b7b125aa0cef
Author: Alexander Jentsch <ajentsch@bnl.gov>
Date:   Wed Aug 11 11:42:19 2021 -0400

    Add crossing angle to gen_particle script.

commit d60f1fa503894f2b4a43ab22573cfe01621d212a
Author: Alexander Jentsch <ajentsch@bnl.gov>
Date:   Wed Aug 11 11:18:21 2021 -0400

    Things seem to work now. Need more events to test operation.

commit 21bab2a4a768d23667782e94c314140d8517b0e5
Author: Alexander Jentsch <ajentsch@bnl.gov>
Date:   Wed Aug 11 11:09:30 2021 -0400

    Another issue.

commit cc86cb9ae1399ccbf5b336716f29ecbec42fad35
Author: Alexander Jentsch <ajentsch@bnl.gov>
Date:   Wed Aug 11 11:03:21 2021 -0400

    Fixed typo(s).

commit e24c272856c28f3675496f0635b0d8694e6d7b6c
Author: Alexander Jentsch <ajentsch@bnl.gov>
Date:   Wed Aug 11 10:49:10 2021 -0400

    Trying to add occupancy plots for the B0.

commit 726efc4f777fa3b09dc1d1c608fa912b1766301f
Author: Alexander Jentsch <ajentsch@bnl.gov>
Date:   Wed Aug 11 09:11:30 2021 -0400

    Fixed typo - still getting used to syntax.

commit a2e0d2955964a9159217a1869b82a3ac094e7ee7
Author: Alexander Jentsch <ajentsch@bnl.gov>
Date:   Wed Aug 11 09:06:41 2021 -0400

    Working on updating particle gun and hits reader for B0.

commit 0bcbf1637ce79168cd80b7351cbe0b14da293830
Author: Alexander Jentsch <ajentsch@bnl.gov>
Date:   Wed Aug 11 11:42:19 2021 -0400

    Add crossing angle to gen_particle script.

commit 2b677af9efff50117e09a96fb579180d24ca8756
Author: Alexander Jentsch <ajentsch@bnl.gov>
Date:   Wed Aug 11 11:18:21 2021 -0400

    Things seem to work now. Need more events to test operation.

commit cab2b05c5aadbc7c96795481623883d55fc80b9a
Author: Alexander Jentsch <ajentsch@bnl.gov>
Date:   Wed Aug 11 11:09:30 2021 -0400

    Another issue.

commit 21eeb20d62eb5d051261053f72568eeb448b0c2f
Author: Alexander Jentsch <ajentsch@bnl.gov>
Date:   Wed Aug 11 11:03:21 2021 -0400

    Fixed typo(s).

commit d3a75da26f10f2a4d8684c9116fe5fbc422fd3f7
Author: Alexander Jentsch <ajentsch@bnl.gov>
Date:   Wed Aug 11 10:49:10 2021 -0400

    Trying to add occupancy plots for the B0.

commit 8984a51fec9bf142b15633601502a212abbbb5ba
Author: Alexander Jentsch <ajentsch@bnl.gov>
Date:   Wed Aug 11 09:11:30 2021 -0400

    Fixed typo - still getting used to syntax.

commit cf469f0ca22a2eae299df84c16d7091a02a5aa81
Author: Alexander Jentsch <ajentsch@bnl.gov>
Date:   Wed Aug 11 09:06:41 2021 -0400

    Working on updating particle gun and hits reader for B0.
---
 benchmarks/b0_tracker/forward_protons.sh      |  2 +-
 .../b0_tracker/scripts/b0_tracker_hits.cxx    | 58 ++++++++++++++++---
 .../scripts/gen_forward_protons.cxx           | 23 ++++----
 3 files changed, 62 insertions(+), 21 deletions(-)

diff --git a/benchmarks/b0_tracker/forward_protons.sh b/benchmarks/b0_tracker/forward_protons.sh
index d78b81b7..e5d3572b 100755
--- a/benchmarks/b0_tracker/forward_protons.sh
+++ b/benchmarks/b0_tracker/forward_protons.sh
@@ -7,7 +7,7 @@ if [[ ! -n  "${JUGGLER_DETECTOR}" ]] ; then
 fi
 
 if [[ ! -n  "${JUGGLER_N_EVENTS}" ]] ; then 
-  export JUGGLER_N_EVENTS=100
+  export JUGGLER_N_EVENTS=1000
 fi
 
 export FILE_NAME_TAG="forward_protons"
diff --git a/benchmarks/b0_tracker/scripts/b0_tracker_hits.cxx b/benchmarks/b0_tracker/scripts/b0_tracker_hits.cxx
index b4984ada..c6dd692c 100644
--- a/benchmarks/b0_tracker/scripts/b0_tracker_hits.cxx
+++ b/benchmarks/b0_tracker/scripts/b0_tracker_hits.cxx
@@ -36,25 +36,67 @@ void b0_tracker_hits(const char* fname = "./sim_output/sim_forward_protons.root"
 
   ROOT::RDataFrame d0(*t);
 
-  auto hits_eta = [&](const std::vector<dd4pod::TrackerHitData>& hits) {
+  auto hits_theta = [&](const std::vector<dd4pod::TrackerHitData>& hits) {
     std::vector<double> result;
     for (const auto& h : hits) {
       ROOT::Math::XYZVector vec(h.position.x,h.position.y,h.position.z);
-      result.push_back(vec.eta());
-      std::cout << vec.eta() << "\n";
+      result.push_back(1000*vec.theta());
+      std::cout << 1000*vec.theta() << "\n";
     }
     return result;
   };
 
-  auto d1 = d0.Define("hits_eta", hits_eta, {"B0TrackerHits"});
+  auto local_position = [&](const std::vector<dd4pod::TrackerHitData>& hits) {
+    std::vector<std::array<double, 2>> result;
+    for (const auto& h : hits) {
+
+      auto pos0 = (h.position);
+	  result.push_back({pos0.x , pos0.y});
+    }
+    return result;
+  };
+
+
+  auto x_pos = [&](const std::vector<std::array<double, 2>>& xypos) {
+    std::vector<double> result;
+    for (const auto& h : xypos) {
+
+		result.push_back(h.at(0));
+    }
+    return result;
+  };
 
-  auto h1 = d1.Histo1D({"h1", "hits_eta", 100, 0,20}, "hits_eta");
+  auto y_pos = [&](const std::vector<std::array<double, 2>>& xypos) {
+    std::vector<double> result;
+    for (const auto& h : xypos) {
+
+        result.push_back(h.at(1));
+    }
+    return result;
+  };
+
+
+  auto d1 = d0.Define("nhits", hits_theta, {"B0TrackerHits"})
+			  .Define("xy_hit_pos", local_position, {"B0TrackerHits"})
+              .Define("x_pos", x_pos, {"xy_hit_pos"})
+              .Define("y_pos", y_pos, {"xy_hit_pos"});
+
+  auto h_local_pos = d1.Histo2D({"h_local_pos", ";x [mm]; y [mm] ", 100,  -100.0, -200.0, 100, -100.0, 100.0}, "x_pos", "y_pos");
+
+  auto d2 = d0.Define("hits_theta", hits_theta, {"B0TrackerHits"});
+
+  auto h1 = d2.Histo1D({"h1", "hits_theta", 100, 0,20}, "hits_theta");
   TCanvas* c = new TCanvas();
   h1->DrawCopy();
-  c->SaveAs("results/b0_tracker_hits_eta.png");
-  c->SaveAs("results/b0_tracker_hits_eta.pdf");
+  c->SaveAs("results/b0_tracker_hits_theta.png");
+  c->SaveAs("results/b0_tracker_hits_theta.pdf");
+  
+  h_local_pos->DrawCopy("colz");
+  c->SaveAs("results/b0_tracker_hits_occupancy_disk_1.png");
+  c->SaveAs("results/b0_tracker_hits_occupancy_disk_1.pdf");
+
   auto n1 = h1->GetMean();
-  std::cout << "Pseudorapidity of hits: " << n1 << std::endl;
+  std::cout << "Polar angle of hits: " << n1 << std::endl;
 
   //if (n1 < 5) {
   //        std::quick_exit(1);
diff --git a/benchmarks/b0_tracker/scripts/gen_forward_protons.cxx b/benchmarks/b0_tracker/scripts/gen_forward_protons.cxx
index 6cb9bd50..9e1a3405 100644
--- a/benchmarks/b0_tracker/scripts/gen_forward_protons.cxx
+++ b/benchmarks/b0_tracker/scripts/gen_forward_protons.cxx
@@ -20,12 +20,19 @@ using namespace HepMC3;
 /** Generate electrons in the central region.
  *  This is for testing detectors in the "barrel" region.
  */
-void gen_forward_protons(int n_events = 100, 
+void gen_forward_protons(int n_events = 1000, 
                      const char* out_fname = "forward_protons.hepmc")
 {
-  double cos_theta_min = std::cos(0.5*(M_PI/180.0));
+
+  double crossingAngle = -0.025; //radiansc
+
+  // generate protons in B0 acceptance - roughly 5 - 20 mrad
+  double cos_theta_min = std::cos(0.02);
   double cos_theta_max = std::cos(0.0*(M_PI/180.0));
 
+  double partEnergyMin = 270.0; // xL 0.98
+  double partEnergyMax = 275.0; // top beam energy
+
   const double M_p = common_bench::particleMap.at(2212).mass;
 
   WriterAscii hepmc_output(out_fname);
@@ -42,12 +49,12 @@ void gen_forward_protons(int n_events = 100,
     // pdgid 111 - pi0
     // pdgid 2212 - proton
     GenParticlePtr p1 =
-        std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 2212, 4);
+        std::make_shared<GenParticle>(FourVector(0.0, 0.0, partEnergyMax, partEnergyMax), 2212, 4);
     GenParticlePtr p2 = std::make_shared<GenParticle>(
         FourVector(0.0, 0.0, 0.0, M_p), 2212, 4);
 
     // Define momentum
-    Double_t p     = r1->Uniform(200.0, 275.0);
+    Double_t p     = r1->Uniform(partEnergyMin, partEnergyMax);
     Double_t phi   = r1->Uniform(0.0, 2.0 * M_PI);
     Double_t costh = r1->Uniform(cos_theta_min, cos_theta_max);
     Double_t th    = std::acos(costh);
@@ -55,20 +62,12 @@ void gen_forward_protons(int n_events = 100,
     Double_t py    = p * std::sin(phi) * std::sin(th);
     Double_t pz    = p * std::cos(th);
 
-
     ROOT::Math::XYZVector p0 = {px,py,pz};
 
     //ROOT::Math::Rotation3D r = (-0.025);
     ROOT::Math::RotationY r(-0.025);
     auto p_rot = r*p0;
 
-
-    // Generates random vectors, uniformly distributed over the surface of a
-    // sphere of given radius, in this case momentum.
-    // r1->Sphere(px, py, pz, p);
-
-    //std::cout << std::sqrt(px*px + py*py + pz*pz) - p << " is zero? \n";
-
     // type 1 is final state
     // pdgid 11 - electron 0.510 MeV/c^2
     GenParticlePtr p3 = std::make_shared<GenParticle>(
-- 
GitLab