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