gwm17/spspy
1
0
Fork 0
mirror of https://github.com/gwm17/spspy.git synced 2024-11-26 12:08:51 -05:00
Code Issues Projects Releases Wiki Activity

Fix a few potential bugs in reaction by better separating UI parameters from calculation parameters.

Browse Source
This commit is contained in:
Gordon McCann 2023-02-03 10:33:58 -05:00
parent 4f622097cd
commit 3971b797d6
4 changed files with 14 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
spspy/SPSPlot.py
View File

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

4
spspy/SPSPlotUI.py
View File

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

18
spspy/SPSReaction.py
View File

@ -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)) 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: class Reaction:
DEG2RAD: float = pi/180.0 #degrees -> radians
C = 299792458 #speed of light m/s C = 299792458 #speed of light m/s
QBRHO2P = 1.0E-9*C #Converts qbrho to momentum (p) (kG*cm -> MeV/c) QBRHO2P = 1.0E-9*C #Converts qbrho to momentum (p) (kG*cm -> MeV/c)
FP_MAGNIFICATION = 0.39 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}" return f"{self.params.target.prettyIsotopicSymbol}({self.params.projectile.prettyIsotopicSymbol},{self.params.ejectile.prettyIsotopicSymbol}){self.residual.prettyIsotopicSymbol}"
def __repr__(self) -> str: 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: 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()}" 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,13 +51,12 @@ class Reaction:
#MeV #MeV
def calculate_ejectile_KE(self, excitation: float) -> float: def calculate_ejectile_KE(self, excitation: float) -> float:
rxnQ = self.Qvalue - excitation 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) 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) threshold = -rxnQ*(self.params.ejectile.mass+self.residual.mass)/(self.params.ejectile.mass + self.residual.mass - self.params.projectile.mass)
if beamRxnEnergy < threshold: if beamRxnEnergy < threshold:
return INVALID_KINETIC_ENERGY return INVALID_KINETIC_ENERGY
term1 = sqrt(self.params.projectile.mass * self.params.ejectile.mass * beamRxnEnergy) / (self.params.ejectile.mass + self.residual.mass) * cos(angleRads) 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) 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: if(term1**2.0 + term2) < 0:
return INVALID_KINETIC_ENERGY return INVALID_KINETIC_ENERGY
@ -72,7 +70,7 @@ class Reaction:
else: else:
ejectileEnergy = ke2**2.0 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 return ejectileEnergy
def convert_ejectile_KE_2_rho(self, ejectileEnergy: float) -> float: def convert_ejectile_KE_2_rho(self, ejectileEnergy: float) -> float:
@ -84,27 +82,27 @@ class Reaction:
return qbrho / (float(self.params.ejectile.Z) * self.params.magneticField) return qbrho / (float(self.params.ejectile.Z) * self.params.magneticField)
def calculate_excitation(self, rho: float) -> float: 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 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 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)) 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) 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)) 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 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 return sqrt(residRxnEnergy**2.0 - residRxnP2) - self.residual.mass
def calculate_focal_plane_offset(self, ejectileEnergy: float) -> float: def calculate_focal_plane_offset(self, ejectileEnergy: float) -> float:
if ejectileEnergy == INVALID_KINETIC_ENERGY: if ejectileEnergy == INVALID_KINETIC_ENERGY:
return 0.0 return 0.0
ejectileRho = self.convert_ejectile_KE_2_rho(ejectileEnergy) 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 = 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 * 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)
return -1.0*k*ejectileRho*self.FP_DISPERSION*self.FP_MAGNIFICATION return -1.0*k*ejectileRho*self.FP_DISPERSION*self.FP_MAGNIFICATION
#(MeV, rad, kG)
def update_parameters(self, beamEnergy: float, spsAngle: float, magenticField: float): def update_parameters(self, beamEnergy: float, spsAngle: float, magenticField: float):
self.params.beamEnergy = beamEnergy self.params.beamEnergy = beamEnergy
self.params.spsAngle = spsAngle self.params.spsAngle = spsAngle

2
spspy/data/NuclearData.py
View File

@ -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 return z*z + z + a if z > a else a*a + z
class NuclearDataMap: class NuclearDataMap:
U2MEV: float = 931.4940954 U2MEV: float = 931.493614838475
ELECTRON_MASS: float = 0.000548579909 ELECTRON_MASS: float = 0.000548579909
def __init__(self): def __init__(self):