/*
 *=====================================================
 *
 *  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"
#include "TRandom.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, fvec *HitTime, fvec *Event
                        )
{

    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];
        
        
        HitTime[i1_V][i1_4]=hitl.time;
        
#ifdef USE_EVENT_NUMBER
        Event[i1_V][i1_4]=hitl.n;
#endif USE_EVENT_NUMBER

        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::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,
                        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( int 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];
        
        
        const fvec &dzli = 1./(zl-targZ);
        
        
        
        StripsToCoor(u, v, xl, yl, stal);
        
        
        const fvec &tx = (xl - targX)*dzli;
        const 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];
        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 == kAllSecEIter) || (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.; T.C44 = MaxInvMom/3.*MaxInvMom/3.;
        
        // #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++){
            if ( ( isec != kAllPrimEIter ) && ( isec != kAllSecEIter ) ) {
#ifdef USE_RL_TABLE
                L1AddMaterial( T, fRadThick[ista].GetRadThick(T.x, T.y), MaxInvMom );
#else
                L1AddMaterial( T, vStations[ista].materialInfo, MaxInvMom );
                
#endif
                if (ista == NMvdStations - 1) L1AddPipeMaterial( T, MaxInvMom );
            }
            else {
#ifdef USE_RL_TABLE
                L1AddMaterial( T, fRadThick[ista].GetRadThick(T.x, T.y), MaxInvMom, 1, 0.000511f*0.000511f );
                
#else
                L1AddMaterial( T, vStations[ista].materialInfo, MaxInvMom, 1, 0.000511f*0.000511f );
                
#endif
                if (ista == NMvdStations - 1) L1AddPipeMaterial( T, MaxInvMom, 1, 0.000511f*0.000511f );
            }
        }
        
        
        // add left hit
        L1Filter( T, stal.frontInfo, u );
        L1Filter( T, stal.backInfo,  v );
#endif
        
        
        if ( ( isec != kAllPrimEIter ) && ( isec != kAllSecEIter ) ) {
#ifdef USE_RL_TABLE
            L1AddMaterial( T, fRadThick[istal].GetRadThick(T.x, T.y), MaxInvMom );
#else
            L1AddMaterial( T, stal.materialInfo, MaxInvMom );
#endif
            if ( (istam >= NMvdStations) && (istal <= NMvdStations - 1) )  L1AddPipeMaterial( T, MaxInvMom );
        }
        else {
#ifdef USE_RL_TABLE
            L1AddMaterial( T, fRadThick[istal].GetRadThick(T.x, T.y), MaxInvMom, 1, 0.000511f*0.000511f );
#else
            L1AddMaterial( T, stal.materialInfo, MaxInvMom, 1, 0.000511f*0.000511f );
#endif
            if ( (istam >= NMvdStations) && (istal <= NMvdStations - 1) )  L1AddPipeMaterial( T, MaxInvMom, 1, 0.000511f*0.000511f );
        }
        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.
inline 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, fvec *HitTime, fvec *Event
//                 ,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;
        L1TrackPar& T1 = T_1[i1_V];

        const 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];
        
        
        L1HitAreaTime areaTime( vGridTime[ &stam - vStations ], T1.x[i1_4]*iz, T1.y[i1_4]*iz, (sqrt(Pick_m22*(T1.C00 + stam.XYInfo.C00))+MaxDZ*fabs(T1.tx))[i1_4]*iz, (sqrt(Pick_m22*(T1.C11 + stam.XYInfo.C11))+MaxDZ*fabs(T1.ty))[i1_4]*iz, HitTime[i1_V][i1_4], 50);

        
        THitI imh = 0;
        
        while( areaTime.GetNext( imh ) ) {
            
            
            
            const L1HitPoint &hitm = vStsHits_m[imh];
            // check y-boundaries
            if (fabs(HitTime[i1_V][i1_4]-hitm.time)>35) continue;
            
#ifdef USE_EVENT_NUMBER
            if ((Event[i1_V][i1_4]!=hitm.n)) continue;
#endif USE_EVENT_NUMBER
            
            // - check whether hit belong to the window ( track position +\- errors ) -
            const fscal &zm = hitm.Z();

            
            fvec y, C11;
            L1ExtrapolateYC11Line( T1, zm, y, C11 );

            const fscal &dy_est2 = Pick_m22[i1_4] * fabs(C11[i1_4] + stam.XYInfo.C11[i1_4]);
            
            const fscal &dy = hitm.Y() - y[i1_4];

            if ( dy*dy > dy_est2 && dy < 0  ) continue;

            //   if ( ( nl != nm ) ) continue;
            
            // check x-boundaries
            fvec x, C00;
            L1ExtrapolateXC00Line( T1, zm, x, C00 );
            const fscal &dx_est2 = Pick_m22[i1_4] * fabs(C00[i1_4] +  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() );
#ifdef DO_NOT_SELECT_TRIPLETS
            if (isec!=TRACKS_FROM_TRIPLETS_ITERATION)
#endif
            if ( chi2[i1_4] > DOUBLET_CHI2_CUT ) continue;
            T1.time[i1_4] = hitm.time;
            
#ifdef USE_EVENT_NUMBER
            T1.n[i1_4] = hitm.n;

#endif USE_EVENT_NUMBER
            L1FilterChi2           ( stam.backInfo,  x, y, C00, C10, C11, chi2, hitm.V() );
            
#ifdef DO_NOT_SELECT_TRIPLETS
            if (isec!=TRACKS_FROM_TRIPLETS_ITERATION)
#endif
            if ( chi2[i1_4] > DOUBLET_CHI2_CUT ) continue;

            
            i1_2.push_back(i1);
#ifdef DOUB_PERFORMANCE
            hitsl_2.push_back(hitsl_1[i1]);
#endif // DOUB_PERFORMANCE
            hitsm_2.push_back(imh);
            
            n2++;
        }
        
    }  // for i1
    
    
    
}


/// 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
                        )
{
    
    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;

    // ---- 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];
                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];
            }  // n2_4
            
            
            
            // add middle hit
            L1ExtrapolateLine( T2, zPos_2 );
            
            L1Filter( T2, stam.frontInfo, u_front_2 );
            
            L1Filter( T2, stam.backInfo,  u_back_2 );
            
            
            if ( ( isec != kAllPrimEIter ) && ( isec != kAllSecEIter ) ) {
#ifdef USE_RL_TABLE
                L1AddMaterial( T2, fRadThick[istam].GetRadThick(T2.x, T2.y), T2.qp );
#else
                L1AddMaterial( T2, stam.materialInfo, T2.qp );
#endif
                if ( (istar >= NMvdStations) && (istam <= NMvdStations - 1) ) L1AddPipeMaterial( T2, T2.qp );
            }
            else {
#ifdef USE_RL_TABLE
                L1AddMaterial( T2, fRadThick[istam].GetRadThick(T2.x, T2.y), T2.qp, 1, 0.000511f*0.000511f );
#else
                L1AddMaterial( T2, stam.materialInfo, T2.qp, 1, 0.000511f*0.000511f );
#endif
                if ( (istar >= NMvdStations) && (istam <= NMvdStations - 1) ) L1AddPipeMaterial( T2, T2.qp, 1, 0.000511f*0.000511f );
            }
            
            // extrapolate to the right hit station
            L1Extrapolate( T2, star.z, T2.qp, f2 );
            
            // ---- Find the triplets(right hit). Reformat data in the portion of triplets. ----
            
            for (Tindex i2_4 = 0; i2_4 < n2_4; ++i2_4){

                
                if ( T2.C00[i2_4] < 0 ||  T2.C11[i2_4] < 0 ||  T2.C22[i2_4] < 0 ||  T2.C33[i2_4] < 0 ||  T2.C44[i2_4] < 0 ) continue;
                
                const fvec &Pick_r22 = (TRIPLET_CHI2_CUT - T2.chi2);
                
                
                // find first possible hit
                
#ifdef DO_NOT_SELECT_TRIPLETS
                if ( isec == TRACKS_FROM_TRIPLETS_ITERATION )
                    Pick_r22 = Pick_r2+1;
#endif
                
                
                const fscal &iz = 1/T2.z[i2_4];
                L1HitAreaTime area( vGridTime[ &star - vStations ], T2.x[i2_4]*iz, T2.y[i2_4]*iz, (sqrt(Pick_r22*(T2.C00 + stam.XYInfo.C00))+MaxDZ*fabs(T2.tx))[i2_4]*iz, (sqrt(Pick_r22*(T2.C11 + stam.XYInfo.C11))+MaxDZ*fabs(T2.ty))[i2_4]*iz, T2.time[i2_4], 50 );
                
                THitI irh = 0;
                
                while( area.GetNext( irh ) ) {

                    
                    const L1HitPoint &hitr = vStsHits_r[irh];
                    if (fabs(T2.time[i2_4]-hitr.time)>35) continue;
                    
#ifdef USE_EVENT_NUMBER
                   if ((T2.n[i2_4]!=hitr.n)) continue;
#endif USE_EVENT_NUMBER
                    const fscal &zr = hitr.Z();
                    const fscal &yr = hitr.Y();
                    
                    // - 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[i2_4]*(fabs(C11[i2_4] + 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;
                    if ( dy2 > dy_est2 && dy < 0 ) continue; // if (yr < y_minus_new[i2_4]) continue;
                    
                    //   cout<<irh<<" irh 2"<<endl;
                    
                    // 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[i2_4]*(fabs(C00[i2_4] + 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;
                    //    cout<<chi2[i2_4]<<" T2"<<endl;
                    L1FilterChi2XYC00C10C11( star.frontInfo, x, y, C00, C10, C11, chi2, hitr.U() );
                    L1FilterChi2           ( star.backInfo,  x, y, C00, C10, C11, chi2, hitr.V() );
#ifdef DO_NOT_SELECT_TRIPLETS
                    if (isec!=TRACKS_FROM_TRIPLETS_ITERATION)
#endif

                        if ( chi2[i2_4] > TRIPLET_CHI2_CUT || C00[i2_4] < 0 || C11[i2_4] < 0 ) 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);
                    
                    
                    
                    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
    }  // 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

    };
    
}


/// 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 = .1;
        T.x  = x[0];
        T.y  = y[0];
        
        fvec dz01 = 1./(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.;
        T.NDF = 2.;
        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 );
#ifdef USE_RL_TABLE
            L1AddMaterial( T, fRadThick[ista[ih]].GetRadThick(T.x, T.y), T.qp );
#else
            L1AddMaterial( T, sta[ih].materialInfo, T.qp );
#endif
            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 );
#ifdef USE_RL_TABLE
                L1AddMaterial( T, fRadThick[ista[ih]].GetRadThick(T.x, T.y), T.qp );
#else
                L1AddMaterial( T, sta[ih].materialInfo, T.qp );
#endif
                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 );
#ifdef USE_RL_TABLE
                L1AddMaterial( T, fRadThick[ista[ih]].GetRadThick(T.x, T.y), T.qp );
#else
                L1AddMaterial( T, sta[ih].materialInfo, T.qp );
#endif
                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
                       
                       
                       )
{
    
     THitI ihitl_priv = 0;
     
     unsigned int Station = 0;
     unsigned int Thread = 0;
     unsigned int Triplet = 0;
     
     unsigned int Location = 0;
     
     unsigned char level = 0;
  
    for( Tindex i3=0; i3<n3; ++i3){
        
        
        const Tindex &i3_V = i3/fvecLen;
        const Tindex &i3_4 = i3%fvecLen;
        L1TrackPar &T3 = T_3[i3_V];
        
        // select
        fscal &chi2 = T3.chi2[i3_4];

        
        Station = istal;
        Thread = omp_get_thread_num();
        Triplet = nTripletsThread[istal][omp_get_thread_num()];
        
       // PackLocation(Location, Triplet, Station, Thread);
        
        Location = Triplet+Station*100000000 +Thread*1000000;
        

        if ( !finite(chi2) || chi2 < 0 || chi2 > TRIPLET_CHI2_CUT )  continue;

        
        // 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;
        const fscal &Cqp = 5.*sqrt(fabs(T3.C44[i3_4]));
        
        const THitI &ihitl = hitsl_3[i3] + StsHitsUnusedStartIndex[istal];
        const THitI &ihitm = hitsm_3[i3] + StsHitsUnusedStartIndex[istam];
        const 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] );

        fscal &time = T3.time[i3_4];
       // int n = T3.n[i3_4];
        
        if (ihitl_priv==0||ihitl_priv!=hitsl_3[i3]) TripForHit[0][hitsl_3[i3] + StsHitsUnusedStartIndex[istal]] = Location;
        TripForHit[1][hitsl_3[i3] + StsHitsUnusedStartIndex[istal]] = Location;

        ihitl_priv=hitsl_3[i3];
        
               
        
        L1Triplet & tr1=TripletsLocal1[Station][Thread][nTripletsThread[istal][omp_get_thread_num()]];

        
        tr1.SetLevel(0);

        tr1.Set( ihitl, ihitm, ihitr,
                   istal, istam, istar,
                   0, qp, chi2, time, Cqp, 0);
        
        
        tr1.first_neighbour=0;
        tr1.last_neighbour=0;

            

            ++nstaltriplets;
            nTripletsThread[istal][omp_get_thread_num()]++;
            
            Triplet=nTripletsThread[istal][omp_get_thread_num()];
            
            Location = Triplet+Station*100000000 +Thread*1000000;
            
          //  PackLocation(Location, Triplet, Station, Thread);
            
            
            if (ihitl_priv==hitsl_3[i3]) TripForHit[1][tr1.GetLHit()] = Location;


        if (istal > (NStations - 4)) continue;
        
        if ( TripForHit[0][ihitm] == TripForHit[1][ihitm]) continue;
        
        level = 0;
 

        
        
      // cout<<TripForHit[1][ihitm]<<" TripForHit[1][ihitm] "<<TripForHit[0][ihitm]<<" TripForHit[0][ihitm] "<<endl;
        
        for (unsigned int iN=0; iN < (TripForHit[1][ihitm]-TripForHit[0][ihitm]); ++iN){
            
            Location = TripForHit[0][ihitm]+iN;

/*            
            UnPackStation(Location, Station);
            UnPackThread(Location, Thread);
            UnPackTriplet(Location, Triplet);*/
            
//             Station = Location/100000000;
//             Thread = (Location -Station*100000000)/100000;
//             Triplet = (Location- Station*100000000-Thread*100000);
            
             Station = Location/100000000;
             Thread = (Location -Station*100000000)/1000000;
             Triplet = (Location- Station*100000000-Thread*1000000);


            L1Triplet &curNeighbour = TripletsLocal1[Station][Thread][Triplet];

            
          
            const fscal &qp2 = curNeighbour.GetQp();
            fscal &Cqp2 = curNeighbour.Cqp;

            
            if ((curNeighbour.GetMHit() != ihitr) )  continue;
            
            if (tr1.first_neighbour==0) {
              tr1.first_neighbour=Location;

            }
            
            tr1.last_neighbour=Location+1;

            // calculate level
            const unsigned char &jlevel = curNeighbour.GetLevel();
            if ( level <= jlevel ) level = jlevel + 1;

            
        }
        
        tr1.SetLevel( level );
        
    }
}



/// ------------------- 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];
        fvec HitTime[Portion/fvecLen];
        fvec Event[Portion/fvecLen];
        
        /// prepare the portion of left hits data
        
        Tindex &n1 = n_g1[ip];
        
        // cout<<n1<<" n1 "<<ip<<endl;
        
        f10(  // input
            (ip - portionStopIndex[istal+1]) * Portion, n1, vStsHits_l,
            // output
            u_front, u_back, zPos,
            hitsl_1, HitTime, Event
            );
        
        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, HitTime, Event
            //             ,            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::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;
        
        nsL1::vector<L1TrackPar>::TSimd &T_3 = fT_3[omp_get_thread_num()];
        vector<THitI> &hitsl_3 = fhitsl_3[omp_get_thread_num()];
        vector<THitI> &hitsm_3 = fhitsm_3[omp_get_thread_num()];
        vector<THitI> &hitsr_3 = fhitsr_3[omp_get_thread_num()];
        nsL1::vector<fvec>::TSimd &u_front3 = fu_front3[omp_get_thread_num()];
        nsL1::vector<fvec>::TSimd &u_back3 = fu_back3[omp_get_thread_num()];
        nsL1::vector<fvec>::TSimd &z_pos3 = fz_pos3[omp_get_thread_num()];
        
        T_3.clear();
        hitsl_3.clear();
        hitsm_3.clear();
        hitsr_3.clear();
        u_front3.clear();
        u_back3.clear();
        z_pos3.clear();
        
        /// Find the triplets(right hit). Reformat data in the portion of triplets.
        // cout<<" f30 "<<endl;
        
        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
            );
        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<<" f31 "<<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,          T_3,         hitsl_3, hitsm_3, hitsr_3, 0 );
        
        
#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
        
    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

    RealIHitP=&RealIHit_v;
    RealIHitPBuf=&RealIHit_v_buf;
    vStsHitsUnused = &vStsDontUsedHits_B;      /// array of hits used on current iteration
    vector< L1StsHit > *vStsHitsUnused_buf = &vStsDontUsedHits_A; /// buffer for copy
    
    vStsHitPointsUnused = &vStsDontUsedHitsxy_B;/// array of info for hits used on current iteration
    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 EventRate = 1e7;
//     float tau = 1.e9/EventRate;  // average time between events in [ns]
//     float T0 = 0;
//     float t0 = 0;
//     vector <float> time_mc;
//     vector <float> time_reco;
//     
//     vector <float> time_mc_ev;
//     vector <float> time_reco_ev;
//     
//     
//     vector <float> time_reco_sts;
//     vector <float> time_mc_sts;
//     
//     vector <float> mass;
//     vector <int> ista;
//     vector <float> reco;
//     vector <float> mc;
//     gRandom->SetSeed(1);
//     
//     
//     for(int iev = 0; iev < 100; ++iev){
//         float deltaT = gRandom->Exp(tau);
//         T0 += deltaT;
//         //  t0 = T0 + gRandom->Gaus(0,5);
//         
//         
//         //  cout<<T0<<" "<<t0<<endl;
//         time_mc_ev.push_back(T0);
//         time_reco_ev.push_back(T0);
//     }
//     
//     float k =0;
//     float speed = 0;
//     float shift =0;
    float lasttime =0;
    
    
    {
        
        for(int ist = 0; ist < NStations; ++ist){
            
            for(THitI ih = StsHitsStartIndex[ist]; ih < StsHitsStopIndex[ist]; ++ih)
            {
                
                
                
                //    (const_cast<float&>((*vStsHits)[ih].t_mc)) = time_mc_ev[(*vStsHits)[ih].n];
                
                //        (*vStsHitsUnused)[ih].t_mc = time_mc_ev[vStsDontUsedHits_Buf[ih].n];
                
                //      (*vStsHitsUnused)[ih].t_reco = time_reco_ev[vStsDontUsedHits_Buf[ih].n];
                
//                time_mc.push_back(time_mc_ev[(*vStsHits)[ih].n]);
//                time_reco.push_back(time_reco_ev[(*vStsHits)[ih].n]);
//                
//                ista.push_back( ist);
                
                
                
                //                 cout<<(*vStsHitsUnused)[ih].t_reco<<" (*vStsHitsUnused)[ih].t_reco"<<endl;
                
                //                 if ((ist!=0)&&(fMcDataHit2[ih]!=0)) k = fMcDataHit[ih]/fMcDataHit2[ih];
                //
                //
                //                if ((ist!=0)&&(fMcDataHit2[ih]!=0)) speed = sqrt (k*k / (1.f + k*k));
                //                shift = 0;
                
                // if ((ist!=0)&&(fMcDataHit2[ih]!=0)) shift =((*vStsZPos)[(*vStsHits)[ih].iz]-vStations[0].z[0])/(speed*29.979245800f);
                
                
                //  if ((ist!=0)&&(fMcDataHit2[ih]!=0)) (const_cast<float&>((*vStsHits)[ih].t_mc)) = (*vStsHits)[ih].t_mc+shift;
                
                //if (shift>100) {cout<<"shift  " << shift << " ih " << ih << " m p " << fMcDataHit[ih] << " " << fMcDataHit2[ih]<<endl; int uuu; cin>>uuu;}
                
                //                 (*vStsHitsUnused)[ih].t_reco = (*vStsHitsUnused)[ih].t_mc+gRandom->Gaus(0,5);
                //    cout<<shift<<endl;
                    //  (const_cast<float&>((*vStsHits)[ih].t_mc)) = time_mc_ev[(*vStsHits)[ih].n];
                    //  (const_cast<float&>((*vStsHits)[ih].t_reco)) = (const_cast<float&>((*vStsHits)[ih].t_mc))+gRandom->Gaus(0,5);
                 if (lasttime<(*vStsHits)[ih].t_reco) lasttime=(*vStsHits)[ih].t_reco;
                //(const_cast<float&>((*vStsHits)[ih].t_mc)) = time_mc_ev[(*vStsHits)[ih].n]+gRandom->Gaus(0,5);
                // vStsDontUsedHitsxy_B[ih].time = (*vStsHitsUnused)[ih].t_reco;
                
                // cout<<(*vStsHits)[ih].t_reco<<endl;
                // if ((*vStsHits)[ih].t_mc<0) { const_cast<float&>((*vStsHits)[ih].t_mc)  = 0;
                
            }
            
            // reco.push_back( (*vStsHits)[ih].t_mc);
            // mc.push_back((*vStsHits)[ih].t_mc);
            
            
            //  time_reco_sts.push_back((*vStsHits)[ih].t_mc);
            // time_mc_sts.push_back((*vStsHits)[ih].t_mc);
            //   mass.push_back(fMcDataHit2[ih]);
            
            //  cout<<(*vStsHitsUnused)[ih].t_reco<<" "<<(*vStsHitsUnused)[ih].t_mc<<endl;
            
            
        }
    }
    
    
    
    
    
#ifdef XXX
    
    
    c_time.Start();
    c_timerG.Start();
    
#endif
    
    
    float yStep = 0.55*3/sqrt( nDontUsedHits ); // empirics. 0.01*sqrt(2374) ~= 0.5
    if (yStep > 0.3) yStep = 0.3;
    float xStep = yStep*3;
  //  const float hitDensity = sqrt( nDontUsedHits );
    
    
    //     float yStep = 0.7*4/hitDensity; // empirics. 0.01*sqrt(2374) ~= 0.5
    //     if (yStep > 0.3)
    yStep = 0.01;
    xStep = yStep*3;
    
    
    
    
    vStsHitsUnused = &vStsDontUsedHits_Buf;
    
    cout<<lasttime<<" lasttime"<<endl;
    
    for( int iS = 0; iS < NStations; ++iS ) {
        
        vGridTime[iS].BuildBins(-1,1,-0.6,0.6, 0, lasttime, xStep,yStep, (lasttime+1)/20);
        
        vGridTime[iS].StoreHits(StsHitsUnusedStopIndex[iS]-StsHitsUnusedStartIndex[iS], &((*vStsHits)[StsHitsUnusedStartIndex[iS]]), iS, *this, StsHitsUnusedStartIndex[iS], &(vStsDontUsedHits_Buf[StsHitsUnusedStartIndex[iS]]), &((*vStsHits)[StsHitsUnusedStartIndex[iS]]), &(RealIHit_v[StsHitsUnusedStartIndex[iS]]));
        
        
    }
    
    
    for(int ist = 0; ist < NStations; ++ist)
      for(THitI ih = StsHitsStartIndex[ist]; ih < StsHitsStopIndex[ist]; ++ih)
      {
        L1StsHit &h = *(const_cast<L1StsHit*>(&((*vStsHits)[ih])));
        SetFUnUsed( const_cast<unsigned char&>((*vSFlag)[h.f]));
        SetFUnUsed( const_cast<unsigned char&>((*vSFlagB)[h.b]));
      }
    
    
    for(int ista = 0; ista < NStations; ++ista){
#pragma omp parallel for schedule(dynamic, 5)
        for(THitI ih = StsHitsStartIndex[ista]; ih < StsHitsStopIndex[ista]; ++ih)
            CreateHitPoint(vStsDontUsedHits_Buf[ih],ista, vStsDontUsedHitsxy_B[ih]);
    }
    
    
    
#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();
    
#endif
    
    
       
TStopwatch c_time1;
c_time1.Start();   
    
    
    for (isec = 0; isec < fNFindIterations; ++isec){ // all finder
      

        
      n_g1.assign(n_g1.size(), Portion);
      
      for (int n=0; n<fNThreads; n++)        
        for (int j=0; j<8; j++)  nTripletsThread[j][n]=0;
        
        TripForHit[0].assign(TripForHit[0].size(), 0);
        TripForHit[1].assign(TripForHit[1].size(), 0);

        
        
        /// isec - number of current iterations, fNFindIterations - number of all iterations
        
#ifdef COUNTERS
        Tindex nSinglets = 0;
#endif

        
        if (isec != 0){
            vector<THitI>* RealIHitPTmp = RealIHitP;
            RealIHitP = RealIHitPBuf;
            RealIHitPBuf = RealIHitPTmp;
            
            vector< L1StsHit > *vStsHitsUnused_temp = vStsHitsUnused;
            vStsHitsUnused = vStsHitsUnused_buf;
            vStsHitsUnused_buf = vStsHitsUnused_temp;
            
            vector< L1HitPoint > *vStsHitsUnused_temp2 = vStsHitPointsUnused;
            vStsHitPointsUnused = vStsHitPointsUnused_buf;
            vStsHitPointsUnused_buf = vStsHitsUnused_temp2;
        }
        
        
        
        {
            
            
            
            
            // #pragma omp  task
            {
                
                // --- SET PARAMETERS FOR THE ITERATION ---
                
                FIRSTCASTATION = 0;
                // if ( (isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter) )
                //   FIRSTCASTATION = 2;
                
                DOUBLET_CHI2_CUT = 2.706; // prob = 0.1
                TRIPLET_CHI2_CUT = 5;
                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 = 6.252;//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 == kAllPrimEIter) || (isec == kAllPrimJumpIter) || (isec == kAllSecIter) || (isec == kAllSecEIter) || (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 = 2.0/0.5;                     // max considered q/p
                if ( (isec == kAllPrimJumpIter) || (isec == kAllSecIter) || (isec == kAllSecJumpIter) ) MaxInvMom =  1.0/0.1;
                if ( (isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllSecEIter) ) MaxInvMom = 1./0.05;
                
                MaxSlope = 1.1;
                if ( // (isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllPrimJumpIter) ||
                    (isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter) ) MaxSlope =  1.5;
                
                // define the target
                targX = 0; targY = 0; targZ = 0;      //  suppose, what target will be at (0,0,0)
                
                float SigmaTargetX = 0, SigmaTargetY = 0; // target constraint [cm]
                if ( (isec == kFastPrimIter) || (isec == kFastPrimIter2) || (isec == kFastPrimJumpIter) ||
                    (isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllPrimJumpIter) ){ // target
                    targB = vtxFieldValue;
                    if ( (isec == kFastPrimIter) || (isec == kAllPrimIter) || (isec == kAllPrimEIter) )
                        SigmaTargetX = SigmaTargetY = 1; // target
                    else
                        SigmaTargetX = SigmaTargetY = 1;
                }
                if ( (isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter) ) { //use outer radius of the 1st station as a constraint
                    L1Station &st = vStations[0];
                    SigmaTargetX = SigmaTargetY = 10;//st.Rmax[0];
                    targZ = 0.;//-1.;
                    st.fieldSlice.GetFieldValue( 0, 0, 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)
                MaxDZ = 0;
                if (  (isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllPrimJumpIter) ||
                    (isec == kAllSecIter ) || (isec == kAllSecEIter)|| (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

            
            /// 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-3; 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
        }
        
        //      c_timer.Stop();
        //      ti[isec]["init  "] = c_timer;
        //      c_timer.Start();
        
        
        ///   stage for triplets creation
        
#ifdef XXX
        TStopwatch c_timer;
        TStopwatch c_time;
        c_timer.Start();
#endif
        
        
        const Tindex vPortion = Portion/fvecLen;
        L1TrackPar T_1[vPortion];
        L1FieldRegion fld_1[vPortion];
        THitI hitsl_1[Portion];
        L1TrackPar TG_1[vPortion];
        L1FieldRegion fldG_1[vPortion];
        THitI hitslG_1[Portion];
        
        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)
        
        
        vector<THitI> hitsmG_2; /// middle hits indexed by number of doublets in portion(i2)
        vector<THitI> i1G_2; /// index in portion of singlets(i1) indexed by index in portion of doublets(i2)
        //      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
        hitsm_2.reserve(3500);
        i1_2.reserve(3500);
        
        hitsmG_2.reserve(800);
        i1G_2.reserve(800);
        
        
        
        for (int istal = NStations-3; istal >= FIRSTCASTATION; istal--){//  //start downstream chambers
          
#pragma omp parallel for firstprivate(T_1, fld_1, hitsl_1, hitsm_2, i1_2, TG_1, fldG_1, hitslG_1, hitsmG_2, i1G_2) schedule(dynamic, 2)
            for(Tindex ip = portionStopIndex[istal+1]; ip < portionStopIndex[istal]; ++ip ){
                
                // 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)
                
                Tindex n_2; /// number of doublets in portion
                
                hitsm_2.clear();
                i1_2.clear();
                
                
                
                DupletsStaPort(
                               istal, istal + 1,
                               ip,
                               n_g1, portionStopIndex,
                               
                               // output
                               T_1,
                               fld_1,
                               hitsl_1,
                               
                               //            lmDuplets[istal],
                               
                               
                               n_2,
                               i1_2,
                               hitsm_2
                               );
                
                
                Tindex 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]
                                );

                
                
                if ( (isec == kFastPrimJumpIter) || (isec == kAllPrimJumpIter) || (isec == kAllSecJumpIter) ) {
                    
                    Tindex nG_2;
                    hitsmG_2.clear();
                    i1G_2.clear();
                    
                    
                    DupletsStaPort(  // input
                                   istal, istal + 2,
                                   ip,
                                   n_g1, portionStopIndex,
                                   // output
                                   TG_1,
                                   fldG_1,
                                   hitslG_1,
                                   
                                   nG_2,
                                   i1G_2,
                                   hitsmG_2
                                   );
                    
                    
                    TripletsStaPort(  // input
                                    istal, istal + 1,  istal + 3,
                                    nstaltriplets,
                                    T_1,
                                    fld_1,
                                    hitsl_1,
                                    
                                    n_2,
                                    i1_2,
                                    hitsm_2,
                                    ip-portionStopIndex[istal+1]
                                    );
                    
                    
                    TripletsStaPort(  // input
                                    istal, istal + 2,  istal + 3,
                                    nstaltriplets,
                                    TG_1,
                                    fldG_1,
                                    hitslG_1,
                                    
                                    nG_2,
                                    i1G_2,
                                    hitsmG_2,
                                    ip-portionStopIndex[istal+1]
                                    );

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



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

        
        
        
        int min_level = 0; // min level for start triplet. So min track length = min_level+3.
        //    if (isec == kFastPrimJumpIter) min_level = 1;
        if ( (isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter) )
            min_level = 1; // only the long low momentum tracks
        
#ifdef TRACKS_FROM_TRIPLETS
        if (isec == TRACKS_FROM_TRIPLETS_ITERATION)
            min_level = 0;
#endif

        //     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
        
        
        L1Branch curr_tr;
        L1Branch new_tr[8];
        L1Branch best_tr;
        fscal curr_chi2;
        
        fscal best_chi2;
        unsigned char best_L;

        unsigned char curr_L = 1;
        int ndf = 1;

        
        // 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

            TStopwatch Time;
            
            
            //  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 
            
            vStripToTrack.assign(StsHitsStopIndex[NStations-1], -1);
            vStripToTrackB.assign(StsHitsStopIndex[NStations-1], -1);
            
            for (int i=0; i<fNThreads; ++i){
                
                numberCandidateThread [i]=0;
            }
            
            for( int istaF = FIRSTCASTATION; istaF <= NStations-3-ilev; ++istaF ){
                
#pragma omp parallel for firstprivate(curr_tr, new_tr, best_tr, curr_chi2, best_chi2, best_L, curr_L, ndf ) // schedule(dynamic, 10)
                   for( Tindex ip=0; ip<fNThreads; ++ip ){
                    
                    
                    
                    

                      for( Tindex itrip=0; itrip<nTripletsThread[istaF][ip]; ++itrip ){
                        
int thread_num = omp_get_thread_num();


                        L1Triplet &first_trip = (TripletsLocal1[istaF][ip][itrip]);

                        
                        // ghost supression !!!
#ifndef FIND_GAPED_TRACKS
                        if( /*(isec == kFastPrimIter) ||*/ (isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter) ) {
#else
                            if( (isec == kFastPrimIter) || (isec == kFastPrimIter2) || (isec == kFastPrimJumpIter) ||
                               (isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllPrimJumpIter) ||
                               (isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter) ) {
#endif
                                { // ghost supression !!!
                                    if ( isec != kFastPrimIter && isec != kAllPrimIter && isec != kAllPrimEIter && isec != kAllSecEIter )
                                        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 ( (ilev == 0) &&
                                        (GetFStation((*vSFlag)[(*vStsHitsUnused)[first_trip.GetLHit()].f]) != 0) ) continue;  // ghost supression // collect only MAPS tracks-triplets
                                }
                                if ( first_trip.GetLevel() < ilev ) continue; // try only triplets, which can start track with ilev+3 length. w\o it have more ghosts, but efficiency either
                            }
                            
                            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[0]=((*RealIHitP)[first_trip.GetLHit()]);
                            

                            (curr_tr).NHits=1;
                            
                            
                            curr_L = 1;
                            curr_chi2 = first_trip.GetChi2();
                            
                            best_tr = (curr_tr);
                            
                            best_chi2 = curr_chi2;
                            best_L = curr_L;
                            
                            
                            CAFindTrack(istaF, best_tr, best_L, best_chi2, &first_trip, (curr_tr), curr_L, curr_chi2, min_best_l, new_tr); /// reqursive func to build a tree of possible track-candidates and choose the best
                            
                            //  if ( best_L < min_best_l ) continue;
                            if ( best_L < ilev + 2 ) continue;      // lose maximum one hit
                            if ( best_L < min_level + 3 ) continue; // should find all hits for min_level
                            
                            ndf = best_L*2-5;
                            best_chi2 = best_chi2/ndf; //normalize
                            
#ifndef TRACKS_FROM_TRIPLETS
                            if( fGhostSuppression ){
                                if( best_L == 3 ){
                                    
                                    if( ( (isec == kAllSecIter) || (isec == kAllSecJumpIter) ) && (istaF != 0) ) continue;
                                    if( (isec != kAllSecIter) && (isec != kAllSecJumpIter) && (best_chi2 > 5.0) ) continue;
                                }
                            }
#endif
                            
                            best_tr.Set( istaF, best_L, best_chi2, first_trip.GetQpOrig());
                            
                            
                            best_tr.CandIndex=numberCandidateThread [thread_num]+thread_num*100000;
                            
                            CandidatesTrack[thread_num][ numberCandidateThread [thread_num]]=best_tr;

                            
                            
                            L1Branch &tr = CandidatesTrack[thread_num][ numberCandidateThread [thread_num]];
                            


                            bool check=1;
                            for (vector<THitI>::iterator phitIt = tr.StsHits.begin();/// used strips are marked
                                 phitIt != tr.StsHits.begin()+tr.NHits; ++phitIt){
                                
                                const L1StsHit &h = (*vStsHits)[*phitIt];
                                
                                  omp_set_lock(&hitToBestTrackB[h.b]);
                                
                                int &strip1 = (vStripToTrackB)[h.b];
                                
                                if (strip1!=-1)
                                {
                                    
                                    int thread = strip1/100000;
                                    int num = (strip1 -thread*100000);
                                    
                                    if (L1Branch::compareCand(tr, CandidatesTrack[thread][num])) {
                                      CandidatesTrack[thread][num].CandIndex=-1;
                                      strip1=tr.CandIndex;
                                      
                                    }
                                    
                                    else {check=0; omp_unset_lock(&hitToBestTrackB[h.b]); break;}
                                    
                                }
                                
                                else strip1 = tr.CandIndex;
                                
                                
                                  omp_unset_lock(&hitToBestTrackB[h.b]);
                                
                               if (check) {omp_set_lock(&hitToBestTrackF[h.f]);
                                
                                
                                int &strip2 = (vStripToTrack)[h.f];
                                
                                if (strip2!=-1)
                                {
                                    
                                    
                                    int thread = strip2/100000;
                                    int num = (strip2 -thread*100000);
                                    
                                    if (L1Branch::compareCand(tr, CandidatesTrack[thread][num])){
                                      CandidatesTrack[thread][num].CandIndex=-1;
                                      strip2=tr.CandIndex;  
                                      
                                    }
                                    else {check=0; omp_unset_lock(&hitToBestTrackF[h.f]); break;}
                                    
                                }
                                
                                else strip2 = tr.CandIndex;
                                
                                omp_unset_lock(&hitToBestTrackF[h.f]);}
                            
                                
                            }
                            
                            if (check) numberCandidateThread [thread_num]++;
                            
                        }  // itrip
                        
                        
                    }
                    
                }

                
                
                
                if (--nlevel == 0) break;

                
                for (int i=0; i<fNThreads; ++i){
                    
                    SavedCand[i]=0;
                    SavedHits[i]=0;
                    
                }
                
                
                
                
                for (int i=0; i<fNThreads; ++i){
                    L1Track t;
                    
                            
                    
#pragma omp parallel for  schedule(dynamic, 10)  firstprivate(t)
                    for ( Tindex iCandidate = 0; iCandidate < numberCandidateThread [i]; ++iCandidate ){
                        
                        
                        L1Branch &tr = CandidatesTrack[i][iCandidate];
                        
                        bool check=1;
                        
                        if (CandidatesTrack[i][iCandidate].CandIndex!=-1){
                          
                        for (vector<THitI>::iterator phIt = tr.StsHits.begin();/// used strips are marked
                             phIt != tr.StsHits.begin()+tr.NHits; ++phIt)
                        {
                            
                            const L1StsHit &h = ((*vStsHits)[*phIt]);
                            
                            if (((vStripToTrackB)[h.b]!=tr.CandIndex)||((vStripToTrack)[h.f]!=tr.CandIndex)) {check=0; break;}
                            
                        }
                        
                        
                        if (check)  {
                            
                            
#ifdef EXTEND_TRACKS
                 if (tr.NHits!=8) BranchExtender(tr);
#endif

                            float  sumTime = 0;
                            
                            
                            for (vector<THitI>::iterator phIt = tr.StsHits.begin();/// used strips are marked
                                 phIt != tr.StsHits.begin()+tr.NHits; ++phIt)
                            {
                                L1StsHit &h = *(const_cast<L1StsHit*>(&((*vStsHits)[*phIt])));
                                
                                
                                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()]++;
                                
                                const   L1StsHit &hit = (*vStsHits)[*phIt];
                                
                                
                                L1HitPoint tempPoint = CreateHitPoint( hit, 0); //TODO take number of station from hit
                                
                                float xcoor = tempPoint.X();
                                float ycoor = tempPoint.Y();
                                float zcoor = tempPoint.Z();
                                float timeFlight = sqrt(xcoor*xcoor+ycoor*ycoor+zcoor*zcoor)/30.f; // c = 30[cm/ns]
                                sumTime += (hit.t_reco - timeFlight);
                            }

                            t.NHits = tr.NHits;
                            t.Momentum = tr.Momentum;
                            t.fTrackTime = sumTime/t.NHits;

                            vTracks_local [omp_get_thread_num()][SavedCand[omp_get_thread_num()]]=(t);
                            SavedCand[omp_get_thread_num()]++;

                        }
                      }   
                        
                    }
                    
                    
                    
                }
                
                
                
                
                Time.Stop();
                ti[isec]["table"] += Time;
                
                Time.Start();
                
                
                
                
                
                
                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]);
                    }
                    
                }
            } //istaf
            
            
            c_timer.Stop();
            ti[isec]["tracks"] = c_timer;
            
            
            
                        c_timer.Start();
            
            
            if (isec<(fNFindIterations-1))
                
            {

                int NHitsOnStation=0;
                
                for(int ista = 0; ista < NStations; ++ista){
                    
                    int Nelements = StsHitsUnusedStopIndex[ista]-StsHitsUnusedStartIndex[ista];
                    int NHitsOnStationTmp=NHitsOnStation;
                    vGridTime[ista].UpdateIterGrid(Nelements, &((*vStsHitsUnused)[StsHitsUnusedStartIndex[ista]]),
                                                   RealIHitPBuf, &((*RealIHitP)[StsHitsUnusedStartIndex[ista]]),
                                                   &RealIHit_v_buf2, vStsHitsUnused_buf, vStsHitPointsUnused_buf,
                                                   &((*vStsHitPointsUnused)[StsHitsUnusedStartIndex[ista]]), NHitsOnStation, ista, *this, vSFlag, vSFlagB );
                    StsHitsUnusedStartIndex[ista] = NHitsOnStationTmp;
                    StsHitsUnusedStopIndex[ista] = NHitsOnStation;
                    
                }
                
#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();
#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,
                                    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
            const THitI &ihitm = curr_trip->GetMHit();
            
            
            if(!GetFUsed( (*vSFlag)[(*vStsHitsUnused)[ihitm].f] | (*vSFlagB)[(*vStsHitsUnused)[ihitm].b] )  ){
                
                //        curr_tr.StsHits.push_back((*RealIHitP)[ihitm]);
                
                curr_tr.StsHits[curr_tr.NHits]=((*RealIHitP)[ihitm]);
                
                curr_tr.NHits++;
                
                curr_L++;
            }
            const THitI &ihitr = curr_trip->GetRHit();
            if( !GetFUsed( (*vSFlag)[(*vStsHitsUnused)[ihitr].f] | (*vSFlagB)[(*vStsHitsUnused)[ihitr].b] ) ){
                
                //curr_tr.StsHits.push_back((*RealIHitP)[ihitr]);
                curr_tr.StsHits[curr_tr.NHits]=((*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   = curr_tr;
                best_chi2 = curr_chi2;
                best_L    = curr_L;
            }
            
            return;
            
        } else { //MEANS level ! = 0
            
            
            
            // L1Triplet* new_trip = curr_trip->neibourBest;
            // try to extend. try all possible triplets
            
            
          //  for (Tindex in = 0; in < (curr_trip->neighbours2.Size()); in++) {
                
     unsigned int Station = 0;
     unsigned int Thread = 0;
     unsigned int Triplet = 0; 
     unsigned int Location = 0;  
          
          
            for (Tindex in = 0; in < (curr_trip->last_neighbour-curr_trip->first_neighbour); in++) {  

                
                Location = curr_trip->first_neighbour+in;
                
//                 const fscal &qp2 = curr_trip->GetQp();
//                 fscal &Cqp2 = curr_trip->Cqp;
            
        //    if (( fabs(qp - qp2) > PickNeighbour * (Cqp + Cqp2) ) )  continue;
                
//                 UnPackStation(Location, Station);
//                 UnPackThread(Location, Thread);
//                 UnPackTriplet(Location, Triplet);
                
                Station = Location/100000000;
                Thread = (Location -Station*100000000)/1000000;
                Triplet = (Location- Station*100000000-Thread*1000000);

                
                const L1Triplet &new_trip = TripletsLocal1[Station][Thread][Triplet];
                
                const fscal &qp1 = curr_trip->GetQp();
                const fscal &qp2 = new_trip.GetQp();
                fscal dqp = fabs(qp1 - qp2);
                fscal Cqp = curr_trip->Cqp;
                Cqp      += new_trip.Cqp;
                
                if  (curr_trip->GetLevel() != new_trip.GetLevel() + 1) continue;
                
                if ( dqp > PickNeighbour * Cqp  ) continue; // bad neighbour // CHECKME why do we need recheck it?? (it really change result)
                
                
                
                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].resize(curr_tr);
                    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].StsHits.push_back((*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 = fabs(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;
                    
                    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, new_tr);
                    
                } // add triplet to track
                
                // break;
                
            } // for neighbours
        } // level = 0
    }