Skip to content
Snippets Groups Projects
Commit 3d618a33 authored by Marshall Scott's avatar Marshall Scott Committed by Sylvester Joosten
Browse files

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

parent 776434d9
Branches
No related tags found
1 merge request!37Ongoing dis_electron Q2 analysis. An issue with Q2 being basically zero for the electrons.
...@@ -6,92 +6,48 @@ ...@@ -6,92 +6,48 @@
#include <util.h> #include <util.h>
#include "ROOT/RDataFrame.hxx" #include "ROOT/RDataFrame.hxx"
#include <TFitResult.h>
#include <TH1D.h>
#include <TRandom3.h>
#include <algorithm>
#include <cmath> #include <cmath>
#include <eicd/ReconstructedParticleData.h>
#include <fmt/color.h> #include <fmt/color.h>
#include <fmt/core.h> #include <fmt/core.h>
#include <fstream> #include <fstream>
#include <iostream> #include <iostream>
#include <nlohmann/json.hpp> #include <nlohmann/json.hpp>
#include <string> #include <string>
#include <vector>
#include <eicd/ReconstructedParticleData.h>
#include <TH1D.h>
#include <TFitResult.h>
#include <TRandom3.h>
#include <algorithm>
#include <utility> #include <utility>
#include <vector>
// Reconstuction functions
// Get a vector of 4-momenta from the reconstructed data. bool mom_sort_recon(eic::ReconstructedParticleData& part1, eic::ReconstructedParticleData& part2)
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)
{ {
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 momenta_from_recon(const std::vector<eic::ReconstructedParticleData>& parts)
inline auto sort_momenta(const std::vector<ROOT::Math::PxPyPzEVector>& mom)
{ {
std::vector <ROOT::Math::PxPyPzEVector> sort_mom = mom; std::vector<eic::ReconstructedParticleData> sort_parts = parts;
sort(sort_mom.begin(), sort_mom.end(), sort_mom_bool); sort(sort_parts.begin(), sort_parts.end(), mom_sort_recon);
return sort_mom; return sort_parts;
} }
// Q2 calculation from 4 Vector inline auto Q2_from_recon(const std::vector<eic::ReconstructedParticleData>& parts)
inline auto Q2(const std::vector<ROOT::Math::PxPyPzEVector>& mom)
{ {
std::vector<double> Q2Vec(mom.size() ); std::vector<double> Q2Vec(parts.size());
ROOT::Math::PxPyPzEVector beamMom = {0, 0, 18, 18}; double beamEnergy = 18;
std::transform(mom.begin(), mom.end(), Q2Vec.begin(), [beamMom](const auto& part) { std::transform(parts.begin(), parts.end(), Q2Vec.begin(), [beamEnergy](const auto& part) {
return -(part - beamMom).M2(); 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; 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) bool mom_sort_sim(dd4pod::Geant4ParticleData& part1, dd4pod::Geant4ParticleData& part2)
{ {
return (part1.psx * part1.psx + part1.psy * part1.psy + part1.psz * part1.psz > return (part1.psx * part1.psx + part1.psy * part1.psy + part1.psz * part1.psz >
...@@ -110,8 +66,10 @@ inline auto Q2_from_sim(const std::vector<dd4pod::Geant4ParticleData>& parts) ...@@ -110,8 +66,10 @@ inline auto Q2_from_sim(const std::vector<dd4pod::Geant4ParticleData>& parts)
std::vector<double> Q2Vec(parts.size()); std::vector<double> Q2Vec(parts.size());
double beamEnergy = 18; double beamEnergy = 18;
std::transform(parts.begin(), parts.end(), Q2Vec.begin(), [beamEnergy](const auto& part) { 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 energy = sqrt(part.psx * part.psx + part.psy * part.psy + part.psz * part.psz +
double q2 = pow(beamEnergy - energy, 2.0) - part.psx*part.psx - part.psy*part.psy - pow(part.psz - beamEnergy, 2.0); 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 -q2;
}); });
return Q2Vec; return Q2Vec;
...@@ -120,14 +78,46 @@ inline auto Q2_from_sim(const std::vector<dd4pod::Geant4ParticleData>& parts) ...@@ -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) inline auto elec_PID_sim(const std::vector<dd4pod::Geant4ParticleData>& parts)
{ {
std::vector<dd4pod::Geant4ParticleData> electrons; std::vector<dd4pod::Geant4ParticleData> electrons;
for (auto part : parts) for (auto part : parts) {
{ if (part.pdgID == 11) {
if (part.pdgID == 11){electrons.push_back(part);} electrons.push_back(part);
}
} }
return electrons; 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) int dis_electrons(const std::string& config_name)
{ {
...@@ -162,8 +152,7 @@ int dis_electrons(const std::string& config_name) ...@@ -162,8 +152,7 @@ int dis_electrons(const std::string& config_name)
// ROOT::EnableImplicitMT(kNumThreads); // ROOT::EnableImplicitMT(kNumThreads);
// Particle number enumeration // Particle number enumeration
enum sidis_particle_ID enum sidis_particle_ID {
{
electron = 11, electron = 11,
pi_0 = 111, pi_0 = 111,
pi_plus = 211, pi_plus = 211,
...@@ -175,7 +164,6 @@ int dis_electrons(const std::string& config_name) ...@@ -175,7 +164,6 @@ int dis_electrons(const std::string& config_name)
neutron = 2112 neutron = 2112
}; };
const double electron_mass = util::get_pdg_mass("electron"); const double electron_mass = util::get_pdg_mass("electron");
// Ensure our output prefix always ends on a dot, a slash or a dash // Ensure our output prefix always ends on a dot, a slash or a dash
...@@ -202,44 +190,34 @@ int dis_electrons(const std::string& config_name) ...@@ -202,44 +190,34 @@ int dis_electrons(const std::string& config_name)
.Define("mom_rec", util::mom, {"p_rec"}); .Define("mom_rec", util::mom, {"p_rec"});
*/ */
auto d0 = d.Define("p_recon", momenta_from_reconstruction, {"DummyReconstructedParticles"}) auto d0 = d.Define("sorted_recon", momenta_from_recon, {"DummyReconstructedParticles"})
.Define("p_recon_sort", sort_momenta, {"p_recon"}) .Define("Q2_recon", Q2_from_recon, {"sorted_recon"})
.Define("Q2_recon", Q2, {"p_recon_sort"}) .Define("elec_Q2_recon", "Q2_recon[0]")
.Define("Q2_recon_rand", randomize, {"Q2_recon"}) .Define("sorted_sim", momenta_from_sim, {"mcparticles2"})
.Define("elec_Q2_recon_rand", "Q2_recon_rand[0]") .Define("Q2_sim", Q2_from_sim, {"sorted_sim"})
.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("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("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]"); .Define("elec_Q2_sim_pid", "Q2_sim_elec_pid[0]");
*/
// Testing script // Testing script
/* /*
auto dis = d0.Display({"p_diff"}); auto dis = d0.Display({"Q2_sim", "elec_Q2_recon", "elec_Q2_sim", "elec_Q2_sim_pid"});
//dis -> Print(); dis -> Print();
cout << *d0.Max<double>("elec_Q2_recon_rand") << " " << *d0.Min<double>("elec_Q2_recon_rand") << endl; cout << *d0.Max<double>("elec_Q2_recon") << " " << *d0.Min<double>("elec_Q2_recon") << endl;
cout << *d0.Max<double>("elec_Q2_sim") << " " << *d0.Min<double>("elec_Q2_sim") << endl; */
/**/
// cout << *d0.Count() << endl; // cout << *d0.Count() << endl;
// Momentum // Momentum
// auto h_mom_sim = d0.Histo1D({"h_mom_sim", "; GeV; counts", 100, 0, 50}, "mom_sim"); // 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_rec = d0.Histo1D({"h_mom_rec", "; GeV; counts", 100, 0, 50}, "mom_rec");
// Q2 // Q2
auto h_Q2_sim = d0.Histo1D({"h_Q2_sim", "; GeV; counts", 100, -5, 25}, "elec_Q2_sim"); 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, -5, 25}, "elec_Q2_recon_rand"); auto h_Q2_rec = d0.Histo1D({"h_Q2_rec", "; GeV; counts", 100, -50, 50}, "elec_Q2_recon");
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_res = (TH1D*)h_Q2_sim->Clone();
TH1D* h_Q2_rec_copy = (TH1D*)h_Q2_rec->Clone(); TH1D* h_Q2_rec_copy = (TH1D*)h_Q2_rec->Clone();
h_Q2_res->Scale(1.0 / h_Q2_res->Integral()); 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()); h_Q2_res->Add(h_Q2_rec_copy, -1.0 / h_Q2_rec_copy->Integral());
*/
TFitResultPtr f1 = h_Q2_res->Fit("gaus", "S"); TFitResultPtr f1 = h_Q2_res->Fit("gaus", "S");
f1->Print("V"); f1->Print("V");
/* /*
...@@ -249,7 +227,6 @@ int dis_electrons(const std::string& config_name) ...@@ -249,7 +227,6 @@ int dis_electrons(const std::string& config_name)
f1 -> GetParameter(2)); f1 -> GetParameter(2));
*/ */
auto c = new TCanvas(); auto c = new TCanvas();
// Plot our histograms. // Plot our histograms.
...@@ -264,6 +241,7 @@ int dis_electrons(const std::string& config_name) ...@@ -264,6 +241,7 @@ int dis_electrons(const std::string& config_name)
// auto& h2 = *h_mom_rec; // auto& h2 = *h_mom_rec;
auto& h1 = *h_Q2_sim; auto& h1 = *h_Q2_sim;
auto& h2 = *h_Q2_rec; auto& h2 = *h_Q2_rec;
h2.Scale(2.0);
// histogram style // histogram style
h1.SetLineColor(plot::kMpBlue); h1.SetLineColor(plot::kMpBlue);
h1.SetLineWidth(2); h1.SetLineWidth(2);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment