#ifndef __JPHYSICS__KM3NET__
#define __JPHYSICS__KM3NET__

#include <map>
#include <cmath>

#include "JTools/JFunction1D_t.hh"

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


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

  using JTOOLS::JPolint1Function1D_t;
  using JTOOLS::JGridPolint1Function1D_t;
  using JPHYSICS::JK40Rates;

  
  /**
   * Get K40 rates.
   *
   * This method returns an object with default values.\n
   * The singles and multiples rates are based on Analysis e-log entry <a href="https://elog.km3net.de/Analysis/597">575</a>.\n
   * A dark count of 700 Hz has been included in the singles rate.\n
   * See also %KM3NeT internal <a href="https://simulation.pages.km3net.de/input_tables/Simulations_Description.pdf">note</a>.
   *
   * \return            K40 rates [Hz]
   */
  inline const JK40Rates& getK40Rates()
  {
    static const JK40Rates rates_Hz(5200, {  568.0,  49.10,   5.48,   0.48});

    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 45.4e-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)
  {
    static JPolint1Function1D_t f1;

    if (f1.empty()) {

      f1[290] =   0.0;
      f1[310] =  11.9;
      f1[330] =  16.4;
      f1[350] =  20.6;
      f1[375] =  29.5;
      f1[412] =  48.5;
      f1[440] =  67.5;
      f1[475] =  59.0;
      f1[488] =  55.1;
      f1[510] =  26.1;
      f1[532] =  19.9;
      f1[555] =  14.7;
      f1[650] =   2.8;
      f1[676] =   2.3;
      f1[715] =   1.0;
      f1[720] =   0.0;

      f1.compile();

      f1.setExceptionHandler(new JPolint1Function1D_t::JDefaultResult(0.0));
    }

    return f1(lambda);
  }


  /**
   * Get scattering length.
   *
   * Use the Kopelevich model for spectral volume scattering functions.\n
   * The model separates the contributions by:
   * - pure_sea:   pure sea water;
   * - small_par: 'small' particles (size < 1 micro meter);
   * - large_par: 'large' particles (size > 1 micro meter).
   * Values are taken from reference C.D. Mobley "Light and Water", ISBN 0-12-502750-8, pag. 119.
   *
   *
   * \param  lambda     wavelength of light [nm]
   * \return            scattering length   [m]
   */
  inline double getScatteringLength(const double lambda)
  {
    static const double Vs = 0.0075;
    static const double Vl = 0.0075;
    static const double bw = 0.0017;
    static const double bs = 1.340;
    static const double bl = 0.312;

    const double x = 550.0/lambda;

    const double pure_sea  = bw *      pow(x, 4.3);
    const double small_par = bs * Vs * pow(x, 1.7);
    const double large_par = bl * Vl * pow(x, 0.3);

    return 1.0 / (pure_sea + small_par + large_par);
  }

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


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

    if (f1.empty()) {

      f1[-1.00] = 1.621;
      f1[-0.95] = 1.346;
      f1[-0.90] = 1.193;
      f1[-0.85] = 1.073;
      f1[-0.80] = 0.973;
      f1[-0.75] = 0.877;
      f1[-0.70] = 0.790;
      f1[-0.65] = 0.711;
      f1[-0.60] = 0.640;
      f1[-0.55] = 0.575;
      f1[-0.50] = 0.517;
      f1[-0.45] = 0.450;
      f1[-0.40] = 0.396;
      f1[-0.35] = 0.341;
      f1[-0.30] = 0.295;
      f1[-0.25] = 0.249;
      f1[-0.20] = 0.207;
      f1[-0.15] = 0.166;
      f1[-0.10] = 0.128;
      f1[-0.05] = 0.095;
      f1[0.00] = 0.065;
      f1[0.05] = 0.038;
      f1[0.10] = 0.017;
      f1[0.15] = 0.006;
      f1[0.20] = 0.003;
      f1[0.25] = 0.002;
      f1[0.30] = 0.001;
      f1[0.35] = 0.001;
      f1[0.40] = 0.001;

      f1.compile();

      f1.setExceptionHandler(new JGridPolint1Function1D_t::JDefaultResult(0.0));
    }

    return f1(x);
  }


  /**
   * Get quantum efficiency of %KM3NeT PMT.
   *
   * Note that the QE includes absorption in glass and gel.
   *
   * \param   lambda  wavelength of photon [nm]
   * \return          quantum efficiency
   */
  inline double getQE(const double lambda)
  {
    static JGridPolint1Function1D_t f1;
    static const double cola    = 0.9;      // collection efficiency

    if (f1.empty()) {

      f1[270] = 0.0;
      f1[275] = 0.0;
      f1[280] = 0.0;
      f1[285] = 0.0;
      f1[290] = 0.0;
      f1[295] = 0.1;
      f1[300] = 0.1;
      f1[305] = 0.9;
      f1[310] = 1.8;
      f1[315] = 4.6;
      f1[320] = 7.4;
      f1[325] = 11.3;
      f1[330] = 15.2;
      f1[335] = 18.1;
      f1[340] = 21.0;
      f1[345] = 22.4;
      f1[350] = 23.7;
      f1[355] = 24.5;
      f1[360] = 25.3;
      f1[365] = 25.8;
      f1[370] = 26.3;
      f1[375] = 26.2;
      f1[380] = 26.0;
      f1[385] = 26.6;
      f1[390] = 27.1;
      f1[395] = 26.9;
      f1[400] = 26.8;
      f1[405] = 26.5;
      f1[410] = 26.3;
      f1[415] = 26.0;
      f1[420] = 25.7;
      f1[425] = 25.6;
      f1[430] = 25.4;
      f1[435] = 25.1;
      f1[440] = 24.9;
      f1[445] = 24.5;
      f1[450] = 24.1;
      f1[455] = 23.5;
      f1[460] = 22.9;
      f1[465] = 22.1;
      f1[470] = 21.3;
      f1[475] = 20.6;
      f1[480] = 19.9;
      f1[485] = 19.5;
      f1[490] = 19.1;
      f1[495] = 18.8;
      f1[500] = 18.5;
      f1[505] = 18.1;
      f1[510] = 17.7;
      f1[515] = 16.8;
      f1[520] = 15.9;
      f1[525] = 14.3;
      f1[530] = 12.7;
      f1[535] = 11.5;
      f1[540] = 10.2;
      f1[545] = 9.5;
      f1[550] = 8.7;
      f1[555] = 8.1;
      f1[560] = 7.6;
      f1[565] = 7.2;
      f1[570] = 6.8;
      f1[575] = 6.4;
      f1[580] = 6.0;
      f1[585] = 5.6;
      f1[590] = 5.1;
      f1[595] = 4.8;
      f1[600] = 4.4;
      f1[605] = 4.0;
      f1[610] = 3.6;
      f1[615] = 3.3;
      f1[620] = 3.0;
      f1[625] = 2.7;
      f1[630] = 2.3;
      f1[635] = 2.1;
      f1[640] = 1.9;
      f1[645] = 1.7;
      f1[650] = 1.5;
      f1[655] = 1.4;
      f1[660] = 1.2;
      f1[665] = 1.0;
      f1[670] = 0.9;
      f1[675] = 0.7;
      f1[680] = 0.6;
      f1[685] = 0.5;
      f1[690] = 0.4;
      f1[695] = 0.3;
      f1[700] = 0.2;
      f1[705] = 0.1;
      f1[710] = 0.0;

      f1.compile();

      f1.setExceptionHandler(new JGridPolint1Function1D_t::JDefaultResult(0.0));
    }

    return f1(lambda) * 0.01 * cola;
  }


  /**
   * 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);
  }
}

/**
 * Name space for %KM3NeT with high-QE PMTs.
 */
namespace KM3NET_HIGHQE {

  using JTOOLS::JGridPolint1Function1D_t;
  using JPHYSICS::JK40Rates;

  
  /**
   * Get K40 rates.
   *
   * This method returns an object with default values.\n
   * The singles and multiples rates are based on Analysis e-log entry <a href="https://elog.km3net.de/Analysis/597">575</a>.\n
   * They have been scaled with the outputs of JRateK40.cc and JMultiplicityK40.cc for the standard and high-QE PMT, respectively using script JPrintK40Rates.sh
   * (see GIT issue <a href="https://git.km3net.de/common/jpp/-/issues/351">351</a>).
   *
   * \return            K40 rates [Hz]
   */
  inline const JK40Rates& getK40Rates()
  {
    static const JK40Rates rates_Hz(8529, { 1014.4, 103.77,  15.44,   1.44});

    return rates_Hz;
  }

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

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


  /**
   * Get absorption length.
   *
   * \param  lambda     wavelength of light [nm]
   * \return            absorption length   [m]
   */
  inline double getAbsorptionLength(const double lambda)
  {
    return KM3NET::getAbsorptionLength(lambda);
  }


  /**
   * Get scattering length.
   *
   * \param  lambda     wavelength of light [nm]
   * \return            scattering length   [m]
   */
  inline double getScatteringLength(const double lambda)
  {
    return KM3NET::getScatteringLength(lambda);
  }

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


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


  /**
   * Get quantum efficiency of %KM3NeT PMT.
   *
   * Note that the QE includes absorption in glass and gel.
   *
   * \param   lambda  wavelength of photon [nm]
   * \return          quantum efficiency
   */
  inline double getQE(const double lambda)
  {
    static JGridPolint1Function1D_t f1;
    static const double cola    = 0.95;     // collection efficiency

    if (f1.empty()) {

      f1[280] = 0.0;
      f1[285] = 0.0;
      f1[290] = 0.0;
      f1[295] = 0.2;
      f1[300] = 0.3;
      f1[305] = 2.3;
      f1[310] = 3.9;
      f1[315] = 9.1;
      f1[320] = 13.7;
      f1[325] = 19.5;
      f1[330] = 25.2;
      f1[335] = 28.9;
      f1[340] = 33.0;
      f1[345] = 34.6;
      f1[350] = 36.6;
      f1[355] = 37.5;
      f1[360] = 38.3;
      f1[365] = 38.2;
      f1[370] = 37.8;
      f1[375] = 36.8;
      f1[380] = 36.3;
      f1[385] = 36.8;
      f1[390] = 37.3;
      f1[395] = 36.6;
      f1[400] = 36.5;
      f1[405] = 36.2;
      f1[410] = 35.5;
      f1[415] = 35.2;
      f1[420] = 34.7;
      f1[425] = 34.1;
      f1[430] = 33.7;
      f1[435] = 33.2;
      f1[440] = 32.6;
      f1[445] = 31.8;
      f1[450] = 31.1;
      f1[455] = 30.1;
      f1[460] = 29.1;
      f1[465] = 28.4;
      f1[470] = 27.6;
      f1[475] = 27.0;
      f1[480] = 26.2;
      f1[485] = 25.8;
      f1[490] = 25.2;
      f1[495] = 24.0;
      f1[500] = 22.8;
      f1[505] = 20.9;
      f1[510] = 19.1;
      f1[515] = 17.1;
      f1[520] = 15.3;
      f1[525] = 14.0;
      f1[530] = 12.9;
      f1[535] = 12.3;
      f1[540] = 11.6;
      f1[545] = 11.3;
      f1[550] = 10.7;
      f1[555] = 10.3;
      f1[560] = 9.6;
      f1[565] = 9.1;
      f1[570] = 8.6;
      f1[575] = 8.0;
      f1[580] = 7.3;
      f1[585] = 6.6;
      f1[590] = 5.9;
      f1[595] = 5.3;
      f1[600] = 4.7;
      f1[605] = 4.0;
      f1[610] = 3.4;
      f1[615] = 2.9;
      f1[620] = 2.5;
      f1[625] = 2.1;
      f1[630] = 1.6;
      f1[635] = 1.3;
      f1[640] = 1.0;
      f1[645] = 0.8;
      f1[650] = 0.6;
      f1[655] = 0.5;
      f1[660] = 0.4;
      f1[665] = 0.3;
      f1[670] = 0.2;
      f1[675] = 0.1;
      f1[680] = 0.1;
      f1[685] = 0.1;
      f1[690] = 0.0;
      f1[695] = 0.0;
      f1[700] = 0.0;
      f1[705] = 0.0;

      f1.compile();

      f1.setExceptionHandler(new JGridPolint1Function1D_t::JDefaultResult(0.0));
    }

    return f1(lambda) * 0.01 * cola;
  }


  /**
   * 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