osc::EarthModel Class Reference

#include <EarthModel.h>

Public Member Functions
	EarthModel (const char *which, double tol)
double	Ne (double r)
double	Density (double r)
double	ZoverA (double r)
void	GetLayers (std::vector< double > &rlo, std::vector< double > &rhi, std::vector< double > &ne)
void	LineProfile (double prodL, double cosQ, double rdet, std::list< double > &Ls, std::list< double > &Ns)

Private Member Functions
void	InitPREM ()
void	InitStacey ()
double	DensityPREM (double r)
double	DensityStacey (double r)
double	AveNe (double r1, double r2, int nsample)
void	MakeLayers (double tol)
int	IntersectLineAndCircle (double x1, double y1, double x2, double y2, double r, double *xa, double *ya, double *xb, double *yb)

Private Attributes
int	fModelID
double	fRouterCore
double	fRearth
std::vector< double >	fRregion
std::vector< double >	fRlo
std::vector< double >	fRhi
std::vector< double >	fNe

Detailed Description

Definition at line 13 of file EarthModel.h.

Constructor & Destructor Documentation

EarthModel::EarthModel	(	const char *	which,
		double	tol
	)

Definition at line 18 of file EarthModel.cxx.

{
  if      (std::string("PREM")  ==which) this->InitPREM();
  else if (std::string("STACEY")==which) this->InitStacey();
  // else if (std::string("ModelX")==which) this->InitModelX();
  else {
    std::cerr << __FILE__ << ":" << __LINE__ 
              << " Model '" << which << "' is not supported." 
              << " Currently only PREM is supported." << std::endl;
    abort();
  }
  this->MakeLayers(tol);
}

Member Function Documentation

double EarthModel::AveNe ( double  r1, double  r2, int  nsample  )

private

Definition at line 266 of file EarthModel.cxx.

{
  double sum = 0.0;
  double r;
  for (int i=0; i<nsample; ++i) {
    r = r1 + (r2-r1)*((float)i+0.5)/(float)nsample;
    sum += this->Ne(r);
  }
  return sum/(float)nsample;
}

double EarthModel::Density

(

double

r

)

Definition at line 39 of file EarthModel.cxx.

{
  switch (fModelID) {
  case kPREM:   return this->DensityPREM(r);
  case kSTACEY: return this->DensityStacey(r);
  default:      abort();
  }
  return 0.0;
}

double EarthModel::DensityPREM ( double  r )

private

Definition at line 69 of file EarthModel.cxx.

{
  double x = r/fRearth;
  
  // The regions are checked in order from largest to smallest
  if (r>1221.5&&r<=3480.0) return 12.5815+x*(-1.2638+x*(-3.6426+x*(-5.5281)));
  if (r>3480.0&&r<=5701.0) return  7.9565+x*(-6.4761+x*(5.5283+x*(-3.0807)));
  if (r<=1221.5)           return 13.0885-8.8381*x*x;
  if (r>5771.0&&r<=5971.0) return 11.2494-8.0298*x;
  if (r>5971.0&&r<=6151.0) return  7.1089-3.8045*x;
  if (r>6151.0&&r<=6291.0) return  2.6910+0.6924*x;
  if (r>5701.0&&r<=5771.0) return  5.3197-1.4836*x;
  if (r>6291.0&&r<=6346.6) return  2.6910+0.6924*x;
  if (r>6346.6&&r<=6356.0) return  2.90;
  if (r>6356.0&&r<=6368.0) return  2.60;
  if (r>6368.0&&r<=6371.0) return  1.020;
  return 0.0;
}

double EarthModel::DensityStacey ( double  r )

private

Definition at line 105 of file EarthModel.cxx.

{
  static double crustD[2] = {       0.0,
                                    15.0};
  static double upperMantleD[6] = { 15.0,
                                    60.0,
                                    100.0,
                                    200.0,
                                    300.0, 
                                    350.0};
  static double middleMantleD[6] = {350.0,
                                    400.0,
                                    413.0,
                                    500.0,
                                    600.0, 
                                    650.0};
  static double lowerMantleD[14] = {650.0,
                                    800.0,
                                    984.0,
                                    1000.0,
                                    1200.0, 
                                    1400.0,
                                    1600.0,
                                    1800.0,
                                    2000.0,
                                    2200.0, 
                                    2400.0,
                                    2600.0,
                                    2800.0,
                                    2878.0};
  static double outerCoreD[14] = {  2878.0,
                                    3000.0,
                                    3200.0,
                                    3400.0,
                                    3600.0, 
                                    3800.0,
                                    4000.0,
                                    4200.0,
                                    4400.0,
                                    4600.0, 
                                    4800.0,
                                    4982.0,
                                    5000.0,
                                    5121.0};
  static double innerCoreD[8] ={    5121.0,
                                    5200.0,
                                    5400.0,
                                    5600.0,
                                    5800.0, 
                                    6000.0,
                                    6200.0,
                                    6371.0};
  static double crustRho[2] =       { 2.840,
                                      2.840};
  static double upperMantleRho[6] = { 3.313,
                                      3.332,
                                      3.348,
                                      3.387,
                                      3.424,
                                      3.441};
  static double middleMantleRho[6] = {3.700,
                                      3.775,
                                      3.795,
                                      3.925,
                                      4.075, 
                                      4.150};
  static double lowerMantleRho[14] = {4.200,
                                      4.380,
                                      4.529,
                                      4.538,
                                      4.655,
                                      4.768,
                                      4.877,
                                      4.983,
                                      5.087,
                                      5.188,
                                      5.288,
                                      5.387,
                                      5.487,
                                      5.527};
  static double outerCoreRho[14] = {  9.927,
                                      10.121,
                                      10.421,
                                      10.697,
                                      10.948, 
                                      11.176,
                                      11.383,
                                      11.570,
                                      11.737,
                                      11.887, 
                                      12.017,
                                      12.121,
                                      12.130,
                                      12.197};
  static double innerCoreRho[7] = {   12.229,
                                      12.301,
                                      12.360,
                                      12.405,
                                      12.437, 
                                      12.455,
                                      12.460};
  // ============ END DENSITY TABLE ============
  
  if (r>fRearth) return 0.0;
  double       d = fRearth-r; // Table uses depth
  register int i;
  
  if (d>=crustD[0] && d<crustD[1]) {
    for (i=0; i<1; ++i) {
      if (d>=crustD[i] && d<crustD[i+1]) {
        return crustRho[i];
      }
    }
  }
  if (d>=upperMantleD[0] && d<upperMantleD[5]) {
    for (i=0; i<5; ++i) {
      if (d>=upperMantleD[i] && d<upperMantleD[i+1]) {
        return upperMantleRho[i];
      }
    }
  }
  if (d>=middleMantleD[0] && d<middleMantleD[5]) {
    for (i=0; i<5; ++i) {
      if (d>=middleMantleD[i] && d<middleMantleD[i+1]) {
        return middleMantleRho[i];
      }
    }
  }
  if (d>=lowerMantleD[0] && d<lowerMantleD[13]) {
    for (i=0; i<13; ++i) {
      if (d>=lowerMantleD[i] && d<lowerMantleD[i+1]) {
        return lowerMantleRho[i];
      }
    }
  }
  if (d>=outerCoreD[0] && d<outerCoreD[13]) {
    for (i=0; i<13; ++i) {
      if (d>=outerCoreD[i] && d<outerCoreD[i+1]) {
        return outerCoreRho[i];
      }
    }
  }
  if (d>=innerCoreD[0] && d<innerCoreD[7]) {
    for (i=0; i<7; ++i) {
      if (d>=innerCoreD[i] && d<innerCoreD[i+1]) {
        return innerCoreRho[i];
      }
    }
  }
  return 0.0;
}

void EarthModel::GetLayers	(	std::vector< double > &	rlo,
		std::vector< double > &	rhi,
		std::vector< double > &	ne
	)

Definition at line 324 of file EarthModel.cxx.

{
  rlo.resize(fRlo.size());
  rhi.resize(fRlo.size());
  ne. resize(fRlo.size());
  for (unsigned int i=0; i<fRlo.size(); ++i) {
    rlo[i] = fRlo[i];
    rhi[i] = fRhi[i];
    ne[i]  = fNe[i];
  }
}

void EarthModel::InitPREM ( )

private

Definition at line 50 of file EarthModel.cxx.

{
  fModelID      = kPREM;
  fRouterCore   = 3480.0;
  fRearth       = 6371.0;

  fRregion.push_back(0.0);
  fRregion.push_back(1221.5); // Inner core
  fRregion.push_back(3480.0); // Outer core
  fRregion.push_back(3630.0); // D''
  fRregion.push_back(5701.0); // Lower mantle
  fRregion.push_back(6151.0); // Transition zone
  fRregion.push_back(6291.0); // Low velocity zone
  fRregion.push_back(6346.6); // LID
  fRregion.push_back(6368.0); // Crust
  fRregion.push_back(6371.0); // Ocean
}

void EarthModel::InitStacey ( )

private

Definition at line 89 of file EarthModel.cxx.

{
  fModelID      = kSTACEY;
  fRouterCore   = 6371.0-2878.0;
  fRearth       = 6371.0;

  fRregion.push_back(0.0);
  fRregion.push_back(6371.0-5121.0); // Inner core
  fRregion.push_back(6371.0-2878.0); // Outer core
  fRregion.push_back(6371.0-650.0);  // Lower mantle
  fRregion.push_back(6371.0-350.0);  // Middle mantle
  fRregion.push_back(6371.0-15.0);   // Upper mantle
  fRregion.push_back(6371.0);        // Crust
}

int EarthModel::IntersectLineAndCircle ( double  x1, double  y1, double  x2, double  y2, double  r, double *  xa, double *  ya, double *  xb, double *  yb  )

private

Definition at line 517 of file EarthModel.cxx.

{
  double dx    = x2-x1;
  double dy    = y2-y1;
  double drsqr = dx*dx + dy*dy;
  double D     = x1*y2 - x2*y1;
  double arg   = r*r*drsqr - D*D;
  if (arg<0.0) return 0;

  double sgndy = 1.0;
  if (dy<0.0) sgndy = -1.0;

  *xa = ( D*dy - sgndy * dx * sqrt(arg) ) / drsqr;
  *ya = (-D*dx - fabs(dy)   * sqrt(arg) ) / drsqr;
  *xb = ( D*dy + sgndy * dx * sqrt(arg) ) / drsqr;
  *yb = (-D*dx + fabs(dy)   * sqrt(arg) ) / drsqr;

  if (arg==0.0) return 1;
  return 2;
}

void EarthModel::LineProfile	(	double	prodL,
		double	cosQ,
		double	rdet,
		std::list< double > &	Ls,
		std::list< double > &	Ns
	)

Definition at line 339 of file EarthModel.cxx.

{
  // A typical number density for the crust
  static const double crustNe = 2.6*0.5;
  
  // Locations are relative to a cicular Earth centered at (0,0)
  double rdet; // Radial location of detector (km)
  double rnu;  // Radial location of neutrino (km)
  double Lnu;  // Total flight distance of neutrino (km)
  double x1;   // Location of detector (km)
  double y1;   // Location of detector (km)
  double x2;   // Location of neutrino production (km)
  double y2;   // Location of neutrino production (km)
  
  rdet = fRearth+adet;
  rnu  = fRearth+anu;
  Lnu  = rdet*(sqrt(util::sqr(rnu/rdet)+cosQ*cosQ-1.0)-cosQ);

  x2 = 0.0;
  y2 = rdet;
  x1 = x2+sqrt(1.0-cosQ*cosQ)*Lnu;
  y1 = y2+cosQ*Lnu;
  
  unsigned int i;
  int n1, n2;
  double L;
  double xina,  yina;
  double xouta, youta;
  double xinb,  yinb;
  double xoutb, youtb;
  Ls.clear();
  Ns.clear();

  //
  // For down-going neutrinos, only consider matter in the top layer
  //
  if (cosQ>=0.0) {
    if (adet>0.0) {
      Ls.push_back(Lnu);
      Ns.push_back(0.0);
    }
    else {
      n1 = this->IntersectLineAndCircle(x1, y1, x2, y2, fRearth,
                                        &xouta, &youta, &xoutb, &youtb);
      L = util::pythag(x1-xoutb, y1-youtb);
      Ls.push_back(L);
      Ns.push_back(0.0);
      
      L = util::pythag(x2-xoutb, y2-youtb);
      Ls.push_back(L);
      Ns.push_back(crustNe); // Typical rock density
    }
    return;
  }
  
  //
  // For up-going neutrinos, find the intersections of the neutrino's
  // path with all layers of the Earth model.
  //
  for (i=0; i<fRhi.size()-1; ++i) {
    n1 = this->IntersectLineAndCircle(x1, y1, x2, y2, fRhi[i],
                                      &xouta, &youta, &xoutb, &youtb);
    
    n2 = this->IntersectLineAndCircle(x1, y1, x2, y2, fRlo[i],
                                      &xina, &yina, &xinb, &yinb);
    //
    // Create segments for each layer we intersect
    //
    if (n1==2 && n2<2) {
      //
      // Intersect outer radius, but not inner radius: 
      // Only one one segment to create.
      //
      L = util::pythag(xouta-xoutb, youta-youtb);
      
      Ls.push_back(L);
      Ns.push_back(fNe[i]);
    }
    else if (n1==2 && n2==2) {
      //
      // Path intersects both inner and outer radii: Two sections to
      // create, one on the way in, one on the way out.
      //
      L = util::pythag(xouta-xina, youta-yina);
      Ls.push_front(L);
      Ns.push_front(fNe[i]);
      
      L = util::pythag(xoutb-xinb, youtb-yinb);
      Ls.push_back(L);
      Ns.push_back(fNe[i]);
    }
  } // Loop on all but last layer
  
  //
  // Handle the last layer specially
  //
  i = fRhi.size()-1;
  n1 = this->IntersectLineAndCircle(x1, y1, x2, y2, fRhi[i],
                                    &xouta, &youta, &xoutb, &youtb);
  n2 = this->IntersectLineAndCircle(x1, y1, x2, y2, fRlo[i],
                                    &xina, &yina, &xinb, &yinb);
  if (n1==2 && n2<2) {
    //
    // Intersect outer radius, but not inner radius: 
    // One one segment to create
    //
    if (youtb>y2) {
      //
      // Don't over shoot the detector location
      //
      L = util::pythag(xouta-x2, youta-y2);
    }
    else {
      L = util::pythag(xouta-xoutb, youta-youtb);
    }
    Ls.push_back(L);
    Ns.push_back(fNe[i]);
  }
  else if (n1==2 && n2==2) {
    //
    // Path intersects both inner and outer radii: Two sections to
    // create, one on the way in, one on the way out.
    //
    L = util::pythag(xouta-xina, youta-yina);
    Ls.push_front(L);
    Ns.push_front(fNe[i]);
    
    if (youtb>y2) {
      //
      // Don't over shoot the detector location
      //
      L = util::pythag(x2-xinb, y2-yinb);
    }
    else {
      L = util::pythag(xoutb-xinb, youtb-yinb);
    }
    Ls.push_back(L);
    Ns.push_back(fNe[i]);
  }  
  //
  // Add segment to get neutrino from production to Earth's surface
  //
  L = util::pythag(xouta-x1, youta-y1);
  Ls.push_front(L);
  Ns.push_front(0.0); // Density of air~=0
  
  //
  // If required, add segment to get from sea level up to the
  // detector's position
  //
  if (youtb<y2) {
    L = util::pythag(xoutb-x2, youtb-y2);
    Ls.push_back(L);
    Ns.push_back(crustNe); // Density of of standard rock (mole/cm^3)
  }

  //
  // The following is a useful piece of debugging - check that we've
  // accounted for all pieces of the neutrino flight distance
  //
  if (false) {
    std::list<double>::iterator itr(Ls.begin());
    std::list<double>::iterator itrEnd(Ls.end());
    double ltot = 0.0;
    for (; itr!=itrEnd; ++itr) {
      ltot += (*itr);
    }
    std::cout << ltot << " " << Lnu << std::endl;
    if (fabs(ltot-Lnu)>1e-6) abort();
  }
}

void EarthModel::MakeLayers ( double  tol )

private

Definition at line 278 of file EarthModel.cxx.

{
  unsigned int nsample = 20;
  for (unsigned int i=0; i<fRregion.size()-1; ++i) {
    // Add layers within each region until all layers are within
    // tolerance
    double rRegionLo = fRregion[i];
    double rRegionHi = fRregion[i+1];
    double r1 = rRegionLo;
    double r2 = rRegionHi;
    double r;
    double ne;
    while (1) {
      double ave  = this->AveNe(r1, r2, nsample);
      bool   isok = true;
      for (unsigned int j=0; j<nsample; ++j) {
        r = r1+(r2-r1)*((float)j+0.5)/(float)nsample;
        ne = this->Ne(r);
        if (fabs(ne-ave)/ne>tol) { isok = false; break; }
      }
      if (isok) {
        // Layer is within tolerance - accept it
        fRlo.push_back(r1);
        fRhi.push_back(r2);
        fNe. push_back(ave);
        if (r2>=rRegionHi) {
          // Finished subdividing this region
          break;
        }
        else {
          // Ready next iteration of subdivision
          r1 = r2;
          r2 = rRegionHi;
        }
      } // if (isok) ...
      else {
        // Layer is not within tolerance - reduce its size
        r2 -= 1.0; // Step back some distance (km) and try again
        if (r2<=r1) r2 = r1+0.5;
      }
    } // making layers within region
  } // loop on regions
}

double EarthModel::Ne

(

double

r

)

Definition at line 33 of file EarthModel.cxx.

{
  return this->ZoverA(r)*this->Density(r);
}

double EarthModel::ZoverA

(

double

r

)

Definition at line 258 of file EarthModel.cxx.

{
  if (r<fRouterCore) return 0.468;
  return 0.497;
}

Member Data Documentation

int osc::EarthModel::fModelID

private

Definition at line 44 of file EarthModel.h.

std::vector<double> osc::EarthModel::fNe

private

Definition at line 51 of file EarthModel.h.

double osc::EarthModel::fRearth

private

Definition at line 46 of file EarthModel.h.

std::vector<double> osc::EarthModel::fRhi

private

Definition at line 50 of file EarthModel.h.

std::vector<double> osc::EarthModel::fRlo

private

Definition at line 49 of file EarthModel.h.

double osc::EarthModel::fRouterCore

private

Definition at line 45 of file EarthModel.h.

std::vector<double> osc::EarthModel::fRregion

private

Definition at line 47 of file EarthModel.h.

---

The documentation for this class was generated from the following files:

• g4lbne/ProductionScripts/OscLib/EarthModel.h
• g4lbne/ProductionScripts/OscLib/EarthModel.cxx