//____________________________________________________________________________
/*!

\class   genie::GPEarth

\brief   A class to evaluate the probability to transmission through the Earth
         The column depth is computed according to the Earth PREM model
         A. M. Dziewonski and D.L. Anderson, Physcis of the Earth and Planetary
         Interiors, 25 (1981) 297, Table 1

\author  Carla Distefano <distefano_c@lns.infn.it>
         LNS-INFN, Catania

\created October 17, 2012


\cpright  Copyright (c) 2012-2019, The KM3NeT Collaboration
          For the full text of the license see $GSEAGEN/LICENSE
*/
//____________________________________________________________________________


#ifndef _PEARTH_H_
#define _PEARTH_H_

#include <string>
#include <map>
#include <vector>
#include <cstdio>
#include <iostream>
#include <fstream>
#include <cmath>




#include "SeaEvent/GBinParam.h"

#define NLAYERS 10

class TGraph;

using namespace std;

struct InterPoint{int iLayer; double T;
                  bool operator < (const InterPoint& point) const {return (T < point.T);} };

struct ExtremPoint{int iLayer; double T1; double T2;};


struct ColumnDepthLayer{int iLayer; double ColumnDepth;};

struct Layers{double r1, r2, rho; string Composition;
                  bool operator < (const Layers& layer) const {return (r1 < layer.r1);}};

namespace genie {

class GPEarth {

public :

   GPEarth(GenParam * GenPar);
  ~GPEarth();

  double GetPEarth              () { return fPEarth;}
  double GetColumnDepth         () { return fColumnDepth;}
  double GetSigmaMean           () { return fSigmaMean;}
  double EarthAbsorption        (double X0, double Y0, double Z0, double l, double m, double n, double Energy, int neutrino_pdg);
  double GetREarth              () { return fREarth;}           // Get the total Earth radius (including ocean)
  double GetSeaBottomRadius (void) { return fRSeaBottom;}       // Get the Earth radius at the sea-bottom level

protected:

  GenParam *       fGenPar;

  double           fR[NLAYERS];                 ///< layer radii
  double           fEarthCoeff[NLAYERS][4];     ///< layer density coeff.
  string           fComp[NLAYERS];              ///< layer composition

  double           fREarth;             ///< Earth Radius
  double           fRSeaBottom;

  double           fPEarth;             ///< probability to transmission through the earth
  double           fColumnDepth;        ///< total column depth in mew
  double           fSigmaMean;          ///< mean cross section: fPEart=exp(-Na*fSigmaMean*fColumnDepth)


  double           fXc,fYc,fZc;         ///< coordinates of the Earth centre
  double           fX0,fY0,fZ0;         ///< coordinates of the neutrino vertex
  double           fParL, fParM, fParN; ///< X=X0+fParL*T; Y=Y0+fParM*T; Z=Z0+fParN*T is the parametric equation of straight line coincident with the neutrino track

  double           fA,fB,fC;

  vector<InterPoint>  fVecInterPoint;   ///< intersection points of the straight line with the layers (circonferences)
  vector<ExtremPoint> fVecExtremPoint;  ///< integration extremes for layers

  vector<ColumnDepthLayer>   fColumnDepthLayer;   ///< ColumnDepth for each crossed layer

  double           fGetXsecTot          (double Energy, int neutrino_pdg, map<int,TGraph> XSecTot);

  void             fCalcColumnDepth     (void);
  void             fExtremes            (void);

  void             fInterCir            (void);

  double           fINTEG0              (double T);
  double           fINTEG1              (double T);
  double           fINTEG2              (double T);
  double           fINTEG3              (double T);

  double           fIntLayer            (int iLayer, double T1, double T2); // integrates the single layer

  void             fBuildPREMEarth      (void); // build geometry for nu propagation through the Earth

};

} // genie namespace

#endif // _PEARTH_H_