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/JGauss.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
    };
 
 
    /**
     * Amplification mode specifiers.
     */
    enum class JAmplificationMode {
      Nominal,                  // no underamplification
      Mixed,                    // underamplification
      UnderAmplified            // complete underamplification
    };
 
 
    /**
     * Analogue pulse domain specifiers.
     */
    enum class JAnaloguePulseDomain {
      Gaussian,                 // Gaussian
      Exponential,              // Gaussian + exponential
      Linear                    // Gaussian + exponential + linear
    };
 
 
    /**
     * 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 getRiseTime(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.
     * 
     * \param NPE             number of photo-electrons
     * \return                gain
     */
    double getGain(int NPE) const
    {
      return ((double) NPE) * gain;
    }
 
 
    /**
     * 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;
      }
    }
 
 
    /**
     * Get amplification mode.
     * 
     * \return                amplification mode  
     */
    JAmplificationMode getAmplificationMode() const
    {
      if (PunderAmplified <= 0.0) {
 
        return JAmplificationMode::Nominal;
 
      } else if (PunderAmplified < 1.0) {
 
        return JAmplificationMode::Mixed;
 
      } else {
 
        return JAmplificationMode::UnderAmplified;
      }
    }
 
 
    /**
     * Get charge domain.
     * 
     * \param  npe            observed number of photo-electrons
     */
    JAnaloguePulseDomain getChargeDomain(const double npe) const
    {
      if        (npe*y1 <= threshold) {
 
        return JAnaloguePulseDomain::Gaussian;
      
      } else if (npe <= getStartOfLinearisation()) {
 
        return JAnaloguePulseDomain::Exponential;
      
      } else {
 
        return JAnaloguePulseDomain::Linear;
      }
    }
 
 
    /**
     * Get time-over-threshold domain.
     * 
     * \param  tot_ns         time-over-threshold (with saturation) [ns]
     */
    JAnaloguePulseDomain getToTDomain(const double tot_ns) const
    {
      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 JAnaloguePulseDomain::Gaussian;
        
      } else if (tot <= TOT) {                                    // Gaussian + Exponential
 
        return JAnaloguePulseDomain::Exponential;
        
      } else {                                                                                  // linear
          
        return JAnaloguePulseDomain::Linear;
      }
    }
 
 
    /**
     * 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) {
 
          switch (getAmplificationMode()) {
 
          case JAmplificationMode::Nominal:
            return getRandomChargeNominal(NPE);
 
          case JAmplificationMode::Mixed:
            return getRandomChargeMixed(NPE);
 
          case JAmplificationMode::UnderAmplified:
            return getRandomChargeUnderAmplified(NPE);
 
          default:
            return 0.0;
          }
        
        } else {
 
          return getGain(NPE);
        }
 
      } else {
 
        return 0.0;
      }
    }
 
 
    /**
     * Get randomised charge according to gain and gain spread with a nominal amplification mode.
     *
     * \param  NPE            number of photo-electrons
     * \return                number of photo-electrons
     */
    double getRandomChargeNominal(const int NPE) const
    {
      double q;
 
      do {
 
        const double mu    = getGain(NPE);
        const double sigma = getGainSpread(NPE);
 
        q = gRandom->Gaus(mu, sigma);
 
      } while (q < 0.0);
      
      return q;
    }
 
 
    /**
     * Get randomised charge according to gain and gain spread with a mixed amplification mode.
     *
     * \param  NPE            number of photo-electrons
     * \return                number of photo-electrons
     */
    double getRandomChargeMixed(const int NPE) const
    {
 
      double q;
 
      do {
        // 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);
 
      } while (q < 0.0);
      
      return q;
    }
 
 
    /**
     * Get randomised charge according to gain and gain spread with an underamplified amplification mode.
     *
     * \param  NPE            number of photo-electrons
     * \return                number of photo-electrons
     */
    double getRandomChargeUnderAmplified(const int NPE) const
    {
      double q;
 
      do {
 
        const double mu    = getUnderamplifiedGain(NPE);
        const double sigma = getUnderamplifiedGainSpread(NPE);
 
        q = gRandom->Gaus(mu, sigma);
 
      } while (q < 0.0);
      
      return q;
    }
 
 
    /**
     * 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)) {
 
        switch (getAmplificationMode()) {
 
        case JAmplificationMode::Nominal:
          return getChargeProbabilityNominal(npe, NPE);
 
        case JAmplificationMode::Mixed:
          return getChargeProbabilityMixed(npe, NPE);
 
        case JAmplificationMode::UnderAmplified:
          return getChargeProbabilityUnderAmplified(npe, NPE);
 
        default:
          return 0.0;
        }
 
      } else {
 
        return 0.0;
      }
    }
 
    /**
     * Get probability density for given charge with a Nominal amplification mode.
     * 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]
     */
    double getChargeProbabilityNominal(const double npe, const int NPE) const
    {
      const double mu    = getGain(NPE);
      const double sigma = getGainSpread(NPE);
 
      if (sigma > 0.0) {
 
        JMATH::JGauss gauss = JMATH::JGauss(mu, sigma);
 
        const double prob = gauss.getValue(npe);
 
        const double normLowerLimit = threshold - thresholdBand;
        const double normUpperLimit = std::numeric_limits<double>::infinity();
        const double norm           = gauss.getIntegral(normLowerLimit, normUpperLimit);
 
        return prob / norm;
      
      } else if (std::fabs(npe - mu) <= 0.0001 * mu) {
 
        return 1.0;
      
      } else {
 
        return 0.0;
      }
    }
 
 
    /**
     * Get probability density for given charge with a mixed amplification mode.
     * 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]
     */
    double getChargeProbabilityMixed(const double npe, const double NPE) const
    {
      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) {
 
          JMATH::JGauss gauss = JMATH::JGauss(mu, sigma);
 
          prob += weight * gauss.getValue(npe);
 
          const double normLowerLimit  = threshold - thresholdBand;
          const double normUpperLimit  = std::numeric_limits<double>::infinity();
          
          norm += weight * gauss.getIntegral(normLowerLimit, normUpperLimit);
        
        } else if (std::fabs(npe - mu) <= 0.0001 * mu) {
 
          prob += weight;
          norm += weight;
        
        } else {
 
          return 0.0;
        }
  
        weight *= ((double) (NPE-k) / ((double) (k+1)))  *  PunderAmplified / (1.0 - PunderAmplified);
      }
 
      return prob / norm;
    }
 
 
    /**
     * Get probability density for given charge with an underamplified amplification mode.
     * 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]
     */
    double getChargeProbabilityUnderAmplified(const double npe, const int NPE) const
    {
      if (gainSpread > 0.0) {
 
        const double mu    = getUnderamplifiedGain(NPE);
        const double sigma = getUnderamplifiedGainSpread(NPE);
        
        JMATH::JGauss gauss = JMATH::JGauss(mu, sigma);
 
        const double prob = gauss.getValue(npe);
 
        const double normLowerLimit = threshold - thresholdBand;
        const double normUpperLimit = std::numeric_limits<double>::infinity();
        const double norm           = gauss.getIntegral(normLowerLimit, normUpperLimit);
 
        return prob / norm;
      
      } 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; }
 
        switch (getAmplificationMode()) {
 
        case JAmplificationMode::Nominal:
          return getIntegralOfChargeProbabilityNominal(zmin, zmax, NPE);
 
        case JAmplificationMode::Mixed:
          return getIntegralOfChargeProbabilityMixed(zmin, zmax, NPE);
 
        case JAmplificationMode::UnderAmplified:
          return getIntegralOfChargeProbabilityUnderAmplified(zmin, zmax, NPE);
 
        default:
          return 0.0;
        }
 
      } else {
 
        return 0.0;
      }
    }
 
 
    /**
     * Get integral of probability with a nominal amplification mode.
     * 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 getIntegralOfChargeProbabilityNominal(const double xmin, const double xmax, const int NPE) const
    {
      const double mu    = getGain(NPE);
      const double sigma = getGainSpread(NPE);
 
      if (sigma > 0.0) {
 
        JMATH::JGauss gauss = JMATH::JGauss(mu, sigma);
        
        const double cumulP = gauss.getIntegral(xmin, xmax);
 
        const double normLowerLimit = threshold - thresholdBand;
        const double normUpperLimit = std::numeric_limits<double>::infinity();
        const double norm           = gauss.getIntegral(normLowerLimit, normUpperLimit);
 
        return cumulP / norm;
 
      } else if ((mu >= xmin) && (mu <= xmax)) {
 
        return 1.0;
      
      } else {
 
        return 0.0;
      }
    }
 
 
    /**
     * Get integral of probability with a mixed amplification mode.
     * 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 getIntegralOfChargeProbabilityMixed(const double xmin, const double xmax, const int NPE) const
    {
      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) {
 
          JMATH::JGauss gauss = JMATH::JGauss(mu, sigma);
 
          const double normLowerLimit = threshold - thresholdBand;
          const double normUpperLimit = std::numeric_limits<double>::infinity();
 
          cumulP += weight * gauss.getIntegral(xmin, xmax);
          norm   += weight * gauss.getIntegral(normLowerLimit, normUpperLimit);
        
        } else if ((mu >= xmin) && (mu <= xmax)) {
 
          cumulP += weight;
          norm   += weight;
        
        } else {
 
          return 0.0;
        }
 
        weight *= ((double) (NPE-k) / ((double) (k+1)))  *  PunderAmplified / (1.0 - PunderAmplified);
      }
 
      return cumulP / norm;
    }
 
 
    /**
     * Get integral of probability with an underamplified amplification mode.
     * 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 getIntegralOfChargeProbabilityUnderAmplified(const double xmin, const double xmax, const int NPE) const
    {
      if (gainSpread > 0.0) {
 
        const double mu    = getUnderamplifiedGain(NPE);
        const double sigma = getUnderamplifiedGainSpread(NPE);
 
        JMATH::JGauss gauss = JMATH::JGauss(mu, sigma);
 
        const double cumulP = gauss.getIntegral(xmin, xmax);
 
        const double normLowerLimit = threshold - thresholdBand;
        const double normUpperLimit = std::numeric_limits<double>::infinity();
        const double norm           = gauss.getIntegral(normLowerLimit, normUpperLimit);
 
        return cumulP / norm;
      
      } 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, std::numeric_limits<double>::infinity(), 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.
     *
     * \param  npe            number of photo-electrons
     * \return                time-over-threshold (with saturation) [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;
 
        switch (getChargeDomain(npe)) {
 
        case JAnaloguePulseDomain::Gaussian:
          tot += getRiseTime(npe, threshold);                           // Gaussian
          tot += getRiseTime(npe, threshold);                           // Gaussian
          break;
 
        case JAnaloguePulseDomain::Exponential:
          tot += getRiseTime (npe, threshold);                          // Gaussian
          tot += getDecayTime(npe, threshold);                          // exponential
          break;
 
        case JAnaloguePulseDomain::Linear:
          tot += getRiseTime (getStartOfLinearisation(), threshold);    // Gaussian
          tot += getDecayTime(getStartOfLinearisation(), threshold);    // exponential
    
          tot += slope * (npe - getStartOfLinearisation());                 // linear
          break;
        }
        
        return applySaturation(tot);
      }
        
      default:
        THROW(JValueOutOfRange, "JPMTAnalogueSignalProcessor::getTimeOverThreshold: Invalid charge " << npe);
      }
    }
 
 
    /**
     * Get derivative of number of photo-electrons to time-over-threshold with saturation.
     *
     * \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;
 
          switch (getChargeDomain(npe)) {
 
          case JAnaloguePulseDomain::Gaussian:
            y += npe / (2.0 * z);
            break;
 
          case JAnaloguePulseDomain::Exponential:
            y += npe / (z + decayTime_ns);
            break;
 
          case JAnaloguePulseDomain::Linear:
            y += 1.0 / slope;
            break;
          }
 
          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) {
 
        switch (getAmplificationMode()) {
 
        case JAmplificationMode::Nominal:
          return getThresholdProbabilityNominal(NPE);
 
        case JAmplificationMode::Mixed:
          return getThresholdProbabilityMixed(NPE);
 
        case JAmplificationMode::UnderAmplified:
          return getThresholdProbabilityUnderAmplified(NPE);
 
        default:
          return 0.0;
        }
 
      } else {
 
        return 0.0;
      }
    }
 
 
    /**
     * Probability that a hit survives the threshold of the PMT with a nominal amplification mode.
     *
     * \param  NPE            true    number of photo-electrons
     * \return                probability
     */
    double getThresholdProbabilityNominal(const int NPE) const
    {
      const double mu    = getGain(NPE);
      const double sigma = getGainSpread(NPE);
 
      if (sigma > 0.0) {
 
        JMATH::JGauss gauss = JMATH::JGauss(mu, sigma);
        
        const double probLowerLimit = threshold - thresholdBand;
        const double probUpperLimit = std::numeric_limits<double>::infinity();
        const double cumulP         = gauss.getIntegral(probLowerLimit, probUpperLimit);
 
        const double normLowerLimit = 0.0;
        const double normUpperLimit = std::numeric_limits<double>::infinity();
        const double norm           = gauss.getIntegral(normLowerLimit, normUpperLimit);
 
        return cumulP / norm;
 
      } else if (mu >= threshold - thresholdBand) {
 
        return 1.0;
      
      } else {
 
        return 0.0;
      }
    }
 
 
    /**
     * Probability that a hit survives the threshold of the PMT with a mixed amplification mode.
     *
     * \param  NPE            true    number of photo-electrons
     * \return                probability
     */
    double getThresholdProbabilityMixed(const int NPE) const
    {
      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) {
 
          JMATH::JGauss gauss = JMATH::JGauss(mu, sigma);
 
          const double probLowerLimit = threshold - thresholdBand;
          const double probUpperLimit = std::numeric_limits<double>::infinity();
          
          cumulP += weight * gauss.getIntegral(probLowerLimit, probUpperLimit);
 
          const double normLowerLimit = 0.0;
          const double normUpperLimit = std::numeric_limits<double>::infinity();
 
          norm   += weight * gauss.getIntegral(normLowerLimit, normUpperLimit);
        
        } else if (mu >= threshold - thresholdBand) {
 
          cumulP += weight;
          norm   += weight;
        
        } else {
 
          return 0.0;
        }
 
        weight *= ((double) (NPE-k) / ((double) (k+1)))  *  PunderAmplified / (1.0 - PunderAmplified);
      }
 
      return cumulP / norm;
    }
 
 
    /**
     * Probability that a hit survives the threshold of the PMT with an underamplified amplification mode.
     *
     * \param  NPE            true    number of photo-electrons
     * \return                probability
     */
    double getThresholdProbabilityUnderAmplified(const int NPE) const
    {
      if (gainSpread > 0.0) {
 
        const double mu    = getUnderamplifiedGain(NPE);
        const double sigma = getUnderamplifiedGainSpread(NPE);
 
        JMATH::JGauss gauss = JMATH::JGauss(mu, sigma);
 
        const double probLowerLimit = threshold - thresholdBand;
        const double probUpperLimit = std::numeric_limits<double>::infinity();
        const double cumulP         = gauss.getIntegral(probLowerLimit, probUpperLimit);
 
        const double normLowerLimit = 0.0;
        const double normUpperLimit = std::numeric_limits<double>::infinity();
        const double norm           = gauss.getIntegral(normLowerLimit, normUpperLimit);
 
        return cumulP / norm;
      
      } 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));
      
      switch (getToTDomain(tot_ns)) {
      
      case JAnaloguePulseDomain::Gaussian:
        return threshold * exp(tot*tot/riseTime_ns/riseTime_ns/8.0);
        
      case JAnaloguePulseDomain::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);
      }
        
      case JAnaloguePulseDomain::Linear:
        return getStartOfLinearisation() + (tot - TOT) / slope;
 
      default:
        return 0.0;
      }
    }
 
 
    /**
     * 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 npe      = getNPE(tot_ns);
        const double y        = getChargeProbability(npe, NPE);
        const double v        = getDerivative(npe);
        const double RaboveTh = y * v;
        
        const double w         = getIntegralOfChargeProbability(ABOVE_THRESHOLD, NPE);
        JMATH::JGauss gauss    = JMATH::JGauss(mean_ns, sigma_ns);
        const double gaussProb = gauss.getValue(tot_ns);
        const double gaussNorm = gauss.getIntegral(0, saturation);
        const double gaussDist = gaussProb / gaussNorm;
        const double RthBand   = (1-w) * gaussDist;
 
        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 IaboveTh = getIntegralOfChargeProbability(getNPE(Xmin), getNPE(Xmax), NPE);
 
      const double w           = getIntegralOfChargeProbability(ABOVE_THRESHOLD, NPE);
      JMATH::JGauss gauss      = JMATH::JGauss(mean_ns, sigma_ns);
      const double gaussCumulP = gauss.getIntegral(Xmin, Xmax);
      const double gaussNorm   = gauss.getIntegral(0, saturation);
      const double gaussDist   = gaussCumulP / gaussNorm;
      const double IthBand     = (1-w) * gaussDist;
 
      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;
  };
 
 
  /**
   * 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