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


#include "CbmTofClusterMC.h"

#include "CbmTofDigiTbPar.h"
#include "CbmTofCounter.h"
#include "CbmTofDigiExp.h"
#include "CbmTofHit.h"
#include "CbmTofWall.h"
#include "CbmTofPoint.h"
#include "CbmMCTrack.h"
#include "CbmMatch.h"
#include "CbmMCDataManager.h"
#include "CbmMCDataArray.h"
#include "CbmEvent.h"
#include "CbmTofDef.h"

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

#include "TClonesArray.h"
#include "TMath.h"
#include "TGeoManager.h"

#include <set>


// ---------------------------------------------------------------------------
CbmTofClusterMC::CbmTofClusterMC() :
  FairTask("CbmTofClusterMC", 0),
  fEvents(NULL),
  fTofPoints(NULL),
  fMCTracks(NULL),
  fTofPointsTB(NULL),
  fTofDigis(NULL),
  fTofCalDigis(NULL),
  fTofHits(NULL),
  fTofHitMatches(NULL),
  fMCDigiMap(),
  fbTimeSortHits(kTRUE),
  fDigiTbParSet(NULL),
  fbPointsInTS(kFALSE),
  fPointClusterMap(),
  fReferencePointMap(),
  fbBuildHitsFromPointClusters(kFALSE),
  fbHitCompatibilityMode(kFALSE),
  fbAlternativeBranchNames(kFALSE),
  fiFileIndex(0),
  fDuplicatedCounterSides(),
  fbBuildHitsFromRefTracks(kFALSE)
{
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
CbmTofClusterMC::~CbmTofClusterMC()
{
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
void CbmTofClusterMC::Exec(Option_t */*option*/)
{
  fTofHits->Clear();
  fTofHitMatches->Delete();

  fPointClusterMap.clear();
  fReferencePointMap.clear();

  // WARNING: Make sure that if a 'return' statement is introduced into
  //          this method in the future, the array pointers are reset correctly
  //          PRIOR to returning!
  TClonesArray* tDigiInputArray(NULL);

  if(fDuplicatedCounterSides.size())
  {
    fTofCalDigis->Delete();

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

      // As the entire downstream 'FairTask' chain is migrated to "TofCalDigi"
      // in this scenario, every single digi needs to be copied to the
      // corresponding output array.
      new((*fTofCalDigis)[fTofCalDigis->GetEntriesFast()]) CbmTofDigiExp(*tDigi);

      Int_t iModuleType   = tDigi->GetType();
      Int_t iModuleIndex  = tDigi->GetSm();
      Int_t iCounterIndex = tDigi->GetRpc();
      Int_t iChannelIndex = tDigi->GetChannel();
      Int_t iCounterSide  = tDigi->GetSide();

      // In addition, some duplicate digis are created on given counter sides.
      auto CounterSideTuple = std::make_tuple(iModuleType, iModuleIndex, iCounterIndex, iCounterSide);

      if(fDuplicatedCounterSides.find(CounterSideTuple) != fDuplicatedCounterSides.end())
      {
        CbmTofDigiExp* tDuplicateDigi = new((*fTofCalDigis)[fTofCalDigis->GetEntriesFast()]) CbmTofDigiExp(*tDigi);

        if(0 == iCounterSide)
        {
          tDuplicateDigi->SetAddress(iModuleIndex, iCounterIndex, iChannelIndex, 1, iModuleType);
        }
        else
        {
          tDuplicateDigi->SetAddress(iModuleIndex, iCounterIndex, iChannelIndex, 0, iModuleType);
        }
      }
    }

    // Digi time sorting after inserting duplicates is not strictly necessary here as
    // MC backlink information (and not spatio-temporal distances between digi pairs)
    // is used for ideal hit building.
    fTofCalDigis->Sort();

    tDigiInputArray = fTofDigis;
    fTofDigis = fTofCalDigis;
  }


  Int_t iTSIndex = FairRootManager::Instance()->GetEntryNr();

  // Create point clusters and the corresponding maps for hit building.
  if(fbBuildHitsFromPointClusters)
  {
    if(fbPointsInTS)
    {
      Double_t dLocalTrackStart[3] = {0., 0., 0.};
      Double_t dGlobalTrackStart[3] = {0., 0., 0.};
      Double_t dLocalMotherTrackStart[3] = {0., 0., 0.};
      Double_t dGlobalMotherTrackStart[3] = {0., 0., 0.};

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

        // Extract the file ID of the associated MC track from 'CbmTofPoint::fUniqueID'
        // before reusing this member variable to store the CbmTofPointType of
        // the MC point for further online analysis (the modified point object
        // is not written to disk).
        Int_t iFileID = tPoint->GetUniqueID();

        // Mark all ToF points in the TS as internal cluster points by default.
        tPoint->SetUniqueID(tof::point_Internal);

        Double_t dPointTime = tPoint->GetTime();

        Int_t iMCEvent = tPoint->GetEventID();
        Int_t iTrackID = tPoint->GetTrackID();
        Int_t iCounterID = tPoint->GetDetectorID();

        CbmTofCounter* tCounter = CbmTofWall::Instance()->GetCounter(iCounterID);
        auto & tCounterPointClusters = fPointClusterMap[iCounterID];
        CbmTofPointCluster* tPointCluster(NULL);

        // beam point or dark point (no secondaries available)
        if(-1 == iTrackID)
        {
          // Only consider beam points here as all dark points would be assigned
          // to a single 'CbmTofPointCluster' due to their (-1, -1) track ID pair.
          // To every counter at most one beam point can be assigned per event
          // which makes beam point clusters (always of size 1) distinguishable.
          // -> special treatment for dark points required later on!
          if(-1 != iMCEvent)
          {
            tPointCluster = &(((tCounterPointClusters.emplace(std::make_tuple(iFileID, iMCEvent, iTrackID), CbmTofPointCluster())).first)->second);
            tPointCluster->SetCharged(kTRUE);
            tPointCluster->AddPoint(fiFileIndex, iTSIndex, iPoint, dPointTime, kTRUE);
          }
          else
          {
            if(!fEvents)
            {
              // As no point clusters are created for dark points, their 'fUniqueID'
              // needs to be set here.
              tPoint->SetUniqueID(tof::point_Dark);
            }
          }
        }
        else
        {
          CbmMCTrack* tTrack = dynamic_cast<CbmMCTrack*>(fMCTracks->Get(iFileID, iMCEvent, iTrackID));

          dGlobalTrackStart[0] = tTrack->GetStartX();
          dGlobalTrackStart[1] = tTrack->GetStartY();
          dGlobalTrackStart[2] = tTrack->GetStartZ();

          tCounter->MasterToLocalBox(dGlobalTrackStart, dLocalTrackStart);

          Int_t iMotherID = tPoint->GetTrackID();
          Int_t iKinshipDegree = 0;
          CbmMCTrack* tMotherTrack = NULL;

          while(-1 != iMotherID)
          {
            tMotherTrack = dynamic_cast<CbmMCTrack*>(fMCTracks->Get(iFileID, iMCEvent, iMotherID));

            dGlobalMotherTrackStart[0] = tMotherTrack->GetStartX();
            dGlobalMotherTrackStart[1] = tMotherTrack->GetStartY();
            dGlobalMotherTrackStart[2] = tMotherTrack->GetStartZ();

            tCounter->MasterToLocalBox(dGlobalMotherTrackStart, dLocalMotherTrackStart);

            // lowest-ranked mother particle created outside the counter box volume
            if(!tCounter->IsContainedInBox(dLocalMotherTrackStart))
            {
              std::tuple<Int_t, Int_t, Int_t> tMotherTag(iFileID, iMCEvent, iMotherID);

              auto itTemp = tCounterPointClusters.find(tMotherTag);

              if(itTemp == tCounterPointClusters.end())
              {
                tPointCluster = &tCounterPointClusters[tMotherTag];

                if(!tMotherTrack->GetCharge())
                {
                  tPointCluster->SetCharged(kFALSE);
                }

                tPointCluster->AddPoint(fiFileIndex, iTSIndex, iPoint, dPointTime, 0 == iKinshipDegree);
              }
              else
              {
                tPointCluster = &tCounterPointClusters.at(tMotherTag);

                if(tCounter->IsContainedInBox(dLocalTrackStart))
                {
                  tPointCluster->AddPoint(fiFileIndex, iTSIndex, iPoint, dPointTime, 0 == iKinshipDegree);
                }
              }

              break;
            }

            iMotherID = tMotherTrack->GetMotherId();
            iKinshipDegree++;
          }
        }

      }

    }
    else
    {
      std::set<std::pair<Int_t, Int_t>> iMCEventSet;

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

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

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

          Int_t iFileID = tLink.GetFile();
          Int_t iMCEvent = tLink.GetEntry();

          if(-1 < iMCEvent)
          {
            iMCEventSet.emplace(iFileID, iMCEvent);
          }

        }
      }

      for(auto const & iFileEventPair : iMCEventSet)
      {
        Int_t iFileID = iFileEventPair.first;
        Int_t iMCEvent = iFileEventPair.second;

        Double_t dLocalTrackStart[3] = {0., 0., 0.};
        Double_t dGlobalTrackStart[3] = {0., 0., 0.};
        Double_t dLocalMotherTrackStart[3] = {0., 0., 0.};
        Double_t dGlobalMotherTrackStart[3] = {0., 0., 0.};

        for(Int_t iPoint = 0; iPoint < fTofPoints->Size(iFileID, iMCEvent); iPoint++)
        {
          CbmTofPoint* tPoint = dynamic_cast<CbmTofPoint*>(fTofPoints->Get(iFileID, iMCEvent, iPoint));

          Double_t dPointTime = tPoint->GetTime();

          Int_t iTrackID = tPoint->GetTrackID();
          Int_t iCounterID = tPoint->GetDetectorID();

          CbmTofCounter* tCounter = CbmTofWall::Instance()->GetCounter(iCounterID);
          auto & tCounterPointClusters = fPointClusterMap[iCounterID];
          CbmTofPointCluster* tPointCluster(NULL);

          // beam point (no secondaries available)
          if(-1 == iTrackID)
          {
            // To every counter at most one beam point can be assigned per event
            // which makes beam point clusters (always of size 1) distinguishable.
            tPointCluster = &(((tCounterPointClusters.emplace(std::make_tuple(iFileID, iMCEvent, iTrackID), CbmTofPointCluster())).first)->second);
            tPointCluster->SetCharged(kTRUE);
            tPointCluster->AddPoint(iFileID, iMCEvent, iPoint, dPointTime, kTRUE);
          }
          else
          {
            CbmMCTrack* tTrack = dynamic_cast<CbmMCTrack*>(fMCTracks->Get(iFileID, iMCEvent, iTrackID));

            dGlobalTrackStart[0] = tTrack->GetStartX();
            dGlobalTrackStart[1] = tTrack->GetStartY();
            dGlobalTrackStart[2] = tTrack->GetStartZ();

            tCounter->MasterToLocalBox(dGlobalTrackStart, dLocalTrackStart);

            Int_t iMotherID = tPoint->GetTrackID();
            Int_t iKinshipDegree = 0;
            CbmMCTrack* tMotherTrack = NULL;

            while(-1 != iMotherID)
            {
              tMotherTrack = dynamic_cast<CbmMCTrack*>(fMCTracks->Get(iFileID, iMCEvent, iMotherID));

              dGlobalMotherTrackStart[0] = tMotherTrack->GetStartX();
              dGlobalMotherTrackStart[1] = tMotherTrack->GetStartY();
              dGlobalMotherTrackStart[2] = tMotherTrack->GetStartZ();

              tCounter->MasterToLocalBox(dGlobalMotherTrackStart, dLocalMotherTrackStart);

              // lowest-ranked mother particle created outside the counter box volume
              if(!tCounter->IsContainedInBox(dLocalMotherTrackStart))
              {
                std::tuple<Int_t, Int_t, Int_t> tMotherTag(iFileID, iMCEvent, iMotherID);

                auto itTemp = tCounterPointClusters.find(tMotherTag);

                if(itTemp == tCounterPointClusters.end())
                {
                  tPointCluster = &tCounterPointClusters[tMotherTag];

                  if(!tMotherTrack->GetCharge())
                  {
                    tPointCluster->SetCharged(kFALSE);
                  }

                  tPointCluster->AddPoint(iFileID, iMCEvent, iPoint, dPointTime, 0 == iKinshipDegree);
                }
                else
                {
                  tPointCluster = &tCounterPointClusters.at(tMotherTag);

                  if(tCounter->IsContainedInBox(dLocalTrackStart))
                  {
                    tPointCluster->AddPoint(iFileID, iMCEvent, iPoint, dPointTime, 0 == iKinshipDegree);
                  }
                }

                break;
              }

              iMotherID = tMotherTrack->GetMotherId();
              iKinshipDegree++;
            }
          }
        }
      }
    }


    for(auto const & itCounter : fPointClusterMap)
    {
      for(auto const & itCluster: itCounter.second)
      {
        const CbmTofPointCluster* tPointCluster = &itCluster.second;

        // TODO: add an option to select direct points as hit building reference only
        fReferencePointMap.emplace(tPointCluster->GetReferencePoint(), tPointCluster->IsCharged() ? tof::point_ExternalCharged : tof::point_ExternalNeutral);
      }
    }


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

        auto itTemp = fReferencePointMap.find(std::make_tuple(fiFileIndex, iTSIndex, iPoint));

        if(itTemp != fReferencePointMap.end())
        {
          tPoint->SetUniqueID(itTemp->second);
        }

        // Consider all dark points for hit building.
        if(tof::point_Dark == tPoint->GetUniqueID())
        {
          fReferencePointMap.emplace(std::make_tuple(fiFileIndex, iTSIndex, iPoint), tof::point_Dark);
        }
      }
    }

  }
  // Find reference track points and mark them for hit building and QA.
  else if(fbBuildHitsFromRefTracks)
  {
    TGeoNode* tNode(NULL);
    TGeoMedium* tMedium(NULL);
    TGeoMaterial* tMaterial(NULL);

    const char* cMaterialName;

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

        // Extract the file ID of the associated MC track from 'CbmTofPoint::fUniqueID'
        // before reusing this member variable to store the CbmTofPointType of
        // the MC point for further online analysis (the modified point object
        // is not written to disk).
        Int_t iFileID = tPoint->GetUniqueID();
        Int_t iMCEvent = tPoint->GetEventID();
        Int_t iTrackID = tPoint->GetTrackID();

        // Mark all ToF points in the TS as internal cluster points by default.
        tPoint->SetUniqueID(tof::point_Internal);

        // MC track object cannot be retrieved.
        if(-1 == iTrackID)
        {
          // Mark beam points for hit building.
          if(-1 != iMCEvent)
          {
            fReferencePointMap.emplace(std::make_tuple(fiFileIndex, iTSIndex, iPoint), tof::point_RefTrack);
            tPoint->SetUniqueID(tof::point_RefTrack);
          }
        }
        // MC track object can be retrieved.
        else
        {
          CbmMCTrack* tTrack = dynamic_cast<CbmMCTrack*>(fMCTracks->Get(iFileID, iMCEvent, iTrackID));

          tNode = gGeoManager->FindNode(tTrack->GetStartX(), tTrack->GetStartY(), tTrack->GetStartZ());
          tMedium = tNode->GetMedium();
          tMaterial = tMedium->GetMaterial();

          tof::GetMaterialName(tMaterial->GetName(), cMaterialName);

          if(0 == std::strcmp("target", cMaterialName) && 3 <= tTrack->GetNPoints(kTof))
          {
            fReferencePointMap.emplace(std::make_tuple(fiFileIndex, iTSIndex, iPoint), tof::point_RefTrack);
            tPoint->SetUniqueID(tof::point_RefTrack);
          }
        }
      }

    }
    else
    {
      std::set<std::pair<Int_t, Int_t>> iMCEventSet;

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

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

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

          Int_t iFileID = tLink.GetFile();
          Int_t iMCEvent = tLink.GetEntry();

          if(-1 < iMCEvent)
          {
            iMCEventSet.emplace(iFileID, iMCEvent);
          }

        }
      }

      for(auto const & iFileEventPair : iMCEventSet)
      {
        Int_t iFileID = iFileEventPair.first;
        Int_t iMCEvent = iFileEventPair.second;

        for(Int_t iPoint = 0; iPoint < fTofPoints->Size(iFileID, iMCEvent); iPoint++)
        {
          CbmTofPoint* tPoint = dynamic_cast<CbmTofPoint*>(fTofPoints->Get(iFileID, iMCEvent, iPoint));

          Int_t iTrackID = tPoint->GetTrackID();

          // MC track object cannot be retrieved (here: beam point).
          if(-1 == iTrackID)
          {
            fReferencePointMap.emplace(std::make_tuple(iFileID, iMCEvent, iPoint), tof::point_RefTrack);
          }
          // MC track object can be retrieved.
          else
          {
            CbmMCTrack* tTrack = dynamic_cast<CbmMCTrack*>(fMCTracks->Get(iFileID, iMCEvent, iTrackID));

            tNode = gGeoManager->FindNode(tTrack->GetStartX(), tTrack->GetStartY(), tTrack->GetStartZ());
            tMedium = tNode->GetMedium();
            tMaterial = tMedium->GetMaterial();

            tof::GetMaterialName(tMaterial->GetName(), cMaterialName);

            if(0 == std::strcmp("target", cMaterialName) && 3 <= tTrack->GetNPoints(kTof))
            {
              fReferencePointMap.emplace(std::make_tuple(iFileID, iMCEvent, iPoint), tof::point_RefTrack);
            }
          }
        }
      }
    }

  }


  // Build hits.
  if(fEvents)
  {
    for(Int_t iEvent = 0; iEvent < fEvents->GetEntriesFast(); iEvent++)
    {
      CbmEvent* tEvent = dynamic_cast<CbmEvent*>(fEvents->At(iEvent));

      for (Int_t iDigi = 0; iDigi < tEvent->GetNofData(kTofDigi); iDigi++)
      {
        Int_t iDigiIndex = static_cast<Int_t>(tEvent->GetIndex(kTofDigi, iDigi));

        CbmTofDigiExp* tDigi = dynamic_cast<CbmTofDigiExp*>(fTofDigis->At(iDigiIndex));

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

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

          // Only consider primary MC signals.
          if(tof::signal_SourcePrimary == tLink.GetUniqueID())
          {
            std::tuple<Int_t, Int_t, Int_t> MCPointID;

            if(fbPointsInTS)
            {
              // MC reference point is not available in the current TS. Thus,
              // digis referring to it will not be clustered.
              if(-1 == tLink.GetIndex())
              {
                continue;
              }
              else
              {
                MCPointID = std::make_tuple(tLink.GetFile(), tLink.GetEntry(), tLink.GetIndex());
              }              
            }
            else
            {
              MCPointID = std::make_tuple(tLink.GetFile(), tLink.GetEntry(), tLink.GetIndex());
            }

            if(!(fbBuildHitsFromPointClusters || fbBuildHitsFromRefTracks) || fReferencePointMap.find(MCPointID) != fReferencePointMap.end())
            {
              auto tMCIterator = fMCDigiMap.find(MCPointID);

              if(tMCIterator != fMCDigiMap.end())
              {
                (tMCIterator->second).push_back(std::make_pair(tDigi, iDigiIndex));
              }
              else
              {
                tMCIterator = (fMCDigiMap.emplace(MCPointID, std::vector<std::pair<CbmTofDigiExp*, Int_t>>())).first;
                (tMCIterator->second).reserve(10);
                (tMCIterator->second).push_back(std::make_pair(tDigi, iDigiIndex));
              }
            }
          }
        }
      }

      for(auto & itMCPoint : fMCDigiMap)
      {
        CbmTofPoint* tTofPoint(NULL);

        if(fbPointsInTS)
        {
          tTofPoint = dynamic_cast<CbmTofPoint*>(fTofPointsTB->At(std::get<2>(itMCPoint.first)));
        }
        else
        {
          tTofPoint = dynamic_cast<CbmTofPoint*>(fTofPoints->Get(std::get<0>(itMCPoint.first), std::get<1>(itMCPoint.first), std::get<2>(itMCPoint.first)));
        }

        CbmTofCounter* tCounter = CbmTofWall::Instance()->GetCounter(tTofPoint->GetDetectorID());
        if(tCounter)
        {
          CbmTofHit* tNewHit = dynamic_cast<CbmTofHit*>(fTofHits->ConstructedAt(fTofHits->GetEntriesFast()));
          tNewHit->SetUniqueID(iEvent); // keep the association of the hit with the CbmEvent of its digis

          tCounter->BuildClusterMC(itMCPoint.second, tNewHit, fbHitCompatibilityMode);

          CbmMatch* tHitMatch(NULL);
          std::set<std::tuple<Int_t, Int_t, Int_t>> MCToHitContributions;

          // Add digi links to the hit match object
          for(auto const & itDigi: itMCPoint.second)
          {
            CbmTofDigiExp* tDigi = itDigi.first;

            if(tDigi)
            {
              if(!tHitMatch)
              {
                tHitMatch = new CbmMatch();
              }

              tHitMatch->AddLink(tDigi->GetTot(), itDigi.second, iTSIndex, fiFileIndex);

              for(Int_t iLink = 0; iLink < tDigi->GetMatch()->GetNofLinks(); iLink++)
              {
                const CbmLink& tLink = tDigi->GetMatch()->GetLink(iLink);

                // FIXME: Multiple dark rate signals that might have contributed to a cluster of digis
                // can only be told apart by their respective 'CbmLink::fEntry' and 'CbmLink::fIndex'
                // if 'CbmMCPointBuffer' continuously indexes dark rate points which are supposed to
                // be stored in a timeslice. Otherwise, each of them would try to insert the same
                // (-1, -1, 0) tuple into this map such that the real dark rate contribution to the hit
                // would possibly be underestimated.
                MCToHitContributions.emplace(tLink.GetFile(), tLink.GetEntry(), tLink.GetIndex());
              }

            }
          }

          // Check if cluster size is at least one (= valid hit)
          if(tHitMatch)
          {
            tNewHit->SetMatch(tHitMatch);
            if(!fbHitCompatibilityMode)
            {
              tNewHit->SetFlag(MCToHitContributions.size());
            }
          }
          else
          {
            fTofHits->RemoveAt(fTofHits->GetEntriesFast() - 1);
          }
          
        }
      }

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

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

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

        // Consider primary source and dark rate signals for clustering when
        // operating on timeslices.
        if(tof::kiPrimarySignalRemainder == tLink.GetUniqueID()%tof::kiPrimarySignalDivisor)
        {
          std::tuple<Int_t, Int_t, Int_t> MCPointID;

          if(fbPointsInTS)
          {
            // MC reference point is not available in the current TS. Thus,
            // digis referring to it will not be clustered.
            if(-1 == tLink.GetIndex())
            {
              continue;
            }
            else
            {
              MCPointID = std::make_tuple(tLink.GetFile(), tLink.GetEntry(), tLink.GetIndex());
            }
          }
          else
          {
            // Dark rate MC information is not available from the MC input chain,
            // i.e. if these points are not written to output timeslices during
            // digitization, digis referring to them will not be clustered.
            if(tof::signal_SourcePrimary == tLink.GetUniqueID())
            {
              MCPointID = std::make_tuple(tLink.GetFile(), tLink.GetEntry(), tLink.GetIndex());
            }
            else
            {
              continue;
            }
          }

          if(!(fbBuildHitsFromPointClusters || fbBuildHitsFromRefTracks) || fReferencePointMap.find(MCPointID) != fReferencePointMap.end())
          {
            auto tMCIterator = fMCDigiMap.find(MCPointID);

            if(tMCIterator != fMCDigiMap.end())
            {
              (tMCIterator->second).push_back(std::make_pair(tDigi, iDigi));
            }
            else
            {
              tMCIterator = (fMCDigiMap.emplace(MCPointID, std::vector<std::pair<CbmTofDigiExp*, Int_t>>())).first;
              (tMCIterator->second).reserve(10);
              (tMCIterator->second).push_back(std::make_pair(tDigi, iDigi));
            }
          }
        }
      }
    }

    for(auto & itMCPoint : fMCDigiMap)
    {
      CbmTofPoint* tTofPoint(NULL);

      if(fbPointsInTS)
      {
        tTofPoint = dynamic_cast<CbmTofPoint*>(fTofPointsTB->At(std::get<2>(itMCPoint.first)));
      }
      else
      {
        tTofPoint = dynamic_cast<CbmTofPoint*>(fTofPoints->Get(std::get<0>(itMCPoint.first), std::get<1>(itMCPoint.first), std::get<2>(itMCPoint.first)));
      }

      CbmTofCounter* tCounter = CbmTofWall::Instance()->GetCounter(tTofPoint->GetDetectorID());
      if(tCounter)
      {
        CbmTofHit* tNewHit = dynamic_cast<CbmTofHit*>(fTofHits->ConstructedAt(fTofHits->GetEntriesFast()));

        tCounter->BuildClusterMC(itMCPoint.second, tNewHit, fbHitCompatibilityMode);

        CbmMatch* tHitMatch(NULL);
        std::set<std::tuple<Int_t, Int_t, Int_t>> MCToHitContributions;

        // Add digi links to the hit match object
        for(auto const & itDigi: itMCPoint.second)
        {
          CbmTofDigiExp* tDigi = itDigi.first;

          if(tDigi)
          {
            if(!tHitMatch)
            {
              tHitMatch = new CbmMatch();
            }

            tHitMatch->AddLink(tDigi->GetTot(), itDigi.second, iTSIndex, fiFileIndex);

            for(Int_t iLink = 0; iLink < tDigi->GetMatch()->GetNofLinks(); iLink++)
            {
              const CbmLink& tLink = tDigi->GetMatch()->GetLink(iLink);

              // FIXME: Multiple dark rate signals that might have contributed to a cluster of digis
              // can only be told apart by their respective 'CbmLink::fEntry' and 'CbmLink::fIndex'
              // if 'CbmMCPointBuffer' continuously indexes dark rate points which are supposed to
              // be stored in a timeslice. Otherwise, each of them would try to insert the same
              // (-1, -1, 0) tuple into this map such that the real dark rate contribution to the hit
              // would possibly be underestimated.
              MCToHitContributions.emplace(tLink.GetFile(), tLink.GetEntry(), tLink.GetIndex());
            }

          }
        }

        // Check if cluster size is at least one (= valid hit)
        if(tHitMatch)
        {
          tNewHit->SetMatch(tHitMatch);
          if(!fbHitCompatibilityMode)
          {
            tNewHit->SetFlag(MCToHitContributions.size());
          }
        }
        else
        {
          fTofHits->RemoveAt(fTofHits->GetEntriesFast() - 1);
        }
        
      }
    }

    fMCDigiMap.clear();
  }


  if(fbTimeSortHits)
  {
    fTofHits->Sort();
  }

  if(fEvents)
  {
    for(Int_t iHit = 0; iHit < fTofHits->GetEntriesFast(); iHit++)
    {
      CbmTofHit* tHit = dynamic_cast<CbmTofHit*>(fTofHits->At(iHit));

      // The hit object and the hit match object are related by equal indices in
      // their respective output array.
			CbmMatch* tMatch = new((*fTofHitMatches)[iHit]) CbmMatch();
      tMatch->AddLinks(*tHit->GetMatch());

      // Although the 'CbmMatch' object pointed to by 'CbmHit::fMatch' is not
      // supposed to be serialized along with the hit but to be stored in a
      // dedicated matching branch, ROOT anyway creates and serializes a default
      // 'CbmMatch' object when the hit is written to file that is created here
      // a priori to overwrite the actual matching object which is deleted by
      // the method call below.
      tHit->SetMatch(new CbmMatch());

      CbmEvent* tEvent = dynamic_cast<CbmEvent*>(fEvents->At(tHit->GetUniqueID()));
      tEvent->AddData(kTofHit, iHit);

      tHit->SetUniqueID(0);
    }
  }
  else
  {
    for(Int_t iHit = 0; iHit < fTofHits->GetEntriesFast(); iHit++)
    {
      CbmTofHit* tHit = dynamic_cast<CbmTofHit*>(fTofHits->At(iHit));

      // The hit object and the hit match object are related by equal indices in
      // their respective output array.
			CbmMatch* tMatch = new((*fTofHitMatches)[iHit]) CbmMatch();
      tMatch->AddLinks(*tHit->GetMatch());

      // Although the 'CbmMatch' object pointed to by 'CbmHit::fMatch' is not
      // supposed to be serialized along with the hit but to be stored in a
      // dedicated matching branch, ROOT anyway creates and serializes a default
      // 'CbmMatch' object when the hit is written to file that is created here
      // a priori to overwrite the actual matching object which is deleted by
      // the method call below.
      tHit->SetMatch(new CbmMatch());
    }
  }


  if(fDuplicatedCounterSides.size())
  {
    fTofDigis = tDigiInputArray;
  }
}
// ---------------------------------------------------------------------------



// ---------------------------------------------------------------------------
InitStatus CbmTofClusterMC::Init()
{
  if(fbBuildHitsFromPointClusters && fbBuildHitsFromRefTracks)
  {
    LOG(FATAL)<<"Cannot request hit building from point clusters and reference tracks simultaneously."<<FairLogger::endl;
  }

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

  // for future use
  for(auto const & itCounter : CbmTofWall::Instance()->GetCounterMap())
  {
    (itCounter.second)->SetFileIndex(fiFileIndex);
  }

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

  // --- Get input array (CbmEvent)
  // If an array named "Event" cannot be retrieved from 'FairRootManager',
  // the class runs in timeslice mode.
  fEvents = dynamic_cast<TClonesArray*>(tRootManager->GetObject("Event"));

  TString tPointBranchName;
  TString tDigiBranchName;
  TString tHitBranchName;
  TString tHitDigiMatchBranchName;

  // If 'CbmEvent' objects are available, there is no need to use alternative
  // branch names.
  if(fbAlternativeBranchNames && !fEvents)
  {
    tPointBranchName        = "ATofPointTB";
    tDigiBranchName         = "ATofDigiExp";
    tHitBranchName          = "ATofHit";
    tHitDigiMatchBranchName = "ATofDigiMatch";
  }
  else
  {
    tPointBranchName        = "TofPointTB";
    tDigiBranchName         = "TofDigiExp";
    tHitBranchName          = "TofHit";
    tHitDigiMatchBranchName = "TofDigiMatch";
  }

  // --- Get input array (CbmTofPoint)
  fTofPointsTB = dynamic_cast<TClonesArray*>(tRootManager->GetObject(tPointBranchName));

  CbmMCDataManager* tMCManager = dynamic_cast<CbmMCDataManager*>(tRootManager->GetObject("MCDataManager"));
  if(!tMCManager)
  {
    LOG(ERROR)<<"Could not retrieve CbmMCDataManager from FairRootManager."<<FairLogger::endl;
    return kFATAL;
  }

  if(fTofPointsTB)
  {
    fbPointsInTS = kTRUE;
  }
  else
  {
    fTofPoints = tMCManager->InitBranch("TofPoint");
    if(!fTofPoints)
    {
      LOG(ERROR)<<"Could not retrieve branch TofPoint from CbmMCDataManager."<<FairLogger::endl;
      return kFATAL;
    }
  }

  fMCTracks = tMCManager->InitBranch("MCTrack");
  if(!fMCTracks)
  {
    LOG(ERROR)<<"Could not retrieve branch MCTrack from CbmMCDataManager."<<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;
    }
  }
  else
  {
    // If a digi branch named "TofCalDigi" is found, an upstream instance of 'CbmTofTestBeamClusterizer'
    // most likely already duplicated digis from counters with one-sided readout. In this case, disable
    // the feature for this class.
    fDuplicatedCounterSides.clear();
  }

  if(fDuplicatedCounterSides.size())
  {
    fTofCalDigis = new TClonesArray("CbmTofDigiExp", 1000);
    tRootManager->Register("TofCalDigi", "TOF cal digis", fTofCalDigis, IsOutputBranchPersistent("TofCalDigi"));
  }

  // --- Register output array (CbmTofHit)
  fTofHits = new TClonesArray("CbmTofHit", 10000);
  tRootManager->Register(tHitBranchName, "CbmTofHit", fTofHits, IsOutputBranchPersistent(tHitBranchName));

  // --- Register output array (CbmMatch)
  fTofHitMatches = new TClonesArray("CbmMatch", 10000);
  tRootManager->Register(tHitDigiMatchBranchName, "CbmMatch", fTofHitMatches, IsOutputBranchPersistent(tHitDigiMatchBranchName));

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



// ---------------------------------------------------------------------------
void CbmTofClusterMC::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 CbmTofClusterMC::SetDuplicateCounterSide(Int_t iModuleType, Int_t iModuleIndex, Int_t iCounterIndex, Int_t iCounterSide)
{
  fDuplicatedCounterSides.emplace(std::make_tuple(iModuleType, iModuleIndex, iCounterIndex, iCounterSide));
}
// ---------------------------------------------------------------------------


ClassImp(CbmTofClusterMC)