diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 79dff3cb2cbf8a0d7a3e704675510136d323343b..a0ba22355bd5e7a684a6e42e32390219ef57eb37 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -64,6 +64,7 @@ common:detector:
include:
- local: 'benchmarks/ecal/config.yml'
+ - local: 'benchmarks/far_forward/config.yml'
- local: 'benchmarks/track_finding/config.yml'
- local: 'benchmarks/track_fitting/config.yml'
- local: 'benchmarks/tracking/config.yml'
diff --git a/benchmarks/far_forward/config.yml b/benchmarks/far_forward/config.yml
new file mode 100644
index 0000000000000000000000000000000000000000..05e79484c57dc0cc1bab6e24d440f26b5bfa8244
--- /dev/null
+++ b/benchmarks/far_forward/config.yml
@@ -0,0 +1,6 @@
+tracking_ff:
+ extends: .rec_benchmark
+ stage: run
+ timeout: 24 hours
+ script:
+ - bash benchmarks/far_forward/dummy_ff.sh
diff --git a/benchmarks/far_forward/dummy_ff.sh b/benchmarks/far_forward/dummy_ff.sh
new file mode 100644
index 0000000000000000000000000000000000000000..cc759747e7919cb86fda1c5039aa0e15c9732dce
--- /dev/null
+++ b/benchmarks/far_forward/dummy_ff.sh
@@ -0,0 +1,135 @@
+#!/bin/bash
+
+function print_the_help {
+ echo "USAGE: ${0} [--rec-only] "
+ echo " OPTIONS: "
+ echo " --rec-only only run gaudi reconstruction part"
+ exit
+}
+
+export PBEAM="100.0"
+export XANGLE="-0.025"
+
+REC_ONLY=
+ANALYSIS_ONLY=
+PARTICLE="proton"
+
+POSITIONAL=()
+while [[ $# -gt 0 ]]
+do
+ key="$1"
+
+ case $key in
+ -h|--help)
+ shift # past argument
+ print_the_help
+ ;;
+ --rec-only)
+ REC_ONLY=1
+ shift # past value
+ ;;
+ --ana-only)
+ ANALYSIS_ONLY=1
+ shift # past value
+ ;;
+ --particle)
+ PARTICLE=$2
+ shift # past key
+ shift # past value
+ ;;
+ *) # 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.
+export DETECTOR_PATH=${DETECTOR_PATH}
+
+if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then
+ export JUGGLER_N_EVENTS=100
+fi
+export JUGGLER_N_EVENTS=$(expr ${JUGGLER_N_EVENTS} \* 1)
+
+
+export JUGGLER_FILE_NAME_TAG="genparts"
+export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
+
+export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root"
+export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
+
+echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
+echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
+
+
+if [[ -z "${REC_ONLY}" && -z "${ANALYSIS_ONLY}" ]] ;
+then
+
+ python benchmarks/far_forward/scripts/gen_particles.py \
+ ${JUGGLER_FILE_NAME_TAG}.hepmc \
+ -n $JUGGLER_N_EVENTS -s 1 \
+ --thetamin 0.0 --thetamax 0.030 \
+ --pmin 40 --pmax 100 \
+ --particles ${PARTICLE}
+ ## generate the input events
+ 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
+
+rootls -t ${JUGGLER_SIM_FILE}
+
+if [[ -z "${ANALYSIS_ONLY}" ]] ;
+then
+ # Need to figure out how to pass file name to juggler from the commandline
+ xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv gaudirun.py benchmarks/far_forward/options/dummy_ff.py
+ if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running juggler"
+ exit 1
+ fi
+fi
+
+mkdir -p results/far_forward
+
+root -b -q "benchmarks/far_forward/scripts/rec_ff.cxx(\"${JUGGLER_REC_FILE}\")"
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running root script"
+ exit 1
+fi
+
+root_filesize=$(stat --format=%s "${JUGGLER_REC_FILE}")
+if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then
+ # file must be less than 10 MB to upload
+ if [[ "${root_filesize}" -lt "10000000" ]] ; then
+ cp ${JUGGLER_REC_FILE} results/.
+ fi
+fi
diff --git a/benchmarks/far_forward/options/dummy_ff.py b/benchmarks/far_forward/options/dummy_ff.py
new file mode 100644
index 0000000000000000000000000000000000000000..17202b8ffb071aa4429d26a88501557e4e44dc0f
--- /dev/null
+++ b/benchmarks/far_forward/options/dummy_ff.py
@@ -0,0 +1,58 @@
+from Gaudi.Configuration import *
+
+from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
+from GaudiKernel import SystemOfUnits as units
+
+detector_name = "athena"
+if "JUGGLER_DETECTOR" in os.environ :
+ detector_name = str(os.environ["JUGGLER_DETECTOR"])
+
+detector_path = ""
+if "DETECTOR_PATH" in os.environ :
+ detector_path = str(os.environ["DETECTOR_PATH"])
+
+# todo add checks
+input_sim_file = str(os.environ["JUGGLER_SIM_FILE"])
+output_rec_file = str(os.environ["JUGGLER_REC_FILE"])
+n_events = str(os.environ["JUGGLER_N_EVENTS"])
+
+podioevent = EICDataSvc("EventDataSvc", inputs=[input_sim_file], OutputLevel=WARNING)
+
+from Configurables import PodioInput
+from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
+
+from Configurables import Jug__Reco__DummyFarForwardParticles as DummyFFRec
+
+algorithms = [ ]
+
+podioinput = PodioInput("PodioReader",
+ collections=["mcparticles"])#, OutputLevel=DEBUG)
+algorithms.append( podioinput )
+
+## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
+copier = MCCopier("MCCopier",
+ inputCollection="mcparticles",
+ outputCollection="mcparticles2")
+algorithms.append( copier )
+
+ff_rec = DummyFFRec("DummyFFRec",
+ inputMCParticles="mcparticles",
+ outputParticles="ReconstructedParticles",
+ OutputLevel=DEBUG)
+algorithms.append( ff_rec )
+
+
+out = PodioOutput("out", filename=output_rec_file)
+out.outputCommands = ["keep *",
+ "drop mcparticles"
+ ]
+algorithms.append(out)
+
+ApplicationMgr(
+ TopAlg = algorithms,
+ EvtSel = 'NONE',
+ EvtMax = n_events,
+ ExtSvc = [podioevent,geo_service],
+ OutputLevel=WARNING
+ )
+
diff --git a/benchmarks/far_forward/scripts/gen_particles.py b/benchmarks/far_forward/scripts/gen_particles.py
new file mode 100644
index 0000000000000000000000000000000000000000..c8afd9acfb2f0ebaa309a06e693614dae275d6a3
--- /dev/null
+++ b/benchmarks/far_forward/scripts/gen_particles.py
@@ -0,0 +1,111 @@
+import os
+from pyHepMC3 import HepMC3 as hm
+import numpy as np
+import argparse
+
+
+PARTICLES = {
+ "pion0": (111, 0.1349766), # pi0
+ "pion+": (211, 0.13957018), # pi+
+ "pion-": (-211, 0.13957018), # pi-
+ "kaon0": (311, 0.497648), # K0
+ "kaon+": (321, 0.493677), # K+
+ "kaon-": (-321, 0.493677), # K-
+ "proton": (2212, 0.938272), # proton
+ "neutron": (2112, 0.939565), # neutron
+ "electron": (11, 0.51099895e-3), # electron
+ "positron": (-11, 0.51099895e-3),# positron
+ "photon": (22, 0), # photon
+}
+
+
+def gen_event(p, theta, phi, pid, mass, xangle):
+ evt = hm.GenEvent(momentum_unit=hm.Units.MomentumUnit.GEV, length_unit=hm.Units.LengthUnit.MM)
+ # final state
+ state = 1
+ e0 = np.sqrt(p*p + mass*mass)
+ pxc = p * np.cos(phi)*np.sin(theta)
+ pyc = p * np.sin(phi)*np.sin(theta)
+ pzc = p * np.cos(theta)
+
+ sa = np.sin(xangle)
+ ca = np.cos(xangle)
+
+ px = pxc * ca + pzc * sa
+ py = pyc
+ pz = - pxc * sa + pzc * ca
+
+ # beam
+ pbeam = hm.GenParticle(hm.FourVector(0, 0, 0, 0.938272), 2212, 4)
+ ebeam = hm.GenParticle(hm.FourVector(0, 0, e0, np.sqrt(e0*e0 + 0.511e-3*0.511e-3)), 11, 4)
+
+ # out particle
+ hout = hm.GenParticle(hm.FourVector(px, py, pz, e0), pid, state)
+
+ # vertex
+ vert = hm.GenVertex()
+ vert.add_particle_in(ebeam)
+ vert.add_particle_in(pbeam)
+ vert.add_particle_out(hout)
+ evt.add_vertex(vert)
+ return evt
+
+
+if __name__ == "__main__":
+ parser = argparse.ArgumentParser()
+
+ parser.add_argument('output', help='path to the output file')
+ parser.add_argument('-n', type=int, default=1000, dest='nev', help='number of events to generate')
+ parser.add_argument('-s', type=int, default=-1, dest='seed', help='seed for random generator')
+ parser.add_argument('--parray', type=str, default="", dest='parray',
+ help='an array of momenta in GeV, separated by \",\"')
+ parser.add_argument('--xangle', type=float, default=-0.025, dest='xangle', help="crossing angle")
+ parser.add_argument('--pmin', type=float, default=8.0, dest='pmin', help='minimum momentum in GeV')
+ parser.add_argument('--pmax', type=float, default=100.0, dest='pmax', help='maximum momentum in GeV')
+ parser.add_argument('--thetamin', type=float, default=-4, dest='thetamin', help='minimum angle (about the ion direction)')
+ parser.add_argument('--thetamax', type=float, default=4, dest='thetamax', help='maximum angle (about the ion direction)')
+ parser.add_argument('--phmin', type=float, default=0.0, dest='phmin', help='minimum angle in degree')
+ parser.add_argument('--phmax', type=float, default=360.0, dest='phmax', help='maximum angle in degree')
+ parser.add_argument('--particles', type=str, default='electron', dest='particles',
+ help='particle names, support {}'.format(list(PARTICLES.keys())))
+
+ args = parser.parse_args()
+
+ if args.seed < 0:
+ args.seed = os.environ.get('SEED', int.from_bytes(os.urandom(4), byteorder='big', signed=False))
+ print("Random seed is {}".format(args.seed))
+ np.random.seed(args.seed)
+
+ output = hm.WriterAscii(args.output);
+ if output.failed():
+ print("Cannot open file \"{}\"".format(args.output))
+ sys.exit(2)
+
+ # build particle info
+ parts = []
+ for pid in args.particles.split(','):
+ pid = pid.strip()
+ if pid not in PARTICLES.keys():
+ print('pid {:d} not found in dictionary, ignored.'.format(pid))
+ continue
+ parts.append(PARTICLES[pid])
+
+ # p values
+ pvals = np.random.uniform(args.pmin, args.pmax, args.nev) if not args.parray else \
+ np.random.choice([float(p.strip()) for p in args.parray.split(',')], args.nev)
+ thvals = np.arccos(np.random.uniform(np.cos(args.thetamin), np.cos(args.thetamax), args.nev))
+ phivals = np.random.uniform(args.phmin, args.phmax, args.nev)/180.*np.pi
+ partvals = [parts[i] for i in np.random.choice(len(parts), args.nev)]
+
+ count = 0
+ for p, theta, phi, (pid, mass) in zip(pvals, thvals, phivals, partvals):
+ if (count % 1000 == 0):
+ print("Generated {} events".format(count), end='\r')
+ evt = gen_event(p, theta, phi, pid, mass, args.xangle)
+ output.write_event(evt)
+ evt.clear()
+ count += 1
+
+ print("Generated {} events".format(args.nev))
+ output.close()
+
diff --git a/benchmarks/far_forward/scripts/rec_ff.cxx b/benchmarks/far_forward/scripts/rec_ff.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..1d29f2cf4659636c9d239f2ac85fe35155088edc
--- /dev/null
+++ b/benchmarks/far_forward/scripts/rec_ff.cxx
@@ -0,0 +1,147 @@
+#include "ROOT/RDataFrame.hxx"
+#include "TCanvas.h"
+#include "TH1D.h"
+#include "TLegend.h"
+#include "TProfile.h"
+#include <Math/Vector4D.h>
+
+R__LOAD_LIBRARY(fmt)
+#include <fmt/format.h>
+
+#include <iostream>
+
+R__LOAD_LIBRARY(libeicd)
+R__LOAD_LIBRARY(libDD4pod)
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "eicd/ReconstructedParticleCollection.h"
+
+using ROOT::RDataFrame;
+using namespace ROOT::VecOps;
+
+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].P());
+ }
+ 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;
+};
+
+std::vector<double> p_rec(const std::vector<eic::ReconstructedParticleData>& parts)
+{
+ std::vector<double> result;
+ for (const auto& part : parts) {
+ result.push_back(part.momentum);
+ }
+ return result;
+}
+std::vector<double> thx(const std::vector<eic::ReconstructedParticleData>& parts)
+{
+ std::vector<double> result;
+ for (const auto& part : parts) {
+ result.push_back(atan(part.p.x / part.p.z) * 1000);
+ }
+ return result;
+}
+std::vector<double> thy(const std::vector<eic::ReconstructedParticleData>& parts)
+{
+ std::vector<double> result;
+ for (const auto& part : parts) {
+ result.push_back(atan(part.p.y / part.p.z) * 1000);
+ }
+ return result;
+}
+std::vector<double> theta_colinear(const std::vector<eic::ReconstructedParticleData>& parts, const double xangle)
+{
+ std::vector<double> result;
+ for (const auto& part : parts) {
+ // undo the crossing angle for this one so things make more sense
+ const double ca = cos(-xangle);
+ const double sa = sin(-xangle);
+ const double pz = -part.p.x * sa + part.p.z * ca;
+ result.push_back(acos(pz / part.momentum) * 1000);
+ }
+ 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;
+};
+
+int rec_ff(const std::string& fname = "athena/rec_ff.root", const std::string& particle = "proton",
+ const double xangle = -0.25)
+{
+
+ ROOT::EnableImplicitMT();
+ ROOT::RDataFrame df("events", fname);
+
+ auto df0 = df.Define("isThrown", "mcparticles2.genStatus == 1")
+ .Define("thrownParticles", "mcparticles2[isThrown]")
+ .Define("xangle", fmt::format("{}", xangle))
+ .Define("thrownP", fourvec, {"thrownParticles"})
+ .Define("p_thrown", momentum, {"thrownP"})
+ .Define("p_rec", p_rec, {"ReconstructedParticles"})
+ .Define("delta_p", delta_p, {"p_rec", "p_thrown"})
+ .Define("delta_p_over_p0", delta_p_over_p, {"p_rec", "p_thrown"})
+ .Define("thx_mrad", thx, {"ReconstructedParticles"})
+ .Define("thy_mrad", thy, {"ReconstructedParticles"})
+ .Define("theta_mrad", theta_colinear, {"ReconstructedParticles", "xangle"});
+
+ auto h_angular = df0.Histo2D({"h_angular", ";#theta_{x} (mrad);#theta_{y} (mrad);", 60, -60, 60, 60, -60, 60},
+ "thx_mrad", "thy_mrad");
+ auto h_thetap =
+ df0.Histo2D({"h_thetap", ";#theta_{with ion} (mrad);P (GeV);", 60, 0, 100, 60, -30, 30}, "p_rec", "theta_mrad");
+
+ auto c = new TCanvas();
+
+ h_angular->DrawCopy("colz");
+ c->SaveAs(fmt::format("results/far_forward/far_forward_angular_{}.png", particle).c_str());
+ c->SaveAs(fmt::format("results/far_forward/far_forward_angular_{}.pdf", particle).c_str());
+
+ // -----------------------------------------------
+ h_thetap->DrawCopy("colz");
+ c->SaveAs(fmt::format("results/far_forward/far_forward_thetap_{}.png", particle).c_str());
+ c->SaveAs(fmt::format("results/far_forward/far_forward_thetap_{}.pdf", particle).c_str());
+
+ return 0;
+}