JPMTSignalProcessorInterface.hh

Go to the documentation of this file.

#ifndef __JDETECTOR__JPMTSIGNALPROCESSORINTERFACE__
#define __JDETECTOR__JPMTSIGNALPROCESSORINTERFACE__

#include <algorithm>
#include <cmath>
#include <limits>

#include "JLang/JCC.hh"
#include "JLang/JException.hh"
#include "JDetector/JPMTSimulator.hh"
#include "JDetector/JPMTTransitTimeGenerator.hh"


/**
 * \author mdejong
 */

namespace JDETECTOR {}
namespace JPP { using namespace JDETECTOR; }

namespace JDETECTOR {


  using JLANG::JValueOutOfRange;


  /**
   * 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 JPMTSignalProcessorInterface::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.
   * All photo-electrons are then time sorted.
   * The photo-electrons which simultaneously arrive are merged.
   * The corresponding condition is defined by member method JPMTSignalProcessorInterface::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.
   * 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 JPMTSignalProcessorInterface::getRandomAmplitude.
   * 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 JPMTSignalProcessorInterface::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.
   * The width of the charge distribution and the two photo-electron resolution are implemented in methods
   * JPMTSignalProcessorInterface::getQmin() and JPMTSignalProcessorInterface::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.
     *
     * \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.

      buffer.clear();
      
      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)));
          }
        }
      }

      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 = getRandomAmplitude(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 always <tt>true</tt>.
     *
     * \return                true
     */
    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.
     *
     * Two (or more) photo-electrons are merged if they are comparable.
     * The default implementation returns uses JPMTSignalProcessorInterface::getTmin() as time window.
     *
     * \param  first          first  photo-electron
     * \param  second         second photo-electron
     * \return                false
     */
    virtual bool compare(const JPhotoElectron& first, const JPhotoElectron& second) const 
    {
      return second.t_ns - first.t_ns <= getTmin();
    }


    /**
     * Get randomised amplitude according gain and gain spread.
     * The default implementation returns the given value.
     *
     * \param  NPE            number of photo-electrons
     * \return                number of photo-electrons
     */
    virtual double getRandomAmplitude(const int NPE) const 
    {
      return (double) NPE;
    }


    /**
     * 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 
    {
      return (fabs(npe - NPE) < 0.5*getQmin() ? 1.0/getQmin() : 0.0);
    }


    /**
     * Apply threshold.
     * The default implementation returns always <tt>true</tt>.
     *
     * \param  npe            number of photo-electrons
     * \return                true if pass; else false
     */
    virtual bool applyThreshold(const double npe) const 
    {
      return true;
    }


    /**
     * Get time to pass 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 returns JPMTSignalProcessorInterface::getTmin().
     *
     * \param  npe            number of photo-electrons
     * \return                ToT [ns]
     */
    virtual double getTimeOverThreshold(const double npe) const 
    {
      return getTmin();
    }


    /**
     * Get derivative of number of photo-electrons to time over threshold.
     * The default implementation returns infinity.
     *
     * \param  npe            number of photo-electrons
     * \return                dnpe/dToT [ns^-1]
     */
    virtual double getDerivative(const double npe) const 
    {
      return std::numeric_limits<double>::max();
    }


    /**
     * Probability that a hit survives the simulation of the PMT.
     * The default implementation returns <tt>1</tt> if given value larger than <tt>0</tt>.
     *
     * \param  NPE            number of photo-electrons
     * \return                probability
     */
    virtual double getSurvivalProbability(const int NPE) const
    {
      if (NPE > 0)
        return 1.0;
      else
        return 0.0;
    }


    /**
     * 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
     */
    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 number of photo-electrons.
     *
     * \param  tot_ns       time over threshold [ns]
     * \param  eps          precision
     * \return              number of photo-electrons
     */
    inline double getNPE(const double tot_ns,
                         const double eps = 1.0e-3) const
    {
      static const int N = 100;

      if (tot_ns > 0.0) {

        double xmin =   0.0;
        double xmax = 100.0;
      
        for (int i = 0; i != N; ++i) {
        
          const double x = 0.5 * (xmin + xmax);
          const double y = this->getTimeOverThreshold(x);
          
          if (fabs(y - tot_ns) < eps) {
            return x;
          }
          
          if (y < tot_ns)
            xmin = x;
          else
            xmax = x;
        }

        return 0.5 * (xmin + xmax);

      } else {

        double xmin = 0.0;
        double xmax = 2.0;
      
        for (int i = 0; i != N; ++i) {
          
          const double x = 0.5 * (xmin + xmax);
          const double y = this->getTimeOverThreshold(x);

          if (y == 0.0)
            xmin = x;
          else
            xmax = x;
        }

        return 0.5 * (xmin + xmax);
      }
    }


    /**
     * Get two photo-electron resolution for time-over-threshold
     *
     * \return              minimal time [ns]
     */
    static double getTmin()
    {
      return 1.0;
    }


    /**
     * Get minimal width of charge distribution. 
     *
     * \return              minimal charge [npe]
     */
    static double getQmin()
    {
      return 1.0e-3;
    }


  private:
    mutable JPMTData<JPhotoElectron> buffer;
  };
}

#endif