From 5c78d656946d2bcc93b5bf6dd4fb26bcec5615bb Mon Sep 17 00:00:00 2001
From: Carl Gwilliam <gwilliam@hep.ph.liv.ac.uk>
Date: Fri, 8 Apr 2022 16:47:49 +0100
Subject: [PATCH] Fix incorrect conversion in fluka E loss

---
 Generators/FlukaReader/python/FlukaReaderAlg.py | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/Generators/FlukaReader/python/FlukaReaderAlg.py b/Generators/FlukaReader/python/FlukaReaderAlg.py
index bee9711a..fb49faba 100644
--- a/Generators/FlukaReader/python/FlukaReaderAlg.py
+++ b/Generators/FlukaReader/python/FlukaReaderAlg.py
@@ -11,6 +11,8 @@ import ROOT
 import numpy as np
 import math
 
+# TODO: correct angle for beam crossing angle: both in angel itself and position
+
 class FlukaReader(EvgenAlg):
     def __init__(self, name="FlukaReader", MCEventKey="BeamTruthEvent", file_name = "", dist = 0, randomSeed = None, nsamples = 1, test = True):
         super(FlukaReader,self).__init__(name=name)
@@ -106,7 +108,7 @@ class FlukaReader(EvgenAlg):
 
     def energy_after_loss_exact(self, e, zcorr):
         "Calculate exact energy after loss in material"
-        return Range.muPropagate(e, zcorr/1000.)
+        return Range.muPropagate(e, zcorr/100.)
 
     def energy_after_loss(self, e, cosThetaX, cosThetaY, zcorr, a = 2e-3, b = 4e-6):
         "Calculate approximate energy after loss in material"
@@ -263,11 +265,13 @@ class FlukaReader(EvgenAlg):
         ROOT.SetOwnership(gv, False)
         evt.add_vertex(gv)
 
+        # TODO: skip event if below a certain energy
+
         # Create HepMC particle
         gp = HepMC.GenParticle()
 
         m = 105.66
-        e = entry["E"] * 1000.
+        e = entry["E"] * 1000.        
         p = np.sqrt(e**2 - m**2)
 
         thetaX = self.angle(entry["cosX"])        
@@ -312,7 +316,7 @@ class FlukaReader(EvgenAlg):
         plt.hist(self.after["E"], bins = ebins, histtype='step', color = "r", fill = False, label = "after")
         plt.gca().set_yscale('log')
         plt.legend()
-        plt.savefig("energy.ong")
+        plt.savefig("energy.png")
 
         plt.figure()
         plt.xlabel("Angle to beam in X dir")        
-- 
GitLab