Final commit, fixed issue 3

benchmarks/barrel_ecal/scripts/emcal_barrel_pi0_analysis.cxx

@@ -30,7 +30,7 @@ using namespace ROOT::VecOps;
 void emcal_barrel_pi0_analysis(const char* input_fname = "sim_output/sim_emcal_barrel_pi0.root")
 {
   //input_fname = "../../../temp_pi0.root";
-  //input_fname = "../sim_output/sim_emcal_barrel_pi0.root";
+  input_fname = "../sim_output/sim_emcal_barrel_pi0.root";
   // Setting for graphs
   gROOT->SetStyle("Plain");
   gStyle->SetOptFit(1);
@@ -50,6 +50,7 @@ void emcal_barrel_pi0_analysis(const char* input_fname = "sim_output/sim_emcal_b
   double samp_frac = 1.0;
 
   // Thrown Energy [GeV]
+  double meanE = 5; // Calculated later
   auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
     return TMath::Sqrt(input[2].psx*input[2].psx + input[2].psy*input[2].psy + input[2].psz*input[2].psz + input[2].mass*input[2].mass);
   };
@@ -110,8 +111,9 @@ void emcal_barrel_pi0_analysis(const char* input_fname = "sim_output/sim_emcal_b
   auto hpid   = d1.Histo1D({"hpid",   "PID; PID; Count",                                   100,  -220,   220}, "pid");
   auto hdau   = d1.Histo1D({"hdau",   "Number of Daughters; Number of Daughters; Count",   10,   0,      10},  "dau");
 
-  // Set sampling Fraction, ideally this will be taken from a json file 
-  samp_frac = hfsam -> GetMean();
+  // Gather Sampling fraction and mean Energy, ideally this will be taken from a json file 
+  samp_frac = hfsam->GetMean();
+  meanE     = hEthr->GetMean(); 
 
   // Event Counts
   auto nevents_thrown = d1.Count();
@@ -183,8 +185,8 @@ void emcal_barrel_pi0_analysis(const char* input_fname = "sim_output/sim_emcal_b
   // Taken from : Initial considerations for EMCal of the EIC detector by A. Bazilevsky
   // sigma_E / E = 12% / E^0.5 convoluted with 2%
   // sigma_E / E = [ (0.12/E^0.5)^2 + 0.02^2]^0.5, with E in [GeV]
-  double thrown_energy    = hEthr->GetMean(); // Current thrown energy, will need to grab from json file
-  double resolutionTarget = TMath::Sqrt(0.12 * 0.12 / thrown_energy + 0.02 * 0.02);
+  
+  double resolutionTarget = TMath::Sqrt(0.12 * 0.12 / meanE + 0.02 * 0.02);
 
   eic::util::Test pi0_energy_resolution{
    {{"name", fmt::format("{}_energy_resolution", test_tag)},
@@ -200,7 +202,7 @@ void emcal_barrel_pi0_analysis(const char* input_fname = "sim_output/sim_emcal_b
   auto hdE_rel      = (d1.Histo1D({"hdE_rel",  "dE Relative; dE Relative; Events", 100, -3.0, 3.0}, "dE_rel"))->DrawCopy();
   hdE->Fit("gaus", "", "", -3.0,  3.0);
   double* res       = hdE->GetFunction("gaus")->GetParameters();
-  double sigmaOverE = res[2] / thrown_energy;
+  double sigmaOverE = res[2] / meanE;
 
   //Pass/Fail
   sigmaOverE <= resolutionTarget ? pi0_energy_resolution.pass(sigmaOverE) : pi0_energy_resolution.fail(sigmaOverE);