start work on DummyFF benchmark

.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'

benchmarks/far_forward/config.yml

+tracking_ff:
+  extends: .rec_benchmark
+  stage: run
+  timeout: 24 hours
+  script:
+    - bash benchmarks/far_forward/dummy_ff.sh

benchmarks/far_forward/dummy_ff.sh

+#!/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

benchmarks/far_forward/options/dummy_ff.py

+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
+ )
+

benchmarks/far_forward/scripts/gen_particles.py

+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()
+

benchmarks/far_forward/scripts/rec_ff.cxx

+#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;
+}