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