gwm17/spspy
1
0
Fork 0
mirror of https://github.com/gwm17/spspy.git synced 2024-11-22 10:18:49 -05:00
Code Issues Projects Releases Wiki Activity

Merge pull request #3 from gwm17/main

Browse Source 
Merge hotfixes to devel
This commit is contained in:
Gordon McCann 2023-02-03 11:30:04 -05:00 committed by GitHub
parent 2562156943 3971b797d6
commit eb092d6c4b
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
5 changed files with 20 additions and 19 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
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:

4
spspy/SPSPlotUI.py
View File

@ -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()

20
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))
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

9
spspy/SPSTarget.py
View File

@ -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

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
class NuclearDataMap:
U2MEV: float = 931.4940954
U2MEV: float = 931.493614838475
ELECTRON_MASS: float = 0.000548579909
def __init__(self):