Repository navigation
Conversation
There was a problem hiding this comment.
🟡 Changes recommended
The PSTAR function derives its maximum from number_of_data_points, a field that actually stores the material enum value, so it returns an incorrect (near-minimum) upper bound for every material.
5 open issues
What changed in this PR
This PR adds two new public helper functions, AT_PSTAR_energy_bounds and AT_ICRU_energy_bounds, that report the minimum and maximum tabulated energies available in the PSTAR and ICRU stopping-power datasets for a given material (and, for ICRU, particle). These are intended as lookup utilities alongside the existing AT_PSTAR_wrapper/AT_ICRU_wrapper interpolation routines in the stopping-power data modules.
Changes:
- Added
AT_PSTAR_energy_bounds(material bounds) in the PSTAR source/header. - Added
AT_ICRU_energy_bounds(material + particle bounds, with a special 250 MeV/u cap for helium) in the ICRU source/header. - The PSTAR implementation derives the upper bound from
number_of_data_points, which does not actually contain the point count.
| File | Description |
|---|---|
| src/AT_StoppingPowerDataPSTAR.c | Implements AT_PSTAR_energy_bounds; incorrectly uses number_of_data_points (holds the material id, not the count) for the max index. |
| src/AT_StoppingPowerDataICRU.c | Implements AT_ICRU_energy_bounds with material/particle validation and He upper-bound cap. |
| include/AT_StoppingPowerDataPSTAR.h | Declares AT_PSTAR_energy_bounds (no Doxygen). |
| include/AT_StoppingPowerDataICRU.h | Declares AT_ICRU_energy_bounds (no Doxygen). |
Review effort: Balanced
… documentation and error handling
|
Thanks for this! I built the branch and ran both functions: ICRU gives 0.025–1000 MeV/u (He capped at 250) for both materials, PSTAR gives 0.001–10000 for materials 1–9, and material 0 is correctly rejected. Before merging, I'd suggest a few changes. 1. API shape: two functions returning
|
…energies directly, with error handling for invalid material numbers
…nergies directly, with error handling for invalid material numbers
grzanka
left a comment
grzanka
left a comment
There was a problem hiding this comment.
Thanks for the update! Several points from my earlier comment are done: scalar double functions returning -1.0, PSTAR bounds converted to MeV/u, the per-material Z limit for Al2O3, and rejecting Z < 1. I found one crash and a few things that are still open. Details and proposed code are in the inline comments.
I tried all the changes below locally on top of 011eead: the library builds with no new warnings, amtrack_test passes with the new checks, and pyamtrack #234 (with the small changes at the end of this comment) builds against it and passes all 121 tests. The full patch is at the bottom if you want to apply it directly.
Must fix
- Segfault for material 0 in
AT_PSTAR_min/max_energy_MeV_u(comment onAT_StoppingPowerDataPSTAR.c).
Still open from my earlier comment
- API through the source number:
AT_Stopping_Power_E_MeV_u_min/max(source, particle, material)(comments onAT_StoppingPower.handAT_StoppingPower.c). - One "has data" check per table, used everywhere: by
AT_Stopping_Power_has_data, the min/max functions and the wrappers (comments on both data files). - 250 MeV/u He limit as one named constant (comment on
AT_StoppingPowerDataICRU.c). - Tests: see below, since
test/C/AT_test.cis not part of this diff.
Smaller
- Argument order
(particle_no, material_no)like the rest of the library (comments on the headers). has_datafor Bethe returns 1 for any input (comment onAT_StoppingPower.c).- 4-space indentation in a file that uses tabs; missing newline at the end of
AT_StoppingPowerDataPSTAR.c. - The PR description still describes the old
*_energy_boundsout-pointer API; please update it.
Tests: test/C/AT_test.c
CI already runs amtrack_test, but it only prints values and always returns success. With a small check() helper and a non-zero exit code on failure, CI catches wrong values, not only crashes:
#include <math.h>
static int failures = 0;
static void check(const int condition, const char* description) {
if (!condition) {
printf("FAILED: %s\n", description);
failures++;
}
}
static void check_stopping_power_energy_range() {
const long proton = PARTICLE_PROTON_NUMBER;
const long helium = 2004;
const long neon = 10020;
const long argon = 18040;
check(AT_Stopping_Power_E_MeV_u_min(ICRU, proton, Water_Liquid) == 0.025, "ICRU min, proton in water");
check(AT_Stopping_Power_E_MeV_u_max(ICRU, proton, Water_Liquid) == 1000.0, "ICRU max, proton in water");
check(AT_Stopping_Power_E_MeV_u_max(ICRU, helium, Water_Liquid) == ICRU_HELIUM_MAX_E_MEV_U, "ICRU max, He-4 in water");
check(AT_Stopping_Power_has_data(ICRU, neon, Aluminum_Oxide), "ICRU data for Ne-20 in Al2O3");
check(!AT_Stopping_Power_has_data(ICRU, argon, Aluminum_Oxide), "no ICRU data for Ar-40 in Al2O3");
check(AT_Stopping_Power_E_MeV_u_max(ICRU, argon, Aluminum_Oxide) == -1.0, "ICRU max, Ar-40 in Al2O3");
check(!AT_Stopping_Power_has_data(PSTAR, proton, User_Defined_Material), "no PSTAR data for material 0");
check(AT_Stopping_Power_E_MeV_u_min(PSTAR, proton, User_Defined_Material) == -1.0, "PSTAR min, material 0");
check(AT_Stopping_Power_E_MeV_u_max(PSTAR, proton, Tungsten) == -1.0, "PSTAR max, material without table");
/* PSTAR tables are in MeV: the MeV/u limit must map back to the 10000 MeV table end */
check(fabs(AT_E_MeV_from_E_MeV_u(AT_Stopping_Power_E_MeV_u_max(PSTAR, proton, Water_Liquid), proton) - 10000.0) < 1e-9,
"PSTAR max, MeV/u conversion");
check(AT_Stopping_Power_E_MeV_u_min(Bethe, proton, Water_Liquid) == BETHE_LOWER_LIMIT_E_MEV_U, "Bethe min, proton in water");
check(!AT_Stopping_Power_has_data(Bethe, proton, 999), "no Bethe data for unknown material");
check(!AT_Stopping_Power_has_data(FromFile, proton, Water_Liquid), "FromFile is never reported as available");
}Call check_stopping_power_energy_range(); at the start of main, and end main with:
if (failures > 0) {
return EXIT_FAILURE;
}
return EXIT_SUCCESS;The "material 0" check would have caught the segfault.
pyamtrack #234
With the source-based functions, pyamtrack needs only AT_StoppingPower.h, and the PSTAR and ICRU branches become one:
extern "C" {
#include "AT_StoppingPower.h"
-#include "AT_StoppingPowerDataICRU.h"
-#include "AT_StoppingPowerDataPSTAR.h"
} case static_cast<long>(StoppingPowerSource::PSTAR):
case static_cast<long>(StoppingPowerSource::ICRU):
if (!AT_Stopping_Power_has_data(source, particle_no, material_no)) {
throw std::invalid_argument(std::string(stopping_power_source_name(source)) +
" data is not available for particle " + std::to_string(particle_no) +
" in material ID " + std::to_string(material_no));
}
return StoppingPowerEnergyBounds{AT_Stopping_Power_E_MeV_u_min(source, particle_no, material_no),
AT_Stopping_Power_E_MeV_u_max(source, particle_no, material_no)};Why:
- The
minimum < 0.0 || maximum < 0.0re-checks are not needed oncehas_datareturned 1. - The current message for Ar-40 in Al2O3 says "ICRU data is not available for material ID 2", which points at the material although the particle is the problem. Naming both makes the cause clear.
- The Bethe branch can stay
std::nullopt, so pyamtrack's behaviour for Bethe does not change.
Two tests for tests/test_stopping.py that pin the new limits (the existing PSTAR tests use 0.000999 and 100000 MeV, far from the table end):
@pytest.mark.parametrize("function", [pyamtrack.stopping.mass_stopping_power, pyamtrack.stopping.stopping_power])
def test_pstar_range_ends_at_the_last_table_energy(function):
# PSTAR tables end at 10000 MeV for protons
assert function(10000.0, particle=PROTON, source="pstar") > 0
with pytest.raises(ValueError, match="outside the PSTAR range"):
function(10000.5, particle=PROTON, source="pstar")
@pytest.mark.parametrize("function", [pyamtrack.stopping.mass_stopping_power, pyamtrack.stopping.stopping_power])
def test_icru_rejects_particles_without_data_in_aluminum_oxide(function):
# ICRU data for Al2O3 cover only Z <= 10
argon = pyamtrack.particles.from_string("40Ar")
with pytest.raises(ValueError, match="ICRU data is not available"):
function(400.0, particle=argon, material=pyamtrack.materials.aluminum_oxide, source="icru")pyamtrack CI fails at the moment only because it builds against library master, which does not have the new functions yet, so this PR needs to be merged first.
Open question
For Bethe I propose min = BETHE_LOWER_LIMIT_E_MEV_U (0.49 MeV/u; below it AT_Bethe_wrapper returns 0) and max = HUGE_VAL (Infinity in JS, inf in Python). The alternative is -1.0 for Bethe, but then has_data says 1 while min/max say "no data". Happy to hear other ideas.
Full patch (applies on top of 011eead with git apply)
diff --git a/include/AT_StoppingPower.h b/include/AT_StoppingPower.h
index 7cd34bd..073d177 100644
--- a/include/AT_StoppingPower.h
+++ b/include/AT_StoppingPower.h
@@ -276,13 +276,41 @@ double AT_Mass_Stopping_Power_E_MeV_u_single( const long stopping_power_source_n
/**
* Checks whether stopping-power data are available for a source,
* particle, and material combination.
+ * For FromFile no file name is given, so 0 is always returned.
*
- * @param[in] stopping_power_source_no stopping-power data-source id
- * @param[in] particle_no particle number
- * @param[in] material_no material number
- * @return 1 if data are available, or 0 otherwise
+ * @param[in] stopping_power_source_no id of the data source
+ * @param[in] particle_no particle number
+ * @param[in] material_no material number
+ * @return 1 if data are available, 0 otherwise
+ */
+int AT_Stopping_Power_has_data( const long stopping_power_source_no,
+ const long particle_no,
+ const long material_no);
+
+/**
+ * Returns the lowest kinetic energy in MeV/u for which the source
+ * has stopping-power data for the given particle and material.
+ *
+ * @param[in] stopping_power_source_no id of the data source
+ * @param[in] particle_no particle number
+ * @param[in] material_no material number
+ * @return minimum kinetic energy in MeV/u, or -1.0 if there are no data
+ */
+double AT_Stopping_Power_E_MeV_u_min( const long stopping_power_source_no,
+ const long particle_no,
+ const long material_no);
+
+/**
+ * Returns the highest kinetic energy in MeV/u for which the source
+ * has stopping-power data for the given particle and material.
+ * For Bethe there is no upper limit and HUGE_VAL is returned.
+ *
+ * @param[in] stopping_power_source_no id of the data source
+ * @param[in] particle_no particle number
+ * @param[in] material_no material number
+ * @return maximum kinetic energy in MeV/u, or -1.0 if there are no data
*/
-int AT_Stopping_Power_has_data(const long stopping_power_source_no,
+double AT_Stopping_Power_E_MeV_u_max( const long stopping_power_source_no,
const long particle_no,
const long material_no);
diff --git a/include/AT_StoppingPowerDataICRU.h b/include/AT_StoppingPowerDataICRU.h
index 7101bb0..422cfa2 100644
--- a/include/AT_StoppingPowerDataICRU.h
+++ b/include/AT_StoppingPowerDataICRU.h
@@ -50,25 +50,38 @@ int AT_ICRU_wrapper( const long n,
const char info[],
double mass_stopping_power_MeV_cm2_g[]);
+/** ICRU 49 data for helium end at 250 MeV/u (1000 MeV for He-4) */
+#define ICRU_HELIUM_MAX_E_MEV_U 250.0
+
+/**
+ * Checks whether ICRU data are available for a particle in a material.
+ *
+ * @param[in] particle_no particle number
+ * @param[in] material_no material number
+ * @return 1 if data are available, 0 otherwise
+ */
+int AT_ICRU_has_data( const long particle_no,
+ const long material_no);
+
/**
- * @brief Returns the minimum kinetic energy covered by the ICRU data.
+ * Returns the lowest kinetic energy covered by the ICRU data.
*
- * @param[in] material_no Material number identifying the material.
- * @param[in] particle_no Particle number identifying the particle.
- * @return Minimum tabulated kinetic energy in MeV/u, or -1.0 if
- * material_no or particle_no has no ICRU data.
+ * @param[in] particle_no particle number
+ * @param[in] material_no material number
+ * @return minimum kinetic energy in MeV/u, or -1.0 if there are no ICRU data
*/
-double AT_ICRU_min_energy_MeV_u(long material_no, long particle_no);
+double AT_ICRU_E_MeV_u_min( const long particle_no,
+ const long material_no);
/**
- * @brief Returns the maximum kinetic energy covered by the ICRU data.
+ * Returns the highest kinetic energy covered by the ICRU data.
*
- * @param[in] material_no Material number identifying the material.
- * @param[in] particle_no Particle number identifying the particle.
- * @return Maximum tabulated kinetic energy in MeV/u, or -1.0 if
- * material_no or particle_no has no ICRU data.
+ * @param[in] particle_no particle number
+ * @param[in] material_no material number
+ * @return maximum kinetic energy in MeV/u, or -1.0 if there are no ICRU data
*/
-double AT_ICRU_max_energy_MeV_u(long material_no, long particle_no);
+double AT_ICRU_E_MeV_u_max( const long particle_no,
+ const long material_no);
/**
* @struct AT_stopping_power_ICRU_table
diff --git a/include/AT_StoppingPowerDataPSTAR.h b/include/AT_StoppingPowerDataPSTAR.h
index 4cd5d94..3668ab6 100644
--- a/include/AT_StoppingPowerDataPSTAR.h
+++ b/include/AT_StoppingPowerDataPSTAR.h
@@ -66,24 +66,41 @@ int AT_PSTAR_wrapper( const long n,
const char info[],
double mass_stopping_power_MeV_cm2_g[]);
+/**
+ * Checks whether PSTAR data are available for a particle in a material.
+ * Particles other than protons are scaled from the proton data, so any
+ * valid particle number is accepted.
+ *
+ * @param[in] particle_no particle number
+ * @param[in] material_no material number
+ * @return 1 if data are available, 0 otherwise
+ */
+int AT_PSTAR_has_data( const long particle_no,
+ const long material_no);
/**
- * @brief Returns the minimum kinetic energy covered by the PSTAR data.
+ * Returns the lowest kinetic energy covered by the PSTAR data.
+ * The PSTAR tables are given in MeV for protons; AT_PSTAR_wrapper looks
+ * them up at the proton energy with the same E_MeV_u, so the limit in MeV/u
+ * is the same for all particles.
*
- * @param[in] material_no Material number identifying the material.
- * @return Minimum tabulated kinetic energy in MeV/u, or -1.0 if
- * material_no is outside the PSTAR material range.
+ * @param[in] particle_no particle number
+ * @param[in] material_no material number
+ * @return minimum kinetic energy in MeV/u, or -1.0 if there are no PSTAR data
*/
-double AT_PSTAR_min_energy_MeV_u(long material_no);
+double AT_PSTAR_E_MeV_u_min( const long particle_no,
+ const long material_no);
/**
- * @brief Returns the maximum kinetic energy covered by the PSTAR data.
+ * Returns the highest kinetic energy covered by the PSTAR data.
+ * See AT_PSTAR_E_MeV_u_min for the unit conversion.
*
- * @param[in] material_no Material number identifying the material.
- * @return Maximum tabulated kinetic energy in MeV/u, or -1.0 if
- * material_no is outside the PSTAR material range.
+ * @param[in] particle_no particle number
+ * @param[in] material_no material number
+ * @return maximum kinetic energy in MeV/u, or -1.0 if there are no PSTAR data
*/
-double AT_PSTAR_max_energy_MeV_u(long material_no);
+double AT_PSTAR_E_MeV_u_max( const long particle_no,
+ const long material_no);
/** PSTAR data downloaded from NIST database: http://www.nist.gov/pml/data/star/index.cfm
* Stopping-Power and Range Tables
diff --git a/src/AT_StoppingPower.c b/src/AT_StoppingPower.c
index 01431c5..14092f0 100644
--- a/src/AT_StoppingPower.c
+++ b/src/AT_StoppingPower.c
@@ -319,34 +319,61 @@ double AT_Mass_Stopping_Power_E_MeV_u_single( const long stopping_power_source_n
return stopping_power_MeV_cm2_g;
}
-int AT_Stopping_Power_has_data(const long stopping_power_source_no,
- const long particle_no,
- const long material_no)
-{
- switch (stopping_power_source_no) {
- case FromFile:
- /*
- * The filename and file contents are not provided here,
- * so availability cannot be checked.
- */
- return 0;
- case Bethe:
- /* Bethe is an analytical model, not a finite data table. */
- return 1;
- case PSTAR:
- if (AT_Z_from_particle_no_single(particle_no) < 1 ||
- AT_A_from_particle_no_single(particle_no) < 1) {
- return 0;
- }
- return material_no >= Water_Liquid &&
- material_no <= Copper;
- case ICRU:
- if (AT_Z_from_particle_no_single(particle_no) < 1 ||
- AT_A_from_particle_no_single(particle_no) < 1) {
- return 0;
- }
- return AT_ICRU_min_energy_MeV_u(material_no, particle_no) >= 0.0;
- default:
- return 0;
- }
+int AT_Stopping_Power_has_data( const long stopping_power_source_no,
+ const long particle_no,
+ const long material_no){
+
+ switch(stopping_power_source_no){
+ case Bethe:
+ /* analytical formula: needs a valid particle and a known material */
+ return AT_Z_from_particle_no_single(particle_no) >= 1 &&
+ AT_material_index_from_material_number(material_no) >= 0;
+ case PSTAR:
+ return AT_PSTAR_has_data(particle_no, material_no);
+ case ICRU:
+ return AT_ICRU_has_data(particle_no, material_no);
+ default:
+ /* FromFile: no file name is given here, so availability cannot be checked */
+ return 0;
+ }
+}
+
+double AT_Stopping_Power_E_MeV_u_min( const long stopping_power_source_no,
+ const long particle_no,
+ const long material_no){
+
+ switch(stopping_power_source_no){
+ case Bethe:
+ if(!AT_Stopping_Power_has_data(Bethe, particle_no, material_no)){
+ return -1.0;
+ }
+ /* below this limit AT_Bethe_wrapper returns zero */
+ return BETHE_LOWER_LIMIT_E_MEV_U;
+ case PSTAR:
+ return AT_PSTAR_E_MeV_u_min(particle_no, material_no);
+ case ICRU:
+ return AT_ICRU_E_MeV_u_min(particle_no, material_no);
+ default:
+ return -1.0;
+ }
+}
+
+double AT_Stopping_Power_E_MeV_u_max( const long stopping_power_source_no,
+ const long particle_no,
+ const long material_no){
+
+ switch(stopping_power_source_no){
+ case Bethe:
+ if(!AT_Stopping_Power_has_data(Bethe, particle_no, material_no)){
+ return -1.0;
+ }
+ /* analytical formula, no upper limit */
+ return HUGE_VAL;
+ case PSTAR:
+ return AT_PSTAR_E_MeV_u_max(particle_no, material_no);
+ case ICRU:
+ return AT_ICRU_E_MeV_u_max(particle_no, material_no);
+ default:
+ return -1.0;
+ }
}
diff --git a/src/AT_StoppingPowerDataICRU.c b/src/AT_StoppingPowerDataICRU.c
index a7f729a..5e7c685 100644
--- a/src/AT_StoppingPowerDataICRU.c
+++ b/src/AT_StoppingPowerDataICRU.c
@@ -139,12 +139,30 @@ static long AT_ICRU_max_atomic_number(const long material_no) {
return 0;
}
-static int AT_ICRU_has_data_for_atomic_number(
- const long material_no,
- const long atomic_number) {
- const long max_atomic_number =
- AT_ICRU_max_atomic_number(material_no);
- return atomic_number >= 1 && atomic_number <= max_atomic_number;
+int AT_ICRU_has_data( const long particle_no,
+ const long material_no){
+ const long Z = AT_Z_from_particle_no_single(particle_no);
+ return Z >= 1 && Z <= AT_ICRU_max_atomic_number(material_no);
+}
+
+double AT_ICRU_E_MeV_u_min( const long particle_no,
+ const long material_no){
+ if(!AT_ICRU_has_data(particle_no, material_no)){
+ return -1.0;
+ }
+ return AT_stopping_power_ICRU_table[material_no-1].energy_and_stopping_power[0][0];
+}
+
+double AT_ICRU_E_MeV_u_max( const long particle_no,
+ const long material_no){
+ if(!AT_ICRU_has_data(particle_no, material_no)){
+ return -1.0;
+ }
+ if(AT_Z_from_particle_no_single(particle_no) == 2){
+ return ICRU_HELIUM_MAX_E_MEV_U;
+ }
+ const AT_stopping_power_ICRU_table_struct * table = &AT_stopping_power_ICRU_table[material_no-1];
+ return table->energy_and_stopping_power[0][table->number_of_data_points - 1];
}
int AT_ICRU_wrapper( const long n,
@@ -155,7 +173,7 @@ int AT_ICRU_wrapper( const long n,
double mass_stopping_power_MeV_cm2_g[]){
long i;
- if((material_no != Water_Liquid) && (material_no != Aluminum_Oxide)){ /* Water and Alox data only */
+ if(AT_ICRU_max_atomic_number(material_no) == 0){ /* Water and Alox data only */
for(i = 0; i < n; i++){
mass_stopping_power_MeV_cm2_g[i] = -1.0;
}
@@ -167,17 +185,13 @@ int AT_ICRU_wrapper( const long n,
for(i = 0; i < n; i++){
- if(!AT_ICRU_has_data_for_atomic_number(material_no, Z[i])){
- mass_stopping_power_MeV_cm2_g[i] = -1.0;
- continue;
- }
-
- if (E_MeV_u[i] < 0.025 || E_MeV_u[i] > 1000){
+ if(!AT_ICRU_has_data(particle_no[i], material_no)){
mass_stopping_power_MeV_cm2_g[i] = -1.0;
continue;
}
- if (Z[i] == 2 && E_MeV_u[i] > 250){ // Data for He only until 250 MeV/u = 1000 MeV !
+ if (E_MeV_u[i] < AT_ICRU_E_MeV_u_min(particle_no[i], material_no) ||
+ E_MeV_u[i] > AT_ICRU_E_MeV_u_max(particle_no[i], material_no)){
mass_stopping_power_MeV_cm2_g[i] = -1.0;
continue;
}
@@ -196,51 +210,3 @@ int AT_ICRU_wrapper( const long n,
free(Z);
return AT_Success;
}
-
-static int AT_ICRU_energy_bounds(long material_no,
- long particle_no,
- double *minimum,
- double *maximum) {
- if (material_no != Water_Liquid &&
- material_no != Aluminum_Oxide) {
- return AT_No_ICRU_Data;
- }
-
- long atomic_number = AT_Z_from_particle_no_single(particle_no);
- if (!AT_ICRU_has_data_for_atomic_number(material_no, atomic_number)) {
- return AT_No_ICRU_Data;
- }
-
- const AT_stopping_power_ICRU_table_struct *table =
- &AT_stopping_power_ICRU_table[material_no - 1];
-
- *minimum = table->energy_and_stopping_power[0][0];
- *maximum = table->energy_and_stopping_power[0]
- [table->number_of_data_points - 1];
-
- if (atomic_number == 2) {
- *maximum = 250.0;
- }
-
- return AT_Success;
-}
-
-double AT_ICRU_min_energy_MeV_u(long material_no, long particle_no) {
- double minimum;
- double maximum;
- if (AT_ICRU_energy_bounds(material_no, particle_no,
- &minimum, &maximum) != AT_Success) {
- return -1.0;
- }
- return minimum;
-}
-
-double AT_ICRU_max_energy_MeV_u(long material_no, long particle_no) {
- double minimum;
- double maximum;
- if (AT_ICRU_energy_bounds(material_no, particle_no,
- &minimum, &maximum) != AT_Success) {
- return -1.0;
- }
- return maximum;
-}
diff --git a/src/AT_StoppingPowerDataPSTAR.c b/src/AT_StoppingPowerDataPSTAR.c
index f26c595..7e1f9e9 100644
--- a/src/AT_StoppingPowerDataPSTAR.c
+++ b/src/AT_StoppingPowerDataPSTAR.c
@@ -1290,6 +1290,38 @@ PSTAR_data_struct PSTAR_data =
&PSTAR_Copper } };
+/* NULL if material_no is out of range or has no PSTAR table (e.g. User_Defined_Material) */
+static const PSTAR_data_for_material_struct * AT_PSTAR_table( const long material_no){
+ if( material_no < 0 || material_no >= N_PSTAR_MATERIALS ){
+ return NULL;
+ }
+ return PSTAR_data.stopping_power_source_data[material_no];
+}
+
+int AT_PSTAR_has_data( const long particle_no,
+ const long material_no){
+ return AT_Z_from_particle_no_single(particle_no) >= 1 &&
+ AT_PSTAR_table(material_no) != NULL;
+}
+
+double AT_PSTAR_E_MeV_u_min( const long particle_no,
+ const long material_no){
+ if(!AT_PSTAR_has_data(particle_no, material_no)){
+ return -1.0;
+ }
+ return AT_E_MeV_u_from_E_MeV(AT_PSTAR_table(material_no)->energy_and_stopping_power[0][0],
+ PARTICLE_PROTON_NUMBER);
+}
+
+double AT_PSTAR_E_MeV_u_max( const long particle_no,
+ const long material_no){
+ if(!AT_PSTAR_has_data(particle_no, material_no)){
+ return -1.0;
+ }
+ return AT_E_MeV_u_from_E_MeV(AT_PSTAR_table(material_no)->energy_and_stopping_power[N_PSTAR_DATAPOINTS - 1][0],
+ PARTICLE_PROTON_NUMBER);
+}
+
int AT_PSTAR_wrapper( const long n,
const double E_MeV_u[],
const long particle_no[],
@@ -1297,7 +1329,7 @@ int AT_PSTAR_wrapper( const long n,
const char info[],
double mass_stopping_power_MeV_cm2_g[]){
- const PSTAR_data_for_material_struct * source_for_given_material = PSTAR_data.stopping_power_source_data[material_no];
+ const PSTAR_data_for_material_struct * source_for_given_material = AT_PSTAR_table(material_no);
if( source_for_given_material == NULL ){
char material_name[MATERIAL_NAME_LENGTH];
@@ -1341,38 +1373,3 @@ int AT_PSTAR_wrapper( const long n,
return AT_Success;
}
-
-static int AT_PSTAR_energy_bounds(long material_no,
- double *minimum,
- double *maximum) {
- if (material_no < 0 || material_no >= N_PSTAR_MATERIALS) {
- return AT_No_PSTAR_Data;
- }
-
- *minimum =
- PSTAR_data.stopping_power_source_data[material_no]
- ->energy_and_stopping_power[0][0];
- *maximum =
- PSTAR_data.stopping_power_source_data[material_no]
- ->energy_and_stopping_power[N_PSTAR_DATAPOINTS - 1][0];
-
- return AT_Success;
-}
-
-double AT_PSTAR_min_energy_MeV_u(long material_no) {
- double minimum;
- double maximum;
- if (AT_PSTAR_energy_bounds(material_no, &minimum, &maximum) != AT_Success) {
- return -1.0;
- }
- return AT_E_MeV_u_from_E_MeV(minimum, PARTICLE_PROTON_NUMBER);
-}
-
-double AT_PSTAR_max_energy_MeV_u(long material_no) {
- double minimum;
- double maximum;
- if (AT_PSTAR_energy_bounds(material_no, &minimum, &maximum) != AT_Success) {
- return -1.0;
- }
- return AT_E_MeV_u_from_E_MeV(maximum, PARTICLE_PROTON_NUMBER);
-}
\ No newline at end of file
diff --git a/test/C/AT_test.c b/test/C/AT_test.c
index 1423aea..ef95974 100644
--- a/test/C/AT_test.c
+++ b/test/C/AT_test.c
@@ -30,12 +30,49 @@
#include <stdio.h>
#include <stdlib.h>
+#include <math.h>
#include "AT_StoppingPower.h"
#include "AT_DataRange.h"
+static int failures = 0;
+
+static void check(const int condition, const char* description) {
+ if (!condition) {
+ printf("FAILED: %s\n", description);
+ failures++;
+ }
+}
+
+static void check_stopping_power_energy_range() {
+ const long proton = PARTICLE_PROTON_NUMBER;
+ const long helium = 2004;
+ const long neon = 10020;
+ const long argon = 18040;
+
+ check(AT_Stopping_Power_E_MeV_u_min(ICRU, proton, Water_Liquid) == 0.025, "ICRU min, proton in water");
+ check(AT_Stopping_Power_E_MeV_u_max(ICRU, proton, Water_Liquid) == 1000.0, "ICRU max, proton in water");
+ check(AT_Stopping_Power_E_MeV_u_max(ICRU, helium, Water_Liquid) == ICRU_HELIUM_MAX_E_MEV_U, "ICRU max, He-4 in water");
+ check(AT_Stopping_Power_has_data(ICRU, neon, Aluminum_Oxide), "ICRU data for Ne-20 in Al2O3");
+ check(!AT_Stopping_Power_has_data(ICRU, argon, Aluminum_Oxide), "no ICRU data for Ar-40 in Al2O3");
+ check(AT_Stopping_Power_E_MeV_u_max(ICRU, argon, Aluminum_Oxide) == -1.0, "ICRU max, Ar-40 in Al2O3");
+
+ check(!AT_Stopping_Power_has_data(PSTAR, proton, User_Defined_Material), "no PSTAR data for material 0");
+ check(AT_Stopping_Power_E_MeV_u_min(PSTAR, proton, User_Defined_Material) == -1.0, "PSTAR min, material 0");
+ check(AT_Stopping_Power_E_MeV_u_max(PSTAR, proton, Tungsten) == -1.0, "PSTAR max, material without table");
+ /* PSTAR tables are in MeV: the MeV/u limit must map back to the 10000 MeV table end */
+ check(fabs(AT_E_MeV_from_E_MeV_u(AT_Stopping_Power_E_MeV_u_max(PSTAR, proton, Water_Liquid), proton) - 10000.0) < 1e-9,
+ "PSTAR max, MeV/u conversion");
+
+ check(AT_Stopping_Power_E_MeV_u_min(Bethe, proton, Water_Liquid) == BETHE_LOWER_LIMIT_E_MEV_U, "Bethe min, proton in water");
+ check(!AT_Stopping_Power_has_data(Bethe, proton, 999), "no Bethe data for unknown material");
+ check(!AT_Stopping_Power_has_data(FromFile, proton, Water_Liquid), "FromFile is never reported as available");
+}
+
int main() {
+ check_stopping_power_energy_range();
+
double Test = AT_CSDA_range_g_cm2_single(270.55,
0.1,
6012,
@@ -73,6 +110,9 @@ int main() {
test);
printf("ICRU: %e, %e, %e\n", test[0], test[1], test[2]);
+ if (failures > 0) {
+ return EXIT_FAILURE;
+ }
return EXIT_SUCCESS;
};
| static int AT_PSTAR_energy_bounds(long material_no, | ||
| double *minimum, | ||
| double *maximum) { | ||
| if (material_no < 0 || material_no >= N_PSTAR_MATERIALS) { |
There was a problem hiding this comment.
Crash for material 0 (User_Defined_Material). This range check lets material_no = 0 through, but slot 0 of PSTAR_data.stopping_power_source_data is NULL, so the next lines dereference a null pointer. AT_PSTAR_min_energy_MeV_u(0) segfaults (I reproduced it). Commit 46131eb had a table == NULL check; b1df4d1 dropped it.
I'd put the table lookup in one small helper before AT_PSTAR_wrapper:
/* NULL if material_no is out of range or has no PSTAR table (e.g. User_Defined_Material) */
static const PSTAR_data_for_material_struct * AT_PSTAR_table( const long material_no){
if( material_no < 0 || material_no >= N_PSTAR_MATERIALS ){
return NULL;
}
return PSTAR_data.stopping_power_source_data[material_no];
}and use it in the wrapper too:
- const PSTAR_data_for_material_struct * source_for_given_material = PSTAR_data.stopping_power_source_data[material_no];
+ const PSTAR_data_for_material_struct * source_for_given_material = AT_PSTAR_table(material_no);Why the wrapper as well: it indexes stopping_power_source_data[material_no] with no range check, so e.g. Tungsten (10) reads past the end of the array. That bug is older than this PR, but with the helper both problems are fixed in one place: the wrapper then returns AT_No_PSTAR_Data for material 10.
| return AT_Success; | ||
| } | ||
|
|
||
| double AT_PSTAR_min_energy_MeV_u(long material_no) { |
There was a problem hiding this comment.
With AT_PSTAR_table, the private AT_PSTAR_energy_bounds helper and its out-pointers are not needed. Each function checks once and returns:
int AT_PSTAR_has_data( const long particle_no,
const long material_no){
return AT_Z_from_particle_no_single(particle_no) >= 1 &&
AT_PSTAR_table(material_no) != NULL;
}
double AT_PSTAR_E_MeV_u_min( const long particle_no,
const long material_no){
if(!AT_PSTAR_has_data(particle_no, material_no)){
return -1.0;
}
return AT_E_MeV_u_from_E_MeV(AT_PSTAR_table(material_no)->energy_and_stopping_power[0][0],
PARTICLE_PROTON_NUMBER);
}
double AT_PSTAR_E_MeV_u_max( const long particle_no,
const long material_no){
if(!AT_PSTAR_has_data(particle_no, material_no)){
return -1.0;
}
return AT_E_MeV_u_from_E_MeV(AT_PSTAR_table(material_no)->energy_and_stopping_power[N_PSTAR_DATAPOINTS - 1][0],
PARTICLE_PROTON_NUMBER);
}Why:
AT_PSTAR_has_databecomes the only place that decides whether PSTAR has data.AT_Stopping_Power_has_datacan call it instead of repeating the material range (see the comment there).- No out-pointers means no unused
minimum/maximumlocals and nothing left uninitialised on the error path.
The MeV to MeV/u conversion itself is right: the max is 9921.62 MeV/u, which the wrapper turns back into exactly the 10000 MeV table end.
| return -1.0; | ||
| } | ||
| return AT_E_MeV_u_from_E_MeV(maximum, PARTICLE_PROTON_NUMBER); | ||
| } No newline at end of file |
There was a problem hiding this comment.
Please add a newline at the end of the file. Without it, the next change to this file shows this line as modified too.
| * @return Minimum tabulated kinetic energy in MeV/u, or -1.0 if | ||
| * material_no is outside the PSTAR material range. | ||
| */ | ||
| double AT_PSTAR_min_energy_MeV_u(long material_no); |
There was a problem hiding this comment.
Please give the PSTAR and ICRU functions the same signature, with particle_no first:
int AT_PSTAR_has_data( const long particle_no,
const long material_no);
double AT_PSTAR_E_MeV_u_min( const long particle_no,
const long material_no);
double AT_PSTAR_E_MeV_u_max( const long particle_no,
const long material_no);Why:
(particle_no, material_no)is the order of every other stopping-power function (AT_PSTAR_wrapper,AT_Mass_Stopping_Power_E_MeV_u_single, ...). The ICRU functions in this PR have it reversed, and since both arelongthe compiler cannot catch a swapped call.- The PSTAR bounds don't depend on the particle, but taking
particle_noletshas_datareject invalid particle numbers, and letsAT_Stopping_Power_E_MeV_u_min/maxcall both tables the same way. E_MeV_u_min/maxmatches the names of the source-based functions (AT_Stopping_Power_E_MeV_u_min/max).
In the docs, please say that the tables are in MeV for protons and that AT_PSTAR_wrapper looks them up at the proton energy for the same E_MeV_u. That is why the MeV/u limit is the same for all particles, which is not obvious from the code.
| * @return Minimum tabulated kinetic energy in MeV/u, or -1.0 if | ||
| * material_no or particle_no has no ICRU data. | ||
| */ | ||
| double AT_ICRU_min_energy_MeV_u(long material_no, long particle_no); |
There was a problem hiding this comment.
Same as for PSTAR: (particle_no, material_no) order and E_MeV_u_min/max names. I'd also declare the "has data" check and the He limit here:
/** ICRU 49 data for helium end at 250 MeV/u (1000 MeV for He-4) */
#define ICRU_HELIUM_MAX_E_MEV_U 250.0
int AT_ICRU_has_data( const long particle_no,
const long material_no);
double AT_ICRU_E_MeV_u_min( const long particle_no,
const long material_no);
double AT_ICRU_E_MeV_u_max( const long particle_no,
const long material_no);Why: the constant replaces the two hard-coded 250s (see the comments in the .c file). AT_ICRU_has_data has to be public so that AT_Stopping_Power_has_data can call it directly, instead of inferring availability from min >= 0.
| [table->number_of_data_points - 1]; | ||
|
|
||
| if (atomic_number == 2) { | ||
| *maximum = 250.0; |
There was a problem hiding this comment.
Please use the named constant ICRU_HELIUM_MAX_E_MEV_U here (see the header comment).
Why: 250 is currently written both here and in AT_ICRU_wrapper. If someone changes one, the reported max and the wrapper disagree. Once the wrapper calls AT_ICRU_E_MeV_u_max (comment on line 170), the constant is used in exactly one place.
| * @param[in] material_no material number | ||
| * @return 1 if data are available, or 0 otherwise | ||
| */ | ||
| int AT_Stopping_Power_has_data(const long stopping_power_source_no, |
There was a problem hiding this comment.
Please declare the source-based min/max here, next to has_data:
/**
* Returns the lowest kinetic energy in MeV/u for which the source
* has stopping-power data for the given particle and material.
*
* @param[in] stopping_power_source_no id of the data source
* @param[in] particle_no particle number
* @param[in] material_no material number
* @return minimum kinetic energy in MeV/u, or -1.0 if there are no data
*/
double AT_Stopping_Power_E_MeV_u_min( const long stopping_power_source_no,
const long particle_no,
const long material_no);
/**
* Returns the highest kinetic energy in MeV/u for which the source
* has stopping-power data for the given particle and material.
* For Bethe there is no upper limit and HUGE_VAL is returned.
*
* (same parameters)
* @return maximum kinetic energy in MeV/u, or -1.0 if there are no data
*/
double AT_Stopping_Power_E_MeV_u_max( const long stopping_power_source_no,
const long particle_no,
const long material_no);Why: this header is the public stopping-power API. With these functions, pyamtrack (and the JS package) needs only this header and the source number it already passes to AT_Mass_Stopping_Power_with_no. At the moment pyamtrack has to include the ICRU and PSTAR data headers and pick the table-specific function itself.
Please also say in the has_data doc that FromFile always returns 0. Otherwise that is a surprise for callers.
| return stopping_power_MeV_cm2_g; | ||
| } | ||
|
|
||
| int AT_Stopping_Power_has_data(const long stopping_power_source_no, |
There was a problem hiding this comment.
I'd make has_data a pure dispatcher to the per-table checks, and add the two min/max functions after it:
int AT_Stopping_Power_has_data( const long stopping_power_source_no,
const long particle_no,
const long material_no){
switch(stopping_power_source_no){
case Bethe:
/* analytical formula: needs a valid particle and a known material */
return AT_Z_from_particle_no_single(particle_no) >= 1 &&
AT_material_index_from_material_number(material_no) >= 0;
case PSTAR:
return AT_PSTAR_has_data(particle_no, material_no);
case ICRU:
return AT_ICRU_has_data(particle_no, material_no);
default:
/* FromFile: no file name is given here, so availability cannot be checked */
return 0;
}
}
double AT_Stopping_Power_E_MeV_u_min( const long stopping_power_source_no,
const long particle_no,
const long material_no){
switch(stopping_power_source_no){
case Bethe:
if(!AT_Stopping_Power_has_data(Bethe, particle_no, material_no)){
return -1.0;
}
/* below this limit AT_Bethe_wrapper returns zero */
return BETHE_LOWER_LIMIT_E_MEV_U;
case PSTAR:
return AT_PSTAR_E_MeV_u_min(particle_no, material_no);
case ICRU:
return AT_ICRU_E_MeV_u_min(particle_no, material_no);
default:
return -1.0;
}
}
double AT_Stopping_Power_E_MeV_u_max( const long stopping_power_source_no,
const long particle_no,
const long material_no){
switch(stopping_power_source_no){
case Bethe:
if(!AT_Stopping_Power_has_data(Bethe, particle_no, material_no)){
return -1.0;
}
/* analytical formula, no upper limit */
return HUGE_VAL;
case PSTAR:
return AT_PSTAR_E_MeV_u_max(particle_no, material_no);
case ICRU:
return AT_ICRU_E_MeV_u_max(particle_no, material_no);
default:
return -1.0;
}
}Why: each table's rules (materials, Z range, energy range) then live only in that table's file, and this file just picks the table by source number. That is the same pattern AT_Mass_Stopping_Power_with_no already uses for the stopping power itself. It also replaces the AT_ICRU_min_energy_MeV_u(...) >= 0.0 check below with a direct AT_ICRU_has_data call, which says what it means.
The A >= 1 check is dropped because none of the three sources use A. The PSTAR and ICRU wrappers take E_MeV_u and use only Z.
| return 0; | ||
| case Bethe: | ||
| /* Bethe is an analytical model, not a finite data table. */ | ||
| return 1; |
There was a problem hiding this comment.
return 1 for any input means has_data(Bethe, 0, 999) is 1. But Bethe needs a valid particle (it uses Z for the effective charge) and a material from the material table (it reads the average Z and A from it). With material 999, AT_Mass_Stopping_Power_E_MeV_u_single(Bethe, 10.0, 1001, 999) aborts on assert( A > 0 ) in a debug build.
Why it matters: callers such as pyamtrack use has_data to decide whether a call can succeed. A clear 0 is better than 1 followed by an abort. The code in the comment above checks both the particle and the material.
| return 0; | ||
| } | ||
| return material_no >= Water_Liquid && | ||
| material_no <= Copper; |
There was a problem hiding this comment.
This repeats the list of PSTAR materials from AT_StoppingPowerDataPSTAR.c.
Why it's a problem: if a PSTAR table is added, or a slot is left NULL, this range and the actual data disagree, and nothing fails at compile time. AT_PSTAR_has_data(particle_no, material_no) looks at the table itself, so the list stays in one place.
This pull request introduces new API functions to retrieve the minimum and maximum kinetic energy bounds from the ICRU and PSTAR stopping-power tables, as well as their implementations and related header documentation. These changes improve the usability of the API by allowing users to easily query the valid energy ranges for materials and particles.
API Additions:
AT_ICRU_energy_boundstoAT_StoppingPowerDataICRU.h, which returns the minimum and maximum kinetic energy per nucleon for a given material and particle, along with detailed documentation.AT_PSTAR_energy_boundstoAT_StoppingPowerDataPSTAR.h, which returns the minimum and maximum kinetic energy for a given material, along with detailed documentation.Implementations:
AT_ICRU_energy_boundsinAT_StoppingPowerDataICRU.c, including input validation and logic to extract the energy bounds from the data tables, with a special case for helium ions.AT_PSTAR_energy_boundsinAT_StoppingPowerDataPSTAR.c, including input validation and logic to extract the energy bounds from the PSTAR data tables.Dependency Updates:
#include "AT_Error.h"to bothAT_StoppingPowerDataICRU.candAT_StoppingPowerDataPSTAR.cto support the new error codes used in the API additions. [1] [2]