Constraints on W^2 for DIS

[pweigel]

I think there should be an additional constraint on the hadronic final state invariant mass ($W^2$) in the kinematics checks for DIS:

SIREN/projects/interactions/private/DISFromSpline.cxx

Lines 331 to 363 in 7d2db08

	double DISFromSpline::DifferentialCrossSection(double energy, double x, double y, double secondary_lepton_mass, double Q2) const {
	double log_energy = log10(energy);
	// check preconditions
	if(log_energy < differential_cross_section_.lower_extent(0)
	|| log_energy>differential_cross_section_.upper_extent(0))
	return 0.0;
	if(x <= 0 || x >= 1)
	return 0.0;
	if(y <= 0 || y >= 1)
	return 0.0;
	// we assume that:
	// the target is stationary so its energy is just its mass
	// the incoming neutrino is massless, so its kinetic energy is its total energy
	if(std::isnan(Q2)) {
	Q2 = 2.0 * energy * target_mass_ * x * y;
	}
	if(Q2 < minimum_Q2_) // cross section not calculated, assumed to be zero
	return 0;
	// cross section should be zero, but this check is missing from the original
	// CSMS calculation, so we must add it here
	if(!kinematicallyAllowed(x, y, energy, target_mass_, secondary_lepton_mass))
	return 0;
	std::array<double,3> coordinates{{log_energy, log10(x), log10(y)}};
	std::array<int,3> centers;
	if(!differential_cross_section_.searchcenters(coordinates.data(), centers.data()))
	return 0;
	double result = pow(10., differential_cross_section_.ndsplineeval(coordinates.data(), centers.data(), 0));
	assert(result >= 0);
	return unit * result;
	}

The typical constraint for DIS is:
$$W^2 = m_N^2 + Q^2 (\frac{1}{x} - 1) > 4~\rm{GeV}^2$$

Other codes include some or all of the shallow inelastic scattering region:
$$(m_N + m_\pi)^2 < W^2 < 4~\rm{GeV}^2$$

See https://arxiv.org/abs/2307.09241 for details about the effects of including/excluding the low $$W^2$$ region.