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/BYPDF.xml b/config/BYPDF.xml index 72ca650251..7d61b01732 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,9 +28,7 @@ Uncorr-PDF-Set alg No Uncorrected PDF model - - genie::LHAPDF5/GRVLO diff --git a/config/BYStrucFunc.xml b/config/BYStrucFunc.xml index 3d379db2f3..51ed66ca70 100644 --- a/config/BYStrucFunc.xml +++ b/config/BYStrucFunc.xml @@ -39,30 +39,32 @@ 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::BYDISNuclearModel/Default diff --git a/config/BYStrucFunc2021.xml b/config/BYStrucFunc2021.xml new file mode 100644 index 0000000000..c02a142fc4 --- /dev/null +++ b/config/BYStrucFunc2021.xml @@ -0,0 +1,96 @@ + + + + + + + + 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::JGDISNuclearModel/Default + + + + genie::BY21DISNuclearModel/Default + + + + genie::BY00DISNuclearModel/Default + + + 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 - @@ -15,5 +25,10 @@ Configuration for the DISXSec cross section algorithm 0.0001 + + 50000000 + 0.00001 + true + diff --git a/config/JGDISNuclearModel.xml b/config/JGDISNuclearModel.xml new file mode 100644 index 0000000000..8a622d7c9b --- /dev/null +++ b/config/JGDISNuclearModel.xml @@ -0,0 +1,45 @@ + + + + + + + + 0.017 + 0.018 + 0.005 + + -0.070 + 2.189 + -24.667 + 145.291 + -497.237 + 1013.129 + -1208.393 + 775.767 + -205.872 + + + + + diff --git a/config/KNOTunedQPMDISPXSec.xml b/config/KNOTunedQPMDISPXSec.xml index 77cc15ad96..25b07b4adb 100644 --- a/config/KNOTunedQPMDISPXSec.xml +++ b/config/KNOTunedQPMDISPXSec.xml @@ -17,12 +17,20 @@ Wcut double no Co - NonResBackground - genie::QPMDISPXSec/Default genie::AGKYLowW2019/Default 1. + genie::DISXSec/Default + false + + + + NonResBackground + genie::DISXSec/GRV982010 + genie::AGKYLowW2019/Default + 1. + genie::QPMDISPXSec/GRV982010 diff --git a/config/QPMDISPXSec.xml b/config/QPMDISPXSec.xml index 265805a20c..975f685d00 100644 --- a/config/QPMDISPXSec.xml +++ b/config/QPMDISPXSec.xml @@ -20,10 +20,18 @@ WeinbergAngle double No genie::BYStrucFunc/Default genie::DISXSec/Default - 1.032 1.032 1. + + WeakInt + genie::DISXSec/GRV982010 + genie::BYStrucFunc2021/Default + 1. + 1. + 1. + + diff --git a/config/QPMDISStrucFunc.xml b/config/QPMDISStrucFunc.xml index 3c25bbde2c..da916fbacf 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,8 +29,7 @@ WeinbergAngle double No CKM,Masses,WeakInt - genie::GRV98LO/Default - + genie::GRV98LO/Default BetheBlochModel.xml 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/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 }; diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc.cxx b/src/Physics/DeepInelastic/XSection/BYStrucFunc.cxx index 01f5caacd8..5480ae3e68 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 @@ -99,12 +99,14 @@ double BYStrucFunc::ScalingVar(const Interaction * interaction) const return xw; } //____________________________________________________________________________ -void BYStrucFunc::KFactors(const Interaction * interaction, - double & kuv, double & kdv, double & kus, double & kds) const +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 3df0c474ef..7cf6e27a7f 100644 --- a/src/Physics/DeepInelastic/XSection/BYStrucFunc.h +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc.h @@ -45,10 +45,9 @@ 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; - void KFactors (const Interaction * i, double & kuv, - double & kdv, double & kus, double & kds) const; - + double ScalingVar (const Interaction * i, double Mf = 0 ) 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 new file mode 100644 index 0000000000..caa0b3073e --- /dev/null +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.cxx @@ -0,0 +1,272 @@ +//____________________________________________________________________________ +/* + 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" +#include "Framework/ParticleData/PDGUtils.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) +{ + // vector mass + 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. + // + GetParam( "BY-A", fA ) ; + GetParam( "BY-B", fB ) ; + GetParam( "BY-CsU", fCsU ) ; + GetParam( "BY-CsD", fCsD ) ; + GetParam( "BY-CvLW", fCvLW ) ; + GetParam( "BY-Cv1U", fCv1U ) ; + GetParam( "BY-Cv2U", fCv2U ) ; + 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-CaLW-nubar" , fCaLW_nu ) ; + GetParam( "BY-CaLW-nubar" , fCaLW_nubar ) ; + GetParam( "BY-H0" , fH0 ) ; + GetParam( "BY-H1" , fH1 ) ; + GetParam( "BY-H2" , fH2 ) ; + GetParam( "BY-H3" , fH3 ) ; +} +//____________________________________________________________________________ +void BYStrucFunc2021::Init(void) +{ + fMv = 0; + fMv2 = 0; + fA = 0; + fB = 0; + fCsU = 0; + fCsD = 0; + fCvLW = 0; + fCv1U = 0; + fCv2U = 0; + fCv1D = 0; + fCv2D = 0; + fCsS = 0; + fPsA = 0; + fPvA = 0; + fCsA = 0; + fCaLW_nu = 0; + fCaLW_nubar = 0; + fH0 = 0; + fH1 = 0; + fH2 = 0; + fH3 = 0; +} +//____________________________________________________________________________ +double BYStrucFunc2021::ScalingVar(const Interaction * interaction, double Mf ) 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 Mf2 = TMath::Power( Mf, 2 ) ; + double xw = 2*x*(myQ2+Mf2+fB) / (myQ2*(1.+TMath::Sqrt(1+a)) + 2*fA*x); + return xw; +} +//____________________________________________________________________________ +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); + + double myQ2 = this->Q2(interaction); + double GD = 1. / TMath::Power(1.+myQ2/fMv2, 2); // p elastic form factor + double GD2 = TMath::Power(GD,2); + + + // 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; + + 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) + kss = myQ2/(myQ2+fCsS); // K - s(sea) +} +//____________________________________________________________________________ +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); + 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); + 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; + // 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 ) { + KLW = ( q02 + fCaLW_nu ) / q02 ; + } else if( is_nubar ) { + KLW = ( q02 + fCaLW_nubar ) / q02 ; + } + } + // 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; + +} + +//____________________________________________________________________________ +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); +} +//____________________________________________________________________________ +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 Q2_int = this->Q2(interaction) ; + if( Q2_int < 0.3 ) { + // Freeze R at Q2 = 0.3 GeV2 + Q2 = 0.3 ; + } + double Q4 = pow(Q2,2); + double Q8 = pow(Q4,2); + + const Kinematics & kinematics = interaction->Kine(); + double x = kinematics.x(); + double x2 = pow(x,2); + double x3 = pow(x,3); + + double Theta = 1 + 12.0 * ( Q2 / (Q2+1.) ) * pow(0.125,2)/(pow(0.125,2)+x2); + + // R1998 is defined as the average of Ra, Rb and Rc, each parameterized to accomodate new data at low x + // arXiv:hep-ex/9808028 + const double a1 = 0.0485; + const double a2 = 0.5470; + const double a3 = 2.0621; + const double a4 = -0.3804; + const double a5 = 0.5090; + const double a6 = -0.0285; + + double Ra = (a1 * Theta / TMath::Log(Q2/0.04) ) ; + Ra += a2 * ( 1 + a4 * x + a5 * x2 ) * pow( x, a6 ) ; + Ra /= pow( Q8 + pow(a3,4), 1./4. ) ; + + const double b1 = 0.0481; + const double b2 = 0.6114; + const double b3 = -0.3509; + const double b4 = -0.4611; + const double b5 = 0.7172; + const double b6 = -0.0317; + + double Rb = b1 * Theta / TMath::Log(Q2/0.04) ; + Rb += (b2 / Q2 + b3 / (Q4 + pow(0.3,2)) ) * (1+b4*x+b5*x2) * pow(x,b6); + + const double c1 = 0.0577; + const double c2 = 0.4644; + const double c3 = 1.8288; + const double c4 = 12.3708; + const double c5 = -43.1043; + const double c6 = 41.7415; + double Q2thr = c4*x + c5*x2 + c6*x3; + double Rc = c1 * Theta / TMath::Log( Q2/0.04 ) ; + Rc += c2 * sqrt( pow(Q2 - Q2thr,2) + pow(c3,2)); + + double R1998 = (Ra + Rb + Rc) / 3. ; + + if( Q2_int < 0.3 ) { + double Q4_int = pow(Q2_int,2); + R1998 *= 3.633 * Q2_int / ( Q4_int + 1 ) ; + } + + return R1998 ; +} + diff --git a/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h new file mode 100644 index 0000000000..21137002b2 --- /dev/null +++ b/src/Physics/DeepInelastic/XSection/BYStrucFunc2021.h @@ -0,0 +1,84 @@ +/*! + +\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, double Mf = 0 ) 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; + // 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 + 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 + 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 fCaLW_nubar; ///< C-axial neutrino LW parameter + double fCaLW_nu; ///< C-axial anti-neutrino LW parameter + double fH0; + double fH1; + double fH2; + double fH3; +}; + +} // genie namespace + +#endif // _BODEK_YANG_STRUCTURE_FUNCTION_MODEL_2021_H_ diff --git a/src/Physics/DeepInelastic/XSection/DISXSec.cxx b/src/Physics/DeepInelastic/XSection/DISXSec.cxx index dd8002194b..89cf1597e6 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,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); } @@ -88,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; @@ -97,85 +98,96 @@ 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) { - 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. - // - 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); + + 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); + delete func; + }//phase space ok? + + LOG("DISXSec", pINFO) << "XSec[DIS] (E = " << Ev << " GeV) = " << xsec; + + delete interaction; + + return xsec; + } //____________________________________________________________________________ void DISXSec::Configure(const Registry & config) @@ -207,19 +219,23 @@ 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( - 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); @@ -250,18 +266,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 @@ -314,9 +330,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(); 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 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; 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); diff --git a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.cxx index 3b1a068830..76020ee7cf 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); @@ -107,11 +107,14 @@ 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 ) ; //-- Use 2016 SF relation corrections - GetParam( "Use2016Corrections", fUse2016Corrections ) ; + GetParam( "Use2016Corrections", fUse2016Corrections, false ) ; //-- Set min for relation between 2xF1 and F2 GetParam( "LowQ2CutoffF1F2", fLowQ2CutoffF1F2 ) ; @@ -119,6 +122,8 @@ void QPMDISStrucFuncBase::LoadConfig(void) //-- turn charm production off? GetParamDef( "Charm-Prod-Off", fCharmOff, false ) ; + fDISNuclCorr = dynamic_cast(this->SubAlg("DISNuclModel")); + //-- weinberg angle double thw ; GetParam( "WeinbergAngle", thw ) ; @@ -129,7 +134,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 } @@ -183,42 +188,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; @@ -268,14 +273,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; @@ -293,33 +298,36 @@ 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; - } - - F2val = 2*(q+qbar); - xF3val = 2*(q-qbar); + 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 = (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 @@ -348,12 +356,14 @@ void QPMDISStrucFuncBase::Calculate(const Interaction * interaction) const double Q2val = this->Q2 (interaction); double x = this->ScalingVar(interaction); - double f = this->NuclMod (interaction); // nuclear modification + 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; + LOG("DISSF", pDEBUG) << "H = " << H; #endif if(fUse2016Corrections) { @@ -370,36 +380,34 @@ 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) } else { double a = TMath::Power(x,2.) / TMath::Max(Q2val, fLowQ2CutoffF1F2); double c = (1. + 4. * kNucleonMass2 * a) / (1.+r); - //double a = TMath::Power(x,2.) / Q2val; - //double c = (1. + 4. * kNucleonMass * a) / (1.+r); - - fF3 = f * xF3val / x; + + 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) - << "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(); @@ -423,68 +431,81 @@ double QPMDISStrucFuncBase::Q2(const Interaction * interaction) const return 0; } //____________________________________________________________________________ -double QPMDISStrucFuncBase::ScalingVar(const Interaction* interaction) const +double QPMDISStrucFuncBase::q0(const Interaction * interaction) const { -// The scaling variable is set to the normal Bjorken x. -// Override DISStructureFuncModel::ScalingVar() to compute corrections + 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, double Mf ) const +{ + // 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) 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 + // 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.; } + //____________________________________________________________________________ -double QPMDISStrucFuncBase::NuclMod(const Interaction * interaction) const +void QPMDISStrucFuncBase::KAxialFactors(const Interaction *, + double & kuv, double & kdv, double & kus, double & kds, double & kss ) 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(); - -// 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; -#endif - } + // 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 - return f; + kuv = 1.; + kdv = 1.; + kus = 1.; + kds = 1.; + kss = 1.; } + //____________________________________________________________________________ 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; @@ -492,6 +513,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 @@ -517,41 +542,58 @@ 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 - double kval_u = 1.; - double kval_d = 1.; - double ksea_u = 1.; - double ksea_d = 1.; + // 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); - this->KFactors(interaction, kval_u, kval_d, ksea_u, ksea_d); + // 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:"; @@ -570,15 +612,15 @@ 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); - 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 ff4782afd3..fc49f52139 100644 --- a/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h +++ b/src/Physics/DeepInelastic/XSection/QPMDISStrucFuncBase.h @@ -1,4 +1,3 @@ -//____________________________________________________________________________ /*! \class genie::QPMDISStrucFuncBase @@ -30,6 +29,7 @@ #include "Physics/DeepInelastic/XSection/DISStructureFuncModelI.h" #include "Framework/Interaction/Interaction.h" #include "Physics/PartonDistributions/PDF.h" +#include "Physics/NuclearState/DISNuclearModelI.h" namespace genie { @@ -63,17 +63,21 @@ class QPMDISStrucFuncBase : public DISStructureFuncModelI { virtual void LoadConfig (void); virtual void InitPDF (void); virtual double Q2 (const Interaction * i) const; - virtual double ScalingVar (const Interaction * i) const; + virtual double q0 (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; - virtual void KFactors (const Interaction * i, double & kuv, - double & kdv, double & kus, double & kds) const; + virtual double H (const Interaction * i) 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 + Universitat de Valencia + + For the class documentation see the corresponding header file. + +*/ +//____________________________________________________________________________ + +#include "Physics/NuclearState/BY00DISNuclearModel.h" +#include "Physics/NuclearState/NuclearUtils.h" + +using std::ostringstream; +using namespace genie; + +//____________________________________________________________________________ + +BY00DISNuclearModel::BY00DISNuclearModel() : + DISNuclearModelI("genie::BY00DISNuclearModel"){} + +//____________________________________________________________________________ + +BY00DISNuclearModel::BY00DISNuclearModel(string config) : + DISNuclearModelI("genie::BY00DISNuclearModel"){} + +//____________________________________________________________________________ + +double BY00DISNuclearModel::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 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; +#endif + + return f; + +} + +//____________________________________________________________________________ + +void BY00DISNuclearModel::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} + +//____________________________________________________________________________ + +void BY00DISNuclearModel::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} + +//____________________________________________________________________________ + +void BY00DISNuclearModel::LoadConfig(void) +{ + // 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/BY00DISNuclearModel.h b/src/Physics/NuclearState/BY00DISNuclearModel.h new file mode 100644 index 0000000000..c90d20a0da --- /dev/null +++ b/src/Physics/NuclearState/BY00DISNuclearModel.h @@ -0,0 +1,43 @@ +//____________________________________________________________________________ +/*! + +\class genie::BY00DISNuclearModel + +\brief Uses DIS Nuclear Model Correction from Bodek-Yang + +\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 _BY00_DIS_NUCLEAR_MODEL_H_ +#define _BY00_DIS_NUCLEAR_MODEL_H_ + +#include + +#include "Physics/NuclearState/DISNuclearModelI.h" + +namespace genie { + +class BY00DISNuclearModel : public DISNuclearModelI { + +public: + BY00DISNuclearModel(); + BY00DISNuclearModel(string config); + virtual ~BY00DISNuclearModel() {}; + + 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_H_ diff --git a/src/Physics/NuclearState/BY21DISNuclearModel.cxx b/src/Physics/NuclearState/BY21DISNuclearModel.cxx new file mode 100644 index 0000000000..200c7e5070 --- /dev/null +++ b/src/Physics/NuclearState/BY21DISNuclearModel.cxx @@ -0,0 +1,138 @@ + //____________________________________________________________________________ +/* + 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/BY21DISNuclearModel.h" + +using std::ostringstream; +using namespace genie; + +//____________________________________________________________________________ + +BY21DISNuclearModel::BY21DISNuclearModel() : + DISNuclearModelI("genie::BY21DISNuclearModel"){} + +//____________________________________________________________________________ + +BY21DISNuclearModel::BY21DISNuclearModel(string config) : + DISNuclearModelI("genie::BY21DISNuclearModel"){} + +//____________________________________________________________________________ + +double BY21DISNuclearModel::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*= f2HScale*( f2Hf0 + f2Hf1*xv + f2Hf2*xv2 + f2Hf3*xv3 + f2Hf4*xv4 + f2Hf5*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 *= ( 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 *= ( 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 /= ( 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__ + LOG("DISSF", pDEBUG) << "Nuclear factor for x of " << x << " = " << f; +#endif + + return f; + +} + +//____________________________________________________________________________ + +void BY21DISNuclearModel::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} + +//____________________________________________________________________________ + +void BY21DISNuclearModel::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} + +//____________________________________________________________________________ + +void BY21DISNuclearModel::LoadConfig(void) +{ + + 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-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/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..d4b40ae69c --- /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..fdb541d338 --- /dev/null +++ b/src/Physics/NuclearState/JGDISNuclearModel.cxx @@ -0,0 +1,103 @@ + //____________________________________________________________________________ +/* + 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 { + 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(); + + 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 ) ; + } + +#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) +{ + 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 ) ; + GetParam( "JG-NuclModel-a8", fAa8 ) ; +} diff --git a/src/Physics/NuclearState/JGDISNuclearModel.h b/src/Physics/NuclearState/JGDISNuclearModel.h new file mode 100644 index 0000000000..c2bf9f3d53 --- /dev/null +++ b/src/Physics/NuclearState/JGDISNuclearModel.h @@ -0,0 +1,58 @@ +//____________________________________________________________________________ +/*! + +\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); + + 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 + double fAa8; //> JG-NuclModel-a8 +}; + +} // genie namespace +#endif // _JG_DIS_NUCLEAR_MODEL_H_ diff --git a/src/Physics/NuclearState/LinkDef.h b/src/Physics/NuclearState/LinkDef.h index c9207f1ca8..c4a74ac043 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 BY00DISNuclearModel; +#pragma link C++ class BY21DISNuclearModel; +#pragma link C++ class JGDISNuclearModel; #endif