diff --git a/config/Messenger.xml b/config/Messenger.xml index 713a15a71f..64268a6aa4 100644 --- a/config/Messenger.xml +++ b/config/Messenger.xml @@ -205,6 +205,7 @@ NOTICE NOTICE INFO + WARN INFO INFO WARN diff --git a/config/SplinePostProcessor.xml b/config/SplinePostProcessor.xml new file mode 100644 index 0000000000..499e2d5032 --- /dev/null +++ b/config/SplinePostProcessor.xml @@ -0,0 +1,35 @@ + + + + + + + + + genie::NievesQELCCPXSec/ZExp + 0.1,0.2,0.5,1.0,2.0,5.0,10.0,50.0,100.0 + 0.005,0.01,0.025,0.05,0.1,0.25,0.47,2.34,4.7 + + 2.8,3.0,4.0,5.0,10.0 + 0.005,0.01,0.03,0.075,0.1 + + true + + + + diff --git a/config/master_config.xml b/config/master_config.xml index 8bfe6aa53f..46a356270f 100644 --- a/config/master_config.xml +++ b/config/master_config.xml @@ -179,6 +179,9 @@ HEDISXSec.xml BostedChristyEMPXSec.xml + + SplinePostProcessor.xml + AhrensNCELPXSec.xml AhrensDMELPXSec.xml diff --git a/src/Apps/gMakeSplines.cxx b/src/Apps/gMakeSplines.cxx index 8685b64324..5a88ef4016 100644 --- a/src/Apps/gMakeSplines.cxx +++ b/src/Apps/gMakeSplines.cxx @@ -373,7 +373,7 @@ PDGCodeList * GetNeutrinoCodes(void) vector::const_iterator iter; for(iter = nuvec.begin(); iter != nuvec.end(); ++iter) { list->push_back( atoi(iter->c_str()) ); - } + } return list; } //____________________________________________________________________________ diff --git a/src/Framework/Numerical/LinkDef.h b/src/Framework/Numerical/LinkDef.h index d03b013a39..c9c239883c 100644 --- a/src/Framework/Numerical/LinkDef.h +++ b/src/Framework/Numerical/LinkDef.h @@ -5,11 +5,15 @@ #pragma link off all functions; #pragma link C++ namespace genie; +#pragma link C++ namespace genie::exceptions; #pragma link C++ namespace genie::utils::math; #pragma link C++ namespace genie::utils::gsl; #pragma link C++ class genie::RandomGen; #pragma link C++ class genie::Spline; +#pragma link C++ class genie::SplinePostProcessorI; +#pragma link C++ class genie::exceptions::SplineProcessingException; +#pragma link C++ class genie::SplinePostProcessor; #pragma link C++ class genie::BLI2DGrid; #pragma link C++ class genie::BLI2DUnifGrid; #pragma link C++ class genie::BLI2DNonUnifGrid; diff --git a/src/Framework/Numerical/SplinePostProcessor.cxx b/src/Framework/Numerical/SplinePostProcessor.cxx new file mode 100644 index 0000000000..aade3e6c33 --- /dev/null +++ b/src/Framework/Numerical/SplinePostProcessor.cxx @@ -0,0 +1,333 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2025, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + John Plows + University of Liverpool +*/ +//____________________________________________________________________________ + +#include "Framework/Numerical/SplinePostProcessor.h" +#include "Framework/Messenger/Messenger.h" + +using std::string; +using std::ostream; +using std::vector; + +using namespace genie; + +//___________________________________________________________________________ +namespace genie::exceptions { + ostream & operator << (ostream & stream, + const SplineProcessingException & exception ) { + exception.Print(stream); + return stream; + } +//___________________________________________________________________________ + SplineProcessingException::SplineProcessingException() { + this->Init(); + } +//___________________________________________________________________________ + SplineProcessingException::SplineProcessingException( const SplineProcessingException & exception ) { + this->Copy(exception); + } +//___________________________________________________________________________ + SplineProcessingException::SplineProcessingException( string reason, bool interrupt ) { + this->Init(); + this->SetReason(reason); + this->SetInterrupt(interrupt); + } +//___________________________________________________________________________ + SplineProcessingException::~SplineProcessingException() { + + } +//___________________________________________________________________________ + void SplineProcessingException::Init(void) { + fReason = ""; + fInterrupt = false; + } +//___________________________________________________________________________ + void SplineProcessingException::Copy(const SplineProcessingException & exception) { + fReason = exception.fReason; + fInterrupt = exception.fInterrupt; + } +//___________________________________________________________________________ + void SplineProcessingException::Print(ostream & stream) const { + stream << "** EXCEPTION in SplinePostProcessor! Reason : " << this->ShowReason(); + if( this->DoInterrupt() ) { + stream << " *~*~*~* INTERRUPTING EXECUTION!! *~*~*~* \n"; + } + } +//___________________________________________________________________________ +} // namespace exceptions + +using namespace genie::exceptions; + +//___________________________________________________________________________ +SplinePostProcessor * SplinePostProcessor::fInstance = 0; +//___________________________________________________________________________ +SplinePostProcessor::SplinePostProcessor() : SplinePostProcessorI("genie::SplinePostProcessor") { + this->Init(); +} +//___________________________________________________________________________ +SplinePostProcessor::SplinePostProcessor(string config) : + SplinePostProcessorI("genie::SplinePostProcessor", config) { + this->Init(); +} +//___________________________________________________________________________ +SplinePostProcessor::SplinePostProcessor(string name, string config) : + SplinePostProcessorI(name, config) { + this->Init(); +} +//___________________________________________________________________________ +SplinePostProcessor::~SplinePostProcessor() { + fStops.clear(); + fWidths.clear(); + fStopMap.clear(); + fWidthMap.clear(); + fAlgs.clear(); +} +//___________________________________________________________________________ +SplinePostProcessor * SplinePostProcessor::Instance() +{ + if(fInstance == 0) fInstance = new SplinePostProcessor; + return fInstance; +} +//___________________________________________________________________________ +void SplinePostProcessor::Init(void) { + if(fAlgs.size() > 0) return; + LOG("SplinePostProcessor", pDEBUG) << "Initialising SplinePostProcessor..."; + this->Configure("Default"); +} +//___________________________________________________________________________ +void SplinePostProcessor::LoadConfig(void) { + fStops.clear(); fWidths.clear(); fAlgs.clear(); fStopMap.clear(); fWidthMap.clear(); + LOG("SplinePostProcessor", pDEBUG) << "Loading the SplinePostProcessor configuration..."; + + vector alg_vect; + this->GetParamVect( "HandledAlgorithms", alg_vect ); + // Drop these into the set for quick lookup + for( string alg : alg_vect ) { + fAlgs.insert(alg); + LOG("SplinePostProcessor", pDEBUG) << "Will handle algorithm named " << alg; + } + + this->GetParamVect( "Stops", fStops ); + this->GetParamVect( "Widths", fWidths ); + fStopMap[0] = fStops; + fWidthMap[0] = fWidths; + + // If there are any more keys, typically delineated with "@Pdg=...", add these to the map + // Similar architecture to NuclearModelMap's implementation + RgIMap entries = (this->GetConfig()).GetItemMap(); + for(RgIMap::const_iterator it = entries.begin(); it != entries.end(); ++it){ + const std::string key = it->first; // Notice this prepends "N" and appends "s" + const std::string sens = "@Pdg="; + if( key.find(sens.c_str()) != std::string::npos + && strcmp(key.substr(0, 1).c_str(), "N") == 0 ) { + // Figure out the size of this thing + int n; + LOG("SplinePostProcessor", pDEBUG) << "Getting N stops from here " << key; + this->GetParam( key, n ); + + // Decide if this is stops or widths + // See https://stackoverflow.com/questions/14265581 + std::string comm_name(key); + comm_name.erase(0, 1); + comm_name.erase(comm_name.size() - 1); + std::string dup_name(comm_name); dup_name.append(","); + LOG("SplinePostProcessor", pDEBUG) << "Key = " << comm_name << " with size " + << comm_name.size(); + if( key.find("Stops") != std::string::npos ) { + std::string subs = dup_name.substr(10); // 5 "Stops" + 5 "@Pdg=" + std::vector pdgs; + size_t pos = 0; + while( (pos = subs.find(",", 0)) != std::string::npos ) { + int pdg = std::stoi(subs.substr(0, pos)); + LOG("SplinePostProcessor", pDEBUG) << "Adding key " << pdg << " to stops"; + subs.erase(0, pos+1); // "," has length 1 + std::vector stops; stops.resize(n); + for( int i = 0; i < n; ++i ) { + RgKey tmp_key = Algorithm::BuildParamVectKey(comm_name, i); + this->GetParam( tmp_key, stops[i] ); + } + fStopMap[pdg] = stops; + } // split strings + + } else if ( key.find("Widths") != std::string::npos ) { + std::string subs = dup_name.substr(11); // 6 "Widths" + 5 "@Pdg=" + std::vector pdgs; + size_t pos = 0; + while( (pos = subs.find(",", 0)) != std::string::npos ) { + int pdg = std::stoi(subs.substr(0, pos)); + LOG("SplinePostProcessor", pDEBUG) << "Adding key " << pdg << " to widths"; + subs.erase(0, pos+1); // "," has length 1 + std::vector widths; widths.resize(n); + for( int i = 0; i < n; ++i ) { + RgKey tmp_key = Algorithm::BuildParamVectKey(comm_name, i); + this->GetParam( tmp_key, widths[i] ); + } + fWidthMap[pdg] = widths; + } // split strings + } // stops or widths + } // extract PDG entries + } // loop over registry entries of this Alg + + // Check for misconfigurations + for ( const auto & mapkey : fStopMap ) { + if( fStopMap[mapkey.first].size() != fWidthMap[mapkey.first].size() ) { + // The Algorithm has been misconfigured. It's early, so complain and ask user to fix. + LOG("SplinePostProcessor", pFATAL) << "Different number of stops and widths! " + << " (" << fStopMap[mapkey.first].size() + << " vs " << fWidthMap[mapkey.first].size() << ")" + << " at PDG key = " << mapkey.first + << " (key 0 means default)"; + throw SplineProcessingException("N stops != N widths", true); + } // misconfiguration check + } // over each key + + this->GetParam( "UsePostProcessing", fUsePostProcessing ); + +} +//___________________________________________________________________________ +void SplinePostProcessor::Configure(const Registry & config) { + Algorithm::Configure(config); + try { + this->LoadConfig(); + } catch (SplineProcessingException exception) { + LOG("SplinePostProcessor", pFATAL) << exception; + exit(1); + } +} +//___________________________________________________________________________ +void SplinePostProcessor::Configure(std::string param_set) { + Algorithm::Configure(param_set); + try { + this->LoadConfig(); + } catch (SplineProcessingException exception) { + LOG("SplinePostProcessor", pFATAL) << exception; + exit(1); + } +} +//___________________________________________________________________________ +bool SplinePostProcessor::IsHandled(std::string alg_name, std::string alg_config) { + std::string alg_key = alg_name + "/" + alg_config; + return (fAlgs.count(alg_key) != 0); +} +//___________________________________________________________________________ +bool SplinePostProcessor::IsHandled(const Algorithm * alg) { + std::string alg_key = alg->Id().Key(); + return (fAlgs.count(alg_key) != 0); +} +//___________________________________________________________________________ +std::vector SplinePostProcessor::ProcessSpline( const std::vector E, + const std::vector xsec, + const int pdg ) const { + + // If we don't want post processing, return out now. + if( !fUsePostProcessing ) { + return xsec; + } // no op if no post processing + + // Set the stops and widths to be what we asked for + int elem = (fStopMap.find(pdg) != fStopMap.end()) ? pdg : 0; + this->SetStops(fStopMap.at(elem)); this->SetWidths(fWidthMap.at(elem)); + + // Are the thresholds sensible? + std::vector stops; + std::vector widths; + try { + this->SetThresholds( stops, widths, E.front(), E.back() ); + } catch (SplineProcessingException exception) { + if(exception.DoInterrupt()) { + LOG("SplinePostProcessor", pFATAL) << exception; + exit(1); + } else { + LOG("SplinePostProcessor", pWARN) << exception; + } + } // Set thresholds + + // Read out the knots from the spline and put them into arrays + const int NKnots = E.size(); + std::vector new_xsec; + + // Apply Gaussian kernel smoother here, with a linear activation function for widths + std::vector::iterator it_stop = stops.begin()+1; + std::vector::iterator it_width = widths.begin()+1; + + double w1 = widths.front(); double w2 = *it_width; + for( int i = 0; i < NKnots; i++ ) { + // Linear activation function between s1 and s2 + auto activation = [&](const double s1, const double s2) { + if(E[i] <= s1) return 0.0; + if(E[i] >= s2) return 1.0; + return (E[i] - s1)/(s2 - s1); + }; + + if(E[i] > *it_stop && E[i] < stops.back()) { // Update iterators + it_stop++; it_width++; + w1 = w2; w2 = *(it_width); + } + + double rat = activation(*(it_stop-1), *(it_stop)); + double width = w1 + (w2 - w1) * rat; // between w1 and w2 + double num = 0.0; double den = 0.0; + for(int j = 0; j < NKnots; j++) { + double w = std::exp(-0.5 * std::pow((E[i] - E[j]) / width, 2.0)); // Gauss + num += w * xsec[j]; // weight * value + den += w; // total weight + } + + // Knot value out + new_xsec.emplace_back( num / den ); + + } // Loop over knots of spline; smooth each one by configurable width + + return new_xsec; +} +//___________________________________________________________________________ +void SplinePostProcessor::SetThresholds( std::vector & stops, + std::vector & widths, + double xmin, double xmax ) const { + // Copy fStops and fWidths to the temp objects, and then start working on them + stops = fStops; widths = fWidths; + if( stops.front() > xmax || stops.back() < xmin ) { // Bad configuration, this will be nonsense. + throw SplineProcessingException("Thresholds out of bounds", true); + return; + } + + int n_rejected_front = 0; int n_rejected_back = 0; + while( stops.front() < xmin ) { // Throw out elements until we're in the range + stops.erase(stops.begin()); + widths.erase(widths.begin()); + n_rejected_front++; + } + while( stops.back() > xmax ) { // Throw out elements until we're in the range + stops.erase(stops.end()-1); + widths.erase(widths.end()-1); + n_rejected_back++; + } + + if( stops.size() == 1 ) { // Too few elements. Stop execution + throw SplineProcessingException("Not enough elements left to interpolate", true); + } + + // If the xmin and/or xmax are outside first / last stop, add some with copies of width + if( xmin < stops.front() ) { + stops.insert(stops.begin(), xmin); + widths.insert(widths.begin(), widths.front()); + } + if( xmax > stops.back() ) { + stops.emplace_back(xmax); + widths.emplace_back(widths.back()); + } + + if(n_rejected_front > 0 || n_rejected_back > 0) { // Complain, but don't stop execution + LOG("SplinePostProcessor", pWARN) << "Removed " << n_rejected_front << " elements from the front" + << " of the stops / widths vectors, and " << n_rejected_back + << " from the back. Take note!"; + throw SplineProcessingException("Using different stops to those configured"); + } + return; +} diff --git a/src/Framework/Numerical/SplinePostProcessor.h b/src/Framework/Numerical/SplinePostProcessor.h new file mode 100644 index 0000000000..ca4df96f9c --- /dev/null +++ b/src/Framework/Numerical/SplinePostProcessor.h @@ -0,0 +1,136 @@ +//____________________________________________________________________________ +/*! + +\class genie::SplinePostProcessor + +\brief Post processor class to smooth out a Spline object. + This acts on the xsec values (not the abscissae, E) of each Spline + by implementing a Gaussian kernel of variable width. + + The width of this kernel (in GeV) is user-configurable, and is linearly interpolated + between stops. + See en.wikipedia.org/wiki/Kernel_smoother#Gaussian_kernel_smoother + +\author John Plows + University of Liverpool + +\created November 3, 2025 + +\cpright Copyright (c) 2003-2025, The GENIE Collaboration + For the full text of the license visit + http://genie-mc.github.io/copyright.html +*/ +//____________________________________________________________________________ + +#ifndef _SPLINE_POST_PROCESSOR_H_ +#define _SPLINE_POST_PROCESSOR_H_ + +#include "Framework/Numerical/SplinePostProcessorI.h" +#include "Framework/Algorithm/AlgConfigPool.h" +#include "Framework/Algorithm/AlgId.h" + +#include +#include +#include +#include +#include +#include + +namespace genie { + + // Add an exception handler for spline post processing + // make it part of genie::exceptions + namespace exceptions { + class SplineProcessingException { + + public: + SplineProcessingException(); + SplineProcessingException(const SplineProcessingException & exception); + SplineProcessingException(std::string reason, bool interrupt = false); + ~SplineProcessingException(); + + void Init(void); + void Copy(const SplineProcessingException & exception); + void Print(std::ostream & stream) const; + + void SetReason(std::string reason) { fReason = reason; } + void SetInterrupt(bool interrupt) { fInterrupt = interrupt; } + string ShowReason (void) const { return fReason; } + bool DoInterrupt (void) const { return fInterrupt; } + + friend std::ostream & operator << ( std::ostream & stream, + const SplineProcessingException & exception ); + + private: + + std::string fReason; + bool fInterrupt; + + }; // class genie::exceptions::SplineProcessingException + } // namespace exceptions + + class Spline; + + class SplinePostProcessor: public SplinePostProcessorI { + + public: + SplinePostProcessor(); + SplinePostProcessor(std::string name); + SplinePostProcessor(std::string name, std::string config); + ~SplinePostProcessor(); + + static SplinePostProcessor * Instance(); + + //-- implement the SplinePostProcessorI interface + std::vector ProcessSpline(const std::vector E, + const std::vector xsec, + const int pdg) const; + + //-- override the Algorithm::Configure methods to load configuration + // data to private data members + void Configure (const Registry & config); + void Configure (std::string param_set); + + //-- check if an algorithm is handled + bool IsHandled (std::string alg_name, std::string alg_config); + bool IsHandled (const Algorithm * alg); + + //-- Check if we've enabled post-processing + bool UsePostProcessing() const { return fUsePostProcessing; } + + //-- Simple getter-setter machinery + std::vector GetStops() const { return fStops; } + void SetStops(std::vector stops ) const { fStops = stops; } + std::vector GetWidths() const { return fWidths; } + void SetWidths(std::vector widths) const { fWidths = widths; } + + private: + + static SplinePostProcessor * fInstance; + + void Init (void); + void LoadConfig (void); + + void SetThresholds( std::vector & stops, std::vector & widths, + double xmin, double xmax ) const; + + //-- private data members + std::unordered_set fAlgs; ///< Algorithms handled by the SplinePostProcessor. + + // Stops and widths to be used when processing a spline + mutable std::vector fStops; ///< Abscissae for the interpolation of the Gaussian kernel smoother + mutable std::vector fWidths; ///< Widths of the Gaussian kernel smoother at the stops + + // Map of the stops and widths to be used, potentially. + // If a key with "@PDG=.." is seen, add this PDG to the map. + // There is always a key "0" that is the default. + std::unordered_map> fStopMap; + std::unordered_map> fWidthMap; + + bool fUsePostProcessing; ///< If true, do post processing. If false, do nothing. + + }; // class SplinePostProcessor + +} // namespace genie + +#endif // #ifndef _SPLINE_POST_PROCESSOR_H_ diff --git a/src/Framework/Numerical/SplinePostProcessorI.cxx b/src/Framework/Numerical/SplinePostProcessorI.cxx new file mode 100644 index 0000000000..e3b0dee894 --- /dev/null +++ b/src/Framework/Numerical/SplinePostProcessorI.cxx @@ -0,0 +1,38 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2025, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + John Plows + University of Liverpool +*/ +//____________________________________________________________________________ + +#include "Framework/Numerical/SplinePostProcessorI.h" + +using namespace genie; + +//____________________________________________________________________________ +SplinePostProcessorI::SplinePostProcessorI() : + Algorithm() +{ + +} +//____________________________________________________________________________ +SplinePostProcessorI::SplinePostProcessorI(string name) : + Algorithm(name) +{ + +} +//____________________________________________________________________________ +SplinePostProcessorI::SplinePostProcessorI(string name, string config) : + Algorithm(name, config) +{ + +} +//____________________________________________________________________________ +SplinePostProcessorI::~SplinePostProcessorI() +{ + +} +//____________________________________________________________________________ diff --git a/src/Framework/Numerical/SplinePostProcessorI.h b/src/Framework/Numerical/SplinePostProcessorI.h new file mode 100644 index 0000000000..11376ecb79 --- /dev/null +++ b/src/Framework/Numerical/SplinePostProcessorI.h @@ -0,0 +1,52 @@ +//____________________________________________________________________________ +/*! + +\class genie::SplinePostProcessorI + +\brief Interface for the SplinePostProcessor. + Concrete implementations of this interface use the 'Visitor' design + to perform an operation on a Spline. + +\author John Plows + University of Liverpool + +\created October 29, 2025 + +\cpright Copyright (c) 2003-2025, The GENIE Collaboration + For the full text of the license visit + http://genie-mc.github.io/copyright.html +*/ +//____________________________________________________________________________ + +#ifndef _SPLINE_POST_PROCESSOR_I_H_ +#define _SPLINE_POST_PROCESSOR_I_H_ + +#include "Framework/Algorithm/Algorithm.h" + +namespace genie { + + class Spline; + + class SplinePostProcessorI : public Algorithm { + + public: + + virtual ~SplinePostProcessorI(); + + //-- define the SplinePostProcessorI interface + + virtual std::vector ProcessSpline(const std::vector E, + const std::vector xsec, + const int pdg) const = 0; + + protected: + + SplinePostProcessorI(); + SplinePostProcessorI(string name); + SplinePostProcessorI(string name, string config); + + }; // class SplinePostProcessorI + +} // namespace genie + +#endif // #ifndef _SPLINE_POST_PROCESSOR_I_H_ diff --git a/src/Framework/Utils/XSecSplineList.cxx b/src/Framework/Utils/XSecSplineList.cxx index 66686cb267..ad5550a5fa 100644 --- a/src/Framework/Utils/XSecSplineList.cxx +++ b/src/Framework/Utils/XSecSplineList.cxx @@ -267,6 +267,24 @@ void XSecSplineList::CreateSpline(const XSecAlgorithmI * alg, } + fPostProcessor = *( SplinePostProcessor::Instance() ); + if( fPostProcessor.UsePostProcessing() ) { + // Check if need to post-process + if( fPostProcessor.IsHandled(alg) ) { + SLOG("XSecSplList", pINFO) << "Post processing spline with key: " << key; + xsec = fPostProcessor.ProcessSpline(E, xsec, + interaction->InitStatePtr()->ProbePdg()); + + // Output so the user can see the new spline + for( size_t i = 0; i < E.size(); i++ ) { + SLOG("XSecSplList", pNOTICE) + << "xsec(E = " << E[i] << ") = " + << (1E+38/units::cm2)*xsec[i] << " x 1E-38 cm^2 post-processing"; + } // output + } // if handled alg + } // if post processing + + // Warn about odd case of decreasing cross section // but allow for small variation due to integration errors const double eps_xsec = 1.0e-5; diff --git a/src/Framework/Utils/XSecSplineList.h b/src/Framework/Utils/XSecSplineList.h index 47604d4b67..3723832702 100644 --- a/src/Framework/Utils/XSecSplineList.h +++ b/src/Framework/Utils/XSecSplineList.h @@ -25,6 +25,7 @@ #include #include "Framework/Conventions/XmlParserStatus.h" +#include "Framework/Numerical/SplinePostProcessor.h" using std::map; using std::set; @@ -38,6 +39,7 @@ namespace genie { class XSecAlgorithmI; class Interaction; class Spline; +class SplinePostProcessor; class XSecSplineList; ostream & operator << (ostream & stream, const XSecSplineList & xsl); @@ -97,6 +99,8 @@ class XSecSplineList { virtual ~XSecSplineList(); static XSecSplineList * fInstance; + + SplinePostProcessor fPostProcessor; bool fUseLogE; int fNKnots;