diff --git a/benchmarks/tof/scripts/gen_tof_forward_hits.cxx b/benchmarks/tof/scripts/gen_tof_forward_hits.cxx new file mode 100644 index 0000000000000000000000000000000000000000..c679611d6c99c8dde8060356329204152dd253c3 --- /dev/null +++ b/benchmarks/tof/scripts/gen_tof_forward_hits.cxx @@ -0,0 +1,100 @@ +#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_tof_forward_hits(int n_events = 100, + const char* out_fname = "tof_forward_hits.hepmc", + int n_parts = 4) +{ + double cos_theta_min = std::cos( 1.0*(M_PI/180.0)); + double cos_theta_max = std::cos(189.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 + Double_t p = r1->Uniform(0, 8); + Double_t costh = -1*cos(12.033/180*3.1415926); //r1->Uniform(cos_theta_min, cos_theta_max); + Double_t th = std::acos(costh); + Double_t phi = r1->Uniform(0.0, 2.0 * M_PI); + 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); + + 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 + phi = r1->Uniform(0.0, 2.0 * M_PI); + px = p * std::cos(phi) * std::sin(th); + py = p * std::sin(phi) * std::sin(th); + pz = p * std::cos(th); + GenParticlePtr p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.000511 * 0.000511))), + ((ip % 2 == 0) ? 11 : -11), 1); + + phi = r1->Uniform(0.0, 2.0 * M_PI); + px = p * std::cos(phi) * std::sin(th); + py = p * std::sin(phi) * std::sin(th); + pz = p * std::cos(th); + GenParticlePtr p4 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.938272 * 0.938272))), + ((ip % 2 == 0) ? 2212 : -2212), 1); + + phi = r1->Uniform(0.0, 2.0 * M_PI); + px = p * std::cos(phi) * std::sin(th); + py = p * std::sin(phi) * std::sin(th); + pz = p * std::cos(th); + GenParticlePtr p5 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.493677 * 0.493677))), + ((ip % 2 == 0) ? 321 : -321), 1); + phi = r1->Uniform(0.0, 2.0 * M_PI); + px = p * std::cos(phi) * std::sin(th); + py = p * std::sin(phi) * std::sin(th); + pz = p * std::cos(th); + GenParticlePtr p6 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.139570 * 0.139570))), + ((ip % 2 == 0) ? 211 : -211), 1); + + GenVertexPtr v1 = std::make_shared<GenVertex>(); + v1->add_particle_in(p1); + v1->add_particle_in(p2); + v1->add_particle_out(p3); + v1->add_particle_out(p4); + v1->add_particle_out(p5); + v1->add_particle_out(p6); + 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/tof/scripts/sim_tof_forward_hits.cxx b/benchmarks/tof/scripts/sim_tof_forward_hits.cxx new file mode 100644 index 0000000000000000000000000000000000000000..dc8f508cebc9094ffb83713cba87fadffaf84a14 --- /dev/null +++ b/benchmarks/tof/scripts/sim_tof_forward_hits.cxx @@ -0,0 +1,203 @@ +#include "ROOT/RDataFrame.hxx" +#include "TCanvas.h" +#include "TH1D.h" +#include "TLegend.h" +#include "TProfile.h" + +#include <cstdlib> +#include <iostream> + +R__LOAD_LIBRARY(libeicd.so) +R__LOAD_LIBRARY(libDD4pod.so) + +#include <fmt/format.h> + +#include "dd4pod/Geant4ParticleCollection.h" +#include "dd4pod/TrackerHitCollection.h" +#include "dd4pod/TrackerHitData.h" +#include "eicd/ClusterCollection.h" +#include "eicd/ClusterData.h" +#include "eicd/TrackParametersCollection.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; +}; + +// Create dataframe with extra definitions for hit multiplicities +ROOT::RDF::RNode add_subsystems(ROOT::RDF::RNode df, std::vector<std::pair<std::string, std::string>> hitcols) { + if (hitcols.empty()) { + return df; + } + const auto [name, collection] = hitcols.back(); + std::cout << " --> registering subsystem " << collection << "\n"; + auto df2 = df.Define("N_" + name, [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size(); }, {collection}); + hitcols.pop_back(); + return add_subsystems(df2, hitcols); +}; + +int sim_tof_forward_hits(const char* fname = "sim_tof_forward_hits.root") { + + ROOT::EnableImplicitMT(); + ROOT::RDataFrame df("events", fname); + + // detect detector setup + std::string detector = "default"; + if (const char* env_detector = std::getenv("JUGGLER_DETECTOR_VERSION")) { + if (detector.find("acadia") != std::string::npos) { + detector = "acadia"; + } else if (detector.find("canyonlands") != std::string::npos) { + detector = "canyonlands"; + } + } + std::cout << "sim_tof_forward_hits: detector set to " << detector << std::endl; + + // minimal hit collection setup + std::vector<std::pair<std::string, std::string>> hitCollections{{"tof_forward", "TOFForwardHits"}}; + //{"vtx_barrel", "VertexBarrelHits"}, + //{"trk_barrel", "TrackerBarrelHits"}, + //{"tof_barrel", "TOFBarrelHits"}, + //{"trk_endcap", "TrackerEndcapHits"}, + //{"gem_endcap", "GEMTrackerEndcapHits"}}; + + // append extra hit collections based on detector setup + if (detector == "acadia") { + hitCollections.push_back({"vtx_endcap", "VertexEndcapHits"}); + } else if (detector == "canyonlands" || detector == "default") { + hitCollections.push_back({"mm_barrel", "MPGDTrackerBarrelHits"}); + } + + 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]"); + auto df1 = add_subsystems(df0, hitCollections); + + // get 1D and 2D histogram definitions + std::vector<decltype(df1.Histo2D({}, "", ""))> h2D; + std::vector<decltype(df1.Histo1D(""))> hnth; + std::vector<decltype(df1.Histo1D(""))> hn; + + for (const auto& [name, col] : hitCollections) { + hnth.push_back( + df1.Histo1D({(name + "_nhits_vs_theta").c_str(), "; #theta [deg.]; ", 20, 0, 180}, "theta0", "N_" + name)); + hn.push_back(df1.Histo1D({(name + "_nhits").c_str(), "; Nhits; #", 20, 0, 20}, "N_" + name)); + h2D.push_back(df1.Histo2D({(name + "_x_vs_y").c_str(), "; x ; y ", 100, -600, 600, 100, -600, 600}, + col + ".position.x", col + ".position.y")); + h2D.push_back(df1.Histo2D({(name + "_x_vs_z").c_str(), "; z ; x ", 100, -2000, 2000, 100, -900, 900}, + col + ".position.z", col + ".position.x")); + h2D.push_back(df1.Histo2D({(name + "_time_vs_length").c_str(), "; time ; length ", 100, 0, 10, 100, 0, 1}, + col + ".truth.time", col + ".length")); + + } + auto hth = df1.Histo1D({"theta0", "; #theta [deg.]", 20, 0, 180}, "theta0"); + + // Draw multiplicity versus theta + { + TCanvas c; + THStack hs{"n_hits", "; #theta [deg.] "}; + int idx = 1; + for (auto& h : hnth) { + h->Divide(&*hth); + h->SetLineColor(idx++); + hs.Add((TH1D*)h->Clone()); + } + hs.Draw("nostack, hist"); + c.BuildLegend(); + c.SaveAs("results/tof/sim_tof_forward_hits_n_hits_vs_theta.png"); + c.SaveAs("results/tof/sim_tof_forward_hits_n_hits_vs_theta.pdf"); + } + // Draw nhits + { + TCanvas c; + THStack hs{"n_hits", "; NHits "}; + int idx = 1; + for (auto& h : hn) { + h->SetLineColor(idx++); + hs.Add((TH1D*)h->Clone()); + } + hs.Draw("nostack, hist"); + c.BuildLegend(); + c.SaveAs("results/tof/sim_tof_forward_hits_n_hits.png"); + c.SaveAs("results/tof/sim_tof_forward_hits_n_hits.pdf"); + } + + // Draw total multiplicity versus theta + { + TCanvas c; + hth->DrawClone(); + c.BuildLegend(); + c.SaveAs("results/tof/sim_tof_forward_hits_theta.png"); + c.SaveAs("results/tof/sim_tof_forward_hits_theta.pdf"); + } + + for (auto& h : h2D) { + TCanvas c; + h->DrawClone("colz"); + c.SaveAs(fmt::format("results/tof/sim_tof_forward_hits_{}.png", h->GetName()).c_str()); + c.SaveAs(fmt::format("results/tof/sim_tof_forward_hits_{}.pdf", h->GetName()).c_str()); + } + + return 0; +} diff --git a/benchmarks/tof/tof_forward_hits.sh b/benchmarks/tof/tof_forward_hits.sh new file mode 100644 index 0000000000000000000000000000000000000000..04af3af842a4fc25feed29d8e63e01ca1b95937b --- /dev/null +++ b/benchmarks/tof/tof_forward_hits.sh @@ -0,0 +1,121 @@ +#!/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_N_EVENTS=100 +export JUGGLER_FILE_NAME_TAG="tof_forward_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/tof/scripts/gen_tof_forward_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/tof + rootls -t ${JUGGLER_SIM_FILE} + root -b -q "benchmarks/tof/scripts/sim_tof_forward_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}") +echo "filesize" $filesize +if [[ "${JUGGLER_N_EVENTS}" -lt "31000" ]] ; then + # file must be less than 10 MB to upload + if [[ "${root_filesize}" -lt "10000000000" ]] ; then + cp ${JUGGLER_SIM_FILE} results/. + fi +fi +ls -ltr ${JUGGLER_SIM_FILE} + +