HLT/global/AliHLTGlobalTrackMerger.cxx

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