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

Add truth inclusive kinematics to reconstruction

parent ef57ccf7
No related branches found
No related tags found
1 merge request!99Add truth inclusive kinematics to reconstruction
Hide whitespace changes
Inline Side-by-side
Showing
with 90 additions and 169 deletions
benchmarks/dvcs/scripts/dvcs_ps_gen.cxx → benchmarks/dvcs/analysis/dvcs_ps_gen.cxx
+ 6
− 0
View file @ a749745a
#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;
......
benchmarks/dvcs/analysis/dvcs_tests.cxx 0 → 100644
+ 54
− 0
View file @ a749745a
#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");
}
benchmarks/dvcs/config.yml
+ 1
− 0
View file @ a749745a
......@@ -5,6 +5,7 @@ dvcs:process:
- s3
needs: ["common:detector"]
script:
- compile_analyses.py dvcs
- bash benchmarks/dvcs/dvcs.sh --all
dvcs:results:
......
benchmarks/dvcs/dvcs.sh
+ 1
− 1
View file @ a749745a
......@@ -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
......
benchmarks/dvcs/scripts/dvcs_tests.cxx deleted 100644 → 0
+ 0
− 145
View file @ ef57ccf7
#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");
}
benchmarks/single/config.yml
+ 1
− 0
View file @ a749745a
......@@ -3,4 +3,5 @@ single:process:
timeout: 24 hours
stage: process
script:
- compile_analyses.py single
- bash benchmarks/single/single.sh e-_1GeV_45to135deg
benchmarks/u_omega/scripts/demo.cxx → benchmarks/u_omega/analysis/demo.cxx
+ 17
− 22
View file @ a749745a
......@@ -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");
}
benchmarks/u_omega/config.yml
+ 1
− 0
View file @ a749745a
......@@ -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:
......
benchmarks/u_omega/u_omega.sh
+ 1
− 1
View file @ a749745a
......@@ -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
......
options/reconstruction.py
+ 8
− 0
View file @ a749745a
......@@ -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,
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment