Skip to content

Fix spline generation issue involving trace isotopes of elements that use a spectral function nuclear model #2

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
Oct 31, 2025
+105 −31
Merged

Fix spline generation issue involving trace isotopes of elements that use a spectral function nuclear model #2

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
7a9b620
Adjust the SpectralFunc class so that a missing spectral function can be
sjgardiner Oct 30, 2025
d577187
Update destructor to protect against double free
kjplows Oct 30, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 5 additions & 5 deletions config/AR25_20i/ModelConfiguration.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,11 +39,11 @@ University of Liverpool

<!--<param type="alg" name="NuclearModel"> genie::FGMBodekRitchie/Default </param>-->

<!-- For 12C, 16O, and 40Ar, use the spectral function -->
<param type="alg" name="NuclearModel@Pdg=1000060120"> genie::SpectralFunc/Default </param>
<!--<param type="alg" name="NuclearModel@Pdg=1000260560"> genie::SpectralFunc/Default </param>-->
<param type="alg" name="NuclearModel@Pdg=1000080160"> genie::SpectralFunc/Default </param>
<param type="alg" name="NuclearModel@Pdg=1000180400"> genie::SpectralFunc/Default </param>
<!-- For C, O, and Ar targets, use the spectral function -->
<param type="alg" name="NuclearModel@Pdg=1000060120"> genie::SpectralFunc/ApproxElements </param>
<!--<param type="alg" name="NuclearModel@Pdg=1000260560"> genie::SpectralFunc/ApproxElements </param>-->
<param type="alg" name="NuclearModel@Pdg=1000080160"> genie::SpectralFunc/ApproxElements </param>
<param type="alg" name="NuclearModel@Pdg=1000180400"> genie::SpectralFunc/ApproxElements </param>

<!--
Example of specific model for specific nuclei
Expand Down
10 changes: 5 additions & 5 deletions config/AR25_22i/ModelConfiguration.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,11 +39,11 @@ University of Liverpool

<!--<param type="alg" name="NuclearModel"> genie::FGMBodekRitchie/Default </param>-->

<!-- For 12C, 16O, and 40Ar, use the spectral function -->
<param type="alg" name="NuclearModel@Pdg=1000060120"> genie::SpectralFunc/Default </param>
<!--<param type="alg" name="NuclearModel@Pdg=1000260560"> genie::SpectralFunc/Default </param>-->
<param type="alg" name="NuclearModel@Pdg=1000080160"> genie::SpectralFunc/Default </param>
<param type="alg" name="NuclearModel@Pdg=1000180400"> genie::SpectralFunc/Default </param>
<!-- For C, O, and Ar targets, use the spectral function -->
<param type="alg" name="NuclearModel@Pdg=1000060120"> genie::SpectralFunc/ApproxElements </param>
<!--<param type="alg" name="NuclearModel@Pdg=1000260560"> genie::SpectralFunc/ApproxElements </param>-->
<param type="alg" name="NuclearModel@Pdg=1000080160"> genie::SpectralFunc/ApproxElements </param>
<param type="alg" name="NuclearModel@Pdg=1000180400"> genie::SpectralFunc/ApproxElements </param>

<!--
Example of specific model for specific nuclei
Expand Down
4 changes: 4 additions & 0 deletions config/SpectralFunc.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,4 +34,8 @@ FermiMomentumTable string No Table of Fermi momentum (kF) constan

</param_set>

<param_set name="ApproxElements">
<param type="bool" name="AllowElementMatch"> true </param>
</param_set>

</alg_conf>
107 changes: 86 additions & 21 deletions src/Physics/NuclearState/SpectralFunc.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,19 +32,9 @@
#include "Framework/Messenger/Messenger.h"
#include "Physics/NuclearState/SpectralFunc.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/ParticleData/PDGUtils.h"
#include "Framework/Numerical/RandomGen.h"

namespace {

// TODO: Replace this with std::to_string when we switch to C++11
std::string replace_with_std_to_string(int an_integer) {
std::ostringstream oss;
oss << an_integer;
return oss.str();
}

}

using namespace genie;
using namespace genie::constants;
using namespace genie::controls;
Expand All @@ -63,13 +53,7 @@ SpectralFunc::SpectralFunc(string config) :
//____________________________________________________________________________
SpectralFunc::~SpectralFunc()
{
// Delete the TH2D objects from the spectral function map
std::map<std::pair<int,int>, TH2D*>::iterator begin = fSpectralFunctionMap.begin();
std::map<std::pair<int,int>, TH2D*>::iterator end = fSpectralFunctionMap.end();
for (std::map<std::pair<int,int>, TH2D*>::iterator iter = begin; iter != end; ++iter) {
TH2D* hist = iter->second;
if ( hist ) delete hist;
}
fSpectralFunctionMap.clear();
}
//____________________________________________________________________________
bool SpectralFunc::GenerateNucleon(const Target& target) const
Expand Down Expand Up @@ -198,6 +182,8 @@ void SpectralFunc::LoadConfig(void)
}

fDataPath = data_path;

this->GetParamDef( "AllowElementMatch", fAllowElementMatch, false );
}
//____________________________________________________________________________
TH2D* SpectralFunc::SelectSpectralFunction(const Target& t) const
Expand All @@ -210,16 +196,95 @@ TH2D* SpectralFunc::SelectSpectralFunction(const Target& t) const
if ( my_iter != fSpectralFunctionMap.end() ) return my_iter->second;

// If not, attempt to build it
std::string target_pdg_string = replace_with_std_to_string( target_pdg );
std::string hitnuc_pdg_string = replace_with_std_to_string( hitnuc_pdg );
std::string target_pdg_string = std::to_string( target_pdg );
std::string hitnuc_pdg_string = std::to_string( hitnuc_pdg );

RgKey sf_key( "SpectFuncTable@Pdg=" + target_pdg_string + "_" + hitnuc_pdg_string );
std::string data_filename;
this->GetParamDef( sf_key, data_filename, std::string() );

if ( data_filename.empty() ) {
LOG("SpectralFunc", pERROR) << "** The spectral function for target "
LOG("SpectralFunc", pNOTICE) << "** The spectral function for target "
<< target_pdg << " and hit nuc PDG " << hitnuc_pdg << " isn't available";

// If element matching is enabled, then check for a spectral function with
// the same proton number Z as the requested nuclide.
if ( fAllowElementMatch ) {

// Get the proton number for the requested target nuclide
int tZ = genie::pdg::IonPdgCodeToZ( target_pdg );

// Not all configured spectral functions may have been built yet, so rely
// on the owned Registry to see if a suitable definition exists in the
// Algorithm configuration
const auto& conf_map = this->GetConfig().GetItemMap();

// Search the configuration for a matching proton number Z and matching
// hit nucleon PDG code. Store the nuclear PDG code for the result if
// we find a match.
int element_match_pdg;
auto iter = std::find_if( conf_map.cbegin(), conf_map.cend(),
[ tZ, hitnuc_pdg, &element_match_pdg ]( const auto& conf_pair ) -> bool {
const std::string prefix = "SpectFuncTable@Pdg=";
const std::size_t prefix_size = prefix.size();

// Check that the current config key begins with the exact prefix above
const auto& key = conf_pair.first;
// If not, then it cannot be a match
if ( key.compare(0, prefix_size, prefix) != 0 ) return false;

// Find an underscore following the prefix
size_t underscore_pos = key.find( '_', prefix_size );
// If there is no underscore, then the current key is not a match
if ( underscore_pos == std::string::npos ) return false;

// Extract the string between the prefix and the underscore. Assume
// it is a string representation of a nuclear PDG code
std::string tgt_pdg_str = key.substr( prefix_size,
underscore_pos - prefix_size );

// Also extract the remainder of the string following the underscore.
// Assume it is a string representation of a nucleon PDG code
std::string nucleon_pdg_str = key.substr( underscore_pos + 1u );

// Convert the string representations to integers
int tgt_pdg = std::stoi( tgt_pdg_str );
int nucleon_pdg = std::stoi( nucleon_pdg_str );

// Extract the proton number for the current nuclide
int target_Z = genie::pdg::IonPdgCodeToZ( tgt_pdg );

if ( hitnuc_pdg == nucleon_pdg && tZ == target_Z ) {
element_match_pdg = tgt_pdg;
return true;
}
return false;
}
);

// If a match was found, use it to approximate the spectral function for
// the current nuclide of interest
if ( iter != conf_map.end() ) {
LOG("SpectralFunc", pNOTICE) << "Approximating the spectral function"
<< " for target " << target_pdg << " with the one for target "
<< element_match_pdg;

// Use recursion to retrieve or build the spectral function for the
// nuclide with the matching proton number
genie::Target temp_tgt( element_match_pdg );
temp_tgt.SetHitNucPdg( hitnuc_pdg );
TH2D* sf_ptr = this->SelectSpectralFunction( temp_tgt );

// Add an entry for the current nuclide to the map so we don't have to do
// the Z-based search again
fSpectralFunctionMap[ std::make_pair( hitnuc_pdg, target_pdg ) ] = sf_ptr;
return sf_ptr;
}
}

// Fail with an error if the user requests a spectral function that
// is not available
LOG("SpectralFunc", pERROR) << "Failed to load spectral function data";
std::exit( 1 );
}

Expand Down
5 changes: 5 additions & 0 deletions src/Physics/NuclearState/SpectralFunc.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,6 +81,11 @@ class SpectralFunc : public NuclearModelI {
/// The number of removal energy bins to use when making spectral function
/// histograms
int fNumYBins;

/// If this flag is true, then a cached spectral function with the same
/// proton number Z will be used when a perfect match is not available
/// for the requested nuclide.
bool fAllowElementMatch = false;
};

} // genie namespace
Expand Down