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),
                           __module.getPMT(pair.second));
    }

  protected:
    const JModule& __module;
    const double   __epsilon;
  };


  /**
   * Coincidence probability of two PMTs within one module.
   *
   * \param  E1          expected counts of PMT 1 in certain time interval
   * \param  E2          expected counts of PMT 2 in certain time interval
   * \param  ED          expected counts on the rest of the DOM, excluding PMT 1 and PMT 2
   * \param  M_min       DOM multiplicity lower limit
   * \param  M_max       DOM multiplicity upper limit
   * \return             coincidence probability of two PMTs in a certain time interval
   */
  inline double coincidenceP(double E1, double E2, double ED, int M_min, int M_max)
  {
    // first,  always calculate the probability for 2-fold multiplicity
    double P = exp(-E1) * E1 * exp(-E2) * E2 * exp(-ED);
      
    // second, calculate the probability of multiplicity at lower limit
    for (int m = 2; m < M_min; m++) {
      P = P * ( ED/(m-1) + E1/m + E2/m );
    }
  
    // third,  add the probabilities of multiplicities upto upper limit
    double P_sum = P;

    for (int m = M_min+1; m <= M_max; m++) {
      P     = P * ( ED/(m-1) + E1/m + E2/m );
      P_sum = P_sum + P;
    }
  
    return P_sum;
  } 
}

#endif