CDCDedxInjectTimeAlgorithm.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 <map>
#include <vector>
#include <iostream>
#include <fstream>
 
#include <TH1D.h>
#include <TH1I.h>
#include <TF1.h>
#include <TCanvas.h>
#include <TLine.h>
#include <TLegend.h>
 
#include <calibration/CalibrationAlgorithm.h>
#include <framework/database/DBObjPtr.h>
#include <cdc/dbobjects/CDCDedxInjectionTime.h>
#include <cdc/utilities/CDCDedxMeanPred.h>
#include <cdc/utilities/CDCDedxSigmaPred.h>
#include <framework/gearbox/Const.h>
 
// namespace constants
namespace numdedx {
  static const int nrings = 2;
  static const int nvectors = 6;
}
 
enum fstatus {fitOK, fitFailed};
 
namespace Belle2 {
  class CDCDedxInjectTimeAlgorithm : public CalibrationAlgorithm {
 
  public:

    CDCDedxInjectTimeAlgorithm();

    virtual ~CDCDedxInjectTimeAlgorithm() {}

    void setMergePayload(bool value = true) {m_isMerge = value;}

    void setMonitoringPlots(bool value = false) {m_ismakePlots = value;}

    void setFitWidth(double value) {m_sigmaR = value;}

    void setMinTracks(int value) {m_thersE = value;}

    void setDedxPars(int value, double min, double max)
    {
      m_dedxBins = value;
      m_dedxMin = min;
      m_dedxMax = max;
    }

    void setChiPars(int value, double min, double max)
    {
      m_chiBins = value;
      m_chiMin = min;
      m_chiMax = max;
    }

    void fitGaussianWRange(TH1D*& temphist, fstatus& status);

    void getExpRunInfo();

    void defineTimeBins();

    void defineHisto(std::array<std::vector<TH1D*>, numdedx::nrings>& htemp,  std::string var);

    void defineTimeHisto(std::array<std::array<TH1D*, 3>, numdedx::nrings>& htemp);

    void checkStatistics(std::array<std::vector<TH1D*>, numdedx::nrings>& hvar);

    void correctBinBias(std::map<int, std::vector<double>>& varcorr, std::map<int, std::vector<double>>& var,
                        std::map<int, std::vector<double>>& time, TH1D*& htimes);

    void createPayload(std::array<double, numdedx::nrings>& scale, std::map<int, std::vector<double>>& vmeans,
                       std::map<int, std::vector<double>>& varscal, std::string svar);

    void getMeanReso(std::array<std::vector<TH1D*>, numdedx::nrings>& hvar,
                     std::map<int, std::vector<double>>& vmeans, std::map<int, std::vector<double>>& vresos);

    void plotEventStats();

    void plotBinLevelDist(std::array<std::vector<TH1D*>, numdedx::nrings>& hvar, std::string var);

    void plotRelConstants(std::map<int, std::vector<double>>& vmeans, std::map<int, std::vector<double>>& vresos,
                          std::map<int, std::vector<double>>& corr, std::string svar);

    void plotTimeStat(std::array<std::vector<TH1D*>, numdedx::nrings>& htime);

    void plotFinalConstants(std::map<int, std::vector<double>>& vmeans, std::map<int, std::vector<double>>& vresos,
                            std::array<double, numdedx::nrings>& scale,  std::array<double, numdedx::nrings>& scale_reso);

    void plotInjectionTime(std::array<std::array<TH1D*, 3>, numdedx::nrings>& hvar);

    void setHistStyle(TH1D*& htemp, const int ic, const int is, const double min, const double max)
    {
      htemp->SetStats(0);
      htemp->LabelsDeflate();
      htemp->SetMarkerColor(ic);
      htemp->SetMarkerStyle(is);
      htemp->GetXaxis()->SetLabelOffset(-0.055);
      htemp->GetYaxis()->SetTitleOffset(0.75);
      htemp->SetMinimum(min);
      htemp->SetMaximum(max);
    }

    std::string getTimeBinLabel(const double& tedges, const int& it)
    {
      std::string label = "";
      if (tedges < 2e4)label = Form("%0.01f-%0.01fK", m_tedges[it] / 1e3, m_tedges[it + 1] / 1e3);
      else if (tedges < 1e5)label = Form("%0.0f-%0.0fK", m_tedges[it] / 1e3, m_tedges[it + 1] / 1e3);
      else label = Form("%0.01f-%0.01fM", m_tedges[it] / 1e6, m_tedges[it + 1] / 1e6);
      return label;
    }

    void deleteHisto(std::array<std::vector<TH1D*>, numdedx::nrings>& htemp)
    {
      for (unsigned int ir = 0; ir < c_rings; ir++) {
        for (unsigned int it = 0; it < m_tbins; it++) {
          delete htemp[ir][it];
        }
      }
    }

    void deleteTimeHisto(std::array<std::array<TH1D*, 3>, numdedx::nrings>& htemp)
    {
      const int tzoom = 3;
      for (unsigned int ir = 0; ir < c_rings; ir++) {
        for (int wt = 0; wt < tzoom; wt++) {
          delete htemp[ir][wt];
        }
      }
    }

    double getCorrection(unsigned int ring, unsigned int time, std::map<int, std::vector<double>>&  vmeans);
 
 
  protected:

    virtual EResult calibrate() override;
 
  private:
 
    static const int c_rings = numdedx::nrings;  
    std::array<std::string, numdedx::nrings> m_sring{"ler", "her"};
 
    std::vector<double> m_vtedges; 
    std::vector<double> m_vtlocaledges; 
 
    double* m_tedges; 
    unsigned int m_tbins;  
    double m_sigmaR; 
 
    int m_dedxBins; 
    double m_dedxMin; 
    double m_dedxMax; 
 
    int m_chiBins; 
    double m_chiMin; 
    double m_chiMax; 
 
    int m_countR; 
    int m_thersE; 
 
    bool m_isminStat; 
    bool m_ismakePlots; 
    bool m_isMerge; 
 
    std::string m_prefix; 
    std::string m_suffix; 
 
    std::vector<std::vector<double>> m_vinjPayload; 
 
    DBObjPtr<CDCDedxInjectionTime> m_DBInjectTime; 
  };

} // namespace Belle2