diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index eeb2c0e7671ec69732bf55e4e9860c3a1f02fd18..251e3cd92cdd1bfa76c2657ced11400ced72f58b 100644 --- a/.gitlab-ci.yml +++ b/.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 diff --git a/benchmarks/full/config.yml b/benchmarks/full/config.yml new file mode 100644 index 0000000000000000000000000000000000000000..3ad6d6ef0bef04368950c78c4172db4cc69b2638 --- /dev/null +++ b/benchmarks/full/config.yml @@ -0,0 +1,20 @@ +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 + diff --git a/benchmarks/full/full_reconstruction.sh b/benchmarks/full/full_reconstruction.sh new file mode 100644 index 0000000000000000000000000000000000000000..c0b57c645b22a63f21a342642c661bbafa6703ee --- /dev/null +++ b/benchmarks/full/full_reconstruction.sh @@ -0,0 +1,184 @@ +#!/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 + diff --git a/benchmarks/full/options/full_reconstruction.py b/benchmarks/full/options/full_reconstruction.py new file mode 100644 index 0000000000000000000000000000000000000000..a40f0eddad05d6c7449bd26c62c45eacf3a3c8d4 --- /dev/null +++ b/benchmarks/full/options/full_reconstruction.py @@ -0,0 +1,513 @@ +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 + ) diff --git a/benchmarks/full/scripts/gen_particles.py b/benchmarks/full/scripts/gen_particles.py new file mode 100644 index 0000000000000000000000000000000000000000..0609714c5bd10b38230e82c6de12be6768be7b08 --- /dev/null +++ b/benchmarks/full/scripts/gen_particles.py @@ -0,0 +1,124 @@ +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() + diff --git a/benchmarks/tracking/config.yml b/benchmarks/tracking/config.yml index 52c3df4d90fc4b082b4890fbcfaf09f2fac26ea7..6e5c5d2f9694c97f7a5773c501424a3b68f387c8 100644 --- a/benchmarks/tracking/config.yml +++ b/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