GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Classes | Public Member Functions | Private Types | Private Member Functions | Private Attributes | List of all members
genie::MKSPPPXSec2020 Class Reference

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: More...

#include <MKSPPPXSec2020.h>

Inheritance diagram for genie::MKSPPPXSec2020:
Inheritance graph
[legend]
Collaboration diagram for genie::MKSPPPXSec2020:
Collaboration graph
[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. More...
 
double Integral (const Interaction *i) const
 
bool ValidProcess (const Interaction *i) const
 Can this cross section algorithm handle the input process? More...
 
bool ValidKinematics (const Interaction *interaction) const
 Is the input kinematical point a physically allowed one? More...
 
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 Registry & GetConfig (void) const
 
Registry * GetOwnedConfig (void)
 
virtual const AlgId & Id (void) const
 Get algorithm ID. More...
 
virtual AlgStatus_t GetStatus (void) const
 Get algorithm status. More...
 
virtual bool AllowReconfig (void) const
 
virtual AlgCmp_t Compare (const Algorithm *alg) const
 Compare with input algorithm. More...
 
virtual void SetId (const AlgId &id)
 Set algorithm ID. More...
 
virtual void SetId (string name, string config)
 
const Algorithm * SubAlg (const RgKey &registry_key) const
 
void AdoptConfig (void)
 
void AdoptSubstructure (void)
 
virtual void Print (ostream &stream) const
 Print algorithm info. More...
 

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 RSHelicityAmplModelI * fHAmplModelCC
 
const RSHelicityAmplModelI * fHAmplModelNCp
 
const RSHelicityAmplModelI * fHAmplModelNCn
 
double fFermiConstant
 
double fCA50
 FKR parameter Zeta. More...
 
double fOmega
 FKR parameter Omega. More...
 
double fMa2
 (axial mass)^2 More...
 
double fMv2
 (vector mass)^2 More...
 
double fCv3
 GV calculation coeffs. More...
 
double fCv4
 
double fCv51
 
double fCv52
 
double fSin2Wein
 sin^2(Weingberg angle) More...
 
double fVud
 |Vud| (magnitude ud-element of CKM-matrix) More...
 
double fXSecScaleCC
 External CC xsec scaling factor. More...
 
double fXSecScaleNC
 External NC xsec scaling factor. More...
 
const ELFormFactorsModelI * fElFFModel
 Elastic form facors model for background contribution. More...
 
const QELFormFactorsModelI * fFormFactorsModel
 Quasielastic form facors model for background contribution. More...
 
const QELFormFactorsModelI * fEMFormFactorsModel
 Electromagnetic form factors model for background contribution. More...
 
string fKFTable
 Table of Fermi momentum (kF) constants for various nuclei. More...
 
bool fUseRFGParametrization
 Use parametrization for fermi momentum insted of table? More...
 
bool fUsePauliBlocking
 Account for Pauli blocking? More...
 
QELFormFactors fFormFactors
 Quasielastic form facors for background contribution. More...
 
QELFormFactors fEMFormFactors
 Electromagnetic form facors for background contribution. More...
 
double f_pi
 Constant for pion-nucleon interaction. More...
 
double FA0
 Axial coupling (value of axial form factor at Q2=0) More...
 
double Frho0
 
double fBkgVWmin
 
double fBkgVWmax
 
double fBkgV3
 
double fBkgV2
 
double fBkgV1
 
double fBkgV0
 
double fRho770Mass
 Mass of rho(770) meson. More...
 
double fWmax
 The value above which the partial cross section is set to zero (if negative then not appliciable) More...
 
bool fUseAuthorCode
 Use author code? More...
 
const XSecIntegratorI * fXSecIntegrator
 
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)
 
Registry * ExtractLocalConfig (const Registry &in) const
 
Registry * ExtractLowerConfig (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. More...
 
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. More...
 
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. More...
 
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 More...
 
int AddLowRegistry (Registry *rp, bool owns=true)
 add registry with lowest priority, also update ownership More...
 
int MergeTopRegistry (const Registry &r)
 
int AddTopRegisties (const vector< Registry * > &rs, bool owns=false)
 Add registries with top priority, also udated Ownerships. More...
 
- Protected Attributes inherited from genie::Algorithm
bool fAllowReconfig
 
bool fOwnsSubstruc
 true if it owns its substructure (sub-algs,...) More...
 
AlgId fID
 algorithm name and configuration set More...
 
vector< Registry * > fConfVect
 
vector< bool > fOwnerships
 ownership for every registry in fConfVect More...
 
AlgStatus_t fStatus
 algorithm execution status More...
 
AlgMap * fOwnedSubAlgMp
 local pool for owned sub-algs (taken out of the factory pool) More...
 

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:
  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:
Nov 12, 2019
License:
Copyright (c) 2003-2024, 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

using genie::MKSPPPXSec2020::BosonPolarizationIterator = Iterator<BosonPolarization, BosonPolarization::LEFT, BosonPolarization::PLUS0>
private

Definition at line 154 of file MKSPPPXSec2020.h.

using genie::MKSPPPXSec2020::CurrentIterator = Iterator<Current, Current::VECTOR, Current::AXIAL>
private

Definition at line 153 of file MKSPPPXSec2020.h.

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.

using genie::MKSPPPXSec2020::NucleonPolarizationIterator = Iterator<NucleonPolarization, NucleonPolarization::MINUS, NucleonPolarization::PLUS>
private

Definition at line 155 of file MKSPPPXSec2020.h.

Member Enumeration Documentation

enum genie::MKSPPPXSec2020::BosonPolarization
private
Enumerator
LEFT 
RIGHT 
MINUS0 
PLUS0 

Definition at line 129 of file MKSPPPXSec2020.h.

129 { LEFT, RIGHT, MINUS0, PLUS0 };
enum genie::MKSPPPXSec2020::Current
private
Enumerator
VECTOR 
AXIAL 

Definition at line 128 of file MKSPPPXSec2020.h.

128 { VECTOR, AXIAL };
enum genie::MKSPPPXSec2020::NucleonPolarization
private
Enumerator
MINUS 
PLUS 

Definition at line 130 of file MKSPPPXSec2020.h.

130 { MINUS, PLUS };

Constructor & Destructor Documentation

MKSPPPXSec2020::MKSPPPXSec2020 ( )

Definition at line 52 of file MKSPPPXSec2020.cxx.

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

Definition at line 58 of file MKSPPPXSec2020.cxx.

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

Definition at line 64 of file MKSPPPXSec2020.cxx.

65 {
66 
67 }

Member Function Documentation

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.

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

1375 {
1376  Algorithm::Configure(config);
1377  this->LoadConfig();
1378 }
genie::Algorithm::Configure
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
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.

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

1381 {
1382  Algorithm::Configure(config);
1383  this->LoadConfig();
1384 }
genie::Algorithm::Configure
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
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.

References fXSecIntegrator, and genie::XSecIntegratorI::Integrate().

1307 {
1308  double xsec = fXSecIntegrator->Integrate(this,interaction);
1309  return xsec;
1310 }
genie::MKSPPPXSec2020::fXSecIntegrator
const XSecIntegratorI * fXSecIntegrator
Definition: MKSPPPXSec2020.h:417
genie::XSecIntegratorI::Integrate
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
int MKSPPPXSec2020::Lambda ( NucleonPolarization  l) const
inlineprivate

Definition at line 1327 of file MKSPPPXSec2020.cxx.

Referenced by XSec().

1328 {
1329  return 2*l - 1;
1330 }
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 }
void MKSPPPXSec2020::LoadConfig ( void  )
private

Definition at line 1386 of file MKSPPPXSec2020.cxx.

References genie::BaryonResList::Clear(), genie::BaryonResList::DecodeFromNameList(), 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(), genie::QELFormFactors::SetModel(), and genie::Algorithm::SubAlg().

Referenced by Configure().

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 }
genie::MKSPPPXSec2020::fXSecIntegrator
const XSecIntegratorI * fXSecIntegrator
Definition: MKSPPPXSec2020.h:417
genie::XSecIntegratorI
Cross Section Integrator Interface.
Definition: XSecIntegratorI.h:29
genie::QELFormFactors::SetModel
void SetModel(const QELFormFactorsModelI *model)
Attach an algorithm.
Definition: QELFormFactors.cxx:43
genie::BaryonResList::DecodeFromNameList
void DecodeFromNameList(string list, string delimiter=",")
Definition: BaryonResList.cxx:101
genie::MKSPPPXSec2020::fSin2Wein
double fSin2Wein
sin^2(Weingberg angle)
Definition: MKSPPPXSec2020.h:384
genie::MKSPPPXSec2020::fHAmplModelCC
const RSHelicityAmplModelI * fHAmplModelCC
Definition: MKSPPPXSec2020.h:369
genie::MKSPPPXSec2020::fUseRFGParametrization
bool fUseRFGParametrization
Use parametrization for fermi momentum insted of table?
Definition: MKSPPPXSec2020.h:393
genie::MKSPPPXSec2020::fWmax
double fWmax
The value above which the partial cross section is set to zero (if negative then not appliciable) ...
Definition: MKSPPPXSec2020.h:413
genie::MKSPPPXSec2020::fMa2
double fMa2
(axial mass)^2
Definition: MKSPPPXSec2020.h:378
genie::MKSPPPXSec2020::fXSecScaleNC
double fXSecScaleNC
External NC xsec scaling factor.
Definition: MKSPPPXSec2020.h:387
genie::MKSPPPXSec2020::fCA50
double fCA50
FKR parameter Zeta.
Definition: MKSPPPXSec2020.h:376
genie::MKSPPPXSec2020::fResList
BaryonResList fResList
Definition: MKSPPPXSec2020.h:419
genie::QELFormFactorsModelI
Pure abstract base class. Defines the QELFormFactorsModelI interface to be implemented by any algorit...
Definition: QELFormFactorsModelI.h:28
genie::MKSPPPXSec2020::fHAmplModelNCp
const RSHelicityAmplModelI * fHAmplModelNCp
Definition: MKSPPPXSec2020.h:370
genie::MKSPPPXSec2020::fRho770Mass
double fRho770Mass
Mass of rho(770) meson.
Definition: MKSPPPXSec2020.h:412
genie::MKSPPPXSec2020::fMv2
double fMv2
(vector mass)^2
Definition: MKSPPPXSec2020.h:379
genie::MKSPPPXSec2020::fHAmplModelNCn
const RSHelicityAmplModelI * fHAmplModelNCn
Definition: MKSPPPXSec2020.h:371
genie::MKSPPPXSec2020::FA0
double FA0
Axial coupling (value of axial form factor at Q2=0)
Definition: MKSPPPXSec2020.h:399
genie::MKSPPPXSec2020::fEMFormFactorsModel
const QELFormFactorsModelI * fEMFormFactorsModel
Electromagnetic form factors model for background contribution.
Definition: MKSPPPXSec2020.h:390
genie::MKSPPPXSec2020::f_pi
double f_pi
Constant for pion-nucleon interaction.
Definition: MKSPPPXSec2020.h:398
genie::RSHelicityAmplModelI
Pure abstract base class. Defines the RSHelicityAmplModelI interface.
Definition: RSHelicityAmplModelI.h:27
genie::MKSPPPXSec2020::fKFTable
string fKFTable
Table of Fermi momentum (kF) constants for various nuclei.
Definition: MKSPPPXSec2020.h:392
genie::MKSPPPXSec2020::fCv3
double fCv3
GV calculation coeffs.
Definition: MKSPPPXSec2020.h:380
genie::MKSPPPXSec2020::fUsePauliBlocking
bool fUsePauliBlocking
Account for Pauli blocking?
Definition: MKSPPPXSec2020.h:394
genie::MKSPPPXSec2020::fOmega
double fOmega
FKR parameter Omega.
Definition: MKSPPPXSec2020.h:377
genie::MKSPPPXSec2020::fVud
double fVud
|Vud| (magnitude ud-element of CKM-matrix)
Definition: MKSPPPXSec2020.h:385
genie::MKSPPPXSec2020::fUseAuthorCode
bool fUseAuthorCode
Use author code?
Definition: MKSPPPXSec2020.h:415
genie::MKSPPPXSec2020::fFormFactors
QELFormFactors fFormFactors
Quasielastic form facors for background contribution.
Definition: MKSPPPXSec2020.h:396
genie::Algorithm::GetParamDef
bool GetParamDef(const RgKey &name, T &p, const T &def) const
genie::Algorithm::GetParam
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
genie::MKSPPPXSec2020::fXSecScaleCC
double fXSecScaleCC
External CC xsec scaling factor.
Definition: MKSPPPXSec2020.h:386
genie::MKSPPPXSec2020::fEMFormFactors
QELFormFactors fEMFormFactors
Electromagnetic form facors for background contribution.
Definition: MKSPPPXSec2020.h:397
genie::MKSPPPXSec2020::fFormFactorsModel
const QELFormFactorsModelI * fFormFactorsModel
Quasielastic form facors model for background contribution.
Definition: MKSPPPXSec2020.h:389
genie::Algorithm::SubAlg
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345
int MKSPPPXSec2020::PhaseFactor ( BosonPolarization  lk,
NucleonPolarization  l1,
NucleonPolarization  l2 
) const
inlineprivate

Definition at line 1337 of file MKSPPPXSec2020.cxx.

Referenced by XSec().

1338 {
1339  return lk*(lk*(4*lk - 21) + 29)/6 - l1 - l2;
1340 }
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.

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

Referenced by XSec().

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 }
genie::Kinematics::W
double W(bool selected=false) const
Definition: Kinematics.cxx:157
genie::utils::kinematics::Q2
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
genie::Range1D_t
A simple [min,max] interval for doubles.
Definition: Range1.h:42
genie::Kinematics
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
genie::KPhaseSpace::Threshold_SPP_iso
double Threshold_SPP_iso(void) const
Energy limit for resonance single pion production on isoscalar nucleon.
Definition: KPhaseSpace.cxx:999
genie::MKSPPPXSec2020::fWmax
double fWmax
The value above which the partial cross section is set to zero (if negative then not appliciable) ...
Definition: MKSPPPXSec2020.h:413
genie::utils::kinematics::W
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
genie::KPhaseSpace::Q2Lim_W_SPP_iso
Range1D_t Q2Lim_W_SPP_iso(void) const
Q2 limits @ fixed W for resonance single pion production on isoscalar nucleon.
Definition: KPhaseSpace.cxx:1110
genie::KPhaseSpace
Kinematical phase space.
Definition: KPhaseSpace.h:33
genie::kISkipKinematicChk
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
genie::Range1D_t::max
double max
Definition: Range1.h:53
genie::KPhaseSpace::WLim_SPP_iso
Range1D_t WLim_SPP_iso(void) const
W limits for resonance single pion production on isoscalar nucleon.
Definition: KPhaseSpace.cxx:1042
genie::Range1D_t::min
double min
Definition: Range1.h:52
genie::Kinematics::Q2
double Q2(bool selected=false) const
Definition: Kinematics.cxx:125
genie::InitialState::ProbeE
double ProbeE(RefFrame_t rf) const
Definition: InitialState.cxx:384
genie::InitialState
Initial State information.
Definition: InitialState.h:48
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.

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

Referenced by XSec().

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 }
genie::SppChannel::FromInteraction
static SppChannel_t FromInteraction(const Interaction *interaction)
Definition: SppChannel.h:402
genie::SppChannel_t
enum genie::ESppChannel SppChannel_t
genie::kISkipProcessChk
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
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.

References genie::Target::A(), genie::units::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(), genie::FKR::B, genie::SppChannel::BranchingRatio(), genie::FKR::C, genie::QELFormFactors::Calculate(), genie::utils::res::Cjsgn_plus(), genie::RSHelicityAmplModelI::Compute(), genie::utils::res::Dsgn(), genie::QELFormFactors::F1V(), f_pi, genie::QELFormFactors::FA(), 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, genie::units::fm3, fMa2, fMv2, fOmega, fResList, Frho0, fRho770Mass, genie::SppChannel::FromInteraction(), fSin2Wein, genie::Interaction::FSPrimLepton(), fUseAuthorCode, fUsePauliBlocking, fUseRFGParametrization, fVud, fXSecScaleCC, fXSecScaleNC, genie::units::g, 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::mec::J(), 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(), genie::FKR::Lamda, LOG, genie::utils::res::Mass(), genie::utils::res::OrbitalAngularMom(), pFATAL, PhaseFactor(), genie::InitialState::ProbeE(), genie::InitialState::ProbePdg(), genie::Kinematics::Q2(), genie::utils::kinematics::Q2(), genie::FKR::R, genie::FKR::Ra, genie::utils::res::ResonanceIndex(), genie::FKR::Rminus, genie::FKR::Rplus, genie::FKR::Rv, genie::FKR::S, genie::units::s, genie::FKR::T, genie::FKR::Ta, genie::InitialState::Tgt(), genie::FKR::Tminus, genie::FKR::Tplus, genie::FKR::Tv, ValidKinematics(), ValidProcess(), genie::Kinematics::W(), genie::utils::kinematics::W(), genie::utils::res::Width(), genie::QELFormFactors::xiF2V(), and genie::Target::Z().

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
222  HelicityBkgAmp<std::complex<double>> Hbkg;
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 
346  Hbkg(Current::VECTOR, BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
347  k1_Sqrt2*SinTheta*CosHalfTheta*(sF_3 + sF_4)*std::complex<double>
348  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
349  TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
350 
351  Hbkg(Current::VECTOR, BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
352  -k1_Sqrt2*SinTheta*SinHalfTheta*(sF_3 - sF_4)*std::complex<double>
353  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
354  TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
355 
356  Hbkg(Current::VECTOR, BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
357  kSqrt2*(CosHalfTheta*(sF_1 - sF_2) - 0.5*SinTheta*SinHalfTheta*(sF_3 - sF_4))*std::complex<double>
358  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
359  TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
360 
361  Hbkg(Current::VECTOR, BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
362  -kSqrt2*(SinHalfTheta*(sF_1 + sF_2) + 0.5*SinTheta*CosHalfTheta*(sF_3 + sF_4))*std::complex<double>
363  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
364  TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
365 
366 
367  Hbkg(Current::VECTOR, BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
368  -kSqrt2*(SinHalfTheta*(sF_1 + sF_2) + 0.5*SinTheta*CosHalfTheta*(sF_3 + sF_4))*std::complex<double>
369  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
370  TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
371 
372  Hbkg(Current::VECTOR, BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
373  -kSqrt2*(CosHalfTheta*(sF_1 - sF_2) - 0.5*SinTheta*SinHalfTheta*(sF_3 - sF_4))*std::complex<double>
374  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
375  TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
376 
377  Hbkg(Current::VECTOR, BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
378  k1_Sqrt2*SinTheta*SinHalfTheta*(sF_3 - sF_4)*std::complex<double>
379  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
380  TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
381 
382  Hbkg(Current::VECTOR, BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
383  k1_Sqrt2*SinTheta*CosHalfTheta*(sF_3 + sF_4)*std::complex<double>
384  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
385  TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
386 
387 
388  if (is_CC || (is_NC && is_nu) )
389  {
390  Hbkg(Current::VECTOR, BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
391  CosHalfTheta*(k_0*eps_z_L - abs_mom_k*eps_zero_L)*(sF_5 + sF_6)/C_S_minus*std::complex<double>
392  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
393  TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
394 
395  Hbkg(Current::VECTOR, BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
396  -SinHalfTheta*(k_0*eps_z_L - abs_mom_k*eps_zero_L)*(sF_5 - sF_6)/C_S_minus*std::complex<double>
397  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
398  TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
399 
400  Hbkg(Current::VECTOR, BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
401  -SinHalfTheta*(k_0*eps_z_L - abs_mom_k*eps_zero_L)*(sF_5 - sF_6)/C_S_minus*std::complex<double>
402  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
403  TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
404 
405  Hbkg(Current::VECTOR, BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
406  -CosHalfTheta*(k_0*eps_z_L - abs_mom_k*eps_zero_L)*(sF_5 + sF_6)/C_S_minus*std::complex<double>
407  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
408  TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
409  }
410  if (is_CC || (is_NC && is_nubar) )
411  {
412  Hbkg(Current::VECTOR, BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
413  CosHalfTheta*(k_0*eps_z_R - abs_mom_k*eps_zero_R)*(sF_5 + sF_6)/C_S_plus*std::complex<double>
414  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
415  TMath::Sin(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
416 
417  Hbkg(Current::VECTOR, BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
418  -SinHalfTheta*(k_0*eps_z_R - abs_mom_k*eps_zero_R)*(sF_5 - sF_6)/C_S_plus*std::complex<double>
419  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
420  TMath::Sin(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
421 
422  Hbkg(Current::VECTOR, BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
423  -SinHalfTheta*(k_0*eps_z_R - abs_mom_k*eps_zero_R)*(sF_5 - sF_6)/C_S_plus*std::complex<double>
424  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
425  TMath::Sin(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
426 
427  Hbkg(Current::VECTOR, BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
428  -CosHalfTheta*(k_0*eps_z_R - abs_mom_k*eps_zero_R)*(sF_5 + sF_6)/C_S_plus*std::complex<double>
429  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
430  TMath::Sin(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
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 
493  HelicityBkgAmp<std::complex<double>> HbkgEM;
494 
495  HbkgEM(Current::VECTOR, BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
496  k1_Sqrt2*SinTheta*CosHalfTheta*(sFem_3 + sFem_4)*std::complex<double>
497  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
498  TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
499 
500  HbkgEM(Current::VECTOR, BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
501  -k1_Sqrt2*SinTheta*SinHalfTheta*(sFem_3 - sFem_4)*std::complex<double>
502  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
503  TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
504 
505  HbkgEM(Current::VECTOR, BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
506  kSqrt2*(CosHalfTheta*(sFem_1 - sFem_2) - 0.5*SinTheta*SinHalfTheta*(sFem_3 - sFem_4))*std::complex<double>
507  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
508  TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
509 
510  HbkgEM(Current::VECTOR, BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
511  -kSqrt2*(SinHalfTheta*(sFem_1 + sFem_2) + 0.5*SinTheta*CosHalfTheta*(sFem_3 + sFem_4))*std::complex<double>
512  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
513  TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
514 
515 
516  HbkgEM(Current::VECTOR, BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
517  -kSqrt2*(SinHalfTheta*(sFem_1 + sFem_2) + 0.5*SinTheta*CosHalfTheta*(sFem_3 + sFem_4))*std::complex<double>
518  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
519  TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
520 
521  HbkgEM(Current::VECTOR, BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
522  -kSqrt2*(CosHalfTheta*(sFem_1 - sFem_2) - 0.5*SinTheta*SinHalfTheta*(sFem_3 - sFem_4))*std::complex<double>
523  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
524  TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
525 
526  HbkgEM(Current::VECTOR, BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
527  k1_Sqrt2*SinTheta*SinHalfTheta*(sFem_3 - sFem_4)*std::complex<double>
528  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
529  TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
530 
531  HbkgEM(Current::VECTOR, BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
532  k1_Sqrt2*SinTheta*CosHalfTheta*(sFem_3 + sFem_4)*std::complex<double>
533  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
534  TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
535 
536 
537  if (is_nu)
538  {
539  HbkgEM(Current::VECTOR, BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
540  CosHalfTheta*(k_0*eps_z_L - abs_mom_k*eps_zero_L)*(sFem_5 + sFem_6)/C_S_minus*std::complex<double>
541  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
542  TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
543 
544  HbkgEM(Current::VECTOR, BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
545  -SinHalfTheta*(k_0*eps_z_L - abs_mom_k*eps_zero_L)*(sFem_5 - sFem_6)/C_S_minus*std::complex<double>
546  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
547  TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
548 
549  HbkgEM(Current::VECTOR, BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
550  -SinHalfTheta*(k_0*eps_z_L - abs_mom_k*eps_zero_L)*(sFem_5 - sFem_6)/C_S_minus*std::complex<double>
551  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
552  TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
553 
554  HbkgEM(Current::VECTOR, BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
555  -CosHalfTheta*(k_0*eps_z_L - abs_mom_k*eps_zero_L)*(sFem_5 + sFem_6)/C_S_minus*std::complex<double>
556  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
557  TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
558  }
559  else
560  {
561  HbkgEM(Current::VECTOR, BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
562  CosHalfTheta*(k_0*eps_z_R - abs_mom_k*eps_zero_R)*(sFem_5 + sFem_6)/C_S_plus*std::complex<double>
563  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
564  TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
565 
566  HbkgEM(Current::VECTOR, BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
567  -SinHalfTheta*(k_0*eps_z_R - abs_mom_k*eps_zero_R)*(sFem_5 - sFem_6)/C_S_plus*std::complex<double>
568  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
569  TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
570 
571  HbkgEM(Current::VECTOR, BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
572  -SinHalfTheta*(k_0*eps_z_R - abs_mom_k*eps_zero_R)*(sFem_5 - sFem_6)/C_S_plus*std::complex<double>
573  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
574  TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
575 
576  HbkgEM(Current::VECTOR, BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
577  -CosHalfTheta*(k_0*eps_z_R - abs_mom_k*eps_zero_R)*(sFem_5 + sFem_6)/C_S_plus*std::complex<double>
578  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
579  TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
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
726  Hbkg(Current::AXIAL, BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
727  k1_Sqrt2*SinTheta*CosHalfTheta*(sG_3 + sG_4)*std::complex<double>
728  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
729  TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
730 
731  Hbkg(Current::AXIAL, BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
732  k1_Sqrt2*SinTheta*SinHalfTheta*(sG_3 - sG_4)*std::complex<double>
733  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
734  TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
735 
736  Hbkg(Current::AXIAL, BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
737  kSqrt2*(CosHalfTheta*(sG_1 - sG_2) - 0.5*SinTheta*SinHalfTheta*(sG_3 - sG_4))*std::complex<double>
738  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
739  TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
740 
741  Hbkg(Current::AXIAL, BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
742  kSqrt2*(SinHalfTheta*(sG_1 + sG_2) + 0.5*SinTheta*CosHalfTheta*(sG_3 + sG_4))*std::complex<double>
743  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
744  TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
745 
746 
747  Hbkg(Current::AXIAL, BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
748  -kSqrt2*(SinHalfTheta*(sG_1 + sG_2) + 0.5*SinTheta*CosHalfTheta*(sG_3 + sG_4))*std::complex<double>
749  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
750  TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
751 
752  Hbkg(Current::AXIAL, BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
753  kSqrt2*(CosHalfTheta*(sG_1 - sG_2) - 0.5*SinTheta*SinHalfTheta*(sG_3 - sG_4))*std::complex<double>
754  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
755  TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
756 
757  Hbkg(Current::AXIAL, BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
758  k1_Sqrt2*SinTheta*SinHalfTheta*(sG_3 - sG_4)*std::complex<double>
759  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
760  TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
761 
762  Hbkg(Current::AXIAL, BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
763  -k1_Sqrt2*SinTheta*CosHalfTheta*(sG_3 + sG_4)*std::complex<double>
764  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
765  TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
766 
767  if (is_CC || (is_NC && is_nu) )
768  {
769  Hbkg(Current::AXIAL, BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
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>
771  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
772  TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
773 
774  Hbkg(Current::AXIAL, BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
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>
776  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
777  TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
778 
779  Hbkg(Current::AXIAL, BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
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>
781  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
782  TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
783 
784  Hbkg(Current::AXIAL, BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
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>
786  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
787  TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
788  }
789  if (is_CC || (is_NC && is_nubar) )
790  {
791  Hbkg(Current::AXIAL, BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
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>
793  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
794  TMath::Sin(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
795 
796  Hbkg(Current::AXIAL, BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
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>
798  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
799  TMath::Sin(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
800 
801  Hbkg(Current::AXIAL, BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
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>
803  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
804  TMath::Sin(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
805 
806  Hbkg(Current::AXIAL, BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
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>
808  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
809  TMath::Sin(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
810  }
811 
812  /* Resonace contribution */
813  HelicityAmpVminusARes<std::complex<double>> Hres;
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}
958  Hres(res, BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
959  JRtSqrt2*Dsgn*Cjsgn_plus*kappa_f_BW*fp3*std::complex<double>
960  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
961  TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
962 
963  Hres(res, BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
964  -JRtSqrt2*Dsgn*Cjsgn_plus*kappa_f_BW*fp3*std::complex<double>
965  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
966  TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
967 
968  Hres(res, BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
969  JRtSqrt2*Dsgn*Cjsgn_plus*kappa_f_BW*fp1*std::complex<double>
970  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
971  TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
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}
976  Hres(res, BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
977  -JRtSqrt2*Dsgn*Cjsgn_plus*kappa_f_BW*fp1*std::complex<double>
978  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
979  TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
980 
981 
982 
983  Hres(res, BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
984  -JRtSqrt2*Dsgn*Cjsgn_plus*kappa_f_BW*fm1*std::complex<double>
985  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
986  TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
987 
988  Hres(res, BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
989  JRtSqrt2*Dsgn*Cjsgn_plus*kappa_f_BW*fm1*std::complex<double>
990  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
991  TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
992 
993  Hres(res, BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
994  -JRtSqrt2*Dsgn*Cjsgn_plus*kappa_f_BW*fm3*std::complex<double>
995  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
996  TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
997 
998  Hres(res, BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
999  JRtSqrt2*Dsgn*Cjsgn_plus*kappa_f_BW*fm3*std::complex<double>
1000  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
1001  TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
1002 
1003 
1004 
1005  Hres(res, BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
1006  -k_sqrtQ2*JRtSqrt2*Dsgn*kappa_f_BW*fm0m*std::complex<double>
1007  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
1008  TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
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}
1012  Hres(res, BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
1013  k_sqrtQ2*JRtSqrt2*Dsgn*kappa_f_BW*fm0m*std::complex<double>
1014  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
1015  TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
1016 
1017  Hres(res, BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
1018  k_sqrtQ2*JRtSqrt2*Dsgn*kappa_f_BW*fm0p*std::complex<double>
1019  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
1020  TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
1021 
1022  Hres(res, BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
1023  -k_sqrtQ2*JRtSqrt2*Dsgn*kappa_f_BW*fm0p*std::complex<double>
1024  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
1025  TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
1026 
1027 
1028 
1029  Hres(res, BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
1030  -k_sqrtQ2*JRtSqrt2*Dsgn*kappa_f_BW*fp0m*std::complex<double>
1031  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
1032  TMath::Sin(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
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}
1036  Hres(res, BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
1037  k_sqrtQ2*JRtSqrt2*Dsgn*kappa_f_BW*fp0m*std::complex<double>
1038  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
1039  TMath::Sin(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
1040 
1041  Hres(res, BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
1042  k_sqrtQ2*JRtSqrt2*Dsgn*kappa_f_BW*fp0p*std::complex<double>
1043  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
1044  TMath::Sin(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
1045 
1046  Hres(res, BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
1047  -k_sqrtQ2*JRtSqrt2*Dsgn*kappa_f_BW*fp0p*std::complex<double>
1048  (TMath::Cos(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
1049  TMath::Sin(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
1050  } //end resonances loop
1051 
1052  // Sum of all helicities
1053  SumHelicityAmpVminusARes<std::complex<double>> sum3;
1054  SumHelicityAmpVminusARes<std::complex<double>> sum1;
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 
1084  for (BosonPolarization lk : BosonPolarizationIterator() )
1085  {
1086  int lambda_k = Lambda(lk);
1087  for (NucleonPolarization l2 : NucleonPolarizationIterator() )
1088  {
1089  int lambda_2 = Lambda(l2);
1090  for (NucleonPolarization l1 : NucleonPolarizationIterator() )
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  {
1159  for (NucleonPolarization l2 : NucleonPolarizationIterator() )
1160  {
1161  for (NucleonPolarization l1 : NucleonPolarizationIterator() )
1162  {
1163  // Eq. 18 ref. 1
1164  xsec +=
1165  TMath::Power(
1166  std::abs(
1167  C_L_minus*(
1168  Hbkg(Current::VECTOR, BosonPolarization::LEFT, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::LEFT, l2, l1) +
1169  C1*sum1( BosonPolarization::LEFT, l2, l1) + C3*sum3( BosonPolarization::LEFT, l2, l1)
1170  ) +
1171  C_R_minus*(
1172  Hbkg(Current::VECTOR, BosonPolarization::RIGHT, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::RIGHT, l2, l1) +
1173  C1*sum1( BosonPolarization::RIGHT, l2, l1) + C3*sum3( BosonPolarization::RIGHT, l2, l1)
1174  ) +
1175  C_S_minus*(
1176  Hbkg(Current::VECTOR, BosonPolarization::MINUS0, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::MINUS0, l2, l1) +
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*(
1185  Hbkg(Current::VECTOR, BosonPolarization::LEFT, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::LEFT, l2, l1) +
1186  C1*sum1( BosonPolarization::LEFT, l2, l1) + C3*sum3( BosonPolarization::LEFT, l2, l1)
1187  ) +
1188  C_R_plus*(
1189  Hbkg(Current::VECTOR, BosonPolarization::RIGHT, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::RIGHT, l2, l1) +
1190  C1*sum1( BosonPolarization::RIGHT, l2, l1) + C3*sum3( BosonPolarization::RIGHT, l2, l1)
1191  ) +
1192  C_S_plus*(
1193  Hbkg(Current::VECTOR, BosonPolarization::PLUS0, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::PLUS0, l2, l1) +
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  {
1203  for (NucleonPolarization l2 : NucleonPolarizationIterator() )
1204  {
1205  for (NucleonPolarization l1 : NucleonPolarizationIterator() )
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
1282  P_Fermi = utils::nuclear::FermiMomentumForIsoscalarNucleonParametrization(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 }
genie::SppChannel::FromInteraction
static SppChannel_t FromInteraction(const Interaction *interaction)
Definition: SppChannel.h:402
genie::Kinematics::W
double W(bool selected=false) const
Definition: Kinematics.cxx:157
genie::RSHelicityAmpl::AmpMinus3
double AmpMinus3(void) const
Definition: RSHelicityAmpl.h:57
genie::ProcessInfo::IsWeakCC
bool IsWeakCC(void) const
Definition: ProcessInfo.cxx:203
genie::pdg::IsNeutrino
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
genie::utils::mec::J
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
genie::units::g
static constexpr double g
Definition: Units.h:144
genie::FKR::Rminus
double Rminus
Definition: FKR.h:50
genie::utils::kinematics::Q2
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
genie::Target::HitNucPdg
int HitNucPdg(void) const
Definition: Target.cxx:304
genie::MKSPPPXSec2020::BosonPolarizationIterator
Iterator< BosonPolarization, BosonPolarization::LEFT, BosonPolarization::PLUS0 > BosonPolarizationIterator
Definition: MKSPPPXSec2020.h:154
genie::FKR::Ra
double Ra
Definition: FKR.h:42
genie::Target::A
int A(void) const
Definition: Target.h:70
genie::RSHelicityAmplModelI::Compute
virtual const RSHelicityAmpl & Compute(Resonance_t res, const FKR &fkr) const =0
genie::RSHelicityAmpl::AmpPlus3
double AmpPlus3(void) const
Definition: RSHelicityAmpl.h:58
pFATAL
#define pFATAL
Definition: Messenger.h:56
genie::units::s
static constexpr double s
Definition: Units.h:95
genie::SppChannel::Isospin1Coefficients
static double Isospin1Coefficients(SppChannel_t channel)
Definition: SppChannel.h:317
genie::FermiMomentumTablePool::Instance
static FermiMomentumTablePool * Instance(void)
Definition: FermiMomentumTablePool.cxx:62
genie::Kinematics
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
genie::MKSPPPXSec2020::fSin2Wein
double fSin2Wein
sin^2(Weingberg angle)
Definition: MKSPPPXSec2020.h:384
genie::FKR::Lamda
double Lamda
Definition: FKR.h:37
genie::utils::res::Mass
double Mass(Resonance_t res)
resonance mass (GeV)
Definition: BaryonResUtils.cxx:436
genie::FKR::R
double R
Definition: FKR.h:45
genie::FermiMomentumTable
A table of Fermi momentum constants.
Definition: FermiMomentumTable.h:33
genie::utils::res::Width
double Width(Resonance_t res)
resonance width (GeV)
Definition: BaryonResUtils.cxx:463
genie::MKSPPPXSec2020::fHAmplModelCC
const RSHelicityAmplModelI * fHAmplModelCC
Definition: MKSPPPXSec2020.h:369
genie::utils::res::AngularMom
int AngularMom(Resonance_t res)
Definition: BaryonResUtils.cxx:622
genie::MKSPPPXSec2020::fUseRFGParametrization
bool fUseRFGParametrization
Use parametrization for fermi momentum insted of table?
Definition: MKSPPPXSec2020.h:393
genie::MKSPPPXSec2020::fMa2
double fMa2
(axial mass)^2
Definition: MKSPPPXSec2020.h:378
genie::MKSPPPXSec2020::fXSecScaleNC
double fXSecScaleNC
External NC xsec scaling factor.
Definition: MKSPPPXSec2020.h:387
genie::MKSPPPXSec2020::fCA50
double fCA50
FKR parameter Zeta.
Definition: MKSPPPXSec2020.h:376
genie::SppChannel_t
enum genie::ESppChannel SppChannel_t
genie::utils::res::Cjsgn_plus
int Cjsgn_plus(Resonance_t res)
Definition: BaryonResUtils.cxx:652
genie::utils::kinematics::W
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
genie::RSHelicityAmpl::Amp0Plus
double Amp0Plus(void) const
Definition: RSHelicityAmpl.h:60
genie::units::fm3
static constexpr double fm3
Definition: Units.h:77
genie::FKR::Tv
double Tv
Definition: FKR.h:38
genie::MKSPPPXSec2020::fResList
BaryonResList fResList
Definition: MKSPPPXSec2020.h:419
genie::RSHelicityAmpl
A class holding the Rein-Sehgal&#39;s helicity amplitudes.
Definition: RSHelicityAmpl.h:40
genie::pdg::IsProton
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
genie::ProcessInfo::IsWeakNC
bool IsWeakNC(void) const
Definition: ProcessInfo.cxx:208
genie::FermiMomentumTablePool
Singleton class to load &amp; serve tables of Fermi momentum constants.
Definition: FermiMomentumTablePool.h:34
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::SppChannel::InitStateNucleon
static int InitStateNucleon(SppChannel_t channel)
Definition: SppChannel.h:103
genie::units::A
static constexpr double A
Definition: Units.h:74
genie::FermiMomentumTablePool::GetTable
const FermiMomentumTable * GetTable(string name)
Definition: FermiMomentumTablePool.cxx:76
genie::KinePhaseSpace::AsString
static string AsString(KinePhaseSpace_t kps)
Definition: KinePhaseSpace.h:85
genie::ProcessInfo
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
genie::MKSPPPXSec2020::Lambda
int Lambda(NucleonPolarization l) const
Definition: MKSPPPXSec2020.cxx:1327
genie::FKR::T
double T
Definition: FKR.h:46
genie::FKR::Rv
double Rv
Definition: FKR.h:39
genie::MKSPPPXSec2020::fHAmplModelNCp
const RSHelicityAmplModelI * fHAmplModelNCp
Definition: MKSPPPXSec2020.h:370
genie::pdg::IsAntiNeutrino
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:118
genie::Target
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
genie::InitialState::ProbePdg
int ProbePdg(void) const
Definition: InitialState.h:64
genie::MKSPPPXSec2020::fRho770Mass
double fRho770Mass
Mass of rho(770) meson.
Definition: MKSPPPXSec2020.h:412
genie::RSHelicityAmpl::AmpMinus1
double AmpMinus1(void) const
return helicity amplitude
Definition: RSHelicityAmpl.h:55
genie::kPdgPiP
const int kPdgPiP
Definition: PDGCodes.h:158
genie::kPdgPi0
const int kPdgPi0
Definition: PDGCodes.h:160
genie::MKSPPPXSec2020::fMv2
double fMv2
(vector mass)^2
Definition: MKSPPPXSec2020.h:379
genie::utils::res::OrbitalAngularMom
int OrbitalAngularMom(Resonance_t res)
orbital angular momentum
Definition: BaryonResUtils.cxx:531
genie::RSHelicityAmpl::Amp0Minus
double Amp0Minus(void) const
Definition: RSHelicityAmpl.h:59
genie::MKSPPPXSec2020::fHAmplModelNCn
const RSHelicityAmplModelI * fHAmplModelNCn
Definition: MKSPPPXSec2020.h:371
genie::Target::Z
int Z(void) const
Definition: Target.h:68
genie::MKSPPPXSec2020::FA0
double FA0
Axial coupling (value of axial form factor at Q2=0)
Definition: MKSPPPXSec2020.h:399
genie::Kinematics::GetKV
double GetKV(KineVar_t kv) const
Definition: Kinematics.cxx:323
genie::MKSPPPXSec2020::f_pi
double f_pi
Constant for pion-nucleon interaction.
Definition: MKSPPPXSec2020.h:398
genie::RSHelicityAmplModelI
Pure abstract base class. Defines the RSHelicityAmplModelI interface.
Definition: RSHelicityAmplModelI.h:27
genie::QELFormFactors::xiF2V
double xiF2V(void) const
Get the computed form factor xi*F2V.
Definition: QELFormFactors.h:56
genie::QELFormFactors::Calculate
void Calculate(const Interaction *interaction)
Compute the form factors for the input interaction using the attached model.
Definition: QELFormFactors.cxx:49
genie::utils::res::Dsgn
int Dsgn(Resonance_t res)
Definition: BaryonResUtils.cxx:685
genie::SppChannel::Isospin3Coefficients
static double Isospin3Coefficients(SppChannel_t channel)
Definition: SppChannel.h:280
genie::MKSPPPXSec2020::fKFTable
string fKFTable
Table of Fermi momentum (kF) constants for various nuclei.
Definition: MKSPPPXSec2020.h:392
genie::FKR::C
double C
Definition: FKR.h:44
genie::utils::nuclear::FermiMomentumForIsoscalarNucleonParametrization
double FermiMomentumForIsoscalarNucleonParametrization(const Target &target)
Definition: NuclearUtils.cxx:501
genie::FKR::Tplus
double Tplus
Definition: FKR.h:47
genie::PDGLibrary::Instance
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
genie::MKSPPPXSec2020::fCv3
double fCv3
GV calculation coeffs.
Definition: MKSPPPXSec2020.h:380
genie::MKSPPPXSec2020::fUsePauliBlocking
bool fUsePauliBlocking
Account for Pauli blocking?
Definition: MKSPPPXSec2020.h:394
genie::FKR::B
double B
Definition: FKR.h:43
genie::PDGLibrary
Singleton class to load &amp; serve a TDatabasePDG.
Definition: PDGLibrary.h:35
genie::FKR::Rplus
double Rplus
Definition: FKR.h:49
genie::kIAssumeFreeNucleon
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
genie::MKSPPPXSec2020::fOmega
double fOmega
FKR parameter Omega.
Definition: MKSPPPXSec2020.h:377
genie::MKSPPPXSec2020::fVud
double fVud
|Vud| (magnitude ud-element of CKM-matrix)
Definition: MKSPPPXSec2020.h:385
genie::FKR::Tminus
double Tminus
Definition: FKR.h:48
genie::RSHelicityAmpl::AmpPlus1
double AmpPlus1(void) const
Definition: RSHelicityAmpl.h:56
genie::pdg::IonPdgCode
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:71
genie::MKSPPPXSec2020::PhaseFactor
int PhaseFactor(BosonPolarization lk, NucleonPolarization l1, NucleonPolarization l2) const
Definition: MKSPPPXSec2020.cxx:1337
genie::kPdgPiM
const int kPdgPiM
Definition: PDGCodes.h:159
genie::utils::kinematics::Jacobian
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:130
genie::MKSPPPXSec2020::fUseAuthorCode
bool fUseAuthorCode
Use author code?
Definition: MKSPPPXSec2020.h:415
genie::kPdgProton
const int kPdgProton
Definition: PDGCodes.h:81
genie::QELFormFactors::F1V
double F1V(void) const
Get the computed form factor F1V.
Definition: QELFormFactors.h:53
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
genie::FermiMomentumTable::FindClosestKF
double FindClosestKF(int target_pdgc, int nucleon_pdgc) const
Definition: FermiMomentumTable.cxx:45
genie::Kinematics::Q2
double Q2(bool selected=false) const
Definition: Kinematics.cxx:125
genie::MKSPPPXSec2020::fFormFactors
QELFormFactors fFormFactors
Quasielastic form facors for background contribution.
Definition: MKSPPPXSec2020.h:396
genie::InitialState::Tgt
const Target & Tgt(void) const
Definition: InitialState.h:66
genie::MKSPPPXSec2020::fXSecScaleCC
double fXSecScaleCC
External CC xsec scaling factor.
Definition: MKSPPPXSec2020.h:386
genie::MKSPPPXSec2020::fEMFormFactors
QELFormFactors fEMFormFactors
Electromagnetic form facors for background contribution.
Definition: MKSPPPXSec2020.h:397
genie::MKSPPPXSec2020::ValidProcess
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
Definition: MKSPPPXSec2020.cxx:1312
genie::SppChannel::BranchingRatio
static double BranchingRatio(Resonance_t res)
Definition: SppChannel.h:357
genie::SppChannel::FinStateIsospin
static int FinStateIsospin(SppChannel_t channel)
Definition: SppChannel.h:211
genie::gAbortingInErr
bool gAbortingInErr
Definition: Messenger.cxx:34
genie::InitialState::ProbeE
double ProbeE(RefFrame_t rf) const
Definition: InitialState.cxx:384
genie::kPdgNeutron
const int kPdgNeutron
Definition: PDGCodes.h:83
genie::QELFormFactors::FA
double FA(void) const
Get the computed form factor FA.
Definition: QELFormFactors.h:59
genie::MKSPPPXSec2020::ValidKinematics
bool ValidKinematics(const Interaction *interaction) const
Is the input kinematical point a physically allowed one?
Definition: MKSPPPXSec2020.cxx:1342
genie::FKR::S
double S
Definition: FKR.h:40
genie::FKR::Ta
double Ta
Definition: FKR.h:41
genie::utils::res::Isospin
int Isospin(Resonance_t res)
Definition: BaryonResUtils.cxx:591
genie::utils::res::ResonanceIndex
int ResonanceIndex(Resonance_t res)
resonance idx, quark model / SU(6)
Definition: BaryonResUtils.cxx:561
genie::InitialState
Initial State information.
Definition: InitialState.h:48
genie::MKSPPPXSec2020::NucleonPolarizationIterator
Iterator< NucleonPolarization, NucleonPolarization::MINUS, NucleonPolarization::PLUS > NucleonPolarizationIterator
Definition: MKSPPPXSec2020.h:155

Member Data Documentation

double genie::MKSPPPXSec2020::f_pi
private

Constant for pion-nucleon interaction.

Definition at line 398 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

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().

double genie::MKSPPPXSec2020::fBkgV0
private

Definition at line 411 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

double genie::MKSPPPXSec2020::fBkgV1
private

Definition at line 410 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

double genie::MKSPPPXSec2020::fBkgV2
private

Definition at line 409 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

double genie::MKSPPPXSec2020::fBkgV3
private

Definition at line 408 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

double genie::MKSPPPXSec2020::fBkgVWmax
private

Definition at line 407 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

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().

double genie::MKSPPPXSec2020::fCA50
private

FKR parameter Zeta.

Definition at line 376 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

double genie::MKSPPPXSec2020::fCv3
private

GV calculation coeffs.

Definition at line 380 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

double genie::MKSPPPXSec2020::fCv4
private

Definition at line 381 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

double genie::MKSPPPXSec2020::fCv51
private

Definition at line 382 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

double genie::MKSPPPXSec2020::fCv52
private

Definition at line 383 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

const ELFormFactorsModelI* genie::MKSPPPXSec2020::fElFFModel
private

Elastic form facors model for background contribution.

Definition at line 388 of file MKSPPPXSec2020.h.

QELFormFactors genie::MKSPPPXSec2020::fEMFormFactors
mutableprivate

Electromagnetic form facors for background contribution.

Definition at line 397 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

const QELFormFactorsModelI* genie::MKSPPPXSec2020::fEMFormFactorsModel
private

Electromagnetic form factors model for background contribution.

Definition at line 390 of file MKSPPPXSec2020.h.

Referenced by LoadConfig().

double genie::MKSPPPXSec2020::fFermiConstant
private

Definition at line 375 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

FKR genie::MKSPPPXSec2020::fFKR
mutableprivate

Definition at line 368 of file MKSPPPXSec2020.h.

Referenced by XSec().

QELFormFactors genie::MKSPPPXSec2020::fFormFactors
mutableprivate

Quasielastic form facors for background contribution.

Definition at line 396 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

const QELFormFactorsModelI* genie::MKSPPPXSec2020::fFormFactorsModel
private

Quasielastic form facors model for background contribution.

Definition at line 389 of file MKSPPPXSec2020.h.

Referenced by LoadConfig().

const RSHelicityAmplModelI* genie::MKSPPPXSec2020::fHAmplModelCC
private

Definition at line 369 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

const RSHelicityAmplModelI* genie::MKSPPPXSec2020::fHAmplModelNCn
private

Definition at line 371 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

const RSHelicityAmplModelI* genie::MKSPPPXSec2020::fHAmplModelNCp
private

Definition at line 370 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

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().

double genie::MKSPPPXSec2020::fMa2
private

(axial mass)^2

Definition at line 378 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

double genie::MKSPPPXSec2020::fMv2
private

(vector mass)^2

Definition at line 379 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

double genie::MKSPPPXSec2020::fOmega
private

FKR parameter Omega.

Definition at line 377 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

BaryonResList genie::MKSPPPXSec2020::fResList
private

Definition at line 419 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

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().

double genie::MKSPPPXSec2020::fRho770Mass
private

Mass of rho(770) meson.

Definition at line 412 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

double genie::MKSPPPXSec2020::fSin2Wein
private

sin^2(Weingberg angle)

Definition at line 384 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

bool genie::MKSPPPXSec2020::fUseAuthorCode
private

Use author code?

Definition at line 415 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

bool genie::MKSPPPXSec2020::fUsePauliBlocking
private

Account for Pauli blocking?

Definition at line 394 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

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().

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().

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().

const XSecIntegratorI* genie::MKSPPPXSec2020::fXSecIntegrator
private

Definition at line 417 of file MKSPPPXSec2020.h.

Referenced by Integral(), and LoadConfig().

double genie::MKSPPPXSec2020::fXSecScaleCC
private

External CC xsec scaling factor.

Definition at line 386 of file MKSPPPXSec2020.h.

Referenced by LoadConfig(), and XSec().

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:
Generated on Wed Jul 17 2024 14:47:02 for GENIEGenerator by   doxygen 1.8.5