Working on updating particle gun and hits reader for B0.

benchmarks/b0_tracker/scripts/b0_tracker_hits.cxx

@@ -40,21 +40,21 @@ void b0_tracker_hits(const char* fname = "./sim_output/sim_forward_protons.root"
     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 d1 = d0.Define("hits_theta", hits_theta, {"B0TrackerHits"});
 
-  auto h1 = d1.Histo1D({"h1", "hits_eta", 100, 0,20}, "hits_eta");
+  auto h1 = d1.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");
   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);

benchmarks/b0_tracker/scripts/gen_forward_protons.cxx

@@ -20,8 +20,13 @@ using namespace HepMC3;
 void gen_forward_protons(int n_events = 100, 
                      const char* out_fname = "forward_protons.hepmc")
 {
-  double cos_theta_min = std::cos(1.0*(M_PI/180.0));
-  double cos_theta_max = std::cos(0.0*(M_PI/180.0));
+  
+  // generate protons in B0 acceptance - roughly 5 - 20 mrad
+  double cos_theta_min = std::cos(0.005*(M_PI/180.0));
+  double cos_theta_max = std::cos(0.020*(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;
 
@@ -39,12 +44,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(1.0, 10.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);