From 576b9679505b03a03f5d312f3b6c485d26130fa3 Mon Sep 17 00:00:00 2001 From: kjplows Date: Mon, 2 Feb 2026 05:54:41 -0600 Subject: [PATCH 1/6] Switch off integral nuclear influence cutoff --- config/NievesQELCCPXSec.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/config/NievesQELCCPXSec.xml b/config/NievesQELCCPXSec.xml index 1a0be1658..43b08d802 100644 --- a/config/NievesQELCCPXSec.xml +++ b/config/NievesQELCCPXSec.xml @@ -40,7 +40,7 @@ RmaxMode string Yes Method to use to comput genie::NewQELXSec/Default genie::NuclearModelMap/Default UseNuclearModel - 10.0 + 1000000.0 true true From 21c08c477fccbd0cdb8de695952a6da9758f0ccb Mon Sep 17 00:00:00 2001 From: kjplows Date: Wed, 18 Feb 2026 08:14:53 -0600 Subject: [PATCH 2/6] Prevent unphysical events from crashing gevgen, part 2 --- src/Physics/QuasiElastic/EventGen/QELEventGenerator.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Physics/QuasiElastic/EventGen/QELEventGenerator.cxx b/src/Physics/QuasiElastic/EventGen/QELEventGenerator.cxx index a04dddbc2..44cc4655d 100644 --- a/src/Physics/QuasiElastic/EventGen/QELEventGenerator.cxx +++ b/src/Physics/QuasiElastic/EventGen/QELEventGenerator.cxx @@ -508,10 +508,10 @@ double QELEventGenerator::ComputeMaxXSec(const Interaction * in) const }// Done with centre-of-mass angles coarsely // Calculate the range for the next layer - costh_range_min = costh_at_xsec_max - costh_increment; - costh_range_max = costh_at_xsec_max + costh_increment; - phi_range_min = phi_at_xsec_max - phi_increment; - phi_range_max = phi_at_xsec_max + phi_increment; + costh_range_min = std::max(-1.0, costh_at_xsec_max - costh_increment); + costh_range_max = std::min( 1.0, costh_at_xsec_max + costh_increment); + phi_range_min = std::max(0.0, phi_at_xsec_max - phi_increment); + phi_range_max = std::max(2.0*TMath::Pi(), phi_at_xsec_max + phi_increment); double improvement_factor = this_nuc_xsec_max/last_layer_xsec_max; if (ilayer && (improvement_factor-1) < acceptable_fraction_of_safety_factor * (fSafetyFactor-1)) { From a516c47992205defc09679598c7af6b02a78d313 Mon Sep 17 00:00:00 2001 From: kjplows Date: Thu, 19 Feb 2026 06:56:42 -0600 Subject: [PATCH 3/6] Revert to explicit clamping for angular range --- src/Physics/QuasiElastic/EventGen/QELEventGenerator.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Physics/QuasiElastic/EventGen/QELEventGenerator.cxx b/src/Physics/QuasiElastic/EventGen/QELEventGenerator.cxx index 44cc4655d..2086ca71b 100644 --- a/src/Physics/QuasiElastic/EventGen/QELEventGenerator.cxx +++ b/src/Physics/QuasiElastic/EventGen/QELEventGenerator.cxx @@ -508,10 +508,10 @@ double QELEventGenerator::ComputeMaxXSec(const Interaction * in) const }// Done with centre-of-mass angles coarsely // Calculate the range for the next layer - costh_range_min = std::max(-1.0, costh_at_xsec_max - costh_increment); - costh_range_max = std::min( 1.0, costh_at_xsec_max + costh_increment); - phi_range_min = std::max(0.0, phi_at_xsec_max - phi_increment); - phi_range_max = std::max(2.0*TMath::Pi(), phi_at_xsec_max + phi_increment); + costh_range_min = costh_at_xsec_max - costh_increment; if(costh_range_min < -1.0) { costh_range_min = -1.0; } + costh_range_max = costh_at_xsec_max + costh_increment; if(costh_range_max > 1.0) { costh_range_max = 1.0; } + phi_range_min = phi_at_xsec_max - phi_increment; if( phi_range_min < -1.0) { phi_range_min = 0.0;} + phi_range_max = phi_at_xsec_max + phi_increment; if( phi_range_max > 2.0 * TMath::Pi()) { phi_range_max = 2.0 * TMath::Pi(); } double improvement_factor = this_nuc_xsec_max/last_layer_xsec_max; if (ilayer && (improvement_factor-1) < acceptable_fraction_of_safety_factor * (fSafetyFactor-1)) { From aa0af1ad7a62f913a1f66996976f9379eb860ec5 Mon Sep 17 00:00:00 2001 From: kjplows Date: Thu, 19 Feb 2026 07:00:08 -0600 Subject: [PATCH 4/6] Phi is never less than zero --- src/Physics/QuasiElastic/EventGen/QELEventGenerator.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Physics/QuasiElastic/EventGen/QELEventGenerator.cxx b/src/Physics/QuasiElastic/EventGen/QELEventGenerator.cxx index 2086ca71b..04f0bb24b 100644 --- a/src/Physics/QuasiElastic/EventGen/QELEventGenerator.cxx +++ b/src/Physics/QuasiElastic/EventGen/QELEventGenerator.cxx @@ -510,7 +510,7 @@ double QELEventGenerator::ComputeMaxXSec(const Interaction * in) const // Calculate the range for the next layer costh_range_min = costh_at_xsec_max - costh_increment; if(costh_range_min < -1.0) { costh_range_min = -1.0; } costh_range_max = costh_at_xsec_max + costh_increment; if(costh_range_max > 1.0) { costh_range_max = 1.0; } - phi_range_min = phi_at_xsec_max - phi_increment; if( phi_range_min < -1.0) { phi_range_min = 0.0;} + phi_range_min = phi_at_xsec_max - phi_increment; if( phi_range_min < 0.0) { phi_range_min = 0.0; } phi_range_max = phi_at_xsec_max + phi_increment; if( phi_range_max > 2.0 * TMath::Pi()) { phi_range_max = 2.0 * TMath::Pi(); } double improvement_factor = this_nuc_xsec_max/last_layer_xsec_max; From 55f715f290430fd7b46e19f895936cb5ad00c625 Mon Sep 17 00:00:00 2001 From: kjplows Date: Thu, 19 Feb 2026 10:06:10 -0600 Subject: [PATCH 5/6] Ad hoc fixes for Fermi motion plus exception handling in InitialState --- src/Framework/Interaction/InitialState.h | 1 + src/Physics/NuclearState/FermiMover.cxx | 18 +++++++++++++++++- src/Physics/NuclearState/FermiMover.h | 1 + src/Physics/NuclearState/NuclearModel.h | 1 + src/Physics/NuclearState/SpectralFunc.cxx | 4 ++-- 5 files changed, 22 insertions(+), 3 deletions(-) diff --git a/src/Framework/Interaction/InitialState.h b/src/Framework/Interaction/InitialState.h index ddc00f64b..91d407e9d 100644 --- a/src/Framework/Interaction/InitialState.h +++ b/src/Framework/Interaction/InitialState.h @@ -14,6 +14,7 @@ Other minor changes / additions and fixes were installed by: Andy Furmanski (Univ. of Manchester) Joe Johnston (Univ of Pittsburgh) + John Plows (Univ. of Liverpool) \created May 02, 2004 diff --git a/src/Physics/NuclearState/FermiMover.cxx b/src/Physics/NuclearState/FermiMover.cxx index 6d5020b3a..2110d8b13 100644 --- a/src/Physics/NuclearState/FermiMover.cxx +++ b/src/Physics/NuclearState/FermiMover.cxx @@ -145,7 +145,9 @@ void FermiMover::KickHitNucleon(GHepRecord * evrec) const FermiMoverInteractionType_t interaction_type = fNuclModel->GetFermiMoverInteractionType(); // EffectiveSF treatment or momentum-dependent removal energy - if (interaction_type == kFermiMoveEffectiveSF1p1h || fMomDepErmv ) { + if (interaction_type == kFermiMoveEffectiveSF1p1h || + interaction_type == kFermiMoveSpectralFunc || + fMomDepErmv ) { EN = nucleon->Mass() - w - pF2 / (2 * (nucleus->Mass() - nucleon->Mass())); } else if (interaction_type == kFermiMoveEffectiveSF2p2h_eject || interaction_type == kFermiMoveEffectiveSF2p2h_noeject) { @@ -197,6 +199,20 @@ void FermiMover::KickHitNucleon(GHepRecord * evrec) const nucleon->SetMomentum(*p4); // update GHEP value + // The below call to KPhaseSpace::IsAboveThreshold() depends on the target p4. + // If the EN < p3.Mag() then this is unphysical and kills code. Handle it here + if( p4->Mag2() < 0.0 ) { + LOG("FermiMover", pNOTICE) + << "Event generates an unphysical initial state target with p4 = " + << utils::print::P4AsShortString(p4) + << ". Skipping event"; + evrec->EventFlags()->SetBitNumber(kBelowThrNRF, true); + genie::exceptions::EVGThreadException exception; + exception.SetReason("Unphysical initial state 4-momentum after FermiMover acted on it"); + exception.SwitchOnFastForward(); + throw exception; + } + // Sometimes, for interactions near threshold, Fermi momentum might bring // the neutrino energy in the nucleon rest frame below threshold (for the // selected interaction). In this case mark the event as unphysical and diff --git a/src/Physics/NuclearState/FermiMover.h b/src/Physics/NuclearState/FermiMover.h index 7c1c94aa7..43846aae6 100644 --- a/src/Physics/NuclearState/FermiMover.h +++ b/src/Physics/NuclearState/FermiMover.h @@ -27,6 +27,7 @@ #include "Framework/Interaction/Target.h" #include "Physics/NuclearState/SRCNuclearRecoil.h" #include "Physics/NuclearState/SecondNucleonEmissionI.h" +#include "Framework/Utils/PrintUtils.h" namespace genie { diff --git a/src/Physics/NuclearState/NuclearModel.h b/src/Physics/NuclearState/NuclearModel.h index 33ac63429..3b805717a 100644 --- a/src/Physics/NuclearState/NuclearModel.h +++ b/src/Physics/NuclearState/NuclearModel.h @@ -40,6 +40,7 @@ typedef enum EFermiMoverInteractionType { kFermiMoveEffectiveSF1p1h, kFermiMoveEffectiveSF2p2h_eject, kFermiMoveEffectiveSF2p2h_noeject, + kFermiMoveSpectralFunc, } FermiMoverInteractionType_t; class NuclearModel { diff --git a/src/Physics/NuclearState/SpectralFunc.cxx b/src/Physics/NuclearState/SpectralFunc.cxx index 10d4395c4..8ae91480f 100644 --- a/src/Physics/NuclearState/SpectralFunc.cxx +++ b/src/Physics/NuclearState/SpectralFunc.cxx @@ -42,13 +42,13 @@ using namespace genie::controls; //____________________________________________________________________________ SpectralFunc::SpectralFunc() : NuclearModelI("genie::SpectralFunc") { - + fFermiMoverInteractionType = kFermiMoveSpectralFunc; } //____________________________________________________________________________ SpectralFunc::SpectralFunc(string config) : NuclearModelI("genie::SpectralFunc", config) { - + fFermiMoverInteractionType = kFermiMoveSpectralFunc; } //____________________________________________________________________________ SpectralFunc::~SpectralFunc() From 6e62ae274b285992a7c945d9ce7eaf7ba18b426e Mon Sep 17 00:00:00 2001 From: kjplows Date: Tue, 31 Mar 2026 10:12:49 -0500 Subject: [PATCH 6/6] Don't change nucleon energy in SpectralFunction --- src/Physics/NuclearState/FermiMover.cxx | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/Physics/NuclearState/FermiMover.cxx b/src/Physics/NuclearState/FermiMover.cxx index 2110d8b13..3cd0f70aa 100644 --- a/src/Physics/NuclearState/FermiMover.cxx +++ b/src/Physics/NuclearState/FermiMover.cxx @@ -145,9 +145,7 @@ void FermiMover::KickHitNucleon(GHepRecord * evrec) const FermiMoverInteractionType_t interaction_type = fNuclModel->GetFermiMoverInteractionType(); // EffectiveSF treatment or momentum-dependent removal energy - if (interaction_type == kFermiMoveEffectiveSF1p1h || - interaction_type == kFermiMoveSpectralFunc || - fMomDepErmv ) { + if (interaction_type == kFermiMoveEffectiveSF1p1h || fMomDepErmv ) { EN = nucleon->Mass() - w - pF2 / (2 * (nucleus->Mass() - nucleon->Mass())); } else if (interaction_type == kFermiMoveEffectiveSF2p2h_eject || interaction_type == kFermiMoveEffectiveSF2p2h_noeject) {