2978901102522706a4cbc2835dae11234154ae6d
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughTrack.cxx
1 //$Id$
2
3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright &copy ASV 
5
6 #include <math.h>
7
8 #include "AliL3Logging.h"
9 #include "AliL3HoughTrack.h"
10 #include "AliL3Transform.h"
11
12 //_____________________________________________________________
13 // AliL3HoughTrack
14 //
15 // Track class for Hough tracklets
16
17 ClassImp(AliL3HoughTrack)
18
19
20 AliL3HoughTrack::AliL3HoughTrack()
21 {
22   //Constructor
23   
24   fWeight = 0;
25   fMinDist=0;
26   fDLine = 0;
27   fPsiLine = 0;
28   fIsHelix = true;
29   fEtaIndex = -1;
30   fEta = 0;
31   
32 }
33
34
35 AliL3HoughTrack::~AliL3HoughTrack()
36 {
37   
38 }
39
40 void AliL3HoughTrack::Set(AliL3Track *track)
41 {
42   
43   AliL3HoughTrack *tpt = (AliL3HoughTrack*)track;
44   SetTrackParameters(tpt->GetKappa(),tpt->GetPsi(),tpt->GetWeight());
45   SetEtaIndex(tpt->GetEtaIndex());
46   SetEta(tpt->GetEta());
47   SetTgl(tpt->GetTgl());
48   SetPsi(tpt->GetPsi());
49   SetCenterX(tpt->GetCenterX());
50   SetCenterY(tpt->GetCenterY());
51   SetFirstPoint(tpt->GetFirstPointX(),tpt->GetFirstPointY(),tpt->GetFirstPointZ());
52   SetLastPoint(tpt->GetLastPointX(),tpt->GetLastPointY(),tpt->GetLastPointZ());
53   SetCharge(tpt->GetCharge());
54   SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
55   SetSlice(tpt->GetSlice());
56   SetNHits(1);
57   return;
58
59   fWeight = tpt->GetWeight();
60   fDLine = tpt->GetDLine();
61   fPsiLine = tpt->GetPsiLine();
62   SetNHits(tpt->GetWeight());
63   SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
64   fIsHelix = false;
65
66
67 }
68
69 Int_t AliL3HoughTrack::Compare(const AliL3Track *tpt) const
70 {
71   AliL3HoughTrack *track = (AliL3HoughTrack*)tpt;
72   if(track->GetWeight() < GetWeight()) return 1;
73   if(track->GetWeight() > GetWeight()) return -1;
74   return 0;
75 }
76
77 void AliL3HoughTrack::SetEta(Double_t f)
78 {
79   //Set eta, and calculate fTanl, which is the tan of dipangle
80
81   fEta = f;
82   Double_t theta = 2*atan(exp(-1.*fEta));
83   Double_t dipangle = AliL3Transform::Pi()/2 - theta;
84   Double_t tgl = tan(dipangle);
85   SetTgl(tgl);
86 }
87
88 void AliL3HoughTrack::CalculateHelix()
89 {
90   return;
91 }
92
93 void AliL3HoughTrack::UpdateToFirstRow()
94 {
95   //Update the track parameters to the point where track cross
96   //its first padrow.
97   
98   //Get the crossing point with the first padrow:
99   Float_t xyz[3];
100   if(!GetCrossingPoint(GetFirstRow(),xyz))
101     LOG(AliL3Log::kWarning,"AliL3HoughTrack::UpdateToFirstRow()","Track parameters")
102       <<AliL3Log::kDec<<"Track does not cross padrow "<<GetFirstRow()<<" centerx "
103       <<GetCenterX()<<" centery "<<GetCenterY()<<" Radius "<<GetRadius()<<" tgl "<<GetTgl()<<ENDLOG;
104   
105   //printf("Track with eta %f tgl %f crosses at x %f y %f z %f on padrow %d\n",GetEta(),GetTgl(),xyz[0],xyz[1],xyz[2],GetFirstRow());
106   //printf("Before: first %f %f %f tgl %f center %f %f charge %d\n",GetFirstPointX(),GetFirstPointY(),GetFirstPointZ(),GetTgl(),GetCenterX(),GetCenterY(),GetCharge());
107   
108   Double_t radius = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
109
110   //Get the track parameters
111   
112   /*
113     Double_t x0    = GetR0() * cos(GetPhi0()) ;
114     Double_t y0    = GetR0() * sin(GetPhi0()) ;
115   */
116   Double_t rc    = GetRadius();//fabs(GetPt()) / ( BFACT * AliL3Transform::GetBField() )  ;
117   Double_t tPhi0 = GetPsi() + GetCharge() * 0.5 * pi / fabs(GetCharge()) ;
118   Double_t xc    = GetCenterX();//x0 - rc * cos(tPhi0) ;
119   Double_t yc    = GetCenterY();//y0 - rc * sin(tPhi0) ;
120   
121   //Check helix and cylinder intersect
122   Double_t fac1 = xc*xc + yc*yc ;
123   Double_t sfac = sqrt( fac1 ) ;
124   
125   if ( fabs(sfac-rc) > radius || fabs(sfac+rc) < radius ) {
126     LOG(AliL3Log::kError,"AliL3HoughTrack::UpdateToFirstRow","Tracks")<<AliL3Log::kDec<<
127       "Track does not intersect"<<ENDLOG;
128     return;
129   }
130   
131   //Find intersection
132   Double_t fac2   = (radius*radius + fac1 - rc*rc) / (2.00 * radius * sfac ) ;
133   Double_t phi    = atan2(yc,xc) + GetCharge()*acos(fac2) ;
134   Double_t td     = atan2(radius*sin(phi) - yc,radius*cos(phi) - xc) ;
135   
136   //Intersection in z
137   if ( td < 0 ) td = td + 2. * pi ;
138   Double_t deltat = fmod((-GetCharge()*td + GetCharge()*tPhi0),2*pi) ;
139   if ( deltat < 0.      ) deltat += 2. * pi ;
140   if ( deltat > 2.*pi ) deltat -= 2. * pi ;
141   Double_t z = GetZ0() + rc * GetTgl() * deltat ;
142   
143   
144   Double_t xExtra = radius * cos(phi) ;
145   Double_t yExtra = radius * sin(phi) ;
146   
147   Double_t tPhi = atan2(yExtra-yc,xExtra-xc);
148   
149   //if ( tPhi < 0 ) tPhi += 2. * M_PI ;
150   
151   Double_t tPsi = tPhi - GetCharge() * 0.5 * pi / fabs(GetCharge()) ;
152   if ( tPsi > 2. * pi ) tPsi -= 2. * pi ;
153   if ( tPsi < 0.        ) tPsi += 2. * pi ;
154   
155   //And finally, update the track parameters
156   SetR0(radius);
157   SetPhi0(phi);
158   SetZ0(z);
159   SetPsi(tPsi);
160   SetFirstPoint(xyz[0],xyz[1],z);
161   //printf("After: first %f %f %f tgl %f center %f %f charge %d\n",GetFirstPointX(),GetFirstPointY(),GetFirstPointZ(),GetTgl(),GetCenterX(),GetCenterY(),GetCharge());
162   
163   //printf("First point set %f %f %f\n",xyz[0],xyz[1],z);
164   
165   //Also, set the coordinates of the point where track crosses last padrow:
166   GetCrossingPoint(GetLastRow(),xyz);
167   SetLastPoint(xyz[0],xyz[1],xyz[2]);
168   //printf("last point %f %f %f\n",xyz[0],xyz[1],xyz[2]);
169 }
170
171 void AliL3HoughTrack::SetTrackParameters(Double_t kappa,Double_t eangle,Int_t weight)
172 {
173
174   fWeight = weight;
175   fMinDist = 100000;
176   SetKappa(kappa);
177   Double_t pt = fabs(BFACT*AliL3Transform::GetBField()/kappa);
178   SetPt(pt);
179   Double_t radius = 1/fabs(kappa);
180   SetRadius(radius);
181   SetFirstPoint(0,0,0);
182   SetPsi(eangle); //Psi = emission angle when first point is vertex
183   SetPhi0(0);     //not defined for vertex reference point
184   SetR0(0);
185   Double_t charge = -1.*kappa;
186   SetCharge((Int_t)copysign(1.,charge));
187   Double_t trackPhi0 = GetPsi() + charge*0.5*AliL3Transform::Pi()/fabs(charge);
188   Double_t xc = GetFirstPointX() - GetRadius() * cos(trackPhi0) ;
189   Double_t yc = GetFirstPointY() - GetRadius() * sin(trackPhi0) ;
190   SetCenterX(xc);
191   SetCenterY(yc);
192   SetNHits(1); //just for the trackarray IO
193   fIsHelix = true;
194 }
195
196 void AliL3HoughTrack::SetLineParameters(Double_t psi,Double_t D,Int_t weight,Int_t *rowrange,Int_t ref_row)
197 {
198   //Initialize a track piece, not yet a track
199   //Used in case of straight line transformation
200
201   //Transform line parameters to coordinate system of slice:
202   
203   /*
204   D = D + fTransform->Row2X(ref_row)*cos(psi);
205
206   fDLine = D;
207   fPsiLine = psi;
208   fWeight = weight;
209   SetNHits(weight);
210   SetRowRange(rowrange[0],rowrange[1]);
211   fIsHelix = false;
212   */
213 }
214
215 void AliL3HoughTrack::SetBestMCid(Int_t mcid,Double_t min_dist)
216 {
217   
218   if(min_dist < fMinDist)
219     {
220       fMinDist = min_dist;
221       SetMCid(mcid);
222     }
223   
224 }
225
226 void AliL3HoughTrack::GetLineCrossingPoint(Int_t padrow,Double_t *xy)
227 {
228   
229   /*
230   if(fIsHelix)
231     {
232       printf("AliL3HoughTrack::GetLineCrossingPoint : Track is not a line\n");
233       return;
234     }
235
236   Double_t xhit = fTransform->Row2X(padrow);
237   Double_t a = -1/tan(fPsiLine);
238   Double_t b = fDLine/sin(fPsiLine);
239   
240   Double_t yhit = a*xhit + b;
241   xy[0] = xhit;
242   xy[1] = yhit;
243   */
244 }
245
246 /*
247 Double_t AliL3HoughTrack::GetCrossingAngle(Int_t padrow)
248 {
249   //Calculate the crossing angle between track and given padrow.
250
251   if(!fIsHelix)
252     {
253       printf("AliL3HoughTrack::GetCrossingAngle : Track is not a helix\n");
254       return 0;
255     }
256
257   if(!IsLocal())
258     {
259       printf("Track is not given in local coordinates\n");
260       return 0;
261     }
262
263   Float_t xyz[3];
264   if(!GetCrossingPoint(padrow,xyz))
265     printf("AliL3HoughTrack::GetCrossingPoint : Track does not cross line!!\n");
266   
267   
268   //Convert center of curvature to local coordinates:
269   //Float_t xyz_coc[3] = {GetCenterX(),GetCenterY(),0};
270   //fTransform->Global2Local(xyz_coc,slice);
271   
272   //Take the dot product of the tangent vector of the track, and
273   //vector perpendicular to the padrow.
274   
275   Double_t tangent[2];
276   //tangent[1] = (xyz[0] - xyz_coc[0])/GetRadius();
277   //tangent[0] = -1.*(xyz[1] - xyz_coc[1])/GetRadius();
278   tangent[1] = (xyz[0] - GetCenterX())/GetRadius();
279   tangent[0] = -1.*(xyz[1] - GetCenterY())/GetRadius();
280
281   Double_t perp_padrow[2] = {1,0}; //locally in slice
282
283   Double_t cos_beta = fabs(tangent[0]*perp_padrow[0] + tangent[1]*perp_padrow[1]);
284   return acos(cos_beta);
285   
286 }
287
288 Bool_t AliL3HoughTrack::GetCrossingPoint(Int_t padrow,Float_t *xyz)
289 {
290   //Assumes the track is given in local coordinates
291
292   if(!fIsHelix)
293     {
294       printf("AliL3HoughTrack::GetCrossingPoint : Track is not a helix\n");
295       return 0;
296     }
297     
298
299   if(!IsLocal())
300     {
301       printf("GetCrossingPoint: Track is given on global coordinates\n");
302       return false;
303     }
304   
305   Double_t xHit = fTransform->Row2X(padrow);
306
307   //xyz[0] = fTransform->Row2X(padrow);
308   xyz[0] = xHit;
309   Double_t aa = (xHit - GetCenterX())*(xHit - GetCenterX());
310   Double_t r2 = GetRadius()*GetRadius();
311   if(aa > r2)
312     return false;
313
314   Double_t aa2 = sqrt(r2 - aa);
315   Double_t y1 = GetCenterY() + aa2;
316   Double_t y2 = GetCenterY() - aa2;
317   xyz[1] = y1;
318   if(fabs(y2) < fabs(y1)) xyz[1] = y2;
319   xyz[2] = 0; //only consider transverse plane
320   
321   return true;
322 }
323
324
325 Bool_t AliL3HoughTrack::GetCrossingPoint(Int_t slice,Int_t padrow,Float_t *xyz)
326 {
327   //Calculate the crossing point in transverse plane of this track and given 
328   //padrow (y = a*x + b). Point is given in local coordinates in slice.
329   //Assumes the track is given in global coordinates.
330
331   if(!fIsHelix)
332     {
333       printf("AliL3HoughTrack::GetCrossingPoint : Track is not a helix\n");
334       return 0;
335     }
336
337   
338   if(IsLocal())
339     {
340       printf("GetCrossingPoint: Track is given in local coordintes!!!\n");
341       return false;
342     }
343
344   Double_t padrowradii = fTransform->Row2X(padrow);
345
346
347   Float_t rotation_angle = (slice*20)*ToRad;
348   
349   Float_t cs,sn;
350   cs = cos(rotation_angle);
351   sn = sin(rotation_angle);
352
353   Double_t a = -1.*cs/sn;
354   Double_t b = padrowradii/sn;
355
356   Double_t ycPrime = GetCenterY() - b ;
357   Double_t aa = ( 1. + a * a ) ;
358   Double_t bb = -2. * ( GetCenterX() + a * ycPrime ) ;
359   Double_t cc = ( GetCenterX() * GetCenterX() + ycPrime * ycPrime - GetRadius() * GetRadius() ) ;
360
361   Double_t racine = bb * bb - 4. * aa * cc ;
362   if ( racine < 0 ) return false ;
363   Double_t rootRacine = sqrt(racine) ;
364   
365   Double_t oneOverA = 1./aa;
366 //
367 //   First solution
368 //
369    Double_t x1 = 0.5 * oneOverA * ( -1. * bb + rootRacine ) ; 
370    Double_t y1 = a * x1 + b ;
371    Double_t r1 = sqrt(x1*x1+y1*y1);
372 //
373 //   Second solution
374 //
375    Double_t x2 = 0.5 * oneOverA * ( -1. * bb - rootRacine ) ; 
376    Double_t y2 = a * x2 + b ;
377    Double_t r2 = sqrt(x2*x2+y2*y2);
378 //
379 //    Choose close to (0,0) 
380 //
381    Double_t xHit ;
382    Double_t yHit ;
383    if ( r1 < r2 ) {
384       xHit = x1 ;
385       yHit = y1 ;
386    }
387    else {
388       xHit = x2 ;
389       yHit = y2 ;
390    }
391   
392    xyz[0] = xHit;
393    xyz[1] = yHit;
394    xyz[2] = 0;
395    
396    fTransform->Global2Local(xyz,slice);
397    
398    return true;
399 }
400 */