From 398969165c7fafbad9e32a207e7e32a60ab03694 Mon Sep 17 00:00:00 2001 From: Sylvester Joosten <sjoosten@anl.gov> Date: Sat, 7 Aug 2021 02:07:28 +0000 Subject: [PATCH] updated for EICD --- .../imaging_ecal/options/hybrid_cluster.py | 22 +++++++++++-------- .../imaging_ecal/options/imaging_2dcluster.py | 6 +++-- .../options/imaging_topocluster.py | 10 +++++---- .../imaging_ecal/options/scfi_cluster.py | 12 +++++----- .../imaging_ecal/scripts/draw_cluster.py | 14 ++++++------ .../scripts/draw_cluster_layers.py | 6 ++--- .../imaging_ecal/scripts/energy_profile.py | 22 +++++++++---------- benchmarks/imaging_ecal/scripts/utils.py | 18 +++++++++------ 8 files changed, 62 insertions(+), 48 deletions(-) diff --git a/benchmarks/imaging_ecal/options/hybrid_cluster.py b/benchmarks/imaging_ecal/options/hybrid_cluster.py index e89cef98..5b2e1ead 100644 --- a/benchmarks/imaging_ecal/options/hybrid_cluster.py +++ b/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']) diff --git a/benchmarks/imaging_ecal/options/imaging_2dcluster.py b/benchmarks/imaging_ecal/options/imaging_2dcluster.py index e58e962f..9baee078 100644 --- a/benchmarks/imaging_ecal/options/imaging_2dcluster.py +++ b/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']) diff --git a/benchmarks/imaging_ecal/options/imaging_topocluster.py b/benchmarks/imaging_ecal/options/imaging_topocluster.py index 572544ef..9c715fcc 100644 --- a/benchmarks/imaging_ecal/options/imaging_topocluster.py +++ b/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 *"] diff --git a/benchmarks/imaging_ecal/options/scfi_cluster.py b/benchmarks/imaging_ecal/options/scfi_cluster.py index 873c3b89..3926326b 100644 --- a/benchmarks/imaging_ecal/options/scfi_cluster.py +++ b/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']) diff --git a/benchmarks/imaging_ecal/scripts/draw_cluster.py b/benchmarks/imaging_ecal/scripts/draw_cluster.py index 88f338d0..3c896f2c 100644 --- a/benchmarks/imaging_ecal/scripts/draw_cluster.py +++ b/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')) diff --git a/benchmarks/imaging_ecal/scripts/draw_cluster_layers.py b/benchmarks/imaging_ecal/scripts/draw_cluster_layers.py index 7baafebb..13a8eb1b 100644 --- a/benchmarks/imaging_ecal/scripts/draw_cluster_layers.py +++ b/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) diff --git a/benchmarks/imaging_ecal/scripts/energy_profile.py b/benchmarks/imaging_ecal/scripts/energy_profile.py index fff925ce..a8406f0c 100644 --- a/benchmarks/imaging_ecal/scripts/energy_profile.py +++ b/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) diff --git a/benchmarks/imaging_ecal/scripts/utils.py b/benchmarks/imaging_ecal/scripts/utils.py index da9954ac..d3847618 100644 --- a/benchmarks/imaging_ecal/scripts/utils.py +++ b/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): -- GitLab