// -----------------------------------------------------------------------
// -----                    CbmL1CATrdTrackFinderSA                  -----
// -----            Created 2/12/2006  by A. Bubak & M. Krauze       -----
// -----------------------------------------------------------------------

#include <iostream>

#include "TClonesArray.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TLinearFitter.h"

#include "CbmRootManager.h"
#include "CbmRunAna.h"
#include "CbmBaseParSet.h"
#include "CbmRuntimeDb.h"
#include "CbmDetector.h"
#include "CbmTrdPoint.h"
#include "CbmTrdHit.h"
#include "CbmStsTrack.h"
#include "CbmTrdTrack.h"
#include "CbmL1CATrdTrackFinderSA.h"
#include "CbmMCTrack.h"
#include "CbmMCPoint.h"
#include "CbmL1.h"
#include "CbmL1TrdTracklet.h"
#include "CbmL1TrdTracklet4.h"
#include "CbmGeoTrdPar.h"

#include "CbmKF.h"
#include "CbmKFTrdHit.h"

// ----------------------- Default constructor ---------------------------
CbmL1CATrdTrackFinderSA::CbmL1CATrdTrackFinderSA()
{
  fArrayStsTrack = NULL;
  fArrayKFTrdHit = new TClonesArray("CbmKFTrdHit");
  fVerbose = 1;

  fNofEvents = 0;

  fTrdSegments = new TClonesArray("CbmTrdTrack");
  fArrayTrdHit = new TClonesArray("CbmTrdHit");
  fArrayTrdTrack = new TClonesArray("CbmTrdTrack");
  trdTrackFitter = new CbmTrdTrackFitterKF();
}
// -----------------------------------------------------------------------

// --------------------------- Destructor --------------------------------
CbmL1CATrdTrackFinderSA::~CbmL1CATrdTrackFinderSA()
{
  delete fArrayKFTrdHit;
  delete fTrdSegments;
  delete fArrayTrdHit;
  delete fArrayTrdTrack;
  delete trdTrackFitter;
}
// -----------------------------------------------------------------------



// ------------------------- Initialisation ------------------------------
void CbmL1CATrdTrackFinderSA::Init()
{
  // Create histogramms
  CreateHistogramms();
  // Activate data branches
  DataBranches();

  CbmRootManager* ioman = CbmRootManager::Instance();
  if (! ioman) {
    cout << "-E- CbmL1CATrdTrackFinderSA::Init: "
	 << "RootManager not instantised!" << endl;
    return;
  }

  // Get MCTrack array
  fMCTrackArray  = (TClonesArray*) ioman->ActivateBranch("MCTrack");
  if ( ! fMCTrackArray) {
    cout << "-E- CbmL1CATrdTrackFinderSA::Init: No MCTrack array!"
	 << endl;
    return;
  }

  // Get MCPoint array
  fMCPointArray  = (TClonesArray*) ioman->ActivateBranch("TRDPoint");
  if ( ! fMCPointArray) {
    cout << "-E- CbmL1CATrdTrackFinderSA::Init: No MCPoint array!"
	 << endl;
    return;
  }

  // Determine the TRD layout
  TrdLayout();
  trdTrackFitter->Init();
}
// -----------------------------------------------------------------------



// -------------------- Algorithm implementation -------------------------
Int_t CbmL1CATrdTrackFinderSA::DoFind(TClonesArray *hitArray,
				      TClonesArray *trackArray)
{

  fArrayTrdHit = hitArray;
  fArrayTrdTrack = trackArray;

  fNofEvents++;

  doFind.Start();

  // Initialize control counters
  //  Int_t nNoMCTrack   = 0;
  //  Int_t nNoTrdPoint  = 0;
  Int_t nNoTrdHit    = 0;
  Int_t trdPlaneID   = 0;
  
  // Create pointers to TRD hit and TrdPoint
  CbmTrdHit*       pHit = NULL;
  CbmTrdPoint*     pMCpt = NULL;
  CbmMCTrack*      pMCtr = NULL;

  //--- Function: FindNeighbour ------------------------------------------
  Double_t
    distPropLongX_FL = 3,       //distance from X-propagated point, around which we look for neighbours
    distPropLongY_FL = 3,       //distance from Y-propagated point, around which we look for neighbours
    distPropLongX_SL = 8,   //the same but for low-momentum tracks
    distPropLongY_SL = 8;

  //--- Function: CreateSpacePoint --------------------------------------
  Double_t
    sigmaA_FL = 3,
    sigmaB_FL = 3,
    sigmaA_SL = 4,
    sigmaB_SL = 4;

  //--- Function: CreateSegments; Data from MC, look inside function ----
  Double_t 
    dY_FL = 0.5,
    dX_FL = 0.5,
    dY_SL = 1.0,
    dX_SL = 1.0;

  //----------------------------------------------------------------------
  Bool_t competition    = true;
  Bool_t removeUsedHits = true;
  Bool_t secondLoop     = true;
  Int_t segment = 0; //0 = normal, 1 = segments only, 2 = unused yet



  Int_t ptIndex = 0;       // MC point index

  sortHits.Start();

  //--- hit sorting procedure ------------------------------------------------
  Int_t nTrdHit = hitArray->GetEntriesFast();
  for (Int_t iHit = 0; iHit<nTrdHit; iHit++) {
    // Loop over hits. Get corresponding MCPoint and MCTrack index  
    pHit = (CbmTrdHit*) hitArray->At(iHit);
    if ( ! pHit ) {
      cout << "-E- CbmL1CATrdTrackFinderSA::DoFind: Empty slot "
	   << "in HitArray at position " << iHit << endl;
      nNoTrdHit++;
      continue;
    }

    ptIndex = pHit->GetRefIndex();    
    //     if (ptIndex < 0) continue;           // fake or background hit
    pMCpt = (CbmTrdPoint*) (fMCPointArray->At(ptIndex));
    //     if (!pMCpt) {
    //       nNoTrdPoint++;
    //       continue;
    //     }
    pMCtr = (CbmMCTrack*) fMCTrackArray->At(pMCpt->GetTrackID());
    //     if ( ! pMCtr ) continue;

    trdPlaneID = pHit->GetPlaneID();
    Int_t trdPlaneIDN = trdPlaneID-1;

    planeHits.mcTrackID = pMCpt->GetTrackID();
    planeHits.hitIndex = iHit;
    planeHits.X  = pHit->GetX();
    planeHits.Y  = pHit->GetY();
    planeHits.Z  = pHit->GetZ();
    planeHits.DX = pHit->GetDx()/1E+4;
    planeHits.DY = pHit->GetDy()/1E+4;
    planeHits.planeID = trdPlaneID;
    fvTrdHitArr[trdPlaneIDN].push_back(planeHits);

    //     ptIndex = pHit->GetRefIndex();
    //     pMCpt = (CbmTrdPoint*) (fMCPointArray->At(ptIndex));

    //     planeHits.hitIndex = ptIndex;
    //     planeHits.X  = pMCpt->GetX();
    //     planeHits.Y  = pMCpt->GetY();
    //     planeHits.Z  = pMCpt->GetZ();
    //     planeHits.DX = 0;
    //     planeHits.DY = 0;;
    //     fvTrdPointArr[trdPlaneIDN].push_back(planeHits);
  }

  for(Int_t i = 0; i < 12; i++){
    sort(fvTrdHitArr[i].begin(),  fvTrdHitArr[i].end(),  CompareY);
  }

  //-----------------------------------------------------------
//   vector <LayerWithHits>::iterator ivT;
//   LayerWithHits pl;
//   for(ivT = fvTrdHitArr[0].begin();
//       ivT != fvTrdHitArr[0].end();
//       ivT++){
//     pl = (*ivT);
//     cout << pl.Y << endl;
//     fPlane1Ydens->Fill(pl.Y);
//   }
//   for(ivT = fvTrdHitArr[5].begin();
//       ivT != fvTrdHitArr[5].end();
//       ivT++){
//     pl = (*ivT);
//     //    cout << pl.Y << endl;
//     fPlane5Ydens->Fill(pl.Y);
//   }
//   for(ivT = fvTrdHitArr[9].begin();
//       ivT != fvTrdHitArr[9].end();
//       ivT++){
//     pl = (*ivT);
//     //    cout << pl.Y << endl;
//     fPlane9Ydens->Fill(pl.Y);
//   }
  //-----------------------------------------------------------

  sortHits.Stop();
  
  if(fVerbose > 2) {    
    cout <<"-I- CbmL1CATrdTrackFinderSA::DoFind " << endl;
    for(Int_t i = 0; i < 12; i++){
      cout <<" Size of vector "<< i <<": "<< fvTrdHitArr[i].size() << endl;
    }
  }

  //create vectors that will hold the tracklet (4-hit) objects

  vector<CbmL1TrdTracklet4*> clTracklets[3];
  vector<CbmL1TrdTracklet4*>::iterator itclTrackletsRight;

  //create vectors that will hold the mini-tracklet (2-hit) objects - Space Points
  vector<CbmL1TrdTracklet*> clSpacePoints[6];

  Int_t noHitsPerTracklet = 4;
  
  //###################################################################################################

  if(noHitsPerTracklet == 4) { //noHitsPerTracklet == 4
    //------ creation of mini tracklet 1-2 -------------------------------------
    cout << endl << "-I- CbmL1CATrdTrackFinderSA::DoFind " << endl
	 << "--------------------------------------------------" << endl
	 << "     ### Start creation of Space Points" <<endl
	 << "--------------------------------------------------" << endl;

    createSPs.Start();    
    for(Int_t i = 0, j = 0; i < 6; i++, j=j+2) {
      CreateSpacePoints(fvTrdHitArr[j], fvTrdHitArr[j+1], 
			clSpacePoints[i], sigmaA_FL, sigmaB_FL);
    }
    createSPs.Stop();

//     Double_t rtime;
//     Double_t ctime;
//     cout << "Create SPs: ";
//     rtime = createSPs.RealTime();
//     ctime = createSPs.CpuTime();
//     printf("RealTime=%4.2f s, CpuTime=%4.2f s",rtime,ctime);

    
    //------- Sorting spacepoints -------------------------------------------------
    sortSPs.Start();            
    for(Int_t i = 0; i < 6; i++){
      sort(clSpacePoints[i].begin(), clSpacePoints[i].end(), CbmL1TrdTracklet::Compare1);
    }
    sortSPs.Stop();  



    // --------------- Creation tracklet from 4 space points ----------------
    cout << endl << "-I- CbmL1CATrdTrackFinderSA::DoFind " << endl
	 << "--------------------------------------------------" << endl
	 << "     ### Start creation tracklets" <<endl
	 << "--------------------------------------------------" << endl;
    createSegments.Start();    
    //joining SP together to create 4-hit tracklets
    for(Int_t i = 0, j = 0; i < 3; i++, j=j+2) {
      CreateSegments(clSpacePoints[j], clSpacePoints[j+1], 
		     clTracklets[i], dX_FL, dY_FL);
    }
    createSegments.Stop();



    //---------- BEGIN: Find friend ------------------------------------------------
    findNeighbour.Start();
    
    cout << "--- Finding neighbour 14-58" << endl;
    FindNeighbour(clTracklets[0],clTracklets[1],
		  distPropLongX_FL, distPropLongY_FL);

    cout << "--- Finding neighbour 58-912" << endl;
    FindNeighbour(clTracklets[1],clTracklets[2],
		  distPropLongX_FL, distPropLongY_FL);
    
    findNeighbour.Stop();

  }//end of noHitsPerTracklet == 4

  //###################################################################################################


  cout << endl << "-I- CbmL1CATrdTrackFinderSA::DoFind " << endl
       << "--------------------------------------------------" << endl
       << "     ### (FL) in Event No. " << fNofEvents << " ###" <<endl
       << "--------------------------------------------------" << endl;

  if(fVerbose > 2){
    cout    << "size of segment vector 14 = "
	    << clTracklets[0].size()
	    << endl
	    << "size of segment vector 58 = "
	    << clTracklets[1].size() 
	    << endl
	    << "size of segment vector 912 = "
	    << clTracklets[2].size() 
	    << endl;
  }
  
  if(fVerbose > 1) 
    cout << "-I- CbmL1CATrdTrackFinderSA::DoFind " << endl
	 << "Tracklet finding phase completed." << endl
	 << "Now constructing tracks from tracklets..."
	 << endl;


  //  CbmL1TrdTracklet4* clSegLeft; //segment farther from the target
  CbmL1TrdTracklet4* clSegRight; //segment nearer to the target
  
  tagTracks.Start();

  if(noHitsPerTracklet == 4) {//noHitsPerTracklet == 4
    //BEGIN: Tagging procedure ------------------------------------------------------
    cout << endl << "-I- CbmL1CATrdTrackFinderSA::DoFind " << endl
	 << "--------------------------------------------------" << endl
	 << "     ### Starting tagging procedure" <<endl
	 << "--------------------------------------------------" << endl;
    
    TagSegments(clTracklets[1], clTracklets[2]);
    cout << "--- Tagging 58-921 done." << endl;

    TagSegments(clTracklets[0], clTracklets[1]);
    cout << "--- Tagging 14-58 done." << endl;  
  } //end of noHitsPerTracklet == 4
 
  tagTracks.Stop();
  //------------------------------------------------------------------------- 


  //creating STL maps to hold used hits, i.e. hits that belong to an already reconstructed track
  map<Int_t, Int_t> usedHits9;
  map<Int_t, Int_t> usedHits12;

  map<Int_t, Int_t> segValues14;
  map<Int_t, Int_t> segValues58;
  map<Int_t, Int_t> segValues912;

  Int_t segValue;

  //--------------------------------------------------------------------------------------
  //    Counting number of segments which have a value of 4, 3, 2, ... etc.

  //for tracklets in the 1st station
  for(itclTrackletsRight = clTracklets[0].begin();
      itclTrackletsRight != clTracklets[0].end();
      itclTrackletsRight++) {
      
    clSegRight = *itclTrackletsRight;
    segValue = clSegRight->GetVal();
    segValues14[segValue]++;
  }

  //for tracklets in the 2nd station
  for(itclTrackletsRight = clTracklets[1].begin();
      itclTrackletsRight != clTracklets[1].end();
      itclTrackletsRight++) {
      
    clSegRight = *itclTrackletsRight;
    segValue = clSegRight->GetVal();
    segValues58[segValue]++;
  }

  //for tracklets in the 3rd station
  for(itclTrackletsRight = clTracklets[2].begin();
      itclTrackletsRight != clTracklets[2].end();
      itclTrackletsRight++) {
      
    clSegRight = *itclTrackletsRight;
    segValue = clSegRight->GetVal();
    segValues912[segValue]++;
  }


  map<Int_t,Int_t>::iterator mIt;

  if(fVerbose > 1) {
    cout << "--- Map no. 14 has: " << endl;
    for(mIt = segValues14.begin();
	mIt != segValues14.end();
	mIt++) {
      cout << mIt->first << " segment has a count of " 
	   << mIt->second << endl;
    }

    cout << "--- Map no. 58 has: " << endl;
    for(mIt = segValues58.begin();
	mIt != segValues58.end();
	mIt++) {
      cout << mIt->first << " segment has a count of " 
	   << mIt->second << endl;
    }
    
    cout << "--- Map no. 912 has: " << endl;
    for(mIt = segValues912.begin();
	mIt != segValues912.end();
	mIt++) {
      cout << mIt->first << " segment has a count of " 
	   << mIt->second << endl;
    }
    
    cout << endl;
  }

  //---------------------------------------------------------------------------------
      
  //vector<CbmTrdTrack*> vTrdTrackCand;
  //vector<CbmTrdTrack*> vTempTrdTrackCand;
  //  vector<CbmTrdTrack*>::iterator ivTempTrdTrackCand;
  //  vector<CbmTrdTrack*> trackCand1;

  //  vector<CbmL1Track> vtrackcandidate;

  map<Int_t, Int_t> globalMapUsedHits;
  globalMapUsedHits.clear();

  //moved from the inside of the if structure
  //vector<CbmTrdTrack*> vTrdTrackTrue;
  //  vector<CbmTrdTrack*>::iterator ivTrdTrackCand;

  //#########################################################################################################

  //for(Int_t i = 0; i < 3; i++) {
  //sort(clTracklets[0].begin(), clTracklets[0].end(), CbmL1TrdTracklet4::Compare3S);
  //  sort(clTracklets[1].begin(), clTracklets[1].end(), CbmL1TrdTracklet4::Compare3S);
  //}

  //   vector<CbmL1TrdTracklet4*>::iterator itsort;
  //   for(itsort =  clTracklets[1].begin();
  //       itsort != clTracklets[1].end();
  //       itsort++) {
  //     cout << (*itsort)->GetVal() << endl;
  //   }


  cout 
    << "--------------------------------------------------" << endl
    << "   ### Starting creating the tracks               " << endl
    << "--------------------------------------------------" << endl;
  
  createTracks.Start();


  if(segment == 0) { //segment == 0 -> we create long tracks
        
    if(fVerbose > 2) {
      cout << "Outer area: " << endl
	   << "--- size of fArrayTrdTrack    = " << fArrayTrdTrack->GetEntriesFast() << endl
	   << "--- size of globalMapUsedHits = " << globalMapUsedHits.size() << endl;
    }

    //create long tracks from 3 sorts of tracklets - station 1,2,3
    CreateTracks(clTracklets[0], clTracklets[1], 
		 clTracklets[2],
		 globalMapUsedHits, 
		 removeUsedHits, competition,
		 fArrayTrdTrack);
    
    createTracks.Stop();
    
    //clearing and removing objects that are not needed anymore
    for(Int_t i = 0; i < 3; i++){
      DeleteTracklets(clTracklets[i]);
      clTracklets[i].clear();
    }
    //----------------------------------------------------

    for(Int_t i = 0; i < 6; i++){
      DeleteTracklets(clSpacePoints[i]); 
      clSpacePoints[i].clear();
    }


    //--- [Second Loop] ------------------------------------------------------
    secondLoopTime.Start();
    if(secondLoop) {    
      
      for(Int_t i = 0; i < 12; i++) {
	fvTrdHitArr[i].clear();
      }

      map<Int_t, Int_t>::iterator iglobalMapUsedHits;

      for (Int_t iHit = 0; iHit<nTrdHit; iHit++) {
	pHit = (CbmTrdHit*) hitArray->At(iHit);
	if ( ! pHit ) {
	  cout << "-E- CbmL1CATrdTrackFinderSA::DoFind: Empty slot "
	       << "in HitArray at position " << iHit << endl;
	  nNoTrdHit++;
	  continue;
	}
	
	iglobalMapUsedHits = globalMapUsedHits.find(iHit);
	if(iglobalMapUsedHits != globalMapUsedHits.end()) continue; 

	trdPlaneID = pHit->GetPlaneID();
	Int_t trdPlaneIDN = trdPlaneID-1;

	//planeHits.mcTrackID = pMCpt->GetTrackID();
	planeHits.hitIndex = iHit;
	planeHits.X  = pHit->GetX();
	planeHits.Y  = pHit->GetY();
	planeHits.Z  = pHit->GetZ();;
	planeHits.DX = pHit->GetDx()/1E+4;;
	planeHits.DY = pHit->GetDy()/1E+4;;
	fvTrdHitArr[trdPlaneIDN].push_back(planeHits);
      }

      for(Int_t i = 0; i < 12; i++){
	sort(fvTrdHitArr[i].begin(),  fvTrdHitArr[i].end(),  CompareY);
      }

      //       for(Int_t i = 0; i < 12; i++){
      // 	cout <<" Size of vector "<< i <<": "<< fvTrdHitArr[i].size() << endl;
      //       }

      //--------------------------------------------------------------------------------------


      createSPs_SL.Start();
      //performing second step - make all previous tasks but with less restrictive conditions
      cout << "[Second Loop]::Creating space points"<< endl;

      for(Int_t i = 0, j = 0; i < 6; i++, j=j+2) {
	CreateSpacePoints(fvTrdHitArr[j], fvTrdHitArr[j+1], 
			  clSpacePoints[i], sigmaA_SL, sigmaB_SL);
      }
      createSPs_SL.Stop();

      for(Int_t i = 0; i < 6; i++){
	sort(clSpacePoints[i].begin(), clSpacePoints[i].end(), CbmL1TrdTracklet::Compare1);
      }

      cout << "[Second Loop]::Creating tracklets" << endl;
      for(Int_t i = 0, j = 0; i < 3; i++, j=j+2) {
	CreateSegments(clSpacePoints[j], clSpacePoints[j+1], 
		       clTracklets[i], dX_SL, dY_SL);
      }


      cout << endl << "-I- CbmL1CATrdTrackFinderSA::DoFind (Second loop)" << endl
	   << "--------------------------------------------------" << endl
	   << "     ### (SL) in Event No. " << fNofEvents << " ###" <<endl
	   << "--------------------------------------------------" << endl
	   << "--- size of segment vector 14 = "
	   << clTracklets[0].size() 
	   << endl
	   << "--- size of segment vector 58 = "
	   << clTracklets[1].size() 
	   << endl
	   << "--- size of segment vector 912 = "
	   << clTracklets[2].size() 
	   << endl;
    
      cout << "[Second Loop]::Finding friends 14-58" << endl;
      FindNeighbour(clTracklets[0],clTracklets[1],
		    distPropLongX_SL, distPropLongY_SL);

      cout << "[Second Loop]::Finding friends 58-912" << endl;
      FindNeighbour(clTracklets[1],clTracklets[2],
		    distPropLongX_SL, distPropLongY_SL);

      cout << "[Second Loop]::Tagging segments 58-912" << endl;
      TagSegments(clTracklets[1], clTracklets[2]);

      //----------------------------------------------------------------------------------------------
      cout << "[Second Loop]::Tagging segments 14-58" << endl;
      TagSegments(clTracklets[0], clTracklets[1]);

      //----------------------------------------------------------------------------------------------

      cout << "[Second Loop]::Creating tracks" << endl;
      CreateTracks(clTracklets[0], clTracklets[1], 
		   clTracklets[2],
		   globalMapUsedHits, 
		   removeUsedHits, competition,
		   fArrayTrdTrack);
      //competition -> selects the best track candidate from each branch    
    }
    secondLoopTime.Stop();

  } //end segment == 0
  //####################################################################################################



  if(segment == 1) { //segment == 1 -> we create only 4-hit, 1-segment tracks!
    
    CreateAndManageSegments(clTracklets[0], clTracklets[1], 
			    clTracklets[2], fArrayTrdTrack);    
  } //koniec segment == 1

  cout << endl << "-I- CbmL1CATrdTrackFinderSA::DoFind " << endl
       << "--------------------------------------------------" << endl
       << "     ### Summary" << endl
       << "--------------------------------------------------" << endl
       << "--- Number of found tracks: " << fArrayTrdTrack->GetEntriesFast() << endl
       << endl;

  delTime.Start();
  
  cout << "\n-I- CbmL1CATrdTrackFinderSA::DoFind " << endl
       << "--------------------------------------------------" << endl
       << "     ### Delete:Clear" <<endl
       << "--------------------------------------------------" << endl;

  for(Int_t i = 0; i < 12; i++){
    fvTrdHitArr[i].clear();
 }

  usedHits9.clear();
  usedHits12.clear();

  //----------------------------------------------------
  for(Int_t i = 0; i < 3; i++){
    DeleteTracklets(clTracklets[i]);
    clTracklets[i].clear();
  }

  //----------------------------------------------------
  for(Int_t i = 0; i < 6; i++){
    DeleteTracklets(clSpacePoints[i]); 
    clSpacePoints[i].clear();
  }
  //----------------------------------------------------
  
  usedHits9.clear();
  usedHits12.clear();
  segValues14.clear();
  segValues58.clear();
  segValues912.clear();

  delTime.Stop();
  doFind.Stop();

  //-----------------------------------------------------------------
  Double_t rtime;
  Double_t ctime;
  Double_t qtime = Double_t(fNofEvents);
  
  cout << endl <<"\n-I- CbmL1CATrdTrackFinderSA::DoFind " << endl
       << "--------------------------------------------------" << endl
       << "     ### Time" <<endl
       << "--------------------------------------------------" << endl;
  cout << "         Do find: ";

  rtime = doFind.RealTime();
  ctime = doFind.CpuTime();
  totDoFind += ctime;
  printf("RealTime=%4.2f s, CpuTime=%4.2f s, TotCpu=%4.2f s, %4.2f/event\n",rtime,ctime,totDoFind,totDoFind/qtime);
  
  cout << "       Sort Hits: ";
  rtime = sortHits.RealTime();
  ctime = sortHits.CpuTime();
  totSortHits += ctime;
  printf("RealTime=%4.2f s, CpuTime=%4.2f s, TotCpu=%4.2f s, %4.2f/event\n",rtime,ctime,totSortHits,totSortHits/qtime);
  
  cout << "      Create SPs: ";
  rtime = createSPs.RealTime();
  ctime = createSPs.CpuTime();
  totCreateSPs += ctime;
  printf("RealTime=%4.2f s, CpuTime=%4.2f s, TotCpu=%4.2f s, %4.2f/event\n",rtime,ctime,totCreateSPs,totCreateSPs/qtime);
  
  cout << "        Sort SPs: ";
  rtime = sortSPs.RealTime();
  ctime = sortSPs.CpuTime();
  totSortSPs += ctime;
  printf("RealTime=%4.2f s, CpuTime=%4.2f s, TotCpu=%4.2f s, %4.2f/event\n",rtime,ctime,totSortSPs,totSortSPs/qtime);

  cout << "Create tracklets: ";
  rtime = createSegments.RealTime();
  ctime = createSegments.CpuTime();
  totCreateSegments += ctime;
  printf("RealTime=%4.2f s, CpuTime=%4.2f s, TotCpu=%4.2f s, %4.2f/event\n",rtime,ctime,totCreateSegments,totCreateSegments/qtime);
  
  cout << "     Find friend: ";
  rtime = findNeighbour.RealTime();
  ctime = findNeighbour.CpuTime();
  totFindNeighbour += ctime;
  printf("RealTime=%4.2f s, CpuTime=%4.2f s, TotCpu=%4.2f s, %4.2f/event\n",rtime,ctime,totFindNeighbour,totFindNeighbour/qtime);

  cout << "   Tag Tracklets: ";
  rtime = tagTracks.RealTime();
  ctime = tagTracks.CpuTime();
  totTagTracks += ctime;
  printf("RealTime=%4.2f s, CpuTime=%4.2f s, TotCpu=%4.2f s, %4.2f/event\n",rtime,ctime,totTagTracks,totTagTracks/qtime);
  
  cout << "   Create Tracks: ";
  rtime = createTracks.RealTime();
  ctime = createTracks.CpuTime();
  totCreateTracks += ctime;
  printf("RealTime=%4.2f s, CpuTime=%4.2f s, TotCpu=%4.2f s, %4.2f/event\n",rtime,ctime,totCreateTracks,totCreateTracks/qtime);

  cout << "     Second Loop: ";
  rtime = secondLoopTime.RealTime();
  ctime = secondLoopTime.CpuTime();
  totSecondLoopTime += ctime;
  printf("RealTime=%4.2f s, CpuTime=%4.2f s, TotCpu=%4.2f s, %4.2f/event\n",rtime,ctime,totSecondLoopTime,totSecondLoopTime/qtime);

  cout << "Deleting Objects: ";
  rtime = delTime.RealTime();
  ctime = delTime.CpuTime();
  totDelTime += ctime;
  printf("RealTime=%4.2f s, CpuTime=%4.2f s, TotCpu=%4.2f s, %4.2f/event\n",rtime,ctime,totDelTime,totDelTime/qtime);

  cout<<"[Second Loop]"<< endl;
  cout << "      Create SPs: ";
  rtime = createSPs_SL.RealTime();
  ctime = createSPs_SL.CpuTime();
  totCreateSPs_SL += ctime;
  printf("RealTime=%4.2f s, CpuTime=%4.2f s, TotCpu=%4.2f s, %4.2f/event\n",rtime,ctime,totCreateSPs_SL,totCreateSPs_SL/qtime);

  cout << endl << "-I- CbmL1CATrdTrackFinderSA::DoFind " << endl
       << "--------------------------------------------------" << endl
       << "     ### END: Track constructing phase completed" <<endl
       << "--------------------------------------------------" << endl;



  fNoEvTime->Fill(fArrayTrdTrack->GetEntriesFast(),doFind.RealTime());


  return 1;
}
// -----------------------------------------------------------------------

void  CbmL1CATrdTrackFinderSA::DeleteTracklets(vector<CbmL1TrdTracklet4*> vect)
{
  vector<CbmL1TrdTracklet4*>::iterator it;
  for(it  = vect.begin();
      it != vect.end();
      it++) {
    delete (*it);
  }
}

// -----------------------------------------------------------------------

void  CbmL1CATrdTrackFinderSA::DeleteTracklets(vector<CbmL1TrdTracklet*> vect)
{
  vector<CbmL1TrdTracklet*>::iterator it;
  for(it  = vect.begin();
      it != vect.end();
      it++) {
    delete (*it);
  }
}


// ---------------------- Create histogramms -----------------------------
void CbmL1CATrdTrackFinderSA::CreateHistogramms()
{
  // Normalized distance to hit
  fh_chi2hit = new TH1F("h_chi2hit",
			"Normalized distance to hit",
			500, 0., 50.);
  fh_chi2hit_plane = new TH2F("h_chi2hit_plane",
			      "Normalized distance to hit vs. plane number",
			      500, 0., 50.,
			      12, 0., 12.);

  fDistLongX   = new TH1F("DistLongX",  "DistLongX",  100,-200,200);
  fDistLongY   = new TH1F("DistLongY",  "DistLongY",  100,-200,200);
  fDistShortX  = new TH1F("DistShortX", "DistShortX", 100,-10,10);
  fDistShortY  = new TH1F("DistShortY", "DistShortY", 100,-10,10);

  fDistLongBX   = new TH1F("DistLongBX",  "DistLongBX",  100,-200,200);
  fDistLongBY   = new TH1F("DistLongBY",  "DistLongBY",  100,-200,200);
  fDistShortBX  = new TH1F("DistShortBX", "DistShortBX", 100,-200,200);
  fDistShortBY  = new TH1F("DistShortBY", "DistShortBY", 100,-200,200);

  fDistY12 = new TH1F("DistY12","DistY12",100,-1,1);

  fMomDistLongPrimaryX = new TH2F("MomDistLongPrimaryX",  "MomDistLongPrimaryX",  100,0,10,100,-30,30);
  fMomDistLongPrimaryY = new TH2F("MomDistLongPrimaryY",  "MomDistLongPrimaryY",  100,0,10,100,-30,30);
  fMomDistExtrapolPrimaryX = new TH2F("MomDistExtrapolPrimaryX",  "MomDistExtrapolPrimaryX",  100,0,10,200,-20,20);
  fMomDistExtrapolPrimaryY = new TH2F("MomDistExtrapolPrimaryY",  "MomDistExtrapolPrimaryY",  100,0,10,200,-20,20);

  fMomDistLongExtraX = new TH2F("MomDistLongExtraX",  "MomDistLongExtraX",  100,0,10,100,-30,30);
  fMomDistLongExtraY = new TH2F("MomDistLongExtraY",  "MomDistLongExtraY",  100,0,10,100,-30,30);
  fMomDistExtrapolExtraX = new TH2F("MomDistExtrapolExtraX",  "MomDistExtrapolExtraX",  100,0,10,200,-20,20);
  fMomDistExtrapolExtraY = new TH2F("MomDistExtrapolExtraY",  "MomDistExtrapolExtraY",  100,0,10,200,-20,20);

  fMomDistShortPrimaryX = new TH2F("MomDistShortPrimaryX",  "MomDistShortPrimaryX",  100,0,10,100,-5,5);
  fMomDistShortPrimaryY = new TH2F("MomDistShortPrimaryY",  "MomDistShortPrimaryY",  100,0,10,100,-5,5);
  fMomDistShortExtraX   = new TH2F("MomDistShortExtraX",    "MomDistShortExtraX",    100,0,10,100,-5,5);
  fMomDistShortExtraY   = new TH2F("MomDistShortExtraY",    "MomDistShortExtraY",    100,0,10,100,-5,5);

  fDistY = new TH1F("DistY","DistY",1000,-10,10);
  fDistX = new TH1F("DistX","DistX",1000,-10,10);

  fPlane1Ydens = new TH1F("Plane1Ydens","Plane1Ydens",1000,-1000,1000);
  fPlane5Ydens = new TH1F("Plane5Ydens","Plane5Ydens",1000,-1000,1000);
  fPlane9Ydens = new TH1F("Plane9Ydens","Plane9Ydens",1000,-1000,1000);

  fSPlengthMC = new TH1F("SPlengthMC","SPlengthMC",200,0,20);
  fSPlength   = new TH1F("SPlength",  "SPlength",  200,0,20);

  fYat0   = new TH1F("Yat0",   "Yat0",   500,-500,500);
  fYat0MC = new TH1F("Yat0MC", "Yat0MC", 500,-500,500);

  fNoEvTime = new TH2F("NoEvTime","NoEvTime",1000,0,1000,1000,0,5);

}
// -----------------------------------------------------------------------



// -------------------- Activates data branches --------------------------
void CbmL1CATrdTrackFinderSA::DataBranches()
{
  // Get pointer to the ROOT manager
  CbmRootManager* rootMgr = CbmRootManager::Instance();
  if(NULL == rootMgr) {
    cout << "-E- CbmL1CATrdTrackFinderSA::DataBranches : "
	 << "ROOT manager is not instantiated" << endl;
    return;
  }

  fArrayTrdHit  = (TClonesArray*) rootMgr->GetObject("TRDHit");
  if ( ! fArrayTrdHit) {
    cout << "-E- CbmL1CATrdTrackFinderSA::Init: No TrdHit array!"
	 << endl;
    return;
  }
}
// -----------------------------------------------------------------------



// ------------------- Determines the TRD layout -------------------------
void CbmL1CATrdTrackFinderSA::TrdLayout()
{
  // Get the pointer to the singleton CbmRunAna object
  CbmRunAna* ana = CbmRunAna::Instance();
  if(NULL == ana) {
    cout << "-E- CbmL1CATrdTrackFinderSA::TrdLayout :"
	 <<" no CbmRunAna object!" << endl;
    return;
  }
  // Get the pointer to run-time data base
  CbmRuntimeDb* rtdb = ana->GetRuntimeDb();
  if(NULL == rtdb) {
    cout << "-E- CbmL1CATrdTrackFinderSA::TrdLayout :"
	 <<" no runtime database!" << endl;
    return;
  }
  // Get the pointer to container of base parameters
  CbmBaseParSet* baseParSet =
    (CbmBaseParSet*) rtdb->getContainer("CbmBaseParSet");
  if(NULL == baseParSet) {
    cout << "-E- CbmL1CATrdTrackFinderSA::TrdLayout :"
	 <<" no container of base parameters!" << endl;
    return;
  }

  
  TrdPar = (CbmGeoTrdPar*) (rtdb->findContainer("CbmGeoTrdPar"));
  TObjArray *Nodes = TrdPar->GetGeoPassiveNodes();
  for( Int_t i=0;i<Nodes->GetEntries(); i++) {
    CbmGeoNode *node = dynamic_cast<CbmGeoNode*> (Nodes->At(i));
    if ( !node ) continue;
    
    TString name = node->getName();
    //TString Sname = node->getShapePointer()->GetName();	
    CbmGeoVector nodeV = node->getLabTransform()->getTranslation(); //in cm

    //     Int_t id = node->getMCid();
    //     cout <<"name: "<< name <<"\tid: "<< id << endl;
    //     cout <<" name: "<< name
    // 	 <<" value: " << nodeV.Z() << endl;	  
    //     if(name.Contains("trd13")) fTrd13_Z = nodeV.Z();
    //     if(name.Contains("trd14")) fTrd14_Z = nodeV.Z();
    //     if(name.Contains("trd21")) fTrd21_Z = nodeV.Z();
    //     if(name.Contains("trd24")) fTrd24_Z = nodeV.Z();
    //     if(name.Contains("trd31")) fTrd31_Z = nodeV.Z();
    
    if(name == "trd13") fTrd13_Z = nodeV.Z();
    if(name == "trd14") fTrd14_Z = nodeV.Z();
    if(name == "trd21") fTrd21_Z = nodeV.Z();
    if(name == "trd24") fTrd24_Z = nodeV.Z();
    if(name == "trd31") fTrd31_Z = nodeV.Z();
  }

  /*
    TrdPar = (CbmGeoTrdPar*) (rtdb->findContainer("CbmGeoTrdPar"));
    //TObjArray *Nodes = TrdPar->GetGeoSensitiveNodes();
    Nodes = TrdPar->GetGeoSensitiveNodes();
    for( Int_t i=0;i<Nodes->GetEntries(); i++) {
    CbmGeoNode *node = dynamic_cast<CbmGeoNode*> (Nodes->At(i));
    if ( !node ) continue;
    
    TString name = node->getName();
    CbmGeoVector nodeV = node->getLabTransform()->getTranslation(); //in cm

    Int_t uid = i+1;
    cout <<"name, uid, Z: "<< name <<", "<< uid <<", "<< nodeV.Z() << endl;
    
    if(uid == 3) fTrd13_Z = nodeV.Z();
    if(uid == 4) fTrd14_Z = nodeV.Z();
    if(uid == 5) fTrd21_Z = nodeV.Z();
    if(uid == 8) fTrd24_Z = nodeV.Z();
    if(uid == 9) fTrd31_Z = nodeV.Z();
    }
  */


  // Get the pointer to detector list
  TObjArray *detList = baseParSet->GetDetList();
  if(NULL == detList) {
    cout << "-E- CbmL1CATrdTrackFinderSA::TrdLayout :"
	 << " no detector list!" << endl;
    return;
  }
  // Find TRD detector
  CbmDetector* trd = (CbmDetector*) detList->FindObject("TRD");
  if(NULL == trd) {
    cout << "-E- CbmL1CATrdTrackFinderSA::TrdLayout :"
	 << " no TRD detector!" << endl;
    return;
  }
  // Determine the geometry version
  TString name = trd->GetGeometryFileName();
  if(name.Contains("9")) {
    cout << "-I- CbmL1CATrdTrackFinderSA::TrdLayout :"
	 << " TRD geometry : 3x3." << endl;
    fNoTrdStations = 3;
    fNoTrdPerStation = 3;
  } else if(name.Contains("12")) {
    cout << "-I- CbmL1CATrdTrackFinderSA::TrdLayout :"
	 << " TRD geometry : 3x4." << endl;
    fNoTrdStations = 3;
    fNoTrdPerStation = 4;
  } else if(name.Contains("6x2")) {
    cout << "-I- CbmL1CATrdTrackFinderSA::TrdLayout :"
	 << " TRD geometry : 6x2." << endl;
    fNoTrdStations = 6;
    fNoTrdPerStation = 2;
  }
}
// -----------------------------------------------------------------------


// ---------------------------- CutTracklet -----------------------------------------------
Double_t CbmL1CATrdTrackFinderSA::DistTrackletCutY(TClonesArray* fHit,
						   Int_t iIndexFirst,   //index of the 1st hit
						   Int_t iIndexSecond)  //index of the 3rd hit
//calculates two points distance in the YZ midplane
//we take Y coordinates of the 1st and 3rd hits in the tracklet
{
  CbmTrdHit
    *pHitFirst,
    *pHitSecond;
  Double_t dist;

  Double_t firstPosY = 0;
  Double_t firstPosZ = 0;


  // Get hits from first and second plane in station
  if(iIndexFirst != -1) {
    pHitFirst  = (CbmTrdHit*) fHit->At(iIndexFirst);
    firstPosY  = fabs(pHitFirst->GetY());
    firstPosZ  = pHitFirst->GetZ();
  }
  
  
  pHitSecond = (CbmTrdHit*) fHit->At(iIndexSecond);

  // Get position Y & Z of hits

  Double_t secondPosY = fabs(pHitSecond->GetY());
  Double_t secondPosZ = pHitSecond->GetZ();

  Double_t Y = (firstPosZ*secondPosY)/secondPosZ;
  dist = fabs(Y-firstPosY);

  //   cout << "From the coords (y,z) = (" 
  //        << firstPosY << ", " << firstPosZ << endl
  //        << "From the coords (y,z) = (" 
  //        << secondPosY << ", " << secondPosZ << endl
  //        << "     we get the value = "
  //        << dist 
  //        << endl;

  return dist;
}
// -----------------------------------------------------------------------


// -----------------------------------------------------------------------
Double_t CbmL1CATrdTrackFinderSA::DistTrackletCutX(TClonesArray* trdHit,
						   Int_t iIndexFirst,    //index of the 2nd hit
						   Int_t iIndexSecond)   //index of the 4th hit
//calculates two points distance in the XZ midplane
//we take X coordinates of the 2nd and 4th hits in the tracklet
{
  CbmTrdHit
    *trdHitFirst,
    *trdHitSecond;
  Double_t dist;

  // Get hits from 2nd and 4th plane in station
  trdHitFirst  = (CbmTrdHit*) trdHit->At(iIndexFirst);
  trdHitSecond = (CbmTrdHit*) trdHit->At(iIndexSecond);

  // Get position X & Z of hits
  Double_t firstPosX  = fabs(trdHitFirst->GetX());
  Double_t firstPosZ  = trdHitFirst->GetZ();

  Double_t secondPosX = fabs(trdHitSecond->GetX());
  Double_t secondPosZ = trdHitSecond->GetZ();

  Double_t X = (firstPosZ*secondPosX)/secondPosZ;
  dist = fabs(X-firstPosX); //!!

  return dist;
}
// -----------------------------------------------------------------------



// ------------------------- Write histogramms ---------------------------
void CbmL1CATrdTrackFinderSA::WriteHistogramms()
{
  fh_chi2hit->Write();
  fh_chi2hit_plane->Write();

  fDistLongBX->Write();
  fDistLongBY->Write();
  fDistShortBX->Write();
  fDistShortBY->Write();

  fDistY12->Write();

  //----------------------------------------------------------
  fMomDistLongPrimaryX->Write();
  fMomDistLongPrimaryY->Write();
  fMomDistLongExtraX->Write();
  fMomDistLongExtraY->Write();

  fMomDistExtrapolPrimaryX->Write();
  fMomDistExtrapolPrimaryY->Write();
  fMomDistExtrapolExtraX->Write();
  fMomDistExtrapolExtraY->Write();

  fMomDistShortPrimaryX->Write();
  fMomDistShortPrimaryY->Write();
  fMomDistShortExtraX->Write();
  fMomDistShortExtraY->Write();

  fDistLongX->Write();
  fDistLongY->Write();
  fDistShortX->Write();
  fDistShortY->Write();

  fDistY->Write();
  fDistX->Write();

  fPlane1Ydens->Write();
  fPlane5Ydens->Write();
  fPlane9Ydens->Write();

  fSPlength->Write();
  fSPlengthMC->Write();

  fYat0->Write();
  fYat0MC->Write();

  fNoEvTime->Write();
}
// -----------------------------------------------------------------------


// -----------------------------------------------------------------------
Double_t CbmL1CATrdTrackFinderSA::FitLinear(CbmTrdTrack *tr, 
					    Int_t var = 1) {
  //fit using a least square method
  

  //  cout << ">>> Trying to fit a track..." << endl;
  Int_t noHits = tr->GetNofTrdHits();
  //  cout << "No Hits in the track : " << noHits << endl;
  Int_t iHit;
  CbmTrdHit *hit;
  
  Double_t C1[12];
  Double_t C2[12];
  Double_t  Z[12];

  Double_t n = 12;

  for(int i=0;i<12;i++) {
    C1[i] = 0;
    C2[i] = 0;
    Z [i] = 0;
  }

  //  Int_t ind = 0;

  for(int i=0;i<noHits;i++) {
    
    iHit = tr->GetTrdHitIndex(i);
    hit = (CbmTrdHit*)fArrayTrdHit->At(iHit);
    
    C1[i] = hit->GetX();
    C2[i] = hit->GetY();
    Z[i]  = hit->GetZ();
  }

  Double_t 
    sumXZ = 0, 
    sumYZ = 0, 
    sumX  = 0, 
    sumY  = 0, 
    sumZx = 0, 
    sumZy = 0, 
    sumX2 = 0,
    sumY2 = 0,
    sumZx2 = 0,
    sumZy2 = 0;  

  for(int i=0;i<12;i+=2) {
    //    cout << "C1[" << i << "]=" << C1[i] << endl; 
    //    cout << "C2[" << i << "]=" << C2[i] << endl; 
    
    //    if(!(i % 2)) 
    {
      sumXZ += C1[i]*Z[ i];
      
      sumX  += C1[i];
      sumX2 += C1[i]*C1[i];

      sumZx  += Z[i];
      sumZx2 += Z[i]*Z[ i];

      
    }
  }

  for(int i=1;i<12;i+=2) {
    //    cout << "C1[" << i << "]=" << C1[i] << endl; 
    //    cout << "C2[" << i << "]=" << C2[i] << endl; 
    
    //    if(!((1i % 2)) 
    {
      sumYZ += C2[i]*Z[ i];
      
      sumY  += C2[i];
      sumY2 += C2[i]*C2[i];
      
      sumZy  += Z[ i];
      sumZy2 += Z[ i]*Z[ i];
    }
  }


  cout << "";
  //   Double_t a = (n*sumXZ - sumX*sumZ)/(n*sumX2 - sumX*sumX);
  //   Double_t b = (sumZ - a*sumX)/n;
  Double_t r1 = (n*sumXZ - sumX*sumZx)/sqrt((n*sumX2-sumX*sumX)*(n*sumZx2-sumZx*sumZx));
  //  Double_t Sa = sqrt((n*(sumY2 -a*sumXZ-b*sumY))/((n-2)*(n*sumX2-sumX*sumX)));
  //  Double_t Sb = Sa*sqrt(sumX2/n);


  Double_t r2 = (n*sumYZ - sumY*sumZy)/sqrt((n*sumY2-sumY*sumY)*(n*sumZy2-sumZy*sumZy));


  //   cout << "Linear coefficients: a= " 
  //        << a << ", Sa= " << Sa << ", b= " << b << ", Sb= " << Sb << endl 
  //        << ", r=" << r << ", chi2=" << tr->GetChi2() << endl;
    

  return sqrt(r1*r1+r2*r2);
}
// -----------------------------------------------------------------------



// -----------------------------------------------------------------------
Double_t CbmL1CATrdTrackFinderSA::DistTwoTrackletsY(TClonesArray* fHit,
						    Int_t iIndexFirst,
						    Int_t iIndexSecond,
						    Double_t zed)
{
  //returns an extrapolated coordinate
  CbmTrdHit
    *pHitFirst,
    *pHitSecond;
  Double_t dist;

  Double_t Y1 = 0;
  Double_t Z1 = 0;

  // Get hits from first and second plane in station
  if(iIndexFirst != -1) {
    pHitFirst  = (CbmTrdHit*) fHit->At(iIndexFirst);
    Y1  = pHitFirst->GetY();
    Z1  = pHitFirst->GetZ();
  }
  
  pHitSecond = (CbmTrdHit*) fHit->At(iIndexSecond);

  // Get position Y & Z of hits

  Double_t Y4 = pHitSecond->GetY();
  Double_t Z4 = pHitSecond->GetZ();

  Double_t a = (Y4 - Y1)/(Z4 - Z1);
  Double_t b = Y1 - a*Z1;
  Double_t Y = a*zed + b ;// - Y1;
  dist = Y;
  //      cout << "From the coords (y1,z1) = (" 
  // 	  << Y1 << ", " << Z1 << endl
  //           << "From the coords (y2,z2) = (" 
  // 	  << Y4 << ", " << Z4 << endl
  //           << "     we get the value = "
  //           << dist 
  //           << endl;
  return dist;
}
// -----------------------------------------------------------------------



// -----------------------------------------------------------------------

Double_t CbmL1CATrdTrackFinderSA::DistTwoTrackletsX(TClonesArray* fHit,
						    Int_t iIndexFirst,
						    Int_t iIndexSecond,
						    Double_t zed)
{
  CbmTrdHit
    *pHitFirst,
    *pHitSecond;
  Double_t dist;

  Double_t X1 = 0;
  Double_t Z1 = 0;

  // Get hits from first and second plane in station
  if(iIndexFirst != -1) {
    pHitFirst  = (CbmTrdHit*) fHit->At(iIndexFirst);
    X1  = pHitFirst->GetX();
    Z1  = pHitFirst->GetZ();
  }
  
  pHitSecond = (CbmTrdHit*) fHit->At(iIndexSecond);


  Double_t X4 = pHitSecond->GetX();
  Double_t Z4 = pHitSecond->GetZ();

  Double_t a = (X4 - X1)/(Z4 - Z1);
  Double_t b = X1 - a*Z1;
  //  Double_t X = (dz)/(dx)*(zed - Z1);
  Double_t X = a*zed + b;
  //    = dy/dz * (zed + Z1) - X1;
  dist = X;
  //       cout << "From the coords (x1,z1) = (" 
  // 	   << X1 << ", " << Z1 << endl
  //            << "From the coords (x4,z4) = (" 
  // 	   << X4 << ", " << Z4 << endl
  // 	   << "and at zed = " << zed << endl
  //            << "     we get the value = "
  //            << dist 
  //            << endl;
  return dist;
}
// -----------------------------------------------------------------------

// -----------------------------------------------------------------------
Bool_t CbmL1CATrdTrackFinderSA::FindNeighbour(vector<CbmL1TrdTracklet4*> &v1,
					      vector<CbmL1TrdTracklet4*> &v2,
					      Double_t dY, Double_t dX) {

  vector<CbmL1TrdTracklet4*>::iterator iv1;
  CbmL1TrdTracklet4 *tr1;

  Double_t
    extY = 0,extX = 0,  //two extrapolated coords from the first tracklet
    mesY = 0,mesX = 0,  ////two measured coords from the second tracklet
    aminY, amaxY,
    aminX, amaxX;
  Int_t indexA, indexB;
  Int_t
    Left = 0,
    Right,
    Mid = 0;

  Int_t il = 0;
  for(iv1 = v1.begin(); iv1 != v1.end(); iv1++) {
    il++;
    tr1 = *iv1;
    extY = tr1->GetExt1();
    extX = tr1->GetExt2();
    indexA = tr1->GetIndex();
    amaxY = extY + dY;
    aminY = extY - dY;
    amaxX = extX + dX;
    aminX = extX - dX;
    
    Left = 0;
    Right = v2.size();
    Mid = 0;

    // ----------- looking by bisection ------------------
    while(Left < Right){
      Mid = (Left+Right)/2;
      mesY = v2[Mid]->GetCoord1();

      if (amaxY < mesY){
	Left = Mid+1;	
      }else{
	Right = Mid-1;
      }
      /*
	cout << il << ": Size: "<< v1.size()
	<<" mesY: "<< v2[Mid]->GetCoord1()
	<<" extY: "<< extY
	<<" Mid: "<< Mid
	<< endl;
      */
    }
    //----------------------------------------------------
    
    Int_t size = v2.size();
    for(Int_t i = Mid; i < size; i++) {      
      mesY = v2[i]->GetCoord1();
      mesX = v2[i]->GetCoord2();
      indexB = v2[i]->GetIndex();
      
      if (mesY > amaxY) continue;
      if (mesY < aminY) break;
      
      if(aminX < mesX && mesX < amaxX) {
	tr1->vAccostRight.push_back(indexB);
	v2[i]->vAccostLeft.push_back(indexA);
      }
    }      
  }   
  return true;
}
// -----------------------------------------------------------------------



// -----------------------------------------------------------------------
Bool_t CbmL1CATrdTrackFinderSA::OverlapsHitsXY(Int_t posA,
					       Int_t posB)
{
  Bool_t overlap = false;
  
  CbmTrdHit
    *pHitA,
    *pHitB;
  
  Double_t
    hitPosA_X,
    hitPosA_Y,
    hitPosB_X,
    hitPosB_Y,
    dMul1 = 3, //sigma multiplier for calculating the range of overlapping hits
    dMul2 = 3; //sigma multiplier for calculating the range of overlapping hits

  Double_t
    hitPosA_DX,
    hitPosA_DY,
    hitPosB_DX,
    hitPosB_DY;
  
  pHitA = (CbmTrdHit*) fArrayTrdHit->At(posA);
  hitPosA_X  = pHitA->GetX();
  hitPosA_Y  = pHitA->GetY();
  hitPosA_DX = pHitA->GetDx()/1E+4;
  hitPosA_DY = pHitA->GetDy()/1E+4;

  pHitB = (CbmTrdHit*) fArrayTrdHit->At(posB);
  hitPosB_X  = pHitB->GetX();
  hitPosB_Y  = pHitB->GetY();
  hitPosB_DX = pHitB->GetDx()/1E+4;
  hitPosB_DY = pHitB->GetDy()/1E+4;

  //   cout 
  //     << "Overlap function...\n"
  //     << "X1=" << hitPosA_X << ", dX1=" << hitPosA_DX << ", Y1=" << hitPosA_Y << ", dY1=" << hitPosA_DY
  //     << "X2=" << hitPosB_X << ", dX2=" << hitPosB_DX << ", Y2=" << hitPosB_Y << ", dY2=" << hitPosB_DY
  //     << endl;

  
  if(((hitPosA_X+dMul1*hitPosA_DX) >= (hitPosB_X-dMul2*hitPosB_DX)) &&
     ((hitPosA_X-dMul1*hitPosA_DX) <= (hitPosB_X+dMul2*hitPosB_DX)) &&
     ((hitPosA_Y+dMul2*hitPosA_DY) >= (hitPosB_Y-dMul1*hitPosB_DY)) &&
     ((hitPosA_Y-dMul2*hitPosA_DY) <= (hitPosB_Y+dMul1*hitPosB_DY))){
    overlap = true;
  }
  
  //   CbmTrdPoint *pt1 = (CbmTrdPoint*)fMCPointArray->At(posA);
  //   CbmTrdPoint *pt2 = (CbmTrdPoint*)fMCPointArray->At(posB);
 
  //   if(overlap && (pt1->GetTrackID() == pt2->GetTrackID())) return true;
  //   else return false;

  return overlap;
}


// -----------------------------------------------------------------------
Bool_t CbmL1CATrdTrackFinderSA::TagSegments(vector <CbmL1TrdTracklet4*>& clTrackletsA,
					    vector <CbmL1TrdTracklet4*>& clTrackletsB) {
  //asigning numbers to each segment, 4- highest, 0 - lowest

  //  cout << "TagSegments: engaging... " << endl;

  //  vector<CbmL1TrdTracklet4*>::iterator itclTrackletsLeft;
  vector<CbmL1TrdTracklet4*>::iterator itclTrackletsRight;
  
  CbmL1TrdTracklet4* clSegRight;
  
  vector <Int_t> vLeft;
  vector <Int_t> vRight;

  vector <Int_t>::iterator ivLeft;
  //  vector <Int_t>::iterator ivRight;

  Int_t
    valLeft  = 0,
    valRight = 0;
    
  //     CbmL1TrdTracklet4* clSegLeft; //segment farther from the target
  //     CbmL1TrdTracklet4* clSegRight; //segment nearer to the target
  

  //  cout << "Engaging Right loop " << endl;
  for(itclTrackletsRight = clTrackletsB.begin();
      itclTrackletsRight != clTrackletsB.end();
      itclTrackletsRight++) {
      
    clSegRight = *itclTrackletsRight;
    //    cout << "---> X of tracklet = " << clSegRight->GetCoord1() << endl;
    vLeft = clSegRight->vAccostLeft;
    //    cout << "Size of Right vector = " << clSegRight->vAccostRight.size() << endl;
    //    cout << "Size of  Left vector = " << clSegRight->vAccostLeft.size() << endl;
    //    iHitRightA = clSegRight->GetInd1();
    //    iHitRightB = clSegRight->GetInd2();

    //      cout << "Tracklets 912: iHitRight:   " << iHitRight << "\t" << endl;
    if(vLeft.size() > 0) {
      //      cout << "Processing Comrade vector of size " << vLeft.size() << endl;
      valRight = clSegRight->GetVal();
      if(valRight == 0) {
	clSegRight->SetVal(1);
	valRight++;
      }
      for(ivLeft = vLeft.begin(); ivLeft != vLeft.end(); ivLeft++) {
	
	valLeft  = (clTrackletsA[*ivLeft])->GetVal();
	if(valLeft <= valRight) (clTrackletsA[*ivLeft])->SetVal(valRight+1);	
      }
    }
    //    cout << "clSegRight = " << clSegRight->GetVal() << endl;
    //    if(isAlone == true) clSegRight->SetIsAlone(true);
  }
  
  return true;
}
// -----------------------------------------------------------------------



	
	    

// -----------------------------------------------------------------------
Bool_t CbmL1CATrdTrackFinderSA::CreateTracks(vector <CbmL1TrdTracklet4*> clTracklets14,
					     vector <CbmL1TrdTracklet4*> clTracklets58,
					     vector <CbmL1TrdTracklet4*> clTracklets912,
					     map<Int_t,Int_t> &globalMapUsedHits,
					     Bool_t removeUsedHits,
					     Bool_t competition,
					     TClonesArray *fArrayTrdTrack) {
  //create long tracks from previously selected segments 

  if(fVerbose > 2) {
    cout << "Inner area: " << endl
	 << "--- size of fArrayTrdTrack    = " << fArrayTrdTrack->GetEntriesFast() << endl
	 << "--- size of globalMapUsedHits = " << globalMapUsedHits.size() << endl;
  }

  CbmL1TrdTracklet4 
    *clTrdSeg2,
    *clTrdSeg1;
  
  vector<CbmL1TrdTracklet4*>::iterator itclSeg3;

  Bool_t isEmpty = true;

  Int_t
    licz = 0,
    counter = 0;
   
  CbmTrdTrack* tr1 = NULL;
  //CbmTrdTrack* tr2 = NULL;
  CbmTrdHit* hit1 = NULL;

  vector<CbmTrdTrack*> vTempTrdTrackCand;
  vector<CbmTrdTrack*>::iterator ivTempTrdTrackCand;

  vector<CbmTrdTrack*> toRemove;  
  //  vector<CbmTrdTrack*>::iterator itoRemove;  

  map<Int_t, Int_t>::iterator iglobalMapUsedHits;
  vector<CbmTrdTrack*>::iterator ivTrdTrackCand;
  
  Int_t iFakeTrack = 0;    
  Int_t iHit = 0;
  Bool_t bTrueTrack = true;

  CbmTrackParam *trParam;
  trParam = new CbmTrackParam();
  trParam->SetQp(1.);

  //CbmTrdTrack *tr[10];
  /** almost: only make sens more than 1 if noHitsAccepted > 0 **/
  Int_t trMax = 1;          //number of tracks from one tree accepted as tracks candidates
  //  Double_t chi2A = 1E+7;    //maximum chi2 value
  //Double_t chi2B = 0;       //minimum chi2 value
  Int_t noHitsAccepted = 0; //number of used hits above which the tracks is recognized as fake; 0 = all tracks are fake.

  vector<TempTrackStruct> auxTrack;
  Int_t trdInd;
  vector <TempTrackStruct>::iterator ikol;

  cout << "--- Engaging main loop..." << endl;

  //   cout << "- No of outer iterations to go: " << clTracklets14.size() << endl;
  //   Int_t vSize = Int_t(clTracklets14.size());
  //   for(int i = 0; i < vSize; i++) {
  //     if( (clTracklets14[i])->GetVal() == 3) licz++;
  //   }
  //   cout << "* No of 3 in vector: " << licz << endl;
  licz = 0;

  Int_t trackArrayInd = fArrayTrdTrack->GetEntriesFast();
  //  cout <<"Array of global tracksCT: " << fArrayTrdTrack->GetEntries() << endl;;

  for(itclSeg3 = clTracklets14.begin();
      itclSeg3 != clTracklets14.end();
      itclSeg3++) { //loop over tracklets with value of 3
    if((*itclSeg3)->GetVal() == 3) {
      
      Int_t iIn2 = 0;
      tempTrack.M[0] = (*itclSeg3)->GetInd1();
      tempTrack.M[1] = (*itclSeg3)->GetInd2();
      tempTrack.M[2] = (*itclSeg3)->GetInd3();
      tempTrack.M[3] = (*itclSeg3)->GetInd4();

      Int_t size2 = Int_t((*itclSeg3)->vAccostRight.size());
      for(Int_t iSeg2 = 0; iSeg2 < size2; iSeg2++) {//loop over clTracklets with value of 2	  
	iIn2 = (*itclSeg3)->vAccostRight[iSeg2];
	// 	  cout << "We have a 2 value, vector index= " 
	// 	       << iIn2 << endl;

	clTrdSeg2 = clTracklets58[iIn2];
	tempTrack.M[4] = clTrdSeg2->GetInd1();
	tempTrack.M[5] = clTrdSeg2->GetInd2();
	tempTrack.M[6] = clTrdSeg2->GetInd3();
	tempTrack.M[7] = clTrdSeg2->GetInd4();

	Int_t iIn1 = 0;

	Int_t size1 = Int_t(clTrdSeg2->vAccostRight.size());
	for(Int_t iSeg1 = 0; iSeg1 < size1; iSeg1++) {//loop over clTracklets with value of 1
	  iIn1 = clTrdSeg2->vAccostRight[iSeg1];

	  clTrdSeg1 = clTracklets912[iIn1];

	  tempTrack.M[8]  = clTrdSeg1->GetInd1();
	  tempTrack.M[9]  = clTrdSeg1->GetInd2();
	  tempTrack.M[10] = clTrdSeg1->GetInd3();
	  tempTrack.M[11] = clTrdSeg1->GetInd4();

	  counter++;
	  isEmpty = false;

// 	  //-------------------------------------	  
// 	  tr2 = new CbmTrdTrack();
// 	  tr2->SetParamFirst(*trParam);
// 	  tr2->SetParamLast(*trParam);
	  
// 	  for(Int_t we = 0; we < 12; we++) {
// 	    hit1 = (CbmTrdHit*)fArrayTrdHit->At(tempTrack.M[we]);
// 	    tr2->AddHit(tempTrack.M[we], hit1);
// 	  }
// 	  tr2->SortHits();
// 	  trdTrackFitter->DoFit(tr2);
// 	  tempTrack.Chi2 = tr2->GetChi2();
// 	  delete tr2;
// 	  //-------------------------------------	  
  
// 	  tr1->SetParamFirst(*trParam);
// 	  tr1->SetParamLast(*trParam);
	    
	  //	    cout << ">>> " << tr1->GetNofTrdHits() << " hits in this track." << endl;
	  //	    tr1->Print();    
	  //     	trParam2 = tr1->GetParamLast();
	  //		cout << "Q/p= " << trParam2->GetQp() <<", tx=" << trParam2->GetTx() << trParam2->GetTy() << endl;
	  //delete trParam;	 	  

	  //chi2B = FitLinear(tempTrack.M, 1);	  
	  //chi2B = FitTLinearFitter(tempTrack.M);
	  tempTrack.Chi2 = Fit(tempTrack.M);;

	  //chi2B = Fit(tr1);
	  //tr1->SetChi2(chi2B);

	  //trdTrackFitter->DoFit(tr1);
	  //chi2B = tr1->GetChi2();

	  //toRemove.push_back(tr1);

	  auxTrack.push_back(tempTrack);

	  // Selection the track with lowest chi2 
	  // 	  if(chi2B < chi2A){
	  // 	    chi2A = chi2B;
	      
	  // 	    Int_t k = 0;
	  // 	    for(Int_t i = 0; i < trMax; i++){
	  // 	      k = i + 1;
	  // 	      tr[k] = tr[i];
	  // 	      tr[i] = tr1;
	  // 	    }
	  // 	  }
	  //---------------------------------------
	      
	  if((counter%1000) == 0) {
	    printf("\rTracks: %i",counter);
	    fflush(stdout);
	  }	      	      
	}//end::loop 1
      }//end::loop 2

      //if(tr2 == NULL) cout << "*** WARNING: tr2 is NULL" << endl;

      //       for(Int_t i = 0; i < trMax; i++){
      // 	if(tr[i] != NULL) vTempTrdTrackCand.push_back(tr[i]);		    
      //       }

      sort(auxTrack.begin(), auxTrack.end(), CompareChi2);
      
      Int_t li = 0;
      for(ikol = auxTrack.begin();
	  ikol != auxTrack.end();
	  ikol++) {
	//cout << (*ikol).Chi2 << endl;
	if(li >= trMax) break;
	
	tr1 = new CbmTrdTrack();

	for(Int_t we = 0; we < 12; we++) {
	  trdInd = (*ikol).M[we];
	  hit1 = (CbmTrdHit*)fArrayTrdHit->At((*ikol).M[we]);
	  tr1->AddHit(trdInd, hit1);
	}	    
	tr1->SortHits();

	tr1->SetParamFirst(*trParam);
	tr1->SetParamLast(*trParam);

	vTempTrdTrackCand.push_back(tr1);

	li++;
      }      
      auxTrack.clear();

      //cout <<"---------------------------------------"<<endl;
    }// loop 3
  }//loop over 3s

  if(competition) {
    cout << "\n--- Refiting by KF...";
    //Refit using KF fitter ---------------------------------------------
    //Int_t iKlops = 0;
    //    cout << "vTempTrdTrackCand.size(): " << vTempTrdTrackCand.size() << endl;
    for(ivTempTrdTrackCand = vTempTrdTrackCand.begin();
     	ivTempTrdTrackCand != vTempTrdTrackCand.end();
     	ivTempTrdTrackCand++) {
      //      iKlops++;
      trdTrackFitter->DoFit(*ivTempTrdTrackCand);
      //cout << iKlops <<": GetNofTrdHits(): " << (*ivTempTrdTrackCand)->GetNofTrdHits() << endl;
    }
    //-------------------------------------------------------------------
    cout << " [OK]" << endl;

    sort(vTempTrdTrackCand.begin(),vTempTrdTrackCand.end(),CbmTrdTrack::CompareChi2);
    
    //     //Erase if chi2 > 20 ------------------------------------------------
    //     CbmTrdTrack *t1;
    //for(ivTempTrdTrackCand = vTempTrdTrackCand.begin();
    // 	ivTempTrdTrackCand != vTempTrdTrackCand.end();
    // 	ivTempTrdTrackCand++) {
    //       t1 = (*ivTempTrdTrackCand);
    //       if(t1->GetChi2() >= 20) {
    // 	vTempTrdTrackCand.erase(ivTempTrdTrackCand, vTempTrdTrackCand.end());
    // 	break;
    //       }
    //     }
    //     //-------------------------------------------------------------------
  }
  
  CbmTrdTrack *trCand;
  //cout << "vTempTrdTrackCand.size(): " << vTempTrdTrackCand.size() << endl;

  Int_t firstHitSunk = 0;
  for(ivTrdTrackCand = vTempTrdTrackCand.begin();
      ivTrdTrackCand != vTempTrdTrackCand.end();
      ivTrdTrackCand++) { //Loop over track candidates; mark used hits
    trCand = (*ivTrdTrackCand);
    //cout << "Track no " << trackInd++ << " has " << (*ivTrdTrackCand)->GetNofTrdHits() << " hits and chi2= " << (*ivTrdTrackCand)->GetChi2() << endl;    

    if(removeUsedHits) {//RemoveHits
      Int_t noHitsA = trCand->GetNofTrdHits();

      bTrueTrack = true; 
      firstHitSunk = 0;

      for(int i = 0; i < noHitsA; i++) {
	//loop over hits
	//cout << "^ Getting index of hit no " << i << endl;
	//iHit = (*ivTrdTrackCand)->GetTrdHitIndex(i);
	iHit = trCand->GetTrdHitIndex(i);
	iglobalMapUsedHits = globalMapUsedHits.find(iHit);
	
	if(iglobalMapUsedHits != globalMapUsedHits.end()) {
	  if(firstHitSunk == noHitsAccepted) {
	    bTrueTrack = false;
	    iFakeTrack++;
	    break;
	  }
	  firstHitSunk++;
	}
      }

      if(bTrueTrack) {
	//cout << "^ We have a true track, adding a track to an array " << endl;
	new ((*fArrayTrdTrack)[trackArrayInd]) CbmTrdTrack(*trCand);
	trackArrayInd++;

	//cout << "^ Track added, size of array: " << trackArray->GetEntriesFast() << ", chi2=" <<tr2->GetChi2()<< endl;
	for(int i = 0; i < noHitsA; i++){
	  //set used hits in correct track
	  iHit = trCand->GetTrdHitIndex(i);
	  globalMapUsedHits.insert(pair<Int_t, Int_t>(iHit, 9));
	}
      }
      delete trCand;

    }else{ // end of RemoveHits	        
      new ((*fArrayTrdTrack)[trackArrayInd]) CbmTrdTrack(*trCand);
      trackArrayInd++;
      delete trCand;
    }
  }
  cout << "Track candidates: " << vTempTrdTrackCand.size() << endl
       << "Fake tracks:      " << iFakeTrack << endl;

  vTempTrdTrackCand.clear();
  
  if(!isEmpty) { licz++; isEmpty = true; }

  //cout << "toRemove.size():" << toRemove.size() << endl;
  //   for(itoRemove =  toRemove.begin();
  //       itoRemove != toRemove.end();
  //       itoRemove++) {
  //     delete (*itoRemove);    
  //   }
  //   toRemove.clear();
 
  // for(itoRemove =  vTempTrdTrackCand.begin();
  //       itoRemove != vTempTrdTrackCand.end();
  //       itoRemove++) {
  //     delete (*itoRemove);    
  //   }
  //  vTempTrdTrackCand.clear();

  delete trParam;

  return true;
}
// -----------------------------------------------------------------------






// -----------------------------------------------------------------------
Bool_t CbmL1CATrdTrackFinderSA::CreateSegments(vector<CbmL1TrdTracklet*> clSpacePointsAB,
					       vector<CbmL1TrdTracklet*> clSpacePointsCD,
					       vector<CbmL1TrdTracklet4*>& clTrackletsAD,
					       Double_t dX, Double_t dY) {
  //--- method to create 4-hit tracklets from Space Points (smaller 2-hit tracklets) ---
  vector<CbmL1TrdTracklet*>::iterator iSpacePtA;
  vector<CbmL1TrdTracklet*>::iterator iSpacePtB;
  vector<CbmL1TrdTracklet*>::iterator iSpacePtStart;

  Int_t
    iIndex = 0,       //index
    indSpacePtA,       //index of 1st miniTracklet
    indSpacePtB,       //index of 2nd miniTracklet
    indSpacePtC,       //index of 3rd miniTracklet
    indSpacePtD,       //index of 4th miniTracklet
    noXSurv = 0,      //number of tracklets that survived the X distance cut
    noYSurv = 0,      //number of tracklets that survived the Y distance cut
    noXYSurv = 0,     //number of tracklets that survived both X & Y distance cuts
    noAllPairs = 0,
    iSegAcc14 = 0;

  CbmL1TrdTracklet4 *clTr;

  CbmL1TrdTracklet
    *trLeft,
    *trRight;
  
  iSpacePtStart  = clSpacePointsCD.begin();

  //---------- BEGIN: Creation tracklet 12-34 ---------------------------------
  for(iSpacePtA  = clSpacePointsAB.begin();
      iSpacePtA != clSpacePointsAB.end();
      iSpacePtA++) {
    trLeft = *iSpacePtA;
    indSpacePtA = trLeft->GetIndLeft();
    indSpacePtB = trLeft->GetIndRight();

    //    cout << "GetCoord1: " << trLeft->GetCoord1() << endl;
    //cout << "--------------------" << clSpacePointsCD.size() << endl;


    //--- extrapolate --------------------------------------------
    Double_t
      y1  = trLeft->GetCoord1(),
      x1  = trLeft->GetCoord2(),
      y1z = trLeft->GetCoord3(),
      x1z = trLeft->GetCoord4(),
      x2  = 0,
      y2  = 0,
      distBetweenLayer = x1z-y1z,    
      y2z = y1z+distBetweenLayer*2,
      x2z = x1z+distBetweenLayer*2;

    y2 = (y1/y1z)*y2z;
    x2 = (x1/x1z)*x2z;
    //--- end: extrapolate --------------------------------------

    Int_t first = 1;

    //for(y = 0; y < trsize; y++){//y
    //trRight = clSpacePointsCD[y];
    //cout << trRight->GetCoord1() << endl;

    for(iSpacePtB  = iSpacePtStart;
	iSpacePtB != clSpacePointsCD.end();
	iSpacePtB++) {
      trRight = *iSpacePtB;

      indSpacePtC = trRight->GetIndLeft();
      indSpacePtD = trRight->GetIndRight();
      noAllPairs++;

      Double_t
	yB  = trRight->GetCoord1(),
	xB  = trRight->GetCoord2(),
	yBz = trRight->GetCoord3(),
	xBz = trRight->GetCoord4();
      
      /*
      //--- for estimation dx and dy ----------------------------
      trdHit1 = (CbmTrdHit*)fArrayTrdHit->At(indSpacePtA);
      Int_t refInd1 = trdHit1->GetRefIndex();
      trdHit1 = (CbmTrdHit*)fArrayTrdHit->At(indSpacePtB);
      Int_t refInd2 = trdHit1->GetRefIndex();
	  
      trdHit1 = (CbmTrdHit*)fArrayTrdHit->At(indSpacePtC);
      Int_t refInd3 = trdHit1->GetRefIndex();
      trdHit1 = (CbmTrdHit*)fArrayTrdHit->At(indSpacePtD);
      Int_t refInd4 = trdHit1->GetRefIndex();

      CbmTrdPoint
      *pPointAy,
      *pPointBx,
      *pPointCy,
      *pPointDx;
	  
      pPointAy = (CbmTrdPoint*) fMCPointArray->At(refInd1);
      pPointBx = (CbmTrdPoint*) fMCPointArray->At(refInd2);
      pPointCy = (CbmTrdPoint*) fMCPointArray->At(refInd3);
      pPointDx = (CbmTrdPoint*) fMCPointArray->At(refInd4);
      
      Int_t trIDAy = pPointAy->GetTrackID();
      Int_t trIDBx = pPointBx->GetTrackID();
      Int_t trIDCy = pPointCy->GetTrackID();
      Int_t trIDDx = pPointDx->GetTrackID();

      CbmMCTrack *pMCtrack;
      Int_t noTRDPoints;

      Double_t
      x1  = pPointBx->GetX(),
      y1  = pPointAy->GetY(),
      z1y = pPointAy->GetZ(),
      z1x = pPointBx->GetZ(),
      x2  = 0,
      y2  = 0,
      z2y = z1y+distBetweenLayer*2,
      z2x = z1x+distBetweenLayer*2;
      
      y2 = (y1/z1y)*z2y;
      x2 = (x1/z1x)*z2x;
      
      if(trIDAy == trIDCy){	
      pMCtrack = (CbmMCTrack*) fMCTrackArray->At(trIDAy);
      noTRDPoints = pMCtrack->GetTrdPoints();

      if(noTRDPoints == 12){
      Double_t fPosY = pPointCy->GetY();
      fDistY->Fill(fabs(fPosY-y2));
      }
      }

      if(trIDBx == trIDDx){	
      pMCtrack = (CbmMCTrack*) fMCTrackArray->At(trIDBx);
      noTRDPoints = pMCtrack->GetTrdPoints();

      if(noTRDPoints == 12){
      Double_t fPosX = pPointDx->GetX();
      fDistX->Fill(fabs(fPosX-x2));
      }
      }
      //--- end: for estimation dx and dy -------------------------------    
      */

      if (trRight->GetCoord1() > y2+dY) continue;
      if (trRight->GetCoord1() < y2-dY) break;

      if (first == 1){
	first = 0;
	iSpacePtStart = iSpacePtB;
      }

      if((trRight->GetCoord2() <= x2+dX) &&
      	 (trRight->GetCoord2() >= x2-dX)){
	
	iSegAcc14++;  //we have one more tracklet
	//we create a new tracklet object to add it into a vector
	clTr = new CbmL1TrdTracklet4();

  	clTr->SetCoord1(y1);
  	clTr->SetCoord2(x1);
  	clTr->SetCoord3(yB);
  	clTr->SetCoord4(xB);

	clTr->M[0] = y1z;
	clTr->M[1] = x1z;
	clTr->M[2] = yBz;
	clTr->M[3] = xBz;

	//zed is a z-axis value to which we extrapolate
	Double_t zed=0;
	  	  
	//cout << trdHit1->GetZ() <<" : "<< trRight->GetCoord4() << endl;
	if(Int_t(trRight->GetCoord4()) == fTrd14_Z) zed = fTrd21_Z; 
	if(Int_t(trRight->GetCoord4()) == fTrd24_Z) zed = fTrd31_Z;

	//if(trRight->GetPlanesID(1) == 4) zed = fTrd21_Z; 
	//if(trRight->GetPlanesID(1) == 8) zed = fTrd31_Z;

	//set the extrapolated coordinates - Y
	clTr->SetExt1(DistTwoTrackletsY(fArrayTrdHit,
					indSpacePtA, indSpacePtC,
					zed));
	//set the extrapolated coordinates - X
	clTr->SetExt2(DistTwoTrackletsX(fArrayTrdHit,
					indSpacePtB, indSpacePtD,
					zed));
	//set the indices of all 4 point in tracklet
	clTr->SetInd1(indSpacePtA);
	clTr->SetInd2(indSpacePtB);
	clTr->SetInd3(indSpacePtC);
	clTr->SetInd4(indSpacePtD);
	//initial value is always 0; it is modified if necessary
	clTr->SetVal(0);

	//position of the tracklet in tracklets vector
	clTr->SetIndex(iIndex);
	iIndex++;

	//add a tracklet to a vector
	clTrackletsAD.push_back(clTr);
	
	noXYSurv++;	
      }
    }
  }
  if(fVerbose > 2) {
    cout << " ----------- Tracklet 12-34 ------------------" << endl;
    cout << " Number of X survivors: " << noXSurv << endl;
    cout << " Number of Y survivors: " << noYSurv << endl;
    cout << "Number of XY survivors: " << noXYSurv << endl;
    cout << "   Number of all pairs: " << noAllPairs << endl;
  }
  return true;
}
// -----------------------------------------------------------------------

void CbmL1CATrdTrackFinderSA::CreateSpacePoints(vector <LayerWithHits> vTrdHitArrayA,
						vector <LayerWithHits> vTrdHitArrayB,
						vector <CbmL1TrdTracklet*>& clSpacePointsAB,
						Double_t sigmaA, Double_t sigmaB) {
  //method to create Space Points (SP) from two individual detector hist

  vector<LayerWithHits>::iterator itA;
  vector<LayerWithHits>::iterator itB;
  vector<LayerWithHits>::iterator itStart;

  Int_t
    noOverlapsAB = 0,
    //    noMCOverlapsAB = 0,
    iInd = 0;

  CbmL1TrdTracklet *clSpacePt;
    
  Int_t
    indA, indB;

  Double_t
    A_X, A_Y, A_Z,
    A_DX, A_DY;

  Double_t
    B_X, B_Y, B_Z,
    B_DX, B_DY;
 
 Int_t
    A_planeID,
    B_planeID;
 
  //   Int_t
  //     A_mcTrID,
  //     B_mcTrID;

  //   Double_t SPlength;
  
  Int_t first = 1;

  /*
    68,0%; <= 1sigma
    95,5%; <= 2sigma
    99,7%, <= 3sigma
  */

  itStart = vTrdHitArrayB.begin();

  for(itA  = vTrdHitArrayA.begin();
      itA != vTrdHitArrayA.end();
      itA++) {
    indA = (*itA).hitIndex;
    A_X  = (*itA).X;
    A_Y  = (*itA).Y;
    A_Z  = (*itA).Z;
    A_DX = (*itA).DX;
    A_DY = (*itA).DY;
    A_planeID = (*itA).planeID;
    //A_mcTrID = (*itA).mcTrackID;

    //     //--- extrapolate --------------------------------------------
    //     Double_t distBetweenLayer = fTrd14_Z-fTrd13_Z;
    //     //    cout << "distBetweenLayer: " << distBetweenLayer << endl;

    //     Double_t
    //       x1  = A_X,
    //       y1  = A_Y,
    //       z1y = A_Z,
    //       z1x = A_Z,
    //       x2  = 0,
    //       y2  = 0,
    //       z2y = z1y+distBetweenLayer, // 12 is a distance beetwen plane in TRD
    //       z2x = z1x+distBetweenLayer; // 12 is a distance beetwen plane in TRD

    //     y2 = (y1/z1y)*z2y;
    //     x2 = (x1/z1x)*z2x;
    //     //--- end: extrapolate --------------------------------------

    for(itB  = itStart;
	itB != vTrdHitArrayB.end();
	itB++) {
      indB = (*itB).hitIndex;
      B_X  = (*itB).X;
      B_Y  = (*itB).Y;
      B_Z  = (*itB).Z;
      B_DX = (*itB).DX;
      B_DY = (*itB).DY;
      B_planeID = (*itB).planeID;
      //B_mcTrID = (*itB).mcTrackID;

      if (B_Y + sigmaB*B_DY < A_Y - sigmaA*A_DY) continue;
      if (B_Y - sigmaB*B_DY > A_Y + sigmaA*A_DY) break;

      //       cout <<"(A_X, A_DX) "<< A_X <<", "<< A_DX << endl;
      //       cout <<"(A_Y, A_DY) "<< A_Y <<", "<< A_DY << endl;
      //       cout <<"(B_X, B_DX) "<< B_X <<", "<< B_DX << endl;
      //       cout <<"(B_Y, B_DY) "<< B_Y <<", "<< B_DY << endl;
      //       cout <<"-----------------------------------------"<< endl;
      //      if (B_Y < y2 - sigmaB*A_DY) continue;
      //      if (B_Y > y2 + sigmaB*A_DY) break;
      //      cout <<"[Extrapol, Real, Diff] "
      //   << y2 <<", "<< B_Y <<", "<< B_Y-y2 << endl;

      if (first == 1){
	first = 0;
	itStart = itB;
      }
      
      if((B_X-sigmaA*B_DX < A_X+sigmaB*A_DX) &&
	  (B_X+sigmaA*B_DX > A_X-sigmaB*A_DX)){	

	noOverlapsAB++;
	clSpacePt = new CbmL1TrdTracklet();
	
	clSpacePt->SetIndLeft(indA);
	clSpacePt->SetIndRight(indB);
	clSpacePt->SetVal(0);
	clSpacePt->SetIndex(iInd);
	clSpacePt->SetCoord1(A_Y);
	clSpacePt->SetCoord2(B_X);
	clSpacePt->SetCoord3(A_Z);
	clSpacePt->SetCoord4(B_Z);
	clSpacePt->SetPlanesID(A_planeID, B_planeID);

	// 	//----------------------------------------
	// 	TVector3 t1, t2, t3;
	// 	t1.SetXYZ(A_X, A_Y, A_Z);
	// 	t2.SetXYZ(B_X, B_Y, B_Z);
	// 	t3 = t2 - t1;
	// 	SPlength = t3.Mag();

	// 	Double_t
	// 	  x = 0,
	// 	  y = 0,
	// 	  z = 0,
	// 	  t = 0;
	
	// 	 t = (z - A_Z)/(B_Z - A_Z);
	// 	 x = A_X + t*(B_X-A_X);
	// 	 y = A_Y + t*(B_Y-A_Y);

	// 	//SPlength = sqrt(pow(B_Y-A_Y,2)+pow(B_Y-A_Y,2));
	// 	 if(Station == 1){
	// 	   if(A_mcTrID == B_mcTrID){
	// 	     noMCOverlapsAB++;
	// 	     fSPlengthMC->Fill(SPlength);
	// 	     fYat0MC->Fill(x);
	// 	   }else{
	// 	     fSPlength->Fill(SPlength);
	// 	     fYat0->Fill(x);
	// 	   }
	// 	 }
	iInd++;
	clSpacePointsAB.push_back(clSpacePt);
      }
    }
  }
  iInd = 0;
  cout <<"--- No Space Points: "<< noOverlapsAB << endl; 
  //   cout <<"--- (No SP, MCSP) "
  //        << noOverlapsAB <<", "
  //        << noMCOverlapsAB <<endl; 
}


// -----------------------------------------------------------------------
Bool_t CbmL1CATrdTrackFinderSA::CreateAndManageSegments(vector <CbmL1TrdTracklet4*> clTracklets14,
							vector <CbmL1TrdTracklet4*> clTracklets58,
							vector <CbmL1TrdTracklet4*> clTracklets912,
							TClonesArray *trackArray) { 
  //CreateAndManageSegments with no long track creation, 1 track = 1 segment
    
  vector<CbmL1TrdTracklet4*>::iterator itclSeg;

  CbmTrdTrack *tr1 = NULL;
  CbmTrdHit* hit1 = NULL;
  CbmTrdHit* hit2 = NULL;
  CbmTrdHit* hit3 = NULL;
  CbmTrdHit* hit4 = NULL;

  CbmL1TrdTracklet4 
    *clTrdSeg;
  
  Int_t
    trackArrayInd = 0;

  Int_t
    indTrdHit1,
    indTrdHit2,
    indTrdHit3,
    indTrdHit4;

  for(itclSeg = clTracklets14.begin();
      itclSeg != clTracklets14.end();
      itclSeg++) {
      
    clTrdSeg = *itclSeg;
    tr1 = new CbmTrdTrack();

    indTrdHit1 = clTrdSeg->GetInd1();
    hit1 = (CbmTrdHit*)fArrayTrdHit->At(indTrdHit1);
    tr1->AddHit(indTrdHit1, hit1);

    indTrdHit2 = clTrdSeg->GetInd2();
    hit2 = (CbmTrdHit*)fArrayTrdHit->At(indTrdHit2);
    tr1->AddHit(indTrdHit2, hit2);
  
    indTrdHit3 = clTrdSeg->GetInd3();
    hit3 = (CbmTrdHit*)fArrayTrdHit->At(indTrdHit3);
    tr1->AddHit(indTrdHit3, hit3);
    
    indTrdHit4 = clTrdSeg->GetInd4();
    hit4 = (CbmTrdHit*)fArrayTrdHit->At(indTrdHit4);
    tr1->AddHit(indTrdHit4, hit4);

    tr1->SortHits();
	
    new ((*fArrayTrdTrack)[trackArrayInd]) CbmTrdTrack(*tr1);

    delete tr1;
    trackArrayInd++;
  }

  for(itclSeg = clTracklets58.begin();
      itclSeg != clTracklets58.end();
      itclSeg++) {

    clTrdSeg = *itclSeg;
    tr1 = new CbmTrdTrack();
    
    indTrdHit1 = clTrdSeg->GetInd1();
    hit1 = (CbmTrdHit*)fArrayTrdHit->At(indTrdHit1);
    tr1->AddHit(indTrdHit1, hit1);

    indTrdHit2 = clTrdSeg->GetInd2();
    hit2 = (CbmTrdHit*)fArrayTrdHit->At(indTrdHit2);
    tr1->AddHit(indTrdHit2, hit2);
  
    indTrdHit3 = clTrdSeg->GetInd3();
    hit3 = (CbmTrdHit*)fArrayTrdHit->At(indTrdHit3);
    tr1->AddHit(indTrdHit3, hit3);
    
    indTrdHit4 = clTrdSeg->GetInd4();
    hit4 = (CbmTrdHit*)fArrayTrdHit->At(indTrdHit4);
    tr1->AddHit(indTrdHit4, hit4);

    tr1->SortHits();

    new ((*fArrayTrdTrack)[trackArrayInd]) CbmTrdTrack(*tr1);

    delete tr1;
    trackArrayInd++;
  }

  for(itclSeg = clTracklets912.begin();
      itclSeg != clTracklets912.end();
      itclSeg++) {

    clTrdSeg = *itclSeg;
    tr1 = new CbmTrdTrack();
    
    indTrdHit1 = clTrdSeg->GetInd1();
    hit1 = (CbmTrdHit*)fArrayTrdHit->At(indTrdHit1);
    tr1->AddHit(indTrdHit1, hit1);

    indTrdHit2 = clTrdSeg->GetInd2();
    hit2 = (CbmTrdHit*)fArrayTrdHit->At(indTrdHit2);
    tr1->AddHit(indTrdHit2, hit2);
  
    indTrdHit3 = clTrdSeg->GetInd3();
    hit3 = (CbmTrdHit*)fArrayTrdHit->At(indTrdHit3);
    tr1->AddHit(indTrdHit3, hit3);
    
    indTrdHit4 = clTrdSeg->GetInd4();
    hit4 = (CbmTrdHit*)fArrayTrdHit->At(indTrdHit4);
    tr1->AddHit(indTrdHit4, hit4);

    tr1->SortHits();

    new ((*fArrayTrdTrack)[trackArrayInd]) CbmTrdTrack(*tr1);

    delete tr1;
    trackArrayInd++;
  }

  return true;
}
// -----------------------------------------------------------------------



// -----------------------------------------------------------------------
Double_t CbmL1CATrdTrackFinderSA::Fit(CbmTrdTrack *tr)
{//fits a track with a straight line and caluculates each point's deviation

  Double_t chi2 = 0;
  
  Int_t iHit;
  CbmTrdHit *pHit[12];

  Double_t
    y1,y2,
    z1,z2,
    yA[12], zA[12],
    yB[12], zB[12];

  Double_t dist1 = 0;
  Double_t dist2 = 0;
  //  Double_t sum1 = 0;
  //  Double_t sum2 = 0;
  Double_t yS1 = 0;
  Double_t yS2 = 0;

  Int_t noHits = tr->GetNofTrdHits();

  for(int i = 0; i < noHits; i++){
    iHit = tr->GetTrdHitIndex(i);

    pHit[i]  = (CbmTrdHit*) fArrayTrdHit->At(iHit);
  }

  //first and last X(orY) pointa in the track
  yA[0]  = pHit[0]->GetY();
  zA[0]  = pHit[0]->GetZ();
  yB[10] = pHit[10]->GetY();
  zB[10] = pHit[10]->GetZ();

  //first and last Y(orX) pointa in the track
  yA[1]  = pHit[1]->GetX();
  zA[1]  = pHit[1]->GetZ();
  yB[11] = pHit[11]->GetX();
  zB[11] = pHit[11]->GetZ();

  Int_t t1 = 0;
  Int_t t2 = 1;

  for(int i = 0; i < 4; i++){

    //for X(orY) coordinates
    t1+=2;
    y1 = fabs(pHit[t1]->GetY());
    z1 = pHit[t1]->GetZ();

    yS1 = fabs((((yB[10]-yA[0])*(z1-zA[0]))/(zB[10]-zA[0]))+yA[0]);
    dist1 = fabs(yS1-y1);

    //for Y(orX) coordinates
    t2+=2;
    y2 = fabs(pHit[t2]->GetX());
    z2 = pHit[t2]->GetZ();

    yS2 = fabs((((yB[11]-yA[1])*(z2-zA[1]))/(zB[11]-zA[1]))+yA[1]);
    dist2 = fabs(yS2-y2);

    chi2 += (dist1+dist2)/2;
  }

  //  cout << "chi2: " << chi2 << endl;


  //   for(int i = 0; i < 4; i++){

  //     //for X(orY) coordinates
  //     t1+=2;
  //     y1 = fabs(pHit[t1]->Y());
  //     z1 = pHit[t1]->Z();

  //     yS1 = fabs((((yB[10]-yA[0])*(z1-zA[0]))/(zB[10]-zA[0]))+yA[0]);
  //     dist1 = fabs(yS1-y1);
  //     sum1+=dist1;    

  //     //for Y(orX) coordinates
  //     t2+=2;
  //     y2 = fabs(pHit[t2]->X());
  //     z2 = pHit[t2]->Z();

  //     yS2 = fabs((((yB[11]-yA[1])*(z2-zA[1]))/(zB[11]-zA[1]))+yA[1]);
  //     dist2 = fabs(yS2-y2);
  //     sum2+=dist2;    
  //   }
  //  chi2=sum1+sum2;



  // Calculate length of the track ---------------
  //   Double_t
  //     length0 = 0,
  //     length1 = 0,
  //     length2 = 0;

  //   t1 = 0;
  //   t2 = 1;

  //   Double_t
  //     x = 0,
  //     y = 0,
  //     z = 0;

  //   for(int i = 0; i < 4; i++){
  //     y = fabs(pHit[t1+2]->Y())-fabs(pHit[t1]->Y());
  //     z = (pHit[t1+2]->Z())-(pHit[t1]->Z());
  //     length1 += sqrt(pow(y,2) + pow(z,2));
  //     t1+=2;

  //     x = fabs(pHit[t2+2]->X())-fabs(pHit[t2]->X());
  //     z = (pHit[t2+2]->Z())-(pHit[t2]->Z());
  //     length2 += sqrt(pow(x,2) + pow(z,2));
  //     t2+=2;
  //   }

  //   width0 = width1 + width2;
  //   //  cout << width0 << endl;

  return chi2;
  //return length0;
}
// -----------------------------------------------------------------------

Double_t CbmL1CATrdTrackFinderSA::Fit(Int_t M[])
{//fits a track with a straight line and caluculates each point's deviation

  Double_t chi2 = 0;

  Int_t iHit;
  CbmTrdHit *pHit[12];

  Double_t
    y1,y2,
    z1,z2,
    yA[12], zA[12],
    yB[12], zB[12];

  Double_t dist1 = 0;
  Double_t dist2 = 0;
  Double_t yS1 = 0;
  Double_t yS2 = 0;

  Int_t noHits = 12;

  for(int i = 0; i < noHits; i++){
    iHit = M[i];
    pHit[i]  = (CbmTrdHit*) fArrayTrdHit->At(iHit);
  }

  //first and last X(orY) pointa in the track
  yA[0]  = pHit[0]->GetY();
  zA[0]  = pHit[0]->GetZ();
  yB[10] = pHit[10]->GetY();
  zB[10] = pHit[10]->GetZ();

  //first and last Y(orX) pointa in the track
  yA[1]  = pHit[1]->GetX();
  zA[1]  = pHit[1]->GetZ();
  yB[11] = pHit[11]->GetX();
  zB[11] = pHit[11]->GetZ();

  Int_t t1 = 0;
  Int_t t2 = 1;

  for(int i = 0; i < 4; i++){

    //for X(orY) coordinates
    t1+=2;
    y1 = fabs(pHit[t1]->GetY());
    z1 = pHit[t1]->GetZ();

    yS1 = fabs((((yB[10]-yA[0])*(z1-zA[0]))/(zB[10]-zA[0]))+yA[0]);
    dist1 = fabs(yS1-y1);

    //for Y(orX) coordinates
    t2+=2;
    y2 = fabs(pHit[t2]->GetX());
    z2 = pHit[t2]->GetZ();

    yS2 = fabs((((yB[11]-yA[1])*(z2-zA[1]))/(zB[11]-zA[1]))+yA[1]);
    dist2 = fabs(yS2-y2);

    chi2 += (dist1+dist2)/2;
  }
  return sqrt(chi2);
}
// -----------------------------------------------------------------------

// -----------------------------------------------------------------------
Double_t CbmL1CATrdTrackFinderSA::FitLinear(Int_t M[], 
					    Int_t var = 1) {
  //fit using a least square method
  

  Int_t noHits = 12;

  Int_t iHit;
  CbmTrdHit *hit;
  
  Double_t C1[12];
  Double_t C2[12];
  Double_t  Z[12];

  Double_t n = 12;

  for(int i=0;i<12;i++) {
    C1[i] = 0;
    C2[i] = 0;
    Z [i] = 0;
  }

  //  Int_t ind = 0;

  for(int i=0;i<noHits;i++) {
    
    iHit = M[i];
    hit = (CbmTrdHit*)fArrayTrdHit->At(iHit);
    
    C1[i] = hit->GetX();
    C2[i] = hit->GetY();
    Z[i]  = hit->GetZ();
  }

  Double_t 
    sumXZ = 0, 
    sumYZ = 0, 
    sumX  = 0, 
    sumY  = 0, 
    sumZx = 0, 
    sumZy = 0, 
    sumX2 = 0,
    sumY2 = 0,
    sumZx2 = 0,
    sumZy2 = 0;  

  for(int i=0;i<12;i+=2) {
    //    cout << "C1[" << i << "]=" << C1[i] << endl; 
    //    cout << "C2[" << i << "]=" << C2[i] << endl; 
    
    //    if(!(i % 2)) 
    {
      sumXZ += C1[i]*Z[ i];
      
      sumX  += C1[i];
      sumX2 += C1[i]*C1[i];

      sumZx  += Z[i];
      sumZx2 += Z[i]*Z[ i];

      
    }
  }

  for(int i=1;i<12;i+=2) {
    //    cout << "C1[" << i << "]=" << C1[i] << endl; 
    //    cout << "C2[" << i << "]=" << C2[i] << endl; 
    
    //    if(!((1i % 2)) 
    {
      sumYZ += C2[i]*Z[ i];
      
      sumY  += C2[i];
      sumY2 += C2[i]*C2[i];
      
      sumZy  += Z[ i];
      sumZy2 += Z[ i]*Z[ i];
    }
  }


  cout << "";
  //   Double_t a = (n*sumXZ - sumX*sumZ)/(n*sumX2 - sumX*sumX);
  //   Double_t b = (sumZ - a*sumX)/n;
  Double_t r1 = (n*sumXZ - sumX*sumZx)/sqrt((n*sumX2-sumX*sumX)*(n*sumZx2-sumZx*sumZx));
  //  Double_t Sa = sqrt((n*(sumY2 -a*sumXZ-b*sumY))/((n-2)*(n*sumX2-sumX*sumX)));
  //  Double_t Sb = Sa*sqrt(sumX2/n);


  Double_t r2 = (n*sumYZ - sumY*sumZy)/sqrt((n*sumY2-sumY*sumY)*(n*sumZy2-sumZy*sumZy));


  //   cout << "Linear coefficients: a= " 
  //        << a << ", Sa= " << Sa << ", b= " << b << ", Sb= " << Sb << endl 
  //        << ", r=" << r << ", chi2=" << tr->GetChi2() << endl;
    

  return sqrt(r1*r1+r2*r2);
}

Double_t CbmL1CATrdTrackFinderSA::FitKF(Int_t M[]) {

  Double_t Chi2 = 0;
  Double_t sigma2 = 100;
  //Double_t dZ;
  Double_t y;
  //Double_t ty;

  CbmTrdHit *trdHit;
  trdHit = (CbmTrdHit*) fArrayTrdHit->At(M[0]);
  y = trdHit->GetY();

  //  ty = 

  TMatrixFSym cMat(2);
  //  for(int i = 0; i < 2; i++) {
  // for(int j = 0; j < 2; i++) {
  //      cMat(i,j);
  //  }
  //}
  cMat(0,1) = sigma2;
  cMat(1,1) = 0;
  cMat(1,0) = 0;
  cMat(1,1) = 1E+10;






  return Chi2;
}

Double_t CbmL1CATrdTrackFinderSA::FitTLinearFitter(Int_t M[]) {

  Double_t chi2 = 0;
  const Int_t nrPnts = 12;
  Int_t k = 0;
  Int_t w = 0;
  Int_t j = 6;
  Double_t ax[6];
  Double_t ay[6];

  CbmTrdHit *trdHit;
  for(int i = 0; i < nrPnts; i++) {
    trdHit = (CbmTrdHit*) fArrayTrdHit->At(M[i]);
    k = i+1;
    if((k%2) == 0) {
      ax[w] = trdHit->GetX();
      w++;
    }else{
      ay[w] = trdHit->GetY();
    }
  }

  TLinearFitter *lf = new TLinearFitter(2);
  lf->SetFormula("hyp2");
  lf->StoreData(0);

  TVectorD x; x.Use(j,ax);
  TVectorD y; y.Use(j,ay);
  lf->AssignData(j,1,ax,ay);
  lf->Eval();
  chi2 = lf->GetChisquare();
  
  delete lf;

  return chi2;
}


ClassImp(CbmL1CATrdTrackFinderSA);