JCalibrateK40.hh

Go to the documentation of this file.

#ifndef __JCALIBRATE_JCALIBRATEK40__
#define __JCALIBRATE_JCALIBRATEK40__
 
#include <cmath>
#include <algorithm>
 
#include "JTools/JCombinatorics.hh"
#include "JDetector/JModule.hh"
#include "JMath/JMathToolkit.hh"
 
/**
 * \author mdejong
 */
 
namespace JCALIBRATE {}
namespace JPP { using namespace JCALIBRATE; }
 
namespace JCALIBRATE {
 
  using JTOOLS::JCombinatorics;
  using JDETECTOR::JModule;
 
 
  /**
   * Histogram naming.
   */
  static const char* const weights_hist_t = "weights_hist";    //!< Name of histogram with normalisation constants.
 
  static const char* const W1_overall_t   = "W1_overall";      //!< Named bin for duration of the run
  static const char* const WS_t           = "WS";              //!< Named bin for time residual bin width
  static const char* const WB_t           = "WB";              //!< Named bin for value of TMax_ns in JCalibrateK40
 
  static const char* const _2S            = ".2S";             //!< Name extension for 2D counts
  static const char* const _1B            = ".1B";             //!< Name extension for 1D background
  static const char* const _1L            = ".1L";             //!< Name extension for 1D live time
  static const char* const _2R            = ".2R";             //!< Name extension for 2D rate measured
  static const char* const _2F            = ".2F";             //!< Name extension for 2F rate fitted
  static const char* const _1D            = ".1D";             //!< Name extension for 1D time offset
  static const char* const _1F            = ".1F";             //!< Name extension for 1D time offset fit
 
  typedef JCombinatorics::pair_type  pair_type;
 
 
  /**
   * Auxiliary class to sort pairs of PMT addresses within optical module.
   */
  struct JPairwiseComparator {
    /**
     * Constructor.
     *
     * \param  module       detector module
     * \param  epsilon      precision
     */
    JPairwiseComparator(const JModule& module,
                        const double   epsilon = 1.0e-5) :
      module (module),
      epsilon(epsilon)
    {}
 
 
    /**
     * Comparison of two pairs of PMT addresses.
     *
     * \param  first        first  pair of PMT addresses
     * \param  second       second pair of PMT addresses
     * \return              true if first pair has larger space angle; else false
     */
    inline bool operator()(const pair_type& first, const pair_type& second) const
    {
      using namespace std;
 
      if      (fabs(this->getDot(first) - this->getDot(second)) > epsilon)
        return this->getDot(first) < this->getDot(second);
      else if (max(first.first, first.second) == max(second.first, second.second))
        return min(first.first, first.second) < min(second.first, second.second);
      else
        return max(first.first, first.second) < max(second.first, second.second);
    }
 
 
    /**
     * Get cosine of space angle between PMT axes.
     *
     * \param  pair         pair of PMT addresses
     * \return              cosine
     */
    inline double getDot(const pair_type& pair) const
    {
      return JMATH::getDot(module.getPMT(pair.first) .getDirection(),
                           module.getPMT(pair.second).getDirection());
    }
 
  protected:
    const JModule& module;
    const double   epsilon;
  };
 
 
  /**
   * PMT combinatorics for optical module.
   */
  struct JCombinatorics_t :
    public JCombinatorics
  {
    /**
     * Default constructor.
     */
    JCombinatorics_t()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  module       detector module
     */
    JCombinatorics_t(const JModule& module)
    {
      this->configure(module.size());
 
      this->sort(JPairwiseComparator(module));
    }
  };
 
 
  /**
   * Get coincidence probability of two PMTs within one module due to random background.
   *
   * \param  E1          expected counts of 1st PMT
   * \param  E2          expected counts of 2nd PMT
   * \param  ED          expected counts of module, excluding PMT 1 and PMT 2
   * \param  M_min       multiplicity lower limit
   * \param  M_max       multiplicity upper limit
   * \return             probability
   */
  inline double getP(const double E1, const double E2, const double ED, const int M_min, const int M_max)
  {
    // 1st, always calculate the probability for twofold multiplicity
 
    double P = E1 * E2 * exp(-(E1 + E2 + ED));
      
    // 2nd, calculate the probability of multiplicity at lower limit
 
    int MD = 1;
 
    for ( ; MD + 2 <= M_min; ++MD) {
      P *= ED/MD;
    }
  
    // 3rd, add the probabilities of multiplicities up to upper limit
 
    double P_sum = P;
 
    for ( ; MD + 2 <= M_max; ++MD) {
      P     *= ED/MD;
      P_sum += P;
    }
  
    return P_sum;
  } 
}
 
#endif