diff --git a/benchmarks/clustering/full_cal_clusters.sh b/benchmarks/clustering/full_cal_clusters.sh index 5cdf14b604b39582ee680481779711df742b624f..0add9a14cb991e08eea8aac0c5bc8b0a5a7f7228 100644 --- a/benchmarks/clustering/full_cal_clusters.sh +++ b/benchmarks/clustering/full_cal_clusters.sh @@ -2,67 +2,137 @@ print_env.sh -## To run the reconstruction, we need the following global variables: -## - JUGGLER_INSTALL_PREFIX: Install prefix for Juggler (simu/recon) -## - JUGGLER_DETECTOR: the detector package we want to use for this benchmark -## - JUGGLER_DETECTOR_VERSION: the detector package we want to use for this benchmark -## - DETECTOR_PATH: full path to the detector definitions -## -## You can ready options/env.sh for more in-depth explanations of the variables -## and how they can be controlled. - -# deprecated -export JUGGLER_DETECTOR_PATH=${DETECTOR_PATH} - -if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then - export JUGGLER_N_EVENTS=100 +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 " --pmin minimum particle momentum (GeV)" + echo " --pmax maximum particle momentum (GeV)" + 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 FULL_CAL_NUMEV="$2" + shift # past argument + shift # past value + ;; + --pmin) + export FULL_CAL_PMIN="$2" + shift + shift + ;; + --pmax) + export FULL_CAL_PMAX="$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 + +export FULL_CAL_COMPACT_PATH=${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml + +if [[ ! -n "${FULL_CAL_NUMEV}" ]] ; then + export FULL_CAL_NUMEV=1000 fi -if [[ ! -n "${CB_EMCAL_SAMP_FRAC}" ]] ; then - export CB_EMCAL_SAMP_FRAC=1.0 +if [[ ! -n "${FULL_CAL_PMIN}" ]] ; then + export FULL_CAL_ENERGY=1.0 fi -export JUGGLER_FILE_NAME_TAG="barrel_clusters" -export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc" +if [[ ! -n "${FULL_CAL_PMAX}" ]] ; then + export FULL_CAL_ENERGY=10.0 +fi + +export FULL_CAL_NAME_TAG="${nametag}" +export FULL_CAL_GEN_FILE="${FULL_CAL_NAME_TAG}.hepmc" + +export FULL_CAL_SIM_FILE="sim_${FULL_CAL_NAME_TAG}.root" +export FULL_CAL_REC_FILE="rec_${FULL_CAL_NAME_TAG}.root" + +echo "FULL_CAL_NUMEV = ${FULL_CAL_NUMEV}" +echo "FULL_CAL_COMPACT_PATH = ${FULL_CAL_COMPACT_PATH}" -export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root" -export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root" +":" << END +# Generate the input events +python benchmarks/imaging_ecal/scripts/gen_particles.py ${FULL_CAL_GEN_FILE} -n ${FULL_CAL_NUMEV}\ + --angmin 60 --angmax 120 --phmin 0 --phmax 360 --pmin ${FULL_CAL_PMIN} --pmax ${FULL_CAL_PMAX}\ + --particles="${particle}" -echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}" -echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}" +if [[ "$?" -ne "0" ]] ; then + echo "ERROR running script: generating input events" + exit 1 +fi -root -b -q "benchmarks/clustering/scripts/gen_central_pions.cxx(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")" +ls -lh ${FULL_CAL_GEN_FILE} -### run geant4 simulations -# part.minimalKineticEnergy sets Ek limit for saving particles in mcparticles branch -# 1*TeV means only saving particles from the event generator +# Run geant4 simulations npsim --runType batch \ - --part.minimalKineticEnergy "1*TeV" \ -v WARNING \ - --numberOfEvents ${JUGGLER_N_EVENTS} \ - --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \ - --inputFiles ${JUGGLER_FILE_NAME_TAG}.hepmc \ - --outputFile ${JUGGLER_SIM_FILE} + --part.minimalKineticEnergy "1*TeV" \ + --numberOfEvents ${FULL_CAL_NUMEV} \ + --compactFile ${FULL_CAL_COMPACT_PATH} \ + --inputFiles ${FULL_CAL_GEN_FILE} \ + --outputFile ${FULL_CAL_SIM_FILE} + if [[ "$?" -ne "0" ]] ; then echo "ERROR running npdet" exit 1 fi +END -# Need to figure out how to pass file name to juggler from the commandline -gaudirun.py benchmarks/clustering/options/fullcalo_clustering.py +rootls -t "${FULL_CAL_SIM_FILE}" + +# Directory for results +mkdir -p results + +FULL_CAL_OPTION_DIR=benchmarks/clustering/options +FULL_CAL_SCRIPT_DIR=benchmarks/clustering/scripts + +# Run Juggler +xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \ + gaudirun.py ${FULL_CAL_OPTION_DIR}/full_cal_clustering.py if [[ "$?" -ne "0" ]] ; then echo "ERROR running juggler" exit 1 fi -mkdir -p results/clustering +# TODO analysis scripts -#root -b -q "benchmarks/clustering/scripts/barrel_clusters.cxx(\"${JUGGLER_REC_FILE}\")" +root_filesize=$(stat --format=%s "${FULL_CAL_REC_FILE}") +if [[ "${FULL_CAL_NUMEV}" -lt "500" ]] ; then + # file must be less than 10 MB to upload + if [[ "${root_filesize}" -lt "10000000" ]] ; then + cp ${FULL_CAL_REC_FILE} results/. + fi +fi -#root_filesize=$(stat --format=%s "${JUGGLER_DETECTOR}/${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/clustering/. -# fi -#fi diff --git a/benchmarks/clustering/options/full_cal_clustering.py b/benchmarks/clustering/options/full_cal_clustering.py index 39bb28eeb213122e307707205631a1a028149c15..49d82f2860bdf6e9528f828c44a5918549b0be0d 100644 --- a/benchmarks/clustering/options/full_cal_clustering.py +++ b/benchmarks/clustering/options/full_cal_clustering.py @@ -4,7 +4,7 @@ import ROOT from GaudiKernel.DataObjectHandleBase import DataObjectHandleBase from Configurables import ApplicationMgr, EICDataSvc, PodioInput, PodioOutput, GeoSvc -from GaudiKernel.SystemOfUnits import MeV, GeV, cm +from GaudiKernel.SystemOfUnits import MeV, GeV, mm, cm, mrad detector_name = str(os.environ.get("JUGGLER_DETECTOR", "athena")) detector_path = str(os.environ.get("DETECTOR_PATH", ".")) @@ -15,14 +15,14 @@ cb_ecal_sf = float(os.environ.get("CB_EMCAL_SAMP_FRAC", 1.0)) cb_hcal_sf = float(os.environ.get("CB_HCAL_SAMP_FRAC", 1.0)) # 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"]) +input_sims = [f.strip() for f in str.split(os.environ["FULL_CAL_SIM_FILE"], ",") if f.strip()] +output_rec = str(os.environ["FULL_CAL_REC_FILE"]) +n_events = int(os.environ["FULL_CAL_NUMEV"]) # geometry service -geo_service = GeoSvc("GeoSvc", detectors=["{}.xml".format(detector_name)]) +geo_service = GeoSvc("GeoSvc", detectors=["{}.xml".format(compact_path)]) # data service -podioevent = EICDataSvc("EventDataSvc", inputs=[input_sims]) +podioevent = EICDataSvc("EventDataSvc", inputs=input_sims) # juggler components @@ -48,11 +48,11 @@ sim_coll = [ # input and output podin = PodioInput("PodioReader", collections=sim_coll) -podout = PodioOutput("out", filename=output_rec_file) - +podout = PodioOutput("out", filename=output_rec) # copier needed to get around input --> output copy bug. So truth (mcparticles) can be saved in output file copier = MCCopier("MCCopier", inputCollection="mcparticles", outputCollection="mcparticles2") + # Crystal Endcap Ecal ce_ecal_digi = CalHitDigi("ce_ecal_digi", inputHitCollection="CrystalEcalHits", @@ -77,8 +77,9 @@ ce_ecal_cl = IslandCluster("ce_ecal_cl", # OutputLevel=DEBUG, inputHitCollection="CrystalEcalHitsReco", outputClusterCollection="CrystalEcalClusters", + splitHitCollection="CrystalEcalHitsSplit", splitCluster=False, - minClusterCenterEdep=30*units.MeV, + minClusterCenterEdep=30*MeV, groupRanges=[2.2*cm, 2.2*cm]) ce_ecal_clreco = RecoCoG("ce_ecal_clreco", @@ -86,7 +87,7 @@ ce_ecal_clreco = RecoCoG("ce_ecal_clreco", logWeightBase=4.6) -# central barrel Ecal +# Central Barrel Ecal (Imaging Cal.) cb_ecal_digi = CalHitDigi("cb_ecal_digi", inputHitCollection="EcalBarrelHits", outputHitCollection="EcalBarrelHitsDigi", @@ -95,7 +96,7 @@ cb_ecal_digi = CalHitDigi("cb_ecal_digi", capacityADC=8192, pedestalMean=400, pedestalSigma=20) # about 6 keV -cb_ecal_reco = ImaCalPixelReco("cb_ecal_reco", +cb_ecal_reco = ImCalPixelReco("cb_ecal_reco", inputHitCollection="EcalBarrelHitsDigi", outputHitCollection="EcalBarrelHitsReco", dynamicRangeADC=3.*MeV, @@ -109,15 +110,16 @@ cb_ecal_reco = ImaCalPixelReco("cb_ecal_reco", cb_ecal_cl = ImagingCluster("cb_ecal_cl", inputHitCollection="EcalBarrelHitsReco", outputClusterCollection="EcalBarrelClusters", - localRanges=[2.*units.mm, 2*units.mm], # same layer - adjLayerRanges=[10*units.mrad, 10*units.mrad], # adjacent layer - adjLayerDiff=2, # id diff for adjacent layer - adjSectorDist=3.*units.cm) # different sector + localRanges=[2.*mm, 2*mm], # same layer + adjLayerRanges=[10*mrad, 10*mrad], # adjacent layer + adjLayerDiff=2, # id diff for adjacent layer + adjSectorDist=3.*cm) # different sector cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco", inputClusterCollection="EcalBarrelClusters", outputLayerCollection="EcalBarrelLayers") -# central barrel hcal + +# Central Barrel Hcal cb_hcal_digi = CalHitDigi("cb_hcal_digi", inputHitCollection="HcalBarrelHits", outputHitCollection="HcalBarrelHitsDigi", @@ -145,25 +147,24 @@ cb_hcal_merger = CalHitsMerger("cb_hcal_merger", cb_hcal_cl = IslandCluster("cb_hcal_cl", inputHitCollection="HcalBarrelHitsRecoXY", outputClusterCollection="HcalBarrelClusters", + splitHitCollection="HcalBarrelHitsSplit", splitCluster=False, minClusterCenterEdep=30.*MeV, - groupRanges=[15.*cm, 15.*cm], - OutputLevel=DEBUG) + groupRanges=[15.*cm, 15.*cm]) cb_hcal_clreco = RecoCoG("cb_hcal_clreco", clusterCollection="HcalBarrelClusters", logWeightBase=6.2, - samplingFraction=cb_hcal_sf, - OutputLevel=DEBUG) + samplingFraction=cb_hcal_sf) -out.outputCommands = ["keep *"] +podout.outputCommands = ["keep *"] ApplicationMgr( - TopAlg = [podioinput, copier, + TopAlg = [podin, copier, ce_ecal_digi, ce_ecal_reco, ce_ecal_cl, ce_ecal_clreco, cb_ecal_digi, cb_ecal_reco, cb_ecal_cl, cb_ecal_clreco, - cb_hcal_digi, cb_hcal_reco, cb_hcal_cl, cb_hcal_clreco, - out], + cb_hcal_digi, cb_hcal_reco, cb_hcal_merger, cb_hcal_cl, cb_hcal_clreco, + podout], EvtSel = 'NONE', EvtMax = n_events, ExtSvc = [podioevent],