//////////////////////////////////////////////////////////////////////////////
//
// File: hrichanalysissim.cc
//
// $Id: hrichanalysis.cc,v 1.37 2009-07-15 11:39:17 halo Exp $
//
//*-- Author : W. Przygoda
//*-- Modified : Tue Feb 15 18:50:45 CET 2005 martin.jurkovic@ph.tum.de
//*-- Modified : Tue Feb 12 16:54:10 CET 2002 teberl@ph.tum.de
//*--            introduced ctr flag to skip events
//*-- Modified Oct. 2000 by W. Koenig (W.Koenig@gsi.de)
//*-- Modified : 13/05/2001 by Laura Fabbietti
//
//_HADES_CLASS_DESCRIPTION 
//////////////////////////////////////////////////////////////////////////////
//
//  HRichAnalysis
//  
//  
//////////////////////////////////////////////////////////////////////////////


#include "hades.h"
#include "hcategory.h"
#include "hdebug.h"
#include "hevent.h"
#include "heventheader.h"
#include "hiterator.h"
#include "hlinearcategory.h"
#include "hparset.h"
#include "hrichanalysis.h"
#include "hrichanalysispar.h"
#include "hrichcal.h"
#include "hrichdetector.h"
#include "hrichgeometrypar.h"
#include "hrichhitheader.h"
#include "hrichpadclean.h"
#include "hrichpadlabel.h"
#include "hrichringfind.h"
#include "hruntimedb.h"
#include "hspectrometer.h"
#include "richdef.h"

#include <iomanip>
#include <iostream> 

using namespace std;

ClassImp(HRichAnalysis)

//----------------------------------------------------------------------------
HRichAnalysis::HRichAnalysis() : HReconstructor()
{

   clear();
 
   pPadClean = new HRichPadClean;
   pPadLabel = new HRichPadLabel;
   pRingFind = new HRichRingFind;

}
//============================================================================

//----------------------------------------------------------------------------
HRichAnalysis::HRichAnalysis(const Text_t* name,const Text_t* title, Bool_t kSkip)
                                  : HReconstructor(name, title)
{
   clear();

   pPadClean = new HRichPadClean;
   pPadLabel = new HRichPadLabel;
   pRingFind = new HRichRingFind;
   kSkipEvtIfNoRing = kSkip;
}
//============================================================================

//----------------------------------------------------------------------------
HRichAnalysis::~HRichAnalysis()
{
   if ( NULL != pPadClean )
   {
      delete pPadClean;
      pPadClean = NULL;
   }
   if ( NULL != pPadLabel )
   {
      delete pPadLabel;
      pPadLabel = NULL;
   }
   if ( NULL != pRingFind )
   {
      delete pRingFind;
      pRingFind = NULL;
   }
   if ( NULL != pLabelArea )
   {
      delete [] pLabelArea;
      pLabelArea = NULL;
   }
   if ( NULL != pPads )
   {
      for (Int_t i = 0; i < 6; i++)
         if (pPads[i]) delete [] pPads[i]; 
      delete [] pPads;
      pPads = NULL;
   }
   if( NULL != pLeftBorder ){
      delete [] pLeftBorder;
      pLeftBorder = NULL;
   }
   if( NULL != pRightBorder )
   {
      delete [] pRightBorder;
      pRightBorder = NULL;
   }
   if ( NULL != fIter )
   {
      delete fIter;
      fIter = NULL;
   }
   if ( NULL != fIterHitHeader )
   {
      delete fIterHitHeader;
      fIterHitHeader = NULL;
   }
}
//============================================================================

//----------------------------------------------------------------------------
void HRichAnalysis::clear()
{

   pLabelArea     = NULL;
   pRings         = NULL;
   pPads          = NULL;
   pSectorPads    = NULL;
   fIter          = NULL;
   fIterHitHeader = NULL;
   pLeftBorder    = NULL;
   pRightBorder   = NULL;
   m_pCalCat      = NULL;
   m_pHitCat      = NULL;
   m_pHitHdrCat   = NULL;
   fpAnalysisPar  = NULL;
   fpGeomPar      = NULL;
   hithdrLoop     = NULL;

   for ( Int_t i = 0; i < 6; ++i )
      fPadFired[i] = 0;

   iClustersCleaned.Reset();

   allPairNrTot      = 0;
   iActiveSector     = 0;
   iClusterCleanedNr = 0;
   iFakeLocalMax4    = 0;
   iFakeLocalMax8    = 0;
   iFakePad          = 0;
   iLabelNr          = 0;
   iPadCleanedNr     = 0;
   iPadFiredNr       = 0;
   iRingNr           = 0;
   iRingNrTot        = 0;
   sectorPairNrTot   = 0;

}
//============================================================================

//----------------------------------------------------------------------------
void
HRichAnalysis::Reset()
{
   iClustersCleaned.Reset();

   iClusterCleanedNr = 0;
   iFakeLocalMax4    = 0;
   iFakeLocalMax8    = 0;
   iFakePad          = 0;
   iLabelNr          = 0;
   iPadCleanedNr     = 0;
   iPadFiredNr       = 0;
   iRingNr           = 0;

   if ( pLabelArea )
   {
      delete [] pLabelArea;
      pLabelArea = NULL; 
   }
}
//============================================================================

//----------------------------------------------------------------------------
HRichAnalysis::HRichAnalysis(const HRichAnalysis& source)
{

   Error("HRichAnalysis", "HRichAnalysis object can not be initialized with values of another object!\nUse instead default constructor");

}
//============================================================================

//----------------------------------------------------------------------------
HRichAnalysis&
HRichAnalysis::operator=(const HRichAnalysis& source)
{
   if (this != &source)
   {
      Error("operator=", "HRichAnalysis object can not be assigned!");
   }
   return *this;
}
//============================================================================

//----------------------------------------------------------------------------
Bool_t
HRichAnalysis::init()
{
// Allocate input/output categories, initialize iterators and parameter containers

   HRichDetector*    pRichDet     = NULL;
   HRuntimeDb*       rtdb         = NULL;
   HRichAnalysisPar* pAnalysisPar = NULL;
   HRichGeometryPar* pGeomPar     = NULL;

   pRichDet = static_cast<HRichDetector*>(gHades->getSetup()
                                          ->getDetector("Rich"));
   if ( NULL == pRichDet )
   {
      Error("init", "Pointer to HRichDetector is NULL");
      return kFALSE;
   }

   m_pCalCat = gHades->getCurrentEvent()->getCategory(catRichCal);
   if ( NULL == m_pCalCat )
   {
      m_pCalCat = pRichDet->buildCategory(catRichCal);
      if ( NULL == m_pCalCat )
      {
         Error("init", "Pointer to HRichCal category is NULL");
         return kFALSE;
      }
      else
         gHades->getCurrentEvent()
         ->addCategory(catRichCal, m_pCalCat, "Rich");
   }
   fIter = static_cast<HIterator*>(m_pCalCat->MakeIterator());
   if ( NULL == fIter )
   {
      Error("init", "Pointer to HRichCalSim iterator is NULL");
      return kFALSE;
   }

   m_pHitCat = gHades->getCurrentEvent()->getCategory(catRichHit);
   if ( NULL == m_pHitCat )
   {
      m_pHitCat = pRichDet->buildCategory(catRichHit);
      if ( NULL == m_pHitCat )
      {
         Error("init", "Pointer to HRichHit category is NULL");
         return kFALSE;
      }
      else gHades->getCurrentEvent()
      ->addCategory(catRichHit, m_pHitCat, "Rich");
   }

   m_pHitHdrCat = gHades->getCurrentEvent()->getCategory(catRichHitHdr);
   if ( NULL == m_pHitHdrCat )
   {
      m_pHitHdrCat = pRichDet->buildCategory(catRichHitHdr);
      if ( NULL == m_pHitHdrCat )
      {
         Error("init", "Pointer to HRichHitHdr category is NULL");
         return kFALSE;
      }
      else gHades->getCurrentEvent()
      ->addCategory(catRichHitHdr, m_pHitHdrCat, "Rich");
   }
   fIterHitHeader = static_cast<HIterator*>(m_pHitHdrCat->MakeIterator("native"));
   if ( NULL == fIterHitHeader )
   {
      Error("init", "Pointer to HRichHitHdr iterator is NULL");
      return kFALSE;
   }


   rtdb = gHades->getRuntimeDb();
   pAnalysisPar = static_cast<HRichAnalysisPar*>(rtdb->
                                                 getContainer("RichAnalysisParameters"));
   if ( NULL == pAnalysisPar )
   {
      Error("init", "Pointer to HRichAnalysisPar parameters is NULL");
      return kFALSE;
   }
   setAnalysisPar(pAnalysisPar);

   pGeomPar = static_cast<HRichGeometryPar*>(rtdb->
                                            getContainer("RichGeometryParameters"));
   if ( NULL == pGeomPar )
   {
      Error("init", "Pointer to RichGeometryParameters parameters is NULL");
      return kFALSE;
   }
   setGeometryPar(pGeomPar);

   if ( kFALSE == initParameters() )
   {
      Error("init", "Parameters not initialized");
      return kFALSE;
   }

   if ( kFALSE == pPadClean->init() )
   {
      Error("init", "HRichPadClean did not initialize");
      return kFALSE;
   }

   return kTRUE;
}
//============================================================================

//----------------------------------------------------------------------------
Bool_t
HRichAnalysis::reinit()
{

   Int_t i,j,k,m,n;
   Int_t iMatrixSurface = 0;
   Int_t iPartSurface   = 0;
   Int_t iMaskSize      = 0;
   Int_t iMaskSizeSquared = 0;

   if ( kFALSE == initParameters() )
   {
      Error("reinit", "Parameters not initialized");
      return kFALSE;
   }

   pRingFind->init(this);


   iPadActive.Set(maxCols * maxRows);
   for ( i = 0; i < maxCols * maxRows; ++i )
      if ( ((HRichGeometryPar*)fpGeomPar)->getPadsPar()->getPad(i)->getPadActive() > 0 )
         iPadActive[i] = 1; 
      else
         iPadActive[i] = 0;


   iMaskSize        = ((HRichAnalysisPar*)fpAnalysisPar)->iRingMaskSize;
   iMaskSizeSquared = iMaskSize*iMaskSize;
   for ( k = 0; k < iMaskSizeSquared; ++k )
      if ( 1 == ((HRichAnalysisPar*)fpAnalysisPar)->iRingMask[k] )
         iMatrixSurface++;


   for ( j = 0; j < maxRows; j++ )
      for ( i = 0; i < maxCols; i++ )
      {
         iPartSurface = 0;
         for ( k = 0; k < iMaskSizeSquared; k++ )
         {
            m = (k % iMaskSize) - iMaskSize/2;
            n = (k / iMaskSize) - iMaskSize/2;
            if ( !IsOut(i,j,m,n) )
               if ( 1 == ((HRichAnalysisPar*)fpAnalysisPar)->iRingMask[k] )
                  iPartSurface++;
         }

         ((HRichGeometryPar*)fpGeomPar)->getPadsPar()->
         getPad(i,j)->setAmplitFraction((Float_t)iPartSurface / (Float_t)iMatrixSurface);
      }

   return kTRUE;
}
//============================================================================

//----------------------------------------------------------------------------
Bool_t
HRichAnalysis::initParameters()
{
// allocate non event by event classes

   HRichGeometryPar *pGeomPar = getGeometryPar();                                  

   iRingNrTot=0;
   allPairNrTot=0;
   sectorPairNrTot=0;
   maxFiredTotalPads  = 3000; // upper limit of number of fired pads

   maxCols = pGeomPar->getColumns();
   maxRows = pGeomPar->getRows();

   if ( pLeftBorder )
      delete [] pLeftBorder;
   pLeftBorder = new Short_t[maxRows];

   if ( pRightBorder )
      delete [] pRightBorder;
   pRightBorder = new Short_t[maxRows];

   for ( Int_t row = 0; row < maxRows; ++row )
   {
      Int_t col = 0;
      Int_t padOffset = row * maxCols;
      while ( !pGeomPar->getPadsPar()->getPad(col+padOffset)->getPadActive() && col<maxCols )
         ++col;
      if( col == maxCols )
      {
         maxRows = row;
         break;
      }
      pLeftBorder[row] = col;
      while ( pGeomPar->getPadsPar()->getPad(col + padOffset)->getPadActive() && col < maxCols )
         ++col;
      pRightBorder[row] = col - 1;
   }
   maxPads = maxRows * maxCols;

   // now creating pads array
   pPads = new HRichPadSignal * [6];
   for ( Int_t j = 0; j < 6; ++j )
   {
      pSectorPads = pPads[j] = new HRichPadSignal[maxPads];
      for ( Int_t i = 0; i < maxPads; pSectorPads[i++].clear() );
   }
   pSectorPads = pPads[iActiveSector];

   return kTRUE;
}
//============================================================================

//----------------------------------------------------------------------------
Bool_t
HRichAnalysis::finalize()
{
   return kTRUE;
}
//============================================================================

//----------------------------------------------------------------------------
Int_t
HRichAnalysis::execute()
{

   HRichCal* pCal      = NULL;
   Bool_t skipSector[] = {0, 0, 0, 0, 0, 0};
   Int_t iRingNrTotEvt = 0;
   Int_t fPadFiredTot  = 0;
   Int_t lastRingNr    = 0;
   Int_t ampl          = -1;
   Int_t sector        = -1;
   Int_t padnr         = -1;
   const Int_t eventNr = gHades->getCurrentEvent()->getHeader()->getEventSeqNumber();
   Int_t j;

   secWithCurrent = -1;
   Reset();

   for ( sector = 0; sector < 6; ++sector )
   {
      if ( fPadFired[sector] > 0 )
      {
         fPadFired[sector] = 0;
         HRichPadSignal* pSecPads = pPads[sector];
         for ( j = 0; j < maxPads; pSecPads[j++].clear() );
      }
   }

   fIter->Reset();
   while ( NULL != (pCal = (HRichCal *)fIter->Next()) )
   {
      if ( (ampl = (Int_t)pCal->getCharge()) > 0)
      {
         sector = pCal->getSector();
         fPadFired[sector]++;
         padnr = pCal->getCol() + pCal->getRow() * maxCols;
         pPads[sector][padnr].setAmplitude(ampl);
      }
      // remember that columns are x and rows are y
   }

   for( sector = 6; sector--; fPadFiredTot += fPadFired[sector] )
   {
      if ( fPadFiredTot > maxFiredTotalPads )
      {
         Warning("execute", "Analysis of event %d skipped: too many pads fired %i/%i",
                 eventNr, fPadFiredTot, maxFiredTotalPads);
         fIter->Reset();
         while( NULL != (pCal = static_cast<HRichCal*>(fIter->Next())) )
         {
            pCal->setIsCleanedSector( kTRUE );
         }
         return 0;
      }
      if ( fPadFired[sector] >= fpAnalysisPar->maxFiredSectorPads )
      {
         Warning("execute", "Event %d: too many fired pads in sector %i. %i/%i",
                 eventNr, sector, fPadFired[sector], fpAnalysisPar->maxFiredSectorPads);
         fIter->Reset();
         while( NULL != (pCal = static_cast<HRichCal*>(fIter->Next())) )
         {
            if ( pCal->getSector() == sector )
               pCal->setIsCleanedSector( kTRUE );
         }
         skipSector[sector] = kTRUE;
      }
   }

  
   // **************************************************************

   // ------- loop over sectors --- begin ---

   lastRingNr = iRingNrTot;
   for ( sector = 0; sector < 6; ++sector )
   {
      // Skip analysis for sector with too many fired pads
      if ( kTRUE == skipSector[sector] )
         continue;

      if ( ((HRichGeometryPar*)fpGeomPar)->getSectorActive(sector) > 0 )
      {
         Reset();
         SetActiveSector(sector);
         iPadFiredNr = fPadFired[sector];
  
         // for each sector:
         //   - cleaning
         //   - labeling
         //   - ring finding
         // procedure are executed

         iPadCleanedNr = pPadClean->Execute(this);
         iLabelNr      = pPadLabel->Execute(this);

         if ( 0 == ((HRichAnalysisPar*)fpAnalysisPar)->isActiveLabelPads )
         {
            iLabelNr = 1;
            pLabelArea = new HRichLabel[1];
            pLabelArea[0].iLeftX = 0;
            pLabelArea[0].iRightX = maxCols;
            pLabelArea[0].iLowerY = 0;
            pLabelArea[0].iUpperY = maxRows;
            pLabelArea[0].iLabeledPadsNr = maxPads;
            pLabelArea[0].iFiredPadsNr = iPadFiredNr;
         }

         iRingNr = pRingFind->Execute(this);
         iRingNrTot    += iRingNr;
         iRingNrTotEvt += iRingNr;

         if (iRingNr > 1)
         {
            ++sectorPairNrTot;
         }

         // Found rings are stored in the hit container 
         updateHeaders(sector, eventNr);
         updateHits(sector);
      }
   }

   // update information about cleaned pads
   fIter->Reset();
   while( NULL != (pCal = (HRichCal *)fIter->Next()) )
   {
      sector = pCal->getSector();
      padnr  = pCal->getCol() + pCal->getRow() * maxCols;
      pCal->setIsCleanedSingle(pPads[sector][padnr].getIsCleanedSingle());
      pCal->setIsCleanedHigh  (pPads[sector][padnr].getIsCleanedHigh());
   }

   // ------- loop over sectors --- end ---
   if( iRingNrTot > lastRingNr + 1)
   {
      ++allPairNrTot;
   }

   // modification to skip event which does not contain any ring
   if ( kSkipEvtIfNoRing && iRingNrTotEvt == 0 )
      return kSkipEvent;
   return 0;
}
//============================================================================

//----------------------------------------------------------------------------
void
HRichAnalysis::updateHits(Int_t nSec)
{

   HRichHit* hit = NULL;
   HLocation loc;

   for (Int_t i = 0; i < iRingNr; i++)
   {
      loc.set(1, nSec);

      hit = (HRichHit*)m_pHitCat->getNewSlot(loc);

      if ( NULL != hit )
      {
         hit = new(hit) HRichHit;

         *hit = pRings[i];
         hit->setEventNr(gHades->getCurrentEvent()->getHeader()->getEventSeqNumber());
         hit->setSector(nSec);
         if ( hit->fX > 1000. )
            Info("updateHits", "fX : %i, %i fY : %i, %i",
                 hit->fX, hit->iRingX, hit->fY, hit->iRingY);

         HRichPad* pad = ((HRichGeometryPar*)fpGeomPar)->getPadsPar()->
         getPad(pRings[i].iRingX,pRings[i].iRingY);

         // now fX and fY is something else - mean cluster x,y
//        pad->getXY(&hit->fX, &hit->fY);
//        pad->getXYZlab(nSec, &hit->fLabX, &hit->fLabY, &hit->fLabZ);
         hit->fPhi   = pad->getPhi(nSec);
         hit->fTheta = pad->getTheta();

      }
   }
}      
//============================================================================

//----------------------------------------------------------------------------
void
HRichAnalysis::updateHeaders(const Int_t nSec,
                             const Int_t nEvNr )
{
   HRichHitHeader* hithdr = NULL;
   HLocation       loc;

   fIterHitHeader->Reset();
   hithdrLoop = NULL;
 
   while ( NULL != (hithdrLoop = (HRichHitHeader *)fIterHitHeader->Next()) )
   {
      secWithCurrent = hithdrLoop->getSector();
      if ( nSec == secWithCurrent )
      {
         hithdr = hithdrLoop;
         break;
      }
   }

   if ( hithdr == NULL )
   {
      loc.set(1, nSec);
      hithdr = (HRichHitHeader*)(m_pHitHdrCat)->getNewSlot(loc);
      if ( NULL != hithdr )
         hithdr = new (hithdr) HRichHitHeader;
   }

   if ( NULL != hithdr )
   {
      hithdr->setSector(nSec);
      hithdr->setEventNr(nEvNr);
      hithdr->setExecutedAnalysis(1);
      hithdr->iPadFiredNr = iPadFiredNr;
      hithdr->iPadCleanedNr = iPadCleanedNr;
      hithdr->iClusterCleanedNr = iClusterCleanedNr;
      hithdr->iClustersCleaned = iClustersCleaned;
      hithdr->iRingNr = iRingNr;
      hithdr->iFakePad = iFakePad;
      hithdr->iFakeLocalMax4 = iFakeLocalMax4;
      hithdr->iFakeLocalMax8 = iFakeLocalMax8;
   }
   else
   {
      Warning("updateHeaders", "No Header Object could be created");
   }
}

//============================================================================
  
//----------------------------------------------------------------------------
void
HRichAnalysis::SetActiveSector(Int_t sectornr)
{
   if ( sectornr == iActiveSector )
      return;
   if ( sectornr > 5 || sectornr < 0 )
   {
      Error("SetActiveSector", "Sector number %d out of range (0..5)!", sectornr);
      return;
   }
   else if ( 0 == ((HRichGeometryPar*)fpGeomPar)->getSectorActive(sectornr) )
   {
      Error("SetActiveSector", "Sector number %d is not present (and cannot be active)!", sectornr);
      return;
   }
   else
   {
      iActiveSector = sectornr;
      pSectorPads = pPads[sectornr];
   } 
}
//============================================================================

//----------------------------------------------------------------------------
Int_t
HRichAnalysis::SetNextSector()
{
   Int_t i = iActiveSector;
   while ( i < 6 )
   {
      i++;
      if ( ((HRichGeometryPar*)fpGeomPar)->getSectorActive(i) > 0 )
      {
         pSectorPads = pPads[i];
         return (iActiveSector = i);
      }
   }
   Error("SetNextSector", "No more sectors (last sector $d)", iActiveSector);
   return iActiveSector;
}
//============================================================================

//----------------------------------------------------------------------------
Bool_t
HRichAnalysis::IsOut(Int_t nowPad, Int_t dx, Int_t dy)
{
   if ( nowPad <= 0 )
      return 1;
   Int_t x = nowPad % maxCols;
   Int_t y = nowPad / maxCols;

   if ( (x + dx) >= 0 && (x + dx) < maxCols &&
        (y + dy) >= 0 && (y + dy) < maxRows &&
        iPadActive[nowPad+dx + maxCols*dy]     )
   {
      return 0;
   }
   else
   {
      return 1;
   }
}
//============================================================================

//----------------------------------------------------------------------------
Bool_t
HRichAnalysis::IsOut(Int_t x, Int_t y, Int_t dx, Int_t dy)
{
   if ( (x + dx) >= 0 && (x + dx) < maxCols &&
        (y + dy) >= 0 && (y + dy) < maxRows &&
        iPadActive[x + dx + maxCols * (y + dy)] )
   {
      return 0;
   }
   else
   {
      return 1;
   }
}
//============================================================================

//----------------------------------------------------------------------------
void
HRichAnalysis::Streamer(TBuffer &R__b)
{
// Stream an object of class HRichAnalysis.

   if (R__b.IsReading())
   {
      Version_t R__v = R__b.ReadVersion(); if (R__v) { }
      R__b >> iPadFiredNr;
      R__b >> iPadCleanedNr;
      R__b >> iRingNr;
      R__b >> pRings;
      R__b >> iClusterCleanedNr;
      iClustersCleaned.Streamer(R__b);
      R__b >> iFakePad;
      R__b >> iFakeLocalMax4;
      R__b >> iFakeLocalMax8;
   }
   else
   {
      R__b.WriteVersion(HRichAnalysis::IsA());
      R__b << iPadFiredNr;
      R__b << iPadCleanedNr;
      R__b << iRingNr;
      R__b << pRings;
      R__b << iClusterCleanedNr;
      iClustersCleaned.Streamer(R__b);
      R__b << iFakePad;
      R__b << iFakeLocalMax4;
      R__b << iFakeLocalMax8;
   }
}
//============================================================================