]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HLT/TPCLib/AliHLTTPCFitter.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
removing DigitReaderPacked and DigitReaderDecoder, all raw data processing goes via...
[u/mrichter/AliRoot.git] / HLT / TPCLib / AliHLTTPCFitter.cxx
CommitLineData
a6c02c85 1// @(#) $Id$
4aa41877 2// Original: AliHLTFitter.cxx,v 1.14 2005/06/14 10:55:21 cvetan
a6c02c85 3
6235cd38 4/**************************************************************************
9be2600f 5 * This file is property of and copyright by the ALICE HLT Project *
6 * ALICE Experiment at CERN, All rights reserved. *
6235cd38 7 * *
9be2600f 8 * Primary Authors: Anders Vestbo, maintained by *
9 * Matthias Richter <Matthias.Richter@ift.uib.no> *
10 * for The ALICE HLT Project. *
6235cd38 11 * *
12 * Permission to use, copy, modify and distribute this software and its *
13 * documentation strictly for non-commercial purposes is hereby granted *
14 * without fee, provided that the above copyright notice appears in all *
15 * copies and that both the copyright notice and this permission notice *
16 * appear in the supporting documentation. The authors make no claims *
17 * about the suitability of this software for any purpose. It is *
18 * provided "as is" without express or implied warranty. *
19 **************************************************************************/
20
21/** @file AliHLTTPCFitter.cxx
22 @author Anders Vestbo, maintained by Matthias Richter
23 @date
24 @brief Fit class HLT for helix
a6c02c85 25*/
26
27#include <math.h>
a6c02c85 28#include "AliHLTTPCLogging.h"
29#include "AliHLTTPCFitter.h"
30#include "AliHLTTPCVertex.h"
31#include "AliHLTTPCTrack.h"
32#include "AliHLTTPCSpacePointData.h"
33#include "AliHLTTPCMemHandler.h"
34#include "AliHLTTPCTransform.h"
a7c32962 35//#include "AliHLTTPC.h"
a6c02c85 36
37#if __GNUC__ >= 3
38using namespace std;
39#endif
40
41ClassImp(AliHLTTPCFitter)
42
43
44AliHLTTPCFitter::AliHLTTPCFitter()
6235cd38 45 :
46 fTrack(NULL),
47 fVertex(NULL),
48 fVertexConstraint(0)
a6c02c85 49{
50 //constructor
e419b223 51 // see header file for class documentation
52 // or
53 // refer to README to build package
54 // or
55 // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
a6c02c85 56 memset(fClusters,0,36*6*sizeof(AliHLTTPCSpacePointData*));
6235cd38 57 memset(fNcl,0,36*6*sizeof(UInt_t));
58}
59
a6c02c85 60AliHLTTPCFitter::AliHLTTPCFitter(AliHLTTPCVertex *vertex,Bool_t vertexconstraint)
e419b223 61 :
62 fTrack(NULL),
63 fVertex(vertex),
64 fVertexConstraint(vertexconstraint)
a6c02c85 65{
66 //constructor
67 fTrack=0;
a6c02c85 68 memset(fClusters,0,36*6*sizeof(AliHLTTPCSpacePointData*));
69}
70
71AliHLTTPCFitter::~AliHLTTPCFitter()
72{
73 //destructor
74 for(Int_t i=0; i<36; i++)
75 {
76 for(Int_t j=0; j<6; j++)
77 {
78 if(fClusters[i][j])
79 delete [] fClusters[i][j];
80 }
81 }
82}
83
84void AliHLTTPCFitter::LoadClusters(Char_t *path,Int_t event,Bool_t sp)
85{
86 //load clusters
379ff4f0 87 const Int_t fnamelen=256;
88 Char_t fname[fnamelen];
a6c02c85 89 AliHLTTPCMemHandler *clusterfile[36][6];
90 for(Int_t s=0; s<=35; s++)
91 {
92 for(Int_t p=0; p<6; p++)
93 {
94 Int_t patch;
95 if(sp==kTRUE)
96 patch=-1;
97 else
98 patch=p;
99 if(fClusters[s][p])
100 delete fClusters[s][p];
101 fClusters[s][p] = 0;
102 clusterfile[s][p] = new AliHLTTPCMemHandler();
379ff4f0 103 snprintf(fname,fnamelen,"%s/points_%d_%d_%d.raw",path,event,s,patch);
a6c02c85 104 if(!clusterfile[s][p]->SetBinaryInput(fname))
105 {
106 delete clusterfile[s][p];
107 clusterfile[s][p] = 0;
108 continue;
109 }
110 fClusters[s][p] = (AliHLTTPCSpacePointData*)clusterfile[s][p]->Allocate();
111 clusterfile[s][p]->Binary2Memory(fNcl[s][p],fClusters[s][p]);
112 clusterfile[s][p]->CloseBinaryInput();
113 if(sp==kTRUE)
114 break;
115 }
116 }
117}
118
119void AliHLTTPCFitter::SortTrackClusters(AliHLTTPCTrack *track) const
120{
121 //Sort the internal cluster list in each track with respect to row numbering.
122 //This may be necessary when no conventional track follower has been
123 //applied, in which the cluster list has been maintained in a more
124 //arbitrary fashion.
125
126 Int_t nhits = track->GetNHits();
127 Int_t *ids = (Int_t*)track->GetHitNumbers();
128 Int_t *origids = new Int_t[nhits];
129 Int_t *mk = new Int_t[nhits];
130 Int_t k;
131
132 for(k=0; k<nhits; k++) {origids[k] = ids[k]; mk[k] = -1;}
133
134 Int_t slice,patch,id,padrow,maxrow,maxk;
135 UInt_t pos;
136 for(Int_t j=0; j<nhits; j++)
137 {
138 maxrow=-1;
139 maxk=200;
140 for(k=0; k<nhits; k++)
141 {
142 id=ids[k];
143 if(id < 0) continue;
a371a266 144 slice = AliHLTTPCSpacePointData::GetSlice(id);
145 patch = AliHLTTPCSpacePointData::GetPatch(id);
146 pos = AliHLTTPCSpacePointData::GetNumber(id);
a6c02c85 147 AliHLTTPCSpacePointData *points = fClusters[slice][patch];
148 padrow = points[pos].fPadRow;
149 if(padrow > maxrow)
150 {
151 maxrow = padrow;
152 maxk=k;
153 }
154 }
155 mk[j]=maxk;
156 ids[maxk]=-1;
157 }
158
159 for(k=0; k<nhits; k++)
160 ids[k] = origids[mk[k]];
161 delete [] origids;
162 delete [] mk;
163}
164
165Int_t AliHLTTPCFitter::FitHelix(AliHLTTPCTrack *track)
166{
167 //fit helix parameters
168 fTrack = track;
169 if(FitCircle())
170 {
171 LOG(AliHLTTPCLog::kError,"AliHLTTPCFitter::FitHelix","TrackFit")<<AliHLTTPCLog::kDec<<
172 "Problems during circle fit"<<ENDLOG;
173 return 1;
174 }
175 if(FitLine())
176 {
177 LOG(AliHLTTPCLog::kError,"AliHLTTPCFitter::FitHelix","TrackFit")<<AliHLTTPCLog::kDec<<
178 "Problems during line fit"<<ENDLOG;
179 return 1;
180 }
181 return 0;
182}
183
184Int_t AliHLTTPCFitter::FitCircle()
185{
186 //-----------------------------------------------------------------
187 //Fits circle parameters using algorithm
188 //described by ChErnov and Oskov in Computer Physics
189 //Communications.
190 //
191 //Written in FORTRAN by Jawluen Tang, Physics department , UT-Austin
192 //Moved to C by Pablo Yepes
193 //Moved to AliROOT by ASV.
194 //------------------------------------------------------------------
379ff4f0 195
196 if (!fTrack) return -1;
a6c02c85 197
198 Double_t wsum = 0.0 ;
199 Double_t xav = 0.0 ;
200 Double_t yav = 0.0 ;
201
202 //
203 // Loop over hits calculating average
204 Double_t * fXYWeight = new Double_t[(fTrack->GetNHits())];
205 UInt_t *hitnum = fTrack->GetHitNumbers();
379ff4f0 206 if (!fXYWeight || !hitnum) {
207 if (fXYWeight) delete [] fXYWeight;
208 return -1;
209 }
210 memset(fXYWeight, 0, fTrack->GetNHits()*sizeof(Double_t));
a6c02c85 211 for(Int_t i=0; i<fTrack->GetNHits(); i++)
212 {
213 UInt_t id = hitnum[i];
a371a266 214 Int_t slice = AliHLTTPCSpacePointData::GetSlice(id);
215 Int_t patch = AliHLTTPCSpacePointData::GetPatch(id);
216 UInt_t pos = AliHLTTPCSpacePointData::GetNumber(id);
a6c02c85 217 AliHLTTPCSpacePointData *points = fClusters[slice][patch];
218 fXYWeight[i] = 1./ (Double_t)(points[pos].fSigmaY2 + points[pos].fSigmaY2);
219 wsum += fXYWeight[i];
220 xav += fXYWeight[i]*points[pos].fX;
221 yav += fXYWeight[i]*points[pos].fY;
222 }
223 if (fVertexConstraint == kTRUE)
224 {
225 wsum += fVertex->GetXYWeight() ;
226 xav += fVertex->GetX() ;
227 yav += fVertex->GetY() ;
228 }
229
230 xav = xav / wsum ;
231 yav = yav / wsum ;
232//
233// CALCULATE <X**2>, <XY>, AND <Y**2> WITH <X> = 0, & <Y> = 0
234//
235 Double_t xxav = 0.0 ;
236 Double_t xyav = 0.0 ;
237 Double_t yyav = 0.0 ;
238 Double_t xi, yi ;
239
240 for(Int_t i=0; i<fTrack->GetNHits(); i++)
241 {
242 UInt_t id = hitnum[i];
a371a266 243 Int_t slice = AliHLTTPCSpacePointData::GetSlice(id);
244 Int_t patch = AliHLTTPCSpacePointData::GetPatch(id);
245 UInt_t pos = AliHLTTPCSpacePointData::GetNumber(id);
a6c02c85 246 AliHLTTPCSpacePointData *points = fClusters[slice][patch];
247
248 xi = points[pos].fX -xav;
249 yi = points[pos].fY - yav ;
250 xxav += xi * xi * fXYWeight[i];
251 xyav += xi * yi * fXYWeight[i];
252 yyav += yi * yi * fXYWeight[i];
253 }
254
255 if (fVertexConstraint == kTRUE)
256 {
257 xi = fVertex->GetX() - xav ;
258 yi = fVertex->GetY() - yav ;
259 xxav += xi * xi * fVertex->GetXYWeight() ;
260 xyav += xi * yi * fVertex->GetXYWeight() ;
261 yyav += yi * yi * fVertex->GetXYWeight() ;
262 }
263 xxav = xxav / wsum ;
264 xyav = xyav / wsum ;
265 yyav = yyav / wsum ;
266//
267//--> ROTATE COORDINATES SO THAT <XY> = 0
268//
269//--> SIGN(C**2 - S**2) = SIGN(XXAV - YYAV) >
270//--> & > ==> NEW : (XXAV-YYAV) > 0
271//--> SIGN(S) = SIGN(XYAV) >
272
273 Double_t a = fabs( xxav - yyav ) ;
274 Double_t b = 4.0 * xyav * xyav ;
275
276 Double_t asqpb = a * a + b ;
277 Double_t rasqpb = sqrt ( asqpb) ;
278
279 Double_t splus = 1.0 + a / rasqpb ;
280 Double_t sminus = b / (asqpb * splus) ;
281
282 splus = sqrt (0.5 * splus ) ;
283 sminus = sqrt (0.5 * sminus) ;
284//
285//-> FIRST REQUIRE : SIGN(C**2 - S**2) = SIGN(XXAV - YYAV)
286//
287 Double_t sinrot, cosrot ;
288 if ( xxav <= yyav ) {
289 cosrot = sminus ;
290 sinrot = splus ;
291 }
292 else {
293 cosrot = splus ;
294 sinrot = sminus ;
295 }
296//
297//-> REQUIRE : SIGN(S) = SIGN(XYAV) * SIGN(C) (ASSUMING SIGN(C) > 0)
298//
299 if ( xyav < 0.0 ) sinrot = - sinrot ;
300//
301//--> WE NOW HAVE THE SMALLEST ANGLE THAT GUARANTEES <X**2> > <Y**2>
302//--> TO GET THE SIGN OF THE CHARGE RIGHT, THE NEW X-AXIS MUST POINT
303//--> OUTWARD FROM THE ORGIN. WE ARE FREE TO CHANGE SIGNS OF BOTH
304//--> COSROT AND SINROT SIMULTANEOUSLY TO ACCOMPLISH THIS.
305//
306//--> CHOOSE SIGN OF C WISELY TO BE ABLE TO GET THE SIGN OF THE CHARGE
307//
308 if ( cosrot*xav+sinrot*yav < 0.0 ) {
309 cosrot = -cosrot ;
310 sinrot = -sinrot ;
311 }
312//
313//-> NOW GET <R**2> AND RSCALE= SQRT(<R**2>)
314//
315 Double_t rrav = xxav + yyav ;
316 Double_t rscale = sqrt(rrav) ;
317
318 xxav = 0.0 ;
319 yyav = 0.0 ;
320 xyav = 0.0 ;
321 Double_t xrrav = 0.0 ;
322 Double_t yrrav = 0.0 ;
323 Double_t rrrrav = 0.0 ;
324
325 Double_t xixi, yiyi, riri, wiriri, xold, yold ;
326
327 for(Int_t i=0; i<fTrack->GetNHits(); i++)
328 {
329 UInt_t id = hitnum[i];
a371a266 330 Int_t slice = AliHLTTPCSpacePointData::GetSlice(id);
331 Int_t patch = AliHLTTPCSpacePointData::GetPatch(id);
332 UInt_t pos = AliHLTTPCSpacePointData::GetNumber(id);
a6c02c85 333 AliHLTTPCSpacePointData *points = fClusters[slice][patch];
334
335 xold = points[pos].fX - xav ;
336 yold = points[pos].fY - yav ;
337 //
338 //--> ROTATE SO THAT <XY> = 0 & DIVIDE BY RSCALE SO THAT <R**2> = 1
339 //
340 xi = ( cosrot * xold + sinrot * yold ) / rscale ;
341 yi = ( -sinrot * xold + cosrot * yold ) / rscale ;
342
343 xixi = xi * xi ;
344 yiyi = yi * yi ;
345 riri = xixi + yiyi ;
346 wiriri = fXYWeight[i] * riri ;
347
348 xyav += fXYWeight[i] * xi * yi ;
349 xxav += fXYWeight[i] * xixi ;
350 yyav += fXYWeight[i] * yiyi ;
351
352 xrrav += wiriri * xi ;
353 yrrav += wiriri * yi ;
354 rrrrav += wiriri * riri ;
355 }
356//
357// Include vertex if required
358//
359 if (fVertexConstraint == kTRUE)
360 {
361 xold = fVertex->GetX() - xav ;
362 yold = fVertex->GetY() - yav ;
363 //
364 //--> ROTATE SO THAT <XY> = 0 & DIVIDE BY RSCALE SO THAT <R**2> = 1
365 //
366 xi = ( cosrot * xold + sinrot * yold ) / rscale ;
367 yi = ( -sinrot * xold + cosrot * yold ) / rscale ;
368
369 xixi = xi * xi ;
370 yiyi = yi * yi ;
371 riri = xixi + yiyi ;
372 wiriri = fVertex->GetXYWeight() * riri ;
373
374 xyav += fVertex->GetXYWeight() * xi * yi ;
375 xxav += fVertex->GetXYWeight() * xixi ;
376 yyav += fVertex->GetXYWeight() * yiyi ;
377
378 xrrav += wiriri * xi ;
379 yrrav += wiriri * yi ;
380 rrrrav += wiriri * riri ;
381 }
382 //
383 //
384 //
385 //--> DIVIDE BY WSUM TO MAKE AVERAGES
386 //
387 xxav = xxav / wsum ;
388 yyav = yyav / wsum ;
389 xrrav = xrrav / wsum ;
390 yrrav = yrrav / wsum ;
391 rrrrav = rrrrav / wsum ;
392 xyav = xyav / wsum ;
393
394 Int_t const kntry = 5 ;
395//
396//--> USE THESE TO GET THE COEFFICIENTS OF THE 4-TH ORDER POLYNIMIAL
397//--> DON'T PANIC - THE THIRD ORDER TERM IS ZERO !
398//
399 Double_t xrrxrr = xrrav * xrrav ;
400 Double_t yrryrr = yrrav * yrrav ;
401 Double_t rrrrm1 = rrrrav - 1.0 ;
402 Double_t xxyy = xxav * yyav ;
403
404 Double_t c0 = rrrrm1*xxyy - xrrxrr*yyav - yrryrr*xxav ;
405 Double_t c1 = - rrrrm1 + xrrxrr + yrryrr - 4.0*xxyy ;
406 Double_t c2 = 4.0 + rrrrm1 - 4.0*xxyy ;
407 Double_t c4 = - 4.0 ;
408//
409//--> COEFFICIENTS OF THE DERIVATIVE - USED IN NEWTON-RAPHSON ITERATIONS
410//
411 Double_t c2d = 2.0 * c2 ;
412 Double_t c4d = 4.0 * c4 ;
413//
414//--> 0'TH VALUE OF LAMDA - LINEAR INTERPOLATION BETWEEN P(0) & P(YYAV)
415//
416// LAMDA = YYAV * C0 / (C0 + YRRSQ * (XXAV-YYAV))
417 Double_t lamda = 0.0 ;
418 Double_t dlamda = 0.0 ;
419//
420 Double_t chiscl = wsum * rscale * rscale ;
421 Double_t dlamax = 0.001 / chiscl ;
422
423 Double_t p, pd ;
424 for ( int itry = 1 ; itry <= kntry ; itry++ ) {
425 p = c0 + lamda * (c1 + lamda * (c2 + lamda * lamda * c4 )) ;
426 pd = (c1 + lamda * (c2d + lamda * lamda * c4d)) ;
427 dlamda = -p / pd ;
428 lamda = lamda + dlamda ;
429 if (fabs(dlamda)< dlamax) break ;
430 }
431
432 //Double_t chi2 = (Double_t)(chiscl * lamda) ;
433 //fTrack->SetChiSq1(chi2);
434 // Double_t dchisq = chiscl * dlamda ;
435 //
436 //--> NOW CALCULATE THE MATRIX ELEMENTS FOR ALPHA, BETA & KAPPA
437 //
438 Double_t h11 = xxav - lamda ;
439 Double_t h14 = xrrav ;
440 Double_t h22 = yyav - lamda ;
441 Double_t h24 = yrrav ;
442 Double_t h34 = 1.0 + 2.0*lamda ;
443 if ( h11 == 0.0 || h22 == 0.0 ){
444 LOG(AliHLTTPCLog::kError,"AliHLTTPCFitter::FitCircle","TrackFit")<<AliHLTTPCLog::kDec<<
445 "Problems fitting circle"<<ENDLOG;
379ff4f0 446 delete [] fXYWeight;
a6c02c85 447 return 1 ;
448 }
449 Double_t rootsq = (h14*h14)/(h11*h11) + 4.0*h34 ;
450
451 Double_t ratio, kappa, beta ;
452 if ( fabs(h22) > fabs(h24) ) {
453 ratio = h24 / h22 ;
454 rootsq = ratio * ratio + rootsq ;
455 kappa = 1.0 / sqrt(rootsq) ;
456 beta = - ratio * kappa ;
457 }
458 else {
459 ratio = h22 / h24 ;
460 rootsq = 1.0 + ratio * ratio * rootsq ;
461 beta = 1.0 / sqrt(rootsq) ;
462 if ( h24 > 0 ) beta = - beta ;
463 kappa = -ratio * beta ;
464 }
465 Double_t alpha = - (h14/h11) * kappa ;
466//
467//--> transform these into the lab coordinate system
468//--> first get kappa and back to real dimensions
469//
470 Double_t kappa1 = kappa / rscale ;
471 Double_t dbro = 0.5 / kappa1 ;
472//
473//--> next rotate alpha and beta and scale
474//
475 Double_t alphar = (cosrot * alpha - sinrot * beta)* dbro ;
476 Double_t betar = (sinrot * alpha + cosrot * beta)* dbro ;
477//
478//--> then translate by (xav,yav)
479//
480 Double_t acent = (double)(xav - alphar) ;
481 Double_t bcent = (double)(yav - betar ) ;
482 Double_t radius = (double)dbro ;
483//
484// Get charge
485//
486 Int_t q = ( ( yrrav < 0 ) ? 1 : -1 ) ;
487 fTrack->SetCharge(q);
488
489 //Set the first point on the track to the space point coordinates of the innermost track
490 //This will be updated to lie on the fit later on (AliHLTTPCTrack::UpdateToFirstPoint).
491 Double_t x0,y0,psi,pt ;
492 Int_t lastid=fTrack->GetNHits()-1;
493 UInt_t id = hitnum[lastid];
a371a266 494 Int_t slice = AliHLTTPCSpacePointData::GetSlice(id);
495 Int_t patch = AliHLTTPCSpacePointData::GetPatch(id);
496 UInt_t pos = AliHLTTPCSpacePointData::GetNumber(id);
a6c02c85 497 AliHLTTPCSpacePointData *points = fClusters[slice][patch];
498 x0 = points[pos].fX;
499 y0 = points[pos].fY;
500 fTrack->SetFirstPoint(x0,y0,0); //Z-value is set in FitLine
501
502 //Set the remaining fit parameters
503 psi = (Double_t)atan2(bcent-y0,acent-x0) ;
504 psi = psi + q * 0.5F * AliHLTTPCTransform::Pi() ;
505 if ( psi < 0 ) psi = psi + 2*AliHLTTPCTransform::Pi();
506
507 pt = (Double_t)(AliHLTTPCTransform::GetBFact() * AliHLTTPCTransform::GetBField() * radius ) ;
508 fTrack->SetPsi(psi);
509 fTrack->SetPt(pt);
510 fTrack->SetRadius(radius);
511 fTrack->SetCenterX(acent);
512 fTrack->SetCenterY(bcent);
513//
514// Get errors from fast fit
515//
516 //if ( getPara()->getErrors ) getErrorsCircleFit ( acent, bcent, radius ) ;
517//
518 delete [] fXYWeight;
519 return 0 ;
520}
521
522//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
523// Fit Line in s-z plane
524//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
525Int_t AliHLTTPCFitter::FitLine ( )
526{
379ff4f0 527 // Fit Line in s-z plane
528 if (!fTrack) return -1;
529
a6c02c85 530 //
531 //Initialization
532 //
533 Double_t sum = 0.F ;
534 Double_t ss = 0.F ;
535 Double_t sz = 0.F ;
536 Double_t sss = 0.F ;
537 Double_t ssz = 0.F ;
538 //
539 //find sum , sums ,sumz, sumss
540 //
541 Double_t dx, dy ;
542 Double_t radius = (Double_t)(fTrack->GetPt() / ( AliHLTTPCTransform::GetBFact() * AliHLTTPCTransform::GetBField() ) ) ;
543
544 Double_t * fS = new Double_t[(fTrack->GetNHits())];
545 Double_t *fZWeight = new Double_t[fTrack->GetNHits()];
546 UInt_t *hitnum = fTrack->GetHitNumbers();
379ff4f0 547 if (!fS || !fZWeight || !hitnum) {
548 if (fS) delete [] fS;
549 if (fZWeight) delete [] fZWeight;
550 return -1;
551 }
552 memset(fS, 0, fTrack->GetNHits()*sizeof(Double_t));
553 memset(fZWeight, 0, fTrack->GetNHits()*sizeof(Double_t));
a6c02c85 554 if (0)//fVertexConstraint==kTRUE)
555 {
556 UInt_t id = hitnum[0];
a371a266 557 Int_t slice = AliHLTTPCSpacePointData::GetSlice(id);
558 Int_t patch = AliHLTTPCSpacePointData::GetPatch(id);
559 UInt_t pos = AliHLTTPCSpacePointData::GetNumber(id);
a6c02c85 560 AliHLTTPCSpacePointData *points = fClusters[slice][patch];
561
562 dx = points[pos].fX - fVertex->GetX();
563 dy = points[pos].fY - fVertex->GetY();
564 }
565 else
566 {
567 UInt_t id = hitnum[0];
a371a266 568 Int_t slice = AliHLTTPCSpacePointData::GetSlice(id);
569 Int_t patch = AliHLTTPCSpacePointData::GetPatch(id);
570 UInt_t posf = AliHLTTPCSpacePointData::GetNumber(id);
a6c02c85 571 AliHLTTPCSpacePointData *pointsf = fClusters[slice][patch];
572 id = hitnum[(fTrack->GetNHits()-1)];
a371a266 573 slice = AliHLTTPCSpacePointData::GetSlice(id);
574 patch = AliHLTTPCSpacePointData::GetPatch(id);
575 UInt_t posl = AliHLTTPCSpacePointData::GetNumber(id);
a6c02c85 576 AliHLTTPCSpacePointData *pointsl = fClusters[slice][patch];
577 dx = pointsf[posf].fX - pointsl[posl].fX;
578 dy = pointsf[posf].fY - pointsl[posl].fY;
579 }
580
581 Double_t localPsi = 0.5F * sqrt ( dx*dx + dy*dy ) / radius ;
582 Double_t totals ;
583
584 if ( fabs(localPsi) < 1. )
585 {
586 totals = 2.0 * radius * asin ( localPsi ) ;
587 }
588 else
589 {
590 totals = 2.0 * radius * AliHLTTPCTransform::Pi() ;
591 }
592
593 Double_t dpsi,s;
594
595 for(Int_t i=0; i<fTrack->GetNHits(); i++)
596 {
597 UInt_t id = hitnum[i];
a371a266 598 Int_t slice = AliHLTTPCSpacePointData::GetSlice(id);
599 Int_t patch = AliHLTTPCSpacePointData::GetPatch(id);
600 UInt_t pos = AliHLTTPCSpacePointData::GetNumber(id);
a6c02c85 601 AliHLTTPCSpacePointData *points = fClusters[slice][patch];
602
603 fZWeight[i] = 1./(Double_t)(points[pos].fSigmaZ2);
604 if(i>0)
605 {
606 id = hitnum[i-1];
a371a266 607 slice = AliHLTTPCSpacePointData::GetSlice(id);
608 patch = AliHLTTPCSpacePointData::GetPatch(id);
609 UInt_t lastpos = AliHLTTPCSpacePointData::GetNumber(id);
a6c02c85 610 AliHLTTPCSpacePointData *lastpoints = fClusters[slice][patch];
611 dx = points[pos].fX -lastpoints[lastpos].fX;
612 dy = points[pos].fY -lastpoints[lastpos].fY;
613 dpsi = 0.5 * (Double_t)sqrt ( dx*dx + dy*dy ) / radius ;
379ff4f0 614 if(fabs(dpsi) > 1) {
615 delete [] fS;
616 delete [] fZWeight;
a6c02c85 617 return 1;
379ff4f0 618 }
a6c02c85 619 fTrack->SetPsierr(dpsi);
620 s = fS[i-1] - 2.0 * radius * (Double_t)asin ( dpsi ) ;
621 fS[i]=s;
622 }
623 else
624 fS[i]=totals;
625
626 sum += fZWeight[i];
627 ss += fZWeight[i] * fS[i];
628 sz += fZWeight[i] * points[pos].fZ;
629 sss += fZWeight[i] * fS[i] * fS[i];
630 ssz += fZWeight[i] * fS[i] * points[pos].fZ;
631
632 }
633
634
635 Double_t chi2,det = sum * sss - ss * ss;
636 if ( fabs(det) < 1e-20)
637 {
638 chi2 = 99999.F ;
639 //fTrack->SetChiSq2(chi2);
379ff4f0 640 delete [] fS;
641 delete [] fZWeight;
a6c02c85 642 return 0 ;
643 }
644
645 //Compute the best fitted parameters A,B
646 Double_t tanl,z0,dtanl,dz0;
647
648 tanl = (Double_t)((sum * ssz - ss * sz ) / det );
649 z0 = (Double_t)((sz * sss - ssz * ss ) / det );
650
651 fTrack->SetTgl(tanl);
652 fTrack->SetZ0(z0);
653
654 //calculate chi-square
655 chi2 = 0.;
656 Double_t r1 ;
657
658 for(Int_t i=0; i<fTrack->GetNHits(); i++)
659 {
660 UInt_t id = hitnum[i];
a371a266 661 Int_t slice = AliHLTTPCSpacePointData::GetSlice(id);
662 Int_t patch = AliHLTTPCSpacePointData::GetPatch(id);
663 UInt_t pos = AliHLTTPCSpacePointData::GetNumber(id);
a6c02c85 664 AliHLTTPCSpacePointData *points = fClusters[slice][patch];
665 r1 = points[pos].fZ - tanl * fS[i] - z0 ;
666 chi2 += (Double_t) ( (Double_t)(fZWeight[i]) * (r1 * r1) );
667 }
668
669 //fTrack->SetChiSq2(chi2);
670 //
671 //calculate estimated variance
672 //varsq=chi/(double(n)-2.)
673 //calculate covariance matrix
674 //siga=sqrt(varsq*sxx/det)
675 //sigb=sqrt(varsq*sum/det)
676 //
677 dtanl = (Double_t) ( sum / det );
678 dz0 = (Double_t) ( sss / det );
679
680 fTrack->SetTglerr(dtanl);
681 fTrack->SetZ0err(dz0);
682 delete [] fZWeight;
683 delete [] fS;
684 return 0 ;
685}