]>
Commit | Line | Data |
---|---|---|
b1886074 | 1 | //$Id$ |
2 | ||
3 | // Author: Anders Vestbo <mailto:vestbo@fi.uib.no> | |
4 | //*-- Copyright © 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 | 18 | ClassImp(AliL3HoughTrack) |
19 | ||
20 | ||
21 | AliL3HoughTrack::AliL3HoughTrack() | |
22 | { | |
23 | //Constructor | |
24 | ||
25 | fWeight = 0; | |
26 | fMinDist=0; | |
27 | fTransform = new AliL3Transform(); | |
28 | fDLine = 0; | |
29 | fPsiLine = 0; | |
30 | fIsHelix = true; | |
31 | fEtaIndex = -1; | |
32 | } | |
33 | ||
34 | ||
35 | AliL3HoughTrack::~AliL3HoughTrack() | |
36 | { | |
37 | //Destructor | |
38 | if(fTransform) | |
39 | delete fTransform; | |
40 | } | |
41 | ||
42 | void AliL3HoughTrack::Set(AliL3Track *track) | |
43 | { | |
44 | ||
45 | AliL3HoughTrack *tpt = (AliL3HoughTrack*)track; | |
46 | SetTrackParameters(tpt->GetKappa(),tpt->GetPhi0(),tpt->GetWeight()); | |
47 | SetEtaIndex(tpt->GetEtaIndex()); | |
48 | SetEta(tpt->GetEta()); | |
49 | SetPsi(tpt->GetPsi()); | |
50 | SetCenterX(tpt->GetCenterX()); | |
51 | SetCenterY(tpt->GetCenterY()); | |
52 | SetFirstPoint(tpt->GetFirstPointX(),tpt->GetFirstPointY(),tpt->GetFirstPointZ()); | |
53 | SetLastPoint(tpt->GetLastPointX(),tpt->GetLastPointY(),tpt->GetLastPointZ()); | |
54 | SetCharge(tpt->GetCharge()); | |
55 | SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow()); | |
b1886074 | 56 | SetSlice(tpt->GetSlice()); |
d712b5b8 | 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 | { | |
b1886074 | 79 | //Set eta, and calculate fTanl, which is the tan of dipangle |
80 | ||
d712b5b8 | 81 | fEta = f; |
82 | Double_t theta = 2*atan(exp(-1.*fEta)); | |
8b007a2d | 83 | Double_t dipangle = Pi/2 - theta; |
84 | Double_t tgl = tan(dipangle); | |
d712b5b8 | 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 | ||
8b007a2d | 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 | ||
d712b5b8 | 108 | Double_t radius = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]); |
109 | ||
110 | //Get the track parameters | |
b1886074 | 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 * BField ) ; | |
d712b5b8 | 117 | Double_t tPhi0 = GetPsi() + GetCharge() * 0.5 * pi / fabs(GetCharge()) ; |
b1886074 | 118 | Double_t xc = GetCenterX();//x0 - rc * cos(tPhi0) ; |
119 | Double_t yc = GetCenterY();//y0 - rc * sin(tPhi0) ; | |
d712b5b8 | 120 | |
121 | //Check helix and cylinder intersect | |
122 | Double_t fac1 = xc*xc + yc*yc ; | |
123 | Double_t sfac = sqrt( fac1 ) ; | |
b1886074 | 124 | |
d712b5b8 | 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 ; | |
8b007a2d | 142 | |
d712b5b8 | 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 | |
d712b5b8 | 156 | SetR0(radius); |
157 | SetPhi0(phi); | |
158 | SetZ0(z); | |
159 | SetPsi(tPsi); | |
160 | SetFirstPoint(xyz[0],xyz[1],z); | |
8b007a2d | 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); | |
d712b5b8 | 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]); | |
8b007a2d | 168 | //printf("last point %f %f %f\n",xyz[0],xyz[1],xyz[2]); |
d712b5b8 | 169 | } |
170 | ||
171 | void AliL3HoughTrack::SetTrackParameters(Double_t kappa,Double_t phi,Int_t weight) | |
172 | { | |
173 | ||
174 | fWeight = weight; | |
175 | fMinDist = 100000; | |
176 | SetKappa(kappa); | |
177 | SetPhi0(phi); | |
b1886074 | 178 | Double_t pt = fabs(BFACT*BField/kappa); |
d712b5b8 | 179 | SetPt(pt); |
180 | Double_t radius = 1/fabs(kappa); | |
181 | SetRadius(radius); | |
182 | SetFirstPoint(0,0,0); | |
8b007a2d | 183 | SetPsi(phi); //Psi = Phi when first point is vertex |
184 | SetR0(0); | |
d712b5b8 | 185 | Double_t charge = -1.*kappa; |
8b007a2d | 186 | SetCharge((Int_t)copysign(1.,charge)); |
187 | ||
188 | Double_t trackPhi0 = GetPsi() + charge*0.5*Pi/fabs(charge); | |
189 | Double_t xc = GetFirstPointX() - GetRadius() * cos(trackPhi0) ; | |
190 | Double_t yc = GetFirstPointY() - GetRadius() * sin(trackPhi0) ; | |
d712b5b8 | 191 | SetCenterX(xc); |
192 | SetCenterY(yc); | |
193 | fIsHelix = true; | |
8b007a2d | 194 | |
d712b5b8 | 195 | } |
196 | ||
197 | void AliL3HoughTrack::SetLineParameters(Double_t psi,Double_t D,Int_t weight,Int_t *rowrange,Int_t ref_row) | |
198 | { | |
199 | //Initialize a track piece, not yet a track | |
200 | //Used in case of straight line transformation | |
201 | ||
202 | //Transform line parameters to coordinate system of slice: | |
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 | */ |