ARSampledNucleus.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 #include "ARSampledNucleus.h"
16 
17 #include <string>
18 #include <iostream>
19 #include <cstdlib>
20 #include <sstream>
21 #include <fstream>
22 #include <istream>
23 #include <complex>
24 
25 #include <TSystem.h>
26 #include <TF1.h>
27 #include <TMath.h>
28 #include <Math/GaussLegendreIntegrator.h>
29 
33 
34 //
35 // Equation/Table numbers refer to:
36 // J. Nieves and E. Oset
37 // A Theoretical approach to pionic atoms and the problem of anomalies
38 // Nuclear Physics A 554 (1993) 509-553
39 //
40 
41 using namespace genie::constants;
42 
43 namespace genie {
44 namespace alvarezruso {
45 
46 double ARSampledNucleus::mean_radius_squared = 0.69; // fm (CA)
47 
48 ARSampledNucleus::ARSampledNucleus(unsigned int ZNumber, unsigned int ANumber, unsigned int fSampling_):
49  fZ(ZNumber),
50  fA(ANumber),
51  fSampling(fSampling_)
52 {
54 
55  fDensities = NULL;
56  fDensitiesOfCentres = NULL;
57  fRadii = NULL;
58  fSample_points_1 = NULL;
59  fSample_points_2 = NULL;
60  fSample_weights_1 = NULL;
61  fSample_weights_2 = NULL;
62 
63  if(fA>20) { // fermi gas
64  if (fA == 27) { fNucRadius = 3.07; fDiffuseness = 0.52; } // aluminum
65  else if (fA == 28) { fNucRadius = 3.07; fDiffuseness = 0.54; } // silicon
66  else if (fA == 40) { fNucRadius = 3.53; fDiffuseness = 0.54; } // argon
67  else if (fA == 56) { fNucRadius = 4.10; fDiffuseness = 0.56; } // iron
68  else if (fA == 208) { fNucRadius = 6.62; fDiffuseness = 0.55; } // lead
69  else {
70  fNucRadius = pow(fA,0.35); fDiffuseness = 0.54;
71  } //others
72  fUseHarmonicOscillator = false;
73  }
74  else {
75  if (fA == 7) { fNucRadius = 1.77 ; fDiffuseness = 0.327;} // lithium
76  else if (fA == 12) { fNucRadius = 1.692 ; fDiffuseness = 1.083;} // carbon
77  else if (fA == 14) { fNucRadius = 1.76 ; fDiffuseness = 1.23; } // nitrogen
78  else if (fA == 16) { fNucRadius = 1.83 ; fDiffuseness = 1.54; } // oxygen
79  else if (fA <= 4 ) { fNucRadius = 1.344 ; fDiffuseness = 0.00; } // helium
80  else {
81  fNucRadius=1.75; fDiffuseness=-0.4+.12*fA;
82  } //others- fDiffuseness=0.08 if A=4
83 
85  }
86 
88 
89  // Calculate number densities (density of 'centres'):
90  double diff2 = fDiffuseness*fDiffuseness;
93  double x = (fDiffuseness * fNucRadiusSq) / ((1 + (1.5*fDiffuseness)) * (fRadiusCentres*fRadiusCentres));
94  fDiffusenessCentres = 2.*x / (2. - 3.*x ) ;
95  }
96  else { //fermi liquid
97  double numerator = 5.0 * mean_radius_squared * fNucRadius;
98  double denominator = (15.0 * (fNucRadiusSq)) + (7.0 * kPi2 * fDiffuseness*fDiffuseness);
99  fRadiusCentres = fNucRadius + (numerator / denominator);
100 
101  numerator = (fNucRadiusSq*fNucRadius) + (kPi2 * diff2 * fRadiusCentres) - (fRadiusCentres*fRadiusCentres*fRadiusCentres);
102  denominator = kPi2 * fRadiusCentres;
103 
104  fDiffusenessCentres = sqrt( numerator / denominator );
105  }
106 
107  this->Fill();
108 }
109 
111 {
112  for(unsigned int i = 0; i != fNDensities; ++i)
113  {
114  if (fDensities && fDensities [i]) delete[] fDensities[i];
116  if (fRadii && fRadii [i]) delete[] fRadii [i];
117  }
118 
119  if (fDensities ) delete[] fDensities;
121  if (fRadii ) delete[] fRadii ;
122  if (fSample_points_1 ) delete[] fSample_points_1;
123  if (fSample_points_2 ) delete[] fSample_points_2;
124  if (fSample_weights_1 ) delete[] fSample_weights_1;
125  if (fSample_weights_2 ) delete[] fSample_weights_2;
126 }
128 {
129 
130  fR_max = 3.0 * TMath::Power(this->A(), (1.0/3.0));
131 
132 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
133  LOG("AR_PiWFunction_Table", pDEBUG)<< "N:: fR_max = " << fR_max
134  << "N:: z = " << fZ
135  << "N:: a = " << fA
136 #endif
137 
138  this->FillSamplePoints();
139  this->FillDensities();
140 }
141 
142 double ARSampledNucleus::Density(const int i, const int j) const
143 {
144  return fDensities[i][j];
145 }
146 
147 double ARSampledNucleus::DensityOfCentres(const int i, const int j) const
148 {
149  return fDensitiesOfCentres[i][j];
150 }
151 
152 double ARSampledNucleus::Radius(const int i, const int j) const
153 {
154  return fRadii[i][j];
155 }
156 
158 {
159  if (fSample_points_1) delete[] fSample_points_1;
160  if (fSample_points_2) delete[] fSample_points_2;
161  if (fSample_weights_1) delete[] fSample_weights_1;
162  if (fSample_weights_2) delete[] fSample_weights_2;
163 
164  fSample_points_1 = new double[fNDensities];
165  fSample_points_2 = new double[fNDensities];
166 
167  fSample_weights_1 = new double[fNDensities];
168  fSample_weights_2 = new double[fNDensities];
169  unsigned int decoy;
170 
173 
174 }
175 
177 {
178  double r;
179 
180  for(unsigned int i = 0; i != fNDensities; ++i)
181  {
182  if (fDensities && fDensities [i]) delete[] fDensities [i];
184  if (fRadii && fRadii [i]) delete[] fRadii [i];
185  }
186  if (fDensities ) delete[] fDensities;
188  if (fRadii ) delete[] fRadii ;
189 
190  fDensities = new double*[fNDensities];
191  fDensitiesOfCentres = new double*[fNDensities];
192  fRadii = new double*[fNDensities];
193 
194  for(unsigned int i = 0; i != fNDensities; ++i)
195  {
196  fDensities [i] = new double[fNDensities];
197  fDensitiesOfCentres[i] = new double[fNDensities];
198  fRadii [i] = new double[fNDensities];
199 
200  for(unsigned int j = 0; j != fNDensities; ++j)
201  {
203  fRadii[i][j] = r;
204  fDensities [i][j]=this->CalcMatterDensity(r);
206  }
207  }
208 }
209 
210 unsigned int ARSampledNucleus::GetSampling (void) const
211 {
212  return fSampling;
213 }
214 unsigned int ARSampledNucleus::GetNDensities(void) const
215 {
216  return fNDensities;
217 }
218 
219 double ARSampledNucleus::CalcDensity(double r, double nuc_rad, double nuc_diff) const
220 {
221  double dens_rel;
223  dens_rel = (1.0 + nuc_diff*r*r/(nuc_rad*nuc_rad)) * exp(-r*r/(nuc_rad*nuc_rad));
224  }
225  else { //fermi liquid
226  dens_rel = 1.0 / (1.0 + exp((r - nuc_rad)/nuc_diff));
227  }
228 
229  //~ double dens_0 = utils::nuclear::Density(nuc_rad, fA);
230  double dens_0 = Density0(fA,nuc_rad,nuc_diff);
231 
232  return dens_0 * dens_rel;
233 }
234 
236  unsigned int number, // A
237  double radius, // R
238  double diffuseness // a
239  ) const
240 {
241  double result = 0.0;
243  {
244  double u = fR_max/radius;
245  double dterm = (3*diffuseness+2);
246  double term1 = TMath::Sqrt(TMath::Pi())*dterm*TMath::Power(radius,3)*TMath::Exp(u*u)*TMath::Erf(u);
247  double term2 = 2*fR_max*(dterm*radius*radius+2*diffuseness*radius*radius);
248  result = (0.5)*TMath::Pi()*TMath::Exp(-u*u)*(term1 - term2);
249  }
250  else
251  {
252  //Probably faster to do this via the integral for the moment as ROOT doesn't have a builtin
253  //PolyLog function
254  TF1* f = this->Density0Function();
255  f->SetParameter(0, diffuseness);
256  f->SetParameter(1, radius);
257  result = f->Integral(0.0,fR_max);
258  delete f;
259  }
260 
261  return ( number / result );
262 }
263 
265 {
266  return (new TF1("density0function", Density0FunctionFermiLiquid, 0.0, fR_max, 2));
267 }
268 
269 Double_t ARSampledNucleus::Density0FunctionFermiLiquid(Double_t* r_, Double_t* parameters)
270 {
271  double r = (*r_);
272  double diffuseness = parameters[0];
273  double radius = parameters[1];
274 
275  double dens = 1.0 / (1.0 + exp((r - radius)/diffuseness));
276 
277  return 4.0*TMath::Pi()*r*r*dens;
278 }
279 
281 {
282  return this->CalcDensity(r,fNucRadius,fDiffuseness);
283 }
284 
286 {
288 }
289 
290 } //namespace alvarezruso
291 } //namespace genie
double Density(const int i, const int j) const
Basic constants.
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
double Radius(const int i, const int j) const
double CalcMatterDensity(double r) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
static Double_t Density0FunctionFermiLiquid(Double_t *r, Double_t *parameters)
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
double Density0(unsigned int number, double diffuseness, double radius) const
const double j
Definition: BetheBloch.cxx:29
Double_t radius
double CalcDensity(double radius, double nuc_rad, double nuc_diff) const
unsigned int GetNDensities(void) const
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
unsigned int GetSampling(void) const
TRandom3 r(0)
double CalcNumberDensity(double r) const
double DensityOfCentres(const int i, const int j) const
void SGNR(const double a, const double b, const unsigned int n, const unsigned int nsamp, double *x, unsigned int &np, double *w)
static const double kPi2
Definition: Constants.h:39
#define pDEBUG
Definition: Messenger.h:64