From a81eb106b10dad089fe4ab5f173e696d3d0a720e Mon Sep 17 00:00:00 2001
From: Carl Gwilliam <gwilliam@hep.ph.liv.ac.uk>
Date: Thu, 28 Apr 2022 17:42:49 +0100
Subject: [PATCH] Add energy cut into Fluka reader
---
.../FlukaReader/python/FlukaReaderAlg.py | 27 +++++++++++++++----
1 file changed, 22 insertions(+), 5 deletions(-)
diff --git a/Generators/FlukaReader/python/FlukaReaderAlg.py b/Generators/FlukaReader/python/FlukaReaderAlg.py
index b54aa07a..3a8728e3 100644
--- a/Generators/FlukaReader/python/FlukaReaderAlg.py
+++ b/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())
--
GitLab