Skip to content

Better handling of SpectalFunction in the Fermi Mover, and switch off integral nuclear influence cutoff energy #14

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

Open
kjplows wants to merge 8 commits into
base: v3_06_02_br
Choose a base branch
from
Open

Better handling of SpectalFunction in the Fermi Mover, and switch off integral nuclear influence cutoff energy #14

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
576b967
Switch off integral nuclear influence cutoff
kjplows Feb 2, 2026
4c21280
Merge tag 'v3_06_02_sbn2p1' into bugfix/remove-unphysical-QEL-costh0
kjplows Feb 16, 2026
21c08c4
Prevent unphysical events from crashing gevgen, part 2
kjplows Feb 18, 2026
a516c47
Revert to explicit clamping for angular range
kjplows Feb 19, 2026
aa0af1a
Phi is never less than zero
kjplows Feb 19, 2026
55f715f
Ad hoc fixes for Fermi motion plus exception handling in InitialState
kjplows Feb 19, 2026
6e62ae2
Don't change nucleon energy in SpectralFunction
kjplows Mar 31, 2026
aa1ba7e
Merge branch 'v3_06_02_br' of https://github.com/SBNSoftware/Generato…
kjplows Apr 1, 2026
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
2 changes: 1 addition & 1 deletion config/NievesQELCCPXSec.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ RmaxMode string Yes Method to use to comput
<param type="alg" name="XSec-Integrator"> genie::NewQELXSec/Default </param>
<param type="alg" name="IntegralNuclearModel"> genie::NuclearModelMap/Default </param>
<param type="string" name="IntegralNucleonBindingMode"> UseNuclearModel </param>
<param type="double" name = "IntegralNuclearInfluenceCutoffEnergy"> 10.0 </param>
<param type="double" name = "IntegralNuclearInfluenceCutoffEnergy"> 1000000.0 </param>

<param type="bool" name="RPA"> true </param>
<param type="bool" name="Coulomb"> true </param>
Expand Down
1 change: 1 addition & 0 deletions src/Framework/Interaction/InitialState.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
14 changes: 14 additions & 0 deletions src/Physics/NuclearState/FermiMover.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -197,6 +197,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
Expand Down
1 change: 1 addition & 0 deletions src/Physics/NuclearState/FermiMover.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 {

Expand Down
1 change: 1 addition & 0 deletions src/Physics/NuclearState/NuclearModel.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,7 @@ typedef enum EFermiMoverInteractionType {
kFermiMoveEffectiveSF1p1h,
kFermiMoveEffectiveSF2p2h_eject,
kFermiMoveEffectiveSF2p2h_noeject,
kFermiMoveSpectralFunc,
} FermiMoverInteractionType_t;

class NuclearModel {
Expand Down
4 changes: 2 additions & 2 deletions src/Physics/NuclearState/SpectralFunc.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down
8 changes: 4 additions & 4 deletions src/Physics/QuasiElastic/EventGen/QELEventGenerator.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 = 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 < 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;
if (ilayer && (improvement_factor-1) < acceptable_fraction_of_safety_factor * (fSafetyFactor-1)) {
Expand Down