SVDCoGCalibrationFunction.h

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
 
#pragma once
 
#include <framework/utilities/MathHelpers.h>
 
#include <TObject.h>
 
#include <cmath>
#include <vector>
 
namespace Belle2 {

  class SVDCoGCalibrationFunction : public TObject {
 
  public:

    typedef double (SVDCoGCalibrationFunction::*cogFunction)(double, int) const;

    typedef double (SVDCoGCalibrationFunction::*cogFunctionErr)(double, double, int) const;

    double calibratedValue(double raw_time, int trigger_bin) const
    {
      cogFunction f = m_implementations[m_current];
      return (this->*f)(raw_time, trigger_bin) ;
    }

    double calibratedValueError(double raw_time, double raw_timeErr, int trigger_bin) const
    {
      cogFunctionErr fErr = m_implementationsErr[m_current];
      return (this->*fErr)(raw_time, raw_timeErr, trigger_bin) ;
    }

    SVDCoGCalibrationFunction()
    {
      // The m_implementations vector is static.
      // We have to initialize it just once.
      if (m_implementations.size() == 0) {
        m_implementations.push_back(&SVDCoGCalibrationFunction::pol1TBdep);
        m_implementations.push_back(&SVDCoGCalibrationFunction::pol3TBindep);
        m_implementations.push_back(&SVDCoGCalibrationFunction::pol5TBindep);
        m_implementations.push_back(&SVDCoGCalibrationFunction::elsTBindep);
        //m_implementations.push_back(
        //  &SVDCoGCalibrationFunction::betterVersion);
      }
 
      // The m_implementationsErr vector is static.
      // We have to initialize it just once.
      if (m_implementationsErr.size() == 0) {
        m_implementationsErr.push_back(&SVDCoGCalibrationFunction::pol1TBdepErr);
        m_implementationsErr.push_back(&SVDCoGCalibrationFunction::pol3TBindepErr);
        m_implementationsErr.push_back(&SVDCoGCalibrationFunction::pol5TBindepErr);
        m_implementationsErr.push_back(&SVDCoGCalibrationFunction::elsTBindepErr);
        //m_implementationsErr.push_back(
        //  &SVDCoGCalibrationFunction::betterVersion);
      }
 
      m_current = m_implementations.size() - 1;
 
    };

    void set_current(int current)
    {
      m_current = current;
    }
 
    //SETTERS FOR function ID = 0 (pol1TBdep)
    void set_bias(double tb0, double tb1, double tb2, double tb3)
    {
      m_bias[0] = tb0;
      m_bias[1] = tb1;
      m_bias[2] = tb2;
      m_bias[3] = tb3;
    }

    void set_scale(double tb0, double tb1, double tb2, double tb3)
    {
      m_scale[0] = tb0;
      m_scale[1] = tb1;
      m_scale[2] = tb2;
      m_scale[3] = tb3;
    }
 
 
    //SETTERS FOR function ID = 1 (pol3TBindep)
    void set_pol3parameters(double a, double b, double c, double d)
    {
      m_par[ 0 ] = a;
      m_par[ 1 ] = b;
      m_par[ 2 ] = c;
      m_par[ 3 ] = d;
    }
 
    //SETTERS FOR function ID = 2 (pol5TBindep)
    void set_pol5parameters(double a, double b, double c, double d, double e, double f)
    {
      m_par[ 0 ] = a;
      m_par[ 1 ] = b;
      m_par[ 2 ] = c;
      m_par[ 3 ] = d;
      m_par[ 4 ] = e;
      m_par[ 5 ] = f;
    }
 
    //SETTERS FOR function ID = 3 (elsTBindep)
    void set_elsparameters(double a, double b, double c, double d)
    {
      m_par[ 0 ] = a;
      m_par[ 1 ] = b;
      m_par[ 2 ] = c;
      m_par[ 3 ] = d;
    }

    SVDCoGCalibrationFunction(const Belle2::SVDCoGCalibrationFunction& a);

    SVDCoGCalibrationFunction& operator=(const Belle2::SVDCoGCalibrationFunction& a);
 
  private:

    static const int nTriggerBins = 4;


    double m_bias[ nTriggerBins ] = {0}; 
    double m_scale[ nTriggerBins ] = {0}; 
    double pol1TBdep(double raw_time, int tb) const
    {
      return raw_time * m_scale[ tb % nTriggerBins] +
             m_bias[ tb % nTriggerBins ];
    };

    double pol1TBdepErr(double, double raw_timeErr, int tb) const
    {
      return raw_timeErr * m_scale[ tb % nTriggerBins];
    };
 
    //data member useful for polynomials
    static const int m_nPar = 6; 
    double m_par[ m_nPar ] = {0}; 

    double pol3TBindep(double raw_time, int) const
    {
      return m_par[0] + m_par[1] * raw_time + m_par[2] * square(raw_time) + m_par[3] * cube(raw_time);
    };

    double pol3TBindepErr(double raw_time, double raw_timeErr, int) const
    {
      return raw_timeErr * (m_par[1] +  2 * m_par[2] * raw_time + 3 * m_par[3] * square(raw_time));
    };
 

    double pol5TBindep(double raw_time, int) const
    {
      return m_par[0] + m_par[1] * raw_time + m_par[2] * square(raw_time) + m_par[3] * cube(raw_time) + m_par[4] * pow4(
               raw_time) + m_par[5] * pow5(raw_time);
    };

    double pol5TBindepErr(double raw_time, double raw_timeErr, int) const
    {
      return raw_timeErr * (m_par[1] +  2 * m_par[2] * raw_time + 3 * m_par[3] * square(raw_time) + 4 * m_par[4] * cube(
                              raw_time) + 5 * m_par[5] * pow4(raw_time));
    };

    double elsTBindep(double raw_time, int) const
    {
      if (raw_time > m_par[3] - sqrt(-m_par[2]) / 4)
        return  std::numeric_limits<float>::quiet_NaN();
 
      return m_par[0] + m_par[1] * raw_time + m_par[2] / (raw_time - m_par[3]);
    };

    double elsTBindepErr(double raw_time, double raw_timeErr, int) const
    {
      if (raw_time > m_par[3] - sqrt(-m_par[2]) / 4)
        return  std::numeric_limits<float>::quiet_NaN();
 
      return raw_timeErr * (m_par[1] - m_par[2] / square(raw_time - m_par[3]));
    };

    int m_current;

    static std::vector < cogFunction > m_implementations; 

    static std::vector < cogFunctionErr > m_implementationsErr; 
 
 
    ClassDef(SVDCoGCalibrationFunction, 6)
  };

}