JPMTAnalogueSignalProcessor.hh

Go to the documentation of this file.

#ifndef __JDETECTOR__JPMTANALOGUESIGNALPROCESSOR__
#define __JDETECTOR__JPMTANALOGUESIGNALPROCESSOR__

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

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

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

  using JLANG::JValueOutOfRange;


  /**
   * Specification for time-over-threshold corresponding to a one photo-electron pulse.
   */
  const double TIME_OVER_THRESHOLD_NS = 26.4;   // [ns]


  /**
   * 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.\n
   *  - If the specified transition time spread (TTS) is negative, 
   *    the transition times are generated according to the measurements (see method JDETECTOR::getTransitionTime).\n
   *  - If the specified TTS is positive,
   *    the transition times are generated according a Gaussian with a sigma equals to the TTS.\n
   *  - If the TTS is zero, the transition times are generated without any spread.
   */
  struct JPMTAnalogueSignalProcessor :
    public JPMTSignalProcessorInterface,
    public JPMTParameters
  {
    /**
     * 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-3;

      // 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 v = slope / getDerivativeOfSaturation(getTimeOverThreshold(1.0));

      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) * v;

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

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

      x1 = 0.5 * (xmin + xmax);
    }


    /**
     * 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 constraint 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 function value with saturation.
     *
     * \param  y              value
     * \return                value
     */
    double getSaturation(const double y) const
    {
      return saturation / sqrt(y*y + saturation*saturation);
    }


    /**
     * Get derivative of function value with saturation.
     *
     * \param  y              value
     * \return                value
     */
    double getDerivativeOfSaturation(const double y) const
    {
      return (saturation*saturation + y*y) * sqrt(saturation*saturation + y*y) / (saturation * saturation * saturation);
    }


    /**
     * 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 getIntegralOfProbability(const double xmin, const double xmax, const int NPE) const
    {
      double zmin = xmin;
      double zmax = xmax;

      if (zmin < threshold) { zmin = threshold; }
      if (zmax < threshold) { zmax = threshold; }
        
      const double mu    = NPE * gain;
      const double sigma = getGainSpread(NPE);

      return (0.5 * TMath::Erfc((zmin - mu) / sqrt(2.0) / sigma)  -
              0.5 * TMath::Erfc((zmax - mu) / sqrt(2.0) / sigma));
    }


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


    /**
     * Apply relative QE.
     *
     * \return                true if accepted; false if rejected
     */
    virtual bool applyQE() const
    {
      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
    {
      if (TTS_ns < 0.0)
        return t_ns + getTransitionTime(gRandom->Rndm());
      else
        return gRandom->Gaus(t_ns, TTS_ns);
    }
    

    /**
     * Compare (arrival times of) photo-electrons.
     *
     * \param  first          first  photo-electron
     * \param  second         second photo-electron
     * \return                true if arrival times of photo-electrons are within rise time of analogue pulses; else false
     */
    virtual bool compare(const JPhotoElectron& first, const JPhotoElectron& second) const
    {
      return second.t_ns < first.t_ns + riseTime_ns;
    }


    /**
     * Get randomised amplitude according gain and gain spread.
     *
     * \param  NPE            number of photo-electrons
     * \return                number of photo-electrons
     */
    virtual double getRandomAmplitude(const int NPE) const
    {
      double x;

      do {
        x = gRandom->Gaus(NPE*gain, getGainSpread(NPE));
      } while (x < 0.0);

      return x;
    }


    /**
     * Get probability for given charge.
     *
     * \param  npe            observed number of photo-electrons
     * \param  NPE            true     number of photo-electrons
     * \return                probability [npe^-1]
     */
    virtual double getProbability(const double npe, const int NPE) const
    {
      if (npe >= threshold) {
        
        const double mu    = NPE * gain;
        const double sigma = getGainSpread(NPE);
        const double V     = 0.5 * TMath::Erfc((threshold - mu) / sqrt(2.0) / sigma);
        
        return TMath::Gaus(npe, mu, sigma, kTRUE) / V;
      }

      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 npe >= 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
    {
      if (slewing)
        return ((getRiseTime(npe, getTH0())  -  getRiseTime(npe, this->threshold)) -
                (getRiseTime(1.0, getTH0())  -  getRiseTime(1.0, this->threshold)));
      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
    {
      if (npe > threshold) {
        
        double tot = 0.0;

        if (npe*y1 > threshold) {

          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
          }

        } else {

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

        return tot * getSaturation(tot);

      } else {

        return 0.0;
      }
    }


    /**
     * 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
    {
      if (npe >= 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
        }
        
        return y * getDerivativeOfSaturation(getTimeOverThreshold(npe));
          
      } else {
        
        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
    {
      if        (QE <= 0.0) {

        return 0.0;

      } else if (QE <= 1.0) {

        double P = 0.0;

        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);
          const double mu    = i * gain;
          const double sigma = getGainSpread(i);
          const double V     = 0.5 * TMath::Erfc((threshold  - mu) / sqrt(2.0) / sigma);
          const double W     = 0.5 * TMath::Erfc((0.0        - mu) / sqrt(2.0) / sigma);

          P += p*V/W;
        }
        
        return P;

      } else {

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


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

#endif