From e66cdaf95470ac3788707c5ef285c4f51e13b1ed Mon Sep 17 00:00:00 2001
From: Whitney Armstrong <warmstrong@anl.gov>
Date: Wed, 21 Oct 2020 13:13:42 -0500
Subject: [PATCH] Tracking electrons in central region

- Looking at "barrel" electrons
---
 .gitlab-ci.yml                             | 43 +----------
 tracking/central_electrons.sh              | 51 +++++++++++++
 tracking/options/tracker_reconstruction.py | 86 ++++++++++++++++++++++
 tracking/scripts/gen_central_electrons.cxx | 83 +++++++++++++++++++++
 tracking/tracking_config.yml               | 14 ++++
 5 files changed, 235 insertions(+), 42 deletions(-)
 create mode 100644 tracking/central_electrons.sh
 create mode 100644 tracking/options/tracker_reconstruction.py
 create mode 100644 tracking/scripts/gen_central_electrons.cxx
 create mode 100644 tracking/tracking_config.yml

diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index aa49676f..1ad1744a 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -6,23 +6,7 @@ default:
   artifacts:
     expire_in: 3 days
     paths:
-      - config/
       - results/
-        # - sim_output/
-        #    exclude:
-        #      - .git/
-        #      - datasets/.git/
-        # before_script:
-    #- git clone --depth 1 https://eicweb.phy.anl.gov/EIC/NPDet.git && mkdir NPDet/build && cd NPDet/build && cmake ../. -DCMAKE_INSTALL_PREFIX=/usr/local && make -j20 install && cd ../..
-    #- git clone --depth 1 https://eicweb.phy.anl.gov/EIC/eicd.git && mkdir eicd/build && cd eicd/build && cmake ../. -DCMAKE_INSTALL_PREFIX=/usr/local && make -j20 install && cd ../..
-    #- git clone --depth 1 https://eicweb.phy.anl.gov/EIC/juggler.git && mkdir juggler/build && cd juggler/build && cmake ../. -DCMAKE_INSTALL_PREFIX=/usr/local && make -j20 install && cd ../..
-    #- git clone --depth 1 https://eicweb.phy.anl.gov/EIC/detectors/topside.git && mkdir topside/build && cd topside/build && cmake ../. -DCMAKE_INSTALL_PREFIX=/usr/local && make -j20 install && cd ../..
-        #    - cd NPDet/build && cmake ../. -DCMAKE_INSTALL_PREFIX=/usr/local && make -j10 && make install
-
-        #before_script:
-        #  - git clone https://eicweb.phy.anl.gov/EIC/NPDet.git
-        #      #    - cd NPDet/build && cmake ../. -DCMAKE_INSTALL_PREFIX=/usr/local && make -j10 && make install
-        #      #    - cd ../.. 
 
 stages:
   - build
@@ -39,35 +23,10 @@ get_data:
     - mkdir -p results
     - mkdir -p config
 
-      #generate_config:
-      #  stage: config
-      #  needs: ["get_data"]
-      #  tags:
-      #    - silicon
-      #  script:
-      #    - mkdir -p config && ./bin/gen_ci_config -p test_ -i dummy > config/dummy_config.yml
-      #    - mkdir -p config && ./bin/gen_ci_config -p clustering_ -i clustering > config/clustering_config.yml
-      #
-      #dummy-pipeline:
-      #  stage: run
-      #  needs: ["generate_config"]
-      #  trigger:
-      #    include:
-      #      - artifact: config/dummy_config.yml
-      #        job: generate_config
-      #    strategy: depend
-      #
-      #clustering-pipeline:
-      #  stage: run
-      #  needs: ["generate_config"]
-      #  trigger:
-      #    include:
-      #      - artifact: config/clustering_config.yml
-      #        job: generate_config
-      #    strategy: depend
 
 include: 
   - local: 'ecal/ecal_config.yml'
+  - local: 'tracking/tracking_config.yml'
 
 
 final_report:
diff --git a/tracking/central_electrons.sh b/tracking/central_electrons.sh
new file mode 100644
index 00000000..2ce49b8d
--- /dev/null
+++ b/tracking/central_electrons.sh
@@ -0,0 +1,51 @@
+#!/bin/bash
+
+if [[ ! -n  "${JUGGLER_DETECTOR}" ]] ; then 
+  export JUGGLER_DETECTOR="topside"
+fi
+
+if [[ ! -n  "${JUGGLER_N_EVENTS}" ]] ; then 
+  export JUGGLER_N_EVENTS=100
+fi
+
+export JUGGLER_FILE_NAME_TAG="central_electrons"
+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
+mkdir ${JUGGLER_DETECTOR}/build
+pushd ${JUGGLER_DETECTOR}/build
+cmake ../. -DCMAKE_INSTALL_PREFIX=/usr/local && make -j30 install
+popd
+
+# generate the input events
+# note datasets is now only used to develop datasets.
+#git clone https://eicweb.phy.anl.gov/EIC/datasets.git datasets
+root -b -q "tracking/scripts/gen_central_electrons.cxx(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+
+pushd ${JUGGLER_DETECTOR}
+ls -l
+# run geant4 simulations
+npsim --runType batch \
+      -v WARNING \
+      --numberOfEvents ${JUGGLER_N_EVENTS} \
+      --compactFile ${JUGGLER_DETECTOR}.xml \
+      --inputFiles ../${JUGGLER_FILE_NAME_TAG}.hepmc \
+      --outputFile  ${JUGGLER_SIM_FILE}
+
+# Need to figure out how to pass file name to juggler from the commandline
+xenv -x /usr/local/Juggler.xenv gaudirun.py ../tracking/options/tracker_reconstruction.py
+ls -l
+popd
+
+pwd
+mkdir -p results
+cp topside/${JUGGLER_REC_FILE} results/.
+
diff --git a/tracking/options/tracker_reconstruction.py b/tracking/options/tracker_reconstruction.py
new file mode 100644
index 00000000..5db1281b
--- /dev/null
+++ b/tracking/options/tracker_reconstruction.py
@@ -0,0 +1,86 @@
+from Gaudi.Configuration import *
+
+from GaudiKernel.DataObjectHandleBase import DataObjectHandleBase
+from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
+
+detector_name = "topside"
+if "JUGGLER_DETECTOR" in os.environ :
+  detector_name = str(os.environ["JUGGLER_DETECTOR"])
+
+# todo add checks
+input_sim_file  = str(os.environ["JUGGLER_SIM_FILE"])
+output_rec_file = str(os.environ["JUGGLER_REC_FILE"])
+n_events = str(os.environ["JUGGLER_N_EVENTS"])
+
+geo_service  = GeoSvc("GeoSvc", detectors=["{}.xml".format(detector_name)])
+podioevent   = EICDataSvc("EventDataSvc", inputs=[input_sim_file], OutputLevel=DEBUG)
+
+from Configurables import PodioInput
+from Configurables import Jug__Digi__ExampleCaloDigi as ExampleCaloDigi
+from Configurables import Jug__Digi__UFSDTrackerDigi as UFSDTrackerDigi
+from Configurables import Jug__Reco__TrackerHitReconstruction as TrackerHitReconstruction
+from Configurables import Jug__Reco__TrackerSourceLinker as TrackerSourceLinker
+from Configurables import Jug__Reco__TrackParamTruthInit as TrackParamTruthInit
+from Configurables import Jug__Reco__TrackFindingAlgorithm as TrackFindingAlgorithm
+from Configurables import Jug__Reco__ParticlesFromTrackFit as ParticlesFromTrackFit
+from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
+from Configurables import Jug__Base__InputCopier_dd4pod__CalorimeterHitCollection_dd4pod__CalorimeterHitCollection_ as CalCopier
+from Configurables import Jug__Base__InputCopier_dd4pod__TrackerHitCollection_dd4pod__TrackerHitCollection_ as TrkCopier
+
+
+podioinput = PodioInput("PodioReader", 
+                        collections=["mcparticles","SiTrackerBarrelHits"])#, OutputLevel=DEBUG)
+
+## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
+copier = MCCopier("MCCopier", inputCollection="mcparticles", outputCollection="mcparticles2",OutputLevel=DEBUG) 
+trkcopier = TrkCopier("TrkCopier", inputCollection="SiTrackerBarrelHits", outputCollection="SiTrackerBarrelHits2",OutputLevel=DEBUG) 
+#caldigi = ExampleCaloDigi(inputHitCollection="FAEC_ShHits",outputHitCollection="RawFAECShowerHits")
+#ufsd_digi = UFSDTrackerDigi("ufsd_digi", inputHitCollection="SiVertexBarrelHits",outputHitCollection="VertexRawHits",timeResolution=8)
+ufsd_digi = UFSDTrackerDigi("ufsd_digi", inputHitCollection="SiTrackerBarrelHits",outputHitCollection="SiTrackerBarrelRawHits",timeResolution=8)
+trackpartruth = TrackParamTruthInit("trk_par_init",inputMCParticles="mcparticles",outputInitialTrackParameters="InitTrackParams")#,OutputLevel=DEBUG)
+trackerhit = TrackerHitReconstruction("trk_hit_reco",inputHitCollection="SiTrackerBarrelRawHits",outputHitCollection="TrackerBarrelRecHits",OutputLevel=DEBUG)
+sourcelinker = TrackerSourceLinker("trk_srclinker",inputHitCollection="TrackerBarrelRecHits",outputSourceLinks="BarrelTrackSourceLinks",OutputLevel=DEBUG)
+trk_find_alg = TrackFindingAlgorithm("trk_find_alg",inputSourceLinks="BarrelTrackSourceLinks",inputInitialTrackParameters= "InitTrackParams", outputTrajectories="trajectories",OutputLevel=DEBUG)
+parts_from_fit = ParticlesFromTrackFit("parts_from_fit",inputTrajectories="trajectories",outputParticles="ReconstructedParticles")#,OutputLevel=DEBUG)
+
+#types = []
+## this printout is useful to check that the type information is passed to python correctly
+#print("---------------------------------------\n")
+#print("---\n# List of input and output types by class")
+#for configurable in sorted([ PodioInput, EICDataSvc, PodioOutput,
+#                             TrackerHitReconstruction,ExampleCaloDigi, 
+#                             UFSDTrackerDigi, TrackerSourceLinker,
+#                             PodioOutput],
+#                           key=lambda c: c.getType()):
+#    print("\"{}\":".format(configurable.getType()))
+#    props = configurable.getDefaultProperties()
+#    for propname, prop in sorted(props.items()):
+#        print(" prop name: {}".format(propname))
+#        if isinstance(prop, DataObjectHandleBase):
+#            types.append(prop.type())
+#            print("  {}: \"{}\"".format(propname, prop.type()))
+#print("---")
+
+out = PodioOutput("out", filename=output_rec_file)
+out.outputCommands = ["keep *", 
+        "drop BarrelTrackSourceLinks", 
+        "drop InitTrackParams",
+        "drop trajectories",
+        "drop outputSourceLinks",
+        "drop outputInitialTrackParameters"
+        ]
+
+ApplicationMgr(
+    TopAlg = [podioinput, 
+              copier, trkcopier,
+              ufsd_digi, trackerhit, 
+              sourcelinker, trackpartruth, 
+              trk_find_alg, parts_from_fit,
+              out
+              ],
+    EvtSel = 'NONE',
+    EvtMax   = n_events,
+    ExtSvc = [podioevent,geo_service],
+    OutputLevel=DEBUG
+ )
+
diff --git a/tracking/scripts/gen_central_electrons.cxx b/tracking/scripts/gen_central_electrons.cxx
new file mode 100644
index 00000000..1f57715a
--- /dev/null
+++ b/tracking/scripts/gen_central_electrons.cxx
@@ -0,0 +1,83 @@
+#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(M_PI * (45 / 180.0));
+  double cos_theta_max = std::cos(125*M_PI);
+
+  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(0.0, 30.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;
+}
diff --git a/tracking/tracking_config.yml b/tracking/tracking_config.yml
new file mode 100644
index 00000000..b139788f
--- /dev/null
+++ b/tracking/tracking_config.yml
@@ -0,0 +1,14 @@
+tracking_central_electrons:
+  image: eicweb.phy.anl.gov:4567/eic/juggler/juggler:latest
+  tags:
+    - silicon
+  needs: ["get_data"]
+  timeout: 24 hours
+  artifacts:
+    expire_in: 20 weeks
+    paths:
+      - results/
+  stage: run
+  script:
+    - bash tracking/central_electrons.sh
+
-- 
GitLab