JK40DefaultSimulatorInterface.hh

Go to the documentation of this file.

#ifndef __JDETECTOR__JK40DEFAULTSIMULATORINTERFACE__
#define __JDETECTOR__JK40DEFAULTSIMULATORINTERFACE__
 
#include <utility>
#include <vector>
#include <set>
 
#include "TRandom3.h"
 
#include "JPhysics/JConstants.hh"
 
#include "JDetector/JK40Simulator.hh"
#include "JDetector/JModuleIdentifier.hh"
#include "JDetector/JPMTIdentifier.hh"
#include "JDetector/JTimeRange.hh"
 
#include "JMath/JMathToolkit.hh"
#include "JLang/JComparator.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
  {
    /**
     * Auxiliary data structure for value with associated probability.
     */
    template<class T>
    struct probability_type {
      T      v;        //!< value
      double P;        //!< probability
    };
 
 
  public:
    /**
     * Auxiliary data structure for selecting values according associated probabilities.
     */
    template<class T>
    struct cdf_type :
      public std::vector< probability_type<T> >
    {
      /**
       * Get random value according associated probabilities.
       *
       * \param  rv             random vavlue [0,1>
       * \return                value
       */
      const T& operator()(const double rv) const
      {
        const double P = this->rbegin()->P * rv;
 
        size_t first = 0;
        size_t count = this->size();
 
        for (size_t i, step; count != 0; ) {
 
          step = count / 2;
          i    = first + step;
 
          if ((*this)[i].P < P) {
            first  = ++i;
            count -= step + 1;
          } else {
            count  = step;
          }
        }
 
        return (*this)[first].v;
      }
 
 
      /**
       * Put value with given probability.
       *
       * \param  v                  value
       * \param  P                  probability
       */
      void put(const T& v, const double P)
      {
        if (this->empty())
          this->push_back({v, P});
        else
          this->push_back({v, this->rbegin()->P + P});
      }
    };
    
 
  protected:
    /**
     * 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;
 
 
    /**
     * Get mixed L0 (lower module) and L1 (upper module) rate [Hz]
     *
     * \return                 rate [Hz]
     */
    virtual double getL01Rate() const = 0;
 
 
    /**
     * Get mixed L1 (lower module) and L0 (upper module) rate [Hz]
     *
     * \return                 rate [Hz]
     */
    virtual double getL10Rate() 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 override
    {
      using namespace std;
      using namespace JPP;
      
      const size_t N = module.size();
 
 
      // generate singles
 
      for (size_t pmt = 0; pmt != N; ++pmt) {
 
        const double rateL0_Hz = getSinglesRate(JPMTIdentifier(module.getID(), pmt));
 
        if (rateL0_Hz > 0.0) {
 
          const bool   asap = (output[pmt].size() != 0);  // sort as soon as possible
 
          const double t_ns = 1.0e9 / rateL0_Hz;          // interval period [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));    // see JPMTSignalProcessorInterface::operator()
          }
 
          if (asap) {
            output[pmt].sort();
          }
        }
      }
 
 
      // generate coincidences
 
      double totalRateL1_Hz = 0.0;
 
      // multiples rate as a function of the multiplicity - index i corresponds to multiplicity M = i + 2
 
      vector<double> rateL1_Hz(N, 0.0);      
 
      for (size_t i = 0; i != N; ++i) {
        totalRateL1_Hz += rateL1_Hz[i] = getMultiplesRate(module.getID(), i + 2);
      }
 
      if (totalRateL1_Hz > 0.0) {
 
        typedef pair<size_t,size_t>  pair_type;
      
        cdf_type<pair_type> probablityL1;
 
        const double t_ns = 1.0e9 / totalRateL1_Hz;       // interval period [ns]
 
        double t1 = period.getLowerLimit() + gRandom->Exp(t_ns);
 
        if (t1 < period.getUpperLimit()) {
 
          // configure pair-wise propabilities
 
          for (size_t pmt1 = 0; pmt1 != N; ++pmt1) {
            for (size_t pmt2 = 0; pmt2 != pmt1; ++pmt2) {
 
              const double ct = getDot(module[pmt1].getDirection(), module[pmt2].getDirection());
              const double p  = getProbability(ct);
 
              probablityL1.put({pmt1, pmt2}, p);
            }
          }
          
          for ( ; t1 < period.getUpperLimit(); t1 += gRandom->Exp(t_ns)) {
 
            try {
 
              // generate two-fold coincidence
 
              const pair_type& pair = probablityL1(gRandom->Rndm());
 
              output[pair.first ].insert(JPMTSignal(gRandom->Gaus(t1, getSigma()), 1));
              output[pair.second].insert(JPMTSignal(gRandom->Gaus(t1, getSigma()), 1));
              
              // generate larger than two-fold coincidences, if any
 
              size_t M = 0;
 
              for (double R = totalRateL1_Hz * gRandom->Rndm(); M != N && (R -= rateL1_Hz[M]) > 0.0; ++M) {}
 
              if (M != 0) {
 
                set<size_t> buffer = { pair.first , pair.second };   // hit PMTs
 
                for ( ; M != 0; --M) {
 
                  vector<double> probability1D(N, 1.0);
 
                  double P = 0.0;
                  
                  for (size_t i = 0; i != N; ++i) {
                    
                    if (buffer.count(i) == 0) {
 
                      for (set<size_t>::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
 
                        const double ct = getDot(module[i].getDirection(), module[*pmt].getDirection());
 
                        probability1D[i] *= getProbability(ct);
                      }
 
                      P += probability1D[i];
 
                    } else {
 
                      probability1D[i] = 0.0;
                    }
                  }
 
                  if (P > 0.0) {
 
                    size_t pmt = 0;
 
                    for (P *= gRandom->Rndm(); pmt != N && (P -= probability1D[pmt]) > 0.0; ++pmt) {}
 
                    if (pmt != N) {
 
                      output[pmt].insert(JPMTSignal(gRandom->Gaus(t1, getSigma()), 1));
 
                      buffer.insert(pmt);
                    }
 
                  } else {
 
                    break;
                  }
                }
              }
            }
            catch (const JNumericalPrecision&) {}
          }
        }
      }
    }
 
 
    /**
     * Generate mixed-L1/L0 hits.
     *
     * \param  input          module pair
     * \param  period         time window [ns]
     * \param  output         background data
     */
    virtual void generateHits(const module_pair& input,
                              const JTimeRange&  period,
                              const module_data& output) const override
    {
      using namespace JPP;
 
      const double rate_Hz = getL01Rate() + getL10Rate();
 
      if (rate_Hz > 0.0) {
 
        if (!input.first.empty() && !input.second.empty()) {
 
          if (neighbours(input.first, input.second)) {
 
            // locate decay in between modules
 
            const double ts = 0.8 * fabs(input.first.getZ() - input.second.getZ()) * getInverseSpeedOfLight() * getIndexOfRefraction();
 
            const JModule&   input_lower  = (input.first.getFloor() == input.second.getFloor() - 1 ? input.first  : input.second);
            const JModule&   input_upper  = (input.first.getFloor() == input.second.getFloor() + 1 ? input.first  : input.second);
 
            JModuleData&     output_lower = (input.first.getFloor() == input.second.getFloor() - 1 ? output.first : output.second);
            JModuleData&     output_upper = (input.first.getFloor() == input.second.getFloor() + 1 ? output.first : output.second);
 
            cdf_type<size_t> probability_lower;
            cdf_type<size_t> probability_upper;
 
            // use angular dependence of coincidence rate with respect to z-axis
 
            for (size_t i = 0; i != input_lower.size(); ++i) { probability_lower.put(i, this->getProbability(+1.0 * input_lower[i].getDZ())); }
            for (size_t i = 0; i != input_upper.size(); ++i) { probability_upper.put(i, this->getProbability(-1.0 * input_upper[i].getDZ())); }
              
            const double t_ns = 1.0e9 / rate_Hz;          // interval period [ns]
 
            for (double t1 = period.getLowerLimit() + gRandom->Exp(t_ns); t1 < period.getUpperLimit(); t1 += gRandom->Exp(t_ns)) {
 
              const bool L01 = gRandom->Uniform(0.0, rate_Hz) < getL01Rate();
 
              JModuleData&      output_L1      = (L01 ? output_upper      : output_lower);
              JModuleData&      output_L0      = (L01 ? output_lower      : output_upper);
 
              cdf_type<size_t>& probability_L1 = (L01 ? probability_upper : probability_lower);
              cdf_type<size_t>& probability_L0 = (L01 ? probability_lower : probability_upper);
 
              output_L1[probability_L1(gRandom->Rndm())].insert(JPMTSignal(gRandom->Gaus(t1, getSigma()), 1));
              output_L1[probability_L1(gRandom->Rndm())].insert(JPMTSignal(gRandom->Gaus(t1, getSigma()), 1));
 
              output_L0[probability_L0(gRandom->Rndm())].insert(JPMTSignal(gRandom->Gaus(t1 + ts, getSigma()), 1));
            }
          }
        }
      }
    }
 
 
    /**
     * Get intrinsic time smearing of K40 coincidences.
     *
     * \return                sigma [ns]
     */
    static double getSigma()
    {
      return get_sigma();
    }
 
 
    /**
     * Set intrinsic time smearing of K40 coincidences.
     *
     * \param  sigma          sigma [ns]
     */
    static void setSigma(const double sigma)
    {
      get_sigma() = sigma;
    }
 
  private:
    /**
     * Get intrinsic time smearing of K40 coincidences.
     *
     * \return                sigma [ns]
     */
    static double& get_sigma()
    {
      static double sigma = 0.5;
 
      return sigma;
    }
  }; 
}
 
#endif