/*
 *=====================================================
 *
 *  CBM Level 1 Reconstruction
 *
 *  Authors: I.Kisel,  S.Gorbunov, I.Kulakov, M.Zyzak
 * Documentation: V.Akishina
 *
 *  e-mail : i.kulakov@gsi.de
 *
 *=====================================================
 *
 *  Finds tracks using the Cellular Automaton algorithm
 *
 */

#include "L1Algo.h"
#include "L1TrackPar.h"
#include "L1Branch.h"
#include "L1Track.h"
#include "L1Extrapolation.h"
#include "L1Filtration.h"
#include "L1AddMaterial.h"
#include "L1HitPoint.h"
//#include "TDHelper.h"
#include "L1Grid.h"
#include "L1HitArea.h"
#include "L1Portion.h"
#include "omp.h"
#include "pthread.h"
#include <parallel/algorithm>
#include "L1HitsSortHelper.h"



#include "L1Timer.h"

#ifdef DRAW
#include "utils/L1AlgoDraw.h"
#endif
#ifdef PULLS
#include "utils/L1AlgoPulls.h"
#endif
#ifdef TRIP_PERFORMANCE
#include "utils/L1AlgoEfficiencyPerformance.h"
#endif
#ifdef DOUB_PERFORMANCE
#include "utils/L1AlgoEfficiencyPerformance.h"
#endif

#ifdef TBB
#include "L1AlgoTBB.h"
#endif

#include <fstream>
#include <iostream>
#include <map>
#include <list>
#include <stdio.h>
#include <algorithm>




// using namespace std;
using std::cout;
using std::endl;
using std::vector;


          /// Prepare the portion of data of left hits of a triplet: all hits except the last and the second last station will be procesesed in the algorythm, 
      /// the data is orginesed in order to be used by SIMD
inline void L1Algo::f10(  // input
                Tindex start_lh, Tindex n1,  L1HitPoint *vStsHits_l,
                    // output
                  fvec *u_front, fvec *u_back, fvec *zPos,
                  THitI* hitsl_1
                  )
{
  const Tindex end_lh = start_lh+n1;
  for (Tindex ilh = start_lh, i1 = 0; ilh < end_lh; ++ilh, ++i1){
    Tindex i1_V= i1/fvecLen;
    Tindex i1_4= i1%fvecLen;
    L1HitPoint &hitl = vStsHits_l[ilh];
      
      
    hitsl_1[i1] = ilh;
    u_front[i1_V][i1_4] = hitl.U();
    u_back [i1_V][i1_4] = hitl.V();
    zPos   [i1_V][i1_4] = hitl.Z();
  }
}


/// Get the field approximation. Add the target to parameters estimation. Propagaete to the middle station of a triplet.
inline void L1Algo::f11SIMD(  /// input 1st stage of singlet search
                        int istal, int istam,  /// indexes of left and middle stations of a triplet
                        Tindex n1_V,           ///
                        
                        fvec *u_front, fvec *u_back,  fvec *zPos,
                        // output
                        
                        
                        L1TrackParV *T_1V,
                       
                        L1FieldRegionV *fld_1V
                        )
{
    L1Station &stal = vStations[istal];
    L1Station &stam = vStations[istam];
    fvec zstal = stal.z;
    fvec zstam = stam.z;
    
    L1FieldRegion fld0;      // by  left   hit, target and "center" (station in-between). Used here for extrapolation to target and to middle hit
    L1FieldValue centB, l_B _fvecalignment; // field for singlet creation
    L1FieldValue m_B _fvecalignment; // field for the next extrapolations
    
    L1FieldRegion fld1; // by  middle hit, left hit and "center" . Will be used for extrapolation to right hit
    
    L1TrackPar T;

    
    for( Tindex i1_V=0; i1_V<n1_V; ++i1_V ){
        
      //  L1TrackParR T(T_1V,(i1_V*fvecLen));
        
        // get the magnetic field along the track.
        // (suppose track is straight line with origin in the target)
        fvec &u = (u_front[i1_V]);
        fvec &v = (u_back [i1_V]);
        fvec xl, yl; // left(1-st) hit coor
        fvec zl = (zPos   [i1_V]);
        
        fvec dzli = 1.f/(zl-targZ);
        StripsToCoor(u, v, xl, yl, stal); //hit coordinates from strips
        
        fvec tx = (xl - targX)*dzli;
        fvec ty = (yl - targY)*dzli;
        
        // estimate field for singlet creation
        int istac = istal/2; // "center" station
        //     if (istal >= NMvdStations) // to make it in the same way for both with and w\o mvd
        //       istac = (istal-NMvdStations)/2+NMvdStations;
        L1Station &stac = vStations[istac]; //estimate target field
        fvec zstac = stac.z;
        stac.fieldSlice.GetFieldValue( targX + tx*(zstac - targZ), targY + ty*(zstac - targZ), centB );
        stal.fieldSlice.GetFieldValue( targX + tx*(zstal - targZ), targY + ty*(zstal - targZ), l_B );
        if( istac != istal  ){
            fld0.Set( l_B, stal.z, centB, stac.z, targB, targZ );
        }
        else{
            fld0.Set( l_B, zstal, targB, targZ );
        }
        // estimate field for the next extrapolations
        stam.fieldSlice.GetFieldValue( targX + tx*(zstam - targZ), targY + ty*(zstam - targZ), m_B );
#define USE_3HITS
#ifndef USE_3HITS
        if( istac != istal ){
            fld1.Set( m_B, zstam, l_B, zstal, centB, zstac  );
        }
        else{
            fld1.Set( m_B, zstam, l_B, zstal, targB, targZ  );
        }
#else // if USE_3HITS  // the best now
        L1FieldValue r_B _fvecalignment;
        L1Station &star = vStations[istam+1];
        fvec zstar = star.z;
        star.fieldSlice.GetFieldValue( targX + tx*(zstar - targZ), targY + ty*(zstar - targZ), r_B );
        fld1.Set( r_B, zstar, m_B, zstam, l_B, zstal);
#endif // USE_3HITS
        
        T.chi2 = 0.;
        T.NDF = 2.;
        if ( (isec == kAllSecIter) || (isec == kAllSecJumpIter) ) T.NDF = 0;
        T.tx = tx;
        T.ty = ty;
        T.qp = 0.;
        T.C20 = T.C21 = 0;
        T.C30 = T.C31 = T.C32 = 0;
        T.C40 = T.C41 = T.C42 = T.C43 = 0;
        T.C22 = T.C33 = MaxSlope*MaxSlope/9.f; T.C44 = MaxInvMom/3.f*MaxInvMom/3.f;
        
        // #define BEGIN_FROM_TARGET
#ifndef BEGIN_FROM_TARGET // the best now
        
        T.x  = xl;
        T.y  = yl;
        T.z  = zl;
        T.C00 = stal.XYInfo.C00;
        T.C10 = stal.XYInfo.C10;
        T.C11 = stal.XYInfo.C11;
        
        //add the target
        {
            fvec eX, eY, J04, J14;
            fvec dz;
            dz = targZ - zl;
            L1ExtrapolateJXY0(T.tx, T.ty, dz, fld0, eX, eY, J04, J14 );
            fvec J[6];
            J[0]= dz; J[1] = 0; J[2]= J04;
            J[3] = 0; J[4]= dz; J[5]= J14;
            L1FilterVtx( T, targX, targY, TargetXYInfo, eX, eY, J );
        }
#else  // TODO: doesn't work. Why?
        
        T.x  = targX;
        T.y  = targY;
        
        T.z  = targZ;
        T.C00 = TargetXYInfo.C00;
        T.C10 = TargetXYInfo.C10;
        T.C11 = TargetXYInfo.C11;
        
        
        // extrapolate to left hit
        L1Extrapolate0( T, zl, fld0 );
        for (int ista = 0; ista <= istal-1; ++ista){
            L1AddMaterial( T, vStations[ista].materialInfo, MaxInvMom );
            if (ista == NMvdStations - 1) L1AddPipeMaterial( T, MaxInvMom );
        }
        // add left hit
        L1Filter( T, stal.frontInfo, u );
        L1Filter( T, stal.backInfo,  v );
#endif
        
        L1AddMaterial( T, stal.materialInfo, MaxInvMom );
        if ( (istam >= NMvdStations) && (istal <= NMvdStations - 1) )  L1AddPipeMaterial( T, MaxInvMom );
        
        L1Extrapolate0( T, zstam, fld0 ); // TODO: fld1 doesn't work!
        //     L1Extrapolate( T, zstam, T.qp, fld1 );
        
        T_1V->Set(i1_V*fvecLen, T);
        fld_1V->Set(i1_V*fvecLen, fld1);
        
    }// i1_V
}

          /// Get the field approximation. Add the target to parameters estimation. Propagaete to the middle station of a triplet.
inline void L1Algo::f11(  /// input 1st stage of singlet search
                int istal, int istam,  /// indexes of left and middle stations of a triplet 
                Tindex n1_V,           /// 

                fvec *u_front, fvec *u_back,  fvec *zPos,
                  // output
                
                L1TrackPar *T_1,
                
                L1FieldRegion *fld_1
               
               )
{
  L1Station &stal = vStations[istal];
  L1Station &stam = vStations[istam];
  fvec zstal = stal.z;
  fvec zstam = stam.z;
  
  L1FieldRegion fld0;      // by  left   hit, target and "center" (station in-between). Used here for extrapolation to target and to middle hit
  L1FieldValue centB, l_B _fvecalignment; // field for singlet creation
  L1FieldValue m_B _fvecalignment; // field for the next extrapolations
  
  for( Tindex i1_V=0; i1_V<n1_V; ++i1_V ){
    L1TrackPar &T = T_1[i1_V];

    L1FieldRegion &fld1 = fld_1[i1_V]; // by  middle hit, left hit and "center" . Will be used for extrapolation to right hit
 
      // get the magnetic field along the track.
      // (suppose track is straight line with origin in the target)
    fvec &u = u_front[i1_V];
    fvec &v = u_back [i1_V];
    fvec xl, yl; // left(1-st) hit coor
    fvec zl = zPos   [i1_V];

    fvec dzli = 1.f/(zl-targZ);
    StripsToCoor(u, v, xl, yl, stal); //hit coordinates from strips

    fvec tx = (xl - targX)*dzli;
    fvec ty = (yl - targY)*dzli;

            // estimate field for singlet creation
    int istac = istal/2; // "center" station
//     if (istal >= NMvdStations) // to make it in the same way for both with and w\o mvd
//       istac = (istal-NMvdStations)/2+NMvdStations;
    L1Station &stac = vStations[istac]; //estimate target field
    fvec zstac = stac.z;
    stac.fieldSlice.GetFieldValue( targX + tx*(zstac - targZ), targY + ty*(zstac - targZ), centB );
    stal.fieldSlice.GetFieldValue( targX + tx*(zstal - targZ), targY + ty*(zstal - targZ), l_B );
    if( istac != istal  ){
      fld0.Set( l_B, stal.z, centB, stac.z, targB, targZ );
    }
    else{
      fld0.Set( l_B, zstal, targB, targZ );
    }
          // estimate field for the next extrapolations
    stam.fieldSlice.GetFieldValue( targX + tx*(zstam - targZ), targY + ty*(zstam - targZ), m_B );
#define USE_3HITS
#ifndef USE_3HITS
    if( istac != istal ){
      fld1.Set( m_B, zstam, l_B, zstal, centB, zstac  );
    }
    else{
      fld1.Set( m_B, zstam, l_B, zstal, targB, targZ  );
    }
#else // if USE_3HITS  // the best now
    L1FieldValue r_B _fvecalignment;
    L1Station &star = vStations[istam+1];
    fvec zstar = star.z;
    star.fieldSlice.GetFieldValue( targX + tx*(zstar - targZ), targY + ty*(zstar - targZ), r_B );
    fld1.Set( r_B, zstar, m_B, zstam, l_B, zstal);
#endif // USE_3HITS

    T.chi2 = 0.;
    T.NDF = 2.;
    if ( (isec == kAllSecIter) || (isec == kAllSecJumpIter) ) T.NDF = 0;
    T.tx = tx;
    T.ty = ty;
    T.qp = 0.;
    T.C20 = T.C21 = 0;
    T.C30 = T.C31 = T.C32 = 0;
    T.C40 = T.C41 = T.C42 = T.C43 = 0;
    T.C22 = T.C33 = MaxSlope*MaxSlope/9.f; T.C44 = MaxInvMom/3.f*MaxInvMom/3.f;
    
// #define BEGIN_FROM_TARGET
#ifndef BEGIN_FROM_TARGET // the best now

    T.x  = xl;
    T.y  = yl;
    T.z  = zl;
    T.C00 = stal.XYInfo.C00;
    T.C10 = stal.XYInfo.C10;
    T.C11 = stal.XYInfo.C11;
    
    //add the target
    {
      fvec eX, eY, J04, J14;
      fvec dz;
      dz = targZ - zl;
      L1ExtrapolateJXY0(T.tx, T.ty, dz, fld0, eX, eY, J04, J14 );
      fvec J[6];
      J[0]= dz; J[1] = 0; J[2]= J04;
      J[3] = 0; J[4]= dz; J[5]= J14;
      L1FilterVtx( T, targX, targY, TargetXYInfo, eX, eY, J );
    }
#else  // TODO: doesn't work. Why? 

    T.x  = targX;
    T.y  = targY;

    T.z  = targZ;
    T.C00 = TargetXYInfo.C00;
    T.C10 = TargetXYInfo.C10;
    T.C11 = TargetXYInfo.C11;

    
      // extrapolate to left hit
    L1Extrapolate0( T, zl, fld0 );
    for (int ista = 0; ista <= istal-1; ++ista){
      L1AddMaterial( T, vStations[ista].materialInfo, MaxInvMom );
      if (ista == NMvdStations - 1) L1AddPipeMaterial( T, MaxInvMom );
    }
      // add left hit
    L1Filter( T, stal.frontInfo, u );
    L1Filter( T, stal.backInfo,  v );
#endif

    L1AddMaterial( T, stal.materialInfo, MaxInvMom );
    if ( (istam >= NMvdStations) && (istal <= NMvdStations - 1) )  L1AddPipeMaterial( T, MaxInvMom );

    L1Extrapolate0( T, zstam, fld0 ); // TODO: fld1 doesn't work!
//     L1Extrapolate( T, zstam, T.qp, fld1 );

   
    
  }// i1_V
}


/// Find the doublets. Reformat data in the portion of doublets.
void L1Algo::f20SIMD(  // input
                 Tindex n1, L1Station &stam,
                 
                
                 
                 L1TrackParV *T_1V,
               L1TrackParV *T_1VBuf,
                      L1TrackParV *T_1VBuf2,
                 nsL1::vector <int>::TSimd &MHits,
                 nsL1::vector <int>::TSimd &LHits,
                 
                 //                 THitI* hitsl_1,
                 
                 // output
                 Tindex &n2,
                     int indexHit,
                     nsL1::vector <int>::TSimd &LHits3
)
{
    n2 = 0; // number of doublet
    
    nsL1::vector <int>::TSimd MiddleHitsIndexes;
    nsL1::vector <int>::TSimd FirstHitsIndexes;
    
    MiddleHitsIndexes.reserve(4000);
    FirstHitsIndexes.reserve(4000);


   Tindex n1_V = (n1+fvecLen-1)/fvecLen;
   
    
   
    for (Tindex i1 = 0; i1 < n1_V-1; i1++)
    {
        Tindex index =i1*fvecLen;
        
        const fvec& X = reinterpret_cast<const fvec&>(T_1V->fP[0][index]);
        const fvec& Y = reinterpret_cast<const fvec&>(T_1V->fP[1][index]);
        const fvec& Z = reinterpret_cast<const fvec&>(T_1V->fP[5][index]);
        const fvec& Tx = reinterpret_cast<const fvec&>(T_1V->fP[2][index]);
        const fvec& Ty = reinterpret_cast<const fvec&>(T_1V->fP[3][index]);
        const fvec& Chi2 = reinterpret_cast<const fvec&>(T_1V->chi2[index]);
        const fvec& C00 = reinterpret_cast<const fvec&>(T_1V->fC[0][index]);
        const fvec& C11 = reinterpret_cast<const fvec&>(T_1V->fC[2][index]);
        
        
        fvec Pick_m22 = (DOUBLET_CHI2_CUT - Chi2);
        
        const fvec iz = 1/Z;
        fvec x = X*iz;
        fvec y = Y*iz;
        fvec dx = (sqrt(Pick_m22*(C00 + stam.XYInfo.C00))+MaxDZ*abs(Tx))*iz;
        fvec dy = (sqrt(Pick_m22*(C11 + stam.XYInfo.C11))+MaxDZ*abs(Ty))*iz;
        
       
        
        L1HitAreaV areaV( vGrid[ &stam - vStations ], x, y, dx, dy);
        
        
        areaV.HitsIndexesMiddle(MiddleHitsIndexes, FirstHitsIndexes, index);
        
        
    }
       Tindex index =(n1_V-1)*fvecLen;

        const fvec& X = reinterpret_cast<const fvec&>(T_1V->fP[0][index]);
        const fvec& Y = reinterpret_cast<const fvec&>(T_1V->fP[1][index]);
        const fvec& Z = reinterpret_cast<const fvec&>(T_1V->fP[5][index]);
        const fvec& Tx = reinterpret_cast<const fvec&>(T_1V->fP[2][index]);
        const fvec& Ty = reinterpret_cast<const fvec&>(T_1V->fP[3][index]);
        const fvec& Chi2 = reinterpret_cast<const fvec&>(T_1V->chi2[index]);
        const fvec& C00 = reinterpret_cast<const fvec&>(T_1V->fC[0][index]);
        const fvec& C11 = reinterpret_cast<const fvec&>(T_1V->fC[2][index]);
        
        
        fvec Pick_m22 = (DOUBLET_CHI2_CUT - Chi2);
        
        const fvec iz = 1/Z;
        fvec x = X*iz;
        fvec y = Y*iz;
        fvec dx = (sqrt(Pick_m22*(C00 + stam.XYInfo.C00))+MaxDZ*abs(Tx))*iz;
        fvec dy = (sqrt(Pick_m22*(C11 + stam.XYInfo.C11))+MaxDZ*abs(Ty))*iz;
        
       
        
        L1HitAreaV areaV( vGrid[ &stam - vStations ], x, y, dx, dy);
        
        
        areaV.HitsIndexesMiddle(MiddleHitsIndexes, FirstHitsIndexes, index, n1);
    
    
        n1_V = (MiddleHitsIndexes.size()+fvecLen-1)/fvecLen;   



    for(int iPair=0; iPair<n1_V; iPair++)
    {


        Tindex index = iPair*fvecLen;
        int_v hitM_V = reinterpret_cast<int_v&>(MiddleHitsIndexes[index]); // vector middle hits indexes in the initial L1HitPonts array
        int_m IsNiceDoublet (true);

        
        
        
        if((index + fvecLen-1)>= MiddleHitsIndexes.size())
            IsNiceDoublet &= int_m( int_v(Vc::IndexesFromZero) < int(MiddleHitsIndexes.size()%fvecLen) );
        
        
        
        hitM_V(!(IsNiceDoublet))=0;
        
        
        const fvec zm(vStsHitPointsUnusedV[ &stam - vStations ].Z(), hitM_V);
        const fvec ym(vStsHitPointsUnusedV[ &stam - vStations ].Y(), hitM_V);

      
        
        int_v hitL_V = reinterpret_cast<int_v&>(FirstHitsIndexes[index]);
        hitL_V(!(IsNiceDoublet))=0;

        T_1V->CopyTrackPar((*T_1VBuf), hitL_V, index);

        
        L1TrackParR T (T_1VBuf, index);
        
        fvec Pick_m22 = (fvec(DOUBLET_CHI2_CUT) - T.chi2);
        
        
        
        fvec y, C11;
        L1ExtrapolateYC11Line( T, zm, y, C11 );
        
        
        const fvec dy_est2 = ( Pick_m22 * abs(C11 + stam.XYInfo.C11) );
        
        const fvec dy = ym - y;
        const fvec dy2 = dy*dy;
        
        
        IsNiceDoublet &=(!(dy2 > dy_est2 && dy < fvec(Vc::Zero)));
        if ((IsNiceDoublet.isEmpty())) continue;
        
        const fvec xm(vStsHitPointsUnusedV[ &stam - vStations ].X(), hitM_V);
        fvec x, C00;
        L1ExtrapolateXC00Line( T, zm, x, C00 );
        const fvec dx_est2 = (Pick_m22 * abs(C00 + stam.XYInfo.C00));
        const fvec dx = xm - x;
        
             IsNiceDoublet &= (dx*dx <= dx_est2);
         if ((IsNiceDoublet.isEmpty())) continue;
        
        // check chi2
        fvec C10;
        L1ExtrapolateC10Line( T, zm, C10 );
        fvec chi2 = T.chi2;
        const fvec um(vStsHitPointsUnusedV[ &stam - vStations ].U(), hitM_V);
        L1FilterChi2XYC00C10C11( stam.frontInfo, x, y, C00, C10, C11, chi2, um );
             IsNiceDoublet &= (chi2 <= DOUBLET_CHI2_CUT);
          if ((IsNiceDoublet.isEmpty())) continue;
        const fvec vm(vStsHitPointsUnusedV[ &stam - vStations ].V(), hitM_V);
        L1FilterChi2           ( stam.backInfo,  x, y, C00, C10, C11, chi2, vm);
             IsNiceDoublet &= (chi2 <= TRIPLET_CHI2_CUT*(T.NDF-fvec(1)) );
        if ((IsNiceDoublet.isEmpty())) continue;


        
        
        for (int iV=0; iV<fvecLen; iV++)
        {
            if (!(IsNiceDoublet[iV])) continue;
          //  LHits.push_back(hitL_V[iV]);
            
            LHits.push_back(index+iV);
            
            MHits.push_back(hitM_V[iV]);
          
            LHits3.push_back(hitL_V[iV]);
           
            // LHits3.push_back(hitL_V[iV]+index);
            n2++;
        }
        }
   
   }



          /// Find the doublets. Reformat data in the portion of doublets.
 void L1Algo::f20(  // input
                Tindex n1, L1Station &stam,

                L1HitPoint *vStsHits_m,
                L1TrackPar *T_1, 
                 
               
//                 THitI* hitsl_1,

                  // output
               Tindex &n2,
                vector<THitI> &i1_2,
#ifdef DOUB_PERFORMANCE
                vector<THitI> &hitsl_2,
#endif // DOUB_PERFORMANCE
                vector<THitI> &hitsm_2
//                 ,vector<bool> &lmDuplets
                )
{
 n2 = 0; // number of doublet


  for (Tindex i1 = 0; i1 < n1; ++i1){ // for each singlet


    const Tindex i1_V = i1/fvecLen;
    const Tindex i1_4 = i1%fvecLen;
    const L1TrackPar& T1 = T_1[i1_V];

//     const int n2Saved = n2;

      //  const THitI nl = vStsHits_l[hitsl_1[i1]].N();
    fvec Pick_m22 = (DOUBLET_CHI2_CUT - T1.chi2); // if make it bigger the found hits will be rejected later because of the chi2 cut.
  // Pick_m22 is not used, search for mean squared, 2nd version
  
    
      // -- collect possible doublets --


      
    const fscal iz = 1/T1.z[i1_4]; 
    L1HitArea area( vGrid[ &stam - vStations ], T1.x[i1_4]*iz, T1.y[i1_4]*iz, (sqrt(Pick_m22*(T1.C00 + stam.XYInfo.C00))+MaxDZ*abs(T1.tx))[i1_4]*iz, (sqrt(Pick_m22*(T1.C11 + stam.XYInfo.C11))+MaxDZ*abs(T1.ty))[i1_4]*iz ); 
    THitI imh = 0;

    while( area.GetNext( imh ) ) { 
     
      const L1HitPoint &hitm = vStsHits_m[imh];


        // check y-boundaries

        // - check whether hit belong to the window ( track position +\- errors ) -
      const fscal &zm = hitm.Z();

        // check lower boundary
      fvec y, C11;
      L1ExtrapolateYC11Line( T1, zm, y, C11 );

      const fscal dy_est2 = ( Pick_m22 * abs(C11 + stam.XYInfo.C11) )[i1_4];

      const fscal dy = hitm.Y() - y[i1_4];
      const fscal dy2 = dy*dy;
      if ( dy2 > dy_est2 && dy < 0 ) continue;

        // check upper boundary
    
        // const Tindex nm = hitm.N();
        // cout<<nl<<" nm"<<nm<<" nl"<<endl;
        //  if ( ( nl != nm ) ) continue;

        // check x-boundaries
      fvec x, C00;
      L1ExtrapolateXC00Line( T1, zm, x, C00 );
      const fscal dx_est2 = (Pick_m22 * abs(C00 + stam.XYInfo.C00))[i1_4];
      const fscal dx = hitm.X() - x[i1_4];
      if ( dx*dx > dx_est2 ) continue;

        // check chi2
      fvec C10;
      L1ExtrapolateC10Line( T1, zm, C10 );
      fvec chi2 = T1.chi2;
      L1FilterChi2XYC00C10C11( stam.frontInfo, x, y, C00, C10, C11, chi2, hitm.U() );
      if ( chi2[i1_4] > DOUBLET_CHI2_CUT ) continue;
      L1FilterChi2           ( stam.backInfo,  x, y, C00, C10, C11, chi2, hitm.V() );
      if ( chi2[i1_4] > TRIPLET_CHI2_CUT*(T1.NDF[i1_4]-1) ) continue; // chi2_doublet < chi2_triplet < CHI2_CUT

      i1_2.push_back(i1);

#ifdef DOUB_PERFORMANCE
//       hitsl_2.push_back(hitsl_1[i1]);
#endif // DOUB_PERFORMANCE
      hitsm_2.push_back(imh);

      n2++;
    }
   
 
//       if (n2Saved >= n2) continue;
//       if ( lmDuplets[hitsl_1[i1]]==1) continue;
//     
//     lmDuplets[hitsl_1[i1]]=1;
     
  }  // for i1

   
// cout<<n2<<" n2 2 start<<"<<endl;

}


          /// Add the middle hits to parameters estimation. Propagate to right station.
          /// Find the triplets(right hit). Reformat data in the portion of triplets.
inline void L1Algo::f30(  // input
                L1HitPoint *vStsHits_r, L1Station &stam, L1Station &star,

                int istam, int istar,
                L1HitPoint *vStsHits_m,
                L1TrackPar *T_1, L1FieldRegion *fld_1,
                THitI *hitsl_1,
               

                Tindex n2,
                vector<THitI> &hitsm_2,
                vector<THitI> &i1_2,

//                 const vector<bool> &mrDuplets,

                  // output
                Tindex &n3,
                nsL1::vector<L1TrackPar>::TSimd &T_3,
                vector<THitI> &hitsl_3,  vector<THitI> &hitsm_3,  vector<THitI> &hitsr_3,
                nsL1::vector<fvec>::TSimd &u_front_3, nsL1::vector<fvec>::TSimd &u_back_3, nsL1::vector<fvec>::TSimd &z_Pos_3
     

                )
{
    // RHits.resize(0);
  THitI hitsl_2[fvecLen];
  THitI hitsm_2_tmp[fvecLen];
  fvec fvec_0;
  L1TrackPar L1TrackPar_0;

 Tindex n3_V = 0, n3_4 = 0;
    
  T_3.push_back(L1TrackPar_0);
  u_front_3.push_back(fvec_0);
  u_back_3.push_back(fvec_0);
  z_Pos_3.push_back(fvec_0);
       L1TrackPar T2;
      L1FieldRegion f2;
      // pack the data
      fvec u_front_2;
      fvec u_back_2;
      fvec zPos_2; 

    
  //  if (grid1>0) return;
    
    

    // ---- Add the middle hits to parameters estimation. Propagate to right station. ----
  if (istar < NStations){


    for (Tindex i2 = 0; i2 < n2;) {



      Tindex n2_4 = 0;
      for (; n2_4 < fvecLen && i2 < n2; ++n2_4, ++i2){

//         if (!mrDuplets[hitsm_2[i2]]) {
//           n2_4--;
//           continue;
//         }

        const Tindex i1 = i1_2[i2];
          
         
        const Tindex i1_V = i1/fvecLen;
        const Tindex i1_4 = i1%fvecLen;

        const L1TrackPar &T1 = T_1[i1_V];
          
        //  cout<<i1_2[i2]<<" hitL_V1"<<endl;
        const L1FieldRegion &f1 = fld_1[i1_V];
        T2.SetOneEntry(n2_4, T1, i1_4);
        f2.SetOneEntry(n2_4, f1, i1_4);
         
        const Tindex imh = hitsm_2[i2];
        const L1HitPoint &hitm = vStsHits_m[imh];
        u_front_2[n2_4] = hitm.U();
        u_back_2 [n2_4] = hitm.V();
        zPos_2   [n2_4] = hitm.Z();
          
         


        hitsl_2[n2_4] = hitsl_1[i1];
        hitsm_2_tmp[n2_4] = hitsm_2[i2];

//cout<<f2.cx1[n2_4]<<" "<<hitsl_2[n2_4]<<" "<<hitsm_2_tmp[n2_4]<<" "<<n2<<" scalar"<<endl;
      }  // n2_4


        // add middle hit
      L1ExtrapolateLine( T2, zPos_2 );
      L1Filter( T2, stam.frontInfo, u_front_2 );
      L1Filter( T2, stam.backInfo,  u_back_2 );
       
      L1AddMaterial( T2, stam.materialInfo, T2.qp );
        

       
      if ( (istar >= NMvdStations) && (istam <= NMvdStations - 1) ) L1AddPipeMaterial( T2, T2.qp );


        // extrapolate to the right hit station
      L1Extrapolate( T2, star.z, T2.qp, f2 );
        
      //  cout<<T2.chi2<<" T2.chi2 i"<<endl;

        // ---- Find the triplets(right hit). Reformat data in the portion of triplets. ----

      for (Tindex i2_4 = 0; i2_4 < n2_4; ++i2_4){
    
        
          //     THitI nm = vStsHits_m[hitsm_2[i2]].n;
        
        if ( T2.C00[i2_4] < 0.f ||  T2.C11[i2_4] < 0.f ||  T2.C22[i2_4] < 0.f ||  T2.C33[i2_4] < 0.f ||  T2.C44[i2_4] < 0.f ) continue;
  
  fvec Pick_r22 = (TRIPLET_CHI2_CUT - T2.chi2);

          // find first possible hit
      

         
      
          const fscal iz = 1/T2.z[i2_4];
          
          
         
         
    L1HitArea area( vGrid[ &star - vStations ], T2.x[i2_4]*iz, T2.y[i2_4]*iz, (sqrt(Pick_r22*(T2.C00 + stam.XYInfo.C00))+MaxDZ*abs(T2.tx))[i2_4]*iz, (sqrt(Pick_r22*(T2.C11 + stam.XYInfo.C11))+MaxDZ*abs(T2.ty))[i2_4]*iz ); 
    THitI irh = 0; 
    while( area.GetNext( irh ) ) {
        
        
     //   cout<<i2<<" i2 hit"<<irh<<endl;
                
          const L1HitPoint &hitr = vStsHits_r[irh];

          const fscal zr = hitr.Z();
          const fscal yr = hitr.Y();
        
     //    cout<<zr<<" zr"<<endl;
        
            // - check whether hit belong to the window ( track position +\- errors ) -
            // check lower boundary
          fvec y, C11;
          L1ExtrapolateYC11Line( T2, zr, y, C11 );
          const fscal dy_est2 = (Pick_r22*(abs(C11 + star.XYInfo.C11)))[i2_4]; // TODO for FastPrim dx < dy - other sort is optimal. But not for doublets
          const fscal dy = yr - y[i2_4];
          const fscal dy2 = dy*dy;
        
     //   cout<<dy_est2<<" dy_est2"<<endl;
        
        
        //cout<<dy<<" dy"<<endl;
        
       
       

        
          if ( dy2 > dy_est2 && dy < 0 ) continue; // if (yr < y_minus_new[i2_4]) continue;
        
        
        
        
            // check upper boundary
          if ( dy2 > dy_est2 ) continue; // if (yr > y_plus_new [i2_4] ) continue;
        
       
        
            // check x-boundaries
          fvec x, C00;
          L1ExtrapolateXC00Line( T2, zr, x, C00 );
          const fscal dx_est2 = (Pick_r22*(abs(C00 + star.XYInfo.C00)))[i2_4];
          const fscal dx = hitr.X() - x[i2_4];
          if ( dx*dx > dx_est2 ) continue;

            // check chi2  // not effective
          fvec C10;
          L1ExtrapolateC10Line( T2, zr, C10 );
          fvec chi2 = T2.chi2;
          L1FilterChi2XYC00C10C11( star.frontInfo, x, y, C00, C10, C11, chi2, hitr.U() );
          L1FilterChi2           ( star.backInfo,  x, y, C00, C10, C11, chi2, hitr.V() );
          if ( chi2[i2_4] > TRIPLET_CHI2_CUT || C00[i2_4] < 0.f || C11[i2_4] < 0.f ) continue; // chi2_triplet < CHI2_CUT
        
        
        
            // pack triplet
          L1TrackPar &T3 = T_3[n3_V];

          hitsl_3.push_back(hitsl_2[i2_4]);
          hitsm_3.push_back(hitsm_2_tmp[i2_4]);
          hitsr_3.push_back(irh);
        
     //   cout<<T2.chi2[i2_4]<<" ch 1"<<endl;

          T3.SetOneEntry(n3_4, T2, i2_4);

          u_front_3[n3_V][n3_4] = hitr.U();
          u_back_3 [n3_V][n3_4] = hitr.V();
          z_Pos_3  [n3_V][n3_4] = zr;

          n3++;               
          n3_V = n3/fvecLen;  
          n3_4 = n3%fvecLen;

          if (!n3_4){
            T_3.push_back(L1TrackPar_0);

            u_front_3.push_back(fvec_0);
            u_back_3.push_back(fvec_0);
            z_Pos_3.push_back(fvec_0);
            
          }
        }
          
          
      } // i2_4
    }   // i2_V

 // cout<<n3<<" n3 scal"<<endl;
  }  // if istar
    
    

  
    
  
}


  /// Add the right hits to parameters estimation.
inline void L1Algo::f31(  // input
                Tindex n3_V,
                L1Station &star,
                nsL1::vector<fvec>::TSimd &u_front, nsL1::vector<fvec>::TSimd &u_back, nsL1::vector<fvec>::TSimd &z_Pos, 
                  // output
                //                L1TrackPar *T_3
                nsL1::vector<L1TrackPar>::TSimd &T_3
               )
{
  for( Tindex i3_V=0; i3_V<n3_V; ++i3_V){
    L1ExtrapolateLine( T_3[i3_V], z_Pos[i3_V] );
    L1Filter( T_3[i3_V], star.frontInfo, u_front[i3_V] );    // 2.1/100 sec
    L1Filter( T_3[i3_V], star.backInfo,  u_back [i3_V] );   // 2.0/100 sec
//     cout<<sqrt(T_3[i3_V].C00)<<endl;
//     cout<<sqrt(T_3[i3_V].C10)<<endl;
//     cout<<sqrt(T_3[i3_V].C11)<<endl;
      

  };
  
}


          /// Refit Triplets.
inline void L1Algo::f32( // input // TODO not updated after gaps introduction
    Tindex n3, int istal,
                nsL1::vector<L1TrackPar>::TSimd &T_3,
                vector<THitI> &hitsl_3,  vector<THitI> &hitsm_3,  vector<THitI> &hitsr_3,
                
                        
                        int nIterations
                       )
{
  const int NHits = 3; // triplets

    // prepare data
  int ista[NHits] = {
    istal,
    istal + 1,
    istal + 2
  };
  
  L1Station sta[3];
  for (int is = 0; is < NHits; ++is){
    sta[is] = vStations[ista[is]];
  };

  for( int i3=0; i3<n3; ++i3){
    int i3_V = i3/fvecLen;
    int i3_4 = i3%fvecLen;

    L1TrackPar &T3 = T_3[i3_V];
    fscal qp0 = T3.qp[i3_4];

      // prepare data
    THitI ihit[NHits] = {
      (*RealIHitP)[hitsl_3[i3] + StsHitsUnusedStartIndex[ista[0]]],
      (*RealIHitP)[hitsm_3[i3] + StsHitsUnusedStartIndex[ista[1]]],
      (*RealIHitP)[hitsr_3[i3] + StsHitsUnusedStartIndex[ista[2]]]
    };

    fscal u[NHits], v[NHits], x[NHits], y[NHits], z[NHits];
    for (int ih = 0; ih < NHits; ++ih){
      const L1StsHit &hit = (*vStsHits)[ihit[ih]];
      u[ih] = (*vStsStrips)[hit.f];
      v[ih] = (*vStsStripsB)[hit.b];
      StripsToCoor(u[ih], v[ih], x[ih], y[ih], sta[ih]);
      z[ih] = (*vStsZPos)[hit.iz];
    };

      // initialize parameters
    L1TrackPar T;
    
    const fvec vINF = .1f;
    T.x  = x[0];
    T.y  = y[0];

    fvec dz01 = 1.f/(z[1]-z[0]);
    T.tx = (x[1]-x[0])*dz01;
    T.ty = (y[1]-y[0])*dz01;

    T.qp = qp0;
    T.z  = z[0];
    T.chi2 = 0.f;
    T.NDF = 2.f;
    T.C00 = sta[0].XYInfo.C00;
    T.C10 = sta[0].XYInfo.C10;
    T.C11 = sta[0].XYInfo.C11;

    T.C20 = T.C21 = 0;
    T.C30 = T.C31 = T.C32 = 0;
    T.C40 = T.C41 = T.C42 = T.C43 = 0;
    T.C22 = T.C33 = vINF;
    T.C44 = 1.;

      // find field along the track
    L1FieldValue B[3] _fvecalignment;
    L1FieldRegion fld _fvecalignment;

    fvec tx[3] = {
      (x[1]-x[0])/(z[1]-z[0]),
      (x[2]-x[0])/(z[2]-z[0]),
      (x[2]-x[1])/(z[2]-z[1])
    };
    fvec ty[3] = {
      (y[1]-y[0])/(z[1]-z[0]),
      (y[2]-y[0])/(z[2]-z[0]),
      (y[2]-y[1])/(z[2]-z[1])
    };
    for (int ih = 0; ih < NHits; ++ih){
      fvec dz = (sta[ih].z - z[ih]);
      sta[ih].fieldSlice.GetFieldValue( x[ih] + tx[ih]*dz, y[ih] + ty[ih]*dz, B[ih] );
    };
    fld.Set( B[0], sta[0].z, B[1], sta[1].z, B[2], sta[2].z );

      // fit
    for( int ih = 1; ih < NHits; ++ih){
      L1Extrapolate( T, z[ih], T.qp, fld );
      L1AddMaterial( T, sta[ih].materialInfo, T.qp );
      if (ista[ih] == NMvdStations - 1) L1AddPipeMaterial( T, T.qp );
      
      L1Filter( T, sta[ih].frontInfo, u[ih] );
      L1Filter( T, sta[ih].backInfo,  v[ih] );
    }

      // repeat several times in order to improve precision
    for (int iiter = 0; iiter < nIterations; ++iiter){
        // fit backward
      // keep tx,ty,q/p
      int ih = NHits-1;
      T.x  = x[ih];
      T.y  = y[ih];
      T.z  = z[ih];
      T.chi2 = 0.;
      T.NDF = 2.;
      T.C00 = sta[ih].XYInfo.C00;
      T.C10 = sta[ih].XYInfo.C10;
      T.C11 = sta[ih].XYInfo.C11;

      T.C20 = T.C21 = 0;
      T.C30 = T.C31 = T.C32 = 0;
      T.C40 = T.C41 = T.C42 = T.C43 = 0;
      T.C22 = T.C33 = vINF;
      T.C44 = 1.;

//       L1Filter( T, sta[ih].frontInfo, u[ih] );
//       L1Filter( T, sta[ih].backInfo,  v[ih] );
      for( ih = NHits-2; ih >= 0; ih--){
        L1Extrapolate( T, z[ih], T.qp, fld );
        L1AddMaterial( T, sta[ih].materialInfo, T.qp );
        if (ista[ih] == NMvdStations) L1AddPipeMaterial( T, T.qp );
        
        L1Filter( T, sta[ih].frontInfo, u[ih] );
        L1Filter( T, sta[ih].backInfo,  v[ih] );
      }
        // fit forward
      ih = 0;
      T.x  = x[ih];
      T.y  = y[ih];
      T.z  = z[ih];
      T.chi2 = 0.;
      T.NDF = 2.;
      T.C00 = sta[ih].XYInfo.C00;
      T.C10 = sta[ih].XYInfo.C10;
      T.C11 = sta[ih].XYInfo.C11;

      T.C20 = T.C21 = 0;
      T.C30 = T.C31 = T.C32 = 0;
      T.C40 = T.C41 = T.C42 = T.C43 = 0;
      T.C22 = T.C33 = vINF;
      T.C44 = 1.;

//       L1Filter( T, sta[ih].frontInfo, u[ih] );
//       L1Filter( T, sta[ih].backInfo,  v[ih] );
      for( ih = 1; ih < NHits; ++ih){
        L1Extrapolate( T, z[ih], T.qp, fld );
        L1AddMaterial( T, sta[ih].materialInfo, T.qp );
        if (ista[ih] == NMvdStations + 1) L1AddPipeMaterial( T, T.qp );
        
        L1Filter( T, sta[ih].frontInfo, u[ih] );
        L1Filter( T, sta[ih].backInfo,  v[ih] );
      }
    } // for iiter

    T3.SetOneEntry(i3_4, T,i3_4);
  }//i3




} // f32


          /// Select triplets. Save them into vTriplets.
inline void L1Algo::f4(  // input
                Tindex n3, int istal, int istam, int istar,
                nsL1::vector<L1TrackPar>::TSimd &T_3,

                vector<THitI> &hitsl_3,  vector<THitI> &hitsm_3,  vector<THitI> &hitsr_3,
                // output
                Tindex &nstaltriplets,
                
                          Tindex ip_cur 
                      
               
               )
{
    
    
  for( Tindex i3=0; i3<n3; ++i3){
    
    
    Tindex i3_V = i3/fvecLen;
    Tindex i3_4 = i3%fvecLen;
    L1TrackPar &T3 = T_3[i3_V];
    
      // select
    fscal chi2 = T3.chi2[i3_4];
     if ( !finite(chi2) || chi2 < 0 || chi2 > TRIPLET_CHI2_CUT )  continue;
    // if ( isec == kAllSecIter && T3.qp[i3_4] +sqrt( T3.C44[i3_4] ) < 1./10. ) continue; // why does it cut so much of ExtraSec ?
    
      // prepare data
//     fscal MaxInvMomS = MaxInvMom[0];
     // fscal qp = MaxInvMomS + T3.qp[i3_4];
    fscal qp = T3.qp[i3_4];
//     if( qp < 0 )            qp = 0;          
//     if( qp > MaxInvMomS*2 )  qp = MaxInvMomS*2;
    fscal Cqp = 5.*sqrt(abs(T3.C44[i3_4]));  


    
      
    THitI ihitl = hitsl_3[i3] + StsHitsUnusedStartIndex[istal];
    THitI ihitm = hitsm_3[i3] + StsHitsUnusedStartIndex[istam];
    THitI ihitr = hitsr_3[i3] + StsHitsUnusedStartIndex[istar];
      
      
      
    L1_ASSERT( ihitl < StsHitsUnusedStopIndex[istal], ihitl << " < " << StsHitsUnusedStopIndex[istal] );
    L1_ASSERT( ihitm < StsHitsUnusedStopIndex[istam], ihitm << " < " << StsHitsUnusedStopIndex[istam] );
    L1_ASSERT( ihitr < StsHitsUnusedStopIndex[istar], ihitr << " < " << StsHitsUnusedStopIndex[istar] );
    
     
    
    L1Triplet & tr=TripletsLocalReal[istal][ip_cur][TripNPortion[istal][ip_cur]];
      TripletsLocalRealAll[istal][ip_cur][TripNPortion[istal][ip_cur]]=&tr;

    if (istal!=FIRSTCASTATION) {
    //    L1Triplet * trip = new L1Triplet;
      //trip[omp_get_thread_num()];

    //  L1Triplet & tr=trip[omp_get_thread_num()];
//       trip->Set( ihitl, ihitm, ihitr,
//                  istal, istam, istar,
//                  0, qp, chi2);
                 tr.Set( ihitl, ihitm, ihitr,
                 istal, istam, istar,
                 0, qp, chi2);
    //trip.Cqp  nthrenthre = static_cast<Tindex>( Cqp );
   // trip->Cqp   = ( Cqp );
    tr.Cqp   = ( Cqp );
          //count statistics
  //  TripletsLocalReal[istal][ip_cur][TripNPortion[istal][ip_cur]]=(tr);
  //  if(TripNPortion[istal][ip_cur]>3000) cout<<" NOOOOOOOOOOOO!!!"<<endl;
   // vTriplets_part.push_back(trip);  // 0.67/100 sec
        
      
   
    TripNPortion[istal][ip_cur]++;
   //   if(ip_cur>1000) cout<<"Q NOOOOOOOOOOOO!!!"<<endl;
    
     // vTriplets_part.push_back(trip);  
      ++nstaltriplets; 
//     #pragma omp critical
//      vTripletsP[istal].push_back(trip);

      TripForHit[tr.GetLHit()].push_back(&TripletsLocalReal[istal][ip_cur][TripNPortion[istal][ip_cur]-1]);
    }

     tr.chi2Chain=tr.GetChi2();
    
    if (istal > (NStations - 4)) continue;
    
    if ( TripForHit[ihitm].size() == 0) continue;
  
    unsigned char level = 0;
    vector<L1Triplet *> neighCands; // save neighbour candidates
    neighCands.reserve(8); // ~average is 1-2 for central, up to 5
          
           

    for (unsigned int iN=0; iN < TripForHit[ihitm].size(); ++iN){
 
      L1Triplet* &curNeighbour = TripForHit[ihitm][iN];

      if (curNeighbour->GetMHit() != ihitr) continue; // neighbours should have 2 common hits
            
      L1Triplet* &tripn = curNeighbour;

      fscal qp2 = tripn->GetQp();
      fscal Cqp2 = tripn->Cqp;
      if ( abs(qp - qp2) > PickNeighbour * (Cqp + Cqp2) )  continue; // neighbours should have same qp
  
              // calculate level
      unsigned char jlevel = tripn->GetLevel();
      if ( level <= jlevel )
        level = jlevel + 1;
        
         //if  (( (fabs(qp - qp2))*(fabs(qp - qp2))/(Cqp*5.*Cqp*5.)) > TRACK_CHI2_CUT * level ) {qit++; continue; };
      
      neighCands.push_back(curNeighbour);
      

    }
   
 
    if (istal==FIRSTCASTATION) {
      if (level==0) continue;

     //  L1Triplet * trip = new L1Triplet;
      

//       trip->Set( ihitl, ihitm, ihitr,
//                  istal, istam, istar,
//                  0, qp, chi2);
                 tr.Set( ihitl, ihitm, ihitr,
                 istal, istam, istar,
                 0, qp, chi2);
    //trip.Cqp  nthrenthre = static_cast<Tindex>( Cqp );
   // trip->Cqp   = ( Cqp );
    tr.Cqp   = ( Cqp );
          //count statistics
  //  TripletsLocalReal[istal][ip_cur][TripNPortion[istal][ip_cur]]=(tr);
   //  if(TripNPortion[istal][ip_cur]>3000) cout<<" NOOOOOOOOOOOO!!!"<<endl;
  //  vTriplets_part.push_back(trip);
        
   
     TripNPortion[istal][ip_cur]++;
 //     if(ip_cur>1000) cout<<"Q NOOOOOOOOOOOO!!!"<<endl;
      ++nstaltriplets; 
//     #pragma omp critical
//      vTripletsP[istal].push_back(trip);

      TripForHit[tr.GetLHit()].push_back(&TripletsLocalReal[istal][ip_cur][TripNPortion[istal][ip_cur]-1]);
      
     
    }
   fscal BestChi2 = 1000000000;
    for (unsigned int in = 0; in < (neighCands.size()); ++in) {
      const Tindex nLevel = neighCands[in]->GetLevel();
 
      if (level == nLevel + 1) {
     //   (vTriplets_part[vTriplets_part.size()-1])->neighbours.push_back(neighCands[in]);
       TripletsLocalReal[istal][ip_cur][TripNPortion[istal][ip_cur]-1].neighbours.push_back(neighCands[in]);
      
      
            if ((neighCands[in]->chi2Chain)<BestChi2)
      {
        BestChi2=(neighCands[in]->chi2Chain);
       
      //  tr.neibourBest=neighCands[in];
      
      }
      
      
      }
       
    }
      
     
     std::sort(TripletsLocalReal[istal][ip_cur][TripNPortion[istal][ip_cur]-1].neighbours.begin(), TripletsLocalReal[istal][ip_cur][TripNPortion[istal][ip_cur]-1].neighbours.end(), L1Triplet::compareChain);
      tr.chi2Chain = (tr.GetChi2())+BestChi2;
  //  (vTriplets_part[vTriplets_part.size()-1])->SetLevel( level );
    TripletsLocalReal[istal][ip_cur][TripNPortion[istal][ip_cur]-1].SetLevel( level );
    
  }
    
    
    
       std::sort(TripletsLocalRealAll[istal][ip_cur].begin(), TripletsLocalRealAll[istal][ip_cur].begin() + TripNPortion[istal][ip_cur], L1Triplet::compareLevel);
} 







/// Add the middle hits to parameters estimation. Propagate to right station.
/// Find the triplets(right hit). Reformat data in the portion of triplets.
inline void L1Algo::f30SIMD(  // input
                        L1Station &stam, L1Station &star,
                        
                        int istam, int istar,
                       
                        L1TrackParV *T_1V, L1TrackParV *T_1VBuf2, L1TrackParV *T_1VBuf, L1FieldRegionV *fld_1V,
                        L1TrackParV &TripletsArray,
                        
                        
                        Tindex n2,
                        
                        
                        //                 const vector<bool> &mrDuplets,
                        
                        // output
                        Tindex &n3,
                      
                      
                        nsL1::vector <int>::TSimd &MHits,
                        nsL1::vector <int>::TSimd &LHits,
                        nsL1::vector <int>::TSimd &RHits,
                        nsL1::vector<fvec>::TSimd &u_front_3V, nsL1::vector<fvec>::TSimd &u_back_3V, nsL1::vector<fvec>::TSimd &z_Pos_3V,
                            nsL1::vector <L1TrackPar>::TSimd &TripletsArrayV,
                            nsL1::vector <int>::TSimd &LHitsV,
                            nsL1::vector <int>::TSimd &MHitsV,
                            nsL1::vector <int>::TSimd &LHits3
                        
                        )
{

    
   
    
//     nsL1::vector <int>::TSimd LastHitsIndexes;
//     nsL1::vector <int>::TSimd FirstHitsIndexes;
//     nsL1::vector <int>::TSimd FirstHitsIndexesLoc;
//     nsL1::vector <int>::TSimd MiddleHitsIndexes;
//     RHits.resize(0);
//     int_v RHitsL, MHitsL, LHitsL;
    
//     L1TrackParV TempTracks;
    L1TrackPar T;
    L1TrackPar T2;
    L1FieldRegion f;
    
    int_v hitR_V, hitM_V, hitL_V3, hitL_V;
    int_v hitR_Vl, hitM_Vl, hitL_V3l, hitL_Vl;
  //  L1TrackPar T2;
    L1TrackPar T3;
    fvec u, v, z;
    int_m IsNiceDoublet (true);
   
    int fill =0;
    int fill2 =0;
    n3 =0;
    
    if (istar < NStations){


      
        
        for (Tindex i2 = 0; i2 < n2; i2+=fvecLen) {
            
           
            IsNiceDoublet=int_m(true);
            
            

                if((i2 + fvecLen)>MHits.size())
                IsNiceDoublet &= int_m( int_v(Vc::IndexesFromZero) < int(MHits.size()%fvecLen) );
            
                int_v hitM_V = reinterpret_cast<int_v&>(MHits[i2]); // vector middle hits indexes in the initial L1HitPonts array
                
                
                
                hitM_V(!(IsNiceDoublet))=0;
                
                
                const fvec zm(vStsHitPointsUnusedV[ &stam - vStations ].Z(), hitM_V);
                const fvec um(vStsHitPointsUnusedV[ &stam - vStations ].U(), hitM_V);
                const fvec vm(vStsHitPointsUnusedV[ &stam - vStations ].V(), hitM_V);
                
                
                
                
                
                
                
                
                int_v hitL_V = reinterpret_cast<int_v&>(LHits[i2]);
                
                int_v hitL_V3 = reinterpret_cast<int_v&>(LHits3[i2]);
                
               
                  hitL_V(!(IsNiceDoublet))=0;
                 hitL_V3(!(IsNiceDoublet))=0;
                
                fld_1V->GetFieldPar(f, hitL_V3);
                T_1VBuf->GetTrackPar(T, hitL_V);

// for (int i=0; i<fvecLen;i++){
// cout<<f.cx1[i]<<" "<<hitL_V[i]<<" "<<hitM_V[i]<<" "<<n2<<" simd"<<endl;
// }
          
                
              //  T_1VBuf2->CopyTrackPar((*T_1VBuf), hitL_V, i2);
               
              //  L1TrackParR T (T_1VBuf, i2);
              
                
                L1ExtrapolateLine( T, zm );
                L1Filter( T, stam.frontInfo, um );
                L1Filter( T, stam.backInfo,  vm );
                
                L1AddMaterial( T, stam.materialInfo, T.qp );
                if ( (istar >= NMvdStations) && (istam <= NMvdStations - 1) ) L1AddPipeMaterial( T, T.qp );
                
                
                // extrapolate to the right hit station
                L1Extrapolate( T, star.z, T.qp, f);
                
                // ---- Find the triplets(right hit). Reformat data in the portion of triplets. ----
                
                
                
                
                  IsNiceDoublet &=(( T.C00 >= fvec(Vc::Zero) &&  T.C11 >= fvec(Vc::Zero) &&  T.C22 >= fvec(Vc::Zero) &&  T.C33 >= fvec(Vc::Zero) &&  T.C44 >= fvec(Vc::Zero) ));
                  if ((IsNiceDoublet.isEmpty())) continue;
                
                
                
                
                fvec Pick_r22 = (TRIPLET_CHI2_CUT - T.chi2);

                
                const fvec iz = 1/T.z;
                
                fvec x = T.x*iz;
                fvec y = T.y*iz;
                fvec dx = (sqrt(Pick_r22*(T.C00 + stam.XYInfo.C00))+MaxDZ*abs(T.tx))*iz;
                fvec dy = (sqrt(Pick_r22*(T.C11 + stam.XYInfo.C11))+MaxDZ*abs(T.ty))*iz;
                
                
                
                L1HitAreaV areaV( vGrid[ &star - vStations ], x, y, dx, dy);


            

                
                areaV.HitsIndexesMiddle(i2, hitL_V, hitM_V, hitL_V3, fill, fill2, T, vStsHitPointsUnusedV[ &star - vStations ], TRIPLET_CHI2_CUT, star, IsNiceDoublet, n3, TripletsArrayV, u_front_3V, u_back_3V, z_Pos_3V,
                                        hitR_V, hitM_V, hitL_V3, hitL_V, hitR_Vl, hitM_Vl, hitL_V3l, hitL_Vl, T2, T3, u, v, z);

           
              //  areaV.HitsIndexesMiddle(LastHitsIndexes, FirstHitsIndexes, MiddleHitsIndexes, FirstHitsIndexesLoc, i2, hitL_V, hitM_V, hitL_V3);
               
                // find first possible hit
     
        }
        

        
        
  
        
        if (fill!=0)  {
            
            int_m IsNiceDoublet (true);
            
            
            IsNiceDoublet&=(int_v(Vc::IndexesFromZero)<int_v(fill));
            
            fvec zr(vStsHitPointsUnusedV[ &star - vStations ].Z(), hitR_V);
            fvec yr(vStsHitPointsUnusedV[ &star - vStations ].Y(), hitR_V);
            fvec xr(vStsHitPointsUnusedV[ &star - vStations ].X(), hitR_V);
            fvec ur(vStsHitPointsUnusedV[ &star - vStations ].U(), hitR_V);
            fvec vr(vStsHitPointsUnusedV[ &star - vStations ].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())) {
            
        
            
            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())) {
            
            // 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())) {
            
            
            
            
            for (int iV=0; iV<fill; 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);
                }
                
                
            }
            
        }
}
}
}
 if (fill2!=0) {TripletsArrayV.push_back(T3);
                    u_front_3V.push_back(u);
                    u_back_3V.push_back(v);
                    z_Pos_3V.push_back(z);
        n3+=fill2;
        
        }


       // cout<<n3<<" n3 simd"<<endl;
 
/*        
        L1TrackPar T2;
        L1TrackPar T_buf;
        int TrackParFilled=0;
        int Triplets=0;
        
        fvec z, u, v = fvec (0.f);
        
        for(int iPair=0; iPair<LastHitsIndexes.size(); iPair+=fvecLen)
        {
            int_v hitR_V = reinterpret_cast<int_v&>(LastHitsIndexes[iPair]); // vector middle hits indexes in the initial L1HitPonts array
            int_m IsNiceDoublet (true);
            
          
            if((iPair + fvecLen)> LastHitsIndexes.size())
                IsNiceDoublet &= int_m( int_v(Vc::IndexesFromZero) < int(LastHitsIndexes.size()%fvecLen) );
            
            int_v hitL_V3 = reinterpret_cast<int_v&>(FirstHitsIndexesLoc[iPair]);
            
            
            hitR_V(!(IsNiceDoublet))=0;
            
            hitL_V3(!(IsNiceDoublet))=0;
            
            
            
            
            fvec zr(vStsHitPointsUnusedV[ &star - vStations ].Z(), hitR_V);
            fvec yr(vStsHitPointsUnusedV[ &star - vStations ].Y(), hitR_V);
            fvec xr(vStsHitPointsUnusedV[ &star - vStations ].X(), hitR_V);
            fvec ur(vStsHitPointsUnusedV[ &star - vStations ].U(), hitR_V);
            fvec vr(vStsHitPointsUnusedV[ &star - vStations ].V(), hitR_V);
            
            
            
            
            
            
           // int_v index = (int_v(Vc::IndexesFromZero)+iPair);
            
            //TempTracks.GetTrackPar(T2, index);
            
      //       T_1VBuf2->GetTrackPar(T2, hitL_V3);
            
            T_1VBuf->GetTrackPar(T2, hitL_V3);
            
  //          T_1V->CopyTrackPar((*T_1VBuf2), hitL_V3, 0);
//            
//            
    //        L1TrackParR T (T_1VBuf2, 0);
            
            
            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())) continue;
            
            
            
            
            
            
            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())) continue;
            
            // 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())) continue;
            
            int_v hitM_V = reinterpret_cast<int_v&>(MiddleHitsIndexes[iPair]);
            
            int_v hitL_V = reinterpret_cast<int_v&>(FirstHitsIndexes[iPair]);
           
            
            
            for (int iV=0; iV<fvecLen; iV++)
            {
                if (!(IsNiceDoublet[iV])) continue;
                
//               
//              //  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++;
//                    
//                }
                
                
            }



            
        }*/
        
//          n3=(TrackParFilled+RHits.size());
//        
//        if (TrackParFilled!=0){
//           // TripletsArray.Set(T_buf);
//        //   TripletsArrayV.push_back(T_buf);
//            RHits.resize(RHits.size()+fvecLen);
//            MHitsV.resize(MHitsV.size()+fvecLen);
//            LHitsV.resize(LHitsV.size()+fvecLen);
//            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&>(MHitsV[Triplets*fvecLen])=MHitsL;
//            reinterpret_cast<int_v&>(LHitsV[Triplets*fvecLen])=LHitsL;
//            Triplets+=TrackParFilled;
//            
//        }
      
        
    }  // if istar
    
   
    
}


/// Add the right hits to parameters estimation.
inline void L1Algo::f31SIMD(  // input
                        Tindex n3_V,
                        L1Station &star,
                            nsL1::vector<fvec>::TSimd &u_front_3, nsL1::vector<fvec>::TSimd &u_back_3, nsL1::vector<fvec>::TSimd &z_Pos_3,
                            // output

                        //                L1TrackPar *T_3
                        L1TrackParV &TripletsArray,
                            
                             nsL1::vector <L1TrackPar>::TSimd &TripletsArrayV
                        )
{
  
    for( Tindex i3_V=0; i3_V<n3_V; ++i3_V){
    L1ExtrapolateLine( TripletsArrayV[i3_V], z_Pos_3[i3_V] );
        L1Filter( TripletsArrayV[i3_V], star.frontInfo, u_front_3[i3_V] );    // 2.1/100 sec
        L1Filter( TripletsArrayV[i3_V], star.backInfo,  u_back_3 [i3_V] );   // 2.0/100 sec
        //     cout<<sqrt(T_3[i3_V].C00)<<endl;
        //     cout<<sqrt(T_3[i3_V].C10)<<endl;
        //     cout<<sqrt(T_3[i3_V].C11)<<endl;
        
       
    };
    
}


/// Refit Triplets.
inline void L1Algo::f32SIMD( // input // TODO not updated after gaps introduction
                        Tindex n3, int istal, L1TrackParV &TripletsArray,
                        
                      
                        
                        nsL1::vector <int>::TSimd &MHits,
                        nsL1::vector <int>::TSimd &LHits,
                        nsL1::vector <int>::TSimd &RHits,
                        int nIterations
                        )
{
//    const int NHits = 3; // triplets
//    
//    // prepare data
//    int ista[NHits] = {
//        istal,
//        istal + 1,
//        istal + 2
//    };
//    
//    L1Station sta[3];
//    for (int is = 0; is < NHits; ++is){
//        sta[is] = vStations[ista[is]];
//    };
//    
//    for( int i3=0; i3<n3; ++i3){
//        int i3_V = i3/fvecLen;
//        int i3_4 = i3%fvecLen;
//        
//        L1TrackPar &T3 = T_3[i3_V];
//        fscal qp0 = T3.qp[i3_4];
//        
//        // prepare data
//        THitI ihit[NHits] = {
//            (*RealIHitP)[hitsl_3[i3] + StsHitsUnusedStartIndex[ista[0]]],
//            (*RealIHitP)[hitsm_3[i3] + StsHitsUnusedStartIndex[ista[1]]],
//            (*RealIHitP)[hitsr_3[i3] + StsHitsUnusedStartIndex[ista[2]]]
//        };
//        
//        fscal u[NHits], v[NHits], x[NHits], y[NHits], z[NHits];
//        for (int ih = 0; ih < NHits; ++ih){
//            const L1StsHit &hit = (*vStsHits)[ihit[ih]];
//            u[ih] = (*vStsStrips)[hit.f];
//            v[ih] = (*vStsStripsB)[hit.b];
//            StripsToCoor(u[ih], v[ih], x[ih], y[ih], sta[ih]);
//            z[ih] = (*vStsZPos)[hit.iz];
//        };
//        
//        // initialize parameters
//        L1TrackPar T;
//        
//        const fvec vINF = .1f;
//        T.x  = x[0];
//        T.y  = y[0];
//        
//        fvec dz01 = 1.f/(z[1]-z[0]);
//        T.tx = (x[1]-x[0])*dz01;
//        T.ty = (y[1]-y[0])*dz01;
//        
//        T.qp = qp0;
//        T.z  = z[0];
//        T.chi2 = 0.f;
//        T.NDF = 2.f;
//        T.C00 = sta[0].XYInfo.C00;
//        T.C10 = sta[0].XYInfo.C10;
//        T.C11 = sta[0].XYInfo.C11;
//        
//        T.C20 = T.C21 = 0;
//        T.C30 = T.C31 = T.C32 = 0;
//        T.C40 = T.C41 = T.C42 = T.C43 = 0;
//        T.C22 = T.C33 = vINF;
//        T.C44 = 1.;
//        
//        // find field along the track
//        L1FieldValue B[3] _fvecalignment;
//        L1FieldRegion fld _fvecalignment;
//        
//        fvec tx[3] = {
//            (x[1]-x[0])/(z[1]-z[0]),
//            (x[2]-x[0])/(z[2]-z[0]),
//            (x[2]-x[1])/(z[2]-z[1])
//        };
//        fvec ty[3] = {
//            (y[1]-y[0])/(z[1]-z[0]),
//            (y[2]-y[0])/(z[2]-z[0]),
//            (y[2]-y[1])/(z[2]-z[1])
//        };
//        for (int ih = 0; ih < NHits; ++ih){
//            fvec dz = (sta[ih].z - z[ih]);
//            sta[ih].fieldSlice.GetFieldValue( x[ih] + tx[ih]*dz, y[ih] + ty[ih]*dz, B[ih] );
//        };
//        fld.Set( B[0], sta[0].z, B[1], sta[1].z, B[2], sta[2].z );
//        
//        // fit
//        for( int ih = 1; ih < NHits; ++ih){
//            L1Extrapolate( T, z[ih], T.qp, fld );
//            L1AddMaterial( T, sta[ih].materialInfo, T.qp );
//            if (ista[ih] == NMvdStations - 1) L1AddPipeMaterial( T, T.qp );
//            
//            L1Filter( T, sta[ih].frontInfo, u[ih] );
//            L1Filter( T, sta[ih].backInfo,  v[ih] );
//        }
//        
//        // repeat several times in order to improve precision
//        for (int iiter = 0; iiter < nIterations; ++iiter){
//            // fit backward
//            // keep tx,ty,q/p
//            int ih = NHits-1;
//            T.x  = x[ih];
//            T.y  = y[ih];
//            T.z  = z[ih];
//            T.chi2 = 0.;
//            T.NDF = 2.;
//            T.C00 = sta[ih].XYInfo.C00;
//            T.C10 = sta[ih].XYInfo.C10;
//            T.C11 = sta[ih].XYInfo.C11;
//            
//            T.C20 = T.C21 = 0;
//            T.C30 = T.C31 = T.C32 = 0;
//            T.C40 = T.C41 = T.C42 = T.C43 = 0;
//            T.C22 = T.C33 = vINF;
//            T.C44 = 1.;
//            
//            //       L1Filter( T, sta[ih].frontInfo, u[ih] );
//            //       L1Filter( T, sta[ih].backInfo,  v[ih] );
//            for( ih = NHits-2; ih >= 0; ih--){
//                L1Extrapolate( T, z[ih], T.qp, fld );
//                L1AddMaterial( T, sta[ih].materialInfo, T.qp );
//                if (ista[ih] == NMvdStations) L1AddPipeMaterial( T, T.qp );
//                
//                L1Filter( T, sta[ih].frontInfo, u[ih] );
//                L1Filter( T, sta[ih].backInfo,  v[ih] );
//            }
//            // fit forward
//            ih = 0;
//            T.x  = x[ih];
//            T.y  = y[ih];
//            T.z  = z[ih];
//            T.chi2 = 0.;
//            T.NDF = 2.;
//            T.C00 = sta[ih].XYInfo.C00;
//            T.C10 = sta[ih].XYInfo.C10;
//            T.C11 = sta[ih].XYInfo.C11;
//            
//            T.C20 = T.C21 = 0;
//            T.C30 = T.C31 = T.C32 = 0;
//            T.C40 = T.C41 = T.C42 = T.C43 = 0;
//            T.C22 = T.C33 = vINF;
//            T.C44 = 1.;
//            
//            //       L1Filter( T, sta[ih].frontInfo, u[ih] );
//            //       L1Filter( T, sta[ih].backInfo,  v[ih] );
//            for( ih = 1; ih < NHits; ++ih){
//                L1Extrapolate( T, z[ih], T.qp, fld );
//                L1AddMaterial( T, sta[ih].materialInfo, T.qp );
//                if (ista[ih] == NMvdStations + 1) L1AddPipeMaterial( T, T.qp );
//                
//                L1Filter( T, sta[ih].frontInfo, u[ih] );
//                L1Filter( T, sta[ih].backInfo,  v[ih] );
//            }
//        } // for iiter
//        
//        T3.SetOneEntry(i3_4, T,i3_4);
//    }//i3
//    
//    
//    //simd here
//    
//    L1TrackPar T3;
//    
//    for( int i3=0; i3<n3; i3+=fvecLen){
//        
//        
//        
//        int_v index = (int_v(Vc::IndexesFromZero)+i3);
//        
//        TripletsArray.GetTrackPar(T3, index);
//        
//        fvec qp0 = T3.qp;
//        
//        
//        int_v hitL_V = reinterpret_cast<int_v&>(LHits[i3]); // vector middle hits indexes in the initial L1HitPonts array
//        int_v hitM_V = reinterpret_cast<int_v&>(MHits[i3]);
//        int_v hitR_V = reinterpret_cast<int_v&>(RHits[i3]);
//        int_m IsNiceDoublet (true);
//        
//        
//        if((i3 + fvecLen)> n3)
//            IsNiceDoublet &= int_m( int_v(Vc::IndexesFromZero) < int(n3%fvecLen) );
//        
//        
//        
//        hitR_V(!(IsNiceDoublet))=0;
//        hitM_V(!(IsNiceDoublet))=0;
//        hitL_V(!(IsNiceDoublet))=0;
//        
//        fvec z[NHits];
//        fvec y[NHits];
//        fvec x[NHits];
//        fvec u[NHits];
//        fvec v[NHits];
//        
//        
//        
//        z[0]=(vStsHitPointsUnusedV[ ista[0] ].Z(), hitR_V);
//        y[0]=(vStsHitPointsUnusedV[ ista[0] ].Y(), hitR_V);
//        x[0]=(vStsHitPointsUnusedV[ ista[0] ].X(), hitR_V);
//        u[0]=(vStsHitPointsUnusedV[ ista[0] ].U(), hitR_V);
//        v[0]=(vStsHitPointsUnusedV[ ista[0] ].V(), hitR_V);
//        
//        z[1]=(vStsHitPointsUnusedV[ ista[1] ].Z(), hitM_V);
//        y[1]=(vStsHitPointsUnusedV[ ista[1] ].Y(), hitM_V);
//        x[1]=(vStsHitPointsUnusedV[ ista[1] ].X(), hitM_V);
//        u[1]=(vStsHitPointsUnusedV[ ista[1] ].U(), hitM_V);
//        v[1]=(vStsHitPointsUnusedV[ ista[1] ].V(), hitM_V);
//        
//        z[2]=(vStsHitPointsUnusedV[ ista[2] ].Z(), hitL_V);
//        y[2]=(vStsHitPointsUnusedV[ ista[2] ].Y(), hitL_V);
//        x[2]=(vStsHitPointsUnusedV[ ista[2] ].X(), hitL_V);
//        u[2]=(vStsHitPointsUnusedV[ ista[2] ].U(), hitL_V);
//        v[2]=(vStsHitPointsUnusedV[ ista[2] ].V(), hitL_V);
//        
//        
//        
//        //        // prepare data
//        //        THitI ihit[NHits] = {
//        //            (*RealIHitP)[hitsl_3[i3] + StsHitsUnusedStartIndex[ista[0]]],
//        //            (*RealIHitP)[hitsm_3[i3] + StsHitsUnusedStartIndex[ista[1]]],
//        //            (*RealIHitP)[hitsr_3[i3] + StsHitsUnusedStartIndex[ista[2]]]
//        //        };
//        //
//        //        fscal u[NHits], v[NHits], x[NHits], y[NHits], z[NHits];
//        //        for (int ih = 0; ih < NHits; ++ih){
//        //            const L1StsHit &hit = (*vStsHits)[ihit[ih]];
//        //            u[ih] = (*vStsStrips)[hit.f];
//        //            v[ih] = (*vStsStripsB)[hit.b];
//        //            StripsToCoor(u[ih], v[ih], x[ih], y[ih], sta[ih]);
//        //            z[ih] = (*vStsZPos)[hit.iz];
//        //        };
//        
//        // initialize parameters
//        L1TrackPar T;
//        
//        const fvec vINF = .1f;
//        T.x  = x[0];
//        T.y  = y[0];
//        
//        fvec dz01 = 1.f/(z[1]-z[0]);
//        T.tx = (x[1]-x[0])*dz01;
//        T.ty = (y[1]-y[0])*dz01;
//        
//        T.qp = qp0;
//        T.z  = z[0];
//        T.chi2 = 0.f;
//        T.NDF = 2.f;
//        T.C00 = sta[0].XYInfo.C00;
//        T.C10 = sta[0].XYInfo.C10;
//        T.C11 = sta[0].XYInfo.C11;
//        
//        T.C20 = T.C21 = 0.f;
//        T.C30 = T.C31 = T.C32 = 0.f;
//        T.C40 = T.C41 = T.C42 = T.C43 = 0.f;
//        T.C22 = T.C33 = vINF;
//        T.C44 = 1.f;
//        
//        // find field along the track
//        L1FieldValue B[3] _fvecalignment;
//        L1FieldRegion fld _fvecalignment;
//        
//        fvec tx[3] = {
//            (x[1]-x[0])/(z[1]-z[0]),
//            (x[2]-x[0])/(z[2]-z[0]),
//            (x[2]-x[1])/(z[2]-z[1])
//        };
//        fvec ty[3] = {
//            (y[1]-y[0])/(z[1]-z[0]),
//            (y[2]-y[0])/(z[2]-z[0]),
//            (y[2]-y[1])/(z[2]-z[1])
//        };
//        for (int ih = 0; ih < NHits; ++ih){
//            fvec dz = (sta[ih].z - z[ih]);
//            sta[ih].fieldSlice.GetFieldValue( x[ih] + tx[ih]*dz, y[ih] + ty[ih]*dz, B[ih] );
//        };
//        fld.Set( B[0], sta[0].z, B[1], sta[1].z, B[2], sta[2].z );
//        
//        // fit
//        for( int ih = 1; ih < NHits; ++ih){
//            L1Extrapolate( T, z[ih], T.qp, fld );
//            L1AddMaterial( T, sta[ih].materialInfo, T.qp );
//            if (ista[ih] == NMvdStations - 1) L1AddPipeMaterial( T, T.qp );
//            
//            L1Filter( T, sta[ih].frontInfo, u[ih] );
//            L1Filter( T, sta[ih].backInfo,  v[ih] );
//        }
//        
//        // repeat several times in order to improve precision
//        for (int iiter = 0; iiter < nIterations; ++iiter){
//            // fit backward
//            // keep tx,ty,q/p
//            int ih = NHits-1;
//            T.x  = x[ih];
//            T.y  = y[ih];
//            T.z  = z[ih];
//            T.chi2 = 0.f;
//            T.NDF = 2.f;
//            T.C00 = sta[ih].XYInfo.C00;
//            T.C10 = sta[ih].XYInfo.C10;
//            T.C11 = sta[ih].XYInfo.C11;
//            
//            T.C20 = T.C21 = 0.f;
//            T.C30 = T.C31 = T.C32 = 0.f;
//            T.C40 = T.C41 = T.C42 = T.C43 = 0.f;
//            T.C22 = T.C33 = vINF;
//            T.C44 = 1.f;
//            
//            //       L1Filter( T, sta[ih].frontInfo, u[ih] );
//            //       L1Filter( T, sta[ih].backInfo,  v[ih] );
//            for( ih = NHits-2; ih >= 0; ih--){
//                L1Extrapolate( T, z[ih], T.qp, fld );
//                L1AddMaterial( T, sta[ih].materialInfo, T.qp );
//                if (ista[ih] == NMvdStations) L1AddPipeMaterial( T, T.qp );
//                
//                L1Filter( T, sta[ih].frontInfo, u[ih] );
//                L1Filter( T, sta[ih].backInfo,  v[ih] );
//            }
//            // fit forward
//            ih = 0;
//            T.x  = x[ih];
//            T.y  = y[ih];
//            T.z  = z[ih];
//            T.chi2 = 0.f;
//            T.NDF = 2.f;
//            T.C00 = sta[ih].XYInfo.C00;
//            T.C10 = sta[ih].XYInfo.C10;
//            T.C11 = sta[ih].XYInfo.C11;
//            
//            T.C20 = T.C21 = 0.f;
//            T.C30 = T.C31 = T.C32 = 0.f;
//            T.C40 = T.C41 = T.C42 = T.C43 = 0.f;
//            T.C22 = T.C33 = vINF;
//            T.C44 = 1.f;
//            
//            //       L1Filter( T, sta[ih].frontInfo, u[ih] );
//            //       L1Filter( T, sta[ih].backInfo,  v[ih] );
//            for( ih = 1; ih < NHits; ++ih){
//                L1Extrapolate( T, z[ih], T.qp, fld );
//                L1AddMaterial( T, sta[ih].materialInfo, T.qp );
//                if (ista[ih] == NMvdStations + 1) L1AddPipeMaterial( T, T.qp );
//                
//                L1Filter( T, sta[ih].frontInfo, u[ih] );
//                L1Filter( T, sta[ih].backInfo,  v[ih] );
//            }
//        } // for iiter
//        
//        TripletsArray.Set(i3, T);
//        
//        /// T3.SetOneEntry(i3_4, T,i3_4); check if it works
//    }//i3
//    
//    
//    
//    
    
    
} // f32


/// Select triplets. Save them into vTriplets.
inline void L1Algo::f4SIMD(  // input
                       Tindex n3, int istal, int istam, int istar, L1TrackParV &TripletsArray,
                       
                       
                                            // output
                       Tindex &nstaltriplets,
                       
                       Tindex ip_cur,
                           nsL1::vector <L1TrackPar>::TSimd &TripletsArrayV,
                           nsL1::vector <int>::TSimd &MHits,
                           nsL1::vector <int>::TSimd &LHits,
                           nsL1::vector <int>::TSimd &RHits,
                           nsL1::vector <int>::TSimd &LHitsV,
                           nsL1::vector <int>::TSimd &MHitsV,
                           nsL1::vector <int>::TSimd &LHits3,
                           Tindex n3S

                       
                       
                       )
{
    
   
   
   
    
    // simd here
     for( Tindex i3_V=0; i3_V<n3; i3_V++)
        {
           
            
           L1TrackPar &T3 = TripletsArrayV[i3_V];

    
            int_m IsNiceTriplet (true);

    
            if(((i3_V+1)*fvecLen)> n3S) IsNiceTriplet &= (int_v(Vc::IndexesFromZero) < int((n3S+fvecLen)%fvecLen));
            
          
    

            fvec chi2 = T3.chi2;
    
    
            IsNiceTriplet &=(!( !isfinite(chi2) || chi2 < 0.f || chi2 > fvec(TRIPLET_CHI2_CUT) ) );
            
            
    
    
            if ((IsNiceTriplet.isEmpty())) continue;

            fvec qp = T3.qp;

            fvec Cqp = 5.f*sqrt(abs(T3.C44));
    
            
            
            
            int_v ihitl = reinterpret_cast<int_v&> (LHits3[i3_V*fvecLen]) + int_v( StsHitsUnusedStartIndex[istal]);
            int_v ihitm = reinterpret_cast<int_v&> (MHitsV[i3_V*fvecLen]) + int_v (StsHitsUnusedStartIndex[istam]);
            int_v ihitr = reinterpret_cast<int_v&> (RHits[i3_V*fvecLen]) + int_v (StsHitsUnusedStartIndex[istar]);
    

//            ihitl(!(IsNiceTriplet))=0;
//            ihitm(!(IsNiceTriplet))=0;
//            ihitr(!(IsNiceTriplet))=0;
//            
//            qp(!(IsNiceTriplet))=0;
//            Cqp(!(IsNiceTriplet))=0;

            
           
    
             L1TripletV  tr;
    
            if (istal!=FIRSTCASTATION) {

                tr.Set( ihitl, ihitm, ihitr,
                       int_v(istal), int_v(istam), int_v(istar),
                       int_v(0.f), qp, chi2);

                tr.Cqp   = ( Cqp );
                tr.chi2Chain=tr.GetChi2();
    


                
    
                
                for(int iV=0; iV<fvecLen; iV++){
                    if (!(IsNiceTriplet[iV])) continue;
                    
                    
                    L1Triplet & trip=TripletsLocalRealV[istal][ip_cur][TripNPortionV[istal][ip_cur]];
                    
                    trip.SetOneEntry(tr, iV);
                    
                    TripletsLocalRealAllV[istal][ip_cur][TripNPortionV[istal][ip_cur]]=&trip;
                    
                     TripForHitV[tr.GetLHit()[iV]].push_back(&trip);
                    
                   
                    TripNPortionV[istal][ip_cur]++;
                   // ++nstaltriplets;
                   

                }
            }
            

        
            if (istal > (NStations - 4)) continue;
            
            int_v level (Vc::Zero);
            
            int_v size (Vc::One);
            
            vector<L1Triplet *> neighCands[fvecLen]; // save neighbour candidates
           for(int iV=0; iV<fvecLen; iV++) neighCands[iV].reserve(8); // ~average is 1-2 for central, up to 5
            
            for(int iV=0; iV<fvecLen; iV++){
                if (!(IsNiceTriplet[iV])) continue;
                
                if ( TripForHitV[ihitm[iV]].size() == 0) {size[iV]=0.f; continue;}
               
                
                for (unsigned int iN=0; iN < TripForHitV[ihitm[iV]].size(); ++iN){

                    L1Triplet* &curNeighbour = TripForHitV[ihitm[iV]][iN];
                    
                    if (curNeighbour->GetMHit() != ihitr[iV]) continue; // neighbours should have 2 common hits
                    
                    L1Triplet* &tripn = curNeighbour;
                    
                    fscal qp2 = tripn->GetQp();
                    fscal Cqp2 = tripn->Cqp;
                    if ( abs(qp - qp2) > PickNeighbour * (Cqp + Cqp2) )  continue; // neighbours should have same qp
                    
                    // calculate level
                    unsigned char jlevel = tripn->GetLevel();
                    if ( level <= jlevel )
                        level[iV] = jlevel + 1;
                    
                  
                    
                    neighCands[iV].push_back(curNeighbour);
                    
                    
                }

                
                

            }
            
            
            
           // IsNiceTriplet &=(size!=0);
            
             if ((IsNiceTriplet.isEmpty())) continue;
            
            
            if (istal==FIRSTCASTATION) {
                
              //  IsNiceTriplet &=(level!=0);
                
                
                if ((IsNiceTriplet.isEmpty())) continue;
               // cout<<ihitl<<" ihitl"<<endl;
                
                tr.Set( ihitl, ihitm, ihitr,
                       int_v(istal), int_v(istam), int_v(istar),
                       int_v(0.f), qp, chi2);
                
                tr.Cqp   = ( Cqp );
                tr.chi2Chain=tr.GetChi2();
                
          
                
                for(int iV=0; iV<fvecLen; iV++){
                    
                    if (!(IsNiceTriplet[iV])) continue;
                    
                    L1Triplet & trip=TripletsLocalRealV[istal][ip_cur][TripNPortionV[istal][ip_cur]];
                    
                    trip.SetOneEntry(tr, iV);
                    
                    TripletsLocalRealAllV[istal][ip_cur][TripNPortionV[istal][ip_cur]]=&trip;
                    
                    TripForHitV[tr.GetLHit()[iV]].push_back(&trip);
                    
                    TripNPortionV[istal][ip_cur]++;
                   

                }
                
                
            }

           
            
            for(int iV=0; iV<fvecLen; iV++){
            fscal BestChi2 = 1000000000;
            for (unsigned int in = 0; in < (neighCands[iV].size()); ++in) {
                const Tindex nLevel = neighCands[iV][in]->GetLevel();
                
                if (level[iV] == nLevel + 1) {
                   
                    TripletsLocalRealV[istal][ip_cur][TripNPortionV[istal][ip_cur]-1+iV].neighbours.push_back(neighCands[iV][in]);

                    
                   if ((neighCands[iV][in]->chi2Chain)<BestChi2)
                    {
                        BestChi2=(neighCands[iV][in]->chi2Chain);
                        
                       
                        
                    }
                    
                    
                }
                
            }
            
              
                
                std::sort(TripletsLocalRealV[istal][ip_cur][TripNPortionV[istal][ip_cur]-1+iV].neighbours.begin(), TripletsLocalRealV[istal][ip_cur][TripNPortionV[istal][ip_cur]-1+iV].neighbours.end(), L1Triplet::compareChain);
                
                tr.chi2Chain = (tr.GetChi2()[iV])+BestChi2;
               
                TripletsLocalRealV[istal][ip_cur][TripNPortionV[istal][ip_cur]-1+iV].SetLevel( level[iV] );
                
            }
            
            

    
            
            
        }
    
    
    
    std::sort(TripletsLocalRealAllV[istal][ip_cur].begin(), TripletsLocalRealAllV[istal][ip_cur].begin() + TripNPortionV[istal][ip_cur], L1Triplet::compareLevel);
}

/// ------------------- doublets on station ----------------------

inline void L1Algo::DupletsStaPortSIMD(  /// creates duplets: input: @istal - start station number, @istam - last station number, @ip - index of portion, @&n_g1 - number of elements in portion, @*portionStopIndex
                                   int istal, int istam,
                                   Tindex ip,
                                   vector<Tindex> &n_g1, Tindex *portionStopIndex,
                                   
                                   /// output:
                                    /// @*T_1 - singlets perameters, @*fld_1 - field aproximation, @*hitsl_1- left hits of triplets, @&lmDuplets - existance of a doublet starting from the left hit,
                                   L1TrackParV *T_1V,
                                   L1TrackParV *T_1VBuf,
                                   L1TrackParV *T_1VBuf2,
                                   ///  @&n_2 - number of douplets,@&i1_2 - index of 1st hit in portion indexed by doublet index, @&hitsm_2 - index of middle hit in hits array indexed by doublet index
                                   L1FieldRegionV *fld_1V,
                            
                                   
                                   //   vector<bool> &lmDuplets,
                                   
                                   Tindex &n_2,
                            
                    
                                   nsL1::vector <int>::TSimd &MHits,
                                   nsL1::vector <int>::TSimd &LHits,
                                       nsL1::vector <int>::TSimd &LHits3
                                   )
{
    if ( istam < NStations ) {
        L1Station &stam = vStations[istam];
      
        
        Tindex n1 = n_g1[ip];

        
        Tindex n1_V = (n1+fvecLen-1)/fvecLen;
        
      
        
        int FirstHitInPortion = (ip - portionStopIndex[istal+1]) * Portion;
        int FirstHitInPortionV = (ip - portionStopIndex[istal+1]) * Portion/fvecLen;
        
        
       
       fvec* u_frontV = &(((vStsHitPointsUnusedV[istal]).u[0]))+FirstHitInPortionV;
       fvec* u_backV =&(((vStsHitPointsUnusedV[istal]).v[0]))+FirstHitInPortionV;
       fvec* zPosV = &(((vStsHitPointsUnusedV[istal]).z[0]))+FirstHitInPortionV;

       
       
        
        f11SIMD(istal, istam,
            n1_V,
            
            u_frontV, u_backV, zPosV,
            // output
            T_1V,fld_1V
            );
        
        
        
        f20SIMD(  // input
            n1, stam,
            T_1V,
             T_1VBuf,
                T_1VBuf2,
            MHits,
            LHits,
            //             hitsl_1,
            // output
            n_2,
            FirstHitInPortion,
            LHits3
            );
        


    }
}


       
                                /// ------------------- doublets on station ----------------------

inline void L1Algo::DupletsStaPort(  /// creates duplets: input: @istal - start station number, @istam - last station number, @ip - index of portion, @&n_g1 - number of elements in portion, @*portionStopIndex
  int istal, int istam,
  Tindex ip,
  vector<Tindex> &n_g1, Tindex *portionStopIndex,
  
    /// output: 
  L1TrackPar *T_1, /// @*T_1 - singlets perameters, @*fld_1 - field aproximation, @*hitsl_1- left hits of triplets, @&lmDuplets - existance of a doublet starting from the left hit, 
 
  L1FieldRegion *fld_1, ///  @&n_2 - number of douplets,@&i1_2 - index of 1st hit in portion indexed by doublet index, @&hitsm_2 - index of middle hit in hits array indexed by doublet index
 
  THitI *hitsl_1,
                       
//   vector<bool> &lmDuplets,
            
  Tindex &n_2,
  vector<THitI> &i1_2,
  vector<THitI> &hitsm_2
  )
{
  if ( istam < NStations ) {
    L1Station &stam = vStations[istam];
  
      // prepare data
    L1HitPoint *vStsHits_l = &((*vStsHitPointsUnused)[0]) + StsHitsUnusedStartIndex[istal];
    L1HitPoint *vStsHits_m = &((*vStsHitPointsUnused)[0]) + StsHitsUnusedStartIndex[istam];



        fvec u_front[Portion/fvecLen], u_back[Portion/fvecLen];
        fvec zPos[Portion/fvecLen];

          /// prepare the portion of left hits data

        Tindex n1 = n_g1[ip];
        
        f10(  // input
            (ip - portionStopIndex[istal+1]) * Portion, n1, vStsHits_l,
              // output
            u_front, u_back, zPos,
            hitsl_1
          );

        for (Tindex i = 0; i < n1; ++i)
          L1_ASSERT(hitsl_1[i] < StsHitsUnusedStopIndex[istal] - StsHitsUnusedStartIndex[istal],
            hitsl_1[i] << " < " << StsHitsUnusedStopIndex[istal] - StsHitsUnusedStartIndex[istal]);
        
        Tindex n1_V = (n1+fvecLen-1)/fvecLen;

          /// Get the field approximation. Add the target to parameters estimation. Propagaete to middle station.

        f11(istal, istam,
            n1_V,

            u_front, u_back, zPos,
              // output
            T_1, fld_1
          );

          /// Find the doublets. Reformat data in the portion of doublets.

#ifdef DOUB_PERFORMANCE
        vector<THitI> hitsl_2;
#endif // DOUB_PERFORMANCE
          
      
      
        f20(  // input
            n1, stam,
             vStsHits_m, 
            T_1,
          
//             hitsl_1,
              // output
            n_2,
            i1_2,
#ifdef DOUB_PERFORMANCE
            hitsl_2,
#endif // DOUB_PERFORMANCE
            hitsm_2
//             ,            lmDuplets
           );
     

        for (Tindex i = 0; i < static_cast<Tindex>(hitsm_2.size()); ++i)
         L1_ASSERT(hitsm_2[i] < StsHitsUnusedStopIndex[istam] - StsHitsUnusedStartIndex[istam], hitsm_2[i] << " " << StsHitsUnusedStopIndex[istam] - StsHitsUnusedStartIndex[istam]); 

#ifdef DOUB_PERFORMANCE
        THitI* RealIHitL = &((*RealIHitP)[StsHitsUnusedStartIndex[istal]]);
        THitI* RealIHitM = &((*RealIHitP)[StsHitsUnusedStartIndex[istam]]);
        for (Tindex i = 0; i < n2; ++i){
          // int i_4 = i%4;
          // int i_V = i/4;
          THitI iHits[2] = {
            RealIHitL[hitsl_2[i]],
            RealIHitM[hitsm_2[i]]
          };
          fL1Eff_doublets->AddOne(iHits);
        }
#endif // DOUB_PERFORMANCE
  }
}



/// ------------------- Triplets on station ----------------------

inline void L1Algo::TripletsStaPortSIMD(  /// creates triplets: input: @istal - start station number, @istam - middle station number, @istar - last station number, @ip - index of portion, @&n_g1 - numer of elements in portion, @*portionStopIndex
                                    int istal, int istam, int istar,
                                    
                                    ///@nstaltriplets - , @*portionStopIndex, @*T_1 - track parameters for singlets, @*fld_1 - field approximation for singlets, @&n_2 - number of doublets in portion
                                    ///  @&n_2 - number of douplets,@&i1_2 - index of 1st hit in portion indexed by doublet index, @&hitsm_2 - index of middle hit in hits array indexed by doublet index
                                    
                                    
                                    Tindex& nstaltriplets,
                                 
                                    L1TrackParV *T_1V,
                                        L1TrackParV *T_1VBuf2,
                                        L1TrackParV *T_1VBuf,
                                   
                                    L1FieldRegionV *fld_1V,
                                   
                                    L1TrackParV &TripletsArray,
                                    
                                    Tindex &n_2,
                                  
                                    
                                    //                             const vector<bool> &mrDuplets,
                                    
                                    /// output: @*vTriplets_part - array of triplets, @*TripStartIndexH, @*TripStopIndexH - start/stop index of a triplet in the array
                                    
                                    Tindex ip_cur,
                                    nsL1::vector <int>::TSimd &MHits,
                                    nsL1::vector <int>::TSimd &LHits,
                                    nsL1::vector <int>::TSimd &RHits,
                                        nsL1::vector <L1TrackPar>::TSimd &TripletsArrayV,
                                        nsL1::vector <int>::TSimd &LHits3
                                    
                                    )
{
    
    if (istar < NStations )
    {
        // prepare data
        L1Station &stam = vStations[istam];
        L1Station &star = vStations[istar];
        
        L1HitPoint *vStsHits_m = &((*vStsHitPointsUnused)[0]) + StsHitsUnusedStartIndex[istam];
        
        L1HitPoint *vStsHits_r = 0;
        vStsHits_r = &((*vStsHitPointsUnused)[0]) + StsHitsUnusedStartIndex[istar];
        
        Tindex n3=0, n3_V;
        
     
        
      
        
     
        
        /// Add the middle hits to parameters estimation. Propagate to right station.
        
        
     
        nsL1::vector<fvec>::TSimd u_front3V, u_back3V, z_pos3V;
        
        nsL1::vector <int>::TSimd MHitsV;
        nsL1::vector <int>::TSimd LHitsV;
        
//        T_3.reserve(MaxPortionTriplets/fvecLen);
//        hitsl_3.reserve(MaxPortionTriplets);
//        hitsm_3.reserve(MaxPortionTriplets);
//        hitsr_3.reserve(MaxPortionTriplets);
//        u_front3.reserve(MaxPortionTriplets/fvecLen);
//        u_back3.reserve(MaxPortionTriplets/fvecLen);
//        z_pos3.reserve(MaxPortionTriplets/fvecLen);
        
        u_front3V.reserve(MaxPortionTriplets/fvecLen);
        u_back3V.reserve(MaxPortionTriplets/fvecLen);
        z_pos3V.reserve(MaxPortionTriplets/fvecLen);

    nsL1::vector<L1TrackPar>::TSimd T_3;
  
    vector<THitI> hitsl_3, hitsm_3, hitsr_3;

      /// Add the middle hits to parameters estimation. Propagate to right station.


        nsL1::vector<fvec>::TSimd u_front3, u_back3, z_pos3;

        T_3.reserve(MaxPortionTriplets/fvecLen);

        
        /// Find the triplets(right hit). Reformat data in the portion of triplets.
        
       f30SIMD(  // input
           stam, star,
           
           istam, istar,
           
           T_1V,  T_1VBuf2, T_1VBuf, fld_1V,
          
           TripletsArray,
           
           n_2,
          
           
           //             mrDuplets,
           // output
           n3,
               
          
           MHits,
           LHits,
           RHits,
           u_front3V, u_back3V, z_pos3V,
           T_3,
           LHitsV,
           MHitsV,
           LHits3
           );
        
   
// 
// for (int i=0; i<T_3.size(); i++) cout<<T_3[i].x<<" simd x"<<endl;


       //  cout<<u_front3V.size()<<" simd"<<endl;
        
        n3_V = (n3+fvecLen-1)/fvecLen;
       
        
        //        if (n3 >= MaxPortionTriplets) cout << "isec: " << isec << " stantion: " << istal << " portion number: " << ip << " CATrackFinder: Warning: Too many Triplets created in portion" << endl;
        
        
        
        
        /// Add the right hits to parameters estimation.
//        f31SIMD(  // input
//            n3_V,
//            star,
//            u_front3V, u_back3V, z_pos3V,
//            
//            // output
//           TripletsArray, TripletsArrayV
//                
//            );
        
        /// refit
        //         f32(  n3, istal, _RealIHit,  TripletsArray,         T_3,         hitsl_3, hitsm_3, hitsr_3, MHits, LHits, RHits );
        
        
#ifdef TRIP_PERFORMANCE
        THitI* RealIHitL = &((*RealIHitP)[StsHitsUnusedStartIndex[istal]]);
        THitI* RealIHitM = &((*RealIHitP)[StsHitsUnusedStartIndex[istam]]);
        THitI* RealIHitR = &((*RealIHitP)[StsHitsUnusedStartIndex[istar]]);
        for (Tindex i = 0; i < n3; ++i){
            Tindex i_4 = i%4;
            Tindex i_V = i/4;
            THitI iHits[3] = {
                RealIHitL[hitsl_3[i]],
                RealIHitM[hitsm_3[i]],
                RealIHitR[hitsr_3[i]]
            };
#ifdef PULLS
            if ( fL1Eff_triplets->AddOne(iHits) )
                fL1Pulls->AddOne(T_3[i_V], i_4, iHits[2]);
#else
            fL1Eff_triplets->AddOne(iHits);
#endif
            if (T_3[i_V].chi2[i_4] < TRIPLET_CHI2_CUT) 
                fL1Eff_triplets2->AddOne(iHits);
        }
#endif // TRIP_PERFORMANCE
        
        /// Fill Triplets.
        
//        f4SIMD( // input
//           n3_V, istal, istam, istar,  TripletsArray,
//         
//          
//           // output
//           nstaltriplets,
//        
//           ip_cur,
//        TripletsArrayV,
//        MHits,
//        LHits,
//               RHits,
//               LHitsV,
//               MHitsV,
//               LHits3,
//               n3
//
//        
//           
//           );
    }
}



                            /// ------------------- Triplets on station ----------------------

inline void L1Algo::TripletsStaPort(  /// creates triplets: input: @istal - start station number, @istam - middle station number, @istar - last station number, @ip - index of portion, @&n_g1 - numer of elements in portion, @*portionStopIndex
                            int istal, int istam, int istar,
              
              ///@nstaltriplets - , @*portionStopIndex, @*T_1 - track parameters for singlets, @*fld_1 - field approximation for singlets, @&n_2 - number of doublets in portion
              ///  @&n_2 - number of douplets,@&i1_2 - index of 1st hit in portion indexed by doublet index, @&hitsm_2 - index of middle hit in hits array indexed by doublet index


                            Tindex& nstaltriplets,
                            L1TrackPar *T_1,
                     
                            L1FieldRegion *fld_1,
                   
                            THitI *hitsl_1,
                  
                            
                            Tindex &n_2,
                            vector<THitI> &i1_2,
                            vector<THitI> &hitsm_2,

//                             const vector<bool> &mrDuplets,
                            
                              /// output: @*vTriplets_part - array of triplets, @*TripStartIndexH, @*TripStopIndexH - start/stop index of a triplet in the array
                           
                                       Tindex ip_cur
                    
                         
  )
{
  
  if (istar < NStations )
  {
      // prepare data
    L1Station &stam = vStations[istam];
    L1Station &star = vStations[istar];
  
    L1HitPoint *vStsHits_m = &((*vStsHitPointsUnused)[0]) + StsHitsUnusedStartIndex[istam];
  
    L1HitPoint *vStsHits_r = 0;
    vStsHits_r = &((*vStsHitPointsUnused)[0]) + StsHitsUnusedStartIndex[istar];
  
    Tindex n3=0, n3_V;
  
    nsL1::vector<L1TrackPar>::TSimd T_3;
  
    vector<THitI> hitsl_3, hitsm_3, hitsr_3;

      /// Add the middle hits to parameters estimation. Propagate to right station.


        nsL1::vector<fvec>::TSimd u_front3, u_back3, z_pos3;

        T_3.reserve(MaxPortionTriplets/fvecLen);
        hitsl_3.reserve(MaxPortionTriplets);
        hitsm_3.reserve(MaxPortionTriplets);
        hitsr_3.reserve(MaxPortionTriplets);
        u_front3.reserve(MaxPortionTriplets/fvecLen);
        u_back3.reserve(MaxPortionTriplets/fvecLen);
        z_pos3.reserve(MaxPortionTriplets/fvecLen);

          /// Find the triplets(right hit). Reformat data in the portion of triplets.

        f30(  // input
            vStsHits_r, stam, star,

            istam, istar,
            vStsHits_m,
            T_1, fld_1,
            hitsl_1,
           

            n_2,
            hitsm_2,
            i1_2,

//             mrDuplets,
              // output
            n3,
            T_3,
            hitsl_3, hitsm_3, hitsr_3,
            u_front3, u_back3, z_pos3
           
           );


// 
// for (int i=0; i<T_3.size(); i++)
// cout<<T_3[i].x<<" scal x"<<endl;
        n3_V = (n3+fvecLen-1)/fvecLen;

        for (Tindex i = 0; i < static_cast<Tindex>(hitsl_3.size()); ++i)
          L1_assert(hitsl_3[i] < StsHitsUnusedStopIndex[istal] - StsHitsUnusedStartIndex[istal]);
        for (Tindex i = 0; i < static_cast<Tindex>(hitsm_3.size()); ++i)
          L1_assert(hitsm_3[i] < StsHitsUnusedStopIndex[istam] - StsHitsUnusedStartIndex[istam]);
        for (Tindex i = 0; i < static_cast<Tindex>(hitsr_3.size()); ++i)
          L1_assert(hitsr_3[i] < StsHitsUnusedStopIndex[istar] - StsHitsUnusedStartIndex[istar]);

//        if (n3 >= MaxPortionTriplets) cout << "isec: " << isec << " stantion: " << istal << " portion number: " << ip << " CATrackFinder: Warning: Too many Triplets created in portion" << endl;


     // cout<<u_front3.size()<<" scalar"<<endl;

          /// Add the right hits to parameters estimation.
//        f31(  // input
//            n3_V,
//            star,
//            u_front3, u_back3, z_pos3,
//           
//              // output
//            T_3
//           );

                  /// refit
      //         f32(  n3, istal, _RealIHit,  TripletsArray,         T_3,         hitsl_3, hitsm_3, hitsr_3, MHits, LHits, RHits );


#ifdef TRIP_PERFORMANCE
        THitI* RealIHitL = &((*RealIHitP)[StsHitsUnusedStartIndex[istal]]);
        THitI* RealIHitM = &((*RealIHitP)[StsHitsUnusedStartIndex[istam]]);
        THitI* RealIHitR = &((*RealIHitP)[StsHitsUnusedStartIndex[istar]]);
        for (Tindex i = 0; i < n3; ++i){
          Tindex i_4 = i%4;
          Tindex i_V = i/4;
          THitI iHits[3] = {
            RealIHitL[hitsl_3[i]],
            RealIHitM[hitsm_3[i]],
            RealIHitR[hitsr_3[i]]
          };
#ifdef PULLS
          if ( fL1Eff_triplets->AddOne(iHits) )
            fL1Pulls->AddOne(T_3[i_V], i_4, iHits[2]);
#else
          fL1Eff_triplets->AddOne(iHits);
#endif
     if (T_3[i_V].chi2[i_4] < TRIPLET_CHI2_CUT) 
            fL1Eff_triplets2->AddOne(iHits);
        }
#endif // TRIP_PERFORMANCE

          /// Fill Triplets.

//        f4( // input
//          n3, istal, istam, istar,
//          T_3,
//          hitsl_3, hitsm_3, hitsr_3,
//            // output
//          nstaltriplets,
//         
//            ip_cur 
//         
//          );
  }
}
  

      
// hitCheck::hitCheck()
// {
//          omp_init_lock(&Occupied); 
//          trackCandidateIndex = -1;
//          UsedByTrack=0;
//          Chi2Track = 100000000;
//          Length = 0;
//          ista = 1000;
// 
// }
// hitCheck::~hitCheck()
// {
// omp_destroy_lock(&Occupied);
// }
    ///**********************************************************************************************************************
    ///**********************************************************************************************************************
    ///**********************************************************************************************************************
    ///*                                                                                                                    *
    ///*                                                                                                                    *
    ///*                                                                                                                    *
    ///*                                          CATrackFinder procedure                                                   *
    ///*                                                                                                                    *
    ///*                                                                                                                    *
    ///*                                                                                                                    *
    ///**********************************************************************************************************************
    ///**********************************************************************************************************************
    ///**********************************************************************************************************************




void L1Algo::CATrackFinder()
{
  omp_set_num_threads(fNThreads);

#ifdef PULLS
  static L1AlgoPulls *l1Pulls_ = new L1AlgoPulls();
  fL1Pulls = l1Pulls_;
  fL1Pulls->Init();
#endif
#ifdef TRIP_PERFORMANCE
  static L1AlgoEfficiencyPerformance<3> *l1Eff_triplets_ = new L1AlgoEfficiencyPerformance<3>();
  fL1Eff_triplets = l1Eff_triplets_;
  fL1Eff_triplets->Init();
  static L1AlgoEfficiencyPerformance<3> *l1Eff_triplets2_ = new L1AlgoEfficiencyPerformance<3>();
  fL1Eff_triplets2 = l1Eff_triplets2_;
  fL1Eff_triplets2->Init();
#endif
#ifdef DOUB_PERFORMANCE
  static L1AlgoEfficiencyPerformance<2> *l1Eff_doublets_ = new L1AlgoEfficiencyPerformance<2>();
  fL1Eff_doublets = l1Eff_doublets_;
  fL1Eff_doublets->Init();
#endif
  
#ifdef DRAW
  static L1AlgoDraw draw;
  draw.InitL1Draw(this);
#endif

  cout.precision(6);

  TStopwatch c_time;   // for performance time
#if defined(XXX) || defined(COUNTERS)
  static unsigned int stat_N = 0; // number of events
  stat_N++;
#endif
  
#ifdef XXX
  TStopwatch c_timerG; // global
  TStopwatch c_timerI; // for iterations
  
  L1CATFIterTimerInfo gti; // global
  gti.Add("init  ");
  gti.Add("iterts");
  gti.Add("merge ");
  
  L1CATFTimerInfo ti;
  ti.SetNIter( fNFindIterations ); // for iterations
  ti.Add("init  ");
  // ti.Add("dblte1");
  // ti.Add("dblte2");
  ti.Add("tripl1");

 
  ti.Add("tracks");
   ti.Add("table");
    ti.Add("save");
     ti.Add("delete");
     ti.Add("copy");
  ti.Add("finish");
  
  static L1CATFIterTimerInfo stat_gti = gti;
  static L1CATFTimerInfo stat_ti = ti;
  

#endif

#ifdef COUNTERS
  static Tindex stat_nStartHits = 0;
  static Tindex stat_nHits[fNFindIterations] = {0};
  
  static Tindex stat_nSinglets[fNFindIterations] = {0};
//  static Tindex stat_nDoublets[fNFindIterations] = {0};
  static Tindex stat_nTriplets[fNFindIterations] = {0};
  
  static Tindex stat_nLevels[MaxNStations-2][fNFindIterations] = {{0}};
  static Tindex stat_nCalls[fNFindIterations] = {0}; // n calls of CAFindTrack
  static Tindex stat_nTrCandidates[fNFindIterations] = {0};
#endif
  
//    vector<hitCheck> hitToBestTrackB[fNThreads];
//    vector<hitCheck> hitToBestTrackF[fNThreads];
    
   
   
//     for (int i=0; i<fNThreads; i++)   
//     {  
//   //  hitToBestTrackB.resize(0);
//     
//     hitToBestTrackB[i].resize((*vSFlagB).size());
//    // hitToBestTrackF.resize(0);
//     hitToBestTrackF[i].resize((*vSFlag).size());
//     }
    
 
/*    

    RealIHit_v.resize(0);
    RealIHit_v.resize(2200000);
    RealIHit_v_buf.resize(0);
    RealIHit_v_buf.resize(2200000);
    RealIHit_v_buf2.resize(0);
    RealIHit_v_buf2.resize(2200000);
    
    
    
    vTracks.resize(0);
    vRecoHits.resize(0);
    vStsDontUsedHits_A.resize(0);
    vStsDontUsedHits_B.resize(0);
    vStsDontUsedHits_Buf.resize(0);
    vStsDontUsedHitsxy_A.resize(0);
    vStsDontUsedHitsxy_buf.resize(0);
    vStsDontUsedHitsxy_B.resize(0);
    
    vTracks.resize(100000);
    vRecoHits.resize(550000);
    vStsDontUsedHits_A.resize(2200000);
    vStsDontUsedHits_B.resize(2200000);
    vStsDontUsedHits_Buf.resize(2200000);
    vStsDontUsedHitsxy_A.resize(2200000);
    vStsDontUsedHitsxy_buf.resize(2200000);
    vStsDontUsedHitsxy_B.resize(2200000);

          for( unsigned int iV=0; iV<(*vSFlag).size(); ++iV)
              SetFUnUsed(const_cast<unsigned char&>((*vSFlag)[iV]));
          for( unsigned int iV=0; iV<(*vSFlagB).size(); ++iV)
              SetFUnUsed(const_cast<unsigned char&>((*vSFlagB)[iV]));
    
    
    for (int i=0; i<2200000; i++) TripForHit[i].resize(0);
    
    
    


    for (int i=0; i<8; i++)
        for (int j=0; j<5000; j++)
        {
                       for (int k=0; k<5000; k++) TripletsLocalReal[i][j][k].neighbours.resize(0);
            
        }




    n_g1.resize(5000, Portion);
    
      
    for (int i=0; i<MaxNStations; i++)
        vGrid[i].Initial1(fNThreads);*/



 vector <L1Branch*> candidatesAll;
 
 vector <L1Branch*> *candidatesAllP;

 candidatesAllP=&candidatesAll;
 
  candidatesAll.resize(30000, 0);

  
//         vector < vector <L1Branch*> > Common[fNThreads];
//         vector <short int> CommonIndex[fNThreads];
//       
//         for (int i=0; i<fNThreads; ++i){
// 
//         Common[i].resize(8000);
//         CommonIndex[i].resize(8000, 0);
//         
//     
//             
//             
//             for ( Tindex iC = 0; iC < 8000; ++iC ){
//                      Common[i][iC].resize(300, 0);
//                 }
//         
//         }

  #ifdef XXX
    
    c_time.Start();
    c_timerG.Start();
    
#endif  
    
/*  
    TripForHit.resize(2000000);
    hitToBestTrackB.resize(2000000);
    hitToBestTrackF.resize(1200000);
    hitToBestTrackB1.resize(2000000);
    hitToBestTrackF1.resize(1200000);

      vTracks.resize(90000); // reconstructed tracks
    vRecoHits.resize(500000);// packed hits of reconstructed tracks
   vStsDontUsedHits_A.resize(2000000);
  vStsDontUsedHits_B.resize(2000000);
  vStsDontUsedHits_Buf.resize(2000000);
  vStsDontUsedHitsxy_A.resize(2000000);
  vStsDontUsedHitsxy_buf.resize(2000000);
 vStsDontUsedHitsxy_B.resize(2000000);
 RealIHit_v.resize(2000000);
    RealIHit_v_buf.resize(2000000);
    RealIHit_v_buf2.resize(2000000);
    */
    
    
    
    
    
      
  RealIHitP=&RealIHit_v;
  RealIHitPBuf=&RealIHit_v_buf;
  vStsHitsUnused = &vStsDontUsedHits_B;      /// array of hits used on current iteration
  std::vector< L1StsHit > *vStsHitsUnused_buf = &vStsDontUsedHits_A; /// buffer for copy

  vStsHitPointsUnused = &vStsDontUsedHitsxy_B;/// array of info for hits used on current iteration
  std::vector< L1HitPoint > *vStsHitPointsUnused_buf = &vStsDontUsedHitsxy_A;
  
     NHitsIsecAll=0;
    NTracksIsecAll=0;
 
  int nDontUsedHits = 0;
   
 // #pragma omp parallel for  reduction(+:nDontUsedHits)
  for(int ista = 0; ista < NStations; ++ista){
    nDontUsedHits += (StsHitsStopIndex[ista]-StsHitsStartIndex[ista]);
     StsHitsUnusedStartIndex[ista] = StsHitsStartIndex[ista];
       StsHitsUnusedStopIndex[ista]  =  StsHitsStopIndex[ista];
  }

  

  
  
#ifndef L1_NO_ASSERT
//     for (int istal = NStations - 1; istal >= 0; istal--){
//       L1_ASSERT( StsHitsStopIndex[istal] <= (*vStsHits).size()             , StsHitsStopIndex[istal] << " <= " << (*vStsHits).size());
//       L1_ASSERT( StsHitsUnusedStopIndex[istal] <= (*vStsHitsUnused).size(), StsHitsUnusedStopIndex[istal] << " <= " << (*vStsHitsUnused).size());
//     }
#endif // L1_NO_ASSERT
  
    float yStep = 0.55/sqrt( nDontUsedHits ); // empirics. 0.01*sqrt(2374) ~= 0.5
    if (yStep > 0.3) yStep = 0.3;
    const float xStep = yStep*8;
    
//    float yStep = 0.5/sqrt( nDontUsedHits ); // empirics. 0.01*sqrt(2374) ~= 0.5
//    if (yStep > 0.3) yStep = 0.3;
//    const float xStep = yStep*3;



    
     #pragma omp parallel
        for( int iS = 0; iS < NStations; ++iS ) {     /// Grid is created for each station with the same step: xStep,yStep
      L1Grid &grid = vGrid[iS];
      grid.CreatePar0(-1,1,-0.6,0.6,xStep,yStep);
            
            
    
    
    }

    
    vStsHitsUnused = &vStsDontUsedHits_Buf;
    

    
    // #pragma omp parallel
//    for( int iS = 0; iS < NStations; ++iS ) {     /// Grid is created for each station with the same step: xStep,yStep
//        L1Grid &grid = vGrid[iS];
//        grid.CreateParSIMD(&((*vStsHitPointsUnused)[StsHitsUnusedStartIndex[iS]]),StsHitsUnusedStopIndex[iS]-StsHitsUnusedStartIndex[iS], &vStsDontUsedHitsxy_buf, &vStsDontUsedHits_Buf, &((*vStsHits)[StsHitsUnusedStartIndex[iS]]), &RealIHit_v, &RealIHit_v_buf, iS, *this, StsHitsUnusedStartIndex[iS], vStsStrips, vStsStripsB, vStsZPos);
//        grid.HitsSortSIMD(&(vStsDontUsedHitsxy_buf[StsHitsUnusedStartIndex[iS]]), &(vStsDontUsedHits_Buf[StsHitsUnusedStartIndex[iS]]), &((*vStsHits)[StsHitsUnusedStartIndex[iS]]), &(RealIHit_v[StsHitsUnusedStartIndex[iS]]), &(RealIHit_v_buf[StsHitsUnusedStartIndex[iS]]), &((*vStsHitPointsUnused)[StsHitsUnusedStartIndex[iS]]), StsHitsUnusedStartIndex[iS], StsHitsUnusedStopIndex[iS]-StsHitsUnusedStartIndex[iS], iS, *this, &vStsHitsUnusedV[iS], vStsZPos, vStsStrips, vStsStripsB);
//        
//        
//    }
    
    
    
  
   // #pragma omp parallel 
    for( int iS = 0; iS < NStations; ++iS ) {     /// Grid is created for each station with the same step: xStep,yStep
      L1Grid &grid = vGrid[iS];
      grid.CreatePar(&((*vStsHitPointsUnused)[StsHitsUnusedStartIndex[iS]]),StsHitsUnusedStopIndex[iS]-StsHitsUnusedStartIndex[iS], &vStsDontUsedHitsxy_buf, &vStsDontUsedHits_Buf, &((*vStsHits)[StsHitsUnusedStartIndex[iS]]), &RealIHit_v, &RealIHit_v_buf, iS, *this, StsHitsUnusedStartIndex[iS]);
         grid.HitsSort(&(vStsDontUsedHitsxy_buf[StsHitsUnusedStartIndex[iS]]), &(vStsDontUsedHits_Buf[StsHitsUnusedStartIndex[iS]]), &((*vStsHits)[StsHitsUnusedStartIndex[iS]]), &(RealIHit_v[StsHitsUnusedStartIndex[iS]]), &(RealIHit_v_buf[StsHitsUnusedStartIndex[iS]]), &((*vStsHitPointsUnused)[StsHitsUnusedStartIndex[iS]]), StsHitsUnusedStartIndex[iS], StsHitsUnusedStopIndex[iS]-StsHitsUnusedStartIndex[iS], iS, *this, &vStsHitsUnusedV[iS], vStsZPos, vStsStrips, vStsStripsB);
    
    
    }
    
    

    


    // SIMD HitPoints

    for(int ista = 0; ista < NStations; ++ista){
        
        L1Station &sta = vStations[ista];
        
       #pragma omp parallel for  schedule(dynamic, 9) 
   
      for(THitI ih = 0; ih < StsHitsStopIndex[ista]-StsHitsStartIndex[ista]; ih+=fvecLen)
      {
        const fvec &f = reinterpret_cast<const fvec&>(vStsHitsUnusedV[ista].f[ih]);
        const fvec &b = reinterpret_cast<const fvec&>(vStsHitsUnusedV[ista].b[ih]);
        const fvec &u = reinterpret_cast<const fvec&>(vStsHitsUnusedV[ista].u[ih]);
        const fvec &v = reinterpret_cast<const fvec&>(vStsHitsUnusedV[ista].v[ih]);
        const fvec &z = reinterpret_cast<const fvec&>(vStsHitsUnusedV[ista].z[ih]);
        int n = ih/fvecLen;

        
        const fvec x = sta.xInfo.sin_phi*u + sta.xInfo.cos_phi*v;
        const fvec y = sta.yInfo.cos_phi*u + sta.yInfo.sin_phi*v;

        vStsHitPointsUnusedV[ista].SetL1HitPointV(x,y,z,v,u, n);

      }
    }

    
    
 

      {

           for(int ista = 0; ista < NStations; ++ista){
                #pragma omp parallel for
      for(THitI ih = StsHitsStartIndex[ista]; ih < StsHitsStopIndex[ista]; ++ih)
      {  vStsDontUsedHitsxy_B[ih] = CreateHitPoint(vStsDontUsedHits_Buf[ih],ista);
          
    }
      }
  }
  
     

//            for(int ista = 0; ista < NStations; ++ista){
//                 
//       for(THitI ih = StsHitsStartIndex[ista]; ih < StsHitsStopIndex[ista]; ++ih)
//       { // if (vStsDontUsedHitsxy_B[ih].Y() != vStsHitPointsUnusedV[ista].Y((ih-StsHitsStartIndex[ista])/fvecLen)[(ih-StsHitsStartIndex[ista])%fvecLen]) 
// cout<<vStsDontUsedHitsxy_B[ih].Y()<<" Y "<<vStsHitPointsUnusedV[ista].Y((ih-StsHitsStartIndex[ista])/fvecLen)[(ih-StsHitsStartIndex[ista])%fvecLen]<<endl;
// cout<<vStsDontUsedHitsxy_B[ih].X()<<" X "<<vStsHitPointsUnusedV[ista].X((ih-StsHitsStartIndex[ista])/fvecLen)[(ih-StsHitsStartIndex[ista])%fvecLen]<<endl;
// cout<<vStsDontUsedHitsxy_B[ih].U()<<" U "<<vStsHitPointsUnusedV[ista].U((ih-StsHitsStartIndex[ista])/fvecLen)[(ih-StsHitsStartIndex[ista])%fvecLen]<<endl;
// cout<<vStsDontUsedHitsxy_B[ih].V()<<" V "<<vStsHitPointsUnusedV[ista].V((ih-StsHitsStartIndex[ista])/fvecLen)[(ih-StsHitsStartIndex[ista])%fvecLen]<<endl;
//           
//     }
//       }
  
  

  
//(const fvec &f, const fvec &b, const fvec &z, char ista, int n, L1HitPointV Point)

      

  

    

    // #pragma omp parallel 
     //   for( int iS = 0; iS < NStations; ++iS ) {     /// Grid is created for each station with the same step: xStep,yStep
     // L1Grid &grid = vGrid[iS];
    //  grid.CreatePar(-1,1,-0.6,0.6,xStep,yStep, &((*vStsHitPointsUnused)[StsHitsUnusedStartIndex[iS]]),StsHitsUnusedStopIndex[iS]-StsHitsUnusedStartIndex[iS], &vStsDontUsedHitsxy_buf, &vStsDontUsedHits_Buf, vStsHitsUnused, &RealIHit_v, &);
      //grid.HitsSort(&vStsDontUsedHitsxy_buf, &vStsDontUsedHits_Buf, vStsHitsUnused, &RealIHit_v, &RealIHit_v_buf, vStsHitPointsUnused, StsHitsUnusedStartIndex[iS], StsHitsUnusedStopIndex[iS]-StsHitsUnusedStartIndex[iS]);
      
//    }

    

  //  sh.Fill(vStsHitsUnused, vStsHitPointsUnused, &RealIHit_v, vGrid, NStations, StsHitsUnusedStartIndex, StsHitsUnusedStopIndex);
    
    

  //  sh.Sort(nDontUsedHits, &vStsDontUsedHitsxy_buf, &vStsDontUsedHits_Buf);



#ifdef COUNTERS
  cout << " Begin " << endl;
  stat_nStartHits += nDontUsedHits;
  cout << " NStartHits = " << stat_nStartHits/stat_N << endl;
#endif 

    
    #ifdef XXX
    
    c_timerG.Stop();
    gti["init  "] = c_timerG;
    c_timerG.Start(1);
    
#endif
    
  for (isec = 0; isec < fNFindIterations; ++isec){ // all finder

  /// isec - number of current iterations, fNFindIterations - number of all iterations 
    
#ifdef COUNTERS
  Tindex nSinglets = 0;
#endif
    
#ifdef XXX
     TStopwatch c_timer;
     TStopwatch c_time;
      c_timer.Start();
#endif

  
{   
     
  //if (isec!=0)
//    for( int iS = 0; iS < NStations; ++iS ) {
//        
//    L1Grid &grid = vGrid[iS]; // reference to array element iS
//    vGrid[iS].FillPar(&((*vStsHitPointsUnused)[StsHitsUnusedStartIndex[iS]]),StsHitsUnusedStopIndex[iS]-StsHitsUnusedStartIndex[iS]); /// input: 1st hit on station, number of hits on station
//         
//    
//      }
    

       c_timer.Stop();
    ti[isec]["init  "] = c_timer;
    c_timer.Start(1);

 
 
  #pragma omp  task
      {
       
     // --- SET PARAMETERS FOR THE ITERATION ---

       // first station used in the triplets finding
     FIRSTCASTATION = 0;
     // if ( (isec == kAllSecIter) || (isec == kAllSecJumpIter) )
     //   FIRSTCASTATION = 2;
      DOUBLET_CHI2_CUT = 2.706; // prob = 0.1 
      switch ( isec ) { 
       case kFastPrimIter: 
       TRIPLET_CHI2_CUT = 7.815;// prob = 0.05 
       break; 
       case kAllPrimIter: 
       TRIPLET_CHI2_CUT = 7.815; // prob = 0.05 
       break; 
       case kAllPrimJumpIter: 
       TRIPLET_CHI2_CUT = 6.252; // prob = 0.1 
       break; 
       case kAllSecIter: 
       TRIPLET_CHI2_CUT = 2.706; // prob = 0.1 
       break; 
  
      }

    Pick_gather = 3.0; /// coefficient for size of region for attach new hits to the created track
    if ( (isec == kAllPrimIter) ||
         (isec == kAllPrimJumpIter) ||
         (isec == kAllSecIter) ||
         (isec == kAllSecJumpIter) )
      Pick_gather = 4.0;

    PickNeighbour = 1.0; // (PickNeighbour < dp/dp_error)  =>  triplets are neighbours
    // if ( (isec == kFastPrimIter) )
    //   PickNeighbour = 0.5; // TODO understand why works with 0.2

    MaxInvMom = 1.0/0.5;                     /// max considered q/p
    if ( (isec == kAllPrimIter) || (isec == kAllPrimJumpIter) || (isec == kAllSecIter) || (isec == kAllSecJumpIter) ) MaxInvMom =  1.0/0.1;

    MaxSlope = 1.1;
    if ( // (isec == kAllPrimIter) || (isec == kAllPrimJumpIter) ||
         (isec == kAllSecIter) || (isec == kAllSecJumpIter) ) MaxSlope =  1.5;
    
      /// define the target position
    targX = 0; targY = 0; targZ = 0;      ///  the target is plased at (0,0,0)
    
    float SigmaTargetX = 0, SigmaTargetY = 0; /// target constraint [cm]
    if ( (isec == kFastPrimIter) || (isec == kFastPrimIter2) || (isec == kFastPrimJumpIter) ||
         (isec == kAllPrimIter) || (isec == kAllPrimJumpIter) ){ // target
      targB = vtxFieldValue;
      if ( (isec == kFastPrimIter) || (isec == kAllPrimIter) )
        SigmaTargetX = SigmaTargetY = 0.01; // target
      else
        SigmaTargetX = SigmaTargetY = 0.1;
    }
    if ( (isec == kAllSecIter) || (isec == kAllSecJumpIter) ) { //to do: use outer radius of the 1st station as a constraint
      L1Station &st = vStations[0];
      SigmaTargetX = SigmaTargetY = 3;//st.Rmax[0];
      targZ = 0.;//-1.;
      st.fieldSlice.GetFieldValue( 0.f, 0.f, targB );
    }

    TargetXYInfo.C00 = SigmaTargetX * SigmaTargetX;
    TargetXYInfo.C10 = 0;
    TargetXYInfo.C11 = SigmaTargetY * SigmaTargetY;

      // Set correction in order to take into account overlaping and iff z.
      // The reason is that low momentum tracks are too curved and goes not from target direction. That's why sort by hit_y/hit_z is not work idealy
      // If sort by y then it is max diff between same station's modules (~0.4cm)
    // ya podumau ob etom
    MaxDZ = 0;
    if (  (isec == kAllPrimIter) || (isec == kAllPrimJumpIter) ||
          (isec == kAllSecIter ) || (isec == kAllSecJumpIter ) ) MaxDZ = 0.1;


    
    if (NStations > MaxNStations) cout << " CATrackFinder: Error: Too many Stantions" << endl;
      }
 
   

    
   

#ifndef L1_NO_ASSERT

    for (int istal = NStations - 1; istal >= 0; istal--){
      L1_ASSERT( StsHitsUnusedStopIndex[istal] >= StsHitsUnusedStartIndex[istal], StsHitsUnusedStopIndex[istal] << " >= " << StsHitsUnusedStartIndex[istal]);
      L1_ASSERT( StsHitsUnusedStopIndex[istal] <= (*vStsHitsUnused).size(), StsHitsUnusedStopIndex[istal] << " <= " << (*vStsHitsUnused).size());
    }
    
//      cout<<" L1_NO_ASSERT"<<endl;
#endif // L1_NO_ASSERT
    Tindex nPortions = 0; ///  number of portions

   int StsHitsUnusedStopIndexEnd[10];
   int StsHitsUnusedStartIndexEnd[10];
   
   for (int i=0; i<NStations;i++)
   {
     StsHitsUnusedStopIndexEnd[i]=StsHitsUnusedStopIndex[i];
    StsHitsUnusedStartIndexEnd[i]=StsHitsUnusedStartIndex[i];
   }
   
   
           /// possible left hits of triplets are splited in portions of 16 (4 SIMDs) to use memory faster
      portionStopIndex[NStations-1] = 0;
      Tindex ip = 0;  //ip - index of curent portion
/// possible left hits of triplets are splited in portions of 16 (4 SIMDs) to use memory faster     
      for (int istal = NStations-2; istal >= FIRSTCASTATION; istal--){  //start downstream chambers
        NHits_l [istal]= StsHitsUnusedStopIndex[istal] - StsHitsUnusedStartIndex[istal];

        NHits_l_P [istal]= NHits_l [istal]/Portion;
       
       
        ip+=NHits_l_P[istal];
     

        n_g1[ip]=(NHits_l[istal] - NHits_l_P[istal]*Portion);
        ip++;
        portionStopIndex[istal] = ip;
       }// lstations
      nPortions = ip;

    
#ifdef COUNTERS
  stat_nSinglets[isec] += nSinglets;
#endif
}


  
    
      ///   stage for triplets creation
      


     const Tindex vPortion = Portion/fvecLen;
     L1TrackPar T_1[vPortion];
     L1TrackParV T_1V;
     L1TrackParV T_1VBuf;
     L1TrackParV T_1VBuf2;
     T_1V.resize(vPortion*fvecLen);
     T_1VBuf.resize(5000);
    // T_1VBuf2.resize(5000);
     L1FieldRegion fld_1[vPortion];
     L1FieldRegionV fld_1V;
     fld_1V.resize(vPortion*fvecLen);
     THitI hitsl_1[Portion];
//      L1TrackPar TG_1[vPortion];
//      L1FieldRegion fldG_1[vPortion];
//      THitI hitslG_1[Portion];
     
//      vector<bool> lmDuplets[MaxNStations]; // is exist a doublet started from indexed by left hit
//      vector<bool> lmDupletsG[MaxNStations]; // is exist a doublet started from indexed by left hit
      
       c_timer.Start();
      
      
      for (int istal = NStations-2; istal >= FIRSTCASTATION; istal--){//  //start downstream chambers
          
          const int num_portions= portionStopIndex[istal]-portionStopIndex[istal+1];
          numPortions[istal]= num_portions;
          //  TripletsLocal[istal] = new vector<L1Triplet*> [num_portions];
//#pragma omp parallel for firstprivate(T_1, T_1V, fld_1, fld_1V, hitsl_1) schedule(dynamic)
      for(Tindex ip = portionStopIndex[istal+1]; ip < portionStopIndex[istal]; ++ip ){
          
         
          TripNPortionV[istal][ip-portionStopIndex[istal+1]]=0;
          
        static int k = 0;
          nsL1::vector <int>::TSimd MHits;
          nsL1::vector <int>::TSimd LHits;
          nsL1::vector <int>::TSimd RHits;
          nsL1::vector <int>::TSimd LHits3;
          
          MHits.reserve(MaxPortionDoublets);
          LHits.reserve(MaxPortionDoublets);
          RHits.reserve(MaxPortionDoublets);
          LHits3.reserve(MaxPortionDoublets);
           Tindex n_2;
          
          L1TrackParV TripletsArray;
          
          
          nsL1::vector <L1TrackPar>::TSimd  TripletsArrayV;
          
          
          
          
          
          
          DupletsStaPortSIMD(
                             istal, istal + 1,
                             ip,
                             n_g1, portionStopIndex,
                             
                             // output
                             
                             &T_1V,
                             &T_1VBuf,
                             &T_1VBuf2,
                             
                             &fld_1V,
                             
                             
                             //            lmDuplets[istal],
                             
                             
                             n_2,
                             
                             
                             MHits,
                             LHits,
                             LHits3
                             );
          Tindex nstaltriplets=0;
          
           TripletsStaPortSIMD(  // input
                               istal, istal + 1,  istal + 2,
                               nstaltriplets,
                              
                               &T_1V,
                               &T_1VBuf2,
                               &T_1VBuf,

                              
                               &fld_1V,
                               
                               
                               TripletsArray,
                               
                               n_2,
                               
                               
                               //            lmDuplets[istal+1],
                               // output
                               
                               ip-portionStopIndex[istal+1],
                               MHits,
                               LHits,
                               RHits,
                               TripletsArrayV,
                                                              LHits3
                               );
          

          
          k++;
// if (k>160) break;
          
          
          
          
          
      }// 
  }

      
#ifdef XXX
      c_timer.Stop();
      ti[isec]["copy"] = c_timer;
      c_timer.Start(1);
#endif
      
      
     for (int istal = NStations-2; istal >= FIRSTCASTATION; istal--){//  //start downstream chambers

       const int num_portions= portionStopIndex[istal]-portionStopIndex[istal+1];   
       numPortions[istal]= num_portions; 
     //  TripletsLocal[istal] = new vector<L1Triplet*> [num_portions];
    
        
       
     //  #pragma omp parallel for firstprivate(T_1, T_1V, fld_1, fld_1V, hitsl_1) schedule(dynamic)
       for(Tindex ip = portionStopIndex[istal+1]; ip < portionStopIndex[istal]; ++ip ){
         
           TripNPortion[istal][ip-portionStopIndex[istal+1]]=0;
           //TripNPortionV[istal][ip-portionStopIndex[istal+1]]=0;

         // nsL1::vector<L1TrackPar>::TSimd T_1; // estimation of parameters
         // nsL1::vector<L1FieldRegion>::TSimd fld_1; // magnetic field
         // vector<THitI> hitsl_1; // left hits indexed by number of singlets in portion(i1)
         vector<THitI> hitsm_2; /// middle hits indexed by number of doublets in portion(i2)
         vector<THitI> i1_2; /// index in portion of singlets(i1) indexed by index in portion of doublets(i2)
         Tindex n_2; /// number of doublets in portion

         hitsm_2.reserve(MaxPortionDoublets);
         i1_2.reserve(MaxPortionDoublets);
//         nsL1::vector <int>::TSimd MHits;
//         nsL1::vector <int>::TSimd LHits;
//         nsL1::vector <int>::TSimd RHits;
//           nsL1::vector <int>::TSimd LHits3;
//           
//           MHits.reserve(MaxPortionDoublets);
//           LHits.reserve(MaxPortionDoublets);
//           RHits.reserve(MaxPortionDoublets);
//           LHits3.reserve(MaxPortionDoublets);
//           
//           
//            L1TrackParV TripletsArray;
//           
//           
//           nsL1::vector <L1TrackPar>::TSimd  TripletsArrayV;
           


static int k1=0;
           
           
//           DupletsStaPortSIMD(
//                          istal, istal + 1,
//                          ip,
//                          n_g1, portionStopIndex,
//                          
//                          // output
//                          
//                          &T_1V,
//                          
//                          &fld_1V,
//                          
//                          
//                          //            lmDuplets[istal],
//                          
//                          
//                          n_2,
//                          
//                        
//                          MHits,
//                          LHits,
//                          LHits3
//                          );
            Tindex nstaltriplets=0;
           
//           TripletsStaPortSIMD(  // input
//                           istal, istal + 1,  istal + 2,
//                           nstaltriplets,
//                           T_1,
//                           &T_1V,
//                           fld_1,
//                           &fld_1V,
//                           hitsl_1,
//                           TripletsArray,
//                           
//                           n_2,
//                           i1_2,
//                           hitsm_2,
//                           
//                           //            lmDuplets[istal+1],
//                           // output
//                           
//                           ip-portionStopIndex[istal+1],
//                           MHits,
//                           LHits,
//                           RHits,
//                           TripletsArrayV,
//                           LHits3
//                           );
      
           
           DupletsStaPort(
                          istal, istal + 1,
                          ip,
                          n_g1, portionStopIndex,
                          
                          // output
                          T_1,
                         
                          fld_1,
                        
                          hitsl_1,
                          
                          //            lmDuplets[istal],
                          
                          
                          n_2,
                          i1_2,
                          hitsm_2
                         
                          );
           

      
           nstaltriplets=0;
  
        TripletsStaPort(  // input
            istal, istal + 1,  istal + 2,
            nstaltriplets,
            T_1,
            fld_1,
            hitsl_1,
   
            n_2,
            i1_2,
            hitsm_2,
   
 //            lmDuplets[istal+1],
               // output
           
                       ip-portionStopIndex[istal+1]
             );
           
           
           
          k1++;
// if (k1>160) break;
           
           
           
       }// 
     }

     
#ifdef XXX
c_timer.Stop();
ti[isec]["tripl1"] = c_timer;
c_timer.Start(1);
#endif
  
      

      

    ///====================================================================
    ///=                                                                  =
    ///=        Collect track candidates. CREATE TRACKS                   =
    ///=                                                                  =
    ///====================================================================

// #ifdef XXX
//     cout<<"---- Collect track candidates. ----"<<endl;
// #endif

 
 

      int size=0;
      int min_level = 1; // min level for start triplet. So min track length = min_level+3.
      if (isec == kAllPrimJumpIter) min_level = 1;
      if ( (isec == kAllSecIter) || (isec == kAllSecJumpIter) ) min_level = 2; // only the long low momentum tracks
      //    if (isec == -1) min_level = NStations-3 - 3; //only the longest tracks

 
   
    int  numberCandidateThread2 [20]={0};
   int  numberCandidateThreadBuf [20]={0};
             int SavedCand[20]={0};    
      int SavedHits[20]={0}; 
    static int count = 0;  
        int UsedF[20]={0};
        int UsedB[20]={0};
    int  numberCandidateThread [20]={0};

c_timer.Start();

      
      // collect consequtive: the longest tracks, shorter, more shorter and so on
    for (int ilev = NStations-3; ilev >= min_level; ilev--){
// choose length in triplets number - ilev - the maximum possible triplet level among all triplets in the searched candidate
      //  cout<<ilev<<endl;
      //array to store candidates
       TStopwatch Time;

     Time.Start(1);
    
//    omp_set_num_threads(fNThreads);


        //  how many levels to check
      int nlevel = (NStations-2)-ilev+1;
    
      
      const unsigned char min_best_l = (ilev > min_level) ? ilev + 2 : min_level + 3; // loose maximum one hit and find all hits for min_level
      // value to check that candidate has at least min_best_l hits (corresponds to ilev or ilev - 1 hit)




      L1Branch curr_tr;
      curr_tr.StsHits.resize(8);
       L1Branch new_tr[8];
         //#pragma omp parallel for collapse(2) schedule(dynamic)  firstprivate(curr_tr, new_tr)
         

     

    
         for (int i=0; i<fNThreads; ++i){
       
         UsedF[i]=0;
         UsedB[i]=0;
         numberCandidateThread [i]=0;
      }
  int check1=0;
  int check2=0;

      for( int istaF = FIRSTCASTATION; istaF <= NStations-3-ilev; ++istaF ){
        

          
            #pragma omp parallel for schedule(dynamic)  firstprivate(curr_tr, new_tr)
             for( Tindex ip=0; ip<numPortions[istaF]; ++ip ){
                 
                
         int thread_num = omp_get_thread_num();
                for( Tindex itrip=0; itrip<TripNPortion[istaF][ip]; ++itrip ){

                
      

   
        L1Triplet &first_trip = *(TripletsLocalRealAll[istaF][ip][itrip]);
                    
              
          { // ghost supression !!!
     //       if ( first_trip->GetLevel() == 0 ) continue; // ghost suppression // find track with 3 hits only if it was created from a chain of triplets, but not from only one triplet
            if ( first_trip.GetLevel() < ilev ) break; // try only triplets, which can start track with ilev+3 length. w\o it have more ghosts, but efficiency either
            if ( (ilev == 0) &&
                 (GetFStation((*vSFlag)[(*vStsHitsUnused)[first_trip.GetLHit()].f]) != 0) ) continue;  // ghost supression // collect only MAPS tracks-triplets
          }
         
          if( GetFUsed( (*vSFlag)[(*vStsHitsUnused)[first_trip.GetLHit()].f] | (*vSFlagB)[(*vStsHitsUnused)[first_trip.GetLHit()].b] ) ) continue;
 
          curr_tr.Momentum = 0.f;
          curr_tr.chi2 = 0.f;
          curr_tr.Lengtha = 0;
          curr_tr.ista = 0;
          //curr_tr.StsHits.resize(8);
         
      
          //L1Branch* curr_tr=&CandidatesTrack[thread_num][ numberCandidateThread [thread_num]];
               (curr_tr).StsHits[0]=((*RealIHitP)[first_trip.GetLHit()]);
               
              
                   
            //  (curr_tr).StsHits.push_back((*RealIHitP)[first_trip->GetLHit()]);
              (curr_tr).NHits=1;
              
      
          unsigned char curr_L = 1;
          fscal curr_chi2 = first_trip.GetChi2();

          L1Branch best_tr = (curr_tr);
           
          fscal best_chi2 = curr_chi2;
          unsigned char best_L = curr_L;
          int NCalls = 1;
             int NHitTrack=1;
  
          CAFindTrack(istaF, best_tr, best_L, best_chi2, &first_trip, (curr_tr), curr_L, curr_chi2, min_best_l, NCalls, new_tr); /// reqursive func to build a tree of possible track-candidates and choose the best
     
            if ( best_L < min_best_l ) continue;
          
          int ndf = best_L*2-5;
          best_chi2 = best_chi2/ndf; //normalize

//             if( fGhostSuppression ){
//               if( best_L == 3 ){
//              
//                 if( ( (isec == kAllSecIter) || (isec == kAllSecJumpIter) ) && (istaF != 0) ) continue; 
//                 if( (isec != kAllSecIter) && (isec != kAllSecJumpIter) && (best_chi2 > 5.0) ) continue;
//               }
//             }
 
          best_tr.Set( istaF, best_L, best_chi2, first_trip.GetQpOrig());
         
            
          best_tr.CandIndex=numberCandidateThread [thread_num]+thread_num*100000;
          best_tr.StsHits.resize(best_tr.NHits);
 
          CandidatesTrack[thread_num][ numberCandidateThread [thread_num]]=best_tr;
          CandidatesTrack_buf2[thread_num][ numberCandidateThread [thread_num]]=&CandidatesTrack[thread_num][ numberCandidateThread [thread_num]];
                 
                    
                    
                    L1Branch &tr = CandidatesTrack[thread_num][ numberCandidateThread [thread_num]];
                    
                   tr.BestCandidate = 1;
                    
                    
                    for (vector<THitI>::iterator phitIt = tr.StsHits.begin();/// used strips are marked
                         phitIt != tr.StsHits.end(); ++phitIt){
                             L1StsHit &h = *(const_cast<L1StsHit*>(&((*vStsHits)[*phitIt])));
                             
                             
                             
                         for ( Tindex in = 0; in < (stripF)[h.f].Cand; ++in ){
                         
                        L1Branch* &tr_neig = (stripF)[h.f].Candidates[in];
                         
                       if (!L1Branch::compareCand(tr_neig, &tr)) {tr_neig->BestCandidate = 0; } else {tr.BestCandidate = 0; }
                         
                       }
                      
                       for ( Tindex in = 0; in < (stripB)[h.b].Cand; ++in ){
                         
                          L1Branch* &tr_neig = (stripB)[h.b].Candidates[in];
                         
                      if (!L1Branch::compareCand(tr_neig, &tr)) tr_neig->BestCandidate = 0; else tr.BestCandidate = 0;
                       }
                             
                             
                    
                      omp_set_lock(&stripB[h.b].Occupied);
            
                        (stripB)[h.b].Candidates[(stripB)[h.b].Cand]=&tr;
                        (stripB)[h.b].Cand++; 
                        
                        
                         if ((stripB)[h.b].Cand==1)  
                        
                      {  stripBUsed[thread_num][UsedB[thread_num]]=h.b;
                        
                         UsedB[thread_num]++;}
                                               
                      omp_unset_lock(&stripB[h.b].Occupied);  
                        
                      omp_set_lock(&stripF[h.f].Occupied);
                      
                         (stripF)[h.f].Candidates[(stripF)[h.f].Cand]=&tr;
                         (stripF)[h.f].Cand++;
                         
                     if ((stripF)[h.f].Cand==1)   { stripFUsed[thread_num][UsedF[thread_num]]=h.f;
                         
                         UsedF[thread_num]++;}
                                       
                      omp_unset_lock(&stripF[h.f].Occupied);   
                   
                      
                      
     
                      
                                                    
                    }
                 
          numberCandidateThread [thread_num]++;

        }  // itrip
        
        
        
       }
        
      }
        

    //    cout<<check1<<" "<<check2<<endl; 
         
        size=0;  
    
          

     for (int i=0; i<fNThreads; ++i){
       
       SavedCand[i]=0;    
     SavedHits[i]=0; 
        
       numberCandidateThreadBuf [i]=size;
      
       size+=numberCandidateThread [i];
      }
      

       
//         int can=0;
//         
//         for (int i=0; i<fNThreads; ++i){
//             
//             
//             
//             for ( Tindex iCandidate = 0; iCandidate < numberCandidateThread [i]; ++iCandidate ){
//                 
//                 L1Branch* &tr = CandidatesTrack_buf2[i][ iCandidate];
// 
//         
//         (*candidatesAllP)[can]=(tr);
//         can++;
//             }
//         }

        



    //  Time.Start();

      // cout <<"Level "<< ilev <<" track candidates "<< ntracks << endl;
      // cout << " total number of track candidates " << vtrackcandidate.size() <<endl;
         
           
           
      if (--nlevel == 0) break;
        // -- make competition between tracks of the same length --
        // sort track candidates
      //vector<L1Branch*> vptrackcandidate[glob_threads];  // vptrackcandidate - array of pointers to vtrackcandidate
  //  TStopwatch cand;
 // cand.Start();


//      #pragma omp parallel for
//        for (int i=0; i<fNThreads; ++i){  
//  vTracks_local[i].resize(0);
// vRecoHits_local[i].resize(0);
// 
//        }

     

/*     
       #pragma omp parallel for 
        for (int i=0; i<fNThreads; ++i){
            
            
            for ( Tindex iC = 0; iC < numberCandidateThread [i]; ++iC ){

                     CommonIndex[i][iC]=0;

                }
            }*/
            
            
            
            
            
           
/*



    #pragma omp parallel for 
        for (int i=0; i<fNThreads; ++i){
            
            
            for ( Tindex iC = 0; iC < numberCandidateThread [i]; ++iC ){
                
                L1Branch* &tr = CandidatesTrack_buf2[i][iC];
                
                 
                
                for (vector<THitI>::iterator phIt = tr->StsHits.begin();
                     phIt != tr->StsHits.end(); ++phIt){
                    
                    const   L1StsHit &h = (*vStsHits)[*phIt];
                    
                    for (int k=0; k<(stripF)[h.f].Cand; ++k){
                    
                     Common[i][iC][CommonIndex[i][iC]]=((stripF)[h.f].Candidates[k]);
                     CommonIndex[i][iC]++;
                    }
                    
                      for (int m=0; m<(stripB)[h.b].Cand; ++m){
                        
                        Common[i][iC][CommonIndex[i][iC]]=((stripB)[h.b].Candidates[m]);
                        CommonIndex[i][iC]++;
                    }
                    
                }
            }
            
            }  */ 
 
       
    
//       #pragma omp parallel for 
//           for (int i=0; i<fNThreads; ++i){
//   
//        for ( Tindex iCandidate = 0; iCandidate < numberCandidateThread [i]; ++iCandidate ){ 
//       
//          L1Branch &tr = *(CandidatesTrack_buf2[i][iCandidate]);
// 
//          for (vector<THitI>::iterator phitIt = tr.StsHits.begin();
//                                       phitIt != tr.StsHits.end(); ++phitIt){
//            
//           L1StsHit &h = *(const_cast<L1StsHit*>(&((*vStsHits)[*phitIt])));
         
           
//                     HitCandidatesSortF[h.f].push_back(&tr);
//                     HitCandidatesSortB[h.b].push_back(&tr);
                    
//                        int &indexF =  ((stripF)[h.f].Cand);
//                        int &indexB =  ((stripB)[h.b].Cand);
                      
                     
//                       ((stripB)[h.b].Candidates[indexB])=(&tr);
//                       ((stripF)[h.f].Candidates[indexF])=(&tr);
/*                     
                      omp_set_lock(&hitToBestTrackB[h.b].Occupied);
                        (stripB)[h.b].Candidates.push_back(&tr);
                        omp_unset_lock(&hitToBestTrackB[h.b].Occupied);  
                        
                        omp_set_lock(&hitToBestTrackF[h.f].Occupied);
                         (stripF)[h.f].Candidates.push_back(&tr);
                        omp_unset_lock(&hitToBestTrackF[h.f].Occupied);  
                     */
//                       indexF++;
//                       indexB++;
//                   h.CandidatesB.push_back(&tr);
//                   h.CandidatesF.push_back(&tr);
//             
           
//          }
//        }
//                         }
                        
 
          
//            #pragma omp parallel for 
//      for( Tindex j=0; j<(*vSFlagB).size(); ++j ) { std::sort((stripB)[j].Candidates.begin(),(stripB)[j].Candidates.end(),L1Branch::compareCand); }
//      #pragma omp parallel for 
//      for( Tindex j=0; j<(*vSFlag).size(); ++j ) std::sort((stripF)[j].Candidates.begin(),(stripF)[j].Candidates.end(),L1Branch::compareCand);
     
        #pragma omp parallel for 
       for (int i=0; i<fNThreads; ++i){
         
         
            for ( Tindex iC = 0; iC < numberCandidateThread [i]; ++iC ){
            
             L1Branch* &tr = CandidatesTrack_buf2[i][iC];

           (*candidatesAllP)[numberCandidateThreadBuf[i]+iC]=(tr);
          }
         
         SavedCand[i]=0;    
         SavedHits[i]=0; 
         
       CandidatesTrack_buf2P[i]  = &CandidatesTrack_buf2[i];
       CandidatesTrack_bufP[i]  = &CandidatesTrack_buf[i];
    
   
//      for( Tindex j=0; j<UsedB[i]; ++j ) {
//      int &index = stripBUsed[i][j];
//        std::sort(&((stripB)[index].Candidates[0]),&((stripB)[index].Candidates[0])+(stripB)[index].Cand,L1Branch::compareCand); }
//      
//      for( Tindex j=0; j<UsedF[i]; ++j ) {
//        
//        int &index = stripFUsed[i][j];
//       
//        std::sort(&((stripF)[index].Candidates[0]),&((stripF)[index].Candidates[0])+(stripF)[index].Cand,L1Branch::compareCand);}
       
       }


      Time.Stop();
      ti[isec]["table"] += Time;
      
   
   

        while (size>0)
        {
       
//cout<<raz<<" raz "<<size<<endl;
             
             
//                  #pragma omp parallel for 
//             for (int i=0; i<fNThreads; ++i){
//             cout<< numberCandidateThread [i]<<endl;
//             
//             for ( Tindex iC = 0; iC < numberCandidateThread [i]; ++iC ){
// 
//                     
//                     L1Branch* &tr =  Common[i][iC][0];
//                     
//                  //   if ((tr->Lengtha)==0) continue;
//                     
//                     int &index1 = (*CandidatesTrack_buf2P)[i][iC]->CandIndex;
//                     int &index2 = tr->CandIndex;
//                     
//                     
// 
//                   
//                      if (index1!=index2) continue;
//                      
//                 
//                       (tr->Lengtha)=-1;
//               
//                     
//                       for ( Tindex j = 1; j <CommonIndex[i][iC]; ++j ){
//                         
//                         L1Branch* &trD =  Common[i][iC][j];
//                         
//                 
//                         
//                       (trD->Lengtha)=-1;
//                   
//                       }
// 
//           
// 
//                               for (vector<THitI>::iterator phIt = tr->StsHits.begin();
//                              phIt != tr->StsHits.end(); ++phIt)
//                         {
// 
//                               const   L1StsHit &h = (*vStsHits)[*phIt];
//                               
//                             //  cout<<h.f<<" h.f"<<" "<< (tr->Lengtha)<<endl;
//                             
//                             SetFUsed( const_cast<unsigned char&>((*vSFlag)[h.f]));
//                             SetFUsed( const_cast<unsigned char&>((*vSFlagB)[h.b]));
// 
//                             
//                             vRecoHits_local [omp_get_thread_num()][SavedHits[omp_get_thread_num()]]=(*phIt);
//                             SavedHits[omp_get_thread_num()]++;
//                           
//                           
//                           
//                         }
//                       
//                         L1Track t;
//                         t.NHits = tr->StsHits.size();
//                         t.Momentum = tr->Momentum;
//                         //vTracks_local [omp_get_thread_num()].push_back(t);
//                         vTracks_local [omp_get_thread_num()][SavedCand[omp_get_thread_num()]]=(t);
//                         SavedCand[omp_get_thread_num()]++;
//                         
//                         
//                        // (tr->NHits)=-1;
//                      
//                     
//                 }
//             }
            // cout<<size<<endl;
            
          
             #pragma omp parallel 
            {
              
              
              
              
              
              
              
              
              
                       #pragma omp for nowait  
       for (int i=0; i<fNThreads; ++i){
    
   
     for( Tindex j=0; j<UsedB[i]; ++j ) {
     int &index = stripBUsed[i][j];
       std::sort(&((stripB)[index].Candidates[0]),&((stripB)[index].Candidates[0])+(stripB)[index].Cand,L1Branch::compareCand); }
       }
       
        #pragma omp for nowait  
       for (int i=0; i<fNThreads; ++i){
     for( Tindex j=0; j<UsedF[i]; ++j ) {
       
       int &index = stripFUsed[i][j];
      
       std::sort(&((stripF)[index].Candidates[0]),&((stripF)[index].Candidates[0])+(stripF)[index].Cand,L1Branch::compareCand);}
       
       }
              
            }       
              
              
              
              
              
               Time.Start(); 
               #pragma omp parallel 
            {
              
              
              
                   #pragma omp for nowait 
                for ( Tindex iCandidate = 0; iCandidate < size; ++iCandidate ){
                  
     
                    
                    L1Branch* &tr = ((*candidatesAllP)[iCandidate]);
                    
                  
                     bool allHitsApproved=1;
                    
                    for (vector<THitI>::iterator phIt = tr->StsHits.begin();
                         phIt != tr->StsHits.end(); ++phIt){
                        
                        const   L1StsHit &h = (*vStsHits)[*phIt];
                        
                       
                        
                   //     if ((stripB)[h.b].Cand==0|| (stripF)[h.f].Cand==0) {allHitsApproved=0; break;};
                        
                         
                        
                      //  for ( Tindex in = 0; in < (stripF)[h.f].Cand; ++in )
                        {
                        //    if ((stripF)[h.f].Candidates[in]->Lengtha==0) continue;
                            if ((stripF)[h.f].Candidates[0]!= tr) allHitsApproved=0;
                      //      break;
                        }
                     //   if (allHitsApproved == 0) break;
                        
                     //  for ( Tindex in = 0; in < (stripB)[h.b].Cand; ++in )
                        {
                         //   if ((stripB)[h.b].Candidates[in]->Lengtha==0) continue;
                            if ((stripB)[h.b].Candidates[0]!= tr) allHitsApproved=0;
                      //      break;
                        }
                     //   if (allHitsApproved == 0) break;
                        
                        
                    }
                    
                    if (allHitsApproved == 1)
                        
                    {
                      
                 //    if (!tr->BestCandidate) count++;
                   //   if (!tr->BestCandidate) cout<<tr->BestCandidate<<endl;
                        
                        for (vector<THitI>::iterator phIt = tr->StsHits.begin();
                             phIt != tr->StsHits.end(); ++phIt)
                        {
                            const   L1StsHit &h = (*vStsHits)[*phIt];
                            
                            
                            
                            for ( Tindex in = 0; in <(stripB)[h.b].Cand; ++in ) {
                                
                              //  if ((stripB)[h.b].Candidates[in]==0)  continue;
                                L1Branch &tr_del = *((stripB)[h.b].Candidates[in]);
                              //   if (tr_del.BestCandidate) countD++;
                                (tr_del.Lengtha)=0;
                               
                             
//                                 for (vector<THitI>::iterator phIt1 = tr_del.StsHits.begin();
//                                      phIt1 != tr_del.StsHits.end(); ++phIt1){
//                                     
//                                     const   L1StsHit &h_del = (*vStsHits)[*phIt1];
//                                     
//                                     for ( Tindex in = 0; in < (stripF)[h_del.f].Cand; ++in ) if ((stripF)[h_del.f].Candidates[in]==&tr_del) {(stripF)[h_del.f].Candidates[in]=0; break;}
//                                     for ( Tindex in = 0; in < (stripB)[h_del.b].Cand; ++in ) if ((stripB)[h_del.b].Candidates[in]==&tr_del) {(stripB)[h_del.b].Candidates[in]=0; break;}
//                                     
//                                 }
                            }
                            
                            
                            for ( Tindex in = 0; in < (stripF)[h.f].Cand; ++in ) {
                                
                             //   if ((stripF)[h.f].Candidates[in]==0) continue;
                                
                                L1Branch &tr_del = *((stripF)[h.f].Candidates[in]);
                                //   if (tr_del.BestCandidate) countD++;
                                (tr_del.Lengtha)=0;
                              
                        
//                                 for (vector<THitI>::iterator phIt2 = tr_del.StsHits.begin();
//                                      phIt2 != tr_del.StsHits.end(); ++phIt2){
//                                     
//                                     const   L1StsHit &h_del = (*vStsHits)[*phIt2];
//                                     
//                                     for ( Tindex in = 0; in < (stripF)[h_del.f].Cand; ++in )
//                                         if ((stripF)[h_del.f].Candidates[in]==&tr_del) {(stripF)[h_del.f].Candidates[in]=0; break;}
//                                     for ( Tindex in = 0; in < (stripB)[h_del.b].Cand; ++in )
//                                         if ((stripB)[h_del.b].Candidates[in]==&tr_del) {(stripB)[h_del.b].Candidates[in]=0; break;}
//                                     
//                                 }
                            }
                            
                            
                            
                            
                            
                            
                            //                   (stripF)[h.f].Candidates.resize(0);
                            //                   (stripB)[h.b].Candidates.resize(0);
                            
                            (stripF)[h.f].Cand=0;
                            (stripB)[h.b].Cand=0;
                            
                            
                            SetFUsed( const_cast<unsigned char&>((*vSFlag)[h.f]));
                            SetFUsed( const_cast<unsigned char&>((*vSFlagB)[h.b]));
                            
                            
                            //vRecoHits_local [omp_get_thread_num()].push_back(*phIt);
                            
                            vRecoHits_local [omp_get_thread_num()][SavedHits[omp_get_thread_num()]]=(*phIt);
                            SavedHits[omp_get_thread_num()]++;
                            
                        }
                        
                        
                        
                        
                        L1Track t;
                        t.NHits = tr->StsHits.size();
                        t.Momentum = tr->Momentum;
                        //vTracks_local [omp_get_thread_num()].push_back(t);
                        vTracks_local [omp_get_thread_num()][SavedCand[omp_get_thread_num()]]=(t);
                        SavedCand[omp_get_thread_num()]++;
                        
                        
                        (tr->Lengtha)=0;
                        
                        
                        
                        
                        
                        
                    } //else { if (!tr->BestCandidate == 0) countD++;}
  
            } //NTreads loop end
            } 
            
   Time.Stop();
      ti[isec]["save"] += Time;
            
            
            
            

            

            
            

          

            



          size=0;

           #pragma omp parallel for
          for (int i=0; i<fNThreads; ++i){
            numberCandidateThreadBuf[i]=0;
       for ( Tindex iCandidate = 0; iCandidate < numberCandidateThread [i]; ++iCandidate ){
           L1Branch* &tr = ((*CandidatesTrack_buf2P)[i][iCandidate]);
           if ((tr->Lengtha)!=0)
           {
             (*CandidatesTrack_bufP)[i][numberCandidateThreadBuf[i]]=tr;
             numberCandidateThreadBuf[i]++;
           }
          
             
         }
          #pragma omp atomic
         size+=numberCandidateThreadBuf[i];
          }

    // cout<<count<<" "<<countD<<" d"<<endl;
      
    #pragma omp parallel for
   for (int i=0; i<fNThreads; ++i){
     
            
      numberCandidateThread[i]=numberCandidateThreadBuf[i];


   vector<L1Branch*>* CandidatesTrackTmp = CandidatesTrack_buf2P[i];
   CandidatesTrack_buf2P[i] = CandidatesTrack_bufP[i];
   CandidatesTrack_bufP[i] = CandidatesTrackTmp;
     
     
   }
   
   
              int index=0;
           for (int i=0; i<fNThreads; ++i){
            numberCandidateThread2[i]=index;
            index+=numberCandidateThread[i];
           }
   
   
   
   
                       #pragma omp parallel for 
       for (int i=0; i<fNThreads; ++i){

            for ( Tindex iC = 0; iC < numberCandidateThread [i]; ++iC ){
               L1Branch* &tr = ((*CandidatesTrack_buf2P)[i][iC]);
             
           (*candidatesAllP)[numberCandidateThread2[i]+iC]=(tr);
                   
          }
            
       } 
   



   
           
        } //end of while
       
 
                   
//      ti[isec]["save"] += Time;  
      

//           int size = 1;
//           bool allHitsApproved;
//               while (size>0)
//      {
//
//        
//       
//        
// #pragma omp parallel for  firstprivate(allHitsApproved)
//                  for (int i=0; i<fNThreads; ++i){
//                    
//
//     
//        for ( Tindex iCandidate = 0; iCandidate < numberCandidateThread [i]; ++iCandidate ){
//       
//          L1Branch &tr = *((*CandidatesTrack_buf2P)[i][iCandidate]);
//          allHitsApproved=1;
//            
//
//            for (vector<THitI>::iterator phIt = tr.StsHits.begin();
//                                        phIt != tr.StsHits.end(); ++phIt){
//
//           const   L1StsHit &h = (*vStsHits)[*phIt];
//           
//           
//            
//               if ((stripB)[h.b].Cand==0|| (stripF)[h.f].Cand==0) {allHitsApproved=0; break;};
//                
//        
//              for ( Tindex in = 0; in < (stripF)[h.f].Cand; ++in )
//              {
//                 if ((stripF)[h.f].Candidates[in]==0) continue;
//                 if ((stripF)[h.f].Candidates[in]!= &tr) allHitsApproved=0;
//                 break;
//              }
//                if (allHitsApproved == 0) break;
//
//                      for ( Tindex in = 0; in < (stripB)[h.b].Cand; ++in )
//                        {
//                         if ((stripB)[h.b].Candidates[in]==0) continue;
//                         if ((stripB)[h.b].Candidates[in]!= &tr) allHitsApproved=0;
//                         break;
//                        }  
//                     if (allHitsApproved == 0) break;
//         
//                        
//            } 
//            
//            if (allHitsApproved == 1)
//           
//            {
//              
//                     for (vector<THitI>::iterator phIt = tr.StsHits.begin();
//                                            phIt != tr.StsHits.end(); ++phIt)
//                 {
//                  const   L1StsHit &h = (*vStsHits)[*phIt];
//             
//                  
//                  
//                  for ( Tindex in = 0; in <(stripB)[h.b].Cand; ++in ) {
//               
//                    if ((stripB)[h.b].Candidates[in]==0)  continue;
//                   L1Branch &tr_del = *((stripB)[h.b].Candidates[in]);
//                   (tr_del.NHits)=-1;
//                   for (vector<THitI>::iterator phIt1 = tr_del.StsHits.begin();
//                                           phIt1 != tr_del.StsHits.end(); ++phIt1){
//
//                   const   L1StsHit &h_del = (*vStsHits)[*phIt1];
//  
//                                    for ( Tindex in = 0; in < (stripF)[h_del.f].Cand; ++in ) if ((stripF)[h_del.f].Candidates[in]==&tr_del) {(stripF)[h_del.f].Candidates[in]=0; break;}
//                                    for ( Tindex in = 0; in < (stripB)[h_del.b].Cand; ++in ) if ((stripB)[h_del.b].Candidates[in]==&tr_del) {(stripB)[h_del.b].Candidates[in]=0; break;}
//                                    
//                   }                    
//                  }   
//
//                  
//                   for ( Tindex in = 0; in < (stripF)[h.f].Cand; ++in ) {
//                   
//                     if ((stripF)[h.f].Candidates[in]==0) continue;
//                     
//                   L1Branch &tr_del = *((stripF)[h.f].Candidates[in]);
//                   (tr_del.NHits)=-1;
//                   for (vector<THitI>::iterator phIt2 = tr_del.StsHits.begin();
//                                           phIt2 != tr_del.StsHits.end(); ++phIt2){
//
//                   const   L1StsHit &h_del = (*vStsHits)[*phIt2];
//  
//                                    for ( Tindex in = 0; in < (stripF)[h_del.f].Cand; ++in ) 
//                                      if ((stripF)[h_del.f].Candidates[in]==&tr_del) {(stripF)[h_del.f].Candidates[in]=0; break;}
//                                    for ( Tindex in = 0; in < (stripB)[h_del.b].Cand; ++in ) 
//                                      if ((stripB)[h_del.b].Candidates[in]==&tr_del) {(stripB)[h_del.b].Candidates[in]=0; break;}
//                                    
//                   }             
//                  }      
//                  
//                  
//  
//
//                                    
//
////                   (stripF)[h.f].Candidates.resize(0);
////                   (stripB)[h.b].Candidates.resize(0);
//                  
//                  (stripF)[h.f].Cand=0;
//                  (stripB)[h.b].Cand=0;
//
//
//                  SetFUsed( const_cast<unsigned char&>((*vSFlag)[h.f]));
//                  SetFUsed( const_cast<unsigned char&>((*vSFlagB)[h.b]));
//
//           
//                  //vRecoHits_local [omp_get_thread_num()].push_back(*phIt);
//                  
//                  vRecoHits_local [i][SavedHits[i]]=(*phIt);
//                  SavedHits[i]++;
//                  
//                }
//              
//              
//              
//              
//              L1Track t;
//              t.NHits = tr.StsHits.size();
//              t.Momentum = tr.Momentum;
//              //vTracks_local [omp_get_thread_num()].push_back(t);
//              vTracks_local [i][SavedCand[i]]=(t);
//              SavedCand[i]++;
//
//            
//              (tr.NHits)=-1;
//              
//              
//
//              
//         
//                
//                          }
//          }
//          
//
//                  } //NTreads loop end
//   
//        
//                   
//                  
//            #pragma omp parallel for
//           for (int i=0; i<fNThreads; ++i){
//             numberCandidateThreadBuf[i]=0;
//        for ( Tindex iCandidate = 0; iCandidate < numberCandidateThread [i]; ++iCandidate ){
//            L1Branch &tr = *((*CandidatesTrack_buf2P)[i][iCandidate]);
//            if ((tr.NHits)!=-1)
//            {
//              (*CandidatesTrack_bufP)[i][numberCandidateThreadBuf[i]]=&tr;
//              numberCandidateThreadBuf[i]++;
//            }
//      
//              
//          }
//           }
// 
//             
//      
//          size=0;
//    
//    for (int i=0; i<fNThreads; ++i){
//      
//             
//       numberCandidateThread[i]=numberCandidateThreadBuf[i];
//
//
//    vector<L1Branch*>* CandidatesTrackTmp = CandidatesTrack_buf2P[i];
//    CandidatesTrack_buf2P[i] = CandidatesTrack_bufP[i];
//    CandidatesTrack_bufP[i] = CandidatesTrackTmp;
//      
//      
//     size+=numberCandidateThread [i];}
//        
//      } //end of while
//    
 
        
         
          vector<int> offset_tracks(fNThreads,NTracksIsecAll);   
          vector<int> offset_hits(fNThreads,NHitsIsecAll);  
          
          NTracksIsecAll+=SavedCand[0];
          NHitsIsecAll+=SavedHits[0];
          

     for (int i=1; i<fNThreads; ++i){
        
       offset_tracks [i]=offset_tracks [i-1]+SavedCand[i-1];
       offset_hits [i]=offset_hits [i-1]+SavedHits[i-1];
       NTracksIsecAll+=SavedCand[i];
       NHitsIsecAll+=SavedHits[i];
      }





#pragma omp parallel for
      for (int i=0; i<fNThreads; ++i){
        
        
            for ( Tindex iC = 0; iC < SavedCand[i]; ++iC ){
          vTracks[ offset_tracks[i]+iC] = (vTracks_local[i][iC]);
        }
        
           for ( Tindex iH = 0; iH < SavedHits[i]; ++iH ){ 
          vRecoHits[offset_hits [i]+iH] = (vRecoHits_local[i][iH]);
                    }
   
     for( Tindex j=0; j<UsedB[i]; ++j ) {
     int &index = stripBUsed[i][j];
       (stripB)[index].Cand=0; }
     
     for( Tindex j=0; j<UsedF[i]; ++j ) {
       
       int &index = stripFUsed[i][j];
      
       (stripF)[index].Cand=0;}

      }


      Time.Stop();
      ti[isec]["save"] += Time;
   
    } //ilev
    


 #if 0
 
   for (int i=0; i<fNThreads; ++i){
       CandidatesTrack_buf2P[i]  = &CandidatesTrack_buf2[i];
       CandidatesTrack_bufP[i]  = &CandidatesTrack_buf[i];
   }
   
      int SavedCand[20]={0};    
      int SavedHits[20]={0};   
      int size =1;
      
      while (size>0) 
      { Time.Start();
      
      cout<<size<<" size"<<endl;
      
         #pragma omp parallel for 
                  for (int i=0; i<fNThreads; ++i){
   
        for ( Tindex iCandidate = 0; iCandidate < numberCandidateThread [i]; ++iCandidate ){
            
           cout<<iCandidate<<endl;
       
          L1Branch &tr = *((*CandidatesTrack_buf2P)[i][iCandidate]);
            
           

          for (vector<THitI>::iterator phitIt = tr.StsHits.begin();/// used strips are marked
                                       phitIt != tr.StsHits.end(); ++phitIt){
            const L1StsHit &h = (*vStsHits)[*phitIt];
       
              omp_set_lock(& hitToBestTrackB[h.b].Occupied); 
            
            {
              if ((hitToBestTrackB[i][h.b].trackCandidateIndex==-1  )||( (tr.Lengtha>hitToBestTrackB[i][h.b].Length) ||
                   ( (tr.Lengtha==hitToBestTrackB[i][h.b].Length)&&
                     ( (tr.ista<hitToBestTrackB[i][h.b].ista) ||
                       ( (tr.ista==hitToBestTrackB[i][h.b].ista) &&
                         (tr.chi2<hitToBestTrackB[i][h.b].Chi2Track  )|| ((tr.chi2==hitToBestTrackB[i][h.b].Chi2Track  )&&
                         (tr.CandIndex>hitToBestTrackB[i][h.b].trackCandidateIndex  )))
                       )
                     )
                   )
                 )// if chi2 considered candidate better than saved one =????
              {
        
                hitToBestTrackB[i][h.b].trackCandidateIndex = (tr.CandIndex);
               
                hitToBestTrackB[i][h.b].Chi2Track = (tr.chi2);
              
                hitToBestTrackB[i][h.b].Length = (tr.Lengtha);
             
                hitToBestTrackB[i][h.b].ista = (tr.ista);
               
              };
            }
            
          omp_unset_lock(&hitToBestTrackB[h.b].Occupied);
              
                  omp_set_lock(&hitToBestTrackF[h.f].Occupied);
            {    
              if ( (hitToBestTrackF[i][h.f].trackCandidateIndex==-1  )||((tr.Lengtha>hitToBestTrackF[i][h.f].Length) ||
                   ( (tr.Lengtha==hitToBestTrackF[i][h.f].Length)&&
                     ( (tr.ista<hitToBestTrackF[i][h.f].ista) ||
                       ( (tr.ista==hitToBestTrackF[i][h.f].ista) &&
                         (tr.chi2<hitToBestTrackF[i][h.f].Chi2Track )|| ((tr.chi2==hitToBestTrackF[i][h.f].Chi2Track  )&&
                         (tr.CandIndex>hitToBestTrackF[i][h.f].trackCandidateIndex  )))
                       )
                     )
                   )
                 )
              {
                { 
                  hitToBestTrackF[i][h.f].trackCandidateIndex = (tr.CandIndex);
                  hitToBestTrackF[i][h.f].Chi2Track = (tr.chi2);
                  hitToBestTrackF[i][h.f].Length = (tr.Lengtha);
                  hitToBestTrackF[i][h.f].ista = (tr.ista);
                 
               
                }
              }
              
            }
           omp_unset_lock(&hitToBestTrackF[h.f].Occupied);  
          }
        }
                  }
                  
                 c_timer.Stop(); 
                  #pragma omp parallel for schedule(dynamic, 40) 
                   for (int k=0; k< (*vSFlag).size(); ++k){
                       for (int i=1; i<fNThreads; ++i){
                     
                     
                                   if ( (hitToBestTrackF[0][k].trackCandidateIndex==-1  )||((hitToBestTrackF[i][k].Length > hitToBestTrackF[0][k].Length) ||
                   ( (hitToBestTrackF[i][k].Length==hitToBestTrackF[0][k].Length)&&
                     ( (hitToBestTrackF[i][k].ista<hitToBestTrackF[0][k].ista) ||
                       ( (hitToBestTrackF[i][k].ista==hitToBestTrackF[0][k].ista) &&
                         (hitToBestTrackF[i][k].Chi2Track<hitToBestTrackF[0][k].Chi2Track )|| ((hitToBestTrackF[i][k].Chi2Track==hitToBestTrackF[0][k].Chi2Track  )&&
                         (hitToBestTrackF[i][k].trackCandidateIndex>hitToBestTrackF[0][k].trackCandidateIndex  )))
                       )
                     )
                   )
                 )
              {
                { 
                     if  (hitToBestTrackF[i][k].trackCandidateIndex==-1  ) continue;
                  hitToBestTrackF[0][k].trackCandidateIndex = hitToBestTrackF[i][k].trackCandidateIndex;
                  hitToBestTrackF[0][k].Chi2Track = hitToBestTrackF[i][k].Chi2Track;
                  hitToBestTrackF[0][k].Length = hitToBestTrackF[i][k].Length;
                  hitToBestTrackF[0][k].ista = hitToBestTrackF[i][k].ista;
                 
               
                }
              }
                     
                     
                       }   
                   }
                     
                     
                     
                     #pragma omp parallel for schedule(dynamic, 40)
                      for (int k=0; k< (*vSFlagB).size(); ++k){
                       for (int i=1; i<fNThreads; ++i){
                     
                     
                                   if ( (hitToBestTrackB[0][k].trackCandidateIndex==-1  )||((hitToBestTrackB[i][k].Length > hitToBestTrackB[0][k].Length) ||
                   ( (hitToBestTrackB[i][k].Length==hitToBestTrackB[0][k].Length)&&
                     ( (hitToBestTrackB[i][k].ista<hitToBestTrackB[0][k].ista) ||
                       ( (hitToBestTrackB[i][k].ista==hitToBestTrackB[0][k].ista) &&
                         (hitToBestTrackB[i][k].Chi2Track<hitToBestTrackB[0][k].Chi2Track )|| ((hitToBestTrackB[i][k].Chi2Track==hitToBestTrackB[0][k].Chi2Track  )&&
                         (hitToBestTrackB[i][k].trackCandidateIndex>hitToBestTrackB[0][k].trackCandidateIndex  )))
                       )
                     )
                   )
                 )
              {
                { 
                   if (hitToBestTrackB[i][k].trackCandidateIndex==-1  ) continue;
                  hitToBestTrackB[0][k].trackCandidateIndex = hitToBestTrackB[i][k].trackCandidateIndex;
                  hitToBestTrackB[0][k].Chi2Track = hitToBestTrackB[i][k].Chi2Track;
                  hitToBestTrackB[0][k].Length = hitToBestTrackB[i][k].Length;
                  hitToBestTrackB[0][k].ista = hitToBestTrackB[i][k].ista;
                 
               
                }
              }
                     
                     
                       }    
                        
                      }
                  
                
                  
                c_timer.Start(0);  
                  
                  
                  
                  
/*      Time.Stop();
      ti[isec]["save"] += Time;
      Time.Start()*/;
                  
             bool allHitsApproved =1;      
 #pragma omp parallel for firstprivate(allHitsApproved)
                  for (int i=0; i<fNThreads; ++i){
     int Saved = SavedCand[i];
     int SavedHit=SavedHits[i];
        for ( Tindex iCandidate = 0; iCandidate < numberCandidateThread [i]; ++iCandidate ){ 
       
          L1Branch &tr = *((*CandidatesTrack_buf2P)[i][iCandidate]);
           

            for (vector<THitI>::iterator phIt = tr.StsHits.begin();
                                        phIt != tr.StsHits.end(); ++phIt){

              const L1StsHit &h = (*vStsHits)[*phIt];

             
   
        
              if ((((hitToBestTrackF[0][h.f].trackCandidateIndex)!= (tr.CandIndex ))) ||
                  (((hitToBestTrackB[0][h.b].trackCandidateIndex)!= (tr.CandIndex ) )))
              {
                allHitsApproved =0;
                 break;
              }
     
             
            } 
            if (allHitsApproved)
           
            {
                
              
                L1Track t;//=vTracks_local [omp_get_thread_num()][tracks_loc[omp_get_thread_num()]];
               tracks_loc[omp_get_thread_num()]++;
              t.NHits = tr.StsHits.size();
              t.Momentum = tr.Momentum;
              for( int i=0; i<6; i++) t.TFirst[i] = t.TLast[i]=0; // CHECKME don't need this
              for( int i=0; i<15; i++ ) t.CFirst[i] = t.CLast[i] = 10;
              t.chi2 = 100;

             vTracks.push_back(t);
             vTracks_local [omp_get_thread_num()].push_back(t);
              vTracks_local [omp_get_thread_num()][Saved]=(t);
              Saved++;

              (tr.FlagRemove)=1;

              {
                for (vector<THitI>::iterator phIt = tr.StsHits.begin();
                                            phIt != tr.StsHits.end(); ++phIt){
                  const L1StsHit &h = (*vStsHits)[*phIt];
                 
                 
                  hitToBestTrackF[0][h.f].UsedByTrack = 1;
                  hitToBestTrackB[0][h.b].UsedByTrack = 1;
                  SetFUsed( const_cast<unsigned char&>((*vSFlag)[h.f]));
                  SetFUsed( const_cast<unsigned char&>((*vSFlagB)[h.b]));

           
                 vRecoHits_local [omp_get_thread_num()].push_back(*phIt);
                  vRecoHits_local [omp_get_thread_num()][SavedHit]=(*phIt);
                  SavedHit++;
             
            
                }
                
              }
            } else allHitsApproved=1;
          }
          
          SavedCand[i]=Saved;
          SavedHits[i]=SavedHit;
                  }

         bool NotAllHitsApproved =0;

       #pragma omp parallel for firstprivate(NotAllHitsApproved)
                   for (int i=0; i<fNThreads; ++i){
     
        for ( Tindex iCandidate = 0; iCandidate < numberCandidateThread [i]; ++iCandidate ){ 
       
          L1Branch &tr = *((*CandidatesTrack_buf2P)[i][iCandidate]);

          
          
          
            for (vector<THitI>::iterator phIt = tr.StsHits.begin();
                                         phIt != tr.StsHits.end(); ++phIt){
              const L1StsHit &h = (*vStsHits)[*phIt];      
         
              if (( hitToBestTrackF[0][h.f].UsedByTrack)||( hitToBestTrackB[0][h.b].UsedByTrack)) 
              {
                NotAllHitsApproved =1;
                break;
              };
            }
            if (NotAllHitsApproved)  
            {
               NotAllHitsApproved =0;
              (tr.FlagRemove)=1;
              
              for (vector<THitI>::iterator phIt = tr.StsHits.begin();
                                           phIt != tr.StsHits.end(); ++phIt){

                const L1StsHit &h = (*vStsHits)[*phIt];
            
               
               omp_set_lock(&hitToBestTrackF[0][h.f].Occupied); 
                
                
                  for (int i=0; i<fNThreads; ++i){
                    hitToBestTrackF[i][h.f].trackCandidateIndex = -1;
                  }
                  
                omp_unset_lock(&hitToBestTrackF[0][h.f].Occupied);
                omp_set_lock(&hitToBestTrackB[0][h.b].Occupied);    
               
                for (int i=0; i<fNThreads; ++i){
                    hitToBestTrackB[i][h.b].trackCandidateIndex = -1;
                }
                
               omp_unset_lock(&hitToBestTrackB[0][h.b].Occupied); 
                
              }
            } 
          }
                   }
                   
          
          
          

          
                  
            #pragma omp parallel for
           for (int i=0; i<fNThreads; ++i){
             numberCandidateThreadBuf[i]=0;
        for ( Tindex iCandidate = 0; iCandidate < numberCandidateThread [i]; ++iCandidate ){
            L1Branch &tr = *((*CandidatesTrack_buf2P)[i][iCandidate]);
            if ((tr.FlagRemove)==0)
            {
              
              (*CandidatesTrack_bufP)[i][numberCandidateThreadBuf[i]]=&tr;
              numberCandidateThreadBuf[i]++;
            }
       
          }
           }
 
             
      
          size=0;
           #pragma omp parallel for reduction(+: size)
    for (int i=0; i<fNThreads; ++i){
      
       numberCandidateThread[i]=numberCandidateThreadBuf[i];


    vector<L1Branch*>* CandidatesTrackTmp = CandidatesTrack_buf2P[i];
    CandidatesTrack_buf2P[i] = CandidatesTrack_bufP[i];
    CandidatesTrack_bufP[i] = CandidatesTrackTmp;
      
      
     size+=numberCandidateThread [i];}

      } //end of while


       
    
       
         vector<int> offset_tracks(fNThreads,NTracksIsecAll);   
          vector<int> offset_hits(fNThreads,NHitsIsecAll);  
          
          NTracksIsecAll+=SavedCand[0];
          NHitsIsecAll+=SavedHits[0];
          
           Time.Stop();
          Time.Start();
     for (int i=1; i<fNThreads; ++i){
        
       offset_tracks [i]=offset_tracks [i-1]+SavedCand[i-1];
       offset_hits [i]=offset_hits [i-1]+SavedHits[i-1];
       NTracksIsecAll+=SavedCand[i];
       NHitsIsecAll+=SavedHits[i];
      }
      

     #pragma omp parallel for  
      for (int i=0; i<fNThreads; ++i){
    
        for ( Tindex iC = 0; iC < SavedCand[i]; ++iC ){
          vTracks[ offset_tracks[i]+iC] = (vTracks_local[i][iC]);
        }

      }
      
      
       #pragma omp parallel for  
      for (int i=0; i<fNThreads; ++i){
        for ( Tindex iH = 0; iH < SavedHits[i]; ++iH ){ 
          vRecoHits[offset_hits [i]+iH] = (vRecoHits_local[i][iH]);
                    }
          
       

      }
        
   

      

//                   Time.Stop();
//       ti[isec]["delete"] += Time;
//       Time.Start();
       
  
#ifdef XXX
/*
for (int i=0; i<glob_threads; i++)
      {
        
       Tindex Bsize=0, nb=0;
      for (vector<L1Branch*>::iterator   trIt = vptrackcandidate.begin();
                                         trIt != vptrackcandidate.end(); ++trIt){
        L1Branch *tr = *trIt;
        Bsize+= sizeof(L1Branch) + sizeof(THitI) * tr->StsHits.size();
        nb++;
      }

      if( Bsize>stat_max_BranchSize ) stat_max_BranchSize = Bsize;
      if( nb>stat_max_n_branches ) stat_max_n_branches = nb;
      }*/
#endif

      //cout<< " track candidates stored: "<< ntracks <<endl;
#ifdef COUNTERS
   for (int i=0; i<omp_get_num_threads(); ++i) stat_nTrCandidates[isec] += vtrackcandidate[i].size();
#endif


    } //ilev
 #endif  
 
//       c_timer.Stop();
//   
//  c_timer.Start(0);
    
      c_timer.Stop();
      ti[isec]["tracks"] = c_timer;
      c_timer.Start(1);
    
          for(int ista=0; ista < NStations; ++ista){
      #pragma omp parallel for   
      for(THitI ih = StsHitsUnusedStartIndex[ista]; ih < StsHitsUnusedStopIndex[ista]; ++ih){

        //L1StsHit &hit = (*vStsHitsUnused)[ih];
        (*vStsHitsUnused)[ih].used = ( GetFUsed( (*vSFlag)[(*vStsHitsUnused)[ih].f] | (*vSFlagB)[(*vStsHitsUnused)[ih].b] ) );
       
      }
 }





      if (isec!=2)

      {
   for (int istal = NStations-2; istal >= FIRSTCASTATION; istal--) 
     n_g1[portionStopIndex[istal]]=Portion;


   
   
          int NHitsOnStation=0;
    
              for(int ista = 0; ista < NStations; ++ista){
              
              int Nelements = 0;
               int NHitsOnStationTmp=NHitsOnStation;
                  vGrid[ista].UpdateIterGrid(Nelements, &((*vStsHitsUnused)[StsHitsUnusedStartIndex[ista]]),
                                             RealIHitPBuf, &((*RealIHitP)[StsHitsUnusedStartIndex[ista]]),
                                             &RealIHit_v_buf2, vStsHitsUnused_buf, vStsHitPointsUnused_buf,
                                             &((*vStsHitPointsUnused)[StsHitsUnusedStartIndex[ista]]), NHitsOnStation);
              StsHitsUnusedStartIndex[ista] = NHitsOnStationTmp;
              StsHitsUnusedStopIndex[ista] =NHitsOnStation;
              
          }

   



    
    vector<THitI>* RealIHitPTmp = RealIHitP;
    RealIHitP = RealIHitPBuf;
    RealIHitPBuf = RealIHitPTmp;
    
    std::vector< L1StsHit > *vStsHitsUnused_temp = vStsHitsUnused;
    vStsHitsUnused = vStsHitsUnused_buf;
    vStsHitsUnused_buf = vStsHitsUnused_temp;
    
    std::vector< L1HitPoint > *vStsHitsUnused_temp2 = vStsHitPointsUnused;
    vStsHitPointsUnused = vStsHitPointsUnused_buf;
    vStsHitPointsUnused_buf = vStsHitsUnused_temp2;
          

     #ifdef XXX
          c_timer.Stop();
          ti[isec]["finish"] = c_timer;
     #endif


#ifdef XXX
//     if( stat_max_n_trip<stat_n_trip ) stat_max_n_trip = vTriplets.size();
//     Tindex tsize = vTripletsP.size()*sizeof(L1Triplet);
//     if( stat_max_trip_size < tsize ) stat_max_trip_size = tsize;
#endif
// #ifdef DRAW
//     draw.ClearVeiw();
// //   draw.DrawInfo();
//     draw.DrawRestHits(StsHitsUnusedStartIndex, StsHitsUnusedStopIndex, RealIHit);
//     draw.DrawRecoTracks();
//     draw.SaveCanvas("Reco_"+isec+"_");
//     draw.DrawAsk();
// #endif

// #ifdef PULLS
//       fL1Pulls->Build(1);
// #endif
#ifdef COUNTERS
    stat_nHits[isec] += vStsHitsUnused->size();
    
    cout << "iter = " << isec << endl;
    cout << " NSinglets = " << stat_nSinglets[isec]/stat_N << endl;
//    cout << " NDoublets = " << stat_nDoublets[isec]/stat_N << endl;
    cout << " NTriplets = " << stat_nTriplets[isec]/stat_N << endl;
    cout << " NHitsUnused = " << stat_nHits[isec]/stat_N << endl;

#endif // COUNTERS


 

  
    
      }
      

    
   
//  trash.Stop();
//   cout<<trash.RealTime()<<" trash 2500 "<<omp_get_num_threads()<<" omp_set_num_threads"<<endl;
  } // for (int isec



 









#ifdef XXX

  c_timerG.Stop();
  gti["iterts"] = c_timerG;
  c_timerG.Start(1);
#endif

#ifdef MERGE_CLONES
 // CAMergeClones();
#endif
  
#ifdef XXX
  c_timerG.Stop();
  gti["merge "] = c_timerG;
#endif
 
  //cout<< "Total number of tracks found " << vTracks.size() <<endl;

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

  c_time.Stop();

//   cout << "End TrackFinder" << endl;
//  CATime = (double(c_time.CpuTime()));
  CATime = (double(c_time.RealTime()));

#ifdef XXX
  

  cout << endl << " --- Timers, ms --- " << endl;
  ti.Calc();
  stat_ti += ti;
  L1CATFTimerInfo tmp_ti = stat_ti/0.001/stat_N; // ms
  
  tmp_ti.PrintReal();
  stat_gti += gti;
  L1CATFIterTimerInfo tmp_gti = stat_gti/0.001/stat_N; // ms
  tmp_gti.PrintReal( 1 );
  fstream filestr;
  filestr.open("speedUp.log", fstream::out | fstream::app);
  float tripl_speed=1000./(tmp_ti.GetTimerAll()["tripl1"].Real());
  filestr << tripl_speed << " ";
  filestr.close();


#if 0
  static long int NTimes =0, NHits=0, HitSize =0, NStrips=0, StripSize =0, NStripsB=0, StripSizeB =0,
    NDup=0, DupSize=0, NTrip=0, TripSize=0, NBranches=0, BranchSize=0, NTracks=0, TrackSize=0 ;

  NTimes++;
  NHits += vStsHitsUnused->size();
  HitSize += vStsHitsUnused->size()*sizeof(L1StsHit);
  NStrips+= vStsStrips.size();
  StripSize += vStsStrips.size()*sizeof(fscal) + (*vSFlag).size()*sizeof(unsigned char);
  NStripsB+= (*vSFlagB).size();
  StripSizeB += vStsStripsB.size()*sizeof(fscal) + (*vSFlagB).size()*sizeof(unsigned char);
  NDup += stat_max_n_dup;
  DupSize += stat_max_n_dup*sizeof(/*short*/ int);
  NTrip += stat_max_n_trip;
  TripSize += stat_max_trip_size;

  NBranches += stat_max_n_branches;
  BranchSize += stat_max_BranchSize;
  NTracks += vTracks.size();
  TrackSize += sizeof(L1Track)*vTracks.size() + sizeof(THitI)*vRecoHits.size();
  int k = 1024*NTimes;

  cout<<"L1 Event size: \n"
      <<HitSize/k<<"kB for "<<NHits/NTimes<<" hits, \n"
      <<StripSize/k<<"kB for "<<NStrips/NTimes<<" strips, \n"
      <<StripSizeB/k<<"kB for "<<NStripsB/NTimes<<" stripsB, \n"
      <<DupSize/k<<"kB for "<<NDup/NTimes<<" doublets, \n"
      <<TripSize/k<<"kB for "<<NTrip/NTimes<<" triplets\n"
      <<BranchSize/k<<"kB for "<<NBranches/NTimes<<" branches, \n"
      <<TrackSize/k<<"kB for "<<NTracks/NTimes<<" tracks. "<<endl;
  cout<<" L1 total event size = "
      <<( HitSize + StripSize +  StripSizeB + DupSize + TripSize + BranchSize + TrackSize )/k
      <<" Kb"<<endl;
#endif // 0
  //cout << "===> ... CA Track Finder    " << endl << endl;
#endif

#ifdef DRAW
  draw.ClearVeiw();
//   draw.DrawInputHits();
//   draw.DrawInfo();
  draw.DrawRecoTracks();
  draw.SaveCanvas("Reco_");
  draw.DrawAsk();
#endif 
#ifdef PULLS
  static int iEvee = 0;
  iEvee++;
  if (iEvee%1 == 0)
    fL1Pulls->Build(1);
#endif
#ifdef DOUB_PERFORMANCE
  fL1Eff_doublets->CalculateEff();
  fL1Eff_doublets->Print("Doublets performance.",1);
#endif
#ifdef TRIP_PERFORMANCE
  fL1Eff_triplets->CalculateEff();
  fL1Eff_triplets->Print("Triplet performance",1);
//   fL1Eff_triplets2->CalculateEff();
//   fL1Eff_triplets2->Print("Triplet performance. After chi2 cut");
#endif


}


  /** *************************************************************
   *                                                              *
   *     The routine performs recursive search for tracks         *
   *                                                              *
   *     I. Kisel                                    06.03.05     *
   *     I.Kulakov                                    2012        *
   *                                                              *
   ****************************************************************/
  
inline void L1Algo::CAFindTrack(int ista,
                         L1Branch &best_tr, unsigned char &best_L, fscal &best_chi2,
                         const L1Triplet* curr_trip,
                         L1Branch &curr_tr, unsigned char &curr_L, fscal &curr_chi2,
                         unsigned char min_best_l,
                         int &NCalls, L1Branch *new_tr)
             /// recursive search for tracks
             /// input: @ista - station index, @&best_tr - best track for the privious call, @&best_L -
                         /// output: @&NCalls - number of function calls             
{
   if (curr_trip->GetLevel() == 0){ // the end of the track -> check and store
    
      // -- finish with current track
      // add rest of hits
    THitI ihitm = curr_trip->GetMHit();
   
    if(!GetFUsed( (*vSFlag)[(*vStsHitsUnused)[ihitm].f] | (*vSFlagB)[(*vStsHitsUnused)[ihitm].b] )  ){
       // curr_tr.StsHits.resize(curr_tr.CandIndex+1);
        curr_tr.StsHits[curr_tr.NHits]=((*RealIHitP)[ihitm]);
        
       // curr_tr.StsHits.push_back((*RealIHitP)[ihitm]);
       
        curr_tr.NHits++;
               
      curr_L++;
    }
     THitI ihitr = curr_trip->GetRHit();
    if( !GetFUsed( (*vSFlag)[(*vStsHitsUnused)[ihitr].f] | (*vSFlagB)[(*vStsHitsUnused)[ihitr].b] ) ){
       // curr_tr.StsHits.resize(curr_tr.CandIndex+1);
      curr_tr.StsHits[curr_tr.NHits]=((*RealIHitP)[ihitr]);
        
        
      //  curr_tr.StsHits.push_back((*RealIHitP)[ihitr]);
     
        curr_tr.NHits++;
             
      curr_L++;
    }
    
    //if( curr_L < min_best_l - 1 ) return; // suppouse that only one hit can be added by extender
    if( curr_chi2 > TRACK_CHI2_CUT * (curr_L*2-5) ) return;

//       // try to find more hits
// #ifdef EXTEND_TRACKS
//     // if( curr_L < min_best_l )
//     if (isec != kFastPrimJumpIter && isec != kAllPrimJumpIter && isec != kAllSecJumpIter && curr_L >= 3){
//       //curr_chi2 = BranchExtender(curr_tr);
//       BranchExtender(curr_tr);
//       curr_L = curr_tr.StsHits.size();
//         //      if( 2*curr_chi2 > (2*(curr_L*2-5) + 1) * 4*4 ) return;
//     }
// #endif // EXTEND_TRACKS
    
      // -- select the best
    if ( (curr_L > best_L ) || ( (curr_L == best_L) && (curr_chi2 < best_chi2) ) ){
     // best_tr.CandIndex   = curr_tr.CandIndex;
        best_tr   = curr_tr;
      best_chi2 = curr_chi2;
      best_L    = curr_L;
    }

  } else { //MEANS level ! = 0

    
    
    
    // L1Triplet* new_trip = curr_trip->neibourBest;
      // try to extend. try all possible triplets
    Tindex NNeighs = curr_trip->neighbours.size();
     
    for (Tindex in = 0; in < NNeighs; in++) {
      
      
    
      
       // if ( GetFUsed( (*vSFlag)[(*vStsHitsUnused)[new_trip->GetLHit()].f] | (*vSFlagB)[(*vStsHitsUnused)[new_trip->GetLHit()].b] )))*/ 
      
/*      Tindex index = curr_trip->neighbours[in];// + offset;*/
      L1Triplet* new_trip = curr_trip->neighbours[in];

    

      if ( GetFUsed( (*vSFlag)[(*vStsHitsUnused)[new_trip->GetLHit()].f] | (*vSFlagB)[(*vStsHitsUnused)[new_trip->GetLHit()].b] )){ //hits are used
          //  no used hits allowed -> compare and store track
        if ( ( curr_L > best_L ) || ( (curr_L == best_L) && (curr_chi2 < best_chi2) ) ){
            best_tr = curr_tr;
           // best_tr.CandIndex = curr_tr.CandIndex;
            
          best_chi2 = curr_chi2;
          best_L    = curr_L;
            
        }
          
      }
      else{ // if hits are used add new triplet to the current track


          
        new_tr[ista] = curr_tr;
          
        unsigned char new_L = curr_L;
        fscal new_chi2 = curr_chi2;

          // add new hit
        //  new_tr.StsHits.resize(new_tr.CandIndex+1);
        new_tr[ista].StsHits[new_tr[ista].NHits]=((*RealIHitP)[new_trip->GetLHit()]);
         
          new_tr[ista].NHits++;
         // NHitTrack++;
         // if (NHitTrack-1!=new_tr.StsHits.size()) cout<<"SARAERSER"<<endl;
        new_L += 1;
       
          //   dqp = dqp/Cqp; // IKu
          //   new_chi2 += dqp*dqp*5.;
      
         
              // check new triplet
      const fscal qp1 = curr_trip->GetQp();
      const fscal qp2 = new_trip->GetQp();
      fscal dqp = abs(qp1 - qp2);
      fscal Cqp = curr_trip->Cqp;
      Cqp      += new_trip->Cqp;
      
       dqp = dqp/Cqp*5.;  // CHECKME: understand 5, why no sqrt(5)?
        new_chi2 += dqp*dqp;

          
          if( new_chi2 > TRACK_CHI2_CUT * new_L ) continue;
                  //  new_tr[ista].StsHits.push_back((*RealIHitP)[new_trip->GetLHit()]);

          
        const int new_ista = ista + new_trip->GetMSta() - new_trip->GetLSta();
        

        
        CAFindTrack(new_ista, best_tr, best_L, best_chi2, new_trip, new_tr[ista], new_L, new_chi2, min_best_l, NCalls, new_tr);
        
        
        NCalls++;
         
        
      } // add triplet to track

      // break;

    } // for neighbours
  } // level = 0
}