diff --git a/config/CommonPhaseSpaceCuts.xml b/config/CommonPhaseSpaceCuts.xml new file mode 100644 index 0000000000..1f2d43780d --- /dev/null +++ b/config/CommonPhaseSpaceCuts.xml @@ -0,0 +1,17 @@ + + + + + + + + 0.02 + + + diff --git a/config/G18_10a/G18_10a_02_11a/CommonPhaseSpaceCuts.xml b/config/G18_10a/G18_10a_02_11a/CommonPhaseSpaceCuts.xml new file mode 100644 index 0000000000..1f2d43780d --- /dev/null +++ b/config/G18_10a/G18_10a_02_11a/CommonPhaseSpaceCuts.xml @@ -0,0 +1,17 @@ + + + + + + + + 0.02 + + + diff --git a/config/GEM21_11a/CommonPhaseSpaceCuts.xml b/config/GEM21_11a/CommonPhaseSpaceCuts.xml new file mode 100644 index 0000000000..1f2d43780d --- /dev/null +++ b/config/GEM21_11a/CommonPhaseSpaceCuts.xml @@ -0,0 +1,17 @@ + + + + + + + + 0.02 + + + diff --git a/config/GEM21_11b/CommonPhaseSpaceCuts.xml b/config/GEM21_11b/CommonPhaseSpaceCuts.xml new file mode 100644 index 0000000000..1f2d43780d --- /dev/null +++ b/config/GEM21_11b/CommonPhaseSpaceCuts.xml @@ -0,0 +1,17 @@ + + + + + + + + 0.02 + + + diff --git a/config/GEM21_11c/CommonPhaseSpaceCuts.xml b/config/GEM21_11c/CommonPhaseSpaceCuts.xml new file mode 100644 index 0000000000..1f2d43780d --- /dev/null +++ b/config/GEM21_11c/CommonPhaseSpaceCuts.xml @@ -0,0 +1,17 @@ + + + + + + + + 0.02 + + + diff --git a/config/GEM21_11d/CommonPhaseSpaceCuts.xml b/config/GEM21_11d/CommonPhaseSpaceCuts.xml new file mode 100644 index 0000000000..1f2d43780d --- /dev/null +++ b/config/GEM21_11d/CommonPhaseSpaceCuts.xml @@ -0,0 +1,17 @@ + + + + + + + + 0.02 + + + diff --git a/src/Apps/gEvGen.cxx b/src/Apps/gEvGen.cxx index 035c31bbd5..371bb1897b 100644 --- a/src/Apps/gEvGen.cxx +++ b/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 #include +#include +#include #include #include #include @@ -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 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"; diff --git a/src/Apps/gMakeSplines.cxx b/src/Apps/gMakeSplines.cxx index 8685b64324..0e5d390356 100644 --- a/src/Apps/gMakeSplines.cxx +++ b/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 +#include #include #include #include @@ -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"; diff --git a/src/Framework/Interaction/KPhaseSpace.cxx b/src/Framework/Interaction/KPhaseSpace.cxx index 989297803e..20e4bf01b0 100644 --- a/src/Framework/Interaction/KPhaseSpace.cxx +++ b/src/Framework/Interaction/KPhaseSpace.cxx @@ -23,6 +23,7 @@ #include "Framework/Conventions/Controls.h" #include "Framework/Interaction/Interaction.h" #include "Framework/Interaction/InteractionException.h" +#include "Framework/Interaction/KPhaseSpaceCuts.h" #include "Framework/Messenger/Messenger.h" #include "Framework/ParticleData/PDGLibrary.h" #include "Framework/ParticleData/PDGUtils.h" @@ -38,6 +39,22 @@ using namespace genie::constants; ClassImp(KPhaseSpace) +namespace { + void ApplyQ2MinCut(const Interaction * interaction, Range1D_t & q2lim, double default_q2min) + { + KPhaseSpaceCuts * cuts = KPhaseSpaceCuts::Instance(); + bool has_cut = cuts->HasQ2MinCut(interaction); + if(q2lim.min < 0. || q2lim.max < 0. || !has_cut) return; + + double q2min = cuts->Q2MinCut(interaction, default_q2min); + if(q2lim.min < q2min) q2lim.min = q2min; + if(q2lim.max < q2lim.min) { + q2lim.min = -1.; + q2lim.max = -1.; + } + } +} + //____________________________________________________________________________ KPhaseSpace::KPhaseSpace(void) : TObject(), fInteraction(NULL) @@ -74,6 +91,11 @@ double KPhaseSpace::GetTMaxDFR() } //___________________________________________________________________________ +double KPhaseSpace::GetQ2MinEM() +{ + return KPhaseSpaceCuts::Instance()->EMQ2MinCut(); +} +//___________________________________________________________________________ void KPhaseSpace::UseInteraction(const Interaction * in) { fInteraction = in; @@ -552,8 +574,10 @@ Range1D_t KPhaseSpace::Q2Lim_W(void) const } else if (is_dme || is_dmdis) { Q2l = kinematics::DarkQ2Lim_W(Ev,M,ml,W); } else { - Q2l = is_em ? kinematics::electromagnetic::InelQ2Lim_W(Ev,ml,M,W) : kinematics::InelQ2Lim_W(Ev,M,ml,W); + Q2l = is_em ? kinematics::electromagnetic::InelQ2Lim_W(Ev,ml,M,W,0.) : kinematics::InelQ2Lim_W(Ev,M,ml,W); } + ApplyQ2MinCut(fInteraction, Q2l, + pi.IsInverseBetaDecay() ? controls::kMinQ2Limit_VLE : (is_em ? 0. : controls::kMinQ2Limit)); return Q2l; } @@ -599,6 +623,7 @@ Range1D_t KPhaseSpace::Q2Lim(void) const if(is_cevns) { double Ev_lab = init_state.ProbeE(kRfLab); Q2l = kinematics::CEvNSQ2Lim(Ev_lab); + ApplyQ2MinCut(fInteraction, Q2l, is_em ? 0. : controls::kMinQ2Limit); return Q2l; } @@ -616,6 +641,7 @@ Range1D_t KPhaseSpace::Q2Lim(void) const } Q2l = kinematics::CohQ2Lim(M, m_other, ml, Ev); + ApplyQ2MinCut(fInteraction, Q2l, is_em ? 0. : controls::kMinQ2Limit); return Q2l; } @@ -632,8 +658,10 @@ Range1D_t KPhaseSpace::Q2Lim(void) const if (pi.IsInverseBetaDecay()) { Q2l = kinematics::InelQ2Lim_W(Ev,M,ml,W,controls::kMinQ2Limit_VLE); } else { - Q2l = is_em ? kinematics::electromagnetic::InelQ2Lim_W(Ev,ml,M,W) : kinematics::InelQ2Lim_W(Ev,M,ml,W); + Q2l = is_em ? kinematics::electromagnetic::InelQ2Lim_W(Ev,ml,M,W,0.) : kinematics::InelQ2Lim_W(Ev,M,ml,W); } + ApplyQ2MinCut(fInteraction, Q2l, + pi.IsInverseBetaDecay() ? controls::kMinQ2Limit_VLE : (is_em ? 0. : controls::kMinQ2Limit)); return Q2l; } @@ -662,9 +690,10 @@ Range1D_t KPhaseSpace::Q2Lim(void) const // TODO: Q2maxConfig if (pi.IsMEC()){ double W = fInteraction->RecoilNucleon()->Mass(); - Q2l = is_em ? kinematics::electromagnetic::InelQ2Lim_W(Ev,ml,M,W) : kinematics::InelQ2Lim_W(Ev,M,ml,W); + Q2l = is_em ? kinematics::electromagnetic::InelQ2Lim_W(Ev,ml,M,W,0.) : kinematics::InelQ2Lim_W(Ev,M,ml,W); double Q2maxConfig = 1.44; // need to pull from config file somehow? if (Q2l.max > Q2maxConfig) Q2l.max = Q2maxConfig; + ApplyQ2MinCut(fInteraction, Q2l, is_em ? 0. : controls::kMinQ2Limit); return Q2l; } @@ -674,7 +703,8 @@ Range1D_t KPhaseSpace::Q2Lim(void) const } // inelastic - Q2l = is_em ? kinematics::electromagnetic::InelQ2Lim(Ev,ml,M) : kinematics::InelQ2Lim(Ev,M,ml); + Q2l = is_em ? kinematics::electromagnetic::InelQ2Lim(Ev,ml,M,0.) : kinematics::InelQ2Lim(Ev,M,ml); + ApplyQ2MinCut(fInteraction, Q2l, is_em ? 0. : controls::kMinQ2Limit); return Q2l; } //____________________________________________________________________________ @@ -1147,4 +1177,3 @@ Range1D_t KPhaseSpace::Q2Lim_W_SPP_iso(void) const return Q2l; } //____________________________________________________________________________ - diff --git a/src/Framework/Interaction/KPhaseSpace.h b/src/Framework/Interaction/KPhaseSpace.h index d8f632a8fe..c6dd24b945 100644 --- a/src/Framework/Interaction/KPhaseSpace.h +++ b/src/Framework/Interaction/KPhaseSpace.h @@ -71,6 +71,7 @@ class KPhaseSpace : public TObject { Range1D_t Q2Lim_W_SPP_iso (void) const; ///< Q2 limits @ fixed W for resonance single pion production on isoscalar nucleon static double GetTMaxDFR(); + static double GetQ2MinEM(); ///< Deprecated; use KPhaseSpaceCuts instead. private: void Init(void); diff --git a/src/Framework/Interaction/KPhaseSpaceCuts.cxx b/src/Framework/Interaction/KPhaseSpaceCuts.cxx new file mode 100644 index 0000000000..349aec59ef --- /dev/null +++ b/src/Framework/Interaction/KPhaseSpaceCuts.cxx @@ -0,0 +1,229 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2025, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org +*/ +//____________________________________________________________________________ + +#include +#include + +#include + +#include "Framework/Algorithm/AlgConfigPool.h" +#include "Framework/Interaction/Interaction.h" +#include "Framework/Interaction/KPhaseSpaceCuts.h" +#include "Framework/Interaction/ProcessInfo.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/ParticleData/PDGUtils.h" +#include "Framework/Registry/Registry.h" + +using namespace genie; + +//____________________________________________________________________________ +KPhaseSpaceCuts * KPhaseSpaceCuts::fInstance = 0; +//____________________________________________________________________________ +KPhaseSpaceCuts::KPhaseSpaceCuts() : +fLoaded(false), +fHasEMQ2Min(false), +fHasWeakQ2Min(false), +fEMQ2Min(-1.), +fWeakQ2Min(-1.), +fHasQ2MinOverride(false), +fQ2MinOverride(-1.) +{ + fInstance = 0; +} +//____________________________________________________________________________ +KPhaseSpaceCuts::~KPhaseSpaceCuts() +{ + fInstance = 0; +} +//____________________________________________________________________________ +KPhaseSpaceCuts * KPhaseSpaceCuts::Instance(void) +{ + if(fInstance == 0) { + static KPhaseSpaceCuts::Cleaner cleaner; + cleaner.DummyMethodAndSilentCompiler(); + fInstance = new KPhaseSpaceCuts; + } + return fInstance; +} +//____________________________________________________________________________ +void KPhaseSpaceCuts::SetQ2MinOverride(double q2min) +{ + fQ2MinOverride = this->ValidateQ2Min("command-line Q2-min", q2min); + fHasQ2MinOverride = true; +} +//____________________________________________________________________________ +bool KPhaseSpaceCuts::HasQ2MinCut(const Interaction * interaction) const +{ + this->LoadConfig(); + + if(!interaction) return fHasQ2MinOverride || fHasEMQ2Min || fHasWeakQ2Min; + + const ProcessInfo & proc = interaction->ProcInfo(); + if(proc.IsEM()) { + if(fHasQ2MinOverride) return true; + if(!fHasEMQ2Min) this->FailMissingEMQ2Min(); + return true; + } + + if(proc.IsWeak()) { + return fHasQ2MinOverride || fHasWeakQ2Min; + } + + return false; +} +//____________________________________________________________________________ +double KPhaseSpaceCuts::Q2MinCut( + const Interaction * interaction, double default_q2min) const +{ + this->LoadConfig(); + + if(!interaction) { + if(fHasQ2MinOverride) return this->RaiseDefault(default_q2min, fQ2MinOverride); + return default_q2min; + } + + const ProcessInfo & proc = interaction->ProcInfo(); + if(proc.IsEM()) { + if(fHasQ2MinOverride) return this->RaiseDefault(default_q2min, fQ2MinOverride); + if(!fHasEMQ2Min) this->FailMissingEMQ2Min(); + return this->RaiseDefault(default_q2min, fEMQ2Min); + } + + if(proc.IsWeak()) { + if(fHasQ2MinOverride) return this->RaiseDefault(default_q2min, fQ2MinOverride); + if(fHasWeakQ2Min) return this->RaiseDefault(default_q2min, fWeakQ2Min); + } + + return default_q2min; +} +//____________________________________________________________________________ +double KPhaseSpaceCuts::EMQ2MinCut(void) const +{ + this->LoadConfig(); + + if(fHasQ2MinOverride) return fQ2MinOverride; + if(!fHasEMQ2Min) this->FailMissingEMQ2Min(); + return fEMQ2Min; +} +//____________________________________________________________________________ +bool KPhaseSpaceCuts::HasSplineQ2MinCut(void) const +{ + this->LoadConfig(); + + return fHasQ2MinOverride || fHasWeakQ2Min || fHasEMQ2Min; +} +//____________________________________________________________________________ +double KPhaseSpaceCuts::SplineQ2MinCut(void) const +{ + this->LoadConfig(); + + if(fHasQ2MinOverride) return fQ2MinOverride; + if(fHasWeakQ2Min) return fWeakQ2Min; + if(fHasEMQ2Min) return fEMQ2Min; + return -1.; +} +//____________________________________________________________________________ +string KPhaseSpaceCuts::SplineQ2MinCutSource(void) const +{ + this->LoadConfig(); + + if(fHasQ2MinOverride) return "command-line --q2-min"; + if(fHasWeakQ2Min || fHasEMQ2Min) return "CommonPhaseSpaceCuts.xml"; + return ""; +} +//____________________________________________________________________________ +bool KPhaseSpaceCuts::HasSplineQ2MinCutForProbe(int probe_pdg) const +{ + this->LoadConfig(); + + if(fHasQ2MinOverride) return true; + if(pdg::IsNeutrino(probe_pdg) || pdg::IsAntiNeutrino(probe_pdg)) { + return fHasWeakQ2Min; + } + if(pdg::IsChargedLepton(probe_pdg)) { + return fHasEMQ2Min; + } + return fHasWeakQ2Min || fHasEMQ2Min; +} +//____________________________________________________________________________ +double KPhaseSpaceCuts::SplineQ2MinCutForProbe(int probe_pdg) const +{ + this->LoadConfig(); + + if(fHasQ2MinOverride) return fQ2MinOverride; + if(pdg::IsNeutrino(probe_pdg) || pdg::IsAntiNeutrino(probe_pdg)) { + if(fHasWeakQ2Min) return fWeakQ2Min; + return -1.; + } + if(pdg::IsChargedLepton(probe_pdg)) { + if(fHasEMQ2Min) return fEMQ2Min; + return -1.; + } + if(fHasWeakQ2Min) return fWeakQ2Min; + if(fHasEMQ2Min) return fEMQ2Min; + return -1.; +} +//____________________________________________________________________________ +string KPhaseSpaceCuts::SplineQ2MinCutSourceForProbe(int probe_pdg) const +{ + this->LoadConfig(); + + if(!this->HasSplineQ2MinCutForProbe(probe_pdg)) return ""; + if(fHasQ2MinOverride) return "command-line --q2-min"; + return "CommonPhaseSpaceCuts.xml"; +} +//____________________________________________________________________________ +void KPhaseSpaceCuts::LoadConfig(void) const +{ + if(fLoaded) return; + + fLoaded = true; + + AlgConfigPool * confp = AlgConfigPool::Instance(); + const Registry * r = confp->CommonList("PhaseSpaceCuts", "Default"); + + if(!r) { + LOG("KPhaseSpaceCuts", pWARN) + << "CommonPhaseSpaceCuts.xml [Default] was not found"; + return; + } + + if(r->Exists("EM-Q2-min")) { + fEMQ2Min = this->ValidateQ2Min("EM-Q2-min", r->GetDouble("EM-Q2-min")); + fHasEMQ2Min = true; + } + + if(r->Exists("Weak-Q2-min")) { + fWeakQ2Min = this->ValidateQ2Min("Weak-Q2-min", r->GetDouble("Weak-Q2-min")); + fHasWeakQ2Min = true; + } +} +//____________________________________________________________________________ +double KPhaseSpaceCuts::ValidateQ2Min(const char * name, double q2min) const +{ + if(!std::isfinite(q2min) || q2min <= 0.) { + LOG("KPhaseSpaceCuts", pFATAL) + << "Invalid " << name << " phase-space cut: " << q2min + << " GeV^2. Q2 minimum cuts must be finite and > 0."; + exit(78); + } + return q2min; +} +//____________________________________________________________________________ +double KPhaseSpaceCuts::RaiseDefault(double default_q2min, double cut_q2min) const +{ + return TMath::Max(default_q2min, cut_q2min); +} +//____________________________________________________________________________ +void KPhaseSpaceCuts::FailMissingEMQ2Min(void) const +{ + LOG("KPhaseSpaceCuts", pFATAL) + << "Electromagnetic interactions require EM-Q2-min in " + << "CommonPhaseSpaceCuts.xml [Default], or a command-line --q2-min override."; + exit(78); +} +//____________________________________________________________________________ diff --git a/src/Framework/Interaction/KPhaseSpaceCuts.h b/src/Framework/Interaction/KPhaseSpaceCuts.h new file mode 100644 index 0000000000..dcdef0a6a8 --- /dev/null +++ b/src/Framework/Interaction/KPhaseSpaceCuts.h @@ -0,0 +1,79 @@ +//____________________________________________________________________________ +/*! + +\class genie::KPhaseSpaceCuts + +\brief Singleton service for configurable kinematic phase-space cuts. + +\author GENIE Collaboration + +\created May 21, 2026 + +\cpright Copyright (c) 2003-2025, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org +*/ +//____________________________________________________________________________ + +#ifndef _KINEMATIC_PHASE_SPACE_CUTS_H_ +#define _KINEMATIC_PHASE_SPACE_CUTS_H_ + +#include + +using std::string; + +namespace genie { + +class Interaction; + +class KPhaseSpaceCuts { + +public: + static KPhaseSpaceCuts * Instance(void); + + void SetQ2MinOverride(double q2min); + bool HasQ2MinCut(const Interaction * interaction) const; + double Q2MinCut(const Interaction * interaction, double default_q2min) const; + double EMQ2MinCut(void) const; + bool HasSplineQ2MinCut(void) const; + double SplineQ2MinCut(void) const; + string SplineQ2MinCutSource(void) const; + bool HasSplineQ2MinCutForProbe(int probe_pdg) const; + double SplineQ2MinCutForProbe(int probe_pdg) const; + string SplineQ2MinCutSourceForProbe(int probe_pdg) const; + +private: + KPhaseSpaceCuts(); + KPhaseSpaceCuts(const KPhaseSpaceCuts & cuts); + virtual ~KPhaseSpaceCuts(); + + void LoadConfig(void) const; + double ValidateQ2Min(const char * name, double q2min) const; + double RaiseDefault(double default_q2min, double cut_q2min) const; + void FailMissingEMQ2Min(void) const; + + static KPhaseSpaceCuts * fInstance; + + mutable bool fLoaded; + mutable bool fHasEMQ2Min; + mutable bool fHasWeakQ2Min; + mutable double fEMQ2Min; + mutable double fWeakQ2Min; + + bool fHasQ2MinOverride; + double fQ2MinOverride; + + struct Cleaner { + void DummyMethodAndSilentCompiler() { } + ~Cleaner() { + if (KPhaseSpaceCuts::fInstance != 0) { + delete KPhaseSpaceCuts::fInstance; + KPhaseSpaceCuts::fInstance = 0; + } + } + }; + friend struct Cleaner; +}; + +} // genie namespace + +#endif // _KINEMATIC_PHASE_SPACE_CUTS_H_ diff --git a/src/Framework/Utils/KineUtils.cxx b/src/Framework/Utils/KineUtils.cxx index 38ad3fa9b0..05ceaed4fd 100644 --- a/src/Framework/Utils/KineUtils.cxx +++ b/src/Framework/Utils/KineUtils.cxx @@ -567,7 +567,7 @@ Range1D_t genie::utils::kinematics::electromagnetic::InelWLim(double El, double } //____________________________________________________________________________ Range1D_t genie::utils::kinematics::electromagnetic::InelQ2Lim_W( - double El, double ml, double M, double W) + double El, double ml, double M, double W, double Q2min_cut) { // Computes Q2 limits (>0) @ the input W for inelastic em interactions @@ -598,18 +598,18 @@ Range1D_t genie::utils::kinematics::electromagnetic::InelQ2Lim_W( Q2.min = TMath::Max(0., Q2.min); // limit the minimum Q2 - if(Q2.min < utils::kinematics::electromagnetic::kMinQ2Limit) {Q2.min = utils::kinematics::electromagnetic::kMinQ2Limit; } // use the relevant threshold for em scattering + if(Q2.min < Q2min_cut) {Q2.min = Q2min_cut; } if(Q2.max < Q2.min ) {Q2.min = -1; Q2.max = -1;} return Q2; } //____________________________________________________________________________ Range1D_t genie::utils::kinematics::electromagnetic::Inelq2Lim_W( - double El, double ml, double M, double W) + double El, double ml, double M, double W, double q2min_cut) { // Computes q2 (<0) limits @ the input W for inelastic em interactions - Range1D_t Q2 = utils::kinematics::electromagnetic::InelQ2Lim_W(El,ml,M,W); + Range1D_t Q2 = utils::kinematics::electromagnetic::InelQ2Lim_W(El,ml,M,W,-1.*q2min_cut); Range1D_t q2; q2.min = - Q2.max; q2.max = - Q2.min; @@ -617,7 +617,7 @@ Range1D_t genie::utils::kinematics::electromagnetic::Inelq2Lim_W( } //____________________________________________________________________________ Range1D_t genie::utils::kinematics::electromagnetic::InelQ2Lim( - double El, double ml, double M) + double El, double ml, double M, double Q2min_cut) { // Computes Q2 (>0) limits irrespective of W for inelastic em interactions @@ -628,16 +628,16 @@ Range1D_t genie::utils::kinematics::electromagnetic::InelQ2Lim( Range1D_t W = utils::kinematics::electromagnetic::InelWLim(El,ml,M); if(W.min<0) return Q2; - Q2 = utils::kinematics::electromagnetic::InelQ2Lim_W(El,ml,M,W.min); + Q2 = utils::kinematics::electromagnetic::InelQ2Lim_W(El,ml,M,W.min,Q2min_cut); return Q2; } //____________________________________________________________________________ Range1D_t genie::utils::kinematics::electromagnetic::Inelq2Lim( - double El, double ml, double M) + double El, double ml, double M, double q2min_cut) { -// Computes Q2 (>0) limits irrespective of W for inelastic em interactions +// Computes q2 (<0) limits irrespective of W for inelastic em interactions - Range1D_t Q2 = utils::kinematics::electromagnetic::InelQ2Lim(El,ml,M); + Range1D_t Q2 = utils::kinematics::electromagnetic::InelQ2Lim(El,ml,M,-1.*q2min_cut); Range1D_t q2; q2.min = - Q2.max; q2.max = - Q2.min; diff --git a/src/Framework/Utils/KineUtils.h b/src/Framework/Utils/KineUtils.h index 74c400e805..ac35060874 100644 --- a/src/Framework/Utils/KineUtils.h +++ b/src/Framework/Utils/KineUtils.h @@ -102,16 +102,18 @@ namespace kinematics namespace electromagnetic { + // Legacy default Q2 threshold for em scattering events (GeV^2). + // Configurable phase-space cuts are managed by KPhaseSpaceCuts. + static const double kMinQ2Limit = 0.02; // GeV^2 + Range1D_t InelWLim (double El, double ml, double M); - Range1D_t InelQ2Lim_W (double El, double ml, double M, double W); - Range1D_t Inelq2Lim_W (double El, double ml, double M, double W); - Range1D_t InelQ2Lim (double El, double ml, double M); - Range1D_t Inelq2Lim (double El, double ml, double M); + Range1D_t InelQ2Lim_W (double El, double ml, double M, double W, double Q2min_cut = kMinQ2Limit); + Range1D_t Inelq2Lim_W (double El, double ml, double M, double W, double q2min_cut = -1*kMinQ2Limit); + Range1D_t InelQ2Lim (double El, double ml, double M, double Q2min_cut = kMinQ2Limit); + Range1D_t Inelq2Lim (double El, double ml, double M, double q2min_cut = -1*kMinQ2Limit); Range1D_t InelXLim (double El, double ml, double M); Range1D_t InelYLim (double El, double ml, double M); Range1D_t InelYLim_X (double El, double ml, double M, double x); - - static const double kMinQ2Limit = 0.02; // GeV^2 // Q2 threshold relevant for em scattering events } } // kinematics namespace diff --git a/src/Framework/Utils/XSecSplineList.cxx b/src/Framework/Utils/XSecSplineList.cxx index 66686cb267..042e383c61 100644 --- a/src/Framework/Utils/XSecSplineList.cxx +++ b/src/Framework/Utils/XSecSplineList.cxx @@ -15,6 +15,7 @@ #include #include +#include #include "libxml/parser.h" #include "libxml/xmlmemory.h" @@ -336,6 +337,38 @@ void XSecSplineList::SetMaxE(double Ev) if(Ev>0) fEmax = Ev; } //____________________________________________________________________________ +void XSecSplineList::SetTuneQ2MinKinematics( + const string & tune, double q2min, const string & source) +{ + if(tune.empty()) return; + + KinematicsMetadata & meta = fKinematicsMetadata[tune]; + meta.has_q2_min = true; + meta.q2_min = q2min; + meta.q2_unit = "GeV^2"; + meta.q2_source = source; +} +//____________________________________________________________________________ +bool XSecSplineList::HasTuneQ2MinKinematics(const string & tune) const +{ + map::const_iterator it = + fKinematicsMetadata.find(tune); + return it != fKinematicsMetadata.end() && it->second.has_q2_min; +} +//____________________________________________________________________________ +bool XSecSplineList::GetTuneQ2MinKinematics( + const string & tune, double & q2min, string * unit, string * source) const +{ + map::const_iterator it = + fKinematicsMetadata.find(tune); + if(it == fKinematicsMetadata.end() || !it->second.has_q2_min) return false; + + q2min = it->second.q2_min; + if(unit) *unit = it->second.q2_unit; + if(source) *source = it->second.q2_source; + return true; +} +//____________________________________________________________________________ void XSecSplineList::SaveAsXml(const string & filename, bool save_init) const { //! Save XSecSplineList to XML file @@ -367,6 +400,18 @@ void XSecSplineList::SaveAsXml(const string & filename, bool save_init) const outxml << " "; outxml << endl << endl; + map::const_iterator kiter = + fKinematicsMetadata.find(tune_name); + if(kiter != fKinematicsMetadata.end() && kiter->second.has_q2_min) { + outxml << " " << endl; + outxml << " second.q2_min + << "\" unit=\"" << kiter->second.q2_unit + << "\" source=\"" << kiter->second.q2_source + << "\"/>" << endl; + outxml << " " << endl << endl; + } + // loop over splines for given tune const map & spl_map_curr_tune = mm_iter->second; map::const_iterator //\/ @@ -411,7 +456,10 @@ XmlParserStatus_t XSecSplineList::LoadFromXml(const string & filename, bool keep << "Option to keep pre-existing splines is switched " << ( (keep) ? "ON" : "OFF" ); - if(!keep) fSplineMap.clear(); + if(!keep) { + fSplineMap.clear(); + fKinematicsMetadata.clear(); + } const int kNodeTypeStartElement = 1; const int kNodeTypeEndElement = 15; @@ -424,6 +472,7 @@ XmlParserStatus_t XSecSplineList::LoadFromXml(const string & filename, bool keep double * E = 0, * xsec = 0; string spline_name = ""; string temp_tune ; + bool in_kinematics = false; reader = xmlNewTextReaderFilename(filename.c_str()); if (reader != NULL) { @@ -464,6 +513,42 @@ XmlParserStatus_t XSecSplineList::LoadFromXml(const string & filename, bool keep xmlFree(xtune); } + if( (!xmlStrcmp(name, (const xmlChar *) "kinematics")) && type==kNodeTypeStartElement) { + in_kinematics = true; + } + + if( in_kinematics && (!xmlStrcmp(name, (const xmlChar *) "q2")) && + type==kNodeTypeStartElement) { + xmlChar * xmin = xmlTextReaderGetAttribute(reader,(const xmlChar*)"min"); + xmlChar * xunit = xmlTextReaderGetAttribute(reader,(const xmlChar*)"unit"); + xmlChar * xsource = xmlTextReaderGetAttribute(reader,(const xmlChar*)"source"); + + if(xmin) { + string smin = utils::str::TrimSpaces((const char *)xmin); + double q2min = atof(smin.c_str()); + if(std::isfinite(q2min) && q2min > 0.) { + KinematicsMetadata & meta = fKinematicsMetadata[temp_tune]; + meta.has_q2_min = true; + meta.q2_min = q2min; + meta.q2_unit = xunit ? + utils::str::TrimSpaces((const char *)xunit) : "GeV^2"; + meta.q2_source = xsource ? + utils::str::TrimSpaces((const char *)xsource) : ""; + SLOG("XSecSplLst", pNOTICE) + << "Loaded Q2 minimum metadata for tune " << temp_tune + << ": " << meta.q2_min << " " << meta.q2_unit; + } else { + SLOG("XSecSplLst", pWARN) + << "Ignoring invalid Q2 minimum metadata in " << filename + << ": " << smin; + } + } + + if(xmin) xmlFree(xmin); + if(xunit) xmlFree(xunit); + if(xsource) xmlFree(xsource); + } + if( (!xmlStrcmp(name, (const xmlChar *) "spline")) && type==kNodeTypeStartElement) { xmlChar * xname = xmlTextReaderGetAttribute(reader,(const xmlChar*)"name"); xmlChar * xnkn = xmlTextReaderGetAttribute(reader,(const xmlChar*)"nknots"); @@ -518,6 +603,9 @@ XmlParserStatus_t XSecSplineList::LoadFromXml(const string & filename, bool keep map::value_type(spline_name, spline) ); fLoadedSplineSet[temp_tune].insert(spline_name); } + if( (!xmlStrcmp(name, (const xmlChar *) "kinematics")) && type==kNodeTypeEndElement) { + in_kinematics = false; + } xmlFree(name); xmlFree(value); ret = xmlTextReaderRead(reader); diff --git a/src/Framework/Utils/XSecSplineList.h b/src/Framework/Utils/XSecSplineList.h index 47604d4b67..1b59883428 100644 --- a/src/Framework/Utils/XSecSplineList.h +++ b/src/Framework/Utils/XSecSplineList.h @@ -62,6 +62,17 @@ class XSecSplineList { void SetCurrentTune (const string & tune) { fCurrentTune = tune; } string CurrentTune (void) const { return fCurrentTune; } bool HasSplineFromTune( const string & tune ) const { return fSplineMap.count(tune) > 0 ; } + void SetTuneQ2MinKinematics(const string & tune, double q2min, const string & source); + bool HasTuneQ2MinKinematics(const string & tune) const; + bool GetTuneQ2MinKinematics(const string & tune, double & q2min, + string * unit = 0, string * source = 0) const; + void SetCurrentTuneQ2MinKinematics(double q2min, const string & source) + { this->SetTuneQ2MinKinematics(fCurrentTune, q2min, source); } + bool HasCurrentTuneQ2MinKinematics(void) const + { return this->HasTuneQ2MinKinematics(fCurrentTune); } + bool GetCurrentTuneQ2MinKinematics(double & q2min, + string * unit = 0, string * source = 0) const + { return this->GetTuneQ2MinKinematics(fCurrentTune, q2min, unit, source); } // Query the existence, access or create a spline // The results of the following methods depend on the current tune setting @@ -105,8 +116,17 @@ class XSecSplineList { string fCurrentTune; ///< The `active' tune, out the many that can co-exist + struct KinematicsMetadata { + KinematicsMetadata() : has_q2_min(false), q2_min(-1.), q2_unit("GeV^2"), q2_source("") { } + bool has_q2_min; + double q2_min; + string q2_unit; + string q2_source; + }; + map > fSplineMap; ///< tune -> { xsec_alg/xsec_config/interaction -> Spline } map > fLoadedSplineSet; ///< tune -> { set of initialy loaded splines } + map fKinematicsMetadata; ///< tune -> optional kinematic cut metadata } struct Cleaner { void DummyMethodAndSilentCompiler() { } diff --git a/src/Physics/DeepInelastic/EventGen/DISKinematicsGenerator.cxx b/src/Physics/DeepInelastic/EventGen/DISKinematicsGenerator.cxx index 82d7459c23..499df17593 100644 --- a/src/Physics/DeepInelastic/EventGen/DISKinematicsGenerator.cxx +++ b/src/Physics/DeepInelastic/EventGen/DISKinematicsGenerator.cxx @@ -23,6 +23,7 @@ #include "Framework/EventGen/RunningThreadInfo.h" #include "Framework/GHEP/GHepRecord.h" #include "Framework/GHEP/GHepFlags.h" +#include "Framework/Interaction/KPhaseSpace.h" #include "Framework/Messenger/Messenger.h" #include "Framework/Numerical/RandomGen.h" #include "Framework/Numerical/MathUtils.h" diff --git a/src/Physics/Multinucleon/EventGen/MECGenerator.cxx b/src/Physics/Multinucleon/EventGen/MECGenerator.cxx index 35907ab57b..9e47a03d19 100644 --- a/src/Physics/Multinucleon/EventGen/MECGenerator.cxx +++ b/src/Physics/Multinucleon/EventGen/MECGenerator.cxx @@ -21,6 +21,7 @@ #include "Framework/Conventions/Constants.h" #include "Framework/Conventions/Controls.h" #include "Framework/EventGen/EVGThreadException.h" +#include "Framework/Interaction/KPhaseSpaceCuts.h" #include "Framework/EventGen/RunningThreadInfo.h" #include "Framework/EventGen/EventGeneratorI.h" #include "Framework/GHEP/GHepStatus.h" @@ -863,9 +864,9 @@ void MECGenerator::SelectSuSALeptonKinematics(GHepRecord* event) const // Choose the appropriate minimum Q^2 value based on the interaction // mode (this is important for EM interactions since the differential // cross section blows up as Q^2 --> 0) - double Q2min = genie::controls::kMinQ2Limit; // CC/NC limit - if ( interaction->ProcInfo().IsEM() ) Q2min = genie::utils::kinematics - ::electromagnetic::kMinQ2Limit; // EM limit + double Q2min = KPhaseSpaceCuts::Instance()->Q2MinCut( + interaction, interaction->ProcInfo().IsEM() ? + 0. : genie::controls::kMinQ2Limit); LOG("MEC", pDEBUG) << "Q2min = " << Q2min; diff --git a/src/Physics/Multinucleon/XSection/EmpiricalMECPXSec2015.cxx b/src/Physics/Multinucleon/XSection/EmpiricalMECPXSec2015.cxx index 1b1dfa24fe..f9a866b478 100644 --- a/src/Physics/Multinucleon/XSection/EmpiricalMECPXSec2015.cxx +++ b/src/Physics/Multinucleon/XSection/EmpiricalMECPXSec2015.cxx @@ -19,6 +19,7 @@ #include "Framework/Conventions/GBuild.h" #include "Framework/Conventions/Units.h" #include "Framework/GHEP/GHepParticle.h" +#include "Framework/Interaction/KPhaseSpaceCuts.h" #include "Framework/Messenger/Messenger.h" #include "Physics/Multinucleon/XSection/EmpiricalMECPXSec2015.h" #include "Framework/ParticleData/PDGCodes.h" @@ -124,8 +125,11 @@ double EmpiricalMECPXSec2015::XSec( {double xsec = 0.; return xsec; } - //use proper Q2 limit from Controls.h - Range1D_t Q2lim = isem ? genie::utils::kinematics::electromagnetic::InelQ2Lim_W(Ev, ml, M2n, W) : genie::utils::kinematics::InelQ2Lim_W (Ev, M2n, ml, W, kMinQ2Limit); + double Q2min = KPhaseSpaceCuts::Instance()->Q2MinCut( + interaction, isem ? 0. : kMinQ2Limit); + Range1D_t Q2lim = isem ? + genie::utils::kinematics::electromagnetic::InelQ2Lim_W(Ev, ml, M2n, W, Q2min) : + genie::utils::kinematics::InelQ2Lim_W (Ev, M2n, ml, W, Q2min); //LOG("MEC", pINFO) << "Q2lim= " << Q2lim.min << " " < Q2lim.max) diff --git a/src/Physics/Multinucleon/XSection/MECUtils.cxx b/src/Physics/Multinucleon/XSection/MECUtils.cxx index 4df7d50439..0362e0c5b5 100644 --- a/src/Physics/Multinucleon/XSection/MECUtils.cxx +++ b/src/Physics/Multinucleon/XSection/MECUtils.cxx @@ -16,6 +16,7 @@ #include "Framework/Conventions/Controls.h" #include "Framework/Utils/KineUtils.h" #include "Framework/Interaction/Interaction.h" +#include "Framework/Interaction/KPhaseSpaceCuts.h" #include "Physics/HadronTensors/HadronTensorModelI.h" #include "Physics/Multinucleon/XSection/MECUtils.h" #include "Framework/EventGen/XSecAlgorithmI.h" @@ -342,9 +343,9 @@ double genie::utils::mec::GetMaxXSecTlctl( const XSecAlgorithmI& xsec_model, // Choose the appropriate minimum Q^2 value based on the interaction // mode (this is important for EM interactions since the differential // cross section blows up as Q^2 --> 0) - double Q2min = genie::controls::kMinQ2Limit; // CC/NC limit - if ( interaction->ProcInfo().IsEM() ) Q2min = genie::utils::kinematics - ::electromagnetic::kMinQ2Limit; // EM limit + double Q2min = KPhaseSpaceCuts::Instance()->Q2MinCut( + interaction, interaction->ProcInfo().IsEM() ? + 0. : genie::controls::kMinQ2Limit); const double Enu = interaction->InitState().ProbeE( kRfLab ); const double ProbeMass = interaction->InitState().Probe()->Mass(); @@ -524,4 +525,3 @@ genie::utils::mec::gsl::d2Xsec_dTCosth::Clone() const new genie::utils::mec::gsl::d2Xsec_dTCosth(fModel,fInteraction, fEnu, fLepMass, fFactor ); } //____________________________________________________________________________ - diff --git a/src/Physics/Multinucleon/XSection/SuSAv2MECPXSec.cxx b/src/Physics/Multinucleon/XSection/SuSAv2MECPXSec.cxx index 80c2dd54bf..6c98957078 100644 --- a/src/Physics/Multinucleon/XSection/SuSAv2MECPXSec.cxx +++ b/src/Physics/Multinucleon/XSection/SuSAv2MECPXSec.cxx @@ -9,6 +9,7 @@ //_________________________________________________________________________ #include "Framework/Algorithm/AlgConfigPool.h" +#include "Framework/Interaction/KPhaseSpaceCuts.h" #include "Framework/Messenger/Messenger.h" #include "Framework/ParticleData/PDGCodes.h" #include "Framework/ParticleData/PDGLibrary.h" @@ -115,9 +116,9 @@ double SuSAv2MECPXSec::XSec(const Interaction* interaction, // Choose the appropriate minimum Q^2 value based on the interaction // mode (this is important for EM interactions since the differential // cross section blows up as Q^2 --> 0) - double Q2min = genie::controls::kMinQ2Limit; // CC/NC limit - if ( interaction->ProcInfo().IsEM() ) Q2min = genie::utils::kinematics - ::electromagnetic::kMinQ2Limit; // EM limit + double Q2min = KPhaseSpaceCuts::Instance()->Q2MinCut( + interaction, interaction->ProcInfo().IsEM() ? + 0. : genie::controls::kMinQ2Limit); // Neglect shift due to binding energy. The cut is on the actual // value of Q^2, not the effective one to use in the tensor contraction. diff --git a/src/Physics/QuasiElastic/EventGen/QELEventGeneratorSuSA.cxx b/src/Physics/QuasiElastic/EventGen/QELEventGeneratorSuSA.cxx index c634b77193..1a4d5605d7 100644 --- a/src/Physics/QuasiElastic/EventGen/QELEventGeneratorSuSA.cxx +++ b/src/Physics/QuasiElastic/EventGen/QELEventGeneratorSuSA.cxx @@ -18,6 +18,7 @@ #include "Framework/Conventions/Constants.h" #include "Framework/Conventions/KineVar.h" #include "Framework/Conventions/KinePhaseSpace.h" +#include "Framework/Interaction/KPhaseSpaceCuts.h" #include "Framework/EventGen/EVGThreadException.h" #include "Framework/EventGen/EventGeneratorI.h" #include "Framework/EventGen/RunningThreadInfo.h" @@ -98,9 +99,9 @@ void QELEventGeneratorSuSA::SelectLeptonKinematics (GHepRecord * event) const // Choose the appropriate minimum Q^2 value based on the interaction // mode (this is important for EM interactions since the differential // cross section blows up as Q^2 --> 0) - double Q2min = genie::controls::kMinQ2Limit; // CC/NC limit - if ( interaction->ProcInfo().IsEM() ) Q2min = genie::utils::kinematics - ::electromagnetic::kMinQ2Limit; // EM limit + double Q2min = KPhaseSpaceCuts::Instance()->Q2MinCut( + interaction, interaction->ProcInfo().IsEM() ? + 0. : genie::controls::kMinQ2Limit); // The SuSA 1p1h model kinematics works in a system where // the whole nuclear target system has no momentum. diff --git a/src/Physics/QuasiElastic/XSection/SuSAv2QELPXSec.cxx b/src/Physics/QuasiElastic/XSection/SuSAv2QELPXSec.cxx index 588f59f157..09454c79b4 100644 --- a/src/Physics/QuasiElastic/XSection/SuSAv2QELPXSec.cxx +++ b/src/Physics/QuasiElastic/XSection/SuSAv2QELPXSec.cxx @@ -10,6 +10,7 @@ #include "Framework/Algorithm/AlgConfigPool.h" #include "Framework/Conventions/Units.h" +#include "Framework/Interaction/KPhaseSpaceCuts.h" #include "Framework/Messenger/Messenger.h" #include "Framework/ParticleData/PDGCodes.h" #include "Framework/ParticleData/PDGLibrary.h" @@ -63,9 +64,9 @@ double SuSAv2QELPXSec::XSec(const Interaction* interaction, // Choose the appropriate minimum Q^2 value based on the interaction // mode (this is important for EM interactions since the differential // cross section blows up as Q^2 --> 0) - double Q2min = genie::controls::kMinQ2Limit; // CC/NC limit - if ( interaction->ProcInfo().IsEM() ) Q2min = genie::utils::kinematics - ::electromagnetic::kMinQ2Limit; // EM limit + double Q2min = KPhaseSpaceCuts::Instance()->Q2MinCut( + interaction, interaction->ProcInfo().IsEM() ? + 0. : genie::controls::kMinQ2Limit); // Neglect shift due to binding energy. The cut is on the actual // value of Q^2, not the effective one to use in the tensor contraction. diff --git a/src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx b/src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx index ce2e53c27b..19212fd2bc 100644 --- a/src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx +++ b/src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx @@ -22,6 +22,8 @@ #include "Framework/EventGen//RunningThreadInfo.h" #include "Framework/GHEP/GHepRecord.h" #include "Framework/GHEP/GHepFlags.h" +#include "Framework/Interaction/KPhaseSpace.h" +#include "Framework/Interaction/KPhaseSpaceCuts.h" #include "Framework/Messenger/Messenger.h" #include "Framework/Numerical/RandomGen.h" #include "Framework/Numerical/MathUtils.h" @@ -285,8 +287,9 @@ double RESKinematicsGenerator::ComputeMaxXSec( const InitialState & init_state = interaction -> InitState(); double E = init_state.ProbeE(kRfHitNucRest); - bool is_em = interaction->ProcInfo().IsEM(); - double Q2Thres = is_em ? utils::kinematics::electromagnetic::kMinQ2Limit : controls::kMinQ2Limit; + double Q2Thres = KPhaseSpaceCuts::Instance()->Q2MinCut( + interaction, interaction->ProcInfo().IsEM() ? + 0. : controls::kMinQ2Limit); double md; if(!interaction->ExclTag().KnownResonance()) md=1.23; diff --git a/src/contrib/test/Makefile b/src/contrib/test/Makefile index 391caf1588..da0657dc26 100644 --- a/src/contrib/test/Makefile +++ b/src/contrib/test/Makefile @@ -14,6 +14,9 @@ include $(GENIE)/src/make/Make.include GENIE_LIBS = $(shell $(GENIE)/src/scripts/setup/genie-config --libs) LIBRARIES := $(GENIE_LIBS) $(LIBRARIES) $(CERN_LIBRARIES) +KPHASESPACE_ROOT_LIBRARIES := $(shell root-config --glibs) -lMinuit -lGeom -lEG -lGenVector -lMathMore +KPHASESPACE_LIBRARIES := -L$(GENIE_LIB_PATH) -lGFwInt -lGFwAlg -lGFwReg -lGFwUtl -lGFwNum -lGFwParDat -lGFwMsg \ + $(KPHASESPACE_ROOT_LIBRARIES) $(XML_LIBRARIES) $(LOG_LIBRARIES) #gtestDecay \ #gtestHadronization \ @@ -182,7 +185,7 @@ gtestResonances: FORCE gtestKPhaseSpace: FORCE $(CXX) $(CXXFLAGS) -c gtestKPhaseSpace.cxx $(CPP_INCLUDES) - $(LD) $(LDFLAGS) gtestKPhaseSpace.o $(LIBRARIES) -o $(GENIE_BIN_PATH)/gtestKPhaseSpace + $(LD) $(LDFLAGS) gtestKPhaseSpace.o $(KPHASESPACE_LIBRARIES) -o $(GENIE_BIN_PATH)/gtestKPhaseSpace gtestROOTGeometry: FORCE ifeq ($(strip $(GOPT_ENABLE_GEOM_DRIVERS)),YES) diff --git a/src/contrib/test/gtestKPhaseSpace.cxx b/src/contrib/test/gtestKPhaseSpace.cxx index 02b94f64cf..fc812878a5 100644 --- a/src/contrib/test/gtestKPhaseSpace.cxx +++ b/src/contrib/test/gtestKPhaseSpace.cxx @@ -19,17 +19,44 @@ #include #include +#include +#include +#include +#include +#include +#include + +#include "Framework/Conventions/Controls.h" #include "Framework/Interaction/Interaction.h" +#include "Framework/Interaction/KPhaseSpaceCuts.h" #include "Framework/Messenger/Messenger.h" #include "Framework/ParticleData/PDGCodes.h" +#include "Framework/Utils/XSecSplineList.h" using namespace genie; void PrintLimits(const Interaction * interaction); +int RunDefaultChecks(const char * self); +int RunMissingEMCutCheck(void); +int RunWeakConfigCutCheck(void); +int RunNoSplineQ2CutCheck(void); +bool Check(bool ok, const char * msg); +bool NearlyEqual(double a, double b); +int RunChild(const char * self, const char * mode, const std::string & xml); +std::string WritePhaseSpaceCutsXML(const std::string & xml_body); +std::string WriteSplineXML(const std::string & xml_body); +void RemoveTempXML(const std::string & xml); //__________________________________________________________________________ -int main(int /*argc*/, char ** /*argv*/) +int main(int argc, char ** argv) { + if(argc >= 2) { + std::string mode = argv[1]; + if(mode == "--expect-missing-em-cut") _exit(RunMissingEMCutCheck()); + if(mode == "--expect-weak-config-cut") _exit(RunWeakConfigCutCheck()); + if(mode == "--expect-no-spline-q2-cut") _exit(RunNoSplineQ2CutCheck()); + } + // -- get a DIS interaction object & access its kinematics int tgt = kPdgTgtFe56; @@ -45,7 +72,7 @@ int main(int /*argc*/, char ** /*argv*/) PrintLimits(rescc); PrintLimits(discc); - return 0; + return RunDefaultChecks(argv[0]); } //__________________________________________________________________________ void PrintLimits(const Interaction * interaction) @@ -65,3 +92,253 @@ void PrintLimits(const Interaction * interaction) LOG("test", pNOTICE) << "W e [" << Wl.min << ", " << Wl.max << "]"; } //__________________________________________________________________________ +int RunDefaultChecks(const char * self) +{ + int failures = 0; + const double em_q2_min = 0.02; + const double weak_default = 0.123; + + Interaction * em = Interaction::QELEM( + kPdgTgtFe56, kPdgProton, kPdgElectron, 3.); + Interaction * weak = Interaction::QELCC( + kPdgTgtFe56, kPdgProton, kPdgNuMu, 3.); + + KPhaseSpaceCuts * cuts = KPhaseSpaceCuts::Instance(); + + failures += !Check(cuts->HasQ2MinCut(em), + "EM interactions have a configured Q2 cut"); + failures += !Check(NearlyEqual(cuts->Q2MinCut(em, 0.), em_q2_min), + "EM Q2 cut is loaded from CommonPhaseSpaceCuts.xml"); + failures += !Check(em->PhaseSpace().Limits(kKVQ2).min >= em_q2_min, + "KPhaseSpace applies the EM Q2 cut to Q2 limits"); + + failures += !Check(!cuts->HasQ2MinCut(weak), + "weak interactions do not use a configurable Q2 cut by default"); + failures += !Check(NearlyEqual(cuts->Q2MinCut(weak, weak_default), weak_default), + "weak interactions keep the caller default without an explicit cut"); + failures += !Check(!cuts->HasSplineQ2MinCutForProbe(kPdgNuMu), + "neutrino spline metadata ignores EM-only Q2 cuts"); + failures += !Check(cuts->HasSplineQ2MinCutForProbe(kPdgElectron), + "electron spline metadata sees EM Q2 cuts"); + failures += !Check(NearlyEqual(cuts->SplineQ2MinCutForProbe(kPdgElectron), em_q2_min), + "electron spline metadata uses EM-Q2-min"); + + cuts->SetQ2MinOverride(0.25); + failures += !Check(cuts->HasQ2MinCut(weak), + "command-line Q2 override opts weak interactions into the cut"); + failures += !Check(NearlyEqual(cuts->Q2MinCut(weak, 0.), 0.25), + "command-line Q2 override is applied to weak interactions"); + failures += !Check(cuts->HasSplineQ2MinCutForProbe(kPdgNuMu), + "command-line Q2 override provides spline metadata"); + failures += !Check(NearlyEqual(cuts->SplineQ2MinCutForProbe(kPdgNuMu), 0.25), + "spline metadata uses the command-line Q2 override"); + + delete em; + delete weak; + + std::string missing_xml = WritePhaseSpaceCutsXML( + "\n" + "\n" + " \n" + " \n" + "\n"); + if(missing_xml.empty()) { + failures += !Check(false, "temporary XML for missing EM-Q2-min was created"); + } else { + int missing_status = RunChild(self, "--expect-missing-em-cut", missing_xml); + failures += !Check(WIFEXITED(missing_status) && + WEXITSTATUS(missing_status) == 78, + "EM interactions fail loudly when EM-Q2-min is missing"); + } + RemoveTempXML(missing_xml); + + std::string weak_xml = WritePhaseSpaceCutsXML( + "\n" + "\n" + " \n" + " 0.02 \n" + " 0.25 \n" + " \n" + "\n"); + if(weak_xml.empty()) { + failures += !Check(false, "temporary XML for Weak-Q2-min was created"); + } else { + int weak_status = RunChild(self, "--expect-weak-config-cut", weak_xml); + failures += !Check(WIFEXITED(weak_status) && + WEXITSTATUS(weak_status) == 0, + "Weak-Q2-min from CommonPhaseSpaceCuts.xml opts weak interactions in"); + } + RemoveTempXML(weak_xml); + + std::string empty_xml = WritePhaseSpaceCutsXML( + "\n" + "\n" + " \n" + " \n" + "\n"); + if(empty_xml.empty()) { + failures += !Check(false, "temporary XML for empty Q2 metadata was created"); + } else { + int empty_status = RunChild(self, "--expect-no-spline-q2-cut", empty_xml); + failures += !Check(WIFEXITED(empty_status) && + WEXITSTATUS(empty_status) == 0, + "empty CommonPhaseSpaceCuts.xml produces no spline Q2 metadata"); + } + RemoveTempXML(empty_xml); + + std::string spline_xml = WriteSplineXML( + "\n" + "\n" + " \n" + " \n" + " \n" + " \n" + " \n" + " 0.10000 1.0000000000e-38 \n" + " 1.00000 2.0000000000e-38 \n" + " 2.00000 3.0000000000e-38 \n" + " \n" + " \n" + "\n"); + if(spline_xml.empty()) { + failures += !Check(false, "temporary spline XML was created"); + } else { + XSecSplineList * xspl = XSecSplineList::Instance(); + failures += !Check(xspl->LoadFromXml(spline_xml) == kXmlOK, + "spline XML with kinematics metadata loads"); + double q2min = -1.; + std::string unit = ""; + std::string source = ""; + failures += !Check(xspl->GetTuneQ2MinKinematics( + "UnitTestTune", q2min, &unit, &source), + "loaded spline XML exposes Q2 metadata"); + failures += !Check(NearlyEqual(q2min, 0.25), + "loaded spline XML preserves Q2 metadata value"); + + std::string saved_spline_xml = spline_xml + ".roundtrip.xml"; + xspl->SaveAsXml(saved_spline_xml); + failures += !Check(xspl->LoadFromXml(saved_spline_xml) == kXmlOK, + "saved spline XML with kinematics metadata reloads"); + q2min = -1.; + failures += !Check(xspl->GetTuneQ2MinKinematics( + "UnitTestTune", q2min, &unit, &source), + "round-tripped spline XML exposes Q2 metadata"); + failures += !Check(NearlyEqual(q2min, 0.25), + "round-tripped spline XML preserves Q2 metadata value"); + unlink(saved_spline_xml.c_str()); + } + RemoveTempXML(spline_xml); + + return failures == 0 ? 0 : 1; +} +//__________________________________________________________________________ +int RunMissingEMCutCheck(void) +{ + Interaction * em = Interaction::QELEM( + kPdgTgtFe56, kPdgProton, kPdgElectron, 3.); + KPhaseSpaceCuts::Instance()->Q2MinCut(em, 0.); + delete em; + return 1; +} +//__________________________________________________________________________ +int RunWeakConfigCutCheck(void) +{ + Interaction * weak = Interaction::QELCC( + kPdgTgtFe56, kPdgProton, kPdgNuMu, 3.); + KPhaseSpaceCuts * cuts = KPhaseSpaceCuts::Instance(); + int failures = 0; + failures += !Check(cuts->HasQ2MinCut(weak), + "weak interaction sees Weak-Q2-min from XML"); + failures += !Check(NearlyEqual(cuts->Q2MinCut(weak, 0.), 0.25), + "weak interaction uses Weak-Q2-min from XML"); + failures += !Check(cuts->HasSplineQ2MinCut(), + "Weak-Q2-min provides spline metadata"); + failures += !Check(NearlyEqual(cuts->SplineQ2MinCut(), 0.25), + "spline metadata uses Weak-Q2-min from XML"); + failures += !Check(cuts->HasSplineQ2MinCutForProbe(kPdgNuMu), + "Weak-Q2-min provides neutrino spline metadata"); + failures += !Check(NearlyEqual(cuts->SplineQ2MinCutForProbe(kPdgNuMu), 0.25), + "neutrino spline metadata uses Weak-Q2-min from XML"); + delete weak; + return failures == 0 ? 0 : 1; +} +//__________________________________________________________________________ +int RunNoSplineQ2CutCheck(void) +{ + KPhaseSpaceCuts * cuts = KPhaseSpaceCuts::Instance(); + int failures = 0; + failures += !Check(!cuts->HasSplineQ2MinCut(), + "empty CommonPhaseSpaceCuts.xml has no spline Q2 metadata"); + failures += !Check(!cuts->HasSplineQ2MinCutForProbe(kPdgNuMu), + "empty CommonPhaseSpaceCuts.xml has no neutrino spline Q2 metadata"); + failures += !Check(!cuts->HasSplineQ2MinCutForProbe(kPdgElectron), + "empty CommonPhaseSpaceCuts.xml has no electron spline Q2 metadata"); + return failures == 0 ? 0 : 1; +} +//__________________________________________________________________________ +bool Check(bool ok, const char * msg) +{ + if(ok) LOG("test", pNOTICE) << "PASS: " << msg; + else LOG("test", pERROR) << "FAIL: " << msg; + return ok; +} +//__________________________________________________________________________ +bool NearlyEqual(double a, double b) +{ + return std::fabs(a - b) < 1e-12; +} +//__________________________________________________________________________ +int RunChild(const char * self, const char * mode, const std::string & xml) +{ + size_t slash = xml.rfind('/'); + std::string dir = slash == std::string::npos ? "." : xml.substr(0, slash); + + pid_t pid = fork(); + if(pid == 0) { + setenv("GXMLPATH", dir.c_str(), 1); + execlp(self, self, mode, static_cast(0)); + _exit(127); + } + + int status = 0; + if(pid < 0 || waitpid(pid, &status, 0) < 0) return -1; + return status; +} +//__________________________________________________________________________ +std::string WritePhaseSpaceCutsXML(const std::string & xml_body) +{ + char tmp[] = "/tmp/genie-kphase-space-cuts-XXXXXX"; + char * dir = mkdtemp(tmp); + if(!dir) return ""; + + std::string xml = std::string(dir) + "/CommonPhaseSpaceCuts.xml"; + std::ofstream out(xml.c_str()); + out << xml_body; + out.close(); + return xml; +} +//__________________________________________________________________________ +std::string WriteSplineXML(const std::string & xml_body) +{ + char tmp[] = "/tmp/genie-xsec-spline-XXXXXX"; + char * dir = mkdtemp(tmp); + if(!dir) return ""; + + std::string xml = std::string(dir) + "/spline.xml"; + std::ofstream out(xml.c_str()); + out << xml_body; + out.close(); + return xml; +} +//__________________________________________________________________________ +void RemoveTempXML(const std::string & xml) +{ + if(xml.empty()) return; + unlink(xml.c_str()); + size_t slash = xml.rfind('/'); + if(slash != std::string::npos) { + std::string dir = xml.substr(0, slash); + rmdir(dir.c_str()); + } +} +//__________________________________________________________________________