CbmDeviceHitBuilderTof.cxx

/* Copyright (C) 2018-2021 PI-UHd, GSI
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Norbert Herrmann [committer] */

/**
 * CbmDeviceHitBuilderTof.cxx
 *
 * @since 2018-05-31
 * @author N. Herrmann
 */

#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;
    /*
    /// Save old global file and folder pointer to avoid messing with FairRoot
    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();
    /// Restore old global file and folder pointer to avoid messing with FairRoot
    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";

  /// Save old global file and folder pointer to avoid messing with FairRoot
  // <= 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();
  /// Restore old global file and folder pointer to avoid messing with FairRoot
  // <= 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);