#ifndef L1HitArea_H
#define L1HitArea_H

#include "L1Grid.h"
#include "L1Algo.h"

#include "CbmL1Def.h"

struct L1Grid;
class L1Row;
class L1SliceData;

class L1HitArea
{
  public:

    L1HitArea( const L1Grid & grid, float y, float z, float dy, float dz );

    /**
     * look up the next hit in the requested area.
     * Sets h to the coordinates and returns the index for the hit data
     */
    bool GetNext( THitI& i );
   

  protected:
    const L1Grid &fGrid;

    unsigned short fBZmax;   // maximal Z bin index
    unsigned short fBDY;     // Y distance of bin indexes
    unsigned int fIndYmin; // minimum index for
    unsigned short fIz;      // current Z bin index (incremented while iterating)
    THitI fHitYlst; // last possible hit index in current z-line
    THitI fIh;      // hit index iterating inside the bins
    int fNy;      // Number of bins in Y direction
};

class L1HitAreaV
{
  public:

    L1HitAreaV( const L1Grid & grid, fvec& y, fvec& z, fvec& dy, fvec& dz );

    /**
     * look up the next hit in the requested area.
     * Sets h to the coordinates and returns the index for the hit data
     */
void HitsIndexesMiddle(nsL1::vector <int>::TSimd &MiddleHitsIndexes, nsL1::vector <int>::TSimd &FirstHitsIndexes, int i1);
    void HitsIndexesMiddle(nsL1::vector <int>::TSimd &MiddleHitsIndexes, nsL1::vector <int>::TSimd &FirstHitsIndexes, int i1, int n1);
       void HitsIndexesMiddle(nsL1::vector <int>::TSimd &MiddleHitsIndexes, nsL1::vector <int>::TSimd &FirstHitsIndexes, int i1, int n1, int hit);
    
     void HitsIndexesMiddle(nsL1::vector <int>::TSimd &MiddleHitsIndexes, nsL1::vector <int>::TSimd &FirstHitsIndexes, int i1, int_v hits, int n1, L1TrackParV &TempTracks, L1TrackPar T);
    
    void HitsIndexesMiddle(nsL1::vector <int>::TSimd &MiddleHitsIndexes, nsL1::vector <int>::TSimd &FirstHitsIndexes, nsL1::vector <int>::TSimd &MiddleHitsIndexes2, int i1, int_v hits, int_v hitsM, int n1, L1TrackParV &TempTracks, L1TrackPar T);
    
    void HitsIndexesMiddle(nsL1::vector <int>::TSimd &MiddleHitsIndexes, nsL1::vector <int>::TSimd &FirstHitsIndexes, nsL1::vector <int>::TSimd &MiddleHitsIndexes2, nsL1::vector <int>::TSimd &FirstHitsIndexes2, int i1, int_v hits, int_v hitsM, int_v hitsL3, int n1, L1TrackParV &TempTracks, L1TrackPar T);
    
    void HitsIndexesMiddle(nsL1::vector <int>::TSimd &MiddleHitsIndexes, nsL1::vector <int>::TSimd &FirstHitsIndexes, nsL1::vector <int>::TSimd &MiddleHitsIndexes2, nsL1::vector <int>::TSimd &FirstHitsIndexes2, int i1, int_v hits, int_v hitsM, int_v hitsL3, int n1);
    
        void HitsIndexesMiddle(nsL1::vector <int>::TSimd &MiddleHitsIndexes, nsL1::vector <int>::TSimd &FirstHitsIndexes, nsL1::vector <int>::TSimd &MiddleHitsIndexes2, nsL1::vector <int>::TSimd &FirstHitsIndexes2, int i1, int_v hits, int_v hitsM, int_v hitsL3);
void HitsIndexesMiddle(int i1, int_v hits, int_v hitsM, int_v hitsL3, int& fill, int& fill2, L1TrackPar& T, L1HitPointV& Points, float& TRIPLET_CHI2_CUT, L1Station &star, int_m& IsNiceDoublet, int& n3, nsL1::vector <L1TrackPar>::TSimd &TripletsArrayV, nsL1::vector<fvec>::TSimd &u_front_3V, nsL1::vector<fvec>::TSimd &u_back_3V, nsL1::vector<fvec>::TSimd &z_Pos_3V, int_v& hitR_V, int_v& hitM_V, int_v& hitL_V3, int_v& hitL_V, int_v& hitR_Vl, int_v& hitM_Vl, int_v& hitL_V3l, int_v& hitL_Vl, L1TrackPar& T2, L1TrackPar& T3, fvec& u, fvec& v, fvec& z);
void HitsIndexesMiddleLast(int i1, int_v hits, int_v hitsM, int_v hitsL3, int& fill, int& fill2, L1TrackPar& T, L1HitPointV& Points, float& TRIPLET_CHI2_CUT, L1Station &star, int_m& IsNiceDoublet, int& n3, nsL1::vector <L1TrackPar>::TSimd &TripletsArrayV, nsL1::vector<fvec>::TSimd &u_front_3V, nsL1::vector<fvec>::TSimd &u_back_3V, nsL1::vector<fvec>::TSimd &z_Pos_3V, int_v& hitR_V, int_v& hitM_V, int_v& hitL_V3, int_v& hitL_V, int_v& hitR_Vl, int_v& hitM_Vl, int_v& hitL_V3l, int_v& hitL_Vl, L1TrackPar& T2, L1TrackPar& T3, fvec& u, fvec& v, fvec& z);

    
    const L1Grid &fGrid;
//     fvec fBZmax;   // maximal Z bin index
//     fvec fBDY;     // Y distance of bin indexes
//     fvec fIndYmin; // minimum index for
//     fvec fIz;      // current Z bin index (incremented while iterating)
//     fvec fHitYlst; // last possible hit index in current z-line
//     fvec fIh;      // hit index iterating inside the bins
//     fvec fNy;      // Number of bins in Y direction

     int_v fFirstBin;
     int_v fNy;
     int_v fNx;
     int_v fGap;
     int_v fLastBin;
//   const fvec minY = y - dy;
//   const fvec maxY = y + dy;
//   const fvec minZ = z - dz;
//   const fvec maxZ = z + dz;
//   fvec bYmin, bZmin, bYmax; // boundary bin indexes
//   fGrid.GetBinBounded( minY, minZ, bYmin, bZmin );
//   fGrid.GetBinBounded( maxY, maxZ, bYmax, fBZmax );
// 
//   fBDY = ( bYmax - bYmin + 1 ); // bin index span in y direction
// 
//   fIndYmin = ( bZmin * fNy + bYmin ); // same as grid.GetBin(fMinY, fMinZ), i.e. the smallest bin index of interest
//   // fIndYmin + fBDY - 1 then is the largest bin index of interest with the same Z
// 
//   fGrid.GetBinBounds( fIndYmin, y, dy, z, dz );
// 
//   fIz = bZmin;
// 
//   fIh = fGrid.FirstHitInBin( fIndYmin );
//   fHitYlst = fGrid.FirstHitInBin( fIndYmin + fBDY );

};

inline L1HitAreaV::L1HitAreaV( const L1Grid & grid, fvec& y, fvec& z, fvec& dy, fvec& dz )
  : fGrid( grid ),
  fFirstBin(0),
  fNx( 0),
  fNy( 0),
  fGap(0),
  fLastBin(0)
{
  const fvec& minX = y - dy;
  const fvec& maxX = y + dy;
  const fvec& minY = z - dz;
  const fvec& maxY = z + dz;

 
  fFirstBin = fGrid.GetBinBounded( minX, minY);
  fNx = fGrid.GetBinBounded( maxX, minY)-fFirstBin+1;
  fGap = fGrid.Ny()-fNx;
  fLastBin=fGrid.GetBinBounded( maxX, maxY);
  fNy = (fLastBin+1+fGap-fFirstBin)/(fGrid.Ny());




}

inline L1HitArea::L1HitArea( const L1Grid & grid, float y, float z, float dy, float dz )
  : fGrid( grid ),
  fHitYlst( 0.f ),
  fIh( 0.f),
  fNy( fGrid.Ny() )
{
  const float minY = y - dy;
  const float maxY = y + dy;
  const float minZ = z - dz;
  const float maxZ = z + dz;
  unsigned short bYmin, bZmin, bYmax; // boundary bin indexes
  fGrid.GetBinBounded( minY, minZ, bYmin, bZmin );
  fGrid.GetBinBounded( maxY, maxZ, bYmax, fBZmax );

  fBDY = ( bYmax - bYmin + 1 ); // bin index span in y direction

  fIndYmin = ( bZmin * fNy + bYmin ); // same as grid.GetBin(fMinY, fMinZ), i.e. the smallest bin index of interest
  // fIndYmin + fBDY - 1 then is the largest bin index of interest with the same Z

//  fGrid.GetBinBounds( fIndYmin, y, dy, z, dz );

  fIz = bZmin;

  fIh = fGrid.FirstHitInBin( fIndYmin );
  fHitYlst = fGrid.FirstHitInBin( fIndYmin + fBDY );


  // const unsigned int hitLst = fGrid.FirstHitInBin( fBZmax * fNy + bZmin + 1 );
  // cout << fIh << " " << hitLst << endl;
  // if ( fIh >= hitLst ) fIz = fBZmax + 1; // finish area in one step

  // cout // << fBZmax << " " << fBDY << " " << fIndYmin << " " << fIz << " " << fHitYlst << " " << fIh << endl; // dbg 
  //  << "HitArea created:\n"
  //  << y << " " << z << endl
  //   << "bYmin:    " << bYmin << "\n"
  //   << "bZmin:    " << bZmin << "\n"
  //   << "bYmax:    " << bYmax << "\n"
  //   << "fBZmax:   " << fBZmax << "\n"
  //   << "fBDY:     " << fBDY << "\n"
  //   << "fIndYmin: " << fIndYmin << "\n"
  //   << "fIz:      " << fIz << "\n"
  //   << "fHitYlst: " << fHitYlst << "\n"
  //   << "fIh:      " << fIh << "\n"
  //   << "fNy:      " << fNy << "\n"
  //   << "Nz " << fGrid.Nz() << std::endl;
  L1_ASSERT ( fIndYmin + fBDY < fGrid.N()+1, fIndYmin + fBDY << " < " << fGrid.N());
}


inline void L1HitAreaV::HitsIndexesMiddle(nsL1::vector <int>::TSimd &MiddleHitsIndexes, nsL1::vector <int>::TSimd &FirstHitsIndexes, int i1)
{
  int nyMax = fNy.max();


 //cout<<fFirstBin<<" firstHitInLine "<<endl;

  for(int iNy=1; iNy<=nyMax; iNy++)
  {
    int_v firstHitInLine = fGrid.FirstHitInBin(fFirstBin);
    int_v LastBinInLine = fFirstBin+fNx;
    int_v lastHitInLine  = fGrid.FirstHitInBin(LastBinInLine);

    for(int iV=0; iV<fvecLen; iV++)
    {
     // if((i1 + iV)>= n1)
       // continue;

      if(iNy<=fNy[iV])
      {
        for(int iHit=firstHitInLine[iV]; iHit<lastHitInLine[iV]; iHit++)
        {


          MiddleHitsIndexes.push_back(iHit);
            
          
          FirstHitsIndexes.push_back(iV+i1);
        }
      }
    }

    fFirstBin(iNy<fNy) += (fGap+fNx);
  }

/*
    float_m IsFullHitVec (true);
    if((i1 + fvecLen)> n1)
      IsFullHitVec &= float_m( int_v(Vc::IndexesFromZero) < int(n1%fvecLen) );*/

}


inline void L1HitAreaV::HitsIndexesMiddle(nsL1::vector <int>::TSimd &MiddleHitsIndexes, nsL1::vector <int>::TSimd &FirstHitsIndexes, int i1, int n1)
{
  int nyMax = fNy.max();


 //cout<<fFirstBin<<" firstHitInLine "<<endl;

  for(int iNy=1; iNy<=nyMax; iNy++)
  {
    int_v firstHitInLine = fGrid.FirstHitInBin(fFirstBin);
    int_v LastBinInLine = fFirstBin+fNx;
    int_v lastHitInLine  = fGrid.FirstHitInBin(LastBinInLine);

    for(int iV=0; iV<fvecLen; iV++)
    {
     if((i1 + iV)>= n1)
        break;

      if(iNy<=fNy[iV])
      {
        for(int iHit=firstHitInLine[iV]; iHit<lastHitInLine[iV]; iHit++)
        {


          MiddleHitsIndexes.push_back(iHit);
            
          
          FirstHitsIndexes.push_back(iV+i1);
        }
      }
    }

    fFirstBin(iNy<fNy) += (fGap+fNx);
  }

/*
    float_m IsFullHitVec (true);
    if((i1 + fvecLen)> n1)
      IsFullHitVec &= float_m( int_v(Vc::IndexesFromZero) < int(n1%fvecLen) );*/

}



inline void L1HitAreaV::HitsIndexesMiddle(nsL1::vector <int>::TSimd &MiddleHitsIndexes, nsL1::vector <int>::TSimd &FirstHitsIndexes, int i1, int n1, int hit)
{
    int nyMax = fNy.max();
    
    
    //cout<<fFirstBin<<" firstHitInLine "<<endl;
    
    for(int iNy=1; iNy<=nyMax; iNy++)
    {
        int_v firstHitInLine = fGrid.FirstHitInBin(fFirstBin);
        int_v LastBinInLine = fFirstBin+fNx;
        int_v lastHitInLine  = fGrid.FirstHitInBin(LastBinInLine);
        
        for(int iV=0; iV<fvecLen; iV++)
        {
            if((i1 + iV)>= n1)
                continue;
            
            if(iNy<=fNy[iV])
            {
                for(int iHit=firstHitInLine[iV]; iHit<lastHitInLine[iV]; iHit++)
                {
                    
                    
                    MiddleHitsIndexes.push_back(iHit);
                    
                    
                    FirstHitsIndexes.push_back(iV+i1+hit);
                }
            }
        }
        
        fFirstBin(iNy<fNy) += (fGap+fNx);
    }
    
    /*
     float_m IsFullHitVec (true);
     if((i1 + fvecLen)> n1)
     IsFullHitVec &= float_m( int_v(Vc::IndexesFromZero) < int(n1%fvecLen) );*/
    
}


inline void L1HitAreaV::HitsIndexesMiddle(nsL1::vector <int>::TSimd &MiddleHitsIndexes, nsL1::vector <int>::TSimd &FirstHitsIndexes, int i1, int_v hits, int n1, L1TrackParV &TempTracks, L1TrackPar T)
{
    int nyMax = fNy.max();
    
    
    //cout<<fFirstBin<<" firstHitInLine "<<endl;
    
    for(int iNy=1; iNy<=nyMax; iNy++)
    {
        int_v firstHitInLine = fGrid.FirstHitInBin(fFirstBin);
        int_v LastBinInLine = fFirstBin+fNx;
        int_v lastHitInLine  = fGrid.FirstHitInBin(LastBinInLine);
        
        for(int iV=0; iV<fvecLen; iV++)
        {
            if((i1 + iV)>= n1)
                continue;
            
            if(iNy<=fNy[iV])
            {
                for(int iHit=firstHitInLine[iV]; iHit<lastHitInLine[iV]; iHit++)
                {
                    
                    TempTracks.Set(MiddleHitsIndexes.size(),T, iV);
                    MiddleHitsIndexes.push_back(iHit);
                    
                   
                    
                    FirstHitsIndexes.push_back(hits[iV]);
                   
                    
                     
                }
            }
        }
        
        fFirstBin(iNy<fNy) += (fGap+fNx);
    }
    
    /*
     float_m IsFullHitVec (true);
     if((i1 + fvecLen)> n1)
     IsFullHitVec &= float_m( int_v(Vc::IndexesFromZero) < int(n1%fvecLen) );*/
    
}



inline void L1HitAreaV::HitsIndexesMiddle(nsL1::vector <int>::TSimd &MiddleHitsIndexes, nsL1::vector <int>::TSimd &FirstHitsIndexes, nsL1::vector <int>::TSimd &MiddleHitsIndexes2, int i1, int_v hits, int_v hitsM, int n1, L1TrackParV &TempTracks, L1TrackPar T)
{
    int nyMax = fNy.max();
    
    
    //cout<<fFirstBin<<" firstHitInLine "<<endl;
    
    for(int iNy=1; iNy<=nyMax; iNy++)
    {
        int_v firstHitInLine = fGrid.FirstHitInBin(fFirstBin);
        int_v LastBinInLine = fFirstBin+fNx;
        int_v lastHitInLine  = fGrid.FirstHitInBin(LastBinInLine);
        
        for(int iV=0; iV<fvecLen; iV++)
        {
            if((i1 + iV)>= n1)
                continue;
            
            if(iNy<=fNy[iV])
            {
                for(int iHit=firstHitInLine[iV]; iHit<lastHitInLine[iV]; iHit++)
                {
                    
                    TempTracks.Set(MiddleHitsIndexes.size(),T, iV);
                    MiddleHitsIndexes.push_back(iHit);
                    MiddleHitsIndexes2.push_back(hitsM[iV]);

                    
                    
                    FirstHitsIndexes.push_back(hits[iV]);
                    
                    
                    
                }
            }
        }
        
        fFirstBin(iNy<fNy) += (fGap+fNx);
    }
    
    /*
     float_m IsFullHitVec (true);
     if((i1 + fvecLen)> n1)
     IsFullHitVec &= float_m( int_v(Vc::IndexesFromZero) < int(n1%fvecLen) );*/
    
}

inline void L1HitAreaV::HitsIndexesMiddle(nsL1::vector <int>::TSimd &MiddleHitsIndexes, nsL1::vector <int>::TSimd &FirstHitsIndexes, nsL1::vector <int>::TSimd &MiddleHitsIndexes2, nsL1::vector <int>::TSimd &FirstHitsIndexes2, int i1, int_v hits, int_v hitsM, int_v hitsL3, int n1, L1TrackParV &TempTracks, L1TrackPar T)
{
    int nyMax = fNy.max();
    
    
    //cout<<fFirstBin<<" firstHitInLine "<<endl;
    
    for(int iNy=1; iNy<=nyMax; iNy++)
    {
        int_v firstHitInLine = fGrid.FirstHitInBin(fFirstBin);
        int_v LastBinInLine = fFirstBin+fNx;
        int_v lastHitInLine  = fGrid.FirstHitInBin(LastBinInLine);
        
        for(int iV=0; iV<fvecLen; iV++)
        {
            if((i1 + iV)>= n1)
                continue;
            
            if(iNy<=fNy[iV])
            {
                for(int iHit=firstHitInLine[iV]; iHit<lastHitInLine[iV]; iHit++)
                {
                    
                    TempTracks.Set(MiddleHitsIndexes.size(),T, iV);
                    MiddleHitsIndexes.push_back(iHit);
                    MiddleHitsIndexes2.push_back(hitsM[iV]);
                    FirstHitsIndexes2.push_back(hitsL3[iV]);
                    
                    
                    FirstHitsIndexes.push_back(hits[iV]);
                    
                    
                    
                }
            }
        }
        
        fFirstBin(iNy<fNy) += (fGap+fNx);
    }
    
    
}

inline void L1HitAreaV::HitsIndexesMiddle(nsL1::vector <int>::TSimd &MiddleHitsIndexes, nsL1::vector <int>::TSimd &FirstHitsIndexes, nsL1::vector <int>::TSimd &MiddleHitsIndexes2, nsL1::vector <int>::TSimd &FirstHitsIndexes2, int i1, int_v hits, int_v hitsM, int_v hitsL3, int n1)
{
    int nyMax = fNy.max();
    
    
    //cout<<fFirstBin<<" firstHitInLine "<<endl;
    
    for(int iNy=1; iNy<=nyMax; iNy++)
    {
        int_v firstHitInLine = fGrid.FirstHitInBin(fFirstBin);
        int_v LastBinInLine = fFirstBin+fNx;
        int_v lastHitInLine  = fGrid.FirstHitInBin(LastBinInLine);
        
        for(int iV=0; iV<fvecLen; iV++)
        {
            if((i1 + iV)>= n1)
                continue;
            
            if(iNy<=fNy[iV])
            {
                for(int iHit=firstHitInLine[iV]; iHit<lastHitInLine[iV]; iHit++)
                {
                    
                  
                    MiddleHitsIndexes.push_back(iHit);
                    MiddleHitsIndexes2.push_back(hitsM[iV]);
                    FirstHitsIndexes2.push_back(hitsL3[iV]);
                    FirstHitsIndexes.push_back(hits[iV]);
                    
                    
                    
                }
            }
        }
        
        fFirstBin(iNy<fNy) += (fGap+fNx);
    }
    
    
}

inline void L1HitAreaV::HitsIndexesMiddle(nsL1::vector <int>::TSimd &MiddleHitsIndexes, nsL1::vector <int>::TSimd &FirstHitsIndexes, nsL1::vector <int>::TSimd &MiddleHitsIndexes2, nsL1::vector <int>::TSimd &FirstHitsIndexes2, int i1, int_v hits, int_v hitsM, int_v hitsL3)
{
    int nyMax = fNy.max();
    
    
    //cout<<fFirstBin<<" firstHitInLine "<<endl;
    
    for(int iNy=1; iNy<=nyMax; iNy++)
    {
        int_v firstHitInLine = fGrid.FirstHitInBin(fFirstBin);
        int_v LastBinInLine = fFirstBin+fNx;
        int_v lastHitInLine  = fGrid.FirstHitInBin(LastBinInLine);
        
        for(int iV=0; iV<fvecLen; iV++)
        {
            
            if(iNy<=fNy[iV])
            {
                for(int iHit=firstHitInLine[iV]; iHit<lastHitInLine[iV]; iHit++)
                {
                    
                    
                    MiddleHitsIndexes.push_back(iHit);
                    MiddleHitsIndexes2.push_back(hitsM[iV]);
                    FirstHitsIndexes2.push_back(hitsL3[iV]);
                    FirstHitsIndexes.push_back(hits[iV]);
                    
                    
                    
                }
            }
        }
        
        fFirstBin(iNy<fNy) += (fGap+fNx);
    }
    
    
}



inline void L1HitAreaV::HitsIndexesMiddle(int i1, int_v hits, int_v hitsM, int_v hitsL3, int& fill, int& fill2, L1TrackPar& T, L1HitPointV& Points, float& TRIPLET_CHI2_CUT, L1Station &star, int_m& IsNiceDoublet, int& n3, nsL1::vector <L1TrackPar>::TSimd &TripletsArrayV, nsL1::vector<fvec>::TSimd &u_front_3V, nsL1::vector<fvec>::TSimd &u_back_3V, nsL1::vector<fvec>::TSimd &z_Pos_3V, int_v& hitR_V, int_v& hitM_V, int_v& hitL_V3, int_v& hitL_V, int_v& hitR_Vl, int_v& hitM_Vl, int_v& hitL_V3l, int_v& hitL_Vl, L1TrackPar& T2, L1TrackPar& T3, fvec& u, fvec& v, fvec& z)
{



    int nyMax = fNy.max();
    
    for(int iNy=1; iNy<=nyMax; iNy++)
    {
        int_v firstHitInLine = fGrid.FirstHitInBin(fFirstBin);
        int_v LastBinInLine = fFirstBin+fNx;
        int_v lastHitInLine  = fGrid.FirstHitInBin(LastBinInLine);
        
        for(int iV=0; iV<fvecLen; iV++)
        {
            
            if(iNy<=fNy[iV])
            {
                for(int iHit=firstHitInLine[iV]; iHit<lastHitInLine[iV]; iHit++)
                {
                    
                    if (!(IsNiceDoublet[iV])) continue;
                    hitR_V[fill]=(iHit);
                    hitM_V[fill]=(hitsM[iV]);
                    hitL_V3[fill]=(hitsL3[iV]);
                    hitL_V[fill]=(hits[iV]);
                    T2.SetOneEntry(fill, T, iV);
                    
                    //cout<<T.x[iV]<<" "<<fill<<" fill o7 "<<iV<<endl;
                    
                    fill++;



             if (fill==fvecLen)  {

        
            
            fill = 0;
            
            fvec zr(Points.Z(), hitR_V);
            fvec yr(Points.Y(), hitR_V);
            fvec xr(Points.X(), hitR_V);
            fvec ur(Points.U(), hitR_V);
            fvec vr(Points.V(), hitR_V);
            
            
            fvec Pick_r22 = (TRIPLET_CHI2_CUT - T2.chi2);
            
            fvec y, C11;
            L1ExtrapolateYC11Line( T2, zr, y, C11 );
            

            const fvec dy_est2 = (Pick_r22*(abs(C11 + star.XYInfo.C11))); // TODO for FastPrim dx < dy - other sort is optimal. But not for doublets
            const fvec dy = yr - y;
            const fvec dy2 = dy*dy;
            
                 int_m IsNiceDoublet1 = int_m(true);
            
            IsNiceDoublet1 &=(!( dy2 > dy_est2 && dy < fvec(Vc::Zero) ));
            
            
            
            
            IsNiceDoublet1 &=(!(dy2 > dy_est2));
           if (!(IsNiceDoublet1.isEmpty())) {
            
            
            
            fvec x, C00;
            L1ExtrapolateXC00Line( T2, zr, x, C00 );
            
            const fvec dx_est2 = (Pick_r22*(abs(C00 + star.XYInfo.C00)));
            const fvec dx = xr - x;
            
            IsNiceDoublet1 &=(dx*dx <= dx_est2);
            
            
            if (!(IsNiceDoublet1.isEmpty())) {
            
            // check chi2  // not effective
            fvec C10;
            L1ExtrapolateC10Line( T2, zr, C10 );
            fvec chi2 = T2.chi2;
            L1FilterChi2XYC00C10C11( star.frontInfo, x, y, C00, C10, C11, chi2, ur );
            L1FilterChi2           ( star.backInfo,  x, y, C00, C10, C11, chi2, vr );
            
            
            
            IsNiceDoublet1 &=(chi2 <= TRIPLET_CHI2_CUT && C00 >= fvec(Vc::Zero) && C11 >= fvec(Vc::Zero));
           if (!(IsNiceDoublet.isEmpty())) {
   
           
            
            
            for (int iV=0; iV<fvecLen; iV++)
            {
                if (!(IsNiceDoublet1[iV])) continue;


                u[fill2]=ur[iV];
                v[fill2]=vr[iV];
                z[fill2]=zr[iV];
                hitR_Vl[fill2]=hitR_V[iV];
                hitL_Vl[fill2]=hitL_V[iV];
                hitM_Vl[fill2]=hitM_V[iV];
                T3.SetOneEntry(fill2, T2, iV);
               //  cout<<T2.x[iV]<<" "<<fill2<<" fill2 o7 "<<iV<<endl;
fill2++;

if (fill2==fvecLen)  {

        n3+=fvecLen;
            
            fill2 = 0;

TripletsArrayV.push_back(T3);
u_front_3V.push_back(u);
u_back_3V.push_back(v);
z_Pos_3V.push_back(z);
}
                
//               
//              //  T_buf.SetOneEntry(TrackParFilled, T2, iV);
//                
//                
//                RHitsL[TrackParFilled]=hitR_V[iV];
//                MHitsL[TrackParFilled]=hitM_V[iV];
//                LHitsL[TrackParFilled]=hitL_V3[iV];
//                u[TrackParFilled]=ur[iV];
//                v[TrackParFilled]=vr[iV];
//                z[TrackParFilled]=zr[iV];
//                TrackParFilled++;
//                
//                if (TrackParFilled==fvecLen) {
//                    
//                    RHits.resize(RHits.size()+fvecLen);
//                    MHitsV.resize(MHitsV.size()+fvecLen);
//                    LHitsV.resize(LHitsV.size()+fvecLen);
//                  //  TripletsArray.Set(T_buf);
//                   
//                   // TripletsArrayV.push_back(T_buf);
//                    u_front_3V.push_back(u);
//                    u_back_3V.push_back(v);
//                    z_Pos_3V.push_back(z);
//                    reinterpret_cast<int_v&>(RHits[Triplets*fvecLen])=RHitsL;
//                    reinterpret_cast<int_v&>(LHitsV[Triplets*fvecLen])=LHitsL;
//                    reinterpret_cast<int_v&>(MHitsV[Triplets*fvecLen])=MHitsL;
//                    
//                    TrackParFilled=0;
//                    Triplets++;
//                    
//                }
                
                
            }

}
}
}







}
                }
            }
        }
        
        fFirstBin(iNy<fNy) += (fGap+fNx);
    }
   
}
inline void L1HitAreaV::HitsIndexesMiddleLast(int i1, int_v hits, int_v hitsM, int_v hitsL3, int& fill, int& fill2, L1TrackPar& T, L1HitPointV& Points, float& TRIPLET_CHI2_CUT, L1Station &star, int_m& IsNiceDoublet, int& n3, nsL1::vector <L1TrackPar>::TSimd &TripletsArrayV, nsL1::vector<fvec>::TSimd &u_front_3V, nsL1::vector<fvec>::TSimd &u_back_3V, nsL1::vector<fvec>::TSimd &z_Pos_3V, int_v& hitR_V, int_v& hitM_V, int_v& hitL_V3, int_v& hitL_V, int_v& hitR_Vl, int_v& hitM_Vl, int_v& hitL_V3l, int_v& hitL_Vl, L1TrackPar& T2, L1TrackPar& T3, fvec& u, fvec& v, fvec& z)
{
    if (fill!=0)  {
        
        
        IsNiceDoublet&=(int_v(Vc::IndexesFromZero)<int_v(fill));
        
        fvec zr(Points.Z(), hitR_V);
        fvec yr(Points.Y(), hitR_V);
        fvec xr(Points.X(), hitR_V);
        fvec ur(Points.U(), hitR_V);
        fvec vr(Points.V(), hitR_V);
        
        
        fvec Pick_r22 = (TRIPLET_CHI2_CUT - T2.chi2);
        
        fvec y, C11;
        L1ExtrapolateYC11Line( T2, zr, y, C11 );
        
        
        const fvec dy_est2 = (Pick_r22*(abs(C11 + star.XYInfo.C11))); // TODO for FastPrim dx < dy - other sort is optimal. But not for doublets
        const fvec dy = yr - y;
        const fvec dy2 = dy*dy;
        
        
        
        IsNiceDoublet &=(!( dy2 > dy_est2 && dy < fvec(Vc::Zero) ));
        
        
        
        
        IsNiceDoublet &=(!(dy2 > dy_est2));
        if ((IsNiceDoublet.isEmpty())) return;
        
        
        
        fvec x, C00;
        L1ExtrapolateXC00Line( T2, zr, x, C00 );
        
        const fvec dx_est2 = (Pick_r22*(abs(C00 + star.XYInfo.C00)));
        const fvec dx = xr - x;
        
        IsNiceDoublet &=(dx*dx <= dx_est2);
        
        
        if ((IsNiceDoublet.isEmpty())) return;
        
        // check chi2  // not effective
        fvec C10;
        L1ExtrapolateC10Line( T2, zr, C10 );
        fvec chi2 = T2.chi2;
        L1FilterChi2XYC00C10C11( star.frontInfo, x, y, C00, C10, C11, chi2, ur );
        L1FilterChi2           ( star.backInfo,  x, y, C00, C10, C11, chi2, vr );
        
        
        
        IsNiceDoublet &=(chi2 <= TRIPLET_CHI2_CUT && C00 >= fvec(Vc::Zero) && C11 >= fvec(Vc::Zero));
        if ((IsNiceDoublet.isEmpty())) return;
        
        
        
        
        for (int iV=0; iV<fvecLen; iV++)
        {
            if (!(IsNiceDoublet[iV])) continue;
            
            
            u[fill2]=ur[iV];
            v[fill2]=vr[iV];
            z[fill2]=zr[iV];
            hitR_Vl[fill2]=hitR_V[iV];
            hitL_Vl[fill2]=hitL_V[iV];
            hitM_Vl[fill2]=hitM_V[iV];
            T3.SetOneEntry(fill2, T2, iV);
            fill2++;
            
            if (fill2==fvecLen)  {
                
                n3+=fvecLen;
                
                fill2 = 0;
                
                TripletsArrayV.push_back(T3);
                u_front_3V.push_back(u);
                u_back_3V.push_back(v);
                z_Pos_3V.push_back(z);
            }
            
            
        }
        
    }

    
    
}


inline bool L1HitArea::GetNext( THitI& i )
{

  bool yIndexOutOfRange = fIh >= fHitYlst;     // current y is not in the area
  bool nextZIndexOutOfRange = (fIz >= fBZmax);   // there isn't any new z-line

  if ( yIndexOutOfRange && nextZIndexOutOfRange ) { // all iterators are over the end
    return false;
  }
  
  // at least one entry in the vector has (fIh >= fHitYlst && fIz < fBZmax)
  bool needNextZ = yIndexOutOfRange && !nextZIndexOutOfRange;
    int num=0;
  // skip as long as fIh is outside of the interesting bin y-index
  while ( ISLIKELY( needNextZ ) ) {   //ISLIKELY to speed the programm and optimise the use of registers
    fIz++;   // get new z-line
   // num++;
      // get next hit
    fIndYmin += fNy;
    fIh = fGrid.FirstHitInBin( fIndYmin ); // get first hit in cell, if z-line is new
    fHitYlst = fGrid.FirstHitInBin( fIndYmin + fBDY );
    yIndexOutOfRange = fIh >= fHitYlst;
    nextZIndexOutOfRange = (fIz >= fBZmax);
    needNextZ = yIndexOutOfRange && !nextZIndexOutOfRange;
    
  }
  //cout<<num<<endl;
  L1_ASSERT ( fIh < fGrid.FirstHitInBin(fGrid.N()-1) || yIndexOutOfRange, fIh << " < " << fGrid.FirstHitInBin(fGrid.N()-1) << " || " <<  yIndexOutOfRange);
  i = fIh; // return
  fIh++; // go to next
  return !yIndexOutOfRange;
}





#endif