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

Compare commits

...

9 Commits

7 changed files with 88 additions and 53 deletions

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)
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(self.params.spsAngle)
return sqrt(residRxnEnergy**2.0 - residRxnP2) - self.residual.mass

View File

@ -13,7 +13,7 @@ class TargetLayer:
thickness: float = 0.0 #ug/cm^2
def __str__(self) -> str:
return "".join([f"{global_nuclear_data.get_data(z, a,).prettyIsotopicSymbol}<sub>{s}<\sub>" for z, a, s in self.compound_list])
return "".join([f"{global_nuclear_data.get_data(z, a,).prettyIsotopicSymbol}<sub>{s}</sub>" for z, a, s in self.compound_list])
#integrate energy loss starting from the final energy and running backwards to initial energy
#catima does not natively provide this type of method
@ -43,10 +43,10 @@ def get_reverse_energyloss(projectile: catima.Projectile, material: catima.Mater
e_step = catima.dedx(projectile, material)*x_step
e_initial += e_step*A_recip
projectile.T(e_initial)
return (e_initial - e_out)*projectile.A()
return (e_initial - e_out)
elif depth == ADAPTIVE_DEPTH_MAX:
return e_out*projectile.A()
return e_out
else:
e_step = catima.dedx(projectile, material)*x_step
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_final -= e_step
projectile.T(e_final)
return (e_in - e_final)*projectile.A()
return (e_in - e_final)
elif depth == ADAPTIVE_DEPTH_MAX:
return e_in*projectile.A()
return e_in
else:
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:
MEV2U: float = 1.0/931.493614838475
UG2G: float = 1.0e-6 #convert ug to g
def __init__(self, layers: list[TargetLayer], name: str = "default"):
self.layer_details = layers
@ -114,31 +115,33 @@ class SPSTarget:
if angle == pi*0.5:
return e_initial
projectile = catima.Projectile(ap, zp)
e_current = e_initial/ap
ap_u = ap * self.MEV2U
projectile = catima.Projectile(ap_u, zp)
e_current = e_initial/ap_u
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
if idx == rxn_layer:
material.thickness(self.layer_details[idx].thickness * self.UG2G / (2.0 * abs(cos(angle))))
e_current -= get_energyloss(projectile, material)
return e_initial - e_current*ap
return e_initial - e_current*ap_u
else:
material.thickness(self.layer_details[idx].thickness * self.UG2G / abs(cos(angle)))
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
def get_outgoing_energyloss(self, zp: int, ap: float, e_initial: float, rxn_layer: int, angle: float) -> float:
if angle == pi*0.5:
return e_initial
projectile = catima.Projectile(ap, zp)
e_current = e_initial/ap
ap_u = ap * self.MEV2U
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):
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
if idx == rxn_layer:
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)))
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)
def get_outgoing_reverse_energyloss(self, zp: int, ap: float, e_final: float, rxn_layer: int, angle: float) -> float:
if angle == pi*0.5:
return 0.0
projectile = catima.Projectile(ap, zp)
e_current = e_final/ap
ap_u = ap * self.MEV2U
projectile = catima.Projectile(ap_u, zp)
e_current = e_final/ap_u
sublist = self.layer_details[rxn_layer:] #only care about rxn_layer -> exit
reveresedRxnLayer = len(sublist) -1 #when reversed rxn_layer is the last layer
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
if idx == reveresedRxnLayer:
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)))
e_current += get_reverse_energyloss(projectile, material)
return e_current*ap - e_final
return e_current*ap_u - e_final

View File

@ -7,6 +7,7 @@ import numpy as np
from enum import Enum
INVALID_PEAK_ID: int = -1
DEG2RAD: float = np.pi / 180.0
class PeakType(Enum):
CALIBRATION = "Calibration"
@ -58,6 +59,7 @@ class Spanc:
print("Cannot create reaction with non-existant target ", targetName)
return
key = f"Rxn{len(self.reactions)}"
params.spsAngle *= DEG2RAD
rxn = Reaction(params, target=self.targets[targetName])
self.reactions[key] = rxn
@ -65,7 +67,7 @@ class Spanc:
if rxnName in self.reactions:
rxn = self.reactions[rxnName]
rxn.params.beamEnergy = beamEnergy
rxn.params.spsAngle = spsAngle
rxn.params.spsAngle = spsAngle * DEG2RAD
rxn.params.magneticField = magneticField
def add_calibration(self, data: Peak) -> None:
@ -77,6 +79,12 @@ class Spanc:
data.peakID = len(self.calibrations)
self.calibrations[data.peakID] = data
return
def remove_calibration(self, data: Peak) -> bool:
if data.peakID not in self.calibrations.keys():
return False
self.calibrations.pop(data.peakID)
return True
def add_output(self, data: Peak) -> None:
if data.peakID == INVALID_PEAK_ID:

View File

@ -1,4 +1,4 @@
from .Spanc import Spanc, PeakType
from .Spanc import Spanc, PeakType, DEG2RAD
from .Fitter import convert_fit_points_to_arrays, convert_resid_points_to_arrays
from .ui.MPLCanvas import MPLCanvas
from .ui.ReactionDialog import ReactionDialog
@ -116,6 +116,7 @@ class SpancGUI(QMainWindow):
self.tablelayout.addWidget(self.targetGroup)
self.targetTable.resizeColumnsToContents()
self.targetTable.cellDoubleClicked.connect(self.handle_update_target)
self.targetTable.setEditTriggers(QTableWidget.EditTrigger.NoEditTriggers)
def create_reaction_table(self) -> None:
self.rxnGroup = QGroupBox("Reactions", self.tableTab)
@ -128,30 +129,33 @@ class SpancGUI(QMainWindow):
self.tablelayout.addWidget(self.rxnGroup)
self.reactionTable.resizeColumnsToContents()
self.reactionTable.cellDoubleClicked.connect(self.handle_update_reaction)
self.reactionTable.setEditTriggers(QTableWidget.EditTrigger.NoEditTriggers)
def create_calibration_table(self) -> None:
self.calGroup = QGroupBox("Calibration Peaks", self.tableTab)
calLayout = QVBoxLayout()
self.calibrationTable = QTableWidget(self.calGroup)
self.calibrationTable.setColumnCount(8)
self.calibrationTable.setHorizontalHeaderLabels(["Reaction","x(mm)","ux stat.(mm)","ux sys.(mm)","rho(cm)","urho(cm)","Ex(MeV)","uEx(MeV)"])
self.calibrationTable.setColumnCount(9)
self.calibrationTable.setHorizontalHeaderLabels(["Peak ID","Reaction","x(mm)","ux stat.(mm)","ux sys.(mm)","rho(cm)","urho(cm)","Ex(MeV)","uEx(MeV)"])
calLayout.addWidget(self.calibrationTable)
self.calGroup.setLayout(calLayout)
self.tablelayout.addWidget(self.calGroup)
self.calibrationTable.resizeColumnsToContents()
self.calibrationTable.cellDoubleClicked.connect(self.handle_update_calibration)
self.calibrationTable.setEditTriggers(QTableWidget.EditTrigger.NoEditTriggers)
def create_output_table(self) -> None:
self.outGroup = QGroupBox("Output Peaks", self.tableTab)
outLayout = QVBoxLayout()
self.outputTable = QTableWidget(self.outGroup)
self.outputTable.setColumnCount(12)
self.outputTable.setHorizontalHeaderLabels(["Reaction","x(mm)","ux stat.(mm)","ux sys.(mm)","rho(cm)","urho(cm)","Ex(MeV)","uEx(MeV)","FWHM(mm)","uFWHM(mm)","FWHM(MeV)","uFWHM(MeV)"])
self.outputTable.setColumnCount(13)
self.outputTable.setHorizontalHeaderLabels(["Peak ID", "Reaction","x(mm)","ux stat.(mm)","ux sys.(mm)","rho(cm)","urho(cm)","Ex(MeV)","uEx(MeV)","FWHM(mm)","uFWHM(mm)","FWHM(MeV)","uFWHM(MeV)"])
outLayout.addWidget(self.outputTable)
self.outGroup.setLayout(outLayout)
self.tablelayout.addWidget(self.outGroup)
self.outputTable.resizeColumnsToContents()
self.outputTable.cellDoubleClicked.connect(self.handle_update_output)
self.outputTable.setEditTriggers(QTableWidget.EditTrigger.NoEditTriggers)
def create_fit_result_text(self) -> None:
self.fitTextGroup = QGroupBox("Fit Results", self.plotTab)
@ -234,11 +238,15 @@ class SpancGUI(QMainWindow):
return
def handle_update_calibration(self, row: int, col: int) -> None:
peakData = self.spanc.calibrations[row]
peakID = int(self.calibrationTable.item(row, 0).text())
peakData = self.spanc.calibrations[peakID]
calDia = PeakDialog(PeakType.CALIBRATION, self.spanc.reactions.keys(), self, peak=peakData)
calDia.new_peak.connect(self.spanc.add_calibration)
calDia.delete_peak.connect(self.spanc.remove_calibration)
if calDia.exec():
self.update_calibration_table()
if self.spanc.isFit == True:
self.handle_run_fit()
return
def handle_new_output(self) -> None:
@ -249,7 +257,8 @@ class SpancGUI(QMainWindow):
return
def handle_update_output(self, row: int, col: int) -> None:
peakData = self.spanc.calibrations[row]
peakID = int(self.calibrationTable.item(row, 0).text())
peakData = self.spanc.outputs[peakID]
outDia = PeakDialog(PeakType.OUTPUT, self.spanc.reactions.keys(), self, peak=peakData)
outDia.new_peak.connect(self.spanc.add_output)
if outDia.exec():
@ -311,7 +320,7 @@ class SpancGUI(QMainWindow):
self.reactionTable.setCellWidget(row, 1, QLabel(str(rxn)))
self.reactionTable.setItem(row, 2, QTableWidgetItem(str(rxn.params.beamEnergy)))
self.reactionTable.setItem(row, 3, QTableWidgetItem(str(rxn.params.magneticField)))
self.reactionTable.setItem(row, 4, QTableWidgetItem(str(rxn.params.spsAngle)))
self.reactionTable.setItem(row, 4, QTableWidgetItem(str(rxn.params.spsAngle / DEG2RAD)))
self.reactionTable.resizeColumnsToContents()
self.reactionTable.resizeRowsToContents()
@ -319,14 +328,15 @@ class SpancGUI(QMainWindow):
self.calibrationTable.setRowCount(len(self.spanc.calibrations))
self.calibrationTable.setVerticalHeaderLabels(self.spanc.calibrations.keys())
for row, peak in enumerate(self.spanc.calibrations.values()):
self.calibrationTable.setItem(row, 0, QTableWidgetItem(peak.rxnName))
self.calibrationTable.setItem(row, 1, QTableWidgetItem(str(peak.position)))
self.calibrationTable.setItem(row, 2, QTableWidgetItem(str(peak.positionErrStat)))
self.calibrationTable.setItem(row, 3, QTableWidgetItem(str(peak.positionErrSys)))
self.calibrationTable.setItem(row, 4, QTableWidgetItem(str(peak.rho)))
self.calibrationTable.setItem(row, 5, QTableWidgetItem(str(peak.rhoErr)))
self.calibrationTable.setItem(row, 6, QTableWidgetItem(str(peak.excitation)))
self.calibrationTable.setItem(row, 7, QTableWidgetItem(str(peak.excitationErr)))
self.calibrationTable.setItem(row, 0, QTableWidgetItem(str(peak.peakID)))
self.calibrationTable.setItem(row, 1, QTableWidgetItem(peak.rxnName))
self.calibrationTable.setItem(row, 2, QTableWidgetItem(str(peak.position)))
self.calibrationTable.setItem(row, 3, QTableWidgetItem(str(peak.positionErrStat)))
self.calibrationTable.setItem(row, 4, QTableWidgetItem(str(peak.positionErrSys)))
self.calibrationTable.setItem(row, 5, QTableWidgetItem(str(peak.rho)))
self.calibrationTable.setItem(row, 6, QTableWidgetItem(str(peak.rhoErr)))
self.calibrationTable.setItem(row, 7, QTableWidgetItem(str(peak.excitation)))
self.calibrationTable.setItem(row, 8, QTableWidgetItem(str(peak.excitationErr)))
self.calibrationTable.resizeColumnsToContents()
self.calibrationTable.resizeRowsToContents()
@ -334,18 +344,19 @@ class SpancGUI(QMainWindow):
self.outputTable.setRowCount(len(self.spanc.outputs))
self.outputTable.setVerticalHeaderLabels(self.spanc.outputs.keys())
for row, peak in enumerate(self.spanc.outputs.values()):
self.outputTable.setItem(row, 0, QTableWidgetItem(peak.rxnName))
self.outputTable.setItem(row, 1, QTableWidgetItem(str(peak.position)))
self.outputTable.setItem(row, 2, QTableWidgetItem(str(peak.positionErrStat)))
self.outputTable.setItem(row, 3, QTableWidgetItem(str(peak.positionErrSys)))
self.outputTable.setItem(row, 4, QTableWidgetItem(str(peak.rho)))
self.outputTable.setItem(row, 5, QTableWidgetItem(str(peak.rhoErr)))
self.outputTable.setItem(row, 6, QTableWidgetItem(str(peak.excitation)))
self.outputTable.setItem(row, 7, QTableWidgetItem(str(peak.excitationErr)))
self.outputTable.setItem(row, 8, QTableWidgetItem(str(peak.positionFWHM)))
self.outputTable.setItem(row, 9, QTableWidgetItem(str(peak.positionFWHMErr)))
self.outputTable.setItem(row, 10, QTableWidgetItem(str(peak.excitationFWHM)))
self.outputTable.setItem(row, 11, QTableWidgetItem(str(peak.excitationFWHMErr)))
self.outputTable.setItem(row, 0, QTableWidgetItem(str(peak.peakID)))
self.outputTable.setItem(row, 1, QTableWidgetItem(peak.rxnName))
self.outputTable.setItem(row, 2, QTableWidgetItem(str(peak.position)))
self.outputTable.setItem(row, 3, QTableWidgetItem(str(peak.positionErrStat)))
self.outputTable.setItem(row, 4, QTableWidgetItem(str(peak.positionErrSys)))
self.outputTable.setItem(row, 5, QTableWidgetItem(str(peak.rho)))
self.outputTable.setItem(row, 6, QTableWidgetItem(str(peak.rhoErr)))
self.outputTable.setItem(row, 7, QTableWidgetItem(str(peak.excitation)))
self.outputTable.setItem(row, 8, QTableWidgetItem(str(peak.excitationErr)))
self.outputTable.setItem(row, 9, QTableWidgetItem(str(peak.positionFWHM)))
self.outputTable.setItem(row, 10, QTableWidgetItem(str(peak.positionFWHMErr)))
self.outputTable.setItem(row, 11, QTableWidgetItem(str(peak.excitationFWHM)))
self.outputTable.setItem(row, 12, QTableWidgetItem(str(peak.excitationFWHMErr)))
self.outputTable.resizeColumnsToContents()
self.outputTable.resizeRowsToContents()

View File

@ -53,7 +53,7 @@ def get_excitations(Z: int, A: int) -> list[float]:
levels = []
text = ''
symbol = global_nuclear_data.get_data(Z, A).isotopicSymbol
site = req.get(f"https://www.nndc.bnl.gov/nudat2/getdatasetClassic.jsp?nucleus={symbol}&unc=nds")
site = req.get(f"https://www.nndc.bnl.gov/nudat3/getdatasetClassic.jsp?nucleus={symbol}&unc=nds")
contents = xhtml.fromstring(site.content)
tables = contents.xpath("//table")
rows = tables[2].xpath("./tr")

View File

@ -2,11 +2,12 @@ from ..Spanc import Peak, PeakType, INVALID_PEAK_ID
from PySide6.QtWidgets import QLabel
from PySide6.QtWidgets import QVBoxLayout, QFormLayout, QGroupBox
from PySide6.QtWidgets import QComboBox, QDoubleSpinBox
from PySide6.QtWidgets import QDialog, QDialogButtonBox
from PySide6.QtWidgets import QDialog, QDialogButtonBox, QPushButton
from PySide6.QtCore import Signal
class PeakDialog(QDialog):
new_peak = Signal(Peak)
delete_peak = Signal(Peak)
def __init__(self, peakType: PeakType, rxnList: list[str], parent=None, peak: Peak=None):
super().__init__(parent)
@ -26,14 +27,19 @@ class PeakDialog(QDialog):
if peakType == PeakType.CALIBRATION:
self.create_calibration_inputs()
if peak is not None:
self.setWindowTitle(f"Update A {peakType.value} Peak")
self.set_calibration_inputs(peak)
self.peakID = peak.peakID
self.buttonBox.accepted.connect(self.send_update_calibration_peak)
self.deleteButton = QPushButton("Delete", self)
self.deleteButton.clicked.connect(self.send_delete_calibration_peak)
self.centralLayout.addWidget(self.deleteButton)
else:
self.buttonBox.accepted.connect(self.send_calibration_peak)
elif peakType == PeakType.OUTPUT:
self.create_output_inputs()
if peak is not None:
self.setWindowTitle(f"Update A {peakType.value} Peak")
self.set_output_inputs(peak)
self.peakID = peak.peakID
self.buttonBox.accepted.connect(self.send_update_output_peak)
@ -58,8 +64,8 @@ class PeakDialog(QDialog):
self.uexInput = QDoubleSpinBox(self.inputGroupBox)
self.uexInput.setDecimals(6)
inputLayout.addRow("Position(mm)", self.xInput)
inputLayout.addRow("Position Stat. Error(mm)", self.uxstatInput)
inputLayout.addRow("Position Sys. Error(mm)", self.uxsysInput)
inputLayout.addRow("Position Stat. Error(mm)", self.uxstatInput)
inputLayout.addRow("Excitation Energy(MeV)", self.exInput)
inputLayout.addRow("Excitation Energy Error(MeV)", self.uexInput)
self.inputGroupBox.setLayout(inputLayout)
@ -82,8 +88,8 @@ class PeakDialog(QDialog):
self.ufwhmInput = QDoubleSpinBox(self.inputGroupBox)
self.ufwhmInput.setDecimals(6)
inputLayout.addRow("Position(mm)", self.xInput)
inputLayout.addRow("Position Stat. Error(mm)", self.uxstatInput)
inputLayout.addRow("Position Sys. Error(mm)", self.uxsysInput)
inputLayout.addRow("Position Stat. Error(mm)", self.uxstatInput)
inputLayout.addRow("Position FWHM(mm)", self.fwhmInput)
inputLayout.addRow("Position FWHM Error(mm)", self.ufwhmInput)
self.inputGroupBox.setLayout(inputLayout)
@ -114,6 +120,12 @@ class PeakDialog(QDialog):
peak = Peak(excitation=self.exInput.value(), excitationErr=self.uexInput.value(), position=self.xInput.value(),
positionErrStat=self.uxstatInput.value(), positionErrSys=self.uxsysInput.value(), rxnName=self.rxnNameBox.currentText(), peakID=self.peakID)
self.new_peak.emit(peak)
def send_delete_calibration_peak(self) -> None:
peak = Peak(excitation=self.exInput.value(), excitationErr=self.uexInput.value(), position=self.xInput.value(),
positionErrStat=self.uxstatInput.value(), positionErrSys=self.uxsysInput.value(), rxnName=self.rxnNameBox.currentText(), peakID=self.peakID)
self.delete_peak.emit(peak)
self.done(3)
def send_output_peak(self) -> None:
peak = Peak(position=self.xInput.value(), positionErrStat=self.uxstatInput.value(), positionErrSys=self.uxsysInput.value(),

View File

@ -5,6 +5,7 @@ from PySide6.QtWidgets import QVBoxLayout, QFormLayout, QGroupBox
from PySide6.QtWidgets import QSpinBox, QDoubleSpinBox, QComboBox
from PySide6.QtWidgets import QDialog, QDialogButtonBox
from PySide6.QtCore import Signal
from numpy import rad2deg
MAXIMUM_NUCLEAR_Z: int = 110
MAXIMUM_NUCLEAR_A: int = 270
@ -116,6 +117,6 @@ class ReactionDialog(QDialog):
self.aeInput.setValue(rxn.params.ejectile.A)
self.aeInput.setEnabled(False)
self.bkeInput.setValue(rxn.params.beamEnergy)
self.thetaInput.setValue(rxn.params.spsAngle)
self.thetaInput.setValue(rad2deg(rxn.params.spsAngle))
self.bfieldInput.setValue(rxn.params.magneticField)