Ongoing dis_electron Q2 analysis. An issue with Q2 being basically zero for the electrons.

So I have made the branch and got the functions necessary, but there is an issue that the electrons Q2 is essentially zero. So the sorted particles within each event with the highest momenta are basically electrons that "missed" the target. Maybe I should sort by highest Q2 and rerun.

benchmarks/dis/analysis/dis_electrons.cxx

@@ -6,7 +6,12 @@
 #include <util.h>
 
 #include "ROOT/RDataFrame.hxx"
+#include <TFitResult.h>
+#include <TH1D.h>
+#include <TRandom3.h>
+#include <algorithm>
 #include <cmath>
+#include <eicd/ReconstructedParticleData.h>
 #include <fmt/color.h>
 #include <fmt/core.h>
 #include <fstream>
@@ -14,120 +19,117 @@
 #include <nlohmann/json.hpp>
 #include <string>
 #include <vector>
-#include <eicd/ReconstructedParticleData.h>
-#include <TH1D.h>
-#include <TFitResult.h>
-#include <TRandom3.h>
-#include <algorithm>
-#include <utility>
 
+// Beam energies
+double eBeamEnergy;
+double pBeamEnergy;
+
+// Gaussian variables
+double gausMean;
+double gausSTD;
+
+// Makes Resolution plots from the input histograms
+void plot_resolutions(ROOT::RDF::RResultPtr<TH1D> h_res, ROOT::RDF::RResultPtr<TH1D> h_rel_res,
+                      std::string varName, std::string detName, std::string outputPREFIX)
+{
+  TCanvas c{"canvas", "canvas", 1200, 1200};
+  c.cd();
+  // gPad->SetLogx(false);
+  gPad->SetLogy(true);
+  auto& h1 = *h_res;
+  auto& h2 = *h_rel_res;
+  // histogram style
+  h1.SetLineColor(plot::kMpBlue);
+  h1.SetLineWidth(2);
+  h2.SetLineColor(plot::kMpOrange);
+  h2.SetLineWidth(2);
+  // axes
+  h1.GetXaxis()->CenterTitle();
+  h1.GetYaxis()->CenterTitle();
+  h2.GetXaxis()->CenterTitle();
+  h2.GetYaxis()->CenterTitle();
+  // draw h1
+  h1.DrawClone("hist");
+  plot::draw_label(eBeamEnergy, pBeamEnergy, detName);
+  TText* tptr1;
+  auto   t1 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
+  t1->SetFillColorAlpha(kWhite, 0);
+  t1->SetTextFont(43);
+  t1->SetTextSize(25);
+  tptr1 = t1->AddText((varName + " Resolution").c_str());
+  tptr1->SetTextColor(plot::kMpBlue);
+  t1->Draw();
+  c.Print(fmt::format("{}{}_resolution.png", outputPREFIX, varName).c_str());
+  t1->Clear();
+  // draw h2
+  h2.DrawClone("hist");
+  plot::draw_label(eBeamEnergy, pBeamEnergy, detName);
+  t1->SetFillColorAlpha(kWhite, 0);
+  t1->SetTextFont(43);
+  t1->SetTextSize(25);
+  tptr1 = t1->AddText((varName + " Relative Resolution").c_str());
+  tptr1->SetTextColor(plot::kMpOrange);
+  t1->Draw();
+  c.Print(fmt::format("{}{}_rel_resolution.png", outputPREFIX, varName).c_str());
+}
 
 // Get a vector of 4-momenta from the reconstructed data.
-inline auto momenta_from_reconstruction(const std::vector<eic::ReconstructedParticleData>& parts)
+inline auto
+momenta_from_reconstruction_pid(const std::vector<eic::ReconstructedParticleData>& parts)
 {
-  std::vector<ROOT::Math::PxPyPzEVector> momenta{parts.size()};
+  std::vector<std::pair<ROOT::Math::PxPyPzEVector, int>> momenta(parts.size());
   // transform our reconstructed particle data into 4-momenta
   std::transform(parts.begin(), parts.end(), momenta.begin(), [](const auto& part) {
-    return ROOT::Math::PxPyPzEVector{part.p.x, part.p.y, part.p.z, part.energy};
+    return std::make_pair(ROOT::Math::PxPyPzEVector{part.p.x, part.p.y, part.p.z, part.energy},
+                          part.pid);
   });
   return momenta;
 }
 
-// Convert PxPyPzMVector to PxPyPzEVector
-inline auto convertMtoE(const std::vector<ROOT::Math::PxPyPzMVector>& mom)
+// Get a vector of 4-momenta from the simulated data.
+inline auto momenta_from_simulation_pid(const std::vector<dd4pod::Geant4ParticleData>& parts)
 {
-  std::vector<ROOT::Math::PxPyPzEVector> momenta{mom.size()};
-  std::transform(mom.begin(), mom.end(), momenta.begin(), [](const auto& part) {
-    return ROOT::Math::PxPyPzEVector{part.Px(), part.Py(), part.Pz(), part.energy()};
+  std::vector<std::pair<ROOT::Math::PxPyPzEVector, int>> momenta(parts.size());
+  // transform our simulation particle data into 4-momenta
+  std::transform(parts.begin(), parts.end(), momenta.begin(), [](const auto& part) {
+    double energy = sqrt(part.psx * part.psx + part.psy * part.psy + part.psz * part.psz +
+                         part.mass * part.mass);
+    return std::make_pair(ROOT::Math::PxPyPzEVector{part.psx, part.psy, part.psz, energy},
+                          part.pdgID);
   });
   return momenta;
 }
 
-// Momentum sorting bool
-bool sort_mom_bool(ROOT::Math::PxPyPzEVector &mom1, ROOT::Math::PxPyPzEVector &mom2)
-{
-  return  mom1.energy() > mom2.energy(); 
-}
-
-// Momentunm sorting function
-inline auto sort_momenta(const std::vector<ROOT::Math::PxPyPzEVector>& mom)
-{
-  std::vector <ROOT::Math::PxPyPzEVector> sort_mom = mom;
-  sort(sort_mom.begin(), sort_mom.end(), sort_mom_bool);
-  return sort_mom;
-}
-
-// Q2 calculation from 4 Vector
-inline auto Q2(const std::vector<ROOT::Math::PxPyPzEVector>& mom)
-{
-  std::vector<double> Q2Vec(mom.size() );
-  ROOT::Math::PxPyPzEVector beamMom = {0, 0, 18, 18};
-  std::transform(mom.begin(), mom.end(), Q2Vec.begin(), [beamMom](const auto& part) {
-    return -(part - beamMom).M2();
-  });
-  return Q2Vec;
-}
-
-// Difference between two 4 Vectors
-inline auto sub(const std::vector<ROOT::Math::PxPyPzEVector>& mom1, const std::vector<ROOT::Math::PxPyPzEVector>& mom2)
+// Finds particle with highest momentum within the event and returns the 4 Vector
+// Currently there is a cut so that only electrons are retained, but that can be changed
+// by chaning the PID
+inline auto find_scat(const std::vector<std::pair<ROOT::Math::PxPyPzEVector, int>>& mom)
 {
-  std::vector<ROOT::Math::PxPyPzEVector> sub(mom1.size() );
-  for (int i = 0; i < sub.size(); i++)
-  {
-    sub[i] = mom1[i] - mom2[i];
+  ROOT::Math::PxPyPzEVector elec_cand = {0, 0, 0, 0};
+  int                       partID    = 11; // For electrons
+  for (auto mom_cand : mom) {
+    if (mom_cand.first.E() > elec_cand.E() && mom_cand.second == partID) {
+      elec_cand = mom_cand.first;
+    }
   }
-  return sub;
-}
-
-// Multiplies a double by a gaussian distributed number 
-inline auto randomize(const std::vector<double>& doubleVec)
-{
-  std::vector<double> outVec(doubleVec.size() );
-  TRandom3 rand;
-  std::transform(doubleVec.begin(), doubleVec.end(), outVec.begin(), [&rand](const auto& indouble) {
-    return indouble * rand.Gaus(1.0, 0.2);
-  });
-  return outVec;
-}
-
-//Simulation functions
-/*
-bool mom_sort_sim(dd4pod::Geant4ParticleData& part1, dd4pod::Geant4ParticleData& part2)
-{
-  return (part1.psx*part1.psx + part1.psy*part1.psy + part1.psz*part1.psz > 
-          part2.psx*part2.psx + part2.psy*part2.psy + part2.psz*part2.psz);
-}
-
-inline auto momenta_from_sim(const std::vector<dd4pod::Geant4ParticleData>& parts)
-{
-  std::vector<dd4pod::Geant4ParticleData> sort_parts = parts;
-  sort(sort_parts.begin(), sort_parts.end(), mom_sort_sim);
-  return sort_parts;
+  return elec_cand;
 }
 
-inline auto Q2_from_sim(const std::vector<dd4pod::Geant4ParticleData>& parts)
+// Q2 calculation from 4 Vector
+// Recall, the electron ebeam is defined as coming from -z and the proton beam from +z
+inline auto Q2(ROOT::Math::PxPyPzEVector& mom)
 {
-  std::vector<double> Q2Vec(parts.size() );
-  double beamEnergy = 18;
-  std::transform(parts.begin(), parts.end(), Q2Vec.begin(), [beamEnergy](const auto& part) {
-    double energy = sqrt(part.psx*part.psx + part.psy*part.psy + part.psz*part.psz + part.mass*part.mass);
-    double q2 = pow(beamEnergy - energy, 2.0) - part.psx*part.psx - part.psy*part.psy - pow(part.psz - beamEnergy, 2.0);
-    return -q2;
-  });
-  return Q2Vec;
+  ROOT::Math::PxPyPzEVector beamMom = {0, 0, -eBeamEnergy, eBeamEnergy};
+  return -(mom - beamMom).M2();
 }
 
-inline auto elec_PID_sim(const std::vector<dd4pod::Geant4ParticleData>& parts)
+// Multiplies a double by a gaussian distributed number
+// The mean and sigma(STD) are set in teh code further down
+inline auto randomize(double& inDouble)
 {
-  std::vector<dd4pod::Geant4ParticleData> electrons;
-  for (auto part : parts)
-  {
-    if (part.pdgID == 11){electrons.push_back(part);}
-  }
-  return electrons;
+  TRandom3 rand(0);
+  return rand.Gaus(gausMean, gausSTD) * inDouble;
 }
-*/
-
 
 int dis_electrons(const std::string& config_name)
 {
@@ -140,6 +142,8 @@ int dis_electrons(const std::string& config_name)
   const std::string detector      = config["detector"];
   std::string       output_prefix = config["output_prefix"];
   const std::string test_tag      = config["test_tag"];
+  eBeamEnergy                     = config["ebeam"];
+  pBeamEnergy                     = config["pbeam"];
 
   fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green),
              "Running DIS electron analysis...\n");
@@ -159,22 +163,22 @@ int dis_electrons(const std::string& config_name)
        {"target", "0.1"}}};
 
   // Run this in multi-threaded mode if desired
-  //ROOT::EnableImplicitMT(kNumThreads);
-
-  //Particle number enumeration
-  enum sidis_particle_ID
-  {
-    electron = 11,
-    pi_0     = 111,
-    pi_plus  = 211,
-    pi_minus = -211,
-    k_0      = 311,
-    k_plus   = 321,
-    k_minus  = -321,
-    proton   = 2212,
-    neutron  = 2112
-  };
-
+  ROOT::EnableImplicitMT(kNumThreads);
+
+  // PIDs for reference
+  // electron = 11
+  // pi_0     = 111
+  // pi_plus  = 211
+  // pi_minus = -211
+  // k_0      = 311
+  // k_plus   = 321
+  // k_minus  = -321
+  // proton   = 2212
+  // neutron  = 2112
+
+  // Gausian variable declarations
+  gausMean = 1.0;
+  gausSTD  = 0.1;
 
   const double electron_mass = util::get_pdg_mass("electron");
 
@@ -193,75 +197,69 @@ int dis_electrons(const std::string& config_name)
         return util::momenta_from_tracking(tracks, electron_mass);
       };
 
-/*
-  //Old dataframe
-  auto d0 = d.Define("p_rec", momenta_from_tracking, {"outputTrackParameters"})
-                .Define("N", "p_rec.size()")
-                .Define("p_sim", util::momenta_from_simulation, {"mcparticles2"})
-                .Define("mom_sim", util::mom, {"p_sim"})
-                .Define("mom_rec", util::mom, {"p_rec"});
-*/
-
-  auto d0 = d.Define("p_recon", momenta_from_reconstruction, {"DummyReconstructedParticles"})
-                .Define("p_recon_sort", sort_momenta, {"p_recon"})
-                .Define("Q2_recon", Q2, {"p_recon_sort"})
-                .Define("Q2_recon_rand", randomize, {"Q2_recon"})
-                .Define("elec_Q2_recon_rand", "Q2_recon_rand[0]")
-                .Define("p_sim_M", util::momenta_from_simulation, {"mcparticles2"})
-                .Define("p_sim", convertMtoE, {"p_sim_M"})
-                .Define("p_sim_sort", sort_momenta, {"p_sim"})
-                .Define("Q2_sim", Q2, {"p_sim_sort"})
-                .Define("elec_Q2_sim", "Q2_sim[0]")
-                .Define("Q2_diff", "(elec_Q2_recon_rand - elec_Q2_sim)/elec_Q2_sim")
-                .Define("p_diff", sub, {"p_recon","p_sim"});
-                /*
-                .Define("electrons_sim", elec_PID_sim, {"sorted_sim"})
-                .Define("Q2_sim_elec_pid", Q2_from_sim, {"electrons_sim"})
-                .Define("elec_Q2_sim_pid", "Q2_sim_elec_pid[0]");
-                */
-  //Testing script
-/*
-  auto dis = d0.Display({"p_diff"});
-  //dis -> Print();
-  cout << *d0.Max<double>("elec_Q2_recon_rand") << " " << *d0.Min<double>("elec_Q2_recon_rand") << endl;
-  cout << *d0.Max<double>("elec_Q2_sim") << " " << *d0.Min<double>("elec_Q2_sim") << endl;
-/**/
-  //cout << *d0.Count() << endl;
-  //Momentum
-  //auto h_mom_sim = d0.Histo1D({"h_mom_sim", "; GeV; counts", 100, 0, 50}, "mom_sim");
-  //auto h_mom_rec = d0.Histo1D({"h_mom_rec", "; GeV; counts", 100, 0, 50}, "mom_rec");
-  //Q2
-  auto h_Q2_sim = d0.Histo1D({"h_Q2_sim", "; GeV; counts", 100, -5, 25}, "elec_Q2_sim");
-  auto h_Q2_rec = d0.Histo1D({"h_Q2_rec", "; GeV; counts", 100, -5, 25}, "elec_Q2_recon_rand");
-  auto h_Q2_res = d0.Histo1D({"h_Q2_res", ";    ; counts", 100, -10, 10}, "Q2_diff");
-/*
-  TH1D *h_Q2_res = (TH1D *)h_Q2_sim -> Clone();
-  TH1D *h_Q2_rec_copy = (TH1D *)h_Q2_rec -> Clone();
-  h_Q2_res -> Scale(1.0 / h_Q2_res -> Integral() );
-  h_Q2_res -> Add(h_Q2_rec_copy, -1.0 / h_Q2_rec_copy -> Integral() );
- */
-  TFitResultPtr f1  = h_Q2_res -> Fit("gaus", "S");
-  f1 -> Print("V");
-/*
-  printf("chisq %f A %f mean %f sigma %f \n", f1 -> Chi2(), 
-                                              f1 -> GetParameter(0), 
-                                              f1 -> GetParameter(1),
-                                              f1 -> GetParameter(2));
-*/
+  /*
+    //Old dataframe
+    auto d0 = d.Define("p_rec", momenta_from_tracking, {"outputTrackParameters"})
+                  .Define("N", "p_rec.size()")
+                  .Define("p_sim", util::momenta_from_simulation, {"mcparticles2"})
+                  .Define("mom_sim", util::mom, {"p_sim"})
+                  .Define("mom_rec", util::mom, {"p_rec"});
+  */
+
+  auto d0 = d.Define("p_recon", momenta_from_reconstruction_pid, {"DummyReconstructedParticles"})
+                .Define("elec_recon", find_scat, {"p_recon"})
+                .Define("elec_recon_Q2", Q2, {"elec_recon"})
+                .Define("elec_recon_Q2_rand", randomize, {"elec_recon_Q2"})
+                .Define("p_sim", momenta_from_simulation_pid, {"mcparticles2"})
+                .Define("elec_sim", find_scat, {"p_sim"})
+                .Define("elec_sim_Q2", Q2, {"elec_sim"})
+                .Define("dQ2", "elec_recon_Q2_rand - elec_sim_Q2")
+                .Define("dQ2_relative", "(elec_recon_Q2_rand - elec_sim_Q2)/elec_sim_Q2");
+  // Testing script
+  /*
+    //auto dis = d0.Display({"pt_recon", "elec_sim_Q2", "elec_recon_Q2"});
+    //dis -> Print();
+    std:vector<string> nameStr = {"elec_recon_Q2_rand", "elec_sim_Q2", "dQ2", "dQ2_relative"};
+    for (auto name : nameStr)
+    {
+      printf("%s Min(%f) Mean(%f) Max(%f)\n", name.c_str(), *d0.Min<double>(name),
+  *d0.Mean<double>(name), *d0.Max<double>(name) );
+    }
+  /**/
+  // Momentum
+  // auto h_mom_sim = d0.Histo1D({"h_mom_sim", "; GeV; counts", 100, 0, 50}, "mom_sim");
+  // auto h_mom_rec = d0.Histo1D({"h_mom_rec", "; GeV; counts", 100, 0, 50}, "mom_rec");
+
+  // Q2 Histograms
+  // auto h_Q2_sim = d0.Histo1D({"h_Q2_sim", "; GeV^{2}; Counts", 50, 0, 25}, "elec_sim_Q2");
+  // auto h_Q2_rec = d0.Histo1D({"h_Q2_rec", "; GeV^{2}; Counts", 50, 0, 25}, "elec_recon_Q2_rand");
+  auto h_Q2_res =
+      d0.Histo1D({"h_Q2_res", "; Q^{2} Resolution (GeV^{2}); Counts", 50, -120, 120}, "dQ2");
+  auto h_Q2_rel_res = d0.Histo1D(
+      {"h_Q2_rel_res", "; Q^{2} Relative Resolution ; Counts", 50, -2, 2}, "dQ2_relative");
+
+  TFitResultPtr f1  = h_Q2_rel_res->Fit("gaus", "S");
+  const double* res = f1->GetParams();
+  if (res[2] <= 0.1) {
+    dis_Q2_resolution.pass(res[2]);
+  } else {
+    dis_Q2_resolution.fail(res[2]);
+  }
+  f1->Print("V");
 
-  
   auto c = new TCanvas();
 
   // Plot our histograms.
   // TODO: to start I'm explicitly plotting the histograms, but want to
   // factorize out the plotting code moving forward.
+  /*
   {
     TCanvas c{"canvas", "canvas", 1200, 1200};
     c.cd();
     // gPad->SetLogx(false);
     gPad->SetLogy(true);
-    //auto& h1 = *h_mom_sim;
-    //auto& h2 = *h_mom_rec;
+    // auto& h1 = *h_mom_sim;
+    // auto& h2 = *h_mom_rec;
     auto& h1 = *h_Q2_sim;
     auto& h2 = *h_Q2_rec;
     // histogram style
@@ -276,7 +274,7 @@ int dis_electrons(const std::string& config_name)
     h1.DrawClone("hist");
     h2.DrawClone("hist same");
     // FIXME hardcoded beam configuration
-    plot::draw_label(18, 275, detector);
+    plot::draw_label(eBeamEnergy, pBeamEnergy, detector);
     TText* tptr1;
     auto   t1 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
     t1->SetFillColorAlpha(kWhite, 0);
@@ -290,6 +288,9 @@ int dis_electrons(const std::string& config_name)
 
     c.Print(fmt::format("{}momentum.png", output_prefix).c_str());
   }
+  */
+  plot_resolutions(h_Q2_res, h_Q2_rel_res, "Q2", detector, output_prefix);
+
   eic::util::write_test({dis_Q2_resolution}, fmt::format("{}dis_electrons.json", output_prefix));
 
   return 0;