genie::ReinSehgalRESPXSec Class Reference

Computes the double differential cross section for resonance electro- or neutrino-production according to the Rein-Sehgal model. More...

#include <ReinSehgalRESPXSec.h>

Inheritance diagram for genie::ReinSehgalRESPXSec:

Collaboration diagram for genie::ReinSehgalRESPXSec:

Public Member Functions
	ReinSehgalRESPXSec ()
	ReinSehgalRESPXSec (string config)
virtual	~ReinSehgalRESPXSec ()
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...
void	Configure (const Registry &config)
void	Configure (string config)
Public Member Functions inherited from genie::XSecAlgorithmI
virtual	~XSecAlgorithmI ()
virtual bool	ValidKinematics (const Interaction *i) const
	Is the input kinematical point a physically allowed one? More...
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 Member Functions
void	LoadConfig (void)

Private Attributes
FKR	fFKR
const RSHelicityAmplModelI *	fHAmplModelCC
const RSHelicityAmplModelI *	fHAmplModelNCp
const RSHelicityAmplModelI *	fHAmplModelNCn
const RSHelicityAmplModelI *	fHAmplModelEMp
const RSHelicityAmplModelI *	fHAmplModelEMn
bool	fWghtBW
	weight with resonance breit-wigner? More...
bool	fNormBW
	normalize resonance breit-wigner to 1? More...
double	fZeta
	FKR parameter Zeta. More...
double	fOmega
	FKR parameter Omega. More...
double	fMa2
	(axial mass)^2 More...
double	fMv2
	(vector mass)^2 More...
double	fSin48w
	sin^4(Weingberg angle) More...
double	fVud2
	|Vud|^2(square of magnitude ud-element of CKM-matrix) More...
bool	fUsingDisResJoin
	use a DIS/RES joining scheme? More...
bool	fUsingNuTauScaling
	use NeuGEN nutau xsec reduction factors? More...
double	fWcut
	apply DIS/RES joining scheme < Wcut More...
double	fN2ResMaxNWidths
	limits allowed phase space for n=2 res More...
double	fN0ResMaxNWidths
	limits allowed phase space for n=0 res More...
double	fGnResMaxNWidths
	limits allowed phase space for other res 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...
Spline *	fNuTauRdSpl
	xsec reduction spline for nu_tau More...
Spline *	fNuTauBarRdSpl
	xsec reduction spline for nu_tau_bar More...
double	fXSecScaleCC
	external CC xsec scaling factor More...
double	fXSecScaleNC
	external NC xsec scaling factor More...
double	fXSecScaleEM
	external EM xsec scaling factor More...
const XSecIntegratorI *	fXSecIntegrator

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

Computes the double differential cross section for resonance electro- or neutrino-production according to the Rein-Sehgal model.

      The computed cross section is the d^2 xsec/ dQ^2 dW \n

      where \n
        \li \c Q^2 : momentum transfer ^ 2
        \li \c W   : invariant mass of the final state hadronic system

      Is a concrete implementation of the XSecAlgorithmI interface.

References:

D.Rein and L.M.Sehgal, Neutrino Excitation of Baryon Resonances and Single Pion Production, Ann.Phys.133, 79 (1981)

Author

Costas Andreopoulos <c.andreopoulos cern.ch> University of Liverpool

Created:

May 05, 2004

License:

Copyright (c) 2003-2024, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org

Definition at line 43 of file ReinSehgalRESPXSec.h.

Constructor & Destructor Documentation

ReinSehgalRESPXSec::ReinSehgalRESPXSec

(

)

Definition at line 42 of file ReinSehgalRESPXSec.cxx.

References fNuTauBarRdSpl, and fNuTauRdSpl.

                                       :
XSecAlgorithmI("genie::ReinSehgalRESPXSec")
{
  fNuTauRdSpl    = 0;
  fNuTauBarRdSpl = 0;
}

ReinSehgalRESPXSec::ReinSehgalRESPXSec

(

string

config

)

Definition at line 49 of file ReinSehgalRESPXSec.cxx.

References fNuTauBarRdSpl, and fNuTauRdSpl.

                                                    :
XSecAlgorithmI("genie::ReinSehgalRESPXSec", config)
{
  fNuTauRdSpl    = 0;
  fNuTauBarRdSpl = 0;
}

ReinSehgalRESPXSec::~ReinSehgalRESPXSec ( )

virtual

Definition at line 56 of file ReinSehgalRESPXSec.cxx.

References fNuTauBarRdSpl, and fNuTauRdSpl.

{
  if(fNuTauRdSpl)    delete fNuTauRdSpl;
  if(fNuTauBarRdSpl) delete fNuTauBarRdSpl;
}

Member Function Documentation

void ReinSehgalRESPXSec::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 414 of file ReinSehgalRESPXSec.cxx.

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

{
  Algorithm::Configure(config);
  this->LoadConfig();
}

void ReinSehgalRESPXSec::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 420 of file ReinSehgalRESPXSec.cxx.

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

{
  Algorithm::Configure(config);
  this->LoadConfig();
}

double ReinSehgalRESPXSec::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 381 of file ReinSehgalRESPXSec.cxx.

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

{
  double xsec = fXSecIntegrator->Integrate(this,interaction);
  return xsec;
}

void ReinSehgalRESPXSec::LoadConfig ( void  )

private

Definition at line 426 of file ReinSehgalRESPXSec.cxx.

References fGnResMaxNWidths, fHAmplModelCC, fHAmplModelEMn, fHAmplModelEMp, fHAmplModelNCn, fHAmplModelNCp, fKFTable, fMa2, fMv2, fN0ResMaxNWidths, fN2ResMaxNWidths, fNormBW, fNuTauBarRdSpl, fNuTauRdSpl, fOmega, fSin48w, fUsePauliBlocking, fUseRFGParametrization, fUsingDisResJoin, fUsingNuTauScaling, fVud2, fWcut, fWghtBW, fXSecIntegrator, fXSecScaleCC, fXSecScaleEM, fXSecScaleNC, fZeta, genie::AlgFactory::GetAlgorithm(), genie::Algorithm::GetParam(), genie::Algorithm::GetParamDef(), genie::AlgFactory::Instance(), LOG, pNOTICE, and genie::Algorithm::SubAlg().

Referenced by Configure().

{
  // Cross section scaling factors
  this->GetParam( "RES-CC-XSecScale", fXSecScaleCC ) ;
  this->GetParam( "RES-NC-XSecScale", fXSecScaleNC ) ;
  this->GetParam( "RES-EM-XSecScale", fXSecScaleEM ) ;

  this->GetParam( "RES-Zeta", fZeta ) ;
  this->GetParam( "RES-Omega", fOmega ) ;

  double ma, mv ;
  this->GetParam( "RES-Ma", ma ) ;
  this->GetParam( "RES-Mv", mv ) ;
  fMa2 = TMath::Power(ma,2);
  fMv2 = TMath::Power(mv,2);

  this->GetParamDef( "BreitWignerWeight", fWghtBW, true ) ;
  this->GetParamDef( "BreitWignerNorm",   fNormBW, true);

  double thw ;
  this->GetParam( "WeinbergAngle", thw ) ;
  fSin48w = TMath::Power( TMath::Sin(thw), 4 );
  double Vud;
  this->GetParam("CKM-Vud", Vud );
  fVud2 = TMath::Power( Vud, 2 );
  this->GetParam("FermiMomentumTable", fKFTable);
  this->GetParam("RFG-UseParametrization", fUseRFGParametrization);
  this->GetParam("UsePauliBlockingForRES", fUsePauliBlocking);

  // Load all the sub-algorithms needed

  fHAmplModelCC     = 0;
  fHAmplModelNCp    = 0;
  fHAmplModelNCn    = 0;
  fHAmplModelEMp    = 0;
  fHAmplModelEMn    = 0;

  AlgFactory * algf = AlgFactory::Instance();

  fHAmplModelCC  = dynamic_cast<const RSHelicityAmplModelI *> (
      algf->GetAlgorithm("genie::RSHelicityAmplModelCC","Default"));
  fHAmplModelNCp = dynamic_cast<const RSHelicityAmplModelI *> (
      algf->GetAlgorithm("genie::RSHelicityAmplModelNCp","Default"));
  fHAmplModelNCn = dynamic_cast<const RSHelicityAmplModelI *> (
      algf->GetAlgorithm("genie::RSHelicityAmplModelNCn","Default"));
  fHAmplModelEMp = dynamic_cast<const RSHelicityAmplModelI *> (
      algf->GetAlgorithm("genie::RSHelicityAmplModelEMp","Default"));
  fHAmplModelEMn = dynamic_cast<const RSHelicityAmplModelI *> (
      algf->GetAlgorithm("genie::RSHelicityAmplModelEMn","Default"));

  assert( fHAmplModelCC  );
  assert( fHAmplModelNCp );
  assert( fHAmplModelNCn );
  assert( fHAmplModelEMp );
  assert( fHAmplModelEMn );

  // Use algorithm within a DIS/RES join scheme. If yes get Wcut
  this->GetParam( "UseDRJoinScheme", fUsingDisResJoin ) ;
  fWcut = 999999;
  if(fUsingDisResJoin) {
    this->GetParam( "Wcut", fWcut ) ;
  }

  // NeuGEN limits in the allowed resonance phase space:
  // W < min{ Wmin(physical), (res mass) + x * (res width) }
  // It limits the integration area around the peak and avoids the
  // problem with huge xsec increase at low Q2 and high W.
  // In correspondence with Hugh, Rein said that the underlying problem
  // are unphysical assumptions in the model.
  this->GetParamDef( "MaxNWidthForN2Res", fN2ResMaxNWidths, 2.0 ) ;
  this->GetParamDef( "MaxNWidthForN0Res", fN0ResMaxNWidths, 6.0 ) ;
  this->GetParamDef( "MaxNWidthForGNRes", fGnResMaxNWidths, 4.0 ) ;

  // NeuGEN reduction factors for nu_tau: a gross estimate of the effect of
  // neglected form factors in the R/S model
  this->GetParamDef( "UseNuTauScalingFactors", fUsingNuTauScaling, true ) ;
  if(fUsingNuTauScaling) {
     if(fNuTauRdSpl)    delete fNuTauRdSpl;
     if(fNuTauBarRdSpl) delete fNuTauBarRdSpl;

     assert( std::getenv( "GENIE") );
     string base = std::getenv( "GENIE") ;

     string filename = base + "/data/evgen/rein_sehgal/res/nutau_xsec_scaling_factors.dat";
     LOG("ReinSehgalRes", pNOTICE)
                << "Loading nu_tau xsec reduction spline from: " << filename;
     fNuTauRdSpl = new Spline(filename);

     filename = base + "/data/evgen/rein_sehgal/res/nutaubar_xsec_scaling_factors.dat";
     LOG("ReinSehgalRes", pNOTICE)
           << "Loading bar{nu_tau} xsec reduction spline from: " << filename;
     fNuTauBarRdSpl = new Spline(filename);
  }

  // load the differential cross section integrator
  fXSecIntegrator =
      dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
  assert(fXSecIntegrator);
}

bool ReinSehgalRESPXSec::ValidProcess ( const Interaction *  i ) const

virtual

Can this cross section algorithm handle the input process?

Implements genie::XSecAlgorithmI.

Definition at line 387 of file ReinSehgalRESPXSec.cxx.

References genie::Interaction::ExclTag(), genie::Target::HitNucPdg(), genie::Interaction::InitState(), genie::pdg::IsChargedLepton(), genie::ProcessInfo::IsEM(), genie::pdg::IsNeutralLepton(), genie::pdg::IsNeutron(), genie::pdg::IsProton(), genie::ProcessInfo::IsResonant(), genie::ProcessInfo::IsWeak(), genie::kISkipProcessChk, genie::XclsTag::KnownResonance(), genie::InitialState::ProbePdg(), genie::Interaction::ProcInfo(), and genie::InitialState::Tgt().

Referenced by XSec().

{
  if(interaction->TestBit(kISkipProcessChk)) return true;

  const InitialState & init_state = interaction->InitState();
  const ProcessInfo &  proc_info  = interaction->ProcInfo();
  const XclsTag &      xcls       = interaction->ExclTag();

  if(!proc_info.IsResonant()) return false;
  if(!xcls.KnownResonance())  return false;

  int  hitnuc = init_state.Tgt().HitNucPdg();
  bool is_pn = (pdg::IsProton(hitnuc) || pdg::IsNeutron(hitnuc));

  if (!is_pn) return false;

  int  probe   = init_state.ProbePdg();
  bool is_weak = proc_info.IsWeak();
  bool is_em   = proc_info.IsEM();
  bool nu_weak = (pdg::IsNeutralLepton(probe) && is_weak);
  bool l_em    = (pdg::IsChargedLepton(probe) && is_em  );

  if (!nu_weak && !l_em) return false;

  return true;
}

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

virtual

Compute the cross section for the input interaction.

Implements genie::XSecAlgorithmI.

Definition at line 62 of file ReinSehgalRESPXSec.cxx.

References genie::Target::A(), genie::units::A, genie::utils::res::AsString(), genie::Interaction::AsString(), genie::FKR::B, genie::utils::bwfunc::BreitWignerL(), genie::utils::bwfunc::BreitWignerLGamma(), genie::utils::res::BWNorm(), genie::FKR::C, genie::RSHelicityAmplModelI::Compute(), genie::Spline::Evaluate(), genie::Interaction::ExclTag(), genie::utils::nuclear::FermiMomentumForIsoscalarNucleonParametrization(), fFKR, fGnResMaxNWidths, fHAmplModelCC, fHAmplModelEMn, fHAmplModelEMp, fHAmplModelNCn, fHAmplModelNCp, genie::FermiMomentumTable::FindClosestKF(), fKFTable, fMa2, fMv2, fN0ResMaxNWidths, fN2ResMaxNWidths, fNormBW, fNuTauBarRdSpl, fNuTauRdSpl, fOmega, fSin48w, fUsePauliBlocking, fUseRFGParametrization, fUsingDisResJoin, fUsingNuTauScaling, fVud2, fWcut, fWghtBW, fXSecScaleCC, fXSecScaleEM, fXSecScaleNC, fZeta, genie::FermiMomentumTablePool::GetTable(), genie::Target::HitNucMass(), genie::Target::HitNucPdg(), genie::FermiMomentumTablePool::Instance(), genie::pdg::IonPdgCode(), genie::pdg::IsAntiNeutrino(), genie::pdg::IsAntiNuTau(), genie::utils::res::IsDelta(), genie::ProcessInfo::IsEM(), genie::pdg::IsNegChargedLepton(), genie::pdg::IsNeutrino(), genie::pdg::IsNeutron(), genie::pdg::IsNuTau(), genie::pdg::IsPosChargedLepton(), genie::pdg::IsProton(), genie::ProcessInfo::IsWeakCC(), genie::ProcessInfo::IsWeakNC(), genie::utils::mec::J(), genie::utils::kinematics::Jacobian(), genie::constants::kAem2, genie::constants::kGF2, genie::kIAssumeFreeNucleon, genie::constants::kPi, genie::constants::kPi2, genie::constants::kPionMass2, genie::kPSWQ2fE, genie::kRfHitNucRest, genie::constants::kSqrt2, genie::FKR::Lamda, LOG, genie::utils::res::Mass(), genie::utils::res::OrbitalAngularMom(), pDEBUG, pINFO, genie::InitialState::ProbeE(), genie::InitialState::ProbePdg(), genie::Kinematics::q2(), genie::utils::kinematics::Q2(), genie::FKR::R, genie::FKR::Ra, genie::XclsTag::Resonance(), genie::utils::res::ResonanceIndex(), genie::FKR::Rminus, genie::FKR::Rplus, genie::FKR::Rv, genie::FKR::S, genie::FKR::T, genie::FKR::Ta, genie::InitialState::Tgt(), genie::FKR::Tminus, genie::FKR::Tplus, genie::FKR::Tv, genie::XSecAlgorithmI::ValidKinematics(), ValidProcess(), genie::Kinematics::W(), genie::utils::kinematics::W(), genie::utils::res::Width(), and genie::Target::Z().

{
  if(! this -> ValidProcess    (interaction) ) return 0.;
  if(! this -> ValidKinematics (interaction) ) return 0.;

  const InitialState & init_state = interaction -> InitState();
  const ProcessInfo &  proc_info  = interaction -> ProcInfo();
  const Target & target = init_state.Tgt();

  // Get kinematical parameters
  const Kinematics & kinematics = interaction -> Kine();
  double W  = kinematics.W();
  double q2 = kinematics.q2();

  // Under the DIS/RES joining scheme, xsec(RES)=0 for W>=Wcut
  if(fUsingDisResJoin) {
    if(W>=fWcut) {
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
       LOG("ReinSehgalRes", pDEBUG)
         << "RES/DIS Join Scheme: XSec[RES, W=" << W
         << " >= Wcut=" << fWcut << "] = 0";
#endif
       return 0;
    }
  }

  // Get the input baryon resonance
  Resonance_t resonance = interaction->ExclTag().Resonance();
  string      resname   = utils::res::AsString(resonance);
  bool        is_delta  = utils::res::IsDelta (resonance);

  // Get the neutrino, hit nucleon & weak current
  int  nucpdgc   = target.HitNucPdg();
  int  probepdgc = init_state.ProbePdg();
  bool is_nu     = pdg::IsNeutrino         (probepdgc);
  bool is_nubar  = pdg::IsAntiNeutrino     (probepdgc);
  bool is_lplus  = pdg::IsPosChargedLepton (probepdgc);
  bool is_lminus = pdg::IsNegChargedLepton (probepdgc);
  bool is_p      = pdg::IsProton  (nucpdgc);
  bool is_n      = pdg::IsNeutron (nucpdgc);
  bool is_CC     = proc_info.IsWeakCC();
  bool is_NC     = proc_info.IsWeakNC();
  bool is_EM     = proc_info.IsEM();

  if(is_CC && !is_delta) {
    if((is_nu && is_p) || (is_nubar && is_n)) return 0;
  }

  // Get baryon resonance parameters
  int    IR  = utils::res::ResonanceIndex    (resonance);
  int    LR  = utils::res::OrbitalAngularMom (resonance);
  double MR  = utils::res::Mass              (resonance);
  double WR  = utils::res::Width             (resonance);
  double NR  = fNormBW?utils::res::BWNorm    (resonance,fN0ResMaxNWidths,fN2ResMaxNWidths,fGnResMaxNWidths):1;

  // Following NeuGEN, avoid problems with underlying unphysical
  // model assumptions by restricting the allowed W phase space
  // around the resonance peak
  if (fNormBW) {
        if      (W > MR + fN0ResMaxNWidths * WR && IR==0) return 0.;
        else if (W > MR + fN2ResMaxNWidths * WR && IR==2) return 0.;
        else if (W > MR + fGnResMaxNWidths * WR)          return 0.;
  }

  // Compute auxiliary & kinematical factors
  double E      = init_state.ProbeE(kRfHitNucRest);
  double Mnuc   = target.HitNucMass();
  double W2     = TMath::Power(W,    2);
  double Mnuc2  = TMath::Power(Mnuc, 2);
  double k      = 0.5 * (W2 - Mnuc2)/Mnuc;
  double v      = k - 0.5 * q2/Mnuc;
  double v2     = TMath::Power(v, 2);
  double Q2     = v2 - q2;
  double Q      = TMath::Sqrt(Q2);
  double Eprime = E - v;
  double U      = 0.5 * (E + Eprime + Q) / E;
  double V      = 0.5 * (E + Eprime - Q) / E;
  double U2     = TMath::Power(U, 2);
  double V2     = TMath::Power(V, 2);
  double UV     = U*V;

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("ReinSehgalRes", pDEBUG)
     << "Kinematical params V = " << V << ", U = " << U;
#endif

  // Calculate the Feynman-Kislinger-Ravndall parameters

  double Go  = TMath::Power(1 - 0.25 * q2/Mnuc2, 0.5-IR);
  double GV  = Go * TMath::Power( 1./(1-q2/fMv2), 2);
  double GA  = Go * TMath::Power( 1./(1-q2/fMa2), 2);

  if(is_EM) {
    GA = 0.; // zero the axial term for EM scattering
  }

  double d      = TMath::Power(W+Mnuc,2.) - q2;
  double sq2omg = TMath::Sqrt(2./fOmega);
  double nomg   = IR * fOmega;
  double mq_w   = Mnuc*Q/W;

  fFKR.Lamda  = sq2omg * mq_w;
  fFKR.Tv     = GV / (3.*W*sq2omg);
  fFKR.Rv     = kSqrt2 * mq_w*(W+Mnuc)*GV / d;
  fFKR.S      = (-q2/Q2) * (3*W*Mnuc + q2 - Mnuc2) * GV / (6*Mnuc2);
  fFKR.Ta     = (2./3.) * (fZeta/sq2omg) * mq_w * GA / d;
  fFKR.Ra     = (kSqrt2/6.) * fZeta * (GA/W) * (W+Mnuc + 2*nomg*W/d );
  fFKR.B      = fZeta/(3.*W*sq2omg) * (1 + (W2-Mnuc2+q2)/ d) * GA;
  fFKR.C      = fZeta/(6.*Q) * (W2 - Mnuc2 + nomg*(W2-Mnuc2+q2)/d) * (GA/Mnuc);
  fFKR.R      = fFKR.Rv;
  fFKR.Rplus  = - (fFKR.Rv + fFKR.Ra);
  fFKR.Rminus = - (fFKR.Rv - fFKR.Ra);
  fFKR.T      = fFKR.Tv;
  fFKR.Tplus  = - (fFKR.Tv + fFKR.Ta);
  fFKR.Tminus = - (fFKR.Tv - fFKR.Ta);

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("FKR", pDEBUG)
     << "FKR params for RES = " << resname << " : " << fFKR;
#endif

  // Calculate the Rein-Sehgal Helicity Amplitudes

  const RSHelicityAmplModelI * hamplmod = 0;
  if(is_CC) {
    hamplmod = fHAmplModelCC;
  }
  else
  if(is_NC) {
    if (is_p) { hamplmod = fHAmplModelNCp;}
    else      { hamplmod = fHAmplModelNCn;}
  }
  else
  if(is_EM) {
    if (is_p) { hamplmod = fHAmplModelEMp;}
    else      { hamplmod = fHAmplModelEMn;}
  }
  assert(hamplmod);

  const RSHelicityAmpl & hampl = hamplmod->Compute(resonance, fFKR);

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("RSHAmpl", pDEBUG)
     << "Helicity Amplitudes for RES = " << resname << " : " << hampl;
#endif

  double g2 = kGF2;
  if(is_CC) g2 = kGF2*fVud2;
  // For EM interaction replace  G_{Fermi} with :
  // a_{em} * pi / ( sqrt(2) * sin^2(theta_weinberg) * Mass_{W}^2 }
  // See C.Quigg, Gauge Theories of the Strong, Weak and E/M Interactions,
  // ISBN 0-8053-6021-2, p.112 (6.3.57)
  // Also, take int account that the photon propagator is 1/p^2 but the
  // W propagator is 1/(p^2-Mass_{W}^2), so weight the EM case with
  // Mass_{W}^4 / q^4
  // So, overall:
  // G_{Fermi}^2 --> a_{em}^2 * pi^2 / (2 * sin^4(theta_weinberg) * q^{4})
  //
  if(is_EM) {
    double q4 = q2*q2;
    g2 = kAem2 * kPi2 / (2.0 * fSin48w * q4);
  }

  // Compute the cross section

  double sig0 = 0.125*(g2/kPi)*(-q2/Q2)*(W/Mnuc);
  double scLR = W/Mnuc;
  double scS  = (Mnuc/W)*(-Q2/q2);
  double sigL = scLR* (hampl.Amp2Plus3 () + hampl.Amp2Plus1 ());
  double sigR = scLR* (hampl.Amp2Minus3() + hampl.Amp2Minus1());
  double sigS = scS * (hampl.Amp20Plus () + hampl.Amp20Minus());

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("ReinSehgalRes", pDEBUG) << "sig_{0} = " << sig0;
  LOG("ReinSehgalRes", pDEBUG) << "sig_{L} = " << sigL;
  LOG("ReinSehgalRes", pDEBUG) << "sig_{R} = " << sigR;
  LOG("ReinSehgalRes", pDEBUG) << "sig_{S} = " << sigS;
#endif

  double xsec = 0.0;
  if (is_nu || is_lminus) {
     xsec = sig0*(V2*sigR + U2*sigL + 2*UV*sigS);
  }
  else
  if (is_nubar || is_lplus) {
     xsec = sig0*(U2*sigR + V2*sigL + 2*UV*sigS);
  }
  xsec = TMath::Max(0.,xsec);

  double mult = 1.0;
  if(is_CC && is_delta) {
    if((is_nu && is_p) || (is_nubar && is_n)) mult=3.0;
  }
  xsec *= mult;

  // Check whether the cross section is to be weighted with a
  // Breit-Wigner distribution (default: true)
  double bw = 1.0;
  if(fWghtBW) {
     //different Delta photon decay branch
     if(is_delta){
     bw = utils::bwfunc::BreitWignerLGamma(W,LR,MR,WR,NR);
     }
     else{
     bw = utils::bwfunc::BreitWignerL(W,LR,MR,WR,NR);
     }
  }
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
     LOG("ReinSehgalRes", pDEBUG)
       << "BreitWigner(RES=" << resname << ", W=" << W << ") = " << bw;
#endif
  xsec *= bw;

  // Apply NeuGEN nutau cross section reduction factors
  double rf = 1.0;
  Spline * spl = 0;
  if (is_CC && fUsingNuTauScaling) {
    if      (pdg::IsNuTau    (probepdgc)) spl = fNuTauRdSpl;
    else if (pdg::IsAntiNuTau(probepdgc)) spl = fNuTauBarRdSpl;

    if(spl) {
      if(E <spl->XMax()) rf = spl->Evaluate(E);
    }
  }
  xsec *= rf;

  // Apply given scaling factor
  double xsec_scale = 1.;
  if      (is_CC) { xsec_scale = fXSecScaleCC; }
  else if (is_NC) { xsec_scale = fXSecScaleNC; }
  else if (is_EM) { xsec_scale = fXSecScaleEM; }
  xsec *= xsec_scale;

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("ReinSehgalRes", pINFO)
    << "\n d2xsec/dQ2dW"  << "[" << interaction->AsString()
          << "](W=" << W << ", q2=" << q2 << ", E=" << E << ") = " << xsec;
#endif

  // The algorithm computes d^2xsec/dWdQ2
  // Check whether variable tranformation is needed
  if(kps!=kPSWQ2fE) {
    double J = utils::kinematics::Jacobian(interaction,kPSWQ2fE,kps);
    xsec *= J;
  }

  // If requested return the free nucleon xsec even for input nuclear tgt
  if( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;


  int Z = target.Z();
  int A = target.A();
  int N = A-Z;

  // Take into account the number of scattering centers in the target
  int NNucl = (is_p) ? Z : N;

  xsec*=NNucl; // nuclear xsec (no nuclear suppression factor)

  if (fUsePauliBlocking && A!=1)
  {
     // Calculation of Pauli blocking according references:
     //
     //     [1] S.L. Adler,  S. Nussinov,  and  E.A.  Paschos,  "Nuclear
     //         charge exchange corrections to leptonic pion  production
     //         in  the (3,3) resonance  region,"  Phys. Rev. D 9 (1974)
     //         2125-2143 [Erratum Phys. Rev. D 10 (1974) 1669].
     //     [2] J.Y. Yu, "Neutrino interactions and  nuclear  effects in
     //         oscillation experiments and the  nonperturbative disper-
     //         sive  sector in strong (quasi-)abelian  fields,"  Ph. D.
     //         Thesis, Dortmund U., Dortmund, 2002 (unpublished).
     //     [3] E.A. Paschos, J.Y. Yu,  and  M. Sakuda,  "Neutrino  pro-
     //         duction  of  resonances,"  Phys. Rev. D 69 (2004) 014013
     //         [arXiv: hep-ph/0308130].

     double P_Fermi = 0.0;

     // Maximum value of Fermi momentum of target nucleon (GeV)
     if (A<6 || !fUseRFGParametrization)
     {
         // Look up the Fermi momentum for this target
         FermiMomentumTablePool * kftp = FermiMomentumTablePool::Instance();
         const FermiMomentumTable * kft = kftp->GetTable(fKFTable);
         P_Fermi = kft->FindClosestKF(pdg::IonPdgCode(A, Z), nucpdgc);
     }
     else {
        // Define the Fermi momentum for this target
        P_Fermi = utils::nuclear::FermiMomentumForIsoscalarNucleonParametrization(target);
        // Correct the Fermi momentum for the struck nucleon
        if(is_p) { P_Fermi *= TMath::Power( 2.*Z/A, 1./3); }
        else     { P_Fermi *= TMath::Power( 2.*N/A, 1./3); }
     }

     double FactorPauli_RES = 1.0;

     double k0 = 0., q = 0., q0 = 0.;

     if (P_Fermi > 0.)
     {
        k0 = (W2-Mnuc2-Q2)/(2*W);
        k = TMath::Sqrt(k0*k0+Q2);                  // previous value of k is overridden
        q0 = (W2-Mnuc2+kPionMass2)/(2*W);
        q = TMath::Sqrt(q0*q0-kPionMass2);
     }

     if (2*P_Fermi < k-q)
        FactorPauli_RES = 1.0;
     if (2*P_Fermi >= k+q)
        FactorPauli_RES = ((3*k*k+q*q)/(2*P_Fermi)-(5*TMath::Power(k,4)+TMath::Power(q,4)+10*k*k*q*q)/(40*TMath::Power(P_Fermi,3)))/(2*k);
     if (2*P_Fermi >= k-q && 2*P_Fermi <= k+q)
        FactorPauli_RES = ((q+k)*(q+k)-4*P_Fermi*P_Fermi/5-TMath::Power(k-q, 3)/(2*P_Fermi)+TMath::Power(k-q, 5)/(40*TMath::Power(P_Fermi, 3)))/(4*q*k);

     xsec *= FactorPauli_RES;
  }

  return xsec;
}

Member Data Documentation

FKR genie::ReinSehgalRESPXSec::fFKR

mutableprivate

Definition at line 64 of file ReinSehgalRESPXSec.h.

Referenced by XSec().

double genie::ReinSehgalRESPXSec::fGnResMaxNWidths

private

limits allowed phase space for other res

Definition at line 86 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

const RSHelicityAmplModelI* genie::ReinSehgalRESPXSec::fHAmplModelCC

private

Definition at line 66 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

const RSHelicityAmplModelI* genie::ReinSehgalRESPXSec::fHAmplModelEMn

private

Definition at line 70 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

const RSHelicityAmplModelI* genie::ReinSehgalRESPXSec::fHAmplModelEMp

private

Definition at line 69 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

const RSHelicityAmplModelI* genie::ReinSehgalRESPXSec::fHAmplModelNCn

private

Definition at line 68 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

const RSHelicityAmplModelI* genie::ReinSehgalRESPXSec::fHAmplModelNCp

private

Definition at line 67 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

string genie::ReinSehgalRESPXSec::fKFTable

private

table of Fermi momentum (kF) constants for various nuclei

Definition at line 87 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::ReinSehgalRESPXSec::fMa2

private

(axial mass)^2

Definition at line 77 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::ReinSehgalRESPXSec::fMv2

private

(vector mass)^2

Definition at line 78 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::ReinSehgalRESPXSec::fN0ResMaxNWidths

private

limits allowed phase space for n=0 res

Definition at line 85 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::ReinSehgalRESPXSec::fN2ResMaxNWidths

private

limits allowed phase space for n=2 res

Definition at line 84 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

bool genie::ReinSehgalRESPXSec::fNormBW

private

normalize resonance breit-wigner to 1?

Definition at line 74 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

Spline* genie::ReinSehgalRESPXSec::fNuTauBarRdSpl

private

xsec reduction spline for nu_tau_bar

Definition at line 91 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), ReinSehgalRESPXSec(), XSec(), and ~ReinSehgalRESPXSec().

Spline* genie::ReinSehgalRESPXSec::fNuTauRdSpl

private

xsec reduction spline for nu_tau

Definition at line 90 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), ReinSehgalRESPXSec(), XSec(), and ~ReinSehgalRESPXSec().

double genie::ReinSehgalRESPXSec::fOmega

private

FKR parameter Omega.

Definition at line 76 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::ReinSehgalRESPXSec::fSin48w

private

sin^4(Weingberg angle)

Definition at line 79 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

bool genie::ReinSehgalRESPXSec::fUsePauliBlocking

private

account for Pauli blocking?

Definition at line 89 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

bool genie::ReinSehgalRESPXSec::fUseRFGParametrization

private

use parametrization for fermi momentum insted of table?

Definition at line 88 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

bool genie::ReinSehgalRESPXSec::fUsingDisResJoin

private

use a DIS/RES joining scheme?

Definition at line 81 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

bool genie::ReinSehgalRESPXSec::fUsingNuTauScaling

private

use NeuGEN nutau xsec reduction factors?

Definition at line 82 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::ReinSehgalRESPXSec::fVud2

private

|Vud|^2(square of magnitude ud-element of CKM-matrix)

Definition at line 80 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::ReinSehgalRESPXSec::fWcut

private

apply DIS/RES joining scheme < Wcut

Definition at line 83 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

bool genie::ReinSehgalRESPXSec::fWghtBW

private

weight with resonance breit-wigner?

Definition at line 73 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

const XSecIntegratorI* genie::ReinSehgalRESPXSec::fXSecIntegrator

private

Definition at line 96 of file ReinSehgalRESPXSec.h.

Referenced by Integral(), and LoadConfig().

double genie::ReinSehgalRESPXSec::fXSecScaleCC

private

external CC xsec scaling factor

Definition at line 92 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::ReinSehgalRESPXSec::fXSecScaleEM

private

external EM xsec scaling factor

Definition at line 94 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::ReinSehgalRESPXSec::fXSecScaleNC

private

external NC xsec scaling factor

Definition at line 93 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::ReinSehgalRESPXSec::fZeta

private

FKR parameter Zeta.

Definition at line 75 of file ReinSehgalRESPXSec.h.

Referenced by LoadConfig(), and XSec().

---

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

• ReinSehgalRESPXSec.h
• ReinSehgalRESPXSec.cxx