diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index af39d9347d77c759a43806401c9a838c4be1eb23..74ed1de8542a4c562e891bdd27e7ca297a6c3811 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -98,6 +98,7 @@ get_data: - ls -lrtha include: + - local: 'benchmarks/tracking_detectors/config.yml' - local: 'benchmarks/barrel_ecal/config.yml' - local: 'benchmarks/barrel_hcal/config.yml' - local: 'benchmarks/roman_pots/config.yml' diff --git a/benchmarks/tracking_detectors/config.yml b/benchmarks/tracking_detectors/config.yml new file mode 100644 index 0000000000000000000000000000000000000000..031cc85039579cca9b0e4e2429adaf778ffc34bd --- /dev/null +++ b/benchmarks/tracking_detectors/config.yml @@ -0,0 +1,13 @@ +tracking_detectors:sim_track_hits: + extends: .det_benchmark + stage: simulate + script: + - bash benchmarks/tracking_detectors/track_hits.sh --sim-only + +tracking_detectors:plot_hits: + extends: .det_benchmark + stage: benchmarks + needs: + - ["tracking_detectors:sim_track_hits"] + script: + - bash benchmarks/tracking_detectors/track_hits.sh --analysis diff --git a/benchmarks/tracking_detectors/scripts/gen_track_hits.cxx b/benchmarks/tracking_detectors/scripts/gen_track_hits.cxx new file mode 100644 index 0000000000000000000000000000000000000000..13bc31e70598c6bb505d8618603b269f930e447e --- /dev/null +++ b/benchmarks/tracking_detectors/scripts/gen_track_hits.cxx @@ -0,0 +1,81 @@ +#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 multiple electrons/positron tracks in the central region. + * This is for testing detectors in the "barrel" region. + */ +void gen_track_hits(int n_events = 100, + const char* out_fname = "track_hits.hepmc", + int n_parts = 2) +{ + double cos_theta_min = std::cos( 10.0*(M_PI/180.0)); + double cos_theta_max = std::cos(170.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 - electron + // pdgid 111 - pi0 + // pdgid 2212 - proton + for (int ip = 0; ip < n_parts; ip++) { + 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 11 - electron 0.510 MeV/c^2 + GenParticlePtr p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.000511 * 0.000511))), + ((ip % 2 == 0) ? 11 : -11), 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_detectors/scripts/sim_track_hits.cxx b/benchmarks/tracking_detectors/scripts/sim_track_hits.cxx new file mode 100644 index 0000000000000000000000000000000000000000..440c40720cc04774ef017470b99fe3fee3ba720b --- /dev/null +++ b/benchmarks/tracking_detectors/scripts/sim_track_hits.cxx @@ -0,0 +1,205 @@ +#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].ps.x * in[i].ps.x + in[i].ps.y * in[i].ps.y)); + } + return result; +} + +auto momentum = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) { + std::vector<double> result; + for (size_t i = 0; i < in.size(); ++i) { + result.push_back(in[i].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].ps.x, in[i].ps.y, in[i].ps.z, 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 sim_track_hits(const char* fname = "sim_track_hits.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("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 hEndcap_x_vs_y = df0.Histo2D({"hEndcap_x_vs_y", "; x ; y ", 100, -900, 900,100, -900, 900 }, "TrackerEndcapHits.position.x", "TrackerEndcapHits.position.y"); + //auto hvtxBarrel_x_vs_y = df0.Histo2D({"hvtxBarrel_x_vs_y", "; x ; y ", 100, -900, 900,100, -900, 900 }, "VertexBarrelHits.position.x", "VertexBarrelHits.position.y"); + //auto hvtxEndcap_x_vs_y = df0.Histo2D({"hvtxEndcap_x_vs_y", "; x ; y ", 100, -900, 900,100, -900, 900 }, "VertexEndcaplHits.position.x", "VertexEndcapHits.position.y"); + //auto hAllHits_x_vs_y = df0.Histo2D({"hAllHitsx_vs_y", "; x ; y ", 100, -900, 900,100, -900, 900 }, "allHits.position.x", "allHits.position.y"); + + auto hBarrel_x_vs_z = df0.Histo2D({"hBarrel_x_vs_z", "; z ; x ", 100, -900, 900,100, -900, 900 }, "TrackerBarrelHits.position.z" , "TrackerBarrelHits.position.y"); + auto hEndcap_x_vs_z = df0.Histo2D({"hEndcap_x_vs_z", "; z ; x ", 100, -900, 900,100, -900, 900 }, "TrackerEndcapHits.position.z", "TrackerEndcapHits.position.y"); + //auto hvtxBarrel_x_vs_z = df0.Histo2D({"hvtxBarrel_x_vs_z", "; z ; x ", 100, -900, 900,100, -900, 900 }, "VertexBarrelHits.position.z" , "VertexBarrelHits.position.y" ); + //auto hvtxEndcap_x_vs_z = df0.Histo2D({"hvtxEndcap_x_vs_z", "; z ; x ", 100, -900, 900,100, -900, 900 }, "VertexEndcaplHits.position.z" , "VertexEndcapHits.position.y" ); + //auto hAllHits_x_vs_z = df0.Histo2D({"hAllHitsx_vs_z","; z ; x ", 100, -900, 900,100, -900, 900 }, "allHits.position.z" , "allHits.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_detectors/sim_track_hits_n_hits_vs_theta.png"); + c->SaveAs("results/tracking_detectors/sim_track_hits_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_detectors/sim_track_hits_theta.png"); + c->SaveAs("results/tracking_detectors/sim_track_hits_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_detectors/sim_track_hits_nhits.png"); + c->SaveAs("results/tracking_detectors/sim_track_hits_nhits.pdf"); + + c = new TCanvas(); + hBarrel_x_vs_y->DrawCopy("colz"); + c->SaveAs("results/tracking_detectors/sim_track_hits_trkBarrel_xy.png"); + c->SaveAs("results/tracking_detectors/sim_track_hits_trkBarrel_xy.pdf"); + + c = new TCanvas(); + hBarrel_x_vs_y->DrawCopy("colz"); + hEndcap_x_vs_y->DrawCopy("colz same"); + //hvtxBarrel_x_vs_y->DrawCopy("colz same"); + //hvtxEndcap_x_vs_y->DrawCopy("colz same"); + + c->SaveAs("results/tracking_detectors/sim_track_hits_Hits_xy.png"); + c->SaveAs("results/tracking_detectors/sim_track_hits_Hits_xy.pdf"); + + //hAllHits_x_vs_z->DrawCopy("colz"); + hBarrel_x_vs_z->DrawCopy("colz"); + hEndcap_x_vs_z->DrawCopy("colz same"); + c->SaveAs("results/tracking_detectors/sim_track_hits_Hits_xz.png"); + c->SaveAs("results/tracking_detectors/sim_track_hits_Hits_xz.pdf"); + return 0; +} diff --git a/benchmarks/tracking_detectors/track_hits.sh b/benchmarks/tracking_detectors/track_hits.sh new file mode 100644 index 0000000000000000000000000000000000000000..a7a54e2e31333c9c6e0d868ae1de3a637ffa7f64 --- /dev/null +++ b/benchmarks/tracking_detectors/track_hits.sh @@ -0,0 +1,118 @@ +#!/bin/bash + +#set -o errexit +#set -o pipefail + +function print_the_help { + echo "USAGE: ${0} [--sim-only] " + echo "OPTIONS: " + echo " --sim-only Only run up to the simulation " + echo " --analysis Only run the analysis scripts " + exit +} + +ANALYSIS_ONLY= +SIM_ONLY= +POSITIONAL=() + +while [[ $# -gt 0 ]] +do + key="$1" + + case $key in + -h|--help) + shift # past argument + print_the_help + ;; + --sim-only) + shift # past argument + SIM_ONLY=1 + ;; + --analysis) + shift # past argument + ANALYSIS_ONLY=1 + ;; + *) # 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. + +if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then + export JUGGLER_N_EVENTS=100 +fi + +export JUGGLER_FILE_NAME_TAG="track_hits" +export JUGGLER_GEN_FILE="${LOCAL_DATA_PATH}/${JUGGLER_FILE_NAME_TAG}.hepmc" + +export JUGGLER_SIM_FILE="${LOCAL_DATA_PATH}/sim_${JUGGLER_FILE_NAME_TAG}.root" + +echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}" +echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}" + + +if [ -z "${ANALYSIS_ONLY}" ] ; then + + echo "Generating Events" + ## generate the input events + root -b -q "benchmarks/tracking_detectors/scripts/gen_track_hits.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 + + +if [ -z "${SIM_ONLY}" ] ; then + + echo "Running analysis scripts" + + mkdir -p results/tracking_detectors + rootls -t ${JUGGLER_SIM_FILE} + root -b -q "benchmarks/tracking_detectors/scripts/sim_track_hits.cxx(\"${JUGGLER_SIM_FILE}\")" + if [[ "$?" -ne "0" ]] ; then + echo "ERROR running root script" + exit 1 + fi + +fi + +root_filesize=$(stat --format=%s "${JUGGLER_SIM_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_SIM_FILE} results/. + fi +fi + +