//////////////////////////////////////////////////////////////////////////////
//
// File: hrichanalysissim.cc
//
// $Id: hrichanalysissim.cc,v 1.27 2009-07-15 11:39:17 halo Exp $
//
//*-- Author : Laura Fabbietti
//*-- Modified : 08/05/2000 L.Fabbietti
//*-- Modified : Oct. 2000 W. Koenig (reduced content of HRichPadSignal)
//*-- Modified : Tue Feb 12 16:54:10 CET 2002 teberl@ph.tum.de
//*--            introduced ctr flag to skip events
//
//_HADES_CLASS_DESCRIPTION 
//////////////////////////////////////////////////////////////////////////////
//
//  HRichAnalysisSim
//
//  This class besides the Rich Analysis contains
//  the necessary functions to retrieve the parent track 
//  number corresponding to each pad that belongs to
//  a Ring structure.
//  The Ring Finding algorithm is executed by the class
//  HRichRingFindSim. This class contains 2 different
//  algorithms to find a ring, once a ring is found
//  the track numbers of the particles that have hit
//  the pads belonging to the ring structure are memorized
//  in an array. (Int_t *iRingTrack). 
//  Once all track numbers are collected, 
//  3 tracks number are saved ( HRichHitSim:track1....)
//  for each Hit with the corresponding weights (cfr. updateHits).
//  (The weight
//  is the number of pads a photon or a Ip has fired.)
//  A flag for each of the three track numbers is saved too,
//  to be able to distinguish between a Chrenkov photon
//  and a direct hit.  
//
//  use kSkip in the constructor to skip writing of 
//  event if no ring was found
//
//////////////////////////////////////////////////////////////////////////////


#include "hades.h"
#include "hevent.h"
#include "heventheader.h"
#include "hiterator.h"
#include "hlinearcategory.h"
#include "hlocation.h"
#include "hmatrixcategory.h"
#include "hparcond.h"
#include "hparset.h"
#include "hrichanalysispar.h"
#include "hrichanalysissim.h"
#include "hrichcalsim.h"
#include "hrichdetector.h"
#include "hrichgeometrypar.h"
#include "hrichhitsim.h"
#include "hrichpadclean.h"
#include "hrichpadlabel.h"
#include "hrichringfindsim.h"
#include "hrichtrack.h"
#include "hruntimedb.h"
#include "hspectrometer.h"
#include "richdef.h"

#include <iomanip>
#include <iostream> 

using namespace std;

ClassImp(HRichAnalysisSim)

//----------------------------------------------------------------------------
HRichAnalysisSim::HRichAnalysisSim(void)
{
   pRingFindSim   = new HRichRingFindSim;
   pRichCalSim    = NULL;
   fIterHitHeader = NULL;
   m_pTrackCat    = NULL;
   pRings         = NULL;

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

//----------------------------------------------------------------------------
HRichAnalysisSim::HRichAnalysisSim(const Text_t*  name,
                                   const Text_t*  title,
                                   Bool_t   kSkip  )
                                  : HRichAnalysis(name, title, kSkip)
{
   pRingFindSim     = new HRichRingFindSim;
   pRichCalSim      = NULL;
   fIterHitHeader   = NULL;
   m_pTrackCat      = NULL;
   pRings           = NULL;
   kSkipEvtIfNoRing = kSkip;
}
//============================================================================

//----------------------------------------------------------------------------
HRichAnalysisSim::~HRichAnalysisSim(void)
{
   if ( NULL != pRingFindSim )
   {
      delete pRingFindSim;
      pRingFindSim = NULL;
   }
   if( NULL != fIterHitHeader )
   {
      delete  fIterHitHeader;
      fIterHitHeader = NULL;
   }
}
//============================================================================

//----------------------------------------------------------------------------
Bool_t
HRichAnalysisSim::init()
{
// allocate input/output categories

   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->buildMatrixCat("HRichCalSim", 1);
      if ( NULL == m_pCalCat )
      {
         Error("init", "Pointer to HRichCalSim 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->buildLinearCat("HRichHitSim");
      if ( NULL == m_pHitCat )
      {
         Error("init", "Pointer to HRichHitSim 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;
   }

   m_pTrackCat = gHades->getCurrentEvent()->getCategory(catRichTrack);
   if ( NULL == m_pTrackCat )
   {
      m_pTrackCat = pRichDet->buildLinearCat("HRichTrack", 1002);
      if ( NULL == m_pTrackCat )
      {
         Error("init", "Pointer to HRichTrack category is NULL");
         return kFALSE;
      }
      else gHades->getCurrentEvent()
      ->addCategory(catRichTrack,m_pTrackCat , "Rich");
   }
     
   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 HRichAnalysisSim::reinit()
{

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

    pRingFindSim->init(this);


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

  Int_t i,j,k,m,n;

    Int_t iMatrixSurface=0, iPartSurface=0;
    Int_t iMaskSize = getAnalysisPar()->iRingMaskSize;
    Int_t iMaskSizeSquared = iMaskSize*iMaskSize;
    for (k = 0; k < iMaskSizeSquared; k++)
       if (getAnalysisPar()->iRingMask[k] == 1) 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 (getAnalysisPar()->iRingMask[k] == 1) iPartSurface++;
        }
        getGeometryPar()->getPadsPar()->
          getPad(i,j)->setAmplitFraction((Float_t)iPartSurface/(Float_t)iMatrixSurface);
      }


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

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

    iRingNrTot=0;
    allPairNrTot=0;
    sectorPairNrTot=0;
    maxFiredTotalPads  = 3000; // upper limit of number of fired pads
    
    maxCols = getGeometryPar()->getColumns();
    maxRows = getGeometryPar()->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(!getGeometryPar()->getPadsPar()->getPad(col+padOffset)->getPadActive() && col<maxCols) ++col;
       if(col == maxCols)
       {
          maxRows=row;
          break;
       }
       pLeftBorder[row]=col;
       while(getGeometryPar()->getPadsPar()->getPad(col+padOffset)->getPadActive() && col<maxCols) ++col;
       pRightBorder[row]=col-1;
    }
    maxPads = maxRows * maxCols;
    
    
    // now creating pads array
    
    Int_t i;
    Int_t j;
    if (pPads) {
     for (j = 0; j < 6; j++) 
        if (pPads[j]) 
           delete [] pPads[j];
     delete [] pPads;
    } 

    pPads = new HRichPadSignal * [6];
    for (j = 0; j < 6; j++) {
     pSectorPads = pPads[j] = new HRichPadSignal[maxPads];
     for (i = 0; i < maxPads; pSectorPads[i++].clear());
    }
    pSectorPads=pPads[iActiveSector];
    
    return kTRUE;
}
//============================================================================

//----------------------------------------------------------------------------
Int_t
HRichAnalysisSim::getPadsNTrack1(Int_t padx, Int_t pady, Int_t sec)
{
// This function returns for each pad the corresponding
// NTrack1 value (cfr. HRichCalSim, HRichDigitizer::execute).
// All track numbers are stored during digitization  in the
// catRichTrack Linear Category.
// NTrack1 is the index in the Track array corresponding to the
// first track for each pad.  
//
// This function is called from HRichRingFindSim::CalcRingParamters.

   HLocation loc1;
   loc1.set(3, sec, pady, padx);

   pRichCalSim = (HRichCalSim*)((HMatrixCategory*)getCalCat())->getObject(loc1);
   if ( NULL != pRichCalSim )
      return pRichCalSim->getNTrack1();   
   return 0;
}
//============================================================================

//----------------------------------------------------------------------------
Int_t
HRichAnalysisSim::getPadsNTrack2()
{
// This functions returns NTrack2, which is the index in 
// the Track Array corresponding to the last track for each pad.
// Pad must be identified by calling getPadsNTrack1(...) first. 

   if ( NULL != pRichCalSim)
      return pRichCalSim->getNTrack2();
   return 0;
}
//============================================================================

//----------------------------------------------------------------------------
Int_t
HRichAnalysisSim::getFlag(Int_t index)
{
// It is called from HRichRingFindSim::LookForTrack().
// This function returns the flag contained in the
// catRichTrack container at the position index.
// This flag is 0 for Chrenkov photons and 1 for IP.

    HRichTrack *trk = NULL;

    trk = (HRichTrack*)((HLinearCategory*)m_pTrackCat)->getObject(index);
    if ( NULL != trk )
       return trk->getFlag();
    return  -1;
}
//============================================================================

//----------------------------------------------------------------------------
Int_t
HRichAnalysisSim::getTrack(Int_t index)
{
//  This function returns the track number contained in the
//  catRichTrack container at the position index.

   HRichTrack *trk = NULL;

   trk = ((HRichTrack*)((HLinearCategory*)m_pTrackCat)->getObject(index));
   if( NULL != trk )
      return trk->getTrack();
   return -10;

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

//----------------------------------------------------------------------------
Int_t
HRichAnalysisSim::execute()
{
 
   HRichCal*  pCal           = NULL;;
   Int_t      iRingNrTotEvt  = 0;
   Int_t      fPadFiredTot   = 0;
   Int_t      sector         = -1;
   Int_t      padnr          = -1;
   Int_t      j;
   Int_t      ampl;
  
   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 = static_cast<HRichCal*>(fIter->Next())) )
   {
      if ( (ampl = (Int_t)pCal->getCharge()) > 0 )
      {
         fPadFired[pCal->getSector()]++;
         padnr = pCal->getCol() + pCal->getRow() * maxCols;
         pPads[pCal->getSector()][padnr].setAmplitude(ampl);  
      }
   }

   fPadFiredTot = 0;
   for ( sector = 0; sector < 6; ++sector )
   {
      fPadFiredTot += fPadFired[sector];
      if ( fPadFired[sector] >= fpAnalysisPar->maxFiredSectorPads )
      {
         Warning("execute", "To many fired pads in sector %i. %i/%i",
                 sector, fPadFired[sector], fpAnalysisPar->maxFiredSectorPads);
         fIter->Reset();
         while( NULL != (pCal = static_cast<HRichCal*>(fIter->Next())) )
         {
            if ( pCal->getSector() == sector )
               pCal->setIsCleanedSector( kTRUE );
         }
         return 0;
      }
   }

   if( fPadFiredTot > maxFiredTotalPads )
   {
      Warning("execute", "Analysis of event skipped: too many pads fired %i/%i",
              fPadFiredTot, maxFiredTotalPads);
      fIter->Reset();
      while( NULL != (pCal = static_cast<HRichCal*>(fIter->Next())) )
      {
         pCal->setIsCleanedSector( kTRUE );
      }
      return 0;
   }

   // **************************************************************
   // ------- loop over sectors --- begin ---
   for ( sector = 0; sector < 6; ++sector ) 
      if ( getGeometryPar()->getSectorActive(sector) > 0 )
      {
	  
         Reset();
         SetActiveSector(sector);
         iPadFiredNr = fPadFired[sector];

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

         if ( 0 == getAnalysisPar()->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;
	    pLabelArea[0].iSignature = 0;
         }
	  
         iRingNr = pRingFindSim->Execute(this);
         iRingNrTot    += iRingNr;
         iRingNrTotEvt += iRingNr;
         updateHeaders( sector, gHades->getCurrentEvent()->getHeader()->getEventSeqNumber() );
         updateHits(sector);
      }
   // ------- loop over sectors --- end ---

   fIter->Reset();
   while( NULL != (pCal = static_cast<HRichCal*>(fIter->Next())) )
   {
      padnr = pCal->getCol() + pCal->getRow() * maxCols;
      pCal->setIsCleanedSingle(pPads[pCal->getSector()][padnr].getIsCleanedSingle());
      pCal->setIsCleanedHigh  (pPads[pCal->getSector()][padnr].getIsCleanedHigh());
   }

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

//----------------------------------------------------------------------------
void
HRichAnalysisSim::updateHits(Int_t nSec)
{
// Each Hit is stored, in addition to real data for simulated one
// the tracks number corresponding to each ring are stored too.
// Maximal 3 tracks are stored for each ring. The track numbers
// are sorted according to their weights.

  HRichHitSim *hit=NULL;
  HLocation    loc;

  for (Int_t i = 0; i < iRingNr; i++) {

    loc.set(1, nSec);

    hit=(HRichHitSim *)m_pHitCat->getNewSlot(loc,NULL);
    if (hit!=NULL) {
      hit=new(hit) HRichHitSim;

      HRichHitSim & ring = pRings[i];
 
      *hit = ring;
      hit->setSector(nSec);

      HRichPad * pad = getGeometryPar()->getPadsPar()->
                         getPad(ring.getRingCenterX(),ring.getRingCenterY());
      if(pad){
	hit->setTheta(pad-> getTheta() );
	hit->setPhi(pad->getPhi(nSec));
      }
      else {
	hit->fPhi = -10;
	hit->fTheta = -10;
      }
      Int_t k = 0;
      while(ring.iRingTrack[k]) {

	if (hit->track1==0) {
	  hit->track1 = ring.iRingTrack[k];
	  hit->flag1 =  ring.iRingFlag[k];
	  (hit->weigTrack1)++;
	} else{
	  if(ring.iRingTrack[k]== hit->track1){
	    (hit->weigTrack1)++;

	  } else {
	    if(hit->track2==0) {
	      hit->track2 = ring.iRingTrack[k];
	      hit->flag2 =  ring.iRingFlag[k];
	      (hit->weigTrack2)++;
	    } else{
	      if(ring.iRingTrack[k]== hit->track2){
		(hit->weigTrack2)++;

	      } else {
		if (hit->track3==0) {
		  hit->track3 = ring.iRingTrack[k];
		  hit->flag3 =  ring.iRingFlag[k];
		  (hit->weigTrack3)++;
		}
		else (hit->weigTrack3)++;
	      }
	    }
	  }  
	}
	k++;

      }//end loop on track array

      sortTracks(hit);

    }

  }//end loop on rings

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

//----------------------------------------------------------------------------
void
HRichAnalysisSim::sortTracks(HRichHitSim *hit)
{
// Sorting tracks according to their weights
// Track with the highest weight is the first one

  Int_t tmp[3];
  tmp[0] = hit->track1;
  tmp[1] = hit->flag1;
  tmp[2] = hit->weigTrack1;

  if (hit->weigTrack1 < hit->weigTrack2) { // 2>1
    if (hit->weigTrack2 > hit->weigTrack3) { // 2 is the largest
      hit->track1     = hit->track2;
      hit->flag1      = hit->flag2;
      hit->weigTrack1 = hit->weigTrack2;
      if (tmp[2] < hit->weigTrack3) { // 2>3>1
	hit->track2     = hit->track3;
	hit->flag2      = hit->flag3;
	hit->weigTrack2 = hit->weigTrack3;
	hit->track3     = tmp[0];
	hit->flag3      = tmp[1];
	hit->weigTrack3 = tmp[2];
      }
      else { // 2>1>=3
	hit->track2     = tmp[0];
	hit->flag2      = tmp[1];
	hit->weigTrack2 = tmp[2];
      }
    }
    else { // 3>= 2; 2>1 
      hit->track1     = hit->track3;
      hit->flag1      = hit->flag3;
      hit->weigTrack1 = hit->weigTrack3;
      hit->track3     = tmp[0];
      hit->flag3      = tmp[1];
      hit->weigTrack3 = tmp[2];
    }
  }
  else { // 1>=2
    if (hit->weigTrack1 < hit->weigTrack3) { // 3>1>=2
      hit->track1     = hit->track3;
      hit->flag1      = hit->flag3;
      hit->weigTrack1 = hit->weigTrack3;
      hit->track3     = hit->track2;
      hit->flag3      = hit->flag2;
      hit->weigTrack3 = hit->weigTrack2;
      hit->track2     = tmp[0];
      hit->flag2      = tmp[1];
      hit->weigTrack2 = tmp[2];
    }
    else { // 1>=2; 1>=3
      if (hit->weigTrack2 < hit->weigTrack3) { // 1>=3>2
	tmp[0] = hit->track2;
	tmp[1] = hit->flag2;
	tmp[2] = hit->weigTrack2;
	hit->track2     = hit->track3;
	hit->flag2      = hit->flag3;
	hit->weigTrack2 = hit->weigTrack3;
	hit->track3     = tmp[0];
	hit->flag3      = tmp[1];
	hit->weigTrack3 = tmp[2];
      }
    } // 1>=2>=3
  }
}