add Forward/Backward TOF

benchmarks/tof/scripts/gen_tof_forward_hits.cxx

+#include "HepMC3/GenEvent.h"
+#include "HepMC3/ReaderAscii.h"
+#include "HepMC3/WriterAscii.h"
+#include "HepMC3/Print.h"
+
+#include <iostream>
+#include <random>
+#include <cmath>
+#include <math.h>
+#include <TMath.h>
+
+using namespace HepMC3;
+
+/** Generate multiple electrons/positron tracks in the central region.
+ *  This is for testing detectors in the "barrel" region.
+ */
+void gen_tof_forward_hits(int n_events = 100,
+                         const char* out_fname = "tof_forward_hits.hepmc",
+                         int n_parts = 4)
+{
+  double cos_theta_min = std::cos( 1.0*(M_PI/180.0));
+  double cos_theta_max = std::cos(189.0*(M_PI/180.0));
+
+  WriterAscii hepmc_output(out_fname);
+  int events_parsed = 0;
+  GenEvent evt(Units::GEV, Units::MM);
+
+  // Random number generator
+  TRandom *r1 = new TRandom();
+
+  for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
+    // FourVector(px,py,pz,e,pdgid,status)
+    // type 4 is beam
+    // pdgid 11 - electron
+    // pdgid 111 - pi0
+    // pdgid 2212 - proton
+    Double_t p     = r1->Uniform(0, 8);
+    Double_t costh = -1*cos(12.033/180*3.1415926); //r1->Uniform(cos_theta_min, cos_theta_max);
+    Double_t th    = std::acos(costh);
+    Double_t phi   = r1->Uniform(0.0, 2.0 * M_PI);
+    Double_t px    = p * std::cos(phi) * std::sin(th);
+    Double_t py    = p * std::sin(phi) * std::sin(th);
+    Double_t pz    = p * std::cos(th);
+
+    for (int ip = 0; ip < n_parts; ip++) {
+      GenParticlePtr p1 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
+      GenParticlePtr p2 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
+
+      // Define momentum
+      phi   = r1->Uniform(0.0, 2.0 * M_PI);
+      px    = p * std::cos(phi) * std::sin(th);
+      py    = p * std::sin(phi) * std::sin(th);
+      pz    = p * std::cos(th);
+      GenParticlePtr p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.000511 * 0.000511))),
+                                                          ((ip % 2 == 0) ? 11 : -11), 1);
+
+      phi   = r1->Uniform(0.0, 2.0 * M_PI);
+      px    = p * std::cos(phi) * std::sin(th);
+      py    = p * std::sin(phi) * std::sin(th);
+      pz    = p * std::cos(th);
+      GenParticlePtr p4 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.938272 * 0.938272))),
+                                                          ((ip % 2 == 0) ? 2212 : -2212), 1);
+
+      phi   = r1->Uniform(0.0, 2.0 * M_PI);
+      px    = p * std::cos(phi) * std::sin(th);
+      py    = p * std::sin(phi) * std::sin(th);
+      pz    = p * std::cos(th);
+      GenParticlePtr p5 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.493677 * 0.493677))),
+                                                          ((ip % 2 == 0) ? 321 : -321), 1);
+      phi   = r1->Uniform(0.0, 2.0 * M_PI);
+      px    = p * std::cos(phi) * std::sin(th);
+      py    = p * std::sin(phi) * std::sin(th);
+      pz    = p * std::cos(th);
+      GenParticlePtr p6 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.139570 * 0.139570))),
+                                                          ((ip % 2 == 0) ? 211 : -211), 1);
+
+      GenVertexPtr v1 = std::make_shared<GenVertex>();
+      v1->add_particle_in(p1);
+      v1->add_particle_in(p2);
+      v1->add_particle_out(p3);
+      v1->add_particle_out(p4);
+      v1->add_particle_out(p5);
+      v1->add_particle_out(p6);
+      evt.add_vertex(v1);
+    }
+
+    if (events_parsed == 0) {
+      std::cout << "First event: " << std::endl;
+      Print::listing(evt);
+    }
+
+    hepmc_output.write_event(evt);
+    if (events_parsed % 10000 == 0) {
+      std::cout << "Event: " << events_parsed << std::endl;
+    }
+    evt.clear();
+  }
+  hepmc_output.close();
+  std::cout << "Events parsed and written: " << events_parsed << std::endl;
+}

benchmarks/tof/scripts/sim_tof_forward_hits.cxx

+#include "ROOT/RDataFrame.hxx"
+#include "TCanvas.h"
+#include "TH1D.h"
+#include "TLegend.h"
+#include "TProfile.h"
+
+#include <cstdlib>
+#include <iostream>
+
+R__LOAD_LIBRARY(libeicd.so)
+R__LOAD_LIBRARY(libDD4pod.so)
+
+#include <fmt/format.h>
+
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "dd4pod/TrackerHitCollection.h"
+#include "dd4pod/TrackerHitData.h"
+#include "eicd/ClusterCollection.h"
+#include "eicd/ClusterData.h"
+#include "eicd/TrackParametersCollection.h"
+#include "eicd/TrackerHitCollection.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;
+};
+
+std::vector<float> 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 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;
+};
+
+auto delta_p_over_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) / p1);
+    }
+  }
+  return res;
+};
+
+// Create dataframe with extra definitions for hit multiplicities
+ROOT::RDF::RNode add_subsystems(ROOT::RDF::RNode df, std::vector<std::pair<std::string, std::string>> hitcols) {
+  if (hitcols.empty()) {
+    return df;
+  }
+  const auto [name, collection] = hitcols.back();
+  std::cout << " --> registering subsystem " << collection << "\n";
+  auto df2 = df.Define("N_" + name, [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size(); }, {collection});
+  hitcols.pop_back();
+  return add_subsystems(df2, hitcols);
+};
+
+int sim_tof_forward_hits(const char* fname = "sim_tof_forward_hits.root") {
+
+  ROOT::EnableImplicitMT();
+  ROOT::RDataFrame df("events", fname);
+
+  // detect detector setup
+  std::string detector = "default";
+  if (const char* env_detector = std::getenv("JUGGLER_DETECTOR_VERSION")) {
+    if (detector.find("acadia") != std::string::npos) {
+      detector = "acadia";
+    } else if (detector.find("canyonlands") != std::string::npos) {
+      detector = "canyonlands";
+    }
+  }
+  std::cout << "sim_tof_forward_hits: detector set to " << detector << std::endl;
+
+  // minimal hit collection setup
+    std::vector<std::pair<std::string, std::string>> hitCollections{{"tof_forward", "TOFForwardHits"}};
+                                                                  //{"vtx_barrel", "VertexBarrelHits"},
+                                                                  //{"trk_barrel", "TrackerBarrelHits"},
+                                                                  //{"tof_barrel", "TOFBarrelHits"},
+                                                                  //{"trk_endcap", "TrackerEndcapHits"},
+                                                                  //{"gem_endcap", "GEMTrackerEndcapHits"}};
+
+  // append extra hit collections based on detector setup
+  if (detector == "acadia") {
+    hitCollections.push_back({"vtx_endcap", "VertexEndcapHits"});
+  } else if (detector == "canyonlands" || detector == "default") {
+    hitCollections.push_back({"mm_barrel", "MPGDTrackerBarrelHits"});
+  }
+
+  auto df0 = df.Define("isThrown", "mcparticles.genStatus == 1")
+                 .Define("thrownParticles", "mcparticles[isThrown]")
+                 .Define("thrownP", fourvec, {"thrownParticles"})
+                 .Define("p_thrown", momentum, {"thrownP"})
+                 .Define("theta_thrown", theta, {"thrownP"})
+                 .Define("theta0", "theta_thrown[0]");
+  auto df1 = add_subsystems(df0, hitCollections);
+
+  // get 1D and 2D histogram definitions
+  std::vector<decltype(df1.Histo2D({}, "", ""))> h2D;
+  std::vector<decltype(df1.Histo1D(""))> hnth;
+  std::vector<decltype(df1.Histo1D(""))> hn;
+
+  for (const auto& [name, col] : hitCollections) {
+    hnth.push_back(
+        df1.Histo1D({(name + "_nhits_vs_theta").c_str(), "; #theta [deg.]; ", 20, 0, 180}, "theta0", "N_" + name));
+    hn.push_back(df1.Histo1D({(name + "_nhits").c_str(), "; Nhits; #", 20, 0, 20}, "N_" + name));
+    h2D.push_back(df1.Histo2D({(name + "_x_vs_y").c_str(), "; x ; y ", 100, -600, 600, 100, -600, 600},
+                              col + ".position.x", col + ".position.y"));
+    h2D.push_back(df1.Histo2D({(name + "_x_vs_z").c_str(), "; z ; x ", 100, -2000, 2000, 100, -900, 900},
+                              col + ".position.z", col + ".position.x"));
+    h2D.push_back(df1.Histo2D({(name + "_time_vs_length").c_str(), "; time ; length ", 100, 0, 10, 100, 0, 1},
+           col + ".truth.time", col + ".length"));
+
+  }
+  auto hth = df1.Histo1D({"theta0", "; #theta [deg.]", 20, 0, 180}, "theta0");
+
+  // Draw multiplicity versus theta
+  {
+    TCanvas c;
+    THStack hs{"n_hits", "; #theta [deg.] "};
+    int idx = 1;
+    for (auto& h : hnth) {
+      h->Divide(&*hth);
+      h->SetLineColor(idx++);
+      hs.Add((TH1D*)h->Clone());
+    }
+    hs.Draw("nostack, hist");
+    c.BuildLegend();
+    c.SaveAs("results/tof/sim_tof_forward_hits_n_hits_vs_theta.png");
+    c.SaveAs("results/tof/sim_tof_forward_hits_n_hits_vs_theta.pdf");
+  }
+  // Draw nhits
+  {
+    TCanvas c;
+    THStack hs{"n_hits", "; NHits "};
+    int idx = 1;
+    for (auto& h : hn) {
+      h->SetLineColor(idx++);
+      hs.Add((TH1D*)h->Clone());
+    }
+    hs.Draw("nostack, hist");
+    c.BuildLegend();
+    c.SaveAs("results/tof/sim_tof_forward_hits_n_hits.png");
+    c.SaveAs("results/tof/sim_tof_forward_hits_n_hits.pdf");
+  }
+
+  // Draw total multiplicity versus theta
+  {
+    TCanvas c;
+    hth->DrawClone();
+    c.BuildLegend();
+    c.SaveAs("results/tof/sim_tof_forward_hits_theta.png");
+    c.SaveAs("results/tof/sim_tof_forward_hits_theta.pdf");
+  }
+
+  for (auto& h : h2D) {
+    TCanvas c;
+    h->DrawClone("colz");
+    c.SaveAs(fmt::format("results/tof/sim_tof_forward_hits_{}.png", h->GetName()).c_str());
+    c.SaveAs(fmt::format("results/tof/sim_tof_forward_hits_{}.pdf", h->GetName()).c_str());
+  }
+
+  return 0;
+}

benchmarks/tof/tof_forward_hits.sh

+#!/bin/bash
+
+#set -o errexit
+#set -o pipefail
+
+function print_the_help {
+  echo "USAGE: ${0}  [--sim-only]  "
+  echo "OPTIONS: "
+  echo "  --sim-only    Only run up to the simulation "
+  echo "  --analysis    Only run the analysis scripts "
+  exit 
+}
+
+ANALYSIS_ONLY=
+SIM_ONLY=
+POSITIONAL=()
+
+while [[ $# -gt 0 ]]
+do
+  key="$1"
+
+  case $key in
+    -h|--help)
+      shift # past argument
+      print_the_help
+      ;;
+    --sim-only)
+      shift # past argument
+      SIM_ONLY=1
+      ;;
+    --analysis)
+      shift # past argument
+      ANALYSIS_ONLY=1
+      ;;
+    *)    # unknown option
+      #POSITIONAL+=("$1") # save it in an array for later
+      echo "unknown option $1"
+      print_the_help
+      shift # past argument
+      ;;
+  esac
+done
+set -- "${POSITIONAL[@]}" # restore positional parameters
+
+print_env.sh
+
+## To run the reconstruction, we need the following global variables:
+## - JUGGLER_INSTALL_PREFIX:   Install prefix for Juggler (simu/recon)
+## - JUGGLER_DETECTOR:         the detector package we want to use for this benchmark
+## - JUGGLER_DETECTOR_VERSION: the detector package we want to use for this benchmark
+## - DETECTOR_PATH:            full path to the detector definitions
+##
+## You can ready options/env.sh for more in-depth explanations of the variables
+## and how they can be controlled.
+
+if [[ ! -n  "${JUGGLER_N_EVENTS}" ]] ; then 
+  export JUGGLER_N_EVENTS=100
+fi
+
+export JUGGLER_N_EVENTS=100
+export JUGGLER_FILE_NAME_TAG="tof_forward_hits"
+export JUGGLER_GEN_FILE="${LOCAL_DATA_PATH}/${JUGGLER_FILE_NAME_TAG}.hepmc"
+
+export JUGGLER_SIM_FILE="${LOCAL_DATA_PATH}/sim_${JUGGLER_FILE_NAME_TAG}.root"
+
+echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
+echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
+
+
+if [ -z "${ANALYSIS_ONLY}" ] ; then
+
+  echo "Generating Events"
+  ## generate the input events
+  root -b -q "benchmarks/tof/scripts/gen_tof_forward_hits.cxx(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+  if [[ "$?" -ne "0" ]] ; then
+    echo "ERROR running script"
+    exit 1
+  fi
+
+  echo "Running geant4 simulation"
+  ## run geant4 simulations
+  npsim --runType batch \
+    --part.minimalKineticEnergy 1000*GeV  \
+    -v WARNING \
+    --numberOfEvents ${JUGGLER_N_EVENTS} \
+    --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
+    --inputFiles  ${JUGGLER_FILE_NAME_TAG}.hepmc \
+    --outputFile  ${JUGGLER_SIM_FILE}
+  if [[ "$?" -ne "0" ]] ; then
+    echo "ERROR running script"
+    exit 1
+  fi
+
+fi
+
+
+if [ -z "${SIM_ONLY}" ] ; then
+
+  echo "Running analysis scripts"
+
+  mkdir -p results/tof
+  rootls -t ${JUGGLER_SIM_FILE}
+  root -b -q "benchmarks/tof/scripts/sim_tof_forward_hits.cxx(\"${JUGGLER_SIM_FILE}\")"
+  if [[ "$?" -ne "0" ]] ; then
+    echo "ERROR running root script"
+    exit 1
+  fi
+
+fi
+
+root_filesize=$(stat --format=%s "${JUGGLER_SIM_FILE}")
+echo "filesize" $filesize
+if [[ "${JUGGLER_N_EVENTS}" -lt "31000" ]] ; then
+  # file must be less than 10 MB to upload
+  if [[ "${root_filesize}" -lt "10000000000" ]] ; then
+    cp ${JUGGLER_SIM_FILE} results/.
+  fi
+fi
+ls -ltr ${JUGGLER_SIM_FILE}
+
+