// -----------------------------------------------------------------------------
// ----- pl_mc_inspection.C                                                -----
// -----                                                                   -----
// ----- created by C. Simon on 2018-04-24                                 -----
// -----                                                                   -----
// -----------------------------------------------------------------------------

//#include "drawStyle.C"

#include "tof/TofTools/CbmTofGeoHandler.h"
#include "CbmTofPoint.h"
#include "data/CbmMCTrack.h"
#include "FairLogger.h"

#include "TGeoManager.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TFile.h"
#include "TMath.h"
#include "TH2F.h"
#include "TClonesArray.h"
#include "TTree.h"
#include "TGeoPhysicalNode.h"
#include "TGaxis.h"
#include "TSystem.h"
#include "TSystemDirectory.h"
#include "TRegexp.h"
#include "TH1F.h"
#include "TF1.h"
#include "TMCProcess.h"
#include "TStyle.h"
#include "TProfile.h"
#include "TRandom.h"

#include "FairGeoLoader.h"
#include "FairGeoInterface.h"
#include "FairGeoBuilder.h"
#include "FairGeoMedia.h"
#include "FairGeoMedium.h"

#include "tof/TofMC/CbmTof.h"
#include "passive/CbmGeoCave.h"

#include <boost/regex.hpp>
#include <boost/lexical_cast.hpp>

#include <tuple>
#include <iostream>


struct TCounter
{
  Int_t fiModuleType;
  Int_t fiModuleIndex;
  Int_t fiCounterIndex;
  TGeoNode* fCounterNode;
  TGeoPhysicalNode* fCentralPN;
  TGeoPhysicalNode* fCounterBoxPN;
  Int_t fiNCells;
  Double_t fdCounterXWidth;
  Double_t fdCounterYWidth;
  Double_t fdCounterZWidth;
  Double_t fdCounterAcrossWidth;
  Double_t fdCounterAlongWidth;
  Double_t fdCellXWidth;
  Double_t fdCellYWidth;
  Double_t fdCellAcrossWidth;
  Double_t fdCellAlongWidth;
  Bool_t fbCellsAlongX;
  Bool_t fbCellIndexingAlongAxis;
  Int_t fiNGlassPlates;
  Int_t fiNGasGaps;
  Double_t fdGlassThickness;
  Double_t fdGapSize;
  Bool_t fbPadReadout;

  TCounter(Int_t iModuleType, Int_t iModuleIndex, Int_t iCounterIndex,
           TGeoNode* tCounterNode, TGeoPhysicalNode* tCentralPN, TGeoPhysicalNode* tCounterBoxPN) :
    fiModuleType(iModuleType),
    fiModuleIndex(iModuleIndex),
    fiCounterIndex(iCounterIndex),
    fCounterNode(tCounterNode),
    fCentralPN(tCentralPN),
    fCounterBoxPN(tCounterBoxPN),
    fiNCells(0),
    fdCounterXWidth(0.),
    fdCounterYWidth(0.),
    fdCounterZWidth(0.),
    fdCounterAcrossWidth(0.),
    fdCounterAlongWidth(0.),
    fdCellXWidth(0.),
    fdCellYWidth(0.),
    fdCellAcrossWidth(0.),
    fdCellAlongWidth(0.),
    fbCellsAlongX(kFALSE),
    fbCellIndexingAlongAxis(kFALSE),
    fiNGlassPlates(0),
    fiNGasGaps(0),
    fdGlassThickness(0.),
    fdGapSize(0.),
    fbPadReadout(kFALSE)
  {
    Double_t dAxisLow(0.);
    Double_t dAxisHigh(0.);

    for(Int_t iComponent = 0; iComponent < fCounterNode->GetNdaughters(); iComponent++)
    {
      if(TString(fCounterNode->GetDaughter(iComponent)->GetName()).Contains("glass", TString::kIgnoreCase))
      {
        TGeoNode* tPlateNode = fCounterNode->GetDaughter(iComponent);

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

        fdCounterZWidth += tPlateShape->GetAxisRange(3, dAxisLow, dAxisHigh);

        if(0 == fiNGlassPlates)
        {
          fdGlassThickness = tPlateShape->GetAxisRange(3, dAxisLow, dAxisHigh);
        }

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

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

        fdCounterZWidth += tGapShape->GetAxisRange(3, dAxisLow, dAxisHigh);

        if(0 == fiNGasGaps)
        {
          fdCounterXWidth = tGapShape->GetAxisRange(1, dAxisLow, dAxisHigh);
          fdCounterYWidth = tGapShape->GetAxisRange(2, dAxisLow, dAxisHigh);
          fdGapSize = tGapShape->GetAxisRange(3, dAxisLow, dAxisHigh);

          fiNCells = tGapNode->GetNdaughters();
          if(!fiNCells)
          {
            cout<<"MRPC gas gap is not divided into readout cells."<<endl;
          }

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

          fdCellXWidth = tFirstCellShape->GetAxisRange(1, dAxisLow, dAxisHigh);
          fdCellYWidth = tFirstCellShape->GetAxisRange(2, dAxisLow, dAxisHigh);


          if(fdCellXWidth == fdCounterXWidth)
          {
            fbCellsAlongX = kTRUE;

            fdCellAcrossWidth = fdCellYWidth;
            fdCellAlongWidth = fdCellXWidth;

            fdCounterAcrossWidth = fdCounterYWidth;
            fdCounterAlongWidth = fdCounterXWidth;
          }
          else
          {
            fdCellAcrossWidth = fdCellXWidth;
            fdCellAlongWidth = fdCellYWidth;

            fdCounterAcrossWidth = fdCounterXWidth;
            fdCounterAlongWidth = fdCounterYWidth;
          }

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

          TGeoNode* tLastCellNode = tGapNode->GetDaughter(fiNCells - 1);
          TGeoVolume* tLastCellVolume = tLastCellNode->GetVolume();
          TGeoShape* tLastCellShape = tLastCellVolume->GetShape();

          Double_t dCellMasterCoordinates[3] = {0.,0.,0.};
          Double_t dCellLocalCoordinates[3] = {0.,0.,0.};

          if(fbCellsAlongX)
          {
            Double_t dFirstCellLocalYLow(0.);
            Double_t dFirstCellLocalYHigh(0.);
            Double_t dFirstCellMasterYLow(0.);
            Double_t dFirstCellMasterYHigh(0.);

            tFirstCellShape->GetAxisRange(2, dFirstCellLocalYLow, dFirstCellLocalYHigh);

            dCellLocalCoordinates[1] = dFirstCellLocalYLow;
            tFirstCellNode->GetMatrix()->LocalToMaster(dCellLocalCoordinates, dCellMasterCoordinates);
            dFirstCellMasterYLow = dCellMasterCoordinates[1];

            dCellLocalCoordinates[1] = dFirstCellLocalYHigh;
            tFirstCellNode->GetMatrix()->LocalToMaster(dCellLocalCoordinates, dCellMasterCoordinates);
            dFirstCellMasterYHigh = dCellMasterCoordinates[1];

            Double_t dLastCellLocalYLow(0.);
            Double_t dLastCellLocalYHigh(0.);
            Double_t dLastCellMasterYLow(0.);
            Double_t dLastCellMasterYHigh(0.);

            tLastCellShape->GetAxisRange(2, dLastCellLocalYLow, dLastCellLocalYHigh);

            dCellLocalCoordinates[1] = dLastCellLocalYLow;
            tLastCellNode->GetMatrix()->LocalToMaster(dCellLocalCoordinates, dCellMasterCoordinates);
            dLastCellMasterYLow = dCellMasterCoordinates[1];

            dCellLocalCoordinates[1] = dLastCellLocalYHigh;
            tLastCellNode->GetMatrix()->LocalToMaster(dCellLocalCoordinates, dCellMasterCoordinates);
            dLastCellMasterYHigh = dCellMasterCoordinates[1];

            if(dFirstCellMasterYLow < 0. && dLastCellMasterYHigh > 0.)
            {
              fbCellIndexingAlongAxis = kTRUE;
            }
            else if(dLastCellMasterYLow < 0. && dFirstCellMasterYHigh > 0.)
            {
              fbCellIndexingAlongAxis = kFALSE;
            }
            else
            {
              cout<<"Incomplete understanding of geometry. Cannot determine direction of cell indexing."<<endl;
            }
          }
          else
          {
            Double_t dFirstCellLocalXLow(0.);
            Double_t dFirstCellLocalXHigh(0.);
            Double_t dFirstCellMasterXLow(0.);
            Double_t dFirstCellMasterXHigh(0.);

            tFirstCellShape->GetAxisRange(1, dFirstCellLocalXLow, dFirstCellLocalXHigh);

            dCellLocalCoordinates[0] = dFirstCellLocalXLow;
            tFirstCellNode->GetMatrix()->LocalToMaster(dCellLocalCoordinates, dCellMasterCoordinates);
            dFirstCellMasterXLow = dCellMasterCoordinates[0];

            dCellLocalCoordinates[0] = dFirstCellLocalXHigh;
            tFirstCellNode->GetMatrix()->LocalToMaster(dCellLocalCoordinates, dCellMasterCoordinates);
            dFirstCellMasterXHigh = dCellMasterCoordinates[0];

            Double_t dLastCellLocalXLow(0.);
            Double_t dLastCellLocalXHigh(0.);
            Double_t dLastCellMasterXLow(0.);
            Double_t dLastCellMasterXHigh(0.);

            tLastCellShape->GetAxisRange(1, dLastCellLocalXLow, dLastCellLocalXHigh);

            dCellLocalCoordinates[0] = dLastCellLocalXLow;
            tLastCellNode->GetMatrix()->LocalToMaster(dCellLocalCoordinates, dCellMasterCoordinates);
            dLastCellMasterXLow = dCellMasterCoordinates[0];

            dCellLocalCoordinates[0] = dLastCellLocalXHigh;
            tLastCellNode->GetMatrix()->LocalToMaster(dCellLocalCoordinates, dCellMasterCoordinates);
            dLastCellMasterXHigh = dCellMasterCoordinates[0];

            if(dFirstCellMasterXLow < 0. && dLastCellMasterXHigh > 0.)
            {
              fbCellIndexingAlongAxis = kTRUE;
            }
            else if(dLastCellMasterXLow < 0. && dFirstCellMasterXHigh > 0.)
            {
              fbCellIndexingAlongAxis = kFALSE;
            }
            else
            {
              cout<<"Incomplete understanding of geometry. Cannot determine direction of cell indexing."<<endl;
            }
          }
        }

        fiNGasGaps++;
      }
    }
  }
};

class TPointCluster
{
  public:
    TPointCluster() :
      fPointIndices(),
      fPointArray(),
      fdSummedPointTime(0.),
      fdSummedLocalPointX(0.),
      fdSummedLocalPointY(0.)
    {
    }

    Double_t GetPointTimeDifference(Int_t iInternalIndex) const
    {
      FairMCPoint* tPoint = dynamic_cast<FairMCPoint*>(fPointArray->At(fPointIndices.at(iInternalIndex)));

      return TMath::Abs(tPoint->GetTime() - GetClusterMeanTime());
    }

    Double_t GetPointPositionDifference(Int_t iInternalIndex) const
    {
      FairMCPoint* tPoint = dynamic_cast<FairMCPoint*>(fPointArray->At(fPointIndices.at(iInternalIndex)));

      Double_t dGlobalPoint[3] = {tPoint->GetX(), tPoint->GetY(), tPoint->GetZ()};
      Double_t dLocalPoint[3]  = {0., 0., 0.};

      fCounter->fCentralPN->GetMatrix()->MasterToLocal(dGlobalPoint, dLocalPoint);

      return TMath::Sqrt(TMath::Power(dLocalPoint[0] - GetClusterMeanPositionX(), 2.) + TMath::Power(dLocalPoint[1] - GetClusterMeanPositionY(), 2.));
    }

    Int_t GetKinshipDegree(Int_t iInternalIndex) const
    {
      FairMCPoint* tPoint = dynamic_cast<FairMCPoint*>(fPointArray->At(fPointIndices.at(iInternalIndex)));

      return tPoint->GetUniqueID();
    }

    Int_t GetNofPoints() const
    {
      return fPointIndices.size();
    }

    Double_t GetClusterMeanTime() const
    {
      return fdSummedPointTime/GetNofPoints();
    }

    Double_t GetClusterMeanPositionX() const
    {
      return fdSummedLocalPointX/GetNofPoints();
    }

    Double_t GetClusterMeanPositionY() const
    {
      return fdSummedLocalPointY/GetNofPoints();
    }

    void AddPoint(Int_t iPointIndex, Int_t iKinshipDegree)
    {
      fPointIndices.push_back(iPointIndex);

      FairMCPoint* tPoint = dynamic_cast<FairMCPoint*>(fPointArray->At(iPointIndex));
      tPoint->SetUniqueID(iKinshipDegree);

      fdSummedPointTime   += tPoint->GetTime();

      Double_t dGlobalPoint[3] = {tPoint->GetX(), tPoint->GetY(), tPoint->GetZ()};
      Double_t dLocalPoint[3]  = {0., 0., 0.};

      fCounter->fCentralPN->GetMatrix()->MasterToLocal(dGlobalPoint, dLocalPoint);

      fdSummedLocalPointX += dLocalPoint[0];
      fdSummedLocalPointY += dLocalPoint[1];
    }

    void SetPointArray(TClonesArray* tPointArray)
    {
      fPointArray = tPointArray;
    }

    void SetCounter(TCounter* tCounter)
    {
      fCounter = tCounter;
    }

  private:
   std::vector<Int_t> fPointIndices; // list of MC point indices
   TClonesArray* fPointArray;        // tree entry array of MC points
   TCounter* fCounter;               // counter containing the cluster
   Double_t fdSummedPointTime;       // summed MC point times
   Double_t fdSummedLocalPointX;     // summed MC point local X positions
   Double_t fdSummedLocalPointY;     // summed MC point local Y positions
};

struct cmp_str
{
   bool operator()(char const *a, char const *b)
   {
      return std::strcmp(a, b) < 0;
   }
};

std::map<Int_t, TCounter*> gCounterMap;
std::vector<std::tuple<Int_t, Int_t, Int_t>> gCounters;
CbmTofGeoHandler* gGeoHandler(NULL);
Int_t giCurrentModuleType(-1);
Int_t giCurrentModuleIndex(-1);
Int_t giCurrentCounterIndex(-1);
TString gCurrentNodePath("");
TString gTofNodePath("");
TString gCurrentBoxNodePath("");
std::map<Int_t, const char*> giPdgMap;
std::map<const char*, const char*, cmp_str> gcMaterialMap;
std::map<const char*, const char*, cmp_str> gcProcessMap;

Int_t CountPads(TVirtualPad *pad);
void AddCounter(Int_t iModuleType, Int_t iModuleIndex, Int_t iCounterIndex);
Bool_t FindTofNode(const char* cGeoFileName);
void ExpandNode(TGeoNode* tMotherNode);
void ReverseXAxis(TH1* h);
void GetPdgName(Int_t iPdgCode, const char*& cPdgName);
void GetMaterialName(const char* cMaterial, const char*& cMaterialName);
void GetProcessName(const char* cProcess, const char*& cProcessName);
Double_t fun_pty_tantheta(Double_t* x, Double_t* par);


void pl_mc_inspection(const char* cFileName, Int_t iNEvents, Double_t dInteractionProbability, Bool_t bBatchMode, const char* cGeoFileName)
{
  gStyle->SetOptStat("emr");
//  GeneralStyle();
  gStyle->SetPaintTextFormat(".2f");
  gErrorIgnoreLevel = kWarning;
  FairLogger::GetLogger()->SetLogScreenLevel("FATAL");
  FairLogger::GetLogger()->SetLogVerbosityLevel("LOW");

//  TLine l;

// feb15

  Double_t dBeamIntensity = 2.37e7;

//  AddCounter(2, 0, 0);
  AddCounter(3, 0, 0);
  AddCounter(9, 0, 0);
  AddCounter(7, 0, 0);
  AddCounter(4, 0, 0);
//  AddCounter(2, 1, 0);

// nov15
/*
//  Double_t dBeamIntensity = 6.7e6;
//  Double_t dBeamIntensity = 1.3e6;
//  Double_t dBeamIntensity = 3.2e6;
//  Double_t dBeamIntensity = 1.0e6;
  Double_t dBeamIntensity = 3.0e6;

//  AddCounter(0, 0, 0);
//  AddCounter(0, 1, 0);

  AddCounter(4, 0, 0);
//  AddCounter(9, 0, 0);
//  AddCounter(9, 0, 1);
  AddCounter(9, 1, 0);
//  AddCounter(9, 1, 1);
//  AddCounter(9, 2, 0);
  AddCounter(9, 2, 1);
  AddCounter(3, 0, 0);
*/

  // mean theta acceptance
  const Double_t kdTyMean = 66.8/367.4;

  // track multiplicity axis (Y1)
  const Double_t kY1Min(0.01);
//  const Double_t kY1Max(20.);
  const Double_t kY1Max(3.);
  const TString kY1Title("track multiplicity []");

  // track flux axis (Y2)
  const Bool_t bFluxAxis = kTRUE;
  const Double_t kY2Min(kY1Min*dBeamIntensity*dInteractionProbability/1000.);
  const Double_t kY2Max(kY1Max*dBeamIntensity*dInteractionProbability/1000.);
  const TString kY2Title(bFluxAxis ? "track flux [kHz/cm^{2}]" : "track rate [kHz]");


  // Gauge bosons
  giPdgMap.emplace( 22, "#gamma");  
  // leptons
  giPdgMap.emplace( 11, "e^{-}");
  giPdgMap.emplace(-11, "e^{+}");
  giPdgMap.emplace( 13, "#mu^{-}");
  giPdgMap.emplace(-13, "#mu^{+}");
  // mesons
  giPdgMap.emplace( 111, "#pi^{0}");
  giPdgMap.emplace( 211, "#pi^{+}");
  giPdgMap.emplace(-211, "#pi^{-}");
  giPdgMap.emplace( 221, "#eta");
  giPdgMap.emplace( 130, "#Kappa^{0}_{L}");
  giPdgMap.emplace( 310, "#Kappa^{0}_{S}");
  giPdgMap.emplace( 321, "#Kappa^{+}");
  giPdgMap.emplace(-321, "#Kappa^{-}");
  // baryons
  giPdgMap.emplace( 2112, "n");
  giPdgMap.emplace(-2112, "#bar{n}");
  giPdgMap.emplace( 2212, "p");
  giPdgMap.emplace(-2212, "#bar{p}");
  giPdgMap.emplace( 3112, "#Sigma^{-}");
  giPdgMap.emplace(-3112, "#Sigma^{+}");
  giPdgMap.emplace( 3122, "#Lambda");
  giPdgMap.emplace(-3122, "#bar{#Lambda}");
  giPdgMap.emplace( 3212, "#Sigma^{0}");
  giPdgMap.emplace(-3212, "#bar{#Sigma^{0}}");
  giPdgMap.emplace( 3222, "#Sigma^{+}");
  giPdgMap.emplace(-3222, "#Sigma^{-}");
  giPdgMap.emplace( 3312, "#Xi^{-}");
  giPdgMap.emplace(-3312, "#Xi^{+}");
  giPdgMap.emplace( 3322, "#Xi^{0}");
  giPdgMap.emplace(-3322, "#bar{#Xi^{0}}");
  giPdgMap.emplace( 3334, "#Omega^{-}");
  giPdgMap.emplace(-3334, "#Omega^{+}");
  // fragments
  giPdgMap.emplace( 1000010020, "d");
  giPdgMap.emplace( 1000010030, "t");
  giPdgMap.emplace( 1000020040, "#alpha");
  // dummy
  giPdgMap.emplace( 0, "Geant");


  gcMaterialMap.emplace("RPCgas", "gas");
  gcMaterialMap.emplace("RPCgas_noact", "gas");
  gcMaterialMap.emplace("RPCglass", "glass");
  gcMaterialMap.emplace("air", "air");
  gcMaterialMap.emplace("aluminium", "box");
  gcMaterialMap.emplace("carbon", "pcb");
  gcMaterialMap.emplace("targetMaterial", "target");
  gcMaterialMap.emplace("silicon", "diamond");
  gcMaterialMap.emplace("Scintillator", "crystal");
/*
  gcProcessMap.emplace("Compton scattering", "Compton");
  gcProcessMap.emplace("Decay", "Decay");
  gcProcessMap.emplace("Delta ray", "Delta");
  gcProcessMap.emplace("Hadronic interaction", "Hadronic");
  gcProcessMap.emplace("Lepton pair production", "Pair");
  gcProcessMap.emplace("Photoelectric effect", "Photo");
  gcProcessMap.emplace("Positron annihilation", "Annil");
  gcProcessMap.emplace("Primary particle emission", "Primary");
  gcProcessMap.emplace("Bremstrahlung", "Brems");
*/
  gcProcessMap.emplace("Compton scattering", "COMP");
  gcProcessMap.emplace("Decay", "DCAY");
  gcProcessMap.emplace("Delta ray", "DRAY");
  gcProcessMap.emplace("Hadronic interaction", "HADR");
  gcProcessMap.emplace("Lepton pair production", "PAIR");
  gcProcessMap.emplace("Photoelectric effect", "PHOT");
  gcProcessMap.emplace("Positron annihilation", "ANNI");
  gcProcessMap.emplace("Primary particle emission", "PRIM");
  gcProcessMap.emplace("Bremstrahlung", "BREM");

  std::map<Int_t, TH1*> tHistoMeanTrackMultiplicities;
  std::map<Int_t, TH1*> tHistoTriggerMeanTrackMultiplicities;
  std::map<Int_t, TH1*> tHistoRefTrackMeanTrackMultiplicities;
  std::map<Int_t, TH1*> tHistoMeanTrackQualityMultiplicities;
  std::map<Int_t, TH1*> tHistoClusterSize;
  std::map<Int_t, TH2*> tHistoEventClusterMulCell;
  std::map<Int_t, TH2*> tHistoTriggerClusterMulCell;
  std::map<Int_t, TH2*> tHistoRefTrackClusterMulCell;
  std::map<Int_t, TH1*> tHistoEventCellClusterOccupancy;
  std::map<Int_t, TH1*> tHistoTriggerCellClusterOccupancy;
  std::map<Int_t, TH1*> tHistoRefTrackCellClusterOccupancy;
  std::map<Int_t, TH1*> tHistoEventCellMultiClusterProbability;
  std::map<Int_t, TH1*> tHistoTriggerCellMultiClusterProbability;
  std::map<Int_t, TH1*> tHistoRefTrackCellMultiClusterProbability;
  std::map<Int_t, TH1*> tHistoClusterKinshipDegree;
  std::map<Int_t, TH1*> tHistoClusterInternalTimeDifference;
  std::map<Int_t, TH1*> tHistoClusterInternalPositionDifference;
  std::map<Int_t, TH1*> tHistoForeignChargedTimeDifference;
  std::map<Int_t, TH1*> tHistoForeignChargedPositionDifference;
  std::map<Int_t, TH1*> tHistoAnyForeignTrackMultiplicityDistribution;
  std::map<Int_t, TH1*> tHistoAnyForeignTrackCreationProcess;
  std::map<Int_t, TH1*> tHistoAnyForeignTrackCreationMaterial;
  std::map<Int_t, TH2*> tHistoAnyForeignTrackCreationProcessAndMaterial;
  std::map<Int_t, TH2*> tHistoAnyForeignTrackCreationProcessAndSpecies;
  std::map<Int_t, TH1*> tHistoAnyDomesticTrackCreationProcess;
  std::map<Int_t, TH1*> tHistoAnyDomesticTrackCreationMaterial;
  std::map<Int_t, TH2*> tHistoAnyDomesticTrackCreationProcessAndMaterial;
  std::map<Int_t, TH2*> tHistoAnyDomesticTrackCreationProcessAndSpecies;
  std::map<Int_t, TH1*> tHistoMeanLocalXClusterFlux;
  std::map<Int_t, TH1*> tHistoMeanLocalYClusterFlux;
  std::map<Int_t, TH2*> tHistoMeanLocalXYClusterFlux;
  std::map<Int_t, TH1*> tHistoMeanClusterTimeInEvent;
  std::map<Int_t, TH1*> tHistoAnyForeignTrackGlobalProductionZ;
  std::map<Int_t, TH1*> tHistoAnyPointLogTimeDistance;
  std::map<Int_t, TH1*> tHistoForeignChargedLogTimeDistance;

  TH1* tHistoTrackableTargetTrackGlobalProductionZ = NULL;
  TH1* tHistoTrackableRandomTrackGlobalProductionZ = NULL;
  TH1* tHistoFullTargetTrackMultiplicityDistribution = NULL;
  TH1* tHistoFullRandomTrackMultiplicityDistribution = NULL;
  TH2* tHistoTrackableTargetTrackAcceptance = NULL;
  TH2* tHistoTrackableTargetTrackCreationProcessAndSpecies = NULL;

  std::map<Int_t, TCanvas*> tCanvas;

  std::map<Int_t, Double_t> dPreviousPointTime;
  std::map<Int_t, Double_t> dPreviousForeignPointTime;
  std::map<Int_t, Double_t> dPreviousForeignPointLocalXPosition;
  std::map<Int_t, Double_t> dPreviousForeignPointLocalYPosition;


  if(!FindTofNode(cGeoFileName))
  {
    cout<<"Could not extract ToF information from geometry file."<<endl;
    return;
  }

  ExpandNode(gGeoManager->GetCurrentNode());

  gGeoManager->CdTop();

  gROOT->cd();

  TH1* tCurrentHistogram = NULL;

  for(auto const & itCounter : gCounterMap)
  {
    Int_t iCounterID = itCounter.first;
    TCounter* tCounter = itCounter.second;
    Int_t iModuleType = tCounter->fiModuleType;
    Int_t iModuleIndex = tCounter->fiModuleIndex;
    Int_t iCounterIndex = tCounter->fiCounterIndex;


    tCurrentHistogram = new TH1F(Form("tHistoMeanTrackMultiplicities_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 Form(" ; ; %s", kY1Title.Data()),
                                 5, 0., 5.);
    tCurrentHistogram->SetCanExtend(TH1::kXaxis);
    tCurrentHistogram->SetStats(0);

    tHistoMeanTrackMultiplicities.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));


    tCurrentHistogram = new TH1F(Form("tHistoTriggerMeanTrackMultiplicities_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 Form(" ; ; %s", kY1Title.Data()),
                                 5, 0., 5.);
    tCurrentHistogram->SetCanExtend(TH1::kXaxis);
    tCurrentHistogram->SetStats(0);

    tHistoTriggerMeanTrackMultiplicities.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));

    tCurrentHistogram = new TH1F(Form("tHistoRefTrackMeanTrackMultiplicities_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 Form(" ; ; %s", kY1Title.Data()),
                                 5, 0., 5.);
    tCurrentHistogram->SetCanExtend(TH1::kXaxis);
    tCurrentHistogram->SetStats(0);

    tHistoRefTrackMeanTrackMultiplicities.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));


    tCurrentHistogram = new TH1F(Form("tHistoMeanTrackQualityMultiplicities_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 Form(" ; ; %s", kY1Title.Data()),
                                 5, 0., 5.);
    tCurrentHistogram->SetCanExtend(TH1::kXaxis);
    tCurrentHistogram->SetStats(0);

    tHistoMeanTrackQualityMultiplicities.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));


    tCurrentHistogram = new TH1F(Form("tHistoClusterSize_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " MC track production; cluster size []; ",
                                 30, -0.5, 29.5);

    tHistoClusterSize.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));

    tCurrentHistogram = new TH2F(Form("tHistoEventClusterMulCell_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " ; readout cell []; MC cluster MUL []",
                                 tCounter->fiNCells, 0, tCounter->fiNCells, 100, -0.5, 99.5);

    tHistoEventClusterMulCell.insert(std::pair<Int_t, TH2*>(iCounterID, dynamic_cast<TH2*>(tCurrentHistogram)));

    tCurrentHistogram = new TH2F(Form("tHistoTriggerClusterMulCell_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " ; readout cell []; MC cluster MUL []",
                                 tCounter->fiNCells, 0, tCounter->fiNCells, 100, -0.5, 99.5);

    tHistoTriggerClusterMulCell.insert(std::pair<Int_t, TH2*>(iCounterID, dynamic_cast<TH2*>(tCurrentHistogram)));

    tCurrentHistogram = new TH2F(Form("tHistoRefTrackClusterMulCell_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " ; readout cell []; MC cluster MUL []",
                                 tCounter->fiNCells, 0, tCounter->fiNCells, 100, -0.5, 99.5);

    tHistoRefTrackClusterMulCell.insert(std::pair<Int_t, TH2*>(iCounterID, dynamic_cast<TH2*>(tCurrentHistogram)));

    tCurrentHistogram = new TH1F(Form("tHistoEventCellClusterOccupancy_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " ; readout cell []; track occupancy [%]",
                                 tCounter->fiNCells, 0, tCounter->fiNCells);

    tHistoEventCellClusterOccupancy.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));

    tCurrentHistogram = new TH1F(Form("tHistoTriggerCellClusterOccupancy_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " ; readout cell []; track occupancy [%]",
                                 tCounter->fiNCells, 0, tCounter->fiNCells);

    tHistoTriggerCellClusterOccupancy.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));

    tCurrentHistogram = new TH1F(Form("tHistoRefTrackCellClusterOccupancy_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " ; readout cell []; track occupancy [%]",
                                 tCounter->fiNCells, 0, tCounter->fiNCells);

    tHistoRefTrackCellClusterOccupancy.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));

    tCurrentHistogram = new TH1F(Form("tHistoEventCellMultiClusterProbability_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " ; readout cell []; multi-track prob. [%]",
                                 tCounter->fiNCells, 0, tCounter->fiNCells);

    tHistoEventCellMultiClusterProbability.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));

    tCurrentHistogram = new TH1F(Form("tHistoTriggerCellMultiClusterProbability_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " ; readout cell []; multi-track prob. [%]",
                                 tCounter->fiNCells, 0, tCounter->fiNCells);

    tHistoTriggerCellMultiClusterProbability.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));

    tCurrentHistogram = new TH1F(Form("tHistoRefTrackCellMultiClusterProbability_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " ; readout cell []; multi-track prob. [%]",
                                 tCounter->fiNCells, 0, tCounter->fiNCells);

    tHistoRefTrackCellMultiClusterProbability.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));

    tCurrentHistogram  = new TH1F(Form("tHistoClusterKinshipDegree_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                  " MC track production; kinship degree []; ",
                                  20, -0.5, 19.5);

    tHistoClusterKinshipDegree.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));


    tCurrentHistogram = new TH1F(Form("tHistoClusterInternalTimeDifference_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " internal MC cluster time; #Deltat to mean [ns]; ",
                                 500, 0., 5.);

    tHistoClusterInternalTimeDifference.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));


    tCurrentHistogram = new TH1F(Form("tHistoClusterInternalPositionDifference_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " internal MC cluster position; #Deltar to mean [cm]; ",
                                 20000, 0., 20.);

    tHistoClusterInternalPositionDifference.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));


    tCurrentHistogram = new TH1F(Form("tHistoForeignChargedTimeDifference_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " two subsequent foreign points; #Deltat [ns]; ",
                                 2000, 0., 20.);

    tHistoForeignChargedTimeDifference.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));


    tCurrentHistogram = new TH1F(Form("tHistoForeignChargedPositionDifference_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " two subsequent foreign points; #Deltar [cm]; ",
                                 500, 0., 50.);

    tHistoForeignChargedPositionDifference.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));


    tCurrentHistogram = new TH1F(Form("tHistoAnyForeignTrackMultiplicityDistribution_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " #color[2]{random}/#color[3]{target} good and all tracks; track multiplicity []; ",
                                 100, -0.5, 99.5);

    tHistoAnyForeignTrackMultiplicityDistribution.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));


    tCurrentHistogram = new TH1F(Form("tHistoAnyForeignTrackCreationProcess_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " foreign creation process; ; ",
                                 5, 0., 5.);
    tCurrentHistogram->SetCanExtend(TH1::kXaxis);
    tCurrentHistogram->SetStats(0);

    tHistoAnyForeignTrackCreationProcess.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));


    tCurrentHistogram = new TH1F(Form("tHistoAnyForeignTrackCreationMaterial_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " foreign creation material; ; ",
                                 5, 0., 5.);
    tCurrentHistogram->SetCanExtend(TH1::kXaxis);
    tCurrentHistogram->SetStats(0);

    tHistoAnyForeignTrackCreationMaterial.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));


    tCurrentHistogram = new TH2F(Form("tHistoAnyForeignTrackCreationProcessAndMaterial_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " foreign creation material vs. process; ; ",
                                 5, 0., 5., 5, 0., 5.);
    tCurrentHistogram->SetCanExtend(TH1::kAllAxes);
    tCurrentHistogram->SetStats(0);

    tHistoAnyForeignTrackCreationProcessAndMaterial.insert(std::pair<Int_t, TH2*>(iCounterID, dynamic_cast<TH2*>(tCurrentHistogram)));


    tCurrentHistogram = new TH2F(Form("tHistoAnyForeignTrackCreationProcessAndSpecies_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " foreign creation species vs. process; ; ",
                                 5, 0., 5., 5, 0., 5.);
    tCurrentHistogram->SetCanExtend(TH1::kAllAxes);
    tCurrentHistogram->SetStats(0);

    tHistoAnyForeignTrackCreationProcessAndSpecies.insert(std::pair<Int_t, TH2*>(iCounterID, dynamic_cast<TH2*>(tCurrentHistogram)));


    tCurrentHistogram = new TH1F(Form("tHistoAnyDomesticTrackCreationProcess_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " domestic creation process; ; ",
                                 5, 0., 5.);
    tCurrentHistogram->SetCanExtend(TH1::kXaxis);
    tCurrentHistogram->SetStats(0);

    tHistoAnyDomesticTrackCreationProcess.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));


    tCurrentHistogram = new TH1F(Form("tHistoAnyDomesticTrackCreationMaterial_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " domestic creation material; ; ",
                                 5, 0., 5.);
    tCurrentHistogram->SetCanExtend(TH1::kXaxis);
    tCurrentHistogram->SetStats(0);

    tHistoAnyDomesticTrackCreationMaterial.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));


    tCurrentHistogram = new TH2F(Form("tHistoAnyDomesticTrackCreationProcessAndMaterial_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " domestic creation material vs. process; ; ",
                                 5, 0., 5., 5, 0., 5.);
    tCurrentHistogram->SetCanExtend(TH1::kAllAxes);
    tCurrentHistogram->SetStats(0);

    tHistoAnyDomesticTrackCreationProcessAndMaterial.insert(std::pair<Int_t, TH2*>(iCounterID, dynamic_cast<TH2*>(tCurrentHistogram)));


    tCurrentHistogram = new TH2F(Form("tHistoAnyDomesticTrackCreationProcessAndSpecies_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " domestic creation species vs. process; ; ",
                                 5, 0., 5., 5, 0., 5.);
    tCurrentHistogram->SetCanExtend(TH1::kAllAxes);
    tCurrentHistogram->SetStats(0);

    tHistoAnyDomesticTrackCreationProcessAndSpecies.insert(std::pair<Int_t, TH2*>(iCounterID, dynamic_cast<TH2*>(tCurrentHistogram)));


    Double_t dXLow  = TMath::Floor(-tCounter->fdCounterXWidth/2.) - 1.;
    Double_t dXHigh = TMath::Ceil(tCounter->fdCounterXWidth/2.) + 1.;
    Int_t iNXBins = /*2*/static_cast<Int_t>(dXHigh - dXLow);

    Double_t dYLow  = TMath::Floor(-tCounter->fdCounterYWidth/2.) - 1.;
    Double_t dYHigh = TMath::Ceil(tCounter->fdCounterYWidth/2.) + 1.;
    Int_t iNYBins = /*2*/static_cast<Int_t>(dYHigh - dYLow);

    tCurrentHistogram = new TH1F(Form("tHistoMeanLocalXClusterFlux_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " MC track cluster flux; local x [cm]",
                                 iNXBins, dXLow, dXHigh);

    tHistoMeanLocalXClusterFlux.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));

    tCurrentHistogram = new TH1F(Form("tHistoMeanLocalYClusterFlux_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " MC track cluster flux; local y [cm]",
                                 iNYBins, dYLow, dYHigh);

    tHistoMeanLocalYClusterFlux.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));

    tCurrentHistogram = new TH2F(Form("tHistoMeanLocalXYClusterFlux_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " MC track cluster flux; local x [cm]; local y [cm]",
                                 iNXBins, dXLow, dXHigh, iNYBins, dYLow, dYHigh);

    tHistoMeanLocalXYClusterFlux.insert(std::pair<Int_t, TH2*>(iCounterID, dynamic_cast<TH2*>(tCurrentHistogram)));


    tCurrentHistogram = new TH1F(Form("tHistoMeanClusterTimeInEvent_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " MC track cluster; time [ns]; ",
                                 1000, 0., 100.);

    tHistoMeanClusterTimeInEvent.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));


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

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

    tCurrentHistogram = new TH1F(Form("tHistoAnyPointLogTimeDistance_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " any two subsequent points; any point #Deltat [ns]; ",
                                 iNHitDistanceBins, dHitDistanceBinEdges);

    tHistoAnyPointLogTimeDistance.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));


    tCurrentHistogram = new TH1F(Form("tHistoForeignChargedLogTimeDistance_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " two subsequent foreign points; #Deltat [ns]; ",
                                 iNHitDistanceBins, dHitDistanceBinEdges);

    tHistoForeignChargedLogTimeDistance.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));



    tCurrentHistogram = new TH1F(Form("tHistoAnyForeignTrackGlobalProductionZ_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                 " #color[2]{random}/#color[3]{target} good and all tracks; start vertex Z [cm]; ",
                                 10100, -10., 1000.);

    tHistoAnyForeignTrackGlobalProductionZ.insert(std::pair<Int_t, TH1*>(iCounterID, tCurrentHistogram));


    dPreviousPointTime.insert(std::pair<Int_t, Double_t>(iCounterID, -1000.));
    dPreviousForeignPointTime.insert(std::pair<Int_t, Double_t>(iCounterID, -1000.));
    dPreviousForeignPointLocalXPosition.insert(std::pair<Int_t, Double_t>(iCounterID, -1000.));
    dPreviousForeignPointLocalYPosition.insert(std::pair<Int_t, Double_t>(iCounterID, -1000.));

  }


  tHistoFullTargetTrackMultiplicityDistribution = new TH1F("tHistoFullTargetTrackMultiplicityDistribution",
                                                           " ; ; ",
                                                           100, -0.5, 99.5);


  tHistoFullRandomTrackMultiplicityDistribution = new TH1F("tHistoFullRandomTrackMultiplicityDistribution",
                                                           " ; ; ",
                                                           100, -0.5, 99.5);


  tHistoTrackableTargetTrackGlobalProductionZ = new TH1F("tHistoTrackableTargetTrackGlobalProductionZ",
                                                         " ; ; ",
                                                         10100, -10., 1000.);


  tHistoTrackableRandomTrackGlobalProductionZ = new TH1F("tHistoTrackableRandomTrackGlobalProductionZ",
                                                         " ; ; ",
                                                         10100, -10., 1000.);

  tHistoTrackableTargetTrackAcceptance = new TH2F("tHistoTrackableTargetTrackAcceptance",
                                                  "ref track acceptance; y_{lab} []; p_{t}/m [1/c]",
                                                  30, -1., 4., 30, 0., 2.5);

  tHistoTrackableTargetTrackCreationProcessAndSpecies = new TH2F("tHistoTrackableTargetTrackCreationProcessAndSpecies",
                                                                 " ref track creation species vs. process; creation process; particle species",
                                                                 5, 0., 5., 5, 0., 5.);
  tHistoTrackableTargetTrackCreationProcessAndSpecies->SetCanExtend(TH1::kAllAxes);
  tHistoTrackableTargetTrackCreationProcessAndSpecies->SetStats(0);


  TClonesArray* tMCTrackArray = new TClonesArray("CbmMCTrack");
  TClonesArray* tTofPointArray = new TClonesArray("CbmTofPoint");

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

  TGeoNode* tNode(NULL);
  TGeoMedium* tMedium(NULL);
  TGeoMaterial* tMaterial(NULL);

  const char* cPdgName;
  const char* cMaterialName;
  const char* cProcessName;

  cout<<endl;

  TFile::Open(cFileName, "read");
  if(!gFile)
  {
    cout<<Form("Could not open MC file %s.", cFileName)<<endl;
    return;
  }

  TTree* tTree = dynamic_cast<TTree*>(gFile->FindObjectAny("cbmsim"));
  if(!tTree)
  {
    cout<<"TTree 'cbmsim' could not be retrieved from MC file."<<endl;
    return;
  }

  TBranch* tMCTrackBranch = tTree->GetBranch("MCTrack");
  TBranch* tTofPointBranch = tTree->GetBranch("TofPoint");

  tMCTrackBranch->SetAddress(&tMCTrackArray);
  tTofPointBranch->SetAddress(&tTofPointArray);

  Int_t iNTreeEntries = tTree->GetEntries();

  if(iNEvents > iNTreeEntries)
  {
    cout<<Form("Tree contains only %d/%d events.", iNTreeEntries, iNEvents)<<endl;
    return;
  }


  Int_t iNTriggerEvents(0);
  Int_t iNRefTrackEvents(0);

  for(Int_t iEntry = 0; iEntry < iNEvents; iEntry++)
  {
    if(0 == (iEntry + 1) % 10000)
    {
      cout<<Form("   event %4d", iEntry + 1)<<endl;
    }

    std::map<Int_t, std::map<Int_t, TPointCluster>> tPointClusterMap;
    std::vector<Int_t> tTrackableTargetTrackVector;
    std::vector<Int_t> tTrackableRandomTrackVector;
    std::map<Int_t, Int_t> tForeignTrackMultiplicity;
    Int_t iTrackableTargetTrackMultiplicity = 0;
    Int_t iTrackableRandomTrackMultiplicity = 0;
    Bool_t bTriggerEvent = kFALSE;

    tMCTrackArray->Clear("C");
    tTofPointArray->Clear("C");

    tTree->GetEntry(iEntry);

    tTofPointArray->Sort();

    Int_t iTofPointArrayEntries = tTofPointArray->GetEntriesFast();
    Int_t iMCTrackArrayEntries = tMCTrackArray->GetEntriesFast();

    for(auto const & itCounter : gCounterMap)
    {
      Int_t iCounterID = itCounter.first;
      tPointClusterMap.emplace(iCounterID, std::map<Int_t, TPointCluster>());
      tForeignTrackMultiplicity.emplace(iCounterID, 0);

      dPreviousPointTime.at(iCounterID) = -1000.;
      dPreviousForeignPointTime.at(iCounterID) = -1000.;
      dPreviousForeignPointLocalXPosition.at(iCounterID) = -1000.;
      dPreviousForeignPointLocalYPosition.at(iCounterID) = -1000.;
    }

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

      CbmMCTrack* tTrack = NULL;

      // necessary to catch diamond points without any MC track reference
      if(-1 != tPoint->GetTrackID())
      {
        tTrack = dynamic_cast<CbmMCTrack*>(tMCTrackArray->At(tPoint->GetTrackID()));
      }

      Int_t iCounterID = tPoint->GetDetectorID();

      auto it = gCounterMap.find(iCounterID);
      if(it != gCounterMap.end())
      {
        TCounter* tCounter = it->second;

        dGlobalPoint[0] = tPoint->GetX();
        dGlobalPoint[1] = tPoint->GetY();
        dGlobalPoint[2] = tPoint->GetZ();

        tCounter->fCentralPN->GetMatrix()->MasterToLocal(dGlobalPoint, dLocalPoint);

        if(-1000. != dPreviousPointTime.at(iCounterID))
        {
          tHistoAnyPointLogTimeDistance.at(iCounterID)->Fill(tPoint->GetTime() - dPreviousPointTime.at(iCounterID));
        }
        dPreviousPointTime.at(iCounterID) = tPoint->GetTime();

        if(tTrack)
        {
          bTriggerEvent = kTRUE;

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

          tNode = gGeoManager->FindNode(dGlobalTrackStart[0], dGlobalTrackStart[1], dGlobalTrackStart[2]);
          tMedium = tNode->GetMedium();
          tMaterial = tMedium->GetMaterial();

          GetPdgName(tTrack->GetPdgCode(), cPdgName);
          GetMaterialName(tMaterial->GetName(), cMaterialName);
          GetProcessName(TMCProcessName[tTrack->GetGeantProcessId()], cProcessName);

          if(it == gCounterMap.begin())
          {
            if(0 == std::strcmp("target", cMaterialName))
            {
              tTrackableTargetTrackVector.push_back(tPoint->GetTrackID());
            }
            else
            {
              tTrackableRandomTrackVector.push_back(tPoint->GetTrackID());
            }
          }

          tCounter->fCounterBoxPN->GetMatrix()->MasterToLocal(dGlobalTrackStart, dLocalTrackStart);

          if(tCounter->fCounterBoxPN->GetVolume()->Contains(dLocalTrackStart))
          {
            tHistoAnyDomesticTrackCreationProcess.at(iCounterID)->Fill(cProcessName, 1.);
            tHistoAnyDomesticTrackCreationMaterial.at(iCounterID)->Fill(cMaterialName, 1.);
            tHistoAnyDomesticTrackCreationProcessAndMaterial.at(iCounterID)->Fill(cProcessName, cMaterialName, 1.);
            tHistoAnyDomesticTrackCreationProcessAndSpecies.at(iCounterID)->Fill(cProcessName, cPdgName, 1.);
          }
          else
          {
            if(-1000. != dPreviousForeignPointTime.at(iCounterID))
            {
              tHistoForeignChargedTimeDifference.at(iCounterID)->Fill(tPoint->GetTime() - dPreviousForeignPointTime.at(iCounterID));
              tHistoForeignChargedLogTimeDistance.at(iCounterID)->Fill(tPoint->GetTime() - dPreviousForeignPointTime.at(iCounterID));
            }
            dPreviousForeignPointTime.at(iCounterID) = tPoint->GetTime();

            if(-1000. != dPreviousForeignPointLocalXPosition.at(iCounterID))
            {
              Double_t dPositionDifference = TMath::Sqrt( TMath::Power(dLocalPoint[0] - dPreviousForeignPointLocalXPosition.at(iCounterID), 2.)
                                                         +TMath::Power(dLocalPoint[1] - dPreviousForeignPointLocalYPosition.at(iCounterID), 2.));
              tHistoForeignChargedPositionDifference.at(iCounterID)->Fill(dPositionDifference);
            }
            dPreviousForeignPointLocalXPosition.at(iCounterID) = dLocalPoint[0];
            dPreviousForeignPointLocalYPosition.at(iCounterID) = dLocalPoint[1];
          }


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

          while(-1 != iMotherID)
          {
            tMotherTrack = dynamic_cast<CbmMCTrack*>(tMCTrackArray->At(iMotherID));

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

            tCounter->fCounterBoxPN->GetMatrix()->MasterToLocal(dGlobalMotherTrackStart, dLocalMotherTrackStart);

            // lowest-ranked mother particle created outside the counter box volume
            if(!tCounter->fCounterBoxPN->GetVolume()->Contains(dLocalMotherTrackStart))
            {
              std::map<Int_t, TPointCluster>& tTempMap = tPointClusterMap.at(iCounterID);

              auto itTemp = tTempMap.find(iMotherID);
              if(itTemp == tTempMap.end())
              {
                TPointCluster& tTempCluster = tTempMap[iMotherID];
                tTempCluster.SetPointArray(tTofPointArray);
                tTempCluster.SetCounter(tCounter);

                // charged particle created inside the counter box volume
                if(tPoint->GetTrackID() != iMotherID)
                {
                  tHistoMeanTrackMultiplicities.at(iCounterID)->Fill("INT C", 1.);
                }

                tTempCluster.AddPoint(iPoint, iKinshipDegree);

                if(!tMotherTrack->GetCharge())
                {
                  tHistoMeanTrackMultiplicities.at(iCounterID)->Fill("EXT U", 1.);
                }
                else
                {
                  tHistoMeanTrackMultiplicities.at(iCounterID)->Fill("EXT C", 1.);
                }

                tNode = gGeoManager->FindNode(dGlobalMotherTrackStart[0], dGlobalMotherTrackStart[1], dGlobalMotherTrackStart[2]);
                tMedium = tNode->GetMedium();
                tMaterial = tMedium->GetMaterial();

                GetPdgName(tMotherTrack->GetPdgCode(), cPdgName);
                GetMaterialName(tMaterial->GetName(), cMaterialName);
                GetProcessName(TMCProcessName[tMotherTrack->GetGeantProcessId()], cProcessName);

                tHistoAnyForeignTrackCreationProcess.at(iCounterID)->Fill(cProcessName, 1.);
                tHistoAnyForeignTrackCreationMaterial.at(iCounterID)->Fill(cMaterialName, 1.);
                tHistoAnyForeignTrackCreationProcessAndMaterial.at(iCounterID)->Fill(cProcessName, cMaterialName, 1.);
                tHistoAnyForeignTrackCreationProcessAndSpecies.at(iCounterID)->Fill(cProcessName, cPdgName, 1.);

                tHistoAnyForeignTrackGlobalProductionZ.at(iCounterID)->Fill(tMotherTrack->GetStartZ());

                tHistoMeanTrackMultiplicities.at(iCounterID)->Fill("EXT T", 1.);

                tForeignTrackMultiplicity.at(iCounterID)++;
              }
              else
              {
                TPointCluster& tTempCluster = tTempMap.at(iMotherID);

                if(tCounter->fCounterBoxPN->GetVolume()->Contains(dLocalTrackStart))
                {
                  tTempCluster.AddPoint(iPoint, iKinshipDegree);

                  tHistoMeanTrackMultiplicities.at(iCounterID)->Fill("INT C", 1.);
                }
              }

              tHistoMeanTrackMultiplicities.at(iCounterID)->Fill("MEM P", 1.);
              break;
            }

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

      }
    }


    // scan for trackable particles (= points in all counters) originating from the target
    for(auto itIndex = tTrackableTargetTrackVector.begin(); itIndex != tTrackableTargetTrackVector.end(); itIndex++)
    {
      Int_t& iTrackableTargetTrackIndex = *itIndex;

      if(1 < gCounterMap.size())
      {
        for(auto itCounter = gCounterMap.begin(); itCounter != gCounterMap.end(); itCounter++)
        {
          if(itCounter != gCounterMap.begin())
          {
            Int_t iCounterID = itCounter->first;
            Bool_t bFoundTrack = kFALSE;

            std::map<Int_t, TPointCluster>& tTempMap = tPointClusterMap.at(iCounterID);

            for(auto const & itCluster : tTempMap)
            {
              Int_t iClusterTrackIndex = itCluster.first;

              if(iClusterTrackIndex == iTrackableTargetTrackIndex)
              {
                bFoundTrack = kTRUE;
              }
            }

            if(!bFoundTrack)
            {
              iTrackableTargetTrackIndex = -1;
              break;
            }

          }
        }
      }

      if(-1 != iTrackableTargetTrackIndex)
      {
        CbmMCTrack* tTargetTrack = dynamic_cast<CbmMCTrack*>(tMCTrackArray->At(iTrackableTargetTrackIndex));

        tHistoTrackableTargetTrackGlobalProductionZ->Fill(tTargetTrack->GetStartZ());

        if(tTargetTrack->GetMass())
        {
          tHistoTrackableTargetTrackAcceptance->Fill(tTargetTrack->GetRapidity(), tTargetTrack->GetPt()/tTargetTrack->GetMass());
        }

        GetPdgName(tTargetTrack->GetPdgCode(), cPdgName);
        GetProcessName(TMCProcessName[tTargetTrack->GetGeantProcessId()], cProcessName);
        tHistoTrackableTargetTrackCreationProcessAndSpecies->Fill(cProcessName, cPdgName, 1.);

        for(auto const & itCounter : gCounterMap)
        {
          Int_t iCounterID = itCounter.first;

          tHistoMeanTrackQualityMultiplicities.at(iCounterID)->Fill("#color[3]{TRGT C}", 1.);
          tHistoMeanTrackQualityMultiplicities.at(iCounterID)->Fill("TRCK C", 1.);
        }

        iTrackableTargetTrackMultiplicity++;
      }
    }

    // scan for random trackable particles (= points in all counters) NOT originating from the target
    for(auto itIndex = tTrackableRandomTrackVector.begin(); itIndex != tTrackableRandomTrackVector.end(); itIndex++)
    {
      Int_t& iTrackableRandomTrackIndex = *itIndex;

      if(1 < gCounterMap.size())
      {
        for(auto itCounter = gCounterMap.begin(); itCounter != gCounterMap.end(); itCounter++)
        {
          if(itCounter != gCounterMap.begin())
          {
            Int_t iCounterID = itCounter->first;
            Bool_t bFoundTrack = kFALSE;

            std::map<Int_t, TPointCluster>& tTempMap = tPointClusterMap.at(iCounterID);

            for(auto const & itCluster : tTempMap)
            {
              Int_t iClusterTrackIndex = itCluster.first;

              if(iClusterTrackIndex == iTrackableRandomTrackIndex)
              {
                bFoundTrack = kTRUE;
              }
            }

            if(!bFoundTrack)
            {
              iTrackableRandomTrackIndex = -1;
              break;
            }

          }
        }
      }

      if(-1 != iTrackableRandomTrackIndex)
      {
        CbmMCTrack* tRandomTrackTrack = dynamic_cast<CbmMCTrack*>(tMCTrackArray->At(iTrackableRandomTrackIndex));

        tHistoTrackableRandomTrackGlobalProductionZ->Fill(tRandomTrackTrack->GetStartZ());

        for(auto const & itCounter : gCounterMap)
        {
          Int_t iCounterID = itCounter.first;

          tHistoMeanTrackQualityMultiplicities.at(iCounterID)->Fill("#color[2]{RNDM C}", 1.);
          tHistoMeanTrackQualityMultiplicities.at(iCounterID)->Fill("TRCK C", 1.);
        }

        iTrackableRandomTrackMultiplicity++;
      }
    }


    for(auto const & itCounter : gCounterMap)
    {
      Int_t iCounterID = itCounter.first;
      TCounter* tCounter = itCounter.second;

      std::map<Int_t, TPointCluster>& tTempMap = tPointClusterMap.at(iCounterID);

      std::map<Int_t, Int_t> ClusterMulCell;

      for(auto const & itCluster : tTempMap)
      {
        const TPointCluster& tCluster = itCluster.second;
        Int_t iNClusterPoints = tCluster.GetNofPoints();
        Int_t iClusterSize = iNClusterPoints;

        for(Int_t iClusterPoint = 0; iClusterPoint < iNClusterPoints; iClusterPoint++)
        {
          tHistoClusterKinshipDegree.at(iCounterID)->Fill(tCluster.GetKinshipDegree(iClusterPoint));

          if(1 < iNClusterPoints)
          {
            tHistoClusterInternalTimeDifference.at(iCounterID)->Fill(tCluster.GetPointTimeDifference(iClusterPoint));
            tHistoClusterInternalPositionDifference.at(iCounterID)->Fill(tCluster.GetPointPositionDifference(iClusterPoint));
          }

          if(0 == tCluster.GetKinshipDegree(iClusterPoint))
          {
            iClusterSize--;
          }
        }

        tHistoClusterSize.at(iCounterID)->Fill(iClusterSize);
        tHistoMeanLocalXClusterFlux.at(iCounterID)->Fill(tCluster.GetClusterMeanPositionX());
        tHistoMeanLocalYClusterFlux.at(iCounterID)->Fill(tCluster.GetClusterMeanPositionY());
        tHistoMeanLocalXYClusterFlux.at(iCounterID)->Fill(tCluster.GetClusterMeanPositionX(), tCluster.GetClusterMeanPositionY());
        tHistoMeanClusterTimeInEvent.at(iCounterID)->Fill(tCluster.GetClusterMeanTime());

        Double_t dLocalClusterPositionAcrossCells(0.);
        if(tCounter->fbCellsAlongX)
        {
          dLocalClusterPositionAcrossCells = tCluster.GetClusterMeanPositionY();
        }   
        else
        {
          dLocalClusterPositionAcrossCells = tCluster.GetClusterMeanPositionX();
        }

        Int_t iClusterCellIndex(-1);
        if(tCounter->fbCellIndexingAlongAxis)
        {
          iClusterCellIndex = static_cast<Int_t>(dLocalClusterPositionAcrossCells/tCounter->fdCellAcrossWidth + static_cast<Double_t>(tCounter->fiNCells)/2.);
        }   
        else
        {
          iClusterCellIndex = -static_cast<Int_t>(dLocalClusterPositionAcrossCells/tCounter->fdCellAcrossWidth - static_cast<Double_t>(tCounter->fiNCells)/2.);
        }

        ClusterMulCell[iClusterCellIndex]++;
      }

      for(Int_t iCell = 0; iCell < tCounter->fiNCells; iCell++)
      {
        tHistoEventClusterMulCell.at(iCounterID)->Fill(iCell, ClusterMulCell[iCell]);

        if(bTriggerEvent)
        {
          tHistoTriggerClusterMulCell.at(iCounterID)->Fill(iCell, ClusterMulCell[iCell]);
        }

        if(iTrackableTargetTrackMultiplicity > 0)
        {
          tHistoRefTrackClusterMulCell.at(iCounterID)->Fill(iCell, ClusterMulCell[iCell]);
        }
      }
    }


    for(auto const & itCounter : gCounterMap)
    {
      Int_t iCounterID = itCounter.first;

      tHistoMeanTrackQualityMultiplicities.at(iCounterID)->Fill("EXT T", tForeignTrackMultiplicity.at(iCounterID));
      tHistoAnyForeignTrackMultiplicityDistribution.at(iCounterID)->Fill(tForeignTrackMultiplicity.at(iCounterID));
    }

    tHistoFullTargetTrackMultiplicityDistribution->Fill(iTrackableTargetTrackMultiplicity);
    tHistoFullRandomTrackMultiplicityDistribution->Fill(iTrackableRandomTrackMultiplicity);


    // repeat point loop to get selected multiplicities
    tPointClusterMap.clear();

    for(auto const & itCounter : gCounterMap)
    {
      Int_t iCounterID = itCounter.first;
      tPointClusterMap.emplace(iCounterID, std::map<Int_t, TPointCluster>());
    }

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

      CbmMCTrack* tTrack = NULL;

      // necessary to catch diamond points without any MC track reference
      if(-1 != tPoint->GetTrackID())
      {
        tTrack = dynamic_cast<CbmMCTrack*>(tMCTrackArray->At(tPoint->GetTrackID()));
      }

      Int_t iCounterID = tPoint->GetDetectorID();

      auto it = gCounterMap.find(iCounterID);
      if(it != gCounterMap.end())
      {
        TCounter* tCounter = it->second;

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

          tCounter->fCounterBoxPN->GetMatrix()->MasterToLocal(dGlobalTrackStart, dLocalTrackStart);

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

          while(-1 != iMotherID)
          {
            tMotherTrack = dynamic_cast<CbmMCTrack*>(tMCTrackArray->At(iMotherID));

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

            tCounter->fCounterBoxPN->GetMatrix()->MasterToLocal(dGlobalMotherTrackStart, dLocalMotherTrackStart);

            // lowest-ranked mother particle created outside the counter box volume
            if(!tCounter->fCounterBoxPN->GetVolume()->Contains(dLocalMotherTrackStart))
            {
              std::map<Int_t, TPointCluster>& tTempMap = tPointClusterMap.at(iCounterID);

              auto itTemp = tTempMap.find(iMotherID);
              if(itTemp == tTempMap.end())
              {
                TPointCluster& tTempCluster = tTempMap[iMotherID];
                tTempCluster.SetPointArray(tTofPointArray);
                tTempCluster.SetCounter(tCounter);

                // charged particle created inside the counter box volume
                if(tPoint->GetTrackID() != iMotherID)
                {
                  if(bTriggerEvent)
                  {
                    tHistoTriggerMeanTrackMultiplicities.at(iCounterID)->Fill("INT C", 1.);
                  }

                  if(iTrackableTargetTrackMultiplicity > 0)
                  {
                    tHistoRefTrackMeanTrackMultiplicities.at(iCounterID)->Fill("INT C", 1.);
                  }
                }

                tTempCluster.AddPoint(iPoint, iKinshipDegree);

                if(!tMotherTrack->GetCharge())
                {
                  if(bTriggerEvent)
                  {
                    tHistoTriggerMeanTrackMultiplicities.at(iCounterID)->Fill("EXT U", 1.);
                  }

                  if(iTrackableTargetTrackMultiplicity > 0)
                  {
                    tHistoRefTrackMeanTrackMultiplicities.at(iCounterID)->Fill("EXT U", 1.);
                  }
                }
                else
                {
                  if(bTriggerEvent)
                  {
                    tHistoTriggerMeanTrackMultiplicities.at(iCounterID)->Fill("EXT C", 1.);
                  }

                  if(iTrackableTargetTrackMultiplicity > 0)
                  {
                    tHistoRefTrackMeanTrackMultiplicities.at(iCounterID)->Fill("EXT C", 1.);
                  }
                }

                if(bTriggerEvent)
                {
                  tHistoTriggerMeanTrackMultiplicities.at(iCounterID)->Fill("EXT T", 1.);
                }

                if(iTrackableTargetTrackMultiplicity > 0)
                {
                  tHistoRefTrackMeanTrackMultiplicities.at(iCounterID)->Fill("EXT T", 1.);
                }
              }
              else
              {
                TPointCluster& tTempCluster = tTempMap.at(iMotherID);

                if(tCounter->fCounterBoxPN->GetVolume()->Contains(dLocalTrackStart))
                {
                  tTempCluster.AddPoint(iPoint, iKinshipDegree);

                  if(bTriggerEvent)
                  {
                    tHistoTriggerMeanTrackMultiplicities.at(iCounterID)->Fill("INT C", 1.);
                  }

                  if(iTrackableTargetTrackMultiplicity > 0)
                  {
                    tHistoRefTrackMeanTrackMultiplicities.at(iCounterID)->Fill("INT C", 1.);
                  }
                }
              }

              if(bTriggerEvent)
              {
                tHistoTriggerMeanTrackMultiplicities.at(iCounterID)->Fill("MEM P", 1.);
              }

              if(iTrackableTargetTrackMultiplicity > 0)
              {
                tHistoRefTrackMeanTrackMultiplicities.at(iCounterID)->Fill("MEM P", 1.);
              }

              break;
            }

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

      }
    }


    if(bTriggerEvent)
    {
      iNTriggerEvents++;
    }

    if(iTrackableTargetTrackMultiplicity > 0)
    {
      iNRefTrackEvents++;
    }

  }


  tMCTrackArray->Clear("C");
  tTofPointArray->Clear("C");

  gFile->Close();
  delete gFile;

  delete tMCTrackArray;
  delete tTofPointArray;


  tHistoTrackableTargetTrackGlobalProductionZ->SetLineColor(kGreen);
  tHistoTrackableRandomTrackGlobalProductionZ->SetLineColor(kRed);

  tHistoFullTargetTrackMultiplicityDistribution->SetFillColor(kGreen);
  tHistoFullRandomTrackMultiplicityDistribution->SetFillColor(kRed);

  for(auto const & itCounter : gCounterMap)
  {
    Int_t iCounterID = itCounter.first;
    TCounter* tCounter = itCounter.second;
    Int_t iModuleType = tCounter->fiModuleType;
    Int_t iModuleIndex = tCounter->fiModuleIndex;
    Int_t iCounterIndex = tCounter->fiCounterIndex;

    Double_t dCounterArea = (tCounter->fdCounterXWidth)*(tCounter->fdCounterYWidth);

    TCanvas* tCurrentCanvas = new TCanvas(Form("tCanvas_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex),
                                          Form("MC inspection of counter %d-%d-%d", iModuleType, iModuleIndex, iCounterIndex),
                                          0, 0, 1750, 1400);

    tCanvas.insert(std::pair<Int_t, TCanvas*>(iCounterID, tCurrentCanvas));

    tCurrentCanvas->Divide(5, 4);


    tHistoMeanTrackMultiplicities.at(iCounterID)->Scale(1./iNEvents);
    tHistoTriggerMeanTrackMultiplicities.at(iCounterID)->Scale(1./iNTriggerEvents);
    tHistoRefTrackMeanTrackMultiplicities.at(iCounterID)->Scale(1./iNRefTrackEvents);
    tHistoMeanTrackQualityMultiplicities.at(iCounterID)->Scale(1./iNEvents);
    tHistoMeanLocalXClusterFlux.at(iCounterID)->Scale(/*2*/dBeamIntensity*dInteractionProbability/iNEvents/tCounter->fdCounterYWidth);
    tHistoMeanLocalYClusterFlux.at(iCounterID)->Scale(/*2*/dBeamIntensity*dInteractionProbability/iNEvents/tCounter->fdCounterXWidth);
    tHistoMeanLocalXYClusterFlux.at(iCounterID)->Scale(/*4*/dBeamIntensity*dInteractionProbability/iNEvents);

    TH1D* tHisto;
    Double_t dNTotalEvents;
    Double_t dNOccupiedEvents;
    Double_t dNMultiEvents;

    for(Int_t iCell = 0; iCell < tCounter->fiNCells; iCell++)
    {
      tHisto = tHistoEventClusterMulCell.at(iCounterID)->ProjectionY("_py", iCell + 1, iCell + 1);
      dNTotalEvents = tHisto->GetEntries();
      dNOccupiedEvents = tHisto->Integral(2, tHisto->GetNbinsX());
      dNMultiEvents = tHisto->Integral(3, tHisto->GetNbinsX());

      if(dNTotalEvents)
      {
        tHistoEventCellClusterOccupancy.at(iCounterID)->SetBinContent(iCell + 1, 100.*dNOccupiedEvents/dNTotalEvents);
      }

      if(dNOccupiedEvents)
      {
        tHistoEventCellMultiClusterProbability.at(iCounterID)->SetBinContent(iCell + 1, 100.*dNMultiEvents/dNOccupiedEvents);
      }

      delete tHisto;


      tHisto = tHistoTriggerClusterMulCell.at(iCounterID)->ProjectionY("_py", iCell + 1, iCell + 1);
      dNTotalEvents = tHisto->GetEntries();
      dNOccupiedEvents = tHisto->Integral(2, tHisto->GetNbinsX());
      dNMultiEvents = tHisto->Integral(3, tHisto->GetNbinsX());

      if(dNTotalEvents)
      {
        tHistoTriggerCellClusterOccupancy.at(iCounterID)->SetBinContent(iCell + 1, 100.*dNOccupiedEvents/dNTotalEvents);
      }

      if(dNOccupiedEvents)
      {
        tHistoTriggerCellMultiClusterProbability.at(iCounterID)->SetBinContent(iCell + 1, 100.*dNMultiEvents/dNOccupiedEvents);
      }

      delete tHisto;


      tHisto = tHistoRefTrackClusterMulCell.at(iCounterID)->ProjectionY("_py", iCell + 1, iCell + 1);
      dNTotalEvents = tHisto->GetEntries();
      dNOccupiedEvents = tHisto->Integral(2, tHisto->GetNbinsX());
      dNMultiEvents = tHisto->Integral(3, tHisto->GetNbinsX());

      if(dNTotalEvents)
      {
        tHistoRefTrackCellClusterOccupancy.at(iCounterID)->SetBinContent(iCell + 1, 100.*dNOccupiedEvents/dNTotalEvents);
      }

      if(dNOccupiedEvents)
      {
        tHistoRefTrackCellMultiClusterProbability.at(iCounterID)->SetBinContent(iCell + 1, 100.*dNMultiEvents/dNOccupiedEvents);
      }

      delete tHisto;
    }


    TAxis* tXaxis = NULL;
    TGaxis* tGY2axis = NULL;

    tCurrentCanvas->cd(1);
    gPad->SetRightMargin(0.20);
    tCurrentHistogram = tHistoMeanTrackMultiplicities.at(iCounterID);
    if(tCurrentHistogram->GetEntries())
    {
      tXaxis = tCurrentHistogram->GetXaxis();
      Double_t dTrackMULExtU = tCurrentHistogram->GetBinContent(tXaxis->FindBin("EXT U"));
      Double_t dTrackMULIntC = tCurrentHistogram->GetBinContent(tXaxis->FindBin("INT C"));
      Double_t dTrackMULExtC = tCurrentHistogram->GetBinContent(tXaxis->FindBin("EXT C"));
      Double_t dTrackMULExtT = tCurrentHistogram->GetBinContent(tXaxis->FindBin("EXT T"));
      Double_t dTrackMULMemP = tCurrentHistogram->GetBinContent(tXaxis->FindBin("MEM P"));

      cout<<endl;
      cout<<Form("counter %d-%d-%d", iModuleType, iModuleIndex, iCounterIndex)<<endl;
      cout<<Form(" EXT U: MUL %8.3f flux %8.3f kHz/cm^2", dTrackMULExtU, dTrackMULExtU*dBeamIntensity*dInteractionProbability/1000./dCounterArea)<<endl;
      cout<<Form(" INT C: MUL %8.3f flux %8.3f kHz/cm^2", dTrackMULIntC, dTrackMULIntC*dBeamIntensity*dInteractionProbability/1000./dCounterArea)<<endl;
      cout<<Form(" EXT C: MUL %8.3f flux %8.3f kHz/cm^2", dTrackMULExtC, dTrackMULExtC*dBeamIntensity*dInteractionProbability/1000./dCounterArea)<<endl;
      cout<<Form(" EXT T: MUL %8.3f flux %8.3f kHz/cm^2", dTrackMULExtT, dTrackMULExtT*dBeamIntensity*dInteractionProbability/1000./dCounterArea)<<endl;
      cout<<Form(" MEM P: MUL %8.3f flux %8.3f kHz/cm^2", dTrackMULMemP, dTrackMULMemP*dBeamIntensity*dInteractionProbability/1000./dCounterArea)<<endl;

      tCurrentHistogram->LabelsDeflate("X");
      tCurrentHistogram->LabelsOption("<d", "X");
      tCurrentHistogram->SetMinimum(kY1Min);
      tCurrentHistogram->SetMaximum(kY1Max);
  //    tCurrentHistogram->SetLabelOffset(0.005);
      tCurrentHistogram->GetXaxis()->SetLabelSize(0.09);
      tCurrentHistogram->GetYaxis()->SetTitleOffset(1.20);

      tCurrentHistogram->Draw("HIST");
      gPad->Update();

      tGY2axis = new TGaxis(gPad->GetUxmax(), gPad->GetUymin(), gPad->GetUxmax(), gPad->GetUymax(),
                            bFluxAxis ? kY2Min/dCounterArea : kY2Min,
                            bFluxAxis ? kY2Max/dCounterArea : kY2Max,
                            510, "+L");

      tGY2axis->SetName(Form("tAxisMeanTrackMultiplicities_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex));
      tGY2axis->SetTitle(kY2Title.Data());
      tGY2axis->SetTitleFont(tCurrentHistogram->GetYaxis()->GetTitleFont());
      tGY2axis->SetTitleSize(tCurrentHistogram->GetYaxis()->GetTitleSize());
      tGY2axis->SetTitleOffset(tCurrentHistogram->GetYaxis()->GetTitleOffset());
      tGY2axis->SetLabelFont(tCurrentHistogram->GetYaxis()->GetLabelFont());
      tGY2axis->SetLabelOffset(tCurrentHistogram->GetYaxis()->GetLabelOffset());
      tGY2axis->SetLabelSize(tCurrentHistogram->GetYaxis()->GetLabelSize());
      tGY2axis->Draw();
    }

    tCurrentHistogram = tHistoTriggerMeanTrackMultiplicities.at(iCounterID);
    if(tCurrentHistogram->GetEntries())
    {
      tXaxis = tCurrentHistogram->GetXaxis();
      Double_t dTrackMULExtU = tCurrentHistogram->GetBinContent(tXaxis->FindBin("EXT U"));
      Double_t dTrackMULIntC = tCurrentHistogram->GetBinContent(tXaxis->FindBin("INT C"));
      Double_t dTrackMULExtC = tCurrentHistogram->GetBinContent(tXaxis->FindBin("EXT C"));
      Double_t dTrackMULExtT = tCurrentHistogram->GetBinContent(tXaxis->FindBin("EXT T"));
      Double_t dTrackMULMemP = tCurrentHistogram->GetBinContent(tXaxis->FindBin("MEM P"));

      cout<<endl;
      cout<<Form(" EXT U: MUL %8.3f trigger", dTrackMULExtU)<<endl;
      cout<<Form(" INT C: MUL %8.3f trigger", dTrackMULIntC)<<endl;
      cout<<Form(" EXT C: MUL %8.3f trigger", dTrackMULExtC)<<endl;
      cout<<Form(" EXT T: MUL %8.3f trigger", dTrackMULExtT)<<endl;
      cout<<Form(" MEM P: MUL %8.3f trigger", dTrackMULMemP)<<endl;
    }

    tCurrentHistogram = tHistoRefTrackMeanTrackMultiplicities.at(iCounterID);
    if(tCurrentHistogram->GetEntries())
    {
      tXaxis = tCurrentHistogram->GetXaxis();
      Double_t dTrackMULExtU = tCurrentHistogram->GetBinContent(tXaxis->FindBin("EXT U"));
      Double_t dTrackMULIntC = tCurrentHistogram->GetBinContent(tXaxis->FindBin("INT C"));
      Double_t dTrackMULExtC = tCurrentHistogram->GetBinContent(tXaxis->FindBin("EXT C"));
      Double_t dTrackMULExtT = tCurrentHistogram->GetBinContent(tXaxis->FindBin("EXT T"));
      Double_t dTrackMULMemP = tCurrentHistogram->GetBinContent(tXaxis->FindBin("MEM P"));

      cout<<endl;
      cout<<Form(" EXT U: MUL %8.3f reftrack", dTrackMULExtU)<<endl;
      cout<<Form(" INT C: MUL %8.3f reftrack", dTrackMULIntC)<<endl;
      cout<<Form(" EXT C: MUL %8.3f reftrack", dTrackMULExtC)<<endl;
      cout<<Form(" EXT T: MUL %8.3f reftrack", dTrackMULExtT)<<endl;
      cout<<Form(" MEM P: MUL %8.3f reftrack", dTrackMULMemP)<<endl;
    }

    tCurrentCanvas->cd(2);
    gPad->SetRightMargin(0.20);
    tCurrentHistogram = tHistoMeanTrackQualityMultiplicities.at(iCounterID);
    if(tCurrentHistogram->GetEntries())
    {
      tXaxis = tCurrentHistogram->GetXaxis();
      Double_t dTrackMULRndmC = tCurrentHistogram->GetBinContent(tXaxis->FindBin("#color[2]{RNDM C}"));
      Double_t dTrackMULTrgtC = tCurrentHistogram->GetBinContent(tXaxis->FindBin("#color[3]{TRGT C}"));
      Double_t dTrackMULTrckC = tCurrentHistogram->GetBinContent(tXaxis->FindBin("TRCK C"));

      cout<<endl;
      cout<<Form("RNDM C: MUL %8.3f flux %8.3f kHz/cm^2", dTrackMULRndmC, dTrackMULRndmC*dBeamIntensity*dInteractionProbability/1000./dCounterArea)<<endl;
      cout<<Form("TRGT C: MUL %8.3f flux %8.3f kHz/cm^2", dTrackMULTrgtC, dTrackMULTrgtC*dBeamIntensity*dInteractionProbability/1000./dCounterArea)<<endl;
      cout<<Form("TRCK C: MUL %8.3f flux %8.3f kHz/cm^2", dTrackMULTrckC, dTrackMULTrckC*dBeamIntensity*dInteractionProbability/1000./dCounterArea)<<endl;
      cout<<endl;

      tCurrentHistogram->LabelsDeflate("X");
      tCurrentHistogram->LabelsOption("<d", "X");
      tCurrentHistogram->SetMinimum(kY1Min);
      tCurrentHistogram->SetMaximum(kY1Max);
  //    tCurrentHistogram->SetLabelOffset(0.005);
      tCurrentHistogram->GetXaxis()->SetLabelSize(0.09);
      tCurrentHistogram->GetYaxis()->SetTitleOffset(1.20);

      tCurrentHistogram->Draw("HIST");
      gPad->Update();

      tGY2axis = new TGaxis(gPad->GetUxmax(), gPad->GetUymin(), gPad->GetUxmax(), gPad->GetUymax(),
                            bFluxAxis ? kY2Min/dCounterArea : kY2Min,
                            bFluxAxis ? kY2Max/dCounterArea : kY2Max,
                            510, "+L");

      tGY2axis->SetName(Form("tAxisMeanTrackQualityMultiplicities_%d_%d_%d", iModuleType, iModuleIndex, iCounterIndex));
      tGY2axis->SetTitle(kY2Title.Data());
      tGY2axis->SetTitleFont(tCurrentHistogram->GetYaxis()->GetTitleFont());
      tGY2axis->SetTitleSize(tCurrentHistogram->GetYaxis()->GetTitleSize());
      tGY2axis->SetTitleOffset(tCurrentHistogram->GetYaxis()->GetTitleOffset());
      tGY2axis->SetLabelFont(tCurrentHistogram->GetYaxis()->GetLabelFont());
      tGY2axis->SetLabelOffset(tCurrentHistogram->GetYaxis()->GetLabelOffset());
      tGY2axis->SetLabelSize(tCurrentHistogram->GetYaxis()->GetLabelSize());
      tGY2axis->Draw();
    }

    tCurrentCanvas->cd(3);
    gPad->SetLogy();

    tHistoAnyForeignTrackMultiplicityDistribution.at(iCounterID)->SetFillColor(kBlack);
    tHistoAnyForeignTrackMultiplicityDistribution.at(iCounterID)->SetMinimum(1.);

    tHistoAnyForeignTrackMultiplicityDistribution.at(iCounterID)->Draw();
    tHistoFullTargetTrackMultiplicityDistribution->Draw("SAME");
    tHistoFullRandomTrackMultiplicityDistribution->Draw("SAME");


    tCurrentCanvas->cd(4);
    gPad->SetLogy();

    tHistoAnyForeignTrackGlobalProductionZ.at(iCounterID)->SetLineColor(kBlack);

    tHistoAnyForeignTrackGlobalProductionZ.at(iCounterID)->Draw();
    tHistoTrackableTargetTrackGlobalProductionZ->Draw("SAME");
    tHistoTrackableRandomTrackGlobalProductionZ->Draw("SAME");


    tCurrentCanvas->cd(5);

    tHistoMeanLocalXYClusterFlux.at(iCounterID)->Draw("colz");
    tHistoMeanLocalXYClusterFlux.at(iCounterID)->GetZaxis()->SetMaxDigits(3);


    tCurrentCanvas->cd(6);
    gPad->SetLogy();

    tHistoForeignChargedPositionDifference.at(iCounterID)->Draw();


    tCurrentCanvas->cd(7);
    gPad->SetLogy();

    tHistoForeignChargedTimeDifference.at(iCounterID)->Draw();


    tCurrentCanvas->cd(8);
    gPad->SetLogx();
    gPad->SetLogy();

    tHistoForeignChargedLogTimeDistance.at(iCounterID)->Draw();


    tCurrentCanvas->cd(9);
    gPad->SetLogx();
    gPad->SetLogy();

    tHistoAnyPointLogTimeDistance.at(iCounterID)->Draw();


    tCurrentCanvas->cd(10);
    gPad->SetLogy();

    tHistoMeanClusterTimeInEvent.at(iCounterID)->Draw();


    tCurrentCanvas->cd(11);
    tCurrentHistogram = tHistoAnyForeignTrackCreationProcess.at(iCounterID);
    if(tCurrentHistogram->GetEntries())
    {
      tCurrentHistogram->LabelsDeflate("X");
      tCurrentHistogram->LabelsOption("<d", "X");
      tCurrentHistogram->GetXaxis()->SetLabelSize(0.08);
      tCurrentHistogram->Scale(1./tCurrentHistogram->GetEntries());
      tCurrentHistogram->Draw("HIST");
    }


    tCurrentCanvas->cd(12);
/*
    tCurrentHistogram = tHistoAnyForeignTrackCreationMaterial.at(iCounterID);
    if(tCurrentHistogram->GetEntries())
    {
      tCurrentHistogram->LabelsDeflate("X");
      tCurrentHistogram->LabelsOption("<d", "X");
      tCurrentHistogram->GetXaxis()->SetLabelSize(0.08);
      tCurrentHistogram->Scale(1./tCurrentHistogram->GetEntries());
      tCurrentHistogram->Draw("HIST");
    }
*/
    tCurrentHistogram = tHistoAnyForeignTrackCreationProcessAndSpecies.at(iCounterID);
    if(tCurrentHistogram->GetEntries())
    {
      tCurrentHistogram->LabelsDeflate("X");
      tCurrentHistogram->LabelsDeflate("Y");
      tCurrentHistogram->LabelsOption("<d", "X");
      tCurrentHistogram->LabelsOption("<v", "Y");
      tCurrentHistogram->GetXaxis()->SetLabelSize(0.08);
      tCurrentHistogram->GetYaxis()->SetLabelSize(0.09);
      tCurrentHistogram->Scale(1./tCurrentHistogram->GetEntries());
      tCurrentHistogram->SetMarkerSize(2);
      tCurrentHistogram->Draw("col text");
      gPad->Update();
      gPad->Modified();
    }

    tCurrentCanvas->cd(13);
    tCurrentHistogram = tHistoAnyForeignTrackCreationProcessAndMaterial.at(iCounterID);
    if(tCurrentHistogram->GetEntries())
    {
      tCurrentHistogram->LabelsDeflate("X");
      tCurrentHistogram->LabelsDeflate("Y");
      tCurrentHistogram->LabelsOption("<d", "X");
      tCurrentHistogram->LabelsOption("<v", "Y");
      tCurrentHistogram->GetXaxis()->SetLabelSize(0.08);
      tCurrentHistogram->GetYaxis()->SetLabelSize(0.09);
      tCurrentHistogram->Scale(1./tCurrentHistogram->GetEntries());
      tCurrentHistogram->SetMarkerSize(2);
      tCurrentHistogram->Draw("col text");
      gPad->Update();
      gPad->Modified();
    }


    tCurrentCanvas->cd(14);
    gPad->SetLogy();

    tHistoClusterSize.at(iCounterID)->Draw();


    tCurrentCanvas->cd(15);
    gPad->SetLogy();

    tHistoClusterKinshipDegree.at(iCounterID)->Draw();


    tCurrentCanvas->cd(16);
    tCurrentHistogram = tHistoAnyDomesticTrackCreationProcess.at(iCounterID);
    if(tCurrentHistogram->GetEntries())
    {
      tCurrentHistogram->LabelsDeflate("X");
      tCurrentHistogram->LabelsOption("<d", "X");
      tCurrentHistogram->GetXaxis()->SetLabelSize(0.08);
      tCurrentHistogram->Scale(1./tCurrentHistogram->GetEntries());
      tCurrentHistogram->Draw("HIST");
    }


    tCurrentCanvas->cd(17);
/*
    tCurrentHistogram = tHistoAnyDomesticTrackCreationMaterial.at(iCounterID);
    if(tCurrentHistogram->GetEntries())
    {
      tCurrentHistogram->LabelsDeflate("X");
      tCurrentHistogram->LabelsOption("<d", "X");
      tCurrentHistogram->GetXaxis()->SetLabelSize(0.08);
      tCurrentHistogram->Scale(1./tCurrentHistogram->GetEntries());
      tCurrentHistogram->Draw("HIST");
    }
*/
    tCurrentHistogram = tHistoAnyDomesticTrackCreationProcessAndSpecies.at(iCounterID);
    if(tCurrentHistogram->GetEntries())
    {
      tCurrentHistogram->LabelsDeflate("X");
      tCurrentHistogram->LabelsDeflate("Y");
      tCurrentHistogram->LabelsOption("<d", "X");
      tCurrentHistogram->LabelsOption("<v", "Y");
      tCurrentHistogram->GetXaxis()->SetLabelSize(0.08);
      tCurrentHistogram->GetYaxis()->SetLabelSize(0.09);
      tCurrentHistogram->Scale(1./tCurrentHistogram->GetEntries());
      tCurrentHistogram->SetMarkerSize(2);
      tCurrentHistogram->Draw("col text");
      gPad->Update();
      gPad->Modified();
    }

    tCurrentCanvas->cd(18);
    tCurrentHistogram = tHistoAnyDomesticTrackCreationProcessAndMaterial.at(iCounterID);
    if(tCurrentHistogram->GetEntries())
    {
      tCurrentHistogram->LabelsDeflate("X");
      tCurrentHistogram->LabelsDeflate("Y");
      tCurrentHistogram->LabelsOption("<d", "X");
      tCurrentHistogram->LabelsOption("<v", "Y");
      tCurrentHistogram->GetXaxis()->SetLabelSize(0.08);
      tCurrentHistogram->GetYaxis()->SetLabelSize(0.09);
      tCurrentHistogram->Scale(1./tCurrentHistogram->GetEntries());
      tCurrentHistogram->SetMarkerSize(2);
      tCurrentHistogram->Draw("col text");
      gPad->Update();
      gPad->Modified();
    }


    tCurrentCanvas->cd(19);
//    gPad->SetLogy();
//    tHistoClusterInternalTimeDifference.at(iCounterID)->Draw();
    tHistoEventCellClusterOccupancy.at(iCounterID)->GetYaxis()->SetRangeUser(0., 40.);
    tHistoEventCellClusterOccupancy.at(iCounterID)->Draw();
    tHistoTriggerCellClusterOccupancy.at(iCounterID)->SetLineColor(kRed);
    tHistoTriggerCellClusterOccupancy.at(iCounterID)->Draw("SAME");
    tHistoRefTrackCellClusterOccupancy.at(iCounterID)->SetLineColor(kGreen);
    tHistoRefTrackCellClusterOccupancy.at(iCounterID)->Draw("SAME");

    tCurrentCanvas->cd(20);
//    gPad->SetLogy();
//    tHistoClusterInternalPositionDifference.at(iCounterID)->Draw();
    tHistoEventCellMultiClusterProbability.at(iCounterID)->GetYaxis()->SetRangeUser(0., 20.);
    tHistoEventCellMultiClusterProbability.at(iCounterID)->Draw();
    tHistoTriggerCellMultiClusterProbability.at(iCounterID)->SetLineColor(kRed);
    tHistoTriggerCellMultiClusterProbability.at(iCounterID)->Draw("SAME");
    tHistoRefTrackCellMultiClusterProbability.at(iCounterID)->SetLineColor(kGreen);
    tHistoRefTrackCellMultiClusterProbability.at(iCounterID)->Draw("SAME");

    cout<<endl;
    cout<<Form("UrQMD      occupancy: %4.1f %%", tHistoEventCellClusterOccupancy.at(iCounterID)->Integral()/tCounter->fiNCells)<<endl;
    cout<<Form("trigger    occupancy: %4.1f %%", tHistoTriggerCellClusterOccupancy.at(iCounterID)->Integral()/tCounter->fiNCells)<<endl;
    cout<<Form("reftrack   occupancy: %4.1f %%", tHistoRefTrackCellClusterOccupancy.at(iCounterID)->Integral()/tCounter->fiNCells)<<endl;
    cout<<Form("UrQMD    multi-track: %4.1f %%", tHistoEventCellMultiClusterProbability.at(iCounterID)->Integral()/tCounter->fiNCells)<<endl;
    cout<<Form("trigger  multi-track: %4.1f %%", tHistoTriggerCellMultiClusterProbability.at(iCounterID)->Integral()/tCounter->fiNCells)<<endl;
    cout<<Form("reftrack multi-track: %4.1f %%", tHistoRefTrackCellMultiClusterProbability.at(iCounterID)->Integral()/tCounter->fiNCells)<<endl;
    cout<<endl;

    if(bBatchMode)
    {
      tCurrentCanvas->SaveAs(Form("pl_mc_inspection_%d_%d_%d.pdf", iModuleType, iModuleIndex, iCounterIndex));

      delete tCurrentCanvas;
    }

  }

 
  TCanvas* tGeneralCanvas = new TCanvas("tGeneralCanvas", "General inspection", 0, 0, 700, 350);
  tGeneralCanvas->Divide(2, 1);

  tGeneralCanvas->cd(1);
  gPad->SetLogz();
  tHistoTrackableTargetTrackAcceptance->Draw("colz");
  TF1* tFunction = new TF1("fun_pty_mean", fun_pty_tantheta, 0., 4., 1);
  tFunction->SetParameter(0, kdTyMean);
  tFunction->SetLineColor(kBlue);
  tFunction->Draw("same");

  tGeneralCanvas->cd(2);
  tCurrentHistogram = tHistoTrackableTargetTrackCreationProcessAndSpecies;
  if(tCurrentHistogram->GetEntries())
  {
    tCurrentHistogram->LabelsDeflate("X");
    tCurrentHistogram->LabelsDeflate("Y");
    tCurrentHistogram->LabelsOption("<d", "X");
    tCurrentHistogram->LabelsOption("<v", "Y");
    tCurrentHistogram->GetXaxis()->SetLabelSize(0.08);
    tCurrentHistogram->GetYaxis()->SetLabelSize(0.09);
    tCurrentHistogram->Scale(1./tCurrentHistogram->GetEntries());
    tCurrentHistogram->SetMarkerSize(2);
    tCurrentHistogram->Draw("col text");
    gPad->Update();
    gPad->Modified();
  }

  if(bBatchMode)
  {
    tGeneralCanvas->SaveAs("pl_mc_inspection_general.pdf");

    delete tGeneralCanvas;
  }


  cout<<endl;
  cout<<Form("beam intensity: %15.3f kHz", dBeamIntensity/1000.)<<endl;
  cout<<Form("collision rate: %15.3f kHz", dBeamIntensity*dInteractionProbability/1000.)<<endl;
  cout<<Form("trigger   rate: %15.3f kHz", dBeamIntensity*dInteractionProbability*iNTriggerEvents/iNEvents/1000.)<<endl;
  cout<<Form("selector  rate: %15.3f kHz", dBeamIntensity*dInteractionProbability*iNRefTrackEvents/iNEvents/1000.)<<endl;
  cout<<endl;

  return;
}


Int_t CountPads(TVirtualPad* tPad)
{
  if(!tPad)
  {
    return 0;
  }

  Int_t iNPads = 0;
  TObject* tObject;

  TIter next(tPad->GetListOfPrimitives());
  while((tObject = next()))
  {
    if(tObject->InheritsFrom(TVirtualPad::Class()))
    {
      iNPads++;
    }
  }

  return iNPads;
}

void AddCounter(Int_t iModuleType, Int_t iModuleIndex, Int_t iCounterIndex)
{
  cout<<Form("Adding counter %d-%d-%d ...", iModuleType, iModuleIndex, iCounterIndex)<<endl;
  gCounters.push_back(std::tuple<Int_t, Int_t, Int_t>(iModuleType, iModuleIndex, iCounterIndex));
}

Bool_t FindTofNode(const char* cGeoFileName)
{
  TFile::Open(cGeoFileName,"read");

  if(!gFile)
  {
    cout<<"Could not open geometry file. Abort macro execution."<<endl;
    return kFALSE;
  }

  gFile->FindObjectAny("FAIRGeom");
  if(!gGeoManager)
  {
    cout<<"Could not retrieve TGeoManager 'FAIRGeom' from geometry file."<<endl;
    return kFALSE;
  }

  gGeoHandler = new CbmTofGeoHandler();

  if(!gGeoHandler->Init(kFALSE))
  {
    cout<<"Failed to initialize 'CbmTofGeoHandler'."<<endl;
    return kFALSE;
  }

  Bool_t bFoundTofNode(kFALSE);
  gGeoManager->CdTop();

  TObjArray* tNodes = gGeoManager->GetCurrentNode()->GetNodes();
  for(Int_t iNode = 0; iNode < tNodes->GetEntriesFast(); iNode++)
  {
    TGeoNode* tNode = dynamic_cast<TGeoNode*>(tNodes->At(iNode));
    if(TString(tNode->GetName()).Contains("tof",TString::kIgnoreCase))
    {
      gGeoManager->CdDown(iNode);
      bFoundTofNode = kTRUE;
      break;
    }
  }

  if(!bFoundTofNode)
  {
    cout<<"Could not retrieve 'tof' node from TGeoManager."<<endl;
    return kFALSE;
  }

  gTofNodePath = gGeoManager->GetPath();

  return kTRUE;
}

void ExpandNode(TGeoNode* tMotherNode)
{
  TObjArray* tDaughterNodesList = tMotherNode->GetNodes();

  for(Int_t iNode = 0; iNode < tDaughterNodesList->GetEntriesFast(); iNode++)
  {
    TGeoNode* tDaughterNode = dynamic_cast<TGeoNode*>(tDaughterNodesList->At(iNode));

    // Extract the current module type and module index from the module node
    if(TString(tDaughterNode->GetName()).Contains("module"))
    {
      giCurrentModuleType = -1;
      giCurrentModuleIndex = -1;
      giCurrentCounterIndex = -1;

      gCurrentNodePath = gTofNodePath + "/" + (TString)tDaughterNode->GetName();

      boost::regex rgx(".*_(\\d+)_.*");
      boost::cmatch match;
      if( boost::regex_search(tDaughterNode->GetName(), match, rgx) )
      {
        giCurrentModuleType = boost::lexical_cast<Int_t>(match[1]);
      }

      giCurrentModuleIndex = tDaughterNode->GetNumber();

      gCurrentBoxNodePath = gCurrentNodePath;
    }
    // Extract the current counter index from the counter node and find the
    // central node in the counter volume (a glass plate OR a gap) to serve as
    // reference coordinate system
    else if(TString(tDaughterNode->GetName()).Contains("counter"))
    {
      gCurrentNodePath += "/" + (TString)tDaughterNode->GetName();

      giCurrentCounterIndex = tDaughterNode->GetNumber();

      TString tCurrentCentralNodePath = gCurrentNodePath;

      Int_t iNGaps(0);

      for(Int_t iDaughter = 0; iDaughter < tDaughterNode->GetNdaughters(); iDaughter++)
      {
        if(TString(tDaughterNode->GetDaughter(iDaughter)->GetName()).Contains("Gap"))
        {
          iNGaps++;
        }
      }

      if(0 == iNGaps%2)
      {
        tCurrentCentralNodePath += TString::Format("/tof_glass_%d",iNGaps/2);
      }
      else
      {
        tCurrentCentralNodePath += TString::Format("/Gap_%d",(iNGaps-1)/2);
      }

      for(auto &itCounter : gCounters)
      {
        if(std::get<0>(itCounter) == giCurrentModuleType &&
           std::get<1>(itCounter) == giCurrentModuleIndex &&
           std::get<2>(itCounter) == giCurrentCounterIndex)
        {
          CbmTofDetectorInfo tCounterInfo(kTof, giCurrentModuleType, giCurrentModuleIndex, giCurrentCounterIndex, 0, 0);
          Int_t iUniqueCounterId = gGeoHandler->GetDetIdPointer()->SetDetectorInfo(tCounterInfo);

          TCounter* tCounter = new TCounter(giCurrentModuleType, giCurrentModuleIndex, giCurrentCounterIndex,
                                            tDaughterNode, new TGeoPhysicalNode(tCurrentCentralNodePath.Data()),
                                            new TGeoPhysicalNode(gCurrentBoxNodePath.Data()));

          gCounterMap.emplace(iUniqueCounterId, tCounter);
        }
      }

    }
    else
    {
      gCurrentNodePath += "/" + (TString)tDaughterNode->GetName();
    }

    // Expand nodes recursively
    if(tDaughterNode->GetNdaughters() > 0) { ExpandNode(tDaughterNode); }

    // Remove a node's name from the current node path upon completing a
    // recursion step
    gCurrentNodePath.ReplaceAll("/"+(TString)tDaughterNode->GetName(), "");
  }
}

void GetPdgName(Int_t iPdgCode, const char*& cPdgName)
{
  auto itPdg = giPdgMap.find(iPdgCode);
  if(itPdg != giPdgMap.end())
  {
    cPdgName = itPdg->second;
  }
  else
  {
    cPdgName = "XXX";
    cout<<"unknown PDG code: "<<iPdgCode<<endl;
  }
}

void GetMaterialName(const char* cMaterial, const char*& cMaterialName)
{
  auto itMaterial = gcMaterialMap.find(cMaterial);
  if(itMaterial != gcMaterialMap.end())
  {
    cMaterialName = itMaterial->second;
  }
  else
  {
    cMaterialName = "XXX";
    cout<<"unknown material: "<<cMaterial<<endl;
  }
}

void GetProcessName(const char* cProcess, const char*& cProcessName)
{
  auto itProcess = gcProcessMap.find(cProcess);
  if(itProcess != gcProcessMap.end())
  {
    cProcessName = itProcess->second;
  }
  else
  {
    cProcessName = "XXX";
    cout<<"unknown process: "<<cProcess<<endl;
  }
}

Double_t fun_pty_tantheta(Double_t* x, Double_t* par)
{
  Double_t dRadicand = 1. - TMath::SinH(x[0])*TMath::SinH(x[0])*TMath::Power(par[0], 2.);

  if(0 < dRadicand)
  {
    return TMath::SinH(x[0])*par[0]/TMath::Sqrt(dRadicand);
  }
  else
  {
    return 2.6;
  }
}