From 3d618a33fecc88fb8490b3f857179317ec4ab873 Mon Sep 17 00:00:00 2001
From: Marshall Scott <mbscott@anl.gov>
Date: Mon, 1 Mar 2021 10:09:30 -0500
Subject: [PATCH] Ongoing dis_electron Q2 analysis. An issue with Q2 being
basically zero for the electrons.
---
benchmarks/dis/analysis/dis_electrons.cxx | 250 ++++++++++------------
1 file changed, 114 insertions(+), 136 deletions(-)
diff --git a/benchmarks/dis/analysis/dis_electrons.cxx b/benchmarks/dis/analysis/dis_electrons.cxx
index f0a19595..63515079 100644
--- a/benchmarks/dis/analysis/dis_electrons.cxx
+++ b/benchmarks/dis/analysis/dis_electrons.cxx
@@ -6,96 +6,52 @@
#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>
#include <iostream>
#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>
+#include <vector>
-
-// Get a vector of 4-momenta from the reconstructed data.
-inline auto momenta_from_reconstruction(const std::vector<eic::ReconstructedParticleData>& parts)
-{
- std::vector<ROOT::Math::PxPyPzEVector> 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 momenta;
-}
-
-// Convert PxPyPzMVector to PxPyPzEVector
-inline auto convertMtoE(const std::vector<ROOT::Math::PxPyPzMVector>& mom)
-{
- 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()};
- });
- return momenta;
-}
-
-// Momentum sorting bool
-bool sort_mom_bool(ROOT::Math::PxPyPzEVector &mom1, ROOT::Math::PxPyPzEVector &mom2)
+// Reconstuction functions
+bool mom_sort_recon(eic::ReconstructedParticleData& part1, eic::ReconstructedParticleData& part2)
{
- return mom1.energy() > mom2.energy();
+ 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);
}
-// Momentunm sorting function
-inline auto sort_momenta(const std::vector<ROOT::Math::PxPyPzEVector>& mom)
+inline auto momenta_from_recon(const std::vector<eic::ReconstructedParticleData>& parts)
{
- std::vector <ROOT::Math::PxPyPzEVector> sort_mom = mom;
- sort(sort_mom.begin(), sort_mom.end(), sort_mom_bool);
- return sort_mom;
+ std::vector<eic::ReconstructedParticleData> sort_parts = parts;
+ sort(sort_parts.begin(), sort_parts.end(), mom_sort_recon);
+ return sort_parts;
}
-// Q2 calculation from 4 Vector
-inline auto Q2(const std::vector<ROOT::Math::PxPyPzEVector>& mom)
+inline auto Q2_from_recon(const std::vector<eic::ReconstructedParticleData>& parts)
{
- 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();
+ 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;
}
-// Difference between two 4 Vectors
-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++)
- {
- sub[i] = mom1[i] - mom2[i];
- }
- 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
-/*
+// 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);
+ 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)
@@ -107,11 +63,13 @@ inline auto momenta_from_sim(const std::vector<dd4pod::Geant4ParticleData>& part
inline auto Q2_from_sim(const std::vector<dd4pod::Geant4ParticleData>& parts)
{
- std::vector<double> Q2Vec(parts.size() );
- double beamEnergy = 18;
+ 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);
+ 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;
@@ -120,14 +78,46 @@ inline auto Q2_from_sim(const std::vector<dd4pod::Geant4ParticleData>& parts)
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);}
+ 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)
{
@@ -159,11 +149,10 @@ 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
- {
+ // Particle number enumeration
+ enum sidis_particle_ID {
electron = 11,
pi_0 = 111,
pi_plus = 211,
@@ -175,7 +164,6 @@ int dis_electrons(const std::string& config_name)
neutron = 2112
};
-
const double electron_mass = util::get_pdg_mass("electron");
// Ensure our output prefix always ends on a dot, a slash or a dash
@@ -193,63 +181,52 @@ 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"})
+ /*
+ //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("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("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({"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));
-*/
+ // 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");
+
+ 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.
@@ -260,10 +237,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);
--
GitLab