New GENIE Hadron FSI Channel + MFP Systematic Weights

config/GSystUncertaintyTable.xml

@@ -217,6 +217,31 @@ values for each Reweight tweak dial. All of these have default values of zero.
    <param type="double" name="FrPiProd_N@PlusOneSigma"> 0.20 </param>
    <param type="double" name="FrPiProd_N@MinusOneSigma"> 0.20 </param>
 
+   <param type="double" name="FrG4_N@PlusOneSigma"> 1.0 </param>
+   <param type="double" name="FrG4_N@MinusOneSigma"> 1.0 </param>
+   <param type="double" name="FrINCL_N@PlusOneSigma"> 1.0 </param>
+   <param type="double" name="FrINCL_N@MinusOneSigma"> 1.0 </param>
+
+   <param type="double" name="FrG4LoE_N@PlusOneSigma"> 1.0 </param>
+   <param type="double" name="FrG4LoE_N@MinusOneSigma"> 1.0 </param>
+   <param type="double" name="FrINCLLoE_N@PlusOneSigma"> 1.0 </param>
+   <param type="double" name="FrINCLLoE_N@MinusOneSigma"> 1.0 </param>
+
+   <param type="double" name="FrG4M1E_N@PlusOneSigma"> 1.0 </param>
+   <param type="double" name="FrG4M1E_N@MinusOneSigma"> 1.0 </param>
+   <param type="double" name="FrINCLM1E_N@PlusOneSigma"> 1.0 </param>
+   <param type="double" name="FrINCLM1E_N@MinusOneSigma"> 1.0 </param>
+
+   <param type="double" name="FrG4M2E_N@PlusOneSigma"> 1.0 </param>
+   <param type="double" name="FrG4M2E_N@MinusOneSigma"> 1.0 </param>
+   <param type="double" name="FrINCLM2E_N@PlusOneSigma"> 1.0 </param>
+   <param type="double" name="FrINCLM2E_N@MinusOneSigma"> 1.0 </param>
+
+   <param type="double" name="FrG4HiE_N@PlusOneSigma"> 1.0 </param>
+   <param type="double" name="FrG4HiE_N@MinusOneSigma"> 1.0 </param>
+   <param type="double" name="FrINCLHiE_N@PlusOneSigma"> 1.0 </param>
+   <param type="double" name="FrINCLHiE_N@MinusOneSigma"> 1.0 </param>
+
    <param type="double" name="RDecBR1gamma@PlusOneSigma"> 0.50 </param>
    <param type="double" name="RDecBR1gamma@MinusOneSigma"> 0.50 </param>
 

data/intranuke/tot_xsec/README.md

@@ -0,0 +1,11 @@
+This folder contains evaluations of energy-dependent nucleon fate fractions in a
+format similar to data/evgen/intranuke/tot_xsec/intranuke-fractions-NA2016.dat
+from the Generator repository. However, the fate fractions given here are
+evaluated according to alternative hadron transport models with an approximate
+mapping between the relevant reaction channels and hA2018 fates:
+
+intranuke-fractions-NA-G4.dat = Geant4 Bertini cascade
+intranuke-fractions-NA-INCL.dat = Liège intranuclear cascade (INCL++)
+
+These alternative model tables were calculated and contributed by Gray Putnam
+as part of a 2025 systematic uncertainty quantification effort for SBN.

data/intranuke/tot_xsec/intranuke-fractions-NA-G4.dat

@@ -0,0 +1,28 @@
+# Provided & maintained by Gray Putnam
+#
+# Cross section fractions computed by G4 for p-Ar interactions.
+# KE is lab energy in MeV, cross sections in lab, units= mb
+# KE     pA tot     pA elas     pA inel      pA cex      pA abs      pA pipro   pA xsec
+#
+#  fate=               2            3           1          4            5
+0	1.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
+50	1.0000	0.0000	0.3510	0.1139	0.4800	0.0000	748.4410
+100	1.0000	0.0000	0.3424	0.0935	0.5015	0.0000	584.0650
+150	1.0000	0.0000	0.3528	0.0969	0.4634	0.0000	519.8320
+200	1.0000	0.0000	0.3501	0.1100	0.4518	0.0001	493.5160
+250	1.0000	0.0000	0.3796	0.1201	0.4433	0.0013	489.0440
+300	1.0000	0.0000	0.3845	0.1172	0.4305	0.0070	483.5990
+350	1.0000	0.0000	0.3685	0.1210	0.4196	0.0195	480.5530
+400	1.0000	0.0000	0.3540	0.1138	0.4146	0.0374	484.1820
+450	1.0000	0.0000	0.3372	0.1013	0.4117	0.0617	495.4610
+500	1.0000	0.0000	0.3163	0.0923	0.4210	0.0935	506.7390
+550	1.0000	0.0000	0.2984	0.0840	0.4391	0.1239	533.2490
+600	1.0000	0.0000	0.2674	0.0730	0.4175	0.1634	541.2860
+650	1.0000	0.0000	0.2622	0.0624	0.4191	0.2020	560.8610
+700	1.0000	0.0000	0.2474	0.0573	0.4204	0.2201	575.5740
+750	1.0000	0.0000	0.2325	0.0479	0.4152	0.2480	577.9080
+800	1.0000	0.0000	0.2225	0.0458	0.4017	0.2784	579.7230
+850	1.0000	0.0000	0.2043	0.0356	0.3837	0.3045	579.2690
+900	1.0000	0.0000	0.2006	0.0339	0.3678	0.3293	579.6580
+950	1.0000	0.0000	0.1890	0.0333	0.3656	0.3429	580.3710
+1000	1.0000	0.0000	0.1897	0.0316	0.3406	0.3575	579.4630

data/intranuke/tot_xsec/intranuke-fractions-NA-INCL.dat

@@ -0,0 +1,28 @@
+# Provided & maintained by Gray Putnam
+#
+# Cross section fractions computed by INCL for p-Ar interactions.
+# KE is lab energy in MeV, cross sections in lab, units= mb
+# KE     pA tot     pA elas     pA inel      pA cex      pA abs      pA pipro   pA xsec
+#
+#  fate=               2            3           1          4            5
+0	1.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
+50	1.0000	0.0000	0.3441	0.2931	0.3628	0.0000	792.0560
+100	1.0000	0.0000	0.2838	0.2302	0.4859	0.0000	623.0340
+150	1.0000	0.0000	0.2706	0.2097	0.5196	0.0000	569.3210
+200	1.0000	0.0000	0.2736	0.2000	0.5262	0.0002	541.1050
+250	1.0000	0.0000	0.2863	0.1912	0.5191	0.0034	517.9060
+300	1.0000	0.0000	0.2957	0.1785	0.5144	0.0115	514.8660
+350	1.0000	0.0000	0.2987	0.1650	0.5045	0.0318	517.4510
+400	1.0000	0.0000	0.2926	0.1478	0.5005	0.0591	526.3530
+450	1.0000	0.0000	0.2817	0.1297	0.4947	0.0939	543.3740
+500	1.0000	0.0000	0.2660	0.1162	0.4877	0.1302	548.9100
+550	1.0000	0.0000	0.2533	0.0997	0.4759	0.1712	564.0840
+600	1.0000	0.0000	0.2483	0.0860	0.4600	0.2057	587.4950
+650	1.0000	0.0000	0.2354	0.0736	0.4527	0.2382	593.1100
+700	1.0000	0.0000	0.2248	0.0704	0.4329	0.2720	608.9550
+750	1.0000	0.0000	0.2219	0.0586	0.4190	0.3005	612.4780
+800	1.0000	0.0000	0.2154	0.0554	0.4049	0.3244	611.2830
+850	1.0000	0.0000	0.2091	0.0526	0.3854	0.3529	624.1370
+900	1.0000	0.0000	0.2100	0.0485	0.3675	0.3740	622.2430
+950	1.0000	0.0000	0.2052	0.0444	0.3596	0.3908	623.5250
+1000	1.0000	0.0000	0.2003	0.0423	0.3377	0.4197	630.3000

src/Apps/gRwght1Param.cxx

@@ -118,6 +118,7 @@
 #include "RwCalculators/GReWeightResonanceDecay.h"
 #include "RwCalculators/GReWeightFZone.h"
 #include "RwCalculators/GReWeightINuke.h"
+#include "RwCalculators/GReWeightINukeExtra.h"
 #include "RwCalculators/GReWeightAGKY.h"
 #include "RwCalculators/GReWeightNuXSecCCQEaxial.h"
 #include "RwCalculators/GReWeightNuXSecCCQEvec.h"
@@ -247,6 +248,7 @@ int main(int argc, char ** argv)
   rw.AdoptWghtCalc( "hadro_res_decay", new GReWeightResonanceDecay  );
   rw.AdoptWghtCalc( "hadro_fzone",     new GReWeightFZone           );
   rw.AdoptWghtCalc( "hadro_intranuke", new GReWeightINuke           );
+  rw.AdoptWghtCalc( "hadro_intranuke_extra", new GReWeightINukeExtra);
   rw.AdoptWghtCalc( "hadro_agky",      new GReWeightAGKY            );
 
   // GReWeightDISNuclMod::CalcWeight() not implemented - don't try to use it ..

src/RwCalculators/GReWeightINuke.cxx

@@ -81,6 +81,8 @@ GReWeightModel("IntraNuke")
   }
 
   fFSIModel->AdoptSubstructure();
+
+  fINukeRwParams = std::make_shared< GReWeightINukeParams >();
 }
 //_______________________________________________________________________________________
 GReWeightINuke::~GReWeightINuke()
@@ -141,19 +143,19 @@ bool GReWeightINuke::AppliesTo(const EventRecord & event) const
 void GReWeightINuke::SetSystematic(GSyst_t syst, double val)
 {
   if(this->IsHandled(syst)) {
-     fINukeRwParams.SetTwkDial(syst, val);
+     fINukeRwParams->SetTwkDial(syst, val);
   }
 }
 //_______________________________________________________________________________________
 void GReWeightINuke::Reset(void)
 {
-  fINukeRwParams.Reset();
+  fINukeRwParams->Reset();
   this->Reconfigure();
 }
 //_______________________________________________________________________________________
 void GReWeightINuke::Reconfigure(void)
 {
-  fINukeRwParams.Reconfigure();
+  fINukeRwParams->Reconfigure();
 }
 //_______________________________________________________________________________________
 double GReWeightINuke::CalcWeight(const EventRecord & event)
@@ -166,7 +168,7 @@ double GReWeightINuke::CalcWeight(const EventRecord & event)
   if (A<=1) return 1.0;
   if (Z<=1) return 1.0;
 
-  fINukeRwParams.SetTargetA( A );
+  fINukeRwParams->SetTargetA( A );
 
   // Get the pre-FSI nuclear remnant. The A and Z values for this particle
   // (distinct from both the post-FSI remnant and the target nucleus)
@@ -274,14 +276,14 @@ double GReWeightINuke::CalcWeight(const EventRecord & event)
      double w_fate = 1.0;
 
      // Check which weights need to be calculated (only if relevant params were tweaked)
-     bool calc_w_mfp  = fINukeRwParams.MeanFreePathParams(pdgc)->IsTweaked();
-     bool calc_w_fate = fINukeRwParams.FateParams(pdgc)->IsTweaked();
+     bool calc_w_mfp  = fINukeRwParams->MeanFreePathParams(pdgc)->IsTweaked();
+     bool calc_w_fate = fINukeRwParams->FateParams(pdgc)->IsTweaked();
 
      // Compute weight to account for changes in the total rescattering probability
      double mfp_scale_factor = 1.;
      if(calc_w_mfp)
      {
-        mfp_scale_factor = fINukeRwParams.MeanFreePathParams(pdgc)->ScaleFactor();
+        mfp_scale_factor = fINukeRwParams->MeanFreePathParams(pdgc)->ScaleFactor(p4);
         w_mfp = utils::rew::MeanFreePathWeight( pdgc, x4, p4, A, Z,
           mfp_scale_factor, interacted, *fFSIModel );
      } // calculate mfp weight?
@@ -290,7 +292,7 @@ double GReWeightINuke::CalcWeight(const EventRecord & event)
      if(calc_w_fate && interacted)
      {
         double fate_fraction_scale_factor =
-             fINukeRwParams.FateParams(pdgc)->ScaleFactor(
+             fINukeRwParams->FateParams(pdgc)->ScaleFactor(
                   GSyst::INukeFate2GSyst((INukeFateHA_t)fsi_code,pdgc), p4);
         w_fate = fate_fraction_scale_factor;
      }

src/RwCalculators/GReWeightINuke.h

@@ -34,6 +34,9 @@
 
 //#define _G_REWEIGHT_INUKE_DEBUG_NTP_
 
+// Standard library includes
+#include <memory>
+
 // GENIE/Reweight includes
 #include "RwCalculators/GReWeightModel.h"
 #include "RwCalculators/GReWeightINukeParams.h"
@@ -56,7 +59,7 @@ namespace rew   {
  {
  public:
    GReWeightINuke();
-  ~GReWeightINuke();
+   virtual ~GReWeightINuke();
 
    // implement the GReWeightI interface
    bool   AppliesTo      (const EventRecord & event) const;
@@ -66,14 +69,14 @@ namespace rew   {
    void   Reconfigure    (void);
    double CalcWeight     (const EventRecord & event);
 
- private:
+ protected:
 
    void CalcDeltaAZ( const EventRecord& event, const GHepParticle& p,
      int& deltaA, int& deltaZ );
 
    void UpdateRemnantAZ( int deltaA, int deltaZ );
 
-   GReWeightINukeParams fINukeRwParams;
+   std::shared_ptr< GReWeightINukeParams > fINukeRwParams;
 
    HAIntranuke2018* fFSIModel;
 

src/RwCalculators/GReWeightINukeExtra.cxx

@@ -0,0 +1,67 @@
+//____________________________________________________________________________
+/*
+ Copyright (c) 2003-2025, The GENIE Collaboration
+ For the full text of the license visit http://copyright.genie-mc.org
+*/
+//____________________________________________________________________________
+
+// GENIE/Generator includes
+#include "Framework/EventGen/EventRecord.h"
+#include "Framework/Messenger/Messenger.h"
+
+// GENIE/Reweight includes
+#include "RwCalculators/GReWeightINukeExtra.h"
+#include "RwCalculators/GReWeightINukeParamsExtra.h"
+
+using namespace genie;
+using namespace genie::rew;
+
+//_______________________________________________________________________________________
+GReWeightINukeExtra::GReWeightINukeExtra() : GReWeightINuke()
+{
+  // Switch to a different name to avoid collisions with an instance
+  // of the base class
+  fName = "IntraNukeExtra";
+
+  // Switch to using the derived "params" object
+  fINukeRwParams = std::make_shared< GReWeightINukeParamsExtra >();
+}
+//_______________________________________________________________________________________
+GReWeightINukeExtra::~GReWeightINukeExtra()
+{
+}
+//_______________________________________________________________________________________
+bool GReWeightINukeExtra::IsHandled(GSyst_t syst) const {
+  bool handle;
+
+  switch (syst) {
+  // Nucleon MFP variations in kinetic energy bins
+  case (kINukeTwkDial_MFPLoE_N):
+  case (kINukeTwkDial_MFPM1E_N):
+  case (kINukeTwkDial_MFPM2E_N):
+  case (kINukeTwkDial_MFPHiE_N):
+
+  // Map hA2018 nucleon fates to either G4 or INCL
+  case (kINukeTwkDial_G4_N):
+  case (kINukeTwkDial_INCL_N):
+
+  // Map hA2018 nucleon fates to G4 in kinetic energy bins
+  case (kINukeTwkDial_G4LoE_N):
+  case (kINukeTwkDial_G4M1E_N):
+  case (kINukeTwkDial_G4M2E_N):
+  case (kINukeTwkDial_G4HiE_N):
+
+  // Map hA2018 nucleon fates to INCL in kinetic energy bins
+  case (kINukeTwkDial_INCLLoE_N):
+  case (kINukeTwkDial_INCLM1E_N):
+  case (kINukeTwkDial_INCLM2E_N):
+  case (kINukeTwkDial_INCLHiE_N):
+    handle = true;
+    break;
+
+  default:
+    handle = false;
+  }
+
+  return handle;
+}

src/RwCalculators/GReWeightINukeExtra.h

@@ -0,0 +1,45 @@
+//____________________________________________________________________________
+/*!
+
+\class    genie::rew::GReWeightINukeExtra
+
+\brief    Adds some additional functionality to reweighting hA2018.
+          Originally developed for use by SBND and ICARUS.
+
+\author   Gray Putnam <gputnam \at fnal.gov>
+          Steven Gardiner <gardiner \at fnal.gov>
+          Fermi National Accelerator Laboratory
+
+\created  Oct 1, 2025
+
+\cpright  Copyright (c) 2003-2025, The GENIE Collaboration
+          For the full text of the license visit http://copyright.genie-mc.org
+*/
+//____________________________________________________________________________
+
+#ifndef _G_REWEIGHT_INUKE_EXTRA_H_
+#define _G_REWEIGHT_INUKE_EXTRA_H_
+
+// GENIE/Reweight includes
+#include "RwCalculators/GReWeightINuke.h"
+#include "RwCalculators/GReWeightINukeParamsExtra.h"
+
+namespace genie {
+
+namespace rew   {
+
+ class GReWeightINukeExtra : public GReWeightINuke
+ {
+ public:
+   GReWeightINukeExtra();
+  ~GReWeightINukeExtra();
+
+   // Adjust some details of the base implementation to accommodate
+   // new dials while keeping the originals untouched
+   virtual bool IsHandled( GSyst_t syst ) const override final;
+ };
+
+} // rew
+} // genie
+
+#endif