// -------------------------------------------------------------------------
// -----                PndRichHitFinder source file                   -----
// -----                created 06 Apr 2017                            -----
// -------------------------------------------------------------------------


#include "TClonesArray.h"
#include "TArrayD.h"
#include "TGeoManager.h"

#include "FairRootManager.h"
#include "PndRichHitFinder.h"
#include "FairRun.h"
#include "FairRuntimeDb.h"
#include "FairGeoNode.h"
#include "FairGeoVector.h"
#include "FairRunAna.h"
#include "FairEventHeader.h"

//#include "PndRichDigi.h"

#include "PndDetectorList.h"

// -----   Default constructor   -------------------------------------------
PndRichHitFinder::PndRichHitFinder() :
FairTask("RichHitFinder", 1)
{
  fPersistency = kTRUE;
  fPixelHits = 0;
  fEventNr = 0;
  fPixelFactor = 1;
  fHitNumber = 0;
  
  if(fVerbose>0) Info("PndRichHitFinder","RICH Hit Finder created, Parameters will be taken from RTDB");  
  fGeoH=NULL;
  fGeo = new PndRichGeo();
  fDigiArray   = NULL;
  fPdHitArray   = NULL;
  fMCEventHeader = NULL;  
  fStopFunctor = NULL;

  fInd = 0;
  fTimeGate = 1;
  fTimeStep = 0.5; //ns
  fIBuffer = 0;
  fIBufferPrev = 1;
  fBufferStartTime.resize(2);
  fBufferStartTime.at(0) = 0;
  fBufferStartTime.at(1) = 0;
  fHitsBuffer.resize(2);
  fHitsBuffer.at(0).resize(1000);
  fHitsBuffer.at(1).resize(1000);
  fBufferNumHits.resize(2);
  fBufferNumHits.at(0) = 0;
  fBufferNumHits.at(1) = 0;
  fThreshold = 5;
  petime = -1;

//  fPixelSize = fGeo->PixelSize();	//pixel size;    
//  fNpix = fGeo->Npixels(); //pixel columns
//  fMcpActiveArea = fGeo->McpActiveArea();
//  fPixelGap      =  (fMcpActiveArea - (Double_t)fNpix*fPixelSize) / ((Double_t)fNpix - 1.);
//  fPixelStep     =  fMcpActiveArea/fNpix;//fPixelSize + 0.5*fPixelGap;

}
// -------------------------------------------------------------------------

PndRichHitFinder::PndRichHitFinder(Int_t iVerbose) :
  FairTask("RichHitFinder", iVerbose)
{
  fPersistency = kTRUE;
  fPixelHits = 0;
  fEventNr = 0; 
  fPixelFactor = 1;
  fGeoH = NULL;
  fGeo = new PndRichGeo();
  fDigiPixelMCInfo = kFALSE;  
  if(fVerbose>0) Info("PndRichHitFinder","RichHitFinder created, Parameters will be taken from RTDB");

  fDigiArray   = NULL;
  fPdHitArray   = NULL;
  fMCEventHeader = NULL;  
  fStopFunctor = NULL;
   
  fInd = 0;
  fTimeGate = 1;
  fTimeStep = 0.5; //ns
  fIBuffer = 0;
  fIBufferPrev = 1;
  fBufferStartTime.resize(2);
  fBufferStartTime.at(0) = 0;
  fBufferStartTime.at(1) = 0;
  fHitsBuffer.resize(2);
  fHitsBuffer.at(0).resize(1000);
  fHitsBuffer.at(1).resize(1000);
  fBufferNumHits.resize(2);
  fBufferNumHits.at(0) = 0;
  fBufferNumHits.at(1) = 0;
  fThreshold = 5;
  petime = -1;

//  fPixelSize = fGeo->PixelSize();	//pixel size;  
//  fNpix = fGeo->Npixels();              //pixel rows in one FE
//  fMcpActiveArea = fGeo->McpActiveArea();
//  fPixelGap      =  (fMcpActiveArea - (Double_t)fNpix*fPixelSize) / ((Double_t)fNpix - 1.);
//  fPixelStep     =  fPixelSize + 0.5*fPixelGap;
}

PndRichHitFinder::PndRichHitFinder(const char* name, Int_t iVerbose) :
FairTask(name, iVerbose)
{
  fPersistency = kTRUE;
  fPixelHits = 0;
  fEventNr = 0;   
  fPixelFactor = 1;
  fGeoH = NULL;
  fGeo = new PndRichGeo();
  fDigiPixelMCInfo = kFALSE;  
  if(fVerbose>0) Info("PndRichHitFinder","%s created, Parameters will be taken from RTDB",name);

  fDigiArray   = NULL;
  fPdHitArray   = NULL;
  fMCEventHeader = NULL;

  fTimeStep = 0.5; //ns
  fIBuffer = 0;
  fHitsBuffer.resize(2);
  fHitsBuffer.at(0).resize(1000);
  fHitsBuffer.at(1).resize(1000);
  petime = -1;

//  fPixelSize = fGeo->PixelSize();	//pixel size;  
//  fNpix = fGeo->Npixels();              //pixel rows in one FE
//  fMcpActiveArea = fGeo->McpActiveArea();
//  fPixelGap      =  (fMcpActiveArea - (Double_t)fNpix*fPixelSize) / ((Double_t)fNpix - 1.);
//  fPixelStep     =  fPixelSize + 0.5*fPixelGap;
}


// -----   Destructor   ----------------------------------------------------
PndRichHitFinder::~PndRichHitFinder(){
 if (fGeo) delete fGeo;
 if (fGeoH) delete fGeoH;	
}
// -------------------------------------------------------------------------

// -----   Initialization  of Parameter Containers -------------------------
void PndRichHitFinder::SetParContainers(){
  if ( fGeoH == NULL ) fGeoH = PndGeoHandling::Instance();
  fGeoH->SetParContainers(); 
  fGeoH->SetVerbose(fVerbose);
  if(fVerbose>1) Info("SetParContainers","done.");
  return;
}

InitStatus PndRichHitFinder::ReInit(){
  SetParContainers();
  return kSUCCESS;
}

// -----   Public method Init   --------------------------------------------
InitStatus PndRichHitFinder::Init(){
  //FairRun* ana = FairRun::Instance(); //[R.K. 01/2017] unused variable?
  FairRootManager* ioman = FairRootManager::Instance();
  if ( ! ioman )
  {
    std::cout << "-E- PndRichHitFinder::Init: "
    << "RootManager not instantiated!" << std::endl;
    return kFATAL;
  }

  if (FairRunAna::Instance()->IsTimeStamp())
      fInBranchName = "RichDigiSorted";
  else
      fInBranchName = "RichDigi";
  // Get input array
  fDigiArray = (TClonesArray*) ioman->GetObject(fInBranchName);  
  if ( ! fDigiArray )
  {
    std::cout << "-W- PndRichHitFinder::Init: "
    << "No RichDigi array!" << std::endl;
    return kERROR;
  }
  
  // Create and register output array
  fPdHitArray	= new TClonesArray("PndRichPDHit");
  ioman->Register("RichPDHit", "Rich", fPdHitArray, fPersistency);
 
  fGapFunctor = new TimeGap();
  fStopFunctor = new StopTime();
  return kSUCCESS;
}

// -----   Public method Exec   --------------------------------------------
void PndRichHitFinder::Exec(Option_t*){
  
  if(fVerbose>3) Info("Exec","Start");  
  if (!fPdHitArray) Fatal("Exec", "No PdHitArray");
  fPdHitArray->Clear();

  Double_t etime = FairRootManager::Instance()->GetEventTime();
  if (FairRunAna::Instance()->IsTimeStamp()){
    //fDigiArray = FairRootManager::Instance()->GetData(fInBranchName, fStopFunctor, etime, fStopFunctor, etime+100);
    fDigiArray = FairRootManager::Instance()->GetData(fInBranchName, fGapFunctor, 10);
  }

  Int_t nDigis = fDigiArray->GetEntriesFast();
  Double_t  timef = ((PndRichDigi*)fDigiArray->At(0))->GetTime();
  Double_t  timel = ((PndRichDigi*)fDigiArray->At(nDigis-1))->GetTime();
  fBufferNumHits.at(fIBuffer) = (timel-fBufferStartTime.at(fIBuffer))/fTimeStep+1;
  Double_t stopTime = fBufferStartTime.at(fIBuffer) + fBufferNumHits.at(fIBuffer)*fTimeStep;
  if ((Int_t)fHitsBuffer.at(fIBuffer).size()<fBufferNumHits.at(fIBuffer))
      fHitsBuffer.at(fIBuffer).resize(fBufferNumHits.at(fIBuffer));
  for(Int_t i=0; i<fBufferNumHits.at(fIBuffer);i++)fHitsBuffer.at(fIBuffer).at(i)=0;
  for (Int_t iDigi = 0; iDigi < nDigis; iDigi++)
   {
      Double_t  hitTime = ((PndRichDigi*)fDigiArray->At(iDigi))->GetTime();
      Int_t ind = (hitTime-fBufferStartTime.at(fIBuffer))/fTimeStep;
      if (hitTime>fBufferStartTime.at(fIBuffer))
      {
         fHitsBuffer.at(fIBuffer).at(ind)++;
      }
      else
      {
         ind += fBufferNumHits.at(fIBufferPrev)-1;
         fHitsBuffer.at(fIBufferPrev).at(ind)++;
      }
   }
  Int_t nbr = fBufferNumHits.at(fIBuffer);
  Int_t nbrp = fBufferNumHits.at(fIBufferPrev);
  Int_t im = nbr + nbrp - 3*fTimeGate;
  Int_t imax_;
  Int_t nhitsmax = 0;
   std::vector< Int_t > nmax;
   std::vector< Int_t > imax;
   std::vector< Int_t > istart;
   std::vector< Int_t > istop;
   Int_t nhitsprev = 0;
   Int_t ilast;
//  for(Int_t i = fInd; i<im; i++)
//   {
//      std::cout << i << " " << (i<nbrp) << " "
//         << ( i<nbrp ?
//              fHitsBuffer.at(fIBufferPrev).at(i) :
//              fHitsBuffer.at(fIBuffer).at(i-nbrp) ) << std::endl;
//   }
  for(Int_t i = fInd; i<im; i++)
   {
      Int_t nhits = 0;
      for(Int_t j = i-fTimeGate; j<=i+fTimeGate; j++)
         if (j>=0) nhits += j<nbrp ?
                            fHitsBuffer.at(fIBufferPrev).at(j) :
                            fHitsBuffer.at(fIBuffer).at(j-nbrp);
//      std::cout << i << " " << (i<nbrp) << " "
//         << ( i<nbrp ?
//              fHitsBuffer.at(fIBufferPrev).at(i) :
 //             fHitsBuffer.at(fIBuffer).at(i-nbrp) ) << " " << nhits << " " << nhitsmax;
      if ((nhits>fThreshold)&&(nhitsprev<=fThreshold)) istart.push_back(i);
      if ((nhits>fThreshold)&&(nhits>nhitsmax))
      {
         nhitsmax = nhits;
         imax_ = i;
      }
//      std::cout << std::endl;
      if ((nhits<=fThreshold)&&(nhitsprev>fThreshold))
      {
         nmax.push_back(nhitsmax);
         imax.push_back(imax_);
         istop.push_back(i-1);
         Double_t t = fBufferStartTime.at(fIBufferPrev) + imax_*fTimeStep;
         std::cout << "peak = " << imax.size() << " " << imax_ << " " << nhitsmax << " "
            << istop.back() << " " << istart.back() << " " << t << std::endl;
         nhitsmax = 0;
         imax_ = -1;
      }
      if (nhits<5) ilast = i;
      nhitsprev = nhits;
   }
  fInd = ilast -  nbrp;
   
  if(fVerbose>1) std::cout<<"-I- PndRichHitFinder: Event # "<< fEventNr<<" has "<<nDigis<<" digis."<< std::endl;
  else if(fVerbose==1 && fEventNr%1000==0) std::cout<<"-I- PndRichHitFinder: Event # "<< fEventNr<<" has "<<nDigis<<" digis."<< std::endl;
   std::cout<<"-I- PndRichHitFinder: Event # "<< fEventNr<<" has "<<nDigis<<" digis.   "
      << FairRunAna::Instance()->IsTimeStamp() << " " << etime << " "
      << timel-timef << " " << etime << " " << timef << " " << timel << " " << petime << std::endl;

   std::cout << "buffer: " << fIBuffer << " " << fBufferNumHits.at(fIBuffer) << " " <<
      fBufferStartTime.at(fIBuffer) << std::endl;
 
  Int_t detID = 0;//, mcpID = 0, pixelID = 0;
  TVector3 HitPosGlobal, HitPosLocal, dPosHit;
  Double_t hitTime = 0.;
  fHitNumber += nDigis;
    
  for (Int_t iDigi = 0; iDigi < nDigis; iDigi++){
    fDigi = (PndRichDigi*) fDigiArray->At(iDigi); 
    detID = -999;//fDigi->GetDetectorId();    
//    pixelID = detID - 100*(Int_t)TMath::Floor((Double_t)detID/100.);
//    mcpID = detID/100;
    hitTime = fDigi->GetTime();
     Int_t sensorId = -999;
        
    // the pixel number shows local coordinates of the hit:
//    HitPosLocal.SetXYZ(fPixelStep*((Double_t)(pixelID % fNpix) - (Double_t)(fNpix/2) + 0.5),
//		       fPixelStep*(TMath::Floor(((Double_t)pixelID)/((Double_t)fNpix)) - (Double_t)(fNpix/2) + 0.5), 0.);
//    Int_t sensorId = fDigi->GetSensorId()/fPixelFactor;
//    if(fPixelFactor==2) { //double pixels
//      sensorId++;
//      if(fDigi->GetSensorId()%2) HitPosLocal.SetX(HitPosLocal.X()+fPixelSize/2.);
//      else  HitPosLocal.SetX(HitPosLocal.X()-fPixelSize/2.);
//    }

    // local coordinates of the hit on the MCP are translated into global ones as the following:
//    HitPosGlobal = fGeoH->LocalToMasterShortId(HitPosLocal, mcpID);   
//    dPosHit.SetXYZ(fPixelSize/2., fPixelSize/2., 0.);   
//    if(fPixelFactor==2)  dPosHit.SetXYZ(fPixelSize, fPixelSize/2., 0.);  
    

    if(fDigi->GetTimeStamp()!=etime+hitTime) hitTime =  fDigi->GetTimeStamp() - etime;
    PndRichPDHit pdhit = PndRichPDHit(iDigi, detID, sensorId , HitPosGlobal, dPosHit, hitTime, 0.);
//    pdhit.SetPdgCode(fDigi->GetPdgCode());
    pdhit.SetTimeStamp(fDigi->GetTimeStamp());
//    pdhit.SetTimeAtBar(fDigi->GetTimeAtBar());
    pdhit.SetLink(fDigi->GetLink(0)); // MCTrack
    pdhit.AddLink(fDigi->GetLink(1)); // RichPDPoint
    pdhit.AddLink(FairLink(-1,fEventNr, "RichDigi", iDigi));
    //((FairMultiLinkedData)pdhit).Print(); std::cout<<std::endl;
    new((*fPdHitArray)[fPdHitArray->GetEntriesFast()]) PndRichPDHit(pdhit); 
  }
  
  fEventNr++;
  if(fVerbose>3) Info("Exec","Loop MC points");

  fPdHitArray->Sort();
  fDigiArray->Delete();

  // switch between buffers
  fIBufferPrev = fIBuffer;
  fIBuffer = fIBuffer ? 0 : 1;
  fBufferStartTime.at(fIBuffer) = stopTime;
  petime = etime;
}


// -------------------------------------------------------------------------
void PndRichHitFinder::FinishEvent()
{
  FinishEvents();
}

// -------------------------------------------------------------------------
void PndRichHitFinder::FinishTask()
{
   std::cout << "-I- PndRichHitFinder: Total number of hits is " << fHitNumber << std::endl;
}

ClassImp(PndRichHitFinder);