/**
 * @file
 * @author Christian Simon <csimon@physi.uni-heidelberg.de>
 * @since 2017-06-19
 */

#include "CbmTofCounter.h"

#include "CbmTofDigiTbPar.h"
#include "CbmTofSignalBuffer.h"
#include "CbmTofPoint.h"
#include "CbmTofDigiExp.h"
#include "CbmTofHit.h"
#include "CbmMatch.h"
#include "CbmTofMemory.h"
#include "CbmDaqBuffer.h"
#include "CbmMCBuffer.h"
#include "CbmTofDef.h"

#include "FairLogger.h"

#include "TGeoManager.h"
#include "TGeoPhysicalNode.h"
#include "TRandom3.h"
#include "TMath.h"

#include <iterator>
#include <set>

#include <gsl/gsl_multimin.h>
#include <gsl/gsl_deriv.h>

// -----   initialization of static member variables   -----------------------
Bool_t CbmTofCounter::fbEnforceTwoSidedReadout = kFALSE;
Bool_t CbmTofCounter::fbIgnoreSignalInterference = kFALSE;
Bool_t CbmTofCounter::fbIgnoreChannelDeadTime = kFALSE;
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
CbmTofCounter::CbmTofCounter(Int_t iModuleType, Int_t iModuleIndex, Int_t iCounterIndex) :
  fiModuleType(iModuleType),
  fiModuleIndex(iModuleIndex),
  fiCounterIndex(iCounterIndex),
  fiCounterAddress(0),
  fCentralNodePath(""),
  fBoxNodePath(""),
  fCounterNode(NULL),
  fCentralPN(NULL),
  fBoxPN(NULL),
  fbPadReadout(kFALSE),
  fiNumberReadoutElectrodes(0),
  fdCellAcrossWidth(0.),
  fdCellAlongWidth(0.),
  fdGapAcrossWidth(0.),
  fdGapAlongWidth(0.),
  fbCellsAlongX(kFALSE),
  fbCellIndexingAlongAxis(kFALSE),
  fiNWalkBinsX(0),
  fdWalkBinWidth(0.),
  fdMinimumToT(0.),
  fdMaximumToT(0.),
  fdSignalTimeOffset(0.),
  fdPropagationVelocity(0.),
  fdDigitizationJitterSigma(0.),
  fdMaximumRuntimeOffset(0.),
  fdToTTargetExpectationValue(0.),
  fdToTTargetStandardDeviation(0.),
  fdSignalEdgeRuntimeOffsetsSigma(0.),
  fiMaximumClusterSize(0),
  fdChargeModus(0.),
  fdChargeScaling(0.),
  fdChargeDistance(0.),
  fdDiscriminationJitterSigma(0.),
  fdSignalTimeConstant(0.),
  fdSignalThreshold(0.),
  fdStripToTOffset(0.),
  fdAvalancheJitterSigma(0.),
  fdToTDistributionScaling(),
  fdToTDistributionOffset(),
  fdRuntimeCorrections(),
  fdSlewingCorrections(),
  fdRuntimeOffsets(),
  fdSignalEdgeRuntimeOffsets(),
  fdAcrossPositionResidualSigma(0.),
  fdAlongPositionResidualSigma(0.),
  fdTimeResidualSigma(0.),
  fdHitDistanceScale(0.),
  fdRelaxationTimeConstant(0.),
  fdChannelDeadTime(0.),
  fdCounterDarkRate(0.),
  fdCalibrationCoefficient(1.),
  fdWorkingCoefficient(1.),
  fdLastDarkPointTime(0.),
  fdMeanDarkPointInterval(0.),
  fdEventTimeOffset(0.),
  fbTwoSidedReadout(kTRUE),
  fbSendChargesToMemory(kFALSE),
  fbCalibrateDigis(kTRUE),
  fdSignalParameters {0., 0., 0., 0.},
  fFindIntersection(NULL),
  fFindSlope(NULL),
  fMinimizer(NULL),
  fRootFinder(NULL),
  fdLandauSignalMPV(0.),
  fdChargeSamplingLimit(100.),
  fSignalBuffer(NULL),
  fdLastLeadingEdgeTime(),
  fdLastTrailingEdgeTime(),
  fbToTWeights(kTRUE),
  fiDigiIndex(0),
  fDigiUnmergedLeadingEdgeToT(),
  fdHitTimeContributions(),
  fiFileIndex(0)
{
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
CbmTofCounter::~CbmTofCounter()
{
  if(fCentralPN)
  {
    delete fCentralPN;
  }

  if(fBoxPN)
  {
    delete fBoxPN;
  }

  if(fFindIntersection)
  {
    delete fFindIntersection;
  }

  if(fFindSlope)
  {
    delete fFindSlope;
  }

  if(fMinimizer)
  {
    delete fMinimizer;
  }

  if(fRootFinder)
  {
    delete fRootFinder;
  }

  if(fSignalBuffer)
  {
    delete fSignalBuffer;
  }
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
void CbmTofCounter::SetEventTimeOffset(Double_t dOffset)
{
  fdEventTimeOffset = dOffset;
  fdLastDarkPointTime += fdEventTimeOffset;
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
void CbmTofCounter::MasterToLocal(const Double_t* dMasterCoordinates, Double_t* dLocalCoordinates) const
{
  if(!fCentralPN)
  {
    LOG(FATAL)<<"Counter must be initialized before any coordinate transformation."<<FairLogger::endl;
  }

  fCentralPN->GetMatrix()->MasterToLocal(dMasterCoordinates, dLocalCoordinates);
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
void CbmTofCounter::LocalToMaster(const Double_t* dLocalCoordinates, Double_t* dMasterCoordinates) const
{
  if(!fCentralPN)
  {
    LOG(FATAL)<<"Counter must be initialized before any coordinate transformation."<<FairLogger::endl;
  }

  fCentralPN->GetMatrix()->LocalToMaster(dLocalCoordinates, dMasterCoordinates);
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
void CbmTofCounter::MasterToLocalBox(const Double_t* dMasterCoordinates, Double_t* dLocalCoordinates) const
{
  if(!fBoxPN)
  {
    LOG(FATAL)<<"Counter must be initialized before any coordinate transformation."<<FairLogger::endl;
  }

  fBoxPN->GetMatrix()->MasterToLocal(dMasterCoordinates, dLocalCoordinates);
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
void CbmTofCounter::LocalToMasterBox(const Double_t* dLocalCoordinates, Double_t* dMasterCoordinates) const
{
  if(!fBoxPN)
  {
    LOG(FATAL)<<"Counter must be initialized before any coordinate transformation."<<FairLogger::endl;
  }

  fBoxPN->GetMatrix()->LocalToMaster(dLocalCoordinates, dMasterCoordinates);
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
Bool_t CbmTofCounter::IsContainedInBox(const Double_t* dLocalCoordinates) const
{
  if(!fBoxPN)
  {
    LOG(FATAL)<<"Counter must be initialized before any coordinate transformation."<<FairLogger::endl;
  }

  return fBoxPN->GetVolume()->Contains(dLocalCoordinates);
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
Bool_t CbmTofCounter::Initialize(const CbmTofDigiTbPar* tTofDigiTbParSet)
{
  if(!tTofDigiTbParSet)
  {
    LOG(ERROR)<<"Pointer to parameter container 'CbmTofDigiTbPar' is invalid."<<FairLogger::endl;
    return kFALSE;
  }

  if(!fCounterNode)
  {
    LOG(ERROR)<<"Pointer to counter node invalid."<<FairLogger::endl;
    return kFALSE;
  }

  if(fCentralNodePath.IsNull() || !gGeoManager || !gGeoManager->CheckPath(fCentralNodePath.Data()))
  {
    LOG(ERROR)<<"Cannot get the central physical node of the counter."<<FairLogger::endl;
    return kFALSE;
  }

  fCentralPN = new TGeoPhysicalNode(fCentralNodePath.Data());

  if(fBoxNodePath.IsNull() || !gGeoManager->CheckPath(fBoxNodePath.Data()))
  {
    LOG(ERROR)<<"Cannot get the physical node of the counter box."<<FairLogger::endl;
    return kFALSE;
  }

  fBoxPN = new TGeoPhysicalNode(fBoxNodePath.Data());

  Int_t iNumberGlassPlates(0);
  Int_t iNumberGasGaps(0);
//  Double_t dGlassThickness(0.);
//  Double_t dGapSize(0.);

  Double_t dCellXWidth(0.);
  Double_t dCellYWidth(0.);
  Double_t dGapXWidth(0.);
  Double_t dGapYWidth(0.);

  Double_t dLocalCoordinates[3] = {0.,0.,0.};
  Double_t dGlobalCoordinates[3] = {0.,0.,0.};

  fCentralPN->GetMatrix()->LocalToMaster(dLocalCoordinates,dGlobalCoordinates);

  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))
    {
      if(0 == iNumberGlassPlates)
      {
        TGeoNode* tPlateNode = fCounterNode->GetDaughter(iComponent);

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

//        dGlassThickness = tPlateShape->GetAxisRange(3,dAxisLow,dAxisHigh);
      }

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

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

        dGapXWidth = tGapShape->GetAxisRange(1,dAxisLow,dAxisHigh);
        dGapYWidth = tGapShape->GetAxisRange(2,dAxisLow,dAxisHigh);
//        dGapSize = tGapShape->GetAxisRange(3,dAxisLow,dAxisHigh);

        fiNumberReadoutElectrodes = tGapNode->GetNdaughters();

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

        dCellXWidth = tFirstCellShape->GetAxisRange(1,dAxisLow,dAxisHigh);
        dCellYWidth = tFirstCellShape->GetAxisRange(2,dAxisLow,dAxisHigh);

        if(dCellXWidth == dGapXWidth)
        {
          fbCellsAlongX = kTRUE;

          fdCellAcrossWidth = dCellYWidth;
          fdCellAlongWidth = dCellXWidth;

          fdGapAcrossWidth = dGapYWidth;
          fdGapAlongWidth = dGapXWidth;
        }
        else
        {
          fdCellAcrossWidth = dCellXWidth;
          fdCellAlongWidth = dCellYWidth;

          fdGapAcrossWidth = dGapXWidth;
          fdGapAlongWidth = dGapYWidth;
        }

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

        TGeoNode* tLastCellNode = tGapNode->GetDaughter(fiNumberReadoutElectrodes-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
          {
            LOG(ERROR)<<"Incomplete understanding of geometry. Cannot determine direction of cell indexing."<<FairLogger::endl;
            return kFALSE;
          }
        }
        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
          {
            LOG(ERROR)<<"Incomplete understanding of geometry. Cannot determine direction of cell indexing."<<FairLogger::endl;
            return kFALSE;
          }
        }
      }

      iNumberGasGaps++;
    }
  }


  if(   !tTofDigiTbParSet->GetNumberReadoutElectrodes(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetCellAcrossWidth(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetCellAlongWidth(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetGapAcrossWidth(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetGapAlongWidth(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetCellsAlongX(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetCellIndexingAlongAxis(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetNWalkBinsX(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetMinimumToT(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetMaximumToT(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetSignalTimeOffset(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetPropagationVelocity(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetDigitizationJitterSigma(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetMaximumRuntimeOffset(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetToTTargetExpectationValue(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetToTTargetStandardDeviation(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetSignalEdgeRuntimeOffsetsSigma(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetMaximumClusterSize(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetChargeModus(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetChargeScaling(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetChargeDistance(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetDiscriminationJitterSigma(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetSignalTimeConstant(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetSignalThreshold(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetStripToTOffset(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetAvalancheJitterSigma(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetToTDistributionScaling(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetToTDistributionOffset(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetRuntimeCorrections(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetSlewingCorrections(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetRuntimeOffsets(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetSignalEdgeRuntimeOffsets(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetAcrossPositionResidualSigma(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetAlongPositionResidualSigma(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetTimeResidualSigma(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetHitDistanceScale(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetRelaxationTimeConstant(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetChannelDeadTime(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetCounterDarkRate(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetTwoSidedReadout(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetCalibrationCoefficient(fiModuleType, fiModuleIndex, fiCounterIndex)
     || !tTofDigiTbParSet->GetWorkingCoefficient(fiModuleType, fiModuleIndex, fiCounterIndex) )
  {
    LOG(ERROR)<<"Incomplete parameter set 'CbmTofDigiTbPar'."<<FairLogger::endl;
    return kFALSE;
  }

  if(   fiNumberReadoutElectrodes != *tTofDigiTbParSet->GetNumberReadoutElectrodes(fiModuleType, fiModuleIndex, fiCounterIndex)
     || fdCellAcrossWidth != *tTofDigiTbParSet->GetCellAcrossWidth(fiModuleType, fiModuleIndex, fiCounterIndex)
     || fdCellAlongWidth != *tTofDigiTbParSet->GetCellAlongWidth(fiModuleType, fiModuleIndex, fiCounterIndex)
     || fdGapAcrossWidth != *tTofDigiTbParSet->GetGapAcrossWidth(fiModuleType, fiModuleIndex, fiCounterIndex)
     || fdGapAlongWidth != *tTofDigiTbParSet->GetGapAlongWidth(fiModuleType, fiModuleIndex, fiCounterIndex)
     || fbCellsAlongX != *tTofDigiTbParSet->GetCellsAlongX(fiModuleType, fiModuleIndex, fiCounterIndex)
     || fbCellIndexingAlongAxis != *tTofDigiTbParSet->GetCellIndexingAlongAxis(fiModuleType, fiModuleIndex, fiCounterIndex) )
  {
    LOG(ERROR)<<"Counter geometry seems to have changed compared to modelling time."<<FairLogger::endl;
    return kFALSE;
  }

  // TODO: Should we run out of memory at some point we could use pointer member variables instead of copy-constructing everything
  fiNWalkBinsX = *tTofDigiTbParSet->GetNWalkBinsX(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdMinimumToT = *tTofDigiTbParSet->GetMinimumToT(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdMaximumToT = *tTofDigiTbParSet->GetMaximumToT(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdSignalTimeOffset = *tTofDigiTbParSet->GetSignalTimeOffset(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdPropagationVelocity = *tTofDigiTbParSet->GetPropagationVelocity(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdDigitizationJitterSigma = *tTofDigiTbParSet->GetDigitizationJitterSigma(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdMaximumRuntimeOffset = *tTofDigiTbParSet->GetMaximumRuntimeOffset(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdToTTargetExpectationValue = *tTofDigiTbParSet->GetToTTargetExpectationValue(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdToTTargetStandardDeviation = *tTofDigiTbParSet->GetToTTargetStandardDeviation(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdSignalEdgeRuntimeOffsetsSigma = *tTofDigiTbParSet->GetSignalEdgeRuntimeOffsetsSigma(fiModuleType, fiModuleIndex, fiCounterIndex);
  fiMaximumClusterSize = *tTofDigiTbParSet->GetMaximumClusterSize(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdChargeModus = *tTofDigiTbParSet->GetChargeModus(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdChargeScaling = *tTofDigiTbParSet->GetChargeScaling(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdChargeDistance = *tTofDigiTbParSet->GetChargeDistance(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdDiscriminationJitterSigma = *tTofDigiTbParSet->GetDiscriminationJitterSigma(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdSignalTimeConstant = *tTofDigiTbParSet->GetSignalTimeConstant(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdSignalThreshold = *tTofDigiTbParSet->GetSignalThreshold(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdStripToTOffset = *tTofDigiTbParSet->GetStripToTOffset(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdAvalancheJitterSigma = *tTofDigiTbParSet->GetAvalancheJitterSigma(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdToTDistributionScaling = *tTofDigiTbParSet->GetToTDistributionScaling(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdToTDistributionOffset = *tTofDigiTbParSet->GetToTDistributionOffset(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdRuntimeCorrections = *tTofDigiTbParSet->GetRuntimeCorrections(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdSlewingCorrections = *tTofDigiTbParSet->GetSlewingCorrections(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdRuntimeOffsets = *tTofDigiTbParSet->GetRuntimeOffsets(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdSignalEdgeRuntimeOffsets = *tTofDigiTbParSet->GetSignalEdgeRuntimeOffsets(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdAcrossPositionResidualSigma = *tTofDigiTbParSet->GetAcrossPositionResidualSigma(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdAlongPositionResidualSigma = *tTofDigiTbParSet->GetAlongPositionResidualSigma(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdTimeResidualSigma = *tTofDigiTbParSet->GetTimeResidualSigma(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdHitDistanceScale = *tTofDigiTbParSet->GetHitDistanceScale(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdRelaxationTimeConstant = *tTofDigiTbParSet->GetRelaxationTimeConstant(fiModuleType, fiModuleIndex, fiCounterIndex)*1e9; // convert to ns
  fdChannelDeadTime = *tTofDigiTbParSet->GetChannelDeadTime(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdCounterDarkRate = *tTofDigiTbParSet->GetCounterDarkRate(fiModuleType, fiModuleIndex, fiCounterIndex);
  fbTwoSidedReadout = *tTofDigiTbParSet->GetTwoSidedReadout(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdCalibrationCoefficient = *tTofDigiTbParSet->GetCalibrationCoefficient(fiModuleType, fiModuleIndex, fiCounterIndex);
  fdWorkingCoefficient = *tTofDigiTbParSet->GetWorkingCoefficient(fiModuleType, fiModuleIndex, fiCounterIndex);

  fdWalkBinWidth = (fdMaximumToT - fdMinimumToT)/static_cast<Double_t>(fiNWalkBinsX);

  // fdSignalParameters[2] holds the strip charge
  fdSignalParameters[0] = fdSignalTimeOffset;
  fdSignalParameters[1] = fdSignalTimeConstant;
  fdSignalParameters[2] = 0.;
  fdSignalParameters[3] = fdSignalThreshold;

  gsl_function find_maximum;
  find_maximum.function = &NegativeLandau;

  fFindIntersection = new gsl_function();
  fFindIntersection->function = &LandauThreshold;

  fFindSlope = new gsl_function();
  fFindSlope->function = &LandauSignal;

  const gsl_min_fminimizer_type* tMinimizerType = gsl_min_fminimizer_brent;
  const gsl_root_fsolver_type* tRootFinderType = gsl_root_fsolver_brent;

  // Calculate the maximum total induced charge sampled from the induced charge
  // distribution to make the cluster size distribution smooth out at a desired size.

  Double_t dMaximumChargeDensityLowerAcrossLimit = (fiMaximumClusterSize/2)*fdCellAcrossWidth;
  Double_t dMaximumChargeDensityUpperAcrossLimit = dMaximumChargeDensityLowerAcrossLimit + fdCellAcrossWidth;
  Double_t dMaximumChargeDensityLowerAlongLimit  = -fdCellAlongWidth/2.;
  Double_t dMaximumChargeDensityUpperAlongLimit  =  fdCellAlongWidth/2.;


  fMinimizer = gsl_min_fminimizer_alloc(tMinimizerType);

  find_maximum.params = fdSignalParameters;

  Int_t iIteration(0);
  Int_t iMinimizationStatus(0);

  Double_t dInitialMaximumEstimate = fdSignalParameters[0];
  Double_t dInitialMaximumLowerLimit = fdSignalParameters[0] - 10.;
  Double_t dInitialMaximumUpperLimit = fdSignalParameters[0] + 10.;

  gsl_min_fminimizer_set(fMinimizer, &find_maximum, dInitialMaximumEstimate, dInitialMaximumLowerLimit, dInitialMaximumUpperLimit);

  do
  {
    iIteration++;
    gsl_min_fminimizer_iterate(fMinimizer);

    iMinimizationStatus = gsl_min_test_interval(fMinimizer->x_lower, fMinimizer->x_upper, 0.0001, 0.0);
  }
  while( GSL_CONTINUE == iMinimizationStatus && 100 > iIteration );

  fdLandauSignalMPV = fMinimizer->x_minimum;

  fdChargeSamplingLimit = 1./(-fMinimizer->f_minimum); // minimize the negative function to maximize it


  fdChargeSamplingLimit *= 2.*fdSignalThreshold*ChargeLimit(fdChargeDistance,dMaximumChargeDensityLowerAlongLimit,dMaximumChargeDensityUpperAlongLimit,
                                                                             dMaximumChargeDensityLowerAcrossLimit,dMaximumChargeDensityUpperAcrossLimit);
  fdChargeSamplingLimit /= fdCalibrationCoefficient;


  fRootFinder = gsl_root_fsolver_alloc(tRootFinderType);

  if(fdCounterDarkRate)
  {
    fdMeanDarkPointInterval = 1e9/(fdCounterDarkRate*fdGapAcrossWidth*fdGapAlongWidth); // convert to ns
  }

  fdLastDarkPointTime += gRandom->Exp(fdMeanDarkPointInterval);

  if(fbEnforceTwoSidedReadout)
  {
    fbTwoSidedReadout = kTRUE;
  }

  // create a signal buffer for all readout channels
  fSignalBuffer = new CbmTofSignalBuffer((static_cast<Int_t>(!fbPadReadout) + 1)*fiNumberReadoutElectrodes);
  fSignalBuffer->SetIgnoreInterference(fbIgnoreSignalInterference);

  fdLastLeadingEdgeTime.resize((static_cast<Int_t>(!fbPadReadout) + 1)*fiNumberReadoutElectrodes, -fdChannelDeadTime);
  fdLastTrailingEdgeTime.resize((static_cast<Int_t>(!fbPadReadout) + 1)*fiNumberReadoutElectrodes, -fdChannelDeadTime);

  fDigiUnmergedLeadingEdgeToT.resize((static_cast<Int_t>(!fbPadReadout) + 1)*fiNumberReadoutElectrodes);

  fdHitTimeContributions.reserve(fiNumberReadoutElectrodes);

  return kTRUE;
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
void CbmTofCounter::ProcessPoint(const CbmTofPoint* tTofPoint)
{
  // TODO: Add support for pad readout.
  //       Requires a new geometrical description of pad MRPCs in CbmRoot.
  //       Currently, single pad rows are considered a counter and several
  //       such counters are combined in software to represent a real pad
  //       MRPC. This digitization approach here, however, requires a coherent
  //       readout plane in which the sampled induced charge can be distributed.

  // IMPORTANT NOTE: To simulate degradation effects in MRPCs due to continuous
  // irradation with heavy-ion reaction products the digitization process cannot
  // be split into many small jobs whose individual results would be merged in
  // the end. The need to keep track of all previously generated charges in the
  // electric field (implemented in 'CbmTofMemory') necessitates simulating
  // (millions of) events in ONE single job.
  // The file source class connects all MC input files of THIS job to a single 'TChain'.
  // In the reconstruction, 'CbmMCDataManager' will allow access to this chain
  // of files by index 'CbmTofCounter::fiFileIndex'.
  // The event ID gives the entry in the chain, NOT in the single MC files.
  // The point ID gives the position of the MC point in its 'TClonesArray' container.
  // The link ID distinguishes source points from dark points and primary signals
  // from afterpulses and reflections (cf. 'CbmTofSignalOrigin' in 'CbmTofDef.h').

  Int_t iMCFileID = fiFileIndex;
  Int_t iMCEventID = tTofPoint->GetEventID();
  Int_t iMCPointID = tTofPoint->GetUniqueID();
  Int_t iMCLinkID = tof::signal_SourcePrimary;

  // Dark rate points do not belong to any event and cannot be read from the
  // input 'TChain'.
  if(-1 == iMCEventID)
  {
    iMCFileID = -1;
    iMCLinkID = tof::signal_DarkPrimary;
  }


  Double_t dLocalHitCoordinates[3] = {0.,0.,0.};
  Double_t dGlobalHitCoordinates[3] = {0.,0.,0.};

  dGlobalHitCoordinates[0] = tTofPoint->GetX();
  dGlobalHitCoordinates[1] = tTofPoint->GetY();
  dGlobalHitCoordinates[2] = tTofPoint->GetZ();

  fCentralPN->GetMatrix()->MasterToLocal(dGlobalHitCoordinates,dLocalHitCoordinates);

  Double_t dMCAcrossStripPosition(0.);
  Double_t dMCAlongStripPosition(0.);

  if(fbCellsAlongX)
  {
    dMCAcrossStripPosition = dLocalHitCoordinates[1];
    dMCAlongStripPosition = dLocalHitCoordinates[0];
  }
  else
  {
    dMCAcrossStripPosition = dLocalHitCoordinates[0];
    dMCAlongStripPosition = dLocalHitCoordinates[1];
  }

  Double_t dMCHitTime = tTofPoint->GetTime();

  Double_t dInducedCharge(0.);

  do
  {
    dInducedCharge = gRandom->Landau(fdChargeModus,fdChargeScaling);
  }
  while(0 > dInducedCharge || fdChargeSamplingLimit < dInducedCharge);

  if(fbSendChargesToMemory && fdRelaxationTimeConstant)
  {
    dInducedCharge *= CbmTofMemory::Instance()->GetMemoryCoefficient(fiCounterAddress, dInducedCharge, dMCHitTime, dMCAcrossStripPosition,
                                                                     dMCAlongStripPosition, fdHitDistanceScale, fdRelaxationTimeConstant,
                                                                     fdChargeSamplingLimit);
  }

  dInducedCharge *= fdWorkingCoefficient;

/*
  Int_t iHitStrip(-1);
  Double_t dHitStripCenterAcross(0.);

  if(fbCellIndexingAlongAxis)
  {
    iHitStrip = static_cast<Int_t>(dMCAcrossStripPosition/fdCellAcrossWidth + static_cast<Double_t>(fiNumberReadoutElectrodes)/2.);
    dHitStripCenterAcross = (static_cast<Double_t>(iHitStrip) - static_cast<Double_t>(fiNumberReadoutElectrodes)/2. + 0.5)*fdCellAcrossWidth;
  }
  else
  {
    iHitStrip = -static_cast<Int_t>(dMCAcrossStripPosition/fdCellAcrossWidth - static_cast<Double_t>(fiNumberReadoutElectrodes)/2.);
    dHitStripCenterAcross = -(static_cast<Double_t>(iHitStrip) - static_cast<Double_t>(fiNumberReadoutElectrodes)/2. + 0.5)*fdCellAcrossWidth;
  }
*/

  Double_t dAvalancheJitter = gRandom->Gaus(0.0, fdAvalancheJitterSigma);

  Double_t dPropagationTimeAlongCell[2] = {0.,0.};

  // negative cell end
  dPropagationTimeAlongCell[0] = (fdCellAlongWidth/2. + dMCAlongStripPosition)/fdPropagationVelocity;
  // positive cell end
  dPropagationTimeAlongCell[1] = (fdCellAlongWidth/2. - dMCAlongStripPosition)/fdPropagationVelocity;


  Double_t dChargeDensityLowerAcrossLimit = -fdGapAcrossWidth/2. - dMCAcrossStripPosition;
  Double_t dChargeDensityUpperAcrossLimit =  fdGapAcrossWidth/2. - dMCAcrossStripPosition;
  Double_t dChargeDensityLowerAlongLimit  = -fdGapAlongWidth/2.  - dMCAlongStripPosition;
  Double_t dChargeDensityUpperAlongLimit  =  fdGapAlongWidth/2.  - dMCAlongStripPosition;

  Double_t dChargeDensityStripLowerAcrossLimit(0.);
  Double_t dChargeDensityStripUpperAcrossLimit(0.);


  for(Int_t iStripIndex = 0; iStripIndex < fiNumberReadoutElectrodes; iStripIndex++)
  {
    if(fbCellIndexingAlongAxis)
    {
      dChargeDensityStripLowerAcrossLimit = dChargeDensityLowerAcrossLimit + iStripIndex*fdCellAcrossWidth;
      dChargeDensityStripUpperAcrossLimit = dChargeDensityLowerAcrossLimit + (iStripIndex+1)*fdCellAcrossWidth;
    }
    else
    {
      dChargeDensityStripLowerAcrossLimit = dChargeDensityUpperAcrossLimit - iStripIndex*fdCellAcrossWidth;
      dChargeDensityStripUpperAcrossLimit = dChargeDensityUpperAcrossLimit - (iStripIndex+1)*fdCellAcrossWidth;
    }

    Double_t dStripCharge = DistributeCharge(dInducedCharge,fdChargeDistance,dChargeDensityLowerAlongLimit,dChargeDensityUpperAlongLimit,dChargeDensityStripLowerAcrossLimit,dChargeDensityStripUpperAcrossLimit);


    for( Int_t iStripSide = 0; iStripSide < static_cast<Int_t>(fbTwoSidedReadout) + 1; iStripSide++ )
    {
      fdSignalParameters[2] = 0.5*dStripCharge;
      fFindIntersection->params = fdSignalParameters;
      fFindSlope->params = fdSignalParameters;

      Double_t dSignalPeak = 0.5*dStripCharge*TMath::Landau(fdLandauSignalMPV,fdSignalTimeOffset,fdSignalTimeConstant,kTRUE);

      if(fdSignalThreshold < dSignalPeak)
      {
        Double_t dStripLeadingEdge(0.);
        Double_t dStripTrailingEdge(0.);

        Double_t dStripLeadingEdgeSlope(0.);
        Double_t dStripTrailingEdgeSlope(0.);


        Int_t iIteration(0);
        Int_t iRootFindingStatus(0);

        // leading edge
        Double_t dInitialRootLowerLimit = -1000.;
        Double_t dInitialRootUpperLimit = fdLandauSignalMPV;

        gsl_root_fsolver_set(fRootFinder, fFindIntersection, dInitialRootLowerLimit, dInitialRootUpperLimit);

        do
        {
          iIteration++;
          gsl_root_fsolver_iterate(fRootFinder);

          iRootFindingStatus = gsl_root_test_interval(fRootFinder->x_lower, fRootFinder->x_upper, 0, 0.00001); // 0.001
        }
        while( GSL_CONTINUE == iRootFindingStatus && 100 > iIteration );

        iIteration = 0;

        dStripLeadingEdge = fRootFinder->root;

        Double_t dStripLeadingEdgeSlopeError(0.);
        gsl_deriv_central(fFindSlope, dStripLeadingEdge, 1.e-8, &dStripLeadingEdgeSlope, &dStripLeadingEdgeSlopeError);

        // trailing edge
        dInitialRootLowerLimit = fdLandauSignalMPV;
        dInitialRootUpperLimit = 1000.;

        gsl_root_fsolver_set(fRootFinder, fFindIntersection, dInitialRootLowerLimit, dInitialRootUpperLimit);

        do
        {
          iIteration++;
          gsl_root_fsolver_iterate(fRootFinder);

          iRootFindingStatus = gsl_root_test_interval(fRootFinder->x_lower, fRootFinder->x_upper, 0, 0.00001); // 0.001
        }
        while( GSL_CONTINUE == iRootFindingStatus && 100 > iIteration );

        dStripTrailingEdge = fRootFinder->root;

        Double_t dStripTrailingEdgeSlopeError(0.);
        gsl_deriv_central(fFindSlope, dStripTrailingEdge, 1.e-8, &dStripTrailingEdgeSlope, &dStripTrailingEdgeSlopeError);


        Double_t dNumericalStripLeadingEdge = dStripLeadingEdge;
        Double_t dNumericalStripTrailingEdge = dStripTrailingEdge;

        Double_t dLeadingEdgeDiscriminationJitter = gRandom->Gaus(0.0, fdDiscriminationJitterSigma/TMath::Abs(dStripLeadingEdgeSlope));
        Double_t dTrailingEdgeDiscriminationJitter = gRandom->Gaus(0.0, fdDiscriminationJitterSigma/TMath::Abs(dStripTrailingEdgeSlope));

        // make sure that the trailing edge is not shifted to before the leading edge discrimination point in time
        if( dNumericalStripLeadingEdge + dLeadingEdgeDiscriminationJitter > dNumericalStripTrailingEdge + dTrailingEdgeDiscriminationJitter )
        {
          dNumericalStripTrailingEdge = dNumericalStripLeadingEdge + dLeadingEdgeDiscriminationJitter;
          dTrailingEdgeDiscriminationJitter = TMath::Abs(dTrailingEdgeDiscriminationJitter);
        }


        // Time and ToT information might be merged for two consecutive events in a single readout channel.
        // Applying the (optional) calibration prior to the merging process would guarantee that the leading edge time
        // information of the result had its slewing corrected for. The (extended) new ToT value, however, were
        // too big for the equalized walk effect.
        // If the calibration were applied after the merging process, the walk effect would be overcorrected
        // for (due to the too big merged ToT) resulting in biased leading edge information.
        // Downscaling ToT spectra and walk correcting leading-edge time spectra before possible signal
        // merging in the RO buffer might cause signals that would have overlapped in the cable without
        // corrections not to overlap anymore when corrected. For this reason, the decision is taken to
        // only calibrate signal information after signal extraction from the RO buffer.
        // Different signal runtime offsets between leading and trailing edge channels throughout the counter
        // can cause the ToT spectra of certain channels to start with negative values. Such negative ToT
        // values cannot be dealt with in the readout buffer. For this reason, the data objects in the buffer
        // will be given a(n uncorrected) ToT with 'fdSignalEdgeRuntimeOffsets' subtracted.
        Double_t dDiscriminatedTime =  dMCHitTime + dAvalancheJitter + dPropagationTimeAlongCell[iStripSide] + dNumericalStripLeadingEdge + dLeadingEdgeDiscriminationJitter
                                     + fdRuntimeOffsets.at(2*iStripIndex+iStripSide);
        Double_t dDiscriminatedToT = (dNumericalStripTrailingEdge + dTrailingEdgeDiscriminationJitter)
                                    -(dNumericalStripLeadingEdge + dLeadingEdgeDiscriminationJitter)
                                    + fdStripToTOffset;


        CbmTofDigiExp* tDiscriminationDigi = new CbmTofDigiExp(fiModuleIndex, fiCounterIndex, iStripIndex, dDiscriminatedTime, dDiscriminatedToT, iStripSide, fiModuleType);
        CbmMatch* tDigiMatch = new CbmMatch();
        tDiscriminationDigi->SetMatch(tDigiMatch);
        // MC point contributions to digi objects are weighted with uncorrected and uncalibrated ToT information
        // because ToT calibration can shift a spectrum partially into the range of negative numbers.
        // Weights are not normalized. This can be done in the reconstruction step if needed.
        tDigiMatch->AddLink(dDiscriminatedToT, iMCPointID, iMCEventID, iMCFileID);

        CbmLink& tLink = tDigiMatch->GetLink(0);
        tLink.SetUniqueID(iMCLinkID);

        auto & DigiToTMap = fDigiUnmergedLeadingEdgeToT.at(2*iStripIndex + iStripSide);
        DigiToTMap.emplace_hint(DigiToTMap.end(), fiDigiIndex, dDiscriminatedToT);
        tDiscriminationDigi->SetUniqueID(fiDigiIndex++);

        fSignalBuffer->Fill(2*iStripIndex + iStripSide, tDiscriminationDigi);
      }
    }
  }
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
Int_t CbmTofCounter::ReadOutChannelBuffers(Double_t dReadoutTime)
{
  Int_t iNDigis(0);
  std::vector<CbmTofDigiExp*> tDigiVector;
  tDigiVector.reserve(10);

  for(Int_t iStripIndex = 0; iStripIndex < fiNumberReadoutElectrodes; iStripIndex++)
  {
    for( Int_t iStripSide = 0; iStripSide < static_cast<Int_t>(fbTwoSidedReadout) + 1; iStripSide++ )
    {
      iNDigis += fSignalBuffer->ReadOutData(2*iStripIndex + iStripSide, dReadoutTime, tDigiVector);

      auto & DigiToTMap = fDigiUnmergedLeadingEdgeToT.at(2*iStripIndex + iStripSide);
      auto & MergedDigiIndices = fSignalBuffer->GetMergedDigiIndices(2*iStripIndex + iStripSide);

      for(auto const & MergedIndex : MergedDigiIndices)
      {
        auto itTemp = DigiToTMap.find(MergedIndex);

        if(itTemp != DigiToTMap.end())
        {
          DigiToTMap.erase(itTemp);
        }
      }

      MergedDigiIndices.clear();

      for(auto & itDigi : tDigiVector)
      {
        CbmTofDigiExp* tCurrentDigi = itDigi;
        Double_t dDigiLeadingEdgeTime = tCurrentDigi->GetTime();
        Double_t dDigiToT = tCurrentDigi->GetTot();
        Double_t dDigiTrailingEdgeTime = dDigiLeadingEdgeTime + dDigiToT;


        // To check if a signal arrives at a currently dead TDC channel both its
        // discriminated leading and its discriminated trailing edge time are
        // checked against the leading and trailing edge times of the last
        // digitized signal, kept track of in 'fdLastLeadingEdgeTime' and in
        // 'fdLastTrailingEdgeTime'. This is done prior to adding digitization
        // jitter as the arrival time of the signal at the delay line matters
        // here and not how the TDC discretizes the information in the end.
        if(fbIgnoreChannelDeadTime ||
           (dDigiLeadingEdgeTime > fdLastLeadingEdgeTime.at(2*iStripIndex + iStripSide) + fdChannelDeadTime
            && dDigiTrailingEdgeTime > fdLastTrailingEdgeTime.at(2*iStripIndex + iStripSide) + fdChannelDeadTime) )
        {
          fdLastLeadingEdgeTime.at(2*iStripIndex + iStripSide) = dDigiLeadingEdgeTime;
          fdLastTrailingEdgeTime.at(2*iStripIndex + iStripSide) = dDigiTrailingEdgeTime;

          Double_t dLeadingEdgeDigitizationJitter = gRandom->Gaus(0.0, fdDigitizationJitterSigma);
          Double_t dTrailingEdgeDigitizationJitter = gRandom->Gaus(0.0, fdDigitizationJitterSigma);

          dDigiLeadingEdgeTime += dLeadingEdgeDigitizationJitter;
          dDigiToT += dTrailingEdgeDigitizationJitter - dLeadingEdgeDigitizationJitter + fdSignalEdgeRuntimeOffsets.at(2*iStripIndex + iStripSide);

          if(fbCalibrateDigis)
          {
            dDigiToT = fdToTDistributionScaling.at(2*iStripIndex + iStripSide)*dDigiToT + fdToTDistributionOffset.at(2*iStripIndex + iStripSide);

            dDigiLeadingEdgeTime -= fdRuntimeCorrections.at(2*iStripIndex + iStripSide);

            Double_t dCalibDigiToT = DigiToTMap.at(tCurrentDigi->GetUniqueID());
            dCalibDigiToT += dTrailingEdgeDigitizationJitter - dLeadingEdgeDigitizationJitter + fdSignalEdgeRuntimeOffsets.at(2*iStripIndex + iStripSide);
            dCalibDigiToT = fdToTDistributionScaling.at(2*iStripIndex + iStripSide)*dCalibDigiToT + fdToTDistributionOffset.at(2*iStripIndex + iStripSide);
            DigiToTMap.erase(tCurrentDigi->GetUniqueID());


            // linear interpolation according to TH1::Interpolate
            if( dCalibDigiToT <= fdMinimumToT + 0.5*fdWalkBinWidth )
            {
              dDigiLeadingEdgeTime -= (fdSlewingCorrections.at(2*iStripIndex + iStripSide)).at(1);
            }
            else if( dCalibDigiToT >= fdMaximumToT - 0.5*fdWalkBinWidth )
            {
              dDigiLeadingEdgeTime -= (fdSlewingCorrections.at(2*iStripIndex + iStripSide)).at(fiNWalkBinsX);
            }
            else
            {
              Int_t iWalkBin = 1 + static_cast<Int_t>((dCalibDigiToT - fdMinimumToT)/fdWalkBinWidth );
              Double_t dWalkBinCenter = fdMinimumToT + (static_cast<Double_t>(iWalkBin) - 0.5)*fdWalkBinWidth;

              Double_t x0,x1,y0,y1;

              if( dCalibDigiToT <= dWalkBinCenter )
              {
                y0 = (fdSlewingCorrections.at(2*iStripIndex + iStripSide)).at(iWalkBin-1);
                x0 = dWalkBinCenter - fdWalkBinWidth;
                y1 = (fdSlewingCorrections.at(2*iStripIndex + iStripSide)).at(iWalkBin);
                x1 = dWalkBinCenter;
              }
              else
              {
                y0 = (fdSlewingCorrections.at(2*iStripIndex + iStripSide)).at(iWalkBin);
                x0 = dWalkBinCenter;
                y1 = (fdSlewingCorrections.at(2*iStripIndex + iStripSide)).at(iWalkBin+1);
                x1 = dWalkBinCenter + fdWalkBinWidth;
              }

              dDigiLeadingEdgeTime -= y0 + (dCalibDigiToT - x0)*((y1 - y0)/(x1 - x0));
            }
          }

          tCurrentDigi->SetTime(dDigiLeadingEdgeTime);
          tCurrentDigi->SetTot(dDigiToT);
          tCurrentDigi->SetUniqueID(0);

          CbmDaqBuffer::Instance()->InsertData(tCurrentDigi);
        }
        else
        {
          DigiToTMap.erase(tCurrentDigi->GetUniqueID());
          delete tCurrentDigi;
        }
      }

      tDigiVector.clear();
    }
  }

  // number of signals extracted from the readout buffer (not the number of actually digitized signals -> dead time!)
  return iNDigis;
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
void CbmTofCounter::BuildClusterMC(std::vector<std::pair<CbmTofDigiExp*, Int_t>>& tDigiVector, CbmTofHit* tHit, Bool_t bCompatibilityMode)
{
  // Sort the digi vector w.r.t. increasing channel indices.
  sort(tDigiVector.begin(), tDigiVector.end(), [] (const std::pair<CbmTofDigiExp*, Int_t>& Digi1, const std::pair<CbmTofDigiExp*, Int_t>& Digi2)
        { return 2*(Digi1.first->GetChannel()) + Digi1.first->GetSide() < 2*(Digi2.first->GetChannel()) + Digi2.first->GetSide();}
      );

  Double_t dClusterPositionAcrossStrips(0.);
  Double_t dClusterPositionAlongStrips(0.);
  Double_t dClusterTime(0.);
  Double_t dClusterTimeError(0.);
  Double_t dClusterSummedToT(0.);
  Int_t iClusterSize(0);

  Double_t dBottomLeadingEdge(0.);
  Double_t dTopLeadingEdge(0.);
  Double_t dBottomToT(0.);
  Double_t dTopToT(0.);
  Int_t iLastChannel(-1);
  Int_t iLastSide(-1);

  fdHitTimeContributions.clear();

  // Assumption: Each readout channel side has at most sent one "good" digi (fUniqueID % tof::kiPrimarySignalDivisor == 0) to DAQ!

  for(auto itDigi = tDigiVector.begin(); itDigi != tDigiVector.end(); ++itDigi)
  {
    CbmTofDigiExp* tCurrentDigi = itDigi->first;
    Int_t iChannel = tCurrentDigi->GetChannel();
    Int_t iSide = tCurrentDigi->GetSide();

    // bottom readout channel
    if(0 == iSide)
    {
      if(iLastSide == iSide)
      {
        // Information from the top readout channel is missing. Mark the PREVIOUS bottom information unused.
        std::prev(itDigi)->first = NULL;
      }

      dBottomLeadingEdge = tCurrentDigi->GetTime();
      dBottomToT = tCurrentDigi->GetTot();

      iLastChannel = iChannel;
      iLastSide = iSide;
    }
    // top readout channel
    else
    {
      if(iLastChannel == iChannel)
      {
        dTopLeadingEdge = tCurrentDigi->GetTime();
        dTopToT = tCurrentDigi->GetTot();

        Double_t dStripPositionAcrossStrips(0.);

        if(fbCellIndexingAlongAxis)
        {
          dStripPositionAcrossStrips = (static_cast<Double_t>(iChannel) - static_cast<Double_t>(fiNumberReadoutElectrodes)/2. + 0.5)*fdCellAcrossWidth;
        }
        else
        {
          dStripPositionAcrossStrips = -(static_cast<Double_t>(iChannel) - static_cast<Double_t>(fiNumberReadoutElectrodes)/2. + 0.5)*fdCellAcrossWidth;
        }

        if(fbToTWeights)
        {
          dClusterPositionAlongStrips += (dBottomToT + dTopToT)*(dBottomLeadingEdge - dTopLeadingEdge);
          dClusterTime += (dBottomToT + dTopToT)*(dBottomLeadingEdge + dTopLeadingEdge);
          dClusterPositionAcrossStrips += (dBottomToT + dTopToT)*dStripPositionAcrossStrips;
          dClusterSummedToT += (dBottomToT + dTopToT);
        }
        else
        {
          dClusterPositionAlongStrips += (dBottomLeadingEdge - dTopLeadingEdge);
          dClusterTime += (dBottomLeadingEdge + dTopLeadingEdge);
          dClusterPositionAcrossStrips += dStripPositionAcrossStrips;
          dClusterSummedToT += 1.;
        }

        fdHitTimeContributions.push_back(0.5*(dBottomLeadingEdge + dTopLeadingEdge));

        iClusterSize++;
      }
      else
      {
        // Information from the bottom readout channel is missing. Mark the CURRENT top information unused.
        itDigi->first = NULL;

        // If previous bottom readout channel information (from a different channel index) exists, mark it unused as well.
        if(iLastChannel != -1 && 0 == iLastSide)
        {
          std::prev(itDigi)->first = NULL;
        }
      }

      iLastSide = iSide;
    }
  }

  // Catch the case that a single bottom digi is the last (or the only) one in the buffer and - if so - mark it unused.
  if(0 == iLastSide)
  {
    (--tDigiVector.end())->first = NULL;
  }


  if(iClusterSize)
  {
    dClusterPositionAcrossStrips /= dClusterSummedToT;
    dClusterPositionAlongStrips *= 0.5*fdPropagationVelocity/dClusterSummedToT;
    dClusterTime *= 0.5/dClusterSummedToT;

    if(1 < iClusterSize)
    {
      for(auto & dStripTime : fdHitTimeContributions)
      {
        dClusterTimeError += TMath::Power(dStripTime - dClusterTime, 2);
      }

      dClusterTimeError = TMath::Sqrt(dClusterTimeError/(iClusterSize - 1));
    }

    // remove the fixed delay of the reconstructed time w.r.t. the MC time due to the finite signal propagation time on the readout cells
    dClusterTime -= 0.5*fdCellAlongWidth/fdPropagationVelocity;

    Int_t iClusterCentralStrip(-1);

    if(fbCellIndexingAlongAxis)
    {
      iClusterCentralStrip = static_cast<Int_t>(dClusterPositionAcrossStrips/fdCellAcrossWidth + static_cast<Double_t>(fiNumberReadoutElectrodes)/2.);
    }
    else
    {
      iClusterCentralStrip = -static_cast<Int_t>(dClusterPositionAcrossStrips/fdCellAcrossWidth - static_cast<Double_t>(fiNumberReadoutElectrodes)/2.);
    }

    Double_t dLocalHitCoordinates[3] = {0.,0.,0.};
    Double_t dGlobalHitCoordinates[3] = {0.,0.,0.};

    if(fbCellsAlongX)
    {
      dLocalHitCoordinates[0] = dClusterPositionAlongStrips;
      dLocalHitCoordinates[1] = dClusterPositionAcrossStrips;
    }
    else
    {
      dLocalHitCoordinates[0] = dClusterPositionAcrossStrips;
      dLocalHitCoordinates[1] = dClusterPositionAlongStrips;
    }

    fCentralPN->GetMatrix()->LocalToMaster(dLocalHitCoordinates, dGlobalHitCoordinates);


    tHit->SetX(dGlobalHitCoordinates[0]);
    tHit->SetY(dGlobalHitCoordinates[1]);
    tHit->SetZ(dGlobalHitCoordinates[2]);
    tHit->SetTime(dClusterTime);

    if(bCompatibilityMode)
    {
      tHit->SetAddress(CbmTofAddress::GetUniqueAddress(fiModuleIndex, fiCounterIndex, iClusterCentralStrip, 0, fiModuleType));
      tHit->SetFlag(2*iClusterSize);
      tHit->SetChannel(10.*dClusterSummedToT);
      tHit->SetDx(0.5);
      tHit->SetDy(0.5);
      tHit->SetDz(0.5);
      tHit->SetTimeError(dClusterTimeError);
    }
    else
    {
      tHit->SetAddress(fiCounterAddress);
      tHit->SetChannel(iClusterCentralStrip);
    }
  }

}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
Int_t CbmTofCounter::GenerateDarkRate(Double_t dStopTime)
{
  LOG(INFO)<<Form("Dark rate stop time: %.3f", dStopTime)<<FairLogger::endl;

  Int_t iNDarkRatePoints(0);
  Double_t dLocalPointCoordinates[3] = {0., 0., 0.};
  Double_t dGlobalPointCoordinates[3] = {0., 0., 0.};
  Double_t dPointPositionAcrossStrips(0.);
  Double_t dPointPositionAlongStrips(0.);
  CbmTofPoint* tPoint(NULL);

  while(dStopTime > fdLastDarkPointTime)
  {
    dPointPositionAcrossStrips = gRandom->Uniform(-0.5*fdGapAcrossWidth, 0.5*fdGapAcrossWidth);
    dPointPositionAlongStrips  = gRandom->Uniform(-0.5*fdGapAlongWidth, 0.5*fdGapAlongWidth);

    if(fbCellsAlongX)
    {
      dLocalPointCoordinates[0] = dPointPositionAlongStrips;
      dLocalPointCoordinates[1] = dPointPositionAcrossStrips;
    }
    else
    {
      dLocalPointCoordinates[0] = dPointPositionAcrossStrips;
      dLocalPointCoordinates[1] = dPointPositionAlongStrips;
    }

    LocalToMaster(dLocalPointCoordinates, dGlobalPointCoordinates);

    tPoint = new CbmTofPoint();
    tPoint->SetTrackID(-1);
    // Although 'FairMCPoint::fEventId' is an unsigned integer, it can be
    // safely assumed that (UInt_t)(-1) = 4,294,967,295 is a unique index
    // in any reasonable simulation scenario to unambiguously identify
    // dark points in 'CbmTofDigitize'.
    tPoint->SetEventID(-1);
    tPoint->SetDetectorID(fiCounterAddress);
    tPoint->SetTime(fdLastDarkPointTime);
    tPoint->SetX(dGlobalPointCoordinates[0]);
    tPoint->SetY(dGlobalPointCoordinates[1]);
    tPoint->SetZ(dGlobalPointCoordinates[2]);

    CbmMCBuffer::Instance()->Fill(tPoint, kTof);

    // By the previous instruction, the 'CbmTofPoint' object was copied into
    // the 'std::multiset' member of 'CbmMCPointBuffer'. For this reason, it
    // can(/needs to) be deleted here.
    delete tPoint;

    fdLastDarkPointTime += gRandom->Exp(fdMeanDarkPointInterval);

    iNDarkRatePoints++;
  }

  return iNDarkRatePoints;
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
Double_t CbmTofCounter::DistributeCharge(const Double_t& dTotalCharge, const Double_t& dChargeDistance, const Double_t& dX1, const Double_t& dX2, const Double_t& dY1, const Double_t& dY2)
{
  Double_t result = dTotalCharge/TMath::TwoPi()*( TMath::ATan(dX2*dY2/dChargeDistance/TMath::Sqrt(TMath::Power(dChargeDistance,2.)+TMath::Power(dX2,2.)+TMath::Power(dY2,2.)))
                                                 -TMath::ATan(dX1*dY2/dChargeDistance/TMath::Sqrt(TMath::Power(dChargeDistance,2.)+TMath::Power(dX1,2.)+TMath::Power(dY2,2.)))
                                                 -TMath::ATan(dX2*dY1/dChargeDistance/TMath::Sqrt(TMath::Power(dChargeDistance,2.)+TMath::Power(dX2,2.)+TMath::Power(dY1,2.)))
                                                 +TMath::ATan(dX1*dY1/dChargeDistance/TMath::Sqrt(TMath::Power(dChargeDistance,2.)+TMath::Power(dX1,2.)+TMath::Power(dY1,2.)))
                                                );
  return result;
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
Double_t CbmTofCounter::ChargeLimit(const Double_t& dChargeDistance, const Double_t& dX1, const Double_t& dX2, const Double_t& dY1, const Double_t& dY2)
{
  Double_t result = TMath::TwoPi()/
                                  ( TMath::ATan(dX2*dY2/dChargeDistance/TMath::Sqrt(TMath::Power(dChargeDistance,2.)+TMath::Power(dX2,2.)+TMath::Power(dY2,2.)))
                                   -TMath::ATan(dX1*dY2/dChargeDistance/TMath::Sqrt(TMath::Power(dChargeDistance,2.)+TMath::Power(dX1,2.)+TMath::Power(dY2,2.)))
                                   -TMath::ATan(dX2*dY1/dChargeDistance/TMath::Sqrt(TMath::Power(dChargeDistance,2.)+TMath::Power(dX2,2.)+TMath::Power(dY1,2.)))
                                   +TMath::ATan(dX1*dY1/dChargeDistance/TMath::Sqrt(TMath::Power(dChargeDistance,2.)+TMath::Power(dX1,2.)+TMath::Power(dY1,2.)))
                                  );
  return result;
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
Double_t CbmTofCounter::NegativeLandau(Double_t x, void* params)
{
  Double_t* p = (Double_t *)params;

  return -TMath::Landau(x,p[0],p[1],kTRUE);
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
Double_t CbmTofCounter::LandauThreshold(Double_t x, void* params)
{
  Double_t* p = (Double_t *)params;

  return p[3]-p[2]*TMath::Landau(x,p[0],p[1],kTRUE);
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
Double_t CbmTofCounter::LandauSignal(Double_t x, void* params)
{
  Double_t* p = (Double_t *)params;

  return p[2]*TMath::Landau(x,p[0],p[1],kTRUE);
}
// ---------------------------------------------------------------------------


ClassImp(CbmTofCounter)