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