diff --git a/benchmarks/imaging_ecal/options/imaging_topocluster.py b/benchmarks/imaging_ecal/options/imaging_topocluster.py
index bdcda122e517fe40ba67afa59c598c8e76b4457a..cdcbfc0238372c79c450035cfda1564532f4d880 100644
--- a/benchmarks/imaging_ecal/options/imaging_topocluster.py
+++ b/benchmarks/imaging_ecal/options/imaging_topocluster.py
@@ -31,7 +31,6 @@ 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(','))
podioevent = EICDataSvc("EventDataSvc", inputs=kwargs['input'].split(','), OutputLevel=DEBUG)
out = PodioOutput("out", filename=kwargs['output'])
@@ -64,13 +63,15 @@ imcalreco = ImagingPixelReco("imcal_reco",
readoutClass="EcalBarrelHits",
layerField="layer",
sectorField="module")
-imcalcluster = ImagingTopoCluster(inputHitCollection="RecoEcalBarrelHits",
+imcalcluster = ImagingTopoCluster("imcal_cluster",
+ inputHitCollection="RecoEcalBarrelHits",
outputClusterCollection="EcalBarrelClusters",
localRanges=[2.*units.mm, 2*units.mm],
adjLayerRanges=[10*units.mrad, 10*units.mrad],
adjLayerDiff=2,
adjSectorDist=3.*units.cm)
-clusterreco = ImagingClusterReco(inputClusterCollection="EcalBarrelClusters",
+clusterreco = ImagingClusterReco("imcal_clreco",
+ inputClusterCollection="EcalBarrelClusters",
outputLayerCollection="EcalBarrelClustersLayers",
samplingFraction=sf,
OutputLevel=DEBUG)
@@ -78,7 +79,7 @@ clusterreco = ImagingClusterReco(inputClusterCollection="EcalBarrelClusters",
out.outputCommands = ["keep *"]
ApplicationMgr(
- TopAlg=[podioinput, copier, calcopier, imcaldigi, imcalreco, imcalcluster, clusterreco, out],
+ TopAlg=[podioinput, copier, calcopier, imcaldigi, imcalreco, imcalcluster, out],
EvtSel='NONE',
EvtMax=kwargs['nev'],
ExtSvc=[podioevent],
diff --git a/benchmarks/imaging_ecal/run_emcal_barrel.sh b/benchmarks/imaging_ecal/run_emcal_barrel.sh
index ad7f7ed89fdea5a3ec743c878e51c4f54cc4ab6c..43cc78c3fcf3523fb8da5241ee62cdf2f0b9b5f2 100644
--- a/benchmarks/imaging_ecal/run_emcal_barrel.sh
+++ b/benchmarks/imaging_ecal/run_emcal_barrel.sh
@@ -5,11 +5,11 @@ print_env.sh
function print_the_help {
echo "USAGE: ${0} -n <nevents> -t <nametag> -p <particle> "
echo " OPTIONS: "
- echo " -n,--nevents Number of eventsj"
+ echo " -n,--nevents Number of events"
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"
- exit
+ exit
}
POSITIONAL=()
@@ -33,9 +33,9 @@ do
shift # past value
;;
-n|--nevents)
- export CB_EMCAL_NUMEV=${OPTARG}
+ export CB_EMCAL_NUMEV="$2"
shift # past argument
- #shift # past value
+ shift # past value
;;
*) # unknown option
#POSITIONAL+=("$1") # save it in an array for later
@@ -54,7 +54,7 @@ if [[ ! -n "${CB_EMCAL_NUMEV}" ]] ; then
fi
if [[ ! -n "${CB_EMCAL_IEV}" ]] ; then
- export CB_EMCAL_IEV="0, 1, 2, 3, 4"
+ export CB_EMCAL_IEV="0, 1"
fi
if [[ ! -n "${CB_EMCAL_ENERGY}" ]] ; then
@@ -87,6 +87,7 @@ 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} \
@@ -108,7 +109,8 @@ CB_EMCAL_OPTION_DIR=benchmarks/imaging_ecal/options
CB_EMCAL_SCRIPT_DIR=benchmarks/imaging_ecal/scripts
# 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
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running juggler"
@@ -125,20 +127,22 @@ for iev in "${ADDR[@]}"; do
fi
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
echo "ERROR running analysis script: draw_cluster_layers"
exit 1
fi
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
echo "ERROR running analysis script: draw_cluster"
exit 1
fi
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 \
${CB_EMCAL_REC_FILE} --type=EM --energy=${CB_EMCAL_ENERGY} -o results/${particle} \
--save=results/profile.csv --color=royalblue
@@ -146,7 +150,7 @@ if [[ "$?" -ne "0" ]] ; then
echo "ERROR running analysis script: energy_profile"
exit 1
fi
-
+COMMENT
root_filesize=$(stat --format=%s "${CB_EMCAL_REC_FILE}")
if [[ "${CB_EMCAL_NUMEV}" -lt "500" ]] ; then
diff --git a/benchmarks/imaging_ecal/scripts/draw_cluster.py b/benchmarks/imaging_ecal/scripts/draw_cluster.py
index 315f689174d497801b285500cad1cf2b249f3faa..d03db476f5270ed0ef3e564cea2fbba1fcf2cc54 100644
--- a/benchmarks/imaging_ecal/scripts/draw_cluster.py
+++ b/benchmarks/imaging_ecal/scripts/draw_cluster.py
@@ -114,10 +114,10 @@ if __name__ == '__main__':
# we can read these values from xml file
desc = compact_constants(args.compact, [
- 'EcalBarrelAstroPix_RMin',
- 'EcalBarrelAstroPix_ReadoutLayerThickness',
- 'EcalBarrelAstroPix_ReadoutLayerNumber',
- 'EcalBarrelAstroPix_Length',
+ 'EcalBarrel_rmin',
+ 'EcalBarrel_ReadoutLayerThickness',
+ 'EcalBarrel_ReadoutLayerNumber',
+ 'EcalBarrel_length',
])
if not len(desc):
# or define Ecal shapes
@@ -152,10 +152,13 @@ if __name__ == '__main__':
# neutral particle, no need to consider magnetic field
if np.isclose(inpart.Charge(), 0., rtol=1e-5):
vec = dfmcp[['px', 'py', 'pz']].values
- # charge particle, use the cluster center
+ # charge particle, use center
+ # TODO implement motion of particles in fields
else:
- flayer = df[df['layer'] == df['layer'].min()]
- vec = flayer[['x', 'y', 'z']].mean().values
+ # filter out outliers
+ 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)
# particle line from (0, 0, 0) to the inner Ecal surface
@@ -168,7 +171,8 @@ if __name__ == '__main__':
fig = plt.figure(figsize=(15, 12), dpi=160)
ax = fig.add_subplot(111, projection='3d')
# 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_hits3d(ax, df[['x', 'y', 'z', 'edep']].values, cmap, s=5.0)
draw_cylinder3d(ax, rmin, length, rstride=10, cstride=10, color='royalblue')
@@ -184,7 +188,8 @@ if __name__ == '__main__':
fig = plt.figure(figsize=(15, 12), dpi=160)
ax = fig.add_subplot(111, projection='3d')
# 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_hits3d(ax, df[['x', 'y', 'z', 'edep']].values, cmap, s=20.0)
# draw_track_fit(ax, df, stop_layer=args.stop,
diff --git a/benchmarks/imaging_ecal/scripts/draw_cluster_layers.py b/benchmarks/imaging_ecal/scripts/draw_cluster_layers.py
index fa4a96c4a028e8ce67a34e86bfba679bd1a6560e..25f102e0296cf9ec98a44c9598b640be6c5bd00b 100644
--- a/benchmarks/imaging_ecal/scripts/draw_cluster_layers.py
+++ b/benchmarks/imaging_ecal/scripts/draw_cluster_layers.py
@@ -68,10 +68,10 @@ if __name__ == '__main__':
# we can read these values from xml file
desc = compact_constants(args.compact, [
- 'EcalBarrelAstroPix_RMin',
- 'EcalBarrelAstroPix_ReadoutLayerThickness',
- 'EcalBarrelAstroPix_ReadoutLayerNumber',
- 'EcalBarrelAstroPix_Length',
+ 'EcalBarrel_rmin',
+ 'EcalBarrel_ReadoutLayerThickness',
+ 'EcalBarrel_ReadoutLayerNumber',
+ 'EcalBarrel_length',
])
if not len(desc):
# or define Ecal shapes
@@ -108,8 +108,9 @@ if __name__ == '__main__':
vec = dfmcp[['px', 'py', 'pz']].values
# charge particle, use the cluster center
else:
- flayer = df[df['layer'] == df['layer'].min()]
- 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)
# particle line from (0, 0, 0) to the inner Ecal surface
@@ -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.
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)))
bpos, bprops = (0.5, 0.95), dict(boxstyle='round', facecolor='wheat', alpha=0.2)
frames = []
diff --git a/benchmarks/imaging_ecal/scripts/energy_profile.py b/benchmarks/imaging_ecal/scripts/energy_profile.py
index c984cfbd572c04cca40a73bedc27b4a72100be10..234e646121657f6e958433f56f64aaa1a796cf62 100644
--- a/benchmarks/imaging_ecal/scripts/energy_profile.py
+++ b/benchmarks/imaging_ecal/scripts/energy_profile.py
@@ -23,7 +23,7 @@ from utils import *
def find_start_layer(grp, min_edep=0.5):
ids, edeps = grp.sort_values('layer')[['layer', 'edep']].values.T
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
@@ -47,7 +47,7 @@ if __name__ == '__main__':
dfe.loc[:, 'total_edep'] = dfe.groupby(['event', 'cluster'])['edep'].transform('sum')
# dfe.loc[:, 'efrac'] = dfe['edep']/dfe['total_edep']
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)