add benchmark for pion0

benchmarks/imaging_ecal/config.yml

@@ -19,3 +19,10 @@ imaging_ecal_pions:
   script:
     - bash benchmarks/imaging_ecal/run_emcal_barrel.sh -t emcal_barrel_pions -p "pion-" -n 100
 
+imaging_ecal_pion0:
+  extends: .rec_benchmark
+  timeout: 48 hours
+  stage: run
+  script:
+    - bash benchmarks/imaging_ecal/run_imcal_pion0.sh -t imcal_barrel_pion0 -p "pion0" -n 100
+

benchmarks/imaging_ecal/options/imaging_2dcluster.py

+import os
+import ROOT
+from Gaudi.Configuration import *
+from GaudiKernel.SystemOfUnits import mm, MeV, rad, ns
+from GaudiKernel.DataObjectHandleBase import DataObjectHandleBase
+from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
+
+from Configurables import PodioInput
+from Configurables import Jug__Digi__CalorimeterHitDigi as CalHitDigi
+from Configurables import Jug__Reco__CalorimeterHitReco as CalHitReco
+from Configurables import Jug__Reco__CalorimeterHitsEtaPhiProjector as CalHitsProj
+from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
+from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
+
+
+# input arguments through environment variables
+kwargs = dict()
+kwargs['sf'] = float(os.environ.get('CB_EMCAL_SCFI, SAMP_FRAC', '0.0134'))
+kwargs['input'] = os.environ.get('CB_EMCAL_SIM_FILE', '../topside/barrel_pion0_5GeV.root')
+kwargs['output'] = os.environ.get('CB_EMCAL_REC_FILE', 'barrel_pion0_5GeV_cluster.root')
+kwargs['compact'] = os.environ.get('CB_EMCAL_COMPACT_PATH', '../topside/test.xml')
+kwargs['nev'] = int(os.environ.get('CB_EMCAL_NUMEV', 100))
+
+if kwargs['nev'] < 1:
+    f = ROOT.TFile(kwargs['input'])
+    kwargs['nev'] = f.events.GetEntries()
+
+print(kwargs)
+
+# get sampling fraction from system environment variable, 1.0 by default
+sf = float(os.environ.get('CB_EMCAL_SAMP_FRAC', '1.0'))
+geo_service = GeoSvc('GeoSvc', detectors=kwargs['compact'].split(','), OutputLevel=INFO)
+podioevent = EICDataSvc('EventDataSvc', inputs=kwargs['input'].split(','), OutputLevel=DEBUG)
+
+podin = PodioInput('PodioReader', collections=['mcparticles', 'EcalBarrelHits'], OutputLevel=DEBUG)
+podout = PodioOutput('out', filename=kwargs['output'])
+
+# use the same daq_setting for digi/reco pair
+imcal_barrel_daq = dict(
+        dynamicRangeADC=3.*MeV,
+        capacityADC=32768,
+        pedestalMean=400,
+        pedestalSigma=50)
+
+imcal_barrel_digi = CalHitDigi('imcal_barrel_digi',
+        inputHitCollection='EcalBarrelHits',
+        outputHitCollection='EcalBarrelHitsDigi',
+        **imcal_barrel_daq)
+
+imcal_barrel_reco = CalHitReco('imcal_barrel_reco',
+        inputHitCollection=imcal_barrel_digi.outputHitCollection,
+        outputHitCollection='EcalBarrelHitsReco',
+        thresholdFactor=5.0,
+        readoutClass='EcalBarrelHits',
+        localDetFields=['system', 'module'], # use local coordinates in each module (stave)
+        **imcal_barrel_daq)
+
+# merge hits to eta, phi bins
+imcal_barrel_merger = CalHitsProj('imcal_barrel_merger',
+        # OutputLevel=DEBUG,
+        inputHitCollection=imcal_barrel_reco.outputHitCollection,
+        outputHitCollection='EcalBarrelHitsMerge',
+        gridSizes=[0.004, 0.004*rad])
+
+imcal_barrel_cl = IslandCluster('imcal_barrel_cl',
+        # OutputLevel=DEBUG,
+        inputHitCollection=imcal_barrel_merger.outputHitCollection,
+        outputHitCollection='EcalBarrelClusterHits',
+        minClusterCenterEdep=1.*MeV,
+        minClusterHitEdep=0.1*MeV,
+        splitCluster=True,
+        globalDistEtaPhi=[0.005, 0.005*rad])
+
+imcal_barrel_clreco = RecoCoG('imcal_barrel_clreco',
+        # OutputLevel=DEBUG,
+        inputHitCollection=imcal_barrel_cl.outputHitCollection,
+        outputClusterCollection='EcalBarrelClusters',
+        logWeightBase=6.2,
+        samplingFraction=kwargs['sf'])
+
+
+podout.outputCommands = ['keep *']
+
+ApplicationMgr(
+    TopAlg=[podin,
+        imcal_barrel_digi, imcal_barrel_reco,
+        imcal_barrel_merger, imcal_barrel_cl, imcal_barrel_clreco,
+        podout],
+    EvtSel='NONE',
+    EvtMax=kwargs['nev'],
+    ExtSvc=[podioevent],
+    OutputLevel=DEBUG
+)
+
+

benchmarks/imaging_ecal/run_imcal_pion0.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_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 --parray ${CB_EMCAL_ENERGY} --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 "1*TeV" \
+      --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
+
+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}/imaging_2dcluster.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 "EcalBarrelClusters"
+
+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
+