diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 9a3563cbc123d34a2cfcaaec4b5d37aaac43959f..f68226e47101879610edfb511f8fb049030297be 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -93,6 +93,7 @@ include: - local: 'benchmarks/crystal_calorimeter/config.yml' - local: 'benchmarks/pid/config.yml' - local: 'benchmarks/timing/config.yml' + - local: 'benchmarks/b0_tracker/config.yml' deploy_results: stage: deploy diff --git a/benchmarks/b0_tracker/config.yml b/benchmarks/b0_tracker/config.yml new file mode 100644 index 0000000000000000000000000000000000000000..5445e7bfe06459c66789ff0351a2f13e0a2290e7 --- /dev/null +++ b/benchmarks/b0_tracker/config.yml @@ -0,0 +1,25 @@ +sim:b0_tracker: + stage: simulate + extends: .det_benchmark + script: + - bash benchmarks/b0_tracker/forward_protons.sh + +bench:b0_tracker: + stage: benchmarks + extends: .det_benchmark + needs: + - ["sim:roman_pot"] + script: + - echo "tracking analysis script here later" + #- root -b -q benchmarks/trackers/simple_tracking.cxx+ + +results:b0_tracker: + extends: .det_benchmark + stage: collect + needs: + - ["bench:b0_tracker"] + script: + - echo "Collecting results for B0 Tracker" + + + diff --git a/benchmarks/b0_tracker/forward_protons.sh b/benchmarks/b0_tracker/forward_protons.sh new file mode 100755 index 0000000000000000000000000000000000000000..fd7db809993d003b60ecf1b72c4849a3df964970 --- /dev/null +++ b/benchmarks/b0_tracker/forward_protons.sh @@ -0,0 +1,51 @@ +#!/bin/bash + +source .local/bin/env.sh + +if [[ ! -n "${JUGGLER_DETECTOR}" ]] ; then + export JUGGLER_DETECTOR="topside" +fi + +if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then + export JUGGLER_N_EVENTS=100 +fi + +export FILE_NAME_TAG="forward_protons" +export JUGGLER_GEN_FILE="${FILE_NAME_TAG}.hepmc" +export JUGGLER_SIM_FILE="sim_${FILE_NAME_TAG}.root" +export JUGGLER_REC_FILE="rec_${FILE_NAME_TAG}.root" + +echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}" +echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}" + +# Generate the input events +root -b -q "benchmarks/b0_tracker/scripts/gen_${FILE_NAME_TAG}.cxx(${JUGGLER_N_EVENTS}, \"${FILE_NAME_TAG}.hepmc\")" +if [[ "$?" -ne "0" ]] ; then + echo "ERROR running script: generating input events" + exit 1 +fi + +# Run geant4 simulations +npsim --runType batch \ + -v WARNING \ + --part.minimalKineticEnergy 0.5*GeV \ + --numberOfEvents ${JUGGLER_N_EVENTS} \ + --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \ + --inputFiles ${FILE_NAME_TAG}.hepmc \ + --outputFile sim_output/${JUGGLER_SIM_FILE} + +if [[ "$?" -ne "0" ]] ; then + echo "ERROR running npdet" + exit 1 +fi + +# Directory for plots +mkdir -p results + +# Plot the input events +root -b -q "benchmarks/b0_tracker/scripts/b0_tracker_hits.cxx(\"sim_output/sim_${FILE_NAME_TAG}.root\")" +if [[ "$?" -ne "0" ]] ; then + echo "ERROR running script: events" + exit 1 +fi + diff --git a/benchmarks/b0_tracker/scripts/b0_tracker_hits.cxx b/benchmarks/b0_tracker/scripts/b0_tracker_hits.cxx new file mode 100644 index 0000000000000000000000000000000000000000..b4984adaa3ef6c31788a38c6ba35f3d0f3e9e466 --- /dev/null +++ b/benchmarks/b0_tracker/scripts/b0_tracker_hits.cxx @@ -0,0 +1,64 @@ +R__LOAD_LIBRARY(libfmt.so) +#include "fmt/core.h" + +#include "ROOT/RDataFrame.hxx" +#include "Math/Vector3D.h" +//#include "Math/Vector4D.h" +//#include "Math/VectorUtil.h" +#include "TCanvas.h" +//#include "TLegend.h" +//#include "TMath.h" +//#include "TRandom3.h" +//#include "TFile.h" +//#include "TH1F.h" +//#include "TH1D.h" +//#include "TTree.h" +#include "TChain.h" +//#include "TF1.h" + +#include "dd4pod/TrackerHitCollection.h" + +#include "common_bench/particles.h" +#include "common_bench/benchmark.h" +#include "common_bench/mt.h" +#include "common_bench/util.h" + + +#include "dd4pod/TrackerHitCollection.h" + +void b0_tracker_hits(const char* fname = "./sim_output/sim_forward_protons.root"){ + + ROOT::EnableImplicitMT(); // Tell ROOT you want to go parallel + double degree = TMath::Pi()/180.0; + + TChain* t = new TChain("events"); + t->Add(fname); + + ROOT::RDataFrame d0(*t); + + auto hits_eta = [&](const std::vector<dd4pod::TrackerHitData>& hits) { + std::vector<double> result; + for (const auto& h : hits) { + ROOT::Math::XYZVector vec(h.position.x,h.position.y,h.position.z); + result.push_back(vec.eta()); + std::cout << vec.eta() << "\n"; + } + return result; + }; + + auto d1 = d0.Define("hits_eta", hits_eta, {"B0TrackerHits"}); + + auto h1 = d1.Histo1D({"h1", "hits_eta", 100, 0,20}, "hits_eta"); + TCanvas* c = new TCanvas(); + h1->DrawCopy(); + c->SaveAs("results/b0_tracker_hits_eta.png"); + c->SaveAs("results/b0_tracker_hits_eta.pdf"); + auto n1 = h1->GetMean(); + std::cout << "Pseudorapidity of hits: " << n1 << std::endl; + + //if (n1 < 5) { + // std::quick_exit(1); + //} + +} + diff --git a/benchmarks/b0_tracker/scripts/gen_forward_protons.cxx b/benchmarks/b0_tracker/scripts/gen_forward_protons.cxx new file mode 100644 index 0000000000000000000000000000000000000000..a8884ca1646b18e1a8851d1e616f92c03df4afd8 --- /dev/null +++ b/benchmarks/b0_tracker/scripts/gen_forward_protons.cxx @@ -0,0 +1,85 @@ +#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" + +#include "common_bench/particles.h" + +using namespace HepMC3; + +/** Generate electrons in the central region. + * This is for testing detectors in the "barrel" region. + */ +void gen_forward_protons(int n_events = 100, + const char* out_fname = "forward_protons.hepmc") +{ + double cos_theta_min = std::cos(1.0*(M_PI/180.0)); + double cos_theta_max = std::cos(0.0*(M_PI/180.0)); + + const double M_p = common_bench::particleMap.at(2212).mass; + + 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 + GenParticlePtr p1 = + std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 2212, 4); + GenParticlePtr p2 = std::make_shared<GenParticle>( + FourVector(0.0, 0.0, 0.0, M_p), 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 + (M_p * M_p))), 2212, 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; +}