Add plot script for e/pi separation performance

• update the json output from epcut and ML so rejection power and uncertainties can be directly read from them
• add a json output for run_benchmark.py to collect all run information
• add a script reading the run info files and plot the e/pi separation power curve

benchmarks/imaging_shower_ML/scripts/draw_imaging_event3d.py

@@ -13,7 +13,7 @@ import json
 import numpy as np
 import pandas as pd
 import argparse
-from utils import flatten_collection, cartesian_to_polar, imcal_info
+from utils import flatten_collection, cartesian_to_polar, imcal_info, center_of_gravity
 
 
 pd.set_option('display.max_rows', 500)
@@ -49,6 +49,28 @@ if __name__ == '__main__':
             default='EcalBarrelImagingHitsReco',
             help='name of data branch (edm4eic::CalorimeterHitCollection).'
             )
+    parser.add_argument(
+            '--eta-ngrid', type=int,
+            dest='eta_ngrid',
+            default=100,
+            help='number of eta grids for drawing.')
+    parser.add_argument(
+            '--eta-grid', type=float,
+            dest='eta_grid',
+            default=0.01,
+            help='eta grid size for drawing.')
+    parser.add_argument(
+            '--phi-ngrid', type=int,
+            dest='phi_ngrid',
+            default=100,
+            help='number of phi grids for drawing.')
+    parser.add_argument(
+            '--phi-grid', type=float,
+            dest='phi_grid',
+            default=0.01,
+            help='phi grid size for drawing (rad).')
+    parser.add_argument('--topo-range', type=float, default=50.0, dest='topo_range',
+                        help='half range for projection plot (mrad)')
     args = parser.parse_args()
 
     # read data and MCParticles
@@ -63,6 +85,12 @@ if __name__ == '__main__':
             'energy',
             'position.x', 'position.y', 'position.z',
             ])
-    print(data)
+    etot = data['energy'].sum()
+    xc, yc, zc = center_of_gravity(*data[['position.x', 'position.y', 'position.z']].values.T,
+                                   np.log(data['energy'].values/etot), 4.6)
+
+    rc, thetac, phic, r0c, etac = cartesian_to_polar(xc, yc, zc)
+    r, theta, phi, rc, eta = cartesian_to_polar(*data[['position.x', 'position.y', 'position.z']].values.T)
 
+    # determine eta and phi windows according to event (cluster) center