// -----------------------------------------------------------------------------
// ----- pl_response_tb.C                                                  -----
// -----                                                                   -----
// ----- created by C. Simon on 2017-09-06                                 -----
// -----                                                                   -----
// -----------------------------------------------------------------------------
#if !defined(__CLING__) || defined(__ROOTCLING__)
#include <Riostream.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TH2.h>
#include <TH3.h>
#include <TF1.h>
#include <TProfile.h>
#include <TEfficiency.h>
#include <TGraphAsymmErrors.h>
#include <TRegexp.h>
#include <TMath.h>
#include <TPaveStats.h>
#include <TSystem.h>
#include <TSystemDirectory.h>
#include <TGeoManager.h>
#include <TKey.h>
#include <TStyle.h>
#include <TGeoPhysicalNode.h>
#include <Math/DistFunc.h>
#include <Math/QuantFuncMathMore.h>
#include <TVectorD.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"
#endif
#include <boost/regex.hpp>
#include <boost/lexical_cast.hpp>

using namespace std;

Int_t giNumberGlassPlates(0);
Int_t giNumberGasGaps(0);
Double_t gdGlassThickness(0.);
Double_t gdGapSize(0.);
Bool_t gbPadReadout(kFALSE);
Int_t giNumberReadoutElectrodes(0);
Double_t gdGapXWidth(0.);
Double_t gdGapYWidth(0.);
Double_t gdCellXWidth(0.);
Double_t gdCellYWidth(0.);
Double_t gdCellAcrossWidth(0.);
Double_t gdCellAlongWidth(0.);
Double_t gdGapAcrossWidth(0.);
Double_t gdGapAlongWidth(0.);
Bool_t gbCellsAlongX(kFALSE);
Bool_t gbCellIndexingAlongAxis(kFALSE);

TString gtCurrentNodePath;
TString gtTofNodePath;
TString gtCurrentCentralNodePath;
TString gtCentralNodePath;
Int_t giCurrentModuleType(-1);
Int_t giCurrentModuleIndex(-1);
Int_t giCurrentCounterIndex(-1);
Int_t giModuleType(-1);
Int_t giModuleIndex(-1);
Int_t giCounterIndex(-1);
TGeoNode* gtCounterNode;
TGeoPhysicalNode* gtCentralPN(NULL);

Int_t countpads(TVirtualPad *pad);
Double_t gpb(Double_t* x, Double_t* p);
Double_t chisqpdf(Double_t* x, Double_t* p);
Double_t chisqcdf(Double_t* x, Double_t* p);
Bool_t FindTofNode(const char* cGeoTag);
Bool_t GetCounterProperties();
void ExpandNode(TGeoNode* tMotherNode);

void pl_response_tb(const char* cFileName, Int_t iNMergedFiles, Bool_t bBatchMode, const char* cGeoTag, Bool_t bSpotResponse)
{

  gStyle->SetOptStat("emr");
  gErrorIgnoreLevel = kWarning;

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

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

  boost::regex rgx("(.*)_(.*)_(.*)");

  for(const auto&& obj: *gFile->GetListOfKeys())
  {
    TString tCounterTag = dynamic_cast<TKey*>(obj)->GetTitle();

    if(!gFile->cd(Form("%s", tCounterTag.Data())))
    {
      cout<<Form("Could not find directory '%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }


    boost::cmatch match;
    if(boost::regex_search(tCounterTag.Data(), match, rgx))
    {
      if(4 == match.size())
      {
        giModuleType   = boost::lexical_cast<Int_t>(match[1]);
        giModuleIndex  = boost::lexical_cast<Int_t>(match[2]);
        giCounterIndex = boost::lexical_cast<Int_t>(match[3]);
      }
      else
      {
        return;
      }
    }
    else
    {
      return;
    }

    if(!GetCounterProperties())
    {
      cout<<Form("Failed to retrieve geometrical properties of counter %s.", tCounterTag.Data())<<endl;
    }


    TH2F* tClusterSizes(NULL);
    TH2F* tStripToTSpectra(NULL);
    TH2F* tStripPositionSpectra(NULL);
    TH2F* tPointPositionSpectra(NULL);
    TH2F* tHitPositionSpectra(NULL);
    TH3F* tStripTimeResidualWalk(NULL);
    TH1F* tHitAcrossPositionResidual(NULL);
    TH1F* tHitAlongPositionResidual(NULL);
    TH1F* tHitTimeResidual(NULL);
    TH1F* tHitChiSquare(NULL);
    TH1F* tCorrectedHitChiSquare(NULL);
    TH2F* tClusterSizeEvolution(NULL);
    TH2F* tAcrossPositionResidualEvolution(NULL);
    TH2F* tAlongPositionResidualEvolution(NULL);
    TH2F* tTimeResidualEvolution(NULL);
    TH2F* tCentralStripToTEvolution(NULL);
    TEfficiency* tIntegratedEfficiency(NULL);
    TEfficiency* tChiSquareEfficiency(NULL);
    TEfficiency* tCorrectedChiSquareEfficiency(NULL);
    TH2F* tClusterSizesSPH(NULL);
    TH2F* tStripToTSpectraSPH(NULL);
    TH2F* tStripPositionSpectraSPH(NULL);
    TH2F* tHitPositionSpectraSPH(NULL);
    TH3F* tStripTimeResidualWalkSPH(NULL);
    TH1F* tHitAcrossPositionResidualSPH(NULL);
    TH1F* tHitAlongPositionResidualSPH(NULL);
    TH1F* tHitTimeResidualSPH(NULL);
    TH1F* tHitChiSquareSPH(NULL);
    TH1F* tCorrectedHitChiSquareSPH(NULL);
    TH2F* tClusterSizeEvolutionSPH(NULL);
    TH2F* tAcrossPositionResidualEvolutionSPH(NULL);
    TH2F* tAlongPositionResidualEvolutionSPH(NULL);
    TH2F* tTimeResidualEvolutionSPH(NULL);
    TH2F* tCentralStripToTEvolutionSPH(NULL);
    TH1F* tDistanceToLastPoint(NULL);
    TH1F* tDistanceToLocalLastPoint(NULL);
    TH1F* tDistanceToLastHitUnbiased(NULL);
    TH1F* tDistanceToLastHitBiased(NULL);
    TEfficiency* tEfficiencyToLastHitUnbiased(NULL);
    TEfficiency* tEfficiencyToLastHitBiased(NULL);
    TH2F* tClusterSizeToLastHitUnbiased(NULL);
    TH2F* tClusterSizeToLastHitBiased(NULL);
    TH2F* tTimeResidualToLastHitUnbiased(NULL);
    TH2F* tTimeResidualToLastHitBiased(NULL);
    TEfficiency* tLocalSurfaceFlux(NULL);
    TEfficiency* tChannelDigiRates(NULL);
    TH1F* tNHitsPerPoint(NULL);
    TH1F* tNPointsPerHit(NULL);

    gDirectory->GetObject(Form("tClusterSizes_%s", tCounterTag.Data()), tClusterSizes);
    gDirectory->GetObject(Form("tStripToTSpectra_%s", tCounterTag.Data()), tStripToTSpectra);
    gDirectory->GetObject(Form("tStripPositionSpectra_%s", tCounterTag.Data()), tStripPositionSpectra);
    gDirectory->GetObject(Form("tPointPositionSpectra_%s", tCounterTag.Data()), tPointPositionSpectra);
    gDirectory->GetObject(Form("tHitPositionSpectra_%s", tCounterTag.Data()), tHitPositionSpectra);
    gDirectory->GetObject(Form("tStripTimeResidualWalk_%s", tCounterTag.Data()), tStripTimeResidualWalk);
    gDirectory->GetObject(Form("tHitAcrossPositionResidual_%s", tCounterTag.Data()), tHitAcrossPositionResidual);
    gDirectory->GetObject(Form("tHitAlongPositionResidual_%s", tCounterTag.Data()), tHitAlongPositionResidual);
    gDirectory->GetObject(Form("tHitTimeResidual_%s", tCounterTag.Data()), tHitTimeResidual);
    gDirectory->GetObject(Form("tHitChiSquare_%s", tCounterTag.Data()), tHitChiSquare);
    gDirectory->GetObject(Form("tCorrectedHitChiSquare_%s", tCounterTag.Data()), tCorrectedHitChiSquare);
    gDirectory->GetObject(Form("tClusterSizeEvolution_%s", tCounterTag.Data()), tClusterSizeEvolution);
    gDirectory->GetObject(Form("tAcrossPositionResidualEvolution_%s", tCounterTag.Data()), tAcrossPositionResidualEvolution);
    gDirectory->GetObject(Form("tAlongPositionResidualEvolution_%s", tCounterTag.Data()), tAlongPositionResidualEvolution);
    gDirectory->GetObject(Form("tTimeResidualEvolution_%s", tCounterTag.Data()), tTimeResidualEvolution);
    gDirectory->GetObject(Form("tCentralStripToTEvolution_%s", tCounterTag.Data()), tCentralStripToTEvolution);
    gDirectory->GetObject(Form("tIntegratedEfficiency_%s", tCounterTag.Data()), tIntegratedEfficiency);
    gDirectory->GetObject(Form("tChiSquareEfficiency_%s", tCounterTag.Data()), tChiSquareEfficiency);
    gDirectory->GetObject(Form("tCorrectedChiSquareEfficiency_%s", tCounterTag.Data()), tCorrectedChiSquareEfficiency);
    gDirectory->GetObject(Form("tClusterSizesSPH_%s", tCounterTag.Data()), tClusterSizesSPH);
    gDirectory->GetObject(Form("tStripToTSpectraSPH_%s", tCounterTag.Data()), tStripToTSpectraSPH);
    gDirectory->GetObject(Form("tStripPositionSpectraSPH_%s", tCounterTag.Data()), tStripPositionSpectraSPH);
    gDirectory->GetObject(Form("tHitPositionSpectraSPH_%s", tCounterTag.Data()), tHitPositionSpectraSPH);
    gDirectory->GetObject(Form("tStripTimeResidualWalkSPH_%s", tCounterTag.Data()), tStripTimeResidualWalkSPH);
    gDirectory->GetObject(Form("tHitAcrossPositionResidualSPH_%s", tCounterTag.Data()), tHitAcrossPositionResidualSPH);
    gDirectory->GetObject(Form("tHitAlongPositionResidualSPH_%s", tCounterTag.Data()), tHitAlongPositionResidualSPH);
    gDirectory->GetObject(Form("tHitTimeResidualSPH_%s", tCounterTag.Data()), tHitTimeResidualSPH);
    gDirectory->GetObject(Form("tHitChiSquareSPH_%s", tCounterTag.Data()), tHitChiSquareSPH);
    gDirectory->GetObject(Form("tCorrectedHitChiSquareSPH_%s", tCounterTag.Data()), tCorrectedHitChiSquareSPH);
    gDirectory->GetObject(Form("tClusterSizeEvolutionSPH_%s", tCounterTag.Data()), tClusterSizeEvolutionSPH);
    gDirectory->GetObject(Form("tAcrossPositionResidualEvolutionSPH_%s", tCounterTag.Data()), tAcrossPositionResidualEvolutionSPH);
    gDirectory->GetObject(Form("tAlongPositionResidualEvolutionSPH_%s", tCounterTag.Data()), tAlongPositionResidualEvolutionSPH);
    gDirectory->GetObject(Form("tTimeResidualEvolutionSPH_%s", tCounterTag.Data()), tTimeResidualEvolutionSPH);
    gDirectory->GetObject(Form("tCentralStripToTEvolutionSPH_%s", tCounterTag.Data()), tCentralStripToTEvolutionSPH);
    gDirectory->GetObject(Form("tDistanceToLastPoint_%s", tCounterTag.Data()), tDistanceToLastPoint);
    gDirectory->GetObject(Form("tDistanceToLocalLastPoint_%s", tCounterTag.Data()), tDistanceToLocalLastPoint);
    gDirectory->GetObject(Form("tDistanceToLastHitUnbiased_%s", tCounterTag.Data()), tDistanceToLastHitUnbiased);
    gDirectory->GetObject(Form("tDistanceToLastHitBiased_%s", tCounterTag.Data()), tDistanceToLastHitBiased);
    gDirectory->GetObject(Form("tEfficiencyToLastHitUnbiased_%s", tCounterTag.Data()), tEfficiencyToLastHitUnbiased);
    gDirectory->GetObject(Form("tEfficiencyToLastHitBiased_%s", tCounterTag.Data()), tEfficiencyToLastHitBiased);
    gDirectory->GetObject(Form("tClusterSizeToLastHitUnbiased_%s", tCounterTag.Data()), tClusterSizeToLastHitUnbiased);
    gDirectory->GetObject(Form("tClusterSizeToLastHitBiased_%s", tCounterTag.Data()), tClusterSizeToLastHitBiased);
    gDirectory->GetObject(Form("tTimeResidualToLastHitUnbiased_%s", tCounterTag.Data()), tTimeResidualToLastHitUnbiased);
    gDirectory->GetObject(Form("tTimeResidualToLastHitBiased_%s", tCounterTag.Data()), tTimeResidualToLastHitBiased);
    gDirectory->GetObject(Form("tLocalSurfaceFlux_%s", tCounterTag.Data()), tLocalSurfaceFlux);
    gDirectory->GetObject(Form("tChannelDigiRates_%s", tCounterTag.Data()), tChannelDigiRates);
    gDirectory->GetObject(Form("tNHitsPerPoint_%s", tCounterTag.Data()), tNHitsPerPoint);
    gDirectory->GetObject(Form("tNPointsPerHit_%s", tCounterTag.Data()), tNPointsPerHit);

    if(!tClusterSizes)
    {
      cout<<Form("Could not find histogram 'tClusterSizes_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tStripToTSpectra)
    {
      cout<<Form("Could not find histogram 'tStripToTSpectra_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tStripPositionSpectra)
    {
      cout<<Form("Could not find histogram 'tStripPositionSpectra_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tPointPositionSpectra)
    {
      cout<<Form("Could not find histogram 'tPointPositionSpectra_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tHitPositionSpectra)
    {
      cout<<Form("Could not find histogram 'tHitPositionSpectra_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tStripTimeResidualWalk)
    {
      cout<<Form("Could not find histogram 'tStripTimeResidualWalk_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tHitAcrossPositionResidual)
    {
      cout<<Form("Could not find histogram 'tHitAcrossPositionResidual_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tHitAlongPositionResidual)
    {
      cout<<Form("Could not find histogram 'tHitAlongPositionResidual_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tHitTimeResidual)
    {
      cout<<Form("Could not find histogram 'tHitTimeResidual_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tHitChiSquare)
    {
      cout<<Form("Could not find histogram 'tHitChiSquare_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tCorrectedHitChiSquare)
    {
      cout<<Form("Could not find histogram 'tCorrectedHitChiSquare_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tClusterSizeEvolution)
    {
      cout<<Form("Could not find histogram 'tClusterSizeEvolution_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tAcrossPositionResidualEvolution)
    {
      cout<<Form("Could not find histogram 'tAcrossPositionResidualEvolution_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tAlongPositionResidualEvolution)
    {
      cout<<Form("Could not find histogram 'tAlongPositionResidualEvolution_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tTimeResidualEvolution)
    {
      cout<<Form("Could not find histogram 'tTimeResidualEvolution_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tCentralStripToTEvolution)
    {
      cout<<Form("Could not find histogram 'tCentralStripToTEvolution_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tIntegratedEfficiency)
    {
      cout<<Form("Could not find histogram 'tIntegratedEfficiency_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tChiSquareEfficiency)
    {
      cout<<Form("Could not find histogram 'tChiSquareEfficiency_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tCorrectedChiSquareEfficiency)
    {
      cout<<Form("Could not find histogram 'tCorrectedChiSquareEfficiency_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tClusterSizesSPH)
    {
      cout<<Form("Could not find histogram 'tClusterSizesSPH_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tStripToTSpectraSPH)
    {
      cout<<Form("Could not find histogram 'tStripToTSpectraSPH_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tStripPositionSpectraSPH)
    {
      cout<<Form("Could not find histogram 'tStripPositionSpectraSPH_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tHitPositionSpectraSPH)
    {
      cout<<Form("Could not find histogram 'tHitPositionSpectraSPH_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tStripTimeResidualWalkSPH)
    {
      cout<<Form("Could not find histogram 'tStripTimeResidualWalkSPH_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tHitAcrossPositionResidualSPH)
    {
      cout<<Form("Could not find histogram 'tHitAcrossPositionResidualSPH_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tHitAlongPositionResidualSPH)
    {
      cout<<Form("Could not find histogram 'tHitAlongPositionResidualSPH_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tHitTimeResidualSPH)
    {
      cout<<Form("Could not find histogram 'tHitTimeResidualSPH_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tHitChiSquareSPH)
    {
      cout<<Form("Could not find histogram 'tHitChiSquareSPH_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tCorrectedHitChiSquareSPH)
    {
      cout<<Form("Could not find histogram 'tCorrectedHitChiSquareSPH_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tClusterSizeEvolutionSPH)
    {
      cout<<Form("Could not find histogram 'tClusterSizeEvolutionSPH_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tAcrossPositionResidualEvolutionSPH)
    {
      cout<<Form("Could not find histogram 'tAcrossPositionResidualEvolutionSPH_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tAlongPositionResidualEvolutionSPH)
    {
      cout<<Form("Could not find histogram 'tAlongPositionResidualEvolutionSPH_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tTimeResidualEvolutionSPH)
    {
      cout<<Form("Could not find histogram 'tTimeResidualEvolutionSPH_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tCentralStripToTEvolutionSPH)
    {
      cout<<Form("Could not find histogram 'tCentralStripToTEvolutionSPH_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tDistanceToLastPoint)
    {
      cout<<Form("Could not find histogram 'tDistanceToLastPoint_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tDistanceToLocalLastPoint)
    {
      cout<<Form("Could not find histogram 'tDistanceToLocalLastPoint_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tDistanceToLastHitUnbiased)
    {
      cout<<Form("Could not find histogram 'tDistanceToLastHitUnbiased_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tDistanceToLastHitBiased)
    {
      cout<<Form("Could not find histogram 'tDistanceToLastHitBiased_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tEfficiencyToLastHitUnbiased)
    {
      cout<<Form("Could not find histogram 'tEfficiencyToLastHitUnbiased_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tEfficiencyToLastHitBiased)
    {
      cout<<Form("Could not find histogram 'tEfficiencyToLastHitBiased_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tClusterSizeToLastHitUnbiased)
    {
      cout<<Form("Could not find histogram 'tClusterSizeToLastHitUnbiased_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tClusterSizeToLastHitBiased)
    {
      cout<<Form("Could not find histogram 'tClusterSizeToLastHitBiased_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tTimeResidualToLastHitUnbiased)
    {
      cout<<Form("Could not find histogram 'tTimeResidualToLastHitUnbiased_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tTimeResidualToLastHitBiased)
    {
      cout<<Form("Could not find histogram 'tTimeResidualToLastHitBiased_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tLocalSurfaceFlux)
    {
      cout<<Form("Could not find histogram 'tLocalSurfaceFlux_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tChannelDigiRates)
    {
      cout<<Form("Could not find histogram 'tChannelDigiRates_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tNHitsPerPoint)
    {
      cout<<Form("Could not find histogram 'tNHitsPerPoint_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }

    if(!tNPointsPerHit)
    {
      cout<<Form("Could not find histogram 'tNPointsPerHit_%s' in histogram file.", tCounterTag.Data())<<endl;
      return;
    }


    if(0 == tClusterSizes->GetEntries())
    {
      cout<<Form("Cluster size histogram of counter %s does not contain any entry.", tCounterTag.Data())<<endl;
      cout<<"Most likely, a one-sided readout was configured. No QA histograms exist for this scenario (yet)."<<endl;
      cout<<"Skip counter."<<endl;

      gFile->cd("../");

      continue;
    }


    const_cast<TH1*>(tLocalSurfaceFlux->GetPassedHistogram())->Scale(1./iNMergedFiles);
    const_cast<TH1*>(tLocalSurfaceFlux->GetTotalHistogram())->Scale(1./iNMergedFiles);
    const_cast<TH1*>(tChannelDigiRates->GetPassedHistogram())->Scale(1./iNMergedFiles);
    const_cast<TH1*>(tChannelDigiRates->GetTotalHistogram())->Scale(1./iNMergedFiles);


    Int_t iCounterMainStripBin = tClusterSizes->GetNbinsX()/2;

    if(bSpotResponse)
    {
      Double_t dLargestStripStatistics = tClusterSizes->Integral(iCounterMainStripBin, iCounterMainStripBin, 1, tClusterSizes->GetNbinsY());

      for(Int_t iStrip = 0; iStrip < tClusterSizes->GetNbinsX(); iStrip++)
      {
        Double_t dStripStatistics = tClusterSizes->Integral(iStrip + 1, iStrip + 1, 1, tClusterSizes->GetNbinsY());

        if(dStripStatistics > dLargestStripStatistics)
        {
          iCounterMainStripBin = iStrip + 1;
          dLargestStripStatistics = dStripStatistics;
        }
      }
    }

    TH1D* tCentralClusterSizeProjection = tClusterSizes->ProjectionY(Form("tCentralClusterSizeProjection_%s", tCounterTag.Data()),
                                                                     iCounterMainStripBin, iCounterMainStripBin);

    TProfile* tStripToTSpectraProfile = tStripToTSpectra->ProfileX(Form("tStripToTSpectraProfile_%s", tCounterTag.Data()),
                                                                   1, -1, "s");
    TH1D* tCentralStripToTSpectrum = tStripToTSpectra->ProjectionY(Form("tCentralStripToTSpectrum_%s", tCounterTag.Data()),
                                                                   2*iCounterMainStripBin, 2*iCounterMainStripBin);


    TProfile* tStripPositionSpectraProfile = tStripPositionSpectra->ProfileX(Form("tStripPositionSpectraProfile_%s", tCounterTag.Data()),
                                                                             1, -1, "s");
    TH1D* tStripPositionSpectraProjection = tStripPositionSpectra->ProjectionY(Form("tStripPositionSpectraProjection_%s", tCounterTag.Data()), 1, tStripPositionSpectra->GetNbinsX());
    tStripPositionSpectraProjection->Scale(1./tClusterSizes->GetMean(2));
    tStripPositionSpectraProjection->SetLineColor(kGreen);

    TProfile* tPointPositionSpectraProfile = tPointPositionSpectra->ProfileX(Form("tPointPositionSpectraProfile_%s", tCounterTag.Data()),
                                                                             1, -1, "s");
    TH1D* tPointPositionSpectraProjection = tPointPositionSpectra->ProjectionY(Form("tPointPositionSpectraProjection_%s", tCounterTag.Data()), 1, tPointPositionSpectra->GetNbinsX());
    tPointPositionSpectraProjection->SetLineColor(kBlack);

    TProfile* tHitPositionSpectraProfile = tHitPositionSpectra->ProfileX(Form("tHitPositionSpectraProfile_%s", tCounterTag.Data()),
                                                                         1, -1, "s");
    TH1D* tHitPositionSpectraProjection = tHitPositionSpectra->ProjectionY(Form("tHitPositionSpectraProjection_%s", tCounterTag.Data()), 1, tHitPositionSpectra->GetNbinsX());
    tHitPositionSpectraProjection->SetLineColor(kRed);

    cout<<Form("%d-%d-%d: %5.1f %% of hits reconstructed on the readout cell.", giModuleType, giModuleIndex, giCounterIndex,
               tHitPositionSpectraProjection->Integral(tHitPositionSpectraProjection->GetXaxis()->FindBin(-gdCellAlongWidth/2.),
                                                       tHitPositionSpectraProjection->GetXaxis()->FindBin(gdCellAlongWidth/2.))
               /tHitPositionSpectraProjection->Integral()*100.)<<endl;

    tStripTimeResidualWalk->GetXaxis()->SetRange(2*iCounterMainStripBin, 2*iCounterMainStripBin);
    TH2D* tCentralStripTopTimeResidualWalk = dynamic_cast<TH2D*>(tStripTimeResidualWalk->Project3D(Form("%s_zy", tCounterTag.Data())));
    tCentralStripTopTimeResidualWalk->SetTitle("");


    TH1F* tStripToTSpectraMean = new TH1F(Form("tStripToTSpectraMean_%s", tCounterTag.Data()), "; ToT mean values [ns]", 50, 0., 10.);

    TH1F* tStripToTSpectraRMS = new TH1F(Form("tStripToTSpectraRMS_%s", tCounterTag.Data()), "; ToT RMS values [ns]", 20, 0., 2.);

    for(Int_t iBin = 1; iBin <= tStripToTSpectraProfile->GetNbinsX(); iBin++)
    {
      tStripToTSpectraMean->Fill(tStripToTSpectraProfile->GetBinContent(iBin));
      tStripToTSpectraRMS->Fill(tStripToTSpectraProfile->GetBinError(iBin));
    }


    TEfficiency* tHitChiSquareAcc = new TEfficiency(Form("tHitChiSquareAcc_%s", tCounterTag.Data()),
                                                    "; #chi^{2}_{3} limit []; #chi^{2}_{3} acceptance []",
                                                    tHitChiSquare->GetNbinsX(), tHitChiSquare->GetXaxis()->GetXmin(), tHitChiSquare->GetXaxis()->GetXmax());

    TEfficiency* tCorrectedHitChiSquareAcc = new TEfficiency(Form("tCorrectedHitChiSquareAcc_%s", tCounterTag.Data()),
                                                             "; #chi^{2}_{3} limit []; #chi^{2}_{3} acceptance []",
                                                             tCorrectedHitChiSquare->GetNbinsX(), tCorrectedHitChiSquare->GetXaxis()->GetXmin(), tCorrectedHitChiSquare->GetXaxis()->GetXmax());

    TVectorD dChiSquareQuantilesTheory(tHitChiSquare->GetNbinsX());
    TVectorD dChiSquareQuantiles(tHitChiSquare->GetNbinsX());
    TVectorD dCorrectedChiSquareQuantilesTheory(tHitChiSquare->GetNbinsX());
    TVectorD dCorrectedChiSquareQuantiles(tHitChiSquare->GetNbinsX());

    for(Int_t iBin = 1; iBin <= tHitChiSquare->GetNbinsX(); iBin++)
    {
      Double_t dChiSqLim = tHitChiSquareAcc->GetTotalHistogram()->GetXaxis()->GetBinUpEdge(iBin);

      tHitChiSquareAcc->SetTotalEvents(iBin + 1, tHitChiSquare->GetEntries());

      Int_t iPassedIntegral = static_cast<Int_t>(tHitChiSquare->Integral(1, tHitChiSquare->FindBin(dChiSqLim) - 1));
      tHitChiSquareAcc->SetPassedEvents(iBin + 1, iPassedIntegral);

      dChiSquareQuantilesTheory[iBin - 1] = ROOT::MathMore::chisquared_quantile(tHitChiSquareAcc->GetEfficiency(iBin), 3.);
      dChiSquareQuantiles[iBin - 1] = dChiSqLim;

      tCorrectedHitChiSquareAcc->SetTotalEvents(iBin + 1, tCorrectedHitChiSquare->GetEntries());

      iPassedIntegral = static_cast<Int_t>(tCorrectedHitChiSquare->Integral(1, tCorrectedHitChiSquare->FindBin(dChiSqLim) - 1));
      tCorrectedHitChiSquareAcc->SetPassedEvents(iBin + 1, iPassedIntegral);

      dCorrectedChiSquareQuantilesTheory[iBin - 1] = ROOT::MathMore::chisquared_quantile(tCorrectedHitChiSquareAcc->GetEfficiency(iBin), 3.);
      dCorrectedChiSquareQuantiles[iBin - 1] = dChiSqLim;
    }


    TH1F* tAcrossPositionResidualEvolutionSigma = new TH1F(Form("tAcrossPositionResidualEvolutionSigma_%s", tCounterTag.Data()),
                                                           Form("; %s; #perp residual sigma [cm]", tAcrossPositionResidualEvolution->GetXaxis()->GetTitle()),
                                                           tAcrossPositionResidualEvolution->GetNbinsX(),
                                                           tAcrossPositionResidualEvolution->GetXaxis()->GetXmin(), tAcrossPositionResidualEvolution->GetXaxis()->GetXmax());

    TH1F* tAlongPositionResidualEvolutionSigma = new TH1F(Form("tAlongPositionResidualEvolutionSigma_%s", tCounterTag.Data()),
                                                          Form("; %s; #parallel residual sigma [cm]", tAlongPositionResidualEvolution->GetXaxis()->GetTitle()),
                                                          tAlongPositionResidualEvolution->GetNbinsX(),
                                                          tAlongPositionResidualEvolution->GetXaxis()->GetXmin(), tAlongPositionResidualEvolution->GetXaxis()->GetXmax());

    TH1F* tTimeResidualEvolutionSigma = new TH1F(Form("tTimeResidualEvolutionSigma_%s", tCounterTag.Data()),
                                                 Form("; %s; residual sigma [ns]", tTimeResidualEvolution->GetXaxis()->GetTitle()),
                                                 tTimeResidualEvolution->GetNbinsX(), tTimeResidualEvolution->GetXaxis()->GetXmin(), tTimeResidualEvolution->GetXaxis()->GetXmax());

    for(Int_t iBin = 1; iBin <= tTimeResidualEvolution->GetNbinsX(); iBin++)
    {
      TH1D* tProjection = tAcrossPositionResidualEvolution->ProjectionY("_py", iBin, iBin);

      if(tProjection->GetEntries())
      {
        tAcrossPositionResidualEvolutionSigma->SetBinContent(iBin, tProjection->GetRMS());
        tAcrossPositionResidualEvolutionSigma->SetBinError(iBin, tProjection->GetRMSError());
      }

      delete tProjection;

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

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

      delete tProjection;

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

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

      delete tProjection;
    }


    TH1D* tCentralClusterSizeProjectionSPH = tClusterSizesSPH->ProjectionY(Form("tCentralClusterSizeProjectionSPH_%s", tCounterTag.Data()),
                                                                           iCounterMainStripBin, iCounterMainStripBin);

    TProfile* tStripToTSpectraProfileSPH = tStripToTSpectraSPH->ProfileX(Form("tStripToTSpectraProfileSPH_%s", tCounterTag.Data()),
                                                                         1, -1, "s");
    TH1D* tCentralStripToTSpectrumSPH = tStripToTSpectraSPH->ProjectionY(Form("tCentralStripToTSpectrumSPH_%s", tCounterTag.Data()),
                                                                         2*iCounterMainStripBin, 2*iCounterMainStripBin);


    TProfile* tStripPositionSpectraProfileSPH = tStripPositionSpectraSPH->ProfileX(Form("tStripPositionSpectraProfileSPH_%s", tCounterTag.Data()),
                                                                                   1, -1, "s");
    TH1D* tStripPositionSpectraProjectionSPH = tStripPositionSpectraSPH->ProjectionY(Form("tStripPositionSpectraProjectionSPH_%s", tCounterTag.Data()), 1, tStripPositionSpectra->GetNbinsX());
    tStripPositionSpectraProjectionSPH->Scale(1./tClusterSizesSPH->GetMean(2));
    tStripPositionSpectraProjectionSPH->SetLineColor(kGreen);

    TProfile* tHitPositionSpectraProfileSPH = tHitPositionSpectraSPH->ProfileX(Form("tHitPositionSpectraProfileSPH_%s", tCounterTag.Data()),
                                                                               1, -1, "s");
    TH1D* tHitPositionSpectraProjectionSPH = tHitPositionSpectraSPH->ProjectionY(Form("tHitPositionSpectraProjectionSPH_%s", tCounterTag.Data()), 1, tHitPositionSpectra->GetNbinsX());
    tHitPositionSpectraProjectionSPH->SetLineColor(kRed);

    tStripTimeResidualWalkSPH->GetXaxis()->SetRange(2*iCounterMainStripBin, 2*iCounterMainStripBin);
    TH2D* tCentralStripTopTimeResidualWalkSPH = dynamic_cast<TH2D*>(tStripTimeResidualWalkSPH->Project3D(Form("SPH%s_zy", tCounterTag.Data())));
    tCentralStripTopTimeResidualWalkSPH->SetTitle("");


    TEfficiency* tHitChiSquareAccSPH = new TEfficiency(Form("tHitChiSquareAccSPH_%s", tCounterTag.Data()),
                                                       "; #chi^{2}_{3} limit []; #chi^{2}_{3} acceptance []",
                                                       tHitChiSquareSPH->GetNbinsX(), tHitChiSquareSPH->GetXaxis()->GetXmin(), tHitChiSquareSPH->GetXaxis()->GetXmax());

    TEfficiency* tCorrectedHitChiSquareAccSPH = new TEfficiency(Form("tCorrectedHitChiSquareAccSPH_%s", tCounterTag.Data()),
                                                                "; #chi^{2}_{3} limit []; #chi^{2}_{3} acceptance []",
                                                                tCorrectedHitChiSquareSPH->GetNbinsX(), tCorrectedHitChiSquareSPH->GetXaxis()->GetXmin(), tCorrectedHitChiSquareSPH->GetXaxis()->GetXmax());

    TVectorD dChiSquareQuantilesTheorySPH(tHitChiSquareSPH->GetNbinsX());
    TVectorD dChiSquareQuantilesSPH(tHitChiSquareSPH->GetNbinsX());
    TVectorD dCorrectedChiSquareQuantilesTheorySPH(tHitChiSquareSPH->GetNbinsX());
    TVectorD dCorrectedChiSquareQuantilesSPH(tHitChiSquareSPH->GetNbinsX());

    for(Int_t iBin = 1; iBin <= tHitChiSquareSPH->GetNbinsX(); iBin++)
    {
      Double_t dChiSqLim = tHitChiSquareAccSPH->GetTotalHistogram()->GetXaxis()->GetBinUpEdge(iBin);

      tHitChiSquareAccSPH->SetTotalEvents(iBin + 1, tHitChiSquareSPH->GetEntries());

      Int_t iPassedIntegral = static_cast<Int_t>(tHitChiSquareSPH->Integral(1, tHitChiSquareSPH->FindBin(dChiSqLim) - 1));
      tHitChiSquareAccSPH->SetPassedEvents(iBin + 1, iPassedIntegral);

      dChiSquareQuantilesTheorySPH[iBin - 1] = ROOT::MathMore::chisquared_quantile(tHitChiSquareAccSPH->GetEfficiency(iBin), 3.);
      dChiSquareQuantilesSPH[iBin - 1] = dChiSqLim;

      tCorrectedHitChiSquareAccSPH->SetTotalEvents(iBin + 1, tCorrectedHitChiSquareSPH->GetEntries());

      iPassedIntegral = static_cast<Int_t>(tCorrectedHitChiSquareSPH->Integral(1, tCorrectedHitChiSquareSPH->FindBin(dChiSqLim) - 1));
      tCorrectedHitChiSquareAccSPH->SetPassedEvents(iBin + 1, iPassedIntegral);

      dCorrectedChiSquareQuantilesTheorySPH[iBin - 1] = ROOT::MathMore::chisquared_quantile(tCorrectedHitChiSquareAccSPH->GetEfficiency(iBin), 3.);
      dCorrectedChiSquareQuantilesSPH[iBin - 1] = dChiSqLim;
    }


    TH1F* tAcrossPositionResidualEvolutionSigmaSPH = new TH1F(Form("tAcrossPositionResidualEvolutionSigmaSPH_%s", tCounterTag.Data()),
                                                              Form("; %s; #perp residual sigma [cm]", tAcrossPositionResidualEvolutionSPH->GetXaxis()->GetTitle()),
                                                              tAcrossPositionResidualEvolutionSPH->GetNbinsX(),
                                                              tAcrossPositionResidualEvolutionSPH->GetXaxis()->GetXmin(), tAcrossPositionResidualEvolutionSPH->GetXaxis()->GetXmax());

    TH1F* tAlongPositionResidualEvolutionSigmaSPH = new TH1F(Form("tAlongPositionResidualEvolutionSigmaSPH_%s", tCounterTag.Data()),
                                                             Form("; %s; #parallel residual sigma [cm]", tAlongPositionResidualEvolutionSPH->GetXaxis()->GetTitle()),
                                                             tAlongPositionResidualEvolutionSPH->GetNbinsX(),
                                                             tAlongPositionResidualEvolutionSPH->GetXaxis()->GetXmin(), tAlongPositionResidualEvolutionSPH->GetXaxis()->GetXmax());

    TH1F* tTimeResidualEvolutionSigmaSPH = new TH1F(Form("tTimeResidualEvolutionSigmaSPH_%s", tCounterTag.Data()),
                                                    Form("; %s; residual sigma [ns]", tTimeResidualEvolutionSPH->GetXaxis()->GetTitle()),
                                                    tTimeResidualEvolutionSPH->GetNbinsX(), tTimeResidualEvolutionSPH->GetXaxis()->GetXmin(), tTimeResidualEvolutionSPH->GetXaxis()->GetXmax());

    for(Int_t iBin = 1; iBin <= tTimeResidualEvolutionSPH->GetNbinsX(); iBin++)
    {
      TH1D* tProjection = tAcrossPositionResidualEvolutionSPH->ProjectionY("_py", iBin, iBin);

      if(tProjection->GetEntries())
      {
        tAcrossPositionResidualEvolutionSigmaSPH->SetBinContent(iBin, tProjection->GetRMS());
        tAcrossPositionResidualEvolutionSigmaSPH->SetBinError(iBin, tProjection->GetRMSError());
      }

      delete tProjection;

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

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

      delete tProjection;

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

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

      delete tProjection;
    }


    TH1F* tTimeResidualToLastHitUnbiasedSigma = new TH1F(Form("tTimeResidualToLastHitUnbiasedSigma_%s", tCounterTag.Data()),
                                                         Form("; %s; residual sigma [ns]", tTimeResidualToLastHitUnbiased->GetXaxis()->GetTitle()),
                                                         tTimeResidualToLastHitUnbiased->GetNbinsX(), tTimeResidualToLastHitUnbiased->GetXaxis()->GetXbins()->GetArray());

    for(Int_t iBin = 1; iBin <= tTimeResidualToLastHitUnbiased->GetNbinsX(); iBin++)
    {
      TH1D* tProjection = tTimeResidualToLastHitUnbiased->ProjectionY("_py", iBin, iBin);

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

      delete tProjection;
    }


    TH1F* tTimeResidualToLastHitBiasedSigma = new TH1F(Form("tTimeResidualToLastHitBiasedSigma_%s", tCounterTag.Data()),
                                                       Form("; %s; residual sigma [ns]", tTimeResidualToLastHitBiased->GetXaxis()->GetTitle()),
                                                       tTimeResidualToLastHitBiased->GetNbinsX(), tTimeResidualToLastHitBiased->GetXaxis()->GetXbins()->GetArray());

    for(Int_t iBin = 1; iBin <= tTimeResidualToLastHitBiased->GetNbinsX(); iBin++)
    {
      TH1D* tProjection = tTimeResidualToLastHitBiased->ProjectionY("_py", iBin, iBin);

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

      delete tProjection;
    }


    TH3D* tLocalSurfaceFluxTotal = dynamic_cast<TH3D*>(tLocalSurfaceFlux->GetCopyTotalHisto());
    TH3D* tLocalSurfaceFluxPassed = dynamic_cast<TH3D*>(tLocalSurfaceFlux->GetCopyPassedHisto());
/*
    TH2D* tMCPointPositionSpectrum = dynamic_cast<TH2D*>(tLocalSurfaceFluxTotal->Project3D("_yxNUFNOF"));
    tMCPointPositionSpectrum->SetName(Form("tMCPointPositionSpectrum_%s", tCounterTag.Data()));
    tMCPointPositionSpectrum->SetTitle("");
    TH1D* tMCPointPositionSpectrumProjection = dynamic_cast<TH1D*>(tLocalSurfaceFluxTotal->Project3D("_yNUFNOF"));
    tMCPointPositionSpectrumProjection->SetName(Form("tMCPointPositionSpectrumProjection_%s", tCounterTag.Data()));
    tMCPointPositionSpectrumProjection->SetTitle("");
    tMCPointPositionSpectrumProjection->SetLineColor(kBlack);
    TProfile* tMCPointPositionSpectrumProfile = tMCPointPositionSpectrum->ProfileX(Form("tMCPointPositionSpectrumProfile_%s", tCounterTag.Data()), 1, -1, "s");
    tMCPointPositionSpectrumProfile->SetTitle("");

    TH2D* tHitPositionSpectrum = dynamic_cast<TH2D*>(tLocalSurfaceFluxPassed->Project3D("_yxNUFNOF"));
    tHitPositionSpectrum->SetName(Form("tHitPositionSpectrum_%s", tCounterTag.Data()));
    tHitPositionSpectrum->SetTitle("");
    TH1D* tHitPositionSpectrumProjection = dynamic_cast<TH1D*>(tLocalSurfaceFluxPassed->Project3D("_yNUFNOF"));
    tHitPositionSpectrumProjection->SetName(Form("tHitPositionSpectrumProjection_%s", tCounterTag.Data()));
    tHitPositionSpectrumProjection->SetTitle("");
    tHitPositionSpectrumProjection->SetLineColor(kRed);
    TProfile* tHitPositionSpectrumProfile = tHitPositionSpectrum->ProfileX(Form("tHitPositionSpectrumProfile_%s", tCounterTag.Data()), 1, -1, "s");
    tHitPositionSpectrumProfile->SetTitle("");
*/

/*
    TCanvas* tCanvas = new TCanvas(Form("tCanvas_global_%s", tCounterTag.Data()),Form("Global response of MRPC %s", tCounterTag.Data()),0,0,1750,1050);
    tCanvas->Divide(5,3);

    TCanvas* tCanvas = new TCanvas(Form("tCanvas_global_%s", tCounterTag.Data()),Form("Global response of MRPC %s", tCounterTag.Data()),0,0,1750,1400);
    tCanvas->Divide(5,4);
*/
    TCanvas* tCanvas = new TCanvas(Form("tCanvas_global_int_%s", tCounterTag.Data()), Form("Global integral response of MRPC %s", tCounterTag.Data()), 0, 0, 1800, 900);
    tCanvas->Divide(6, 3);

    tCanvas->cd(1);
    gPad->SetLogz();
    tClusterSizes->Draw("colz");
/*
    tCanvas->cd(3);
    tStripToTSpectraMean->Draw();

    tCanvas->cd(4);
    tStripToTSpectraRMS->Draw();
*/
    tCanvas->cd(2);
    gPad->SetLogy();
    tCentralClusterSizeProjection->Draw();

    tCanvas->cd(3);
    gPad->SetLogz();
    tStripToTSpectra->Draw("colz");
    tStripToTSpectraProfile->Draw("same");

    tCanvas->cd(4);
    gPad->SetLogy();
    tCentralStripToTSpectrum->Draw();

    tCanvas->cd(5);
    gPad->SetLogz();
    tStripPositionSpectra->Draw("colz");
    tStripPositionSpectraProfile->Draw("same");

    tCanvas->cd(6);
    gPad->SetLogz();
    tCentralStripTopTimeResidualWalk ->Draw("colz");
    tCentralStripTopTimeResidualWalk ->ProfileX()->Draw("sameE1");

    tCanvas->cd(7);
    gPad->SetLogy();
    tHitAcrossPositionResidual->Draw();
    tHitAcrossPositionResidual->Fit("gaus","Q");
    gPad->Update();
    TPaveStats* tStatsBoxHitAcrossPositionResidual = dynamic_cast<TPaveStats*>(tHitAcrossPositionResidual->FindObject("stats"));
    tStatsBoxHitAcrossPositionResidual->SetOptFit(10001);
    gPad->Modified();

    tCanvas->cd(8);
    gPad->SetLogy();
    tHitAlongPositionResidual->Draw();
    tHitAlongPositionResidual->Fit("gaus","Q");
    gPad->Update();
    TPaveStats* tStatsBoxHitAlongPositionResidual = dynamic_cast<TPaveStats*>(tHitAlongPositionResidual->FindObject("stats"));
    tStatsBoxHitAlongPositionResidual->SetOptFit(10001);
    gPad->Modified();

    tCanvas->cd(9);
    gPad->SetLogy();
    tHitTimeResidual->Fit("gaus","Q0");
    const Int_t kNotDraw = 1 << 9;
    TF1* tFitFunction = tHitTimeResidual->GetFunction("gaus");
    if(tFitFunction)
    {
      tFitFunction->ResetBit(kNotDraw);
    }
    tHitTimeResidual->Draw();
    gPad->Update();
    TPaveStats* tStatsBoxHitTimeResidual = dynamic_cast<TPaveStats*>(tHitTimeResidual->FindObject("stats"));
    tStatsBoxHitTimeResidual->SetOptFit(10001);
    gPad->Modified();

    tCanvas->cd(10);
    TF1* tChiSquareTheory = new TF1(Form("tChiSquareTheory_%s", tCounterTag.Data()), chisqpdf, 0., 30., 2);
    tChiSquareTheory->SetParameter(0, 3);
    tChiSquareTheory->SetParameter(1, tHitChiSquare->Integral("width"));
    tChiSquareTheory->SetNpx(1000000);
    tChiSquareTheory->SetLineColor(kRed);
    tChiSquareTheory->SetLineWidth(2);
    gPad->SetLogy();
    tCorrectedHitChiSquare->SetLineColor(kGreen + 2);
    tCorrectedHitChiSquare->SetLineWidth(2);
    tCorrectedHitChiSquare->Draw();
    tHitChiSquare->SetLineColor(kBlue);
    tHitChiSquare->SetLineWidth(2);
    tHitChiSquare->Draw("SAME");
    tChiSquareTheory->Draw("SAME");

    tCanvas->cd(11);
    TF1* tChiSquareAccTheory = new TF1(Form("tChiSquareAccTheory_%s", tCounterTag.Data()), chisqcdf, 0., 30., 1);
    tChiSquareAccTheory->SetParameter(0, 3);
    tChiSquareAccTheory->SetNpx(1000000);
    tChiSquareAccTheory->SetLineColor(kRed);
    tChiSquareAccTheory->SetLineWidth(2);
    gPad->SetGridx();
    gPad->SetGridy();
    tHitChiSquareAcc->SetLineColor(kBlue);
    tHitChiSquareAcc->SetMarkerColor(kBlue);
//    tHitChiSquareAcc->SetMarkerStyle(20);
//    tHitChiSquareAcc->SetMarkerSize(0.5);
    tHitChiSquareAcc->Draw("AP");
    gPad->Update();
    gPad->Modified();
    tHitChiSquareAcc->GetPaintedGraph()->GetYaxis()->SetRangeUser(0.5, 1.);
    tHitChiSquareAcc->GetPaintedGraph()->GetXaxis()->SetLimits(0., 30.);
    gPad->Modified();
    tCorrectedHitChiSquareAcc->SetLineColor(kGreen + 2);
    tCorrectedHitChiSquareAcc->SetMarkerColor(kGreen + 2);
//    tCorrectedHitChiSquareAcc->SetMarkerStyle(20);
//    tCorrectedHitChiSquareAcc->SetMarkerSize(0.5);
    tCorrectedHitChiSquareAcc->Draw("P SAME");
    tChiSquareAccTheory->Draw("SAME");

    tCanvas->cd(12);
    TGraph* tChiSquareQQPlot = new TGraph(dChiSquareQuantilesTheory, dChiSquareQuantiles);
    tChiSquareQQPlot->SetTitle("");
    tChiSquareQQPlot->GetXaxis()->SetTitle("#chi^{2}_{3} quantiles []");
    tChiSquareQQPlot->GetYaxis()->SetTitle("cluster quantiles []");
    tChiSquareQQPlot->SetMarkerStyle(21);
    tChiSquareQQPlot->SetMarkerSize(0.5);
    tChiSquareQQPlot->SetMarkerColor(kBlue);
    tChiSquareQQPlot->Draw("AP");
    gPad->Update();
    tChiSquareQQPlot->SetMinimum(0.);
    tChiSquareQQPlot->SetMaximum(30.);
    tChiSquareQQPlot->GetXaxis()->SetLimits(0., 30.);
    gPad->Update();
    gPad->Modified();
    TGraph* tCorrectedChiSquareQQPlot = new TGraph(dCorrectedChiSquareQuantilesTheory, dCorrectedChiSquareQuantiles);
    tCorrectedChiSquareQQPlot->SetTitle("");
    tCorrectedChiSquareQQPlot->GetXaxis()->SetTitle("#chi^{2}_{3} quantiles []");
    tCorrectedChiSquareQQPlot->GetYaxis()->SetTitle("cluster quantiles []");
    tCorrectedChiSquareQQPlot->SetMarkerStyle(21);
    tCorrectedChiSquareQQPlot->SetMarkerSize(0.5);
    tCorrectedChiSquareQQPlot->SetMarkerColor(kGreen + 2);
    tCorrectedChiSquareQQPlot->Draw("P SAME");
    TLine* tLine = new TLine(0., 0., 30., 30.);
    tLine->SetLineColor(kRed);
    tLine->SetLineWidth(2);
    tLine->SetLineStyle(1);
    tLine->Draw("SAME");

    tCanvas->cd(13);
    tPointPositionSpectra->Draw("colz");
    tPointPositionSpectraProfile->Draw("SAME");

    tCanvas->cd(14);
    tHitPositionSpectra->Draw("colz");
    tHitPositionSpectraProfile->Draw("SAME");

    tCanvas->cd(15);
    tPointPositionSpectraProjection->Draw("HIST");
    tHitPositionSpectraProjection->Draw("HIST SAME");
    tStripPositionSpectraProjection->Draw("SAME");

    tCanvas->cd(16);
    Double_t dStripResolution = 1.;
    TF1* tResidualFunction = tHitAlongPositionResidual->GetFunction("gaus");
    if(tResidualFunction)
    {
      dStripResolution = tResidualFunction->GetParameter(2);
    }
    TF1* tFunctionGPB = new TF1(Form("tFunctionGPB_%s", tCounterTag.Data()), gpb, -2.*gdCellAlongWidth, 2.*gdCellAlongWidth, 4);
    tFunctionGPB->SetParameter(0, 0.);
    tFunctionGPB->SetParameter(1, dStripResolution);
    if(bSpotResponse)
    {
      tFunctionGPB->SetParameter(2, 1.);
    }
    else
    {
      tFunctionGPB->SetParameter(2, gdCellAlongWidth);
    }
    tFunctionGPB->SetParameter(3, tHitPositionSpectraProjection->Integral("width"));
    tFunctionGPB->SetNpx(1000000);
    tFunctionGPB->SetLineColor(kBlue);
    tFunctionGPB->SetLineStyle(9);
    tHitPositionSpectraProjection->Draw();
    TF1* tFunctionGaus(NULL);
    if(5 == giModuleType && !bSpotResponse)
    {
      tFunctionGaus = new TF1(*dynamic_cast<TF1*>(gROOT->GetFunction("gaus")));
      tFunctionGaus->SetParameter(0, tHitPositionSpectraProjection->Integral("width"));
      tFunctionGaus->SetParameter(1, 0.);
      tFunctionGaus->SetParameter(2, TMath::Sqrt(TMath::Power(gdCellAlongWidth/10., 2.) + TMath::Power(dStripResolution, 2.)));
      tFunctionGaus->SetNpx(1000000);
      tFunctionGaus->SetLineColor(kBlue);
      tFunctionGaus->SetLineStyle(9);
      tFunctionGaus->Draw("SAME");
    }
    else
    {
      tFunctionGPB->Draw("SAME");
    }

    tCanvas->cd(17);
    gPad->SetLogy();
    tNPointsPerHit->Draw();

    tCanvas->cd(18);
    gPad->SetLogy();
    tNHitsPerPoint->Draw();

    if(bBatchMode)
    {
      tCanvas->SaveAs(Form("pl_global_integral_response_%s.pdf", tCounterTag.Data()));

      delete tCanvas;
    }


    TCanvas* tCanvas6 = new TCanvas(Form("tCanvas_global_diff_%s", tCounterTag.Data()), Form("Global differential response of MRPC %s", tCounterTag.Data()), 0, 0, 1800, 900);
    tCanvas6->Divide(6, 3);

    tCanvas6->cd(1);
    tIntegratedEfficiency->SetLineColor(kBlack);
    tIntegratedEfficiency->SetMarkerColor(kBlack);
    tIntegratedEfficiency->SetMarkerStyle(20);
    tIntegratedEfficiency->SetMarkerSize(0.5);
    tIntegratedEfficiency->Draw("AP");
    gPad->Update();
    tIntegratedEfficiency->GetPaintedGraph()->SetMinimum(0.);
    tIntegratedEfficiency->GetPaintedGraph()->SetMaximum(1.1);
    tIntegratedEfficiency->GetPaintedGraph()->GetXaxis()->SetLimits(tIntegratedEfficiency->GetTotalHistogram()->GetXaxis()->GetXmin(),
                                                                    tIntegratedEfficiency->GetTotalHistogram()->GetXaxis()->GetXmax());
    gPad->Modified();

    cout<<Form("       %5.1f %% mean efficiency", 100.*tIntegratedEfficiency->GetPassedHistogram()->Integral()/tIntegratedEfficiency->GetTotalHistogram()->Integral()) << endl;

    tCanvas6->cd(2);
    tIntegratedEfficiency->GetPaintedGraph()->Draw("AP");
    gPad->Update();
    gPad->Modified();
    tCorrectedChiSquareEfficiency->SetLineColor(kGreen + 2);
    tCorrectedChiSquareEfficiency->SetMarkerColor(kGreen + 2);
    tCorrectedChiSquareEfficiency->SetMarkerStyle(20);
    tCorrectedChiSquareEfficiency->SetMarkerSize(0.5);
    tCorrectedChiSquareEfficiency->Draw("P SAME");
    tChiSquareEfficiency->SetLineColor(kBlue);
    tChiSquareEfficiency->SetMarkerColor(kBlue);
    tChiSquareEfficiency->SetMarkerStyle(20);
    tChiSquareEfficiency->SetMarkerSize(0.5);
    tChiSquareEfficiency->Draw("P SAME");

    tCanvas6->cd(3);
    gPad->SetLogz();
    tClusterSizeEvolution->Draw("colz");

    tCanvas6->cd(4);
    TProfile* tClusterSizeEvolutionProfile = tClusterSizeEvolution->ProfileX();
    tClusterSizeEvolutionProfile->GetYaxis()->SetTitle("mean cluster size [strips]");
    tClusterSizeEvolutionProfile->Draw("PE1");
    tClusterSizeEvolutionProfile->GetYaxis()->SetRangeUser(1., 3.);

    tCanvas6->cd(5);
    gPad->SetLogz();
    tCentralStripToTEvolution->Draw("colz");

    tCanvas6->cd(6);
    TProfile* tCentralStripToTEvolutionProfile = tCentralStripToTEvolution->ProfileX();
    tCentralStripToTEvolutionProfile->GetYaxis()->SetTitle("central ToT mean [ns]");
    tCentralStripToTEvolutionProfile->Draw("PE1");
    tCentralStripToTEvolutionProfile->GetYaxis()->SetRangeUser(1., 3.);

    tCanvas6->cd(7);
    gPad->SetLogz();
    tAcrossPositionResidualEvolution->Draw("colz");

    tCanvas6->cd(8);
    tAcrossPositionResidualEvolutionSigma->Draw("PE1");
    tAcrossPositionResidualEvolutionSigma->GetYaxis()->SetRangeUser(0., 0.5);

    tCanvas6->cd(9);
    gPad->SetLogz();
    tAlongPositionResidualEvolution->Draw("colz");

    tCanvas6->cd(10);
    tAlongPositionResidualEvolutionSigma->Draw("PE1");
    tAlongPositionResidualEvolutionSigma->GetYaxis()->SetRangeUser(0., 3.);

    tCanvas6->cd(11);
    gPad->SetLogz();
    tTimeResidualEvolution->Draw("colz");
//    tTimeResidualEvolution->ProfileX("_pfx", 1, -1, "s")->Draw("sameE1");

    tCanvas6->cd(12);
    tTimeResidualEvolutionSigma->Draw("PE1");
    tTimeResidualEvolutionSigma->GetYaxis()->SetRangeUser(0., 0.2);


    if(bBatchMode)
    {
      tCanvas6->SaveAs(Form("pl_global_differential_response_%s.pdf", tCounterTag.Data()));

      delete tCanvas6;
    }


    TCanvas* tCanvas1 = new TCanvas(Form("tCanvas_global_int_sph_%s", tCounterTag.Data()), Form("Global integral SPH response of MRPC %s", tCounterTag.Data()), 0, 0, 1800, 900);
    tCanvas1->Divide(6, 3);

    tCanvas1->cd(1);
    gPad->SetLogz();
    tClusterSizesSPH->Draw("colz");

    tCanvas1->cd(2);
    gPad->SetLogy();
    tCentralClusterSizeProjectionSPH->Draw();

    tCanvas1->cd(3);
    gPad->SetLogz();
    tStripToTSpectraSPH->Draw("colz");
    tStripToTSpectraProfileSPH->Draw("same");

    tCanvas1->cd(4);
    gPad->SetLogy();
    tCentralStripToTSpectrumSPH->Draw();

    tCanvas1->cd(5);
    gPad->SetLogz();
    tStripPositionSpectraSPH->Draw("colz");
    tStripPositionSpectraProfileSPH->Draw("same");

    tCanvas1->cd(6);
    gPad->SetLogz();
    tCentralStripTopTimeResidualWalkSPH->Draw("colz");
    tCentralStripTopTimeResidualWalkSPH->ProfileX()->Draw("sameE1");

    tCanvas1->cd(7);
    gPad->SetLogy();
    tHitAcrossPositionResidualSPH->Draw();
    tHitAcrossPositionResidualSPH->Fit("gaus","Q");
    gPad->Update();
    TPaveStats* tStatsBoxHitAcrossPositionResidualSPH = dynamic_cast<TPaveStats*>(tHitAcrossPositionResidualSPH->FindObject("stats"));
    tStatsBoxHitAcrossPositionResidualSPH->SetOptFit(10001);
    gPad->Modified();

    tCanvas1->cd(8);
    gPad->SetLogy();
    tHitAlongPositionResidualSPH->Draw();
    tHitAlongPositionResidualSPH->Fit("gaus","Q");
    gPad->Update();
    TPaveStats* tStatsBoxHitAlongPositionResidualSPH = dynamic_cast<TPaveStats*>(tHitAlongPositionResidualSPH->FindObject("stats"));
    tStatsBoxHitAlongPositionResidualSPH->SetOptFit(10001);
    gPad->Modified();

    tCanvas1->cd(9);
    gPad->SetLogy();
    tHitTimeResidualSPH->Fit("gaus","Q0");
    tFitFunction = tHitTimeResidualSPH->GetFunction("gaus");
    if(tFitFunction)
    {
      tFitFunction->ResetBit(kNotDraw);
    }
    tHitTimeResidualSPH->Draw();
    gPad->Update();
    TPaveStats* tStatsBoxHitTimeResidualSPH = dynamic_cast<TPaveStats*>(tHitTimeResidualSPH->FindObject("stats"));
    tStatsBoxHitTimeResidualSPH->SetOptFit(10001);
    gPad->Modified();

    tCanvas1->cd(10);
    TF1* tChiSquareTheorySPH = new TF1(Form("tChiSquareTheorySPH_%s", tCounterTag.Data()), chisqpdf, 0., 30., 2);
    tChiSquareTheorySPH->SetParameter(0, 3);
    tChiSquareTheorySPH->SetParameter(1, tHitChiSquareSPH->Integral("width"));
    tChiSquareTheorySPH->SetNpx(1000000);
    tChiSquareTheorySPH->SetLineColor(kRed);
    tChiSquareTheorySPH->SetLineWidth(2);
    gPad->SetLogy();
    tCorrectedHitChiSquareSPH->SetLineColor(kGreen + 2);
    tCorrectedHitChiSquareSPH->SetLineWidth(2);
    tCorrectedHitChiSquareSPH->Draw();
    tHitChiSquareSPH->SetLineColor(kBlue);
    tHitChiSquareSPH->SetLineWidth(2);
    tHitChiSquareSPH->Draw("SAME");
    tChiSquareTheorySPH->Draw("SAME");

    tCanvas1->cd(11);
    TF1* tChiSquareAccTheorySPH = new TF1(Form("tChiSquareAccTheorySPH_%s", tCounterTag.Data()), chisqcdf, 0., 30., 1);
    tChiSquareAccTheorySPH->SetParameter(0, 3);
    tChiSquareAccTheorySPH->SetNpx(1000000);
    tChiSquareAccTheorySPH->SetLineColor(kRed);
    tChiSquareAccTheorySPH->SetLineWidth(2);
    gPad->SetGridx();
    gPad->SetGridy();
    tHitChiSquareAccSPH->SetLineColor(kBlue);
    tHitChiSquareAccSPH->SetMarkerColor(kBlue);
//    tHitChiSquareAccSPH->SetMarkerStyle(20);
//    tHitChiSquareAccSPH->SetMarkerSize(0.5);
    tHitChiSquareAccSPH->Draw("AP");
    gPad->Update();
    gPad->Modified();
    tHitChiSquareAccSPH->GetPaintedGraph()->GetYaxis()->SetRangeUser(0.5, 1.);
    tHitChiSquareAccSPH->GetPaintedGraph()->GetXaxis()->SetLimits(0., 30.);
    gPad->Modified();
    tCorrectedHitChiSquareAccSPH->SetLineColor(kGreen + 2);
    tCorrectedHitChiSquareAccSPH->SetMarkerColor(kGreen + 2);
//    tCorrectedHitChiSquareAccSPH->SetMarkerStyle(20);
//    tCorrectedHitChiSquareAccSPH->SetMarkerSize(0.5);
    tCorrectedHitChiSquareAccSPH->Draw("P SAME");
    tChiSquareAccTheorySPH->Draw("SAME");

    tCanvas1->cd(12);
    TGraph* tChiSquareQQPlotSPH = new TGraph(dChiSquareQuantilesTheorySPH, dChiSquareQuantilesSPH);
    tChiSquareQQPlotSPH->SetTitle("");
    tChiSquareQQPlotSPH->GetXaxis()->SetTitle("#chi^{2}_{3} quantiles []");
    tChiSquareQQPlotSPH->GetYaxis()->SetTitle("cluster quantiles []");
    tChiSquareQQPlotSPH->SetMarkerStyle(21);
    tChiSquareQQPlotSPH->SetMarkerSize(0.5);
    tChiSquareQQPlotSPH->SetMarkerColor(kBlue);
    tChiSquareQQPlotSPH->Draw("AP");
    gPad->Update();
    tChiSquareQQPlotSPH->SetMinimum(0.);
    tChiSquareQQPlotSPH->SetMaximum(30.);
    tChiSquareQQPlotSPH->GetXaxis()->SetLimits(0., 30.);
    gPad->Update();
    gPad->Modified();
    TGraph* tCorrectedChiSquareQQPlotSPH = new TGraph(dCorrectedChiSquareQuantilesTheorySPH, dCorrectedChiSquareQuantilesSPH);
    tCorrectedChiSquareQQPlotSPH->SetTitle("");
    tCorrectedChiSquareQQPlotSPH->GetXaxis()->SetTitle("#chi^{2}_{3} quantiles []");
    tCorrectedChiSquareQQPlotSPH->GetYaxis()->SetTitle("cluster quantiles []");
    tCorrectedChiSquareQQPlotSPH->SetMarkerStyle(21);
    tCorrectedChiSquareQQPlotSPH->SetMarkerSize(0.5);
    tCorrectedChiSquareQQPlotSPH->SetMarkerColor(kGreen + 2);
    tCorrectedChiSquareQQPlotSPH->Draw("P SAME");
    tLine = new TLine(0., 0., 30., 30.);
    tLine->SetLineColor(kRed);
    tLine->SetLineWidth(2);
    tLine->SetLineStyle(1);
    tLine->Draw("SAME");

    tCanvas1->cd(13);
    tPointPositionSpectra->Draw("colz");
    tPointPositionSpectraProfile->Draw("SAME");

    tCanvas1->cd(14);
    tHitPositionSpectraSPH->Draw("colz");
    tHitPositionSpectraProfileSPH->Draw("SAME");

    tCanvas1->cd(15);
    tPointPositionSpectraProjection->Draw("HIST");
    tHitPositionSpectraProjectionSPH->Draw("HIST SAME");
    tStripPositionSpectraProjectionSPH->Draw("SAME");

    tCanvas1->cd(16);
    dStripResolution = 1.;
    tResidualFunction = tHitAlongPositionResidualSPH->GetFunction("gaus");
    if(tResidualFunction)
    {
      dStripResolution = tResidualFunction->GetParameter(2);
    }
    tFunctionGPB->SetParameter(0, 0.);
    tFunctionGPB->SetParameter(1, dStripResolution);
    if(bSpotResponse)
    {
      tFunctionGPB->SetParameter(2, 1.);
    }
    else
    {
      tFunctionGPB->SetParameter(2, gdCellAlongWidth);
    }
    tFunctionGPB->SetParameter(3, tHitPositionSpectraProjectionSPH->Integral("width"));
    tFunctionGPB->SetNpx(1000000);
    tFunctionGPB->SetLineColor(kBlue);
    tFunctionGPB->SetLineStyle(9);
    tHitPositionSpectraProjectionSPH->Draw();
    if(5 == giModuleType && !bSpotResponse)
    {
      tFunctionGaus->SetParameter(0, tHitPositionSpectraProjectionSPH->Integral("width"));
      tFunctionGaus->SetParameter(1, 0.);
      tFunctionGaus->SetParameter(2, TMath::Sqrt(TMath::Power(gdCellAlongWidth/10., 2.) + TMath::Power(dStripResolution, 2.)));
      tFunctionGaus->SetNpx(1000000);
      tFunctionGaus->SetLineColor(kBlue);
      tFunctionGaus->SetLineStyle(9);
      tFunctionGaus->Draw("SAME");
    }
    else
    {
      tFunctionGPB->Draw("SAME");
    }

    if(bBatchMode)
    {
      tCanvas1->SaveAs(Form("pl_global_integral_sph_response_%s.pdf", tCounterTag.Data()));

      delete tCanvas1;
    }


    TCanvas* tCanvas7 = new TCanvas(Form("tCanvas_global_diff_sph_%s", tCounterTag.Data()), Form("Global differential SPH response of MRPC %s", tCounterTag.Data()), 0, 0, 1800, 900);
    tCanvas7->Divide(6, 3);

    tCanvas7->cd(3);
    gPad->SetLogz();
    tClusterSizeEvolutionSPH->Draw("colz");

    tCanvas7->cd(4);
    TProfile* tClusterSizeEvolutionProfileSPH = tClusterSizeEvolutionSPH->ProfileX();
    tClusterSizeEvolutionProfileSPH->GetYaxis()->SetTitle("mean cluster size [strips]");
    tClusterSizeEvolutionProfileSPH->Draw("PE1");
    tClusterSizeEvolutionProfileSPH->GetYaxis()->SetRangeUser(1., 3.);

    tCanvas7->cd(5);
    gPad->SetLogz();
    tCentralStripToTEvolutionSPH->Draw("colz");

    tCanvas7->cd(6);
    TProfile* tCentralStripToTEvolutionProfileSPH = tCentralStripToTEvolutionSPH->ProfileX();
    tCentralStripToTEvolutionProfileSPH->GetYaxis()->SetTitle("central ToT mean [ns]");
    tCentralStripToTEvolutionProfileSPH->Draw("PE1");
    tCentralStripToTEvolutionProfileSPH->GetYaxis()->SetRangeUser(1., 3.);

    tCanvas7->cd(7);
    gPad->SetLogz();
    tAcrossPositionResidualEvolutionSPH->Draw("colz");

    tCanvas7->cd(8);
    tAcrossPositionResidualEvolutionSigmaSPH->Draw("PE1");
    tAcrossPositionResidualEvolutionSigmaSPH->GetYaxis()->SetRangeUser(0., 0.5);

    tCanvas7->cd(9);
    gPad->SetLogz();
    tAlongPositionResidualEvolutionSPH->Draw("colz");

    tCanvas7->cd(10);
    tAlongPositionResidualEvolutionSigmaSPH->Draw("PE1");
    tAlongPositionResidualEvolutionSigmaSPH->GetYaxis()->SetRangeUser(0., 3.);

    tCanvas7->cd(11);
    gPad->SetLogz();
    tTimeResidualEvolutionSPH->Draw("colz");
//    tTimeResidualEvolutionSPH->ProfileX("_pfx", 1, -1, "s")->Draw("sameE1");

    tCanvas7->cd(12);
    tTimeResidualEvolutionSigmaSPH->Draw("PE1");
    tTimeResidualEvolutionSigmaSPH->GetYaxis()->SetRangeUser(0., 0.2);


    if(bBatchMode)
    {
      tCanvas7->SaveAs(Form("pl_global_differential_sph_response_%s.pdf", tCounterTag.Data()));

      delete tCanvas7;
    }


    TCanvas* tCanvas2 = new TCanvas(Form("tCanvas_local_%s", tCounterTag.Data()),Form("Local response of MRPC %s", tCounterTag.Data()),0,0,1750,1050);
    tCanvas2->Divide(5,3);

    tCanvas2->cd(1);
    gPad->SetLogx();
    gPad->SetLogy();
    tDistanceToLastHitUnbiased->Draw();

    tCanvas2->cd(2);
    gPad->SetLogx();
    tEfficiencyToLastHitUnbiased->Draw("AP");
    gPad->Update();
    tEfficiencyToLastHitUnbiased->GetPaintedGraph()->GetXaxis()->SetLimits(tEfficiencyToLastHitUnbiased->GetTotalHistogram()->GetXaxis()->GetXmin(),
                                                                           tEfficiencyToLastHitUnbiased->GetTotalHistogram()->GetXaxis()->GetXmax());
    gPad->Modified();

    tCanvas2->cd(4);
    gPad->SetLogx();
    gPad->SetLogy();
    tDistanceToLastHitBiased->Draw();

    tCanvas2->cd(5);
    gPad->SetLogx();
    tEfficiencyToLastHitBiased->Draw("AP");
    gPad->Update();
    tEfficiencyToLastHitBiased->GetPaintedGraph()->GetXaxis()->SetLimits(tEfficiencyToLastHitBiased->GetTotalHistogram()->GetXaxis()->GetXmin(),
                                                                         tEfficiencyToLastHitBiased->GetTotalHistogram()->GetXaxis()->GetXmax());
    gPad->Modified();

    tCanvas2->cd(6);
    gPad->SetLogx();
    gPad->SetLogz();
    tClusterSizeToLastHitUnbiased->Draw("colz");
    TProfile* tClusterSizeToLastHitUnbiasedProfile = tClusterSizeToLastHitUnbiased->ProfileX();
    tClusterSizeToLastHitUnbiasedProfile->Draw("sameE1");

    tCanvas2->cd(7);
    gPad->SetLogx();
    tClusterSizeToLastHitUnbiasedProfile->Draw("E1");

    tCanvas2->cd(8);
    gPad->SetLogx();
    gPad->SetLogy();
    tDistanceToLastPoint->Draw();

    tCanvas2->cd(9);
    gPad->SetLogx();
    gPad->SetLogz();
    tClusterSizeToLastHitBiased->Draw("colz");
    TProfile* tClusterSizeToLastHitBiasedProfile = tClusterSizeToLastHitBiased->ProfileX();
    tClusterSizeToLastHitBiasedProfile->Draw("sameE1");

    tCanvas2->cd(10);
    gPad->SetLogx();
    tClusterSizeToLastHitBiasedProfile->Draw("E1");

    tCanvas2->cd(11);
    gPad->SetLogx();
    gPad->SetLogz();
    tTimeResidualToLastHitUnbiased->Draw("colz");
    tTimeResidualToLastHitUnbiased->ProfileX("_pfx", 1, -1, "s")->Draw("sameE1");

    tCanvas2->cd(12);
    gPad->SetLogx();
    tTimeResidualToLastHitUnbiasedSigma->Draw("PE1");

    tCanvas2->cd(13);
    gPad->SetLogx();
    gPad->SetLogy();
    tDistanceToLocalLastPoint->Draw();

    tCanvas2->cd(14);
    gPad->SetLogx();
    gPad->SetLogz();
    tTimeResidualToLastHitBiased->Draw("colz");
    tTimeResidualToLastHitBiased->ProfileX("_pfx", 1, -1, "s")->Draw("sameE1");

    tCanvas2->cd(15);
    gPad->SetLogx();
    tTimeResidualToLastHitBiasedSigma->Draw("PE1");

    if(bBatchMode)
    {
      tCanvas2->SaveAs(Form("pl_local_response_tb_%s.pdf", tCounterTag.Data()));

      delete tCanvas2;
    }


    TCanvas* tCanvas3 = new TCanvas(Form("tCanvas_mc_flux_%s", tCounterTag.Data()),Form("MC point flux on MRPC %s", tCounterTag.Data()),0,0,1750,1050);
    tCanvas3->Divide(5,3);
    Int_t iNMaxTimeBins3 = countpads(tCanvas3);

//    TH3D* tLocalSurfaceFluxTotal = dynamic_cast<TH3D*>(tLocalSurfaceFlux->GetCopyTotalHisto());
    // histogram name: "tLocalSurfaceFlux_%d_%d_%d_total"
    tLocalSurfaceFluxTotal->SetName(Form("%s_copy", tLocalSurfaceFluxTotal->GetName()));
    // histogram name: "tLocalSurfaceFlux_%d_%d_%d_total_copy"

    for(Int_t iTimeBin = 1; iTimeBin <= iNMaxTimeBins3; iTimeBin++)
    {
      if(iTimeBin > tLocalSurfaceFluxTotal->GetNbinsZ())
      {
        break;
      }

      tLocalSurfaceFluxTotal->GetZaxis()->SetRange(iTimeBin, iTimeBin);
      TH2D* tLocalSurfaceFluxTotalProjection = dynamic_cast<TH2D*>(tLocalSurfaceFluxTotal->Project3D(Form("%d_yx", iTimeBin)));
      // histogram name: "tLocalSurfaceFlux_%d_%d_%d_total_copy_%d_yx"
      tLocalSurfaceFluxTotalProjection->SetTitle(Form("time window: [%2d,%2d] s", iTimeBin - 1, iTimeBin));
      tLocalSurfaceFluxTotalProjection->SetZTitle("MC point flux [Hz/cm^2]");

      tLocalSurfaceFluxTotalProjection->Clone(Form("tLocalSurfaceEfficiency%d_%s", iTimeBin, tCounterTag.Data()));

      tCanvas3->cd(iTimeBin);
      tLocalSurfaceFluxTotalProjection->SetStats(0);
      tLocalSurfaceFluxTotalProjection->Draw("colz");
    }

    if(bBatchMode)
    {
      tCanvas3->SaveAs(Form("pl_point_flux_tb_%s.pdf", tCounterTag.Data()));

      delete tCanvas3;
    }


    TCanvas* tCanvas4 = new TCanvas(Form("tCanvas_hit_flux_%s", tCounterTag.Data()),Form("Hit flux on MRPC %s", tCounterTag.Data()),0,0,1750,1050);
    tCanvas4->Divide(5,3);
    Int_t iNMaxTimeBins4 = countpads(tCanvas4);

//    TH3D* tLocalSurfaceFluxPassed = dynamic_cast<TH3D*>(tLocalSurfaceFlux->GetCopyPassedHisto());
    // histogram name: "tLocalSurfaceFlux_%d_%d_%d_passed"
    tLocalSurfaceFluxPassed->SetName(Form("%s_copy", tLocalSurfaceFluxPassed->GetName()));
    // histogram name: "tLocalSurfaceFlux_%d_%d_%d_passed_copy"

    for(Int_t iTimeBin = 1; iTimeBin <= iNMaxTimeBins4; iTimeBin++)
    {
      if(iTimeBin > tLocalSurfaceFluxPassed->GetNbinsZ())
      {
        break;
      }

      tLocalSurfaceFluxPassed->GetZaxis()->SetRange(iTimeBin, iTimeBin);
      TH2D* tLocalSurfaceFluxPassedProjection = dynamic_cast<TH2D*>(tLocalSurfaceFluxPassed->Project3D(Form("%d_yx", iTimeBin)));
      // histogram name: "tLocalSurfaceFlux_%d_%d_%d_passed_copy_%d_yx"
      tLocalSurfaceFluxPassedProjection->SetTitle(Form("time window: [%2d,%2d] s", iTimeBin - 1, iTimeBin));
      tLocalSurfaceFluxPassedProjection->SetZTitle("hit flux [Hz/cm^2]");

      tCanvas4->cd(iTimeBin);
      tLocalSurfaceFluxPassedProjection->SetStats(0);
      tLocalSurfaceFluxPassedProjection->Draw("colz");
    }

    if(bBatchMode)
    {
      tCanvas4->SaveAs(Form("pl_hit_flux_tb_%s.pdf", tCounterTag.Data()));

      delete tCanvas4;
    }


    TCanvas* tCanvas5 = new TCanvas(Form("tCanvas_efficiency_%s", tCounterTag.Data()),Form("Efficiency on MRPC %s", tCounterTag.Data()),0,0,1750,1050);
    tCanvas5->Divide(5,3);
    Int_t iNMaxTimeBins5 = countpads(tCanvas5);

    for(Int_t iTimeBin = 1; iTimeBin <= iNMaxTimeBins5; iTimeBin++)
    {
      if(iTimeBin > tLocalSurfaceFluxPassed->GetNbinsZ())
      {
        break;
      }

      TH2D* tLocalSurfaceEfficiency = dynamic_cast<TH2D*>(gFile->FindObjectAny(Form("tLocalSurfaceEfficiency%d_%s", iTimeBin, tCounterTag.Data())));
      TH2D* tLocalSurfacePassed = dynamic_cast<TH2D*>(gFile->FindObjectAny(Form("tLocalSurfaceFlux_%s_passed_copy_%d_yx", tCounterTag.Data(), iTimeBin)));
      TH2D* tLocalSurfaceTotal = dynamic_cast<TH2D*>(gFile->FindObjectAny(Form("tLocalSurfaceFlux_%s_total_copy_%d_yx", tCounterTag.Data(), iTimeBin)));

      tLocalSurfaceEfficiency->Divide(tLocalSurfacePassed, tLocalSurfaceTotal);

      tLocalSurfaceEfficiency->SetTitle(Form("time window: [%2d,%2d] s", iTimeBin - 1, iTimeBin));
      tLocalSurfaceEfficiency->SetZTitle("efficiency []");

      tCanvas5->cd(iTimeBin);
      tLocalSurfaceEfficiency->SetStats(0);
      tLocalSurfaceEfficiency->Draw("colz");
    }

    if(bBatchMode)
    {
      tCanvas5->SaveAs(Form("pl_efficiency_tb_%s.pdf", tCounterTag.Data()));

      delete tCanvas5;
    }

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

}


Int_t countpads(TVirtualPad *pad)
{
   //count the number of pads in pad
   if (!pad) return 0;
   Int_t npads = 0;
   TObject *obj;
   TIter next(pad->GetListOfPrimitives());
   while ((obj = next())) {
      if (obj->InheritsFrom(TVirtualPad::Class())) npads++;
   }
   return npads;
}

Double_t gpb(Double_t* x, Double_t* p)
{
  // p[0] : location
  // p[1] : scale
  // p[2] : shape
  // p[3] : normalization

  Double_t dReturnValue = p[3]*0.5/p[2]*( TMath::Erf((x[0] - p[0] + p[2]/2.)/p[1]/TMath::Sqrt(2.))
                                         -TMath::Erf((x[0] - p[0] - p[2]/2.)/p[1]/TMath::Sqrt(2.)));

  return dReturnValue;
}

Double_t chisqpdf(Double_t* x, Double_t* p)
{
  // p[0] : ndof
  // p[1] : normalization

  return p[1]*ROOT::Math::chisquared_pdf(x[0], p[0]);
}

Double_t chisqcdf(Double_t* x, Double_t* p)
{
  // p[0] : ndof

  return ROOT::Math::chisquared_cdf(x[0], p[0]);
}

Bool_t FindTofNode(const char* cGeoTag)
{
  // Locate the alphabetically first geometry file in 'geometry/tof' that contains
  // the given cGeoTag argument
  TString tDirectoryName = (TString)(gSystem->Getenv("VMCWORKDIR")) + "/geometry/tof/";
  FileStat_t tFileStat;
  if(1 == gSystem->GetPathInfo(tDirectoryName.Data(), tFileStat))
  {
    cout<<Form("Geometry directory %s does not exist.", tDirectoryName.Data())<<endl;
  }

  TRegexp* tRegexp = new TRegexp(((TString)("*")+(TString)(cGeoTag)+"."+"*root").Data(), kTRUE);

  TSystemDirectory* tSystemDirectory = new TSystemDirectory("dir", tDirectoryName.Data());

  TString tTargetDirectoryName(tDirectoryName);

  if(!tTargetDirectoryName.EndsWith("/"))
  {
      tTargetDirectoryName += "/";
  }

  TList* tList = tSystemDirectory->GetListOfFiles();
  tList->Sort();

  TIterator* tIter = tList->MakeIterator();
  TSystemFile* tSystemFile;

  TString tGeoFileName;
  Bool_t bFoundGeoFile(kFALSE);

  while(NULL != (tSystemFile = (TSystemFile*)tIter->Next()))
  {
    tGeoFileName = tSystemFile->GetName();

    if(tGeoFileName.Contains(*tRegexp))
    {
      tGeoFileName = tTargetDirectoryName + tGeoFileName;
      bFoundGeoFile = kTRUE;
      break;
    }
  }

  tList->Delete();

  if(!bFoundGeoFile)
  {
    cout<<"Could not find geometry file."<<endl;
    return kFALSE;
  }

  TFile::Open(tGeoFileName.Data(),"read");
  if(!gFile)
  {
    cout<<"Could not open geometry file."<<endl;
    return kFALSE;
  }


  gFile->FindObjectAny("FAIRGeom");
  if(!gGeoManager)
  {
    cout<<"Could not retrieve TGeoManager 'FAIRGeom' from geometry file."<<endl;
    cout<<"Create a TGeoManager and register the top volume in the file with it."<<endl;
    cout<<endl;

    // Use the FairRoot geometry interface to load the media
    FairGeoLoader* tGeoLoad = new FairGeoLoader("TGeo", "FairGeoLoader");
    FairGeoInterface* tGeoFace = tGeoLoad->getGeoInterface();
    FairGeoBuilder* tGeoBuild = tGeoLoad->getGeoBuilder();

    tGeoFace->setMediaFile((TString)(gSystem->Getenv("VMCWORKDIR"))+"/geometry/media.geo");
    tGeoFace->readMedia();

    // Create the media in the TGeoManager
    FairGeoMedia* tGeoMedia = tGeoFace->getMedia();

    TIter next(tGeoMedia->getListOfMedia());
    while(FairGeoMedium* tGeoMedium = dynamic_cast<FairGeoMedium*>(next()))
    {
      tGeoBuild->createMedium(tGeoMedium);
//      cout<<Form("Created medium %s.",tGeoMedium->GetName())<<endl;
    }

    gFile->Close();
    delete gFile;

    CbmGeoCave* tCave = new CbmGeoCave();
    tCave->setGeomFile((TString)(gSystem->Getenv("VMCWORKDIR"))+"/geometry/cave.geo");
    tGeoFace->addGeoModule(tCave);
    tGeoFace->readSet(tCave);
    tCave->create(tGeoBuild);

    FairModule* tModule = new CbmTof("TOF", kTRUE);
    tModule->SetGeometryFileName(tGeoFileName.Data());
    tModule->ConstructRootGeometry();
  }
  else
  {
    gFile->Close();
    delete gFile;
  }

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

  // Just in case a geometry file contains a TGeoManager with a 'tof' top node
  if(!TString(gGeoManager->GetCurrentNode()->GetName()).Contains("tof",TString::kIgnoreCase))
  {
    // Normally, the top node should be a 'cave' by which the 'tof' node is contained
    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 geometry file."<<endl;
      return kFALSE;
    }
  }

  gtTofNodePath = gGeoManager->GetPath();

  return kTRUE;
}

Bool_t GetCounterProperties()
{
  gGeoManager->cd(gtTofNodePath.Data());

  ExpandNode(gGeoManager->GetCurrentNode());

  gGeoManager->CdTop();

  if(gtCentralNodePath.IsNull() || !gGeoManager->CheckPath(gtCentralNodePath.Data()))
  {
    cout<<Form("Counter %d-%d-%d not found in geometry file.", giModuleType, giModuleIndex, giCounterIndex)<<endl;
    return kFALSE;
  }

  gtCentralPN = new TGeoPhysicalNode(gtCentralNodePath.Data());

  giNumberGlassPlates = 0;
  giNumberGasGaps = 0;

  Double_t dLocalCoordinates[3] = {0.,0.,0.};
  Double_t dGlobalCoordinates[3] = {0.,0.,0.};

  gtCentralPN->GetMatrix()->LocalToMaster(dLocalCoordinates,dGlobalCoordinates);

  Double_t dAxisLow(0.);
  Double_t dAxisHigh(0.);

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

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

        gdGlassThickness = tPlateShape->GetAxisRange(3,dAxisLow,dAxisHigh);
      }

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

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

        gdGapXWidth = tGapShape->GetAxisRange(1,dAxisLow,dAxisHigh);
        gdGapYWidth = tGapShape->GetAxisRange(2,dAxisLow,dAxisHigh);
        gdGapSize = tGapShape->GetAxisRange(3,dAxisLow,dAxisHigh);

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

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

        gdCellXWidth = tFirstCellShape->GetAxisRange(1,dAxisLow,dAxisHigh);
        gdCellYWidth = tFirstCellShape->GetAxisRange(2,dAxisLow,dAxisHigh);


        if(gdCellXWidth == gdGapXWidth)
        {
          gbCellsAlongX = kTRUE;

          gdCellAcrossWidth = gdCellYWidth;
          gdCellAlongWidth = gdCellXWidth;

          gdGapAcrossWidth = gdGapYWidth;
          gdGapAlongWidth = gdGapXWidth;
        }
        else
        {
          gdCellAcrossWidth = gdCellXWidth;
          gdCellAlongWidth = gdCellYWidth;

          gdGapAcrossWidth = gdGapXWidth;
          gdGapAlongWidth = gdGapYWidth;
        }

        if(gdCellXWidth <= gdCellYWidth)
        {
          if(gdCellXWidth/gdCellYWidth >= 0.5)
          {
            gbPadReadout = kTRUE;
          }
        }
        else
        {
          if(gdCellXWidth/gdCellYWidth < 2.)
          {
            gbPadReadout = kTRUE;
          }
        }

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

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

        if(gbCellsAlongX)
        {
          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.)
          {
            gbCellIndexingAlongAxis = kTRUE;
          }
          else if(dLastCellMasterYLow < 0. && dFirstCellMasterYHigh > 0.)
          {
            gbCellIndexingAlongAxis = kFALSE;
          }
          else
          {
            cout<<"Incomplete understanding of geometry. Cannot determine direction of cell indexing."<<endl;
            return kFALSE;
          }
        }
        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.)
          {
            gbCellIndexingAlongAxis = kTRUE;
          }
          else if(dLastCellMasterXLow < 0. && dFirstCellMasterXHigh > 0.)
          {
            gbCellIndexingAlongAxis = kFALSE;
          }
          else
          {
            cout<<"Incomplete understanding of geometry. Cannot determine direction of cell indexing."<<endl;
            return kFALSE;
          }
        }
      }

      giNumberGasGaps++;
    }
  }

  delete gtCentralPN;

  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;

      gtCurrentNodePath = gtTofNodePath + "/" + (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();
    }
    // 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 during (optional) beam point creation
    else if(TString(tDaughterNode->GetName()).Contains("counter"))
    {
      gtCurrentNodePath += "/" + (TString)tDaughterNode->GetName();

      giCurrentCounterIndex = tDaughterNode->GetNumber();

      gtCurrentCentralNodePath = gtCurrentNodePath;

      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)
      {
        gtCurrentCentralNodePath += TString::Format("/tof_glass_%d",iNGaps/2);
      }
      else
      {
        gtCurrentCentralNodePath += TString::Format("/Gap_%d",(iNGaps-1)/2);
      }

      if(giModuleType == giCurrentModuleType && giModuleIndex == giCurrentModuleIndex && giCounterIndex == giCurrentCounterIndex)
      {
        gtCentralNodePath = gtCurrentCentralNodePath;
        gtCounterNode = tDaughterNode;
      }
    }
    else
    {
      gtCurrentNodePath += "/" + (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
    gtCurrentNodePath.ReplaceAll("/"+(TString)tDaughterNode->GetName(), "");
  }
}