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

Fixed two bugs related to energy loss units in SPSTarget. Result is change in final energies by around 20 keV>

Browse Source
This commit is contained in:
Gordon McCann 2023-03-31 16:26:43 -04:00
parent 1d6af4966a
commit b0a4cc3609
2 changed files with 21 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
spspy/SPSReaction.py
View File

@ -89,7 +89,6 @@ class Reaction:
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(self.params.spsAngle) 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

38
spspy/SPSTarget.py
View File

@ -43,10 +43,10 @@ def get_reverse_energyloss(projectile: catima.Projectile, material: catima.Mater
e_step = catima.dedx(projectile, material)*x_step e_step = catima.dedx(projectile, material)*x_step
e_initial += e_step*A_recip e_initial += e_step*A_recip
projectile.T(e_initial) projectile.T(e_initial)
return (e_initial - e_out)*projectile.A() return (e_initial - e_out)
elif depth == ADAPTIVE_DEPTH_MAX: elif depth == ADAPTIVE_DEPTH_MAX:
return e_out*projectile.A() return e_out
else: else:
e_step = catima.dedx(projectile, material)*x_step e_step = catima.dedx(projectile, material)*x_step
e_initial += e_step*A_recip e_initial += e_step*A_recip
@ -80,10 +80,10 @@ def get_energyloss(projectile: catima.Projectile, material: catima.Material) ->
e_step = catima.dedx(projectile, material)*x_step*A_recip e_step = catima.dedx(projectile, material)*x_step*A_recip
e_final -= e_step e_final -= e_step
projectile.T(e_final) projectile.T(e_final)
return (e_in - e_final)*projectile.A() return (e_in - e_final)
elif depth == ADAPTIVE_DEPTH_MAX: elif depth == ADAPTIVE_DEPTH_MAX:
return e_in*projectile.A() return e_in
else: else:
e_step = catima.dedx(projectile, material)*x_step*A_recip e_step = catima.dedx(projectile, material)*x_step*A_recip
@ -94,6 +94,7 @@ def get_energyloss(projectile: catima.Projectile, material: catima.Material) ->
class SPSTarget: class SPSTarget:
MEV2U: float = 1.0/931.493614838475
UG2G: float = 1.0e-6 #convert ug to g UG2G: float = 1.0e-6 #convert ug to g
def __init__(self, layers: list[TargetLayer], name: str = "default"): def __init__(self, layers: list[TargetLayer], name: str = "default"):
self.layer_details = layers self.layer_details = layers
@ -114,31 +115,33 @@ class SPSTarget:
if angle == pi*0.5: if angle == pi*0.5:
return e_initial return e_initial
projectile = catima.Projectile(ap, zp) ap_u = ap * self.MEV2U
e_current = e_initial/ap projectile = catima.Projectile(ap_u, zp)
e_current = e_initial/ap_u
for (idx, layer) in enumerate(self.layer_details): 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]) material = catima.Material([(global_nuclear_data.get_data(z, a).mass * self.MEV2U, z, float(s)) for (z, a, s) in layer.compound_list])
projectile.T(e_current) #catima wants MeV/u projectile.T(e_current) #catima wants MeV/u
if idx == rxn_layer: if idx == rxn_layer:
material.thickness(self.layer_details[idx].thickness * self.UG2G / (2.0 * abs(cos(angle)))) material.thickness(self.layer_details[idx].thickness * self.UG2G / (2.0 * abs(cos(angle))))
e_current -= get_energyloss(projectile, material) e_current -= get_energyloss(projectile, material)
return e_initial - e_current*ap return e_initial - e_current*ap_u
else: else:
material.thickness(self.layer_details[idx].thickness * self.UG2G / abs(cos(angle))) material.thickness(self.layer_details[idx].thickness * self.UG2G / abs(cos(angle)))
e_current -= get_energyloss(projectile, material) e_current -= get_energyloss(projectile, material)
return e_initial - e_current*ap return e_initial - e_current*ap_u
#Calculate energy loss for a particle leaving the target, from rxn layer (halfway through rxn layer) to end #Calculate energy loss for a particle leaving the target, from rxn layer (halfway through rxn layer) to end
def get_outgoing_energyloss(self, zp: int, ap: float, e_initial: float, rxn_layer: int, angle: float) -> float: def get_outgoing_energyloss(self, zp: int, ap: float, e_initial: float, rxn_layer: int, angle: float) -> float:
if angle == pi*0.5: if angle == pi*0.5:
return e_initial return e_initial
projectile = catima.Projectile(ap, zp) ap_u = ap * self.MEV2U
e_current = e_initial/ap projectile = catima.Projectile(ap_u, zp)
e_current = e_initial/ap_u
for (idx, layer) in enumerate(self.layer_details[rxn_layer:], start=rxn_layer): for (idx, layer) in enumerate(self.layer_details[rxn_layer:], start=rxn_layer):
material = catima.Material([(global_nuclear_data.get_data(z, a).mass, z, float(s)) for (z, a, s) in layer.compound_list]) material = catima.Material([(global_nuclear_data.get_data(z, a).mass * self.MEV2U, z, float(s)) for (z, a, s) in layer.compound_list])
projectile.T(e_current) #catima wants MeV/u projectile.T(e_current) #catima wants MeV/u
if idx == rxn_layer: if idx == rxn_layer:
material.thickness(self.layer_details[idx].thickness * self.UG2G / (2.0 * abs(cos(angle)))) material.thickness(self.layer_details[idx].thickness * self.UG2G / (2.0 * abs(cos(angle))))
@ -146,19 +149,20 @@ class SPSTarget:
material.thickness(self.layer_details[idx].thickness * self.UG2G / abs(cos(angle))) material.thickness(self.layer_details[idx].thickness * self.UG2G / abs(cos(angle)))
e_current -= get_energyloss(projectile, material) e_current -= get_energyloss(projectile, material)
return e_initial - e_current*ap return e_initial - e_current*ap_u
#Calculate reverse energy loss (energy gain) for a particle that left the target after a reaction (end -> rxn_layer) #Calculate reverse energy loss (energy gain) for a particle that left the target after a reaction (end -> rxn_layer)
def get_outgoing_reverse_energyloss(self, zp: int, ap: float, e_final: float, rxn_layer: int, angle: float) -> float: def get_outgoing_reverse_energyloss(self, zp: int, ap: float, e_final: float, rxn_layer: int, angle: float) -> float:
if angle == pi*0.5: if angle == pi*0.5:
return 0.0 return 0.0
projectile = catima.Projectile(ap, zp) ap_u = ap * self.MEV2U
e_current = e_final/ap projectile = catima.Projectile(ap_u, zp)
e_current = e_final/ap_u
sublist = self.layer_details[rxn_layer:] #only care about rxn_layer -> exit sublist = self.layer_details[rxn_layer:] #only care about rxn_layer -> exit
reveresedRxnLayer = len(sublist) -1 #when reversed rxn_layer is the last layer reveresedRxnLayer = len(sublist) -1 #when reversed rxn_layer is the last layer
for (idx, layer) in reversed(list(enumerate(sublist))): for (idx, layer) in reversed(list(enumerate(sublist))):
material = catima.Material([(global_nuclear_data.get_data(z, a).mass, z, float(s)) for (z, a, s) in layer.compound_list]) material = catima.Material([(global_nuclear_data.get_data(z, a).mass * self.MEV2U, z, float(s)) for (z, a, s) in layer.compound_list])
projectile.T(e_current) #catima wants MeV/u projectile.T(e_current) #catima wants MeV/u
if idx == reveresedRxnLayer: if idx == reveresedRxnLayer:
material.thickness(self.layer_details[idx].thickness * self.UG2G / (2.0 * abs(cos(angle)))) material.thickness(self.layer_details[idx].thickness * self.UG2G / (2.0 * abs(cos(angle))))
@ -166,5 +170,5 @@ class SPSTarget:
material.thickness(self.layer_details[idx].thickness * self.UG2G / abs(cos(angle))) material.thickness(self.layer_details[idx].thickness * self.UG2G / abs(cos(angle)))
e_current += get_reverse_energyloss(projectile, material) e_current += get_reverse_energyloss(projectile, material)
return e_current*ap - e_final return e_current*ap_u - e_final