JK40Rates.hh

Go to the documentation of this file.

#ifndef __JPHYSICS__JK40RATES___
#define __JPHYSICS__JK40RATES___

#include <vector>

#include "JMath/JMathSupportkit.hh"


/**
 * \author mdejong
 */

namespace JPHYSICS {}
namespace JPP { using namespace JPHYSICS; }

namespace JPHYSICS {

  /**
   * Type definition of singles rate [Hz].
   */
  typedef double                 JRateL0_t;
  
  /**
   * Type definition of count rate as a function of multiplicty [Hz]
   * The multiples rate start counting at two-fold coincidences.
   */
  typedef std::vector<double>    JRateL1_t;


  /**
   * Type definition of multiplicity.
   */
  typedef size_t                 multiplicity_type;

  
  /**
   * Auxiliary class for K40 rates.
   * The singles rate refers to the counting rate of a PMT and
   * the multiples rate to the rate of genuine coincidences due to K40 decays.
   */
  struct JK40Rates {
    /**
     * Default constructor.
     */
    JK40Rates() :
      rateL0(0.0),
      rateL1()
    {}


    /**
     * Constructor.
     *
     * The multiples rates start counting at two-fold coincidences.
     *
     * \param  rateL0_Hz       singles   rate  [Hz]
     * \param  rateL1_Hz       multiples rates [Hz]
     */
    JK40Rates(const JRateL0_t  rateL0_Hz,
              const JRateL1_t& rateL1_Hz = JRateL1_t()) :
      rateL0(rateL0_Hz),
      rateL1(rateL1_Hz)
    {}


    /**
     * Get reference to unique instance of this class object.
     *
     * This method returns an object with default values.\n
     * The singles and multiples rates are taken from 
     * %KM3NeT internal <a href="https://simulation.pages.km3net.de/input_tables/Simulations_Description.pdf">note</a>
     *
     * \return                 reference to this class object
     */
    static JK40Rates& getInstance()
    {
      static JK40Rates rates_Hz;

      if (rates_Hz.rateL1.empty()) {

        rates_Hz.rateL0 = 5.0e3;             // [Hz]

        rates_Hz.rateL1.push_back(500.0);    // [Hz]
        rates_Hz.rateL1.push_back( 50.0);
        rates_Hz.rateL1.push_back(  5.0);
      }

      return rates_Hz;
    }


    /**
     * Get singles rate.
     *
     * \return                 rate [Hz]
     */
    double getSinglesRate() const
    {
      return rateL0;
    }


    /**
     * Get multiples rate.
     *
     * \param  M               multiplicity (M >= JK40Rates::LOWER_L1_MULTIPLICITY)
     * \return                 rate [Hz]
     */
    double getMultiplesRate(const multiplicity_type M) const
    {
      if (M >= LOWER_L1_MULTIPLICITY && M - LOWER_L1_MULTIPLICITY < (multiplicity_type) rateL1.size())
        return rateL1[M - LOWER_L1_MULTIPLICITY];
      else
        return 0.0;
    }


    /**
     * Get lower multiplicty.
     *
     * \return                 lower multiplicity
     */
    multiplicity_type getLowerL1Multiplicity() const
    {
      return LOWER_L1_MULTIPLICITY;
    }


    /**
     * Get upper multiplicty.
     *
     * \return                 upper multiplicity
     */
    multiplicity_type getUpperL1Multiplicity() const
    {
      return rateL1.size() + 1;
    }

    
    /**
     * Correct rates for global efficiency,
     *
     * \param  QE              global efficiency
     */
    void correct(const double QE)
    {
      if (QE > 0.0) {
      
        rateL0 /= QE;

        JRateL1_t buffer = rateL1;
        
        for (multiplicity_type M = getUpperL1Multiplicity(); M >= getLowerL1Multiplicity(); --M) {

          // determine contribution from higher multiplicities
          
          double R = 0.0;
          
          for (multiplicity_type i = M + 1; i <= getUpperL1Multiplicity(); ++i) {
            R += buffer[i - LOWER_L1_MULTIPLICITY] * JMATH::binomial(i, M) * pow(QE, M) * pow(1.0 - QE, i - M);
          }

          if (getMultiplesRate(M) > R)
            buffer[M - LOWER_L1_MULTIPLICITY] = (getMultiplesRate(M) - R) / pow(QE, M);
          else
            buffer[M - LOWER_L1_MULTIPLICITY] =  0.0;
        }

        rateL1 = buffer;
        
      } else {

        rateL0 = 0.0;

        for (JRateL1_t::iterator i = rateL1.begin(); i != rateL1.end(); ++i) {
          *i = 0.0;
        }
      }
    }
    
    
    /**
     * Read K40 rates from input.
     *
     * \param  in              input stream
     * \param  object          K40 rates
     * \return                 input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JK40Rates& object)
    {
      const double rateL0 = object.rateL0;

      if (in >> object.rateL0) {

        object.rateL1.clear();

        for (double x; in >> x; ) {
          object.rateL1.push_back(x);
        }

      } else {

        object.rateL0 = rateL0;
      }

      return in;
    }


    /**
     * Write K40 rates to output.
     *
     * \param  out             output stream
     * \param  object          K40 rates
     * \return                 output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JK40Rates& object)
    {
      out << object.rateL0;

      for (JRateL1_t::const_iterator i = object.rateL1.begin(); i != object.rateL1.end(); ++i) {
        out << ' ' << *i;
      }

      return out;
    }

    
    /**
     * Lower L1 multiplicity
     */
    static const multiplicity_type  LOWER_L1_MULTIPLICITY  =  2;

    
  protected:
    JRateL0_t rateL0;          //!< singles   rate  [Hz]
    JRateL1_t rateL1;          //!< multiples rates [Hz]
  };
}

#endif