#ifndef __JPHYSICS__ANTARES__
#define __JPHYSICS__ANTARES__

#include <map>
#include <cmath>

#include "JPhysics/JPhysicsSupportkit.hh"
#include "JPhysics/JK40Rates.hh"


/**
 * \file
 *
 * Properties of %Antares PMT and deep-sea water.
 * \author mdejong
 */
namespace ANTARES {

  using JPHYSICS::JK40Rates;

  
  /**
   * Get K40 rates.
   *
   * \return            K40 rates [Hz]
   */
  inline const JK40Rates& getK40Rates()
  {
    static const JK40Rates rates_Hz(70e3, { 25.0 });

    return rates_Hz;
  }

  
  /**
   * Get ambient pressure.
   *
   * \return            pressure [Atm]
   */
  inline double getAmbientPressure()
  {
    return 240.0;    // ambient pressure [Atm]
  }
  

  /**
   * Get photo-cathode area of PMT.
   *
   * \return            photo-cathode area [m^2]
   */
  inline double getPhotocathodeArea()
  {
    return 440e-4;   // photo-cathode area [m^2]
  }


  /**
   * Get absorption length.
   *
   * \param  lambda     wavelength of light [nm]
   * \return            absorption length   [m]
   */
  inline double getAbsorptionLength(const double lambda)
  {
    typedef std::map<double, double>     map_t;

    static map_t zmap;

    if (zmap.empty()) {
      /**
       * values taken from:
       *
       * afs/in2p3.fr/home/throng/antares/src/km3/v3r6/src/hit-ini_optic.f
       */
      //static const double a0 = 0.971;
      static const double a0 = 1.0;
      
      zmap[620.0] = 0.0;
      zmap[610.0] = a0/0.2890;
      zmap[600.0] = a0/0.2440;
      zmap[590.0] = a0/0.1570;
      zmap[580.0] = a0/0.1080;
      zmap[570.0] = a0/0.0799;
      zmap[560.0] = a0/0.0708;
      zmap[550.0] = a0/0.0638;
      zmap[540.0] = a0/0.0558;
      zmap[530.0] = a0/0.0507;
      zmap[520.0] = a0/0.0477;
      zmap[510.0] = a0/0.0357;
      zmap[500.0] = a0/0.0257;
      zmap[490.0] = a0/0.0196;
      zmap[480.0] = a0/0.0182;
      zmap[470.0] = a0/0.0182;
      zmap[460.0] = a0/0.0191;
      zmap[450.0] = a0/0.0200;
      zmap[440.0] = a0/0.0218;
      zmap[430.0] = a0/0.0237;
      zmap[420.0] = a0/0.0255;
      zmap[410.0] = a0/0.0291;
      zmap[400.0] = a0/0.0325;
      zmap[390.0] = a0/0.0363;
      zmap[380.0] = a0/0.0415;
      zmap[370.0] = a0/0.0473;
      zmap[360.0] = a0/0.0528;
      zmap[350.0] = a0/0.0629;
      zmap[340.0] = a0/0.0710;
      zmap[330.0] = a0/0.0792;
      zmap[320.0] = a0/0.0946;
      zmap[310.0] = a0/0.1090;
      zmap[300.0] = a0/0.1390;
      zmap[290.0] = 0.0;
    }

    const double x = lambda;

    if (x > zmap.begin()->first && x < zmap.rbegin()->first) {

      map_t::const_iterator i = zmap.lower_bound(x);
      map_t::const_iterator j = i;

      --j;

      if (i == zmap.begin()) {
	++i;
	++j;
      }	

      const double x1 = i->first;
      const double x2 = j->first;
    
      const double y1 = i->second;
      const double y2 = j->second;

      return y1  +  (y2 - y1) * (x - x1) / (x2 - x1);

    } else

      return 0.0;
  }


  /**
   * Get scattering length.
   *
   * \param  lambda     wavelength of light [nm]
   * \return            scattering length   [m]
   */
  inline double getScatteringLength(const double lambda)
  {
    typedef std::map<double, double>     map_t;

    static map_t zmap;

    if (zmap.empty()) {

      zmap[200.0] =    0.0;
      zmap[307.0] =   19.7132;
      zmap[330.0] =   23.8280;
      zmap[350.0] =   27.6075;
      zmap[370.0] =   31.5448;
      zmap[390.0] =   35.6151;
      zmap[410.0] =   39.7971;
      zmap[430.0] =   44.0726;
      zmap[450.0] =   48.4261;
      zmap[470.0] =   52.8447;
      zmap[490.0] =   57.3172;
      zmap[510.0] =   61.8344;
      zmap[530.0] =   66.3884;
      zmap[550.0] =   70.9723;
      zmap[570.0] =   75.5804;
      zmap[590.0] =   80.2077;
      zmap[610.0] =   84.8497;
      zmap[650.0] =   95.0;
    }

    const double x = lambda;

    if (x > zmap.begin()->first && x < zmap.rbegin()->first) {

      map_t::const_iterator i = zmap.lower_bound(x);
      map_t::const_iterator j = i;

      --j;

      if (i == zmap.begin()) {
	++i;
	++j;
      }	

      const double x1 = i->first;
      const double x2 = j->first;
    
      const double y1 = i->second;
      const double y2 = j->second;

      return y1  +  (y2 - y1) * (x - x1) / (x2 - x1);

    } else

      return 0.0;
  }

  
  /**
   * Function to describe light scattering in water.
   *
   * \param  x      cosine scattering angle
   * \return        probability
   */
  inline double getScatteringProbability(const double x)
  {
    return JPHYSICS::p00075(x);
  }


  /**
   * Angular acceptance of PMT (Genova).
   *
   * \param  x      cosine of angle of incidence
   * \return        probability
   */
  inline double genova(const double x)
  {
    static const double a0 =  0.3265;
    static const double a1 =  0.6144;
    static const double a2 = -0.0343;
    static const double a3 = -0.0641;
    static const double a4 =  0.2988;
    static const double a5 = -0.1422;
    
    const double z = -x;

    double y = 0.0;

    if (z < -0.65)
      y = 0.0;
    else
      y = a0 + z*(a1 + z*(a2 + z*(a3 + z*(a4 + z*a5))));

    return y;
  }


  /**
   * Angular acceptance of PMT (Gamelle).
   *
   * \param  x      cosine angle of incidence
   * \return        probability
   */
  inline double gamelle(const double x)
  {    
    static const double a0 = 59.115;
    static const double a1 =  0.52258;
    static const double a2 =  0.60944E-02;
    static const double a3 = -0.16955E-03;
    static const double a4 =  0.60929E-06;

    double y = 0.0;

    if      (x > +0.36)
      y = 0.0;
    else if (x < -1.00)
      y = 1.0;
    else {

      const double z = acos(-x)*57.29578 + 57.75;

      y = (a0 + z*(a1 + z*(a2 + z*(a3 + z*a4)))) / 84.0;
    }

    return y;
  }


  /**
   * Get angular acceptance of PMT.
   *
   * \param  x      cosine of angle of incidence
   * \return        probability
   */
  inline double getAngularAcceptance(const double x)
  {
    return genova(x);
  }


  /**
   * Get quantum efficiency of PMT.
   *
   * \param   lambda  wavelength of photon [nm]
   * \param   option  include absorption in glass and gel (if true, otherwise not)
   * \return          quantum efficiency
   */
  inline double getQE(const double lambda, const bool option)
  {
    class tuple {
    public:
      tuple(const double __QE,
	    const double __l_gel,
	    const double __l_glass) :
	QE     (__QE),
	l_gel  (__l_gel),
	l_glass(__l_glass)
      {}

      double QE;         // Quantum efficiiency
      double l_gel;      // gel   absorption length [cm]
      double l_glass;    // glass absorption length [cm]
    };

    static const tuple ntuple[] = {
      tuple(0.000e-2,   0.00,   0.00),
      tuple(1.988e-2, 100.81, 148.37),
      tuple(2.714e-2,  99.94, 142.87),
      tuple(3.496e-2,  99.89, 135.64),
      tuple(4.347e-2,  96.90, 134.58),
      tuple(5.166e-2,  96.42, 138.27),
      tuple(6.004e-2,  94.36, 142.40),
      tuple(6.885e-2,  89.09, 147.16),
      tuple(8.105e-2,  90.10, 151.80),
      tuple(10.13e-2,  86.95, 150.88),
      tuple(13.03e-2,  85.88, 145.68),
      tuple(15.29e-2,  84.49, 139.70),
      tuple(16.37e-2,  81.08, 126.55),
      tuple(17.11e-2,  78.18, 118.86),
      tuple(17.86e-2,  76.48, 113.90),
      tuple(18.95e-2,  74.55, 116.08),
      tuple(20.22e-2,  72.31, 109.23),
      tuple(21.26e-2,  68.05,  81.63),
      tuple(22.10e-2,  66.91,  65.66),
      tuple(22.65e-2,  64.48,  77.30),
      tuple(23.07e-2,  62.53,  73.02),
      tuple(23.14e-2,  59.38,  81.25),
      tuple(23.34e-2,  56.64, 128.04),
      tuple(22.95e-2,  53.29,  61.84),
      tuple(22.95e-2,  48.96,  19.23),
      tuple(22.74e-2,  45.71,  27.21),
      tuple(23.48e-2,  41.88,  18.09),
      tuple(22.59e-2,  37.14,   8.41),
      tuple(20.61e-2,  30.49,   3.92),
      tuple(17.68e-2,  23.08,   1.82),
      tuple(13.18e-2,  15.60,   0.84),
      tuple(7.443e-2,   8.00,   0.39),
      tuple(2.526e-2,   0.00,   0.17),
      tuple(0.000e-2,   0.00,   0.00)
    };

    static const double cola    = 0.9;      // collection efficiency
    static const double x_glass = 1.5;      // glass thickness [cm]
    static const double x_gel   = 1.0;      // gel   thickness [cm]

    // ntuple 
    static const int    N    = sizeof(ntuple) / sizeof(ntuple[0])  -  1;

    static const double xmax = 620.0;       // maximal wavelength [nm] (tuple[ 0 ])
    static const double xmin = 290.0;       // minimal wavelength [nm] (tuple[N-1])

    const double x = lambda;

    double y = 0.0;

    if (x > xmin && x < xmax) {

      const int i = (int) (N * (x - xmax) / (xmin - xmax));
      const int j = (i == N ? i - 1 : i + 1); 
    
      const double x1 = xmax  +  i * (xmin - xmax) / N;
      const double x2 = xmax  +  j * (xmin - xmax) / N;

      const double dx = (x - x1) / (x1 - x2);

      const double QE      = ntuple[i].QE       +  (ntuple[i].QE      - ntuple[j].QE     ) * dx;
      const double l_gel   = ntuple[i].l_gel    +  (ntuple[i].l_gel   - ntuple[j].l_gel  ) * dx;
      const double l_glass = ntuple[i].l_glass  +  (ntuple[i].l_glass - ntuple[j].l_glass) * dx;

      y = cola * QE;

      if (option) {

	if (l_glass > 0.0 && l_gel > 0.0)
	  y *= exp(-x_glass/l_glass) * exp(-x_gel/l_gel);
	else
	  y = 0.0;
      }
    } 

    return y;
  }


  /**
   * Get quantum efficiency of PMT.
   *
   * \param   lambda  wavelength of photon [nm]
   * \return          quantum efficiency
   */
  inline double getQE(const double lambda)
  {
    return getQE(lambda, true);
  }


  /**
   * Get effective photo-cathode area of PMT.
   *
   * \param  x          cosine of angle of incidence
   * \param  lambda     wavelength of photon [nm]
   * \return            photo-cathode area [m^2]
   */
  inline double getPhotocathodeArea2D(const double x, const double lambda)
  {
    return getPhotocathodeArea() * getQE(lambda) * getAngularAcceptance(x);
  }
}

#endif