/********************************************************************************
 *    Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH    *
 *                                                                              *
 *              This software is distributed under the terms of the             *
 *              GNU Lesser General Public Licence (LGPL) version 3,             *
 *                  copied verbatim in the file "LICENSE"                       *
 ********************************************************************************/
#include "CbmMcbm2019CheckPulser.h"

#include "CbmStsDigi.h"
#include "CbmTofDigiExp.h"
#include "CbmMuchBeamTimeDigi.h"

#include "FairLogger.h"
#include "FairRootManager.h"
#include "FairRunOnline.h"

#include "TClonesArray.h"
#include "TH1.h"
#include "TH2.h"
#include "TProfile.h"
#include "THttpServer.h"
#include <TFile.h>

#include <iomanip>
#include <iostream>
using std::fixed;
using std::setprecision;

// ---- Default constructor -------------------------------------------
CbmMcbm2019CheckPulser::CbmMcbm2019CheckPulser()
  : FairTask("CbmMcbm2019CheckPulser")
{
}

// ---- Destructor ----------------------------------------------------
CbmMcbm2019CheckPulser::~CbmMcbm2019CheckPulser()
{
}

// ----  Initialisation  ----------------------------------------------
void CbmMcbm2019CheckPulser::SetParContainers()
{
  // Load all necessary parameter containers from the runtime data base
  /*
  FairRunAna* ana = FairRunAna::Instance();
  FairRuntimeDb* rtdb=ana->GetRuntimeDb();

  <CbmMcbm2019CheckPulserDataMember> = (<ClassPointer>*)
    (rtdb->getContainer("<ContainerName>"));
  */
}

// ---- Init ----------------------------------------------------------
InitStatus CbmMcbm2019CheckPulser::Init()
{

  // Get a handle from the IO manager
  FairRootManager* ioman = FairRootManager::Instance();

  // Get a pointer to the previous already existing data level

  fT0Digis = static_cast<TClonesArray*>(ioman->GetObject("CbmT0Digi"));
  if ( ! fT0Digis ) {
    LOG(fatal) << "No TClonesArray with T0 digis found.";
  }

  fStsDigis = static_cast<TClonesArray*>(ioman->GetObject("CbmStsDigi"));
  if ( ! fStsDigis ) {
    LOG(info) << "No TClonesArray with STS digis found.";
  }

  fMuchDigis = static_cast<TClonesArray*>(ioman->GetObject("CbmMuchBeamTimeDigi"));
  if ( ! fMuchDigis ) {
    LOG(info) << "No TClonesArray with MUCH digis found.";
  }

  fTofDigis = static_cast<TClonesArray*>(ioman->GetObject("CbmTofDigi"));
  if ( ! fTofDigis ) {
    LOG(info) << "No TClonesArray with TOF digis found.";
  }

  fRichDigis = static_cast<TClonesArray*>(ioman->GetObject("CbmRichDigi"));
  if ( ! fRichDigis ) {
    LOG(info) << "No TClonesArray with RICH digis found.";
  }

  fPsdDigis = static_cast<TClonesArray*>(ioman->GetObject("CbmPsdDigi"));
  if ( ! fPsdDigis ) {
    LOG(info) << "No TClonesArray with PSD digis found.";
  }

  CreateHistos();

  return kSUCCESS;
}

Int_t CbmMcbm2019CheckPulser::CalcNrBins(Int_t offsetRange) {

  if (offsetRange < 1001) {
    fBinWidth = 5;
  } else if (offsetRange < 10001) {
    fBinWidth = 10;
  } else if (offsetRange < 100001) {
    fBinWidth = 100;
  } else {
    fBinWidth = 100;
  }

  return (offsetRange/fBinWidth *2);
}

void CbmMcbm2019CheckPulser::CreateHistos()
{
   /// Logarithmic bining for self time comparison
   // Number of log bins =
   //      9 for the sub-unit decade
   //    + 9 for each unit of each decade * 10 for the subdecade range
   //    + 1 for the closing bin top edge
   const Int_t iNbDecadesLog     = 9;
   const Int_t iNbStepsDecade    = 9;
   const Int_t iNbSubStepsInStep = 1;
   const Int_t iNbBinsLog = iNbStepsDecade
                           + iNbStepsDecade * iNbSubStepsInStep * iNbDecadesLog
                           + 1;
   Double_t dBinsLog[iNbBinsLog];
      // First fill sub-unit decade
   for( Int_t iSubU = 0; iSubU < iNbStepsDecade; iSubU ++ )
      dBinsLog[ iSubU ] = 0.1 * ( 1 + iSubU );
   std::cout << std::endl;
      // Then fill the main decades
   Double_t dSubstepSize = 1.0 / iNbSubStepsInStep;
   for( Int_t iDecade = 0; iDecade < iNbDecadesLog; iDecade ++)
   {
      Double_t dBase = std::pow( 10, iDecade );
      Int_t iDecadeIdx = iNbStepsDecade
                       + iDecade * iNbStepsDecade * iNbSubStepsInStep;
      for( Int_t iStep = 0; iStep < iNbStepsDecade; iStep++ )
      {
         Int_t iStepIdx = iDecadeIdx + iStep * iNbSubStepsInStep;
         for( Int_t iSubStep = 0; iSubStep < iNbSubStepsInStep; iSubStep++ )
         {
            dBinsLog[ iStepIdx + iSubStep ] = dBase * (1 + iStep)
                                             + dBase * dSubstepSize * iSubStep;
         } // for( Int_t iSubStep = 0; iSubStep < iNbSubStepsInStep; iSubStep++ )
      } // for( Int_t iStep = 0; iStep < iNbStepsDecade; iStep++ )
   } // for( Int_t iDecade = 0; iDecade < iNbDecadesLog; iDecade ++)
   dBinsLog[ iNbBinsLog - 1 ] = std::pow( 10, iNbDecadesLog );

  Int_t nrOfBinsSts =  CalcNrBins(fStsOffsetRange);
  // T0 vs. Sts
  fT0StsDiff = new TH1F("fT0StsDiff","T0-Sts;time diff [ns];Counts",
			nrOfBinsSts, -fStsOffsetRange, fStsOffsetRange);

  fT0StsDiffEvo = new TH2F("fT0StsDiffEvo","T0-Sts;TS; time diff [ns];Counts",
       1000, 0, 10000,
			nrOfBinsSts, -fStsOffsetRange, fStsOffsetRange);


  Int_t nrOfBinsMuch =  CalcNrBins(fMuchOffsetRange);
  // T0 vs. Much
  fT0MuchDiff = new TH1F("fT0MuchDiff","T0-Much;time diff [ns];Counts",
			 nrOfBinsMuch, -fMuchOffsetRange, fMuchOffsetRange);

  fT0MuchDiffEvo = new TH2F("fT0MuchDiffEvo","T0-Much;TS; time diff [ns];Counts",
       1000, 0, 10000,
			 nrOfBinsMuch, -fMuchOffsetRange, fMuchOffsetRange);


  Int_t nrOfBinsTof =  CalcNrBins(fTofOffsetRange);
  // To vs. Tof
  fT0TofDiff = new TH1F("fT0TofDiff","T0-Tof;time diff [ns];Counts",
			nrOfBinsTof, -fTofOffsetRange, fTofOffsetRange);

  fT0TofDiffEvo = new TH2F("fT0TofDiffEvo","T0-Tof;TS; time diff [ns];Counts",
       1000, 0, 10000,
			nrOfBinsTof, -fTofOffsetRange, fTofOffsetRange);


  Int_t nrOfBinsRich =  CalcNrBins(fRichOffsetRange);
  // To vs. Rich
  fT0RichDiff = new TH1F("fT0RichDiff","T0-Rich;time diff [ns];Counts",
			nrOfBinsRich, -fRichOffsetRange, fRichOffsetRange);

  fT0RichDiffEvo = new TH2F("fT0RichDiffEvo","T0-Rich;TS; time diff [ns];Counts",
       1000, 0, 10000,
			nrOfBinsRich, -fRichOffsetRange, fRichOffsetRange);

  Int_t nrOfBinsPsd =  CalcNrBins(fPsdOffsetRange);
  // To vs. Rich
  fT0PsdDiff = new TH1F("fT0PsdDiff","T0-Psd;time diff [ns];Counts",
			nrOfBinsPsd, -fPsdOffsetRange, fPsdOffsetRange);

  fT0PsdDiffEvo = new TH2F("fT0PsdDiffEvo","T0-Psd;TS; time diff [ns];Counts",
       1000, 0, 10000,
			nrOfBinsPsd, -fPsdOffsetRange, fPsdOffsetRange);

  fT0PsdDiffCharge = new TH2F("fT0PsdDiffCharge","T0-Psd;time diff [ns]; Charge [a.u]; ;Counts",
			nrOfBinsPsd, -fPsdOffsetRange, fPsdOffsetRange,
      7000, 0, 70000
      );


  // T0 vs. Sts
  fT0StsDiffEvoLong = new TH2F("fT0StsDiffEvoLong","T0-Sts;TS; time diff [ns];Counts",
       1800, 0, 180000,
			nrOfBinsSts, -fStsOffsetRange, fStsOffsetRange);
  // T0 vs. Much
  fT0MuchDiffEvoLong = new TH2F("fT0MuchDiffEvoLong","T0-Much;TS; time diff [ns];Counts",
       1800, 0, 180000,
			 nrOfBinsMuch, -fMuchOffsetRange, fMuchOffsetRange);
  // To vs. Tof
  fT0TofDiffEvoLong = new TH2F("fT0TofDiffEvoLong","T0-Tof;TS; time diff [ns];Counts",
       1800, 0, 180000,
			nrOfBinsTof, -fTofOffsetRange, fTofOffsetRange);
  // To vs. Rich
  fT0RichDiffEvoLong = new TH2F("fT0RichDiffEvoLong","T0-Rich;TS; time diff [ns];Counts",
       1800, 0, 180000,
			nrOfBinsRich, -fRichOffsetRange, fRichOffsetRange);

  // To vs. Psd
  fT0PsdDiffEvoLong = new TH2F("fT0PsdDiffEvoLong","T0-Psd;TS; time diff [ns];Counts",
       1800, 0, 180000,
			nrOfBinsPsd, -fPsdOffsetRange, fPsdOffsetRange);


  // T0 vs. Sts
  fT0StsMeanEvo = new TProfile("fT0StsMeanEvo","T0-Sts; time in run [s]; Mean time diff [ns]",
       4320, 0, 4320 );
  // T0 vs. Much
  fT0MuchMeanEvo = new TProfile("fT0MuchMeanEvo","T0-Much; time in run [s]; Mean time diff [ns]",
       4320, 0, 4320 );
  // To vs. Tof
  fT0TofMeanEvo = new TProfile("fT0TofMeanEvo","T0-Tof; time in run [s]; Mean time diff [ns]",
       4320, 0, 4320 );
  // To vs. Rich
  fT0RichMeanEvo = new TProfile("fT0RichMeanEvo","T0-Rich; time in run [s]; Mean time diff [ns]",
       4320, 0, 4320);
  // To vs. Psd
  fT0PsdMeanEvo = new TProfile("fT0PsdMeanEvo","T0-Psd; time in run [s]; Mean time diff [ns]",
       4320, 0, 4320);
//       4320, 0, 259200);


  // T0 vs. STS for the different DPBs
  fT0StsDpbDiff = new TH2F("fT0StsDpbDiff","T0-Much;DPB; time diff [ns];Counts",
       2, -0.5, 1.5, nrOfBinsSts, -fStsOffsetRange, fStsOffsetRange);

  for( UInt_t uStsDpb = 0; uStsDpb < kuMaxNbStsDpbs; ++uStsDpb )
  {
    fT0StsDpbDiffEvo[uStsDpb] = new TH2F( Form("fT0StsDpbDiffEvo%02u", uStsDpb),
      Form( "T0-STS DPB %02u;TS; time diff [ns];Counts", uStsDpb),
       1800, 0, 180000,
			nrOfBinsSts, -fStsOffsetRange, fStsOffsetRange);
    fStsDpbCntsEvo[uStsDpb] = new TH1F( Form("fStsDpbCntsEvo%02u", uStsDpb),
      Form( "Time STS DPB %02u;TS; Hit Counts", uStsDpb),
       1800, 0, 180000);
  } // for( UInt_t uStsDpb = 0; uStsDpb < kuMaxNbStsDpbs; ++uStsDpb )

  // T0 vs. Much for the different DPBs/AFCK
  fT0MuchRocDiff = new TH2F("fT0MuchRocDiff","T0-Much;AFCK; time diff [ns];Counts",
       kuMaxNbMuchDpbs, -0.5, kuMaxNbMuchDpbs - 0.5, nrOfBinsMuch, -fMuchOffsetRange, fMuchOffsetRange);

  // T0 vs. Much for the different ASICs
  fT0MuchAsicDiff = new TH2F("fT0MuchAsicDiff","T0-Much;ASIC; time diff [ns];Counts",
       kuMaxNbMuchAsics, -0.5, kuMaxNbMuchAsics - 0.5, nrOfBinsMuch, -fMuchOffsetRange, fMuchOffsetRange);

  for( UInt_t uMuchAsic = 0; uMuchAsic < kuMaxNbMuchAsics; ++uMuchAsic )
    fT0MuchAsicDiffEvo[uMuchAsic] = new TH2F( Form("fT0MuchAsicDiffEvo%02u", uMuchAsic),
      Form( "T0-Much ASIC %02u;TS; time diff [ns];Counts", uMuchAsic),
       1800, 0, 180000,
			nrOfBinsMuch, -fMuchOffsetRange, fMuchOffsetRange);

  // T0 vs. T0
  fT0T0Diff = new TH1F("fT0T0Diff","T0-T0_prev;time diff [ns];Counts",
		       iNbBinsLog - 1, dBinsLog);
  // sts vs. Sts
  fStsStsDiff = new TH1F("fStsStsDiff","Sts-Sts_prev;time diff [ns];Counts",
		       iNbBinsLog - 1, dBinsLog);
  // Much vs. Much
  fMuchMuchDiff = new TH1F("fMuchMuchDiff","Much-Much_prev;time diff [ns];Counts",
		       iNbBinsLog - 1, dBinsLog);
  // Tof vs. Tof
  fTofTofDiff = new TH1F("fTofTofDiff","Tof-Tof_prev;time diff [ns];Counts",
		       iNbBinsLog - 1, dBinsLog);
  // Rich vs. Rich
  fRichRichDiff = new TH1F("fRichRichDiff","Rich-Rich_prev;time diff [ns];Counts",
		       iNbBinsLog - 1, dBinsLog);
  // Psd vs. Psd
  fPsdPsdDiff = new TH1F("fPsdPsdDiff","Psd-Psd_prev;time diff [ns];Counts",
		       iNbBinsLog - 1, dBinsLog);

  fT0Address = new TH1F("fT0Address","T0 address;address;Counts",
                        1000000, 0, 1000000.);

  fT0Channel = new TH1F("fT0Channel","T0 channel;channel nr;Counts",
                        100, -0.5, 99.5);

  /// Register the histos in the HTTP server
  FairRunOnline* run = FairRunOnline::Instance();
  if (run) {
    THttpServer* server = run->GetHttpServer();
    if( nullptr != server ) {
      server->Register("/CheckTiming", fT0StsDiff);
      server->Register("/CheckTiming", fT0MuchDiff);
      server->Register("/CheckTiming", fT0TofDiff);
      server->Register("/CheckTiming", fT0RichDiff);
      server->Register("/CheckTiming", fT0PsdDiff);
      server->Register("/CheckTiming", fT0StsDiffEvo);
      server->Register("/CheckTiming", fT0MuchDiffEvo);
      server->Register("/CheckTiming", fT0TofDiffEvo);
      server->Register("/CheckTiming", fT0RichDiffEvo);
      server->Register("/CheckTiming", fT0PsdDiffEvo);
      server->Register("/CheckTiming", fT0StsDiffEvoLong);
      server->Register("/CheckTiming", fT0MuchDiffEvoLong);
      server->Register("/CheckTiming", fT0TofDiffEvoLong);
      server->Register("/CheckTiming", fT0RichDiffEvoLong);
      server->Register("/CheckTiming", fT0PsdDiffEvoLong);

      server->Register("/CheckTiming", fT0StsMeanEvo);
      server->Register("/CheckTiming", fT0MuchMeanEvo);
      server->Register("/CheckTiming", fT0TofMeanEvo);
      server->Register("/CheckTiming", fT0RichMeanEvo);
      server->Register("/CheckTiming", fT0PsdMeanEvo);

      server->Register("/CheckTiming", fT0T0Diff);
      server->Register("/CheckTiming", fStsStsDiff);
      server->Register("/CheckTiming", fMuchMuchDiff);
      server->Register("/CheckTiming", fTofTofDiff);
      server->Register("/CheckTiming", fRichRichDiff);
      server->Register("/CheckTiming", fPsdPsdDiff);

      server->Register("/CheckTiming", fT0StsDpbDiff);
      for( UInt_t uStsDpb = 0; uStsDpb < kuMaxNbStsDpbs; ++uStsDpb )
        server->Register("/CheckTiming/STS", fT0StsDpbDiffEvo[uStsDpb]);

      server->Register("/CheckTiming", fT0MuchRocDiff);
      server->Register("/CheckTiming", fT0MuchAsicDiff);
      for( UInt_t uMuchAsic = 0; uMuchAsic < kuMaxNbMuchAsics; ++uMuchAsic )
        server->Register("/CheckTiming/MUCH", fT0MuchAsicDiffEvo[uMuchAsic]);
    }
  }
}
// ---- ReInit  -------------------------------------------------------
InitStatus CbmMcbm2019CheckPulser::ReInit()
{
  return kSUCCESS;
}

// ---- Exec ----------------------------------------------------------
void CbmMcbm2019CheckPulser::Exec(Option_t* /*option*/)
{
  LOG(debug) << "executing TS " << fNrTs;

  if( 0 < fNrTs && 0 == fNrTs % 1000 )
    LOG(info) << Form( "Processing TS %6d", fNrTs );


  CheckInterSystemOffset();
  fNrTs++;

  if( 0 < fNrTs && 0 == fNrTs % 90000 )
    WriteHistos();
}

void CbmMcbm2019CheckPulser::CheckInterSystemOffset()
{
  LOG(debug) <<"Begin";
  Int_t nrT0Digis=fT0Digis->GetEntriesFast();
  LOG(debug) << "T0Digis: " << nrT0Digis;

  Int_t nrStsDigis{0};
  if (fStsDigis) {
    nrStsDigis=fStsDigis->GetEntriesFast();
    LOG(debug) << "StsDigis: " << nrStsDigis;
  }

  Int_t nrMuchDigis{0};
  if (fMuchDigis) {
    nrMuchDigis=fMuchDigis->GetEntriesFast();
    LOG(debug) << "MuchDigis: " << nrMuchDigis;
  }

  Int_t nrTofDigis{0};
  if (fTofDigis) {
    nrTofDigis=fTofDigis->GetEntriesFast();
    LOG(debug) << "TofDigis: " << nrTofDigis;
  }

  Int_t nrRichDigis{0};
  if (fRichDigis) {
    nrRichDigis=fRichDigis->GetEntriesFast();
    LOG(debug) << "RichDigis: " << nrRichDigis;
  }

  Int_t nrPsdDigis{0};
  if (fPsdDigis) {
    nrPsdDigis=fPsdDigis->GetEntriesFast();
    LOG(debug) << "PsdDigis: " << nrPsdDigis;
  }

//  if (nrT0Digis < 100000) {
  if (nrT0Digis < 1000000) {
    /// Re-initialize array references
    fPrevT0FirstDigiSts = 0.;
    fPrevT0FirstDigiMuch = 0.;
    fPrevT0FirstDigiTof = 0.;
    fPrevT0FirstDigiRich = 0.;
    fPrevT0FirstDigiPsd = 0.;

    for (Int_t iT0 = 0; iT0 < nrT0Digis; ++iT0) {

      if (iT0%1000 == 0) LOG(debug) << "Executing entry " << iT0;

      CbmDigi* T0Digi = static_cast<CbmDigi*>(fT0Digis->At(iT0));

      /// Keep only pulser Digis
      if( fuMaxTotPulserT0 < T0Digi->GetCharge() || T0Digi->GetCharge() < fuMinTotPulserT0 )
        continue;

      Double_t T0Time = T0Digi->GetTime();
      Int_t T0Address = T0Digi->GetAddress();

      /// Keep only pulser Digis
      if( 0x00005006 != T0Address && 0x04005006 != T0Address )
        continue;

      fT0Address->Fill(T0Address);
/*
      std::cout << Form( "T0 pulser in TS %5d: address 0x%08X T0 time %12.0f dt %12.0f",
                         fNrTs, T0Address, T0Time, T0Time - fdLastT0DigiPulser )
                << std::endl;
*/
      fT0Channel->Fill(static_cast<CbmTofDigiExp*>(T0Digi)->GetChannel());

      if (nrStsDigis > 0 && nrStsDigis < 1000000 && fuMinAdcPulserSts < fuMaxAdcPulserSts )
        fPrevT0FirstDigiSts = FillSystemOffsetHistos( fStsDigis,
                                                      fT0StsDiff, fT0StsDiffEvo,
                                                      fT0StsDiffEvoLong, fT0StsMeanEvo, fT0StsDpbDiff,
                                                      T0Time, fStsOffsetRange, fPrevT0FirstDigiSts,
                                                      kSts);
      if (nrMuchDigis > 0 && nrMuchDigis < 1000000 && fuMinAdcPulserMuch < fuMaxAdcPulserMuch )
        fPrevT0FirstDigiMuch = FillSystemOffsetHistos(fMuchDigis,
                                                      fT0MuchDiff, fT0MuchDiffEvo,
                                                      fT0MuchDiffEvoLong, fT0MuchMeanEvo, fT0MuchRocDiff,
                                                      T0Time, fMuchOffsetRange, fPrevT0FirstDigiMuch,
                                                      kMuch);
      if (nrTofDigis > 0 && nrTofDigis < 1000000 && fuMinTotPulserTof < fuMaxTotPulserTof )
        fPrevT0FirstDigiTof = FillSystemOffsetHistos(fTofDigis,
                                                      fT0TofDiff, fT0TofDiffEvo,
                                                      fT0TofDiffEvoLong, fT0TofMeanEvo, nullptr,
                                                      T0Time, fTofOffsetRange, fPrevT0FirstDigiTof,
                                                      kTof);
      if (nrRichDigis > 0 && nrRichDigis < 1000000 && fuMinTotPulserRich < fuMaxTotPulserRich )
        fPrevT0FirstDigiRich = FillSystemOffsetHistos(fRichDigis,
                                                      fT0RichDiff, fT0RichDiffEvo,
                                                      fT0RichDiffEvoLong, fT0RichMeanEvo, nullptr,
                                                      T0Time, fRichOffsetRange, fPrevT0FirstDigiRich,
                                                      kRich);
      if (nrPsdDigis > 0 && nrPsdDigis < 1000000 && fuMinAdcPulserPsd < fuMaxAdcPulserPsd )
        fPrevT0FirstDigiPsd = FillSystemOffsetHistos(fPsdDigis,
                                                      fT0PsdDiff, fT0PsdDiffEvo,
                                                      fT0PsdDiffEvoLong, fT0PsdMeanEvo, nullptr,
                                                      T0Time, fPsdOffsetRange, fPrevT0FirstDigiPsd,
                                                      kPsd);
    }
  }
}

Int_t CbmMcbm2019CheckPulser::FillSystemOffsetHistos(TClonesArray* array,
                                             TH1* histo,
                                             TH2* histoEvo, TH2* histoEvoLong, TProfile * profMeanEvo,
                                             TH2* histoAFCK,
                                             const Double_t T0Time, const Int_t offsetRange,
                                             Int_t iStartDigi,
                                             ECbmModuleId iDetId )
{
  Int_t nrDigis=array->GetEntriesFast();
  Int_t iFirstDigiInWin = iStartDigi;

  if( kLastModule == iDetId )
  {
  } // if( kLastModule == iDetId )

  for (Int_t i = iStartDigi; i < nrDigis; ++i) {

    CbmDigi* Digi = static_cast<CbmDigi*>(array->At(i));

    switch( iDetId )
    {
      case kSts: ///< Silicon Tracking System
      {
        if( fuMaxAdcPulserSts < Digi->GetCharge() || Digi->GetCharge() < fuMinAdcPulserSts )
          continue;
        break;
      } // case kSts:
      case kMuch:        ///< Muon detection system
      {

        UInt_t uAsic = static_cast<CbmMuchBeamTimeDigi*>(Digi)->GetNxId();
        UInt_t uChan = static_cast<CbmMuchBeamTimeDigi*>(Digi)->GetNxCh();
        if( 9 != uAsic || uChan < 63 )
          continue;

/*
        UInt_t uChan = static_cast<CbmMuchBeamTimeDigi*>(Digi)->GetNxCh();
        if( 32 != uChan )
          continue;
*/
        if( fuMaxAdcPulserMuch < Digi->GetCharge() || Digi->GetCharge() < fuMinAdcPulserMuch )
          continue;
        break;
      } // case kMuch:
      case kTof:         ///< Time-of-flight Detector
      {
        if( fuMaxTotPulserTof < Digi->GetCharge() || Digi->GetCharge() < fuMinTotPulserTof )
          continue;
        break;
      } // case kTof:
      case kRich:        ///< Ring-Imaging Cherenkov Detector
      {
        if( fuMaxTotPulserRich < Digi->GetCharge() || Digi->GetCharge() < fuMinTotPulserRich )
          continue;
        break;
      } // case kRich:
      case kPsd:         ///< Projectile spectator detector
      {
        if( fuMaxAdcPulserPsd < Digi->GetCharge() || Digi->GetCharge() < fuMinAdcPulserPsd )
          continue;
        break;
      } // case kPsd:
      default:
        return 0;
    } // switch( iDetId )

    Double_t diffTime = T0Time - Digi->GetTime();
//    profMeanEvo->Fill( T0Time * 1e-9 - fdStartTime, diffTime );

    if (diffTime > offsetRange) {
      ++ iFirstDigiInWin; // Update Index of first digi in Win to next digi
      continue; // not yet in interesting range
    } // if (diffTime > offsetRange)
    if (diffTime < -offsetRange) break;     // already past interesting range
    histo->Fill(diffTime);
    histoEvo->Fill(fNrTs, diffTime);
    histoEvoLong->Fill(fNrTs, diffTime);

    if( -1 == fdStartTime )
      fdStartTime = T0Time * 1e-9;
    profMeanEvo->Fill( T0Time * 1e-9 - fdStartTime, diffTime );
/*
    if (kMuch == iDetId )
      std::cout << Form( "MUCH coinc in TS %5d: %7.2f T0 time %12.0f", fNrTs, diffTime, T0Time )
                << std::endl;
*/
    if( kPsd == iDetId )
      fT0PsdDiffCharge->Fill( diffTime, Digi->GetCharge() );

    /// STS DPB mapping: ladder 1 is in DPB 1
    if (kSts == iDetId && histoAFCK) {
      UInt_t uDPB = ( 0 < (Digi->GetAddress() & 0x00000400 ) );
      histoAFCK->Fill(uDPB, diffTime);
      if( uDPB < kuMaxNbStsDpbs )
        fT0StsDpbDiffEvo[ uDPB ]->Fill(fNrTs, diffTime);
    } // if (kSts == iDetId && histoAFCK)

    /// MUCH DPB mapping
    if (kMuch == iDetId && histoAFCK) {
      UInt_t afck = static_cast<CbmMuchBeamTimeDigi*>(Digi)->GetRocId();
      UInt_t asic = static_cast<CbmMuchBeamTimeDigi*>(Digi)->GetNxId();
      histoAFCK->Fill(afck, diffTime);
      fT0MuchAsicDiff->Fill(asic, diffTime);
      if( asic < kuMaxNbMuchAsics )
        fT0MuchAsicDiffEvo[ asic ]->Fill(fNrTs, diffTime);
    } // if (kMuch == iDetId && histoAFCK)
  }

  return iFirstDigiInWin;
}

// ---- Finish --------------------------------------------------------
void CbmMcbm2019CheckPulser::Finish()
{
  WriteHistos();

}

void CbmMcbm2019CheckPulser::WriteHistos()
{
  TFile* old = gFile;
  TFile* outfile = TFile::Open(fOutFileName,"RECREATE");

  fT0StsDiff->Write();
  fT0MuchDiff->Write();
  fT0TofDiff->Write();
  fT0RichDiff->Write();
  fT0PsdDiff->Write();
  fT0PsdDiffCharge->Write();

  fT0StsDiffEvo->Write();
  fT0MuchDiffEvo->Write();
  fT0TofDiffEvo->Write();
  fT0RichDiffEvo->Write();
  fT0PsdDiffEvo->Write();

  fT0StsDiffEvoLong->Write();
  fT0MuchDiffEvoLong->Write();
  fT0TofDiffEvoLong->Write();
  fT0RichDiffEvoLong->Write();
  fT0PsdDiffEvoLong->Write();

  fT0StsMeanEvo->Write();
  fT0MuchMeanEvo->Write();
  fT0TofMeanEvo->Write();
  fT0RichMeanEvo->Write();
  fT0PsdMeanEvo->Write();

  fT0T0Diff->Write();
  fStsStsDiff->Write();
  fMuchMuchDiff->Write();
  fTofTofDiff->Write();
  fRichRichDiff->Write();
  fPsdPsdDiff->Write();

  fT0Address->Write();
  fT0Channel->Write();

  fT0StsDpbDiff->Write();
/*
  for( UInt_t uStsDpb = 0; uStsDpb < kuMaxNbStsDpbs; ++uStsDpb )
  {
    fT0StsDpbDiffEvo[uStsDpb]->Write();
    fStsDpbCntsEvo[uStsDpb]->Write();
  }
*/
  fT0MuchRocDiff->Write();
  fT0MuchAsicDiff->Write();
/*
  for( UInt_t uMuchAsic = 0; uMuchAsic < kuMaxNbMuchAsics; ++uMuchAsic )
    fT0MuchAsicDiffEvo[uMuchAsic]->Write();
*/
  outfile->Close();
  delete outfile;

  gFile = old;
}

ClassImp(CbmMcbm2019CheckPulser)