HLT/TPCLib/AliHLTTPCConfMapper.cxx

   1 // @(#) $Id$
   2 // Original: AliHLTConfMapper.cxx,v 1.26 2005/06/14 10:55:21 cvetan Exp $
   3
   4 /** \class AliHLTTPCConfMapper
   5 <pre>
   6 //_____________________________________________________________
   7 // AliHLTTPCConfMapper
   8 //
   9 // Conformal mapping base class
  10 //
  11 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
  12 // Copyright &copy ALICE HLT Group
  13 </pre>
  14 */
  15
  16 #include <cassert>
  17 #include <sys/time.h>
  18
  19 #include "AliHLTTPCRootTypes.h"
  20 #include "AliHLTTPCSpacePointData.h"
  21 #include "AliHLTTPCLogging.h"
  22 #include "AliHLTTPCVertex.h"
  23 #include "AliHLTTPCConfMapTrack.h"
  24 #include "AliHLTTPCConfMapPoint.h"
  25 #include "AliHLTTPCTrackArray.h"
  26 #include "AliHLTTPCTransform.h"
  27 #include "AliHLTTPCConfMapper.h"
  28
  29 #if __GNUC__ >= 3
  30 using namespace std;
  31 #endif
  32
  33 ClassImp(AliHLTTPCConfMapper)
  34
  35 AliHLTTPCConfMapper::AliHLTTPCConfMapper()
  36   :
  37   fBench(kTRUE),
  38   fNTracks(0),
  39   fVertex(NULL),
  40   fVertexFinder(kFALSE),
  41   fHit(),
  42   fTrack(NULL),
  43   fMaxDca(0.0), // no clue whether this is reasonable, but at least better than without initialization
  44   fVolume(NULL),
  45   fRow(NULL),
  46   fNumRowSegment(0),
  47   fNumPhiSegment(0),
  48   fNumEtaSegment(0),
  49   fNumRowSegmentPlusOne(0),
  50   fNumPhiSegmentPlusOne(0),
  51   fNumEtaSegmentPlusOne(0),
  52   fNumPhiEtaSegmentPlusOne(0),
  53   fBounds(0),
  54   fPhiHitsOutOfRange(0),
  55   fEtaHitsOutOfRange(0),
  56   fPhiMin(0),
  57   fPhiMax(0),
  58   fEtaMin(0),
  59   fEtaMax(0),
  60   fRowMin(0),
  61   fRowMax(0),
  62   fVertexConstraint(kTRUE),
  63   fGoodDist(0.0),
  64   fMaxPhi(0.0),
  65   fMaxEta(0.0),
  66   fMainVertexTracks(0),
  67   fClustersUnused(0)
  68
  69 {
  70   //Default constructor
  71   fParamSet[0]=0;
  72   fParamSet[1]=0;
  73 }
  74
  75 AliHLTTPCConfMapper::~AliHLTTPCConfMapper()
  76 {
  77   // Destructor.
  78   if(fRow) {
  79     delete [] fRow;
  80   }
  81   if(fTrack) {
  82     delete fTrack;
  83   }
  84 }
  85
  86 void AliHLTTPCConfMapper::InitVolumes()
  87 {
  88   //Data organization.
  89   //Allocate volumes, set conformal coordinates and pointers.
  90
  91   //Should be done after setting the track parameters
  92
  93   fNumRowSegmentPlusOne = AliHLTTPCTransform::GetNRows();//NumRows[0]; //Maximum 32.
  94   fNumPhiSegmentPlusOne = fNumPhiSegment+1;
  95   fNumEtaSegmentPlusOne = fNumEtaSegment+1;
  96   fNumPhiEtaSegmentPlusOne = fNumPhiSegmentPlusOne*fNumEtaSegmentPlusOne;
  97   fBounds = fNumRowSegmentPlusOne * fNumPhiSegmentPlusOne * fNumEtaSegmentPlusOne;
  98
  99   Reset();
 100
 101   fTrack = new AliHLTTPCTrackArray("AliHLTTPCConfMapTrack",10);
 102 }
 103
 104 void AliHLTTPCConfMapper::Reset()
 105 {
 106   if(fVolume) delete [] fVolume;
 107   fVolume=NULL;
 108   if(fRow) delete [] fRow;
 109   fRow=NULL;
 110
 111   fClustersUnused=0;
 112   fHit.clear();
 113 }
 114
 115 void AliHLTTPCConfMapper::InitSector(Int_t sector,Int_t *rowrange,Float_t *etarange)
 116 { //sector means slice here
 117   //Initialize tracker for tracking in a given sector.
 118   //Resets track and hit arrays.
 119   //Here it is also possible to specify a subsector, by defining
 120   //rowrange[0]=innermost row;
 121   //rowrange[1]=outermostrow;
 122   //Finally you can specify etaslices to save time (assuming a good seed from TRD...)
 123
 124   //Define tracking area:
 125   if(rowrange)
 126     {
 127       fRowMin = rowrange[0];
 128       fRowMax = rowrange[1];
 129     }
 130   else //complete sector
 131     {
 132       fRowMin = 0;
 133       fRowMax = AliHLTTPCTransform::GetNRows() - 1;
 134     }
 135   if(etarange)
 136     {
 137       fEtaMin = etarange[0];
 138       fEtaMax = sector < 18 ? etarange[1] : -etarange[1];
 139     }
 140   else
 141     {
 142       fEtaMin = 0;
 143       fEtaMax = sector < 18 ? 0.9 : -0.9;
 144     }
 145
 146   //Set the angles to sector 2:
 147   fPhiMin = -10*AliHLTTPCTransform::ToRad();//fParam->GetAngle(sector) - 10/todeg;
 148   fPhiMax =  10*AliHLTTPCTransform::ToRad();//fParam->GetAngle(sector) + 10/todeg;
 149
 150   fNTracks=0;
 151   fMainVertexTracks = 0;
 152   fClustersUnused = 0;
 153   fEtaHitsOutOfRange=0;
 154   fPhiHitsOutOfRange=0;
 155
 156   fNumRowSegment = fRowMax - fRowMin; //number of rows to be considered by tracker
 157   LOG(AliHLTTPCLog::kInformational,"AliHLTTPCConfMapper::InitSector","B-field")
 158     <<"Tracker initializing with a magnetic field of "<<AliHLTTPCTransform::GetBField()<<ENDLOG;
 159
 160   fTrack->Reset();
 161 }
 162
 163 Bool_t AliHLTTPCConfMapper::ReadHits(UInt_t count, AliHLTTPCSpacePointData* hits )
 164 {
 165   //read hits
 166   if (fHit.size()<fClustersUnused+count) fHit.resize(fClustersUnused+count);
 167   assert(fHit.size()>=fClustersUnused+count);
 168   for (Int_t i=0;(UInt_t)i<count;i++)
 169     {
 170       fHit[i+fClustersUnused].Reset();
 171       fHit[i+fClustersUnused].Read(hits[i]);
 172     }
 173   fClustersUnused += count;
 174
 175   LOG(AliHLTTPCLog::kDebug,"AliHLTTPCConfMapper::ReadHits","#hits")
 176     <<AliHLTTPCLog::kDec<<"#hits: "<<count<<" total: "<<fClustersUnused<<ENDLOG;
 177
 178   return true;
 179 }
 180
 181 void AliHLTTPCConfMapper::SetPointers()
 182 {
 183   //Check if there are not enough clusters to make a track in this sector
 184   //Can happen in pp events.
 185
 186   if(fClustersUnused < fMinPoints[fVertexConstraint])
 187     return;
 188
 189   //Allocate detector volumes
 190   if (fVolume==NULL) {
 191     LOG(AliHLTTPCLog::kDebug,"AliHLTTPCConfMapper::InitVolumes","Memory")<<AliHLTTPCLog::kDec<<
 192       "Allocating "<<fBounds*sizeof(AliHLTTPCConfMapContainer)<<" Bytes to fVolume"<<ENDLOG;
 193     fVolume = new AliHLTTPCConfMapContainer[fBounds];
 194   }
 195
 196   if (fRow==NULL) {
 197     LOG(AliHLTTPCLog::kDebug,"AliHLTTPCConfMapper::InitVolumes","Memory")<<AliHLTTPCLog::kDec<<
 198       "Allocating "<<fNumRowSegmentPlusOne*sizeof(AliHLTTPCConfMapContainer)<<" Bytes to fRow"<<ENDLOG;
 199     fRow = new AliHLTTPCConfMapContainer[fNumRowSegmentPlusOne];
 200   }
 201
 202   memset(fVolume,0,fBounds*sizeof(AliHLTTPCConfMapContainer));
 203   memset(fRow,0,fNumRowSegmentPlusOne*sizeof(AliHLTTPCConfMapContainer));
 204
 205   Float_t phiSlice = (fPhiMax-fPhiMin)/fNumPhiSegment;
 206   Float_t etaSlice = (fEtaMax-fEtaMin)/fNumEtaSegment;
 207
 208   Int_t volumeIndex;
 209   Int_t localcounter=0;
 210   assert((int)fHit.size()>=fClustersUnused);
 211   for(Int_t j=0; j<fClustersUnused; j++)
 212     {
 213       AliHLTTPCConfMapPoint *thisHit = &(fHit[j]);
 214
 215       thisHit->Setup(fVertex);
 216
 217       Int_t localrow = thisHit->GetPadRow();
 218
 219       if(localrow < fRowMin || localrow > fRowMax)
 220         continue;
 221
 222       //Get indexes:
 223       thisHit->SetPhiIndex((Int_t)((thisHit->GetPhi()-fPhiMin)/phiSlice +1));
 224
 225       if(thisHit->GetPhiIndex()<1 || thisHit->GetPhiIndex()>fNumPhiSegment)
 226         {
 227           //cout << "Phiindex: " << thisHit->phiIndex << " " << thisHit->GetPhi() << endl;
 228           fPhiHitsOutOfRange++;
 229           continue;
 230         }
 231
 232       thisHit->SetEtaIndex((Int_t)((thisHit->GetEta()-fEtaMin)/etaSlice + 1));
 233       if(thisHit->GetEtaIndex()<1 || thisHit->GetEtaIndex()>fNumEtaSegment)
 234         {
 235           //cout << "Etaindex: " << thisHit->etaIndex << " " << thisHit->GetEta() << endl;
 236           fEtaHitsOutOfRange++;
 237           continue;
 238         }
 239       localcounter++;
 240
 241       volumeIndex = (localrow-fRowMin)*fNumPhiEtaSegmentPlusOne +
 242                     thisHit->GetPhiIndex()*fNumEtaSegmentPlusOne+thisHit->GetEtaIndex();
 243
 244       if(fVolume[volumeIndex].first == NULL)
 245         fVolume[volumeIndex].first = (void *)thisHit;
 246       else
 247         ((AliHLTTPCConfMapPoint *)fVolume[volumeIndex].last)->SetNextVolumeHit(thisHit);
 248       fVolume[volumeIndex].last = (void *)thisHit;
 249
 250
 251       //set row pointers
 252       if(fRow[(localrow-fRowMin)].first == NULL)
 253         fRow[(localrow-fRowMin)].first = (void *)thisHit;
 254       else
 255         ((AliHLTTPCConfMapPoint *)(fRow[(localrow-fRowMin)].last))->SetNextRowHit(thisHit);
 256         fRow[(localrow-fRowMin)].last = (void *)thisHit;
 257     }
 258
 259   if(fClustersUnused>0 && localcounter==0)
 260     LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapper::SetPointers","Parameters")
 261       <<AliHLTTPCLog::kDec<<"No points passed to track finder, hits out of range: "
 262       <<fEtaHitsOutOfRange+fPhiHitsOutOfRange<<ENDLOG;
 263
 264   Int_t hits_accepted=fClustersUnused-(fEtaHitsOutOfRange+fPhiHitsOutOfRange);
 265   LOG(AliHLTTPCLog::kDebug,"AliHLTTPCConfMapper::SetPointers","Setup")
 266     <<"Setup finished, hits out of range: "<<fEtaHitsOutOfRange+fPhiHitsOutOfRange
 267     <<" hits accepted "<<hits_accepted<<ENDLOG;
 268 }
 269
 270 void AliHLTTPCConfMapper::MainVertexTrackingA()
 271 {
 272   //Tracking with vertex constraint.
 273
 274   if(!fParamSet[(Int_t)kTRUE])
 275     {
 276       LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapper::MainVertexTracking","Parameters")<<AliHLTTPCLog::kDec<<
 277         "Tracking parameters not set!"<<ENDLOG;
 278       return;
 279     }
 280
 281   Double_t initCpuTime,cpuTime;
 282   initCpuTime = CpuTime();
 283
 284   SetPointers();
 285   SetVertexConstraint(true);
 286   cpuTime = CpuTime() - initCpuTime;
 287   if(fBench)
 288     LOG(AliHLTTPCLog::kBenchmark,"AliHLTTPCConfMapper::MainVertexTrackingA","Timing")
 289       <<AliHLTTPCLog::kDec<<"Setup finished in "<<cpuTime*1000<<" ms"<<ENDLOG;
 290
 291 }
 292
 293 void AliHLTTPCConfMapper::MainVertexTrackingB()
 294 {
 295   //Tracking with vertex constraint.
 296
 297   if(!fParamSet[(Int_t)kTRUE])
 298     {
 299       LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapper::MainVertexTracking","Parameters")<<AliHLTTPCLog::kDec<<
 300         "Tracking parameters not set!"<<ENDLOG;
 301       return;
 302     }
 303   Double_t initCpuTime,cpuTime;
 304   initCpuTime = CpuTime();
 305
 306   ClusterLoop();
 307
 308   cpuTime = CpuTime() - initCpuTime;
 309   if(fBench)
 310     LOG(AliHLTTPCLog::kBenchmark,"AliHLTTPCConfMapper::MainVertexTrackingB","Timing")
 311       <<AliHLTTPCLog::kDec<<"Main Tracking finished in "<<cpuTime*1000<<" ms"<<ENDLOG;
 312 }
 313
 314 void AliHLTTPCConfMapper::MainVertexTracking()
 315 {
 316   //Tracking with vertex constraint.
 317
 318   if(!fParamSet[(Int_t)kTRUE])
 319     {
 320       LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapper::MainVertexTracking","Parameters")<<AliHLTTPCLog::kDec<<
 321         "Tracking parameters not set!"<<ENDLOG;
 322       return;
 323     }
 324
 325   Double_t initCpuTime,cpuTime;
 326   initCpuTime = CpuTime();
 327
 328   SetPointers();
 329
 330   SetVertexConstraint(true);
 331
 332   ClusterLoop();
 333
 334   cpuTime = CpuTime() - initCpuTime;
 335   if(fBench)
 336     LOG(AliHLTTPCLog::kBenchmark,"AliHLTTPCConfMapper::MainVertexTracking","Timing")<<AliHLTTPCLog::kDec<<
 337       "Tracking finished in "<<cpuTime*1000<<" ms"<<ENDLOG;
 338
 339   return;
 340 }
 341
 342 void AliHLTTPCConfMapper::NonVertexTracking()
 343 {
 344   //Tracking with no vertex constraint. This should be called after doing MainVertexTracking,
 345   //in order to do tracking on the remaining clusters.
 346   //The conformal mapping is now done with respect to the first cluster
 347   //assosciated with this track.
 348
 349   if(!fParamSet[(Int_t)kFALSE])
 350     {
 351       LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapper::NonVertexTracking","Parameters")<<AliHLTTPCLog::kDec<<
 352         "Tracking parameters not set!"<<ENDLOG;
 353       return;
 354     }
 355
 356   SetVertexConstraint(false);
 357
 358   SetPointers(); //To be able to do only nonvertextracking (more testing)
 359
 360   ClusterLoop();
 361   LOG(AliHLTTPCLog::kInformational,"AliHLTTPCConfMapper::NonVertexTracking","ntracks")<<AliHLTTPCLog::kDec<<
 362     "Number of nonvertex tracks found: "<<(fNTracks-fMainVertexTracks)<<ENDLOG;
 363   return;
 364 }
 365
 366 void AliHLTTPCConfMapper::MainVertexSettings(Int_t trackletlength, Int_t tracklength,
 367                                          Int_t rowscopetracklet, Int_t rowscopetrack,
 368                                          Double_t maxphi,Double_t maxeta)
 369 {
 370   //Settings for main vertex tracking. The cuts are:
 371   //TrackletLength:      #hits on segment, before trying to build a track
 372   //TrackLength:         Minimum hits on a track
 373   //RowScopeTracklet:    Row search range for segments
 374   //RowScopeTrack:       Row search range for tracks
 375
 376   SetTrackletLength(trackletlength,(Bool_t)true);
 377   SetRowScopeTracklet(rowscopetracklet, (Bool_t) true);
 378   SetRowScopeTrack(rowscopetrack, (Bool_t) true);
 379   SetMinPoints(tracklength,(Bool_t)true);
 380   fMaxPhi=maxphi;
 381   fMaxEta=maxeta;
 382   SetParamDone(kTRUE);
 383 }
 384
 385 void AliHLTTPCConfMapper::NonVertexSettings(Int_t trackletlength, Int_t tracklength,
 386                                         Int_t rowscopetracklet, Int_t rowscopetrack)
 387 {
 388   //set parameters for non-vertex tracking
 389   SetTrackletLength(trackletlength,(Bool_t)false);
 390   SetRowScopeTracklet(rowscopetracklet, (Bool_t)false);
 391   SetRowScopeTrack(rowscopetrack, (Bool_t)false);
 392   SetMinPoints(tracklength,(Bool_t)false);
 393   SetParamDone(kFALSE);
 394 }
 395
 396 void AliHLTTPCConfMapper::SetTrackCuts(Double_t hitChi2Cut, Double_t goodHitChi2, Double_t trackChi2Cut,Int_t maxdist,Bool_t vertexconstraint)
 397 {
 398   //Settings for tracks. The cuts are:
 399   //HitChi2Cut:     Maximum hit chi2
 400   //goodHitChi2:    Chi2 to stop look for next hit
 401   //trackChi2Cut:   Maximum track chi2
 402   //maxdist:        Maximum distance between two clusters when forming segments
 403
 404   SetHitChi2Cut(hitChi2Cut,vertexconstraint);
 405   SetGoodHitChi2(goodHitChi2,vertexconstraint);
 406   SetTrackChi2Cut(trackChi2Cut,vertexconstraint);
 407   SetMaxDist(maxdist,vertexconstraint);
 408 }
 409
 410 void AliHLTTPCConfMapper::SetTrackletCuts(Double_t maxangle,Double_t goodDist, Bool_t vc)
 411 {
 412   //Sets cuts of tracklets. Right now this is only:
 413   //maxangle:  Maximum angle when forming segments (if trackletlength > 2)
 414
 415   fGoodDist=goodDist;
 416   SetMaxAngleTracklet(maxangle, vc);
 417 }
 418
 419 void AliHLTTPCConfMapper::ClusterLoop()
 420 {
 421   //Loop over hits, starting at outermost padrow, and trying to build segments.
 422
 423   //Check if there are not enough clusters to make a track in this sector
 424   //Can happen in pp events.
 425   if(fClustersUnused < fMinPoints[fVertexConstraint])
 426     return;
 427
 428   Int_t rowsegm,lastrow = fRowMin + fMinPoints[fVertexConstraint];
 429   AliHLTTPCConfMapPoint *hit;
 430
 431   //Loop over rows, and try to create tracks from the hits.
 432   //Starts at the outermost row, and loops as long as a track can be build, due to length.
 433
 434   for(rowsegm = fRowMax; rowsegm >= lastrow; rowsegm--)
 435     {
 436       if(fRow[(rowsegm-fRowMin)].first && ((AliHLTTPCConfMapPoint*)fRow[(rowsegm-fRowMin)].first)->GetPadRow() < fRowMin + 1)
 437         break;
 438
 439       for(hit = (AliHLTTPCConfMapPoint*)fRow[(rowsegm-fRowMin)].first; hit!=0; hit=hit->GetNextRowHit())
 440         {
 441           if(hit->GetUsage() == true)
 442             continue;
 443           else
 444             CreateTrack(hit);
 445         }
 446     }
 447
 448   return;
 449 }
 450
 451
 452 void AliHLTTPCConfMapper::CreateTrack(AliHLTTPCConfMapPoint *hit)
 453 {
 454   //Tries to create a track from the initial hit given by ClusterLoop()
 455
 456   AliHLTTPCConfMapPoint *closesthit = NULL;
 457   AliHLTTPCConfMapTrack *track = NULL;
 458
 459   Int_t point;
 460   Int_t tracks = fNTracks;
 461   fNTracks++;
 462
 463   track = (AliHLTTPCConfMapTrack*)fTrack->NextTrack();
 464
 465   //reset hit parameters:
 466   track->Reset();
 467
 468   UInt_t *trackhitnumber = track->GetHitNumbers();
 469
 470   //set conformal coordinates if we are looking for non vertex tracks
 471   if(!fVertexConstraint)
 472     {
 473       hit->SetAllCoord(hit);
 474     }
 475
 476   //fill fit parameters of initial track:
 477   track->UpdateParam(hit); //here the number of hits is incremented.
 478   trackhitnumber[track->GetNumberOfPoints()-1] = hit->GetHitNumber();
 479
 480   Double_t dx,dy;
 481   //create tracklets:
 482
 483   for(point=1; point<fTrackletLength[fVertexConstraint]; point++)
 484     {
 485       if((closesthit = GetNextNeighbor(hit)))
 486         {//closest hit exist
 487
 488           //   Calculate track length in sz plane
 489           dx = ((AliHLTTPCConfMapPoint*)closesthit)->GetX() - ((AliHLTTPCConfMapPoint*)hit)->GetX();
 490           dy = ((AliHLTTPCConfMapPoint*)closesthit)->GetY() - ((AliHLTTPCConfMapPoint*)hit)->GetY();
 491           //track->fLength += (Double_t)sqrt ( dx * dx + dy * dy ) ;
 492           Double_t length = track->GetLength()+(Double_t)sqrt ( dx * dx + dy * dy );
 493           track->SetLength(length);
 494
 495           //closesthit->SetS(track->fLength);
 496           closesthit->SetS(track->GetLength());
 497
 498           //update fit parameters
 499           track->UpdateParam(closesthit);
 500           trackhitnumber[track->GetNumberOfPoints()-1] = closesthit->GetHitNumber();
 501
 502           hit = closesthit;
 503         }
 504       else
 505         {
 506           //closest hit does not exist:
 507           track->DeleteCandidate();
 508           fTrack->RemoveLast();
 509           fNTracks--;
 510           point = fTrackletLength[fVertexConstraint];
 511         }
 512     }
 513
 514   //tracklet is long enough to be extended to a track
 515   if(track->GetNumberOfPoints() == fTrackletLength[fVertexConstraint])
 516     {
 517
 518       track->SetProperties(true);
 519
 520       if(TrackletAngle(track) > fMaxAngleTracklet[fVertexConstraint])
 521         {//proof if the first points seem to be a beginning of a track
 522           track->SetProperties(false);
 523           track->DeleteCandidate();
 524           fTrack->RemoveLast();
 525           fNTracks--;
 526         }
 527
 528       else//good tracklet ->proceed, follow the trackfit
 529         {
 530           tracks++;
 531
 532           //define variables to keep the total chi:
 533           Double_t xyChi2 = track->GetChiSq1();
 534           Double_t szChi2 = track->GetChiSq2();
 535
 536           for(point = fTrackletLength[fVertexConstraint]; point <= fNumRowSegment; point++)
 537             {
 538               track->SetChiSq1(fHitChi2Cut[fVertexConstraint]);
 539               closesthit = GetNextNeighbor((AliHLTTPCConfMapPoint*)track->GetLastHit(),track);
 540
 541               if(closesthit)
 542                 {
 543                   //keep total chi:
 544                   Double_t lxyChi2 = track->GetChiSq1()-track->GetChiSq2();
 545                   xyChi2 += lxyChi2;
 546                   closesthit->SetXYChi2(lxyChi2);
 547
 548                   //update track length:
 549                   track->SetLength(closesthit->GetS());
 550                   szChi2 += track->GetChiSq2();
 551                   closesthit->SetSZChi2(track->GetChiSq2());
 552
 553                   track->UpdateParam(closesthit);
 554                   trackhitnumber[track->GetNumberOfPoints()-1] = closesthit->GetHitNumber();
 555
 556                   //add closest hit to track
 557                   closesthit->SetUsage(true);
 558                   closesthit->SetTrackNumber(tracks-1);
 559
 560                 }//closesthit
 561
 562               else
 563                 {
 564                   //closest hit does not exist
 565                   point = fNumRowSegment; //continue with next hit in segment
 566                 }//else
 567
 568             }//create tracks
 569
 570           //store track chi2:
 571           track->SetChiSq1(xyChi2);
 572           track->SetChiSq2(szChi2);
 573           Double_t normalizedchi2 = (track->GetChiSq1()+track->GetChiSq2())/track->GetNumberOfPoints();
 574
 575           //remove tracks with not enough points already now
 576           if(track->GetNumberOfPoints() < fMinPoints[fVertexConstraint] || normalizedchi2 > fTrackChi2Cut[fVertexConstraint])
 577             {
 578               track->SetProperties(false);
 579               fNTracks--;
 580               track->DeleteCandidate();
 581               fTrack->RemoveLast();
 582               tracks--;
 583             }
 584
 585           else
 586             {
 587               fClustersUnused -= track->GetNumberOfPoints();
 588               track->ComesFromMainVertex(fVertexConstraint);
 589               //mark track as main vertex track or not
 590               track->SetSector(2); //only needed for testing purposes.
 591               track->SetRowRange(fRowMin,fRowMax);
 592
 593               if(fVertexConstraint)
 594                 fMainVertexTracks++;
 595             }
 596
 597         }//good tracklet
 598
 599     }
 600
 601   return;
 602 }
 603
 604 AliHLTTPCConfMapPoint *AliHLTTPCConfMapper::GetNextNeighbor(AliHLTTPCConfMapPoint *starthit,
 605                                           AliHLTTPCConfMapTrack *track)
 606 {
 607   //When forming segments: Finds closest hit to input hit
 608   //When forming tracks: Find closest hit to track fit.
 609
 610   Double_t dist,closestdist = fMaxDist[fVertexConstraint];
 611
 612   AliHLTTPCConfMapPoint *hit = NULL;
 613   AliHLTTPCConfMapPoint *closesthit = NULL;
 614
 615   Int_t subrowsegm;
 616   Int_t subphisegm;
 617   Int_t subetasegm;
 618   Int_t volumeIndex;
 619   Int_t testhit;
 620
 621   Int_t maxrow = starthit->GetPadRow()-1;
 622   Int_t minrow;
 623
 624   if(track) //finding hit close to trackfit
 625     {
 626       minrow = starthit->GetPadRow()-fRowScopeTrack[fVertexConstraint];
 627     }
 628   else
 629     {
 630       minrow = starthit->GetPadRow()-fRowScopeTracklet[fVertexConstraint];
 631     }
 632
 633   //make a smart loop
 634   Int_t loopeta[25] = {0,0,0,-1,-1,-1,1,1,1, 0,0,-1,-1,1,1,-2,-2,-2,-2,-2,2,2,2,2,2};
 635   Int_t loopphi[25] = {0,-1,1,0,-1,1,0,-1,1, -2,2,-2,2,-2,2,-2,-1,0,1,2,-2,-1,0,1,2};
 636
 637   if(minrow < fRowMin)
 638     minrow = fRowMin;
 639   if(maxrow < fRowMin)
 640     return 0;  //reached the last padrow under consideration
 641
 642   else
 643     {
 644       //loop over sub rows
 645       for(subrowsegm=maxrow; subrowsegm>=minrow; subrowsegm--)
 646         {
 647           //loop over subsegments, in the order defined above.
 648           for(Int_t i=0; i<9; i++)
 649             {
 650               subphisegm = starthit->GetPhiIndex() + loopphi[i];
 651
 652               if(subphisegm < 0 || subphisegm >= fNumPhiSegment)
 653                 continue;
 654               /*
 655                 if(subphisegm<0)
 656                 subphisegm += fNumPhiSegment;
 657
 658                 else if(subphisegm >=fNumPhiSegment)
 659                 subphisegm -= fNumPhiSegment;
 660               */
 661               //loop over sub eta segments
 662
 663               subetasegm = starthit->GetEtaIndex() + loopeta[i];
 664
 665               if(subetasegm < 0 || subetasegm >=fNumEtaSegment)
 666                 continue;//segment exceeds bounds->skip it
 667
 668               //loop over hits in this sub segment:
 669               volumeIndex=(subrowsegm-fRowMin)*fNumPhiEtaSegmentPlusOne +
 670                 subphisegm*fNumEtaSegmentPlusOne + subetasegm;
 671
 672               if(volumeIndex<0)
 673                 {//debugging
 674                   LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapper::GetNextNeighbor","Memory")<<AliHLTTPCLog::kDec<<
 675                     "VolumeIndex error "<<volumeIndex<<ENDLOG;
 676                 }
 677
 678               assert(fVolume!=NULL);
 679               for(hit = (AliHLTTPCConfMapPoint*)fVolume[volumeIndex].first;
 680                   hit!=0; hit = hit->GetNextVolumeHit())
 681                 {
 682
 683                   if(!hit->GetUsage())
 684                     {//hit was not used before
 685
 686                       //set conformal mapping if looking for nonvertex tracks:
 687                       if(!fVertexConstraint)
 688                         {
 689                           hit->SetAllCoord(starthit);
 690                         }
 691
 692                       if(track)//track search - look for nearest neighbor to extrapolated track
 693                         {
 694                           if (fVertexConstraint) {
 695                             if(!VerifyRange(starthit,hit))
 696                               continue;
 697                           }
 698                           testhit = EvaluateHit(starthit,hit,track);
 699
 700                           if(testhit == 0)//chi2 not good enough, keep looking
 701                             continue;
 702                           else if(testhit==2)//chi2 good enough, return it
 703                             return hit;
 704                           else
 705                             closesthit = hit;//chi2 acceptable, but keep looking
 706
 707                         }//track search
 708
 709                       else //tracklet search, look for nearest neighbor
 710                         {
 711
 712                           if((dist=CalcDistance(starthit,hit)) < closestdist)
 713                             {
 714                               if (fVertexConstraint) {
 715                                 if(!VerifyRange(starthit,hit))
 716                                   continue;
 717                               }
 718                               closestdist = dist;
 719                               closesthit = hit;
 720
 721                               //if this hit is good enough, return it:
 722                               if(closestdist < fGoodDist)
 723                                 return closesthit;
 724                             }
 725                           else
 726                             continue;//sub hit was farther away than a hit before
 727
 728                         }//tracklet search
 729
 730                     }//hit not used before
 731
 732                   else continue; //sub hit was used before
 733
 734                 }//loop over hits in sub segment
 735
 736             }//loop over sub segments
 737
 738         }//loop over subrows
 739
 740     }//else
 741
 742   //closest hit found:
 743   if(closesthit)// && closestdist < mMaxDist)
 744     return closesthit;
 745   else
 746     return 0;
 747 }
 748
 749 Int_t AliHLTTPCConfMapper::EvaluateHit(AliHLTTPCConfMapPoint *starthit,AliHLTTPCConfMapPoint *hit,AliHLTTPCConfMapTrack *track)
 750 {
 751   //Check if space point gives a fit with acceptable chi2.
 752
 753   Double_t temp,dxy,lchi2,dx,dy,slocal,dsz,lszChi2;
 754   temp = (track->GetA2Xy()*hit->GetXprime()-hit->GetYprime()+track->GetA1Xy());
 755   dxy = temp*temp/(track->GetA2Xy()*track->GetA2Xy() + 1.);
 756
 757   //Calculate chi2
 758   lchi2 = (dxy*hit->GetXYWeight());
 759
 760   if(lchi2 > track->GetChiSq1())//chi2 was worse than before.
 761     return 0;
 762
 763   //calculate s and the distance hit-line
 764   dx = starthit->GetX()-hit->GetX();
 765   dy = starthit->GetY()-hit->GetY();
 766   //slocal = track->fLength+sqrt(dx*dx+dy*dy);
 767   slocal = track->GetLength()+sqrt(dx*dx+dy*dy);
 768
 769   temp = (track->GetA2Sz()*slocal-hit->GetZ()+track->GetA1Sz());
 770   dsz = temp*temp/(track->GetA2Sz()*track->GetA2Sz()+1);
 771
 772   //calculate chi2
 773   lszChi2 = dsz*hit->GetZWeight();
 774   lchi2 += lszChi2;
 775
 776
 777   //check whether chi2 is better than previous one:
 778   if(lchi2 < track->GetChiSq1())
 779     {
 780       track->SetChiSq1(lchi2);
 781       track->SetChiSq2(lszChi2);
 782
 783       hit->SetS(slocal);
 784
 785       //if chi2 good enough, stop here:
 786       if(lchi2 < fGoodHitChi2[fVertexConstraint])
 787         return 2;
 788
 789       return 1;
 790     }
 791
 792   return 0;
 793
 794 }
 795
 796 Double_t AliHLTTPCConfMapper::CalcDistance(const AliHLTTPCConfMapPoint *hit1,const AliHLTTPCConfMapPoint *hit2) const
 797 {
 798   //Return distance between two clusters, defined by Pablo
 799
 800   Double_t phidiff = fabs( hit1->GetPhi() - hit2->GetPhi() );
 801   if (phidiff > AliHLTTPCTransform::Pi()) phidiff = AliHLTTPCTransform::TwoPi() - phidiff;
 802
 803   return AliHLTTPCTransform::ToDeg()*fabs((Float_t)((hit1->GetPadRow() - hit2->GetPadRow()) *
 804          (phidiff + fabs( hit1->GetEta() - hit2->GetEta()))));
 805 }
 806
 807 Bool_t AliHLTTPCConfMapper::VerifyRange(const AliHLTTPCConfMapPoint *hit1,const AliHLTTPCConfMapPoint *hit2) const
 808 {
 809   //Check if the hit are within reasonable range in phi and eta
 810   Double_t dphi,deta;//maxphi=0.1,maxeta=0.1;
 811   dphi = fabs(hit1->GetPhi() - hit2->GetPhi());
 812   if(dphi > AliHLTTPCTransform::Pi()) dphi = fabs(AliHLTTPCTransform::TwoPi() - dphi);
 813   if(dphi > fMaxPhi) return false;
 814
 815   deta = fabs(hit1->GetEta() - hit2->GetEta());
 816   if(deta > fMaxEta) return false;
 817
 818   return true;
 819
 820 }
 821
 822 Double_t AliHLTTPCConfMapper::TrackletAngle(AliHLTTPCConfMapTrack *track,Int_t n) const
 823 {
 824   // Returns the angle 'between' the last three points (started at point number n) on this track.
 825
 826   if(n > track->GetNumberOfPoints())
 827     n = track->GetNumberOfPoints();
 828
 829   if(n<3)
 830     return 0;
 831
 832   Double_t x1[2];
 833   Double_t x2[2];
 834   Double_t x3[2];
 835   Double_t angle1,angle2;
 836   Int_t counter=0;
 837   for(track->StartLoop(); track->LoopDone(); track->GetNextHit())
 838     {
 839       AliHLTTPCConfMapPoint *p = (AliHLTTPCConfMapPoint*)track->GetCurrentHit();
 840       if( (n-1) == counter)
 841         {
 842           x1[0] = p->GetX();
 843           x1[1] = p->GetY();
 844         }
 845       else if( (n-2) == counter)
 846         {
 847           x2[0] = p->GetX();
 848           x2[1] = p->GetY();
 849         }
 850       else if( (n-3) == counter)
 851         {
 852           x3[0] = p->GetX();
 853           x3[1] = p->GetY();
 854         }
 855       counter++;
 856     }
 857
 858   angle1 = atan2(x2[1]-x3[1],x2[0]-x3[0]);
 859   angle2 = atan2(x1[1]-x2[1],x1[0]-x2[0]);
 860
 861   return fabs(angle1-angle2);
 862
 863   /*
 864     Double_t x1[2];
 865   Double_t x2[2];
 866   Double_t angle1,angle2;
 867   TObjArray *hits = track->GetHits();
 868
 869   if (n > track->GetNumberOfPoints()) {
 870     n = track->GetNumberOfPoints();
 871   }
 872
 873   if (n<3)
 874     return 0;
 875
 876
 877   x1[0] = ((AliHLTTPCConfMapPoint *)hits->At(n-2))->GetX() - ((AliHLTTPCConfMapPoint *)hits->At(n-3))->GetX();
 878   x1[1] = ((AliHLTTPCConfMapPoint *)hits->At(n-2))->GetY() - ((AliHLTTPCConfMapPoint *)hits->At(n-3))->GetY();
 879
 880   x2[0] = ((AliHLTTPCConfMapPoint *)hits->At(n-1))->GetX() - ((AliHLTTPCConfMapPoint *)hits->At(n-2))->GetX();
 881   x2[1] = ((AliHLTTPCConfMapPoint *)hits->At(n-1))->GetY() - ((AliHLTTPCConfMapPoint *)hits->At(n-2))->GetY();
 882
 883   angle1 = atan2(x1[1],x1[0]);
 884   angle2 = atan2(x2[1],x1[0]);
 885   return fabs(angle1-angle2);
 886   */
 887 }
 888
 889 Int_t AliHLTTPCConfMapper::FillTracks()
 890 {
 891   //Fill track parameters. Which basically means do a fit of helix in real space,
 892   //which should be done in order to get nice tracks.
 893
 894   Int_t numoftracks = fNTracks;
 895   if(fNTracks == 0)
 896     {
 897       LOG(AliHLTTPCLog::kDebug,"AliHLTTPCConfMapper::FillTracks","fNTracks")<<AliHLTTPCLog::kDec<<
 898         "No tracks found!!"<<ENDLOG;
 899       return 0;
 900     }
 901
 902   LOG(AliHLTTPCLog::kInformational,"AliHLTTPCConfMapper::FillTracks","fNTracks")<<AliHLTTPCLog::kDec<<
 903     "Number of found tracks: "<<fNTracks<<ENDLOG;
 904
 905   //  fTrack->Sort();
 906   for(Int_t i=0; i<numoftracks; i++)
 907     {
 908       AliHLTTPCConfMapTrack *track = (AliHLTTPCConfMapTrack*)fTrack->GetTrack(i);
 909       track->Fill(fVertex,fMaxDca);
 910     }
 911   return 1;
 912 }
 913
 914 Double_t AliHLTTPCConfMapper::CpuTime()
 915 {
 916   //Return the Cputime in seconds.
 917  struct timeval tv;
 918  gettimeofday( &tv, NULL );
 919  return tv.tv_sec+(((Double_t)tv.tv_usec)/1000000.);
 920 }