/*
 *=====================================================
 *  
 *  CBM Level 1 Reconstruction 
 *   
 *  Authors: I.Kisel,  S.Gorbunov
 *
 *  e-mail : ikisel@kip.uni-heidelberg.de 
 *
 *=====================================================
 *
 *  Finds tracks using the Cellular Automaton algorithm
 *
 */
// iklm v8.2.1.end

// #define XXX  
// #define TBB 




		

	
#include "TStopwatch.h" 

/* there are next parameters in next functions:
f11(...):
	fvec vINF = 10.0;
	fvec qp1 = MaxInvMom; // used for scattering on stantion, then we don't know momentum
	if( isec>=2 ) qp1 = 1;
	fvec Infqp = MaxInvMom/3.*MaxInvMom/3.; // used as C44, then we don't know momentum

	
*/



#include "L1Algo.h"
#include "L1TrackPar.h"
#include "L1Branch.h"
#include "L1Extrapolation.h"
#include "L1Filtration.h"
#include "L1AddMaterial.h"
#include "L1HitPoint.h"

#include "L1Portion.h"

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

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



// for f
// 	unsigned char inline GetFStation( unsigned char flag ){ return flag/4; }
//   bool inline GetFUsed   ( unsigned char flag ){ return (flag&0x02)!=0; }
//   bool inline GetFUsedD  ( unsigned char flag ){ return (flag&0x01)!=0; }
// 	
// 	void inline SetFStation ( unsigned char &flag, unsigned iStation ){ flag = iStation*4 + (flag%4); }
//   void inline SetFUsed    ( unsigned char &flag ){ flag |= 0x02; }
//   void inline SetFUsedD   ( unsigned char &flag ){ flag |= 0x01; }
//   void inline SetFUnUsed  ( unsigned char &flag ){ flag &= 0xFC; }
//   void inline SetFUnUsedD ( unsigned char &flag ){ flag &= 0xFE; }

	
					/// Prepare the portion of left hits data
inline void L1Algo::f10(	// input
                int start_lh, int n1,  L1HitPoint *vStsHits_l, 
									// output
								fvec *u_front, fvec *u_back, 
								vector<unsigned /*short*/ int> &hitsl_1
								)
{
	int i1 = 0, i1_V = 0, i1_4 = 0;
	
	int end_lh = start_lh+n1;
#ifdef HAVE_VC_SSE
// cout << "b1" << endl;
	fvec::Memory temp_front,temp_back;
	for (int ilh = start_lh; ilh < end_lh; ){
		for (int i=0; ((i < 4) && (ilh < end_lh)); i++,ilh++){

			L1HitPoint &hitl = vStsHits_l[ilh];

// 			hitsl_1[i1] = ilh;
      hitsl_1.push_back(ilh);
			temp_front[i] = hitl.x;
			temp_back[i] = hitl.v;

			i1++;
		}

		fvec vtempf(temp_front);
		fvec vtempb(temp_back);
		u_front[i1_V] = vtempf;
		u_back[i1_V]  = vtempb;
		i1_V++;
	}

	int n1_4 = n1%fvecLen;
	if( n1_4 > 0 ){
		i1_V--;
		for( int ii1_4=n1_4-1; n1_4<fvecLen; n1_4++ ){
			temp_front[n1_4] = temp_front[ii1_4];
			temp_back[n1_4] = temp_back[ii1_4];
		}	
		fvec vtempf(temp_front);
		fvec vtempb(temp_back);
		u_front[i1_V] = vtempf;
		u_back[i1_V] = vtempb;
	}
// cout << "e1" << endl;
#else // HAVE_VC_SSE	
	for (int ilh = start_lh; ilh < end_lh; ilh++){

		L1HitPoint &hitl = vStsHits_l[ilh];

//    hitsl_1[i1] = ilh;
    hitsl_1.push_back(ilh);
		u_front[i1_V][i1_4] = hitl.x;
		u_back [i1_V][i1_4] = hitl.v;

		i1++;
		i1_V= i1/fvecLen;
		i1_4= i1%fvecLen;
	}

	if( i1_4>0 ){
		for( int ii1_4=i1_4-1; i1_4<fvecLen; i1_4++ ){
			u_front[i1_V][i1_4] = u_front[i1_V][ii1_4];
			u_back [i1_V][i1_4] = u_back [i1_V][ii1_4];
		}	
	}
#endif // HAVE_VC_SSE	

// 	if(0){ // check
// 		for (int i=0; i<n1/fvecLen; i++){
// 			cout << "front = " << u_front[i] << "\t back = " << u_back[i] << endl;
// 		}
// 	int kk; std::cin >> kk;
// 	}

}


					/// Get the field approximation. Add the target to parameters estimation. Propagaete to middle station.
inline void L1Algo::f11(	// input
								int isec, int istal, L1Station *vStations,
								int n1_V,
                double MaxInvMom, double Pick_m,
								fvec targX, fvec targY, fvec targZ, L1FieldValue targB, L1XYMeasurementInfo TargetXYInfo,
                fvec *u_front, fvec *u_back, 								
									// output
// 								L1TrackPar *T_1,
                nsL1::vector<L1TrackPar>::TSimd &T_1,
                vector<L1FieldRegion> &fld_1,
								fvec* x_minusV, fvec* x_plusV, fvec* y_minusV, fvec* y_plusV
							 )
{	
	fvec vINF = 10.0;
		
	fvec qp1 = MaxInvMom; // used for scattering on stantion, then we don't know momentum
	if( isec>=2 ) qp1 = 1;

	fvec Infqp = MaxInvMom/3.*MaxInvMom/3.; // used as C44, then we don't know momentum
		
	L1Station &stal = vStations[istal];
	L1Station &stam = vStations[istal+1];
	
	fvec zl = stal.z[0];
	fvec dzli = 1./(zl-targZ[0]);	
	for( int i1_V=0; i1_V<n1_V; i1_V++ ){
	
		fvec &x = u_front[i1_V];
		fvec &v = u_back [i1_V];
		fvec y  = stal.yInfo.cos_phi*x + stal.yInfo.sin_phi*v;

		fvec tx = (x-targX)*dzli;
		fvec ty = (y-targY)*dzli;
		
		L1FieldRegion fld; // by  hit, targ and centr
    L1FieldRegion fld_1new; // by  hit, targ and centr

				// find the fieldRegions
		L1FieldValue l_B, m_B _fvecalignment;					
		stal.fieldSlice.GetFieldValue( x, y, l_B );
		stam.fieldSlice.GetFieldValue( tx*stam.z, ty*stam.z, m_B );
		if( istal>0 ){
			L1FieldValue centB;
			L1Station &st = vStations[istal/2];
			st.fieldSlice.GetFieldValue( tx*st.z, ty*st.z, centB );
			fld.Set( l_B, stal.z, centB, st.z, targB, targZ );
			fld_1new.Set( m_B, stam.z, l_B, stal.z, centB, st.z  );
		}
		else{
			fld.Set( l_B, stal.z, targB, targZ );
			fld_1new.Set( m_B, stam.z, l_B, stal.z, targB, targZ  );
		}
		fld_1.push_back(fld_1new);
	
// 		L1TrackPar &T = T_1[i1_V];
    L1TrackPar T;
	
		T.chi2 = 0.;
		T.NDF = 2.;
		if( isec==2) T.NDF = 0;
		T.x  = x;
		T.y  = y;
		T.tx = tx;
		T.ty = ty;
		T.qp = 0.;
		T.z  = zl;
		T.C00 = stal.XYInfo.C00;
		T.C10 = stal.XYInfo.C10;
		T.C11 = stal.XYInfo.C11;
		T.C22 = T.C33 = vINF; T.C44 = Infqp;
		
		//add the target
		{
			fvec eX, eY, J04, J14;
			fvec dz; 
			dz = targZ - zl;
			L1ExtrapolateJXY0(T.tx, T.ty, dz, fld, 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 );
		}
	
//----
		
		L1AddMaterial( T, stal.materialInfo, qp1 );	  
		L1Extrapolate0( T, stam.z, fld );
	
		fvec dxm_est = Pick_m*sqrt(fabs(T.C00+stam.XYInfo.C00));
		fvec dym_est = Pick_m*sqrt(fabs(T.C11+stam.XYInfo.C11));
	
		x_minusV[i1_V] = T.x - dxm_est;
		x_plusV [i1_V] = T.x + dxm_est;
		y_minusV[i1_V] = T.y - dym_est;
		y_plusV [i1_V] = T.y + dym_est;
		
//----	
		
		T_1.push_back(T);
	}// i1_V
}				

	
					/// Find the doublets. Reformat data in the portion of doublets.
inline void L1Algo::f20(	// input
								int n1, int istar, int NStations, L1Station &stal, L1Station &stam,
                L1HitPoint *vStsHits_l, L1HitPoint *vStsHits_m, int NHits_m,
								fscal *y_minus, fscal *x_minus, fscal *y_plus, fscal *x_plus,								
								nsL1::vector<L1TrackPar>::TSimd &T_1, vector<L1FieldRegion> &fld_1,
								vector<unsigned /*short*/ int> &hitsl_1,
								map<unsigned /*short*/ int, unsigned /*short*/ int> &mrDuplets_start,
									// output
								int &n2,
                vector<unsigned /*short*/ int> &i1_2,
                int &start_mhit,
                vector<unsigned /*short*/ int> &hitsm_2,
								map<unsigned /*short*/ int, unsigned /*short*/ int> &lmDuplets_start, vector<unsigned /*short*/ int> &lmDuplets_hits, unsigned int &nDuplets_lm,
								std::vector<unsigned char> &vSFlag, std::vector<unsigned char> &vSFlagB
								)
{
	
	n2 = 0;
	for (int i1 = 0; i1 < n1; i1++){
// 		cout << lmDuplets_hits.size() << " " <<  nDuplets_lm << endl; // dbg 12
		lmDuplets_start[hitsl_1[i1]] = nDuplets_lm;	

    unsigned /*short*/ int nl = vStsHits_l[hitsl_1[i1]].n;
			// find first need mhit
		for (; start_mhit < NHits_m; start_mhit++){	    // 1.0/1000 sec
			fscal ym = vStsHits_m[start_mhit].y;
			if ( ym >= y_minus[i1]) break;
		}

		for (int imh = start_mhit; imh < NHits_m; imh++){	 // 6.3/100 sec

// 			if( n2 >= MaxPortionDoublets ) break; // dbg 0
			
			L1HitPoint &hitm = vStsHits_m[imh];
			fscal xm = hitm.x;
			fscal ym = hitm.y;
      unsigned /*short*/ int nm =  hitm.n;
			
			if ( ym > y_plus[i1] ) break;
      if ( (xm < x_minus[i1] ) || (xm > x_plus[i1]) || ( nl != nm)) continue; // dbg 1
	
// 			lmDuplets_hits[nDuplets_lm++]=imh;
      lmDuplets_hits.push_back(imh);
      nDuplets_lm++;
			
// 			i1_2[n2] = i1;		
//       hitsm_2[n2] = imh;
      i1_2.push_back(i1);
      hitsm_2.push_back(imh);
						
			n2++; 
		}
		
		lmDuplets_start[hitsl_1[i1]+1] = nDuplets_lm; // mark the end

	}//i1
}

					
					/// Add the middle hits to parameters estimation. Propagate to right station.
/*void inline f21(	// input
								int n2_V, 
								L1Station &stam, L1Station &star,
								fvec *u_front, fvec *u_back,
								L1FieldRegion *fld_2,
								double Pick_r,
									// output
								L1TrackPar *T_2,
								fvec* x_minusV, fvec* x_plusV, fvec* y_minusV, fvec* y_plusV
							 )
{				
	for( int i2_V=0; i2_V<n2_V; i2_V++){
	
		L1TrackPar &T2 = T_2[i2_V];
		L1FieldRegion &f2 = fld_2[i2_V];
		L1Filter( T2, stam.frontInfo, u_front[i2_V] );
		L1Filter( T2, stam.backInfo,  u_back[i2_V] );
		
		L1AddMaterial( T2, stam.materialInfo, T2.qp );
		
		L1Extrapolate( T2, star.z, T2.qp, f2 );	
		fvec dxm_est = Pick_r*sqrt(fabs(T2.C00+star.XYInfo.C00));
		fvec dym_est = Pick_r*sqrt(fabs(T2.C11+star.XYInfo.C11));	
		x_minusV[i2_V] = T2.x - dxm_est;
		x_plusV [i2_V] = T2.x + dxm_est;	  
		y_minusV[i2_V] = T2.y - dym_est;
		y_plusV [i2_V] = T2.y + dym_est;
		
	}
} 
*/
					
					/// 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 istar, int NStations, int n1, 
								L1HitPoint *vStsHits_m,
								nsL1::vector<L1TrackPar>::TSimd &T_1, vector<L1FieldRegion> &fld_1,
                vector<unsigned /*short*/ int> &hitsl_1,
								map<unsigned /*short*/ int, unsigned /*short*/ int> &lmDuplets_start, vector<unsigned /*short*/ int> &lmDuplets_hits,
								double Pick_r,

								int n2,
                vector<unsigned /*short*/ int> &hitsm_2,
                vector<unsigned /*short*/ int> &i1_2,
																
								map<unsigned /*short*/ int, unsigned /*short*/ int> &mrDuplets_start, vector<unsigned /*short*/ int> &mrDuplets_hits,

									// output
								int &n3,
                nsL1::vector<L1TrackPar>::TSimd &T_3,
                vector<unsigned /*short*/ int> &hitsl_3,  vector<unsigned /*short*/ int> &hitsm_3,  vector<unsigned /*short*/ int> &hitsr_3,
                nsL1::vector_fvec &u_front_3, nsL1::vector_fvec &u_back_3
								)
{
  nsL1::vector<L1TrackPar>::TSimd T_2;
  vector<L1FieldRegion> fld_2;
  vector<unsigned /*short*/ int> hitsl_2;
  fvec fvec_0 _fvecalignment;
  L1TrackPar L1TrackPar_0 _fvecalignment;
  nsL1::vector_fvec u_front_2;
  nsL1::vector_fvec u_back_2;
  nsL1::vector_fvec x_minusV_2;
  nsL1::vector_fvec x_plusV_2;
  nsL1::vector_fvec y_minusV_2;
  nsL1::vector_fvec y_plusV_2;

  T_2.reserve(MaxPortionDoublets/fvecLen);
  fld_2.reserve(MaxPortionDoublets/fvecLen);
  hitsl_2.reserve(MaxPortionDoublets);
  u_front_2.reserve(MaxPortionDoublets/fvecLen);
  u_back_2.reserve(MaxPortionDoublets/fvecLen);
  x_minusV_2.reserve(MaxPortionDoublets/fvecLen);
  x_plusV_2.reserve(MaxPortionDoublets/fvecLen);
  y_minusV_2.reserve(MaxPortionDoublets/fvecLen);
  y_plusV_2.reserve(MaxPortionDoublets/fvecLen);
  


		/// f20.2
	
	int n2_V = n2/fvecLen;
	if (istar < NStations){
		for (int i2_V = 0; i2_V < n2_V; i2_V++){
		
// 			L1TrackPar &T2 = T_2[i2_V];
// 			L1FieldRegion &f2 = fld_2[i2_V];
      L1TrackPar T2;
      L1FieldRegion f2;

      u_front_2.push_back(fvec_0);
      u_back_2.push_back(fvec_0);
			for (int i2_4 = 0; i2_4 < fvecLen; i2_4++){	
				
				int i2 = i2_V*fvecLen+i2_4;
				int i1 = i1_2[i2];
				int i1_V = i1/fvecLen;
				int i1_4 = i1%fvecLen;
				
				L1TrackPar &T1 = T_1[i1_V];
				L1FieldRegion &f1 = fld_1[i1_V];
					
				int imh = hitsm_2[i2];
// 		cout << i1 << " " << i2 << " " << imh << endl;
				
				L1HitPoint &hitm = vStsHits_m[imh];
				fscal &xm = hitm.x;
				fscal &vm = hitm.v;	
					//-- copy 
#ifdef HAVE_VC_SSE
				fscal *temp;
				{			

					temp = (fscal *)(&T2.x);
					temp[i2_4] = T1.x[i1_4];
					temp[i2_4+4] = T1.y[i1_4];
					temp[i2_4+8] = T1.tx[i1_4];
					temp[i2_4+12] = T1.ty[i1_4];
					temp[i2_4+16] = T1.qp[i1_4];
					temp[i2_4+20] = T1.z[i1_4];
					temp[i2_4+24] = T1.C00[i1_4];
					temp[i2_4+28] = T1.C10[i1_4];
					temp[i2_4+32] = T1.C11[i1_4];
					temp[i2_4+36] = T1.C20[i1_4];
					temp[i2_4+40] = T1.C21[i1_4];
					temp[i2_4+44] = T1.C22[i1_4];
					temp[i2_4+48] = T1.C30[i1_4];
					temp[i2_4+52] = T1.C31[i1_4];
					temp[i2_4+56] = T1.C32[i1_4];
					temp[i2_4+60] = T1.C33[i1_4];
					temp[i2_4+64] = T1.C40[i1_4];
					temp[i2_4+68] = T1.C41[i1_4];
					temp[i2_4+72] = T1.C42[i1_4];
					temp[i2_4+76] = T1.C43[i1_4];
					temp[i2_4+80] = T1.C44[i1_4];
					temp[i2_4+84] = T1.chi2[i1_4];
					temp[i2_4+88] = T1.NDF[i1_4];

					temp = (fscal *)(&f2.cx0);
					temp[i2_4] = f1.cx0[i1_4];
					temp[i2_4+4] = f1.cx1[i1_4];
					temp[i2_4+8] = f1.cx2[i1_4];
					temp[i2_4+12] = f1.cy0[i1_4];
					temp[i2_4+16] = f1.cy1[i1_4];
					temp[i2_4+20] = f1.cy2[i1_4];
					temp[i2_4+24] = f1.cz0[i1_4];
					temp[i2_4+28] = f1.cz1[i1_4];
					temp[i2_4+32] = f1.cz2[i1_4];
					temp[i2_4+36] = f1.z0[i1_4];
				}
				temp = (fscal *)(&u_front_2[i2_V]);
				temp[i2_4] = xm;
				temp = (fscal *)(&u_back_2[i2_V]);
				temp[i2_4] = vm;
#else
				{														//	  1.83/100 sec
					T2.x[i2_4] = T1.x[i1_4];
					T2.y[i2_4] = T1.y[i1_4];
					T2.tx[i2_4] = T1.tx[i1_4];
					T2.ty[i2_4] = T1.ty[i1_4];
					T2.qp[i2_4] = T1.qp[i1_4];
					T2.z[i2_4] = T1.z[i1_4];
					T2.C00[i2_4] = T1.C00[i1_4];
					T2.C10[i2_4] = T1.C10[i1_4];
					T2.C11[i2_4] = T1.C11[i1_4];
					T2.C20[i2_4] = T1.C20[i1_4];
					T2.C21[i2_4] = T1.C21[i1_4];
					T2.C22[i2_4] = T1.C22[i1_4];
					T2.C30[i2_4] = T1.C30[i1_4];
					T2.C31[i2_4] = T1.C31[i1_4];
					T2.C32[i2_4] = T1.C32[i1_4];
					T2.C33[i2_4] = T1.C33[i1_4];
					T2.C40[i2_4] = T1.C40[i1_4];
					T2.C41[i2_4] = T1.C41[i1_4];
					T2.C42[i2_4] = T1.C42[i1_4];
					T2.C43[i2_4] = T1.C43[i1_4];
					T2.C44[i2_4] = T1.C44[i1_4];
					T2.chi2[i2_4] = T1.chi2[i1_4];
					T2.NDF[i2_4] = T1.NDF[i1_4];
	
					f2.cx0[i2_4] = f1.cx0[i1_4];
					f2.cx1[i2_4] = f1.cx1[i1_4];
					f2.cx2[i2_4] = f1.cx2[i1_4];
					f2.cy0[i2_4] = f1.cy0[i1_4];
					f2.cy1[i2_4] = f1.cy1[i1_4];
					f2.cy2[i2_4] = f1.cy2[i1_4];
					f2.cz0[i2_4] = f1.cz0[i1_4];
					f2.cz1[i2_4] = f1.cz1[i1_4];
					f2.cz2[i2_4] = f1.cz2[i1_4];
					f2.z0[i2_4] = f1.z0[i1_4];

				}
	
				u_front_2[i2_V][i2_4] = xm;
				u_back_2 [i2_V][i2_4] = vm;
#endif // HAVE_VC_SSE		
// 				if(0){ //check
// 					cout << T2.x << endl;
// 					cout << T2.y << endl;
// 					cout << T2.NDF << endl;
// 					cout << f2.cx0 << endl;
// 					cout << f2.cx1 << endl;
// 					cout << f2.cx2 << endl;
// 					cout << f2.cy0 << endl;
// 					cout << f2.cy1 << endl;
// 					cout << f2.cy2 << endl;
// 					cout << f2.cz0 << endl;
// 					cout << f2.cz1 << endl;
// 					cout << f2.cz2 << endl;
// 					cout << f2.z0 << endl;
// 
// 					cout << u_front_2[i2_V] << endl;
// 					cout << u_back_2 [i2_V] << endl;
// 					
// 					cout  << endl;
// 					int kk; std::cin >> kk;
// 				}
	
// 				hitsl_2[i2] = hitsl_1[i1];
        hitsl_2.push_back(hitsl_1[i1]);

			}  // i2_4
			
			L1Filter( T2, stam.frontInfo, u_front_2[i2_V] );
			L1Filter( T2, stam.backInfo,  u_back_2[i2_V] );
			
			L1AddMaterial( T2, stam.materialInfo, T2.qp );
			
			L1Extrapolate( T2, star.z, T2.qp, f2 );	
			fvec dxm_est = Pick_r*sqrt(fabs(T2.C00+star.XYInfo.C00));
			fvec dym_est = Pick_r*sqrt(fabs(T2.C11+star.XYInfo.C11));	
// 			x_minusV_2[i2_V] = T2.x - dxm_est;
// 			x_plusV_2 [i2_V] = T2.x + dxm_est;
// 			y_minusV_2[i2_V] = T2.y - dym_est;
// 			y_plusV_2 [i2_V] = T2.y + dym_est;
      x_minusV_2.push_back(T2.x - dxm_est);
      x_plusV_2 .push_back(T2.x + dxm_est);
      y_minusV_2.push_back(T2.y - dym_est);
      y_plusV_2 .push_back(T2.y + dym_est);

      T_2.push_back(T2);
      fld_2.push_back(f2);

		}	 // i2_V 

		int n2_4 = n2%fvecLen;
		if (n2_4 != 0){
		
// 			L1TrackPar &T2 = T_2[n2_V];
// 			L1FieldRegion &f2 = fld_2[n2_V];	
      L1TrackPar T2;
      L1FieldRegion f2;

      u_front_2.push_back(fvec_0);
      u_back_2.push_back(fvec_0);
			for (int i2_4 = 0; i2_4 < n2_4; i2_4++){

				int i2_V = n2_V;
				int i2 = i2_V*fvecLen+i2_4;			
				int i1 = i1_2[i2];
				int i1_V = i1/fvecLen;
				int i1_4 = i1%fvecLen;
				
				L1TrackPar &T1 = T_1[i1_V];
				L1FieldRegion &f1 = fld_1[i1_V];
					
				int imh = hitsm_2[i2];
				
				L1HitPoint &hitm = vStsHits_m[imh];
				fscal &xm = hitm.x;
				fscal &vm = hitm.v;	
					//-- copy 

#ifdef HAVE_VC_SSE
				fscal *temp;
				{			

					temp = (fscal *)(&T2.x);
					temp[i2_4] = T1.x[i1_4];
					temp[i2_4+4] = T1.y[i1_4];
					temp[i2_4+8] = T1.tx[i1_4];
					temp[i2_4+12] = T1.ty[i1_4];
					temp[i2_4+16] = T1.qp[i1_4];
					temp[i2_4+20] = T1.z[i1_4];
					temp[i2_4+24] = T1.C00[i1_4];
					temp[i2_4+28] = T1.C10[i1_4];
					temp[i2_4+32] = T1.C11[i1_4];
					temp[i2_4+36] = T1.C20[i1_4];
					temp[i2_4+40] = T1.C21[i1_4];
					temp[i2_4+44] = T1.C22[i1_4];
					temp[i2_4+48] = T1.C30[i1_4];
					temp[i2_4+52] = T1.C31[i1_4];
					temp[i2_4+56] = T1.C32[i1_4];
					temp[i2_4+60] = T1.C33[i1_4];
					temp[i2_4+64] = T1.C40[i1_4];
					temp[i2_4+68] = T1.C41[i1_4];
					temp[i2_4+72] = T1.C42[i1_4];
					temp[i2_4+76] = T1.C43[i1_4];
					temp[i2_4+80] = T1.C44[i1_4];
					temp[i2_4+84] = T1.chi2[i1_4];
					temp[i2_4+88] = T1.NDF[i1_4];

					temp = (fscal *)(&f2.cx0);
					temp[i2_4] = f1.cx0[i1_4];
					temp[i2_4+4] = f1.cx1[i1_4];
					temp[i2_4+8] = f1.cx2[i1_4];
					temp[i2_4+12] = f1.cy0[i1_4];
					temp[i2_4+16] = f1.cy1[i1_4];
					temp[i2_4+20] = f1.cy2[i1_4];
					temp[i2_4+24] = f1.cz0[i1_4];
					temp[i2_4+28] = f1.cz1[i1_4];
					temp[i2_4+32] = f1.cz2[i1_4];
					temp[i2_4+36] = f1.z0[i1_4];
				}
				temp = (fscal *)(&u_front_2[i2_V]);
				temp[i2_4] = xm;
				temp = (fscal *)(&u_back_2[i2_V]);
				temp[i2_4] = vm;
#else
				{														//	  1.83/100 sec
					T2.x[i2_4] = T1.x[i1_4];
					T2.y[i2_4] = T1.y[i1_4];
					T2.tx[i2_4] = T1.tx[i1_4];
					T2.ty[i2_4] = T1.ty[i1_4];
					T2.qp[i2_4] = T1.qp[i1_4];
					T2.z[i2_4] = T1.z[i1_4];
					T2.C00[i2_4] = T1.C00[i1_4];
					T2.C10[i2_4] = T1.C10[i1_4];
					T2.C11[i2_4] = T1.C11[i1_4];
					T2.C20[i2_4] = T1.C20[i1_4];
					T2.C21[i2_4] = T1.C21[i1_4];
					T2.C22[i2_4] = T1.C22[i1_4];
					T2.C30[i2_4] = T1.C30[i1_4];
					T2.C31[i2_4] = T1.C31[i1_4];
					T2.C32[i2_4] = T1.C32[i1_4];
					T2.C33[i2_4] = T1.C33[i1_4];
					T2.C40[i2_4] = T1.C40[i1_4];
					T2.C41[i2_4] = T1.C41[i1_4];
					T2.C42[i2_4] = T1.C42[i1_4];
					T2.C43[i2_4] = T1.C43[i1_4];
					T2.C44[i2_4] = T1.C44[i1_4];
					T2.chi2[i2_4] = T1.chi2[i1_4];
					T2.NDF[i2_4] = T1.NDF[i1_4];
	
					f2.cx0[i2_4] = f1.cx0[i1_4];
					f2.cx1[i2_4] = f1.cx1[i1_4];
					f2.cx2[i2_4] = f1.cx2[i1_4];
					f2.cy0[i2_4] = f1.cy0[i1_4];
					f2.cy1[i2_4] = f1.cy1[i1_4];
					f2.cy2[i2_4] = f1.cy2[i1_4];
					f2.cz0[i2_4] = f1.cz0[i1_4];
					f2.cz1[i2_4] = f1.cz1[i1_4];
					f2.cz2[i2_4] = f1.cz2[i1_4];
					f2.z0[i2_4] = f1.z0[i1_4];

				}
	
				u_front_2[i2_V][i2_4] = xm;
				u_back_2 [i2_V][i2_4] = vm;
#endif // HAVE_VC_SSE
				
// 				hitsl_2[i2] = hitsl_1[i1];
        hitsl_2.push_back(hitsl_1[i1]);
			}
			
			for( int i2_4=n2_4-1; n2_4<fvecLen; n2_4++ ){

#ifdef HAVE_VC_SSE
				fscal *temp;
				temp = (fscal *)(&T2.x);
				for( int i=0; i<=88; i+=4 ) temp[i+n2_4]=temp[i+i2_4];
				temp = (fscal *)(&f2.cx0);
				for( int i=0; i<=36; i+=4 ) temp[i+n2_4]=temp[i+i2_4];
				temp = (fscal *)(&u_front_2[n2_V]);
				temp[n2_4] = temp[i2_4];
				temp = (fscal *)(&u_back_2[n2_V]);
				temp[n2_4] = temp[i2_4];
#else	
				for( fvec *i=&T2.x; i<=&T2.NDF; i++ ) (*i)[n2_4]=(*i)[i2_4];
				for( fvec *i=&f2.cx0; i<=&f2.z0; i++ ) (*i)[n2_4]=(*i)[i2_4];
				
				u_front_2[n2_V][n2_4] = u_front_2[n2_V][i2_4];
				u_back_2 [n2_V][n2_4] = u_back_2 [n2_V][i2_4];
#endif // HAVE_VC_SSE				
			}	

			L1Filter( T2, stam.frontInfo, u_front_2[n2_V] );
			L1Filter( T2, stam.backInfo,  u_back_2[n2_V] );
			
			L1AddMaterial( T2, stam.materialInfo, T2.qp );
			
			L1Extrapolate( T2, star.z, T2.qp, f2 );	
			fvec dxm_est = Pick_r*sqrt(fabs(T2.C00+star.XYInfo.C00));
			fvec dym_est = Pick_r*sqrt(fabs(T2.C11+star.XYInfo.C11));	
// 			x_minusV_2[n2_V] = T2.x - dxm_est;
// 			x_plusV_2 [n2_V] = T2.x + dxm_est;	  
// 			y_minusV_2[n2_V] = T2.y - dym_est;
// 			y_plusV_2 [n2_V] = T2.y + dym_est;		
      x_minusV_2.push_back(T2.x - dxm_est);
      x_plusV_2 .push_back(T2.x + dxm_est);
      y_minusV_2.push_back(T2.y - dym_est);
      y_plusV_2 .push_back(T2.y + dym_est);

      T_2.push_back(T2);
      fld_2.push_back(f2);
		}	// if (n2_V*fvecLen < n2)
			
	}	// if

  fscal *x_minus  = (fscal*) (&(x_minusV_2[0]));
  fscal *x_plus  =  (fscal*) (&(x_plusV_2 [0]));
  fscal *y_minus =  (fscal*) (&(y_minusV_2[0]));
  fscal *y_plus  =  (fscal*) (&(y_plusV_2 [0]));

  
	int n3_V = 0, n3_4 = 0;	
	
	for( int i2 = 0; i2 < n2; i2++){  
		int i2_V = i2/fvecLen;
		int i2_4 = i2%fvecLen;
				
		unsigned /*short*/ int duplet_b = mrDuplets_start[hitsm_2[i2]];			// < 0.001/100000 sec (cash!)
		unsigned /*short*/ int duplet_e = mrDuplets_start[hitsm_2[i2]+1];

//     unsigned /*short*/ int nm = vStsHits_m[hitsm_2[i2]].n;
		L1TrackPar &T2 = T_2[i2_V];

    T_3.push_back(L1TrackPar_0);
    u_front_3.push_back(fvec_0);
    u_back_3.push_back(fvec_0);
		for (int irh_index = duplet_b; irh_index < duplet_e; irh_index++){ //  2.1/10 sec

// 			if( n3 >= MaxPortionTriplets ) break;  // dbg 0
			int irh = mrDuplets_hits[irh_index];

			L1HitPoint &hitr = vStsHits_r[irh];

			fscal xr = hitr.x;   // 0.000/1000000  sec
			fscal vr = hitr.v;
			fscal yr = hitr.y;
//       unsigned /*short*/ int nr = hitr.n;

			if (yr < y_minus[i2]) continue;   // < 0.001/1000
		
			if (yr > y_plus [i2] ) break;	      
			if ((xr < x_minus[i2]) || (xr > x_plus[i2]) /*|| ( nr != nm )*/ ) continue; // dbg 1
																						  
			L1TrackPar &T3 = T_3[n3_V];

// 			hitsl_3[n3] = hitsl_2[i2];   // 1/3  2.3/1000 sec
// 			hitsm_3[n3] = hitsm_2[i2];	 // 1/3  2.3/1000 sec
// 			hitsr_3[n3] = irh;					 // 1/3  2.3/1000 sec
      hitsl_3.push_back(hitsl_2[i2]);   // 1/3  2.3/1000 sec
      hitsm_3.push_back(hitsm_2[i2]);   // 1/3  2.3/1000 sec
      hitsr_3.push_back(irh);           // 1/3  2.3/1000 sec

#ifdef HAVE_VC_SSE		
			fscal *temp;
			{			

				temp = (fscal *)(&T3.x);
				temp[n3_4] = T2.x[i2_4];
				temp[n3_4+4] = T2.y[i2_4];
				temp[n3_4+8] = T2.tx[i2_4];
				temp[n3_4+12] = T2.ty[i2_4];
				temp[n3_4+16] = T2.qp[i2_4];
				temp[n3_4+20] = T2.z[i2_4];
				temp[n3_4+24] = T2.C00[i2_4];
				temp[n3_4+28] = T2.C10[i2_4];
				temp[n3_4+32] = T2.C11[i2_4];
				temp[n3_4+36] = T2.C20[i2_4];
				temp[n3_4+40] = T2.C21[i2_4];
				temp[n3_4+44] = T2.C22[i2_4];
				temp[n3_4+48] = T2.C30[i2_4];
				temp[n3_4+52] = T2.C31[i2_4];
				temp[n3_4+56] = T2.C32[i2_4];
				temp[n3_4+60] = T2.C33[i2_4];
				temp[n3_4+64] = T2.C40[i2_4];
				temp[n3_4+68] = T2.C41[i2_4];
				temp[n3_4+72] = T2.C42[i2_4];
				temp[n3_4+76] = T2.C43[i2_4];
				temp[n3_4+80] = T2.C44[i2_4];
				temp[n3_4+84] = T2.chi2[i2_4];
				temp[n3_4+88] = T2.NDF[i2_4];

			}
			temp = (fscal *)(&u_front_3[n3_V]);
			temp[n3_4] = xr;
			temp = (fscal *)(&u_back_3[n3_V]);
			temp[n3_4] = vr;

#else // HAVE_VC_SSE	

			{														//	2.0/100 sec
				T3.x[n3_4] = T2.x[i2_4];
				T3.y[n3_4] = T2.y[i2_4];
				T3.tx[n3_4] = T2.tx[i2_4];
				T3.ty[n3_4] = T2.ty[i2_4];
				T3.qp[n3_4] = T2.qp[i2_4];
				T3.z[n3_4] = T2.z[i2_4];
				T3.C00[n3_4] = T2.C00[i2_4];
				T3.C10[n3_4] = T2.C10[i2_4];
				T3.C11[n3_4] = T2.C11[i2_4];
				T3.C20[n3_4] = T2.C20[i2_4];
				T3.C21[n3_4] = T2.C21[i2_4];
				T3.C22[n3_4] = T2.C22[i2_4];
				T3.C30[n3_4] = T2.C30[i2_4];
				T3.C31[n3_4] = T2.C31[i2_4];
				T3.C32[n3_4] = T2.C32[i2_4];
				T3.C33[n3_4] = T2.C33[i2_4];
				T3.C40[n3_4] = T2.C40[i2_4];
				T3.C41[n3_4] = T2.C41[i2_4];
				T3.C42[n3_4] = T2.C42[i2_4];
				T3.C43[n3_4] = T2.C43[i2_4];
				T3.C44[n3_4] = T2.C44[i2_4];
				
				T3.chi2[n3_4] = T2.chi2[i2_4];
				T3.NDF[n3_4] = T2.NDF[i2_4];
			}

			u_front_3[n3_V][n3_4] = xr;
			u_back_3 [n3_V][n3_4] = vr; 
#endif // HAVE_VC_SSE			
														
			n3++;										// < 0.001/1000000
			n3_V = n3/fvecLen;			// < 0.001/1000000
			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); }
		}
	}// i2
	
	

	if( n3_4 > 0 ){   // < 0.004/1000
		for( int i3_4=n3_4-1; n3_4<fvecLen; n3_4++ ){
			L1TrackPar &T3 = T_3[n3_V];

#ifdef HAVE_VC_SSE
				fscal *temp;
				temp = (fscal *)(&T3.x);
				for( int i=0; i<=88; i+=4 ) temp[i+n3_4]=temp[i+i3_4];
				temp = (fscal *)(&u_front_2[n3_V]);
				temp[n3_4] = temp[i3_4];
				temp = (fscal *)(&u_back_2[n3_V]);
				temp[n3_4] = temp[i3_4];
#else	
			for( fvec *i=&T3.x, *j=&T3.x; i<=&T3.NDF; i++,j++ )  (*i)[n3_4]=(*j)[i3_4];
			
			u_front_3[n3_V][n3_4] = u_front_3[n3_V][i3_4];
			u_back_3 [n3_V][n3_4] = u_back_3 [n3_V][i3_4];
#endif // HAVE_VC_SSE		
		}	
	}
}

					
					/// Add the right hits to parameters estimation.
inline void L1Algo::f31(	// input
								int n3_V,  
								L1Station &star, 
                nsL1::vector_fvec &u_front, nsL1::vector_fvec &u_back,
									// output
                //                L1TrackPar *T_3
                nsL1::vector<L1TrackPar>::TSimd &T_3
							 )
{	
	for( int i3_V=0; i3_V<n3_V; 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
	};
}



					/// Fill vTriplets.
inline void L1Algo::f4(	// input
								int n3, int istal,
// 								L1TrackPar *T_3,
                nsL1::vector<L1TrackPar>::TSimd &T_3,
								double TRACK_CHI2_CUT, fscal MaxInvMom,
								int *StsHitsStartIndex,  
                vector<unsigned /*short*/ int> &hitsl_3,  vector<unsigned /*short*/ int> &hitsm_3,  vector<unsigned /*short*/ int> &hitsr_3,
								// output
								unsigned &nstaltriplets,
								std::vector<L1Triplet> &vTriplets,
								unsigned *TripStartIndexH, unsigned *TripStopIndexH
// #ifdef XXX								
// 								,unsigned int &stat_n_trip			
// #endif
							 )

{	
				int istam = istal + 1;
				int istar = istal + 2;

				for( int i3=0; i3<n3; i3++){
					int i3_V = i3/fvecLen;
					int i3_4 = i3%fvecLen;
					
					L1TrackPar &T3 = T_3[i3_V];
					
					fscal chi2 = T3.chi2[i3_4];
					if ( !finite(chi2) || chi2<0 ||chi2 > /*TRACK_CHI2_CUT*/ 5. * (T3.NDF[i3_4]-5))	continue;  // < 0.001/100000 sec
// for(int ii = 0; ii < 100000; ii++){
// 
// }	
				
					L1Triplet trip;
				
// 					fscal tx = fabs(T3.tx[i3_4]);
// 					if( tx>1 )	continue;										// not effectiv. just kill secondary tracks.
// 					fscal ty = fabs(T3.ty[i3_4]);
// 					if( ty>1 )	continue;
					
					fscal qp = MaxInvMom + T3.qp[i3_4];						// 1/4 3.0/100000 sec
					if( qp < 0 )						qp = 0;								// 1/4 3.0/100000 sec
					if( qp > MaxInvMom*2 )	qp = MaxInvMom*2;     // 1/4 3.0/100000 sec
					fscal Cqp = 5.*sqrt(fabs(T3.C44[i3_4]));      // 1/4 3.0/100000 sec
					//cout<<520<<endl;
					fscal scale = 255/(MaxInvMom*2);
					//cout<<qp*scale<<" "<<Cqp*scale<<endl;
					qp    = ((unsigned int)(qp*scale) )%256;
					Cqp   = ((unsigned int)(Cqp*scale))%256;
					Cqp+= 1;
					//cout<<525<<endl;
					if( Cqp<0   ) Cqp= 0;
					//if( Cqp>255 ) Cqp= 255;
					if( Cqp>20 ) Cqp= 20;
					qp = (unsigned char) qp;
					//cout<<"qp "<<T3.qp[i3_4]<<" Cqp "<<T3.C44[i3_4]
					//<<", qqp="<<qp<<", Cqqp="<<Cqp<<endl;
					//if( Cqp>30 ) cout<<"track: "<<
					//T3.x[i3_4]<<" "<<T3.y[i3_4]<<" "<<T3.z[i3_4]<<" "<<T3.tx[i3_4]<<" "<<T3.ty[i3_4]<<" "<<T3.qp[i3_4]<<endl;
					//cout<<530<<endl;

					unsigned /*short*/ int ihitl = hitsl_3[i3] + StsHitsStartIndex[istal];
					unsigned /*short*/ int ihitm = hitsm_3[i3] + StsHitsStartIndex[istam];
					unsigned /*short*/ int ihitr = hitsr_3[i3] + StsHitsStartIndex[istar];
				
					
					trip.Set( ihitl, ihitm, ihitr, 
							0, 0, 
							0, (unsigned char) qp, chi2);		// 1.2/1000 sec

					//cout << "Triplet Qp " << trip.GetQpOrig(MaxInvMom) << " mom " << 1.0/trip.GetQpOrig(MaxInvMom) << endl;
				
					//if( trip.GetLHit()<0 || trip.GetLHit()>60000) cout<<"trip.GetLHit()="<<trip.GetLHit()<<endl;
					trip.Cqp   = (unsigned char) Cqp;
						
					//count statistics
					
					vTriplets.push_back(trip);  // 0.67/1000 sec
					
					if ( TripStopIndexH[ihitl] == 0 )
						TripStartIndexH[ihitl] = nstaltriplets;
					TripStopIndexH[ihitl] = ++nstaltriplets;	// < 0.001/100000

// #ifdef XXX				
// 					stat_n_trip++;  
// #endif
				}//i3
}				

					/// Find neighbours of triplets.
					
inline void L1Algo::f5(	// input
							 int NStations, int *TripStartIndex, int *TripStopIndex, 
							 	// output
							 std::vector<L1Triplet> &vTriplets,
							 unsigned *TripStartIndexH, unsigned *TripStopIndexH,
							 int *nlevel
							 )
{	
// 		cout << "Start  Finding levels of triplets "  << endl;

	for	(int istal = NStations-4; istal >=0; istal--){
		int istam = istal + 1;
		int istar = istal + 2;	

					
		unsigned offset_m = TripStartIndex[istam];
		unsigned offset_r = TripStartIndex[istar];

		for (int it = TripStartIndex[istal]; it < TripStopIndex[istal]; it++){
		

			L1Triplet *trip = &(vTriplets[it]);

			
			unsigned char level = 0;
			float  chi2 = trip->GetChi2();
			unsigned char  qp = trip->GetQp();
									
			unsigned /*short*/ int ihitl = trip->GetLHit();
			unsigned /*short*/ int ihitm = trip->GetMHit();
			unsigned /*short*/ int ihitr = trip->GetRHit();

			unsigned /*short*/ int first_neighbour = 0;
			unsigned /*short*/ int end_neighbour = 0;

				// have common 2 hits 
			unsigned int first_triplet;
			unsigned int last_triplet;
			unsigned itn = TripStartIndexH[ihitm];
			unsigned en = TripStopIndexH[ihitm]; 	// end posible neighbour
			
				// find first triplets with 2 same hits
			for (; (itn < en); ++itn){
				unsigned jhitm;
				if( vTriplets[itn].GetLevel()==0 ) 
					jhitm = vTriplets[itn].GetMHit();	      
				else
					jhitm = vTriplets[offset_r+vTriplets[itn].GetFirstNeighbour()].GetLHit();	     
				if (jhitm == ihitr) break;
			};
			
			if (itn < en) { // if exist any one trip
				first_triplet = itn;
					// find last triplets with 2 same hits
				for (; itn < en; ++itn){
					if (vTriplets[itn].GetLHit() != ihitm) break;
					unsigned jhitm;
					if( vTriplets[itn].GetLevel()==0 ) 
							jhitm = vTriplets[itn].GetMHit();	      
						else 
							jhitm = vTriplets[offset_r+vTriplets[itn].GetFirstNeighbour()].GetLHit();	     
					if (jhitm != ihitr) break;
				};
				last_triplet = itn-1;

					// have same qp
				for (itn = first_triplet; itn <= last_triplet; ++itn){ 
					unsigned int index = itn - offset_m;
					L1Triplet *tripn = &vTriplets[itn];
		
					fscal qp2 = tripn->GetQp();
					fscal Cqp1 = trip->Cqp;
					fscal Cqp2 = tripn->Cqp;
					if ( fabs(qp - qp2) > ( Cqp1 + Cqp2 ) )	continue;

					
					if( end_neighbour == 0 ){
						if( index > tripn->GetMaxFirstNeighbour() )	break;
						first_neighbour = index;
					}
					if( index + 1 - first_neighbour >= tripn->GetMaxNNeighbours() )	break;
					end_neighbour = index + 1;	  
					unsigned char jlevel = tripn->GetLevel();
					if( level<=jlevel) level = jlevel + 1 ;		

				}
			
			}; // if exist any one trip
			trip->Set( ihitl, ihitm, ihitr, 
					first_neighbour, end_neighbour-first_neighbour, 
					level, qp, chi2);

			nlevel[level]++;

		}// vTriplets
	} // istal
}	

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

inline void L1Algo::DupletsStaPort(	// input
											int isec,
											int istal, L1Station *vStations, int NStations,
											int *StsHitsStartIndex, int *StsHitsStopIndex,
											vector<unsigned char> &vSFlag, vector<unsigned char> &vSFlagB,
											std::vector<L1HitPoint> &vStsHits,
											double Pick_r,
											double Pick_m, double MaxInvMom,
											fvec targX, fvec targY, fvec targZ, L1FieldValue &targB, L1XYMeasurementInfo &TargetXYInfo,
											
//                      int *n_g1, unsigned *portionStopIndex,
//                      L1TrackPar *T_g1,
//                      L1FieldRegion *fld_g1,
//                      unsigned /*short*/ int *hitsl_g1,

                      vector<int> &n_g1, unsigned *portionStopIndex,
                      L1Portion<L1TrackPar> &T_g1,
                      L1Portion<L1FieldRegion> &fld_g1,
                      L1Portion<unsigned /*short*/ int> &hitsl_g1,
						
												// output
											map<unsigned /*short*/ int,unsigned /*short*/ int> *Duplets_start, vector<unsigned /*short*/ int> *Duplets_hits,
											
//                      int *n_g2,
//                      unsigned /*short*/ int *i1_g2,
//                      unsigned /*short*/ int *hitsm_g2

                      vector<int> &n_g2,
                      L1Portion<unsigned /*short*/ int> &i1_g2,
                      L1Portion<unsigned /*short*/ int> &hitsm_g2
					
											)
{
			int istam = istal+1;
			int istar = istal+2;

			L1Station &stal = vStations[istal];
			L1Station &stam = vStations[istam];
	
						// set lm and mr doublets place for first iteration of loop
						
			map<unsigned /*short*/ int, unsigned /*short*/ int>
			&mrDuplets_start = Duplets_start[istal+1],// mrDuplets_start - first right hit in mrDuplets_hits, array paralel to middle hits
			&lmDuplets_start = Duplets_start[istal];// lmDuplets_start -
      vector<unsigned /*short*/ int>
//       &mrDuplets_hits = Duplets_hits[istal+1], // mrDuplets_hits - right hits of middle-right doublets
      &lmDuplets_hits = Duplets_hits[istal]; // lmDuplets_hits - same for left-middle doublets

			L1HitPoint *vStsHits_l = &(vStsHits[0]) + StsHitsStartIndex[istal];
			L1HitPoint *vStsHits_m = &(vStsHits[0]) + StsHitsStartIndex[istam];
			L1HitPoint *vStsHits_r = 0;

			int NHits_m = StsHitsStopIndex[istam] - StsHitsStartIndex[istam] + 1;
			int NHits_r = 0;
			if( istar < NStations ){
				vStsHits_r = &(vStsHits[0]) + StsHitsStartIndex[istar];
				NHits_r = StsHitsStopIndex[istar] - StsHitsStartIndex[istar] + 1;
			}
	
		
//-----
           
			unsigned int nDuplets_lm = 0;  // number of created doublets on this stantion
			int start_mhit = 0;      //hit on the middle stantion to start

			for(unsigned ip = portionStopIndex[istal+1]; ip < portionStopIndex[istal]; ip++ ){


					
				fvec u_front[Portion/fvecLen], u_back[Portion/fvecLen];				
				
					/// prepare the portion of left hits data
		
				
// 				unsigned /*short*/ int *hitsl_1;
// 				hitsl_1 = &(hitsl_g1[ip*Portion+0]);
				hitsl_g1.add_void(); vector<unsigned /*short*/ int> &hitsl_1 = hitsl_g1[ip];

				int n1 = n_g1[ip];

// 				cout << "Run f10. n1 = " << n1 << endl;

				f10(	// input
						(ip-portionStopIndex[istal+1]) * Portion, n1, vStsHits_l, 
							// output
						u_front, u_back, 
						hitsl_1								
					);
				
				int n1_V = (n1+fvecLen-1)/fvecLen;
			
					/// Get the field approximation. Add the target to parameters estimation. Propagaete to middle station.
		
// 				L1TrackPar *T_1;			// parameters estimation

// 				L1FieldRegion *fld_1;// field approximation on track

// 				T_1 = &(T_g1[ip*Portion/fvecLen+0]);
        T_g1.add_void(); nsL1::vector<L1TrackPar>::TSimd &T_1 = T_g1[ip];
// 				fld_1 = &(fld_g1[ip*Portion/fvecLen+0]);
        fld_g1.add_void(); vector<L1FieldRegion> &fld_1 = fld_g1[ip];
				
				fvec x_minusV_1[Portion/fvecLen];		// region on next station to add hit
				fvec x_plusV_1 [Portion/fvecLen];
				fvec y_minusV_1[Portion/fvecLen];
				fvec y_plusV_1 [Portion/fvecLen];

//         cout << "Run f11" << endl;

				f11(isec, istal, vStations,
						n1_V, 
						MaxInvMom, Pick_m,
						targX, targY, targZ, targB, TargetXYInfo,
						u_front, u_back, 								
							// output
						T_1, fld_1,
						x_minusV_1, x_plusV_1, y_minusV_1, y_plusV_1
					);

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

				fscal *x_minus;
				fscal *x_plus;
				fscal *y_minus;
				fscal *y_plus;
				x_minus  = (fscal*) (&(x_minusV_1[0]));
				x_plus   = (fscal*) (&(x_plusV_1 [0]));
				y_minus  = (fscal*) (&(y_minusV_1[0]));
				y_plus   = (fscal*) (&(y_plusV_1 [0]));

// 				unsigned /*short*/ int *hitsm_2;
// 				hitsm_2 = &(hitsm_g2[ip*MaxPortionDoublets+0]);
        hitsm_g2.add_void(); vector<unsigned /*short*/ int> &hitsm_2 = hitsm_g2[ip];

// 				unsigned /*short*/ int *i1_2;
// 				i1_2 = &(i1_g2[ip*MaxPortionDoublets+0]);											
        i1_g2.add_void(); vector<unsigned /*short*/ int> &i1_2 = i1_g2[ip];

 				int n2;


				f20(	// input
						n1, istar, NStations, stal, stam,
						vStsHits_l, vStsHits_m, NHits_m,
						y_minus, x_minus, y_plus, x_plus,								
						T_1, fld_1,								
						hitsl_1,
						mrDuplets_start,
							// output
						n2,
						i1_2,
						start_mhit,
						hitsm_2,
						lmDuplets_start, lmDuplets_hits, nDuplets_lm, 
						vSFlag, vSFlagB
					 );


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

//  				int n2_V = (n2+fvecLen-1)/fvecLen;
				
// 				fvec *x_minusV_2;
// 				fvec *x_plusV_2;
// 				fvec *y_minusV_2;
// 				fvec *y_plusV_2;		
// 				x_minusV_2 = (&(x_minusV_g2[ip*MaxPortionDoublets/fvecLen+0]));
// 				x_plusV_2  = (&(x_plusV_g2 [ip*MaxPortionDoublets/fvecLen+0]));
// 				y_minusV_2 = (&(y_minusV_g2[ip*MaxPortionDoublets/fvecLen+0]));
// 				y_plusV_2  = (&(y_plusV_g2 [ip*MaxPortionDoublets/fvecLen+0]));	
							
// 				f21(	// input
// 						n2_V, 
// 						stam, star,
// 						u_front, u_back,
// 						fld_2,
// 						Pick_r,
// 							// output
// 						T_2,
// 						x_minusV_2, x_plusV_2, y_minusV_2, y_plusV_2
// 					 );
	
		 
// 				if (n2 >= MaxPortionDoublets) cout << "isec: " << isec << " stantion: " << istal << " portion number: " << ip << " CATrackFinder: Warning: Too many doublets created in portion" << endl;	
					 					 
// 				n_g2[ip] = n2;
        n_g2.push_back(n2);

			} // portion of left hits 

						
// 		if (nDuplets_lm > MaxArrSize) cout << "isec: " << isec << " stantion: " << istal << " CATrackFinder: Error: Too many doublets on stantion was created" << endl;


}


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

inline void L1Algo::TripletsStaPort(	// input
														int isec,
														int istal, L1Station *vStations, int NStations,
														int *StsHitsStartIndex, int *StsHitsStopIndex,
														std::vector<L1HitPoint> &vStsHits,
														double Pick_r, double TRACK_CHI2_CUT, double MaxInvMom, 

// 														int *n_g1,
// 														L1TrackPar *T_g1, 
// 														L1FieldRegion *fld_g1,								
// 														unsigned /*short*/ int *hitsl_g1,

                            vector<int> &n_g1,
                            L1Portion<L1TrackPar> &T_g1,
                            L1Portion<L1FieldRegion> &fld_g1,
                            L1Portion<unsigned /*short*/ int> &hitsl_g1,
																	
	                          vector<int> &n_g2, unsigned *portionStopIndex,
                            L1Portion<unsigned /*short*/ int> &i1_g2,
                            L1Portion<unsigned /*short*/ int> &hitsm_g2,
															
															// output
														map<unsigned /*short*/ int,unsigned /*short*/ int> *Duplets_start, vector<unsigned /*short*/ int>  *Duplets_hits,
														std::vector<L1Triplet> *vTriplets_part,
														unsigned *TripStartIndexH, unsigned *TripStopIndexH
														)
{

      map<unsigned /*short*/ int, unsigned /*short*/ int>
      &mrDuplets_start = Duplets_start[istal+1],// mrDuplets_start - first right hit in mrDuplets_hits, array paralel to middle hits
      &lmDuplets_start = Duplets_start[istal];// lmDuplets_start -
      vector<unsigned /*short*/ int>
      &mrDuplets_hits = Duplets_hits[istal+1], // mrDuplets_hits - right hits of middle-right doublets
      &lmDuplets_hits = Duplets_hits[istal]; // lmDuplets_hits - same for left-middle doublets

			int istam = istal+1;
			int istar = istal+2;
			L1Station &stam = vStations[istam];
			L1Station &star = vStations[istar];
			
			L1HitPoint *vStsHits_m = &(vStsHits[0]) + StsHitsStartIndex[istam];
 			
			L1HitPoint *vStsHits_r = 0;
			int NHits_r = 0;
			if( istar < NStations ){
				vStsHits_r = &(vStsHits[0]) + StsHitsStartIndex[istar];
				NHits_r = StsHitsStopIndex[istar] - StsHitsStartIndex[istar] + 1;
			}

			unsigned nstaltriplets = 0; // find triplets begin from this stantion

			for(unsigned ip = portionStopIndex[istal+1]; ip < portionStopIndex[istal]; ip++ ){		

	
				int n3=0, n3_V;
				
// 				L1TrackPar *T_1;								
// 				T_1 = &(T_g1[ip*Portion/fvecLen+0]);
        nsL1::vector<L1TrackPar>::TSimd &T_1 = T_g1[ip];
// 				L1FieldRegion *fld_1;
//         fld_1 = &(fld_g1[ip*Portion/fvecLen+0]);        
        vector<L1FieldRegion> &fld_1 = fld_g1[ip];
// 				unsigned /*short*/ int *hitsl_1;
// 				hitsl_1 = &(hitsl_g1[ip*Portion+0]);
        vector<unsigned /*short*/ int> &hitsl_1 = hitsl_g1[ip];

				int n1 = n_g1[ip];

// 				unsigned /*short*/ int *hitsm_2;
// 				hitsm_2 = &(hitsm_g2[ip*MaxPortionDoublets+0]);
        vector<unsigned /*short*/ int> &hitsm_2 = hitsm_g2[ip];

// 				unsigned /*short*/ int *i1_2;
// 				i1_2 = &(i1_g2[ip*MaxPortionDoublets+0]);
        vector<unsigned /*short*/ int> &i1_2 = i1_g2[ip];
				int n2 = n_g2[ip];

        
// 				L1TrackPar T_3[MaxPortionTriplets/fvecLen];
        nsL1::vector<L1TrackPar>::TSimd T_3;

// 				unsigned /*short*/ int hitsl_3[MaxPortionTriplets];
// 				unsigned /*short*/ int hitsm_3[MaxPortionTriplets];
// 				unsigned /*short*/ int hitsr_3[MaxPortionTriplets];
        vector<unsigned /*short*/ int> hitsl_3, hitsm_3, hitsr_3;

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

			
// 				fvec u_front3[MaxPortionTriplets/fvecLen], u_back3[MaxPortionTriplets/fvecLen];
        nsL1::vector_fvec u_front3, u_back3;
        u_front3.push_back(fvec(0));
        T_3.reserve(MaxPortionTriplets/fvecLen);
        hitsl_3.reserve(MaxPortionTriplets);
        hitsm_3.reserve(MaxPortionTriplets);
        hitsr_3.reserve(MaxPortionTriplets);
        u_front3.reserve(MaxPortionTriplets/fvecLen);
        u_back3.reserve(MaxPortionTriplets/fvecLen);

					/// Find the triplets(right hit). Reformat data in the portion of triplets.
					
//         cout << mrDuplets_hits.size() << endl; // dbg 12
				f30(	// input
						vStsHits_r, stam, star,

						istar, NStations, n1, 
						vStsHits_m,
						T_1, fld_1,								
						hitsl_1,
						lmDuplets_start, lmDuplets_hits,	
						Pick_r,						
						
						n2, 
						hitsm_2,
						i1_2,
												
						mrDuplets_start, mrDuplets_hits,
// 						MaxPortionTriplets,
							// output
						n3,
						T_3,
						hitsl_3, hitsm_3, hitsr_3,
						u_front3, u_back3
					 );
 				n3_V = (n3+fvecLen-1)/fvecLen;
		
// 				if (n3 >= MaxPortionTriplets) cout << "isec: " << isec << " stantion: " << istal << " portion number: " << ip << " CATrackFinder: Warning: Too many Triplets created in portion" << endl;	
					 					 

				
			
					/// Add the right hits to parameters estimation.

				f31(	// input
						n3_V,  
						star, 
						u_front3, u_back3,
							// output
						T_3
					 );
					 

					/// Fill Triplets.

				f4(	// input
					n3, istal,
					T_3, 
					TRACK_CHI2_CUT, MaxInvMom,
					StsHitsStartIndex,  
					hitsl_3, hitsm_3, hitsr_3,								
						// output
					nstaltriplets,
					vTriplets_part[istal],
					TripStartIndexH, TripStopIndexH
					);
			
			} // start_lh - portion of left hits 

}


  /// 					--------------------------------------------	ITBB part ----------------------------
	
#ifdef TBB 
	
#include "tbb/task_scheduler_init.h"
#include "tbb/blocked_range.h"
#include "tbb/parallel_for.h"
using namespace tbb;

class ParalleledDup {
		// input
	int isec;
	L1Station *vStations; int NStations;
	int *StsHitsStartIndex; int *StsHitsStopIndex;
	vector<unsigned char> &vSFlag; vector<unsigned char> &vSFlagB;
	std::vector<L1HitPoint> &vStsHits;
	double Pick_r;
	double Pick_m; double MaxInvMom;
	fvec targX; fvec targY; fvec targZ; L1FieldValue &targB; L1XYMeasurementInfo &TargetXYInfo;
		
	int *n_g1; unsigned *portionStopIndex;
	L1TrackPar *T_g1;
	L1FieldRegion *fld_g1;
	unsigned /*short*/ int *hitsl_g1;
// 	fvec *x_minusV_g1; fvec *x_plusV_g1; fvec *y_minusV_g1; fvec *y_plusV_g1;											
		// output
	unsigned /*short*/ int *Duplets_start; unsigned /*short*/ int *Duplets_hits;
	
	int *n_g2;
	unsigned /*short*/ int *i1_g2;
// 	L1TrackPar *T_g2;
// 	L1FieldRegion *fld_g2; 
// 	fvec *u_front_g2; fvec *u_back_g2;			
// 	unsigned /*short*/ int *hitsl_g2;
	unsigned /*short*/ int *hitsm_g2;
// 	fvec *x_minusV_g2; fvec *x_plusV_g2; fvec *y_minusV_g2; fvec *y_plusV_g2;		
  
 public:
  void operator()( const blocked_range<size_t>& r ) const;
	
  ParalleledDup(   
									// input
								int isec_,
								L1Station *vStations_, int NStations_,
								int *StsHitsStartIndex_, int *StsHitsStopIndex_,
								vector<unsigned char> &vSFlag_, vector<unsigned char> &vSFlagB_,
								std::vector<L1HitPoint> &vStsHits_,
								double Pick_r_,
								double Pick_m_, double MaxInvMom_,
								fvec targX_, fvec targY_, fvec targZ_, L1FieldValue &targB_, L1XYMeasurementInfo &TargetXYInfo_,			
																
								int *n_g1_, unsigned *portionStopIndex_,
								L1TrackPar *T_g1_,
								L1FieldRegion *fld_g1_,
								unsigned /*short*/ int *hitsl_g1_,
// 								fvec *x_minusV_g1_, fvec *x_plusV_g1_, fvec *y_minusV_g1_, fvec *y_plusV_g1_,											
									// output
								unsigned /*short*/ int *Duplets_start_, unsigned /*short*/ int *Duplets_hits_,
								
								int *n_g2_,
								unsigned /*short*/ int *i1_g2_,
// 								L1TrackPar *T_g2_,
// 								L1FieldRegion *fld_g2_, 
// 								fvec *u_front_g2_, fvec *u_back_g2_,										
// 								unsigned /*short*/ int *hitsl_g2_,
								unsigned /*short*/ int *hitsm_g2_
// 								fvec *x_minusV_g2_, fvec *x_plusV_g2_, fvec *y_minusV_g2_, fvec *y_plusV_g2_		
					 ): 
											isec(isec_),
											vStations(vStations_), NStations(NStations_),
											StsHitsStartIndex(StsHitsStartIndex_), StsHitsStopIndex(StsHitsStopIndex_),
											vSFlag(vSFlag_), vSFlagB(vSFlagB_),
											vStsHits(vStsHits_),
											Pick_r(Pick_r_),
											Pick_m(Pick_m_), MaxInvMom(MaxInvMom_),
											targX(targX_), targY(targY_), targZ(targZ_), targB(targB_), TargetXYInfo(TargetXYInfo_),
																						
											n_g1(n_g1_), portionStopIndex(portionStopIndex_),
											T_g1(T_g1_),
											fld_g1(fld_g1_),
											hitsl_g1(hitsl_g1_),
// 											x_minusV_g1(x_minusV_g1_), x_plusV_g1(x_plusV_g1_), y_minusV_g1(y_minusV_g1_), y_plusV_g1(y_plusV_g1_),											
												// output
											Duplets_start(Duplets_start_), Duplets_hits(Duplets_hits_),
											
											n_g2(n_g2_),
											i1_g2(i1_g2_),
// 											T_g2(T_g2_),
// 											fld_g2(fld_g2_), 
// 											u_front_g2(u_front_g2_), u_back_g2(u_back_g2_),													
// 											hitsl_g2(hitsl_g2_),	
											hitsm_g2(hitsm_g2_)
// 											x_minusV_g2(x_minusV_g2_), x_plusV_g2(x_plusV_g2_), y_minusV_g2(y_minusV_g2_), y_plusV_g2(y_plusV_g2_)	
													{};
		
  ~ParalleledDup(){};
};

void ParalleledDup::operator()( const blocked_range<size_t>& r ) const {

    for(int istal = r.begin(); istal < (int)r.end(); ++istal ){
// 		cout<< " ParalleledDup::operator(). Station: " << istal << endl;		
			DupletsStaPort(	// input
											isec,
											istal, vStations, NStations,
											StsHitsStartIndex, StsHitsStopIndex,
											vSFlag, vSFlagB,
											vStsHits,
											Pick_r,
											Pick_m, MaxInvMom,
											targX, targY, targZ, targB, TargetXYInfo,
																					
											n_g1, portionStopIndex,
											T_g1,
											fld_g1,
											hitsl_g1,
// 											x_minusV_g1, x_plusV_g1, y_minusV_g1, y_plusV_g1,											
												// output
											Duplets_start, Duplets_hits,
											
											n_g2,
											i1_g2,
// 											T_g2,
// 											fld_g2, 
// 											u_front_g2, u_back_g2,														
// 											hitsl_g2,	
											hitsm_g2
// 											x_minusV_g2, x_plusV_g2, y_minusV_g2, y_plusV_g2											
											);

    }
}


class ParalleledTrip {
// input
	int isec;
	L1Station *vStations; int NStations;
	int *StsHitsStartIndex; int *StsHitsStopIndex;
	std::vector<L1HitPoint> &vStsHits;
	double Pick_r; double TRACK_CHI2_CUT; double MaxInvMom; 

														int *n_g1;
														L1TrackPar *T_g1; 
														L1FieldRegion *fld_g1;					
														unsigned /*short*/ int *hitsl_g1;
	
	int *n_g2; unsigned *portionStopIndex;
	unsigned /*short*/ int *i1_g2;
/*	L1FieldRegion *fld_g2; 
	fvec *u_front_g2; fvec *u_back_g2;	
	L1TrackPar *T_g2;*/
// 	unsigned /*short*/ int *hitsl_g2;
	unsigned /*short*/ int *hitsm_g2;
// 	fvec *x_minusV_g2; fvec *x_plusV_g2; fvec *y_minusV_g2; fvec *y_plusV_g2;										
		// output
	unsigned /*short*/ int *Duplets_start; unsigned /*short*/ int *Duplets_hits;
	std::vector<L1Triplet> *vTriplets_part;
	unsigned *TripStartIndexH; unsigned *TripStopIndexH;
  
 public:
  void operator()( const blocked_range<size_t>& r ) const;
	
  ParalleledTrip(   
// input
								int isec_,
								L1Station *vStations_, int NStations_,
								int *StsHitsStartIndex_, int *StsHitsStopIndex_,
								std::vector<L1HitPoint> &vStsHits_,
								double Pick_r_, double TRACK_CHI2_CUT_, double MaxInvMom_, 

														int *n_g1_,
														L1TrackPar *T_g1_, 
														L1FieldRegion *fld_g1_,								
														unsigned /*short*/ int *hitsl_g1_,
													
								int *n_g2_, unsigned *portionStopIndex_,
								unsigned /*short*/ int *i1_g2_,
								/*L1FieldRegion *fld_g2_, 
								fvec *u_front_g2_, fvec *u_back_g2_,
								L1TrackPar *T_g2_,*/
// 								unsigned /*short*/ int *hitsl_g2_,
								unsigned /*short*/ int *hitsm_g2_,
// 								fvec *x_minusV_g2_, fvec *x_plusV_g2_, fvec *y_minusV_g2_, fvec *y_plusV_g2_,										
									// output
								unsigned /*short*/ int *Duplets_start_, unsigned /*short*/ int *Duplets_hits_,
								std::vector<L1Triplet> *vTriplets_part_,
								unsigned *TripStartIndexH_, unsigned *TripStopIndexH_
					 ): 
												// input
											isec(isec_),
											vStations(vStations_), NStations(NStations_),
											StsHitsStartIndex(StsHitsStartIndex_), StsHitsStopIndex(StsHitsStopIndex_),
											vStsHits(vStsHits_),
											Pick_r(Pick_r_), TRACK_CHI2_CUT(TRACK_CHI2_CUT_), MaxInvMom(MaxInvMom_), 
											
														n_g1(n_g1_),
														T_g1(T_g1_), 
														fld_g1(fld_g1_),								
														hitsl_g1(hitsl_g1_),
														
											n_g2(n_g2_), portionStopIndex(portionStopIndex_),
											i1_g2(i1_g2_),
/*											fld_g2(fld_g2_), 
											u_front_g2(u_front_g2_), u_back_g2(u_back_g2_),											
											T_g2(T_g2_),*/
// 											hitsl_g2(hitsl_g2_),	
											hitsm_g2(hitsm_g2_),
// 											x_minusV_g2(x_minusV_g2_), x_plusV_g2(x_plusV_g2_), y_minusV_g2(y_minusV_g2_), y_plusV_g2(y_plusV_g2_),										
												// output
											Duplets_start(Duplets_start_), Duplets_hits(Duplets_hits_),
											vTriplets_part(vTriplets_part_),
											TripStartIndexH(TripStartIndexH_), TripStopIndexH(TripStopIndexH_)
					{};
		
  ~ParalleledTrip(){};
};

void ParalleledTrip::operator()( const blocked_range<size_t>& r ) const {

    for(int istal = r.begin(); istal < (int)r.end(); ++istal ){
// 		cout<< " ParalleledTrip::operator(). Station: " << istal << endl;		
			TripletsStaPort(	// input
											isec,
											istal, vStations, NStations,
											StsHitsStartIndex, StsHitsStopIndex,
											vStsHits,
											Pick_r, TRACK_CHI2_CUT, MaxInvMom, 

														n_g1,
														T_g1, 
														fld_g1,								
														hitsl_g1,											
											
											n_g2, portionStopIndex,
											i1_g2,
/*											fld_g2, 
											u_front_g2, u_back_g2,											
											T_g2,*/
// 											hitsl_g2,	
											hitsm_g2,
// 											x_minusV_g2, x_plusV_g2, y_minusV_g2, y_plusV_g2,										
												// output
											Duplets_start, Duplets_hits,
											vTriplets_part,
											TripStartIndexH, TripStopIndexH
											);

    }
}
#endif // TBB 
  
	/// 										----------------------	end of ITBB part ----------------------------

	
	
	
    ///**********************************************************************************************************************
    ///**********************************************************************************************************************
		///**********************************************************************************************************************
    ///*                                                                                                                    *
		///*                                                                                                                    *
		///*                                                                                                                    *
    ///*                                          CATrackFinder procedure                                                   *
		///*                                                                                                                    *
		///*                                                                                                                    *
    ///*                                                                                                                    *
    ///**********************************************************************************************************************		
    ///**********************************************************************************************************************
		///**********************************************************************************************************************			
	
	


void L1Algo::CATrackFinder()
{

//   cout << "Start TrackFinder" << endl;



// 	cout<<" total STS hits "<<vStsHits.size()<<endl;
		//cout<<" total STS hits, strips, back strips "
		//<<vStsHits.size()<<" "<<vStsStrips.size()<<" "<<vStsStripsB.size()<<endl;
		
#ifdef DRAW
  InitL1Draw();
#endif

#ifdef DEBUG_CA
  InitL1Histo();
#endif
	
	cout.precision(6);
					
  static unsigned int stat_N = 0; // number of events
		
  const int Nisec = 3; //number of main loops		
#ifdef XXX  
		// timers for different part of CATrackFinder
	const int ntimers = 6; //number of timers + 1
	
  static double 
    stat_time = 0.0, stat_time_real = 0.0, stat_time_trip[Nisec], stat_time_fit[Nisec],  stat_time_fit_bra[Nisec], stat_time_fit_fin[Nisec];
		
	if (stat_N == 0)		
		for(int i1 = 0; i1 < Nisec; i1++){
			stat_time_trip[i1] = .0;
			stat_time_fit[i1] = .0;
			stat_time_fit_bra[i1] = .0;
			stat_time_fit_fin[i1] = .0;
		}	
	
  double 
  time_trip[Nisec], time_fit[Nisec], time_fit_bra[Nisec], time_fit_fin[Nisec];

	for(int i1 = 0; i1 < Nisec; i1++){
		time_trip[i1] = .0;
		time_fit[i1] = .0;
    time_fit_bra[i1] = .0;
		time_fit_fin[i1] = .0;
	}		
			
	static double stat_time_find[ntimers][2][Nisec];
	if (stat_N == 0)
		for(int i1 = 0; i1 < Nisec; i1++)
			for(int i2 = 0; i2 < ntimers; i2++)
				for(int i3 = 0; i3 < 2; i3++)
					stat_time_find[i2][i3][i1] = 0.;
	
	double stat_time_find_all[ntimers][2]; 		
	for(int i2 = 0; i2 < ntimers; i2++)
		for(int i3 = 0; i3 < 2; i3++)
			stat_time_find_all[i2][i3] = 0.;		

	
	double time_find[ntimers][2][Nisec];
	for(int i1 = 0; i1 < Nisec; i1++)
		for(int i2 = 0; i2 < ntimers; i2++)
			for(int i3 = 0; i3 < 2; i3++)
				time_find[i2][i3][i1] = 0.;	
		
	TStopwatch c_time_find[ntimers][2];			


#endif


  TStopwatch c_time;	// all CATrackFinder work time
	c_time.Start();		

#ifdef XXX 			
		c_time_find[0][0].Start();
#endif
  
	vTracks.clear();
	vRecoHits.clear();
  
		// set all hits as usable
	for( vector<unsigned char>::iterator is= vSFlag.begin(); is!=vSFlag.end(); ++is)
		SetFUnUsed(*is);
	for( vector<unsigned char>::iterator is= vSFlagB.begin(); is!=vSFlagB.end(); ++is) 
		SetFUnUsed(*is);
		

#ifdef XXX 
	unsigned int 
		stat_max_n_trip = 0,
		stat_max_trip_size = 0,
		stat_max_n_dup = 0,
		stat_max_BranchSize = 0,
		stat_max_n_branches = 0;
#endif
  
#ifdef DRAW
	Draw1();
#endif

	vector< L1StsHit > vStsDontUsedHits_A;	// arrays for copy from vStsHits only don't uses in created tracks hits 
	vector< L1StsHit > vStsDontUsedHits_B;

	vector< L1HitPoint > vStsDontUsedHitsxy_A;	// arrays for contain prepared information for don't uses in created tracks hits 
	vector< L1HitPoint > vStsDontUsedHitsxy_B;	
		
	std::vector< L1StsHit > *svStsHits = &vStsDontUsedHits_B;			// array of hits uses on current iteration of isec-loop
	std::vector< L1StsHit > *svStsHits_buf = &vStsDontUsedHits_A; // buffer for copy
	
	std::vector< L1HitPoint > *svStsHits2 = &vStsDontUsedHitsxy_B;			// array of info for hits uses on current iteration of isec-loop
	std::vector< L1HitPoint > *svStsHits_buf2 = &vStsDontUsedHitsxy_A;
		
	svStsHits_buf->clear();
	svStsHits_buf2->clear();
// 	unsigned int RealIHit[MaxNPortion*Portion*MaxNStations]; // index of hit in vStsHits indexed by index of hit in svStsHits;
  vector<unsigned int> RealIHit_v;
  RealIHit_v.reserve(vStsHits.size());
  unsigned int *RealIHit = &(RealIHit_v[0]);

	
		// full arrays for first iteration
	for(int ista = 0; ista < NStations; ista++){
		for(int ih = StsHitsStartIndex[ista]; ih <= StsHitsStopIndex[ista]; ih++){
		
			RealIHit[ih] = ih;
			vStsDontUsedHits_B.push_back(vStsHits[ih]);
			
			L1Station &sta = vStations[ista];
			L1StsHit &hit = vStsHits[ih];
			L1Strip &x = vStsStrips[hit.f];
			L1Strip &v = vStsStripsB[hit.b];
			fscal y  = (sta.yInfo.cos_phi*x + sta.yInfo.sin_phi*v)[0];
			
      L1HitPoint hitxy(x,y,v,hit.n);
			vStsDontUsedHitsxy_B.push_back(hitxy);
		}
	}
	
#ifdef XXX 			
				c_time_find[0][0].Stop();    
				time_find[0][0][0] += double(c_time_find[0][0].CpuTime()); //ms
				time_find[0][1][0] += double(c_time_find[0][0].RealTime()); //ms
#endif		
		
		// itetation of finding:
		// isec == 0 - primary fast track
    // isec == 1 - primary all track  
    // isec == 2 - secondary all track 
	for (int isec = 0; isec <= fNFindIterations-1; isec++){ // all finder

// 		for( vector<unsigned char>::iterator is= vSFlag.begin(); is!=vSFlag.end(); ++is)
// 			SetFUnUsedD(*is);
// 		for( vector<unsigned char>::iterator is= vSFlagB.begin(); is!=vSFlagB.end(); ++is)
// 			SetFUnUsedD(*is);

#ifdef XXX
//     cout << " Begin of iteration " << isec << endl;
		TStopwatch c_time_trip;  // time for find all triplets
#endif   

		TRACK_CHI2_CUT = 10.0;
 
		Pick_m = 3.0;	// coefficient for size of region on middle station for add middle hits in triplets: Dx = Pick*sigma_x Dy = Pick*sigma_y
		Pick_r = 5.0;	// coefficient for size of region on right  station for add right  hits in triplets
		//double Pick_gather = 5.0; //!!!
		if (isec == 0){
			Pick_m =  2.0;
		}



		//double MaxInvMom = 1.0/1.;		
		//if (isec == 1) MaxInvMom =  1.0/0.2;
		//if (isec == 2) MaxInvMom =  1.0/0.2;		
		MaxInvMom = 1.0/0.5;		 								// max considered q/p  		
		if (isec == 1) MaxInvMom =  1.0/0.1;
		if (isec == 2) MaxInvMom =  1.0/0.1;
		if (isec == fTrackingLevel )  MaxInvMom = 1./fMomentumCutOff;
 

			// define the target
		static L1FieldValue targB _fvecalignment;
		fvec targX= 0, targY= 0, targZ= 0;		  //  suppose, what target will be at (0,0,0)

		if (isec <= 1){ // target
			targB = vtxFieldValue;
			if (isec ==-1)
				SigmaTargetX = SigmaTargetY = 0.01; // target
			else
				SigmaTargetX = SigmaTargetY = 0.1;
		}
		else{ //use outer radius of the 1st station as a constraint
			L1Station &st = vStations[0];
			SigmaTargetX = SigmaTargetY = 3;//st.Rmax[0];
			targZ = 0.;//-1.;
			st.fieldSlice.GetFieldValue( 0, 0, targB );
		}

		static L1XYMeasurementInfo TargetXYInfo _fvecalignment;
		TargetXYInfo.C00 = SigmaTargetX * SigmaTargetX;
		TargetXYInfo.C10 = 0;
		TargetXYInfo.C11 = SigmaTargetY * SigmaTargetY;


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

// 		static unsigned /*short*/ int Duplets_hits[MaxArrSize*MaxNStations];	// right hits of doublets(left-right)
    vector<unsigned /*short*/ int> Duplets_hits[MaxNStations];
// 		static unsigned /*short*/ int Duplets_start[MaxArrSize*MaxNStations];	// index of first right hits of doublets in Duplets_hits array   indexed by left hit
    map<unsigned /*short*/ int, unsigned /*short*/ int> Duplets_start[MaxNStations];
// 					// set 0 for all doublets begin from last station - they don't exist.

//     Duplets_hits[NStations-1].push_back(0);
// 		for (unsigned i = MaxArrSize*(NStations-1); i < MaxArrSize*NStations; i++){
// 			Duplets_hits[i]  = 0;
// 			Duplets_start[i] = 0;
// 		};


// 		static int n_g1[MaxNPortion];	// number of singlets in same portion
    vector<int> n_g1;
// 		static L1TrackPar T_g1[MaxNPortion*Portion/fvecLen];	// estimation of parameters
    L1Portion<L1TrackPar> T_g1(MaxNPortion, Portion/fvecLen);

// 		static L1FieldRegion fld_g1[MaxNPortion*Portion/fvecLen];	// magnetic field
    L1Portion<L1FieldRegion> fld_g1(MaxNPortion, Portion/fvecLen);
// 		static unsigned /*short*/ int hitsl_g1[MaxNPortion*Portion];  // left and middle hits indexed by number of doublets in portion(i2)
    L1Portion<unsigned /*short*/ int> hitsl_g1(MaxNPortion, Portion);
	
// 		static int n_g2[MaxNPortion];	// number of doublets in same portion
    vector<int> n_g2;

// 		static unsigned /*short*/ int i1_g2[MaxNPortion*MaxPortionDoublets]; 	// index in portion of singlets(i1) indexed by index in portion of doublets(i2)
    L1Portion<unsigned /*short*/ int> i1_g2(MaxNPortion, MaxPortionDoublets);

// 		static unsigned /*short*/ int hitsm_g2[MaxNPortion*MaxPortionDoublets];// left and middle hits indexed by number of doublets in portion(i2)
    L1Portion<unsigned /*short*/ int> hitsm_g2(MaxNPortion, MaxPortionDoublets);

    for (int i = 0; i < MaxNStations; i++){
      Duplets_hits[i].reserve(MaxArrSize);
//       Duplets_start[i].reserve(MaxArrSize);
    }
    n_g1.reserve(MaxNPortion);
//     T_g1.reserve();
//     fld_g1.reserve();
//     hitsl_g1.reserve(MaxNPortion*Portion);
    n_g2.reserve(MaxNPortion);
//     i1_g2.reserve(MaxNPortion*MaxPortionDoublets);
//     hitsm_g2.reserve(MaxNPortion*MaxPortionDoublets);
    

    

    
		/*static*/ vector <L1Triplet> vTriplets_part[MaxNStations];  // temporary arrays for parallelizing save triplets
																
		unsigned int portionStopIndex[MaxNStations];	//number of portion on each station

// 		static unsigned int TripStartIndexH[MaxArrSize], TripStopIndexH[MaxArrSize]; // index region for triplets, which begins from hit, indexed by index of lhit.
    vector<unsigned int> TripStartIndexH_v, TripStopIndexH_v;
    TripStartIndexH_v.reserve(svStsHits->size());
    TripStopIndexH_v.reserve(svStsHits->size());
    unsigned int *TripStartIndexH = &(TripStartIndexH_v[0]), *TripStopIndexH = &(TripStopIndexH_v[0]);
		for(unsigned i = 0; i < svStsHits->size(); i++)
		{
// 			TripStartIndexH[i] = 0;
			TripStopIndexH[i] = 0;
		}

		{		// split left hits on portions
			portionStopIndex[NStations-1] = 0;
			unsigned int ip = 0;	//index of curent portion	
			
			for (int istal = NStations-2; istal >= FIRSTCASTATION; istal--){  //start downstream chambers
				int NHits_l = StsHitsStopIndex[istal] - StsHitsStartIndex[istal] + 1;
				
				int NHits_l_P = NHits_l/Portion;
				
				for( int ipp = 0; ipp < NHits_l_P; ipp++ ){
// 					n_g1[ip++] = Portion;
          n_g1.push_back(Portion);
          ip++;
				} // start_lh - portion of left hits
				
// 				n_g1[ip++] = NHits_l - NHits_l_P*Portion;
        n_g1.push_back(NHits_l - NHits_l_P*Portion);
				ip++;
				portionStopIndex[istal]=ip;
			}// lstations		
		}

    	///   CREATE doublets       

#ifdef XXX
    cout << "CREATE doublets" << endl;
#endif


#ifdef XXX 			
		c_time_find[1][0].Start();
#endif

#ifdef TBB
		{
			task_scheduler_init init(nthreads);
			static affinity_partitioner ap;
			ParalleledDup af(	// input
											isec,
											vStations, NStations,
											StsHitsStartIndex, StsHitsStopIndex,
											vSFlag, vSFlagB,
											*svStsHits2,
											Pick_r,
											Pick_m, MaxInvMom,
											targX, targY, targZ, targB, TargetXYInfo,
																						
											n_g1, portionStopIndex,
											T_g1,
											fld_g1,
											hitsl_g1,
	
												// output
											Duplets_start, Duplets_hits,
											
											n_g2,
											i1_g2,
											hitsm_g2
											);
			parallel_for(blocked_range<size_t>(FIRSTCASTATION,NStations-2+1,
			                                   nblocks
																				), af,ap);
		}
#else // TBB 
		for (int istal = NStations-2; istal >= FIRSTCASTATION; istal--){//  //start downstream chambers
			DupletsStaPort(	// input
											isec,
											istal, vStations, NStations,
											StsHitsStartIndex, StsHitsStopIndex,
											vSFlag, vSFlagB,
											*svStsHits2,
											Pick_r,
											Pick_m, MaxInvMom,
											targX, targY, targZ, targB, TargetXYInfo,
																						
											n_g1, portionStopIndex,
											T_g1,
											fld_g1,
											hitsl_g1,
									
												// output
											Duplets_start, Duplets_hits,
											
											n_g2,
											i1_g2,
											hitsm_g2
											);
    }// lstations
#endif // TBB
	
#ifdef XXX 			
				c_time_find[1][0].Stop();    
				time_find[1][0][isec] += double(c_time_find[1][0].CpuTime()); //ms
				time_find[1][1][isec] += double(c_time_find[1][0].RealTime()); //ms
#endif		

#ifdef XXX
		int istat_n_triplet_c=0; 
		for(unsigned ip = 0; ip < portionStopIndex[FIRSTCASTATION]; ip++)
			istat_n_triplet_c += n_g2[ip];
#endif

    	///   CREATE triplets      

#ifdef XXX
    cout << " CREATE TRIPLETS " << endl;
#endif



#ifdef XXX 			
		c_time_find[1][0].Start();
#endif

#ifdef TBB 
		{
			task_scheduler_init init(nthreads);
			static affinity_partitioner ap;
			ParalleledTrip af(	// input
											isec,
											vStations, NStations,
											StsHitsStartIndex, StsHitsStopIndex,
											*svStsHits2,
											Pick_r, TRACK_CHI2_CUT, MaxInvMom, 
											n_g1,
											T_g1, 
											fld_g1,								
											hitsl_g1,				
											n_g2, portionStopIndex,
											i1_g2,
											hitsm_g2,
												// output
											Duplets_start, Duplets_hits,
											vTriplets_part,
											TripStartIndexH, TripStopIndexH
											);
			parallel_for(blocked_range<size_t>(FIRSTCASTATION,NStations-2+1,
																				nblocks
																				), af,ap);
		}
#else // TBB 		
		for (int istal = NStations-2; istal >= FIRSTCASTATION; istal--){//  //start downstream chambers
			TripletsStaPort(	// input
											isec,
											istal, vStations, NStations,
											StsHitsStartIndex, StsHitsStopIndex,
											*svStsHits2,
											Pick_r, TRACK_CHI2_CUT, MaxInvMom, 

											n_g1,
											T_g1, 
											fld_g1,								
											hitsl_g1,											
											
											n_g2, portionStopIndex,
											i1_g2,
											hitsm_g2,
						
												// output
											Duplets_start, Duplets_hits,
											vTriplets_part,
											TripStartIndexH, TripStopIndexH
											);
    }// l-stations	
#endif // TBB
#ifdef XXX 			
				c_time_find[1][0].Stop();    
				time_find[2][0][isec] += double(c_time_find[1][0].CpuTime()); //ms
				time_find[2][1][isec] += double(c_time_find[1][0].RealTime()); //ms
#endif					
		
// #ifdef XXX		
// 		cout<<"isec: " << isec <<  "  n hits, dup, tr_c, trip: "
// 				<<istat_nhits<<" "<<istat_nduplets<<" "<<istat_n_triplet_c<<" "<<istat_n_triplets<<endl;
// #endif

					/// copy triplets in vTriplets
#ifdef XXX
    cout << " copy triplets in vTriplets " << endl;
#endif
#ifdef XXX 
		c_time_find[2][0].Start();
#endif
					
		vTriplets.clear();		
		
		for (int istal = NStations-2; istal >= FIRSTCASTATION; istal--){
			TripStartIndex[istal] = vTriplets.size();
			
			for (unsigned i = 0; i < vTriplets_part[istal].size(); i++){
				vTriplets.push_back(vTriplets_part[istal][i]);
			}

			for (int i = StsHitsStartIndex[istal]; i <=StsHitsStopIndex[istal]; i++){
				if (TripStopIndexH[i] > 0){
					TripStartIndexH[i]+=TripStartIndex[istal];
					TripStopIndexH[i] +=TripStartIndex[istal];
				};
			}
			
			TripStopIndex[istal] = vTriplets.size();
		}

#ifdef XXX  
    { // calculate occupied memory size

//       int sizeD, sizeN, sizeT, sizeF, sizeH, sizeTrip;
//       sizeD = sizeof(unsigned /*short*/ int)*MaxArrSize*MaxNStations*2; // Duplets_*,
//       sizeN = sizeof(int)*MaxNPortion*2; // n_g*
//       sizeT = sizeof(L1TrackPar)*MaxNPortion*Portion/fvecLen; // T_g1
//       sizeF = sizeof(L1FieldRegion)*MaxNPortion*Portion/fvecLen; // fld_g1
//       sizeH = sizeof(unsigned /*short*/ int)*MaxNPortion*(Portion + MaxPortionDoublets*2); // hitsl_g1, hitsm_g2,i1_g2
//       sizeTrip = sizeof(L1Triplet)*vTriplets.size();

      int sizeD2 = 0, sizeN2, sizeT2, sizeF2, sizeH2, sizeTrip2;
      int ndoublets = 0;
      for(int i = 0; i < NStations; i++){
        sizeD2 += Duplets_hits[i].size();
        sizeD2 += Duplets_start[i].size();
        ndoublets += Duplets_hits[i].size();
        sizeD2*=sizeof(unsigned /*short*/ int);
      }

      cout << "number of doublets: " << ndoublets << endl;
      cout << "number of triplets: " << vTriplets.size() << endl;
/*
      sizeN2 = ( n_g1.size() + n_g2.size() )*sizeof(int);
      sizeT2 = T_g1.CalcSize();
      sizeF2 = fld_g1.CalcSize();
      sizeH2 = ( hitsl_g1.CalcSize() + hitsm_g2.CalcSize() + i1_g2.CalcSize());

      cout << sizeD/1000 << " D " << sizeD2/1000 << " [1000elem]" << endl;
      cout << sizeN/1000 << " N " << sizeN2/1000 << " [1000elem]" << endl;
      cout << sizeT/1000 << " T " << sizeT2/1000 << " [1000elem]" << endl;
      cout << sizeF/1000 << " F " << sizeF2/1000 << " [1000elem]" << endl;
      cout << sizeH/1000 << " H " << sizeH2/1000 << " [1000elem]" << endl;

      long int size = sizeD2 + sizeN2 + sizeT2+ sizeF2 + sizeH2 + sizeTrip;
      cout << "Size: " << size/1000 << " KB" << endl;*/
    }
#endif

#ifdef XXX 	
		c_time_find[2][0].Stop();    
		time_find[2][0][isec] += double(c_time_find[2][0].CpuTime()); //ms
#endif							
   
					/// Find neighbours of triplets.
#ifdef XXX
    cout << " Find neighbours of triplets. " << endl;
#endif
#ifdef XXX 
		c_time_find[4][0].Start();
#endif
 
		int nlevel[NStations];	// number of triplets with some number of neighbours.
		for (int il = 0; il < NStations; ++il) nlevel[il] = 0;   
		
 		f5( 	// input 
				NStations, 
				TripStartIndex, TripStopIndex, 
					// output
				vTriplets,
				TripStartIndexH, TripStopIndexH,
				nlevel
			);
		
#ifdef XXX 	
		c_time_find[4][0].Stop();    
		time_find[4][0][isec] += double(c_time_find[4][0].CpuTime()); //ms
#endif		

#ifdef XXX    
    c_time_trip.Stop();    
    time_trip[isec] += double(c_time_trip.CpuTime()); //ms
#endif 

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

#ifdef XXX
    cout<<"---- Collect track candidates. ----"<<endl;
#endif   
#ifdef XXX    
    TStopwatch c_time_fit;
#endif    
    int min_level = 0;//1; //remove pile-up in MAPS //!!!
    if (isec == 1) min_level = 0; 
    if (isec == 2) min_level = 2; //only the longest tracks
    if (isec == -1) min_level = NStations-3 - 3; //only the longest tracks
    //min_level = 1; //!!!
    
    for (int ilev = NStations-3; ilev >= min_level; ilev--){

      //for (int ilev = NStations-3; ilev >= 0; ilev--){
#ifdef XXX 
      TStopwatch c_time_fit_bra;
#endif
      static vector<L1Branch> vtrackcandidate; vtrackcandidate.clear();
      
      //how many levels to check
      int nlevel = (NStations-2)-ilev+1;
      int ntracks = 0;

      for( int istaF = 0; istaF<= NStations-3-ilev; istaF++ ){
				
				//cout<<"istaF="<<istaF<<", max= "<<NStations-3-ilev<<endl;
				
				if( isec==2 && istaF>=4 ) break; // ghost supression !!!
				
				int sta_first = TripStartIndex[istaF];
				int sta_end   = TripStopIndex[istaF];
				//cout<<" level "<<ilev<<", station "<<istaF<<", triplets "<<sta_end-sta_first
				//<<", trip indices ("<<sta_first<<"->"<<sta_end<<")"<<endl;
				//cout<<svStsHits->size()<<endl;
				for( int itrip=sta_first; itrip<sta_end; itrip++ ){	  	
					L1Triplet *first_trip = &vTriplets[itrip];
					int lv=first_trip->GetLevel();
					int jst= GetFStation( vSFlag[(*svStsHits)[first_trip->GetLHit()].f] );
					int fl = GetFUsed( vSFlag[(*svStsHits)[first_trip->GetLHit()].f] | vSFlagB[(*svStsHits)[first_trip->GetLHit()].b] );
			
					if ( lv != ilev) continue;
					if (ilev == 0){ //remove pile-up -> collect only MAPS triplets
						if ( jst != 0) continue;
					}
					if( fl ) continue; // if used	  
					
					/*
					//check the level	  
					if (first_trip->GetLevel() != ilev) continue;
					if (ilev == 0){ //remove pile-up -> collect only MAPS triplets
						if ( GetFStation(vSFlag[svStsHits[first_trip->GetLHit()].f]) != 0) continue;
					}
			
					//check if all hits are not used => clear start
					
					if( GetFUsed( vSFlag[svStsHits[first_trip->GetLHit()].f] | vSFlagB[svStsHits[first_trip->GetLHit()].b] )
							) continue; // if used	  
					*/
			
					//cout<<itrip<<" "<<sta_end<<endl;
					
					//OK, search the best track candidate in the branches starting from the first_trip
			
					L1Branch curr_tr;
					curr_tr.StsHits.push_back(RealIHit[first_trip->GetLHit()]);	  
					unsigned char curr_L = 1;	
					fscal curr_chi2 = first_trip->GetChi2();
					//fscal curr_Qp = first_trip->GetQpOrig(MaxInvMom);
					//cout << "Triplet Qp " << trip.GetQpOrig(MaxInvMom) << " mom " << 1.0/trip.GetQpOrig(MaxInvMom) << endl;
			
					L1Branch best_tr = curr_tr;
					fscal best_chi2 = curr_chi2;
					unsigned char best_L = curr_L;
					int NCalls = 1;
					//cout<<"CAFindTrack, triplet N "<<itrip-sta_first<<" of "<<sta_end-sta_first<<endl;
					CAFindTrack((*svStsHits), RealIHit, istaF, first_trip, best_tr, best_L, best_chi2, curr_tr, curr_L, curr_chi2, NCalls);
					//cout<<"CAFindTrack ok, L="<<best_L<<", chi2="<<best_chi2/(best_L*2-5)
					//<<", Ncalls = "<<NCalls<<endl;
					//if(sta_first==34218 && ilev==1) cout<<4<<endl;
					if (best_L < ilev+2) continue; // maximum one hit missing
			
					int ndf = best_L*2-5;
					best_chi2 = best_chi2/ndf; //normalize
				
					if (best_chi2 > TRACK_CHI2_CUT) continue;
			    
					// fit track candidate
			
					if( best_L <3 ) continue;
			
					static L1TrackPar T _fvecalignment;
					//int nGathered = 0;
					
/*					
					if(0){ // fit forward
			
						int nHits = best_tr.StsHits.size();
			
						//T.qp = first_trip->GetQpOrig(MaxInvMom);
						//FitTrack( T, &(best_tr.StsHits[0]), best_tr.StsHits.size(), 1 );
						
						L1StsHit &hit0 = vStsHits[best_tr.StsHits[0]];
						L1StsHit &hit1 = vStsHits[best_tr.StsHits[1]];
						L1StsHit &hit2 = vStsHits[best_tr.StsHits[2]];
			
						int ista0 = GetFStation( vSFlag[hit0.f] );
						int ista1 = GetFStation( vSFlag[hit1.f] );
						int ista2 = GetFStation( vSFlag[hit2.f] );
			
						L1Station &sta0 = vStations[ista0];
						L1Station &sta1 = vStations[ista1];
						L1Station &sta2 = vStations[ista2];
			
						fvec x0  = vStsStrips[hit0.f];
						fvec v0  = vStsStripsB[hit0.b];
						fvec y0  = sta0.yInfo.cos_phi*x0 + sta0.yInfo.sin_phi*v0;
						fvec z0 = sta0.z;
			
						fvec x1  = vStsStrips[hit1.f];
						fvec v1  = vStsStripsB[hit1.b];
						fvec y1  = sta1.yInfo.cos_phi*x1 + sta1.yInfo.sin_phi*v1;
						fvec z1 = sta1.z;
			
						fvec x2  = vStsStrips[hit2.f];
						fvec v2  = vStsStripsB[hit2.b];
						fvec y2  = sta2.yInfo.cos_phi*x2 + sta2.yInfo.sin_phi*v2;
						fvec z2 = sta2.z;
					
						fvec dzi = 1./(z1-z0);
			
						fvec qp0 = first_trip->GetQpOrig(MaxInvMom);
			
						T.chi2 = 0.;
						T.NDF = 0.;
						T.x  = x0;
						T.y  = y0;
						T.tx = (x1-x0)*dzi;
						T.ty = (y1-y0)*dzi;
						T.qp = qp0;
						T.z  = z0;
						T.C00 = sta0.XYInfo.C00;
						T.C10 = sta0.XYInfo.C10;
						T.C11 = sta0.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.;
			
						static L1FieldValue fB0, fB1, fB2 _fvecalignment;
						static L1FieldRegion fld _fvecalignment;
						fvec fz0 = sta1.z;
						fvec fz1 = sta2.z;
						fvec fz2 = sta0.z;
						sta1.fieldSlice.GetFieldValue( x1, y1, fB0 );
						sta2.fieldSlice.GetFieldValue( x2, y2, fB1 );
						sta0.fieldSlice.GetFieldValue( x0, y0, fB2 );
			
						fld.Set( fB2, fz2, fB1, fz1, fB0, fz0 );
						int ista = ista2;
			
						for( int i=1; i<nHits; i++){
							L1StsHit &hit = vStsHits[best_tr.StsHits[i]];	      
							ista = vSFlag[hit.f]/4;
							L1Station &sta = vStations[ista];
							fvec x = vStsStrips[hit.f];
							fvec v = vStsStripsB[hit.b];
							fvec y  = sta.yInfo.cos_phi*x + sta.yInfo.sin_phi*v;
							L1Extrapolate( T, sta.z, qp0, fld );	
							L1Filter( T, sta.frontInfo, x );
							L1Filter( T, sta.backInfo,  v );
							L1AddMaterial( T, sta.materialInfo, qp0 );
							fB0 = fB1;
							fB1 = fB2;
							fz0 = fz1;
							fz1 = fz2;
							sta.fieldSlice.GetFieldValue( x, y, fB2 );
							fz2 = sta.z;
							fld.Set( fB2, fz2, fB1, fz1, fB0, fz0 );
						}
					
						// gather missed hits
						
						if(1){
							for( ista+=2; ista<NStations; ista++ ){
								L1Station &sta = vStations[ista];
								L1Extrapolate( T, sta.z, qp0, fld );
				
								fscal dxm_est = ( Pick_gather*sqrt(fabs(T.C00+sta.XYInfo.C00)) )[0];
								fscal dym_est = ( Pick_gather*sqrt(fabs(T.C11+sta.XYInfo.C11)) )[0];	
								fscal x_minus = T.x[0] - dxm_est;
								fscal x_plus  = T.x[0] + dxm_est;	  
								fscal y_minus = T.y[0] - dym_est;
								fscal y_plus  = T.y[0] + dym_est;
				
								double D2 = 1.e8;
								int ibest = -1;
								for( int ih = StsHitsStartIndex[ista]; ih<=StsHitsStopIndex[ista]; ih++ ){
							
									L1StsHit &hit = vStsHits[ih];
									if( GetFUsed( vSFlag[hit.f] | vSFlagB[hit.b] ) ) continue; // if used
							
									fscal x = vStsStrips[hit.f];
									fscal v = vStsStripsB[hit.b];
									fscal y  = (sta.yInfo.cos_phi*x + sta.yInfo.sin_phi*v)[0];
									if (y < y_minus) continue;
									if (y > y_plus ) break;	      
									if ((x < x_minus) || (x > x_plus)) continue;
									double dx = x-T.x[0];
									double dy = y-T.y[0];
									double d2 = dx*dx + dy*dy;
									if( d2>D2 ) continue;
									D2 = d2;
									ibest = ih;
								}
								if( ibest <0 ) break;
								
								best_tr.StsHits.push_back(ibest);
								best_L++;
								nGathered ++;
								L1StsHit &hit = vStsHits[ibest];
								fvec x = vStsStrips[hit.f];
								fvec v = vStsStripsB[hit.b];
								fvec y  = sta.yInfo.cos_phi*x + sta.yInfo.sin_phi*v;
								L1Filter( T, sta.frontInfo, x );
								L1Filter( T, sta.backInfo,  v );
								L1AddMaterial( T, sta.materialInfo, qp0 );
								fB0 = fB1;
								fB1 = fB2;
								fz0 = fz1;
								fz1 = fz2;
								sta.fieldSlice.GetFieldValue( x, y, fB2 );
								fz2 = sta.z;
								fld.Set( fB2, fz2, fB1, fz1, fB0, fz0 );
							}
						} // if
					}
			
					if(0){ // fit backward
						cout<<" back 1"<<endl;
						int nHits = best_tr.StsHits.size();
			
						L1StsHit &hit0 = vStsHits[best_tr.StsHits[nHits-1]];
						L1StsHit &hit1 = vStsHits[best_tr.StsHits[nHits-2]];
						L1StsHit &hit2 = vStsHits[best_tr.StsHits[nHits-3]];
			
						int ista0 = vSFlag[hit0.f]/4;
						int ista1 = vSFlag[hit1.f]/4;
						int ista2 = vSFlag[hit2.f]/4;
			
						cout<<"back: ista012="<<ista0<<" "<<ista1<<" "<<ista2<<endl;
						L1Station &sta0 = vStations[ista0];
						L1Station &sta1 = vStations[ista1];
						L1Station &sta2 = vStations[ista2];
			
						fvec x0  = vStsStrips[hit0.f];
						fvec v0  = vStsStripsB[hit0.b];
						fvec y0  = sta0.yInfo.cos_phi*x0 + sta0.yInfo.sin_phi*v0;
						fvec z0 = sta0.z;
			
						fvec x1  = vStsStrips[hit1.f];
						fvec v1  = vStsStripsB[hit1.b];
						fvec y1  = sta1.yInfo.cos_phi*x1 + sta1.yInfo.sin_phi*v1;
						fvec z1 = sta1.z;
			
						fvec x2  = vStsStrips[hit2.f];
						fvec v2  = vStsStripsB[hit2.b];
						fvec y2  = sta2.yInfo.cos_phi*x2 + sta2.yInfo.sin_phi*v2;
						fvec z2 = sta2.z;
					
						fvec dzi = 1./(z1-z0);
			
						cout<<"qp0 ="<<T.qp<<endl;
						fvec qp0 = first_trip->GetQpOrig(MaxInvMom);
			
						T.x  = x0;
						T.y  = y0;
						T.tx = (x1-x0)*dzi;
						T.ty = (y1-y0)*dzi;
						T.qp = qp0;
						T.z  = z0;
						T.chi2 = 0.;
						T.NDF = 0.;
						T.C00 = sta0.XYInfo.C00;
						T.C10 = sta0.XYInfo.C10;
						T.C11 = sta0.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.;
			
						static L1FieldValue fB0, fB1, fB2 _fvecalignment;
						static L1FieldRegion fld _fvecalignment;
						fvec fz0 = sta1.z;
						fvec fz1 = sta2.z;
						fvec fz2 = sta0.z;
						sta1.fieldSlice.GetFieldValue( x1, y1, fB0 );
						sta2.fieldSlice.GetFieldValue( x2, y2, fB1 );
						sta0.fieldSlice.GetFieldValue( x0, y0, fB2 );
			
			
						fld.Set( fB2, fz2, fB1, fz1, fB0, fz0 );
						cout<<" back 2"<<endl;
						int ista = ista2;
						for( int i=nHits-2; i>=0; i--){
							if( fabs(T.qp[0])>2. ) break;
							L1StsHit &hit = vStsHits[best_tr.StsHits[i]];
							ista = vSFlag[hit.f]/4;
							cout<<i<<" :ista="<<ista<<endl;
							L1Station &sta = vStations[ista];
							fvec x = vStsStrips[hit.f];
							fvec v = vStsStripsB[hit.b];
							fvec y  = sta.yInfo.cos_phi*x + sta.yInfo.sin_phi*v;
							cout<<i<<" :"<<x<<" "<<y<<" "<<v<<endl;
							L1Extrapolate( T, sta.z, qp0, fld );	
							cout<<i<<":1 "<<endl;
							L1Filter( T, sta.frontInfo, x );
							cout<<i<<":2 "<<endl;
							L1Filter( T, sta.backInfo,  v );
							cout<<i<<":3 "<<T.tx<<" "<<T.ty<<endl;
							L1AddMaterial( T, sta.materialInfo, qp0 );
							cout<<i<<":4  "<<endl;
							fB0 = fB1;
							fB1 = fB2;
							fz0 = fz1;
							fz1 = fz2;
							sta.fieldSlice.GetFieldValue( x, y, fB2 );
							fz2 = sta.z;
							fld.Set( fB2, fz2, fB1, fz1, fB0, fz0 );
							cout<<i<<" ok "<<endl;
						} // i
					
						cout<<" back 3"<<endl;
						// gather missed hits
						
						if(1){
							for( ista-=2; ista>=0; ista-- ){
				
								L1Station &sta = vStations[ista];
				
								L1Extrapolate( T, sta.z, qp0, fld );
				
								fscal dxm_est = ( Pick_gather*sqrt(fabs(T.C00+sta.XYInfo.C00)) )[0];
								fscal dym_est = ( Pick_gather*sqrt(fabs(T.C11+sta.XYInfo.C11)) )[0];	
								fscal x_minus = T.x[0] - dxm_est;
								fscal x_plus  = T.x[0] + dxm_est;	  
								fscal y_minus = T.y[0] - dym_est;
								fscal y_plus  = T.y[0] + dym_est;
				
								double D2 = 1.e8;
								int ibest = -1;
								for( int ih = StsHitsStartIndex[ista]; ih<=StsHitsStopIndex[ista]; ih++ ){
							
									L1StsHit &hit = vStsHits[ih];
									if( GetFUsed( vSFlag[hit.f] | vSFlagB[hit.b] ) ) continue; // if used
							
									fscal x = vStsStrips[hit.f];
									fscal v = vStsStripsB[hit.b];
									fscal y  = (sta.yInfo.cos_phi*x + sta.yInfo.sin_phi*v)[0];
									if (y < y_minus) continue;
									if (y > y_plus ) break;	      
									if ((x < x_minus) || (x > x_plus)) continue;
									double dx = x-T.x[0];
									double dy = y-T.y[0];
									double d2 = dx*dx + dy*dy;
									if( d2>D2 ) continue;
									D2 = d2;
									ibest = ih;
								}
								if( ibest <0 ) break;
				
								best_tr.StsHits.push_back(ibest);
								best_L++;
								nGathered ++;
								L1StsHit &hit = vStsHits[ibest];
								fvec x = vStsStrips[hit.f];
								fvec v = vStsStripsB[hit.b];
								fvec y  = sta.yInfo.cos_phi*x + sta.yInfo.sin_phi*v;
								L1Filter( T, sta.frontInfo, x );
								L1Filter( T, sta.backInfo,  v );
								L1AddMaterial( T, sta.materialInfo, qp0 );
								fB0 = fB1;
								fB1 = fB2;
								fz0 = fz1;
								fz1 = fz2;
								sta.fieldSlice.GetFieldValue( x, y, fB2 );
								fz2 = sta.z;
								fld.Set( fB2, fz2, fB1, fz1, fB0, fz0 );
							} // ista
						} // if 
					}
	*/				
		
					//if(  nGathered>0 ) cout<<nGathered<<" "<<ndf<<" "<<T.chi2[0]<<endl;
					//ndf = best_L*2-5;
					//if( !finite(T.chi2[0]) ) continue;
					//best_chi2 = T.chi2[0]/ndf; //normalize
					//if ( T.chi2[0] > TRACK_CHI2_CUT *ndf ) continue;
					//if ( fabs(T.qp[0]) > MaxInvMom ) continue;
					
					if( best_L < min_level+3 ) continue;
						
					if( fGhostSuppression ){//suppress ghost	    
						if( best_L <3 ) continue;
						if( best_L == 3 ){
							//if( isec == 2 ) continue; // too /*short*/ secondary track
							if( isec==2 && istaF != 0 ) continue; // too /*short*/ non-MAPS track
							if( (isec<2)&&(best_chi2>5.0) ) continue; 
						} 
					}
					// track candidate is OK, store
			
					best_tr.Set( istaF, best_L, best_chi2, first_trip->GetQpOrig(MaxInvMom));

					vtrackcandidate.push_back(best_tr);
					ntracks++;	
				}  // itrip
				//cout<<" level "<<ilev<<", station "<<ista<<" ok"<<endl;
      } // istaF
      
      //cout <<"Level "<< ilev <<" track candidates "<< ntracks << ", ";
      //cout << " total number of track candidates " << vtrackcandidate.size() <<endl;
      
      if (--nlevel == 0) break;    
#ifdef XXX      
      c_time_fit_bra.Stop();    
      time_fit_bra[isec] += double(c_time_fit_bra.CpuTime());
#endif
      //make competition between tracks of the same length
#ifdef XXX 
      TStopwatch c_time_fit_fin;
#endif
      // 1. sort tracks
      static vector<L1Branch*> vptrackcandidate;  // vptrackcandidate - array of pointers to vtrackcandidate
			vptrackcandidate.clear();

      for (vector<L1Branch>::iterator trIt = vtrackcandidate.begin();
	   																	trIt != vtrackcandidate.end(); ++trIt){
				vptrackcandidate.push_back(&(*trIt));
      }
      //if (isec <= 1) 
			sort(vptrackcandidate.begin(), vptrackcandidate.end(), L1Branch::comparePChi2);
			//else
			//sort(vptrackcandidate.begin(), vptrackcandidate.end(), L1Branch::comparePMomentum);
				
      ntracks = 0;

      for (vector<L1Branch*>::iterator 	trIt = vptrackcandidate.begin();
																				trIt != vptrackcandidate.end(); ++trIt){
				
				//count number of used hits
				
				L1Branch *tr = *trIt;
			
				//cout << "Branch isec " << isec << " Nhits " << tr->StsHits.size() << " Quality " << tr->Quality  << " chi2 " << tr-> chi2 << " momentum " << tr->Momentum << endl;
				int nused = 0;
				for (vector<unsigned /*short*/ int>::iterator phIt = tr->StsHits.begin();
						phIt != tr->StsHits.end(); ++phIt){
					L1StsHit &h = vStsHits[*phIt];
					if ( GetFUsed( vSFlag[h.f] | vSFlagB[h.b] ) ) nused++; 
				}
				
				if (nused != 0) continue;
				
				//=======================================================
				// gather MAPS hits using the Kalman filter
				// was here, skipped for a moment
				//=======================================================
				
				//store track
				for (vector<unsigned /*short*/ int>::iterator phitIt = tr->StsHits.begin();
																									phitIt != tr->StsHits.end(); ++phitIt){	  
					L1StsHit &h = vStsHits[*phitIt];
					SetFUsed( vSFlag[h.f] );
					SetFUsed( vSFlagB[h.b] );
					vRecoHits.push_back(*phitIt);
				}
				
				//store track
				L1Track t;
				t.NHits = tr->StsHits.size();
				t.Momentum = tr->Momentum;
				for( int i=0; i<6; i++) t.TFirst[i] = t.TLast[i]=0; 
				for( int i=0; i<15; i++ ) t.CFirst[i] = t.CLast[i] = 10;
				t.chi2 = 100;
			
				vTracks.push_back(t);
				ntracks++;   

      } //trIt
#ifdef XXX    
      c_time_fit_fin.Stop();    
      time_fit_fin[isec] += double(c_time_fit_fin.CpuTime());
#endif  
      unsigned int Bsize=0, nb=0;
      for (vector<L1Branch*>::iterator 	trIt = vptrackcandidate.begin();
	   																		trIt != vptrackcandidate.end(); ++trIt){
				L1Branch *tr = *trIt;
				Bsize+= sizeof(L1Branch) + sizeof(unsigned /*short*/ int) * tr->StsHits.size();
				nb++;
      }
#ifdef XXX 			
      if( Bsize>stat_max_BranchSize ) stat_max_BranchSize = Bsize;       
      if( nb>stat_max_n_branches ) stat_max_n_branches = nb;
#endif      
      //for (vector<L1StsHit>::iterator ih = svStsHits->begin(); ih != svStsHits->end(); ++ih){
      //if( ih->strip_f->used || ih->strip_b->used) ih->used = 1;
      //}
    
      //cout<< " track candidates stored: "<< ntracks <<endl;
    } //ilev

      // clear hits array
#ifdef XXX
    cout << " clear of hits arrays " << endl;
#endif

		int StsHitsStartIndex_temp;
		int ista = 0;
		int nDontUsedHits = 0;
		svStsHits_buf->clear();
		svStsHits_buf2->clear();
		for(; ista < NStations; ista++){
		
			StsHitsStartIndex_temp = StsHitsStartIndex[ista];
			StsHitsStartIndex[ista] = nDontUsedHits;
			for(int ih = StsHitsStartIndex_temp; ih <= StsHitsStopIndex[ista]; ih++){
				int rih = RealIHit[ih];
				
				L1StsHit hit = vStsHits[rih];
				if( GetFUsed( vSFlag[hit.f] | vSFlagB[hit.b] ) ){ continue;} // if used
				svStsHits_buf->push_back(hit);
				svStsHits_buf2->push_back((*svStsHits2)[ih]);
				
// cout << RealIHit[ih]  << " " << ih  << " " << nDontUsedHits << endl;
				RealIHit[nDontUsedHits] = rih;
				nDontUsedHits++;
			}
			StsHitsStopIndex[ista] = nDontUsedHits-1;
// 			int kk; std::cin >> kk;
		}
		std::vector< L1StsHit > *svStsHits_temp = svStsHits;
		svStsHits = svStsHits_buf;
		svStsHits_buf = svStsHits_temp;
		std::vector< L1HitPoint > *svStsHits_temp2 = svStsHits2;
		svStsHits2 = svStsHits_buf2;
		svStsHits_buf2 = svStsHits_temp2;				
// 		cout << vStsHits.size() << " " << svStsHits->size() << " " << svStsHits_buf->size() << endl;		
		svStsHits_buf->clear();
		svStsHits_buf2->clear();		
// cout << " clear hits array end " << endl;


//     cout<<"End of collecting track candidates "<<endl;
#ifdef XXX
    c_time_fit.Stop();    
    time_fit[isec] += double(c_time_fit.CpuTime());
#endif
#ifdef XXX
//     if( stat_max_n_trip<stat_n_trip ) stat_max_n_trip = vTriplets.size();
    unsigned int tsize = vTriplets.size()*sizeof(L1Triplet);
    if( stat_max_trip_size < tsize ) stat_max_trip_size = tsize;
#endif
  } // for (int isec

  //cout<< "Total number of tracks found " << vTracks.size() <<endl;
  
  //==================================

  c_time.Stop();						

// 	cout << "End TrackFinder" << endl; 
  CATime = (double(c_time.CpuTime()));
	double CATimeReal = (double(c_time.RealTime()));
	
#ifdef XXX
  /*
  cout<<"CA Track Finder: " << (double(c_time.CpuTime()))
      <<" (" << time_trip <<" [" 
      << time_find[1][0] << "/"<< time_find[1][1]<< "+"
      << time_find[1][0] << "/"<< time_find[1][1]<< "+"
      << time_find[2][0] << "/"<< time_find[2][1]<< "+"
      << time_find[4][0] << "+"
      << time_find[4][0] << "+"
      << "]"
      << " + " << time_fit <<" [" << time_fit_bra << "/" << time_fit_fin << "]"
      << " ) s/ev"<<endl;
  */
  stat_N++;
  stat_time     += CATime;	
	stat_time_real += CATimeReal;

	
	cout.precision(3);
	cout<<"CA stat time: " << (stat_time/stat_N) << "/" << stat_time_real/stat_N << "(k=" << stat_time/stat_time_real <<") s/ev"<<endl;
	
	for (int isec = 0; isec < Nisec; isec++){
		stat_time_trip[isec] += time_trip[isec];
		
		for (int i1 = 0; i1 < ntimers; i1++)
			for (int i2 = 0; i2 < 2; i2++)
				stat_time_find[i1][i2][isec] += time_find[i1][i2][isec];
	
		stat_time_fit[isec] += time_fit[isec];
		stat_time_fit_bra[isec] += time_fit_bra[isec];
		stat_time_fit_fin[isec] += time_fit_fin[isec];
	
  
		
	cout  << isec << " isec : 	"
				<<" (" << stat_time_trip[isec]/stat_N <<" [" 
					<< stat_time_find[1][0][isec]/stat_N  <<"/" << stat_time_find[1][1][isec]/stat_N 
					<<	"(k=" << stat_time_find[1][0][isec]/stat_time_find[1][1][isec] <<")+"  
					<< stat_time_find[2][0][isec]/stat_N  <<"/" << stat_time_find[2][1][isec]/stat_N
					<<	"(k=" << stat_time_find[1][0][isec]/stat_time_find[1][1][isec] <<")+" 
					<< stat_time_find[2][0][isec]/stat_N  <<"/" << stat_time_find[2][1][isec]/stat_N <<" + "
					<< stat_time_find[3][0][isec]/stat_N <<" + "
					<< stat_time_find[4][0][isec]/stat_N 
					<< "]"
				<< " + " << stat_time_fit[isec]/stat_N 
					<<" [" << stat_time_fit_bra[isec]/stat_N << "\\" << stat_time_fit_fin[isec]/stat_N << "]"
				<< ") s/ev"<<endl;
				
	stat_time_find_all[1][0] += stat_time_find[1][0][isec];
	stat_time_find_all[1][1] += stat_time_find[1][1][isec];
	stat_time_find_all[2][0] += stat_time_find[2][0][isec];
	stat_time_find_all[2][1] += stat_time_find[2][1][isec];
	stat_time_find_all[3][0] += stat_time_find[3][0][isec];	
	stat_time_find_all[3][1] += stat_time_find[4][0][isec];
	stat_time_find_all[4][0] += stat_time_fit[isec];
  } // isec
	
	for(int i2 = 0; i2 < ntimers; i2++)
		for(int i3 = 0; i3 < 2; i3++)
			stat_time_find_all[i2][i3] /= stat_N;		
	
	cout  << "all : 	"
				<< stat_time_find[0][0] << "+"
				<<" (" << "	[" 
					<< stat_time_find_all[1][0]  <<"/" << 	stat_time_find_all[1][1] << "(k=" << stat_time_find_all[1][0]/stat_time_find_all[1][1] <<")+" 
					<< stat_time_find_all[2][0]  <<"/" << 	stat_time_find_all[2][1] << "(k=" << stat_time_find_all[2][0]/stat_time_find_all[2][1] <<")+" 
					<< stat_time_find_all[3][0]  <<"+"
					<< stat_time_find_all[3][1] 
					<< "]"
				<< " + " << stat_time_find_all[4][0] 
				<< ") s/ev"<<endl;  
	
	
  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 += svStsHits->size();
  HitSize += svStsHits->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(unsigned /*short*/ int)*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;
  
  //cout << "===> ... CA Track Finder    " << endl << endl;
#endif

#ifdef DRAW
  DrawAsk();
#endif 
  

}

//1234567890123456789012345678901234567890123456789012345678901234567890


void L1Algo::CAFindTrack(std::vector< L1StsHit > &svStsHits, unsigned int *RealIHit, int ista, const L1Triplet* curr_trip, 
			  L1Branch& best_tr, unsigned char &best_L, fscal &best_chi2,
			  L1Branch &curr_tr, unsigned char &curr_L, fscal &curr_chi2,
			  int &NCalls ){

  /****************************************************************
   *                                                              *
   *     The routine performs recusrsive search for tracks        *
   *                                                              *
   *     I. Kisel                                    06.03.05     *
   *                                                              *
   ***************************************************************/

  //  cout << "CA track finder " << curr_trip->level << endl;

  if (curr_trip->GetLevel() == 0){ // the last triplet -> check and store    
    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;
    } 
  } 
  else {// curr_trip->level != 0
    
    //store hits & triplets of the current track
    L1Branch temp_tr = curr_tr;
    unsigned char temp_L = curr_L;
    fscal temp_chi2 = curr_chi2;
    int offset = TripStartIndex[ista+1];
    int first_neighbour = offset + curr_trip->GetFirstNeighbour();
    int end_neighbour = first_neighbour + curr_trip->GetNNeighbours();

    for(int index = first_neighbour; index < end_neighbour; ++index){

      L1Triplet *new_trip = &vTriplets[index];      

      if (curr_trip->GetLevel() != new_trip->GetLevel()+1) continue;
      fscal qp1 = curr_trip->GetQp();
      fscal qp2 = new_trip->GetQp();
      fscal dqp = fabs(qp1 - qp2);
      fscal Cqp = curr_trip->Cqp; 
      Cqp      += new_trip->Cqp;
      if ( dqp > Cqp  ) continue; //bad neighbour
      
      //no used hits allowed
      if ( GetFUsed( vSFlag[svStsHits[new_trip->GetLHit()].f] | vSFlagB[svStsHits[new_trip->GetLHit()].b]
		     ) )
	{
	  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;
	  } //compare and store track
	}
      else{//add this triplet to the current track

	// restore current track 

	curr_tr = temp_tr;
	curr_chi2 = temp_chi2;
	curr_L = temp_L;
	curr_tr.StsHits.push_back(RealIHit[new_trip->GetLHit()]); 
	curr_L+= 1;
	dqp = dqp/Cqp*5.;
	curr_chi2 += dqp*dqp;
	if( new_trip->GetLevel()==0 ){
	  unsigned /*short*/ int ihitm = new_trip->GetMHit();
	  unsigned /*short*/ int ihitr = new_trip->GetRHit();
	  if( !GetFUsed( vSFlag[svStsHits[ihitm].f] |  vSFlagB[svStsHits[ihitm].b] )  ){
	    curr_tr.StsHits.push_back(RealIHit[ihitm]);
	    curr_L+= 1;
	  }
	  if( !GetFUsed( vSFlag[svStsHits[ihitr].f] | vSFlagB[svStsHits[ihitr].b] ) ){
	    curr_tr.StsHits.push_back(RealIHit[ihitr]);
	    curr_L+= 1;
	  }
	}

	if( curr_chi2 <= TRACK_CHI2_CUT * (curr_L) ){ // go further      
	  CAFindTrack(svStsHits, RealIHit, ista+1,new_trip, best_tr, best_L, best_chi2, curr_tr, curr_L, curr_chi2,NCalls);
	  NCalls++;
	}
      }//good triplet added
    }//first_triplet --- last_triplet
  }//level != 0

  //cout << "Track finder finished " << endl;

  return;
}