TOPTimeBaseCalibratorModule.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/core/Module.h>
#include <framework/datastore/StoreArray.h>
#include <top/dataobjects/TOPDigit.h>
 
#include <string>
#include <vector>
#include <TMatrixDSym.h>
#include <TH1F.h>
#include <TH2F.h>
 
namespace Belle2 {
  struct Hit {
    double time = 0; 
    double timeErr = 0; 
    int pulseHeight = 0; 
    double pulseWidth = 0; 
    Hit(double t, double et, int height, double width)
    {
      time = t;
      timeErr = et;
      pulseHeight = height;
      pulseWidth = width;
    }
  };
 
  struct TwoTimes {
    float t1 = 0; 
    float t2 = 0; 
    float sigma = 0; 
    bool good = false; 
    TwoTimes(double time1, double time2, double sig)
    {
      if (time1 < time2) {
        t1 = time1;
        t2 = time2;
      } else {
        t1 = time2;
        t2 = time1;
      }
      sigma = sig;
      good = true;
    }
  };
 
  class TOPTimeBaseCalibratorModule : public Module {
 
  public:
 
    TOPTimeBaseCalibratorModule();
 
    virtual ~TOPTimeBaseCalibratorModule();
 
    virtual void initialize() override;
 
    virtual void beginRun() override;
 
    virtual void event() override;
 
    virtual void endRun() override;
 
    virtual void terminate() override;
 
  private:
 
    enum {c_NumChannels = 512, 
          c_NumBoardstacks = 4, 
          c_NumScrodChannels =  c_NumChannels / c_NumBoardstacks,
          c_WindowSize = 64,   
          c_TimeAxisSize = c_WindowSize * 4 
         };
 
    bool calibrateChannel(std::vector<TwoTimes>& ntuple,
                          unsigned scrodID, unsigned scrodChannel,
                          TH1F& Hchi2, TH1F& Hndf, TH1F& HDeltaT);
 
    bool matrixInversion(const std::vector<TwoTimes>& ntuple,
                         unsigned scrodID, unsigned scrodChannel,
                         double meanTimeDifference,
                         TH1F& Hchi2, TH1F& Hndf, TH1F& HDeltaT);
 
    bool iterativeTBC(const std::vector<TwoTimes>& ntuple,
                      unsigned scrodID, unsigned scrodChannel,
                      double meanTimeDifference,
                      TH1F& Hchi2, TH1F& Hndf, TH1F& HDeltaT);
 
    void Iteration(const std::vector<TwoTimes>& ntuple, std::vector<double>& xval);
 
    double Chisq(const std::vector<TwoTimes>& ntuple, const std::vector<double>& xxval);
 
    void saveAsHistogram(const std::vector<double>& vec,
                         const std::string& name,
                         const std::string& title,
                         const std::string& xTitle = "",
                         const std::string& yTitle = "") const;
 
    void saveAsHistogram(const std::vector<double>& vec,
                         const std::vector<double>& err,
                         const std::string& name,
                         const std::string& title,
                         const std::string& xTitle = "",
                         const std::string& yTitle = "") const;
 
    void saveAsHistogram(const TMatrixDSym& M,
                         const std::string& name,
                         const std::string& title) const;
 
 
    int m_moduleID = 0;    
    double m_minTimeDiff = 0;  
    double m_maxTimeDiff = 0;  
    unsigned m_minHits = 0;    
    unsigned m_numIterations = 100;    
    std::string m_directoryName; 
    unsigned m_method = 0; 
    bool m_useFallingEdge = false; 
    bool m_saveMatrix = false; 
    std::vector<TwoTimes> m_ntuples[c_NumChannels]; 
    double m_syncTimeBase = 0; 
    StoreArray<TOPDigit> m_digits; 
    TH2F m_goodHits;       
    TH2F m_calPulseFirst;  
    TH2F m_calPulseSecond; 
    // parameters for iTBC
    int m_good = 0;   
    double m_dt_min = 20.0; 
    double m_dt_max = 24.0; 
    double m_dev_step = 0.001; 
    double m_xstep = 0.020; 
    double m_dchi2dxv = 0.0;
    double m_change_xstep = 0.015;
    double m_new_xstep = 2.0 * m_xstep; 
    double m_min_binwidth = 0.05; 
    double m_max_binwidth = 2.0; 
    double m_dchi2_min = 0.2; 
    unsigned   m_conv_iter = 100; 
    double m_deltamin = 0.01 * 0.01; 
    double  m_DtSigma = 0.0424; 
    double m_sigm2_exp = m_DtSigma * m_DtSigma;
  };
 
} // Belle2 namespace