JK40DefaultSimulatorInterface.hh

Go to the documentation of this file.

#ifndef __JDETECTOR__JK40DEFAULTSIMULATORINTERFACE__
#define __JDETECTOR__JK40DEFAULTSIMULATORINTERFACE__
 
#include "TRandom3.h"
 
#include "JDetector/JK40Simulator.hh"
#include "JDetector/JModuleIdentifier.hh"
#include "JDetector/JPMTIdentifier.hh"
#include "JDetector/JTimeRange.hh"
 
#include "JMath/JMathToolkit.hh"
#include "JLang/JException.hh"
#include "JLang/JCC.hh"
 
 
/**
 * \author mdejong
 */
 
namespace JDETECTOR {}
namespace JPP { using namespace JDETECTOR; }
 
namespace JDETECTOR {
 
  using JLANG::JNumericalPrecision;
 
 
  /**
   * Default K40 simulator interface.
   *
   * This class provides for a default implementation of the JK40Simulator interface
   * which is based on a set of virtual methods.
   * These methods constitute a user interface to the K40 background simulation.
   */
  class JK40DefaultSimulatorInterface :
    public JK40Simulator
  {
  protected:
    /**
     * 1D-linear buffer.
     */
    struct JBuffer1D_t : 
      public std::vector<double>
    {};
 
 
    /**
     * 2D-square buffer.
     */
    struct JBuffer2D_t : 
      public std::vector<JBuffer1D_t>
    {
      /**
       * Resize.
       *
       * \param  size    size
       */
      void resize(size_t size)
      {
        std::vector<JBuffer1D_t>::resize(size);
 
        for (iterator i = begin(); i != end(); ++i) {
          i->resize(size);
        }
      }
    };
 
 
    /**
     * Default constructor.
     */
    JK40DefaultSimulatorInterface()
    {}
 
 
  public:
    /**
     * Get singles rate as a function of PMT.
     *
     * \param  pmt            PMT identifier
     * \return                rate [Hz]
     */
    virtual double getSinglesRate(const JPMTIdentifier& pmt) const = 0;
 
    
    /**
     * Get multiples rate as a function of optical module.
     *
     * \param  module         optical module identifier
     * \param  M              multiplicity (M >= 2)
     * \return                rate [Hz]
     */
    virtual double getMultiplesRate(const JModuleIdentifier& module, const int M) const = 0;
 
 
    /**
     * Get probability of coincidence.
     *
     * \param  ct             cosine space angle between PMT axes
     * \return                probability
     */
    virtual double getProbability(const double ct) const = 0;
 
 
    /**
     * Generate hits.
     *
     * \param  module         module
     * \param  period         time window [ns]
     * \param  output         background data
     */
    virtual void generateHits(const JModule&    module,
                              const JTimeRange& period,
                              JModuleData&      output) const
    {
 
      // resize internal buffers
      
      const int N = module.size();
 
      probability2D.resize(N);
      probability1D.resize(N);
      probabilityND.resize(N);
 
      rateL1_Hz.resize(N);
 
 
      // generate singles
 
      for (int pmt = 0; pmt != N; ++pmt) {
 
        const double rateL0_Hz = getSinglesRate(JPMTIdentifier(module.getID(), pmt));
 
        if (rateL0_Hz > 0.0) {
              
          const double t_ns = 1.0e9 / rateL0_Hz;          // [ns]
              
          for (double t1 = period.getLowerLimit() + gRandom->Exp(t_ns); t1 < period.getUpperLimit(); t1 += gRandom->Exp(t_ns)) {
            output[pmt].push_back(JPMTSignal(t1, 1));
          }
        }
      }
 
 
      // generate coincidences
 
      double totalRateL1_Hz = 0.0;
 
      for (int i = 0; i != N; ++i) {
        totalRateL1_Hz += rateL1_Hz[i] = getMultiplesRate(module.getID(), i + 2);
      }
 
      if (totalRateL1_Hz > 0.0) {
 
        const double t_ns = 1.0e9 / totalRateL1_Hz;       // [ns]
 
        double t1 = period.getLowerLimit() + gRandom->Exp(t_ns);
 
        if (t1 < period.getUpperLimit()) {
 
 
          // configure correlation matrix
 
          for (int i = 0; i != N; ++i) {
            
            probability2D[i][i] = 0.0;
 
            for (int j = i + 1; j != N; ++j) {
              
              const double ct = JMATH::getDot(module[i], module[j]);
              const double p  = getProbability(ct);
              
              probability2D[i][j] = p;
              probability2D[j][i] = p;
            }
          }
 
 
          // determine probability of a coincidence as a function of the PMT number
          
          for (int i = 0; i != N; ++i) {
 
            probability1D[i] = 0.0;
 
            for (int j = 0; j != N; ++j) {
              probability1D[i] += probability2D[i][j];
            }
          }
 
 
          for ( ; t1 < period.getUpperLimit(); t1 += gRandom->Exp(t_ns)) {
 
 
            // generate two-fold coincidence
 
            const unsigned int pmt1 = getRandomIndex(probability1D,       gRandom->Rndm());
            const unsigned int pmt2 = getRandomIndex(probability2D[pmt1], gRandom->Rndm());
 
            output[pmt1].insert(JPMTSignal(gRandom->Gaus(t1, getSigma()), 1));
            output[pmt2].insert(JPMTSignal(gRandom->Gaus(t1, getSigma()), 1));
 
 
            try {
 
              // generate larger than two-fold coincidences, if any
 
              unsigned int M = getRandomIndex(rateL1_Hz, gRandom->Rndm());
 
              if (M != 0) {
 
                probabilityND = probability2D[pmt1];
 
                for (unsigned int pmtN = pmt2; M != 0; --M) {
 
                  for (int i = 0; i != N; ++i) {
                    probabilityND[i] *= probability2D[pmtN][i];
                  }
 
                  probabilityND[pmtN] = 0.0;
                  
                  pmtN = getRandomIndex(probabilityND, gRandom->Rndm());
                  
                  output[pmtN].insert(JPMTSignal(gRandom->Gaus(t1, getSigma()), 1));
                }
              }
            }
            catch (const JNumericalPrecision&) {}
          }
        }
      }
    }
 
 
    /**
     * Get index based on random value.
     *
     * It is assumed that the values of the input buffer monotonously decrease or increase.
     * This method throws an exception when the summed values in the input buffer is zero.
     *
     * \param  buffer         input  values
     * \param  random         random value <0,1]
     * \return                index
     */
    static inline unsigned int getRandomIndex(const JBuffer1D_t& buffer, const double random)
    {
      double x = 0.0;
 
      for (JBuffer1D_t::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
        x += *i;
      }
 
      if (x > 0.0) {
 
        x *= random;
 
        unsigned int index = 0;
 
        for (JBuffer1D_t::const_iterator i = buffer.begin(); i != buffer.end() && (x -= *i) > 0.0; ++i, ++index) {}
 
        if (index == buffer.size()) {
          --index;
        }
 
        return index;
 
      } else {
 
        THROW(JNumericalPrecision, "getRandomIndex(): zero or negative probability.");
      }
    }
 
 
    /**
     * Get intrinsic time smearing of K40 coincidences.
     *
     * \return                sigma [ns]
     */
    static double getSigma()
    {
      return 0.5;
    }
 
  private:
    /**
     * This correlation matrix is a two-dimensional array in which element [i][j] 
     * corresponds to the probability of a genuine coincidence due to K40 decays, 
     * where i and j refer to the indices of the PMTs in the optical module.
     */
    mutable JBuffer2D_t probability2D;
 
    /**
     * This probability vector is a one-dimensional array in which element [i] 
     * corresponds to the probability of a genuine coincidence due to K40 decays, 
     * where i refers to the index of the PMT in the optical module.
     */
    mutable JBuffer1D_t probability1D;
 
    /**
     * This probability vector is a one-dimensional array in which element [i] 
     * corresponds to the probability of an additional hit,
     * where i refers to the index of the PMT in the optical module.
     */
    mutable JBuffer1D_t probabilityND;
 
    /**
     * Multiples rate as a function of the multiplicity.
     * The index i corresponds to multiplicity M = i + 2.
     */
    mutable JBuffer1D_t rateL1_Hz;
  }; 
}
 
#endif