Add configurable Q2 minimum cuts and spline consistency metadata

config/CommonPhaseSpaceCuts.xml

@@ -0,0 +1,17 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+
+<!--
+***************************************************************************************************
+This file defines configurable kinematic phase-space cuts that are shared by the generator.
+Physical limits remain in KineUtils/KPhaseSpace; these values are user-configurable cuts layered
+on top of those physical limits.
+***************************************************************************************************
+-->
+
+<common_PhaseSpaceCuts_list>
+
+  <param_set name="Default">
+    <param type="double" name="EM-Q2-min"> 0.02 </param>
+  </param_set>
+
+</common_PhaseSpaceCuts_list>

config/G18_10a/G18_10a_02_11a/CommonPhaseSpaceCuts.xml

@@ -0,0 +1,17 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+
+<!--
+***************************************************************************************************
+This file defines configurable kinematic phase-space cuts that are shared by the generator.
+Physical limits remain in KineUtils/KPhaseSpace; these values are user-configurable cuts layered
+on top of those physical limits.
+***************************************************************************************************
+-->
+
+<common_PhaseSpaceCuts_list>
+
+  <param_set name="Default">
+    <param type="double" name="EM-Q2-min"> 0.02 </param>
+  </param_set>
+
+</common_PhaseSpaceCuts_list>

config/GEM21_11a/CommonPhaseSpaceCuts.xml

@@ -0,0 +1,17 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+
+<!--
+***************************************************************************************************
+This file defines configurable kinematic phase-space cuts that are shared by the generator.
+Physical limits remain in KineUtils/KPhaseSpace; these values are user-configurable cuts layered
+on top of those physical limits.
+***************************************************************************************************
+-->
+
+<common_PhaseSpaceCuts_list>
+
+  <param_set name="Default">
+    <param type="double" name="EM-Q2-min"> 0.02 </param>
+  </param_set>
+
+</common_PhaseSpaceCuts_list>

config/GEM21_11b/CommonPhaseSpaceCuts.xml

@@ -0,0 +1,17 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+
+<!--
+***************************************************************************************************
+This file defines configurable kinematic phase-space cuts that are shared by the generator.
+Physical limits remain in KineUtils/KPhaseSpace; these values are user-configurable cuts layered
+on top of those physical limits.
+***************************************************************************************************
+-->
+
+<common_PhaseSpaceCuts_list>
+
+  <param_set name="Default">
+    <param type="double" name="EM-Q2-min"> 0.02 </param>
+  </param_set>
+
+</common_PhaseSpaceCuts_list>

config/GEM21_11c/CommonPhaseSpaceCuts.xml

@@ -0,0 +1,17 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+
+<!--
+***************************************************************************************************
+This file defines configurable kinematic phase-space cuts that are shared by the generator.
+Physical limits remain in KineUtils/KPhaseSpace; these values are user-configurable cuts layered
+on top of those physical limits.
+***************************************************************************************************
+-->
+
+<common_PhaseSpaceCuts_list>
+
+  <param_set name="Default">
+    <param type="double" name="EM-Q2-min"> 0.02 </param>
+  </param_set>
+
+</common_PhaseSpaceCuts_list>

config/GEM21_11d/CommonPhaseSpaceCuts.xml

@@ -0,0 +1,17 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+
+<!--
+***************************************************************************************************
+This file defines configurable kinematic phase-space cuts that are shared by the generator.
+Physical limits remain in KineUtils/KPhaseSpace; these values are user-configurable cuts layered
+on top of those physical limits.
+***************************************************************************************************
+-->
+
+<common_PhaseSpaceCuts_list>
+
+  <param_set name="Default">
+    <param type="double" name="EM-Q2-min"> 0.02 </param>
+  </param_set>
+
+</common_PhaseSpaceCuts_list>

src/Apps/gEvGen.cxx

@@ -31,6 +31,7 @@
              [--force-flux-ray-interaction]
              [--seed random_number_seed]
              [--cross-sections xml_file]
+             [--q2-min q2_value]
 
              // command line args handled by RunOpt:
              [--event-generator-list list_name] // default "Default"
@@ -105,6 +106,9 @@
            --cross-sections
               Name (incl. full path) of an XML file with pre-computed
               cross-section values used for constructing splines.
+           --q2-min
+              Override the minimum Q^2 phase-space cut (in GeV^2). For weak
+              interactions, this explicitly opts the run into the Q2 cut.
 
            --event-generator-list
               List of event generators to load in event generation drivers.
@@ -146,6 +150,8 @@
 
 #include <cstdlib>
 #include <cassert>
+#include <algorithm>
+#include <cmath>
 #include <sstream>
 #include <string>
 #include <vector>
@@ -171,6 +177,7 @@
 #include "Framework/EventGen/GMCJDriver.h"
 #include "Framework/EventGen/GMCJMonitor.h"
 #include "Framework/Interaction/Interaction.h"
+#include "Framework/Interaction/KPhaseSpaceCuts.h"
 #include "Framework/Messenger/Messenger.h"
 #include "Framework/Ntuple/NtpWriter.h"
 #include "Framework/Ntuple/NtpMCFormat.h"
@@ -207,6 +214,7 @@ using namespace genie::controls;
 void GetCommandLineArgs (int argc, char ** argv);
 void Initialize         (void);
 void PrintSyntax        (void);
+void ValidateQ2MinAgainstSplineMetadata(void);
 
 #ifdef __CAN_GENERATE_EVENTS_USING_A_FLUX_OR_TGTMIX__
 void            GenerateEventsUsingFluxOrTgtMix();
@@ -240,6 +248,8 @@ long int        gOptRanSeed;      // random number seed
 string          gOptInpXSecFile;  // cross-section splines
 string          gOptOutFileName;  // Optional outfile name
 string          gOptStatFileName; // Status file name, set if gOptOutFileName was set.
+double          gOptQ2Min;        // Q2 minimum override value
+bool            gOptQ2MinSet;     // whether --q2-min was specified
 
 //____________________________________________________________________________
 int main(int argc, char ** argv)
@@ -279,17 +289,75 @@ void Initialize()
   }
   RunOpt::Instance()->BuildTune();
 
+  // Apply Q2 minimum override before physics code queries phase-space cuts.
+  if(gOptQ2MinSet) {
+    KPhaseSpaceCuts::Instance()->SetQ2MinOverride(gOptQ2Min);
+    LOG("gevgen", pNOTICE)
+      << "Overriding Q2 minimum phase-space cut from command line: "
+      << gOptQ2Min << " GeV^2";
+  }
+
   // Initialization of random number generators, cross-section table,
   // messenger thresholds, cache file
   utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
   utils::app_init::CacheFile(RunOpt::Instance()->CacheFile());
   utils::app_init::RandGen(gOptRanSeed);
   utils::app_init::XSecTable(gOptInpXSecFile, false);
+  ValidateQ2MinAgainstSplineMetadata();
 
   // Set GHEP print level
   GHepRecord::SetPrintLevel(RunOpt::Instance()->EventRecordPrintLevel());
 }
 //____________________________________________________________________________
+void ValidateQ2MinAgainstSplineMetadata(void)
+{
+  if(gOptInpXSecFile.empty()) return;
+
+  KPhaseSpaceCuts * cuts = KPhaseSpaceCuts::Instance();
+  if(!cuts->HasSplineQ2MinCutForProbe(gOptNuPdgCode)) return;
+
+  double configured_q2_min = cuts->SplineQ2MinCutForProbe(gOptNuPdgCode);
+  string configured_source =
+    cuts->SplineQ2MinCutSourceForProbe(gOptNuPdgCode);
+
+  XSecSplineList * xspl = XSecSplineList::Instance();
+  double spline_q2_min = -1.;
+  string unit = "";
+  string source = "";
+
+  if(!xspl->GetCurrentTuneQ2MinKinematics(spline_q2_min, &unit, &source)) {
+    LOG("gevgen", pFATAL)
+      << "Input cross-section spline file [" << gOptInpXSecFile
+      << "] does not declare <kinematics><q2 min=\"...\"/> metadata "
+      << "for tune " << RunOpt::Instance()->Tune()->Name()
+      << ". Refusing to run because the configured Q2 cut from "
+      << configured_source << " cannot be checked against the spline.";
+    gAbortingInErr = true;
+    exit(1);
+  }
+
+  const double abs_tol = 1E-12;
+  const double rel_tol = 1E-12;
+  const double diff = std::fabs(configured_q2_min - spline_q2_min);
+  const double scale = std::max(std::fabs(configured_q2_min), std::fabs(spline_q2_min));
+  if(diff > abs_tol && diff > rel_tol * scale) {
+    LOG("gevgen", pFATAL)
+      << "Q2 minimum cut mismatch between the configured run cut and "
+      << "input spline file [" << gOptInpXSecFile << "]. "
+      << configured_source << " sets " << configured_q2_min
+      << " GeV^2, but the spline metadata records " << spline_q2_min
+      << " " << unit << " from " << source << ".";
+    gAbortingInErr = true;
+    exit(1);
+  }
+
+  LOG("gevgen", pNOTICE)
+    << "Validated configured Q2 minimum " << configured_q2_min
+    << " GeV^2 from " << configured_source
+    << " against spline metadata for tune "
+    << RunOpt::Instance()->Tune()->Name();
+}
+//____________________________________________________________________________
 void GenerateEventsAtFixedInitState(void)
 {
   int neutrino = gOptNuPdgCode;
@@ -816,6 +884,14 @@ void GetCommandLineArgs(int argc, char ** argv)
     gOptInpXSecFile = "";
   }
 
+  // Q2 minimum override
+  gOptQ2MinSet = false;
+  if(parser.OptionExists("q2-min")) {
+    LOG("gevgen", pINFO) << "Reading Q2 minimum cut override";
+    gOptQ2Min = parser.ArgAsDouble("q2-min");
+    gOptQ2MinSet = true;
+  }
+
   //
   // print-out the command line options
   //
@@ -865,6 +941,10 @@ void GetCommandLineArgs(int argc, char ** argv)
       LOG("gevgen", pNOTICE)
           << " >> " <<  tgtpdgc << " (weight fraction = " << wgt << ")";
   }
+  if(gOptQ2MinSet) {
+     LOG("gevgen", pNOTICE)
+       << "Q2 minimum override: " << gOptQ2Min << " GeV^2";
+  }
   LOG("gevgen", pNOTICE) << "\n";
 
   LOG("gevgen", pNOTICE) << *RunOpt::Instance();
@@ -887,6 +967,7 @@ void PrintSyntax(void)
     << "\n              [--force-flux-ray-interaction]"
     << "\n              [--seed random_number_seed]"
     << "\n              [--cross-sections xml_file]"
+    << "\n              [--q2-min q2_value]"
     << RunOpt::RunOptSyntaxString(true)
     << "\n";
 

src/Apps/gMakeSplines.cxx

@@ -19,6 +19,7 @@
                   [--no-copy]
                   [--seed seed_number]
                   [--input-cross-sections xml_file]
+                  [--q2-min q2_value]
 
                   // command line args handled by RunOpt:
                   [--event-generator-list list_name] // default "Default"
@@ -59,6 +60,9 @@
               Name (incl. full path) of an XML file with pre-computed
               free-nucleon cross-section values. If loaded, it can speed-up
               cross-section calculation for nuclear targets.
+           --q2-min
+              Override the minimum Q^2 phase-space cut (in GeV^2). For weak
+              interactions, this explicitly opts the run into the Q2 cut.
 
            --event-generator-list
               List of event generators to load in event generation drivers.
@@ -87,6 +91,7 @@
 //____________________________________________________________________________
 
 #include <cassert>
+#include <cmath>
 #include <cstdlib>
 #include <string>
 #include <vector>
@@ -100,6 +105,7 @@
 #include "Framework/Conventions/GBuild.h"
 #include "Framework/EventGen/GEVGDriver.h"
 #include "Framework/Interaction/Interaction.h"
+#include "Framework/Interaction/KPhaseSpaceCuts.h"
 #include "Framework/Messenger/Messenger.h"
 #include "Framework/Numerical/RandomGen.h"
 #include "Framework/ParticleData/PDGCodeList.h"
@@ -140,6 +146,8 @@ bool     gOptNoCopy         = false;
 long int gOptRanSeed        = -1;   // random number seed
 string   gOptInpXSecFile    = "";   // input cross-section file
 string   gOptOutXSecFile    = "";   // output cross-section file
+double   gOptQ2Min;                 // Q2 minimum override value
+bool     gOptQ2MinSet       = false; // whether --q2-min was specified
 
 //____________________________________________________________________________
 int main(int argc, char ** argv)
@@ -153,6 +161,14 @@ int main(int argc, char ** argv)
   }
   RunOpt::Instance()->BuildTune();
 
+  // Apply Q2 minimum override before physics code queries phase-space cuts.
+  if(gOptQ2MinSet) {
+    KPhaseSpaceCuts::Instance()->SetQ2MinOverride(gOptQ2Min);
+    LOG("gmkspl", pNOTICE)
+      << "Overriding Q2 minimum phase-space cut from command line: "
+      << gOptQ2Min << " GeV^2";
+  }
+
   // throw on NaNs and Infs...
 #if defined(HAVE_FENV_H) && defined(HAVE_FEENABLEEXCEPT)
   feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
@@ -202,6 +218,31 @@ int main(int argc, char ** argv)
 
   // Save the splines at the requested XML file
   XSecSplineList * xspl = XSecSplineList::Instance();
+  KPhaseSpaceCuts * cuts = KPhaseSpaceCuts::Instance();
+  bool has_q2_metadata = false;
+  double q2_metadata = -1.;
+  string q2_metadata_source = "";
+  for(nuiter = neutrinos->begin(); nuiter != neutrinos->end(); ++nuiter) {
+    int nupdgc = *nuiter;
+    if(!cuts->HasSplineQ2MinCutForProbe(nupdgc)) continue;
+
+    double probe_q2_metadata = cuts->SplineQ2MinCutForProbe(nupdgc);
+    string probe_q2_metadata_source =
+      cuts->SplineQ2MinCutSourceForProbe(nupdgc);
+    if(!has_q2_metadata) {
+      has_q2_metadata = true;
+      q2_metadata = probe_q2_metadata;
+      q2_metadata_source = probe_q2_metadata_source;
+    } else if(std::fabs(q2_metadata - probe_q2_metadata) > 1E-12) {
+      LOG("gmkspl", pFATAL)
+        << "Cannot write one spline Q2 metadata value for mixed probes "
+        << "with different configured Q2 cuts.";
+      exit(1);
+    }
+  }
+  if(has_q2_metadata) {
+    xspl->SetCurrentTuneQ2MinKinematics(q2_metadata, q2_metadata_source);
+  }
   bool save_init = !gOptNoCopy;
   xspl->SaveAsXml(gOptOutXSecFile, save_init);
 
@@ -329,6 +370,14 @@ void GetCommandLineArgs(int argc, char ** argv)
     gOptInpXSecFile = "";
   }
 
+  // Q2 minimum override
+  gOptQ2MinSet = false;
+  if(parser.OptionExists("q2-min")) {
+    LOG("gmkspl", pINFO) << "Reading Q2 minimum cut override";
+    gOptQ2Min = parser.ArgAsDouble("q2-min");
+    gOptQ2MinSet = true;
+  }
+
   //
   // print the command-line options
   //
@@ -340,8 +389,12 @@ void GetCommandLineArgs(int argc, char ** argv)
      << "\n Input ROOT geometry : " << gOptGeomFilename
      << "\n Output cross-section file : " << gOptOutXSecFile
      << "\n Input cross-section file : " << gOptInpXSecFile
-     << "\n Random number seed : " << gOptRanSeed
-     << "\n";
+     << "\n Random number seed : " << gOptRanSeed;
+  if(gOptQ2MinSet) {
+     LOG("gmkspl", pNOTICE)
+       << "\n Q2 minimum override : " << gOptQ2Min << " GeV^2";
+  }
+  LOG("gmkspl", pNOTICE) << "\n";
 
   LOG("gmkspl", pNOTICE) << *RunOpt::Instance();
 }
@@ -358,6 +411,7 @@ void PrintSyntax(void)
     << "\n    [--no-copy]"
     << "\n    [--seed seed_number]"
     << "\n    [--input-cross-sections xml_file]"
+    << "\n    [--q2-min q2_value]"
     << RunOpt::RunOptSyntaxString(false)
     << "\n";