// ------------------------------------------------------------------
// -----                     CbmTofAnaTestbeam                        -----
// -----              Created 12/04/2014 by nh                 -----
// ------------------------------------------------------------------

#include "CbmTofAnaTestbeam.h"

// TOF Classes and includes
#include "CbmTofDigi.h"       // in cbmdata/tof
#include "CbmTofDigiExp.h"    // in cbmdata/tof
#include "CbmTofHit.h"        // in cbmdata/tof
#include "CbmTofGeoHandler.h" // in tof/TofTools
#include "CbmTofDetectorId_v12b.h" // in cbmdata/tof
#include "CbmTofDetectorId_v14a.h" // in cbmdata/tof
#include "CbmTofCell.h"       // in tof/TofData
#include "CbmTofDigiPar.h"    // in tof/TofParam
#include "CbmTofDigiBdfPar.h" // in tof/TofParam
#include "CbmTofAddress.h"    // in cbmdata/tof
#include "CbmMatch.h"

#include "CbmTofTestBeamClusterizer.h"

#include "TMbsMappingTofPar.h"

// CBMroot classes and includes

// FAIR classes and includes
#include "FairRunAna.h"
#include "FairRuntimeDb.h"
#include "FairRootManager.h"
#include "FairLogger.h"

// ROOT Classes and includes
#include "Riostream.h"
#include "TClonesArray.h"
#include "TH1.h"
#include "TH2.h"
#include "TProfile.h"
#include "TString.h"
#include "TFile.h"
#include "TMath.h"
#include "TRandom.h"
#include "TROOT.h"
#include "TDirectory.h"
#include "TGeoManager.h"

const Int_t DetMask = 4194303;
const Double_t DTDMAX=6000.;  // diamond inspection range in ps  

Double_t dTDia;
Double_t dDTD4Min=1.E8;

//___________________________________________________________________
//
// CbmTofAnaTestbeam
//
// Task for analysis of Testbeam data
//
// ------------------------------------------------------------------
CbmTofAnaTestbeam::CbmTofAnaTestbeam()
  : FairTask("HadronAnalysis"),
    fEvents(0),
    fGeoHandler(new CbmTofGeoHandler()),
    fTofId(NULL),
    fChannelInfo(NULL),
    iNbSmTot(0),
    fvTypeSmOffs(),
    iNbRpcTot(0),
    fvSmRpcOffs(),
    iNbChTot(0),
    fvRpcChOffs(),
    fDigiPar(NULL),
    fDigiBdfPar(NULL),
    fTofDigisColl(NULL),
    fTofDigiMatchColl(NULL),
    fTofHitsColl(NULL),
    fStart(),
    fStop()
{
  cout << "CbmTofTests: Task started " << endl;
}
// ------------------------------------------------------------------

// ------------------------------------------------------------------
CbmTofAnaTestbeam::CbmTofAnaTestbeam(const char* name, Int_t verbose)
  : FairTask(name, verbose),
    fEvents(0),
    fGeoHandler(new CbmTofGeoHandler()),
    fTofId(NULL),
    fChannelInfo(NULL),
    iNbSmTot(0),
    fvTypeSmOffs(),
    iNbRpcTot(0),
    fvSmRpcOffs(),
    iNbChTot(0),
    fvRpcChOffs(),
    fDigiPar(NULL),
    fDigiBdfPar(NULL),
    fTofDigisColl(NULL),
    fTofDigiMatchColl(NULL),
    fTofHitsColl(NULL),
    fStart(),
    fStop()
{
}
// ------------------------------------------------------------------

// ------------------------------------------------------------------
CbmTofAnaTestbeam::~CbmTofAnaTestbeam()
{
    // Destructor
}
// ------------------------------------------------------------------
/************************************************************************************/
// FairTasks inherited functions
InitStatus CbmTofAnaTestbeam::Init()
{
   if( kFALSE == RegisterInputs() )
      return kFATAL;

   //   fTofId = new ( CbmTofDetectorId )CbmTofDetectorId_v14a();
   if( kFALSE == InitParameters() )
      return kFATAL;

   if( kFALSE == LoadCalParameter() )
      return kFATAL;

   if( kFALSE == CreateHistos() )
      return kFATAL;

   return kSUCCESS;
}

void CbmTofAnaTestbeam::SetParContainers()
{
   LOG(INFO)<<" CbmTofAnaTestbeam => Get the digi parameters for tof"<<FairLogger::endl;
   return;
   // Get Base Container
   FairRunAna* ana = FairRunAna::Instance();
   FairRuntimeDb* rtdb=ana->GetRuntimeDb();

   fDigiPar = (CbmTofDigiPar*) (rtdb->getContainer("CbmTofDigiPar"));

   fDigiBdfPar = (CbmTofDigiBdfPar*) (rtdb->getContainer("CbmTofDigiBdfPar"));
}

void CbmTofAnaTestbeam::Exec(Option_t * option)
{
   // Task execution

   LOG(DEBUG)<<" CbmTofAnaTestbeam => New event"<<FairLogger::endl;

   fStart.Set();
   FillHistos();
   fStop.Set();


   if( 0 == ( fEvents%100000 ) && 0 < fEvents )
   {
      cout << "-I- CbmTofAnaTestbeam::Exec : "
           << "event " << fEvents << " processed." << endl;
   }
   fEvents += 1;
}

void CbmTofAnaTestbeam::Finish()
{
   // Normalisations
   cout << "CbmTofAnaTestbeam::Finish up with " << fEvents << " analyzed events " << endl;

   WriteHistos();
   // Prevent them from being sucked in by the CbmHadronAnalysis WriteHistograms method
   DeleteHistos();
}

/************************************************************************************/
// Functions common for all clusters approximations
Bool_t   CbmTofAnaTestbeam::RegisterInputs()
{
   FairRootManager *fManager = FairRootManager::Instance();

   fTofDigisColl   = (TClonesArray *) fManager->GetObject("CbmTofDigiExp");

   if( NULL == fTofDigisColl)
      fTofDigisColl = (TClonesArray *) fManager->GetObject("CbmTofDigi");

   if( NULL == fTofDigisColl)
   {
      LOG(ERROR)<<"CbmTofAnaTestbeam::RegisterInputs => Could not get the TofDigi TClonesArray!!!"<<FairLogger::endl;
      return kFALSE;
   } // if( NULL == fTofDigisColl)

   fTofHitsColl   = (TClonesArray *) fManager->GetObject("TofHit");
   if( NULL == fTofHitsColl)
   {
      LOG(ERROR)<<"CbmTofAnaTestbeam::RegisterInputs => Could not get the TofHit TClonesArray!!!"<<FairLogger::endl;
      return kFALSE;
   } // if( NULL == fTofHitsColl)

   fTofDigiMatchColl= (TClonesArray *) fManager->GetObject("TofDigiMatch");
   if( NULL == fTofDigiMatchColl)
   {
      LOG(ERROR)<<"CbmTofAnaTestbeam::RegisterInputs => Could not get the Match TClonesArray!!!"<<FairLogger::endl;
      return kFALSE;
   } // if( NULL == fTofDigiMatchColl)

   return kTRUE;
}
/************************************************************************************/
Bool_t   CbmTofAnaTestbeam::InitParameters()
{

   // Initialize the TOF GeoHandler
   Bool_t isSimulation=kFALSE;
   LOG(INFO)<<"CbmTofAnaTestbeam::InitParameters - Geometry, Mapping, ... "
             <<FairLogger::endl;
   // Get Base Container
   FairRun* ana = FairRun::Instance();
   FairRuntimeDb* rtdb=ana->GetRuntimeDb();   

   Int_t iGeoVersion = fGeoHandler->Init(isSimulation);
   if( k14a > iGeoVersion )
   {

   }

   fDigiPar = (CbmTofDigiPar*) (rtdb->getContainer("CbmTofDigiPar"));
   if( 0 == fDigiPar )
   {
      LOG(ERROR)<<"CbmTofAnaTestbeam::InitParameters => Could not obtain the CbmTofDigiPar "<<FairLogger::endl;
      return kFALSE; 
   }

   fDigiBdfPar = (CbmTofDigiBdfPar*) (rtdb->getContainer("CbmTofDigiBdfPar"));
   if( 0 == fDigiBdfPar )
   {
      LOG(ERROR)<<"CbmTofAnaTestbeam::InitParameters => Could not obtain the CbmTofDigiBdfPar "<<FairLogger::endl;
      return kFALSE; 
   }

   // Mapping parameter
   fMbsMappingPar = (TMbsMappingTofPar*) (rtdb->getContainer("TMbsMappingTofPar"));
   if( 0 == fMbsMappingPar )
   {
      LOG(ERROR)<<"CbmTofAnaTestBeam::InitParameters => Could not obtain the TMbsMappingTofPar "<<FairLogger::endl;
      return kFALSE; 
   }

   rtdb->initContainers(  ana->GetRunId() );

   LOG(INFO)<<"CbmTofAnaTestbeam::InitParameter: currently " 
            << fDigiPar->GetNrOfModules() << " digi cells " <<FairLogger::endl;

   return kTRUE;
}
/************************************************************************************/
Bool_t   CbmTofAnaTestbeam::LoadCalParameter()
{
    if(fCalParFileName.IsNull()) return kTRUE;

    fCalParFile = new TFile(fCalParFileName,"");
    if(NULL == fCalParFile) {
      LOG(ERROR) << "CbmTofAnaTestBeam::LoadCalParameter: "
		 << "file " << fCalParFileName << " does not exist!" << FairLogger::endl;
      return kTRUE;
    }

    TProfile *fhtmp=(TProfile *) gDirectory->FindObjectAny( Form("hDTD4DT04D4best_pfx_px"));
    if (NULL == fhtmp) {
      LOG(INFO)<<" Histo " << Form("hDTD4DT04D4best_pfx_px") << " not found. "
             <<FairLogger::endl;
    }

    TProfile *fhtmpx=(TProfile *) gDirectory->FindObjectAny( Form("hDTX4D4best_pfx_px"));
    if (NULL == fhtmpx) {
      LOG(INFO)<<" Histo " << Form("hDTX4D4best_pfx_px") << " not found. "
             <<FairLogger::endl;
    }

    TProfile *fhtmpy=(TProfile *) gDirectory->FindObjectAny( Form("hDTY4D4best_pfx_px"));
    if (NULL == fhtmpy) {
      LOG(INFO)<<" Histo " << Form("hDTY4D4best_pfx_px") << " not found. "
             <<FairLogger::endl;
    }

    TProfile *fhtmpt=(TProfile *) gDirectory->FindObjectAny( Form("hTexpDT04D4best_pfx_px"));
    if (NULL == fhtmpt) {
      LOG(INFO)<<" Histo " << Form("hTexpDT04D4best_pfx_px") << " not found. "
             <<FairLogger::endl;
    }

    gROOT->cd();
    if(NULL != fhtmp)  fhDTD4DT04D4Off=(TH1D *)fhtmp->Clone();
    if(NULL != fhtmpx) fhDTX4D4Off=(TH1D *)fhtmpx->Clone();
    if(NULL != fhtmpy) fhDTY4D4Off=(TH1D *)fhtmpy->Clone();
    if(NULL != fhtmpt) fhDTTexpD4Off=(TH1D *)fhtmpt->Clone();

    fCalParFile->Close();
    //    fhDTD4DT04D4Off->Draw();
    
    if (fdDCh4Sel==0.)  fdDCh4Sel=1000.;  // open default window
    if (fdPosY4Sel==0.) fdPosY4Sel=10.;   // open default window

    return kTRUE;
}
/************************************************************************************/
// ------------------------------------------------------------------
Bool_t CbmTofAnaTestbeam::CreateHistos()
{
  // Create histogramms
   Double_t XDMAX=100.;
   Double_t YDMAX=100.;
   Double_t TDMAX=100000.;
   Double_t DTMAX=1000.;
   Double_t DXMAX=10.;
   Double_t DYMAX=10.;
   Double_t XMAX=100.;
   Double_t YMAX=100.;

   TDirectory * oldir = gDirectory; // <= To prevent histos from being sucked in by the param file of the TRootManager!
   gROOT->cd(); // <= To prevent histos from being sucked in by the param file of the TRootManager !

   // define histos here
   Int_t iNbDet=fMbsMappingPar->GetNbMappedDet();
   fhXYPos.resize( iNbDet  );
   for(Int_t iDet=0; iDet<iNbDet; iDet++){
     fhXYPos[iDet] =  new TH2F( Form("hXY_SmT%d",iDet),
				Form("XY Position correlation of SmT %d; X[cm]; Y [cm]",iDet),
			       100, -YDMAX, YDMAX, 100, -YDMAX, YDMAX);       
   }

    fhXX2 =  new TH2F( 
          Form("hXX2"),
          Form("Plastic position correlation; XPla0; XPla1"),
          100, -XDMAX, XDMAX, 100, -XDMAX, XDMAX);
    fhYY2 =  new TH2F( 
          Form("hYY2"),
          Form("Plastic position correlation; YPla0; YPla1"),
          100, -YDMAX, YDMAX, 100, -YDMAX, YDMAX);
    fhDT2 =  new TH1F( 
          Form("hDT2"),
          Form("Plastic time difference; TPla0 - TPla1 (ps)"),
          100, -1000., 1000.);
 
   for (Int_t iDet=0; iDet<2; iDet++)
   {
     fhXX02[iDet] =  new TH2F( Form("hXX02_%d",iDet),Form("X Position correlation; X0; X2"),
			       100, -YDMAX, YDMAX, 100, -YDMAX, YDMAX);      
     fhYY02[iDet] =  new TH2F( Form("hYY02_%d",iDet),Form("Y Position correlation; Y0; Y2"),
			       100, -YDMAX, YDMAX, 100, -YDMAX, YDMAX); 
   }

     fhXX04 =  new TH2F( Form("hXX04"),Form("X Position correlation; X0; X4"),
			       200, -YDMAX, YDMAX, 200, -YDMAX, YDMAX);      
     fhYY04 =  new TH2F( Form("hYY04"),Form("Y Position correlation; Y0; Y4"),
			       200, -YDMAX, YDMAX, 200, -YDMAX, YDMAX); 
     fhTT04 =  new TH2F( Form("hTT04"),Form("Time  correlation; T0; T4"),
			       100, -TDMAX, TDMAX, 100, -TDMAX, TDMAX); 
     fhDXDY04 = new TH2F( Form("hDXDY04"),Form("position correlation; #Delta x [cm]; #DeltaY [cm]"),
			       100, -DXMAX, DXMAX, 100, -DYMAX, DYMAX); 
     fhDXDT04 = new TH2F( Form("hDXDT04"),Form("Time - position correlation; #Delta x [cm]; #DeltaT [ps]"),
			       100, -DXMAX, DXMAX, 100, -DTMAX, DTMAX); 
     fhDYDT04 = new TH2F( Form("hDYDT04"),Form("Time - position correlation; #Delta y [cm]; #DeltaT [ps]"),
			       100, -DYMAX, DYMAX, 100, -DTMAX, DTMAX); 
     fhChi04 =  new TH1F( Form("hChi04"),Form("Matching Chi2; #chi; Nhits"),
			  100, 0., 10.); 

     fhXY24   = new TH2F( Form("hXY24"),Form("position of matches; x [cm]; Y [cm]"),
			       100, -XMAX, XMAX, 100, -YMAX, YMAX); 
     fhDXDY24 = new TH2F( Form("hDXDY24"),Form("position correlation; #Delta x [cm]; #DeltaY [cm]"),
			       100, -DXMAX, DXMAX, 100, -DYMAX, DYMAX); 
     fhDXDT24 = new TH2F( Form("hDXDT24"),Form("Time - position correlation; #Delta x [cm]; #DeltaT [ps]"),
			       100, -DXMAX, DXMAX, 100, -DTMAX*5., DTMAX*5.); 
     fhDYDT24 = new TH2F( Form("hDYDT24"),Form("Time - position correlation; #Delta y [cm]; #DeltaT [ps]"),
			       100, -DYMAX, DYMAX, 100, -DTMAX*5., DTMAX*5.); 
     fhChi24 =  new TH1F( Form("hChi24"),Form("Matching Chi2; #chi; Nhits"),
			  100, 0., 10.); 

     fhDXDY04best = new TH2F( Form("hDXDY04best"),Form("position correlation; #Delta x [cm]; #DeltaY [cm]"),
			       100, -DXMAX, DXMAX, 100, -DYMAX, DYMAX); 
     fhDXDT04best = new TH2F( Form("hDXDT04best"),Form("time - position correlation; #Delta x [cm]; #DeltaT [ps]"),
			       100, -DXMAX, DXMAX, 100, -DTMAX, DTMAX); 
     fhDYDT04best = new TH2F( Form("hDYDT04best"),Form("time - position correlation; #Delta y [cm]; #DeltaT [ps]"),
			       100, -DYMAX, DYMAX, 100, -DTMAX, DTMAX); 

     fhX0DT04best = new TH2F( Form("hX0DT04best"),Form("time - position correlation; #Delta x [cm]; #DeltaT [ps]"),
			       100, -50., 50., 100, -DTMAX, DTMAX); 
     fhY0DT04best = new TH2F( Form("hY0DT04best"),Form("time - position correlation; #Delta y [cm]; #DeltaT [ps]"),
			       100, -50., 50., 100, -DTMAX, DTMAX); 

     fhChi04best  =  new TH1F( Form("hChi04best"),Form("matching chi2; #chi; Nhits"),
			  100, 0., 100.);      
     fhDigiMul0best =  new TH1F( Form("hDigiMul0best"),Form("Number of digis in cluster; N_{digi}; "),
			  20, 0., 20.); 
     fhDigiMul4best =  new TH1F( Form("hDigiMul4best"),Form("Number of digis in cluster; N_{digi}; "),
			  20, 0., 20.); 

     fhChiDT04best = new TH2F( Form("hChiDT04best"),Form("Time - Chi correlation; #chi; #DeltaT [ps]"),
			       100, 0., 100., 100, -DTMAX, DTMAX);
     fhDT24DT04best = new TH2F( Form("hDT24DT04best"),
			    Form("Time - position correlation; #DeltaT24 [ps]; #DeltaT04 [ps]"),
			    100, -DTMAX*3., DTMAX*3., 100, -DTMAX, DTMAX); 
     fhDTD4DT04best = new TH2F( Form("hDTD4DT04best"),
			    Form("Time - velocity correlation; #DeltaTD4 [ps]; #DeltaT04 [ps]"),
			    100, -DTMAX*6., DTMAX*6., 100, -DTMAX, DTMAX); 

     fhDXDY04D4best = new TH2F( Form("hDXDY04D4best"),Form("position correlation; #Delta x [cm]; #DeltaY [cm]"),
			       100, -DXMAX, DXMAX, 100, -DYMAX, DYMAX); 
     fhDXDT04D4best = new TH2F( Form("hDXDT04D4best"),Form("time - position correlation; #Delta x [cm]; #DeltaT [ps]"),
			       100, -DXMAX, DXMAX, 100, -DTMAX, DTMAX); 
     fhDYDT04D4best = new TH2F( Form("hDYDT04D4best"),Form("time - position correlation; #Delta y [cm]; #DeltaT [ps]"),
			       100, -DYMAX, DYMAX, 100, -DTMAX, DTMAX); 
     fhDistDT04D4best = new TH2F( Form("hDistDT04D4best"),Form("time - distance correlation; Dist [cm]; #DeltaT [ps]"),
			       100, fdHitDistMin, fdHitDistMin + 5., 100, -DTMAX, DTMAX); 
     fhTexpDT04D4best = new TH2F( Form("hTexpDT04D4best"),Form("measured - expected time - correlation; Texp [ps]; #DeltaT [ps]"),
			       100, 0., 3000., 100, -DTMAX, DTMAX); 
     fhCluSize0DT04D4best = new TH2F( Form("hCluSize0DT04D4best"),
	 			     Form("time - CluSize correlation; N_{strips} ; #DeltaT [ps]"),
				     15, 1, 16, 100, -DTMAX, DTMAX); 
     fhCluSize4DT04D4best = new TH2F( Form("hCluSize4DT04D4best"),
	 			     Form("time - CluSize correlation; N_{strips} ; #DeltaT [ps]"),
				     15, 1, 16, 100, -DTMAX, DTMAX);
     fhTot0DT04D4best = new TH2F( Form("hTot0DT04D4best"),
	 			     Form("time - Tot correlation; TOT ; #DeltaT [ps]"),
				     100, 0., 10000., 100, -DTMAX, DTMAX);  
     fhTot4DT04D4best = new TH2F( Form("hTot4DT04D4best"),
	 			     Form("time - Tot correlation; TOT ; #DeltaT [ps]"),
				     100, 0., 10000., 100, -DTMAX, DTMAX);  

     fhX0DT04D4best = new TH2F( Form("hX0DT04D4best"),Form("time - position correlation; #Delta x [cm]; #DeltaT [ps]"),
			       100, -50., 50., 100, -DTMAX, DTMAX); 
     fhY0DT04D4best = new TH2F( Form("hY0DT04D4best"),Form("time - position correlation; #Delta y [cm]; #DeltaT [ps]"),
			       100, -50., 50., 100, -DTMAX, DTMAX); 

     fhChi04D4best  =  new TH1F( Form("hChi04D4best"),Form("matching chi2; #chi; Nhits"),
			  100, 0., 10.);
     fhTofD4best  =  new TH1F( Form("hTofD4best"),Form("tof D4; t [ps]; Counts"),
			  100, 0., 50000.); 
     fhVelD4best  =  new TH1F( Form("hVelD4best"),Form("vel D4; v [cm/ns]; Counts"),
			  100, 0., 50.); 

     fhDigiMul0D4best =  new TH1F( Form("hDigiMul0D4best"),Form("Number of digis in cluster; N_{digi}; "),
			  20, 0., 20.); 
     fhDigiMul4D4best =  new TH1F( Form("hDigiMul4D4best"),Form("Number of digis in cluster; N_{digi}; "),
			  20, 0., 20.); 
     fhCluMul04D4best =  new TH2F( Form("hCluMul04D4best"),Form("cluster multiplicity ; Mul0; Mul4"),
				   10, 0., 10., 10, 0., 10.); 


     fhChiDT04D4best = new TH2F( Form("hChiDT04D4best"),Form("Time - position correlation; #chi; #DeltaT [ps]"),
			       100, 0., 100., 100, -DTMAX, DTMAX);
     fhDT24DT04D4best = new TH2F( Form("hDT24DT04D4best"),
			    Form("Time - position correlation; #DeltaT24 [ps]; #DeltaT04 [ps]"),
			    100, -DTMAX*3., DTMAX*3., 100, -DTMAX, DTMAX); 
     fhDTD4DT04D4best = new TH2F( Form("hDTD4DT04D4best"),
			    Form("Time - velocity correlation; #DeltaTD4 [ps]; #DeltaT04 [ps]"),
			    100, -DTMAX*6., DTMAX*6., 100, -DTMAX, DTMAX); 

     Double_t dXMAX=30.;
     Double_t dYMAX=20.; 
     fhDTMul4D4best = new TH2F( Form("hDTMul4D4best"),
			    Form("Time - Multiplicity correlation; Mul4 ; #DeltaT04 [ps]"),
			    10, 0., 10., 100, -DTMAX, DTMAX);
     fhDTX4D4best = new TH2F( Form("hDTX4D4best"),
			    Form("Time - position correlation; X4 [cm]; #DeltaT04 [ps]"),
			    50, -15., 15., 100, -DTMAX, DTMAX); 
     fhDTY4D4best = new TH2F( Form("hDTY4D4best"),
			    Form("Time - position correlation; Y4 [cm]; #DeltaT04 [ps]"),
			    50, -5., 5., 100, -DTMAX, DTMAX); 
     fhDXX4D4best = new TH2F( Form("hDXX4D4best"),
			    Form("DX - position correlation; X4 [cm]; #DeltaX04 [cm]"),
			    50, -15., 15., 50, -10., 10.); 
     fhDXY4D4best = new TH2F( Form("hDXY4D4best"),
			    Form("DX - position correlation; Y4 [cm]; #DeltaX04 [cm]"),
			    50, -5., 5., 50, -10., 10.); 
     fhDYX4D4best = new TH2F( Form("hDYX4D4best"),
			    Form("DY - position correlation; X4 [cm]; #DeltaY04 [cm]"),
			    50, -15., 15., 50, -10., 10.); 
     fhDYY4D4best = new TH2F( Form("hDYY4D4best"),
			    Form("DY - position correlation; Y4 [cm]; #DeltaY04 [cm]"),
			    50, -5., 5., 50, -10., 10.); 

     fhDTMul0D4best = new TH2F( Form("hDTMul0D4best"),
			    Form("Time - Multiplicity correlation; Mul0 ; #DeltaT04 [ps]"),
			    10, 0., 10., 100, -DTMAX, DTMAX);
     fhDTX0D4best = new TH2F( Form("hDTX0D4best"),
			    Form("Time - position correlation; X0 [cm]; #DeltaT04 [ps]"),
			    50, -25., 25., 100, -DTMAX, DTMAX); 
     fhDTY0D4best = new TH2F( Form("hDTY0D4best"),
			    Form("Time - position correlation; Y0 [cm]; #DeltaT04 [ps]"),
			    50, -25., 25., 100, -DTMAX, DTMAX); 
     fhDXX0D4best = new TH2F( Form("hDXX0D4best"),
			    Form("DX - position correlation; X0 [cm]; #DeltaX04 [cm]"),
			    50, -25., 25., 50, -15., 15.); 
     fhDXY0D4best = new TH2F( Form("hDXY0D4best"),
			    Form("DX - position correlation; Y0 [cm]; #DeltaX04 [cm]"),
			    50, -25., 25., 50, -15., 15.); 
     fhDYX0D4best = new TH2F( Form("hDYX0D4best"),
			    Form("DY - position correlation; X0 [cm]; #DeltaY04 [cm]"),
			    50, -25., 25., 50, -15., 15.); 
     fhDYY0D4best = new TH2F( Form("hDYY0D4best"),
			    Form("DY - position correlation; Y0 [cm]; #DeltaY04 [cm]"),
			    50, -25., 25., 50, -15., 15.); 

     fhNMatch04 =  new TH1F( Form("hNMatch04"),Form("Number of Matched Hit pairs 0-4; NMatched"),
			  50, 0., 50.); 

     fhNMatch04sel =  new TH1F( Form("hNMatch04sel"),Form("Number of Matched Hit pairs 0-4; NMatched"),
			  50, 0., 50.); 

     fhNMatchD4sel =  new TH1F( Form("hNMatchD4sel"),Form("Number of Matched Hit pairs 0-4; NMatched"),
			  50, 0., 50.); 

     fhNMatch24 =  new TH1F( Form("hNMatch24"),Form("Number of Matched Hit pairs 2-4; NMatched"),
			  20, 0., 20.); 

     fhNMatch24sel =  new TH1F( Form("hNMatch24sel"),
                          Form("Number of Matched Hit pairs 0-4 uder 2-4 selection; NMatched"),
			  20, 0., 20.); 

     fhDT24sel =  new TH1F( Form("hDT24sel"),Form("Time difference Pla - HdRef; #DeltaT24 (ps)"),
			  100, -DTMAX*3., DTMAX*3.); 

     fhDTD4sel =  new TH1F( Form("hDTD4sel"),Form("Time difference Dia - HdRef; #DeltaTD4 (ps)"),
			  100, -DTDMAX, DTDMAX); 



   gDirectory->cd( oldir->GetPath() ); // <= To prevent histos from being sucked in by the param file of the TRootManager!

   return kTRUE;
}

// ------------------------------------------------------------------
Bool_t CbmTofAnaTestbeam::FillHistos()
{
   // Constants, TODO => put as parameter !!!

  Int_t kTOF=6;

   // Declare variables outside the loop
   CbmTofHit   *pHit;
   CbmTofHit   *pHit1;
   CbmTofHit   *pHit2;
   CbmTofHit   *pHitRef;
   CbmTofCell  *fChannelInfo1;
   CbmTofCell  *fChannelInfo2;

   Int_t iNbTofDigis, iNbTofHits;


   iNbTofDigis   = fTofDigisColl->GetEntriesFast();
   iNbTofHits    = fTofHitsColl->GetEntriesFast();
   /*
   LOG(INFO)<<Form("CbmTofAnaTestbeam::FillHistos for %d digis and %d Tof hits",iNbTofDigis,iNbTofHits)
	    <<FairLogger::endl;
   */
   // Digis info
   /*
   Double_t dTotalDataSize = 0;
   if( kTRUE == fDigiBdfPar->UseExpandedDigi() )
   {
      CbmTofDigiExp *pDigi;
      for( Int_t iDigInd = 0; iDigInd < iNbTofDigis; iDigInd++ )
      {
         pDigi = (CbmTofDigiExp*) fTofDigisColl->At( iDigInd );

         Int_t iSmType = pDigi->GetType();
         Int_t iSm     = pDigi->GetSm();
         Int_t iRpc    = pDigi->GetRpc();
         Int_t iCh     = pDigi->GetChannel();

      } // for( Int_t iDigInd = 0; iDigInd < iNbTofDigis; iDigInd++ )
   } // if( kTRUE == fDigiBdfPar->UseExpandedDigi() )
   */
   // Hits info
   Int_t iNbMatchedHits = 0;
   Int_t iNbMaxMatch=100;
   Double_t Zref=200.;
   Double_t Chi2MatchMin=1.E8;
   Double_t Chi2List[iNbMaxMatch];
   CbmTofHit *pChi2Hit1[iNbMaxMatch];
   CbmTofHit *pChi2Hit2[iNbMaxMatch];
   Chi2List[0]=1.E8;
   pChi2Hit1[0]=NULL;
   pChi2Hit2[0]=NULL;
   Int_t iNSel=2;
   Bool_t  BSel[iNSel];
   for(Int_t iSel=0;iSel<iNSel;iSel++){BSel[iSel]=kFALSE;}
   Double_t dDT24=0.;
   Double_t dDTD4=0.;
   Int_t iMat24=0;
   Double_t Chi2Mat24Min=1.E8;
   dTDia=1.E8;
   dDTD4Min=1.E8;
   Double_t dMul0=0.;
   Double_t dMul4=0.;
   Double_t dMulD=0.;
   Double_t dM4Max=1; // modify
   if (fdMul4Max>0) dM4Max=fdMul4Max;
   Double_t dM0Max=100; // modify
   if (fdMul0Max>0) dM0Max=fdMul0Max;
   Double_t dMDMax=1; 
   if (fdMulDMax>0) dMDMax=fdMulDMax;
   Double_t hitpos[3];
   Double_t hitpos_local[3];

   for( Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++)
   {
      pHit = (CbmTofHit*) fTofHitsColl->At( iHitInd );
      if(NULL == pHit) continue;
      Int_t iDetId = (pHit->GetAddress() & DetMask);
      Int_t iChId = pHit->GetAddress();
      fChannelInfo = fDigiPar->GetCell( iChId );
      Int_t iSmType=CbmTofAddress::GetSmType( iDetId );

      LOG(DEBUG)<<Form("CbmTofAnaTestbeam::FillHistos: process %d.(%d) Tof hit 0x%08x, x = %6.1f, y = %6.1f, z=%6.1f, t=%10.1f",
		      iHitInd, iNbTofHits, iChId, 
		      pHit->GetX(), pHit->GetY(), pHit->GetZ(), pHit->GetTime())
	    <<FairLogger::endl;
      
      if(NULL == fChannelInfo){
        LOG(FATAL) << "CbmTofAnaTestbeam::FillHistos: NULL Channel Pointer for ChId "
		   << Form(" 0x%08x ",iChId)
		   <<FairLogger::endl;
      }

      if(iSmType<fhXYPos.size())
      fhXYPos[iSmType]->Fill(pHit->GetX(),pHit->GetY());

      switch(iSmType){

      case 2 :   // Plastic hit
	if (0 != CbmTofAddress::GetSmId( iDetId ) - fiPlaSelect) break;
	// get here for Plastic SmId == 0 only! 
	for( Int_t iHitInd2 = iHitInd; iHitInd2 < iNbTofHits; iHitInd2++)
	{
	  pHit2 = (CbmTofHit*) fTofHitsColl->At( iHitInd2 );
	  if(pHit2==NULL) continue;
	  Int_t iDetId2 = (pHit2->GetAddress() & DetMask);
	  Int_t iChId2 = pHit2->GetAddress();
          fChannelInfo2 = fDigiPar->GetCell( iChId2 );
      
	  if( 2 == CbmTofAddress::GetSmType( iDetId2 )){   // Plastic hit correlation 
	   if (0==CbmTofAddress::GetSmId( iDetId ) - fiPlaSelect ){
	    if (1==CbmTofAddress::GetSmId( iDetId2 ) - fiPlaSelect) {
	      fhDT2->Fill(pHit->GetTime()-pHit2->GetTime());
	      fhXX2->Fill(pHit->GetX(),pHit2->GetX());
	      fhYY2->Fill(pHit->GetY(),pHit2->GetY());
	    }
	   }
	   else if (1==CbmTofAddress::GetSmId( iDetId ) - fiPlaSelect){ // should never be reached 
	     LOG(WARNING) << "CbmTofAnaTestbeam::FillHistos: inconsistent Plastic analysis, check!"
			  <<FairLogger::endl;
	    if (0==CbmTofAddress::GetSmId( iDetId2 )){
	      fhXX2->Fill(pHit2->GetX(),pHit->GetX());
	      fhYY2->Fill(pHit2->GetY(),pHit->GetY());
	    }	  
	   }
	  }

	  if( fiMrpcRef == CbmTofAddress::GetSmType( iDetId2 )){   // HdRef hit
	    Double_t xPos=Zref/pHit->GetZ()*pHit->GetX();
	    Double_t yPos=Zref/pHit->GetZ()*pHit->GetY();
	    Double_t tof=pHit->GetTime();
            Double_t xPos2=Zref/pHit2->GetZ()*pHit2->GetX();
	    Double_t yPos2=Zref/pHit2->GetZ()*pHit2->GetY();
            Double_t tof2=pHit2->GetTime();
            Double_t Chi2Match =TMath::Power((xPos-xPos2)/fdDXWidth,2.)
	                       +TMath::Power((yPos-yPos2)/fdDYWidth,2.)
	                       +TMath::Power((tof-tof2)/fdDTWidth,2.);

	    if(   TMath::Abs(xPos-xPos2)<10.
	       && TMath::Abs(yPos-yPos2)<10.
	       && TMath::Abs(tof-tof2)<2500.
	       && TMath::Abs(pHit2->GetY()-fChannelInfo2->GetY())<fChannelInfo2->GetSizey() )
	    {
	      iMat24++;
	      BSel[0]=kTRUE;                              // first selection scheme 
	      fhChi24->Fill(Chi2Match);
	      fhXY24->Fill(xPos2,yPos2);
	      fhDXDY24->Fill(xPos-xPos2,yPos-yPos2);
	      fhDXDT24->Fill(xPos-xPos2,tof-tof2);
	      fhDYDT24->Fill(yPos-yPos2,tof-tof2);
	      if(Chi2Match<Chi2Mat24Min){
		Chi2Mat24Min=Chi2Match;
		dDT24=tof-tof2;
	      }
	    }
	  }
	}
	break;

      case 0:   // found HdBig
      case 3:   // found Hd-P2 
      case 6:   // found Buc2013 
      case 8:   // found THU Pad 
	{
	if(fiDut != iSmType) continue;
	dMul0++;
	Double_t xPos=Zref/pHit->GetZ()*pHit->GetX();
	Double_t yPos=Zref/pHit->GetZ()*pHit->GetY();
        Double_t tof=pHit->GetTime();

	for( Int_t iHitInd2 = 0; iHitInd2 < iNbTofHits; iHitInd2++)
	{
          if(iHitInd2 != iHitInd){
	   pHit2 = (CbmTofHit*) fTofHitsColl->At( iHitInd2 );
	   if (NULL== pHit2) continue;

	   Int_t iDetId2 = (pHit2->GetAddress() & DetMask);
	   Int_t iChId2 = pHit2->GetAddress();
           fChannelInfo2 = fDigiPar->GetCell( iChId2 );

	   if (fiMrpcRef==CbmTofAddress::GetSmType( iDetId2 )){  //HdRef -> MrpcRef

             Double_t xPos2=Zref/pHit2->GetZ()*pHit2->GetX();
	     Double_t yPos2=Zref/pHit2->GetZ()*pHit2->GetY();
             Double_t tof2=pHit2->GetTime();
             Double_t Chi2Match =TMath::Power((xPos-xPos2-fdDXMean)/fdDXWidth,2.)
	                        +TMath::Power((yPos-yPos2-fdDYMean)/fdDYWidth,2.)
	                        +TMath::Power((tof-tof2-fdDTMean)/fdDTWidth,2.);

             if (Chi2Match > 1.E8) continue;

	     iNbMatchedHits++;
	     if(iNbMatchedHits==iNbMaxMatch) iNbMatchedHits=iNbMaxMatch-1; //prevent array overflow

 	     LOG(DEBUG)
	       <<Form("CbmTofAnaTestbeam::FillHistos: match %d:  %d.-%d. Tof hit 0x%08x of SmT %d, 0x%08x at 0x%08x,0x%08x",
		      iNbMatchedHits, iHitInd, iHitInd2, iDetId2, CbmTofAddress::GetSmType( iDetId2 ), iChId2, pHit, pHit2)
	       <<FairLogger::endl;

	     fhXX04->Fill(xPos,xPos2);
	     fhYY04->Fill(yPos,yPos2);
	     fhTT04->Fill(tof,tof2);
	     fhDXDY04->Fill(xPos-xPos2,yPos-yPos2);
	     fhDXDT04->Fill(xPos-xPos2,tof-tof2);
	     fhDYDT04->Fill(yPos-yPos2,tof-tof2);
             fhChi04->Fill(Chi2Match);

	     for(Int_t iM=0; iM<iNbMatchedHits; iM++)
	     {
              if(Chi2Match < Chi2List[iM]){

	       LOG(DEBUG)<<Form(" fill Chi2 %d (%d), %6.1f, %6.1f: HitInd %d, %d, Ids 0x%08x, 0x%08x",  
				iM, iNbMatchedHits, Chi2Match, Chi2List[iM], iHitInd,iHitInd2,pHit,pHit2)
			 <<FairLogger::endl;

               for(Int_t iMM=iNbMatchedHits; iMM>=iM; iMM--){
		 Chi2List[iMM]= Chi2List[iMM-1];
		 pChi2Hit1[iMM]= pChi2Hit1[iMM-1];
		 pChi2Hit2[iMM]= pChi2Hit2[iMM-1];
	       }

 	       Chi2List[iM]=Chi2Match;
	       pChi2Hit1[iM]=pHit;
	       pChi2Hit2[iM]=pHit2;
	       Chi2List[iNbMatchedHits]=1.E8;
               if(iM>0){
		 if(Chi2Match == Chi2List[iM-1]){
		   LOG(DEBUG)<<Form("same Chi2?: M %d, Mul0 %3.0f, HitInd %d, %d, Ids 0x%08x, 0x%08x - 0x%08x, 0x%08x",  
				   iM,dMul0,iHitInd,iHitInd2,pHit,pHit2,pChi2Hit1[iM-1],pChi2Hit2[iM-1])
			    <<FairLogger::endl;
		   LOG(DEBUG) << pHit->ToString()           <<FairLogger::endl;
		   LOG(DEBUG) << pChi2Hit1[iM-1]->ToString()<<FairLogger::endl;
		   LOG(DEBUG) << pHit2->ToString()          <<FairLogger::endl;
		 
		 }
	       }
	       break;
	      }
	     }
	   }
	   if (2==CbmTofAddress::GetSmType( iDetId2 )){  // Plastic 
	     Int_t iSm=CbmTofAddress::GetSmId( iDetId2 );
             if(iSm<2){
	      fhXX02[iSm]->Fill(Zref/pHit->GetZ()*pHit->GetX(),Zref/pHit2->GetZ()*pHit2->GetX());
	      fhYY02[iSm]->Fill(Zref/pHit->GetZ()*pHit->GetY(),Zref/pHit2->GetZ()*pHit2->GetY());
	     }
	   }
	  }
	}
	}
	break;

      case 1:
      case 4:
        if(fiMrpcRef != iSmType) break;
	dMul4++;
	break;

      case 5:
        if( 1 != CbmTofAddress::GetSmId( iDetId )) break;
	dMulD++;
	dTDia = pHit->GetTime();
	if( fdDTDia>0. )
	{
	  Double_t dDDia=0.;
	  for( Int_t iHitInd1 = 0; iHitInd1 < iNbTofHits; iHitInd1++) 
	  if( iHitInd1!=iHitInd)
	  {
	    pHit1 = (CbmTofHit*) fTofHitsColl->At( iHitInd1 );
	    if(pHit1==NULL) continue;
	    Int_t iDetId1 = (pHit1->GetAddress() & DetMask);
	    Int_t iChId1 = pHit1->GetAddress();
            fChannelInfo1 = fDigiPar->GetCell( iChId1 );
	    pHit1 = (CbmTofHit*) fTofHitsColl->At( iHitInd1 );
	    if(  5 == CbmTofAddress::GetSmType( iDetId1 ) 
	      && 1 != CbmTofAddress::GetSmId  ( iDetId1 )
	      && TMath::Abs( pHit1->GetTime()-dTDia)<fdDTDia){   // second diamond fired
	      dDDia= pHit1->GetTime()-dTDia;
	    }
	  }
	  LOG(DEBUG)<<Form("CbmTofAnaTestbeam:FillHisto: DDia %f",dDDia)<<FairLogger::endl; 
	  if (dDDia==0.) break;
	}

	for( Int_t iHitInd2 = 0; iHitInd2 < iNbTofHits; iHitInd2++) 
	  if(iHitInd2!=iHitInd)
	  {
	   pHit2 = (CbmTofHit*) fTofHitsColl->At( iHitInd2 );
	   if(pHit2==NULL) continue;

	   Int_t iDetId2 = (pHit2->GetAddress() & DetMask);
	   Int_t iChId2 = pHit2->GetAddress();
	   fChannelInfo2 = fDigiPar->GetCell( iChId2 );
	   if(NULL == fChannelInfo2){
	     LOG(ERROR) << "CbmTofAnaTestbeam::FillHistos: Invalid Channel Pointer for ChId2 "
			<< Form(" 0x%08x ",iChId2)
			<<FairLogger::endl;
	     continue;
	   }
	   LOG(DEBUG2)<<Form("CbmTofAnaTestbeam:FillHisto: TDia %f, THit %f",dTDia,pHit2->GetTime())
		      <<FairLogger::endl;
	   if( fiMrpcRef == CbmTofAddress::GetSmType( iDetId2 )){   // Reference RPC hit
             pHitRef=pHit2;
	     dDTD4=pHit2->GetTime()-dTDia;
	     LOG(DEBUG1)<<Form("CbmTofAnaTestbeam:FillHisto: dDTD4 %f, min: %f",dDTD4,dDTD4Min)
			<<FairLogger::endl;

	     if(TMath::Abs(dDTD4)<DTDMAX                            // second selection scheme, selector 1
	      &&TMath::Abs(pHit2->GetY()-fChannelInfo2->GetY())<fdPosY4Sel*fChannelInfo2->GetSizey()
	      &&TMath::Abs(CbmTofAddress::GetChannelId( iChId2 ) - fdCh4Sel) < fdDCh4Sel
	      //	       && CbmTofAddress::GetChannelId( iChId2 ) < 15
	       )
	       {
		 BSel[1]=kTRUE; 
		 if(TMath::Abs(dDTD4)<TMath::Abs(dDTD4Min)){
		   dDTD4Min=dDTD4;
		 }
	       }
	    }
	   }	//LOG(INFO)<<" TDia="<<dTDia<<FairLogger::endl;
      break;
      default:
	;
      }
   } // for( Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++)

   if(dMul4>dM4Max || dMulD>dMDMax || dMul0>dM0Max)
     {
       BSel[0]=kFALSE;
       BSel[1]=kFALSE;
     }
   //  normalisation distributions 
   fhNMatch24->Fill(iMat24);
   fhNMatch04->Fill(iNbMatchedHits);
   LOG(DEBUG)<<Form(" Matches: %d with first chi2s = %6.1f, %6.1f, %6.1f",iNbMatchedHits,
		   Chi2List[0],Chi2List[1],Chi2List[2])
 	     <<FairLogger::endl;
   // selector 0 distributions 
   if(BSel[0]){
    fhNMatch24sel->Fill(iNbMatchedHits);  //use as normalisation
    fhDT24sel ->Fill(dDT24);              //  refined normalisation 

    if(iNbMatchedHits>0){
     // best match
     LOG(DEBUG)<<Form(" best match (%d): 0x%08x, 0x%8x, %f",iNbMatchedHits,pChi2Hit1[0],pChi2Hit2[0], Chi2List[0])
	      <<FairLogger::endl;

     pHit1=pChi2Hit1[0];
     pHit2=pChi2Hit2[0];

     Int_t iDetId1 = (pHit1->GetAddress() & DetMask);
     Int_t iChId1 = pHit1->GetAddress();
     fChannelInfo1 = fDigiPar->GetCell( iChId1 );

     Int_t iDetId2 = (pHit2->GetAddress() & DetMask);
     Int_t iChId2 = pHit2->GetAddress();
     fChannelInfo2 = fDigiPar->GetCell( iChId2 );

     Double_t xPos=Zref/pChi2Hit1[0]->GetZ()*pChi2Hit1[0]->GetX();
     Double_t yPos=Zref/pChi2Hit1[0]->GetZ()*pChi2Hit1[0]->GetY();
     Double_t tof=pChi2Hit1[0]->GetTime();

     Double_t xPos2=Zref/pChi2Hit2[0]->GetZ()*pChi2Hit2[0]->GetX();
     Double_t yPos2=Zref/pChi2Hit2[0]->GetZ()*pChi2Hit2[0]->GetY();
     Double_t tof2=pChi2Hit2[0]->GetTime();

     fhDXDY04best->Fill(xPos-xPos2,yPos-yPos2);
     fhDXDT04best->Fill(xPos-xPos2,tof-tof2);
     fhDYDT04best->Fill(yPos-yPos2,tof-tof2);
     fhChi04best->Fill(Chi2List[0]);

     fhX0DT04best->Fill(xPos-fChannelInfo1->GetX(),tof-tof2);
     fhY0DT04best->Fill(yPos-fChannelInfo1->GetY(),tof-tof2);

     CbmMatch* digiMatch1=(CbmMatch *)fTofDigiMatchColl->At(pHit1->GetRefId());
     fhDigiMul0best->Fill(digiMatch1->GetNofLinks()/2.);

     fhChiDT04best->Fill(Chi2List[0],tof-tof2);
     fhDT24DT04best->Fill(dDT24,tof-tof2);
     Double_t dTcor=0.;
     if(fhDTD4DT04D4Off != NULL) dTcor=(Double_t)fhDTD4DT04D4Off->GetBinContent(fhDTD4DT04D4Off->FindBin(-dDTD4));
     fhDTD4DT04best->Fill(-dDTD4,tof-tof2-dTcor);
    }
   }

   // selector 1 distributions 
   if(BSel[1]){
    /* 
    for( Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++)
    {
      pHit = (CbmTofHit*) fTofHitsColl->At( iHitInd );
      Int_t iDetId = (pHit->GetAddress() & DetMask);
      Int_t iChId = pHit->GetAddress();
      fChannelInfo = fDigiPar->GetCell( iChId );
      Int_t iSmType=CbmTofAddress::GetSmType( iDetId );

    }
    */
    fhNMatchD4sel->Fill(iNbMatchedHits);  // use as normalisation
    fhDTD4sel->Fill(-dDTD4Min);            //  general normalisation  

    if(iNbMatchedHits>0){
     // best match
     LOG(DEBUG)<<Form(" best match (%d): 0x%08x, 0x%8x, %f",iNbMatchedHits,pChi2Hit1[0],pChi2Hit2[0], Chi2List[0])
	      <<FairLogger::endl;

     pHit1=pChi2Hit1[0];
     pHit2=pChi2Hit2[0];

     Int_t iDetId1 = (pHit1->GetAddress() & DetMask);
     Int_t iChId1 = pHit1->GetAddress();
     fChannelInfo1 = fDigiPar->GetCell( iChId1 );

     Int_t iDetId2 = (pHit2->GetAddress() & DetMask);
     Int_t iChId2 = pHit2->GetAddress();
     fChannelInfo2 = fDigiPar->GetCell( iChId2 );

     Double_t xPos=Zref/pChi2Hit1[0]->GetZ()*pChi2Hit1[0]->GetX();
     Double_t yPos=Zref/pChi2Hit1[0]->GetZ()*pChi2Hit1[0]->GetY();
     Double_t tof=pChi2Hit1[0]->GetTime();

     Double_t xPos2=Zref/pChi2Hit2[0]->GetZ()*pChi2Hit2[0]->GetX();
     Double_t yPos2=Zref/pChi2Hit2[0]->GetZ()*pChi2Hit2[0]->GetY();
     Double_t tof2=pChi2Hit2[0]->GetTime();

     Double_t dDist=TMath::Sqrt(
				TMath::Power(pHit1->GetX()-pHit2->GetX(),2)
			       +TMath::Power(pHit1->GetY()-pHit2->GetY(),2)
			       +TMath::Power(pHit1->GetZ()-pHit2->GetZ(),2)
				);

     fhChi04D4best->Fill(Chi2List[0]);

     CbmMatch* digiMatch1=(CbmMatch *)fTofDigiMatchColl->At(pHit1->GetRefId());
     fhDigiMul0D4best->Fill(digiMatch1->GetNofLinks()/2.);
     fhDigiMul4D4best->Fill(digiMatch1->GetNofLinks()/2.);

     fhCluMul04D4best->Fill(dMul0,dMul4);

     // check for dependence in counter reference frame
     TGeoNode *fNode=        // prepare global->local trafo
     gGeoManager->FindNode(fChannelInfo2->GetX(),fChannelInfo2->GetY(),fChannelInfo2->GetZ());

     hitpos[0]=pChi2Hit2[0]->GetX();
     hitpos[1]=pChi2Hit2[0]->GetY();
     hitpos[2]=pChi2Hit2[0]->GetZ();
     TGeoNode* cNode= gGeoManager->GetCurrentNode();
     gGeoManager->MasterToLocal(hitpos, hitpos_local);

     Double_t dTofD4  = fdTOffD4 + dDTD4;
     Double_t dInvVel = dTofD4/pHitRef->GetR(); // in ps/cm
     Double_t dDTexp  = dDist*dInvVel;

     Double_t dTcor=0.;
     if(fhDTD4DT04D4Off != NULL) dTcor=(Double_t)fhDTD4DT04D4Off->GetBinContent(fhDTD4DT04D4Off->FindBin(-dDTD4));
     if(fhDTX4D4Off     != NULL) dTcor+=(Double_t)fhDTX4D4Off->GetBinContent(fhDTX4D4Off->FindBin(hitpos_local[0]));
     if(fhDTY4D4Off     != NULL) dTcor+=(Double_t)fhDTY4D4Off->GetBinContent(fhDTY4D4Off->FindBin(hitpos_local[1]));
     if(fhDTTexpD4Off   != NULL) dTcor+=(Double_t)fhDTTexpD4Off->GetBinContent(fhDTTexpD4Off->FindBin(dDTexp));

     //     LOG(INFO) << "dTcor for "<<-dDTD4<<" from "<<fhDTD4DT04D4Off<<": "<<dTcor<<FairLogger::endl;



     fhTofD4best->Fill(dTofD4);
     if(dInvVel>0.) fhVelD4best->Fill(1000./dInvVel);

     fhChiDT04D4best->Fill(Chi2List[0],tof-tof2-dTcor);
     if(dDT24 != 0.) fhDT24DT04D4best->Fill(dDT24,tof-tof2-dTcor);

     fhDTD4DT04D4best->Fill(-dDTD4,tof-tof2-dTcor);
     fhDTMul4D4best->Fill(dMul4,tof-tof2-dTcor);

     fhDTX4D4best->Fill(hitpos_local[0],tof-tof2-dTcor);
     fhDTY4D4best->Fill(hitpos_local[1],tof-tof2-dTcor);
     fhDXX4D4best->Fill(hitpos_local[0],xPos-xPos2);
     fhDXY4D4best->Fill(hitpos_local[1],xPos-xPos2);
     fhDYX4D4best->Fill(hitpos_local[0],yPos-yPos2);
     fhDYY4D4best->Fill(hitpos_local[1],yPos-yPos2);

     CbmMatch* digiMatch4=(CbmMatch *)fTofDigiMatchColl->At(pChi2Hit2[0]->GetRefId());
     fhCluSize4DT04D4best->Fill(digiMatch4->GetNofLinks()/2.,tof-tof2-dTcor);
     Double_t dTot = 0.;
     for (Int_t iLink=0; iLink<digiMatch4->GetNofLinks(); iLink++){  // loop over digis
         CbmLink L0 = digiMatch4->GetLink(iLink);  
	 Int_t iDigInd0=L0.GetIndex();
         if (iDigInd0 < fTofDigisColl->GetEntries()){
            CbmTofDigiExp *pDig0 = (CbmTofDigiExp*) (fTofDigisColl->At(iDigInd0));
	    dTot += pDig0->GetTot();
	 }
     } 
     dTot /= digiMatch4->GetNofLinks();  // average time over threshold
     fhTot4DT04D4best->Fill(dTot,tof-tof2-dTcor);

     CbmMatch* digiMatch0=(CbmMatch *)fTofDigiMatchColl->At(pChi2Hit1[0]->GetRefId());
     fhCluSize0DT04D4best->Fill(digiMatch0->GetNofLinks()/2.,tof-tof2-dTcor);
     dTot = 0.;
     for (Int_t iLink=0; iLink<digiMatch0->GetNofLinks(); iLink++){  // loop over digis
         CbmLink L0 = digiMatch0->GetLink(iLink);  
	 Int_t iDigInd0=L0.GetIndex();
         if (iDigInd0 < fTofDigisColl->GetEntries()){
            CbmTofDigiExp *pDig0 = (CbmTofDigiExp*) (fTofDigisColl->At(iDigInd0));
	    dTot += pDig0->GetTot();
	 }
     } 
     dTot /= digiMatch0->GetNofLinks();  // average time over threshold
     fhTot0DT04D4best->Fill(dTot,tof-tof2-dTcor);

     fhDTMul0D4best->Fill(dMul0,tof-tof2-dTcor);


     // check for dependence in counter reference frame
     TGeoNode *fNode1=        // prepare global->local trafo
     gGeoManager->FindNode(fChannelInfo1->GetX(),fChannelInfo1->GetY(),fChannelInfo1->GetZ());

     hitpos[0]=pChi2Hit1[0]->GetX();
     hitpos[1]=pChi2Hit1[0]->GetY();
     hitpos[2]=pChi2Hit1[0]->GetZ();

     TGeoNode* cNode1= gGeoManager->GetCurrentNode();
     gGeoManager->MasterToLocal(hitpos, hitpos_local);

     fhDTX0D4best->Fill(hitpos_local[0],tof-tof2-dTcor);
     fhDTY0D4best->Fill(hitpos_local[1],tof-tof2-dTcor);
     fhDXX0D4best->Fill(hitpos_local[0],xPos-xPos2);
     fhDXY0D4best->Fill(hitpos_local[1],xPos-xPos2);
     fhDYX0D4best->Fill(hitpos_local[0],yPos-yPos2);
     fhDYY0D4best->Fill(hitpos_local[1],yPos-yPos2);

     fhDXDY04D4best->Fill(xPos-xPos2,yPos-yPos2);
     fhDXDT04D4best->Fill(xPos-xPos2,tof-tof2-dTcor);
     fhDYDT04D4best->Fill(yPos-yPos2,tof-tof2-dTcor);
     fhDistDT04D4best->Fill(dDist,tof-tof2-dTcor);
     fhTexpDT04D4best->Fill(dDTexp,tof-tof2-dTcor);
     fhX0DT04D4best->Fill(hitpos_local[0],tof-tof2-dTcor);
     fhY0DT04D4best->Fill(hitpos_local[1],tof-tof2)-dTcor;

    }
   }
   return kTRUE;
}
// ------------------------------------------------------------------

Bool_t CbmTofAnaTestbeam::WriteHistos()
{

   LOG(INFO)<<"CbmTofAnaTestbeam::WriteHistos: ./tofAnaTestBeam.hst.root, mode = "
	    << fiCorMode
	    << FairLogger::endl;

   // Write histogramms to the file
   TDirectory * oldir = gDirectory;
   TFile *fHist = new TFile("./tofAnaTestBeam.hst.root","RECREATE");
   fHist->cd();

   switch(fiCorMode){
   case 1 : {
     //   fhDT24DT04D4best->ProfileX();
     //   fhDT24DT04D4best->Write();
     TProfile *htmp=fhDTD4DT04D4best->ProfileX();
     TH1D *htmp1D=htmp->ProjectionX();

     //     htmp1D->Draw();

     if(fhDTD4DT04D4Off != NULL){
       // fhDTD4DT04D4Off->Draw("same");
       // fhDTD4DT04D4Off->Write();
       //LOG(INFO)<<"Update hDTD4DT04D4best"<<FairLogger::endl;
       Double_t nx=htmp1D->GetNbinsX();
       for (Int_t ix=0; ix<nx; ix++){
	 Double_t dVal=htmp1D->GetBinContent(ix) + fhDTD4DT04D4Off->GetBinContent(ix);
	 // Double_t dVal=fhDTD4DT04D4Off->GetBinContent(ix);
	 LOG(DEBUG1)<<"Update hDTD4DT04D4best "<<ix<<": "<<htmp1D->GetBinContent(ix)<<" + "
		  << fhDTD4DT04D4Off->GetBinContent(ix) << " -> " << dVal << FairLogger::endl;
	 htmp1D->SetBinContent(ix,dVal);
       }
     }else
     {
       LOG(WARNING)<<"CbmTofAnaTestbeam::WriteHistos: fhDTD4DT04D4Off not found "
		   << FairLogger::endl;
     }
   //   fhDTD4DT04D4best->Write();
     htmp1D->Write();
     if(fhDTX4D4Off != NULL)    fhDTX4D4Off->Write();
     if(fhDTY4D4Off != NULL)    fhDTY4D4Off->Write();
     if(fhDTTexpD4Off != NULL)  fhDTTexpD4Off->Write();
     }
     break;
   case 2 :{
     TProfile *htmpx=fhDTX4D4best->ProfileX();
     TH1D *htmpx1D=htmpx->ProjectionX();
     if(fhDTX4D4Off != NULL){
       Double_t nx=htmpx1D->GetNbinsX();
       for (Int_t ix=0; ix<nx; ix++){
	 Double_t dVal=htmpx1D->GetBinContent(ix) + fhDTX4D4Off->GetBinContent(ix);
	 LOG(DEBUG1)<<"Update hDTX4D4best "<<ix<<": "<<htmpx1D->GetBinContent(ix)<<" + "
		    << fhDTX4D4Off->GetBinContent(ix) << " -> " << dVal << FairLogger::endl;
	 htmpx1D->SetBinContent(ix,dVal);
       }
     }else
     {
       LOG(WARNING)<<"CbmTofAnaTestbeam::WriteHistos: fhDTX4D4Off not found "
		   << FairLogger::endl;
     }
     htmpx1D->Write();
     if(fhDTD4DT04D4Off != NULL) fhDTD4DT04D4Off->Write();
     if(fhDTY4D4Off != NULL)     fhDTY4D4Off->Write();
     if(fhDTTexpD4Off != NULL)   fhDTTexpD4Off->Write();
     }
     break;
   case 3 :{
     TProfile *htmpx=fhDTY4D4best->ProfileX();
     TH1D *htmpx1D=htmpx->ProjectionX();
     if(fhDTY4D4Off != NULL){
       Double_t nx=htmpx1D->GetNbinsX();
       for (Int_t ix=0; ix<nx; ix++){
	 Double_t dVal=htmpx1D->GetBinContent(ix) + fhDTY4D4Off->GetBinContent(ix);
	 LOG(DEBUG1)<<"Update hDTY4D4best "<<ix<<": "<<htmpx1D->GetBinContent(ix)<<" + "
		    << fhDTY4D4Off->GetBinContent(ix) << " -> " << dVal << FairLogger::endl;
	 htmpx1D->SetBinContent(ix,dVal);
       }
     }else
     {
       LOG(WARNING)<<"CbmTofAnaTestbeam::WriteHistos: fhDTY4D4Off not found "
		   << FairLogger::endl;
     }
     htmpx1D->Write();
     if(fhDTD4DT04D4Off != NULL) fhDTD4DT04D4Off->Write();
     if(fhDTX4D4Off != NULL)     fhDTX4D4Off->Write();
     if(fhDTTexpD4Off != NULL)   fhDTTexpD4Off->Write();
     }
     break;

   case 4 :{
     TProfile *htmpx=fhTexpDT04D4best->ProfileX();
     TH1D *htmpx1D=htmpx->ProjectionX();
     if(fhDTTexpD4Off != NULL){
       Double_t nx=htmpx1D->GetNbinsX();
       for (Int_t ix=0; ix<nx; ix++){
	 Double_t dVal=htmpx1D->GetBinContent(ix) + fhDTTexpD4Off->GetBinContent(ix);
	 LOG(DEBUG1)<<"Update hDTTexpD4best "<<ix<<": "<<htmpx1D->GetBinContent(ix)<<" + "
		    << fhDTTexpD4Off->GetBinContent(ix) << " -> " << dVal << FairLogger::endl;
	 htmpx1D->SetBinContent(ix,dVal);
       }
     }else
     {
       LOG(WARNING)<<"CbmTofAnaTestbeam::WriteHistos: fhDTTexpD4Off not found "
		   << FairLogger::endl;
     }
     htmpx1D->Write();
     if(fhDTD4DT04D4Off != NULL) fhDTD4DT04D4Off->Write();
     if(fhDTX4D4Off != NULL)     fhDTX4D4Off->Write();
     if(fhDTY4D4Off != NULL)     fhDTY4D4Off->Write();
     }
     break;

   default:
     ;
   }

   //   fHist->Write();
   if(0) {
   fhXX2->Write();
   fhYY2->Write();
   for (Int_t iDet=0; iDet<2; iDet++)
   {
     fhXX02[iDet]->Write();
     fhYY02[iDet]->Write(); 
   }
   fhXX04->Write();
   fhYY04->Write();
   }
   gDirectory->cd( oldir->GetPath() );

   fHist->Close();
   return kTRUE;
}

Bool_t   CbmTofAnaTestbeam::DeleteHistos()
{
   // Test class performance

   // Mapping


   return kTRUE;
}

ClassImp(CbmTofAnaTestbeam);