README.md

@@ -11,21 +11,20 @@ Reconstruction Benchmarks for the EIC
 Here we setup to use our local build of the `juggler` library.
 First set some environment variables.
 ```
-export JUGGLER_INSTALL_PREFIX=$HOME/stow/juggler
+export JUGGLER_INSTALL_PREFIX=/usr/local
 export JUGGLER_DETECTOR=athena   # athena is the default
 export BEAMLINE_CONFIG=ip6       # ip6 is the default
 ```
 
 
+
 ```
-git@eicweb.phy.anl.gov:EIC/benchmarks/reconstruction_benchmarks.git && cd reconstruction_benchmarks
+git clone https://eicweb.phy.anl.gov/EIC/benchmarks/reconstruction_benchmarks.git && cd reconstruction_benchmarks
 git clone https://eicweb.phy.anl.gov/EIC/benchmarks/common_bench.git setup
 source setup/bin/env.sh && ./setup/bin/install_common.sh
 source .local/bin/env.sh && build_detector.sh
 mkdir_local_data_link sim_output
-mkdir -p results
-mkdir -p config
-
+mkdir -p results config
 ```
 
 

benchmarks/imaging_ecal/config.yml

@@ -26,3 +26,24 @@ imaging_ecal_pion0:
   script:
     - bash benchmarks/imaging_ecal/run_imcal_pion0.sh -t imcal_barrel_pion0 -p "pion0" -n 100
 
+## SJJ removed as these are causing the CI to fail
+#imaging_ecal_pion_rejection:
+#  extends: .rec_benchmark
+#  timeout: 48 hours
+#  stage: run
+#  script:
+#    - compile_analyses.py imaging_ecal
+#    - bash benchmarks/imaging_ecal/run_emcal_barrel_pion_rej.sh -t emcal_barrel_pion_rej_electron -p "electron" -n 100
+#    - bash benchmarks/imaging_ecal/run_emcal_barrel_pion_rej.sh -t emcal_barrel_pion_rej_piminus -p "pion-" -n 100
+
+#imaging_ecal_pion_rejection:bench:
+#  extends: .rec_benchmark
+#  timeout: 1 hours
+#  stage: benchmarks2
+#  needs:
+#    - ["imaging_ecal_pion_rejection"]
+#  script:
+#    - ls -lhtR
+#    - compile_analyses.py imaging_ecal
+#    - root -b -q benchmarks/imaging_ecal/analysis/emcal_barrel_pion_rejection_analysis.cxx+
+#    #- bash run_pion_rej_analysis.sh

benchmarks/imaging_ecal/options/hybrid_cluster.py

@@ -19,11 +19,13 @@ from Configurables import Jug__Reco__ImagingClusterReco as ImagingClusterReco
 with open('config/emcal_barrel_calibration.json') as f:
     calib_data = json.load(f)['electron']
 
-print(calib_data)
+#print(calib_data)
 
 kwargs = dict()
 kwargs['img_sf'] = float(calib_data['sampling_fraction_img'])
 kwargs['scfi_sf'] = float(calib_data['sampling_fraction_scfi'])
+#kwargs['img_sf'] = 0.0060795
+#kwargs['scfi_sf'] = 0.10343
 
 # input arguments through environment variables
 kwargs['input'] = os.environ['CB_EMCAL_SIM_FILE']
@@ -56,14 +58,15 @@ imcaldaq = dict(
         pedestalMean=400,
         pedestalSigma=50)   # 50/32767*3 MeV ~ 5 keV
 
+
 imcaldigi = CalHitDigi('imcal_digi',
-        # OutputLevel=DEBUG,
+        OutputLevel=DEBUG,
         inputHitCollection='EcalBarrelHits',
         outputHitCollection='DigiEcalBarrelImagingHits',
         energyResolutions=[0., 0.02, 0.],
         **imcaldaq)
 imcalreco = ImagingPixelReco('imcal_reco',
-        # OutputLevel=DEBUG,
+        #OutputLevel=DEBUG,
         inputHitCollection=imcaldigi.outputHitCollection,
         outputHitCollection='RecoEcalBarrelImagingHits',
         readoutClass='EcalBarrelHits',
@@ -72,7 +75,7 @@ imcalreco = ImagingPixelReco('imcal_reco',
         **imcaldaq)
 
 imcalcluster = ImagingTopoCluster('imcal_cluster',
-        # OutputLevel=DEBUG,
+        #OutputLevel=DEBUG,
         inputHitCollection=imcalreco.outputHitCollection,
         outputProtoClusterCollection='EcalBarrelImagingProtoClusters',
         localDistXY=[2.*mm, 2*mm],
@@ -80,7 +83,7 @@ imcalcluster = ImagingTopoCluster('imcal_cluster',
         neighbourLayersRange=2,
         sectorDist=3.*cm)
 clusterreco = ImagingClusterReco('imcal_clreco',
-        # OutputLevel=DEBUG,
+        #OutputLevel=DEBUG,
         inputHitCollection=imcalcluster.inputHitCollection,
         inputProtoClusterCollection=imcalcluster.outputProtoClusterCollection,
         outputLayerCollection='EcalBarrelImagingClustersLayers',
@@ -112,7 +115,7 @@ scfi_barrel_reco = CalHitReco("scfi_barrel_reco",
 
 # merge hits in different layer (projection to local x-y plane)
 scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger",
-        # OutputLevel=DEBUG,
+        #OutputLevel=DEBUG,
         inputHitCollection=scfi_barrel_reco.outputHitCollection,
         outputHitCollection="EcalBarrelScFiGridReco",
         fields=["fiber"],
@@ -120,7 +123,7 @@ scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger",
         readoutClass="EcalBarrelScFiHits")
 
 scfi_barrel_cl = IslandCluster("scfi_barrel_cl",
-        # OutputLevel=DEBUG,
+        #OutputLevel=DEBUG,
         inputHitCollection=scfi_barrel_merger.outputHitCollection,
         outputProtoClusterCollection="EcalBarrelScFiProtoClusters",
         splitCluster=False,

benchmarks/imaging_ecal/run_emcal_barrel_pion_rej.sh

+#!/bin/bash
+
+print_env.sh
+
+function print_the_help {
+  echo "USAGE: ${0} -n <nevents> -t <nametag> -p <particle> "
+  echo "  OPTIONS: "
+  echo "    -n,--nevents     Number of events"
+  echo "    -t,--nametag     name tag"
+  echo "    -p,--particle    particle type"
+  echo "                     allowed types: pion0, pion+, pion-, kaon0, kaon+, kaon-, proton, neutron, electron, positron, photon"
+  exit
+}
+
+POSITIONAL=()
+while [[ $# -gt 0 ]]
+do
+  key="$1"
+
+  case $key in
+    -h|--help)
+      shift # past argument
+      print_the_help
+      ;;
+    -t|--nametag)
+      nametag="$2"
+      shift # past argument
+      shift # past value
+      ;;
+    -p|--particle)
+      particle="$2"
+      shift # past argument
+      shift # past value
+      ;;
+    -n|--nevents)
+      export CB_EMCAL_NUMEV="$2"
+      shift # past argument
+      shift # past value
+      ;;
+    *)    # unknown option
+      #POSITIONAL+=("$1") # save it in an array for later
+      echo "unknown option $1"
+      print_the_help
+      shift # past argument
+      ;;
+  esac
+done
+set -- "${POSITIONAL[@]}" # restore positional parameters
+
+export CB_EMCAL_COMPACT_PATH=${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml
+
+if [[ ! -n  "${CB_EMCAL_NUMEV}" ]] ; then
+  export CB_EMCAL_NUMEV=1000
+fi
+
+if [[ ! -n "${CB_EMCAL_IEV}" ]] ; then
+  export CB_EMCAL_IEV="0, 1"
+fi
+
+if [[ ! -n  "${CB_EMCAL_ENERGY}" ]] ; then
+  export CB_EMCAL_ENERGY=5.0
+fi
+
+if [[ ! -n "${CB_EMCAL_SAMP_FRAC}" ]] ; then
+  export CB_EMCAL_SAMP_FRAC=0.014
+fi
+
+export CB_EMCAL_NAME_TAG="${nametag}"
+export CB_EMCAL_GEN_FILE="${CB_EMCAL_NAME_TAG}.hepmc"
+
+export CB_EMCAL_SIM_FILE="sim_${CB_EMCAL_NAME_TAG}.root"
+export CB_EMCAL_REC_FILE="rec_${CB_EMCAL_NAME_TAG}.root"
+
+echo "CB_EMCAL_NUMEV = ${CB_EMCAL_NUMEV}"
+echo "CB_EMCAL_COMPACT_PATH = ${CB_EMCAL_COMPACT_PATH}"
+
+# Generate the input events
+python benchmarks/imaging_ecal/scripts/gen_particles.py ${CB_EMCAL_GEN_FILE} -n ${CB_EMCAL_NUMEV}\
+    --angmin 60 --angmax 120 --pmin 0.05 --pmax 18.0 --particles="${particle}"
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running script: generating input events"
+  exit 1
+fi
+
+ls -lh ${CB_EMCAL_GEN_FILE}
+
+# Run geant4 simulations
+npsim --runType batch \
+      -v WARNING \
+      --part.minimalKineticEnergy "0.5*MeV" \
+      --numberOfEvents ${CB_EMCAL_NUMEV} \
+      --compactFile ${CB_EMCAL_COMPACT_PATH} \
+      --inputFiles ${CB_EMCAL_GEN_FILE} \
+      --outputFile ${CB_EMCAL_SIM_FILE}
+
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running npdet"
+  exit 1
+fi
+
+echo "Simulation file"
+rootls -t "${CB_EMCAL_SIM_FILE}"
+
+# Directory for plots
+mkdir -p results
+
+CB_EMCAL_OPTION_DIR=benchmarks/imaging_ecal/options
+# Run Juggler
+xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
+    gaudirun.py ${CB_EMCAL_OPTION_DIR}/hybrid_cluster.py
+# gaudirun.py ${CB_EMCAL_OPTION_DIR}/imaging_topocluster.py
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running juggler"
+  exit 1
+fi
+
+# Plot clusters first
+FULL_CAL_SCRIPT_DIR=benchmarks/clustering/scripts
+python ${FULL_CAL_SCRIPT_DIR}/cluster_plots.py ${CB_EMCAL_SIM_FILE} ${CB_EMCAL_REC_FILE} -o results \
+    --collections "EcalBarrelImagingClusters, EcalBarrelScFiClusters"
+
+# check required python modules
+python -m pip install -r benchmarks/imaging_ecal/requirements.txt
+CB_EMCAL_SCRIPT_DIR=benchmarks/imaging_ecal/scripts
+IFS=',' read -ra ADDR <<< "$CB_EMCAL_IEV"
+for iev in "${ADDR[@]}"; do
+    if [[ $iev -ge "${CB_EMCAL_NUMEV}" ]] ; then
+        continue
+    fi
+
+    python ${CB_EMCAL_SCRIPT_DIR}/draw_cluster_layers.py \
+        ${CB_EMCAL_REC_FILE} -e ${iev} --topo-size=1.0 --compact=${CB_EMCAL_COMPACT_PATH} -o results/${particle}
+    if [[ "$?" -ne "0" ]] ; then
+      echo "ERROR running analysis script: draw_cluster_layers"
+      exit 1
+    fi
+
+    python ${CB_EMCAL_SCRIPT_DIR}/draw_cluster.py \
+        ${CB_EMCAL_REC_FILE} -e ${iev} --topo-size=2.0 --compact=${CB_EMCAL_COMPACT_PATH} -o results/${particle}
+    if [[ "$?" -ne "0" ]] ; then
+      echo "ERROR running analysis script: draw_cluster"
+      exit 1
+    fi
+done
+
+python ${CB_EMCAL_SCRIPT_DIR}/energy_profile.py \
+    ${CB_EMCAL_REC_FILE} --type=EM --energy=${CB_EMCAL_ENERGY} -o results/${particle} \
+    --save=results/profile.csv --color=royalblue
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running analysis script: energy_profile"
+  exit 1
+fi
+
+echo "Reconstruction File"
+rootls -t "${CB_EMCAL_REC_FILE}"
+
+root_filesize=$(stat --format=%s "${CB_EMCAL_REC_FILE}")
+if [[ "${CB_EMCAL_NUMEV}" -lt "500" ]] ; then
+  # file must be less than 10 MB to upload
+  if [[ "${root_filesize}" -lt "10000000" ]] ; then
+    cp ${CB_EMCAL_REC_FILE} results/.
+  fi
+fi
+
+pwd
+echo "Current directory"
+ls 
+
+#echo "Results directory"
+#ls benchmarks/imaging_ecal/results/
+
+
+

benchmarks/imaging_ecal/scripts/draw_cluster_layers.py

@@ -110,7 +110,7 @@ if __name__ == '__main__':
     else:
         dfc = df[(df['eta'] <= df['eta'].quantile(0.95)) & (df['eta'] >= df['eta'].quantile(0.05)) &
                  (df['phi'] <= df['phi'].quantile(0.95)) & (df['phi'] >= df['phi'].quantile(0.05))]
-        vec = np.average(dfc[['x', 'y', 'z']].values, axis=0, weights=dfc['energy'].values)
+        vec = np.average(dfc[['x', 'y', 'z']].values, axis=0, weights= dfc['energy'].values if np.sum(dfc['energy'].values) > 0 else None)
     vec = vec/np.linalg.norm(vec)
 
     # particle line from (0, 0, 0) to the inner Ecal surface

benchmarks/track_finding/options/track_reconstruction.py

@@ -5,11 +5,17 @@ from GaudiKernel import SystemOfUnits as units
 
 detector_name = "topside"
 if "JUGGLER_DETECTOR" in os.environ :
-  detector_name = str(os.environ["JUGGLER_DETECTOR"])
+    detector_name = str(os.environ["JUGGLER_DETECTOR"])
 
 detector_path = ""
 if "DETECTOR_PATH" in os.environ :
-  detector_path = str(os.environ["DETECTOR_PATH"])
+    detector_path = str(os.environ["DETECTOR_PATH"])
+
+detector_version = 'default'
+if "JUGGLER_DETECTOR_VERSION" in os.environ:
+    env_version = str(os.environ["JUGGLER_DETECTOR_VERSION"])
+    if 'acadia' in env_version:
+        detector_version = 'acadia'
 
 # todo add checks
 input_sim_file  = str(os.environ["JUGGLER_SIM_FILE"])
@@ -40,8 +46,13 @@ from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
 
 algorithms = [ ]
 
+input_collections = ['mcparticles','TrackerEndcapHits','TrackerBarrelHits','VertexBarrelHits','GEMTrackerEndcapHits']
+if 'acadia' in detector_version:
+    input_collections.append('VertexEndcapHits')
+else:
+    input_collections.append('MPGDTrackerBarrelHits')
 podioinput = PodioInput("PodioReader", 
-                        collections=["mcparticles","TrackerEndcapHits","TrackerBarrelHits","VertexBarrelHits","VertexEndcapHits","GEMTrackerEndcapHits"])
+                        collections=input_collections)#, OutputLevel=DEBUG)
 algorithms.append( podioinput )
 
 trk_b_digi = TrackerDigi("trk_b_digi", 
@@ -61,11 +72,18 @@ vtx_b_digi = TrackerDigi("vtx_b_digi",
         timeResolution=8)
 algorithms.append( vtx_b_digi )
 
-vtx_ec_digi = TrackerDigi("vtx_ec_digi", 
-        inputHitCollection="VertexEndcapHits",
-        outputHitCollection="VertexEndcapRawHits",
-        timeResolution=8)
-algorithms.append( vtx_ec_digi )
+if 'acadia' in detector_version:
+    vtx_ec_digi = TrackerDigi("vtx_ec_digi", 
+            inputHitCollection="VertexEndcapHits",
+            outputHitCollection="VertexEndcapRawHits",
+            timeResolution=8)
+    algorithms.append( vtx_ec_digi )
+else:
+    mm_b_digi = TrackerDigi("mm_b_digi", 
+            inputHitCollection="MPGDTrackerBarrelHits",
+            outputHitCollection="MPGDTrackerBarrelRawHits",
+            timeResolution=8)
+    algorithms.append( mm_b_digi )
 
 gem_ec_digi = TrackerDigi("gem_ec_digi",
         inputHitCollection="GEMTrackerEndcapHits",
@@ -89,23 +107,35 @@ vtx_b_reco = TrackerHitReconstruction("vtx_b_reco",
         outputHitCollection="VertexBarrelRecHits")
 algorithms.append( vtx_b_reco )
 
-vtx_ec_reco = TrackerHitReconstruction("vtx_ec_reco",
-        inputHitCollection = vtx_ec_digi.outputHitCollection,
-        outputHitCollection="VertexEndcapRecHits")
-algorithms.append( vtx_ec_reco )
+if 'acadia' in detector_version:
+    vtx_ec_reco = TrackerHitReconstruction("vtx_ec_reco",
+            inputHitCollection = vtx_ec_digi.outputHitCollection,
+            outputHitCollection="VertexEndcapRecHits")
+    algorithms.append( vtx_ec_reco )
+else:
+    mm_b_reco = TrackerHitReconstruction("mm_b_reco",
+            inputHitCollection = mm_b_digi.outputHitCollection,
+            outputHitCollection="MPGDTrackerBarrelRecHits")
+    algorithms.append( mm_b_reco )
+
 
 gem_ec_reco = TrackerHitReconstruction("gem_ec_reco",
         inputHitCollection=gem_ec_digi.outputHitCollection,
         outputHitCollection="GEMTrackerEndcapRecHits")
 algorithms.append(gem_ec_reco)
 
+input_tracking_hits = [
+    str(trk_b_reco.outputHitCollection),
+    str(trk_ec_reco.outputHitCollection),
+    str(vtx_b_reco.outputHitCollection),
+    str(gem_ec_reco.outputHitCollection) ]
+if 'acadia' in detector_version:
+    input_tracking_hits.append(str(vtx_ec_reco.outputHitCollection))
+else:
+    input_tracking_hits.append(str(mm_b_reco.outputHitCollection))
+
 trk_hit_col = TrackingHitsCollector("trk_hit_col",
-        inputTrackingHits=[
-            str(trk_b_reco.outputHitCollection),
-            str(trk_ec_reco.outputHitCollection),
-            str(vtx_b_reco.outputHitCollection),
-            str(vtx_ec_reco.outputHitCollection),
-            str(gem_ec_reco.outputHitCollection) ],
+        inputTrackingHits=input_tracking_hits,
         trackingHits="trackingHits",
         OutputLevel=DEBUG)
 algorithms.append( trk_hit_col )

benchmarks/track_fitting/options/track_reconstruction.py

@@ -5,11 +5,17 @@ from GaudiKernel import SystemOfUnits as units
 
 detector_name = "athena"
 if "JUGGLER_DETECTOR" in os.environ :
-  detector_name = str(os.environ["JUGGLER_DETECTOR"])
+    detector_name = str(os.environ["JUGGLER_DETECTOR"])
 
 detector_path = ""
 if "DETECTOR_PATH" in os.environ :
-  detector_path = str(os.environ["DETECTOR_PATH"])
+    detector_path = str(os.environ["DETECTOR_PATH"])
+
+detector_version = 'default'
+if "JUGGLER_DETECTOR_VERSION" in os.environ:
+    env_version = str(os.environ["JUGGLER_DETECTOR_VERSION"])
+    if 'acadia' in env_version:
+        detector_version = 'acadia'
 
 # todo add checks
 input_sim_file  = str(os.environ["JUGGLER_SIM_FILE"])
@@ -40,8 +46,13 @@ from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
 
 algorithms = [ ]
 
+input_collections = ['mcparticles','TrackerEndcapHits','TrackerBarrelHits','VertexBarrelHits','GEMTrackerEndcapHits']
+if 'acadia' in detector_version:
+    input_collections.append('VertexEndcapHits')
+else:
+    input_collections.append('MPGDTrackerBarrelHits')
 podioinput = PodioInput("PodioReader", 
-                        collections=["mcparticles","TrackerEndcapHits","TrackerBarrelHits","VertexBarrelHits","VertexEndcapHits","GEMTrackerEndcapHits"])#, OutputLevel=DEBUG)
+                        collections=input_collections)#, OutputLevel=DEBUG)
 algorithms.append( podioinput )
 
 trk_b_digi = TrackerDigi("trk_b_digi", 
@@ -61,11 +72,18 @@ vtx_b_digi = TrackerDigi("vtx_b_digi",
         timeResolution=8)
 algorithms.append( vtx_b_digi )
 
-vtx_ec_digi = TrackerDigi("vtx_ec_digi", 
-        inputHitCollection="VertexEndcapHits",
-        outputHitCollection="VertexEndcapRawHits",
-        timeResolution=8)
-algorithms.append( vtx_ec_digi )
+if 'acadia' in detector_version:
+    vtx_ec_digi = TrackerDigi("vtx_ec_digi", 
+            inputHitCollection="VertexEndcapHits",
+            outputHitCollection="VertexEndcapRawHits",
+            timeResolution=8)
+    algorithms.append( vtx_ec_digi )
+else:
+    mm_b_digi = TrackerDigi("mm_b_digi", 
+            inputHitCollection="MPGDTrackerBarrelHits",
+            outputHitCollection="MPGDTrackerBarrelRawHits",
+            timeResolution=8)
+    algorithms.append( mm_b_digi )
 
 gem_ec_digi = TrackerDigi("gem_ec_digi",
         inputHitCollection="GEMTrackerEndcapHits",
@@ -89,24 +107,36 @@ vtx_b_reco = TrackerHitReconstruction("vtx_b_reco",
         outputHitCollection="VertexBarrelRecHits")
 algorithms.append( vtx_b_reco )
 
-vtx_ec_reco = TrackerHitReconstruction("vtx_ec_reco",
-        inputHitCollection = vtx_ec_digi.outputHitCollection,
-        outputHitCollection="VertexEndcapRecHits")
-algorithms.append( vtx_ec_reco )
+if 'acadia' in detector_version:
+    vtx_ec_reco = TrackerHitReconstruction("vtx_ec_reco",
+            inputHitCollection = vtx_ec_digi.outputHitCollection,
+            outputHitCollection="VertexEndcapRecHits")
+    algorithms.append( vtx_ec_reco )
+else:
+    mm_b_reco = TrackerHitReconstruction("mm_b_reco",
+            inputHitCollection = mm_b_digi.outputHitCollection,
+            outputHitCollection="MPGDTrackerBarrelRecHits")
+    algorithms.append( mm_b_reco )
 
 gem_ec_reco = TrackerHitReconstruction("gem_ec_reco",
         inputHitCollection=gem_ec_digi.outputHitCollection,
         outputHitCollection="GEMTrackerEndcapRecHits")
 algorithms.append(gem_ec_reco)
 
+input_tracking_hits = [
+    str(trk_b_reco.outputHitCollection),
+    str(trk_ec_reco.outputHitCollection),
+    str(vtx_b_reco.outputHitCollection),
+    str(gem_ec_reco.outputHitCollection) ]
+if 'acadia' in detector_version:
+    input_tracking_hits.append(str(vtx_ec_reco.outputHitCollection))
+else:
+    input_tracking_hits.append(str(mm_b_reco.outputHitCollection))
+
 trk_hit_col = TrackingHitsCollector("trk_hit_col",
-        inputTrackingHits=[
-            str(vtx_b_reco.outputHitCollection),
-            str(vtx_ec_reco.outputHitCollection),
-            str(trk_b_reco.outputHitCollection),
-            str(trk_ec_reco.outputHitCollection),
-            str(gem_ec_reco.outputHitCollection) ],
-        trackingHits="trackingHits")
+        inputTrackingHits=input_tracking_hits,
+        trackingHits="trackingHits",
+        OutputLevel=DEBUG)
 algorithms.append( trk_hit_col )
 
 # Hit Source linker 

benchmarks/tracking/options/track_reconstruction.py

@@ -5,18 +5,36 @@ from GaudiKernel import SystemOfUnits as units
 
 detector_name = "topside"
 if "JUGGLER_DETECTOR" in os.environ :
-  detector_name = str(os.environ["JUGGLER_DETECTOR"])
+    detector_name = str(os.environ["JUGGLER_DETECTOR"])
 
 detector_path = ""
 if "DETECTOR_PATH" in os.environ :
-  detector_path = str(os.environ["DETECTOR_PATH"])
+    detector_path = str(os.environ["DETECTOR_PATH"])
+
+detector_version = 'default'
+if "JUGGLER_DETECTOR_VERSION" in os.environ:
+    env_version = str(os.environ["JUGGLER_DETECTOR_VERSION"])
+    if 'acadia' in env_version:
+        detector_version = 'acadia'
 
 # 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_path,detector_name)], OutputLevel=WARNING)
+## note: old version of material map is called material-maps.XXX, new version is materials-map.XXX
+##       these names are somewhat inconsistent, and should probably all be renamed to 'material-map.XXX' 
+##       FIXME
+if detector_version == 'acadia':
+    geo_service = GeoSvc("GeoSvc",
+                         detectors=["{}/{}.xml".format(detector_path,detector_name)],
+                         materials="config/material-maps.json",
+                         OutputLevel=WARNING)
+else:
+    geo_service = GeoSvc("GeoSvc",
+                         detectors=["{}/{}.xml".format(detector_path,detector_name)],
+                         materials="calibrations/materials-map.cbor",
+                         OutputLevel=WARNING)
 podioevent   = EICDataSvc("EventDataSvc", inputs=[input_sim_file], OutputLevel=WARNING)
 
 from Configurables import PodioInput
@@ -33,15 +51,20 @@ from Configurables import Jug__Reco__TrackParamVertexClusterInit as TrackParamVe
 
 from Configurables import Jug__Reco__TrackFindingAlgorithm as TrackFindingAlgorithm
 from Configurables import Jug__Reco__ParticlesFromTrackFit as ParticlesFromTrackFit
-from Configurables import Jug__Reco__TrajectoryFromTrackFit as TrajectoryFromTrackFit
+#from Configurables import Jug__Reco__TrajectoryFromTrackFit as TrajectoryFromTrackFit
 
 from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
 
 
 algorithms = [ ]
 
+input_collections = ['mcparticles','TrackerEndcapHits','TrackerBarrelHits','VertexBarrelHits','GEMTrackerEndcapHits']
+if 'acadia' in detector_version:
+    input_collections.append('VertexEndcapHits')
+else:
+    input_collections.append('MPGDTrackerBarrelHits')
 podioinput = PodioInput("PodioReader", 
-                        collections=["mcparticles","TrackerEndcapHits","TrackerBarrelHits","VertexBarrelHits","VertexEndcapHits","GEMTrackerEndcapHits"])#, OutputLevel=DEBUG)
+                        collections=input_collections)#, OutputLevel=DEBUG)
 algorithms.append( podioinput )
 
 trk_b_digi = TrackerDigi("trk_b_digi", 
@@ -61,11 +84,18 @@ vtx_b_digi = TrackerDigi("vtx_b_digi",
         timeResolution=8)
 algorithms.append( vtx_b_digi )
 
-vtx_ec_digi = TrackerDigi("vtx_ec_digi", 
-        inputHitCollection="VertexEndcapHits",
-        outputHitCollection="VertexEndcapRawHits",
-        timeResolution=8)
-algorithms.append( vtx_ec_digi )
+if 'acadia' in detector_version:
+    vtx_ec_digi = TrackerDigi("vtx_ec_digi", 
+            inputHitCollection="VertexEndcapHits",
+            outputHitCollection="VertexEndcapRawHits",
+            timeResolution=8)
+    algorithms.append( vtx_ec_digi )
+else:
+    mm_b_digi = TrackerDigi("mm_b_digi", 
+            inputHitCollection="MPGDTrackerBarrelHits",
+            outputHitCollection="MPGDTrackerBarrelRawHits",
+            timeResolution=8)
+    algorithms.append( mm_b_digi )
 
 gem_ec_digi = TrackerDigi("gem_ec_digi",
         inputHitCollection="GEMTrackerEndcapHits",
@@ -89,25 +119,35 @@ vtx_b_reco = TrackerHitReconstruction("vtx_b_reco",
         outputHitCollection="VertexBarrelRecHits")
 algorithms.append( vtx_b_reco )
 
-vtx_ec_reco = TrackerHitReconstruction("vtx_ec_reco",
-        inputHitCollection = vtx_ec_digi.outputHitCollection,
-        outputHitCollection="VertexEndcapRecHits")
-algorithms.append( vtx_ec_reco )
+if 'acadia' in detector_version:
+    vtx_ec_reco = TrackerHitReconstruction("vtx_ec_reco",
+            inputHitCollection = vtx_ec_digi.outputHitCollection,
+            outputHitCollection="VertexEndcapRecHits")
+    algorithms.append( vtx_ec_reco )
+else:
+    mm_b_reco = TrackerHitReconstruction("mm_b_reco",
+            inputHitCollection = mm_b_digi.outputHitCollection,
+            outputHitCollection="MPGDTrackerBarrelRecHits")
+    algorithms.append( mm_b_reco )
 
 gem_ec_reco = TrackerHitReconstruction("gem_ec_reco",
         inputHitCollection=gem_ec_digi.outputHitCollection,
         outputHitCollection="GEMTrackerEndcapRecHits")
 algorithms.append(gem_ec_reco)
 
+input_tracking_hits = [
+    str(trk_b_reco.outputHitCollection),
+    str(trk_ec_reco.outputHitCollection),
+    str(vtx_b_reco.outputHitCollection),
+    str(gem_ec_reco.outputHitCollection) ]
+if 'acadia' in detector_version:
+    input_tracking_hits.append(str(vtx_ec_reco.outputHitCollection))
+else:
+    input_tracking_hits.append(str(mm_b_reco.outputHitCollection))
+
 trk_hit_col = TrackingHitsCollector("trk_hit_col",
-        inputTrackingHits=[
-            str(trk_b_reco.outputHitCollection),
-            str(trk_ec_reco.outputHitCollection),
-            str(vtx_b_reco.outputHitCollection),
-            str(vtx_ec_reco.outputHitCollection),
-            str(gem_ec_reco.outputHitCollection) ],
-        trackingHits="trackingHits",
-        OutputLevel=DEBUG)
+        inputTrackingHits=input_tracking_hits,
+        trackingHits="trackingHits")
 algorithms.append( trk_hit_col )
 
 # Hit Source linker 
@@ -115,23 +155,21 @@ sourcelinker = TrackerSourceLinker("sourcelinker",
         inputHitCollection=trk_hit_col.trackingHits,
         outputSourceLinks="TrackSourceLinks",
         outputMeasurements="TrackMeasurements",
-        OutputLevel=DEBUG)
+        OutputLevel=VERBOSE)
 algorithms.append( sourcelinker )
 
 ## Track param init
 truth_trk_init = TrackParamTruthInit("truth_trk_init",
         inputMCParticles="mcparticles",
         outputInitialTrackParameters="InitTrackParams")
-        #OutputLevel=DEBUG)
 algorithms.append( truth_trk_init )
 
 # Tracking algorithms
 trk_find_alg = TrackFindingAlgorithm("trk_find_alg",
         inputSourceLinks = sourcelinker.outputSourceLinks,
         inputMeasurements = sourcelinker.outputMeasurements,
-        inputInitialTrackParameters= "InitTrackParams",#"InitTrackParamsFromClusters", 
+        inputInitialTrackParameters= "InitTrackParams",
         outputTrajectories="trajectories")
-        #OutputLevel=DEBUG)
 algorithms.append( trk_find_alg )
 
 parts_from_fit = ParticlesFromTrackFit("parts_from_fit",
@@ -141,29 +179,11 @@ parts_from_fit = ParticlesFromTrackFit("parts_from_fit",
         #OutputLevel=DEBUG)
 algorithms.append( parts_from_fit )
 
-trajs_from_fit = TrajectoryFromTrackFit("trajs_from_fit",
-inputTrajectories = trk_find_alg.outputTrajectories,
-outputTrajectoryParameters = "outputTrajectoryParameters")
+#trajs_from_fit = TrajectoryFromTrackFit("trajs_from_fit",
+#inputTrajectories = trk_find_alg.outputTrajectories,
+#outputTrajectoryParameters = "outputTrajectoryParameters")
         #OutputLevel=DEBUG)
-algorithms.append(trajs_from_fit)
-
-#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, DataHandleBase):
-#            types.append(prop.type())
-#            print("  {}: \"{}\"".format(propname, prop.type()))
-#print("---")
+#algorithms.append(trajs_from_fit)
 
 out = PodioOutput("out", filename=output_rec_file)
 out.outputCommands = ["keep *", 
@@ -183,4 +203,3 @@ ApplicationMgr(
     ExtSvc = [podioevent,geo_service],
     OutputLevel=WARNING
  )
-

benchmarks/tracking/scripts/tracking_performance.py

@@ -186,7 +186,7 @@ if __name__ == '__main__':
                              weights=np.repeat(1./float(rec_p.shape[0]), rec_p.shape[0]), ec='k')
     nbins = hbins.shape[0] - 1
     ax.set_ylabel('Normalized Counts / {:d} Bins'.format(nbins), fontsize=20)
-    ax.set_xlabel('$\delta p$ (%)', fontsize=20)
+    ax.set_xlabel(r'$\delta$ p (%)', fontsize=20)
 
     # theta resolution
     ax = axs.flat[2]

bin/get_calibrations

@@ -4,7 +4,7 @@ branch=${1:-master}
 detector_benchmarks=https://eicweb.phy.anl.gov/EIC/benchmarks/detector_benchmarks/-/jobs/artifacts/${branch}/raw/
 
 mkdir -p config
-for i in results/emcal_barrel_calibration.json ; do
+for i in results/emcal_barrel_calibration.json results/material-maps.json ; do
   curl --fail -sL ${detector_benchmarks}/${i}?job=deploy_results --output config/$(basename ${i})
   if [[ "$?" -ne "0" ]] ; then
     if find ${DETECTOR_PATH} -name $(basename $i) ; then