From a749745a4ad866f773dd1198951cae12ae20ef42 Mon Sep 17 00:00:00 2001
From: Wouter Deconinck <wouter.deconinck@umanitoba.ca>
Date: Mon, 11 Oct 2021 21:34:23 +0000
Subject: [PATCH] Add truth inclusive kinematics to reconstruction
---
.../{scripts => analysis}/dvcs_ps_gen.cxx | 6 +
benchmarks/dvcs/analysis/dvcs_tests.cxx | 54 +++++++
benchmarks/dvcs/config.yml | 1 +
benchmarks/dvcs/dvcs.sh | 2 +-
benchmarks/dvcs/scripts/dvcs_tests.cxx | 145 ------------------
benchmarks/single/config.yml | 1 +
.../u_omega/{scripts => analysis}/demo.cxx | 39 ++---
benchmarks/u_omega/config.yml | 1 +
benchmarks/u_omega/u_omega.sh | 2 +-
options/reconstruction.py | 8 +
10 files changed, 90 insertions(+), 169 deletions(-)
rename benchmarks/dvcs/{scripts => analysis}/dvcs_ps_gen.cxx (91%)
create mode 100644 benchmarks/dvcs/analysis/dvcs_tests.cxx
delete mode 100644 benchmarks/dvcs/scripts/dvcs_tests.cxx
rename benchmarks/u_omega/{scripts => analysis}/demo.cxx (74%)
diff --git a/benchmarks/dvcs/scripts/dvcs_ps_gen.cxx b/benchmarks/dvcs/analysis/dvcs_ps_gen.cxx
similarity index 91%
rename from benchmarks/dvcs/scripts/dvcs_ps_gen.cxx
rename to benchmarks/dvcs/analysis/dvcs_ps_gen.cxx
index c0642a7c..1631d930 100644
--- a/benchmarks/dvcs/scripts/dvcs_ps_gen.cxx
+++ b/benchmarks/dvcs/analysis/dvcs_ps_gen.cxx
@@ -1,3 +1,9 @@
+#include "TH2F.h"
+#include "TLorentzVector.h"
+#include "TGenPhaseSpace.h"
+
+#include <iostream>
+
void dvcs_ps_gen() {
double E_p = 100.0;
double M_p = 0.938;
diff --git a/benchmarks/dvcs/analysis/dvcs_tests.cxx b/benchmarks/dvcs/analysis/dvcs_tests.cxx
new file mode 100644
index 00000000..43695195
--- /dev/null
+++ b/benchmarks/dvcs/analysis/dvcs_tests.cxx
@@ -0,0 +1,54 @@
+#include <cmath>
+#include <iostream>
+#include <string>
+#include <vector>
+
+#include "ROOT/RDataFrame.hxx"
+#include "Math/Vector4D.h"
+#include "TCanvas.h"
+
+#include <nlohmann/json.hpp>
+using json = nlohmann::json;
+
+R__LOAD_LIBRARY(libfmt.so)
+#include "fmt/core.h"
+#include "fmt/color.h"
+
+R__LOAD_LIBRARY(libeicd.so)
+R__LOAD_LIBRARY(libDD4pod.so)
+
+#include "eicd/InclusiveKinematicsCollection.h"
+#include "eicd/ReconstructedParticleCollection.h"
+
+void dvcs_tests(const char* fname = "rec_dvcs.root"){
+
+ fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green), "Running DVCS analysis...\n");
+
+ // Run this in multi-threaded mode if desired
+ ROOT::EnableImplicitMT();
+ ROOT::RDataFrame df("events", fname);
+
+ auto df0 = df.Define("n_parts", "ReconstructedParticles.size()")
+ .Define("isQ2gt1", "InclusiveKinematicsTruth.Q2 > 1.0")
+ .Define("n_Q2gt1", "isQ2gt1.size()");
+
+ auto h_n_parts = df0.Histo1D({"h_n_parts", "; h_n_parts n", 10, 0, 10}, "n_parts");
+ auto h_Q2 = df0.Histo1D({"h_Q2", "; Q^{2} [GeV^{2}/c^{2}]", 100, 0, 30}, "InclusiveKinematicsTruth.Q2");
+ auto n_Q2gt1 = df0.Mean("n_Q2gt1");
+ auto n_parts = df0.Mean("n_parts");
+
+ // ---------------------------
+ // Do evaluation
+
+ auto c = new TCanvas();
+ h_Q2->DrawCopy();
+ c->SaveAs("results/dvcs/Q2.png");
+ c->SaveAs("results/dvcs/Q2.pdf");
+ fmt::print("{} DVCS events Q2>1\n",*n_Q2gt1);
+ fmt::print("{} tracks per event\n",*n_parts);
+
+ c = new TCanvas();
+ h_n_parts->DrawCopy();
+ c->SaveAs("results/dvcs/n_parts.png");
+ c->SaveAs("results/dvcs/n_parts.pdf");
+}
diff --git a/benchmarks/dvcs/config.yml b/benchmarks/dvcs/config.yml
index ee77a0f3..9b00dfeb 100644
--- a/benchmarks/dvcs/config.yml
+++ b/benchmarks/dvcs/config.yml
@@ -5,6 +5,7 @@ dvcs:process:
- s3
needs: ["common:detector"]
script:
+ - compile_analyses.py dvcs
- bash benchmarks/dvcs/dvcs.sh --all
dvcs:results:
diff --git a/benchmarks/dvcs/dvcs.sh b/benchmarks/dvcs/dvcs.sh
index 7a47ac14..c05c3486 100644
--- a/benchmarks/dvcs/dvcs.sh
+++ b/benchmarks/dvcs/dvcs.sh
@@ -145,7 +145,7 @@ if [[ -n "${DO_ANALYSIS}" || -n "${DO_ALL}" ]] ; then
mkdir -p results/dvcs
# here you can add as many scripts as you want.
- root -b -q "benchmarks/dvcs/scripts/dvcs_tests.cxx(\"${JUGGLER_REC_FILE}\")"
+ root -b -q "benchmarks/dvcs/analysis/dvcs_tests.cxx+(\"${JUGGLER_REC_FILE}\")"
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running root script"
exit 1
diff --git a/benchmarks/dvcs/scripts/dvcs_tests.cxx b/benchmarks/dvcs/scripts/dvcs_tests.cxx
deleted file mode 100644
index 2ebf9f45..00000000
--- a/benchmarks/dvcs/scripts/dvcs_tests.cxx
+++ /dev/null
@@ -1,145 +0,0 @@
-#include <cmath>
-#include <iostream>
-#include <string>
-#include <vector>
-
-#include "ROOT/RDataFrame.hxx"
-
-#include <nlohmann/json.hpp>
-using json = nlohmann::json;
-
-R__LOAD_LIBRARY(libfmt.so)
-#include "fmt/core.h"
-#include "fmt/color.h"
-
-R__LOAD_LIBRARY(libeicd.so)
-R__LOAD_LIBRARY(libDD4pod.so)
-
-#include "dd4pod/Geant4ParticleCollection.h"
-#include "eicd/TrackParametersCollection.h"
-#include "eicd/ClusterCollection.h"
-#include "eicd/ReconstructedParticleCollection.h"
-#include "eicd/BasicParticleCollection.h"
-
-using ROOT::RDataFrame;
-using namespace ROOT::VecOps;
-
-auto p_track = [](std::vector<eic::TrackParametersData> const& in) {
- std::vector<double> result;
- for (size_t i = 0; i < in.size(); ++i) {
- result.push_back(std::abs(1.0/(in[i].qOverP)));
- }
- return result;
-};
-
-
-auto pt = [](std::vector<dd4pod::Geant4ParticleData> const& in){
- std::vector<float> result;
- for (size_t i = 0; i < in.size(); ++i) {
- result.push_back(std::sqrt(in[i].ps.x * in[i].ps.x + in[i].ps.y * in[i].ps.y));
- }
- return result;
-};
-
-auto momentum = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
- std::vector<double> result;
- for (size_t i = 0; i < in.size(); ++i) {
- result.push_back(in[i].E());
- }
- return result;
-};
-auto theta = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
- std::vector<double> result;
- for (size_t i = 0; i < in.size(); ++i) {
- result.push_back(in[i].Theta()*180/M_PI);
- }
- return result;
-};
-auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
- std::vector<ROOT::Math::PxPyPzMVector> result;
- ROOT::Math::PxPyPzMVector lv;
- for (size_t i = 0; i < in.size(); ++i) {
- lv.SetCoordinates(in[i].ps.x, in[i].ps.y, in[i].ps.z, in[i].mass);
- result.push_back(lv);
- }
- return result;
-};
-auto recfourvec = [](ROOT::VecOps::RVec<eic::ReconstructedParticleData> const& in) {
- std::vector<ROOT::Math::PxPyPzMVector> result;
- ROOT::Math::PxPyPzMVector lv;
- for (size_t i = 0; i < in.size(); ++i) {
- lv.SetCoordinates(in[i].p.x, in[i].p.y, in[i].p.z, in[i].mass);
- result.push_back(lv);
- }
- return result;
-};
-
-auto delta_p = [](const std::vector<double>& tracks, const std::vector<double>& thrown) {
- std::vector<double> res;
- for (const auto& p1 : thrown) {
- for (const auto& p2 : tracks) {
- res.push_back(p1 - p2);
- }
- }
- return res;
-};
-
-
-void dvcs_tests(const char* fname = "rec_dvcs.root"){
-
- fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green), "Running DVCS analysis...\n");
-
- // Run this in multi-threaded mode if desired
- ROOT::EnableImplicitMT();
- ROOT::RDataFrame df("events", fname);
-
- using ROOT::Math::PxPyPzMVector;
- PxPyPzMVector p_ebeam = {0,0,-10, 0.000511};
- PxPyPzMVector p_pbeam = {0,0,275, 0.938 };
-
- auto eprime = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
- for(const auto& p : in){
- if(p.pdgID == 11 ) {
- return PxPyPzMVector(p.ps.x,p.ps.y,p.ps.z,p.mass);
- }
- }
- return PxPyPzMVector(0,0,0,0);
- };
- auto q_vec = [=](PxPyPzMVector const& p) {
- return p_ebeam - p;
- };
-
- auto df0 = df.Define("isThrown", "mcparticles.genStatus == 1")
- .Define("thrownParticles", "mcparticles[isThrown]")
- .Define("thrownP", fourvec, {"thrownParticles"})
- .Define("recP", recfourvec, {"ReconstructedParticles"})
- .Define("NPart", "recP.size()")
- .Define("p_thrown", momentum, {"thrownP"})
- .Define("nTracks", "outputTrackParameters.size()")
- .Define("p_track", p_track, {"outputTrackParameters"})
- .Define("delta_p",delta_p, {"p_track", "p_thrown"})
- .Define("eprime", eprime, {"thrownParticles"})
- .Define("q", q_vec, {"eprime"})
- .Define("Q2", "-1.0*(q.Dot(q))");
-
- auto h_n_dummy = df0.Histo1D({"h_n_part", "; h_n_part n", 10, 0, 10}, "NPart");
- auto h_Q2 = df0.Histo1D({"h_Q2", "; Q^{2} [GeV^{2}/c^{2}]", 100, 0, 30}, "Q2");
- auto n_Q2 = df0.Filter("Q2>1").Count();
- auto n_tracks = df0.Mean("nTracks");
-
- // ---------------------------
- // Do evaluation
-
- auto c = new TCanvas();
- h_Q2->DrawCopy();
- c->SaveAs("results/dvcs/Q2.png");
- c->SaveAs("results/dvcs/Q2.pdf");
- fmt::print("{} DVCS events\n",*n_Q2);
- fmt::print("{} tracks per event\n",*n_tracks);
-
- c = new TCanvas();
- h_n_dummy->DrawCopy();
- c->SaveAs("results/dvcs/n_dummy.png");
- //c->SaveAs("results/dvcs/n_dummy.pdf");
-
-}
diff --git a/benchmarks/single/config.yml b/benchmarks/single/config.yml
index 35ed2a0d..a84da3f8 100644
--- a/benchmarks/single/config.yml
+++ b/benchmarks/single/config.yml
@@ -3,4 +3,5 @@ single:process:
timeout: 24 hours
stage: process
script:
+ - compile_analyses.py single
- bash benchmarks/single/single.sh e-_1GeV_45to135deg
diff --git a/benchmarks/u_omega/scripts/demo.cxx b/benchmarks/u_omega/analysis/demo.cxx
similarity index 74%
rename from benchmarks/u_omega/scripts/demo.cxx
rename to benchmarks/u_omega/analysis/demo.cxx
index afdc0d9a..ad63c895 100644
--- a/benchmarks/u_omega/scripts/demo.cxx
+++ b/benchmarks/u_omega/analysis/demo.cxx
@@ -4,6 +4,8 @@
#include <vector>
#include "ROOT/RDataFrame.hxx"
+#include "Math/Vector4D.h"
+#include "TCanvas.h"
#include <nlohmann/json.hpp>
using json = nlohmann::json;
@@ -20,6 +22,7 @@ R__LOAD_LIBRARY(libDD4pod.so)
#include "eicd/ClusterCollection.h"
#include "eicd/ReconstructedParticleCollection.h"
#include "eicd/BasicParticleCollection.h"
+#include "eicd/InclusiveKinematicsCollection.h"
using ROOT::RDataFrame;
using namespace ROOT::VecOps;
@@ -109,23 +112,15 @@ void demo(const char* fname = "rec_dvcs.root"){
return p_ebeam - p;
};
- auto df0 = df.Define("isThrown", "mcparticles.genStatus == 1")
- .Define("thrownParticles", "mcparticles[isThrown]")
- .Define("thrownP", fourvec, {"thrownParticles"})
- .Define("recP", recfourvec, {"ReconstructedParticles"})
- .Define("NPart", "recP.size()")
- .Define("p_thrown", momentum, {"thrownP"})
- .Define("nTracks", "outputTrackParameters.size()")
- .Define("p_track", p_track, {"outputTrackParameters"})
- .Define("delta_p",delta_p, {"p_track", "p_thrown"})
- .Define("eprime", eprime, {"thrownParticles"})
- .Define("q", q_vec, {"eprime"})
- .Define("Q2", "-1.0*(q.Dot(q))");
-
- auto h_n_dummy = df0.Histo1D({"h_n_part", "; h_n_part n", 10, 0, 10}, "NPart");
- auto h_Q2 = df0.Histo1D({"h_Q2", "; Q^{2} [GeV^{2}/c^{2}]", 100, 0, 30}, "Q2");
- auto n_Q2 = df0.Filter("Q2>1").Count();
- auto n_tracks = df0.Mean("nTracks");
+ auto df0 = df.Define("n_parts", "ReconstructedParticles.size()")
+ .Define("isQ2gt1", "InclusiveKinematicsTruth.Q2 > 1.0")
+ .Define("n_Q2gt1", "isQ2gt1.size()");
+
+
+ auto h_n_parts = df0.Histo1D({"h_n_parts", "; h_n_parts n", 10, 0, 10}, "n_parts");
+ auto h_Q2 = df0.Histo1D({"h_Q2", "; Q^{2} [GeV^{2}/c^{2}]", 100, 0, 30}, "InclusiveKinematicsTruth.Q2");
+ auto n_Q2gt1 = df0.Mean("n_Q2gt1");
+ auto n_parts = df0.Mean("n_parts");
// ---------------------------
// Do evaluation
@@ -134,12 +129,12 @@ void demo(const char* fname = "rec_dvcs.root"){
h_Q2->DrawCopy();
c->SaveAs("results/u_omega/Q2.png");
c->SaveAs("results/u_omega/Q2.pdf");
- fmt::print("{} u_omega events\n",*n_Q2);
- fmt::print("{} tracks per event\n",*n_tracks);
+ fmt::print("{} u_omega events Q2>1\n",*n_Q2gt1);
+ fmt::print("{} tracks per event\n",*n_parts);
c = new TCanvas();
- h_n_dummy->DrawCopy();
- c->SaveAs("results/u_omega/n_dummy.png");
- //c->SaveAs("results/u_omega/n_dummy.pdf");
+ h_n_parts->DrawCopy();
+ c->SaveAs("results/u_omega/n_parts.png");
+ c->SaveAs("results/u_omega/n_parts.pdf");
}
diff --git a/benchmarks/u_omega/config.yml b/benchmarks/u_omega/config.yml
index 9e4925b6..40d5f4be 100644
--- a/benchmarks/u_omega/config.yml
+++ b/benchmarks/u_omega/config.yml
@@ -5,6 +5,7 @@ u_omega:process:
- s3
needs: ["common:detector"]
script:
+ - compile_analyses.py u_omega
- bash benchmarks/u_omega/u_omega.sh --all
u_omega:results:
diff --git a/benchmarks/u_omega/u_omega.sh b/benchmarks/u_omega/u_omega.sh
index 6d051589..23bbf403 100644
--- a/benchmarks/u_omega/u_omega.sh
+++ b/benchmarks/u_omega/u_omega.sh
@@ -144,7 +144,7 @@ if [[ -n "${DO_ANALYSIS}" || -n "${DO_ALL}" ]] ; then
mkdir -p "results/${FILE_NAME_TAG}"
# here you can add as many scripts as you want.
- root -b -q "benchmarks/${FILE_NAME_TAG}/scripts/demo.cxx(\"${JUGGLER_REC_FILE}\")"
+ root -b -q "benchmarks/${FILE_NAME_TAG}/analysis/demo.cxx+(\"${JUGGLER_REC_FILE}\")"
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running root script"
exit 1
diff --git a/options/reconstruction.py b/options/reconstruction.py
index 66b449dc..927a8ffc 100644
--- a/options/reconstruction.py
+++ b/options/reconstruction.py
@@ -62,6 +62,7 @@ from Configurables import Jug__Fast__SmearedFarForwardParticles as FFSmearedPart
from Configurables import Jug__Fast__MatchClusters as MatchClusters
from Configurables import Jug__Fast__ClusterMerger as ClusterMerger
from Configurables import Jug__Fast__TruthEnergyPositionClusterMerger as EnergyPositionClusterMerger
+from Configurables import Jug__Fast__InclusiveKinematicsTruth as InclusiveKinematicsTruth
from Configurables import Jug__Digi__PhotoMultiplierDigi as PhotoMultiplierDigi
from Configurables import Jug__Digi__CalorimeterHitDigi as CalHitDigi
@@ -129,6 +130,13 @@ dummy = MC2DummyParticle("dummy",
smearing=0)
algorithms.append(dummy)
+# Truth level kinematics
+truth_incl_kin = InclusiveKinematicsTruth("truth_incl_kin",
+ inputMCParticles="mcparticles",
+ outputData="InclusiveKinematicsTruth"
+)
+algorithms.append(truth_incl_kin)
+
# Crystal Endcap Ecal
ce_ecal_daq = dict(
dynamicRangeADC=5.*units.GeV,
--
GitLab