Skip to content
Snippets Groups Projects
Commit 1a929908 authored by Sylvester Joosten's avatar Sylvester Joosten
Browse files

restructuring in progress, still have to do dvcs

parent c8780170
No related branches found
No related tags found
1 merge request!10Update dis
This commit is part of merge request !10. Comments created here will be created in the context of that merge request.
Show whitespace changes
Inline Side-by-side
Showing
with 647 additions and 453 deletions
.clang-format
+ 1
− 1
View file @ 1a929908
...@@ -44,7 +44,7 @@ BreakConstructorInitializersBeforeComma: false ...@@ -44,7 +44,7 @@ BreakConstructorInitializersBeforeComma: false
BreakConstructorInitializers: BeforeColon BreakConstructorInitializers: BeforeColon
BreakAfterJavaFieldAnnotations: false BreakAfterJavaFieldAnnotations: false
BreakStringLiterals: true BreakStringLiterals: true
ColumnLimit: 120 ColumnLimit: 100
CommentPragmas: '^ IWYU pragma:' CommentPragmas: '^ IWYU pragma:'
CompactNamespaces: false CompactNamespaces: false
ConstructorInitializerAllOnOneLineOrOnePerLine: false ConstructorInitializerAllOnOneLineOrOnePerLine: false
......
.rootlogon.C
+ 1
− 0
View file @ 1a929908
...@@ -7,6 +7,7 @@ ...@@ -7,6 +7,7 @@
// setup a local build directory so we don't polute our source code with // setup a local build directory so we don't polute our source code with
// ROOT dictionaries etc. // ROOT dictionaries etc.
gSystem->Exec("mkdir -p /tmp/root_build");
gSystem->SetBuildDir("/tmp/root_build"); gSystem->SetBuildDir("/tmp/root_build");
// style definition based off the ATLAS style // style definition based off the ATLAS style
......
accelerator @ f3ff428e
Subproject commit f3ff428e3b926a41e95beaa984d8dc05cec37cc7
benchmarks/dis/analysis/dis.h 0 → 100644
+ 26
− 0
View file @ 1a929908
#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
benchmarks/dis/analysis/dis_electrons.cxx 0 → 100644
+ 116
− 0
View file @ 1a929908
#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;
}
benchmarks/dis/analysis/rec_dis_electrons.cxx deleted 100644 → 0
+ 0
− 113
View file @ c8780170
#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;
}
benchmarks/dis/dis.sh 100644 → 100755
+ 82
− 41
View file @ 1a929908
#!/bin/bash #!/bin/bash
## ============================================================================= ## =============================================================================
## Run the DVMP benchmarks in 5 steps: ## Run the DIS benchmarks in 5 steps:
## 1. Build/install detector package ## 1. Parse the command line and setup environment
## 2. Detector simulation through npsim ## 2. Detector simulation through npsim
## 3. Digitization and reconstruction through Juggler ## 3. Digitization and reconstruction through Juggler
## 4. Root-based Physics analyses ## 4. Root-based Physics analyses
## 5. Finalize ## 5. Finalize
## ============================================================================= ## =============================================================================
echo "Running the DIS benchmarks"
## make sure we launch this script from the project root directory ## 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} pushd ${PROJECT_ROOT}
echo "Running the DIS benchmarks"
## ============================================================================= ## =============================================================================
## Load the environment variables. To build the detector we need the following ## Step 1: Setup the environment variables
## 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_INSTALL_PREFIX: Install prefix for Juggler (simu/recon)
## - JUGGLER_DETECTOR: the detector package we want to use for this benchmark ## - JUGGLER_DETECTOR: the detector package we want to use for this benchmark
## - DETECTOR_PATH: full path to the detector definitions ## - 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. ## and how they can be controlled.
source config/env.sh source options/env.sh
## Extra environment variables for DVMP: ## We also need the following benchmark-specific variables:
## file tag for these tests ##
JUGGLER_FILE_NAME_TAG="dis" ## - BENCHMARK_TAG: Unique identified for this benchmark process.
# TODO use the input file name, as we will be generating a lot of these ## - BEAM_TAG: Identifier for the chosen beam configuration
# in the future... ## - INPUT_PATH: Path for generator-level input to the benchmarks
# FIXME Generator file hardcoded for now ## - TMP_PATH: Path for temporary data (not exported as artifacts)
## note: these variables need to be exported to be accessible from ## - RESULTS_PATH: Path for benchmark output figures and files
## the juggler options.py. We should really work on a dedicated ##
## juggler launcher to get rid of these "magic" variables. FIXME ## You can read dvmp/env.sh for more in-depth explanations of the variables.
export JUGGLER_GEN_FILE="results/dis/${JUGGLER_FILE_NAME_TAG}.hepmc" source benchmarks/dis/env.sh
export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root"
export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
## 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
## Step 1: Build/install the desired detector package SIM_LOG=${TMP_PATH}/sim-${CONFIG}.log
## 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 ## Step 2: Run the simulation
...@@ -55,8 +65,8 @@ npsim --runType batch \ ...@@ -55,8 +65,8 @@ npsim --runType batch \
-v WARNING \ -v WARNING \
--numberOfEvents ${JUGGLER_N_EVENTS} \ --numberOfEvents ${JUGGLER_N_EVENTS} \
--compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \ --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
--inputFiles ${JUGGLER_GEN_FILE} \ --inputFiles ${GEN_FILE} \
--outputFile ${JUGGLER_SIM_FILE} --outputFile ${SIM_FILE}
if [ "$?" -ne "0" ] ; then if [ "$?" -ne "0" ] ; then
echo "ERROR running npsim" echo "ERROR running npsim"
exit 1 exit 1
...@@ -65,36 +75,67 @@ fi ...@@ -65,36 +75,67 @@ fi
## ============================================================================= ## =============================================================================
## Step 3: Run digitization & reconstruction ## Step 3: Run digitization & reconstruction
echo "Running the digitization and 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 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 \ 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 if [ "$?" -ne "0" ] ; then
echo "ERROR running juggler" echo "ERROR running juggler, both attempts failed"
exit 1 exit 1
fi fi
ls -l fi
## ============================================================================= ## =============================================================================
## Step 4: Analysis ## 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 if [[ "$?" -ne "0" ]] ; then
echo "ERROR running root script" echo "ERROR running rec_dis_electron script"
exit 1 exit 1
fi fi
## ============================================================================= ## =============================================================================
## Step 5: finalize ## 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 ## too many events
if [ "${JUGGLER_N_EVENTS}" -lt "500" ] ; then if [ "${JUGGLER_N_EVENTS}" -lt "500" ] ; then
cp ${JUGGLER_REC_FILE} results/dis/. cp ${REC_FILE} ${RESULTS_PATH}
fi fi
## cleanup output files ## Always move over log files to the results path
rm ${JUGGLER_REC_FILE} ${JUGGLER_SIM_FILE} cp ${REC_LOG} ${RESULTS_PATH}
## ============================================================================= ## =============================================================================
## All done! ## All done!
echo "${JUGGLER_FILE_NAME_TAG} benchmarks complete" echo "DIS benchmarks complete"
benchmarks/dis/env.sh 0 → 100644
+ 50
− 0
View file @ 1a929908
#!/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."
benchmarks/dis/gen.sh 100644 → 100755
+ 56
− 25
View file @ 1a929908
...@@ -3,38 +3,63 @@ ...@@ -3,38 +3,63 @@
## ============================================================================= ## =============================================================================
## Standin for a proper pythia generation process, similar to how we ## Standin for a proper pythia generation process, similar to how we
## generate events for DVMP ## 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 ## 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} pushd ${PROJECT_ROOT}
## ============================================================================= ## =============================================================================
## Load the environment variables. To build the detector we need the following ## Step 1: Setup the environment variables
## 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) ## - LOCAL_PREFIX: Place to cache local packages and data
## - JUGGLER_DETECTOR: the detector package we want to use for this benchmark ## - JUGGLER_N_EVENTS: Number of events to process
## - DETECTOR_PATH: full path to the detector definitions ## - 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. ## 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 ## Get a unique file name prefix based on the configuration options
export DATA_PATH=results/dis GEN_TAG=gen-${CONFIG}_${JUGGLER_N_EVENTS} ## Generic file prefix
## Extra environment variables for DVMP: ## =============================================================================
## file tag for these tests ## Step 2: Check if we really need to run, or can use the cache.
JUGGLER_FILE_NAME_TAG="dis" 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 ## Step 3: Build generator exe
## TODO better file name that encodes the actual configuration we're running ## TODO: need to configurability to the generator exe
echo "Compiling benchmarks/dis/generator/pythia_dis.cxx ..." echo "Compiling benchmarks/dis/generator/pythia_dis.cxx ..."
g++ benchmarks/dis/generator/pythia_dis.cxx -o pythia_dis \ g++ benchmarks/dis/generator/pythia_dis.cxx -o pythia_dis \
-I/usr/local/include -Iinclude \ -I/usr/local/include -Iinclude \
...@@ -46,20 +71,26 @@ if [[ "$?" -ne "0" ]] ; then ...@@ -46,20 +71,26 @@ if [[ "$?" -ne "0" ]] ; then
exit 1 exit 1
fi fi
echo "done" 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 if [[ "$?" -ne "0" ]] ; then
echo "ERROR running pythia" echo "ERROR running pythia"
exit 1 exit 1
fi fi
## =============================================================================
## Step 5: Finally, move relevant output into the artifacts directory and clean up
## ============================================================================= ## =============================================================================
## Step 2: finalize echo "Moving generator output into ${INPUT_PATH}"
echo "Moving event generator output into ${DATA_PATH}" mv ${TMP_PATH}/${GEN_TAG}.hepmc ${INPUT_PATH}/${GEN_TAG}.hepmc
#mv .local/${JUGGLER_FILE_NAME_TAG}.hepmc ${DATA_PATH}/${JUGGLER_FILE_NAME_TAG}.hepmc ## this step only matters for local execution
echo "Cleaning up"
## does nothing
## ============================================================================= ## =============================================================================
## All done! ## All done!
echo "dis event generation complete" echo "$BENCHMARK_TAG event generation complete"
benchmarks/dis/generator/pythia_dis.cxx
+ 12
− 12
View file @ 1a929908
#include "Pythia8/Pythia.h" #include "Pythia8/Pythia.h"
#include "Pythia8Plugins/HepMC3.h" #include "Pythia8Plugins/HepMC3.h"
#include <unistd.h>
#include <cassert> #include <cassert>
#include <unistd.h>
using namespace Pythia8; using namespace Pythia8;
...@@ -13,7 +13,7 @@ using std::string; ...@@ -13,7 +13,7 @@ using std::string;
enum class mode { none, help, list, part }; enum class mode { none, help, list, part };
struct settings { struct settings {
double E_electron = 10.0; // GeV double E_electron = 18.0; // GeV
double E_ion = 275.0; // GeV double E_ion = 275.0; // GeV
std::string outfile = "dis.hepmc"; std::string outfile = "dis.hepmc";
int ion_PID = 2212; int ion_PID = 2212;
...@@ -63,16 +63,15 @@ void print_man_page(T cli, const char* argv0) ...@@ -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; settings s;
auto lastOpt = auto lastOpt = " options:" % (option("-h", "--help").set(s.selected, mode::help) % "show help",
" options:" % (option("-h", "--help").set(s.selected, mode::help) % "show help",
value("file", s.outfile).if_missing([] { value("file", s.outfile).if_missing([] {
std::cout << "You need to provide an output filename argument!\n"; std::cout << "You need to provide an output filename argument!\n";
}) % "output file"); }) % "output file");
auto cli = auto cli = (command("help").set(s.selected, mode::help) | lastOpt);
(command("help").set(s.selected, mode::help) | lastOpt);
assert(cli.flags_are_prefix_free()); assert(cli.flags_are_prefix_free());
...@@ -95,7 +94,8 @@ settings cmdline_settings(int argc, char* argv[]) { ...@@ -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); settings s = cmdline_settings(argc, argv);
if (!s.success) { if (!s.success) {
......
benchmarks/dvcs/dvcs.sh
+ 3
− 3
View file @ 1a929908
...@@ -15,9 +15,9 @@ echo "JUGGLER_FILE_NAME_TAG = ${JUGGLER_FILE_NAME_TAG}" ...@@ -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 ## - JUGGLER_DETECTOR: the detector package we want to use for this benchmark
## - DETECTOR_PATH: full path to the detector definitions ## - 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. ## 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" 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 ...@@ -43,7 +43,7 @@ fi
# Need to figure out how to pass file name to juggler from the commandline # Need to figure out how to pass file name to juggler from the commandline
xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \ xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
gaudirun.py config/tracker_reconstruction.py gaudirun.py options/tracker_reconstruction.py
if [[ "$?" -ne "0" ]] ; then if [[ "$?" -ne "0" ]] ; then
echo "ERROR running juggler" echo "ERROR running juggler"
exit 1 exit 1
......
benchmarks/dvmp/analysis/dvmp.h 0 → 100644
+ 56
− 0
View file @ 1a929908
#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
benchmarks/dvmp/analysis/plot.h deleted 100644 → 0
+ 0
− 40
View file @ c8780170
#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
benchmarks/dvmp/analysis/vm_invar.cxx
+ 14
− 12
View file @ 1a929908
#include "dvmp.h"
#include "plot.h" #include "plot.h"
#include <benchmark.h> #include <benchmark.h>
...@@ -18,8 +19,6 @@ ...@@ -18,8 +19,6 @@
// a desired vector meson (e.g. jpsi) and a desired decay particle (e.g. muon) // 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 // Output figures are written to our output prefix (which includes the output
// file prefix), and labeled with our detector name. // 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 // FIXME: MC does not trace back into particle history. Need to fix that
int vm_invar(const std::string& config_name) int vm_invar(const std::string& config_name)
{ {
...@@ -35,17 +34,20 @@ 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"]; std::string output_prefix = config["output_prefix"];
const std::string test_tag = config["test_tag"]; 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(" - Vector meson: {}\n", vm_name);
fmt::print(" - Decay particle: {}\n", decay_name); fmt::print(" - Decay particle: {}\n", decay_name);
fmt::print(" - Detector package: {}\n", detector); fmt::print(" - Detector package: {}\n", detector);
fmt::print(" - input file: {}\n", rec_file);
fmt::print(" - output prefix: {}\n", output_prefix); fmt::print(" - output prefix: {}\n", output_prefix);
// create our test definition // create our test definition
// test_tag // test_tag
eic::util::Test vm_mass_resolution_test{ eic::util::Test Q2_resolution_test{
{{"name", fmt::format("{}_{}_{}_mass_resolution", test_tag, vm_name, decay_name)}, {{"name", fmt::format("{}_{}_{}_Q2_resolution", test_tag, vm_name, decay_name)},
{"title", fmt::format("{} -> {} Invariant Mass Resolution", 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 " {"description", "Invariant Mass Resolution calculated from raw "
"tracking data using a Gaussian fit."}, "tracking data using a Gaussian fit."},
{"quantity", "resolution"}, {"quantity", "resolution"},
...@@ -114,7 +116,7 @@ int vm_invar(const std::string& config_name) ...@@ -114,7 +116,7 @@ int vm_invar(const std::string& config_name)
// draw everything // draw everything
hnu.DrawClone("hist"); hnu.DrawClone("hist");
// FIXME hardcoded beam configuration // FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector, vm_name, "#nu"); plot::draw_label(10, 100, detector);
TText* tptr21; TText* tptr21;
auto t21 = new TPaveText(.6, .8417, .9, .925, "NB NDC"); auto t21 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t21->SetFillColorAlpha(kWhite, 0); t21->SetFillColorAlpha(kWhite, 0);
...@@ -139,7 +141,7 @@ int vm_invar(const std::string& config_name) ...@@ -139,7 +141,7 @@ int vm_invar(const std::string& config_name)
// draw everything // draw everything
hQ2.DrawClone("hist"); hQ2.DrawClone("hist");
// FIXME hardcoded beam configuration // FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector, vm_name, "Q^{2}"); plot::draw_label(10, 100, detector);
TText* tptr22; TText* tptr22;
auto t22 = new TPaveText(.6, .8417, .9, .925, "NB NDC"); auto t22 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t22->SetFillColorAlpha(kWhite, 0); t22->SetFillColorAlpha(kWhite, 0);
...@@ -164,7 +166,7 @@ int vm_invar(const std::string& config_name) ...@@ -164,7 +166,7 @@ int vm_invar(const std::string& config_name)
// draw everything // draw everything
hx.DrawClone("hist"); hx.DrawClone("hist");
// FIXME hardcoded beam configuration // FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector, vm_name, "x"); plot::draw_label(10, 100, detector);
TText* tptr23; TText* tptr23;
auto t23 = new TPaveText(.6, .8417, .9, .925, "NB NDC"); auto t23 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t23->SetFillColorAlpha(kWhite, 0); t23->SetFillColorAlpha(kWhite, 0);
...@@ -179,12 +181,12 @@ int vm_invar(const std::string& config_name) ...@@ -179,12 +181,12 @@ int vm_invar(const std::string& config_name)
c.Print(fmt::format("{}InvariantQuantities.png", output_prefix).c_str()); 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 // error for the test result
vm_mass_resolution_test.error(-1); Q2_resolution_test.error(-1);
// write out our test data // 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! // That's all!
return 0; return 0;
......
benchmarks/dvmp/analysis/vm_mass.cxx
+ 32
− 23
View file @ 1a929908
#include "dvmp.h"
#include "plot.h" #include "plot.h"
#include <benchmark.h> #include <benchmark.h>
...@@ -35,17 +36,20 @@ int vm_mass(const std::string& config_name) ...@@ -35,17 +36,20 @@ int vm_mass(const std::string& config_name)
std::string output_prefix = config["output_prefix"]; std::string output_prefix = config["output_prefix"];
const std::string test_tag = config["test_tag"]; 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(" - Vector meson: {}\n", vm_name);
fmt::print(" - Decay particle: {}\n", decay_name); fmt::print(" - Decay particle: {}\n", decay_name);
fmt::print(" - Detector package: {}\n", detector); fmt::print(" - Detector package: {}\n", detector);
fmt::print(" - input file: {}\n", rec_file);
fmt::print(" - output prefix: {}\n", output_prefix); fmt::print(" - output prefix: {}\n", output_prefix);
// create our test definition // create our test definition
// test_tag // 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)}, {{"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 " {"description", "Invariant Mass Resolution calculated from raw "
"tracking data using a Gaussian fit."}, "tracking data using a Gaussian fit."},
{"quantity", "resolution"}, {"quantity", "resolution"},
...@@ -82,18 +86,23 @@ int vm_mass(const std::string& config_name) ...@@ -82,18 +86,23 @@ int vm_mass(const std::string& config_name)
.Define("p_sim", util::momenta_from_simulation, {"mcparticles2"}) .Define("p_sim", util::momenta_from_simulation, {"mcparticles2"})
.Define("decay_pair_rec", find_decay_pair, {"p_rec"}) .Define("decay_pair_rec", find_decay_pair, {"p_rec"})
.Define("decay_pair_sim", find_decay_pair, {"p_sim"}) .Define("decay_pair_sim", find_decay_pair, {"p_sim"})
.Define("mass_rec", util::get_im, {"decay_pair_rec"}) .Define("p_vm_rec", "decay_pair_rec.first + decay_pair_rec.second")
.Define("mass_sim", util::get_im, {"decay_pair_sim"}) .Define("p_vm_sim", "decay_pair_sim.first + decay_pair_sim.second")
.Define("pt_rec", util::get_pt, {"decay_pair_rec"}) //.Define("p_vm_sim", util::get_sum, {"decay_pair_sim"})
.Define("pt_sim", util::get_pt, {"decay_pair_sim"}) .Define("mass_rec", "p_vm_rec.M()")
.Define("phi_rec", util::get_phi, {"decay_pair_rec"}) .Define("mass_sim", "p_vm_sim.M()")
.Define("phi_sim", util::get_phi, {"decay_pair_sim"}) .Define("pt_rec", "p_vm_rec.pt()")
.Define("rapidity_rec", util::get_y, {"decay_pair_rec"}) .Define("pt_sim", "p_vm_sim.pt()")
.Define("rapidity_sim", util::get_y, {"decay_pair_sim"}); .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 // 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_rec =
auto h_im_sim = d_im.Histo1D({"h_im_sim", ";m_{ll'} (GeV/c^{2});#", 100, -1.1, vm_mass + 5}, "mass_sim"); 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_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"); 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) ...@@ -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_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_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_eta_rec = d_im.Histo1D({"h_eta_rec", ";#eta_{ll'};#", 1000, -5., 5.}, "eta_rec");
auto h_y_sim = d_im.Histo1D({"h_y_sim", ";y_{ll'};#", 1000, -5., 5.}, "rapidity_sim"); auto h_eta_sim = d_im.Histo1D({"h_eta_sim", ";#eta_{ll'};#", 1000, -5., 5.}, "eta_sim");
// Plot our histograms. // Plot our histograms.
// TODO: to start I'm explicitly plotting the histograms, but want to // TODO: to start I'm explicitly plotting the histograms, but want to
...@@ -128,7 +137,7 @@ int vm_mass(const std::string& config_name) ...@@ -128,7 +137,7 @@ int vm_mass(const std::string& config_name)
h11.DrawClone("hist"); h11.DrawClone("hist");
h12.DrawClone("hist same"); h12.DrawClone("hist same");
// FIXME hardcoded beam configuration // FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector, vm_name, "Invariant mass"); plot::draw_label(10, 100, detector);
TText* tptr1; TText* tptr1;
auto t1 = new TPaveText(.6, .8417, .9, .925, "NB NDC"); auto t1 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t1->SetFillColorAlpha(kWhite, 0); t1->SetFillColorAlpha(kWhite, 0);
...@@ -158,7 +167,7 @@ int vm_mass(const std::string& config_name) ...@@ -158,7 +167,7 @@ int vm_mass(const std::string& config_name)
h21.DrawClone("hist"); h21.DrawClone("hist");
h22.DrawClone("hist same"); h22.DrawClone("hist same");
// FIXME hardcoded beam configuration // FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector, vm_name, "Transverse Momentum"); plot::draw_label(10, 100, detector);
TText* tptr2; TText* tptr2;
auto t2 = new TPaveText(.6, .8417, .9, .925, "NB NDC"); auto t2 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t2->SetFillColorAlpha(kWhite, 0); t2->SetFillColorAlpha(kWhite, 0);
...@@ -188,7 +197,7 @@ int vm_mass(const std::string& config_name) ...@@ -188,7 +197,7 @@ int vm_mass(const std::string& config_name)
h31.DrawClone("hist"); h31.DrawClone("hist");
h32.DrawClone("hist same"); h32.DrawClone("hist same");
// FIXME hardcoded beam configuration // FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector, vm_name, "#phi"); plot::draw_label(10, 100, detector);
TText* tptr3; TText* tptr3;
auto t3 = new TPaveText(.6, .8417, .9, .925, "NB NDC"); auto t3 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t3->SetFillColorAlpha(kWhite, 0); t3->SetFillColorAlpha(kWhite, 0);
...@@ -204,8 +213,8 @@ int vm_mass(const std::string& config_name) ...@@ -204,8 +213,8 @@ int vm_mass(const std::string& config_name)
c.cd(4); c.cd(4);
// gPad->SetLogx(false); // gPad->SetLogx(false);
// gPad->SetLogy(false); // gPad->SetLogy(false);
auto& h41 = *h_y_sim; auto& h41 = *h_eta_sim;
auto& h42 = *h_y_rec; auto& h42 = *h_eta_rec;
// histogram style // histogram style
h41.SetLineColor(plot::kMpBlue); h41.SetLineColor(plot::kMpBlue);
h41.SetLineWidth(2); h41.SetLineWidth(2);
...@@ -218,7 +227,7 @@ int vm_mass(const std::string& config_name) ...@@ -218,7 +227,7 @@ int vm_mass(const std::string& config_name)
h41.DrawClone("hist"); h41.DrawClone("hist");
h42.DrawClone("hist same"); h42.DrawClone("hist same");
// FIXME hardcoded beam configuration // FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector, vm_name, "Rapidity"); plot::draw_label(10, 100, detector);
TText* tptr4; TText* tptr4;
auto t4 = new TPaveText(.6, .8417, .9, .925, "NB NDC"); auto t4 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t4->SetFillColorAlpha(kWhite, 0); t4->SetFillColorAlpha(kWhite, 0);
...@@ -235,10 +244,10 @@ int vm_mass(const std::string& config_name) ...@@ -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 // TODO we're not actually doing an IM fit yet, so for now just return an
// error for the test result // error for the test result
vm_mass_resolution_test.error(-1); mass_resolution_test.error(-1);
// write out our test data // 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! // That's all!
return 0; return 0;
......
benchmarks/dvmp/dvmp.sh
+ 13
− 10
View file @ 1a929908
...@@ -35,9 +35,9 @@ source util/parse_cmd.sh $@ ...@@ -35,9 +35,9 @@ source util/parse_cmd.sh $@
## - JUGGLER_DETECTOR: the detector package we want to use for this benchmark ## - JUGGLER_DETECTOR: the detector package we want to use for this benchmark
## - DETECTOR_PATH: full path to the detector definitions ## - 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. ## and how they can be controlled.
source config/env.sh source options/env.sh
## We also need the following benchmark-specific variables: ## We also need the following benchmark-specific variables:
## ##
...@@ -93,21 +93,20 @@ export JUGGLER_SIM_FILE=${SIM_FILE} ...@@ -93,21 +93,20 @@ export JUGGLER_SIM_FILE=${SIM_FILE}
export JUGGLER_REC_FILE=${REC_FILE} export JUGGLER_REC_FILE=${REC_FILE}
export JUGGLER_DETECTOR_PATH=${DETECTOR_PATH} export JUGGLER_DETECTOR_PATH=${DETECTOR_PATH}
xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \ xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
gaudirun.py config/tracker_reconstruction.py \ gaudirun.py options/tracker_reconstruction.py \
2>&1 > ${REC_LOG} 2>&1 > ${REC_LOG}
## on-error, first retry running juggler again as there is still a random ## on-error, first retry running juggler again as there is still a random
## crash we need to address FIXME ## crash we need to address FIXME
if [ "$?" -ne "0" ] ; then if [ "$?" -ne "0" ] ; then
echo "Juggler crashed, retrying..." echo "Juggler crashed, retrying..."
xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \ xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
gaudirun.py config/tracker_reconstruction.py \ gaudirun.py options/tracker_reconstruction.py \
2>&1 > ${REC_LOG} 2>&1 > ${REC_LOG}
if [ "$?" -ne "0" ] ; then if [ "$?" -ne "0" ] ; then
echo "ERROR running juggler, both attempts failed" echo "ERROR running juggler, both attempts failed"
exit 1 exit 1
fi fi
fi fi
#ls -l
## ============================================================================= ## =============================================================================
## Step 4: Analysis ## Step 4: Analysis
...@@ -127,10 +126,14 @@ EOF ...@@ -127,10 +126,14 @@ EOF
#cat ${CONFIG} #cat ${CONFIG}
## run the analysis script with this configuration ## run the analysis script with this configuration
root -b -q "benchmarks/dvmp/analysis/vm_mass.cxx(\"${CONFIG}\")" root -b -q "benchmarks/dvmp/analysis/vm_mass.cxx+(\"${CONFIG}\")"
root -b -q "benchmarks/dvmp/analysis/vm_invar.cxx(\"${CONFIG}\")"
if [ "$?" -ne "0" ] ; then 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 exit 1
fi fi
...@@ -145,11 +148,11 @@ if [ "${JUGGLER_N_EVENTS}" -lt "500" ] ; then ...@@ -145,11 +148,11 @@ if [ "${JUGGLER_N_EVENTS}" -lt "500" ] ; then
fi fi
## Always move over log files to the results path ## Always move over log files to the results path
mv ${SIM_LOG} ${REC_LOG} ${RESULTS_PATH} mv ${REC_LOG} ${RESULTS_PATH}
## cleanup output files ## cleanup output files
rm -f ${REC_FILE} ${SIM_FILE} #rm -f ${REC_FILE} ${SIM_FILE} ## --> not needed for CI
## ============================================================================= ## =============================================================================
## All done! ## All done!
......
benchmarks/dvmp/gen.sh
+ 4
− 2
View file @ 1a929908
...@@ -33,9 +33,9 @@ source util/parse_cmd.sh $@ ...@@ -33,9 +33,9 @@ source util/parse_cmd.sh $@
## - JUGGLER_N_EVENTS: Number of events to process ## - JUGGLER_N_EVENTS: Number of events to process
## - JUGGLER_RNG_SEED: Random seed for event generation. ## - 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. ## and how they can be controlled.
source config/env.sh source options/env.sh
## We also need the following benchmark-specific variables: ## We also need the following benchmark-specific variables:
## ##
...@@ -107,4 +107,6 @@ done ...@@ -107,4 +107,6 @@ done
echo "Cleaning up" echo "Cleaning up"
rm ${TMP_PATH}/${GEN_TAG}.json rm ${TMP_PATH}/${GEN_TAG}.json
## =============================================================================
## All done! ## All done!
echo "$BENCHMARK_TAG event generation complete"
include/benchmark.h
+ 9
− 3
View file @ 1a929908
...@@ -4,6 +4,7 @@ ...@@ -4,6 +4,7 @@
#include "exception.h" #include "exception.h"
#include <fmt/core.h> #include <fmt/core.h>
#include <fstream> #include <fstream>
#include <iomanip>
#include <iostream> #include <iostream>
#include <nlohmann/json.hpp> #include <nlohmann/json.hpp>
#include <string> #include <string>
...@@ -67,14 +68,19 @@ namespace eic::util { ...@@ -67,14 +68,19 @@ namespace eic::util {
// - weight: Weight for this test (this is defaulted to 1.0 if not specified) // - weight: Weight for this test (this is defaulted to 1.0 if not specified)
// - result: pass/fail/error // - result: pass/fail/error
struct Test { 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) // initialize with error (as we don't have a value yet)
error(); error();
// Check that all required fields are present // 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()) { 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 // Default "weight" to 1 if not set
......
include/plot.h 0 → 100644
+ 42
− 0
View file @ 1a929908
#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
include/util.h
+ 129
− 168
View file @ 1a929908
#ifndef UTIL_H #ifndef UTIL_H
#define UTIL_H #define UTIL_H
// TODO: should probably be moved to a global benchmark utility library
#include <algorithm> #include <algorithm>
#include <cmath> #include <cmath>
#include <exception> #include <exception>
...@@ -22,7 +24,8 @@ namespace util { ...@@ -22,7 +24,8 @@ namespace util {
class unknown_particle_error : public std::exception { class unknown_particle_error : public std::exception {
public: public:
unknown_particle_error(std::string_view particle) : m_particle{particle} {} unknown_particle_error(std::string_view particle) : m_particle{particle} {}
virtual const char* what() const throw() { virtual const char* what() const throw()
{
return fmt::format("Unknown particle type: {}", m_particle).c_str(); return fmt::format("Unknown particle type: {}", m_particle).c_str();
} }
virtual const char* type() const throw() { return "unknown_particle_error"; } virtual const char* type() const throw() { return "unknown_particle_error"; }
...@@ -35,7 +38,8 @@ private: ...@@ -35,7 +38,8 @@ private:
// we care about for this process. // we care about for this process.
// FIXME: consider something more robust (maybe based on hepPDT) to the // FIXME: consider something more robust (maybe based on hepPDT) to the
// analysis utility library // analysis utility library
inline double get_pdg_mass(std::string_view part) { inline double get_pdg_mass(std::string_view part)
{
if (part == "electron") { if (part == "electron") {
return 0.0005109989461; return 0.0005109989461;
} else if (part == "muon") { } else if (part == "muon") {
...@@ -51,14 +55,12 @@ inline double get_pdg_mass(std::string_view part) { ...@@ -51,14 +55,12 @@ inline double get_pdg_mass(std::string_view part) {
// Get a vector of 4-momenta from raw tracking info, using an externally // Get a vector of 4-momenta from raw tracking info, using an externally
// provided particle mass assumption. // provided particle mass assumption.
// FIXME: should be part of utility library inline auto momenta_from_tracking(const std::vector<eic::TrackParametersData>& tracks,
inline auto const double mass)
momenta_from_tracking(const std::vector<eic::TrackParametersData>& tracks, {
const double mass) {
std::vector<ROOT::Math::PxPyPzMVector> momenta{tracks.size()}; std::vector<ROOT::Math::PxPyPzMVector> momenta{tracks.size()};
// transform our raw tracker info into proper 4-momenta // transform our raw tracker info into proper 4-momenta
std::transform(tracks.begin(), tracks.end(), momenta.begin(), std::transform(tracks.begin(), tracks.end(), momenta.begin(), [mass](const auto& track) {
[mass](const auto& track) {
// make sure we don't divide by zero // make sure we don't divide by zero
if (fabs(track.qOverP) < 1e-9) { if (fabs(track.qOverP) < 1e-9) {
return ROOT::Math::PxPyPzMVector{}; return ROOT::Math::PxPyPzMVector{};
...@@ -73,27 +75,22 @@ momenta_from_tracking(const std::vector<eic::TrackParametersData>& tracks, ...@@ -73,27 +75,22 @@ momenta_from_tracking(const std::vector<eic::TrackParametersData>& tracks,
} }
// Get a vector of 4-momenta from the simulation data. // Get a vector of 4-momenta from the simulation data.
// FIXME: should be part of utility library
// TODO: Add PID selector (maybe using ranges?) // TODO: Add PID selector (maybe using ranges?)
inline auto inline auto momenta_from_simulation(const std::vector<dd4pod::Geant4ParticleData>& parts)
momenta_from_simulation(const std::vector<dd4pod::Geant4ParticleData>& parts) { {
std::vector<ROOT::Math::PxPyPzMVector> momenta{parts.size()}; std::vector<ROOT::Math::PxPyPzMVector> momenta{parts.size()};
// transform our simulation particle data into 4-momenta // transform our simulation particle data into 4-momenta
std::transform(parts.begin(), parts.end(), momenta.begin(), std::transform(parts.begin(), parts.end(), momenta.begin(), [](const auto& part) {
[](const auto& part) { return ROOT::Math::PxPyPzMVector{part.psx, part.psy, part.psz, part.mass};
return ROOT::Math::PxPyPzMVector{part.psx, part.psy,
part.psz, part.mass};
}); });
return momenta; return momenta;
} }
// Find the decay pair candidates from a vector of particles (parts), // Find the decay pair candidates from a vector of particles (parts),
// with invariant mass closest to a desired value (pdg_mass) // 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> inline std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>
find_decay_pair(const std::vector<ROOT::Math::PxPyPzMVector>& parts, find_decay_pair(const std::vector<ROOT::Math::PxPyPzMVector>& parts, const double pdg_mass)
const double pdg_mass) { {
int first = -1; int first = -1;
int second = -1; int second = -1;
double best_mass = -1; double best_mass = -1;
...@@ -101,8 +98,8 @@ find_decay_pair(const std::vector<ROOT::Math::PxPyPzMVector>& parts, ...@@ -101,8 +98,8 @@ find_decay_pair(const std::vector<ROOT::Math::PxPyPzMVector>& parts,
// go through all particle combinatorics, calculate the invariant mass // go through all particle combinatorics, calculate the invariant mass
// for each combination, and remember which combination is the closest // for each combination, and remember which combination is the closest
// to the desired pdg_mass // to the desired pdg_mass
for (int i = 0; i < parts.size(); ++i) { for (size_t i = 0; i < parts.size(); ++i) {
for (int j = i + 1; j < parts.size(); ++j) { for (size_t j = i + 1; j < parts.size(); ++j) {
const double new_mass{(parts[i] + parts[j]).mass()}; const double new_mass{(parts[i] + parts[j]).mass()};
if (fabs(new_mass - pdg_mass) < fabs(best_mass - pdg_mass)) { if (fabs(new_mass - pdg_mass) < fabs(best_mass - pdg_mass)) {
first = i; first = i;
...@@ -117,82 +114,46 @@ find_decay_pair(const std::vector<ROOT::Math::PxPyPzMVector>& parts, ...@@ -117,82 +114,46 @@ find_decay_pair(const std::vector<ROOT::Math::PxPyPzMVector>& parts,
return {parts[first], parts[second]}; return {parts[first], parts[second]};
} }
// Calculate the invariant mass of a given pair of particles // Calculate the magnitude of the momentum of a vector of 4-vectors
inline double inline auto mom(const std::vector<ROOT::Math::PxPyPzMVector>& momenta)
get_im(const std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>& {
particle_pair) { std::vector<double> P(momenta.size());
return (particle_pair.first + particle_pair.second).mass(); // transform our raw tracker info into proper 4-momenta
} std::transform(momenta.begin(), momenta.end(), P.begin(),
[](const auto& mom) { return mom.P(); });
// Calculate the transverse momentum of a given pair of particles return P;
inline double }
get_pt(const std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>& // Calculate the transverse momentum of a vector of 4-vectors
particle_pair) { inline auto pt(const std::vector<ROOT::Math::PxPyPzMVector>& momenta)
double px_pair = (particle_pair.first + particle_pair.second).px(); {
double py_pair = (particle_pair.first + particle_pair.second).py(); std::vector<double> pt(momenta.size());
return sqrt(px_pair*px_pair + py_pair*py_pair); // transform our raw tracker info into proper 4-momenta
} std::transform(momenta.begin(), momenta.end(), pt.begin(),
[](const auto& mom) { return mom.pt(); });
// Calculate the azimuthal angle of a given pair of particles return pt;
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) { // Calculate the azimuthal angle phi of a vector of 4-vectors
return quantities.nu/1000.; inline auto phi(const std::vector<ROOT::Math::PxPyPzMVector>& momenta)
} {
inline double get_Q2_simu(inv_quant quantities) { std::vector<double> phi(momenta.size());
return quantities.Q2; // transform our raw tracker info into proper 4-momenta
} std::transform(momenta.begin(), momenta.end(), phi.begin(),
inline double get_x_simu(inv_quant quantities) { [](const auto& mom) { return mom.phi(); });
return quantities.x; 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;
} }
//for tracking, add later
//========================================================================================================= //=========================================================================================================
} // namespace util } // namespace util
#endif #endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment