//////////////////////////////////////////////////////////////////////////////
//
// File: hrichpadclean.cc
//
// $Id: hrichpadclean.cc,v 1.16 2009-04-09 16:34:56 jurkovic Exp $
//
//*-- Author : Witold Przygoda (przygoda@psja1.if.uj.edu.pl)
//*-- Modified : 1999/12/04 by Witold Przygoda (przygoda@psja1.if.uj.edu.pl)
//*-- Modified : Oct. 2000 W. Koenig
//*-- Modified : Mar 12 18:50:45 CET 2005 martin.jurkovic@ph.tum
//
//_HADES_CLASS_DESCRIPTION 
//////////////////////////////////////////////////////////////////////////////
//
//  HRichPadClean
//
//  This class cleans pad plane provided by HRichAnalysis. This class 
//  isn't any task but it is called for each event from
//  the HRichAnalysis::execute() function.
//
//////////////////////////////////////////////////////////////////////////////


#include "hades.h"
#include "hevent.h"
#include "hmatrixcategory.h"
#include "hrichanalysispar.h"
#include "hrichdetector.h"
#include "hrichgeometrypar.h"
#include "hrichdirclus.h"
#include "hrichpadclean.h"
#include "hspectrometer.h"
#include "richdef.h"

#include <iomanip>
#include <iostream> 
#include <stdlib.h>

using namespace std;

ClassImp(HRichPadClean)

//----------------------------------------------------------------------------
HRichPadClean::HRichPadClean()
{
   iCount       = 0;
   minX         = 200;
   maxX         = 200;
   minY         = 0;
   maxY         = 0;
   deltaPhi     = 0;
   fRichClusCat = NULL;
   iTempCluster.Set(100);
}
//============================================================================

//----------------------------------------------------------------------------
HRichPadClean::~HRichPadClean()
{

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

//----------------------------------------------------------------------------
HRichPadClean::HRichPadClean(const HRichPadClean& source)
{
   iCount = source.iCount;
   iTempCluster = source.iTempCluster;
}
//============================================================================

//----------------------------------------------------------------------------
HRichPadClean&
HRichPadClean::operator=(const HRichPadClean& source)
{
   if (this != &source)
   {
      iCount = source.iCount;
      iTempCluster = source.iTempCluster;
   }
   return *this;
}
//============================================================================

//----------------------------------------------------------------------------
Bool_t
HRichPadClean::init()
{
   fRichClusCat = gHades->getCurrentEvent()->getCategory(catRichDirClus);
   if ( NULL == fRichClusCat )
   {
      fRichClusCat = gHades->getSetup()->getDetector("Rich")->buildCategory(catRichDirClus);
      if ( NULL == fRichClusCat )
      {
         Error("init", "Category catRichDirClus could not be created");
         return kFALSE;
      }
      else
         gHades->getCurrentEvent()->addCategory(catRichDirClus, fRichClusCat, "Rich");
   }
   return kTRUE;
}
//============================================================================

//----------------------------------------------------------------------------
Int_t
HRichPadClean::Execute(HRichAnalysis *giveMe)
{
// Cleaning algorithm works in two steps:
//   1. First the charged particle clusters are identified and cleaned. 
//      Cleaning is triggered by any pad, which charge is higher than
//      the value set in HRichAnalysisPar::iCleanHighPulseUpperThreshold.
//      If such a pad is found, all "continuously" surrounding area is 
//      cleaned as well. The limit for "continuous" is defined by the 
//      parameter HRichAnalysisPar::iCleanHighPulseBorder
//   2. Afterwards all stand alone pads with charge below certain threshold
//      defined in HRichAnalysisPar::iCleanAlonePadLowerThreshold are cleaned.
//      Condition for stand alone pad is defined by
//      HRichAnalysisPar::iCleanAlonePadBorder.

   if ( giveMe->iPadFiredNr > 2 )
   {
      // what is the reason for this check???? MJ
      Int_t iReducedNr = 0;
      HRichAnalysisPar *analParams = giveMe->getAnalysisPar(); 

      if (analParams->isActiveCleanHighPulse) 
         iReducedNr += CleanHighPulse(giveMe,analParams->iCleanHighPulseBorder,
                                      analParams->iCleanHighPulseUpperThreshold);

      if (analParams->isActiveCleanAlonePad) 
         iReducedNr += CleanAlonePad(giveMe,analParams->iCleanAlonePadBorder,
                                     analParams->iCleanAlonePadLowerThreshold);

      return (iReducedNr);
   } else
      return giveMe->iPadFiredNr; //either 0, 1 or 2
}

//============================================================================

//----------------------------------------------------------------------------
Int_t 
HRichPadClean::CleanHighPulse(HRichAnalysis* showMe, 
                              Int_t          border, 
                              Int_t          upperThr) {
//
// This method deletes all the pads above upper threshold (upperThr)
// and the cluster connected to each pad.
//

   if ( 0 == border )
   {
      Error("CleanAlonePad", "Argument \'border\' in function CleanAlonePad must be larger than 0!");
      return -1;
   }

#ifdef HRICH_DEBUGMODE
   cout << "RICH DEBUG MODE: DeletePulse cleans the following pads: \n";
#endif

   const Int_t maxRows = showMe->GetPadsYNr();
   const Int_t maxCols = showMe->GetPadsXNr();
   Int_t iCountBefore  = 0;
   Int_t nCleaned      = 0;
   Int_t offset        = 0;
   Int_t leftBorder    = 0;
   Int_t rightBorder   = 0;
   HRichDirClus* clus  = NULL;
   HLocation loc;

   Int_t old_row = 0;
   Int_t old_col = 0;

   iCount = 0;

   for( Int_t row = 0; row < maxRows; ++row )
   {
      offset      = row * maxCols;
      leftBorder  = showMe->pLeftBorder[row];
      rightBorder = showMe->pRightBorder[row];
      for( Int_t col = showMe->pLeftBorder[row]; col <= rightBorder; ++col )
      {
         if ( showMe->GetPad(col + offset)->getAmplitude() > upperThr )
         {
            chargeTot = 0; // fix MJ
            minX = 200;
            maxX = 0;
            minY = 200; 
            maxY = 0;
            deltaPhi = 0;
            iCountBefore = iCount;
            old_row = row;
            old_col = col;
            DeletePulse(showMe, border, col, row);
            if ( old_row != row || old_col != col )
               Error("CleanHighPulse", "Numbers changed: old = [%d,%d], new = old = [%d,%d]", old_row, old_col, row, col);
            if ( (iCount - iCountBefore) > 0 && nCleaned < 100 )
            {
               iTempCluster[nCleaned++] = iCount - iCountBefore;

               // Store the cluster and its parameters to the HRichDirClus category
               clus = NULL;
               loc.set(1, 0);
               clus = static_cast<HRichDirClus*>( static_cast<HMatrixCategory*>(fRichClusCat)->getNewSlot(loc) );
   
               if( NULL == clus ) 
                  Error("CleanHighPulse","Error no slot free in HRichDirClus !!!!!!!!!!!  ");
               else 
               {
                  clus = new(clus) HRichDirClus;
                  clus->setSector(showMe->GetActiveSector());
                  clus->setX((Int_t)(maxX + minX)/2);
                  clus->setY((Int_t)(maxY + minY)/2);
                  clus->setXYDim(TMath::Abs(minX - maxX),
                                 TMath::Abs(minY - maxY));
                  clus->setTotalCharge(chargeTot);
                  clus->setNrPad(iCount - iCountBefore);
                  clus->setPhiDiff(calculateDPhi(showMe, minX, minY, maxX, maxY));
                  Float_t theta = -1;
                  if ( NULL != ((HRichGeometryPar*)showMe->getGeometryPar())->getPadsPar()->
                                 getPad( ((Int_t)(maxX + minX)/2), ((Int_t)(maxY + minY)/2)) )
                  {
                     theta = ((HRichGeometryPar*)showMe->getGeometryPar())->getPadsPar()->
                             getPad( ((Int_t)(maxX + minX)/2), ((Int_t)(maxY + minY)/2))->getTheta();
                  }
                  clus->setTheta( theta );
               }
            }
         }
      }
   }

   showMe->iClusterCleanedNr = nCleaned;
   showMe->iClustersCleaned.Set(nCleaned);

   for (Int_t i = 0; i < nCleaned; i++)
      showMe->iClustersCleaned[i] = iTempCluster[i];

#ifdef HRICH_DEBUGMODE
   if (!iCount) cout << "None. \n" << "RICH DEBUG MODE: Total number "
                     << "of pads cleaned by CleanHighPulse = 0\n";
   else cout << "\nRICH DEBUG MODE: Total number of pads cleaned "
             << "by CleanHighPulse = " << iCount << "\n";
#endif

   return iCount;
}
//============================================================================

//----------------------------------------------------------------------------
Int_t
HRichPadClean::CleanAlonePad(HRichAnalysis *showMe,
                             Int_t          border,
                             Int_t          lowerThr)
{
// This method deletes all single pad under threshold (lowerThr).
// A single pad is defined as pad whose first fired neighbour
// has at least a distance of "border" from the pad


   const Int_t     maxCols      = showMe->GetPadsXNr();
   const Int_t     maxRows      = showMe->GetPadsYNr();
   Int_t           leftBorder;
   Int_t           rightBorder;
   Int_t           offset;
   HRichPadSignal* pad          = NULL;

   iCount = 0;

   if ( 0 == border )
   {
      Error("CleanAlonePad", "Argument \'border\' in function CleanAlonePad must be larger than 0!");
      return -1;
   }

#ifdef HRICH_DEBUGMODE
   cout << "RICH DEBUG MODE: CleanAlonePad cleans the following pads: \n";
#endif
 

   for( Int_t row = 0; row < maxRows; ++row )
   {
      leftBorder  = showMe->pLeftBorder[row];
      rightBorder = showMe->pRightBorder[row];
      pad = showMe->GetPad(leftBorder + row * maxCols);
      for( Int_t col = leftBorder; col <= rightBorder; ++col )
      {
         if (pad->getAmplitude() > 0)
         {
            for (Int_t k = row - border; k <= row + border; ++k)
            {
               offset = k * maxCols;
               for ( Int_t l = col - border; l <= col + border; ++l )
                  if ( ! (l == col && k == row))
                     if ( !showMe->IsOut(l, k) )
                        if (showMe->GetPad(l + offset)->getAmplitude() > 0)
                           goto nextPad;
            }

            if (pad->getAmplitude() < lowerThr)
            {
               pad->setAmplitude(0);
               pad->setIsCleanedSingle(kTRUE);
               iCount++;

#ifdef HRICH_DEBUGMODE
               cout << " (" << colStart+border << "," << rowStart+border << ")... ";
#endif

            }
         }
        nextPad:
         ++pad;
      }
   }

#ifdef HRICH_DEBUGMODE
   if ( 0 == iCount) cout << "None. \n";
   else cout << "\nRICH DEBUG MODE: Total number of pads cleaned by "
   "CleanAlonePad = " << iCount << "\n";
#endif

   return iCount;
}
//============================================================================

//----------------------------------------------------------------------------
void
HRichPadClean::DeletePulse(HRichAnalysis* showYou,
                           Int_t          border, 
                           Int_t          col, 
                           Int_t          row)
{
// This method deletes a given pad and all the neighbours within
// a distance of "border"

   const Int_t maxCols = showYou->GetPadsXNr();
   Int_t       offset;

   chargeTot += showYou->GetPad(col + maxCols * row)->getAmplitude();

#ifdef HRICH_DEBUGMODE
   cout << " charge  " << showYou->GetPad(col + maxCols * row)->getAmplitude()
        << "  chargeTot  " <<chargeTot << endl;
#endif

   showYou->GetPad(col + maxCols * row)->setAmplitude(0);
   showYou->GetPad(col + maxCols * row)->setIsCleanedHigh(kTRUE);


#ifdef HRICH_DEBUGMODE
   cout << ((HRichGeometryPar*)showYou->getGeometryPar())->getPadsPar()->getPad(col, row)
        << "  first pad  with the second method" << endl;
   cout << "  ==================== it is assignment time!========= " << endl;
   cout << " minX " <<   minX    << " maxY  " <<   maxY    << "  col  " << col << " minY  "
        <<   minY   << " row   " <<   row     << " maxY  " << maxY << endl;
#endif

   minX = minX >= col ? col:minX;
   minY = minY >= row ? row:minY;
   maxX = maxX <= col ? col:maxX;
   maxY = maxY <= row ? row:maxY;

#ifdef HRICH_DEBUGMODE
   cout << "  the harder you try the dumber you look " << endl;
   cout << " minX " << minX << " maxX " << maxX << "  col  " << col
        << " minY " << minY << " row  " << row  << " maxY  " << maxY<<endl;
#endif

   iCount++;

#ifdef HRICH_DEBUGMODE
   cout << " (" << colStart+border << "," << rowStart+border << ")... ";
#endif

   for (Int_t i = row - border; i <= row + border; ++i) 
   {
      offset = maxCols * i;
      for ( Int_t j = col - border; j <= col + border; ++j )
         if ( !showYou->IsOut(j, i) )
            if( !(i == row && j == col) )
               if ( showYou->GetPad(j + offset)->getAmplitude() > 0 )
                  DeletePulse(showYou, border, j, i);
   }

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

//----------------------------------------------------------------------------
Float_t
HRichPadClean::calculateDPhi(HRichAnalysis *showMe,
                             Int_t          xmin,
                             Int_t          ymin,
                             Int_t          xmax,
                             Int_t          ymax  )
{
   Float_t phi1 = -1;
   Float_t phi2 = -1;

#ifdef HRICH_DEBUGMODE
   cout <<" xmin " <<xmin<<" ymin   "<< ymin<<" xmax  "<<xmax <<"   ymax "<<ymax<<endl;
#endif
  
   Int_t maxCols = showMe->GetPadsXNr();

#ifdef HRICH_DEBUGMODE
   cout <<((HRichGeometryPar*)showMe->getGeometryPar())->getPadsPar()->getPad(xmin,ymin)<<"  first pad  "<<endl;
   cout <<((HRichGeometryPar*)showMe->getGeometryPar())->getPadsPar()->getPad(xmin+maxCols*ymin)<<"  first pad  "<<endl;
   cout <<((HRichGeometryPar*)showMe->getGeometryPar())->getPadsPar()->getPad(xmax+maxCols*ymax)<<"  second pad  "<<endl;
   cout << showMe->GetPad(xmin+maxCols*ymin)<<" first pad, second method  "<<endl; 
#endif
   
   if( showMe->GetPad(xmin+maxCols*ymin) )
      phi1 = ((HRichGeometryPar*)showMe->getGeometryPar())->getPadsPar()->getPad(xmin+maxCols*ymin)->getPhi(showMe->GetActiveSector());
   if ( showMe->GetPad(xmax+maxCols*ymax) )
      phi2 = ((HRichGeometryPar*)showMe->getGeometryPar())->getPadsPar()->getPad(xmax+maxCols*ymax)->getPhi(showMe->GetActiveSector());


#ifdef HRICH_DEBUGMODE
   cout<<"  delta Phi "<<TMath::Abs(phi1-phi2)<<endl;
#endif


   if ( -1 == phi1 && -1 == phi2 )
      return -1;
  
   return TMath::Abs(phi1-phi2);
}
//============================================================================