Some additional changes related to the previous changes. AliL3Transform
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughTrack.cxx
CommitLineData
b1886074 1//$Id$
2
3// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4//*-- Copyright &copy ASV
d712b5b8 5
6#include <math.h>
7
8#include "AliL3Logging.h"
9#include "AliL3Transform.h"
10#include "AliL3Defs.h"
11#include "AliL3HoughTrack.h"
12
b1886074 13//_____________________________________________________________
14// AliL3HoughTrack
15//
16// Track class for Hough tracklets
17
d712b5b8 18ClassImp(AliL3HoughTrack)
19
20
21AliL3HoughTrack::AliL3HoughTrack()
22{
23 //Constructor
24
25 fWeight = 0;
26 fMinDist=0;
d712b5b8 27 fDLine = 0;
28 fPsiLine = 0;
29 fIsHelix = true;
30 fEtaIndex = -1;
31}
32
33
34AliL3HoughTrack::~AliL3HoughTrack()
35{
237d3f5c 36
d712b5b8 37}
38
39void 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());
ae4d5136 46 SetTgl(tpt->GetTgl());
d712b5b8 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());
b1886074 54 SetSlice(tpt->GetSlice());
242c143e 55 SetNHits(1);
d712b5b8 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
68Int_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
76void AliL3HoughTrack::SetEta(Double_t f)
77{
b1886074 78 //Set eta, and calculate fTanl, which is the tan of dipangle
79
d712b5b8 80 fEta = f;
81 Double_t theta = 2*atan(exp(-1.*fEta));
8b007a2d 82 Double_t dipangle = Pi/2 - theta;
83 Double_t tgl = tan(dipangle);
d712b5b8 84 SetTgl(tgl);
85}
86
87void AliL3HoughTrack::CalculateHelix()
88{
89 return;
90}
91
92void 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
8b007a2d 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
d712b5b8 107 Double_t radius = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
108
109 //Get the track parameters
b1886074 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 ) ;
d712b5b8 116 Double_t tPhi0 = GetPsi() + GetCharge() * 0.5 * pi / fabs(GetCharge()) ;
b1886074 117 Double_t xc = GetCenterX();//x0 - rc * cos(tPhi0) ;
118 Double_t yc = GetCenterY();//y0 - rc * sin(tPhi0) ;
d712b5b8 119
120 //Check helix and cylinder intersect
121 Double_t fac1 = xc*xc + yc*yc ;
122 Double_t sfac = sqrt( fac1 ) ;
b1886074 123
d712b5b8 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 ;
8b007a2d 141
d712b5b8 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
d712b5b8 155 SetR0(radius);
156 SetPhi0(phi);
157 SetZ0(z);
158 SetPsi(tPsi);
159 SetFirstPoint(xyz[0],xyz[1],z);
8b007a2d 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);
d712b5b8 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]);
8b007a2d 167 //printf("last point %f %f %f\n",xyz[0],xyz[1],xyz[2]);
d712b5b8 168}
169
170void 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);
b1886074 177 Double_t pt = fabs(BFACT*BField/kappa);
d712b5b8 178 SetPt(pt);
179 Double_t radius = 1/fabs(kappa);
180 SetRadius(radius);
181 SetFirstPoint(0,0,0);
8b007a2d 182 SetPsi(phi); //Psi = Phi when first point is vertex
183 SetR0(0);
d712b5b8 184 Double_t charge = -1.*kappa;
8b007a2d 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) ;
d712b5b8 190 SetCenterX(xc);
191 SetCenterY(yc);
242c143e 192 SetNHits(1); //just for the trackarray IO
d712b5b8 193 fIsHelix = true;
8b007a2d 194
242c143e 195
d712b5b8 196}
197
198void 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
237d3f5c 205 /*
d712b5b8 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;
237d3f5c 214 */
d712b5b8 215}
216
217void 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
228void AliL3HoughTrack::GetLineCrossingPoint(Int_t padrow,Double_t *xy)
229{
230
237d3f5c 231 /*
d712b5b8 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;
237d3f5c 245 */
d712b5b8 246}
247
248/*
249Double_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
290Bool_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
327Bool_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*/