updated for EICD

benchmarks/imaging_ecal/options/hybrid_cluster.py

@@ -65,7 +65,7 @@ imcaldigi = CalHitDigi('imcal_digi',
         **imcaldaq)
 imcalreco = ImagingPixelReco('imcal_reco',
         # OutputLevel=DEBUG,
-        inputHitCollection='DigiEcalBarrelHits',
+        inputHitCollection=imcaldigi.outputHitCollection,
         outputHitCollection='RecoEcalBarrelHits',
         readoutClass='EcalBarrelHits',
         layerField='layer',
@@ -74,17 +74,19 @@ imcalreco = ImagingPixelReco('imcal_reco',
 
 imcalcluster = ImagingTopoCluster('imcal_cluster',
         # OutputLevel=DEBUG,
-        inputHitCollection='RecoEcalBarrelHits',
-        outputHitCollection='EcalBarrelClusterHits',
+        inputHitCollection=imcalreco.outputHitCollection,
+        outputProtoClusterCollection='EcalBarrelProtoClusters',
         localDistXY=[2.*mm, 2*mm],
         layerDistEtaPhi=[10*mrad, 10*mrad],
         neighbourLayersRange=2,
         sectorDist=3.*cm)
 clusterreco = ImagingClusterReco('imcal_clreco',
         # OutputLevel=DEBUG,
-        inputHitCollection='EcalBarrelClusterHits',
+        inputHitCollection=imcalcluster.inputHitCollection,
+        inputProtoClusterCollection=imcalcluster.outputProtoClusterCollection,
         outputLayerCollection='EcalBarrelClustersLayers',
         outputClusterCollection='EcalBarrelClusters',
+        outputInfoCollection='EcalBarrelClustersInfo',
         samplingFraction=kwargs['img_sf'])
 
 # scfi layers
@@ -100,7 +102,7 @@ scfi_barrel_digi = CalHitDigi("scfi_barrel_digi",
         **scfi_barrel_daq)
 
 scfi_barrel_reco = CalHitReco("scfi_barrel_reco",
-        inputHitCollection="EcalBarrelScFiHitsDigi",
+        inputHitCollection=scfi_barrel_digi.outputHitCollection,
         outputHitCollection="EcalBarrelScFiHitsReco",
         thresholdFactor=5.0,
         readoutClass="EcalBarrelScFiHits",
@@ -112,7 +114,7 @@ scfi_barrel_reco = CalHitReco("scfi_barrel_reco",
 # merge hits in different layer (projection to local x-y plane)
 scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger",
         # OutputLevel=DEBUG,
-        inputHitCollection="EcalBarrelScFiHitsReco",
+        inputHitCollection=scfi_barrel_reco.outputHitCollection,
         outputHitCollection="EcalBarrelScFiGridReco",
         fields=["fiber"],
         fieldRefNumbers=[1],
@@ -120,16 +122,18 @@ scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger",
 
 scfi_barrel_cl = IslandCluster("scfi_barrel_cl",
         # OutputLevel=DEBUG,
-        inputHitCollection="EcalBarrelScFiGridReco",
-        outputHitCollection="EcalBarrelScFiClusterHits",
+        inputHitCollection=scfi_barrel_merger.outputHitCollection,
+        outputProtoClusterCollection="EcalBarrelScFiProtoClusters",
         splitCluster=False,
         minClusterCenterEdep=10.*MeV,
         localDistXZ=[30*mm, 30*mm],
         sectorDist=5.*cm)
 
 scfi_barrel_clreco = RecoCoG("scfi_barrel_clreco",
-       inputHitCollection="EcalBarrelScFiClusterHits",
+       inputHitCollection=scfi_barrel_cl.inputHitCollection,
+       inputProtoClusterCollection=scfi_barrel_cl.outputProtoClusterCollection,
        outputClusterCollection="EcalBarrelScFiClusters",
+       outputInfoCollection="EcalBarrelScFiClustersInfo",
        logWeightBase=6.2,
        samplingFraction=kwargs['scfi_sf'])
 

benchmarks/imaging_ecal/options/imaging_2dcluster.py

@@ -64,7 +64,7 @@ imcal_barrel_merger = CalHitsProj('imcal_barrel_merger',
 imcal_barrel_cl = IslandCluster('imcal_barrel_cl',
         # OutputLevel=DEBUG,
         inputHitCollection=imcal_barrel_merger.outputHitCollection,
-        outputHitCollection='EcalBarrelClusterHits',
+        outputProtoClusterCollection='EcalBarrelProtoClusters',
         minClusterCenterEdep=1.*MeV,
         minClusterHitEdep=0.1*MeV,
         splitCluster=True,
@@ -72,8 +72,10 @@ imcal_barrel_cl = IslandCluster('imcal_barrel_cl',
 
 imcal_barrel_clreco = RecoCoG('imcal_barrel_clreco',
         # OutputLevel=DEBUG,
-        inputHitCollection=imcal_barrel_cl.outputHitCollection,
+        inputHitCollection=imcal_barrel_cl.inputHitCollection,
+        inputProtoClusterCollection = imcal_barrel_cl.outputProtoClusterCollection,
         outputClusterCollection='EcalBarrelClusters',
+        outputInfoCollection='EcalBarrelClustersInfo',
         logWeightBase=6.2,
         samplingFraction=kwargs['sf'])
 

benchmarks/imaging_ecal/options/imaging_topocluster.py

@@ -59,7 +59,7 @@ imcaldigi = CalorimeterHitDigi("imcal_digi",
         **daq_setting)
 imcalreco = ImagingPixelReco("imcal_reco",
         # OutputLevel=DEBUG,
-        inputHitCollection="DigiEcalBarrelHits",
+        inputHitCollection=imcaldigi.outputHitCollection,
         outputHitCollection="RecoEcalBarrelHits",
         readoutClass="EcalBarrelHits",
         layerField="layer",
@@ -68,8 +68,8 @@ imcalreco = ImagingPixelReco("imcal_reco",
 
 imcalcluster = ImagingTopoCluster("imcal_cluster",
         # OutputLevel=DEBUG,
-        inputHitCollection="RecoEcalBarrelHits",
-        outputHitCollection="EcalBarrelClusterHits",
+        inputHitCollection=imcalreco.outputHitCollection,
+        outputProtoClusterCollection="EcalBarrelProtoClusters",
         localDistXY=[2.*units.mm, 2*units.mm],
         layerDistEtaPhi=[10*units.mrad, 10*units.mrad],
         neighbourLayersRange=2,
@@ -77,9 +77,11 @@ imcalcluster = ImagingTopoCluster("imcal_cluster",
         minClusterNhits=5)
 clusterreco = ImagingClusterReco("imcal_clreco",
         # OutputLevel=DEBUG,
-        inputHitCollection="EcalBarrelClusterHits",
+        inputHitCollection=imcalcluster.inputHitCollection,
+        inputProtoClusterCollection=imcalcluster.outputProtoClusterCollection,
         outputLayerCollection="EcalBarrelClustersLayers",
         outputClusterCollection="EcalBarrelClusters",
+        outputInfoCollection="EcalBarrelClustersInfo",
         samplingFraction=sf)
 
 out.outputCommands = ["keep *"]

benchmarks/imaging_ecal/options/scfi_cluster.py

@@ -47,7 +47,7 @@ scfi_barrel_digi = CalHitDigi("scfi_barrel_digi",
         **scfi_barrel_daq)
 
 scfi_barrel_reco = CalHitReco("scfi_barrel_reco",
-        inputHitCollection="EcalBarrelScFiHitsDigi",
+        inputHitCollection=scfi_barrel_digi.outputHitCollection,
         outputHitCollection="EcalBarrelScFiHitsReco",
         thresholdFactor=5.0,
         readoutClass="EcalBarrelScFiHits",
@@ -59,7 +59,7 @@ scfi_barrel_reco = CalHitReco("scfi_barrel_reco",
 # merge hits in different layer (projection to local x-y plane)
 scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger",
         # OutputLevel=DEBUG,
-        inputHitCollection="EcalBarrelScFiHitsReco",
+        inputHitCollection=scfi_barrel_reco.outputHitCollection,
         outputHitCollection="EcalBarrelScFiGridReco",
         fields=["fiber"],
         fieldRefNumbers=[1],
@@ -67,15 +67,17 @@ scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger",
 
 scfi_barrel_cl = IslandCluster("scfi_barrel_cl",
         # OutputLevel=DEBUG,
-        inputHitCollection="EcalBarrelScFiGridReco",
-        outputHitCollection="EcalBarrelScFiClusterHits",
+        inputHitCollection=scfi_barrel_reco.outputHitCollection,
+        outputProtoClusterCollection="EcalBarrelScFiProtoClusters",
         splitCluster=False,
         minClusterCenterEdep=10.*MeV,
         localDistXZ=[30*mm, 30*mm])
 
 scfi_barrel_clreco = RecoCoG("scfi_barrel_clreco",
-       inputHitCollection="EcalBarrelScFiClusterHits",
+       inputHitCollection=scfi_barrel_cl.inputHitCollection,
+       inputProtoClusterCollection=scfi_barrel.outputProtoClusterCollection,
        outputClusterCollection="EcalBarrelScFiClusters",
+       outputInfoCollection="EcalBarrelScFiClustersInfo",
        logWeightBase=6.2,
        samplingFraction=kwargs['sf'])
 

benchmarks/imaging_ecal/scripts/draw_cluster.py

@@ -26,8 +26,8 @@ import sys
 # note z and x axes are switched
 def draw_hits3d(axis, data, cmap, units=('mm', 'mm', 'mm', 'MeV'), fontsize=24, **kwargs):
     # normalize energy to get colors
-    x, y, z, edep = np.transpose(data)
-    cvals = edep - min(edep) / (max(edep) - min(edep))
+    x, y, z, energy = np.transpose(data)
+    cvals = energy - min(energy) / (max(energy) - min(energy))
     cvals[np.isnan(cvals)] = 1.0
     colors = cmap(cvals)
 
@@ -37,7 +37,7 @@ def draw_hits3d(axis, data, cmap, units=('mm', 'mm', 'mm', 'MeV'), fontsize=24,
     axis.set_zlabel('x ({})'.format(units[2]), fontsize=fontsize + 2, labelpad=fontsize)
     axis.set_ylabel('y ({})'.format(units[1]), fontsize=fontsize + 2, labelpad=fontsize)
     axis.set_xlabel('z ({})'.format(units[0]), fontsize=fontsize + 2, labelpad=fontsize)
-    cb = plt.colorbar(cm.ScalarMappable(norm=matplotlib.colors.Normalize(vmin=min(edep), vmax=max(edep)), cmap=cmap),
+    cb = plt.colorbar(cm.ScalarMappable(norm=matplotlib.colors.Normalize(vmin=min(energy), vmax=max(energy)), cmap=cmap),
                       ax=axis, shrink=0.85)
     cb.ax.tick_params(labelsize=fontsize)
     cb.ax.get_yaxis().labelpad = fontsize
@@ -158,7 +158,7 @@ if __name__ == '__main__':
         # 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 = np.average(dfc[['x', 'y', 'z']].values, axis=0, weights=dfc['energy'].values)
     vec = vec/np.linalg.norm(vec)
 
     # particle line from (0, 0, 0) to the inner Ecal surface
@@ -174,7 +174,7 @@ if __name__ == '__main__':
     # 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_hits3d(ax, df[['x', 'y', 'z', 'energy']].values, cmap, s=5.0)
     draw_cylinder3d(ax, rmin, length, rstride=10, cstride=10, color='royalblue')
     draw_cylinder3d(ax, rmin + thickness, length, rstride=10, cstride=10, color='forestgreen')
     ax.set_zlim(-(rmin + thickness), rmin + thickness)
@@ -191,7 +191,7 @@ if __name__ == '__main__':
     # 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_hits3d(ax, df[['x', 'y', 'z', 'energy']].values, cmap, s=20.0)
     # draw_track_fit(ax, df, stop_layer=args.stop,
     #                scat_kw=dict(color='k', s=50.0), line_kw=dict(linestyle=':', color='k', lw=3))
     # view range
@@ -213,7 +213,7 @@ 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.
 
     fig, axs = plt.subplots(1, 2, figsize=(13, 12), dpi=160, gridspec_kw={'wspace':0., 'width_ratios': [12, 1]})
-    ax, sm = draw_heatmap(axs[0], df['eta'].values, df['phi'].values, weights=df['edep'].values,
+    ax, sm = draw_heatmap(axs[0], df['eta'].values, df['phi'].values, weights=df['energy'].values,
                           bins=(np.arange(*eta_rg, step=args.topo_size/1000.), np.arange(*phi_rg, step=args.topo_size)),
                           cmap=cmap, cmin=0., pc_kw=dict(alpha=0.8, edgecolor='k'))
 

benchmarks/imaging_ecal/scripts/draw_cluster_layers.py

@@ -110,7 +110,7 @@ if __name__ == '__main__':
     else:
         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 = np.average(dfc[['x', 'y', 'z']].values, axis=0, weights=dfc['energy'].values)
     vec = vec/np.linalg.norm(vec)
 
     # particle line from (0, 0, 0) to the inner Ecal surface
@@ -133,7 +133,7 @@ if __name__ == '__main__':
     for i in np.arange(1, df['layer'].max() + 1, dtype=int):
         data = df[df['layer'] == i]
         fig, axs = plt.subplots(1, 2, figsize=(17, 16), dpi=160, gridspec_kw={'wspace':0., 'width_ratios': [16, 1]})
-        ax, sm = draw_heatmap(axs[0], data['eta'].values, data['phi'].values, weights=data['edep'].values,
+        ax, sm = draw_heatmap(axs[0], data['eta'].values, data['phi'].values, weights=data['energy'].values,
                               bins=(np.arange(*eta_rg, step=args.topo_size/1000.), np.arange(*phi_rg, step=args.topo_size)),
                               cmap=cmap, vmin=0.0, vmax=1.0, pc_kw=dict(alpha=0.8, edgecolor='k'))
         ax.set_ylabel(r'$\phi$ (mrad)', fontsize=28)
@@ -170,7 +170,7 @@ if __name__ == '__main__':
         data = df[df['layer'] == i]
         fig, axs = plt.subplots(1, 2, figsize=(17, 16), dpi=160, gridspec_kw={'wspace':0., 'width_ratios': [16, 1]})
         ax = axs[0]
-        colors = cmap(data['edep']/1.0)
+        colors = cmap(data['energy']/1.0)
         ax.scatter(data['eta'].values, data['phi'].values, c=colors, marker='s', s=15.0)
         ax.set_xlim(*eta_rg)
         ax.set_ylim(*phi_rg)

benchmarks/imaging_ecal/scripts/energy_profile.py

@@ -20,10 +20,10 @@ import sys
 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 -1
+def find_start_layer(grp, min_energy=0.5):
+    ids, energys = grp.sort_values('layer')[['layer', 'energy']].values.T
+    energys = np.cumsum(energys)
+    return min(ids[energys > min_energy]) if sum(energys > min_energy) > 0 else -1
 
 
 # execute this script
@@ -44,9 +44,9 @@ if __name__ == '__main__':
     load_root_macros(args.macros)
     # prepare data
     dfe = get_layers_data(args.file, branch=args.branch)
-    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.loc[:, 'total_energy'] = dfe.groupby(['event', 'cluster'])['energy'].transform('sum')
+    # dfe.loc[:, 'efrac'] = dfe['energy']/dfe['total_energy']
+    dfe.loc[:, 'efrac'] = dfe['energy']/(args.energy*1000.)
     # dfe = dfe[dfe['cluster'] == 0]
 
     os.makedirs(args.outdir, exist_ok=True)
@@ -54,7 +54,7 @@ if __name__ == '__main__':
     slayer = dfe.groupby('event').apply(lambda x: find_start_layer(x, 1.0)).astype(int)
     dfe = dfe.merge(slayer.to_frame(name='slayer'), on='event')
     # dfe.loc[:, 'eff_layer'] = dfe['layer'] - dfe['slayer']
-    # prof = dfe[dfe['eff_layer'] > 0].groupby('eff_layer')['edep'].describe()
+    # prof = dfe[dfe['eff_layer'] > 0].groupby('eff_layer')['energy'].describe()
     prof = dfe.groupby('layer')['efrac'].describe()
 
     # print(prof['mean'])
@@ -73,7 +73,7 @@ if __name__ == '__main__':
     ax.set_ylabel('Normalized Counts', fontsize=26)
     ax.text(*bpos, 'Mininum Edep\n' + '{:.1f} MeV'.format(1.0),
             transform=ax.transAxes, fontsize=26, verticalalignment='top', bbox=bprops)
-    fig.savefig(os.path.join(args.outdir, 'edep_start_{}_{}.png'.format(args.type, int(args.energy))))
+    fig.savefig(os.path.join(args.outdir, 'energy_start_{}_{}.png'.format(args.type, int(args.energy))))
 
 
     fig, ax = plt.subplots(figsize=(16, 9), dpi=160)
@@ -90,7 +90,7 @@ if __name__ == '__main__':
     ax.set_ylabel('Energy Deposit Percentage', fontsize=26)
     fig.savefig(os.path.join(args.outdir, 'efrac_{}_{}.png'.format(args.type, int(args.energy))))
 
-    # edep fraction by layers
+    # energy fraction by layers
     layers = np.asarray([
         [1, 5, 8,],
         [10, 15, 20],
@@ -111,7 +111,7 @@ if __name__ == '__main__':
         ax.grid(linestyle=':', which='both')
         # ax.set_xlabel('Energy Deposit (MeV)', fontsize=26)
         # ax.set_ylabel('Normalized Counts', fontsize=26)
-        mean, std = data.describe().loc[['mean', 'std'], 'edep'].values
+        mean, std = data.describe().loc[['mean', 'std'], 'energy'].values
         label = 'Layer {}'.format(layer)
     #            + '\n' + r'$\mu = {:.3f}$ MeV'.format(mean)
     #            + '\n' + r'$\sigma = {:.3f}$ MeV'.format(std)

benchmarks/imaging_ecal/scripts/utils.py

@@ -98,10 +98,12 @@ def get_hits_data(path, evnums=None, branch='RecoEcalBarrelHits'):
 
         events.GetEntry(iev)
         for hit in getattr(events, branch):
-            dbuf[idb] = (iev, hit.clusterID, hit.layerID, hit.position.x, hit.position.y, hit.position.z, hit.edep*1000.)
+            dbuf[idb] = (iev, hit.ID.value, hit.layer, hit.position.x, hit.position.y,
+                    hit.position.z, hit.energy*1000.)
             idb += 1
 
-    return pd.DataFrame(data=dbuf[:idb], columns=['event', 'cluster', 'layer', 'x', 'y', 'z', 'edep'])
+    return pd.DataFrame(data=dbuf[:idb], columns=['event', 'cluster', 'layer', 'x', 'y',
+        'z', 'energy'])
 
 
 # read layers data from root file
@@ -122,11 +124,13 @@ def get_layers_data(path, evnums=None, branch="EcalBarrelClustersLayers"):
 
         events.GetEntry(iev)
         for layer in getattr(events, branch):
-            dbuf[idb] = (iev, layer.clusterID, layer.layerID,
-                         layer.position.x, layer.position.y, layer.position.z, layer.edep*1000.)
+            dbuf[idb] = (iev, layer.clusterID.value, layer.layer,
+                         layer.position.x, layer.position.y, layer.position.z,
+                         layer.energy*1000.)
             idb += 1
 
-    return pd.DataFrame(data=dbuf[:idb], columns=['event', 'cluster', 'layer', 'x', 'y', 'z', 'edep'])
+    return pd.DataFrame(data=dbuf[:idb], columns=['event', 'cluster', 'layer', 'x', 'y',
+        'z', 'energy'])
 
 
 # read clusters data from root file
@@ -147,10 +151,10 @@ def get_clusters_data(path, evnums=None, branch='EcalBarrelClustersReco'):
 
         events.GetEntry(iev)
         for k, cl in enumerate(getattr(events, branch)):
-            dbuf[idb] = (iev, k, cl.nhits, cl.edep*1000., cl.cl_theta, cl.cl_phi)
+            dbuf[idb] = (iev, k, cl.nhits, cl.energy*1000., cl.cl_theta, cl.cl_phi)
             idb += 1
 
-    return pd.DataFrame(data=dbuf[:idb], columns=['event', 'cluster', 'nhits', 'edep', 'cl_theta', 'cl_phi'])
+    return pd.DataFrame(data=dbuf[:idb], columns=['event', 'cluster', 'nhits', 'energy', 'cl_theta', 'cl_phi'])
 
 
 def compact_constants(path, names):