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

Fix benchmarks for Imaging Cal.

parent 31e0c968
No related branches found
No related tags found
1 merge request!106Fix benchmarks for Imaging Cal.
Hide whitespace changes
Inline Side-by-side
benchmarks/imaging_ecal/options/imaging_topocluster.py
+ 5
− 4
View file @ a9d56427
...@@ -31,7 +31,6 @@ print(kwargs) ...@@ -31,7 +31,6 @@ print(kwargs)
# get sampling fraction from system environment variable, 1.0 by default # get sampling fraction from system environment variable, 1.0 by default
sf = float(os.environ.get('CB_EMCAL_SAMP_FRAC', '1.0')) sf = float(os.environ.get('CB_EMCAL_SAMP_FRAC', '1.0'))
geo_service = GeoSvc("GeoSvc", detectors=kwargs['compact'].split(',')) geo_service = GeoSvc("GeoSvc", detectors=kwargs['compact'].split(','))
podioevent = EICDataSvc("EventDataSvc", inputs=kwargs['input'].split(','), OutputLevel=DEBUG) podioevent = EICDataSvc("EventDataSvc", inputs=kwargs['input'].split(','), OutputLevel=DEBUG)
out = PodioOutput("out", filename=kwargs['output']) out = PodioOutput("out", filename=kwargs['output'])
...@@ -64,13 +63,15 @@ imcalreco = ImagingPixelReco("imcal_reco", ...@@ -64,13 +63,15 @@ imcalreco = ImagingPixelReco("imcal_reco",
readoutClass="EcalBarrelHits", readoutClass="EcalBarrelHits",
layerField="layer", layerField="layer",
sectorField="module") sectorField="module")
imcalcluster = ImagingTopoCluster(inputHitCollection="RecoEcalBarrelHits", imcalcluster = ImagingTopoCluster("imcal_cluster",
inputHitCollection="RecoEcalBarrelHits",
outputClusterCollection="EcalBarrelClusters", outputClusterCollection="EcalBarrelClusters",
localRanges=[2.*units.mm, 2*units.mm], localRanges=[2.*units.mm, 2*units.mm],
adjLayerRanges=[10*units.mrad, 10*units.mrad], adjLayerRanges=[10*units.mrad, 10*units.mrad],
adjLayerDiff=2, adjLayerDiff=2,
adjSectorDist=3.*units.cm) adjSectorDist=3.*units.cm)
clusterreco = ImagingClusterReco(inputClusterCollection="EcalBarrelClusters", clusterreco = ImagingClusterReco("imcal_clreco",
inputClusterCollection="EcalBarrelClusters",
outputLayerCollection="EcalBarrelClustersLayers", outputLayerCollection="EcalBarrelClustersLayers",
samplingFraction=sf, samplingFraction=sf,
OutputLevel=DEBUG) OutputLevel=DEBUG)
...@@ -78,7 +79,7 @@ clusterreco = ImagingClusterReco(inputClusterCollection="EcalBarrelClusters", ...@@ -78,7 +79,7 @@ clusterreco = ImagingClusterReco(inputClusterCollection="EcalBarrelClusters",
out.outputCommands = ["keep *"] out.outputCommands = ["keep *"]
ApplicationMgr( ApplicationMgr(
TopAlg=[podioinput, copier, calcopier, imcaldigi, imcalreco, imcalcluster, clusterreco, out], TopAlg=[podioinput, copier, calcopier, imcaldigi, imcalreco, imcalcluster, out],
EvtSel='NONE', EvtSel='NONE',
EvtMax=kwargs['nev'], EvtMax=kwargs['nev'],
ExtSvc=[podioevent], ExtSvc=[podioevent],
......
benchmarks/imaging_ecal/run_emcal_barrel.sh
+ 14
− 10
View file @ a9d56427
...@@ -5,11 +5,11 @@ print_env.sh ...@@ -5,11 +5,11 @@ print_env.sh
function print_the_help { function print_the_help {
echo "USAGE: ${0} -n <nevents> -t <nametag> -p <particle> " echo "USAGE: ${0} -n <nevents> -t <nametag> -p <particle> "
echo " OPTIONS: " echo " OPTIONS: "
echo " -n,--nevents Number of eventsj" echo " -n,--nevents Number of events"
echo " -t,--nametag name tag" echo " -t,--nametag name tag"
echo " -p,--particle particle type" echo " -p,--particle particle type"
echo " allowed types: pion0, pion+, pion-, kaon0, kaon+, kaon-, proton, neutron, electron, positron, photon" echo " allowed types: pion0, pion+, pion-, kaon0, kaon+, kaon-, proton, neutron, electron, positron, photon"
exit exit
} }
POSITIONAL=() POSITIONAL=()
...@@ -33,9 +33,9 @@ do ...@@ -33,9 +33,9 @@ do
shift # past value shift # past value
;; ;;
-n|--nevents) -n|--nevents)
export CB_EMCAL_NUMEV=${OPTARG} export CB_EMCAL_NUMEV="$2"
shift # past argument shift # past argument
#shift # past value shift # past value
;; ;;
*) # unknown option *) # unknown option
#POSITIONAL+=("$1") # save it in an array for later #POSITIONAL+=("$1") # save it in an array for later
...@@ -54,7 +54,7 @@ if [[ ! -n "${CB_EMCAL_NUMEV}" ]] ; then ...@@ -54,7 +54,7 @@ if [[ ! -n "${CB_EMCAL_NUMEV}" ]] ; then
fi fi
if [[ ! -n "${CB_EMCAL_IEV}" ]] ; then if [[ ! -n "${CB_EMCAL_IEV}" ]] ; then
export CB_EMCAL_IEV="0, 1, 2, 3, 4" export CB_EMCAL_IEV="0, 1"
fi fi
if [[ ! -n "${CB_EMCAL_ENERGY}" ]] ; then if [[ ! -n "${CB_EMCAL_ENERGY}" ]] ; then
...@@ -87,6 +87,7 @@ ls -lh ${CB_EMCAL_GEN_FILE} ...@@ -87,6 +87,7 @@ ls -lh ${CB_EMCAL_GEN_FILE}
# Run geant4 simulations # Run geant4 simulations
npsim --runType batch \ npsim --runType batch \
-v WARNING \ -v WARNING \
--part.minimalKineticEnergy "1*TeV" \
--numberOfEvents ${CB_EMCAL_NUMEV} \ --numberOfEvents ${CB_EMCAL_NUMEV} \
--compactFile ${CB_EMCAL_COMPACT_PATH} \ --compactFile ${CB_EMCAL_COMPACT_PATH} \
--inputFiles ${CB_EMCAL_GEN_FILE} \ --inputFiles ${CB_EMCAL_GEN_FILE} \
...@@ -108,7 +109,8 @@ CB_EMCAL_OPTION_DIR=benchmarks/imaging_ecal/options ...@@ -108,7 +109,8 @@ CB_EMCAL_OPTION_DIR=benchmarks/imaging_ecal/options
CB_EMCAL_SCRIPT_DIR=benchmarks/imaging_ecal/scripts CB_EMCAL_SCRIPT_DIR=benchmarks/imaging_ecal/scripts
# Run Juggler # Run Juggler
gaudirun.py ${CB_EMCAL_OPTION_DIR}/imaging_topocluster.py xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
gaudirun.py ${CB_EMCAL_OPTION_DIR}/imaging_topocluster.py
# gaudirun.py ${CB_EMCAL_OPTION_DIR}/imaging_topocluster.py # gaudirun.py ${CB_EMCAL_OPTION_DIR}/imaging_topocluster.py
if [[ "$?" -ne "0" ]] ; then if [[ "$?" -ne "0" ]] ; then
echo "ERROR running juggler" echo "ERROR running juggler"
...@@ -125,20 +127,22 @@ for iev in "${ADDR[@]}"; do ...@@ -125,20 +127,22 @@ for iev in "${ADDR[@]}"; do
fi fi
python ${CB_EMCAL_SCRIPT_DIR}/draw_cluster_layers.py \ python ${CB_EMCAL_SCRIPT_DIR}/draw_cluster_layers.py \
${CB_EMCAL_REC_FILE} -e iev --topo-size=1.0 --compact=${CB_EMCAL_COMPACT_PATH} -o results/${particle} ${CB_EMCAL_REC_FILE} -e ${iev} --topo-size=1.0 --compact=${CB_EMCAL_COMPACT_PATH} -o results/${particle}
if [[ "$?" -ne "0" ]] ; then if [[ "$?" -ne "0" ]] ; then
echo "ERROR running analysis script: draw_cluster_layers" echo "ERROR running analysis script: draw_cluster_layers"
exit 1 exit 1
fi fi
python ${CB_EMCAL_SCRIPT_DIR}/draw_cluster.py \ python ${CB_EMCAL_SCRIPT_DIR}/draw_cluster.py \
${CB_EMCAL_REC_FILE} -e iev --topo-size=2.0 --compact=${CB_EMCAL_COMPACT_PATH} -o results/${particle} ${CB_EMCAL_REC_FILE} -e ${iev} --topo-size=2.0 --compact=${CB_EMCAL_COMPACT_PATH} -o results/${particle}
if [[ "$?" -ne "0" ]] ; then if [[ "$?" -ne "0" ]] ; then
echo "ERROR running analysis script: draw_cluster" echo "ERROR running analysis script: draw_cluster"
exit 1 exit 1
fi fi
done done
# TODO add this part back (solve mismatched EICD issue and add clusterreco back in option file)
":" << COMMENT
python ${CB_EMCAL_SCRIPT_DIR}/energy_profile.py \ python ${CB_EMCAL_SCRIPT_DIR}/energy_profile.py \
${CB_EMCAL_REC_FILE} --type=EM --energy=${CB_EMCAL_ENERGY} -o results/${particle} \ ${CB_EMCAL_REC_FILE} --type=EM --energy=${CB_EMCAL_ENERGY} -o results/${particle} \
--save=results/profile.csv --color=royalblue --save=results/profile.csv --color=royalblue
...@@ -146,7 +150,7 @@ if [[ "$?" -ne "0" ]] ; then ...@@ -146,7 +150,7 @@ if [[ "$?" -ne "0" ]] ; then
echo "ERROR running analysis script: energy_profile" echo "ERROR running analysis script: energy_profile"
exit 1 exit 1
fi fi
COMMENT
root_filesize=$(stat --format=%s "${CB_EMCAL_REC_FILE}") root_filesize=$(stat --format=%s "${CB_EMCAL_REC_FILE}")
if [[ "${CB_EMCAL_NUMEV}" -lt "500" ]] ; then if [[ "${CB_EMCAL_NUMEV}" -lt "500" ]] ; then
......
benchmarks/imaging_ecal/scripts/draw_cluster.py
+ 14
− 9
View file @ a9d56427
...@@ -114,10 +114,10 @@ if __name__ == '__main__': ...@@ -114,10 +114,10 @@ if __name__ == '__main__':
# we can read these values from xml file # we can read these values from xml file
desc = compact_constants(args.compact, [ desc = compact_constants(args.compact, [
'EcalBarrelAstroPix_RMin', 'EcalBarrel_rmin',
'EcalBarrelAstroPix_ReadoutLayerThickness', 'EcalBarrel_ReadoutLayerThickness',
'EcalBarrelAstroPix_ReadoutLayerNumber', 'EcalBarrel_ReadoutLayerNumber',
'EcalBarrelAstroPix_Length', 'EcalBarrel_length',
]) ])
if not len(desc): if not len(desc):
# or define Ecal shapes # or define Ecal shapes
...@@ -152,10 +152,13 @@ if __name__ == '__main__': ...@@ -152,10 +152,13 @@ if __name__ == '__main__':
# neutral particle, no need to consider magnetic field # neutral particle, no need to consider magnetic field
if np.isclose(inpart.Charge(), 0., rtol=1e-5): if np.isclose(inpart.Charge(), 0., rtol=1e-5):
vec = dfmcp[['px', 'py', 'pz']].values vec = dfmcp[['px', 'py', 'pz']].values
# charge particle, use the cluster center # charge particle, use center
# TODO implement motion of particles in fields
else: else:
flayer = df[df['layer'] == df['layer'].min()] # filter out outliers
vec = flayer[['x', 'y', 'z']].mean().values dfc = df[(df['eta'] <= df['eta'].quantile(0.95)) & (df['eta'] >= df['eta'].quantile(0.05)) &
(df['phi'] <= df['phi'].quantile(0.95)) & (df['phi'] >= df['phi'].quantile(0.05))]
vec = np.average(dfc[['x', 'y', 'z']].values, axis=0, weights=dfc['edep'].values)
vec = vec/np.linalg.norm(vec) vec = vec/np.linalg.norm(vec)
# particle line from (0, 0, 0) to the inner Ecal surface # particle line from (0, 0, 0) to the inner Ecal surface
...@@ -168,7 +171,8 @@ if __name__ == '__main__': ...@@ -168,7 +171,8 @@ if __name__ == '__main__':
fig = plt.figure(figsize=(15, 12), dpi=160) fig = plt.figure(figsize=(15, 12), dpi=160)
ax = fig.add_subplot(111, projection='3d') ax = fig.add_subplot(111, projection='3d')
# draw particle line # draw particle line
ax.plot(*pline[[2, 1]], '--', zs=pline[0], color='green') # TODO need to implement motion of particles in a field
# ax.plot(*pline[[2, 1]], '--', zs=pline[0], color='green')
# draw hits # draw hits
draw_hits3d(ax, df[['x', 'y', 'z', 'edep']].values, cmap, s=5.0) draw_hits3d(ax, df[['x', 'y', 'z', 'edep']].values, cmap, s=5.0)
draw_cylinder3d(ax, rmin, length, rstride=10, cstride=10, color='royalblue') draw_cylinder3d(ax, rmin, length, rstride=10, cstride=10, color='royalblue')
...@@ -184,7 +188,8 @@ if __name__ == '__main__': ...@@ -184,7 +188,8 @@ if __name__ == '__main__':
fig = plt.figure(figsize=(15, 12), dpi=160) fig = plt.figure(figsize=(15, 12), dpi=160)
ax = fig.add_subplot(111, projection='3d') ax = fig.add_subplot(111, projection='3d')
# draw particle line # draw particle line
ax.plot(*pline[[2, 1]], '--', zs=pline[0], color='green') # TODO need to implement motion of particles in a field
# ax.plot(*pline[[2, 1]], '--', zs=pline[0], color='green')
# draw hits # draw hits
draw_hits3d(ax, df[['x', 'y', 'z', 'edep']].values, cmap, s=20.0) draw_hits3d(ax, df[['x', 'y', 'z', 'edep']].values, cmap, s=20.0)
# draw_track_fit(ax, df, stop_layer=args.stop, # draw_track_fit(ax, df, stop_layer=args.stop,
......
benchmarks/imaging_ecal/scripts/draw_cluster_layers.py
+ 9
− 7
View file @ a9d56427
...@@ -68,10 +68,10 @@ if __name__ == '__main__': ...@@ -68,10 +68,10 @@ if __name__ == '__main__':
# we can read these values from xml file # we can read these values from xml file
desc = compact_constants(args.compact, [ desc = compact_constants(args.compact, [
'EcalBarrelAstroPix_RMin', 'EcalBarrel_rmin',
'EcalBarrelAstroPix_ReadoutLayerThickness', 'EcalBarrel_ReadoutLayerThickness',
'EcalBarrelAstroPix_ReadoutLayerNumber', 'EcalBarrel_ReadoutLayerNumber',
'EcalBarrelAstroPix_Length', 'EcalBarrel_length',
]) ])
if not len(desc): if not len(desc):
# or define Ecal shapes # or define Ecal shapes
...@@ -108,8 +108,9 @@ if __name__ == '__main__': ...@@ -108,8 +108,9 @@ if __name__ == '__main__':
vec = dfmcp[['px', 'py', 'pz']].values vec = dfmcp[['px', 'py', 'pz']].values
# charge particle, use the cluster center # charge particle, use the cluster center
else: else:
flayer = df[df['layer'] == df['layer'].min()] dfc = df[(df['eta'] <= df['eta'].quantile(0.95)) & (df['eta'] >= df['eta'].quantile(0.05)) &
vec = flayer[['x', 'y', 'z']].mean().values (df['phi'] <= df['phi'].quantile(0.95)) & (df['phi'] >= df['phi'].quantile(0.05))]
vec = np.average(dfc[['x', 'y', 'z']].values, axis=0, weights=dfc['edep'].values)
vec = vec/np.linalg.norm(vec) vec = vec/np.linalg.norm(vec)
# particle line from (0, 0, 0) to the inner Ecal surface # particle line from (0, 0, 0) to the inner Ecal surface
...@@ -124,7 +125,8 @@ if __name__ == '__main__': ...@@ -124,7 +125,8 @@ if __name__ == '__main__':
eta_rg = np.resize(-np.log(np.tan(vecp[0]/1000./2.)), 2) + np.asarray([-args.topo_range, args.topo_range])/1000. eta_rg = np.resize(-np.log(np.tan(vecp[0]/1000./2.)), 2) + np.asarray([-args.topo_range, args.topo_range])/1000.
os.makedirs(args.outdir, exist_ok=True) os.makedirs(args.outdir, exist_ok=True)
# cluster plot by layers (rebinned)
# cluster plot by layers (rebinned grids)
pdf = PdfPages(os.path.join(args.outdir, 'e{}_c{}_layers.pdf'.format(args.iev, args.icl))) pdf = PdfPages(os.path.join(args.outdir, 'e{}_c{}_layers.pdf'.format(args.iev, args.icl)))
bpos, bprops = (0.5, 0.95), dict(boxstyle='round', facecolor='wheat', alpha=0.2) bpos, bprops = (0.5, 0.95), dict(boxstyle='round', facecolor='wheat', alpha=0.2)
frames = [] frames = []
......
benchmarks/imaging_ecal/scripts/energy_profile.py
+ 2
− 2
View file @ a9d56427
...@@ -23,7 +23,7 @@ from utils import * ...@@ -23,7 +23,7 @@ from utils import *
def find_start_layer(grp, min_edep=0.5): def find_start_layer(grp, min_edep=0.5):
ids, edeps = grp.sort_values('layer')[['layer', 'edep']].values.T ids, edeps = grp.sort_values('layer')[['layer', 'edep']].values.T
edeps = np.cumsum(edeps) edeps = np.cumsum(edeps)
return min(ids[edeps > min_edep]) if sum(edeps > min_edep) > 0 else 20 return min(ids[edeps > min_edep]) if sum(edeps > min_edep) > 0 else -1
# execute this script # execute this script
...@@ -47,7 +47,7 @@ if __name__ == '__main__': ...@@ -47,7 +47,7 @@ if __name__ == '__main__':
dfe.loc[:, 'total_edep'] = dfe.groupby(['event', 'cluster'])['edep'].transform('sum') dfe.loc[:, 'total_edep'] = dfe.groupby(['event', 'cluster'])['edep'].transform('sum')
# dfe.loc[:, 'efrac'] = dfe['edep']/dfe['total_edep'] # dfe.loc[:, 'efrac'] = dfe['edep']/dfe['total_edep']
dfe.loc[:, 'efrac'] = dfe['edep']/(args.energy*1000.) dfe.loc[:, 'efrac'] = dfe['edep']/(args.energy*1000.)
dfe = dfe[dfe['cluster'] == 0] dfe = dfe[dfe['cluster'] == 1]
os.makedirs(args.outdir, exist_ok=True) os.makedirs(args.outdir, exist_ok=True)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment