Re-enable DIS

benchmarks/dis/analysis/dis_electrons.cxx

 #include "common_bench/benchmark.h"
 #include "common_bench/mt.h"
 #include "common_bench/util.h"
+#include "common_bench/plot.h"
 
 #include <cmath>
 #include <fstream>
@@ -14,6 +15,7 @@
 #include <TH1D.h>
 #include <TFitResult.h>
 #include <TRandom3.h>
+#include <TCanvas.h>
 
 #include "fmt/color.h"
 #include "fmt/core.h"
@@ -21,6 +23,16 @@
 #include "nlohmann/json.hpp"
 #include "eicd/ReconstructedParticleData.h"
 
+// Select only electrons
+inline auto select_electrons(const std::vector<eic::ReconstructedParticleData>& parts)
+{
+  std::vector<eic::ReconstructedParticleData> electrons;
+  for (auto &part: parts)
+    if (part.pid == 11)
+      electrons.push_back(part);
+  return electrons;
+}
+
 // Get a vector of 4-momenta from the reconstructed data.
 inline auto momenta_from_reconstruction(const std::vector<eic::ReconstructedParticleData>& parts)
 {
@@ -42,6 +54,15 @@ inline auto convertMtoE(const std::vector<ROOT::Math::PxPyPzMVector>& mom)
   return momenta;
 }
 
+// Momentum
+auto momentum(std::vector<ROOT::Math::PxPyPzEVector> const& in) {
+  std::vector<double> result;
+  for (size_t i = 0; i < in.size(); ++i) {
+   result.push_back(in[i].P());
+  }
+  return result;
+};
+
 // Momentum sorting bool
 bool sort_mom_bool(ROOT::Math::PxPyPzEVector &mom1, ROOT::Math::PxPyPzEVector &mom2)
 {
@@ -60,7 +81,7 @@ inline auto sort_momenta(const std::vector<ROOT::Math::PxPyPzEVector>& mom)
 inline auto Q2(const std::vector<ROOT::Math::PxPyPzEVector>& mom)
 {
   std::vector<double> Q2Vec(mom.size() );
-  ROOT::Math::PxPyPzEVector beamMom = {0, 0, 18, 18};
+  ROOT::Math::PxPyPzEVector beamMom = {0, 0, -18, 18};
   std::transform(mom.begin(), mom.end(), Q2Vec.begin(), [beamMom](const auto& part) {
     return -(part - beamMom).M2();
   });
@@ -71,7 +92,7 @@ inline auto Q2(const std::vector<ROOT::Math::PxPyPzEVector>& mom)
 inline auto sub(const std::vector<ROOT::Math::PxPyPzEVector>& mom1, const std::vector<ROOT::Math::PxPyPzEVector>& mom2)
 {
   std::vector<ROOT::Math::PxPyPzEVector> sub(mom1.size() );
-  for (int i = 0; i < sub.size(); i++)
+  for (size_t i = 0; i < sub.size(); i++)
   {
     sub[i] = mom1[i] - mom2[i];
   }
@@ -158,7 +179,7 @@ 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
@@ -175,7 +196,7 @@ int dis_electrons(const std::string& config_name)
   };
 
 
-  const double electron_mass = util::get_pdg_mass("electron");
+  const double electron_mass = common_bench::get_pdg_mass("electron");
 
   // Ensure our output prefix always ends on a dot, a slash or a dash
   // Necessary when generating output plots
@@ -189,84 +210,94 @@ int dis_electrons(const std::string& config_name)
   // (as we have no PID yet)
   auto momenta_from_tracking =
       [electron_mass](const std::vector<eic::TrackParametersData>& tracks) {
-        return util::momenta_from_tracking(tracks, electron_mass);
+        return common_bench::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, {"mcparticles"})
-                .Define("mom_sim", util::mom, {"p_sim"})
-                .Define("mom_rec", util::mom, {"p_rec"});
-*/
+  auto d0 = d   .Define("electrons_rec", select_electrons, {"ReconstructedParticles"})
+                .Define("p_rec", momenta_from_reconstruction, {"electrons_rec"})
+                .Define("p_rec_sort", sort_momenta, {"p_rec"})
+                .Define("Q2_rec", Q2, {"p_rec_sort"})
+                .Define("Q2_rec_rand", randomize, {"Q2_rec"})
+                .Define("elec_Q2_rec_rand", "Q2_rec_rand[0]")
 
-  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, {"mcparticles"})
-                .Define("p_sim", convertMtoE, {"p_sim_M"})
+                .Define("electrons_sim", select_electrons, {"GeneratedParticles"})
+                .Define("p_sim", momenta_from_reconstruction, {"electrons_sim"})
                 .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("Q2_diff", "(elec_Q2_rec_rand - elec_Q2_sim)/elec_Q2_sim")
+                .Define("p_diff", sub, {"p_rec_sort","p_sim_sort"})
+                .Define("mom_rec", momentum, {"p_rec_sort"})
+                .Define("mom_sim", momentum, {"p_sim_sort"})
+                .Define("mom_diff", "(mom_rec[0]-mom_sim[0])/mom_sim[0]")
                 /*
                 .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]");
+                .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");
+  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");
+  auto h_mom_res = d0.Histo1D({"h_mom_res", ";    ; counts", 100, -10, 10}, "mom_diff");
   //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_rec = d0.Histo1D({"h_Q2_rec", "; GeV; counts", 100, -5, 25}, "elec_Q2_rec_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));
-*/
 
-  
-  auto c = new TCanvas();
+  //TFitResultPtr f1 = h_Q2_res -> Fit("gaus", "S");
+  //f1 -> Print("V");
 
   // Plot our histograms.
   // TODO: to start I'm explicitly plotting the histograms, but want to
   // factorize out the plotting code moving forward.
+
+  // Momentum
   {
-    TCanvas c{"canvas", "canvas", 1200, 1200};
+    TCanvas c{"c", "c", 1200, 1200};
     c.cd();
-    // gPad->SetLogx(false);
-    gPad->SetLogy(true);
-    //auto& h1 = *h_mom_sim;
-    //auto& h2 = *h_mom_rec;
+    //gPad->SetLogy(true);
+    auto& h1 = *h_mom_sim;
+    auto& h2 = *h_mom_rec;
+    // histogram style
+    h1.SetLineColor(common_bench::plot::kMpBlue);
+    h1.SetLineWidth(2);
+    h2.SetLineColor(common_bench::plot::kMpOrange);
+    h2.SetLineWidth(2);
+    // axes
+    h1.GetXaxis()->CenterTitle();
+    h1.GetYaxis()->CenterTitle();
+    // draw everything
+    h1.DrawClone("hist");
+    h2.DrawClone("hist same");
+    // FIXME hardcoded beam configuration
+    common_bench::plot::draw_label(18, 275, detector);
+    TText* tptr1;
+    TPaveText t1(.6, .8417, .9, .925, "NB NDC");
+    t1.SetFillColorAlpha(kWhite, 0);
+    t1.SetTextFont(43);
+    t1.SetTextSize(25);
+    tptr1 = t1.AddText("simulated");
+    tptr1->SetTextColor(common_bench::plot::kMpBlue);
+    tptr1 = t1.AddText("reconstructed");
+    tptr1->SetTextColor(common_bench::plot::kMpOrange);
+    t1.Draw();
+    c.Print(fmt::format("{}mom.png", output_prefix).c_str());
+  }
+
+  // Q2
+  {
+    TCanvas c{"c", "c", 1200, 1200};
+    c.cd();
+    //gPad->SetLogy(true);
     auto& h1 = *h_Q2_sim;
     auto& h2 = *h_Q2_rec;
     // histogram style
-    h1.SetLineColor(plot::kMpBlue);
+    h1.SetLineColor(common_bench::plot::kMpBlue);
     h1.SetLineWidth(2);
-    h2.SetLineColor(plot::kMpOrange);
+    h2.SetLineColor(common_bench::plot::kMpOrange);
     h2.SetLineWidth(2);
     // axes
     h1.GetXaxis()->CenterTitle();
@@ -275,20 +306,20 @@ 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);
+    common_bench::plot::draw_label(18, 275, detector);
     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("simulated");
-    tptr1->SetTextColor(plot::kMpBlue);
-    tptr1 = t1->AddText("reconstructed");
-    tptr1->SetTextColor(plot::kMpOrange);
-    t1->Draw();
-
-    c.Print(fmt::format("{}momentum.png", output_prefix).c_str());
+    TPaveText t1(.6, .8417, .9, .925, "NB NDC");
+    t1.SetFillColorAlpha(kWhite, 0);
+    t1.SetTextFont(43);
+    t1.SetTextSize(25);
+    tptr1 = t1.AddText("simulated");
+    tptr1->SetTextColor(common_bench::plot::kMpBlue);
+    tptr1 = t1.AddText("reconstructed");
+    tptr1->SetTextColor(common_bench::plot::kMpOrange);
+    t1.Draw();
+    c.Print(fmt::format("{}Q2.png", output_prefix).c_str());
   }
+
   common_bench::write_test({dis_Q2_resolution}, fmt::format("{}dis_electrons.json", output_prefix));
 
   return 0;