diff --git a/README.md b/README.md
index 362982095dfd51819cd74cbc588635e7936d4458..5e22ee369a9db98f6743773ec97adbe6bb3bf13c 100644
--- a/README.md
+++ b/README.md
@@ -3,6 +3,31 @@ Physics Benchmarks for the EIC
+## Running Locally
+
+### Local development example
+
+Here we setup to use our local build of the `juggler` library.
+First set some environment variables.
+```
+export JUGGLER_INSTALL_PREFIX=$HOME/stow/juggler
+export JUGGLER_DETECTOR=athena # athena is the default
+export BEAMLINE_CONFIG=ip6 # ip6 is the default
+```
+
+
+```
+git@eicweb.phy.anl.gov:EIC/benchmarks/physics_benchmarks.git && cd physics_benchmarks
+git clone https://eicweb.phy.anl.gov/EIC/benchmarks/common_bench.git setup
+source setup/bin/env.sh && ./setup/bin/install_common.sh
+source .local/bin/env.sh && build_detector.sh
+mkdir_local_data_link sim_output
+mkdir -p results
+mkdir -p config
+
+```
+
+
## Common bench
See [common_bench](https://eicweb.phy.anl.gov/EIC/benchmarks/common_bench) for details.
diff --git a/benchmarks/dvcs/config.yml b/benchmarks/dvcs/config.yml
index 63d78140a8fee8660eb340b3fc9114a036136d82..ee77a0f3b6498bfc30da81457fbff4d7598a9d9e 100644
--- a/benchmarks/dvcs/config.yml
+++ b/benchmarks/dvcs/config.yml
@@ -5,8 +5,7 @@ dvcs:process:
- s3
needs: ["common:detector"]
script:
- - bash benchmarks/dvcs/dvcs.sh --data-init
- - bash benchmarks/dvcs/dvcs.sh
+ - bash benchmarks/dvcs/dvcs.sh --all
dvcs:results:
stage: collect
diff --git a/benchmarks/dvcs/dvcs.sh b/benchmarks/dvcs/dvcs.sh
index bcbfa38cabc7c20c0ff89070183e4d0ebf9c15f0..1dd20afb1aa1e8121ef9404c418e6c2bd0582e79 100644
--- a/benchmarks/dvcs/dvcs.sh
+++ b/benchmarks/dvcs/dvcs.sh
@@ -1,14 +1,22 @@
#!/bin/bash
function print_the_help {
- echo "USAGE: ${0} [--data-init] "
- echo " OPTIONS: "
+ echo "USAGE: ${0} [--rec] [--sim] [--analysis] [--all] "
+ echo " The default options are to run all steps (sim,rec,analysis) "
+ echo "OPTIONS: "
+ echo " --data-init download the input event data"
+ echo " --sim,-s Runs the Geant4 simulation"
+ echo " --rec,-r Run the juggler reconstruction"
+ echo " --analysis,-a Run the analysis scripts"
+ echo " --all (default) Do all steps. Argument is included so usage can convey intent."
exit
}
+DO_ALL=1
DATA_INIT=
-REC_ONLY=
-ANALYSIS_ONLY=
+DO_SIM=
+DO_REC=
+DO_ANALYSIS=
POSITIONAL=()
while [[ $# -gt 0 ]]
@@ -20,8 +28,33 @@ do
shift # past argument
print_the_help
;;
+ --all)
+ DO_ALL=2
+ if [[ ! "${DO_REC}${DO_SIM}${DO_ANALYSIS}" -eq "" ]] ; then
+ echo "Error: cannot use --all with other arguments." 1>&2
+ print_the_help
+ exit 1
+ fi
+ shift # past value
+ ;;
+ -s|--sim)
+ DO_SIM=1
+ DO_ALL=
+ shift # past value
+ ;;
--data-init)
DATA_INIT=1
+ DO_ALL=
+ shift # past value
+ ;;
+ -r|--rec)
+ DO_REC=1
+ DO_ALL=
+ shift # past value
+ ;;
+ -a|--analysis)
+ DO_ANALYSIS=1
+ DO_ALL=
shift # past value
;;
*) # unknown option
@@ -33,73 +66,86 @@ do
esac
done
set -- "${POSITIONAL[@]}" # restore positional parameters
-# these variables might not need exported.
-export FILE_NAME_TAG="dvcs"
-export JUGGLER_SIM_FILE="sim_${FILE_NAME_TAG}.root"
-export JUGGLER_REC_FILE="rec_${FILE_NAME_TAG}.root"
+# assuming something like .local/bin/env.sh has already been sourced.
+print_env.sh
-echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
-echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
-echo "FILE_NAME_TAG = ${FILE_NAME_TAG}"
+FILE_NAME_TAG="dvcs"
+DATA_URL="S3/eictest/ATHENA/EVGEN/DVCS/DVCS_10x100_2M/DVCS.1.hepmc"
+
+export JUGGLER_MC_FILE="${LOCAL_DATA_PATH}/mc_${FILE_NAME_TAG}.hepmc"
+export JUGGLER_SIM_FILE="${LOCAL_DATA_PATH}/sim_${FILE_NAME_TAG}.root"
+export JUGGLER_REC_FILE="${LOCAL_DATA_PATH}/rec_${FILE_NAME_TAG}.root"
+
+echo "FILE_NAME_TAG = ${FILE_NAME_TAG}"
+echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
+echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
-print_env.sh
## To run the reconstruction, we need the following global variables:
## - JUGGLER_INSTALL_PREFIX: Install prefix for Juggler (simu/recon)
## - JUGGLER_DETECTOR: the detector package we want to use for this benchmark
## - DETECTOR_PATH: full path to the detector definitions
-if [[ -n "${DATA_INIT}" ]] ; then
+## Step 1. Get the data
+if [[ -n "${DATA_INIT}" || -n "${DO_ALL}" ]] ; then
mc -C . config host add S3 https://dtn01.sdcc.bnl.gov:9000 $S3_ACCESS_KEY $S3_SECRET_KEY
- mc -C . cat --insecure S3/eictest/ATHENA/EVGEN/DVCS/DVCS_10x100_2M/DVCS.1.hepmc | head -n 1004 > "${LOCAL_DATA_PATH}/dvcs_test.hepmc"
+ mc -C . cat --insecure ${DATA_URL} | head -n 1004 > "${JUGGLER_MC_FILE}"
if [[ "$?" -ne "0" ]] ; then
echo "Failed to download hepmc file"
exit 1
fi
- exit
fi
-
-#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"
-
-
-## run geant4 simulations
-npsim --runType batch \
- --part.minimalKineticEnergy 1000*GeV \
- -v ERROR \
- --numberOfEvents ${JUGGLER_N_EVENTS} \
- --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
- --inputFiles "${LOCAL_DATA_PATH}/dvcs_test.hepmc" \
- --outputFile ${JUGGLER_SIM_FILE}
-if [[ "$?" -ne "0" ]] ; then
- echo "ERROR running npsim"
- exit 1
-fi
-
-# Need to figure out how to pass file name to juggler from the commandline
-xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
- gaudirun.py options/reconstruction.py
-if [[ "$?" -ne "0" ]] ; then
- echo "ERROR running juggler"
- exit 1
+### Step 2. Run the simulation (geant4)
+if [[ -n "${DO_SIM}" || -n "${DO_ALL}" ]] ; then
+ ## run geant4 simulations
+ npsim --runType batch \
+ --part.minimalKineticEnergy 1000*GeV \
+ -v ERROR \
+ --numberOfEvents ${JUGGLER_N_EVENTS} \
+ --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
+ --inputFiles "${JUGGLER_MC_FILE}" \
+ --outputFile ${JUGGLER_SIM_FILE}
+ if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running npsim"
+ exit 1
+ fi
fi
-mkdir -p results/dvcs
-rootls -t ${JUGGLER_SIM_FILE}
+### Step 3. Run the reconstruction (juggler)
+if [[ -n "${DO_REC}" || -n "${DO_ALL}" ]] ; then
+ xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
+ gaudirun.py options/reconstruction.py
+ if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running juggler"
+ exit 1
+ fi
-root -b -q "benchmarks/dvcs/scripts/dvcs_tests.cxx(\"${JUGGLER_REC_FILE}\")"
-if [[ "$?" -ne "0" ]] ; then
- echo "ERROR running root script"
- exit 1
+ root_filesize=$(stat --format=%s "${JUGGLER_REC_FILE}")
+ if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then
+ # file must be less than 10 MB to upload
+ if [[ "${root_filesize}" -lt "10000000" ]] ; then
+ cp ${JUGGLER_REC_FILE} results/.
+ fi
+ fi
fi
-# copy data if it is not too big
-if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then
-cp ${JUGGLER_REC_FILE} results/dvcs/.
-fi
+### Step 4. Run the analysis code
+if [[ -n "${DO_ANALYSIS}" || -n "${DO_ALL}" ]] ; then
+ echo "Running analysis scripts"
+ rootls -t ${JUGGLER_REC_FILE}
+ # Store all plots here (preferribly png and pdf files)
+ mkdir -p results/dvcs
+ # here you can add as many scripts as you want.
+ root -b -q "benchmarks/dvcs/scripts/dvcs_tests.cxx(\"${JUGGLER_REC_FILE}\")"
+ if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running root script"
+ exit 1
+ fi
+fi
diff --git a/benchmarks/dvcs/scripts/dvcs_tests.cxx b/benchmarks/dvcs/scripts/dvcs_tests.cxx
index d7d0df8dade34e1b77d98ad8a72d43302392f1cd..e801b4f78f5d4d1e8cca1b981af4a936a01816c2 100644
--- a/benchmarks/dvcs/scripts/dvcs_tests.cxx
+++ b/benchmarks/dvcs/scripts/dvcs_tests.cxx
@@ -1,19 +1,25 @@
-#include <ROOT/RDataFrame.hxx>
#include <cmath>
-#include "fmt/color.h"
-#include "fmt/core.h"
#include <iostream>
#include <string>
#include <vector>
+
+#include "ROOT/RDataFrame.hxx"
+
#include <nlohmann/json.hpp>
using json = nlohmann::json;
+R__LOAD_LIBRARY(libfmt.so)
+#include "fmt/core.h"
+#include "fmt/color.h"
+
+R__LOAD_LIBRARY(libeicd.so)
+R__LOAD_LIBRARY(libDD4pod.so)
+
#include "dd4pod/Geant4ParticleCollection.h"
#include "eicd/TrackParametersCollection.h"
#include "eicd/ClusterCollection.h"
-#include "eicd/ClusterData.h"
#include "eicd/ReconstructedParticleCollection.h"
-#include "eicd/ReconstructedParticleData.h"
+#include "eicd/BasicParticleCollection.h"
using ROOT::RDataFrame;
using namespace ROOT::VecOps;
@@ -58,7 +64,7 @@ auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
}
return result;
};
-auto dumfourvec = [](ROOT::VecOps::RVec<eic::ReconstructedParticleData> const& in) {
+auto recfourvec = [](ROOT::VecOps::RVec<eic::ReconstructedParticleData> const& in) {
std::vector<ROOT::Math::PxPyPzMVector> result;
ROOT::Math::PxPyPzMVector lv;
for (size_t i = 0; i < in.size(); ++i) {
@@ -106,34 +112,34 @@ void dvcs_tests(const char* fname = "rec_dvcs.root"){
auto df0 = df.Define("isThrown", "mcparticles2.genStatus == 1")
.Define("thrownParticles", "mcparticles2[isThrown]")
.Define("thrownP", fourvec, {"thrownParticles"})
- .Define("dumRec", dumfourvec, {"ReconstructedParticles"})
- .Define("dumNPart", "dumRec.size()")
+ .Define("recP", recfourvec, {"ReconstructedParticles"})
+ .Define("NPart", "recP.size()")
.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"})
.Define("eprime", eprime, {"thrownParticles"})
.Define("q", q_vec, {"eprime"})
.Define("Q2", "-1.0*(q.Dot(q))");
- auto h_n_dummy = df0.Histo1D({"h_n_dummy", "; h_n_dummy n", 10, 0, 10}, "dumNPart");
- auto h_Q2 = df0.Histo1D({"h_Q2", "; Q^{2} [GeV^{2}/c^{2}]", 100, 0, 30}, "Q2");
- auto n_Q2 = df0.Filter("Q2>1").Count();
+ auto h_n_dummy = df0.Histo1D({"h_n_part", "; h_n_part n", 10, 0, 10}, "NPart");
+ auto h_Q2 = df0.Histo1D({"h_Q2", "; Q^{2} [GeV^{2}/c^{2}]", 100, 0, 30}, "Q2");
+ auto n_Q2 = df0.Filter("Q2>1").Count();
+ auto n_tracks = df0.Mean("nTracks");
+
+ // ---------------------------
+ // Do evaluation
auto c = new TCanvas();
h_Q2->DrawCopy();
c->SaveAs("results/dvcs/Q2.png");
c->SaveAs("results/dvcs/Q2.pdf");
fmt::print("{} DVCS events\n",*n_Q2);
+ fmt::print("{} tracks per event\n",*n_tracks);
c = new TCanvas();
h_n_dummy->DrawCopy();
c->SaveAs("results/dvcs/n_dummy.png");
//c->SaveAs("results/dvcs/n_dummy.pdf");
-
}
diff --git a/options/reconstruction.py b/options/reconstruction.py
index ca61cbaec805b5ebc55e06f6272cb2464b399598..ac1173277ccac65a720903e46e01e6827cc8d88f 100644
--- a/options/reconstruction.py
+++ b/options/reconstruction.py
@@ -52,6 +52,8 @@ from Configurables import Jug__Base__InputCopier_dd4pod__CalorimeterHitCollectio
from Configurables import Jug__Base__InputCopier_dd4pod__TrackerHitCollection_dd4pod__TrackerHitCollection_ as TrkCopier
from Configurables import Jug__Base__InputCopier_dd4pod__PhotoMultiplierHitCollection_dd4pod__PhotoMultiplierHitCollection_ as PMTCopier
+from Configurables import Jug__Base__MC2DummyParticle as MC2DummyParticle
+
from Configurables import Jug__Digi__PhotoMultiplierDigi as PhotoMultiplierDigi
from Configurables import Jug__Digi__CalorimeterHitDigi as CalHitDigi
from Configurables import Jug__Digi__UFSDTrackerDigi as TrackerDigi
@@ -66,6 +68,7 @@ from Configurables import Jug__Reco__TrackParamClusterInit as TrackParamClusterI
from Configurables import Jug__Reco__TrackParamVertexClusterInit as TrackParamVertexClusterInit
from Configurables import Jug__Reco__TrackFindingAlgorithm as TrackFindingAlgorithm
from Configurables import Jug__Reco__ParticlesFromTrackFit as ParticlesFromTrackFit
+from Configurables import Jug__Reco__ParticlesWithTruthPID as ParticlesWithTruthPID
from Configurables import Jug__Reco__CalorimeterHitReco as CalHitReco
from Configurables import Jug__Reco__CalorimeterHitsMerger as CalHitsMerger
@@ -114,15 +117,12 @@ copier = MCCopier("MCCopier",
outputCollection="mcparticles2")
algorithms.append(copier)
-trkcopier = TrkCopier("TrkCopier",
- inputCollection="TrackerBarrelHits",
- outputCollection="TrackerBarrelHits2")
-algorithms.append(trkcopier)
-
-pmtcopier = PMTCopier("PMTCopier",
- inputCollection="DRICHHits",
- outputCollection="DRICHHits2")
-algorithms.append(pmtcopier)
+# Generated particles
+dummy = MC2DummyParticle("dummy",
+ inputCollection="mcparticles",
+ outputCollection="GeneratedParticles",
+ smearing=0)
+algorithms.append(dummy)
# Crystal Endcap Ecal
ce_ecal_daq = dict(
@@ -161,6 +161,7 @@ ce_ecal_clreco = RecoCoG("ce_ecal_clreco",
inputHitCollection=ce_ecal_cl.inputHitCollection,
inputProtoClusterCollection=ce_ecal_cl.outputProtoClusterCollection,
outputClusterCollection="EcalEndcapNClusters",
+ outputInfoCollection="EcalEndcapNClusterInfo",
samplingFraction=0.998, # this accounts for a small fraction of leakage
logWeightBase=4.6)
algorithms.append(ce_ecal_clreco)
@@ -188,7 +189,7 @@ algorithms.append(ci_ecal_reco)
# merge hits in different layer (projection to local x-y plane)
ci_ecal_merger = CalHitsMerger("ci_ecal_merger",
inputHitCollection=ci_ecal_reco.outputHitCollection,
- outputHitCollection="EcalEndcapPRecHitsXY",
+ outputHitCollection="EcalEndcapPRecMergedHits",
fields=["layer", "slice"],
fieldRefNumbers=[1, 0],
readoutClass="EcalEndcapPHits")
@@ -280,7 +281,7 @@ algorithms.append(scfi_barrel_reco)
# merge hits in different layer (projection to local x-y plane)
scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger",
inputHitCollection=scfi_barrel_reco.outputHitCollection,
- outputHitCollection="EcalBarrelScFiGridReco",
+ outputHitCollection="EcalBarrelScFiMergedHits",
fields=["fiber"],
fieldRefNumbers=[1],
readoutClass="EcalBarrelScFiHits")
@@ -328,7 +329,7 @@ algorithms.append(cb_hcal_reco)
cb_hcal_merger = CalHitsMerger("cb_hcal_merger",
inputHitCollection=cb_hcal_reco.outputHitCollection,
- outputHitCollection="HcalBarrelRecHitsXY",
+ outputHitCollection="HcalBarrelMergedHits",
readoutClass="HcalBarrelHits",
fields=["layer", "slice"],
fieldRefNumbers=[1, 0])
@@ -373,7 +374,7 @@ algorithms.append(ci_hcal_reco)
ci_hcal_merger = CalHitsMerger("ci_hcal_merger",
inputHitCollection=ci_hcal_reco.outputHitCollection,
- outputHitCollection="HcalEndcapPRecHitsXY",
+ outputHitCollection="HcalEndcapPMergedHits",
readoutClass="HcalEndcapPHits",
fields=["layer", "slice"],
fieldRefNumbers=[1, 0])
@@ -418,7 +419,7 @@ algorithms.append(ce_hcal_reco)
ce_hcal_merger = CalHitsMerger("ce_hcal_merger",
inputHitCollection=ce_hcal_reco.outputHitCollection,
- outputHitCollection="HcalEndcapNRecHitsXY",
+ outputHitCollection="HcalEndcapNMergedHits",
readoutClass="HcalEndcapNHits",
fields=["layer", "slice"],
fieldRefNumbers=[1, 0])
@@ -532,10 +533,17 @@ algorithms.append(trk_find_alg)
parts_from_fit = ParticlesFromTrackFit("parts_from_fit",
inputTrajectories = trk_find_alg.outputTrajectories,
- outputParticles = "ReconstructedParticles",
+ outputParticles = "outputParticles",
outputTrackParameters = "outputTrackParameters")
algorithms.append(parts_from_fit)
+# Event building
+parts_with_truth_pid = ParticlesWithTruthPID("parts_with_truth_pid",
+ inputMCParticles = "mcparticles",
+ inputTrackParameters = parts_from_fit.outputTrackParameters,
+ outputParticles = "ReconstructedParticles")
+algorithms.append(parts_with_truth_pid)
+
# DRICH
drich_digi = PhotoMultiplierDigi("drich_digi",
inputHitCollection="DRICHHits",
@@ -576,11 +584,11 @@ algorithms.append(mrich_reco)
podout = PodioOutput("out", filename=output_rec)
podout.outputCommands = [
"keep *",
- "drop *Digi",
- "keep *Reco*",
- "keep *ClusterHits",
- "keep *Clusters",
+ "drop *Hits",
"keep *Layers",
+ "keep *Clusters",
+ "drop *ProtoClusters",
+ "drop outputParticles",
"drop InitTrackParams",
] + [ "drop " + c for c in sim_coll]
algorithms.append(podout)