diff --git a/benchmarks/tracking/central_pions.sh b/benchmarks/tracking/central_pions.sh
new file mode 100644
index 0000000000000000000000000000000000000000..214913fb215ac319a5ad6a25138f5a4b5768f247
--- /dev/null
+++ b/benchmarks/tracking/central_pions.sh
@@ -0,0 +1,127 @@
+#!/bin/bash
+
+function print_the_help {
+ echo "USAGE: ${0} [--rec-only] "
+ echo " OPTIONS: "
+ echo " --rec-only only run gaudi reconstruction part"
+ exit
+}
+
+REC_ONLY=
+ANALYSIS_ONLY=
+
+POSITIONAL=()
+while [[ $# -gt 0 ]]
+do
+ key="$1"
+
+ case $key in
+ -h|--help)
+ shift # past argument
+ print_the_help
+ ;;
+ --rec-only)
+ REC_ONLY=1
+ shift # past value
+ ;;
+ --ana-only)
+ ANALYSIS_ONLY=1
+ shift # past value
+ ;;
+ *) # unknown option
+ #POSITIONAL+=("$1") # save it in an array for later
+ echo "unknown option $1"
+ print_the_help
+ shift # past argument
+ ;;
+ esac
+done
+set -- "${POSITIONAL[@]}" # restore positional parameters
+
+
+print_env.sh
+
+## To run the reconstruction, we need the following global variables:
+## - JUGGLER_INSTALL_PREFIX: Install prefix for Juggler (simu/recon)
+## - JUGGLER_DETECTOR: the detector package we want to use for this benchmark
+## - JUGGLER_DETECTOR_VERSION: the detector package we want to use for this benchmark
+## - DETECTOR_PATH: full path to the detector definitions
+##
+## You can ready options/env.sh for more in-depth explanations of the variables
+## and how they can be controlled.
+export DETECTOR_PATH=${DETECTOR_PATH}
+
+if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then
+ export JUGGLER_N_EVENTS=100
+fi
+
+
+export JUGGLER_FILE_NAME_TAG="central_pions"
+export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
+
+export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root"
+export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
+
+echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
+echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
+
+
+if [[ -z "${REC_ONLY}" && -z "${ANALYSIS_ONLY}" ]] ;
+then
+
+ ## generate the input events
+ root -b -q "benchmarks/tracking/scripts/gen_central_pions.cxx(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+ if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running script"
+ exit 1
+ fi
+
+ echo "Running geant4 simulation"
+ ## run geant4 simulations
+ npsim --runType batch \
+ --part.minimalKineticEnergy 1000*GeV \
+ -v WARNING \
+ --numberOfEvents ${JUGGLER_N_EVENTS} \
+ --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
+ --inputFiles ${JUGGLER_FILE_NAME_TAG}.hepmc \
+ --outputFile ${JUGGLER_SIM_FILE}
+ if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running script"
+ exit 1
+ fi
+fi
+
+rootls -t ${JUGGLER_SIM_FILE}
+
+if [[ -z "${ANALYSIS_ONLY}" ]] ;
+then
+ # Need to figure out how to pass file name to juggler from the commandline
+ xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv gaudirun.py benchmarks/tracking/options/tracker_reconstruction.py
+ if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running juggler"
+ exit 1
+ fi
+fi
+
+mkdir -p results/tracking
+
+root -b -q "benchmarks/tracking/scripts/rec_central_pions.cxx(\"${JUGGLER_REC_FILE}\")"
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running root script"
+ exit 1
+fi
+
+root -b -q "benchmarks/tracking/scripts/hits_central_pions.cxx(\"${JUGGLER_SIM_FILE}\")"
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running root script"
+ exit 1
+fi
+
+root_filesize=$(stat --format=%s "${JUGGLER_REC_FILE}")
+if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then
+ # file must be less than 10 MB to upload
+ if [[ "${root_filesize}" -lt "10000000" ]] ; then
+ cp ${JUGGLER_REC_FILE} results/.
+ fi
+fi
+
diff --git a/benchmarks/tracking/config.yml b/benchmarks/tracking/config.yml
index b707a2eb2143c42bbe95e26633393a0ab9394e29..070ff3a842e4dc4a4d470a46e2bd2b96ffec046d 100644
--- a/benchmarks/tracking/config.yml
+++ b/benchmarks/tracking/config.yml
@@ -6,10 +6,24 @@ tracking_central_electrons:
- bash benchmarks/tracking/central_electrons.sh
#allow_failure: true
+tracking_central_pions:
+ extends: .rec_benchmark
+ stage: run
+ timeout: 24 hours
+ script:
+ - bash benchmarks/tracking/central_pions.sh
+ #allow_failure: true
+
tracking_truth_init_electrons:
extends: .rec_benchmark
stage: run
timeout: 24 hours
script:
- - python benchmarks/tracking/run_tracking_benchmarks.py --etamin=-3 --etamax=3 -n 100
+ - python benchmarks/tracking/run_tracking_benchmarks.py --nametag=electron --particle=electron --etamin=-3 --etamax=3 -n 100
+tracking_truth_init_pions:
+ extends: .rec_benchmark
+ stage: run
+ timeout: 24 hours
+ script:
+ - python benchmarks/tracking/run_tracking_benchmarks.py --nametag=pion --particle=pion+ --etamin=-3 --etamax=3 -n 100
diff --git a/benchmarks/tracking/run_tracking_benchmarks.py b/benchmarks/tracking/run_tracking_benchmarks.py
index 97078593eb4708eabbd9fc88df0df417b4681692..4929f0e6daa6b1b26ea05017a0b37e31c5d3a894 100755
--- a/benchmarks/tracking/run_tracking_benchmarks.py
+++ b/benchmarks/tracking/run_tracking_benchmarks.py
@@ -94,7 +94,7 @@ if 'ana' in procs:
ana_cmd = ['python', analysis_script, rec_file,
'--mc-collection', 'mcparticles2',
'--tracking-collection', 'outputTrackParameters',
- '-o', 'results']
+ '-o', 'results', '-t', args.nametag]
return_code = subprocess.run(ana_cmd).returncode
if return_code is not None and return_code < 0:
print('ERROR running analysis ({})!'.format(ana))
diff --git a/benchmarks/tracking/scripts/gen_central_pions.cxx b/benchmarks/tracking/scripts/gen_central_pions.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..abf691b5a984ad23758d94d900b9fa74e438b874
--- /dev/null
+++ b/benchmarks/tracking/scripts/gen_central_pions.cxx
@@ -0,0 +1,83 @@
+#include "HepMC3/GenEvent.h"
+#include "HepMC3/ReaderAscii.h"
+#include "HepMC3/WriterAscii.h"
+#include "HepMC3/Print.h"
+
+#include <iostream>
+#include<random>
+#include<cmath>
+#include <math.h>
+#include <TMath.h>
+
+using namespace HepMC3;
+
+/** Generate pions in the central region.
+ * This is for testing detectors in the "barrel" region.
+ */
+void gen_central_pions(int n_events = 100,
+ const char* out_fname = "central_pions.hepmc")
+{
+ double cos_theta_min = std::cos( 80.0*(M_PI/180.0));
+ double cos_theta_max = std::cos(100.0*(M_PI/180.0));
+
+ WriterAscii hepmc_output(out_fname);
+ int events_parsed = 0;
+ GenEvent evt(Units::GEV, Units::MM);
+
+ // Random number generator
+ TRandom *r1 = new TRandom();
+
+ for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
+ // FourVector(px,py,pz,e,pdgid,status)
+ // type 4 is beam
+ // pdgid 11 - pion
+ // pdgid 111 - pi0
+ // pdgid 2212 - proton
+ GenParticlePtr p1 =
+ std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
+ GenParticlePtr p2 = std::make_shared<GenParticle>(
+ FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
+
+ // Define momentum
+ Double_t p = r1->Uniform(1.0, 10.0);
+ Double_t phi = r1->Uniform(0.0, 2.0 * M_PI);
+ Double_t costh = r1->Uniform(cos_theta_min, cos_theta_max);
+ Double_t th = std::acos(costh);
+ Double_t px = p * std::cos(phi) * std::sin(th);
+ Double_t py = p * std::sin(phi) * std::sin(th);
+ Double_t pz = p * std::cos(th);
+ // Generates random vectors, uniformly distributed over the surface of a
+ // sphere of given radius, in this case momentum.
+ // r1->Sphere(px, py, pz, p);
+
+ //std::cout << std::sqrt(px*px + py*py + pz*pz) - p << " is zero? \n";
+
+ // type 1 is final state
+ // pdgid 211 - pion 139.57039 MeV/c^2
+ GenParticlePtr p3 = std::make_shared<GenParticle>(
+ FourVector(
+ px, py, pz,
+ sqrt(p*p + (0.13957039 * 0.13957039))),
+ 211, 1);
+
+ GenVertexPtr v1 = std::make_shared<GenVertex>();
+ v1->add_particle_in(p1);
+ v1->add_particle_in(p2);
+
+ v1->add_particle_out(p3);
+ evt.add_vertex(v1);
+
+ if (events_parsed == 0) {
+ std::cout << "First event: " << std::endl;
+ Print::listing(evt);
+ }
+
+ hepmc_output.write_event(evt);
+ if (events_parsed % 10000 == 0) {
+ std::cout << "Event: " << events_parsed << std::endl;
+ }
+ evt.clear();
+ }
+ hepmc_output.close();
+ std::cout << "Events parsed and written: " << events_parsed << std::endl;
+}
diff --git a/benchmarks/tracking/scripts/hits_central_pions.cxx b/benchmarks/tracking/scripts/hits_central_pions.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..2e97e0dfd3bb7659a3cd4b1dd02e86275bf3a5e9
--- /dev/null
+++ b/benchmarks/tracking/scripts/hits_central_pions.cxx
@@ -0,0 +1,190 @@
+#include "ROOT/RDataFrame.hxx"
+#include "TCanvas.h"
+#include "TLegend.h"
+#include "TH1D.h"
+#include "TProfile.h"
+
+#include <iostream>
+
+R__LOAD_LIBRARY(libeicd.so)
+R__LOAD_LIBRARY(libDD4pod.so)
+#include "dd4pod/TrackerHitCollection.h"
+#include "dd4pod/TrackerHitData.h"
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "eicd/TrackParametersCollection.h"
+#include "eicd/ClusterCollection.h"
+#include "eicd/ClusterData.h"
+#include "eicd/TrackerHitCollection.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;
+};
+
+auto delta_p_over_p = [](const std::vector<double>& tracks, const std::vector<double>& thrown) {
+ std::vector<double> res;
+ for (const auto& p1 : thrown) {
+ for (const auto& p2 : tracks) {
+ res.push_back((p1 - p2)/p1);
+ }
+ }
+ return res;
+};
+
+int hits_central_pions(const char* fname = "sim_central_pions.root")
+{
+
+ ROOT::EnableImplicitMT();
+ ROOT::RDataFrame df("events", fname);
+
+ auto df0 = df.Define("isThrown", "mcparticles.genStatus == 1")
+ .Define("thrownParticles", "mcparticles[isThrown]")
+ .Define("thrownP", fourvec, {"thrownParticles"})
+ .Define("p_thrown", momentum, {"thrownP"})
+ .Define("theta_thrown", theta, {"thrownP"})
+ .Define("theta0", "theta_thrown[0]")
+ //.Define("nTracks", "outputTrackParameters.size()")
+ //.Define("p_track", p_track, {"outputTrackParameters"})
+ //.Define("p_track1", p_track, {"outputTrackParameters1"})
+ //.Define("p_track2", p_track, {"outputTrackParameters2"})
+ //.Define("delta_p0",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("delta_p_over_p0",delta_p_over_p, {"p_track", "p_thrown"})
+ //.Define("delta_p_over_p1",delta_p_over_p, {"p_track1", "p_thrown"})
+ //.Define("delta_p_over_p2",delta_p_over_p, {"p_track2", "p_thrown"})
+ //.Define("N_VtxBarrelHits",[](std::vector<eic::TrackerHitData> hits) { return hits.size();},{"VertexBarrelRecHits"})
+ .Define("N_BarrelHits", [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"TrackerBarrelHits"})
+ .Define("N_EndcapHits", [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"TrackerEndcapHits"})
+ ;
+
+ auto hBarrel_x_vs_y = df0.Histo2D({"hBarrel_x_vs_y", "; x ; y ", 100, -900, 900,100, -900, 900 }, "TrackerBarrelHits.position.x", "TrackerBarrelHits.position.y");
+
+ auto hBarrel_N_vs_theta = df0.Histo1D({"hBarrel_N_vs_theta", "; #theta [deg.]", 20, 0, 180 }, "theta0", "N_BarrelHits");
+ auto hEndcap_N_vs_theta = df0.Histo1D({"hEndcap_N_vs_theta", "; #theta [deg.]", 20, 0, 180 }, "theta0", "N_EndcapHits");
+ //auto hVtxBarrel_N_vs_theta = df0.Histo1D({"hVtxBarrel_N_vs_theta", "; #theta [deg.]", 20, 0, 180 }, "theta0", "N_VtxBarrelHits");
+
+ auto hBarrel_Nhits = df0.Histo1D({"hBarrel_Nhits", "; #theta [deg.]", 20, 0, 20 }, "N_BarrelHits");
+ auto hEndcap_Nhits = df0.Histo1D({"hEndcap_Nhits", "; #theta [deg.]", 20, 0, 20 }, "N_EndcapHits");
+ //auto hVtxBarrel_Nhits = df0.Histo1D({"hVtxBarrel_Nhits", "; #theta [deg.]", 20, 0, 20 }, "N_VtxBarrelHits");
+
+ auto hBarrel_Ntheta = df0.Histo1D({"hBarrel_Ntheta", "; #theta [deg.]", 20, 0, 180 }, "theta0");
+ auto hEndcap_Ntheta = df0.Histo1D({"hEndcap_Ntheta", "; #theta [deg.]", 20, 0, 180 }, "theta0");
+ //auto hVtxBarrel_Ntheta = df0.Histo1D({"hVtxBarrel_Ntheta", "; #theta [deg.]", 20, 0, 180 }, "theta0");
+
+ auto c = new TCanvas();
+ auto hs = new THStack("n_hits","; #theta ");
+ auto h1 = (TH1D*) hBarrel_N_vs_theta->Clone();
+ auto h2 = (TH1D*) hBarrel_Ntheta->Clone();
+ h1->Divide(h2);
+ hs->Add(h1);
+ h1 = (TH1D*) hEndcap_N_vs_theta->Clone();
+ h2 = (TH1D*) hEndcap_Ntheta->Clone();
+ h1->Divide(h2);
+ h1->SetLineColor(2);
+ hs->Add(h1);
+ //h1 = (TH1D*) hVtxBarrel_vs_theta->Clone();
+ //h1->SetLineColor(4);
+ //h1->SetFillStyle(3001);
+ //h1->SetFillColor(4);
+ //hs->Add(h1);
+ hs->Draw("nostack, hist");
+ c->BuildLegend();
+ c->SaveAs("results/tracking/hits_central_pions_n_hits_vs_theta.png");
+ c->SaveAs("results/tracking/hits_central_pions_n_hits_vs_theta.pdf");
+
+ c = new TCanvas();
+ hs = new THStack("theta","; #theta ");
+ h1 = (TH1D*) hBarrel_N_vs_theta->Clone();
+ h2 = (TH1D*) hBarrel_Ntheta->Clone();
+ //h1->Divide(h2);
+ hs->Add(h2);
+ h1 = (TH1D*) hEndcap_N_vs_theta->Clone();
+ h2 = (TH1D*) hEndcap_Ntheta->Clone();
+ //h1->Divide(h2);
+ h1->SetLineColor(2);
+ h2->SetLineColor(2);
+ hs->Add(h2);
+ //h1 = (TH1D*) hVtxBarrel_vs_theta->Clone();
+ //h1->SetLineColor(4);
+ //h1->SetFillStyle(3001);
+ //h1->SetFillColor(4);
+ //hs->Add(h1);
+ hs->Draw("nostack hist");
+ c->BuildLegend();
+ c->SaveAs("results/tracking/hits_central_pions_theta.png");
+ c->SaveAs("results/tracking/hits_central_pions_theta.pdf");
+
+ c = new TCanvas();
+ hs = new THStack("hits","; hits ");
+ h1 = (TH1D*) hBarrel_Nhits->Clone();
+ hs->Add(h1);
+ h1 = (TH1D*) hEndcap_Nhits->Clone();
+ h1->SetLineColor(2);
+ h2->SetLineColor(2);
+ hs->Add(h2);
+ //h1 = (TH1D*) hVtxBarrel_Nhits->Clone();
+ //h1->SetLineColor(4);
+ //h1->SetFillStyle(3001);
+ //h1->SetFillColor(4);
+ //hs->Add(h1);
+ //hs->Draw("nostack hist");
+ c->BuildLegend();
+ c->SaveAs("results/tracking/hits_central_pions_nhits.png");
+ c->SaveAs("results/tracking/hits_central_pions_nhits.pdf");
+
+ c = new TCanvas();
+ hBarrel_x_vs_y->DrawCopy("colz");
+ c->SaveAs("results/tracking/hits_central_pions_xy.png");
+ c->SaveAs("results/tracking/hits_central_pions_xy.pdf");
+ return 0;
+}
diff --git a/benchmarks/tracking/scripts/rec_central_electrons.cxx b/benchmarks/tracking/scripts/rec_central_electrons.cxx
index 962054a2a8b4d6a18479dbc01e76d1fb253fb380..2dc0077cc9e71436308c541bbf55729066dab09b 100644
--- a/benchmarks/tracking/scripts/rec_central_electrons.cxx
+++ b/benchmarks/tracking/scripts/rec_central_electrons.cxx
@@ -37,7 +37,7 @@ std::vector<float> pt (std::vector<dd4pod::Geant4ParticleData> const& in){
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());
+ result.push_back(in[i].P());
}
return result;
};
diff --git a/benchmarks/tracking/scripts/rec_central_pions.cxx b/benchmarks/tracking/scripts/rec_central_pions.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..24be4f35dabb6577658f8708c4f23e2a5d36988f
--- /dev/null
+++ b/benchmarks/tracking/scripts/rec_central_pions.cxx
@@ -0,0 +1,238 @@
+#include "ROOT/RDataFrame.hxx"
+#include "TCanvas.h"
+#include "TLegend.h"
+#include "TH1D.h"
+#include "TProfile.h"
+
+#include <iostream>
+
+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/TrackerHitCollection.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].P());
+ }
+ return result;
+};
+auto theta = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
+ std::vector<double> result;
+ for (size_t i = 0; i < in.size(); ++i) {
+ result.push_back(in[i].Theta()*180/M_PI);
+ }
+ return result;
+};
+auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
+ std::vector<ROOT::Math::PxPyPzMVector> result;
+ ROOT::Math::PxPyPzMVector lv;
+ for (size_t i = 0; i < in.size(); ++i) {
+ lv.SetCoordinates(in[i].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;
+};
+
+auto delta_p_over_p = [](const std::vector<double>& tracks, const std::vector<double>& thrown) {
+ std::vector<double> res;
+ for (const auto& p1 : thrown) {
+ for (const auto& p2 : tracks) {
+ res.push_back((p1 - p2)/p1);
+ }
+ }
+ return res;
+};
+
+int rec_central_pions(const char* fname = "topside/rec_central_pions.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("theta_thrown", theta, {"thrownP"})
+ .Define("theta0", "theta_thrown[0]")
+ .Define("nTracks", "outputTrackParameters.size()")
+ .Define("p_track", p_track, {"outputTrackParameters"})
+ .Define("p_track1", p_track, {"outputTrackParameters1"})
+ //.Define("p_track2", p_track, {"outputTrackParameters2"})
+ .Define("delta_p0",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("delta_p_over_p0",delta_p_over_p, {"p_track", "p_thrown"})
+ .Define("delta_p_over_p1",delta_p_over_p, {"p_track1", "p_thrown"})
+ //.Define("delta_p_over_p2",delta_p_over_p, {"p_track2", "p_thrown"})
+ //.Define("N_VtxBarrelHits",[](std::vector<eic::TrackerHitData> hits) { return hits.size();},{"VertexBarrelRecHits"})
+ .Define("N_SiBarrelHits", [](std::vector<eic::TrackerHitData> hits) { return hits.size();}, {"TrackerBarrelRecHits"})
+ .Define("N_SiEndcapHits", [](std::vector<eic::TrackerHitData> hits) { return hits.size();}, {"TrackerEndcapRecHits"})
+ ;
+
+ 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_p0 = df0.Histo1D({"h_delta_p0", "Truth Track Init; GeV/c ", 100, -10, 10}, "delta_p0");
+ auto h_delta_p1 = df0.Histo1D({"h_delta_p1", "Ecal Cluster Init; GeV/c ", 100, -10, 10}, "delta_p1");
+ //auto h_delta_p2 = df0.Histo1D({"h_delta_p2", "Ecal Cluster, innner vtx hit Init; GeV/c ", 100, -10, 10}, "delta_p2");
+
+ auto h_delta_p0_over_p = df0.Histo1D({"h_delta_p0_over_p", "Truth Track Init; delta p/p ", 100, -0.1, 0.1}, "delta_p_over_p0");
+ auto h_delta_p1_over_p = df0.Histo1D({"h_delta_p1_over_p", "Ecal Cluster Init; delta p/p ", 100, -0.1, 0.1}, "delta_p_over_p1");
+ //auto h_delta_p2_over_p = df0.Histo1D({"h_delta_p2_over_p", "Ecal Cluster, innner vtx hit Init; delta p/p ", 100, -0.1, 0.1}, "delta_p_over_p2");
+
+ auto hSiBarrel_N_vs_theta = df0.Histo1D({"hSiBarrel_N_vs_theta", "; #theta [deg.]", 20, 0, 180 }, "theta0", "N_SiBarrelHits");
+ auto hSiEndcap_N_vs_theta = df0.Histo1D({"hSiEndcap_N_vs_theta", "; #theta [deg.]", 20, 0, 180 }, "theta0", "N_SiEndcapHits");
+ //auto hVtxBarrel_N_vs_theta = df0.Histo1D({"hVtxBarrel_N_vs_theta", "; #theta [deg.]", 20, 0, 180 }, "theta0", "N_VtxBarrelHits");
+
+ auto hSiBarrel_Nhits = df0.Histo1D({"hSiBarrel_Nhits", "; #theta [deg.]", 20, 0, 20 }, "N_SiBarrelHits");
+ auto hSiEndcap_Nhits = df0.Histo1D({"hSiEndcap_Nhits", "; #theta [deg.]", 20, 0, 20 }, "N_SiEndcapHits");
+ //auto hVtxBarrel_Nhits = df0.Histo1D({"hVtxBarrel_Nhits", "; #theta [deg.]", 20, 0, 20 }, "N_VtxBarrelHits");
+
+ auto hSiBarrel_Ntheta = df0.Histo1D({"hSiBarrel_Ntheta", "; #theta [deg.]", 20, 0, 180 }, "theta0");
+ auto hSiEndcap_Ntheta = df0.Histo1D({"hSiEndcap_Ntheta", "; #theta [deg.]", 20, 0, 180 }, "theta0");
+ //auto hVtxBarrel_Ntheta = df0.Histo1D({"hVtxBarrel_Ntheta", "; #theta [deg.]", 20, 0, 180 }, "theta0");
+
+ auto c = new TCanvas();
+
+ h_nTracks->DrawCopy();
+ c->SaveAs("results/tracking/rec_central_pions_nTracks.png");
+ c->SaveAs("results/tracking/rec_central_pions_nTracks.pdf");
+
+ h_pTracks->DrawCopy();
+ c->SaveAs("results/tracking/rec_central_pions_pTracks.png");
+ c->SaveAs("results/tracking/rec_central_pions_pTracks.pdf");
+
+ c = new TCanvas();
+ THStack * hs = new THStack("hs_delta_p","; GeV/c ");
+ TH1D* h1 = (TH1D*) h_delta_p0->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->BuildLegend();
+ c->SaveAs("results/tracking/rec_central_pions_delta_p.png");
+ c->SaveAs("results/tracking/rec_central_pions_delta_p.pdf");
+
+ c = new TCanvas();
+ hs = new THStack("hs_delta_p_over_p","; delta p/p ");
+ h1 = (TH1D*) h_delta_p0_over_p->Clone();
+ hs->Add(h1);
+ h1 = (TH1D*) h_delta_p1_over_p->Clone();
+ h1->SetLineColor(2);
+ hs->Add(h1);
+ //h1 = (TH1D*) h_delta_p2_over_p->Clone();
+ //h1->SetLineColor(4);
+ //h1->SetFillStyle(3001);
+ //h1->SetFillColor(4);
+ //hs->Add(h1);
+ hs->Draw("nostack");
+ c->BuildLegend();
+ c->SaveAs("results/tracking/rec_central_pions_delta_p_over_p.png");
+ c->SaveAs("results/tracking/rec_central_pions_delta_p_over_p.pdf");
+
+ c = new TCanvas();
+ hs = new THStack("n_hits","; #theta ");
+ h1 = (TH1D*) hSiBarrel_N_vs_theta->Clone();
+ auto h2 = (TH1D*) hSiBarrel_Ntheta->Clone();
+ h1->Divide(h2);
+ hs->Add(h1);
+ h1 = (TH1D*) hSiEndcap_N_vs_theta->Clone();
+ h2 = (TH1D*) hSiEndcap_Ntheta->Clone();
+ h1->Divide(h2);
+ h1->SetLineColor(2);
+ hs->Add(h1);
+ //h1 = (TH1D*) hVtxBarrel_vs_theta->Clone();
+ //h1->SetLineColor(4);
+ //h1->SetFillStyle(3001);
+ //h1->SetFillColor(4);
+ //hs->Add(h1);
+ hs->Draw("nostack, hist");
+ c->BuildLegend();
+ c->SaveAs("results/tracking/rec_central_pions_n_hits_vs_theta.png");
+ c->SaveAs("results/tracking/rec_central_pions_n_hits_vs_theta.pdf");
+
+ c = new TCanvas();
+ hs = new THStack("theta","; #theta ");
+ h1 = (TH1D*) hSiBarrel_N_vs_theta->Clone();
+ h2 = (TH1D*) hSiBarrel_Ntheta->Clone();
+ //h1->Divide(h2);
+ hs->Add(h2);
+ h1 = (TH1D*) hSiEndcap_N_vs_theta->Clone();
+ h2 = (TH1D*) hSiEndcap_Ntheta->Clone();
+ //h1->Divide(h2);
+ h1->SetLineColor(2);
+ h2->SetLineColor(2);
+ hs->Add(h2);
+ //h1 = (TH1D*) hVtxBarrel_vs_theta->Clone();
+ //h1->SetLineColor(4);
+ //h1->SetFillStyle(3001);
+ //h1->SetFillColor(4);
+ //hs->Add(h1);
+ hs->Draw("nostack hist");
+ c->BuildLegend();
+ c->SaveAs("results/tracking/rec_central_pions_theta.png");
+ c->SaveAs("results/tracking/rec_central_pions_theta.pdf");
+
+ c = new TCanvas();
+ hs = new THStack("hits","; hits ");
+ h1 = (TH1D*) hSiBarrel_Nhits->Clone();
+ hs->Add(h1);
+ h1 = (TH1D*) hSiEndcap_Nhits->Clone();
+ h1->SetLineColor(2);
+ h2->SetLineColor(2);
+ hs->Add(h2);
+ //h1 = (TH1D*) hVtxBarrel_Nhits->Clone();
+ //h1->SetLineColor(4);
+ //h1->SetFillStyle(3001);
+ //h1->SetFillColor(4);
+ //hs->Add(h1);
+ //hs->Draw("nostack hist");
+ c->BuildLegend();
+ c->SaveAs("results/tracking/rec_central_pions_nhits.png");
+ c->SaveAs("results/tracking/rec_central_pions_nhits.pdf");
+
+
+ return 0;
+}
diff --git a/benchmarks/tracking/scripts/tracking_performance.py b/benchmarks/tracking/scripts/tracking_performance.py
index 5f7c414c651b4c63882cdb0848d42e60506db393..9c34738339696cdf528e140365969ba1f1f6c9d5 100644
--- a/benchmarks/tracking/scripts/tracking_performance.py
+++ b/benchmarks/tracking/scripts/tracking_performance.py
@@ -110,6 +110,7 @@ if __name__ == '__main__':
parser = argparse.ArgumentParser()
parser.add_argument('rec_file', help='Path to reconstruction output file.')
parser.add_argument('-o', dest='outdir', default='.', help='Output directory.')
+ parser.add_argument('-t', '--nametag', type=str, default='tracking', help='Name tag for output files.')
parser.add_argument('-m', '--macros', dest='macros', default='rootlogon.C',
help='Macro files to be loaded by root, separated by \",\".')
parser.add_argument('--mc-collection', dest='mc', default='mcparticles',
@@ -133,7 +134,7 @@ if __name__ == '__main__':
load_root_macros(args.macros)
rdf_rec = ROOT.RDataFrame('events', args.rec_file)
- dfm = thrown_particles_figure(rdf_rec, save=os.path.join(args.outdir, 'thrown_particles.png'), mcbranch=args.mc)
+ dfm = thrown_particles_figure(rdf_rec, save=os.path.join(args.outdir, '{}_thrown_particles.png'.format(args.nametag)), mcbranch=args.mc)
df = flatten_collection(rdf_rec, args.coll)
df.rename(columns={c: c.replace(args.coll + '.', '') for c in df.columns}, inplace=True)
@@ -197,7 +198,6 @@ if __name__ == '__main__':
nbins = hbins.shape[0] - 1
ax.set_ylabel('Normalized Counts / {:d} Bins'.format(nbins), fontsize=20)
ax.set_xlabel(r'$d\theta$ (rad)', fontsize=20)
- fig.savefig('test.png')
# phi resolution
ax = axs.flat[3]
@@ -212,5 +212,5 @@ if __name__ == '__main__':
ax.set_xlabel(r'$d\phi$ (rad)', fontsize=20)
fig.text(0.5, 0.95, 'Barrel Tracker Benchmark (Truth Init.)', fontsize=22, ha='center')
- fig.savefig(os.path.join(args.outdir, 'tracking_performance.png'))
+ fig.savefig(os.path.join(args.outdir, '{}_performance.png'.format(args.nametag)))