JLED.hh

Go to the documentation of this file.

#ifndef __JPHYSICS__JLED__
#define __JPHYSICS__JLED__

#include "JPhysics/JAbstractPMT.hh"
#include "JPhysics/JAbstractMedium.hh"
#include "JPhysics/JDispersionInterface.hh"
#include "JPhysics/JDispersion.hh"
#include "JPhysics/JConstants.hh"
#include "JTools/JElement.hh"
#include "JTools/JQuadrature.hh"


/**
 * \author mdejong
 */

namespace JPHYSICS {}
namespace JPP { using namespace JPHYSICS; }

namespace JPHYSICS {


  using JTOOLS::JQuadrature;
  using JTOOLS::JGaussLegendre;
  using JTOOLS::JCotangent;

  typedef JTOOLS::JElement2D<double, double> JElement2D_t;
  typedef JTOOLS::JElement3D<double, double> JElement3D_t;


  /**
   * Interface for emission profile from LED.
   */
  class JAbstractLED {
  public:

    /**
     * Virtual destructor.
     */
    virtual ~JAbstractLED()
    {}


    /**
     * Light yield from LED (number of p.e. per unit solid angle per unit time).
     *
     * \param  ct         zenith  angle of emission
     * \param  phi        azimuth angle of emission
     * \param  dt         time of emission [ns]
     * \return            d^2P / dOmega dt [npe/ns/sr]
     */
    virtual double getLightFromLED(const double ct, 
                                   const double phi,
                                   const double dt) const = 0;
  };


  /**
   * Probability Density Functions of the time response of a PMT.
   */
  class JLED : 
    public virtual JDispersionInterface, 
    public JAbstractPMT, 
    public JAbstractLED,
    public JAbstractMedium
  {
  public:
    /**
     * Constructor.
     *
     * \param  lambda               wavelength photon [nm]
     * \param  Tmin_ns              minimal time of emmision [ns]
     * \param  Tmax_ns              minimal time of emmision [ns]
     * \param  engine               scattering angle integrator
     * \param  numberOfPoints       number    of points for integration
     * \param  epsilon              precision of points for integration
     */
    JLED(const double lambda,
         const double Tmin_ns,
         const double Tmax_ns,
         const JQuadrature& engine         = JCotangent(20),
         const int          numberOfPoints = 20,
         const double       epsilon        = 1e-12) :
      wavelength(lambda),
      tmin(Tmin_ns),
      tmax(Tmax_ns),
      main_engine(JGaussLegendre(numberOfPoints,epsilon)),
      beta_engine(engine)
    {
      const double dp = PI / numberOfPoints;

      for (double x = 0.5*dp; x < PI; x += dp) {
        phi_engine.push_back(JElement3D_t(cos(x),sin(x), dp));
      }
    }


    /**
     * Probability density function for direct light from LED.
     *
     * \param  D_m        distance between LED and PMT [m]
     * \param  cd         cosine angle LED orientation and LED - PMT position
     * \param  theta      zenith  angle orientation PMT
     * \param  phi        azimuth angle orientation PMT
     * \param  t_ns       time difference relative to direct light [ns]
     * \return            dP/dt [npe/ns]
     */
    double getDirectLightFromLED(const double D_m,
                                 const double cd,
                                 const double theta,
                                 const double phi,
                                 const double t_ns) const
    {
      const double sd =  sqrt((1.0 + cd)*(1.0 - cd));
      const double d  =  D_m;                                  // photon path [m]

      const double A  =  getPhotocathodeArea();

      const double px =  sin(theta)*cos(phi);
      //const double py =  sin(theta)*sin(phi);
      const double pz =  cos(theta);

      const double ct = sd*px + cd*pz;                         // cosine angle of incidence on PMT

      const double qe    = getQE(wavelength);
      const double l_abs = getAbsorptionLength(wavelength);
      const double ls    = getScatteringLength(wavelength);

      const double U  = getAngularAcceptance(ct);
      const double V  = exp(-d/l_abs) * exp(-d/ls);            // absorption & scattering
      const double W  = A/(d*d);                               // solid angle
  
      return qe * getLightFromLED(cd, 0.0, t_ns) * U * V * W;
    }

                                 
    /**
     * Probability density function for scattered light from LED.
     *
     * \param  D_m        distance between LED and PMT [m]
     * \param  cd         cosine angle LED orientation and LED - PMT position
     * \param  theta      zenith  angle orientation PMT
     * \param  phi        azimuth angle orientation PMT
     * \param  t_ns       time difference relative to direct light [ns]
     * \return            dP/dt [npe/ns]
     */
    double getScatteredLightFromLED(const double D_m,
                                    const double cd,
                                    const double theta,
                                    const double phi,
                                    const double t_ns) const
    {
      using namespace std;

      double value = 0;

      if (t_ns >= tmin) {

        const double ng =  getIndexOfRefractionGroup(wavelength);

        const double sd =  sqrt((1.0 + cd)*(1.0 - cd));
        const double l  =  D_m;                                    // distance [m]

        const double A  =  getPhotocathodeArea();
        //const double sr =  sqrt(A/PI) / D_m;
        //const double cr =  sqrt((1.0 + sr)*(1.0 - sr));

        const double px =  sin(theta)*cos(phi);
        const double py =  sin(theta)*sin(phi);
        const double pz =  cos(theta);

        const double qx = cd*px +  0  - sd*pz;
        const double qy =         1*py;
        const double qz = sd*px +  0  + cd*pz;

        const double qe    = getQE(wavelength);
        const double l_abs = getAbsorptionLength(wavelength);
        const double ls    = getScatteringLength(wavelength);

        const double Jc = 1.0 / ls;                                // dN/dx

        const double xmin = tmin;
        const double xmax = std::min(tmax, t_ns);

        JQuadrature buffer;
        
        if (xmax > xmin) {

          for (JQuadrature::const_iterator i = main_engine.begin(); i != main_engine.end(); ++i) {

            const double t  = 0.5 * (xmax + xmin)  +  i->getX() * 0.5 * (xmax - xmin);
            const double dt = i->getY() * 0.5 * (xmax - xmin);

            buffer[t] = dt;
          }

        } else {

          buffer[tmin] = 1.0;
        }

        for (JQuadrature::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {

          const double t  = i->getX();
          const double dt = i->getY();

          const double d  = l + C*(t_ns - t)/ng;                   // photon path
          const double V  = exp(-d/l_abs);                         // absorption

          for (JQuadrature::const_iterator j = beta_engine.begin(); j != beta_engine.end(); ++j) {

            const double cb = j->getX();
            const double dc = j->getY();
            
            const double sb = sqrt((1.0 + cb)*(1.0 - cb));
            
            const double v  = 0.5 * (d + l) * (d - l) / (d - l*cb);
            const double u  = d - v;
            
            if (u <= 0) continue;
            if (v <= 0) continue;
            
            const double cts = (l*cb - v) / u;                   // cosine scattering angle
            
            const double W  = min(A/(v*v), 2.0*PI);              // solid angle
            const double Ja = getScatteringProbability(cts);     // P(cos(),phi)
            const double Jd = ng * (1.0 - cts) / C;              // dt/du
            
            const double ca = (l - v*cb) / u;
            const double sa =      v*sb  / u;
            
            for (std::vector<JElement3D_t>::const_iterator k = phi_engine.begin(); k != phi_engine.end(); ++k) {
              
              const double cp = k->getX();
              const double sp = k->getY();
              const double dp = k->getZ();

              const double dom = dc * dp * v*v / (u*u);
              
              const double ct0 = cd*ca - sd*sa*cp;
              const double ph0 = atan2(sa*sp, cd*sa*cp + sd*ca);
              
              const double vx  = -sb*cp * qx;
              const double vy  = -sb*sp * qy;
              const double vz  =  cb    * qz;
              
              const double U  =
                getAngularAcceptance(vx + vy + vz) * getLightFromLED(ct0, +ph0, t) +
                getAngularAcceptance(vx - vy + vz) * getLightFromLED(ct0, -ph0, t);
              
              value += qe * dt * dom * Ja * Jc * U * V * W / Jd;
            }
          }
        }
      }

      return value; 
    }                            

  protected:
    double wavelength;                    // wavelength photon [nm]

    double tmin;                          // minimal time of emmision [ns]
    double tmax;                          // maximal time of emmision [ns]

    JQuadrature main_engine;              // numerical integration
    JQuadrature beta_engine;              // numerical integration of scattering angle
    std::vector<JElement3D_t> phi_engine; // numerical integration of phi angle
  };


  /**
   * Probability Density Functions of the time response of a PMT (C-like interface)
   */
  class JLED_C : 
    public JLED,
    public JDispersion 
  {
  public:
    /**
     * Constructor.
     *
     * \param  Area            photo-cathode area [m^2]
     * \param  LED             pointer to interface for emission profile from LED
     * \param  QE              pointer to function for quantum efficiency of PMT
     * \param  Pmt             pointer to function for angular acceptance of PMT
     * \param  L_abs           absorption length [m]
     * \param  L_s             scattering length [m]
     * \param  Km3             pointer to model specific function to describe light scattering in water
     * \param  P_atm           ambient pressure [atm]
     * \param  wavelength      wavelength of light [nm]
     * \param  Tmin_ns         minimal time of emmision [ns]
     * \param  Tmax_ns         minimal time of emmision [ns]
     * \param  engine          scattering angle integrator
     * \param  numberOfPoints  number    of points for integration
     * \param  epsilon         precision of points for integration
     */
    JLED_C(const double Area,
           //
           // pointers to static functions
           //
           const JAbstractLED* LED,
           double (*QE)  (const double),
           double (*Pmt) (const double),
           const double L_abs,
           const double L_s,
           double (*Km3) (const double),
           //
           // parameters for base class
           //
           const double P_atm,
           const double wavelength,
           const double Tmin_ns,
           const double Tmax_ns,
           const JQuadrature& engine   = JCotangent(20),
           const int    numberOfPoints = 20,
           const double epsilon        = 1e-12) :
      JLED(wavelength, Tmin_ns, Tmax_ns, engine, numberOfPoints, epsilon),
      JDispersion(P_atm),
      A    (Area ),
      led  (LED  ),
      qe   (QE   ),
      l_abs(L_abs),
      ls   (L_s  ),
      pmt  (Pmt  ),
      km3  (Km3  )
    {}


    /**
     * Light yield from LED (number of p.e. per unit solid angle per unit time).
     *
     * \param  ct         zenith  angle of emission
     * \param  phi        azimuth angle of emission
     * \param  dt         time of emission [ns]
     * \return            d^2P / dOmega dt [npe/ns/sr]
     */
    virtual double getLightFromLED(const double ct, 
                                   const double phi, 
                                   const double dt) const
    { 
      return led->getLightFromLED(ct, phi, dt);
    }


    /**
     * Photo-cathode area of PMT.
     *
     * \return            photo-cathode area [m^2]
     */
    virtual double getPhotocathodeArea() const override 
    { 
      return A;
    }


    /**
     * Quantum efficiency of PMT (incl. absorption in glass, gel, etc.).
     *
     * \param  lambda     wavelenth [nm]
     * \return            QE
     */
    virtual double getQE(const double lambda) const override 
    { 
      return qe(lambda);
    }

 
    /**
     * Angular acceptance of PMT.
     *
     * \param  ct         cosine angle of incidence
     * \return            relative efficiency
     */
    virtual double getAngularAcceptance(const double ct) const override 
    { 
      return pmt(ct);
    }


    /**
     * Absorption length.
     *
     * \param  lambda     wavelenth [nm]
     * \return            absorption length [m]
     */
    virtual double getAbsorptionLength(const double lambda) const override 
    { 
      return l_abs;
    }


    /**
     * Scattering length.
     *
     * \param  lambda     wavelenth [nm]
     * \return            scattering length [m]
     */
    virtual double getScatteringLength(const double lambda) const override 
    { 
      return ls;
    }


    /**
     * Model specific function to describe light scattering in water
     * (integral over full solid angle normalised to one).
     *
     * \param  ct         cosine scattering angle
     * \return            probability
     */
    virtual double getScatteringProbability(const double ct) const override 
    { 
      return km3(ct);
    }

  protected:
    /**
     * photo-cathode area [m2]
     */
    const double A;

    /**
     * Pointer to interface for emission profile from LED
     */
    const JAbstractLED* led;

    /**
     * Quantum efficiency of PMT (incl. absorption in glass, gel, etc.)
     *
     * \param  lambda     wavelenth [nm]
     * \return            QE
     */
    double (*qe)(const double lambda);

    /**
     * Absorption length
     */
    double l_abs;

    /**
     * Scattering length
     */
    double ls;

    /**
     * Angular acceptance of PMT
     *
     * \param  ct     cosine angle of incidence
     * \return        relative efficiency
     */
    double (*pmt)(const double ct);

    /**
     * Model specific function to describe light scattering in water
     *
     * \param  ct     cosine scattering angle
     * \return        probability
     */
    double (*km3)(const double ct);
   }; 
}

#endif