B0tracker basics

.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

benchmarks/b0_tracker/config.yml

+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"
+
+
+

benchmarks/b0_tracker/forward_protons.sh

+#!/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
+

benchmarks/b0_tracker/scripts/b0_tracker_hits.cxx

+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);
+  //}
+
+}
+

benchmarks/b0_tracker/scripts/gen_forward_protons.cxx

+#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;
+}