// -----------------------------------------------------------------------------
// ----- calib_digis.C                                                     -----
// -----                                                                   -----
// ----- adapted by C. Simon on 2020-02-12                                 -----
// -----                                                                   -----
// -----------------------------------------------------------------------------

#include <iostream>
#include <map>
#include <tuple>
#include <set>

#include <boost/regex.hpp>
#include <boost/lexical_cast.hpp>

#include "TStyle.h"
#include "TSystem.h"
#include "TFile.h"
#include "TTree.h"
#include "TClonesArray.h"
#include "TH1.h"
#include "TH2.h"
#include "TProfile.h"
#include "TCanvas.h"
#include "TMath.h"
#include "TF1.h"
#include "TGeoManager.h"
#include "TGeoPhysicalNode.h"
#include "TStopwatch.h"
#include "TPaveStats.h"
#include "TRandom3.h"

#include "FairLogger.h"
#include "FairRootFileSink.h"
#include "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"

#include "CbmTofDigiExp.h"
#include "tof/TofTools/CbmTofGeoHandler.h"


struct TCounter;


const char* gcInputPath;
const char* gcInputFileName;
const char* gcOutputFileName;
const char* gcCalibFileName;
const char* gcGeoVersion;
Int_t giBRefModuleType;
Int_t giBRefModuleIndex;
Int_t giBRefCounterIndex;
Bool_t gbBRefPadReadout;
Int_t giBRefCalibRefModuleType;
Int_t giBRefCalibRefModuleIndex;
Int_t giBRefCalibRefCounterIndex;
UInt_t guBRefAddress;
UInt_t guBRefCalibRefAddress;
std::map<UInt_t, TCounter*> gAllCounterMap;
std::map<UInt_t, TCounter*> gCalibCounterMap;
std::set<UInt_t> gBRefAddressSet;
CbmTofGeoHandler* gGeoHandler(NULL);
Int_t giCurrentModuleType(-1);
Int_t giCurrentModuleIndex(-1);
Int_t giCurrentCounterIndex(-1);
TString gCurrentNodePath("");
TString gTofNodePath("");
TString gCurrentBoxNodePath("");
TClonesArray* gInputDigiArray(NULL);
TBranch* gInputDigiBranch(NULL);
Int_t giNEvents(0);
Bool_t gbNanoSeconds(kTRUE);
Bool_t gbPreCalibration(kFALSE);

const Double_t gdBRefMeanPosition = 0.;
const Double_t gdDefaultSignalVelocity = 18.;
const Double_t gdToTExpectationValue = 2.;
const Double_t gdToTStandardDeviation = 0.5;
const Double_t gdTriggerWindow = 600.;
const Int_t giMaxNCells = 32;
const Int_t giNWalkBinsX = 1000;
const Double_t gdMinimumToT = 0.001;
const Double_t gdMaximumToT = 100.;
const Double_t gdDigiToTOffset = 20.;
const Double_t gdWalkBinWidth = (gdMaximumToT - gdMinimumToT)/static_cast<Double_t>(giNWalkBinsX);
const Double_t gdMeanTimeResidual = 0.;
const Double_t gdMinimumTimeResidual = -5.;
const Double_t gdMaximumTimeResidual = 15.;
const Int_t giNTimeResidualBins = 211; // 2111; TODO
const Int_t giNBRefWalkCorrections = 45;
const Int_t giNWalkCorrections = 30;

const Bool_t gbToTScalingBeforeCalibration = kFALSE;
const Bool_t gbParallelSlewingCorrections = kTRUE;
const Bool_t gbUseInputSignalPropagationVelocity = kTRUE;
const Bool_t gbSkipCounterSlewingCorrections = kTRUE;
const Bool_t gbDebugFlag = kTRUE; // not in use
const Bool_t gbPreserveHistograms = kTRUE;

TF1* gTimeDifferenceFitNH(NULL);


std::map<std::tuple<Int_t, Int_t, Int_t>, std::tuple<Int_t, UInt_t, Bool_t, Bool_t, Double_t>> gCounterInformation
{
//                                           ref   10987654321098765432109876543210  walk?  BRef?
  {std::make_tuple(5, 2, 0), std::make_tuple(10, 0b00000000000000000000000000000000, kTRUE, kTRUE, 18.000)},
  {std::make_tuple(5, 1, 0), std::make_tuple( 9, 0b00000000000000000000000000000000, kTRUE, kTRUE, 18.000)},
  {std::make_tuple(3, 0, 0), std::make_tuple( 1, 0b10000000000000000000000000000001, kFALSE, kFALSE, 16.684/*17.275*/)},
  {std::make_tuple(9, 0, 0), std::make_tuple( 1, 0b00000000100000000000000000000000, kFALSE, kFALSE, 16.800/*16.980*/)},
  {std::make_tuple(7, 0, 0), std::make_tuple( 1, 0b00000000000000001000000000000001, kFALSE, kFALSE, 16.195)},
  {std::make_tuple(4, 0, 0), std::make_tuple( 7, 0b00000000000000000000001000000001, kFALSE, kFALSE, 14.539)}
};


struct TCounter
{
  Int_t fiModuleType;
  Int_t fiModuleIndex;
  Int_t fiCounterIndex;
  TGeoNode* fCounterNode;
  TGeoPhysicalNode* fCentralPN;
  TGeoPhysicalNode* fCounterBoxPN;
  Int_t fiNCells;
  Double_t fdCounterXWidth;
  Double_t fdCounterYWidth;
  Double_t fdCounterZWidth;
  Double_t fdCounterAcrossWidth;
  Double_t fdCounterAlongWidth;
  Double_t fdCellXWidth;
  Double_t fdCellYWidth;
  Double_t fdCellAcrossWidth;
  Double_t fdCellAlongWidth;
  Bool_t fbCellsAlongX;
  Bool_t fbCellIndexingAlongAxis;
  Int_t fiNGlassPlates;
  Int_t fiNGasGaps;
  Double_t fdGlassThickness;
  Double_t fdGapSize;
  Double_t fbPadReadout;
  Double_t fdSignalPropagationDelay[giMaxNCells][2];
  Double_t fdSlewingCorrection[giMaxNCells][2][giNWalkBinsX + 2];
  Double_t fdSignalPropagationVelocity[giMaxNCells];
  Double_t fdToTScaling[giMaxNCells][2];
  Double_t fdToTOffset[giMaxNCells][2];
  Double_t fdMeanSignalPropagationVelocity;
  Double_t fdInputSignalPropagationVelocity;
  Bool_t fbDeadCell[giMaxNCells];
  Bool_t fbIsBRef;
  Int_t fiReferenceCell;
  Bool_t fbCorrectSlewing;
  Bool_t fbDigiAvailable[giMaxNCells][2];
  CbmTofDigiExp* fDigi[giMaxNCells][2];
  Double_t fdDigiSlewing[giMaxNCells][2];
  TH1* fTimeResidual[giMaxNCells];
  TH1* fTimeDifference[giMaxNCells];
  TH2* fTimeResidualWalk[giMaxNCells][2];
  TH1* fTimeMeanCorrelation[giMaxNCells];
  TH1* fTimeDifferenceCorrelation[giMaxNCells];
  TF1* fTimeDifferenceFitFunction[giMaxNCells];
  TH2* fAllTimeResidual;
  TH2* fAllTimeMeanCorrelation;
  TH2* fAllTimeDifferenceCorrelation;
  TH2* fAllTimeDifference;
  TH2* fAllLocalPosition;
  TH2* fAllLocalPositionMeanVelocity;
  TProfile* fSignalPropagationVelocity;
  TH2* fToTSpectra[2];
  TH2* fToTCorrelation[giMaxNCells];
  TH2* fAllToTSpectra;

  TCounter(Int_t iModuleType, Int_t iModuleIndex, Int_t iCounterIndex,
           TGeoNode* tCounterNode, TGeoPhysicalNode* tCentralPN, TGeoPhysicalNode* tCounterBoxPN) :
    fiModuleType(iModuleType),
    fiModuleIndex(iModuleIndex),
    fiCounterIndex(iCounterIndex),
    fCounterNode(tCounterNode),
    fCentralPN(tCentralPN),
    fCounterBoxPN(tCounterBoxPN),
    fiNCells(0),
    fdCounterXWidth(0.),
    fdCounterYWidth(0.),
    fdCounterZWidth(0.),
    fdCounterAcrossWidth(0.),
    fdCounterAlongWidth(0.),
    fdCellXWidth(0.),
    fdCellYWidth(0.),
    fdCellAcrossWidth(0.),
    fdCellAlongWidth(0.),
    fbCellsAlongX(kFALSE),
    fbCellIndexingAlongAxis(kFALSE),
    fiNGlassPlates(0),
    fiNGasGaps(0),
    fdGlassThickness(0.),
    fdGapSize(0.),
    fbPadReadout(kFALSE),
    fdSignalPropagationDelay{0.},
    fdSlewingCorrection{0.},
    fdSignalPropagationVelocity{0.},
    fdToTScaling{1.},
    fdToTOffset{0.},
    fdMeanSignalPropagationVelocity(gdDefaultSignalVelocity),
    fdInputSignalPropagationVelocity(gdDefaultSignalVelocity),
    fbDeadCell{kFALSE},
    fbIsBRef(kFALSE),
    fiReferenceCell(-1),
    fbCorrectSlewing(kTRUE),
    fbDigiAvailable{kFALSE},
    fDigi{NULL},
    fdDigiSlewing{0.},
    fTimeResidual{NULL},
    fTimeDifference{NULL},
    fTimeResidualWalk{NULL},
    fTimeMeanCorrelation{NULL},
    fTimeDifferenceCorrelation{NULL},
    fTimeDifferenceFitFunction{NULL},
    fAllTimeResidual(NULL),
    fAllTimeMeanCorrelation(NULL),
    fAllTimeDifferenceCorrelation(NULL),
    fAllTimeDifference(NULL),
    fAllLocalPosition(NULL),
    fAllLocalPositionMeanVelocity(NULL),
    fSignalPropagationVelocity(NULL),
    fToTSpectra{NULL},
    fToTCorrelation{NULL},
    fAllToTSpectra(NULL)
  {
    Double_t dAxisLow(0.);
    Double_t dAxisHigh(0.);

    for(Int_t iComponent = 0; iComponent < fCounterNode->GetNdaughters(); iComponent++)
    {
      if(TString(fCounterNode->GetDaughter(iComponent)->GetName()).Contains("glass", TString::kIgnoreCase))
      {
        TGeoNode* tPlateNode = fCounterNode->GetDaughter(iComponent);

        TGeoVolume* tPlateVolume = tPlateNode->GetVolume();
        TGeoShape* tPlateShape = tPlateVolume->GetShape();

        fdCounterZWidth += tPlateShape->GetAxisRange(3, dAxisLow, dAxisHigh);

        if(0 == fiNGlassPlates)
        {
          fdGlassThickness = tPlateShape->GetAxisRange(3, dAxisLow, dAxisHigh);
        }

        fiNGlassPlates++;
      }
      else if(TString(fCounterNode->GetDaughter(iComponent)->GetName()).Contains("gap", TString::kIgnoreCase))
      {
        TGeoNode* tGapNode = fCounterNode->GetDaughter(iComponent);

        TGeoVolume* tGapVolume = tGapNode->GetVolume();
        TGeoShape* tGapShape = tGapVolume->GetShape();

        fdCounterZWidth += tGapShape->GetAxisRange(3, dAxisLow, dAxisHigh);

        if(0 == fiNGasGaps)
        {
          fdCounterXWidth = tGapShape->GetAxisRange(1, dAxisLow, dAxisHigh);
          fdCounterYWidth = tGapShape->GetAxisRange(2, dAxisLow, dAxisHigh);
          fdGapSize = tGapShape->GetAxisRange(3, dAxisLow, dAxisHigh);

          fiNCells = tGapNode->GetNdaughters();
          if(!fiNCells)
          {
            cout<<"MRPC gas gap is not divided into readout cells."<<endl;
          }

          TGeoNode* tFirstCellNode = tGapNode->GetDaughter(0);
          TGeoVolume* tFirstCellVolume = tFirstCellNode->GetVolume();
          TGeoShape* tFirstCellShape = tFirstCellVolume->GetShape();

          fdCellXWidth = tFirstCellShape->GetAxisRange(1, dAxisLow, dAxisHigh);
          fdCellYWidth = tFirstCellShape->GetAxisRange(2, dAxisLow, dAxisHigh);


          if(fdCellXWidth == fdCounterXWidth)
          {
            fbCellsAlongX = kTRUE;

            fdCellAcrossWidth = fdCellYWidth;
            fdCellAlongWidth = fdCellXWidth;

            fdCounterAcrossWidth = fdCounterYWidth;
            fdCounterAlongWidth = fdCounterXWidth;
          }
          else
          {
            fdCellAcrossWidth = fdCellXWidth;
            fdCellAlongWidth = fdCellYWidth;

            fdCounterAcrossWidth = fdCounterXWidth;
            fdCounterAlongWidth = fdCounterYWidth;
          }

          if(fdCellXWidth <= fdCellYWidth)
          {
            if(fdCellXWidth/fdCellYWidth >= 0.5)
            {
              fbPadReadout = kTRUE;
            }
          }
          else
          {
            if(fdCellXWidth/fdCellYWidth < 2.)
            {
              fbPadReadout = kTRUE;
            }
          }

          TGeoNode* tLastCellNode = tGapNode->GetDaughter(fiNCells - 1);
          TGeoVolume* tLastCellVolume = tLastCellNode->GetVolume();
          TGeoShape* tLastCellShape = tLastCellVolume->GetShape();

          Double_t dCellMasterCoordinates[3] = {0.,0.,0.};
          Double_t dCellLocalCoordinates[3] = {0.,0.,0.};

          if(fbCellsAlongX)
          {
            Double_t dFirstCellLocalYLow(0.);
            Double_t dFirstCellLocalYHigh(0.);
            Double_t dFirstCellMasterYLow(0.);
            Double_t dFirstCellMasterYHigh(0.);

            tFirstCellShape->GetAxisRange(2, dFirstCellLocalYLow, dFirstCellLocalYHigh);

            dCellLocalCoordinates[1] = dFirstCellLocalYLow;
            tFirstCellNode->GetMatrix()->LocalToMaster(dCellLocalCoordinates, dCellMasterCoordinates);
            dFirstCellMasterYLow = dCellMasterCoordinates[1];

            dCellLocalCoordinates[1] = dFirstCellLocalYHigh;
            tFirstCellNode->GetMatrix()->LocalToMaster(dCellLocalCoordinates, dCellMasterCoordinates);
            dFirstCellMasterYHigh = dCellMasterCoordinates[1];

            Double_t dLastCellLocalYLow(0.);
            Double_t dLastCellLocalYHigh(0.);
            Double_t dLastCellMasterYLow(0.);
            Double_t dLastCellMasterYHigh(0.);

            tLastCellShape->GetAxisRange(2, dLastCellLocalYLow, dLastCellLocalYHigh);

            dCellLocalCoordinates[1] = dLastCellLocalYLow;
            tLastCellNode->GetMatrix()->LocalToMaster(dCellLocalCoordinates, dCellMasterCoordinates);
            dLastCellMasterYLow = dCellMasterCoordinates[1];

            dCellLocalCoordinates[1] = dLastCellLocalYHigh;
            tLastCellNode->GetMatrix()->LocalToMaster(dCellLocalCoordinates, dCellMasterCoordinates);
            dLastCellMasterYHigh = dCellMasterCoordinates[1];
            
            if(dFirstCellMasterYLow < 0. && dLastCellMasterYHigh > 0.)
            {
              fbCellIndexingAlongAxis = kTRUE;
            }
            else if(dLastCellMasterYLow < 0. && dFirstCellMasterYHigh > 0.)
            {
              fbCellIndexingAlongAxis = kFALSE;
            }
            else
            {
              cout<<"Incomplete understanding of geometry. Cannot determine direction of cell indexing."<<endl;
            }
          }
          else
          {
            Double_t dFirstCellLocalXLow(0.);
            Double_t dFirstCellLocalXHigh(0.);
            Double_t dFirstCellMasterXLow(0.);
            Double_t dFirstCellMasterXHigh(0.);

            tFirstCellShape->GetAxisRange(1, dFirstCellLocalXLow, dFirstCellLocalXHigh);

            dCellLocalCoordinates[0] = dFirstCellLocalXLow;
            tFirstCellNode->GetMatrix()->LocalToMaster(dCellLocalCoordinates, dCellMasterCoordinates);
            dFirstCellMasterXLow = dCellMasterCoordinates[0];

            dCellLocalCoordinates[0] = dFirstCellLocalXHigh;
            tFirstCellNode->GetMatrix()->LocalToMaster(dCellLocalCoordinates, dCellMasterCoordinates);
            dFirstCellMasterXHigh = dCellMasterCoordinates[0];

            Double_t dLastCellLocalXLow(0.);
            Double_t dLastCellLocalXHigh(0.);
            Double_t dLastCellMasterXLow(0.);
            Double_t dLastCellMasterXHigh(0.);

            tLastCellShape->GetAxisRange(1, dLastCellLocalXLow, dLastCellLocalXHigh);

            dCellLocalCoordinates[0] = dLastCellLocalXLow;
            tLastCellNode->GetMatrix()->LocalToMaster(dCellLocalCoordinates, dCellMasterCoordinates);
            dLastCellMasterXLow = dCellMasterCoordinates[0];

            dCellLocalCoordinates[0] = dLastCellLocalXHigh;
            tLastCellNode->GetMatrix()->LocalToMaster(dCellLocalCoordinates, dCellMasterCoordinates);
            dLastCellMasterXHigh = dCellMasterCoordinates[0];

            if(dFirstCellMasterXLow < 0. && dLastCellMasterXHigh > 0.)
            {
              fbCellIndexingAlongAxis = kTRUE;
            }
            else if(dLastCellMasterXLow < 0. && dFirstCellMasterXHigh > 0.)
            {
              fbCellIndexingAlongAxis = kFALSE;
            }
            else
            {
              cout<<"Incomplete understanding of geometry. Cannot determine direction of cell indexing."<<endl;
            }
          }
        }

        fiNGasGaps++;
      }
    }

    CreateHistograms();


    for(Int_t iCell = 0; iCell < fiNCells; iCell++)
    {
      for(Int_t iSide = 0; iSide < 2; iSide++)
      {
        fdSignalPropagationDelay[iCell][iSide] = 0.;

        for(Int_t iWalkBin = 0; iWalkBin < giNWalkBinsX + 2; iWalkBin++)
        {
          fdSlewingCorrection[iCell][iSide][iWalkBin] = 0.;
        }

        fdToTScaling[iCell][iSide] = 1.;
        fdToTOffset[iCell][iSide] = 0.;
      }

      fdSignalPropagationVelocity[iCell] = gdDefaultSignalVelocity;
      fbDeadCell[iCell] = kFALSE;
    }
  }


  void CreateHistograms()
  {
    for(Int_t iCell = 0; iCell < fiNCells; iCell++)
    {
      fTimeResidual[iCell] = new TH1F(Form("tTimeResidual_%d_%d_%d_%d", fiModuleType, fiModuleIndex, fiCounterIndex, iCell), Form("; time residual %d [ns]; ", iCell),
//                                      201*gdTriggerWindow, -gdTriggerWindow, gdTriggerWindow); // TODO: avoid hard coding
                                      (4 == fiModuleType ? 21 : 201)*gdTriggerWindow, -gdTriggerWindow, gdTriggerWindow);
      fTimeResidual[iCell]->SetDirectory(gROOT);

      fTimeDifference[iCell] = new TH1F(Form("tTimeDifference_%d_%d_%d_%d", fiModuleType, fiModuleIndex, fiCounterIndex, iCell), Form("; time difference %d [ns]; ", iCell),
                                             101*gdTriggerWindow, -gdTriggerWindow/2., gdTriggerWindow/2.);
      fTimeDifference[iCell]->SetDirectory(gROOT);

      for(Int_t iSide = 0; iSide < 2; iSide++)
      {
        fTimeResidualWalk[iCell][iSide] = new TH2F(Form("tTimeResidualWalk_%d_%d_%d_%d_%d", fiModuleType, fiModuleIndex, fiCounterIndex, iCell, iSide), Form("; ToT %d-%d [ns]; time residual [ns]", iCell, iSide),
                                                   giNWalkBinsX, gdMinimumToT, gdMaximumToT,
                                                   giNTimeResidualBins, gdMinimumTimeResidual, gdMaximumTimeResidual);
        fTimeResidualWalk[iCell][iSide]->SetDirectory(gROOT);
      }

      fTimeMeanCorrelation[iCell] = new TH1F(Form("tTimeMeanCorrelation_%d_%d_%d_%d", fiModuleType, fiModuleIndex, fiCounterIndex, iCell), Form("; time mean corr %d [ns]; ", iCell),
                                             1001, -5., 5.);
      fTimeMeanCorrelation[iCell]->SetDirectory(gROOT);

      fTimeDifferenceCorrelation[iCell] = new TH1F(Form("tTimeDifferenceCorrelation_%d_%d_%d_%d", fiModuleType, fiModuleIndex, fiCounterIndex, iCell), Form("; time difference corr %d [ns]; ", iCell),
                                                   1001, -5., 5.);
      fTimeDifferenceCorrelation[iCell]->SetDirectory(gROOT);
    }


    fAllTimeResidual = new TH2F(Form("tAllTimeResidual_%d_%d_%d", fiModuleType, fiModuleIndex, fiCounterIndex),
                                  "; cell []; time residual [ns]",
                                  fiNCells, 0., static_cast<Double_t>(fiNCells),
                                  (4 == fiModuleType ? 11 : /*101*/51)*200, -100., 100.);
    fAllTimeResidual->SetDirectory(gROOT);

    fAllTimeMeanCorrelation = new TH2F(Form("tAllTimeMeanCorrelation_%d_%d_%d", fiModuleType, fiModuleIndex, fiCounterIndex),
                                  "; cell []; time mean corr [ns]",
                                  fiNCells, 0., static_cast<Double_t>(fiNCells),
                                  1501/*1001*/, -5., 5.);
    fAllTimeMeanCorrelation->SetDirectory(gROOT);

    fAllTimeDifferenceCorrelation = new TH2F(Form("tAllTimeDifferenceCorrelation_%d_%d_%d", fiModuleType, fiModuleIndex, fiCounterIndex),
                                  "; cell []; time diff corr [ns]",
                                  fiNCells, 0., static_cast<Double_t>(fiNCells),
                                  1501/*1001*/, -5., 5.);
    fAllTimeDifferenceCorrelation->SetDirectory(gROOT);

    fAllTimeDifference = new TH2F(Form("tAllTimeDifference_%d_%d_%d", fiModuleType, fiModuleIndex, fiCounterIndex),
                                  "; cell []; time difference [ns]",
                                  fiNCells, 0., static_cast<Double_t>(fiNCells),
                                  751, -7.5, 7.5);
    fAllTimeDifference->SetDirectory(gROOT);

    fAllLocalPosition = new TH2F(Form("tAllLocalPosition_%d_%d_%d", fiModuleType, fiModuleIndex, fiCounterIndex),
                                  "; cell []; position w/ v_{cell} [cm]",
                                  fiNCells, 0., static_cast<Double_t>(fiNCells),
                                  /*105*/27*TMath::CeilNint(fdCellAlongWidth), -5.*TMath::Ceil(fdCellAlongWidth), 5.*TMath::Ceil(fdCellAlongWidth));
    fAllLocalPosition->SetDirectory(gROOT);

    fAllLocalPositionMeanVelocity = new TH2F(Form("tAllLocalPositionMeanVelocity_%d_%d_%d", fiModuleType, fiModuleIndex, fiCounterIndex),
                                             "; cell []; position w/ v_{mean} [cm]",
                                             fiNCells, 0., static_cast<Double_t>(fiNCells),
                                             /*105*/27*TMath::CeilNint(fdCellAlongWidth), -5.*TMath::Ceil(fdCellAlongWidth), 5.*TMath::Ceil(fdCellAlongWidth));
    fAllLocalPositionMeanVelocity->SetDirectory(gROOT);

    fSignalPropagationVelocity = new TProfile(Form("tSignalPropagationVelocity_%d_%d_%d", fiModuleType, fiModuleIndex, fiCounterIndex),
                                              "; cell []; v_{sig} [cm/ns]",
                                              fiNCells, 0., static_cast<Double_t>(fiNCells),
                                              10., 20.);
    fSignalPropagationVelocity->SetDirectory(gROOT);

    for(Int_t iSide = 0; iSide < 2; iSide++)
    {
      fToTSpectra[iSide] = new TH2F(Form("tToTSpectra_%d_%d_%d_%d", fiModuleType, fiModuleIndex, fiCounterIndex, iSide),
                                    Form("; cell []; ToT side %d [ns]", iSide),
                                    fiNCells, 0., static_cast<Double_t>(fiNCells),
                                    giNWalkBinsX, gdMinimumToT, gdMaximumToT);
      fToTSpectra[iSide]->SetDirectory(gROOT);
    }

    for(Int_t iCell = 0; iCell < fiNCells; iCell++)
    {
      fToTCorrelation[iCell] = new TH2F(Form("tToTCorrelation_%d_%d_%d_%d", fiModuleType, fiModuleIndex, fiCounterIndex, iCell),
                                        Form("; ToT side %d-0 [ns]; ToT side %d-1 [ns]", iCell, iCell),
                                        giNWalkBinsX, gdMinimumToT, gdMaximumToT,
                                        giNWalkBinsX, gdMinimumToT, gdMaximumToT);
      fToTCorrelation[iCell]->SetDirectory(gROOT);
    }

    fAllToTSpectra = new TH2F(Form("tAllToTSpectra_%d_%d_%d", fiModuleType, fiModuleIndex, fiCounterIndex),
                              "; 2*cell + side []; ToT [ns]",
                              2*fiNCells, 0., 2.*static_cast<Double_t>(fiNCells),
                              giNWalkBinsX, gdMinimumToT, gdMaximumToT);
    fAllToTSpectra->SetDirectory(gROOT);
  }


  void ResetHistograms()
  {
    for(Int_t iCell = 0; iCell < fiNCells; iCell++)
    {
      fTimeResidual[iCell]->Reset();
      fTimeResidual[iCell]->GetXaxis()->SetRangeUser(-gdTriggerWindow, gdTriggerWindow);
      
      fTimeDifference[iCell]->Reset();
      fTimeDifference[iCell]->GetXaxis()->SetRangeUser(-gdTriggerWindow/2., gdTriggerWindow/2.);

      for(Int_t iSide = 0; iSide < 2; iSide++)
      {
        fTimeResidualWalk[iCell][iSide]->Reset();
        fTimeResidualWalk[iCell][iSide]->GetYaxis()->SetRangeUser(gdMinimumTimeResidual, gdMaximumTimeResidual);
      }

      fTimeMeanCorrelation[iCell]->Reset();
      fTimeMeanCorrelation[iCell]->GetXaxis()->SetRangeUser(-5., 5.);

      fTimeDifferenceCorrelation[iCell]->Reset();
      fTimeDifferenceCorrelation[iCell]->GetXaxis()->SetRangeUser(-5., 5.);
    }
  }


  void ResetAllHistograms()
  {
    fAllTimeResidual->Reset();
    fAllTimeResidual->GetYaxis()->SetRangeUser(-100., 100.);

    fAllTimeMeanCorrelation->Reset();
    fAllTimeMeanCorrelation->GetYaxis()->SetRangeUser(-5., 5.);

    fAllTimeDifferenceCorrelation->Reset();
    fAllTimeDifferenceCorrelation->GetYaxis()->SetRangeUser(-5., 5.);

    fAllTimeDifference->Reset();
    if(4 == fiModuleType)
    {
      fAllTimeDifference->GetYaxis()->SetRangeUser(-2.5, 2.5);
    }
    else
    {
      fAllTimeDifference->GetYaxis()->SetRangeUser(-7.5, 7.5);
    }

    fAllLocalPosition->Reset();
    fAllLocalPosition->GetYaxis()->SetRangeUser(-5.*TMath::Ceil(fdCellAlongWidth), 5.*TMath::Ceil(fdCellAlongWidth));

    fAllLocalPositionMeanVelocity->Reset();
    fAllLocalPositionMeanVelocity->GetYaxis()->SetRangeUser(-5.*TMath::Ceil(fdCellAlongWidth), 5.*TMath::Ceil(fdCellAlongWidth));
  }


  void ResetToTSpectra()
  {
    for(Int_t iSide = 0; iSide < 2; iSide++)
    {
      fToTSpectra[iSide]->Reset();
      fToTSpectra[iSide]->GetYaxis()->SetRangeUser(gdMinimumToT, gdMaximumToT);
    }

    for(Int_t iCell = 0; iCell < fiNCells; iCell++)
    {
      fToTCorrelation[iCell]->Reset();
      fToTCorrelation[iCell]->GetXaxis()->SetRangeUser(gdMinimumToT, gdMaximumToT);
      fToTCorrelation[iCell]->GetYaxis()->SetRangeUser(gdMinimumToT, gdMaximumToT);
    }

    fAllToTSpectra->Reset();
    fAllToTSpectra->GetYaxis()->SetRangeUser(gdMinimumToT, gdMaximumToT);
  }


  void PlotHistograms(const char* cTag)
  {
    TCanvas* tCanvasTimeResidual = new TCanvas(Form("tCanvas%sTimeResidual_%d_%d_%d", cTag, fiModuleType, fiModuleIndex, fiCounterIndex),
                                               Form("%s time residual %d-%d-%d", cTag, fiModuleType, fiModuleIndex, fiCounterIndex),
                                               0, 0, 1400, 1400);
    tCanvasTimeResidual->Divide(6, 6);

    TCanvas* tCanvasTimeDifference = new TCanvas(Form("tCanvas%sTimeDifference_%d_%d_%d", cTag, fiModuleType, fiModuleIndex, fiCounterIndex),
                                                 Form("%s time difference %d-%d-%d", cTag, fiModuleType, fiModuleIndex, fiCounterIndex),
                                                 0, 0, 1400, 1400);
    tCanvasTimeDifference->Divide(6, 6);

    TCanvas* tCanvasTimeResidualWalk[2];

    for(Int_t iSide = 0; iSide < 2; iSide++)
    {
      tCanvasTimeResidualWalk[iSide] = new TCanvas(Form("tCanvas%sTimeResidualWalk%d_%d_%d_%d", cTag, iSide, fiModuleType, fiModuleIndex, fiCounterIndex),
                                                   Form("%s time residual walk %d_%d-%d-%d", cTag, iSide, fiModuleType, fiModuleIndex, fiCounterIndex),
                                                   0, 0, 1400, 1400);
      tCanvasTimeResidualWalk[iSide]->Divide(6, 6);
    }

    TCanvas* tCanvasTimeMeanCorrelation = new TCanvas(Form("tCanvas%sTimeMeanCorrelation_%d_%d_%d", cTag, fiModuleType, fiModuleIndex, fiCounterIndex),
                                                      Form("%s time mean correlation %d-%d-%d", cTag, fiModuleType, fiModuleIndex, fiCounterIndex),
                                                      0, 0, 1400, 1400);
    tCanvasTimeMeanCorrelation->Divide(6, 6);

    TCanvas* tCanvasTimeDifferenceCorrelation = new TCanvas(Form("tCanvas%sTimeDifferenceCorrelation_%d_%d_%d", cTag, fiModuleType, fiModuleIndex, fiCounterIndex),
                                                            Form("%s time difference correlation %d-%d-%d", cTag, fiModuleType, fiModuleIndex, fiCounterIndex),
                                                            0, 0, 1400, 1400);
    tCanvasTimeDifferenceCorrelation->Divide(6, 6);


    TH1F* tTimeResidual[giMaxNCells];
    TH1F* tTimeDifference[giMaxNCells];
    TH2F* tTimeResidualWalk[giMaxNCells][2];
    TProfile* tTimeResidualWalkProfile[giMaxNCells][2];
    TH1F* tTimeMeanCorrelation[giMaxNCells];
    TH1F* tTimeDifferenceCorrelation[giMaxNCells];
    TF1* tFitFunction(NULL);

    for(Int_t iCell = 0; iCell < fiNCells; iCell++)
    {
      tCanvasTimeResidual->cd(iCell + 1);
      tTimeResidual[iCell] = (TH1F*)fTimeResidual[iCell]->Clone(TString(fTimeResidual[iCell]->GetName()).Insert(1, cTag));
      tTimeResidual[iCell]->SetDirectory(gROOT);
      gPad->SetLogy();
      tTimeResidual[iCell]->Draw();
      tTimeResidual[iCell]->GetXaxis()->SetRangeUser(-2., 5.);

      tCanvasTimeDifference->cd(iCell + 1);
      tTimeDifference[iCell] = (TH1F*)fTimeDifference[iCell]->Clone(TString(fTimeDifference[iCell]->GetName()).Insert(1, cTag));
      tTimeDifference[iCell]->SetDirectory(gROOT);
      gPad->SetLogy();
      tTimeDifference[iCell]->Draw();
      if(4 == fiModuleType)
      {
        tTimeDifference[iCell]->GetXaxis()->SetRangeUser(-1., 1.);
      }
      else
      {
        tTimeDifference[iCell]->GetXaxis()->SetRangeUser(-3., 3.);
      }
      if(fTimeDifferenceFitFunction[iCell])
      {
        if(4 == fiModuleType)
        {
          fTimeDifferenceFitFunction[iCell]->SetRange(-0.25, 0.25); // TODO: avoid hard coding
        }
        else
        {
          fTimeDifferenceFitFunction[iCell]->SetRange(-1., 1.); // TODO: avoid hard coding
        }

        fTimeDifferenceFitFunction[iCell]->Draw("SAME");
      }

      for(Int_t iSide = 0; iSide < 2; iSide++)
      {
        tCanvasTimeResidualWalk[iSide]->cd(iCell + 1);
        tTimeResidualWalk[iCell][iSide] = (TH2F*)fTimeResidualWalk[iCell][iSide]->Clone(TString(fTimeResidualWalk[iCell][iSide]->GetName()).Insert(1, cTag));
        tTimeResidualWalk[iCell][iSide]->SetDirectory(gROOT);
        gPad->SetLogz();
        tTimeResidualWalk[iCell][iSide]->Draw("colz");
        if(gbToTScalingBeforeCalibration)
        {
          tTimeResidualWalk[iCell][iSide]->GetXaxis()->SetRangeUser(0., 5.);
        }
        else
        {
          tTimeResidualWalk[iCell][iSide]->GetXaxis()->SetRangeUser(10., 20.);
        }
        tTimeResidualWalkProfile[iCell][iSide] = tTimeResidualWalk[iCell][iSide]->ProfileX();
        tTimeResidualWalkProfile[iCell][iSide]->SetDirectory(gROOT);
        tTimeResidualWalkProfile[iCell][iSide]->Draw("SAME");
      }

      tCanvasTimeMeanCorrelation->cd(iCell + 1);
      tTimeMeanCorrelation[iCell] = (TH1F*)fTimeMeanCorrelation[iCell]->Clone(TString(fTimeMeanCorrelation[iCell]->GetName()).Insert(1, cTag));
      tTimeMeanCorrelation[iCell]->SetDirectory(gROOT);
      gPad->SetLogy();
      tTimeMeanCorrelation[iCell]->Draw();
      tTimeMeanCorrelation[iCell]->GetXaxis()->SetRangeUser(-3., 3.);

      tCanvasTimeDifferenceCorrelation->cd(iCell + 1);
      tTimeDifferenceCorrelation[iCell] = (TH1F*)fTimeDifferenceCorrelation[iCell]->Clone(TString(fTimeDifferenceCorrelation[iCell]->GetName()).Insert(1, cTag));
      tTimeDifferenceCorrelation[iCell]->SetDirectory(gROOT);
      gPad->SetLogy();
      tTimeDifferenceCorrelation[iCell]->Draw();
    }


    tCanvasTimeResidual->SaveAs(Form("plots/pl_%s_time_residual_%d-%d-%d.pdf", cTag, fiModuleType, fiModuleIndex, fiCounterIndex));
    tCanvasTimeDifference->SaveAs(Form("plots/pl_%s_time_difference_%d-%d-%d.pdf", cTag, fiModuleType, fiModuleIndex, fiCounterIndex));

    for(Int_t iSide = 0; iSide < 2; iSide++)
    {
      tCanvasTimeResidualWalk[iSide]->SaveAs(Form("plots/pl_%s_time_residual_walk_%d-%d-%d-%d.pdf", cTag, iSide, fiModuleType, fiModuleIndex, fiCounterIndex));
    }

    tCanvasTimeMeanCorrelation->SaveAs(Form("plots/pl_%s_time_mean_correlation_%d-%d-%d.pdf", cTag, fiModuleType, fiModuleIndex, fiCounterIndex));
    tCanvasTimeDifferenceCorrelation->SaveAs(Form("plots/pl_%s_time_difference_correlation_%d-%d-%d.pdf", cTag, fiModuleType, fiModuleIndex, fiCounterIndex));


    if(gROOT->IsBatch() && !gbPreserveHistograms)
    {
      delete tCanvasTimeResidual;
      delete tCanvasTimeDifference;

      for(Int_t iSide = 0; iSide < 2; iSide++)
      {
        delete tCanvasTimeResidualWalk[iSide];
      }

      delete tCanvasTimeMeanCorrelation;
      delete tCanvasTimeDifferenceCorrelation;


      for(Int_t iCell = 0; iCell < fiNCells; iCell++)
      {
        delete tTimeResidual[iCell];
        delete tTimeDifference[iCell];

        for(Int_t iSide = 0; iSide < 2; iSide++)
        {
          delete tTimeResidualWalk[iCell][iSide];
          delete tTimeResidualWalkProfile[iCell][iSide];
        }

        delete tTimeMeanCorrelation[iCell];
        delete tTimeDifferenceCorrelation[iCell];
      }
    }

  }


  void PlotOverview(const char* cTag)
  {
    TCanvas* tCanvasCounterOverview = new TCanvas(Form("tCanvas%sCounterOverview_%d_%d_%d", cTag, fiModuleType, fiModuleIndex, fiCounterIndex),
                                                  Form("%s counter overview %d-%d-%d", cTag, fiModuleType, fiModuleIndex, fiCounterIndex),
                                                  0, 0, 1400, 700);

    tCanvasCounterOverview->Divide(4, 2);


    TH2F* tAllTimeResidual(NULL);
    TH2F* tAllTimeMeanCorrelation(NULL);
    TH2F* tAllTimeDifferenceCorrelation(NULL);
    TH2F* tAllTimeDifference(NULL);
    TH2F* tAllLocalPosition(NULL);
    TH2F* tAllLocalPositionMeanVelocity(NULL);
    TProfile* tAllTimeResidualProfile(NULL);
    TProfile* tAllTimeMeanCorrelationProfile(NULL);
    TProfile* tAllTimeDifferenceCorrelationProfile(NULL);
    TProfile* tAllTimeDifferenceProfile(NULL);
    TProfile* tAllLocalPositionProfile(NULL);
    TProfile* tAllLocalPositionMeanVelocityProfile(NULL);
    TProfile* tSignalPropagationVelocity(NULL);
    TLine* tMeanVelocityLine(NULL);

    tCanvasCounterOverview->cd(1);
    tAllTimeResidual = (TH2F*)fAllTimeResidual->Clone(TString(fAllTimeResidual->GetName()).Insert(1, cTag));
    tAllTimeResidual->SetDirectory(gROOT);
    gPad->SetLogz();
    tAllTimeResidual->Draw("colz");
    tAllTimeResidualProfile = tAllTimeResidual->ProfileX();
    tAllTimeResidualProfile->SetDirectory(gROOT);
    tAllTimeResidualProfile->Draw("SAME");

    tCanvasCounterOverview->cd(2);
    tAllTimeMeanCorrelation = (TH2F*)fAllTimeMeanCorrelation->Clone(TString(fAllTimeMeanCorrelation->GetName()).Insert(1, cTag));
    tAllTimeMeanCorrelation->SetDirectory(gROOT);
    gPad->SetLogz();
    tAllTimeMeanCorrelation->Draw("colz");
    tAllTimeMeanCorrelationProfile = tAllTimeMeanCorrelation->ProfileX();
    tAllTimeMeanCorrelationProfile->SetDirectory(gROOT);
    tAllTimeMeanCorrelationProfile->Draw("SAME");

    tCanvasCounterOverview->cd(3);
    tAllTimeDifferenceCorrelation = (TH2F*)fAllTimeDifferenceCorrelation->Clone(TString(fAllTimeDifferenceCorrelation->GetName()).Insert(1, cTag));
    tAllTimeDifferenceCorrelation->SetDirectory(gROOT);
    gPad->SetLogz();
    tAllTimeDifferenceCorrelation->Draw("colz");
    tAllTimeDifferenceCorrelationProfile = tAllTimeDifferenceCorrelation->ProfileX();
    tAllTimeDifferenceCorrelationProfile->SetDirectory(gROOT);
    tAllTimeDifferenceCorrelationProfile->Draw("SAME");

    tCanvasCounterOverview->cd(5);
    tAllTimeDifference = (TH2F*)fAllTimeDifference->Clone(TString(fAllTimeDifference->GetName()).Insert(1, cTag));
    tAllTimeDifference->SetDirectory(gROOT);
    gPad->SetLogz();
    tAllTimeDifference->Draw("colz");
    tAllTimeDifferenceProfile = tAllTimeDifference->ProfileX();
    tAllTimeDifferenceProfile->SetDirectory(gROOT);
    tAllTimeDifferenceProfile->Draw("SAME");

    tCanvasCounterOverview->cd(6);
    tAllLocalPosition = (TH2F*)fAllLocalPosition->Clone(TString(fAllLocalPosition->GetName()).Insert(1, cTag));
    tAllLocalPosition->SetDirectory(gROOT);
    gPad->SetLogz();
    tAllLocalPosition->Draw("colz");
    tAllLocalPositionProfile = tAllLocalPosition->ProfileX();
    tAllLocalPositionProfile->SetDirectory(gROOT);
    tAllLocalPositionProfile->Draw("SAME");

    tCanvasCounterOverview->cd(7);
    tAllLocalPositionMeanVelocity = (TH2F*)fAllLocalPositionMeanVelocity->Clone(TString(fAllLocalPositionMeanVelocity->GetName()).Insert(1, cTag));
    tAllLocalPositionMeanVelocity->SetDirectory(gROOT);
    gPad->SetLogz();
    tAllLocalPositionMeanVelocity->Draw("colz");
    tAllLocalPositionMeanVelocityProfile = tAllLocalPositionMeanVelocity->ProfileX();
    tAllLocalPositionMeanVelocityProfile->SetDirectory(gROOT);
    tAllLocalPositionMeanVelocityProfile->Draw("SAME");

    tCanvasCounterOverview->cd(8);
    if(fSignalPropagationVelocity->GetEntries())
    {
      tSignalPropagationVelocity = (TProfile*)fSignalPropagationVelocity->Clone(TString(fSignalPropagationVelocity->GetName()).Insert(1, cTag));
      tSignalPropagationVelocity->SetDirectory(gROOT);
      tSignalPropagationVelocity->SetMarkerStyle(34);
      tSignalPropagationVelocity->SetMarkerSize(1);
      tSignalPropagationVelocity->GetYaxis()->SetRangeUser(TMath::Nint(tSignalPropagationVelocity->GetMean(2)) - 1., TMath::Nint(tSignalPropagationVelocity->GetMean(2)) + 1.);
      tSignalPropagationVelocity->Draw("P");
  /*
      gPad->Update();
      dynamic_cast<TPaveStats*>(tSignalPropagationVelocity->FindObject("stats"))->SetOptStat(1110);
      gPad->Modified();
  */
      tMeanVelocityLine = new TLine(0., tSignalPropagationVelocity->GetMean(2), static_cast<Double_t>(fiNCells), tSignalPropagationVelocity->GetMean(2));
      tMeanVelocityLine->SetLineColor(kBlue);
      tMeanVelocityLine->SetLineWidth(2);
      tMeanVelocityLine->SetLineStyle(2);
      tMeanVelocityLine->Draw("SAME");
    }


    tCanvasCounterOverview->SaveAs(Form("plots/pl_%s_counter_overview_%d-%d-%d.pdf", cTag, fiModuleType, fiModuleIndex, fiCounterIndex));


    if(gROOT->IsBatch() && !gbPreserveHistograms)
    {
      delete tCanvasCounterOverview;

      delete tAllTimeResidual;
      delete tAllTimeMeanCorrelation;
      delete tAllTimeDifferenceCorrelation;
      delete tAllTimeDifference;
      delete tAllLocalPosition;
      delete tAllLocalPositionMeanVelocity;
      delete tAllTimeResidualProfile;
      delete tAllTimeMeanCorrelationProfile;
      delete tAllTimeDifferenceCorrelationProfile;
      delete tAllTimeDifferenceProfile;
      delete tAllLocalPositionProfile;
      delete tAllLocalPositionMeanVelocityProfile;
      delete tSignalPropagationVelocity;
      delete tMeanVelocityLine;
    }
  }


  void PlotToTSpectra(const char* cTag, Bool_t bFinalAlignment)
  {
    TCanvas* tCanvasToTSpectra = new TCanvas(Form("tCanvas%sToTSpectra%d_%d_%d", cTag, fiModuleType, fiModuleIndex, fiCounterIndex),
                                             Form("%s ToT spectra %d-%d-%d", cTag, fiModuleType, fiModuleIndex, fiCounterIndex),
                                             0, 0, 1050, 350);

    tCanvasToTSpectra->Divide(3, 1);

    TCanvas* tCanvasToTCellSpectra[2];

    for(Int_t iSide = 0; iSide < 2; iSide++)
    {
      tCanvasToTCellSpectra[iSide] = new TCanvas(Form("tCanvas%sToTCellSpectra%d_%d_%d_%d", cTag, iSide, fiModuleType, fiModuleIndex, fiCounterIndex),
                                                 Form("%s ToT cell spectra %d_%d-%d-%d", cTag, iSide, fiModuleType, fiModuleIndex, fiCounterIndex),
                                                 0, 0, 1400, 1400);
      tCanvasToTCellSpectra[iSide]->Divide(6, 6);
    }

    TCanvas* tCanvasToTCorrelation = new TCanvas(Form("tCanvas%sToTCorrelation%d_%d_%d", cTag, fiModuleType, fiModuleIndex, fiCounterIndex),
                                                 Form("%s ToT cell correlation %d-%d-%d", cTag, fiModuleType, fiModuleIndex, fiCounterIndex),
                                                 0, 0, 1400, 1400);
    tCanvasToTCorrelation->Divide(6, 6);


    TH2F* tToTSpectra[2];
    TProfile* tToTSpectraProfile[2];
    TH1D* tToTCellSpectra[giMaxNCells][2];
    TH2F* tToTCorrelation[giMaxNCells];
    TProfile* tToTCorrelationProfile[giMaxNCells];
    TLine* tToTCorrelationLine[giMaxNCells];
    TH2F* tAllToTSpectra;
    TProfile* tAllToTSpectraProfile;

    for(Int_t iSide = 0; iSide < 2; iSide++)
    {
      tCanvasToTSpectra->cd(iSide + 1);
      tToTSpectra[iSide] = (TH2F*)fToTSpectra[iSide]->Clone(TString(fToTSpectra[iSide]->GetName()).Insert(1, cTag));
      tToTSpectra[iSide]->SetDirectory(gROOT);
      gPad->SetLogz();
      tToTSpectra[iSide]->Draw("colz");
      tToTSpectraProfile[iSide] = tToTSpectra[iSide]->ProfileX("_pfx", 1, -1, "s");
      tToTSpectraProfile[iSide]->SetDirectory(gROOT);
      tToTSpectraProfile[iSide]->Draw("SAME");
    }
    tCanvasToTSpectra->cd(3);
    tAllToTSpectra = (TH2F*)fAllToTSpectra->Clone(TString(fAllToTSpectra->GetName()).Insert(1, cTag));
    tAllToTSpectra->SetDirectory(gROOT);
    gPad->SetLogz();
    tAllToTSpectra->Draw("colz");
    tAllToTSpectraProfile = tAllToTSpectra->ProfileX("_pfx", 1, -1, "s");
    tAllToTSpectraProfile->SetDirectory(gROOT);
    tAllToTSpectraProfile->Draw("SAME");

    for(Int_t iCell = 0; iCell < fiNCells; iCell++)
    {
      for(Int_t iSide = 0; iSide < 2; iSide++)
      {
        tCanvasToTCellSpectra[iSide]->cd(iCell + 1);
        tToTCellSpectra[iCell][iSide] = fToTSpectra[iSide]->ProjectionY(Form("t%sToTSpectra_%d_%d_%d_%d_%d", cTag, fiModuleType, fiModuleIndex, fiCounterIndex, iCell, iSide), iCell + 1, iCell + 1);
        tToTCellSpectra[iCell][iSide]->SetDirectory(gROOT);
        gPad->SetLogy();
        tToTCellSpectra[iCell][iSide]->Draw();
      }

      tCanvasToTCorrelation->cd(iCell + 1);
      tToTCorrelation[iCell] = (TH2F*)fToTCorrelation[iCell]->Clone(TString(fToTCorrelation[iCell]->GetName()).Insert(1, cTag));
      tToTCorrelation[iCell]->SetDirectory(gROOT);
      gPad->SetLogz();
      tToTCorrelation[iCell]->Draw("colz");
      if(bFinalAlignment)
      {
        tToTCorrelation[iCell]->GetXaxis()->SetRangeUser(0., 5.);
        tToTCorrelation[iCell]->GetYaxis()->SetRangeUser(0., 5.);
      }
      else
      {
        tToTCorrelation[iCell]->GetXaxis()->SetRangeUser(10., 20.);
        tToTCorrelation[iCell]->GetYaxis()->SetRangeUser(10., 20.);
      }

      tToTCorrelationProfile[iCell] = tToTCorrelation[iCell]->ProfileX("_pfx", 1, -1, "s");
      tToTCorrelationProfile[iCell]->SetDirectory(gROOT);
      tToTCorrelationProfile[iCell]->SetMarkerStyle(22);
      tToTCorrelationProfile[iCell]->SetMarkerSize(0.5);
      tToTCorrelationProfile[iCell]->Draw("HIST P SAME");

      tToTCorrelationLine[iCell] = new TLine(gPad->GetLeftMargin(), gPad->GetBottomMargin(), 1. - gPad->GetRightMargin(), 1. - gPad->GetTopMargin());
      tToTCorrelationLine[iCell]->SetNDC(kTRUE);
      tToTCorrelationLine[iCell]->SetLineColor(kMagenta);
      tToTCorrelationLine[iCell]->SetLineWidth(1);
      tToTCorrelationLine[iCell]->SetLineStyle(1);
      tToTCorrelationLine[iCell]->Draw("SAME");

      gPad->Update();
      TPaveStats* tStats = dynamic_cast<TPaveStats*>(tToTCorrelation[iCell]->FindObject("stats"));
      gPad->GetListOfPrimitives()->Remove(tStats);
      gPad->Modified();
      gPad->Update();
      tStats->Draw();
    }

    tCanvasToTSpectra->SaveAs(Form("plots/pl_%s_tot_spectra_%d-%d-%d.pdf", cTag, fiModuleType, fiModuleIndex, fiCounterIndex));

    for(Int_t iSide = 0; iSide < 2; iSide++)
    {
      tCanvasToTCellSpectra[iSide]->SaveAs(Form("plots/pl_%s_tot_cell_spectra_%d-%d-%d-%d.pdf", cTag, iSide, fiModuleType, fiModuleIndex, fiCounterIndex));
    }

    tCanvasToTCorrelation->SaveAs(Form("plots/pl_%s_tot_correlation_%d-%d-%d.pdf", cTag, fiModuleType, fiModuleIndex, fiCounterIndex));


    if(gROOT->IsBatch() && !gbPreserveHistograms)
    {
      delete tCanvasToTSpectra;

      for(Int_t iSide = 0; iSide < 2; iSide++)
      {
        delete tCanvasToTCellSpectra[iSide];
      }

      delete tCanvasToTCorrelation;

      for(Int_t iSide = 0; iSide < 2; iSide++)
      {
        delete tToTSpectra[iSide];
        delete tToTSpectraProfile[iSide];
      }

      for(Int_t iCell = 0; iCell < fiNCells; iCell++)
      {
        for(Int_t iSide = 0; iSide < 2; iSide++)
        {
          delete tToTCellSpectra[iCell][iSide];
        }

        delete tToTCorrelation[iCell];
        delete tToTCorrelationProfile[iCell];
        delete tToTCorrelationLine[iCell];
      }

      delete tAllToTSpectra;
      delete tAllToTSpectraProfile;
    }
  }


  void ZoomHistograms()
  {
    for(Int_t iCell = 0; iCell < fiNCells; iCell++)
    {
      fTimeResidual[iCell]->GetXaxis()->SetRangeUser(gdMeanTimeResidual - 1., gdMeanTimeResidual + 1.);
      fTimeDifference[iCell]->GetXaxis()->SetRangeUser(-3., 3.);

      for(Int_t iSide = 0; iSide < 2; iSide++)
      {
        fTimeResidualWalk[iCell][iSide]->GetYaxis()->SetRangeUser(gdMeanTimeResidual - 1., gdMeanTimeResidual + 1.);
      }

// TODO
//      fTimeMeanCorrelation[iCell]->GetXaxis()->SetRangeUser(?, ?);
//      fTimeDifferenceCorrelation[iCell]->GetXaxis()->SetRangeUser(?, ?);
    }
  }


  void ZoomAllHistograms()
  {
    fAllTimeResidual->GetYaxis()->SetRangeUser(-2., 5.);
    fAllTimeMeanCorrelation->GetYaxis()->SetRangeUser(-3., 3.);
    fAllTimeDifferenceCorrelation->GetYaxis()->SetRangeUser(-5., 5.);
    if(4 == fiModuleType)
    {
      fAllTimeDifference->GetYaxis()->SetRangeUser(-0.5, 0.5);
    }
    else
    {
      fAllTimeDifference->GetYaxis()->SetRangeUser(-1.5, 1.5);
    }
    fAllLocalPosition->GetYaxis()->SetRangeUser(-TMath::Ceil(fdCellAlongWidth), TMath::Ceil(fdCellAlongWidth));
    fAllLocalPositionMeanVelocity->GetYaxis()->SetRangeUser(-TMath::Ceil(fdCellAlongWidth), TMath::Ceil(fdCellAlongWidth));
  }


  void ZoomToTSpectra()
  {
    for(Int_t iSide = 0; iSide < 2; iSide++)
    {
      fToTSpectra[iSide]->GetYaxis()->SetRangeUser(gdMinimumToT, 10.);
    }

    for(Int_t iCell = 0; iCell < fiNCells; iCell++)
    {
      fToTCorrelation[iCell]->GetXaxis()->SetRangeUser(gdMinimumToT, 10.);
      fToTCorrelation[iCell]->GetYaxis()->SetRangeUser(gdMinimumToT, 10.);
    }

    fAllToTSpectra->GetYaxis()->SetRangeUser(gdMinimumToT, 10.);
  }


  void FillToTSpectra()
  {
    for(Int_t iCell = 0; iCell < fiNCells; iCell++)
    {
      if(!fbDeadCell[iCell])
      {
        if(fbIsBRef)
        {
          if(gbBRefPadReadout)
          {
            if(fbDigiAvailable[iCell][0])
            {
              if(0. < fDigi[iCell][0]->GetTot())
              {
                fToTSpectra[0]->Fill(iCell, fDigi[iCell][0]->GetTot());
                fAllToTSpectra->Fill(2*iCell + 0, fDigi[iCell][0]->GetTot());
              }
            }
          }
          else
          {
            if(fbDigiAvailable[iCell][0] && fbDigiAvailable[iCell][1])
            {
              if(0. < fDigi[iCell][0]->GetTot() && 0. < fDigi[iCell][1]->GetTot())
              {
                for(Int_t iSide = 0; iSide < 2; iSide++)
                {
                  fToTSpectra[iSide]->Fill(iCell, fDigi[iCell][iSide]->GetTot());
                  fAllToTSpectra->Fill(2*iCell + iSide, fDigi[iCell][iSide]->GetTot());
                }

                fToTCorrelation[iCell]->Fill(fDigi[iCell][0]->GetTot(), fDigi[iCell][1]->GetTot());
              }
            }
          }
        }
        else
        {
          if(fbDigiAvailable[iCell][0] && fbDigiAvailable[iCell][1])
          {
            if(0. < fDigi[iCell][0]->GetTot() && 0. < fDigi[iCell][1]->GetTot())
            {
              for(Int_t iSide = 0; iSide < 2; iSide++)
              {
                fToTSpectra[iSide]->Fill(iCell, fDigi[iCell][iSide]->GetTot());
                fAllToTSpectra->Fill(2*iCell + iSide, fDigi[iCell][iSide]->GetTot());
              }

              fToTCorrelation[iCell]->Fill(fDigi[iCell][0]->GetTot(), fDigi[iCell][1]->GetTot());
            }
          }
        }
      }
    }
  }


  void ClearDigiArrays()
  {
    for(Int_t iCell = 0; iCell < fiNCells; iCell++)
    {
      for(Int_t iSide = 0; iSide < 2; iSide++)
      {
        fbDigiAvailable[iCell][iSide] = kFALSE;
        fDigi[iCell][iSide] = NULL;
      }
    }
  }


  void ApplyToTScaling()
  {
    for(Int_t iCell = 0; iCell < fiNCells; iCell++)
    {
      for(Int_t iSide = 0; iSide < 2; iSide++)
      {
        if(fbDigiAvailable[iCell][iSide])
        {
          Double_t dDigiToT = fDigi[iCell][iSide]->GetTot();

          if(0. < dDigiToT)
          {
            fDigi[iCell][iSide]->SetTot(fdToTScaling[iCell][iSide]*dDigiToT + fdToTOffset[iCell][iSide]);
          }
        }
      }
    }
  }


  void ComputeDigiWalk()
  {
    for(Int_t iCell = 0; iCell < fiNCells; iCell++)
    {
      for(Int_t iSide = 0; iSide < 2; iSide++)
      {
        fdDigiSlewing[iCell][iSide] = 0.;

        if(fbDigiAvailable[iCell][iSide])
        {
          Double_t dDigiToT = fDigi[iCell][iSide]->GetTot();

          // linear interpolation according to TH1::Interpolate
          // N.B.: corrects the BRef counter for ToT values that are out of range
          //       (necessary to get rid of additional time peaks in 5-2-0 spectra)
          if(fbIsBRef && dDigiToT < gdMinimumToT)
          {
            fdDigiSlewing[iCell][iSide] = fdSlewingCorrection[iCell][iSide][0];
          }
          else if(dDigiToT <= gdMinimumToT + 0.5*gdWalkBinWidth)
          {
            fdDigiSlewing[iCell][iSide] = fdSlewingCorrection[iCell][iSide][1];
          }
          else if(fbIsBRef && dDigiToT >= gdMaximumToT)
          {
            fdDigiSlewing[iCell][iSide] = fdSlewingCorrection[iCell][iSide][giNWalkBinsX + 1];
          }
          else if(dDigiToT >= gdMaximumToT - 0.5*gdWalkBinWidth)
          {
            fdDigiSlewing[iCell][iSide] = fdSlewingCorrection[iCell][iSide][giNWalkBinsX];
          }
          else
          {
            Int_t iWalkBin = 1 + static_cast<Int_t>((dDigiToT - gdMinimumToT)/gdWalkBinWidth );
            Double_t dWalkBinCenter = gdMinimumToT + (static_cast<Double_t>(iWalkBin) - 0.5)*gdWalkBinWidth;

            Double_t x0, x1, y0, y1;

            if(dDigiToT <= dWalkBinCenter)
            {
              y0 = fdSlewingCorrection[iCell][iSide][iWalkBin - 1];
              x0 = dWalkBinCenter - gdWalkBinWidth;
              y1 = fdSlewingCorrection[iCell][iSide][iWalkBin];
              x1 = dWalkBinCenter;
            }
            else
            {
              y0 = fdSlewingCorrection[iCell][iSide][iWalkBin];
              x0 = dWalkBinCenter;
              y1 = fdSlewingCorrection[iCell][iSide][iWalkBin + 1];
              x1 = dWalkBinCenter + gdWalkBinWidth;
            }

            fdDigiSlewing[iCell][iSide] = y0 + (dDigiToT - x0)*((y1 - y0)/(x1 - x0));
          }

          if(fbIsBRef)
          {
            fdDigiSlewing[iCell][iSide] *= -1.;

            if(!gbBRefPadReadout)
            {
              if(!gbParallelSlewingCorrections)
              {
                fdDigiSlewing[iCell][iSide] *= 2.;
              }
            }
          }
          else
          {
            if(!gbParallelSlewingCorrections)
            {
              fdDigiSlewing[iCell][iSide] *= 2.;
            }
          }
        }
      }
    }
  }


  void UpdateDelayOffsets(Bool_t bTimePeak, Bool_t bFitTimeDifference, Bool_t bLastFit)
  {
    if(bFitTimeDifference && bLastFit)
    {
      for(Int_t iCell = 0; iCell < fiNCells; iCell++)
      {
        gTimeDifferenceFitNH->SetRange(-3., 3.);
        gTimeDifferenceFitNH->SetParameter(0, fTimeDifference[iCell]->GetMean());
        gTimeDifferenceFitNH->SetParameter(1, 0.03);
        gTimeDifferenceFitNH->SetParameter(2, 1.5);
        gTimeDifferenceFitNH->SetParameter(3, fTimeDifference[iCell]->Integral("width"));
        gTimeDifferenceFitNH->SetParameter(4, -0.5);
        gTimeDifferenceFitNH->SetParameter(5, 0.);

        gTimeDifferenceFitNH->SetParLimits(1, 0.01, 0.05);
        gTimeDifferenceFitNH->SetParLimits(2, 1.3, 1.7);
//        gTimeDifferenceFitNH->SetParLimits(4, -1.5, 1.5);
//        gTimeDifferenceFitNH->SetParLimits(5, -3., 3.);

        if(4 == fiModuleType)
        {
          gTimeDifferenceFitNH->SetRange(-1., 1.);
          gTimeDifferenceFitNH->SetParameter(2, 0.28);
          gTimeDifferenceFitNH->SetParLimits(2, 0.2, 0.4);
//          gTimeDifferenceFitNH->SetParLimits(4, -3., 3.);

          fTimeDifference[iCell]->Fit("gTimeDifferenceFitNH", "NQ", "", -0.25, 0.25); // TODO: avoid hard coding
        }
        else
        {
          fTimeDifference[iCell]->Fit("gTimeDifferenceFitNH", "NQ", "", -1., 1.);
        }

        // TODO: identify and repeat bad fits
        fTimeDifferenceFitFunction[iCell] = new TF1(*gTimeDifferenceFitNH);
        fTimeDifferenceFitFunction[iCell]->SetName(Form("tTimeDifferenceFitFunction_%d_%d_%d_%d", fiModuleType, fiModuleIndex, fiCounterIndex, iCell));
        fTimeDifferenceFitFunction[iCell]->SetParameter(0, 0.);
        fTimeDifferenceFitFunction[iCell]->SetParameter(5, gTimeDifferenceFitNH->GetParameter(5) - gTimeDifferenceFitNH->GetParameter(0));

        if(0. != gTimeDifferenceFitNH->GetParameter(2))
        {
          fdSignalPropagationVelocity[iCell] = fdCellAlongWidth/gTimeDifferenceFitNH->GetParameter(2);

          if(!fbDeadCell[iCell])
          {
            fSignalPropagationVelocity->Fill(iCell, fdSignalPropagationVelocity[iCell]);
          }
        }

        cout<<Form("%d-%d-%d-%d: velocity %.3f cm/ns; location %.3f ns; scale %.3f ns; shape %.3f ns; norm %.3f; p-scale %.3f 1/ns^2; p-vertex %.3f ns",
                   fiModuleType, fiModuleIndex, fiCounterIndex, iCell,
                   fdSignalPropagationVelocity[iCell],
                   fTimeDifferenceFitFunction[iCell]->GetParameter(0),
                   fTimeDifferenceFitFunction[iCell]->GetParameter(1),
                   fTimeDifferenceFitFunction[iCell]->GetParameter(2),
                   fTimeDifferenceFitFunction[iCell]->GetParameter(3),
                   fTimeDifferenceFitFunction[iCell]->GetParameter(4),
                   fTimeDifferenceFitFunction[iCell]->GetParameter(5))<<endl;
      }

      fdMeanSignalPropagationVelocity = fSignalPropagationVelocity->GetMean(2);

      cout<<Form("%d-%d-%d: mean velocity %.3f cm/ns", fiModuleType, fiModuleIndex, fiCounterIndex, fdMeanSignalPropagationVelocity)<<endl;

      if(gbUseInputSignalPropagationVelocity)
      {
        fdMeanSignalPropagationVelocity = fdInputSignalPropagationVelocity;
      }

      cout<<Form("%d-%d-%d: appl velocity %.3f cm/ns", fiModuleType, fiModuleIndex, fiCounterIndex, fdMeanSignalPropagationVelocity)<<endl;
    }


    for(Int_t iCell = 0; iCell < fiNCells; iCell++)
    {
      Double_t dTimeResidual = fTimeResidual[iCell]->GetMean();
      Double_t dTimeDifference = fTimeDifference[iCell]->GetMean();

      if(bTimePeak)
      {
        dTimeResidual = fTimeResidual[iCell]->GetBinCenter(fTimeResidual[iCell]->GetMaximumBin());

        if(7 != fiModuleType) // TODO: avoid counter-specific exception handling
        {
          dTimeDifference = fTimeDifference[iCell]->GetBinCenter(fTimeDifference[iCell]->GetMaximumBin());
        }
      }

      if(bFitTimeDifference)
      {
        if(bLastFit)
        {
          dTimeResidual -= 0.5*fdCellAlongWidth/fdMeanSignalPropagationVelocity;

          dTimeDifference = fTimeDifferenceFitFunction[iCell]->GetParameter(0);
        }
        else
        {
          gTimeDifferenceFitNH->SetRange(-3., 3.);
          gTimeDifferenceFitNH->SetParameter(0, fTimeDifference[iCell]->GetMean());
          gTimeDifferenceFitNH->SetParameter(1, 0.03);
          gTimeDifferenceFitNH->SetParameter(2, 1.5);
          gTimeDifferenceFitNH->SetParameter(3, fTimeDifference[iCell]->Integral("width"));
          gTimeDifferenceFitNH->SetParameter(4, -0.5);
          gTimeDifferenceFitNH->SetParameter(5, 0.);

          gTimeDifferenceFitNH->SetParLimits(1, 0.01, 0.05);
          gTimeDifferenceFitNH->SetParLimits(2, 1.3, 1.7);
//          gTimeDifferenceFitNH->SetParLimits(4, -1.5, 1.5);
//          gTimeDifferenceFitNH->SetParLimits(5, -3., 3.);

          if(4 == fiModuleType)
          {
            gTimeDifferenceFitNH->SetRange(-1., 1.);
            gTimeDifferenceFitNH->SetParameter(2, 0.28);
            gTimeDifferenceFitNH->SetParLimits(2, 0.2, 0.4);
//            gTimeDifferenceFitNH->SetParLimits(4, -3., 3.);
          }

          fTimeDifference[iCell]->Fit("gTimeDifferenceFitNH", "RNQ");
          // TODO: identify and repeat bad fits

          dTimeDifference = gTimeDifferenceFitNH->GetParameter(0);
        }
      }

      if(fbIsBRef)
      {
        dTimeResidual *= -1.;

        if(gbBRefPadReadout)
        {
          dTimeDifference = 0.;
        }
      }

      for(Int_t iSide = 0; iSide < 2; iSide++)
      {
        fdSignalPropagationDelay[iCell][iSide] += dTimeResidual + (0 == iSide ? 1. : -1.)*dTimeDifference;
      }
    }
  }


  void UpdateSlewingOffsets(Bool_t bShiftPeak, Int_t iSideToUpdate)
  {
    if(bShiftPeak)
    {
      for(Int_t iCell = 0; iCell < fiNCells; iCell++)
      {
        for(Int_t iSide = 0; iSide < 2; iSide++)
        {
          if(iSideToUpdate != iSide)
          {
            if(!gbParallelSlewingCorrections || (fbIsBRef && gbBRefPadReadout))
            {
              continue;
            }
          }

          for(Int_t iWalkBin = 0; iWalkBin < giNWalkBinsX + 2; iWalkBin++)
          {
            TH1D* tSlewingProjection = fTimeResidualWalk[iCell][iSide]->ProjectionY("_py", iWalkBin, iWalkBin);
            Int_t iMaximumBin = tSlewingProjection->GetMaximumBin();

            if( 0. != tSlewingProjection->GetBinContent(iMaximumBin))
            {
              fdSlewingCorrection[iCell][iSide][iWalkBin] += tSlewingProjection->GetBinCenter(iMaximumBin) - gdMeanTimeResidual;
            }

            delete tSlewingProjection;
          }
        }
      }
    }
    else
    {
      for(Int_t iCell = 0; iCell < fiNCells; iCell++)
      {
        for(Int_t iSide = 0; iSide < 2; iSide++)
        {
          if(iSideToUpdate != iSide)
          {
            if(!gbParallelSlewingCorrections || (fbIsBRef && gbBRefPadReadout))
            {
              continue;
            }
          }

          TProfile* tSlewingProfile = fTimeResidualWalk[iCell][iSide]->ProfileX();

          for(Int_t iWalkBin = 0; iWalkBin < giNWalkBinsX + 2; iWalkBin++)
          {
            // transition from a minimum entries criterion to a standard error of the mean (SEM) criterion
            if( 0. != tSlewingProfile->GetBinError(iWalkBin) )// && 0.01 > tSlewingProfile->GetBinError(iWalkBin) )
            {
              fdSlewingCorrection[iCell][iSide][iWalkBin] += tSlewingProfile->GetBinContent(iWalkBin) - gdMeanTimeResidual;
            }
          }

          delete tSlewingProfile;
        }
      }
    }
  }


  void UpdateCorrelationOffsets()
  {
    Double_t dTimeMeanCorrelation[giMaxNCells] = {0.};
    Double_t dTimeDifferenceCorrelation[giMaxNCells] = {0.};

    for(Int_t iCell = fiReferenceCell + 1; iCell < fiNCells; iCell++)
    {
      dTimeMeanCorrelation[iCell] = fTimeMeanCorrelation[iCell]->GetMean();
      dTimeDifferenceCorrelation[iCell] = fTimeDifferenceCorrelation[iCell]->GetMean();

      Double_t dTimeMean(0.);
      Double_t dTimeDifference(0.);

      for(Int_t iPreviousCell = fiReferenceCell + 1; iPreviousCell <= iCell; iPreviousCell++)
      {
        dTimeMean += dTimeMeanCorrelation[iPreviousCell];
        dTimeDifference += dTimeDifferenceCorrelation[iPreviousCell];
      }

      for(Int_t iSide = 0; iSide < 2; iSide++)
      {
        fdSignalPropagationDelay[iCell][iSide] += dTimeMean + (0 == iSide ? 1. : -1.)*dTimeDifference;
      }
    }

    for(Int_t iCell = fiReferenceCell - 1; iCell > - 1; iCell--)
    {
      dTimeMeanCorrelation[iCell] = fTimeMeanCorrelation[iCell]->GetMean();
      dTimeDifferenceCorrelation[iCell] = fTimeDifferenceCorrelation[iCell]->GetMean();

      Double_t dTimeMean(0.);
      Double_t dTimeDifference(0.);

      for(Int_t iPreviousCell = fiReferenceCell - 1; iPreviousCell >= iCell; iPreviousCell--)
      {
        dTimeMean += dTimeMeanCorrelation[iPreviousCell];
        dTimeDifference += dTimeDifferenceCorrelation[iPreviousCell];
      }

      for(Int_t iSide = 0; iSide < 2; iSide++)
      {
        fdSignalPropagationDelay[iCell][iSide] += dTimeMean + (0 == iSide ? 1. : -1.)*dTimeDifference;
      }
    }
  }


  void UpdateToTAlignment()
  {
    TProfile* tToTSpectraProfile[2];

    for(Int_t iSide = 0; iSide < 2; iSide++)
    {
      tToTSpectraProfile[iSide] = fToTSpectra[iSide]->ProfileX("_pfx", 1, -1, "s");
    }

    for(Int_t iCell = 0; iCell < fiNCells; iCell++)
    {
      if(!fbDeadCell[iCell])
      {
        for(Int_t iSide = 0; iSide < 2; iSide++)
        {
          if(0. != tToTSpectraProfile[iSide]->GetBinError(iCell + 1))
          {
            fdToTScaling[iCell][iSide] = gdToTStandardDeviation/tToTSpectraProfile[iSide]->GetBinError(iCell + 1);
            fdToTOffset[iCell][iSide]  = gdToTExpectationValue - fdToTScaling[iCell][iSide]*tToTSpectraProfile[iSide]->GetBinContent(iCell + 1);
          }
        }
      }
    }

    for(Int_t iSide = 0; iSide < 2; iSide++)
    {
      delete tToTSpectraProfile[iSide];
    }
  }


  void ClearCalibration()
  {
    for(Int_t iCell = 0; iCell < fiNCells; iCell++)
    {
      for(Int_t iSide = 0; iSide < 2; iSide++)
      {
        fdSignalPropagationDelay[iCell][iSide] = 0.;

        for(Int_t iWalkBin = 0; iWalkBin < giNWalkBinsX + 2; iWalkBin++)
        {
          fdSlewingCorrection[iCell][iSide][iWalkBin] = 0.;
        }
      }

      fdSignalPropagationVelocity[iCell] = gdDefaultSignalVelocity;
    }
  }


  void CalibrateDigi(CbmTofDigiExp* tDigi)
  {
    Int_t iDigiCell = static_cast<Int_t>(tDigi->GetChannel());
    Int_t iDigiSide = static_cast<Int_t>(tDigi->GetSide());
    Double_t dDigiTime = tDigi->GetTime();
    Double_t dDigiToT = tDigi->GetTot();

    Double_t dDigiSlewing(0.);

    // linear interpolation according to TH1::Interpolate
    // N.B.: corrects the BRef counter for ToT values that are out of range
    //       (necessary to get rid of additional time peaks in 5-2-0 spectra)
    if(fbIsBRef && dDigiToT < gdMinimumToT)
    {
      dDigiSlewing = fdSlewingCorrection[iDigiCell][iDigiSide][0];
    }
    else if(dDigiToT <= gdMinimumToT + 0.5*gdWalkBinWidth)
    {
      dDigiSlewing = fdSlewingCorrection[iDigiCell][iDigiSide][1];
    }
    else if(fbIsBRef && dDigiToT >= gdMaximumToT)
    {
      dDigiSlewing = fdSlewingCorrection[iDigiCell][iDigiSide][giNWalkBinsX + 1];
    }
    else if(dDigiToT >= gdMaximumToT - 0.5*gdWalkBinWidth)
    {
      dDigiSlewing = fdSlewingCorrection[iDigiCell][iDigiSide][giNWalkBinsX];
    }
    else
    {
      Int_t iWalkBin = 1 + static_cast<Int_t>((dDigiToT - gdMinimumToT)/gdWalkBinWidth );
      Double_t dWalkBinCenter = gdMinimumToT + (static_cast<Double_t>(iWalkBin) - 0.5)*gdWalkBinWidth;

      Double_t x0, x1, y0, y1;

      if( dDigiToT <= dWalkBinCenter )
      {
        y0 = fdSlewingCorrection[iDigiCell][iDigiSide][iWalkBin - 1];
        x0 = dWalkBinCenter - gdWalkBinWidth;
        y1 = fdSlewingCorrection[iDigiCell][iDigiSide][iWalkBin];
        x1 = dWalkBinCenter;
      }
      else
      {
        y0 = fdSlewingCorrection[iDigiCell][iDigiSide][iWalkBin];
        x0 = dWalkBinCenter;
        y1 = fdSlewingCorrection[iDigiCell][iDigiSide][iWalkBin + 1];
        x1 = dWalkBinCenter + gdWalkBinWidth;
      }

      dDigiSlewing = y0 + (dDigiToT - x0)*((y1 - y0)/(x1 - x0));
    }

    if(fbIsBRef)
    {
      dDigiSlewing *= -1.;

      if(!gbBRefPadReadout)
      {
        if(!gbParallelSlewingCorrections)
        {
          dDigiSlewing *= 2.;
        }
      }
    }
    else
    {
      if(!gbParallelSlewingCorrections)
      {
        dDigiSlewing *= 2.;
      }
    }

    tDigi->SetTime(dDigiTime - fdSignalPropagationDelay[iDigiCell][iDigiSide] - dDigiSlewing);


    if(0. < dDigiToT && !gbPreCalibration)
    {
      tDigi->SetTot(fdToTScaling[iDigiCell][iDigiSide]*dDigiToT + fdToTOffset[iDigiCell][iDigiSide]);
    }
  }


};


void CalibrateBRef(UInt_t uBRefAddress, Bool_t bCalibrateBRefCalibRef);
void CalibrateCounters();
void CalibrateToT();
void WriteCalibratedDigis();
void FillDigiArrays();
void CalibrateDigi(CbmTofDigiExp* tDigi, TCounter* tCounter);
void FillHistograms(TCounter* tCounter, Bool_t bBRefVelocityKnown, Bool_t bCounterVelocityKnown, Bool_t bFillAllHistograms);
Bool_t AddCounter(Int_t iModuleType, Int_t iModuleIndex, Int_t iCounterIndex);
Bool_t FindTofNode();
void ExpandNode(TGeoNode* tMotherNode);
Double_t gpb(Double_t* x, Double_t* p);
Double_t nh(Double_t* x, Double_t* p);


void calib_digis(Int_t iMRef, Int_t iBRef, const TString& tUnpackDir, const TString& tTofGeoVersion, Bool_t bNanoSeconds, Bool_t bPreCalibration)
{
  gStyle->SetOptStat("emr");
//  GeneralStyle();
  gStyle->SetPaintTextFormat(".2f");
  gErrorIgnoreLevel = kWarning;
  FairLogger::GetLogger()->SetLogScreenLevel("FATAL");
  FairLogger::GetLogger()->SetLogVerbosityLevel("LOW");

  if(gRandom)
  {
    delete gRandom;
  }
  gRandom = new TRandom3(4357);

  gcInputPath      = tUnpackDir.Data();
  gcInputFileName  = "data/unpack.out.root";
  gcOutputFileName = "data/digi.in.root";
//  gcCalibFileName  = "data/digi.calib.root"; TODO: implement functionality to store and retrieve calibration histograms from file
  gcGeoVersion = tTofGeoVersion.Data();

  gbNanoSeconds = bNanoSeconds;
  gbPreCalibration = bPreCalibration;

  Int_t iMRefCounter = iMRef%10;
  Int_t iMRefModule  = ((iMRef - iMRefCounter)/10)%10;
  Int_t iMRefType    = (iMRef - iMRefCounter - 10*iMRefModule)/100;

  Int_t iBRefCounter = iBRef%10;
  Int_t iBRefModule  = ((iBRef - iBRefCounter)/10)%10;
  Int_t iBRefType    = (iBRef - iBRefCounter - 10*iBRefModule)/100;

  giBRefModuleType = iBRefType;
  giBRefModuleIndex = iBRefModule;
  giBRefCounterIndex = iBRefCounter;

  if(5 == giBRefModuleType)
  {
    gbBRefPadReadout = kTRUE;
  }

  giBRefCalibRefModuleType = iMRefType;
  giBRefCalibRefModuleIndex = iMRefModule;
  giBRefCalibRefCounterIndex = iMRefCounter;

  guBRefAddress = static_cast<UInt_t>(CbmTofAddress::GetModFullId(CbmTofAddress::GetUniqueAddress(giBRefModuleIndex, giBRefCounterIndex, 0, 0, giBRefModuleType)));
  guBRefCalibRefAddress = static_cast<UInt_t>(CbmTofAddress::GetModFullId(CbmTofAddress::GetUniqueAddress(giBRefCalibRefModuleIndex, giBRefCalibRefCounterIndex, 0, 0, giBRefCalibRefModuleType)));


  if(gCounterInformation.find(std::make_tuple(giBRefModuleType, giBRefModuleIndex, giBRefCounterIndex)) == gCounterInformation.end())
  {
    cout<<Form("BRef counter %d-%d-%d is not in the calibration set!", giBRefModuleType, giBRefModuleIndex, giBRefCounterIndex)<<endl;
    return;
  }

  if(gCounterInformation.find(std::make_tuple(giBRefCalibRefModuleType, giBRefCalibRefModuleIndex, giBRefCalibRefCounterIndex)) == gCounterInformation.end())
  {
    cout<<Form("MRef counter %d-%d-%d is not in the calibration set!", giBRefCalibRefModuleType, giBRefCalibRefModuleIndex, giBRefCalibRefCounterIndex)<<endl;
    return;
  }


  if(!FindTofNode())
  {
    return;
  }
  else
  {
    ExpandNode(gGeoManager->GetCurrentNode());
    gGeoManager->CdTop();
  }


  for(auto const & Counter: gCounterInformation)
  {
    if(!AddCounter(std::get<0>(Counter.first), std::get<1>(Counter.first), std::get<2>(Counter.first)))
    {
      return;
    }
  }


  TFile::Open(TString(gcInputPath) + "/" + gcInputFileName, "read");
  if(!gFile)
  {
    cout<<Form("Could not open input file %s.", gcInputFileName)<<endl;
    return;
  }


  TTree* tTree = dynamic_cast<TTree*>(gFile->FindObjectAny("cbmsim"));
  if(!tTree)
  {
    cout<<"TTree 'cbmsim' could not be retrieved from input file."<<endl;
    return;
  }

  giNEvents = tTree->GetEntries();


  gInputDigiArray = new TClonesArray("CbmTofDigiExp");

  gInputDigiBranch = tTree->GetBranch("CbmTofDigi");
  gInputDigiBranch->SetAddress(&gInputDigiArray);

  gTimeDifferenceFitNH = new TF1("gTimeDifferenceFitNH", nh, -3., 3., 6);
  gTimeDifferenceFitNH->SetNpx(1000000);
  gTimeDifferenceFitNH->SetLineColor(kRed);

  TStopwatch timer;
  timer.Start();


  if(gbToTScalingBeforeCalibration)
  {
    CalibrateToT();
  }

  // calibrate MRef against main BRef
  CalibrateBRef(guBRefAddress, kTRUE);

  // calibrate BRefs against MRef (previously calibrated against the main BRef)
  for(auto const & uAddress: gBRefAddressSet)
  {
    CalibrateBRef(uAddress, kFALSE);
  }

  CalibrateCounters();

  if(!gbToTScalingBeforeCalibration)
  {
    CalibrateToT();
  }

  WriteCalibratedDigis();


  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();

  cout<<endl;
  cout<<" real time " << rtime << " s, CPU time " << ctime << " s" << endl;
  cout<<endl;
}


void CalibrateBRef(UInt_t uBRefAddress, Bool_t bCalibrateBRefCalibRef)
{
  cout<<Form("CalibrateBRef: initial delay correction...")<<endl;

  for(auto const & MapElement : gCalibCounterMap)
  {
    UInt_t uCounterAddress = MapElement.first;
    TCounter* tCounter = MapElement.second;

    tCounter->ResetHistograms();
    if(!bCalibrateBRefCalibRef)
    {
      tCounter->ResetAllHistograms();
    }

    if(uBRefAddress == uCounterAddress && !bCalibrateBRefCalibRef)
    {
      tCounter->ClearCalibration();
    }
  }


  for(Int_t iEvent = 0; iEvent < giNEvents; iEvent++)
  {
    if(0 == (iEvent + 1) % 100000)
    {
      cout<<Form("   event %4d", iEvent + 1)<<endl;
    }

    gInputDigiArray->Delete();
    gInputDigiBranch->GetEntry(iEvent);

    FillDigiArrays();

    for(auto const & MapElement : gCalibCounterMap)
    {
      UInt_t uCounterAddress = MapElement.first;
      TCounter* tCounter = MapElement.second;

      Bool_t bBRefVelocityKnown = kFALSE;
      Bool_t bCounterVelocityKnown = kFALSE;
      Bool_t bFillAllHistograms = (!bCalibrateBRefCalibRef ? kTRUE : kFALSE);

      if(uBRefAddress == uCounterAddress)
      {
        FillHistograms(tCounter, bBRefVelocityKnown, bCounterVelocityKnown, bFillAllHistograms);
      }
    }
  }


  for(auto const & MapElement : gCalibCounterMap)
  {
    UInt_t uCounterAddress = MapElement.first;
    TCounter* tCounter = MapElement.second;

    Bool_t bTimePeak = kTRUE;
    Bool_t bFitTimeDifference = kFALSE;
    Bool_t bLastFit = kFALSE;

    if(uBRefAddress == uCounterAddress)
    {
      tCounter->UpdateDelayOffsets(bTimePeak, bFitTimeDifference, bLastFit);
      if(!bCalibrateBRefCalibRef)
      {
        tCounter->PlotHistograms("OriginalBRefCalib");
        tCounter->PlotOverview("OriginalBRefCalib");
      }
    }
  }


  cout<<Form("CalibrateBRef: walk corrections...")<<endl;

  for( Int_t iCorrection = 0; iCorrection < giNBRefWalkCorrections + 1; iCorrection++)
  {
    for(auto const & MapElement : gCalibCounterMap)
    {
      TCounter* tCounter = MapElement.second;
      tCounter->ResetHistograms();
      if(0 < iCorrection)
      {
        tCounter->ZoomHistograms();
      }
      if(!bCalibrateBRefCalibRef && giNBRefWalkCorrections == iCorrection)
      {
        tCounter->ResetAllHistograms();
        tCounter->ZoomAllHistograms();
      }
    }


    for(Int_t iEvent = 0; iEvent < giNEvents; iEvent++)
    {
      if(0 == (iEvent + 1) % 100000)
      {
        cout<<Form("   event %4d", iEvent + 1)<<endl;
      }

      gInputDigiArray->Delete();
      gInputDigiBranch->GetEntry(iEvent);

      FillDigiArrays();

      for(auto const & MapElement : gCalibCounterMap)
      {
        UInt_t uCounterAddress = MapElement.first;
        TCounter* tCounter = MapElement.second;

        Bool_t bBRefVelocityKnown = kFALSE;
        Bool_t bCounterVelocityKnown = kFALSE;
        Bool_t bFillAllHistograms = ((!bCalibrateBRefCalibRef && giNBRefWalkCorrections == iCorrection) ? kTRUE : kFALSE);

        if(uBRefAddress == uCounterAddress || (bCalibrateBRefCalibRef && guBRefCalibRefAddress == uCounterAddress))
        {
          FillHistograms(tCounter, bBRefVelocityKnown, bCounterVelocityKnown, bFillAllHistograms);
        }
      }
    }


    for(auto const & MapElement : gCalibCounterMap)
    {
      UInt_t uCounterAddress = MapElement.first;
      TCounter* tCounter = MapElement.second;

      // TODO: make slewing corrections for BRef and BRefCalibRef counters optional

      if(gbBRefPadReadout)
      {
        if((0 == iCorrection%3 && uBRefAddress == uCounterAddress) ||
           (bCalibrateBRefCalibRef && (0 != iCorrection%3 && guBRefCalibRefAddress == uCounterAddress)))
        {
          Bool_t bShiftPeak = kFALSE;

          if(3 > iCorrection)
          {
            bShiftPeak = kTRUE;
          }

          // Do not update slewing offsets in the final correction step
          if(giNBRefWalkCorrections != iCorrection)
          {
            tCounter->UpdateSlewingOffsets(bShiftPeak, (uBRefAddress == uCounterAddress ? iCorrection%3 : iCorrection%3 - 1));
          }
        }
      }
      else
      {
        if(((0 == iCorrection%4 || 1 == iCorrection%4) && uBRefAddress == uCounterAddress) ||
           (bCalibrateBRefCalibRef && ((2 == iCorrection%4 || 3 == iCorrection%4) && guBRefCalibRefAddress == uCounterAddress)))
        {
          Bool_t bShiftPeak = kFALSE;

          if(4 > iCorrection)
          {
            bShiftPeak = kTRUE;
          }

          // Do not update slewing offsets in the final correction step
          if(giNBRefWalkCorrections != iCorrection)
          {
            tCounter->UpdateSlewingOffsets(bShiftPeak, (uBRefAddress == uCounterAddress ? iCorrection%4 : iCorrection%4 - 2));
          }
        }
      }

      if((uBRefAddress == uCounterAddress && !bCalibrateBRefCalibRef) || (bCalibrateBRefCalibRef && guBRefCalibRefAddress == uCounterAddress))
      {
        if(0 == iCorrection)
        {
//          tCounter->PlotHistograms("OriginalBRefWalk");
        }
        else if(giNBRefWalkCorrections == iCorrection)
        {
          tCounter->PlotHistograms("FinalBRefWalk");
          tCounter->PlotOverview("FinalBRefWalk");
        }
      }
    }

  }


  if(!gbBRefPadReadout)
  {
    cout<<Form("CalibrateBRef: initial time difference fit correction...")<<endl;

    for(auto const & MapElement : gCalibCounterMap)
    {
      UInt_t uCounterAddress = MapElement.first;
      TCounter* tCounter = MapElement.second;

      tCounter->ResetHistograms();
      tCounter->ZoomHistograms();
    }


    for(Int_t iEvent = 0; iEvent < giNEvents; iEvent++)
    {
      if(0 == (iEvent + 1) % 100000)
      {
        cout<<Form("   event %4d", iEvent + 1)<<endl;
      }

      gInputDigiArray->Delete();
      gInputDigiBranch->GetEntry(iEvent);

      FillDigiArrays();

      for(auto const & MapElement : gCalibCounterMap)
      {
        UInt_t uCounterAddress = MapElement.first;
        TCounter* tCounter = MapElement.second;

        Bool_t bBRefVelocityKnown = kFALSE;
        Bool_t bCounterVelocityKnown = kFALSE;
        Bool_t bFillAllHistograms = kFALSE;

        if(uBRefAddress == uCounterAddress)
        {
          FillHistograms(tCounter, bBRefVelocityKnown, bCounterVelocityKnown, bFillAllHistograms);
        }
      }
    }


    for(auto const & MapElement : gCalibCounterMap)
    {
      UInt_t uCounterAddress = MapElement.first;
      TCounter* tCounter = MapElement.second;

      Bool_t bTimePeak = kFALSE;
      Bool_t bFitTimeDifference = kTRUE;
      Bool_t bLastFit = kFALSE;

      if(uBRefAddress == uCounterAddress)
      {
        tCounter->UpdateDelayOffsets(bTimePeak, bFitTimeDifference, bLastFit);
      }
    }


    cout<<Form("CalibrateBRef: final time difference fit correction...")<<endl;

    for(auto const & MapElement : gCalibCounterMap)
    {
      UInt_t uCounterAddress = MapElement.first;
      TCounter* tCounter = MapElement.second;

      tCounter->ResetHistograms();
      tCounter->ZoomHistograms();
    }


    for(Int_t iEvent = 0; iEvent < giNEvents; iEvent++)
    {
      if(0 == (iEvent + 1) % 100000)
      {
        cout<<Form("   event %4d", iEvent + 1)<<endl;
      }

      gInputDigiArray->Delete();
      gInputDigiBranch->GetEntry(iEvent);

      FillDigiArrays();

      for(auto const & MapElement : gCalibCounterMap)
      {
        UInt_t uCounterAddress = MapElement.first;
        TCounter* tCounter = MapElement.second;

        Bool_t bBRefVelocityKnown = kFALSE;
        Bool_t bCounterVelocityKnown = kFALSE;
        Bool_t bFillAllHistograms = kFALSE;

        if(uBRefAddress == uCounterAddress)
        {
          FillHistograms(tCounter, bBRefVelocityKnown, bCounterVelocityKnown, bFillAllHistograms);
        }
      }
    }


    for(auto const & MapElement : gCalibCounterMap)
    {
      UInt_t uCounterAddress = MapElement.first;
      TCounter* tCounter = MapElement.second;

      Bool_t bTimePeak = kFALSE;
      Bool_t bFitTimeDifference = kTRUE;
      Bool_t bLastFit = kTRUE;

      if(uBRefAddress == uCounterAddress)
      {
        tCounter->UpdateDelayOffsets(bTimePeak, bFitTimeDifference, bLastFit);
      }
    }


    cout<<Form("CalibrateBRef: display iteration...")<<endl;


    for(auto const & MapElement : gCalibCounterMap)
    {
      UInt_t uCounterAddress = MapElement.first;
      TCounter* tCounter = MapElement.second;

      tCounter->ResetHistograms();
      tCounter->ZoomHistograms();

      if(!bCalibrateBRefCalibRef)
      {
        tCounter->ResetAllHistograms();
        tCounter->ZoomAllHistograms();
      }
    }


    for(Int_t iEvent = 0; iEvent < giNEvents; iEvent++)
    {
      if(0 == (iEvent + 1) % 100000)
      {
        cout<<Form("   event %4d", iEvent + 1)<<endl;
      }

      gInputDigiArray->Delete();
      gInputDigiBranch->GetEntry(iEvent);

      FillDigiArrays();

      for(auto const & MapElement : gCalibCounterMap)
      {
        UInt_t uCounterAddress = MapElement.first;
        TCounter* tCounter = MapElement.second;

        Bool_t bBRefVelocityKnown = kTRUE;
        Bool_t bCounterVelocityKnown = kFALSE;
        Bool_t bFillAllHistograms = (!bCalibrateBRefCalibRef ? kTRUE : kFALSE);

        if(uBRefAddress == uCounterAddress)
        {
          FillHistograms(tCounter, bBRefVelocityKnown, bCounterVelocityKnown, bFillAllHistograms);
        }
      }
    }


    for(auto const & MapElement : gCalibCounterMap)
    {
      UInt_t uCounterAddress = MapElement.first;
      TCounter* tCounter = MapElement.second;

      if(uBRefAddress == uCounterAddress)
      {
        if(!bCalibrateBRefCalibRef)
        {
          tCounter->PlotHistograms("FinalBRefCalib");
          tCounter->PlotOverview("FinalBRefCalib");
        }
      }
    }
  }

  // TODO: Call TCounter::UpdateCorrelationOffsets here?
  //       The difference between the TOF from the BRef counter to the center of
  //       counter 4-0-0, for example, and the TOF from the BRef counter to the
  //       periphery of counter 4-0-0 amounts to about 18 ps (16 cm flight path
  //       difference on 400 cm). This value is smaller than the measurement precision
  //       of the counter. For this reason, it does not make much sense to correct
  //       for slightly different flight pathes from the BRef reference cell to
  //       individually calibrated counter cells by summing up correlation offsets
  //       between neighboring cells until the counter reference cell is reached.
  //       Additionally, if the counter reference cell is positioned in the counter
  //       center, multiple error-prone correlation offsets between neighboring cells
  //       need to be summed up if a cell in the counter periphery is to be
  //       corrected this way.

}


void CalibrateCounters()
{
  cout<<Form("CalibrateCounters: initial delay correction...")<<endl;

  for(auto const & MapElement : gCalibCounterMap)
  {
    UInt_t uCounterAddress = MapElement.first;
    TCounter* tCounter = MapElement.second;

    if(gBRefAddressSet.find(uCounterAddress) == gBRefAddressSet.end())
    {
      tCounter->ClearCalibration();
      tCounter->ResetHistograms();
      tCounter->ResetAllHistograms();
    }
  }


  for(Int_t iEvent = 0; iEvent < giNEvents; iEvent++)
  {
    if(0 == (iEvent + 1) % 100000)
    {
      cout<<Form("   event %4d", iEvent + 1)<<endl;
    }

    gInputDigiArray->Delete();
    gInputDigiBranch->GetEntry(iEvent);

    FillDigiArrays();

    for(auto const & MapElement : gCalibCounterMap)
    {
      UInt_t uCounterAddress = MapElement.first;
      TCounter* tCounter = MapElement.second;

      Bool_t bBRefVelocityKnown = kTRUE;
      Bool_t bCounterVelocityKnown = kFALSE;
      Bool_t bFillAllHistograms = kTRUE;

      if(gBRefAddressSet.find(uCounterAddress) == gBRefAddressSet.end())
      {
        FillHistograms(tCounter, bBRefVelocityKnown, bCounterVelocityKnown, bFillAllHistograms);
      }
    }
  }


  for(auto const & MapElement : gCalibCounterMap)
  {
    UInt_t uCounterAddress = MapElement.first;
    TCounter* tCounter = MapElement.second;

    Bool_t bTimePeak = kTRUE;
    Bool_t bFitTimeDifference = kFALSE;
    Bool_t bLastFit = kFALSE;

    if(gBRefAddressSet.find(uCounterAddress) == gBRefAddressSet.end())
    {
      tCounter->UpdateDelayOffsets(bTimePeak, bFitTimeDifference, bLastFit);
      tCounter->PlotHistograms("OriginalCounterCalib");
      tCounter->PlotOverview("OriginalCounterCalib");
    }
  }


  if(!gbSkipCounterSlewingCorrections)
  {
    cout<<Form("CalibrateCounters: walk corrections...")<<endl;

    for(Int_t iCorrection = 0; iCorrection < giNWalkCorrections + 1; iCorrection++)
    {
      for(auto const & MapElement : gCalibCounterMap)
      {
        UInt_t uCounterAddress = MapElement.first;
        TCounter* tCounter = MapElement.second;

        if(gBRefAddressSet.find(uCounterAddress) == gBRefAddressSet.end())
        {
          tCounter->ResetHistograms();
          if(0 < iCorrection)
          {
            tCounter->ZoomHistograms();
          }
        }
      }


      for(Int_t iEvent = 0; iEvent < giNEvents; iEvent++)
      {
        if(0 == (iEvent + 1) % 100000)
        {
          cout<<Form("   event %4d", iEvent + 1)<<endl;
        }

        gInputDigiArray->Delete();
        gInputDigiBranch->GetEntry(iEvent);

        FillDigiArrays();

        for(auto const & MapElement : gCalibCounterMap)
        {
          UInt_t uCounterAddress = MapElement.first;
          TCounter* tCounter = MapElement.second;

          Bool_t bBRefVelocityKnown = kTRUE;
          Bool_t bCounterVelocityKnown = kFALSE;
          Bool_t bFillAllHistograms = kFALSE;

          if(gBRefAddressSet.find(uCounterAddress) == gBRefAddressSet.end())
          {
            FillHistograms(tCounter, bBRefVelocityKnown, bCounterVelocityKnown, bFillAllHistograms);
          }
        }
      }


      for(auto const & MapElement : gCalibCounterMap)
      {
        UInt_t uCounterAddress = MapElement.first;
        TCounter* tCounter = MapElement.second;

        if(gBRefAddressSet.find(uCounterAddress) == gBRefAddressSet.end() && tCounter->fbCorrectSlewing)
        {
          Bool_t bShiftPeak = kFALSE;

          if(0 == iCorrection)
          {
            bShiftPeak = kTRUE;
          }

          // Do not update slewing offsets in the final correction step
          if(giNWalkCorrections != iCorrection)
          {
            tCounter->UpdateSlewingOffsets(bShiftPeak, iCorrection%2);
          }

          if(0 == iCorrection)
          {
            tCounter->PlotHistograms("OriginalCounterWalk");
          }
          else if(giNWalkCorrections == iCorrection)
          {
            tCounter->PlotHistograms("FinalCounterWalk");
          }
        }
      }

    }
  }


  cout<<Form("CalibrateCounters: initial time difference fit correction...")<<endl;

  for(auto const & MapElement : gCalibCounterMap)
  {
    UInt_t uCounterAddress = MapElement.first;
    TCounter* tCounter = MapElement.second;

    if(gBRefAddressSet.find(uCounterAddress) == gBRefAddressSet.end())
    {
      tCounter->ResetHistograms();
      tCounter->ZoomHistograms();
    }
  }


  for(Int_t iEvent = 0; iEvent < giNEvents; iEvent++)
  {
    if(0 == (iEvent + 1) % 100000)
    {
      cout<<Form("   event %4d", iEvent + 1)<<endl;
    }

    gInputDigiArray->Delete();
    gInputDigiBranch->GetEntry(iEvent);

    FillDigiArrays();

    for(auto const & MapElement : gCalibCounterMap)
    {
      UInt_t uCounterAddress = MapElement.first;
      TCounter* tCounter = MapElement.second;

      Bool_t bBRefVelocityKnown = kTRUE;
      Bool_t bCounterVelocityKnown = kFALSE;
      Bool_t bFillAllHistograms = kFALSE;

      if(gBRefAddressSet.find(uCounterAddress) == gBRefAddressSet.end())
      {
        FillHistograms(tCounter, bBRefVelocityKnown, bCounterVelocityKnown, bFillAllHistograms);
      }
    }
  }


  for(auto const & MapElement : gCalibCounterMap)
  {
    UInt_t uCounterAddress = MapElement.first;
    TCounter* tCounter = MapElement.second;

    Bool_t bTimePeak = kFALSE;
    Bool_t bFitTimeDifference = kTRUE;
    Bool_t bLastFit = kFALSE;

    if(gBRefAddressSet.find(uCounterAddress) == gBRefAddressSet.end())
    {
      tCounter->UpdateDelayOffsets(bTimePeak, bFitTimeDifference, bLastFit);
    }
  }


  cout<<Form("CalibrateCounters: final time difference fit correction...")<<endl;

  for(auto const & MapElement : gCalibCounterMap)
  {
    UInt_t uCounterAddress = MapElement.first;
    TCounter* tCounter = MapElement.second;

    if(gBRefAddressSet.find(uCounterAddress) == gBRefAddressSet.end())
    {
      tCounter->ResetHistograms();
      tCounter->ZoomHistograms();
    }
  }


  for(Int_t iEvent = 0; iEvent < giNEvents; iEvent++)
  {
    if(0 == (iEvent + 1) % 100000)
    {
      cout<<Form("   event %4d", iEvent + 1)<<endl;
    }

    gInputDigiArray->Delete();
    gInputDigiBranch->GetEntry(iEvent);

    FillDigiArrays();

    for(auto const & MapElement : gCalibCounterMap)
    {
      UInt_t uCounterAddress = MapElement.first;
      TCounter* tCounter = MapElement.second;

      Bool_t bBRefVelocityKnown = kTRUE;
      Bool_t bCounterVelocityKnown = kFALSE;
      Bool_t bFillAllHistograms = kFALSE;

      if(gBRefAddressSet.find(uCounterAddress) == gBRefAddressSet.end())
      {
        FillHistograms(tCounter, bBRefVelocityKnown, bCounterVelocityKnown, bFillAllHistograms);
      }
    }
  }


  for(auto const & MapElement : gCalibCounterMap)
  {
    UInt_t uCounterAddress = MapElement.first;
    TCounter* tCounter = MapElement.second;

    Bool_t bTimePeak = kFALSE;
    Bool_t bFitTimeDifference = kTRUE;
    Bool_t bLastFit = kTRUE;

    if(gBRefAddressSet.find(uCounterAddress) == gBRefAddressSet.end())
    {
      tCounter->UpdateDelayOffsets(bTimePeak, bFitTimeDifference, bLastFit);
    }
  }


  cout<<Form("CalibrateCounters: display iteration...")<<endl;

  for(auto const & MapElement : gCalibCounterMap)
  {
    UInt_t uCounterAddress = MapElement.first;
    TCounter* tCounter = MapElement.second;

    if(gBRefAddressSet.find(uCounterAddress) == gBRefAddressSet.end())
    {
      tCounter->ResetHistograms();
      tCounter->ZoomHistograms();

      tCounter->ResetAllHistograms();
      tCounter->ZoomAllHistograms();
    }
  }


  for(Int_t iEvent = 0; iEvent < giNEvents; iEvent++)
  {
    if(0 == (iEvent + 1) % 100000)
    {
      cout<<Form("   event %4d", iEvent + 1)<<endl;
    }

    gInputDigiArray->Delete();
    gInputDigiBranch->GetEntry(iEvent);

    FillDigiArrays();

    for(auto const & MapElement : gCalibCounterMap)
    {
      UInt_t uCounterAddress = MapElement.first;
      TCounter* tCounter = MapElement.second;

      Bool_t bBRefVelocityKnown = kTRUE;
      Bool_t bCounterVelocityKnown = kTRUE;
      Bool_t bFillAllHistograms = kTRUE;

      if(gBRefAddressSet.find(uCounterAddress) == gBRefAddressSet.end())
      {
        FillHistograms(tCounter, bBRefVelocityKnown, bCounterVelocityKnown, bFillAllHistograms);
      }
    }
  }


  for(auto const & MapElement : gCalibCounterMap)
  {
    UInt_t uCounterAddress = MapElement.first;
    TCounter* tCounter = MapElement.second;

    if(gBRefAddressSet.find(uCounterAddress) == gBRefAddressSet.end())
    {
      tCounter->PlotHistograms("FinalCounterCalib");
      tCounter->PlotOverview("FinalCounterCalib");
    }
  }


  // TODO: correlation offsets? velocity offsets?
}


void CalibrateToT()
{
  cout<<Form("CalibrateToT: align ToT spectra...")<<endl;

  for(Int_t iEvent = 0; iEvent < giNEvents; iEvent++)
  {
    if(0 == (iEvent + 1) % 100000)
    {
      cout<<Form("   event %4d", iEvent + 1)<<endl;
    }

    gInputDigiArray->Delete();
    gInputDigiBranch->GetEntry(iEvent);

    FillDigiArrays();

    for(auto const & MapElement : gCalibCounterMap)
    {
      UInt_t uCounterAddress = MapElement.first;
      TCounter* tCounter = MapElement.second;

      tCounter->FillToTSpectra();
    }
  }


  for(auto const & MapElement : gCalibCounterMap)
  {
    UInt_t uCounterAddress = MapElement.first;
    TCounter* tCounter = MapElement.second;

    Bool_t bFinalAlignment = kFALSE;

    tCounter->UpdateToTAlignment();
    tCounter->PlotToTSpectra("InitialToTAlign", bFinalAlignment);
  }


  cout<<Form("CalibrateToT: display ToT alignment...")<<endl;

  for(auto const & MapElement : gCalibCounterMap)
  {
    UInt_t uCounterAddress = MapElement.first;
    TCounter* tCounter = MapElement.second;

    tCounter->ResetToTSpectra();
    tCounter->ZoomToTSpectra();
  }


  for(Int_t iEvent = 0; iEvent < giNEvents; iEvent++)
  {
    if(0 == (iEvent + 1) % 100000)
    {
      cout<<Form("   event %4d", iEvent + 1)<<endl;
    }

    gInputDigiArray->Delete();
    gInputDigiBranch->GetEntry(iEvent);

    FillDigiArrays();

    for(auto const & MapElement : gCalibCounterMap)
    {
      UInt_t uCounterAddress = MapElement.first;
      TCounter* tCounter = MapElement.second;

      tCounter->FillToTSpectra();
    }
  }


  for(auto const & MapElement : gCalibCounterMap)
  {
    UInt_t uCounterAddress = MapElement.first;
    TCounter* tCounter = MapElement.second;

    Bool_t bFinalAlignment = kTRUE;

    tCounter->PlotToTSpectra("FinalToTAlign", bFinalAlignment);
  }
}


void WriteCalibratedDigis()
{
  FairRootFileSink* tFileSink = new FairRootFileSink(TString(gcOutputFileName));
  FairRootManager* tRootManager = FairRootManager::Instance();
  FairRunAna* tRun = new FairRunAna();

  tRootManager->SetSink(tFileSink);
  tRootManager->InitSink();

  TClonesArray* tOutputDigiArray = new TClonesArray("CbmTofDigiExp", 1000);
  tRootManager->Register("CbmTofDigiExp", "Tof", tOutputDigiArray, kTRUE);

  tRootManager->WriteFolder();

  for(Int_t iEvent = 0; iEvent < giNEvents; iEvent++)
  {
    if(0 == (iEvent + 1) % 100000)
    {
      cout<<Form("   event %4d", iEvent + 1)<<endl;
    }

    gInputDigiArray->Delete();
    gInputDigiBranch->GetEntry(iEvent);

    Int_t iNDigis = gInputDigiArray->GetEntriesFast();

    for(Int_t iDigi = 0; iDigi < iNDigis; iDigi++)
    {
      CbmTofDigiExp* tDigi = dynamic_cast<CbmTofDigiExp*>(gInputDigiArray->At(iDigi));

      auto itCounter = gCalibCounterMap.find(static_cast<UInt_t>(CbmTofAddress::GetModFullId(static_cast<UInt_t>(tDigi->GetAddress()))));

      if(itCounter != gCalibCounterMap.end())
      {
        TCounter* tCounter = itCounter->second;
        CalibrateDigi(tDigi, tCounter);

        CbmTofDigiExp* tCalibratedDigi = new((*tOutputDigiArray)[iDigi]) CbmTofDigiExp(*tDigi);
      }
    }

    tOutputDigiArray->Sort(iNDigis);

    tRootManager->Fill();

    tOutputDigiArray->Delete();
  }

  tRootManager->Write();
  tFileSink->Close();

/*
  // TODO: code fragments kept here to allow for writing a 'FairFileHeader' object
  //       and a 'FairEventHeader' branch to the output file if ever required

  FairRunIdGenerator genid;
  fRunId = genid.generateId();

  FairEventHeader* tEvtHeader = new FairEventHeader();
  feh->SetEventTime(fEventTime);
  feh->SetRunId(fEvtHeader->GetRunId());
  feh->SetMCEntryNumber(fEvtHeader->GetMCEntryNumber());
  feh->SetInputFileId(0);

  fEvtHeader->Register(fStoreEventHeader);
  fRootManager->FillEventHeader(fEvtHeader);
  fFileHeader->SetRunId(fRunId);
  fRootManager->WriteFileHeader(fFileHeader);
*/
}


void FillDigiArrays()
{
  for(auto const & MapElement : gCalibCounterMap)
  {
    TCounter* tCounter = MapElement.second;
    tCounter->ClearDigiArrays();
  }

  Int_t iNDigis = gInputDigiArray->GetEntriesFast();

  for(Int_t iDigi = 0; iDigi < iNDigis; iDigi++)
  {
    CbmTofDigiExp* tDigi = dynamic_cast<CbmTofDigiExp*>(gInputDigiArray->At(iDigi));

    if(!gbNanoSeconds)
    {
      tDigi->SetTime(tDigi->GetTime()/1000.);
      tDigi->SetTot(tDigi->GetTot()/1000.);
    }

    tDigi->SetTot(tDigi->GetTot() + gdDigiToTOffset);

    auto itCounter = gCalibCounterMap.find(static_cast<UInt_t>(CbmTofAddress::GetModFullId(static_cast<UInt_t>(tDigi->GetAddress()))));

    if(itCounter != gCalibCounterMap.end())
    {
      TCounter* tCounter = itCounter->second;
      if(!tCounter->fbDigiAvailable[static_cast<Int_t>(tDigi->GetChannel())][static_cast<Int_t>(tDigi->GetSide())])
      {
        tCounter->fbDigiAvailable[static_cast<Int_t>(tDigi->GetChannel())][static_cast<Int_t>(tDigi->GetSide())] = kTRUE;
        tCounter->fDigi[static_cast<Int_t>(tDigi->GetChannel())][static_cast<Int_t>(tDigi->GetSide())] = tDigi;
      }
    }
  }

  for(auto const & MapElement : gCalibCounterMap)
  {
    TCounter* tCounter = MapElement.second;
    tCounter->ApplyToTScaling();
    tCounter->ComputeDigiWalk();
  }
}


void CalibrateDigi(CbmTofDigiExp* tDigi, TCounter* tCounter)
{
  if(!gbNanoSeconds)
  {
    tDigi->SetTime(tDigi->GetTime()/1000.);
    tDigi->SetTot(tDigi->GetTot()/1000.);
  }

  tDigi->SetTot(tDigi->GetTot() + gdDigiToTOffset);

  tCounter->CalibrateDigi(tDigi);

  if(!gbNanoSeconds)
  {
    tDigi->SetTime(tDigi->GetTime()*1000.);
    tDigi->SetTot(tDigi->GetTot()*1000.);
  }
}


void FillHistograms(TCounter* tCounter, Bool_t bBRefVelocityKnown, Bool_t bCounterVelocityKnown, Bool_t bFillAllHistograms)
{
  if(tCounter->fbIsBRef)
  {
    TCounter* tCalibCounter = gCalibCounterMap.at(static_cast<UInt_t>(CbmTofAddress::GetModFullId(CbmTofAddress::GetUniqueAddress(giBRefCalibRefModuleIndex, giBRefCalibRefCounterIndex, 0, 0, giBRefCalibRefModuleType))));
    Int_t iCalibCell = tCalibCounter->fiReferenceCell;

    if(tCalibCounter->fbDigiAvailable[iCalibCell][0] && tCalibCounter->fbDigiAvailable[iCalibCell][1])
    {
      for(Int_t iCell = 0; iCell < tCounter->fiNCells; iCell++)
      {
        if(!tCounter->fbDeadCell[iCell])
        {
          if(gbBRefPadReadout)
          {
            if(tCounter->fbDigiAvailable[iCell][0])
            {
              Double_t dTimeResidual = 0.5*( tCalibCounter->fDigi[iCalibCell][0]->GetTime()
                                            -tCalibCounter->fdSignalPropagationDelay[iCalibCell][0]
                                            -tCalibCounter->fdDigiSlewing[iCalibCell][0]

                                            +tCalibCounter->fDigi[iCalibCell][1]->GetTime()
                                            -tCalibCounter->fdSignalPropagationDelay[iCalibCell][1]
                                            -tCalibCounter->fdDigiSlewing[iCalibCell][1]
                                           )

                                        // The term is commented out here as the BRef counter just needs to be positioned somewhere
                                        // in time with respect to a stable reference cell of another counter. For this positioning,
                                        // it is not required that the signal runtime on the reference cell is taken into account.
                                        // In a later step, the counter that provides the reference cell is placed with respect
                                        // to the BRef counter CONSIDERING the runtime on its cells.
//                                        -0.5*tCalibCounter->fdCellAlongWidth/tCalibCounter->fdMeanSignalPropagationVelocity

                                      -( tCounter->fDigi[iCell][0]->GetTime()
                                        -tCounter->fdSignalPropagationDelay[iCell][0]
                                        -tCounter->fdDigiSlewing[iCell][0]

                                        // The (fixed) values of the parameters 'gdBRefMeanPosition' and 'fdMeanSignalPropagationVelocity'
                                        // cannot be derived from the data for a BRef counter which is read out on one strip end only.
                                        // Anyway, the term is commented out for the reason given above.
//                                        -(0.5*tCounter->fdCellAlongWidth + gdBRefMeanPosition)/tCounter->fdMeanSignalPropagationVelocity
                                       );


              Double_t dTimeDifference = 0.;

              tCounter->fTimeResidual[iCell]->Fill(dTimeResidual);
              tCounter->fTimeDifference[iCell]->Fill(dTimeDifference);
              tCounter->fTimeResidualWalk[iCell][0]->Fill(tCounter->fDigi[iCell][0]->GetTot(), dTimeResidual);

              if(bFillAllHistograms)
              {
                tCounter->fAllTimeResidual->Fill(iCell, dTimeResidual);
                tCounter->fAllTimeDifference->Fill(iCell, dTimeDifference);
                tCounter->fAllLocalPosition->Fill(iCell, dTimeDifference*tCounter->fdSignalPropagationVelocity[iCell]);
                tCounter->fAllLocalPositionMeanVelocity->Fill(iCell, dTimeDifference*tCounter->fdMeanSignalPropagationVelocity);
              }


              Int_t iNeighboringCell = -1;

              if(iCell < tCounter->fiReferenceCell)
              {
                iNeighboringCell = iCell + 1;

                while(iNeighboringCell < tCounter->fiReferenceCell)
                {
                  if(!tCounter->fbDeadCell[iNeighboringCell])
                  {
                    break;
                  }

                  iNeighboringCell++;
                }
              }
              else if(iCell > tCounter->fiReferenceCell)
              {
                iNeighboringCell = iCell - 1;

                while(iNeighboringCell > tCounter->fiReferenceCell)
                {
                  if(!tCounter->fbDeadCell[iNeighboringCell])
                  {
                    break;
                  }

                  iNeighboringCell--;
                }
              }

              if(-1 < iNeighboringCell && tCounter->fbDigiAvailable[iNeighboringCell][0])
              {
                Double_t dTimeMeanCorrelation =  tCounter->fDigi[iCell][0]->GetTime()
                                                -tCounter->fdSignalPropagationDelay[iCell][0]
                                                -tCounter->fdDigiSlewing[iCell][0]

                                                // The (fixed) values of the parameters 'gdBRefMeanPosition' and 'fdMeanSignalPropagationVelocity'
                                                // cannot be derived from the data for a BRef counter which is read out on one strip end only.
                                                // Thus, the mean signal runtime can unscrupulously be subtracted here.
                                                -(0.5*tCounter->fdCellAlongWidth + gdBRefMeanPosition)/tCounter->fdMeanSignalPropagationVelocity

                                                -( tCounter->fDigi[iNeighboringCell][0]->GetTime()
                                                  -tCounter->fdSignalPropagationDelay[iNeighboringCell][0]
                                                  -tCounter->fdDigiSlewing[iNeighboringCell][0]

                                                  // The (fixed) values of the parameters 'gdBRefMeanPosition' and 'fdMeanSignalPropagationVelocity'
                                                  // cannot be derived from the data for a BRef counter which is read out on one strip end only.
                                                  // Thus, the mean signal runtime can unscrupulously be subtracted here.
                                                  -(0.5*tCounter->fdCellAlongWidth + gdBRefMeanPosition)/tCounter->fdMeanSignalPropagationVelocity
                                                 );

                Double_t dTimeDifferenceCorrelation = 0.;

                tCounter->fTimeMeanCorrelation[iCell]->Fill(dTimeMeanCorrelation);
                tCounter->fTimeDifferenceCorrelation[iCell]->Fill(dTimeDifferenceCorrelation);

                if(bFillAllHistograms)
                {
                  tCounter->fAllTimeMeanCorrelation->Fill(iCell, dTimeMeanCorrelation);
                  tCounter->fAllTimeDifferenceCorrelation->Fill(iCell, dTimeDifferenceCorrelation);
                }
              }
            }
          }
          else
          {
            if(tCounter->fbDigiAvailable[iCell][0] && tCounter->fbDigiAvailable[iCell][1])
            {
              Double_t dTimeResidual = 0.5*( tCalibCounter->fDigi[iCalibCell][0]->GetTime()
                                            -tCalibCounter->fdSignalPropagationDelay[iCalibCell][0]
                                            -tCalibCounter->fdDigiSlewing[iCalibCell][0]

                                            +tCalibCounter->fDigi[iCalibCell][1]->GetTime()
                                            -tCalibCounter->fdSignalPropagationDelay[iCalibCell][1]
                                            -tCalibCounter->fdDigiSlewing[iCalibCell][1]
                                           )

                                        // The term is commented out here as the BRef counter just needs to be positioned somewhere
                                        // in time with respect to a stable reference cell of another counter. For this positioning,
                                        // it is not required that the signal runtime on the reference cell is taken into account.
                                        // In a later step, the counter that provides the reference cell is placed with respect
                                        // to the BRef counter CONSIDERING the runtime on its cells.
//                                        -0.5*tCalibCounter->fdCellAlongWidth/tCalibCounter->fdMeanSignalPropagationVelocity

                                      -( 0.5*( tCounter->fDigi[iCell][0]->GetTime()
                                              -tCounter->fdSignalPropagationDelay[iCell][0]
                                              -tCounter->fdDigiSlewing[iCell][0]

                                              +tCounter->fDigi[iCell][1]->GetTime()
                                              -tCounter->fdSignalPropagationDelay[iCell][1]
                                              -tCounter->fdDigiSlewing[iCell][1]
                                             )
                                        // The term is commented out for the reason given above.
//                                        -(bBRefVelocityKnown ? 1. : 0.)*0.5*tCounter->fdCellAlongWidth/tCounter->fdMeanSignalPropagationVelocity
                                       );

              Double_t dTimeDifference = 0.5*( tCounter->fDigi[iCell][0]->GetTime()
                                              -tCounter->fdSignalPropagationDelay[iCell][0]
                                              -tCounter->fdDigiSlewing[iCell][0]

                                              -( tCounter->fDigi[iCell][1]->GetTime()
                                                -tCounter->fdSignalPropagationDelay[iCell][1]
                                                -tCounter->fdDigiSlewing[iCell][1]
                                               )
                                             );

              tCounter->fTimeResidual[iCell]->Fill(dTimeResidual);
              tCounter->fTimeDifference[iCell]->Fill(dTimeDifference);
              tCounter->fTimeResidualWalk[iCell][0]->Fill(tCounter->fDigi[iCell][0]->GetTot(), dTimeResidual);
              tCounter->fTimeResidualWalk[iCell][1]->Fill(tCounter->fDigi[iCell][1]->GetTot(), dTimeResidual);

              if(bFillAllHistograms)
              {
                tCounter->fAllTimeResidual->Fill(iCell, dTimeResidual);
                tCounter->fAllTimeDifference->Fill(iCell, dTimeDifference);
                tCounter->fAllLocalPosition->Fill(iCell, dTimeDifference*tCounter->fdSignalPropagationVelocity[iCell]);
                tCounter->fAllLocalPositionMeanVelocity->Fill(iCell, dTimeDifference*tCounter->fdMeanSignalPropagationVelocity);
              }


              Int_t iNeighboringCell = -1;

              if(iCell < tCounter->fiReferenceCell)
              {
                iNeighboringCell = iCell + 1;

                while(iNeighboringCell < tCounter->fiReferenceCell)
                {
                  if(!tCounter->fbDeadCell[iNeighboringCell])
                  {
                    break;
                  }

                  iNeighboringCell++;
                }
              }
              else if(iCell > tCounter->fiReferenceCell)
              {
                iNeighboringCell = iCell - 1;

                while(iNeighboringCell > tCounter->fiReferenceCell)
                {
                  if(!tCounter->fbDeadCell[iNeighboringCell])
                  {
                    break;
                  }

                  iNeighboringCell--;
                }
              }

              if(-1 < iNeighboringCell && tCounter->fbDigiAvailable[iNeighboringCell][0] && tCounter->fbDigiAvailable[iNeighboringCell][1])
              {
                Double_t dTimeMeanCorrelation = 0.5*( tCounter->fDigi[iCell][0]->GetTime()
                                                     -tCounter->fdSignalPropagationDelay[iCell][0]
                                                     -tCounter->fdDigiSlewing[iCell][0]

                                                    +tCounter->fDigi[iCell][1]->GetTime()
                                                    -tCounter->fdSignalPropagationDelay[iCell][1]
                                                    -tCounter->fdDigiSlewing[iCell][1]
                                                    )

                                               -(bBRefVelocityKnown ? 1. : 0.)*0.5*tCounter->fdCellAlongWidth/tCounter->fdMeanSignalPropagationVelocity

                                               -( 0.5*( tCounter->fDigi[iNeighboringCell][0]->GetTime()
                                                       -tCounter->fdSignalPropagationDelay[iNeighboringCell][0]
                                                       -tCounter->fdDigiSlewing[iNeighboringCell][0]

                                                       +tCounter->fDigi[iNeighboringCell][1]->GetTime()
                                                       -tCounter->fdSignalPropagationDelay[iNeighboringCell][1]
                                                       -tCounter->fdDigiSlewing[iNeighboringCell][1]
                                                      )

                                                 -(bBRefVelocityKnown ? 1. : 0.)*0.5*tCounter->fdCellAlongWidth/tCounter->fdMeanSignalPropagationVelocity
                                                );


                Double_t dTimeDifferenceCorrelation = 0.5*( tCounter->fDigi[iCell][0]->GetTime()
                                                           -tCounter->fdSignalPropagationDelay[iCell][0]
                                                           -tCounter->fdDigiSlewing[iCell][0]

                                                           -( tCounter->fDigi[iCell][1]->GetTime()
                                                             -tCounter->fdSignalPropagationDelay[iCell][1]
                                                             -tCounter->fdDigiSlewing[iCell][1]
                                                            )
                                                          )

                                                     -0.5*( tCounter->fDigi[iNeighboringCell][0]->GetTime()
                                                           -tCounter->fdSignalPropagationDelay[iNeighboringCell][0]
                                                           -tCounter->fdDigiSlewing[iNeighboringCell][0]

                                                           -( tCounter->fDigi[iNeighboringCell][1]->GetTime()
                                                             -tCounter->fdSignalPropagationDelay[iNeighboringCell][1]
                                                             -tCounter->fdDigiSlewing[iNeighboringCell][1]
                                                            )
                                                          );

                tCounter->fTimeMeanCorrelation[iCell]->Fill(dTimeMeanCorrelation);
                tCounter->fTimeDifferenceCorrelation[iCell]->Fill(dTimeDifferenceCorrelation);

                if(bFillAllHistograms)
                {
                  tCounter->fAllTimeMeanCorrelation->Fill(iCell, dTimeMeanCorrelation);
                  tCounter->fAllTimeDifferenceCorrelation->Fill(iCell, dTimeDifferenceCorrelation);
                }
              }
            }
          }
        }
      }
    }
  }
  else
  {
    TCounter* tCalibCounter = gCalibCounterMap.at(static_cast<UInt_t>(CbmTofAddress::GetModFullId(CbmTofAddress::GetUniqueAddress(giBRefModuleIndex, giBRefCounterIndex, 0, 0, giBRefModuleType))));
    Int_t iCalibCell = tCalibCounter->fiReferenceCell;

    if(gbBRefPadReadout)
    {
      if(tCalibCounter->fbDigiAvailable[iCalibCell][0])
      {
        for(Int_t iCell = 0; iCell < tCounter->fiNCells; iCell++)
        {
          if(!tCounter->fbDeadCell[iCell])
          {
            if(tCounter->fbDigiAvailable[iCell][0] && tCounter->fbDigiAvailable[iCell][1])
            {
              Double_t dTimeResidual = 0.5*( tCounter->fDigi[iCell][0]->GetTime()
                                            -tCounter->fdSignalPropagationDelay[iCell][0]
                                            -tCounter->fdDigiSlewing[iCell][0]

                                            +tCounter->fDigi[iCell][1]->GetTime()
                                            -tCounter->fdSignalPropagationDelay[iCell][1]
                                            -tCounter->fdDigiSlewing[iCell][1]
                                           )

                                      -(bCounterVelocityKnown ? 1. : 0.)*0.5*tCounter->fdCellAlongWidth/tCounter->fdMeanSignalPropagationVelocity

                                      -( tCalibCounter->fDigi[iCalibCell][0]->GetTime()
                                        -tCalibCounter->fdSignalPropagationDelay[iCalibCell][0]
                                        -tCalibCounter->fdDigiSlewing[iCalibCell][0]

                                        -(0.5*tCalibCounter->fdCellAlongWidth + gdBRefMeanPosition)/tCalibCounter->fdMeanSignalPropagationVelocity
                                       );

              Double_t dTimeDifference = 0.5*( tCounter->fDigi[iCell][0]->GetTime()
                                              -tCounter->fdSignalPropagationDelay[iCell][0]
                                              -tCounter->fdDigiSlewing[iCell][0]

                                              -( tCounter->fDigi[iCell][1]->GetTime()
                                                -tCounter->fdSignalPropagationDelay[iCell][1]
                                                -tCounter->fdDigiSlewing[iCell][1]
                                               )
                                             );

              tCounter->fTimeResidual[iCell]->Fill(dTimeResidual);
              tCounter->fTimeDifference[iCell]->Fill(dTimeDifference);
              tCounter->fTimeResidualWalk[iCell][0]->Fill(tCounter->fDigi[iCell][0]->GetTot(), dTimeResidual);
              tCounter->fTimeResidualWalk[iCell][1]->Fill(tCounter->fDigi[iCell][1]->GetTot(), dTimeResidual);

              if(bFillAllHistograms)
              {
                tCounter->fAllTimeResidual->Fill(iCell, dTimeResidual);
                tCounter->fAllTimeDifference->Fill(iCell, dTimeDifference);
                tCounter->fAllLocalPosition->Fill(iCell, dTimeDifference*tCounter->fdSignalPropagationVelocity[iCell]);
                tCounter->fAllLocalPositionMeanVelocity->Fill(iCell, dTimeDifference*tCounter->fdMeanSignalPropagationVelocity);
              }


              Int_t iNeighboringCell = -1;

              if(iCell < tCounter->fiReferenceCell)
              {
                iNeighboringCell = iCell + 1;

                while(iNeighboringCell < tCounter->fiReferenceCell)
                {
                  if(!tCounter->fbDeadCell[iNeighboringCell])
                  {
                    break;
                  }

                  iNeighboringCell++;
                }
              }
              else if(iCell > tCounter->fiReferenceCell)
              {
                iNeighboringCell = iCell - 1;

                while(iNeighboringCell > tCounter->fiReferenceCell)
                {
                  if(!tCounter->fbDeadCell[iNeighboringCell])
                  {
                    break;
                  }

                  iNeighboringCell--;
                }
              }

              if(-1 < iNeighboringCell && tCounter->fbDigiAvailable[iNeighboringCell][0] && tCounter->fbDigiAvailable[iNeighboringCell][1])
              {
                Double_t dTimeMeanCorrelation = 0.5*( tCounter->fDigi[iCell][0]->GetTime()
                                                     -tCounter->fdSignalPropagationDelay[iCell][0]
                                                     -tCounter->fdDigiSlewing[iCell][0]

                                                    +tCounter->fDigi[iCell][1]->GetTime()
                                                    -tCounter->fdSignalPropagationDelay[iCell][1]
                                                    -tCounter->fdDigiSlewing[iCell][1]
                                                    )

                                               -0.5*tCounter->fdCellAlongWidth/tCounter->fdMeanSignalPropagationVelocity

                                               -( 0.5*( tCounter->fDigi[iNeighboringCell][0]->GetTime()
                                                       -tCounter->fdSignalPropagationDelay[iNeighboringCell][0]
                                                       -tCounter->fdDigiSlewing[iNeighboringCell][0]

                                                       +tCounter->fDigi[iNeighboringCell][1]->GetTime()
                                                       -tCounter->fdSignalPropagationDelay[iNeighboringCell][1]
                                                       -tCounter->fdDigiSlewing[iNeighboringCell][1]
                                                      )

                                                 -0.5*tCounter->fdCellAlongWidth/tCounter->fdMeanSignalPropagationVelocity
                                                );


                Double_t dTimeDifferenceCorrelation = 0.5*( tCounter->fDigi[iCell][0]->GetTime()
                                                           -tCounter->fdSignalPropagationDelay[iCell][0]
                                                           -tCounter->fdDigiSlewing[iCell][0]

                                                           -( tCounter->fDigi[iCell][1]->GetTime()
                                                             -tCounter->fdSignalPropagationDelay[iCell][1]
                                                             -tCounter->fdDigiSlewing[iCell][1]
                                                            )
                                                          )

                                                     -0.5*( tCounter->fDigi[iNeighboringCell][0]->GetTime()
                                                           -tCounter->fdSignalPropagationDelay[iNeighboringCell][0]
                                                           -tCounter->fdDigiSlewing[iNeighboringCell][0]

                                                           -( tCounter->fDigi[iNeighboringCell][1]->GetTime()
                                                             -tCounter->fdSignalPropagationDelay[iNeighboringCell][1]
                                                             -tCounter->fdDigiSlewing[iNeighboringCell][1]
                                                            )
                                                          );

                tCounter->fTimeMeanCorrelation[iCell]->Fill(dTimeMeanCorrelation);
                tCounter->fTimeDifferenceCorrelation[iCell]->Fill(dTimeDifferenceCorrelation);

                if(bFillAllHistograms)
                {
                  tCounter->fAllTimeMeanCorrelation->Fill(iCell, dTimeMeanCorrelation);
                  tCounter->fAllTimeDifferenceCorrelation->Fill(iCell, dTimeDifferenceCorrelation);
                }
              }
            }
          }
        }
      }
    }
    else
    {
      if(tCalibCounter->fbDigiAvailable[iCalibCell][0] && tCounter->fbDigiAvailable[iCalibCell][1])
      {
        for(Int_t iCell = 0; iCell < tCounter->fiNCells; iCell++)
        {
          if(!tCounter->fbDeadCell[iCell])
          {
            if(tCounter->fbDigiAvailable[iCell][0] && tCounter->fbDigiAvailable[iCell][1])
            {
              Double_t dTimeResidual = 0.5*( tCounter->fDigi[iCell][0]->GetTime()
                                            -tCounter->fdSignalPropagationDelay[iCell][0]
                                            -tCounter->fdDigiSlewing[iCell][0]

                                            +tCounter->fDigi[iCell][1]->GetTime()
                                            -tCounter->fdSignalPropagationDelay[iCell][1]
                                            -tCounter->fdDigiSlewing[iCell][1]
                                           )

                                      -(bCounterVelocityKnown ? 1. : 0.)*0.5*tCounter->fdCellAlongWidth/tCounter->fdMeanSignalPropagationVelocity

                                      -( 0.5*( tCalibCounter->fDigi[iCalibCell][0]->GetTime()
                                              -tCalibCounter->fdSignalPropagationDelay[iCalibCell][0]
                                              -tCalibCounter->fdDigiSlewing[iCalibCell][0]

                                              +tCalibCounter->fDigi[iCalibCell][1]->GetTime()
                                              -tCalibCounter->fdSignalPropagationDelay[iCalibCell][1]
                                              -tCalibCounter->fdDigiSlewing[iCalibCell][1]
                                             )

                                        -(bBRefVelocityKnown ? 1. : 0.)*0.5*tCalibCounter->fdCellAlongWidth/tCalibCounter->fdMeanSignalPropagationVelocity
                                       );


              Double_t dTimeDifference = 0.5*( tCounter->fDigi[iCell][0]->GetTime()
                                              -tCounter->fdSignalPropagationDelay[iCell][0]
                                              -tCounter->fdDigiSlewing[iCell][0]

                                              -( tCounter->fDigi[iCell][1]->GetTime()
                                                -tCounter->fdSignalPropagationDelay[iCell][1]
                                                -tCounter->fdDigiSlewing[iCell][1]
                                               )
                                             );

              tCounter->fTimeResidual[iCell]->Fill(dTimeResidual);
              tCounter->fTimeDifference[iCell]->Fill(dTimeDifference);
              tCounter->fTimeResidualWalk[iCell][0]->Fill(tCounter->fDigi[iCell][0]->GetTot(), dTimeResidual);
              tCounter->fTimeResidualWalk[iCell][1]->Fill(tCounter->fDigi[iCell][1]->GetTot(), dTimeResidual);

              if(bFillAllHistograms)
              {
                tCounter->fAllTimeResidual->Fill(iCell, dTimeResidual);
                tCounter->fAllTimeDifference->Fill(iCell, dTimeDifference);
                tCounter->fAllLocalPosition->Fill(iCell, dTimeDifference*tCounter->fdSignalPropagationVelocity[iCell]);
                tCounter->fAllLocalPositionMeanVelocity->Fill(iCell, dTimeDifference*tCounter->fdMeanSignalPropagationVelocity);
              }


              Int_t iNeighboringCell = -1;

              if(iCell < tCounter->fiReferenceCell)
              {
                iNeighboringCell = iCell + 1;

                while(iNeighboringCell < tCounter->fiReferenceCell)
                {
                  if(!tCounter->fbDeadCell[iNeighboringCell])
                  {
                    break;
                  }

                  iNeighboringCell++;
                }
              }
              else if(iCell > tCounter->fiReferenceCell)
              {
                iNeighboringCell = iCell - 1;

                while(iNeighboringCell > tCounter->fiReferenceCell)
                {
                  if(!tCounter->fbDeadCell[iNeighboringCell])
                  {
                    break;
                  }

                  iNeighboringCell--;
                }
              }

              if(-1 < iNeighboringCell && tCounter->fbDigiAvailable[iNeighboringCell][0] && tCounter->fbDigiAvailable[iNeighboringCell][1])
              {
                Double_t dTimeMeanCorrelation = 0.5*( tCounter->fDigi[iCell][0]->GetTime()
                                                     -tCounter->fdSignalPropagationDelay[iCell][0]
                                                     -tCounter->fdDigiSlewing[iCell][0]

                                                    +tCounter->fDigi[iCell][1]->GetTime()
                                                    -tCounter->fdSignalPropagationDelay[iCell][1]
                                                    -tCounter->fdDigiSlewing[iCell][1]
                                                    )

                                               -0.5*tCounter->fdCellAlongWidth/tCounter->fdMeanSignalPropagationVelocity

                                               -( 0.5*( tCounter->fDigi[iNeighboringCell][0]->GetTime()
                                                       -tCounter->fdSignalPropagationDelay[iNeighboringCell][0]
                                                       -tCounter->fdDigiSlewing[iNeighboringCell][0]

                                                       +tCounter->fDigi[iNeighboringCell][1]->GetTime()
                                                       -tCounter->fdSignalPropagationDelay[iNeighboringCell][1]
                                                       -tCounter->fdDigiSlewing[iNeighboringCell][1]
                                                      )

                                                 -0.5*tCounter->fdCellAlongWidth/tCounter->fdMeanSignalPropagationVelocity
                                                );


                Double_t dTimeDifferenceCorrelation = 0.5*( tCounter->fDigi[iCell][0]->GetTime()
                                                           -tCounter->fdSignalPropagationDelay[iCell][0]
                                                           -tCounter->fdDigiSlewing[iCell][0]

                                                           -( tCounter->fDigi[iCell][1]->GetTime()
                                                             -tCounter->fdSignalPropagationDelay[iCell][1]
                                                             -tCounter->fdDigiSlewing[iCell][1]
                                                            )
                                                          )

                                                     -0.5*( tCounter->fDigi[iNeighboringCell][0]->GetTime()
                                                           -tCounter->fdSignalPropagationDelay[iNeighboringCell][0]
                                                           -tCounter->fdDigiSlewing[iNeighboringCell][0]

                                                           -( tCounter->fDigi[iNeighboringCell][1]->GetTime()
                                                             -tCounter->fdSignalPropagationDelay[iNeighboringCell][1]
                                                             -tCounter->fdDigiSlewing[iNeighboringCell][1]
                                                            )
                                                          );

                tCounter->fTimeMeanCorrelation[iCell]->Fill(dTimeMeanCorrelation);
                tCounter->fTimeDifferenceCorrelation[iCell]->Fill(dTimeDifferenceCorrelation);

                if(bFillAllHistograms)
                {
                  tCounter->fAllTimeMeanCorrelation->Fill(iCell, dTimeMeanCorrelation);
                  tCounter->fAllTimeDifferenceCorrelation->Fill(iCell, dTimeDifferenceCorrelation);
                }
              }
            }
          }
        }
      }
    }
  }
}


Bool_t AddCounter(Int_t iModuleType, Int_t iModuleIndex, Int_t iCounterIndex)
{
  cout<<Form("Adding counter %d-%d-%d ...", iModuleType, iModuleIndex, iCounterIndex)<<endl;

  UInt_t uCounterAddress = CbmTofAddress::GetUniqueAddress(iModuleIndex, iCounterIndex, 0, 0, iModuleType);

  auto itCounter = gAllCounterMap.find(uCounterAddress);
  auto itInformation = gCounterInformation.find(std::make_tuple(iModuleType, iModuleIndex, iCounterIndex));

  if(itCounter != gAllCounterMap.end() && itInformation != gCounterInformation.end())
  {
    TCounter* tCounter = itCounter->second;

    Int_t iReferenceCell = std::get<0>(itInformation->second);
    UInt_t uDeadCells = std::get<1>(itInformation->second);
    Bool_t bCorrectSlewing = std::get<2>(itInformation->second);
    Bool_t bIsBRef = std::get<3>(itInformation->second);
    Double_t dInputSignalPropagationVelocity = std::get<4>(itInformation->second);

    tCounter->fiReferenceCell = iReferenceCell;
    tCounter->fbCorrectSlewing = bCorrectSlewing;
    tCounter->fbIsBRef = bIsBRef;
    tCounter->fdInputSignalPropagationVelocity = dInputSignalPropagationVelocity;

    for(Int_t iCell = 0; iCell < tCounter->fiNCells; iCell++)
    {
      tCounter->fbDeadCell[iCell] = (uDeadCells >> iCell ) & 0x1;
    }

    if(guBRefAddress == uCounterAddress && !bIsBRef)
    {
      cout<<Form("BRef counter %d-%d-%d not marked as BRef in 'gCounterInformation'!", iModuleType, iModuleIndex, iCounterIndex)<<endl;

      return kFALSE;
    }

    if(guBRefCalibRefAddress == uCounterAddress && bIsBRef)
    {
      cout<<Form("MRef counter %d-%d-%d marked as BRef in 'gCounterInformation'!", iModuleType, iModuleIndex, iCounterIndex)<<endl;

      return kFALSE;
    }

    gCalibCounterMap.emplace(uCounterAddress, tCounter);

    if(bIsBRef)
    {
      gBRefAddressSet.emplace(uCounterAddress);
    }
  }
  else
  {
    if(itCounter == gAllCounterMap.end())
    {
      cout<<Form("Counter %d-%d-%d not found in geometry file!", iModuleType, iModuleIndex, iCounterIndex)<<endl;
    }

    if(itInformation == gCounterInformation.end())
    {
      cout<<Form("No cell information available for counter %d-%d-%d!", iModuleType, iModuleIndex, iCounterIndex)<<endl;
    }

    return kFALSE;
  }

  return kTRUE;
}


Bool_t FindTofNode()
{
  TFile::Open(TString(gSystem->Getenv("VMCWORKDIR")) + "/geometry" + "/tof/geofile_tof_" + gcGeoVersion + ".root", "read");

  if(!gFile)
  {
    cout<<"Could not open geometry file. Abort macro execution."<<endl;
    return kFALSE;
  }

  gFile->FindObjectAny("FAIRGeom");
  if(!gGeoManager)
  {
    cout<<"Could not retrieve TGeoManager 'FAIRGeom' from geometry file."<<endl;
    return kFALSE;
  }

  gGeoHandler = new CbmTofGeoHandler();

  if(!gGeoHandler->Init(kFALSE))
  {
    cout<<"Failed to initialize 'CbmTofGeoHandler'."<<endl;
    return kFALSE;
  }

  Bool_t bFoundTofNode(kFALSE);
  gGeoManager->CdTop();

  TObjArray* tNodes = gGeoManager->GetCurrentNode()->GetNodes();
  for(Int_t iNode = 0; iNode < tNodes->GetEntriesFast(); iNode++)
  {
    TGeoNode* tNode = dynamic_cast<TGeoNode*>(tNodes->At(iNode));
    if(TString(tNode->GetName()).Contains("tof",TString::kIgnoreCase))
    {
      gGeoManager->CdDown(iNode);
      bFoundTofNode = kTRUE;
      break;
    }
  }

  if(!bFoundTofNode)
  {
    cout<<"Could not retrieve 'tof' node from TGeoManager."<<endl;
    return kFALSE;
  }

  gTofNodePath = gGeoManager->GetPath();

  return kTRUE;
}


void ExpandNode(TGeoNode* tMotherNode)
{
  TObjArray* tDaughterNodesList = tMotherNode->GetNodes();

  for(Int_t iNode = 0; iNode < tDaughterNodesList->GetEntriesFast(); iNode++)
  {
    TGeoNode* tDaughterNode = dynamic_cast<TGeoNode*>(tDaughterNodesList->At(iNode));

    // Extract the current module type and module index from the module node
    if(TString(tDaughterNode->GetName()).Contains("module"))
    {
      giCurrentModuleType = -1;
      giCurrentModuleIndex = -1;
      giCurrentCounterIndex = -1;

      gCurrentNodePath = gTofNodePath + "/" + (TString)tDaughterNode->GetName();

      boost::regex rgx(".*_(\\d+)_.*");
      boost::cmatch match;
      if( boost::regex_search(tDaughterNode->GetName(), match, rgx) )
      {
        giCurrentModuleType = boost::lexical_cast<Int_t>(match[1]);
      }

      giCurrentModuleIndex = tDaughterNode->GetNumber();

      gCurrentBoxNodePath = gCurrentNodePath;
    }
    // Extract the current counter index from the counter node and find the
    // central node in the counter volume (a glass plate OR a gap) to serve as
    // reference coordinate system
    else if(TString(tDaughterNode->GetName()).Contains("counter"))
    {
      gCurrentNodePath += "/" + (TString)tDaughterNode->GetName();

      giCurrentCounterIndex = tDaughterNode->GetNumber();

      TString tCurrentCentralNodePath = gCurrentNodePath;

      Int_t iNGaps(0);

      for(Int_t iDaughter = 0; iDaughter < tDaughterNode->GetNdaughters(); iDaughter++)
      {
        if(TString(tDaughterNode->GetDaughter(iDaughter)->GetName()).Contains("Gap"))
        {
          iNGaps++;
        }
      }

      if(0 == iNGaps%2)
      {
        tCurrentCentralNodePath += TString::Format("/tof_glass_%d",iNGaps/2);
      }
      else
      {
        tCurrentCentralNodePath += TString::Format("/Gap_%d",(iNGaps-1)/2);
      }

      UInt_t uCounterAddress = CbmTofAddress::GetUniqueAddress(giCurrentModuleIndex, giCurrentCounterIndex, 0, 0, giCurrentModuleType);

      TCounter* tCounter = new TCounter(giCurrentModuleType, giCurrentModuleIndex, giCurrentCounterIndex,
                                        tDaughterNode, new TGeoPhysicalNode(tCurrentCentralNodePath.Data()),
                                        new TGeoPhysicalNode(gCurrentBoxNodePath.Data()));

      gAllCounterMap.emplace(uCounterAddress, tCounter);

    }
    else
    {
      gCurrentNodePath += "/" + (TString)tDaughterNode->GetName();
    }

    // Expand nodes recursively
    if(tDaughterNode->GetNdaughters() > 0) { ExpandNode(tDaughterNode); }

    // Remove a node's name from the current node path upon completing a
    // recursion step
    gCurrentNodePath.ReplaceAll("/"+(TString)tDaughterNode->GetName(), "");
  }
}


Double_t gpb(Double_t* x, Double_t* p)
{
  // p[0] : location
  // p[1] : scale
  // p[2] : shape
  // p[3] : normalization

  Double_t dReturnValue = p[3]*0.5/p[2]*( TMath::Erf((x[0] - p[0] + p[2]/2.)/p[1]/TMath::Sqrt(2.))
                                         -TMath::Erf((x[0] - p[0] - p[2]/2.)/p[1]/TMath::Sqrt(2.)));

  return dReturnValue;
}


Double_t nh(Double_t* x, Double_t* p)
{
  // p[0] : location
  // p[1] : scale
  // p[2] : shape
  // p[3] : normalization
  // p[4] : parabola scaling
  // p[5] : parabola vertex (p[5], 1)

  Double_t dReturnValue = p[3]*0.5/p[2]*( TMath::Erf((x[0] - p[0] + p[2]/2.)/p[1]/TMath::Sqrt(2.))
                                         -TMath::Erf((x[0] - p[0] - p[2]/2.)/p[1]/TMath::Sqrt(2.)))
                                       *(1. + p[4]*TMath::Power(x[0] - p[5], 2));

  return dReturnValue;
}