Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
genie::BSKLNBaseRESPXSec2014 Class Reference

Base class for the Berger-Sehgal and the Kuzmin, Lyubushkin, Naumov resonance models, implemented as modifications to the Rein-Sehgal model. More...

#include <BSKLNBaseRESPXSec2014.h>

Inheritance diagram for genie::BSKLNBaseRESPXSec2014:
genie::XSecAlgorithmI genie::Algorithm genie::BergerSehgalRESPXSec2014 genie::KuzminLyubushkinNaumovRESPXSec2014

Public Member Functions

virtual ~BSKLNBaseRESPXSec2014 ()
 
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 RegistryGetConfig (void) const
 
RegistryGetOwnedConfig (void)
 
virtual const AlgIdId (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 AlgorithmSubAlg (const RgKey &registry_key) const
 
void AdoptConfig (void)
 
void AdoptSubstructure (void)
 
virtual void Print (ostream &stream) const
 Print algorithm info. More...
 

Protected Member Functions

 BSKLNBaseRESPXSec2014 (string name)
 
 BSKLNBaseRESPXSec2014 (string name, string config)
 
void LoadConfig (void)
 
- Protected Member Functions inherited from genie::XSecAlgorithmI
 XSecAlgorithmI ()
 
 XSecAlgorithmI (string name)
 
 XSecAlgorithmI (string name, string config)
 
- Protected Member Functions inherited from genie::Algorithm
 Algorithm ()
 
 Algorithm (string name)
 
 Algorithm (string name, string config)
 
void Initialize (void)
 
void DeleteConfig (void)
 
void DeleteSubstructure (void)
 
RegistryExtractLocalConfig (const Registry &in) const
 
RegistryExtractLowerConfig (const Registry &in, const string &alg_key) const
 Split an incoming configuration Registry into a block valid for the sub-algo identified by alg_key. 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
 
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

FKR fFKR
 
const RSHelicityAmplModelIfHAmplModelCC
 
const RSHelicityAmplModelIfHAmplModelNCp
 
const RSHelicityAmplModelIfHAmplModelNCn
 
const RSHelicityAmplModelIfHAmplModelEMp
 
const RSHelicityAmplModelIfHAmplModelEMn
 
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...
 
double fXSecScaleCC
 external CC xsec scaling factor More...
 
double fXSecScaleNC
 external NC xsec scaling factor More...
 
bool fKLN
 
bool fBRS
 
bool fGA
 
bool fGV
 
const XSecIntegratorIfXSecIntegrator
 
- 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< boolfOwnerships
 ownership for every registry in fConfVect More...
 
AlgStatus_t fStatus
 algorithm execution status More...
 
AlgMapfOwnedSubAlgMp
 local pool for owned sub-algs (taken out of the factory pool) More...
 

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)
 

Detailed Description

Base class for the Berger-Sehgal and the Kuzmin, Lyubushkin, Naumov resonance models, implemented as modifications to the Rein-Sehgal model.

Berger, Sehgal Phys. Rev. D76, 113004 (2007)
Kuzmin, Lyubushkin, Naumov Mod. Phys. Lett. A19 (2004) 2815
D.Rein and L.M.Sehgal, Neutrino Excitation of Baryon Resonances and Single Pion Production, Ann.Phys.133, 79 (1981)
Modifications based on a MiniBooNE tune courtesy of J. Nowak, S.Dytman

Author
Steve Dytman University of Pittsburgh

Jarek Nowak University of Lancaster

Gabe Perdue Fermilab

Costas Andreopoulos <constantinos.andreopoulos cern.ch> University of Liverpool & STFC Rutherford Appleton Laboratory

Sep 15, 2015

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

Definition at line 48 of file BSKLNBaseRESPXSec2014.h.

Constructor & Destructor Documentation

BSKLNBaseRESPXSec2014::~BSKLNBaseRESPXSec2014 ( )
virtual

Definition at line 69 of file BSKLNBaseRESPXSec2014.cxx.

70 {
71 
72 }
BSKLNBaseRESPXSec2014::BSKLNBaseRESPXSec2014 ( string  name)
protected

Definition at line 57 of file BSKLNBaseRESPXSec2014.cxx.

57  :
59 {
60 
61 }
static QCString name
Definition: declinfo.cpp:673
BSKLNBaseRESPXSec2014::BSKLNBaseRESPXSec2014 ( string  name,
string  config 
)
protected

Definition at line 63 of file BSKLNBaseRESPXSec2014.cxx.

63  :
65 {
66 
67 }
static QCString name
Definition: declinfo.cpp:673
static Config * config
Definition: config.cpp:1054

Member Function Documentation

void BSKLNBaseRESPXSec2014::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 749 of file BSKLNBaseRESPXSec2014.cxx.

750 {
751  Algorithm::Configure(config);
752  this->LoadConfig();
753 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void BSKLNBaseRESPXSec2014::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 755 of file BSKLNBaseRESPXSec2014.cxx.

756 {
758  this->LoadConfig();
759 }
static Config * config
Definition: config.cpp:1054
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
double BSKLNBaseRESPXSec2014::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 716 of file BSKLNBaseRESPXSec2014.cxx.

717 {
718  double xsec = fXSecIntegrator->Integrate(this,interaction);
719  return xsec;
720 }
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
const XSecIntegratorI * fXSecIntegrator
void BSKLNBaseRESPXSec2014::LoadConfig ( void  )
protected

Definition at line 761 of file BSKLNBaseRESPXSec2014.cxx.

762 {
763  // Cross section scaling factors
764  this->GetParam( "RES-CC-XSecScale", fXSecScaleCC ) ;
765  this->GetParam( "RES-NC-XSecScale", fXSecScaleNC ) ;
766 
767  // Load all configuration data or set defaults
768 
769  this->GetParam( "RES-Zeta" , fZeta ) ;
770  this->GetParam( "RES-Omega" , fOmega ) ;
771  this->GetParam( "minibooneGA", fGA ) ;
772  this->GetParam( "minibooneGV", fGV ) ;
773 
774  double ma, mv ;
775  this->GetParam( "RES-Ma", ma ) ;
776  this->GetParam( "RES-Mv", mv ) ;
777  fMa2 = TMath::Power(ma,2);
778  fMv2 = TMath::Power(mv,2);
779 
780  this->GetParamDef( "BreitWignerWeight", fWghtBW, true ) ;
781  this->GetParamDef( "BreitWignerNorm", fNormBW, true);
782  double thw ;
783  this->GetParam( "WeinbergAngle", thw ) ;
784  fSin48w = TMath::Power( TMath::Sin(thw), 4 );
785  double Vud;
786  this->GetParam("CKM-Vud", Vud );
787  fVud2 = TMath::Power( Vud, 2 );
788  this->GetParam("FermiMomentumTable", fKFTable);
789  this->GetParam("RFG-UseParametrization", fUseRFGParametrization);
790  this->GetParam("UsePauliBlockingForRES", fUsePauliBlocking);
791 
792  // Load all the sub-algorithms needed
793 
794  fHAmplModelCC = 0;
795  fHAmplModelNCp = 0;
796  fHAmplModelNCn = 0;
797  fHAmplModelEMp = 0;
798  fHAmplModelEMn = 0;
799 
800  AlgFactory * algf = AlgFactory::Instance();
801 
802  fHAmplModelCC = dynamic_cast<const RSHelicityAmplModelI *> (
803  algf->GetAlgorithm("genie::RSHelicityAmplModelCC","Default"));
804  fHAmplModelNCp = dynamic_cast<const RSHelicityAmplModelI *> (
805  algf->GetAlgorithm("genie::RSHelicityAmplModelNCp","Default"));
806  fHAmplModelNCn = dynamic_cast<const RSHelicityAmplModelI *> (
807  algf->GetAlgorithm("genie::RSHelicityAmplModelNCn","Default"));
808  fHAmplModelEMp = dynamic_cast<const RSHelicityAmplModelI *> (
809  algf->GetAlgorithm("genie::RSHelicityAmplModelEMp","Default"));
810  fHAmplModelEMn = dynamic_cast<const RSHelicityAmplModelI *> (
811  algf->GetAlgorithm("genie::RSHelicityAmplModelEMn","Default"));
812 
813  assert( fHAmplModelCC );
814  assert( fHAmplModelNCp );
815  assert( fHAmplModelNCn );
816  assert( fHAmplModelEMp );
817  assert( fHAmplModelEMn );
818 
819  // Use algorithm within a DIS/RES join scheme. If yes get Wcut
820  this->GetParam( "UseDRJoinScheme", fUsingDisResJoin ) ;
821  fWcut = 999999;
822  if(fUsingDisResJoin) {
823  this->GetParam( "Wcut", fWcut ) ;
824  }
825 
826  // NeuGEN limits in the allowed resonance phase space:
827  // W < min{ Wmin(physical), (res mass) + x * (res width) }
828  // It limits the integration area around the peak and avoids the
829  // problem with huge xsec increase at low Q2 and high W.
830  // In correspondence with Hugh, Rein said that the underlying problem
831  // are unphysical assumptions in the model.
832  this->GetParamDef( "MaxNWidthForN2Res", fN2ResMaxNWidths, 2.0 ) ;
833  this->GetParamDef( "MaxNWidthForN0Res", fN0ResMaxNWidths, 6.0 ) ;
834  this->GetParamDef( "MaxNWidthForGNRes", fGnResMaxNWidths, 4.0 ) ;
835 
836  // Load the differential cross section integrator
838  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
839  assert(fXSecIntegrator);
840 }
bool fNormBW
normalize resonance breit-wigner to 1?
string fKFTable
table of Fermi momentum (kF) constants for various nuclei
double fOmega
FKR parameter Omega.
bool fUsingDisResJoin
use a DIS/RES joining scheme?
double fXSecScaleNC
external NC xsec scaling factor
Cross Section Integrator Interface.
double fN0ResMaxNWidths
limits allowed phase space for n=0 res
const RSHelicityAmplModelI * fHAmplModelEMp
const RSHelicityAmplModelI * fHAmplModelCC
double fVud2
|Vud|^2(square of magnitude ud-element of CKM-matrix)
bool fWghtBW
weight with resonance breit-wigner?
double fXSecScaleCC
external CC xsec scaling factor
const Algorithm * GetAlgorithm(const AlgId &algid)
Definition: AlgFactory.cxx:75
bool fUsePauliBlocking
account for Pauli blocking?
double fWcut
apply DIS/RES joining scheme < Wcut
const RSHelicityAmplModelI * fHAmplModelEMn
Pure abstract base class. Defines the RSHelicityAmplModelI interface.
double fGnResMaxNWidths
limits allowed phase space for other res
bool fUseRFGParametrization
use parametrization for fermi momentum insted of table?
const RSHelicityAmplModelI * fHAmplModelNCp
static AlgFactory * Instance()
Definition: AlgFactory.cxx:64
double fSin48w
sin^4(Weingberg angle)
const RSHelicityAmplModelI * fHAmplModelNCn
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
The GENIE Algorithm Factory.
Definition: AlgFactory.h:39
double fN2ResMaxNWidths
limits allowed phase space for n=2 res
double fZeta
FKR parameter Zeta.
const XSecIntegratorI * fXSecIntegrator
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345
bool BSKLNBaseRESPXSec2014::ValidProcess ( const Interaction i) const
virtual

Can this cross section algorithm handle the input process?

Implements genie::XSecAlgorithmI.

Definition at line 722 of file BSKLNBaseRESPXSec2014.cxx.

723 {
724  if(interaction->TestBit(kISkipProcessChk)) return true;
725 
726  const InitialState & init_state = interaction->InitState();
727  const ProcessInfo & proc_info = interaction->ProcInfo();
728  const XclsTag & xcls = interaction->ExclTag();
729 
730  if(!proc_info.IsResonant()) return false;
731  if(!xcls.KnownResonance()) return false;
732 
733  int hitnuc = init_state.Tgt().HitNucPdg();
734  bool is_pn = (pdg::IsProton(hitnuc) || pdg::IsNeutron(hitnuc));
735 
736  if (!is_pn) return false;
737 
738  int probe = init_state.ProbePdg();
739  bool is_weak = proc_info.IsWeak();
740  bool is_em = proc_info.IsEM();
741  bool nu_weak = (pdg::IsNeutralLepton(probe) && is_weak);
742  bool l_em = (pdg::IsChargedLepton(probe) && is_em );
743 
744  if (!nu_weak && !l_em) return false;
745 
746  return true;
747 }
bool IsResonant(void) const
Definition: ProcessInfo.cxx:94
bool IsWeak(void) const
int HitNucPdg(void) const
Definition: Target.cxx:304
bool KnownResonance(void) const
Definition: XclsTag.h:68
bool IsChargedLepton(int pdgc)
Definition: PDGUtils.cxx:98
Contains minimal information for tagging exclusive processes.
Definition: XclsTag.h:39
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:338
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:333
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
int ProbePdg(void) const
Definition: InitialState.h:64
bool IsEM(void) const
bool IsNeutralLepton(int pdgc)
Definition: PDGUtils.cxx:92
const Target & Tgt(void) const
Definition: InitialState.h:66
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
Initial State information.
Definition: InitialState.h:48
double BSKLNBaseRESPXSec2014::XSec ( const Interaction i,
KinePhaseSpace_t  k 
) const
virtual

Compute the cross section for the input interaction.

Implements genie::XSecAlgorithmI.

Definition at line 74 of file BSKLNBaseRESPXSec2014.cxx.

76 {
77  if(! this -> ValidProcess (interaction) ) return 0.;
78  if(! this -> ValidKinematics (interaction) ) return 0.;
79 
80  const InitialState & init_state = interaction -> InitState();
81  const ProcessInfo & proc_info = interaction -> ProcInfo();
82  const Target & target = init_state.Tgt();
83 
84  // Get kinematical parameters
85  const Kinematics & kinematics = interaction -> Kine();
86  double W = kinematics.W();
87  double q2 = kinematics.q2();
88  double costh = kinematics.FSLeptonP4().CosTheta();
89 
90  // Under the DIS/RES joining scheme, xsec(RES)=0 for W>=Wcut
91  if(fUsingDisResJoin) {
92  if(W>=fWcut) {
93 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
94  LOG("BSKLNBaseRESPXSec2014", pDEBUG)
95  << "RES/DIS Join Scheme: XSec[RES, W=" << W
96  << " >= Wcut=" << fWcut << "] = 0";
97 #endif
98  return 0;
99  }
100  }
101 
102  // Get the input baryon resonance
103  Resonance_t resonance = interaction->ExclTag().Resonance();
104  string resname = utils::res::AsString(resonance);
105  bool is_delta = utils::res::IsDelta (resonance);
106 
107  // Get the neutrino, hit nucleon & weak current
108  int nucpdgc = target.HitNucPdg();
109  int probepdgc = init_state.ProbePdg();
110  bool is_nu = pdg::IsNeutrino (probepdgc);
111  bool is_nubar = pdg::IsAntiNeutrino (probepdgc);
112  bool is_lplus = pdg::IsPosChargedLepton (probepdgc);
113  bool is_lminus = pdg::IsNegChargedLepton (probepdgc);
114  bool is_p = pdg::IsProton (nucpdgc);
115  bool is_n = pdg::IsNeutron (nucpdgc);
116  bool is_CC = proc_info.IsWeakCC();
117  bool is_NC = proc_info.IsWeakNC();
118  bool is_EM = proc_info.IsEM();
119 
120  // bool new_GV = fGA; //JN
121  // bool new_GA = fGV; //JN
122 
123 
124  if(is_CC && !is_delta) {
125  if((is_nu && is_p) || (is_nubar && is_n)) return 0;
126  }
127 
128  // Get baryon resonance parameters
129  int IR = utils::res::ResonanceIndex (resonance);
130  int LR = utils::res::OrbitalAngularMom (resonance);
131  double MR = utils::res::Mass (resonance);
132  double WR = utils::res::Width (resonance);
134 
135  // Following NeuGEN, avoid problems with underlying unphysical
136  // model assumptions by restricting the allowed W phase space
137  // around the resonance peak
138  if (fNormBW) {
139  if (W > MR + fN0ResMaxNWidths * WR && IR==0) return 0.;
140  else if (W > MR + fN2ResMaxNWidths * WR && IR==2) return 0.;
141  else if (W > MR + fGnResMaxNWidths * WR) return 0.;
142  }
143 
144  // Compute auxiliary & kinematical factors
145  double E = init_state.ProbeE(kRfHitNucRest);
146  double Mnuc = target.HitNucMass();
147  double W2 = TMath::Power(W, 2);
148  double Mnuc2 = TMath::Power(Mnuc, 2);
149  double k = 0.5 * (W2 - Mnuc2)/Mnuc;
150  double v = k - 0.5 * q2/Mnuc;
151  double v2 = TMath::Power(v, 2);
152  double Q2 = v2 - q2;
153  double Q = TMath::Sqrt(Q2);
154  double Eprime = E - v;
155  double U = 0.5 * (E + Eprime + Q) / E;
156  double V = 0.5 * (E + Eprime - Q) / E;
157  double U2 = TMath::Power(U, 2);
158  double V2 = TMath::Power(V, 2);
159  double UV = U*V;
160 
161 
162  //JN parameter from the KUZMIN et al.
163 
164  // bool is_RS = true;
165  bool is_KLN = false;
166  if(fKLN && is_CC) is_KLN=true;
167 
168  bool is_BRS = false;
169  if(fBRS && is_CC) is_BRS=true;
170 
171  double ml = interaction->FSPrimLepton()->Mass();
172  double Pl = TMath::Sqrt(Eprime*Eprime - ml*ml);
173 
174  double vstar = (Mnuc*v + q2)/W; //missing W
175  double Qstar = TMath::Sqrt(-q2 + vstar*vstar);
176  double sqrtq2 = TMath::Sqrt(-q2);
177  double a = 1. + 0.5*(W2-q2+Mnuc2)/Mnuc/W;
178 
179  double KNL_Alambda_plus = 0;
180  double KNL_Alambda_minus = 0;
181  double KNL_j0_plus = 0;
182  double KNL_j0_minus = 0;
183  double KNL_jx_plus = 0;
184  double KNL_jx_minus = 0;
185  double KNL_jy_plus = 0;
186  double KNL_jy_minus = 0;
187  double KNL_jz_plus = 0;
188  double KNL_jz_minus = 0;
189  double KNL_Qstar_plus =0;
190  double KNL_Qstar_minus =0;
191 
192  double KNL_K = Q/E/TMath::Sqrt(2*(-q2));
193 
194  double KNL_cL_plus = 0;
195  double KNL_cL_minus = 0;
196 
197  double KNL_cR_plus = 0;
198  double KNL_cR_minus = 0;
199 
200  double KNL_cS_plus = 0;
201  double KNL_cS_minus = 0;
202 
203  double KNL_vstar_plus = 0;
204  double KNL_vstar_minus = 0;
205 
206  if(is_CC && (is_KLN || is_BRS)){
207 
208  LOG("BSKLNBaseRESPXSec2014",pINFO) "costh1="<<costh;
209  costh = (q2 - ml*ml + 2.*E*Eprime)/2./E/Pl;
210  //ml=0;
211  LOG("BSKLNBaseRESPXSec2014",pINFO) "q2="<<q2<< "m2="<<ml*ml<<" 2.*E*Eprime="<<2.*E*Eprime<<" nom="<< (q2 - ml*ml + 2.*E*Eprime)<<" den="<<2.*E*Pl;
212  LOG("BSKLNBaseRESPXSec2014",pINFO) "costh2="<<costh;
213 
214  KNL_Alambda_plus = TMath::Sqrt(E*(Eprime - Pl));
215  KNL_Alambda_minus = TMath::Sqrt(E*(Eprime + Pl));
216  LOG("BSKLNBaseRESPXSec2014",pINFO)
217  << "\n+++++++++++++++++++++++ \n"
218  << "E="<<E << " K= "<<KNL_K << "\n"
219  << "El="<<Eprime<<" Pl="<<Pl<<" ml="<<ml << "\n"
220  << "W="<<W<<" Q="<<Q<<" q2="<<q2 << "\n"
221  << "A-="<<KNL_Alambda_minus<<" A+="<<KNL_Alambda_plus << "\n"
222  << "xxxxxxxxxxxxxxxxxxxxxxx";
223 
224  KNL_j0_plus = KNL_Alambda_plus /W * TMath::Sqrt(1 - costh) * (Mnuc - Eprime - Pl);
225  KNL_j0_minus = KNL_Alambda_minus/W * TMath::Sqrt(1 + costh) * (Mnuc - Eprime + Pl);
226 
227  KNL_jx_plus = KNL_Alambda_plus/ Q * TMath::Sqrt(1 + costh) * (Pl - E);
228  KNL_jx_minus = KNL_Alambda_minus/Q * TMath::Sqrt(1 - costh) * (Pl + E);
229 
230  KNL_jy_plus = KNL_Alambda_plus * TMath::Sqrt(1 + costh);
231  KNL_jy_minus = -KNL_Alambda_minus * TMath::Sqrt(1 - costh);
232 
233  KNL_jz_plus = KNL_Alambda_plus /W/Q * TMath::Sqrt(1 - costh) * ( (E + Pl)*(Mnuc -Eprime) + Pl*( E + 2*E*costh -Pl) );
234  KNL_jz_minus = KNL_Alambda_minus/W/Q * TMath::Sqrt(1 + costh) * ( (E - Pl)*(Mnuc -Eprime) + Pl*( -E + 2*E*costh -Pl) );
235 
236  if (is_nu || is_lminus) {
237  KNL_Qstar_plus = sqrtq2 * KNL_j0_plus / TMath::Sqrt(TMath::Abs(KNL_j0_plus*KNL_j0_plus - KNL_jz_plus*KNL_jz_plus) );
238  KNL_Qstar_minus = sqrtq2 * KNL_j0_minus / TMath::Sqrt(TMath::Abs(KNL_j0_minus*KNL_j0_minus - KNL_jz_minus*KNL_jz_minus) );
239  }
240 
241  else if (is_nubar || is_lplus){
242  KNL_Qstar_plus = sqrtq2 * KNL_j0_minus / TMath::Sqrt(TMath::Abs(KNL_j0_minus*KNL_j0_minus - KNL_jz_minus*KNL_jz_minus) );
243  KNL_Qstar_minus = sqrtq2 * KNL_j0_plus / TMath::Sqrt(TMath::Abs(KNL_j0_plus*KNL_j0_plus - KNL_jz_plus*KNL_jz_plus) );
244  }
245 
246  if (is_nu || is_lminus) {
247  KNL_vstar_plus = sqrtq2 * KNL_jz_plus / TMath::Sqrt(TMath::Abs(KNL_j0_plus*KNL_j0_plus - KNL_jz_plus*KNL_jz_plus) );
248  KNL_vstar_minus = sqrtq2 * KNL_jz_minus / TMath::Sqrt(TMath::Abs(KNL_j0_minus*KNL_j0_minus - KNL_jz_minus*KNL_jz_minus) );
249  }
250  else if (is_nubar || is_lplus) {
251  KNL_vstar_minus = sqrtq2 * KNL_jz_plus / TMath::Sqrt(TMath::Abs(KNL_j0_plus*KNL_j0_plus - KNL_jz_plus*KNL_jz_plus) );
252  KNL_vstar_plus = sqrtq2 * KNL_jz_minus / TMath::Sqrt(TMath::Abs(KNL_j0_minus*KNL_j0_minus - KNL_jz_minus*KNL_jz_minus) );
253  }
254 
255  if(is_nu || is_lminus){
256  KNL_cL_plus = TMath::Sqrt(0.5)* KNL_K * (KNL_jx_plus - KNL_jy_plus);
257  KNL_cL_minus = TMath::Sqrt(0.5)* KNL_K * (KNL_jx_minus - KNL_jy_minus);
258 
259  KNL_cR_plus = -TMath::Sqrt(0.5)* KNL_K * (KNL_jx_plus + KNL_jy_plus);
260  KNL_cR_minus = -TMath::Sqrt(0.5)* KNL_K * (KNL_jx_minus + KNL_jy_minus);
261 
262  KNL_cS_plus = KNL_K * TMath::Sqrt(TMath::Abs(KNL_j0_plus *KNL_j0_plus - KNL_jz_plus *KNL_jz_plus ) );
263  KNL_cS_minus = KNL_K * TMath::Sqrt(TMath::Abs(KNL_j0_minus*KNL_j0_minus - KNL_jz_minus*KNL_jz_minus) );
264  }
265 
266  if (is_nubar || is_lplus) {
267  KNL_cL_plus = -1 * TMath::Sqrt(0.5)* KNL_K * (KNL_jx_minus + KNL_jy_minus);
268  KNL_cL_minus = 1 * TMath::Sqrt(0.5)* KNL_K * (KNL_jx_plus + KNL_jy_plus);
269 
270  KNL_cR_plus = 1 * TMath::Sqrt(0.5)* KNL_K * (KNL_jx_minus - KNL_jy_minus);
271  KNL_cR_minus = -1 * TMath::Sqrt(0.5)* KNL_K * (KNL_jx_plus - KNL_jy_plus);
272 
273  KNL_cS_plus = -1 * KNL_K * TMath::Sqrt(TMath::Abs(KNL_j0_minus*KNL_j0_minus - KNL_jz_minus*KNL_jz_minus) );
274  KNL_cS_minus = 1 * KNL_K * TMath::Sqrt(TMath::Abs(KNL_j0_plus*KNL_j0_plus - KNL_jz_plus*KNL_jz_plus) );
275  }
276  }
277 
278  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"j0-="<<KNL_j0_minus<<" j0+="<<KNL_j0_plus;
279  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"jx-="<<KNL_jx_minus<<" jx+="<<KNL_jx_plus;
280  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"jy-="<<KNL_jy_minus<<" jy+="<<KNL_jy_plus;
281  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"jz-="<<KNL_jz_minus<<" jz+="<<KNL_jz_plus;
282 
283  LOG("BSKLNBaseRESPXSec2014",pINFO) "sqrt2="<<sqrtq2<<" jz+=:"<<KNL_jz_plus<<" j0+="<<KNL_j0_plus<<" denom="<<TMath::Sqrt(TMath::Abs(KNL_j0_plus*KNL_j0_plus - KNL_jz_plus*KNL_jz_plus) );
284 
285  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"vstar-="<<KNL_vstar_minus<<" vstar+="<<KNL_vstar_plus;
286  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"Qstar-="<<KNL_Qstar_minus<<" Qstar+="<<KNL_Qstar_plus;
287 
288 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
289  LOG("BSKLNBaseRESPXSec2014", pDEBUG)
290  << "Kinematical params V = " << V << ", U = " << U;
291 #endif
292 
293  // Calculate the Feynman-Kislinger-Ravndall parameters
294 
295  double Go = TMath::Power(1 - 0.25 * q2/Mnuc2, 0.5-IR);
296  double GV = Go * TMath::Power( 1./(1-q2/fMv2), 2);
297  double GA = Go * TMath::Power( 1./(1-q2/fMa2), 2);
298 
299  if(fGV){
300 
301  LOG("BSKLNBaseRESPXSec2014",pDEBUG) <<"Using new GV";
302  double CV0 = 1./(1-q2/fMv2/4.);
303  double CV3 = 2.13 * CV0 * TMath::Power( 1-q2/fMv2,-2);
304  double CV4 = -1.51 * CV0 * TMath::Power( 1-q2/fMv2,-2);
305  double CV5 = 0.48 * CV0 * TMath::Power( 1-q2/fMv2/0.766, -2);
306 
307  double GV3 = 0.5 / TMath::Sqrt(3) * ( CV3 * (W + Mnuc)/Mnuc
308  + CV4 * (W2 + q2 -Mnuc2)/2./Mnuc2
309  + CV5 * (W2 - q2 -Mnuc2)/2./Mnuc2 );
310 
311  double GV1 = - 0.5 / TMath::Sqrt(3) * ( CV3 * (Mnuc2 -q2 +Mnuc*W)/W/Mnuc
312  + CV4 * (W2 +q2 - Mnuc2)/2./Mnuc2
313  + CV5 * (W2 -q2 - Mnuc2)/2./Mnuc2 );
314 
315  GV = 0.5 * TMath::Power( 1 - q2/(Mnuc + W)/(Mnuc + W), 0.5-IR)
316  * TMath::Sqrt( 3 * GV3*GV3 + GV1*GV1);
317  }
318 
319  if(fGA){
320  LOG("BSKLNBaseRESPXSec2014",pDEBUG) << "Using new GA";
321 
322  double CA5_0 = 1.2;
323  double CA5 = CA5_0 * TMath::Power( 1./(1-q2/fMa2), 2);
324  // GA = 0.5 * TMath::Sqrt(3.) * TMath::Power( 1 - q2/(Mnuc + W)/(Mnuc + W), 0.5-IR) * (1- (W2 +q2 -Mnuc2)/8./Mnuc2) * CA5/fZeta;
325  GA = 0.5 * TMath::Sqrt(3.) * TMath::Power( 1 - q2/(Mnuc + W)/(Mnuc + W), 0.5-IR) * (1- (W2 +q2 -Mnuc2)/8./Mnuc2) * CA5;
326 
327  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"GA= " <<GA << " C5A= " <<CA5;
328  }
329  //JN end of new form factors code
330 
331  if(is_EM) {
332  GA = 0.; // zero the axial term for EM scattering
333  }
334 
335  double d = TMath::Power(W+Mnuc,2.) - q2;
336  double sq2omg = TMath::Sqrt(2./fOmega);
337  double nomg = IR * fOmega;
338  double mq_w = Mnuc*Q/W;
339 
340  fFKR.Lamda = sq2omg * mq_w;
341  fFKR.Tv = GV / (3.*W*sq2omg);
342  fFKR.Rv = kSqrt2 * mq_w*(W+Mnuc)*GV / d;
343  fFKR.S = (-q2/Q2) * (3*W*Mnuc + q2 - Mnuc2) * GV / (6*Mnuc2);
344  fFKR.Ta = (2./3.) * (fZeta/sq2omg) * mq_w * GA / d;
345  fFKR.Ra = (kSqrt2/6.) * fZeta * (GA/W) * (W+Mnuc + 2*nomg*W/d );
346  fFKR.B = fZeta/(3.*W*sq2omg) * (1 + (W2-Mnuc2+q2)/ d) * GA;
347  fFKR.C = fZeta/(6.*Q) * (W2 - Mnuc2 + nomg*(W2-Mnuc2+q2)/d) * (GA/Mnuc);
348  fFKR.R = fFKR.Rv;
349  fFKR.Rplus = - (fFKR.Rv + fFKR.Ra);
350  fFKR.Rminus = - (fFKR.Rv - fFKR.Ra);
351  fFKR.T = fFKR.Tv;
352  fFKR.Tplus = - (fFKR.Tv + fFKR.Ta);
353  fFKR.Tminus = - (fFKR.Tv - fFKR.Ta);
354 
355  //JN KNL
356  double KNL_S_plus = 0;
357  double KNL_S_minus = 0;
358  double KNL_B_plus = 0;
359  double KNL_B_minus = 0;
360  double KNL_C_plus = 0;
361  double KNL_C_minus = 0;
362 
363  if(is_CC && is_KLN){
364  KNL_S_plus = (KNL_vstar_plus*vstar - KNL_Qstar_plus *Qstar )* (Mnuc2 -q2 - 3*W*Mnuc ) * GV / (6*Mnuc2)/Q2; //possibly missing minus sign ()
365  KNL_S_minus = (KNL_vstar_minus*vstar - KNL_Qstar_minus*Qstar )* (Mnuc2 -q2 - 3*W*Mnuc ) * GV / (6*Mnuc2)/Q2;
366 
367  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"KNL S= " <<KNL_S_plus<<"\t"<<KNL_S_minus<<"\t"<<fFKR.S;
368 
369  KNL_B_plus = fZeta/(3.*W*sq2omg)/Qstar * (KNL_Qstar_plus + KNL_vstar_plus *Qstar/a/Mnuc ) * GA;
370  KNL_B_minus = fZeta/(3.*W*sq2omg)/Qstar * (KNL_Qstar_minus + KNL_vstar_minus*Qstar/a/Mnuc ) * GA;
371  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"KNL B= " <<KNL_B_plus<<"\t"<<KNL_B_minus<<"\t"<<fFKR.B;
372 
373  KNL_C_plus = ( (KNL_Qstar_plus*Qstar - KNL_vstar_plus*vstar ) * ( 1./3. + vstar/a/Mnuc)
374  + KNL_vstar_plus*(2./3.*W +q2/a/Mnuc + nomg/3./a/Mnuc) )* fZeta * (GA/2./W/Qstar);
375 
376  KNL_C_minus = ( (KNL_Qstar_minus*Qstar - KNL_vstar_minus*vstar ) * ( 1./3. + vstar/a/Mnuc)
377  + KNL_vstar_minus*(2./3.*W +q2/a/Mnuc + nomg/3./a/Mnuc) )* fZeta * (GA/2./W/Qstar);
378 
379  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"KNL C= "<<KNL_C_plus<<"\t"<<KNL_C_minus<<"\t"<<fFKR.C;
380  }
381  double BRS_S_plus = 0;
382  double BRS_S_minus = 0;
383  double BRS_B_plus = 0;
384  double BRS_B_minus = 0;
385  double BRS_C_plus = 0;
386  double BRS_C_minus = 0;
387 
388 
389  if(is_CC && is_BRS){
390 
391  KNL_S_plus = (KNL_vstar_plus*vstar - KNL_Qstar_plus *Qstar )* (Mnuc2 -q2 - 3*W*Mnuc ) * GV / (6*Mnuc2)/Q2;
392  KNL_S_minus = (KNL_vstar_minus*vstar - KNL_Qstar_minus*Qstar )* (Mnuc2 -q2 - 3*W*Mnuc ) * GV / (6*Mnuc2)/Q2;
393 
394 
395  KNL_B_plus = fZeta/(3.*W*sq2omg)/Qstar * (KNL_Qstar_plus + KNL_vstar_plus *Qstar/a/Mnuc ) * GA;
396  KNL_B_minus = fZeta/(3.*W*sq2omg)/Qstar * (KNL_Qstar_minus + KNL_vstar_minus*Qstar/a/Mnuc ) * GA;
397 
398 
399  KNL_C_plus = ( (KNL_Qstar_plus*Qstar - KNL_vstar_plus*vstar ) * ( 1./3. + vstar/a/Mnuc)
400  + KNL_vstar_plus*(2./3.*W +q2/a/Mnuc + nomg/3./a/Mnuc) )* fZeta * (GA/2./W/Qstar);
401 
402  KNL_C_minus = ( (KNL_Qstar_minus*Qstar - KNL_vstar_minus*vstar ) * ( 1./3. + vstar/a/Mnuc)
403  + KNL_vstar_minus*(2./3.*W +q2/a/Mnuc + nomg/3./a/Mnuc) )* fZeta * (GA/2./W/Qstar);
404 
405  BRS_S_plus = KNL_S_plus;
406  BRS_S_minus = KNL_S_minus;
407  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"BRS S= " <<KNL_S_plus<<"\t"<<KNL_S_minus<<"\t"<<fFKR.S;
408 
409  BRS_B_plus = KNL_B_plus + fZeta*GA/2./W/Qstar*( KNL_Qstar_plus*vstar - KNL_vstar_plus*Qstar)
410  *( 2./3 /sq2omg *(vstar + Qstar*Qstar/Mnuc/a))/(kPionMass2 -q2);
411 
412  BRS_B_minus = KNL_B_minus + fZeta*GA/2./W/Qstar*( KNL_Qstar_minus*vstar - KNL_vstar_minus*Qstar)
413  *( 2./3 /sq2omg *(vstar + Qstar*Qstar/Mnuc/a))/(kPionMass2 -q2);
414  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"BRS B= " <<KNL_B_plus<<"\t"<<KNL_B_minus<<"\t"<<fFKR.B;
415 
416  BRS_C_plus = KNL_C_plus + fZeta*GA/2./W/Qstar*( KNL_Qstar_plus*vstar - KNL_vstar_plus*Qstar)
417  * Qstar*(2./3.*W +q2/Mnuc/a +nomg/3./a/Mnuc)/(kPionMass2 -q2);
418 
419  BRS_C_minus = KNL_C_minus + fZeta*GA/2./W/Qstar*( KNL_Qstar_minus*vstar - KNL_vstar_minus*Qstar)
420  * Qstar*(2./3.*W +q2/Mnuc/a +nomg/3./a/Mnuc)/(kPionMass2 -q2);
421  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"BRS C= " <<KNL_C_plus<<"\t"<<KNL_C_minus<<"\t"<<fFKR.C;
422  }
423 
424 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
425  LOG("FKR", pDEBUG)
426  << "FKR params for RES = " << resname << " : " << fFKR;
427 #endif
428 
429  // Calculate the Rein-Sehgal Helicity Amplitudes
430  double sigL_minus = 0;
431  double sigR_minus = 0;
432  double sigS_minus = 0;
433 
434  double sigL_plus = 0;
435  double sigR_plus = 0;
436  double sigS_plus = 0;
437 
438  const RSHelicityAmplModelI * hamplmod = 0;
439  const RSHelicityAmplModelI * hamplmod_KNL_minus = 0;
440  const RSHelicityAmplModelI * hamplmod_KNL_plus = 0;
441  const RSHelicityAmplModelI * hamplmod_BRS_minus = 0;
442  const RSHelicityAmplModelI * hamplmod_BRS_plus = 0;
443 
444  // These lines were ~ 100 lines below, which means that, for EM interactions, the coefficients below were still calculated using the weak coupling constant - Afro
445  double g2 = kGF2;
446 
447  // For EM interaction replace G_{Fermi} with :
448  // a_{em} * pi / ( sqrt(2) * sin^2(theta_weinberg) * Mass_{W}^2 }
449  // See C.Quigg, Gauge Theories of the Strong, Weak and E/M Interactions,
450  // ISBN 0-8053-6021-2, p.112 (6.3.57)
451  // Also, take int account that the photon propagator is 1/p^2 but the
452  // W propagator is 1/(p^2-Mass_{W}^2), so weight the EM case with
453  // Mass_{W}^4 / q^4
454  // So, overall:
455  // G_{Fermi}^2 --> a_{em}^2 * pi^2 / (2 * sin^4(theta_weinberg) * q^{4})
456  //
457 
458  if(is_EM) {
459  double q4 = q2*q2;
460  g2 = kAem2 * kPi2 / (2.0 * fSin48w * q4);
461  }
462 
463  if(is_CC) g2 = kGF2*fVud2;
464 
465  double sig0 = 0.125*(g2/kPi)*(-q2/Q2)*(W/Mnuc);
466  double scLR = W/Mnuc;
467  double scS = (Mnuc/W)*(-Q2/q2);
468 
469  double sigL =0;
470  double sigR =0;
471  double sigS =0;
472 
473  double sigRSL =0;
474  double sigRSR =0;
475  double sigRSS =0;
476 
477  if(is_CC && !(is_KLN || is_BRS) ) {
478 
479  hamplmod = fHAmplModelCC;
480  }
481  else
482  if(is_NC) {
483  if (is_p) { hamplmod = fHAmplModelNCp;}
484  else { hamplmod = fHAmplModelNCn;}
485  }
486  else
487  if(is_EM) {
488  if (is_p) { hamplmod = fHAmplModelEMp;}
489  else { hamplmod = fHAmplModelEMn;}
490  }
491  else
492  if(is_CC && is_KLN ){
493  fFKR.S = KNL_S_minus; //2 times fFKR.S?
494  fFKR.B = KNL_B_minus;
495  fFKR.C = KNL_C_minus;
496 
497  hamplmod_KNL_minus = fHAmplModelCC;
498 
499  assert(hamplmod_KNL_minus);
500 
501  const RSHelicityAmpl & hampl_KNL_minus = hamplmod_KNL_minus->Compute(resonance, fFKR);
502 
503  sigL_minus = (hampl_KNL_minus.Amp2Plus3 () + hampl_KNL_minus.Amp2Plus1 ());
504  sigR_minus = (hampl_KNL_minus.Amp2Minus3() + hampl_KNL_minus.Amp2Minus1());
505  sigS_minus = (hampl_KNL_minus.Amp20Plus () + hampl_KNL_minus.Amp20Minus());
506 
507 
508  fFKR.S = KNL_S_plus;
509  fFKR.B = KNL_B_plus;
510  fFKR.C = KNL_C_plus;
511  hamplmod_KNL_plus = fHAmplModelCC;
512  assert(hamplmod_KNL_plus);
513 
514  const RSHelicityAmpl & hampl_KNL_plus = hamplmod_KNL_plus->Compute(resonance, fFKR);
515 
516  sigL_plus = (hampl_KNL_plus.Amp2Plus3 () + hampl_KNL_plus.Amp2Plus1 ());
517  sigR_plus = (hampl_KNL_plus.Amp2Minus3() + hampl_KNL_plus.Amp2Minus1());
518  sigS_plus = (hampl_KNL_plus.Amp20Plus () + hampl_KNL_plus.Amp20Minus());
519 
520  }
521  else
522  if(is_CC && is_BRS ){
523  fFKR.S = BRS_S_minus;
524  fFKR.B = BRS_B_minus;
525  fFKR.C = BRS_C_minus;
526 
527  hamplmod_BRS_minus = fHAmplModelCC;
528  assert(hamplmod_BRS_minus);
529 
530  const RSHelicityAmpl & hampl_BRS_minus = hamplmod_BRS_minus->Compute(resonance, fFKR);
531 
532  sigL_minus = (hampl_BRS_minus.Amp2Plus3 () + hampl_BRS_minus.Amp2Plus1 ());
533  sigR_minus = (hampl_BRS_minus.Amp2Minus3() + hampl_BRS_minus.Amp2Minus1());
534  sigS_minus = (hampl_BRS_minus.Amp20Plus () + hampl_BRS_minus.Amp20Minus());
535 
536  fFKR.S = BRS_S_plus;
537  fFKR.B = BRS_B_plus;
538  fFKR.C = BRS_C_plus;
539  hamplmod_BRS_plus = fHAmplModelCC;
540  assert(hamplmod_BRS_plus);
541 
542  const RSHelicityAmpl & hampl_BRS_plus = hamplmod_BRS_plus->Compute(resonance, fFKR);
543 
544  sigL_plus = (hampl_BRS_plus.Amp2Plus3 () + hampl_BRS_plus.Amp2Plus1 ());
545  sigR_plus = (hampl_BRS_plus.Amp2Minus3() + hampl_BRS_plus.Amp2Minus1());
546  sigS_plus = (hampl_BRS_plus.Amp20Plus () + hampl_BRS_plus.Amp20Minus());
547  }
548 
549  // Compute the cross section
550  if(is_KLN || is_BRS) {
551 
552  sigL_minus *= scLR;
553  sigR_minus *= scLR;
554  sigS_minus *= scS;
555  sigL_plus *= scLR;
556  sigR_plus *= scLR;
557  sigS_plus *= scS;
558 
559  LOG("BSKLNBaseRESPXSec2014", pINFO)
560  << "sL,R,S minus = " << sigL_minus << "," << sigR_minus << "," << sigS_minus;
561  LOG("BSKLNBaseRESPXSec2014", pINFO)
562  << "sL,R,S plus = " << sigL_plus << "," << sigR_plus << "," << sigS_plus;
563  }
564  else {
565  assert(hamplmod);
566 
567  const RSHelicityAmpl & hampl = hamplmod->Compute(resonance, fFKR);
568 
569  sigL = scLR* (hampl.Amp2Plus3 () + hampl.Amp2Plus1 ());
570  sigR = scLR* (hampl.Amp2Minus3() + hampl.Amp2Minus1());
571  sigS = scS * (hampl.Amp20Plus () + hampl.Amp20Minus());
572  }
573 
574 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
575  LOG("BSKLNBaseRESPXSec2014", pDEBUG) << "sig_{0} = " << sig0;
576  LOG("BSKLNBaseRESPXSec2014", pDEBUG) << "sig_{L} = " << sigL;
577  LOG("BSKLNBaseRESPXSec2014", pDEBUG) << "sig_{R} = " << sigR;
578  LOG("BSKLNBaseRESPXSec2014", pDEBUG) << "sig_{S} = " << sigS;
579 #endif
580 
581  double xsec = 0.0;
582 
583  if(is_KLN || is_BRS) {
584  xsec = TMath::Power(KNL_cL_minus,2)*sigL_minus + TMath::Power(KNL_cL_plus,2)*sigL_plus
585  + TMath::Power(KNL_cR_minus,2)*sigR_minus + TMath::Power(KNL_cR_plus,2)*sigR_plus
586  + TMath::Power(KNL_cS_minus,2)*sigS_minus + TMath::Power(KNL_cS_plus,2)*sigS_plus;
587  xsec *=sig0;
588 
589  LOG("BSKLNBaseRESPXSec2014",pINFO) << "A-="<<KNL_Alambda_minus<<" A+="<<KNL_Alambda_plus;
590  // protect against sigRSR=sigRSL=sigRSS=0
591  LOG("BSKLNBaseRESPXSec2014",pINFO) <<q2<<"\t"<<xsec<<"\t"<<sig0*(V2*sigR + U2*sigL + 2*UV*sigS)<<"\t"<<xsec/TMath::Max(sig0*(V2*sigRSR + U2*sigRSL + 2*UV*sigRSS),1.0e-100);
592  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"fFKR.B="<<fFKR.B<<" fFKR.C="<<fFKR.C<<" fFKR.S="<<fFKR.S;
593  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"CL-="<<TMath::Power(KNL_cL_minus,2)<<" CL+="<<TMath::Power(KNL_cL_plus,2)<<" U2="<<U2;
594  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"SL-="<<sigL_minus<<" SL+="<<sigL_plus<<" SL="<<sigRSL;
595 
596  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"CR-="<<TMath::Power(KNL_cR_minus,2)<<" CR+="<<TMath::Power(KNL_cR_plus,2)<<" V2="<<V2;
597  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"SR-="<<sigR_minus<<" SR+="<<sigR_plus<<" sR="<<sigRSR;
598 
599  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"CS-="<<TMath::Power(KNL_cS_minus,2)<<" CS+="<<TMath::Power(KNL_cS_plus,2)<<" UV="<<UV;
600  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"SS-="<<sigL_minus<<" SS+="<<sigS_plus<<" sS="<<sigRSS;
601  }
602  else {
603  if (is_nu || is_lminus) {
604  xsec = sig0*(V2*sigR + U2*sigL + 2*UV*sigS);
605  }
606  else
607  if (is_nubar || is_lplus) {
608  xsec = sig0*(U2*sigR + V2*sigL + 2*UV*sigS);
609  }
610  xsec = TMath::Max(0.,xsec);
611  }
612  double mult = 1.0;
613  if ( is_CC && is_delta ) {
614  if ( (is_nu && is_p) || (is_nubar && is_n) ) mult=3.0;
615  }
616  xsec *= mult;
617 
618  // Check whether the cross section is to be weighted with a Breit-Wigner distribution
619  // (default: true)
620  double bw = 1.0;
621  if ( fWghtBW ) {
622  bw = utils::bwfunc::BreitWignerL(W,LR,MR,WR,NR);
623  }
624 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
625  LOG("BSKLNBaseRESPXSec2014", pDEBUG)
626  << "BreitWigner(RES=" << resname << ", W=" << W << ") = " << bw;
627 #endif
628  xsec *= bw;
629 
630 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
631  LOG("BSKLNBaseRESPXSec2014", pINFO)
632  << "\n d2xsec/dQ2dW" << "[" << interaction->AsString()
633  << "](W=" << W << ", q2=" << q2 << ", E=" << E << ") = " << xsec;
634 #endif
635 
636  // The algorithm computes d^2xsec/dWdQ2
637  // Check whether variable tranformation is needed
638  if ( kps != kPSWQ2fE ) {
640  xsec *= J;
641  }
642 
643  // Apply given scaling factor
644  if (is_CC) { xsec *= fXSecScaleCC; }
645  else if (is_NC) { xsec *= fXSecScaleNC; }
646 
647  // If requested return the free nucleon xsec even for input nuclear tgt
648  if ( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;
649 
650  int Z = target.Z();
651  int A = target.A();
652  int N = A-Z;
653 
654  // Take into account the number of scattering centers in the target
655  int NNucl = (is_p) ? Z : N;
656  xsec*=NNucl; // nuclear xsec (no nuclear suppression factor)
657 
658  if ( fUsePauliBlocking && A!=1 )
659  {
660  // Calculation of Pauli blocking according references:
661  //
662  // [1] S.L. Adler, S. Nussinov, and E.A. Paschos, "Nuclear
663  // charge exchange corrections to leptonic pion production
664  // in the (3,3) resonance region," Phys. Rev. D 9 (1974)
665  // 2125-2143 [Erratum Phys. Rev. D 10 (1974) 1669].
666  // [2] J.Y. Yu, "Neutrino interactions and nuclear effects in
667  // oscillation experiments and the nonperturbative disper-
668  // sive sector in strong (quasi-)abelian fields," Ph. D.
669  // Thesis, Dortmund U., Dortmund, 2002 (unpublished).
670  // [3] E.A. Paschos, J.Y. Yu, and M. Sakuda, "Neutrino pro-
671  // duction of resonances," Phys. Rev. D 69 (2004) 014013
672  // [arXiv: hep-ph/0308130].
673 
674  double P_Fermi = 0.0;
675 
676  // Maximum value of Fermi momentum of target nucleon (GeV)
677  if ( A<6 || ! fUseRFGParametrization )
678  {
679  // look up the Fermi momentum for this target
681  const FermiMomentumTable * kft = kftp->GetTable(fKFTable);
682  P_Fermi = kft->FindClosestKF(pdg::IonPdgCode(A, Z), nucpdgc);
683  }
684  else {
685  // define the Fermi momentum for this target
687  // correct the Fermi momentum for the struck nucleon
688  if(is_p) { P_Fermi *= TMath::Power( 2.*Z/A, 1./3); }
689  else { P_Fermi *= TMath::Power( 2.*N/A, 1./3); }
690  }
691 
692  double FactorPauli_RES = 1.0;
693 
694  double k0 = 0., q = 0., q0 = 0.;
695 
696  if (P_Fermi > 0.)
697  {
698  k0 = (W2-Mnuc2-Q2)/(2*W);
699  k = TMath::Sqrt(k0*k0+Q2); // previous value of k is overridden
700  q0 = (W2-Mnuc2+kPionMass2)/(2*W);
701  q = TMath::Sqrt(q0*q0-kPionMass2);
702  }
703 
704  if ( 2*P_Fermi < k-q )
705  FactorPauli_RES = 1.0;
706  if ( 2*P_Fermi >= k+q )
707  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);
708  if ( 2*P_Fermi >= k-q && 2*P_Fermi <= k+q )
709  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);
710 
711  xsec *= FactorPauli_RES;
712  }
713  return xsec;
714 }
bool IsDelta(Resonance_t res)
is it a Delta resonance?
bool fNormBW
normalize resonance breit-wigner to 1?
virtual const RSHelicityAmpl & Compute(Resonance_t res, const FKR &fkr) const =0
string fKFTable
table of Fermi momentum (kF) constants for various nuclei
double fOmega
FKR parameter Omega.
double W(bool selected=false) const
Definition: Kinematics.cxx:157
bool IsWeakCC(void) const
static const double kSqrt2
Definition: Constants.h:115
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:107
bool fUsingDisResJoin
use a DIS/RES joining scheme?
double fXSecScaleNC
external NC xsec scaling factor
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
double Rminus
Definition: FKR.h:50
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1064
int HitNucPdg(void) const
Definition: Target.cxx:304
double Ra
Definition: FKR.h:42
double Amp2Plus3(void) const
int A(void) const
Definition: Target.h:70
double Amp2Minus3(void) const
double HitNucMass(void) const
Definition: Target.cxx:233
double fN0ResMaxNWidths
limits allowed phase space for n=0 res
static FermiMomentumTablePool * Instance(void)
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
double Lamda
Definition: FKR.h:37
double Mass(Resonance_t res)
resonance mass (GeV)
double R
Definition: FKR.h:45
A table of Fermi momentum constants.
double Width(Resonance_t res)
resonance width (GeV)
double Amp2Plus1(void) const
double Amp2Minus1(void) const
return |helicity amplitude|^2
double BreitWignerL(double W, int L, double mass, double width0, double norm)
Definition: BWFunc.cxx:99
double BWNorm(Resonance_t res, double N0ResMaxNWidths=6, double N2ResMaxNWidths=2, double GnResMaxNWidths=4)
breit-wigner normalization factor
enum genie::EResonance Resonance_t
const RSHelicityAmplModelI * fHAmplModelEMp
const RSHelicityAmplModelI * fHAmplModelCC
double fVud2
|Vud|^2(square of magnitude ud-element of CKM-matrix)
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:338
bool IsPosChargedLepton(int pdgc)
Definition: PDGUtils.cxx:145
double Tv
Definition: FKR.h:38
virtual bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
double q2(bool selected=false) const
Definition: Kinematics.cxx:141
A class holding the Rein-Sehgal&#39;s helicity amplitudes.
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:333
bool IsWeakNC(void) const
const TLorentzVector & FSLeptonP4(void) const
Definition: Kinematics.h:65
Singleton class to load & serve tables of Fermi momentum constants.
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
bool fWghtBW
weight with resonance breit-wigner?
const FermiMomentumTable * GetTable(string name)
static const double kAem2
Definition: Constants.h:57
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
const double a
double T
Definition: FKR.h:46
double Rv
Definition: FKR.h:39
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:115
double fXSecScaleCC
external CC xsec scaling factor
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
int ProbePdg(void) const
Definition: InitialState.h:64
bool fUsePauliBlocking
account for Pauli blocking?
double fWcut
apply DIS/RES joining scheme < Wcut
int OrbitalAngularMom(Resonance_t res)
orbital angular momentum
const RSHelicityAmplModelI * fHAmplModelEMn
int Z(void) const
Definition: Target.h:68
#define pINFO
Definition: Messenger.h:62
Pure abstract base class. Defines the RSHelicityAmplModelI interface.
double Amp20Minus(void) const
bool IsEM(void) const
double fGnResMaxNWidths
limits allowed phase space for other res
double C
Definition: FKR.h:44
double FermiMomentumForIsoscalarNucleonParametrization(const Target &target)
bool fUseRFGParametrization
use parametrization for fermi momentum insted of table?
const RSHelicityAmplModelI * fHAmplModelNCp
double Tplus
Definition: FKR.h:47
double B
Definition: FKR.h:43
double Rplus
Definition: FKR.h:49
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
E
Definition: 018_def.c:13
double Tminus
Definition: FKR.h:48
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:68
double fSin48w
sin^4(Weingberg angle)
#define A
Definition: memgrp.cpp:38
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:130
double Amp20Plus(void) const
const char * AsString(Resonance_t res)
resonance id -> string
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
double FindClosestKF(int target_pdgc, int nucleon_pdgc) const
const RSHelicityAmplModelI * fHAmplModelNCn
const Target & Tgt(void) const
Definition: InitialState.h:66
static const double kGF2
Definition: Constants.h:59
double fN2ResMaxNWidths
limits allowed phase space for n=2 res
double fZeta
FKR parameter Zeta.
double ProbeE(RefFrame_t rf) const
static const double kPi
Definition: Constants.h:37
bool IsNegChargedLepton(int pdgc)
Definition: PDGUtils.cxx:136
static const double kPi2
Definition: Constants.h:38
double S
Definition: FKR.h:40
double Ta
Definition: FKR.h:41
int ResonanceIndex(Resonance_t res)
resonance idx, quark model / SU(6)
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63
static const double kPionMass2
Definition: Constants.h:86

Member Data Documentation

bool genie::BSKLNBaseRESPXSec2014::fBRS
protected

Definition at line 101 of file BSKLNBaseRESPXSec2014.h.

FKR genie::BSKLNBaseRESPXSec2014::fFKR
mutableprotected

Definition at line 70 of file BSKLNBaseRESPXSec2014.h.

bool genie::BSKLNBaseRESPXSec2014::fGA
protected

Definition at line 103 of file BSKLNBaseRESPXSec2014.h.

double genie::BSKLNBaseRESPXSec2014::fGnResMaxNWidths
protected

limits allowed phase space for other res

Definition at line 92 of file BSKLNBaseRESPXSec2014.h.

bool genie::BSKLNBaseRESPXSec2014::fGV
protected

Definition at line 104 of file BSKLNBaseRESPXSec2014.h.

const RSHelicityAmplModelI* genie::BSKLNBaseRESPXSec2014::fHAmplModelCC
protected

Definition at line 72 of file BSKLNBaseRESPXSec2014.h.

const RSHelicityAmplModelI* genie::BSKLNBaseRESPXSec2014::fHAmplModelEMn
protected

Definition at line 76 of file BSKLNBaseRESPXSec2014.h.

const RSHelicityAmplModelI* genie::BSKLNBaseRESPXSec2014::fHAmplModelEMp
protected

Definition at line 75 of file BSKLNBaseRESPXSec2014.h.

const RSHelicityAmplModelI* genie::BSKLNBaseRESPXSec2014::fHAmplModelNCn
protected

Definition at line 74 of file BSKLNBaseRESPXSec2014.h.

const RSHelicityAmplModelI* genie::BSKLNBaseRESPXSec2014::fHAmplModelNCp
protected

Definition at line 73 of file BSKLNBaseRESPXSec2014.h.

string genie::BSKLNBaseRESPXSec2014::fKFTable
protected

table of Fermi momentum (kF) constants for various nuclei

Definition at line 93 of file BSKLNBaseRESPXSec2014.h.

bool genie::BSKLNBaseRESPXSec2014::fKLN
protected

Definition at line 100 of file BSKLNBaseRESPXSec2014.h.

double genie::BSKLNBaseRESPXSec2014::fMa2
protected

(axial mass)^2

Definition at line 83 of file BSKLNBaseRESPXSec2014.h.

double genie::BSKLNBaseRESPXSec2014::fMv2
protected

(vector mass)^2

Definition at line 84 of file BSKLNBaseRESPXSec2014.h.

double genie::BSKLNBaseRESPXSec2014::fN0ResMaxNWidths
protected

limits allowed phase space for n=0 res

Definition at line 91 of file BSKLNBaseRESPXSec2014.h.

double genie::BSKLNBaseRESPXSec2014::fN2ResMaxNWidths
protected

limits allowed phase space for n=2 res

Definition at line 90 of file BSKLNBaseRESPXSec2014.h.

bool genie::BSKLNBaseRESPXSec2014::fNormBW
protected

normalize resonance breit-wigner to 1?

Definition at line 80 of file BSKLNBaseRESPXSec2014.h.

double genie::BSKLNBaseRESPXSec2014::fOmega
protected

FKR parameter Omega.

Definition at line 82 of file BSKLNBaseRESPXSec2014.h.

double genie::BSKLNBaseRESPXSec2014::fSin48w
protected

sin^4(Weingberg angle)

Definition at line 85 of file BSKLNBaseRESPXSec2014.h.

bool genie::BSKLNBaseRESPXSec2014::fUsePauliBlocking
protected

account for Pauli blocking?

Definition at line 95 of file BSKLNBaseRESPXSec2014.h.

bool genie::BSKLNBaseRESPXSec2014::fUseRFGParametrization
protected

use parametrization for fermi momentum insted of table?

Definition at line 94 of file BSKLNBaseRESPXSec2014.h.

bool genie::BSKLNBaseRESPXSec2014::fUsingDisResJoin
protected

use a DIS/RES joining scheme?

Definition at line 87 of file BSKLNBaseRESPXSec2014.h.

bool genie::BSKLNBaseRESPXSec2014::fUsingNuTauScaling
protected

use NeuGEN nutau xsec reduction factors?

Definition at line 88 of file BSKLNBaseRESPXSec2014.h.

double genie::BSKLNBaseRESPXSec2014::fVud2
protected

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

Definition at line 86 of file BSKLNBaseRESPXSec2014.h.

double genie::BSKLNBaseRESPXSec2014::fWcut
protected

apply DIS/RES joining scheme < Wcut

Definition at line 89 of file BSKLNBaseRESPXSec2014.h.

bool genie::BSKLNBaseRESPXSec2014::fWghtBW
protected

weight with resonance breit-wigner?

Definition at line 79 of file BSKLNBaseRESPXSec2014.h.

const XSecIntegratorI* genie::BSKLNBaseRESPXSec2014::fXSecIntegrator
protected

Definition at line 106 of file BSKLNBaseRESPXSec2014.h.

double genie::BSKLNBaseRESPXSec2014::fXSecScaleCC
protected

external CC xsec scaling factor

Definition at line 97 of file BSKLNBaseRESPXSec2014.h.

double genie::BSKLNBaseRESPXSec2014::fXSecScaleNC
protected

external NC xsec scaling factor

Definition at line 98 of file BSKLNBaseRESPXSec2014.h.

double genie::BSKLNBaseRESPXSec2014::fZeta
protected

FKR parameter Zeta.

Definition at line 81 of file BSKLNBaseRESPXSec2014.h.


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