CbmDeviceHitBuilderTof.cxx

Go to the documentation of this file.

/* Copyright (C) 2018-2021 PI-UHd, GSI
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Norbert Herrmann [committer] */
 
#include "CbmDeviceHitBuilderTof.h"
 
// CBM Classes and includes
#include "CbmDigiManager.h"
 
// TOF Classes and includes
#include "CbmMatch.h"
#include "CbmTofAddress.h"  // in cbmdata/tof
#include "CbmTofCell.h"     // in tof/TofData
#include "CbmTofClusterizersDef.h"
#include "CbmTofCreateDigiPar.h"    // in tof/TofTools
#include "CbmTofDetectorId_v12b.h"  // in cbmdata/tof
#include "CbmTofDetectorId_v14a.h"  // in cbmdata/tof
#include "CbmTofDetectorId_v21a.h"  // in cbmdata/tof
#include "CbmTofDigi.h"             // in cbmdata/tof
#include "CbmTofDigiBdfPar.h"       // in tof/TofParam
#include "CbmTofDigiPar.h"          // in tof/TofParam
#include "CbmTofGeoHandler.h"       // in tof/TofTools
#include "CbmTofHit.h"              // in cbmdata/tof
#include "CbmTofPoint.h"            // in cbmdata/tof
 
#include "FairEventHeader.h"
#include "FairFileHeader.h"
#include "FairGeoParSet.h"
#include "FairMQLogger.h"
#include "FairMQProgOptions.h"  // device->fConfig
#include "FairRootFileSink.h"
#include "FairRootManager.h"
#include "FairRunOnline.h"
#include "FairRuntimeDb.h"
#include "FairSource.h"
 
// ROOT Classes and includes
#include "TClonesArray.h"
#include "TDirectory.h"
#include "TF1.h"
#include "TF2.h"
#include "TGeoManager.h"
#include "TH1.h"
#include "TH2.h"
#include "TLine.h"
#include "TMath.h"
#include "TMinuit.h"
#include "TProfile.h"
#include "TROOT.h"
#include "TRandom3.h"
#include "TVector3.h"
 
#include <boost/archive/binary_iarchive.hpp>
#include <boost/archive/binary_oarchive.hpp>
#include <boost/serialization/vector.hpp>
 
#include <chrono>
#include <iomanip>
#include <stdexcept>
#include <string>
#include <thread>  // this_thread::sleep_for
struct InitTaskError : std::runtime_error {
  using std::runtime_error::runtime_error;
};
 
using namespace std;
 
// Constants definitions
static Int_t iMess                = 0;
static Int_t iIndexDut            = 0;
static Double_t StartAnalysisTime = 0.;
//const Double_t cLight             = 29.9792;  // in cm/ns
static FairRootManager* rootMgr = NULL;
static Int_t iRunId             = 1;
static Int_t SelMask            = DetMask;
static Double_t dTstart         = 0.;
static Double_t dTmax           = 0.;
 
CbmTofDigi* pRef;
CbmTofDigi* pRefCal;
CbmDigiManager* fDigiMan;
 
CbmDeviceHitBuilderTof::CbmDeviceHitBuilderTof()
  : fNumMessages(0)
  , fGeoMan(NULL)
  , fGeoHandler(new CbmTofGeoHandler())
  , fTofId(NULL)
  , fDigiPar(NULL)
  , fChannelInfo(NULL)
  , fDigiBdfPar(NULL)
  , fiNDigiIn(0)
  , fvDigiIn()
  , fEventHeader()
  , fEvtHeader(NULL)
  , fTofHitsColl(NULL)
  , fTofDigiMatchColl(NULL)
  , fTofHitsCollOut(NULL)
  , fTofDigiMatchCollOut(NULL)
  , fiNbHits(0)
  , fiNevtBuild(0)
  , fiMsgCnt(100)
  , fdTOTMax(50.)
  , fdTOTMin(0.)
  , fdTTotMean(2.)
  , fdMaxTimeDist(0.)
  , fdMaxSpaceDist(0.)
  , fdEvent(0)
  , fiMaxEvent(-1)
  , fiRunId(111)
  , fiOutputTreeEntry(0)
  , fiFileIndex(0)
  , fStorDigi()
  , fStorDigiInd()
  , vDigiIndRef()
  , fviClusterMul()
  , fviClusterSize()
  , fviTrkMul()
  , fvdX()
  , fvdY()
  , fvdDifX()
  , fvdDifY()
  , fvdDifCh()
  , fvCPDelTof()
  , fvCPTOff()
  , fvCPTotGain()
  , fvCPTotOff()
  , fvCPWalk()
  , fvLastHits()
  , fvDeadStrips()
  , fvPulserOffset()
  , fvPulserTimes()
  , fhEvDetMul(NULL)
  , fhEvDigiMul(NULL)
  , fhEvRateIn(NULL)
  , fhEvRateOut(NULL)
  , fhPulMul(NULL)
  , fhDigiTdif(NULL)
  , fhPulserTimesRaw()
  , fhPulserTimeRawEvo()
  , fhPulserTimesCor()
  , fhDigiTimesRaw()
  , fhDigiTimesCor()
  , fhRpcDigiTot()
  , fhRpcDigiCor()
  , fhRpcCluMul()
  , fhRpcCluRate()
  , fhRpcCluPosition()
  , fhRpcCluDelPos()
  , fhRpcCluDelMatPos()
  , fhRpcCluTOff()
  , fhRpcCluDelTOff()
  , fhRpcCluDelMatTOff()
  , fhRpcCluTrms()
  , fhRpcCluTot()
  , fhRpcCluSize()
  , fhRpcCluAvWalk()
  , fhRpcCluAvLnWalk()
  , fhRpcCluWalk()
  , fhSmCluPosition()
  , fhSmCluTOff()
  , fhSmCluSvel()
  , fhSmCluFpar()
  , fhRpcDTLastHits()
  , fhRpcDTLastHits_Tot()
  , fhRpcDTLastHits_CluSize()
  , fhTRpcCluMul()
  , fhTRpcCluPosition()
  , fhTRpcCluTOff()
  , fhTRpcCluTot()
  , fhTRpcCluSize()
  , fhTRpcCluAvWalk()
  , fhTRpcCluDelTof()
  , fhTRpcCludXdY()
  , fhTRpcCluWalk()
  , fhTSmCluPosition()
  , fhTSmCluTOff()
  , fhTSmCluTRun()
  , fhTRpcCluTOffDTLastHits()
  , fhTRpcCluTotDTLastHits()
  , fhTRpcCluSizeDTLastHits()
  , fhTRpcCluMemMulDTLastHits()
  , fhSeldT()
  , dTRef(0.)
  , fCalMode(0)
  , fdCaldXdYMax(0.)
  , fiCluMulMax(0)
  , fTRefMode(0)
  , fTRefHits(0)
  , fDutId(0)
  , fDutSm(0)
  , fDutRpc(0)
  , fDutAddr(0)
  , fSelId(0)
  , fSelSm(0)
  , fSelRpc(0)
  , fSelAddr(0)
  , fSel2Id(0)
  , fSel2Sm(0)
  , fSel2Rpc(0)
  , fSel2Addr(0)
  , fiMode(0)
  , fiPulserMode(0)
  , fiPulMulMin(0)
  , fiPulDetRef(0)
  , fiPulTotMin(0)
  , fiPulTotMax(1000)
  , fDetIdIndexMap()
  , fviDetId()
  , fPosYMaxScal(0.)
  , fTRefDifMax(0.)
  , fTotMax(0.)
  , fTotMin(0.)
  , fTotMean(0.)
  , fdDelTofMax(60.)
  , fMaxTimeDist(0.)
  , fdChannelDeadtime(0.)
  , fdMemoryTime(0.)
  , fEnableMatchPosScaling(kTRUE)
  , fbPs2Ns(kFALSE)
  , fCalParFileName("")
  , fOutHstFileName("")
  , fOutRootFileName("")
  , fCalParFile(NULL)
  , fOutRootFile(NULL)
{
}
 
CbmDeviceHitBuilderTof::~CbmDeviceHitBuilderTof()
{
  if (NULL != fOutRootFile) {
    LOG(info) << "Finally close root file " << fOutRootFile->GetName();
    fOutRootFile->cd();
    rootMgr->LastFill();
    rootMgr->Write();
    WriteHistograms();
    fOutRootFile->Write();
    fOutRootFile->Close();
  }
}
 
void CbmDeviceHitBuilderTof::InitTask()
try {
  // Get the information about created channels from the device
  // Check if the defined channels from the topology (by name)
  // are in the list of channels which are possible/allowed
  // for the device
  // The idea is to check at initilization if the devices are
  // properly connected. For the time beeing this is done with a
  // nameing convention. It is not avoided that someone sends other
  // data on this channel.
  int noChannel = fChannels.size();
  LOG(info) << "Number of defined input channels: " << noChannel;
  for (auto const& entry : fChannels) {
    LOG(info) << "Channel name: " << entry.first;
    if (!IsChannelNameAllowed(entry.first)) throw InitTaskError("Channel name does not match.");
    if (entry.first != "syscmd") OnData(entry.first, &CbmDeviceHitBuilderTof::HandleData);
    else
      OnData(entry.first, &CbmDeviceHitBuilderTof::HandleMessage);
  }
  InitWorkspace();
  InitContainers();
  LoadGeometry();
  InitRootOutput();
}
catch (InitTaskError& e) {
  LOG(error) << e.what();
  ChangeState(fair::mq::Transition::ErrorFound);
}
 
bool CbmDeviceHitBuilderTof::IsChannelNameAllowed(std::string channelName)
{
  for (auto const& entry : fAllowedChannels) {
    std::size_t pos1 = channelName.find(entry);
    if (pos1 != std::string::npos) {
      const vector<std::string>::const_iterator pos =
        std::find(fAllowedChannels.begin(), fAllowedChannels.end(), entry);
      const vector<std::string>::size_type idx = pos - fAllowedChannels.begin();
      LOG(info) << "Found " << entry << " in " << channelName;
      LOG(info) << "Channel name " << channelName << " found in list of allowed channel names at position " << idx;
      return true;
    }
  }
  LOG(info) << "Channel name " << channelName << " not found in list of allowed channel names.";
  LOG(error) << "Stop device.";
  return false;
}
 
Bool_t CbmDeviceHitBuilderTof::InitWorkspace()
{
  LOG(info) << "Init work space for CbmDeviceHitBuilderTof.";
  fOutRootFileName = fConfig->GetValue<string>("OutRootFile");
  fiMaxEvent       = fConfig->GetValue<int64_t>("MaxEvent");
  LOG(info) << "Max number of events to be processed: " << fiMaxEvent;
  fiRunId      = fConfig->GetValue<int64_t>("RunId");
  fiMode       = fConfig->GetValue<int64_t>("Mode");
  fiPulserMode = fConfig->GetValue<int64_t>("PulserMode");
  fiPulMulMin  = fConfig->GetValue<uint64_t>("PulMulMin");
  fiPulDetRef  = fConfig->GetValue<uint64_t>("PulDetRef");
  fiPulTotMin  = fConfig->GetValue<uint64_t>("PulTotMin");
  fiPulTotMax  = fConfig->GetValue<uint64_t>("PulTotMax");
  dTmax        = (Double_t) fConfig->GetValue<uint64_t>("ReqTint");
 
  //fTofCalDigisColl     = new TClonesArray("CbmTofDigi", 100);
  //fTofCalDigisCollOut  = new TClonesArray("CbmTofDigi", 100);
  fTofCalDigiVec = new std::vector<CbmTofDigi>();
 
  fTofHitsColl         = new TClonesArray("CbmTofHit", 100);
  fTofHitsCollOut      = new TClonesArray("CbmTofHit", 100);
  fTofDigiMatchColl    = new TClonesArray("CbmMatch", 100);
  fTofDigiMatchCollOut = new TClonesArray("CbmMatch", 100);
 
  /*
  fDigiMan = CbmDigiManager::Instance();
  fDigiMan->Init();
  if (!fDigiMan->IsPresent(ECbmModuleId::kTof)) {
    LOG(error) << "HitBuilder: No digi input!";
    //return kFALSE;
  }
  */
  if (fOutRootFileName != "") {  // prepare root output
 
    FairRunOnline* fRun = new FairRunOnline(0);
    rootMgr             = FairRootManager::Instance();
    //fOutRootFile = rootMgr->OpenOutFile(fOutRootFileName);
    if (rootMgr->InitSink()) {
      fRun->SetSink(new FairRootFileSink(fOutRootFileName));
      fOutRootFile = rootMgr->GetOutFile();
      if (NULL == fOutRootFile) LOG(fatal) << "could not open root file";
    }
    else
      LOG(fatal) << "could not init Sink";
  }
 
  // steering variables
  fDutId  = fConfig->GetValue<uint64_t>("DutType");
  fDutSm  = fConfig->GetValue<uint64_t>("DutSm");
  fDutRpc = fConfig->GetValue<uint64_t>("DutRpc");
 
  fSelId  = fConfig->GetValue<uint64_t>("SelType");
  fSelSm  = fConfig->GetValue<uint64_t>("SelSm");
  fSelRpc = fConfig->GetValue<uint64_t>("SelRpc");
 
  fSel2Id  = fConfig->GetValue<uint64_t>("Sel2Type");
  fSel2Sm  = fConfig->GetValue<uint64_t>("Sel2Sm");
  fSel2Rpc = fConfig->GetValue<uint64_t>("Sel2Rpc");
 
  fiBeamRefType = fConfig->GetValue<uint64_t>("BRefType");
  fiBeamRefSm   = fConfig->GetValue<uint64_t>("BRefSm");
  fiBeamRefDet  = fConfig->GetValue<uint64_t>("BRefDet");
 
  return kTRUE;
}
 
Bool_t CbmDeviceHitBuilderTof::InitRootOutput()
{
  if (NULL != fOutRootFile) {
    LOG(info) << "Init Root Output to " << fOutRootFile->GetName();
 
    /*
    fFileHeader->SetRunId(iRunId);
    rootMgr->WriteFileHeader(fFileHeader);
    */
    rootMgr->InitSink();
    fEvtHeader = new FairEventHeader();
    fEvtHeader->SetRunId(iRunId);
    rootMgr->Register("EventHeader.", "Event", fEvtHeader, kTRUE);
    auto source = rootMgr->GetSource();
    if (source) { source->FillEventHeader(fEvtHeader); }
 
    // rootMgr->Register("CbmTofDigi", "Tof raw Digi", fTofCalDigisColl, kTRUE);
    //    fOutRootFile->cd();
    // rootMgr->RegisterAny("TofCalDigi", fTofCalDigiVec, kTRUE);
    rootMgr->RegisterAny("TofDigi", fTofCalDigiVec, kTRUE);
 
    TTree* outTree = new TTree(FairRootManager::GetTreeName(), "/cbmout", 99);
    LOG(info) << "define Tree " << outTree->GetName();
    //rootMgr->TruncateBranchNames(outTree, "cbmout");
    //rootMgr->SetOutTree(outTree);
    rootMgr->GetSink()->SetOutTree(outTree);
    rootMgr->WriteFolder();
    LOG(info) << "Initialized outTree with rootMgr at " << rootMgr;
    /*
    TFile* oldFile     = gFile;
    TDirectory* oldDir = gDirectory;
    fOutRootFile = new TFile(fOutRootFileName,"recreate");
    fRootEvent   = new TTree("CbmEvent","Cbm Event");
    fRootEvent->Branch("CbmDigi",fTofCalDigisColl);
    LOG(info)<<"Open Root Output file " << fOutRootFileName;
    fRootEvent->Write();
    gFile      = oldFile;
    gDirectory = oldDir;
    */
  }
  return kTRUE;
}
 
 
Bool_t CbmDeviceHitBuilderTof::InitContainers()
{
  LOG(info) << "Init parameter containers for CbmDeviceHitBuilderTof.";
 
  FairRuntimeDb* fRtdb = FairRuntimeDb::instance();
  if (NULL == fRtdb) LOG(error) << "No FairRuntimeDb found";
 
  // NewSimpleMessage creates a copy of the data and takes care of its destruction (after the transfer takes place).
  // Should only be used for small data because of the cost of an additional copy
 
  // Int_t fiRunId=1535700811;  // from *geo*.par.root file
  Int_t NSet = 3;
  std::string parSet[NSet];
  parSet[0]           = "CbmTofDigiPar";
  parSet[1]           = "CbmTofDigiBdfPar";
  parSet[2]           = "FairGeoParSet";
  std::string Channel = "parameters";
 
  Bool_t isSimulation = kFALSE;
  Int_t iGeoVersion;
  FairParSet* cont;
 
  for (Int_t iSet = 1; iSet < NSet; iSet++) {
    std::string message = parSet[iSet] + "," + to_string(fiRunId);
    LOG(info) << "Requesting parameter container, sending message: " << message;
 
    FairMQMessagePtr req(NewSimpleMessage(message));
    //FairMQMessagePtr req(NewSimpleMessage( "CbmTofDigiBdfPar,111" )); //original format
    FairMQMessagePtr rep(NewMessage());
    CbmTofCreateDigiPar* tofDigiPar = NULL;
 
    if (Send(req, Channel) > 0) {
      if (Receive(rep, Channel) >= 0) {
        if (rep->GetSize() != 0) {
          CbmMqTMessage tmsg(rep->GetData(), rep->GetSize());
          switch (iSet) {
            case 0:
              fDigiPar = static_cast<CbmTofDigiPar*>(tmsg.ReadObject(tmsg.GetClass()));
              //fDigiPar->print();
              break;
            case 1:
              fDigiBdfPar = static_cast<CbmTofDigiBdfPar*>(tmsg.ReadObject(tmsg.GetClass()));
              //fDigiBdfPar->print();
              LOG(info) << "Calib data file: " << fDigiBdfPar->GetCalibFileName();
              fdMaxTimeDist  = fDigiBdfPar->GetMaxTimeDist();     // in ns
              fdMaxSpaceDist = fDigiBdfPar->GetMaxDistAlongCh();  // in cm
 
              if (fMaxTimeDist != fdMaxTimeDist) {
                fdMaxTimeDist  = fMaxTimeDist;  // modify default
                fdMaxSpaceDist = fdMaxTimeDist * fDigiBdfPar->GetSignalSpeed()
                                 * 0.5;  // cut consistently on positions (with default signal velocity)
              }
 
              break;
            case 2:  // Geometry container
              cont = static_cast<FairParSet*>(tmsg.ReadObject(tmsg.GetClass()));
              cont->init();
              cont->Print();
              //fRtdb->InitContainer(parSet[iSet]);
              if (NULL == fGeoMan) fGeoMan = (TGeoManager*) ((FairGeoParSet*) cont)->GetGeometry();  //crashes
              LOG(info) << "GeoMan: " << fGeoMan << " " << gGeoManager;
              iGeoVersion = fGeoHandler->Init(isSimulation);
              if (k21a > iGeoVersion) {
                LOG(error) << "Incompatible geometry !!!";
                ChangeState(fair::mq::Transition::Stop);
              }
              switch (iGeoVersion) {
                case k14a: fTofId = new CbmTofDetectorId_v14a(); break;
                case k21a: fTofId = new CbmTofDetectorId_v21a(); break;
              }
              if (NULL == fDigiPar) {
                if (NULL == fRtdb) {
                  LOG(info) << "Instantiate FairRunDb";
                  fRtdb = FairRuntimeDb::instance();
                  assert(fRtdb);
                }
                // create digitization parameters from geometry file
                tofDigiPar = new CbmTofCreateDigiPar("TOF Digi Producer", "TOF task");
                LOG(info) << "Create DigiPar ";
                tofDigiPar->Init();
                fDigiPar = (CbmTofDigiPar*) (fRtdb->getContainer("CbmTofDigiPar"));
                if (NULL == fDigiPar) {
                  LOG(error) << "CbmTofEventClusterizer::InitParameters => Could not obtain CbmTofDigiPar";
                  return kFALSE;
                }
              }
              gGeoManager->Export("HitBuilder.geo.root");
              break;
 
            case 3:  // Calib
              break;
            default: LOG(warn) << "Parameter Set " << iSet << " not implemented ";
          }
          LOG(info) << "Received parameter from server:";
        }
        else {
          LOG(warn) << "Received empty reply. Parameter not available";
        }
      }
    }
  }
  Bool_t initOK = ReInitContainers();
 
  CreateHistograms();
 
  if (!InitCalibParameter()) return kFALSE;  // ChangeState(PAUSE); // for debugging
 
  fDutAddr  = CbmTofAddress::GetUniqueAddress(fDutSm, fDutRpc, 0, 0, fDutId);
  fSelAddr  = CbmTofAddress::GetUniqueAddress(fSelSm, fSelRpc, 0, 0, fSelId);
  fSel2Addr = CbmTofAddress::GetUniqueAddress(fSel2Sm, fSel2Rpc, 0, 0, fSel2Id);
  if (fiBeamRefType > -1) {
    if (fiBeamRefDet > -1) {
      fiBeamRefAddr = CbmTofAddress::GetUniqueAddress(fiBeamRefSm, fiBeamRefDet, 0, 0, fiBeamRefType);
    }
    else {
      SelMask       = ModMask;
      fiBeamRefAddr = CbmTofAddress::GetUniqueAddress(fiBeamRefSm, 0, 0, 0, fiBeamRefType);
    }
  }
  else
    fiBeamRefAddr = -1;
 
  iIndexDut = fDigiBdfPar->GetDetInd(fDutAddr);
  LOG(info) << Form("Use Dut 0x%08x, Sel 0x%08x, Sel2 0x%08x, BRef 0x%08x", fDutAddr, fSelAddr, fSel2Addr,
                    fiBeamRefAddr);
  return initOK;
}
 
Bool_t CbmDeviceHitBuilderTof::ReInitContainers()
{
  LOG(info) << "ReInit parameter containers for CbmDeviceHitBuilderTof.";
 
  return kTRUE;
}
 
// handler is called whenever a message arrives on "data", with a reference to the message and a sub-channel index (here 0)
//bool CbmDeviceHitBuilderTof::HandleData(FairMQMessagePtr& msg, int /*index*/)
bool CbmDeviceHitBuilderTof::HandleData(FairMQParts& parts, int /*index*/)
{
  // Don't do anything with the data
  // Maybe add an message counter which counts the incomming messages and add
  // an output
  fNumMessages++;
  LOG(debug) << "Received message " << fNumMessages << " with " << parts.Size() << " parts"
             << ", size0: " << parts.At(0)->GetSize();
 
  std::string msgStrE(static_cast<char*>(parts.At(0)->GetData()), (parts.At(0))->GetSize());
  std::istringstream issE(msgStrE);
  boost::archive::binary_iarchive inputArchiveE(issE);
  inputArchiveE >> fEventHeader;
  LOG(debug) << "EventHeader: " << fEventHeader[0] << " " << fEventHeader[1] << " " << fEventHeader[2] << " "
             << fEventHeader[3] << " " << fEventHeader[4];
 
  fiNDigiIn = 0;
  //  LOG(debug) << "Received message # "<<  fNumMessages
  //         << " with size " << msg->GetSize()<<" at "<< msg->GetData();
 
  //std::string msgStr(static_cast<char*>(msg->GetData()), msg->GetSize());
  std::string msgStr(static_cast<char*>(parts.At(1)->GetData()), (parts.At(1))->GetSize());
  std::istringstream iss(msgStr);
  boost::archive::binary_iarchive inputArchive(iss);
 
  std::vector<CbmTofDigi*> vdigi;
  inputArchive >> vdigi;
 
  /*  ---- for debugging ----------------
  int *pData = static_cast <int *>(msg->GetData());
  for (int iData=0; iData<msg->GetSize()/NBytes; iData++) {
    LOG(info) << Form(" ind %d, poi %p, data: 0x%08x",iData,pData,*pData++);
  }
  */
 
  //vector descriptor and data separated -> transfer of vectors does not work reliably
  //std::vector<CbmTofDigi>* vdigi = static_cast<std::vector<CbmTofDigi>*>(msg->GetData());
  //  (*vdigi).resize(fiNDigiIn);
  LOG(debug) << "vdigi vector at " << vdigi.data() << " with size " << vdigi.size();
 
  for (UInt_t iDigi = 0; iDigi < vdigi.size(); iDigi++) {
    LOG(debug) << "#" << iDigi << " " << vdigi[iDigi]->ToString();
  }
 
  /*
  const Int_t iNDigiIn=100;
  std::array<CbmTofDigi,iNDigiIn> *aTofDigi = static_cast<std::array<CbmTofDigi,iNDigiIn>*>(msg->GetData());
  for (int iDigi=0; iDigi<fiNDigiIn; iDigi++) {
    LOG(info) << "#" << iDigi << " " <<(*aTofDigi)[iDigi].ToString();
  }
 
 
  pDigiIn=static_cast<CbmTofDigi*> (msg->GetData());
  CbmTofDigi*  pDigi=pDigiIn;
  CbmTofDigi  aTofDigi[fiNDigiIn];
 
 
  for (int iDigi=0; iDigi<fiNDigiIn; iDigi++) {
  //aTofDigi[iDigi] = *pDigi++;
    aTofDigi[iDigi] = *pDigi;
    fvDigiIn[iDigi] = *pDigi;
    LOG(info) << "#" << iDigi << " at "<<pDigi<< " " <<aTofDigi[iDigi].ToString();
    // LOG(info) << "#" << iDigi << " at "<<pDigi<< " " <<pDigi->ToString();   // does not work ???
    pDigi++;
  }
  */
 
  fiNDigiIn = vdigi.size();
  fvDigiIn.resize(fiNDigiIn);
  for (int iDigi = 0; iDigi < fiNDigiIn; iDigi++) {
    fvDigiIn[iDigi] = *vdigi[iDigi];
 
    // Remap Digis for v21a_cosmicHD
    if (1) {
      int iSmType       = 8;
      int iSm           = -1;
      int iRpc          = 0;
      int iDetId        = 0;
      CbmTofDigi* pDigi = &fvDigiIn[iDigi];
 
      if (pDigi->GetType() == 6 && pDigi->GetRpc() == 1) {
        switch ((int) (pDigi->GetChannel() * 2 + pDigi->GetSide())) {
          case 62:  //800
            iSm = 0;
            break;
          case 46:  //810
            iSm = 1;
            break;
          case 22:  //820
            iSm = 2;
            break;
          case 2:  //830
            iSm = 3;
            break;
          default:;
        }
        if (iSm > -1) {
          iDetId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, 0, 0, iSmType);
          pDigi->SetAddress(iDetId);
        }
      }
    }  //remapping end
 
    delete vdigi[iDigi];
  }
  vdigi.clear();
 
  //  for( Int_t i = 0; i < fTofCalDigisColl->GetEntriesFast(); i++ ) ((CbmTofDigi*) fTofCalDigisColl->At(i))->Delete();
  // fTofCalDigisColl->Clear("C");
 
  //  for( Int_t i = 0; i < fTofHitsColl->GetEntriesFast(); i++ ) ((CbmTofHit*) fTofHitsColl->At(i))->Delete();
  fTofHitsColl->Clear("C");
  //fTofHitsColl->Delete(); // Are the elements deleted?
  //fTofHitsColl         = new TClonesArray("CbmTofHit",100);
 
  //for( Int_t i = 0; i < fTofDigiMatchColl->GetEntriesFast(); i++ ) ((CbmMatch*) fTofDigiMatchColl->At(i))->Delete();
  fTofDigiMatchColl->Clear("C");
 
  fiNbHits = 0;
  if (fNumMessages % 10000 == 0) LOG(info) << "Processed " << fNumMessages << " messages";
  if (fEventHeader.size() > 3) {
    fhPulMul->Fill((Double_t) fEventHeader[3]);
    if (fEventHeader[3] > 0) {
      //LOG(debug) << "Pulser event found, Mul "<< fEventHeader[3];
      if (!MonitorPulser()) return kFALSE;
      return kTRUE;  // separate events from pulser
    }
  }
 
  LOG(debug) << " Process msg " << fNumMessages << " at evt " << fdEvent << ", PulMode " << fiPulserMode;
 
  if (fiPulserMode > 0) {  // don't process events without valid pulser correction
    if (fvPulserTimes[fiPulDetRef][0].size() == 0) return kTRUE;
  }
 
  fdEvent++;
  fhEvDetMul->Fill((Double_t) fEventHeader[1]);
  fhEvDigiMul->Fill(fiNDigiIn);
  if (fiNDigiIn > 0) {
    if (dTstart == 0) {
      dTstart = fvDigiIn[0].GetTime() + fEventHeader[4];
      LOG(info) << "Start time set to " << dTstart / 1.E9 << " sec ";
    }
    fhEvRateIn->Fill((fvDigiIn[0].GetTime() + fEventHeader[4] - dTstart) / 1.E9);
  }
  if (!InspectRawDigis()) return kFALSE;
 
  if (fiPulserMode > 0)
    if (!ApplyPulserCorrection()) return kFALSE;
 
  if (!BuildClusters()) return kFALSE;
  //if(!MergeClusters())       return kFALSE;
 
  if (NULL != fOutRootFile) {  // CbmEvent output to root file
    fEvtHeader->SetEventTime((double) fEventHeader[4]);
    auto source = rootMgr->GetSource();
    if (source) { source->FillEventHeader(fEvtHeader); }
    //LOG(info) << "Fill WriteOutBuffer with rootMgr at " << rootMgr;
    fOutRootFile->cd();
    rootMgr->Fill();
    rootMgr->StoreWriteoutBufferData(rootMgr->GetEventTime());
    fhEvRateOut->Fill((fvDigiIn[0].GetTime() + fEventHeader[4] - dTstart) / 1.E9);
    //rootMgr->StoreAllWriteoutBufferData();
    rootMgr->DeleteOldWriteoutBufferData();
    if ((Int_t) fdEvent == fiMaxEvent) {
      rootMgr->Write();
      WriteHistograms();
      fOutRootFile->Close();
      LOG(info) << "File closed after " << fdEvent << " events. ";
      ChangeState(fair::mq::Transition::Stop);
    }
  }
  if (!FillHistos()) return kTRUE;  // event not selected for histogramming, skip sending it
  //if(!SendHits())      return kFALSE;
  if (!SendAll()) return kFALSE;
 
  return kTRUE;
}
 
/************************************************************************************/
 
bool CbmDeviceHitBuilderTof::HandleMessage(FairMQMessagePtr& msg, int /*index*/)
{
  const char* cmd    = (char*) (msg->GetData());
  const char cmda[4] = {*cmd};
  LOG(info) << "Handle message " << cmd << ", " << cmd[0];
  //LOG(info) << "Current State: " << fair::mq::StateMachine::GetCurrentStateName();
 
  // only one implemented so far "Stop"
  if (NULL != fOutRootFile) {
    LOG(info) << "Close root file " << fOutRootFile->GetName();
    fOutRootFile->cd();
    rootMgr->LastFill();
    rootMgr->Write();
    WriteHistograms();
    fOutRootFile->Write();
    fOutRootFile->Close();
  }
 
  if (strcmp(cmda, "STOP")) {
    LOG(info) << "STOP";
    /*  Invalid syntax
    ChangeState(internal_READY);
    LOG(info) << "Current State: " <<  FairMQStateMachine::GetCurrentStateName();
    ChangeState(internal_DEVICE_READY);
    LOG(info) << "Current State: " <<  FairMQStateMachine::GetCurrentStateName();
    ChangeState(internal_IDLE);
    LOG(info) << "Current State: " <<  FairMQStateMachine::GetCurrentStateName();
    ChangeState(END);
    LOG(info) << "Current State: " <<  FairMQStateMachine::GetCurrentStateName();
    */
    ChangeState(fair::mq::Transition::Stop);
    std::this_thread::sleep_for(std::chrono::milliseconds(1000));
    ChangeState(fair::mq::Transition::End);
  }
 
  return true;
}
 
/************************************************************************************/
Bool_t CbmDeviceHitBuilderTof::InitCalibParameter()
{
  // dimension and initialize calib parameter
  // Int_t iNbDet     = fDigiBdfPar->GetNbDet();
  Int_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
 
  fTotMean = 1.;
  fCalMode = -1;
 
  if (fTotMean != 0.) fdTTotMean = fTotMean;  // adjust target mean for TTT
 
  fvCPTOff.resize(iNbSmTypes);
  fvCPTotGain.resize(iNbSmTypes);
  fvCPTotOff.resize(iNbSmTypes);
  fvCPWalk.resize(iNbSmTypes);
  fvCPDelTof.resize(iNbSmTypes);
  for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
    Int_t iNbSm  = fDigiBdfPar->GetNbSm(iSmType);
    Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
    fvCPTOff[iSmType].resize(iNbSm * iNbRpc);
    fvCPTotGain[iSmType].resize(iNbSm * iNbRpc);
    fvCPTotOff[iSmType].resize(iNbSm * iNbRpc);
    fvCPWalk[iSmType].resize(iNbSm * iNbRpc);
    fvCPDelTof[iSmType].resize(iNbSm * iNbRpc);
    for (Int_t iSm = 0; iSm < iNbSm; iSm++) {
      for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
        //LOG(info)<<Form(" fvCPDelTof resize for SmT %d, R %d, B %d ",iSmType,iNbSm*iNbRpc,nbClDelTofBinX);
        fvCPDelTof[iSmType][iSm * iNbRpc + iRpc].resize(nbClDelTofBinX);
        for (Int_t iBx = 0; iBx < nbClDelTofBinX; iBx++) {
          // LOG(info)<<Form(" fvCPDelTof for SmT %d, R %d, B %d",iSmType,iSm*iNbRpc+iRpc,iBx);
          fvCPDelTof[iSmType][iSm * iNbRpc + iRpc][iBx].resize(iNSel);
          for (Int_t iSel = 0; iSel < iNSel; iSel++)
            fvCPDelTof[iSmType][iSm * iNbRpc + iRpc][iBx][iSel] = 0.;  // initialize
        }
 
        Int_t iNbChan = fDigiBdfPar->GetNbChan(iSmType, iRpc);
        fvCPTOff[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
        fvCPTotGain[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
        fvCPTotOff[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
        fvCPWalk[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
        Int_t nbSide = 2 - fDigiBdfPar->GetChanType(iSmType, iRpc);
        for (Int_t iCh = 0; iCh < iNbChan; iCh++) {
          fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh].resize(nbSide);
          fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh].resize(nbSide);
          fvCPTotOff[iSmType][iSm * iNbRpc + iRpc][iCh].resize(nbSide);
          fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh].resize(nbSide);
          for (Int_t iSide = 0; iSide < nbSide; iSide++) {
            fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][iSide]    = 0.;  //initialize
            fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] = 1.;  //initialize
            fvCPTotOff[iSmType][iSm * iNbRpc + iRpc][iCh][iSide]  = 0.;  //initialize
            fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][iSide].resize(nbClWalkBinX);
            for (Int_t iWx = 0; iWx < nbClWalkBinX; iWx++) {
              fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][iSide][iWx] = 0.;
            }
          }
        }
      }
    }
  }
  LOG(info) << "defaults set";
 
  // <= To prevent histos from being sucked in by the param file of the TRootManager!
  TFile* oldFile     = gFile;
  TDirectory* oldDir = gDirectory;
 
  if (0 < fCalMode) {
    fCalParFileName = fDigiBdfPar->GetCalibFileName();
    LOG(info) << "InitCalibParameter: read histos from "
              << "file " << fCalParFileName;
 
    // read parameter from histos
    if (fCalParFileName.IsNull()) return kTRUE;
 
    fCalParFile = new TFile(fCalParFileName, "");
    if (NULL == fCalParFile) {
      LOG(error) << "InitCalibParameter: "
                 << "file " << fCalParFileName << " does not exist!";
      ChangeState(fair::mq::Transition::Stop);
    }
    /*
    gDirectory->Print();
    fCalParFile->cd();
    fCalParFile->ls();
    */
    for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
      Int_t iNbSm     = fDigiBdfPar->GetNbSm(iSmType);
      Int_t iNbRpc    = fDigiBdfPar->GetNbRpc(iSmType);
      TProfile* hSvel = (TProfile*) gDirectory->FindObjectAny(Form("cl_SmT%01d_Svel", iSmType));
 
      // copy Histo to memory
      TDirectory* curdir = gDirectory;
      if (NULL != hSvel) {
        gDirectory->cd(oldDir->GetPath());
        //TProfile *hSvelmem =
        //(TProfile*) hSvel->Clone();
        gDirectory->cd(curdir->GetPath());
      }
      else {
        LOG(info) << "Svel histogram not found for module type " << iSmType;
      }
 
      for (Int_t iPar = 0; iPar < 4; iPar++) {
        TProfile* hFparcur = (TProfile*) gDirectory->FindObjectAny(Form("cl_SmT%01d_Fpar%1d", iSmType, iPar));
        if (NULL != hFparcur) {
          gDirectory->cd(oldDir->GetPath());
          //TProfile *hFparmem =
          //(TProfile*) hFparcur->Clone();
          gDirectory->cd(curdir->GetPath());
        }
      }
 
      for (Int_t iSm = 0; iSm < iNbSm; iSm++)
        for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
          // update default parameter
          if (NULL != hSvel) {
            Double_t Vscal = 1.;  //hSvel->GetBinContent(iSm*iNbRpc+iRpc+1);
            if (Vscal == 0.) Vscal = 1.;
            fDigiBdfPar->SetSigVel(iSmType, iSm, iRpc, fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc) * Vscal);
            LOG(info) << "Modify " << iSmType << iSm << iRpc << " Svel by " << Vscal << " to "
                      << fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc);
          }
          TH2F* htempPos_pfx =
            (TH2F*) gDirectory->FindObjectAny(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Pos_pfx", iSmType, iSm, iRpc));
          TH2F* htempTOff_pfx =
            (TH2F*) gDirectory->FindObjectAny(Form("cl_CorSmT%01d_sm%03d_rpc%03d_TOff_pfx", iSmType, iSm, iRpc));
          TH1D* htempTot_Mean =
            (TH1D*) gDirectory->FindObjectAny(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Tot_Mean", iSmType, iSm, iRpc));
          TH1D* htempTot_Off =
            (TH1D*) gDirectory->FindObjectAny(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Tot_Off", iSmType, iSm, iRpc));
          if (NULL != htempPos_pfx && NULL != htempTOff_pfx && NULL != htempTot_Mean && NULL != htempTot_Off) {
            Int_t iNbCh    = fDigiBdfPar->GetNbChan(iSmType, iRpc);
            Int_t iNbinTot = htempTot_Mean->GetNbinsX();
            for (Int_t iCh = 0; iCh < iNbCh; iCh++) {
              Double_t YMean = ((TProfile*) htempPos_pfx)->GetBinContent(iCh + 1);  //nh +1 empirical(?)
              Double_t TMean = ((TProfile*) htempTOff_pfx)->GetBinContent(iCh + 1);
              //Double_t dTYOff=YMean/fDigiBdfPar->GetSignalSpeed() ;
              Double_t dTYOff = YMean / fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc);
              fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][0] += -dTYOff + TMean;
              fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][1] += +dTYOff + TMean;
 
              for (Int_t iSide = 0; iSide < 2; iSide++) {
                Double_t TotMean = htempTot_Mean->GetBinContent(iCh * 2 + 1 + iSide);  //nh +1 empirical(?)
                if (0.001 < TotMean) { fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] *= fdTTotMean / TotMean; }
                fvCPTotOff[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] = htempTot_Off->GetBinContent(iCh * 2 + 1 + iSide);
              }
 
              if (5 == iSmType || 8 == iSmType) {
                fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][1]    = fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][0];
                fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][1] = fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][0];
                fvCPTotOff[iSmType][iSm * iNbRpc + iRpc][iCh][1]  = fvCPTotOff[iSmType][iSm * iNbRpc + iRpc][iCh][0];
              }
 
              LOG(debug) << "InitCalibParameter:"
                         << " SmT " << iSmType << " Sm " << iSm << " Rpc " << iRpc << " Ch " << iCh
                         << Form(": YMean %f, TMean %f", YMean, TMean) << " -> "
                         << Form(" %f, %f, %f, %f ", fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][0],
                                 fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][1],
                                 fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][0],
                                 fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][1])
                         << ", NbinTot " << iNbinTot;
              TH1D* htempWalk0 = (TH1D*) gDirectory->FindObjectAny(
                Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S0_Walk_px", iSmType, iSm, iRpc, iCh));
              TH1D* htempWalk1 = (TH1D*) gDirectory->FindObjectAny(
                Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S1_Walk_px", iSmType, iSm, iRpc, iCh));
              if (NULL != htempWalk0 && NULL != htempWalk1) {  // reinitialize Walk array
                LOG(debug) << "Initialize Walk correction for "
                           << Form(" SmT%01d_sm%03d_rpc%03d_Ch%03d", iSmType, iSm, iRpc, iCh);
                if (htempWalk0->GetNbinsX() != nbClWalkBinX)
                  LOG(error) << "InitCalibParameter: Inconsistent Walk histograms";
                for (Int_t iBin = 0; iBin < nbClWalkBinX; iBin++) {
                  fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][0][iBin] = htempWalk0->GetBinContent(iBin + 1);
                  fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][1][iBin] = htempWalk1->GetBinContent(iBin + 1);
                  //LOG(debug)<<Form(" SmT%01d_sm%03d_rpc%03d_Ch%03d bin %d walk %f ",iSmType, iSm, iRpc, iCh, iBin,
                  //                 fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][0][iBin]);
                  if (5 == iSmType || 8 == iSmType) {  // Pad structure
                    fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][1][iBin] =
                      fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][0][iBin];
                  }
                }
              }
            }
          }
          else {
            LOG(warn) << " Calibration histos " << Form("cl_SmT%01d_sm%03d_rpc%03d_XXX", iSmType, iSm, iRpc)
                      << " not found. ";
          }
          for (Int_t iSel = 0; iSel < iNSel; iSel++) {
            TH1D* htmpDelTof = (TH1D*) gDirectory->FindObjectAny(
              Form("cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof", iSmType, iSm, iRpc, iSel));
            if (NULL == htmpDelTof) {
              LOG(debug) << " Histos " << Form("cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof", iSmType, iSm, iRpc, iSel)
                         << " not found. ";
              continue;
            }
            LOG(debug) << " Load DelTof from histos "
                       << Form("cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof", iSmType, iSm, iRpc, iSel) << ".";
            for (Int_t iBx = 0; iBx < nbClDelTofBinX; iBx++) {
              fvCPDelTof[iSmType][iSm * iNbRpc + iRpc][iBx][iSel] += htmpDelTof->GetBinContent(iBx + 1);
            }
 
            // copy Histo to memory
            // TDirectory * curdir = gDirectory;
            gDirectory->cd(oldDir->GetPath());
            TH1D* h1DelTof =
              (TH1D*) htmpDelTof->Clone(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof", iSmType, iSm, iRpc, iSel));
 
            LOG(debug) << " copy histo " << h1DelTof->GetName() << " to directory " << oldDir->GetName();
            gDirectory->cd(curdir->GetPath());
          }
        }
    }
  }
  //   fCalParFile->Delete();
  // <= To prevent histos from being sucked in by the param file of the TRootManager!
  gFile      = oldFile;
  gDirectory = oldDir;
 
  LOG(info) << "InitCalibParameter: initialization done";
  return kTRUE;
}
 
static TH2F* fhBucDigiCor = NULL;
/************************************************************************************/
// Histogram definitions
void CbmDeviceHitBuilderTof::CreateHistograms()
{
  TDirectory* oldir = gDirectory;  // <= To prevent histos from being sucked in by the param file of the TRootManager!
  gROOT->cd();                     // <= To prevent histos from being sucked in by the param file of the TRootManager !
  // process event header info
  fhEvDetMul  = new TH1F("hEvDetMul", "Detector multiplicity of Selector; Mul", 50, 0, 50);
  fhEvDigiMul = new TH1F("hEvDigiMul", "Digi multiplicity of Selector; Mul", 1000, 0, 10000);
  fhEvRateIn  = new TH1F("hEvRateIn", "Incoming Event rate; time (s); rate (Hz)", 10000, 0, 10000);
  fhEvRateOut = new TH1F("hEvRateOut", "Outgoing Event rate; time (s); rate (Hz)", 10000, 0, 10000);
  fhPulMul    = new TH1F("hPulMul", "Pulser multiplicity; Mul", 110, 0, 110);
  fhDigiTdif  = new TH1F("hDigiTdif", "Digi time difference; #Deltat ns)", 8000, -20, 20);
 
  Int_t iNDet = 0;
  dTmax *= 2;  // cover full range of selector
 
  for (Int_t iModTyp = 0; iModTyp < 10; iModTyp++)
    iNDet += fDigiBdfPar->GetNbSm(iModTyp) * fDigiBdfPar->GetNbRpc(iModTyp);
 
  fhPulserTimesRaw = new TH2F(Form("hPulserTimesRaw"), Form("Pulser Times uncorrected; Det# []; t - t0 [ns]"),
                              iNDet * 2, 0, iNDet * 2, 999, -dTmax, dTmax);
 
  fhPulserTimeRawEvo.resize(iNDet * 2);
  for (Int_t iDetSide = 0; iDetSide < iNDet * 2; iDetSide++) {
    fhPulserTimeRawEvo[iDetSide] = new TProfile(
      Form("hPulserTimeRawEvo_%d", iDetSide),
      Form("Raw Pulser TimeEvolution  %d; PulserEvent# ; DeltaT [ns] ", iDetSide), 1000, 0., 1.E5, -dTmax, dTmax);
  }
 
  fhPulserTimesCor = new TH2F(Form("hPulserTimesCor"), Form("Pulser Times corrected; Det# []; t - t0 [ns]"), iNDet * 2,
                              0, iNDet * 2, 999, -10., 10.);
 
  fhDigiTimesRaw = new TH2F(Form("hDigiTimesRaw"), Form("Digi Times uncorrected; Det# []; t - t0 [ns]"), iNDet * 2, 0,
                            iNDet * 2, 999, -dTmax, dTmax);
 
  fhDigiTimesCor = new TH2F(Form("hDigiTimesCor"), Form("Digi Times corrected; Det# []; t - t0 [ns]"), iNDet * 2, 0,
                            iNDet * 2, 999, -dTmax, dTmax);
 
  fhBucDigiCor = new TH2F(Form("hBucDigiCor"), Form("Buc Digi Correlation; ch1 []; ch2 []"), 128, 0, 128, 128, 0, 128);
 
  Int_t iNbDet = fDigiBdfPar->GetNbDet();
  fDetIdIndexMap.clear();
  fhRpcDigiTot.resize(iNbDet);
  fhRpcDigiCor.resize(iNbDet);
  fviDetId.resize(iNbDet);
  for (Int_t iDetIndx = 0; iDetIndx < iNbDet; iDetIndx++) {
    Int_t iUniqueId           = fDigiBdfPar->GetDetUId(iDetIndx);
    fDetIdIndexMap[iUniqueId] = iDetIndx;
    fviDetId[iDetIndx]        = iUniqueId;
    Int_t iSmType             = CbmTofAddress::GetSmType(iUniqueId);
    Int_t iSmId               = CbmTofAddress::GetSmId(iUniqueId);
    Int_t iRpcId              = CbmTofAddress::GetRpcId(iUniqueId);
    LOG(info) << "Index map entry " << iDetIndx << ": TSR " << iSmType << iSmId << iRpcId
              << Form(", addr 0x%08x", iUniqueId);
 
    fhRpcDigiTot[iDetIndx] = new TH2F(
      Form("hDigiTot_SmT%01d_sm%03d_rpc%03d", iSmType, iSmId, iRpcId),
      Form("Digi Tot of Rpc #%03d in Sm %03d of type %d; digi 0; digi 1", iRpcId, iSmId, iSmType),
      2 * fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, 2 * fDigiBdfPar->GetNbChan(iSmType, iRpcId), 256, 0, 256);
 
    fhRpcDigiCor[iDetIndx] =
      new TH2F(Form("cl_SmT%01d_sm%03d_rpc%03d_DigiCor", iSmType, iSmId, iRpcId),
               Form("Digi Correlation of Rpc #%03d in Sm %03d of type %d; digi 0; digi 1", iRpcId, iSmId, iSmType),
               fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId),
               fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId));
  }
 
  if (0 == fiMode) return;  // no cluster histograms needed
 
  // Sm related distributions
  fhSmCluPosition.resize(fDigiBdfPar->GetNbSmTypes());
  fhSmCluTOff.resize(fDigiBdfPar->GetNbSmTypes());
  fhSmCluSvel.resize(fDigiBdfPar->GetNbSmTypes());
  fhSmCluFpar.resize(fDigiBdfPar->GetNbSmTypes());
  fhTSmCluPosition.resize(fDigiBdfPar->GetNbSmTypes());
  fhTSmCluTOff.resize(fDigiBdfPar->GetNbSmTypes());
  fhTSmCluTRun.resize(fDigiBdfPar->GetNbSmTypes());
 
  for (Int_t iS = 0; iS < fDigiBdfPar->GetNbSmTypes(); iS++) {
    Double_t YSCAL = 50.;
    if (fPosYMaxScal != 0.) YSCAL = fPosYMaxScal;
 
    Int_t iUCellId(0);
    fChannelInfo = NULL;
 
    // Cover the case that the geometry file does not contain the module
    // indexed with 0 of a certain module type BUT modules with higher indices.
    for (Int_t iSM = 0; iSM < fDigiBdfPar->GetNbSm(iS); iSM++) {
      iUCellId     = CbmTofAddress::GetUniqueAddress(iSM, 0, 0, 0, iS);
      fChannelInfo = fDigiPar->GetCell(iUCellId);
 
      // Retrieve geometry information from the first module of a certain
      // module type that is found in the geometry file.
      if (NULL != fChannelInfo) { break; }
    }
 
    if (NULL == fChannelInfo) {
      LOG(warn) << "No DigiPar for SmType " << Form("%d, 0x%08x", iS, iUCellId);
      continue;
    }
    Double_t YDMAX = TMath::Max(2., fChannelInfo->GetSizey()) * YSCAL;
 
    fhSmCluPosition[iS] =
      new TH2F(Form("cl_SmT%01d_Pos", iS), Form("Clu position of SmType %d; Sm+Rpc# []; ypos [cm]", iS),
               fDigiBdfPar->GetNbSm(iS) * fDigiBdfPar->GetNbRpc(iS), 0,
               fDigiBdfPar->GetNbSm(iS) * fDigiBdfPar->GetNbRpc(iS), 99, -YDMAX, YDMAX);
 
    Double_t TSumMax = 1.E3;
    if (fTRefDifMax != 0.) TSumMax = fTRefDifMax;
    fhSmCluTOff[iS] =
      new TH2F(Form("cl_SmT%01d_TOff", iS), Form("Clu TimeZero in SmType %d; Sm+Rpc# []; TOff [ns]", iS),
               fDigiBdfPar->GetNbSm(iS) * fDigiBdfPar->GetNbRpc(iS), 0,
               fDigiBdfPar->GetNbSm(iS) * fDigiBdfPar->GetNbRpc(iS), 99, -TSumMax, TSumMax);
 
    TProfile* hSvelcur = (TProfile*) gDirectory->FindObjectAny(Form("cl_SmT%01d_Svel", iS));
    if (NULL == hSvelcur) {
      fhSmCluSvel[iS] =
        new TProfile(Form("cl_SmT%01d_Svel", iS), Form("Clu Svel in SmType %d; Sm+Rpc# []; v/v_{nominal} ", iS),
                     fDigiBdfPar->GetNbSm(iS) * fDigiBdfPar->GetNbRpc(iS), 0,
                     fDigiBdfPar->GetNbSm(iS) * fDigiBdfPar->GetNbRpc(iS), 0.8, 1.2);
      fhSmCluSvel[iS]->Sumw2();
    }
    else
      fhSmCluSvel[iS] = (TProfile*) hSvelcur->Clone();
 
    fhSmCluFpar[iS].resize(4);
    for (Int_t iPar = 0; iPar < 4; iPar++) {
      TProfile* hFparcur = (TProfile*) gDirectory->FindObjectAny(Form("cl_SmT%01d_Fpar%1d", iS, iPar));
      if (NULL == hFparcur) {
        LOG(info) << "Histo " << Form("cl_SmT%01d_Fpar%1d", iS, iPar) << " not found, recreate ...";
        fhSmCluFpar[iS][iPar] = new TProfile(Form("cl_SmT%01d_Fpar%1d", iS, iPar),
                                             Form("Clu Fpar %d in SmType %d; Sm+Rpc# []; value ", iPar, iS),
                                             fDigiBdfPar->GetNbSm(iS) * fDigiBdfPar->GetNbRpc(iS), 0,
                                             fDigiBdfPar->GetNbSm(iS) * fDigiBdfPar->GetNbRpc(iS), -dTmax, dTmax);
      }
      else
        fhSmCluFpar[iS][iPar] = (TProfile*) hFparcur->Clone();
    }
 
    fhTSmCluPosition[iS].resize(iNSel);
    fhTSmCluTOff[iS].resize(iNSel);
    fhTSmCluTRun[iS].resize(iNSel);
    for (Int_t iSel = 0; iSel < iNSel; iSel++) {  // Loop over selectors
      fhTSmCluPosition[iS][iSel] = new TH2F(Form("cl_TSmT%01d_Sel%02d_Pos", iS, iSel),
                                            Form("Clu position of SmType %d under Selector %02d; Sm+Rpc# "
                                                 "[]; ypos [cm]",
                                                 iS, iSel),
                                            fDigiBdfPar->GetNbSm(iS) * fDigiBdfPar->GetNbRpc(iS), 0,
                                            fDigiBdfPar->GetNbSm(iS) * fDigiBdfPar->GetNbRpc(iS), 99, -YDMAX, YDMAX);
      fhTSmCluTOff[iS][iSel]     = new TH2F(Form("cl_TSmT%01d_Sel%02d_TOff", iS, iSel),
                                        Form("Clu TimeZero in SmType %d under Selector %02d; Sm+Rpc# "
                                             "[]; TOff [ns]",
                                             iS, iSel),
                                        fDigiBdfPar->GetNbSm(iS) * fDigiBdfPar->GetNbRpc(iS), 0,
                                        fDigiBdfPar->GetNbSm(iS) * fDigiBdfPar->GetNbRpc(iS), 99, -TSumMax, TSumMax);
      fhTSmCluTRun[iS][iSel]     = new TH2F(Form("cl_TSmT%01d_Sel%02d_TRun", iS, iSel),
                                        Form("Clu TimeZero in SmType %d under Selector %02d; Event# "
                                             "[]; TMean [ns]",
                                             iS, iSel),
                                        100, 0, MaxNbEvent, 99, -TSumMax, TSumMax);
    }
  }
 
  // RPC related distributions
  LOG(info) << " Define Clusterizer histos for " << iNbDet << " detectors ";
 
  fhRpcCluMul.resize(iNbDet);
  fhRpcCluRate.resize(iNbDet);
  fhRpcCluPosition.resize(iNbDet);
  fhRpcCluDelPos.resize(iNbDet);
  fhRpcCluDelMatPos.resize(iNbDet);
  fhRpcCluTOff.resize(iNbDet);
  fhRpcCluDelTOff.resize(iNbDet);
  fhRpcCluDelMatTOff.resize(iNbDet);
  fhRpcCluTrms.resize(iNbDet);
  fhRpcCluTot.resize(iNbDet);
  fhRpcCluSize.resize(iNbDet);
  fhRpcCluAvWalk.resize(iNbDet);
  fhRpcCluAvLnWalk.resize(iNbDet);
  fhRpcCluWalk.resize(iNbDet);
  fhRpcDTLastHits.resize(iNbDet);
  fhRpcDTLastHits_Tot.resize(iNbDet);
  fhRpcDTLastHits_CluSize.resize(iNbDet);
 
 
  if (fTotMax != 0.) fdTOTMax = fTotMax;
  if (fTotMin != 0.) fdTOTMin = fTotMin;
 
  for (Int_t iDetIndx = 0; iDetIndx < iNbDet; iDetIndx++) {
    Int_t iUniqueId = fDigiBdfPar->GetDetUId(iDetIndx);
    Int_t iSmType   = CbmTofAddress::GetSmType(iUniqueId);
    Int_t iSmId     = CbmTofAddress::GetSmId(iUniqueId);
    Int_t iRpcId    = CbmTofAddress::GetRpcId(iUniqueId);
    Int_t iUCellId  = CbmTofAddress::GetUniqueAddress(iSmId, iRpcId, 0, 0, iSmType);
    fChannelInfo    = fDigiPar->GetCell(iUCellId);
    if (NULL == fChannelInfo) {
      LOG(warn) << "No DigiPar for Det " << Form("0x%08x", iUCellId);
      continue;
    }
    LOG(info) << "DetIndx " << iDetIndx << ", SmType " << iSmType << ", SmId " << iSmId << ", RpcId " << iRpcId
              << " => UniqueId " << Form("(0x%08x, 0x%08x)", iUniqueId, iUCellId) << ", dx " << fChannelInfo->GetSizex()
              << ", dy " << fChannelInfo->GetSizey() << Form(" ChPoi: %p ", fChannelInfo) << ", nbCh "
              << fDigiBdfPar->GetNbChan(iSmType, 0);
 
    // check access to all channel infos
    for (Int_t iCh = 0; iCh < fDigiBdfPar->GetNbChan(iSmType, iRpcId); iCh++) {
      Int_t iCCellId = CbmTofAddress::GetUniqueAddress(iSmId, iRpcId, iCh, 0, iSmType);
      fChannelInfo   = fDigiPar->GetCell(iCCellId);
      if (NULL == fChannelInfo) LOG(warn) << Form("missing ChannelInfo for ch %d addr 0x%08x", iCh, iCCellId);
    }
 
    fChannelInfo = fDigiPar->GetCell(iUCellId);
 
    fhRpcCluMul[iDetIndx] =
      new TH1F(Form("cl_SmT%01d_sm%03d_rpc%03d_Mul", iSmType, iSmId, iRpcId),
               Form("Clu multiplicity of Rpc #%03d in Sm %03d of type %d; M []; Cnts", iRpcId, iSmId, iSmType),
               2. + 2. * fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, 2. + 2. * fDigiBdfPar->GetNbChan(iSmType, iRpcId));
 
    fhRpcCluRate[iDetIndx] =
      new TH1F(Form("cl_SmT%01d_sm%03d_rpc%03d_rate", iSmType, iSmId, iRpcId),
               Form("Clu rate of Rpc #%03d in Sm %03d of type %d; Time (s); Rate (Hz)", iRpcId, iSmId, iSmType), 3600.,
               0., 3600.);
 
    fhRpcDTLastHits[iDetIndx] = new TH1F(Form("cl_SmT%01d_sm%03d_rpc%03d_DTLastHits", iSmType, iSmId, iRpcId),
                                         Form("Clu #DeltaT to last hits  of Rpc #%03d in Sm %03d of type %d; log( "
                                              "#DeltaT (ns)); counts",
                                              iRpcId, iSmId, iSmType),
                                         100., 0., 10.);
 
    fhRpcDTLastHits_Tot[iDetIndx] = new TH2F(Form("cl_SmT%01d_sm%03d_rpc%03d_Tot_DTLH", iSmType, iSmId, iRpcId),
                                             Form("Clu Tot of Rpc #%03d in Sm %03d of type %d; log(#DeltaT (ns)); TOT "
                                                  "[ns]",
                                                  iRpcId, iSmId, iSmType),
                                             100, 0., 10., 100, fdTOTMin, 4. * fdTOTMax);
 
    fhRpcDTLastHits_CluSize[iDetIndx] = new TH2F(Form("cl_SmT%01d_sm%03d_rpc%03d_CluSize_DTLH", iSmType, iSmId, iRpcId),
                                                 Form("Clu Size of Rpc #%03d in Sm %03d of type %d; log(#DeltaT (ns)); "
                                                      "CluSize []",
                                                      iRpcId, iSmId, iSmType),
                                                 100, 0., 10., 16, 0.5, 16.5);
 
    Double_t YSCAL = 50.;
    if (fPosYMaxScal != 0.) YSCAL = fPosYMaxScal;
    Double_t YDMAX = TMath::Max(2., fChannelInfo->GetSizey()) * YSCAL;
    fhRpcCluPosition[iDetIndx] =
      new TH2F(Form("cl_SmT%01d_sm%03d_rpc%03d_Pos", iSmType, iSmId, iRpcId),
               Form("Clu position of Rpc #%03d in Sm %03d of type %d; Strip []; ypos [cm]", iRpcId, iSmId, iSmType),
               fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 99, -YDMAX, YDMAX);
 
    fhRpcCluDelPos[iDetIndx] =
      new TH2F(Form("cl_SmT%01d_sm%03d_rpc%03d_DelPos", iSmType, iSmId, iRpcId),
               Form("Clu position difference of Rpc #%03d in Sm %03d of type "
                    "%d; Strip []; #Deltaypos(clu) [cm]",
                    iRpcId, iSmId, iSmType),
               fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 99, -10., 10.);
 
    fhRpcCluDelMatPos[iDetIndx] =
      new TH2F(Form("cl_SmT%01d_sm%03d_rpc%03d_DelMatPos", iSmType, iSmId, iRpcId),
               Form("Matched Clu position difference of Rpc #%03d in Sm %03d of type "
                    "%d; Strip []; #Deltaypos(mat) [cm]",
                    iRpcId, iSmId, iSmType),
               fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 99, -5., 5.);
 
    Double_t TSumMax = 1.E3;
    if (fTRefDifMax != 0.) TSumMax = fTRefDifMax;
    fhRpcCluTOff[iDetIndx] = new TH2F(
      Form("cl_SmT%01d_sm%03d_rpc%03d_TOff", iSmType, iSmId, iRpcId),
      Form("Clu TimeZero of Rpc #%03d in Sm %03d of type %d; Strip []; TOff [ns]", iRpcId, iSmId, iSmType),
      fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 99, -TSumMax, TSumMax);
 
    fhRpcCluDelTOff[iDetIndx] =
      new TH2F(Form("cl_SmT%01d_sm%03d_rpc%03d_DelTOff", iSmType, iSmId, iRpcId),
               Form("Clu TimeZero Difference of Rpc #%03d in Sm %03d of type %d; Strip "
                    "[]; #DeltaTOff(clu) [ns]",
                    iRpcId, iSmId, iSmType),
               fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 99, -0.6, 0.6);
 
    fhRpcCluDelMatTOff[iDetIndx] =
      new TH2F(Form("cl_SmT%01d_sm%03d_rpc%03d_DelMatTOff", iSmType, iSmId, iRpcId),
               Form("Clu TimeZero Difference of Rpc #%03d in Sm %03d of type %d; Strip "
                    "[]; #DeltaTOff(mat) [ns]",
                    iRpcId, iSmId, iSmType),
               fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 99, -0.6, 0.6);
 
    fhRpcCluTrms[iDetIndx] =
      new TH2F(Form("cl_SmT%01d_sm%03d_rpc%03d_Trms", iSmType, iSmId, iRpcId),
               Form("Clu Time RMS of Rpc #%03d in Sm %03d of type %d; Strip []; Trms [ns]", iRpcId, iSmId, iSmType),
               fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 99, 0., 0.5);
 
    fhRpcCluTot[iDetIndx] =
      new TH2F(Form("cl_SmT%01d_sm%03d_rpc%03d_Tot", iSmType, iSmId, iRpcId),
               Form("Clu Tot of Rpc #%03d in Sm %03d of type %d; StripSide []; TOT [ns]", iRpcId, iSmId, iSmType),
               2 * fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, 2 * fDigiBdfPar->GetNbChan(iSmType, iRpcId), 100,
               fdTOTMin, fdTOTMax);
 
    fhRpcCluSize[iDetIndx] =
      new TH2F(Form("cl_SmT%01d_sm%03d_rpc%03d_Size", iSmType, iSmId, iRpcId),
               Form("Clu size of Rpc #%03d in Sm %03d of type %d; Strip []; size [strips]", iRpcId, iSmId, iSmType),
               fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 16, 0.5, 16.5);
 
    // Walk histos
    fhRpcCluAvWalk[iDetIndx] = new TH2F(Form("cl_SmT%01d_sm%03d_rpc%03d_AvWalk", iSmType, iSmId, iRpcId),
                                        Form("Walk in SmT%01d_sm%03d_rpc%03d_AvWalk", iSmType, iSmId, iRpcId),
                                        nbClWalkBinX, fdTOTMin, fdTOTMax, nbClWalkBinY, -TSumMax, TSumMax);
 
    fhRpcCluAvLnWalk[iDetIndx] =
      new TH2D(Form("cl_SmT%01d_sm%03d_rpc%03d_AvLnWalk", iSmType, iSmId, iRpcId),
               Form("Walk in SmT%01d_sm%03d_rpc%03d_AvLnWalk", iSmType, iSmId, iRpcId), nbClWalkBinX,
               TMath::Log10(fdTOTMax / 50.), TMath::Log10(fdTOTMax), nbClWalkBinY, -TSumMax, TSumMax);
 
    fhRpcCluWalk[iDetIndx].resize(fDigiBdfPar->GetNbChan(iSmType, iRpcId));
    for (Int_t iCh = 0; iCh < fDigiBdfPar->GetNbChan(iSmType, iRpcId); iCh++) {
      fhRpcCluWalk[iDetIndx][iCh].resize(2);
      for (Int_t iSide = 0; iSide < 2; iSide++) {
        fhRpcCluWalk[iDetIndx][iCh][iSide] =
          new TH2D(Form("cl_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Walk", iSmType, iSmId, iRpcId, iCh, iSide),
                   Form("Walk in SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Walk", iSmType, iSmId, iRpcId, iCh, iSide),
                   nbClWalkBinX, fdTOTMin, fdTOTMax, nbClWalkBinY, -TSumMax, TSumMax);
      }
      /*
         (fhRpcCluWalk[iDetIndx]).push_back( hTemp );
       */
    }
  }
 
  // Trigger selected histograms
  if (0 < iNSel) {
 
    fhSeldT.resize(iNSel);
    for (Int_t iSel = 0; iSel < iNSel; iSel++) {
      fhSeldT[iSel] = new TH2F(Form("cl_dt_Sel%02d", iSel), Form("Selector time %02d; dT [ns]", iSel), 99,
                               -fdDelTofMax * 10., fdDelTofMax * 10., 16, -0.5, 15.5);
    }
 
    fhTRpcCluMul.resize(iNbDet);
    fhTRpcCluPosition.resize(iNbDet);
    fhTRpcCluTOff.resize(iNbDet);
    fhTRpcCluTot.resize(iNbDet);
    fhTRpcCluSize.resize(iNbDet);
    fhTRpcCluAvWalk.resize(iNbDet);
    fhTRpcCluDelTof.resize(iNbDet);
    fhTRpcCludXdY.resize(iNbDet);
    fhTRpcCluWalk.resize(iNbDet);
    fhTRpcCluTOffDTLastHits.resize(iNbDet);
    fhTRpcCluTotDTLastHits.resize(iNbDet);
    fhTRpcCluSizeDTLastHits.resize(iNbDet);
    fhTRpcCluMemMulDTLastHits.resize(iNbDet);
 
    for (Int_t iDetIndx = 0; iDetIndx < iNbDet; iDetIndx++) {
      Int_t iUniqueId = fDigiBdfPar->GetDetUId(iDetIndx);
      Int_t iSmType   = CbmTofAddress::GetSmType(iUniqueId);
      Int_t iSmId     = CbmTofAddress::GetSmId(iUniqueId);
      Int_t iRpcId    = CbmTofAddress::GetRpcId(iUniqueId);
      Int_t iUCellId  = CbmTofAddress::GetUniqueAddress(iSmId, iRpcId, 0, 0, iSmType);
      fChannelInfo    = fDigiPar->GetCell(iUCellId);
      if (NULL == fChannelInfo) {
        LOG(warn) << "No DigiPar for Det " << Form("0x%08x", iUCellId);
        continue;
      }
      LOG(debug) << "DetIndx " << iDetIndx << ", SmType " << iSmType << ", SmId " << iSmId << ", RpcId " << iRpcId
                 << " => UniqueId " << Form("(0x%08x, 0x%08x)", iUniqueId, iUCellId) << ", dx "
                 << fChannelInfo->GetSizex() << ", dy " << fChannelInfo->GetSizey() << Form(" poi: 0x%p ", fChannelInfo)
                 << ", nbCh " << fDigiBdfPar->GetNbChan(iSmType, iRpcId);
 
      fhTRpcCluMul[iDetIndx].resize(iNSel);
      fhTRpcCluPosition[iDetIndx].resize(iNSel);
      fhTRpcCluTOff[iDetIndx].resize(iNSel);
      fhTRpcCluTot[iDetIndx].resize(iNSel);
      fhTRpcCluSize[iDetIndx].resize(iNSel);
      fhTRpcCluAvWalk[iDetIndx].resize(iNSel);
      fhTRpcCluDelTof[iDetIndx].resize(iNSel);
      fhTRpcCludXdY[iDetIndx].resize(iNSel);
      fhTRpcCluWalk[iDetIndx].resize(iNSel);
      fhTRpcCluTOffDTLastHits[iDetIndx].resize(iNSel);
      fhTRpcCluTotDTLastHits[iDetIndx].resize(iNSel);
      fhTRpcCluSizeDTLastHits[iDetIndx].resize(iNSel);
      fhTRpcCluMemMulDTLastHits[iDetIndx].resize(iNSel);
 
      for (Int_t iSel = 0; iSel < iNSel; iSel++) {
        fhTRpcCluMul[iDetIndx][iSel] =
          new TH1F(Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_Mul", iSmType, iSmId, iRpcId, iSel),
                   Form("Clu multiplicity of Rpc #%03d in Sm %03d of type %d "
                        "under Selector %02d; M []; cnts",
                        iRpcId, iSmId, iSmType, iSel),
                   fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0., fDigiBdfPar->GetNbChan(iSmType, iRpcId));
 
        if (NULL == fhTRpcCluMul[iDetIndx][iSel]) LOG(error) << " Histo not generated !";
 
        Double_t YSCAL = 50.;
        if (fPosYMaxScal != 0.) YSCAL = fPosYMaxScal;
        Double_t YDMAX                    = TMath::Max(2., fChannelInfo->GetSizey()) * YSCAL;
        fhTRpcCluPosition[iDetIndx][iSel] = new TH2F(
          Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_Pos", iSmType, iSmId, iRpcId, iSel),
          Form("Clu position of Rpc #%03d in Sm %03d of type %d under "
               "Selector %02d; Strip []; ypos [cm]",
               iRpcId, iSmId, iSmType, iSel),
          fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 100, -YDMAX, YDMAX);
 
        Double_t TSumMax = 1.E4;
        if (fTRefDifMax != 0.) TSumMax = fTRefDifMax;
        fhTRpcCluTOff[iDetIndx][iSel] = new TH2F(
          Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_TOff", iSmType, iSmId, iRpcId, iSel),
          Form("Clu TimeZero of Rpc #%03d in Sm %03d of type %d under "
               "Selector %02d; Strip []; TOff [ns]",
               iRpcId, iSmId, iSmType, iSel),
          fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 99, -TSumMax, TSumMax);
 
        if (fTotMax != 0.) fdTOTMax = fTotMax;
        fhTRpcCluTot[iDetIndx][iSel] =
          new TH2F(Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_Tot", iSmType, iSmId, iRpcId, iSel),
                   Form("Clu Tot of Rpc #%03d in Sm %03d of type %d under "
                        "Selector %02d; StripSide []; TOT [ns]",
                        iRpcId, iSmId, iSmType, iSel),
                   fDigiBdfPar->GetNbChan(iSmType, iRpcId) * 2, 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId) * 2, 100,
                   fdTOTMin, fdTOTMax);
 
        fhTRpcCluSize[iDetIndx][iSel] =
          new TH2F(Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_Size", iSmType, iSmId, iRpcId, iSel),
                   Form("Clu size of Rpc #%03d in Sm %03d of type %d under "
                        "Selector %02d; Strip []; size [strips]",
                        iRpcId, iSmId, iSmType, iSel),
                   fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 16, 0.5, 16.5);
 
        // Walk histos
        fhTRpcCluAvWalk[iDetIndx][iSel] =
          new TH2D(Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_AvWalk", iSmType, iSmId, iRpcId, iSel),
                   Form("Walk in SmT%01d_sm%03d_rpc%03d_Sel%02d_AvWalk; TOT; T-TSel", iSmType, iSmId, iRpcId, iSel),
                   nbClWalkBinX, fdTOTMin, fdTOTMax, nbClWalkBinY, -TSumMax, TSumMax);
 
        // Tof Histos
        fhTRpcCluDelTof[iDetIndx][iSel] =
          new TH2F(Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof", iSmType, iSmId, iRpcId, iSel),
                   Form("SmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof; TRef-TSel; T-TSel", iSmType, iSmId, iRpcId, iSel),
                   nbClDelTofBinX, -fdDelTofMax, fdDelTofMax, nbClDelTofBinY, -TSumMax, TSumMax);
 
        // Position deviation histos
        fhTRpcCludXdY[iDetIndx][iSel] =
          new TH2F(Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_dXdY", iSmType, iSmId, iRpcId, iSel),
                   Form("SmT%01d_sm%03d_rpc%03d_Sel%02d_dXdY; #Delta x [cm]; "
                        "#Delta y [cm];",
                        iSmType, iSmId, iRpcId, iSel),
                   nbCldXdYBinX, -dXdYMax, dXdYMax, nbCldXdYBinY, -dXdYMax, dXdYMax);
 
        fhTRpcCluWalk[iDetIndx][iSel].resize(fDigiBdfPar->GetNbChan(iSmType, iRpcId));
        for (Int_t iCh = 0; iCh < fDigiBdfPar->GetNbChan(iSmType, iRpcId); iCh++) {
          fhTRpcCluWalk[iDetIndx][iSel][iCh].resize(2);
          for (Int_t iSide = 0; iSide < 2; iSide++) {
            fhTRpcCluWalk[iDetIndx][iSel][iCh][iSide] = new TH2D(
              Form("cl_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Sel%02d_Walk", iSmType, iSmId, iRpcId, iCh, iSide, iSel),
              Form("Walk in SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Sel%02d_Walk", iSmType, iSmId, iRpcId, iCh, iSide,
                   iSel),
              nbClWalkBinX, fdTOTMin, fdTOTMax, nbClWalkBinY, -TSumMax, TSumMax);
          }
        }
 
        fhTRpcCluTOffDTLastHits[iDetIndx][iSel] =
          new TH2F(Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_TOff_DTLH", iSmType, iSmId, iRpcId, iSel),
                   Form("Clu TimeZero of Rpc #%03d in Sm %03d of type %d under "
                        "Selector %02d; log(#DeltaT (ns)); TOff [ns]",
                        iRpcId, iSmId, iSmType, iSel),
                   100, 0., 10., 99, -TSumMax, TSumMax);
 
        fhTRpcCluTotDTLastHits[iDetIndx][iSel] =
          new TH2F(Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_Tot_DTLH", iSmType, iSmId, iRpcId, iSel),
                   Form("Clu Tot of Rpc #%03d in Sm %03d of type %d under "
                        "Selector %02d; log(#DeltaT (ns)); TOT [ns]",
                        iRpcId, iSmId, iSmType, iSel),
                   100, 0., 10., 100, fdTOTMin, fdTOTMax);
 
        fhTRpcCluSizeDTLastHits[iDetIndx][iSel] =
          new TH2F(Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_Size_DTLH", iSmType, iSmId, iRpcId, iSel),
                   Form("Clu size of Rpc #%03d in Sm %03d of type %d under "
                        "Selector %02d; log(#DeltaT (ns)); size [strips]",
                        iRpcId, iSmId, iSmType, iSel),
                   100, 0., 10., 10, 0.5, 10.5);
 
        fhTRpcCluMemMulDTLastHits[iDetIndx][iSel] =
          new TH2F(Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_MemMul_DTLH", iSmType, iSmId, iRpcId, iSel),
                   Form("Clu Memorized Multiplicity of Rpc #%03d in Sm %03d of type %d "
                        "under Selector %02d; log(#DeltaT (ns)); TOff [ns]",
                        iRpcId, iSmId, iSmType, iSel),
                   100, 0., 10., 10, 0, 10);
      }
    }
  }
 
  gDirectory->cd(oldir->GetPath());  // <= To prevent histos from being sucked in by the param file of the TRootManager!
 
  return;
}
/************************************************************************************/
void CbmDeviceHitBuilderTof::WriteHistograms()
{
  TList* tHistoList(NULL);
  tHistoList = gROOT->GetList();
 
  TIter next(tHistoList);
  // Write histogramms to the file
  fOutRootFile->cd();
  {
    TH1* h;
    TObject* obj;
    while ((obj = (TObject*) next())) {
      if (obj->InheritsFrom(TH1::Class())) {
        h = (TH1*) obj;
        //         cout << "Write histo " << h->GetTitle() << endl;
        h->Write();
      }
    }
  }
}
 
/************************************************************************************/
Bool_t CbmDeviceHitBuilderTof::BuildClusters()
{
  fiNevtBuild++;
  if (!CalibRawDigis()) return kFALSE;
  if (0 == fiMode) return kTRUE;  // no customer yet
  if (!FillDigiStor()) return kFALSE;
  if (!BuildHits()) return kFALSE;
  return kTRUE;
}
 
/************************************************************************************/
Bool_t CbmDeviceHitBuilderTof::InspectRawDigis()
{
  Int_t iNbTofDigi = fiNDigiIn;
  pRef             = NULL;
  for (Int_t iDigInd = 0; iDigInd < fiNDigiIn; iDigInd++) {
    CbmTofDigi* pDigi = &fvDigiIn[iDigInd];
    //LOG(debug)<<iDigInd<<" "<<pDigi;
    Int_t iAddr    = pDigi->GetAddress() & DetMask;
    Int_t iDetIndx = fDigiBdfPar->GetDetInd(iAddr);
    /*
    LOG(info)<<iDigInd
              //<<Form(" Address : 0x%08x ",pDigi->GetAddress())
              <<Form(" TSRCS %d%d%d%02d%d ", (int)pDigi->GetType(),
                          (int)pDigi->GetSm(),(int)pDigi->GetRpc(),(int)pDigi->GetChannel(),(int)pDigi->GetSide())
          <<" : " << Form("Ind %2d, T %15.3f, Tot %5.1f",iDetIndx,pDigi->GetTime(),pDigi->GetTot());
    */
 
    if (fiBeamRefAddr < 0 || iAddr == fiBeamRefAddr) {
      //LOG(debug) << Form("Ref digi found for 0x%08x, Mask  0x%08x ", fiBeamRefAddr, DetMask);
      if (NULL == pRef) pRef = pDigi;
      else {
        if (pDigi->GetTime() < pRef->GetTime()) pRef = pDigi;
      }
    }
 
    if (fDigiBdfPar->GetNbDet() - 1 < iDetIndx || iDetIndx < 0) {
      LOG(warn) << Form(" Wrong DetIndx %d >< %d,0 ", iDetIndx, fDigiBdfPar->GetNbDet());
      break;
    }
 
    fhRpcDigiTot[iDetIndx]->Fill(2 * pDigi->GetChannel() + pDigi->GetSide(), pDigi->GetTot());
 
    if (NULL == fhRpcDigiCor[iDetIndx]) {
      if (100 < iMess++)
        LOG(warn) << Form(" DigiCor Histo for  DetIndx %d derived from 0x%08x not found", iDetIndx,
                          pDigi->GetAddress());
      continue;
    }
 
    Double_t dTDifMin     = dDoubleMax;
    CbmTofDigi* pDigi2Min = NULL;
    //       for (Int_t iDigI2 =iDigInd+1; iDigI2<iNbTofDigi;iDigI2++){
    for (Int_t iDigI2 = 0; iDigI2 < iNbTofDigi; iDigI2++) {
      CbmTofDigi* pDigi2 = &fvDigiIn[iDigI2];
      if (pDigi != pDigi2) { fhDigiTdif->Fill(pDigi->GetTime() - pDigi2->GetTime()); }
      if (pDigi->GetAddress() != pDigi2->GetAddress())
        if (pDigi->GetType() == 6 && pDigi2->GetType() == 6)
          fhBucDigiCor->Fill(pDigi->GetRpc() * 64 + pDigi->GetChannel() * 2 + pDigi->GetSide(),
                             pDigi2->GetRpc() * 64 + pDigi2->GetChannel() * 2 + pDigi2->GetSide());
 
      if (iDetIndx == fDigiBdfPar->GetDetInd(pDigi2->GetAddress())) {
        /*
        LOG(info) <<" Correlate DetIndx "<<iDetIndx
                      <<Form(", TSRCS %d%d%d%02d%d ", (int)pDigi->GetType(),
                                  (int)pDigi->GetSm(),(int)pDigi->GetRpc(),(int)pDigi->GetChannel(),(int)pDigi->GetSide())
                              <<Form(" with %d%d%d%02d%d ", (int)pDigi2->GetType(),
                                    (int)pDigi2->GetSm(),(int)pDigi2->GetRpc(),(int)pDigi2->GetChannel(),(int)pDigi2->GetSide());
        */
        if (0. == pDigi->GetSide() && 1. == pDigi2->GetSide()) {
          fhRpcDigiCor[iDetIndx]->Fill(pDigi->GetChannel(), pDigi2->GetChannel());
        }
        else {
          if (1. == pDigi->GetSide() && 0. == pDigi2->GetSide()) {
            fhRpcDigiCor[iDetIndx]->Fill(pDigi2->GetChannel(), pDigi->GetChannel());
          }
        }
        if (pDigi->GetSide() != pDigi2->GetSide()) {
          if (pDigi->GetChannel() == pDigi2->GetChannel()) {
            Double_t dTDif = TMath::Abs(pDigi->GetTime() - pDigi2->GetTime());
            if (dTDif < dTDifMin) {
              dTDifMin  = dTDif;
              pDigi2Min = pDigi2;
              //LOG(debug) << "Digi2 found at "<<iDigI2<<" with TDif = "<<dTDifMin<<" ns";
            }
          }
          else if (TMath::Abs(pDigi->GetChannel() - pDigi2->GetChannel())
                   == 1) {  // opposite side missing, neighbouring channel has hit on opposite side // FIXME
            // check that same side digi of neighbouring channel is absent
            Int_t iDigI3 = 0;
            for (; iDigI3 < iNbTofDigi; iDigI3++) {
              CbmTofDigi* pDigi3 = &fvDigiIn[iDigI3];
              if (pDigi3->GetSide() == pDigi->GetSide() && pDigi2->GetChannel() == pDigi3->GetChannel()) break;
            }
            if (iDigI3 == iNbTofDigi) {  // same side neighbour did not fire
              Int_t iCorMode = 0;        // Missing hit correction mode
              switch (iCorMode) {
                case 0:  // no action
                  break;
                case 1:  // shift found hit
                  LOG(debug) << Form("shift channel %d%d%d%d%d and  ", (Int_t) pDigi->GetType(), (Int_t) pDigi->GetSm(),
                                     (Int_t) pDigi->GetRpc(), (Int_t) pDigi->GetChannel(), (Int_t) pDigi->GetSide())
                             << Form(" %d%d%d%d%d ", (Int_t) pDigi2->GetType(), (Int_t) pDigi2->GetSm(),
                                     (Int_t) pDigi2->GetRpc(), (Int_t) pDigi2->GetChannel(), (Int_t) pDigi2->GetSide());
                  //if(pDigi->GetTime() < pDigi2->GetTime())
                  if (pDigi->GetSide() == 0)
                    pDigi2->SetAddress(pDigi->GetSm(), pDigi->GetRpc(), pDigi->GetChannel(), 1 - pDigi->GetSide(),
                                       pDigi->GetType());
                  else
                    pDigi->SetAddress(pDigi2->GetSm(), pDigi2->GetRpc(), pDigi2->GetChannel(), 1 - pDigi2->GetSide(),
                                      pDigi2->GetType());
 
                  LOG(debug) << Form("resultchannel %d%d%d%d%d and  ", (Int_t) pDigi->GetType(), (Int_t) pDigi->GetSm(),
                                     (Int_t) pDigi->GetRpc(), (Int_t) pDigi->GetChannel(), (Int_t) pDigi->GetSide())
                             << Form(" %d%d%d%d%d ", (Int_t) pDigi2->GetType(), (Int_t) pDigi2->GetSm(),
                                     (Int_t) pDigi2->GetRpc(), (Int_t) pDigi2->GetChannel(), (Int_t) pDigi2->GetSide());
                  break;
                case 2:  // insert missing hits
                  CbmTofDigi* pDigiN = new CbmTofDigi(*pDigi);
                  pDigiN->SetAddress(pDigi->GetSm(), pDigi->GetRpc(), pDigi2->GetChannel(), pDigi->GetSide(),
                                     pDigi->GetType());
                  pDigiN->SetTot(pDigi2->GetTot());
                  fvDigiIn.push_back(*pDigiN);
 
                  CbmTofDigi* pDigiN2 = new CbmTofDigi(*pDigi2);
                  pDigiN2->SetAddress(pDigi2->GetSm(), pDigi2->GetRpc(), pDigi->GetChannel(), pDigi2->GetSide(),
                                      pDigi2->GetType());
                  pDigiN2->SetTot(pDigi->GetTot());
                  fvDigiIn.push_back(*pDigiN2);
 
                  break;
              }
            }
          }
        }
      }
    }
    if (pDigi2Min != NULL) {
      CbmTofDetectorInfo xDetInfo(ECbmModuleId::kTof, pDigi->GetType(), pDigi->GetSm(), pDigi->GetRpc(), 0,
                                  pDigi->GetChannel());
      Int_t iChId  = fTofId->SetDetectorInfo(xDetInfo);
      fChannelInfo = fDigiPar->GetCell(iChId);
      if (NULL == fChannelInfo) {
        LOG(warn) << Form("Invalid ChannelInfo for 0x%08x, 0x%08x", iChId, pDigi2Min->GetAddress()) << " TSRC "
                  << pDigi->GetType() << pDigi->GetSm() << pDigi->GetRpc() << pDigi->GetChannel();
        continue;
      }
      if (fDigiBdfPar->GetSigVel(pDigi->GetType(), pDigi->GetSm(), pDigi->GetRpc()) * dTDifMin * 0.5
          < fPosYMaxScal * fChannelInfo->GetSizey()) {
        //check consistency
        if (8 == pDigi->GetType() || 5 == pDigi->GetType()) {
          if (pDigi->GetTime() != pDigi2Min->GetTime()) {
            if (fiMsgCnt-- > 0) {
              LOG(warn) << "Inconsistent duplicated digis in event " << fiNevtBuild << ", Ind: " << iDigInd;
              LOG(warn) << "   " << pDigi->ToString();
              LOG(warn) << "   " << pDigi2Min->ToString();
            }
            pDigi2Min->SetTot(pDigi->GetTot());
            pDigi2Min->SetTime(pDigi->GetTime());
          }
        }
      }
    }
  }
 
  if (NULL != pRef) {
    // LOG(debug) << Form("pRef from 0x%08x ",pRef->GetAddress());
 
    for (Int_t iDigInd = 0; iDigInd < fiNDigiIn; iDigInd++) {
      CbmTofDigi* pDigi = &fvDigiIn[iDigInd];
      Int_t iAddr       = pDigi->GetAddress() & DetMask;
      Int_t iDet        = fDetIdIndexMap[iAddr];  // Detector Index
      Int_t iSide       = pDigi->GetSide();
      fhDigiTimesRaw->Fill(iDet * 2 + iSide, pDigi->GetTime() - pRef->GetTime());
    }
  }
  return kTRUE;
}
 
/************************************************************************************/
Bool_t CbmDeviceHitBuilderTof::CalibRawDigis()
{
  CbmTofDigi* pDigi;
  CbmTofDigi* pCalDigi = NULL;
  Int_t iDigIndCal     = -1;
  // channel deadtime map
  std::map<Int_t, Double_t> mChannelDeadTime;
  fTofCalDigiVec->clear();
 
  Int_t iNbTofDigi = fvDigiIn.size();
  for (Int_t iDigInd = 0; iDigInd < iNbTofDigi; iDigInd++) {
    pDigi       = (CbmTofDigi*) &fvDigiIn[iDigInd];
    Int_t iAddr = pDigi->GetAddress();
    /*
    LOG(debug)<<"BC "  // Before Calibration
              <<Form("0x%08x",pDigi->GetAddress())<<" TSRC "
              <<pDigi->GetType()
              <<pDigi->GetSm()
              <<pDigi->GetRpc()
              <<Form("%2d",(Int_t)pDigi->GetChannel())<<" "
              <<pDigi->GetSide()<<" "
              <<Form("%f",pDigi->GetTime())<<" "
              <<pDigi->GetTot();
    */
    if (pDigi->GetType() == 5 || pDigi->GetType() == 8)  // for Pad counters generate fake digi to mockup a strip
      if (pDigi->GetSide() == 1) continue;               // skip one side to avoid double entries
 
    Bool_t bValid = kTRUE;
    std::map<Int_t, Double_t>::iterator it;
    it = mChannelDeadTime.find(iAddr);
    if (it != mChannelDeadTime.end()) {
      /*
      LOG(debug)<<"CCT found valid ChannelDeadtime entry "<<mChannelDeadTime[iAddr]
                <<", DeltaT "<<pDigi->GetTime()-mChannelDeadTime[iAddr];
      */
      if ((bValid = (pDigi->GetTime() > mChannelDeadTime[iAddr] + fdChannelDeadtime))) {
        //  pCalDigi =
        //      new ((*fTofCalDigisColl)[++iDigIndCal]) CbmTofDigi(*pDigi);
        fTofCalDigiVec->push_back(CbmTofDigi(*pDigi));
        pCalDigi = &(fTofCalDigiVec->back());
        iDigIndCal++;
      }
    }
    else {
      // pCalDigi = new ((*fTofCalDigisColl)[++iDigIndCal]) CbmTofDigi(*pDigi);
      fTofCalDigiVec->push_back(CbmTofDigi(*pDigi));
      pCalDigi = &(fTofCalDigiVec->back());
      iDigIndCal++;
    }
    mChannelDeadTime[iAddr] = pDigi->GetTime();
    if (!bValid) continue;
    if (pRef != NULL)
      if (pDigi == pRef) pRefCal = pCalDigi;
    /*
    LOG(debug)<<"DC "  // After deadtime check. before Calibration
              <<Form("0x%08x",pDigi->GetAddress())<<" TSRC "
              <<pDigi->GetType()
              <<pDigi->GetSm()
              <<pDigi->GetRpc()
              <<Form("%2d",(Int_t)pDigi->GetChannel())<<" "
              <<pDigi->GetSide()<<" "
              <<Form("%f",pDigi->GetTime())<<" "
              <<pDigi->GetTot();
    */
 
    if (fbPs2Ns) {
      pCalDigi->SetTime(pCalDigi->GetTime() / 1000.);  // for backward compatibility
      pCalDigi->SetTot(pCalDigi->GetTot() / 1000.);    // for backward compatibility
    }
    if (fDigiBdfPar->GetNbSmTypes() > pDigi->GetType()  // prevent crash due to misconfiguration
        && fDigiBdfPar->GetNbSm(pDigi->GetType()) > pDigi->GetSm()
        && fDigiBdfPar->GetNbRpc(pDigi->GetType()) > pDigi->GetRpc()
        && fDigiBdfPar->GetNbChan(pDigi->GetType(), pDigi->GetRpc()) > pDigi->GetChannel()) {
      // apply calibration vectors
      pCalDigi->SetTime(pCalDigi->GetTime() -  // calibrate Digi Time
                        fvCPTOff[pDigi->GetType()][pDigi->GetSm() * fDigiBdfPar->GetNbRpc(pDigi->GetType())
                                                   + pDigi->GetRpc()][pDigi->GetChannel()][pDigi->GetSide()]);
      //      LOG(debug)<<" CluCal-TOff: "<<pCalDigi->ToString();
 
      Double_t dTot = pCalDigi->GetTot() -  // subtract Offset
                      fvCPTotOff[pDigi->GetType()][pDigi->GetSm() * fDigiBdfPar->GetNbRpc(pDigi->GetType())
                                                   + pDigi->GetRpc()][pDigi->GetChannel()][pDigi->GetSide()];
      if (dTot < 0.001) dTot = 0.001;
      pCalDigi->SetTot(dTot *  // calibrate Digi ToT
                       fvCPTotGain[pDigi->GetType()][pDigi->GetSm() * fDigiBdfPar->GetNbRpc(pDigi->GetType())
                                                     + pDigi->GetRpc()][pDigi->GetChannel()][pDigi->GetSide()]);
 
      // walk correction
      Double_t dTotBinSize = (fdTOTMax - fdTOTMin) / nbClWalkBinX;
      Int_t iWx            = (Int_t)((pCalDigi->GetTot() - fdTOTMin) / dTotBinSize);
      if (0 > iWx) iWx = 0;
      if (iWx >= nbClWalkBinX) iWx = nbClWalkBinX - 1;
      Double_t dDTot = (pCalDigi->GetTot() - fdTOTMin) / dTotBinSize - (Double_t) iWx - 0.5;
      Double_t dWT =
        fvCPWalk[pCalDigi->GetType()][pCalDigi->GetSm() * fDigiBdfPar->GetNbRpc(pCalDigi->GetType())
                                      + pCalDigi->GetRpc()][pCalDigi->GetChannel()][pCalDigi->GetSide()][iWx];
      if (dDTot > 0) {                 // linear interpolation to next bin
        if (iWx < nbClWalkBinX - 1) {  // linear interpolation to next bin
 
          dWT +=
            dDTot
            * (fvCPWalk[pCalDigi->GetType()][pCalDigi->GetSm() * fDigiBdfPar->GetNbRpc(pCalDigi->GetType())
                                             + pCalDigi->GetRpc()][pCalDigi->GetChannel()][pCalDigi->GetSide()][iWx + 1]
               - fvCPWalk[pCalDigi->GetType()]
                         [pCalDigi->GetSm() * fDigiBdfPar->GetNbRpc(pCalDigi->GetType()) + pCalDigi->GetRpc()]
                         [pCalDigi->GetChannel()][pCalDigi->GetSide()][iWx]);  //memory leak???
        }
      }
      else  // dDTot < 0,  linear interpolation to next bin
      {
        if (0 < iWx) {  // linear interpolation to next bin
          dWT -=
            dDTot
            * (fvCPWalk[pCalDigi->GetType()][pCalDigi->GetSm() * fDigiBdfPar->GetNbRpc(pCalDigi->GetType())
                                             + pCalDigi->GetRpc()][pCalDigi->GetChannel()][pCalDigi->GetSide()][iWx - 1]
               - fvCPWalk[pCalDigi->GetType()]
                         [pCalDigi->GetSm() * fDigiBdfPar->GetNbRpc(pCalDigi->GetType()) + pCalDigi->GetRpc()]
                         [pCalDigi->GetChannel()][pCalDigi->GetSide()][iWx]);  //memory leak???
        }
      }
      pCalDigi->SetTime(pCalDigi->GetTime() - dWT);  // calibrate Digi Time
      // LOG(debug)<<" CluCal-Walk: "<<pCalDigi->ToString();
    }
    else {
      LOG(info) << "Skip1 Digi "
                << " Type " << pDigi->GetType() << " " << fDigiBdfPar->GetNbSmTypes() << " Sm " << pDigi->GetSm() << " "
                << fDigiBdfPar->GetNbSm(pDigi->GetType()) << " Rpc " << pDigi->GetRpc() << " "
                << fDigiBdfPar->GetNbRpc(pDigi->GetType()) << " Ch " << pDigi->GetChannel() << " "
                << fDigiBdfPar->GetNbChan(pDigi->GetType(), 0);
    }
    if (pCalDigi->GetType() == 5
        || pCalDigi->GetType() == 8) {  // for Pad counters generate fake digi to mockup a strip
      //CbmTofDigi* pCalDigi2 =
      //  new ((*fTofCalDigisColl)[++iDigIndCal]) CbmTofDigi(*pCalDigi);
      fTofCalDigiVec->push_back(CbmTofDigi(*pCalDigi));
      CbmTofDigi* pCalDigi2 = &(fTofCalDigiVec->back());
      iDigIndCal++;
      if (pCalDigi->GetSide() == 0)
        pCalDigi2->SetAddress(pCalDigi->GetSm(), pCalDigi->GetRpc(), pCalDigi->GetChannel(), 1, pCalDigi->GetType());
      else
        pCalDigi2->SetAddress(pCalDigi->GetSm(), pCalDigi->GetRpc(), pCalDigi->GetChannel(), 0, pCalDigi->GetType());
    }
  }  // for( Int_t iDigInd = 0; iDigInd < nTofDigi; iDigInd++ )
 
  iNbTofDigi = fTofCalDigiVec->size();
  //  fTofCalDigisColl->GetEntriesFast();  // update because of added duplicted digis
  LOG(debug) << "CbmTofHitBuilder: Sort " << fTofCalDigiVec->size() << " calibrated digis ";
  if (iNbTofDigi > 1) {
    //    fTofCalDigisColl->Sort(iNbTofDigi); // Time order again, in case modified by the calibration
    std::sort(fTofCalDigiVec->begin(), fTofCalDigiVec->end(),
              [](const CbmTofDigi& a, const CbmTofDigi& b) -> bool { return a.GetTime() < b.GetTime(); });
    //    std::sort(fTofCalDigiVec->begin(), fTofCalDigiVec->end());
    //    if(!fTofCalDigisColl->IsSorted()){
    //    if ( ! std::is_sorted(fTofCalDigiVec->begin(), fTofCalDigiVec->end()))
    if (!std::is_sorted(fTofCalDigiVec->begin(), fTofCalDigiVec->end(),
                        [](const CbmTofDigi& a, const CbmTofDigi& b) -> bool { return a.GetTime() < b.GetTime(); }))
      LOG(warning) << "CbmTofHitBuilder: Digi sorting not successful ";
  }
 
  if (NULL != pRef) {
    // LOG(debug) << Form("pRef from 0x%08x ",pRef->GetAddress());
 
    for (Int_t iDigInd = 0; iDigInd < iNbTofDigi; iDigInd++) {
      // pDigi       = (CbmTofDigi*) fTofCalDigisColl->At(iDigInd);
      pDigi       = &(fTofCalDigiVec->at(iDigInd));
      Int_t iAddr = pDigi->GetAddress() & DetMask;
      Int_t iDet  = fDetIdIndexMap[iAddr];  // Detector Index
      Int_t iSide = pDigi->GetSide();
      fhDigiTimesCor->Fill(iDet * 2 + iSide, pDigi->GetTime() - pRefCal->GetTime());
    }
  }
  return kTRUE;
}
 
/************************************************************************************/
Bool_t CbmDeviceHitBuilderTof::BuildHits()
{
  Int_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
  // Hit variables
  Double_t dWeightedTime = 0.0;
  Double_t dWeightedPosX = 0.0;
  Double_t dWeightedPosY = 0.0;
  Double_t dWeightedPosZ = 0.0;
  Double_t dWeightsSum   = 0.0;
  //vPtsRef.clear();
  vDigiIndRef.clear();
  //CbmTofCell *fTrafoCell=NULL;
  //Int_t iTrafoCell=-1;
  Int_t iNbChanInHit = 0;
  // Last Channel Temp variables
  Int_t iLastChan    = -1;
  Double_t dLastPosX = 0.0;  // -> Comment to remove warning because set but never used
  Double_t dLastPosY = 0.0;
  Double_t dLastTime = 0.0;
  // Channel Temp variables
  Double_t dPosX     = 0.0;
  Double_t dPosY     = 0.0;
  Double_t dPosZ     = 0.0;
  Double_t dTime     = 0.0;
  Double_t dTimeDif  = 0.0;
  Double_t dTotS     = 0.0;
  Int_t fiNbSameSide = 0;
 
  //  gGeoManager->SetTopVolume( gGeoManager->FindVolumeFast("tof_v14a") );
  gGeoManager->SetTopVolume(gGeoManager->FindVolumeFast("cave"));
  gGeoManager->CdTop();
 
  if (kTRUE == fDigiBdfPar->UseExpandedDigi()) {
    for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
      Int_t iNbSm  = fDigiBdfPar->GetNbSm(iSmType);
      Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
      for (Int_t iSm = 0; iSm < iNbSm; iSm++)
        for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
          Int_t iNbCh   = fDigiBdfPar->GetNbChan(iSmType, iRpc);
          Int_t iChType = fDigiBdfPar->GetChanType(iSmType, iRpc);
          /*
          LOG(debug)<<"RPC - Loop  "
                    << Form(" %3d %3d %3d %3d ",iSmType,iSm,iRpc,iChType);
          */
          fviClusterMul[iSmType][iSm][iRpc] = 0;
          Int_t iChId                       = 0;
          Int_t iDetId                      = CbmTofAddress::GetUniqueAddress(iSm, iRpc, 0, 0, iSmType);
          ;
          Int_t iDetIndx = fDetIdIndexMap[iDetId];  // Detector Index
          if (0 == iChType) {
            // Don't spread clusters over RPCs!!!
            dWeightedTime = 0.0;
            dWeightedPosX = 0.0;
            dWeightedPosY = 0.0;
            dWeightedPosZ = 0.0;
            dWeightsSum   = 0.0;
            iNbChanInHit  = 0;
            //vPtsRef.clear();
            // For safety reinitialize everything
            iLastChan = -1;
            //                  dLastPosX = 0.0; // -> Comment to remove warning because set but never used
            dLastPosY = 0.0;
            dLastTime = 0.0;
            //LOG(debug2)<<"ChanOrient "
            //          << Form(" %3d %3d %3d %3d %3d ",iSmType,iSm,iRpc,fDigiBdfPar->GetChanOrient( iSmType, iRpc ),iNbCh);
 
            if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc)) {
              // Horizontal strips => X comes from left right time difference
            }       // if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
            else {  // Vertical strips => Y comes from bottom top time difference
              for (Int_t iCh = 0; iCh < iNbCh; iCh++) {
                //LOG(debug3)<<"VDigisize "
                //          << Form(" T %3d Sm %3d R %3d Ch %3d Size %3lu ",
                //                  iSmType,iSm,iRpc,iCh,fStorStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh].size());
                if (0 == fStorDigi[iSmType][iSm * iNbRpc + iRpc].size()) continue;
                if (fvDeadStrips[iDetIndx] & (1 << iCh)) continue;  // skip over dead channels
                //if( 0 < fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh].size() )
                //  fhNbDigiPerChan->Fill( fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh].size() );
 
                while (1 < fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size()) {
 
                  while ((fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][0])->GetSide()
                         == (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][1])->GetSide()) {
                    // Not one Digi of each end!
                    fiNbSameSide++;
                    if (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size() > 2) {
                      LOG(debug) << "SameSide Digis! on TSRC " << iSmType << iSm << iRpc << iCh << ", Times: "
                                 << Form("%f", (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][0])->GetTime()) << ", "
                                 << Form("%f", (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][1])->GetTime())
                                 << ", DeltaT "
                                 << (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][1])->GetTime()
                                      - (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][0])->GetTime()
                                 << ", array size: " << fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size();
                      if (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][2]->GetSide()
                          == fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][0]->GetSide()) {
                        LOG(debug) << "3 consecutive SameSide Digis! on TSRC " << iSmType << iSm << iRpc << iCh
                                   << ", Times: " << (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][0])->GetTime()
                                   << ", " << (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][1])->GetTime()
                                   << ", DeltaT "
                                   << (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][1])->GetTime()
                                        - (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][0])->GetTime()
                                   << ", array size: " << fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size();
                        fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                          fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
                        fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                          fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
                      }
                      else {
                        if (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][2]->GetTime()
                              - fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][0]->GetTime()
                            > fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][2]->GetTime()
                                - fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][1]->GetTime()) {
                          fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                            fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
                          fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                            fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
                        }
                        else {
                          LOG(debug) << Form("Ev %8.0f, digis not properly time ordered, TSRCS "
                                             "%d%d%d%d%d ",
                                             fdEvent, iSmType, iSm, iRpc, iCh,
                                             (Int_t) fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][0]->GetSide());
                          fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                            fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin() + 1);
                          fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                            fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin() + 1);
                        }
                      }
                    }
                    else {
                      LOG(debug) << "SameSide Erase fStor entries(d) " << iSmType << ", SR " << iSm * iNbRpc + iRpc
                                 << ", Ch" << iCh;
                      fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                        fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
                      fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                        fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
                    }
                    if (2 > fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size()) break;
                    continue;
                  }  // same condition side end
                  LOG(debug) << "digis processing for "
                             << Form(" SmT %3d Sm %3d Rpc %3d Ch %3d # %3lu ", iSmType, iSm, iRpc, iCh,
                                     fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size());
                  if (2 > fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size()) {
                    LOG(debug) << Form("Leaving digi processing for TSRC %d%d%d%d, size  %3lu", iSmType, iSm, iRpc, iCh,
                                       fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size());
                    break;
                  }
                  /* Int_t iLastChId = iChId; // Save Last hit channel*/
 
                  // 2 Digis = both sides present
                  CbmTofDetectorInfo xDetInfo(ECbmModuleId::kTof, iSmType, iSm, iRpc, 0, iCh);
                  iChId          = fTofId->SetDetectorInfo(xDetInfo);
                  Int_t iUCellId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, iCh, 0, iSmType);
                  LOG(debug) << Form("TSRC %d%d%d%d size %3lu ", iSmType, iSm, iRpc, iCh,
                                     fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size())
                             << Form(" ChId: 0x%08x 0x%08x ", iChId, iUCellId);
                  fChannelInfo = fDigiPar->GetCell(iChId);
 
                  if (NULL == fChannelInfo) {
                    LOG(error) << "BuildHits: no geometry info! "
                               << Form(" %3d %3d %3d %3d 0x%08x 0x%08x ", iSmType, iSm, iRpc, iCh, iChId, iUCellId);
                    break;
                  }
 
                  TGeoNode* fNode =  // prepare local->global trafo
                    gGeoManager->FindNode(fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ());
                  //          fNode->Print();
                  if (NULL == fNode) {  // Transformation matrix not available !!!??? - Check
                    LOG(error) << Form("Node at (%6.1f,%6.1f,%6.1f) : %p", fChannelInfo->GetX(), fChannelInfo->GetY(),
                                       fChannelInfo->GetZ(), fNode);
                    ChangeState(fair::mq::Transition::Stop);
                  }
 
                  CbmTofDigi* xDigiA = fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][0];
                  CbmTofDigi* xDigiB = fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][1];
 
                  dTimeDif = (xDigiA->GetTime() - xDigiB->GetTime());
                  if (5 == iSmType && dTimeDif != 0.) {
                    // FIXME -> Overflow treatment in calib/tdc/TMbsCalibTdcTof.cxx
                    LOG(debug) << "BuildHits: Diamond hit in " << iSm << " with inconsistent digits "
                               << xDigiA->GetTime() << ", " << xDigiB->GetTime() << " -> " << dTimeDif;
                    /*
                    LOG(debug) << "    "<<xDigiA->ToString();
                    LOG(debug) << "    "<<xDigiB->ToString();
                    */
                  }
                  if (1 == xDigiA->GetSide())
                    // 0 is the top side, 1 is the bottom side
                    dPosY = fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc) * dTimeDif * 0.5;
                  else
                    // 0 is the bottom side, 1 is the top side
                    dPosY = -fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc) * dTimeDif * 0.5;
 
                  if (TMath::Abs(dPosY) > fChannelInfo->GetSizey()
                      && fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size() > 2) {
                    LOG(debug) << "Hit candidate outside correlation window, check for "
                                  "better possible digis, "
                               << " mul " << fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size();
 
                    CbmTofDigi* xDigiC = fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][2];
                    Double_t dPosYN    = 0.;
                    Double_t dTimeDifN = 0;
                    if (xDigiC->GetSide() == xDigiA->GetSide()) dTimeDifN = xDigiC->GetTime() - xDigiB->GetTime();
                    else
                      dTimeDifN = xDigiA->GetTime() - xDigiC->GetTime();
 
                    if (1 == xDigiA->GetSide()) dPosYN = fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc) * dTimeDifN * 0.5;
                    else
                      dPosYN = -fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc) * dTimeDifN * 0.5;
 
                    if (TMath::Abs(dPosYN) < TMath::Abs(dPosY)) {
                      LOG(debug) << "Replace digi on side " << xDigiC->GetSide() << ", yPosNext " << dPosYN
                                 << " old: " << dPosY;
                      dTimeDif = dTimeDifN;
                      dPosY    = dPosYN;
                      if (xDigiC->GetSide() == xDigiA->GetSide()) {
                        xDigiA = xDigiC;
                        fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                          fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
                        fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                          fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
                      }
                      else {
                        xDigiB = xDigiC;
                        fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                          ++(fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin() + 1));
                        fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                          ++(fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin() + 1));
                      }
                    }
                  }
                  if (xDigiA->GetSide() == xDigiB->GetSide()) {
                    LOG(error) << "Wrong combinations of digis " << fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][0]
                               << "," << fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][1];
                  }
                  // The "Strip" time is the mean time between each end
                  dTime = 0.5 * (xDigiA->GetTime() + xDigiB->GetTime());
 
                  // Weight is the total charge => sum of both ends ToT
                  dTotS = xDigiA->GetTot() + xDigiB->GetTot();
 
                  // use local coordinates, (0,0,0) is in the center of counter  ?
                  dPosX = ((Double_t)(-iNbCh / 2 + iCh) + 0.5) * fChannelInfo->GetSizex();
                  dPosZ = 0.;
 
                  /*
                  LOG(debug)
                    <<"NbChanInHit  "
                    << Form(" %3d %3d %3d %3d %3d 0x%p %1.0f Time %f PosX %f PosY %f Svel %f ",
                            iNbChanInHit,iSmType,iRpc,iCh,iLastChan,xDigiA,xDigiA->GetSide(),
                            dTime,dPosX,dPosY,fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc))
                    //  << Form( ", Offs %f, %f ",fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0],
                    //                            fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1])
                    ;
                  */
                  // Now check if a hit/cluster is already started
                  if (0 < iNbChanInHit) {
                    if (iLastChan == iCh - 1) {
                      /*
                       fhDigTimeDifClust->Fill( dTime - dLastTime );
                       fhDigSpacDifClust->Fill( dPosY - dLastPosY );
                       fhDigDistClust->Fill( dPosY - dLastPosY,
                                             dTime - dLastTime );
                      */
                    }
                    // if( iLastChan == iCh - 1 )
                    // a cluster is already started => check distance in space/time
                    // For simplicity, just check along strip direction for now
                    // and break cluster when a not fired strip is found
                    if (TMath::Abs(dTime - dLastTime) < fdMaxTimeDist && iLastChan == iCh - 1
                        && TMath::Abs(dPosY - dLastPosY) < fdMaxSpaceDist) {
                      // Add to cluster/hit
                      dWeightedTime += dTime * dTotS;
                      dWeightedPosX += dPosX * dTotS;
                      dWeightedPosY += dPosY * dTotS;
                      dWeightedPosZ += dPosZ * dTotS;
                      dWeightsSum += dTotS;
                      iNbChanInHit += 1;
 
                      vDigiIndRef.push_back((Int_t)(fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][0]));
                      vDigiIndRef.push_back((Int_t)(fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][1]));
 
                      LOG(debug) << " Add Digi and erase fStor entries(a): NbChanInHit " << iNbChanInHit << ", "
                                 << iSmType << ", SR " << iSm * iNbRpc + iRpc << ", Ch" << iCh;
 
                      fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                        fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
                      fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                        fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
                      fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                        fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
                      fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                        fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
 
                    }  // if current Digis compatible with last fired chan
                    else {
                      // Save Hit
                      dWeightedTime /= dWeightsSum;
                      dWeightedPosX /= dWeightsSum;
                      dWeightedPosY /= dWeightsSum;
                      dWeightedPosZ /= dWeightsSum;
                      //  TVector3 hitPosLocal(dWeightedPosX, dWeightedPosY, dWeightedPosZ);
                      //TVector3 hitPos;
                      Double_t hitpos_local[3];
                      hitpos_local[0] = dWeightedPosX;
                      hitpos_local[1] = dWeightedPosY;
                      hitpos_local[2] = dWeightedPosZ;
 
                      Double_t hitpos[3];
                      //TGeoNode*         cNode   =
                      gGeoManager->GetCurrentNode();
                      /*TGeoHMatrix* cMatrix =*/gGeoManager->GetCurrentMatrix();
                      //cNode->Print();
                      //cMatrix->Print();
                      gGeoManager->LocalToMaster(hitpos_local, hitpos);
                      /*
                      LOG(debug)<<
                        Form("LocalToMaster for node %p: (%6.1f,%6.1f,%6.1f) ->(%6.1f,%6.1f,%6.1f)",
                             cNode, hitpos_local[0], hitpos_local[1], hitpos_local[2],
                             hitpos[0], hitpos[1], hitpos[2]);
                      */
                      TVector3 hitPos(hitpos[0], hitpos[1], hitpos[2]);
 
                      // Simple errors, not properly done at all for now
                      // Right way of doing it should take into account the weight distribution
                      // and real system time resolution
                      TVector3 hitPosErr(0.5, 0.5, 0.5);  // including positioning uncertainty
                      /*
                        TVector3 hitPosErr( fChannelInfo->GetSizex()/TMath::Sqrt(12.0),   // Single strips approximation
                        0.5, // Use generic value
                        1.);
 
                      */
                      //fDigiBdfPar->GetFeeTimeRes() * fDigiBdfPar->GetSigVel(iSmType,iRpc), // Use the electronics resolution
                      //fDigiBdfPar->GetNbGaps( iSmType, iRpc)*
                      //fDigiBdfPar->GetGapSize( iSmType, iRpc)/ //10.0 / // Change gap size in cm
                      //TMath::Sqrt(12.0) ); // Use full RPC thickness as "Channel" Z size
 
                      // Int_t iDetId = vPtsRef[0]->GetDetectorID();// detID = pt->GetDetectorID() <= from TofPoint
                      // calc mean ch from dPosX=((Double_t)(-iNbCh/2 + iCh)+0.5)*fChannelInfo->GetSizex();
 
                      Int_t iChm = floor(dWeightedPosX / fChannelInfo->GetSizex()) + iNbCh / 2;
                      if (iChm < 0) iChm = 0;
                      if (iChm > iNbCh - 1) iChm = iNbCh - 1;
                      iDetId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, iChm, 0, iSmType);
                      //Int_t iRefId = 0; // Index of the correspondng TofPoint
                      //                       if(NULL != fTofPointsColl) {
                      //iRefId = fTofPointsColl->IndexOf( vPtsRef[0] );
                      //}
                      /*
                      LOG(debug)<<"Save Hit  "
                                << Form(" %3d %3d 0x%08x %3d %3d %3d %f %f",
                                        fiNbHits,iNbChanInHit,iDetId,iChm,iLastChan,iRefId,
                                        dWeightedTime,dWeightedPosY)
                                <<", DigiSize: "<<vDigiIndRef.size()
                                <<", DigiInds: ";
                      for (UInt_t i=0; i<vDigiIndRef.size();i++){
                        LOG(debug)<<" "<<vDigiIndRef.at(i)<<"(M"<<fviClusterMul[iSmType][iSm][iRpc]<<")";
                      }
                      */
 
                      fviClusterMul[iSmType][iSm][iRpc]++;
                      if (vDigiIndRef.size() < 2) {
                        LOG(warn) << "Digi refs for Hit " << fiNbHits << ": " << vDigiIndRef.size();
                      }
                      if (fiNbHits > 0) {
                        CbmTofHit* pHitL = (CbmTofHit*) fTofHitsColl->At(fiNbHits - 1);
                        if (iDetId == pHitL->GetAddress() && dWeightedTime == pHitL->GetTime()) {
                          LOG(debug) << "Store Hit twice? "
                                     << " fiNbHits " << fiNbHits << ", " << Form("0x%08x", iDetId);
                        }
                      }
                      CbmTofHit* pHit =
                        new CbmTofHit(iDetId, hitPos,
                                      hitPosErr,  //local detector coordinates
                                      fiNbHits,   // this number is used as reference!!
                                      dWeightedTime,
                                      vDigiIndRef.size(),  // number of linked digis =  2*CluSize
                                      //vPtsRef.size(), // flag  = number of TofPoints generating the cluster
                                      Int_t(dWeightsSum * 10.));  //channel -> Tot
                      //0) ; //channel
                      // output hit
                      new ((*fTofHitsColl)[fiNbHits]) CbmTofHit(*pHit);
                      // memorize hit
                      if (fdMemoryTime > 0.) { LH_store(iSmType, iSm, iRpc, iChm, pHit); }
                      else {
                        pHit->Delete();
                      }
                      /*
                        new((*fTofDigiMatchColl)[fiNbHits]) CbmMatch();
                        CbmMatch* digiMatch = (CbmMatch *)fTofDigiMatchColl->At(fiNbHits);
                      */
                      CbmMatch* digiMatch = new ((*fTofDigiMatchColl)[fiNbHits]) CbmMatch();
                      for (UInt_t i = 0; i < vDigiIndRef.size(); i++) {
                        Double_t dTot = ((CbmTofDigi*) (&(fTofCalDigiVec->at(vDigiIndRef.at(i)))))->GetTot();
                        digiMatch->AddLink(CbmLink(dTot, vDigiIndRef.at(i), fiOutputTreeEntry, fiFileIndex));
                      }
 
                      fiNbHits++;
                      // For Histogramming
                      fviClusterSize[iSmType][iRpc].push_back(iNbChanInHit);
                      //fviTrkMul[iSmType][iRpc].push_back( vPtsRef.size() );
                      fvdX[iSmType][iRpc].push_back(dWeightedPosX);
                      fvdY[iSmType][iRpc].push_back(dWeightedPosY);
                      /*  no TofPoint available for data!
                          fvdDifX[iSmType][iRpc].push_back( vPtsRef[0]->GetX() - dWeightedPosX);
                          fvdDifY[iSmType][iRpc].push_back( vPtsRef[0]->GetY() - dWeightedPosY);
                          fvdDifCh[iSmType][iRpc].push_back( fGeoHandler->GetCell( vPtsRef[0]->GetDetectorID() ) -1 -iLastChan );
                      */
                      //vPtsRef.clear();
                      vDigiIndRef.clear();
 
                      // Start a new hit
                      dWeightedTime = dTime * dTotS;
                      dWeightedPosX = dPosX * dTotS;
                      dWeightedPosY = dPosY * dTotS;
                      dWeightedPosZ = dPosZ * dTotS;
                      dWeightsSum   = dTotS;
                      iNbChanInHit  = 1;
                      // Save pointer on CbmTofPoint
                      // vPtsRef.push_back( (CbmTofPoint*)(xDigiA->GetLinks()) );
                      // Save next digi address
 
                      vDigiIndRef.push_back((Int_t)(fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][0]));
                      vDigiIndRef.push_back((Int_t)(fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][1]));
                      //LOG(debug2)<<"Erase fStor entries(b) "<<iSmType<<", SR "<<iSm*iNbRpc+iRpc<<", Ch"<<iCh;
 
                      fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                        fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
                      fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                        fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
                      fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                        fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
                      fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                        fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
                    }  // else of if current Digis compatible with last fired chan
                  }    // if( 0 < iNbChanInHit)
                  else {
                    LOG(debug) << Form("1.Hit on channel %d, time: %f, PosY %f", iCh, dTime, dPosY);
 
                    // first fired strip in this RPC
                    dWeightedTime = dTime * dTotS;
                    dWeightedPosX = dPosX * dTotS;
                    dWeightedPosY = dPosY * dTotS;
                    dWeightedPosZ = dPosZ * dTotS;
                    dWeightsSum   = dTotS;
                    iNbChanInHit  = 1;
                    // Save pointer on CbmTofPoint
                    //if(NULL != fTofPointsColl)
                    //                                    vPtsRef.push_back( (CbmTofPoint*)(xDigiA->GetLinks()) );
                    vDigiIndRef.push_back((Int_t)(fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][0]));
                    vDigiIndRef.push_back((Int_t)(fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][1]));
 
                    //LOG(debug)<<"Erase fStor entries(c) "<<iSmType<<", SR "<<iSm*iNbRpc+iRpc<<", Ch"<<iCh;
 
                    fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                      fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
                    fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                      fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
                    fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                      fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
                    fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                      fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
 
                  }  // else of if( 0 < iNbChanInHit)
                  iLastChan = iCh;
                  dLastPosX = dPosX;
                  dLastPosY = dPosY;
                  dLastTime = dTime;
                  if (AddNextChan(iSmType, iSm, iRpc, iLastChan, dLastPosX, dLastPosY, dLastTime, dWeightsSum)) {
                    iNbChanInHit = 0;  // cluster already stored
                  }
                }  // while( 1 < fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh].size() )
                fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].clear();
                fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].clear();
              }  // for( Int_t iCh = 0; iCh < iNbCh; iCh++ )
              //LOG(debug2)<<"finished V-RPC"
              //                  << Form(" %3d %3d %3d %d %f %fx",iSmType,iSm,iRpc,fTofHitsColl->GetEntriesFast(),dLastPosX,dLastPosY);
            }  // else of if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
          }    // if( 0 == iChType)
          else {
            LOG(error) << "=> Cluster building "
                       << "from digis to hits not implemented for pads, Sm type " << iSmType << " Rpc " << iRpc;
            return kFALSE;
          }  // else of if( 0 == iChType)
 
          // Now check if another hit/cluster is started
          // and save it if it's the case
          if (0 < iNbChanInHit) {
            /*
            LOG(debug)<<"Process cluster "
                      <<iNbChanInHit;
            */
            // Check orientation to properly assign errors
            if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc)) {
              // LOG(debug1)<<"H-Hit ";
            }  // if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
            else {
              //LOG(debug2)<<"V-Hit ";
              // Save Hit
              dWeightedTime /= dWeightsSum;
              dWeightedPosX /= dWeightsSum;
              dWeightedPosY /= dWeightsSum;
              dWeightedPosZ /= dWeightsSum;
              //TVector3 hitPos(dWeightedPosX, dWeightedPosY, dWeightedPosZ);
 
              Double_t hitpos_local[3] = {3 * 0.};
              hitpos_local[0]          = dWeightedPosX;
              hitpos_local[1]          = dWeightedPosY;
              hitpos_local[2]          = dWeightedPosZ;
 
              Double_t hitpos[3];
              //TGeoNode*        cNode=
              gGeoManager->GetCurrentNode();
              /*TGeoHMatrix* cMatrix =*/gGeoManager->GetCurrentMatrix();
              //cNode->Print();
              //cMatrix->Print();
 
              gGeoManager->LocalToMaster(hitpos_local, hitpos);
              /*
              LOG(debug)<<
                Form(" LocalToMaster for V-node %p: (%6.1f,%6.1f,%6.1f) ->(%6.1f,%6.1f,%6.1f)",
                     cNode, hitpos_local[0], hitpos_local[1], hitpos_local[2],
                     hitpos[0], hitpos[1], hitpos[2])
                ;
              */
              TVector3 hitPos(hitpos[0], hitpos[1], hitpos[2]);
              // Event errors, not properly done at all for now
              // Right way of doing it should take into account the weight distribution
              // and real system time resolution
              TVector3 hitPosErr(0.5, 0.5, 0.5);  // including positioning uncertainty
              /*
                TVector3 hitPosErr( fChannelInfo->GetSizex()/TMath::Sqrt(12.0),   // Single strips approximation
                0.5, // Use generic value
                1.);
              */
              Int_t iChm = floor(dWeightedPosX / fChannelInfo->GetSizex()) + iNbCh / 2;
              if (iChm < 0) iChm = 0;
              if (iChm > iNbCh - 1) iChm = iNbCh - 1;
              iDetId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, iChm, 0, iSmType);
              //Int_t iRefId = 0; // Index of the correspondng TofPoint
              //if(NULL != fTofPointsColl) iRefId = fTofPointsColl->IndexOf( vPtsRef[0] );
              /*
              LOG(debug)<<"Save V-Hit  "
                        << Form(" %3d %3d 0x%08x %3d 0x%08x", // %3d %3d
                                fiNbHits,iNbChanInHit,iDetId,iLastChan,iRefId) //vPtsRef.size(),vPtsRef[0])
                //   dWeightedTime,dWeightedPosY)
                        <<", DigiSize: "<<vDigiIndRef.size();
              for (UInt_t i=0; i<vDigiIndRef.size();i++){
                LOG(debug)<<"DigiIndRef "<<i<<" "<<vDigiIndRef.at(i)<<"(M"<<fviClusterMul[iSmType][iSm][iRpc]<<")";
              }
              */
              fviClusterMul[iSmType][iSm][iRpc]++;
              if (vDigiIndRef.size() < 2) {
                LOG(warn) << "Digi refs for Hit " << fiNbHits << ":  " << vDigiIndRef.size();
              }
              if (fiNbHits > 0) {
                CbmTofHit* pHitL = (CbmTofHit*) fTofHitsColl->At(fiNbHits - 1);
                if (iDetId == pHitL->GetAddress() && dWeightedTime == pHitL->GetTime())
                  LOG(debug) << "Store Hit twice? "
                             << " fiNbHits " << fiNbHits << ", " << Form("0x%08x", iDetId);
              }
 
              CbmTofHit* pHit = new CbmTofHit(iDetId, hitPos,
                                              hitPosErr,  //local detector coordinates
                                              fiNbHits,   // this number is used as reference!!
                                              dWeightedTime,
                                              vDigiIndRef.size(),  // number of linked digis =  2*CluSize
                                              //vPtsRef.size(), // flag  = number of TofPoints generating the cluster
                                              Int_t(dWeightsSum * 10.));  //channel -> Tot
              //                0) ; //channel
              //                vDigiIndRef);
              // output hit
              new ((*fTofHitsColl)[fiNbHits]) CbmTofHit(*pHit);
              // memorize hit
              if (fdMemoryTime > 0.) { LH_store(iSmType, iSm, iRpc, iChm, pHit); }
              else {
                pHit->Delete();
              }
              /*
                new((*fTofDigiMatchColl)[fiNbHits]) CbmMatch();
                CbmMatch* digiMatch = (CbmMatch *)fTofDigiMatchColl->At(fiNbHits);
              */
              CbmMatch* digiMatch = new ((*fTofDigiMatchColl)[fiNbHits]) CbmMatch();
              for (UInt_t i = 0; i < vDigiIndRef.size(); i++) {
                Double_t dTot = ((CbmTofDigi*) (&(fTofCalDigiVec->at(vDigiIndRef.at(i)))))->GetTot();
                digiMatch->AddLink(CbmLink(dTot, vDigiIndRef.at(i), fiOutputTreeEntry, fiFileIndex));
              }
 
              fiNbHits++;
              // For Histogramming
              fviClusterSize[iSmType][iRpc].push_back(iNbChanInHit);
              //fviTrkMul[iSmType][iRpc].push_back( vPtsRef.size() );
              fvdX[iSmType][iRpc].push_back(dWeightedPosX);
              fvdY[iSmType][iRpc].push_back(dWeightedPosY);
              /*
                fvdDifX[iSmType][iRpc].push_back( vPtsRef[0]->GetX() - dWeightedPosX);
                fvdDifY[iSmType][iRpc].push_back( vPtsRef[0]->GetY() - dWeightedPosY);
                fvdDifCh[iSmType][iRpc].push_back( fGeoHandler->GetCell( vPtsRef[0]->GetDetectorID() ) -1 -iLastChan );
              */
              //vPtsRef.clear();
              vDigiIndRef.clear();
            }  // else of if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
          }    // if( 0 < iNbChanInHit)
        }      // for each sm/rpc pair
    }          // for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ )
  }
 
  return kTRUE;
}
 
/************************************************************************************/
Bool_t CbmDeviceHitBuilderTof::MergeClusters() { return kTRUE; }
 
void CbmDeviceHitBuilderTof::LH_store(Int_t iSmType, Int_t iSm, Int_t iRpc, Int_t iChm, CbmTofHit* pHit)
{
 
  if (fvLastHits[iSmType][iSm][iRpc][iChm].size() == 0) fvLastHits[iSmType][iSm][iRpc][iChm].push_back(pHit);
  else {
    Double_t dLastTime = pHit->GetTime();
    if (dLastTime >= fvLastHits[iSmType][iSm][iRpc][iChm].back()->GetTime()) {
      fvLastHits[iSmType][iSm][iRpc][iChm].push_back(pHit);
      LOG(debug) << Form(" Store LH from Ev  %8.0f for TSRC %d%d%d%d, size %lu, addr 0x%08x, "
                         "time %f, dt %f",
                         fdEvent, iSmType, iSm, iRpc, iChm, fvLastHits[iSmType][iSm][iRpc][iChm].size(),
                         pHit->GetAddress(), dLastTime,
                         dLastTime - fvLastHits[iSmType][iSm][iRpc][iChm].front()->GetTime());
    }
    else {
      if (dLastTime
          >= fvLastHits[iSmType][iSm][iRpc][iChm].front()->GetTime()) {  // hit has to be inserted in the proper place
        std::list<CbmTofHit*>::iterator it;
        for (it = fvLastHits[iSmType][iSm][iRpc][iChm].begin(); it != fvLastHits[iSmType][iSm][iRpc][iChm].end(); ++it)
          if ((*it)->GetTime() > dLastTime) break;
        fvLastHits[iSmType][iSm][iRpc][iChm].insert(--it, pHit);
        Double_t deltaTime = dLastTime - (*it)->GetTime();
        LOG(debug) << Form("Hit inserted into LH from Ev  %8.0f for TSRC "
                           "%d%d%d%d, size %lu, addr 0x%08x, delta time %f  ",
                           fdEvent, iSmType, iSm, iRpc, iChm, fvLastHits[iSmType][iSm][iRpc][iChm].size(),
                           pHit->GetAddress(), deltaTime);
      }
      else {  // this hit is first
        Double_t deltaTime = dLastTime - fvLastHits[iSmType][iSm][iRpc][iChm].front()->GetTime();
        LOG(debug) << Form("first LH from Ev  %8.0f for TSRC %d%d%d%d, size "
                           "%lu, addr 0x%08x, delta time %f ",
                           fdEvent, iSmType, iSm, iRpc, iChm, fvLastHits[iSmType][iSm][iRpc][iChm].size(),
                           pHit->GetAddress(), deltaTime);
        if (deltaTime == 0.) {
          // remove hit, otherwise double entry?
          pHit->Delete();
        }
        else {
          fvLastHits[iSmType][iSm][iRpc][iChm].push_front(pHit);
        }
      }
    }
  }
}
 
void CbmDeviceHitBuilderTof::CheckLHMemory()
{
  if ((Int_t) fvLastHits.size() != fDigiBdfPar->GetNbSmTypes())
    LOG(error) << Form("Inconsistent LH Smtype size %lu, %d ", fvLastHits.size(), fDigiBdfPar->GetNbSmTypes());
 
  for (Int_t iSmType = 0; iSmType < fDigiBdfPar->GetNbSmTypes(); iSmType++) {
    if ((Int_t) fvLastHits[iSmType].size() != fDigiBdfPar->GetNbSm(iSmType))
      LOG(error) << Form("Inconsistent LH Sm size %lu, %d T %d", fvLastHits[iSmType].size(),
                         fDigiBdfPar->GetNbSm(iSmType), iSmType);
    for (Int_t iSm = 0; iSm < fDigiBdfPar->GetNbSm(iSmType); iSm++) {
      if ((Int_t) fvLastHits[iSmType][iSm].size() != fDigiBdfPar->GetNbRpc(iSmType))
        LOG(error) << Form("Inconsistent LH Rpc size %lu, %d TS %d%d ", fvLastHits[iSmType][iSm].size(),
                           fDigiBdfPar->GetNbRpc(iSmType), iSmType, iSm);
      for (Int_t iRpc = 0; iRpc < fDigiBdfPar->GetNbRpc(iSmType); iRpc++) {
        if ((Int_t) fvLastHits[iSmType][iSm][iRpc].size() != fDigiBdfPar->GetNbChan(iSmType, iRpc))
          LOG(error) << Form("Inconsistent LH RpcChannel size %lu, %d TSR %d%d%d ",
                             fvLastHits[iSmType][iSm][iRpc].size(), fDigiBdfPar->GetNbChan(iSmType, iRpc), iSmType, iSm,
                             iRpc);
        for (Int_t iCh = 0; iCh < fDigiBdfPar->GetNbChan(iSmType, iRpc); iCh++)
          if (fvLastHits[iSmType][iSm][iRpc][iCh].size() > 0) {
            CbmTofDetectorInfo xDetInfo(ECbmModuleId::kTof, iSmType, iSm, iRpc, 0, iCh);
            Int_t iAddr = fTofId->SetDetectorInfo(xDetInfo);
            if (fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetAddress() != iAddr)
              LOG(error) << Form("Inconsistent address for Ev %8.0f in list of size %lu for "
                                 "TSRC %d%d%d%d: 0x%08x, time  %f",
                                 fdEvent, fvLastHits[iSmType][iSm][iRpc][iCh].size(), iSmType, iSm, iRpc, iCh,
                                 fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetAddress(),
                                 fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetTime());
          }
      }
    }
  }
  LOG(debug) << Form("LH check passed for event %8.0f", fdEvent);
}
 
void CbmDeviceHitBuilderTof::CleanLHMemory()
{
  if ((Int_t) fvLastHits.size() != fDigiBdfPar->GetNbSmTypes())
    LOG(error) << Form("Inconsistent LH Smtype size %lu, %d ", fvLastHits.size(), fDigiBdfPar->GetNbSmTypes());
  for (Int_t iSmType = 0; iSmType < fDigiBdfPar->GetNbSmTypes(); iSmType++) {
    if ((Int_t) fvLastHits[iSmType].size() != fDigiBdfPar->GetNbSm(iSmType))
      LOG(error) << Form("Inconsistent LH Sm size %lu, %d T %d", fvLastHits[iSmType].size(),
                         fDigiBdfPar->GetNbSm(iSmType), iSmType);
    for (Int_t iSm = 0; iSm < fDigiBdfPar->GetNbSm(iSmType); iSm++) {
      if ((Int_t) fvLastHits[iSmType][iSm].size() != fDigiBdfPar->GetNbRpc(iSmType))
        LOG(error) << Form("Inconsistent LH Rpc size %lu, %d TS %d%d ", fvLastHits[iSmType][iSm].size(),
                           fDigiBdfPar->GetNbRpc(iSmType), iSmType, iSm);
      for (Int_t iRpc = 0; iRpc < fDigiBdfPar->GetNbRpc(iSmType); iRpc++) {
        if ((Int_t) fvLastHits[iSmType][iSm][iRpc].size() != fDigiBdfPar->GetNbChan(iSmType, iRpc))
          LOG(error) << Form("Inconsistent LH RpcChannel size %lu, %d TSR %d%d%d ",
                             fvLastHits[iSmType][iSm][iRpc].size(), fDigiBdfPar->GetNbChan(iSmType, iRpc), iSmType, iSm,
                             iRpc);
        for (Int_t iCh = 0; iCh < fDigiBdfPar->GetNbChan(iSmType, iRpc); iCh++)
          while (fvLastHits[iSmType][iSm][iRpc][iCh].size() > 0) {
            CbmTofDetectorInfo xDetInfo(ECbmModuleId::kTof, iSmType, iSm, iRpc, 0, iCh);
            Int_t iAddr = fTofId->SetDetectorInfo(xDetInfo);
            if (fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetAddress() != iAddr)
              LOG(error) << Form("Inconsistent address for Ev %8.0f in list of size %lu for "
                                 "TSRC %d%d%d%d: 0x%08x, time  %f",
                                 fdEvent, fvLastHits[iSmType][iSm][iRpc][iCh].size(), iSmType, iSm, iRpc, iCh,
                                 fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetAddress(),
                                 fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetTime());
            fvLastHits[iSmType][iSm][iRpc][iCh].front()->Delete();
            fvLastHits[iSmType][iSm][iRpc][iCh].pop_front();
          }
      }
    }
  }
  LOG(info) << Form("LH cleaning done after %8.0f events", fdEvent);
}
 
Bool_t CbmDeviceHitBuilderTof::AddNextChan(Int_t iSmType, Int_t iSm, Int_t iRpc, Int_t iLastChan, Double_t dLastPosX,
                                           Double_t dLastPosY, Double_t dLastTime, Double_t dLastTotS)
{
  //Int_t iNbSm  = fDigiBdfPar->GetNbSm(  iSmType);
  Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
  Int_t iNbCh  = fDigiBdfPar->GetNbChan(iSmType, iRpc);
  //  Int_t iChType = fDigiBdfPar->GetChanType( iSmType, iRpc );
 
  Int_t iCh = iLastChan + 1;
  LOG(debug) << Form("Inspect channel TSRC %d%d%d%d at time %f, pos %f, size ", iSmType, iSm, iRpc, iCh, dLastTime,
                     dLastPosY)
             << fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size();
  if (iCh == iNbCh) return kFALSE;
  if (0 == fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size()) return kFALSE;
  //if( 0 < fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh].size() )
  //   fhNbDigiPerChan->Fill( fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh].size() );
  if (1 < fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size()) {
    Bool_t AddedHit = kFALSE;
    for (Int_t i1 = 0; i1 < (Int_t) fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size() - 1; i1++) {
      if (AddedHit) break;
      Int_t i2 = i1 + 1;
      while (!AddedHit && i2 < (Int_t) fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size()) {
        // LOG(debug)<<"check digi pair "<<i1<<","<<i2<<" with size "<<fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh].size();
 
        if ((fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][i1])->GetSide()
            == (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][i2])->GetSide()) {
          i2++;
          continue;
        }  // endif same side
           // 2 Digis, both sides present
        CbmTofDigi* xDigiA = fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][i1];
        CbmTofDigi* xDigiB = fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][i2];
        Double_t dTime     = 0.5 * (xDigiA->GetTime() + xDigiB->GetTime());
        if (TMath::Abs(dTime - dLastTime) < fdMaxTimeDist) {
          CbmTofDetectorInfo xDetInfo(ECbmModuleId::kTof, iSmType, iSm, iRpc, 0, iCh);
          Int_t iChId  = fTofId->SetDetectorInfo(xDetInfo);
          fChannelInfo = fDigiPar->GetCell(iChId);
          gGeoManager->FindNode(fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ());
 
          Double_t dTimeDif = xDigiA->GetTime() - xDigiB->GetTime();
          Double_t dPosY    = 0.;
          if (1 == xDigiA->GetSide()) dPosY = fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc) * dTimeDif * 0.5;
          else
            dPosY = -fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc) * dTimeDif * 0.5;
 
          if (TMath::Abs(dPosY - dLastPosY) < fdMaxSpaceDist) {  // append digi pair to current cluster
 
            Double_t dNClHits = (Double_t)(vDigiIndRef.size() / 2);
            Double_t dPosX    = ((Double_t)(-iNbCh / 2 + iCh) + 0.5) * fChannelInfo->GetSizex();
            Double_t dTotS    = xDigiA->GetTot() + xDigiB->GetTot();
            Double_t dNewTotS = (dLastTotS + dTotS);
            dLastPosX         = (dLastPosX * dLastTotS + dPosX * dTotS) / dNewTotS;
            dLastPosY         = (dLastPosY * dLastTotS + dPosY * dTotS) / dNewTotS;
            dLastTime         = (dLastTime * dLastTotS + dTime * dTotS) / dNewTotS;
            dLastTotS         = dNewTotS;
            // attach selected digis from pool
            Int_t Ind1 = fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][i1];
            Int_t Ind2 = fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][i2];
            vDigiIndRef.push_back(Ind1);
            vDigiIndRef.push_back(Ind2);
            // remove selected digis from pool
            fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin()
                                                               + i1);
            fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
              fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin() + i1);
 
            std::vector<int>::iterator it;
            it = find(fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin(),
                      fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].end(), Ind2);
            if (it != fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].end()) {
              auto ipos = it - fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin();
              LOG(debug) << "Found i2 " << i2 << " with Ind2 " << Ind2 << " at position " << ipos;
              fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin()
                                                                 + ipos);
              fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
                fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin() + ipos);
            }
            else {
              LOG(error) << " Did not find  i2 " << i2 << " with Ind2 " << Ind2;
            }
 
            //if(iCh == iNbCh-1) break;  //Last strip reached
            if (iCh != (iNbCh - 1)
                && AddNextChan(iSmType, iSm, iRpc, iCh, dLastPosX, dLastPosY, dLastTime, dLastTotS)) {
              LOG(debug) << "Added Strip " << iCh << " to cluster of size " << dNClHits;
              return kTRUE;  // signal hit was already added
            }
            AddedHit = kTRUE;
          }  //TMath::Abs(dPosY - dLastPosY) < fdMaxSpaceDist
        }    //TMath::Abs(dTime-dLastTime)<fdMaxTimeDist)
        i2++;
      }  //  while(i2 < fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh].size()-1 )
    }    // end for i1
  }      // end if size
  Double_t hitpos_local[3] = {3 * 0.};
  hitpos_local[0]          = dLastPosX;
  hitpos_local[1]          = dLastPosY;
  hitpos_local[2]          = 0.;
  Double_t hitpos[3];
 
  TGeoNode* cNode = gGeoManager->GetCurrentNode();
  if (NULL == cNode) {  // Transformation matrix not available !!!??? - Check
    LOG(error) << Form("Node at (%6.1f,%6.1f,%6.1f) : %p", fChannelInfo->GetX(), fChannelInfo->GetY(),
                       fChannelInfo->GetZ(), cNode);
    ChangeState(fair::mq::Transition::Stop);
  }
 
  /*TGeoHMatrix* cMatrix = */ gGeoManager->GetCurrentMatrix();
  gGeoManager->LocalToMaster(hitpos_local, hitpos);
  TVector3 hitPos(hitpos[0], hitpos[1], hitpos[2]);
  TVector3 hitPosErr(0.5, 0.5, 0.5);  // FIXME including positioning uncertainty
  Int_t iChm = floor(dLastPosX / fChannelInfo->GetSizex()) + iNbCh / 2;
  if (iChm < 0) iChm = 0;
  if (iChm > iNbCh - 1) iChm = iNbCh - 1;
  Int_t iDetId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, iChm, 0, iSmType);
 
  //  Int_t iNbChanInHit=vDigiIndRef.size()/2;
  fviClusterMul[iSmType][iSm][iRpc]++;
  /*
  LOG(debug)<<"Save A-Hit "
            << Form("%2d %2d 0x%08x %3d t %f, y %f ",
                    fiNbHits,iNbChanInHit,iDetId,iLastChan,dLastTime,dLastPosY)
            <<", DigiSize: "<<vDigiIndRef.size();
  for (UInt_t i=0; i<vDigiIndRef.size();i++){
    LOG(debug)<<"DigiIndRef "<<i<<" "<<vDigiIndRef.at(i)<<"(M"<<fviClusterMul[iSmType][iSm][iRpc]<<")";
  }
  */
  CbmTofHit* pHit = new CbmTofHit(iDetId, hitPos,
                                  hitPosErr,  //local detector coordinates
                                  fiNbHits,   // this number is used as reference!!
                                  dLastTime,
                                  vDigiIndRef.size(),  // number of linked digis =  2*CluSize
                                  //vPtsRef.size(), // flag  = number of TofPoints generating the cluster
                                  Int_t(dLastTotS * 10.));  //channel -> Tot
  // output hit
  new ((*fTofHitsColl)[fiNbHits]) CbmTofHit(*pHit);
  if (fdMemoryTime > 0.) {  // memorize hit
    LH_store(iSmType, iSm, iRpc, iChm, pHit);
  }
  else {
    pHit->Delete();
  }
  CbmMatch* digiMatch = new ((*fTofDigiMatchColl)[fiNbHits]) CbmMatch();
  for (Int_t i = 0; i < (Int_t) vDigiIndRef.size(); i++) {
    Double_t dTot = ((CbmTofDigi*) (&(fTofCalDigiVec->at(vDigiIndRef.at(i)))))->GetTot();
    digiMatch->AddLink(CbmLink(dTot, vDigiIndRef.at(i), fiOutputTreeEntry, fiFileIndex));
  }
  fiNbHits++;
  vDigiIndRef.clear();
 
  return kTRUE;
}
 
static Double_t f1_xboxe(double* x, double* par)
{
  double xx    = x[0];
  double wx    = 1. - par[4] * TMath::Power(xx + par[5], 2);
  double xboxe = par[0] * 0.25 * (1. + TMath::Erf((xx + par[1] - par[3]) / par[2]))
                 * (1. + TMath::Erf((-xx + par[1] + par[3]) / par[2]));
  return xboxe * wx;
}
 
void CbmDeviceHitBuilderTof::fit_ybox(const char* hname)
{
  TH1* h1;
  h1 = (TH1*) gROOT->FindObjectAny(hname);
  if (NULL != h1) { fit_ybox(h1, 0.); }
}
 
void CbmDeviceHitBuilderTof::fit_ybox(TH1* h1, Double_t ysize)
{
  Double_t* fpar = NULL;
  fit_ybox(h1, ysize, fpar);
}
 
void CbmDeviceHitBuilderTof::fit_ybox(TH1* h1, Double_t ysize, Double_t* fpar = NULL)
{
  TAxis* xaxis  = h1->GetXaxis();
  Double_t Ymin = xaxis->GetXmin();
  Double_t Ymax = xaxis->GetXmax();
  TF1* f1       = new TF1("YBox", f1_xboxe, Ymin, Ymax, 6);
  Double_t yini = (h1->GetMaximum() + h1->GetMinimum()) * 0.5;
  if (ysize == 0.) ysize = Ymax * 0.8;
  f1->SetParameters(yini, ysize * 0.5, 1., 0., 0., 0.);
  //  f1->SetParLimits(1,ysize*0.8,ysize*1.2);
  f1->SetParLimits(2, 0.2, 3.);
  f1->SetParLimits(3, -4., 4.);
  if (fpar != NULL) {
    Double_t fp[4];
    for (Int_t i = 0; i < 4; i++)
      fp[i] = *fpar++;
    for (Int_t i = 0; i < 4; i++)
      f1->SetParameter(2 + i, fp[i]);
    LOG(debug) << "Ini Fpar for " << h1->GetName() << " with "
               << Form(" %6.3f %6.3f %6.3f %6.3f ", fp[0], fp[1], fp[2], fp[3]);
  }
 
  h1->Fit("YBox", "Q");
 
  double res[10];
  double err[10];
  res[9] = f1->GetChisquare();
 
  for (int i = 0; i < 6; i++) {
    res[i] = f1->GetParameter(i);
    err[i] = f1->GetParError(i);
    //cout << " FPar "<< i << ": " << res[i] << ", " << err[i] << endl;
  }
  LOG(debug) << "YBox Fit of " << h1->GetName() << " ended with chi2 = " << res[9]
             << Form(", strip length %7.2f +/- %5.2f, position resolution "
                     "%7.2f +/- %5.2f at y_cen = %7.2f +/- %5.2f",
                     2. * res[1], 2. * err[1], res[2], err[2], res[3], err[3]);
}
 
Bool_t CbmDeviceHitBuilderTof::LoadGeometry()
{
  LOG(info) << "LoadGeometry starting for  " << fDigiBdfPar->GetNbDet() << " described detectors, "
            << fDigiPar->GetNrOfModules() << " geometrically known modules ";
 
  //gGeoManager->Export("HitBuilder.loadgeo.root");  // write current status to file
 
  Int_t iNrOfCells = fDigiPar->GetNrOfModules();
  LOG(info) << "Digi Parameter container contains " << iNrOfCells << " cells.";
 
  for (Int_t icell = 0; icell < iNrOfCells; ++icell) {
 
    Int_t cellId = fDigiPar->GetCellId(icell);  // cellId is assigned in CbmTofCreateDigiPar
    fChannelInfo = fDigiPar->GetCell(cellId);
 
    Int_t smtype  = fGeoHandler->GetSMType(cellId);
    Int_t smodule = fGeoHandler->GetSModule(cellId);
    Int_t module  = fGeoHandler->GetCounter(cellId);
    Int_t cell    = fGeoHandler->GetCell(cellId);
 
    Double_t x  = fChannelInfo->GetX();
    Double_t y  = fChannelInfo->GetY();
    Double_t z  = fChannelInfo->GetZ();
    Double_t dx = fChannelInfo->GetSizex();
    Double_t dy = fChannelInfo->GetSizey();
    LOG(debug) << "-I- InitPar " << icell << " Id: " << Form("0x%08x", cellId) << " " << cell << " tmcs: " << smtype
               << " " << smodule << " " << module << " " << cell << " x=" << Form("%6.2f", x)
               << " y=" << Form("%6.2f", y) << " z=" << Form("%6.2f", z) << " dx=" << dx << " dy=" << dy;
 
    TGeoNode* fNode =  // prepare local->global trafo
      gGeoManager->FindNode(fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ());
 
    if (NULL == fNode) {  // Transformation matrix not available !!!??? - Check
      LOG(error) << Form("Node at (%6.1f,%6.1f,%6.1f) : %p", fChannelInfo->GetX(), fChannelInfo->GetY(),
                         fChannelInfo->GetZ(), fNode);
      ChangeState(fair::mq::Transition::Stop);
    }
    if (icell == 0) {
      TGeoHMatrix* cMatrix = gGeoManager->GetCurrentMatrix();
      //cNode->Print();
      cMatrix->Print();
    }
  }
 
  Int_t iNbDet = fDigiBdfPar->GetNbDet();
  fvDeadStrips.resize(iNbDet);
  fvPulserOffset.resize(iNbDet);
  fvPulserTimes.resize(iNbDet);
 
  for (Int_t iDetIndx = 0; iDetIndx < iNbDet; iDetIndx++) {
    Int_t iUniqueId = fDigiBdfPar->GetDetUId(iDetIndx);
    Int_t iSmType   = CbmTofAddress::GetSmType(iUniqueId);
    Int_t iSmId     = CbmTofAddress::GetSmId(iUniqueId);
    Int_t iRpcId    = CbmTofAddress::GetRpcId(iUniqueId);
    LOG(info) << " DetIndx " << iDetIndx << "(" << iNbDet << "), SmType " << iSmType << ", SmId " << iSmId << ", RpcId "
              << iRpcId << " => UniqueId " << Form("0x%08x ", iUniqueId)
              << Form(" Svel %6.6f, DeadStrips 0x%08x ", fDigiBdfPar->GetSigVel(iSmType, iSmId, iRpcId),
                      fvDeadStrips[iDetIndx]);
 
    Int_t iCell = -1;
    while (kTRUE) {
      Int_t iUCellId = CbmTofAddress::GetUniqueAddress(iSmId, iRpcId, ++iCell, 0, iSmType);
      fChannelInfo   = fDigiPar->GetCell(iUCellId);
      if (NULL == fChannelInfo) break;
    }
 
    fvPulserOffset[iDetIndx].resize(2);  // provide vector for both sides
    for (Int_t i = 0; i < 2; i++)
      fvPulserOffset[iDetIndx][i] = 0.;  // initialize
    fvPulserTimes[iDetIndx].resize(2);   // provide vector for both sides
  }
 
  //   return kTRUE;
 
  Int_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
 
  if (kTRUE == fDigiBdfPar->UseExpandedDigi()) {
    fStorDigi.resize(iNbSmTypes);
    fStorDigiInd.resize(iNbSmTypes);
    fviClusterSize.resize(iNbSmTypes);
    fviTrkMul.resize(iNbSmTypes);
    fvdX.resize(iNbSmTypes);
    fvdY.resize(iNbSmTypes);
    fvdDifX.resize(iNbSmTypes);
    fvdDifY.resize(iNbSmTypes);
    fvdDifCh.resize(iNbSmTypes);
    fviClusterMul.resize(iNbSmTypes);
    fvLastHits.resize(iNbSmTypes);
 
    for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
      Int_t iNbSm  = fDigiBdfPar->GetNbSm(iSmType);
      Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
      fStorDigi[iSmType].resize(iNbSm * iNbRpc);
      fStorDigiInd[iSmType].resize(iNbSm * iNbRpc);
      fviClusterSize[iSmType].resize(iNbRpc);
      fviTrkMul[iSmType].resize(iNbRpc);
      fvdX[iSmType].resize(iNbRpc);
      fvdY[iSmType].resize(iNbRpc);
      fvdDifX[iSmType].resize(iNbRpc);
      fvdDifY[iSmType].resize(iNbRpc);
      fvdDifCh[iSmType].resize(iNbRpc);
      for (Int_t iSm = 0; iSm < iNbSm; iSm++) {
        fviClusterMul[iSmType].resize(iNbSm);
        fvLastHits[iSmType].resize(iNbSm);
        for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
          fviClusterMul[iSmType][iSm].resize(iNbRpc);
          fvLastHits[iSmType][iSm].resize(iNbRpc);
          Int_t iNbChan = fDigiBdfPar->GetNbChan(iSmType, iRpc);
          if (iNbChan == 0) {
            LOG(warn) << "LoadGeometry: StoreDigi without channels "
                      << Form("SmTy %3d, Sm %3d, NbRpc %3d, Rpc, %3d ", iSmType, iSm, iNbRpc, iRpc);
          }
          fStorDigi[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
          fStorDigiInd[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
          fvLastHits[iSmType][iSm][iRpc].resize(iNbChan);
        }  // for( Int_t iRpc = 0; iRpc < iNbRpc; iRpc++ )
      }    // for( Int_t iSm = 0; iSm < iNbSm; iSm++ )
    }      // for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ )
  }        // if( kTRUE == fDigiBdfPar->UseExpandedDigi() )
  else {
    fStorDigi.resize(iNbSmTypes);
    fStorDigiInd.resize(iNbSmTypes);
    fviClusterSize.resize(iNbSmTypes);
    fviTrkMul.resize(iNbSmTypes);
    fvdX.resize(iNbSmTypes);
    fvdY.resize(iNbSmTypes);
    fvdDifX.resize(iNbSmTypes);
    fvdDifY.resize(iNbSmTypes);
    fvdDifCh.resize(iNbSmTypes);
    for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
      Int_t iNbSm  = fDigiBdfPar->GetNbSm(iSmType);
      Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
      fStorDigi[iSmType].resize(iNbSm * iNbRpc);
      fStorDigiInd[iSmType].resize(iNbSm * iNbRpc);
      fviClusterSize[iSmType].resize(iNbRpc);
      fviTrkMul[iSmType].resize(iNbRpc);
      fvdX[iSmType].resize(iNbRpc);
      fvdY[iSmType].resize(iNbRpc);
      fvdDifX[iSmType].resize(iNbRpc);
      fvdDifY[iSmType].resize(iNbRpc);
      fvdDifCh[iSmType].resize(iNbRpc);
      for (Int_t iSm = 0; iSm < iNbSm; iSm++) {
        for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
          Int_t iNbChan = fDigiBdfPar->GetNbChan(iSmType, iRpc);
          fStorDigi[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
          fStorDigiInd[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
        }  // for( Int_t iRpc = 0; iRpc < iNbRpc; iRpc++ )
      }    // for( Int_t iSm = 0; iSm < iNbSm; iSm++ )
    }      // for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ )
  }        // else of if( kTRUE == fDigiBdfPar->UseExpandedDigi() )
  return kTRUE;
}
 
Bool_t CbmDeviceHitBuilderTof::FillDigiStor()
{
 
  // Loop over the digis array and store the Digis in separate vectors for
  // each RPC modules
 
  CbmTofDigi* pDigi;
 
  Int_t iNbTofDigi = (Int_t) fTofCalDigiVec->size();
  //fTofCalDigisColl->GetEntriesFast();
  for (Int_t iDigInd = 0; iDigInd < iNbTofDigi; iDigInd++) {
    //pDigi = (CbmTofDigi*) fTofCalDigisColl->At(iDigInd);
    pDigi = &(fTofCalDigiVec->at(iDigInd));
 
    /*
    LOG(info)<<"AC " // After Calibration
              <<Form("0x%08x",pDigi->GetAddress())<<" TSRC "
              <<pDigi->GetType()
              <<pDigi->GetSm()
              <<pDigi->GetRpc()
              <<Form("%2d",(Int_t)pDigi->GetChannel())<<" "
              <<pDigi->GetSide()<<" "
              <<Form("%f",pDigi->GetTime())<<" "
              <<pDigi->GetTot();
    */
    if (fDigiBdfPar->GetNbSmTypes() > pDigi->GetType()  // prevent crash due to misconfiguration
        && fDigiBdfPar->GetNbSm(pDigi->GetType()) > pDigi->GetSm()
        && fDigiBdfPar->GetNbRpc(pDigi->GetType()) > pDigi->GetRpc()
        && fDigiBdfPar->GetNbChan(pDigi->GetType(), pDigi->GetRpc()) > pDigi->GetChannel()) {
      fStorDigi[pDigi->GetType()][pDigi->GetSm() * fDigiBdfPar->GetNbRpc(pDigi->GetType()) + pDigi->GetRpc()]
               [pDigi->GetChannel()]
                 .push_back(pDigi);
      fStorDigiInd[pDigi->GetType()][pDigi->GetSm() * fDigiBdfPar->GetNbRpc(pDigi->GetType()) + pDigi->GetRpc()]
                  [pDigi->GetChannel()]
                    .push_back(iDigInd);
    }
    else {
      LOG(info) << "Skip2 Digi "
                << " Type " << pDigi->GetType() << " " << fDigiBdfPar->GetNbSmTypes() << " Sm " << pDigi->GetSm() << " "
                << fDigiBdfPar->GetNbSm(pDigi->GetType()) << " Rpc " << pDigi->GetRpc() << " "
                << fDigiBdfPar->GetNbRpc(pDigi->GetType()) << " Ch " << pDigi->GetChannel() << " "
                << fDigiBdfPar->GetNbChan(pDigi->GetType(), 0);
    }
  }  // for( Int_t iDigInd = 0; iDigInd < nTofDigi; iDigInd++ )
  return kTRUE;
}
 
Bool_t CbmDeviceHitBuilderTof::SendHits()
{
  //Output Log
  for (Int_t iHit = 0; iHit < fiNbHits; iHit++) {
    CbmTofHit* pHit = (CbmTofHit*) fTofHitsColl->At(iHit);
    LOG(debug) << Form("Found Hit %d, addr 0x%08x, X %6.2f, Y %6.2f Z %6.2f T %6.2f CLS %d", iHit, pHit->GetAddress(),
                       pHit->GetX(), pHit->GetY(), pHit->GetZ(), pHit->GetTime(), pHit->GetFlag());
  }
  // prepare output hit vector
  std::vector<CbmTofHit*> vhit;
  vhit.resize(fiNbHits);
  for (Int_t iHit = 0; iHit < fiNbHits; iHit++) {
    CbmTofHit* pHit = (CbmTofHit*) fTofHitsColl->At(iHit);
    vhit[iHit]      = pHit;
  }
 
  // prepare output string streams
  std::stringstream ossE;
  boost::archive::binary_oarchive oaE(ossE);
  oaE << fEventHeader;
  std::string* strMsgE = new std::string(ossE.str());
 
  std::stringstream oss;
  boost::archive::binary_oarchive oa(oss);
  oa << vhit;
  std::string* strMsg = new std::string(oss.str());
 
  FairMQParts parts;
  parts.AddPart(NewMessage(
    const_cast<char*>(strMsgE->c_str()),  // data
    strMsgE->length(),                    // size
    [](void*, void* object) { delete static_cast<std::string*>(object); },
    strMsgE));  // object that manages the data
 
  parts.AddPart(NewMessage(
    const_cast<char*>(strMsg->c_str()),  // data
    strMsg->length(),                    // size
    [](void*, void* object) { delete static_cast<std::string*>(object); },
    strMsg));  // object that manages the data
 
  if (Send(parts, "tofhits")) {
    LOG(error) << "Problem sending data ";
    return false;
  }
 
  return kTRUE;
}
 
Bool_t CbmDeviceHitBuilderTof::SendAll() { return kTRUE; }
 
Bool_t CbmDeviceHitBuilderTof::FillHistos()
{
  Int_t iNbTofHits = fTofHitsColl->GetEntriesFast();
  CbmTofHit* pHit;
  //gGeoManager->SetTopVolume( gGeoManager->FindVolumeFast("tof_v14a") );
  gGeoManager->CdTop();
 
  if (0 < iNbTofHits) {
    Bool_t BSel[iNSel];
    Double_t dTTrig[iNSel];
    CbmTofHit* pTrig[iNSel];
    Double_t ddXdZ[iNSel];
    Double_t ddYdZ[iNSel];
    Double_t dSel2dXdYMin[iNSel];
 
    Int_t iBeamRefMul    = 0;
    Int_t iBeamAddRefMul = 0;
 
    if (0 < iNSel) {  // check software triggers
 
      LOG(debug) << "FillHistos() for " << iNSel << " triggers"
                 << ", Dut " << fDutId << fDutSm << fDutRpc << Form(", 0x%08x", fDutAddr) << ", Sel " << fSelId
                 << fSelSm << fSelRpc << Form(", 0x%08x", fSelAddr) << ", Sel2 " << fSel2Id << fSel2Sm << fSel2Rpc
                 << Form(", 0x%08x", fSel2Addr);
      /*
      LOG(debug) <<"FillHistos: Muls: "
                 <<fviClusterMul[fDutId][fDutSm][fDutRpc]
                 <<", "<<fviClusterMul[fSelId][fSelSm][fSelRpc]
        ;
      */
      // monitor multiplicities
      Int_t iNbDet = fDigiBdfPar->GetNbDet();
      for (Int_t iDetIndx = 0; iDetIndx < iNbDet; iDetIndx++) {
        Int_t iDetId  = fviDetId[iDetIndx];
        Int_t iSmType = CbmTofAddress::GetSmType(iDetId);
        Int_t iSm     = CbmTofAddress::GetSmId(iDetId);
        Int_t iRpc    = CbmTofAddress::GetRpcId(iDetId);
        //LOG(info) << Form(" indx %d, Id 0x%08x, TSR %d %d %d", iDetIndx, iDetId, iSmType, iSm, iRpc)
        //          ;
        if (NULL != fhRpcCluMul[iDetIndx]) fhRpcCluMul[iDetIndx]->Fill(fviClusterMul[iSmType][iSm][iRpc]);  //
      }
 
      // do input distributions first
      for (Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++) {
        pHit = (CbmTofHit*) fTofHitsColl->At(iHitInd);
        if (NULL == pHit) continue;
        if (StartAnalysisTime == 0.) {
          StartAnalysisTime = pHit->GetTime();
          LOG(info) << "StartAnalysisTime set to " << StartAnalysisTime << " ns. ";
        }
        Int_t iDetId                          = (pHit->GetAddress() & DetMask);
        std::map<UInt_t, UInt_t>::iterator it = fDetIdIndexMap.find(iDetId);
        if (it == fDetIdIndexMap.end()) continue;  // continue for invalid detector index
        Int_t iDetIndx = it->second;               //fDetIdIndexMap[iDetId];
 
        Int_t iSmType = CbmTofAddress::GetSmType(iDetId);
        Int_t iSm     = CbmTofAddress::GetSmId(iDetId);
        Int_t iRpc    = CbmTofAddress::GetRpcId(iDetId);
        Int_t iCh     = CbmTofAddress::GetChannelId(pHit->GetAddress());
 
        Double_t dTimeAna = (pHit->GetTime() - StartAnalysisTime) / 1.E9;
        //LOG(debug)<<"TimeAna "<<StartAnalysisTime<<", "<< pHit->GetTime()<<", "<<dTimeAna;
        fhRpcCluRate[iDetIndx]->Fill(dTimeAna);
 
        if (fdMemoryTime > 0. && fvLastHits[iSmType][iSm][iRpc][iCh].size() == 0)
          LOG(error) << Form(" <E> hit not stored in memory for TSRC %d%d%d%d", iSmType, iSm, iRpc, iCh);
        //CheckLHMemory();
 
        if (fvLastHits[iSmType][iSm][iRpc][iCh].size() > 1) {  // check for outdated hits
          //std::list<CbmTofHit *>::iterator it0=fvLastHits[iSmType][iSm][iRpc][iCh].begin();
          //std::list<CbmTofHit *>::iterator itL=fvLastHits[iSmType][iSm][iRpc][iCh].end();
          //CbmTofHit* pH0 = *it0;
          //CbmTofHit* pHL = *(--itL);
          CbmTofHit* pH0 = fvLastHits[iSmType][iSm][iRpc][iCh].front();
          CbmTofHit* pHL = fvLastHits[iSmType][iSm][iRpc][iCh].back();
          if (pH0->GetTime() > pHL->GetTime())
            LOG(warn) << Form("Invalid time ordering in ev %8.0f in list of "
                              "size %lu for TSRC %d%d%d%d: Delta t %f  ",
                              fdEvent, fvLastHits[iSmType][iSm][iRpc][iCh].size(), iSmType, iSm, iRpc, iCh,
                              pHL->GetTime() - pH0->GetTime());
          //while( (*((std::list<CbmTofHit *>::iterator) fvLastHits[iSmType][iSm][iRpc][iCh].begin()))->GetTime()+fdMemoryTime < pHit->GetTime()
          while (fvLastHits[iSmType][iSm][iRpc][iCh].size() > 2.
                 || fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetTime() + fdMemoryTime < pHit->GetTime()) {
            LOG(debug) << " pop from list size " << fvLastHits[iSmType][iSm][iRpc][iCh].size()
                       << Form(" outdated hits for ev %8.0f in TSRC %d%d%d%d", fdEvent, iSmType, iSm, iRpc, iCh)
                       << Form(" with tHit - tLast %f ",
                               pHit->GetTime() - fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetTime())
              //(*((std::list<CbmTofHit *>::iterator) fvLastHits[iSmType][iSm][iRpc][iCh].begin()))->GetTime())
              ;
            if (fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetAddress() != pHit->GetAddress())
              LOG(error) << Form("Inconsistent address in list of size %lu for TSRC %d%d%d%d: "
                                 "0x%08x, time  %f",
                                 fvLastHits[iSmType][iSm][iRpc][iCh].size(), iSmType, iSm, iRpc, iCh,
                                 fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetAddress(),
                                 fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetTime());
            fvLastHits[iSmType][iSm][iRpc][iCh].front()->Delete();
            fvLastHits[iSmType][iSm][iRpc][iCh].pop_front();
          }
        }  //fvLastHits[iSmType][iSm][iRpc][iCh].size()>1)
 
        // plot remaining time difference to previous hits
        if (fvLastHits[iSmType][iSm][iRpc][iCh].size() > 1) {  // check for previous hits in memory time interval
          CbmMatch* digiMatch = (CbmMatch*) fTofDigiMatchColl->At(iHitInd);
          Double_t dTotSum    = 0.;
          for (Int_t iLink = 0; iLink < digiMatch->GetNofLinks(); iLink += 2) {  // loop over digis
            CbmLink L0        = digiMatch->GetLink(iLink);
            Int_t iDigInd0    = L0.GetIndex();
            Int_t iDigInd1    = (digiMatch->GetLink(iLink + 1)).GetIndex();
            CbmTofDigi* pDig0 = (CbmTofDigi*) (&(fTofCalDigiVec->at(iDigInd0)));
            CbmTofDigi* pDig1 = (CbmTofDigi*) (&(fTofCalDigiVec->at(iDigInd1)));
            dTotSum += pDig0->GetTot() + pDig1->GetTot();
          }
 
          std::list<CbmTofHit*>::iterator itL = fvLastHits[iSmType][iSm][iRpc][iCh].end();
          itL--;
          for (UInt_t iH = 0; iH < fvLastHits[iSmType][iSm][iRpc][iCh].size() - 1; iH++) {
            itL--;
            fhRpcDTLastHits[iDetIndx]->Fill(TMath::Log10(pHit->GetTime() - (*itL)->GetTime()));
            fhRpcDTLastHits_CluSize[iDetIndx]->Fill(TMath::Log10(pHit->GetTime() - (*itL)->GetTime()),
                                                    digiMatch->GetNofLinks() / 2.);
            fhRpcDTLastHits_Tot[iDetIndx]->Fill(TMath::Log10(pHit->GetTime() - (*itL)->GetTime()), dTotSum);
          }
        }
      }  // iHitInd loop end
 
      // do reference first
      dTRef     = dDoubleMax;
      fTRefHits = 0;
      for (Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++) {
        pHit = (CbmTofHit*) fTofHitsColl->At(iHitInd);
        if (NULL == pHit) continue;
        Int_t iDetId = (pHit->GetAddress() & DetMask);
 
        if (fiBeamRefAddr == iDetId) {
          if (fviClusterMul[fiBeamRefType][fiBeamRefSm][fiBeamRefDet] > fiBeamRefMulMax) break;
          // Check Tot
          CbmMatch* digiMatch = (CbmMatch*) fTofDigiMatchColl->At(iHitInd);
          Double_t TotSum     = 0.;
          for (Int_t iLink = 0; iLink < digiMatch->GetNofLinks(); iLink += 2) {  // loop over digis
            CbmLink L0     = digiMatch->GetLink(iLink);                          //vDigish.at(ivDigInd);
            Int_t iDigInd0 = L0.GetIndex();
            if (iDigInd0 < (Int_t) fTofCalDigiVec->size()) {
              CbmTofDigi* pDig0 = &(fTofCalDigiVec->at(iDigInd0));
              TotSum += pDig0->GetTot();
            }
          }
          TotSum /= (0.5 * digiMatch->GetNofLinks());
          if (TotSum > fhRpcCluTot[iIndexDut]->GetYaxis()->GetXmax()) continue;  // ignore too large clusters
 
          fTRefHits = 1;
          if (pHit->GetTime() < dTRef) { dTRef = pHit->GetTime(); }
          iBeamRefMul++;
        }
        else {  //additional reference type multiplicity
          if (fiBeamRefType == CbmTofAddress::GetSmType(iDetId)) iBeamAddRefMul++;
        }
      }
      LOG(debug) << "FillHistos: BRefMul: " << iBeamRefMul << ", " << iBeamAddRefMul;
      if (iBeamRefMul == 0) return kFALSE;                  // don't fill histos without reference time
      if (iBeamAddRefMul < fiBeamAddRefMul) return kFALSE;  // ask for confirmation by other beam counters
 
      for (Int_t iSel = 0; iSel < iNSel; iSel++) {
        BSel[iSel]         = kFALSE;
        pTrig[iSel]        = NULL;
        Int_t iDutMul      = 0;
        Int_t iRefMul      = 0;
        Int_t iR0          = 0;
        Int_t iRl          = 0;
        ddXdZ[iSel]        = 0.;
        ddYdZ[iSel]        = 0.;
        dSel2dXdYMin[iSel] = 1.E300;
 
        switch (iSel) {
          case 0:  //  Detector under Test (Dut) && Diamonds,BeamRef
            iRl = fviClusterMul[fDutId][fDutSm].size();
            if (fDutRpc > -1) {
              iR0 = fDutRpc;
              iRl = fDutRpc + 1;
            }
            for (Int_t iRpc = iR0; iRpc < iRl; iRpc++)
              iDutMul += fviClusterMul[fDutId][fDutSm][iRpc];
            LOG(debug) << "Selector 0: DutMul " << fviClusterMul[fDutId][fDutSm][fDutRpc] << ", " << iDutMul
                       << ", BRefMul " << iBeamRefMul << " TRef: " << dTRef << ", BeamAddRefMul " << iBeamAddRefMul
                       << ", " << fiBeamAddRefMul;
            if (iDutMul > 0 && iDutMul < fiCluMulMax) {
              dTTrig[iSel] = dDoubleMax;
              // LOG(debug1)<<"Found selector 0, NbHits  "<<iNbTofHits;
              for (Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++) {
                pHit = (CbmTofHit*) fTofHitsColl->At(iHitInd);
                if (NULL == pHit) continue;
                Int_t iDetId = (pHit->GetAddress() & DetMask);
                // LOG(debug1)<<Form(" Det 0x%08x, Dut 0x%08x, T %f, TTrig %f",
                //              iDetId, fDutAddr, pHit->GetTime(),  dTTrig[iSel])
                //               ;
                //if( fDutId == CbmTofAddress::GetSmType( iDetId ))
                if (fDutAddr == iDetId) {
                  if (pHit->GetTime() < dTTrig[iSel]) {
                    dTTrig[iSel] = pHit->GetTime();
                    pTrig[iSel]  = pHit;
                    BSel[iSel]   = kTRUE;
                  }
                }
              }
              LOG(debug) << Form("Found selector 0 with mul %d from 0x%08x at %f ", iDutMul, pTrig[iSel]->GetAddress(),
                                 dTTrig[iSel]);
            }
            break;
 
          case 1:  // MRef & BRef
            iRl = fviClusterMul[fSelId][fSelSm].size();
            if (fSelRpc > -1) {
              iR0 = fSelRpc;
              iRl = fSelRpc + 1;
            }
            for (Int_t iRpc = iR0; iRpc < iRl; iRpc++)
              iRefMul += fviClusterMul[fSelId][fSelSm][iRpc];
            LOG(debug) << "FillHistos(): selector 1: RefMul " << fviClusterMul[fSelId][fSelSm][fSelRpc] << ", "
                       << iRefMul << ", BRefMul " << iBeamRefMul;
            if (iRefMul > 0 && iRefMul < fiCluMulMax) {
              // LOG(debug1)<<"FillHistos(): Found selector 1";
              dTTrig[iSel] = dDoubleMax;
              for (Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++) {
                pHit = (CbmTofHit*) fTofHitsColl->At(iHitInd);
                if (NULL == pHit) continue;
 
                Int_t iDetId = (pHit->GetAddress() & DetMask);
                if (fSelAddr == iDetId) {
                  if (pHit->GetTime() < dTTrig[iSel]) {
                    dTTrig[iSel] = pHit->GetTime();
                    pTrig[iSel]  = pHit;
                    BSel[iSel]   = kTRUE;
                  }
                }
              }
              LOG(debug) << Form("Found selector 1 with mul %d from 0x%08x at %f ", iRefMul, pTrig[iSel]->GetAddress(),
                                 dTTrig[iSel]);
            }
            break;
 
          default: LOG(info) << "FillHistos: selection not implemented " << iSel; ;
        }  // switch end
        if (fTRefMode > 10) { dTTrig[iSel] = dTRef; }
      }  // iSel - loop end
 
      if (fSel2Id > 0) {  // confirm selector by independent match
        for (Int_t iSel = 0; iSel < iNSel; iSel++) {
          if (BSel[iSel]) {
            BSel[iSel] = kFALSE;
            if (fviClusterMul[fSel2Id][fSel2Sm][fSel2Rpc] > 0
                && fviClusterMul[fSel2Id][fSel2Sm][fSel2Rpc] < fiCluMulMax)
              for (Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++) {
                pHit = (CbmTofHit*) fTofHitsColl->At(iHitInd);
                if (NULL == pHit) continue;
                Int_t iDetId = (pHit->GetAddress() & DetMask);
                if (fSel2Addr == iDetId) {
                  Double_t dzscal = 1.;
                  if (fEnableMatchPosScaling) dzscal = pHit->GetZ() / pTrig[iSel]->GetZ();
                  Double_t dSEl2dXdz = (pHit->GetX() - pTrig[iSel]->GetX()) / (pHit->GetZ() - pTrig[iSel]->GetZ());
                  Double_t dSEl2dYdz = (pHit->GetY() - pTrig[iSel]->GetY()) / (pHit->GetZ() - pTrig[iSel]->GetZ());
 
                  if (TMath::Sqrt(TMath::Power(pHit->GetX() - dzscal * pTrig[iSel]->GetX(), 2.)
                                  + TMath::Power(pHit->GetY() - dzscal * pTrig[iSel]->GetY(), 2.))
                      < fdCaldXdYMax) {
                    BSel[iSel]     = kTRUE;
                    Double_t dX2Y2 = TMath::Sqrt(dSEl2dXdz * dSEl2dXdz + dSEl2dYdz * dSEl2dYdz);
                    if (dX2Y2 < dSel2dXdYMin[iSel]) {
                      ddXdZ[iSel]        = dSEl2dXdz;
                      ddYdZ[iSel]        = dSEl2dYdz;
                      dSel2dXdYMin[iSel] = dX2Y2;
                    }
                    break;
                  }
                }
              }
          }  // BSel condition end
        }    // iSel lopp end
      }      // Sel2Id condition end
      UInt_t uTriggerPattern = 1;
 
      for (Int_t iSel = 0; iSel < iNSel; iSel++)
        if (BSel[iSel]) {
          uTriggerPattern |= (0x1 << (iSel * 3 + CbmTofAddress::GetRpcId(pTrig[iSel]->GetAddress() & DetMask)));
        }
 
      for (Int_t iSel = 0; iSel < iNSel; iSel++) {
        if (BSel[iSel]) {
          if (dTRef != 0. && fTRefHits > 0) {
            for (UInt_t uChannel = 0; uChannel < 16; uChannel++) {
              if (uTriggerPattern & (0x1 << uChannel)) { fhSeldT[iSel]->Fill(dTTrig[iSel] - dTRef, uChannel); }
            }
          }
        }
      }
    }  // 0<iNSel software triffer check end
 
    for (Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++) {
      pHit = (CbmTofHit*) fTofHitsColl->At(iHitInd);
      if (NULL == pHit) continue;
 
      Int_t iDetId                          = (pHit->GetAddress() & DetMask);
      std::map<UInt_t, UInt_t>::iterator it = fDetIdIndexMap.find(iDetId);
      if (it == fDetIdIndexMap.end()) continue;  // continue for invalid detector index
      Int_t iDetIndx = it->second;               //fDetIdIndexMap[iDetId];
 
      Int_t iSmType = CbmTofAddress::GetSmType(iDetId);
      Int_t iSm     = CbmTofAddress::GetSmId(iDetId);
      Int_t iRpc    = CbmTofAddress::GetRpcId(iDetId);
      Int_t iNbRpc  = fDigiBdfPar->GetNbRpc(iSmType);
      if (-1 < fviClusterMul[iSmType][iSm][iRpc]) {
        for (Int_t iSel = 0; iSel < iNSel; iSel++)
          if (BSel[iSel]) {
            Double_t w = fviClusterMul[iSmType][iSm][iRpc];
            if (w == 0.) w = 1.;
            else
              w = 1. / w;
            fhTRpcCluMul[iDetIndx][iSel]->Fill(fviClusterMul[iSmType][iSm][iRpc], w);
          }
      }
 
      if (fviClusterMul[iSmType][iSm][iRpc] > fiCluMulMax) continue;  // skip this event
      if (iBeamRefMul == 0) break;
 
      Int_t iChId  = pHit->GetAddress();
      fChannelInfo = fDigiPar->GetCell(iChId);
      Int_t iCh    = CbmTofAddress::GetChannelId(iChId);
      if (NULL == fChannelInfo) {
        LOG(error) << "Invalid Channel Pointer for ChId " << Form(" 0x%08x ", iChId) << ", Ch " << iCh;
        continue;
      }
      /*TGeoNode *fNode=*/  // prepare global->local trafo
      gGeoManager->FindNode(fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ());
      /*
        LOG(debug1)<<"Hit info: "
        <<Form(" 0x%08x %d %f %f %f %f %f %d",iChId,iCh,
        pHit->GetX(),pHit->GetY(),pHit->GetTime(),fChannelInfo->GetX(),fChannelInfo->GetY(), iHitInd )
        ;
      */
      Double_t hitpos[3];
      hitpos[0] = pHit->GetX();
      hitpos[1] = pHit->GetY();
      hitpos[2] = pHit->GetZ();
      Double_t hitpos_local[3];
      //TGeoNode* cNode=
      gGeoManager->GetCurrentNode();
      gGeoManager->MasterToLocal(hitpos, hitpos_local);
      /*
        LOG(debug1)<< Form(" MasterToLocal for %d, %d%d%d, node %p: (%6.1f,%6.1f,%6.1f) ->(%6.1f,%6.1f,%6.1f)",
                      iDetIndx, iSmType, iSm, iRpc,
                 cNode, hitpos[0], hitpos[1], hitpos[2],
                 hitpos_local[0], hitpos_local[1], hitpos_local[2])
              ;
      */
      fhRpcCluPosition[iDetIndx]->Fill((Double_t) iCh, hitpos_local[1]);  //pHit->GetY()-fChannelInfo->GetY());
      fhSmCluPosition[iSmType]->Fill((Double_t)(iSm * iNbRpc + iRpc), hitpos_local[1]);
 
      for (Int_t iSel = 0; iSel < iNSel; iSel++)
        if (BSel[iSel]) {
          fhTRpcCluPosition[iDetIndx][iSel]->Fill((Double_t) iCh,
                                                  hitpos_local[1]);  //pHit->GetY()-fChannelInfo->GetY());
          fhTSmCluPosition[iSmType][iSel]->Fill((Double_t)(iSm * iNbRpc + iRpc), hitpos_local[1]);
        }
 
      if (TMath::Abs(hitpos_local[1]) > fChannelInfo->GetSizey() * fPosYMaxScal) continue;
      /*
        LOG(debug1)<<" TofDigiMatchColl entries:"
                <<fTofDigiMatchColl->GetEntriesFast()
                ;
      */
      if (iHitInd > fTofDigiMatchColl->GetEntriesFast()) {
        LOG(error) << " Inconsistent DigiMatches for Hitind " << iHitInd
                   << ", TClonesArraySize: " << fTofDigiMatchColl->GetEntriesFast();
      }
      CbmMatch* digiMatch = (CbmMatch*) fTofDigiMatchColl->At(iHitInd);
      /*
        LOG(debug1)<<" got "
                <<digiMatch->GetNofLinks()<< " matches for iCh "<<iCh<<" at iHitInd "<<iHitInd
                ;
      */
      fhRpcCluSize[iDetIndx]->Fill((Double_t) iCh, digiMatch->GetNofLinks() / 2.);
 
      for (Int_t iSel = 0; iSel < iNSel; iSel++)
        if (BSel[iSel]) {
          fhTRpcCluSize[iDetIndx][iSel]->Fill((Double_t) iCh, digiMatch->GetNofLinks() / 2.);
          if (fvLastHits[iSmType][iSm][iRpc][iCh].size() > 1) {  // check for previous hits in memory time interval
            std::list<CbmTofHit*>::iterator itL = fvLastHits[iSmType][iSm][iRpc][iCh].end();
            itL--;
            for (UInt_t iH = 0; iH < fvLastHits[iSmType][iSm][iRpc][iCh].size() - 1; iH++) {
              itL--;
              fhTRpcCluSizeDTLastHits[iDetIndx][iSel]->Fill(TMath::Log10(pHit->GetTime() - (*itL)->GetTime()),
                                                            digiMatch->GetNofLinks() / 2.);
            }
          }
        }
 
      Double_t TotSum = 0.;
      for (Int_t iLink = 0; iLink < digiMatch->GetNofLinks(); iLink++) {  // loop over digis
        CbmLink L0     = digiMatch->GetLink(iLink);                       //vDigish.at(ivDigInd);
        Int_t iDigInd0 = L0.GetIndex();
        if (iDigInd0 < (Int_t) fTofCalDigiVec->size()) {
          CbmTofDigi* pDig0 = &(fTofCalDigiVec->at(iDigInd0));
          TotSum += pDig0->GetTot();
        }
      }
      Double_t dMeanTimeSquared = 0.;
      Double_t dNstrips         = 0.;
 
      Double_t dDelTof = 0.;
      Double_t dTcor[iNSel];
      Double_t dTTcor[iNSel];
      Double_t dZsign[iNSel];
      Double_t dzscal = 1.;
      //Double_t dDist=0.;
 
      for (Int_t iLink = 0; iLink < digiMatch->GetNofLinks(); iLink += 2) {  // loop over digis
        CbmLink L0     = digiMatch->GetLink(iLink);                          //vDigish.at(ivDigInd);
        Int_t iDigInd0 = L0.GetIndex();
        Int_t iDigInd1 = (digiMatch->GetLink(iLink + 1)).GetIndex();  //vDigish.at(ivDigInd+1);
        //LOG(debug1)<<" " << iDigInd0<<", "<<iDigInd1;
 
        if (iDigInd0 < (Int_t) fTofCalDigiVec->size() && iDigInd1 < (Int_t) fTofCalDigiVec->size()) {
          CbmTofDigi* pDig0 = &(fTofCalDigiVec->at(iDigInd0));
          CbmTofDigi* pDig1 = &(fTofCalDigiVec->at(iDigInd1));
          if ((Int_t) pDig0->GetType() != iSmType) {
            LOG(error) << Form(" Wrong Digi SmType for Tofhit %d in iDetIndx "
                               "%d, Ch %d with %3.0f strips at Indx %d, %d",
                               iHitInd, iDetIndx, iCh, dNstrips, iDigInd0, iDigInd1);
          }
          /*
            LOG(debug1)<<" fhRpcCluTot:  Digi 0 "<<iDigInd0<<": Ch "<<pDig0->GetChannel()<<", Side "<<pDig0->GetSide()
                       <<", StripSide "<<(Double_t)iCh*2.+pDig0->GetSide()
                       <<" Digi 1 "<<iDigInd1<<": Ch "<<pDig1->GetChannel()<<", Side "<<pDig1->GetSide()
                       <<" , StripSide "<<(Double_t)iCh*2.+pDig1->GetSide()
                       <<", Tot0 " << pDig0->GetTot() <<", Tot1 "<<pDig1->GetTot();
          */
          fhRpcCluTot[iDetIndx]->Fill(pDig0->GetChannel() * 2. + pDig0->GetSide(), pDig0->GetTot());
          fhRpcCluTot[iDetIndx]->Fill(pDig1->GetChannel() * 2. + pDig1->GetSide(), pDig1->GetTot());
 
          Int_t iCh0 = pDig0->GetChannel();
          Int_t iCh1 = pDig1->GetChannel();
          Int_t iS0  = pDig0->GetSide();
          Int_t iS1  = pDig1->GetSide();
          if (iCh0 != iCh1 || iS0 == iS1) {
            LOG(error) << Form(" MT2 for Tofhit %d in iDetIndx %d, Ch %d from %3.0f strips: ", iHitInd, iDetIndx, iCh,
                               dNstrips)
                       << Form(" Dig0: Ind %d, Ch %d, Side %d, T: %6.1f ", iDigInd0, iCh0, iS0, pDig0->GetTime())
                       << Form(" Dig1: Ind %d, Ch %d, Side %d, T: %6.1f ", iDigInd1, iCh1, iS1, pDig1->GetTime());
            continue;
          }
 
          if (0 > iCh0 || fDigiBdfPar->GetNbChan(iSmType, iRpc) <= iCh0) {
            LOG(error) << Form(" Wrong Digi for Tofhit %d in iDetIndx %d, Ch %d at Indx %d, %d "
                               "from %3.0f strips:  %d, %d, %d, %d",
                               iHitInd, iDetIndx, iCh, iDigInd0, iDigInd1, dNstrips, iCh0, iCh1, iS0, iS1);
            continue;
          }
 
          if (digiMatch->GetNofLinks() > 2)  //&& digiMatch->GetNofLinks()<8 ) // FIXME: hardcoded limits on CluSize
          {
            dNstrips += 1.;
            dMeanTimeSquared += TMath::Power(0.5 * (pDig0->GetTime() + pDig1->GetTime()) - pHit->GetTime(), 2);
            //             fhRpcCluAvWalk[iDetIndx]->Fill(0.5*(pDig0->GetTot()+pDig1->GetTot()),
            //                        0.5*(pDig0->GetTime()+pDig1->GetTime())-pHit->GetTime());
            fhRpcCluAvLnWalk[iDetIndx]->Fill(TMath::Log10(0.5 * (pDig0->GetTot() + pDig1->GetTot())),
                                             0.5 * (pDig0->GetTime() + pDig1->GetTime()) - pHit->GetTime());
 
            Double_t dTotWeigth = (pDig0->GetTot() + pDig1->GetTot()) / TotSum;
            Double_t dCorWeigth = 1. - dTotWeigth;
 
            fhRpcCluDelTOff[iDetIndx]->Fill(
              pDig0->GetChannel(), dCorWeigth * (0.5 * (pDig0->GetTime() + pDig1->GetTime()) - pHit->GetTime()));
 
            Double_t dDelPos = 0.5 * (pDig0->GetTime() - pDig1->GetTime()) * fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc);
            if (0 == pDig0->GetSide()) dDelPos *= -1.;
            fhRpcCluDelPos[iDetIndx]->Fill(pDig0->GetChannel(), dCorWeigth * (dDelPos - hitpos_local[1]));
 
            fhRpcCluWalk[iDetIndx][iCh0][iS0]->Fill(
              pDig0->GetTot(),
              pDig0->GetTime()
                - (pHit->GetTime()
                   - (1. - 2. * pDig0->GetSide()) * hitpos_local[1] / fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc)));
            fhRpcCluWalk[iDetIndx][iCh1][iS1]->Fill(
              pDig1->GetTot(),
              pDig1->GetTime()
                - (pHit->GetTime()
                   - (1. - 2. * pDig1->GetSide()) * hitpos_local[1] / fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc)));
 
            fhRpcCluAvWalk[iDetIndx]->Fill(pDig0->GetTot(), pDig0->GetTime()
                                                              - (pHit->GetTime()
                                                                 - (1. - 2. * pDig0->GetSide()) * hitpos_local[1]
                                                                     / fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc)));
            fhRpcCluAvWalk[iDetIndx]->Fill(pDig1->GetTot(), pDig1->GetTime()
                                                              - (pHit->GetTime()
                                                                 - (1. - 2. * pDig1->GetSide()) * hitpos_local[1]
                                                                     / fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc)));
          }  // end of Clustersize > 1 condition
          /*
            LOG(debug1)<<" fhTRpcCluTot: Digi 0 "<<iDigInd0<<": Ch "<<pDig0->GetChannel()<<", Side "<<pDig0->GetSide()
                    <<", StripSide "<<(Double_t)iCh*2.+pDig0->GetSide()
                    <<" Digi 1 "<<iDigInd1<<": Ch "<<pDig1->GetChannel()<<", Side "<<pDig1->GetSide()
                    <<", StripSide "<<(Double_t)iCh*2.+pDig1->GetSide()
                    ;
          */
          for (Int_t iSel = 0; iSel < iNSel; iSel++)
            if (BSel[iSel]) {
              if (NULL == pHit || NULL == pTrig[iSel]) {
                LOG(info) << " invalid pHit, iSel " << iSel << ", iDetIndx " << iDetIndx;
                break;
              }
              if (pHit->GetAddress() == pTrig[iSel]->GetAddress()) continue;
 
              fhTRpcCluTot[iDetIndx][iSel]->Fill(pDig0->GetChannel() * 2. + pDig0->GetSide(), pDig0->GetTot());
              fhTRpcCluTot[iDetIndx][iSel]->Fill(pDig1->GetChannel() * 2. + pDig1->GetSide(), pDig1->GetTot());
              if (fvLastHits[iSmType][iSm][iRpc][iCh].size() > 1) {  // check for previous hits in memory time interval
                std::list<CbmTofHit*>::iterator itL = fvLastHits[iSmType][iSm][iRpc][iCh].end();
                itL--;
                for (UInt_t iH = 0; iH < fvLastHits[iSmType][iSm][iRpc][iCh].size() - 1; iH++) {
                  itL--;
                  fhTRpcCluTotDTLastHits[iDetIndx][iSel]->Fill(TMath::Log10(pHit->GetTime() - (*itL)->GetTime()),
                                                               pDig0->GetTot());
                  fhTRpcCluTotDTLastHits[iDetIndx][iSel]->Fill(TMath::Log10(pHit->GetTime() - (*itL)->GetTime()),
                                                               pDig1->GetTot());
                }
              }
              if (iLink == 0) {  // Fill histo only once (for 1. digi entry)
                if (fEnableMatchPosScaling) dzscal = pHit->GetZ() / pTrig[iSel]->GetZ();
                fhTRpcCludXdY[iDetIndx][iSel]->Fill(pHit->GetX() - dzscal * pTrig[iSel]->GetX(),
                                                    pHit->GetY() - dzscal * pTrig[iSel]->GetY());
                dZsign[iSel] = 1.;
                if (pHit->GetZ() < pTrig[iSel]->GetZ()) dZsign[iSel] = -1.;
              }
              // look for geometrical match  with selector hit
              //if(  iSmType==fiBeamRefType      // to get entries in diamond/BeamRef histos
              if (iSmType == 5  // FIXME, to get entries in diamond histos
                  || TMath::Sqrt(TMath::Power(pHit->GetX() - dzscal * pTrig[iSel]->GetX(), 2.)
                                 + TMath::Power(pHit->GetY() - dzscal * pTrig[iSel]->GetY(), 2.))
                       < fdCaldXdYMax) {
                if (!fEnableMatchPosScaling && dSel2dXdYMin[iSel] < 1.E300)
                  if (TMath::Sqrt(
                        TMath::Power(pHit->GetX()
                                       - (pTrig[iSel]->GetX() + ddXdZ[iSel] * (pHit->GetZ() - (pTrig[iSel]->GetZ()))),
                                     2.)
                        + TMath::Power(pHit->GetY()
                                         - (pTrig[iSel]->GetY() + ddYdZ[iSel] * (pHit->GetZ() - (pTrig[iSel]->GetZ()))),
                                       2.))
                      > 0.5 * fdCaldXdYMax)
                    continue;      // refine position selection cut in cosmic measurement
                dTcor[iSel] = 0.;  // precaution
                if (dTRef != 0.
                    && TMath::Abs(dTRef - dTTrig[iSel]) < fdDelTofMax) {  // correct times for DelTof - velocity spread
                  if (iLink == 0) {  // do calculations only once (at 1. digi entry) // interpolate!
                    // calculate spatial distance to trigger hit
                    /*
                 dDist=TMath::Sqrt(TMath::Power(pHit->GetX()-pTrig[iSel]->GetX(),2.)
                                  +TMath::Power(pHit->GetY()-pTrig[iSel]->GetY(),2.)
                                  +TMath::Power(pHit->GetZ()-pTrig[iSel]->GetZ(),2.));
                 */
                    // determine correction value
                    //if(fiBeamRefAddr  != iDetId) // do not do this for reference counter itself
                    if (fTRefMode < 11)  // do not do this for trigger counter itself
                    {
                      Double_t dTentry = dTRef - dTTrig[iSel] + fdDelTofMax;
                      Int_t iBx        = dTentry / 2. / fdDelTofMax * nbClDelTofBinX;
                      if (iBx < 0) iBx = 0;
                      if (iBx > nbClDelTofBinX - 1) iBx = nbClDelTofBinX - 1;
                      Double_t dBinWidth = 2. * fdDelTofMax / nbClDelTofBinX;
                      Double_t dDTentry  = dTentry - ((Double_t) iBx) * dBinWidth;
                      Int_t iBx1         = 0;
                      dDTentry < 0 ? iBx1 = iBx - 1 : iBx1 = iBx + 1;
                      Double_t w0 = 1. - TMath::Abs(dDTentry) / dBinWidth;
                      Double_t w1 = 1. - w0;
                      if (iBx1 < 0) iBx1 = 0;
                      if (iBx1 > nbClDelTofBinX - 1) iBx1 = nbClDelTofBinX - 1;
                      dDelTof = fvCPDelTof[iSmType][iSm * iNbRpc + iRpc][iBx][iSel] * w0
                                + fvCPDelTof[iSmType][iSm * iNbRpc + iRpc][iBx1][iSel] * w1;
                      //dDelTof *= dDist; // has to be consistent with fhTRpcCluDelTof filling
                      /*
                        LOG(debug1)<<Form(" DelTof for SmT %d, Sm %d, R %d, T %d, dTRef %6.1f, Bx %d, Bx1 %d, DTe %6.1f -> DelT %6.1f",
                                     iSmType, iSm, iRpc, iSel, dTRef-dTTrig[iSel], iBx, iBx1, dDTentry, dDelTof)
                              ;
                      */
                    }
                    dTTcor[iSel] = dDelTof * dZsign[iSel];
                    dTcor[iSel]  = pHit->GetTime() - dDelTof - dTTrig[iSel];
                    //              Double_t dAvTot=0.5*(pDig0->GetTot()+pDig1->GetTot());
                  }  // if(iLink==0)
 
                  LOG(debug) << Form(" TRpcCluWalk for Ev %d, Link %d(%d), Sel %d, TSR %d%d%d, "
                                     "Ch %d,%d, S %d,%d T %f, DelTof %6.1f, W-ent:  %6.0f,%6.0f",
                                     fiNevtBuild, iLink, (Int_t) digiMatch->GetNofLinks(), iSel, iSmType, iSm, iRpc,
                                     iCh0, iCh1, iS0, iS1, dTTrig[iSel], dDelTof,
                                     fhTRpcCluWalk[iDetIndx][iSel][iCh0][iS0]->GetEntries(),
                                     fhTRpcCluWalk[iDetIndx][iSel][iCh1][iS1]->GetEntries());
 
                  if (fhTRpcCluWalk[iDetIndx][iSel][iCh0][iS0]->GetEntries()
                      != fhTRpcCluWalk[iDetIndx][iSel][iCh1][iS1]->GetEntries())
                    LOG(error) << Form(" Inconsistent walk histograms -> debugging "
                                       "necessary ... for %d, %d, %d, %d, %d, %d, %d ",
                                       fiNevtBuild, iDetIndx, iSel, iCh0, iCh1, iS0, iS1);
                  /*
                    LOG(debug1)<<Form(" TRpcCluWalk values side %d: %f, %f, side %d:  %f, %f ",
                                 iS0,pDig0->GetTot(),pDig0->GetTime()
                                 +((1.-2.*pDig0->GetSide())*hitpos_local[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc))-dTTcor[iSel]-dTTrig[iSel],
                                 iS1,pDig1->GetTot(),pDig1->GetTime()
                                 +((1.-2.*pDig1->GetSide())*hitpos_local[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc))-dTTcor[iSel]-dTTrig[iSel])
                          ;
                  */
                  fhTRpcCluWalk[iDetIndx][iSel][iCh0][iS0]->Fill(
                    pDig0->GetTot(),
                    //(pDig0->GetTime()+((1.-2.*pDig0->GetSide())*hitpos_local[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc))-dTTrig[iSel])-dTTcor[iSel]);
                    //dTcor[iSel]+(1.-2.*pDig0->GetSide())*hitpos_local[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc));
                    dTcor[iSel]);
                  fhTRpcCluWalk[iDetIndx][iSel][iCh1][iS1]->Fill(
                    pDig1->GetTot(),
                    //(pDig1->GetTime()+((1.-2.*pDig1->GetSide())*hitpos_local[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc))-dTTrig[iSel])-dTTcor[iSel]);
                    //dTcor[iSel]+(1.-2.*pDig1->GetSide())*hitpos_local[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc));
                    dTcor[iSel]);
 
                  fhTRpcCluAvWalk[iDetIndx][iSel]->Fill(
                    pDig0->GetTot(),
                    //(pDig0->GetTime()+((1.-2.*pDig0->GetSide())*hitpos_local[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc))-dTTrig[iSel])-dTTcor[iSel]);
                    //dTcor[iSel]+(1.-2.*pDig0->GetSide())*hitpos_local[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc));
                    dTcor[iSel]);
                  fhTRpcCluAvWalk[iDetIndx][iSel]->Fill(
                    pDig1->GetTot(),
                    //(pDig1->GetTime()+((1.-2.*pDig1->GetSide())*hitpos_local[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc))-dTTrig[iSel])-dTTcor[iSel]);
                    //dTcor[iSel]+(1.-2.*pDig1->GetSide())*hitpos_local[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc));
                    dTcor[iSel]);
 
                  if (iLink == 0) {  // Fill histo only once (for 1. digi entry)
                    //fhTRpcCluDelTof[iDetIndx][iSel]->Fill(dTRef-dTTrig[iSel],dTcor[iSel]/dDist);
                    fhTRpcCluDelTof[iDetIndx][iSel]->Fill(dTRef - dTTrig[iSel], dTcor[iSel]);
                    fhTSmCluTOff[iSmType][iSel]->Fill((Double_t)(iSm * iNbRpc + iRpc), dTcor[iSel]);
                    fhTSmCluTRun[iSmType][iSel]->Fill(fdEvent, dTcor[iSel]);
                    if (iDetId
                        != (pTrig[iSel]->GetAddress()
                            & DetMask)) {  // transform matched hit-pair back into detector frame
                      hitpos[0] = pHit->GetX() - dzscal * pTrig[iSel]->GetX() + fChannelInfo->GetX();
                      hitpos[1] = pHit->GetY() - dzscal * pTrig[iSel]->GetY() + fChannelInfo->GetY();
                      hitpos[2] = pHit->GetZ();
                      gGeoManager->MasterToLocal(hitpos, hitpos_local);  //  transform into local frame
                      fhRpcCluDelMatPos[iDetIndx]->Fill((Double_t) iCh, hitpos_local[1]);
                      fhRpcCluDelMatTOff[iDetIndx]->Fill((Double_t) iCh,
                                                         (pHit->GetTime() - dTTrig[iSel]) - dTTcor[iSel]);
                    }
                  }  // iLink==0 condition end
                }    // position condition end
              }      // Match condition end
            }        // closing of selector loop
        }
        else {
          LOG(error) << "FillHistos: invalid digi index " << iDetIndx << " digi0,1" << iDigInd0 << ", " << iDigInd1
                     << " - max:" << fTofCalDigiVec->size()
            //                       << " in event " << XXX
            ;
        }
      }  // iLink digi loop end;
      if (1 < dNstrips) {
        //           Double_t dVar=dMeanTimeSquared/dNstrips - TMath::Power(pHit->GetTime(),2);
        Double_t dVar = dMeanTimeSquared / (dNstrips - 1);
        //if(dVar<0.) dVar=0.;
        Double_t dTrms = TMath::Sqrt(dVar);
        LOG(debug) << Form(" Trms for Tofhit %d in iDetIndx %d, Ch %d from "
                           "%3.0f strips: %6.3f ns",
                           iHitInd, iDetIndx, iCh, dNstrips, dTrms);
        fhRpcCluTrms[iDetIndx]->Fill((Double_t) iCh, dTrms);
        pHit->SetTimeError(dTrms);
      }
      /*
        LOG(debug1)<<" Fill Time of iDetIndx "<<iDetIndx<<", hitAddr "
                <<Form(" %08x, y = %5.2f",pHit->GetAddress(),hitpos_local[1])
                <<" for |y| <"
                <<fhRpcCluPosition[iDetIndx]->GetYaxis()->GetXmax()
                ;
      */
      if (TMath::Abs(hitpos_local[1]) < (fhRpcCluPosition[iDetIndx]->GetYaxis()->GetXmax())) {
        if (dTRef != 0. && fTRefHits == 1) {
          fhRpcCluTOff[iDetIndx]->Fill((Double_t) iCh, pHit->GetTime() - dTRef);
          fhSmCluTOff[iSmType]->Fill((Double_t)(iSm * iNbRpc + iRpc), pHit->GetTime() - dTRef);
 
          for (Int_t iSel = 0; iSel < iNSel; iSel++)
            if (BSel[iSel]) {
              /*
                LOG(debug1)<<" TRpcCluTOff "<< iDetIndx <<", Sel "<< iSel
                           <<Form(", Dt %7.3f, LHsize %lu ",
                             pHit->GetTime()-dTTrig[iSel],fvLastHits[iSmType][iSm][iRpc][iCh].size());
             */
              if (pHit->GetAddress() == pTrig[iSel]->GetAddress()) continue;
 
              if (fvLastHits[iSmType][iSm][iRpc][iCh].size() > 1) {  // check for previous hits in memory time interval
                std::list<CbmTofHit*>::iterator itL = fvLastHits[iSmType][iSm][iRpc][iCh].end();
                itL--;
                for (UInt_t iH = 0; iH < fvLastHits[iSmType][iSm][iRpc][iCh].size() - 1; iH++) {
                  itL--;
                  //LOG(debug1)<<Form(" %f,",pHit->GetTime()-(*itL)->GetTime());
                }
              }
 
              // fill Time Offset histograms without velocity spread (DelTof) correction
              fhTRpcCluTOff[iDetIndx][iSel]->Fill((Double_t) iCh,
                                                  pHit->GetTime()
                                                    - dTTrig[iSel]);  // -dTTcor[iSel] only valid for matches
              if (fvLastHits[iSmType][iSm][iRpc][iCh].size() > 1) {   // check for previous hits in memory time interval
                std::list<CbmTofHit*>::iterator itL = fvLastHits[iSmType][iSm][iRpc][iCh].end();
                itL--;
                for (Int_t iH = 0; iH < 1; iH++) {  // use only last hit
                  //  for(Int_t iH=0; iH<fvLastHits[iSmType][iSm][iRpc][iCh].size()-1; iH++){//fill for all memorized hits
                  itL--;
                  Double_t dTsinceLast = pHit->GetTime() - (*itL)->GetTime();
                  if (dTsinceLast > fdMemoryTime)
                    LOG(error) << Form("Invalid Time since last hit on channel "
                                       "TSRC %d%d%d%d: %f > %f",
                                       iSmType, iSm, iRpc, iCh, dTsinceLast, fdMemoryTime);
 
                  fhTRpcCluTOffDTLastHits[iDetIndx][iSel]->Fill(TMath::Log10(dTsinceLast),
                                                                pHit->GetTime() - dTTrig[iSel]);
                  fhTRpcCluMemMulDTLastHits[iDetIndx][iSel]->Fill(TMath::Log10(dTsinceLast),
                                                                  fvLastHits[iSmType][iSm][iRpc][iCh].size() - 1);
                }
              }
            }
        }
      }
    }
  }
  return kTRUE;
}
 
static Int_t iNPulserFound = 0;
Bool_t CbmDeviceHitBuilderTof::MonitorPulser()
{
  iNPulserFound++;
  const Int_t iDet0   = fiPulDetRef;  // Define reference detector
  const Double_t Tlim = 0.5;
  Int_t iDigi0        = -1;
  switch (fiPulserMode) {
    case 1:                                             // mcbm :
      if ((UInt_t) fiNDigiIn != fiPulMulMin * 2 + 2) {  // 2 * working RPCs + 1 Diamond
        LOG(debug) << "Incomplete or distorted pulser event " << iNPulserFound << " with " << fiNDigiIn
                   << " digis instead of " << fiPulMulMin * 2 + 2;
        if ((UInt_t) fiNDigiIn < fiPulMulMin * 2) return kTRUE;
      }
      break;
    case 2:  // micro CBM cosmic
      if ((UInt_t) fiNDigiIn < fiPulMulMin * 2) return kTRUE;
      break;
      ;
    default:;
  }
 
  for (Int_t iDigi = 0; iDigi < (Int_t) fiNDigiIn; iDigi++) {
    Int_t iCh = fvDigiIn[iDigi].GetChannel();
    if (iCh != 0 && iCh != 31) continue;
 
    Int_t iAddr = fvDigiIn[iDigi].GetAddress() & DetMask;
    Int_t iDet  = fDetIdIndexMap[iAddr];
    if (iDet == iDet0) {
      Int_t iSide = fvDigiIn[iDigi].GetSide();
      if (iSide == 1) {  // for pulser mode 1
        iDigi0 = iDigi;
        break;
      }
    }
  }
 
  if (iDigi0 == -1) {
    LOG(debug) << Form("Pulser reference %d not found in pulser event %d of mul %d, return", iDet0, iNPulserFound,
                       fiNDigiIn);
    return kTRUE;
  }
 
  if (fvPulserTimes[iDet0][0].size() > 0)
    LOG(debug) << fiNDigiIn << " pulser digis with ref in "
               << Form("0x%08x at %d with deltaT %12.3f", fvDigiIn[iDigi0].GetAddress(), iDigi0,
                       fvDigiIn[iDigi0].GetTime() - fvPulserTimes[iDet0][0].back());
 
  for (Int_t iDigi = 0; iDigi < fiNDigiIn; iDigi++) {
    Int_t iCh = fvDigiIn[iDigi].GetChannel();
 
    Int_t iDetIndx = fDigiBdfPar->GetDetInd(fvDigiIn[iDigi].GetAddress() & DetMask);
    fhRpcDigiTot[iDetIndx]->Fill(2 * iCh + fvDigiIn[iDigi].GetSide(), fvDigiIn[iDigi].GetTot());
 
    Int_t iAddr = fvDigiIn[iDigi].GetAddress() & DetMask;
    Int_t iDet  = fDetIdIndexMap[iAddr];
    Int_t iSide = fvDigiIn[iDigi].GetSide();
 
    switch (fiPulserMode) {
      case 0:
      case 1:  // mCBM
        if (fvDigiIn[iDigi].GetType() != 5) {
          if (iCh != 0 && iCh != 31) continue;
          if (fvDigiIn[iDigi].GetTot() > fiPulTotMax || fvDigiIn[iDigi].GetTot() < fiPulTotMin) continue;
        }
        else {
          LOG(debug) << "Found PulHit of Type 5 in Ch " << iCh;
          //    if (iCh == 5)    iSide=1; // Special case for diamond, pulser for channels 5-8 in channel 5, stored as side 1
          //    if (iCh !=0 && iCh != 5) continue;
          if (iCh > 39) iSide = 1;  // Special case for diamond in Mar2019
          if (iCh != 0 && iCh != 40) continue;
        }
        break;
 
      case 2:
        if (fvDigiIn[iDigi].GetType() == 8)  // ceramic RPC
          if (fvDigiIn[iDigi].GetRpc() != 7) continue;
        if (iCh != 0 && iCh != 31) continue;
        break;
    }
 
    if (fvPulserTimes[iDet][iSide].size() == (UInt_t) NPulserTimes) fvPulserTimes[iDet][iSide].pop_front();
    if (iDet == iDet0 && 1 == iSide) {
      // check consistency of latest hit
      if (fvPulserTimes[iDet][iSide].size() > 1) {
        Double_t TDif = (fvPulserTimes[iDet][iSide].back() - fvPulserTimes[iDet][iSide].front())
                        / (fvPulserTimes[iDet][iSide].size() - 1);
        if (TMath::Abs(fvDigiIn[iDigi].GetTime() - fvPulserTimes[iDet][iSide].back() - TDif)
            > 1.E3) {  // probably due to parallel processing
 
          // action, reset
          fvPulserTimes[iDet][iSide].clear();
        }
        else {
          if (TMath::Abs(fvDigiIn[iDigi].GetTime() - fvPulserTimes[iDet][iSide].back() - TDif) > 10.) {
            LOG(info) << Form("Unexpected Pulser Time for event %d, pulser %d: "
                              "%15.1f instead %15.1f, TDif: %10.1f != %10.1f",
                              (int) fdEvent, iNPulserFound, fvDigiIn[iDigi].GetTime(),
                              fvPulserTimes[iDet][iSide].back() + TDif,
                              fvDigiIn[iDigi].GetTime() - fvPulserTimes[iDet][iSide].back(), TDif);
            // append dummy entry
            fvPulserTimes[iDet][iSide].push_back((Double_t)(fvPulserTimes[iDet][iSide].back() + TDif));
            continue;
          }
        }
      }
      // append new entry
      fvPulserTimes[iDet][iSide].push_back((Double_t)(fvDigiIn[iDigi].GetTime()));
    }
    else {
      Double_t dDeltaT0 = (Double_t)(fvDigiIn[iDigi].GetTime() - fvDigiIn[iDigi0].GetTime());
 
      // check consistency of latest hit
      if (fvPulserOffset[iDet][iSide] != 0.)
        if (TMath::Abs(dDeltaT0 - fvPulserOffset[iDet][iSide]) > Tlim) {
          LOG(info) << "Unexpected pulser offset at ev " << fdEvent << ", pulcnt " << iNPulserFound << " for Det "
                    << iDet << Form(", addr 0x%08x", iAddr) << ", side " << iSide << ": "
                    << dDeltaT0 - fvPulserOffset[iDet][iSide] << " ( " << fvPulserTimes[iDet][iSide].size() << " ) ";
          fvPulserTimes[iDet][iSide].clear();
          // skip this event
          // continue;
        }
 
      // fill monitoring histo
      if (fvPulserOffset[iDet][iSide] != 0.) {
        fhPulserTimesRaw->Fill(iDet * 2 + iSide, dDeltaT0 - fvPulserOffset[iDet][iSide]);
        fhPulserTimeRawEvo[iDet * 2 + iSide]->Fill(iNPulserFound, dDeltaT0 - fvPulserOffset[iDet][iSide]);
      }
      // append new entry
      fvPulserTimes[iDet][iSide].push_back(dDeltaT0);
    }
  }
 
  for (Int_t iDet = 0; iDet < (Int_t) fvPulserTimes.size(); iDet++) {
    for (Int_t iSide = 0; iSide < 2; iSide++) {
      if (iDet == iDet0 && iSide == 1) continue;  // skip reference counter
      if (fvPulserTimes[iDet][iSide].size() > 0) {
        Double_t Tmean = 0.;
        std::list<Double_t>::iterator it;
        for (it = fvPulserTimes[iDet][iSide].begin(); it != fvPulserTimes[iDet][iSide].end(); ++it)
          Tmean += *it;
        Tmean /= fvPulserTimes[iDet][iSide].size();
 
        if (TMath::Abs(Tmean - fvPulserOffset[iDet][iSide]) > Tlim)
          LOG(debug) << "New pulser offset at ev " << fdEvent << ", pulcnt " << iNPulserFound << " for Det " << iDet
                     << ", side " << iSide << ": " << Tmean << " ( " << fvPulserTimes[iDet][iSide].size() << " ) ";
 
        fvPulserOffset[iDet][iSide] = Tmean;
      }
    }
  }
 
  for (Int_t iDigi = 0; iDigi < fiNDigiIn; iDigi++) {
 
    Int_t iCh   = fvDigiIn[iDigi].GetChannel();
    Int_t iAddr = fvDigiIn[iDigi].GetAddress() & DetMask;
    Int_t iDet  = fDetIdIndexMap[iAddr];
    Int_t iSide = fvDigiIn[iDigi].GetSide();
 
    switch (fiPulserMode) {
      case 1:
        if (fvDigiIn[iDigi].GetType() != 5) {
          if (iCh != 0 && iCh != 31) continue;
        }
        else {
          //if (iCh == 5)    iSide=1; // Special case for diamond, pulser for channels 5-8 in channel 5, stored as side 1
          if (iCh > 39) iSide = 1;  // Special case for diamond in Mar2019
          //    if(iCh !=0 && iCh != 5) continue;
          if (iCh != 0 && iCh != 40) continue;
        }
        break;
      case 2:
        if (fvDigiIn[iDigi].GetType() == 8)  // ceramic RPC
          if (fvDigiIn[iDigi].GetRpc() != 7) continue;
        if (iCh != 0 && iCh != 31) continue;
        break;
    }
 
    Double_t dDeltaT0 =
      (Double_t)(fvDigiIn[iDigi].GetTime() - fvDigiIn[iDigi0].GetTime()) - fvPulserOffset[iDet][iSide];
 
    // fill monitoring histo
    fhPulserTimesCor->Fill(iDet * 2 + iSide, dDeltaT0);
  }
  return kTRUE;
}
 
Bool_t CbmDeviceHitBuilderTof::ApplyPulserCorrection()
{
  for (Int_t iDigi = 0; iDigi < fiNDigiIn; iDigi++) {
    Int_t iAddr = fvDigiIn[iDigi].GetAddress() & DetMask;
    Int_t iDet  = fDetIdIndexMap[iAddr];
    Int_t iSide = fvDigiIn[iDigi].GetSide();
    switch (fiPulserMode) {
      case 1:  // diamond has pulser in ch 0 and 5
        if (5 == fvDigiIn[iDigi].GetType()) {
          continue;  // no pulser correction for diamond, mapping to be resolved
          Int_t iCh = fvDigiIn[iDigi].GetChannel();
          if (iCh > 4) iSide = 1;
        }
        break;
      case 2:  // correct all cer pads by same rpc/side 0
        if (8 == fvDigiIn[iDigi].GetType() || 5 == fvDigiIn[iDigi].GetType()) {
          const Int_t iRefAddr = 0x00078006;
          iDet                 = fDetIdIndexMap[iRefAddr];
        }
        break;
    }
    fvDigiIn[iDigi].SetTime(fvDigiIn[iDigi].GetTime() - fvPulserOffset[iDet][iSide]);
  }
  return kTRUE;
}