Draft: Resolve "Add Roman Pots Reconstruction Matrix"

Closes #70

benchmarks/far_forward/scripts/gen_far_forward_protons.cxx

@@ -14,11 +14,17 @@ using namespace HepMC3;
 /** Generate protons in the far forward region.
  *  This is for testing detectors in the B0 tracker.
  */
-void gen_far_forward_protons(int n_events = 100, 
+void gen_far_forward_protons(int n_events = 1000, 
                      const char* out_fname = "far_forward_protons.hepmc")
 {
-  double cos_theta_min = std::cos(0.3*(M_PI/180.0));
-  double cos_theta_max = std::cos(2.0*(M_PI/180.0));
+
+  double crossingAngle = -0.025; //-25 mrad crossing angle
+
+  double minMomentum = 273.0; //GeV
+  double maxMomentum = 275.0; //GeV
+  
+  double theta_min = 0.0; //we need MRAD for the units - syntax is degrees*Pi/180
+  double theta_max = 0.005; //1.15 degrees = 20mrad FF region is 0.0 to 20.0 mrad
 
   WriterAscii hepmc_output(out_fname);
   int events_parsed = 0;
@@ -39,10 +45,10 @@ void gen_far_forward_protons(int n_events = 100,
         FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
 
     // Define momentum
-    Double_t p     = r1->Uniform(1.0, 10.0);
+    Double_t p     = r1->Uniform(minMomentum, maxMomentum);
     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);
+    Double_t th    = r1->Uniform(theta_min, theta_max);
+    //Double_t th    = std::acos(costh);
     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);
@@ -52,6 +58,15 @@ void gen_far_forward_protons(int n_events = 100,
 
     //std::cout << std::sqrt(px*px + py*py + pz*pz) - p << " is zero? \n";
 
+    //Rotate vector by -25mrad crossing angle
+
+    Double_t s = std::sin(crossingAngle);
+    Double_t c = std::cos(crossingAngle);
+    Double_t zz = pz;
+    Double_t xx = px;
+    pz = c*zz - s*xx;
+    px = s*zz + c*xx;
+
     // type 1 is final state
     // pdgid 11 - proton 938.272 MeV/c^2
     GenParticlePtr p3 = std::make_shared<GenParticle>(