EarthModel.cxx
Go to the documentation of this file.
1 ///
2 /// \file EarthModel.cxx
3 /// \brief Provide interface to different Earth models
4 /// \author messier@indiana.edu
5 /// \version $Id: EarthModel.cxx,v 1.2 2012-07-01 22:21:59 gsdavies Exp $
6 ///
7 #include "OscLib/EarthModel.h"
8 
9 #include <cstdlib>
10 #include <iostream>
11 #include <cmath>
12 #include <list>
13 using namespace osc;
14 
15 static const int kPREM = 0;
16 static const int kSTACEY = 1;
17 
18 namespace{
19  namespace util{
20  template<class T> T sqr(T x){return x*x;}
21  template<class T> T pythag(T x, T y){return sqrt(sqr(x)+sqr(y));}
22  }
23 }
24 
25 EarthModel::EarthModel(const char* which, double tol)
26 {
27  if (std::string("PREM") ==which) this->InitPREM();
28  else if (std::string("STACEY")==which) this->InitStacey();
29  // else if (std::string("ModelX")==which) this->InitModelX();
30  else {
31  std::cerr << __FILE__ << ":" << __LINE__
32  << " Model '" << which << "' is not supported."
33  << " Currently only PREM is supported." << std::endl;
34  abort();
35  }
36  this->MakeLayers(tol);
37 }
38 
39 //......................................................................
40 ///
41 /// Return the electron number density at radius r
42 ///
43 /// \param r - radius in units of km
44 /// \returns electron densty in units of moles/cm^3
45 ///
46 double EarthModel::Ne(double r)
47 {
48  return this->ZoverA(r)*this->Density(r);
49 }
50 
51 //......................................................................
52 ///
53 /// Earth density in g/cm^3 as a function of radius in km
54 ///
55 double EarthModel::Density(double r)
56 {
57  switch (fModelID) {
58  case kPREM: return this->DensityPREM(r);
59  case kSTACEY: return this->DensityStacey(r);
60  default: abort();
61  }
62  return 0.0;
63 }
64 
65 //......................................................................
66 ///
67 /// Initialize the Earth model to "PREM"
68 ///
69 /// [*] Dziewonski and Anderson, Physics of the Earth and Planetary
70 /// Interiors, 25 (1981) 297-356. Table I.
71 ///
73 {
74  fModelID = kPREM;
75  fRouterCore = 3480.0;
76  fRearth = 6371.0;
77 
78  fRregion.push_back(0.0);
79  fRregion.push_back(1221.5); // Inner core
80  fRregion.push_back(3480.0); // Outer core
81  fRregion.push_back(3630.0); // D''
82  fRregion.push_back(5701.0); // Lower mantle
83  fRregion.push_back(6151.0); // Transition zone
84  fRregion.push_back(6291.0); // Low velocity zone
85  fRregion.push_back(6346.6); // LID
86  fRregion.push_back(6368.0); // Crust
87  fRregion.push_back(6371.0); // Ocean
88 }
89 
90 //......................................................................
91 ///
92 /// The PREM Earth density profile in units of g/cm^3 as a function of
93 /// radius r (km)
94 ///
95 /// [*] Dziewonski and Anderson, Physics of the Earth and Planetary
96 /// Interiors, 25 (1981) 297-356. Table I.
97 ///
98 /// \param r - radius in km
99 /// \returns density in g/cm^3
100 ///
101 double EarthModel::DensityPREM(double r)
102 {
103  double x = r/fRearth;
104 
105  // The regions are checked in order from largest to smallest
106  if (r>1221.5&&r<=3480.0) return 12.5815+x*(-1.2638+x*(-3.6426+x*(-5.5281)));
107  if (r>3480.0&&r<=5701.0) return 7.9565+x*(-6.4761+x*(5.5283+x*(-3.0807)));
108  if (r<=1221.5) return 13.0885-8.8381*x*x;
109  if (r>5771.0&&r<=5971.0) return 11.2494-8.0298*x;
110  if (r>5971.0&&r<=6151.0) return 7.1089-3.8045*x;
111  if (r>6151.0&&r<=6291.0) return 2.6910+0.6924*x;
112  if (r>5701.0&&r<=5771.0) return 5.3197-1.4836*x;
113  if (r>6291.0&&r<=6346.6) return 2.6910+0.6924*x;
114  if (r>6346.6&&r<=6356.0) return 2.90;
115  if (r>6356.0&&r<=6368.0) return 2.60;
116  if (r>6368.0&&r<=6371.0) return 1.020;
117  return 0.0;
118 }
119 
120 //......................................................................
121 ///
122 /// Initialize the Earth model to Stacey:
123 ///
124 /// [*] F.D. Stacey, Physics of the Earth (Wiley, New York, 1969)
125 ///
126 /// \param r - radius in km
127 /// \returns density in g/cm^3
128 ///
130 {
131  fModelID = kSTACEY;
132  fRouterCore = 6371.0-2878.0;
133  fRearth = 6371.0;
134 
135  fRregion.push_back(0.0);
136  fRregion.push_back(6371.0-5121.0); // Inner core
137  fRregion.push_back(6371.0-2878.0); // Outer core
138  fRregion.push_back(6371.0-650.0); // Lower mantle
139  fRregion.push_back(6371.0-350.0); // Middle mantle
140  fRregion.push_back(6371.0-15.0); // Upper mantle
141  fRregion.push_back(6371.0); // Crust
142 }
143 
144 //......................................................................
145 ///
146 /// Earth density profile according to
147 ///
148 /// [*] F.D. Stacey, Physics of the Earth (Wiley, New York, 1969)
149 ///
151 {
152  static double crustD[2] = { 0.0,
153  15.0};
154  static double upperMantleD[6] = { 15.0,
155  60.0,
156  100.0,
157  200.0,
158  300.0,
159  350.0};
160  static double middleMantleD[6] = {350.0,
161  400.0,
162  413.0,
163  500.0,
164  600.0,
165  650.0};
166  static double lowerMantleD[14] = {650.0,
167  800.0,
168  984.0,
169  1000.0,
170  1200.0,
171  1400.0,
172  1600.0,
173  1800.0,
174  2000.0,
175  2200.0,
176  2400.0,
177  2600.0,
178  2800.0,
179  2878.0};
180  static double outerCoreD[14] = { 2878.0,
181  3000.0,
182  3200.0,
183  3400.0,
184  3600.0,
185  3800.0,
186  4000.0,
187  4200.0,
188  4400.0,
189  4600.0,
190  4800.0,
191  4982.0,
192  5000.0,
193  5121.0};
194  static double innerCoreD[8] ={ 5121.0,
195  5200.0,
196  5400.0,
197  5600.0,
198  5800.0,
199  6000.0,
200  6200.0,
201  6371.0};
202  static double crustRho[2] = { 2.840,
203  2.840};
204  static double upperMantleRho[6] = { 3.313,
205  3.332,
206  3.348,
207  3.387,
208  3.424,
209  3.441};
210  static double middleMantleRho[6] = {3.700,
211  3.775,
212  3.795,
213  3.925,
214  4.075,
215  4.150};
216  static double lowerMantleRho[14] = {4.200,
217  4.380,
218  4.529,
219  4.538,
220  4.655,
221  4.768,
222  4.877,
223  4.983,
224  5.087,
225  5.188,
226  5.288,
227  5.387,
228  5.487,
229  5.527};
230  static double outerCoreRho[14] = { 9.927,
231  10.121,
232  10.421,
233  10.697,
234  10.948,
235  11.176,
236  11.383,
237  11.570,
238  11.737,
239  11.887,
240  12.017,
241  12.121,
242  12.130,
243  12.197};
244  static double innerCoreRho[7] = { 12.229,
245  12.301,
246  12.360,
247  12.405,
248  12.437,
249  12.455,
250  12.460};
251  // ============ END DENSITY TABLE ============
252 
253  if (r>fRearth) return 0.0;
254  double d = fRearth-r; // Table uses depth
255  int i;
256 
257  if (d>=crustD[0] && d<crustD[1]) {
258  for (i=0; i<1; ++i) {
259  if (d>=crustD[i] && d<crustD[i+1]) {
260  return crustRho[i];
261  }
262  }
263  }
264  if (d>=upperMantleD[0] && d<upperMantleD[5]) {
265  for (i=0; i<5; ++i) {
266  if (d>=upperMantleD[i] && d<upperMantleD[i+1]) {
267  return upperMantleRho[i];
268  }
269  }
270  }
271  if (d>=middleMantleD[0] && d<middleMantleD[5]) {
272  for (i=0; i<5; ++i) {
273  if (d>=middleMantleD[i] && d<middleMantleD[i+1]) {
274  return middleMantleRho[i];
275  }
276  }
277  }
278  if (d>=lowerMantleD[0] && d<lowerMantleD[13]) {
279  for (i=0; i<13; ++i) {
280  if (d>=lowerMantleD[i] && d<lowerMantleD[i+1]) {
281  return lowerMantleRho[i];
282  }
283  }
284  }
285  if (d>=outerCoreD[0] && d<outerCoreD[13]) {
286  for (i=0; i<13; ++i) {
287  if (d>=outerCoreD[i] && d<outerCoreD[i+1]) {
288  return outerCoreRho[i];
289  }
290  }
291  }
292  if (d>=innerCoreD[0] && d<innerCoreD[7]) {
293  for (i=0; i<7; ++i) {
294  if (d>=innerCoreD[i] && d<innerCoreD[i+1]) {
295  return innerCoreRho[i];
296  }
297  }
298  }
299  return 0.0;
300 }
301 
302 //......................................................................
303 ///
304 /// Mean ratio of atomic number to atomic mass for different segments
305 /// of the Earth. These are taken from:
306 ///
307 /// [*] Bahcall and Krastev, Phys. Rev. C56, p.2839 (1997)
308 ///
309 /// \param r - radius in km
310 /// \returns Z/A at this radius
311 ///
312 double EarthModel::ZoverA(double r)
313 {
314  if (r<fRouterCore) return 0.468;
315  return 0.497;
316 }
317 
318 //......................................................................
319 
320 double EarthModel::AveNe(double r1, double r2, int nsample)
321 {
322  double sum = 0.0;
323  double r;
324  for (int i=0; i<nsample; ++i) {
325  r = r1 + (r2-r1)*((float)i+0.5)/(float)nsample;
326  sum += this->Ne(r);
327  }
328  return sum/(float)nsample;
329 }
330 
331 //......................................................................
332 ///
333 /// Make layers of constant density keeping the electron density
334 /// profile to within a fractional tolerance
335 ///
336 void EarthModel::MakeLayers(double tol)
337 {
338  unsigned int nsample = 20;
339  for (unsigned int i=0; i<fRregion.size()-1; ++i) {
340  // Add layers within each region until all layers are within
341  // tolerance
342  double rRegionLo = fRregion[i];
343  double rRegionHi = fRregion[i+1];
344  double r1 = rRegionLo;
345  double r2 = rRegionHi;
346  double r;
347  double ne;
348  while (1) {
349  double ave = this->AveNe(r1, r2, nsample);
350  bool isok = true;
351  for (unsigned int j=0; j<nsample; ++j) {
352  r = r1+(r2-r1)*((float)j+0.5)/(float)nsample;
353  ne = this->Ne(r);
354  if (fabs(ne-ave)/ne>tol) { isok = false; break; }
355  }
356  if (isok) {
357  // Layer is within tolerance - accept it
358  fRlo.push_back(r1);
359  fRhi.push_back(r2);
360  fNe. push_back(ave);
361  if (r2>=rRegionHi) {
362  // Finished subdividing this region
363  break;
364  }
365  else {
366  // Ready next iteration of subdivision
367  r1 = r2;
368  r2 = rRegionHi;
369  }
370  } // if (isok) ...
371  else {
372  // Layer is not within tolerance - reduce its size
373  r2 -= 1.0; // Step back some distance (km) and try again
374  if (r2<=r1) r2 = r1+0.5;
375  }
376  } // making layers within region
377  } // loop on regions
378 }
379 
380 //......................................................................
381 
382 void EarthModel::GetLayers(std::vector<double>& rlo,
383  std::vector<double>& rhi,
384  std::vector<double>& ne)
385 {
386  rlo.resize(fRlo.size());
387  rhi.resize(fRlo.size());
388  ne. resize(fRlo.size());
389  for (unsigned int i=0; i<fRlo.size(); ++i) {
390  rlo[i] = fRlo[i];
391  rhi[i] = fRhi[i];
392  ne[i] = fNe[i];
393  }
394 }
395 
396 //......................................................................
397 ///
398 /// Find the list of distances a neutrino travels through each layer
399 /// of the Earth model
400 ///
401 /// \param anu - Altitude of neutrino production (km above sea level)
402 /// \param cosQ - cosine of neutrnio zenith angle (+1=down, -1=up)
403 /// \param adet - Altitude of detector (km above sea level)
404 ///
405 void EarthModel::LineProfile(double anu,
406  double cosQ,
407  double adet,
408  std::list<double>& Ls,
409  std::list<double>& Ns)
410 {
411  // A typical number density for the crust
412  static const double crustNe = 2.6*0.5;
413 
414  // Locations are relative to a cicular Earth centered at (0,0)
415  double rdet; // Radial location of detector (km)
416  double rnu; // Radial location of neutrino (km)
417  double Lnu; // Total flight distance of neutrino (km)
418  double x1; // Location of detector (km)
419  double y1; // Location of detector (km)
420  double x2; // Location of neutrino production (km)
421  double y2; // Location of neutrino production (km)
422 
423  rdet = fRearth+adet;
424  rnu = fRearth+anu;
425  Lnu = rdet*(sqrt(util::sqr(rnu/rdet)+cosQ*cosQ-1.0)-cosQ);
426 
427  x2 = 0.0;
428  y2 = rdet;
429  x1 = x2+sqrt(1.0-cosQ*cosQ)*Lnu;
430  y1 = y2+cosQ*Lnu;
431 
432  unsigned int i;
433  int n1, n2;
434  double L;
435  double xina, yina;
436  double xouta, youta;
437  double xinb, yinb;
438  double xoutb, youtb;
439  Ls.clear();
440  Ns.clear();
441 
442  //
443  // For down-going neutrinos, only consider matter in the top layer
444  //
445  if (cosQ>=0.0) {
446  if (adet>0.0) {
447  Ls.push_back(Lnu);
448  Ns.push_back(0.0);
449  }
450  else {
451  n1 = this->IntersectLineAndCircle(x1, y1, x2, y2, fRearth,
452  &xouta, &youta, &xoutb, &youtb);
453  L = util::pythag(x1-xoutb, y1-youtb);
454  Ls.push_back(L);
455  Ns.push_back(0.0);
456 
457  L = util::pythag(x2-xoutb, y2-youtb);
458  Ls.push_back(L);
459  Ns.push_back(crustNe); // Typical rock density
460  }
461  return;
462  }
463 
464  //
465  // For up-going neutrinos, find the intersections of the neutrino's
466  // path with all layers of the Earth model.
467  //
468  for (i=0; i<fRhi.size()-1; ++i) {
469  n1 = this->IntersectLineAndCircle(x1, y1, x2, y2, fRhi[i],
470  &xouta, &youta, &xoutb, &youtb);
471 
472  n2 = this->IntersectLineAndCircle(x1, y1, x2, y2, fRlo[i],
473  &xina, &yina, &xinb, &yinb);
474  //
475  // Create segments for each layer we intersect
476  //
477  if (n1==2 && n2<2) {
478  //
479  // Intersect outer radius, but not inner radius:
480  // Only one one segment to create.
481  //
482  L = util::pythag(xouta-xoutb, youta-youtb);
483 
484  Ls.push_back(L);
485  Ns.push_back(fNe[i]);
486  }
487  else if (n1==2 && n2==2) {
488  //
489  // Path intersects both inner and outer radii: Two sections to
490  // create, one on the way in, one on the way out.
491  //
492  L = util::pythag(xouta-xina, youta-yina);
493  Ls.push_front(L);
494  Ns.push_front(fNe[i]);
495 
496  L = util::pythag(xoutb-xinb, youtb-yinb);
497  Ls.push_back(L);
498  Ns.push_back(fNe[i]);
499  }
500  } // Loop on all but last layer
501 
502  //
503  // Handle the last layer specially
504  //
505  i = fRhi.size()-1;
506  n1 = this->IntersectLineAndCircle(x1, y1, x2, y2, fRhi[i],
507  &xouta, &youta, &xoutb, &youtb);
508  n2 = this->IntersectLineAndCircle(x1, y1, x2, y2, fRlo[i],
509  &xina, &yina, &xinb, &yinb);
510  if (n1==2 && n2<2) {
511  //
512  // Intersect outer radius, but not inner radius:
513  // One one segment to create
514  //
515  if (youtb>y2) {
516  //
517  // Don't over shoot the detector location
518  //
519  L = util::pythag(xouta-x2, youta-y2);
520  }
521  else {
522  L = util::pythag(xouta-xoutb, youta-youtb);
523  }
524  Ls.push_back(L);
525  Ns.push_back(fNe[i]);
526  }
527  else if (n1==2 && n2==2) {
528  //
529  // Path intersects both inner and outer radii: Two sections to
530  // create, one on the way in, one on the way out.
531  //
532  L = util::pythag(xouta-xina, youta-yina);
533  Ls.push_front(L);
534  Ns.push_front(fNe[i]);
535 
536  if (youtb>y2) {
537  //
538  // Don't over shoot the detector location
539  //
540  L = util::pythag(x2-xinb, y2-yinb);
541  }
542  else {
543  L = util::pythag(xoutb-xinb, youtb-yinb);
544  }
545  Ls.push_back(L);
546  Ns.push_back(fNe[i]);
547  }
548  //
549  // Add segment to get neutrino from production to Earth's surface
550  //
551  L = util::pythag(xouta-x1, youta-y1);
552  Ls.push_front(L);
553  Ns.push_front(0.0); // Density of air~=0
554 
555  //
556  // If required, add segment to get from sea level up to the
557  // detector's position
558  //
559  if (youtb<y2) {
560  L = util::pythag(xoutb-x2, youtb-y2);
561  Ls.push_back(L);
562  Ns.push_back(crustNe); // Density of of standard rock (mole/cm^3)
563  }
564 
565  //
566  // The following is a useful piece of debugging - check that we've
567  // accounted for all pieces of the neutrino flight distance
568  //
569  if (false) {
570  std::list<double>::iterator itr(Ls.begin());
571  std::list<double>::iterator itrEnd(Ls.end());
572  double ltot = 0.0;
573  for (; itr!=itrEnd; ++itr) {
574  ltot += (*itr);
575  }
576  std::cout << ltot << " " << Lnu << std::endl;
577  if (fabs(ltot-Lnu)>1e-6) abort();
578  }
579 }
580 
581 //......................................................................
582 
584  double x2, double y2,
585  double r,
586  double* xa, double* ya,
587  double* xb, double* yb)
588 {
589  double dx = x2-x1;
590  double dy = y2-y1;
591  double drsqr = dx*dx + dy*dy;
592  double D = x1*y2 - x2*y1;
593  double arg = r*r*drsqr - D*D;
594  if (arg<0.0) return 0;
595 
596  double sgndy = 1.0;
597  if (dy<0.0) sgndy = -1.0;
598 
599  *xa = ( D*dy - sgndy * dx * sqrt(arg) ) / drsqr;
600  *ya = (-D*dx - fabs(dy) * sqrt(arg) ) / drsqr;
601  *xb = ( D*dy + sgndy * dx * sqrt(arg) ) / drsqr;
602  *yb = (-D*dx + fabs(dy) * sqrt(arg) ) / drsqr;
603 
604  if (arg==0.0) return 1;
605  return 2;
606 }
607 
608 ////////////////////////////////////////////////////////////////////////
EarthModel(const char *which, double tol)
Definition: EarthModel.cxx:25
Filter events based on their run/event numbers.
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
Float_t y1[n_points_granero]
Definition: compare.C:5
Float_t x1[n_points_granero]
Definition: compare.C:5
T sqrt(T number)
Definition: d0nt_math.hpp:156
OStream cerr
Definition: OStream.cxx:7
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
void MakeLayers(double tol)
Definition: EarthModel.cxx:336
Provide and interface to different Earth models.
double DensityStacey(double r)
Definition: EarthModel.cxx:150
void GetLayers(std::vector< double > &rlo, std::vector< double > &rhi, std::vector< double > &ne)
Definition: EarthModel.cxx:382
int IntersectLineAndCircle(double x1, double y1, double x2, double y2, double r, double *xa, double *ya, double *xb, double *yb)
Definition: EarthModel.cxx:583
double dy[NP][NC]
double dx[NP][NC]
static constexpr double L
double Ne(double r)
Definition: EarthModel.cxx:46
void resize(T &x, std::vector< int > dims)
Definition: resize.hpp:41
base_types push_back(int_type())
T sqr(T x)
void LineProfile(double prodL, double cosQ, double rdet, std::list< double > &Ls, std::list< double > &Ns)
Definition: EarthModel.cxx:405
Float_t d
Definition: plot.C:236
const double j
Definition: BetheBloch.cxx:29
double AveNe(double r1, double r2, int nsample)
Definition: EarthModel.cxx:320
Oscillation probability calculators.
Definition: Calcs.h:5
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
double DensityPREM(double r)
Definition: EarthModel.cxx:101
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
double ZoverA(double r)
Definition: EarthModel.cxx:312
TRandom3 r(0)
double Density(double r)
Definition: EarthModel.cxx:55
static const int kPREM
Definition: EarthModel.cxx:15
double T
Definition: Xdiff_gwt.C:5
const cmplx< T > & y
Double_t sum
Definition: plot.C:31
Float_t e
Definition: plot.C:35
static const int kSTACEY
Definition: EarthModel.cxx:16
double Density(double r)
Definition: PREM.cxx:26