Adding dis benchmark.

- Currently throwing (central) electrons.
- This will be replace with pythia8 input file.
  This should be added to the TODO list

dis/config.yml

-dis:dummy_test:
+dis:run_test:
   stage: analyze
   timeout: 1 hours
   script:
-    - dis/dummy_test.sh
-
-dis:dummy_test2:
-  stage: analyze
-  timeout: 1 hours
-  script:
-    - dis/dummy_test2.sh
-
-dis:dummy_fail_test:
-  stage: analyze
-  timeout: 1 hours
-  allow_failure: true
-  script:
-    - dis/dummy_fail_test.sh
+    - bash dis/dis.sh
 
 dis:results:
   stage: collect
-  needs: ["dis:dummy_test", "dis:dummy_test2", "dis:dummy_fail_test"]
+  needs: ["dis:run_test"]
   script:
     - echo "All DIS benchmarks successful"

dis/dis.sh

+#!/bin/bash
+
+if [[ ! -n  "${JUGGLER_DETECTOR}" ]] ; then 
+  export JUGGLER_DETECTOR="topside"
+fi
+
+if [[ ! -n  "${JUGGLER_N_EVENTS}" ]] ; then 
+  export JUGGLER_N_EVENTS=100
+fi
+
+if [[ ! -n  "${JUGGLER_INSTALL_PREFIX}" ]] ; then 
+  export JUGGLER_INSTALL_PREFIX="/usr/local"
+fi
+
+export JUGGLER_FILE_NAME_TAG="dis"
+export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
+
+export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root"
+export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
+
+echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
+echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
+
+
+### Build the detector constructors.
+git clone https://eicweb.phy.anl.gov/EIC/detectors/${JUGGLER_DETECTOR}.git
+git clone https://eicweb.phy.anl.gov/EIC/detectors/accelerator.git
+pushd ${JUGGLER_DETECTOR}
+ln -s ../accelerator/eic
+popd
+mkdir ${JUGGLER_DETECTOR}/build
+pushd ${JUGGLER_DETECTOR}/build
+cmake ../. -DCMAKE_INSTALL_PREFIX=/usr/local && make -j30 install
+popd
+
+# generate the input events
+# temporary standin until hepmc output from pythia is generated.
+root -b -q "dis/scripts/gen_central_electrons.cxx(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running script"
+  exit 1
+fi
+
+#
+pushd ${JUGGLER_DETECTOR}
+
+## run geant4 simulations
+npsim --runType batch \
+      --part.minimalKineticEnergy 1000*GeV  \
+      -v WARNING \
+      --numberOfEvents ${JUGGLER_N_EVENTS} \
+      --compactFile ${JUGGLER_DETECTOR}.xml \
+      --inputFiles ../${JUGGLER_FILE_NAME_TAG}.hepmc \
+      --outputFile  ${JUGGLER_SIM_FILE}
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running script"
+  exit 1
+fi
+
+# Need to figure out how to pass file name to juggler from the commandline
+xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
+  gaudirun.py ../options/tracker_reconstruction.py
+
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running juggler"
+  exit 1
+fi
+ls -l
+popd
+
+pwd
+mkdir -p results/dis
+
+root -b -q "dis/scripts/rec_dis_electrons.cxx(\"${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}\")"
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running root script"
+  exit 1
+fi
+
+if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then 
+cp ${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE} results/.
+fi
+

dis/dummy_fail_test.sh

-#!/bin/bash
-
-echo "Dummy Test..."
-echo "..."
-echo "Fails!"
-
-exit 1

dis/dummy_test.sh

-#!/bin/bash
-
-echo "Dummy Test..."
-echo "..."
-echo "Passes!"
-
-#exit 1

dis/dummy_test2.sh

-#!/bin/bash
-
-echo "Dummy Test number 2..."
-echo "..."
-echo "Passes!"
-
-#exit 1

dis/scripts/gen_central_electrons.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>
+
+using namespace HepMC3;
+
+/** Generate electrons in the central region.
+ *  This is for testing detectors in the "barrel" region.
+ */
+void gen_central_electrons(int n_events = 100, 
+                     const char* out_fname = "central_electrons.hepmc")
+{
+  double cos_theta_min = std::cos( 50.0*(M_PI/180.0));
+  double cos_theta_max = std::cos(130.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
+    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))),
+        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;
+}

dis/scripts/rec_dis_electrons.cxx

+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "eicd/TrackParametersCollection.h"
+#include "eicd/ClusterCollection.h"
+#include "eicd/ClusterData.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].psx * in[i].psx + in[i].psy * in[i].psy));
+  }
+  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].psx, in[i].psy, in[i].psz, 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;
+};
+
+int rec_dis_electrons(const char* fname = "topside/rec_central_electrons.root")
+{
+
+  ROOT::EnableImplicitMT();
+  ROOT::RDataFrame df("events", fname);
+
+  auto df0 = df.Define("isThrown", "mcparticles2.genStatus == 1")
+                 .Define("thrownParticles", "mcparticles2[isThrown]")
+                 .Define("thrownP", fourvec, {"thrownParticles"})
+                 .Define("p_thrown", momentum, {"thrownP"})
+                 .Define("nTracks", "outputTrackParameters.size()")
+                 .Define("p_track", p_track, {"outputTrackParameters"})
+                 .Define("p_track1", p_track, {"outputTrackParameters1"})
+                 .Define("p_track2", p_track, {"outputTrackParameters2"})
+                 .Define("delta_p",delta_p, {"p_track", "p_thrown"})
+                 .Define("delta_p1",delta_p, {"p_track1", "p_thrown"})
+                 .Define("delta_p2",delta_p, {"p_track2", "p_thrown"});
+
+  auto h_nTracks = df0.Histo1D({"h_nTracks", "; N tracks ", 10, 0, 10}, "nTracks");
+  auto h_pTracks = df0.Histo1D({"h_pTracks", "; GeV/c ", 100, 0, 10}, "p_track");
+
+  auto h_delta_p  = df0.Histo1D({"h_delta_p", "; GeV/c ",  100, -10, 10}, "delta_p");
+  auto h_delta_p1 = df0.Histo1D({"h_delta_p1", "; GeV/c ", 100, -10, 10}, "delta_p1");
+  auto h_delta_p2 = df0.Histo1D({"h_delta_p2", "; GeV/c ", 100, -10, 10}, "delta_p2");
+
+  auto c = new TCanvas();
+
+  h_nTracks->DrawCopy();
+  c->SaveAs("results/dis/rec_central_electrons_nTracks.png");
+  c->SaveAs("results/dis/rec_central_electrons_nTracks.pdf");
+
+  h_pTracks->DrawCopy();
+  c->SaveAs("results/dis/rec_central_electrons_pTracks.png");
+  c->SaveAs("results/dis/rec_central_electrons_pTracks.pdf");
+
+  THStack * hs = new THStack("hs_delta_p","; GeV/c "); 
+  TH1D* h1 = (TH1D*) h_delta_p->Clone();
+  hs->Add(h1);
+  h1 = (TH1D*) h_delta_p1->Clone();
+  h1->SetLineColor(2);
+  hs->Add(h1);
+  h1 = (TH1D*) h_delta_p2->Clone();
+  h1->SetLineColor(4);
+  h1->SetFillStyle(3001);
+  h1->SetFillColor(4);
+  hs->Add(h1);
+  hs->Draw("nostack");
+  c->SaveAs("results/dis/rec_central_electrons_delta_p.png");
+  c->SaveAs("results/dis/rec_central_electrons_delta_p.pdf");
+
+  return 0;
+}