Public Member Functions | Private Member Functions | Private Attributes | List of all members
genie::alvarezruso::AlvarezRusoCOHPiPDXSec Class Reference

#include <AlvarezRusoCOHPiPDXSec.h>

Public Member Functions

 AlvarezRusoCOHPiPDXSec (unsigned int Z_, unsigned int A_, const current_t current_, const flavour_t flavour_=kE, const nutype_t nutype=kNu, const formfactors_t ff_=kNieves)
 
 ~AlvarezRusoCOHPiPDXSec ()
 
double DXSec (const double E_nu_, const double E_l_, const double theta_l_, const double phi_l_, const double theta_pi_, const double phi_pi_)
 
void SetDebug (bool debug)
 
ARConstantsGetConstants (void)
 
ARSampledNucleusGetNucleus (void)
 
int GetSampling () const
 
double GetPiMass () const
 
double GetLeptonMass () const
 

Private Member Functions

void SetKinematics ()
 
void SetFlavour ()
 
void SetCurrent ()
 
std::complex< double > DeltaCouplingInMed (ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > delta_momentum, ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > pion_momentum, double density_cent)
 
double PiDecayVertex (ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > pion_momentum, double mass)
 
std::complex< double > DeltaPropagatorInMed (ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > delta_momentum)
 
double DeltaWidthPauliBlocked (ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > delta_momentum, double density)
 
double DeltaWidthFree (ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > delta_momentum)
 
std::complex< double > H (unsigned int i, unsigned int j) const
 
double DifferentialCrossSection ()
 
double PionMomentumCM (ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > delta_momentum)
 
double PNVertexFactor (ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > momentum, double mass)
 
double DeltaSelfEnergyRe (double density)
 
double DeltaSelfEnergyIm (double density)
 
double DeltaSelfEnergyConstant (double a, double b, double c, double E)
 
std::complex< double > NucleonPropagator (ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > nucleon_momentum)
 
void NuclearCurrent (ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > q, ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > pdir, ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > pcrs, ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > ppi, std::complex< double > *jPtr)
 
void SolveWavefunctions ()
 

Private Attributes

bool debug_
 
unsigned int fZ
 
unsigned int fA
 
unsigned int fSampling
 
current_t current
 
flavour_t flavour
 
nutype_t nutype
 
formfactors_t formfactors
 
ARConstantsfConstants
 
ARSampledNucleusfNucleus
 
ARWFSolutionfWfsolution
 
double fE_nu
 
double fE_l
 
double fTheta_l
 
double fTheta_pi
 
double fPhi
 
double fLastE_pi
 
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fQ
 
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fP_nu
 
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fP_l
 
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fP_pi
 
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fP_n_i
 
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fP_n_o
 
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fP_direct
 
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fP_cross
 
double fM_pi
 
double fM_l
 
double fg_factor
 
double fF_direct_delta
 
double fF_direct_nucleon
 
double fF_cross_delta
 
double fF_cross_nucleon
 
ARWavefunctionfUwave
 
ARWavefunctionfUwaveDr
 
ARWavefunctionfUwaveDtheta
 
std::complex< double > fJ_hadronic [4]
 

Detailed Description

Definition at line 46 of file AlvarezRusoCOHPiPDXSec.h.

Constructor & Destructor Documentation

genie::alvarezruso::AlvarezRusoCOHPiPDXSec::AlvarezRusoCOHPiPDXSec ( unsigned int  Z_,
unsigned int  A_,
const current_t  current_,
const flavour_t  flavour_ = kE,
const nutype_t  nutype = kNu,
const formfactors_t  ff_ = kNieves 
)

Definition at line 49 of file AlvarezRusoCOHPiPDXSec.cxx.

51  : debug_(false),
52  fZ(Z_),
53  fA(A_),
54  //~ fSampling( A_ >= 56 ? 48 : 20),
55  fSampling(20),
56  current( current_ ),
57  flavour( flavour_ ),
58  nutype( nutype_ ),
59  formfactors( ff_ ),
60  fConstants ( new ARConstants() ),
61  fNucleus ( new ARSampledNucleus(fZ, fA, fSampling) ),
62  fWfsolution ( new AREikonalSolution(debug_, this) ),
63  fLastE_pi (-9999999.),
64  fUwave ( new ARWavefunction(fSampling, debug_) ),
65  fUwaveDr ( new ARWavefunction(fSampling, debug_) ),
66  fUwaveDtheta( new ARWavefunction(fSampling, debug_) )
67 {
68  SetCurrent();
69  SetFlavour();
70 
71 }
genie::alvarezruso::AlvarezRusoCOHPiPDXSec::~AlvarezRusoCOHPiPDXSec ( )

Definition at line 73 of file AlvarezRusoCOHPiPDXSec.cxx.

74 {
75  delete this->fWfsolution;
76  delete this->fUwave;
77  delete this->fUwaveDr;
78  delete this->fUwaveDtheta;
79  delete this->fNucleus;
80  delete this->fConstants;
81 }

Member Function Documentation

cdouble genie::alvarezruso::AlvarezRusoCOHPiPDXSec::DeltaCouplingInMed ( ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > >  delta_momentum,
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > >  pion_momentum,
double  density_cent 
)
private

Definition at line 973 of file AlvarezRusoCOHPiPDXSec.cxx.

974 {
975  cdouble gdmed;
976  cdouble s_delta (delta_momentum.mag2(),0);
977  cdouble I(0,1);
978  if( real(s_delta) < ((fConstants->NucleonMass() + fM_pi)*(fConstants->NucleonMass() + fM_pi)) )
979  {
980  gdmed = 1.0 / ( s_delta - fConstants->DeltaPMass()*fConstants->DeltaPMass() ) ;
981  }
982  else if( delta_momentum.E() < 0.0 )
983  {
984  gdmed = 0.0;
985  }
986  else
987  {
988  cdouble sqrt_delta = sqrt(s_delta);
989  double gamdpb = DeltaWidthPauliBlocked(delta_momentum, density);
990  double real = DeltaSelfEnergyRe(density);
991  double imaginary = DeltaSelfEnergyIm(density);
992  double ofshel = PiDecayVertex(pion_momentum, fConstants->DeltaPMass());
993 
994  cdouble part_1 = sqrt_delta - fConstants->DeltaPMass() + (ofshel*ofshel*(I*gamdpb)/2.0) - real - I*imaginary;
995  cdouble part_2 = sqrt_delta + fConstants->DeltaPMass();
996 
997  gdmed = 1.0 / (part_1 * part_2);
998  }
999 
1000  return gdmed;
1001 }
double DeltaWidthPauliBlocked(ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > delta_momentum, double density)
double PiDecayVertex(ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > pion_momentum, double mass)
std::complex< double > cdouble
cdouble genie::alvarezruso::AlvarezRusoCOHPiPDXSec::DeltaPropagatorInMed ( ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > >  delta_momentum)
private

Definition at line 358 of file AlvarezRusoCOHPiPDXSec.cxx.

359 {
360  //Energy dependent in-medium Delta propagator
361  double W = TMath::Abs( delta_momentum.mag() );
362  double width = DeltaWidthPauliBlocked(delta_momentum, 0.0);
363  double imSigma = DeltaSelfEnergyIm(0.0);
364  double reSigma = DeltaSelfEnergyRe(0.0);
365 
366  cdouble denom1( (W + fConstants->DeltaPMass() + reSigma), 0.0);
367  cdouble denom2( (W - fConstants->DeltaPMass() - reSigma), ((width/2.0) - imSigma) );
368 
369  cdouble result = 1.0 / (denom1 * denom2);
370  return result;
371 }
static QCString result
double DeltaWidthPauliBlocked(ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > delta_momentum, double density)
std::complex< double > cdouble
double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::DeltaSelfEnergyConstant ( double  a,
double  b,
double  c,
double  E 
)
private

Definition at line 547 of file AlvarezRusoCOHPiPDXSec.cxx.

548 {
549  double x = (E / fM_pi) - 1.0;
550  return (a*x*x + b*x + c);
551 }
const double a
static bool * b
Definition: config.cpp:1043
list x
Definition: train.py:276
double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::DeltaSelfEnergyIm ( double  density)
private

Definition at line 502 of file AlvarezRusoCOHPiPDXSec.cxx.

503 {
504  // Oset, Salcedo, NPA 468(87)631
505  // Using eq. (3.5) to relate the energy of the delta with the pion
506  // energy used in the parametrization
507 
508  double E = fP_pi.E();
509 
510  // The parameterization is valid for 85 MeV < tpi < 315
511  // above which we take a contant values
512 
513  if( E >= (fM_pi + (0.315/fConstants->HBar())) )
514  {
515  E = fM_pi + 0.315/fConstants->HBar();
516  }
517 
518  double Cq = DeltaSelfEnergyConstant(-5.19, 15.35, 2.06, E) / (1000.0 * fConstants->HBar());
519  double Ca2 = DeltaSelfEnergyConstant(1.06, -6.64, 22.66, E) / (1000.0 * fConstants->HBar());
520 
521  // Ca3 extrapolated to zero at low kin. energies
522  double Ca3;
523  if( E <= (fM_pi + (0.085/fConstants->HBar())) )
524  {
525  Ca3 = DeltaSelfEnergyConstant(-13.46, 46.17 , -20.34, (fM_pi + 0.085/(1000.0*fConstants->HBar())));
526  Ca3 /= 85.0;
527  Ca3 *= (E - fM_pi);
528  }
529  else
530  {
531  Ca3 = DeltaSelfEnergyConstant(-13.46, 46.17, -20.34, E) / (1000.0 * fConstants->HBar());
532  }
533  double alpha = DeltaSelfEnergyConstant(0.382, -1.322, 1.466, E);
534  double beta = DeltaSelfEnergyConstant(-0.038, 0.204, 0.613, E);
535  double gamma = 2.0*beta;
536 
537  double ratio = density / fConstants->Rho0();
538 
539  double result = Cq * TMath::Power(ratio, alpha);
540  result += Ca2 * TMath::Power(ratio, beta);
541  result += Ca3 * TMath::Power(ratio, gamma);
542  result *= -1.0;
543 
544  return result;
545 }
double beta(double KE, const simb::MCParticle *part)
static QCString result
double alpha
Definition: doAna.cpp:15
double gamma(double KE, const simb::MCParticle *part)
double DeltaSelfEnergyConstant(double a, double b, double c, double E)
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fP_pi
double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::DeltaSelfEnergyRe ( double  density)
private

Definition at line 496 of file AlvarezRusoCOHPiPDXSec.cxx.

497 {
498  double result = (0.04/fConstants->HBar())*(density/fConstants->Rho0());
499  return result;
500 }
static QCString result
double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::DeltaWidthFree ( ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > >  delta_momentum)
private

Definition at line 431 of file AlvarezRusoCOHPiPDXSec.cxx.

432 {
433  // Free Delta width
434  double pre_factor_1 = 1.0 / (6.0 * kPi );
435  double pre_factor_2 = fConstants->DeltaNCoupling()*fConstants->DeltaNCoupling()/(fM_pi*fM_pi);
436 
437  double qcm = PionMomentumCM(delta_momentum);
438  double qcm3 = qcm*qcm*qcm;
439  double mn = fConstants->NucleonMass();
440  // Luis' code has the next pre-factor equal to
441  // double pre_factor_3 = fConstants->NucleonMass() / TMath::Sqrt(TMath::Abs( delta_momentum.mag2() ));
442  // but the paper uses the Delta mass?
443  double p = sqrt(delta_momentum.mag2());
444 
445  if (not(p>=fM_pi+mn)) {
446  LOG("AlvarezRusoCOHPiPDXSec",pWARN)<<"DeltaWidthFree >> Delta invariant mass less than pion mass plus nucleon mass: "<<p;
447  }
448 
449  double pre_factor_3 = mn/p;
450 
451  double delta_width = pre_factor_1 * pre_factor_2 * pre_factor_3 * qcm3;
452 
453  return delta_width;
454 }
double PionMomentumCM(ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > delta_momentum)
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
p
Definition: test.py:223
#define pWARN
Definition: Messenger.h:60
static const double kPi
Definition: Constants.h:37
double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::DeltaWidthPauliBlocked ( ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > >  delta_momentum,
double  density 
)
private

Definition at line 373 of file AlvarezRusoCOHPiPDXSec.cxx.

374 {
375  //In-medium Delta width including Pauli blocking
376 
377  double width;
378  double free_width = DeltaWidthFree(delta_momentum);
379 
380  if(free_width == 0.0)
381  {
382  width = 0.0;
383  }
384  else
385  {
386  double f = 0.0;
387  double p_f = TMath::Power( ((3.0/2.0)*constants::kPi2*density), (1.0/3.0) ); // Fermi-momentum
388  double p_cm = PionMomentumCM(delta_momentum); // nucleon (and pion) momentum in CoM
389 
390  if(formfactors == kNieves)
391  {
392  // Use the approximation from Nieves et al. NPA 554(93)554
393  double r = p_cm / p_f;
394  if(r > 1.0)
395  {
396  double r2 = r*r;
397  f = 1.0 + ( (-2.0)/(5.0*r2) + (9.0)/(35.0*r2*r2) - (2.0)/(21.0*r2*r2*r2) );
398  }
399  else
400  {
401  f = (34.0/35.0)*r - (22.0/105.0)*r*r*r;
402  }
403  }
404  else if(formfactors == kGarcia)
405  {
406  //Use the approximation from Garcia-Recio, NPA 526(91)685
407  // Delta inv. mass
408  double wd = TMath::Abs( delta_momentum.M());
409  // Fermi-energy
410  double E_f = TMath::Sqrt( fConstants->NucleonMassSq() + p_f*p_f );
411  // Delta 3-momentum in Lab.
412  double kd = delta_momentum.R();
413  // Nucleon energy in CoM
414  double E_n = TMath::Sqrt( fConstants->NucleonMassSq() + p_cm*p_cm );
415  f = (kd*p_cm + delta_momentum.E()*E_n - E_f*wd) / (2.0*kd*p_cm);
416 
417  if(f < 0.0) f = 0;
418  else if(f > 1.0) f = 1.0;
419  }
420  else
421  {
422  LOG("AlvarezRusoCOHPiPDXSec",pERROR) << "[ERROR] Choice of form-factor approximation not properly made";
423  exit(1);
424  }
425 
426  width = free_width*f;
427  }
428  return width;
429 }
#define pERROR
Definition: Messenger.h:59
double DeltaWidthFree(ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > delta_momentum)
double PionMomentumCM(ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > delta_momentum)
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static const double kPi2
Definition: Constants.h:38
double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::DifferentialCrossSection ( )
private

Definition at line 151 of file AlvarezRusoCOHPiPDXSec.cxx.

152 {
153  const cdouble i(0,1);
154 
155  cdouble term1,term2,term3,term4;
156  if (nutype == kNu) {
157  term1 = fQ.X() * ( H(0,1) - i*H(0,2) + H(1,0) - H(1,3) + i*H(2,0) - i*H(2,3) -H(3,1) + i*H(3,2) );
158  term2 = fQ.Z() * ( -H(0,0) + H(0,3) + H(1,1) - i*H(1,2) + i*H(2,1) + H(2,2) + H(3,0) - H(3,3) );
159  term3 = 2.0 * fP_nu.E() * ( H(0,0) - H(0,3) - H(3,0) + H(3,3) );
160  term4 = -fQ.E() * ( H(0,0) - H(0,3) + H(1,1) - i*H(1,2) + i*H(2,1) + H(2,2) - H(3,0) + H(3,3) );
161  } else {
162  term1 = fQ.X() * ( H(0,1) + i*H(0,2) + H(1,0) - H(1,3) - i*H(2,0) + i*H(2,3) -H(3,1) - i*H(3,2) );
163  term2 = fQ.Z() * ( -H(0,0) + H(0,3) + H(1,1) + i*H(1,2) - i*H(2,1) + H(2,2) + H(3,0) - H(3,3) );
164  term3 = 2.0 * fP_nu.E() * ( H(0,0) - H(0,3) - H(3,0) + H(3,3) );
165  term4 = -fQ.E() * ( H(0,0) - H(0,3) + H(1,1) + i*H(1,2) - i*H(2,1) + H(2,2) - H(3,0) + H(3,3) );
166  }
167 
168  cdouble complex_amp2 = 8.0 * fP_nu.E() * (term1 + term2 + term3 + term4);
169 
170  double amp2 = real(complex_amp2);
171 
172  double d5 = fg_factor / 8.0 * fP_l.P() * fP_pi.P() / fP_nu.E() / (32 * kPi4 * kPi ) *
173  amp2 * fConstants->cm38Conversion() / fConstants->HBar();
174 
175  return d5;
176 }
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fP_nu
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fQ
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fP_l
std::complex< double > H(unsigned int i, unsigned int j) const
std::complex< double > cdouble
static const double kPi4
Definition: Constants.h:40
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fP_pi
static const double kPi
Definition: Constants.h:37
double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::DXSec ( const double  E_nu_,
const double  E_l_,
const double  theta_l_,
const double  phi_l_,
const double  theta_pi_,
const double  phi_pi_ 
)

Definition at line 84 of file AlvarezRusoCOHPiPDXSec.cxx.

85 {
86 
87  fE_nu = E_nu_;
88  fE_l = E_l_;
89  fTheta_l = theta_l_;
90  fTheta_pi = theta_pi_;
91  fPhi = phi_pi_ - phi_l_;
92 
93  if( (fE_nu/fConstants->HBar()) < (fM_pi + fM_l) )
94  {
95  return 0.0;
96  }
97  else if( (fE_l /fConstants->HBar()) < fM_l )
98  {
99  return 0.0;
100  }
101 
102  SetKinematics();
103 
104  if( fP_pi.E() < fM_pi )
105  {
106  LOG("AlvarezRusoCOHPiPDXSec",pERROR)<<"Pion energy smaller than mass: "<<fP_pi.E();
107  return 0.0;
108  }
109 
110  if( TMath::Abs((fQ - fP_pi).M()) > (2.0 / fConstants->HBar()) )
111  {
112  // Comment from original fortran:
113  // This is to eliminate very high momentum transfers to the nucleus, which
114  // have negligible impact on observables but might create numerical instabilities
115  return 0.0;
116  }
117 
122 
123  // Only need to resolve wave funtions if Epi changes
124  if ( TMath::Abs(fLastE_pi-fP_pi.E()) > 1E-10 ){
126  }
127 
128  LorentzVector pni = fP_pi - fQ;
129  pni *= 0.5;
130  pni.SetE( fConstants->NucleonMass() );
131  fP_direct = fQ + pni;
132  fP_cross = pni - fP_pi;
134 
135  double dxsec = DifferentialCrossSection();
136 
137  fLastE_pi = fP_pi.E();
138 
139  return dxsec;
140 }
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fQ
#define pERROR
Definition: Messenger.h:59
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > LorentzVector
void NuclearCurrent(ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > q, ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > pdir, ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > pcrs, ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > ppi, std::complex< double > *jPtr)
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
double PiDecayVertex(ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > pion_momentum, double mass)
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fP_cross
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fP_direct
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fP_pi
ARConstants & genie::alvarezruso::AlvarezRusoCOHPiPDXSec::GetConstants ( void  )

Definition at line 1014 of file AlvarezRusoCOHPiPDXSec.cxx.

1015 {
1016  return *fConstants;
1017 }
double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::GetLeptonMass ( ) const
inline

Definition at line 71 of file AlvarezRusoCOHPiPDXSec.h.

ARSampledNucleus & genie::alvarezruso::AlvarezRusoCOHPiPDXSec::GetNucleus ( void  )

Definition at line 1019 of file AlvarezRusoCOHPiPDXSec.cxx.

1020 {
1021  return *fNucleus;
1022 }
double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::GetPiMass ( ) const
inline

Definition at line 68 of file AlvarezRusoCOHPiPDXSec.h.

int genie::alvarezruso::AlvarezRusoCOHPiPDXSec::GetSampling ( ) const
inline

Definition at line 64 of file AlvarezRusoCOHPiPDXSec.h.

64  {
65  return fSampling;
66  }
cdouble genie::alvarezruso::AlvarezRusoCOHPiPDXSec::H ( unsigned int  i,
unsigned int  j 
) const
private

Definition at line 144 of file AlvarezRusoCOHPiPDXSec.cxx.

145 {
146  cdouble H_ = ( conj(fJ_hadronic[i]) * fJ_hadronic[j] );
147  return H_;
148 }
std::complex< double > cdouble
void genie::alvarezruso::AlvarezRusoCOHPiPDXSec::NuclearCurrent ( ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > >  q,
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > >  pdir,
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > >  pcrs,
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > >  ppi,
std::complex< double > *  jPtr 
)
private

Definition at line 553 of file AlvarezRusoCOHPiPDXSec.cxx.

555 {
556  CVector j1;
557  CVector j2;
558  CVector j3;
559  CVector j4;
560 
561  //double ga = fConstants->GAxial();
562  //double fpi = fConstants->PiDecayConst();
563  double mn = fConstants->NucleonMass();
564  double mn2 = mn*mn;
565  double mn3 = mn*mn*mn;
566  double mdel = fConstants->DeltaPMass();
567  double mdel2 = mdel*mdel;
568  double mdel3 = mdel*mdel*mdel;
569  double mpi = fM_pi;
570  double q0 = q.E();
571  double q02 = q0*q0;
572  double q03 = q0*q0*q0;
573  double q04 = q0*q0*q0*q0;
574  double q05 = q0*q0*q0*q0*q0;
575  double q1 = q.X();
576  double q12 = q1*q1;
577  double q13 = q1*q1*q1;
578  double q14 = q1*q1*q1*q1;
579  double q3 = q.Z();
580  double q32 = q3*q3;
581  double q33 = q3*q3*q3;
582  double q34 = q3*q3*q3*q3;
583  double pi = constants::kPi;
584  double ppim = ppi.P();
585  double ppim2 = ppim*ppim;
586  double ppi1 = ppi.X();
587  double ppi12 = ppi1*ppi1;
588  double ppi2 = ppi.Y();
589  double ppi22 = ppi2*ppi2;
590  double ppi3 = ppi.Z();
591  double ppi32 = ppi3*ppi3;
592  double ppitr = TMath::Sqrt( ppi12 + ppi22 );
593  double ppitr2 = ppitr*ppitr;
594  double fs = fConstants->DeltaNCoupling();
595  double rmax = fNucleus->RadiusMax();
596 
597  cdouble I(0,1);
598  cdouble twoI(0,2);
599 
600  double hbarsq = fConstants->HBar()*fConstants->HBar();
601 
602  double alp;
603  if( current == kCC )
604  {
605  alp = 3.0;
606  }
607  else
608  {
609  alp = 1.0;
610  }
611 
612  double mod;
613  if(current == kCC)
614  {
615  mod = 1.0;
616  }
617  else
618  {
619  mod = constants::kSqrt2;
620  }
621 
622  double t = q.mag2() * hbarsq; // in GeV
623  double Ma2_Delta = fConstants->Ma_Delta() * fConstants->Ma_Delta();
624  double Mv2_Delta = fConstants->Mv_Delta() * fConstants->Mv_Delta();
625 
626  double C3v = 2.05 / ( (1.0 - (t/Mv2_Delta))*(1.0 - (t/Mv2_Delta)) );
627  double C4v = (-fConstants->NucleonMass() / fConstants->DeltaPMass() ) * C3v;
628  double C5v = 0.0;
629 
630  if( current == kNC )
631  {
632  double nc_factor = fConstants->NCFactor();
633  C3v *= nc_factor;
634  C4v *= nc_factor;
635  C5v *= nc_factor;
636  }
637 
638  double C5a = 1.2 * ( 1.0 - ((fConstants->CA5_A() * t)/(fConstants->CA5_B() - t)) ) /
639  ( ( 1.0 - (t / Ma2_Delta))*( 1.0 - (t / Ma2_Delta)) );
640  double C4a = -C5a / 4.0;
641  double C6a = (C5a * mn*mn) / ((mpi*mpi) - (t / hbarsq));
642 
643  // QE Form Factors
644  double F1;
645  double F2;
646  double FA;
647  double FP;
648  {
649  double mun = -1.913;
650  double mup = 2.793;
651 
652  double MNucleon = fConstants->NucleonMass() * fConstants->HBar();
653  double MPion = fM_pi * fConstants->HBar();
654  double Mv2_Nucleon = fConstants->Mv_Nucleon()*fConstants->Mv_Nucleon();
655  double Ma2_Nucleon = fConstants->Ma_Nucleon()*fConstants->Ma_Nucleon();
656 
657  double Q_s = -t;
658  double Q_s2 = Q_s* Q_s;
659  double Q_s3 = Q_s2*Q_s;
660  double Q_s4 = Q_s3*Q_s;
661  double Q_s5 = Q_s4*Q_s;
662  double Q_s6 = Q_s5*Q_s;
663 
664  double tau = Q_s / (4.0 * MNucleon*MNucleon);
665 
666  // parametrization by Budd, Bodek, Arrington (hep-ex/0308005) - BBA2003 formfactors
667  // valid up to t = 6 GeV**2
668 
669  double GEp = 1.0 / (1.0 + 3.253*Q_s + 1.422*Q_s2 + 0.08582*Q_s3 + 0.3318*Q_s4 - 0.09371*Q_s5 + 0.01076*Q_s6);
670  double GMp = mup / (1.0 + 3.104*Q_s + 1.428*Q_s2 + 0.1112*Q_s3 - 0.006981*Q_s4 + 0.0003705*Q_s5 - 7.063e-6*Q_s6);
671  double GEn = ((-mun * 0.942 * tau) / (1.0 + 4.61*tau)) / ( (1.0 + Q_s/Mv2_Nucleon) * (1.0 + Q_s/Mv2_Nucleon) );
672 
673  // parametrization of Krutov (hep-ph/0202183)
674 
675  double GMn = mun / (1.+3.043*Q_s + 0.8548*Q_s2 + 0.6806*Q_s3 - 0.1287*Q_s4 + 0.008912*Q_s5);
676  F1 = ((GEp - GEn) + tau*(GMp - GMn)) / (1.0 + tau);
677  F2 = ((GMp - GMn) - (GEp - GEn)) / (1.0 + tau);
678 
679  FA = 1.0 + ( Q_s / Ma2_Nucleon );
680  FA *= FA;
681  FA = fConstants->GAxial() / FA;
682 
683  FP = (2.0 * MNucleon*MNucleon);
684  FP /= ( MPion*MPion + Q_s );
685  FP *= FA;
686  FP /= fConstants->NucleonMass();
687 
688  if( current == kNC )
689  {
690  double nc_factor = fConstants->NCFactor();
691  F1 *= nc_factor;
692  F2 *= nc_factor;
693  }
694  }
695 
696  // Get q momentum component perpendicular to pion momentum
697  double qper2 = q.P2();
698  double tot = ppi.P2();
699  double dot = q.Vect().Dot(ppi.Vect());
700  if (tot > 0.) qper2 -=dot*dot/tot;
701  qper2 = TMath::Max(qper2,0.);
702 
703  double qper = TMath::Sqrt(qper2);
704  double qpar = dot / ppim;
705 
706  unsigned int n = this->GetNucleus().GetNDensities();
707 
708  std::vector<cdouble > empty_row(n, cdouble(0.0,0.0));
709  std::vector< std::vector<cdouble > > ordez (4, empty_row);
710  std::vector< std::vector<cdouble > > ordez1(4, empty_row);
711  std::vector< std::vector<cdouble > > ordez2(4, empty_row);
712  std::vector< std::vector<cdouble > > ordez3(4, empty_row);
713  std::vector< std::vector<cdouble > > ordez4(4, empty_row);
714 
715  std::vector< std::vector<cdouble > > ordeb2(4, empty_row);
716  // IMPORTANT !!!
717  // ORDEB HAS ITS INDICES REVERSED W.R.T THE FORTRAN
718  std::vector<cdouble > empty_row_backwards(4, cdouble(0.0,0.0));
719  std::vector< std::vector<cdouble > > ordeb(n, empty_row_backwards);
720 
721  cdouble ppi1d, ppi2d, ppi3d;
722 
723  std::vector<cdouble > jnuclear(4);
724 
725 
726  for(unsigned int i = 0; i != n; ++i)
727  {
728  double be = fNucleus->SamplePoint1(i);
729  double bej0 = TMath::BesselJ0( qper*be );
730  double bej1 = TMath::BesselJ1( qper*be );
731 
732  for(unsigned int l = 0; l != n; ++l)
733  {
734  double za = fNucleus->SamplePoint2(l);
735  double r = TMath::Sqrt(za*za + be*be);
736  double r2 = r*r;
737 
738  //double dens = fNucleus->Density(i,l);
739  double dens_cent = fNucleus->DensityOfCentres(i,l);
740  double dens_p_cent = dens_cent * fZ / fA ;
741  double dens_n_cent = dens_cent * (fA-fZ)/fA;
742 
743  cdouble exp_i_qpar_za = exp(I*qpar*za);
744 
745  const cdouble & uwavefunc = (*fUwave)[i][l];
746  const cdouble & uwavedr = (*fUwaveDr)[i][l];
747  const cdouble & uwavedtheta = (*fUwaveDtheta)[i][l];
748 
749  cdouble A = I * ( uwavedr - (za/r2) * uwavedtheta ) * (be/r);
750  cdouble B = I * ( uwavedr*za + (be*be/r2)*uwavedtheta ) / r;
751 
752  // Calculate distorted pion momentum components
753  if( (qper == 0.0) && ppitr != 0.0)
754  {
755  ppi1d = ( q1 * ((ppi12*ppi32)+(ppi22*ppim2)) - q3*ppi1*ppi3*ppitr2 ) /
756  (ppim2*ppitr2)*A*(I*be/2.0) + ppi1/ppim *B*bej0;
757  ppi2d = -ppi2*(ppi1*q1+ppi3*q3)/ppim2*A*(I*be/2.)+ppi2/ppim*B*bej0;
758  ppi3d = -(q1*ppi1*ppi3-q3*ppitr2)/ppim2*A*(I*be/2.)+ppi3/ppim*B*bej0 ;
759  }
760  else if( (qper != 0.0) && (ppitr == 0.0) )
761  {
762  ppi1d = (q1/qper)*A*(I*bej1);
763  ppi2d = 0.0;
764  ppi3d = B*bej0;
765  }
766  else if( (qper == 0.0) && (ppitr == 0.0) )
767  {
768  ppi1d = q1*A*(I*be/2.0);
769  ppi2d = 0.0;
770  ppi3d = B*bej0;
771  }
772  else
773  {
774  ppi1d=(q1*((ppi12*ppi32)+(ppi22*ppim2))-q3*ppi1*ppi3*ppitr2)/(ppim2*ppitr2)/
775  qper*A*(I*bej1)+ppi1/ppim*B*bej0;
776  ppi2d = -ppi2 * (ppi1*q1 + ppi3*q3) / ppim2/qper*A*(I*bej1)+ppi2/ppim*B*bej0;
777  ppi3d=-(q1*ppi1*ppi3-q3*ppitr2)/ppim2/qper*A*(I*bej1)+ppi3/ppim*B*bej0;
778  }
779 
780  // Calculate the current for four different processes (See Fig 1 of Alvarez-Ruso et al,
781  // "Neutral current coherent pion production", arXiv:0707.2172
782  // j1 : current for direct delta production
783  // j2 : current for Crossed delta production
784  // j3 : current for direct nucleon production
785  // j4 : current for crossed nucleon production
786  j1[0] = -4.*(mdel + mn + q0)*(
787  (C5a*mdel2*mn2*q0 + C6a*mdel2*q03 + C5a*mn2*(-mn - q0)*q0*(mn + q0) -
788  C4a*mdel2*q0*q12 - C4a*mdel2*q0*q32 - C6a*q02*(mn + q0)*(q0*(mn + q0) - q12 - q32))*bej0*uwavefunc +
789  ppi1d*(-(C5a*mn2*(-mn - q0)*q1) - C6a*mdel2*q0*q1 + C4a*mdel2*(mn + q0)*q1 +
790  C6a*q0*q1*(q0*(mn + q0) - q12 - q32)) +
791  ppi3d*(-(C5a*mn2*(-mn - q0)*q3) -
792  C6a*mdel2*q0*q3 + C4a*mdel2*(mn + q0)*q3 + C6a*q0*q3*(q0*(mn + q0) - q12 - q32))
793  );
794 
795  j1[1] = (-4.*C6a*mdel3*q02*q1 - 4.*C6a*mdel2*mn*q02*q1 + 4.*C6a*mdel*mn2*q02*q1 + 4.*C6a*mn3*q02*q1 -
796  4.*C6a*mdel2*q03*q1 + 8.*C6a*mdel*mn*q03*q1 + 12.*C6a*mn2*q03*q1 + 4.*C6a*mdel*q04*q1 +
797  12.*C6a*mn*q04*q1 + 4.*C6a*q05*q1 + 4.*C4a*mdel2*q02*(mdel + mn + q0)*q1 +
798  4.*C5a*mn2*q0*(mn + q0)*(mdel + mn + q0)*q1 - 4.*C6a*mdel*mn*q0*q13 - 4.*C6a*mn2*q0*q13 -
799  4.*C6a*mdel*q02*q13 - 8.*C6a*mn*q02*q13 - 4.*C6a*q03*q13 -
800  4.*C6a*q0*(mn + q0)*(mdel + mn + q0)*q1*q32)*bej0*uwavefunc +
801  ppi1d*(-4.*C4a*mdel2*q0*(mn + q0)*(mdel + mn + q0) + 4.*C6a*mdel3*q12 + 4.*C6a*mdel2*mn*q12 +
802  4.*C6a*mdel2*q0*q12 - 4.*C6a*mdel*mn*q0*q12 - 4.*C6a*mn2*q0*q12 - 4.*C6a*mdel*q02*q12 -
803  8.*C6a*mn*q02*q12 - 4.*C6a*q03*q12 + 4.*C6a*mdel*q14 + 4.*C6a*mn*q14 + 4.*C6a*q0*q14 -
804  4.*C5a*mn2*(mdel + mn + q0)*(mdel2 + q12) + 4.*C4a*mdel2*(mdel + mn + q0)*q32 +
805  4.*C6a*(mdel + mn + q0)*q12*q32) +
806  ppi2d*(twoI*C4v*mdel2*(q0*(mn + q0) - q12)*q3 +
807  twoI*C3v*mdel*mn*(2*mdel2 + 2.*mdel*mn - q0*(mn + q0) + q12)*q3 -
808  twoI*mdel*(C4v*mdel - C3v*mn)*q33) +
809  ppi3d*(-4.*C4a*mdel2*(mdel + mn + q0)*q1*q3 -
810  4.*C5a*mn2*(mdel + mn + q0)*q1*q3 + 4.*C6a*(mdel + mn + q0)*q1*(mdel2 - q0*(mn + q0) + q12)*q3 +
811  4.*C6a*(mdel + mn + q0)*q1*q33) +
812  ppi2d*twoI*C5v*mdel2*mn*q0*q3;
813 
814  j1[2] = -2.*I*(-2.*I*C5a*mdel*mn2*ppi2d*(mdel + mn + q0) -
815  twoI*C4a*mdel*ppi2d*(mdel + mn + q0)*(q0*(mn + q0) - q12 - q32) -
816  (ppi3d*q1 - ppi1d*q3)*(C4v*mdel*(q0*(mn + q0) - q12 - q32) +
817  C3v*mn*(2*mdel2 + 2*mdel*mn - q0*(mn + q0) + q12 + q32)))*mdel +
818  twoI*C5v*mdel*mn*q0*(ppi3d*q1 - ppi1d*q3)*mdel;
819 
820  j1[3] = (4.*C4a*mdel2*q02*(mdel + mn + q0)*q3 - 4.*C6a*mdel2*q02*(mdel + mn + q0)*q3 +
821  4.*C5a*mn2*q0*(mn + q0)*(mdel + mn + q0)*q3 + 4.*C6a*q0*(mn + q0)*(mdel + mn + q0)*
822  (q0*(mn + q0) - q12)*q3 - 4.*C6a*q0*(mn + q0)*(mdel + mn + q0)*q33)*bej0*uwavefunc +
823  ppi2d*(-twoI*mdel*q1*(C4v*mdel*(q0*(mn + q0) - q12) +
824  C3v*mn*(2.*mdel2 + 2*mdel*mn - q0*(mn + q0) + q12)) + twoI*mdel*
825  (C4v*mdel - C3v*mn)*q1*q32) +
826  ppi1d*(-4.*C4a*mdel2*(mdel + mn + q0)*q1*q3 +
827  4.*C6a*mdel2*(mdel + mn + q0)*q1*q3 - 4.*C5a*mn2*(mdel + mn + q0)*q1*q3 -
828  4.*C6a*(mdel + mn + q0)*q1*(q0*(mn + q0) - q12)*q3 + 4.*C6a*(mdel + mn + q0)*q1*q33) +
829  ppi3d*(-4.*C4a*mdel2*mn*q0*(mdel + mn + q0) - 4.*C4a*mdel2*(mdel + mn + q0)*(q0 - q1)*(q0 + q1) +
830  4.*C6a*(mdel + mn + q0)*(mdel2 - q0*(mn + q0) + q12)*q32 + 4.*C6a*(mdel + mn + q0)*q34 -
831  4.*C5a*mn2*(mdel + mn + q0)*(mdel2 + q32)) +
832  -twoI*C5v*mdel2*mn*ppi2d*q0*q1;
833 
834  // Crossed Delta
835 
836  j2[0]=-4.*(mdel + mn - q0)*(
837  (C5a*mdel2*mn2*q0 - C5a*mn2*(mn - q0)*(mn - q0)*q0 + C6a*mdel2*q03 -
838  C4a*mdel2*q0*q12 - C4a*mdel2*q0*q32 - C6a*(mn - q0)*q02*((mn - q0)*q0 + q12 + q32))*bej0*uwavefunc +
839  ppi1d*(-(C4a*mdel2*mn*q1) - C5a*mn2*(mn - q0)*q1 + C4a*mdel2*q0*q1 -
840  C6a*mdel2*q0*q1 - C6a*q0*q1*((mn - q0)*q0 + q12 + q32)) +
841  ppi3d*(-(C4a*mdel2*mn*q3) - C5a*mn2*(mn - q0)*q3 + C4a*mdel2*q0*q3 -
842  C6a*mdel2*q0*q3 - C6a*q0*q3*((mn - q0)*q0 + q12 + q32)));
843 
844  j2[1]=(-4.*C5a*mn2*(mn - q0)*(mdel + mn - q0)*q0*q1 - 4.*C6a*mdel3*q02*q1 -
845  4.*C6a*mdel2*mn*q02*q1 + 4.*C6a*mdel*mn2*q02*q1 + 4.*C6a*mn3*q02*q1 +
846  4.*C4a*mdel2*(mdel + mn - q0)*q02*q1 + 4.*C6a*mdel2*q03*q1 -
847  8.*C6a*mdel*mn*q03*q1 - 12.*C6a*mn2*q03*q1 + 4.*C6a*mdel*q04*q1 +
848  12.*C6a*mn*q04*q1 - 4.*C6a*q05*q1 + 4.*C6a*mdel*mn*q0*q13 +
849  4.*C6a*mn2*q0*q13 - 4.*C6a*mdel*q02*q13 - 8.*C6a*mn*q02*q13 +
850  4.*C6a*q03*q13 + 4.*C6a*(mn - q0)*(mdel + mn - q0)*q0*q1*q32)*bej0*uwavefunc +
851  ppi1d*(-4.*C5a*mdel2*mn2*(mdel + mn - q0) + 4.*C4a*mdel2*mn*(mdel + mn - q0)*q0 -
852  4.*C4a*mdel2*(mdel + mn - q0)*q02 + 4.*C6a*mdel3*q12 +
853  4.*C6a*mdel2*mn*q12 - 4.*C5a*mn2*(mdel + mn - q0)*q12 -
854  4.*C6a*mdel2*q0*q12 + 4.*C6a*mdel*mn*q0*q12 + 4.*C6a*mn2*q0*q12 -
855  4.*C6a*mdel*q02*q12 - 8.*C6a*mn*q02*q12 + 4.*C6a*q03*q12 +
856  4.*C6a*mdel*q14 + 4.*C6a*mn*q14 - 4.*C6a*q0*q14 +
857  4.*C4a*mdel2*(mdel + mn - q0)*q32 + 4.*C6a*(mdel + mn - q0)*q12*q32) +
858  ppi2d*(twoI*mdel*(mdel*q0*(-((C4v + C5v)*mn) + C4v*q0) +
859  C3v*mn*(2*mdel*(mdel + mn) + mn*q0 - q02))*q3 +
860  twoI*mdel*(-(C4v*mdel) + C3v*mn)*q12*q3 - twoI*mdel*(C4v*mdel - C3v*mn)*q33) +
861  ppi3d*(-4.*C4a*mdel2*(mdel + mn - q0)*q1*q3 -
862  4.*C5a*mn2*(mdel + mn - q0)*q1*q3 + 4.*C6a*(mdel + mn - q0)*(mdel2 + (mn - q0)*q0)*q1*q3 +
863  4.*C6a*(mdel + mn - q0)*q13*q3 + 4.*C6a*(mdel + mn - q0)*q1*q33);
864 
865  j2[2]=-(twoI*ppi2d*(-twoI*C5a*mdel*mn2*(mdel + mn - q0) +
866  twoI*C4a*mdel*(mdel + mn - q0)*((mn - q0)*q0 + q12 + q32)) -
867  ppi3d*twoI*q1*(C3v*mn*(2*mdel*(mdel + mn) + mn*q0 - q02 + q12 + q32) -
868  mdel*(C5v*mn*q0 + C4v*((mn - q0)*q0 + q12 + q32))) +
869  ppi1d*twoI*q3*(C3v*mn*(2*mdel*(mdel + mn) + mn*q0 - q02 + q12 + q32) -
870  mdel*(C5v*mn*q0 + C4v*((mn - q0)*q0 + q12 + q32)))
871  )*mdel;
872 
873  j2[3] = (-4.*C5a*mn2*(mn - q0)*(mdel + mn - q0)*q0*q3 +
874  4.*C4a*mdel2*(mdel + mn - q0)*q02*q3 - 4.*C6a*mdel2*(mdel + mn - q0)*q02*q3 +
875  4.*C6a*(mn - q0)*(mdel + mn - q0)*q0*((mn - q0)*q0 + q12)*q3 +
876  4.*C6a*(mn - q0)*(mdel + mn - q0)*q0*q33)*bej0*uwavefunc +
877  ppi2d*(-twoI*mdel*q1*(C3v*mn*(2*mdel*(mdel + mn) + mn*q0 - q02 + q12) -
878  mdel*(q0*((C4v + C5v)*mn - C4v*q0) + C4v*q12)) +
879  twoI*mdel*(C4v*mdel - C3v*mn)*q1*q32) +
880  ppi1d*(-4.*C4a*mdel2*(mdel + mn - q0)*q1*q3 + 4.*C6a*mdel2*(mdel + mn - q0)*q1*q3 -
881  4.*C5a*mn2*(mdel + mn - q0)*q1*q3 +
882  4.*C6a*(mdel + mn - q0)*q1*((mn - q0)*q0 + q12)*q3 +
883  4.*C6a*(mdel + mn - q0)*q1*q33) +
884  ppi3d*(-4.*C5a*mdel2*mn2*(mdel + mn - q0) +
885  4.*C4a*mdel2*(mdel + mn - q0)*((mn - q0)*q0 + q12) -
886  4.*C5a*mn2*(mdel + mn - q0)*q32 +
887  4.*C6a*(mdel + mn - q0)*(mdel2 + (mn - q0)*q0 + q12)*q32 +
888  4.*C6a*(mdel + mn - q0)*q34);
889 
890  //
891  // Direct Nucleon
892 
893  j3[0]=-2.0*(FA - FP*q0)*(q02*bej0*uwavefunc - ppi1d*q1 - ppi3d*q3);
894 
895  j3[1] = 2.0*(-(FA*q0*q1) + FP*q02*q1) * bej0 * uwavefunc +
896  2.0*ppi1d*(2.0*FA*mn + FA*q0 - FP*q12) -
897  twoI*(F1 + F2)*ppi2d*q3 - 2.*FP*ppi3d*q1*q3;
898 
899  j3[2]=twoI*(-I*FA*ppi2d*(2.0*mn + q0) - (F1 + F2)*(ppi3d*q1 - ppi1d*q3));
900 
901  j3[3]=twoI*(F1 + F2)*ppi2d*q1 - 2.0*FP*ppi1d*q1*q3 +
902  2.0*(-(FA*q0*q3) + FP*q02*q3)*bej0*uwavefunc +
903  2.0*ppi3d*(2.0*FA*mn + FA*q0 - FP*q32);
904 
905  //
906  // Crossed Nucleon
907 
908  j4[0] = 2.0*(FA + FP*q0)*(q02 *bej0*uwavefunc - ppi1d*q1 - ppi3d*q3);
909 
910  j4[1] = -2.0*(-(FA*q0*q1) - FP*q02*q1)*bej0*uwavefunc - 2.0*ppi1d*(-2.0*FA*mn + FA*q0 + FP*q12) -
911  twoI*(F1 + F2)*ppi2d*q3 - 2.0*FP*ppi3d*q1*q3;
912 
913  j4[2] = twoI*(-I*FA*ppi2d*(2.0*mn - q0) - (F1 + F2)*(ppi3d*q1 - ppi1d*q3));
914 
915  j4[3] = twoI*(F1 + F2)*ppi2d*q1 - 2.0*FP*ppi1d*q1*q3 + 2.0*(FA*q0*q3 + FP*q02*q3)*bej0*uwavefunc +
916  2.0*ppi3d*(2.0*FA*mn - FA*q0 - FP*q32);
917 
918  cdouble pre_factor_1 = mod * I * (fs/mpi) / constants::kSqrt3 *
919  exp_i_qpar_za *
920  (alp*dens_p_cent+dens_n_cent) *
921  DeltaCouplingInMed(pdir,ppi,dens_cent) *
922  pi/(3.0*mn2*mdel2)*fF_direct_delta;
923 
924  cdouble pre_factor_2 = mod * I * (fs/mpi) / constants::kSqrt3 *exp_i_qpar_za *
925  (dens_p_cent+ alp*dens_n_cent)*pi/(3.*mn2*mdel2)*fF_cross_delta *
926  DeltaCouplingInMed(pcrs,ppi,dens_cent);
927 
928  double PreFacMult = (fConstants->GAxial()/constants::kSqrt2/fConstants->PiDecayConst());
929  cdouble pre_factor_3 = 1./mod*(-I)*PreFacMult*exp_i_qpar_za*
930  dens_n_cent*NucleonPropagator(pdir)*pi*fF_direct_nucleon;
931  cdouble pre_factor_4 = 1./mod*(-I)*PreFacMult*exp_i_qpar_za*
932  dens_p_cent*NucleonPropagator(pcrs)*pi*fF_cross_nucleon;
933 
934  for(int m = 0; m != 4; ++m)
935  {
936  ordez1[m][l] = pre_factor_1 * j1[m];
937  ordez2[m][l] = pre_factor_2 * j2[m];
938  ordez3[m][l] = pre_factor_3 * j3[m];
939  ordez4[m][l] = pre_factor_4 * j4[m];
940 
941  ordez[m][l] = ordez1[m][l] + ordez2[m][l] + ordez3[m][l] + ordez4[m][l];
942  }
943  } //l
944 
945  // IMPORTANT !!!
946  // ORDEB HAS ITS INDICES REVERSED W.R.T THE FORTRAN
947 
948  integrationtools::RGN2D(-rmax, rmax, 2, 0, 3, ordez, fSampling, ordeb[i]);
949 
950  for(unsigned int z = 0; z != 4; ++z)
951  {
952  ordeb[i][z] *= be;
953  }
954  }// fSampling point 1 loop (i)
955 
956  for(unsigned int z = 0; z != n; ++z)
957  {
958  for(unsigned int y = 0; y != 4; ++y)
959  {
960  ordeb2[y][z] = ordeb[z][y];
961  }
962  }
963 
964  integrationtools::RGN2D(0.0, rmax, 2, 0, 3, ordeb2, fSampling, jnuclear);
965 
966  for (int i=0; i<4; i++) {
967  cdouble & jn = jnuclear[i];
968  *(jHadCurrent+i) = cdouble(jn);
969  }
970 }
static const double kSqrt3
Definition: Constants.h:116
#define F2(x, y, z)
Definition: md5.c:176
static const double kSqrt2
Definition: Constants.h:115
ROOT::Math::SVector< cdouble, 4 > CVector
std::complex< double > cdouble
#define F1(x, y, z)
Definition: md5.c:175
void RGN2D(const double a, const double b, unsigned int n, unsigned int l, unsigned int m, std::vector< std::vector< std::complex< double > > > &cf, const unsigned int nsamp, std::vector< std::complex< double > > &cres)
static constexpr double fs
Definition: Units.h:100
static QStrList * l
Definition: config.cpp:1044
std::complex< double > DeltaCouplingInMed(ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > delta_momentum, ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > pion_momentum, double density_cent)
double SamplePoint2(const unsigned int i) const
std::void_t< T > n
unsigned int GetNDensities(void) const
std::complex< double > cdouble
#define A
Definition: memgrp.cpp:38
double DensityOfCentres(const int i, const int j) const
float pi
Definition: units.py:11
static const double kPi
Definition: Constants.h:37
double SamplePoint1(const unsigned int i) const
std::complex< double > NucleonPropagator(ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > nucleon_momentum)
cdouble genie::alvarezruso::AlvarezRusoCOHPiPDXSec::NucleonPropagator ( ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > >  nucleon_momentum)
private

Definition at line 1003 of file AlvarezRusoCOHPiPDXSec.cxx.

1004 {
1005  // relativistic nucleon propagator (its denominator)
1006 
1007  cdouble gn( nucleon_momentum.mag2() - fConstants->NucleonMassSq() ,
1008  ( fConstants->NucleonMass()*10.0 / (1000.0*fConstants->HBar()) ) );
1009  gn = 1.0 / gn;
1010 
1011  return gn;
1012 }
std::complex< double > cdouble
double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::PiDecayVertex ( ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > >  pion_momentum,
double  mass 
)
private

Definition at line 183 of file AlvarezRusoCOHPiPDXSec.cxx.

184 {
185  // s-channel form factor for the piNDelta vertex as in M. Post thesis (pg 35)
186  // Calculated for a NUCLEON AT REST
187 
188  double Lam = 1.0 / fConstants->HBar();
189  double Lam4 = Lam*Lam*Lam*Lam;
190 
191  double mass2 = mass*mass;
192 
193  LorentzVector decaying_momentum(momentum.x(),momentum.y(),momentum.z(), (momentum.E()+fConstants->NucleonMass()));
194 
195  double factor = decaying_momentum.mag2() - mass2;
196  factor *= factor;
197 
198  double ofshel = Lam4 / ( Lam4 + factor );
199  return ofshel;
200 }
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > LorentzVector
def momentum(x1, x2, x3, scale=1.)
double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::PionMomentumCM ( ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > >  delta_momentum)
private

Definition at line 461 of file AlvarezRusoCOHPiPDXSec.cxx.

462 {
463  double m_pi2 = fM_pi*fM_pi;
464  double m_n2 = fConstants->NucleonMassSq();
465  double s = delta_momentum.mag2();
466  assert(s>=0);
467 
468  double p_pi_cm = ((s-m_pi2-m_n2)*(s-m_pi2-m_n2)) - 4.0*m_pi2*m_n2;
469 
470  assert(p_pi_cm > 0.);
471 
472  p_pi_cm = TMath::Sqrt(p_pi_cm / s) / 2.0;
473  return p_pi_cm;
474 }
static QCString * s
Definition: config.cpp:1042
double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::PNVertexFactor ( ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > >  momentum,
double  mass 
)
private

Definition at line 476 of file AlvarezRusoCOHPiPDXSec.cxx.

477 {
478  // s-channel form factor for the piNDelta/piNN vertex as in M. Post thesis (pg 35)
479  // Calculated for a NUCLEON AT REST
480 
481  double lambda = 1.0 / fConstants->HBar();
482  double lambda4 = lambda*lambda*lambda*lambda;
483 
484  double mass2 = mass*mass;
485 
486  LorentzVector decaying_momentum(momentum.x(), momentum.y(),momentum.z(), (momentum.E()+fConstants->NucleonMass()));
487 
488  double factor = decaying_momentum.mag2() - mass2;
489  factor *= factor;
490 
491  double result = lambda4 / (lambda4 + factor);
492 
493  return result;
494 }
static QCString result
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > LorentzVector
def momentum(x1, x2, x3, scale=1.)
void genie::alvarezruso::AlvarezRusoCOHPiPDXSec::SetCurrent ( )
private

Definition at line 271 of file AlvarezRusoCOHPiPDXSec.cxx.

272 {
273  switch(this->current)
274  {
275  case kCC:
276  fM_pi = fConstants->PiPMass();
278  break;
279  case kNC:
280  fM_pi = fConstants->Pi0Mass();
282  break;
283  default:
284  LOG("AlvarezRusoCOHPiPDXSec",pERROR) << "[ERROR] Unknown current type";
285  exit(1);
286  }
287 }
#define pERROR
Definition: Messenger.h:59
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
void genie::alvarezruso::AlvarezRusoCOHPiPDXSec::SetDebug ( bool  debug)
inline
void genie::alvarezruso::AlvarezRusoCOHPiPDXSec::SetFlavour ( )
private

Definition at line 233 of file AlvarezRusoCOHPiPDXSec.cxx.

234 {
235  if(current == kNC)
236  {
237  fM_l = 0.0;
238  }
239  else if(current == kCC)
240  {
241  switch(flavour)
242  {
243  case kE:
245  break;
246  case kMu:
247  fM_l = fConstants->MuonMass();
248  break;
249  case kTau:
250  fM_l = fConstants->TauMass();
251  break;
252  default:
253  LOG("AlvarezRusoCOHPiPDXSec",pERROR) << "[ERROR] Unknown lepton flavour";
254  exit(1);
255  }
256  }
257  else
258  {
259  LOG("AlvarezRusoCOHPiPDXSec",pERROR) << "[ERROR] Unknown current";
260  exit(1);
261  }
262 }
#define pERROR
Definition: Messenger.h:59
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
void genie::alvarezruso::AlvarezRusoCOHPiPDXSec::SetKinematics ( )
private

Definition at line 206 of file AlvarezRusoCOHPiPDXSec.cxx.

207 {
208  // Initial neutrino momentum
210 
211  // Final lepton momentum
212  double mod_p_l = TMath::Sqrt( (fE_l/fConstants->HBar())*(fE_l/fConstants->HBar()) - fM_l*fM_l );
213  double p_l_x = mod_p_l * TMath::Sin(fTheta_l);
214  double p_l_z = mod_p_l * TMath::Cos(fTheta_l);
215  fP_l = LorentzVector(p_l_x, 0, p_l_z, (fE_l/fConstants->HBar()));
216 
217  // Momentum transfer
218  fQ = fP_nu - fP_l;
219 
220  // Pion momentum
221  double E_pi = fQ.E();
222  double mod_fP_pi = TMath::Sqrt( E_pi*E_pi - fM_pi*fM_pi );
223  double p_pi_x = mod_fP_pi * TMath::Sin(fTheta_pi) * TMath::Cos(fPhi);
224  double p_pi_y = mod_fP_pi * TMath::Sin(fTheta_pi) * TMath::Sin(fPhi);
225  double p_pi_z = mod_fP_pi * TMath::Cos(fTheta_pi);
226  fP_pi = LorentzVector(p_pi_x, p_pi_y, p_pi_z, E_pi);
227 }
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fP_nu
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fQ
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > LorentzVector
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fP_l
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fP_pi
void genie::alvarezruso::AlvarezRusoCOHPiPDXSec::SolveWavefunctions ( )
private

This is only a function of the nucleus and pion momentum/energy so if neither of those have changed there is no need to re-calculate the wavefunction values. Such a caching has not been implemented here yet!

Definition at line 301 of file AlvarezRusoCOHPiPDXSec.cxx.

302 {
303  unsigned int n_points = fNucleus->GetNDensities();
304 
305  double x2;
306  double radius;
307  double cosine_rz; // angle w.r.t the pion momentum
308 
309  // for calculating derivatives
310  double delta_r;
311  double delta_c;
312  cdouble uwave_plus;
313  cdouble uwave_minus;
314 
315  // Loop over grid of points in the nuclear potential
316  for(unsigned int i = 0; i != n_points; ++i)
317  {
318  for(unsigned int j = 0; j != n_points; ++j)
319  {
320  //double x1 = fNucleus->SamplePoint1(i); // unused
321  x2 = fNucleus->SamplePoint2(j);
322 
323  // radius of position in potential from centre
324  radius = fNucleus->Radius(i,j);
325  // angle of sampling point wrt to neutrino direction
326  cosine_rz = x2 / radius;
327 
328  // Calculate wavefunction
329  fUwave->set(i, j, fWfsolution->Element(radius, -cosine_rz,
330  fP_pi.E()));
331  delta_r = 0.0001;
332  if( radius < delta_r ) delta_r = radius;
333 
334  // Calculate derivative of wavefunction in the radial direction
335  uwave_plus = fWfsolution->Element( (radius+delta_r), -cosine_rz,
336  fP_pi.E());
337  uwave_minus = fWfsolution->Element( (radius-delta_r), -cosine_rz,
338  fP_pi.E());
339 
340  fUwaveDr->set(i, j, (uwave_plus - uwave_minus) / (2.0 * delta_r) );
341 
342  // Calculate derivative of wavefunction in the angle space
343  delta_c = 0.0001;
344  if ( (cosine_rz - delta_c) <= -1.0 ) delta_c = cosine_rz + 1.0 - 1E-12;
345  else if( (cosine_rz + delta_c) >= 1.0 ) delta_c = 1.0 - cosine_rz - 1E-12;
346 
347  uwave_plus = fWfsolution->Element(radius, -(cosine_rz+delta_c),
348  fP_pi.E());
349  uwave_minus = fWfsolution->Element(radius, -(cosine_rz-delta_c),
350  fP_pi.E());
351  fUwaveDtheta->set( i, j, (uwave_plus - uwave_minus) / (2.0 * delta_c) );
352 
353  }
354  }
355 
356 }
double Radius(const int i, const int j) const
double SamplePoint2(const unsigned int i) const
virtual std::complex< double > Element(const double radius, const double cosine_rz, const double e_pion)=0
unsigned int GetNDensities(void) const
std::complex< double > cdouble
void set(unsigned int i, unsigned int j, const std::complex< double > &value)
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fP_pi

Member Data Documentation

current_t genie::alvarezruso::AlvarezRusoCOHPiPDXSec::current
private

Definition at line 114 of file AlvarezRusoCOHPiPDXSec.h.

bool genie::alvarezruso::AlvarezRusoCOHPiPDXSec::debug_
private

Definition at line 108 of file AlvarezRusoCOHPiPDXSec.h.

unsigned int genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fA
private

Definition at line 111 of file AlvarezRusoCOHPiPDXSec.h.

ARConstants* genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fConstants
private

Definition at line 122 of file AlvarezRusoCOHPiPDXSec.h.

double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fE_l
private

Definition at line 130 of file AlvarezRusoCOHPiPDXSec.h.

double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fE_nu
private

Definition at line 129 of file AlvarezRusoCOHPiPDXSec.h.

double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fF_cross_delta
private

Definition at line 156 of file AlvarezRusoCOHPiPDXSec.h.

double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fF_cross_nucleon
private

Definition at line 157 of file AlvarezRusoCOHPiPDXSec.h.

double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fF_direct_delta
private

Definition at line 154 of file AlvarezRusoCOHPiPDXSec.h.

double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fF_direct_nucleon
private

Definition at line 155 of file AlvarezRusoCOHPiPDXSec.h.

double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fg_factor
private

Definition at line 151 of file AlvarezRusoCOHPiPDXSec.h.

std::complex<double> genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fJ_hadronic[4]
private

Definition at line 164 of file AlvarezRusoCOHPiPDXSec.h.

double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fLastE_pi
private

Definition at line 135 of file AlvarezRusoCOHPiPDXSec.h.

flavour_t genie::alvarezruso::AlvarezRusoCOHPiPDXSec::flavour
private

Definition at line 116 of file AlvarezRusoCOHPiPDXSec.h.

double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fM_l
private

Definition at line 150 of file AlvarezRusoCOHPiPDXSec.h.

double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fM_pi
private

Definition at line 149 of file AlvarezRusoCOHPiPDXSec.h.

ARSampledNucleus* genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fNucleus
private

Definition at line 124 of file AlvarezRusoCOHPiPDXSec.h.

formfactors_t genie::alvarezruso::AlvarezRusoCOHPiPDXSec::formfactors
private

Definition at line 120 of file AlvarezRusoCOHPiPDXSec.h.

ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fP_cross
private

Definition at line 145 of file AlvarezRusoCOHPiPDXSec.h.

ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fP_direct
private

Definition at line 144 of file AlvarezRusoCOHPiPDXSec.h.

ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fP_l
private

Definition at line 140 of file AlvarezRusoCOHPiPDXSec.h.

ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fP_n_i
private

Definition at line 142 of file AlvarezRusoCOHPiPDXSec.h.

ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fP_n_o
private

Definition at line 143 of file AlvarezRusoCOHPiPDXSec.h.

ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fP_nu
private

Definition at line 139 of file AlvarezRusoCOHPiPDXSec.h.

ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fP_pi
private

Definition at line 141 of file AlvarezRusoCOHPiPDXSec.h.

double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fPhi
private

Definition at line 133 of file AlvarezRusoCOHPiPDXSec.h.

ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fQ
private

Definition at line 138 of file AlvarezRusoCOHPiPDXSec.h.

unsigned int genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fSampling
private

Definition at line 112 of file AlvarezRusoCOHPiPDXSec.h.

double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fTheta_l
private

Definition at line 131 of file AlvarezRusoCOHPiPDXSec.h.

double genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fTheta_pi
private

Definition at line 132 of file AlvarezRusoCOHPiPDXSec.h.

ARWavefunction* genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fUwave
private

Definition at line 160 of file AlvarezRusoCOHPiPDXSec.h.

ARWavefunction* genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fUwaveDr
private

Definition at line 161 of file AlvarezRusoCOHPiPDXSec.h.

ARWavefunction* genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fUwaveDtheta
private

Definition at line 162 of file AlvarezRusoCOHPiPDXSec.h.

ARWFSolution* genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fWfsolution
private

Definition at line 126 of file AlvarezRusoCOHPiPDXSec.h.

unsigned int genie::alvarezruso::AlvarezRusoCOHPiPDXSec::fZ
private

Definition at line 110 of file AlvarezRusoCOHPiPDXSec.h.

nutype_t genie::alvarezruso::AlvarezRusoCOHPiPDXSec::nutype
private

Definition at line 118 of file AlvarezRusoCOHPiPDXSec.h.


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