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

#include "CbmTofRecoMCQa.h"

#include "CbmMatch.h"
#include "CbmTofPoint.h"
#include "CbmTofDigiExp.h"
#include "CbmTofHit.h"
#include "CbmTofWall.h"
#include "CbmTofCounter.h"
#include "CbmTofDigiTbPar.h"
#include "CbmTofGeoHandler.h"
#include "CbmTofDef.h"

#include "FairLogger.h"
#include "FairRootManager.h"
#include "FairRun.h"
#include "FairRuntimeDb.h"
#include "FairFileSource.h"

#include "TClonesArray.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TEfficiency.h"
#include "TFile.h"
#include "TMath.h"

#include <cmath>

// ---------------------------------------------------------------------------
CbmTofRecoMCQa::CbmTofRecoMCQa() :
	FairTask("TofRecoMCQa"),
  fFileSource(NULL),
  fTofPointsTB(NULL),
  fTofDigis(NULL),
  fTofHits(NULL),
  fTofHitDigiMatches(NULL),
  fTofHitPointMatches(NULL),
  fTofPointMatches(NULL),
  fDigiTbParSet(NULL),
  fiNTimeAxisBins(10),
  fdTimeAxisBinWidth(0.1),
  fdTimeAxisRange(1.),
  fClusterSizes(),
  fStripToTSpectra(),
  fStripPositionSpectra(),
  fPointPositionSpectra(),
  fHitPositionSpectra(),
  fStripTimeResidualWalk(),
  fHitAcrossPositionResidual(),
  fHitAlongPositionResidual(),
  fHitTimeResidual(),
  fCentralClusterSizeEvolution(),
  fCentralTimeResidualEvolution(),
  fHitChiSquare(),
  fCorrectedHitChiSquare(),
  fClusterSizeEvolution(),
  fAcrossPositionResidualEvolution(),
  fAlongPositionResidualEvolution(),
  fTimeResidualEvolution(),
  fCentralStripToTEvolution(),
  fIntegratedEfficiency(),
  fChiSquareEfficiency(),
  fCorrectedChiSquareEfficiency(),
  fClusterSizesSPH(),
  fStripToTSpectraSPH(),
  fStripPositionSpectraSPH(),
  fHitPositionSpectraSPH(),
  fStripTimeResidualWalkSPH(),
  fHitAcrossPositionResidualSPH(),
  fHitAlongPositionResidualSPH(),
  fHitTimeResidualSPH(),
  fHitChiSquareSPH(),
  fCorrectedHitChiSquareSPH(),
  fClusterSizeEvolutionSPH(),
  fAcrossPositionResidualEvolutionSPH(),
  fAlongPositionResidualEvolutionSPH(),
  fTimeResidualEvolutionSPH(),
  fCentralStripToTEvolutionSPH(),
  fDistanceToLastPoint(),
  fDistanceToLocalLastPoint(),
  fDistanceToLastHitUnbiased(),
  fDistanceToLastHitBiased(),
  fEfficiencyToLastHitUnbiased(),
  fEfficiencyToLastHitBiased(),
  fClusterSizeToLastHitUnbiased(),
  fClusterSizeToLastHitBiased(),
  fTimeResidualToLastHitUnbiased(),
  fTimeResidualToLastHitBiased(),
  fLocalSurfaceFlux(),
  fChannelDigiRates(),
  fNHitsPerPoint(),
  fNPointsPerHit(),
  fAcrossPositionResidualEvolutionSigma(),
  fAlongPositionResidualEvolutionSigma(),
  fTimeResidualEvolutionSigma(),
  fAcrossPositionResidualEvolutionSigmaSPH(),
  fAlongPositionResidualEvolutionSigmaSPH(),
  fTimeResidualEvolutionSigmaSPH(),
  fHistogramFile(NULL),
  fHistogramFileName(),
  fiTotalNPoints(0),
  fiTotalNMergedPoints(0),
  fbAverageSpills(kTRUE),
  fdSpillLength(1.),
  fdSpillBreakLength(0.),
  fdDataTimeOffset(0.),
  fdLocalRadius(0.5),
  fbFillLocalHistos(kFALSE),
  fiNSerialSpills(1),
  fdPreviousPointTime(),
  fPointLinkNbSet(),
  fbHitCompatibilityMode(kFALSE),
  fbAlternativeBranchNames(kFALSE),
  fbConsiderPointClusters(kFALSE),
  fbRealClustering(kFALSE),
  fbConsiderRefTracks(kFALSE),
  fProcessedHits(),
  fProcessedDigis(),
  fdChiSquareLimit(TMath::Limits<Double_t>::Max())
{
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
CbmTofRecoMCQa::~CbmTofRecoMCQa()
{
  fPointLinkNbSet.clear();
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
void CbmTofRecoMCQa::Exec(Option_t*)
{
  // HINT: The case '1 < tPointMatch->GetNofLinks()' which covers efficient
  //       MC points that contribute to multiple ToF hits by digi merging
  //       might not be ignorable for efficiency histograms etc. at higher
  //       particle fluxes.

  Int_t iInputTSIndex = FairRootManager::Instance()->GetEntryNr();
  fProcessedHits.clear();
  fProcessedDigis.clear();

  for (Int_t iPoint = 0; iPoint < fTofPointsTB->GetEntriesFast(); iPoint++)
  {
    CbmTofPoint* tPoint = dynamic_cast<CbmTofPoint*>(fTofPointsTB->At(iPoint));

    if(fbConsiderPointClusters && tof::point_Internal == tPoint->GetUniqueID())
    {
      continue;
    }
    else if(fbConsiderRefTracks && tof::point_RefTrack != tPoint->GetUniqueID())
    {
      continue;
    }

    Int_t iCounterAddress = tPoint->GetDetectorID();
    Double_t dPointTime = tPoint->GetTime() - fdDataTimeOffset;

    fDistanceToLastPoint.at(iCounterAddress)->Fill(dPointTime - fdPreviousPointTime.at(iCounterAddress));
    fdPreviousPointTime.at(iCounterAddress) = dPointTime;

    CbmTofCounter* tCounter = CbmTofWall::Instance()->GetCounter(iCounterAddress);

    Double_t dGlobalPointCoordinates[3] = {tPoint->GetX(), tPoint->GetY(), tPoint->GetZ()};
    Double_t dLocalPointCoordinates[3] = {0., 0., 0.};
    tCounter->MasterToLocal(dGlobalPointCoordinates, dLocalPointCoordinates);

    if(tCounter->GetCellsAlongX())
    {
      fPointPositionSpectra.at(iCounterAddress)->Fill(dLocalPointCoordinates[1], dLocalPointCoordinates[0]);
    }
    else
    {
      fPointPositionSpectra.at(iCounterAddress)->Fill(dLocalPointCoordinates[0], dLocalPointCoordinates[1]);
    }

    CbmMatch* tPointMatch = dynamic_cast<CbmMatch*>(fTofPointMatches->At(iPoint));

    Int_t iNHitsPerPoint = tPointMatch->GetNofLinks();

    fNHitsPerPoint.at(iCounterAddress)->Fill(iNHitsPerPoint);

    if(0 == iNHitsPerPoint)
    {
      fIntegratedEfficiency.at(iCounterAddress)->Fill(kFALSE, fmod(dPointTime/1e9, fdTimeAxisRange));
      fChiSquareEfficiency.at(iCounterAddress)->Fill(kFALSE, fmod(dPointTime/1e9, fdTimeAxisRange));
      fCorrectedChiSquareEfficiency.at(iCounterAddress)->Fill(kFALSE, fmod(dPointTime/1e9, fdTimeAxisRange));

      if(tCounter->GetCellsAlongX())
      {
        fLocalSurfaceFlux.at(iCounterAddress)->Fill(kFALSE, dLocalPointCoordinates[1], dLocalPointCoordinates[0], fmod(dPointTime/1e9, fdTimeAxisRange));
      }
      else
      {
        fLocalSurfaceFlux.at(iCounterAddress)->Fill(kFALSE, dLocalPointCoordinates[0], dLocalPointCoordinates[1], fmod(dPointTime/1e9, fdTimeAxisRange));
      }

      if(fbFillLocalHistos)
      {
        // Scan - if necessary - the entire "TofPointTB" branch in reverse order to find
        // the closest efficient point backwards in time.
        Int_t iIndexPrevPoint(-1);
        Int_t iTimeSlicePrevPoint(-1);
        Int_t iIndexPrevEffPoint(-1);
        Int_t iTimeSlicePrevEffPoint(-1);

        Double_t dGlobalPrevPointCoordinates[3] = {0., 0., 0.};
        Double_t dLocalPrevPointCoordinates[3] = {0., 0., 0.};

        // Start form the current time slice.
        for(Int_t iPrevPoint = iPoint - 1; iPrevPoint >= 0; iPrevPoint--)
        {
          CbmTofPoint* tPrevPoint = dynamic_cast<CbmTofPoint*>(fTofPointsTB->At(iPrevPoint));
          Int_t iPrevCounterAddress = tPrevPoint->GetDetectorID();

          if(iCounterAddress == iPrevCounterAddress)
          {
            dGlobalPrevPointCoordinates[0] = tPrevPoint->GetX();
            dGlobalPrevPointCoordinates[1] = tPrevPoint->GetY();
            dGlobalPrevPointCoordinates[2] = tPrevPoint->GetZ();
            tCounter->MasterToLocal(dGlobalPrevPointCoordinates, dLocalPrevPointCoordinates);

            if(TMath::Power(fdLocalRadius, 2.) > TMath::Power(dLocalPointCoordinates[0] - dLocalPrevPointCoordinates[0], 2.)
                                                +TMath::Power(dLocalPointCoordinates[1] - dLocalPrevPointCoordinates[1], 2.))
            {
              if(-1 == iIndexPrevPoint)
              {
                iIndexPrevPoint = iPrevPoint;
                iTimeSlicePrevPoint = iInputTSIndex;

                fDistanceToLocalLastPoint.at(iCounterAddress)->Fill(dPointTime - (tPrevPoint->GetTime() - fdDataTimeOffset));
              }

              CbmMatch* tPrevPointMatch = dynamic_cast<CbmMatch*>(fTofPointMatches->At(iPrevPoint));
              if(0 < tPrevPointMatch->GetNofLinks())
              {
                iIndexPrevEffPoint = iPrevPoint;
                break;
              }
            }
          }
        }

        // Found previous efficient point in the current time slice.
        if(-1 != iIndexPrevEffPoint)
        {
          CbmTofPoint* tPrevEffPoint = dynamic_cast<CbmTofPoint*>(fTofPointsTB->At(iIndexPrevEffPoint));
          Double_t dTimeDifference = dPointTime - (tPrevEffPoint->GetTime() - fdDataTimeOffset);

          // No inefficient point between two efficient ones.
          if(iIndexPrevPoint == iIndexPrevEffPoint)
          {
            fDistanceToLastHitUnbiased.at(iCounterAddress)->Fill(dTimeDifference);
            fEfficiencyToLastHitUnbiased.at(iCounterAddress)->Fill(kFALSE, dTimeDifference);
          }

          // Possibly some inefficient points between two efficient ones.
          fDistanceToLastHitBiased.at(iCounterAddress)->Fill(dTimeDifference);
          fEfficiencyToLastHitBiased.at(iCounterAddress)->Fill(kFALSE, dTimeDifference);

        }
        // Need to scan earlier time slices for the previous efficient point.
        else
        {
          TChain* tChain = fFileSource->GetInChain();

          for(Int_t iTimeSlice = iInputTSIndex - 1; iTimeSlice >= 0; iTimeSlice--)
          {
            // Load only the required "TofPointTB.*" branch entries into memory!
            tChain->SetBranchStatus("*", kFALSE);
            tChain->SetBranchStatus("TofPointTB.fTime", kTRUE);
            tChain->SetBranchStatus("TofPointTB.fDetectorID", kTRUE);
            tChain->SetBranchStatus("TofPointTB.fX", kTRUE);
            tChain->SetBranchStatus("TofPointTB.fY", kTRUE);
            tChain->GetEntry(iTimeSlice);

            for(Int_t iPrevPoint = fTofPointsTB->GetEntriesFast() - 1; iPrevPoint >= 0; iPrevPoint--)
            {
              CbmTofPoint* tPrevPoint = dynamic_cast<CbmTofPoint*>(fTofPointsTB->At(iPrevPoint));
              Int_t iPrevCounterAddress = tPrevPoint->GetDetectorID();

              if(iCounterAddress == iPrevCounterAddress)
              {
                dGlobalPrevPointCoordinates[0] = tPrevPoint->GetX();
                dGlobalPrevPointCoordinates[1] = tPrevPoint->GetY();
                dGlobalPrevPointCoordinates[2] = tPrevPoint->GetZ();
                tCounter->MasterToLocal(dGlobalPrevPointCoordinates, dLocalPrevPointCoordinates);

                if(TMath::Power(fdLocalRadius, 2.) > TMath::Power(dLocalPointCoordinates[0] - dLocalPrevPointCoordinates[0], 2.)
                                                    +TMath::Power(dLocalPointCoordinates[1] - dLocalPrevPointCoordinates[1], 2.))
                {
                  if(-1 == iIndexPrevPoint)
                  {
                    iIndexPrevPoint = iPrevPoint;
                    iTimeSlicePrevPoint = iTimeSlice;

                    fDistanceToLocalLastPoint.at(iCounterAddress)->Fill(dPointTime - (tPrevPoint->GetTime() - fdDataTimeOffset));
                  }

                  if(fPointLinkNbSet.count(std::make_pair(iTimeSlice, iPrevPoint)))
                  {
                    iIndexPrevEffPoint = iPrevPoint;
                    iTimeSlicePrevEffPoint = iTimeSlice;
                    break;
                  }
                }
              }
            }

            if(-1 != iIndexPrevEffPoint)
            {
              break;
            }
          }

          // Found previous efficient point.
          if(-1 != iIndexPrevEffPoint)
          {
            CbmTofPoint* tPrevEffPoint = dynamic_cast<CbmTofPoint*>(fTofPointsTB->At(iIndexPrevEffPoint));
            Double_t dTimeDifference = dPointTime - (tPrevEffPoint->GetTime() - fdDataTimeOffset);

            if(iIndexPrevPoint == iIndexPrevEffPoint)
            {
              // No inefficient point between two efficient ones.
              if(iTimeSlicePrevPoint == iTimeSlicePrevEffPoint)
              {
                fDistanceToLastHitUnbiased.at(iCounterAddress)->Fill(dTimeDifference);
                fEfficiencyToLastHitUnbiased.at(iCounterAddress)->Fill(kFALSE, dTimeDifference);
              }
            }

            // Possibly some inefficient points between two efficient ones.
            fDistanceToLastHitBiased.at(iCounterAddress)->Fill(dTimeDifference);
            fEfficiencyToLastHitBiased.at(iCounterAddress)->Fill(kFALSE, dTimeDifference);

          }

          // Reload the current "TofPointTB" branch entry into memory.
          if(0 < iInputTSIndex)
          {
            tChain->GetEntry(iInputTSIndex);
            tChain->SetBranchStatus("*", kTRUE);
          }

        }

      }
    }
    else if(0 < iNHitsPerPoint)
    {
      fIntegratedEfficiency.at(iCounterAddress)->Fill(kTRUE, fmod(dPointTime/1e9, fdTimeAxisRange));

      if(tCounter->GetCellsAlongX())
      {
        fLocalSurfaceFlux.at(iCounterAddress)->Fill(kTRUE, dLocalPointCoordinates[1], dLocalPointCoordinates[0], fmod(dPointTime/1e9, fdTimeAxisRange));
      }
      else
      {
        fLocalSurfaceFlux.at(iCounterAddress)->Fill(kTRUE, dLocalPointCoordinates[0], dLocalPointCoordinates[1], fmod(dPointTime/1e9, fdTimeAxisRange));
      }

      // IMPORTANT NOTE: Information from a single MC point can contribute to multiple
      // MC clusters because of signal interference during digitization. Here, we 
      // assume that the hit which contains most information from a point,
      // i.e. which has the largest 'CbmLink::fWeight' among all point to hit forward
      // references, is the corresponding one that 'CbmTofClusterMC' uniquely built from
      // all digis containing primary MC backlinks to the point.
      const CbmLink& tHitLink = tPointMatch->GetMatchedLink(); // LEGACY: tPointMatch->GetLink(0);
      CbmTofHit* tHit = dynamic_cast<CbmTofHit*>(fTofHits->At(tHitLink.GetIndex()));

      // Make sure to not fill information from a single hit into QA histograms multiple times
      // if this hit is referenced by several MC points as best match.
      Bool_t bHitNotYetProcessed = (fProcessedHits.emplace(tHitLink.GetFile(), tHitLink.GetEntry(), tHitLink.GetIndex())).second;

      CbmMatch* tHitPointMatch = dynamic_cast<CbmMatch*>(fTofHitPointMatches->At(tHitLink.GetIndex()));
      Bool_t bSinglePointHit = (1 == tHitPointMatch->GetNofLinks());

      if(bHitNotYetProcessed)
      {
        fNPointsPerHit.at(iCounterAddress)->Fill(tHitPointMatch->GetNofLinks());
      }

      CbmMatch* tHitDigiMatch = dynamic_cast<CbmMatch*>(fTofHitDigiMatches->At(tHitLink.GetIndex()));

      Int_t iNDigiLinks = tHitDigiMatch->GetNofLinks();
      Int_t iClusterSize = iNDigiLinks/2;
      Double_t dHitTime = tHit->GetTime() - fdDataTimeOffset;
      if(fbRealClustering)
      {
        dHitTime -= 0.5*tCounter->GetCellAlongWidth()/tCounter->GetPropagationVelocity();
      }

      Int_t iHitMainChannel(0);
      if(fbHitCompatibilityMode)
      {
        iHitMainChannel = CbmTofAddress::GetChannelId(tHit->GetAddress());
      }
      else
      {
        iHitMainChannel = tHit->GetCh();
      }

      for(Int_t iDigiLink = 0; iDigiLink < iNDigiLinks; iDigiLink += 2)
      {
        // 'CbmTofDigi' objects referred to by the link objects should be sorted with respect to
        // increasing channel side indices by 'CbmTofCounter::BuildClusterMC'.
        // In 'CbmTofTestBeamClusterizer', however, digi link pairs are sorted with respect to
        // time.
        const CbmLink& tEvenDigiLink = tHitDigiMatch->GetLink(iDigiLink);
        const CbmLink& tOddDigiLink  = tHitDigiMatch->GetLink(iDigiLink + 1);

        Bool_t bEvenDigiNotYetProcessed = (fProcessedDigis.emplace(tEvenDigiLink.GetFile(), tEvenDigiLink.GetEntry(), tEvenDigiLink.GetIndex())).second;
        Bool_t bOddDigiNotYetProcessed = (fProcessedDigis.emplace(tOddDigiLink.GetFile(), tOddDigiLink.GetEntry(), tOddDigiLink.GetIndex())).second;

        CbmTofDigiExp* tEvenDigi = dynamic_cast<CbmTofDigiExp*>(fTofDigis->At(tEvenDigiLink.GetIndex()));
        CbmTofDigiExp* tOddDigi  = dynamic_cast<CbmTofDigiExp*>(fTofDigis->At(tOddDigiLink.GetIndex()));

        CbmTofDigiExp* tBottomDigi(NULL);
        CbmTofDigiExp* tTopDigi(NULL);
        Bool_t bBottomDigiNotYetProcessed(kTRUE);
        Bool_t bTopDigiNotYetProcessed(kTRUE);

        if(0 == static_cast<Int_t>(tEvenDigi->GetSide()))
        {
          if(fbRealClustering)
          {
            tBottomDigi = tOddDigi;
            tTopDigi    = tEvenDigi;

            bBottomDigiNotYetProcessed = bOddDigiNotYetProcessed;
            bTopDigiNotYetProcessed = bEvenDigiNotYetProcessed;
          }
          else
          {
            tBottomDigi = tEvenDigi;
            tTopDigi    = tOddDigi;

            bBottomDigiNotYetProcessed = bEvenDigiNotYetProcessed;
            bTopDigiNotYetProcessed = bOddDigiNotYetProcessed;
          }
        }
        else
        {
          if(fbRealClustering)
          {
            tBottomDigi = tEvenDigi;
            tTopDigi    = tOddDigi;

            bBottomDigiNotYetProcessed = bEvenDigiNotYetProcessed;
            bTopDigiNotYetProcessed = bOddDigiNotYetProcessed;
          }
          else
          {
            tBottomDigi = tOddDigi;
            tTopDigi    = tEvenDigi;

            bBottomDigiNotYetProcessed = bOddDigiNotYetProcessed;
            bTopDigiNotYetProcessed = bEvenDigiNotYetProcessed;
          }
        }

        // remove the fixed delay of the reconstructed strip time w.r.t. the MC point time
        // due to the finite signal propagation time on the readout cells
        Double_t dStripTime = 0.5*(tBottomDigi->GetTime() - fdDataTimeOffset + tTopDigi->GetTime() - fdDataTimeOffset)
                             -0.5*tCounter->GetCellAlongWidth()/tCounter->GetPropagationVelocity();
        Double_t dStripAlongPosition = 0.5*tCounter->GetPropagationVelocity()*(tBottomDigi->GetTime() - fdDataTimeOffset - (tTopDigi->GetTime() - fdDataTimeOffset));

        if(bHitNotYetProcessed)
        {
          // In case of multiple MC points close in time and space creating multiple non-identical ideal hits
          // on the counter which overlap by at least one digi pair, respectively, the digi information that
          // contributes to more than one hit would enter the ToT histograms multiple times.
          if(bBottomDigiNotYetProcessed)
          {
            fStripToTSpectra.at(iCounterAddress)->Fill(2*tBottomDigi->GetChannel() + tBottomDigi->GetSide(), tBottomDigi->GetTot());

            if(tBottomDigi->GetChannel() == tCounter->GetNumberReadoutElectrodes()/2 - 1)
            {
              fCentralStripToTEvolution.at(iCounterAddress)->Fill(fmod(dPointTime/1e9, fdTimeAxisRange), tBottomDigi->GetTot());
            }
          }

          if(bTopDigiNotYetProcessed)
          {
            fStripToTSpectra.at(iCounterAddress)->Fill(2*tTopDigi->GetChannel() + tTopDigi->GetSide(), tTopDigi->GetTot());
          }

          // TODO: Here, one would need to check for the exact same digi pair (bottom and top) not to
          //       have entered the position histogram previously due to digi pair sharing among neighboring hits.
          fStripPositionSpectra.at(iCounterAddress)->Fill(tBottomDigi->GetChannel(), dStripAlongPosition);

          if(bSinglePointHit)
          {
            fStripToTSpectraSPH.at(iCounterAddress)->Fill(2*tBottomDigi->GetChannel() + tBottomDigi->GetSide(), tBottomDigi->GetTot());
            fStripToTSpectraSPH.at(iCounterAddress)->Fill(2*tTopDigi->GetChannel() + tTopDigi->GetSide(), tTopDigi->GetTot());

            if(tBottomDigi->GetChannel() == tCounter->GetNumberReadoutElectrodes()/2 - 1)
            {
              fCentralStripToTEvolutionSPH.at(iCounterAddress)->Fill(fmod(dPointTime/1e9, fdTimeAxisRange), tBottomDigi->GetTot());
            }

            fStripPositionSpectraSPH.at(iCounterAddress)->Fill(tBottomDigi->GetChannel(), dStripAlongPosition);
          }
        }

        fStripTimeResidualWalk.at(iCounterAddress)->Fill(2*tBottomDigi->GetChannel() + tBottomDigi->GetSide(), tBottomDigi->GetTot(), dStripTime - dPointTime);
        fStripTimeResidualWalk.at(iCounterAddress)->Fill(2*tTopDigi->GetChannel() + tTopDigi->GetSide(), tTopDigi->GetTot(), dStripTime - dPointTime);

        if(bSinglePointHit)
        {
          fStripTimeResidualWalkSPH.at(iCounterAddress)->Fill(2*tBottomDigi->GetChannel() + tBottomDigi->GetSide(), tBottomDigi->GetTot(), dStripTime - dPointTime);
          fStripTimeResidualWalkSPH.at(iCounterAddress)->Fill(2*tTopDigi->GetChannel() + tTopDigi->GetSide(), tTopDigi->GetTot(), dStripTime - dPointTime);
        }
      }

      if(bHitNotYetProcessed)
      {
        fClusterSizeEvolution.at(iCounterAddress)->Fill(fmod(dPointTime/1e9, fdTimeAxisRange), iClusterSize);

        if(bSinglePointHit)
        {
          fClusterSizeEvolutionSPH.at(iCounterAddress)->Fill(fmod(dPointTime/1e9, fdTimeAxisRange), iClusterSize);
        }
      }
      fTimeResidualEvolution.at(iCounterAddress)->Fill(fmod(dPointTime/1e9, fdTimeAxisRange), dHitTime - dPointTime);

      if(bSinglePointHit)
      {
        fTimeResidualEvolutionSPH.at(iCounterAddress)->Fill(fmod(dPointTime/1e9, fdTimeAxisRange), dHitTime - dPointTime);
      }

      if(iHitMainChannel == (tCounter->GetNumberReadoutElectrodes()/2 - 1))
      {
        if(bHitNotYetProcessed)
        {
          fCentralClusterSizeEvolution.at(iCounterAddress)->Fill(fmod(dPointTime/1e9, fdTimeAxisRange), iClusterSize);
        }
        fCentralTimeResidualEvolution.at(iCounterAddress)->Fill(fmod(dPointTime/1e9, fdTimeAxisRange), dHitTime - dPointTime);
      }

      if(bHitNotYetProcessed)
      {
        fClusterSizes.at(iCounterAddress)->Fill(iHitMainChannel, iClusterSize);

        if(bSinglePointHit)
        {
          fClusterSizesSPH.at(iCounterAddress)->Fill(iHitMainChannel, iClusterSize);
        }
      }


      Double_t dGlobalHitCoordinates[3] = {tHit->GetX(), tHit->GetY(), tHit->GetZ()};
      Double_t dLocalHitCoordinates[3] = {0., 0., 0.};
      tCounter->MasterToLocal(dGlobalHitCoordinates, dLocalHitCoordinates);

      Double_t dHitChiSquare(TMath::Limits<Double_t>::Max());
      Double_t dCorrectedHitChiSquare(TMath::Limits<Double_t>::Max());
      Double_t dHitChiSquareSPH(TMath::Limits<Double_t>::Max());
      Double_t dCorrectedHitChiSquareSPH(TMath::Limits<Double_t>::Max());

      if(tCounter->GetCellsAlongX())
      {
        if(bHitNotYetProcessed)
        {
          fHitPositionSpectra.at(iCounterAddress)->Fill(dLocalHitCoordinates[1], dLocalHitCoordinates[0]);

          if(bSinglePointHit)
          {
            fHitPositionSpectraSPH.at(iCounterAddress)->Fill(dLocalHitCoordinates[1], dLocalHitCoordinates[0]);
          }
        }
        fHitAcrossPositionResidual.at(iCounterAddress)->Fill(dLocalHitCoordinates[1] - dLocalPointCoordinates[1]);
        fHitAlongPositionResidual.at(iCounterAddress)->Fill(dLocalHitCoordinates[0] - dLocalPointCoordinates[0]);
        fAcrossPositionResidualEvolution.at(iCounterAddress)->Fill(fmod(dPointTime/1e9, fdTimeAxisRange), dLocalHitCoordinates[1] - dLocalPointCoordinates[1]);
        fAlongPositionResidualEvolution.at(iCounterAddress)->Fill(fmod(dPointTime/1e9, fdTimeAxisRange), dLocalHitCoordinates[0] - dLocalPointCoordinates[0]);

        // TODO: In case of multiple points contributing to a hit the residual sigmas derived during parametrization against a MC reference
        //       do not accurately describe the (larger) sigmas that result from interfering point information in hits. To compensate for this
        //       effect one could use the residual sigmas obtained from the first TIS bin for the entire run which comes at the cost of
        //       relatively low statistics for the residual fits, though.
        dHitChiSquare = TMath::Power((dLocalHitCoordinates[1] - dLocalPointCoordinates[1])/tCounter->GetAcrossPositionResidualSigma(), 2.)
                       +TMath::Power((dLocalHitCoordinates[0] - dLocalPointCoordinates[0])/tCounter->GetAlongPositionResidualSigma(), 2.)
                       +TMath::Power((dHitTime - dPointTime)/tCounter->GetTimeResidualSigma(), 2.);

        dCorrectedHitChiSquare = TMath::Power((dLocalHitCoordinates[1] - dLocalPointCoordinates[1])/GetAcrossPositionResidualSigma(iCounterAddress, fmod(dPointTime/1e9, fdTimeAxisRange)), 2.)
                                +TMath::Power((dLocalHitCoordinates[0] - dLocalPointCoordinates[0])/GetAlongPositionResidualSigma(iCounterAddress, fmod(dPointTime/1e9, fdTimeAxisRange)), 2.)
                                +TMath::Power((dHitTime - dPointTime)/GetTimeResidualSigma(iCounterAddress, fmod(dPointTime/1e9, fdTimeAxisRange)), 2.);

        dHitChiSquareSPH = TMath::Power((dLocalHitCoordinates[1] - dLocalPointCoordinates[1])/tCounter->GetAcrossPositionResidualSigma(), 2.)
                          +TMath::Power((dLocalHitCoordinates[0] - dLocalPointCoordinates[0])/tCounter->GetAlongPositionResidualSigma(), 2.)
                          +TMath::Power((dHitTime - dPointTime)/tCounter->GetTimeResidualSigma(), 2.);

        dCorrectedHitChiSquareSPH = TMath::Power((dLocalHitCoordinates[1] - dLocalPointCoordinates[1])/GetAcrossPositionResidualSigmaSPH(iCounterAddress, fmod(dPointTime/1e9, fdTimeAxisRange)), 2.)
                                   +TMath::Power((dLocalHitCoordinates[0] - dLocalPointCoordinates[0])/GetAlongPositionResidualSigmaSPH(iCounterAddress, fmod(dPointTime/1e9, fdTimeAxisRange)), 2.)
                                   +TMath::Power((dHitTime - dPointTime)/GetTimeResidualSigmaSPH(iCounterAddress, fmod(dPointTime/1e9, fdTimeAxisRange)), 2.);

        if(bSinglePointHit)
        {
          fHitAcrossPositionResidualSPH.at(iCounterAddress)->Fill(dLocalHitCoordinates[1] - dLocalPointCoordinates[1]);
          fHitAlongPositionResidualSPH.at(iCounterAddress)->Fill(dLocalHitCoordinates[0] - dLocalPointCoordinates[0]);
          fAcrossPositionResidualEvolutionSPH.at(iCounterAddress)->Fill(fmod(dPointTime/1e9, fdTimeAxisRange), dLocalHitCoordinates[1] - dLocalPointCoordinates[1]);
          fAlongPositionResidualEvolutionSPH.at(iCounterAddress)->Fill(fmod(dPointTime/1e9, fdTimeAxisRange), dLocalHitCoordinates[0] - dLocalPointCoordinates[0]);
        }
      }
      else
      {
        if(bHitNotYetProcessed)
        {
          fHitPositionSpectra.at(iCounterAddress)->Fill(dLocalHitCoordinates[0], dLocalHitCoordinates[1]);

          if(bSinglePointHit)
          {
            fHitPositionSpectraSPH.at(iCounterAddress)->Fill(dLocalHitCoordinates[0], dLocalHitCoordinates[1]);
          }
        }
        fHitAcrossPositionResidual.at(iCounterAddress)->Fill(dLocalHitCoordinates[0] - dLocalPointCoordinates[0]);
        fHitAlongPositionResidual.at(iCounterAddress)->Fill(dLocalHitCoordinates[1] - dLocalPointCoordinates[1]);
        fAcrossPositionResidualEvolution.at(iCounterAddress)->Fill(fmod(dPointTime/1e9, fdTimeAxisRange), dLocalHitCoordinates[0] - dLocalPointCoordinates[0]);
        fAlongPositionResidualEvolution.at(iCounterAddress)->Fill(fmod(dPointTime/1e9, fdTimeAxisRange), dLocalHitCoordinates[1] - dLocalPointCoordinates[1]);

        // TODO: In case of multiple points contributing to a hit the residual sigmas derived during parametrization against a MC reference
        //       do not accurately describe the (larger) sigmas that result from interfering point information in hits. To compensate for this
        //       effect one could use the residual sigmas obtained from the first TIS bin for the entire run which comes at the cost of
        //       relatively low statistics for the residual fits, though.
        dHitChiSquare = TMath::Power((dLocalHitCoordinates[0] - dLocalPointCoordinates[0])/tCounter->GetAcrossPositionResidualSigma(), 2.)
                       +TMath::Power((dLocalHitCoordinates[1] - dLocalPointCoordinates[1])/tCounter->GetAlongPositionResidualSigma(), 2.)
                       +TMath::Power((dHitTime - dPointTime)/tCounter->GetTimeResidualSigma(), 2.);

        dCorrectedHitChiSquare = TMath::Power((dLocalHitCoordinates[0] - dLocalPointCoordinates[0])/GetAcrossPositionResidualSigma(iCounterAddress, fmod(dPointTime/1e9, fdTimeAxisRange)), 2.)
                                +TMath::Power((dLocalHitCoordinates[1] - dLocalPointCoordinates[1])/GetAlongPositionResidualSigma(iCounterAddress, fmod(dPointTime/1e9, fdTimeAxisRange)), 2.)
                                +TMath::Power((dHitTime - dPointTime)/GetTimeResidualSigma(iCounterAddress, fmod(dPointTime/1e9, fdTimeAxisRange)), 2.);

        dHitChiSquareSPH = TMath::Power((dLocalHitCoordinates[0] - dLocalPointCoordinates[0])/tCounter->GetAcrossPositionResidualSigma(), 2.)
                          +TMath::Power((dLocalHitCoordinates[1] - dLocalPointCoordinates[1])/tCounter->GetAlongPositionResidualSigma(), 2.)
                          +TMath::Power((dHitTime - dPointTime)/tCounter->GetTimeResidualSigma(), 2.);

        dCorrectedHitChiSquareSPH = TMath::Power((dLocalHitCoordinates[0] - dLocalPointCoordinates[0])/GetAcrossPositionResidualSigmaSPH(iCounterAddress, fmod(dPointTime/1e9, fdTimeAxisRange)), 2.)
                                   +TMath::Power((dLocalHitCoordinates[1] - dLocalPointCoordinates[1])/GetAlongPositionResidualSigmaSPH(iCounterAddress, fmod(dPointTime/1e9, fdTimeAxisRange)), 2.)
                                   +TMath::Power((dHitTime - dPointTime)/GetTimeResidualSigmaSPH(iCounterAddress, fmod(dPointTime/1e9, fdTimeAxisRange)), 2.);

        if(bSinglePointHit)
        {
          fHitAcrossPositionResidualSPH.at(iCounterAddress)->Fill(dLocalHitCoordinates[0] - dLocalPointCoordinates[0]);
          fHitAlongPositionResidualSPH.at(iCounterAddress)->Fill(dLocalHitCoordinates[1] - dLocalPointCoordinates[1]);
          fAcrossPositionResidualEvolutionSPH.at(iCounterAddress)->Fill(fmod(dPointTime/1e9, fdTimeAxisRange), dLocalHitCoordinates[0] - dLocalPointCoordinates[0]);
          fAlongPositionResidualEvolutionSPH.at(iCounterAddress)->Fill(fmod(dPointTime/1e9, fdTimeAxisRange), dLocalHitCoordinates[1] - dLocalPointCoordinates[1]);
        }
      }

      fHitChiSquare.at(iCounterAddress)->Fill(dHitChiSquare);
      fCorrectedHitChiSquare.at(iCounterAddress)->Fill(dCorrectedHitChiSquare);
      if(bSinglePointHit)
      {
        fHitChiSquareSPH.at(iCounterAddress)->Fill(dHitChiSquareSPH);
        fCorrectedHitChiSquareSPH.at(iCounterAddress)->Fill(dCorrectedHitChiSquareSPH);
      }

      if(fdChiSquareLimit >= dHitChiSquare)
      {
        fChiSquareEfficiency.at(iCounterAddress)->Fill(kTRUE, fmod(dPointTime/1e9, fdTimeAxisRange));
      }
      else
      {
        fChiSquareEfficiency.at(iCounterAddress)->Fill(kFALSE, fmod(dPointTime/1e9, fdTimeAxisRange));
      }

      if(fdChiSquareLimit >= dCorrectedHitChiSquare)
      {
        fCorrectedChiSquareEfficiency.at(iCounterAddress)->Fill(kTRUE, fmod(dPointTime/1e9, fdTimeAxisRange));
      }
      else
      {
        fCorrectedChiSquareEfficiency.at(iCounterAddress)->Fill(kFALSE, fmod(dPointTime/1e9, fdTimeAxisRange));
      }

      fHitTimeResidual.at(iCounterAddress)->Fill(dHitTime - dPointTime);

      if(bSinglePointHit)
      {
        fHitTimeResidualSPH.at(iCounterAddress)->Fill(dHitTime - dPointTime);
      }

      if(fbFillLocalHistos)
      {
        // Scan - if necessary - the entire "TofPointTB" branch in reverse order to find
        // the closest efficient point backwards in time.
        Int_t iIndexPrevPoint(-1);
        Int_t iTimeSlicePrevPoint(-1);
        Int_t iIndexPrevEffPoint(-1);
        Int_t iTimeSlicePrevEffPoint(-1);

        Double_t dGlobalPrevPointCoordinates[3] = {0., 0., 0.};
        Double_t dLocalPrevPointCoordinates[3] = {0., 0., 0.};

        // Start form the current time slice.
        for(Int_t iPrevPoint = iPoint - 1; iPrevPoint >= 0; iPrevPoint--)
        {
          CbmTofPoint* tPrevPoint = dynamic_cast<CbmTofPoint*>(fTofPointsTB->At(iPrevPoint));
          Int_t iPrevCounterAddress = tPrevPoint->GetDetectorID();

          if(iCounterAddress == iPrevCounterAddress)
          {
            dGlobalPrevPointCoordinates[0] = tPrevPoint->GetX();
            dGlobalPrevPointCoordinates[1] = tPrevPoint->GetY();
            dGlobalPrevPointCoordinates[2] = tPrevPoint->GetZ();
            tCounter->MasterToLocal(dGlobalPrevPointCoordinates, dLocalPrevPointCoordinates);

            if(TMath::Power(fdLocalRadius, 2.) > TMath::Power(dLocalPointCoordinates[0] - dLocalPrevPointCoordinates[0], 2.)
                                                +TMath::Power(dLocalPointCoordinates[1] - dLocalPrevPointCoordinates[1], 2.))
            {
              if(-1 == iIndexPrevPoint)
              {
                iIndexPrevPoint = iPrevPoint;
                iTimeSlicePrevPoint = iInputTSIndex;

                fDistanceToLocalLastPoint.at(iCounterAddress)->Fill(dPointTime - (tPrevPoint->GetTime() - fdDataTimeOffset));
              }

              CbmMatch* tPrevPointMatch = dynamic_cast<CbmMatch*>(fTofPointMatches->At(iPrevPoint));
              if(0 < tPrevPointMatch->GetNofLinks())
              {
                iIndexPrevEffPoint = iPrevPoint;
                break;
              }
            }
          }
        }

        // Found previous efficient point in the current time slice.
        if(-1 != iIndexPrevEffPoint)
        {
          CbmTofPoint* tPrevEffPoint = dynamic_cast<CbmTofPoint*>(fTofPointsTB->At(iIndexPrevEffPoint));
          Double_t dTimeDifference = dPointTime - (tPrevEffPoint->GetTime() - fdDataTimeOffset);

          // No inefficient point between two efficient ones.
          if(iIndexPrevPoint == iIndexPrevEffPoint)
          {
            fDistanceToLastHitUnbiased.at(iCounterAddress)->Fill(dTimeDifference);
            fEfficiencyToLastHitUnbiased.at(iCounterAddress)->Fill(kTRUE, dTimeDifference);
            fClusterSizeToLastHitUnbiased.at(iCounterAddress)->Fill(dTimeDifference, iClusterSize);
            fTimeResidualToLastHitUnbiased.at(iCounterAddress)->Fill(dTimeDifference, dHitTime - dPointTime);
          }

          // Possibly some inefficient points between two efficient ones.
          fDistanceToLastHitBiased.at(iCounterAddress)->Fill(dTimeDifference);
          fEfficiencyToLastHitBiased.at(iCounterAddress)->Fill(kTRUE, dTimeDifference);
          fClusterSizeToLastHitBiased.at(iCounterAddress)->Fill(dTimeDifference, iClusterSize);
          fTimeResidualToLastHitBiased.at(iCounterAddress)->Fill(dTimeDifference, dHitTime - dPointTime);

        }
        // Need to scan earlier time slices for the previous efficient point.
        else
        {
          TChain* tChain = fFileSource->GetInChain();

          for(Int_t iTimeSlice = iInputTSIndex - 1; iTimeSlice >= 0; iTimeSlice--)
          {
            // Load only the "TofPointTB" branch entry into memory!
            tChain->SetBranchStatus("*", kFALSE);
            tChain->SetBranchStatus("TofPointTB.fTime", kTRUE);
            tChain->SetBranchStatus("TofPointTB.fDetectorID", kTRUE);
            tChain->SetBranchStatus("TofPointTB.fX", kTRUE);
            tChain->SetBranchStatus("TofPointTB.fY", kTRUE);
            tChain->GetEntry(iTimeSlice);

            for(Int_t iPrevPoint = fTofPointsTB->GetEntriesFast() - 1; iPrevPoint >= 0; iPrevPoint--)
            {
              CbmTofPoint* tPrevPoint = dynamic_cast<CbmTofPoint*>(fTofPointsTB->At(iPrevPoint));
              Int_t iPrevCounterAddress = tPrevPoint->GetDetectorID();

              if(iCounterAddress == iPrevCounterAddress)
              {
                dGlobalPrevPointCoordinates[0] = tPrevPoint->GetX();
                dGlobalPrevPointCoordinates[1] = tPrevPoint->GetY();
                dGlobalPrevPointCoordinates[2] = tPrevPoint->GetZ();
                tCounter->MasterToLocal(dGlobalPrevPointCoordinates, dLocalPrevPointCoordinates);

                if(TMath::Power(fdLocalRadius, 2.) > TMath::Power(dLocalPointCoordinates[0] - dLocalPrevPointCoordinates[0], 2.)
                                                    +TMath::Power(dLocalPointCoordinates[1] - dLocalPrevPointCoordinates[1], 2.))
                {
                  if(-1 == iIndexPrevPoint)
                  {
                    iIndexPrevPoint = iPrevPoint;
                    iTimeSlicePrevPoint = iTimeSlice;

                    fDistanceToLocalLastPoint.at(iCounterAddress)->Fill(dPointTime - (tPrevPoint->GetTime() - fdDataTimeOffset));
                  }

                  if(fPointLinkNbSet.count(std::make_pair(iTimeSlice, iPrevPoint)))
                  {
                    iIndexPrevEffPoint = iPrevPoint;
                    iTimeSlicePrevEffPoint = iTimeSlice;
                    break;
                  }
                }
              }
            }

            if(-1 != iIndexPrevEffPoint)
            {
              break;
            }
          }

          // Found previous efficient point.
          if(-1 != iIndexPrevEffPoint)
          {
            CbmTofPoint* tPrevEffPoint = dynamic_cast<CbmTofPoint*>(fTofPointsTB->At(iIndexPrevEffPoint));
            Double_t dTimeDifference = dPointTime - (tPrevEffPoint->GetTime() - fdDataTimeOffset);

            if(iIndexPrevPoint == iIndexPrevEffPoint)
            {
              // No inefficient point between two efficient ones.
              if(iTimeSlicePrevPoint == iTimeSlicePrevEffPoint)
              {
                fDistanceToLastHitUnbiased.at(iCounterAddress)->Fill(dTimeDifference);
                fEfficiencyToLastHitUnbiased.at(iCounterAddress)->Fill(kTRUE, dTimeDifference);
                fClusterSizeToLastHitUnbiased.at(iCounterAddress)->Fill(dTimeDifference, iClusterSize);
                fTimeResidualToLastHitUnbiased.at(iCounterAddress)->Fill(dTimeDifference, dHitTime - dPointTime);
              }
            }

            // Possibly some inefficient points between two efficient ones.
            fDistanceToLastHitBiased.at(iCounterAddress)->Fill(dTimeDifference);
            fEfficiencyToLastHitBiased.at(iCounterAddress)->Fill(kTRUE, dTimeDifference);
            fClusterSizeToLastHitBiased.at(iCounterAddress)->Fill(dTimeDifference, iClusterSize);
            fTimeResidualToLastHitBiased.at(iCounterAddress)->Fill(dTimeDifference, dHitTime - dPointTime);

          }

          // Reload the current "TofPointTB" branch entry into memory.
          if(0 < iInputTSIndex)
          {
            tChain->GetEntry(iInputTSIndex);
            tChain->SetBranchStatus("*", kTRUE);
          }

        }

        fPointLinkNbSet.emplace(std::make_pair(iInputTSIndex, iPoint));

      }
    }

    if(1 < iNHitsPerPoint)
    {
      // MC point signals close in time might have been merged in CbmTofSignalBuffer.
      // LOG(INFO)<<Form("Primary MC point %d information appears in more than one (%d) ideal ToF hit.", iPoint, iNHitsPerPoint)<<FairLogger::endl;

      fiTotalNMergedPoints++;
    }

    fiTotalNPoints++;
  }


  CbmTofDetectorInfo* tCounterInfo = new CbmTofDetectorInfo(kTof, 0, 0, 0, 0, 0);

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

    tCounterInfo->fSMtype = tDigi->GetType();
    tCounterInfo->fSModule = tDigi->GetSm();
    tCounterInfo->fCounter = tDigi->GetRpc();

    Int_t iCounterAddress = CbmTofWall::Instance()->GetGeoHandler()->GetDetIdPointer()->SetDetectorInfo(*tCounterInfo);

    Int_t iNMCLinks = tDigi->GetMatch()->GetNofLinks();
    Int_t iLink(0);

    for(; iLink < iNMCLinks; iLink++)
    {
      const CbmLink& tLink = tDigi->GetMatch()->GetLink(iLink);
      Int_t iPointIndex = tLink.GetIndex();

      // MC reference point is not available.
      if(-1 == iPointIndex)
      {
        continue;
      }

      CbmTofPoint* tTofPoint = dynamic_cast<CbmTofPoint*>(fTofPointsTB->At(iPointIndex));
      CbmMatch* tPointMatch = dynamic_cast<CbmMatch*>(fTofPointMatches->At(iPointIndex));

      // MC clustering paradigm: If the MC point that created a digi created a hit, then
      // the digi contributes to the hit.
      if(0 < tPointMatch->GetNofLinks())
      {
        fChannelDigiRates.at(iCounterAddress)->Fill(kTRUE, fmod((tDigi->GetTime() - fdDataTimeOffset)/1e9, fdTimeAxisRange), 2*tDigi->GetChannel() + tDigi->GetSide());
        break;
      }
    }

    // The digi does not contribute to any hit.
    if(iLink == iNMCLinks)
    {
      fChannelDigiRates.at(iCounterAddress)->Fill(kFALSE, fmod((tDigi->GetTime() - fdDataTimeOffset)/1e9, fdTimeAxisRange), 2*tDigi->GetChannel() + tDigi->GetSide());
    }
  }

  delete tCounterInfo;
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
InitStatus CbmTofRecoMCQa::Init()
{
  if(fbConsiderPointClusters && fbConsiderRefTracks)
  {
    LOG(FATAL)<<"Cannot request point cluster and reference track selection simultaneously."<<FairLogger::endl;
  }

  if(!(CbmTofWall::Instance()->Initialize(fDigiTbParSet)))
  {
    return kFATAL;
  }

	// --- Get FairRootManager instance
  FairRootManager* tRootManager = FairRootManager::Instance();
  if(!tRootManager)
  {
    LOG(ERROR)<<"FairRootManager not found."<<FairLogger::endl;
    return kFATAL;
  }

  fFileSource = dynamic_cast<FairFileSource*>(tRootManager->GetSource());
  if(!fFileSource)
  {
    LOG(ERROR)<<"Could not get FairFileSource from FairRootManager."<<FairLogger::endl;
    return kFATAL;
  }

  TString tPointBranchName;
  TString tDigiBranchName;
  TString tHitBranchName;
  TString tHitDigiMatchBranchName;
  TString tHitPointMatchBranchName;
  TString tPointHitMatchBranchName;

  if(fbAlternativeBranchNames)
  {
    tPointBranchName         = "ATofPointTB"; // TODO: take into account if 'fbFillLocalHistos' is true
    tDigiBranchName          = "ATofDigiExp";
    tHitBranchName           = "ATofHit";
    tHitDigiMatchBranchName  = "ATofDigiMatch";
    tHitPointMatchBranchName = "ATofHitMatch";
    tPointHitMatchBranchName = "ATofPointMatch";
  }
  else
  {
    tPointBranchName         = "TofPointTB";
    tDigiBranchName          = "TofDigiExp";
    tHitBranchName           = "TofHit";
    tHitDigiMatchBranchName  = "TofDigiMatch";
    tHitPointMatchBranchName = "TofHitMatch";
    tPointHitMatchBranchName = "TofPointMatch";
  }

  // --- Get input array (CbmTofPoint)
  fTofPointsTB = dynamic_cast<TClonesArray*>(tRootManager->GetObject(tPointBranchName));
  if(!fTofPointsTB)
  {
    LOG(ERROR)<<Form("Could not retrieve branch %s from FairRootManager.", tPointBranchName.Data())<<FairLogger::endl;
    return kFATAL;
  }

  // --- Get input array (CbmTofDigiExp)
  // Ensure compatibility with realistically calibrated digis created by 'CbmTofTestBeamClusterizer'.
  fTofDigis = dynamic_cast<TClonesArray*>(tRootManager->GetObject("TofCalDigi"));
  if(!fTofDigis)
  {
    fTofDigis = dynamic_cast<TClonesArray*>(tRootManager->GetObject(tDigiBranchName));
    if(!fTofDigis)
    {
      LOG(ERROR)<<Form("Could not retrieve branch %s from FairRootManager.", tDigiBranchName.Data())<<FairLogger::endl;
      return kFATAL;
    }
  }

  // --- Get input array (CbmTofHit)
  fTofHits = dynamic_cast<TClonesArray*>(tRootManager->GetObject(tHitBranchName));
  if(!fTofHits)
  {
    LOG(ERROR)<<Form("Could not retrieve branch %s from FairRootManager.", tHitBranchName.Data())<<FairLogger::endl;
    return kFATAL;
  }

  // --- Get input array (CbmMatch)
  fTofHitDigiMatches = dynamic_cast<TClonesArray*>(tRootManager->GetObject(tHitDigiMatchBranchName));
  if(!fTofHitDigiMatches)
  {
    LOG(ERROR)<<Form("Could not retrieve branch %s from FairRootManager.", tHitDigiMatchBranchName.Data())<<FairLogger::endl;
    return kFATAL;
  }

  // --- Get input array (CbmMatch)
  fTofHitPointMatches = dynamic_cast<TClonesArray*>(tRootManager->GetObject(tHitPointMatchBranchName));
  if(!fTofHitPointMatches)
  {
    LOG(ERROR)<<Form("Could not retrieve branch %s from FairRootManager.", tHitPointMatchBranchName.Data())<<FairLogger::endl;
    return kFATAL;
  }

  fTofPointMatches = dynamic_cast<TClonesArray*>(tRootManager->GetObject(tPointHitMatchBranchName));
  if(!fTofPointMatches)
  {
    LOG(ERROR)<<Form("Could not retrieve branch %s from FairRootManager.", tPointHitMatchBranchName.Data())<<FairLogger::endl;
    return kFATAL;
  }


  if(fbAverageSpills)
  {
    fdTimeAxisRange = fdSpillLength + fdSpillBreakLength;
  }
  else
  {
    fdTimeAxisRange = fiNSerialSpills*(fdSpillLength + fdSpillBreakLength);
  }

  fiNTimeAxisBins = static_cast<Int_t>(fdTimeAxisRange/fdTimeAxisBinWidth);

  CreateSigmaHistograms();
  CreateHistograms();

  for(auto const & itCounter : CbmTofWall::Instance()->GetCounterMap())
  {
    Int_t iCounterAddress = itCounter.first;
    fdPreviousPointTime.emplace(iCounterAddress, 0.);
  }

  return kSUCCESS;
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
void CbmTofRecoMCQa::SetParContainers()
{
  FairRun* tRun = FairRun::Instance();
  if(!tRun)
  {
    LOG(FATAL)<<Form("FairRun instance not found.")<<FairLogger::endl;
  }

  FairRuntimeDb* tDataBase = tRun->GetRuntimeDb();
  if(!tDataBase)
  {
    LOG(FATAL)<<Form("FairRuntimeDb instance not found.")<<FairLogger::endl;
  }

  // request initialization of parameter containers
  // The actual parameter container initialization happens only upon calling
  // FairRuntimeDb::initContainers in the FairRun instance
  fDigiTbParSet = dynamic_cast<CbmTofDigiTbPar*>(tDataBase->getContainer("CbmTofDigiTbPar"));
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
void CbmTofRecoMCQa::Finish()
{
  // Downscale the histograms displaying rate according to the number of serial spills.
  if(fbAverageSpills)
  {
    for(auto const & itCounter : CbmTofWall::Instance()->GetCounterMap())
    {
      Int_t iCounterAddress = itCounter.first;

      // Should be safe.
      const_cast<TH1*>(fLocalSurfaceFlux.at(iCounterAddress)->GetPassedHistogram())->Scale(1./fiNSerialSpills);
      const_cast<TH1*>(fLocalSurfaceFlux.at(iCounterAddress)->GetTotalHistogram())->Scale(1./fiNSerialSpills);
      const_cast<TH1*>(fChannelDigiRates.at(iCounterAddress)->GetPassedHistogram())->Scale(1./fiNSerialSpills);
      const_cast<TH1*>(fChannelDigiRates.at(iCounterAddress)->GetTotalHistogram())->Scale(1./fiNSerialSpills);
    }
  }

  WriteHistograms();

  LOG(INFO)<<FairLogger::endl;
  LOG(INFO)<<Form("Total number of MC points scanned: %d", fiTotalNPoints)<<FairLogger::endl;
  LOG(INFO)<<Form("Number of merged MC points skipped: %d", fiTotalNMergedPoints)<<FairLogger::endl;
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
void CbmTofRecoMCQa::CreateSigmaHistograms()
{
  TDirectory* tOldDirectory = gDirectory;
  gROOT->cd();

  for(auto const & itCounter : CbmTofWall::Instance()->GetCounterMap())
  {
    Int_t iCounterAddress = itCounter.first;
    CbmTofCounter* tCounter = itCounter.second;

    Int_t iModuleType = tCounter->GetModuleType();
    Int_t iModuleIndex = tCounter->GetModuleIndex();
    Int_t iCounterIndex = tCounter->GetCounterIndex();

    fAcrossPositionResidualEvolutionSigma.emplace(iCounterAddress, new TH1F(Form("tAcrossPositionResidualEvolutionSigma_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                            "; MC point time [s]; across residual sigma [cm]",
                                                                            fiNTimeAxisBins, 0., fdTimeAxisRange));

    fAlongPositionResidualEvolutionSigma.emplace(iCounterAddress, new TH1F(Form("tAlongPositionResidualEvolutionSigma_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                           "; MC point time [s]; along residual sigma [cm]",
                                                                           fiNTimeAxisBins, 0., fdTimeAxisRange));

    fTimeResidualEvolutionSigma.emplace(iCounterAddress, new TH1F(Form("tTimeResidualEvolutionSigma_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                  "; MC point time [s]; time residual sigma [ns]",
                                                                  fiNTimeAxisBins, 0., fdTimeAxisRange));

    fAcrossPositionResidualEvolutionSigmaSPH.emplace(iCounterAddress, new TH1F(Form("tAcrossPositionResidualEvolutionSigmaSPH_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                               "; MC point time [s]; across residual sigma [cm]",
                                                                               fiNTimeAxisBins, 0., fdTimeAxisRange));

    fAlongPositionResidualEvolutionSigmaSPH.emplace(iCounterAddress, new TH1F(Form("tAlongPositionResidualEvolutionSigmaSPH_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                              "; MC point time [s]; along residual sigma [cm]",
                                                                              fiNTimeAxisBins, 0., fdTimeAxisRange));

    fTimeResidualEvolutionSigmaSPH.emplace(iCounterAddress, new TH1F(Form("tTimeResidualEvolutionSigmaSPH_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                     "; MC point time [s]; time residual sigma [ns]",
                                                                     fiNTimeAxisBins, 0., fdTimeAxisRange));

    for(Int_t iBin = 1; iBin <= fiNTimeAxisBins; iBin++)
    {
      (fAcrossPositionResidualEvolutionSigma.at(iCounterAddress))->SetBinContent(iBin, 1.);
      (fAlongPositionResidualEvolutionSigma.at(iCounterAddress))->SetBinContent(iBin, 1.);
      (fTimeResidualEvolutionSigma.at(iCounterAddress))->SetBinContent(iBin, 1.);

      (fAcrossPositionResidualEvolutionSigmaSPH.at(iCounterAddress))->SetBinContent(iBin, 1.);
      (fAlongPositionResidualEvolutionSigmaSPH.at(iCounterAddress))->SetBinContent(iBin, 1.);
      (fTimeResidualEvolutionSigmaSPH.at(iCounterAddress))->SetBinContent(iBin, 1.);
    }
  }

  if(fHistogramFileName.IsNull())
  {
    LOG(FATAL)<<"No histogram file name specified!"<<FairLogger::endl;
  }

  TFile* tHistogramFile = new TFile(fHistogramFileName.Data(), "READ");
  if(tHistogramFile->IsZombie())
  {
    tHistogramFile->Close();
    delete tHistogramFile;

    LOG(INFO)<<"Could not open histogram file to extract residual widthes."<<FairLogger::endl;
  }
  else
  {
    for(auto const & itCounter : CbmTofWall::Instance()->GetCounterMap())
    {
      Int_t iCounterAddress = itCounter.first;
      CbmTofCounter* tCounter = itCounter.second;

      Int_t iModuleType = tCounter->GetModuleType();
      Int_t iModuleIndex = tCounter->GetModuleIndex();
      Int_t iCounterIndex = tCounter->GetCounterIndex();

      if(!gDirectory->cd(Form("%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex)))
      {
        LOG(FATAL)<<Form("Could not find directory '%d_%d_%d' in histogram file.", iModuleType, iModuleIndex, iCounterIndex)<<FairLogger::endl;
      }

      TH2F* tAcrossPositionResidualEvolution(NULL);
      TH2F* tAlongPositionResidualEvolution(NULL);
      TH2F* tTimeResidualEvolution(NULL);
      TH2F* tAcrossPositionResidualEvolutionSPH(NULL);
      TH2F* tAlongPositionResidualEvolutionSPH(NULL);
      TH2F* tTimeResidualEvolutionSPH(NULL);

      gDirectory->GetObject(Form("tAcrossPositionResidualEvolution_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex), tAcrossPositionResidualEvolution);
      gDirectory->GetObject(Form("tAlongPositionResidualEvolution_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex), tAlongPositionResidualEvolution);
      gDirectory->GetObject(Form("tTimeResidualEvolution_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex), tTimeResidualEvolution);
      gDirectory->GetObject(Form("tAcrossPositionResidualEvolutionSPH_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex), tAcrossPositionResidualEvolutionSPH);
      gDirectory->GetObject(Form("tAlongPositionResidualEvolutionSPH_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex), tAlongPositionResidualEvolutionSPH);
      gDirectory->GetObject(Form("tTimeResidualEvolutionSPH_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex), tTimeResidualEvolutionSPH);


      if(!tAcrossPositionResidualEvolution)
      {
        LOG(FATAL)<<Form("Could not find histogram 'tAcrossPositionResidualEvolution_%d_%d_%d' in histogram file.", iModuleType, iModuleIndex, iCounterIndex)<<FairLogger::endl;
      }

      if(!tAlongPositionResidualEvolution)
      {
        LOG(FATAL)<<Form("Could not find histogram 'tAlongPositionResidualEvolution_%d_%d_%d' in histogram file.", iModuleType, iModuleIndex, iCounterIndex)<<FairLogger::endl;
      }

      if(!tTimeResidualEvolution)
      {
        LOG(FATAL)<<Form("Could not find histogram 'tTimeResidualEvolution_%d_%d_%d' in histogram file.", iModuleType, iModuleIndex, iCounterIndex)<<FairLogger::endl;
      }

      if(!tAcrossPositionResidualEvolutionSPH)
      {
        LOG(FATAL)<<Form("Could not find histogram 'tAcrossPositionResidualEvolutionSPH_%d_%d_%d' in histogram file.", iModuleType, iModuleIndex, iCounterIndex)<<FairLogger::endl;
      }

      if(!tAlongPositionResidualEvolutionSPH)
      {
        LOG(FATAL)<<Form("Could not find histogram 'tAlongPositionResidualEvolutionSPH_%d_%d_%d' in histogram file.", iModuleType, iModuleIndex, iCounterIndex)<<FairLogger::endl;
      }

      if(!tTimeResidualEvolutionSPH)
      {
        LOG(FATAL)<<Form("Could not find histogram 'tTimeResidualEvolutionSPH_%d_%d_%d' in histogram file.", iModuleType, iModuleIndex, iCounterIndex)<<FairLogger::endl;
      }


      TH1D* tProjection(NULL);

      for(Int_t iBin = 1; iBin <= fiNTimeAxisBins; iBin++)
      {
        tProjection = tAcrossPositionResidualEvolution->ProjectionY("_py", iBin, iBin);

        if(tProjection->GetEntries())
        {
          (fAcrossPositionResidualEvolutionSigma.at(iCounterAddress))->SetBinContent(iBin, tProjection->GetRMS());
          (fAcrossPositionResidualEvolutionSigma.at(iCounterAddress))->SetBinError(iBin, tProjection->GetRMSError());
        }

        delete tProjection;

        tProjection = tAlongPositionResidualEvolution->ProjectionY("_py", iBin, iBin);

        if(tProjection->GetEntries())
        {
          if(!static_cast<Int_t>(tProjection->Fit("gaus","Q0")))
          {
            (fAlongPositionResidualEvolutionSigma.at(iCounterAddress))->SetBinContent(iBin, (tProjection->GetFunction("gaus"))->GetParameter("Sigma"));
            (fAlongPositionResidualEvolutionSigma.at(iCounterAddress))->SetBinError(iBin, (tProjection->GetFunction("gaus"))->GetParError(2));
          }
        }

        delete tProjection;

        tProjection = tTimeResidualEvolution->ProjectionY("_py", iBin, iBin);

        if(tProjection->GetEntries())
        {
          if(!static_cast<Int_t>(tProjection->Fit("gaus", "Q0")))
          {
            (fTimeResidualEvolutionSigma.at(iCounterAddress))->SetBinContent(iBin, (tProjection->GetFunction("gaus"))->GetParameter("Sigma"));
            (fTimeResidualEvolutionSigma.at(iCounterAddress))->SetBinError(iBin, (tProjection->GetFunction("gaus"))->GetParError(2));
          }
        }

        delete tProjection;


        tProjection = tAcrossPositionResidualEvolutionSPH->ProjectionY("_py", iBin, iBin);

        if(tProjection->GetEntries())
        {
          (fAcrossPositionResidualEvolutionSigmaSPH.at(iCounterAddress))->SetBinContent(iBin, tProjection->GetRMS());
          (fAcrossPositionResidualEvolutionSigmaSPH.at(iCounterAddress))->SetBinError(iBin, tProjection->GetRMSError());
        }

        delete tProjection;

        tProjection = tAlongPositionResidualEvolutionSPH->ProjectionY("_py", iBin, iBin);

        if(tProjection->GetEntries())
        {
          if(!static_cast<Int_t>(tProjection->Fit("gaus","Q0")))
          {
            (fAlongPositionResidualEvolutionSigmaSPH.at(iCounterAddress))->SetBinContent(iBin, (tProjection->GetFunction("gaus"))->GetParameter("Sigma"));
            (fAlongPositionResidualEvolutionSigmaSPH.at(iCounterAddress))->SetBinError(iBin, (tProjection->GetFunction("gaus"))->GetParError(2));
          }
        }

        delete tProjection;

        tProjection = tTimeResidualEvolutionSPH->ProjectionY("_py", iBin, iBin);

        if(tProjection->GetEntries())
        {
          if(!static_cast<Int_t>(tProjection->Fit("gaus", "Q0")))
          {
            (fTimeResidualEvolutionSigmaSPH.at(iCounterAddress))->SetBinContent(iBin, (tProjection->GetFunction("gaus"))->GetParameter("Sigma"));
            (fTimeResidualEvolutionSigmaSPH.at(iCounterAddress))->SetBinError(iBin, (tProjection->GetFunction("gaus"))->GetParError(2));
          }
        }

        delete tProjection;
      }

      gDirectory->cd("../");
    }

    tHistogramFile->Close();
    delete tHistogramFile;
  }

  tOldDirectory->cd();
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
Double_t CbmTofRecoMCQa::GetAcrossPositionResidualSigma(Int_t iCounterAddress, Double_t dPointTime)
{
  return (fAcrossPositionResidualEvolutionSigma.at(iCounterAddress))->GetBinContent((fAcrossPositionResidualEvolutionSigma.at(iCounterAddress))->FindBin(dPointTime));
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
Double_t CbmTofRecoMCQa::GetAlongPositionResidualSigma(Int_t iCounterAddress, Double_t dPointTime)
{
  return (fAlongPositionResidualEvolutionSigma.at(iCounterAddress))->GetBinContent((fAlongPositionResidualEvolutionSigma.at(iCounterAddress))->FindBin(dPointTime));
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
Double_t CbmTofRecoMCQa::GetTimeResidualSigma(Int_t iCounterAddress, Double_t dPointTime)
{
  return (fTimeResidualEvolutionSigma.at(iCounterAddress))->GetBinContent((fTimeResidualEvolutionSigma.at(iCounterAddress))->FindBin(dPointTime));
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
Double_t CbmTofRecoMCQa::GetAcrossPositionResidualSigmaSPH(Int_t iCounterAddress, Double_t dPointTime)
{
  return (fAcrossPositionResidualEvolutionSigmaSPH.at(iCounterAddress))->GetBinContent((fAcrossPositionResidualEvolutionSigmaSPH.at(iCounterAddress))->FindBin(dPointTime));
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
Double_t CbmTofRecoMCQa::GetAlongPositionResidualSigmaSPH(Int_t iCounterAddress, Double_t dPointTime)
{
  return (fAlongPositionResidualEvolutionSigmaSPH.at(iCounterAddress))->GetBinContent((fAlongPositionResidualEvolutionSigmaSPH.at(iCounterAddress))->FindBin(dPointTime));
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
Double_t CbmTofRecoMCQa::GetTimeResidualSigmaSPH(Int_t iCounterAddress, Double_t dPointTime)
{
  return (fTimeResidualEvolutionSigmaSPH.at(iCounterAddress))->GetBinContent((fTimeResidualEvolutionSigmaSPH.at(iCounterAddress))->FindBin(dPointTime));
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
void CbmTofRecoMCQa::CreateHistograms()
{
  TDirectory* tOldDirectory = gDirectory;

  if(fHistogramFileName.IsNull())
  {
    LOG(FATAL)<<"No histogram file name specified!"<<FairLogger::endl;
  }

  fHistogramFile = new TFile(fHistogramFileName.Data(),"RECREATE");

  for(auto const & itCounter : CbmTofWall::Instance()->GetCounterMap())
  {
    Int_t iCounterAddress = itCounter.first;
    CbmTofCounter* tCounter = itCounter.second;

    Int_t iModuleType = tCounter->GetModuleType();
    Int_t iModuleIndex = tCounter->GetModuleIndex();
    Int_t iCounterIndex = tCounter->GetCounterIndex();

    Int_t iNumberReadoutElectrodes = tCounter->GetNumberReadoutElectrodes();
    Int_t iNWalkBinsX = tCounter->GetNWalkBinsX();
    Double_t dMinimumToT = tCounter->GetMinimumToT();
    Double_t dMaximumToT = tCounter->GetMaximumToT();
    Double_t dCellAlongWidth = tCounter->GetCellAlongWidth();
    Double_t dCellAcrossWidth = tCounter->GetCellAcrossWidth();
    Bool_t bCellsAlongX = tCounter->GetCellsAlongX();

    TDirectory* tCounterDirectory = gDirectory->mkdir(Form("%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex));
    tCounterDirectory->cd();

    fClusterSizes.emplace(iCounterAddress, new TH2F(Form("tClusterSizes_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                    "; hit main strip []; cluster size [strips]",
                                                    iNumberReadoutElectrodes, 0, iNumberReadoutElectrodes,
                                                    iNumberReadoutElectrodes, 0.5, iNumberReadoutElectrodes + 0.5));
//                                                    iNumberReadoutElectrodes/2, 0.5, iNumberReadoutElectrodes/2+0.5));

    fStripToTSpectra.emplace(iCounterAddress, new TH2F(Form("tStripToTSpectra_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                       "; strip side [2*index]; strip side ToT [ns]",
                                                       2*iNumberReadoutElectrodes, 0, 2*iNumberReadoutElectrodes,
                                                       iNWalkBinsX, dMinimumToT, dMaximumToT));

    fStripPositionSpectra.emplace(iCounterAddress, new TH2F(Form("tStripPositionSpectra_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                            "; strip [index]; strip along position [cm]",
                                                            iNumberReadoutElectrodes, 0, iNumberReadoutElectrodes,
                                                            20*TMath::CeilNint(dCellAlongWidth), -TMath::Ceil(dCellAlongWidth), TMath::Ceil(dCellAlongWidth)));

    fPointPositionSpectra.emplace(iCounterAddress, new TH2F(Form("tPointPositionSpectra_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                            "; point across position [cm]; point along position [cm]",
                                                            2*TMath::CeilNint(iNumberReadoutElectrodes/2.*dCellAcrossWidth),
                                                            -TMath::Ceil(iNumberReadoutElectrodes/2.*dCellAcrossWidth),
                                                            TMath::Ceil(iNumberReadoutElectrodes/2.*dCellAcrossWidth),
                                                            20*TMath::CeilNint(dCellAlongWidth), -TMath::Ceil(dCellAlongWidth), TMath::Ceil(dCellAlongWidth)));

    fHitPositionSpectra.emplace(iCounterAddress, new TH2F(Form("tHitPositionSpectra_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                          "; hit across position [cm]; hit along position [cm]",
                                                          2*TMath::CeilNint(iNumberReadoutElectrodes/2.*dCellAcrossWidth),
                                                          -TMath::Ceil(iNumberReadoutElectrodes/2.*dCellAcrossWidth),
                                                          TMath::Ceil(iNumberReadoutElectrodes/2.*dCellAcrossWidth),
                                                          20*TMath::CeilNint(dCellAlongWidth), -TMath::Ceil(dCellAlongWidth), TMath::Ceil(dCellAlongWidth)));

    fStripTimeResidualWalk.emplace(iCounterAddress, new TH3F(Form("tStripTimeResidualWalk_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                             "; strip side [2*index]; strip side ToT [ns]; strip time residual [ns]",
                                                             2*iNumberReadoutElectrodes, 0, 2*iNumberReadoutElectrodes,
                                                             iNWalkBinsX, dMinimumToT, dMaximumToT,
                                                             811, -2., 2.));

    fHitAcrossPositionResidual.emplace(iCounterAddress, new TH1F(Form("tHitAcrossPositionResidual_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                 "; hit across residual [cm]",
                                                                 1001, -5., 5.));

    fHitAlongPositionResidual.emplace(iCounterAddress, new TH1F(Form("tHitAlongPositionResidual_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                "; hit along residual [cm]",
                                                                3001, -15., 15.));

    fHitTimeResidual.emplace(iCounterAddress, new TH1F(Form("tHitTimeResidual_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                       "; hit time residual [ns]",
                                                       3001, -1.5, 1.5));

    fCentralClusterSizeEvolution.emplace(iCounterAddress, new TH2F(Form("tCentralClusterSizeEvolution_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                   "; MC point time [s]; cluster size [strips]",
                                                                   fiNTimeAxisBins, 0., fdTimeAxisRange,
                                                                   iNumberReadoutElectrodes, 0.5, iNumberReadoutElectrodes + 0.5));
//                                                                   iNumberReadoutElectrodes/2, 0.5, iNumberReadoutElectrodes/2+0.5));

    fCentralTimeResidualEvolution.emplace(iCounterAddress, new TH2F(Form("tCentralTimeResidualEvolution_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                    "; MC point time [s]; hit time residual [ns]",
                                                                    fiNTimeAxisBins, 0., fdTimeAxisRange,
                                                                    3001, -1.5, 1.5));
//                                                                    2001, -1., 1.));

    fHitChiSquare.emplace(iCounterAddress, new TH1F(Form("tHitChiSquare_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                    "; hit #chi^{2}_{3} []",
                                                    300, 0., 30.));

    fCorrectedHitChiSquare.emplace(iCounterAddress, new TH1F(Form("tCorrectedHitChiSquare_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                             "; hit #chi^{2}_{3} []",
                                                             300, 0., 30.));

    fClusterSizeEvolution.emplace(iCounterAddress, new TH2F(Form("tClusterSizeEvolution_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                            "; MC point time [s]; cluster size [strips]",
                                                            fiNTimeAxisBins, 0., fdTimeAxisRange,
                                                            iNumberReadoutElectrodes, 0.5, iNumberReadoutElectrodes + 0.5));
//                                                            iNumberReadoutElectrodes/2, 0.5, iNumberReadoutElectrodes/2+0.5));

    fAcrossPositionResidualEvolution.emplace(iCounterAddress, new TH2F(Form("tAcrossPositionResidualEvolution_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                       "; MC point time [s]; hit across residual [cm]",
                                                                       fiNTimeAxisBins, 0., fdTimeAxisRange,
                                                                       1001, -5., 5.));

    fAlongPositionResidualEvolution.emplace(iCounterAddress, new TH2F(Form("tAlongPositionResidualEvolution_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                      "; MC point time [s]; hit along residual [cm]",
                                                                      fiNTimeAxisBins, 0., fdTimeAxisRange,
                                                                      3001, -15., 15.));

    fTimeResidualEvolution.emplace(iCounterAddress, new TH2F(Form("tTimeResidualEvolution_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                             "; MC point time [s]; hit time residual [ns]",
                                                             fiNTimeAxisBins, 0., fdTimeAxisRange,
                                                             3001, -1.5, 1.5));
//                                                             2001, -1., 1.));

    fCentralStripToTEvolution.emplace(iCounterAddress, new TH2F(Form("tCentralStripToTEvolution_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                "; MC point time [s]; strip side ToT [ns]",
                                                                fiNTimeAxisBins, 0., fdTimeAxisRange,
                                                                iNWalkBinsX, dMinimumToT, dMaximumToT));

    fIntegratedEfficiency.emplace(iCounterAddress, new TEfficiency(Form("tIntegratedEfficiency_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                   "; MC point time [s]; counter efficiency []",
                                                                   fiNTimeAxisBins, 0., fdTimeAxisRange));

    fChiSquareEfficiency.emplace(iCounterAddress, new TEfficiency(Form("tChiSquareEfficiency_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                  "; MC point time [s]; #chi^{2}_{3} efficiency []",
                                                                  fiNTimeAxisBins, 0., fdTimeAxisRange));

    fCorrectedChiSquareEfficiency.emplace(iCounterAddress, new TEfficiency(Form("tCorrectedChiSquareEfficiency_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                           "; MC point time [s]; #chi^{2}_{3} efficiency []",
                                                                           fiNTimeAxisBins, 0., fdTimeAxisRange));


    // QA histograms for MC points referring to a hit which in turn refers to this MC point only (SPH := single-point hit)
    fClusterSizesSPH.emplace(iCounterAddress, new TH2F(Form("tClusterSizesSPH_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                       "; hit main strip []; cluster size [strips]",
                                                       iNumberReadoutElectrodes, 0, iNumberReadoutElectrodes,
                                                       iNumberReadoutElectrodes, 0.5, iNumberReadoutElectrodes + 0.5));
//                                                       iNumberReadoutElectrodes/2, 0.5, iNumberReadoutElectrodes/2+0.5));

    fStripToTSpectraSPH.emplace(iCounterAddress, new TH2F(Form("tStripToTSpectraSPH_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                          "; strip side [2*index]; strip side ToT [ns]",
                                                          2*iNumberReadoutElectrodes, 0, 2*iNumberReadoutElectrodes,
                                                          iNWalkBinsX, dMinimumToT, dMaximumToT));

    fStripPositionSpectraSPH.emplace(iCounterAddress, new TH2F(Form("tStripPositionSpectraSPH_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                               "; strip [index]; strip along position [cm]",
                                                               iNumberReadoutElectrodes, 0, iNumberReadoutElectrodes,
                                                               20*TMath::CeilNint(dCellAlongWidth), -TMath::Ceil(dCellAlongWidth), TMath::Ceil(dCellAlongWidth)));

    fHitPositionSpectraSPH.emplace(iCounterAddress, new TH2F(Form("tHitPositionSpectraSPH_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                             "; hit across position [cm]; hit along position [cm]",
                                                             2*TMath::CeilNint(iNumberReadoutElectrodes/2.*dCellAcrossWidth),
                                                             -TMath::Ceil(iNumberReadoutElectrodes/2.*dCellAcrossWidth),
                                                             TMath::Ceil(iNumberReadoutElectrodes/2.*dCellAcrossWidth),
                                                             20*TMath::CeilNint(dCellAlongWidth), -TMath::Ceil(dCellAlongWidth), TMath::Ceil(dCellAlongWidth)));

    fStripTimeResidualWalkSPH.emplace(iCounterAddress, new TH3F(Form("tStripTimeResidualWalkSPH_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                "; strip side [2*index]; strip side ToT [ns]; strip time residual [ns]",
                                                                2*iNumberReadoutElectrodes, 0, 2*iNumberReadoutElectrodes,
                                                                iNWalkBinsX, dMinimumToT, dMaximumToT,
                                                                811, -2., 2.));

    fHitAcrossPositionResidualSPH.emplace(iCounterAddress, new TH1F(Form("tHitAcrossPositionResidualSPH_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                    "; hit across residual [cm]",
                                                                    1001, -5., 5.));

    fHitAlongPositionResidualSPH.emplace(iCounterAddress, new TH1F(Form("tHitAlongPositionResidualSPH_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                   "; hit along residual [cm]",
                                                                   3001, -15., 15.));

    fHitTimeResidualSPH.emplace(iCounterAddress, new TH1F(Form("tHitTimeResidualSPH_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                          "; hit time residual [ns]",
                                                          3001, -1.5, 1.5));

    fHitChiSquareSPH.emplace(iCounterAddress, new TH1F(Form("tHitChiSquareSPH_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                       "; hit #chi^{2}_{3} []",
                                                       300, 0., 30.));

    fCorrectedHitChiSquareSPH.emplace(iCounterAddress, new TH1F(Form("tCorrectedHitChiSquareSPH_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                "; hit #chi^{2}_{3} []",
                                                                300, 0., 30.));

    fClusterSizeEvolutionSPH.emplace(iCounterAddress, new TH2F(Form("tClusterSizeEvolutionSPH_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                               "; MC point time [s]; cluster size [strips]",
                                                               fiNTimeAxisBins, 0., fdTimeAxisRange,
                                                               iNumberReadoutElectrodes, 0.5, iNumberReadoutElectrodes + 0.5));
//                                                               iNumberReadoutElectrodes/2, 0.5, iNumberReadoutElectrodes/2+0.5));

    fAcrossPositionResidualEvolutionSPH.emplace(iCounterAddress, new TH2F(Form("tAcrossPositionResidualEvolutionSPH_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                          "; MC point time [s]; hit across residual [cm]",
                                                                          fiNTimeAxisBins, 0., fdTimeAxisRange,
                                                                          1001, -5., 5.));

    fAlongPositionResidualEvolutionSPH.emplace(iCounterAddress, new TH2F(Form("tAlongPositionResidualEvolutionSPH_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                         "; MC point time [s]; hit along residual [cm]",
                                                                         fiNTimeAxisBins, 0., fdTimeAxisRange,
                                                                         3001, -15., 15.));

    fTimeResidualEvolutionSPH.emplace(iCounterAddress, new TH2F(Form("tTimeResidualEvolutionSPH_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                "; MC point time [s]; hit time residual [ns]",
                                                                fiNTimeAxisBins, 0., fdTimeAxisRange,
                                                                3001, -1.5, 1.5));
//                                                                2001, -1., 1.));

    fCentralStripToTEvolutionSPH.emplace(iCounterAddress, new TH2F(Form("tCentralStripToTEvolutionSPH_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                   "; MC point time [s]; strip side ToT [ns]",
                                                                   fiNTimeAxisBins, 0., fdTimeAxisRange,
                                                                   iNWalkBinsX, dMinimumToT, dMaximumToT));


    Int_t iMinDecimalPower = -9; // [1 ps]
    Int_t iMaxDecimalPower = 11; // [100 s]
    Double_t dHitDistanceBinWidths = 0.1;

    Int_t iNHitDistanceBins = (iMaxDecimalPower - iMinDecimalPower)/dHitDistanceBinWidths;
    Double_t dHitDistanceBinEdges[iNHitDistanceBins + 1];
    Double_t dLog10 = TMath::Log(10);
    for(Int_t iEdge = 0; iEdge <= iNHitDistanceBins; iEdge++)
    {
      dHitDistanceBinEdges[iEdge] = TMath::Exp(dLog10*(iEdge*dHitDistanceBinWidths + iMinDecimalPower));
    }


    fDistanceToLastPoint.emplace(iCounterAddress, new TH1F(Form("tDistanceToLastPoint_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                           "; overall TTLP [ns]",
                                                           iNHitDistanceBins, dHitDistanceBinEdges));

    fDistanceToLocalLastPoint.emplace(iCounterAddress, new TH1F(Form("tDistanceToLocalLastPoint_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                "; local TTLP [ns]",
                                                                iNHitDistanceBins, dHitDistanceBinEdges));

    fDistanceToLastHitUnbiased.emplace(iCounterAddress, new TH1F(Form("tDistanceToLastHitUnbiased_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                 "; time since last hit [ns]",
                                                                 iNHitDistanceBins, dHitDistanceBinEdges));

    fDistanceToLastHitBiased.emplace(iCounterAddress, new TH1F(Form("tDistanceToLastHitBiased_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                               "; time since last hit [ns]",
                                                               iNHitDistanceBins, dHitDistanceBinEdges));

    fEfficiencyToLastHitUnbiased.emplace(iCounterAddress, new TEfficiency(Form("tEfficiencyToLastHitUnbiased_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                          "; time since last hit [ns]; local efficiency []",
                                                                          iNHitDistanceBins, dHitDistanceBinEdges));

    fEfficiencyToLastHitBiased.emplace(iCounterAddress, new TEfficiency(Form("tEfficiencyToLastHitBiased_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                        "; time since last hit [ns]; local efficiency []",
                                                                        iNHitDistanceBins, dHitDistanceBinEdges));

    fClusterSizeToLastHitUnbiased.emplace(iCounterAddress, new TH2F(Form("tClusterSizeToLastHitUnbiased_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                    "; time since last hit [ns]; cluster size [strips]",
                                                                    iNHitDistanceBins, dHitDistanceBinEdges,
                                                                    iNumberReadoutElectrodes, 0.5, iNumberReadoutElectrodes + 0.5));
//                                                                    iNumberReadoutElectrodes/2, 0.5, iNumberReadoutElectrodes/2+0.5));

    fClusterSizeToLastHitBiased.emplace(iCounterAddress, new TH2F(Form("tClusterSizeToLastHitBiased_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                  "; time since last hit [ns]; cluster size [strips]",
                                                                  iNHitDistanceBins, dHitDistanceBinEdges,
                                                                  iNumberReadoutElectrodes, 0.5, iNumberReadoutElectrodes + 0.5));
//                                                                  iNumberReadoutElectrodes/2, 0.5, iNumberReadoutElectrodes/2+0.5));

    fTimeResidualToLastHitUnbiased.emplace(iCounterAddress, new TH2F(Form("tTimeResidualToLastHitUnbiased_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                     "; time since last hit [ns]; hit time residual [ns]",
                                                                     iNHitDistanceBins, dHitDistanceBinEdges,
                                                                     3001, -1.5, 1.5));
//                                                                     2001, -1., 1.));

    fTimeResidualToLastHitBiased.emplace(iCounterAddress, new TH2F(Form("tTimeResidualToLastHitBiased_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                                   "; time since last hit [ns]; hit time residual [ns]",
                                                                   iNHitDistanceBins, dHitDistanceBinEdges,
                                                                   3001, -1.5, 1.5));
//                                                                   2001, -1., 1.));


    fLocalSurfaceFlux.emplace(iCounterAddress, new TEfficiency(Form("tLocalSurfaceFlux_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                               "; strip across position [cm]; strip along position [cm]; MC point time [s]",
                                                               2*TMath::CeilNint(iNumberReadoutElectrodes/2.*dCellAcrossWidth),
                                                               -TMath::Ceil(iNumberReadoutElectrodes/2.*dCellAcrossWidth),
                                                               TMath::Ceil(iNumberReadoutElectrodes/2.*dCellAcrossWidth),
                                                               2*TMath::CeilNint(dCellAlongWidth), -TMath::Ceil(dCellAlongWidth), TMath::Ceil(dCellAlongWidth),
                                                               TMath::CeilNint(fdTimeAxisRange), 0., TMath::Ceil(fdTimeAxisRange)));

    fChannelDigiRates.emplace(iCounterAddress, new TEfficiency(Form("tChannelDigiRates_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                               "; digi time [s]; strip [index]; digi rate [Hz]",
                                                               fiNTimeAxisBins, 0., fdTimeAxisRange,
                                                               2*iNumberReadoutElectrodes, 0, 2*iNumberReadoutElectrodes));


    const Int_t iNMaxHitsPerPoint = 30;
    const Int_t iNMaxPointsPerHit = 100;

    fNHitsPerPoint.emplace(iCounterAddress, new TH1F(Form("tNHitsPerPoint_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                     "; #hits per point []",
                                                     iNMaxHitsPerPoint/2, 0.5, iNMaxHitsPerPoint/2+0.5));

    fNPointsPerHit.emplace(iCounterAddress, new TH1F(Form("tNPointsPerHit_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                                     "; #points per hit []",
                                                     iNMaxPointsPerHit/2, 0.5, iNMaxPointsPerHit/2+0.5));

    tCounterDirectory->cd("../");
  }

  tOldDirectory->cd();
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
void CbmTofRecoMCQa::WriteHistograms()
{
  TDirectory* tOldDirectory = gDirectory;

  fHistogramFile->cd();

  for(auto const & itCounter : CbmTofWall::Instance()->GetCounterMap())
  {
    Int_t iCounterAddress = itCounter.first;
    CbmTofCounter* tCounter = itCounter.second;

    Int_t iModuleType = tCounter->GetModuleType();
    Int_t iModuleIndex = tCounter->GetModuleIndex();
    Int_t iCounterIndex = tCounter->GetCounterIndex();

    gDirectory->cd(Form("%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex));

    (fClusterSizes.at(iCounterAddress))->Write();
    (fStripToTSpectra.at(iCounterAddress))->Write();
    (fStripPositionSpectra.at(iCounterAddress))->Write();
    (fPointPositionSpectra.at(iCounterAddress))->Write();
    (fHitPositionSpectra.at(iCounterAddress))->Write();
    (fStripTimeResidualWalk.at(iCounterAddress))->Write();
    (fHitAcrossPositionResidual.at(iCounterAddress))->Write();
    (fHitAlongPositionResidual.at(iCounterAddress))->Write();
    (fHitTimeResidual.at(iCounterAddress))->Write();
    (fCentralClusterSizeEvolution.at(iCounterAddress))->Write();
    (fCentralTimeResidualEvolution.at(iCounterAddress))->Write();
    (fHitChiSquare.at(iCounterAddress))->Write();
    (fCorrectedHitChiSquare.at(iCounterAddress))->Write();
    (fClusterSizeEvolution.at(iCounterAddress))->Write();
    (fAcrossPositionResidualEvolution.at(iCounterAddress))->Write();
    (fAlongPositionResidualEvolution.at(iCounterAddress))->Write();
    (fTimeResidualEvolution.at(iCounterAddress))->Write();
    (fCentralStripToTEvolution.at(iCounterAddress))->Write();
    (fIntegratedEfficiency.at(iCounterAddress))->Write();
    (fChiSquareEfficiency.at(iCounterAddress))->Write();
    (fCorrectedChiSquareEfficiency.at(iCounterAddress))->Write();
    (fClusterSizesSPH.at(iCounterAddress))->Write();
    (fStripToTSpectraSPH.at(iCounterAddress))->Write();
    (fStripPositionSpectraSPH.at(iCounterAddress))->Write();
    (fHitPositionSpectraSPH.at(iCounterAddress))->Write();
    (fStripTimeResidualWalkSPH.at(iCounterAddress))->Write();
    (fHitAcrossPositionResidualSPH.at(iCounterAddress))->Write();
    (fHitAlongPositionResidualSPH.at(iCounterAddress))->Write();
    (fHitTimeResidualSPH.at(iCounterAddress))->Write();
    (fHitChiSquareSPH.at(iCounterAddress))->Write();
    (fCorrectedHitChiSquareSPH.at(iCounterAddress))->Write();
    (fClusterSizeEvolutionSPH.at(iCounterAddress))->Write();
    (fAcrossPositionResidualEvolutionSPH.at(iCounterAddress))->Write();
    (fAlongPositionResidualEvolutionSPH.at(iCounterAddress))->Write();
    (fTimeResidualEvolutionSPH.at(iCounterAddress))->Write();
    (fCentralStripToTEvolutionSPH.at(iCounterAddress))->Write();
    (fDistanceToLastPoint.at(iCounterAddress))->Write();
    (fDistanceToLocalLastPoint.at(iCounterAddress))->Write();
    (fDistanceToLastHitUnbiased.at(iCounterAddress))->Write();
    (fDistanceToLastHitBiased.at(iCounterAddress))->Write();
    (fEfficiencyToLastHitUnbiased.at(iCounterAddress))->Write();
    (fEfficiencyToLastHitBiased.at(iCounterAddress))->Write();
    (fClusterSizeToLastHitUnbiased.at(iCounterAddress))->Write();
    (fClusterSizeToLastHitBiased.at(iCounterAddress))->Write();
    (fTimeResidualToLastHitUnbiased.at(iCounterAddress))->Write();
    (fTimeResidualToLastHitBiased.at(iCounterAddress))->Write();
    (fLocalSurfaceFlux.at(iCounterAddress))->Write();
    (fChannelDigiRates.at(iCounterAddress))->Write();
    (fNHitsPerPoint.at(iCounterAddress))->Write();
    (fNPointsPerHit.at(iCounterAddress))->Write();

    delete fClusterSizes.at(iCounterAddress);
    delete fStripToTSpectra.at(iCounterAddress);
    delete fStripPositionSpectra.at(iCounterAddress);
    delete fPointPositionSpectra.at(iCounterAddress);
    delete fHitPositionSpectra.at(iCounterAddress);
    delete fStripTimeResidualWalk.at(iCounterAddress);
    delete fHitAcrossPositionResidual.at(iCounterAddress);
    delete fHitAlongPositionResidual.at(iCounterAddress);
    delete fHitTimeResidual.at(iCounterAddress);
    delete fCentralClusterSizeEvolution.at(iCounterAddress);
    delete fCentralTimeResidualEvolution.at(iCounterAddress);
    delete fHitChiSquare.at(iCounterAddress);
    delete fCorrectedHitChiSquare.at(iCounterAddress);
    delete fClusterSizeEvolution.at(iCounterAddress);
    delete fAcrossPositionResidualEvolution.at(iCounterAddress);
    delete fAlongPositionResidualEvolution.at(iCounterAddress);
    delete fTimeResidualEvolution.at(iCounterAddress);
    delete fCentralStripToTEvolution.at(iCounterAddress);
    delete fIntegratedEfficiency.at(iCounterAddress);
    delete fChiSquareEfficiency.at(iCounterAddress);
    delete fCorrectedChiSquareEfficiency.at(iCounterAddress);
    delete fClusterSizesSPH.at(iCounterAddress);
    delete fStripToTSpectraSPH.at(iCounterAddress);
    delete fStripPositionSpectraSPH.at(iCounterAddress);
    delete fHitPositionSpectraSPH.at(iCounterAddress);
    delete fStripTimeResidualWalkSPH.at(iCounterAddress);
    delete fHitAcrossPositionResidualSPH.at(iCounterAddress);
    delete fHitAlongPositionResidualSPH.at(iCounterAddress);
    delete fHitTimeResidualSPH.at(iCounterAddress);
    delete fHitChiSquareSPH.at(iCounterAddress);
    delete fCorrectedHitChiSquareSPH.at(iCounterAddress);
    delete fClusterSizeEvolutionSPH.at(iCounterAddress);
    delete fAcrossPositionResidualEvolutionSPH.at(iCounterAddress);
    delete fAlongPositionResidualEvolutionSPH.at(iCounterAddress);
    delete fTimeResidualEvolutionSPH.at(iCounterAddress);
    delete fCentralStripToTEvolutionSPH.at(iCounterAddress);
    delete fDistanceToLastPoint.at(iCounterAddress);
    delete fDistanceToLocalLastPoint.at(iCounterAddress);
    delete fDistanceToLastHitUnbiased.at(iCounterAddress);
    delete fDistanceToLastHitBiased.at(iCounterAddress);
    delete fEfficiencyToLastHitUnbiased.at(iCounterAddress);
    delete fEfficiencyToLastHitBiased.at(iCounterAddress);
    delete fClusterSizeToLastHitUnbiased.at(iCounterAddress);
    delete fClusterSizeToLastHitBiased.at(iCounterAddress);
    delete fTimeResidualToLastHitUnbiased.at(iCounterAddress);
    delete fTimeResidualToLastHitBiased.at(iCounterAddress);
    delete fLocalSurfaceFlux.at(iCounterAddress);
    delete fChannelDigiRates.at(iCounterAddress);
    delete fNHitsPerPoint.at(iCounterAddress);
    delete fNPointsPerHit.at(iCounterAddress);

    // The following histograms belong to directory 'gROOT' and can be safely
    // deleted here.
    delete fAcrossPositionResidualEvolutionSigma.at(iCounterAddress);
    delete fAlongPositionResidualEvolutionSigma.at(iCounterAddress);
    delete fTimeResidualEvolutionSigma.at(iCounterAddress);
    delete fAcrossPositionResidualEvolutionSigmaSPH.at(iCounterAddress);
    delete fAlongPositionResidualEvolutionSigmaSPH.at(iCounterAddress);
    delete fTimeResidualEvolutionSigmaSPH.at(iCounterAddress);

    gDirectory->cd("../");
  }

  fHistogramFile->Close();
  delete fHistogramFile;

  tOldDirectory->cd();
}
// ---------------------------------------------------------------------------

ClassImp(CbmTofRecoMCQa)