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