[u/mrichter/AliRoot.git] / HLT / TPCLib / tracking-ca / AliHLTTPCCATracker.cxx

// @(#) $Id$
//*************************************************************************
// This file is property of and copyright by the ALICE HLT Project        * 
// ALICE Experiment at CERN, All rights reserved.                         *
//                                                                        *
// Primary Authors: Jochen Thaeder <thaeder@kip.uni-heidelberg.de>        *
//                  Ivan Kisel <kisel@kip.uni-heidelberg.de>              *
//                  for The ALICE HLT Project.                            *
//                                                                        *
// Permission to use, copy, modify and distribute this software and its   *
// documentation strictly for non-commercial purposes is hereby granted   *
// without fee, provided that the above copyright notice appears in all   *
// copies and that both the copyright notice and this permission notice   *
// appear in the supporting documentation. The authors make no claims     *
// about the suitability of this software for any purpose. It is          *
// provided "as is" without express or implied warranty.                  *
//*************************************************************************

#include "AliHLTTPCCATracker.h"

#include "AliHLTTPCCAHit.h"
#include "AliHLTTPCCACell.h"
#include "AliHLTTPCCAOutTrack.h"

#include "TMath.h"
//#include "Riostream.h"
#include <vector>
#ifndef __SUNPRO_CC
#include <algo.h>
#endif

//#define DRAW

#ifdef DRAW
#include "AliHLTTPCCADisplay.h"
#include "TApplication.h"
#include "TROOT.h"
#include "TSystem.h"
#endif //DRAW

ClassImp(AliHLTTPCCATracker);


AliHLTTPCCATracker::AliHLTTPCCATracker()
  :fParam(),fRows(0),fOutTrackHits(0),fNOutTrackHits(0),fOutTracks(0),fNOutTracks(0),fTrackCells(0),fNHitsTotal(0),fTracks(0),fNTracks(0)
{
  // constructor
  fRows = new AliHLTTPCCARow[fParam.NRows()];
  Initialize( fParam );
}

AliHLTTPCCATracker::AliHLTTPCCATracker( const AliHLTTPCCATracker& )
  :fParam(),fRows(0),fOutTrackHits(0),fNOutTrackHits(0),fOutTracks(0),fNOutTracks(0),fTrackCells(0),fNHitsTotal(0),fTracks(0),fNTracks(0)
{
  // dummy
}

AliHLTTPCCATracker &AliHLTTPCCATracker::operator=( const AliHLTTPCCATracker& )
{
  // dummy
  return *this;
}

AliHLTTPCCATracker::~AliHLTTPCCATracker()
{
  // destructor
  StartEvent();
  delete[] fRows;
}

// ----------------------------------------------------------------------------------
void AliHLTTPCCATracker::Initialize( AliHLTTPCCAParam &param )
{
  // initialosation
  StartEvent();
  delete[] fRows;
  fRows = 0;
  fParam = param;
  fParam.Update();
  fRows = new AliHLTTPCCARow[fParam.NRows()];
  for( Int_t irow=0; irow<fParam.NRows(); irow++ ){
    fRows[irow].X() = fParam.RowXFirst() + irow*fParam.RowXStep();
  }
  StartEvent();
}

void AliHLTTPCCATracker::StartEvent()
{
  // start new event and fresh the memory  
  delete[] fTracks;
  delete[] fTrackCells;
  delete[] fOutTrackHits;
  delete[] fOutTracks;
  fTracks = 0;
  fTrackCells = 0;
  fOutTrackHits = 0;
  fOutTracks = 0;
  fNTracks = 0;
  fNOutTrackHits = 0;
  fNOutTracks = 0;
  fNHitsTotal = 0;
  for( Int_t irow=0; irow<fParam.NRows(); irow++ ){
    fRows[irow].Clear();
  }
}


void AliHLTTPCCATracker::ReadHitRow( Int_t iRow, AliHLTTPCCAHit *Row, Int_t NHits )
{
  // read row of hits
  AliHLTTPCCARow &row = fRows[iRow];
  row.Hits() = new AliHLTTPCCAHit[NHits];
  for( Int_t i=0; i<NHits; i++ ){ 
    row.Hits()[i]=Row[i];
    row.Hits()[i].ErrY()*= fParam.YErrorCorrection();
    row.Hits()[i].ErrZ()*= fParam.ZErrorCorrection();
  }
  row.NHits() = NHits;
  fNHitsTotal += NHits;
}

void AliHLTTPCCATracker::Reconstruct()
{
  // reconstruction of event
  FindCells();
  FindTracks();
}


void AliHLTTPCCATracker::FindCells()
{
  // cell finder - neighbouring hits are grouped to cells

  for( Int_t irow=0; irow<fParam.NRows(); irow++ ){
    AliHLTTPCCARow &row=fRows[irow];
    Int_t nHits = row.NHits();
    if( nHits<1 ) continue;
    row.CellHitPointers() = new Int_t[ nHits ];
    Int_t nPointers = 0;
    std::vector<AliHLTTPCCACell> vCells; 
    Bool_t vUsed[10000];
    for (Int_t ih = 0; ih<nHits; ih++) vUsed[ih] = 0;
  
    for (Int_t ih = 0; ih<nHits; ih++){    
      if( vUsed[ih] ) continue;
      // cell start
      AliHLTTPCCACell cell;
      cell.FirstHitRef() = nPointers;
      cell.NHits() = 1;
      cell.IDown() = -1;
      cell.IUp() = -1;
      cell.IUsed() = 0;
      row.CellHitPointers()[nPointers++] = ih;
      vUsed[ih] = 1;
      Int_t jLast = ih;
      while( jLast<nHits-1 ){
	AliHLTTPCCAHit &h  = row.Hits()[jLast];    
	Double_t d2min = 1.e10;
	Int_t jBest = -1;
	for (Int_t j = jLast+1; j<nHits; j++){    
	  if( vUsed[j] ) continue;
	  AliHLTTPCCAHit &h1  = row.Hits()[j];
	  Double_t dy = TMath::Abs(h.Y() - h1.Y() );
	  if( dy>(h.ErrY()+h1.ErrY())*4. ) break;
	  Double_t dz = TMath::Abs(h.Z() - h1.Z() );
	  if( dz>(h.ErrZ()+h1.ErrZ())*4. ) continue;
	  Double_t d2 = dz*dz+dy*dy;
	  if( d2<d2min ){
	    d2min = d2;
	    jBest = j;
	  }
	}
	if( jBest<0 ) break;
	row.CellHitPointers()[nPointers++] = jBest;
	cell.NHits()++;
	vUsed[jBest] = 1;     
	jLast = jBest;
      }

      AliHLTTPCCAHit &h = row.GetCellHit(cell,0);
      AliHLTTPCCAHit &h1= row.GetCellHit(cell,cell.NHits()-1);
    
      cell.Y() = .5*(h.Y() + h1.Y());
      cell.Z() = .5*(h.Z() + h1.Z());
      cell.ErrY() = .5*(TMath::Abs(h.Y() - h1.Y())/3 + h.ErrY() + h1.ErrY());
      cell.ErrZ() = .5*(TMath::Abs(h.Z() - h1.Z())/3 + h.ErrZ() + h1.ErrZ());
      vCells.push_back(cell);
    }
    
    row.Cells() = new AliHLTTPCCACell[vCells.size()];
    row.NCells() = vCells.size();
    for( Int_t i=0; i<vCells.size(); i++ ) row.Cells()[i]=vCells[i];  
  }
}


void AliHLTTPCCATracker::FindTracks()
{
  // the Cellular Automaton track finder
  
  if( fNHitsTotal < 1 ) return;

#ifdef DRAW
  if( !gApplication ){
    TApplication *myapp = new TApplication("myapp",0,0);
  }
  AliHLTTPCCADisplay::Instance().SetCurrentSector( this );
  AliHLTTPCCADisplay::Instance().DrawSector( this );
  //for (Int_t i = 0; i<fHits.size(); i++) AliHLTTPCCADisplay::Instance().DrawHit( i );
  //cout<<"hits"<<endl;
  //AliHLTTPCCADisplay::Instance().Ask();
  //AliHLTTPCCADisplay::Instance().Clear();
  //AliHLTTPCCADisplay::Instance().DrawSector( this );
  for( Int_t iRow=0; iRow<fParam.NRows(); iRow++ )
    for (Int_t i = 0; i<fRows[iRow].NCells(); i++) 
      AliHLTTPCCADisplay::Instance().DrawCell( iRow, i );
  //cout<<"cells"<<endl;
  AliHLTTPCCADisplay::Instance().Ask();
  Int_t nConnectedCells = 0;
#endif 

  std::vector<AliHLTTPCCATrack> vTracks; 
  std::vector<Int_t> vTrackCells; 
  fTrackCells = new Int_t[2*fNHitsTotal];

  for( Int_t iRow1=0; iRow1<fParam.NRows()-1; iRow1++ ){
    AliHLTTPCCARow &row1 = fRows[iRow1];
    Int_t lastRow2 = iRow1+3;
    if( lastRow2>=fParam.NRows() ) lastRow2 = fParam.NRows()-1;    
    for( Int_t iRow2=iRow1+1; iRow2<=lastRow2; iRow2++ ){
      AliHLTTPCCARow &row2 = fRows[iRow2];
      for (Int_t i1 = 0; i1<row1.NCells(); i1++){
	AliHLTTPCCACell *c1  = &(row1.Cells()[i1]);
	if( c1->IUp()>=0 ) continue;
	Double_t sy1 = c1->ErrY()*c1->ErrY();
	Double_t sz1 = c1->ErrZ()*c1->ErrZ();
	for (Int_t i2 = 0; i2<row2.NCells(); i2++){
	  AliHLTTPCCACell *c2  = &(row2.Cells()[i2]);
	  Double_t sy2 = c2->ErrY()*c2->ErrY();
	  Double_t sz2 = c2->ErrZ()*c2->ErrZ();
	  Double_t dist1 = sqrt( (c1->Y()-c2->Y())*(c1->Y()-c2->Y())/(sy1+sy2) );
	  Double_t dist2 = sqrt( (c1->Z()-c2->Z())*(c1->Z()-c2->Z())/(sz1+sz2) );
	  if( dist1>3.5*fParam.CellConnectionFactor() ) continue;
	  if( dist2>3.5*fParam.CellConnectionFactor() ) continue;
	  if( c1->IUp() ==-1 ) c1->IUp() = (i2<<8)+iRow2;
	  else c1->IUp() = -2;
	  if( c2->IDown() ==-1 ) c2->IDown() = (i1<<8)+iRow1;
	  else c2->IDown() = -2;
	}
      }
    }
  }

  Int_t nOutTrackHits = 0;
  Int_t nTrackCells = 0;
  for( Int_t iRow1=0; iRow1<fParam.NRows(); iRow1++ ){
    AliHLTTPCCARow &row1 = fRows[iRow1];
    for (Int_t i1 = 0; i1<row1.NCells(); i1++){ 
      AliHLTTPCCACell *c1  = &(row1.Cells()[i1]);
      if( c1->IDown()==-2 || c1->IUp()==-2 ) continue;
      if( c1->IUsed()>0 ) continue;
      c1->IUsed() = 1;
      AliHLTTPCCATrack track;
      track.Used() = 0;
      track.NCells() = 1;
      track.IFirstCell() = nTrackCells;
      fTrackCells[nTrackCells++] = (i1<<8)+iRow1;
      AliHLTTPCCACell *last = c1;
      Int_t lastRow = iRow1;
      while( last->IUp() >=0 ){
	Int_t iRow2 = last->IUp()%256;
	AliHLTTPCCARow &row2 = fRows[iRow2];
	AliHLTTPCCACell *next = &(row2.Cells()[last->IUp()>>8]);
	if( next->IDown()==-2 || next->IUp()==-2 ) break;
#ifdef DRAW
	AliHLTTPCCADisplay::Instance().ConnectCells( lastRow,*last,iRow2,*next );
	nConnectedCells++;
#endif 
	next->IUsed() = 1;      
	fTrackCells[nTrackCells++] = last->IUp();
	track.NCells()++;
	last = next;
	lastRow = iRow2;
      } 
      vTracks.push_back(track);
    }
  }
  
  Int_t nTracks = vTracks.size();
  std::sort( vTracks.begin(), vTracks.end(), AliHLTTPCCATrack::CompareSize);
  
  fTracks = new AliHLTTPCCATrack[nTracks];
  fNTracks = 0;
  vTrackCells.clear();
 
  Int_t * vMatchedTracks = new Int_t[nTracks];
  for( Int_t itr=0; itr<nTracks; itr++ ){
    AliHLTTPCCATrack &iTrack = vTracks[itr];
    if( iTrack.Used() ) continue;    
    FitTrack( iTrack );
    if( iTrack.Param().Chi2() > fParam.TrackChi2Cut()*iTrack.Param().NDF() ) continue;    

    Int_t iFirstRow =  GetTrackCellIRow( iTrack, 0);
    Int_t iLastRow  =  GetTrackCellIRow( iTrack, iTrack.NCells()-1);
    AliHLTTPCCACell *iFirstCell = &GetTrackCell( iTrack, 0);
    AliHLTTPCCACell *iLastCell = &GetTrackCell( iTrack, iTrack.NCells()-1);
    
    Bool_t updated = 1;
    Int_t nMatched = 0;
    std::vector<Int_t> vMatchedCells;
    while( updated ){
      updated = 0;
      for( Int_t jtr=1; jtr<nTracks; jtr++ ){
	AliHLTTPCCATrack &jTrack = vTracks[jtr];
	if( jTrack.Used() ) continue;
	Int_t jFirstRow =  GetTrackCellIRow( jTrack, 0);
	Int_t jLastRow =  GetTrackCellIRow( jTrack, jTrack.NCells()-1);
	AliHLTTPCCACell *jFirstCell = &GetTrackCell( jTrack, 0);
	AliHLTTPCCACell *jLastCell = &GetTrackCell( jTrack, jTrack.NCells()-1);

	Int_t dFirstRow1 = TMath::Abs(iFirstRow-jLastRow);
	Int_t dFirstRow2 = TMath::Abs(iFirstRow-jFirstRow);
	Int_t dLastRow1 = TMath::Abs(iLastRow-jLastRow);
	Int_t dLastRow2 = TMath::Abs(iLastRow-jFirstRow);
	if( dFirstRow1 > fParam.MaxTrackMatchDRow() && 
	    dFirstRow2 > fParam.MaxTrackMatchDRow() &&
	    dLastRow1  > fParam.MaxTrackMatchDRow() && 
	    dLastRow2  > fParam.MaxTrackMatchDRow()    ) continue;
	Int_t iCase=0;
	AliHLTTPCCACell *iC, *jC;
	if( dFirstRow1<dFirstRow2 && dFirstRow1<dLastRow1 && dFirstRow1<dLastRow2 ){
	  iCase = 0;
	  iC = iFirstCell;
	  jC = jLastCell;
	}else if( dFirstRow2<dLastRow1 && dFirstRow2<dLastRow2 ){
	  iCase = 1;
	  iC = iFirstCell;
	  jC = jFirstCell;
	}else if( dLastRow1<dLastRow2 ){
	  iCase = 2;
	  iC = iLastCell; 
	  jC = jLastCell;
	}else{
	  iCase = 3;
	  iC = iLastCell; 
	  jC = jFirstCell;
	}
	{
	  Double_t dy = TMath::Abs(iC->Y() - jC->Y());
	  Double_t dz = TMath::Abs(iC->Z() - jC->Z());
	  Double_t sy1 = iC->ErrY()*iC->ErrY();
	  Double_t sz1 = iC->ErrZ()*iC->ErrZ();
	  Double_t sy2 = jC->ErrY()*jC->ErrY();
	  Double_t sz2 = jC->ErrZ()*jC->ErrZ();
	  Double_t dist1 = sqrt( (dy)/(sy1+sy2) );
	  Double_t dist2 = sqrt( (dz)/(sz1+sz2) );
	  if( dist1>3.5*5*fParam.CellConnectionFactor() ) continue;
	  if( dist2>3.5*5*fParam.CellConnectionFactor() ) continue;
	}
	AliHLTTPCCATrackPar t = iTrack.Param();
	//t.Chi2() = 0;
	//t.NDF() = 0;
	for( Int_t i=0; i<jTrack.NCells() ; i++){
	  AliHLTTPCCACell &c = GetTrackCell(jTrack,i);	
	  AliHLTTPCCARow &row = GetTrackCellRow(jTrack,i);
	  for( Int_t j=0; j<c.NHits(); j++){
	    AliHLTTPCCAHit &h = row.GetCellHit(c,j);
	    Double_t m[3] = {row.X(), h.Y(), h.Z() };
	    Double_t mV[6] = {fParam.ErrX()*fParam.ErrX(), 0, h.ErrY()*h.ErrY(), 0, 0, h.ErrZ()*h.ErrZ() };
	    Double_t mV1[6];
	    t.TransportBz(fParam.Bz(),m);
	    t.GetConnectionMatrix(fParam.Bz(),m, mV1);
	    t.Filter(m, mV, mV1);
	  }
	}
	if( t.Chi2() > fParam.TrackChi2Cut()*t.NDF() ) continue;
	if( iCase==0 ){
	  iFirstRow = jFirstRow;
	  iFirstCell = jFirstCell;
	}else if( iCase ==1 ){
	  iFirstRow = jLastRow;
	  iFirstCell = jLastCell;
	}else if( iCase == 2 ){
	  iLastRow = jFirstRow;
	  iLastCell = jFirstCell;
	}else{
	  iLastRow = jLastRow;
	  iLastCell = jLastCell;
	}
	t.Normalize();

	for( Int_t i=0; i<jTrack.NCells(); i++){
	  vMatchedCells.push_back(fTrackCells[jTrack.IFirstCell()+i]);
	}
	//t.NDF()+= iTrack.Param().NDF();
	//t.Chi2()+= iTrack.Param().Chi2();
	iTrack.Param() = t;
	vMatchedTracks[nMatched++] = jtr;
	jTrack.Used()=1;
	updated = 1;
	break;
      }
    }

    if(0){
      Double_t t0[7];
      for( Int_t i=0; i<7; i++ ) t0[i] = iTrack.Param().Par()[i];
      iTrack.Param().Init();
      for( Int_t i=0; i<7; i++ ) iTrack.Param().Par()[i] = t0[i];
      for( Int_t i=0; i<iTrack.NCells() ; i++){
	AliHLTTPCCACell &c = GetTrackCell( iTrack, i);
	AliHLTTPCCARow &row = GetTrackCellRow( iTrack, i);
	for( Int_t j=0; j<c.NHits(); j++){
	  AliHLTTPCCAHit &h = row.GetCellHit(c,j);
	  Double_t m[3] = {row.X(), h.Y(), h.Z() };
	  Double_t mV[6] = {0, 0, h.ErrY()*h.ErrY(), 0, 0, h.ErrZ()*h.ErrZ() };
	  Double_t mV1[6];
	  iTrack.Param().TransportBz(fParam.Bz(), m, t0);
	  iTrack.Param().GetConnectionMatrix(fParam.Bz(), m, mV1, t0);
	  iTrack.Param().Filter(m, mV, mV1);
	}
      }
      iTrack.Param().Normalize();
    }
    //FitTrack(iTrack,5);
    Int_t nHits = 0;
    for( Int_t iCell=0; iCell<iTrack.NCells(); iCell++){
      Int_t ind = fTrackCells[iTrack.IFirstCell()+iCell];
      AliHLTTPCCARow &row = fRows[ind%256];
      AliHLTTPCCACell &c = row.Cells()[ind>>8];
      nHits+=c.NHits();
    }
    for( UInt_t i=0; i<vMatchedCells.size(); i++){
      Int_t ind = vMatchedCells[i];
      AliHLTTPCCARow &row = fRows[ind%256];
      AliHLTTPCCACell &c = row.Cells()[ind>>8];
      nHits+=c.NHits();
    }

    if( nHits<5 ){
      for( Int_t i=0; i<nMatched; i++ ) vTracks[vMatchedTracks[i]].Used()=0;
      continue;
    }
    iTrack.Used() = 1;
    Int_t oldNCells = vTrackCells.size();
    for( Int_t i=0; i<iTrack.NCells() ; i++){
      vTrackCells.push_back(fTrackCells[iTrack.IFirstCell()+i]);
    }      
    iTrack.IFirstCell() = oldNCells;
    for( UInt_t i=0; i<vMatchedCells.size(); i++){
      vTrackCells.push_back(vMatchedCells[i]);
      iTrack.NCells()++;
    }      
    fTracks[fNTracks++] = iTrack;  
    nOutTrackHits+= nHits;  
  }
  delete [] vMatchedTracks;

  //fTrackCells = new Int_t[vTrackCells.size()];
  for( UInt_t i=0; i<vTrackCells.size(); i++ ) fTrackCells[i] = vTrackCells[i];
  
#ifdef DRAW
  if( nConnectedCells>0 ) AliHLTTPCCADisplay::Instance().Ask();
#endif 


  fOutTrackHits = new Int_t[nOutTrackHits];
  fNOutTrackHits = 0;
  fNOutTracks = fNTracks;
  fOutTracks = new AliHLTTPCCAOutTrack[fNOutTracks];
  for( Int_t itr=0; itr<fNOutTracks; itr++ ){
    AliHLTTPCCATrack &t = fTracks[itr];
#ifdef DRAW
  AliHLTTPCCADisplay::Instance().DrawTrack( t );
#endif 
    AliHLTTPCCAOutTrack &tmp = fOutTracks[itr];
    tmp.FirstHitRef() = fNOutTrackHits;
    tmp.NHits() = 0;
    tmp.Param() = t.Param();
    for( Int_t iCell=0; iCell<t.NCells(); iCell++){
      AliHLTTPCCACell &cell = GetTrackCell(t,iCell);
      AliHLTTPCCARow &row = GetTrackCellRow(t,iCell);
      for( Int_t iHit=0; iHit<cell.NHits(); iHit++ ){
	AliHLTTPCCAHit &hit = row.GetCellHit(cell,iHit);
	fOutTrackHits[fNOutTrackHits] = hit.ID();
	fNOutTrackHits++;
	tmp.NHits()++;
      }
    }
  }  
  
#ifdef DRAW
  AliHLTTPCCADisplay::Instance().Ask();
  //AliHLTTPCCADisplay::Instance().DrawMCTracks(fParam.fISec);
  //AliHLTTPCCADisplay::Instance().Update();
  //AliHLTTPCCADisplay::Instance().Ask();
#endif 
}


void AliHLTTPCCATracker::FitTrack( AliHLTTPCCATrack &track, Int_t nIter )
{    
  // fit the track with nIter iterations

  AliHLTTPCCATrackPar &t = track.Param();
  t.Init();
  
  AliHLTTPCCACell &c1 = GetTrackCell(track,0);
  AliHLTTPCCACell &c2 = GetTrackCell(track,track.NCells()-1);
  AliHLTTPCCARow &row1 = GetTrackCellRow(track,0);
  AliHLTTPCCARow &row2 = GetTrackCellRow(track,track.NCells()-1);
  Double_t t0[7];
  t0[0]=row1.X();
  t0[1]=c1.Y();
  t0[2]=c1.Z();
  t0[3]= row2.X() - row1.X();
  t0[4]= c2.Y() - c1.Y();
  t0[5]= c2.Z() - c1.Z();
  Double_t tt = sqrt(t0[3]*t0[3]+t0[4]*t0[4]+t0[5]*t0[5]);
  if( TMath::Abs(tt)>1.e-4 ){
    t0[3]/=tt;
    t0[4]/=tt;
    t0[5]/=tt;
  }else{
    t0[4]=1;
  }
  t0[6] = 0;

  for( Int_t iter=0; iter<nIter; iter++ ){
    t.Init();
    for( Int_t i=0; i<7; i++) t.Par()[i] = t0[i];
    {
      Double_t m[3] = {row1.X(), c1.Y(), c1.Z() };
      t.TransportBz(fParam.Bz(),m,t0);
    }
    t.Init();
    for( Int_t i=0; i<7; i++ ) t.Par()[i] = t0[i];
    
    for( Int_t i=0; i<track.NCells() ; i++){
      AliHLTTPCCACell &c = GetTrackCell(track,i);
      AliHLTTPCCARow &row = GetTrackCellRow(track,i);
      for( Int_t j=0; j<c.NHits(); j++){
	AliHLTTPCCAHit &h = row.GetCellHit(c,j);
	Double_t m[3] = {row.X(), h.Y(), h.Z() };
	Double_t mV[6] = {fParam.ErrX()*fParam.ErrX(), 0, h.ErrY()*h.ErrY(), 0, 0, h.ErrZ()*h.ErrZ() };
	Double_t mV1[6];
	t.TransportBz(fParam.Bz(),m, t0);
	t.GetConnectionMatrix(fParam.Bz(),m, mV1, t0);
	t.Filter(m, mV, mV1);
      }
    }
    t.Normalize();
    for( Int_t i=0; i<7; i++ ) t0[i] = t.Par()[i];
  }
}
Commit	Line	Data
326c2d4b	1	// @(#) $Id$
	2	//*************************************************************************
	3	// This file is property of and copyright by the ALICE HLT Project *
	4	// ALICE Experiment at CERN, All rights reserved. *
	5	// *
	6	// Primary Authors: Jochen Thaeder <thaeder@kip.uni-heidelberg.de> *
	7	// Ivan Kisel <kisel@kip.uni-heidelberg.de> *
	8	// for The ALICE HLT Project. *
	9	// *
	10	// Permission to use, copy, modify and distribute this software and its *
	11	// documentation strictly for non-commercial purposes is hereby granted *
	12	// without fee, provided that the above copyright notice appears in all *
	13	// copies and that both the copyright notice and this permission notice *
	14	// appear in the supporting documentation. The authors make no claims *
	15	// about the suitability of this software for any purpose. It is *
	16	// provided "as is" without express or implied warranty. *
	17	//*************************************************************************
	18
	19	#include "AliHLTTPCCATracker.h"
	20
	21	#include "AliHLTTPCCAHit.h"
	22	#include "AliHLTTPCCACell.h"
	23	#include "AliHLTTPCCAOutTrack.h"
	24
	25	#include "TMath.h"
	26	//#include "Riostream.h"
	27	#include <vector>
ce622827	28	#ifndef __SUNPRO_CC
326c2d4b	29	#include <algo.h>
ce622827	30	#endif
326c2d4b	31
	32	//#define DRAW
	33
	34	#ifdef DRAW
	35	#include "AliHLTTPCCADisplay.h"
	36	#include "TApplication.h"
	37	#include "TROOT.h"
	38	#include "TSystem.h"
	39	#endif //DRAW
	40
	41	ClassImp(AliHLTTPCCATracker);
	42
	43
	44	AliHLTTPCCATracker::AliHLTTPCCATracker()
	45	:fParam(),fRows(0),fOutTrackHits(0),fNOutTrackHits(0),fOutTracks(0),fNOutTracks(0),fTrackCells(0),fNHitsTotal(0),fTracks(0),fNTracks(0)
	46	{
	47	// constructor
	48	fRows = new AliHLTTPCCARow[fParam.NRows()];
	49	Initialize( fParam );
	50	}
	51
	52	AliHLTTPCCATracker::AliHLTTPCCATracker( const AliHLTTPCCATracker& )
	53	:fParam(),fRows(0),fOutTrackHits(0),fNOutTrackHits(0),fOutTracks(0),fNOutTracks(0),fTrackCells(0),fNHitsTotal(0),fTracks(0),fNTracks(0)
	54	{
	55	// dummy
	56	}
	57
	58	AliHLTTPCCATracker &AliHLTTPCCATracker::operator=( const AliHLTTPCCATracker& )
	59	{
	60	// dummy
ce622827	61	return *this;
326c2d4b	62	}
	63
	64	AliHLTTPCCATracker::~AliHLTTPCCATracker()
	65	{
	66	// destructor
	67	StartEvent();
	68	delete[] fRows;
	69	}
	70
	71	// ----------------------------------------------------------------------------------
	72	void AliHLTTPCCATracker::Initialize( AliHLTTPCCAParam &param )
	73	{
	74	// initialosation
	75	StartEvent();
	76	delete[] fRows;
	77	fRows = 0;
	78	fParam = param;
	79	fParam.Update();
	80	fRows = new AliHLTTPCCARow[fParam.NRows()];
	81	for( Int_t irow=0; irow<fParam.NRows(); irow++ ){
	82	fRows[irow].X() = fParam.RowXFirst() + irow*fParam.RowXStep();
	83	}
	84	StartEvent();
	85	}
	86
	87	void AliHLTTPCCATracker::StartEvent()
	88	{
	89	// start new event and fresh the memory
	90	delete[] fTracks;
	91	delete[] fTrackCells;
	92	delete[] fOutTrackHits;
	93	delete[] fOutTracks;
	94	fTracks = 0;
	95	fTrackCells = 0;
	96	fOutTrackHits = 0;
	97	fOutTracks = 0;
	98	fNTracks = 0;
	99	fNOutTrackHits = 0;
	100	fNOutTracks = 0;
	101	fNHitsTotal = 0;
	102	for( Int_t irow=0; irow<fParam.NRows(); irow++ ){
	103	fRows[irow].Clear();
	104	}
	105	}
	106
	107
	108	void AliHLTTPCCATracker::ReadHitRow( Int_t iRow, AliHLTTPCCAHit *Row, Int_t NHits )
	109	{
	110	// read row of hits
	111	AliHLTTPCCARow &row = fRows[iRow];
	112	row.Hits() = new AliHLTTPCCAHit[NHits];
	113	for( Int_t i=0; i<NHits; i++ ){
	114	row.Hits()[i]=Row[i];
	115	row.Hits()[i].ErrY()*= fParam.YErrorCorrection();
	116	row.Hits()[i].ErrZ()*= fParam.ZErrorCorrection();
	117	}
	118	row.NHits() = NHits;
	119	fNHitsTotal += NHits;
	120	}
	121
	122	void AliHLTTPCCATracker::Reconstruct()
	123	{
	124	// reconstruction of event
	125	FindCells();
126	FindTracks();
127	}
128
129
130	void AliHLTTPCCATracker::FindCells()
131	{
132	// cell finder - neighbouring hits are grouped to cells
133
134	for( Int_t irow=0; irow<fParam.NRows(); irow++ ){
135	AliHLTTPCCARow &row=fRows[irow];
136	Int_t nHits = row.NHits();
137	if( nHits<1 ) continue;
138	row.CellHitPointers() = new Int_t[ nHits ];
139	Int_t nPointers = 0;
140	std::vector<AliHLTTPCCACell> vCells;
141	Bool_t vUsed[10000];
142	for (Int_t ih = 0; ih<nHits; ih++) vUsed[ih] = 0;
143
144	for (Int_t ih = 0; ih<nHits; ih++){
145	if( vUsed[ih] ) continue;
146	// cell start
147	AliHLTTPCCACell cell;
148	cell.FirstHitRef() = nPointers;
149	cell.NHits() = 1;
150	cell.IDown() = -1;
151	cell.IUp() = -1;
152	cell.IUsed() = 0;
153	row.CellHitPointers()[nPointers++] = ih;
154	vUsed[ih] = 1;
155	Int_t jLast = ih;
156	while( jLast<nHits-1 ){
157	AliHLTTPCCAHit &h = row.Hits()[jLast];
158	Double_t d2min = 1.e10;
159	Int_t jBest = -1;
160	for (Int_t j = jLast+1; j<nHits; j++){
161	if( vUsed[j] ) continue;
162	AliHLTTPCCAHit &h1 = row.Hits()[j];
163	Double_t dy = TMath::Abs(h.Y() - h1.Y() );
164	if( dy>(h.ErrY()+h1.ErrY())*4. ) break;
165	Double_t dz = TMath::Abs(h.Z() - h1.Z() );
166	if( dz>(h.ErrZ()+h1.ErrZ())*4. ) continue;
167	Double_t d2 = dzdz+dydy;
168	if( d2<d2min ){
169	d2min = d2;
170	jBest = j;
171	}
172	}
173	if( jBest<0 ) break;
174	row.CellHitPointers()[nPointers++] = jBest;
175	cell.NHits()++;
176	vUsed[jBest] = 1;
177	jLast = jBest;
178	}
179
180	AliHLTTPCCAHit &h = row.GetCellHit(cell,0);
181	AliHLTTPCCAHit &h1= row.GetCellHit(cell,cell.NHits()-1);
182
183	cell.Y() = .5*(h.Y() + h1.Y());
184	cell.Z() = .5*(h.Z() + h1.Z());
185	cell.ErrY() = .5*(TMath::Abs(h.Y() - h1.Y())/3 + h.ErrY() + h1.ErrY());
186	cell.ErrZ() = .5*(TMath::Abs(h.Z() - h1.Z())/3 + h.ErrZ() + h1.ErrZ());
187	vCells.push_back(cell);
188	}
189
190	row.Cells() = new AliHLTTPCCACell[vCells.size()];
191	row.NCells() = vCells.size();
192	for( Int_t i=0; i<vCells.size(); i++ ) row.Cells()[i]=vCells[i];
193	}
194	}
195
196
197	void AliHLTTPCCATracker::FindTracks()
198	{
199	// the Cellular Automaton track finder
200
201	if( fNHitsTotal < 1 ) return;
202
203	#ifdef DRAW
204	if( !gApplication ){
205	TApplication *myapp = new TApplication("myapp",0,0);
206	}
207	AliHLTTPCCADisplay::Instance().SetCurrentSector( this );
208	AliHLTTPCCADisplay::Instance().DrawSector( this );
209	//for (Int_t i = 0; i<fHits.size(); i++) AliHLTTPCCADisplay::Instance().DrawHit( i );
210	//cout<<"hits"<<endl;
211	//AliHLTTPCCADisplay::Instance().Ask();
212	//AliHLTTPCCADisplay::Instance().Clear();
213	//AliHLTTPCCADisplay::Instance().DrawSector( this );
214	for( Int_t iRow=0; iRow<fParam.NRows(); iRow++ )
215	for (Int_t i = 0; i<fRows[iRow].NCells(); i++)
216	AliHLTTPCCADisplay::Instance().DrawCell( iRow, i );
217	//cout<<"cells"<<endl;
218	AliHLTTPCCADisplay::Instance().Ask();
219	Int_t nConnectedCells = 0;
220	#endif
221
222	std::vector<AliHLTTPCCATrack> vTracks;
223	std::vector<Int_t> vTrackCells;
224	fTrackCells = new Int_t[2*fNHitsTotal];
225
226	for( Int_t iRow1=0; iRow1<fParam.NRows()-1; iRow1++ ){
227	AliHLTTPCCARow &row1 = fRows[iRow1];
228	Int_t lastRow2 = iRow1+3;
229	if( lastRow2>=fParam.NRows() ) lastRow2 = fParam.NRows()-1;
230	for( Int_t iRow2=iRow1+1; iRow2<=lastRow2; iRow2++ ){
231	AliHLTTPCCARow &row2 = fRows[iRow2];
232	for (Int_t i1 = 0; i1<row1.NCells(); i1++){
233	AliHLTTPCCACell *c1 = &(row1.Cells()[i1]);
234	if( c1->IUp()>=0 ) continue;
235	Double_t sy1 = c1->ErrY()*c1->ErrY();
236	Double_t sz1 = c1->ErrZ()*c1->ErrZ();
237	for (Int_t i2 = 0; i2<row2.NCells(); i2++){
238	AliHLTTPCCACell *c2 = &(row2.Cells()[i2]);
239	Double_t sy2 = c2->ErrY()*c2->ErrY();
240	Double_t sz2 = c2->ErrZ()*c2->ErrZ();
241	Double_t dist1 = sqrt( (c1->Y()-c2->Y())*(c1->Y()-c2->Y())/(sy1+sy2) );
242	Double_t dist2 = sqrt( (c1->Z()-c2->Z())*(c1->Z()-c2->Z())/(sz1+sz2) );
243	if( dist1>3.5*fParam.CellConnectionFactor() ) continue;
244	if( dist2>3.5*fParam.CellConnectionFactor() ) continue;
245	if( c1->IUp() ==-1 ) c1->IUp() = (i2<<8)+iRow2;
246	else c1->IUp() = -2;
247	if( c2->IDown() ==-1 ) c2->IDown() = (i1<<8)+iRow1;
248	else c2->IDown() = -2;
249	}
250	}
251	}
252	}
253
254	Int_t nOutTrackHits = 0;
255	Int_t nTrackCells = 0;
256	for( Int_t iRow1=0; iRow1<fParam.NRows(); iRow1++ ){
257	AliHLTTPCCARow &row1 = fRows[iRow1];
258	for (Int_t i1 = 0; i1<row1.NCells(); i1++){
259	AliHLTTPCCACell *c1 = &(row1.Cells()[i1]);
260	if( c1->IDown()==-2 \|\| c1->IUp()==-2 ) continue;
261	if( c1->IUsed()>0 ) continue;
262	c1->IUsed() = 1;
263	AliHLTTPCCATrack track;
264	track.Used() = 0;
265	track.NCells() = 1;
266	track.IFirstCell() = nTrackCells;
267	fTrackCells[nTrackCells++] = (i1<<8)+iRow1;
268	AliHLTTPCCACell *last = c1;
269	Int_t lastRow = iRow1;
270	while( last->IUp() >=0 ){
271	Int_t iRow2 = last->IUp()%256;
272	AliHLTTPCCARow &row2 = fRows[iRow2];
273	AliHLTTPCCACell *next = &(row2.Cells()[last->IUp()>>8]);
274	if( next->IDown()==-2 \|\| next->IUp()==-2 ) break;
275	#ifdef DRAW
276	AliHLTTPCCADisplay::Instance().ConnectCells( lastRow,last,iRow2,next );
277	nConnectedCells++;
278	#endif
279	next->IUsed() = 1;
280	fTrackCells[nTrackCells++] = last->IUp();
281	track.NCells()++;
282	last = next;
283	lastRow = iRow2;
284	}
285	vTracks.push_back(track);
286	}
287	}
288
289	Int_t nTracks = vTracks.size();
290	std::sort( vTracks.begin(), vTracks.end(), AliHLTTPCCATrack::CompareSize);
291
292	fTracks = new AliHLTTPCCATrack[nTracks];
293	fNTracks = 0;
294	vTrackCells.clear();
295
ce622827	296	Int_t * vMatchedTracks = new Int_t[nTracks];
326c2d4b	297	for( Int_t itr=0; itr<nTracks; itr++ ){
	298	AliHLTTPCCATrack &iTrack = vTracks[itr];
	299	if( iTrack.Used() ) continue;
	300	FitTrack( iTrack );
	301	if( iTrack.Param().Chi2() > fParam.TrackChi2Cut()*iTrack.Param().NDF() ) continue;
	302
	303	Int_t iFirstRow = GetTrackCellIRow( iTrack, 0);
	304	Int_t iLastRow = GetTrackCellIRow( iTrack, iTrack.NCells()-1);
	305	AliHLTTPCCACell *iFirstCell = &GetTrackCell( iTrack, 0);
	306	AliHLTTPCCACell *iLastCell = &GetTrackCell( iTrack, iTrack.NCells()-1);
	307
	308	Bool_t updated = 1;
	309	Int_t nMatched = 0;
	310	std::vector<Int_t> vMatchedCells;
	311	while( updated ){
	312	updated = 0;
	313	for( Int_t jtr=1; jtr<nTracks; jtr++ ){
	314	AliHLTTPCCATrack &jTrack = vTracks[jtr];
	315	if( jTrack.Used() ) continue;
	316	Int_t jFirstRow = GetTrackCellIRow( jTrack, 0);
	317	Int_t jLastRow = GetTrackCellIRow( jTrack, jTrack.NCells()-1);
	318	AliHLTTPCCACell *jFirstCell = &GetTrackCell( jTrack, 0);
	319	AliHLTTPCCACell *jLastCell = &GetTrackCell( jTrack, jTrack.NCells()-1);
	320
	321	Int_t dFirstRow1 = TMath::Abs(iFirstRow-jLastRow);
	322	Int_t dFirstRow2 = TMath::Abs(iFirstRow-jFirstRow);
	323	Int_t dLastRow1 = TMath::Abs(iLastRow-jLastRow);
	324	Int_t dLastRow2 = TMath::Abs(iLastRow-jFirstRow);
	325	if( dFirstRow1 > fParam.MaxTrackMatchDRow() &&
	326	dFirstRow2 > fParam.MaxTrackMatchDRow() &&
	327	dLastRow1 > fParam.MaxTrackMatchDRow() &&
	328	dLastRow2 > fParam.MaxTrackMatchDRow() ) continue;
	329	Int_t iCase=0;
	330	AliHLTTPCCACell iC, jC;
	331	if( dFirstRow1<dFirstRow2 && dFirstRow1<dLastRow1 && dFirstRow1<dLastRow2 ){
	332	iCase = 0;
	333	iC = iFirstCell;
	334	jC = jLastCell;
	335	}else if( dFirstRow2<dLastRow1 && dFirstRow2<dLastRow2 ){
	336	iCase = 1;
	337	iC = iFirstCell;
	338	jC = jFirstCell;
	339	}else if( dLastRow1<dLastRow2 ){
	340	iCase = 2;
	341	iC = iLastCell;
	342	jC = jLastCell;
	343	}else{
	344	iCase = 3;
	345	iC = iLastCell;
	346	jC = jFirstCell;
	347	}
	348	{
	349	Double_t dy = TMath::Abs(iC->Y() - jC->Y());
	350	Double_t dz = TMath::Abs(iC->Z() - jC->Z());
	351	Double_t sy1 = iC->ErrY()*iC->ErrY();
	352	Double_t sz1 = iC->ErrZ()*iC->ErrZ();
	353	Double_t sy2 = jC->ErrY()*jC->ErrY();
	354	Double_t sz2 = jC->ErrZ()*jC->ErrZ();
	355	Double_t dist1 = sqrt( (dy)/(sy1+sy2) );
	356	Double_t dist2 = sqrt( (dz)/(sz1+sz2) );
	357	if( dist1>3.55fParam.CellConnectionFactor() ) continue;
	358	if( dist2>3.55fParam.CellConnectionFactor() ) continue;
	359	}
	360	AliHLTTPCCATrackPar t = iTrack.Param();
361	//t.Chi2() = 0;
362	//t.NDF() = 0;
363	for( Int_t i=0; i<jTrack.NCells() ; i++){
364	AliHLTTPCCACell &c = GetTrackCell(jTrack,i);
365	AliHLTTPCCARow &row = GetTrackCellRow(jTrack,i);
366	for( Int_t j=0; j<c.NHits(); j++){
367	AliHLTTPCCAHit &h = row.GetCellHit(c,j);
368	Double_t m[3] = {row.X(), h.Y(), h.Z() };
369	Double_t mV[6] = {fParam.ErrX()fParam.ErrX(), 0, h.ErrY()h.ErrY(), 0, 0, h.ErrZ()*h.ErrZ() };
370	Double_t mV1[6];
371	t.TransportBz(fParam.Bz(),m);
372	t.GetConnectionMatrix(fParam.Bz(),m, mV1);
373	t.Filter(m, mV, mV1);
374	}
375	}
376	if( t.Chi2() > fParam.TrackChi2Cut()*t.NDF() ) continue;
377	if( iCase==0 ){
378	iFirstRow = jFirstRow;
379	iFirstCell = jFirstCell;
380	}else if( iCase ==1 ){
381	iFirstRow = jLastRow;
382	iFirstCell = jLastCell;
383	}else if( iCase == 2 ){
384	iLastRow = jFirstRow;
385	iLastCell = jFirstCell;
386	}else{
387	iLastRow = jLastRow;
388	iLastCell = jLastCell;
389	}
390	t.Normalize();
391
392	for( Int_t i=0; i<jTrack.NCells(); i++){
393	vMatchedCells.push_back(fTrackCells[jTrack.IFirstCell()+i]);
394	}
395	//t.NDF()+= iTrack.Param().NDF();
396	//t.Chi2()+= iTrack.Param().Chi2();
397	iTrack.Param() = t;
398	vMatchedTracks[nMatched++] = jtr;
399	jTrack.Used()=1;
400	updated = 1;
401	break;
402	}
403	}
404
405	if(0){
406	Double_t t0[7];
407	for( Int_t i=0; i<7; i++ ) t0[i] = iTrack.Param().Par()[i];
408	iTrack.Param().Init();
409	for( Int_t i=0; i<7; i++ ) iTrack.Param().Par()[i] = t0[i];
410	for( Int_t i=0; i<iTrack.NCells() ; i++){
411	AliHLTTPCCACell &c = GetTrackCell( iTrack, i);
412	AliHLTTPCCARow &row = GetTrackCellRow( iTrack, i);
413	for( Int_t j=0; j<c.NHits(); j++){
414	AliHLTTPCCAHit &h = row.GetCellHit(c,j);
415	Double_t m[3] = {row.X(), h.Y(), h.Z() };
416	Double_t mV[6] = {0, 0, h.ErrY()h.ErrY(), 0, 0, h.ErrZ()h.ErrZ() };
417	Double_t mV1[6];
418	iTrack.Param().TransportBz(fParam.Bz(), m, t0);
419	iTrack.Param().GetConnectionMatrix(fParam.Bz(), m, mV1, t0);
420	iTrack.Param().Filter(m, mV, mV1);
421	}
422	}
423	iTrack.Param().Normalize();
424	}
425	//FitTrack(iTrack,5);
426	Int_t nHits = 0;
427	for( Int_t iCell=0; iCell<iTrack.NCells(); iCell++){
428	Int_t ind = fTrackCells[iTrack.IFirstCell()+iCell];
429	AliHLTTPCCARow &row = fRows[ind%256];
430	AliHLTTPCCACell &c = row.Cells()[ind>>8];
431	nHits+=c.NHits();
432	}
433	for( UInt_t i=0; i<vMatchedCells.size(); i++){
434	Int_t ind = vMatchedCells[i];
435	AliHLTTPCCARow &row = fRows[ind%256];
436	AliHLTTPCCACell &c = row.Cells()[ind>>8];
437	nHits+=c.NHits();
438	}
439
440	if( nHits<5 ){
441	for( Int_t i=0; i<nMatched; i++ ) vTracks[vMatchedTracks[i]].Used()=0;
442	continue;
443	}
444	iTrack.Used() = 1;
445	Int_t oldNCells = vTrackCells.size();
446	for( Int_t i=0; i<iTrack.NCells() ; i++){
447	vTrackCells.push_back(fTrackCells[iTrack.IFirstCell()+i]);
448	}
449	iTrack.IFirstCell() = oldNCells;
450	for( UInt_t i=0; i<vMatchedCells.size(); i++){
451	vTrackCells.push_back(vMatchedCells[i]);
452	iTrack.NCells()++;
453	}
454	fTracks[fNTracks++] = iTrack;
455	nOutTrackHits+= nHits;
456	}
ce622827	457	delete [] vMatchedTracks;
326c2d4b	458
	459	//fTrackCells = new Int_t[vTrackCells.size()];
	460	for( UInt_t i=0; i<vTrackCells.size(); i++ ) fTrackCells[i] = vTrackCells[i];
	461
	462	#ifdef DRAW
	463	if( nConnectedCells>0 ) AliHLTTPCCADisplay::Instance().Ask();
	464	#endif
	465
	466
	467	fOutTrackHits = new Int_t[nOutTrackHits];
	468	fNOutTrackHits = 0;
	469	fNOutTracks = fNTracks;
	470	fOutTracks = new AliHLTTPCCAOutTrack[fNOutTracks];
	471	for( Int_t itr=0; itr<fNOutTracks; itr++ ){
	472	AliHLTTPCCATrack &t = fTracks[itr];
	473	#ifdef DRAW
	474	AliHLTTPCCADisplay::Instance().DrawTrack( t );
	475	#endif
	476	AliHLTTPCCAOutTrack &tmp = fOutTracks[itr];
	477	tmp.FirstHitRef() = fNOutTrackHits;
	478	tmp.NHits() = 0;
	479	tmp.Param() = t.Param();
	480	for( Int_t iCell=0; iCell<t.NCells(); iCell++){
	481	AliHLTTPCCACell &cell = GetTrackCell(t,iCell);
	482	AliHLTTPCCARow &row = GetTrackCellRow(t,iCell);
	483	for( Int_t iHit=0; iHit<cell.NHits(); iHit++ ){
	484	AliHLTTPCCAHit &hit = row.GetCellHit(cell,iHit);
	485	fOutTrackHits[fNOutTrackHits] = hit.ID();
	486	fNOutTrackHits++;
	487	tmp.NHits()++;
	488	}
	489	}
	490	}
	491
	492	#ifdef DRAW
	493	AliHLTTPCCADisplay::Instance().Ask();
	494	//AliHLTTPCCADisplay::Instance().DrawMCTracks(fParam.fISec);
	495	//AliHLTTPCCADisplay::Instance().Update();
	496	//AliHLTTPCCADisplay::Instance().Ask();
	497	#endif
	498	}
	499
	500
	501
	502	void AliHLTTPCCATracker::FitTrack( AliHLTTPCCATrack &track, Int_t nIter )
	503	{
	504	// fit the track with nIter iterations
	505
	506	AliHLTTPCCATrackPar &t = track.Param();
	507	t.Init();
	508
	509	AliHLTTPCCACell &c1 = GetTrackCell(track,0);
	510	AliHLTTPCCACell &c2 = GetTrackCell(track,track.NCells()-1);
	511	AliHLTTPCCARow &row1 = GetTrackCellRow(track,0);
	512	AliHLTTPCCARow &row2 = GetTrackCellRow(track,track.NCells()-1);
	513	Double_t t0[7];
	514	t0[0]=row1.X();
	515	t0[1]=c1.Y();
	516	t0[2]=c1.Z();
	517	t0[3]= row2.X() - row1.X();
	518	t0[4]= c2.Y() - c1.Y();
	519	t0[5]= c2.Z() - c1.Z();
	520	Double_t tt = sqrt(t0[3]t0[3]+t0[4]t0[4]+t0[5]*t0[5]);
dd94bda2	521	if( TMath::Abs(tt)>1.e-4 ){
326c2d4b	522	t0[3]/=tt;
	523	t0[4]/=tt;
	524	t0[5]/=tt;
	525	}else{
	526	t0[4]=1;
	527	}
	528	t0[6] = 0;
	529
	530	for( Int_t iter=0; iter<nIter; iter++ ){
	531	t.Init();
	532	for( Int_t i=0; i<7; i++) t.Par()[i] = t0[i];
	533	{
	534	Double_t m[3] = {row1.X(), c1.Y(), c1.Z() };
	535	t.TransportBz(fParam.Bz(),m,t0);
	536	}
	537	t.Init();
	538	for( Int_t i=0; i<7; i++ ) t.Par()[i] = t0[i];
	539
	540	for( Int_t i=0; i<track.NCells() ; i++){
	541	AliHLTTPCCACell &c = GetTrackCell(track,i);
	542	AliHLTTPCCARow &row = GetTrackCellRow(track,i);
	543	for( Int_t j=0; j<c.NHits(); j++){
	544	AliHLTTPCCAHit &h = row.GetCellHit(c,j);
	545	Double_t m[3] = {row.X(), h.Y(), h.Z() };
	546	Double_t mV[6] = {fParam.ErrX()fParam.ErrX(), 0, h.ErrY()h.ErrY(), 0, 0, h.ErrZ()*h.ErrZ() };
	547	Double_t mV1[6];
	548	t.TransportBz(fParam.Bz(),m, t0);
	549	t.GetConnectionMatrix(fParam.Bz(),m, mV1, t0);
	550	t.Filter(m, mV, mV1);
	551	}
	552	}
	553	t.Normalize();
	554	for( Int_t i=0; i<7; i++ ) t0[i] = t.Par()[i];
	555	}
	556	}