SpaceResolutionCalibrationAlgorithm.h

/**************************************************************************
 * BASF2 (Belle Analysis Framework 2)                                     *
 * Copyright(C) 2010 - Belle II Collaboration                             *
 *                                                                        *
 * Author: The Belle II Collaboration                                     *
 * Contributors:  Makoto Uchida                                           *
 *                                                                        *
 * This software is provided "as is" without any warranty.                *
 **************************************************************************/
 
#pragma once
 
#include <TH1F.h>
#include <TH2F.h>
#include <TGraphErrors.h>
#include <calibration/CalibrationAlgorithm.h>
#include <cdc/dbobjects/CDCGeometry.h>
#include <framework/database/DBObjPtr.h>
 
namespace Belle2 {
  namespace CDC {
 
    class SpaceResolutionCalibrationAlgorithm :  public CalibrationAlgorithm {
      //      typedef std::array<float, 3> array3; /**< angle bin info. */
    public:
 
      SpaceResolutionCalibrationAlgorithm();
 
      ~SpaceResolutionCalibrationAlgorithm() {}
 
      void setDebug(bool debug = false) {m_debug = debug; }
 
      void setMinimumNDF(double ndf) {m_minNdf = ndf;}
 
      void setMinimumPval(double pval) {m_minPval = pval;}
 
      void setBinWidth(double bw) {m_binWidth = bw;}
 
      void setBField(bool bfield) {m_bField = bfield;}
 
      void setStoreHisto(bool storeHist = false) {m_storeHisto = storeHist;}
 
      void enableTextOutput(bool output = true) {m_textOutput = output;}
 
      void setOutputFileName(std::string outputname) {m_outputFileName.assign(outputname);}
 
      void setHistFileName(const std::string& name) {m_histName = "histSigma_" + name + ".root";}
 
      void setThreshold(double th = 0.6) {m_threshold = th;}
 
    protected:
      EResult calibrate() override;
      void createHisto();
      void storeHisto();
      void write();
      void prepare();
 
      double getUpperBoundaryForFit(TGraphErrors* graph)
      {
        double ymax = 0;
        double xmax = 0;
        int imax = 0;
        double x, y;
        int unCount = floor(0.05 / m_binWidth);
        int N = graph->GetN();
        int Nstart = floor(0.5 * (N - unCount));
        int Nend = N - unCount;
        for (int i  = Nstart; i < Nend; ++i) {
          graph->GetPoint(i, x, y);
          if (graph->GetErrorY(i) > 0.06E-3) continue;
          if (y > ymax) {
            xmax = x; ymax = y;
            imax = i;
          }
        }
        if (imax <= Nstart) {
          graph->GetPoint(Nend, x, y);
          xmax = x;
        }
        return xmax;
      }
 
    private:
      static const int Max_nalpha = 18; 
      static const int Max_ntheta = 7; 
      static const unsigned short Max_np = 40; 
      double m_minNdf = 5; 
      double m_minPval = 0.; 
      double m_binWidth = 0.05; 
      bool m_debug = false;   
      bool m_storeHisto = false; 
      bool m_bField = true;                   
      double m_threshold = 0.6 ; 
      double m_sigma[56][2][18][7][8]; 
      TGraphErrors* m_gFit[56][2][18][7];  
      TGraphErrors* m_graph[56][2][18][7];    
      TH2F* m_hBiased[56][2][Max_nalpha][Max_ntheta]; 
      TH2F* m_hUnbiased[56][2][Max_nalpha][Max_ntheta]; 
      TH1F* m_hMeanBiased[56][2][Max_nalpha][Max_ntheta]; 
      TH1F* m_hSigmaBiased[56][2][Max_nalpha][Max_ntheta]; 
      TH1F* m_hMeanUnbiased[56][2][Max_nalpha][Max_ntheta]; 
      TH1F* m_hSigmaUnbiased[56][2][Max_nalpha][Max_ntheta]; 
      int m_fitStatus[56][2][Max_nalpha][Max_ntheta] = {{{{0}}}} ; 
      int m_nAlphaBins; 
      int m_nThetaBins;
      float m_lowerAlpha[18];
      float m_upperAlpha[18];
      float m_iAlpha[18]; 
      float m_lowerTheta[7]; 
      float m_upperTheta[7];
      float m_iTheta[7]; 
      unsigned short m_sigmaParamMode = 0; 
      double m_sigmaPost[56][2][18][7][8]; 
      unsigned short m_sigmaParamModePost; 
      bool  m_textOutput = false; 
      std::string m_outputFileName = "sigma_new.dat"; 
      std::string m_histName = "histSigma.root"; 
      DBObjPtr<CDCGeometry> m_cdcGeo; 
    };
  }
}