From 8d3a8bd85a04943b3eb28d26b5addead8d6506cc Mon Sep 17 00:00:00 2001
From: Wouter Deconinck <wdconinc@gmail.com>
Date: Thu, 13 Jan 2022 01:03:31 +0000
Subject: [PATCH] tcs_tests.cxx: plot theta in mrad for various FF detectors
---
benchmarks/tcs/analysis/tcs_tests.cxx | 49 ++++++++++++++++++++++++---
benchmarks/tcs/tcs.sh | 3 +-
2 files changed, 46 insertions(+), 6 deletions(-)
diff --git a/benchmarks/tcs/analysis/tcs_tests.cxx b/benchmarks/tcs/analysis/tcs_tests.cxx
index 5fe49b7b..f833b43a 100644
--- a/benchmarks/tcs/analysis/tcs_tests.cxx
+++ b/benchmarks/tcs/analysis/tcs_tests.cxx
@@ -6,6 +6,7 @@
#include "ROOT/RDataFrame.hxx"
#include "Math/Vector4D.h"
#include "TCanvas.h"
+#include "TLegend.h"
#include <nlohmann/json.hpp>
using json = nlohmann::json;
@@ -28,13 +29,33 @@ void tcs_tests(const char* fname = "rec_tcs.root"){
ROOT::EnableImplicitMT();
ROOT::RDataFrame df("events", fname);
+ auto ff_theta_mrad = [] (
+ const std::vector<eic::ReconstructedParticleData>& v,
+ const int16_t status
+ ) -> std::vector<float> {
+ std::vector<float> theta;
+ for (const auto& p: v)
+ if (p.status == status)
+ theta.push_back(1000. * p.direction.theta);
+ return theta;
+ };
+
auto df0 = df.Define("n_parts", "ReconstructedParticles.size()")
.Define("isQ2gt1", "InclusiveKinematicsTruth.Q2 > 1.0")
- .Define("n_Q2gt1", "isQ2gt1.size()");
+ .Define("n_Q2gt1", "isQ2gt1.size()")
+ .Define("ff_theta_mrad_B0", [&](const std::vector<eic::ReconstructedParticleData>& v){return ff_theta_mrad(v,1);}, {"ReconstructedFFParticles"})
+ .Define("ff_theta_mrad_RP", [&](const std::vector<eic::ReconstructedParticleData>& v){return ff_theta_mrad(v,2);}, {"ReconstructedFFParticles"})
+ .Define("ff_theta_mrad_OMD",[&](const std::vector<eic::ReconstructedParticleData>& v){return ff_theta_mrad(v,3);}, {"ReconstructedFFParticles"})
+ .Define("ff_theta_mrad_ZDC",[&](const std::vector<eic::ReconstructedParticleData>& v){return ff_theta_mrad(v,4);}, {"ReconstructedFFParticles"})
+ ;
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 h_ff_status = df0.Histo1D({"h_ff_status", "; FF status", 10, -0.5, 9.5}, "ReconstructedFFParticles.status");
+ auto h_FF = df0.Histo1D({"h_FF", "; FF status", 10, -0.5, 9.5}, "ReconstructedFFParticles.status");
+ auto h_FF_B0 = df0.Histo1D({"h_FF_B0", "; FF B0 Theta [mrad]", 100, 0.0, 25.0}, "ff_theta_mrad_B0");
+ auto h_FF_RP = df0.Histo1D({"h_FF_RP", "; FF RP Theta [mrad]", 100, 0.0, 25.0}, "ff_theta_mrad_RP");
+ auto h_FF_OMD = df0.Histo1D({"h_FF_OMD", "; FF OMD Theta [mrad]", 100, 0.0, 25.0}, "ff_theta_mrad_OMD");
+ auto h_FF_ZDC = df0.Histo1D({"h_FF_ZDC", "; FF ZDC Theta [mrad]", 100, 0.0, 25.0}, "ff_theta_mrad_ZDC");
auto n_Q2gt1 = df0.Mean("n_Q2gt1");
auto n_parts = df0.Mean("n_parts");
@@ -54,7 +75,25 @@ void tcs_tests(const char* fname = "rec_tcs.root"){
c->SaveAs("results/tcs/n_parts.pdf");
c = new TCanvas();
- h_ff_status->DrawCopy();
- c->SaveAs("results/tcs/ff_status.png");
- c->SaveAs("results/tcs/ff_status.pdf");
+ h_FF->DrawCopy();
+ c->SaveAs("results/tcs/ff.png");
+ c->SaveAs("results/tcs/ff.pdf");
+
+ c = new TCanvas();
+ h_FF_B0->SetLineColor(1);
+ h_FF_B0->DrawCopy();
+ h_FF_RP->SetLineColor(2);
+ h_FF_RP->DrawCopy("same");
+ h_FF_OMD->SetLineColor(3);
+ h_FF_OMD->DrawCopy("same");
+ h_FF_ZDC->SetLineColor(4);
+ h_FF_ZDC->DrawCopy("same");
+ auto legend = new TLegend(0.1,0.7,0.48,0.9);
+ legend->AddEntry(h_FF_B0.GetPtr(), "B0", "L");
+ legend->AddEntry(h_FF_RP.GetPtr(), "RP", "L");
+ legend->AddEntry(h_FF_OMD.GetPtr(), "OMD", "L");
+ legend->AddEntry(h_FF_ZDC.GetPtr(), "ZDC", "L");
+ legend->Draw();
+ c->SaveAs("results/tcs/ff_theta.png");
+ c->SaveAs("results/tcs/ff_theta.pdf");
}
diff --git a/benchmarks/tcs/tcs.sh b/benchmarks/tcs/tcs.sh
index d2c6e879..6c858947 100644
--- a/benchmarks/tcs/tcs.sh
+++ b/benchmarks/tcs/tcs.sh
@@ -111,7 +111,7 @@ echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
## Step 1. Get the data
if [[ -n "${DATA_INIT}" || -n "${DO_ALL}" ]] ; then
mc -C . config host add S3 https://dtn01.sdcc.bnl.gov:9000 $S3_ACCESS_KEY $S3_SECRET_KEY
- mc -C . cat --insecure ${DATA_URL} | gunzip -c | head -n 1004 | sanitize_hepmc3 > "${JUGGLER_MC_FILE}"
+ mc -C . cat --insecure ${DATA_URL} | gunzip -c | head -n $((20+10*JUGGLER_N_EVENTS)) | sanitize_hepmc3 > "${JUGGLER_MC_FILE}"
if [[ "$?" -ne "0" ]] ; then
echo "Failed to download hepmc file"
exit 1
@@ -135,6 +135,7 @@ if [[ -n "${DO_SIM}" || -n "${DO_ALL}" ]] ; then
fi
### Step 3. Run the reconstruction (juggler)
+export PBEAM
if [[ -n "${DO_REC}" || -n "${DO_ALL}" ]] ; then
for rec in options/*.py ; do
unset tag
--
GitLab