Skip to content
Snippets Groups Projects
Commit 22b07fd9 authored by Chao Peng's avatar Chao Peng
Browse files

remove compact files from barrel emcal benchmarks

parent 6f768f56
No related branches found
No related tags found
1 merge request!81remove compact files from barrel emcal benchmarks
Hide whitespace changes
Inline Side-by-side
benchmarks/sampling_ecal/compact/barrel_ecal_test.xml deleted 100644 → 0
+ 0
− 129
View file @ 6f768f56
<?xml version="1.0" encoding="UTF-8"?>
<lccdd>
<info
name="topside"
title="Time-of-flight Optimized PID Silicon Detector for EIC (TOPSiDE)"
author="Whitney R. Armstrong"
url=""
status="development"
version="$Id: 1$">
<comment> </comment>
</info>
<define>
<include ref="topside/topside_defs.xml" />
<include ref="eic/eic_defs.xml" />
</define>
<properties>
<matrix name="RINDEX__Vacuum" coldim="2" values="
1.0*eV 1.0
5.1*eV 1.0
"/>
<matrix name="RINDEX__Air" coldim="2" values="
1.0*eV 1.00029
5.1*eV 1.00029
"/>
<matrix name="RINDEX__Quartz" coldim="2" values="
1.0*eV 1.46
5.1*eV 1.46
"/>
<matrix name="RINDEX__N2" coldim="2" values="
1.0*eV 1.00033
4.0*eV 1.00033
5.1*eV 1.00033
"/>
<matrix name="RINDEX__Pyrex" coldim="2" values="
1.0*eV 1.5
4.0*eV 1.5
5.1*eV 1.5
"/>
<matrix name= "REFLECTIVITY_mirror" coldim="2" values="
1.0*eV 0.9
4.0*eV 0.9
5.1*eV 0.9
"/>
</properties>
<includes>
<gdmlFile ref="topside/elements.xml"/>
<gdmlFile ref="topside/materials.xml"/>
</includes>
<materials>
<material name="AirOptical">
<D type="density" unit="g/cm3" value="0.0012"/>
<fraction n="0.754" ref="N"/>
<fraction n="0.234" ref="O"/>
<fraction n="0.012" ref="Ar"/>
<property name="RINDEX" ref="RINDEX__Vacuum"/>
</material>
<material name="N2cherenkov">
<D type="density" value="0.00125" unit="g/cm3"/>
<composite n="1" ref="N"/>
<property name="RINDEX" ref="RINDEX__N2"/>
</material>
<material name="PyrexGlass">
<D type="density" value="2.23" unit="g/cm3"/>
<fraction n="0.806" ref="SiliconOxide"/>
<fraction n="0.130" ref="BoronOxide"/>
<fraction n="0.040" ref="SodiumOxide"/>
<fraction n="0.023" ref="AluminumOxide"/>
<property name="RINDEX" ref="RINDEX__Pyrex"/>
</material>
</materials>
<surfaces>
<comment> For the values of "finish", model and type, see TGeoOpticalSurface.h !
</comment>
<opticalsurface finish="polished" model="glisur" name="MirrorOpticalSurface" type="dielectric_metal" value="0">
<property name="REFLECTIVITY" ref="REFLECTIVITY_mirror"/>
<property name="RINDEX" coldim="2" values="1.034*eV 1.5 4.136*eV 1.5"/>
<!--<property name="EFFICIENCY" ref="EFFICIENCY0x8b77240"/>-->
</opticalsurface>
<opticalsurface name="mirror2" finish="polished" model="glisur" type="dielectric_dielectric">
<property name="REFLECTIVITY" coldim="2" values="1.034*eV 0.8 4.136*eV 0.9"/>
<property name="EFFICIENCY" coldim="2" values="2.034*eV 0.8 4.136*eV 1.0"/>
<property name="RINDEX" coldim="2" values="1.034*eV 1.5 4.136*eV 1.5"/>
</opticalsurface>
</surfaces>
<limits>
<limitset name="EICBeamlineLimits">
<limit name="step_length_max" particles="*" value="1.0" unit="mm" />
<limit name="track_length_max" particles="*" value="1.0" unit="mm" />
<limit name="time_max" particles="*" value="0.1" unit="ns" />
<limit name="ekin_min" particles="*" value="0.001" unit="MeV" />
<limit name="range_min" particles="*" value="0.1" unit="mm" />
</limitset>
<limitset name="cal_limits">
<limit name="step_length_max" particles="*" value="0.1" unit="mm"/>
</limitset>
</limits>
<display>
<include ref="topside/display.xml" />
</display>
<!--
<include ref="topside/vertex_tracker.xml"/>
-->
<include ref="topside/beampipe.xml"/>
<include ref="topside/silicon_tracker.xml"/>
<!--<include ref="topside/ecal.xml"/>--> <comment> old version of em barrel - SiW sampling design </comment>
<include ref="topside/ecal_wAstroPixSiW.xml"/> <comment> new version of em barrel - SiW AstroPix sampling design </comment>
<!--
<include ref="topside/hcal.xml"/>
<include ref="topside/forward_rich.xml"/>
<include ref="topside/solenoid.xml"/>
<include ref="topside/roman_pots.xml"/>
<include ref="eic/forward_ion_beamline.xml"/>
-->
<detectors>
</detectors>
<readouts>
</readouts>
</lccdd>
benchmarks/sampling_ecal/run_emcal_barrel.sh
+ 1
− 7
View file @ 22b07fd9
...@@ -8,13 +8,7 @@ do ...@@ -8,13 +8,7 @@ do
esac esac
done done
if [[ ! -n "${CB_EMCAL_TEST_DETECTOR}" ]] ; then export CB_EMCAL_COMPACT_PATH=${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml
export CB_EMCAL_TEST_DETECTOR="barrel_ecal_test"
fi
# copy the compact file to detector path
cp benchmarks/sampling_ecal/compact/${CB_EMCAL_TEST_DETECTOR}.xml ${DETECTOR_PATH}
export CB_EMCAL_COMPACT_PATH=${DETECTOR_PATH}/${CB_EMCAL_TEST_DETECTOR}.xml
if [[ ! -n "${CB_EMCAL_NUMEV}" ]] ; then if [[ ! -n "${CB_EMCAL_NUMEV}" ]] ; then
export CB_EMCAL_NUMEV=1000 export CB_EMCAL_NUMEV=1000
......
benchmarks/sampling_ecal/scripts/draw_cluster.py
+ 41
− 26
View file @ 22b07fd9
...@@ -107,6 +107,8 @@ if __name__ == '__main__': ...@@ -107,6 +107,8 @@ if __name__ == '__main__':
help='bin size for projection plot (mrad)') help='bin size for projection plot (mrad)')
parser.add_argument('--topo-range', type=float, default=50.0, dest='topo_range', parser.add_argument('--topo-range', type=float, default=50.0, dest='topo_range',
help='half range for projection plot (mrad)') help='half range for projection plot (mrad)')
parser.add_argument('--zoom-factor', type=float, default=1.0, dest='zoom_factor',
help='factor for zoom-in')
args = parser.parse_args() args = parser.parse_args()
...@@ -115,7 +117,7 @@ if __name__ == '__main__': ...@@ -115,7 +117,7 @@ if __name__ == '__main__':
'EcalBarrelAstroPix_RMin', 'EcalBarrelAstroPix_RMin',
'EcalBarrelAstroPix_ReadoutLayerThickness', 'EcalBarrelAstroPix_ReadoutLayerThickness',
'EcalBarrelAstroPix_ReadoutLayerNumber', 'EcalBarrelAstroPix_ReadoutLayerNumber',
'EcalBarrelAstroPix_Length' 'EcalBarrelAstroPix_Length',
]) ])
if not len(desc): if not len(desc):
# or define Ecal shapes # or define Ecal shapes
...@@ -130,22 +132,39 @@ if __name__ == '__main__': ...@@ -130,22 +132,39 @@ if __name__ == '__main__':
# read data # read data
load_root_macros(args.macros) load_root_macros(args.macros)
df = get_hits_data(args.file, args.iev, branch=args.branch) df = get_hits_data(args.file, args.iev, branch=args.branch)
dfmcp = get_mcp_data(args.file, args.iev, 'mcparticles2')
vec = dfmcp.loc[dfmcp['status'] == 24578, ['px', 'py', 'pz']].iloc[0].values
vec = vec/np.linalg.norm(vec)
df = df[df['cluster'] == args.icl] df = df[df['cluster'] == args.icl]
if not len(df): if not len(df):
print("Error: do not find any hits for cluster {:d} in event {:d}".format(args.icl, args.iev)) print("Error: do not find any hits for cluster {:d} in event {:d}".format(args.icl, args.iev))
exit(-1) exit(-1)
# convert to polar coordinates (mrad), and stack all r values
df['r'] = np.sqrt(df['x'].values**2 + df['y'].values**2 + df['z'].values**2)
df['phi'] = np.arctan2(df['y'].values, df['x'].values)*1000.
df['theta'] = np.arccos(df['z'].values/df['r'].values)*1000.
df['eta'] = -np.log(np.tan(df['theta'].values/1000./2.))
# truth
dfmcp = get_mcp_simple(args.file, args.iev, 'mcparticles2').iloc[0]
pdgbase = ROOT.TDatabasePDG()
inpart = pdgbase.GetParticle(int(dfmcp['pid']))
print("Incoming particle = {}, pdgcode = {}, charge = {}, mass = {}"\
.format(inpart.GetName(), inpart.PdgCode(), inpart.Charge(), inpart.Mass()))
# 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
else:
flayer = df[df['layer'] == df['layer'].min()]
vec = flayer[['x', 'y', 'z']].mean().values
vec = vec/np.linalg.norm(vec)
# particle line from (0, 0, 0) to the inner Ecal surface # particle line from (0, 0, 0) to the inner Ecal surface
length = rmin/np.sqrt(vec[0]**2 + vec[1]**2) length = rmin/np.sqrt(vec[0]**2 + vec[1]**2)
pline = np.transpose(vec*np.mgrid[0:length:2j][:, np.newaxis]) pline = np.transpose(vec*np.mgrid[0:length:2j][:, np.newaxis])
cmap = truncate_colormap(plt.get_cmap('jet'), 0.1, 0.9)
os.makedirs(args.outdir, exist_ok=True) os.makedirs(args.outdir, exist_ok=True)
# cluster plot # cluster plot
cmap = truncate_colormap(plt.get_cmap('jet'), 0.1, 0.9) fig = plt.figure(figsize=(15, 12), dpi=160)
fig = plt.figure(figsize=(20, 16), dpi=160)
ax = fig.add_subplot(111, projection='3d') ax = fig.add_subplot(111, projection='3d')
# draw particle line # draw particle line
ax.plot(*pline[[2, 1]], '--', zs=pline[0], color='green') ax.plot(*pline[[2, 1]], '--', zs=pline[0], color='green')
...@@ -161,17 +180,17 @@ if __name__ == '__main__': ...@@ -161,17 +180,17 @@ if __name__ == '__main__':
# zoomed-in cluster plot # zoomed-in cluster plot
fig = plt.figure(figsize=(20, 16), dpi=160) fig = plt.figure(figsize=(15, 12), dpi=160)
ax = fig.add_subplot(111, projection='3d') ax = fig.add_subplot(111, projection='3d')
# draw particle line # draw particle line
ax.plot(*pline[[2, 1]], '--', zs=pline[0], color='green') ax.plot(*pline[[2, 1]], '--', zs=pline[0], color='green')
# draw hits # draw hits
draw_hits3d(ax, df[['x', 'y', 'z', 'edep']].values, cmap, s=20.0) draw_hits3d(ax, df[['x', 'y', 'z', 'edep']].values, cmap, s=20.0)
draw_track_fit(ax, df, stop_layer=args.stop, # 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)) # scat_kw=dict(color='k', s=50.0), line_kw=dict(linestyle=':', color='k', lw=3))
# view range # view range
center = (length + thickness/2.)*vec center = (length + thickness/2.)*vec
ranges = np.vstack([center - thickness/4., center + thickness/4.]).T ranges = np.vstack([center - thickness/args.zoom_factor, center + thickness/args.zoom_factor]).T
ax.set_zlim(*ranges[0]) ax.set_zlim(*ranges[0])
ax.set_ylim(*ranges[1]) ax.set_ylim(*ranges[1])
ax.set_xlim(*ranges[2]) ax.set_xlim(*ranges[2])
...@@ -181,31 +200,27 @@ if __name__ == '__main__': ...@@ -181,31 +200,27 @@ if __name__ == '__main__':
# projection plot # projection plot
# convert to polar coordinates (mrad), and stack all r values
df['r'] = np.sqrt(df['x'].values**2 + df['y'].values**2 + df['z'].values**2)
df['phi'] = np.arctan2(df['y'].values, df['x'].values)*1000.
df['theta'] = np.arccos(df['z'].values/df['r'].values)*1000.
# convert to mrad # convert to mrad
vecp = np.asarray([np.arccos(vec[2]), np.arctan2(vec[1], vec[0])])*1000. vecp = np.asarray([np.arccos(vec[2]), np.arctan2(vec[1], vec[0])])*1000.
phi_rg = np.asarray([vecp[1] - args.topo_range, vecp[1] + args.topo_range]) phi_rg = np.asarray([vecp[1] - args.topo_range, vecp[1] + args.topo_range])
th_rg = np.asarray([vecp[0] - args.topo_range, vecp[0] + args.topo_range]) th_rg = np.asarray([vecp[0] - args.topo_range, vecp[0] + args.topo_range])
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=(17, 16), dpi=160, gridspec_kw={'wspace':0., 'width_ratios': [16, 1]}) 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['theta'].values, df['phi'].values, weights=df['edep'].values, ax, sm = draw_heatmap(axs[0], df['eta'].values, df['phi'].values, weights=df['edep'].values,
bins=(np.arange(*th_rg, step=args.topo_size), np.arange(*phi_rg, step=args.topo_size)), 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')) cmap=cmap, cmin=0., pc_kw=dict(alpha=0.8, edgecolor='k'))
ax.set_xlabel(r'$\theta$ (mrad)', fontsize=28) ax.set_ylabel(r'$\phi$ (mrad)', fontsize=32)
ax.set_ylabel(r'$\phi$ (mrad)', fontsize=28) ax.set_xlabel(r'$\eta$', fontsize=32)
ax.tick_params(labelsize=24) ax.tick_params(labelsize=28)
ax.xaxis.set_minor_locator(MultipleLocator(5)) ax.xaxis.set_minor_locator(MultipleLocator(5))
ax.yaxis.set_minor_locator(MultipleLocator(5)) ax.yaxis.set_minor_locator(MultipleLocator(5))
ax.grid(linestyle=':', which='both') ax.grid(linestyle=':', which='both')
ax.set_axisbelow(True) ax.set_axisbelow(True)
cb = plt.colorbar(sm, cax=axs[1], shrink=0.85) cb = plt.colorbar(sm, cax=axs[1], shrink=0.85, aspect=1.2*20)
cb.ax.tick_params(labelsize=24) cb.ax.tick_params(labelsize=28)
cb.ax.get_yaxis().labelpad = 24 cb.ax.get_yaxis().labelpad = 10
cb.ax.set_ylabel('Energy Deposit (MeV)', rotation=90, fontsize=28) cb.ax.set_ylabel('Energy Deposit (MeV)', rotation=90, fontsize=32)
fig.savefig(os.path.join(args.outdir, 'e{}_topo.png'.format(args.iev))) fig.savefig(os.path.join(args.outdir, 'e{}_topo.png'.format(args.iev)))
benchmarks/sampling_ecal/scripts/draw_cluster_layers.py
+ 28
− 15
View file @ 22b07fd9
...@@ -71,9 +71,8 @@ if __name__ == '__main__': ...@@ -71,9 +71,8 @@ if __name__ == '__main__':
'EcalBarrelAstroPix_RMin', 'EcalBarrelAstroPix_RMin',
'EcalBarrelAstroPix_ReadoutLayerThickness', 'EcalBarrelAstroPix_ReadoutLayerThickness',
'EcalBarrelAstroPix_ReadoutLayerNumber', 'EcalBarrelAstroPix_ReadoutLayerNumber',
'EcalBarrelAstroPix_Length' 'EcalBarrelAstroPix_Length',
]) ])
if not len(desc): if not len(desc):
# or define Ecal shapes # or define Ecal shapes
rmin, thickness, length = 890, 20*(10. + 1.65), 860*2+500 rmin, thickness, length = 890, 20*(10. + 1.65), 860*2+500
...@@ -87,27 +86,41 @@ if __name__ == '__main__': ...@@ -87,27 +86,41 @@ if __name__ == '__main__':
# read data # read data
load_root_macros(args.macros) load_root_macros(args.macros)
df = get_hits_data(args.file, args.iev, args.branch) df = get_hits_data(args.file, args.iev, args.branch)
dfmcp = get_mcp_data(args.file, args.iev, 'mcparticles2')
vec = dfmcp.loc[dfmcp['status'] == 24578, ['px', 'py', 'pz']].iloc[0].values
vec = vec/np.linalg.norm(vec)
df = df[df['cluster'] == args.icl] df = df[df['cluster'] == args.icl]
if not len(df):
print("Error: do not find any hits for cluster {:d} in event {:d}".format(args.icl, args.iev))
exit(-1)
# convert to polar coordinates (mrad), and stack all r values # convert to polar coordinates (mrad), and stack all r values
df['r'] = np.sqrt(df['x'].values**2 + df['y'].values**2 + df['z'].values**2) df['r'] = np.sqrt(df['x'].values**2 + df['y'].values**2 + df['z'].values**2)
df['phi'] = np.arctan2(df['y'].values, df['x'].values)*1000. df['phi'] = np.arctan2(df['y'].values, df['x'].values)*1000.
df['theta'] = np.arccos(df['z'].values/df['r'].values)*1000. df['theta'] = np.arccos(df['z'].values/df['r'].values)*1000.
if not len(df): df['eta'] = -np.log(np.tan(df['theta'].values/1000./2.))
print("Error: do not find any hits for cluster {:d} in event {:d}".format(args.icl, args.iev))
exit(-1) # truth
dfmcp = get_mcp_simple(args.file, args.iev, 'mcparticles2').iloc[0]
pdgbase = ROOT.TDatabasePDG()
inpart = pdgbase.GetParticle(int(dfmcp['pid']))
print("Incoming particle = {}, pdgcode = {}, charge = {}, mass = {}"\
.format(inpart.GetName(), inpart.PdgCode(), inpart.Charge(), inpart.Mass()))
# 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
else:
flayer = df[df['layer'] == df['layer'].min()]
vec = flayer[['x', 'y', 'z']].mean().values
vec = vec/np.linalg.norm(vec)
# particle line from (0, 0, 0) to the inner Ecal surface # particle line from (0, 0, 0) to the inner Ecal surface
length = rmin/np.sqrt(vec[0]**2 + vec[1]**2) length = rmin/np.sqrt(vec[0]**2 + vec[1]**2)
pline = np.transpose(vec*np.mgrid[0:length:2j][:, np.newaxis]) pline = np.transpose(vec*np.mgrid[0:length:2j][:, np.newaxis])
cmap = truncate_colormap(plt.get_cmap('jet'), 0.1, 0.9) cmap = truncate_colormap(plt.get_cmap('jet'), 0.1, 0.9)
# convert truth to mrad # convert truth to mrad
vecp = np.asarray([np.arccos(vec[2]), np.arctan2(vec[1], vec[0])])*1000. vecp = np.asarray([np.arccos(vec[2]), np.arctan2(vec[1], vec[0])])*1000.
phi_rg = np.asarray([vecp[1] - args.topo_range, vecp[1] + args.topo_range]) phi_rg = np.asarray([vecp[1] - args.topo_range, vecp[1] + args.topo_range])
th_rg = np.asarray([vecp[0] - args.topo_range, vecp[0] + args.topo_range]) th_rg = np.asarray([vecp[0] - args.topo_range, vecp[0] + args.topo_range])
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) os.makedirs(args.outdir, exist_ok=True)
# cluster plot by layers (rebinned) # cluster plot by layers (rebinned)
...@@ -117,11 +130,11 @@ if __name__ == '__main__': ...@@ -117,11 +130,11 @@ if __name__ == '__main__':
for i in np.arange(1, df['layer'].max() + 1, dtype=int): for i in np.arange(1, df['layer'].max() + 1, dtype=int):
data = df[df['layer'] == i] data = df[df['layer'] == i]
fig, axs = plt.subplots(1, 2, figsize=(17, 16), dpi=160, gridspec_kw={'wspace':0., 'width_ratios': [16, 1]}) 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['theta'].values, data['phi'].values, weights=data['edep'].values, ax, sm = draw_heatmap(axs[0], data['eta'].values, data['phi'].values, weights=data['edep'].values,
bins=(np.arange(*th_rg, step=args.topo_size), np.arange(*phi_rg, step=args.topo_size)), 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')) cmap=cmap, vmin=0.0, vmax=1.0, pc_kw=dict(alpha=0.8, edgecolor='k'))
ax.set_ylabel(r'$\phi$ (mrad)', fontsize=28) ax.set_ylabel(r'$\phi$ (mrad)', fontsize=28)
ax.set_xlabel(r'$\theta$ (mrad)', fontsize=28) ax.set_xlabel(r'$\eta$', fontsize=28)
ax.tick_params(labelsize=24) ax.tick_params(labelsize=24)
ax.xaxis.set_minor_locator(MultipleLocator(5)) ax.xaxis.set_minor_locator(MultipleLocator(5))
ax.yaxis.set_minor_locator(MultipleLocator(5)) ax.yaxis.set_minor_locator(MultipleLocator(5))
...@@ -155,11 +168,11 @@ if __name__ == '__main__': ...@@ -155,11 +168,11 @@ if __name__ == '__main__':
fig, axs = plt.subplots(1, 2, figsize=(17, 16), dpi=160, gridspec_kw={'wspace':0., 'width_ratios': [16, 1]}) fig, axs = plt.subplots(1, 2, figsize=(17, 16), dpi=160, gridspec_kw={'wspace':0., 'width_ratios': [16, 1]})
ax = axs[0] ax = axs[0]
colors = cmap(data['edep']/1.0) colors = cmap(data['edep']/1.0)
ax.scatter(data['theta'].values, data['phi'].values, c=colors, marker='s', s=15.0) ax.scatter(data['eta'].values, data['phi'].values, c=colors, marker='s', s=15.0)
ax.set_xlim(*th_rg) ax.set_xlim(*eta_rg)
ax.set_ylim(*phi_rg) ax.set_ylim(*phi_rg)
ax.set_ylabel(r'$\phi$ (mrad)', fontsize=28) ax.set_ylabel(r'$\phi$ (mrad)', fontsize=28)
ax.set_xlabel(r'$\theta$ (mrad)', fontsize=28) ax.set_xlabel(r'$\eta$', fontsize=28)
ax.tick_params(labelsize=24) ax.tick_params(labelsize=24)
ax.xaxis.set_minor_locator(MultipleLocator(5)) ax.xaxis.set_minor_locator(MultipleLocator(5))
ax.yaxis.set_minor_locator(MultipleLocator(5)) ax.yaxis.set_minor_locator(MultipleLocator(5))
......
benchmarks/sampling_ecal/scripts/utils.py
+ 51
− 1
View file @ 22b07fd9
...@@ -41,7 +41,7 @@ def get_mcp_data(path, evnums=None, branch='mcparticles2'): ...@@ -41,7 +41,7 @@ def get_mcp_data(path, evnums=None, branch='mcparticles2'):
elif isinstance(evnums, int): elif isinstance(evnums, int):
evnums = [evnums] evnums = [evnums]
dbuf = np.zeros(shape=(2000*len(evnums), 6)) dbuf = np.zeros(shape=(5000*len(evnums), 6))
idb = 0 idb = 0
for iev in evnums: for iev in evnums:
if iev >= events.GetEntries(): if iev >= events.GetEntries():
...@@ -122,3 +122,53 @@ def compact_constants(path, names): ...@@ -122,3 +122,53 @@ def compact_constants(path, names):
kernel.terminate() kernel.terminate()
return vals return vals
# read clusters data from root file
def get_clusters_data(path, evnums=None, branch='EcalBarrelClustersReco'):
f = ROOT.TFile(path)
events = f.events
if evnums is None:
evnums = np.arange(events.GetEntries())
elif isinstance(evnums, int):
evnums = [evnums]
dbuf = np.zeros(shape=(20*len(evnums), 6))
idb = 0
for iev in evnums:
if iev >= events.GetEntries():
print('Error: event {:d} is out of range (0 - {:d})'.format(iev, events.GetEntries() - 1))
continue
events.GetEntry(iev)
for k, cl in enumerate(getattr(events, branch)):
dbuf[idb] = (iev, k, cl.nhits, cl.edep, cl.cl_theta, cl.cl_phi)
idb += 1
return pd.DataFrame(data=dbuf[:idb], columns=['event', 'cluster', 'nhits', 'edep', 'cl_theta', 'cl_phi'])
# read mc particles from root file
def get_mcp_simple(path, evnums=None, branch='mcparticles2'):
f = ROOT.TFile(path)
events = f.events
if evnums is None:
evnums = np.arange(events.GetEntries())
elif isinstance(evnums, int):
evnums = [evnums]
dbuf = np.zeros(shape=(len(evnums), 6))
idb = 0
for iev in evnums:
if iev >= events.GetEntries():
print('Error: event {:d} is out of range (0 - {:d})'.format(iev, events.GetEntries() - 1))
continue
events.GetEntry(iev)
# extract full mc particle data
part = getattr(events, branch)[2]
dbuf[idb] = (iev, part.psx, part.psy, part.psz, part.pdgID, part.status)
idb += 1
return pd.DataFrame(data=dbuf[:idb], columns=['event', 'px', 'py', 'pz', 'pid', 'status'])
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment