Fix bug in SPSTarget where catima projectile was given (z, a) when should have been given (a, z)

spspy/SPSReaction.py

@@ -58,8 +58,10 @@ class Reaction:
         if beamRxnEnergy < threshold:
             return INVALID_KINETIC_ENERGY
         
-        term1 = sqrt(self.params.projectile.mass * self.params.ejectile.mass * beamRxnEnergy) / (self.params.ejectile.mass + self.residual.mass) * cos(angleRads * self.DEG2RAD)
+        term1 = sqrt(self.params.projectile.mass * self.params.ejectile.mass * beamRxnEnergy) / (self.params.ejectile.mass + self.residual.mass) * cos(angleRads)
         term2 = (beamRxnEnergy * (self.residual.mass - self.params.projectile.mass) + self.residual.mass * rxnQ) / (self.params.ejectile.mass + self.residual.mass)
+        if(term1**2.0 + term2) < 0:
+            return INVALID_KINETIC_ENERGY
 
         ke1 = term1 + sqrt(term1**2.0 + term2)
         ke2 = term1 + sqrt(term1**2.0 + term2)

spspy/SPSTarget.py

@@ -104,7 +104,7 @@ class SPSTarget:
 
     def get_rxn_layer(self, zt: int, at: int) -> int:
         for idx, layer in enumerate(self.layer_details):
-            for (a, z, s) in layer.compound_list:
+            for (z, a, s) in layer.compound_list:
                 if at == a and zt == z:
                     return idx
         return INVALID_RXN_LAYER
@@ -114,19 +114,21 @@ class SPSTarget:
         if angle == pi*0.5:
             return e_initial
 
-        projectile = catima.Projectile(zp, ap)
+        print(f"z{zp} a{ap}")
+        projectile = catima.Projectile(ap, zp)
         e_current = e_initial/ap
-
         for (idx, layer) in enumerate(self.layer_details):
             material = catima.Material([(global_nuclear_data.get_data(z, a).mass, z, float(s)) for (z, a, s) in layer.compound_list])
             projectile.T(e_current) #catima wants MeV/u
             if idx == rxn_layer:
                 material.thickness(self.layer_details[idx].thickness * self.UG2G / (2.0 * abs(cos(angle))))
                 e_current -= get_energyloss(projectile, material)
+                print("e_current: ", e_current*ap)
                 return e_initial - e_current*ap
             else:
                 material.thickness(self.layer_details[idx].thickness * self.UG2G / abs(cos(angle)))
                 e_current -= get_energyloss(projectile, material)
+                print("e_current2: ", e_current*ap)
 
         return e_initial - e_current*ap
 
@@ -135,7 +137,7 @@ class SPSTarget:
         if angle == pi*0.5:
             return e_initial
 
-        projectile = catima.Projectile(zp, ap)
+        projectile = catima.Projectile(ap, zp)
         e_current = e_initial/ap
 
         for (idx, layer) in enumerate(self.layer_details[rxn_layer:], start=rxn_layer):
@@ -154,7 +156,7 @@ class SPSTarget:
         if angle == pi*0.5:
             return 0.0
 
-        projectile = catima.Projectile(zp, ap)
+        projectile = catima.Projectile(ap, zp)
         e_current = e_final/ap
         sublist = self.layer_details[rxn_layer:] #only care about rxn_layer -> exit
         reveresedRxnLayer = len(sublist) -1 #when reversed rxn_layer is the last layer