Analysis script for checking merged hits for imaging barrel ecal

Add an analysis script for checking merged hits for imaging barrel ecal

Update some scripts to make sure they are working with the latest changes in detectors/algorithms.

benchmarks/imaging_ecal/scripts/check_merged_hits.py

@@ -42,12 +42,11 @@ if __name__ == '__main__':
     print('{:20s} {:#066b}'.format('Merging Mask', mmask))
 
     # build dataframes from reconstructed root files
-    '''
     # get hits
     columns = ['cellID', 'position.x', 'position.y', 'position.z', 'energy']
     rdf_rec = ROOT.RDataFrame('events', args.rec_file)
     # raw hits
-    dfr = flatten_collection(rdf_rec, args.rcoll, columns)
+    dfr = flatten_collection(rdf_rec.Range(100), args.rcoll, columns)
     dfr.rename(columns={c: c.replace(args.rcoll + '.', '') for c in dfr.columns}, inplace=True)
     # merged hits
     dfm = flatten_collection(rdf_rec, args.mcoll, columns)
@@ -65,11 +64,11 @@ if __name__ == '__main__':
     print(dfr[['masked_id', 'merged_id', 'raw_merge_num']])
     print(dfm[['masked_id', 'merged_id', 'merge_raw_num']])
 
-    dfr.to_csv('raw_hits.csv')
-    dfm.to_csv('merged_hits.csv')
-    '''
+    # dfr.to_csv('raw_hits.csv')
+    # dfm.to_csv('merged_hits.csv')
 
-    dfm = pd.read_csv('merged_hits.csv', index_col=0)
+    # merged hits (grid)
+    # dfm = pd.read_csv('merged_hits.csv', index_col=0)
     ev_grpm = dfm.groupby('event')
     # energy scaling (assume mean at 5 GeV)
     escale = 5000./np.mean(ev_grpm['energy'].sum())
@@ -105,12 +104,14 @@ if __name__ == '__main__':
     fig.savefig('signals_per_grid.png')
 
 
-    dfr = pd.read_csv('raw_hits.csv', index_col=0)
+    # raw hits (fiber)
+    # dfr = pd.read_csv('raw_hits.csv', index_col=0)
+    dfr['energy'] = dfr['energy']*escale
     ev_grpr = dfr.groupby('event')
     ev_nfibers = ev_grpr['energy'].size()
     nev = len(ev_nfibers)
 
-    # number of grids per event
+    # number of fiber per event
     fig, ax = plt.subplots(figsize=(12, 9), dpi=160)
     ax.hist(ev_nfibers, ec='k', bins=30, rwidth=1.0, weights=np.repeat(1./nev, nev))
     ax.tick_params(labelsize=24)
@@ -119,12 +120,12 @@ if __name__ == '__main__':
     ax.set_ylabel('Normalized Counts', fontsize=24)
     fig.savefig('fiber_per_event.png')
 
-    # energy deposit per grid
+    # energy deposit per fiber
     fig, ax = plt.subplots(figsize=(12, 9), dpi=160)
     ax.hist(dfr['energy'], ec='k', bins=30, rwidth=1.0, weights=np.repeat(1./len(dfr), len(dfr)))
     ax.tick_params(labelsize=24)
     ax.set_title('5.0 GeV $e^-$', fontsize=26)
-    ax.set_xlabel('Grid $E_{dep}$ (MeV)', fontsize=24)
+    ax.set_xlabel('Fiber $E_{dep}$ (MeV)', fontsize=24)
     ax.set_ylabel('Normalized Counts', fontsize=24)
     fig.savefig('edep_per_fiber.png')
 
@@ -134,7 +135,23 @@ if __name__ == '__main__':
             rwidth=1.0, weights=np.repeat(1./len(dfr), len(dfr)))
     ax.tick_params(labelsize=24)
     ax.set_title('5.0 GeV $e^-$', fontsize=26)
-    ax.set_xlabel('# of Fiber-Signals per Grid', fontsize=24)
+    ax.set_xlabel('# of Grids per Fiber-Signal', fontsize=24)
     ax.set_ylabel('Normalized Counts', fontsize=24)
     fig.savefig('grids_per_fiber.png')
 
+    # build relations between two datasets
+    dfe = dfr.groupby(['event', 'masked_id'])['energy'].sum().reset_index()
+    print(dfe)
+    dfe.rename(columns={'energy':'energy_fibers'}, inplace=True)
+    dfem = dfm[['event', 'merged_id', 'energy']].rename(columns={'merged_id': 'masked_id'})
+    dfe = dfe.merge(dfem, on=['event', 'masked_id'])
+
+    fig, ax = plt.subplots(figsize=(12, 9), dpi=160)
+    ax.scatter(*dfe[['energy_fibers', 'energy']].values.T)
+    ax.tick_params(labelsize=24)
+    ax.set_title('5.0 GeV $e^-$ (100 events)', fontsize=26)
+    ax.set_xlabel('Fiber Group Energy Sum (MeV)', fontsize=24)
+    ax.set_ylabel('Grid Energy (MeV)', fontsize=24)
+    fig.savefig('fiber_grid_energy.png')
+
+