KLMTimeAlgorithm.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
 
/* KLM headers. */
#include <klm/bklm/geometry/GeometryPar.h>
#include <klm/dataobjects/KLMChannelIndex.h>
#include <klm/dataobjects/KLMElementNumbers.h>
#include <klm/dbobjects/KLMTimeCableDelay.h>
#include <klm/dbobjects/KLMTimeConstants.h>
#include <klm/dbobjects/KLMTimeResolution.h>
#include <klm/dbobjects/KLMEventT0HitResolution.h>
#include <klm/eklm/geometry/GeometryData.h>
 
/* Basf2 headers. */
#include <calibration/CalibrationAlgorithm.h>
 
/* ROOT headers. */
#include <Math/MinimizerOptions.h>
#include <TF1.h>
#include <TGraphErrors.h>
#include <TH1F.h>
#include <TH1I.h>
#include <TH2F.h>
#include <TProfile.h>
 
/* C++ STL headers. */
#include <functional>
#include <map>
#include <utility>
#include <vector>
 
class TFile;  // forward declaration
 
namespace Belle2 {

  class KLMTimeAlgorithm : public CalibrationAlgorithm {
 
  public:

    struct Event {

      double Run = 0;

      double Events = 0;

      double nTrack = 0;

      double Track_Charge = 0;

      double t0 = 0;

      double t0_uc = 0;

      double flyTime = 0;

      double recTime = 0;

      double dist = 0;

      double diffDistX = 0;

      double diffDistY = 0;

      double diffDistZ = 0;

      double eDep = 0;

      double nPE = 0;

      int channelId = 0;

      bool inRPC = 0;

      bool isFlipped = 0;

      bool isGood = 0;

      uint16_t getADCcount = 0;

      double time() const
      {
        return recTime - flyTime;
      }
 
    };

    enum ChannelCalibrationStatus {
 
      /* Not enough data. */
      c_NotEnoughData = 0,
 
      /* Failed fit. */
      c_FailedFit = 1,
 
      /* Successful calibration. */
      c_SuccessfulCalibration = 2,
 
    };

    KLMTimeAlgorithm();

    ~KLMTimeAlgorithm();

    void setDebug()
    {
      m_debug = true;
    }

    void setMC(bool mc)
    {
      m_mc = mc;
    }

    void useEvtT0()
    {
      m_useEventT0 = true;
    }

    void setMinimalDigitNumber(int minimalDigitNumber)
    {
      m_MinimalDigitNumber = minimalDigitNumber;
    }

    void setLowerLimit(int counts)
    {
      m_lower_limit_counts = counts;
    }

    void setApplyChargeRestriction(bool apply)
    {
      m_applyChargeRestriction = apply;
    }

    void setSaveAllPlots(bool on)
    {
      m_saveAllPlots = on;
    }

    void setSaveChannelHists(bool on)
    {
      m_saveChannelHists = on;
    }

    void saveHist();

    double esti_timeShift(const KLMChannelIndex& klmChannel);

    std::pair<int, double> tS_upperStrip(const KLMChannelIndex& klmChannel);

    std::pair<int, double> tS_lowerStrip(const KLMChannelIndex& klmChannel);

    double esti_timeRes(const KLMChannelIndex& klmChannel);

    std::pair<int, double> tR_upperStrip(const KLMChannelIndex& klmChannel);

    std::pair<int, double> tR_lowerStrip(const KLMChannelIndex& klmChannel);
 
 
  protected:

    virtual EResult calibrate() override;
 
  private:

    void setupDatabase();

    CalibrationAlgorithm::EResult readCalibrationData();

    void readCalibrationDataCounts(std::map<KLMChannelNumber, unsigned int>& eventCounts);

    void readCalibrationDataFor2DFit(
      const std::vector<std::pair<KLMChannelNumber, unsigned int>>& channelsBKLM,
      const std::vector<std::pair<KLMChannelNumber, unsigned int>>& channelsEKLM);

    void readCalibrationDataBatch(std::function<bool(const KLMChannelIndex&)> channelFilter);

    void createHistograms();

    void writeThenDelete_(TH1* h, bool write, TDirectory* dir = nullptr);

    void writeThenDelete_(TH2* h, bool write, TDirectory* dir = nullptr);

    bool passesADCCut(const Event& event, int subdetector, int layer) const;

    void fillTimeDistanceProfiles(
      TProfile* profileRpcPhi, TProfile* profileRpcZ,
      TProfile* profileBKLMScintillatorPhi, TProfile* profileBKLMScintillatorZ,
      TProfile* profileEKLMScintillatorPlane1,
      TProfile* profileEKLMScintillatorPlane2, bool fill2dHistograms);

    void timeDistance2dFit(
      const std::vector< std::pair<KLMChannelNumber, unsigned int> >& channels,
      double& delay, double& delayError);

    std::map<KLMChannelNumber, std::vector<struct Event> > m_evts;

    std::map<KLMChannelNumber, int> m_cFlag;

    std::map<KLMChannelNumber, double> m_timeShift;

    std::map<KLMChannelNumber, double> m_timeRes;

    std::map<KLMChannelNumber, double> m_time_channel;

    std::map<KLMChannelNumber, double> m_etime_channel;

    std::map<KLMChannelNumber, double> m_ctime_channel;

    std::map<KLMChannelNumber, double> mc_etime_channel;

    double m_LowerTimeBoundaryRPC = -10.0;

    double m_UpperTimeBoundaryRPC = 10.0;

    double m_LowerTimeBoundaryScintillatorsBKLM = 20.0;

    double m_UpperTimeBoundaryScintillatorsBKLM = 70.0;

    double m_LowerTimeBoundaryScintillatorsEKLM = 20.0;

    double m_UpperTimeBoundaryScintillatorsEKLM = 70.0;

    double m_LowerTimeBoundaryCalibratedRPC = -40.0;

    double m_UpperTimeBoundaryCalibratedRPC = 40.0;

    bool m_useFixedRPCDelay = false;

    double m_fixedRPCDelay = 0.0667;

    double m_LowerTimeBoundaryCalibratedScintillatorsBKLM = -40.0;

    double m_UpperTimeBoundaryCalibratedScintillatorsBKLM = 40.0;

    double m_LowerTimeBoundaryCalibratedScintillatorsEKLM = -40.0;

    double m_UpperTimeBoundaryCalibratedScintillatorsEKLM = 40.0;

    double m_time_channelAvg_rpc = 0.0;

    double m_etime_channelAvg_rpc = 0.0;

    double m_time_channelAvg_scint = 0.0;

    double m_etime_channelAvg_scint = 0.0;

    double m_time_channelAvg_scint_end = 0.0;

    double m_etime_channelAvg_scint_end = 0.0;

    double m_ctime_channelAvg_rpc = 0.0;

    double mc_etime_channelAvg_rpc = 0.0;

    double m_ctime_channelAvg_scint = 0.0;

    double mc_etime_channelAvg_scint = 0.0;

    double m_ctime_channelAvg_scint_end = 0.0;

    double mc_etime_channelAvg_scint_end = 0.0;

    int m_MinimalDigitNumber = 100000000;

    int m_lower_limit_counts = 50;

    const KLMElementNumbers* m_ElementNumbers;

    const bklm::GeometryPar* m_BKLMGeometry = nullptr;

    const EKLM::GeometryData* m_EKLMGeometry = nullptr;

    KLMChannelIndex m_klmChannels;

    ROOT::Math::MinimizerOptions m_minimizerOptions;

    KLMTimeConstants* m_timeConstants = nullptr;

    KLMTimeCableDelay* m_timeCableDelay = nullptr;

    KLMTimeResolution* m_timeResolution = nullptr;

    KLMEventT0HitResolution* m_eventT0HitResolution = nullptr;

    bool m_debug = false;

    bool m_mc = false;

    bool m_useEventT0 = true;

    bool m_applyChargeRestriction = false;

    TH1I* h_calibrated = nullptr;

    TH1I* hc_calibrated = nullptr;

    TH1F* h_diff = nullptr;
 
    /* Monitor graphs of peak value of time distribution for each channel. */

    TGraphErrors* gre_time_channel_rpc = nullptr;

    TGraphErrors* gre_time_channel_scint = nullptr;

    TGraphErrors* gre_time_channel_scint_end = nullptr;
 
    /* Monitor graphs of peak value of calibrated time distribution for each channel. */

    TGraphErrors* gre_ctime_channel_rpc = nullptr;

    TGraphErrors* gre_ctime_channel_scint = nullptr;

    TGraphErrors* gre_ctime_channel_scint_end = nullptr;
 
    /* Monitor graphs of calibration constant value of each channel. */

    TGraph* gr_timeShift_channel_rpc = nullptr;

    TGraph* gr_timeShift_channel_scint = nullptr;

    TGraph* gr_timeShift_channel_scint_end = nullptr;
 
    /* Monitor graphs of calibrated time resolution value of each channel. */

    TGraph* gr_timeRes_channel_rpc = nullptr;

    TGraph* gr_timeRes_channel_scint = nullptr;

    TGraph* gr_timeRes_channel_scint_end = nullptr;
 
    /* Profiles used for effective light speed estimation. */

    TProfile* m_ProfileRpcPhi = nullptr;

    TProfile* m_ProfileRpcZ = nullptr;

    TProfile* m_ProfileBKLMScintillatorPhi = nullptr;

    TProfile* m_ProfileBKLMScintillatorZ = nullptr;

    TProfile* m_ProfileEKLMScintillatorPlane1 = nullptr;

    TProfile* m_ProfileEKLMScintillatorPlane2 = nullptr;
 
    /* Profiles of time versus distance (after fit). */

    TProfile* m_Profile2RpcPhi = nullptr;

    TProfile* m_Profile2RpcZ = nullptr;

    TProfile* m_Profile2BKLMScintillatorPhi = nullptr;

    TProfile* m_Profile2BKLMScintillatorZ = nullptr;

    TProfile* m_Profile2EKLMScintillatorPlane1 = nullptr;

    TProfile* m_Profile2EKLMScintillatorPlane2 = nullptr;
 
    /*
     * Histograms of global time distribution used for effective light speed
     * estimation.
     */

    TH1F* h_time_rpc_tc = nullptr;

    TH1F* h_time_scint_tc = nullptr;

    TH1F* h_time_scint_tc_end = nullptr;
 
    /* Histograms of global time distribution before calibration. */

    TH1F* h_time_rpc = nullptr;

    TH1F* h_time_scint = nullptr;

    TH1F* h_time_scint_end = nullptr;
 
    /* Histograms of global time distribution after calibration. */

    TH1F* hc_time_rpc = nullptr;

    TH1F* hc_time_scint = nullptr;

    TH1F* hc_time_scint_end = nullptr;
 
    /*
     * Histograms of time distribution for forward (backward)
     * before calibration.
     */

    TH1F* h_timeF_rpc[2] = {nullptr};

    TH1F* h_timeF_scint[2] = {nullptr};

    TH1F* h_timeF_scint_end[2] = {nullptr};
 
    /*
     * Histograms of time distribution for forward (backward)
     * after calibration.
     */

    TH1F* hc_timeF_rpc[2] = {nullptr};

    TH1F* hc_timeF_scint[2] = {nullptr};

    TH1F* hc_timeF_scint_end[2] = {nullptr};
 
    /*
     * Histograms of time dependent on sector for forward (backward)
     * before calibration.
     */

    TH2F* h2_timeF_rpc[2] = {nullptr};

    TH2F* h2_timeF_scint[2] = {nullptr};

    TH2F* h2_timeF_scint_end[2] = {nullptr};
 
    /*
     * Histograms of time dependent on sector for forward (backward)
     * after calibration.
     */

    TH2F* h2c_timeF_rpc[2] = {nullptr};

    TH2F* h2c_timeF_scint[2] = {nullptr};

    TH2F* h2c_timeF_scint_end[2] = {nullptr};
 
    /* Histograms of time distribution for sectors before calibration. */

    TH1F* h_timeFS_rpc[2][8] = {{nullptr}};

    TH1F* h_timeFS_scint[2][8] = {{nullptr}};

    TH1F* h_timeFS_scint_end[2][4] = {{nullptr}};
 
    /* Histograms of time distribution for sectors after calibration. */

    TH1F* hc_timeFS_rpc[2][8] = {{nullptr}};

    TH1F* hc_timeFS_scint[2][8] = {{nullptr}};

    TH1F* hc_timeFS_scint_end[2][4] = {{nullptr}};
 
    /*
     * Histograms of time distribution dependent on layer of sectors
     * before calibration.
     */

    TH2F* h2_timeFS[2][8] = {{nullptr}};

    TH2F* h2_timeFS_end[2][4] = {{nullptr}};
 
    /*
     * Histograms of time distribution dependent on layer of sectors
     * after calibration.
     */

    TH2F* h2c_timeFS[2][8] = {{nullptr}};

    TH2F* h2c_timeFS_end[2][4] = {{nullptr}};
 
    /* Histograms of time distribution of one layer before calibration. */

    TH1F* h_timeFSL[2][8][15] = {{{nullptr}}};

    TH1F* h_timeFSL_end[2][4][14] = {{{nullptr}}};
 
    /* Histograms of time distribution of one layer after calibration. */

    TH1F* hc_timeFSL[2][8][15] = {{{nullptr}}};

    TH1F* hc_timeFSL_end[2][4][14] = {{{nullptr}}};
 
    /* Histograms of time distribution of one plane before calibration. */

    TH1F* h_timeFSLP[2][8][15][2] = {{{{nullptr}}}};

    TH1F* h_timeFSLP_end[2][4][14][2] = {{{{nullptr}}}};
 
    /* Histograms of time distribution of one plane after calibration. */

    TH1F* hc_timeFSLP[2][8][15][2] = {{{{nullptr}}}};

    TH1F* hc_timeFSLP_end[2][4][14][2] = {{{{nullptr}}}};
 
    /*
     * Histograms of time distribution dependent on channels
     * before calibration.
     */

    TH2F* h2_timeFSLP[2][8][15][2] = {{{{nullptr}}}};

    TH2F* h2_timeFSLP_end[2][4][14][2] = {{{{nullptr}}}};
 
    /*
     * Histograms of time distribution dependent on channels
     * after calibration.
     */

    TH2F* h2c_timeFSLP[2][8][15][2] = {{{{nullptr}}}};

    TH2F* h2c_timeFSLP_end[2][4][14][2] = {{{{nullptr}}}};
 
    /* Histograms of time distribution of each channel before calibration. */

    TH1F* h_timeFSLPC_tc[2][8][15][2][54] = {{{{{nullptr}}}}};

    TH1F* h_timeFSLPC[2][8][15][2][54] = {{{{{nullptr}}}}};

    TH2F* m_HistTimeLengthBKLM[2][8][15][2][54] = {{{{{nullptr}}}}};

    TH1F* h_timeFSLPC_tc_end[2][4][14][2][75] = {{{{{nullptr}}}}};

    TH1F* h_timeFSLPC_end[2][4][14][2][75] = {{{{{nullptr}}}}};

    TH2F* m_HistTimeLengthEKLM[2][4][14][2][75] = {{{{{nullptr}}}}};
 
    /* Histograms of time distribution of each channel after calibration. */

    TH1F* hc_timeFSLPC[2][8][15][2][54] = {{{{{nullptr}}}}};

    TH1F* hc_timeFSLPC_end[2][4][14][2][75] = {{{{{nullptr}}}}};

    TH1F* h_eventT0_rpc = nullptr;

    TH1F* h_eventT0_scint = nullptr;

    TH1F* h_eventT0_scint_end = nullptr;

    TH1F* hc_eventT0_rpc = nullptr;

    TH1F* hc_eventT0_scint = nullptr;

    TH1F* hc_eventT0_scint_end = nullptr;

    TH1F* h_nHits_plus_rpc;

    TH1F* h_nHits_minus_rpc;

    TH1F* h_nHits_plus_scint;

    TH1F* h_nHits_minus_scint;

    TH1F* h_nHits_plus_scint_end;

    TH1F* h_nHits_minus_scint_end;

    TH2F* h2_deltaT0_vs_v_rpc;

    TH2F* h2_deltaT0_vs_v_scint;

    TH2F* h2_deltaT0_vs_v_scint_end;

    TProfile* prof_deltaT0_rms_vs_v_rpc;

    TProfile* prof_deltaT0_rms_vs_v_scint;

    TProfile* prof_deltaT0_rms_vs_v_scint_end;

    TH2F* h2_deltaT0_vs_nhits_rpc;

    TH2F* h2_deltaT0_vs_nhits_scint;

    TH2F* h2_deltaT0_vs_nhits_scint_end;

    TH1F* hc_eventT0_rpc_lowN;

    TH1F* hc_eventT0_rpc_midN;

    TH1F* hc_eventT0_rpc_highN;

    TH1F* hc_eventT0_scint_lowN;

    TH1F* hc_eventT0_scint_midN;

    TH1F* hc_eventT0_scint_highN;

    TH1F* hc_eventT0_scint_end_lowN;

    TH1F* hc_eventT0_scint_end_midN;

    TH1F* hc_eventT0_scint_end_highN;
 
    /* Formulas used for fitting. */

    TF1* fcn_pol1 = nullptr;

    TF1* fcn_const = nullptr;

    TF1* fcn_gaus = nullptr;

    TF1* fcn_land = nullptr;

    TFile* m_outFile = nullptr;

    bool m_saveAllPlots = false;

    bool m_saveChannelHists = false;

    TDirectory* m_channelHistDir_BKLM[2][8][15][2] = {};

    TDirectory* m_channelHistDir_EKLM[2][4][14][2] = {};
 
  };

}