SVDCoGCalibrationFunction.h

/**************************************************************************
\ * BASF2 (Belle Analysis Framework 2)                                     *
 * Copyright(C) 2017 - Belle II Collaboration                             *
 *                                                                        *
 * Author: The Belle II Collaboration                                     *
 * Contributors: Eugenio Paoloni, Giulia Casarosa                         *
 *                                                                        *
 * This software is provided "as is" without any warranty.                *
 **************************************************************************/
 
#pragma once
 
#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);
 
 
  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 usefule for polinomials
    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] * pow(raw_time, 2) + m_par[3] * pow(raw_time, 3);
    };
    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] * pow(raw_time, 2));
    };
 
 
    double pol5TBindep(double raw_time, int) const
    {
      return m_par[0] + m_par[1] * raw_time + m_par[2] * pow(raw_time, 2) + m_par[3] * pow(raw_time, 3) + m_par[4] * pow(raw_time,
             4) + m_par[5] * pow(raw_time, 5);
    };
    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] * pow(raw_time, 2) + 4 * m_par[4] * pow(raw_time,
                            3) + 5 * m_par[5] * pow(raw_time, 4));
    };
 
    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] / pow(raw_time - m_par[3], 2));
    };
 
    int m_current;
 
    static std::vector < cogFunction > m_implementations; 
 
    static std::vector < cogFunctionErr > m_implementationsErr; 
 
 
    ClassDef(SVDCoGCalibrationFunction, 6)
  };
 
}