Skip to content

Handling isomers produced beyond (n,n*) reactions#229

Open
eitan-weinstein wants to merge 48 commits intosvalinn:mainfrom
eitan-weinstein:non_n_isomer_rxns
Open

Handling isomers produced beyond (n,n*) reactions#229
eitan-weinstein wants to merge 48 commits intosvalinn:mainfrom
eitan-weinstein:non_n_isomer_rxns

Conversation

@eitan-weinstein
Copy link
Copy Markdown
Contributor

This PR expands on the isomer handling methods for ALARAJOYWrapper developed over the past month to properly address isomer daughters produced from pathways other than (n,n*) reactions. I came across this issue when doing inventory benchmarking against FISPACT-II and the JAEA FNS experimental data and noticed certain isomers appearing to be stable when they obviously were not, just like before. The old method failed to include cases where the parent itself is an isomer, which can subsequently produce excited daughters by normal reaction pathways.

Copy link
Copy Markdown
Member

@gonuke gonuke left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

a few questions and comments

"""

return M == 0
return M == 0 and dKZA % 10 == 0
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should we have set M correctly before this in these circumstances?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The way that we define M in reaction_data.py determines if a daughter is left in an excited state by specific emitted particle tags from the MT string formatting. However, that formatting does not account for cases where the parent is already excited, meaning the del_KZA value would calculate an excited daughter as well. M is generally parent-ambivalent, so I think I'll need to do a parent-by-parent recalculation of M at the start of the MT for loop to account for excited parents.

_is_ground_state_nuc(rxn['isomer']) or
_is_isomer_with_decay_data(dKZA, radionucs, rxn['isomer'])
_is_ground_state_daughter(rxn['isomer'], dKZA) or
_is_isomer_with_decay_data(dKZA, pKZA, radionucs, rxn['isomer'])
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

do we care if its an isomer anymore? given this or we only care that it has decay data, right?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think we can only go decay data because then we would be excluding all stable nuclides from the standard processing and directing them incorrectly to the special isomer handling in the else block.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But all stable nuclides are in the ground state, so they should be handled by the first condition.

I guess my point is, we've already captured all the ground state nuclides in the first condition, so we know that any nuclides we are testing must be an isomer

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, I agree with what you are saying about stable nuclides all being handled by the first condition, but I'm not sure I understand what change you are requesting.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

_is_isomer_with_decay_data can just be refactored to be _has_decay_data and doesn't need to check if it's an isomer

else:
if to_ground:
decay_KZA = f'{pKZA // 10}*'
div_KZA = pKZA if rxn['emitted'] == 'n' else dKZA
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

cant we always use dkza? for an emitted neutron, I assume pkza == dkza

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I think you're right. Although, after closer inspection, I realized the assignment of decay_KZA was reversed, based on the actual meaning of the to_ground Boolean, which I'll have updated in the next commit.

return sigma_list[::-1]

def _is_ground_state_nuc(M):
def _is_ground_state_daughter(M):
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't this it's correct to name this about daughters. You are only passing in M which could theoretically come from anywhere. The fact that you may only use it for daughters should not affect its name.

_is_ground_state_nuc(rxn['isomer']) or
_is_isomer_with_decay_data(dKZA, radionucs, rxn['isomer'])
_is_ground_state_daughter(rxn['isomer'], dKZA) or
_is_isomer_with_decay_data(dKZA, pKZA, radionucs, rxn['isomer'])
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But all stable nuclides are in the ground state, so they should be handled by the first condition.

I guess my point is, we've already captured all the ground state nuclides in the first condition, so we know that any nuclides we are testing must be an isomer

Copy link
Copy Markdown
Member

@gonuke gonuke left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

continued cleanup of the isomer checking

"""

return dKZA in radionucs and M in range(1,10)
return dKZA in radionucs and ((pKZA % 10) + M in range(1,10))
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

do we still need to test the isomeric state of the parent?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In principle, yes, although I don't know if we'd ever actually run into a case where this would come up, which would be having a parent whose excitation level is 10 and somehow still has decay data. This is obviously very unlikely to happen, given that I don't think there is decay data beyond any M=2, but technically this covers a physical case that could exist.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But didn't we just decide that this method only cares whether there is decay data or not?

This is called in one place: to decide whether we keep the cross-section as-is or collapse it. It seems to me that we should keep all cross-sections where the daughter is ground state (you have a function for that) or where the daughter has decay data (this function). It shouldn't matter whether it's an isomer or not.

Copy link
Copy Markdown
Contributor Author

@eitan-weinstein eitan-weinstein Feb 5, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe I'm not defining the function clearly enough, but I think the problem that I am trying to address here is that the KZA arithmetic allows for the production of daughters that could appear to have decay data, but only as a result of the limitations of the KZA convention, not because the actual daughter produced is represented in the radionucs dictionary. The simplest example of this would an (n,n10) reaction of a ground state parent, like what I have in the _has_decay_data() docstring: 6Li(n,n10) -> 6,m10Li + n, but in KZA notation it would be 30060 + n10 -> 300610 + n. The problem here is that the KZA of the daughter does correspond to a nuclide with decay data, just not the correct daughter because 300610 = 61Zn. Perhaps my implementation is not capturing this accurately, but what I am trying to capture the parameters of the problem of not allowing any double digit excited states because the isomeric properties inescapably affect the legitimacy of a KZA value, and by extension our ability to use that KZA to cross-reference the decay data.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we should be creating some clear metadata about these cases sooner then, perhaps? As soon as we determine what the daughter KZA is, we should also set some metadata to capture this.

else:
if to_ground:
decay_KZA = f'{pKZA // 10}*'
decay_KZA = (dKZA // 10 - rxn['isomer'] // 10) * 10
Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I realized that the version of this formula that I had in previous commits for this PR was incorrect.

Previously it was:

decay_KZA = (dKZA // 10) * 10

When benchmarking the FNS data, I was noticing a number of non-physical reaction pathways represented (specifically it was 51V(n,n*)->53V which caught my attention as being obviously incorrect). I realized that the reason for this was that for double digit excited states, this rounding-down to 10 method would not account for the dKZA already being corrupted by the unallowed KZA operation with a double-digit M. Therefore, I changed it to:

decay_KZA = (dKZA // 10 - rxn['isomer'] // 10) * 10

to bring the tens-digit back to appropriate value after any isomeric corruption may occur (unaffected by any M < 10)

Copy link
Copy Markdown
Member

@gonuke gonuke left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The suggestion we discussed today

Copy link
Copy Markdown
Member

@gonuke gonuke left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you provide a handful of examples of the kinds of reactions that need the special treatment of the isomer data? It seems like we keep ratcheting up the value of or isomeric information.

Comment on lines +279 to +280
if parent_excitation > 0 and MT not in EXCITATION_REACTIONS:
M += parent_excitation
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This surprises me a little - that we're adding the parent excitation to the isomeric state that we already had in the MT table? Or does this only occur when rxn['isomer'] == 0?


# Force isomeric daughters with high excitation levels or without known
# half-lives to their ground state, and flag the change
if pKZA % 10 + M >= 10 or dKZA not in radionucs:
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Below, we've already added the parent excitation to M and it seems like we do it again here?

@eitan-weinstein
Copy link
Copy Markdown
Contributor Author

eitan-weinstein commented Feb 10, 2026

Per my Thread in the Slack channel, I think on top of figuring out whether isomeric parents' reaction residuals default to the ground state or the equivalent isomeric state, I also think that I need to recalibrate the forced-ground accumulation in general, because I think that I was overlooking that MTs 4 and 103-107 were already accumulated, so to calculate just the ground state, I would have to subtract from the accumulated isomeric cross-sections instead.

Copy link
Copy Markdown
Member

@gonuke gonuke left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The cognitive burden of our processing of these excitation reactions is getting hard to follow. I've requested that we pull some parts out into separate functions to make it clearer.

if MT in excitation_pathways[pKZA]:
return excitation_pathways[pKZA][MT], excitation_pathways
else:
return 0, excitation_pathways
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(Magic number)
Let's define a variable

default_state  = 0

And return that instead of just returning 0 so that we can easily change it later.

# Parse File 10 for allowed daughter excitations from a given reaction
excitation_pathways[pKZA] = {}

za = f'{(pKZA // 10) * 10 + delKZA}'[:-1]
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I need to be reminded (again?) when we every produce pKZA values that need this treatment...

I suppose that we sometimes have a parent nuclide in an isomeric state and we want to block that isomeric state when we apply the MT-specific change in KZA. This is particularly true given recent discussion of not assuming that reaction products are in the same isomeric state as their parents.

for line in section.content.split('\n'):
if za in line:
excitation_levels.append(int(
line.split(za)[1].strip().split(' ')[0]
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we add a comment somewhere with the format of the line we are parsing here? This is compact/concise, but a little obtuse.

Comment on lines +337 to +338
# Reset deLKZA for (n,n*) excitation reactions
if delKZA >= 0 and MT in EXCITATION_REACTIONS:
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The EXCITATION_REACTIONS are no longer just (n,n*) so don't some of them still have a non-zero delKZA?

return excitation_pathways[pKZA][MT], excitation_pathways

def iterate_MTs(MTs, file_obj, mt_dict, pKZA, all_rxns, radionucs, to_ground):
def iterate_MTs(MTs, file_obj, mt_dict, pKZA, all_rxns, eaf_nucs, endf_path):
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it will lower cognitive burden if we break some of the code in this method into different functions.

Originally, the purpose of this function was just to iterate through the data and extract simple cross-sections for writing to the final file.

Now we do a lot of cleanup and cross-checking. It would be great if those steps were done in their own functions to make the overall flow of this method more clear and isolate where we perform those operations.

We might end up passing around a lot of data (which makes me feel like we should probably think about making this object oriented some day but not now).

Copy link
Copy Markdown
Member

@gonuke gonuke left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've noted a few places where it may be helpful to put things into functions to make it more readable/lower cognitive burden

This preprocessor is designed to update ALARA's data input capabilities for updated FENDL3.2x data sets formatted as TENDL (ENDF-6) files. Unlike previous versions of [FENDL](https://www-nds.iaea.org/fendl_library/websites/fendl32b/) (Fusion Evaluated Nuclear Data Library), which have available groupwise cross-section data for neutron activation, as is required for ALARA's functionality, FENDL3 data requires conversion to that format.

This preprocessor uses [NJOY 2016](https://github.com/njoy/NJOY2016) Nuclear Data Processing System to produce requisite pointwise data (PENDF) to subsequntly convert to the Vitamin-J 175 energy group groupwise format (GENDF) for activation cross-sections and to handle the processed data to be fed back to ALARA. For pointwise conversion, ALARAJOYWrapper uses the NJOY modules MODER, RECONR, BROADR, UNRESR, and GASPR, and for the ultimate conversion to the groupwise format, GROUPR.
This preprocessor uses [NJOY 2016](https://github.com/njoy/NJOY2016) Nuclear Data Processing System to produce requisite pointwise data (PENDF) to subsequntly convert to the Vitamin-J 175 energy group groupwise format (GENDF) for activation cross-sections and to handle the processed data to be fed back to ALARA. For pointwise conversion, ALARAJOYWrapper uses the NJOY modules MODER, RECONR, BROADR, UNRESR, and GASPR, and for the ultimate conversion to the groupwise format, GROUPR. Note that in order to properly handle isomeric reactions, NJOY 2016 has to be built with the `develop` branch, which contains recently developed capacity for holistic isomer handling.
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you include a specific git hash here? At some point it will no longer be necessary to use the develop branch if/when they issue a release, and it would be good to try to capture that here. The best way is to note a specific SHA and then when there is a new release, we can ensure that the new capability is in that release and change this reference to that version.

Comment on lines +320 to +327
# Stable nuclide processing
if 'STABLE NUCLIDE' in line and prev_line is not None:
iso = prev_line.rstrip().split('DECAY')[0]

if 'FILE' not in iso:
Z, _, A = iso.split('-')
kza = (int(Z) * 1000 + int(A)) * 10
radionucs.setdefault(kza, 0.0)
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you need to add these in this way? Is there a reason that this dictionary can't only contain radioactive nuclides and the absence of a nuclide implies that it's stable?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The reason I chose to do it this way would be to reduce ambiguity later on when checking a KZA's validity for inclusion in the DSV without special handling. In that case, it does not matter whether a nuclide is radioactive or stable, but rather whether it is referenceable in the library because otherwise, we would have to create an additional reference for all stable nuclides, which I think would make less sense than the current approach.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But this assumes that the EAF decay library is an authoritative reference on all stable nuclides? Seems outside of its scope/purpose...

Comment on lines +308 to +324
pKZA = tp.extract_pendf_pkza(pendf_path)
isomer_dict = tp.determine_all_excitations(
TAPE20, MTs, pKZA, mt_dict
)
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could we put all of this into a function that returns the few things that we ultimately need when it's done? This will help readability.

Comment on lines 327 to 354
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perhaps this also all goes into another function that returns all_rxns?

@@ -326,29 +318,29 @@ def main():
_, pendf_MTs, _ = tp.extract_endf_specs(pendf_path)
gas_MTs = set(pendf_MTs) & set(rxd.GAS_DF['total_mt'])
MTs |= {int(gas_MT) for gas_MT in gas_MTs}
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Which of these are not already integers? Should we promote them to integers sooner? or does that break something else?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm actually not sure why I have that int setting here. They all are already, so I'll just take it out.

Comment on lines +170 to +190
za = str((pKZA + mt_dict[MT]['delKZA']) // 10).zfill(6)
has_isomers = False
# Isomer pathways contained either in MF 9 ("Multiplicities for
# Production of Radioactive Nuclides") and MF 10 ("Cross Sections
# for Production of Radioactive Nuclides").
for MF in [9,10]:
if isomer_dict and len(isomer_dict[MT][MF]) > 1:
has_isomers = True
for M in isomer_dict[MT][MF]:
if M < 10:
card9_lines.append(
f'{3 if MF == 9 else 4}{za}{M} ' \
f'{MT} "{mtname}, M={M}" /'
)
else:
# Special handling for high LFS
# Requires develop branch of NJOY2016
card9_lines.append(
'-1 / \n' \
f' {MF} {MT} {za} {M} "{mtname}, M={M}" /'
)
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe this belongs in a function for readability?

Copy link
Copy Markdown
Member

@gonuke gonuke left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Some cleanup from removing the stables from eaf_dict

f.readline()

_in_decay_block = False
prev_line = None
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think you need prev_line anymore

_in_decay_block = False

# Advance to next line
prev_line = line
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Don't need prev_line anymore

Copy link
Copy Markdown
Member

@gonuke gonuke left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Mostly some questions about the "Card 9" stuff...

This preprocessor is designed to update ALARA's data input capabilities for updated FENDL3.2x data sets formatted as TENDL (ENDF-6) files. Unlike previous versions of [FENDL](https://www-nds.iaea.org/fendl_library/websites/fendl32b/) (Fusion Evaluated Nuclear Data Library), which have available groupwise cross-section data for neutron activation, as is required for ALARA's functionality, FENDL3 data requires conversion to that format.

This preprocessor uses [NJOY 2016](https://github.com/njoy/NJOY2016) Nuclear Data Processing System to produce requisite pointwise data (PENDF) to subsequntly convert to the Vitamin-J 175 energy group groupwise format (GENDF) for activation cross-sections and to handle the processed data to be fed back to ALARA. For pointwise conversion, ALARAJOYWrapper uses the NJOY modules MODER, RECONR, BROADR, UNRESR, and GASPR, and for the ultimate conversion to the groupwise format, GROUPR.
This preprocessor uses [NJOY 2016](https://github.com/njoy/NJOY2016) Nuclear Data Processing System to produce requisite pointwise data (PENDF) to subsequntly convert to the Vitamin-J 175 energy group groupwise format (GENDF) for activation cross-sections and to handle the processed data to be fed back to ALARA. For pointwise conversion, ALARAJOYWrapper uses the NJOY modules MODER, RECONR, BROADR, UNRESR, and GASPR, and for the ultimate conversion to the groupwise format, GROUPR. Note that in order to properly handle isomeric reactions, NJOY 2016 has to be built with from the commit `db71977593d084ae5bbb9e5c88a926541718d313` of the `develop` branch, which includes the necessary support for holistic isomer handling.
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
This preprocessor uses [NJOY 2016](https://github.com/njoy/NJOY2016) Nuclear Data Processing System to produce requisite pointwise data (PENDF) to subsequntly convert to the Vitamin-J 175 energy group groupwise format (GENDF) for activation cross-sections and to handle the processed data to be fed back to ALARA. For pointwise conversion, ALARAJOYWrapper uses the NJOY modules MODER, RECONR, BROADR, UNRESR, and GASPR, and for the ultimate conversion to the groupwise format, GROUPR. Note that in order to properly handle isomeric reactions, NJOY 2016 has to be built with from the commit `db71977593d084ae5bbb9e5c88a926541718d313` of the `develop` branch, which includes the necessary support for holistic isomer handling.
This preprocessor uses [NJOY 2016](https://github.com/njoy/NJOY2016) Nuclear Data Processing System to produce requisite pointwise data (PENDF) to subsequntly convert to the Vitamin-J 175 energy group groupwise format (GENDF) for activation cross-sections and to handle the processed data to be fed back to ALARA. For pointwise conversion, ALARAJOYWrapper uses the NJOY modules MODER, RECONR, BROADR, UNRESR, and GASPR, and for the ultimate conversion to the groupwise format, GROUPR. Note that in order to properly handle isomeric reactions, NJOY 2016 has to be built with from the commit `db71977593d084ae5bbb9e5c88a926541718d313` or newer of the `develop` branch, which includes the necessary support for holistic isomer handling.

has_isomers = False

for MF in [9,10]:
if isomer_dict and len(isomer_dict[MT][MF]) > 1:
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This seems like a strange place to check for isomer_dict. If it isn't provided, then the whole method can be skipped, right? So we can either check for it immediately and exit if it's not there, or we can only call the method if/when it is defined.

for M in isomer_dict[MT][MF]:
if M < 10:
card9_lines.append(
f'{3 if MF == 9 else 4}{za}{M} ' \
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

While I don't understand the NJOY syntax, can't we do this with math?

Suggested change
f'{3 if MF == 9 else 4}{za}{M} ' \
f'{MF-6}{za}{M} ' \

f'{3 if MF == 9 else 4}{za}{M} ' \
f'{MT} "{mtname}, M={M}" /'
)
else:
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is the new capability only valid for M >= 10? or only necessary for M >= 10? If it is valid for M < 10, then why not use it for all values of M?

"""

za = str((pKZA + mt_dict[MT]['delKZA']) // 10).zfill(6)
has_isomers = False
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Whether or not has_isomers is true appears to be based entirely on isomer_dict which is an input to this method - can't we determine this outside of this method?

Then, this could just return a set of cards to return and then append/extend them outside this method.

Copy link
Copy Markdown
Member

@gonuke gonuke left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is a lot in this change and I wish I could better grapple with it all in my head. Here are a few simplifications, but it is still a little hard to read the whole thing.

Comment on lines +225 to +235
if pKZA:
if has_isomers:
card9_lines = write_card9_special_isomer_reactions(
pKZA, mt_dict, isomer_dict, MT, card9_lines, mtname
)

else:
card9_lines.append(f'{MFD} {MT} "{mtname}" /')

else:
card9_lines.append(f'{MFD} {MT} "{mtname}" /')
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this is the same...

Suggested change
if pKZA:
if has_isomers:
card9_lines = write_card9_special_isomer_reactions(
pKZA, mt_dict, isomer_dict, MT, card9_lines, mtname
)
else:
card9_lines.append(f'{MFD} {MT} "{mtname}" /')
else:
card9_lines.append(f'{MFD} {MT} "{mtname}" /')
if pKZA and has_isomers:
card9_lines = write_card9_special_isomer_reactions(
pKZA, mt_dict, isomer_dict, MT, card9_lines, mtname
)
else:
card9_lines.append(f'{MFD} {MT} "{mtname}" /')

Comment on lines +222 to +235
has_isomers = (
MT in isomer_dict and isomer_dict[MT].keys() != {MFD}
)
if pKZA:
if has_isomers:
card9_lines = write_card9_special_isomer_reactions(
pKZA, mt_dict, isomer_dict, MT, card9_lines, mtname
)

else:
card9_lines.append(f'{MFD} {MT} "{mtname}" /')

else:
card9_lines.append(f'{MFD} {MT} "{mtname}" /')
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I feel as though I may be thrashing back and forth in my guidance, but can we do all of this in the special card9 routine?

card9_lines.append(write_card9_special_isomer_reactions(pKZA, MT, mt_dict[MT], isomer_dict[MT], mtname, MFD))

Copy link
Copy Markdown
Contributor Author

@eitan-weinstein eitan-weinstein Mar 24, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To clarify the "specialness" of the isomer handling, but still keeping all of the card9_lines inclusions in one-line, what if I went with something like this instead:

card9_lines.extend(write_card9_special_isomer_reactions(
      pKZA, mt_dict[MT], isomer_dict[MT], MT, mtname
) if pKZA and has_isomers else [f'{MFD} {MT} "{mtname}" /'])

(Note, I'm using extend here instead of append following the changes to write_card9_special_isomer_reactions() from your earlier comment on restructuring the function itself. If we're not operating on the original card9_lines list within the function, I think it makes most sense to be extending the list with the return of the function).

Copy link
Copy Markdown
Member

@gonuke gonuke Mar 25, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's not so much about wanting it to be a single line as wanting there to be more clarity about what's happening. In theory that clarity could be in comments, but with a lot of if/then logic, comments get a little confusing.

The current structure has the question of "whether there are special isomers" separate from "what to do when there are special isomers".

I'd like to contain the "whether there are special isomers" in some kind of method so that it communicates what's going on. My first suggestion was to put it in the card9 handling method so that becomes a method that generates card9, including the "whether" question and the "what to do" question.

Comment on lines +163 to +171
za = str((pKZA + mt_dict[MT]['delKZA']) // 10).zfill(6)
for MF in isomer_dict[MT].keys():
if len(isomer_dict[MT][MF]) > 1:
for M in isomer_dict[MT][MF]:
card9_lines.append(
f'-1 / \n {MF} {MT} {za} {M} "{mtname}, M={M}" /'
)

return card9_lines
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If we pass in mt_dict[MT] as mt_data and isomer_dict[MT] as isomer_data
and return a fresh set of card9_lines to be appended by the function that calls this one.

Suggested change
za = str((pKZA + mt_dict[MT]['delKZA']) // 10).zfill(6)
for MF in isomer_dict[MT].keys():
if len(isomer_dict[MT][MF]) > 1:
for M in isomer_dict[MT][MF]:
card9_lines.append(
f'-1 / \n {MF} {MT} {za} {M} "{mtname}, M={M}" /'
)
return card9_lines
card9_lines = []
za = str((pKZA + mt_data['delKZA']) // 10).zfill(6)
for MF in isomer_data:
for M in isomer_data[MF]:
card9_lines.append(
f'-1 / \n {MF} {MT} {za} {M} "{mtname}, M={M}" /'
)
return card9_lines

Comment on lines +300 to +311
elif mf == 3 and mt == 0:
# SEND record; end of MT section
if current_section_lines:
sections.append(current_section_lines)
current_section_lines = []
line_count = 0

else:
if current_section_lines:
sections.append(current_section_lines)
current_section_lines = []
line_count = 0
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Don't these two blocks do the same thing? Can we collapse them to one? or should thet be different?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good catch. I don't think they need to be separated like this. I'll put them both under the else

Copy link
Copy Markdown
Member

@gonuke gonuke left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A few more thoughts on the card9 clarity

]
elements = dict(zip(elements, range(1, len(elements)+1)))

def write_card9_special_isomer_reactions(
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

While we're here... maybe we can drop the write_ from the beginning here since it doesn't actually write anything.

In may cases, we like methods to be verbs, but this one is more like a "getter" function and could get away with just being card9_special_isomer_reactions()

Comment on lines +223 to +225
has_isomers = (
MT in isomer_dict and isomer_dict[MT].keys() != {MFD}
)
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What if we changed this to be

Suggested change
has_isomers = (
MT in isomer_dict and isomer_dict[MT].keys() != {MFD}
)
has_special_isomers = pKZA is not None and (
MT in isomer_dict and isomer_dict[MT].keys() != {MFD}
)

Copy link
Copy Markdown
Member

@gonuke gonuke left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I did a full read through everything according to the program flow so see how it all fits together now that we're reaching the end and found a set of small polishing suggestions.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This function is only used once and only the values() of this dict are used. Maybe it can just return a list of dictionaries instead of a dictionary of dictionaries? Does it need to assemble this as a dictionary or could it just assemble a list internally? (I could imagine that the dictionary structure could provide some value vis-a-vis ensuring only one entry per ZA... but maybe that's not a real problem?)

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I could restructure this to return a list and just print out a warning in case of a duplicate ZA, which would probably be useful to have anyways.


material_id, MTs, endftk_file_obj = tp.extract_endf_specs(TAPE20)
material_id, MTs = tp.extract_endf_specs(TAPE20)
MTs = set(MTs).intersection(mt_dict.keys())
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should we issue a warning if there are MTs in the file that we ignore because of this?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure that we should, at least with the way that mt_dict is defined as is. When we construct mt_dict in reaction_data.py, it already filters out certain MTs that are not relevant, such as MTs 1-3. That said, I can take a closer look to for other MTs that would always be removed and then check for differences beyond them, although I'm not expecting to see so many mismatches in that case because I doubt that we would be encountering MTs in any authentic TENDL files that are not present in the ENDF-6 Manual, Appendix B, where we get the MT data from in the first place.

Comment on lines +120 to +122
MTs |= {
gas_MT for gas_MT in set(pendf_MTs) & set(rxd.GAS_DF['total_mt'])
}
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For readability and knowledge transfer, I think it might be good to separate out this line and provide some semantic meaning.

  1. I'm not sure where the logical & happens in the order of operations
  2. I'm not sure, semantically, what this new set represents

Is this what's happening??

Suggested change
MTs |= {
gas_MT for gas_MT in set(pendf_MTs) & set(rxd.GAS_DF['total_mt'])
}
known_gas_MTs = set(pendf_MTs)
known_gas_MTs &= set(rxd.GAS_DF['total_mt'])
MTs |= known_gas_MTs

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, your interpretation is correct, although I don't know if I would want to make the initial assignment of known_gas_MTs to be the full set(pendf_MTs), given that includes all present MTs, including but not limited to the gases.

How would this look for some greater clarity, while avoiding a premature assignment of known_gas_MTs:

known_gas_MTs = set(pendf_MTs).intersection(set(rxd.GAS_DF['total_mt']))
MTs |= known_gas_MTs

Do you think that would be reasonably clear enough?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's fine, but then maybe it's as simple as:

known_gas_MTs = set(pendf_MTs) & set(rxd.GAS_DF['total_mt'])

Copy link
Copy Markdown
Member

@gonuke gonuke Mar 27, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Any maybe it doesn't actually need to be changed that much... maybe we're just realizing that the for loop in the list comprehension isn't necessary and this can just be:

MTs |= (set(pendf_MTs) & set(rxd.GAS_DF['total_mt']))

MTs |= {
gas_MT for gas_MT in set(pendf_MTs) & set(rxd.GAS_DF['total_mt'])
}
pKZA = tp.extract_pendf_pkza(pendf_path)
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm still a little confused about why we don't already know the pKZA? This method creates a PENDF file for a particular combination of element and A - doesn't that define the pKZA without having to consult the file?

It appears to be checking to see if the parent is an isomer, but shouldn't we know that already? Isn't that part of the TENDL file name that got us here in the first place?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's a good point. We really shouldn't need to be opening the PENDF file at all because once the main loop has extracted and defined element and A, we already have all of the necessary information to produce the pKZA independently of any NJOY-related steps. I'll modify this logic to be used in a new function that can be called right after element and A are defined.

if gendf_path:
# Extract MT values again from GENDF file as there may be some
# difference from the original MT values in the ENDF/PENDF files
MTs = tp.update_gendf_MTs(gendf_path)
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If this list of MTs doesn't match the original list, should we warn the user?

reaction pathways for the given parent and its MTs.
"""

for MT in MTs:
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since you only seem to do something when it's NOT in EXCITATION_REACTIONS, maybe you can do some set operation like this and avoid the if statement?
(syntax untested)

Suggested change
for MT in MTs:
for MT in (MTs - EXCITATION_REACTIONS):

Returns:
is_ground_state (bool): True if the nuclide is in its ground state,
False if excited.
sigma_dict (dict): Dictionary keyed by tuples of (MT, M) with values
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The only use of this method never uses the MT number in the tuple that is the dictionary key - could we just use only M as the key? Or return two variables being M and the cross-section?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I suppose that's true because it's already being looped over for a single MT when it is called.

Comment on lines +434 to +437
dKZA = (
GAS_DF.loc[GAS_DF['gas'] == gas, 'kza'].iat[0] if gas
else (((pKZA // 10) * 10 + rxn['delKZA']) // 10) * 10 + M
)
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

readability - I'm realizing these inline if statements should probably be used more sparingly? (maybe it's just me....)

Suggested change
dKZA = (
GAS_DF.loc[GAS_DF['gas'] == gas, 'kza'].iat[0] if gas
else (((pKZA // 10) * 10 + rxn['delKZA']) // 10) * 10 + M
)
dKZA = (((pKZA // 10) * 10 + rxn['delKZA']) // 10) * 10 + M
if gas:
dKZA = GAS_DF.loc[GAS_DF['gas'] == gas, 'kza'].iat[0]

Comment on lines +368 to +371
for lowered_M in range(min(M - 1, 9), 0, -1):
trial_KZA = dKZA + lowered_M
if trial_KZA in eaf_nucs:
return trial_KZA
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
for lowered_M in range(min(M - 1, 9), 0, -1):
trial_KZA = dKZA + lowered_M
if trial_KZA in eaf_nucs:
return trial_KZA
trial_M = min(M-1, 9)
while (dKZA + trial_M) not in eaf_nucs and trial_M>0:
trial_M -= 1
return dKZA + trial_M

M, dKZA, eaf_nucs
)

if dKZA:
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Based on the original implementation of incrementally_deexcite_isomer() I think (?) dKZA will be None if and only if none of the isomers above ground state are in eaf_nucs. If that's true, does that mean we've decided to not include this cross-section at all, rather than including it in the ground state or our special MT?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I believe that your understanding is correct here. To fix this, following the while loop suggestion that you offered above for incrementally_deexcite_isomer(), I've modified it slightly to go all the way down to ground, rather than the first excited state (it will still be saved with a "special MT" even as it goes to ground).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants