GENIEGenerator
Loading...
Searching...
No Matches
genie::MKSPPPXSec2020 Class Reference


Class calculate differental cross-sections More...

#include <MKSPPPXSec2020.h>

Inheritance diagram for genie::MKSPPPXSec2020:
[legend]
Collaboration diagram for genie::MKSPPPXSec2020:
[legend]

Classes

class  HelicityAmpVminusARes
class  HelicityBkgAmp
struct  is_complex
struct  is_complex< std::complex< T > >
class  Iterator
class  SumHelicityAmpVminusARes

Public Member Functions

 MKSPPPXSec2020 ()
 MKSPPPXSec2020 (string config)
virtual ~MKSPPPXSec2020 ()
double XSec (const Interaction *i, KinePhaseSpace_t k) const
 Compute the cross section for the input interaction.
double Integral (const Interaction *i) const
bool ValidProcess (const Interaction *i) const
 Can this cross section algorithm handle the input process?
bool ValidKinematics (const Interaction *interaction) const
 Is the input kinematical point a physically allowed one?
void Configure (const Registry &config)
void Configure (string config)
Public Member Functions inherited from genie::XSecAlgorithmI
virtual ~XSecAlgorithmI ()
Public Member Functions inherited from genie::Algorithm
virtual ~Algorithm ()
virtual void FindConfig (void)
virtual const RegistryGetConfig (void) const
RegistryGetOwnedConfig (void)
virtual const AlgIdId (void) const
 Get algorithm ID.
virtual AlgStatus_t GetStatus (void) const
 Get algorithm status.
virtual bool AllowReconfig (void) const
virtual AlgCmp_t Compare (const Algorithm *alg) const
 Compare with input algorithm.
virtual void SetId (const AlgId &id)
 Set algorithm ID.
virtual void SetId (string name, string config)
const AlgorithmSubAlg (const RgKey &registry_key) const
void AdoptConfig (void)
void AdoptSubstructure (void)
virtual void Print (ostream &stream) const
 Print algorithm info.

Private Types

enum  Current { VECTOR , AXIAL }
enum  BosonPolarization { LEFT , RIGHT , MINUS0 , PLUS0 }
enum  NucleonPolarization { MINUS , PLUS }
using CurrentIterator = Iterator<Current, Current::VECTOR, Current::AXIAL>
using BosonPolarizationIterator = Iterator<BosonPolarization, BosonPolarization::LEFT, BosonPolarization::PLUS0>
using NucleonPolarizationIterator = Iterator<NucleonPolarization, NucleonPolarization::MINUS, NucleonPolarization::PLUS>
template<bool C, typename T = void>
using enable_if_t = typename std::enable_if<C, T>::type

Private Member Functions

int Lambda (NucleonPolarization l) const
int Lambda (BosonPolarization l) const
int PhaseFactor (BosonPolarization lk, NucleonPolarization l1, NucleonPolarization l2) const
void LoadConfig (void)

Private Attributes

FKR fFKR
const RSHelicityAmplModelIfHAmplModelCC
const RSHelicityAmplModelIfHAmplModelNCp
const RSHelicityAmplModelIfHAmplModelNCn
double fFermiConstant
double fCA50
 FKR parameter Zeta.
double fOmega
 FKR parameter Omega.
double fMa2
 (axial mass)^2
double fMv2
 (vector mass)^2
double fCv3
 GV calculation coeffs.
double fCv4
double fCv51
double fCv52
double fSin2Wein
 sin^2(Weingberg angle)
double fVud
 |Vud| (magnitude ud-element of CKM-matrix)
double fXSecScaleCC
 External CC xsec scaling factor.
double fXSecScaleNC
 External NC xsec scaling factor.
const ELFormFactorsModelIfElFFModel
 Elastic form facors model for background contribution.
const QELFormFactorsModelIfFormFactorsModel
 Quasielastic form facors model for background contribution.
const QELFormFactorsModelIfEMFormFactorsModel
 Electromagnetic form factors model for background contribution.
string fKFTable
 Table of Fermi momentum (kF) constants for various nuclei.
bool fUseRFGParametrization
 Use parametrization for fermi momentum insted of table?
bool fUsePauliBlocking
 Account for Pauli blocking?
QELFormFactors fFormFactors
 Quasielastic form facors for background contribution.
QELFormFactors fEMFormFactors
 Electromagnetic form facors for background contribution.
double f_pi
 Constant for pion-nucleon interaction.
double FA0
 Axial coupling (value of axial form factor at Q2=0)
double Frho0
double fBkgVWmin
double fBkgVWmax
double fBkgV3
double fBkgV2
double fBkgV1
double fBkgV0
double fRho770Mass
 Mass of rho(770) meson.
double fWmax
 The value above which the partial cross section is set to zero (if negative then not appliciable)
bool fUseAuthorCode
 Use author code?
const XSecIntegratorIfXSecIntegrator
BaryonResList fResList

Additional Inherited Members

Static Public Member Functions inherited from genie::Algorithm
static string BuildParamVectKey (const std::string &comm_name, unsigned int i)
static string BuildParamVectSizeKey (const std::string &comm_name)
static string BuildParamMatKey (const std::string &comm_name, unsigned int i, unsigned int j)
static string BuildParamMatRowSizeKey (const std::string &comm_name)
static string BuildParamMatColSizeKey (const std::string &comm_name)
Protected Member Functions inherited from genie::XSecAlgorithmI
 XSecAlgorithmI ()
 XSecAlgorithmI (string name)
 XSecAlgorithmI (string name, string config)
Protected Member Functions inherited from genie::Algorithm
 Algorithm ()
 Algorithm (string name)
 Algorithm (string name, string config)
void Initialize (void)
void DeleteConfig (void)
void DeleteSubstructure (void)
RegistryExtractLocalConfig (const Registry &in) const
RegistryExtractLowerConfig (const Registry &in, const string &alg_key) const
 Split an incoming configuration Registry into a block valid for the sub-algo identified by alg_key.
template<class T>
bool GetParam (const RgKey &name, T &p, bool is_top_call=true) const
template<class T>
bool GetParamDef (const RgKey &name, T &p, const T &def) const
template<class T>
int GetParamVect (const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
 Handle to load vectors of parameters.
int GetParamVectKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const
template<class T>
int GetParamMat (const std::string &comm_name, TMatrixT< T > &mat, bool is_top_call=true) const
 Handle to load matrix of parameters.
template<class T>
int GetParamMatSym (const std::string &comm_name, TMatrixTSym< T > &mat, bool is_top_call=true) const
int GetParamMatKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const
int AddTopRegistry (Registry *rp, bool owns=true)
 add registry with top priority, also update ownership
int AddLowRegistry (Registry *rp, bool owns=true)
 add registry with lowest priority, also update ownership
int MergeTopRegistry (const Registry &r)
int AddTopRegisties (const vector< Registry * > &rs, bool owns=false)
 Add registries with top priority, also udated Ownerships.
Protected Attributes inherited from genie::Algorithm
bool fAllowReconfig
bool fOwnsSubstruc
 true if it owns its substructure (sub-algs,...)
AlgId fID
 algorithm name and configuration set
vector< Registry * > fConfVect
vector< bool > fOwnerships
 ownership for every registry in fConfVect
AlgStatus_t fStatus
 algorithm execution status
AlgMapfOwnedSubAlgMp
 local pool for owned sub-algs (taken out of the factory pool)

Detailed Description


Class calculate differental cross-sections

\[\frac{d^4\sigma}{dQ^2dWd\cos\theta_\pi d\phi_\pi}
\]

or

\[\frac{d^3\sigma}{dQ^2dWd\cos\theta_\pi}
\]

for specific neutrino energy (in lab frame), where:

Variable Description
$W$ Invariant mass
$Q^2$ Sqaured 4-momentum transfer
$\cos\theta_\pi$ Cosine of pion polar angle in $\pi$N rest frame
$\phi_\pi$ Pion azimuthal angle in $\pi$N rest frame

for the following channels:

  1. $\nu            + p \to \ell^-          + p + \pi^+$
  2. $\nu            + n \to \ell^-          + p + \pi^0$
  3. $\nu            + n \to \ell^-          + n + \pi^+$
  4. $\overline{\nu} + n \to \ell^+          + n + \pi^-$
  5. $\overline{\nu} + p \to \ell^+          + n + \pi^0$
  6. $\overline{\nu} + p \to \ell^+          + p + \pi^-$
  7. $\nu            + p \to \nu             + p + \pi^0$
  8. $\nu            + p \to \nu             + n + \pi^+$
  9. $\nu            + n \to \nu             + n + \pi^0$
  10. $\nu            + n \to \nu             + p + \pi^-$
  11. $\overline{\nu} + p \to \overline{\nu}  + p + \pi^0$
  12. $\overline{\nu} + p \to \overline{\nu}  + n + \pi^+$
  13. $\overline{\nu} + n \to \overline{\nu}  + n + \pi^0$
  14. $\overline{\nu} + n \to \overline{\nu}  + p + \pi^-$
References:\n

  1. Kabirnezhad M., Phys.Rev.D 97 (2018) 013002 (Erratim: arXiv:1711.02403)
  2. Kabirnezhad M., Ph. D. Thesis (https://www.ncbj.gov.pl/sites/default/files/m.kabirnezhad_thesis_0.pdf , https://www.ncbj.gov.pl/sites/default/files/m.kabirnezhad-thesis_final_0.pdf)
  3. Rein D., Sehgal L., Ann. of Phys. 133 (1981) 79-153
  4. Rein D., Z.Phys.C 35 (1987) 43-64
  5. Adler S.L., Ann. Phys. 50 (1968) 189
  6. Graczyk K., Sobczyk J., Phys.Rev.D 77 (2008) 053001 [Erratum: ibid.D 79 (2009) 079903]
  7. Jacob M., Wick G.C., Ann. of Phys. 7 (1959) 404-428
  8. Hernandez E., Nieves J., Valverde M., Phys.Rev.D 76 (2007) 033005
  9. Adler S.L., Nussinov S., Paschos E.A., Phys. Rev. D 9 (1974) 2125-2143 [Erratum: ibid D 10 (1974) 1669]
  10. Paschos E.A., Yu J.Y., Sakuda M., Phys. Rev. D 69 (2004) 014013 [arXiv: hep-ph/0308130]
  11. Yu J.Y., "Neutrino interactions and nuclear effects in oscillation experiments and the nonperturbative dispersive sector in strong (quasi-)abelian fields", Ph. D.Thesis, Dortmund U., Dortmund, 2002 (unpublished)
  12. Kakorin I., Kuzmin K., Naumov V. "Report on implementation of the MK-model for resonance single-pion production into GENIE" (https://genie-docdb.pp.rl.ac.uk/cgi-bin/private/ShowDocument?docid=181, http://theor.jinr.ru/NeutrinoOscillations/Papers/Report_MK_implementation_GENIE.pdf)
Authors
Igor Kakorin kakor.nosp@m.in@j.nosp@m.inr.r.nosp@m.u, Joint Institute for Nuclear Research
Vadim Naumov vnaum.nosp@m.ov@t.nosp@m.heor..nosp@m.jinr.nosp@m..ru, Joint Institute for Nuclear Research
adapted from code provided by
Minoo Kabirnezhad minoo.nosp@m..kab.nosp@m.irnez.nosp@m.had@.nosp@m.physi.nosp@m.cs.o.nosp@m.x.ac..nosp@m.uk University of Oxford, Department of Physics
based on code of
Costas Andreopoulos c.and.nosp@m.reop.nosp@m.oulos.nosp@m.@cer.nosp@m.n.ch University of Liverpool
Created:\n Nov 12, 2019
License:\n Copyright (c) 2003-2025, The GENIE Collaboration
For the full text of the license visit http://copyright.genie-mc.org or see $GENIE/LICENSE

Definition at line 99 of file MKSPPPXSec2020.h.

Member Typedef Documentation

◆ BosonPolarizationIterator

◆ CurrentIterator

◆ enable_if_t

template<bool C, typename T = void>
using genie::MKSPPPXSec2020::enable_if_t = typename std::enable_if<C, T>::type
private

Definition at line 164 of file MKSPPPXSec2020.h.

◆ NucleonPolarizationIterator

Member Enumeration Documentation

◆ BosonPolarization

◆ Current

Enumerator
VECTOR 
AXIAL 

Definition at line 128 of file MKSPPPXSec2020.h.

◆ NucleonPolarization

Enumerator
MINUS 
PLUS 

Definition at line 130 of file MKSPPPXSec2020.h.

Constructor & Destructor Documentation

◆ MKSPPPXSec2020() [1/2]

MKSPPPXSec2020::MKSPPPXSec2020 ( )

Definition at line 52 of file MKSPPPXSec2020.cxx.

52 :
53XSecAlgorithmI("genie::MKSPPPXSec2020")
54{
55
56}

References genie::XSecAlgorithmI::XSecAlgorithmI().

◆ MKSPPPXSec2020() [2/2]

MKSPPPXSec2020::MKSPPPXSec2020 ( string config)

Definition at line 58 of file MKSPPPXSec2020.cxx.

58 :
59XSecAlgorithmI("genie::MKSPPPXSec2020", config)
60{
61
62}

References genie::XSecAlgorithmI::XSecAlgorithmI().

◆ ~MKSPPPXSec2020()

MKSPPPXSec2020::~MKSPPPXSec2020 ( )
virtual

Definition at line 64 of file MKSPPPXSec2020.cxx.

65{
66
67}

Member Function Documentation

◆ Configure() [1/2]

void MKSPPPXSec2020::Configure ( const Registry & config)
virtual

Configure the algorithm with an external registry The registry is merged with the top level registry if it is owned, Otherwise a copy of it is added with the highest priority

Reimplemented from genie::Algorithm.

Definition at line 1374 of file MKSPPPXSec2020.cxx.

1375{
1376 Algorithm::Configure(config);
1377 this->LoadConfig();
1378}
virtual void Configure(const Registry &config)
Definition Algorithm.cxx:62

References genie::Algorithm::Configure(), and LoadConfig().

◆ Configure() [2/2]

void MKSPPPXSec2020::Configure ( string config)
virtual

Configure the algorithm from the AlgoConfigPool based on param_set string given in input An algorithm contains a vector of registries coming from different xml configuration files, which are loaded according a very precise prioriy This methods will load a number registries in order of priority: 1) "Tunable" parameter set from CommonParametes. This is loaded with the highest prioriry and it is designed to be used for tuning procedure Usage not expected from the user. 2) For every string defined in "CommonParame" the corresponding parameter set will be loaded from CommonParameter.xml 3) parameter set specified by the config string and defined in the xml file of the algorithm 4) if config is not "Default" also the Default parameter set from the same xml file will be loaded Effectively this avoids the repetion of a parameter when it is not changed in the requested configuration

Reimplemented from genie::Algorithm.

Definition at line 1380 of file MKSPPPXSec2020.cxx.

1381{
1382 Algorithm::Configure(config);
1383 this->LoadConfig();
1384}

References genie::Algorithm::Configure(), and LoadConfig().

◆ Integral()

double MKSPPPXSec2020::Integral ( const Interaction * i) const
virtual

Integrate the model over the kinematic phase space available to the input interaction (kinematical cuts can be included)

Implements genie::XSecAlgorithmI.

Definition at line 1306 of file MKSPPPXSec2020.cxx.

1307{
1308 double xsec = fXSecIntegrator->Integrate(this,interaction);
1309 return xsec;
1310}
const XSecIntegratorI * fXSecIntegrator

References fXSecIntegrator.

◆ Lambda() [1/2]

int MKSPPPXSec2020::Lambda ( BosonPolarization l) const
inlineprivate

Definition at line 1332 of file MKSPPPXSec2020.cxx.

1333{
1334 return (l*(l*(4*l-21)+29)-6)/3;
1335}

◆ Lambda() [2/2]

int MKSPPPXSec2020::Lambda ( NucleonPolarization l) const
inlineprivate

Definition at line 1327 of file MKSPPPXSec2020.cxx.

1328{
1329 return 2*l - 1;
1330}

Referenced by XSec().

◆ LoadConfig()

void MKSPPPXSec2020::LoadConfig ( void )
private

Definition at line 1386 of file MKSPPPXSec2020.cxx.

1387{
1388
1389 fResList.Clear();
1390 string resonances ;
1391 GetParam( "ResonanceNameList", resonances ) ;
1392 fResList.DecodeFromNameList(resonances);
1393
1394 // Cross section scaling symmetry_factors
1395 this->GetParam( "RES-CC-XSecScale", fXSecScaleCC );
1396 this->GetParam( "RES-NC-XSecScale", fXSecScaleNC );
1397
1398 // Load all configuration data or set defaults
1399 this->GetParamDef( "RES-Omega" , fOmega, 1.05);
1400
1401 double ma, mv;
1402 this->GetParam( "RES-Ma" , ma);
1403 this->GetParam( "RES-Mv" , mv);
1404 fMa2 = ma*ma;
1405 fMv2 = mv*mv;
1406 this->GetParam( "RES-CA50" , fCA50) ;
1407
1408 this->GetParam( "GVCAL-Cv3" , fCv3) ;
1409 this->GetParam( "GVCAL-Cv4" , fCv4) ;
1410 this->GetParam( "GVCAL-Cv51" , fCv51) ;
1411 this->GetParam( "GVCAL-Cv52" , fCv52) ;
1412
1413 this->GetParam( "FermiConstant", fFermiConstant );
1414
1415 double thw;
1416 this->GetParam( "WeinbergAngle", thw );
1417 fSin2Wein = TMath::Power( TMath::Sin(thw), 2 );
1418
1419 this->GetParam("CKM-Vud", fVud );
1420
1421 // Load all the sub-algorithms needed
1422 fHAmplModelCC = 0;
1423 fHAmplModelNCp = 0;
1424 fHAmplModelNCn = 0;
1425
1426 fHAmplModelCC = dynamic_cast<const RSHelicityAmplModelI *> ( this -> SubAlg( "HelictyAmplCCAlg" ) );
1427 fHAmplModelNCp = dynamic_cast<const RSHelicityAmplModelI *> ( this -> SubAlg( "HelictyAmplNCpAlg" ));
1428 fHAmplModelNCn = dynamic_cast<const RSHelicityAmplModelI *> ( this -> SubAlg( "HelictyAmplNCnAlg" ));
1429
1430 assert(fHAmplModelCC);
1431 assert(fHAmplModelNCp);
1432 assert(fHAmplModelNCn);
1433
1434 // Parameters for calculation of background contribution
1435 fFormFactorsModel = dynamic_cast<const QELFormFactorsModelI *> (this->SubAlg("CCFormFactorsAlg"));
1436 assert(fFormFactorsModel);
1437 fFormFactors.SetModel(fFormFactorsModel); // <-- attach algorithm
1438
1439 fEMFormFactorsModel = dynamic_cast<const QELFormFactorsModelI *> (this->SubAlg("EMFormFactorsAlg"));
1440 assert(fEMFormFactorsModel);
1441 fEMFormFactors.SetModel(fEMFormFactorsModel); // <-- attach algorithm
1442
1443 this->GetParam("FermiMomentumTable", fKFTable);
1444 this->GetParam("RFG-UseParametrization", fUseRFGParametrization);
1445 this->GetParam("UsePauliBlockingForRES", fUsePauliBlocking);
1446
1447
1448 this->GetParamDef("MK-fPi", f_pi, 0.093 );
1449 this->GetParamDef("QEL-FA0", FA0, -1.26 );
1450 this->GetParamDef("MK-Frho0", Frho0, 1.0 );
1451
1452 this->GetParamDef("MK-NonResBkg-VWmin", fBkgVWmin, 1.30 );
1453 this->GetParamDef("MK-NonResBkg-VWmax", fBkgVWmax, 1.60 );
1454 this->GetParamDef("MK-NonResBkg-V3", fBkgV3, 8.08497 );
1455 this->GetParamDef("MK-NonResBkg-V2", fBkgV2, -41.6767 );
1456 this->GetParamDef("MK-NonResBkg-V1", fBkgV1, 66.3419 );
1457 this->GetParamDef("MK-NonResBkg-V0", fBkgV0, -32.5733 );
1458 this->GetParamDef("Mass-rho770", fRho770Mass, 0.7758 );
1459 this->GetParamDef("MK-WMax", fWmax, -1.0 );
1460
1461 this->GetParamDef("UseAuthorCode", fUseAuthorCode, false );
1462
1463 // Load the differential cross section integrator
1464 fXSecIntegrator = dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
1465 assert(fXSecIntegrator);
1466
1467}
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
bool GetParamDef(const RgKey &name, T &p, const T &def) const
const Algorithm * SubAlg(const RgKey &registry_key) const
double fRho770Mass
Mass of rho(770) meson.
const RSHelicityAmplModelI * fHAmplModelNCp
double fSin2Wein
sin^2(Weingberg angle)
bool fUseRFGParametrization
Use parametrization for fermi momentum insted of table?
const RSHelicityAmplModelI * fHAmplModelNCn
bool fUsePauliBlocking
Account for Pauli blocking?
double fCA50
FKR parameter Zeta.
double fMa2
(axial mass)^2
double fMv2
(vector mass)^2
const RSHelicityAmplModelI * fHAmplModelCC
const QELFormFactorsModelI * fFormFactorsModel
Quasielastic form facors model for background contribution.
const QELFormFactorsModelI * fEMFormFactorsModel
Electromagnetic form factors model for background contribution.
double fWmax
The value above which the partial cross section is set to zero (if negative then not appliciable)
double fXSecScaleNC
External NC xsec scaling factor.
QELFormFactors fEMFormFactors
Electromagnetic form facors for background contribution.
double fVud
|Vud| (magnitude ud-element of CKM-matrix)
double fCv3
GV calculation coeffs.
double fXSecScaleCC
External CC xsec scaling factor.
double FA0
Axial coupling (value of axial form factor at Q2=0)
string fKFTable
Table of Fermi momentum (kF) constants for various nuclei.
bool fUseAuthorCode
Use author code?
double f_pi
Constant for pion-nucleon interaction.
double fOmega
FKR parameter Omega.
QELFormFactors fFormFactors
Quasielastic form facors for background contribution.

References f_pi, FA0, fBkgV0, fBkgV1, fBkgV2, fBkgV3, fBkgVWmax, fBkgVWmin, fCA50, fCv3, fCv4, fCv51, fCv52, fEMFormFactors, fEMFormFactorsModel, fFermiConstant, fFormFactors, fFormFactorsModel, fHAmplModelCC, fHAmplModelNCn, fHAmplModelNCp, fKFTable, fMa2, fMv2, fOmega, fResList, Frho0, fRho770Mass, fSin2Wein, fUseAuthorCode, fUsePauliBlocking, fUseRFGParametrization, fVud, fWmax, fXSecIntegrator, fXSecScaleCC, fXSecScaleNC, genie::Algorithm::GetParam(), genie::Algorithm::GetParamDef(), and genie::Algorithm::SubAlg().

Referenced by Configure(), and Configure().

◆ PhaseFactor()

int MKSPPPXSec2020::PhaseFactor ( BosonPolarization lk,
NucleonPolarization l1,
NucleonPolarization l2 ) const
inlineprivate

Definition at line 1337 of file MKSPPPXSec2020.cxx.

1338{
1339 return lk*(lk*(4*lk - 21) + 29)/6 - l1 - l2;
1340}

Referenced by XSec().

◆ ValidKinematics()

bool MKSPPPXSec2020::ValidKinematics ( const Interaction * i) const
virtual

Is the input kinematical point a physically allowed one?

Reimplemented from genie::XSecAlgorithmI.

Definition at line 1342 of file MKSPPPXSec2020.cxx.

1343{
1344 // call only after ValidProcess
1345 if ( interaction->TestBit(kISkipKinematicChk) ) return true;
1346
1347 const KPhaseSpace& kps = interaction->PhaseSpace();
1348
1349 // Get kinematical parameters
1350 const InitialState & init_state = interaction -> InitState();
1351 const Kinematics & kinematics = interaction -> Kine();
1352 double Enu = init_state.ProbeE(kRfHitNucRest);
1353 double W = kinematics.W();
1354 double Q2 = kinematics.Q2();
1355
1356 if (Enu < kps.Threshold_SPP_iso())
1357 return false;
1358
1359 Range1D_t Wl = kps.WLim_SPP_iso();
1360 if (fWmax >= Wl.min)
1361 Wl.max = TMath::Min(fWmax, Wl.max);
1362 Range1D_t Q2l = kps.Q2Lim_W_SPP_iso();
1363
1364 if (W < Wl.min || W > Wl.max)
1365 return false;
1366
1367 if (Q2 < Q2l.min || Q2 > Q2l.max)
1368 return false;
1369
1370 return true;
1371
1372}
double ProbeE(RefFrame_t rf) const
double Threshold_SPP_iso(void) const
Energy limit for resonance single pion production on isoscalar nucleon.
Range1D_t WLim_SPP_iso(void) const
W limits for resonance single pion production on isoscalar nucleon.
Range1D_t Q2Lim_W_SPP_iso(void) const
Q2 limits @ fixed W for resonance single pion production on isoscalar nucleon.
double Q2(bool selected=false) const
double W(bool selected=false) const
double W(const Interaction *const i)
double Q2(const Interaction *const i)
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition Interaction.h:48
@ kRfHitNucRest
Definition RefFrame.h:30

References fWmax, genie::kISkipKinematicChk, genie::kRfHitNucRest, genie::Range1D_t::max, genie::Range1D_t::min, genie::Interaction::PhaseSpace(), genie::InitialState::ProbeE(), genie::Kinematics::Q2(), genie::KPhaseSpace::Q2Lim_W_SPP_iso(), genie::KPhaseSpace::Threshold_SPP_iso(), genie::Kinematics::W(), and genie::KPhaseSpace::WLim_SPP_iso().

Referenced by XSec().

◆ ValidProcess()

bool MKSPPPXSec2020::ValidProcess ( const Interaction * i) const
virtual

Can this cross section algorithm handle the input process?

Implements genie::XSecAlgorithmI.

Definition at line 1312 of file MKSPPPXSec2020.cxx.

1313{
1314
1315 if(interaction->TestBit(kISkipProcessChk)) return true;
1316
1317 //-- Get the requested SPP channel
1318 SppChannel_t spp_channel = SppChannel::FromInteraction(interaction);
1319 if( spp_channel == kSppNull ) {
1320 return false;
1321 }
1322
1323 return true;
1324
1325}
static SppChannel_t FromInteraction(const Interaction *interaction)
Definition SppChannel.h:402
@ kSppNull
Definition SppChannel.h:46
enum genie::ESppChannel SppChannel_t
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition Interaction.h:47

References genie::SppChannel::FromInteraction(), genie::kISkipProcessChk, and genie::kSppNull.

Referenced by XSec().

◆ XSec()

double MKSPPPXSec2020::XSec ( const Interaction * i,
KinePhaseSpace_t k ) const
virtual

Compute the cross section for the input interaction.

Implements genie::XSecAlgorithmI.

Definition at line 69 of file MKSPPPXSec2020.cxx.

70{
71 if(! this -> ValidProcess (interaction) ) return 0;
72 if(! this -> ValidKinematics (interaction) ) return 0;
73
74 const InitialState & init_state = interaction -> InitState();
75 const ProcessInfo & proc_info = interaction -> ProcInfo();
76 const Target & target = init_state.Tgt();
77
78 double xsec = 0;
79 //-- Get 1pi exclusive channel
80 SppChannel_t spp_channel = SppChannel::FromInteraction(interaction);
81
82 int nucpdgc = target.HitNucPdg();
83 int probepdgc = init_state.ProbePdg();
84 bool is_nu = pdg::IsNeutrino (probepdgc);
85 bool is_nubar = pdg::IsAntiNeutrino (probepdgc);
86 bool is_CC = proc_info.IsWeakCC();
87 bool is_NC = proc_info.IsWeakNC();
88 bool is_p = pdg::IsProton (nucpdgc);
89
90
91 // Get kinematical parameters
92 const Kinematics & kinematics = interaction -> Kine();
93 double E = init_state.ProbeE(kRfHitNucRest);
94 double ml = interaction->FSPrimLepton()->Mass();
95 double ml2 = ml*ml;
96 double Q2 = kinematics.Q2();
97 double W = kinematics.W();
98 double W2 = W*W;
99 double Wt2 = W*2;
100
101
102 // dimension of kine phase space
103 std::string s = KinePhaseSpace::AsString(kps);
104 int kpsdim = s!="<|E>"?1 + std::count(s.begin(), s.begin()+s.find('}'), ','):0;
105 if (kpsdim < 3 || kpsdim > 4) return 0;
106 // Pion angles should be given in Adler frame
107 double CosTheta = kinematics.GetKV(kKVctp);
108 double Phi = 0;
109 if (kpsdim == 4)
110 Phi = kinematics.GetKV(kKVphip);
111
112 double SinTheta = TMath::Sqrt(1 - CosTheta*CosTheta);
113 double SinHalfTheta = TMath::Sqrt((1 - CosTheta)/2);
114 double CosHalfTheta = TMath::Sqrt((1 + CosTheta)/2);
115
116 PDGLibrary * pdglib = PDGLibrary::Instance();
117
118 // imply isospin symmetry
119 double m_pi = (pdglib->Find(kPdgPiP)->Mass() + pdglib->Find(kPdgPi0)->Mass() + pdglib->Find(kPdgPiM)->Mass())/3;
120 double m_pi2 = m_pi*m_pi;
121
122 double M = (pdglib->Find(kPdgProton)->Mass() + pdglib->Find(kPdgNeutron)->Mass())/2;
123 double M2 = M*M;
124 double Mt2 = M*2;
125
126 double q_0 = (W2 - M2 + m_pi2)/Wt2;
127 double q_02 = q_0*q_0;
128 double k_0 = (W2 - M2 - Q2)/Wt2;
129 double abs_mom_q = TMath::Sqrt(TMath::Max(0., q_02 - m_pi2));
130 double abs_mom_k = TMath::Sqrt(k_0*k_0 + Q2);
131 //double E_2L = (M2 - W2 - Q2 + Mt2*E)/Mt2;
132 double abs_mom_k_L = W*abs_mom_k/M;
133 double abs_mom_k_L2 = abs_mom_k_L*abs_mom_k_L;
134 double k_2 = (M2 + Mt2*E - W2 - ml2)/Wt2;
135 double k_1 = k_2 + k_0;
136 double p_10 = (W2 + M2 + Q2)/Wt2;
137 double qk = q_0*k_0 - abs_mom_k*abs_mom_q*CosTheta;
138 //double k_2L = TMath::Sqrt(E_2L*E_2L - ml2); //magnitude of lepton momentum in lab frame
139
140 // There is no check of positivity under root expression in the original code
141 double k_2_iso = TMath::Sqrt(TMath::Max(0., k_2*k_2 - ml2)); //magnitude of lepton momentum in isobaric frame
142 double cos_theta = k_2_iso>0?(2*k_1*k_2 - Q2 - ml2)/2/k_1/k_2_iso:0;
143 // There are no checks presented below in the original code
144 if (cos_theta > 1)
145 cos_theta = 1;
146 if (cos_theta < -1)
147 cos_theta = -1;
148
149 // Eq. 7 of ref. 1
150 double A_plus = TMath::Sqrt( k_1*(k_2 - k_2_iso) );
151 double A_minus = TMath::Sqrt( k_1*(k_2 + k_2_iso) );
152 //Eq. 6 of ref. 1
153 double eps_1_plus = 2*A_plus *(k_1 - k_2_iso)/abs_mom_k*TMath::Sqrt(1 + cos_theta);
154 double eps_1_minus = -2*A_minus*(k_1 + k_2_iso)/abs_mom_k*TMath::Sqrt(1 - cos_theta);
155 double eps_2_plus = 2*A_plus *TMath::Sqrt(1 + cos_theta);
156 double eps_2_minus = 2*A_minus*TMath::Sqrt(1 - cos_theta);
157 //Eq. 9 of ref. 1
158 double eps_zero_L = -2*A_minus*TMath::Sqrt(1 + cos_theta); // L->lambda = -1
159 double eps_zero_R = 2*A_plus *TMath::Sqrt(1 - cos_theta); // R->lambda = +1
160 double eps_z_L = -2*A_minus*(k_1 - k_2_iso)/abs_mom_k*TMath::Sqrt(1 + cos_theta);
161 double eps_z_R = 2*A_plus *(k_1 + k_2_iso)/abs_mom_k*TMath::Sqrt(1 - cos_theta);
162
163 if (is_nubar)
164 {
165 if (fUseAuthorCode)
166 {
167 //This ``recipe'' of transition from neutrino to antineutrino case is promoted by Minoo Kabirnezhad.
168 //However, it is not correct in our opinion. All details can be found in Ref. [12], see
169 //section "Problem with transition from neutrino to antineutrino case".
170 Phi = -Phi;
171 eps_1_plus = -2*A_minus *(k_1 + k_2_iso)/abs_mom_k*TMath::Sqrt(1 - cos_theta);
172 eps_1_minus = -2*A_plus*(k_1 - k_2_iso)/abs_mom_k*TMath::Sqrt(1 + cos_theta);
173 eps_2_plus = -2*A_minus *TMath::Sqrt(1 - cos_theta);
174 eps_2_minus = 2*A_plus*TMath::Sqrt(1 + cos_theta);
175 eps_zero_L = 2*A_plus*TMath::Sqrt(1 - cos_theta); // L->lambda = -1
176 eps_zero_R = 2*A_minus *TMath::Sqrt(1 + cos_theta); // R->lambda = +1
177 eps_z_L = 2*A_plus*(k_1 + k_2_iso)/abs_mom_k*TMath::Sqrt(1 - cos_theta);
178 eps_z_R = 2*A_minus *(k_1 - k_2_iso)/abs_mom_k*TMath::Sqrt(1 + cos_theta);
179 }
180 else
181 {
182 eps_1_plus = 2*A_minus *(k_1 + k_2_iso)/abs_mom_k*TMath::Sqrt(1 - cos_theta);
183 eps_1_minus = 2*A_plus*(k_1 - k_2_iso)/abs_mom_k*TMath::Sqrt(1 + cos_theta);
184 eps_2_plus = 2*A_minus *TMath::Sqrt(1 - cos_theta);
185 eps_2_minus = -2*A_plus*TMath::Sqrt(1 + cos_theta);
186 eps_zero_L = 2*A_plus*TMath::Sqrt(1 - cos_theta); // L->lambda = -1
187 eps_zero_R = 2*A_minus *TMath::Sqrt(1 + cos_theta); // R->lambda = +1
188 eps_z_L = 2*A_plus*(k_1 + k_2_iso)/abs_mom_k*TMath::Sqrt(1 - cos_theta);
189 eps_z_R = 2*A_minus *(k_1 - k_2_iso)/abs_mom_k*TMath::Sqrt(1 + cos_theta);
190 }
191 }
192 //Eq. 10 of ref. 1
193 double C_L_plus = k1_Sqrt2*(eps_1_plus - eps_2_plus);
194 double C_L_minus = k1_Sqrt2*(eps_1_minus - eps_2_minus);
195 double C_R_plus = -k1_Sqrt2*(eps_1_plus + eps_2_plus);
196 double C_R_minus = -k1_Sqrt2*(eps_1_minus + eps_2_minus);
197 double C_S_plus = TMath::Sqrt(TMath::Abs(eps_zero_R*eps_zero_R - eps_z_R*eps_z_R));
198 double C_S_minus = TMath::Sqrt(TMath::Abs(eps_zero_L*eps_zero_L - eps_z_L*eps_z_L));
199
200 // if C_S_plus or C_S_minus are equal to zero, then the singularity appers
201 // in the expressions for Hbkg and Hres at BosonPolarization=PLUS0 or BosonPolarization=MINUS0
202 // which further eliminated by corresponding factors in the sum for calculating of the cross section.
203 // Therefore we can just set them equal to 1 in this case. In our opinion it is simplier to eliminate
204 // singularity in such way, because it allows to keep intact original formulas from paper.
205 C_S_plus = C_S_plus==0?1:C_S_plus;
206 C_S_minus = C_S_minus==0?1:C_S_minus;
207
208
209 /* Bkg Contribution */
210 //Auxiliary functions
211 double W_plus = W + M;
212 double W_minus = W - M;
213 double W_plus2 = W_plus*W_plus;
214 double W_minus2 = W_minus*W_minus;
215 double s_minus_M2 = W_plus*W_minus;
216 double u_minus_M2 = m_pi2 - 2*(q_0*p_10 + abs_mom_q*abs_mom_k*CosTheta);
217 double t_minus_mpi2 = -(Q2 + 2*qk);
218 double mpi2_minus_k2 = Q2 + m_pi2;
219 double Q = TMath::Sqrt(Q2);
220
221 // Helicity amplitudes for bkg
223
224 //Vector_helicity amplituds for bkg
225 // vector cut
226 double FV_cut = 0;
227 //virtual form symmetry_factor to kill the bkg smothly
228 if (W < fBkgVWmin)
229 FV_cut = 1;
230 else if (W >= fBkgVWmin && W < fBkgVWmax)
231 FV_cut = fBkgV0 + W*(fBkgV1 + W*(fBkgV2 + W*fBkgV3));
232
233 fFormFactors.Calculate(interaction);
234 double g_A = -FA0;
235 double FF_pi = fFormFactors.F1V();
236 double FF_1 = 0.5*FF_pi;
237 double FF_2 = 0.5*fFormFactors.xiF2V();
238
239
240
241 // Eq. 38 and Table 5 of ref. 1
242 double V_1 = 0, V_2 = 0, V_3 = 0, V_4 = 0, V_5 = 0;
243
244 // [p,n,pi+,pi0,pi-]
245 switch (spp_channel) {
246
247 //CC
248 case (kSpp_vp_cc_10100) :
249 case (kSpp_vbn_cc_01001):
250 V_1 = g_A*kSqrt2*FV_cut/f_pi*(Mt2*FF_1/u_minus_M2 + FF_2/M);
251 V_2 = g_A*kSqrt2*FV_cut/f_pi*Mt2*FF_1/u_minus_M2/qk;
252 V_3 = g_A*kSqrt2*FV_cut/f_pi*2*FF_2/u_minus_M2;
253 V_4 = -V_3;
254 V_5 = g_A*kSqrt2*FV_cut/f_pi*M*FF_pi/t_minus_mpi2/qk;
255 break;
256
257 case (kSpp_vn_cc_01100) :
258 case (kSpp_vbp_cc_10001):
259 V_1 = g_A*kSqrt2*FV_cut/f_pi*(Mt2*FF_1/s_minus_M2 + FF_2/M);
260 V_2 = g_A*kSqrt2*FV_cut/f_pi*Mt2*FF_1/s_minus_M2/qk;
261 V_3 = -g_A*kSqrt2*FV_cut/f_pi*2*FF_2/s_minus_M2;
262 V_4 = V_3;
263 V_5 = -g_A*kSqrt2*FV_cut/f_pi*M*FF_pi/t_minus_mpi2/qk;
264 break;
265
266 case (kSpp_vn_cc_10010) :
267 case (kSpp_vbp_cc_01010):
268 V_1 = g_A*FV_cut/f_pi*Mt2*FF_1*(1/s_minus_M2 - 1/u_minus_M2);
269 V_2 = g_A*FV_cut/f_pi*Mt2*FF_1/qk*(1/s_minus_M2 - 1/u_minus_M2);
270 V_3 = -g_A*FV_cut/f_pi*2*FF_2*(1/s_minus_M2 + 1/u_minus_M2);
271 V_4 = -g_A*FV_cut/f_pi*2*FF_2*(1/s_minus_M2 - 1/u_minus_M2);;
272 V_5 = -g_A*FV_cut/f_pi*Mt2*FF_pi/t_minus_mpi2/qk;
273 break;
274
275 //NC
276 case (kSpp_vp_nc_10010) :
277 case (kSpp_vbp_nc_10010):
278 case (kSpp_vn_nc_01010) :
279 case (kSpp_vbn_nc_01010):
280 V_1 = g_A*FV_cut*M/f_pi*(FF_1*(1/u_minus_M2 + 1/s_minus_M2) + FF_2/M2);
281 V_2 = g_A*FV_cut*M/f_pi*FF_1*(1/u_minus_M2 + 1/s_minus_M2)/qk;
282 V_3 = g_A*FV_cut/2/f_pi*FF_2*(1/u_minus_M2 - 1/s_minus_M2);
283 V_4 = -g_A*FV_cut/2/f_pi*FF_2*(1/u_minus_M2 + 1/s_minus_M2);
284 break;
285
286 case (kSpp_vp_nc_01100) :
287 case (kSpp_vbp_nc_01100):
288 V_1 = -k1_Sqrt2*g_A*FV_cut*Mt2/f_pi*FF_1*(1/u_minus_M2 - 1/s_minus_M2);
289 V_2 = -k1_Sqrt2*g_A*FV_cut*Mt2/f_pi*FF_1*(1/u_minus_M2 - 1/s_minus_M2)/qk;
290 V_3 = -k1_Sqrt2*g_A*FV_cut/f_pi*FF_2*(1/u_minus_M2 + 1/s_minus_M2);
291 V_4 = k1_Sqrt2*g_A*FV_cut/f_pi*FF_2*(1/u_minus_M2 - 1/s_minus_M2);
292 V_5 = -k1_Sqrt2*g_A*FV_cut*Mt2/f_pi*FF_pi/t_minus_mpi2/qk;
293 break;
294
295 case (kSpp_vn_nc_10001) :
296 case (kSpp_vbn_nc_10001):
297 V_1 = k1_Sqrt2*g_A*FV_cut*Mt2/f_pi*FF_1*(1/u_minus_M2 - 1/s_minus_M2);
298 V_2 = k1_Sqrt2*g_A*FV_cut*Mt2/f_pi*FF_1*(1/u_minus_M2 - 1/s_minus_M2)/qk;
299 V_3 = k1_Sqrt2*g_A*FV_cut/f_pi*FF_2*(1/u_minus_M2 + 1/s_minus_M2);
300 V_4 = -k1_Sqrt2*g_A*FV_cut/f_pi*FF_2*(1/u_minus_M2 - 1/s_minus_M2);
301 V_5 = k1_Sqrt2*g_A*FV_cut*Mt2/f_pi*FF_pi/t_minus_mpi2/qk;
302 break;
303 default:
304 // should not be here - meaningless to return anything
305 gAbortingInErr = true;
306 LOG("MKSPPPXSec2020", pFATAL) << "Unknown resonance production channel: " << spp_channel;
307 exit(1);
308
309
310 }
311 double V_25 = (W_plus*W_minus)*V_2 - Q2*V_5;
312
313
314 double O_1_plus = TMath::Sqrt((W_plus2 + Q2)*( W_plus2 - m_pi2 ))/Wt2;
315 // There is no check of positivity under root expression in the original code
316 double O_1_minus = TMath::Sqrt(TMath::Max(0., (W_minus2 + Q2)*(W_minus2 - m_pi2)))/Wt2;
317 double O_2_plus = TMath::Sqrt((W_plus2 + Q2)/(W_plus2 - m_pi2));
318 // There is no check of positivity under root expression in the original code
319 double O_2_minus = W_minus2>m_pi2?TMath::Sqrt((W_minus2 + Q2)/(W_minus2 - m_pi2)):0;
320
321 double K_1_V = W_minus*O_1_plus;
322 double K_2_V = W_plus*O_1_minus;
323 // There is no check of positivity under root expression in the original code
324 double K_3_V = W_minus2>m_pi2?abs_mom_q*abs_mom_q*W_plus*O_2_minus:0;
325 double K_4_V = abs_mom_q*abs_mom_q*W_minus*O_2_plus;
326 double K_5_V = 1/O_2_plus;
327 // the of K_6_V = 1/O_2_minus differ from original code
328 double K_6_V = TMath::Sqrt(TMath::Max(0., (W_minus2 - m_pi2))/(W_minus2 + Q2));
329
330 double F_1 = V_1 + qk*(V_3 - V_4)/W_minus + W_minus*V_4;
331 double F_2 = -V_1 + qk*(V_3 - V_4)/W_plus + W_plus *V_4;
332 double F_3 = (V_3 - V_4) + V_25/W_plus;
333 double F_4 = (V_3 - V_4) - V_25/W_minus;
334 double F_5 = (W_plus2 + Q2)*(V_1 + W_minus*V_4)/Wt2 + (W_plus*q_0 - qk)*(V_3 - V_4) + q_0*V_25 - k_0*qk*V_5 -
335 qk*(W_plus2 + Q2 + 2*W*W_minus)*V_2/Wt2;
336 double F_6 =-(W_minus2 + Q2)*(V_1 - W_plus *V_4)/Wt2 + (W_minus*q_0 - qk)*(V_3 - V_4) - q_0*V_25 + k_0*qk*V_5 +
337 qk*(W_plus2 + Q2 + 2*W*W_minus)*V_2/Wt2;
338
339 double sF_1 = K_1_V*F_1/Mt2;
340 double sF_2 = K_2_V*F_2/Mt2;
341 double sF_3 = K_3_V*F_3/Mt2;
342 double sF_4 = K_4_V*F_4/Mt2;
343 double sF_5 = K_5_V*F_5/Mt2;
344 double sF_6 = K_6_V*F_6/Mt2;
345
347 k1_Sqrt2*SinTheta*CosHalfTheta*(sF_3 + sF_4)*std::complex<double>
350
352 -k1_Sqrt2*SinTheta*SinHalfTheta*(sF_3 - sF_4)*std::complex<double>
355
357 kSqrt2*(CosHalfTheta*(sF_1 - sF_2) - 0.5*SinTheta*SinHalfTheta*(sF_3 - sF_4))*std::complex<double>
360
362 -kSqrt2*(SinHalfTheta*(sF_1 + sF_2) + 0.5*SinTheta*CosHalfTheta*(sF_3 + sF_4))*std::complex<double>
365
366
368 -kSqrt2*(SinHalfTheta*(sF_1 + sF_2) + 0.5*SinTheta*CosHalfTheta*(sF_3 + sF_4))*std::complex<double>
371
373 -kSqrt2*(CosHalfTheta*(sF_1 - sF_2) - 0.5*SinTheta*SinHalfTheta*(sF_3 - sF_4))*std::complex<double>
376
378 k1_Sqrt2*SinTheta*SinHalfTheta*(sF_3 - sF_4)*std::complex<double>
381
383 k1_Sqrt2*SinTheta*CosHalfTheta*(sF_3 + sF_4)*std::complex<double>
386
387
388 if (is_CC || (is_NC && is_nu) )
389 {
391 CosHalfTheta*(k_0*eps_z_L - abs_mom_k*eps_zero_L)*(sF_5 + sF_6)/C_S_minus*std::complex<double>
394
396 -SinHalfTheta*(k_0*eps_z_L - abs_mom_k*eps_zero_L)*(sF_5 - sF_6)/C_S_minus*std::complex<double>
399
401 -SinHalfTheta*(k_0*eps_z_L - abs_mom_k*eps_zero_L)*(sF_5 - sF_6)/C_S_minus*std::complex<double>
404
406 -CosHalfTheta*(k_0*eps_z_L - abs_mom_k*eps_zero_L)*(sF_5 + sF_6)/C_S_minus*std::complex<double>
409 }
410 if (is_CC || (is_NC && is_nubar) )
411 {
413 CosHalfTheta*(k_0*eps_z_R - abs_mom_k*eps_zero_R)*(sF_5 + sF_6)/C_S_plus*std::complex<double>
416
418 -SinHalfTheta*(k_0*eps_z_R - abs_mom_k*eps_zero_R)*(sF_5 - sF_6)/C_S_plus*std::complex<double>
421
423 -SinHalfTheta*(k_0*eps_z_R - abs_mom_k*eps_zero_R)*(sF_5 - sF_6)/C_S_plus*std::complex<double>
426
428 -CosHalfTheta*(k_0*eps_z_R - abs_mom_k*eps_zero_R)*(sF_5 + sF_6)/C_S_plus*std::complex<double>
431 }
432
433 if (is_NC)
434 {
435 // isoscalar electromagnetic amplitudes eq.68 of ref. 8
436 fEMFormFactors.Calculate(interaction);
437 double FF_em_1 = 0.5*fEMFormFactors.F1V();
438 double FF_em_2 = 0.5*fEMFormFactors.xiF2V();
439
440 double Vem_1 = 0, Vem_2 = 0, Vem_3 = 0, Vem_4 = 0, Vem_5 = 0;
441
442 // [p,n,pi+,pi0,pi-]
443 switch (spp_channel) {
444 //NC
445 case (kSpp_vp_nc_10010) :
446 case (kSpp_vbp_nc_10010):
447 Vem_1 = -k1_Sqrt2*g_A*FV_cut*M/f_pi*(FF_em_1*(1/u_minus_M2 + 1/s_minus_M2) + FF_em_2/M2);
448 Vem_2 = -k1_Sqrt2*g_A*FV_cut*M/f_pi*FF_em_1*(1/u_minus_M2 + 1/s_minus_M2)/qk;
449 Vem_3 = -k1_Sqrt2*g_A*FV_cut/2/f_pi*FF_em_2*(1/u_minus_M2 - 1/s_minus_M2);
450 Vem_4 = k1_Sqrt2*g_A*FV_cut/2/f_pi*FF_em_2*(1/u_minus_M2 + 1/s_minus_M2);
451 break;
452 case (kSpp_vn_nc_01010) :
453 case (kSpp_vbn_nc_01010):
454 Vem_1 = k1_Sqrt2*g_A*FV_cut*M/f_pi*(FF_em_1*(1/u_minus_M2 + 1/s_minus_M2) + FF_em_2/M2);
455 Vem_2 = k1_Sqrt2*g_A*FV_cut*M/f_pi*FF_em_1*(1/u_minus_M2 + 1/s_minus_M2)/qk;
456 Vem_3 = k1_Sqrt2*g_A*FV_cut/2/f_pi*FF_em_2*(1/u_minus_M2 - 1/s_minus_M2);
457 Vem_4 =-k1_Sqrt2*g_A*FV_cut/2/f_pi*FF_em_2*(1/u_minus_M2 + 1/s_minus_M2);
458 break;
459
460 case (kSpp_vp_nc_01100) :
461 case (kSpp_vbp_nc_01100):
462 case (kSpp_vn_nc_10001) :
463 case (kSpp_vbn_nc_10001):
464 Vem_1 = g_A*FV_cut*M/f_pi*(FF_em_1*(1/u_minus_M2 + 1/s_minus_M2) + FF_em_2/M2/2);
465 Vem_2 = g_A*FV_cut*M/f_pi*FF_em_1*(1/u_minus_M2 + 1/s_minus_M2)/qk;
466 Vem_3 = g_A*FV_cut/2/f_pi*FF_em_2*(1/u_minus_M2 - 1/s_minus_M2);
467 Vem_4 =-g_A*FV_cut/2/f_pi*FF_em_2*(1/u_minus_M2 + 1/s_minus_M2);
468 break;
469 default:
470 // should not be here - meaningless to return anything
471 gAbortingInErr = true;
472 LOG("MKSPPPXSec2020", pFATAL) << "CC channel in NC mode " << spp_channel;
473 exit(1);
474 }
475 double Vem_25 = (W_plus*W_minus)*Vem_2 + Vem_5;
476
477 double Fem_1 = Vem_1 + qk*(Vem_3 - Vem_4)/W_minus + W_minus*Vem_4;
478 double Fem_2 = -Vem_1 + qk*(Vem_3 - Vem_4)/W_plus + W_plus*Vem_4;
479 double Fem_3 = (Vem_3 - Vem_4) + Vem_25/W_plus;
480 double Fem_4 = (Vem_3 - Vem_4) - Vem_25/W_minus;
481 double Fem_5 = (W_plus2 + Q2)*(Vem_1 + W_minus*Vem_4)/Wt2 + (W_plus*q_0 - qk)*(Vem_3 - Vem_4) + q_0*Vem_25 -
482 qk*(W_plus2 + Q2 + 2*W*W_minus)*Vem_2/Wt2;
483 double Fem_6 =-(W_minus2 + Q2)*(Vem_1 - W_plus*Vem_4)/Wt2 + (W_minus*q_0 - qk)*(Vem_3 - Vem_4) - q_0*Vem_25 +
484 qk*(W_plus2 + Q2 + 2*W*W_minus)*Vem_2/Wt2;
485
486 double sFem_1 = K_1_V*Fem_1/Mt2;
487 double sFem_2 = K_2_V*Fem_2/Mt2;
488 double sFem_3 = K_3_V*Fem_3/Mt2;
489 double sFem_4 = K_4_V*Fem_4/Mt2;
490 double sFem_5 = K_5_V*Fem_5/Mt2;
491 double sFem_6 = K_6_V*Fem_6/Mt2;
492
494
496 k1_Sqrt2*SinTheta*CosHalfTheta*(sFem_3 + sFem_4)*std::complex<double>
499
501 -k1_Sqrt2*SinTheta*SinHalfTheta*(sFem_3 - sFem_4)*std::complex<double>
504
506 kSqrt2*(CosHalfTheta*(sFem_1 - sFem_2) - 0.5*SinTheta*SinHalfTheta*(sFem_3 - sFem_4))*std::complex<double>
509
511 -kSqrt2*(SinHalfTheta*(sFem_1 + sFem_2) + 0.5*SinTheta*CosHalfTheta*(sFem_3 + sFem_4))*std::complex<double>
514
515
517 -kSqrt2*(SinHalfTheta*(sFem_1 + sFem_2) + 0.5*SinTheta*CosHalfTheta*(sFem_3 + sFem_4))*std::complex<double>
520
522 -kSqrt2*(CosHalfTheta*(sFem_1 - sFem_2) - 0.5*SinTheta*SinHalfTheta*(sFem_3 - sFem_4))*std::complex<double>
525
527 k1_Sqrt2*SinTheta*SinHalfTheta*(sFem_3 - sFem_4)*std::complex<double>
530
532 k1_Sqrt2*SinTheta*CosHalfTheta*(sFem_3 + sFem_4)*std::complex<double>
535
536
537 if (is_nu)
538 {
540 CosHalfTheta*(k_0*eps_z_L - abs_mom_k*eps_zero_L)*(sFem_5 + sFem_6)/C_S_minus*std::complex<double>
543
545 -SinHalfTheta*(k_0*eps_z_L - abs_mom_k*eps_zero_L)*(sFem_5 - sFem_6)/C_S_minus*std::complex<double>
548
550 -SinHalfTheta*(k_0*eps_z_L - abs_mom_k*eps_zero_L)*(sFem_5 - sFem_6)/C_S_minus*std::complex<double>
553
555 -CosHalfTheta*(k_0*eps_z_L - abs_mom_k*eps_zero_L)*(sFem_5 + sFem_6)/C_S_minus*std::complex<double>
558 }
559 else
560 {
562 CosHalfTheta*(k_0*eps_z_R - abs_mom_k*eps_zero_R)*(sFem_5 + sFem_6)/C_S_plus*std::complex<double>
565
567 -SinHalfTheta*(k_0*eps_z_R - abs_mom_k*eps_zero_R)*(sFem_5 - sFem_6)/C_S_plus*std::complex<double>
570
572 -SinHalfTheta*(k_0*eps_z_R - abs_mom_k*eps_zero_R)*(sFem_5 - sFem_6)/C_S_plus*std::complex<double>
575
577 -CosHalfTheta*(k_0*eps_z_R - abs_mom_k*eps_zero_R)*(sFem_5 + sFem_6)/C_S_plus*std::complex<double>
580 }
581
582 Hbkg = (1 - 2*fSin2Wein)*Hbkg - 4*fSin2Wein*HbkgEM;
583 }
584
585 //Axial helicity amp for bkg
586 //Axial cut is same as FV_cut in the latest version
587 double FA_cut = FV_cut;
588
589 // FF_A here is equal to -g_A*FF_A of original code
590 double FF_A = fFormFactors.FA();
591 double F_rho = Frho0/(1 - (t_minus_mpi2 + m_pi2 )/fRho770Mass/fRho770Mass);
592
593 // Eq. 38 and Table 5 of ref. 1
594 double A_1 = 0, A_3 = 0, A_4 = 0, A_7 = 0, A_8 = 0;
595
596 // [p,n,pi+,pi0,pi-]
597 switch (spp_channel) {
598
599 //CC
600 case (kSpp_vp_cc_10100) :
601 case (kSpp_vbn_cc_01001):
602 A_1 = -kSqrt2*FA_cut*M*FF_A/f_pi/u_minus_M2;
603 A_3 = -A_1;
604 A_4 = k1_Sqrt2*FA_cut/f_pi*(FF_A + F_rho)/M;
605 A_7 = kSqrt2*FA_cut/f_pi*M*FF_A/mpi2_minus_k2;
606 A_8 = -k1_Sqrt2*FA_cut/f_pi*(FF_A*(1 + 4*M2/u_minus_M2) + F_rho)/mpi2_minus_k2;
607 break;
608
609 case (kSpp_vn_cc_01100) :
610 case (kSpp_vbp_cc_10001):
611 A_1 = kSqrt2*FA_cut/f_pi*M*FF_A/s_minus_M2;
612 A_3 = A_1;
613 A_4 = -k1_Sqrt2*FA_cut/f_pi/M*(FF_A + F_rho);
614 A_7 = kSqrt2*FA_cut/f_pi*M*FF_A/mpi2_minus_k2;
615 A_8 = k1_Sqrt2/f_pi*FA_cut*(FF_A/mpi2_minus_k2*(1 + 4*M2/s_minus_M2) + F_rho/mpi2_minus_k2);
616 break;
617
618 case (kSpp_vn_cc_10010) :
619 case (kSpp_vbp_cc_01010):
620 A_1 = FA_cut/f_pi*M*FF_A*(1/u_minus_M2 + 1/s_minus_M2);
621 A_3 = -FA_cut/f_pi*M*FF_A*(1/u_minus_M2 - 1/s_minus_M2) ;
622 A_4 = -FA_cut/f_pi/M*(FF_A + F_rho);
623 A_8 = FA_cut/f_pi*(F_rho + FF_A*(1 + 2*M2*(1/u_minus_M2 + 1/s_minus_M2)))/mpi2_minus_k2;
624 break;
625
626 //NC
627 case (kSpp_vp_nc_10010) :
628 case (kSpp_vbp_nc_10010):
629 case (kSpp_vn_nc_01010) :
630 case (kSpp_vbn_nc_01010):
631 A_1 = -FA_cut/f_pi*M*FF_A*(1/u_minus_M2 - 1/s_minus_M2)/2;
632 A_3 = FA_cut/f_pi*M*FF_A*(1/u_minus_M2 + 1/s_minus_M2)/2;
633 A_7 = FA_cut/f_pi*M*FF_A/mpi2_minus_k2;
634 break;
635
636 case (kSpp_vp_nc_01100) :
637 case (kSpp_vbp_nc_01100):
638 A_1 = k1_Sqrt2*FA_cut/f_pi*M*FF_A*(1/u_minus_M2 + 1/s_minus_M2);
639 A_3 = -k1_Sqrt2*FA_cut/f_pi*M*FF_A*(1/u_minus_M2 - 1/s_minus_M2);
640 A_4 = -k1_Sqrt2*FA_cut/f_pi/M*(FF_A + F_rho);
641 A_8 = k1_Sqrt2*FA_cut/f_pi/mpi2_minus_k2*(F_rho + FF_A*(1 + 2*M2*(1/u_minus_M2 + 1/s_minus_M2)));
642 break;
643
644 case (kSpp_vn_nc_10001) :
645 case (kSpp_vbn_nc_10001):
646 A_1 = -k1_Sqrt2*FA_cut/f_pi*M*FF_A*(1/u_minus_M2 + 1/s_minus_M2);
647 A_3 = k1_Sqrt2*FA_cut/f_pi*M*FF_A*(1/u_minus_M2 - 1/s_minus_M2);
648 A_4 = k1_Sqrt2*FA_cut/f_pi/M*(FF_A + F_rho);
649 A_8 = k1_Sqrt2*FA_cut/f_pi/mpi2_minus_k2*(FF_A*(1 + 2*M2*(1/u_minus_M2 + 1/s_minus_M2)) + F_rho);
650 break;
651 default:
652 // should not be here - meaningless to return anything
653 gAbortingInErr = true;
654 LOG("MKSPPPXSec2020", pFATAL) << "Unknown resonance production channel: " << spp_channel;
655 exit(1);
656 }
657
658 double K_1_A = abs_mom_q*O_2_plus;
659 // There is no check of positivity under root expression in the original code
660 double K_2_A = W_minus2>m_pi2?abs_mom_q*O_2_minus:1;
661 double K_3_A = abs_mom_q*O_1_minus;
662 double K_4_A = abs_mom_q* O_1_plus;
663 double K_5_A = O_1_minus;
664 double K_6_A = O_1_plus;
665 double abs_mom_k2 = abs_mom_k*abs_mom_k;
666 double Delta = k_0*(q_0*k_0 - qk)/abs_mom_k2;
667
668 double G_1 = W_plus* A_1 - M*A_4;
669 double G_2 = -W_minus*A_1 - M*A_4;
670 double G_3 = A_1 - A_3;
671 double G_4 = -G_3;
672 double G_5 = (Delta + (W_plus2 - m_pi2)/Wt2 + Wt2*k_0*W_plus/(W_minus2 + Q2))*A_1 + (q_0 - Delta)*A_3 - M*W_minus*A_4/(p_10 - M);
673 double G_6 =-(Delta + (W_minus2 - m_pi2)/Wt2 + Wt2*k_0*W_minus/(W_plus2 + Q2))*A_1 - (q_0 - Delta)*A_3 - M*W_plus*A_4 /(p_10 + M);
674 double G_7= (W_plus2 - m_pi2)*A_1/Wt2 + q_0*A_3 - M*A_4 + k_0*A_7 + W_plus *k_0*A_8;
675 double G_8 = -(W_minus2 - m_pi2)*A_1/Wt2 - q_0*A_3 - M*A_4 - k_0*A_7 + W_minus*k_0*A_8;
676
677 //Isobaric amplitud (Scribal G)
678 double sG_1 = K_1_A*G_1/Mt2;
679 double sG_2 = K_2_A*G_2/Mt2;
680 double sG_3 = K_3_A*G_3/Mt2;
681 double sG_4 = K_4_A*G_4/Mt2;
682 double sG_5 = K_5_A*G_5/Mt2;
683 double sG_6 = K_6_A*G_6/Mt2;
684 double sG_7 = K_5_A*G_7/Mt2;
685 double sG_8 = K_6_A*G_8/Mt2;
686
687 // scalar part eq. 71 of ref.8
688 if (is_NC)
689 {
690 double g_s = 0.15;
691 if (spp_channel == kSpp_vp_nc_01100 || spp_channel == kSpp_vbp_nc_01100 || spp_channel == kSpp_vn_nc_10001 || spp_channel == kSpp_vbn_nc_10001)
692 g_s *= kSqrt2;
693
694 double Gs_A = g_s/g_A*FF_A;
695 double As_1 = k1_Sqrt2*g_A*FA_cut/f_pi*M*Gs_A*(1/u_minus_M2 - 1/s_minus_M2);
696 double As_3 =-k1_Sqrt2*g_A*FA_cut/f_pi*M*Gs_A*(1/u_minus_M2 + 1/s_minus_M2);
697 double As_4 = 0;
698 double As_7 = 0;
699 double As_8 = 0;
700 if (spp_channel == kSpp_vn_nc_01010 || spp_channel == kSpp_vbn_nc_01010)
701 {
702 As_1 *= -1;
703 As_3 *= -1;
704 }
705
706 double Gs_1 = W_plus* As_1 - M*As_4;
707 double Gs_2 = -W_minus*As_1 - M*As_4;
708 double Gs_3 = As_1 - As_3;
709 double Gs_4 = -Gs_3;
710 double Gs_5 = (Delta + (W_plus2 - m_pi2)/Wt2 + 2*W*k_0*W_plus /(W_minus2 + Q2))*As_1 + (q_0 - Delta)*As_3 - M*W_minus*As_4/(p_10-M);
711 double Gs_6 =-(Delta + (W_minus2 - m_pi2)/Wt2 + 2*W*k_0*W_minus/(W_plus2 + Q2)) *As_1 - (q_0 - Delta)*As_3 - M*W_plus*As_4/(p_10 + M);
712 double Gs_7 = (W_plus2 - m_pi2)*As_1/Wt2 + q_0*As_3 - M*As_4 + k_0*As_7 + W_plus* k_0*As_8;
713 double Gs_8 =-(W_minus2 - m_pi2)*As_1/Wt2 - q_0*As_3 - M*As_4 - k_0*As_7 + W_minus* k_0*As_8;
714
715 sG_1 -= K_1_A*Gs_1/Mt2;
716 sG_2 -= K_2_A*Gs_2/Mt2;
717 sG_3 -= K_3_A*Gs_3/Mt2;
718 sG_4 -= K_4_A*Gs_4/Mt2;
719 sG_5 -= K_5_A*Gs_5/Mt2;
720 sG_6 -= K_6_A*Gs_6/Mt2;
721 sG_7 -= K_5_A*Gs_7/Mt2;
722 sG_8 -= K_6_A*Gs_8/Mt2;
723 }
724
725 // helicity amplitudes
727 k1_Sqrt2*SinTheta*CosHalfTheta*(sG_3 + sG_4)*std::complex<double>
730
732 k1_Sqrt2*SinTheta*SinHalfTheta*(sG_3 - sG_4)*std::complex<double>
735
737 kSqrt2*(CosHalfTheta*(sG_1 - sG_2) - 0.5*SinTheta*SinHalfTheta*(sG_3 - sG_4))*std::complex<double>
740
742 kSqrt2*(SinHalfTheta*(sG_1 + sG_2) + 0.5*SinTheta*CosHalfTheta*(sG_3 + sG_4))*std::complex<double>
745
746
748 -kSqrt2*(SinHalfTheta*(sG_1 + sG_2) + 0.5*SinTheta*CosHalfTheta*(sG_3 + sG_4))*std::complex<double>
751
753 kSqrt2*(CosHalfTheta*(sG_1 - sG_2) - 0.5*SinTheta*SinHalfTheta*(sG_3 - sG_4))*std::complex<double>
756
758 k1_Sqrt2*SinTheta*SinHalfTheta*(sG_3 - sG_4)*std::complex<double>
761
763 -k1_Sqrt2*SinTheta*CosHalfTheta*(sG_3 + sG_4)*std::complex<double>
766
767 if (is_CC || (is_NC && is_nu) )
768 {
770 CosHalfTheta*(abs_mom_k*eps_z_L*(sG_5 + sG_6) + (k_0*eps_zero_L - abs_mom_k*eps_z_L)*(sG_7 + sG_8))/k_0/C_S_minus*std::complex<double>
773
775 SinHalfTheta*(abs_mom_k*eps_z_L*(sG_5 - sG_6) + (k_0*eps_zero_L - abs_mom_k*eps_z_L)*(sG_7 - sG_8))/k_0/C_S_minus*std::complex<double>
778
780 -SinHalfTheta*(abs_mom_k*eps_z_L*(sG_5 - sG_6) + (k_0*eps_zero_L - abs_mom_k*eps_z_L)*(sG_7 - sG_8))/k_0/C_S_minus*std::complex<double>
783
785 CosHalfTheta*(abs_mom_k*eps_z_L*(sG_5 + sG_6) + (k_0*eps_zero_L - abs_mom_k*eps_z_L)*(sG_7 + sG_8))/k_0/C_S_minus*std::complex<double>
788 }
789 if (is_CC || (is_NC && is_nubar) )
790 {
792 CosHalfTheta*(abs_mom_k*eps_z_R*(sG_5 + sG_6) + (k_0*eps_zero_R - abs_mom_k*eps_z_R)*(sG_7 + sG_8))/k_0/C_S_plus*std::complex<double>
795
797 SinHalfTheta*(abs_mom_k*eps_z_R*(sG_5 - sG_6) + (k_0*eps_zero_R - abs_mom_k*eps_z_R)*(sG_7 - sG_8))/k_0/C_S_plus*std::complex<double>
800
802 -SinHalfTheta*(abs_mom_k*eps_z_R*(sG_5 - sG_6) + (k_0*eps_zero_R - abs_mom_k*eps_z_R)*(sG_7 - sG_8))/k_0/C_S_plus*std::complex<double>
805
807 CosHalfTheta*(abs_mom_k*eps_z_R*(sG_5 + sG_6) + (k_0*eps_zero_R - abs_mom_k*eps_z_R)*(sG_7 + sG_8))/k_0/C_S_plus*std::complex<double>
810 }
811
812 /* Resonace contribution */
814 //f_BW_A is same as f_BW_V in the latest version, so rename f_BW_V = f_BW_A -> f_BW
815 std::complex<double> kappa_f_BW;
816
817 for (auto res : fResList)
818 {
819
820 if (utils::res::Isospin(res) == 1 && SppChannel::FinStateIsospin(spp_channel) == 3) // skip resonances with I=1/2 if isospin of final state is 3/2
821 continue;
822
823 int NR = utils::res::ResonanceIndex (res);
824 int LR = utils::res::OrbitalAngularMom (res);
825 int JR = (utils::res::AngularMom (res) + 1)/2;
826 double MR = utils::res::Mass (res);
827 double WR = utils::res::Width (res);
828 double Cjsgn_plus = utils::res::Cjsgn_plus (res);
829 double Dsgn = utils::res::Dsgn (res);
830 double BR = SppChannel::BranchingRatio (res);
831
832
833 double d = W_plus2 + Q2;
834 double sq2omg = TMath::Sqrt(2/fOmega);
835 double nomg = NR*fOmega;
836
837 //Graczyk and Sobczyk vector form-factors
838 double CV_factor = 1/(1 + Q2/fMv2/4);
839 // Eq. 29 of ref. 6
840 double CV3 = fCv3*CV_factor/(1 + Q2/fMv2)/(1 + Q2/fMv2);
841 // Eq. 30 of ref. 6
842 double CV4 = -1. * fCv4 / fCv3 * CV3;
843 // Eq. 31 of ref. 6
844 double CV5 = fCv51*CV_factor/(1 + Q2/fMv2/fCv52)/(1 + Q2/fMv2/fCv52);
845
846 // Eq. 38 of ref. 6
847 double GV3 = 0.5*k1_Sqrt3*(CV4*(W2 - M2 - Q2)/2/M2 + CV5*(W2 - M2 + Q2)/2/M2 + CV3*W_plus/M);
848 // Eq. 39 of ref. 6
849 double GV1 = -0.5*k1_Sqrt3*(CV4*(W2 - M2 - Q2)/2/M2 + CV5*(W2 - M2 + Q2)/2/M2 - CV3*(W_plus*M + Q2)/W/M);
850 // Eq. 36 of ref. 6, which is implied to use for EM-production
851 double GV = 0.5*TMath::Sqrt(1 + Q2/W_plus2)/TMath::Power(1 + Q2/4/M2, 0.5*NR)*TMath::Sqrt(3*GV3*GV3 + GV1*GV1);
852 // Eq. 37 of ref. 6, which is implied to use for neutrino-production
853 // double GV = 0.5*TMath::Sqrt(1 + Q2/W_plus2)/TMath::Power(1 + Q2/4/M2, NR)*TMath::Sqrt(3*GV3*GV3 + GV1*GV1);
854
855
856 //Graczyk and Sobczyk axial form-symmetry_factor
857 // Eq. 52 of ref. 6
858 double CA5 = fCA50/(1 + Q2/fMa2)/(1 + Q2/fMa2);
859
860 // The form is the same like in Eq. 54 of ref. 6, but differ from it by index, which in ref. 6 is equal to NR.
861 double GA = 0.5*kSqrt3*TMath::Sqrt(1 + Q2/W_plus2)*(1 - (W2 - Q2 -M2)/8/M2)*CA5/TMath::Power(1+ Q2/4/M2, 0.5*NR);
862
863 double qMR_0 = (MR*MR - M2 + m_pi2)/(2*MR);
864 double abs_mom_qMR = TMath::Sqrt( qMR_0*qMR_0 - m_pi2);
865 double Gamma = WR*TMath::Power((abs_mom_q/abs_mom_qMR), 2*LR + 1);
866
867 // denominator of Breit-Wigner function
868 std::complex<double> denom(W - MR, Gamma/2);
869 // Breit-Wigner amplitude multiplied by kappa*sqrt(BR) to avoid singularity at abs_mom_q=0, where BR = chi_E (see eq. 25 and 27 of ref. 1)
870 // double kappa = kPi*W*TMath::Sqrt(2/JR/abs_mom_q)/M;
871 // f_BW = TMath::Sqrt(BR*Gamma/2/kPi)/denom;
872 kappa_f_BW = W*TMath::Sqrt(kPi*BR*WR/JR/abs_mom_qMR)*TMath::Power((abs_mom_q/abs_mom_qMR), LR)/denom/M;
873
874 fFKR.Lamda = sq2omg*abs_mom_k;
875 fFKR.Tv = GV/3/W/sq2omg;
876 fFKR.Ta = 2./3/sq2omg*abs_mom_k*GA/d;
877 fFKR.Rv = kSqrt2*abs_mom_k*W_plus*GV/d;
878 fFKR.Ra = kSqrt2/6*(W_plus + 2*nomg*W/d)*GA/W;
879 fFKR.R = fFKR.Rv;
880 fFKR.T = fFKR.Tv;
881 fFKR.Rplus = - (fFKR.Rv + fFKR.Ra);
882 fFKR.Rminus = - (fFKR.Rv - fFKR.Ra);
883 fFKR.Tplus = - (fFKR.Tv + fFKR.Ta);
884 fFKR.Tminus = - (fFKR.Tv - fFKR.Ta);
885
886 double a_aux = 1 + ((W2 + Q2 + M2)/(Mt2*W));
887 // if Q2=Q=sqrt(Q2)=0 then the singularity appears in k_sqrtQ2
888 // which is eliminated when in Hres at BosonPolarization=PLUS0 or BosonPolarization=MINUS0
889 // when k_sqrtQ2 is multiplied by C, B and S. Therefore in this case we put Q equal to 1.
890 // In our opinion it is simplier to eliminate singularity in such way, because it allows to
891 // keep intact original formulas from paper.
892 if (Q == 0) Q = 1;
893
894 double C_plus = Q/C_S_plus*((eps_zero_R*abs_mom_k - eps_z_R*k_0)*(1./3 + k_0/a_aux/M) +
895 (Wt2/3 - Q2/a_aux/M + nomg/a_aux/M/3)*(eps_z_R + (eps_zero_R*k_0 - eps_z_R*abs_mom_k)*abs_mom_k/mpi2_minus_k2))*GA/Wt2/abs_mom_k;
896 double C_minus = Q/C_S_minus*((eps_zero_L*abs_mom_k - eps_z_L*k_0)*(1./3 + k_0/a_aux/M) +
897 (Wt2/3 - Q2/a_aux/M + nomg/a_aux/M/3)*(eps_z_L + (eps_zero_L*k_0 - eps_z_L*abs_mom_k)*abs_mom_k/mpi2_minus_k2))*GA/Wt2/abs_mom_k;
898
899 double B_plus = Q/C_S_plus *(eps_zero_R + eps_z_R*abs_mom_k/a_aux/M +
900 (eps_zero_R*k_0 - eps_z_R*abs_mom_k)*(k_0 + abs_mom_k*abs_mom_k/M/a_aux)/mpi2_minus_k2)*GA/W/3/sq2omg/abs_mom_k;
901 double B_minus = Q/C_S_minus*(eps_zero_L + eps_z_L*abs_mom_k/a_aux/M +
902 (eps_zero_L*k_0 - eps_z_L*abs_mom_k)*(k_0 + abs_mom_k*abs_mom_k/M/a_aux)/mpi2_minus_k2)*GA/W/3/sq2omg/abs_mom_k;
903
904 double S_plus = Q/C_S_plus *(eps_z_R*k_0 - eps_zero_R*abs_mom_k)*(1 + Q2/M2 - 3*W/M)*GV/abs_mom_k_L2/6;
905 double S_minus = Q/C_S_minus*(eps_z_L*k_0 - eps_zero_L*abs_mom_k)*(1 + Q2/M2 - 3*W/M)*GV/abs_mom_k_L2/6;
906 double k_sqrtQ2 = abs_mom_k/Q;
907
908 const RSHelicityAmplModelI * hamplmod = 0;
909
910 if (is_CC)
911 hamplmod = fHAmplModelCC;
912 else if(is_NC)
913 {
914 if (is_p)
915 hamplmod = fHAmplModelNCp;
916 else
917 hamplmod = fHAmplModelNCn;
918 }
919
920 fFKR.S = S_minus;
921 fFKR.B = B_minus;
922 fFKR.C = C_minus;
923 const RSHelicityAmpl & hampl = hamplmod->Compute(res, fFKR);
924 double fp3 = hampl.AmpPlus3();
925 double fp1 = hampl.AmpPlus1();
926 double fm3 = hampl.AmpMinus3();
927 double fm1 = hampl.AmpMinus1();
928 double fm0m = 0, fm0p = 0, fp0m = 0, fp0p = 0;
929 if (is_CC || (is_NC && is_nu) )
930 {
931 fm0m = hampl.Amp0Minus();
932 fm0p = hampl.Amp0Plus();
933 }
934 if (is_CC || (is_NC && is_nubar) )
935 {
936 fFKR.S = S_plus;
937 fFKR.B = B_plus;
938 fFKR.C = C_plus;
939 const RSHelicityAmpl & hampl_plus = hamplmod->Compute(res, fFKR);
940 fp0m = hampl_plus.Amp0Minus();
941 fp0p = hampl_plus.Amp0Plus();
942 }
943
944 double JRtSqrt2 = kSqrt2*JR;
945
946 // The signs of Hres are not correct. Now we consider these signs are exactly equal to ones in the latest version
947 // of code provided by Minoo Kabirnezhad, however pay attention to Cjsgn_plus
948 // which differ from what stated in refs. 1 and 2.
949 // More details about other problems can be found in Ref. 12, see section "Ambiguity in calculation of signs of the amplitudes"
950 // (most important conclusion in subsection "Paradox and probable explanation").
951
952 // Helicity amplitudes V-A, eq. 23 - 25 and Table 3 of ref. 1
953 // Cjsgn_plus are C^j_{lS2z} for S2z=1/2 and given in Table 9 of ref. 4
954 // Cjsgn_minus are C^j_{lS2z} for S2z=-1/2 and all equal to 1 (see Table 7 of ref. 4)
955
956 // The sign of the following amplitude is opposite to one from original code, because of it will be change
957 // when it will be multiplied by Wigner functions d^j_{\lambda\mu}
959 JRtSqrt2*Dsgn*Cjsgn_plus*kappa_f_BW*fp3*std::complex<double>
962
964 -JRtSqrt2*Dsgn*Cjsgn_plus*kappa_f_BW*fp3*std::complex<double>
967
969 JRtSqrt2*Dsgn*Cjsgn_plus*kappa_f_BW*fp1*std::complex<double>
972
973
974 // The sign of the following amplitude is opposite to one from original code, because of it will be change
975 // when it will be multiplied by Wigner functions d^j_{\lambda\mu}
977 -JRtSqrt2*Dsgn*Cjsgn_plus*kappa_f_BW*fp1*std::complex<double>
980
981
982
984 -JRtSqrt2*Dsgn*Cjsgn_plus*kappa_f_BW*fm1*std::complex<double>
987
989 JRtSqrt2*Dsgn*Cjsgn_plus*kappa_f_BW*fm1*std::complex<double>
992
994 -JRtSqrt2*Dsgn*Cjsgn_plus*kappa_f_BW*fm3*std::complex<double>
997
999 JRtSqrt2*Dsgn*Cjsgn_plus*kappa_f_BW*fm3*std::complex<double>
1002
1003
1004
1006 -k_sqrtQ2*JRtSqrt2*Dsgn*kappa_f_BW*fm0m*std::complex<double>
1009
1010 // The sign of the following amplitude is opposite to one from original code, because of it will be change
1011 // when it will be multiplied by Wigner functions d^j_{\lambda\mu}
1013 k_sqrtQ2*JRtSqrt2*Dsgn*kappa_f_BW*fm0m*std::complex<double>
1016
1018 k_sqrtQ2*JRtSqrt2*Dsgn*kappa_f_BW*fm0p*std::complex<double>
1021
1023 -k_sqrtQ2*JRtSqrt2*Dsgn*kappa_f_BW*fm0p*std::complex<double>
1026
1027
1028
1030 -k_sqrtQ2*JRtSqrt2*Dsgn*kappa_f_BW*fp0m*std::complex<double>
1033
1034 // The sign of the following amplitude is opposite to one from original code, because of it will be change
1035 // when it will be multiplied by Wigner functions d^j_{\lambda\mu}
1037 k_sqrtQ2*JRtSqrt2*Dsgn*kappa_f_BW*fp0m*std::complex<double>
1040
1042 k_sqrtQ2*JRtSqrt2*Dsgn*kappa_f_BW*fp0p*std::complex<double>
1045
1047 -k_sqrtQ2*JRtSqrt2*Dsgn*kappa_f_BW*fp0p*std::complex<double>
1050 } //end resonances loop
1051
1052 // Sum of all helicities
1055
1056 double pt_1 = 1;
1057 double pt_2 = 3*CosTheta;
1058 double pt_3 = 15./2*CosTheta*CosTheta - 3./2;
1059 double pt_4 = 35./2*CosTheta*CosTheta*CosTheta - 15./2*CosTheta;
1060
1061 // Wigner functions in terms of pion polar angle (eq. B4 of ref. 1)
1062 double d[4][2][2];
1063 d[0][0][1] = CosHalfTheta;
1064 d[0][0][0] =-SinHalfTheta;
1065 d[0][1][1] = 0;
1066 d[0][1][0] = 0;
1067
1068 d[1][0][1] = CosHalfTheta*(pt_2 - pt_1)/2;
1069 d[1][0][0] =-SinHalfTheta*(pt_2 + pt_1)/2;
1070 d[1][1][1] =-SinHalfTheta*(k1_Sqrt3*pt_2 + kSqrt3*pt_1)/2;
1071 d[1][1][0] = CosHalfTheta*(-k1_Sqrt3*pt_2 + kSqrt3*pt_1)/2;
1072
1073 d[2][0][1] = CosHalfTheta*(pt_3 - pt_2)/3;
1074 d[2][0][0] =-SinHalfTheta*(pt_3 + pt_2)/3;
1075 d[2][1][1] =-SinHalfTheta*(k1_Sqrt2*pt_3 + kSqrt2*pt_2)/3;
1076 d[2][1][0] = CosHalfTheta*(-k1_Sqrt2*pt_3 + kSqrt2*pt_2)/3;
1077
1078 d[3][0][1] = CosHalfTheta*(pt_4 - pt_3)/4;
1079 d[3][0][0] =-SinHalfTheta*(pt_4 + pt_3)/4;
1080 d[3][1][1] =-SinHalfTheta*(kSqrt3_5*pt_4 + kSqrt5_3*pt_3)/4;
1081 d[3][1][0] = CosHalfTheta*(-kSqrt3_5*pt_4 + kSqrt5_3*pt_3)/4;
1082
1083
1085 {
1086 int lambda_k = Lambda(lk);
1088 {
1089 int lambda_2 = Lambda(l2);
1091 {
1092 int lambda_1 = Lambda(l1);
1093 int lambda = lambda_k - lambda_1;
1094 int mu = -lambda_2;
1095 // Symmetry properties of Wigner d-functions, see eq. A1 from ref. 7
1096 int symmetry_factor = 1;
1097 int itemp;
1098 if (lambda < 0)
1099 {
1100 if (TMath::Abs(lambda)>TMath::Abs(mu)) //swap 1 + swap 2
1101 {
1102 symmetry_factor = TMath::Power(-1, (lambda - mu)/2);
1103 lambda*=-1;
1104 mu*=-1;
1105 }
1106 else
1107 {
1108 if (mu<0) //swap 1
1109 {
1110 itemp = lambda;
1111 lambda = -mu;
1112 mu = -itemp;
1113 }
1114 else //swap 2
1115 {
1116 symmetry_factor = TMath::Power(-1, (lambda - mu)/2);
1117 itemp = lambda;
1118 lambda = mu;
1119 mu = itemp;
1120 }
1121 }
1122 }
1123
1124 for (auto res : fResList)
1125 {
1126 // Get baryon resonance parameters
1127 int IR = utils::res::Isospin (res);
1128 int JR = (utils::res::AngularMom (res) + 1)/2;
1129 if (SppChannel::FinStateIsospin(spp_channel) == 3 && IR == 1) // skip resonances with I=1/2 if isospin of final state is 3/2
1130 continue;
1131 // Eq. 24 of ref. 1
1132 if (IR == 3)
1133 sum3(lk, l2, l1) += symmetry_factor*d[JR - 1][(lambda - 1)/2][(mu + 1)/2]*Hres(res, lk, l2, l1);
1134
1135
1136 if (IR == 1)
1137 sum1(lk, l2, l1) += symmetry_factor*d[JR - 1][(lambda - 1)/2][(mu + 1)/2]*Hres(res, lk, l2, l1);
1138 }
1139 }
1140 }
1141 }
1142
1143 // Isospin Clebsch–Gordan coefficients to sum amplitudes for I=1/2 and I=3/2, see eq.25 and Table 2 of ref. 1
1144 double C1 = SppChannel::Isospin1Coefficients(spp_channel);
1145 double C3 = SppChannel::Isospin3Coefficients(spp_channel);
1146
1147
1148 double g = fFermiConstant ;
1149 if(is_CC) g *= fVud;
1150
1151 double Lcoeff= abs_mom_k_L/2/kSqrt2/E;
1152 // Eq. 3.59 of ref. 2, which is multiplied by Q to avoid singularity at Q2=0
1153 double xsec0 = TMath::Power(g/8/kPi/kPi, 2)*abs_mom_q*Lcoeff*Lcoeff/abs_mom_k_L2/2;
1154 // We divide xsec0 by Q2 due to redifinition of Lcoeff
1155
1156
1157 if (kpsdim == 4)
1158 {
1160 {
1162 {
1163 // Eq. 18 ref. 1
1164 xsec +=
1165 TMath::Power(
1166 std::abs(
1167 C_L_minus*(
1169 C1*sum1( BosonPolarization::LEFT, l2, l1) + C3*sum3( BosonPolarization::LEFT, l2, l1)
1170 ) +
1171 C_R_minus*(
1173 C1*sum1( BosonPolarization::RIGHT, l2, l1) + C3*sum3( BosonPolarization::RIGHT, l2, l1)
1174 ) +
1175 C_S_minus*(
1177 C1*sum1( BosonPolarization::MINUS0, l2, l1) + C3*sum3( BosonPolarization::MINUS0, l2, l1)
1178 )
1179 )
1180 , 2)
1181 +
1182 TMath::Power(
1183 std::abs(
1184 C_L_plus*(
1186 C1*sum1( BosonPolarization::LEFT, l2, l1) + C3*sum3( BosonPolarization::LEFT, l2, l1)
1187 ) +
1188 C_R_plus*(
1190 C1*sum1( BosonPolarization::RIGHT, l2, l1) + C3*sum3( BosonPolarization::RIGHT, l2, l1)
1191 ) +
1192 C_S_plus*(
1194 C1*sum1( BosonPolarization::PLUS0, l2, l1) + C3*sum3( BosonPolarization::PLUS0, l2, l1)
1195 )
1196 )
1197 , 2);
1198 }
1199 }
1200 }
1201 else if (kpsdim == 3)
1202 {
1204 {
1206 {
1207 xsec +=
1208 (C_L_minus*C_L_minus + C_L_plus*C_L_plus)*(
1209 TMath::Power((Hbkg(Current::VECTOR, BosonPolarization::LEFT, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::LEFT, l2, l1)).real(), 2) +
1210 TMath::Power(std::abs(C1*sum1( BosonPolarization::LEFT, l2, l1) + C3*sum3( BosonPolarization::LEFT, l2, l1)), 2) +
1211 2*( Hbkg(Current::VECTOR, BosonPolarization::LEFT, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::LEFT, l2, l1)).real()*
1212 (C1*sum1( BosonPolarization::LEFT, l2, l1) + C3*sum3( BosonPolarization::LEFT, l2, l1)).real()
1213 )+
1214 (C_R_minus*C_R_minus + C_R_plus*C_R_plus)*(
1215 TMath::Power((Hbkg(Current::VECTOR, BosonPolarization::RIGHT, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::RIGHT, l2, l1)).real(), 2) +
1216 TMath::Power(std::abs(C1*sum1( BosonPolarization::RIGHT, l2, l1) + C3*sum3( BosonPolarization::RIGHT, l2, l1)), 2) +
1217 2*( Hbkg(Current::VECTOR, BosonPolarization::RIGHT, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::RIGHT, l2, l1)).real()*
1218 (C1*sum1( BosonPolarization::RIGHT, l2, l1) + C3*sum3( BosonPolarization::RIGHT, l2, l1)).real()
1219 )+
1220 (C_S_minus*C_S_minus)*(
1221 TMath::Power((Hbkg(Current::VECTOR, BosonPolarization::MINUS0, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::MINUS0, l2, l1)).real(), 2) +
1222 TMath::Power(std::abs(C1*sum1( BosonPolarization::MINUS0, l2, l1) + C3*sum3( BosonPolarization::MINUS0, l2, l1)), 2) +
1223 2*( Hbkg(Current::VECTOR, BosonPolarization::MINUS0, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::MINUS0, l2, l1)).real()*
1224 (C1*sum1( BosonPolarization::MINUS0, l2, l1) + C3*sum3( BosonPolarization::MINUS0, l2, l1)).real()
1225 )+
1226 (C_S_plus*C_S_plus)*(
1227 TMath::Power((Hbkg(Current::VECTOR, BosonPolarization::PLUS0, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::PLUS0, l2, l1)).real(), 2) +
1228 TMath::Power(std::abs(C1*sum1( BosonPolarization::PLUS0, l2, l1) + C3*sum3( BosonPolarization::PLUS0, l2, l1)), 2) +
1229 2*( Hbkg(Current::VECTOR, BosonPolarization::PLUS0, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::PLUS0, l2, l1)).real()*
1230 (C1*sum1( BosonPolarization::PLUS0, l2, l1) + C3*sum3( BosonPolarization::PLUS0, l2, l1)).real()
1231 );
1232 }
1233 }
1234 xsec*= 2*kPi;
1235 }
1236
1237 xsec*=xsec0;
1238
1239
1240 // The algorithm computes d^4xsec/dWdQ2dCosTheta_pidPhi_pi or d^3xsec/dWdQ2dCosTheta_pi
1241 // Check whether variable tranformation is needed
1242 if ( kps != kPSWQ2ctpphipfE && kps != kPSWQ2ctpfE )
1243 {
1244 double J = 1;
1245 if (kpsdim == 3)
1246 J = utils::kinematics::Jacobian(interaction, kPSWQ2ctpfE, kps);
1247 else if (kpsdim == 4)
1248 J = utils::kinematics::Jacobian(interaction, kPSWQ2ctpphipfE, kps);
1249 xsec *= J;
1250 }
1251
1252 // Apply given scaling factor
1253 if (is_CC) { xsec *= fXSecScaleCC; }
1254 else if (is_NC) { xsec *= fXSecScaleNC; }
1255
1256 // If requested return the free nucleon xsec even for input nuclear tgt
1257 if ( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;
1258
1259 int Z = target.Z();
1260 int A = target.A();
1261 int N = A-Z;
1262
1263 // Take into account the number of scattering centers in the target
1264 int NNucl = (SppChannel::InitStateNucleon(spp_channel) == kPdgProton) ? Z : N;
1265 xsec*=NNucl; // nuclear xsec (no nuclear suppression symmetry_factor)
1266
1267 if ( fUsePauliBlocking && A!=1 && kps == kPSWQ2ctpfE )
1268 {
1269 // Calculation of Pauli blocking according to refs. 9-11
1270 double P_Fermi = 0;
1271
1272 // Maximum value of Fermi momentum of target nucleon (GeV)
1273 if ( A<6 || ! fUseRFGParametrization )
1274 {
1275 // look up the Fermi momentum for this target
1276 FermiMomentumTablePool * kftp = FermiMomentumTablePool::Instance();
1277 const FermiMomentumTable * kft = kftp->GetTable(fKFTable);
1278 P_Fermi = kft->FindClosestKF(pdg::IonPdgCode(A, Z), nucpdgc);
1279 }
1280 else {
1281 // define the Fermi momentum for this target
1283 // correct the Fermi momentum for the struck nucleon
1284 if(is_p) { P_Fermi *= TMath::Power( 2.*Z/A, 1./3); }
1285 else { P_Fermi *= TMath::Power( 2.*N/A, 1./3); }
1286 }
1287
1288 double FactorPauli_RES = 1;
1289
1290 if (P_Fermi > 0)
1291 {
1292 if ( 2*P_Fermi < abs_mom_k-abs_mom_q )
1293 FactorPauli_RES = 1;
1294 if ( 2*P_Fermi >= abs_mom_k+abs_mom_q )
1295 FactorPauli_RES = ((3*abs_mom_k*abs_mom_k + abs_mom_q*abs_mom_q)/(2*P_Fermi) - (5*TMath::Power(abs_mom_k,4) + TMath::Power(abs_mom_q,4) + 10*abs_mom_k*abs_mom_k*abs_mom_q*abs_mom_q)/(40*TMath::Power(P_Fermi,3)))/(2*abs_mom_k);
1296 if ( 2*P_Fermi >= abs_mom_k-abs_mom_q && 2*P_Fermi <= abs_mom_k+abs_mom_q )
1297 FactorPauli_RES = ((abs_mom_q + abs_mom_k)*(abs_mom_q + abs_mom_k) - 4*P_Fermi*P_Fermi/5 - TMath::Power(abs_mom_k - abs_mom_q, 3)/(2*P_Fermi)+TMath::Power(abs_mom_k - abs_mom_q, 5)/(40*TMath::Power(P_Fermi, 3)))/(4*abs_mom_q*abs_mom_k);
1298 }
1299 xsec *= FactorPauli_RES;
1300 }
1301 return xsec;
1302
1303}
#define pFATAL
Definition Messenger.h:56
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE,...
Definition Messenger.h:96
const FermiMomentumTable * GetTable(string name)
static FermiMomentumTablePool * Instance(void)
double FindClosestKF(int target_pdgc, int nucleon_pdgc) const
const Target & Tgt(void) const
int ProbePdg(void) const
static string AsString(KinePhaseSpace_t kps)
double GetKV(KineVar_t kv) const
bool ValidKinematics(const Interaction *interaction) const
Is the input kinematical point a physically allowed one?
Iterator< NucleonPolarization, NucleonPolarization::MINUS, NucleonPolarization::PLUS > NucleonPolarizationIterator
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
int Lambda(NucleonPolarization l) const
Iterator< BosonPolarization, BosonPolarization::LEFT, BosonPolarization::PLUS0 > BosonPolarizationIterator
int PhaseFactor(BosonPolarization lk, NucleonPolarization l1, NucleonPolarization l2) const
static PDGLibrary * Instance(void)
TParticlePDG * Find(int pdgc, bool must_exist=true)
bool IsWeakNC(void) const
bool IsWeakCC(void) const
virtual const RSHelicityAmpl & Compute(Resonance_t res, const FKR &fkr) const =0
double Amp0Minus(void) const
double AmpPlus3(void) const
double Amp0Plus(void) const
double AmpPlus1(void) const
double AmpMinus3(void) const
double AmpMinus1(void) const
return helicity amplitude
static double Isospin1Coefficients(SppChannel_t channel)
Definition SppChannel.h:317
static double Isospin3Coefficients(SppChannel_t channel)
Definition SppChannel.h:280
static int FinStateIsospin(SppChannel_t channel)
Definition SppChannel.h:211
static int InitStateNucleon(SppChannel_t channel)
Definition SppChannel.h:103
static double BranchingRatio(Resonance_t res)
Definition SppChannel.h:357
int HitNucPdg(void) const
Definition Target.cxx:304
int Z(void) const
Definition Target.h:68
int A(void) const
Definition Target.h:70
bool IsNeutrino(int pdgc)
Definition PDGUtils.cxx:110
int IonPdgCode(int A, int Z)
Definition PDGUtils.cxx:71
bool IsProton(int pdgc)
Definition PDGUtils.cxx:336
bool IsAntiNeutrino(int pdgc)
Definition PDGUtils.cxx:118
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
double J(double q0, double q3, double Enu, double ml)
Definition MECUtils.cxx:147
double FermiMomentumForIsoscalarNucleonParametrization(const Target &target)
int AngularMom(Resonance_t res)
int Cjsgn_plus(Resonance_t res)
double Width(Resonance_t res)
resonance width (GeV)
int Dsgn(Resonance_t res)
int Isospin(Resonance_t res)
double Mass(Resonance_t res)
resonance mass (GeV)
int OrbitalAngularMom(Resonance_t res)
orbital angular momentum
int ResonanceIndex(Resonance_t res)
resonance idx, quark model / SU(6)
const int kPdgPiM
Definition PDGCodes.h:159
bool gAbortingInErr
Definition Messenger.cxx:34
const int kPdgProton
Definition PDGCodes.h:81
const int kPdgPi0
Definition PDGCodes.h:160
@ kPSWQ2ctpphipfE
const int kPdgNeutron
Definition PDGCodes.h:83
@ kKVphip
Definition KineVar.h:54
@ kKVctp
Definition KineVar.h:55
@ kSpp_vp_cc_10100
Definition SppChannel.h:50
@ kSpp_vbn_nc_10001
Definition SppChannel.h:66
@ kSpp_vbp_nc_10010
Definition SppChannel.h:63
@ kSpp_vp_nc_01100
Definition SppChannel.h:55
@ kSpp_vn_nc_01010
Definition SppChannel.h:56
@ kSpp_vbn_cc_01001
Definition SppChannel.h:59
@ kSpp_vn_cc_01100
Definition SppChannel.h:52
@ kSpp_vn_nc_10001
Definition SppChannel.h:57
@ kSpp_vbp_nc_01100
Definition SppChannel.h:64
@ kSpp_vn_cc_10010
Definition SppChannel.h:51
@ kSpp_vbp_cc_01010
Definition SppChannel.h:60
@ kSpp_vbp_cc_10001
Definition SppChannel.h:61
@ kSpp_vbn_nc_01010
Definition SppChannel.h:65
@ kSpp_vp_nc_10010
Definition SppChannel.h:54
const int kPdgPiP
Definition PDGCodes.h:158
const UInt_t kIAssumeFreeNucleon
Definition Interaction.h:49

References genie::Target::A(), genie::RSHelicityAmpl::Amp0Minus(), genie::RSHelicityAmpl::Amp0Plus(), genie::RSHelicityAmpl::AmpMinus1(), genie::RSHelicityAmpl::AmpMinus3(), genie::RSHelicityAmpl::AmpPlus1(), genie::RSHelicityAmpl::AmpPlus3(), genie::utils::res::AngularMom(), genie::KinePhaseSpace::AsString(), AXIAL, genie::SppChannel::BranchingRatio(), genie::utils::res::Cjsgn_plus(), genie::RSHelicityAmplModelI::Compute(), genie::utils::res::Dsgn(), f_pi, FA0, fBkgV0, fBkgV1, fBkgV2, fBkgV3, fBkgVWmax, fBkgVWmin, fCA50, fCv3, fCv4, fCv51, fCv52, fEMFormFactors, genie::utils::nuclear::FermiMomentumForIsoscalarNucleonParametrization(), fFermiConstant, fFKR, fFormFactors, fHAmplModelCC, fHAmplModelNCn, fHAmplModelNCp, genie::PDGLibrary::Find(), genie::FermiMomentumTable::FindClosestKF(), genie::SppChannel::FinStateIsospin(), fKFTable, fMa2, fMv2, fOmega, fResList, Frho0, fRho770Mass, genie::SppChannel::FromInteraction(), fSin2Wein, genie::Interaction::FSPrimLepton(), fUseAuthorCode, fUsePauliBlocking, fUseRFGParametrization, fVud, fXSecScaleCC, fXSecScaleNC, genie::gAbortingInErr, genie::Kinematics::GetKV(), genie::FermiMomentumTablePool::GetTable(), genie::Target::HitNucPdg(), genie::SppChannel::InitStateNucleon(), genie::FermiMomentumTablePool::Instance(), genie::PDGLibrary::Instance(), genie::pdg::IonPdgCode(), genie::pdg::IsAntiNeutrino(), genie::pdg::IsNeutrino(), genie::utils::res::Isospin(), genie::SppChannel::Isospin1Coefficients(), genie::SppChannel::Isospin3Coefficients(), genie::pdg::IsProton(), genie::ProcessInfo::IsWeakCC(), genie::ProcessInfo::IsWeakNC(), genie::utils::kinematics::Jacobian(), genie::constants::k1_Sqrt2, genie::constants::k1_Sqrt3, genie::kIAssumeFreeNucleon, genie::kKVctp, genie::kKVphip, genie::kPdgNeutron, genie::kPdgPi0, genie::kPdgPiM, genie::kPdgPiP, genie::kPdgProton, genie::constants::kPi, genie::kPSWQ2ctpfE, genie::kPSWQ2ctpphipfE, genie::kRfHitNucRest, genie::kSpp_vbn_cc_01001, genie::kSpp_vbn_nc_01010, genie::kSpp_vbn_nc_10001, genie::kSpp_vbp_cc_01010, genie::kSpp_vbp_cc_10001, genie::kSpp_vbp_nc_01100, genie::kSpp_vbp_nc_10010, genie::kSpp_vn_cc_01100, genie::kSpp_vn_cc_10010, genie::kSpp_vn_nc_01010, genie::kSpp_vn_nc_10001, genie::kSpp_vp_cc_10100, genie::kSpp_vp_nc_01100, genie::kSpp_vp_nc_10010, genie::constants::kSqrt2, genie::constants::kSqrt3, genie::constants::kSqrt3_5, genie::constants::kSqrt5_3, Lambda(), LEFT, LOG, genie::utils::res::Mass(), MINUS, MINUS0, genie::utils::res::OrbitalAngularMom(), pFATAL, PhaseFactor(), PLUS, PLUS0, genie::InitialState::ProbeE(), genie::InitialState::ProbePdg(), genie::Kinematics::Q2(), genie::utils::res::ResonanceIndex(), RIGHT, genie::InitialState::Tgt(), ValidKinematics(), ValidProcess(), VECTOR, genie::Kinematics::W(), genie::utils::res::Width(), and genie::Target::Z().

Member Data Documentation

◆ f_pi

double genie::MKSPPPXSec2020::f_pi
private

Constant for pion-nucleon interaction.

Definition at line 398 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ FA0

double genie::MKSPPPXSec2020::FA0
private

Axial coupling (value of axial form factor at Q2=0)

Definition at line 399 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fBkgV0

double genie::MKSPPPXSec2020::fBkgV0
private

Definition at line 411 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fBkgV1

double genie::MKSPPPXSec2020::fBkgV1
private

Definition at line 410 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fBkgV2

double genie::MKSPPPXSec2020::fBkgV2
private

Definition at line 409 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fBkgV3

double genie::MKSPPPXSec2020::fBkgV3
private

Definition at line 408 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fBkgVWmax

double genie::MKSPPPXSec2020::fBkgVWmax
private

Definition at line 407 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fBkgVWmin

double genie::MKSPPPXSec2020::fBkgVWmin
private

Parameters for vector virtual form factor for background contribution, which equal to:
1, W<VWmin
V3*W^3+V2*W^2+V1*W+V0 VWmin<W<VWmax 0 W>VWmax

Definition at line 406 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fCA50

double genie::MKSPPPXSec2020::fCA50
private

FKR parameter Zeta.

Definition at line 376 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fCv3

double genie::MKSPPPXSec2020::fCv3
private

GV calculation coeffs.

Definition at line 380 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fCv4

double genie::MKSPPPXSec2020::fCv4
private

Definition at line 381 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fCv51

double genie::MKSPPPXSec2020::fCv51
private

Definition at line 382 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fCv52

double genie::MKSPPPXSec2020::fCv52
private

Definition at line 383 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fElFFModel

const ELFormFactorsModelI* genie::MKSPPPXSec2020::fElFFModel
private

Elastic form facors model for background contribution.

Definition at line 388 of file MKSPPPXSec2020.h.

◆ fEMFormFactors

QELFormFactors genie::MKSPPPXSec2020::fEMFormFactors
mutableprivate

Electromagnetic form facors for background contribution.

Definition at line 397 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fEMFormFactorsModel

const QELFormFactorsModelI* genie::MKSPPPXSec2020::fEMFormFactorsModel
private

Electromagnetic form factors model for background contribution.

Definition at line 390 of file MKSPPPXSec2020.h.

Referenced by LoadConfig().

◆ fFermiConstant

double genie::MKSPPPXSec2020::fFermiConstant
private

Definition at line 375 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fFKR

FKR genie::MKSPPPXSec2020::fFKR
mutableprivate

Definition at line 368 of file MKSPPPXSec2020.h.

Referenced by XSec().

◆ fFormFactors

QELFormFactors genie::MKSPPPXSec2020::fFormFactors
mutableprivate

Quasielastic form facors for background contribution.

Definition at line 396 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fFormFactorsModel

const QELFormFactorsModelI* genie::MKSPPPXSec2020::fFormFactorsModel
private

Quasielastic form facors model for background contribution.

Definition at line 389 of file MKSPPPXSec2020.h.

Referenced by LoadConfig().

◆ fHAmplModelCC

const RSHelicityAmplModelI* genie::MKSPPPXSec2020::fHAmplModelCC
private

Definition at line 369 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fHAmplModelNCn

const RSHelicityAmplModelI* genie::MKSPPPXSec2020::fHAmplModelNCn
private

Definition at line 371 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fHAmplModelNCp

const RSHelicityAmplModelI* genie::MKSPPPXSec2020::fHAmplModelNCp
private

Definition at line 370 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fKFTable

string genie::MKSPPPXSec2020::fKFTable
private

Table of Fermi momentum (kF) constants for various nuclei.

Definition at line 392 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fMa2

double genie::MKSPPPXSec2020::fMa2
private

(axial mass)^2

Definition at line 378 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fMv2

double genie::MKSPPPXSec2020::fMv2
private

(vector mass)^2

Definition at line 379 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fOmega

double genie::MKSPPPXSec2020::fOmega
private

FKR parameter Omega.

Definition at line 377 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fResList

BaryonResList genie::MKSPPPXSec2020::fResList
private

Definition at line 419 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ Frho0

double genie::MKSPPPXSec2020::Frho0
private

Value of form factor F_rho at t=0

Definition at line 400 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fRho770Mass

double genie::MKSPPPXSec2020::fRho770Mass
private

Mass of rho(770) meson.

Definition at line 412 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fSin2Wein

double genie::MKSPPPXSec2020::fSin2Wein
private

sin^2(Weingberg angle)

Definition at line 384 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fUseAuthorCode

bool genie::MKSPPPXSec2020::fUseAuthorCode
private

Use author code?

Definition at line 415 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fUsePauliBlocking

bool genie::MKSPPPXSec2020::fUsePauliBlocking
private

Account for Pauli blocking?

Definition at line 394 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fUseRFGParametrization

bool genie::MKSPPPXSec2020::fUseRFGParametrization
private

Use parametrization for fermi momentum insted of table?

Definition at line 393 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fVud

double genie::MKSPPPXSec2020::fVud
private

|Vud| (magnitude ud-element of CKM-matrix)

Definition at line 385 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fWmax

double genie::MKSPPPXSec2020::fWmax
private

The value above which the partial cross section is set to zero (if negative then not appliciable)

Definition at line 413 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and ValidKinematics().

◆ fXSecIntegrator

const XSecIntegratorI* genie::MKSPPPXSec2020::fXSecIntegrator
private

Definition at line 417 of file MKSPPPXSec2020.h.

Referenced by Integral(), and LoadConfig().

◆ fXSecScaleCC

double genie::MKSPPPXSec2020::fXSecScaleCC
private

External CC xsec scaling factor.

Definition at line 386 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

◆ fXSecScaleNC

double genie::MKSPPPXSec2020::fXSecScaleNC
private

External NC xsec scaling factor.

Definition at line 387 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().


The documentation for this class was generated from the following files: