Bugfix in SetTrackParameters, wrong charge...
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughTrack.cxx
CommitLineData
d712b5b8 1
2#include <math.h>
3
4#include "AliL3Logging.h"
5#include "AliL3Transform.h"
6#include "AliL3Defs.h"
7#include "AliL3HoughTrack.h"
8
9ClassImp(AliL3HoughTrack)
10
11
12AliL3HoughTrack::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
26AliL3HoughTrack::~AliL3HoughTrack()
27{
28 //Destructor
29 if(fTransform)
30 delete fTransform;
31}
32
33void 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
59Int_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
67void AliL3HoughTrack::SetEta(Double_t f)
68{
69
70 fEta = f;
71 Double_t theta = 2*atan(exp(-1.*fEta));
8b007a2d 72 Double_t dipangle = Pi/2 - theta;
73 Double_t tgl = tan(dipangle);
d712b5b8 74 SetTgl(tgl);
75}
76
77void AliL3HoughTrack::CalculateHelix()
78{
79 return;
80}
81
82void 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
8b007a2d 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
d712b5b8 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 ;
8b007a2d 128
d712b5b8 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
d712b5b8 142 SetR0(radius);
143 SetPhi0(phi);
144 SetZ0(z);
145 SetPsi(tPsi);
146 SetFirstPoint(xyz[0],xyz[1],z);
8b007a2d 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);
d712b5b8 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]);
8b007a2d 154 //printf("last point %f %f %f\n",xyz[0],xyz[1],xyz[2]);
d712b5b8 155}
156
157void 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);
8b007a2d 169 SetPsi(phi); //Psi = Phi when first point is vertex
170 SetR0(0);
d712b5b8 171 Double_t charge = -1.*kappa;
8b007a2d 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) ;
d712b5b8 177 SetCenterX(xc);
178 SetCenterY(yc);
179 fIsHelix = true;
8b007a2d 180
d712b5b8 181}
182
183void 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
201void 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
212void 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/*
233Double_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
274Bool_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
311Bool_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*/