CDCDedxADCNonLinearity.h

/**************************************************************************
 * BASF2 (Belle Analysis Framework 2)                                     *
 * Copyright(C) 2020 - Belle II Collaboration                             *
 
 * Author: The Belle II Collaboration                                     *
 
 * This software is provided "as is" without any warranty.                *
 **************************************************************************/
 
#pragma once
#include <TObject.h>
#include <vector>
#include <cmath>
#include <framework/logging/Logger.h>
#include <iostream>
 
namespace Belle2 {
  class CDCDedxADCNonLinearity: public TObject {
  public:
 
    CDCDedxADCNonLinearity(): m_nonlADC() {};
 
    explicit CDCDedxADCNonLinearity(const std::vector<std::vector<double>>& nonlinearity): m_nonlADC(nonlinearity) {};
 
    ~CDCDedxADCNonLinearity() {};
 
    unsigned int getSize(int layer, int axis) const
    {
 
      int mylayer = 0;
      if (layer < 0 || layer > 55) {
        B2ERROR("CDCDedxADCNonLinearity:CDC layer out of range, choose between 0-55");
        return 0;
      } else {
        if (layer >= 8) mylayer = 2;
      }
 
      if (axis < 0 || axis > 1) {
        B2ERROR("CDCDedxADCNonLinearity:Problem with X/Y axis: choose 0/1 for X/Y");
        return 0;
      } else {
        mylayer = mylayer + axis;
      }
 
      return m_nonlADC[mylayer].size();
    };
 
 
    double getCorrectedADC(const double& ADC, const int layer)const
    {
 
      int mylayer = 0;
      if (layer < 0 || layer > 55) {
        B2WARNING("CDCDedxADCNonLinearity:returning uncorrected ADC as input layer is out of range: must be 0-55");
        return ADC;
      } else {
        if (layer >= 8) mylayer = 2;
      }
 
      std::vector<double> tempX = m_nonlADC[mylayer];//inner or outer X array
      std::vector<double> tempY = m_nonlADC[mylayer + 1];//inner or outer Y array
 
      if (tempX.size() != tempY.size()) {
        B2WARNING("CDCDedxADCNonLinearity:returning uncorrected ADC as parameters range don't match: X=Y in bins");
        return ADC;
      }
 
      //Find bin for ADC correction
      unsigned int ibin = TMath::BinarySearch(tempY.size(), tempY.data(), double(ADC));
 
      if (ibin >= tempY.size() - 1)ibin = tempY.size() - 2; //overflow to last bin
      if (ibin >= tempY.size() - 1) {
        B2WARNING("CDCDedxADCNonLinearity:returning uncorrected ADC as bins are not in range");
        return ADC;
      }
 
      double slope = (tempY[ibin + 1] - tempY[ibin]) / (tempX[ibin + 1] - tempX[ibin]);
      return  std::round(tempX[ibin] + (ADC - tempY[ibin]) / slope);
    }
 
    double getNonLinearityPar(int layer, int axis, unsigned int par) const
    {
 
      int mylayer = 0;
      if (layer < 0 || layer > 55) {
        B2ERROR("CDCDedxADCNonLinearity:CDC layer out of range, choose between 0-55");
        return -99.0;
      } else {
        if (layer >= 8) mylayer = 2;
      }
 
      if (axis < 0 || axis > 1) {
        B2ERROR("CDCDedxADCNonLinearity:Problem with X/Y axis: choose 0/1 for X/Y");
        return -99.0;
      } else {
        mylayer = mylayer + axis;
      }
 
      if (par >= m_nonlADC[mylayer].size()) {
        B2ERROR("CDCDedxADCNonLinearity:Problem with par index: choose 0 and " << m_nonlADC[mylayer].size()); //
        return -99.0;
      }
 
      return m_nonlADC[mylayer][par];
    };
 
    void printNonLinearityPars(int layer, int axis) const
    {
 
      int mylayer = 0;
      if (layer < 0 || layer > 55) {
        B2ERROR("CDCDedxADCNonLinearity:CDC layer out of range, choose between 0-55");
        return;
      } else {
        if (layer >= 8) mylayer = 2;
      }
 
      if (axis < 0 || axis > 1) {
        B2ERROR("CDCDedxADCNonLinearity:Problem with X/Y axis: choose 0/1 for X/Y");
        return;
      } else {
        mylayer = mylayer + axis;
      }
 
      B2INFO("Printing parameters for layer: " << layer << ", axis: " << axis << ", nPars:" << m_nonlADC[mylayer].size());
      for (unsigned int iPar = 0; iPar < m_nonlADC[mylayer].size(); iPar++)
        std::cout << "Par # " << iPar << ": " << m_nonlADC[mylayer][iPar] << std::endl;
 
    };
 
    void setNonLinearityPar(unsigned int layer, unsigned int axis, unsigned int par, double value)
    {
      int mylayer = 0;
      if (layer > 55) {
        B2ERROR("CDCDedxADCNonLinearity:CDC layer out of range, choose between 0-55");
        return;
      } else {
        if (layer >= 8) mylayer = 2;
      }
 
      if (axis > 1) {
        B2ERROR("CDCDedxADCNonLinearity:Problem with X/Y axis: choose 0/1 for X/Y");
        return;
      } else {
        mylayer = mylayer + axis;
      }
 
      if (par >= m_nonlADC[mylayer].size()) {
        B2ERROR("CDCDedxADCNonLinearity:Problem with par index: choose 0 and " << m_nonlADC[mylayer].size()); //
      }
 
      //here set parameter to requested value
      m_nonlADC[mylayer][par] = value;
 
    };
 
 
  private:
    std::vector<std::vector<double>> m_nonlADC; 
    ClassDef(CDCDedxADCNonLinearity, 1); 
  };
} // end namespace Belle2