Skip to content
Snippets Groups Projects
Commit 09566ddd authored by Chao Peng's avatar Chao Peng
Browse files

update scripts

parent d8496373
No related branches found
No related tags found
1 merge request!107modify clustering to use the latest components
This commit is part of merge request !107. Comments created here will be created in the context of that merge request.
...@@ -2,67 +2,137 @@ ...@@ -2,67 +2,137 @@
print_env.sh print_env.sh
## To run the reconstruction, we need the following global variables: function print_the_help {
## - JUGGLER_INSTALL_PREFIX: Install prefix for Juggler (simu/recon) echo "USAGE: ${0} -n <nevents> -t <nametag> -p <particle> "
## - JUGGLER_DETECTOR: the detector package we want to use for this benchmark echo " OPTIONS: "
## - JUGGLER_DETECTOR_VERSION: the detector package we want to use for this benchmark echo " -n,--nevents Number of events"
## - DETECTOR_PATH: full path to the detector definitions echo " -t,--nametag name tag"
## echo " -p,--particle particle type"
## You can ready options/env.sh for more in-depth explanations of the variables echo " --pmin minimum particle momentum (GeV)"
## and how they can be controlled. echo " --pmax maximum particle momentum (GeV)"
echo " allowed types: pion0, pion+, pion-, kaon0, kaon+, kaon-, proton, neutron, electron, positron, photon"
# deprecated exit
export JUGGLER_DETECTOR_PATH=${DETECTOR_PATH} }
if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then POSITIONAL=()
export JUGGLER_N_EVENTS=100 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 fi
if [[ ! -n "${CB_EMCAL_SAMP_FRAC}" ]] ; then if [[ ! -n "${FULL_CAL_PMIN}" ]] ; then
export CB_EMCAL_SAMP_FRAC=1.0 export FULL_CAL_ENERGY=1.0
fi fi
export JUGGLER_FILE_NAME_TAG="barrel_clusters" if [[ ! -n "${FULL_CAL_PMAX}" ]] ; then
export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc" 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" ":" << END
export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root" # 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}" if [[ "$?" -ne "0" ]] ; then
echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}" 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 # Run geant4 simulations
# part.minimalKineticEnergy sets Ek limit for saving particles in mcparticles branch
# 1*TeV means only saving particles from the event generator
npsim --runType batch \ npsim --runType batch \
--part.minimalKineticEnergy "1*TeV" \
-v WARNING \ -v WARNING \
--numberOfEvents ${JUGGLER_N_EVENTS} \ --part.minimalKineticEnergy "1*TeV" \
--compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \ --numberOfEvents ${FULL_CAL_NUMEV} \
--inputFiles ${JUGGLER_FILE_NAME_TAG}.hepmc \ --compactFile ${FULL_CAL_COMPACT_PATH} \
--outputFile ${JUGGLER_SIM_FILE} --inputFiles ${FULL_CAL_GEN_FILE} \
--outputFile ${FULL_CAL_SIM_FILE}
if [[ "$?" -ne "0" ]] ; then if [[ "$?" -ne "0" ]] ; then
echo "ERROR running npdet" echo "ERROR running npdet"
exit 1 exit 1
fi fi
END
# Need to figure out how to pass file name to juggler from the commandline rootls -t "${FULL_CAL_SIM_FILE}"
gaudirun.py benchmarks/clustering/options/fullcalo_clustering.py
# 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 if [[ "$?" -ne "0" ]] ; then
echo "ERROR running juggler" echo "ERROR running juggler"
exit 1 exit 1
fi 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
...@@ -4,7 +4,7 @@ import ROOT ...@@ -4,7 +4,7 @@ import ROOT
from GaudiKernel.DataObjectHandleBase import DataObjectHandleBase from GaudiKernel.DataObjectHandleBase import DataObjectHandleBase
from Configurables import ApplicationMgr, EICDataSvc, PodioInput, PodioOutput, GeoSvc 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_name = str(os.environ.get("JUGGLER_DETECTOR", "athena"))
detector_path = str(os.environ.get("DETECTOR_PATH", ".")) 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)) ...@@ -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)) cb_hcal_sf = float(os.environ.get("CB_HCAL_SAMP_FRAC", 1.0))
# input and output # input and output
input_sims = [f.strip() for f in str.split(os.environ["JUGGLER_SIM_FILE"], ",") if f.strip()] input_sims = [f.strip() for f in str.split(os.environ["FULL_CAL_SIM_FILE"], ",") if f.strip()]
output_rec = str(os.environ["JUGGLER_REC_FILE"]) output_rec = str(os.environ["FULL_CAL_REC_FILE"])
n_events = int(os.environ["JUGGLER_N_EVENTS"]) n_events = int(os.environ["FULL_CAL_NUMEV"])
# geometry service # geometry service
geo_service = GeoSvc("GeoSvc", detectors=["{}.xml".format(detector_name)]) geo_service = GeoSvc("GeoSvc", detectors=["{}.xml".format(compact_path)])
# data service # data service
podioevent = EICDataSvc("EventDataSvc", inputs=[input_sims]) podioevent = EICDataSvc("EventDataSvc", inputs=input_sims)
# juggler components # juggler components
...@@ -48,11 +48,11 @@ sim_coll = [ ...@@ -48,11 +48,11 @@ sim_coll = [
# input and output # input and output
podin = PodioInput("PodioReader", collections=sim_coll) 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 needed to get around input --> output copy bug. So truth (mcparticles) can be saved in output file
copier = MCCopier("MCCopier", inputCollection="mcparticles", outputCollection="mcparticles2") copier = MCCopier("MCCopier", inputCollection="mcparticles", outputCollection="mcparticles2")
# Crystal Endcap Ecal # Crystal Endcap Ecal
ce_ecal_digi = CalHitDigi("ce_ecal_digi", ce_ecal_digi = CalHitDigi("ce_ecal_digi",
inputHitCollection="CrystalEcalHits", inputHitCollection="CrystalEcalHits",
...@@ -77,8 +77,9 @@ ce_ecal_cl = IslandCluster("ce_ecal_cl", ...@@ -77,8 +77,9 @@ ce_ecal_cl = IslandCluster("ce_ecal_cl",
# OutputLevel=DEBUG, # OutputLevel=DEBUG,
inputHitCollection="CrystalEcalHitsReco", inputHitCollection="CrystalEcalHitsReco",
outputClusterCollection="CrystalEcalClusters", outputClusterCollection="CrystalEcalClusters",
splitHitCollection="CrystalEcalHitsSplit",
splitCluster=False, splitCluster=False,
minClusterCenterEdep=30*units.MeV, minClusterCenterEdep=30*MeV,
groupRanges=[2.2*cm, 2.2*cm]) groupRanges=[2.2*cm, 2.2*cm])
ce_ecal_clreco = RecoCoG("ce_ecal_clreco", ce_ecal_clreco = RecoCoG("ce_ecal_clreco",
...@@ -86,7 +87,7 @@ ce_ecal_clreco = RecoCoG("ce_ecal_clreco", ...@@ -86,7 +87,7 @@ ce_ecal_clreco = RecoCoG("ce_ecal_clreco",
logWeightBase=4.6) logWeightBase=4.6)
# central barrel Ecal # Central Barrel Ecal (Imaging Cal.)
cb_ecal_digi = CalHitDigi("cb_ecal_digi", cb_ecal_digi = CalHitDigi("cb_ecal_digi",
inputHitCollection="EcalBarrelHits", inputHitCollection="EcalBarrelHits",
outputHitCollection="EcalBarrelHitsDigi", outputHitCollection="EcalBarrelHitsDigi",
...@@ -95,7 +96,7 @@ cb_ecal_digi = CalHitDigi("cb_ecal_digi", ...@@ -95,7 +96,7 @@ cb_ecal_digi = CalHitDigi("cb_ecal_digi",
capacityADC=8192, capacityADC=8192,
pedestalMean=400, pedestalMean=400,
pedestalSigma=20) # about 6 keV pedestalSigma=20) # about 6 keV
cb_ecal_reco = ImaCalPixelReco("cb_ecal_reco", cb_ecal_reco = ImCalPixelReco("cb_ecal_reco",
inputHitCollection="EcalBarrelHitsDigi", inputHitCollection="EcalBarrelHitsDigi",
outputHitCollection="EcalBarrelHitsReco", outputHitCollection="EcalBarrelHitsReco",
dynamicRangeADC=3.*MeV, dynamicRangeADC=3.*MeV,
...@@ -109,15 +110,16 @@ cb_ecal_reco = ImaCalPixelReco("cb_ecal_reco", ...@@ -109,15 +110,16 @@ cb_ecal_reco = ImaCalPixelReco("cb_ecal_reco",
cb_ecal_cl = ImagingCluster("cb_ecal_cl", cb_ecal_cl = ImagingCluster("cb_ecal_cl",
inputHitCollection="EcalBarrelHitsReco", inputHitCollection="EcalBarrelHitsReco",
outputClusterCollection="EcalBarrelClusters", outputClusterCollection="EcalBarrelClusters",
localRanges=[2.*units.mm, 2*units.mm], # same layer localRanges=[2.*mm, 2*mm], # same layer
adjLayerRanges=[10*units.mrad, 10*units.mrad], # adjacent layer adjLayerRanges=[10*mrad, 10*mrad], # adjacent layer
adjLayerDiff=2, # id diff for adjacent layer adjLayerDiff=2, # id diff for adjacent layer
adjSectorDist=3.*units.cm) # different sector adjSectorDist=3.*cm) # different sector
cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco", cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco",
inputClusterCollection="EcalBarrelClusters", inputClusterCollection="EcalBarrelClusters",
outputLayerCollection="EcalBarrelLayers") outputLayerCollection="EcalBarrelLayers")
# central barrel hcal
# Central Barrel Hcal
cb_hcal_digi = CalHitDigi("cb_hcal_digi", cb_hcal_digi = CalHitDigi("cb_hcal_digi",
inputHitCollection="HcalBarrelHits", inputHitCollection="HcalBarrelHits",
outputHitCollection="HcalBarrelHitsDigi", outputHitCollection="HcalBarrelHitsDigi",
...@@ -145,25 +147,24 @@ cb_hcal_merger = CalHitsMerger("cb_hcal_merger", ...@@ -145,25 +147,24 @@ cb_hcal_merger = CalHitsMerger("cb_hcal_merger",
cb_hcal_cl = IslandCluster("cb_hcal_cl", cb_hcal_cl = IslandCluster("cb_hcal_cl",
inputHitCollection="HcalBarrelHitsRecoXY", inputHitCollection="HcalBarrelHitsRecoXY",
outputClusterCollection="HcalBarrelClusters", outputClusterCollection="HcalBarrelClusters",
splitHitCollection="HcalBarrelHitsSplit",
splitCluster=False, splitCluster=False,
minClusterCenterEdep=30.*MeV, minClusterCenterEdep=30.*MeV,
groupRanges=[15.*cm, 15.*cm], groupRanges=[15.*cm, 15.*cm])
OutputLevel=DEBUG)
cb_hcal_clreco = RecoCoG("cb_hcal_clreco", cb_hcal_clreco = RecoCoG("cb_hcal_clreco",
clusterCollection="HcalBarrelClusters", clusterCollection="HcalBarrelClusters",
logWeightBase=6.2, logWeightBase=6.2,
samplingFraction=cb_hcal_sf, samplingFraction=cb_hcal_sf)
OutputLevel=DEBUG)
out.outputCommands = ["keep *"] podout.outputCommands = ["keep *"]
ApplicationMgr( ApplicationMgr(
TopAlg = [podioinput, copier, TopAlg = [podin, copier,
ce_ecal_digi, ce_ecal_reco, ce_ecal_cl, ce_ecal_clreco, 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_ecal_digi, cb_ecal_reco, cb_ecal_cl, cb_ecal_clreco,
cb_hcal_digi, cb_hcal_reco, cb_hcal_cl, cb_hcal_clreco, cb_hcal_digi, cb_hcal_reco, cb_hcal_merger, cb_hcal_cl, cb_hcal_clreco,
out], podout],
EvtSel = 'NONE', EvtSel = 'NONE',
EvtMax = n_events, EvtMax = n_events,
ExtSvc = [podioevent], ExtSvc = [podioevent],
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment