Alter dis script

benchmarks/dis/analysis-only.sh

+#!/bin/bash
+
+## =============================================================================
+## Run the DIS benchmarks in 5 steps:
+## 1. Parse the command line and setup environment
+## 2. Detector simulation through npsim
+## 3. Digitization and reconstruction through Juggler
+## 4. Root-based Physics analyses
+## 5. Finalize
+## =============================================================================
+
+## make sure we launch this script from the project root directory
+PROJECT_ROOT="$( cd "$(dirname "$0")" >/dev/null 2>&1 ; pwd -P )"/../..
+pushd ${PROJECT_ROOT}
+
+echo "Running the DIS benchmarks"
+
+## =============================================================================
+## Step 1: Setup the environment variables
+##
+## First parse the command line flags.
+## This sets the following environment variables:
+## - CONFIG:   The specific generator configuration
+## - EBEAM:    The electron beam energy
+## - PBEAM:    The ion beam energy
+source util/parse_cmd.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
+## - 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.
+source options/env.sh
+
+## We also need the following benchmark-specific variables:
+##
+## - BENCHMARK_TAG: Unique identified for this benchmark process.
+## - BEAM_TAG:      Identifier for the chosen beam configuration
+## - INPUT_PATH:    Path for generator-level input to the benchmarks
+## - TMP_PATH:      Path for temporary data (not exported as artifacts)
+## - RESULTS_PATH:  Path for benchmark output figures and files
+##
+## You can read dvmp/env.sh for more in-depth explanations of the variables.
+source benchmarks/dis/env.sh
+
+## Get a unique file names based on the configuration options
+GEN_FILE=${INPUT_PATH}/gen-${CONFIG}_${JUGGLER_N_EVENTS}.hepmc
+
+SIM_FILE=${TMP_PATH}/sim-${CONFIG}.root
+SIM_LOG=${TMP_PATH}/sim-${CONFIG}.log
+
+
+REC_FILE=${TMP_PATH}/rec-${CONFIG}.root
+REC_LOG=${TMP_PATH}/sim-${CONFIG}.log
+
+PLOT_TAG=${CONFIG}
+
+## =============================================================================
+## Step 2: Run the simulation
+echo "Running Geant4 simulation"
+#npsim --runType batch \
+#      --part.minimalKineticEnergy 1000*GeV  \
+#      -v WARNING \
+#      --numberOfEvents ${JUGGLER_N_EVENTS} \
+#      --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
+#      --inputFiles ${GEN_FILE} \
+#      --outputFile ${SIM_FILE}
+#if [ "$?" -ne "0" ] ; then
+#  echo "ERROR running npsim"
+#  exit 1
+#fi
+
+## =============================================================================
+## Step 3: Run digitization & reconstruction
+#echo "Running the digitization and reconstruction"
+## FIXME Need to figure out how to pass file name to juggler from the commandline
+## the tracker_reconstruction.py options file uses the following environment
+## variables:
+## - JUGGLER_SIM_FILE:    input detector simulation
+## - JUGGLER_REC_FILE:    output reconstructed data
+## - JUGGLER_DETECTOR_PATH: Location of the detector geometry
+## - JUGGLER_N_EVENTS:    number of events to process (part of global environment)
+## - JUGGLER_DETECTOR:    detector package (part of global environment)
+#export JUGGLER_SIM_FILE=${SIM_FILE}
+#export JUGGLER_REC_FILE=${REC_FILE}
+#export JUGGLER_DETECTOR_PATH=${DETECTOR_PATH}
+#xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
+#  gaudirun.py options/tracker_reconstruction.py 
+## on-error, first retry running juggler again as there is still a random
+## crash we need to address FIXME
+#if [ "$?" -ne "0" ] ; then
+#  echo "Juggler crashed, retrying..."
+#  xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
+#    gaudirun.py options/tracker_reconstruction.py \
+#    2>&1 > ${REC_LOG}
+#  if [ "$?" -ne "0" ] ; then
+#    echo "ERROR running juggler, both attempts failed"
+#    exit 1
+#  fi
+#fi
+
+## =============================================================================
+## Step 4: Analysis
+## write a temporary configuration file for the analysis script
+echo "Running analysis"
+CONFIG="${TMP_PATH}/${PLOT_TAG}.json"
+cat << EOF > ${CONFIG}
+{
+  "rec_file": "${REC_FILE}",
+  "detector": "${JUGGLER_DETECTOR}",
+  "output_prefix": "${RESULTS_PATH}/${PLOT_TAG}",
+  "test_tag": "${BEAM_TAG}"
+}
+EOF
+#cat ${CONFIG}
+root -b -q "benchmarks/dis/analysis/dis_electrons.cxx+(\"${CONFIG}\")"
+#root -b -q "benchmarks/dis/analysis/dis_electrons.cxx(\"${CONFIG}\")"
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running rec_dis_electron script"
+  exit 1
+fi
+
+## =============================================================================
+## Step 5: finalize
+#echo "Finalizing DIS benchmark"
+
+## Move over reconsturction artifacts as long as we don't have
+## too many events
+#if [ "${JUGGLER_N_EVENTS}" -lt "500" ] ; then 
+#  cp ${REC_FILE} ${RESULTS_PATH}
+#fi
+
+## Always move over log files to the results path
+#cp ${REC_LOG} ${RESULTS_PATH}
+
+## =============================================================================
+## All done!
+echo "DIS benchmarks complete"

benchmarks/dis/analysis/dis.h

@@ -17,6 +17,9 @@
 // promoted to the top-level util library
 namespace util {
 
+  //The functions below were copied from dvmp.h
+
+
   // ADD EXTRA DIS UTILTIY FUNCTIONS HERE
 
   //=========================================================================================================

benchmarks/dis/analysis/dis_electrons.cxx

@@ -14,6 +14,120 @@
 #include <nlohmann/json.hpp>
 #include <string>
 #include <vector>
+#include <eicd/ReconstructedParticleData.h>
+#include <TH1D.h>
+#include <TFitResult.h>
+#include <TRandom3.h>
+#include <algorithm>
+#include <utility>
+
+
+// Get a vector of 4-momenta from the reconstructed data.
+inline auto momenta_from_reconstruction(const std::vector<eic::ReconstructedParticleData>& parts)
+{
+  std::vector<ROOT::Math::PxPyPzEVector> momenta{parts.size()};
+  // transform our reconstructed particle data into 4-momenta
+  std::transform(parts.begin(), parts.end(), momenta.begin(), [](const auto& part) {
+    return ROOT::Math::PxPyPzEVector{part.p.x, part.p.y, part.p.z, part.energy};
+  });
+  return momenta;
+}
+
+// Convert PxPyPzMVector to PxPyPzEVector
+inline auto convertMtoE(const std::vector<ROOT::Math::PxPyPzMVector>& mom)
+{
+  std::vector<ROOT::Math::PxPyPzEVector> momenta{mom.size()};
+  std::transform(mom.begin(), mom.end(), momenta.begin(), [](const auto& part) {
+    return ROOT::Math::PxPyPzEVector{part.Px(), part.Py(), part.Pz(), part.energy()};
+  });
+  return momenta;
+}
+
+// Momentum sorting bool
+bool sort_mom_bool(ROOT::Math::PxPyPzEVector &mom1, ROOT::Math::PxPyPzEVector &mom2)
+{
+  return  mom1.energy() > mom2.energy(); 
+}
+
+// Momentunm sorting function
+inline auto sort_momenta(const std::vector<ROOT::Math::PxPyPzEVector>& mom)
+{
+  std::vector <ROOT::Math::PxPyPzEVector> sort_mom = mom;
+  sort(sort_mom.begin(), sort_mom.end(), sort_mom_bool);
+  return sort_mom;
+}
+
+// Q2 calculation from 4 Vector
+inline auto Q2(const std::vector<ROOT::Math::PxPyPzEVector>& mom)
+{
+  std::vector<double> Q2Vec(mom.size() );
+  ROOT::Math::PxPyPzEVector beamMom = {0, 0, 18, 18};
+  std::transform(mom.begin(), mom.end(), Q2Vec.begin(), [beamMom](const auto& part) {
+    return -(part - beamMom).M2();
+  });
+  return Q2Vec;
+}
+
+// Difference between two 4 Vectors
+inline auto sub(const std::vector<ROOT::Math::PxPyPzEVector>& mom1, const std::vector<ROOT::Math::PxPyPzEVector>& mom2)
+{
+  std::vector<ROOT::Math::PxPyPzEVector> sub(mom1.size() );
+  for (int i = 0; i < sub.size(); i++)
+  {
+    sub[i] = mom1[i] - mom2[i];
+  }
+  return sub;
+}
+
+// Multiplies a double by a gaussian distributed number 
+inline auto randomize(const std::vector<double>& doubleVec)
+{
+  std::vector<double> outVec(doubleVec.size() );
+  TRandom3 rand;
+  std::transform(doubleVec.begin(), doubleVec.end(), outVec.begin(), [&rand](const auto& indouble) {
+    return indouble * rand.Gaus(1.0, 0.2);
+  });
+  return outVec;
+}
+
+//Simulation functions
+/*
+bool mom_sort_sim(dd4pod::Geant4ParticleData& part1, dd4pod::Geant4ParticleData& part2)
+{
+  return (part1.psx*part1.psx + part1.psy*part1.psy + part1.psz*part1.psz > 
+          part2.psx*part2.psx + part2.psy*part2.psy + part2.psz*part2.psz);
+}
+
+inline auto momenta_from_sim(const std::vector<dd4pod::Geant4ParticleData>& parts)
+{
+  std::vector<dd4pod::Geant4ParticleData> sort_parts = parts;
+  sort(sort_parts.begin(), sort_parts.end(), mom_sort_sim);
+  return sort_parts;
+}
+
+inline auto Q2_from_sim(const std::vector<dd4pod::Geant4ParticleData>& parts)
+{
+  std::vector<double> Q2Vec(parts.size() );
+  double beamEnergy = 18;
+  std::transform(parts.begin(), parts.end(), Q2Vec.begin(), [beamEnergy](const auto& part) {
+    double energy = sqrt(part.psx*part.psx + part.psy*part.psy + part.psz*part.psz + part.mass*part.mass);
+    double q2 = pow(beamEnergy - energy, 2.0) - part.psx*part.psx - part.psy*part.psy - pow(part.psz - beamEnergy, 2.0);
+    return -q2;
+  });
+  return Q2Vec;
+}
+
+inline auto elec_PID_sim(const std::vector<dd4pod::Geant4ParticleData>& parts)
+{
+  std::vector<dd4pod::Geant4ParticleData> electrons;
+  for (auto part : parts)
+  {
+    if (part.pdgID == 11){electrons.push_back(part);}
+  }
+  return electrons;
+}
+*/
+
 
 int dis_electrons(const std::string& config_name)
 {
@@ -45,7 +159,22 @@ int dis_electrons(const std::string& config_name)
        {"target", "0.1"}}};
 
   // Run this in multi-threaded mode if desired
-  ROOT::EnableImplicitMT(kNumThreads);
+  //ROOT::EnableImplicitMT(kNumThreads);
+
+  //Particle number enumeration
+  enum sidis_particle_ID
+  {
+    electron = 11,
+    pi_0     = 111,
+    pi_plus  = 211,
+    pi_minus = -211,
+    k_0      = 311,
+    k_plus   = 321,
+    k_minus  = -321,
+    proton   = 2212,
+    neutron  = 2112
+  };
+
 
   const double electron_mass = util::get_pdg_mass("electron");
 
@@ -64,15 +193,63 @@ int dis_electrons(const std::string& config_name)
         return util::momenta_from_tracking(tracks, electron_mass);
       };
 
+/*
+  //Old dataframe
   auto d0 = d.Define("p_rec", momenta_from_tracking, {"outputTrackParameters"})
                 .Define("N", "p_rec.size()")
                 .Define("p_sim", util::momenta_from_simulation, {"mcparticles2"})
                 .Define("mom_sim", util::mom, {"p_sim"})
                 .Define("mom_rec", util::mom, {"p_rec"});
+*/
 
-  auto h_mom_sim = d0.Histo1D({"h_mom_sim", "; GeV; counts", 100, 0, 50}, "mom_sim");
-  auto h_mom_rec = d0.Histo1D({"h_mom_rec", "; GeV; counts", 100, 0, 50}, "mom_rec");
+  auto d0 = d.Define("p_recon", momenta_from_reconstruction, {"DummyReconstructedParticles"})
+                .Define("p_recon_sort", sort_momenta, {"p_recon"})
+                .Define("Q2_recon", Q2, {"p_recon_sort"})
+                .Define("Q2_recon_rand", randomize, {"Q2_recon"})
+                .Define("elec_Q2_recon_rand", "Q2_recon_rand[0]")
+                .Define("p_sim_M", util::momenta_from_simulation, {"mcparticles2"})
+                .Define("p_sim", convertMtoE, {"p_sim_M"})
+                .Define("p_sim_sort", sort_momenta, {"p_sim"})
+                .Define("Q2_sim", Q2, {"p_sim_sort"})
+                .Define("elec_Q2_sim", "Q2_sim[0]")
+                .Define("Q2_diff", "(elec_Q2_recon_rand - elec_Q2_sim)/elec_Q2_sim")
+                .Define("p_diff", sub, {"p_recon","p_sim"});
+                /*
+                .Define("electrons_sim", elec_PID_sim, {"sorted_sim"})
+                .Define("Q2_sim_elec_pid", Q2_from_sim, {"electrons_sim"})
+                .Define("elec_Q2_sim_pid", "Q2_sim_elec_pid[0]");
+                */
+  //Testing script
+/*
+  auto dis = d0.Display({"p_diff"});
+  //dis -> Print();
+  cout << *d0.Max<double>("elec_Q2_recon_rand") << " " << *d0.Min<double>("elec_Q2_recon_rand") << endl;
+  cout << *d0.Max<double>("elec_Q2_sim") << " " << *d0.Min<double>("elec_Q2_sim") << endl;
+/**/
+  //cout << *d0.Count() << endl;
+  //Momentum
+  //auto h_mom_sim = d0.Histo1D({"h_mom_sim", "; GeV; counts", 100, 0, 50}, "mom_sim");
+  //auto h_mom_rec = d0.Histo1D({"h_mom_rec", "; GeV; counts", 100, 0, 50}, "mom_rec");
+  //Q2
+  auto h_Q2_sim = d0.Histo1D({"h_Q2_sim", "; GeV; counts", 100, -5, 25}, "elec_Q2_sim");
+  auto h_Q2_rec = d0.Histo1D({"h_Q2_rec", "; GeV; counts", 100, -5, 25}, "elec_Q2_recon_rand");
+  auto h_Q2_res = d0.Histo1D({"h_Q2_res", ";    ; counts", 100, -10, 10}, "Q2_diff");
+/*
+  TH1D *h_Q2_res = (TH1D *)h_Q2_sim -> Clone();
+  TH1D *h_Q2_rec_copy = (TH1D *)h_Q2_rec -> Clone();
+  h_Q2_res -> Scale(1.0 / h_Q2_res -> Integral() );
+  h_Q2_res -> Add(h_Q2_rec_copy, -1.0 / h_Q2_rec_copy -> Integral() );
+ */
+  TFitResultPtr f1  = h_Q2_res -> Fit("gaus", "S");
+  f1 -> Print("V");
+/*
+  printf("chisq %f A %f mean %f sigma %f \n", f1 -> Chi2(), 
+                                              f1 -> GetParameter(0), 
+                                              f1 -> GetParameter(1),
+                                              f1 -> GetParameter(2));
+*/
 
+  
   auto c = new TCanvas();
 
   // Plot our histograms.
@@ -83,8 +260,10 @@ int dis_electrons(const std::string& config_name)
     c.cd();
     // gPad->SetLogx(false);
     gPad->SetLogy(true);
-    auto& h1 = *h_mom_sim;
-    auto& h2 = *h_mom_rec;
+    //auto& h1 = *h_mom_sim;
+    //auto& h2 = *h_mom_rec;
+    auto& h1 = *h_Q2_sim;
+    auto& h2 = *h_Q2_rec;
     // histogram style
     h1.SetLineColor(plot::kMpBlue);
     h1.SetLineWidth(2);

benchmarks/dis/dis.sh

@@ -111,6 +111,8 @@ cat << EOF > ${CONFIG}
   "rec_file": "${REC_FILE}",
   "detector": "${JUGGLER_DETECTOR}",
   "output_prefix": "${RESULTS_PATH}/${PLOT_TAG}",
+  "ebeam": ${EBEAM},
+  "pbeam": ${PBEAM},
   "test_tag": "${BEAM_TAG}"
 }
 EOF