// Gaudi
#include "GaudiKernel/PhysicalConstants.h"
#include "GaudiUtils/HistoLabels.h"
#include "GaudiUtils/Aida2ROOT.h"

// Tb/TbEvent
#include "Event/TbTrack.h"
#include "Event/TbCluster.h"

// Tb/TbKernel
#include "TbKernel/TbConstants.h"
#include "TbKernel/TbModule.h"

// Local
#include "TbTimewalkMonitor.h"

// ROOT
#include "TH1D.h"

using namespace Gaudi::Utils::Histos;

DECLARE_ALGORITHM_FACTORY(TbTimewalkMonitor)

  //=============================================================================
  // Standard constructor
  //=============================================================================
TbTimewalkMonitor::TbTimewalkMonitor(const std::string& name, 
    ISvcLocator* pSvcLocator)
  : TbAlgorithm(name, pSvcLocator) {

    declareProperty("TrackLocation",
        m_trackLocation = LHCb::TbTrackLocation::Default);
    declareProperty("WidthMin", m_widthMin = 37 );
    declareProperty("WidthMax", m_widthMax = 42 );
  }

//=============================================================================
// Initialization
//=============================================================================
StatusCode TbTimewalkMonitor::initialize() {

  // Initialise the base class.
  StatusCode sc = TbAlgorithm::initialize();
  if (sc.isFailure()) return sc;

  // Book the histograms.
  for (unsigned int i = 0; i < m_nPlanes; ++i) {
    const std::string plane = std::to_string(i);
    const std::string title = geomSvc()->modules().at(i)->id();
    m_space.push_back(bookProfile1D("SpatialResidual/Plane"+plane,title,-0.5,255.5,256));

    m_spd2D.push_back(book2D("spd2D/Plane"+plane,title,-0.5,255.5,256,-100,100,20));
    m_twd.push_back(book2D("twd/Plane"+plane,title,1.,101.,100,0.,62.5,40));
    m_twdQ.push_back(book2D("twdQ/Plane"+plane,title,1000,10000,100,0.,62.5,40));
    m_twdQL.push_back(book2D("twdQL/Plane"+plane,title,1000,10000,100,0.,62.5,40));
    m_twdQR.push_back(book2D("twdQR/Plane"+plane,title,1000,10000,100,0.,62.5,40));

    setAxisLabels( m_twd[i], "ToT", "#Delta t [ns]");

    m_twdL.push_back(book2D("twdL/Plane"+plane,title,0.,100.,100,0.,62.5,40));
    m_twdR.push_back(book2D("twdR/Plane"+plane,title,0.,100.,100,0.,62.5,40));
    m_LRSYNC.emplace_back(bookProfile1D("LR_SYNC/Plane"+plane,title,-0.5,255.5,256));
    m_UDSYNC.emplace_back(bookProfile1D("UD_SYNC/Plane"+plane,title,-0.5,255.5,256));

    m_quad.emplace_back(bookProfile1D("QUAD/Plane"+plane,title,-0.5,255.5,256));

    m_inscol.emplace_back(bookProfile1D("Scol/Plane"+plane,title,-0.5,127.5,128));
    m_interscol.emplace_back(bookProfile1D("InterScol/Plane"+plane,title,-0.5,127.5,128));

    m_timewalk.push_back(bookProfile1D("Timewalk/Plane"+plane,title, 1., 401., 400));
    m_dtDist.push_back( book1D( "dtDist/Plane"+plane,title, 0., 90., 58 ) );
    m_cDist.push_back( book1D( "cDist/Plane"+plane,title,0.,90.,58 ) );
    m_timewalkQ.push_back(bookProfile1D("TimewalkQ/Plane"+plane,title,0,20000,200));
    setAxisLabels(m_timewalkQ[i], "charge [e^{-}]","#Delta t [ns]");
    setAxisLabels(m_timewalk[i], "ToT", "#Delta t [ns]");
  }

  return StatusCode::SUCCESS;
}

StatusCode TbTimewalkMonitor::finalize() {

  for( unsigned int i = 0 ; i < m_nPlanes; ++i ){
    TH1D* tdist = Gaudi::Utils::Aida2ROOT::aida2root( m_dtDist[i] );
    AIDA::IHistogram1D* cdist = m_cDist[i];
    int total = 0;    
    for (int i=0;i<tdist->GetNbinsX();i++)
    {
      total += tdist->GetBinContent(i);
      cdist->fill( tdist->GetBinCenter(i) , total /tdist->GetEntries() ) ;
    }    
  }
  return StatusCode::SUCCESS; 
}

//=============================================================================
// Main execution
//=============================================================================
StatusCode TbTimewalkMonitor::execute() {

  // Grab the tracks.

  const LHCb::TbTracks* tracks = getIfExists<LHCb::TbTracks>(m_trackLocation);
  if (!tracks) return StatusCode::SUCCESS; 
  for (const LHCb::TbTrack* track : *tracks) {
 //   info() << track << endmsg; 
    const SmartRefVector<LHCb::TbCluster>& clusters = track->clusters();
    const auto ref = std::min_element(
        clusters.begin(), clusters.end(),
        [](const LHCb::TbCluster* h1,
          const LHCb::TbCluster* h2) { return h1->htime() < h2->htime(); });
    const double tref = (*ref)->htime();
    const double syncRef = (*clusters.rbegin())->htime();
    unsigned int p = (*clusters.rbegin())->plane();

    for (auto it = clusters.cbegin(), end = clusters.cend(); it != end; ++it){

      auto hits = (*it)->hits();
      for (const auto hit : hits) {
        if( p==5 ){
          m_space[(*it)->plane()]->fill( hit->col(), hit->htime() - syncRef );
          m_spd2D[(*it)->plane()]->fill( hit->col(), hit->htime() - syncRef );
        }
        m_timewalk[(*it)->plane()]->fill(hit->ToT(), hit->htime() - tref);
        m_twd[(*it)->plane()]->fill(hit->ToT(), hit->htime() - tref);
        m_timewalkQ[(*it)->plane()]->fill(hit->charge(), hit->htime() - tref );
        m_twdQ[(*it)->plane()]->fill(hit->charge(), hit->htime() - tref );
        m_dtDist[(*it)->plane()]->fill( hit->htime() - tref );   
      }
      if( (*it)->size() == 2 ){
        auto minmax = std::minmax_element(hits.begin(), hits.end(),[](LHCb::TbHit* h1, LHCb::TbHit* h2) { return h1->col() < h2->col(); });
        if( (*minmax.first)->col()  != (*minmax.second)->col() ){
          m_twdL[(*it)->plane()]->fill( (*minmax.first)->ToT(), (*minmax.first)->htime() - tref );
          m_twdR[(*it)->plane()]->fill( (*minmax.second)->ToT(), (*minmax.second)->htime() - tref );
          m_twdQL[(*it)->plane()]->fill( (*minmax.first)->charge(), (*minmax.first)->htime() - tref );
          m_twdQR[(*it)->plane()]->fill( (*minmax.second)->charge(), (*minmax.second)->htime() - tref );
        }
      }
    }
    const SmartRefVector<LHCb::TbCluster>& associatedClusters  = track->associatedClusters();

    for (auto it = associatedClusters.cbegin(), end = associatedClusters.cend(); it != end; ++it){
      auto hits = (*it)->hits();
      //info() << "Filling associated clusters for " << track << endmsg; 
      for (const auto hit : (*it)->hits()) {
        m_timewalk[(*it)->plane()]->fill(hit->ToT(), hit->htime() - tref);
        m_twd[(*it)->plane()]->fill(hit->ToT(), hit->htime() - tref);
        //m_twdQ[(*it)->plane()]->fill(hit->charge(), hit->htime() - tref );
        //m_dtDist[(*it)->plane()]->fill( hit->htime() - tref );
        m_timewalkQ[(*it)->plane()]->fill(hit->charge(), hit->htime() - tref );
        //if( (*it)->size() == 1 ) m_timewalkOneHit[(*it)->plane()]->fill(hit->ToT(), ( hit->htime() - tref ));
        //if( (*it)->size() == 2 ) m_timewalkTwoHit[(*it)->plane()]->fill(hit->ToT(), ( hit->htime() - tref ));
      }
      /*
      if( (*it)->size() == 2 ){
        auto minmax = std::minmax_element(hits.begin(), hits.end(),[](LHCb::TbHit* h1, LHCb::TbHit* h2) { return h1->col() < h2->col(); });
        if( (*minmax.first)->col()  != (*minmax.second)->col() ){
          m_twdL[(*it)->plane()]->fill( (*minmax.first)->ToT(), (*minmax.first)->htime() - tref );
          m_twdR[(*it)->plane()]->fill( (*minmax.second)->ToT(), (*minmax.second)->htime() - tref );
          m_twdQL[(*it)->plane()]->fill( (*minmax.first)->charge(), (*minmax.first)->htime() - tref );
          m_twdQR[(*it)->plane()]->fill( (*minmax.second)->charge(), (*minmax.second)->htime() - tref );
        }
      }
      */
      //info() << "Filled associated clusters " << endmsg; 
    }
  }

  return StatusCode::SUCCESS;
}
/*
   for(unsigned int plane=0;plane<m_nPlanes; ++plane ){
   const LHCb::TbClusters* clusters = 
   getIfExists<LHCb::TbClusters>(LHCb::TbClusterLocation::Default+std::to_string(plane));
   for( auto& cluster : *clusters ){
   if( cluster->size() == 2 ){
   auto hits = cluster->hits();
   auto minmax = std::minmax_element(hits.begin(), hits.end(),
   [](LHCb::TbHit* h1, LHCb::TbHit* h2) { return h1->col() < h2->col(); });
   const LHCb::TbHit* left = (*minmax.first);
   const LHCb::TbHit* right = (*minmax.second);
   if( left->col() != right->col() ){
   m_LRSYNC[plane]->fill( left->col() , left->htime() - right->htime() );
   if( int(left->col() /2 ) == int(right->col() /2 ) ) 
   m_inscol[plane]->fill( int(left->col()/2) , left->htime() - right->htime() );
   else 
   m_interscol[plane]->fill( int(left->col()/2) , left->htime() - right->htime() );
   }
   else {
   minmax = std::minmax_element(hits.begin(), hits.end(),
   [](LHCb::TbHit* h1, LHCb::TbHit* h2) { return h1->row() < h2->row(); });
   m_UDSYNC[plane]->fill( (*minmax.first)->row(), 
   (*minmax.first)->htime() - (*minmax.second)->htime());
   }
   }
   if( cluster->size() == 4 ){
   auto hits = cluster->hits();
   auto minmax = std::minmax_element(hits.begin(), hits.end(),
   [](LHCb::TbHit* h1, LHCb::TbHit* h2) { return h1->col() < h2->col(); });
   const LHCb::TbHit* left = (*minmax.first);
   const LHCb::TbHit* right = (*minmax.second);
   if( left->col() == right->col() + 4 ){
   m_quad[plane]->fill( left->col() , left->htime() - right->htime() );
   }
   }

   }
   }
   return StatusCode::SUCCESS;
   }
   */