diff --git a/spspy/SPSReaction.py b/spspy/SPSReaction.py index 20cf052..d4c845d 100644 --- a/spspy/SPSReaction.py +++ b/spspy/SPSReaction.py @@ -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 diff --git a/spspy/SPSTarget.py b/spspy/SPSTarget.py index 6b456c4..f66858f 100644 --- a/spspy/SPSTarget.py +++ b/spspy/SPSTarget.py @@ -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}{s}<\sub>" for z, a, s in self.compound_list]) + return "".join([f"{global_nuclear_data.get_data(z, a,).prettyIsotopicSymbol}{s}" 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 diff --git a/spspy/Spanc.py b/spspy/Spanc.py index debe303..30e5b2d 100644 --- a/spspy/Spanc.py +++ b/spspy/Spanc.py @@ -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: diff --git a/spspy/SpancUI.py b/spspy/SpancUI.py index 7d9a3da..0c5584e 100644 --- a/spspy/SpancUI.py +++ b/spspy/SpancUI.py @@ -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() diff --git a/spspy/data/NuclearData.py b/spspy/data/NuclearData.py index 5e27655..87c69da 100644 --- a/spspy/data/NuclearData.py +++ b/spspy/data/NuclearData.py @@ -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") diff --git a/spspy/ui/PeakDialog.py b/spspy/ui/PeakDialog.py index 782be62..f91e2c0 100644 --- a/spspy/ui/PeakDialog.py +++ b/spspy/ui/PeakDialog.py @@ -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(), diff --git a/spspy/ui/ReactionDialog.py b/spspy/ui/ReactionDialog.py index 2d505d4..a5fde00 100644 --- a/spspy/ui/ReactionDialog.py +++ b/spspy/ui/ReactionDialog.py @@ -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) \ No newline at end of file