[u/mrichter/AliRoot.git] / HLT / global / AliHLTGlobalTrackMerger.cxx

//$Id$

/**************************************************************************
 * This file is property of and copyright by the ALICE HLT Project        * 
 * ALICE Experiment at CERN, All rights reserved.                         *
 *                                                                        *
 * Primary Authors: Jacek Otwinowski (Jacek.Otwinowski@gsi.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.                  *
 **************************************************************************/

//  @file   AliHLTGlobalTrackMerger.cxx
//  @author Jacek Otwinowski
//  @date   
//  @brief  The HLT global merger base class
// 

//#include "AliTPCReconstructor.h"

#include "AliESDEvent.h"
#include "AliESDVertex.h"
#include "AliESDtrack.h"
#include "AliTracker.h"
#include "TTreeStream.h"

#include <TClonesArray.h>

#include "AliHLTGlobalTrackMerger.h"

#if __GNUC__>= 3
using namespace std;
#endif

ClassImp(AliHLTGlobalTrackMerger)

AliHLTGlobalTrackMerger::AliHLTGlobalTrackMerger()
  :
  fMaxY(0.0),
  fMaxZ(0.0),
  fMaxSnp(0.0),
  fMaxTgl(0.0),
  fMaxSigned1Pt(0.0),
  fVertex(0),
  fDebugStreamer(0)
{
  //Default constructor

  // standard vertex settings at the moment
  // V(0.,0.,0.), sigmaVx=sigmaVy=5.e-3 [cm], sigmaVz=5.3 [cm]    
  fVertex = new AliESDVertex;

  // if (AliTPCReconstructor::StreamLevel()>0) {
  //   fDebugStreamer = new TTreeSRedirector("debugGlobalMerger.root");
  // }

}

//_____________________________________________________________________________
AliHLTGlobalTrackMerger::~AliHLTGlobalTrackMerger()
{
  //Destructor
  if(fVertex) delete fVertex; fVertex =0;
  if(fDebugStreamer) delete fDebugStreamer; fDebugStreamer =0;
}

//_____________________________________________________________________________
Bool_t AliHLTGlobalTrackMerger::Merge(AliESDEvent* esdEvent)
{
  // merge TPC and TRD tracks
  // 1. propagate TPC track to the radius between TPC and TRD
  // 2. propagate TRD track to the same radius between TPC and TRD
  // 3. matches TPC and TRD tracks at the radius
  // 4. propagate matched TRD track to the merging radius (first measured TPC point - x coordinate)
  // 5. merge TPC and TRD track parameters at the merging radius 
  // 6. create AliESDtrack from merged tracks
  // 7. add AliESDtrack to AliESDEvent

  if(!esdEvent) return kFALSE;

  const Double_t kMaxStep     = 10.0;    // [cm] track propagation step
  const Double_t kMatchRadius = 285.0;   // [cm] matching at radius between TPC and TRD
  Double_t kMergeRadius = 0.0;
  Bool_t isOk = kFALSE;
  Bool_t isMatched = kFALSE;
  AliESDtrack *track=0;
  AliExternalTrackParam *extTPCTrack = 0;

  Int_t nTracks = esdEvent->GetNumberOfTracks(); 
  HLTWarning("nTracks %d",nTracks);

  Int_t nTracksTPC =0;
  for(Int_t iTrack = 0; iTrack<nTracks; ++iTrack) 
  {
    track = esdEvent->GetTrack(iTrack); 
    if(!track) continue;

    // TPC tracks
    if((track->GetStatus()&AliESDtrack::kTPCin)==0) continue;
    AliESDtrack *tpcTrack = track;
    nTracksTPC++;

    // create external tpc track param (needed to propagate to matching radius) 
    if ((extTPCTrack = new AliExternalTrackParam(*tpcTrack)) == 0) continue; 

    kMergeRadius = tpcTrack->GetTPCPoints(0); // [cm] merging at first measured TPC point

    HLTInfo("-------------------------------------------------------------------------------------");
    //HLTInfo("-----tpc track before matching: alpha %f, x %f, y, %f, z %f, snp %f, tgl %f, 1pt %f",tpcTrack->GetAlpha(),tpcTrack->GetX(),tpcTrack->GetY(),tpcTrack->GetZ(),tpcTrack->GetSnp(),tpcTrack->GetTgl(),tpcTrack->GetSigned1Pt());
    HLTInfo("-----ext tpc track before matching: alpha %f, x %f, y, %f, z %f, snp %f, tgl %f, 1pt %f",extTPCTrack->GetAlpha(),extTPCTrack->GetX(),extTPCTrack->GetY(),extTPCTrack->GetZ(),extTPCTrack->GetSnp(),extTPCTrack->GetTgl(),extTPCTrack->GetSigned1Pt());

    // propagate tracks to the matching radius 
    //isOk = AliTracker::PropagateTrackTo(tpcTrack,kMatchRadius,tpcTrack->GetMass(),kMaxStep,kFALSE);
    isOk = AliTracker::PropagateTrackTo(extTPCTrack,kMatchRadius,tpcTrack->GetMass(),kMaxStep,kFALSE);
    if(!isOk) continue;

    HLTInfo("-----ext tpc track at matching radius: alpha %f, x %f, y, %f, z %f, snp %f, tgl %f, 1pt %f",extTPCTrack->GetAlpha(),extTPCTrack->GetX(),extTPCTrack->GetY(),extTPCTrack->GetZ(),extTPCTrack->GetSnp(),extTPCTrack->GetTgl(),extTPCTrack->GetSigned1Pt());

    Int_t nTracksTRD =0;
    for(Int_t jTrack = 0; jTrack<nTracks; ++jTrack) 
    {
      track = esdEvent->GetTrack(jTrack); 
      if(!track) continue;

      // TRD tracks
      if((track->GetStatus()&AliESDtrack::kTRDin)==0) continue;
      AliESDtrack *trdTrack = track;
      nTracksTRD++;

      HLTInfo("-----trd track before matching: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",trdTrack->GetAlpha(),trdTrack->GetX(),trdTrack->GetY(),trdTrack->GetZ(),trdTrack->GetSnp(),trdTrack->GetTgl(),trdTrack->GetSigned1Pt());

      isOk = AliTracker::PropagateTrackTo(trdTrack,kMatchRadius,trdTrack->GetMass(),kMaxStep,kFALSE);
      if(!isOk) continue;

      HLTInfo("-----trd track at matching radius: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",trdTrack->GetAlpha(),trdTrack->GetX(),trdTrack->GetY(),trdTrack->GetZ(),trdTrack->GetSnp(),trdTrack->GetTgl(),trdTrack->GetSigned1Pt());

      // match TPC and TRD tracks
      //isMatched = MatchTracks(tpcTrack,trdTrack);
      isMatched = MatchTracks(extTPCTrack,trdTrack);
      if(!isMatched) continue;

      //HLTInfo("-----tpc track after matching: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",tpcTrack->GetAlpha(),tpcTrack->GetX(),tpcTrack->GetY(),tpcTrack->GetZ(),tpcTrack->GetSnp(),tpcTrack->GetTgl(),tpcTrack->GetSigned1Pt());
      HLTInfo("-----ext tpc track after matching: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",extTPCTrack->GetAlpha(),extTPCTrack->GetX(),extTPCTrack->GetY(),extTPCTrack->GetZ(),extTPCTrack->GetSnp(),extTPCTrack->GetTgl(),extTPCTrack->GetSigned1Pt());
      HLTInfo("-----trd track after matching: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",trdTrack->GetAlpha(),trdTrack->GetX(),trdTrack->GetY(),trdTrack->GetZ(),trdTrack->GetSnp(),trdTrack->GetTgl(),trdTrack->GetSigned1Pt());

      // propagate TRD tracks to the merging radius 
      //isOk = AliTracker::PropagateTrackTo(tpcTrack,kMergeRadius,tpcTrack->GetMass(),kMaxStep,kFALSE);
      //if(!isOk) continue;
      isOk = AliTracker::PropagateTrackTo(trdTrack,kMergeRadius,trdTrack->GetMass(),kMaxStep,kFALSE);
      if(!isOk) continue;

      HLTInfo("-----tpc before merging: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",tpcTrack->GetAlpha(),tpcTrack->GetX(),tpcTrack->GetY(),tpcTrack->GetZ(),tpcTrack->GetSnp(),tpcTrack->GetTgl(),tpcTrack->GetSigned1Pt());
      HLTInfo("-----trd before merging: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",trdTrack->GetAlpha(),trdTrack->GetX(),trdTrack->GetY(),trdTrack->GetZ(),trdTrack->GetSnp(),trdTrack->GetTgl(),trdTrack->GetSigned1Pt());

      // merge TPC and TRD tracks
      // create AliESDtrack and add it to AliESDevent
      Bool_t isMerged = MergeTracks(tpcTrack,trdTrack,esdEvent);
      if(!isMerged) 
        HLTInfo("No merged tracks");
    }
    HLTInfo("nTracksTRD %d",nTracksTRD);

    // delete external TPC track
    if(extTPCTrack) delete extTPCTrack; extTPCTrack=0;
  }
  HLTInfo("nTracksTPC %d",nTracksTPC);

return kTRUE;
}

//_____________________________________________________________________________
Bool_t AliHLTGlobalTrackMerger::MergeTracks(AliESDtrack *tpcTrack, AliESDtrack* trdTrack, AliESDEvent *esdEvent)
{
  // merge TPC and TRD track parameters
  // create new AliESDtrack with TPC+TRD merged track parameters
  // add AliESDtrack to AliESDEvent

  if(!tpcTrack) return kFALSE;
  if(!trdTrack) return kFALSE;

  Double_t trackParam[5], trackCovar[15]; 
  Double_t trackChi2;
  Int_t trackNDF;

  // calculate merged track parameters
  Bool_t isNotOK = SmoothTracks(tpcTrack->GetParameter(), tpcTrack->GetCovariance(), tpcTrack->GetTPCchi2(), 5,
                                trdTrack->GetParameter(), trdTrack->GetCovariance(), trdTrack->GetTRDchi2(), 5,
                                trackParam, trackCovar, trackChi2, trackNDF,5);

  if(isNotOK) 
      return kFALSE;

  //
  // create AliESDtrack
  // merged TPC+TRD information
  // 

  AliESDtrack track;
  //track.UpdateTrackParams(tpcTrack, AliESDtrack::kTPCrefit);
  track.SetStatus(AliESDtrack::kGlobalMerge);
  track.SetLabel(tpcTrack->GetLabel());
  track.Set(tpcTrack->GetX(),tpcTrack->GetAlpha(),trackParam,trackCovar);
  track.SetGlobalChi2(trackChi2);

  //track.SetTPCLabel(tpcTrack->GetLabel());
  Double32_t tpcPID[AliPID::kSPECIES];
  tpcTrack->GetTPCpid(tpcPID);
  track.SetTPCpid(tpcPID);
  //fTPCncls=t->GetNumberOfClusters();  // no cluster on HLT
  //fTPCchi2=t->GetChi2();

  //track.SetTRDLabel(trdTrack->GetLabel());
  Double32_t trdPID[AliPID::kSPECIES];
  trdTrack->GetTRDpid(trdPID);
  track.SetTRDpid(trdPID);
  //fTRDchi2  = t->GetChi2();
  //fTRDncls  = t->GetNumberOfClusters();
  //for (Int_t i=0;i<6;i++) index[i]=t->GetTrackletIndex(i);

  // add track to AliESDEvent
  esdEvent->AddTrack(&track);

return kTRUE;
}

//_____________________________________________________________________________
void AliHLTGlobalTrackMerger::SetParameter(Double_t maxy, Double_t maxz, Double_t maxsnp, Double_t maxtgl, Double_t signed1Pt)
{ 
  //set parameters for merger
  fMaxY = maxy;
  fMaxZ = maxz;
  fMaxSnp = maxsnp;
  fMaxTgl = maxtgl;
  fMaxSigned1Pt = signed1Pt;
}

//_____________________________________________________________________________
Bool_t AliHLTGlobalTrackMerger::MatchTracks(AliExternalTrackParam *trackTPC, const AliESDtrack *trackTRD)
{ 
  // match TPC and TRD tracks 
  // return kTRUE in case of matching
 
  if(!trackTPC) return kFALSE;
  if(!trackTRD) return kFALSE;

  Double_t xTpc=trackTPC->GetX();
  Double_t yTpc=trackTPC->GetY();
  Double_t zTpc=trackTPC->GetZ();
  Double_t snpTpc=trackTPC->GetSnp();
  Double_t tglTpc=trackTPC->GetTgl();
  Double_t signed1PtTpc=trackTPC->GetSigned1Pt();

  Double_t xTrd=trackTRD->GetX();
  Double_t yTrd=trackTRD->GetY();
  Double_t zTrd=trackTRD->GetZ();
  Double_t snpTrd=trackTRD->GetSnp();
  Double_t tglTrd=trackTRD->GetTgl();
  Double_t signed1PtTrd=trackTRD->GetSigned1Pt();

  // debug stream
  // if (AliTPCReconstructor::StreamLevel()>0) {
  // //TTreeSRedirector &cstream = *fDebugStreamer;
  // *fDebugStreamer<<"match"<<
  // "xTpc="<<xTpc<<
  // "yTpc="<<yTpc<<
  // "zTpc="<<zTpc<<
  // "snpTpc="<<snpTpc<<
  // "tglTpc="<<tglTpc<<
  // "signed1PtTpc="<<signed1PtTpc<<
  // "xTrd="<<xTrd<<
  // "yTrd="<<yTrd<<
  // "zTrd="<<zTrd<<
  // "snpTrd="<<snpTrd<<
  // "tglTrd="<<tglTrd<<
  // "signed1PtTrd="<<signed1PtTrd<<
  // "\n";
  // }

  if (TMath::Abs(xTpc-xTrd) > 0) {/* get rid of warning*/;}
  if (TMath::Abs(yTpc-yTrd) > fMaxY) return kFALSE;
  if (TMath::Abs(zTpc-zTrd) > fMaxZ) return kFALSE;
  if (TMath::Abs(snpTpc-snpTrd) > fMaxSnp) return kFALSE;
  if (TMath::Abs(tglTpc-tglTrd) > fMaxTgl) return kFALSE;
  if (TMath::Abs(signed1PtTpc-signed1PtTrd) > fMaxSigned1Pt) return kFALSE;

return kTRUE;
}

//_____________________________________________________________________________
Bool_t AliHLTGlobalTrackMerger::SmoothTracks( const Double_t T1[], const Double_t C1[], Double_t Chi21, Int_t NDF1,
				    const  Double_t T2[], const Double_t C2[], Double_t Chi22, Int_t NDF2,
				     Double_t T [], Double_t C [], Double_t &Chi2, Int_t &NDF,
				     Int_t N  )
{
  //* Smooth two tracks with parameter vectors of size N
  //*
  //* Input:
  //*
  //* T1[N], T2[N] - tracks
  //* C1[N*(N+1)/2], C2[N*(N+1)/2] - covariance matrices in low-diagonal form:
  //* C = { c00, 
  //*       c10, c11, 
  //*       c20, c21, c22, 
  //*       ...             };
  //* Chi2{1,2}, NDF{1,2} - \Chi^2 and "Number of Degrees of Freedom" values for both tracks
  //* Output: 
  //*
  //* T[N], C[N] ( can be aqual to {T1,C1}, or {T2,C2} )
  //* Chi2, NDF
  //*
  //* returns error flag (0 means OK, 1 not OK )
    
  Int_t M = N*(N+1)/2;
  
  Double_t A[M];
  Double_t K[N*N];
  
  for(Int_t k=0; k<M; k++) A[k] = C1[k] + C2[k]; 
  Bool_t err = InvertS(A,N);
  if( err ) return 1;

  Chi2 = Chi21 + Chi22;
  NDF = NDF1 + NDF2;
  
  MultSSQ( C1, A, K, N);	
  Double_t r[N];
  for( Int_t k=0; k<N;k++) r[k] = T1[k] - T2[k]; 
  for( Int_t k=0; k<N;k++ )
    for( Int_t l=0;l<N;l++) T[k] = T1[k] - K[k*N+l]*r[l];

  for( Int_t ind=0,i=0; i<N; i++ ){
    for( Int_t j=0; j<i; j++ ) Chi2+= 2*r[i]*r[j]*A[ind++];
    Chi2+= r[i]*r[i]*A[ind++];
  }
  NDF+=N;

  for( Int_t l=0; l<N; l++ ) K[ (N+1)*l ] -= 1;
  
  for( Int_t ind = 0, l=0; l<N; ++l ){
    for( Int_t j=0; j<=l; ++j, ind++ ){
      A[ind] = 0;
      for( Int_t k=0; k<N; ++k ) A[ind] -= K[l*N+k] * C1[IndexS(j,k)];
    }
  }
  for( Int_t l=0; l<N; l++ ) C[l] = A[l];
  return 0;
}

//_____________________________________________________________________________
void AliHLTGlobalTrackMerger::MultSSQ( const Double_t *A, const Double_t *B, Double_t *C, Int_t N )
{
  // no clue
  for( Int_t ind=0, i=0; i<N; ++i ){
    for( Int_t j=0; j<N; ++j, ++ind ){
      C[ind] = 0;
      for( Int_t k=0; k<N; ++k ) C[ind] += A[IndexS(i,k)] * B[IndexS(k,j)];
    }
  }
}

//_____________________________________________________________________________
Bool_t AliHLTGlobalTrackMerger::InvertS( Double_t A[], Int_t N )
{
  //* input: simmetric > 0 NxN matrix A = {a11,a21,a22,a31..a33,..}  
  //* output: inverse A, in case of problems fill zero and return 1
  //*  
  //* A->low triangular Anew : A = Anew x Anew^T
  //* method:
  //* for(j=1,N) for(i=j,N) Aij=(Aii-sum_{k=1}^{j-1}Aik*Ajk )/Ajj
  //*   

  Bool_t ret = 0;
  
  const Double_t ZERO = 1.E-20;
    
  {
    Double_t *j1 = A, *jj = A;
    for( Int_t j=1; j<=N; j1+=j++, jj+=j ){
      Double_t *ik = j1, x = 0;
      while( ik!=jj ){
	x -= (*ik) * (*ik);
	ik++;
      }
      x += *ik;
      if( x > ZERO ){
	x = sqrt(x);
	*ik = x;
	ik++;
	x = 1 / x;
	for( Int_t step=1; step<=N-j; ik+=++step ){ // ik==Ai1
	  Double_t sum = 0;
	  for( Double_t *jk=j1; jk!=jj; sum += (*(jk++)) * (*(ik++)) ) {}
	  *ik = (*ik - sum) * x; // ik == Aij
	}
      }else{
	Double_t *ji=jj;
	for( Int_t i=j; i<N; i++ ) *(ji+=i) = 0.;
	ret = 1;
      }   
    }
  }
  
  //* A -> Ainv
  //* method : 
  //* for(i=1,N){ 
  //*   Aii = 1/Aii; 
  //*   for(j=1,i-1) Aij=-(sum_{k=j}^{i-1} Aik * Akj) / Aii ;
  //* }
  
  {
    Double_t *ii=A,*ij=A;
    for( Int_t i = 1; i<=N; ij=ii+1, ii+=++i ){
      if( *ii > ZERO ){
	Double_t x = -(*ii = 1./ *ii);
	{ 
	  Double_t *jj = A;
	  for( Int_t j=1; j<i; jj+=++j, ij++ ){
	    Double_t *ik = ij, *kj = jj, sum = 0.;
	    for( Int_t k=j; ik!=ii; kj+=k++, ik++ ){
	      sum += *ik * *kj;
	    }
	    *kj = sum * x;
	  }
	}
      }else{      
	for( Double_t *ik = ij; ik!=ii+1; ik++ ){
	  *ik = 0.;
	}
	ret = 1;
      }
    }
  }
  
  //* A -> A^T x A
  //* method: 
  //* Aij = sum_{k=i}^N Aki * Akj
  
  {
    Double_t *ii=A, *ij=A;
    for( Int_t i=1; i<=N; ii+=++i ){
      do{ 
	Double_t *ki = ii, *kj = ij, sum = 0.;
	for( Int_t k=i; k<=N; ki+=k, kj+=k++ ) sum += (*ki) * (*kj);
	*ij = sum;
      }while( (ij++)!=ii );
    }    
  }
  return ret;    
}

//_____________________________________________________________________________
void  AliHLTGlobalTrackMerger::PropagateTracksToDCA(const AliESDEvent *esdEvent)
{
  // try to propagate all tracks to DCA to primary vertex
  if(!esdEvent) return;

  const Double_t kBz = esdEvent->GetMagneticField();
  const Double_t kSmallRadius  = 2.8; // [cm] something less than the beam pipe radius
  const Double_t kMaxStep  = 10.0;    // [cm] track propagation step
  Bool_t isOK = kFALSE;

  Int_t nTracks = esdEvent->GetNumberOfTracks(); 
  for(Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
    AliESDtrack *track = esdEvent->GetTrack(iTrack); 
    if(!track) continue;

    // propagate to small radius (material budget included)
    isOK = AliTracker::PropagateTrackTo(track,kSmallRadius,track->GetMass(),kMaxStep,kFALSE);

    // relate tracks to DCA to primary vertex
    if(isOK) 
    {
      track->RelateToVertex(fVertex, kBz, kVeryBig);
      HLTInfo("1-------: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",track->GetAlpha(),track->GetX(),track->GetY(),track->GetZ(),track->GetSnp(),track->GetTgl(),track->GetSigned1Pt());
    }

    //
    // the same procedure must be repeated for TPCinner (TPC only) tracks
    //
    AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
    if(!tpcTrack) continue;

    // propagate tpcTracks to small radius (material budget included)
    isOK = AliTracker::PropagateTrackTo(tpcTrack,kSmallRadius,track->GetMass(),kMaxStep,kFALSE);

    // relate tracks to DCA to primary vertex
    if(isOK) 
    {
      Double_t par[2], cov[3];
      tpcTrack->PropagateToDCA(fVertex, kBz, kVeryBig,par,cov);
      HLTInfo("2-------: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",tpcTrack->GetAlpha(),tpcTrack->GetX(),tpcTrack->GetY(),tpcTrack->GetZ(),tpcTrack->GetSnp(),tpcTrack->GetTgl(),tpcTrack->GetSigned1Pt());
    }
  }
}
Commit	Line	Data
ec6160d5	1	//$Id$
	2
	3	/**************************************************************************
	4	* This file is property of and copyright by the ALICE HLT Project *
	5	* ALICE Experiment at CERN, All rights reserved. *
	6	* *
	7	* Primary Authors: Jacek Otwinowski (Jacek.Otwinowski@gsi.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
a0395f81	19	// @file AliHLTGlobalTrackMerger.cxx
	20	// @author Jacek Otwinowski
	21	// @date
	22	// @brief The HLT global merger base class
	23	//
ec6160d5	24
18ada816	25	//#include "AliTPCReconstructor.h"
ec6160d5	26
ec6160d5	27	#include "AliESDEvent.h"
ec6160d5	28	#include "AliESDVertex.h"
a0395f81	29	#include "AliESDtrack.h"
ec6160d5	30	#include "AliTracker.h"
a9f060e5	31	#include "TTreeStream.h"
ec6160d5	32
	33	#include <TClonesArray.h>
	34
	35	#include "AliHLTGlobalTrackMerger.h"
	36
	37	#if __GNUC__>= 3
	38	using namespace std;
	39	#endif
	40
	41	ClassImp(AliHLTGlobalTrackMerger)
	42
	43	AliHLTGlobalTrackMerger::AliHLTGlobalTrackMerger()
	44	:
	45	fMaxY(0.0),
	46	fMaxZ(0.0),
	47	fMaxSnp(0.0),
	48	fMaxTgl(0.0),
	49	fMaxSigned1Pt(0.0),
a9f060e5	50	fVertex(0),
a9f060e5	51	fDebugStreamer(0)
ec6160d5	52	{
	53	//Default constructor
	54
	55	// standard vertex settings at the moment
	56	// V(0.,0.,0.), sigmaVx=sigmaVy=5.e-3 [cm], sigmaVz=5.3 [cm]
	57	fVertex = new AliESDVertex;
a9f060e5	58
18ada816	59	// if (AliTPCReconstructor::StreamLevel()>0) {
	60	// fDebugStreamer = new TTreeSRedirector("debugGlobalMerger.root");
	61	// }
a9f060e5	62
ec6160d5	63	}
	64
	65	//_____________________________________________________________________________
	66	AliHLTGlobalTrackMerger::~AliHLTGlobalTrackMerger()
	67	{
	68	//Destructor
	69	if(fVertex) delete fVertex; fVertex =0;
a9f060e5	70	if(fDebugStreamer) delete fDebugStreamer; fDebugStreamer =0;
ec6160d5	71	}
ec6160d5	72
ec6160d5	73	//_____________________________________________________________________________
	74	Bool_t AliHLTGlobalTrackMerger::Merge(AliESDEvent* esdEvent)
	75	{
	76	// merge TPC and TRD tracks
	77	// 1. propagate TPC track to the radius between TPC and TRD
	78	// 2. propagate TRD track to the same radius between TPC and TRD
	79	// 3. matches TPC and TRD tracks at the radius
4b6d0cc6	80	// 4. propagate matched TRD track to the merging radius (first measured TPC point - x coordinate)
ec6160d5	81	// 5. merge TPC and TRD track parameters at the merging radius
	82	// 6. create AliESDtrack from merged tracks
	83	// 7. add AliESDtrack to AliESDEvent
	84
	85	if(!esdEvent) return kFALSE;
	86
	87	const Double_t kMaxStep = 10.0; // [cm] track propagation step
	88	const Double_t kMatchRadius = 285.0; // [cm] matching at radius between TPC and TRD
	89	Double_t kMergeRadius = 0.0;
	90	Bool_t isOk = kFALSE;
	91	Bool_t isMatched = kFALSE;
a9f060e5	92	AliESDtrack *track=0;
4b6d0cc6	93	AliExternalTrackParam *extTPCTrack = 0;
ec6160d5	94
ec6160d5	95	Int_t nTracks = esdEvent->GetNumberOfTracks();
ec6160d5	96	HLTWarning("nTracks %d",nTracks);
ec6160d5	97
a9f060e5	98	Int_t nTracksTPC =0;
ec6160d5	99	for(Int_t iTrack = 0; iTrack<nTracks; ++iTrack)
ec6160d5	100	{
a9f060e5	101	track = esdEvent->GetTrack(iTrack);
ec6160d5	102	if(!track) continue;
	103
	104	// TPC tracks
	105	if((track->GetStatus()&AliESDtrack::kTPCin)==0) continue;
	106	AliESDtrack *tpcTrack = track;
a9f060e5	107	nTracksTPC++;
ec6160d5	108
4b6d0cc6	109	// create external tpc track param (needed to propagate to matching radius)
	110	if ((extTPCTrack = new AliExternalTrackParam(*tpcTrack)) == 0) continue;
	111
a9f060e5	112	kMergeRadius = tpcTrack->GetTPCPoints(0); // [cm] merging at first measured TPC point
ec6160d5	113
a9f060e5	114	HLTInfo("-------------------------------------------------------------------------------------");
4b6d0cc6	115	//HLTInfo("-----tpc track before matching: alpha %f, x %f, y, %f, z %f, snp %f, tgl %f, 1pt %f",tpcTrack->GetAlpha(),tpcTrack->GetX(),tpcTrack->GetY(),tpcTrack->GetZ(),tpcTrack->GetSnp(),tpcTrack->GetTgl(),tpcTrack->GetSigned1Pt());
4b6d0cc6	116	HLTInfo("-----ext tpc track before matching: alpha %f, x %f, y, %f, z %f, snp %f, tgl %f, 1pt %f",extTPCTrack->GetAlpha(),extTPCTrack->GetX(),extTPCTrack->GetY(),extTPCTrack->GetZ(),extTPCTrack->GetSnp(),extTPCTrack->GetTgl(),extTPCTrack->GetSigned1Pt());
ec6160d5	117
ec6160d5	118	// propagate tracks to the matching radius
4b6d0cc6	119	//isOk = AliTracker::PropagateTrackTo(tpcTrack,kMatchRadius,tpcTrack->GetMass(),kMaxStep,kFALSE);
4b6d0cc6	120	isOk = AliTracker::PropagateTrackTo(extTPCTrack,kMatchRadius,tpcTrack->GetMass(),kMaxStep,kFALSE);
ec6160d5	121	if(!isOk) continue;
ec6160d5	122
4b6d0cc6	123	HLTInfo("-----ext tpc track at matching radius: alpha %f, x %f, y, %f, z %f, snp %f, tgl %f, 1pt %f",extTPCTrack->GetAlpha(),extTPCTrack->GetX(),extTPCTrack->GetY(),extTPCTrack->GetZ(),extTPCTrack->GetSnp(),extTPCTrack->GetTgl(),extTPCTrack->GetSigned1Pt());
ec6160d5	124
a9f060e5	125	Int_t nTracksTRD =0;
ec6160d5	126	for(Int_t jTrack = 0; jTrack<nTracks; ++jTrack)
ec6160d5	127	{
a9f060e5	128	track = esdEvent->GetTrack(jTrack);
ec6160d5	129	if(!track) continue;
	130
	131	// TRD tracks
	132	if((track->GetStatus()&AliESDtrack::kTRDin)==0) continue;
	133	AliESDtrack *trdTrack = track;
a9f060e5	134	nTracksTRD++;
a9f060e5	135
4b6d0cc6	136	HLTInfo("-----trd track before matching: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",trdTrack->GetAlpha(),trdTrack->GetX(),trdTrack->GetY(),trdTrack->GetZ(),trdTrack->GetSnp(),trdTrack->GetTgl(),trdTrack->GetSigned1Pt());
ec6160d5	137
	138	isOk = AliTracker::PropagateTrackTo(trdTrack,kMatchRadius,trdTrack->GetMass(),kMaxStep,kFALSE);
	139	if(!isOk) continue;
	140
4b6d0cc6	141	HLTInfo("-----trd track at matching radius: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",trdTrack->GetAlpha(),trdTrack->GetX(),trdTrack->GetY(),trdTrack->GetZ(),trdTrack->GetSnp(),trdTrack->GetTgl(),trdTrack->GetSigned1Pt());
ec6160d5	142
ec6160d5	143	// match TPC and TRD tracks
4b6d0cc6	144	//isMatched = MatchTracks(tpcTrack,trdTrack);
4b6d0cc6	145	isMatched = MatchTracks(extTPCTrack,trdTrack);
ec6160d5	146	if(!isMatched) continue;
ec6160d5	147
4b6d0cc6	148	//HLTInfo("-----tpc track after matching: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",tpcTrack->GetAlpha(),tpcTrack->GetX(),tpcTrack->GetY(),tpcTrack->GetZ(),tpcTrack->GetSnp(),tpcTrack->GetTgl(),tpcTrack->GetSigned1Pt());
	149	HLTInfo("-----ext tpc track after matching: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",extTPCTrack->GetAlpha(),extTPCTrack->GetX(),extTPCTrack->GetY(),extTPCTrack->GetZ(),extTPCTrack->GetSnp(),extTPCTrack->GetTgl(),extTPCTrack->GetSigned1Pt());
	150	HLTInfo("-----trd track after matching: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",trdTrack->GetAlpha(),trdTrack->GetX(),trdTrack->GetY(),trdTrack->GetZ(),trdTrack->GetSnp(),trdTrack->GetTgl(),trdTrack->GetSigned1Pt());
ec6160d5	151
4b6d0cc6	152	// propagate TRD tracks to the merging radius
	153	//isOk = AliTracker::PropagateTrackTo(tpcTrack,kMergeRadius,tpcTrack->GetMass(),kMaxStep,kFALSE);
	154	//if(!isOk) continue;
ec6160d5	155	isOk = AliTracker::PropagateTrackTo(trdTrack,kMergeRadius,trdTrack->GetMass(),kMaxStep,kFALSE);
	156	if(!isOk) continue;
	157
4b6d0cc6	158	HLTInfo("-----tpc before merging: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",tpcTrack->GetAlpha(),tpcTrack->GetX(),tpcTrack->GetY(),tpcTrack->GetZ(),tpcTrack->GetSnp(),tpcTrack->GetTgl(),tpcTrack->GetSigned1Pt());
4b6d0cc6	159	HLTInfo("-----trd before merging: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",trdTrack->GetAlpha(),trdTrack->GetX(),trdTrack->GetY(),trdTrack->GetZ(),trdTrack->GetSnp(),trdTrack->GetTgl(),trdTrack->GetSigned1Pt());
ec6160d5	160
	161	// merge TPC and TRD tracks
	162	// create AliESDtrack and add it to AliESDevent
	163	Bool_t isMerged = MergeTracks(tpcTrack,trdTrack,esdEvent);
a9f060e5	164	if(!isMerged)
a9f060e5	165	HLTInfo("No merged tracks");
ec6160d5	166	}
a9f060e5	167	HLTInfo("nTracksTRD %d",nTracksTRD);
4b6d0cc6	168
	169	// delete external TPC track
	170	if(extTPCTrack) delete extTPCTrack; extTPCTrack=0;
ec6160d5	171	}
a9f060e5	172	HLTInfo("nTracksTPC %d",nTracksTPC);
ec6160d5	173
	174	return kTRUE;
	175	}
	176
	177	//_____________________________________________________________________________
	178	Bool_t AliHLTGlobalTrackMerger::MergeTracks(AliESDtrack tpcTrack, AliESDtrack trdTrack, AliESDEvent *esdEvent)
	179	{
	180	// merge TPC and TRD track parameters
	181	// create new AliESDtrack with TPC+TRD merged track parameters
	182	// add AliESDtrack to AliESDEvent
	183
	184	if(!tpcTrack) return kFALSE;
	185	if(!trdTrack) return kFALSE;
	186
ec6160d5	187	Double_t trackParam[5], trackCovar[15];
	188	Double_t trackChi2;
	189	Int_t trackNDF;
	190
	191	// calculate merged track parameters
	192	Bool_t isNotOK = SmoothTracks(tpcTrack->GetParameter(), tpcTrack->GetCovariance(), tpcTrack->GetTPCchi2(), 5,
	193	trdTrack->GetParameter(), trdTrack->GetCovariance(), trdTrack->GetTRDchi2(), 5,
	194	trackParam, trackCovar, trackChi2, trackNDF,5);
	195
	196	if(isNotOK)
	197	return kFALSE;
	198
	199	//
	200	// create AliESDtrack
	201	// merged TPC+TRD information
	202	//
	203
	204	AliESDtrack track;
	205	//track.UpdateTrackParams(tpcTrack, AliESDtrack::kTPCrefit);
	206	track.SetStatus(AliESDtrack::kGlobalMerge);
	207	track.SetLabel(tpcTrack->GetLabel());
	208	track.Set(tpcTrack->GetX(),tpcTrack->GetAlpha(),trackParam,trackCovar);
	209	track.SetGlobalChi2(trackChi2);
	210
	211	//track.SetTPCLabel(tpcTrack->GetLabel());
	212	Double32_t tpcPID[AliPID::kSPECIES];
	213	tpcTrack->GetTPCpid(tpcPID);
	214	track.SetTPCpid(tpcPID);
	215	//fTPCncls=t->GetNumberOfClusters(); // no cluster on HLT
	216	//fTPCchi2=t->GetChi2();
	217
	218	//track.SetTRDLabel(trdTrack->GetLabel());
	219	Double32_t trdPID[AliPID::kSPECIES];
	220	trdTrack->GetTRDpid(trdPID);
	221	track.SetTRDpid(trdPID);
	222	//fTRDchi2 = t->GetChi2();
	223	//fTRDncls = t->GetNumberOfClusters();
	224	//for (Int_t i=0;i<6;i++) index[i]=t->GetTrackletIndex(i);
	225
	226	// add track to AliESDEvent
	227	esdEvent->AddTrack(&track);
	228
	229	return kTRUE;
	230	}
	231
	232	//_____________________________________________________________________________
	233	void AliHLTGlobalTrackMerger::SetParameter(Double_t maxy, Double_t maxz, Double_t maxsnp, Double_t maxtgl, Double_t signed1Pt)
	234	{
	235	//set parameters for merger
	236	fMaxY = maxy;
	237	fMaxZ = maxz;
	238	fMaxSnp = maxsnp;
	239	fMaxTgl = maxtgl;
	240	fMaxSigned1Pt = signed1Pt;
	241	}
	242
	243	//_____________________________________________________________________________
a0395f81	244	Bool_t AliHLTGlobalTrackMerger::MatchTracks(AliExternalTrackParam trackTPC, const AliESDtrack trackTRD)
ec6160d5	245	{
	246	// match TPC and TRD tracks
	247	// return kTRUE in case of matching
	248
	249	if(!trackTPC) return kFALSE;
	250	if(!trackTRD) return kFALSE;
	251
a0395f81	252	Double_t xTpc=trackTPC->GetX();
	253	Double_t yTpc=trackTPC->GetY();
	254	Double_t zTpc=trackTPC->GetZ();
	255	Double_t snpTpc=trackTPC->GetSnp();
	256	Double_t tglTpc=trackTPC->GetTgl();
	257	Double_t signed1PtTpc=trackTPC->GetSigned1Pt();
a9f060e5	258
a0395f81	259	Double_t xTrd=trackTRD->GetX();
	260	Double_t yTrd=trackTRD->GetY();
	261	Double_t zTrd=trackTRD->GetZ();
	262	Double_t snpTrd=trackTRD->GetSnp();
	263	Double_t tglTrd=trackTRD->GetTgl();
	264	Double_t signed1PtTrd=trackTRD->GetSigned1Pt();
a9f060e5	265
a9f060e5	266	// debug stream
18ada816	267	// if (AliTPCReconstructor::StreamLevel()>0) {
	268	// //TTreeSRedirector &cstream = *fDebugStreamer;
	269	// *fDebugStreamer<<"match"<<
a0395f81	270	// "xTpc="<<xTpc<<
	271	// "yTpc="<<yTpc<<
	272	// "zTpc="<<zTpc<<
	273	// "snpTpc="<<snpTpc<<
	274	// "tglTpc="<<tglTpc<<
	275	// "signed1PtTpc="<<signed1PtTpc<<
	276	// "xTrd="<<xTrd<<
	277	// "yTrd="<<yTrd<<
	278	// "zTrd="<<zTrd<<
	279	// "snpTrd="<<snpTrd<<
	280	// "tglTrd="<<tglTrd<<
	281	// "signed1PtTrd="<<signed1PtTrd<<
18ada816	282	// "\n";
18ada816	283	// }
a9f060e5	284
a0395f81	285	if (TMath::Abs(xTpc-xTrd) > 0) {/* get rid of warning*/;}
	286	if (TMath::Abs(yTpc-yTrd) > fMaxY) return kFALSE;
	287	if (TMath::Abs(zTpc-zTrd) > fMaxZ) return kFALSE;
	288	if (TMath::Abs(snpTpc-snpTrd) > fMaxSnp) return kFALSE;
	289	if (TMath::Abs(tglTpc-tglTrd) > fMaxTgl) return kFALSE;
	290	if (TMath::Abs(signed1PtTpc-signed1PtTrd) > fMaxSigned1Pt) return kFALSE;
ec6160d5	291
	292	return kTRUE;
	293	}
	294
	295	//_____________________________________________________________________________
	296	Bool_t AliHLTGlobalTrackMerger::SmoothTracks( const Double_t T1[], const Double_t C1[], Double_t Chi21, Int_t NDF1,
	297	const Double_t T2[], const Double_t C2[], Double_t Chi22, Int_t NDF2,
	298	Double_t T [], Double_t C [], Double_t &Chi2, Int_t &NDF,
	299	Int_t N )
	300	{
	301	//* Smooth two tracks with parameter vectors of size N
	302	//*
	303	//* Input:
	304	//*
	305	//* T1[N], T2[N] - tracks
	306	//* C1[N(N+1)/2], C2[N(N+1)/2] - covariance matrices in low-diagonal form:
	307	//* C = { c00,
	308	//* c10, c11,
	309	//* c20, c21, c22,
	310	//* ... };
	311	//* Chi2{1,2}, NDF{1,2} - \Chi^2 and "Number of Degrees of Freedom" values for both tracks
	312	//* Output:
	313	//*
	314	//* T[N], C[N] ( can be aqual to {T1,C1}, or {T2,C2} )
	315	//* Chi2, NDF
	316	//*
	317	//* returns error flag (0 means OK, 1 not OK )
	318
	319	Int_t M = N*(N+1)/2;
	320
	321	Double_t A[M];
	322	Double_t K[N*N];
	323
	324	for(Int_t k=0; k<M; k++) A[k] = C1[k] + C2[k];
	325	Bool_t err = InvertS(A,N);
	326	if( err ) return 1;
	327
	328	Chi2 = Chi21 + Chi22;
	329	NDF = NDF1 + NDF2;
	330
	331	MultSSQ( C1, A, K, N);
	332	Double_t r[N];
	333	for( Int_t k=0; k<N;k++) r[k] = T1[k] - T2[k];
	334	for( Int_t k=0; k<N;k++ )
	335	for( Int_t l=0;l<N;l++) T[k] = T1[k] - K[kN+l]r[l];
	336
	337	for( Int_t ind=0,i=0; i<N; i++ ){
	338	for( Int_t j=0; j<i; j++ ) Chi2+= 2r[i]r[j]*A[ind++];
	339	Chi2+= r[i]r[i]A[ind++];
	340	}
	341	NDF+=N;
	342
	343	for( Int_t l=0; l<N; l++ ) K[ (N+1)*l ] -= 1;
	344
	345	for( Int_t ind = 0, l=0; l<N; ++l ){
	346	for( Int_t j=0; j<=l; ++j, ind++ ){
	347	A[ind] = 0;
	348	for( Int_t k=0; k<N; ++k ) A[ind] -= K[lN+k] C1[IndexS(j,k)];
	349	}
	350	}
	351	for( Int_t l=0; l<N; l++ ) C[l] = A[l];
	352	return 0;
	353	}
	354
355	//_____________________________________________________________________________
356	void AliHLTGlobalTrackMerger::MultSSQ( const Double_t A, const Double_t B, Double_t *C, Int_t N )
357	{
a0395f81	358	// no clue
ec6160d5	359	for( Int_t ind=0, i=0; i<N; ++i ){
	360	for( Int_t j=0; j<N; ++j, ++ind ){
	361	C[ind] = 0;
	362	for( Int_t k=0; k<N; ++k ) C[ind] += A[IndexS(i,k)] * B[IndexS(k,j)];
	363	}
	364	}
	365	}
	366
	367	//_____________________________________________________________________________
	368	Bool_t AliHLTGlobalTrackMerger::InvertS( Double_t A[], Int_t N )
	369	{
	370	//* input: simmetric > 0 NxN matrix A = {a11,a21,a22,a31..a33,..}
	371	//* output: inverse A, in case of problems fill zero and return 1
	372	//*
	373	//* A->low triangular Anew : A = Anew x Anew^T
	374	//* method:
	375	//* for(j=1,N) for(i=j,N) Aij=(Aii-sum_{k=1}^{j-1}Aik*Ajk )/Ajj
	376	//*
	377
	378	Bool_t ret = 0;
	379
	380	const Double_t ZERO = 1.E-20;
	381
	382	{
	383	Double_t j1 = A, jj = A;
	384	for( Int_t j=1; j<=N; j1+=j++, jj+=j ){
	385	Double_t *ik = j1, x = 0;
	386	while( ik!=jj ){
	387	x -= (ik) (*ik);
	388	ik++;
	389	}
	390	x += *ik;
	391	if( x > ZERO ){
	392	x = sqrt(x);
	393	*ik = x;
	394	ik++;
	395	x = 1 / x;
	396	for( Int_t step=1; step<=N-j; ik+=++step ){ // ik==Ai1
	397	Double_t sum = 0;
bb7774d6	398	for( Double_t jk=j1; jk!=jj; sum += ((jk++)) * (*(ik++)) ) {}
ec6160d5	399	ik = (ik - sum) * x; // ik == Aij
	400	}
	401	}else{
	402	Double_t *ji=jj;
	403	for( Int_t i=j; i<N; i++ ) *(ji+=i) = 0.;
	404	ret = 1;
	405	}
	406	}
	407	}
	408
	409	//* A -> Ainv
	410	//* method :
	411	//* for(i=1,N){
	412	//* Aii = 1/Aii;
	413	//* for(j=1,i-1) Aij=-(sum_{k=j}^{i-1} Aik * Akj) / Aii ;
	414	//* }
	415
	416	{
	417	Double_t ii=A,ij=A;
	418	for( Int_t i = 1; i<=N; ij=ii+1, ii+=++i ){
	419	if( *ii > ZERO ){
	420	Double_t x = -(ii = 1./ ii);
	421	{
	422	Double_t *jj = A;
	423	for( Int_t j=1; j<i; jj+=++j, ij++ ){
	424	Double_t ik = ij, kj = jj, sum = 0.;
	425	for( Int_t k=j; ik!=ii; kj+=k++, ik++ ){
	426	sum += ik *kj;
	427	}
	428	kj = sum x;
	429	}
	430	}
	431	}else{
	432	for( Double_t *ik = ij; ik!=ii+1; ik++ ){
	433	*ik = 0.;
	434	}
	435	ret = 1;
	436	}
	437	}
	438	}
	439
	440	//* A -> A^T x A
	441	//* method:
	442	//* Aij = sum_{k=i}^N Aki * Akj
	443
	444	{
	445	Double_t ii=A, ij=A;
	446	for( Int_t i=1; i<=N; ii+=++i ){
	447	do{
	448	Double_t ki = ii, kj = ij, sum = 0.;
	449	for( Int_t k=i; k<=N; ki+=k, kj+=k++ ) sum += (ki) (*kj);
	450	*ij = sum;
	451	}while( (ij++)!=ii );
	452	}
	453	}
	454	return ret;
	455	}
	456
	457	//_____________________________________________________________________________
a0395f81	458	void AliHLTGlobalTrackMerger::PropagateTracksToDCA(const AliESDEvent *esdEvent)
ec6160d5	459	{
	460	// try to propagate all tracks to DCA to primary vertex
	461	if(!esdEvent) return;
	462
	463	const Double_t kBz = esdEvent->GetMagneticField();
	464	const Double_t kSmallRadius = 2.8; // [cm] something less than the beam pipe radius
	465	const Double_t kMaxStep = 10.0; // [cm] track propagation step
	466	Bool_t isOK = kFALSE;
	467
	468	Int_t nTracks = esdEvent->GetNumberOfTracks();
	469	for(Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
	470	AliESDtrack *track = esdEvent->GetTrack(iTrack);
	471	if(!track) continue;
	472
ec6160d5	473	// propagate to small radius (material budget included)
	474	isOK = AliTracker::PropagateTrackTo(track,kSmallRadius,track->GetMass(),kMaxStep,kFALSE);
	475
ec6160d5	476	// relate tracks to DCA to primary vertex
	477	if(isOK)
	478	{
	479	track->RelateToVertex(fVertex, kBz, kVeryBig);
a9f060e5	480	HLTInfo("1-------: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",track->GetAlpha(),track->GetX(),track->GetY(),track->GetZ(),track->GetSnp(),track->GetTgl(),track->GetSigned1Pt());
	481	}
	482
	483	//
	484	// the same procedure must be repeated for TPCinner (TPC only) tracks
	485	//
	486	AliExternalTrackParam tpcTrack = (AliExternalTrackParam )track->GetTPCInnerParam();
	487	if(!tpcTrack) continue;
	488
	489	// propagate tpcTracks to small radius (material budget included)
	490	isOK = AliTracker::PropagateTrackTo(tpcTrack,kSmallRadius,track->GetMass(),kMaxStep,kFALSE);
	491
	492	// relate tracks to DCA to primary vertex
	493	if(isOK)
	494	{
	495	Double_t par[2], cov[3];
	496	tpcTrack->PropagateToDCA(fVertex, kBz, kVeryBig,par,cov);
	497	HLTInfo("2-------: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",tpcTrack->GetAlpha(),tpcTrack->GetX(),tpcTrack->GetY(),tpcTrack->GetZ(),tpcTrack->GetSnp(),tpcTrack->GetTgl(),tpcTrack->GetSigned1Pt());
ec6160d5	498	}
	499	}
	500	}