TOPLocalCalFitter.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 <array>
#include <string>
#include <calibration/CalibrationAlgorithm.h>
 
class TH1;
 
namespace Belle2 {
  namespace TOP {

    class TOPLocalCalFitter : public CalibrationAlgorithm {
    public:

      TOPLocalCalFitter();

      ~TOPLocalCalFitter() override;

      void setMinEntries(int minEntries)
      {
        m_minEntries = minEntries;
      }

      void setOutputFileName(const std::string& output)
      {
        m_output = output;
      }

      void setFitConstraintsFileName(const std::string& fitConstraints)
      {
        m_fitConstraints = fitConstraints;
      }

      void setTTSFileName(const std::string& TTSData)
      {
        m_TTSData = TTSData;
      }

      void setFitMode(const std::string& fitterMode)
      {
        if (fitterMode == "calibration")
          B2INFO("Fitter set to calibration mode");
        else if (fitterMode == "monitoring")
          B2INFO("Fitter set to monitoring mode");
        else if (fitterMode == "MC")
          B2INFO("Fitter set to MC mode");
        else
          B2ERROR("Unknown fitter type " << fitterMode  << ". The valid options are calibration, monitoring or MC");
 
        m_fitterMode = fitterMode;
      }

      void fitInAmpliduteBins(bool isFitInAmplitudeBins)
      {
        m_isFitInAmplitudeBins = isFitInAmplitudeBins;
      }
 
    protected:

      void setupOutputTreeAndFile();

      void loadMCInfoTrees();

      void determineFitStatus();

      void fitChannel(short slot, short channel, TH1* h);

      void fitChannel(short slot, short channel, TH1* h, bool inBins, double frac);

      void fitPulser(TH1*, TH1*);

      void calculateChannelT0();

      EResult calibrate() override;
 
    private:
 
      int m_minEntries = 50; 
      std::string m_output = "laserFitResult.root"; 
      std::string m_fitConstraints =
        "/group/belle2/group/detector/TOP/calibration/MCreferences/LaserMCParameters.root"; 
      std::string m_TTSData =
        "/group/belle2/group/detector/TOP/calibration/MCreferences/TTSParametrization.root"; 
      std::string m_fitterMode = "calibration";
      bool m_isFitInAmplitudeBins = false; 
      std::vector<float> m_binEdges = {50, 100, 150, 200, 250, 300, 350, 400, 500, 600, 800, 1000, 1500, 2000}; 
      TFile* m_inputTTS = nullptr; 
      TFile* m_inputConstraints = nullptr; 
      TTree* m_treeTTS = nullptr; 
      TTree* m_treeConstraints =
        nullptr; 
      TFile* m_histFile = nullptr; 
      TTree* m_fitTree = nullptr; 
      TTree* m_timewalkTree =
        nullptr; 
 
      bool m_detectCrosstalk = false; 
 
      TTree* m_crosstalkTree   = nullptr; 
      TTree* m_fitTree_noXtalk = nullptr; 
 
      // Per-candidate (crosstalk) variables
      short m_sl0 = -1; 
      short m_sl1 = -1; 
      short m_ch0 = -1; 
      short m_ch1 = -1; 
      float m_ht0  = NAN; 
      float m_ht1  = NAN; 
      float m_a0   = NAN; 
      float m_a1   = NAN; 
      float m_w0   = NAN; 
      float m_w1   = NAN; 
      float m_q0   = NAN; 
      float m_q1   = NAN; 
      float m_f_q0 = NAN; 
 
      // Variables for the TTS parametrization tree
      float m_mean2 = 0;  
      float m_sigma1 = 0; 
      float m_sigma2 = 0; 
      float m_f1 = 0; 
      float m_f2 = 0; 
      short m_pixelRow = 0; 
      short m_pixelCol = 0; 
 
      // Variables for the MC truth infos
      float m_peakTimeConstraints = 0; 
      float m_deltaTConstraints = 0; 
      float m_fractionConstraints = 0; 
      float m_timeExtraConstraints = 0; 
      float m_sigmaExtraConstraints = 0; 
      float m_alphaExtraConstraints = 0.; 
      float m_nExtraConstraints = 0.; 
      float m_timeBackgroundConstraints = 0.; 
      float m_sigmaBackgroundConstraints = 0.; 
 
 
      // Variables for the output tree
      float m_binLowerEdge = 0; 
      float m_binUpperEdge = 0; 
      short m_channel = 0; 
      short m_slot = 0; 
      short m_row = 0; 
      short m_col = 0; 
      float m_peakTime = 0; 
      float m_deltaT =
        0; 
      float m_sigma = 0.; 
      float m_fraction = 0.; 
      float m_yieldLaser = 0.; 
      float m_histoIntegral = 0.; 
 
      float m_peakTimeErr = 0; 
      float m_deltaTErr = 0; 
      float m_sigmaErr = 0.; 
      float m_fractionErr = 0.; 
      float m_yieldLaserErr = 0.; 
 
      float m_timeExtra = 0.; 
      float m_sigmaExtra = 0.; 
      float m_yieldLaserExtra = 0.; 
      float m_alphaExtra = 0.; 
      float m_nExtra = 0.; 
 
      float m_timeBackground = 0.; 
      float m_sigmaBackground = 0.; 
      float m_yieldLaserBackground = 0.; 
 
      float m_fractionMC = 0.; 
      float m_deltaTMC = 0.; 
      float m_peakTimeMC =
        0.; 
 
      float m_chi2 = 0; 
      float m_rms = 0; 
 
      float m_channelT0 =
        0.; 
      float m_channelT0Err = 0.; 
 
      float m_firstPulserTime = 0.; 
      float m_firstPulserSigma = 0.; 
 
      float m_secondPulserTime =
        0.; 
      float m_secondPulserSigma = 0.;  
 
      short m_fitStatus = 1; 
 
      double m_width = 0; 
      double m_amplitude = 0; 
 
      // Row/Col per (slot, channel). Slots are 0..15, channels 0..511 (0-based inside code).
      std::array<std::array<short, 512>, 16> m_rowOf{}; 
      std::array<std::array<short, 512>, 16> m_colOf{}; 
      bool m_hasChannelMaps{false}; 

      inline short rowOf(short slot, short ch) const noexcept
      {
        return (slot >= 0 && slot < 16 && ch >= 0 && ch < 512) ? m_rowOf[slot][ch] : short(-1);
      }

      inline short colOf(short slot, short ch) const noexcept
      {
        return (slot >= 0 && slot < 16 && ch >= 0 && ch < 512) ? m_colOf[slot][ch] : short(-1);
      }

      inline bool areNeighbors(short slot, short a, short b, int drMax = 1, int dcMax = 1) const noexcept
      {
        const int dr = std::abs(rowOf(slot, a) - rowOf(slot, b));
        const int dc = std::abs(colOf(slot, a) - colOf(slot, b));
        return (dr + dc > 0) && (dr <= drMax) && (dc <= dcMax);
      }

      void buildChannelMaps();
 
    };
  } // namespace TOP
} // namespace Belle2