Skip to content

Conversation

@sjgardiner
Copy link
Member

Add in new systematic weights associated with hadron FSI channel and mean-free-path. Add in "model switch" uncertainties to the FSI fates from GENIE hA2018 to both INCL-like and G4 Bertini Cascade-like. Split up MFP systematic uncertainty by the nucleon kinetic energy.

This PR combines and supersedes Reweight #38 and GENIE-MC/Generator#450 originally submitted by @gputnam.

gputnam and others added 2 commits September 22, 2025 15:10
…mean-free-path. Add in model switch uncertainties to the FSI channel from GENIE hA to both INCL-like and G4 Bertini Cascade-like. Split up MFP systematic uncertainty by the nucleon kinetic energy.
…ata.

This clarifies that its sole purpose is to manage fate fraction splines
for reweighting from hA2018 to a G4-like or INCL-like model.
@sjgardiner sjgardiner force-pushed the feature/gputnam-fsi-model-variations branch from 7fe6054 to 40af3dd Compare October 8, 2025 23:25
kinetic-energy-dependent MFP dials so that they live in derived classes.
This leaves the original set of MFP and fate fraction dials untouched
apart from adding dummy arguments to interfaces.
@sjgardiner sjgardiner force-pushed the feature/gputnam-fsi-model-variations branch from 40af3dd to 69a4784 Compare October 8, 2025 23:34
@sjgardiner sjgardiner requested a review from nusense October 9, 2025 17:47
Copy link
Member

@nusense nusense left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There are a few comments / questions included here. But because it's factorized to not disturb the original code and is limited in its use, I'll approve without requiring changes.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The second column seems to be labelled "pA tot" and always 1, the 7th column is labelled "pA xsec". I infer that columns 3-6 are the individual fate fractions. The thing that bothers me is that the sum doesn't add up to 1 for all rows. Is that expected? What's the purpose of "pA tot" (or is it just historical)?

Spline *fINCLFracNA_PiPro;

// constants
static constexpr double fMinKinEnergy = 0.;
Copy link
Member

@nusense nusense Oct 10, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These energies ... are in MeV or GeV? Maybe note that.

std::string datafile_NAG4 = data_dir + "intranuke-fractions-NA-G4.dat";
std::string datafile_NAINCL = data_dir + "intranuke-fractions-NA-INCL.dat";

// verify files exist
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When compiled in optimized (ie. not debug) mode assert() can be a no-op and these wouldn't fail. I might be better to explicitly use the message service to print a warning if either is missing and then explicitly use exit(N). Not a show stopper but a nicety.

virtual void SetTwkDial(GSyst_t s, double val) override;

ModelSwitch_t fModelSwitch;
double fSystKELow; ///< Lower limit in kinetic energy range for this systematic
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

With at least some of the values kicked around in this code being in MeV vs GENIE's native GeV, it might be best to specify units for fSystKE* here as documentation


//.........................................................................
//
// nested class: MFP
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What does MFP stand for?

void GReWeightINukeParamsExtra::Fates::SetSystKERange(GSyst_t syst) {
if (syst == kINukeTwkDial_INCLLoE_N || syst == kINukeTwkDial_G4LoE_N) {
fSystKELow = 0;
fSystKEHigh = 0.15 / units::GeV;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are you using 0.15 / units::GeV to covert from "native units" to GeV? Usually I think when expressing a value one gives it a 0.15 * unit::GeV if the 0.15 is in GeV to convert it to "native units". I think this works here because GeV is 1.0. At least that's the way I remember it being.

}
else if (syst == kINukeTwkDial_INCLHiE_N || syst == kINukeTwkDial_G4HiE_N) {
fSystKELow = 0.6 / units::GeV;
fSystKEHigh = -1;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see that in IsInSystKERange() the -1 values is handled as "unspecified". One might note it here.

{
// Leave inelastic as cushion term to absorb the residual fraction in the systematic.
// Since the fractions sum to 1, this will set the value correctly.
if (syst != kINukeTwkDial_FrInel_N) {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does this mean that kINukeTwkDial_FrInel_N is always the cushion term? Or is this term somehow special? Perhaps it can be found a different way?


// set all the syst values to 1 -- these will be scaled by the model difference
int i=0;
while( (syst = (fHadType == kRwINukePion) ?
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This loop is a bit tricky-- I'm not sure I understand the intent, but I'll go along with it assuming it does what's desired.

void GReWeightINukeParamsExtra::MFP::SetSystKERange(GSyst_t syst) {
if (syst == kINukeTwkDial_MFPLoE_N) {
fSystKELow = 0;
fSystKEHigh = 0.15 / units::GeV;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same question as above about unit conversion scheme

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants