JPMTSignalProcessorInterface.hh

Go to the documentation of this file.

#ifndef __JDETECTOR__JPMTSIGNALPROCESSORINTERFACE__
#define __JDETECTOR__JPMTSIGNALPROCESSORINTERFACE__
 
#include <algorithm>
#include <cmath>
#include <limits>
 
#include "JDetector/JPMTSimulator.hh"
#include "JDetector/JCalibration.hh"
 
 
/**
 * \author mdejong
 */
 
namespace JDETECTOR {}
namespace JPP { using namespace JDETECTOR; }
 
namespace JDETECTOR {
 
  /**
   * PMT signal processor interface.
   *
   * This class supports the implementation of the PMT simulator interface using 
   * an alternative set of virtual methods.
   * These methods constitute a user interface to the expected performance of a PMT.
   *
   * Each photon is converted to a photo-electron using member method getRandomTime().
   * For this, the data structure JPhotoElectron is used.
   * Note that the quantum efficiency (QE) of the PMT actually is included in the simulation of the detector response.
   * A relative QE is applied via method applyQE().
   * All photo-electrons are then time sorted.
   * The photo-electrons which simultaneously arrive are merged.
   * The corresponding condition is defined by member method compare().
   * The time of the combined signal is determined by the time of the first photo-electron and
   * the rise time of the analogue pulse to the threshold of the discriminator via method getRiseTime().
   * In this, the actual amplitude of the combined analogue signal and the calibration of the PMT are taken into account.
   * The amplitude of the combined analogue signal is simulated using member method getRandomCharge().
   * For this, the data structure JPMTSignal is used.
   *
   * The analogue signals of the PMT are processed by a discriminator.
   * This discriminator produces a time-over-threshold signal for each analogue signal that passes a preset threshold.
   * The output signal is described by the time of the leading edge and the length of the time-over-threshold signal.
   * For this, the data structure JPMTPulse is used.
   * The determination of the time of the leading edge and the length of the time-over-threshold signal
   * require a designated model.
   * The class JPMTAnalogueSignalProcessor provides for an implementation of such a model.
   *
   * Overlapping time-over-threshold pulses are merged.
   * The time of the leading edge is then set to the earliest leading edge and 
   * the time-over-threshold to the difference between
   * the latest trailing edge and the earliest leading edge.
   * The merging of hits is implemented in member method merge().
   *
   * The default implementation of these virtual methods corresponds to no smearing 
   * and a time-over-threshold value equal to a fixed two photo-electron resolution times the number of photo-electrons.
   * The width of the charge distribution and the two photo-electron resolution are implemented in methods
   * getQmin() and getTmin(), respectively.
   *
   * For a realistic PMT simulation, the derived class should provide for 
   * an implementation of each of these virtual methods.
   */
  class JPMTSignalProcessorInterface {
  public:
    /**
     * Default constructor.
     */
    JPMTSignalProcessorInterface()
    {}
 
 
    /**
     * Virtual destructor.
     */
    virtual ~JPMTSignalProcessorInterface()
    {}
 
 
    /**
     * Process hits.
     *
     * Two (or more) photo-electrons are combined if they are comparable according method compare.\n
     * Two (or more) consecutive hits hits maybe merged (according method merge).\n
     * A PMT signal with negative number of photo-electrons is not subjected to the quantum efficiency and transit-time distribution.\n
     * This feature is used to simulate the singles rate of a PMT when the measured QE is zero. 
     *
     * \param  calibration    PMT calibration
     * \param  input          PMT signals
     * \param  output         PMT hits
     */
    void operator()(const JCalibration&         calibration,
                    const JPMTData<JPMTSignal>& input,
                    JPMTData<JPMTPulse>&        output) const
    {
      // apply transition time distribution to each photo-electron.
 
      JPMTData<JPhotoElectron> buffer;
      
      for (JPMTData<JPMTSignal>::const_iterator hit = input.begin(); hit != input.end(); ++hit) {
 
        for (int i = 0; i < hit->npe; ++i) {
          if (applyQE()) {
            buffer.push_back(JPhotoElectron(getRandomTime(hit->t_ns)));
          }
        }
 
        for (int i = 0; i > hit->npe; --i) {
          buffer.push_back(JPhotoElectron(hit->t_ns));
        }
      }
 
      if (!buffer.empty()) {
      
        buffer.push_back(JPhotoElectron::getEndMarker());
      
        buffer.sort();
      
      
        // generate PMT hits from time sequence of photo-electrons.
        
        for (JPMTData<JPhotoElectron>::const_iterator q = buffer.begin(), p = q++; q != buffer.end(); ++q) {
          
          while (compare(*p,*q)) {
            ++q;
          }
          
          const double npe = getRandomCharge(distance(p,q));
 
          if (applyThreshold(npe)) {
            output.push_back(JPMTPulse(putTime(p->t_ns + getRiseTime(npe), calibration), getTimeOverThreshold(npe)));
          }
 
          p = q;
        }
        
        // merge overlapping PMT hits.
      
        merge(output);
      }
    }
 
 
    /**
     * Apply relative QE.
     * The default implementation returns <tt>true</tt>.
     *
     * \return                true if accepted; false if rejected
     */
    virtual bool applyQE() const 
    {
      return true;
    }
 
 
    /**
     * Get randomised time according transition time distribution.
     * The default implementation returns the given value.
     *
     * \param  t_ns           time [ns]
     * \return                time [ns]
     */
    virtual double getRandomTime(const double t_ns) const
    {
      return t_ns;
    }
 
                        
    /**
     * Compare arrival times of photo-electrons.
     * The default implementation uses method getTmin as two photo-electron resolution.
     *
     * \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 
    {
      return second.t_ns - first.t_ns <= getTmin();
    }
 
 
    /**
     * Get randomised charge according to gain and gain spread.
     * The default implementation returns the given value.
     *
     * \param  NPE            number of photo-electrons
     * \return                number of photo-electrons
     */
    virtual double getRandomCharge(const int NPE) const 
    {
      return (double) NPE;
    }
 
 
    /**
     * Get probability density for given charge.
     *
     * \param  npe            observed number of photo-electrons
     * \param  NPE            true     number of photo-electrons
     * \return                probability [npe^-1]
     */
    virtual double getChargeProbability(const double npe, const int NPE) const 
    {
      return (fabs(npe - NPE) <= 0.5 * getQmin() ? 1.0 / getQmin() : 0.0);
    }
 
 
    /**
     * Apply threshold.
     * The default implementation returns <tt>true</tt>.
     *
     * \param  npe            number of photo-electrons
     * \return                true if pass; else false
     */
    virtual bool applyThreshold(const double npe) const 
    {
      return (npe > 0.0);
    }
 
 
    /**
     * Get time to reach threshold.
     * The default implementation returns <tt>0</tt>.
     *
     * \param  npe            number of photo-electrons
     * \return                time [ns]
     */
    virtual double getRiseTime(const double npe) const 
    {
      return 0.0;
    }
 
 
    /**
     * Get time-over-threshold (ToT).
     * The default implementation corresponds to a linear relation between the number of photo-electrons and the time-over-threshold.
     *
     * \param  npe            number of photo-electrons
     * \return                ToT [ns]
     */
    virtual double getTimeOverThreshold(const double npe) const 
    {
      return TIME_OVER_THRESHOLD_NS  +  getTmin() * (round(npe) - 1.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 
    {
      return 1.0 / getTmin();
    }
 
 
    /**
     * Probability that a hit survives the simulation of the PMT.
     * The default implementation returns probability to pass the PMT threshold.
     *
     * \param  NPE            number of photo-electrons
     * \return                probability
     */
    virtual double getSurvivalProbability(const int NPE) const
    {
      return getThresholdProbability(NPE);
    }
 
 
    /**
     * Probability that a hit survives the threshold of the PMT.
     * The default implementation returns one (zero) if given number of photo-electrons is larger than (less or equal to) zero.
     *
     * \param  NPE            number of photo-electrons
     * \return                probability
     */
    virtual double getThresholdProbability(const int NPE) const
    {
      if (NPE > 0)
        return 1.0;
      else
        return 0.0;
    }
 
 
    /**
     * Get number of photo-electrons.
     * The default implementation corresponds to a linear relation between the number of photo-electrons and the time-over-threshold.
     *
     * \param  tot_ns       time over threshold [ns]
     * \return              number of photo-electrons
     */
    virtual double getNPE(const double tot_ns) const
    {
      return 1.0  +  (tot_ns - TIME_OVER_THRESHOLD_NS) / getTmin();
    }
 
 
    /**
     * Merging of PMT hits.
     *
     * Hits with overlapping time-over-threshold signals should -de facto- be combined.
     * In this, the leading edge is maintained and the time-over-threshold is
     * set to the difference between the overall trailing and leading edges.
     * As a result, the number of PMT hits may be reduced.
     *
     * \param  data           PMT hits (I/O)
     */
    virtual void merge(JPMTData<JPMTPulse>& data) const
    {
      using namespace std;
      
      JPMTData<JPMTPulse>::iterator out = data.begin();
      
      for (JPMTData<JPMTPulse>::iterator i = data.begin(); i != data.end(); ) {
        
        double t1 = i->t_ns;
        double t2 = i->t_ns + i->tot_ns;
        
        while (++i != data.end() && i->t_ns < t2 + getTmin()) {
          t2 = max(t2, i->t_ns + i->tot_ns);
        }
        
        out->t_ns   = t1;
        out->tot_ns = t2 - t1;
        
        ++out;
      }
      
      data.resize(distance(data.begin(), out));
    }
 
 
    /**
     * Get two photo-electron resolution for time-over-threshold
     *
     * \return              minimal time [ns]
     */
    static double getTmin()
    {
      return 1.0;
    }
 
 
    /**
     * Get width of charge distribution. 
     *
     * \return              width charge distribution [npe]
     */
    static double getQmin()
    {
      return 1.0e-3;
    }
  };
 
 
  /**
   * Get charge probability.
   *
   * \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;
  }
}
 
#endif