GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AlvarezRusoCOHPiPDXSec.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Daniel Scully ( d.i.scully \at warwick.ac.uk)
7  University of Warwick
8 
9  Based on Fortran code by Luis Alvarez-Ruso.
10 */
11 //____________________________________________________________________________
12 
13 // C++
14 #include <iostream>
15 #include <iomanip>
16 #include <sstream>
17 #include <string>
18 #include <cstdlib>
19 #include <complex>
20 
21 // Root
22 #include <TVector3.h>
23 #include <TMath.h>
24 #include <Math/SMatrix.h>
25 #include <Math/SVector.h>
26 #include <Math/LorentzVector.h>
27 
28 //Genie
32 
33 //AlvarezRuso
39 
40 using namespace genie::constants;
41 
42 typedef std::complex<double> cdouble;
43 typedef ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > LorentzVector;
44 typedef ROOT::Math::SVector< cdouble , 4> CVector;
45 
46 namespace genie {
47 namespace alvarezruso {
48 
49 AlvarezRusoCOHPiPDXSec::AlvarezRusoCOHPiPDXSec(unsigned int Z_, unsigned int A_, const current_t current_,
50  const flavour_t flavour_, const nutype_t nutype_,const formfactors_t ff_)
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 }
72 
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 }
82 
83 
84 double 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_)
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 }
141 
142 
143 
144 cdouble AlvarezRusoCOHPiPDXSec::H(unsigned int i, unsigned int j) const
145 {
146  cdouble H_ = ( conj(fJ_hadronic[i]) * fJ_hadronic[j] );
147  return H_;
148 }
149 
150 
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 }
177 
178 
179 /* *********************************************************************
180  * Vertex factor for Delta/Nucleon decaying to Pi + Nucleon
181  * formerly offshell()
182  */
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 }
201 
202 /* *********************************************************************
203  * Fill the LorentzVectors with values based on the values from the
204  * kinematics
205  */
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 }
228 
229 
230 /* *********************************************************************
231  *
232  */
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 }
263 
264 
265 
266 
267 
268 /* *********************************************************************
269  * Fill values based on the current
270  */
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 }
288 
289 
290 
291 
292 /*
293  * Solve the wavefunctions
294  */
295 
296 /// This is only a function of the nucleus and pion momentum/energy
297 /// so if neither of those have changed there is no need to re-calculate
298 /// the wavefunction values.
299 /// Such a caching has not been implemented here yet!
300 
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 }
357 
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 }
372 
373 double AlvarezRusoCOHPiPDXSec::DeltaWidthPauliBlocked(LorentzVector delta_momentum, double density)
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 }
430 
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 }
455 
456 // 1.1.1.1
457 /*
458  * Calculate the three-momentum of a pion coming from the decay of a
459  * Delta resonance.
460  */
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 }
475 
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 }
495 
497 {
498  double result = (0.04/fConstants->HBar())*(density/fConstants->Rho0());
499  return result;
500 }
501 
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 }
546 
547 double AlvarezRusoCOHPiPDXSec::DeltaSelfEnergyConstant(double a, double b, double c, double E)
548 {
549  double x = (E / fM_pi) - 1.0;
550  return (a*x*x + b*x + c);
551 }
552 
554 //calculates the nuclear current
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 }
971 
972 
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 }
1002 
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 }
1013 
1015 {
1016  return *fConstants;
1017 }
1018 
1020 {
1021  return *fNucleus;
1022 }
1023 
1024 } // namespace alvarezruso
1025 } // namespace genie
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fP_nu
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fQ
Wave function class for AlvarezRuso Coherent pion production xsec.
ROOT::Math::SVector< cdouble, 4 > CVector
std::complex< double > cdouble
#define pERROR
Definition: Messenger.h:59
double DeltaWidthFree(ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > delta_momentum)
static constexpr double s
Definition: Units.h:95
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)
std::complex< double > DeltaPropagatorInMed(ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > delta_momentum)
static constexpr double fs
Definition: Units.h:100
double Radius(const int i, const int j) const
double DeltaWidthPauliBlocked(ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > delta_momentum, double density)
double PionMomentumCM(ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > delta_momentum)
Eikonal wavefunction solution for Alvarez-Ruso Coherent Pion Production xsec.
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > LorentzVector
static constexpr double b
Definition: Units.h:78
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
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)
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)
double PNVertexFactor(ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > momentum, double mass)
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_)
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static constexpr double A
Definition: Units.h:74
double SamplePoint2(const unsigned int i) const
const double a
virtual std::complex< double > Element(const double radius, const double cosine_rz, const double e_pion)=0
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fP_l
#define pWARN
Definition: Messenger.h:60
double DeltaSelfEnergyConstant(double a, double b, double c, double E)
std::complex< double > H(unsigned int i, unsigned int j) const
Nucleus class for Alvarez-Ruso Coherent Pion Production xsec.
unsigned int GetNDensities(void) const
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
void set(unsigned int i, unsigned int j, const std::complex< double > &value)
double DensityOfCentres(const int i, const int j) const
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > fP_pi
static constexpr double m
Definition: Units.h:71
double SamplePoint1(const unsigned int i) const
std::complex< double > NucleonPropagator(ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > nucleon_momentum)