Add energy cut into Fluka reader

Generators/FlukaReader/python/FlukaReaderAlg.py

@@ -36,7 +36,8 @@ class FlukaReader(EvgenAlg):
 
         if self.test:
             self.before = dict(zip(self.columns, [[] for i in range(len(self.columns))]))
-            self.after = dict(zip(self.columns, [[] for i in range(len(self.columns))]))                           
+            self.after = dict(zip(self.columns, [[] for i in range(len(self.columns))]))
+
                
         return
 
@@ -153,6 +154,8 @@ class FlukaReader(EvgenAlg):
         thetaX = self.angle(cosThetaX)
         
         xout = np.copy(x)
+        if xout.ndim == 0:
+            xout = float(xout)
         
         # Add displacement due to initial angle
         #xang = z * np.tan(thetaX)
@@ -275,7 +278,16 @@ class FlukaReader(EvgenAlg):
         gp = HepMC.GenParticle()
 
         m = 105.66
-        e = newentry["E"] * 1000.  #MeV      
+        e = newentry["E"] * 1000.  #MeV
+
+        # If the energy is less than mass then skip the event
+        if e < m:
+            self.setFilterPassed(False)
+            self.msg.debug("Event failed energy cut")
+            return False
+        else:
+            self.setFilterPassed(True)
+        
         p = np.sqrt(e**2 - m**2)
 
         thetaX = self.angle(newentry["cosX"])        
@@ -288,22 +300,27 @@ class FlukaReader(EvgenAlg):
         #phi = np.arctan(newentry["cosY"] / newentry["cosX"])
         if phi < 0: phi += 2*np.pi
         if phi == 2*np.pi: phi = 0
+
+        #self.msg.debug(f"INPUT: {e}, {m}, {p}, {theta}, {phi}, {np.sin(theta)}, {np.cos(theta)}, {np.sin(phi)}, {np.cos(phi)}")
                     
         px = p * np.sin(theta) * np.cos(phi)
         py = p * np.sin(theta) * np.sin(phi)        
         pz = p * np.cos(theta) 
-        
+
         mom = HepMC.FourVector(px, py, pz, e)
 
         gp.set_momentum(mom)
         gp.set_generated_mass(m)
         gp.set_pdg_id(self.pid(newentry["type"]))
         gp.set_status(1)
+
+        #self.msg.debug(f"HEPMC:{px, py, pz, e}")
+        #gp.print()
         
         ROOT.SetOwnership(gp, False)
         gv.add_particle_out(gp)
 
-        return
+        return True
 
     def plot(self):
         "Plot entries before and after propagation/smeating for tests"
@@ -434,6 +451,6 @@ if __name__ == "__main__":
     itemList = [ "EventInfo#McEventInfo", "McEventCollection#*" ]
     from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
     cfg.merge(OutputStreamCfg(ConfigFlags, "EVNT", itemList, disableEventTag = True))
-
+    cfg.getEventAlgo("OutputStreamEVNT").AcceptAlgs = ["FlukaReader"]
     sc = cfg.run(maxEvents = getNEvents(args.file, args.nevents) * args.nsamples)
     sys.exit(not sc.isSuccess())