Skip to content
Snippets Groups Projects
Commit f1cc8329 authored by Wouter Deconinck's avatar Wouter Deconinck Committed by Wouter Deconinck
Browse files

Select electrons only; use GeneratedParticles

parent e9a7de19
No related branches found
No related tags found
1 merge request!81Re-enable DIS
This commit is part of merge request !81. Comments created here will be created in the context of that merge request.
...@@ -23,6 +23,16 @@ ...@@ -23,6 +23,16 @@
#include "nlohmann/json.hpp" #include "nlohmann/json.hpp"
#include "eicd/ReconstructedParticleData.h" #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 (parts.pid == 11)
electrons.push_back(parts);
return electrons;
}
// Get a vector of 4-momenta from the reconstructed data. // 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(const std::vector<eic::ReconstructedParticleData>& parts)
{ {
...@@ -62,7 +72,7 @@ inline auto sort_momenta(const std::vector<ROOT::Math::PxPyPzEVector>& mom) ...@@ -62,7 +72,7 @@ inline auto sort_momenta(const std::vector<ROOT::Math::PxPyPzEVector>& mom)
inline auto Q2(const std::vector<ROOT::Math::PxPyPzEVector>& mom) inline auto Q2(const std::vector<ROOT::Math::PxPyPzEVector>& mom)
{ {
std::vector<double> Q2Vec(mom.size() ); 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) { std::transform(mom.begin(), mom.end(), Q2Vec.begin(), [beamMom](const auto& part) {
return -(part - beamMom).M2(); return -(part - beamMom).M2();
}); });
...@@ -203,23 +213,29 @@ int dis_electrons(const std::string& config_name) ...@@ -203,23 +213,29 @@ int dis_electrons(const std::string& config_name)
.Define("mom_rec", common_bench::mom, {"p_rec"}); .Define("mom_rec", common_bench::mom, {"p_rec"});
*/ */
auto d0 = d.Define("p_recon", momenta_from_reconstruction, {"ReconstructedParticles"}) auto d0 = d .Define("electrons_rec", select_electrons, {"ReconstructedParticles"})
.Define("p_recon_sort", sort_momenta, {"p_recon"}) .Define("p_rec", momenta_from_reconstruction, {"electrons_rec"})
.Define("Q2_recon", Q2, {"p_recon_sort"}) .Define("p_rec_sort", sort_momenta, {"p_rec"})
.Define("Q2_recon_rand", randomize, {"Q2_recon"}) .Define("Q2_rec", Q2, {"p_rec_sort"})
.Define("elec_Q2_recon_rand", "Q2_recon_rand[0]") .Define("Q2_rec_rand", randomize, {"Q2_rec"})
.Define("p_sim_M", common_bench::momenta_from_simulation, {"mcparticles2"}) .Define("elec_Q2_rec_rand", "Q2_rec_rand[0]")
.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("p_sim_sort", sort_momenta, {"p_sim"})
.Define("Q2_sim", Q2, {"p_sim_sort"}) .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("Q2_diff", "(elec_Q2_rec_rand - elec_Q2_sim)/elec_Q2_sim")
.Define("p_diff", sub, {"p_rec","p_sim"})
.Define("mom_rec", common_bench::mom, {"p_rec"})
.Define("mom_sim", common_bench::mom, {"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, {"electrons_sim"})
.Define("elec_Q2_sim_pid", "Q2_sim_elec_pid[0]"); .Define("elec_Q2_sim_pid", "Q2_sim_elec_pid[0]");
*/ */
;
//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");
...@@ -227,7 +243,7 @@ int dis_electrons(const std::string& config_name) ...@@ -227,7 +243,7 @@ int dis_electrons(const std::string& config_name)
auto h_mom_res = d0.Histo1D({"h_mom_res", "; ; counts", 100, -10, 10}, "p_diff"); auto h_mom_res = d0.Histo1D({"h_mom_res", "; ; counts", 100, -10, 10}, "p_diff");
//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, -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"); auto h_Q2_res = d0.Histo1D({"h_Q2_res", "; ; counts", 100, -10, 10}, "Q2_diff");
/* /*
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment