JPMTAnalogueSignalProcessor.hh

Go to the documentation of this file.

#ifndef __JDETECTOR__JPMTANALOGUESIGNALPROCESSOR__
#define __JDETECTOR__JPMTANALOGUESIGNALPROCESSOR__

#include <istream>
#include <cmath>
#include <limits>

#include "TRandom3.h"
#include "TMath.h"

#include "JLang/JException.hh"
#include "JDetector/JCalibration.hh"
#include "JDetector/JPMTParameters.hh"
#include "JDetector/JPMTSignalProcessorInterface.hh"
#include "JDetector/JPMTTransitTimeGenerator.hh"

/**
 * \file
 *
 * PMT analogue signal processor.
 * \author mdejong
 */
namespace JDETECTOR {

  using JLANG::JValueOutOfRange;


  /**
   * PMT analogue signal processor.
   *
   * This class provides for an implementation of the JDETECTOR::JPMTSignalProcessorInterface
   * using a specific model for the analogue pulse of the PMT.\n 
   * In this, the leading edge of the analogue pulse from the PMT is assumed to be a Gaussian and the tail an exponential.\n
   * The width of the Gaussian is referred to as the rise time and 
   * the inverse slope of the exponential to the decay time.\n
   * The two functions are matched at a point where the values and first derivatives are identical.\n
   *
   * Note that the decay time is related to the rise time via the specification JDETECTOR::TIME_OVER_THRESHOLD_NS.
   *
   * The charge distribution is assumed to be a Gaussian which is centered at the specified gain
   * and truncated by the specified threshold.
   *
   * The transition times are generated according the specified spread as follows.
   *  - If the specified transition-time spread (TTS) is negative, 
   *    the transition times are generated according to measurements (see method JDETECTOR::getTransitionTime).\n
   *    In this, the negated integral value of the TTS corresponds to the option 
   *    which in turn corresponds to the detector identifier of the measurements.
   *  - If the specified TTS is positive,
   *    the transition times are generated according a Gaussian with a sigma equals to the given TTS.
   *  - If the TTS is zero, the transition times are generated without any spread.
   */
  struct JPMTAnalogueSignalProcessor :
    public JPMTSignalProcessorInterface,
    public JPMTParameters
  {
    /** 
     * Threshold domain specifiers
     */
    enum JThresholdDomain {
      BELOW_THRESHOLD  = -1,     //!< below  threshold
      THRESHOLDBAND    =  0,     //!< inside threshold band
      ABOVE_THRESHOLD  =  2      //!< above  threshold
    };


    /**
     * Constructor.
     *
     * \param  parameters     PMT parameters
     */
    JPMTAnalogueSignalProcessor(const JPMTParameters& parameters = JPMTParameters()) :
      JPMTSignalProcessorInterface(),
      JPMTParameters(parameters),
      decayTime_ns(0.0),
      t1(0.0),
      y1(0.0),
      x1(std::numeric_limits<double>::max())
    {
      configure();
    }

    
    /**
     * Configure internal parameters.
     *
     * This method provides the implementations for
     *  - matching of the leading edge of the analogue pulse (Gaussian) and the tail (exponential); and
     *  - determination of number of photo-electrons above which the time-over-threshold 
     *    linearly depends on the number of photo-electrons (apart from saturation).
     *
     * Note that this method will throw an error if the value of the rise time (i.e.\ width of the Gaussian)
     * is too large with respect to the specification JDETECTOR::TIME_OVER_THRESHOLD_NS.
     */ 
    void configure()
    {
      static const int    N         = 100;
      static const double precision = 1.0e-4;

      // check thresholdband

      if (threshold - thresholdBand < getTH0()) {
        THROW(JValueOutOfRange, "JPMTAnalogueSignalProcessor::configure(): Invalid thresholdband [npe] " << thresholdBand);
      }

      // check rise time      

      if (riseTime_ns > getMaximalRiseTime(threshold)) {
        THROW(JValueOutOfRange, "JPMTAnalogueSignalProcessor::configure(): Invalid rise time [ns] " << riseTime_ns);
      } 

      // decay time

      const double y  = -log(threshold);

      const double a  =  y;
      const double b  =  riseTime_ns * sqrt(2.0*y) - TIME_OVER_THRESHOLD_NS;
      const double c  =  0.5*riseTime_ns*riseTime_ns;
      const double Q  =  b*b - 4.0*a*c;

      if (Q > 0.0) 
        decayTime_ns = (-b + sqrt(Q)) / (2.0*a);
      else
        decayTime_ns =  -b / (2.0*a);

      // fix matching of Gaussian and exponential

      const double x  =  riseTime_ns / decayTime_ns;

      t1 = riseTime_ns*x;
      y1 = exp(-0.5*x*x);

      // determine transition point to linear dependence of time-over-threshold as a function of number of photo-electrons
      
      const double xs = saturation;              // disable saturation

      saturation = 1.0e50;

      x1 = std::numeric_limits<double>::max();   // disable linearisation

      double xmin = 1.0;
      double xmax = 1.0 / (getDerivative(1.0) * slope);

      for (int i = 0; i != N; ++i) {

        const double x = 0.5 * (xmin + xmax);
        const double u = getDerivative(x) * slope;

        if (fabs(1.0 - u) < precision) {
          break;
        }

        if (u < 1.0)
          xmin = x;
        else
          xmax = x;
      }

      x1 = 0.5 * (xmin + xmax);

      saturation = xs;                           // restore saturation
    }


    /**
     * Get decay time.
     *
     * \return                decay time [ns]
     */
    double getDecayTime() const
    {
      return decayTime_ns;
    }


    /**
     * Get time at transition point from Gaussian to exponential.
     *
     * \return                time [ns]
     */
    double getT1() const
    {
      return t1;
    }


    /**
     * Get amplitude at transition point from Gaussian to exponential.
     *
     * \return                amplitude [npe]
     */
    double getY1() const
    {
      return y1;
    }


    /**
     * Get transition point from a model-dependent to linear relation between time-over-threshold and number of photo-electrons.
     *
     * \return                number of photo-electrons [npe]
     */
    double getStartOfLinearisation() const
    {
      return x1;
    }


    /**
     * Get amplitude at given time for a one photo-electron pulse.
     *
     * \param  t1_ns          time [ns]
     * \return                amplitude [npe]
     */
    double getAmplitude(const double t1_ns) const
    {
      if (t1_ns < t1) {

        const double x = t1_ns / riseTime_ns;
        
        return exp(-0.5*x*x);                                // Gaussian

      } else {

        const double x = t1_ns / decayTime_ns;

        return exp(-x) / y1;                                 // exponential
      }
    }


    /**
     * Get time to pass from threshold to top of analogue pulse.\n
     * In this, the leading edge of the analogue pulse is assumed to be Gaussian.
     *
     * \param  npe            number of photo-electrons
     * \param  th             threshold [npe]
     * \return                time [ns]
     */
    double getRiseTime(const double npe, const double th) const
    {
      return riseTime_ns * sqrt(2.0*log(npe/th));            // Gaussian
    }


    /**
     * Get time to pass from top of analogue pulse to threshold.\n
     * In this, the trailing edge of the analogue pulse is assumed to be exponential.
     *
     * \param  npe            number of photo-electrons
     * \param  th             threshold [npe]
     * \return                time [ns]
     */
    double getDecayTime(const double npe, const double th) const
    {
      if (npe*y1 > th)
        return decayTime_ns * (log(npe/th) - log(y1));       // exponential
      else
        return riseTime_ns * sqrt(2.0*log(npe/th));          // Gaussian
    }


    /**
     * Get maximal rise time for given threshold.
     *
     * Note that the rise time is entirely constrained by the specification JDETECTOR::TIME_OVER_THRESHOLD_NS.
     *
     * \param  th             threshold [npe]
     * \return                rise time [ns]
     */
    static double getMaximalRiseTime(const double th)
    {
      if (th > 0.0 && th < 1.0)
        return 0.5 * TIME_OVER_THRESHOLD_NS / sqrt(-2.0*log(th));
      else
        THROW(JValueOutOfRange, "JPMTAnalogueSignalProcessor::getMaximalRiseTime(): Invalid threshold " << th);
    }


    /**
     * Get time-over-threshold with saturation.
     *
     * \param  tot_ns         time-over-threshold without saturation 
     * \return                time-over-threshold with saturation 
     */
    double applySaturation(const double tot_ns) const
    {
      return saturation / sqrt(tot_ns*tot_ns + saturation*saturation) * tot_ns;
    }


    /**
     * Get time-over-threshold without saturation.
     *
     * \param  tot_ns         time-over-threshold with saturation 
     * \return                time-over-threshold without saturation 
     */
    double removeSaturation(const double tot_ns) const
    {
      if (tot_ns < saturation)
        return saturation / sqrt(saturation*saturation - tot_ns*tot_ns) * tot_ns;
      else
        return std::numeric_limits<double>::max();
        //THROW(JValueOutOfRange, "Time-over-threshold exceeds saturation " << tot_ns << " >= " << saturation);
    }


    /**
     * Get derivative of saturation factor.
     *
     * \param  tot_ns         time-over-threshold without saturation 
     * \return                derivative of saturation factor
     */
    double getDerivativeOfSaturation(const double tot_ns) const
    {
      return saturation * saturation * saturation / ((saturation*saturation + tot_ns*tot_ns) * sqrt(saturation*saturation + tot_ns*tot_ns));
    }


    /**
     * Get gain spread for given number of photo-electrons.
     *
     * \param  NPE            number of photo-electrons
     * \return                gain spread
     */
    double getGainSpread(int NPE) const
    {
      return sqrt((double) NPE * gain) * gainSpread;
    }


    /**
     * Get integral of probability.
     *
     * \param  xmin           minimum number of photo-electrons
     * \param  xmax           maximum number of photo-electrons
     * \param  NPE            true    number of photo-electrons
     * \return                probability
     */
    double getIntegralOfChargeProbability(const double xmin, const double xmax, const int NPE) const
    {
      double zmin = xmin;
      double zmax = xmax;

      const double th = threshold - thresholdBand;
      const double f  = gainSpread * gainSpread;

      if (zmin < th) { zmin = th; }
      if (zmax < th) { zmax = th; }
      
      double norm   = 0.0;
      double cumulP = 0.0;
      double weight = (PunderAmplified > 0.0 ? pow(PunderAmplified, NPE) : 1.0);

      for (int k = (PunderAmplified > 0.0 ? 0 : NPE); k <= NPE; ++k) {

        const double mu    = ((NPE-k) * f * gain) + (k * gain);
        const double sigma = sqrt((NPE-k) * f + k) * getGainSpread(1);

        norm   += weight * (0.5 * TMath::Erfc((th   - mu) / sqrt(2.0) / sigma));
        cumulP += weight * (0.5 * TMath::Erfc((zmin - mu) / sqrt(2.0) / sigma) -
                            0.5 * TMath::Erfc((zmax - mu) / sqrt(2.0) / sigma));

        weight *= ( (NPE - k) / ((double) (k+1)) *
                    (1.0 - PunderAmplified) / PunderAmplified);
      }

      return cumulP / norm;
    }

    
    /**
     * Get integral of probability in specific threshold domain
     *
     * \param  domain         threshold domain
     * \param  NPE            true number of photo-electrons
     * \return                probability
     */
    double getIntegralOfChargeProbability(const JThresholdDomain domain, const int NPE) const
    {
      switch (domain) {

      case ABOVE_THRESHOLD:
        return 1.0 - getIntegralOfChargeProbability(threshold-thresholdBand, threshold, NPE);

      case THRESHOLDBAND:
        return getIntegralOfChargeProbability(threshold-thresholdBand, threshold, NPE);

      default:
        return 0.0;
      }
    }


    /**
     * Set PMT parameters.
     *
     * \param  parameters              PMT parameters
     */
    void setPMTParameters(const JPMTParameters& parameters)
    {
      static_cast<JPMTParameters&>(*this).setPMTParameters(parameters);

      configure();
    }


    /**
     * Read PMT signal from input.
     *
     * \param  in             input stream
     * \param  object         PMT signal
     * \return                input stream
     */    
    friend std::istream& operator>>(std::istream& in, JPMTAnalogueSignalProcessor& object)
    {
      in >> static_cast<JPMTParameters&>(object);
      
      object.configure();

      return in;
    }


    /**
     * Get threshold domain.
     *
     * \param  npe            number of photo-electrons
     * \return                threshold domain
     */
    JThresholdDomain getThresholdDomain(const double npe) const
    {
      if (npe > threshold) {

        return ABOVE_THRESHOLD;
        
      } else if (npe > threshold - thresholdBand) {

        return THRESHOLDBAND;
        
      } else {

        return BELOW_THRESHOLD;
      }
    }


    /**
     * Apply relative QE.
     *
     * \return                true if accepted; false if rejected
     */
    virtual bool applyQE() const override 
    {
      if      (QE <= 0.0)
        return false;
      else if (QE <  1.0)
        return gRandom->Rndm() < QE;
      else
        return true;
    }


    /**
     * Get randomised time according transition time distribution.
     *
     * \param  t_ns           time [ns]
     * \return                time [ns]
     */
    virtual double getRandomTime(const double t_ns) const override 
    {
      if      (TTS_ns < 0.0)
        return t_ns + getTransitionTime(gRandom->Rndm(), -lrint(TTS_ns));
      else if (TTS_ns > 0.0)
        return gRandom->Gaus(t_ns, TTS_ns);
      else
        return t_ns;
    }
    

    /**
     * Compare arrival times of photo-electrons.
     * This implementation uses the internal rise time as two photo-electron resolution.
     *
     * Two (or more) photo-electrons are merged if they are comparable.
     *
     * \param  first          first  photo-electron
     * \param  second         second photo-electron
     * \return                true if arrival times of photo-electrons are within two photo-electron resolution; else false
     */
    virtual bool compare(const JPhotoElectron& first, const JPhotoElectron& second) const override 
    {
      return second.t_ns < first.t_ns + riseTime_ns;
    }


    /**
     * Get randomised charge according to gain and gain spread.
     *
     * \param  NPE            number of photo-electrons
     * \return                number of photo-electrons
     */
    virtual double getRandomCharge(const int NPE) const override 
    {
      double q;

      do {

        // Determine which contribution to sample from
        int k = NPE;

        if (PunderAmplified > 0.0) {

          const double X = gRandom->Uniform();
          double weight  = pow(1.0 - PunderAmplified, NPE);
        
          for (double sum_p = 0.0; k > 0; --k) {
            
            sum_p += weight;
            if (sum_p > X) { break; }
            
            weight *= ( k / ((double) (NPE - (k-1))) *
                        PunderAmplified / (1.0 - PunderAmplified) );
          }
        } 

        // Sample from chosen contribution
        const double f     = gainSpread * gainSpread;
        const double mu    = ((NPE-k) * f * gain) + (k * gain);
        const double sigma = sqrt((NPE-k) * f + k) * getGainSpread(1);

        q = gRandom->Gaus(mu,sigma);

      } while (q < 0.0);
      
      return q;
    }


    /**
     * Get probability density for given charge.
     *
     * \param  npe            observed number of photo-electrons
     * \param  NPE            true     number of photo-electrons
     * \return                probability [npe^-1]
     */
    virtual double getChargeProbability(const double npe, const int NPE) const override 
    {
      if (getThresholdDomain(npe) > BELOW_THRESHOLD) {

        const double f = gainSpread * gainSpread;

        double norm   = 0.0;
        double prob   = 0.0;
        double weight = (PunderAmplified > 0.0 ? pow(PunderAmplified, NPE) : 1.0);

        for (int k = (PunderAmplified > 0.0 ? 0 : NPE); k <= NPE; ++k) {
          
          const double mu    = ((NPE-k) * f * gain) + (k * gain);
          const double sigma = sqrt((NPE-k) * f + k) * getGainSpread(1);
          
          norm += weight * (0.5 * TMath::Erfc(((threshold-thresholdBand) - mu) / sqrt(2.0) / sigma));
          prob += weight * TMath::Gaus(npe, mu, sigma, kTRUE);
          
          weight *= ( (NPE - k) / ((double) (k+1)) *
                      (1.0 - PunderAmplified) / PunderAmplified );
        }
        
        return prob / norm;
      }

      return 0.0;
    }


    /**
     * Apply threshold.
     *
     * \param  npe            number of photo-electrons
     * \return                true if pass; else false
     */
    virtual bool applyThreshold(const double npe) const 
    {
      return getThresholdDomain(npe) > BELOW_THRESHOLD;
    }


    /**
     * Get time to reach threshold.
     *
     * Note that the rise time is defined to be zero for a one photo-electron signal.
     *
     * \param  npe            number of photo-electrons
     * \return                time [ns]
     */
    virtual double getRiseTime(const double npe) const override 
    {
      if (slewing) {

        switch (getThresholdDomain(npe)) {

        case THRESHOLDBAND:
          return ((getRiseTime(npe, getTH0()) - getRiseTime(npe, threshold-thresholdBand)) -
                  (getRiseTime(1.0, getTH0()) - getRiseTime(1.0, threshold-thresholdBand))) + this->mean_ns;

        case ABOVE_THRESHOLD: 
          return ((getRiseTime(npe, getTH0()) - getRiseTime(npe, threshold)) -
                  (getRiseTime(1.0, getTH0()) - getRiseTime(1.0, threshold)));

        default:
          THROW(JValueOutOfRange, "JPMTAnalogueSignalProcessor::getRiseTime: Invalid charge " << npe);
        }

      } else {

        return 0.0;
      }
    }


    /**
     * Get time-over-threshold (ToT).
     *
     * \param  npe            number of photo-electrons
     * \return                ToT [ns]
     */
    virtual double getTimeOverThreshold(const double npe) const override 
    {
      switch (getThresholdDomain(npe)) {

      case THRESHOLDBAND: {
        
        return gRandom->Gaus(mean_ns, sigma_ns);
      }

      case ABOVE_THRESHOLD: {
        
        double tot = 0.0;

        if (npe*y1 <= threshold) { 
          
          tot += getRiseTime(npe, threshold);                           // Gaussian
          tot += getRiseTime(npe, threshold);                           // Gaussian

        } else if (npe <= getStartOfLinearisation()) {

          tot += getRiseTime (npe, threshold);                          // Gaussian
          tot += getDecayTime(npe, threshold);                          // exponential

        } else {

          tot += getRiseTime (getStartOfLinearisation(), threshold);    // Gaussian
          tot += getDecayTime(getStartOfLinearisation(), threshold);    // exponential
          
          tot += slope * (npe - getStartOfLinearisation());             // linear
        }
        
        return applySaturation(tot);
      }
        
      default: {

        THROW(JValueOutOfRange, "JPMTAnalogueSignalProcessor::getTimeOverThreshold: Invalid charge " << npe);
      }}
    }


    /**
     * Get derivative of number of photo-electrons to time-over-threshold.
     *
     * \param  npe            number of photo-electrons 
     * \return                dnpe/dToT [ns^-1]
     */
    virtual double getDerivative(const double npe) const override 
    {
      switch (getThresholdDomain(npe)) {

      case ABOVE_THRESHOLD: {

        const double z  = riseTime_ns / sqrt(2.0 * log(npe/threshold));

        double y = 0.0;

        if (npe*y1 > threshold) {
          
          if (npe <= getStartOfLinearisation())
            y = npe / (z + decayTime_ns);           // Gaussian + exponential
          else
            y = 1.0 / slope;                        // linear

        } else {

          y = npe / (2.0 * z);                      // Gaussian + Gaussian
        }

        const double tot_ns = getTimeOverThreshold(npe);

        return y / getDerivativeOfSaturation(removeSaturation(tot_ns));
      }

      default:
        return 0.0;
      }
    }


    /**
     * Probability that a hit survives the simulation of the PMT.
     *
     * \param  NPE            number of photo-electrons
     * \return                probability
     */
    virtual double getSurvivalProbability(const int NPE) const override 
    {
      if        (QE <= 0.0) {

        return 0.0;

      } else if (QE <= 1.0) {

        double P = 0.0;
        
        const double f = gainSpread * gainSpread;

        for (int i = 1; i <= NPE; ++i) {

          const double p = (TMath::Binomial(NPE, i) *
                            TMath::Power(QE, i) * TMath::Power(1.0 - QE, NPE - i));
          
          double Ptotal = 0.0;
          double Pabove = 0.0;
          double weight = (PunderAmplified > 0.0 ? pow(PunderAmplified, i) : 1.0);

          for (int k = (PunderAmplified > 0.0 ? 0 : NPE); k <= NPE; ++k) {

            const double mu    = ((i-k) * f * gain) + (k * gain);
            const double sigma = sqrt((i-k) * f + k) * getGainSpread(1);

            Ptotal += weight * (0.5 * TMath::Erfc((0.0                       - mu) / sqrt(2.0) / sigma));
            Pabove += weight * (0.5 * TMath::Erfc((threshold - thresholdBand - mu) / sqrt(2.0) / sigma));

            weight *= ( (i - k) / ((double) (k+1)) *
                        (1.0 - PunderAmplified) / PunderAmplified );
          }

          P += p*Pabove/Ptotal;
        }
        
        return P;

      } else {

        THROW(JValueOutOfRange, "JPMTAnalogueSignalProcessor::getSurvivalProbability: Invalid QE " << QE);
      }
    }


    /**
     * Get number of photo-electrons.
     *
     * \param  tot_ns       time over threshold [ns]
     * \return              number of photo-electrons
     */
    virtual double getNPE(const double tot_ns) const override
    {
      double tot = std::numeric_limits<double>::max();
      
      if (tot_ns < saturation) {
        tot = removeSaturation(tot_ns);
      }

      const double TOT = (getRiseTime (getStartOfLinearisation(), threshold)  + 
                          getDecayTime(getStartOfLinearisation(), threshold));
      
      if        (tot <= 2*getRiseTime(threshold/y1,threshold)) {              // Gaussian + Gaussian
        
        return threshold * exp(tot*tot/riseTime_ns/riseTime_ns/8.0);
        
      } else if (tot <= TOT) {                                                // Gaussian + Exponential

        const double a = decayTime_ns;
        const double b = sqrt(2.0) * riseTime_ns;
        const double c = -(decayTime_ns*log(y1) + tot);
        const double z = (-b + sqrt(b*b - 4*a*c)) / (2*a);

        return threshold * exp(z*z);
        
      } else {                                                                // linear
          
        return getStartOfLinearisation() + (tot - TOT) / slope;
      }
    }


    /**
     * Get probability of having a pulse with specific time-over-threshold
     *
     * \param  tot_ns           time-over-threshold with saturation [ns]
     * \param  NPE              true number of photo-electrons
     * \return                  probability [ns^-1]
     */
    double getTimeOverThresholdProbability(const double tot_ns, const int NPE) const
    {
      const double PthBand  = getIntegralOfChargeProbability(THRESHOLDBAND, NPE);

      const double npe      = getNPE(tot_ns);
      const double y        = getChargeProbability(npe, NPE);
      const double v        = getDerivative(npe);
      
      const double RthBand  = PthBand * TMath::Gaus(tot_ns, mean_ns, sigma_ns, kTRUE);
      const double RaboveTh = y * v;

      return RthBand + RaboveTh;
    }


    /**
     * Get cumulative probability of time-over-threshold distribution
     *
     * \param  Tmin             minimum time-over-threshold (with saturation) [ns]
     * \param  Tmax             maximum time-over-threshold (with saturation) [ns]
     * \param  NPE              true number of photo-electrons
     * \return                  probability [ns^-1]
     */
    double getIntegralOfTimeOverThresholdProbability(const double Tmin, const double Tmax, const int NPE) const
    {

      const double PthBand  = getIntegralOfChargeProbability(THRESHOLDBAND, NPE);

      const double IthBand  = PthBand * (0.5 * TMath::Erfc((Tmin - mean_ns) / sqrt(2.0) / sigma_ns) -
                                         0.5 * TMath::Erfc((Tmax - mean_ns) / sqrt(2.0) / sigma_ns)); 
      const double IaboveTh = getIntegralOfChargeProbability(getNPE(Tmin), getNPE(Tmax), NPE);

      return IthBand + IaboveTh;
    }


    /**
     * Get lower threshold for rise time evaluation.
     *
     * \return                threshold [npe]
     */
    static double getTH0()
    {
      return 0.1;
    }


    /**
     * Get upper threshold for rise time evaluation.
     *
     * \return                threshold [npe]
     */
    static double getTH1()
    {
      return 0.9;
    }
 
  protected:

    double decayTime_ns;                //!< decay time               [ns]
    double t1;                          //!< time      at match point [ns]
    double y1;                          //!< amplitude at match point [npe]
    /** 
     * Transition point from a logarithmic to a linear relation 
     * between time-over-threshold and number of photo-electrons.\n
     * Measurements by B. Schermer and R. Bruijn at Nikhef.
     */
    double x1;
  };


  /**
   * Get time-over-threshold probability.
   *
   * \param  pmt             PMT signal processor
   * \param  tot_ns          time-over-threshold [ns]
   * \param  NPE             expected number of photo-electrons
   * \param  precision       precision
   * \return                 probability
   */
  inline double getTimeOverThresholdProbability(const JPMTAnalogueSignalProcessor& pmt,
                                                const double tot_ns,
                                                const double NPE,
                                                const double precision = 1.0e-4)
  {
    int i = (int) (NPE - 5.0 * sqrt(NPE) - 0.5);
    int M = (int) (NPE + 5.0 * sqrt(NPE) + 0.5);

    if (i < 1) { i = 1; }
    if (M < i) { M = i; }

    double p = NPE * exp(-NPE) / (double) 1;

    for (int __i = 1; __i != i; ++__i) {
      p *= NPE / __i;
    }

    double P = 0.0;

    for (double p0 = 0.0; (p >= p0 || p > precision) && i != M; ++i, p0 = p, p *= NPE / (double) i) {
      P += pmt.getTimeOverThresholdProbability(tot_ns, i) * p;
    }

    return P;
  }
}

#endif