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 "JLang/JException.hh"
#include "JMath/JMathSupportkit.hh"
#include "JDetector/JCalibration.hh"
#include "JDetector/JPMTParameters.hh"
#include "JDetector/JPMTSignalProcessorInterface.hh"
#include "JDetector/JPMTTransitTimeGenerator.hh"
 
/**
 * \file
 *
 * PMT analogue signal processor.
 * \author mdejong & bjjung & evancampenhout
 */
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 transit times are generated according the specified spread as follows.
   *  - If the specified transit-time spread (TTS) is negative, 
   *    the transit times are generated according to measurements (see method JDETECTOR::getTransitTime).\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 transit times are generated according a Gaussian with a sigma equals to the given TTS.
   *  - If the TTS is zero, the transit 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);
      } 
 
      // check probability for under-amplified signal
 
      if (PunderAmplified < 0.0) {
        THROW(JValueOutOfRange, "JPMTAnalogueSignalProcessor::configure(): Invalid probability of under-amplification " << PunderAmplified);
      }
 
      // 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();
    }
 
 
    /**
     * 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.
     *
     * \return                gain spread
     */
    double getGainSpread() const
    {
      return sqrt(gain) * gainSpread;
    }
 
 
    /**
     * 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) * getGainSpread();
    }
 
 
    /**
     * Get gain of underamplified pulses.
     * 
     * \return                gain of underamplified pulses
     */
    double getUnderamplifiedGain() const
    {
      return gain * gainSpread * gainSpread;
    }
 
 
    /**
     * Get gain of underamplified pulses.
     * 
     * \param NPE             number of photo-electrons
     * \return                gain of underamplified pulses
     */
    double getUnderamplifiedGain(int NPE) const
    {
      return ((double) NPE) * getUnderamplifiedGain();
    }
 
 
    /**
     * Get gain spread of underamplified pulses.
     * 
     * \return                gain spread of underamplified pulses
     */
    double getUnderamplifiedGainSpread() const
    {
      return sqrt(gain) * gainSpread * gainSpread;
    }
 
 
    /**
     * Get gain spread of underamplified pulses.
     * 
     * \param NPE             number of photo-electrons
     * \return                gain spread of underamplified pulses
     */
    double getUnderamplifiedGainSpread(int NPE) const
    {
      return sqrt((double) NPE) * getUnderamplifiedGainSpread();
    }
 
 
    /**
     * 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 transit-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 + getTransitTime(gRandom->Rndm(), getType());
      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 
    {
      if (NPE > 0) {
 
        if (gainSpread > 0.0) {
 
          double q;
 
          do {
 
            if        (PunderAmplified <= 0.0) {
 
              const double mu    = NPE * gain;
              const double sigma = getGainSpread(NPE);
 
              q = gRandom->Gaus(mu,sigma);
 
            } else if (PunderAmplified < 1.0) {
 
              // Determine which contribution to sample from
              // Method uses inverse transform sampling for a binomial distribution
 
              const double X = gRandom->Uniform();
              double sum_p   = 0.0;
              double weight  = pow(1-PunderAmplified, NPE);
              int k          = 0;  // k underamplified photo-electrons
      
              for (; k <= NPE; ++k) {
          
                sum_p += weight;
                if (sum_p > X) { break; }
          
                weight *= ((double) (NPE-k) / ((double) (k+1)))  *  PunderAmplified / (1.0 - PunderAmplified);
              }
 
              // Sample from chosen contribution
              const double mu    = (NPE-k) * gain  +  (k) * getUnderamplifiedGain();
              const double sigma = gainSpread * sqrt(mu);
 
              q = gRandom->Gaus(mu,sigma);
 
            } else {
 
              const double mu    = getUnderamplifiedGain(NPE);
              const double sigma = getUnderamplifiedGainSpread(NPE);
 
              q = gRandom->Gaus(mu,sigma);
            }
 
          } while (q < 0.0);
          
          return q;
        
        } else {
 
          return NPE * gain;
        }
 
      } else {
 
        return 0.0;
      }
    }
 
 
    /**
     * Get probability density for given charge.
     * The normalisation is from (threshold-thresholdband) to infinity.
     *
     * \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 ((NPE > 0) && (getThresholdDomain(npe) > BELOW_THRESHOLD)) {
 
        if        (PunderAmplified <= 0.0) {
 
          const double mu    = NPE * gain;
          const double sigma = getGainSpread(NPE);
 
          if (sigma > 0.0) {
 
            const double prob = JMATH::Gauss(npe, mu, sigma);
            const double norm = 0.5 * erfc(((threshold-thresholdBand) - mu) / sqrt(2.0) / sigma);
 
            return prob / norm;
          
          } else if (mu >= threshold - thresholdBand) {
 
            return 1.0;
          
          } else {
 
            return 0.0;
          }
        
        } else if (PunderAmplified < 1.0) {
 
          double norm   = 0.0;
          double prob   = 0.0;
          double weight = pow(1-PunderAmplified, NPE);
 
          for (int k = 0; k <= NPE; ++k) {  // k underamplified photo-electrons
      
            const double mu    = (NPE-k) * gain  +  (k) * getUnderamplifiedGain();
            const double sigma = gainSpread * sqrt(mu);
 
            if (sigma > 0.0) {
      
              prob += weight * JMATH::Gauss(npe, mu, sigma);
              norm += weight * (0.5 * erfc(((threshold-thresholdBand) - mu) / sqrt(2.0) / sigma));
            
            } else if (mu >= threshold - thresholdBand) {
 
              prob += weight;
              norm += weight;
            
            } else {
 
              return 0.0;
            }
      
            weight *= ((double) (NPE-k) / ((double) (k+1)))  *  PunderAmplified / (1.0 - PunderAmplified);
          }
    
          return prob / norm;
        
        } else {
 
          const double mu    = getUnderamplifiedGain(NPE);
          const double sigma = getUnderamplifiedGainSpread(NPE);
 
          if (sigma > 0.0) {
 
            const double prob = JMATH::Gauss(npe, mu, sigma);
            const double norm = 0.5 * erfc(((threshold-thresholdBand) - mu) / sqrt(2.0) / sigma);
 
            return prob / norm;
          
          } else if (mu >= threshold - thresholdBand) {
 
            return 1.0;
          
          } else {
 
            return 0.0;
          }
        }
 
      } else {
 
        return 0.0;
      }
    }
 
    /**
     * Get probability density for given charge.
     * The normalisation is from a configurable threshold to infinity.
     *
     * \param  npe            observed number of photo-electrons
     * \param  NPE            true     number of photo-electrons
     * \param  normTh         normalisation threshold
     * \return                probability [npe^-1]
     */
    double getChargeProbability(const double npe, const int NPE, const double normTh) const
    {
      if ((NPE > 0) && (getThresholdDomain(npe) > BELOW_THRESHOLD)) {
 
        double minCharge = normTh;
        const double th = threshold - thresholdBand;
 
        if (minCharge < th) {minCharge = th;}
 
        if        (PunderAmplified <= 0.0) {
 
          const double mu    = NPE * gain;
          const double sigma = getGainSpread(NPE);
 
          if (sigma > 0.0) {
 
            const double prob = JMATH::Gauss(npe, mu, sigma);
            const double norm = 0.5 * erfc((normTh - mu) / sqrt(2.0) / sigma);
 
            return prob / norm;
 
          } else if (mu >= normTh) {
 
            return 1.0;
          
          } else {
 
            return 0.0;
          }
        
        } else if (PunderAmplified < 1.0) {
 
          double norm   = 0.0;
          double prob   = 0.0;
          double weight = pow(1-PunderAmplified, NPE);
 
          for (int k = 0; k <= NPE; ++k) {  // k underamplified photo-electrons
      
            const double mu    = (NPE-k) * gain  +  (k) * getUnderamplifiedGain();
            const double sigma = gainSpread * sqrt(mu);
 
            if (sigma > 0.0) {
      
              prob += weight * JMATH::Gauss(npe, mu, sigma);
              norm += weight * (0.5 * erfc((normTh - mu) / sqrt(2.0) / sigma));
            
            } else if (mu >= normTh) {
 
              prob += weight;
              norm += weight;
            
            } else {
 
              return 0.0;
            }
      
            weight *= ((double) (NPE-k) / ((double) (k+1)))  *  PunderAmplified / (1.0 - PunderAmplified);
          }
    
          return prob / norm;
        
        } else {
 
          const double mu    = getUnderamplifiedGain(NPE);
          const double sigma = getUnderamplifiedGainSpread(NPE);
 
          if (sigma > 0.0) {
 
            const double prob = JMATH::Gauss(npe, mu, sigma);
            const double norm = 0.5 * erfc((normTh - mu) / sqrt(2.0) / sigma);
 
            return prob / norm;
 
          } else if (mu >= normTh) {
 
            return 1.0;
          
          } else {
 
            return 0.0;
          }
        }
      
      } else {
 
        return 0.0;
      }
    }
 
 
    /**
     * Get integral of probability.
     * The normalisation is from (threshold-thresholdband) to infinity.
     *
     * \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
    {
      if (NPE > 0) {
 
        double zmin = xmin;
        double zmax = xmax;
 
        const double th = threshold - thresholdBand;
 
        if (zmin < th) { zmin = th; }
        if (zmax < th) { zmax = th; }
 
        if        (PunderAmplified <= 0.0) {
 
          const double mu    = NPE * gain;
          const double sigma = getGainSpread(NPE);
 
          if (sigma > 0.0) {
 
            const double cumulP = (0.5 * erfc((zmin - mu) / sqrt(2.0) / sigma) -
                                   0.5 * erfc((zmax - mu) / sqrt(2.0) / sigma));
            const double norm   =  0.5 * erfc((th   - mu) / sqrt(2.0) / sigma);
 
            return cumulP / norm;
 
          } else if ((mu >= zmin) && (mu <= zmax)) {
 
            return 1.0;
          
          } else {
 
            return 0.0;
          }
 
        } else if (PunderAmplified < 1.0) {
 
          double cumulP = 0.0;
          double norm   = 0.0;
          double weight = pow(1.0 - PunderAmplified, NPE);
 
          for (int k = 0; k <= NPE; ++k) {  // k underamplified photo-electrons
 
            const double mu    = (NPE-k) * gain  +  (k) * getUnderamplifiedGain();
            const double sigma = gainSpread * sqrt(mu);
 
            if (sigma > 0.0) {
 
              cumulP += weight * (0.5 * erfc((zmin - mu) / sqrt(2.0) / sigma) -
                                  0.5 * erfc((zmax - mu) / sqrt(2.0) / sigma));
              norm   += weight * (0.5 * erfc((th   - mu) / sqrt(2.0) / sigma));
            
            } else if ((mu >= zmin) && (mu <= zmax)) {
 
              cumulP += weight;
              norm   += weight;
            
            } else {
 
              return 0.0;
            }
 
            weight *= ((double) (NPE-k) / ((double) (k+1)))  *  PunderAmplified / (1.0 - PunderAmplified);
          }
 
          return cumulP / norm;
        
        } else {
 
          const double mu    = getUnderamplifiedGain(NPE);
          const double sigma = getUnderamplifiedGainSpread(NPE);
 
          if (sigma > 0.0) {
 
            const double cumulP = (0.5 * erfc((zmin - mu) / sqrt(2.0) / sigma) -
                                   0.5 * erfc((zmax - mu) / sqrt(2.0) / sigma));
            const double norm   =  0.5 * erfc((th   - mu) / sqrt(2.0) / sigma);
 
            return cumulP / norm;
 
          } else if ((mu >= zmin) && (mu <= zmax)) {
 
            return 1.0;
          
          } else {
 
            return 0.0;
          }
        }
 
      } else {
 
        return 0.0;
      }
    }
 
 
    /**
     * Get integral of probability.
     * The normalisation is from a configurable threshold to infinity.
     *
     * \param  xmin           minimum number of photo-electrons
     * \param  xmax           maximum number of photo-electrons
     * \param  NPE            true    number of photo-electrons
     * \param  normTh         normalisation threshold
     * \return                probability
     */
    double getIntegralOfChargeProbability(const double xmin, const double xmax, const int NPE, const double normTh) const
    {
      if (NPE > 0) {
 
        double minCharge = normTh;
        const double th = threshold - thresholdBand;
 
        if (minCharge < th) {minCharge = th;}
 
        double zmin = xmin;
        double zmax = xmax;
 
        if (zmin < minCharge) { zmin = minCharge; }
        if (zmax < minCharge) { zmax = minCharge; }
 
        if        (PunderAmplified <= 0.0) {
 
          const double mu    = NPE * gain;
          const double sigma = getGainSpread(NPE);
 
          if (sigma > 0.0) {
 
            const double cumulP = (0.5 * erfc((zmin - mu) / sqrt(2.0) / sigma) -
                                   0.5 * erfc((zmax - mu) / sqrt(2.0) / sigma));
            const double norm   =  0.5 * erfc((normTh   - mu) / sqrt(2.0) / sigma);
 
            return cumulP / norm;
          
          } else if ((mu >= zmin) && (mu <= zmax)) {
 
            return 1.0;
          
          } else {
 
            return 0.0;
          }
 
        } else if (PunderAmplified < 1.0) {
 
          double cumulP = 0.0;
          double norm   = 0.0;
          double weight = pow(1.0 - PunderAmplified, NPE);
 
          for (int k = 0; k <= NPE; ++k) {  // k underamplified photo-electrons
 
            const double mu    = (NPE-k) * gain  +  (k) * getUnderamplifiedGain();
            const double sigma = gainSpread * sqrt(mu);
 
            if (sigma > 0.0) {
 
              cumulP += weight * (0.5 * erfc((zmin - mu) / sqrt(2.0) / sigma) -
                                  0.5 * erfc((zmax - mu) / sqrt(2.0) / sigma));
              norm   += weight * (0.5 * erfc((normTh   - mu) / sqrt(2.0) / sigma));
 
            } else if ((mu >= zmin) && (mu <= zmax)) {
 
              cumulP += weight;
              norm   += weight;
            
            } else {
 
              return 0.0;
            }
 
            weight *= ((double) (NPE-k) / ((double) (k+1)))  *  PunderAmplified / (1.0 - PunderAmplified);
          }
 
          return cumulP / norm;
        
        } else {
 
          const double mu    = getUnderamplifiedGain(NPE);
          const double sigma = getUnderamplifiedGainSpread(NPE);
 
          if (sigma > 0.0) {
 
            const double cumulP = (0.5 * erfc((zmin - mu) / sqrt(2.0) / sigma) -
                                   0.5 * erfc((zmax - mu) / sqrt(2.0) / sigma));
            const double norm   =  0.5 * erfc((normTh   - mu) / sqrt(2.0) / sigma);
 
            return cumulP / norm;
          
          } else if ((mu >= zmin) && (mu <= zmax)) {
 
            return 1.0;
          
          } else {
 
            return 0.0;
          }
        }
 
      } else {
 
        return 0.0;
      }
    }
 
 
    /**
     * Get integral of probability.
     * The maximum number of photo-electrons is infinity.
     * The normalisation is from (threshold-thresholdband) to infinity.
     *
     * \param  xmin           minimum number of photo-electrons
     * \param  NPE            true    number of photo-electrons
     * \return                probability
     */
    double getIntegralOfChargeProbability(const double xmin, const int NPE) const
    {
      if (NPE > 0) {
 
        double zmin = xmin;
 
        const double th = threshold - thresholdBand;
 
        if (zmin < th) { zmin = th; }
 
        if        (PunderAmplified <= 0.0) {
 
          const double mu    = NPE * gain;
          const double sigma = getGainSpread(NPE);
 
          if (sigma > 0.0) {
 
            const double cumulP =  0.5 * erfc((zmin - mu) / sqrt(2.0) / sigma);
            const double norm   =  0.5 * erfc((th   - mu) / sqrt(2.0) / sigma);
 
            return cumulP / norm;
 
          } else if (mu >= zmin) {
 
            return 1.0;
          
          } else {
 
            return 0.0;
          }
 
        } else if (PunderAmplified < 1.0) {
 
          double cumulP = 0.0;
          double norm   = 0.0;
          double weight = pow(1.0 - PunderAmplified, NPE);
 
          for (int k = 0; k <= NPE; ++k) {  // k underamplified photo-electrons
 
            const double mu    = (NPE-k) * gain  +  (k) * getUnderamplifiedGain();
            const double sigma = gainSpread * sqrt(mu);
 
            if (sigma > 0.0) {
 
              cumulP += weight * (0.5 * erfc((zmin - mu) / sqrt(2.0) / sigma));
              norm   += weight * (0.5 * erfc((th   - mu) / sqrt(2.0) / sigma));
 
            } else if (mu >= zmin) {
 
              cumulP += weight;
              norm   += weight;
            
            } else {
 
              return 0.0;
            }
 
            weight *= ((double) (NPE-k) / ((double) (k+1)))  *  PunderAmplified / (1.0 - PunderAmplified);
          }
 
          return cumulP / norm;
        
        } else {
 
          const double mu    = getUnderamplifiedGain(NPE);
          const double sigma = getGainSpread(NPE);
 
          if (sigma > 0.0) {
 
            const double cumulP =  0.5 * erfc((zmin - mu) / sqrt(2.0) / sigma);
            const double norm   =  0.5 * erfc((th   - mu) / sqrt(2.0) / sigma);
 
            return cumulP / norm;
 
          } else if (mu >= zmin) {
 
            return 1.0;
          
          } else {
 
            return 0.0;
          }
        }
 
      } else {
 
        return 0.0;
      }
    }
 
    
    /**
     * Get integral of probability in specific threshold domain.
     * The normalisation is from (threshold-thresholdband) to infinity.
     *
     * \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 getIntegralOfChargeProbability(threshold, NPE);
 
      case THRESHOLDBAND:
        return getIntegralOfChargeProbability(threshold-thresholdBand, threshold, NPE);
 
      default:
        return 0.0;
      }
    }
 
 
    /**
     * Apply threshold.
     *
     * \param  npe            number of photo-electrons
     * \return                true if pass; else false
     */
    virtual bool applyThreshold(const double npe) const override
    {
      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 tot_ns = getTimeOverThreshold(npe);
 
        if (tot_ns < saturation) {
 
          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
          }
 
          return y / getDerivativeOfSaturation(removeSaturation(tot_ns));
        
        } else {
 
          return 0.0;
        }
      }
 
      default:
        return 0.0;
      }
    }
 
 
    /**
     * Probability that a hit survives the simulation of the PMT.
     * The survival probability takes into account the number of photo-electrons, the analogue signal of the PMT and the threshold of the discriminator.
     *
     * \param  NPE            number of photo-electrons
     * \return                probability
     */
    virtual double getSurvivalProbability(const int NPE) const override
    {
      if (NPE > 0) {
 
        if        (QE <= 0.0) {
 
          return 0.0;
 
        } else if (QE <  1.0) {
 
          double P = 0.0;
    
          for (int i = 1; i <= NPE; ++i) {                                             // i corresponds to number of photo-electrons passing relative QE
 
            const double p = JMATH::binomial(NPE, i) * pow(QE, i) * pow(1.0 - QE, NPE - i);
 
            P += p * getThresholdProbability(NPE);
          }
 
          return P;
 
        } else {
 
          return getThresholdProbability(NPE);
        }
 
      } else {
 
        return 0.0;
      }
    }
 
 
    /**
     * Probability that a hit survives the threshold of the PMT.
     *
     * \param  NPE            number of photo-electrons
     * \return                probability
     */
    virtual double getThresholdProbability(const int NPE) const override
    {
      if (NPE > 0) {
 
        if        (PunderAmplified <= 0.0) {
 
          const double mu    = NPE * gain;
          const double sigma = getGainSpread(NPE);
 
          if (sigma > 0.0) {
 
            const double Ptotal = 0.5 * erfc((0.0                       - mu) / (sqrt(2.0) * sigma));
            const double Pabove = 0.5 * erfc((threshold - thresholdBand - mu) / (sqrt(2.0) * sigma));
 
            return Pabove / Ptotal;
 
          } else if (mu >= threshold - thresholdBand) {
 
            return 1.0;
 
          } else {
 
            return 0.0;
          }
 
        } else if (PunderAmplified < 1.0) {
 
          double Ptotal = 0.0;
          double Pabove = 0.0;
          double weight = pow(1.0 - PunderAmplified, NPE);                         // binomial weight - start with all photo-electrons nominal
 
          for (int k = 0; k <= NPE; ++k) {                                         // NPE-k nominal photo-electrons and k underamplified photo-electrons
 
            const double mu    = (NPE-k) * gain  +  (k) * getUnderamplifiedGain();
            const double sigma = gainSpread * sqrt(mu);
 
            if (sigma > 0.0) {
 
              Ptotal += weight * 0.5 * erfc((0.0                       - mu) / (sqrt(2.0) * sigma));
              Pabove += weight * 0.5 * erfc((threshold - thresholdBand - mu) / (sqrt(2.0) * sigma));
 
            } else if (mu >= threshold - thresholdBand) {
 
              Ptotal += weight;
              Pabove += weight;
            }
 
            weight *= ((double) (NPE-k) / ((double) (k+1)))  *  PunderAmplified / (1.0 - PunderAmplified);
          }
 
          return Pabove / Ptotal;
 
        } else {
 
          const double mu    = getUnderamplifiedGain(NPE);
          const double sigma = getUnderamplifiedGainSpread(NPE);
 
          if (sigma > 0.0) {
 
            const double Ptotal = 0.5 * erfc((0.0                       - mu) / (sqrt(2.0) * sigma));
            const double Pabove = 0.5 * erfc((threshold - thresholdBand - mu) / (sqrt(2.0) * sigma));
 
            return Pabove / Ptotal;
 
          } else if (mu >= threshold - thresholdBand) {
 
            return 1.0;
 
          } else {
 
            return 0.0;
          }
        }
 
      } else {
 
        return 0.0;
      }
    }
 
 
    /**
     * Get number of photo-electrons.
     *
     * \param  tot_ns       time-over-threshold with saturation [ns]
     * \return              number of photo-electrons
     */
    virtual double getNPE(const double tot_ns) const override
    {
      if (tot_ns >= saturation) {
        return std::numeric_limits<double>::max();
      }
 
      const double 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
    {
      if (tot_ns < saturation) {
 
        const double PaboveTh = getIntegralOfChargeProbability(threshold, NPE);
        const double npe      = getNPE(tot_ns);
        const double y        = getChargeProbability(npe, NPE, threshold);
        const double v        = getDerivative(npe);
        const double RaboveTh = PaboveTh * y * v;
        
        const double PthBand  = getIntegralOfChargeProbability(THRESHOLDBAND, NPE);
        const double RthBand  = PthBand * getNormalisedGauss(tot_ns, mean_ns, sigma_ns, 0, saturation);
 
        return RaboveTh + RthBand;
 
      } else {
 
        return 0.0;
      }
    }
 
 
    /**
     * 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
    {
      double Xmin = Tmin;
      double Xmax = Tmax;
 
      if (Tmin > saturation) {Xmin = saturation;}
      if (Tmax > saturation) {Xmax = saturation;}
 
      const double PaboveTh = getIntegralOfChargeProbability(ABOVE_THRESHOLD, NPE);
      const double IaboveTh = PaboveTh * getIntegralOfChargeProbability(getNPE(Xmin), getNPE(Xmax), NPE, threshold);
 
      const double PthBand  = getIntegralOfChargeProbability(THRESHOLDBAND,   NPE);
      const double IthBand  = PthBand * getIntegralOfNormalisedGauss(Xmin, Xmax, mean_ns, sigma_ns, 0, saturation);
 
      return IaboveTh + IthBand;
    }
 
 
    /**
     * 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;
 
  private:
 
    /**
     * Gauss function normalised to a configurable domain.
     *
     * \param  x              x
     * \param  x0             central value
     * \param  sigma          sigma
     * \param  normXmin       normalisation lower limit
     * \param  normXmax       normalisation upper limit
     * \return                probability
     */
    double getNormalisedGauss(const double x, const double x0, const double sigma, const double normXmin, const double normXmax) const
    {
      const double prob = JMATH::Gauss(x, x0, sigma);
      const double norm = (0.5 * erfc((normXmin - x0) / sqrt(2.0) / sigma) -
                           0.5 * erfc((normXmax - x0) / sqrt(2.0) / sigma));
      return prob / norm;
    }
 
 
    /**
     * Integral of Gauss function normalised to a configurable domain.
     *
     * \param  xmin           integral lower limit
     * \param  xmax           integral upper limit
     * \param  x0             central value
     * \param  sigma          sigma
     * \param  normXmin       normalisation lower limit
     * \param  normXmax       normalisation upper limit
     * \return                probability
     */
    double getIntegralOfNormalisedGauss(const double xmin, const double xmax, const double x0, const double sigma, const double normXmin, const double normXmax) const
    {
      const double numerator   = (0.5 * erfc((xmin     - x0) / sqrt(2.0) / sigma_ns) -
                                  0.5 * erfc((xmax     - x0) / sqrt(2.0) / sigma_ns));
      const double denominator = (0.5 * erfc((normXmin - x0) / sqrt(2.0) / sigma_ns) - 
                                  0.5 * erfc((normXmax - x0) / sqrt(2.0) / sigma_ns));
      return numerator / denominator;
    }
  };
 
 
  /**
   * Get charge probability convoluted with a Poisson distribution.
   *
   * \param  pmt             PMT signal processor
   * \param  npe             measured number of photo-electrons
   * \param  NPE             expected number of photo-electrons
   * \param  precision       precision
   * \return                 probability
   */
  inline double getChargeProbability(const JPMTSignalProcessorInterface& pmt,
                                     const double npe,
                                     const double NPE,
                                     const double precision = 1.0e-4)
  {
    int i = (int) (NPE - 5.0 * sqrt(NPE));
 
    if (i < 1) {
      i = 1;
    }
 
    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, p0 = p, p *= NPE / (double) i) {
      P += pmt.getChargeProbability(npe, i) * p;
    }
 
    return P;
  }
 
  
  /**
   * Get time-over-threshold probability convoluted with a Poisson distribution.
   *
   * \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