tcs_tests.cxx: plot theta in mrad for various FF detectors

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");
 }

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