Update Analysis Script

Adding analysis script

• Put a cut on the energy difference, | E_{\mathrm{cluster}} - E_{\mathrm{thrown}} | < \delta_{E} or something similar.
• acceptance as a function of P_{T} , \eta, etc., for bins in different bins in E_{\mathrm{thrown}}
• resolution
• mass reconstruction

ecal/scripts/emcal_electrons.cxx

@@ -19,12 +19,8 @@
 
 using namespace HepMC3;
 
-void emcal_electrons(int n_events = 1e2, double e_start = 1.0, double e_end = 1.0,
-                     const char* out_fname = "./data/emcal_electron_0GeVto30GeV_100kEvt.hepmc")
+void emcal_electrons(int n_events = 1e6, double e_start = 0.0, double e_end = 30.0, const char* out_fname = "./data/emcal_electron_0GeVto30GeV_1000kEvt.hepmc")
 {
-  double cos_theta_min = std::cos(M_PI * (155.0 / 180.0));
-  double cos_theta_max = std::cos(M_PI);
-
   WriterAscii hepmc_output(out_fname);
   int events_parsed = 0;
   GenEvent evt(Units::GEV, Units::MM);
@@ -32,6 +28,12 @@ void emcal_electrons(int n_events = 1e2, double e_start = 1.0, double e_end = 1.
   // Random number generator
   TRandom *r1 = new TRandom();
 
+  // Constraining the solid angle, but larger than that subtended by the 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 * (155.0 / 180.0));
+  double cos_theta_max = std::cos(M_PI);
+
   for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
     // FourVector(px,py,pz,e,pdgid,status)
     // type 4 is beam
@@ -44,18 +46,18 @@ void emcal_electrons(int n_events = 1e2, double e_start = 1.0, double e_end = 1.
         FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
 
     // Define momentum
-    Double_t p     = r1->Uniform(0.0, 30.0);
-    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 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);
+    // Momentum starting at 0 GeV/c to 30 GeV/c
+    Double_t p        = r1->Uniform(e_start, e_end);
+    Double_t phi      = r1->Uniform(0.0, 2.0 * M_PI);
+    Double_t costheta = r1->Uniform(cos_theta_min, cos_theta_max);
+    Double_t theta    = std::acos(costheta);
+    Double_t px       = p * std::cos(phi) * std::sin(theta);
+    Double_t py       = p * std::sin(phi) * std::sin(theta);
+    Double_t pz       = p * std::cos(theta);
+
     // 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";
+    //r1->Sphere(px, py, pz, p);
 
     // type 1 is final state
     // pdgid 11 - electron 0.510 MeV/c^2