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 "JLang/JException.hh"

#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;

  static const char         FITTOT_SEPARATOR      = '.';
  static const std::string  FITTOT_SUFFIX         = "1ToT";
  static const std::string  FITTOT_FNAME          = "fittot";

  static const double       FITTOT_GAIN_MIN       =  0.25;   //!< Default minimal gain 
  static const double       FITTOT_GAIN_MAX       =  2.00;   //!< Default maximal gain

  static const double       FITTOT_GAINSPREAD_MIN =  0.05;   //!< Default minimal gain spread 
  static const double       FITTOT_GAINSPREAD_MAX =  1.00;   //!< Default maximal gain spread
  

  /**
   * 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),
      normalization     (1.0)
    {}


    /**
     * Copy constructor.
     *
     * \param  data         data
     */
    JFitToTParameters(const Double_t* data) :
      gain              (0.0),
      gainSpread        (0.0),
      normalization     (1.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
    Double_t normalization;
  };


  /**
   * 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 TF1,
    public JFitToTParameters
  {
    
    /**
     * Constructor.
     *
     * \param  parameters   parameters
     * \param  range        abscissa range
     */
    JFitToT(const JPMTParameters& parameters, const JRange<double>& range) :

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

      JFitToTParameters(parameters),
      
      cpu(parameters)
    {}


    /**
     * Access method for the analogue signal processor
     *
     * \return              reference to analogue signal processor
     */
    const JPMTAnalogueSignalProcessor& getCPU() const
    {
      return cpu;
    }

    
    /**
     * 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  option       fit option
     * \return              fit result
     */
    TFitResultPtr operator()(TH1& h1, const std::string& option) 
    {
      using namespace std;
      using namespace JPP;

      // Set initial gain and gain-spread
      const Int_t Ngain          = this->getModelParameter(&JFitToT::JFitToTParameters::gain);
      const Int_t NgainSpread    = this->getModelParameter(&JFitToT::JFitToTParameters::gainSpread);
      const Int_t Nnormalization = this->getModelParameter(&JFitToT::JFitToTParameters::normalization);      
      
      SetParName(Ngain,          MAKE_CSTRING("gain"));
      SetParName(NgainSpread,    MAKE_CSTRING("gainSpread"));
      SetParName(Nnormalization, MAKE_CSTRING("normalization"));

      normalization = h1.Integral(h1.FindBin(this->GetXmin()),
                                  h1.FindBin(this->GetXmax()));

      SetParameters(getModelParameters());

      // Set parameter and fit ranges
      setParLimits(*this, Ngain,       FITTOT_GAIN_MIN,       FITTOT_GAIN_MAX);
      setParLimits(*this, NgainSpread, FITTOT_GAINSPREAD_MIN, FITTOT_GAINSPREAD_MAX);
      
      FixParameter(Nnormalization, normalization);

      // Fit histogram
      const TFitResultPtr result = h1.Fit(this, option.c_str());
        
      this->setModelParameters(this->GetParameters());
      
      cpu.gain       = getModelParameter(Ngain);
      cpu.gainSpread = getModelParameter(NgainSpread);
        
      return result;
    }
    

    /**
     * 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(const Double_t* x, const Double_t* par)
    {
      const double tot_ns = x[0];

      // Set new parameter values
      cpu.gain            = par[0]; 
      cpu.gainSpread      = par[1];

      // Determine normalization factor
      const int    NPE    = 1;
      const double Whist  = par[2];
      const double Wpdf   = cpu.getIntegralOfTimeOverThresholdProbability(GetXmin(), GetXmax(), NPE);
      
      return (Wpdf > 0.0 ? (Whist / Wpdf) : Whist) * cpu.getTimeOverThresholdProbability(tot_ns, NPE);
    }

  private:

    JPMTAnalogueSignalProcessor cpu;
  };
}

#endif