diff --git a/.clang-format b/.clang-format
index 05b10dc8ce7802fa782322aca13d6d67d14c6cf8..c43c6165c7c67310628bd74e3b38b2c0a5a072a7 100644
--- a/.clang-format
+++ b/.clang-format
@@ -44,7 +44,7 @@ BreakConstructorInitializersBeforeComma: false
BreakConstructorInitializers: BeforeColon
BreakAfterJavaFieldAnnotations: false
BreakStringLiterals: true
-ColumnLimit: 120
+ColumnLimit: 100
CommentPragmas: '^ IWYU pragma:'
CompactNamespaces: false
ConstructorInitializerAllOnOneLineOrOnePerLine: false
diff --git a/.rootlogon.C b/.rootlogon.C
index 808058da0a709843cf4aa848d950cd2bbb64e6f1..7055e9129f770742f8288bd5196505582a2ec46c 100644
--- a/.rootlogon.C
+++ b/.rootlogon.C
@@ -7,6 +7,7 @@
// setup a local build directory so we don't polute our source code with
// ROOT dictionaries etc.
+ gSystem->Exec("mkdir -p /tmp/root_build");
gSystem->SetBuildDir("/tmp/root_build");
// style definition based off the ATLAS style
diff --git a/accelerator b/accelerator
new file mode 160000
index 0000000000000000000000000000000000000000..f3ff428e3b926a41e95beaa984d8dc05cec37cc7
--- /dev/null
+++ b/accelerator
@@ -0,0 +1 @@
+Subproject commit f3ff428e3b926a41e95beaa984d8dc05cec37cc7
diff --git a/benchmarks/dis/analysis/dis.h b/benchmarks/dis/analysis/dis.h
new file mode 100644
index 0000000000000000000000000000000000000000..0e88a01a008726cd6c6f2301fa07aef9d0da03ce
--- /dev/null
+++ b/benchmarks/dis/analysis/dis.h
@@ -0,0 +1,26 @@
+#ifndef DVMP_H
+#define DVMP_H
+
+#include <util.h>
+
+#include <algorithm>
+#include <cmath>
+#include <exception>
+#include <fmt/core.h>
+#include <limits>
+#include <string>
+#include <vector>
+
+#include <Math/Vector4D.h>
+
+// Additional utility functions for DVMP benchmarks. Where useful, these can be
+// promoted to the top-level util library
+namespace util {
+
+ // ADD EXTRA DIS UTILTIY FUNCTIONS HERE
+
+ //=========================================================================================================
+
+} // namespace util
+
+#endif
diff --git a/benchmarks/dis/analysis/dis_electrons.cxx b/benchmarks/dis/analysis/dis_electrons.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..02dc7b0dac383580c671aa98591f0a64cca3671b
--- /dev/null
+++ b/benchmarks/dis/analysis/dis_electrons.cxx
@@ -0,0 +1,116 @@
+#include "dis.h"
+#include "plot.h"
+
+#include <benchmark.h>
+#include <mt.h>
+#include <util.h>
+
+#include "ROOT/RDataFrame.hxx"
+#include <cmath>
+#include <fmt/color.h>
+#include <fmt/core.h>
+#include <fstream>
+#include <iostream>
+#include <nlohmann/json.hpp>
+#include <string>
+#include <vector>
+
+int dis_electrons(const std::string& config_name)
+{
+ // read our configuration
+ std::ifstream config_file{config_name};
+ nlohmann::json config;
+ config_file >> config;
+
+ const std::string rec_file = config["rec_file"];
+ const std::string detector = config["detector"];
+ std::string output_prefix = config["output_prefix"];
+ const std::string test_tag = config["test_tag"];
+
+ fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green),
+ "Running DIS electron analysis...\n");
+ fmt::print(" - Detector package: {}\n", detector);
+ fmt::print(" - input file: {}\n", rec_file);
+ fmt::print(" - output prefix: {}\n", output_prefix);
+ fmt::print(" - test tag: {}\n", test_tag);
+
+ // create our test definition
+ // test_tag
+ eic::util::Test dis_Q2_resolution{
+ {{"name", fmt::format("{}_Q2_resolution", test_tag)},
+ {"title", "DIS Q2 resolution"},
+ {"description",
+ fmt::format("DIS Q2 resolution with {}, estimated using a Gaussian fit.", detector)},
+ {"quantity", "resolution (in %)"},
+ {"target", "0.1"}}};
+
+ // Run this in multi-threaded mode if desired
+ ROOT::EnableImplicitMT(kNumThreads);
+
+ const double electron_mass = util::get_pdg_mass("electron");
+
+ // Ensure our output prefix always ends on a dot, a slash or a dash
+ // Necessary when generating output plots
+ if (output_prefix.back() != '.' && output_prefix.back() != '/' && output_prefix.back() != '-') {
+ output_prefix += "-";
+ }
+
+ ROOT::RDataFrame d("events", rec_file);
+
+ // utility lambda functions to bind the reconstructed particle type
+ // (as we have no PID yet)
+ auto momenta_from_tracking =
+ [electron_mass](const std::vector<eic::TrackParametersData>& tracks) {
+ return util::momenta_from_tracking(tracks, electron_mass);
+ };
+
+ auto d0 = d.Define("p_rec", momenta_from_tracking, {"outputTrackParameters"})
+ .Define("N", "p_rec.size()")
+ .Define("p_sim", util::momenta_from_simulation, {"mcparticles2"})
+ .Define("mom_sim", util::mom, {"p_sim"})
+ .Define("mom_rec", util::mom, {"p_rec"});
+
+ auto h_mom_sim = d0.Histo1D({"h_mom_sim", "; GeV; counts", 100, 0, 50}, "mom_sim");
+ auto h_mom_rec = d0.Histo1D({"h_mom_rec", "; GeV; counts", 100, 0, 50}, "mom_rec");
+
+ auto c = new TCanvas();
+
+ // Plot our histograms.
+ // TODO: to start I'm explicitly plotting the histograms, but want to
+ // factorize out the plotting code moving forward.
+ {
+ TCanvas c{"canvas", "canvas", 1200, 1200};
+ c.cd();
+ // gPad->SetLogx(false);
+ gPad->SetLogy(true);
+ auto& h1 = *h_mom_sim;
+ auto& h2 = *h_mom_rec;
+ // histogram style
+ h1.SetLineColor(plot::kMpBlue);
+ h1.SetLineWidth(2);
+ h2.SetLineColor(plot::kMpOrange);
+ h2.SetLineWidth(2);
+ // axes
+ h1.GetXaxis()->CenterTitle();
+ h1.GetYaxis()->CenterTitle();
+ // draw everything
+ h1.DrawClone("hist");
+ h2.DrawClone("hist same");
+ // FIXME hardcoded beam configuration
+ plot::draw_label(18, 275, detector);
+ TText* tptr1;
+ auto t1 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
+ t1->SetFillColorAlpha(kWhite, 0);
+ t1->SetTextFont(43);
+ t1->SetTextSize(25);
+ tptr1 = t1->AddText("simulated");
+ tptr1->SetTextColor(plot::kMpBlue);
+ tptr1 = t1->AddText("reconstructed");
+ tptr1->SetTextColor(plot::kMpOrange);
+ t1->Draw();
+
+ c.Print(fmt::format("{}momentum.png", output_prefix).c_str());
+ }
+
+ return 0;
+}
diff --git a/benchmarks/dis/analysis/rec_dis_electrons.cxx b/benchmarks/dis/analysis/rec_dis_electrons.cxx
deleted file mode 100644
index b7227a1ddb33c2564fd160e2124da774a9711161..0000000000000000000000000000000000000000
--- a/benchmarks/dis/analysis/rec_dis_electrons.cxx
+++ /dev/null
@@ -1,113 +0,0 @@
-#include "ROOT/RDataFrame.hxx"
-#include <iostream>
-
-#include "dd4pod/Geant4ParticleCollection.h"
-#include "eicd/TrackParametersCollection.h"
-#include "eicd/ClusterCollection.h"
-#include "eicd/ClusterData.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].psx * in[i].psx + in[i].psy * in[i].psy));
- }
- 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].psx, in[i].psy, in[i].psz, 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;
-};
-
-int rec_dis_electrons(const char* fname = "topside/rec_central_electrons.root")
-{
-
- ROOT::EnableImplicitMT();
- ROOT::RDataFrame df("events", fname);
-
- auto df0 = df.Define("isThrown", "mcparticles2.genStatus == 1")
- .Define("thrownParticles", "mcparticles2[isThrown]")
- .Define("thrownP", fourvec, {"thrownParticles"})
- .Define("p_thrown", momentum, {"thrownP"})
- .Define("nTracks", "outputTrackParameters.size()")
- .Define("p_track", p_track, {"outputTrackParameters"})
- .Define("p_track1", p_track, {"outputTrackParameters1"})
- .Define("p_track2", p_track, {"outputTrackParameters2"})
- .Define("delta_p",delta_p, {"p_track", "p_thrown"})
- .Define("delta_p1",delta_p, {"p_track1", "p_thrown"})
- .Define("delta_p2",delta_p, {"p_track2", "p_thrown"});
-
- auto h_nTracks = df0.Histo1D({"h_nTracks", "; N tracks ", 10, 0, 10}, "nTracks");
- auto h_pTracks = df0.Histo1D({"h_pTracks", "; GeV/c ", 100, 0, 10}, "p_track");
-
- auto h_delta_p = df0.Histo1D({"h_delta_p", "; GeV/c ", 100, -10, 10}, "delta_p");
- auto h_delta_p1 = df0.Histo1D({"h_delta_p1", "; GeV/c ", 100, -10, 10}, "delta_p1");
- auto h_delta_p2 = df0.Histo1D({"h_delta_p2", "; GeV/c ", 100, -10, 10}, "delta_p2");
-
- auto c = new TCanvas();
-
- h_nTracks->DrawCopy();
- c->SaveAs("results/dis/rec_central_electrons_nTracks.png");
- c->SaveAs("results/dis/rec_central_electrons_nTracks.pdf");
-
- h_pTracks->DrawCopy();
- c->SaveAs("results/dis/rec_central_electrons_pTracks.png");
- c->SaveAs("results/dis/rec_central_electrons_pTracks.pdf");
-
- THStack * hs = new THStack("hs_delta_p","; GeV/c ");
- TH1D* h1 = (TH1D*) h_delta_p->Clone();
- hs->Add(h1);
- h1 = (TH1D*) h_delta_p1->Clone();
- h1->SetLineColor(2);
- hs->Add(h1);
- h1 = (TH1D*) h_delta_p2->Clone();
- h1->SetLineColor(4);
- h1->SetFillStyle(3001);
- h1->SetFillColor(4);
- hs->Add(h1);
- hs->Draw("nostack");
- c->SaveAs("results/dis/rec_central_electrons_delta_p.png");
- c->SaveAs("results/dis/rec_central_electrons_delta_p.pdf");
-
- return 0;
-}
diff --git a/benchmarks/dis/dis.sh b/benchmarks/dis/dis.sh
old mode 100644
new mode 100755
index a863c7c5c66bcbccc6c97c4c99fd615e2bb95e1a..77fb85a577d6726e413166a153f2893b81aab1ef
--- a/benchmarks/dis/dis.sh
+++ b/benchmarks/dis/dis.sh
@@ -1,51 +1,61 @@
#!/bin/bash
## =============================================================================
-## Run the DVMP benchmarks in 5 steps:
-## 1. Build/install detector package
+## Run the DIS benchmarks in 5 steps:
+## 1. Parse the command line and setup environment
## 2. Detector simulation through npsim
## 3. Digitization and reconstruction through Juggler
## 4. Root-based Physics analyses
## 5. Finalize
## =============================================================================
-echo "Running the DIS benchmarks"
-
## make sure we launch this script from the project root directory
-PROJECT_ROOT="$( cd "$(dirname "$0")" >/dev/null 2>&1 ; pwd -P )"/..
+PROJECT_ROOT="$( cd "$(dirname "$0")" >/dev/null 2>&1 ; pwd -P )"/../..
pushd ${PROJECT_ROOT}
+echo "Running the DIS benchmarks"
+
## =============================================================================
-## Load the environment variables. To build the detector we need the following
-## variables:
+## Step 1: Setup the environment variables
##
+## First parse the command line flags.
+## This sets the following environment variables:
+## - CONFIG: The specific generator configuration
+## - EBEAM: The electron beam energy
+## - PBEAM: The ion beam energy
+source util/parse_cmd.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
## - DETECTOR_PATH: full path to the detector definitions
##
-## You can ready config/env.sh for more in-depth explanations of the variables
+## You can ready options/env.sh for more in-depth explanations of the variables
## and how they can be controlled.
-source config/env.sh
-
-## Extra environment variables for DVMP:
-## file tag for these tests
-JUGGLER_FILE_NAME_TAG="dis"
-# TODO use the input file name, as we will be generating a lot of these
-# in the future...
-# FIXME Generator file hardcoded for now
-## note: these variables need to be exported to be accessible from
-## the juggler options.py. We should really work on a dedicated
-## juggler launcher to get rid of these "magic" variables. FIXME
-export JUGGLER_GEN_FILE="results/dis/${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"
+source options/env.sh
+
+## We also need the following benchmark-specific variables:
+##
+## - BENCHMARK_TAG: Unique identified for this benchmark process.
+## - BEAM_TAG: Identifier for the chosen beam configuration
+## - INPUT_PATH: Path for generator-level input to the benchmarks
+## - TMP_PATH: Path for temporary data (not exported as artifacts)
+## - RESULTS_PATH: Path for benchmark output figures and files
+##
+## You can read dvmp/env.sh for more in-depth explanations of the variables.
+source benchmarks/dis/env.sh
+## Get a unique file names based on the configuration options
+GEN_FILE=${INPUT_PATH}/gen-${CONFIG}_${JUGGLER_N_EVENTS}.hepmc
+
+SIM_FILE=${TMP_PATH}/sim-${CONFIG}.root
+SIM_LOG=${TMP_PATH}/sim-${CONFIG}.log
-## =============================================================================
-## Step 1: Build/install the desired detector package
-## TODO remove this
-#bash util/build_detector.sh
+REC_FILE=${TMP_PATH}/rec-${CONFIG}.root
+REC_LOG=${TMP_PATH}/sim-${CONFIG}.log
+
+PLOT_TAG=${CONFIG}
## =============================================================================
## Step 2: Run the simulation
@@ -55,8 +65,8 @@ npsim --runType batch \
-v WARNING \
--numberOfEvents ${JUGGLER_N_EVENTS} \
--compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
- --inputFiles ${JUGGLER_GEN_FILE} \
- --outputFile ${JUGGLER_SIM_FILE}
+ --inputFiles ${GEN_FILE} \
+ --outputFile ${SIM_FILE}
if [ "$?" -ne "0" ] ; then
echo "ERROR running npsim"
exit 1
@@ -65,36 +75,67 @@ fi
## =============================================================================
## Step 3: Run digitization & reconstruction
echo "Running the digitization and reconstruction"
-# FIXME Need to figure out how to pass file name to juggler from the commandline
+## FIXME Need to figure out how to pass file name to juggler from the commandline
+## the tracker_reconstruction.py options file uses the following environment
+## variables:
+## - JUGGLER_SIM_FILE: input detector simulation
+## - JUGGLER_REC_FILE: output reconstructed data
+## - JUGGLER_DETECTOR_PATH: Location of the detector geometry
+## - JUGGLER_N_EVENTS: number of events to process (part of global environment)
+## - JUGGLER_DETECTOR: detector package (part of global environment)
+export JUGGLER_SIM_FILE=${SIM_FILE}
+export JUGGLER_REC_FILE=${REC_FILE}
+export JUGGLER_DETECTOR_PATH=${DETECTOR_PATH}
xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
- gaudirun.py config/tracker_reconstruction.py
+ gaudirun.py options/tracker_reconstruction.py \
+ 2>&1 > ${REC_LOG}
+## on-error, first retry running juggler again as there is still a random
+## crash we need to address FIXME
if [ "$?" -ne "0" ] ; then
- echo "ERROR running juggler"
- exit 1
+ echo "Juggler crashed, retrying..."
+ xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
+ gaudirun.py options/tracker_reconstruction.py \
+ 2>&1 > ${REC_LOG}
+ if [ "$?" -ne "0" ] ; then
+ echo "ERROR running juggler, both attempts failed"
+ exit 1
+ fi
fi
-ls -l
## =============================================================================
## Step 4: Analysis
-root -b -q "benchmarks/dis/analysis/rec_dis_electrons.cxx(\"${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}\")"
+## write a temporary configuration file for the analysis script
+echo "Running analysis"
+CONFIG="${TMP_PATH}/${PLOT_TAG}.json"
+cat << EOF > ${CONFIG}
+{
+ "rec_file": "${REC_FILE}",
+ "detector": "${JUGGLER_DETECTOR}",
+ "output_prefix": "${RESULTS_PATH}/${PLOT_TAG}",
+ "test_tag": "${BEAM_TAG}"
+}
+EOF
+#cat ${CONFIG}
+root -b -q "benchmarks/dis/analysis/dis_electrons.cxx+(\"${CONFIG}\")"
+#root -b -q "benchmarks/dis/analysis/dis_electrons.cxx(\"${CONFIG}\")"
if [[ "$?" -ne "0" ]] ; then
- echo "ERROR running root script"
+ echo "ERROR running rec_dis_electron script"
exit 1
fi
## =============================================================================
## Step 5: finalize
-echo "Finalizing ${JUGGLER_FILE_NAME_TAG} benchmark"
+echo "Finalizing DIS benchmark"
-## Copy over reconsturction artifacts as long as we don't have
+## Move over reconsturction artifacts as long as we don't have
## too many events
if [ "${JUGGLER_N_EVENTS}" -lt "500" ] ; then
- cp ${JUGGLER_REC_FILE} results/dis/.
+ cp ${REC_FILE} ${RESULTS_PATH}
fi
-## cleanup output files
-rm ${JUGGLER_REC_FILE} ${JUGGLER_SIM_FILE}
+## Always move over log files to the results path
+cp ${REC_LOG} ${RESULTS_PATH}
## =============================================================================
## All done!
-echo "${JUGGLER_FILE_NAME_TAG} benchmarks complete"
+echo "DIS benchmarks complete"
diff --git a/benchmarks/dis/env.sh b/benchmarks/dis/env.sh
new file mode 100644
index 0000000000000000000000000000000000000000..1a5b153f74f55b5a764019826ad061c7a6d196ba
--- /dev/null
+++ b/benchmarks/dis/env.sh
@@ -0,0 +1,50 @@
+#!/bin/bash
+
+## =============================================================================
+## Local configuration variables for this particular benchmark
+## It defines the following additional variables:
+##
+## - BENCHMARK_TAG: Tag to identify this particular benchmark
+## - BEAM_TAG Tag to identify the beam configuration
+## - INPUT_PATH: Path for generator-level input to the benchmarks
+## - TMP_PATH: Path for temporary data (not exported as artifacts)
+## - RESULTS_PATH: Path for benchmark output figures and files
+##
+## This script assumes that EBEAM and PBEAM are set as part of the
+## calling script (usually as command line argument).
+##
+## =============================================================================
+
+## Tag for the local benchmark. Should be the same as the directory name for
+## this particular benchmark set (for clarity).
+## This tag is used for the output artifacts directory (results/${JUGGLER_TAG})
+## and a tag in some of the output files.
+export BENCHMARK_TAG="dis"
+echo "Setting up the local environment for the ${BENCHMARK_TAG^^} benchmarks"
+
+## Extra beam tag to identify the desired beam configuration
+export BEAM_TAG="${EBEAM}on${PBEAM}"
+
+## Data path for input data (generator-level hepmc file)
+INPUT_PATH="input/${BENCHMARK_TAG}/${BEAM_TAG}"
+mkdir -p ${INPUT_PATH}
+export INPUT_PATH=`realpath ${INPUT_PATH}`
+echo "INPUT_PATH: ${INPUT_PATH}"
+
+
+## Data path for temporary data (not exported as artifacts)
+TMP_PATH=${LOCAL_PREFIX}/tmp/${BEAM_TAG}
+mkdir -p ${TMP_PATH}
+export TMP_PATH=`realpath ${TMP_PATH}`
+echo "TMP_PATH: ${TMP_PATH}"
+
+## Data path for benchmark output (plots and reconstructed files
+## if not too big).
+RESULTS_PATH="results/${BENCHMARK_TAG}/${BEAM_TAG}"
+mkdir -p ${RESULTS_PATH}
+export RESULTS_PATH=`realpath ${RESULTS_PATH}`
+echo "RESULTS_PATH: ${RESULTS_PATH}"
+
+## =============================================================================
+## That's all!
+echo "Local environment setup complete."
diff --git a/benchmarks/dis/gen.sh b/benchmarks/dis/gen.sh
old mode 100644
new mode 100755
index e09f2c4c7bf79f6e38e9ee09bdd0deacf3124177..4a42081ae6f437d542de677613eef22f4c6c9516
--- a/benchmarks/dis/gen.sh
+++ b/benchmarks/dis/gen.sh
@@ -3,38 +3,63 @@
## =============================================================================
## Standin for a proper pythia generation process, similar to how we
## generate events for DVMP
+## Runs in 5 steps:
+## 1. Parse the command line and setup the environment
+## 2. Check if we can load the requested file from the cache
+## 3. Build generator exe
+## 4. Run the actual generator
+## 5. Finalize
+## =============================================================================
## =============================================================================
-
-## TODO: use JUGGLER_FILE_NAME_TAG instead of explicitly refering to dis
-
-echo "Running the DIS benchmarks"
## make sure we launch this script from the project root directory
PROJECT_ROOT="$( cd "$(dirname "$0")" >/dev/null 2>&1 ; pwd -P )"/../..
pushd ${PROJECT_ROOT}
## =============================================================================
-## Load the environment variables. To build the detector we need the following
-## variables:
+## Step 1: Setup the environment variables
+## First parse the command line flags.
+## This sets the following environment variables:
+## - CONFIG: The specific generator configuration --> not currenlty used FIXME
+## - EBEAM: The electron beam energy --> not currently used FIXME
+## - PBEAM: The ion beam energy --> not currently used FIXME
+source util/parse_cmd.sh $@
+
+## To run the generator, 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
-## - DETECTOR_PATH: full path to the detector definitions
+## - LOCAL_PREFIX: Place to cache local packages and data
+## - JUGGLER_N_EVENTS: Number of events to process
+## - JUGGLER_RNG_SEED: Random seed for event generation.
##
-## You can ready config/env.sh for more in-depth explanations of the variables
+## You can read options/env.sh for more in-depth explanations of the variables
## and how they can be controlled.
-source config/env.sh
+source options/env.sh
+
+## We also need the following benchmark-specific variables:
+##
+## - BENCHMARK_TAG: Unique identified for this benchmark process.
+## - INPUT_PATH: Path for generator-level input to the benchmarks
+## - TMP_PATH: Path for temporary data (not exported as artifacts)
+##
+## You can read dvmp/env.sh for more in-depth explanations of the variables.
+source benchmarks/dis/env.sh
-## Setup local environment
-export DATA_PATH=results/dis
+## Get a unique file name prefix based on the configuration options
+GEN_TAG=gen-${CONFIG}_${JUGGLER_N_EVENTS} ## Generic file prefix
-## Extra environment variables for DVMP:
-## file tag for these tests
-JUGGLER_FILE_NAME_TAG="dis"
+## =============================================================================
+## Step 2: Check if we really need to run, or can use the cache.
+if [ -f "${INPUT_PATH}/${GEN_TAG}.hepmc" ]; then
+ echo "Found cached generator output for $GEN_TAG, no need to rerun"
+ exit 0
+fi
+
+echo "Generator output for $GEN_TAG not found in cache, need to run generator"
## =============================================================================
-## Step 1: Dummy event generator
-## TODO better file name that encodes the actual configuration we're running
+## Step 3: Build generator exe
+## TODO: need to configurability to the generator exe
+
echo "Compiling benchmarks/dis/generator/pythia_dis.cxx ..."
g++ benchmarks/dis/generator/pythia_dis.cxx -o pythia_dis \
-I/usr/local/include -Iinclude \
@@ -46,20 +71,26 @@ if [[ "$?" -ne "0" ]] ; then
exit 1
fi
echo "done"
-pwd
-ls -lrth
-./pythia_dis dis.hepmc
+## =============================================================================
+## Step 4: Run the event generator
+echo "Running the generator"
+./pythia_dis ${TMP_PATH}/${GEN_TAG}.hepmc
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running pythia"
exit 1
fi
+
+## =============================================================================
+## Step 5: Finally, move relevant output into the artifacts directory and clean up
## =============================================================================
-## Step 2: finalize
-echo "Moving event generator output into ${DATA_PATH}"
-#mv .local/${JUGGLER_FILE_NAME_TAG}.hepmc ${DATA_PATH}/${JUGGLER_FILE_NAME_TAG}.hepmc
+echo "Moving generator output into ${INPUT_PATH}"
+mv ${TMP_PATH}/${GEN_TAG}.hepmc ${INPUT_PATH}/${GEN_TAG}.hepmc
+## this step only matters for local execution
+echo "Cleaning up"
+## does nothing
## =============================================================================
## All done!
-echo "dis event generation complete"
+echo "$BENCHMARK_TAG event generation complete"
diff --git a/benchmarks/dis/generator/pythia_dis.cxx b/benchmarks/dis/generator/pythia_dis.cxx
index ff819a7e3243373539aace192d0990003ca3368d..a48679683baf6ae692825275e713f79b8317f685 100644
--- a/benchmarks/dis/generator/pythia_dis.cxx
+++ b/benchmarks/dis/generator/pythia_dis.cxx
@@ -1,7 +1,7 @@
#include "Pythia8/Pythia.h"
#include "Pythia8Plugins/HepMC3.h"
-#include <unistd.h>
#include <cassert>
+#include <unistd.h>
using namespace Pythia8;
@@ -13,7 +13,7 @@ using std::string;
enum class mode { none, help, list, part };
struct settings {
- double E_electron = 10.0; // GeV
+ double E_electron = 18.0; // GeV
double E_ion = 275.0; // GeV
std::string outfile = "dis.hepmc";
int ion_PID = 2212;
@@ -22,7 +22,7 @@ struct settings {
bool success = false;
double Q2_min = 4.0;
int N_events = 1000;
- mode selected = mode::none;
+ mode selected = mode::none;
};
template <typename T>
@@ -63,16 +63,15 @@ void print_man_page(T cli, const char* argv0)
}
//______________________________________________________________________________
-settings cmdline_settings(int argc, char* argv[]) {
+settings cmdline_settings(int argc, char* argv[])
+{
settings s;
- auto lastOpt =
- " options:" % (option("-h", "--help").set(s.selected, mode::help) % "show help",
- value("file", s.outfile).if_missing([] {
- std::cout << "You need to provide an output filename argument!\n";
- }) % "output file");
+ auto lastOpt = " options:" % (option("-h", "--help").set(s.selected, mode::help) % "show help",
+ value("file", s.outfile).if_missing([] {
+ std::cout << "You need to provide an output filename argument!\n";
+ }) % "output file");
- auto cli =
- (command("help").set(s.selected, mode::help) | lastOpt);
+ auto cli = (command("help").set(s.selected, mode::help) | lastOpt);
assert(cli.flags_are_prefix_free());
@@ -95,7 +94,8 @@ settings cmdline_settings(int argc, char* argv[]) {
}
//______________________________________________________________________________
-int main(int argc, char* argv[]) {
+int main(int argc, char* argv[])
+{
settings s = cmdline_settings(argc, argv);
if (!s.success) {
diff --git a/benchmarks/dvcs/dvcs.sh b/benchmarks/dvcs/dvcs.sh
index 39143562a090fdde535a899ae5bf7e1420f7a0a5..35bbcafeddfe3e2e218ba286f209464cf8a297d7 100644
--- a/benchmarks/dvcs/dvcs.sh
+++ b/benchmarks/dvcs/dvcs.sh
@@ -15,9 +15,9 @@ echo "JUGGLER_FILE_NAME_TAG = ${JUGGLER_FILE_NAME_TAG}"
## - JUGGLER_DETECTOR: the detector package we want to use for this benchmark
## - DETECTOR_PATH: full path to the detector definitions
##
-## You can ready config/env.sh for more in-depth explanations of the variables
+## You can ready options/env.sh for more in-depth explanations of the variables
## and how they can be controlled.
-source config/env.sh
+source options/env.sh
curl -o test_proton_dvcs_eic.hepmc "https://eicweb.phy.anl.gov/api/v4/projects/345/jobs/artifacts/master/raw/data/test_proton_dvcs_eic.hepmc?job=compile"
@@ -43,7 +43,7 @@ fi
# Need to figure out how to pass file name to juggler from the commandline
xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
- gaudirun.py config/tracker_reconstruction.py
+ gaudirun.py options/tracker_reconstruction.py
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running juggler"
exit 1
diff --git a/benchmarks/dvmp/analysis/dvmp.h b/benchmarks/dvmp/analysis/dvmp.h
new file mode 100644
index 0000000000000000000000000000000000000000..371b454f8d92f4c0cd5b17edd4a8f5e1ecc77e77
--- /dev/null
+++ b/benchmarks/dvmp/analysis/dvmp.h
@@ -0,0 +1,56 @@
+#ifndef DVMP_H
+#define DVMP_H
+
+#include <util.h>
+
+#include <algorithm>
+#include <cmath>
+#include <exception>
+#include <fmt/core.h>
+#include <limits>
+#include <string>
+#include <vector>
+
+#include <Math/Vector4D.h>
+
+// Additional utility functions for DVMP benchmarks. Where useful, these can be
+// promoted to the top-level util library
+namespace util {
+
+ // Calculate the 4-vector sum of a given pair of particles
+ inline ROOT::Math::PxPyPzMVector
+ get_sum(const std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>& particle_pair)
+ {
+ return (particle_pair.first + particle_pair.second);
+ }
+
+ //========================================================================================================
+ // for structure functions
+
+ struct inv_quant { // add more when needed
+ double nu, Q2, x;
+ };
+
+ // for simu
+ inline inv_quant calc_inv_quant_simu(const std::vector<ROOT::Math::PxPyPzMVector>& parts)
+ {
+ ROOT::Math::PxPyPzMVector q(parts[0] - parts[2]);
+ ROOT::Math::PxPyPzMVector P(parts[3]);
+
+ double nu = q.Dot(P) / P.mass();
+ double Q2 = -q.Dot(q);
+ inv_quant quantities = {nu, Q2, Q2 / 2. / P.mass() / nu};
+ return quantities;
+ }
+
+ inline double get_nu_simu(inv_quant quantities) { return quantities.nu / 1000.; }
+ inline double get_Q2_simu(inv_quant quantities) { return quantities.Q2; }
+ inline double get_x_simu(inv_quant quantities) { return quantities.x; }
+
+ // for tracking, add later
+
+ //=========================================================================================================
+
+} // namespace util
+
+#endif
diff --git a/benchmarks/dvmp/analysis/plot.h b/benchmarks/dvmp/analysis/plot.h
deleted file mode 100644
index 69ebba341af239541496a410c66f78eb33e8adac..0000000000000000000000000000000000000000
--- a/benchmarks/dvmp/analysis/plot.h
+++ /dev/null
@@ -1,40 +0,0 @@
-#ifndef PLOT_H
-#define PLOT_H
-
-#include <TColor.h>
-#include <fmt/core.h>
-#include <vector>
-
-namespace plot {
-
-const int kMpBlue = TColor::GetColor(0x1f, 0x77, 0xb4);
-const int kMpOrange = TColor::GetColor(0xff, 0x7f, 0x0e);
-const int kMpGreen = TColor::GetColor(0x2c, 0xa0, 0x2c);
-const int kMpRed = TColor::GetColor(0xd6, 0x27, 0x28);
-const int kMpPurple = TColor::GetColor(0x94, 0x67, 0xbd);
-const int kMpBrown = TColor::GetColor(0x8c, 0x56, 0x4b);
-const int kMpPink = TColor::GetColor(0xe3, 0x77, 0xc2);
-const int kMpGrey = TColor::GetColor(0x7f, 0x7f, 0x7f);
-const int kMpMoss = TColor::GetColor(0xbc, 0xbd, 0x22);
-const int kMpCyan = TColor::GetColor(0x17, 0xbe, 0xcf);
-
-const std::vector<int> kPalette = {kMpBlue, kMpOrange, kMpGreen, kMpRed,
- kMpPurple, kMpBrown, kMpPink, kMpGrey,
- kMpMoss, kMpCyan};
-
-void draw_label(int ebeam, int pbeam, const std::string_view detector,
- std::string_view vm, std::string_view what) {
- auto t = new TPaveText(.15, 0.800, .7, .925, "NB NDC");
- t->SetFillColorAlpha(kWhite, 0.4);
- t->SetTextFont(43);
- t->SetTextSize(25);
- t->AddText(fmt::format("#bf{{{} }}SIMULATION", detector).c_str());
- t->AddText(fmt::format("{} for {} DVMP.", what, vm).c_str());
- t->AddText(fmt::format("{} GeV on {} GeV", ebeam, pbeam, what).c_str());
- t->SetTextAlign(12);
- t->Draw();
-}
-
-} // namespace plot
-
-#endif
diff --git a/benchmarks/dvmp/analysis/vm_invar.cxx b/benchmarks/dvmp/analysis/vm_invar.cxx
index d6d0804f73b0e18d2a90c2ec478d26168d214f21..03349825902a7e1210251fab0b4d8b467ddfae9d 100644
--- a/benchmarks/dvmp/analysis/vm_invar.cxx
+++ b/benchmarks/dvmp/analysis/vm_invar.cxx
@@ -1,3 +1,4 @@
+#include "dvmp.h"
#include "plot.h"
#include <benchmark.h>
@@ -18,8 +19,6 @@
// a desired vector meson (e.g. jpsi) and a desired decay particle (e.g. muon)
// Output figures are written to our output prefix (which includes the output
// file prefix), and labeled with our detector name.
-// TODO: I think it would be better to pass small json configuration file to
-// the test, instead of this ever-expanding list of function arguments.
// FIXME: MC does not trace back into particle history. Need to fix that
int vm_invar(const std::string& config_name)
{
@@ -35,17 +34,20 @@ int vm_invar(const std::string& config_name)
std::string output_prefix = config["output_prefix"];
const std::string test_tag = config["test_tag"];
- fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green), "Running VM invariant mass analysis...\n");
+ fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green),
+ "Running VM invariant mass analysis...\n");
fmt::print(" - Vector meson: {}\n", vm_name);
fmt::print(" - Decay particle: {}\n", decay_name);
fmt::print(" - Detector package: {}\n", detector);
+ fmt::print(" - input file: {}\n", rec_file);
fmt::print(" - output prefix: {}\n", output_prefix);
// create our test definition
// test_tag
- eic::util::Test vm_mass_resolution_test{
- {{"name", fmt::format("{}_{}_{}_mass_resolution", test_tag, vm_name, decay_name)},
- {"title", fmt::format("{} -> {} Invariant Mass Resolution", vm_name, decay_name)},
+ eic::util::Test Q2_resolution_test{
+ {{"name", fmt::format("{}_{}_{}_Q2_resolution", test_tag, vm_name, decay_name)},
+ {"title",
+ fmt::format("Q^2 Resolution for {} -> {} events with {}", vm_name, decay_name, detector)},
{"description", "Invariant Mass Resolution calculated from raw "
"tracking data using a Gaussian fit."},
{"quantity", "resolution"},
@@ -114,7 +116,7 @@ int vm_invar(const std::string& config_name)
// draw everything
hnu.DrawClone("hist");
// FIXME hardcoded beam configuration
- plot::draw_label(10, 100, detector, vm_name, "#nu");
+ plot::draw_label(10, 100, detector);
TText* tptr21;
auto t21 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t21->SetFillColorAlpha(kWhite, 0);
@@ -139,7 +141,7 @@ int vm_invar(const std::string& config_name)
// draw everything
hQ2.DrawClone("hist");
// FIXME hardcoded beam configuration
- plot::draw_label(10, 100, detector, vm_name, "Q^{2}");
+ plot::draw_label(10, 100, detector);
TText* tptr22;
auto t22 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t22->SetFillColorAlpha(kWhite, 0);
@@ -164,7 +166,7 @@ int vm_invar(const std::string& config_name)
// draw everything
hx.DrawClone("hist");
// FIXME hardcoded beam configuration
- plot::draw_label(10, 100, detector, vm_name, "x");
+ plot::draw_label(10, 100, detector);
TText* tptr23;
auto t23 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t23->SetFillColorAlpha(kWhite, 0);
@@ -179,12 +181,12 @@ int vm_invar(const std::string& config_name)
c.Print(fmt::format("{}InvariantQuantities.png", output_prefix).c_str());
}
- // TODO we're not actually doing an IM fit yet, so for now just return an
+ // TODO we're not actually getting the resolutions yet
// error for the test result
- vm_mass_resolution_test.error(-1);
+ Q2_resolution_test.error(-1);
// write out our test data
- eic::util::write_test(vm_mass_resolution_test, fmt::format("{}vm_invar.json", output_prefix));
+ eic::util::write_test(Q2_resolution_test, fmt::format("{}vm_invar.json", output_prefix));
// That's all!
return 0;
diff --git a/benchmarks/dvmp/analysis/vm_mass.cxx b/benchmarks/dvmp/analysis/vm_mass.cxx
index 6420e1fe1b55d2304d0817dfe6b4b51a7623c1c2..82abba6f77d8343df259e89d823e0b7d14419534 100644
--- a/benchmarks/dvmp/analysis/vm_mass.cxx
+++ b/benchmarks/dvmp/analysis/vm_mass.cxx
@@ -1,3 +1,4 @@
+#include "dvmp.h"
#include "plot.h"
#include <benchmark.h>
@@ -35,17 +36,20 @@ int vm_mass(const std::string& config_name)
std::string output_prefix = config["output_prefix"];
const std::string test_tag = config["test_tag"];
- fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green), "Running VM invariant mass analysis...\n");
+ fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green),
+ "Running VM invariant mass analysis...\n");
fmt::print(" - Vector meson: {}\n", vm_name);
fmt::print(" - Decay particle: {}\n", decay_name);
fmt::print(" - Detector package: {}\n", detector);
+ fmt::print(" - input file: {}\n", rec_file);
fmt::print(" - output prefix: {}\n", output_prefix);
// create our test definition
// test_tag
- eic::util::Test vm_mass_resolution_test{
+ eic::util::Test mass_resolution_test{
{{"name", fmt::format("{}_{}_{}_mass_resolution", test_tag, vm_name, decay_name)},
- {"title", fmt::format("{} -> {} Invariant Mass Resolution", vm_name, decay_name)},
+ {"title", fmt::format("{} Invariant Mass Resolution for {} -> {} with {}", vm_name, vm_name,
+ decay_name, detector)},
{"description", "Invariant Mass Resolution calculated from raw "
"tracking data using a Gaussian fit."},
{"quantity", "resolution"},
@@ -82,18 +86,23 @@ int vm_mass(const std::string& config_name)
.Define("p_sim", util::momenta_from_simulation, {"mcparticles2"})
.Define("decay_pair_rec", find_decay_pair, {"p_rec"})
.Define("decay_pair_sim", find_decay_pair, {"p_sim"})
- .Define("mass_rec", util::get_im, {"decay_pair_rec"})
- .Define("mass_sim", util::get_im, {"decay_pair_sim"})
- .Define("pt_rec", util::get_pt, {"decay_pair_rec"})
- .Define("pt_sim", util::get_pt, {"decay_pair_sim"})
- .Define("phi_rec", util::get_phi, {"decay_pair_rec"})
- .Define("phi_sim", util::get_phi, {"decay_pair_sim"})
- .Define("rapidity_rec", util::get_y, {"decay_pair_rec"})
- .Define("rapidity_sim", util::get_y, {"decay_pair_sim"});
+ .Define("p_vm_rec", "decay_pair_rec.first + decay_pair_rec.second")
+ .Define("p_vm_sim", "decay_pair_sim.first + decay_pair_sim.second")
+ //.Define("p_vm_sim", util::get_sum, {"decay_pair_sim"})
+ .Define("mass_rec", "p_vm_rec.M()")
+ .Define("mass_sim", "p_vm_sim.M()")
+ .Define("pt_rec", "p_vm_rec.pt()")
+ .Define("pt_sim", "p_vm_sim.pt()")
+ .Define("phi_rec", "p_vm_rec.phi()")
+ .Define("phi_sim", "p_vm_sim.phi()")
+ .Define("eta_rec", "p_vm_rec.eta()")
+ .Define("eta_sim", "p_vm_sim.eta()");
// Define output histograms
- auto h_im_rec = d_im.Histo1D({"h_im_rec", ";m_{ll'} (GeV/c^{2});#", 100, -1.1, vm_mass + 5}, "mass_rec");
- auto h_im_sim = d_im.Histo1D({"h_im_sim", ";m_{ll'} (GeV/c^{2});#", 100, -1.1, vm_mass + 5}, "mass_sim");
+ auto h_im_rec =
+ d_im.Histo1D({"h_im_rec", ";m_{ll'} (GeV/c^{2});#", 100, -1.1, vm_mass + 5}, "mass_rec");
+ auto h_im_sim =
+ d_im.Histo1D({"h_im_sim", ";m_{ll'} (GeV/c^{2});#", 100, -1.1, vm_mass + 5}, "mass_sim");
auto h_pt_rec = d_im.Histo1D({"h_pt_rec", ";p_{T} (GeV/c);#", 400, 0., 40.}, "pt_rec");
auto h_pt_sim = d_im.Histo1D({"h_pt_sim", ";p_{T} (GeV/c);#", 400, 0., 40.}, "pt_sim");
@@ -101,8 +110,8 @@ int vm_mass(const std::string& config_name)
auto h_phi_rec = d_im.Histo1D({"h_phi_rec", ";#phi_{ll'};#", 90, -M_PI, M_PI}, "phi_rec");
auto h_phi_sim = d_im.Histo1D({"h_phi_sim", ";#phi_{ll'};#", 90, -M_PI, M_PI}, "phi_sim");
- auto h_y_rec = d_im.Histo1D({"h_y_rec", ";y_{ll'};#", 1000, -5., 5.}, "rapidity_rec");
- auto h_y_sim = d_im.Histo1D({"h_y_sim", ";y_{ll'};#", 1000, -5., 5.}, "rapidity_sim");
+ auto h_eta_rec = d_im.Histo1D({"h_eta_rec", ";#eta_{ll'};#", 1000, -5., 5.}, "eta_rec");
+ auto h_eta_sim = d_im.Histo1D({"h_eta_sim", ";#eta_{ll'};#", 1000, -5., 5.}, "eta_sim");
// Plot our histograms.
// TODO: to start I'm explicitly plotting the histograms, but want to
@@ -128,7 +137,7 @@ int vm_mass(const std::string& config_name)
h11.DrawClone("hist");
h12.DrawClone("hist same");
// FIXME hardcoded beam configuration
- plot::draw_label(10, 100, detector, vm_name, "Invariant mass");
+ plot::draw_label(10, 100, detector);
TText* tptr1;
auto t1 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t1->SetFillColorAlpha(kWhite, 0);
@@ -158,7 +167,7 @@ int vm_mass(const std::string& config_name)
h21.DrawClone("hist");
h22.DrawClone("hist same");
// FIXME hardcoded beam configuration
- plot::draw_label(10, 100, detector, vm_name, "Transverse Momentum");
+ plot::draw_label(10, 100, detector);
TText* tptr2;
auto t2 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t2->SetFillColorAlpha(kWhite, 0);
@@ -188,7 +197,7 @@ int vm_mass(const std::string& config_name)
h31.DrawClone("hist");
h32.DrawClone("hist same");
// FIXME hardcoded beam configuration
- plot::draw_label(10, 100, detector, vm_name, "#phi");
+ plot::draw_label(10, 100, detector);
TText* tptr3;
auto t3 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t3->SetFillColorAlpha(kWhite, 0);
@@ -204,8 +213,8 @@ int vm_mass(const std::string& config_name)
c.cd(4);
// gPad->SetLogx(false);
// gPad->SetLogy(false);
- auto& h41 = *h_y_sim;
- auto& h42 = *h_y_rec;
+ auto& h41 = *h_eta_sim;
+ auto& h42 = *h_eta_rec;
// histogram style
h41.SetLineColor(plot::kMpBlue);
h41.SetLineWidth(2);
@@ -218,7 +227,7 @@ int vm_mass(const std::string& config_name)
h41.DrawClone("hist");
h42.DrawClone("hist same");
// FIXME hardcoded beam configuration
- plot::draw_label(10, 100, detector, vm_name, "Rapidity");
+ plot::draw_label(10, 100, detector);
TText* tptr4;
auto t4 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t4->SetFillColorAlpha(kWhite, 0);
@@ -235,10 +244,10 @@ int vm_mass(const std::string& config_name)
// TODO we're not actually doing an IM fit yet, so for now just return an
// error for the test result
- vm_mass_resolution_test.error(-1);
+ mass_resolution_test.error(-1);
// write out our test data
- eic::util::write_test(vm_mass_resolution_test, fmt::format("{}vm_mass.json", output_prefix));
+ eic::util::write_test(mass_resolution_test, fmt::format("{}_mass.json", output_prefix));
// That's all!
return 0;
diff --git a/benchmarks/dvmp/dvmp.sh b/benchmarks/dvmp/dvmp.sh
index 87a0cef4e895a75fe73f111e2da92726fb12dfde..1746a8600c1a3e20a8968b43db51884d2311e2fc 100755
--- a/benchmarks/dvmp/dvmp.sh
+++ b/benchmarks/dvmp/dvmp.sh
@@ -35,9 +35,9 @@ source util/parse_cmd.sh $@
## - JUGGLER_DETECTOR: the detector package we want to use for this benchmark
## - DETECTOR_PATH: full path to the detector definitions
##
-## You can ready config/env.sh for more in-depth explanations of the variables
+## You can ready options/env.sh for more in-depth explanations of the variables
## and how they can be controlled.
-source config/env.sh
+source options/env.sh
## We also need the following benchmark-specific variables:
##
@@ -93,21 +93,20 @@ export JUGGLER_SIM_FILE=${SIM_FILE}
export JUGGLER_REC_FILE=${REC_FILE}
export JUGGLER_DETECTOR_PATH=${DETECTOR_PATH}
xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
- gaudirun.py config/tracker_reconstruction.py \
+ gaudirun.py options/tracker_reconstruction.py \
2>&1 > ${REC_LOG}
## on-error, first retry running juggler again as there is still a random
## crash we need to address FIXME
if [ "$?" -ne "0" ] ; then
echo "Juggler crashed, retrying..."
xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
- gaudirun.py config/tracker_reconstruction.py \
+ gaudirun.py options/tracker_reconstruction.py \
2>&1 > ${REC_LOG}
if [ "$?" -ne "0" ] ; then
echo "ERROR running juggler, both attempts failed"
exit 1
fi
fi
-#ls -l
## =============================================================================
## Step 4: Analysis
@@ -127,10 +126,14 @@ EOF
#cat ${CONFIG}
## run the analysis script with this configuration
-root -b -q "benchmarks/dvmp/analysis/vm_mass.cxx(\"${CONFIG}\")"
-root -b -q "benchmarks/dvmp/analysis/vm_invar.cxx(\"${CONFIG}\")"
+root -b -q "benchmarks/dvmp/analysis/vm_mass.cxx+(\"${CONFIG}\")"
if [ "$?" -ne "0" ] ; then
- echo "ERROR running root script"
+ echo "ERROR running vm_mass script"
+ exit 1
+fi
+root -b -q "benchmarks/dvmp/analysis/vm_invar.cxx+(\"${CONFIG}\")"
+if [ "$?" -ne "0" ] ; then
+ echo "ERROR running vm_invar script"
exit 1
fi
@@ -145,11 +148,11 @@ if [ "${JUGGLER_N_EVENTS}" -lt "500" ] ; then
fi
## Always move over log files to the results path
-mv ${SIM_LOG} ${REC_LOG} ${RESULTS_PATH}
+mv ${REC_LOG} ${RESULTS_PATH}
## cleanup output files
-rm -f ${REC_FILE} ${SIM_FILE}
+#rm -f ${REC_FILE} ${SIM_FILE} ## --> not needed for CI
## =============================================================================
## All done!
diff --git a/benchmarks/dvmp/gen.sh b/benchmarks/dvmp/gen.sh
index cfb64aa079b5a266b89b89a05cc082232e9a48f2..8bef0a116f14014869d134c38720e7717ad5e49d 100755
--- a/benchmarks/dvmp/gen.sh
+++ b/benchmarks/dvmp/gen.sh
@@ -33,9 +33,9 @@ source util/parse_cmd.sh $@
## - JUGGLER_N_EVENTS: Number of events to process
## - JUGGLER_RNG_SEED: Random seed for event generation.
##
-## You can read config/env.sh for more in-depth explanations of the variables
+## You can read options/env.sh for more in-depth explanations of the variables
## and how they can be controlled.
-source config/env.sh
+source options/env.sh
## We also need the following benchmark-specific variables:
##
@@ -107,4 +107,6 @@ done
echo "Cleaning up"
rm ${TMP_PATH}/${GEN_TAG}.json
+## =============================================================================
## All done!
+echo "$BENCHMARK_TAG event generation complete"
diff --git a/include/benchmark.h b/include/benchmark.h
index 48c9694dc5491bc32fa6b89bba8f672ec4cc3eae..0634be0a1494739c4e7337f99e35ffbac0be7f7f 100644
--- a/include/benchmark.h
+++ b/include/benchmark.h
@@ -4,6 +4,7 @@
#include "exception.h"
#include <fmt/core.h>
#include <fstream>
+#include <iomanip>
#include <iostream>
#include <nlohmann/json.hpp>
#include <string>
@@ -67,14 +68,19 @@ namespace eic::util {
// - weight: Weight for this test (this is defaulted to 1.0 if not specified)
// - result: pass/fail/error
struct Test {
- Test(nlohmann::json definition) : json{std::move(definition)}
+ // note: round braces for the json constructor, as else it will pick the wrong
+ // initializer-list constructur (it will put everything in an array)
+ Test(const std::map<std::string, std::string>& definition) : json(definition)
{
+ // std::cout << json.dump() << std::endl;
// initialize with error (as we don't have a value yet)
error();
// Check that all required fields are present
- for (const auto& field : {"name", "title", "description", "quantity", "target", "value", "result"}) {
+ for (const auto& field :
+ {"name", "title", "description", "quantity", "target", "value", "result"}) {
if (json.find(field) == json.end()) {
- throw TestDefinitionError{fmt::format("Error in test definition: field '{}' missing", field)};
+ throw TestDefinitionError{
+ fmt::format("Error in test definition: field '{}' missing", field)};
}
}
// Default "weight" to 1 if not set
diff --git a/include/plot.h b/include/plot.h
new file mode 100644
index 0000000000000000000000000000000000000000..c198616325d54d73ae5e4c08f5ba7113ae77f817
--- /dev/null
+++ b/include/plot.h
@@ -0,0 +1,42 @@
+#ifndef PLOT_H
+#define PLOT_H
+
+#include <TCanvas.h>
+#include <TColor.h>
+#include <TPad.h>
+#include <TPaveText.h>
+#include <TStyle.h>
+#include <fmt/core.h>
+#include <vector>
+
+namespace plot {
+
+ const int kMpBlue = TColor::GetColor(0x1f, 0x77, 0xb4);
+ const int kMpOrange = TColor::GetColor(0xff, 0x7f, 0x0e);
+ const int kMpGreen = TColor::GetColor(0x2c, 0xa0, 0x2c);
+ const int kMpRed = TColor::GetColor(0xd6, 0x27, 0x28);
+ const int kMpPurple = TColor::GetColor(0x94, 0x67, 0xbd);
+ const int kMpBrown = TColor::GetColor(0x8c, 0x56, 0x4b);
+ const int kMpPink = TColor::GetColor(0xe3, 0x77, 0xc2);
+ const int kMpGrey = TColor::GetColor(0x7f, 0x7f, 0x7f);
+ const int kMpMoss = TColor::GetColor(0xbc, 0xbd, 0x22);
+ const int kMpCyan = TColor::GetColor(0x17, 0xbe, 0xcf);
+
+ const std::vector<int> kPalette = {kMpBlue, kMpOrange, kMpGreen, kMpRed, kMpPurple,
+ kMpBrown, kMpPink, kMpGrey, kMpMoss, kMpCyan};
+
+ void draw_label(int ebeam, int pbeam, const std::string_view detector)
+ {
+ auto t = new TPaveText(.15, 0.800, .7, .925, "NB NDC");
+ t->SetFillColorAlpha(kWhite, 0.4);
+ t->SetTextFont(43);
+ t->SetTextSize(25);
+ t->AddText(fmt::format("#bf{{{} }}SIMULATION", detector).c_str());
+ t->AddText(fmt::format("{} GeV on {} GeV", ebeam, pbeam).c_str());
+ t->SetTextAlign(12);
+ t->Draw();
+ }
+
+} // namespace plot
+
+#endif
diff --git a/include/util.h b/include/util.h
index dad7048e016d7c5d03cca9eb37481472ea4e5a74..56fb12893787e1d1fff00bd01653ad23b149f48e 100644
--- a/include/util.h
+++ b/include/util.h
@@ -1,6 +1,8 @@
#ifndef UTIL_H
#define UTIL_H
+// TODO: should probably be moved to a global benchmark utility library
+
#include <algorithm>
#include <cmath>
#include <exception>
@@ -16,182 +18,141 @@
namespace util {
-// Exception definition for unknown particle errors
-// FIXME: A utility exception base class should be included in the analysis
-// utility library, so we can skip most of this boilerplate
-class unknown_particle_error : public std::exception {
-public:
- unknown_particle_error(std::string_view particle) : m_particle{particle} {}
- virtual const char* what() const throw() {
- return fmt::format("Unknown particle type: {}", m_particle).c_str();
- }
- virtual const char* type() const throw() { return "unknown_particle_error"; }
-
-private:
- const std::string m_particle;
-};
-
-// Simple function to return the appropriate PDG mass for the particles
-// we care about for this process.
-// FIXME: consider something more robust (maybe based on hepPDT) to the
-// analysis utility library
-inline double get_pdg_mass(std::string_view part) {
- if (part == "electron") {
- return 0.0005109989461;
- } else if (part == "muon") {
- return .1056583745;
- } else if (part == "jpsi") {
- return 3.0969;
- } else if (part == "upsilon") {
- return 9.49630;
- } else {
- throw unknown_particle_error{part};
- }
-}
-
-// Get a vector of 4-momenta from raw tracking info, using an externally
-// provided particle mass assumption.
-// FIXME: should be part of utility library
-inline auto
-momenta_from_tracking(const std::vector<eic::TrackParametersData>& tracks,
- const double mass) {
- std::vector<ROOT::Math::PxPyPzMVector> momenta{tracks.size()};
- // transform our raw tracker info into proper 4-momenta
- std::transform(tracks.begin(), tracks.end(), momenta.begin(),
- [mass](const auto& track) {
- // make sure we don't divide by zero
- if (fabs(track.qOverP) < 1e-9) {
- return ROOT::Math::PxPyPzMVector{};
- }
- const double p = fabs(1. / track.qOverP);
- const double px = p * cos(track.phi) * sin(track.theta);
- const double py = p * sin(track.phi) * sin(track.theta);
- const double pz = p * cos(track.theta);
- return ROOT::Math::PxPyPzMVector{px, py, pz, mass};
- });
- return momenta;
-}
-
-// Get a vector of 4-momenta from the simulation data.
-// FIXME: should be part of utility library
-// TODO: Add PID selector (maybe using ranges?)
-inline auto
-momenta_from_simulation(const std::vector<dd4pod::Geant4ParticleData>& parts) {
- std::vector<ROOT::Math::PxPyPzMVector> momenta{parts.size()};
- // transform our simulation particle data into 4-momenta
- std::transform(parts.begin(), parts.end(), momenta.begin(),
- [](const auto& part) {
- return ROOT::Math::PxPyPzMVector{part.psx, part.psy,
- part.psz, part.mass};
- });
- return momenta;
-}
-
-// Find the decay pair candidates from a vector of particles (parts),
-// with invariant mass closest to a desired value (pdg_mass)
-// FIXME: not sure if this belongs here, or in the utility library. Probably the
-// utility library
-inline std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>
-find_decay_pair(const std::vector<ROOT::Math::PxPyPzMVector>& parts,
- const double pdg_mass) {
- int first = -1;
- int second = -1;
- double best_mass = -1;
-
- // go through all particle combinatorics, calculate the invariant mass
- // for each combination, and remember which combination is the closest
- // to the desired pdg_mass
- for (int i = 0; i < parts.size(); ++i) {
- for (int j = i + 1; j < parts.size(); ++j) {
- const double new_mass{(parts[i] + parts[j]).mass()};
- if (fabs(new_mass - pdg_mass) < fabs(best_mass - pdg_mass)) {
- first = i;
- second = j;
- best_mass = new_mass;
- }
+ // Exception definition for unknown particle errors
+ // FIXME: A utility exception base class should be included in the analysis
+ // utility library, so we can skip most of this boilerplate
+ class unknown_particle_error : public std::exception {
+ public:
+ unknown_particle_error(std::string_view particle) : m_particle{particle} {}
+ virtual const char* what() const throw()
+ {
+ return fmt::format("Unknown particle type: {}", m_particle).c_str();
+ }
+ virtual const char* type() const throw() { return "unknown_particle_error"; }
+
+ private:
+ const std::string m_particle;
+ };
+
+ // Simple function to return the appropriate PDG mass for the particles
+ // we care about for this process.
+ // FIXME: consider something more robust (maybe based on hepPDT) to the
+ // analysis utility library
+ inline double get_pdg_mass(std::string_view part)
+ {
+ if (part == "electron") {
+ return 0.0005109989461;
+ } else if (part == "muon") {
+ return .1056583745;
+ } else if (part == "jpsi") {
+ return 3.0969;
+ } else if (part == "upsilon") {
+ return 9.49630;
+ } else {
+ throw unknown_particle_error{part};
}
}
- if (first < 0) {
- return {{}, {}};
+
+ // Get a vector of 4-momenta from raw tracking info, using an externally
+ // provided particle mass assumption.
+ inline auto momenta_from_tracking(const std::vector<eic::TrackParametersData>& tracks,
+ const double mass)
+ {
+ std::vector<ROOT::Math::PxPyPzMVector> momenta{tracks.size()};
+ // transform our raw tracker info into proper 4-momenta
+ std::transform(tracks.begin(), tracks.end(), momenta.begin(), [mass](const auto& track) {
+ // make sure we don't divide by zero
+ if (fabs(track.qOverP) < 1e-9) {
+ return ROOT::Math::PxPyPzMVector{};
+ }
+ const double p = fabs(1. / track.qOverP);
+ const double px = p * cos(track.phi) * sin(track.theta);
+ const double py = p * sin(track.phi) * sin(track.theta);
+ const double pz = p * cos(track.theta);
+ return ROOT::Math::PxPyPzMVector{px, py, pz, mass};
+ });
+ return momenta;
}
- return {parts[first], parts[second]};
-}
-
-// Calculate the invariant mass of a given pair of particles
-inline double
-get_im(const std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>&
- particle_pair) {
- return (particle_pair.first + particle_pair.second).mass();
-}
-
-// Calculate the transverse momentum of a given pair of particles
-inline double
-get_pt(const std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>&
- particle_pair) {
- double px_pair = (particle_pair.first + particle_pair.second).px();
- double py_pair = (particle_pair.first + particle_pair.second).py();
- return sqrt(px_pair*px_pair + py_pair*py_pair);
-}
-
-// Calculate the azimuthal angle of a given pair of particles
-inline double
-get_phi(const std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>&
- particle_pair) {
- double px_pair = (particle_pair.first + particle_pair.second).px();
- double py_pair = (particle_pair.first + particle_pair.second).py();
- double phi_pair = std::atan2(py_pair,px_pair);
- //if(py_pair <= 0.) phi_pair = - phi_pair;
- return phi_pair;
-}
-
-// Calculate the rapidity of a given pair of particles
-inline double
-get_y(const std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>&
- particle_pair) {
- double px_pair = (particle_pair.first + particle_pair.second).px();
- double py_pair = (particle_pair.first + particle_pair.second).py();
- double pz_pair = (particle_pair.first + particle_pair.second).pz();
- double mass_pair = (particle_pair.first + particle_pair.second).mass();
- double energy_pair = sqrt(mass_pair*mass_pair + px_pair*px_pair + py_pair*py_pair + pz_pair*pz_pair);
- return 0.5*log((energy_pair + pz_pair)/(energy_pair - pz_pair));
-}
-
-//========================================================================================================
-//for structure functions
-
-struct inv_quant{ //add more when needed
- double nu, Q2, x;
-};
-
-//for simu
-inline inv_quant calc_inv_quant_simu(const std::vector<ROOT::Math::PxPyPzMVector>& parts){
- ROOT::Math::PxPyPzMVector q(parts[0] - parts[2]);
- ROOT::Math::PxPyPzMVector P(parts[3]);
-
- double nu = q.Dot(P) / P.mass();
- double Q2 = - q.Dot(q);
- inv_quant quantities = {nu, Q2, Q2/2./P.mass()/nu};
- return quantities;
-}
-
-inline double get_nu_simu(inv_quant quantities) {
- return quantities.nu/1000.;
-}
-inline double get_Q2_simu(inv_quant quantities) {
- return quantities.Q2;
-}
-inline double get_x_simu(inv_quant quantities) {
- return quantities.x;
-}
-
-//for tracking, add later
-
-//=========================================================================================================
+ // Get a vector of 4-momenta from the simulation data.
+ // TODO: Add PID selector (maybe using ranges?)
+ inline auto momenta_from_simulation(const std::vector<dd4pod::Geant4ParticleData>& parts)
+ {
+ std::vector<ROOT::Math::PxPyPzMVector> momenta{parts.size()};
+ // transform our simulation particle data into 4-momenta
+ std::transform(parts.begin(), parts.end(), momenta.begin(), [](const auto& part) {
+ return ROOT::Math::PxPyPzMVector{part.psx, part.psy, part.psz, part.mass};
+ });
+ return momenta;
+ }
+ // Find the decay pair candidates from a vector of particles (parts),
+ // with invariant mass closest to a desired value (pdg_mass)
+ inline std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>
+ find_decay_pair(const std::vector<ROOT::Math::PxPyPzMVector>& parts, const double pdg_mass)
+ {
+ int first = -1;
+ int second = -1;
+ double best_mass = -1;
+
+ // go through all particle combinatorics, calculate the invariant mass
+ // for each combination, and remember which combination is the closest
+ // to the desired pdg_mass
+ for (size_t i = 0; i < parts.size(); ++i) {
+ for (size_t j = i + 1; j < parts.size(); ++j) {
+ const double new_mass{(parts[i] + parts[j]).mass()};
+ if (fabs(new_mass - pdg_mass) < fabs(best_mass - pdg_mass)) {
+ first = i;
+ second = j;
+ best_mass = new_mass;
+ }
+ }
+ }
+ if (first < 0) {
+ return {{}, {}};
+ }
+ return {parts[first], parts[second]};
+ }
+ // Calculate the magnitude of the momentum of a vector of 4-vectors
+ inline auto mom(const std::vector<ROOT::Math::PxPyPzMVector>& momenta)
+ {
+ std::vector<double> P(momenta.size());
+ // transform our raw tracker info into proper 4-momenta
+ std::transform(momenta.begin(), momenta.end(), P.begin(),
+ [](const auto& mom) { return mom.P(); });
+ return P;
+ }
+ // Calculate the transverse momentum of a vector of 4-vectors
+ inline auto pt(const std::vector<ROOT::Math::PxPyPzMVector>& momenta)
+ {
+ std::vector<double> pt(momenta.size());
+ // transform our raw tracker info into proper 4-momenta
+ std::transform(momenta.begin(), momenta.end(), pt.begin(),
+ [](const auto& mom) { return mom.pt(); });
+ return pt;
+ }
+ // Calculate the azimuthal angle phi of a vector of 4-vectors
+ inline auto phi(const std::vector<ROOT::Math::PxPyPzMVector>& momenta)
+ {
+ std::vector<double> phi(momenta.size());
+ // transform our raw tracker info into proper 4-momenta
+ std::transform(momenta.begin(), momenta.end(), phi.begin(),
+ [](const auto& mom) { return mom.phi(); });
+ return phi;
+ }
+ // Calculate the pseudo-rapidity of a vector of particles
+ inline auto eta(const std::vector<ROOT::Math::PxPyPzMVector>& momenta)
+ {
+ std::vector<double> eta(momenta.size());
+ // transform our raw tracker info into proper 4-momenta
+ std::transform(momenta.begin(), momenta.end(), eta.begin(),
+ [](const auto& mom) { return mom.eta(); });
+ return eta;
+ }
+ //=========================================================================================================
} // namespace util
diff --git a/config/env.sh b/options/env.sh
similarity index 100%
rename from config/env.sh
rename to options/env.sh
diff --git a/config/tracker_reconstruction.py b/options/tracker_reconstruction.py
similarity index 100%
rename from config/tracker_reconstruction.py
rename to options/tracker_reconstruction.py
diff --git a/pythia_dis b/pythia_dis
new file mode 100755
index 0000000000000000000000000000000000000000..91a83f279189632193eaa87181405aaad17b8378
Binary files /dev/null and b/pythia_dis differ
diff --git a/tools/start_dev_shell.sh b/tools/dev-shell
similarity index 100%
rename from tools/start_dev_shell.sh
rename to tools/dev-shell
diff --git a/util/build_detector.sh b/util/build_detector.sh
index 03641f48151b95011700f80ffd1c0814417b47aa..02313585d56e0ca28af161bc19a7f7b3ae5f53c7 100755
--- a/util/build_detector.sh
+++ b/util/build_detector.sh
@@ -18,9 +18,9 @@ pushd ${PROJECT_ROOT}
## - DETECTOR_PATH: full path for the detector definitions
## this is the same as ${DETECTOR_PREFIX}/${JUGGLER_DETECTOR}
##
-## You can read config/env.sh for more in-depth explanations of the variables
+## You can read options/env.sh for more in-depth explanations of the variables
## and how they can be controlled.
-source config/env.sh
+source options/env.sh
## =============================================================================
## Step 1: download/update the detector definitions (if needed)