JK40.hh

Go to the documentation of this file.

#ifndef __JFIT__JK40__
#define __JFIT__JK40__

#include <cmath>
#include <vector>

#include "JPhysics/JK40Rates.hh"

#include "JFit/JTimeRange.hh"
#include "JFit/JFitToolkit.hh"


/**
 * \author mdejong
 */

namespace JFIT {}
namespace JPP { using namespace JFIT; }

namespace JFIT {

  using JPHYSICS::JK40Rates;
  
  
  /**
   * Auxiliary class for converting various rates to expectation values
   * of the number of hits within a given time interval.
   */
  struct JK40 {
    /**
     * Type definition of expectation values of number of multiple hits.
     */
    typedef std::vector<double>  JL1_t;

    
    /**
     * Default constructor.
     */
    JK40() :
      __y0(0.0),
      __y1()
    {}


    /**
     * Constructor for single PMT.
     *
     * \param  T_ns               time range [ns]
     * \param  R_Hz               K40 singles rate [Hz]
     */
    JK40(const JTimeRange& T_ns,
         const double      R_Hz) :
      __y0(0.0),
      __y1()
    {
      const double T_s = (T_ns.getUpperLimit() - T_ns.getLowerLimit()) * 1.0e-9;   // period [s]
      
      __y0 = R_Hz * T_s;
    }
    

    /**
     * Constructor for modules.
     *
     * \param  T_ns               time range [ns]
     * \param  numberOfPMTs       number of PMTs in module
     * \param  rates_Hz           K40 singles and multiples rates [Hz]
     */
    JK40(const JTimeRange& T_ns,
         const size_t      numberOfPMTs,
         const JK40Rates&  rates_Hz) :
      __y0(0.0),
      __y1()
    {
      const double T_s = (T_ns.getUpperLimit() - T_ns.getLowerLimit())* 1.0e-9;   // period [s]
      
      __y0 = numberOfPMTs * rates_Hz.getSinglesRate() * T_s;

      // genuine coincidences
      
      for (size_t M = rates_Hz.getLowerL1Multiplicity(); M <= rates_Hz.getUpperL1Multiplicity(); ++M) {
        __y1.push_back(rates_Hz.getMultiplesRate(M) * T_s);
      }

      // random coincidences
      
      if (__y1.empty()) {
        __y1.push_back(0.0);
      }

      __y1[0] +=  __y0 * (1.0 - exp(-__y0));
    }


    /**
     * Get expectation value for number of single hits.
     *
     * \return                    expectation value
     */
    double getY0() const
    {
      return __y0;
    }


    /**
     * Get expectation value for number of multiple hits.
     *
     * \param  M                  multiplicity (M >= 2)
     * \return                    expectation value
     */
    double getY1(const size_t M) const
    {
      if (M >= 2 && M - 2 < __y1.size()) 
        return __y1[M-2];
      else
        return 0.0;
    }
    

    /**
     * Get probability to observe a hit with given multiplicity.
     * Note that if M = 0, the probability corresponds to observing no hits.
     *
     * \param  M                  multiplicity (M >= 0)
     * \return                    probability
     */
    double getP(size_t M) const
    {
      double expval = 0.0;
      
      if (M < 2)
        expval = getY0();
      else
        expval = getY1(M);

      return JFIT::getP(expval, M != 0);
    }

  protected:
    double __y0;           //!< expectation value  of number of single   hits
    JL1_t  __y1;           //!< expectation values of number of multiple hits
  };
}

#endif