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 "JDetector/JPMTParameters.hh"
#include "JDetector/JPMTAnalogueSignalProcessor.hh"
#include "JROOT/JRootToolkit.hh"


/**
 * \author mkarel
 */

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

namespace JCALIBRATE {

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

  static const double FITTOT_TOT_MIN_NS      = 10.0;   //!< Minimal time-over-threshold [ns]
  static const double FITTOT_GAIN_MIN        =  0.0;   //!< Minimal gain [unit]
  static const double FITTOT_GAIN_MAX        = 10.0;   //!< Maximal gain [unit]
  static const double FITTOT_GAINSPREAD_MIN  =  0.0;   //!< Minimal gain spread [unit]
  static const double FITTOT_GAINSPREAD_MAX  = 10.0;   //!< Maximal gain spread [unit]

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


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


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


    /** 
     * Set model parameters.
     *
     * \param  data         pointer to parameters
     */
    void setModelParameters(const Double_t* data)
    {
      for (Int_t i = 0; i != getNumberOfModelParameters(); ++i) {
        getModelParameters()[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 W;                //!< normalisation
    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  xmin         minimal x
     * \param  xmax         maximal x
     */
    JFitToT(const JPMTParameters& parameters,
            const double          xmin,
            const double          xmax) :

      JFitToTParameters(parameters),

      TF1("fittot",
          JFitToT::getValue,
          xmin, xmax,
          JFitToT::getNumberOfModelParameters()),

      Xmin(xmin),
      Xmax(xmax)
    {
      getInstance().setPMTParameters(parameters);
    }


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

      Double_t x0 = numeric_limits<Double_t>::max();
      Double_t y0 = 0.0;
      this->W     = 0.0;
      
      for (Int_t i = 1; i <= h1.GetXaxis()->GetNbins(); ++i) {

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

        if (x >= this->Xmin && 
            x <= this->Xmax) {
          this->W += y;
        }

        if (x >= 0.5*TIME_OVER_THRESHOLD_NS && 
            x <= 2.0*TIME_OVER_THRESHOLD_NS) {

          if (y >= y0) {
            x0 = x;
            y0 = y;
          }
        }
      }

      double xmin = Xmin + x0 - TIME_OVER_THRESHOLD_NS;
      double xmax = Xmax + x0 - TIME_OVER_THRESHOLD_NS;
      
      if (xmin < FITTOT_TOT_MIN_NS) {
        xmin = FITTOT_TOT_MIN_NS;
      }

      this->SetRange(xmin, xmax);

      this->gain = getInstance().getNPE(x0);

      setParameter(*this, this->getModelParameter(&JFitToT::gain));
      setParameter(*this, this->getModelParameter(&JFitToT::gainSpread));

      setParLimits(*this, this->getModelParameter(&JFitToT::gain),       FITTOT_GAIN_MIN,       FITTOT_GAIN_MAX);
      setParLimits(*this, this->getModelParameter(&JFitToT::gainSpread), FITTOT_GAINSPREAD_MIN, FITTOT_GAINSPREAD_MAX);

      this->SetParameters(this->getModelParameters());

      getInstance().gain       = this->gain;
      getInstance().gainSpread = this->gainSpread;

      const TFitResultPtr result = h1.Fit(this, option.c_str());

      this->setModelParameters(this->GetParameters());

      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]
     */
    static Double_t getValue(const Double_t* x, const Double_t* par)
    {
      static const int NPE = 1;
      
      const Double_t W         = par[0];
      getInstance().gain       = par[1];
      getInstance().gainSpread = par[2];
      
      const Double_t tot_ns = x[0];
      const Double_t npe    = getInstance().getNPE(tot_ns);
      
      return W * getInstance().getProbability(npe, NPE) * getInstance().getDerivative(npe);
    }


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

      return object;
    }

    double Xmin;
    double Xmax;
  };
}

#endif