diff --git a/benchmarks/imaging_ecal/scripts/draw_cluster.py b/benchmarks/imaging_ecal/scripts/draw_cluster.py
index 69b641ecb7fee402e10c2a9b9036dc759a82cdf2..35dbf0c97deeb0490540edf613f1ac6d019cd8f9 100644
--- a/benchmarks/imaging_ecal/scripts/draw_cluster.py
+++ b/benchmarks/imaging_ecal/scripts/draw_cluster.py
@@ -151,17 +151,18 @@ if __name__ == '__main__':
     dfallmcp = get_all_mcp(args.file, args.iev, 'mcparticles2')
     # Select decaying particles
     dftemp = dfallmcp[dfallmcp['g4Parent'] == 1.0]
-    dfdecaymcp = dftemp.copy()
-    pdgbase = ROOT.TDatabasePDG()
-    for iptl in [0, len(dfdecaymcp) - 1]:
-        infoptl = pdgbase.GetParticle(int(dfdecaymcp['pid'].iloc[iptl]))
-        print("{} Decaying particle = {}, pdgcode = {}, charge = {}, mass = {}"\
-            .format(iptl, infoptl.GetName(), infoptl.PdgCode(), infoptl.Charge(), infoptl.Mass()))
-    # Calculate geometric variables of decaying particles
-    dfdecaymcp['r'] = np.sqrt(dfdecaymcp['vex'].values**2 + dfdecaymcp['vey'].values**2 + dfdecaymcp['vez'].values**2)
-    dfdecaymcp['phi'] = np.arctan2(dfdecaymcp['vey'].values, dfdecaymcp['vex'].values)*1000.
-    dfdecaymcp['theta'] = np.arccos(dfdecaymcp['vez'].values/dfdecaymcp['r'].values)*1000.
-    dfdecaymcp['eta'] = -np.log(np.tan(dfdecaymcp['theta'].values/1000./2.))
+    if len(dfdecaymcp) > 0:
+        dfdecaymcp = dftemp.copy()
+        pdgbase = ROOT.TDatabasePDG()
+        for iptl in [0, len(dfdecaymcp) - 1]:
+            infoptl = pdgbase.GetParticle(int(dfdecaymcp['pid'].iloc[iptl]))
+            print("{} Decaying particle = {}, pdgcode = {}, charge = {}, mass = {}"\
+                .format(iptl, infoptl.GetName(), infoptl.PdgCode(), infoptl.Charge(), infoptl.Mass()))
+        # Calculate geometric variables of decaying particles
+        dfdecaymcp['r'] = np.sqrt(dfdecaymcp['vex'].values**2 + dfdecaymcp['vey'].values**2 + dfdecaymcp['vez'].values**2)
+        dfdecaymcp['phi'] = np.arctan2(dfdecaymcp['vey'].values, dfdecaymcp['vex'].values)*1000.
+        dfdecaymcp['theta'] = np.arccos(dfdecaymcp['vez'].values/dfdecaymcp['r'].values)*1000.
+        dfdecaymcp['eta'] = -np.log(np.tan(dfdecaymcp['theta'].values/1000./2.))
 
     # truth
     dfmcp = get_mcp_simple(args.file, args.iev, 'mcparticles2').iloc[0]