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