//-*- Mode: C++ -*-
// $Id: AliHLTTPCCAGrid.h,v 1.1.1.1 2010/07/26 20:55:38 ikulakov Exp $
// ************************************************************************
// This file is property of and copyright by the ALICE HLT Project        *
// ALICE Experiment at CERN, All rights reserved.                         *
// See cxx source for full Copyright notice                               *
//                                                                        *
//*************************************************************************

#ifndef L1GRID_H
#define L1GRID_H

#include "CbmL1Def.h"
#include <algorithm>
#include <string.h>

#include <assert.h>
#include <cstdio>
#include "omp.h"
#include "L1StsHit.h"
#include "L1HitPoint.h"

using namespace std;

class L1Algo;
//class L1HitPoint;


/**
 * @class L1Grid
 *
 * 2-dimensional grid of pointers.
 * pointers to (y,z)-like objects are assigned to the corresponding grid bin
 * used by L1Tracker to speed-up the hit operations
 * grid axis are named Z,Y to be similar to TPC row coordinates.
 */
class L1Grid {
 public:
 L1Grid():
 fN(0), fNy(0), fNz(0), fNThreads(0){}
 L1Grid(const L1Grid& grid):
 fN(grid.N()), fNy(grid.Ny()), fNz(grid.Nz()), fNThreads(0){  }
 ~L1Grid(){ //if ( fFirstHitInBin ) delete[] fFirstHitInBin;
     
     for (int i=0; i<8000000; i++)
         omp_destroy_lock(&(lock[i]));
 
   for( int i = 0; i < fNThreads; i++ ) {
     if  (fFirstHitInBinArray[i]) delete[] fFirstHitInBinArray[i];
     if  (fFirstHitInBin[i]) delete[] fFirstHitInBin[i];
   }
 }
 
  void CreatePar( L1HitPoint* points, THitI nhits, vector <L1HitPoint>* pointsBuf, vector <L1StsHit>* hitsBuf, const L1StsHit* hits, vector <THitI>* indices, vector <THitI>* indicesBuf, char iS, L1Algo &Algo, THitI n);
    
      void CreateParSIMD( L1HitPoint* points, THitI nhits, vector <L1HitPoint>* pointsBuf, vector <L1StsHit>* hitsBuf, const L1StsHit* hits, vector <THitI>* indices, vector <THitI>* indicesBuf, char iS, L1Algo &Algo, THitI n, const vector< L1Strip > *vStsStrips, const vector< L1Strip > *vStsStripsB, const vector< fscal >   *vStsZPos);
  
    void CreatePar0( float yMin, float yMax, float zMin, float zMax, float sy, float sz);
    
     void HitsSort(L1HitPoint* pointsBuf, L1StsHit* hitsBuf, const L1StsHit* hits, THitI* indices, THitI* indicesBuf, L1HitPoint* points, THitI n, THitI nhits, char iS, L1Algo &Algo,  L1StsHitV* vStsHitsUnusedV, const vector< fscal >   *vStsZPos, const vector< L1Strip >  *vStsStrips, const vector< L1Strip >  *vStsStripsB);
    
    void HitsSortSIMD(L1HitPoint* pointsBuf, L1StsHit* hitsBuf, const L1StsHit* hits, THitI* indices, THitI* indicesBuf, L1HitPoint* points, THitI n, THitI nhits, char iS, L1Algo &Algo,  L1StsHitV* vStsHitsUnusedV, const vector< fscal >   *vStsZPos, const vector< L1Strip >  *vStsStrips, const vector< L1Strip >  *vStsStripsB);
  
  void Initial1( int NThreads );
    
    void Create( float yMin, float yMax, float zMin, float zMax, float sy, float sz  );
    
  void Fill( const L1HitPoint* points, THitI n ); // call after sort
    void FillPar( const L1HitPoint* points, THitI n);
    
  
  int GetBin( float Y, float Z ) const;
  
  unsigned int GetBinBounded( const float &Y, const float &Z ) const;
  int_v GetBinBounded( const fvec &Y, const fvec &Z ) const;
  void GetBinBounded( const float &Y, const float &Z, unsigned short &bY, unsigned short &bZ ) const;
  void GetBinBounded( const fvec &Y, const fvec &Z, int_v &bY, int_v &bZ ) const;

    // static unsigned short GetBinBounded( const L1Grid *array, const unsigned short &indexes, const float &Y, const float &Z );
    // static void GetBinBounded( const L1Grid *array, const unsigned short &indexes, const float &Y, const float &Z, unsigned short *bY, unsigned short *bZ );
    // static unsigned short Ny( const L1Grid *array, const unsigned short &indexes ) { return unsigned short( array, &L1Grid::fNy, indexes ); }

  void GetBinBounds( unsigned int iBin, float &Ymin, float &Ymax, float &Zmin, float &Zmax) const;
  void GetBinBounds( fvec iBin, fvec &Ymin, fvec &Ymax, fvec &Zmin, fvec &Zmax) const;
    
  unsigned int   N()        const { return fN;  }
  unsigned short Ny()       const { return fNy; }
  unsigned short Nz()       const { return fNz; }

  THitI FirstHitInBin(int i) const  {
   
    return fFirstHitInBinArray[i/fBinInGrid][i%fBinInGrid];
  }
  
    int_v FirstHitInBin(int_v &i) const  {
        
      
      int_v tempo(-1);
      
      for (int m=0; m<fvecLen; m++)
      {
        int k = int(i[m]);
     //   cout<<k<<endl;
     if (k<0) continue;
       
      tempo[m]=fFirstHitInBinArray[k/fBinInGrid][k%fBinInGrid];
      }
    //  tempo = fvec (0.f);
      
      return tempo;
    
  }

//     inline void UpdateGrid (fscal x, fscal y)
//     {
//       //  #pragma omp atomic
//         fNumberHitsInBin[GetBinBounded(x, y)]++;
// 
//     };
    

    
void UpdateIterGrid(int Nelements, L1StsHit* hits,  vector <THitI>* indicesBuf, THitI* indices,  vector <THitI>* indicesBuf2, vector <L1StsHit>* hits2, vector <L1HitPoint>* pointsBuf, L1HitPoint* points,  int &NHitsOnStation);
void fill(int Nelements, L1StsHit* hits,  vector <THitI>* indicesBuf);

  
 
 private:

  unsigned int   fN;       //* total N bins
  unsigned short fNy;      //* N bins in Y
  unsigned short fNz;      //* N bins in Z
  float fYMinOverStep;     //* minimal Y value * fStepYInv
  float fZMinOverStep;     //* minimal Z value * fStepZInv
  float fStepYInv; //* inverse bin size in X
  float fStepZInv; //* inverse bin size in Y
  int fBinInGrid;
  
   vector <THitI*> fFirstHitInBin;
   vector <THitI*> fFirstHitInBinArray;
   vector <THitI> fOffsets;
   vector <THitI> fNumberHitsInBin;
    
   unsigned short fNThreads;
   
   vector <omp_lock_t> lock;

};

// inline unsigned short L1Grid::GetBinBounded( const L1Grid *array, const unsigned short &indexes, const float &Y, const float &Z )
// {
  
//   const float fZMinOverStep( array, &L1Grid::fZMinOverStep, indexes );
//   const float fStepZInv( array, &L1Grid::fStepZInv, indexes );
//   const unsigned short fNz( array, &L1Grid::fNz, indexes );
//   short_v zBin = ( Z * fStepZInv - fZMinOverStep ).staticCast<short_v>();
//   unsigned short zBin2 = CAMath::Max( short_v( Vc::Zero ), CAMath::Min( short_v( fNz - 1 ), zBin ) ).staticCast<unsigned short>();

//   const float fYMinOverStep( array, &L1Grid::fYMinOverStep, indexes );
//   const float fStepYInv( array, &L1Grid::fStepYInv, indexes );
//   const unsigned short fNy( array, &L1Grid::fNy, indexes );
//   short_v yBin = ( Y * fStepYInv - fYMinOverStep ).staticCast<short_v>();
//   unsigned short yBin2 = CAMath::Max( short_v( Vc::Zero ), CAMath::Min( short_v( fNy - 1 ), yBin ) ).staticCast<unsigned short>();
//   return zBin2 * fNy + yBin2;
// }

// inline void L1Grid::GetBinBounded( const L1Grid *array, const unsigned short &indexes, const float &Y, const float &Z, unsigned short *bY, unsigned short *bZ )
// {
//   const float fYMinOverStep( array, &L1Grid::fYMinOverStep, indexes );
//   const float fStepYInv( array, &L1Grid::fStepYInv, indexes );
//   const unsigned short fNy( array, &L1Grid::fNy, indexes );
//   const short fNy2 = fNy.staticCast<short_v>();
//   const short &yBin =  static_cast<short>( Y * fStepYInv - fYMinOverStep );
//   *bY = static_cast<unsigned short>(max( 0, min( fNy2 - 1, yBin ) ));

//   const float fZMinOverStep( array, &L1Grid::fZMinOverStep, indexes );
//   const float fStepZInv( array, &L1Grid::fStepZInv, indexes );
//   const unsigned short fNz( array, &L1Grid::fNz, indexes );
//   const short_v fNz2 = fNz.staticCast<short_v>();
//   const short_v &zBin = ( Z * fStepZInv - fZMinOverStep ).staticCast<short_v>();
//   *bZ = CAMath::Max( short_v( Vc::Zero ), CAMath::Min( fNz2 - 1, zBin ) ).staticCast<unsigned short>();
// }
inline int L1Grid::GetBin( float Y, float Z ) const
{
  //* get the bin pointer
  const int& yBin = static_cast<int>( Y * fStepYInv - fYMinOverStep );
  const int& zBin = static_cast<int>( Z * fStepZInv - fZMinOverStep );
  assert( yBin >= 0 );
  assert( zBin >= 0 );
  assert( yBin < static_cast<int>( fNy ) );
  assert( zBin < static_cast<int>( fNz ) );
  const int& bin = zBin * fNy + yBin;
  return bin;
}

inline void L1Grid::GetBinBounds( unsigned int iBin, float &Ymin, float &Ymax, float &Zmin, float &Zmax) const
{
  int zBin = iBin / fNy;
  int yBin = iBin % fNy;
  Ymin = (fYMinOverStep + yBin)/fStepYInv;
  Zmin = (fZMinOverStep + zBin)/fStepZInv;
  Ymax = Ymin + 1. / fStepYInv;
  Zmax = Zmin + 1. / fStepZInv;
}

inline void L1Grid::GetBinBounds( fvec iBin, fvec &Ymin, fvec &Ymax, fvec &Zmin, fvec &Zmax) const
{
  
  
  int_v zBin = int_v(iBin)/int_v(int(fNy));
  int_v yBin = int_v(iBin) - zBin*((int_v)fNy);
  Ymin = (fYMinOverStep + yBin)/fStepYInv;
  Zmin = (fZMinOverStep + zBin)/fStepZInv;
  Ymax = Ymin + 1.f / fStepYInv;
  Zmax = Zmin + 1.f / fStepZInv;
}

inline unsigned int L1Grid::GetBinBounded( const float &Y, const float &Z ) const
{
  //* get the bin pointer at 

  unsigned short yBin, zBin;
  GetBinBounded( Y, Z, yBin, zBin );
  return (unsigned int)zBin * (unsigned int)fNy + (unsigned int)yBin;
}

inline void L1Grid::GetBinBounded( const float &Y, const float &Z, unsigned short &bY, unsigned short &bZ ) const
{
  const short &yBin = ( Y * fStepYInv - fYMinOverStep );
  const short &zBin = ( Z * fStepZInv - fZMinOverStep );

  bY = max( short(0), min( short( fNy - 1 ), yBin ) );
  bZ = max( short(0), min( short( fNz - 1 ), zBin ) );
}

inline void L1Grid::GetBinBounded( const fvec &Y, const fvec &Z, int_v &bY, int_v &bZ ) const
{
  const fvec &yBin = ( Y * fStepYInv - fYMinOverStep );
  const fvec &zBin = ( Z * fStepZInv - fZMinOverStep );
  
//   for (int i=0; i<fvecLen; i++)
//   {
//   bY[i] = max( short(0), min( short( fNy - 1 ), short(yBin[i]) ) );
//   bZ[i] = max( short(0), min( short( fNz - 1 ), short(zBin[i]) ) );
//       
// 
//   }

  bY = max( int_v(Vc::Zero), (min( int_v( fNy - 1 ), int_v(yBin) ) ));
  bZ = max( int_v(Vc::Zero), (min( int_v( fNz - 1 ), int_v(zBin) ) ));

}

inline int_v L1Grid::GetBinBounded( const fvec &Y, const fvec &Z ) const
{
  //* get the bin pointer at 

  int_v yBin, zBin;
  GetBinBounded( Y, Z, yBin, zBin );
  return zBin * int_v(fNy) + yBin;
}
#endif