Antares.hh

Go to the documentation of this file.

#ifndef __JPHYSICS__ANTARES__
#define __JPHYSICS__ANTARES__

#include <map>
#include <cmath>

#include "JTools/JConstants.hh"
#include "JTools/JFunction1D_t.hh"


/**
 * \file
 *
 * Optical properties of Antares deep-sea site.
 * \author mdejong
 */
namespace ANTARES {

  using JTOOLS::PI;
  using JTOOLS::JPolint1Function1D_t;
  using JTOOLS::JGridPolint1Function1D_t;

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

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


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


  /**
   * Scattering length
   *
   * \param  lambda     wavelength of light [nm]
   * \return            absorption 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;
  }

  
  /**
   * Auxiliary method to describe light scattering in water (Heneyey-Greenstein)
   *
   * \param  g      angular dependence parameter
   * \param  x      cosine scattering angle
   * \return        probability
   */
  inline double henyey_greenstein(const double g,
                                  const double x)
    
  {
    const double a0 = (1.0 - g*g) / (4*PI);
    const double y  =  1.0 + g*g - 2.0*g*x;
    
    return a0 / (y*sqrt(y));
  }
  

  /**
   * Auxiliary method to describe light scattering in water (Heneyey-Greenstein)
   *
   * \param  x      cosine scattering angle
   * \return        probability
   */
  inline double henyey_greenstein(const double x)
  {
    static const double g = 0.924;
    
    return henyey_greenstein(g, x);
  }


  /**
   * Auxiliary method to describe light scattering in water (Rayleigh)
   *
   * \param  a      angular dependence parameter
   * \param  x      cosine scattering angle
   * \return        probability
   */
  inline double rayleigh(const double a,
                         const double x)

  {
    const double a0 = 1.0 / (1.0 + a/3.0) / (4*PI);

    return a0 * (1.0 + a*x*x); 
  }


  /**
   * Auxiliary method to describe light scattering in water (Rayleigh)
   *
   * \param  x      cosine scattering angle
   * \return        probability
   */
  inline double rayleigh(const double x)
  {
    return rayleigh(0.853, x);
  }


  /**
   * Model specific function to describe light scattering in water (f4)
   *
   * \param  x      cosine scattering angle
   * \return        probability
   */
  inline double f4(const double x)
  {
    static const double g1 = 0.77;
    static const double g2 = 0.00;
    static const double f  = 1.00;

    return f * henyey_greenstein(g1,x)  +  (1.0 - f) * henyey_greenstein(g2,x);
  }


  /**
   * Model specific function to describe light scattering in water (p00075)
   *
   * \param  x      cosine scattering angle
   * \return        probability
   */
  inline double p00075(const double x)
  {
    static const double g = 0.924;
    static const double f = 0.17;

    return f * rayleigh(x)  +  (1.0 - f) * henyey_greenstein(g,x);
  }

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


  /**
   * Angular acceptence of Antares 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 acceptence of Antares 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;
  }


  /**
   * Angular acceptence of Antares PMT
   *
   * \param  x      cosine of angle of incidence
   * \return        probability
   */
  inline double getAngularAcceptance(const double x)
  {
    return genova(x);
  }


  /**
   * Quantum efficiency of 10-inch Hamamatsu 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;
  }


  /**
   * Quantum efficiency of 10-inch Hamamatsu PMT
   *
   * \param   lambda  wavelength of photon [nm]
   * \return          quantum efficiency
   */
  inline double getQE(const double lambda)
  {
    return getQE(lambda, true);
  }


  /**
   * Photo-cathode area 10 inch 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() * getAngularAcceptance(x);
  }
}

#endif