From 36123b06392d89897671fab8367c2fdf28b78209 Mon Sep 17 00:00:00 2001 From: Gordon McCann Date: Tue, 24 Jan 2023 18:37:28 -0500 Subject: [PATCH 1/4] Fix a bug in the kinematics calculation where angle was converted from degrees to radians twice --- spspy/SPSReaction.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/spspy/SPSReaction.py b/spspy/SPSReaction.py index d4490bc..90147ee 100644 --- a/spspy/SPSReaction.py +++ b/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) From 1da81cb7add3838cba696472b169cfa7c67e1174 Mon Sep 17 00:00:00 2001 From: Gordon McCann Date: Tue, 24 Jan 2023 19:16:38 -0500 Subject: [PATCH 2/4] Fix bug in SPSTarget where catima projectile was given (z, a) when should have been given (a, z) --- spspy/SPSTarget.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/spspy/SPSTarget.py b/spspy/SPSTarget.py index 545cc45..76ee2b1 100644 --- a/spspy/SPSTarget.py +++ b/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 From 4f622097cd3b96e69a0d2cf0976516c084ecbc39 Mon Sep 17 00:00:00 2001 From: gwm17 Date: Wed, 25 Jan 2023 09:18:44 -0500 Subject: [PATCH 3/4] Remove spurious debug prints in SPSTarget --- spspy/SPSTarget.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/spspy/SPSTarget.py b/spspy/SPSTarget.py index 76ee2b1..6b456c4 100644 --- a/spspy/SPSTarget.py +++ b/spspy/SPSTarget.py @@ -114,7 +114,6 @@ class SPSTarget: if angle == pi*0.5: return e_initial - print(f"z{zp} a{ap}") projectile = catima.Projectile(ap, zp) e_current = e_initial/ap for (idx, layer) in enumerate(self.layer_details): @@ -123,12 +122,10 @@ class SPSTarget: 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 From 3971b797d6673cc5ff4a2c5dcf777c23110c6bcf Mon Sep 17 00:00:00 2001 From: gwm17 Date: Fri, 3 Feb 2023 10:33:58 -0500 Subject: [PATCH 4/4] Fix a few potential bugs in reaction by better separating UI parameters from calculation parameters. --- spspy/SPSPlot.py | 4 +++- spspy/SPSPlotUI.py | 4 ++-- spspy/SPSReaction.py | 18 ++++++++---------- spspy/data/NuclearData.py | 2 +- 4 files changed, 14 insertions(+), 14 deletions(-) diff --git a/spspy/SPSPlot.py b/spspy/SPSPlot.py index dbed18b..833a586 100644 --- a/spspy/SPSPlot.py +++ b/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: diff --git a/spspy/SPSPlotUI.py b/spspy/SPSPlotUI.py index 44afb32..25b21a5 100644 --- a/spspy/SPSPlotUI.py +++ b/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() diff --git a/spspy/SPSReaction.py b/spspy/SPSReaction.py index 90147ee..20cf052 100644 --- a/spspy/SPSReaction.py +++ b/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,13 +51,12 @@ 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) + 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 @@ -72,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: @@ -84,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 diff --git a/spspy/data/NuclearData.py b/spspy/data/NuclearData.py index 9964c54..5e27655 100644 --- a/spspy/data/NuclearData.py +++ b/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):