Merge pull request #3 from gwm17/main

Merge hotfixes to devel

spspy/SPSPlot.py

@@ -4,6 +4,8 @@ from .data.NuclearData import *
 from dataclasses import dataclass, field
 import csv
 
+DEG2RAD: float = np.pi / 180.0
+
 @dataclass
 class Excitation:
     excitation: float = 0.0 #MeV
@@ -45,7 +47,7 @@ class SPSPlot:
         
     def update_reactions(self) -> None:
         for datum in self.data.values():
-            datum.rxn.update_parameters(self.beamEnergy, self.spsAngle, self.magneticField)
+            datum.rxn.update_parameters(self.beamEnergy, self.spsAngle * DEG2RAD, self.magneticField)
             if datum.rxn.targetMaterial.name in self.targets:
                 datum.rxn.targetMaterial = self.targets[datum.rxn.targetMaterial.name]
             for ex in datum.excitations:

spspy/SPSPlotUI.py

@@ -1,4 +1,4 @@
-from .SPSPlot import SPSPlot
+from .SPSPlot import SPSPlot, DEG2RAD
 from .SPSReaction import RxnParameters
 from .ui.MPLCanvas import MPLCanvas
 from .ui.ReactionDialog import ReactionDialog
@@ -236,7 +236,7 @@ class SPSPlotGUI(QMainWindow):
 
     def add_reaction(self, rxnParams: RxnParameters, targName: str) -> None:
         rxnParams.beamEnergy = self.bkeInput.value()
-        rxnParams.spsAngle = self.angleInput.value()
+        rxnParams.spsAngle = self.angleInput.value() * DEG2RAD
         rxnParams.magneticField = self.bfieldInput.value()
         self.sps.add_reaction(rxnParams, targName)
         self.update_plot()

spspy/SPSReaction.py

@@ -18,7 +18,6 @@ def create_reaction_parameters(zt: int, at: int, zp: int, ap: int, ze: int, ae:
     return RxnParameters(global_nuclear_data.get_data(zt, at), global_nuclear_data.get_data(zp, ap), global_nuclear_data.get_data(ze, ae))
 
 class Reaction:
-    DEG2RAD: float = pi/180.0 #degrees -> radians
     C = 299792458 #speed of light m/s
     QBRHO2P = 1.0E-9*C #Converts qbrho to momentum (p) (kG*cm -> MeV/c)
     FP_MAGNIFICATION = 0.39
@@ -44,7 +43,7 @@ class Reaction:
         return f"{self.params.target.prettyIsotopicSymbol}({self.params.projectile.prettyIsotopicSymbol},{self.params.ejectile.prettyIsotopicSymbol}){self.residual.prettyIsotopicSymbol}"
 
     def __repr__(self) -> str:
-        return f"{self.params.target.isotopicSymbol}({self.params.projectile.isotopicSymbol},{self.params.ejectile.isotopicSymbol}){self.residual.isotopicSymbol}_{self.params.beamEnergy}MeV_{self.params.spsAngle}deg_{self.params.magneticField}kG"
+        return f"{self.params.target.isotopicSymbol}({self.params.projectile.isotopicSymbol},{self.params.ejectile.isotopicSymbol}){self.residual.isotopicSymbol}_{self.params.beamEnergy}MeV_{self.params.spsAngle}rad_{self.params.magneticField}kG"
 
     def get_latex_rep(self) -> str:
         return f"{self.params.target.get_latex_rep()}({self.params.projectile.get_latex_rep()},{self.params.ejectile.get_latex_rep()}){self.residual.get_latex_rep()}"
@@ -52,14 +51,15 @@ class Reaction:
     #MeV
     def calculate_ejectile_KE(self, excitation: float) -> float:
         rxnQ = self.Qvalue - excitation
-        angleRads = self.params.spsAngle * self.DEG2RAD
         beamRxnEnergy = self.params.beamEnergy - self.targetMaterial.get_incoming_energyloss(self.params.projectile.Z, self.params.projectile.mass, self.params.beamEnergy, self.rxnLayer, 0.0)
         threshold = -rxnQ*(self.params.ejectile.mass+self.residual.mass)/(self.params.ejectile.mass + self.residual.mass - self.params.projectile.mass)
         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(self.params.spsAngle)
         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)
@@ -70,7 +70,7 @@ class Reaction:
         else:
             ejectileEnergy = ke2**2.0
 
-        ejectileEnergy -= self.targetMaterial.get_outgoing_energyloss(self.params.ejectile.Z, self.params.ejectile.mass, ejectileEnergy, self.rxnLayer, angleRads)
+        ejectileEnergy -= self.targetMaterial.get_outgoing_energyloss(self.params.ejectile.Z, self.params.ejectile.mass, ejectileEnergy, self.rxnLayer, self.params.spsAngle)
         return ejectileEnergy
 
     def convert_ejectile_KE_2_rho(self, ejectileEnergy: float) -> float:
@@ -82,27 +82,27 @@ class Reaction:
         return qbrho / (float(self.params.ejectile.Z) * self.params.magneticField)
 
     def calculate_excitation(self, rho: float) -> float:
-        angleRads = self.params.spsAngle * self.DEG2RAD
         ejectileP = rho * float(self.params.ejectile.Z) * self.params.magneticField * self.QBRHO2P
         ejectileEnergy  = sqrt(ejectileP**2.0 + self.params.ejectile.mass**2.0) - self.params.ejectile.mass
-        ejectileRxnEnergy = ejectileEnergy +  self.targetMaterial.get_outgoing_reverse_energyloss(self.params.ejectile.Z, self.params.ejectile.mass, ejectileEnergy, self.rxnLayer, angleRads)
+        ejectileRxnEnergy = ejectileEnergy +  self.targetMaterial.get_outgoing_reverse_energyloss(self.params.ejectile.Z, self.params.ejectile.mass, ejectileEnergy, self.rxnLayer, self.params.spsAngle)
         ejectileRxnP = sqrt(ejectileRxnEnergy * (ejectileRxnEnergy + 2.0 * self.params.ejectile.mass))
         beamRxnEnergy = self.params.beamEnergy - self.targetMaterial.get_incoming_energyloss(self.params.projectile.Z, self.params.projectile.mass, self.params.beamEnergy, self.rxnLayer, 0.0)
         beamRxnP = sqrt(beamRxnEnergy * (beamRxnEnergy + 2.0 * self.params.projectile.mass))
 
 
         residRxnEnergy = beamRxnEnergy + self.params.projectile.mass + self.params.target.mass - ejectileRxnEnergy - self.params.ejectile.mass
-        residRxnP2 = beamRxnP**2.0 + ejectileRxnP**2.0 - 2.0 * ejectileRxnP * beamRxnP * cos(angleRads)
+        residRxnP2 = beamRxnP**2.0 + ejectileRxnP**2.0 - 2.0 * ejectileRxnP * beamRxnP * cos(self.params.spsAngle)
         return sqrt(residRxnEnergy**2.0 - residRxnP2) - self.residual.mass
 
     def calculate_focal_plane_offset(self, ejectileEnergy: float) -> float:
         if ejectileEnergy == INVALID_KINETIC_ENERGY:
             return 0.0
         ejectileRho = self.convert_ejectile_KE_2_rho(ejectileEnergy)
-        k = sqrt(self.params.projectile.mass * self.params.ejectile.mass * self.params.beamEnergy / ejectileEnergy) * sin(self.params.spsAngle * self.DEG2RAD)
-        k /= self.params.ejectile.mass + self.residual.mass - sqrt(self.params.projectile.mass * self.params.ejectile.mass * self.params.beamEnergy/ejectileEnergy) * cos(self.params.spsAngle * self.DEG2RAD)
+        k = sqrt(self.params.projectile.mass * self.params.ejectile.mass * self.params.beamEnergy / ejectileEnergy) * sin(self.params.spsAngle)
+        k /= self.params.ejectile.mass + self.residual.mass - sqrt(self.params.projectile.mass * self.params.ejectile.mass * self.params.beamEnergy/ejectileEnergy) * cos(self.params.spsAngle)
         return -1.0*k*ejectileRho*self.FP_DISPERSION*self.FP_MAGNIFICATION
 
+    #(MeV, rad, kG)
     def update_parameters(self, beamEnergy: float, spsAngle: float, magenticField: float):
         self.params.beamEnergy = beamEnergy
         self.params.spsAngle = spsAngle

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,9 +114,8 @@ 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):
             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
@@ -135,7 +134,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 +153,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

spspy/data/NuclearData.py

@@ -24,7 +24,7 @@ def generate_nucleus_id(z: np.uint32, a: np.uint32) -> np.uint32 :
     return z*z + z + a if z > a else a*a + z
 
 class NuclearDataMap:
-    U2MEV: float = 931.4940954
+    U2MEV: float = 931.493614838475
     ELECTRON_MASS: float = 0.000548579909
 
     def __init__(self):