diff --git a/benchmarks/dis/analysis/dis_electrons.cxx b/benchmarks/dis/analysis/dis_electrons.cxx
index c3ce781f2e804e36593ad65fdd1d5c27ba4c19c8..6a4d9896210442546ce282fabf492db64f7792f6 100644
--- a/benchmarks/dis/analysis/dis_electrons.cxx
+++ b/benchmarks/dis/analysis/dis_electrons.cxx
@@ -1,6 +1,7 @@
#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,120 +15,15 @@
#include <TH1D.h>
#include <TFitResult.h>
#include <TRandom3.h>
+#include <TCanvas.h>
#include "fmt/color.h"
#include "fmt/core.h"
#include "nlohmann/json.hpp"
+#include "eicd/InclusiveKinematicsData.h"
#include "eicd/ReconstructedParticleData.h"
-// 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)
-{
- return mom1.energy() > mom2.energy();
-}
-
-// Momentunm sorting function
-inline auto sort_momenta(const std::vector<ROOT::Math::PxPyPzEVector>& mom)
-{
- std::vector <ROOT::Math::PxPyPzEVector> sort_mom = mom;
- sort(sort_mom.begin(), sort_mom.end(), sort_mom_bool);
- return sort_mom;
-}
-
-// Q2 calculation from 4 Vector
-inline auto Q2(const std::vector<ROOT::Math::PxPyPzEVector>& mom)
-{
- 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();
- });
- 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
-/*
-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;
-}
-*/
-
-
int dis_electrons(const std::string& config_name)
{
// read our configuration
@@ -137,7 +33,7 @@ int dis_electrons(const std::string& config_name)
const std::string rec_file = config["rec_file"];
const std::string detector = config["detector"];
- std::string output_prefix = config["output_prefix"];
+ const std::string output_prefix = config["output_prefix"];
const std::string test_tag = config["test_tag"];
fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green),
@@ -158,115 +54,62 @@ int dis_electrons(const std::string& config_name)
{"target", "0.1"}}};
// Run this in multi-threaded mode if desired
- //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");
-
- // Ensure our output prefix always ends on a dot, a slash or a dash
- // Necessary when generating output plots
- if (output_prefix.back() != '.' && output_prefix.back() != '/' && output_prefix.back() != '-') {
- output_prefix += "-";
- }
-
+ ROOT::EnableImplicitMT(kNumThreads);
ROOT::RDataFrame d("events", rec_file);
- // utility lambda functions to bind the reconstructed particle type
- // (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);
- };
+ auto combinatorial_diff_ratio = [] (
+ const ROOT::VecOps::RVec<float>& v1,
+ const ROOT::VecOps::RVec<float>& v2
+ ) {
+ std::vector<float> v;
+ for (auto& i1: v1) {
+ for (auto& i2: v2) {
+ if (i1 != 0) {
+ v.push_back((i1-i2)/i1);
+ }
+ }
+ }
+ return v;
+ };
-/*
- //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("Q2_sim", "InclusiveKinematicsTruth.Q2")
+ .Define("Q2_rec", "InclusiveKinematicsElectron.Q2")
+ .Define("Q2_res", combinatorial_diff_ratio, {"Q2_sim", "Q2_rec"})
+ .Define("x_sim", "InclusiveKinematicsTruth.x")
+ .Define("x_rec", "InclusiveKinematicsElectron.x")
+ .Define("x_res", combinatorial_diff_ratio, {"x_sim", "x_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, {"mcparticles"})
- .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("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("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));
-*/
-
-
- auto c = new TCanvas();
+ auto h_Q2_sim = d0.Histo1D({"h_Q2_sim", "; GeV^2; counts", 100, -5, 25}, "Q2_sim");
+ auto h_Q2_rec = d0.Histo1D({"h_Q2_rec", "; GeV^2; counts", 100, -5, 25}, "Q2_rec");
+ auto h_Q2_res = d0.Histo1D({"h_Q2_res", "; ; counts", 100, -10, 10}, "Q2_res");
+ //x
+ auto h_x_sim = d0.Histo1D({"h_x_sim", "; ; counts", 100, 0, +1}, "x_sim");
+ auto h_x_rec = d0.Histo1D({"h_x_rec", "; ; counts", 100, 0, +1}, "x_rec");
+ auto h_x_res = d0.Histo1D({"h_x_res", "; ; counts", 100, -1, 1}, "x_res");
+
+ TFitResultPtr f_Q2_res = h_Q2_res->Fit("gaus", "S");
+ if (f_Q2_res == 0) f_Q2_res->Print("V");
+ TFitResultPtr f_x_res = h_x_res->Fit("gaus", "S");
+ if (f_x_res == 0) f_x_res->Print("V");
// Plot our histograms.
// TODO: to start I'm explicitly plotting the histograms, but want to
// factorize out the plotting code moving forward.
+
+ // Q2 comparison
{
- TCanvas c{"canvas", "canvas", 1200, 1200};
+ TCanvas c("c", "c", 1200, 1200);
c.cd();
- // gPad->SetLogx(false);
+ gPad->SetLogx(false);
gPad->SetLogy(true);
- //auto& h1 = *h_mom_sim;
- //auto& h2 = *h_mom_rec;
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 +118,94 @@ 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();
+ 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());
+ }
- c.Print(fmt::format("{}momentum.png", output_prefix).c_str());
+ // Q2 resolution
+ {
+ TCanvas c("c", "c", 1200, 1200);
+ c.cd();
+ gPad->SetLogx(false);
+ gPad->SetLogy(true);
+ auto& h1 = *h_Q2_res;
+ // histogram style
+ h1.SetLineColor(common_bench::plot::kMpBlue);
+ h1.SetLineWidth(2);
+ // axes
+ h1.GetXaxis()->CenterTitle();
+ h1.GetYaxis()->CenterTitle();
+ // draw everything
+ h1.DrawClone("hist");
+ // FIXME hardcoded beam configuration
+ common_bench::plot::draw_label(18, 275, detector);
+ c.Print(fmt::format("{}Q2resolution.png", output_prefix).c_str());
}
+
+ // x comparison
+ {
+ TCanvas c("c", "c", 1200, 1200);
+ c.cd();
+ gPad->SetLogx(true);
+ gPad->SetLogy(true);
+ auto& h1 = *h_x_sim;
+ auto& h2 = *h_x_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("{}x.png", output_prefix).c_str());
+ }
+
+ // x resolution
+ {
+ TCanvas c("c", "c", 1200, 1200);
+ c.cd();
+ gPad->SetLogx(true);
+ gPad->SetLogy(true);
+ auto& h1 = *h_x_res;
+ // histogram style
+ h1.SetLineColor(common_bench::plot::kMpBlue);
+ h1.SetLineWidth(2);
+ // axes
+ h1.GetXaxis()->CenterTitle();
+ h1.GetYaxis()->CenterTitle();
+ // draw everything
+ h1.DrawClone("hist");
+ // FIXME hardcoded beam configuration
+ common_bench::plot::draw_label(18, 275, detector);
+ c.Print(fmt::format("{}xresolution.png", output_prefix).c_str());
+ }
+
common_bench::write_test({dis_Q2_resolution}, fmt::format("{}dis_electrons.json", output_prefix));
return 0;
diff --git a/benchmarks/dis/config.yml b/benchmarks/dis/config.yml
index a3bb11959f757b7598afec95a8f9a8048baa4afa..4b50d5e6379b3fc5dd60f9cc68d18307fb826a43 100644
--- a/benchmarks/dis/config.yml
+++ b/benchmarks/dis/config.yml
@@ -1,5 +1,6 @@
dis:generate:
stage: generate
+ extends: .phy_benchmark
needs: ["common:detector"]
timeout: 1 hours
script:
@@ -7,12 +8,12 @@ dis:generate:
dis:process:
stage: process
+ extends: .phy_benchmark
needs: ["common:detector", "dis:generate"]
timeout: 2 hour
script:
- - echo "Temporarily disabling!!!"
- #- compile_analyses.py dis
- #- ./benchmarks/dis/dis.sh --config barrel --ebeam 18 --pbeam 275
+ - compile_analyses.py dis
+ - bash benchmarks/dis/dis.sh --config barrel --ebeam 18 --pbeam 275
retry:
max: 2
when:
diff --git a/benchmarks/dis/gen.sh b/benchmarks/dis/gen.sh
index f0e74c99e0cf97d70a68caba5a5c62a29e258cde..2714b1ab224bbd7df5206d7b948c6d247a1ea0ce 100755
--- a/benchmarks/dis/gen.sh
+++ b/benchmarks/dis/gen.sh
@@ -65,8 +65,8 @@ echo "Compiling benchmarks/dis/generator/pythia_dis.cxx ..."
g++ benchmarks/dis/generator/pythia_dis.cxx -o ${TMP_PATH}/pythia_dis \
-I/usr/local/include -I${LOCAL_PREFIX}/include \
-O2 -std=c++11 -pedantic -W -Wall -Wshadow -fPIC \
- -L/usr/local/lib -Wl,-rpath,/usr/local/lib -lpythia8 -ldl \
- -L/usr/local/lib -Wl,-rpath,/usr/local/lib -lHepMC3
+ $(pythia8-config --cxxflags --ldflags) \
+ $(HepMC3-config --cxxflags --ldflags)
if [[ "$?" -ne "0" ]] ; then
echo "ERROR compiling pythia"
exit 1