Full/reference/best-effort reconstruction

.gitlab-ci.yml

@@ -53,6 +53,7 @@ common:detector:
     - ln -s "${LOCAL_DATA_PATH}/sim_output" sim_output
     - ln -s "${LOCAL_DATA_PATH}/datasets/data" data
     - ls -lrtha
+  interruptible: true
 
 include: 
   - local: 'benchmarks/ecal/config.yml'
@@ -61,11 +62,12 @@ include:
   - local: 'benchmarks/rich/config.yml'
   - local: 'benchmarks/imaging_ecal/config.yml'
   - local: 'benchmarks/imaging_shower_ML/config.yml'
+  - local: 'benchmarks/full/config.yml'
 
 
 final_report:
   stage: finish
-  needs: ["ecal_collect", "tracking_central_electrons","clustering:results"]
+  needs: ["ecal_collect", "tracking_central_electrons", "clustering:results", "full:results"]
   script:
     # disabled while we address ACTS issues
     #- mkdir -p results/views && cd results/views && bash ../../bin/download_views

benchmarks/full/config.yml

+full:reco:
+  extends: .rec_benchmark
+  stage: run
+  timeout: 24 hours
+  script:
+    - bash benchmarks/full/full_reconstruction.sh -t full_${REGIONS}_${PARTICLES} -n 100 -p "${PARTICLES}" --pmin 1 --pmax 20 --region "${REGIONS}"
+  parallel:
+    matrix:
+      - REGIONS: ["forward", "barrel", "backward"]
+        PARTICLES: ["electron", "pion+", "pion-", "kaon+", "kaon-", "proton", "neutron", "positron", "photon"]
+        #PARTICLES: ["electron", "pion0", "pion+", "pion-", "kaon0S", "kaon0L", "kaon+", "kaon-", "proton", "neutron", "positron", "photon"]
+
+full:results:
+  extends: .rec_benchmark
+  stage: collect
+  needs: ["full:reco"]
+  script:
+    - ls -lrth
+      #- python dvcs/scripts/merge_results.py
+

benchmarks/full/full_reconstruction.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"
+  echo "    --pmin           minimum particle momentum (GeV)"
+  echo "    --pmax           maximum particle momentum (GeV)"
+  echo "    --angmin         minimum particle angle (deg)"
+  echo "    --angmax         maximum particle angle (deg)"
+  echo "    --region         detector region"
+  echo "                     allowed regions: farbackward, backward, barrel, forward, farforward"
+  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 JUGGLER_N_EVENTS="$2"
+      shift # past argument
+      shift # past value
+      ;;
+    --pmin)
+      export PMIN="$2"
+      shift
+      shift
+      ;;
+    --pmax)
+      export PMAX="$2"
+      shift
+      shift
+      ;;
+    --angmin)
+      export ANGMIN="$2"
+      shift
+      shift
+      ;;
+    --angmax)
+      export ANGMAX="$2"
+      shift
+      shift
+      ;;
+    --region)
+      export REGION="$2"
+      shift
+      shift
+      ;;
+    *)    # 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
+
+if [[ ! -n  "${JUGGLER_N_EVENTS}" ]] ; then
+  export JUGGLER_N_EVENTS=1000
+fi
+
+if [[ ! -n  "${PMIN}" ]] ; then
+  export PMIN=1.0
+fi
+
+if [[ ! -n  "${PMAX}" ]] ; then
+  export PMAX=10.0
+fi
+
+if [[ ! -n  "${ANGMIN}" ]] ; then
+  export ANGMIN=3
+fi
+
+if [[ ! -n  "${ANGMAX}" ]] ; then
+  export ANGMAX=177
+fi
+
+if [[ ! -n  "${REGION}" ]] ; then
+  export REGION="ignore"
+fi
+
+if [[ ! -n "${JUGGLER_DETECTOR}" ]] ; then
+  export JUGGLER_DETECTOR="athena"
+fi
+
+if [[ ! -n "${JUGGLER_DETECTOR_PATH}" ]] ; then
+  export JUGGLER_DETECTOR_PATH="/opt/detector/share/athena"
+fi
+
+if [[ ! -n "${JUGGLER_INSTALL_PREFIX}" ]] ; then
+  export JUGGLER_INSTALL_PREFIX="/usr/local"
+fi
+
+compact_path=${JUGGLER_DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml
+export JUGGLER_FILE_NAME_TAG="${nametag}"
+export JUGGLER_GEN_FILE="gen_${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_PATH = ${JUGGLER_DETECTOR_PATH}"
+echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
+
+# Generate the input events
+if [ ! -f ${JUGGLER_GEN_FILE} ] ; then
+  python benchmarks/full/scripts/gen_particles.py ${JUGGLER_GEN_FILE} -n ${JUGGLER_N_EVENTS}\
+    --angmin ${ANGMIN} --angmax ${ANGMAX} --phmin 0 --phmax 360 --pmin ${PMIN} --pmax ${PMAX}\
+    --regions="${REGION}" --particles="${particle}"
+fi
+
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running script: generating input events"
+  exit 1
+fi
+
+ls -lh ${JUGGLER_GEN_FILE}
+
+# Run geant4 simulations
+if [ ! -f ${JUGGLER_SIM_FILE} ] ; then
+  npsim --runType batch \
+        -v WARNING \
+        --part.minimalKineticEnergy "1*TeV" \
+        --numberOfEvents ${JUGGLER_N_EVENTS} \
+        --compactFile ${compact_path} \
+        --inputFiles ${JUGGLER_GEN_FILE} \
+        --outputFile ${JUGGLER_SIM_FILE}
+fi
+
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running npdet"
+  exit 1
+fi
+
+rootls -t "${JUGGLER_SIM_FILE}"
+
+# Directory for results
+mkdir -p results
+
+OPTION_DIR=benchmarks/full/options
+SCRIPT_DIR=benchmarks/full/scripts
+
+# Run Juggler
+xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
+    gaudirun.py ${OPTION_DIR}/full_reconstruction.py
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running reconstruction (juggler)"
+  exit 1
+fi
+
+# Run analysis scripts
+#python ${SCRIPT_DIR}/full_reconstruction.py ${JUGGLER_SIM_FILE} ${JUGGLER_REC_FILE} -o results
+
+root_filesize=$(stat --format=%s "${JUGGLER_REC_FILE}")
+if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then
+  # file must be less than 10 MB to upload
+  if [[ "${root_filesize}" -lt "10000000" ]] ; then
+    cp ${JUGGLER_REC_FILE} results/.
+  fi
+fi
+

benchmarks/full/options/full_reconstruction.py

+from Gaudi.Configuration import *
+
+from GaudiKernel.DataObjectHandleBase import DataObjectHandleBase
+from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
+from GaudiKernel import SystemOfUnits as units
+from GaudiKernel.SystemOfUnits import MeV, GeV, mm, cm, mrad
+
+detector_name = "athena"
+if "JUGGLER_DETECTOR" in os.environ :
+  detector_name = str(os.environ["JUGGLER_DETECTOR"])
+
+detector_path = ""
+if "JUGGLER_DETECTOR_PATH" in os.environ :
+  detector_path = str(os.environ["JUGGLER_DETECTOR_PATH"])
+
+compact_path = os.path.join(detector_path, detector_name)
+
+# RICH reconstruction
+qe_data = [(1.0, 0.25), (7.5, 0.25),]
+
+# CAL reconstruction
+# get sampling fractions from system environment variable, 1.0 by default
+ci_ecal_sf = float(os.environ.get("CI_ECAL_SAMP_FRAC", 0.253))
+cb_ecal_sf = float(os.environ.get("CB_ECAL_SAMP_FRAC", 0.01324))
+cb_hcal_sf = float(os.environ.get("CB_HCAL_SAMP_FRAC", 0.038))
+ci_hcal_sf = float(os.environ.get("CI_HCAL_SAMP_FRAC", 0.025))
+ce_hcal_sf = float(os.environ.get("CE_HCAL_SAMP_FRAC", 0.025))
+
+# input and output
+input_sims = [f.strip() for f in str.split(os.environ["JUGGLER_SIM_FILE"], ",") if f.strip()]
+output_rec = str(os.environ["JUGGLER_REC_FILE"])
+n_events = int(os.environ["JUGGLER_N_EVENTS"])
+
+# geometry service
+geo_service = GeoSvc("GeoSvc", detectors=["{}.xml".format(compact_path)], OutputLevel=WARNING)
+# data service
+podioevent = EICDataSvc("EventDataSvc", inputs=input_sims, OutputLevel=WARNING)
+
+# juggler components
+from Configurables import PodioInput
+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
+from Configurables import Jug__Base__InputCopier_dd4pod__PhotoMultiplierHitCollection_dd4pod__PhotoMultiplierHitCollection_ as PMTCopier
+
+from Configurables import Jug__Base__MC2DummyParticle as MC2DummyParticle
+
+from Configurables import Jug__Digi__PhotoMultiplierDigi as PhotoMultiplierDigi
+from Configurables import Jug__Digi__CalorimeterHitDigi as CalHitDigi
+from Configurables import Jug__Digi__ExampleCaloDigi as ExampleCaloDigi
+from Configurables import Jug__Digi__UFSDTrackerDigi as TrackerDigi
+from Configurables import Jug__Digi__EMCalorimeterDigi as EMCalorimeterDigi
+
+from Configurables import Jug__Reco__TrackerHitReconstruction as TrackerHitReconstruction
+from Configurables import Jug__Reco__TrackerSourceLinker as TrackerSourceLinker
+from Configurables import Jug__Reco__TrackerSourcesLinker as TrackerSourcesLinker
+#from Configurables import Jug__Reco__TrackingHitsSourceLinker as TrackingHitsSourceLinker
+from Configurables import Jug__Reco__TrackParamTruthInit as TrackParamTruthInit
+from Configurables import Jug__Reco__TrackParamClusterInit as TrackParamClusterInit
+from Configurables import Jug__Reco__TrackParamVertexClusterInit as TrackParamVertexClusterInit
+from Configurables import Jug__Reco__TrackFindingAlgorithm as TrackFindingAlgorithm
+from Configurables import Jug__Reco__ParticlesFromTrackFit as ParticlesFromTrackFit
+
+from Configurables import Jug__Reco__CalorimeterHitReco as CalHitReco
+from Configurables import Jug__Reco__CalorimeterHitsMerger as CalHitsMerger
+from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
+
+from Configurables import Jug__Reco__ImagingPixelReco as ImCalPixelReco
+from Configurables import Jug__Reco__ImagingTopoCluster as ImagingCluster
+
+from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
+from Configurables import Jug__Reco__ImagingClusterReco as ImagingClusterReco
+
+from Configurables import Jug__Reco__PhotoMultiplierReco as PhotoMultiplierReco
+from Configurables import Jug__Reco__PhotoRingClusters as PhotoRingClusters
+
+from Configurables import Jug__Reco__EMCalReconstruction as EMCalReconstruction
+
+from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
+
+# branches needed from simulation root file
+sim_coll = [
+    "mcparticles",
+    "EcalEndcapNHits",
+    "EcalEndcapPHits",
+    "EcalBarrelHits",
+    "HcalBarrelHits",
+    "HcalHadronEndcapHits",
+    "HcalElectronEndcapHits",
+    "TrackerEndcapHits",
+    "TrackerBarrelHits",
+    "VertexBarrelHits",
+    "VertexEndcapHits",
+    "ForwardRICHHits"
+]
+
+# input and output
+podin = PodioInput("PodioReader", collections=sim_coll)
+podout = PodioOutput("out", filename=output_rec)
+
+## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
+copier = MCCopier("MCCopier", 
+        inputCollection="mcparticles", 
+        outputCollection="mcparticles2") 
+trkcopier = TrkCopier("TrkCopier", 
+        inputCollection="TrackerBarrelHits", 
+        outputCollection="TrackerBarrelHits2") 
+pmtcopier = PMTCopier("PMTCopier",
+        inputCollection="ForwardRICHHits",
+        outputCollection="ForwardRICHHits2")
+
+# Dummy reconstruction
+dummy = MC2DummyParticle("MC2Dummy",
+        inputCollection="mcparticles",
+        outputCollection="ReconstructedParticles",
+        smearing = 0.01)
+
+# Simple clusters
+ecal_digi = EMCalorimeterDigi("ecal_digi",
+        inputHitCollection="EcalBarrelHits",
+        outputHitCollection="EcalBarrelHitsSimpleDigi")
+
+ecal_reco = EMCalReconstruction("ecal_reco",
+        inputHitCollection="EcalBarrelHitsSimpleDigi",
+        outputHitCollection="EcalBarrelHitsSimpleReco",
+        minModuleEdep=0.0*units.MeV)
+
+simple_cluster = SimpleClustering("simple_cluster",
+        inputHitCollection="EcalBarrelHitsSimpleReco",
+        outputClusters="SimpleClusters",
+        minModuleEdep=1.0*units.MeV,
+        maxDistance=50.0*units.cm)
+
+# Crystal Endcap Ecal
+ce_ecal_daq = dict(
+        dynamicRangeADC=5.*units.GeV,
+        capacityADC=32768,
+        pedestalMean=400,
+        pedestalSigma=3)
+
+ce_ecal_digi = CalHitDigi("ce_ecal_digi",
+        inputHitCollection="EcalEndcapNHits",
+        outputHitCollection="EcalEndcapNHitsDigi",
+        energyResolutions=[0., 0.02, 0.],
+        **ce_ecal_daq)
+
+ce_ecal_reco = CalHitReco("ce_ecal_reco",
+        inputHitCollection="EcalEndcapNHitsDigi",
+        outputHitCollection="EcalEndcapNHitsReco",
+        thresholdFactor=4,          # 4 sigma cut on pedestal sigma
+        readoutClass="EcalEndcapNHits",
+        sectorField="sector",
+        **ce_ecal_daq)
+
+ce_ecal_cl = IslandCluster("ce_ecal_cl",
+        inputHitCollection="EcalEndcapNHitsReco",
+        outputHitCollection="EcalEndcapNClusterHits",
+        splitCluster=False,
+        minClusterHitEdep=1.0*units.MeV,  # discard low energy hits
+        minClusterCenterEdep=30*units.MeV,
+        sectorDist=5.0*units.cm,
+        dimScaledLocalDistXY=[1.8, 1.8]) # dimension scaled dist is good for hybrid sectors with different module size
+
+ce_ecal_clreco = RecoCoG("ce_ecal_clreco",
+        inputHitCollection="EcalEndcapNClusterHits",
+        outputClusterCollection="EcalEndcapNClusters",
+        samplingFraction=0.998,      # this accounts for a small fraction of leakage
+        logWeightBase=4.6)
+
+# Endcap Sampling Ecal
+ci_ecal_daq = dict(
+        dynamicRangeADC=50.*units.MeV,
+        capacityADC=32768,
+        pedestalMean=400,
+        pedestalSigma=10)
+
+ci_ecal_digi = CalHitDigi("ci_ecal_digi",
+        inputHitCollection="EcalEndcapPHits",
+        outputHitCollection="EcalEndcapPHitsDigi",
+        **ci_ecal_daq)
+
+ci_ecal_reco = CalHitReco("ci_ecal_reco",
+        inputHitCollection="EcalEndcapPHitsDigi",
+        outputHitCollection="EcalEndcapPHitsReco",
+        thresholdFactor=5.0,
+        **ci_ecal_daq)
+
+# merge hits in different layer (projection to local x-y plane)
+ci_ecal_merger = CalHitsMerger("ci_ecal_merger",
+        inputHitCollection="EcalEndcapPHitsReco",
+        outputHitCollection="EcalEndcapPHitsRecoXY",
+        fields=["layer", "slice"],
+        fieldRefNumbers=[1, 0],
+        readoutClass="EcalEndcapPHits")
+
+ci_ecal_cl = IslandCluster("ci_ecal_cl",
+        inputHitCollection="EcalEndcapPHitsRecoXY",
+        outputHitCollection="EcalEndcapPClusterHits",
+        splitCluster=False,
+        minClusterCenterEdep=10.*units.MeV,
+        localDistXY=[10*units.mm, 10*units.mm])
+
+ci_ecal_clreco = RecoCoG("ci_ecal_clreco",
+        inputHitCollection="EcalEndcapPClusterHits",
+        outputClusterCollection="EcalEndcapPClusters",
+        logWeightBase=6.2,
+        samplingFraction=ci_ecal_sf)
+
+# Central Barrel Ecal (Imaging Cal.)
+cb_ecal_daq = dict(
+        dynamicRangeADC=3*units.MeV,
+        capacityADC=8192,
+        pedestalMean=400,
+        pedestalSigma=20)   # about 6 keV
+
+cb_ecal_digi = CalHitDigi("cb_ecal_digi",
+        inputHitCollection="EcalBarrelHits",
+        outputHitCollection="EcalBarrelHitsDigi",
+        energyResolutions=[0., 0.02, 0.],   # 2% flat resolution
+        **cb_ecal_daq)
+
+cb_ecal_reco = ImCalPixelReco("cb_ecal_reco",
+        inputHitCollection="EcalBarrelHitsDigi",
+        outputHitCollection="EcalBarrelHitsReco",
+        thresholdFactor=3,  # about 20 keV
+        readoutClass="EcalBarrelHits",  # readout class
+        layerField="layer",             # field to get layer id
+        sectorField="module",           # field to get sector id
+        **cb_ecal_daq)
+
+cb_ecal_cl = ImagingCluster("cb_ecal_cl",
+        inputHitCollection="EcalBarrelHitsReco",
+        outputHitCollection="EcalBarrelClusterHits",
+        localDistXY=[2.*units.mm, 2*units.mm],  # same layer
+        layerDistEtaPhi=[10*units.mrad, 10*units.mrad],     # adjacent layer
+        neighbourLayersRange=2,                 # id diff for adjacent layer
+        sectorDist=3.*units.cm)                       # different sector
+
+cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco",
+        samplingFraction=cb_ecal_sf,
+        inputHitCollection="EcalBarrelClusterHits",
+        outputClusterCollection="EcalBarrelClusters",
+        outputLayerCollection="EcalBarrelLayers")
+
+# Central Barrel Hcal
+cb_hcal_daq = dict(
+         dynamicRangeADC=50.*units.MeV,
+         capacityADC=32768,
+         pedestalMean=400,
+         pedestalSigma=10)
+
+cb_hcal_digi = CalHitDigi("cb_hcal_digi",
+         inputHitCollection="HcalBarrelHits",
+         outputHitCollection="HcalBarrelHitsDigi",
+         **cb_hcal_daq)
+
+cb_hcal_reco = CalHitReco("cb_hcal_reco",
+        inputHitCollection="HcalBarrelHitsDigi",
+        outputHitCollection="HcalBarrelHitsReco",
+        thresholdFactor=5.0,
+        readoutClass="HcalBarrelHits",
+        layerField="layer",
+        sectorField="module",
+        **cb_hcal_daq)
+
+cb_hcal_merger = CalHitsMerger("cb_hcal_merger",
+        inputHitCollection="HcalBarrelHitsReco",
+        outputHitCollection="HcalBarrelHitsRecoXY",
+        readoutClass="HcalBarrelHits",
+        fields=["layer", "slice"],
+        fieldRefNumbers=[1, 0])
+
+cb_hcal_cl = IslandCluster("cb_hcal_cl",
+        inputHitCollection="HcalBarrelHitsRecoXY",
+        outputHitCollection="HcalBarrelClusterHits",
+        splitCluster=False,
+        minClusterCenterEdep=30.*units.MeV,
+        localDistXY=[15.*units.cm, 15.*units.cm])
+
+cb_hcal_clreco = RecoCoG("cb_hcal_clreco",
+        inputHitCollection="HcalBarrelClusterHits",
+        outputClusterCollection="HcalBarrelClusters",
+        logWeightBase=6.2,
+        samplingFraction=cb_hcal_sf)
+
+
+# Hcal Hadron Endcap
+ci_hcal_daq = dict(
+         dynamicRangeADC=50.*units.MeV,
+         capacityADC=32768,
+         pedestalMean=400,
+         pedestalSigma=10)
+ci_hcal_digi = CalHitDigi("ci_hcal_digi",
+         inputHitCollection="HcalHadronEndcapHits",
+         outputHitCollection="HcalHadronEndcapHitsDigi",
+         **ci_hcal_daq)
+
+ci_hcal_reco = CalHitReco("ci_hcal_reco",
+        inputHitCollection="HcalHadronEndcapHitsDigi",
+        outputHitCollection="HcalHadronEndcapHitsReco",
+        thresholdFactor=5.0,
+        **ci_hcal_daq)
+
+ci_hcal_merger = CalHitsMerger("ci_hcal_merger",
+        inputHitCollection="HcalHadronEndcapHitsReco",
+        outputHitCollection="HcalHadronEndcapHitsRecoXY",
+        readoutClass="HcalHadronEndcapHits",
+        fields=["layer", "slice"],
+        fieldRefNumbers=[1, 0])
+
+ci_hcal_cl = IslandCluster("ci_hcal_cl",
+        inputHitCollection="HcalHadronEndcapHitsRecoXY",
+        outputHitCollection="HcalHadronEndcapClusterHits",
+        splitCluster=False,
+        minClusterCenterEdep=30.*units.MeV,
+        localDistXY=[15.*units.cm, 15.*units.cm])
+
+ci_hcal_clreco = RecoCoG("ci_hcal_clreco",
+        inputHitCollection="HcalHadronEndcapClusterHits",
+        outputClusterCollection="HcalHadronEndcapClusters",
+        logWeightBase=6.2,
+        samplingFraction=ci_hcal_sf)
+
+# Hcal Electron Endcap
+ce_hcal_daq = dict(
+        dynamicRangeADC=50.*units.MeV,
+        capacityADC=32768,
+        pedestalMean=400,
+        pedestalSigma=10)
+
+ce_hcal_digi = CalHitDigi("ce_hcal_digi",
+        inputHitCollection="HcalElectronEndcapHits",
+        outputHitCollection="HcalElectronEndcapHitsDigi",
+        **ce_hcal_daq)
+
+ce_hcal_reco = CalHitReco("ce_hcal_reco",
+        inputHitCollection="HcalElectronEndcapHitsDigi",
+        outputHitCollection="HcalElectronEndcapHitsReco",
+        thresholdFactor=5.0,
+        **ce_hcal_daq)
+
+ce_hcal_merger = CalHitsMerger("ce_hcal_merger",
+        inputHitCollection="HcalElectronEndcapHitsReco",
+        outputHitCollection="HcalElectronEndcapHitsRecoXY",
+        readoutClass="HcalElectronEndcapHits",
+        fields=["layer", "slice"],
+        fieldRefNumbers=[1, 0])
+
+ce_hcal_cl = IslandCluster("ce_hcal_cl",
+        inputHitCollection="HcalElectronEndcapHitsRecoXY",
+        outputHitCollection="HcalElectronEndcapClusterHits",
+        splitCluster=False,
+        minClusterCenterEdep=30.*units.MeV,
+        localDistXY=[15.*units.cm, 15.*units.cm])
+
+ce_hcal_clreco = RecoCoG("ce_hcal_clreco",
+        inputHitCollection="HcalElectronEndcapClusterHits",
+        outputClusterCollection="HcalElectronEndcapClusters",
+        logWeightBase=6.2,
+        samplingFraction=ce_hcal_sf)
+
+
+# Tracking
+trk_b_digi = TrackerDigi("trk_b_digi",
+        inputHitCollection="TrackerBarrelHits",
+        outputHitCollection="TrackerBarrelRawHits",
+        timeResolution=8)
+trk_ec_digi = TrackerDigi("trk_ec_digi", 
+        inputHitCollection="TrackerEndcapHits",
+        outputHitCollection="TrackerEndcapRawHits",
+        timeResolution=8)
+
+vtx_b_digi = TrackerDigi("vtx_b_digi", 
+        inputHitCollection="VertexBarrelHits",
+        outputHitCollection="VertexBarrelRawHits",
+        timeResolution=8)
+
+vtx_ec_digi = TrackerDigi("vtx_ec_digi", 
+        inputHitCollection="VertexEndcapHits",
+        outputHitCollection="VertexEndcapRawHits",
+        timeResolution=8)
+
+# Tracker and vertex reconstruction
+trk_b_reco = TrackerHitReconstruction("trk_b_reco",
+        inputHitCollection = trk_b_digi.outputHitCollection,
+        outputHitCollection="TrackerBarrelRecHits")
+
+trk_ec_reco = TrackerHitReconstruction("trk_ec_reco",
+        inputHitCollection = trk_ec_digi.outputHitCollection,
+        outputHitCollection="TrackerEndcapRecHits")
+
+vtx_b_reco = TrackerHitReconstruction("vtx_b_reco",
+        inputHitCollection = vtx_b_digi.outputHitCollection,
+        outputHitCollection="VertexBarrelRecHits")
+
+vtx_ec_reco = TrackerHitReconstruction("vtx_ec_reco",
+        inputHitCollection = vtx_ec_digi.outputHitCollection,
+        outputHitCollection="VertexEndcapRecHits")
+
+# Hit Source linker 
+sourcelinker = TrackerSourceLinker("sourcelinker",
+        inputHitCollection=trk_b_reco.outputHitCollection,
+        outputSourceLinks="BarrelTrackSourceLinks",
+        outputMeasurements="BarrelTrackMeasurements")
+
+#trk_hits_srclnkr = TrackerSourcesLinker("trk_srcslnkr",
+#        ITrackerBarrelHits = vtx_b_reco.outputHitCollection,
+#        ITrackerEndcapHits = vtx_ec_reco.outputHitCollection,
+#        OTrackerBarrelHits = trk_b_reco.outputHitCollection,
+#        OTrackerEndcapHits = trk_ec_reco.outputHitCollection,
+#        outputSourceLinks="TrackerMeasurements")
+
+## Track param init
+truth_trk_init = TrackParamTruthInit("truth_trk_init",
+        inputMCParticles="mcparticles",
+        outputInitialTrackParameters="InitTrackParamsFromTruth")
+
+#clust_trk_init = TrackParamClusterInit("clust_trk_init",
+#        inputClusters="SimpleClusters",
+#        outputInitialTrackParameters="InitTrackParamsFromClusters")
+
+#vtxcluster_trk_init = TrackParamVertexClusterInit("vtxcluster_trk_init",
+#        inputVertexHits="VertexBarrelRecHits",
+#        inputClusters="EcalBarrelClusters",
+#        outputInitialTrackParameters="InitTrackParamsFromVtxClusters",
+#        maxHitRadius=40.0*units.mm)
+
+# Tracking algorithms
+trk_find_alg = TrackFindingAlgorithm("trk_find_alg",
+        inputSourceLinks = sourcelinker.outputSourceLinks,
+        inputMeasurements = sourcelinker.outputMeasurements,
+        inputInitialTrackParameters= "InitTrackParamsFromTruth",
+        outputTrajectories="trajectories")
+
+parts_from_fit = ParticlesFromTrackFit("parts_from_fit",
+        inputTrajectories="trajectories",
+        outputParticles="ReconstructedParticlesInitFromTruth",
+        outputTrackParameters="outputTrackParameters")
+
+#trk_find_alg1 = TrackFindingAlgorithm("trk_find_alg1",
+#        inputSourceLinks = sourcelinker.outputSourceLinks,
+#        inputMeasurements = sourcelinker.outputMeasurements,
+#        inputInitialTrackParameters= "InitTrackParamsFromClusters", 
+#        outputTrajectories="trajectories1")
+
+#parts_from_fit1 = ParticlesFromTrackFit("parts_from_fit1",
+#        inputTrajectories="trajectories1",
+#        outputParticles="ReconstructedParticlesInitFromClusters",
+#        outputTrackParameters="outputTrackParameters1")
+
+#trk_find_alg2 = TrackFindingAlgorithm("trk_find_alg2",
+#        inputSourceLinks = trk_hits_srclnkr.outputSourceLinks,
+#        inputMeasurements = trk_hits_srclnkr.outputMeasurements,
+#        inputInitialTrackParameters= "InitTrackParamsFromVtxClusters", 
+#        outputTrajectories="trajectories2")
+#parts_from_fit2 = ParticlesFromTrackFit("parts_from_fit2",
+#        inputTrajectories="trajectories2",
+#        outputParticles="ReconstructedParticles2",
+#        outputTrackParameters="outputTrackParameters2")
+
+
+pmtdigi = PhotoMultiplierDigi("pmtdigi",
+        inputHitCollection="ForwardRICHHits",
+        outputHitCollection="ForwardRICHHitsDigi",
+        quantumEfficiency=[(a*units.eV, b) for a, b in qe_data])
+
+pmtreco = PhotoMultiplierReco("pmtreco",
+        inputHitCollection="ForwardRICHHitsDigi",
+        outputHitCollection="ForwardRICHHitsReco")
+
+# FIXME
+#richcluster = PhotoRingClusters("richcluster",
+#        inputHitCollection="ForwardRICHHitsReco",
+#        #inputTrackCollection="ReconstructedParticles",
+#        outputClusterCollection="ForwardRICHClusters")
+
+# Output
+podout.outputCommands = ["keep *",
+        "keep *Digi",
+        "keep *Reco*",
+        "keep *ClusterHits",
+        "keep *Clusters",
+        "keep *Layers",
+        "drop mcparticles"
+        ]
+
+ApplicationMgr(
+    TopAlg = [podin, copier, trkcopier,
+              dummy,
+              ecal_digi, ecal_reco, #simple_cluster,
+              ce_ecal_digi, ce_ecal_reco, ce_ecal_cl, ce_ecal_clreco,
+              ci_ecal_digi, ci_ecal_reco, ci_ecal_merger, ci_ecal_cl, ci_ecal_clreco,
+              cb_ecal_digi, cb_ecal_reco, cb_ecal_cl, cb_ecal_clreco,
+              cb_hcal_digi, cb_hcal_reco, cb_hcal_merger, cb_hcal_cl, cb_hcal_clreco,
+              ce_hcal_digi, ce_hcal_reco, ce_hcal_merger, ce_hcal_cl, ce_hcal_clreco,
+              ci_hcal_digi, ci_hcal_reco, ci_hcal_merger, ci_hcal_cl, ci_hcal_clreco,
+              trk_b_digi, trk_ec_digi, vtx_b_digi, vtx_ec_digi,
+              trk_b_reco, trk_ec_reco, vtx_b_reco, vtx_ec_reco, 
+              sourcelinker, #trk_hits_srclnkr,
+              #clust_trk_init,
+              truth_trk_init,
+              #vtxcluster_trk_init, 
+              trk_find_alg, parts_from_fit,
+              #trk_find_alg1, parts_from_fit1,
+              #trk_find_alg2, parts_from_fit2,
+              pmtcopier, pmtdigi, pmtreco, #richcluster,
+              podout
+              ],
+    EvtSel = 'NONE',
+    EvtMax   = n_events,
+    ExtSvc = [podioevent,geo_service],
+    OutputLevel=WARNING
+ )

benchmarks/full/scripts/gen_particles.py

+import os
+from pyHepMC3 import HepMC3 as hm
+import numpy as np
+import argparse
+
+
+PARTICLES = {
+    "pion0": (111, 0.1349766),       # pi0
+    "pion+": (211, 0.13957018),      # pi+
+    "pion-": (-211, 0.13957018),     # pi-
+    "kaon0L": (130, 0.497648),       # K0L
+    "kaon0S": (310, 0.497648),       # K0S
+    "kaon0": (311, 0.497648),        # K0
+    "kaon+": (321, 0.493677),        # K+
+    "kaon-": (-321, 0.493677),       # K-
+    "proton": (2212, 0.938272),      # proton
+    "neutron": (2112, 0.939565),     # neutron
+    "electron": (11, 0.51099895e-3), # electron
+    "positron": (-11, 0.51099895e-3),# positron
+    "photon": (22, 0),               # photon
+}
+
+REGIONS = {
+    "farforward":  (1, 4),
+    "forward":     (3, 50),
+    "barrel":      (45, 135),
+    "backward":    (130, 177),
+    "farbackward": (176, 179),
+}
+
+def gen_event(p, theta, phi, pid, mass):
+    evt = hm.GenEvent(momentum_unit=hm.Units.MomentumUnit.GEV, length_unit=hm.Units.LengthUnit.MM)
+    # final state
+    state = 1
+    e0 = np.sqrt(p*p + mass*mass)
+    px = np.cos(phi)*np.sin(theta)
+    py = np.sin(phi)*np.sin(theta)
+    pz = np.cos(theta)
+
+    # beam
+    pbeam = hm.GenParticle(hm.FourVector(0, 0, 0, 0.938272), 2212, 4)
+    ebeam = hm.GenParticle(hm.FourVector(0, 0, e0, np.sqrt(e0*e0 + 0.511e-3*0.511e-3)), 11, 4)
+
+    # out particle
+    hout = hm.GenParticle(hm.FourVector(px*p, py*p, pz*p, e0), pid, state)
+
+    # vertex
+    vert = hm.GenVertex()
+    vert.add_particle_in(ebeam)
+    vert.add_particle_in(pbeam)
+    vert.add_particle_out(hout)
+    evt.add_vertex(vert)
+    return evt
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser()
+
+    parser.add_argument('output', help='path to the output file')
+    parser.add_argument('-n', type=int, default=1000, dest='nev', help='number of events to generate')
+    parser.add_argument('-s', type=int, default=-1, dest='seed', help='seed for random generator')
+    parser.add_argument('--parray', type=str, default="", dest='parray',
+                        help='an array of momenta in GeV, separated by \",\"')
+    parser.add_argument('--pmin', type=float, default=8.0, dest='pmin', help='minimum momentum in GeV')
+    parser.add_argument('--pmax', type=float, default=100.0, dest='pmax', help='maximum momentum in GeV')
+    parser.add_argument('--angmin', type=float, default=0.0, dest='angmin', help='minimum angle in degree')
+    parser.add_argument('--angmax', type=float, default=20.0, dest='angmax', help='maximum angle in degree')
+    parser.add_argument('--phmin', type=float, default=0.0, dest='phmin', help='minimum angle in degree')
+    parser.add_argument('--phmax', type=float, default=360.0, dest='phmax', help='maximum angle in degree')
+    parser.add_argument('--particles', type=str, default='electron', dest='particles',
+                        help='particle names, support {}'.format(list(PARTICLES.keys())))
+    parser.add_argument('--regions', type=str, default='ignored', dest='regions',
+                        help='detector regions, support {}'.format(list(REGIONS.keys())))
+
+    args = parser.parse_args()
+
+    # random seed (< 0 will get it from enviroment variable 'SEED', or a system random number)
+    if args.seed < 0:
+        args.seed = os.environ.get('SEED', int.from_bytes(os.urandom(4), byteorder='big', signed=False))
+    print("Random seed is {}".format(args.seed))
+    np.random.seed(args.seed)
+
+    output = hm.WriterAscii(args.output);
+    if output.failed():
+        print("Cannot open file \"{}\"".format(args.output))
+        sys.exit(2)
+
+    # build particle info
+    parts = []
+    for pid in args.particles.split(','):
+        pid = pid.strip()
+        if pid not in PARTICLES.keys():
+            print('pid {:d} not found in dictionary, ignored.'.format(pid))
+            continue
+        parts.append(PARTICLES[pid])
+
+    # interpret region
+    (angmin, angmax) = (args.angmin, args.angmax)
+    for region in args.regions.split(','):
+        region = region.strip()
+        if region not in REGIONS.keys():
+            print('region {:d} not found in dictionary, ignored.'.format(pid))
+            continue
+        (angmin, angmax) = REGIONS[region]
+
+    # p values
+    pvals = np.random.uniform(args.pmin, args.pmax, args.nev) if not args.parray else \
+            np.random.choice([float(p.strip()) for p in args.parray.split(',')], args.nev)
+    thvals = np.random.uniform(angmin, angmax, args.nev)/180.*np.pi
+    phivals = np.random.uniform(args.phmin, args.phmax, args.nev)/180.*np.pi
+    partvals = [parts[i] for i in np.random.choice(len(parts), args.nev)]
+
+    count = 0
+    for p, theta, phi, (pid, mass) in zip(pvals, thvals, phivals, partvals):
+        if (count % 1000 == 0):
+            print("Generated {} events".format(count), end='\r')
+        evt = gen_event(p, theta, phi, pid, mass)
+        output.write_event(evt)
+        evt.clear()
+        count += 1
+
+    print("Generated {} events".format(args.nev))
+    output.close()
+

benchmarks/tracking/config.yml

@@ -4,5 +4,4 @@ tracking_central_electrons:
   timeout: 24 hours
   script:
     - bash benchmarks/tracking/central_electrons.sh
-      #allow_failure: true
-
+      #allow_failure: true
\ No newline at end of file