JFitToT.hh

Go to the documentation of this file.

#ifndef __JCALIBRATE__JFITTOT__
#define __JCALIBRATE__JFITTOT__

#include <limits>

#include "TMath.h"
#include "TString.h"
#include "TF1.h"
#include "TH1.h"
#include "TFitResult.h"

#include "JTools/JRange.hh"

#include "JDetector/JPMTParameters.hh"
#include "JDetector/JPMTAnalogueSignalProcessor.hh"

#include "JROOT/JRootToolkit.hh"


/**
 * \author mkarel, bjung
 */

namespace JCALIBRATE {}
namespace JPP { using namespace JCALIBRATE; }

namespace JCALIBRATE {

  using JTOOLS::JRange;
  using JROOT::JFitParameter_t;
  using JDETECTOR::JPMTParameters;
  using JDETECTOR::JPMTAnalogueSignalProcessor;

  typedef JRange<Double_t>   JRange_t;

  static const double FITTOT_GAIN_MIN        =  0.25;   //!< Minimal gain 
  static const double FITTOT_GAIN_MAX        =  2.00;   //!< Maximal gain

  static const double FITTOT_GAINSPREAD_MIN  =  0.05;   //!< Minimal gain spread 
  static const double FITTOT_GAINSPREAD_MAX  =  1.00;   //!< Maximal gain spread

  static const double FITTOT_TOT_LMARGIN     = 15.00;   //!< Fit-region left  of peak [ns]
  static const double FITTOT_TOT_RMARGIN     = 10.00;   //!< Fit-region right of peak [ns]

  static const std::string  FITTOT_FNAME     = "fittot";

  /**
   * Fit parameters for two-fold coincidence rate due to K40.
   */
  struct JFitToTParameters {
    /**
     * Constructor.
     *
     * \param  parameters   PMT parameters
     */
    JFitToTParameters(const JPMTParameters& parameters) :
      gain              (parameters.gain),
      gainSpread        (parameters.gainSpread)
    {}


    /**
     * Copy constructor.
     *
     * \param  data         data
     */
    JFitToTParameters(const Double_t* data) :
      gain              (0.0),
      gainSpread        (0.0)
    {
      if (data != NULL) {
        setModelParameters(data);
      }
    }


    /** 
     * Get number of model parameters.
     *
     * \return              number of parameters
     */
    static Int_t getNumberOfModelParameters()
    {
      return sizeof(JFitToTParameters) / sizeof(Double_t);
    }


    /** 
     * Get model parameters.
     *
     * \return              pointer to parameters
     */
    const Double_t* getModelParameters() const
    {
      return &this->gain; 
    }


    /** 
     * Get model parameters.
     *
     * \return              pointer to parameters
     */
    Double_t* getModelParameters()
    {
      return &this->gain; 
    }


    /** 
     * Set model parameter.
     *
     * \param  i            parameter index
     * \param  value        parameter value
     */
    void setModelParameter(const int i, const Double_t value)
    {
      getModelParameters()[i] = value;
    }


    /** 
     * Set model parameters.
     *
     * \param  data         pointer to parameters
     */
    void setModelParameters(const Double_t* data)
    {
      for (Int_t i = 0; i != getNumberOfModelParameters(); ++i) {
        setModelParameter(i, data[i]);
      }
    }


    /** 
     * Get model parameter.
     *
     * \param  i            parameter index
     * \return              parameter value
     */
    const Double_t getModelParameter(const int i) const
    {
      return getModelParameters()[i];
    }
   

    /**
     * Get model parameter.
     *
     * \param  p            pointer to data member
     * \return              parameter index and value
     */
    JFitParameter_t getModelParameter(Double_t JFitToTParameters::*p) const
    {
      const Int_t i = &(this->*p) - getModelParameters();

      return JFitParameter_t(i, getModelParameter(i));
    }


    // fit parameters

    Double_t gain;             //!< PMT gain
    Double_t gainSpread;       //!< PMT gain spread
  };


  /**
   * Parametrisation of time-over-threshold distribution.
   *
   * Note that for use in ROOT fit operations, the member method JFitToT::getValue is static.\n
   */
  struct JFitToT :
    public JFitToTParameters,
    public TF1
  {
    
    /**
     * Constructor.
     *
     * \param  parameters   parameters
     * \param  range        abscissa range
     */
    JFitToT(const JPMTParameters& parameters, const JRange_t& range) :

      JFitToTParameters(parameters),

      TF1              (FITTOT_FNAME.c_str(),
                        this,
                        &JFitToT::getValue,
                        range.getLowerLimit(),
                        range.getUpperLimit(),
                        JFitToT::getNumberOfModelParameters()),

      fitRange         (range),
      WfitRange        (1.0),
      
      NPE              (1)
    {
      getInstance().setPMTParameters(parameters);

      init_gain       = getInstance().gain;
      init_gainSpread = getInstance().gainSpread; 
    }


    /**
     * Fit initializer
     *
     * \param  h1                   ROOT 1D-histogram
     * \param  gradientThreshold    treshold gradient for search of maximum
     */
    void fitInit(TH1& h1, const double gradientThreshold = -0.005)
    {
      using namespace std;
      using namespace JPP;
      
      // Find gain and gainspread initial values
      const Int_t N = h1.GetXaxis()->GetNbins();

      double max    = h1.GetBinContent(N);
      double mode   = h1.GetBinCenter(N);

      for (int i = h1.GetBinCenter(N-1);
               i > h1.GetBinCenter(h1.FindBin(getInstance().mean_ns));
             --i) {

        const double x = h1.GetBinCenter (i);
        const double y = h1.GetBinContent(i);

        const double gradient = ( (h1.GetBinContent(i-1) - h1.GetBinContent(i+1)) /
                                  (h1.GetBinCenter (i+1) - h1.GetBinCenter (i-1)) );

        if (y > max) {

          mode = x;
          max  = y;
          
        } else if (gradient < gradientThreshold) {

          break;
        }
      }

      init_gain         = getInstance().getNPE(mode);
      init_gainSpread   = ( init_gain / 3.0 > FITTOT_GAINSPREAD_MIN ?
                            init_gain / 3.0 : FITTOT_GAINSPREAD_MIN ); // Range rule of thumb

      
      // Set parameter names and initial values
      const Int_t Ngain        = this->getModelParameter(&JFitToT::gain);
      const Int_t NgainSpread  = this->getModelParameter(&JFitToT::gainSpread);
      
      this->SetParName(Ngain,       MAKE_CSTRING("gain"));
      this->SetParName(NgainSpread, MAKE_CSTRING("gainSpread"));
      
      getInstance().gain       = init_gain;
      getInstance().gainSpread = init_gainSpread;
      
      setModelParameter(Ngain,       init_gain);
      setModelParameter(NgainSpread, init_gainSpread);
      
      this->SetParameters(this->getModelParameters());


      // Set parameter ranges
      setParLimits(*this, Ngain,       FITTOT_GAIN_MIN,       FITTOT_GAIN_MAX);
      setParLimits(*this, NgainSpread, FITTOT_GAINSPREAD_MIN, FITTOT_GAINSPREAD_MAX);


      // Set fit range and compute corresponding histogram weight
      if (!fitRange.is_valid()) {
        fitRange.setRange(mode - FITTOT_TOT_LMARGIN,
                          mode + FITTOT_TOT_RMARGIN);
      }

      fitRange.setRange(std::max(h1.GetXaxis()->GetXmin(), fitRange.getLowerLimit()),
                        std::min(h1.GetXaxis()->GetXmax(), fitRange.getUpperLimit()));

      WfitRange = h1.Integral(h1.FindBin(fitRange.getLowerLimit()),
                              h1.FindBin(fitRange.getUpperLimit()));
      
      this->SetRange(fitRange.getLowerLimit(), fitRange.getUpperLimit());
    }
    

    /**
     * Fit histogram.
     *
     * Note that the PMT parameters which are part of the model are reset before the fit according the status of each PMT and
     * the obtained fit parameters are copied back to the model parameters after the fit.
     *
     * \param  h1           ROOT 1D-histogram
     * \param  minWeight    minimum histogram weight in fit-range
     * \param  option       fit option
     */
    TFitResultPtr operator()(TH1& h1, const double minWeight, const std::string& option) 
    {
      using namespace std;
      using namespace JPP;

      fitInit(h1);

      if (WfitRange >  minWeight &&
          init_gain >= FITTOT_GAIN_MIN && init_gain <= FITTOT_GAIN_MAX) {

        const TFitResultPtr result = h1.Fit(this, option.c_str());
        
        this->setModelParameters(this->GetParameters());

        const double ToTcheck = getInstance().getTimeOverThreshold(gain);

        return (fitRange(ToTcheck) ? result : TFitResultPtr());
        
      } else {

        return TFitResultPtr();
      }
    }
    

    /**
     * Get rate as a function of the fit parameters.
     *
     * \param  x            pointer to abscissa  values
     * \param  par          pointer to parameter values
     * \return              rate distribution at specified time-over-threshold [Hz/ns]
     */ 
    Double_t getValue(Double_t* x, Double_t* par)
    {
      const Double_t tot_ns = x[0];

      getInstance().gain        = par[0];
      getInstance().gainSpread  = par[1];

      // Compute total weight (= histogram weight within fit range
      //                         plus pdf weight outside fit range)
      const double Wtot = (WfitRange +
                           1.0 - getInstance().getIntegralOfTimeOverThresholdProbability(fitRange.getLowerLimit(),
                                                                                         fitRange.getUpperLimit(),
                                                                                         NPE));
      
      return Wtot * getInstance().getTimeOverThresholdProbability(tot_ns, NPE);
    }


    /**
     * Retrieve initial gain
     *
     *\return initial gain setting
     */
    const double getInitGain() const {
      return init_gain;
    }


    /**
     * Retrieve initial gainspread
     *
     *\return initial gainspread setting
     */
    const double getInitGainSpread() const {
      return init_gainSpread;
    }


  private:

    /**
     * Get unique instance of PMT analogue signal processor.
     *
     * \return              reference to PMT analogue signal processor.
     */
    static JPMTAnalogueSignalProcessor& getInstance()
    {
      static JPMTAnalogueSignalProcessor object;

      return object;
    }

    JRange_t fitRange;          //!< fit-range [ns]
    double   WfitRange;         //!< Cumulative weight within fit-range

    double   init_gain;         //!< Initial input value for gain
    double   init_gainSpread;   //!< Initial input value for gainSpread

    const int NPE;              //!< True number of photo-electrons
  };
}

#endif