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

@@ -14,6 +14,112 @@
 #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>
+
+
+//Reconstuction functions
+bool mom_sort_recon(eic::ReconstructedParticleData &part1, eic::ReconstructedParticleData &part2)
+{
+  return (part1.p.x*part1.p.x + part1.p.y*part1.p.y + part1.p.z*part1.p.z > 
+          part2.p.x*part2.p.x + part2.p.y*part2.p.y + part2.p.z*part2.p.z);
+}
+
+inline auto momenta_from_recon(const std::vector<eic::ReconstructedParticleData>& parts)
+{
+  std::vector <eic::ReconstructedParticleData> sort_parts = parts;
+  sort(sort_parts.begin(), sort_parts.end(), mom_sort_recon);
+  return sort_parts;
+}
+
+inline auto Q2_from_recon(const std::vector<eic::ReconstructedParticleData>& parts)
+{
+  std::vector<double> Q2Vec(parts.size() );
+  double beamEnergy = 18;
+  std::transform(parts.begin(), parts.end(), Q2Vec.begin(), [beamEnergy](const auto& part) {
+    double q2 = pow(beamEnergy - part.energy, 2.0) - part.p.x*part.p.x - part.p.y*part.p.y - pow(beamEnergy - part.p.z, 2.0);
+    return -q2;
+  });
+  return Q2Vec;
+}
+
+//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;
+}
+
+inline auto Q2_from_sim(const std::vector<dd4pod::Geant4ParticleData>& parts)
+{
+  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;
+}
+
+inline auto elec_PID_sim(const std::vector<dd4pod::Geant4ParticleData>& parts)
+{
+  std::vector<dd4pod::Geant4ParticleData> electrons;
+  for (auto part : parts)
+  {
+    if (part.pdgID == 11){electrons.push_back(part);}
+  }
+  return electrons;
+}
+
+inline auto mass_from_recon(const std::vector<eic::ReconstructedParticleData>& parts)
+{
+  std::vector<double> mass(parts.size() );
+  std::transform(parts.begin(), parts.end(), mass.begin(), [](const auto& part) {
+    return part.mass;
+  });
+  return mass;
+}
+
+inline auto mass_random(const std::vector<double>& massVec)
+{
+  std::vector<double> mass(massVec.size() );
+  TRandom3 rand;
+  std::transform(massVec.begin(), massVec.end(), mass.begin(), [&rand](const auto& inmass) {
+    return inmass * rand.Gaus(1.0, 0.2);
+  });
+  return mass;
+}
+
+inline auto mass_from_simulation(const std::vector<ROOT::Math::PxPyPzMVector>& momenta)
+{
+  std::vector<double> mass(momenta.size() );
+  std::transform(momenta.begin(), momenta.end(), mass.begin(), [](const auto& mom) {
+    return mom.mass();
+  });
+  return mass;
+}
+
+inline auto getSimPID(const std::vector<dd4pod::Geant4ParticleData>& parts)
+{
+  std::vector<int> pid(parts.size() );
+  std::transform(parts.begin(), parts.end(), pid.begin(), [](const auto& part) {
+    return part.pdgID;
+  });
+  return pid;
+}
+
 
 int dis_electrons(const std::string& config_name)
 {
@@ -45,7 +151,22 @@ int dis_electrons(const std::string& config_name)
        {"target", "0.1"}}};
 
   // Run this in multi-threaded mode if desired
-  ROOT::EnableImplicitMT(kNumThreads);
+  //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
+  };
+
 
   const double electron_mass = util::get_pdg_mass("electron");
 
@@ -64,15 +185,53 @@ 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("sorted_recon", momenta_from_recon, {"DummyReconstructedParticles"})
+                .Define("Q2_recon", Q2_from_recon, {"sorted_recon"})
+                .Define("elec_Q2_recon", "Q2_recon[0]")
+                .Define("sorted_sim", momenta_from_sim, {"mcparticles2"})
+                .Define("Q2_sim", Q2_from_sim, {"sorted_sim"})
+                .Define("elec_Q2_sim", "Q2_sim[0]")
+                .Define("electrons_sim", elec_PID_sim, {"sorted_sim"})
+                .Define("Q2_sim_elec_pid", Q2_from_sim, {"sorted_sim"})
+                .Define("elec_Q2_sim_pid", "Q2_sim_elec_pid[0]");
+  //Testing script
+/*
+  auto dis = d0.Display({"Q2_sim", "elec_Q2_recon", "elec_Q2_sim", "elec_Q2_sim_pid"});
+  dis -> Print();
+  cout << *d0.Max<double>("elec_Q2_recon") << " " << *d0.Min<double>("elec_Q2_recon") << 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, -50, 50}, "elec_Q2_sim");
+  auto h_Q2_rec = d0.Histo1D({"h_Q2_rec", "; GeV; counts", 100, -50, 50}, "elec_Q2_recon");
 
-  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");
+  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));
+*/
 
+  
   auto c = new TCanvas();
 
   // Plot our histograms.
@@ -83,8 +242,11 @@ int dis_electrons(const std::string& config_name)
     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;
+    h2.Scale(2.0);
     // histogram style
     h1.SetLineColor(plot::kMpBlue);
     h1.SetLineWidth(2);