HLT/TPCLib/AliHLTTPCCATracker.cxx

   1 /** \class AliHLTTPCCATracker
   2 <pre>
   3 //_____________________________________________________________
   4 // AliHLTTPCCATracker
   5 //
   6 //
   7 // Author: Ivan Kisel
   8 // Copyright &copy ALICE HLT Group
   9 </pre>
  10 */
  11
  12 #define WRITETRACKS 1
  13
  14 #include <sys/time.h>
  15
  16
  17 #include "AliHLTTPCRootTypes.h"
  18 #include "AliHLTTPCSpacePointData.h"
  19 #include "AliHLTTPCLogging.h"
  20 #include "AliHLTTPCVertex.h"
  21 #include "AliHLTTPCConfMapTrack.h"
  22 #include "AliHLTTPCConfMapPoint.h"
  23 #include "AliHLTTPCTrackArray.h"
  24 #include "AliHLTTPCTransform.h"
  25 #include "AliHLTTPCCATracker.h"
  26 #include "AliHLTTPCTrackSegmentData.h"
  27
  28 #include "TApplication.h"
  29 #include "TStyle.h"
  30 #include "TCanvas.h"
  31 #include "TFrame.h"
  32 #include "TMarker.h"
  33 #include "TPaveText.h"
  34 #include "TEllipse.h"
  35 #include "TText.h"
  36 #include "TROOT.h"
  37
  38 #include <iostream>
  39
  40 #if __GNUC__ >= 3
  41 using namespace std;
  42 #endif
  43
  44 static TList *listHisto;
  45 static TH1F *h_NClusters;
  46 static TH2F *h_ClustersXY[10];
  47
  48 static Int_t h_Event;
  49
  50 //#define DRAW
  51
  52 #ifdef DRAW
  53
  54 static TApplication *myapp;
  55 static TCanvas *YX, *ZX;
  56
  57 static bool ask = false;
  58
  59 #endif //DRAW
  60
  61 ClassImp(AliHLTTPCCATracker)
  62
  63 // ----------------------------------------------------------------------------------
  64 AliHLTTPCCATracker::AliHLTTPCCATracker()
  65 {
  66   //Default constructor
  67   fVertex = NULL;
  68   fTrack = NULL;
  69   fHit = NULL;
  70   fVolume = NULL;
  71   fRow = NULL;
  72   fBench = (Bool_t)true;
  73   fVertexConstraint = (Bool_t)true;
  74   fParamSet[0]=0;
  75   fParamSet[1]=0;
  76   //JMT 2006/11/13
  77   fOutputPtr= NULL;
  78   fOutputNTracks=0;
  79   fOutputSize=0;
  80 }
  81
  82 // ----------------------------------------------------------------------------------
  83 AliHLTTPCCATracker::~AliHLTTPCCATracker()
  84 {
  85   // Destructor.
  86   if(fVolume) {
  87     delete [] fVolume;
  88   }
  89   if(fRow) {
  90     delete [] fRow;
  91   }
  92   if(fHit) {
  93     delete [] fHit;
  94   }
  95   if(fTrack) {
  96     delete fTrack;
  97   }
  98   //JMT 2006/11/13
  99   if(fOutputPtr){
 100     fOutputPtr = NULL;
 101     delete fOutputPtr;
 102   }
 103
 104 }
 105
 106 // ----------------------------------------------------------------------------------
 107 void AliHLTTPCCATracker::CACreateHistos()
 108 {
 109
 110     // IK histogramming
 111     static bool first_call = true;
 112     if (first_call){
 113       first_call = false;
 114
 115       h_Event = 0;
 116
 117       TDirectory *curdir = gDirectory;
 118       TDirectory *histodir = gROOT->mkdir("histodir");
 119       histodir->cd();
 120
 121       h_NClusters  = new TH1F("h_NClusters", "Number of clusters in event", 100, 0.0, 1000.0);
 122
 123       h_ClustersXY[0] = new TH2F("h_ClustersXY_0", "Clusters in XY plane Ev 0", 100, 50.0, 300.0, 100, -100.0, 100.0);
 124       h_ClustersXY[1] = new TH2F("h_ClustersXY_1", "Clusters in XY plane Ev 1", 100, 50.0, 300.0, 100, -100.0, 100.0);
 125       h_ClustersXY[2] = new TH2F("h_ClustersXY_2", "Clusters in XY plane Ev 2", 100, 50.0, 300.0, 100, -100.0, 100.0);
 126       h_ClustersXY[3] = new TH2F("h_ClustersXY_3", "Clusters in XY plane Ev 3", 100, 50.0, 300.0, 100, -100.0, 100.0);
 127       h_ClustersXY[4] = new TH2F("h_ClustersXY_4", "Clusters in XY plane Ev 4", 100, 50.0, 300.0, 100, -100.0, 100.0);
 128       h_ClustersXY[5] = new TH2F("h_ClustersXY_5", "Clusters in XY plane Ev 5", 100, 50.0, 300.0, 100, -100.0, 100.0);
 129       h_ClustersXY[6] = new TH2F("h_ClustersXY_6", "Clusters in XY plane Ev 6", 100, 50.0, 300.0, 100, -100.0, 100.0);
 130       h_ClustersXY[7] = new TH2F("h_ClustersXY_7", "Clusters in XY plane Ev 7", 100, 50.0, 300.0, 100, -100.0, 100.0);
 131       h_ClustersXY[8] = new TH2F("h_ClustersXY_8", "Clusters in XY plane Ev 8", 100, 50.0, 300.0, 100, -100.0, 100.0);
 132       h_ClustersXY[9] = new TH2F("h_ClustersXY_9", "Clusters in XY plane Ev 9", 100, 50.0, 300.0, 100, -100.0, 100.0);
 133
 134       listHisto = gDirectory->GetList();
 135       curdir->cd();
 136
 137 #ifdef DRAW
 138
 139       myapp = new TApplication("myapp", 0, 0);
 140
 141       gStyle->SetCanvasBorderMode(0);
 142       gStyle->SetCanvasBorderSize(1);
 143       gStyle->SetCanvasColor(0);
 144
 145
 146       YX = new TCanvas ("YX", "YX (Pad-Row) window", -10, 0, 680, 400);
 147       YX->Range(-50.0, 80.0, 50.0, 250.0);
 148       YX->Clear();
 149       YX->Draw();
 150
 151       ZX = new TCanvas ("ZX", "ZX window", -700, 0, 680, 400);
 152       ZX->Range(-250.0, 80.0, 250.0, 250.0);
 153       ZX->Clear();
 154       ZX->Draw();
 155
 156 #ifdef XXX
 157       char symbol;
 158       if (ask){
 159         do{
 160           std::cin.get(symbol);
 161           if (symbol == 'r')
 162             ask = false;
 163         } while (symbol != '\n');
 164       }
 165 #endif //XXX
 166
 167 #endif //DRAW
 168
 169     }
 170 }
 171
 172 // ----------------------------------------------------------------------------------
 173 void AliHLTTPCCATracker::CAWriteHistos()
 174 {
 175     // IK Open the output file and write histogramms
 176
 177     TFile *outfile = new TFile("ConfMapper_histo.root", "RECREATE");
 178     TIter hiter(listHisto);
 179     while (TObject *obj = hiter()) obj->Write();
 180     outfile->Close();
 181 }
 182
 183 // ----------------------------------------------------------------------------------
 184 Bool_t AliHLTTPCCATracker::ReadHits(UInt_t count, AliHLTTPCSpacePointData* hits )
 185 {
 186
 187   //read hits
 188   Int_t nhit=(Int_t)count;
 189
 190   vector<CAHit> vec_hits;
 191   vec_hits.clear();
 192
 193   for (Int_t i=0;i<nhit;i++)
 194     {
 195       fHit[i+fClustersUnused].Reset();
 196       fHit[i+fClustersUnused].ReadHits(&(hits[i]));
 197     }
 198
 199   fClustersUnused += nhit;
 200
 201   return true;
 202 }
 203
 204 // ----------------------------------------------------------------------------------
 205 void AliHLTTPCCATracker::CAInitialize()
 206 {
 207   vec_hits.clear();
 208   vec_patch_tracks.clear();
 209   vec_slice_tracks.clear();
 210
 211   return;
 212 }
 213
 214 // ----------------------------------------------------------------------------------
 215 void AliHLTTPCCATracker::CAReadPatchHits(Int_t patch, UInt_t count, AliHLTTPCSpacePointData* hits )
 216 {
 217   //read hits
 218   Int_t nhit=(Int_t)count;
 219
 220   patch_first_hit_ind = vec_hits.size();
 221   patch_last_hit_ind  = patch_first_hit_ind + nhit - 1;
 222
 223   for (Int_t i=0;i<nhit;i++)
 224     {
 225       AliHLTTPCSpacePointData* pSP = &(hits[i]);
 226
 227       CAHit hit;
 228       hit.x = pSP->fX;
 229       hit.y = pSP->fY;
 230       hit.z = pSP->fZ;
 231
 232       hit.errx = 1.0;//pSP->fSigmaY2;
 233       if (patch <= 3)
 234         hit.erry = 10.0;//1.0;//pSP->fSigmaY2;
 235       else
 236         hit.erry = 10.0;//pSP->fSigmaY2;
 237       hit.errz = 1.0;//pSP->fSigmaZ2;
 238
 239       hit.index = i;
 240       hit.counter = 0;
 241
 242       vec_hits.push_back(hit);
 243     }
 244
 245 #ifdef DRAW
 246
 247   TMarker *mhit = new TMarker(0.0, 0.0, 6);
 248   mhit->SetMarkerColor(kBlue);
 249
 250   for (Int_t i = patch_first_hit_ind; i <= patch_last_hit_ind; i++)
 251     {
 252       YX->cd();
 253       mhit->DrawMarker(vec_hits[i].y, vec_hits[i].x);
 254
 255       ZX->cd();
 256       mhit->DrawMarker(vec_hits[i].z, vec_hits[i].x);
 257     }
 258
 259   delete mhit;
 260
 261   YX->cd(); YX->Update(); YX->Draw();
 262   ZX->cd(); ZX->Update(); ZX->Draw();
 263
 264 #endif //DRAW
 265
 266   return;
 267 }
 268
 269 // ----------------------------------------------------------------------------------
 270 void AliHLTTPCCATracker::CAFindPatchTracks(Int_t patch)
 271 {
 272
 273 #ifdef DRAW
 274
 275   TLine *line = new TLine();
 276   line->SetLineColor(kBlack);
 277
 278 #endif //DRAW
 279
 280   Double_t dist01, dist12, dist02, mindist, chi2, minchi2;
 281   CAHit *ph0, *ph1, *ph2, *phbest, *next;
 282
 283   Double_t x, y, z, ty, tz, cov_y, cov_ty, cov_yty, cov_z, cov_tz, cov_ztz, ye, ze, dx, dy, dz, sy2, sz2;
 284
 285   for (Int_t i1 = patch_first_hit_ind; i1 < patch_last_hit_ind; i1++) // <nhit-1
 286     {
 287       ph1  = &(vec_hits[i1]);
 288       next = &(vec_hits[i1+1]);
 289
 290       mindist = 20.0;
 291       minchi2 = 1.0;
 292       phbest = NULL;
 293
 294       if (ph1->counter != 0){ // the previous hit exists
 295         x = ph1->x;
 296         y = ph1->y;
 297         z = ph1->z;
 298
 299         dx = (ph1->x-ph0->x); // assume different x !!!
 300         if (abs(dx) < 0.0001) dx = 0.0001;
 301         ty = (ph1->y-ph0->y)/dx;
 302         tz = (ph1->z-ph0->z)/dx;
 303
 304         cov_y   = ph1->erry*ph1->erry;
 305         cov_ty  = (ph0->erry*ph0->erry + ph1->erry*ph1->erry)/(dx*dx);
 306         cov_yty = (ph1->erry*ph1->erry)/dx;
 307
 308         cov_z   = ph1->errz*ph1->errz;
 309         cov_tz  = (ph0->errz*ph0->errz + ph1->errz*ph1->errz)/(dx*dx);
 310         cov_ztz = (ph1->errz*ph1->errz)/dx;
 311
 312       }
 313
 314       for (Int_t i2 = i1+1; i2 <= patch_last_hit_ind; i2++)
 315         {
 316           ph2 = &(vec_hits[i2]);
 317           if (ph2->x < ph1->x) continue; // no backward direction!
 318
 319           dist12 = sqrt((ph1->x-ph2->x)*(ph1->x-ph2->x) + (ph1->y-ph2->y)*(ph1->y-ph2->y) + (ph1->z-ph2->z)*(ph1->z-ph2->z));
 320           if (dist12 < mindist){
 321             if (ph1->counter == 0){
 322               mindist = dist12;
 323               phbest  = ph2;
 324             }
 325             else{
 326               // calculate chi2 distance between the hit and the line
 327               dx = (ph2->x-x); // assume different x !!!
 328               if (abs(dx) < 0.0001) dx = 0.0001;
 329
 330               ye = y + ty*dx;
 331               ze = z + tz*dx;
 332
 333               sy2 = cov_y + 2.0*cov_yty*dx + cov_ty*dx*dx;
 334               sz2 = cov_z + 2.0*cov_ztz*dx + cov_tz*dx*dx;
 335
 336               dy = ph2->y - ye;
 337               dz = ph2->z - ze;
 338
 339               chi2 = (dy*dy)/(sy2+ph2->erry*ph2->erry) + (dz*dz)/(sz2+ph2->errz*ph2->errz);
 340
 341               // check the straight line model
 342               if (chi2 < minchi2){
 343                 mindist = dist12;
 344                 phbest  = ph2;
 345               }
 346
 347             }
 348           }
 349         }
 350
 351       if (phbest != NULL){// closest hit found
 352
 353         // exchange jbest with the next hit
 354         Double_t x_tmp = next->x;
 355         Double_t y_tmp = next->y;
 356         Double_t z_tmp = next->z;
 357
 358         Double_t errx_tmp = next->errx;
 359         Double_t erry_tmp = next->erry;
 360         Double_t errz_tmp = next->errz;
 361
 362         Int_t index_tmp = next->index;
 363         Int_t counter_tmp = next->counter;
 364
 365         next->x = phbest->x;
 366         next->y = phbest->y;
 367         next->z = phbest->z;
 368
 369         next->errx = phbest->errx;
 370         next->erry = phbest->erry;
 371         next->errz = phbest->errz;
 372
 373         next->index = phbest->index;
 374         next->counter = phbest->counter;
 375
 376         phbest->x = x_tmp;
 377         phbest->y = y_tmp;
 378         phbest->z = z_tmp;
 379
 380         phbest->errx = errx_tmp;
 381         phbest->erry = erry_tmp;
 382         phbest->errz = errz_tmp;
 383
 384         phbest->index = index_tmp;
 385         phbest->counter = counter_tmp;
 386
 387         if (ph1->counter == 0)
 388           ph1->counter = 1;
 389         next->counter = ph1->counter+1;
 390
 391         ph0 = ph1;
 392         dist01 = mindist;
 393
 394 #ifdef DRAW
 395
 396         YX->cd();
 397         line->DrawLine(ph1->y, ph1->x, next->y, next->x);
 398
 399         ZX->cd();
 400         line->DrawLine(ph1->z, ph1->x, next->z, next->x);
 401
 402 #endif //DRAW
 403
 404       }
 405
 406     }
 407
 408 #ifdef DRAW
 409
 410   delete line;
 411
 412   YX->cd(); YX->Update(); YX->Draw();
 413   ZX->cd(); ZX->Update(); ZX->Draw();
 414
 415 #endif //DRAW
 416
 417   CATrack tr;
 418   tr.nhits = 0;
 419   tr.chi2  = 0.0;
 420   tr.ndf   = -4;
 421   tr.vec_ihits.clear();
 422
 423   patch_first_track_ind = vec_patch_tracks.size();
 424   patch_last_track_ind  = patch_first_track_ind - 1;
 425
 426   CAHit *my_phf;
 427
 428   // collect tracks
 429   for (Int_t ihr = patch_last_hit_ind; ihr >= patch_first_hit_ind; ihr--){
 430
 431     CAHit* phit = &(vec_hits[ihr]);
 432
 433     if ((phit->counter < 5)&&(tr.nhits == 0))
 434       continue; // do not collect short tracks
 435
 436     if ((phit->counter != 0)&&(tr.nhits == 0))// store the first hit
 437         my_phf = phit;
 438
 439     if (phit->counter != 0){//add hit to the track
 440       tr.nhits++;
 441       tr.vec_ihits.push_back(ihr);
 442
 443       if (tr.nhits == 2){ // calculate initial track parameters
 444         CAHit* phit0 = &(vec_hits[tr.vec_ihits[0]]);
 445         CAHit* phit1 = &(vec_hits[tr.vec_ihits[1]]);
 446
 447         tr.x = phit1->x;
 448         tr.y = phit1->y;
 449         tr.z = phit1->z;
 450
 451         Double_t dx = (phit1->x-phit0->x); // assume different x !!!
 452         if (abs(dx) < 0.0001) dx = 0.0001;
 453         tr.ty = (phit1->y-phit0->y)/dx;
 454         tr.tz = (phit1->z-phit0->z)/dx;
 455
 456         tr.cov_y   = phit1->erry*phit1->erry;
 457         tr.cov_ty  = (phit0->erry*phit0->erry + phit1->erry*phit1->erry)/(dx*dx);
 458         tr.cov_yty = (phit1->erry*phit1->erry)/dx;
 459
 460         tr.cov_z   = phit1->errz*phit1->errz;
 461         tr.cov_tz  = (phit0->errz*phit0->errz + phit1->errz*phit1->errz)/(dx*dx);
 462         tr.cov_ztz = (phit1->errz*phit1->errz)/dx;
 463       }
 464
 465       if (tr.nhits > 2){
 466         // propagate the track
 467         Double_t dx = (phit->x-tr.x); // assume different x !!!
 468         if (abs(dx) < 0.0001) dx = 0.0001;
 469
 470         Double_t ye = tr.y + tr.ty*dx;
 471         Double_t ze = tr.z + tr.tz*dx;
 472
 473         Double_t cov_y   = tr.cov_y + 2.0*tr.cov_yty*dx + tr.cov_ty*dx*dx;
 474         Double_t cov_yty = tr.cov_yty + tr.cov_ty*dx;
 475         Double_t cov_ty  = tr.cov_ty;
 476
 477         Double_t cov_z   = tr.cov_z + 2.0*tr.cov_ztz*dx + tr.cov_tz*dx*dx;
 478         Double_t cov_ztz = tr.cov_ztz + tr.cov_tz*dx;
 479         Double_t cov_tz  = tr.cov_tz;
 480
 481
 482         Double_t dy = phit->y - ye;
 483         Double_t dz = phit->z - ze;
 484
 485         // add the measurement
 486
 487         Double_t w = 1.0/(phit->erry*phit->erry + cov_y);
 488
 489         tr.y  = ye + cov_y*dy*w;
 490         tr.ty = tr.ty + cov_yty*dy*w;
 491
 492         tr.cov_y   = cov_y - cov_y*cov_y*w;
 493         tr.cov_yty = cov_yty - cov_y*cov_yty*w;
 494         tr.cov_ty  = cov_ty -cov_yty*cov_yty*w;
 495
 496         tr.chi2 += dy*dy*w;
 497         tr.ndf  += 1;
 498
 499         w = 1.0/(phit->errz*phit->errz + cov_z);
 500
 501         tr.z  = ze + cov_z*dz*w;
 502         tr.tz = tr.tz + cov_ztz*dz*w;
 503
 504         tr.cov_z   = cov_z - cov_z*cov_z*w;
 505         tr.cov_ztz = cov_ztz - cov_z*cov_ztz*w;
 506         tr.cov_tz  = cov_tz -cov_ztz*cov_ztz*w;
 507
 508         tr.x += dx;
 509
 510         tr.chi2 += dz*dz*w;
 511         tr.ndf  += 1;
 512       }
 513     }
 514     if (phit->counter == 1){// store the track
 515
 516       tr.good  =  1;
 517       tr.used  =  0;
 518       tr.next  = -1;
 519       tr.patch = patch;
 520
 521       Double_t trdist = sqrt((my_phf->x-phit->x)*(my_phf->x-phit->x)+
 522                              (my_phf->y-phit->y)*(my_phf->y-phit->y)+
 523                              (my_phf->z-phit->z)*(my_phf->z-phit->z));
 524
 525       if ((trdist > 1.0)&&(tr.chi2/tr.ndf < 10.0)){
 526         vec_patch_tracks.push_back(tr);
 527         patch_last_track_ind++;
 528       }
 529
 530       // reset the track
 531       tr.nhits = 0;
 532       tr.chi2  = 0.0;
 533       tr.ndf   = -4;
 534       tr.vec_ihits.clear();
 535     }
 536   }
 537
 538   // sort tracks according (currently) to number of hits
 539   //sort(vec_patch_tracks.begin(), vec_patch_tracks.end(), compareNhits); //don't need to sort here
 540
 541
 542 #ifdef DRAW
 543
 544   TLine *trline = new TLine();
 545   trline->SetLineColor(kGreen);
 546
 547   for (Int_t itr = patch_first_track_ind; itr <= patch_last_track_ind; itr++){
 548
 549     CATrack *ptr = &(vec_patch_tracks[itr]);
 550
 551     CAHit* phitf = &(vec_hits[ptr->vec_ihits.back()]);
 552     CAHit* phitl = &(vec_hits[ptr->vec_ihits.front()]);
 553
 554     YX->cd();
 555     Double_t yf = ptr->y + ptr->ty*(phitf->x-ptr->x);
 556     Double_t yl = ptr->y + ptr->ty*(phitl->x-ptr->x);
 557     trline->DrawLine(yf, phitf->x, yl, phitl->x);
 558
 559     ZX->cd();
 560     Double_t zf = ptr->z + ptr->tz*(phitf->x-ptr->x);
 561     Double_t zl = ptr->z + ptr->tz*(phitl->x-ptr->x);
 562     trline->DrawLine(zf, phitf->x, zl, phitl->x);
 563
 564   }
 565
 566   delete trline;
 567
 568   YX->cd(); YX->Update(); YX->Draw();
 569   ZX->cd(); ZX->Update(); ZX->Draw();
 570
 571 #endif //DRAW
 572
 573   return;
 574 }
 575
 576 // ----------------------------------------------------------------------------------
 577 void AliHLTTPCCATracker::CAFindSliceTracks()
 578 {
 579
 580   Int_t ntracks = vec_patch_tracks.size();
 581
 582   // merge patch tracks into slice tracks
 583   // start with the longest tracks in the first patches
 584
 585   // sort tracks according (currently) to patch and within patch - number of hits
 586   sort(vec_patch_tracks.begin(), vec_patch_tracks.end(), compareCATracks);
 587
 588   for (Int_t itr1 = 0; itr1 < ntracks-1; itr1++){
 589
 590     CATrack *ptr1 = &(vec_patch_tracks[itr1]);
 591
 592     Double_t maxdisty = 10.0;
 593     Double_t maxdistz =  5.0;
 594
 595     Double_t mindist = 10.0;
 596
 597     Int_t next = -1;
 598 #ifdef XXX
 599     for (Int_t itr2 = itr1+1; itr2 < ntracks; itr2++){
 600
 601       CATrack *ptr2 = &(vec_patch_tracks[itr2]);
 602
 603       if (ptr2->patch == ptr1->patch) continue; // the same patch - no prolongation
 604       if (ptr2->patch >  ptr1->patch+1) break;  // sorted tracks  - no prolongation skippeng over patches
 605
 606       if (ptr2->used == 1) continue; // already matched
 607
 608       CAHit* phitf2 = &(vec_hits[ptr2->vec_ihits.back()]);
 609       CAHit* phitl2 = &(vec_hits[ptr2->vec_ihits.front()]);
 610
 611       // extrapolate tr1 to the both ends of tr2
 612
 613       Double_t dyf = ptr1->y + ptr1->ty*(phitf2->x-ptr1->x) - phitf2->y;
 614       Double_t dyl = ptr1->y + ptr1->ty*(phitl2->x-ptr1->x) - phitl2->y;
 615
 616       Double_t dzf = ptr1->z + ptr1->tz*(phitf2->x-ptr1->x) - phitf2->z;
 617       Double_t dzl = ptr1->z + ptr1->tz*(phitl2->x-ptr1->x) - phitl2->z;
 618
 619       // roughly similar tracks ?
 620       if ((abs(dyf) > maxdisty)||(abs(dyl) > maxdisty)||(abs(dzf) > maxdistz)||(abs(dzl) > maxdistz)) continue;
 621
 622       // roughly parallel tracks ?
 623       Double_t dist = sqrt((dyf-dyl)*(dyf-dyl)+(dzf-dzl)*(dzf-dzl));
 624       if (dist > mindist) continue;
 625
 626       mindist = dist;
 627       next    = itr2;
 628
 629     }
 630 #endif
 631
 632     // found track prolongations?
 633     if (next != -1){
 634       CATrack *ptr2 = &(vec_patch_tracks[next]);
 635
 636       ptr1->next = next;
 637       ptr2->used = 1;
 638
 639
 640     }
 641
 642   }
 643
 644   CATrack tr;
 645   tr.nhits = 0;
 646   tr.chi2  = 0.0;
 647   tr.ndf   = -4;
 648   tr.vec_ihits.clear();
 649
 650
 651   //collect tracks
 652   for (Int_t itr = 0; itr < ntracks; itr++){
 653
 654     CATrack *ptr = &(vec_patch_tracks[itr]);
 655     if (ptr->used) continue; // start with a track not used in prolongations
 656
 657     Int_t first = 1;
 658     do{
 659       if (first == 1){
 660         first = 0;
 661       }
 662       else{
 663         ptr = &(vec_patch_tracks[ptr->next]);
 664       }
 665
 666       for (Int_t ih = 0; ih < ptr->nhits; ih++){
 667         tr.vec_ihits.push_back(ptr->vec_ihits[ih]);
 668         tr.nhits++;
 669       }
 670     }while(ptr->next != -1);
 671
 672     //sort hits according to increasing x
 673     sort(tr.vec_ihits.begin(), tr.vec_ihits.end(), compareCAHitsX);
 674
 675     vec_slice_tracks.push_back(tr);
 676
 677     tr.nhits = 0;
 678     tr.chi2  = 0.0;
 679     tr.ndf   = -4;
 680     tr.vec_ihits.clear();
 681
 682   }
 683
 684 #ifdef DRAW
 685
 686   TLine *trline = new TLine();
 687   trline->SetLineColor(kRed);
 688
 689 #endif //DRAW
 690
 691 #if WRITETRACKS
 692   UInt_t size = 0;
 693   UInt_t nTracks = 0;
 694 #endif
 695
 696
 697   for (vector<CATrack>::iterator trIt = vec_slice_tracks.begin(); trIt != vec_slice_tracks.end(); trIt++){
 698
 699     CAHit* phit0 = &(vec_hits[trIt->vec_ihits[trIt->nhits-1]]);
 700     CAHit* phit1 = &(vec_hits[trIt->vec_ihits[trIt->nhits-2]]);
 701
 702     trIt->x = phit1->x;
 703     trIt->y = phit1->y;
 704     trIt->z = phit1->z;
 705
 706     Double_t dx = (phit1->x-phit0->x); // assume different x !!!
 707     if (abs(dx) < 0.0001) dx = 0.0001;
 708     trIt->ty = (phit1->y-phit0->y)/dx;
 709     trIt->tz = (phit1->z-phit0->z)/dx;
 710
 711     trIt->cov_y   = phit1->erry*phit1->erry;
 712     trIt->cov_ty  = (phit0->erry*phit0->erry + phit1->erry*phit1->erry)/(dx*dx);
 713     trIt->cov_yty = (phit1->erry*phit1->erry)/dx;
 714
 715     trIt->cov_z   = phit1->errz*phit1->errz;
 716     trIt->cov_tz  = (phit0->errz*phit0->errz + phit1->errz*phit1->errz)/(dx*dx);
 717     trIt->cov_ztz = (phit1->errz*phit1->errz)/dx;
 718
 719     for (Int_t ih = trIt->nhits-3; ih >= 0; ih--){
 720
 721       CAHit* phit = &(vec_hits[trIt->vec_ihits[ih]]);
 722
 723       // propagate the track
 724       Double_t dx = (phit->x-trIt->x); // assume different x !!!
 725       if (abs(dx) < 0.0001) dx = 0.0001;
 726
 727       Double_t ye = trIt->y + trIt->ty*dx;
 728       Double_t ze = trIt->z + trIt->tz*dx;
 729
 730       Double_t cov_y   = trIt->cov_y + 2.0*trIt->cov_yty*dx + trIt->cov_ty*dx*dx;
 731       Double_t cov_yty = trIt->cov_yty + trIt->cov_ty*dx;
 732       Double_t cov_ty  = trIt->cov_ty;
 733
 734       Double_t cov_z   = trIt->cov_z + 2.0*trIt->cov_ztz*dx + trIt->cov_tz*dx*dx;
 735       Double_t cov_ztz = trIt->cov_ztz + trIt->cov_tz*dx;
 736       Double_t cov_tz  = trIt->cov_tz;
 737
 738       Double_t dy = phit->y - ye;
 739       Double_t dz = phit->z - ze;
 740
 741       // add the measurement
 742
 743       Double_t w = 1.0/(phit->erry*phit->erry + cov_y);
 744
 745       trIt->y  = ye + cov_y*dy*w;
 746       trIt->ty = trIt->ty + cov_yty*dy*w;
 747
 748       trIt->cov_y   = cov_y - cov_y*cov_y*w;
 749       trIt->cov_yty = cov_yty - cov_y*cov_yty*w;
 750       trIt->cov_ty  = cov_ty -cov_yty*cov_yty*w;
 751
 752       trIt->chi2 += dy*dy*w;
 753       trIt->ndf  += 1;
 754
 755       w = 1.0/(phit->errz*phit->errz + cov_z);
 756
 757       trIt->z  = ze + cov_z*dz*w;
 758       trIt->tz = trIt->tz + cov_ztz*dz*w;
 759
 760       trIt->cov_z   = cov_z - cov_z*cov_z*w;
 761       trIt->cov_ztz = cov_ztz - cov_z*cov_ztz*w;
 762       trIt->cov_tz  = cov_tz -cov_ztz*cov_ztz*w;
 763
 764       trIt->chi2 += dz*dz*w;
 765       trIt->ndf  += 1;
 766
 767       trIt->x    += dx;
 768     }
 769
 770     trIt->good  =  1;
 771     trIt->used  =  0;
 772     trIt->next  = -1;
 773     trIt->patch = -1;
 774
 775     //if (trIt->chi2/trIt->ndf < 1.0)
 776     //trIt->good  =  0;
 777
 778
 779     // JMT 2006/11/13 Write Tracks to container
 780 #if WRITETRACKS
 781     CAHit* firstHit = &(vec_hits[trIt->vec_ihits[0]]);
 782     CAHit* lastHit = &(vec_hits[trIt->vec_ihits[trIt->nhits-1]]);
 783
 784     Float_t xFirst = (Float_t) firstHit->x;
 785     Float_t yFirst = (Float_t) trIt->y + trIt->ty*(xFirst-trIt->x);
 786     float_t zFirst = (Float_t) trIt->z + trIt->tz*(xFirst-trIt->x);
 787
 788     Float_t xLast = (Float_t) lastHit->x;
 789     Float_t yLast = (Float_t) trIt->y + trIt->ty*(xLast-trIt->x);
 790     Float_t zLast = (Float_t) trIt->z + trIt->tz*(xLast-trIt->x);
 791
 792     fOutputPtr->fX = xFirst;
 793     fOutputPtr->fY = yFirst;
 794     fOutputPtr->fZ = zFirst;
 795     fOutputPtr->fPt = -9876.0;
 796     fOutputPtr->fPterr = 1;
 797     fOutputPtr->fLastX = xLast;
 798     fOutputPtr->fLastY = yLast;
 799     fOutputPtr->fLastZ = zLast;
 800     fOutputPtr->fPsi = atan(trIt->ty);
 801     fOutputPtr->fTgl = trIt->tz;
 802     fOutputPtr->fPsierr = 1;
 803     fOutputPtr->fTglerr = 1;
 804     fOutputPtr->fCharge = 1;
 805     fOutputPtr->fNPoints = 0;
 806
 807     Byte_t *tmpP = (Byte_t *)fOutputPtr;
 808
 809     tmpP += sizeof(AliHLTTPCTrackSegmentData); //+fOutputPtr->fNPoints*sizeof(UInt_t);
 810     size += sizeof(AliHLTTPCTrackSegmentData); //+fOutputPtr->fNPoints*sizeof(UInt_t);
 811     fOutputPtr = (AliHLTTPCTrackSegmentData*)tmpP;
 812     nTracks++;
 813 #endif
 814
 815 #ifdef DRAW
 816     CAHit* phitf = &(vec_hits[trIt->vec_ihits[0]]);
 817     CAHit* phitl = &(vec_hits[trIt->vec_ihits[trIt->nhits-1]]);
 818
 819     Double_t xf = phitf->x;
 820     Double_t yf = trIt->y + trIt->ty*(xf-trIt->x);
 821     Double_t zf = trIt->z + trIt->tz*(xf-trIt->x);
 822
 823     Double_t xl = phitl->x;
 824     Double_t yl = trIt->y + trIt->ty*(xl-trIt->x);
 825     Double_t zl = trIt->z + trIt->tz*(xl-trIt->x);
 826
 827     YX->cd();
 828     trline->DrawLine(yf, xf, yl, xl);
 829
 830     ZX->cd();
 831     trline->DrawLine(zf, xf, zl, xl);
 832
 833 #endif //DRAW
 834
 835   }
 836
 837 #if WRITETRACKS
 838   fOutputNTracks = nTracks;
 839   fOutputSize = size;
 840   cout << "NTRACKS=" << nTracks << endl;
 841   cout << "SIZEoF=" <<  sizeof(AliHLTTPCTrackSegmentData) << endl;
 842   cout << "SIZE=" << fOutputSize << endl;
 843 #endif
 844
 845 #ifdef DRAW
 846
 847   delete trline;
 848
 849   YX->cd(); YX->Update(); YX->Draw();
 850   ZX->cd(); ZX->Update(); ZX->Draw();
 851
 852 #endif //DRAW
 853
 854   return;
 855 }