GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AREikonalSolution.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 //____________________________________________________________________________
10 
11 
12 #include <TMath.h>
13 
14 #include <cstdlib>
15 
19 
20 typedef std::complex<double> cdouble;
21 
22 namespace genie {
23 namespace alvarezruso {
24 
25 cdouble AREikonalSolution::Element(const double radius, const double cosine_rz,
26  const double e_pion)
27 {
28 
29  const double mpik = this->Parent()->GetPiMass();
30  const double mpi = this->Con()->PiPMass();
31  const double hb = this->Con()->HBar() * 1000.0;
32 
33  const double r = (radius);
34 
35  const double ekin = (e_pion - mpik) * hb;
36 
37  const double cosa = cosine_rz;
38 
39  const double rmax = this->Nucleus()->RadiusMax();
40 
41  const unsigned int nz = 1;
42 
43  const double za = r * cosa;
44  const double be = r * TMath::Sqrt(1.0 - cosa*cosa);
45  const double omepi = ekin / hb + mpi;
46  const double ppim = TMath::Sqrt(omepi*omepi - mpi*mpi);
47  unsigned int sampling = (this->Nucleus())->GetSampling();
48 
49  double absiz[sampling];
50  double decoy[sampling];
51 
52  unsigned int junk;
53 
54  integrationtools::SGNR(za, rmax, nz, sampling, absiz, junk, decoy);
55 
56  //do i=1,nzs
57 // cdouble ordez[sampling];
58  cdouble * ordez = new cdouble[sampling]; // CA
59  double zp, rp;
60  cdouble piself;
61 
62  unsigned int A = fNucleus->A();
63  unsigned int Z = fNucleus->Z();
64 
65  for(unsigned int i = 0; i != sampling; ++i)
66  {
67  // Sample point in nucleus
68  zp = absiz[i];
69 
70  // Radius in nucleus
71  rp = TMath::Sqrt( be*be + zp*zp );
72 
73  // Get nuclear densities
74  double dens_cent = fNucleus->CalcNumberDensity(rp);
75  double dens_p_cent = dens_cent * Z / A ;
76  double dens_n_cent = dens_cent * (A-Z)/A;
77 
78  // Calculate pion self energy
79  piself = this->PionSelfEnergy(dens_p_cent, dens_n_cent, omepi, ppim);
80 
81  // Optical potential at each point in the nucleus
82  ordez[i] = piself / 2.0 / ppim;
83 
84  }
85 
86  //Integrate the optical potential through the nucleus
87  cdouble resu = integrationtools::RGN1D(za, rmax, nz, sampling, ordez);
88 
89  // Eikonal approximation to the wave function
90  cdouble uwaveik = exp( - cdouble(0,1) * ( ppim*za + resu ) );
91 
92  delete [] ordez; // CA
93 
94  return uwaveik;
95 }
96 
97 
98 cdouble AREikonalSolution::PionSelfEnergy(const double rhop_cent, const double rhon_cent, const double omepi, const double ppim)
99 {
100  const double rho0 = this->Con()->Rho0();
101  const double mn = this->Con()->NucleonMass();
102  const double hb = this->Con()->HBar() * 1000.0;
103  const double pi = constants::kPi;
104  const double fs = this->Con()->DeltaNCoupling();
105  const double mdel = this->Con()->DeltaPMass();
106  const double mpi = this->Con()->PiPMass();
107  const double fs_mpi2 = fs*fs/(mpi*mpi);
108 
109 
110  const cdouble ui(0,1);
111 
112  const double resig = -53.0/hb;
113  const double rho = rhop_cent + rhon_cent;
114  const double rat = rho / rho0;
115  const double pf = TMath::Power( (3.*pi*pi/2.*rho), (1.0/3.0) );
116 
117  const double sdel = mn*mn + mpi*mpi + 2.*omepi*(mn+3./5.*pf*pf/2./mn);
118  const double sqsdel = TMath::Sqrt(sdel);
119 
120  double gamdpb, imsig;
121  this->Deltamed(sdel, pf, rat, gamdpb, imsig, ppim, omepi);
122 
123  const cdouble pe = -1./6./pi*fs_mpi2*
124  ( rhop_cent/(sqsdel-mdel-resig*(2.*rhon_cent/rho0)+ui*(gamdpb/2.-imsig)) +
125  1./3.*rhon_cent/(sqsdel-mdel-resig*2./3.*(2.*rhon_cent+rhop_cent)/rho0+ui*(gamdpb/2.-imsig)) +
126  rhon_cent/(-sqsdel-mdel+2.*mn-resig*(2.*rhop_cent/rho0))+
127  1./3.*rhop_cent/(-sqsdel-mdel+2.*mn-resig*2./3.*(2.*rhop_cent+rhon_cent)/rho0) );
128 
129  const cdouble efe = 4.*pi*mn*mn/sdel*pe/(1.+4.*pi*0.63*pe);
130 
131  const cdouble piself = -1.0*efe*(1.-1./2.*omepi/mn)*ppim*ppim/(1.+efe*(1.-1./2.*omepi/mn));
132 
133  return piself;
134 }
135 
136 
137 void AREikonalSolution::Deltamed(const double sdel, const double pf, const double rat, double& gamdpb, double& imsig, const double ppim, const double omepi)
138 {
139  unsigned int iapr = 1; // approximation chosen to calculate gamdpb
140 
141  // Calculation of the pauli-blocked width
142  double gamdfree = this->Gamd(sdel);
143 
144  if( gamdfree == 0.0)
145  {
146  gamdpb = 0.0;
147  }
148  else
149  {
150  double f;
151  if( iapr == 1 )
152  {
153  // Approximation from Nieves et al. NPA 554(93)554
154  const double r = this->Qcm(sdel) / pf;
155  if( r > 1.0 )
156  {
157  // f=1.+(-2./5./r**2+9./35./r**4-2./21./r**6)
158  f = 1.0 + ( -2.0 / 5.0 / (r*r) + 9.0 / 35.0 / (r*r*r*r) - 2.0 / 21.0 / (r*r*r*r*r*r) );
159  }
160  else
161  {
162  //f=34./35.*r-22./105*r**3
163  f = 34.0 / 35.0 * r - 22.0 / 105 * (r*r*r);
164  }
165  }
166  else
167  {
168  //Approximation from Garcia-Recio, NPA 526(91)685
169 
170  const double mn = this->Con()->NucleonMass();
171  const double wd = TMath::Sqrt(sdel); // Delta inv. mass
172  const double ef = TMath::Sqrt(mn*mn + pf*pf); // Fermi energy
173  const double kd = ppim; // modulus of the Delta 3-momentum in Lab.
174  const double pn = this->Qcm(sdel); // nucleon(and pion) momentum in C.M.
175  const double en = TMath::Sqrt(mn*mn + pn*pn);
176 
177  f = ( kd * pn + TMath::Sqrt(sdel+kd*kd) * en - ef * wd ) / (2.0 * kd * pn);
178  if (f < 0.0) f = 0.0;
179  if (f > 1.0) f = 1.0;
180  }
181  gamdpb = gamdfree * f;
182  }
183 
184  //Calculation of the delta selfenergy
185 
186  // Imaginary part: using Oset, Salcedo, NPA 468(87)631
187  // Using eq. (3.5) to relate the energy of the delta with the pion energy used
188  // in the parametrization
189 
190  // Prescriptions for the effective pion energy
191  // nucleon at rest
192  // ! ome=p(0)-mn
193  // ! ome=(sdel-mn**2-mpi**2)/2./mn
194  // nucleon with an average momentum
195  // ! ekp=3./5.*pf**2/2./mn
196  // ! ome=p(0)-mn-ekp
197  // ! ome=(sdel-mn**2-mpi**2-ekp**2)/2./(mn+ekp)
198  double ome = omepi;
199  double mpi = this->Parent()->GetPiMass();
200 
201  // The parameterization is valid for 85 MeV < tpi < 315. outside we take a contant values
202  const double hb = this->Con()->HBar() * 1000.0;
203  if( ome <= (mpi + 85.0 / hb) ) ome = mpi + 85.0 / hb;
204  if( ome >= (mpi + 315.0 / hb) ) ome = mpi + 315.0 / hb;
205 
206  // The parameterization of Oset, Salcedo, with ca3 extrapolated to zero at low kin. energies
207  const double cq = this->Cc(-5.19,15.35,2.06,ome)/hb;
208  const double ca2 = this->Cc(1.06,-6.64,22.66,ome)/hb;
209  double ca3 = this->Cc(-13.46,46.17,-20.34,ome)/hb;
210  const double alpha= this->Cc(0.382,-1.322,1.466,ome);
211  const double beta = this->Cc(-0.038,0.204,0.613,ome);
212  const double gamma=2.*beta;
213  if( ome <= (mpi + 85.0/hb) ) ca3 = this->Cc(-13.46,46.17,-20.34,(mpi+85./hb))/85.*(ome-mpi);
214 
215  imsig = - ( cq*TMath::Power(rat,alpha) + ca2*TMath::Power(rat, beta) + ca3*TMath::Power(rat,gamma) );
216 }
217 
218 
219 double AREikonalSolution::Cc(const double a, const double b, const double c, const double ome)
220 {
221  double mpi = this->Parent()->GetPiMass();
222  const double x = ome / mpi - 1.0;
223  return (a*x*x + b*x + c);
224 }
225 
226 
227 double AREikonalSolution::Gamd(const double s)
228 {
229  // Delta -> N pi
230  const double mpi = this->Con()->PiPMass();
231  const double mn = this->Con()->NucleonMass();
232  if( s <= (mn+mpi)*(mn+mpi) )
233  {
234  return 0.0;
235  }
236  else
237  {
238  double fs_mpi2 = this->Con()->DeltaNCoupling()/mpi;
239  fs_mpi2 *= fs_mpi2;
240  const double qcm = this->Qcm(s);
241  return 1.0 / (6.0*constants::kPi)*fs_mpi2*mn*(qcm*qcm*qcm)/TMath::Sqrt(s);
242  }
243 }
244 
245 
246 double AREikonalSolution::Qcm(const double s)
247 {
248  // Returns the 3-momentum of the pion formed after the decay of a
249  // resonance (R->N pi) of inv. mass s in the rest frame of the resonance
250  const double mpi = this->Con()->PiPMass();;
251  const double mn = this->Con()->NucleonMass();
252  return TMath::Sqrt((s-mpi*mpi-mn*mn)*(s-mpi*mpi-mn*mn) - 4.*mpi*mpi*mn*mn)/2.0/TMath::Sqrt(s);
253 }
254 
255 AREikonalSolution::AREikonalSolution(bool debug, AlvarezRusoCOHPiPDXSec* parent): ARWFSolution(debug), parent_(parent)
256 {
257  if( debug_ ) std::cerr << "AREikonalSolution::AREikonalSolution" << std::endl;
258  this->fNucleus = &(this->Parent()->GetNucleus());
259  this->constants_ = &(this->Parent()->GetConstants());
260  owns_constants = false;
261 }
262 
263 AREikonalSolution::AREikonalSolution(bool debug, ARSampledNucleus* nucl): ARWFSolution(debug), fNucleus(nucl)
264 {
265  if( debug_ ) std::cerr << "AREikonalSolution::AREikonalSolution" << std::endl;
266  this->constants_ = new ARConstants();
267  owns_constants = true;
268 }
269 
271  if (owns_constants) delete constants_;
272 }
273 
275 {
276  if(false) std::cout << "Hi!" << std::endl;
277 }
278 
279 } //namespace alvarezruso
280 } //namespace genie
AREikonalSolution(bool debug, AlvarezRusoCOHPiPDXSec *parent)
std::complex< double > cdouble
std::complex< double > PionSelfEnergy(const double rhop_cent, const double rhon_cent, const double omepi, const double ppim)
static constexpr double s
Definition: Units.h:95
static constexpr double fs
Definition: Units.h:100
static constexpr double b
Definition: Units.h:78
Abstract base class for Alvarez-Ruso wavefunction solution.
Definition: ARWFSolution.h:32
static constexpr double A
Definition: Units.h:74
const double a
double Cc(const double a, const double b, const double c, const double ome)
Nucleus class for Alvarez-Ruso Coherent Pion Production xsec.
AlvarezRusoCOHPiPDXSec * Parent()
virtual std::complex< double > Element(const double radius, const double cosine_rz, const double e_pion)
double CalcNumberDensity(double r) const
std::complex< double > RGN1D(const double A, const double B, const unsigned int N, const unsigned int nsamp, const std::complex< double > CF[])
void SGNR(const double a, const double b, const unsigned int n, const unsigned int nsamp, double *x, unsigned int &np, double *w)
void Deltamed(const double sdel, const double pf, const double rat, double &gamdpb, double &imsig, const double ppim, const double omepi)