From 96725594ca85a74dc587874c8f5504f7906e431a Mon Sep 17 00:00:00 2001 From: candreop Date: Fri, 20 Oct 2023 12:12:01 +0100 Subject: [PATCH 01/68] add skeleton for BY2021 (copy of BY) to start tweaking w/o affecting the BY model --- config/master_config.xml | 1 + .../XSection/BYStrucFunc2021.cxx | 121 ++++++++++++++++++ .../DeepInelastic/XSection/BYStrucFunc2021.h | 69 ++++++++++ src/Physics/DeepInelastic/XSection/LinkDef.h | 1 + 4 files changed, 192 insertions(+) create mode 100644 src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx create mode 100644 src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h diff --git a/config/master_config.xml b/config/master_config.xml index 048508c3f1..1e74c8c02a 100644 --- a/config/master_config.xml +++ b/config/master_config.xml @@ -141,6 +141,7 @@ LwlynSmithFFNC.xml QPMDISStrucFunc.xml BYStrucFunc.xml + BYStrucFunc2021.xml RSHelicityAmplModelCC.xml RSHelicityAmplModelNCp.xml RSHelicityAmplModelNCn.xml diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx new file mode 100644 index 0000000000..7ca87bd7a8 --- /dev/null +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx @@ -0,0 +1,121 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + Júlia Tena Vidal + Tel Aviv University + + Costas Andreopoulos + University of Liverpool + +*/ +//____________________________________________________________________________ + +#include + +#include "Framework/Algorithm/AlgFactory.h" +#include "Framework/Algorithm/AlgConfigPool.h" +#include "Physics/DeepInelastic/XSection/BYStrucFunc2021.h" +#include "Framework/Conventions/Constants.h" +#include "Framework/Messenger/Messenger.h" + +using namespace genie; +using namespace genie::constants; + +//____________________________________________________________________________ +BYStrucFunc2021::BYStrucFunc2021() : +QPMDISStrucFuncBase("genie::BYStrucFunc2021") +{ + this->Init(); +} +//____________________________________________________________________________ +BYStrucFunc2021::BYStrucFunc2021(string config): +QPMDISStrucFuncBase("genie::BYStrucFunc2021", config) +{ + this->Init(); +} +//____________________________________________________________________________ +BYStrucFunc2021::~BYStrucFunc2021() +{ + +} +//____________________________________________________________________________ +void BYStrucFunc2021::Configure(const Registry & config) +{ +// Overload Algorithm::Configure() to read the config. registry and set +// private data members. +// QPMDISStrucFuncBase::Configure() creates the owned PDF object that gets +// configured with the specified PDFModelI +// For the ReadBYParams() method see below + + QPMDISStrucFuncBase::Configure(config); + this->ReadBYParams(); +} +//____________________________________________________________________________ +void BYStrucFunc2021::Configure(string param_set) +{ + QPMDISStrucFuncBase::Configure(param_set); + this->ReadBYParams(); +} +//____________________________________________________________________________ +void BYStrucFunc2021::ReadBYParams(void) +{ +// Get the Bodek-Yang model parameters A,B,Csea,Cv1,Cv2 from the config. +// registry and set some private data members so as not to accessing the +// registry at every calculation. +// + GetParam( "BY-A", fA ) ; + GetParam( "BY-B", fB ) ; + GetParam( "BY-CsU", fCsU ) ; + GetParam( "BY-CsD", fCsD ) ; + GetParam( "BY-Cv1U", fCv1U ) ; + GetParam( "BY-Cv2U", fCv2U ) ; + GetParam( "BY-Cv1D", fCv1D ) ; + GetParam( "BY-Cv2D", fCv2D ) ; + +} +//____________________________________________________________________________ +void BYStrucFunc2021::Init(void) +{ + fA = 0; + fB = 0; + fCsU = 0; + fCsD = 0; + fCv1U = 0; + fCv2U = 0; + fCv1D = 0; + fCv2D = 0; +} +//____________________________________________________________________________ +double BYStrucFunc2021::ScalingVar(const Interaction * interaction) const +{ +// Overrides QPMDISStrucFuncBase::ScalingVar() to compute the BY scaling var + + const Kinematics & kine = interaction->Kine(); + double x = kine.x(); + double myQ2 = this->Q2(interaction); + //myQ2 = TMath::Max(Q2,fQ2min); + LOG("BodekYang", pDEBUG) << "Q2 at scaling var calculation = " << myQ2; + + double a = TMath::Power( 2*kProtonMass*x, 2 ) / myQ2; + double xw = 2*x*(myQ2+fB) / (myQ2*(1.+TMath::Sqrt(1+a)) + 2*fA*x); + return xw; +} +//____________________________________________________________________________ +void BYStrucFunc2021::KFactors(const Interaction * interaction, + double & kuv, double & kdv, double & kus, double & kds) const +{ +// Overrides QPMDISStrucFuncBase::KFactors() to compute the BY K factors for +// u(valence), d(valence), u(sea), d(sea); + + double myQ2 = this->Q2(interaction); + double GD = 1. / TMath::Power(1.+myQ2/0.71, 2); // p elastic form factor + double GD2 = TMath::Power(GD,2); + + kuv = (1.-GD2)*(myQ2+fCv2U)/(myQ2+fCv1U); // K - u(valence) + kdv = (1.-GD2)*(myQ2+fCv2D)/(myQ2+fCv1D); // K - d(valence) + kus = myQ2/(myQ2+fCsU); // K - u(sea) + kds = myQ2/(myQ2+fCsD); // K - d(sea) +} +//____________________________________________________________________________ diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h new file mode 100644 index 0000000000..b21a69dd39 --- /dev/null +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h @@ -0,0 +1,69 @@ +//____________________________________________________________________________ +/*! + +\class genie::BYStrucFunc2021 + +\brief 2021 update of the Bodek and Yang structure function model + +\ref arXiv:2108.09240v2 [hep-ph] + +\author Júlia Tena Vidal + Tel Aviv University + + Costas Andreopoulos + University of Liverpool + +\created October 20, 2023 + +\cpright Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org +*/ +//____________________________________________________________________________ + +#ifndef _BODEK_YANG_STRUCTURE_FUNCTION_MODEL_2021_H_ +#define _BODEK_YANG_STRUCTURE_FUNCTION_MODEL_2021_H_ + +#include "Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h" +#include "Framework/Interaction/Interaction.h" +#include "Physics/PartonDistributions/PDFModelI.h" + +namespace genie { + +class BYStrucFunc2021 : public QPMDISStrucFuncBase { + +public: + BYStrucFunc2021(); + BYStrucFunc2021(string config); + virtual ~BYStrucFunc2021(); + + // overload Algorithm::Configure() to read the config. registry + // at the algorithm initialization and set private data members + void Configure (const Registry & config); + void Configure (string param_set); + +protected: + + void Init (void); + void ReadBYParams (void); + + // override part of the DISStructureFuncModel implementation + // to compute all the corrections applied by the Bodek-Yang model. + double ScalingVar (const Interaction * i) const; + void KFactors (const Interaction * i, double & kuv, + double & kdv, double & kus, double & kds) const; + + // Bodek-Yang model-specific parameters + + double fA; ///< better scaling var parameter A + double fB; ///< better scaling var parameter B + double fCsU; ///< U-sea K factor parameter + double fCsD; ///< D-sea K factor parameter + double fCv1U; ///< U-val K factor parameter + double fCv2U; ///< U-val K factor parameter + double fCv1D; ///< D-val K factor parameter + double fCv2D; ///< D-val K factor parameter +}; + +} // genie namespace + +#endif // _BODEK_YANG_STRUCTURE_FUNCTION_MODEL_2021_H_ diff --git a/src/Physics/DeepInelastic/XSection/LinkDef.h b/src/Physics/DeepInelastic/XSection/LinkDef.h index f171f01773..97eb51a6c4 100644 --- a/src/Physics/DeepInelastic/XSection/LinkDef.h +++ b/src/Physics/DeepInelastic/XSection/LinkDef.h @@ -14,6 +14,7 @@ #pragma link C++ class genie::BYPDF; #pragma link C++ class genie::BYStrucFunc; +#pragma link C++ class genie::BYStrucFunc2021; #pragma link C++ class genie::DISStructureFunc; #pragma link C++ class genie::DISStructureFuncModelI; From 66032e0844be65e1dc198d88b257b9c78bd99500 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Fri, 17 Oct 2025 14:40:17 +0200 Subject: [PATCH 02/68] adding updated parameters --- config/BYStrucFunc.xml | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/config/BYStrucFunc.xml b/config/BYStrucFunc.xml index 3d379db2f3..d73233e9cc 100644 --- a/config/BYStrucFunc.xml +++ b/config/BYStrucFunc.xml @@ -36,7 +36,7 @@ WeinbergAngle double No - + CKM,Masses,WeakInt genie::BYPDF/Default @@ -63,7 +63,21 @@ WeinbergAngle double No false 0.8 + + + + + 0.621 + 0.380 + 0.369 + 0.561 + 0.417 + 0.264 + 0.341 + 0.323 + + From f6fd906d05a52ad60ce4dff65e04e8366006f5e6 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Fri, 17 Oct 2025 17:58:43 +0200 Subject: [PATCH 03/68] adding updated BY parameters --- config/KNOTunedQPMDISPXSec.xml | 2 +- config/QPMDISPXSec.xml | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/config/KNOTunedQPMDISPXSec.xml b/config/KNOTunedQPMDISPXSec.xml index 77cc15ad96..fef4d7b460 100644 --- a/config/KNOTunedQPMDISPXSec.xml +++ b/config/KNOTunedQPMDISPXSec.xml @@ -19,7 +19,7 @@ Wcut double no Co NonResBackground - + 0.1 genie::QPMDISPXSec/Default genie::AGKYLowW2019/Default 1. diff --git a/config/QPMDISPXSec.xml b/config/QPMDISPXSec.xml index 265805a20c..e72ff32e52 100644 --- a/config/QPMDISPXSec.xml +++ b/config/QPMDISPXSec.xml @@ -26,4 +26,7 @@ WeinbergAngle double No 1. + + genie::BYStrucFunc/GRV982010 + From 52b6b410c8fa9d61fd9abd6bc8991f19fab3bb28 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Fri, 17 Oct 2025 18:06:07 +0200 Subject: [PATCH 04/68] update BY model --- config/KNOTunedQPMDISPXSec.xml | 1 - 1 file changed, 1 deletion(-) diff --git a/config/KNOTunedQPMDISPXSec.xml b/config/KNOTunedQPMDISPXSec.xml index fef4d7b460..1abf6169e8 100644 --- a/config/KNOTunedQPMDISPXSec.xml +++ b/config/KNOTunedQPMDISPXSec.xml @@ -19,7 +19,6 @@ Wcut double no Co NonResBackground - 0.1 genie::QPMDISPXSec/Default genie::AGKYLowW2019/Default 1. From ef624fe3beefc12ae70113c5471fc359f1e13993 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Mon, 20 Oct 2025 14:15:15 +0200 Subject: [PATCH 05/68] adding missing parameters and fixing K factors --- .../DeepInelastic/XSection/BYStrucFunc.cxx | 2 +- .../DeepInelastic/XSection/BYStrucFunc.h | 2 +- .../XSection/BYStrucFunc2021.cxx | 20 +++++++++++++------ .../DeepInelastic/XSection/BYStrucFunc2021.h | 4 +++- .../XSection/QPMDISStrucFuncBase.cxx | 3 ++- .../XSection/QPMDISStrucFuncBase.h | 2 +- 6 files changed, 22 insertions(+), 11 deletions(-) diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc.cxx index 2f58e1a497..b81b3719f7 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc.cxx @@ -100,7 +100,7 @@ double BYStrucFunc::ScalingVar(const Interaction * interaction) const } //____________________________________________________________________________ void BYStrucFunc::KFactors(const Interaction * interaction, - double & kuv, double & kdv, double & kus, double & kds) const + double & kuv, double & kdv, double & kus, double & kds, double & kss ) const { // Overrides QPMDISStrucFuncBase::KFactors() to compute the BY K factors for // u(valence), d(valence), u(sea), d(sea); diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc.h b/src/Physics/DeepInelastic/XSection/BYStrucFunc.h index afc58163cb..f000572923 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc.h +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc.h @@ -47,7 +47,7 @@ class BYStrucFunc : public QPMDISStrucFuncBase { // to compute all the corrections applied by the Bodek-Yang model. double ScalingVar (const Interaction * i) const; void KFactors (const Interaction * i, double & kuv, - double & kdv, double & kus, double & kds) const; + double & kdv, double & kus, double & kds, double & kss ) const; // Bodek-Yang model-specific parameters diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx index 7ca87bd7a8..fbc86c3738 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx @@ -69,11 +69,12 @@ void BYStrucFunc2021::ReadBYParams(void) GetParam( "BY-B", fB ) ; GetParam( "BY-CsU", fCsU ) ; GetParam( "BY-CsD", fCsD ) ; + GetParam( "BY-CsS", fCsD ) ; GetParam( "BY-Cv1U", fCv1U ) ; GetParam( "BY-Cv2U", fCv2U ) ; GetParam( "BY-Cv1D", fCv1D ) ; GetParam( "BY-Cv2D", fCv2D ) ; - + GetParam( "BY-CsS" , fCsS ) ; } //____________________________________________________________________________ void BYStrucFunc2021::Init(void) @@ -82,10 +83,12 @@ void BYStrucFunc2021::Init(void) fB = 0; fCsU = 0; fCsD = 0; + fCvLW = 0; fCv1U = 0; fCv2U = 0; fCv1D = 0; fCv2D = 0; + fCsS = 0; } //____________________________________________________________________________ double BYStrucFunc2021::ScalingVar(const Interaction * interaction) const @@ -104,18 +107,23 @@ double BYStrucFunc2021::ScalingVar(const Interaction * interaction) const } //____________________________________________________________________________ void BYStrucFunc2021::KFactors(const Interaction * interaction, - double & kuv, double & kdv, double & kus, double & kds) const + double & kuv, double & kdv, double & kus, double & kds, double & kss ) const { // Overrides QPMDISStrucFuncBase::KFactors() to compute the BY K factors for -// u(valence), d(valence), u(sea), d(sea); +// u(valence), d(valence), u(sea), d(sea), s(sea); double myQ2 = this->Q2(interaction); double GD = 1. / TMath::Power(1.+myQ2/0.71, 2); // p elastic form factor double GD2 = TMath::Power(GD,2); - - kuv = (1.-GD2)*(myQ2+fCv2U)/(myQ2+fCv1U); // K - u(valence) - kdv = (1.-GD2)*(myQ2+fCv2D)/(myQ2+fCv1D); // K - d(valence) + // double Ev = init_state.ProbeE(kRfHitNucRest); + double q0 = 0; + double q02 = TMath::Power(q0,2); + double KLW = q02 / ( q02 + fCvLW ); + + kuv = KLW * (1.-GD2)*(myQ2+fCv2U)/(myQ2+fCv1U); // K - u(valence) + kdv = KLW * (1.-GD2)*(myQ2+fCv2D)/(myQ2+fCv1D); // K - d(valence) kus = myQ2/(myQ2+fCsU); // K - u(sea) kds = myQ2/(myQ2+fCsD); // K - d(sea) + double kss = myQ2/(myQ2+fCsS); // K - s(sea) } //____________________________________________________________________________ diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h index b21a69dd39..1bcf158297 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h @@ -50,14 +50,16 @@ class BYStrucFunc2021 : public QPMDISStrucFuncBase { // to compute all the corrections applied by the Bodek-Yang model. double ScalingVar (const Interaction * i) const; void KFactors (const Interaction * i, double & kuv, - double & kdv, double & kus, double & kds) const; + double & kdv, double & kus, double & kds, double & kss ) const; // Bodek-Yang model-specific parameters double fA; ///< better scaling var parameter A double fB; ///< better scaling var parameter B + double fCvLW; ///< C low-nu vector paramter double fCsU; ///< U-sea K factor parameter double fCsD; ///< D-sea K factor parameter + double fCsS; ///< S-sea K factor parameter double fCv1U; ///< U-val K factor parameter double fCv2U; ///< U-val K factor parameter double fCv1D; ///< D-val K factor parameter diff --git a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx index 4f5cef8dc6..567fc437e3 100644 --- a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx +++ b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx @@ -432,7 +432,7 @@ double QPMDISStrucFuncBase::ScalingVar(const Interaction* interaction) const } //____________________________________________________________________________ void QPMDISStrucFuncBase::KFactors(const Interaction *, - double & kuv, double & kdv, double & kus, double & kds) const + double & kuv, double & kdv, double & kus, double & kds, double & kds ) const { // This is an abstract class: no model-specific correction // The PDF scaling variables are set to 1 @@ -442,6 +442,7 @@ void QPMDISStrucFuncBase::KFactors(const Interaction *, kdv = 1.; kus = 1.; kds = 1.; + kss = 1.; } //____________________________________________________________________________ double QPMDISStrucFuncBase::NuclMod(const Interaction * interaction) const diff --git a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h index ba58c99e85..54fa18c08d 100644 --- a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h +++ b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h @@ -68,7 +68,7 @@ class QPMDISStrucFuncBase : public DISStructureFuncModelI { virtual double NuclMod (const Interaction * i) const; virtual double R (const Interaction * i) const; virtual void KFactors (const Interaction * i, double & kuv, - double & kdv, double & kus, double & kds) const; + double & kdv, double & kus, double & kds, double & kss ) const; // configuration // double fQ2min; ///< min Q^2 allowed for PDFs: PDF(Q2 Date: Mon, 20 Oct 2025 14:22:16 +0200 Subject: [PATCH 06/68] adding GRV98 Second iteration parameters --- config/BYStrucFunc2021.xml | 71 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) create mode 100644 config/BYStrucFunc2021.xml diff --git a/config/BYStrucFunc2021.xml b/config/BYStrucFunc2021.xml new file mode 100644 index 0000000000..92af26824d --- /dev/null +++ b/config/BYStrucFunc2021.xml @@ -0,0 +1,71 @@ + + + + + + + + CKM,Masses,WeakInt + + genie::BYPDF/Default + + 0.621 + 0.380 + 0.218 + 0.369 + 0.561 + 0.417 + 0.264 + 0.341 + 0.323 + 0.561 + + + true + true + false + 0.8 + + + + + From 906a85beacacee05692907697289752499f814a1 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Mon, 20 Oct 2025 14:24:15 +0200 Subject: [PATCH 07/68] fixing bug in name --- src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx index fbc86c3738..320620e9b3 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx @@ -69,7 +69,7 @@ void BYStrucFunc2021::ReadBYParams(void) GetParam( "BY-B", fB ) ; GetParam( "BY-CsU", fCsU ) ; GetParam( "BY-CsD", fCsD ) ; - GetParam( "BY-CsS", fCsD ) ; + GetParam( "BY-CvLW", fCvLW ) ; GetParam( "BY-Cv1U", fCv1U ) ; GetParam( "BY-Cv2U", fCv2U ) ; GetParam( "BY-Cv1D", fCv1D ) ; From 982f6d0b66df72b9ae976a90fd32f92b639072d7 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Mon, 20 Oct 2025 14:49:00 +0200 Subject: [PATCH 08/68] fixing compilation error --- src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx index 320620e9b3..f27baf8e57 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx @@ -122,8 +122,8 @@ void BYStrucFunc2021::KFactors(const Interaction * interaction, kuv = KLW * (1.-GD2)*(myQ2+fCv2U)/(myQ2+fCv1U); // K - u(valence) kdv = KLW * (1.-GD2)*(myQ2+fCv2D)/(myQ2+fCv1D); // K - d(valence) - kus = myQ2/(myQ2+fCsU); // K - u(sea) - kds = myQ2/(myQ2+fCsD); // K - d(sea) - double kss = myQ2/(myQ2+fCsS); // K - s(sea) + kus = myQ2/(myQ2+fCsU); // K - u(sea) + kds = myQ2/(myQ2+fCsD); // K - d(sea) + kss = myQ2/(myQ2+fCsS); // K - s(sea) } //____________________________________________________________________________ From 0a7d42d6c9a263dd806f7f43a6bf6d2a0fad7140 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Tue, 21 Oct 2025 11:56:41 +0200 Subject: [PATCH 09/68] adding q0 definition --- .../XSection/BYStrucFunc2021.cxx | 4 +-- .../XSection/QPMDISStrucFuncBase.cxx | 27 ++++++++++++++++--- .../XSection/QPMDISStrucFuncBase.h | 1 + 3 files changed, 26 insertions(+), 6 deletions(-) diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx index f27baf8e57..4938225996 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx @@ -115,8 +115,8 @@ void BYStrucFunc2021::KFactors(const Interaction * interaction, double myQ2 = this->Q2(interaction); double GD = 1. / TMath::Power(1.+myQ2/0.71, 2); // p elastic form factor double GD2 = TMath::Power(GD,2); - // double Ev = init_state.ProbeE(kRfHitNucRest); - double q0 = 0; + + double q0 = this->q0(interaction); double q02 = TMath::Power(q0,2); double KLW = q02 / ( q02 + fCvLW ); diff --git a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx index 33890e6e18..f54a0ee168 100644 --- a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx +++ b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx @@ -423,6 +423,24 @@ double QPMDISStrucFuncBase::Q2(const Interaction * interaction) const return 0; } //____________________________________________________________________________ +double QPMDISStrucFuncBase::q0(const Interaction * interaction) const +{ + const Kinematics & kinematics = interaction->Kine(); + + // Compute from Q2 and x + if (kinematics.KVSet(kKVQ2) || kinematics.KVSet(kKVq2)) { + const Kinematics & kinematics = interaction->Kine(); + const InitialState & init_state = interaction->InitState(); + double Mn = init_state.Tgt().HitNucP4Ptr()->M(); // could be off-shell + double Q2val = kinematics.Q2(); + double x = kinematics.x(); + double q0 = Q2val / ( 2 * Mn * x ) ; + return q0; + } + LOG("DISSF", pERROR) << "Could not compute Q2!"; + return 0; +} +//____________________________________________________________________________ double QPMDISStrucFuncBase::ScalingVar(const Interaction* interaction) const { // The scaling variable is set to the normal Bjorken x. @@ -432,7 +450,7 @@ double QPMDISStrucFuncBase::ScalingVar(const Interaction* interaction) const } //____________________________________________________________________________ void QPMDISStrucFuncBase::KFactors(const Interaction *, - double & kuv, double & kdv, double & kus, double & kds, double & kds ) const + double & kuv, double & kdv, double & kus, double & kds, double & kss ) const { // This is an abstract class: no model-specific correction // The PDF scaling variables are set to 1 @@ -551,8 +569,9 @@ void QPMDISStrucFuncBase::CalcPDFs(const Interaction * interaction) const double kval_d = 1.; double ksea_u = 1.; double ksea_d = 1.; - - this->KFactors(interaction, kval_u, kval_d, ksea_u, ksea_d); + double ksea_s = 1.; + + this->KFactors(interaction, kval_u, kval_d, ksea_u, ksea_d, ksea_s); #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ LOG("DISSF", pDEBUG) << "K-Factors:"; @@ -571,7 +590,7 @@ void QPMDISStrucFuncBase::CalcPDFs(const Interaction * interaction) const fPDF->ScaleDownValence (kval_d); fPDF->ScaleUpSea (ksea_u); fPDF->ScaleDownSea (ksea_d); - fPDF->ScaleStrange (ksea_d); + fPDF->ScaleStrange (ksea_s); fPDF->ScaleCharm (ksea_u); if(above_charm) { fPDFc->ScaleUpValence (kval_u); diff --git a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h index 56c2d3c485..a88cbd06a1 100644 --- a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h +++ b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h @@ -63,6 +63,7 @@ class QPMDISStrucFuncBase : public DISStructureFuncModelI { virtual void LoadConfig (void); virtual void InitPDF (void); virtual double Q2 (const Interaction * i) const; + virtual double q0 (const Interaction * i) const; virtual double ScalingVar (const Interaction * i) const; virtual void CalcPDFs (const Interaction * i) const; virtual double NuclMod (const Interaction * i) const; From 678c11e29ea8ee3988e15aa7167451346131a3c1 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Wed, 22 Oct 2025 10:46:21 +0200 Subject: [PATCH 10/68] adding Higher Order QCD modifications for F3 --- config/BYStrucFunc2021.xml | 10 +++++++++- config/QPMDISPXSec.xml | 4 +++- .../DeepInelastic/XSection/BYStrucFunc2021.cxx | 12 ++++++++++++ .../DeepInelastic/XSection/BYStrucFunc2021.h | 5 +++++ .../DeepInelastic/XSection/QPMDISStrucFuncBase.cxx | 14 +++++++++++++- .../DeepInelastic/XSection/QPMDISStrucFuncBase.h | 2 ++ 6 files changed, 44 insertions(+), 3 deletions(-) diff --git a/config/BYStrucFunc2021.xml b/config/BYStrucFunc2021.xml index 92af26824d..520d614741 100644 --- a/config/BYStrucFunc2021.xml +++ b/config/BYStrucFunc2021.xml @@ -15,6 +15,10 @@ BY-Cv1U double No Bodek-Yang u val K factor param BY-Cv2U double No Bodek-Yang u val K factor param BY-Cv1D double No Bodek-Yang d val K factor param BY-Cv2D double No Bodek-Yang d val K factor param +BY-H0 double No Bodek-Yang high order QCD parameter +BY-H1 double No Bodek-Yang high order QCD parameter +BY-H2 double No Bodek-Yang high order QCD parameter +BY-H3 double No Bodek-Yang high order QCD parameter +++ This model inherits from QPMDISStrucFunc and calls its configure +++ The following parameters have to be the same as in QPMDISStrucFunc.xml @@ -51,7 +55,10 @@ WeinbergAngle double No 0.341 0.323 0.561 - + 0.914 + 0.296 + -0.374 + 0.165 true + true true false 0.8 diff --git a/config/QPMDISPXSec.xml b/config/QPMDISPXSec.xml index e72ff32e52..32b873f0a0 100644 --- a/config/QPMDISPXSec.xml +++ b/config/QPMDISPXSec.xml @@ -27,6 +27,8 @@ WeinbergAngle double No - genie::BYStrucFunc/GRV982010 + WeakInt + 1. + genie::BYStrucFunc2021/Default diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx index 4938225996..1a244404d0 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx @@ -75,6 +75,10 @@ void BYStrucFunc2021::ReadBYParams(void) GetParam( "BY-Cv1D", fCv1D ) ; GetParam( "BY-Cv2D", fCv2D ) ; GetParam( "BY-CsS" , fCsS ) ; + GetParam( "BY-H0" , fH0 ) ; + GetParam( "BY-H1" , fH1 ) ; + GetParam( "BY-H2" , fH2 ) ; + GetParam( "BY-H3" , fH3 ) ; } //____________________________________________________________________________ void BYStrucFunc2021::Init(void) @@ -127,3 +131,11 @@ void BYStrucFunc2021::KFactors(const Interaction * interaction, kss = myQ2/(myQ2+fCsS); // K - s(sea) } //____________________________________________________________________________ +double BYStrucFunc2021::H(const Interaction * interaction) const { + // Overrides QPMDISStrucFuncBase::H() function to compute the correction of the BY 2021 update + // The correction is given by Eq. 34 of https://arxiv.org/pdf/2108.09240 + const Kinematics & kinematics = interaction->Kine(); + double bjx = kinematics.x(); + return fH0 + fH1 * bjx + fH2 * pow(bjx,2) + fH3 * pow(bjx,3); +} + diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h index 1bcf158297..038f83dc96 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h @@ -51,6 +51,7 @@ class BYStrucFunc2021 : public QPMDISStrucFuncBase { double ScalingVar (const Interaction * i) const; void KFactors (const Interaction * i, double & kuv, double & kdv, double & kus, double & kds, double & kss ) const; + double H(const Interaction * interaction) const ; // overrides QPMDISStrucFuncBase implementation // Bodek-Yang model-specific parameters @@ -64,6 +65,10 @@ class BYStrucFunc2021 : public QPMDISStrucFuncBase { double fCv2U; ///< U-val K factor parameter double fCv1D; ///< D-val K factor parameter double fCv2D; ///< D-val K factor parameter + double fH0; + double fH1; + double fH2; + double fH3; }; } // genie namespace diff --git a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx index f54a0ee168..98016ea9a0 100644 --- a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx +++ b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx @@ -107,6 +107,9 @@ void QPMDISStrucFuncBase::LoadConfig(void) //-- include R (~FL)? GetParam( "IncludeR", fIncludeR ) ; + //-- include H? + GetParam( "IncludeH", fIncludeH, false ) ; + //-- include nuclear factor (shadowing / anti-shadowing / ...) GetParam( "IncludeNuclMod", fIncludeNuclMod ) ; @@ -350,10 +353,12 @@ void QPMDISStrucFuncBase::Calculate(const Interaction * interaction) const double x = this->ScalingVar(interaction); double f = this->NuclMod (interaction); // nuclear modification double r = this->R (interaction); // R ~ FL - + double H = fIncludeH ? this->H(interaction) : 1; + #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ LOG("DISSF", pDEBUG) << "Nucl. mod = " << f; LOG("DISSF", pDEBUG) << "R(=FL/2xF1) = " << r; + LOG("DISSF", pDEBUG) << "H = " << H; #endif if(fUse2016Corrections) { @@ -389,6 +394,9 @@ void QPMDISStrucFuncBase::Calculate(const Interaction * interaction) const fF4 = 0.; // Nucl.Phys.B 84, 467 (1975) } + // Apply H(x,Q2) factor for xF3: + fF3 *= H; + #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ LOG("DISSF", pDEBUG) << "F1-F5 = " @@ -511,6 +519,10 @@ double QPMDISStrucFuncBase::R(const Interaction * interaction) const return 0; } //____________________________________________________________________________ +double QPMDISStrucFuncBase::H(const Interaction * interaction) const { + return 1 ; +} +//____________________________________________________________________________ void QPMDISStrucFuncBase::CalcPDFs(const Interaction * interaction) const { // Clean-up previous calculation diff --git a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h index a88cbd06a1..9d8f0402d0 100644 --- a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h +++ b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h @@ -68,6 +68,7 @@ class QPMDISStrucFuncBase : public DISStructureFuncModelI { virtual void CalcPDFs (const Interaction * i) const; virtual double NuclMod (const Interaction * i) const; virtual double R (const Interaction * i) const; + virtual double H (const Interaction * i) const; virtual void KFactors (const Interaction * i, double & kuv, double & kdv, double & kus, double & kds, double & kss ) const; // configuration @@ -75,6 +76,7 @@ class QPMDISStrucFuncBase : public DISStructureFuncModelI { double fQ2min; ///< min Q^2 allowed for PDFs: PDF(Q2 Date: Wed, 22 Oct 2025 11:00:50 +0200 Subject: [PATCH 11/68] adding modification of R correction factor at low Q2 --- .../DeepInelastic/XSection/BYStrucFunc2021.cxx | 14 ++++++++++++++ .../DeepInelastic/XSection/BYStrucFunc2021.h | 1 + 2 files changed, 15 insertions(+) diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx index 1a244404d0..5d1aaf509b 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx @@ -138,4 +138,18 @@ double BYStrucFunc2021::H(const Interaction * interaction) const { double bjx = kinematics.x(); return fH0 + fH1 * bjx + fH2 * pow(bjx,2) + fH3 * pow(bjx,3); } +//____________________________________________________________________________ +double BYStrucFunc2021::R(const Interaction * interaction) const { + + // Evaluate correction for Q2 < 0.3 GeV2/c4 according to Sec. 7 of https://arxiv.org/pdf/2108.09240 + double Q2 = this->Q2(interaction); + double R1998 = QPMDISStrucFuncBase::R(interaction); + if( Q2 < 0.3 ) { + Q2 = 0.3 ; + double Q4 = pow(Q2,2); + R1998 *= 3.633 * Q2 / ( Q4 + 1 ) ; + } + + return R1998 ; +} diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h index 038f83dc96..69ad76bab7 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h @@ -51,6 +51,7 @@ class BYStrucFunc2021 : public QPMDISStrucFuncBase { double ScalingVar (const Interaction * i) const; void KFactors (const Interaction * i, double & kuv, double & kdv, double & kus, double & kds, double & kss ) const; + double R(const Interaction * interaction) const ; // overrides QPMDISStrucFuncBase implementation double H(const Interaction * interaction) const ; // overrides QPMDISStrucFuncBase implementation // Bodek-Yang model-specific parameters From fdb33964e0eb2887a2616bc9561fc7aa91dd6421 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Wed, 22 Oct 2025 11:29:01 +0200 Subject: [PATCH 12/68] adding updates on nuclear corrections --- .../XSection/BYStrucFunc2021.cxx | 51 +++++++++++++++++++ .../DeepInelastic/XSection/BYStrucFunc2021.h | 2 +- 2 files changed, 52 insertions(+), 1 deletion(-) diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx index 5d1aaf509b..71243b5bdd 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx @@ -153,3 +153,54 @@ double BYStrucFunc2021::R(const Interaction * interaction) const { return R1998 ; } +//____________________________________________________________________________ +double BYStrucFunc2021::NuclMod(const Interaction * interaction) const { +// Nuclear modification to Fi +// The scaling variable can be overwritten to include corrections + + if( interaction->TestBit(kIAssumeFreeNucleon) ) return 1.0; + if( interaction->TestBit(kINoNuclearCorrection) ) return 1.0; + + double f = 1.; + if(fIncludeNuclMod) { + const Target & tgt = interaction->InitState().Tgt(); + const Kinematics & kinematics = interaction->Kine(); + double x = kinematics.x(); + int A = tgt.A(); + + double xv = TMath::Min(0.75, TMath::Max(0.05, x)); + double xv2 = xv * xv; + double xv3 = xv2 * xv; + double xv4 = xv3 * xv; + double xv5 = xv4 * xv; + + double f = 1.; + + // first factor goes from free nucleons to deuterium + if(A >= 2) { + f= 0.985*(1.+0.422*xv - 2.745*xv2 + 7.570*xv3 - 10.335*xv4 + 5.422*xv5); + } + + double chiTM = this->ScalingVar(interaction); + // 2nd factor goes from deuterium to iso-scalar iron + if(A > 2) { + f *= (1.096 - 0.38*chiTM - 0.3 * TMath::Exp(-23*chiTM) + 8 * pow(chiTM,15) ) ; + } + + if( A == 79 || A == 82 ) { // Gold nd Lead F2(Au,Pb)/F2(Fe) + f *= (0.932 + 2.461 * chiTM - 24.23*pow(chiTM,2) + 101.03*pow(chiTM,3) - 203.47*pow(chiTM,4) + 193.85*pow(chiTM,5) - 69.82*pow(chiTM,6)); + } + + if( A == 12 ) { // F2(Fe)/F2(C) + f /= ( 0.919 + 1.844*chiTM - 12.73*pow(chiTM,2) + 36.89*pow(chiTM,3) - 46.77*pow(chiTM,4) + 21.22*pow(chiTM,5)); + } + + return f; + +#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ + LOG("DISSF", pDEBUG) << "Nuclear factor for x of " << x << " = " << f; +#endif + } + + return f; +} diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h index 69ad76bab7..a4a611b079 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h @@ -53,7 +53,7 @@ class BYStrucFunc2021 : public QPMDISStrucFuncBase { double & kdv, double & kus, double & kds, double & kss ) const; double R(const Interaction * interaction) const ; // overrides QPMDISStrucFuncBase implementation double H(const Interaction * interaction) const ; // overrides QPMDISStrucFuncBase implementation - + double NuclMod(const Interaction * interaction) const; // Bodek-Yang model-specific parameters double fA; ///< better scaling var parameter A From d92bf6e71b02195e32f289599ac4ec067dee0808 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Wed, 22 Oct 2025 11:44:08 +0200 Subject: [PATCH 13/68] adding K factors parameterization for axial components --- config/BYStrucFunc2021.xml | 6 ++++++ .../DeepInelastic/XSection/BYStrucFunc2021.cxx | 17 +++++++++++++++++ .../DeepInelastic/XSection/BYStrucFunc2021.h | 7 +++++-- 3 files changed, 28 insertions(+), 2 deletions(-) diff --git a/config/BYStrucFunc2021.xml b/config/BYStrucFunc2021.xml index 520d614741..39c17ad953 100644 --- a/config/BYStrucFunc2021.xml +++ b/config/BYStrucFunc2021.xml @@ -15,6 +15,9 @@ BY-Cv1U double No Bodek-Yang u val K factor param BY-Cv2U double No Bodek-Yang u val K factor param BY-Cv1D double No Bodek-Yang d val K factor param BY-Cv2D double No Bodek-Yang d val K factor param +BY-PsA double No Bodek-Yang P sea axial param +BY-PvA double No Bodek-Yang P sea valance param +BY-PsC double No Bodek-Yang C sea axial param BY-H0 double No Bodek-Yang high order QCD parameter BY-H1 double No Bodek-Yang high order QCD parameter BY-H2 double No Bodek-Yang high order QCD parameter @@ -55,6 +58,9 @@ WeinbergAngle double No 0.341 0.323 0.561 + 0.55 + 0.30 + 0.75 0.914 0.296 -0.374 diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx index 71243b5bdd..81ca7d4858 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx @@ -75,6 +75,9 @@ void BYStrucFunc2021::ReadBYParams(void) GetParam( "BY-Cv1D", fCv1D ) ; GetParam( "BY-Cv2D", fCv2D ) ; GetParam( "BY-CsS" , fCsS ) ; + GetParam( "BY-PsA" , fPsA ) ; + GetParam( "BY-PvA" , fPvA ) ; + GetParam( "BY-CsA" , fCsA ) ; GetParam( "BY-H0" , fH0 ) ; GetParam( "BY-H1" , fH1 ) ; GetParam( "BY-H2" , fH2 ) ; @@ -93,6 +96,13 @@ void BYStrucFunc2021::Init(void) fCv1D = 0; fCv2D = 0; fCsS = 0; + fPsA = 0; + fPvA = 0; + fCsA = 0; + fH0 = 0; + fH1 = 0; + fH2 = 0; + fH3 = 0; } //____________________________________________________________________________ double BYStrucFunc2021::ScalingVar(const Interaction * interaction) const @@ -131,6 +141,13 @@ void BYStrucFunc2021::KFactors(const Interaction * interaction, kss = myQ2/(myQ2+fCsS); // K - s(sea) } //____________________________________________________________________________ +void BYStrucFunc2021::KAxialFactors (const Interaction * i, double & ksea, double & kvalance ) const { + // https://arxiv.org/pdf/2108.09240 Sec 11.2 + double myQ2 = this->Q2(interaction); + ksea = ( myQ2 + kPsA*kCsA ) / ( myQ2 + kCsA ) ; + kvalance = ( myQ2 + kPvA * 0.18 ) / ( myQ2 + 0.18 ) ; +} +//____________________________________________________________________________ double BYStrucFunc2021::H(const Interaction * interaction) const { // Overrides QPMDISStrucFuncBase::H() function to compute the correction of the BY 2021 update // The correction is given by Eq. 34 of https://arxiv.org/pdf/2108.09240 diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h index a4a611b079..14ec11b8e3 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h @@ -49,8 +49,8 @@ class BYStrucFunc2021 : public QPMDISStrucFuncBase { // override part of the DISStructureFuncModel implementation // to compute all the corrections applied by the Bodek-Yang model. double ScalingVar (const Interaction * i) const; - void KFactors (const Interaction * i, double & kuv, - double & kdv, double & kus, double & kds, double & kss ) const; + void KFactors (const Interaction * i, double & kuv, double & kdv, double & kus, double & kds, double & kss ) const; + void KAxialFactors (const Interaction * i, double & ksea, double & kvalance ) const; double R(const Interaction * interaction) const ; // overrides QPMDISStrucFuncBase implementation double H(const Interaction * interaction) const ; // overrides QPMDISStrucFuncBase implementation double NuclMod(const Interaction * interaction) const; @@ -66,6 +66,9 @@ class BYStrucFunc2021 : public QPMDISStrucFuncBase { double fCv2U; ///< U-val K factor parameter double fCv1D; ///< D-val K factor parameter double fCv2D; ///< D-val K factor parameter + double fPsA; ///< P-axial sea parameter + double fPvA; ///< P-axial valance paramter + double fCsA; ///< C-axial sea parameter double fH0; double fH1; double fH2; From 9f4ff2cd92f2faeb87ca3a09c3af1a0d19ed21e8 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Wed, 22 Oct 2025 13:11:22 +0200 Subject: [PATCH 14/68] adding KAxial factors --- config/BYStrucFunc2021.xml | 4 +++ .../XSection/BYStrucFunc2021.cxx | 28 +++++++++++++++++-- .../DeepInelastic/XSection/BYStrucFunc2021.h | 2 ++ 3 files changed, 31 insertions(+), 3 deletions(-) diff --git a/config/BYStrucFunc2021.xml b/config/BYStrucFunc2021.xml index 39c17ad953..ed2b664931 100644 --- a/config/BYStrucFunc2021.xml +++ b/config/BYStrucFunc2021.xml @@ -18,6 +18,8 @@ BY-Cv2D double No Bodek-Yang d val K factor param BY-PsA double No Bodek-Yang P sea axial param BY-PvA double No Bodek-Yang P sea valance param BY-PsC double No Bodek-Yang C sea axial param +BY-CaLW-nu double No Low-q0 axial neutrino parameter +BY-CaLW-nubar double No Low-q0 axial anti-neutrino parameter BY-H0 double No Bodek-Yang high order QCD parameter BY-H1 double No Bodek-Yang high order QCD parameter BY-H2 double No Bodek-Yang high order QCD parameter @@ -61,6 +63,8 @@ WeinbergAngle double No 0.55 0.30 0.75 + 0.436 + 0.654 0.914 0.296 -0.374 diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx index 81ca7d4858..4abdfbe1fa 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx @@ -19,6 +19,7 @@ #include "Physics/DeepInelastic/XSection/BYStrucFunc2021.h" #include "Framework/Conventions/Constants.h" #include "Framework/Messenger/Messenger.h" +#include "Framework/ParticleData/PDGUtils.h" using namespace genie; using namespace genie::constants; @@ -78,6 +79,8 @@ void BYStrucFunc2021::ReadBYParams(void) GetParam( "BY-PsA" , fPsA ) ; GetParam( "BY-PvA" , fPvA ) ; GetParam( "BY-CsA" , fCsA ) ; + GetParam( "BY-CaLW-nubar" , fCaLW_nu ) ; + GetParam( "BY-CaLW-nubar" , fCaLW_nubar ) ; GetParam( "BY-H0" , fH0 ) ; GetParam( "BY-H1" , fH1 ) ; GetParam( "BY-H2" , fH2 ) ; @@ -99,6 +102,8 @@ void BYStrucFunc2021::Init(void) fPsA = 0; fPvA = 0; fCsA = 0; + fCaLW_nu = 0; + fCaLW_nubar = 0; fH0 = 0; fH1 = 0; fH2 = 0; @@ -141,11 +146,28 @@ void BYStrucFunc2021::KFactors(const Interaction * interaction, kss = myQ2/(myQ2+fCsS); // K - s(sea) } //____________________________________________________________________________ -void BYStrucFunc2021::KAxialFactors (const Interaction * i, double & ksea, double & kvalance ) const { +void BYStrucFunc2021::KAxialFactors (const Interaction * interaction, double & ksea, double & kvalance ) const { // https://arxiv.org/pdf/2108.09240 Sec 11.2 double myQ2 = this->Q2(interaction); - ksea = ( myQ2 + kPsA*kCsA ) / ( myQ2 + kCsA ) ; - kvalance = ( myQ2 + kPvA * 0.18 ) / ( myQ2 + 0.18 ) ; + ksea = ( myQ2 + fPsA*fCsA ) / ( myQ2 + fCsA ) ; + kvalance = ( myQ2 + fPvA * 0.18 ) / ( myQ2 + 0.18 ) ; + + // Compute Low-q0 modificaion factor for neutrinos and anti-neutrinos + double q0 = this->q0(interaction); + double q02 = TMath::Power(q0,2); + const ProcessInfo & proc_info = interaction->ProcInfo(); + const InitialState & init_state = interaction->InitState(); + int probe_pdgc = init_state.ProbePdg(); + bool is_nu = pdg::IsNeutrino ( probe_pdgc ); + bool is_nubar = pdg::IsAntiNeutrino ( probe_pdgc ); + double KLW = 1; + if( is_nu ) { + KLW = ( q02 + fCaLW_nu ) / q02 ; + } else if( is_nubar ) { + KLW = ( q02 + fCaLW_nubar ) / q02 ; + } + // double check it only affects valance factors: + kvalance *= KLW ; } //____________________________________________________________________________ double BYStrucFunc2021::H(const Interaction * interaction) const { diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h index 14ec11b8e3..7ed42b640a 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h @@ -69,6 +69,8 @@ class BYStrucFunc2021 : public QPMDISStrucFuncBase { double fPsA; ///< P-axial sea parameter double fPvA; ///< P-axial valance paramter double fCsA; ///< C-axial sea parameter + double fCaLW_nubar; ///< C-axial neutrino LW parameter + double fCaLW_nu; ///< C-axial anti-neutrino LW parameter double fH0; double fH1; double fH2; From 3fa55305d3ef9225aa5a369a22aa3396daed069f Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Wed, 22 Oct 2025 14:07:36 +0200 Subject: [PATCH 15/68] fixing chi calculation for charm production --- .../XSection/DMBYStrucFunc.cxx | 2 +- .../XSection/DMBYStrucFunc.h | 2 +- .../DeepInelastic/XSection/BYStrucFunc.cxx | 2 +- .../DeepInelastic/XSection/BYStrucFunc.h | 2 +- .../XSection/BYStrucFunc2021.cxx | 5 +- .../DeepInelastic/XSection/BYStrucFunc2021.h | 2 +- .../XSection/QPMDISStrucFuncBase.cxx | 268 +++++++++--------- .../XSection/QPMDISStrucFuncBase.h | 2 +- 8 files changed, 146 insertions(+), 139 deletions(-) diff --git a/src/Physics/BoostedDarkMatter/XSection/DMBYStrucFunc.cxx b/src/Physics/BoostedDarkMatter/XSection/DMBYStrucFunc.cxx index 5d04f67326..bda8701828 100644 --- a/src/Physics/BoostedDarkMatter/XSection/DMBYStrucFunc.cxx +++ b/src/Physics/BoostedDarkMatter/XSection/DMBYStrucFunc.cxx @@ -92,7 +92,7 @@ void DMBYStrucFunc::Init(void) fCv2D = 0; } //____________________________________________________________________________ -double DMBYStrucFunc::ScalingVar(const Interaction * interaction) const +double DMBYStrucFunc::ScalingVar(const Interaction * interaction, double Mf ) const { // Overrides QPMDMDISStrucFuncBase::ScalingVar() to compute the BY scaling var diff --git a/src/Physics/BoostedDarkMatter/XSection/DMBYStrucFunc.h b/src/Physics/BoostedDarkMatter/XSection/DMBYStrucFunc.h index 2b1eee9f41..d3e7c02da4 100644 --- a/src/Physics/BoostedDarkMatter/XSection/DMBYStrucFunc.h +++ b/src/Physics/BoostedDarkMatter/XSection/DMBYStrucFunc.h @@ -43,7 +43,7 @@ class DMBYStrucFunc : public QPMDMDISStrucFuncBase { // override part of the DISStructureFuncModel implementation // to compute all the corrections applied by the Bodek-Yang model. - double ScalingVar (const Interaction * i) const; + double ScalingVar (const Interaction * i, double Mf = 0) const; void KFactors (const Interaction * i, double & kuv, double & kdv, double & kus, double & kds) const; diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc.cxx index 310538cf9b..e5f370da92 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc.cxx @@ -84,7 +84,7 @@ void BYStrucFunc::Init(void) fCv2D = 0; } //____________________________________________________________________________ -double BYStrucFunc::ScalingVar(const Interaction * interaction) const +double BYStrucFunc::ScalingVar(const Interaction * interaction, double Mf ) const { // Overrides QPMDISStrucFuncBase::ScalingVar() to compute the BY scaling var diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc.h b/src/Physics/DeepInelastic/XSection/BYStrucFunc.h index 80ccd5cf9a..8e00d64d22 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc.h +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc.h @@ -45,7 +45,7 @@ class BYStrucFunc : public QPMDISStrucFuncBase { // override part of the DISStructureFuncModel implementation // to compute all the corrections applied by the Bodek-Yang model. - double ScalingVar (const Interaction * i) const; + double ScalingVar (const Interaction * i, double Mf = 0 ) const; void KFactors (const Interaction * i, double & kuv, double & kdv, double & kus, double & kds, double & kss ) const; diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx index 4abdfbe1fa..0874a865b8 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx @@ -110,7 +110,7 @@ void BYStrucFunc2021::Init(void) fH3 = 0; } //____________________________________________________________________________ -double BYStrucFunc2021::ScalingVar(const Interaction * interaction) const +double BYStrucFunc2021::ScalingVar(const Interaction * interaction, double Mf ) const { // Overrides QPMDISStrucFuncBase::ScalingVar() to compute the BY scaling var @@ -121,7 +121,8 @@ double BYStrucFunc2021::ScalingVar(const Interaction * interaction) const LOG("BodekYang", pDEBUG) << "Q2 at scaling var calculation = " << myQ2; double a = TMath::Power( 2*kProtonMass*x, 2 ) / myQ2; - double xw = 2*x*(myQ2+fB) / (myQ2*(1.+TMath::Sqrt(1+a)) + 2*fA*x); + double Mf2 = TMath::Power( Mf, 2 ) ; + double xw = 2*x*(myQ2+Mf2+fB) / (myQ2*(1.+TMath::Sqrt(1+a)) + 2*fA*x); return xw; } //____________________________________________________________________________ diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h index 7ed42b640a..269445da58 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h @@ -48,7 +48,7 @@ class BYStrucFunc2021 : public QPMDISStrucFuncBase { // override part of the DISStructureFuncModel implementation // to compute all the corrections applied by the Bodek-Yang model. - double ScalingVar (const Interaction * i) const; + double ScalingVar (const Interaction * i, double Mf = 0 ) const; void KFactors (const Interaction * i, double & kuv, double & kdv, double & kus, double & kds, double & kss ) const; void KAxialFactors (const Interaction * i, double & ksea, double & kvalance ) const; double R(const Interaction * interaction) const ; // overrides QPMDISStrucFuncBase implementation diff --git a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx index 98016ea9a0..24569b98b7 100644 --- a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx +++ b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx @@ -1,25 +1,25 @@ //____________________________________________________________________________ /* - Copyright (c) 2003-2025, The GENIE Collaboration - For the full text of the license visit http://copyright.genie-mc.org + Copyright (c) 2003-2025, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org - Costas Andreopoulos - University of Liverpool + Costas Andreopoulos + University of Liverpool - This GENIE code was adapted from the neugen3 code co-authored by Donna Naples - (Pittsburgh U.), Hugh Gallagher (Tufts U), and Costas Andreopoulos (RAL) + This GENIE code was adapted from the neugen3 code co-authored by Donna Naples + (Pittsburgh U.), Hugh Gallagher (Tufts U), and Costas Andreopoulos (RAL) - A fix was installed (Aug 12, 2014) by Brian Tice (Rochester) so that - the nuclear modification to the pdf should be calculated in terms - of the experimental x, not the rescaled x. The same goes for R(x,Q2). + A fix was installed (Aug 12, 2014) by Brian Tice (Rochester) so that + the nuclear modification to the pdf should be calculated in terms + of the experimental x, not the rescaled x. The same goes for R(x,Q2). - A fix of the scaling variable used for the relations between structure - functions was installed by C. Bronner and J. Morrison Jun 06, 2016 - after it was confirmed by A. Bodek that x and not the modified - scaling variable should be used there. + A fix of the scaling variable used for the relations between structure + functions was installed by C. Bronner and J. Morrison Jun 06, 2016 + after it was confirmed by A. Bodek that x and not the modified + scaling variable should be used there. - Changes required to implement the GENIE Boosted Dark Matter module - were installed by Josh Berger (Univ. of Wisconsin) + Changes required to implement the GENIE Boosted Dark Matter module + were installed by Josh Berger (Univ. of Wisconsin) */ //____________________________________________________________________________ @@ -42,19 +42,19 @@ using namespace genie::constants; //____________________________________________________________________________ QPMDISStrucFuncBase::QPMDISStrucFuncBase() : -DISStructureFuncModelI() + DISStructureFuncModelI() { this->InitPDF(); } //____________________________________________________________________________ QPMDISStrucFuncBase::QPMDISStrucFuncBase(string name) : -DISStructureFuncModelI(name) + DISStructureFuncModelI(name) { this->InitPDF(); } //____________________________________________________________________________ QPMDISStrucFuncBase::QPMDISStrucFuncBase(string name, string config): -DISStructureFuncModelI(name, config) + DISStructureFuncModelI(name, config) { this->InitPDF(); } @@ -83,7 +83,7 @@ void QPMDISStrucFuncBase::LoadConfig(void) //-- pdf const PDFModelI * pdf_model = - dynamic_cast (this->SubAlg("PDF-Set")); + dynamic_cast (this->SubAlg("PDF-Set")); fPDF -> SetModel(pdf_model); fPDFc -> SetModel(pdf_model); @@ -132,7 +132,7 @@ void QPMDISStrucFuncBase::LoadConfig(void) //____________________________________________________________________________ void QPMDISStrucFuncBase::InitPDF(void) { - // evaluated at: + // evaluated at: fPDF = new PDF(); // x = computed (+/-corrections) scaling var, Q2 fPDFc = new PDF(); // x = computed charm slow re-scaling var, Q2 } @@ -186,42 +186,42 @@ void QPMDISStrucFuncBase::Calculate(const Interaction * interaction) const if(tgt.HitQrkIsSet()) { - switch_uv = 0.; - switch_us = 0.; - switch_ubar = 0.; - switch_dv = 0.; - switch_ds = 0.; - switch_dbar = 0.; - switch_s = 0.; - switch_sbar = 0.; - switch_c = 0.; - switch_cbar = 0.; - - int qpdg = tgt.HitQrkPdg(); - bool sea = tgt.HitSeaQrk(); - - bool is_u = pdg::IsUQuark (qpdg); - bool is_ubar = pdg::IsAntiUQuark (qpdg); - bool is_d = pdg::IsDQuark (qpdg); - bool is_dbar = pdg::IsAntiDQuark (qpdg); - bool is_s = pdg::IsSQuark (qpdg); - bool is_sbar = pdg::IsAntiSQuark (qpdg); - bool is_c = pdg::IsCQuark (qpdg); - bool is_cbar = pdg::IsAntiCQuark (qpdg); - - if (!sea && is_u ) { switch_uv = 1; } - else if ( sea && is_u ) { switch_us = 1; } - else if ( sea && is_ubar) { switch_ubar = 1; } - else if (!sea && is_d ) { switch_dv = 1; } - else if ( sea && is_d ) { switch_ds = 1; } - else if ( sea && is_dbar) { switch_dbar = 1; } - else if ( sea && is_s ) { switch_s = 1; } - else if ( sea && is_sbar) { switch_sbar = 1; } - else if ( sea && is_c ) { switch_c = 1; } - else if ( sea && is_cbar) { switch_cbar = 1; } - else return; - - // make sure user inputs make sense + switch_uv = 0.; + switch_us = 0.; + switch_ubar = 0.; + switch_dv = 0.; + switch_ds = 0.; + switch_dbar = 0.; + switch_s = 0.; + switch_sbar = 0.; + switch_c = 0.; + switch_cbar = 0.; + + int qpdg = tgt.HitQrkPdg(); + bool sea = tgt.HitSeaQrk(); + + bool is_u = pdg::IsUQuark (qpdg); + bool is_ubar = pdg::IsAntiUQuark (qpdg); + bool is_d = pdg::IsDQuark (qpdg); + bool is_dbar = pdg::IsAntiDQuark (qpdg); + bool is_s = pdg::IsSQuark (qpdg); + bool is_sbar = pdg::IsAntiSQuark (qpdg); + bool is_c = pdg::IsCQuark (qpdg); + bool is_cbar = pdg::IsAntiCQuark (qpdg); + + if (!sea && is_u ) { switch_uv = 1; } + else if ( sea && is_u ) { switch_us = 1; } + else if ( sea && is_ubar) { switch_ubar = 1; } + else if (!sea && is_d ) { switch_dv = 1; } + else if ( sea && is_d ) { switch_ds = 1; } + else if ( sea && is_dbar) { switch_dbar = 1; } + else if ( sea && is_s ) { switch_s = 1; } + else if ( sea && is_sbar) { switch_sbar = 1; } + else if ( sea && is_c ) { switch_c = 1; } + else if ( sea && is_cbar) { switch_cbar = 1; } + else return; + + // make sure user inputs make sense if(is_nu && is_CC && is_u ) return; if(is_nu && is_CC && is_c ) return; if(is_nu && is_CC && is_dbar) return; @@ -271,14 +271,14 @@ void QPMDISStrucFuncBase::Calculate(const Interaction * interaction) const double gad2 = TMath::Power(gad, 2.); double q2 = (switch_uv * fuv + switch_us * fus + switch_c * fc) * (gvu2+gau2) + - (switch_dv * fdv + switch_ds * fds + switch_s * fs) * (gvd2+gad2); + (switch_dv * fdv + switch_ds * fds + switch_s * fs) * (gvd2+gad2); double q3 = (switch_uv * fuv + switch_us * fus + switch_c * fc) * (2*gvu*gau) + - (switch_dv * fdv + switch_ds * fds + switch_s * fs) * (2*gvd*gad); + (switch_dv * fdv + switch_ds * fds + switch_s * fs) * (2*gvd*gad); double qb2 = (switch_ubar * fus + switch_cbar * fc) * (gvu2+gau2) + - (switch_dbar * fds + switch_sbar * fs) * (gvd2+gad2); + (switch_dbar * fds + switch_sbar * fs) * (gvd2+gad2); double qb3 = (switch_ubar * fus + switch_cbar * fc) * (2*gvu*gau) + - (switch_dbar * fds + switch_sbar * fs) * (2*gvd*gad); + (switch_dbar * fds + switch_sbar * fs) * (2*gvd*gad); #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ LOG("DISSF", pINFO) << "f2 : q = " << q2 << ", bar{q} = " << qb2; @@ -296,30 +296,30 @@ void QPMDISStrucFuncBase::Calculate(const Interaction * interaction) const if (is_nu) { q = ( switch_dv * fdv + switch_ds * fds ) * fVud2 + - ( switch_s * fs ) * fVus2 + - ( switch_dv * fdv_c + switch_ds * fds_c ) * fVcd2 + - ( switch_s * fs_c ) * fVcs2; + ( switch_s * fs ) * fVus2 + + ( switch_dv * fdv_c + switch_ds * fds_c ) * fVcd2 + + ( switch_s * fs_c ) * fVcs2; qbar = ( switch_ubar * fus ) * fVud2 + - ( switch_ubar * fus ) * fVus2 + - ( switch_cbar * fc_c ) * fVcd2 + - ( switch_cbar * fc_c ) * fVcs2; + ( switch_ubar * fus ) * fVus2 + + ( switch_cbar * fc_c ) * fVcd2 + + ( switch_cbar * fc_c ) * fVcs2; } else - if (is_nubar) { - q = ( switch_uv * fuv + switch_us * fus ) * fVud2 + - ( switch_uv * fuv + switch_us * fus ) * fVus2 + - ( switch_c * fc_c ) * fVcd2 + - ( switch_c * fc_c ) * fVcs2; - - qbar = ( switch_dbar * fds_c ) * fVcd2 + - ( switch_dbar * fds ) * fVud2 + - ( switch_sbar * fs ) * fVus2 + - ( switch_sbar * fs_c ) * fVcs2; - } - else { - return; - } + if (is_nubar) { + q = ( switch_uv * fuv + switch_us * fus ) * fVud2 + + ( switch_uv * fuv + switch_us * fus ) * fVus2 + + ( switch_c * fc_c ) * fVcd2 + + ( switch_c * fc_c ) * fVcs2; + + qbar = ( switch_dbar * fds_c ) * fVcd2 + + ( switch_dbar * fds ) * fVud2 + + ( switch_sbar * fs ) * fVus2 + + ( switch_sbar * fs_c ) * fVcs2; + } + else { + return; + } F2val = 2*(q+qbar); xF3val = 2*(q-qbar); @@ -399,15 +399,15 @@ void QPMDISStrucFuncBase::Calculate(const Interaction * interaction) const #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ LOG("DISSF", pDEBUG) - << "F1-F5 = " - << fF1 << ", " << fF2 << ", " << fF3 << ", " << fF4 << ", " << fF5; + << "F1-F5 = " + << fF1 << ", " << fF2 << ", " << fF3 << ", " << fF4 << ", " << fF5; #endif } //____________________________________________________________________________ double QPMDISStrucFuncBase::Q2(const Interaction * interaction) const { -// Return Q2 from the kinematics or, if not set, compute it from x,y -// The x might be corrected + // Return Q2 from the kinematics or, if not set, compute it from x,y + // The x might be corrected const Kinematics & kinematics = interaction->Kine(); @@ -449,20 +449,26 @@ double QPMDISStrucFuncBase::q0(const Interaction * interaction) const return 0; } //____________________________________________________________________________ -double QPMDISStrucFuncBase::ScalingVar(const Interaction* interaction) const +double QPMDISStrucFuncBase::ScalingVar(const Interaction* interaction, double Mf ) const { -// The scaling variable is set to the normal Bjorken x. -// Override DISStructureFuncModel::ScalingVar() to compute corrections - + // The scaling variable is set to the normal Bjorken x. + // Override DISStructureFuncModel::ScalingVar() to compute corrections + if( Mf != 0 ) { // For Charm production only + const Target & tgt = interaction->InitState().Tgt(); + double x = this->ScalingVar(interaction); + double Q2val = this->Q2(interaction); + double M = tgt.HitNucP4().M(); + return utils::kinematics::SlowRescalingVar(x, Q2val, M, fMc); + } return interaction->Kine().x(); } //____________________________________________________________________________ void QPMDISStrucFuncBase::KFactors(const Interaction *, - double & kuv, double & kdv, double & kus, double & kds, double & kss ) const + double & kuv, double & kdv, double & kus, double & kds, double & kss ) const { -// This is an abstract class: no model-specific correction -// The PDF scaling variables are set to 1 -// Override this method to compute model-dependent corrections + // This is an abstract class: no model-specific correction + // The PDF scaling variables are set to 1 + // Override this method to compute model-dependent corrections kuv = 1.; kdv = 1.; @@ -473,26 +479,26 @@ void QPMDISStrucFuncBase::KFactors(const Interaction *, //____________________________________________________________________________ double QPMDISStrucFuncBase::NuclMod(const Interaction * interaction) const { -// Nuclear modification to Fi -// The scaling variable can be overwritten to include corrections + // Nuclear modification to Fi + // The scaling variable can be overwritten to include corrections if( interaction->TestBit(kIAssumeFreeNucleon) ) return 1.0; if( interaction->TestBit(kINoNuclearCorrection) ) return 1.0; double f = 1.; if(fIncludeNuclMod) { - const Target & tgt = interaction->InitState().Tgt(); - -// The x used for computing the DIS Nuclear correction factor should be the -// experimental x, not the rescaled x or off-shell-rest-frame version of x -// (i.e. selected x). Since we do not have access to experimental x at this -// point in the calculation, just use selected x. - const Kinematics & kine = interaction->Kine(); - double x = kine.x(); - int A = tgt.A(); - f = utils::nuclear::DISNuclFactor(x,A); + const Target & tgt = interaction->InitState().Tgt(); + + // The x used for computing the DIS Nuclear correction factor should be the + // experimental x, not the rescaled x or off-shell-rest-frame version of x + // (i.e. selected x). Since we do not have access to experimental x at this + // point in the calculation, just use selected x. + const Kinematics & kine = interaction->Kine(); + double x = kine.x(); + int A = tgt.A(); + f = utils::nuclear::DISNuclFactor(x,A); #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ - LOG("DISSF", pDEBUG) << "Nuclear factor for x of " << x << " = " << f; + LOG("DISSF", pDEBUG) << "Nuclear factor for x of " << x << " = " << f; #endif } @@ -501,17 +507,17 @@ double QPMDISStrucFuncBase::NuclMod(const Interaction * interaction) const //____________________________________________________________________________ double QPMDISStrucFuncBase::R(const Interaction * interaction) const { -// Computes R ( ~ longitudinal structure function FL = R * 2xF1) -// The scaling variable can be overwritten to include corrections + // Computes R ( ~ longitudinal structure function FL = R * 2xF1) + // The scaling variable can be overwritten to include corrections -// The x used for computing the DIS Nuclear correction factor should be the -// experimental x, not the rescaled x or off-shell-rest-frame version of x -// (i.e. selected x). Since we do not have access to experimental x at this -// point in the calculation, just use selected x. + // The x used for computing the DIS Nuclear correction factor should be the + // experimental x, not the rescaled x or off-shell-rest-frame version of x + // (i.e. selected x). Since we do not have access to experimental x at this + // point in the calculation, just use selected x. if(fIncludeR) { const Kinematics & kine = interaction->Kine(); double x = kine.x(); -// double x = this->ScalingVar(interaction); + // double x = this->ScalingVar(interaction); double Q2val = this->Q2(interaction); double Rval = utils::phys::RWhitlow(x, Q2val); return Rval; @@ -548,32 +554,32 @@ void QPMDISStrucFuncBase::CalcPDFs(const Interaction * interaction) const // Check whether it is above charm threshold bool above_charm = - utils::kinematics::IsAboveCharmThreshold(x, Q2val, M, fMc); + utils::kinematics::IsAboveCharmThreshold(x, Q2val, M, fMc); if(above_charm) { #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ LOG("DISSF", pDEBUG) << "The event is above the charm threshold (mcharm = " << fMc << ")"; #endif if(fCharmOff) { - LOG("DISSF", pINFO) << "Charm production is turned off"; + LOG("DISSF", pINFO) << "Charm production is turned off"; } else { - // compute the slow rescaling var - double xc = utils::kinematics::SlowRescalingVar(x, Q2val, M, fMc); - if(xc<0 || xc>1) { - LOG("DISSF", pINFO) << "Unphys. slow rescaling var: xc = " << xc; - } else { - // compute PDFs at (xc,Q2) + // compute the slow rescaling var + double xc = ScalingVar(interaction, fMc); + if(xc<0 || xc>1) { + LOG("DISSF", pINFO) << "Unphys. slow rescaling var: xc = " << xc; + } else { + // compute PDFs at (xc,Q2) #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ - LOG("DISSF", pDEBUG) - << "Calculating PDFs @ xc (slow rescaling) = " << x << ", Q2 = " << Q2val; + LOG("DISSF", pDEBUG) + << "Calculating PDFs @ xc (slow rescaling) = " << x << ", Q2 = " << Q2val; #endif - fPDFc->Calculate(xc, Q2pdf); - } + fPDFc->Calculate(xc, Q2pdf); + } }// charm off? }//above charm thr? else { LOG("DISSF", pDEBUG) - << "The event is below the charm threshold (mcharm = " << fMc << ")"; + << "The event is below the charm threshold (mcharm = " << fMc << ")"; } // Compute the K factors @@ -605,12 +611,12 @@ void QPMDISStrucFuncBase::CalcPDFs(const Interaction * interaction) const fPDF->ScaleStrange (ksea_s); fPDF->ScaleCharm (ksea_u); if(above_charm) { - fPDFc->ScaleUpValence (kval_u); - fPDFc->ScaleDownValence (kval_d); - fPDFc->ScaleUpSea (ksea_u); - fPDFc->ScaleDownSea (ksea_d); - fPDFc->ScaleStrange (ksea_d); - fPDFc->ScaleCharm (ksea_u); + fPDFc->ScaleUpValence (kval_u); + fPDFc->ScaleDownValence (kval_d); + fPDFc->ScaleUpSea (ksea_u); + fPDFc->ScaleDownSea (ksea_d); + fPDFc->ScaleStrange (ksea_d); + fPDFc->ScaleCharm (ksea_u); } // Rules of thumb diff --git a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h index 9d8f0402d0..86d4138df7 100644 --- a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h +++ b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h @@ -64,7 +64,7 @@ class QPMDISStrucFuncBase : public DISStructureFuncModelI { virtual void InitPDF (void); virtual double Q2 (const Interaction * i) const; virtual double q0 (const Interaction * i) const; - virtual double ScalingVar (const Interaction * i) const; + virtual double ScalingVar (const Interaction * i, double Mf = 0 ) const; virtual void CalcPDFs (const Interaction * i) const; virtual double NuclMod (const Interaction * i) const; virtual double R (const Interaction * i) const; From d3a8fe4acd82192be6cfc0ad872b203f028f7d4b Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Wed, 22 Oct 2025 14:20:10 +0200 Subject: [PATCH 16/68] avoid double counting with chiTM cut --- src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx index 0874a865b8..8a8b5c4153 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx @@ -224,6 +224,8 @@ double BYStrucFunc2021::NuclMod(const Interaction * interaction) const { double chiTM = this->ScalingVar(interaction); // 2nd factor goes from deuterium to iso-scalar iron if(A > 2) { + // Fermi motion is already accounted for at high chiTM. To avoid double-counting, we limit the chiTM range: + if(! interaction->TestBit(kIAssumeFreeNucleon) && chiTM > 0.65 ) chiTM = 0.65; f *= (1.096 - 0.38*chiTM - 0.3 * TMath::Exp(-23*chiTM) + 8 * pow(chiTM,15) ) ; } From 02dada612bda24a403e810d6c03b4d9f58cb170c Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Wed, 22 Oct 2025 14:31:00 +0200 Subject: [PATCH 17/68] adding new configuration in KNOTuned model --- config/KNOTunedQPMDISPXSec.xml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/config/KNOTunedQPMDISPXSec.xml b/config/KNOTunedQPMDISPXSec.xml index 1abf6169e8..d51bfe1a45 100644 --- a/config/KNOTunedQPMDISPXSec.xml +++ b/config/KNOTunedQPMDISPXSec.xml @@ -24,4 +24,8 @@ Wcut double no Co 1. + + genie::QPMDISPXSec/GRV982010 + + From e8cd9cda6c5645caca200b9e1884dec27d83dac9 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Thu, 23 Oct 2025 11:03:34 +0200 Subject: [PATCH 18/68] adding changes --- config/BergerSehgalRESPXSec2014.xml | 1 - config/KNOTunedQPMDISPXSec.xml | 6 +++++- config/QPMDISPXSec.xml | 6 ++++-- 3 files changed, 9 insertions(+), 4 deletions(-) diff --git a/config/BergerSehgalRESPXSec2014.xml b/config/BergerSehgalRESPXSec2014.xml index f4d9b1b814..a0fd33270b 100644 --- a/config/BergerSehgalRESPXSec2014.xml +++ b/config/BergerSehgalRESPXSec2014.xml @@ -48,7 +48,6 @@ RES-Mb double Yes Mb parameter used for the GA calculati WeakInt,NonResBackground,CKM,FermiGas - @@ -80,7 +79,6 @@ WeinbergAngle double No true true true - false 0.8 diff --git a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx index 26fd60ee60..ed5b17ac14 100644 --- a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx +++ b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx @@ -114,7 +114,7 @@ void QPMDISStrucFuncBase::LoadConfig(void) GetParam( "IncludeNuclMod", fIncludeNuclMod ) ; //-- Use 2016 SF relation corrections - GetParam( "Use2016Corrections", fUse2016Corrections ) ; + GetParam( "Use2016Corrections", fUse2016Corrections, false ) ; //-- Set min for relation between 2xF1 and F2 GetParam( "LowQ2CutoffF1F2", fLowQ2CutoffF1F2 ) ; @@ -360,15 +360,37 @@ void QPMDISStrucFuncBase::Calculate(const Interaction * interaction) const LOG("DISSF", pDEBUG) << "R(=FL/2xF1) = " << r; LOG("DISSF", pDEBUG) << "H = " << H; #endif - - double a = TMath::Power(x,2.) / TMath::Max(Q2val, fLowQ2CutoffF1F2); - double c = (1. + 4. * kNucleonMass2 * a) / (1.+r); - - fF3 = f * H * xF3val / x; - fF2 = f * F2val; - fF1 = fF2 * 0.5 * c / x; - fF5 = fF2 / x; // Albright-Jarlskog relation - fF4 = 0.; // Nucl.Phys.B 84, 467 (1975) + + if(fUse2016Corrections) { + //It was confirmed by A.Bodek that the modified scaling variable + //should just be used to compute the strucure functions F2 and xF3, + //but that the usual Bjorken x should be used for the relations + //between the structure functions. + //For the same reason remove the freezing of Q2 at 0.8 for those relations, + //although it has not been explicitly asked to A.Bodek if it should be done. + + const Kinematics & kinematics = interaction->Kine(); + double bjx = kinematics.x(); + + double a = TMath::Power(bjx,2.) / TMath::Max(Q2val, fLowQ2CutoffF1F2); + double c = (1. + 4. * kNucleonMass2 * a) / (1.+r); + + fF3 = f * xF3val/bjx; + fF2 = f * F2val; + fF1 = fF2 * 0.5*c/bjx; + fF5 = fF2/bjx; // Albright-Jarlskog relation + fF4 = 0.; // Nucl.Phys.B 84, 467 (1975) + } + else { + double a = TMath::Power(x,2.) / TMath::Max(Q2val, fLowQ2CutoffF1F2); + double c = (1. + 4. * kNucleonMass2 * a) / (1.+r); + + fF3 = f * H * xF3val / x; + fF2 = f * F2val; + fF1 = fF2 * 0.5 * c / x; + fF5 = fF2 / x; // Albright-Jarlskog relation + fF4 = 0.; // Nucl.Phys.B 84, 467 (1975) + } #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ LOG("DISSF", pDEBUG) From 6221061948a135b2dcbb30cbac245eff6c92c24f Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Wed, 3 Dec 2025 15:23:56 +0100 Subject: [PATCH 27/68] setting chi modified only for F1 F2 calculations --- config/BYStrucFunc2021.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/config/BYStrucFunc2021.xml b/config/BYStrucFunc2021.xml index 605c87249f..0b885082f9 100644 --- a/config/BYStrucFunc2021.xml +++ b/config/BYStrucFunc2021.xml @@ -80,7 +80,7 @@ WeinbergAngle double No true true 0.8 - + true From 149bd0685d08188b41857965da365d533ac36b41 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Wed, 10 Dec 2025 16:13:57 +0100 Subject: [PATCH 28/68] commenting out low energy transfer factor --- .../DeepInelastic/XSection/BYStrucFunc2021.cxx | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx index 212ce76dc0..e3be402911 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx @@ -138,10 +138,13 @@ void BYStrucFunc2021::KFactors(const Interaction * interaction, double q0 = this->q0(interaction); double q02 = TMath::Power(q0,2); - double KLW = 1 ; + double KLW = 1. ; + + // The KLW is only important in the SIS region (W<2GeV). As we scale it altogether with the RES region, we do not need this factor. + // It also causes issues at low-W. After discussing with A. Bodek, we agreed to comment this part out. // The K factor blows up at q0 = 0. A. Bodek recomends to use it until W = 1.1 GeV - double W = interaction->Kine().W(); - //if ( W > 1.1 ) KLW = ( q02 + fCvLW ) / q02; + // double W = interaction->Kine().W(); + // if ( W > 1.1 ) KLW = ( q02 + fCvLW ) / q02; kuv = KLW * (1.-GD2)*(myQ2+fCv2U)/(myQ2+fCv1U); // K - u(valence) kdv = KLW * (1.-GD2)*(myQ2+fCv2D)/(myQ2+fCv1D); // K - d(valence) @@ -165,6 +168,9 @@ void BYStrucFunc2021::KAxialFactors (const Interaction * interaction, double & k bool is_nu = pdg::IsNeutrino ( probe_pdgc ); bool is_nubar = pdg::IsAntiNeutrino ( probe_pdgc ); double KLW = 1; + // The KLW is only important in the SIS region (W<2GeV). As we scale it altogether with the RES region, we do not need this factor. + // It also causes issues at low-W. After discussing with A. Bodek, we agreed to comment this part out. + /* if( q02 > 0 ) { // TO DECIDE ! It breaks the e- implementation if( is_nu ) { @@ -172,10 +178,12 @@ void BYStrucFunc2021::KAxialFactors (const Interaction * interaction, double & k } else if( is_nubar ) { KLW = ( q02 + fCaLW_nubar ) / q02 ; } - } + } // double check it only affects valance factors: //kvalance *= KLW ; + */ } + //____________________________________________________________________________ double BYStrucFunc2021::H(const Interaction * interaction) const { // Overrides QPMDISStrucFuncBase::H() function to compute the correction of the BY 2021 update From d77ded80b9d9544e1835c7ec9576b540f2207b52 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Thu, 11 Dec 2025 10:57:58 +0100 Subject: [PATCH 29/68] adding scaling of sea and valance quark contributions --- src/Physics/DeepInelastic/XSection/BYPDF.cxx | 10 +++++++++- src/Physics/DeepInelastic/XSection/BYPDF.h | 2 ++ 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/src/Physics/DeepInelastic/XSection/BYPDF.cxx b/src/Physics/DeepInelastic/XSection/BYPDF.cxx index 0af663fea6..a2c0811282 100644 --- a/src/Physics/DeepInelastic/XSection/BYPDF.cxx +++ b/src/Physics/DeepInelastic/XSection/BYPDF.cxx @@ -95,6 +95,13 @@ PDF_t BYPDF::AllPDFs(double x, double q2) const double dv = uncorrected_pdfs.dval; double ds = uncorrected_pdfs.dsea; + // we increase the up and down quark sea by 5%, and decrease the up and down valence quarks + // such that the sum of quark and antiquark distributions remain the same. + uv -= 2 * fUpScale * uncorrected_pdfs.uval; + dv -= 2 * fDownScale * uncorrected_pdfs.dval; + us *= ( 1 + fUpScale ) ; + ds *= ( 1 + fUpScale ) ; + // compute correction factor delta(d/u) double delta = this->DeltaDU(x); #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ @@ -168,7 +175,8 @@ void BYPDF::LoadConfig(void) GetParam( "BY-X0", fX0 ) ; GetParam( "BY-X1", fX1 ) ; GetParam( "BY-X2", fX2 ) ; - + GetParamDef( "BY-UpScale", fUpScale, 0. ); + GetParamDef( "BY-DownScale", fDownScale, 0. ); GetParam( "PDF-Q2min", fQ2min ) ; // get the base PDF model (typically GRV9* LO) diff --git a/src/Physics/DeepInelastic/XSection/BYPDF.h b/src/Physics/DeepInelastic/XSection/BYPDF.h index 6061def401..127f67f9f5 100644 --- a/src/Physics/DeepInelastic/XSection/BYPDF.h +++ b/src/Physics/DeepInelastic/XSection/BYPDF.h @@ -61,6 +61,8 @@ class BYPDF : public PDFModelI { double fX0; ///< correction param X0 double fX1; ///< correction param X1 double fX2; ///< correction param X2 + double fUpScale; ///< Scaling parameter to scale up contributions. Default 1. + double fDownScale; ///< Scaling parameter to scale down contributions. Default 1. double fQ2min; ///< min. Q2 for PDF evaluation }; From 8ec34a9ca096d30f4181413c8825bbe5ec37569d Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Thu, 11 Dec 2025 11:07:25 +0100 Subject: [PATCH 30/68] adding parameter in xml configuration --- config/BYPDF.xml | 7 +++++++ config/BYStrucFunc2021.xml | 3 ++- 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/config/BYPDF.xml b/config/BYPDF.xml index 72ca650251..2d75efb54b 100644 --- a/config/BYPDF.xml +++ b/config/BYPDF.xml @@ -26,6 +26,13 @@ Uncorr-PDF-Set alg No Uncorrected PDF model + + + 0.05 + 0.05 + + + diff --git a/config/BYStrucFunc2021.xml b/config/BYStrucFunc2021.xml index 0b885082f9..7581d6c583 100644 --- a/config/BYStrucFunc2021.xml +++ b/config/BYStrucFunc2021.xml @@ -47,7 +47,7 @@ WeinbergAngle double No CKM,Masses,WeakInt - genie::BYPDF/Default + genie::BYPDF/BY2021 0.621 0.380 @@ -81,6 +81,7 @@ WeinbergAngle double No true 0.8 true + From 2f87eef2a8dbc9d28eac3c567c8743fcb6e84c30 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Mon, 5 Jan 2026 12:05:15 +0100 Subject: [PATCH 31/68] adding vector mass as a free parameter --- src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx | 8 +++++++- src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h | 2 ++ 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx index e3be402911..96119a81b1 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx @@ -62,6 +62,10 @@ void BYStrucFunc2021::Configure(string param_set) //____________________________________________________________________________ void BYStrucFunc2021::ReadBYParams(void) { + // vector mass + GetParam( "EL-Mv",fMv ) ; + fMv2 = TMath::Power(fMv,2); + // Get the Bodek-Yang model parameters A,B,Csea,Cv1,Cv2 from the config. // registry and set some private data members so as not to accessing the // registry at every calculation. @@ -89,6 +93,8 @@ void BYStrucFunc2021::ReadBYParams(void) //____________________________________________________________________________ void BYStrucFunc2021::Init(void) { + fMv = 0; + fMv2 = 0; fA = 0; fB = 0; fCsU = 0; @@ -133,7 +139,7 @@ void BYStrucFunc2021::KFactors(const Interaction * interaction, // u(valence), d(valence), u(sea), d(sea), s(sea); double myQ2 = this->Q2(interaction); - double GD = 1. / TMath::Power(1.+myQ2/0.71, 2); // p elastic form factor + double GD = 1. / TMath::Power(1.+myQ2/fMv2, 2); // p elastic form factor double GD2 = TMath::Power(GD,2); double q0 = this->q0(interaction); diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h index 269445da58..3714b13b8b 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h @@ -56,6 +56,8 @@ class BYStrucFunc2021 : public QPMDISStrucFuncBase { double NuclMod(const Interaction * interaction) const; // Bodek-Yang model-specific parameters + double fMv; ///< Vector Mass + double fMv2; ///< Vector Mass Squared double fA; ///< better scaling var parameter A double fB; ///< better scaling var parameter B double fCvLW; ///< C low-nu vector paramter From 054ed58a82c35f5d8a20da41a109b989bea8bbe3 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Mon, 5 Jan 2026 12:22:40 +0100 Subject: [PATCH 32/68] removing scale up and down in dedicated configuration --- config/BYPDF.xml | 10 ++-------- config/BYStrucFunc2021.xml | 6 +++++- config/QPMDISPXSec.xml | 10 ++++++++++ src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx | 2 +- 4 files changed, 18 insertions(+), 10 deletions(-) diff --git a/config/BYPDF.xml b/config/BYPDF.xml index 2d75efb54b..cbda33ab15 100644 --- a/config/BYPDF.xml +++ b/config/BYPDF.xml @@ -14,6 +14,8 @@ BY-X1 double No Bodek-Yang param X1 BY-X2 double No Bodek-Yang param X2 PDF-Q2min double Yes Min Q2 for PDF evaluation from "Uncorr-PDF-Set" register Uncorr-PDF-Set alg No Uncorrected PDF model +BY-UpScale double Yes Scaling the up distribution. Default 0. +By-DownScale double Yes Scaling the down distribution. Default 0. --> @@ -26,14 +28,6 @@ Uncorr-PDF-Set alg No Uncorrected PDF model - - - 0.05 - 0.05 - - - - genie::LHAPDF5/GRVLO diff --git a/config/BYStrucFunc2021.xml b/config/BYStrucFunc2021.xml index 7581d6c583..81ecf5fb1c 100644 --- a/config/BYStrucFunc2021.xml +++ b/config/BYStrucFunc2021.xml @@ -83,6 +83,10 @@ WeinbergAngle double No true - + + + genie::BYPDF/Default + + diff --git a/config/QPMDISPXSec.xml b/config/QPMDISPXSec.xml index 4c1b4d94b4..ef571b6fa1 100644 --- a/config/QPMDISPXSec.xml +++ b/config/QPMDISPXSec.xml @@ -33,4 +33,14 @@ WeinbergAngle double No 1. 1. + + + WeakInt + genie::DISXSec/Default + genie::BYStrucFunc2021/Test + 1. + 1. + 1. + + diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx index 96119a81b1..f2833aa3ca 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx @@ -63,7 +63,7 @@ void BYStrucFunc2021::Configure(string param_set) void BYStrucFunc2021::ReadBYParams(void) { // vector mass - GetParam( "EL-Mv",fMv ) ; + GetParamDef( "EL-Mv",fMv, 0.84 ) ; fMv2 = TMath::Power(fMv,2); // Get the Bodek-Yang model parameters A,B,Csea,Cv1,Cv2 from the config. From d866f2b3e4c81259d771f8120fac368730b9a176 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Tue, 6 Jan 2026 21:24:49 +0100 Subject: [PATCH 33/68] removing unused configuration --- config/QPMDISPXSec.xml | 9 --------- src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx | 6 +++--- 2 files changed, 3 insertions(+), 12 deletions(-) diff --git a/config/QPMDISPXSec.xml b/config/QPMDISPXSec.xml index ef571b6fa1..2d72da2e6f 100644 --- a/config/QPMDISPXSec.xml +++ b/config/QPMDISPXSec.xml @@ -34,13 +34,4 @@ WeinbergAngle double No 1. - - WeakInt - genie::DISXSec/Default - genie::BYStrucFunc2021/Test - 1. - 1. - 1. - - diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx index f2833aa3ca..7fe1c3aa51 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx @@ -142,13 +142,13 @@ void BYStrucFunc2021::KFactors(const Interaction * interaction, double GD = 1. / TMath::Power(1.+myQ2/fMv2, 2); // p elastic form factor double GD2 = TMath::Power(GD,2); - double q0 = this->q0(interaction); - double q02 = TMath::Power(q0,2); - double KLW = 1. ; // The KLW is only important in the SIS region (W<2GeV). As we scale it altogether with the RES region, we do not need this factor. // It also causes issues at low-W. After discussing with A. Bodek, we agreed to comment this part out. // The K factor blows up at q0 = 0. A. Bodek recomends to use it until W = 1.1 GeV + double KLW = 1. ; + //double q0 = this->q0(interaction); + //double q02 = TMath::Power(q0,2); // double W = interaction->Kine().W(); // if ( W > 1.1 ) KLW = ( q02 + fCvLW ) / q02; From f1e477159f825f65956e0793386a5d467748bd42 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Tue, 13 Jan 2026 14:50:22 +0100 Subject: [PATCH 34/68] adding correct target-mass-corrected scaling variable --- .../DeepInelastic/XSection/BYStrucFunc2021.cxx | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx index 7fe1c3aa51..1329dd0a94 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx @@ -280,15 +280,20 @@ double BYStrucFunc2021::NuclMod(const Interaction * interaction) const { double xv3 = xv2 * xv; double xv4 = xv3 * xv; double xv5 = xv4 * xv; - - double f = 1.; // first factor goes from free nucleons to deuterium if(A >= 2) { f= 0.985*(1.+0.422*xv - 2.745*xv2 + 7.570*xv3 - 10.335*xv4 + 5.422*xv5); } - double chiTM = this->ScalingVar(interaction); + + // Computing target-mass-corrected scaling variable + double Q2 = kinematics.Q2(); + double M = kNucleonMass; + double M2 = pow( M, 2); + double chi_TM = 2 * xv / (1+sqrt(1+4*xv2*M2/Q2) ); + double f = 1.; + // 2nd factor goes from deuterium to iso-scalar iron if(A > 2) { // Fermi motion is already accounted for at high chiTM. To avoid double-counting, we limit the chiTM range: From bbd6d268137d515a64e80f0268903df572b6fec4 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Tue, 13 Jan 2026 15:08:14 +0100 Subject: [PATCH 35/68] fixing typo --- config/BYPDF.xml | 1 - config/BYStrucFunc2021.xml | 79 +++++++++---------- config/QPMDISPXSec.xml | 9 +++ .../XSection/BYStrucFunc2021.cxx | 2 +- 4 files changed, 48 insertions(+), 43 deletions(-) diff --git a/config/BYPDF.xml b/config/BYPDF.xml index cbda33ab15..7d61b01732 100644 --- a/config/BYPDF.xml +++ b/config/BYPDF.xml @@ -29,7 +29,6 @@ By-DownScale double Yes Scaling the down distribution. Default - genie::LHAPDF5/GRVLO diff --git a/config/BYStrucFunc2021.xml b/config/BYStrucFunc2021.xml index 81ecf5fb1c..86449125aa 100644 --- a/config/BYStrucFunc2021.xml +++ b/config/BYStrucFunc2021.xml @@ -44,49 +44,46 @@ WeinbergAngle double No - + CKM,Masses,WeakInt + genie::BYPDF/Default + 0.621 + 0.380 + 0.218 + 0.369 + 0.561 + 0.417 + 0.264 + 0.341 + 0.323 + 0.561 + 0.55 + 0.30 + 0.75 + 0.436 + 0.654 + 0.914 + 0.296 + -0.374 + 0.165 + + true + true + true + 0.8 + true - genie::BYPDF/BY2021 - - 0.621 - 0.380 - 0.218 - 0.369 - 0.561 - 0.417 - 0.264 - 0.341 - 0.323 - 0.561 - 0.55 - 0.30 - 0.75 - 0.436 - 0.654 - 0.914 - 0.296 - -0.374 - 0.165 - - true - true - true - 0.8 - true - - - - - genie::BYPDF/Default - + + + false + diff --git a/config/QPMDISPXSec.xml b/config/QPMDISPXSec.xml index 2d72da2e6f..804a7ef374 100644 --- a/config/QPMDISPXSec.xml +++ b/config/QPMDISPXSec.xml @@ -33,5 +33,14 @@ WeinbergAngle double No 1. 1. + + + WeakInt + genie::DISXSec/Default + genie::BYStrucFunc2021/Test + 1. + 1. + 1. + diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx index 1329dd0a94..e6cd3e0ab2 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx @@ -291,7 +291,7 @@ double BYStrucFunc2021::NuclMod(const Interaction * interaction) const { double Q2 = kinematics.Q2(); double M = kNucleonMass; double M2 = pow( M, 2); - double chi_TM = 2 * xv / (1+sqrt(1+4*xv2*M2/Q2) ); + double chiTM = 2 * xv / (1+sqrt(1+4*xv2*M2/Q2) ); double f = 1.; // 2nd factor goes from deuterium to iso-scalar iron From 3efbd810f7088e5ab4872e291ab0bb4dfcc8861d Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Wed, 14 Jan 2026 10:54:20 +0100 Subject: [PATCH 36/68] Nuclear correction returns correct value for deuterium --- .../XSection/BYStrucFunc2021.cxx | 126 +++++++++--------- 1 file changed, 61 insertions(+), 65 deletions(-) diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx index e6cd3e0ab2..544c3c495d 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx @@ -1,13 +1,13 @@ //____________________________________________________________________________ /* - Copyright (c) 2003-2023, The GENIE Collaboration - For the full text of the license visit http://copyright.genie-mc.org + Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org - Júlia Tena Vidal - Tel Aviv University + Júlia Tena Vidal + Tel Aviv University - Costas Andreopoulos - University of Liverpool + Costas Andreopoulos + University of Liverpool */ //____________________________________________________________________________ @@ -26,13 +26,13 @@ using namespace genie::constants; //____________________________________________________________________________ BYStrucFunc2021::BYStrucFunc2021() : -QPMDISStrucFuncBase("genie::BYStrucFunc2021") + QPMDISStrucFuncBase("genie::BYStrucFunc2021") { this->Init(); } //____________________________________________________________________________ BYStrucFunc2021::BYStrucFunc2021(string config): -QPMDISStrucFuncBase("genie::BYStrucFunc2021", config) + QPMDISStrucFuncBase("genie::BYStrucFunc2021", config) { this->Init(); } @@ -44,11 +44,11 @@ BYStrucFunc2021::~BYStrucFunc2021() //____________________________________________________________________________ void BYStrucFunc2021::Configure(const Registry & config) { -// Overload Algorithm::Configure() to read the config. registry and set -// private data members. -// QPMDISStrucFuncBase::Configure() creates the owned PDF object that gets -// configured with the specified PDFModelI -// For the ReadBYParams() method see below + // Overload Algorithm::Configure() to read the config. registry and set + // private data members. + // QPMDISStrucFuncBase::Configure() creates the owned PDF object that gets + // configured with the specified PDFModelI + // For the ReadBYParams() method see below QPMDISStrucFuncBase::Configure(config); this->ReadBYParams(); @@ -66,10 +66,10 @@ void BYStrucFunc2021::ReadBYParams(void) GetParamDef( "EL-Mv",fMv, 0.84 ) ; fMv2 = TMath::Power(fMv,2); -// Get the Bodek-Yang model parameters A,B,Csea,Cv1,Cv2 from the config. -// registry and set some private data members so as not to accessing the -// registry at every calculation. -// + // Get the Bodek-Yang model parameters A,B,Csea,Cv1,Cv2 from the config. + // registry and set some private data members so as not to accessing the + // registry at every calculation. + // GetParam( "BY-A", fA ) ; GetParam( "BY-B", fB ) ; GetParam( "BY-CsU", fCsU ) ; @@ -118,7 +118,7 @@ void BYStrucFunc2021::Init(void) //____________________________________________________________________________ double BYStrucFunc2021::ScalingVar(const Interaction * interaction, double Mf ) const { -// Overrides QPMDISStrucFuncBase::ScalingVar() to compute the BY scaling var + // Overrides QPMDISStrucFuncBase::ScalingVar() to compute the BY scaling var const Kinematics & kine = interaction->Kine(); double x = kine.x(); @@ -133,10 +133,10 @@ double BYStrucFunc2021::ScalingVar(const Interaction * interaction, double Mf ) } //____________________________________________________________________________ void BYStrucFunc2021::KFactors(const Interaction * interaction, - double & kuv, double & kdv, double & kus, double & kds, double & kss ) const + double & kuv, double & kdv, double & kus, double & kds, double & kss ) const { -// Overrides QPMDISStrucFuncBase::KFactors() to compute the BY K factors for -// u(valence), d(valence), u(sea), d(sea), s(sea); + // Overrides QPMDISStrucFuncBase::KFactors() to compute the BY K factors for + // u(valence), d(valence), u(sea), d(sea), s(sea); double myQ2 = this->Q2(interaction); double GD = 1. / TMath::Power(1.+myQ2/fMv2, 2); // p elastic form factor @@ -177,17 +177,17 @@ void BYStrucFunc2021::KAxialFactors (const Interaction * interaction, double & k // The KLW is only important in the SIS region (W<2GeV). As we scale it altogether with the RES region, we do not need this factor. // It also causes issues at low-W. After discussing with A. Bodek, we agreed to comment this part out. /* - if( q02 > 0 ) { + if( q02 > 0 ) { // TO DECIDE ! It breaks the e- implementation if( is_nu ) { - KLW = ( q02 + fCaLW_nu ) / q02 ; + KLW = ( q02 + fCaLW_nu ) / q02 ; } else if( is_nubar ) { - KLW = ( q02 + fCaLW_nubar ) / q02 ; + KLW = ( q02 + fCaLW_nubar ) / q02 ; } } - // double check it only affects valance factors: - //kvalance *= KLW ; - */ + // double check it only affects valance factors: + //kvalance *= KLW ; + */ } //____________________________________________________________________________ @@ -262,8 +262,8 @@ double BYStrucFunc2021::R(const Interaction * interaction) const { } //____________________________________________________________________________ double BYStrucFunc2021::NuclMod(const Interaction * interaction) const { -// Nuclear modification to Fi -// The scaling variable can be overwritten to include corrections + // Nuclear modification to Fi + // The scaling variable can be overwritten to include corrections if( interaction->TestBit(kIAssumeFreeNucleon) ) return 1.0; if( interaction->TestBit(kINoNuclearCorrection) ) return 1.0; @@ -271,48 +271,44 @@ double BYStrucFunc2021::NuclMod(const Interaction * interaction) const { double f = 1.; if(fIncludeNuclMod) { const Target & tgt = interaction->InitState().Tgt(); - const Kinematics & kinematics = interaction->Kine(); - double x = kinematics.x(); - int A = tgt.A(); + const Kinematics & kinematics = interaction->Kine(); + double x = kinematics.x(); + int A = tgt.A(); - double xv = TMath::Min(0.75, TMath::Max(0.05, x)); - double xv2 = xv * xv; - double xv3 = xv2 * xv; - double xv4 = xv3 * xv; - double xv5 = xv4 * xv; + double xv = TMath::Min(0.75, TMath::Max(0.05, x)); + double xv2 = xv * xv; + double xv3 = xv2 * xv; + double xv4 = xv3 * xv; + double xv5 = xv4 * xv; - // first factor goes from free nucleons to deuterium - if(A >= 2) { - f= 0.985*(1.+0.422*xv - 2.745*xv2 + 7.570*xv3 - 10.335*xv4 + 5.422*xv5); - } - - - // Computing target-mass-corrected scaling variable - double Q2 = kinematics.Q2(); - double M = kNucleonMass; - double M2 = pow( M, 2); - double chiTM = 2 * xv / (1+sqrt(1+4*xv2*M2/Q2) ); - double f = 1.; - - // 2nd factor goes from deuterium to iso-scalar iron - if(A > 2) { - // Fermi motion is already accounted for at high chiTM. To avoid double-counting, we limit the chiTM range: - if(! interaction->TestBit(kIAssumeFreeNucleon) && chiTM > 0.65 ) chiTM = 0.65; - f *= (1.096 - 0.38*chiTM - 0.3 * TMath::Exp(-23*chiTM) + 8 * pow(chiTM,15) ) ; - } + // first factor goes from free nucleons to deuterium + if(A >= 2) { + f*= 0.985*(1.+0.422*xv - 2.745*xv2 + 7.570*xv3 - 10.335*xv4 + 5.422*xv5); + } + + // Computing target-mass-corrected scaling variable + double Q2 = kinematics.Q2(); + double M = kNucleonMass; + double M2 = pow( M, 2); + double chiTM = 2 * xv / (1+sqrt(1+4*xv2*M2/Q2) ); + + // 2nd factor goes from deuterium to iso-scalar iron + if(A > 2) { + // Fermi motion is already accounted for at high chiTM. To avoid double-counting, we limit the chiTM range: + if(! interaction->TestBit(kIAssumeFreeNucleon) && chiTM > 0.65 ) chiTM = 0.65; + f *= (1.096 - 0.38*chiTM - 0.3 * TMath::Exp(-23*chiTM) + 8 * pow(chiTM,15) ) ; + } - if( A == 79 || A == 82 ) { // Gold nd Lead F2(Au,Pb)/F2(Fe) - f *= (0.932 + 2.461 * chiTM - 24.23*pow(chiTM,2) + 101.03*pow(chiTM,3) - 203.47*pow(chiTM,4) + 193.85*pow(chiTM,5) - 69.82*pow(chiTM,6)); - } + if( A == 79 || A == 82 ) { // Gold nd Lead F2(Au,Pb)/F2(Fe) + f *= (0.932 + 2.461 * chiTM - 24.23*pow(chiTM,2) + 101.03*pow(chiTM,3) - 203.47*pow(chiTM,4) + 193.85*pow(chiTM,5) - 69.82*pow(chiTM,6)); + } - if( A == 12 ) { // F2(Fe)/F2(C) - f /= ( 0.919 + 1.844*chiTM - 12.73*pow(chiTM,2) + 36.89*pow(chiTM,3) - 46.77*pow(chiTM,4) + 21.22*pow(chiTM,5)); - } + if( A == 12 ) { // F2(Fe)/F2(C) + f /= ( 0.919 + 1.844*chiTM - 12.73*pow(chiTM,2) + 36.89*pow(chiTM,3) - 46.77*pow(chiTM,4) + 21.22*pow(chiTM,5)); + } - return f; - #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ - LOG("DISSF", pDEBUG) << "Nuclear factor for x of " << x << " = " << f; + LOG("DISSF", pDEBUG) << "Nuclear factor for x of " << x << " = " << f; #endif } From 33a9b814b11ccbf969ae16e4372ea8f85e287dd2 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Thu, 15 Jan 2026 10:51:38 +0100 Subject: [PATCH 37/68] fixing A for Lead and gold --- .../DeepInelastic/XSection/BYStrucFunc2021.cxx | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx index 544c3c495d..e476544833 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx @@ -291,22 +291,23 @@ double BYStrucFunc2021::NuclMod(const Interaction * interaction) const { double M = kNucleonMass; double M2 = pow( M, 2); double chiTM = 2 * xv / (1+sqrt(1+4*xv2*M2/Q2) ); - + // Fermi motion is already accounted for at high chiTM. To avoid double-counting, we limit the chiTM range: + if( chiTM > 0.65 ) chiTM = 0.65; + // 2nd factor goes from deuterium to iso-scalar iron if(A > 2) { - // Fermi motion is already accounted for at high chiTM. To avoid double-counting, we limit the chiTM range: - if(! interaction->TestBit(kIAssumeFreeNucleon) && chiTM > 0.65 ) chiTM = 0.65; - f *= (1.096 - 0.38*chiTM - 0.3 * TMath::Exp(-23*chiTM) + 8 * pow(chiTM,15) ) ; + f *= (1.096 - 0.38*chiTM - 0.3 * TMath::Exp(-23*chiTM) + 8 * pow(chiTM,15) ) ; } - if( A == 79 || A == 82 ) { // Gold nd Lead F2(Au,Pb)/F2(Fe) + if( A == 197 || A == 208 ) { // Gold nd Lead F2(Au,Pb)/F2(Fe) f *= (0.932 + 2.461 * chiTM - 24.23*pow(chiTM,2) + 101.03*pow(chiTM,3) - 203.47*pow(chiTM,4) + 193.85*pow(chiTM,5) - 69.82*pow(chiTM,6)); + std::cout << " correction for gold " << (0.932 + 2.461 * chiTM - 24.23*pow(chiTM,2) + 101.03*pow(chiTM,3) - 203.47*pow(chiTM,4) + 193.85*pow(chiTM,5) - 69.82*pow(chiTM,6)) << std::endl; } if( A == 12 ) { // F2(Fe)/F2(C) f /= ( 0.919 + 1.844*chiTM - 12.73*pow(chiTM,2) + 36.89*pow(chiTM,3) - 46.77*pow(chiTM,4) + 21.22*pow(chiTM,5)); } - + #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ LOG("DISSF", pDEBUG) << "Nuclear factor for x of " << x << " = " << f; #endif From 38c166d7e6fc48ba21d9deaa0145cc1a0114061b Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Fri, 16 Jan 2026 11:33:12 +0100 Subject: [PATCH 38/68] adding H factor correction in F3 --- .../DeepInelastic/XSection/BYStrucFunc2021.cxx | 18 +++++++++++++++--- .../XSection/QPMDISStrucFuncBase.cxx | 4 ++-- 2 files changed, 17 insertions(+), 5 deletions(-) diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx index e476544833..f65d89ba2b 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx @@ -285,11 +285,23 @@ double BYStrucFunc2021::NuclMod(const Interaction * interaction) const { if(A >= 2) { f*= 0.985*(1.+0.422*xv - 2.745*xv2 + 7.570*xv3 - 10.335*xv4 + 5.422*xv5); } - + + // Implementing parameterization as a function of A from PhysRevD.49.4348 + double xv6 = xv5 * xv; + double xv7 = xv6 * xv; + double xv8 = xv7 * xv; + double alpha = -0.070 + 2.189 * xv - 24.667 * xv2 + 145.291 * xv3 - 497.237 * xv4 +1013.129 * xv5 - 1208.393 * xv6 +775.767 * xv7 - 205.872 * xv8; + double lnC = 0.017 +0.018 * TMath::Log(xv) + 0.005 * pow( TMath::Log(xv),2 ); + double C = TMath::Exp( lnC ) ; + f *= C * pow( A, alpha ) ; + + /* // Computing target-mass-corrected scaling variable double Q2 = kinematics.Q2(); double M = kNucleonMass; double M2 = pow( M, 2); + + double chiTM = 2 * xv / (1+sqrt(1+4*xv2*M2/Q2) ); // Fermi motion is already accounted for at high chiTM. To avoid double-counting, we limit the chiTM range: if( chiTM > 0.65 ) chiTM = 0.65; @@ -301,13 +313,13 @@ double BYStrucFunc2021::NuclMod(const Interaction * interaction) const { if( A == 197 || A == 208 ) { // Gold nd Lead F2(Au,Pb)/F2(Fe) f *= (0.932 + 2.461 * chiTM - 24.23*pow(chiTM,2) + 101.03*pow(chiTM,3) - 203.47*pow(chiTM,4) + 193.85*pow(chiTM,5) - 69.82*pow(chiTM,6)); - std::cout << " correction for gold " << (0.932 + 2.461 * chiTM - 24.23*pow(chiTM,2) + 101.03*pow(chiTM,3) - 203.47*pow(chiTM,4) + 193.85*pow(chiTM,5) - 69.82*pow(chiTM,6)) << std::endl; } if( A == 12 ) { // F2(Fe)/F2(C) f /= ( 0.919 + 1.844*chiTM - 12.73*pow(chiTM,2) + 36.89*pow(chiTM,3) - 46.77*pow(chiTM,4) + 21.22*pow(chiTM,5)); } - + */ + #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ LOG("DISSF", pDEBUG) << "Nuclear factor for x of " << x << " = " << f; #endif diff --git a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx index ed5b17ac14..5698be6141 100644 --- a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx +++ b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx @@ -375,9 +375,9 @@ void QPMDISStrucFuncBase::Calculate(const Interaction * interaction) const double a = TMath::Power(bjx,2.) / TMath::Max(Q2val, fLowQ2CutoffF1F2); double c = (1. + 4. * kNucleonMass2 * a) / (1.+r); - fF3 = f * xF3val/bjx; + fF3 = f * H * xF3val / bjx; fF2 = f * F2val; - fF1 = fF2 * 0.5*c/bjx; + fF1 = fF2 * 0.5 * c / bjx; fF5 = fF2/bjx; // Albright-Jarlskog relation fF4 = 0.; // Nucl.Phys.B 84, 467 (1975) } From 3dbe6814b0eb9f2a6ff4edcab3350a21a743e4a6 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Fri, 16 Jan 2026 12:14:52 +0100 Subject: [PATCH 39/68] using BY paper param of EMC effect --- src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx index f65d89ba2b..935da70170 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx @@ -286,6 +286,7 @@ double BYStrucFunc2021::NuclMod(const Interaction * interaction) const { f*= 0.985*(1.+0.422*xv - 2.745*xv2 + 7.570*xv3 - 10.335*xv4 + 5.422*xv5); } + /* // Implementing parameterization as a function of A from PhysRevD.49.4348 double xv6 = xv5 * xv; double xv7 = xv6 * xv; @@ -294,8 +295,8 @@ double BYStrucFunc2021::NuclMod(const Interaction * interaction) const { double lnC = 0.017 +0.018 * TMath::Log(xv) + 0.005 * pow( TMath::Log(xv),2 ); double C = TMath::Exp( lnC ) ; f *= C * pow( A, alpha ) ; + */ - /* // Computing target-mass-corrected scaling variable double Q2 = kinematics.Q2(); double M = kNucleonMass; @@ -318,7 +319,6 @@ double BYStrucFunc2021::NuclMod(const Interaction * interaction) const { if( A == 12 ) { // F2(Fe)/F2(C) f /= ( 0.919 + 1.844*chiTM - 12.73*pow(chiTM,2) + 36.89*pow(chiTM,3) - 46.77*pow(chiTM,4) + 21.22*pow(chiTM,5)); } - */ #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ LOG("DISSF", pDEBUG) << "Nuclear factor for x of " << x << " = " << f; From 9c985d7535593e1a9e2109e9308e63052213d635 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Thu, 5 Feb 2026 10:42:09 +0100 Subject: [PATCH 40/68] adding new algorithms for DIS nuclear modification --- .../NuclearState/BYDISNuclearModel.cxx | 110 ++++++++++++++++++ src/Physics/NuclearState/BYDISNuclearModel.h | 44 +++++++ src/Physics/NuclearState/DISNuclearModelI.cxx | 33 ++++++ src/Physics/NuclearState/DISNuclearModelI.h | 55 +++++++++ .../NuclearState/JGDISNuclearModel.cxx | 90 ++++++++++++++ 5 files changed, 332 insertions(+) create mode 100644 src/Physics/NuclearState/BYDISNuclearModel.cxx create mode 100644 src/Physics/NuclearState/BYDISNuclearModel.h create mode 100644 src/Physics/NuclearState/DISNuclearModelI.cxx create mode 100644 src/Physics/NuclearState/DISNuclearModelI.h create mode 100644 src/Physics/NuclearState/JGDISNuclearModel.cxx diff --git a/src/Physics/NuclearState/BYDISNuclearModel.cxx b/src/Physics/NuclearState/BYDISNuclearModel.cxx new file mode 100644 index 0000000000..fcd0bf576a --- /dev/null +++ b/src/Physics/NuclearState/BYDISNuclearModel.cxx @@ -0,0 +1,110 @@ + //____________________________________________________________________________ +/* + Copyright (c) 2003-2025, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + + J. Tena Vidal + Universitat de Valencia + + For the class documentation see the corresponding header file. + +*/ +//____________________________________________________________________________ + +#include "Physics/NuclearState/BYDISNuclearModel.h" + +using std::ostringstream; +using namespace genie; + +//____________________________________________________________________________ + +BYDISNuclearModel::BYDISNuclearModel() : + DISNuclearModelI("genie::BYDISNuclearModel"){} + +//____________________________________________________________________________ + +BYDISNuclearModel::BYDISNuclearModel(string config) : + DISNuclearModelI("genie::BYDISNuclearModel"){} + +//____________________________________________________________________________ + +double BYDISNuclearModel::DISACorrection (const Interaction & interaction) const { + double f = 1.; + + // Nuclear modification to Fi + // The scaling variable can be overwritten to include corrections + + if( interaction.TestBit(kIAssumeFreeNucleon) ) return f; + if( interaction.TestBit(kINoNuclearCorrection) ) return f; + + const Target & tgt = interaction.InitState().Tgt(); + const Kinematics & kinematics = interaction.Kine(); + double x = kinematics.x(); + int A = tgt.A(); + + double xv = TMath::Min(0.75, TMath::Max(0.05, x)); + double xv2 = xv * xv; + double xv3 = xv2 * xv; + double xv4 = xv3 * xv; + double xv5 = xv4 * xv; + + // first factor goes from free nucleons to deuterium + if(A >= 2) { + f*= 0.985*(1.+0.422*xv - 2.745*xv2 + 7.570*xv3 - 10.335*xv4 + 5.422*xv5); + } + + // Computing target-mass-corrected scaling variable + double Q2 = kinematics.Q2(); + double M = kNucleonMass; + double M2 = pow( M, 2); + + double chiTM = 2 * xv / (1+sqrt(1+4*xv2*M2/Q2) ); + // Fermi motion is already accounted for at high chiTM. To avoid double-counting, we limit the chiTM range: + if( chiTM > 0.65 ) chiTM = 0.65; + + // 2nd factor goes from deuterium to iso-scalar iron + if(A > 2) { + f *= (1.096 - 0.38*chiTM - 0.3 * TMath::Exp(-23*chiTM) + 8 * pow(chiTM,15) ) ; + } + + if( A == 197 || A == 208 ) { // Gold nd Lead F2(Au,Pb)/F2(Fe) + f *= (0.932 + 2.461 * chiTM - 24.23*pow(chiTM,2) + 101.03*pow(chiTM,3) - 203.47*pow(chiTM,4) + 193.85*pow(chiTM,5) - 69.82*pow(chiTM,6)); + } + + if( A == 12 ) { // F2(Fe)/F2(C) + f /= ( 0.919 + 1.844*chiTM - 12.73*pow(chiTM,2) + 36.89*pow(chiTM,3) - 46.77*pow(chiTM,4) + 21.22*pow(chiTM,5)); + } + +#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ + LOG("DISSF", pDEBUG) << "Nuclear factor for x of " << x << " = " << f; +#endif + + return f; + +} + +//____________________________________________________________________________ + +void BYDISNuclearModel::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} + +//____________________________________________________________________________ + +void BYDISNuclearModel::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} + +//____________________________________________________________________________ + +void BYDISNuclearModel::LoadConfig(void) +{ + //this->GetParam("SetName", fSetName ); + // this->GetParam("MemberID", fMemberID); + +} diff --git a/src/Physics/NuclearState/BYDISNuclearModel.h b/src/Physics/NuclearState/BYDISNuclearModel.h new file mode 100644 index 0000000000..e365435d09 --- /dev/null +++ b/src/Physics/NuclearState/BYDISNuclearModel.h @@ -0,0 +1,44 @@ +//____________________________________________________________________________ +/*! + +\class genie::BYDISNuclearModel + +\brief Uses DIS Nuclear Model Correction from Bodek-Yang + Reference: arxiv.org/pdf/2108.09240 (Section 9) + +\author J. Tena Vidal + Universitat de Valencia + +\created February, 2026 + +\cpright Copyright (c) 2003-2027, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org +*/ +//____________________________________________________________________________ + +#ifndef _BY_DIS_NUCLEAR_MODEL_I_H_ +#define _BY_DIS_NUCLEAR_MODEL_I_H_ + +#include + +#include "Physics/NuclearState/DISNuclearModelI.h" + +namespace genie { + +class BYDISNuclearModel : public DISNuclearModelI { + +public: + BYDISNuclearModel(); + BYDISNuclearModel(string config); + virtual ~BYDISNuclearModel() {}; + + double DISACorrection (const Interaction & interaction) const ; + void Configure (const Registry & config); + void Configure (string config); + +private: + void LoadConfig (void); +}; + +} // genie namespace +#endif // _BY_DIS_NUCLEAR_MODEL_I_H_ diff --git a/src/Physics/NuclearState/DISNuclearModelI.cxx b/src/Physics/NuclearState/DISNuclearModelI.cxx new file mode 100644 index 0000000000..232a628c00 --- /dev/null +++ b/src/Physics/NuclearState/DISNuclearModelI.cxx @@ -0,0 +1,33 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2025, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + J. Tena Vidal + Universitat de Valencia +*/ +//____________________________________________________________________________ + +#include "Physics/NuclearState/DISNuclearModelI.h" + +using namespace genie; + +//____________________________________________________________________________ +DISNuclearModelI::DISNuclearModelI() : +Algorithm() +{ + +} +//____________________________________________________________________________ +DISNuclearModelI::DISNuclearModelI(string name) : +Algorithm(name) +{ + +} +//____________________________________________________________________________ +DISNuclearModelI::DISNuclearModelI(string name, string config) : +Algorithm(name, config) +{ + +} +//____________________________________________________________________________ diff --git a/src/Physics/NuclearState/DISNuclearModelI.h b/src/Physics/NuclearState/DISNuclearModelI.h new file mode 100644 index 0000000000..21d7312f0e --- /dev/null +++ b/src/Physics/NuclearState/DISNuclearModelI.h @@ -0,0 +1,55 @@ +//____________________________________________________________________________ +/*! + +\class genie::DISNuclearModelI + +\brief Pure abstract base class. + Defines the EMCModelI interface to be implemented by any physics + model describing the ratio of l-A DIS cross-sections with respect to + the free nucleon calculation + +\author J. Tena Vidal + Universitat de Valencia + +\created February, 2026 + +\cpright Copyright (c) 2003-2027, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org +*/ +//____________________________________________________________________________ + +#ifndef _DIS_NUCLEAR_MODEL_I_H_ +#define _DIS_NUCLEAR_MODEL_I_H_ + +#include +#include "Framework/Algorithm/AlgConfigPool.h" +#include "Framework/Conventions/Constants.h" +#include "Framework/Conventions/Controls.h" +#include "Framework/Conventions/Units.h" +#include "Framework/Interaction/Interaction.h" +#include "Framework/Interaction/Target.h" +#include "Framework/Interaction/Kinematics.h" + +using std::ostringstream; +using namespace genie; +using namespace genie::constants; +using namespace genie::controls; + +namespace genie { + +class DISNuclearModelI : public Algorithm { + +public: + virtual ~DISNuclearModelI() {}; + virtual double DISACorrection (const Interaction & interaction) const = 0; + +protected: + + DISNuclearModelI(); + DISNuclearModelI(string name); + DISNuclearModelI(string name, string config); + +}; + +} // genie namespace +#endif // _DIS_NUCLEAR_MODEL_I_H_ diff --git a/src/Physics/NuclearState/JGDISNuclearModel.cxx b/src/Physics/NuclearState/JGDISNuclearModel.cxx new file mode 100644 index 0000000000..d7e312ae86 --- /dev/null +++ b/src/Physics/NuclearState/JGDISNuclearModel.cxx @@ -0,0 +1,90 @@ + //____________________________________________________________________________ +/* + Copyright (c) 2003-2025, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + + J. Tena Vidal + Universitat de Valencia + + For the class documentation see the corresponding header file. + +*/ +//____________________________________________________________________________ + +#include "Physics/NuclearState/JGDISNuclearModel.h" + +using std::ostringstream; +using namespace genie; + +//____________________________________________________________________________ + +JGDISNuclearModel::JGDISNuclearModel() : + DISNuclearModelI("genie::JGDISNuclearModel"){} + +//____________________________________________________________________________ + +JGDISNuclearModel::JGDISNuclearModel(string config) : + DISNuclearModelI("genie::JGDISNuclearModel"){} + +//____________________________________________________________________________ + +double JGDISNuclearModel::DISACorrection (const Interaction & interaction) const { + double f = 1.; + + // Nuclear modification to Fi + // The scaling variable can be overwritten to include corrections + + if( interaction.TestBit(kIAssumeFreeNucleon) ) return f; + if( interaction.TestBit(kINoNuclearCorrection) ) return f; + + const Target & tgt = interaction.InitState().Tgt(); + const Kinematics & kinematics = interaction.Kine(); + double x = kinematics.x(); + int A = tgt.A(); + + double xv = TMath::Min(0.75, TMath::Max(0.05, x)); + double xv2 = xv * xv; + double xv3 = xv2 * xv; + double xv4 = xv3 * xv; + double xv5 = xv4 * xv; + double xv6 = xv5 * xv; + double xv7 = xv6 * xv; + double xv8 = xv7 * xv; + double alpha = -0.070 + 2.189 * xv - 24.667 * xv2 + 145.291 * xv3 - 497.237 * xv4 + 1013.129 * xv5 - 1208.393 * xv6 +775.767 * xv7 - 205.872 * xv8; + double lnC = 0.017 +0.018 * TMath::Log(xv) + 0.005 * pow( TMath::Log(xv),2 ); + double C = TMath::Exp( lnC ) ; + f *= C * pow( A, alpha ) ; + +#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ + LOG("DISSF", pDEBUG) << "J. Gomez et. al. Nuclear factor for x of " << x << " = " << f; +#endif + + return f; + +} + +//____________________________________________________________________________ + +void JGDISNuclearModel::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} + +//____________________________________________________________________________ + +void JGDISNuclearModel::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} + +//____________________________________________________________________________ + +void JGDISNuclearModel::LoadConfig(void) +{ + //this->GetParam("SetName", fSetName ); + // this->GetParam("MemberID", fMemberID); + +} From ccf046abc612d5ff0a6f36a6880d572ba9beaf37 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Thu, 5 Feb 2026 10:55:34 +0100 Subject: [PATCH 41/68] adding xml files for new DIS nuclear models --- config/BYDISNuclearModel.xml | 20 ++++++++++++++++++++ config/JGDISNuclearModel.xml | 20 ++++++++++++++++++++ 2 files changed, 40 insertions(+) create mode 100644 config/BYDISNuclearModel.xml create mode 100644 config/JGDISNuclearModel.xml diff --git a/config/BYDISNuclearModel.xml b/config/BYDISNuclearModel.xml new file mode 100644 index 0000000000..4cd56a4f35 --- /dev/null +++ b/config/BYDISNuclearModel.xml @@ -0,0 +1,20 @@ + + + + + + + + + + + + + diff --git a/config/JGDISNuclearModel.xml b/config/JGDISNuclearModel.xml new file mode 100644 index 0000000000..410c1de3fa --- /dev/null +++ b/config/JGDISNuclearModel.xml @@ -0,0 +1,20 @@ + + + + + + + + + + + + + From 97a5e6fc812e60631183d97256088d1119b4d212 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Thu, 5 Feb 2026 11:10:49 +0100 Subject: [PATCH 42/68] adding nuclear model in xml --- config/QPMDISStrucFunc.xml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/config/QPMDISStrucFunc.xml b/config/QPMDISStrucFunc.xml index 3c25bbde2c..e2da40f561 100644 --- a/config/QPMDISStrucFunc.xml +++ b/config/QPMDISStrucFunc.xml @@ -8,6 +8,7 @@ Algorithm Configurable Parameters: Name Type Opt Comment Default .................................................................................................... PDF-Set alg No PDF model +DISNuclModel alg No DIS Nuclear Model Charm-Mass double No charm mass CommonParam[Masses] Charm-Prod-Off bool Yes charm production is turned off false PDF-Q2min double No min Q2 for PDF evaluation From "PDF-Set" @@ -28,7 +29,8 @@ WeinbergAngle double No CKM,Masses,WeakInt - genie::GRV98LO/Default + genie::GRV98LO/Default + genie::BYDISNuclearModel/Default BetheBlochModel.xml From 62cd8b55093f984cd66d8c83726dfdc1dccbec32 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Thu, 5 Feb 2026 14:28:44 +0100 Subject: [PATCH 45/68] update structure function with nuclear model --- config/BYStrucFunc.xml | 18 +----- config/BYStrucFunc2021.xml | 5 +- config/QPMDISPXSec.xml | 1 - .../NuclearState/BYDISNuclearModel.cxx | 59 ++++++------------- src/Physics/NuclearState/BYDISNuclearModel.h | 9 ++- src/Physics/NuclearState/DISNuclearModelI.h | 2 +- .../NuclearState/JGDISNuclearModel.cxx | 11 ++-- 7 files changed, 34 insertions(+), 71 deletions(-) diff --git a/config/BYStrucFunc.xml b/config/BYStrucFunc.xml index d73233e9cc..2d4cfa58ea 100644 --- a/config/BYStrucFunc.xml +++ b/config/BYStrucFunc.xml @@ -62,22 +62,10 @@ WeinbergAngle double No true false 0.8 - - - - - - 0.621 - 0.380 - 0.369 - 0.561 - 0.417 - 0.264 - 0.341 - 0.323 - - + + genie::BYDISNuclearModel21/Default + diff --git a/config/BYStrucFunc2021.xml b/config/BYStrucFunc2021.xml index 86449125aa..08ba63a907 100644 --- a/config/BYStrucFunc2021.xml +++ b/config/BYStrucFunc2021.xml @@ -78,8 +78,9 @@ WeinbergAngle double No true true 0.8 - true - + true + + genie::BYDISNuclearModel21/Default diff --git a/config/QPMDISPXSec.xml b/config/QPMDISPXSec.xml index f7662dd802..9bd788c6c8 100644 --- a/config/QPMDISPXSec.xml +++ b/config/QPMDISPXSec.xml @@ -30,7 +30,6 @@ WeinbergAngle double No WeakInt genie::DISXSec/Default genie::BYStrucFunc2021/Default - genie::BYDISNuclearModel21/Default 1. 1. 1. diff --git a/src/Physics/NuclearState/BYDISNuclearModel.cxx b/src/Physics/NuclearState/BYDISNuclearModel.cxx index fcd0bf576a..eee565672e 100644 --- a/src/Physics/NuclearState/BYDISNuclearModel.cxx +++ b/src/Physics/NuclearState/BYDISNuclearModel.cxx @@ -13,6 +13,7 @@ //____________________________________________________________________________ #include "Physics/NuclearState/BYDISNuclearModel.h" +#include "Physics/NuclearState/NuclearUtils.h" using std::ostringstream; using namespace genie; @@ -29,55 +30,29 @@ BYDISNuclearModel::BYDISNuclearModel(string config) : //____________________________________________________________________________ -double BYDISNuclearModel::DISACorrection (const Interaction & interaction) const { +double BYDISNuclearModel::DISACorrection (const Interaction * interaction) const { + if ( !interaction ) return 0; double f = 1.; // Nuclear modification to Fi // The scaling variable can be overwritten to include corrections - if( interaction.TestBit(kIAssumeFreeNucleon) ) return f; - if( interaction.TestBit(kINoNuclearCorrection) ) return f; - - const Target & tgt = interaction.InitState().Tgt(); - const Kinematics & kinematics = interaction.Kine(); - double x = kinematics.x(); - int A = tgt.A(); - - double xv = TMath::Min(0.75, TMath::Max(0.05, x)); - double xv2 = xv * xv; - double xv3 = xv2 * xv; - double xv4 = xv3 * xv; - double xv5 = xv4 * xv; - - // first factor goes from free nucleons to deuterium - if(A >= 2) { - f*= 0.985*(1.+0.422*xv - 2.745*xv2 + 7.570*xv3 - 10.335*xv4 + 5.422*xv5); - } - - // Computing target-mass-corrected scaling variable - double Q2 = kinematics.Q2(); - double M = kNucleonMass; - double M2 = pow( M, 2); - - double chiTM = 2 * xv / (1+sqrt(1+4*xv2*M2/Q2) ); - // Fermi motion is already accounted for at high chiTM. To avoid double-counting, we limit the chiTM range: - if( chiTM > 0.65 ) chiTM = 0.65; - - // 2nd factor goes from deuterium to iso-scalar iron - if(A > 2) { - f *= (1.096 - 0.38*chiTM - 0.3 * TMath::Exp(-23*chiTM) + 8 * pow(chiTM,15) ) ; - } - - if( A == 197 || A == 208 ) { // Gold nd Lead F2(Au,Pb)/F2(Fe) - f *= (0.932 + 2.461 * chiTM - 24.23*pow(chiTM,2) + 101.03*pow(chiTM,3) - 203.47*pow(chiTM,4) + 193.85*pow(chiTM,5) - 69.82*pow(chiTM,6)); - } - - if( A == 12 ) { // F2(Fe)/F2(C) - f /= ( 0.919 + 1.844*chiTM - 12.73*pow(chiTM,2) + 36.89*pow(chiTM,3) - 46.77*pow(chiTM,4) + 21.22*pow(chiTM,5)); - } + if( interaction->TestBit(kIAssumeFreeNucleon) ) return 1.0; + if( interaction->TestBit(kINoNuclearCorrection) ) return 1.0; + const Target & tgt = interaction->InitState().Tgt(); + + // The x used for computing the DIS Nuclear correction factor should be the + // experimental x, not the rescaled x or off-shell-rest-frame version of x + // (i.e. selected x). Since we do not have access to experimental x at this + // point in the calculation, just use selected x. + const Kinematics & kine = interaction->Kine(); + double x = kine.x(); + int A = tgt.A(); + f = utils::nuclear::DISNuclFactor(x,A); + #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ - LOG("DISSF", pDEBUG) << "Nuclear factor for x of " << x << " = " << f; + LOG("DISSF", pDEBUG) << "Nuclear factor for x of " << x << " = " << f; #endif return f; diff --git a/src/Physics/NuclearState/BYDISNuclearModel.h b/src/Physics/NuclearState/BYDISNuclearModel.h index e365435d09..07a968fe97 100644 --- a/src/Physics/NuclearState/BYDISNuclearModel.h +++ b/src/Physics/NuclearState/BYDISNuclearModel.h @@ -4,7 +4,6 @@ \class genie::BYDISNuclearModel \brief Uses DIS Nuclear Model Correction from Bodek-Yang - Reference: arxiv.org/pdf/2108.09240 (Section 9) \author J. Tena Vidal Universitat de Valencia @@ -16,8 +15,8 @@ */ //____________________________________________________________________________ -#ifndef _BY_DIS_NUCLEAR_MODEL_I_H_ -#define _BY_DIS_NUCLEAR_MODEL_I_H_ +#ifndef _BY_DIS_NUCLEAR_MODEL_H_ +#define _BY_DIS_NUCLEAR_MODEL_H_ #include @@ -32,7 +31,7 @@ class BYDISNuclearModel : public DISNuclearModelI { BYDISNuclearModel(string config); virtual ~BYDISNuclearModel() {}; - double DISACorrection (const Interaction & interaction) const ; + double DISACorrection (const Interaction * interaction) const ; void Configure (const Registry & config); void Configure (string config); @@ -41,4 +40,4 @@ class BYDISNuclearModel : public DISNuclearModelI { }; } // genie namespace -#endif // _BY_DIS_NUCLEAR_MODEL_I_H_ +#endif // _BY_DIS_NUCLEAR_MODEL_H_ diff --git a/src/Physics/NuclearState/DISNuclearModelI.h b/src/Physics/NuclearState/DISNuclearModelI.h index 21d7312f0e..d4b40ae69c 100644 --- a/src/Physics/NuclearState/DISNuclearModelI.h +++ b/src/Physics/NuclearState/DISNuclearModelI.h @@ -41,7 +41,7 @@ class DISNuclearModelI : public Algorithm { public: virtual ~DISNuclearModelI() {}; - virtual double DISACorrection (const Interaction & interaction) const = 0; + virtual double DISACorrection (const Interaction * interaction) const = 0; protected: diff --git a/src/Physics/NuclearState/JGDISNuclearModel.cxx b/src/Physics/NuclearState/JGDISNuclearModel.cxx index d7e312ae86..04d9ac6291 100644 --- a/src/Physics/NuclearState/JGDISNuclearModel.cxx +++ b/src/Physics/NuclearState/JGDISNuclearModel.cxx @@ -29,17 +29,18 @@ JGDISNuclearModel::JGDISNuclearModel(string config) : //____________________________________________________________________________ -double JGDISNuclearModel::DISACorrection (const Interaction & interaction) const { +double JGDISNuclearModel::DISACorrection (const Interaction * interaction) const { + if ( !interaction ) return 0; double f = 1.; // Nuclear modification to Fi // The scaling variable can be overwritten to include corrections - if( interaction.TestBit(kIAssumeFreeNucleon) ) return f; - if( interaction.TestBit(kINoNuclearCorrection) ) return f; + if( interaction->TestBit(kIAssumeFreeNucleon) ) return f; + if( interaction->TestBit(kINoNuclearCorrection) ) return f; - const Target & tgt = interaction.InitState().Tgt(); - const Kinematics & kinematics = interaction.Kine(); + const Target & tgt = interaction->InitState().Tgt(); + const Kinematics & kinematics = interaction->Kine(); double x = kinematics.x(); int A = tgt.A(); From 753027424d89cf9ada2ae31c1abac10bd1f62eda Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Thu, 5 Feb 2026 14:30:25 +0100 Subject: [PATCH 46/68] update xml --- config/BYStrucFunc.xml | 50 +++++++++++++++++++++--------------------- 1 file changed, 25 insertions(+), 25 deletions(-) diff --git a/config/BYStrucFunc.xml b/config/BYStrucFunc.xml index 2d4cfa58ea..51ed66ca70 100644 --- a/config/BYStrucFunc.xml +++ b/config/BYStrucFunc.xml @@ -36,36 +36,36 @@ WeinbergAngle double No - + CKM,Masses,WeakInt - genie::BYPDF/Default + genie::BYPDF/Default - 0.538 - 0.305 - 0.363 - 0.621 - 0.291 - 0.189 - 0.202 - 0.255 + 0.538 + 0.305 + 0.363 + 0.621 + 0.291 + 0.189 + 0.202 + 0.255 - - true - true - false - 0.8 + + true + true + false + 0.8 - - genie::BYDISNuclearModel21/Default - + + genie::BYDISNuclearModel/Default + From 28ee254dce907a45c7e5581304b489c8af5b738f Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Thu, 5 Feb 2026 14:31:50 +0100 Subject: [PATCH 47/68] adding structure --- config/BYDISNuclearModel21.xml | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) create mode 100644 config/BYDISNuclearModel21.xml diff --git a/config/BYDISNuclearModel21.xml b/config/BYDISNuclearModel21.xml new file mode 100644 index 0000000000..4cd56a4f35 --- /dev/null +++ b/config/BYDISNuclearModel21.xml @@ -0,0 +1,20 @@ + + + + + + + + + + + + + From 6e3afcb61e8e2f162dce117dc70edcf72a88a656 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Thu, 5 Feb 2026 14:57:58 +0100 Subject: [PATCH 48/68] editing link file --- src/Physics/NuclearState/LinkDef.h | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/Physics/NuclearState/LinkDef.h b/src/Physics/NuclearState/LinkDef.h index c9207f1ca8..742f427860 100644 --- a/src/Physics/NuclearState/LinkDef.h +++ b/src/Physics/NuclearState/LinkDef.h @@ -25,5 +25,8 @@ #pragma link C++ class genie::SRCNuclearRecoil; #pragma link C++ class genie::SecondNucleonEmissionI; #pragma link C++ class genie::SpectralFunction2p2h; +#pragma link C++ class BYDISNuclearModel; +#pragma link C++ class BYDISNuclearModel21; +#pragma link C++ class JGDISNuclearModel; #endif From 5b0ca9bc9adc4af58af6a44230814197a022ad07 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Fri, 13 Feb 2026 16:04:25 +0100 Subject: [PATCH 49/68] adding missing files --- config/QPMDISPXSec.xml | 11 -- .../NuclearState/BYDISNuclearModel21.cxx | 111 ++++++++++++++++++ .../NuclearState/BYDISNuclearModel21.h | 44 +++++++ src/Physics/NuclearState/JGDISNuclearModel.h | 44 +++++++ 4 files changed, 199 insertions(+), 11 deletions(-) create mode 100644 src/Physics/NuclearState/BYDISNuclearModel21.cxx create mode 100644 src/Physics/NuclearState/BYDISNuclearModel21.h create mode 100644 src/Physics/NuclearState/JGDISNuclearModel.h diff --git a/config/QPMDISPXSec.xml b/config/QPMDISPXSec.xml index 9bd788c6c8..6cfccbb50c 100644 --- a/config/QPMDISPXSec.xml +++ b/config/QPMDISPXSec.xml @@ -20,7 +20,6 @@ WeinbergAngle double No genie::BYStrucFunc/Default genie::DISXSec/Default - genie::BYDISNuclearModel/Default 1.032 1.032 1. @@ -35,14 +34,4 @@ WeinbergAngle double No 1. - - WeakInt - genie::DISXSec/Default - genie::BYStrucFunc2021/Test - genie::JGDISNuclearModel/Default - 1. - 1. - 1. - - diff --git a/src/Physics/NuclearState/BYDISNuclearModel21.cxx b/src/Physics/NuclearState/BYDISNuclearModel21.cxx new file mode 100644 index 0000000000..d69e6ce0f4 --- /dev/null +++ b/src/Physics/NuclearState/BYDISNuclearModel21.cxx @@ -0,0 +1,111 @@ + //____________________________________________________________________________ +/* + Copyright (c) 2003-2025, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + + J. Tena Vidal + Universitat de Valencia + + For the class documentation see the corresponding header file. + +*/ +//____________________________________________________________________________ + +#include "Physics/NuclearState/BYDISNuclearModel21.h" + +using std::ostringstream; +using namespace genie; + +//____________________________________________________________________________ + +BYDISNuclearModel21::BYDISNuclearModel21() : + DISNuclearModelI("genie::BYDISNuclearModel21"){} + +//____________________________________________________________________________ + +BYDISNuclearModel21::BYDISNuclearModel21(string config) : + DISNuclearModelI("genie::BYDISNuclearModel21"){} + +//____________________________________________________________________________ + +double BYDISNuclearModel21::DISACorrection (const Interaction * interaction) const { + if ( !interaction ) return 0; + double f = 1.; + + // Nuclear modification to Fi + // The scaling variable can be overwritten to include corrections + + if( interaction->TestBit(kIAssumeFreeNucleon) ) return f; + if( interaction->TestBit(kINoNuclearCorrection) ) return f; + + const Target & tgt = interaction->InitState().Tgt(); + const Kinematics & kinematics = interaction->Kine(); + double x = kinematics.x(); + int A = tgt.A(); + + double xv = TMath::Min(0.75, TMath::Max(0.05, x)); + double xv2 = xv * xv; + double xv3 = xv2 * xv; + double xv4 = xv3 * xv; + double xv5 = xv4 * xv; + + // first factor goes from free nucleons to deuterium + if(A >= 2) { + f*= 0.985*(1.+0.422*xv - 2.745*xv2 + 7.570*xv3 - 10.335*xv4 + 5.422*xv5); + } + + // Computing target-mass-corrected scaling variable + double Q2 = kinematics.Q2(); + double M = kNucleonMass; + double M2 = pow( M, 2); + + double chiTM = 2 * xv / (1+sqrt(1+4*xv2*M2/Q2) ); + // Fermi motion is already accounted for at high chiTM. To avoid double-counting, we limit the chiTM range: + if( chiTM > 0.65 ) chiTM = 0.65; + + // 2nd factor goes from deuterium to iso-scalar iron + if(A > 2) { + f *= (1.096 - 0.38*chiTM - 0.3 * TMath::Exp(-23*chiTM) + 8 * pow(chiTM,15) ) ; + } + + if( A == 197 || A == 208 ) { // Gold nd Lead F2(Au,Pb)/F2(Fe) + f *= (0.932 + 2.461 * chiTM - 24.23*pow(chiTM,2) + 101.03*pow(chiTM,3) - 203.47*pow(chiTM,4) + 193.85*pow(chiTM,5) - 69.82*pow(chiTM,6)); + } + + if( A == 12 ) { // F2(Fe)/F2(C) + f /= ( 0.919 + 1.844*chiTM - 12.73*pow(chiTM,2) + 36.89*pow(chiTM,3) - 46.77*pow(chiTM,4) + 21.22*pow(chiTM,5)); + } + +#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ + LOG("DISSF", pDEBUG) << "Nuclear factor for x of " << x << " = " << f; +#endif + + return f; + +} + +//____________________________________________________________________________ + +void BYDISNuclearModel21::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} + +//____________________________________________________________________________ + +void BYDISNuclearModel21::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} + +//____________________________________________________________________________ + +void BYDISNuclearModel21::LoadConfig(void) +{ + //this->GetParam("SetName", fSetName ); + // this->GetParam("MemberID", fMemberID); + +} diff --git a/src/Physics/NuclearState/BYDISNuclearModel21.h b/src/Physics/NuclearState/BYDISNuclearModel21.h new file mode 100644 index 0000000000..bf3a60a69a --- /dev/null +++ b/src/Physics/NuclearState/BYDISNuclearModel21.h @@ -0,0 +1,44 @@ +//____________________________________________________________________________ +/*! + +\class genie::BYDISNuclearModel21 + +\brief Uses DIS Nuclear Model Correction from the Updated Bodek-Yang Model + Reference: arxiv.org/pdf/2108.09240 (Section 9) + +\author J. Tena Vidal + Universitat de Valencia + +\created February, 2026 + +\cpright Copyright (c) 2003-2027, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org +*/ +//____________________________________________________________________________ + +#ifndef _BY_DIS_NUCLEAR_MODEL_21_H_ +#define _BY_DIS_NUCLEAR_MODEL_21_H_ + +#include + +#include "Physics/NuclearState/DISNuclearModelI.h" + +namespace genie { + +class BYDISNuclearModel21 : public DISNuclearModelI { + +public: + BYDISNuclearModel21(); + BYDISNuclearModel21(string config); + virtual ~BYDISNuclearModel21() {}; + + double DISACorrection (const Interaction * interaction) const ; + void Configure (const Registry & config); + void Configure (string config); + +private: + void LoadConfig (void); +}; + +} // genie namespace +#endif // _BY_DIS_NUCLEAR_MODEL_21_H_ diff --git a/src/Physics/NuclearState/JGDISNuclearModel.h b/src/Physics/NuclearState/JGDISNuclearModel.h new file mode 100644 index 0000000000..375cda7999 --- /dev/null +++ b/src/Physics/NuclearState/JGDISNuclearModel.h @@ -0,0 +1,44 @@ +//____________________________________________________________________________ +/*! + +\class genie::BYDISNuclearModel + +\brief Uses DIS Nuclear Model Correction from J.Gomez et. al. + Reference: journals.aps.org/prd/pdf/10.1103/PhysRevD.49.4348 + +\author J. Tena Vidal + Universitat de Valencia + +\created February, 2026 + +\cpright Copyright (c) 2003-2027, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org +*/ +//____________________________________________________________________________ + +#ifndef _JG_DIS_NUCLEAR_MODEL_H_ +#define _JG_DIS_NUCLEAR_MODEL_H_ + +#include + +#include "Physics/NuclearState/DISNuclearModelI.h" + +namespace genie { + +class JGDISNuclearModel : public DISNuclearModelI { + +public: + JGDISNuclearModel(); + JGDISNuclearModel(string config); + virtual ~JGDISNuclearModel() {}; + + double DISACorrection (const Interaction * interaction) const ; + void Configure (const Registry & config); + void Configure (string config); + +private: + void LoadConfig (void); +}; + +} // genie namespace +#endif // _JG_DIS_NUCLEAR_MODEL_H_ From a95aecd0200352cf0ba0aa3fc189545cd2e8aab5 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Tue, 17 Feb 2026 14:23:39 +0100 Subject: [PATCH 50/68] adding nuclear correction in xml file --- config/BYStrucFunc2021.xml | 12 +++++++++--- config/DISXSec.xml | 15 +++++++++++++++ config/QPMDISPXSec.xml | 4 ++-- 3 files changed, 26 insertions(+), 5 deletions(-) diff --git a/config/BYStrucFunc2021.xml b/config/BYStrucFunc2021.xml index 08ba63a907..e17448194a 100644 --- a/config/BYStrucFunc2021.xml +++ b/config/BYStrucFunc2021.xml @@ -79,12 +79,18 @@ WeinbergAngle double No true 0.8 true - + + genie::JGDISNuclearModel/Default + + + genie::BYDISNuclearModel21/Default - - false + + genie::BYDISNuclearModel/Default diff --git a/config/DISXSec.xml b/config/DISXSec.xml index 61bcd392c1..78a27b14ac 100644 --- a/config/DISXSec.xml +++ b/config/DISXSec.xml @@ -2,6 +2,16 @@ @@ -15,5 +25,10 @@ Configuration for the DISXSec cross section algorithm 0.0001 + + 50000000 + 0.00001 + true + diff --git a/config/QPMDISPXSec.xml b/config/QPMDISPXSec.xml index 6cfccbb50c..975f685d00 100644 --- a/config/QPMDISPXSec.xml +++ b/config/QPMDISPXSec.xml @@ -27,8 +27,8 @@ WeinbergAngle double No WeakInt - genie::DISXSec/Default - genie::BYStrucFunc2021/Default + genie::DISXSec/GRV982010 + genie::BYStrucFunc2021/Default 1. 1. 1. From 8d2a40c27b3c7a3bdb890273cdb6e10dc607743e Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Tue, 17 Feb 2026 14:42:21 +0100 Subject: [PATCH 51/68] adding changes to make DIS Nuclear correction optional --- src/Physics/DeepInelastic/XSection/DISXSec.cxx | 11 ++++++++--- src/Physics/DeepInelastic/XSection/DISXSec.h | 1 + 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/src/Physics/DeepInelastic/XSection/DISXSec.cxx b/src/Physics/DeepInelastic/XSection/DISXSec.cxx index dd8002194b..9d29f8a99d 100644 --- a/src/Physics/DeepInelastic/XSection/DISXSec.cxx +++ b/src/Physics/DeepInelastic/XSection/DISXSec.cxx @@ -77,7 +77,8 @@ double DISXSec::Integrate( // If yes, calculate the nuclear cross section based on that value. // XSecSplineList * xsl = XSecSplineList::Instance(); - if(init_state.Tgt().IsNucleus() && !xsl->IsEmpty() ) { + if(init_state.Tgt().IsNucleus() && !xsl->IsEmpty() && !fDISNuclCorr ) { + // Computes xsec from free nucleon calculation wo nuclear effects Interaction * interaction = new Interaction(*in); Target * target = interaction->InitStatePtr()->TgtPtr(); if(pdg::IsProton(nucpdgc)) { target->SetId(kPdgTgtFreeP); } @@ -104,7 +105,7 @@ double DISXSec::Integrate( // at any subsequent call. // bool precalc_bare_xsec = RunOpt::Instance()->BareXSecPreCalc(); - if(precalc_bare_xsec) { + if(precalc_bare_xsec && !fDISNuclCorr) { ///// CHECK LATER. I THINK NOT NEEDED THE CHANGE HERE Cache * cache = Cache::Instance(); Interaction * interaction = new Interaction(*in); string key = this->CacheBranchName(model,interaction); @@ -140,7 +141,7 @@ double DISXSec::Integrate( // Since nuclear corrections don't need to be included at this stage, all the // nuclear cross sections can be trivially built from the free nucleon ones. // - interaction->SetBit(kINoNuclearCorrection); + if(!fDISNuclCorr) interaction->SetBit(kINoNuclearCorrection); Range1D_t Wl = kps.WLim(); Range1D_t Q2l = kps.Q2Lim(); @@ -207,6 +208,10 @@ void DISXSec::LoadConfig(void) GetParam( "GVLD-Emin", fVldEmin) ; GetParam( "GVLD-Emax", fVldEmax) ; + // This is set to false by default for historical reasons. + // The new option set this to true which distorts sigma/E + // Set this to false to keep the historical behavior flat sigma/E + GetParamDef("DIS-NuclCorr", fDISNuclCorr, false); } //____________________________________________________________________________ void DISXSec::CacheFreeNucleonXSec( diff --git a/src/Physics/DeepInelastic/XSection/DISXSec.h b/src/Physics/DeepInelastic/XSection/DISXSec.h index e4e4d97992..1e2c8105ad 100644 --- a/src/Physics/DeepInelastic/XSection/DISXSec.h +++ b/src/Physics/DeepInelastic/XSection/DISXSec.h @@ -46,6 +46,7 @@ class DISXSec : public XSecIntegratorI { double fVldEmin; double fVldEmax; + bool fDISNuclCorr ; }; } // genie namespace From 80a8d4ebbddbfa0490d4b13916ae112e10a76340 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Tue, 17 Feb 2026 14:47:21 +0100 Subject: [PATCH 52/68] changing name of variable --- src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx | 4 ++-- src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx index 39321da9f0..adf1b28eea 100644 --- a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx +++ b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx @@ -122,7 +122,7 @@ void QPMDISStrucFuncBase::LoadConfig(void) //-- turn charm production off? GetParamDef( "Charm-Prod-Off", fCharmOff, false ) ; - fNuclMod = dynamic_cast(this->SubAlg("DISNuclModel")); + fDISNuclCorr = dynamic_cast(this->SubAlg("DISNuclModel")); //-- weinberg angle double thw ; @@ -353,7 +353,7 @@ void QPMDISStrucFuncBase::Calculate(const Interaction * interaction) const double Q2val = this->Q2 (interaction); double x = this->ScalingVar(interaction); - double f = fIncludeNuclMod ? fNuclMod->DISACorrection(interaction) : 1 ; + double f = fIncludeNuclMod ? fDISNuclCorr->DISACorrection(interaction) : 1 ; double r = this->R (interaction); // R ~ FL double H = fIncludeH ? this->H(interaction) : 1; diff --git a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h index 88350dd186..1a3384b988 100644 --- a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h +++ b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h @@ -91,7 +91,7 @@ class QPMDISStrucFuncBase : public DISStructureFuncModelI { bool fUse2016Corrections;///< Use 2016 SF relation corrections double fLowQ2CutoffF1F2; ///< Set min for relation between 2xF1 and F2 - const DISNuclearModelI * fNuclMod ; ///< model for nuclear factors (shadowing, anti-shadowing,...) + const DISNuclearModelI * fDISNuclCorr ; ///< model for nuclear factors (shadowing, anti-shadowing,...) mutable double fF1; mutable double fF2; From a99a4e2dc7abc4a1e8e5375bbf4e885990d3522b Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Tue, 17 Feb 2026 15:18:38 +0100 Subject: [PATCH 53/68] changing integrator configuration --- config/DISXSec.xml | 8 ++++---- config/KNOTunedQPMDISPXSec.xml | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/config/DISXSec.xml b/config/DISXSec.xml index 78a27b14ac..58ec4be45c 100644 --- a/config/DISXSec.xml +++ b/config/DISXSec.xml @@ -7,10 +7,10 @@ Algorithm Configurable Parameters: .......................................................................................................................... Name Type Opt Comment Default .......................................................................................................................... -gsl-integration-type string No -gsl-max-eval int No -gsl-min-eval int No -gsl-relative-tolerance double No +gsl-integration-type string Yes adaptative +gsl-max-eval int Yes 500000 +gsl-min-eval int Yes 10000 +gsl-relative-tolerance double Yes 1E-2 DIS-NuclCorr bool Yes Apply DIS Nuclear Correction false --> diff --git a/config/KNOTunedQPMDISPXSec.xml b/config/KNOTunedQPMDISPXSec.xml index 54c0d7af2d..25b07b4adb 100644 --- a/config/KNOTunedQPMDISPXSec.xml +++ b/config/KNOTunedQPMDISPXSec.xml @@ -27,7 +27,7 @@ Wcut double no Co NonResBackground - genie::DISXSec/Default + genie::DISXSec/GRV982010 genie::AGKYLowW2019/Default 1. genie::QPMDISPXSec/GRV982010 From e8e24bd53a0844df45b0552e1d88a45e49c40da0 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Wed, 25 Feb 2026 15:11:54 +0100 Subject: [PATCH 54/68] adding free parameters in DIS Nuclear Model --- config/BY00DISNuclearModel.xml | 17 ++++++++ config/BY21DISNuclearModel.xml | 78 ++++++++++++++++++++++++++++++++++ config/BYDISNuclearModel.xml | 20 --------- config/BYDISNuclearModel21.xml | 20 --------- 4 files changed, 95 insertions(+), 40 deletions(-) create mode 100644 config/BY00DISNuclearModel.xml create mode 100644 config/BY21DISNuclearModel.xml delete mode 100644 config/BYDISNuclearModel.xml delete mode 100644 config/BYDISNuclearModel21.xml diff --git a/config/BY00DISNuclearModel.xml b/config/BY00DISNuclearModel.xml new file mode 100644 index 0000000000..d48e28aedd --- /dev/null +++ b/config/BY00DISNuclearModel.xml @@ -0,0 +1,17 @@ + + + + + + + + + + + + + diff --git a/config/BY21DISNuclearModel.xml b/config/BY21DISNuclearModel.xml new file mode 100644 index 0000000000..ee05e38eac --- /dev/null +++ b/config/BY21DISNuclearModel.xml @@ -0,0 +1,78 @@ + + + + + + + + 0.985 + 1 + 0.422 + -2.745 + 7.570 + -10.335 + 5.422 + + 1.096 + -0.38 + -0.30 + -23.0 + 8.0 + + 0.919 + 1.844 + -12.73 + 36.89 + -46.77 + 21.22 + + 0.932 + 2.461 + -24.23 + 101.03 + -203.47 + 193.85 + -69.82 + + + + + + diff --git a/config/BYDISNuclearModel.xml b/config/BYDISNuclearModel.xml deleted file mode 100644 index 4cd56a4f35..0000000000 --- a/config/BYDISNuclearModel.xml +++ /dev/null @@ -1,20 +0,0 @@ - - - - - - - - - - - - - diff --git a/config/BYDISNuclearModel21.xml b/config/BYDISNuclearModel21.xml deleted file mode 100644 index 4cd56a4f35..0000000000 --- a/config/BYDISNuclearModel21.xml +++ /dev/null @@ -1,20 +0,0 @@ - - - - - - - - - - - - - From 943f5896d39c9a75645316b89735fbd734cf5918 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Wed, 25 Feb 2026 15:12:36 +0100 Subject: [PATCH 55/68] adding parameters instead of hardcoded values --- ...clearModel.cxx => BY00DISNuclearModel.cxx} | 24 +++--- ...ISNuclearModel.h => BY00DISNuclearModel.h} | 14 ++-- ...earModel21.cxx => BY21DISNuclearModel.cxx} | 58 +++++++++++---- .../NuclearState/BY21DISNuclearModel.h | 74 +++++++++++++++++++ .../NuclearState/BYDISNuclearModel21.h | 44 ----------- 5 files changed, 136 insertions(+), 78 deletions(-) rename src/Physics/NuclearState/{BYDISNuclearModel.cxx => BY00DISNuclearModel.cxx} (75%) rename src/Physics/NuclearState/{BYDISNuclearModel.h => BY00DISNuclearModel.h} (76%) rename src/Physics/NuclearState/{BYDISNuclearModel21.cxx => BY21DISNuclearModel.cxx} (52%) create mode 100644 src/Physics/NuclearState/BY21DISNuclearModel.h delete mode 100644 src/Physics/NuclearState/BYDISNuclearModel21.h diff --git a/src/Physics/NuclearState/BYDISNuclearModel.cxx b/src/Physics/NuclearState/BY00DISNuclearModel.cxx similarity index 75% rename from src/Physics/NuclearState/BYDISNuclearModel.cxx rename to src/Physics/NuclearState/BY00DISNuclearModel.cxx index eee565672e..e2c6bbc94c 100644 --- a/src/Physics/NuclearState/BYDISNuclearModel.cxx +++ b/src/Physics/NuclearState/BY00DISNuclearModel.cxx @@ -12,7 +12,7 @@ */ //____________________________________________________________________________ -#include "Physics/NuclearState/BYDISNuclearModel.h" +#include "Physics/NuclearState/BY00DISNuclearModel.h" #include "Physics/NuclearState/NuclearUtils.h" using std::ostringstream; @@ -20,17 +20,17 @@ using namespace genie; //____________________________________________________________________________ -BYDISNuclearModel::BYDISNuclearModel() : - DISNuclearModelI("genie::BYDISNuclearModel"){} +BY00DISNuclearModel::BY00DISNuclearModel() : + DISNuclearModelI("genie::BY00DISNuclearModel"){} //____________________________________________________________________________ -BYDISNuclearModel::BYDISNuclearModel(string config) : - DISNuclearModelI("genie::BYDISNuclearModel"){} +BY00DISNuclearModel::BY00DISNuclearModel(string config) : + DISNuclearModelI("genie::BY00DISNuclearModel"){} //____________________________________________________________________________ -double BYDISNuclearModel::DISACorrection (const Interaction * interaction) const { +double BY00DISNuclearModel::DISACorrection (const Interaction * interaction) const { if ( !interaction ) return 0; double f = 1.; @@ -61,7 +61,7 @@ double BYDISNuclearModel::DISACorrection (const Interaction * interaction) const //____________________________________________________________________________ -void BYDISNuclearModel::Configure(const Registry & config) +void BY00DISNuclearModel::Configure(const Registry & config) { Algorithm::Configure(config); this->LoadConfig(); @@ -69,7 +69,7 @@ void BYDISNuclearModel::Configure(const Registry & config) //____________________________________________________________________________ -void BYDISNuclearModel::Configure(string config) +void BY00DISNuclearModel::Configure(string config) { Algorithm::Configure(config); this->LoadConfig(); @@ -77,9 +77,9 @@ void BYDISNuclearModel::Configure(string config) //____________________________________________________________________________ -void BYDISNuclearModel::LoadConfig(void) +void BY00DISNuclearModel::LoadConfig(void) { - //this->GetParam("SetName", fSetName ); - // this->GetParam("MemberID", fMemberID); - + // In this version, parameters are hardcoded. + // This model is not advised, as a new version of this parameterization is available + // in BY21DISNuclearModel } diff --git a/src/Physics/NuclearState/BYDISNuclearModel.h b/src/Physics/NuclearState/BY00DISNuclearModel.h similarity index 76% rename from src/Physics/NuclearState/BYDISNuclearModel.h rename to src/Physics/NuclearState/BY00DISNuclearModel.h index 07a968fe97..c90d20a0da 100644 --- a/src/Physics/NuclearState/BYDISNuclearModel.h +++ b/src/Physics/NuclearState/BY00DISNuclearModel.h @@ -1,7 +1,7 @@ //____________________________________________________________________________ /*! -\class genie::BYDISNuclearModel +\class genie::BY00DISNuclearModel \brief Uses DIS Nuclear Model Correction from Bodek-Yang @@ -15,8 +15,8 @@ */ //____________________________________________________________________________ -#ifndef _BY_DIS_NUCLEAR_MODEL_H_ -#define _BY_DIS_NUCLEAR_MODEL_H_ +#ifndef _BY00_DIS_NUCLEAR_MODEL_H_ +#define _BY00_DIS_NUCLEAR_MODEL_H_ #include @@ -24,12 +24,12 @@ namespace genie { -class BYDISNuclearModel : public DISNuclearModelI { +class BY00DISNuclearModel : public DISNuclearModelI { public: - BYDISNuclearModel(); - BYDISNuclearModel(string config); - virtual ~BYDISNuclearModel() {}; + BY00DISNuclearModel(); + BY00DISNuclearModel(string config); + virtual ~BY00DISNuclearModel() {}; double DISACorrection (const Interaction * interaction) const ; void Configure (const Registry & config); diff --git a/src/Physics/NuclearState/BYDISNuclearModel21.cxx b/src/Physics/NuclearState/BY21DISNuclearModel.cxx similarity index 52% rename from src/Physics/NuclearState/BYDISNuclearModel21.cxx rename to src/Physics/NuclearState/BY21DISNuclearModel.cxx index d69e6ce0f4..82c30ecf5f 100644 --- a/src/Physics/NuclearState/BYDISNuclearModel21.cxx +++ b/src/Physics/NuclearState/BY21DISNuclearModel.cxx @@ -12,24 +12,24 @@ */ //____________________________________________________________________________ -#include "Physics/NuclearState/BYDISNuclearModel21.h" +#include "Physics/NuclearState/BY21DISNuclearModel.h" using std::ostringstream; using namespace genie; //____________________________________________________________________________ -BYDISNuclearModel21::BYDISNuclearModel21() : - DISNuclearModelI("genie::BYDISNuclearModel21"){} +BY21DISNuclearModel::BY21DISNuclearModel() : + DISNuclearModelI("genie::BY21DISNuclearModel"){} //____________________________________________________________________________ -BYDISNuclearModel21::BYDISNuclearModel21(string config) : - DISNuclearModelI("genie::BYDISNuclearModel21"){} +BY21DISNuclearModel::BY21DISNuclearModel(string config) : + DISNuclearModelI("genie::BY21DISNuclearModel"){} //____________________________________________________________________________ -double BYDISNuclearModel21::DISACorrection (const Interaction * interaction) const { +double BY21DISNuclearModel::DISACorrection (const Interaction * interaction) const { if ( !interaction ) return 0; double f = 1.; @@ -52,7 +52,7 @@ double BYDISNuclearModel21::DISACorrection (const Interaction * interaction) con // first factor goes from free nucleons to deuterium if(A >= 2) { - f*= 0.985*(1.+0.422*xv - 2.745*xv2 + 7.570*xv3 - 10.335*xv4 + 5.422*xv5); + f*= f2HScale*( f2Hf0 + f2Hf1*xv + f2Hf2*xv2 + f2Hf3*xv3 + f2Hf4*xv4 + f2Hf5*xv5); } // Computing target-mass-corrected scaling variable @@ -66,15 +66,15 @@ double BYDISNuclearModel21::DISACorrection (const Interaction * interaction) con // 2nd factor goes from deuterium to iso-scalar iron if(A > 2) { - f *= (1.096 - 0.38*chiTM - 0.3 * TMath::Exp(-23*chiTM) + 8 * pow(chiTM,15) ) ; + f *= ( fFef0 + fFef1 * chiTM + fFef2 * TMath::Exp( fFef3*chiTM ) + fFef4 * pow(chiTM,15) ) ; } if( A == 197 || A == 208 ) { // Gold nd Lead F2(Au,Pb)/F2(Fe) - f *= (0.932 + 2.461 * chiTM - 24.23*pow(chiTM,2) + 101.03*pow(chiTM,3) - 203.47*pow(chiTM,4) + 193.85*pow(chiTM,5) - 69.82*pow(chiTM,6)); + f *= ( fAuf0 + fAuf1 * chiTM + fAuf2 *pow(chiTM,2) + fAuf3*pow(chiTM,3) + fAuf4*pow(chiTM,4) + fAuf5*pow(chiTM,5) + fAuf6*pow(chiTM,6)); } if( A == 12 ) { // F2(Fe)/F2(C) - f /= ( 0.919 + 1.844*chiTM - 12.73*pow(chiTM,2) + 36.89*pow(chiTM,3) - 46.77*pow(chiTM,4) + 21.22*pow(chiTM,5)); + f /= ( fCf0 + fCf1*chiTM + fCf2*pow(chiTM,2) + fCf3*pow(chiTM,3) + fCf4*pow(chiTM,4) + fCf5*pow(chiTM,5)); } #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ @@ -87,7 +87,7 @@ double BYDISNuclearModel21::DISACorrection (const Interaction * interaction) con //____________________________________________________________________________ -void BYDISNuclearModel21::Configure(const Registry & config) +void BY21DISNuclearModel::Configure(const Registry & config) { Algorithm::Configure(config); this->LoadConfig(); @@ -95,7 +95,7 @@ void BYDISNuclearModel21::Configure(const Registry & config) //____________________________________________________________________________ -void BYDISNuclearModel21::Configure(string config) +void BY21DISNuclearModel::Configure(string config) { Algorithm::Configure(config); this->LoadConfig(); @@ -103,9 +103,37 @@ void BYDISNuclearModel21::Configure(string config) //____________________________________________________________________________ -void BYDISNuclearModel21::LoadConfig(void) +void BY21DISNuclearModel::LoadConfig(void) { - //this->GetParam("SetName", fSetName ); - // this->GetParam("MemberID", fMemberID); + GetParam( "BY21-NuclModel-2H-Scale", f2HScale ) ; + GetParam( "BY21-NuclModel-2H-f0", f2Hf0 ) ; + GetParam( "BY21-NuclModel-2H-f1", f2Hf1 ) ; + GetParam( "BY21-NuclModel-2H-f2", f2Hf2 ) ; + GetParam( "BY21-NuclModel-2H-f3", f2Hf3 ) ; + GetParam( "BY21-NuclModel-2H-f4", f2Hf4 ) ; + GetParam( "BY21-NuclModel-2H-f5", f2Hf5 ) ; + + GetParam( "BY21-NuclModel-Fe-f0", fFef0 ) ; + GetParam( "BY21-NuclModel-Fe-f1", fFef1 ) ; + GetParam( "BY21-NuclModel-Fe-f2", fFef2 ) ; + GetParam( "BY21-NuclModel-Fe-f3", fFef3 ) ; + GetParam( "BY21-NuclModel-Fe-f4", fFef4 ) ; + + GetParam( "BY21-NuclModel-C-f0", fCf0 ) ; + GetParam( "BY21-NuclModel-C-f1", fCf1 ) ; + GetParam( "BY21-NuclModel-C-f2", fCf2 ) ; + GetParam( "BY21-NuclModel-C-f3", fCf3 ) ; + GetParam( "BY21-NuclModel-C-f4", fCf4 ) ; + GetParam( "BY21-NuclModel-C-f5", fCf5 ) ; + GetParam( "BY21-NuclModel-C-f6", fCf6 ) ; + + GetParam( "BY21-NuclModel-Au-f0", fAuf0 ) ; + GetParam( "BY21-NuclModel-Au-f1", fAuf1 ) ; + GetParam( "BY21-NuclModel-Au-f2", fAuf2 ) ; + GetParam( "BY21-NuclModel-Au-f3", fAuf3 ) ; + GetParam( "BY21-NuclModel-Au-f4", fAuf4 ) ; + GetParam( "BY21-NuclModel-Au-f5", fAuf5 ) ; + GetParam( "BY21-NuclModel-Au-f6", fAuf6 ) ; + } diff --git a/src/Physics/NuclearState/BY21DISNuclearModel.h b/src/Physics/NuclearState/BY21DISNuclearModel.h new file mode 100644 index 0000000000..9f9c30905e --- /dev/null +++ b/src/Physics/NuclearState/BY21DISNuclearModel.h @@ -0,0 +1,74 @@ +//____________________________________________________________________________ +/*! + +\class genie::BY21DISNuclearModel + +\brief Uses DIS Nuclear Model Correction from the Updated Bodek-Yang Model + Reference: arxiv.org/pdf/2108.09240 (Section 9) + +\author J. Tena Vidal + Universitat de Valencia + +\created February, 2026 + +\cpright Copyright (c) 2003-2027, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org +*/ +//____________________________________________________________________________ + +#ifndef _BY_DIS_NUCLEAR_MODEL_21_H_ +#define _BY_DIS_NUCLEAR_MODEL_21_H_ + +#include + +#include "Physics/NuclearState/DISNuclearModelI.h" + +namespace genie { + +class BY21DISNuclearModel : public DISNuclearModelI { + +public: + BY21DISNuclearModel(); + BY21DISNuclearModel(string config); + virtual ~BY21DISNuclearModel() {}; + + double DISACorrection (const Interaction * interaction) const ; + void Configure (const Registry & config); + void Configure (string config); + +private: + void LoadConfig (void); + + double f2HScale ; //> BY21-NuclModel-2H-Scale + double f2Hf0 ; //> BY21-NuclModel-2H-f0 + double f2Hf1 ; //> BY21-NuclModel-2H-f1 + double f2Hf2 ; //> BY21-NuclModel-2H-f2 + double f2Hf3 ; //> BY21-NuclModel-2H-f3 + double f2Hf4 ; //> BY21-NuclModel-2H-f4 + double f2Hf5 ; //> BY21-NuclModel-2H-f5 + + double fFef0 ; //> BY21-NuclModel-Fe-f0 + double fFef1 ; //> BY21-NuclModel-Fe-f1 + double fFef2 ; //> BY21-NuclModel-Fe-f2 + double fFef3 ; //> BY21-NuclModel-Fe-f3 + double fFef4 ; //> BY21-NuclModel-Fe-f4 + + double fCf0 ; //> BY21-NuclModel-C-f0 + double fCf1 ; //> BY21-NuclModel-C-f1 + double fCf2 ; //> BY21-NuclModel-C-f2 + double fCf3 ; //> BY21-NuclModel-C-f3 + double fCf4 ; //> BY21-NuclModel-C-f4 + double fCf5 ; //> BY21-NuclModel-C-f5 + + double fAuf0 ; //> BY21-NuclModel-Au-f0 + double fAuf1 ; //> BY21-NuclModel-Au-f1 + double fAuf2 ; //> BY21-NuclModel-Au-f2 + double fAuf3 ; //> BY21-NuclModel-Au-f3 + double fAuf4 ; //> BY21-NuclModel-Au-f4 + double fAuf5 ; //> BY21-NuclModel-Au-f5 + double fAuf6 ; //> BY21-NuclModel-Au-f6 + +}; + +} // genie namespace +#endif // _BY_DIS_NUCLEAR_MODEL_21_H_ diff --git a/src/Physics/NuclearState/BYDISNuclearModel21.h b/src/Physics/NuclearState/BYDISNuclearModel21.h deleted file mode 100644 index bf3a60a69a..0000000000 --- a/src/Physics/NuclearState/BYDISNuclearModel21.h +++ /dev/null @@ -1,44 +0,0 @@ -//____________________________________________________________________________ -/*! - -\class genie::BYDISNuclearModel21 - -\brief Uses DIS Nuclear Model Correction from the Updated Bodek-Yang Model - Reference: arxiv.org/pdf/2108.09240 (Section 9) - -\author J. Tena Vidal - Universitat de Valencia - -\created February, 2026 - -\cpright Copyright (c) 2003-2027, The GENIE Collaboration - For the full text of the license visit http://copyright.genie-mc.org -*/ -//____________________________________________________________________________ - -#ifndef _BY_DIS_NUCLEAR_MODEL_21_H_ -#define _BY_DIS_NUCLEAR_MODEL_21_H_ - -#include - -#include "Physics/NuclearState/DISNuclearModelI.h" - -namespace genie { - -class BYDISNuclearModel21 : public DISNuclearModelI { - -public: - BYDISNuclearModel21(); - BYDISNuclearModel21(string config); - virtual ~BYDISNuclearModel21() {}; - - double DISACorrection (const Interaction * interaction) const ; - void Configure (const Registry & config); - void Configure (string config); - -private: - void LoadConfig (void); -}; - -} // genie namespace -#endif // _BY_DIS_NUCLEAR_MODEL_21_H_ From c2bdc9099669ac8e4a9e988e85eeac2f973f9fd2 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Wed, 25 Feb 2026 15:25:50 +0100 Subject: [PATCH 56/68] adding free parameters for JG --- config/JGDISNuclearModel.xml | 35 ++++++++++++++++++++++++++++++----- 1 file changed, 30 insertions(+), 5 deletions(-) diff --git a/config/JGDISNuclearModel.xml b/config/JGDISNuclearModel.xml index 410c1de3fa..8a622d7c9b 100644 --- a/config/JGDISNuclearModel.xml +++ b/config/JGDISNuclearModel.xml @@ -6,14 +6,39 @@ Configuration for the J. Gomez et. al. DIS Nuclear model journals.aps.org/prd/pdf/10.1103/PhysRevD.49.4348 Configurable Parameters: -......................................................................................... -Name Type Optional Comment Default -......................................................................................... +......................................................................................................... +Name Type Optional Comment +......................................................................................................... +JG-NuclModel-c0 double No +JG-NuclModel-c1 double No +JG-NuclModel-c2 double No + +JG-NuclModel-a0 double No +JG-NuclModel-a1 double No +JG-NuclModel-a2 double No +JG-NuclModel-a3 double No +JG-NuclModel-a4 double No +JG-NuclModel-a5 double No +JG-NuclModel-a6 double No +JG-NuclModel-a7 double No +JG-NuclModel-a8 double No --> - - + 0.017 + 0.018 + 0.005 + + -0.070 + 2.189 + -24.667 + 145.291 + -497.237 + 1013.129 + -1208.393 + 775.767 + -205.872 + From cda31b2e82b8af5c9f56a8b5935a036570e3ee71 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Wed, 25 Feb 2026 15:26:44 +0100 Subject: [PATCH 57/68] parameters not hardcoded --- .../NuclearState/JGDISNuclearModel.cxx | 19 ++++++++++++++----- src/Physics/NuclearState/JGDISNuclearModel.h | 13 +++++++++++++ 2 files changed, 27 insertions(+), 5 deletions(-) diff --git a/src/Physics/NuclearState/JGDISNuclearModel.cxx b/src/Physics/NuclearState/JGDISNuclearModel.cxx index 04d9ac6291..622b0dbe2b 100644 --- a/src/Physics/NuclearState/JGDISNuclearModel.cxx +++ b/src/Physics/NuclearState/JGDISNuclearModel.cxx @@ -52,8 +52,8 @@ double JGDISNuclearModel::DISACorrection (const Interaction * interaction) const double xv6 = xv5 * xv; double xv7 = xv6 * xv; double xv8 = xv7 * xv; - double alpha = -0.070 + 2.189 * xv - 24.667 * xv2 + 145.291 * xv3 - 497.237 * xv4 + 1013.129 * xv5 - 1208.393 * xv6 +775.767 * xv7 - 205.872 * xv8; - double lnC = 0.017 +0.018 * TMath::Log(xv) + 0.005 * pow( TMath::Log(xv),2 ); + double alpha = fAa0 + fAa1 * xv + fAa2 * xv2 + fAa3 * xv3 + fAa4 * xv4 + fAa5 * xv5 + fAa6 * xv6 + fAa7 * xv7 + fAa8 * xv8; + double lnC = fAc0 + fAc1 * TMath::Log(xv) + fAc2 * pow( TMath::Log(xv),2 ); double C = TMath::Exp( lnC ) ; f *= C * pow( A, alpha ) ; @@ -85,7 +85,16 @@ void JGDISNuclearModel::Configure(string config) void JGDISNuclearModel::LoadConfig(void) { - //this->GetParam("SetName", fSetName ); - // this->GetParam("MemberID", fMemberID); - + GetParam( "JG-NuclModel-c0", fAc0 ) ; + GetParam( "JG-NuclModel-c1", fAc1 ) ; + GetParam( "JG-NuclModel-c2", fAc2 ) ; + + GetParam( "JG-NuclModel-a0", fAa0 ) ; + GetParam( "JG-NuclModel-a1", fAa1 ) ; + GetParam( "JG-NuclModel-a2", fAa2 ) ; + GetParam( "JG-NuclModel-a3", fAa3 ) ; + GetParam( "JG-NuclModel-a4", fAa4 ) ; + GetParam( "JG-NuclModel-a5", fAa5 ) ; + GetParam( "JG-NuclModel-a6", fAa6 ) ; + GetParam( "JG-NuclModel-a7", fAa7 ) ; } diff --git a/src/Physics/NuclearState/JGDISNuclearModel.h b/src/Physics/NuclearState/JGDISNuclearModel.h index 375cda7999..69e765704d 100644 --- a/src/Physics/NuclearState/JGDISNuclearModel.h +++ b/src/Physics/NuclearState/JGDISNuclearModel.h @@ -38,6 +38,19 @@ class JGDISNuclearModel : public DISNuclearModelI { private: void LoadConfig (void); + + double fAc0; //> JG-NuclModel-c0 + double fAc1; //> JG-NuclModel-c1 + double fAc2; //> JG-NuclModel-c2 + + double fAa0; //> JG-NuclModel-a0 + double fAa1; //> JG-NuclModel-a1 + double fAa2; //> JG-NuclModel-a2 + double fAa3; //> JG-NuclModel-a3 + double fAa4; //> JG-NuclModel-a4 + double fAa5; //> JG-NuclModel-a5 + double fAa6; //> JG-NuclModel-a6 + double fAa7; //> JG-NuclModel-a7 }; } // genie namespace From 37f5b61aed90a2dc425b7bea159f9fe169ae407f Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Wed, 25 Feb 2026 15:38:49 +0100 Subject: [PATCH 58/68] fixing bugs --- src/Physics/DeepInelastic/XSection/DISXSec.cxx | 1 + src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx | 2 +- src/Physics/NuclearState/BY21DISNuclearModel.cxx | 1 - src/Physics/NuclearState/JGDISNuclearModel.cxx | 1 + src/Physics/NuclearState/JGDISNuclearModel.h | 1 + 5 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/Physics/DeepInelastic/XSection/DISXSec.cxx b/src/Physics/DeepInelastic/XSection/DISXSec.cxx index 9d29f8a99d..a65d2bcaf5 100644 --- a/src/Physics/DeepInelastic/XSection/DISXSec.cxx +++ b/src/Physics/DeepInelastic/XSection/DISXSec.cxx @@ -126,6 +126,7 @@ double DISXSec::Integrate( return xsec; } else { + // Just go ahead and integrate the input differential cross section for the // specified interaction. // diff --git a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx index adf1b28eea..14cbb2df6c 100644 --- a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx +++ b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx @@ -356,7 +356,7 @@ void QPMDISStrucFuncBase::Calculate(const Interaction * interaction) const double f = fIncludeNuclMod ? fDISNuclCorr->DISACorrection(interaction) : 1 ; double r = this->R (interaction); // R ~ FL double H = fIncludeH ? this->H(interaction) : 1; - + #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ LOG("DISSF", pDEBUG) << "Nucl. mod = " << f; LOG("DISSF", pDEBUG) << "R(=FL/2xF1) = " << r; diff --git a/src/Physics/NuclearState/BY21DISNuclearModel.cxx b/src/Physics/NuclearState/BY21DISNuclearModel.cxx index 82c30ecf5f..b46e9c8dcd 100644 --- a/src/Physics/NuclearState/BY21DISNuclearModel.cxx +++ b/src/Physics/NuclearState/BY21DISNuclearModel.cxx @@ -126,7 +126,6 @@ void BY21DISNuclearModel::LoadConfig(void) GetParam( "BY21-NuclModel-C-f3", fCf3 ) ; GetParam( "BY21-NuclModel-C-f4", fCf4 ) ; GetParam( "BY21-NuclModel-C-f5", fCf5 ) ; - GetParam( "BY21-NuclModel-C-f6", fCf6 ) ; GetParam( "BY21-NuclModel-Au-f0", fAuf0 ) ; GetParam( "BY21-NuclModel-Au-f1", fAuf1 ) ; diff --git a/src/Physics/NuclearState/JGDISNuclearModel.cxx b/src/Physics/NuclearState/JGDISNuclearModel.cxx index 622b0dbe2b..6babfeacc7 100644 --- a/src/Physics/NuclearState/JGDISNuclearModel.cxx +++ b/src/Physics/NuclearState/JGDISNuclearModel.cxx @@ -97,4 +97,5 @@ void JGDISNuclearModel::LoadConfig(void) GetParam( "JG-NuclModel-a5", fAa5 ) ; GetParam( "JG-NuclModel-a6", fAa6 ) ; GetParam( "JG-NuclModel-a7", fAa7 ) ; + GetParam( "JG-NuclModel-a8", fAa8 ) ; } diff --git a/src/Physics/NuclearState/JGDISNuclearModel.h b/src/Physics/NuclearState/JGDISNuclearModel.h index 69e765704d..c2bf9f3d53 100644 --- a/src/Physics/NuclearState/JGDISNuclearModel.h +++ b/src/Physics/NuclearState/JGDISNuclearModel.h @@ -51,6 +51,7 @@ class JGDISNuclearModel : public DISNuclearModelI { double fAa5; //> JG-NuclModel-a5 double fAa6; //> JG-NuclModel-a6 double fAa7; //> JG-NuclModel-a7 + double fAa8; //> JG-NuclModel-a8 }; } // genie namespace From b23828d8a9602ca838f085cc22dc93afda8fa7fd Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Wed, 25 Feb 2026 16:19:34 +0100 Subject: [PATCH 59/68] changing x range for JG model --- src/Physics/NuclearState/JGDISNuclearModel.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Physics/NuclearState/JGDISNuclearModel.cxx b/src/Physics/NuclearState/JGDISNuclearModel.cxx index 6babfeacc7..885297b783 100644 --- a/src/Physics/NuclearState/JGDISNuclearModel.cxx +++ b/src/Physics/NuclearState/JGDISNuclearModel.cxx @@ -44,7 +44,7 @@ double JGDISNuclearModel::DISACorrection (const Interaction * interaction) const double x = kinematics.x(); int A = tgt.A(); - double xv = TMath::Min(0.75, TMath::Max(0.05, x)); + double xv = TMath::Min(0.65, TMath::Max(0.05, x)); double xv2 = xv * xv; double xv3 = xv2 * xv; double xv4 = xv3 * xv; From a9fe9913241dd0da19ebbd360710a3c18b8869f6 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Thu, 26 Feb 2026 09:11:38 +0100 Subject: [PATCH 60/68] fix potential bug. If free nucleon is asked, scaling factor is lost. It should be applied at the free nucleon level --- src/Physics/DeepInelastic/XSection/QPMDISPXSec.cxx | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/Physics/DeepInelastic/XSection/QPMDISPXSec.cxx b/src/Physics/DeepInelastic/XSection/QPMDISPXSec.cxx index 9585cfaa8a..1d49596f38 100644 --- a/src/Physics/DeepInelastic/XSection/QPMDISPXSec.cxx +++ b/src/Physics/DeepInelastic/XSection/QPMDISPXSec.cxx @@ -147,6 +147,12 @@ double QPMDISPXSec::XSec( double xsec = front_factor * (term1 + term2 + term3 + term4 + term5); xsec = TMath::Max(xsec,0.); + + // Apply scaling / if required to reach well known asymmptotic value + if( proc_info.IsWeakCC() ) xsec *= fCCScale; + else if( proc_info.IsWeakNC() ) xsec *= fEMScale; + else if( proc_info.IsEM() ) xsec *= fEMScale; + #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ LOG("DISPXSec", pINFO) << "d2xsec/dxdy[FreeN] (E= " << E @@ -169,12 +175,7 @@ double QPMDISPXSec::XSec( int nucpdgc = target.HitNucPdg(); int NNucl = (pdg::IsProton(nucpdgc)) ? target.Z() : target.N(); xsec *= NNucl; - - // Apply scaling / if required to reach well known asymmptotic value - if( proc_info.IsWeakCC() ) xsec *= fCCScale; - else if( proc_info.IsWeakNC() ) xsec *= fEMScale; - else if( proc_info.IsEM() ) xsec *= fEMScale; - + // Subtract the inclusive charm production cross section interaction->ExclTagPtr()->SetCharm(); double xsec_charm = fCharmProdModel->XSec(interaction,kps); From 8ef639da8e01bad2f78c72e547c3ef9ef66f218b Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Thu, 26 Feb 2026 09:48:22 +0100 Subject: [PATCH 61/68] update code so it does not apply correction for neutron and proton targets when specifying hydrogen --- .../NuclearState/BY21DISNuclearModel.cxx | 2 +- .../NuclearState/JGDISNuclearModel.cxx | 28 ++++++++++--------- 2 files changed, 16 insertions(+), 14 deletions(-) diff --git a/src/Physics/NuclearState/BY21DISNuclearModel.cxx b/src/Physics/NuclearState/BY21DISNuclearModel.cxx index b46e9c8dcd..200c7e5070 100644 --- a/src/Physics/NuclearState/BY21DISNuclearModel.cxx +++ b/src/Physics/NuclearState/BY21DISNuclearModel.cxx @@ -43,7 +43,7 @@ double BY21DISNuclearModel::DISACorrection (const Interaction * interaction) con const Kinematics & kinematics = interaction->Kine(); double x = kinematics.x(); int A = tgt.A(); - + double xv = TMath::Min(0.75, TMath::Max(0.05, x)); double xv2 = xv * xv; double xv3 = xv2 * xv; diff --git a/src/Physics/NuclearState/JGDISNuclearModel.cxx b/src/Physics/NuclearState/JGDISNuclearModel.cxx index 885297b783..fdb541d338 100644 --- a/src/Physics/NuclearState/JGDISNuclearModel.cxx +++ b/src/Physics/NuclearState/JGDISNuclearModel.cxx @@ -43,20 +43,22 @@ double JGDISNuclearModel::DISACorrection (const Interaction * interaction) const const Kinematics & kinematics = interaction->Kine(); double x = kinematics.x(); int A = tgt.A(); + + if ( A >= 2 ) { + double xv = TMath::Min(0.65, TMath::Max(0.05, x)); + double xv2 = xv * xv; + double xv3 = xv2 * xv; + double xv4 = xv3 * xv; + double xv5 = xv4 * xv; + double xv6 = xv5 * xv; + double xv7 = xv6 * xv; + double xv8 = xv7 * xv; + double alpha = fAa0 + fAa1 * xv + fAa2 * xv2 + fAa3 * xv3 + fAa4 * xv4 + fAa5 * xv5 + fAa6 * xv6 + fAa7 * xv7 + fAa8 * xv8; + double lnC = fAc0 + fAc1 * TMath::Log(xv) + fAc2 * pow( TMath::Log(xv),2 ); + double C = TMath::Exp( lnC ) ; + f *= C * pow( A, alpha ) ; + } - double xv = TMath::Min(0.65, TMath::Max(0.05, x)); - double xv2 = xv * xv; - double xv3 = xv2 * xv; - double xv4 = xv3 * xv; - double xv5 = xv4 * xv; - double xv6 = xv5 * xv; - double xv7 = xv6 * xv; - double xv8 = xv7 * xv; - double alpha = fAa0 + fAa1 * xv + fAa2 * xv2 + fAa3 * xv3 + fAa4 * xv4 + fAa5 * xv5 + fAa6 * xv6 + fAa7 * xv7 + fAa8 * xv8; - double lnC = fAc0 + fAc1 * TMath::Log(xv) + fAc2 * pow( TMath::Log(xv),2 ); - double C = TMath::Exp( lnC ) ; - f *= C * pow( A, alpha ) ; - #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ LOG("DISSF", pDEBUG) << "J. Gomez et. al. Nuclear factor for x of " << x << " = " << f; #endif From 5c93a8078e4a124937317f541cc26682d4976553 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Thu, 26 Feb 2026 09:57:32 +0100 Subject: [PATCH 62/68] updating xml file --- config/BYStrucFunc2021.xml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/config/BYStrucFunc2021.xml b/config/BYStrucFunc2021.xml index e17448194a..c02a142fc4 100644 --- a/config/BYStrucFunc2021.xml +++ b/config/BYStrucFunc2021.xml @@ -85,12 +85,12 @@ WeinbergAngle double No genie::JGDISNuclearModel/Default - - genie::BYDISNuclearModel21/Default + + genie::BY21DISNuclearModel/Default - - genie::BYDISNuclearModel/Default + + genie::BY00DISNuclearModel/Default From 46b09fe0819ceec5b8046d6ac4f403f0a6a51dd3 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Thu, 26 Feb 2026 10:00:24 +0100 Subject: [PATCH 63/68] updating name of class --- src/Physics/NuclearState/LinkDef.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Physics/NuclearState/LinkDef.h b/src/Physics/NuclearState/LinkDef.h index 742f427860..c4a74ac043 100644 --- a/src/Physics/NuclearState/LinkDef.h +++ b/src/Physics/NuclearState/LinkDef.h @@ -25,8 +25,8 @@ #pragma link C++ class genie::SRCNuclearRecoil; #pragma link C++ class genie::SecondNucleonEmissionI; #pragma link C++ class genie::SpectralFunction2p2h; -#pragma link C++ class BYDISNuclearModel; -#pragma link C++ class BYDISNuclearModel21; +#pragma link C++ class BY00DISNuclearModel; +#pragma link C++ class BY21DISNuclearModel; #pragma link C++ class JGDISNuclearModel; #endif From 98fb06271e20b760e6437df3ad26e6a49d02d477 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Thu, 26 Feb 2026 10:09:35 +0100 Subject: [PATCH 64/68] updating config in algorithm pool --- config/master_config.xml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/config/master_config.xml b/config/master_config.xml index 72d1486b24..2cd3a1e9dc 100644 --- a/config/master_config.xml +++ b/config/master_config.xml @@ -247,8 +247,8 @@ SpectralFunc1d.xml SpectralFunc.xml SmithMonizUtils.xml - BYDISNuclearModel.xml - BYDISNuclearModel21.xml + BY00DISNuclearModel.xml + BY21DISNuclearModel.xml JGDISNuclearModel.xml From e71a1323bf25e6f28d8f1e2c59d81e22baa4fffb Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Thu, 26 Feb 2026 10:40:02 +0100 Subject: [PATCH 65/68] commenting out second int. method as it takes a very long time. Computing from first principles is faster --- .../DeepInelastic/XSection/DISXSec.cxx | 252 +++++++++--------- 1 file changed, 127 insertions(+), 125 deletions(-) diff --git a/src/Physics/DeepInelastic/XSection/DISXSec.cxx b/src/Physics/DeepInelastic/XSection/DISXSec.cxx index a65d2bcaf5..90babd0277 100644 --- a/src/Physics/DeepInelastic/XSection/DISXSec.cxx +++ b/src/Physics/DeepInelastic/XSection/DISXSec.cxx @@ -1,10 +1,10 @@ //____________________________________________________________________________ /* - Copyright (c) 2003-2025, The GENIE Collaboration - For the full text of the license visit http://copyright.genie-mc.org + Copyright (c) 2003-2025, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org - Costas Andreopoulos - University of Liverpool + Costas Andreopoulos + University of Liverpool */ //____________________________________________________________________________ @@ -37,13 +37,13 @@ using namespace genie::constants; //____________________________________________________________________________ DISXSec::DISXSec() : -XSecIntegratorI("genie::DISXSec") + XSecIntegratorI("genie::DISXSec") { } //____________________________________________________________________________ DISXSec::DISXSec(string config) : -XSecIntegratorI("genie::DISXSec", config) + XSecIntegratorI("genie::DISXSec", config) { } @@ -54,14 +54,14 @@ DISXSec::~DISXSec() } //____________________________________________________________________________ double DISXSec::Integrate( - const XSecAlgorithmI * model, const Interaction * in) const + const XSecAlgorithmI * model, const Interaction * in) const { if(! model->ValidProcess(in) ) return 0.; const KPhaseSpace & kps = in->PhaseSpace(); if(!kps.IsAboveThreshold()) { - LOG("DISXSec", pDEBUG) << "*** Below energy threshold"; - return 0; + LOG("DISXSec", pDEBUG) << "*** Below energy threshold"; + return 0; } const InitialState & init_state = in->InitState(); @@ -69,7 +69,7 @@ double DISXSec::Integrate( int nucpdgc = init_state.Tgt().HitNucPdg(); int NNucl = (pdg::IsProton(nucpdgc)) ? - init_state.Tgt().Z() : init_state.Tgt().N(); + init_state.Tgt().Z() : init_state.Tgt().N(); // If the input interaction is off a nuclear target, then chek whether // the corresponding free nucleon cross section already exists at the @@ -77,7 +77,7 @@ double DISXSec::Integrate( // If yes, calculate the nuclear cross section based on that value. // XSecSplineList * xsl = XSecSplineList::Instance(); - if(init_state.Tgt().IsNucleus() && !xsl->IsEmpty() && !fDISNuclCorr ) { + if( init_state.Tgt().IsNucleus() && !xsl->IsEmpty() && !fDISNuclCorr ) { // Computes xsec from free nucleon calculation wo nuclear effects Interaction * interaction = new Interaction(*in); Target * target = interaction->InitStatePtr()->TgtPtr(); @@ -89,8 +89,8 @@ double DISXSec::Integrate( LOG("DISXSec", pINFO) << "From XSecSplineList: XSec[DIS,free nucleon] (E = " << Ev << " GeV) = " << xsec; if(! interaction->TestBit(kIAssumeFreeNucleon) ) { - xsec *= NNucl; - LOG("DISXSec", pINFO) << "XSec[DIS] (E = " << Ev << " GeV) = " << xsec; + xsec *= NNucl; + LOG("DISXSec", pINFO) << "XSec[DIS] (E = " << Ev << " GeV) = " << xsec; } delete interaction; return xsec; @@ -98,86 +98,88 @@ double DISXSec::Integrate( delete interaction; } + ///// THIS TAKES FAR TOO LONG! Removed. // There was no corresponding free nucleon spline saved in XSecSplineList that // could be used to speed up this calculation. // Check whether local caching of free nucleon cross sections is allowed. // If yes, store free nucleon cross sections at a cache branch and use those // at any subsequent call. + // If DIS nuclear effects are used, precompute from scratch. Cannot reuse free nucleon spline. + /* + bool precalc_bare_xsec = RunOpt::Instance()->BareXSecPreCalc(); + if(precalc_bare_xsec && !fDISNuclCorr) { + std::cout << " ******** IN INT METHOD 2 " << std::endl; + Cache * cache = Cache::Instance(); + Interaction * interaction = new Interaction(*in); + string key = this->CacheBranchName(model,interaction); + LOG("DISXSec", pINFO) << "Finding cache branch with key: " << key; + CacheBranchFx * cache_branch = + dynamic_cast (cache->FindCacheBranch(key)); + if(!cache_branch) { + this->CacheFreeNucleonXSec(model,interaction); + cache_branch = + dynamic_cast (cache->FindCacheBranch(key)); + assert(cache_branch); + } + const CacheBranchFx & cb = (*cache_branch); + double xsec = cb(Ev); + if(! interaction->TestBit(kIAssumeFreeNucleon) ) { xsec *= NNucl; } + LOG("DISXSec", pINFO) << "XSec[DIS] (E = " << Ev << " GeV) = " << xsec; + delete interaction; + return xsec; + } + */ + + // Just go ahead and integrate the input differential cross section for the + // specified interaction. // - bool precalc_bare_xsec = RunOpt::Instance()->BareXSecPreCalc(); - if(precalc_bare_xsec && !fDISNuclCorr) { ///// CHECK LATER. I THINK NOT NEEDED THE CHANGE HERE - Cache * cache = Cache::Instance(); - Interaction * interaction = new Interaction(*in); - string key = this->CacheBranchName(model,interaction); - LOG("DISXSec", pINFO) << "Finding cache branch with key: " << key; - CacheBranchFx * cache_branch = - dynamic_cast (cache->FindCacheBranch(key)); - if(!cache_branch) { - this->CacheFreeNucleonXSec(model,interaction); - cache_branch = - dynamic_cast (cache->FindCacheBranch(key)); - assert(cache_branch); - } - const CacheBranchFx & cb = (*cache_branch); - double xsec = cb(Ev); - if(! interaction->TestBit(kIAssumeFreeNucleon) ) { xsec *= NNucl; } - LOG("DISXSec", pINFO) << "XSec[DIS] (E = " << Ev << " GeV) = " << xsec; - delete interaction; - return xsec; - } - else { - - // Just go ahead and integrate the input differential cross section for the - // specified interaction. - // - Interaction * interaction = new Interaction(*in); - interaction->SetBit(kISkipProcessChk); -// interaction->SetBit(kISkipKinematicChk); - - // **Important note** - // Based on discussions with Hugh at the GENIE mini-workshop / RAL - July '07 - // The DIS nuclear corrections re-distribute the strength in x,y but do not - // affect the total cross-section They should be disabled at this step. - // But they should be enabled at the DIS thread's kinematical selection. - // Since nuclear corrections don't need to be included at this stage, all the - // nuclear cross sections can be trivially built from the free nucleon ones. - // - if(!fDISNuclCorr) interaction->SetBit(kINoNuclearCorrection); - - Range1D_t Wl = kps.WLim(); - Range1D_t Q2l = kps.Q2Lim(); - LOG("DISXSec", pINFO) - << "W integration range = [" << Wl.min << ", " << Wl.max << "]"; - LOG("DISXSec", pINFO) - << "Q2 integration range = [" << Q2l.min << ", " << Q2l.max << "]"; - - bool phsp_ok = - (Q2l.min >= 0. && Q2l.max >= 0. && Q2l.max >= Q2l.min && - Wl.min >= 0. && Wl.max >= 0. && Wl.max >= Wl.min); - - double xsec = 0.; - - if(phsp_ok) { - ROOT::Math::IBaseFunctionMultiDim * func = - new utils::gsl::d2XSec_dWdQ2_E(model, interaction); - ROOT::Math::IntegrationMultiDim::Type ig_type = - utils::gsl::IntegrationNDimTypeFromString(fGSLIntgType); - - double abstol = 1; //We mostly care about relative tolerance. - ROOT::Math::IntegratorMultiDim ig(*func, ig_type, abstol, fGSLRelTol, fGSLMaxEval); - double kine_min[2] = { Wl.min, Q2l.min }; - double kine_max[2] = { Wl.max, Q2l.max }; - xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2); - delete func; - }//phase space ok? - - LOG("DISXSec", pINFO) << "XSec[DIS] (E = " << Ev << " GeV) = " << xsec; - - delete interaction; - - return xsec; - } - return 0; + Interaction * interaction = new Interaction(*in); + interaction->SetBit(kISkipProcessChk); + // interaction->SetBit(kISkipKinematicChk); + + // **Important note** + // Based on discussions with Hugh at the GENIE mini-workshop / RAL - July '07 + // The DIS nuclear corrections re-distribute the strength in x,y but do not + // affect the total cross-section They should be disabled at this step. + // But they should be enabled at the DIS thread's kinematical selection. + // Since nuclear corrections don't need to be included at this stage, all the + // nuclear cross sections can be trivially built from the free nucleon ones. + // + if(!fDISNuclCorr) interaction->SetBit(kINoNuclearCorrection); + + Range1D_t Wl = kps.WLim(); + Range1D_t Q2l = kps.Q2Lim(); + LOG("DISXSec", pINFO) + << "W integration range = [" << Wl.min << ", " << Wl.max << "]"; + LOG("DISXSec", pINFO) + << "Q2 integration range = [" << Q2l.min << ", " << Q2l.max << "]"; + + bool phsp_ok = + (Q2l.min >= 0. && Q2l.max >= 0. && Q2l.max >= Q2l.min && + Wl.min >= 0. && Wl.max >= 0. && Wl.max >= Wl.min); + + double xsec = 0.; + + if(phsp_ok) { + ROOT::Math::IBaseFunctionMultiDim * func = + new utils::gsl::d2XSec_dWdQ2_E(model, interaction); + ROOT::Math::IntegrationMultiDim::Type ig_type = + utils::gsl::IntegrationNDimTypeFromString(fGSLIntgType); + + double abstol = 1; //We mostly care about relative tolerance. + ROOT::Math::IntegratorMultiDim ig(*func, ig_type, abstol, fGSLRelTol, fGSLMaxEval); + double kine_min[2] = { Wl.min, Q2l.min }; + double kine_max[2] = { Wl.max, Q2l.max }; + xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2); + delete func; + }//phase space ok? + + LOG("DISXSec", pINFO) << "XSec[DIS] (E = " << Ev << " GeV) = " << xsec; + + delete interaction; + + return xsec; + } //____________________________________________________________________________ void DISXSec::Configure(const Registry & config) @@ -216,16 +218,16 @@ void DISXSec::LoadConfig(void) } //____________________________________________________________________________ void DISXSec::CacheFreeNucleonXSec( - const XSecAlgorithmI * model, const Interaction * interaction) const + const XSecAlgorithmI * model, const Interaction * interaction) const { LOG("DISXSec", pWARN) - << "Wait while computing/caching free nucleon DIS xsections first..."; + << "Wait while computing/caching free nucleon DIS xsections first..."; // Create the cache branch Cache * cache = Cache::Instance(); string key = this->CacheBranchName(model,interaction); CacheBranchFx * cache_branch = - dynamic_cast (cache->FindCacheBranch(key)); + dynamic_cast (cache->FindCacheBranch(key)); assert(!cache_branch); cache_branch = new CacheBranchFx("DIS XSec"); cache->AddCacheBranch(key, cache_branch); @@ -256,18 +258,18 @@ void DISXSec::CacheFreeNucleonXSec( // knots < energy threshold double dEb = (Ethr>Emin) ? (Ethr - Emin) / nkb : 0; for(int i=0; i= energy threshold double E0 = TMath::Max(Ethr,Emin); double dEa = (TMath::Log10(Emax) - TMath::Log10(E0)) /(nka-1); for(int i=0; iInitStatePtr()->SetProbeP4(p4); double xsec = 0.; if(Ev>Ethr+kASmallNum) { - Range1D_t Wl = kps.WLim(); - Range1D_t Q2l = kps.Q2Lim(); - LOG("DISXSec", pINFO) - << "W integration range = [" << Wl.min << ", " << Wl.max << "]"; - LOG("DISXSec", pINFO) - << "Q2 integration range = [" << Q2l.min << ", " << Q2l.max << "]"; - - bool phsp_ok = - (Q2l.min >= 0. && Q2l.max >= 0. && Q2l.max >= Q2l.min && - Wl.min >= 0. && Wl.max >= 0. && Wl.max >= Wl.min); - - if(phsp_ok) { - ROOT::Math::IntegrationMultiDim::Type ig_type = - utils::gsl::IntegrationNDimTypeFromString(fGSLIntgType); - double abstol = 1; //We mostly care about relative tolerance. - ROOT::Math::IntegratorMultiDim ig(*func, ig_type, abstol, fGSLRelTol, fGSLMaxEval); - - if (ig_type == ROOT::Math::IntegrationMultiDim::kADAPTIVE) { - ROOT::Math::AdaptiveIntegratorMultiDim * cast = - dynamic_cast( ig.GetIntegrator() ); - assert(cast); - cast->SetMinPts(fGSLMinEval); - } - - double kine_min[2] = { Wl.min, Q2l.min }; - double kine_max[2] = { Wl.max, Q2l.max }; - xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2); - }// phase space limits ok? + Range1D_t Wl = kps.WLim(); + Range1D_t Q2l = kps.Q2Lim(); + LOG("DISXSec", pINFO) + << "W integration range = [" << Wl.min << ", " << Wl.max << "]"; + LOG("DISXSec", pINFO) + << "Q2 integration range = [" << Q2l.min << ", " << Q2l.max << "]"; + + bool phsp_ok = + (Q2l.min >= 0. && Q2l.max >= 0. && Q2l.max >= Q2l.min && + Wl.min >= 0. && Wl.max >= 0. && Wl.max >= Wl.min); + + if(phsp_ok) { + ROOT::Math::IntegrationMultiDim::Type ig_type = + utils::gsl::IntegrationNDimTypeFromString(fGSLIntgType); + double abstol = 1; //We mostly care about relative tolerance. + ROOT::Math::IntegratorMultiDim ig(*func, ig_type, abstol, fGSLRelTol, fGSLMaxEval); + + if (ig_type == ROOT::Math::IntegrationMultiDim::kADAPTIVE) { + ROOT::Math::AdaptiveIntegratorMultiDim * cast = + dynamic_cast( ig.GetIntegrator() ); + assert(cast); + cast->SetMinPts(fGSLMinEval); + } + + double kine_min[2] = { Wl.min, Q2l.min }; + double kine_max[2] = { Wl.max, Q2l.max }; + xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2); + }// phase space limits ok? }//Ev>threshold LOG("DISXSec", pNOTICE) - << "Caching: XSec[DIS] (E = " << Ev << " GeV) = " - << xsec / (1E-38 * units::cm2) << " x 1E-38 cm^2"; + << "Caching: XSec[DIS] (E = " << Ev << " GeV) = " + << xsec / (1E-38 * units::cm2) << " x 1E-38 cm^2"; cache_branch->AddValues(Ev,xsec); }//ie @@ -320,9 +322,9 @@ void DISXSec::CacheFreeNucleonXSec( } //____________________________________________________________________________ string DISXSec::CacheBranchName( - const XSecAlgorithmI * model, const Interaction * interaction) const + const XSecAlgorithmI * model, const Interaction * interaction) const { -// Build a unique name for the cache branch + // Build a unique name for the cache branch Cache * cache = Cache::Instance(); From 50f729f043aa5cef9253ac9ee4a9a68b9a5f24f2 Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Thu, 26 Feb 2026 13:00:23 +0100 Subject: [PATCH 66/68] adding minimum No points --- src/Physics/DeepInelastic/XSection/DISXSec.cxx | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/Physics/DeepInelastic/XSection/DISXSec.cxx b/src/Physics/DeepInelastic/XSection/DISXSec.cxx index 90babd0277..89cf1597e6 100644 --- a/src/Physics/DeepInelastic/XSection/DISXSec.cxx +++ b/src/Physics/DeepInelastic/XSection/DISXSec.cxx @@ -168,6 +168,14 @@ double DISXSec::Integrate( double abstol = 1; //We mostly care about relative tolerance. ROOT::Math::IntegratorMultiDim ig(*func, ig_type, abstol, fGSLRelTol, fGSLMaxEval); + + if (ig_type == ROOT::Math::IntegrationMultiDim::kADAPTIVE) { + ROOT::Math::AdaptiveIntegratorMultiDim * cast = + dynamic_cast( ig.GetIntegrator() ); + assert(cast); + cast->SetMinPts(fGSLMinEval); + } + double kine_min[2] = { Wl.min, Q2l.min }; double kine_max[2] = { Wl.max, Q2l.max }; xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2); From 02933478b45d09abe23fc8580e34ce7cc605032e Mon Sep 17 00:00:00 2001 From: jtenavidal Date: Mon, 2 Mar 2026 15:10:31 +0100 Subject: [PATCH 67/68] splitting functions for KAxial and KVector Factors --- config/DISHadronicSystemGenerator.xml | 2 +- config/UnstableParticleDecayer.xml | 6 +++--- .../DeepInelastic/XSection/BYStrucFunc.cxx | 15 ++++++++++--- .../DeepInelastic/XSection/BYStrucFunc.h | 5 ++--- .../XSection/BYStrucFunc2021.cxx | 19 ++++++++++++----- .../DeepInelastic/XSection/BYStrucFunc2021.h | 6 +++--- .../XSection/QPMDISStrucFuncBase.cxx | 21 ++++++++++++++++--- .../XSection/QPMDISStrucFuncBase.h | 7 ++++--- 8 files changed, 57 insertions(+), 24 deletions(-) diff --git a/config/DISHadronicSystemGenerator.xml b/config/DISHadronicSystemGenerator.xml index ce78b8229b..925e2f915b 100644 --- a/config/DISHadronicSystemGenerator.xml +++ b/config/DISHadronicSystemGenerator.xml @@ -40,7 +40,7 @@ NUCL-NR double No - genie::AGCharmPythia6Hadro2023/Default + genie::AGCharmPythia8Hadro2023/Default genie::UnstableParticleDecayer/BeforeHadronTransport false diff --git a/config/UnstableParticleDecayer.xml b/config/UnstableParticleDecayer.xml index 49af36c185..08237bab1f 100644 --- a/config/UnstableParticleDecayer.xml +++ b/config/UnstableParticleDecayer.xml @@ -13,19 +13,19 @@ Name Type Optional Comment Default 2 - genie::Pythia6Decayer2023/BeforeHadronTransport + genie::Pythia8Decayer2023/BeforeHadronTransport genie::BaryonResonanceDecayer/BeforeHadronTransport 2 - genie::Pythia6Decayer2023/AfterHadronTransport + genie::Pythia8Decayer2023/AfterHadronTransport genie::BaryonResonanceDecayer/AfterHadronTransport 3 - genie::Pythia6Decayer2023/AfterHadronTransport + genie::Pythia8Decayer2023/AfterHadronTransport genie::BaryonResonanceDecayer/AfterHadronTransport genie::DarkSectorDecayer/Default diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc.cxx index e5f370da92..5480ae3e68 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc.cxx @@ -99,12 +99,14 @@ double BYStrucFunc::ScalingVar(const Interaction * interaction, double Mf ) cons return xw; } //____________________________________________________________________________ -void BYStrucFunc::KFactors(const Interaction * interaction, +void BYStrucFunc::KVectorFactors(const Interaction * interaction, double & kuv, double & kdv, double & kus, double & kds, double & kss ) const { -// Overrides QPMDISStrucFuncBase::KFactors() to compute the BY K factors for + +// Overrides QPMDISStrucFuncBase::KVectorFactors() to compute the BY K factors for // u(valence), d(valence), u(sea), d(sea); - +// In this version, Vector K-Factors are used for the axial component + double myQ2 = this->Q2(interaction); double GD = 1. / TMath::Power(1.+myQ2/0.71, 2); // p elastic form factor double GD2 = TMath::Power(GD,2); @@ -115,3 +117,10 @@ void BYStrucFunc::KFactors(const Interaction * interaction, kds = myQ2/(myQ2+fCsD); // K - d(sea) } //____________________________________________________________________________ +void BYStrucFunc::KAxialFactors(const Interaction * interaction, + double & kuv, double & kdv, double & kus, double & kds, double & kss ) const { + // In this version, Vector K-Factors are used for the axial component + KVectorFactors( interaction, kuv, kdv, kus, kds, kss ); +} + +//____________________________________________________________________________ diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc.h b/src/Physics/DeepInelastic/XSection/BYStrucFunc.h index 8e00d64d22..7cf6e27a7f 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc.h +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc.h @@ -46,9 +46,8 @@ class BYStrucFunc : public QPMDISStrucFuncBase { // override part of the DISStructureFuncModel implementation // to compute all the corrections applied by the Bodek-Yang model. double ScalingVar (const Interaction * i, double Mf = 0 ) const; - void KFactors (const Interaction * i, double & kuv, - double & kdv, double & kus, double & kds, double & kss ) const; - + void KVectorFactors (const Interaction * i, double & kuv, double & kdv, double & kus, double & kds, double & kss ) const; + void KAxialFactors (const Interaction * i, double & kuv, double & kdv, double & kus, double & kds, double & kss ) const ; // Bodek-Yang model-specific parameters double fA; ///< better scaling var parameter A diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx index 4b7bd595bc..caa0b3073e 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx @@ -132,8 +132,8 @@ double BYStrucFunc2021::ScalingVar(const Interaction * interaction, double Mf ) return xw; } //____________________________________________________________________________ -void BYStrucFunc2021::KFactors(const Interaction * interaction, - double & kuv, double & kdv, double & kus, double & kds, double & kss ) const +void BYStrucFunc2021::KVectorFactors(const Interaction * interaction, + double & kuv, double & kdv, double & kus, double & kds, double & kss ) const { // Overrides QPMDISStrucFuncBase::KFactors() to compute the BY K factors for // u(valence), d(valence), u(sea), d(sea), s(sea); @@ -159,11 +159,13 @@ void BYStrucFunc2021::KFactors(const Interaction * interaction, kss = myQ2/(myQ2+fCsS); // K - s(sea) } //____________________________________________________________________________ -void BYStrucFunc2021::KAxialFactors (const Interaction * interaction, double & ksea, double & kvalance ) const { +void BYStrucFunc2021::KAxialFactors(const Interaction * interaction, + double & kuv, double & kdv, double & kus, double & kds, double & kss ) const { + // https://arxiv.org/pdf/2108.09240 Sec 11.2 double myQ2 = this->Q2(interaction); - ksea = ( myQ2 + fPsA*fCsA ) / ( myQ2 + fCsA ) ; - kvalance = ( myQ2 + fPvA * 0.18 ) / ( myQ2 + 0.18 ) ; + double ksea = ( myQ2 + fPsA*fCsA ) / ( myQ2 + fCsA ) ; + double kvalance = ( myQ2 + fPvA * 0.18 ) / ( myQ2 + 0.18 ) ; // Compute Low-q0 modificaion factor for neutrinos and anti-neutrinos double q0 = this->q0(interaction); @@ -188,6 +190,13 @@ void BYStrucFunc2021::KAxialFactors (const Interaction * interaction, double & k // double check it only affects valance factors: //kvalance *= KLW ; */ + // The 2021 BY model uses the same factors for up and down valance quarks, and for u, d, s sea quarks. + kuv = kvalance; + kdv = kvalance; + kus = ksea; + kds = ksea; + kss = ksea; + } //____________________________________________________________________________ diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h index 3714b13b8b..21137002b2 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h @@ -1,4 +1,3 @@ -//____________________________________________________________________________ /*! \class genie::BYStrucFunc2021 @@ -49,8 +48,9 @@ class BYStrucFunc2021 : public QPMDISStrucFuncBase { // override part of the DISStructureFuncModel implementation // to compute all the corrections applied by the Bodek-Yang model. double ScalingVar (const Interaction * i, double Mf = 0 ) const; - void KFactors (const Interaction * i, double & kuv, double & kdv, double & kus, double & kds, double & kss ) const; - void KAxialFactors (const Interaction * i, double & ksea, double & kvalance ) const; + void KVectorFactors (const Interaction * i, double & kuv, double & kdv, double & kus, double & kds, double & kss ) const; + void KAxialFactors(const Interaction * i, double & kuv, double & kdv, double & kus, double & kds, double & kss ) const ; + double R(const Interaction * interaction) const ; // overrides QPMDISStrucFuncBase implementation double H(const Interaction * interaction) const ; // overrides QPMDISStrucFuncBase implementation double NuclMod(const Interaction * interaction) const; diff --git a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx index 14cbb2df6c..e1a18e59d0 100644 --- a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx +++ b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx @@ -460,8 +460,23 @@ double QPMDISStrucFuncBase::ScalingVar(const Interaction* interaction, double Mf return interaction->Kine().x(); } //____________________________________________________________________________ -void QPMDISStrucFuncBase::KFactors(const Interaction *, - double & kuv, double & kdv, double & kus, double & kds, double & kss ) const +void QPMDISStrucFuncBase::KVectorFactors(const Interaction *, + double & kuv, double & kdv, double & kus, double & kds, double & kss ) const +{ + // This is an abstract class: no model-specific correction + // The PDF scaling variables are set to 1 + // Override this method to compute model-dependent corrections + + kuv = 1.; + kdv = 1.; + kus = 1.; + kds = 1.; + kss = 1.; +} + +//____________________________________________________________________________ +void QPMDISStrucFuncBase::KAxialFactors(const Interaction *, + double & kuv, double & kdv, double & kus, double & kds, double & kss ) const { // This is an abstract class: no model-specific correction // The PDF scaling variables are set to 1 @@ -559,7 +574,7 @@ void QPMDISStrucFuncBase::CalcPDFs(const Interaction * interaction) const double ksea_d = 1.; double ksea_s = 1.; - this->KFactors(interaction, kval_u, kval_d, ksea_u, ksea_d, ksea_s); + this->KVectorFactors(interaction, kval_u, kval_d, ksea_u, ksea_d, ksea_s); #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ LOG("DISSF", pDEBUG) << "K-Factors:"; diff --git a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h index 1a3384b988..fc49f52139 100644 --- a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h +++ b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h @@ -1,4 +1,3 @@ -//____________________________________________________________________________ /*! \class genie::QPMDISStrucFuncBase @@ -69,8 +68,10 @@ class QPMDISStrucFuncBase : public DISStructureFuncModelI { virtual void CalcPDFs (const Interaction * i) const; virtual double R (const Interaction * i) const; virtual double H (const Interaction * i) const; - virtual void KFactors (const Interaction * i, double & kuv, - double & kdv, double & kus, double & kds, double & kss ) const; + virtual void KVectorFactors (const Interaction * i, double & kuv, + double & kdv, double & kus, double & kds, double & kss ) const; + virtual void KAxialFactors (const Interaction * i, double & kuv, + double & kdv, double & kus, double & kds, double & kss ) const; // configuration // double fQ2min; ///< min Q^2 allowed for PDFs: PDF(Q2 Date: Fri, 13 Mar 2026 11:52:21 +0100 Subject: [PATCH 68/68] splitting axial and vector contributions for CC F2 and CCF3 only --- .../XSection/QPMDISStrucFuncBase.cxx | 39 ++++++++++++++----- 1 file changed, 29 insertions(+), 10 deletions(-) diff --git a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx index e1a18e59d0..76020ee7cf 100644 --- a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx +++ b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx @@ -323,8 +323,11 @@ void QPMDISStrucFuncBase::Calculate(const Interaction * interaction) const return; } - F2val = 2*(q+qbar); - xF3val = 2*(q-qbar); + F2val = (q+qbar); + xF3val = (q-qbar); + // Removing factor 2, at low Q2 we have the split between axial and vector + // F2val = 2*(q+qbar); + // xF3val = 2*(q-qbar); } // *** ELECTROMAGNETIC @@ -567,14 +570,30 @@ void QPMDISStrucFuncBase::CalcPDFs(const Interaction * interaction) const << "The event is below the charm threshold (mcharm = " << fMc << ")"; } - // Compute the K factors - double kval_u = 1.; - double kval_d = 1.; - double ksea_u = 1.; - double ksea_d = 1.; - double ksea_s = 1.; - - this->KVectorFactors(interaction, kval_u, kval_d, ksea_u, ksea_d, ksea_s); + // for vector K Factors + double kval_u_V = 1.; + double kval_d_V = 1.; + double ksea_u_V = 1.; + double ksea_d_V = 1.; + double ksea_s_V = 1.; + + this->KVectorFactors(interaction, kval_u_V, kval_d_V, ksea_u_V, ksea_d_V, ksea_s_V); + + // For axial K factors + double kval_u_A = 1.; + double kval_d_A = 1.; + double ksea_u_A = 1.; + double ksea_d_A = 1.; + double ksea_s_A = 1.; + + this->KAxialFactors(interaction, kval_u_A, kval_d_A, ksea_u_A, ksea_d_A, ksea_s_A); + + // Compute the final (V+A) K factors + double kval_u = kval_u_V + kval_u_A ; + double kval_d = kval_d_V + kval_d_A ; + double ksea_u = ksea_u_V + ksea_u_A ; + double ksea_d = ksea_d_V + ksea_d_A ; + double ksea_s = ksea_s_V + ksea_s_A ; #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ LOG("DISSF", pDEBUG) << "K-Factors:";