]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSTPArrayFit.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
correcting printf format specifiers
[u/mrichter/AliRoot.git] / ITS / AliITSTPArrayFit.cxx
CommitLineData
6be22b3f 1/**************************************************************************
2 * Copyright(c) 2009-2011, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15/* $Id$ */
16///////////////////////////////////////////////////////////////////////////////////////////////
17// //
18// The line is defined by equations (1) //
19// a0*z+a1*x-a0*a1=0 and //
20// b0*z+b1*y-b0*b1=0 //
21// where x,y,z are NOT the lab axes but z is the lab axis along which the track //
22// has the largest lever arm and x,y are the remaining 2 axis in //
23// the order of fgkAxisID[z][0], fgkAxisID[z][1] //
24// The parameters are fParams[kA0,kB0,kA1,kB1] and the axis chosen as the independent //
25// var. is fParAxis (i.e. if fParAxis==kZ, then a0=ax,b0=bx, a1=ay,b1=by) //
26// //
27// //
28// The helix is defined by the equations (2) //
29// X(t) = (dr+R)*cos(phi0) - (R+sum{dRi})*cos(t+phi0) + sum{dRi*cos(phi0+ti)} //
30// Y(t) = (dr+R)*sin(phi0) - (R+sum{dRi})*sin(t+phi0) + sum{dRi*sin(phi0+ti)} //
31// Z(t) = dz - (R+sum{dRi})*t*tg(dip) + sum{dRi*ti}*tg(dip) //
32// where dRi is the change of the radius due to the ELoss at parameter ti //
33// //
34// Author: ruben.shahoyan@cern.ch //
35// //
36///////////////////////////////////////////////////////////////////////////////////////////////
37
38#include "AliITSTPArrayFit.h"
39#include "AliExternalTrackParam.h"
40#include "AliSymMatrix.h"
41#include "AliLog.h"
42#include "AliParamSolver.h"
43#include "AliGeomManager.h"
44#include "AliITSgeomTGeo.h"
45#include "AliTracker.h"
46#include <TRandom.h>
ef24eb3b 47#include <TArrayD.h>
6be22b3f 48
49ClassImp(AliITSTPArrayFit)
50
51const Int_t AliITSTPArrayFit::fgkAxisID[3][3] = {
52 {AliITSTPArrayFit::kY,AliITSTPArrayFit::kZ,AliITSTPArrayFit::kX},
53 {AliITSTPArrayFit::kZ,AliITSTPArrayFit::kX,AliITSTPArrayFit::kY},
54 {AliITSTPArrayFit::kX,AliITSTPArrayFit::kY,AliITSTPArrayFit::kZ} };
55
56const Int_t AliITSTPArrayFit::fgkAxisCID[3][6] = {
57 {AliITSTPArrayFit::kYY,AliITSTPArrayFit::kYZ,AliITSTPArrayFit::kXY,
58 AliITSTPArrayFit::kZZ,AliITSTPArrayFit::kXZ,AliITSTPArrayFit::kXX},
59 //
60 {AliITSTPArrayFit::kZZ,AliITSTPArrayFit::kXZ,AliITSTPArrayFit::kYZ,
61 AliITSTPArrayFit::kXX,AliITSTPArrayFit::kYX,AliITSTPArrayFit::kYY},
62 //
63 {AliITSTPArrayFit::kXX,AliITSTPArrayFit::kXY,AliITSTPArrayFit::kXZ,
64 AliITSTPArrayFit::kYY,AliITSTPArrayFit::kYZ,AliITSTPArrayFit::kZZ}
65};
66//
67
68const Double_t AliITSTPArrayFit::fgkAlmostZero = 1E-55;
69const Double_t AliITSTPArrayFit::fgkCQConv = 0.299792458e-3;// R = PT/Bz/fgkCQConv with GeV,kGauss,cm
70const Double_t AliITSTPArrayFit::fgkZSpanITS[AliITSTPArrayFit::kMaxLrITS] = {
71 36. ,14.1,14.1, 38., 22.2,29.7, 51. ,43.1,48.9};
72
73const Double_t AliITSTPArrayFit::fgkRLayITS[AliITSTPArrayFit::kMaxLrITS] = {
74 2.94, 3.9,7.6, 11.04, 15.0,23.9, 29.44 ,38.0,43.0};
75
76const Int_t AliITSTPArrayFit::fgkPassivLrITS[3] =
77 {AliITSTPArrayFit::kLrBeamPime,AliITSTPArrayFit::kLrShield1,AliITSTPArrayFit::kLrShield2};
78
79const Int_t AliITSTPArrayFit::fgkActiveLrITS[6] =
80 {AliITSTPArrayFit::kLrSPD1,AliITSTPArrayFit::kLrSPD2,
81 AliITSTPArrayFit::kLrSDD1,AliITSTPArrayFit::kLrSDD2,
82 AliITSTPArrayFit::kLrSSD1,AliITSTPArrayFit::kLrSSD2};
83
84Double_t AliITSTPArrayFit::fgRhoLITS[AliITSTPArrayFit::kMaxLrITS] = {
85 1.48e-01, 2.48e-01,2.57e-01, 1.34e-01, 3.34e-01,3.50e-01, 2.22e-01, 2.38e-01,2.25e-01};
86
87//____________________________________________________
88AliITSTPArrayFit::AliITSTPArrayFit() :
24391cd5 89 fkPoints(0),fParSol(0),fBz(0),fCharge(0),fPntFirst(-1),
6be22b3f 90 fPntLast(-1),fNPBooked(0),fParAxis(-1),fCovI(0),fChi2NDF(0),
ef24eb3b 91 fMaxIter(20),fIter(0),fEps(1e-6),fMass(0),fSwitch2Line(kFALSE),fMaxRforHelix(6.e5),
1d06ac63 92 fkAxID(0),fkAxCID(0),fCurT(0),
6be22b3f 93 fFirstPosT(0),fNElsPnt(0),fElsId(0),fElsDR(0)
94{
95 // default constructor
96 for (int i=kMaxParam;i--;) fParams[i] = 0;
97 for (int i=kMaxParamSq;i--;) fParamsCov[i] = 0;
98 SetMass();
99}
100
101//____________________________________________________
102AliITSTPArrayFit::AliITSTPArrayFit(Int_t np) :
24391cd5 103 fkPoints(0),fParSol(0),fBz(0),fCharge(0),fPntFirst(-1),
6be22b3f 104 fPntLast(-1),fNPBooked(np),fParAxis(-1),fCovI(0),fChi2NDF(0),
ef24eb3b 105 fMaxIter(20),fIter(0),fEps(1e-6),fMass(0),fSwitch2Line(kFALSE),fMaxRforHelix(6.e5),
1d06ac63 106 fkAxID(0),fkAxCID(0),fCurT(0),fFirstPosT(0),fNElsPnt(0),fElsId(0),fElsDR(0)
6be22b3f 107{
108 // constructor with booking of np points
109 for (int i=kMaxParam;i--;) fParams[i] = 0;
110 for (int i=kMaxParamSq;i--;) fParamsCov[i] = 0;
111 InitAux();
112 SetEps();
113 SetMass();
114 SetMaxIterations();
115}
116
117//____________________________________________________
118AliITSTPArrayFit::AliITSTPArrayFit(const AliITSTPArrayFit &src) :
24391cd5 119 TObject(src),fkPoints(src.fkPoints),fParSol(0),fBz(src.fBz),
6be22b3f 120 fCharge(src.fCharge),fPntFirst(src.fPntFirst),fPntLast(src.fPntLast),fNPBooked(src.fNPBooked),
121 fParAxis(src.fParAxis),fCovI(0),fChi2NDF(0),fMaxIter(20),fIter(0),fEps(0),fMass(src.fMass),
1d06ac63 122 fSwitch2Line(src.fSwitch2Line),fMaxRforHelix(src.fMaxRforHelix),fkAxID(0),fkAxCID(0),fCurT(0),
123 fFirstPosT(0),fNElsPnt(0),fElsId(0),fElsDR(0)
6be22b3f 124{
125 // copy constructor
126 InitAux();
127 memcpy(fCovI,src.fCovI,fNPBooked*6*sizeof(Double_t));
128 for (int i=kMaxParam;i--;) fParams[i] = src.fParams[i];
129 for (int i=kMaxParamSq;i--;) fParamsCov[i] = src.fParamsCov[i];
130 memcpy(fCurT,src.fCurT,fNPBooked*sizeof(Double_t));
131 SetEps(src.fEps);
132 SetMaxIterations(src.fMaxIter);
133 //
134}
135
136//____________________________________________________
137AliITSTPArrayFit &AliITSTPArrayFit::operator =(const AliITSTPArrayFit& src)
138{
139 // assignment operator
140 if (this==&src) return *this;
141 ((TObject*)this)->operator=(src);
24391cd5 142 fkPoints = src.fkPoints;
6be22b3f 143 if (!fParSol) fParSol = new AliParamSolver(*src.fParSol);
144 else *fParSol = *src.fParSol;
145 fBz = src.fBz;
146 fCharge = src.fCharge;
147 fNPBooked = src.fNPBooked;
148 fPntFirst = src.fPntFirst;
149 fPntLast = src.fPntLast;
150 InitAux();
151 memcpy(fCovI,src.fCovI,fNPBooked*6*sizeof(Double_t));
152 for (int i=kMaxParam;i--;) fParams[i] = src.fParams[i];
153 for (int i=kMaxParamSq;i--;) fParamsCov[i] = src.fParamsCov[i];
154 SetParAxis(src.fParAxis);
155 fNElsPnt = src.fNElsPnt;
156 fFirstPosT = src.fFirstPosT;
157 memcpy(fCurT ,src.fCurT ,fNPBooked*sizeof(Double_t));
158 memcpy(fElsId ,src.fElsId ,fNPBooked*sizeof(Int_t));
159 memcpy(fElsDR ,src.fElsDR ,fNPBooked*sizeof(Double_t));
160 memcpy(fCurT ,src.fCurT ,fNPBooked*sizeof(Double_t));
161 SetEps(src.fEps);
162 SetMaxIterations(src.fMaxIter);
163 //
164 return *this;
165 //
166}
167
168//____________________________________________________
169AliITSTPArrayFit::~AliITSTPArrayFit()
170{
171 // destructor
172 delete fParSol;
173 delete[] fCovI;
174 delete[] fCurT;
175 delete[] fElsId;
176 delete[] fElsDR;
177}
178
179//____________________________________________________
180void AliITSTPArrayFit::Reset()
181{
182 // reset to process new track
183 if (fParSol) fParSol->Clear();
24391cd5 184 fkPoints=0;
6be22b3f 185 fNElsPnt = 0;
186 fFirstPosT = 0;
187 // fBz = 0;
188 fCharge = 0;
189 fIter = 0;
190 fPntFirst=fPntLast=-1;
191 SetParAxis(-1);
1d06ac63 192 fSwitch2Line = kFALSE;
6be22b3f 193 ResetBit(kFitDoneBit|kCovInvBit);
194}
195
196//____________________________________________________
197void AliITSTPArrayFit::AttachPoints(const AliTrackPointArray* points, Int_t pfirst,Int_t plast)
198{
199 // create from piece of AliTrackPointArray
200 Reset();
24391cd5 201 fkPoints = points;
6be22b3f 202 int np = points->GetNPoints();
203 if (fNPBooked<np) {
204 fNPBooked = np;
205 InitAux();
206 }
207 fPntFirst = pfirst<0 ? 0 : pfirst;
208 fPntLast = plast<fPntFirst ? np-1 : plast;
209 //
210 for (int i=kMaxParam;i--;) fParams[i] = 0;
211 for (int i=kMaxParamSq;i--;) fParamsCov[i] = 0;
212 //
213 InvertPointsCovMat();
214 //
215}
216
217//____________________________________________________
218Bool_t AliITSTPArrayFit::SetFirstLast(Int_t pfirst,Int_t plast)
219{
220 // set first and last point to fit
24391cd5 221 const AliTrackPointArray* pnts = fkPoints;
6be22b3f 222 if (!pnts) {AliError("TrackPointArray is not attached yet"); return kFALSE;}
223 AttachPoints(pnts,pfirst,plast);
224 return kTRUE;
225 //
226}
227
228//____________________________________________________
229Bool_t AliITSTPArrayFit::InvertPointsCovMat()
230{
231 // invert the cov.matrices of the points
232 for (int i=fPntFirst;i<=fPntLast;i++) {
233 //
24391cd5 234 float *cov = (float*)fkPoints->GetCov() + i*6; // pointer on cov.matrix
6be22b3f 235 //
236 Double_t t0 = cov[kYY]*cov[kZZ] - cov[kYZ]*cov[kYZ];
237 Double_t t1 = cov[kXY]*cov[kZZ] - cov[kXZ]*cov[kYZ];
238 Double_t t2 = cov[kXY]*cov[kYZ] - cov[kXZ]*cov[kYY];
239 Double_t det = cov[kXX]*t0 - cov[kXY]*t1 + cov[kXZ]*t2;
66214d86 240 if (IsZero(det,1e-18)) { // one of errors is 0, fix this
241 double norm[3];
242 TGeoHMatrix hcov;
243 TGeoRotation hrot,hrotI;
244 GetNormal(norm,cov);
245 double phi = TMath::ATan2(norm[1],norm[0]);
246 hrot.SetAngles(-phi*TMath::RadToDeg(),0.,0.);
247 hrotI.SetAngles(phi*TMath::RadToDeg(),0.,0.);
248 //
249 Double_t hcovel[9];
250 hcovel[0] = cov[kXX];
251 hcovel[1] = cov[kXY];
252 hcovel[2] = cov[kXZ];
253 hcovel[3] = cov[kXY];
254 hcovel[4] = cov[kYY];
255 hcovel[5] = cov[kYZ];
256 hcovel[6] = cov[kXZ];
257 hcovel[7] = cov[kYZ];
258 hcovel[8] = cov[kZZ];
259 hcov.SetRotation(hcovel);
260 //
261 Double_t *hcovscl = hcov.GetRotationMatrix();
262 // printf(">> %+e %+e %+e\n %+e %+e %+e\n %+e %+e %+e\n\n",hcovscl[0],hcovscl[1],hcovscl[2],hcovscl[3],hcovscl[4],hcovscl[5],hcovscl[6],hcovscl[7],hcovscl[8]);
263 // printf("Rot by %+.e (%+.e %+.e %+.e)\n",phi, norm[0],norm[1],norm[2]);
264 hcov.Multiply(&hrotI);
265 hcov.MultiplyLeft(&hrot);
266 // printf("|| %+e %+e %+e\n %+e %+e %+e\n %+e %+e %+e\n\n",hcovscl[0],hcovscl[1],hcovscl[2],hcovscl[3],hcovscl[4],hcovscl[5],hcovscl[6],hcovscl[7],hcovscl[8]);
267 if (hcovscl[0]<1e-8) hcovscl[0] = 1e-8;
268 if (hcovscl[4]<1e-8) hcovscl[4] = 1e-8;
269 if (hcovscl[8]<1e-8) hcovscl[8] = 1e-8;
270 // printf("** %+e %+e %+e\n %+e %+e %+e\n %+e %+e %+e\n\n",hcovscl[0],hcovscl[1],hcovscl[2],hcovscl[3],hcovscl[4],hcovscl[5],hcovscl[6],hcovscl[7],hcovscl[8]);
271 hcov.Multiply(&hrot);
272 hcov.MultiplyLeft(&hrotI);
273 // printf("^^ %+e %+e %+e\n %+e %+e %+e\n %+e %+e %+e\n\n",hcovscl[0],hcovscl[1],hcovscl[2],hcovscl[3],hcovscl[4],hcovscl[5],hcovscl[6],hcovscl[7],hcovscl[8]);
274 cov[kXX] = hcovscl[0];
275 cov[kXY] = hcovscl[1];
276 cov[kXZ] = hcovscl[2];
277 cov[kYY] = hcovscl[4];
278 cov[kYZ] = hcovscl[5];
279 cov[kZZ] = hcovscl[8];
280 }
281 t0 = cov[kYY]*cov[kZZ] - cov[kYZ]*cov[kYZ];
282 t1 = cov[kXY]*cov[kZZ] - cov[kXZ]*cov[kYZ];
283 t2 = cov[kXY]*cov[kYZ] - cov[kXZ]*cov[kYY];
284 det = cov[kXX]*t0 - cov[kXY]*t1 + cov[kXZ]*t2;
285 //
abd7ef79 286 AliDebug(2,Form("%+.4e %+.4e %+.4e -> %+.4e",t0,t1,t2,det));
66214d86 287 if (IsZero(det,fgkAlmostZero)) {
6be22b3f 288 AliInfo(Form("Cov.Matrix for point %d is singular",i));
289 return kFALSE;
290 }
291 //
292 Double_t *covI = GetCovI(i);
293 covI[kXX] = t0/det;
294 covI[kXY] = -t1/det;
295 covI[kXZ] = t2/det;
296 covI[kYY] = (cov[kXX]*cov[kZZ] - cov[kXZ]*cov[kXZ])/det;
297 covI[kYZ] = (cov[kXY]*cov[kXZ] - cov[kXX]*cov[kYZ])/det;
298 covI[kZZ] = (cov[kXX]*cov[kYY] - cov[kXY]*cov[kXY])/det;
299 //
300 }
301 SetCovInv();
302 return kTRUE;
303}
304
305//____________________________________________________
306void AliITSTPArrayFit::InitAux()
307{
308 // init auxiliary space
309 if (fCovI) delete[] fCovI;
310 if (fCurT) delete[] fCurT;
311 //
312 fCovI = new Double_t[6*fNPBooked];
313 fCurT = new Double_t[fNPBooked+kMaxLrITS];
314 fElsId = new Int_t[fNPBooked+kMaxLrITS];
315 fElsDR = new Double_t[fNPBooked+kMaxLrITS];
316 memset(fElsDR,0,(fNPBooked+kMaxLrITS)*sizeof(Double_t));
317 memset(fCovI,0,fNPBooked*6*sizeof(Double_t));
318 //
319}
320
321//____________________________________________________
322Bool_t AliITSTPArrayFit::FitLineCrude()
323{
324 // perform linear fit w/o accounting the errors
325 // fit is done in the parameterization
326 // x = res[0] + res[1]*z
327 // y = res[2] + res[3]*z
328 // where x,y,z are NOT the lab axes but z is the lab axis along which the track
329 // has the largest lever arm and x,y are the remaining 2 axis in
330 // the order of fgkAxisID[z][0], fgkAxisID[z][1]
331 //
332 int np = fPntLast - fPntFirst + 1;
333 if (np<2) {
334 AliError("At least 2 points are needed for straight line fit");
335 return kFALSE;
336 }
337 //
338 if (fParAxis<0) SetParAxis(ChoseParAxis());
339 Double_t sZ=0,sZZ=0,sY=0,sYZ=0,sX=0,sXZ=0,det=0;
340 //
24391cd5 341 const float *coord[3] = {fkPoints->GetX(),fkPoints->GetY(),fkPoints->GetZ()};
6be22b3f 342 const Float_t *varZ = coord[ fParAxis ];
343 const Float_t *varX = coord[ fkAxID[kX] ];
344 const Float_t *varY = coord[ fkAxID[kY] ];
345 //
346 for (int i=fPntFirst;i<=fPntLast;i++) {
347 sZ += varZ[i];
348 sZZ += varZ[i]*varZ[i];
349 //
350 sX += varX[i];
351 sXZ += varX[i]*varZ[i];
352 //
353 sY += varY[i];
354 sYZ += varY[i]*varZ[i];
355 }
356 det = sZZ*np-sZ*sZ;
357 if (TMath::Abs(det)<fgkAlmostZero) return kFALSE;
358 fParams[0] = (sX*sZZ-sZ*sXZ)/det;
359 fParams[1] = (sXZ*np-sZ*sX)/det;
360 //
361 fParams[2] = (sY*sZZ-sZ*sYZ)/det;
362 fParams[3] = (sYZ*np-sZ*sY)/det;
363 //
364 return kTRUE;
365 //
366}
367
368//____________________________________________________
369void AliITSTPArrayFit::SetParAxis(Int_t ax)
370{
371 // select the axis which will be used as a parameter for the line: longest baseline
372 if (ax>kZ) {
373 AliInfo(Form("Wrong axis choice: %d",ax));
374 fParAxis = -1;
375 }
376 fParAxis = ax;
377 if (ax>=0) {
378 fkAxID = fgkAxisID[ax];
379 fkAxCID = fgkAxisCID[ax];
380 }
381 else {
382 fkAxID = fkAxCID = 0;
383 }
384 //
385}
386
387//____________________________________________________
388Int_t AliITSTPArrayFit::ChoseParAxis() const
389{
390 // select the variable with largest base as a parameter
391 Double_t cmn[3]={1.e9,1.e9,1.e9},cmx[3]={-1.e9,-1.e9,-1.e9};
392 //
24391cd5 393 const float *coord[3] = {fkPoints->GetX(),fkPoints->GetY(),fkPoints->GetZ()};
6be22b3f 394 for (int i=fPntFirst;i<=fPntLast;i++) {
395 for (int j=3;j--;) {
396 Double_t val = coord[j][i];
397 if (cmn[j]>val) cmn[j] = val;
398 if (cmx[j]<val) cmx[j] = val;
399 }
400 }
401 //
402 int axis = kZ;
403 if (cmx[axis]-cmn[axis] < cmx[kX]-cmn[kX]) axis = kX;
404 if (cmx[axis]-cmn[axis] < cmx[kY]-cmn[kY]) axis = kY;
405 return axis;
406 //
407}
408
409//____________________________________________________
410Double_t AliITSTPArrayFit::GetPosition(Double_t *xyzPCA, const Double_t *xyz, const Double_t *covI) const
411{
412 // calculate the position of the track at PCA to xyz
413 Double_t t = GetParPCA(xyz,covI);
414 GetPosition(xyzPCA,t);
415 return t;
416}
417
418//____________________________________________________
ef24eb3b 419Double_t AliITSTPArrayFit::GetPosition(Double_t *xyzPCA, const AliTrackPoint *pntCovInv, Bool_t useErr) const
6be22b3f 420{
421 // calculate the position of the track at PCA to pntCovInv
422 // NOTE: the covariance matrix of the point must be inverted
ef24eb3b 423 double t;
424 double xyz[3] = {pntCovInv->GetX(),pntCovInv->GetY(),pntCovInv->GetZ()};
425 if (useErr) {
426 Double_t covI[6];;
427 for (int i=6;i--;) covI[i] = pntCovInv->GetCov()[i];
428 t = GetParPCA(xyz,covI);
429 }
430 else t = GetParPCA(xyz);
6be22b3f 431 GetPosition(xyzPCA,t);
432 return t;
433}
434
435//____________________________________________________
ef24eb3b 436void AliITSTPArrayFit::GetResiduals(Double_t *resPCA, const AliTrackPoint *pntCovInv, Bool_t useErr) const
6be22b3f 437{
438 // calculate the residuals of the track at PCA to pntCovInv
439 // NOTE: the covariance matrix of the point must be inverted
ef24eb3b 440 GetPosition(resPCA,pntCovInv,useErr);
6be22b3f 441 resPCA[0] -= pntCovInv->GetX();
442 resPCA[1] -= pntCovInv->GetY();
443 resPCA[2] -= pntCovInv->GetZ();
444}
445
446//____________________________________________________
447void AliITSTPArrayFit::GetResiduals(Double_t *resPCA, const Double_t *xyz, const Double_t *covI) const
448{
449 // calculate the residuals of the track at PCA to xyz
450 GetPosition(resPCA,xyz,covI);
451 resPCA[kX] -= xyz[kX];
452 resPCA[kY] -= xyz[kY];
453 resPCA[kZ] -= xyz[kZ];
454}
455
456//____________________________________________________
457Double_t AliITSTPArrayFit::GetParPCALine(const Double_t *xyz, const Double_t *covI) const
458{
459 // get parameter for the point with least weighted distance to the point
460 //
461 Double_t rhs,denom;
462 Double_t dx = fParams[kA0]-xyz[ fkAxID[kX] ];
463 Double_t dy = fParams[kA1]-xyz[ fkAxID[kY] ];
464 Double_t dz = -xyz[ fkAxID[kZ] ];
465 //
466 if (covI) {
467 Double_t tx = fParams[kB0]*covI[ fkAxCID[kXX] ] + fParams[kB1]*covI[ fkAxCID[kXY] ] + covI[ fkAxCID[kXZ] ];
468 Double_t ty = fParams[kB0]*covI[ fkAxCID[kXY] ] + fParams[kB1]*covI[ fkAxCID[kYY] ] + covI[ fkAxCID[kYZ] ];
469 Double_t tz = fParams[kB0]*covI[ fkAxCID[kXZ] ] + fParams[kB1]*covI[ fkAxCID[kYZ] ] + covI[ fkAxCID[kZZ] ];
470 rhs = tx*dx + ty*dy + tz*dz;
471 denom = -(fParams[kB0]*(covI[ fkAxCID[kXZ] ] + tx) + fParams[kB1]*(covI[ fkAxCID[kYZ] ] + ty) + covI[ fkAxCID[kZZ] ]);
472 }
473 else {
474 rhs = fParams[kB0]*dx + fParams[kB1]*dy + dz;
475 denom = -(fParams[kB0]*fParams[kB0] + fParams[kB1]*fParams[kB1] + 1);
476 }
477 //
478 return rhs/denom;
479 //
480}
481
482//____________________________________________________
483void AliITSTPArrayFit::GetDResDPosLine(Double_t *dXYZdP, /*const Double_t *xyz,*/ const Double_t *covI) const
484{
485 // calculate detivative of the PCA residuals vs point position and fill in user provide
486 // array in the format {dXdXp,dY/dXp,dZdXp, ... dXdZp,dYdZp,dZdZp}
487 //
488 Double_t dTdP[3];
489 GetDtDPosLine(dTdP, /*xyz,*/ covI); // derivative of the t-param over point position
490 //
491 for (int i=3;i--;) {
492 int var = fkAxID[i];
493 Double_t *curd = dXYZdP + var*3; // d/dCoord_i
494 curd[ fkAxID[kX] ] = fParams[kB0]*dTdP[var];
495 curd[ fkAxID[kY] ] = fParams[kB1]*dTdP[var];
496 curd[ fkAxID[kZ] ] = dTdP[var];
497 curd[ var ]-= 1.;
498 }
499 //
500}
501
502//____________________________________________________
503void AliITSTPArrayFit::GetDResDParamsLine(Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI) const
504{
505 // calculate detivative of the PCA residuals vs line parameters and fill in user provide
506 // array in the format {dXdP0,dYdP0,dZdP0, ... dXdPn,dYdPn,dZdPn}
507 //
508 Double_t dTdP[4];
509 Double_t t = GetDtDParamsLine(dTdP, xyz, covI); // derivative of the t-param over line params
510 //
511 Double_t *curd = dXYZdP + kA0*3; // d/dA0
512 curd[ fkAxID[kX] ] = fParams[kB0]*dTdP[kA0] + 1.;
513 curd[ fkAxID[kY] ] = fParams[kB1]*dTdP[kA0];
514 curd[ fkAxID[kZ] ] = dTdP[kA0];
515 //
516 curd = dXYZdP + kB0*3; // d/dB0
517 curd[ fkAxID[kX] ] = fParams[kB0]*dTdP[kB0] + t;
518 curd[ fkAxID[kY] ] = fParams[kB1]*dTdP[kB0];
519 curd[ fkAxID[kZ] ] = dTdP[kB0];
520 //
521 curd = dXYZdP + kA1*3; // d/dA1
522 curd[ fkAxID[kX] ] = fParams[kB0]*dTdP[kA1];
523 curd[ fkAxID[kY] ] = fParams[kB1]*dTdP[kA1] + 1.;
524 curd[ fkAxID[kZ] ] = dTdP[kA1];
525 //
526 curd = dXYZdP + kB1*3; // d/dB1
527 curd[ fkAxID[kX] ] = fParams[kB0]*dTdP[kB1];
528 curd[ fkAxID[kY] ] = fParams[kB1]*dTdP[kB1] + t;
529 curd[ fkAxID[kZ] ] = dTdP[kB1];
530 //
531}
532
533//____________________________________________________
534Double_t AliITSTPArrayFit::GetDtDParamsLine(Double_t *dtparam,const Double_t *xyz, const Double_t *covI) const
535{
536 // get t-param detivative over the parameters for the point with least weighted distance to the point
537 //
538 Double_t rhs,denom;
539 Double_t dx = fParams[kA0]-xyz[ fkAxID[kX] ];
540 Double_t dy = fParams[kA1]-xyz[ fkAxID[kY] ];
541 Double_t dz = -xyz[ fkAxID[kZ] ];
542 Double_t rhsDA0,rhsDA1,rhsDB0,rhsDB1,denDB0,denDB1;
543 //
544 if (covI) {
545 Double_t tx = fParams[kB0]*covI[ fkAxCID[kXX] ] + fParams[kB1]*covI[ fkAxCID[kXY] ] + covI[ fkAxCID[kXZ] ];
546 Double_t ty = fParams[kB0]*covI[ fkAxCID[kXY] ] + fParams[kB1]*covI[ fkAxCID[kYY] ] + covI[ fkAxCID[kYZ] ];
547 Double_t tz = fParams[kB0]*covI[ fkAxCID[kXZ] ] + fParams[kB1]*covI[ fkAxCID[kYZ] ] + covI[ fkAxCID[kZZ] ];
548 rhs = tx*dx + ty*dy + tz*dz;
549 denom = -(fParams[kB0]*(covI[ fkAxCID[kXZ] ] + tx) + fParams[kB1]*(covI[ fkAxCID[kYZ] ] + ty) + covI[ fkAxCID[kZZ] ]);
550 //
551 rhsDA0 = tx;
552 rhsDA1 = ty;
553 rhsDB0 = covI[ fkAxCID[kXX] ]*dx + covI[ fkAxCID[kXY] ]*dy + covI[ fkAxCID[kXZ] ]*dz;
554 rhsDB1 = covI[ fkAxCID[kXY] ]*dx + covI[ fkAxCID[kYY] ]*dy + covI[ fkAxCID[kYZ] ]*dz;
555 //
556 denDB0 = -(tx + tx);
557 denDB1 = -(ty + ty);
558 }
559 else {
560 rhs = fParams[kB0]*dx + fParams[kB1]*dy + dz;
561 denom = -(fParams[kB0]*fParams[kB0] + fParams[kB1]*fParams[kB1] + 1);
562 //
563 rhsDA0 = fParams[kB0];
564 rhsDB0 = dx;
565 rhsDA1 = fParams[kB1];
566 rhsDB1 = dy;
567 //
568 denDB0 = -(fParams[kB0]+fParams[kB0]);
569 denDB1 = -(fParams[kB1]+fParams[kB1]);
570 //
571 }
572 //
573 Double_t denom2 = denom*denom;
574 dtparam[kA0] = rhsDA0/denom; // denom does not depend on A0,A1
575 dtparam[kA1] = rhsDA1/denom;
576 dtparam[kB0] = rhsDB0/denom - rhs/denom2 * denDB0;
577 dtparam[kB1] = rhsDB1/denom - rhs/denom2 * denDB1;
578 //
579 return rhs/denom;
580}
581
582//____________________________________________________
583void AliITSTPArrayFit::GetDtDPosLine(Double_t *dtpos,/*const Double_t *xyz,*/ const Double_t *covI) const
584{
585 // get t-param detivative over the parameters for the point with least weighted distance to the point
586 //
587 // Double_t rhs;
588 // Double_t dx = fParams[kA0]-xyz[ fkAxID[kX] ];
589 // Double_t dy = fParams[kA1]-xyz[ fkAxID[kY] ];
590 // Double_t dz = -xyz[ fkAxID[kZ] ];
591 Double_t denom;
592 Double_t rhsDX,rhsDY,rhsDZ;
593 //
594 if (covI) {
595 Double_t tx = fParams[kB0]*covI[ fkAxCID[kXX] ] + fParams[kB1]*covI[ fkAxCID[kXY] ] + covI[ fkAxCID[kXZ] ];
596 Double_t ty = fParams[kB0]*covI[ fkAxCID[kXY] ] + fParams[kB1]*covI[ fkAxCID[kYY] ] + covI[ fkAxCID[kYZ] ];
597 Double_t tz = fParams[kB0]*covI[ fkAxCID[kXZ] ] + fParams[kB1]*covI[ fkAxCID[kYZ] ] + covI[ fkAxCID[kZZ] ];
598 // rhs = tx*dx + ty*dy + tz*dz;
599 denom = -(fParams[kB0]*(covI[ fkAxCID[kXZ] ] + tx) + fParams[kB1]*(covI[ fkAxCID[kYZ] ] + ty) + covI[ fkAxCID[kZZ] ]);
600 //
601 rhsDX = -tx;
602 rhsDY = -ty;
603 rhsDZ = -tz;
604 }
605 else {
606 // rhs = fParams[kB0]*dx + fParams[kB1]*dy + dz;
607 denom = -(fParams[kB0]*fParams[kB0] + fParams[kB1]*fParams[kB1] + 1);
608 //
609 rhsDX = -fParams[kB0];
610 rhsDY = -fParams[kB1];
611 rhsDZ = -1;
612 //
613 }
614 //
615 dtpos[ fkAxID[kX] ] = rhsDX/denom;
616 dtpos[ fkAxID[kY] ] = rhsDY/denom;
617 dtpos[ fkAxID[kZ] ] = rhsDZ/denom;
618 //
619 // return rhs/denom;
620}
621
622//____________________________________________________
623void AliITSTPArrayFit::GetDResDParamsLine(Double_t *dXYZdP, Int_t ipnt) const
624{
625 // calculate detivative of the PCA residuals vs line parameters and fill in user provide
626 // array in the format {dXdP0,dYdP0,dZdP0, ... dXdPn,dYdPn,dZdPn}
627 //
628 if (ipnt<fPntFirst || ipnt>fPntLast) {
629 AliError(Form("Attempt to access the point %d not in the fitted points [%d:%d]",ipnt,fPntFirst,fPntLast));
630 return;
631 }
632 GetDResDParamsLine(dXYZdP, GetPoint(ipnt) , IsCovIgnored() ? 0 : GetCovI(ipnt));
633}
634
635//____________________________________________________
636void AliITSTPArrayFit::GetDResDPosLine(Double_t *dXYZdP, Int_t ipnt) const
637{
638 // calculate detivative of the PCA residuals vs point position and fill in user provide
639 // array in the format {dXdXp,dY/dXp,dZdXp, ... dXdZp,dYdZp,dZdZp}
640 //
641 if (ipnt<fPntFirst || ipnt>fPntLast) {
642 AliError(Form("Attempt to access the point %d not in the fitted points [%d:%d]",ipnt,fPntFirst,fPntLast));
643 return;
644 }
645 GetDResDPosLine(dXYZdP,IsCovIgnored() ? 0 : GetCovI(ipnt));
646}
647
648//____________________________________________________
649void AliITSTPArrayFit::GetDResDParams(Double_t *dXYZdP, Int_t ipnt)
650{
651 // calculate detivative of the PCA residuals vs track parameters and fill in user provide
652 // array in the format {dXdP0,dYdP0,dZdP0, ... dXdPn,dYdPn,dZdPn}
653 //
654 if (ipnt<fPntFirst || ipnt>fPntLast) {
655 AliError(Form("Attempt to access the point %d not in the fitted points [%d:%d]",ipnt,fPntFirst,fPntLast));
656 return;
657 }
658 GetDResDParams(dXYZdP, GetPoint(ipnt) , IsCovIgnored() ? 0 : GetCovI(ipnt));
659}
660
661//____________________________________________________
662void AliITSTPArrayFit::GetDResDPos(Double_t *dXYZdP, Int_t ipnt)
663{
664 // calculate detivative of the PCA residuals vs point position and fill in user provide
665 // array in the format {dXdXp,dY/dXp,dZdXp, ... dXdZp,dYdZp,dZdZp}
666 //
667 if (ipnt<fPntFirst || ipnt>fPntLast) {
668 AliError(Form("Attempt to access the point %d not in the fitted points [%d:%d]",ipnt,fPntFirst,fPntLast));
669 return;
670 }
671 GetDResDPos(dXYZdP, GetPoint(ipnt), IsCovIgnored() ? 0 : GetCovI(ipnt));
672}
673
674//____________________________________________________
675void AliITSTPArrayFit::GetDResDParams(Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI)
676{
677 // get residual detivatives over the track parameters for the point with least weighted distance to the point
678 //
ef24eb3b 679 if (!IsHelix()) { // for the straight line calculate analytically
6be22b3f 680 GetDResDParamsLine(dXYZdP, xyz, covI);
681 return;
682 }
683 //
684 // calculate derivative numerically
685 const Double_t delta = 0.01;
686 Double_t xyzVar[4][3];
687 //
688 for (int ipar = 5;ipar--;) {
689 double sav = fParams[ipar];
690 fParams[ipar] -= delta;
691 GetPosition(xyzVar[0],xyz,covI);
692 fParams[ipar] += delta/2;
693 GetPosition(xyzVar[1],xyz,covI);
694 fParams[ipar] += delta;
695 GetPosition(xyzVar[2],xyz,covI);
696 fParams[ipar] += delta/2;
697 GetPosition(xyzVar[3],xyz,covI);
698 fParams[ipar] = sav; // restore
699 //
700 double *curd = dXYZdP + 3*ipar;
701 for (int i=3;i--;) curd[i] = (8.*(xyzVar[2][i]-xyzVar[1][i]) - (xyzVar[3][i]-xyzVar[0][i]))/6./delta;
702 }
703 //
704}
705
706
707//____________________________________________________
708void AliITSTPArrayFit::GetDResDPos(Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI)
709{
710 // get residuals detivative over the point position for the point with least weighted distance to the point
711 //
712
ef24eb3b 713 if (!IsHelix()) { // for the straight line calculate analytically
6be22b3f 714 GetDResDPosLine(dXYZdP, /*xyz,*/ covI);
715 return;
716 }
717 //
718 // calculate derivative numerically
719 const Double_t delta = 0.005;
720 Double_t xyzVar[4][3];
721 Double_t xyzv[3] = {xyz[0],xyz[1],xyz[2]};
722 //
723 for (int ipar = 3;ipar--;) {
724 double sav = xyzv[ipar];
725 xyzv[ipar] -= delta;
726 GetPosition(xyzVar[0],xyzv,covI);
727 xyzv[ipar] += delta/2;
728 GetPosition(xyzVar[1],xyzv,covI);
729 xyzv[ipar] += delta;
730 GetPosition(xyzVar[2],xyzv,covI);
731 xyzv[ipar] += delta/2;
732 GetPosition(xyzVar[3],xyzv,covI);
733 xyzv[ipar] = sav; // restore
734 //
735 double *curd = dXYZdP + 3*ipar;
736 for (int i=3;i--;) curd[i] = (8.*(xyzVar[2][i]-xyzVar[1][i]) - (xyzVar[3][i]-xyzVar[0][i]))/6./delta;
737 curd[ipar] -= 1.;
738 }
739 //
740}
741
742//________________________________________________________________________________________________________
743Double_t AliITSTPArrayFit::GetParPCAHelix(const Double_t* xyz, const Double_t* covI) const
744{
745 // find track parameter t (eq.2) corresponding to point of closest approach to xyz
746 //
747 Double_t phi = GetParPCACircle(xyz[kX],xyz[kY]);
748 Double_t cs = TMath::Cos(fParams[kPhi0]);
749 Double_t sn = TMath::Sin(fParams[kPhi0]);
750 Double_t xc = (fParams[kD0]+fParams[kR0])*cs;
751 Double_t yc = (fParams[kD0]+fParams[kR0])*sn;
752 Double_t dchi2,ddchi2;
753 //
754 Double_t dzD = -fParams[kR0]*fParams[kDip];
755 Double_t dphi = 0;
756 //
ef24eb3b 757 double rEps = 1e-5/TMath::Abs(fParams[kR0]); // dphi corresponding to 0.1 micron
758 if (rEps>fEps) rEps = fEps;
759 //
6be22b3f 760 int it=0;
761 do {
762 cs = TMath::Cos(phi + fParams[kPhi0]);
763 sn = TMath::Sin(phi + fParams[kPhi0]);
764 //
765 Double_t dxD = fParams[kR0]*sn;
766 Double_t dyD = -fParams[kR0]*cs;
767 Double_t dxDD = -dyD;
768 Double_t dyDD = dxD;
769 //
770 Double_t dx = xc - fParams[kR0]*cs - xyz[kX];
771 Double_t dy = yc - fParams[kR0]*sn - xyz[kY];
772 Double_t dz = fParams[kDZ] + dzD*phi- xyz[kZ];
773 //
774 if (covI) {
775 Double_t tx = dx*covI[kXX] + dy*covI[kXY] + dz*covI[kXZ];
776 Double_t ty = dx*covI[kXY] + dy*covI[kYY] + dz*covI[kYZ];
777 Double_t tz = dx*covI[kXZ] + dy*covI[kYZ] + dz*covI[kZZ];
778 //
779 Double_t ttx = dxD*covI[kXX] + dyD*covI[kXY] + dzD*covI[kXZ];
780 Double_t tty = dxD*covI[kXY] + dyD*covI[kYY] + dzD*covI[kYZ];
781 Double_t ttz = dxD*covI[kXZ] + dyD*covI[kYZ] + dzD*covI[kZZ];
782 //
783 // chi2 = dx*tx + dy*ty + dz*tz;
784 dchi2 = dxD*tx + dyD*ty + dzD*tz;
785 ddchi2 = dxDD*tx + dyDD*ty + dxD *ttx + dyD *tty + dzD *ttz;
786 //
787 }
788 else {
789 // chi2 = dx*dx + dy*dy + dz*dz;
790 dchi2 = dxD*dx + dyD*dy + dzD*dz;
791 ddchi2 = dxDD*dx + dyDD*dy + + dxD*dxD + dyD*dyD + dzD*dzD;
792 }
793 //
ef24eb3b 794 if (TMath::Abs(ddchi2)<fgkAlmostZero || TMath::Abs(dphi=dchi2/ddchi2)<rEps) break;
6be22b3f 795 phi -= dphi;
796 } while(++it<fMaxIter);
ef24eb3b 797
6be22b3f 798 //
799 return phi;
800}
801
802//________________________________________________________________________________________________________
803Double_t AliITSTPArrayFit::GetParPCACircle(Double_t x,Double_t y) const
804{
805 // find track parameter t (eq.2) corresponding to point on the circle with closest approach to x,y
806 //
807 Double_t r = fParams[kD0]+fParams[kR0];
808 Double_t t = TMath::ATan2( r*TMath::Sin(fParams[kPhi0])-y, r*TMath::Cos(fParams[kPhi0])-x ) - fParams[kPhi0];
809 if (fParams[kR0] < 0) t += TMath::Pi();
810 if (t > TMath::Pi()) t -= TMath::Pi()*2;
811 if (t <-TMath::Pi()) t += TMath::Pi()*2;
812 return t;
813}
814
815//________________________________________________________________________________________________________
816Double_t AliITSTPArrayFit::GetHelixParAtR(Double_t r) const
817{
818 // find helix parameter t (eq.2) corresponding to point on the circle of radius t
819 //
820 double gam = 1. - (r-fParams[kD0])*(r+fParams[kD0])/fParams[kR0]/(fParams[kD0]+fParams[kR0])/2.;
821 return (TMath::Abs(gam)>1) ? -1e9 : TMath::ACos(gam);
822}
823
824//________________________________________________________________________________________________________
825Double_t AliITSTPArrayFit::CalcChi2NDF() const
826{
827 // calculate fit chi2/ndf
828 Double_t chi2 = 0;
829 Double_t dr[3]; // residuals
830 //if (!IsFitDone()) return -1;
831 for (int ipnt=fPntFirst;ipnt<=fPntLast;ipnt++) {
832 GetResiduals(dr,ipnt);
833 Double_t* covI = GetCovI(ipnt);
834 chi2 += dr[kX]*(dr[kX]*covI[ kXX ]+dr[kY]*covI[ kXY ]+dr[kZ]*covI[ kXZ ])
835 + dr[kY]*(dr[kX]*covI[ kXY ]+dr[kY]*covI[ kYY ]+dr[kZ]*covI[ kYZ ])
836 + dr[kZ]*(dr[kX]*covI[ kXZ ]+dr[kY]*covI[ kYZ ]+dr[kZ]*covI[ kZZ ]);
837 }
838 int ndf = (fPntLast-fPntFirst+1)*3 - GetNParams();
839 chi2 /= ndf;
840 return chi2;
841}
842
843//________________________________________________________________________________________________________
844void AliITSTPArrayFit::GetResiduals(Double_t *res,Int_t ipnt) const
845{
846 // calculate residuals at point
847 if (ipnt<fPntFirst || ipnt>fPntLast) {
848 AliError(Form("Attempt to access the point %d not in the fitted points [%d:%d]",ipnt,fPntFirst,fPntLast));
849 return;
850 }
851 GetPosition(res,fCurT[ipnt]);
24391cd5 852 res[kX] -= fkPoints->GetX()[ipnt];
853 res[kY] -= fkPoints->GetY()[ipnt];
854 res[kZ] -= fkPoints->GetZ()[ipnt];
6be22b3f 855}
856
857//________________________________________________________________________________________________________
858void AliITSTPArrayFit::GetPosition(Double_t *xyz, Double_t t) const
859{
860 // calculate track position for parameter value t
ef24eb3b 861 if (IsHelix()) {
6be22b3f 862 //
863 Double_t rrho = fParams[kD0]+fParams[kR0];
864 Double_t xc = rrho*TMath::Cos(fParams[kPhi0]);
865 Double_t yc = rrho*TMath::Sin(fParams[kPhi0]);
866 Double_t r = fParams[kR0];
867 Double_t ze = 0;
868 //
869 if (IsELossON()) {
870 if (t>0) {
871 for (int i=fFirstPosT;i<fNElsPnt;i++) { // along the track direction
872 int indE = fElsId[i];
873 if ( t<fCurT[indE] ) break; // does not reach this layer on its way to t
874 xc += fElsDR[indE] * TMath::Cos(fParams[kPhi0] + fCurT[indE]);
875 yc += fElsDR[indE] * TMath::Sin(fParams[kPhi0] + fCurT[indE]);
876 ze += fElsDR[indE] * fCurT[indE];
877 r += fElsDR[indE];
878 //printf("ELoss@ %+.2e r:%+.3e got %+.3e\n",fCurT[indE],r,fElsDR[indE]);
879 }
880 } else {
881 for (int i=fFirstPosT;i--;) { // against the track direction
882 int indE = fElsId[i];
883 if ( t>=fCurT[indE] ) break; // does not reach this layer on its way to t
884 xc += fElsDR[indE] * TMath::Cos(fParams[kPhi0] + fCurT[indE]);
885 yc += fElsDR[indE] * TMath::Sin(fParams[kPhi0] + fCurT[indE]);
886 ze += fElsDR[indE] * fCurT[indE];
887 r += fElsDR[indE];
888 //printf("ELoss@ %+.2e r:%+.3e got %+.3e\n",fCurT[indE],r,fElsDR[indE]);
889 }
890 }
891 }
892 //
893 xyz[kZ] = fParams[kDZ] - fParams[kDip]*(t*r - ze);
894 //
895 t += fParams[kPhi0];
896 xyz[kX] = xc - r*TMath::Cos(t);
897 xyz[kY] = yc - r*TMath::Sin(t);
898 // printf("t: %+.3e xyz:%+.2e %+.2e %+.2e | R %+.6e -> %+.6e | sign %d\n",t-fParams[kPhi0],xyz[0],xyz[1],xyz[2],fParams[kR0],r,GetSignQB());
899 }
900 else {
901 xyz[ fkAxID[kX] ] = fParams[kA0] + fParams[kB0]*t;
902 xyz[ fkAxID[kY] ] = fParams[kA1] + fParams[kB1]*t;
903 xyz[ fParAxis ] = t;
904 }
905}
906
907//________________________________________________________________________________________________________
908void AliITSTPArrayFit::GetDirCos(Double_t *dircos, Double_t t) const
909{
910 // calculate track direction cosines for parameter value t
ef24eb3b 911 if (IsHelix()) {
6be22b3f 912 dircos[kZ] = -fParams[kDip];
913 t += fParams[kPhi0];
914 dircos[kX] = TMath::Sin(t);
915 dircos[kY] =-TMath::Cos(t);
916 double gam = TMath::Sign(1/TMath::Sqrt(dircos[kZ]*dircos[kZ]+dircos[kY]*dircos[kY]+dircos[kX]*dircos[kX]),fParams[kR0]);
917 for (int i=3;i--;) dircos[i] *= gam;
66214d86 918 if (GetSignQB()>0) for (int i=3;i--;) dircos[i] = -dircos[i]; // positive tracks move along decreasing t
6be22b3f 919 }
920 else {
921 double gam = 1/TMath::Sqrt( fParams[kB0]*fParams[kB0] + fParams[kB1]*fParams[kB1] + 1.);
922 dircos[ fkAxID[kX] ] = fParams[kB0]*gam;
923 dircos[ fkAxID[kY] ] = fParams[kB1]*gam;
924 dircos[ fParAxis ] = gam;
66214d86 925 // decide direction
926 if (IsTypeCollision()) {
927 static double xyzF[3],xyzL[3];
928 GetPosition(xyzF,fPntFirst);
929 GetPosition(xyzL,fPntLast);
930 double dif = fCurT[fPntLast] - fCurT[fPntFirst];
931 double dr = (xyzL[kX]-xyzF[kX])*(xyzL[kX]+xyzF[kX]) + (xyzL[kY]-xyzF[kY])*(xyzL[kY]+xyzF[kY]);
932 if (dr*dif<0) for (int i=3;i--;) dircos[i] = -dircos[i]; // with increasing t the tracks comes closer to origin
933 }
934 else if (dircos[kY]>0) for (int i=3;i--;) dircos[i] = -dircos[i]; // cosmic tracks have negative angle to Y axis
6be22b3f 935 }
66214d86 936 //
6be22b3f 937}
938
939//________________________________________________________________________________________________________
940Double_t AliITSTPArrayFit::GetMachinePrec()
941{
942 // estimate machine precision
943 Double_t eps=1.0,a;
944 do { a = 1. + (eps=eps/2.0); } while(a>1.);
945 return TMath::Abs(2.*eps);
946}
947
ef24eb3b 948//________________________________________________________________________________________________________
949Bool_t AliITSTPArrayFit::FitHelixCrude(Int_t extQ)
950{
951 // crude estimate of helix parameters, w/o errors and Eloss.
952 // Fast Riemann fit: Comp.Phy.Comm.131 (2000) 95
953 //
954 // if charge is not imposed (extQ==0) then it will be determined from the collision type
955 //
956 static TArrayD arrU,arrV,arrW;
957 double *parrW,*parrU,*parrV;
958 Bool_t eloss = IsELossON();
959 //
960 int np = fPntLast - fPntFirst + 1;
961 if (np<3) { AliError("At least 3 points are needed for helix fit"); return kFALSE; }
962 //
963 const float *x=fkPoints->GetX(),*y=fkPoints->GetY(),*z=fkPoints->GetZ(),*cov=fkPoints->GetCov();
964 //
965 if (fPntLast>arrU.GetSize()) {
966 arrU.Set(2*fPntLast);
967 arrV.Set(2*fPntLast);
968 arrW.Set(2*fPntLast);
969 }
970 parrU = arrU.GetArray();
971 parrV = arrV.GetArray();
972 parrW = arrW.GetArray();
973 //
974 double uav=0,vav=0,wav=0,muu=0,muv=0,muw=0,mvv=0,mvw=0,mww=0;
975 int minRId = fPntFirst;
976 //
977 // get points span
978 double xmn=1e9,xmx=-1e9, ymn=1e9,ymx=-1e9;
979 for (int i=fPntFirst;i<=fPntLast;i++) {
980 parrW[i] = x[i]*x[i]+y[i]*y[i];
981 if (parrW[i]<parrW[minRId]) minRId = i; // point closest to origin
982 if (xmn>x[i]) xmn = x[i];
983 if (xmx<x[i]) xmx = x[i];
984 if (ymn>y[i]) ymn = y[i];
985 if (ymx<y[i]) ymx = y[i];
986 }
987 int minRId1 = minRId>fPntFirst ? fPntFirst:fPntFirst+1;
988 for (int i=fPntFirst;i<=fPntLast;i++) if (parrW[i]<parrW[minRId1] && i!=minRId) minRId1 = i;
989 //
990 double xshift = (xmx+xmn)/2 + 10*(ymx-ymn); // shift origin to have uniform weights
991 double yshift = (ymx+ymn)/2 - 10*(xmx-xmn);
992 // printf("X: %+e %+e Y: %+e %+e | shift: %+e %+e\n",xmn,xmx,ymn,ymx,xshift,yshift);
993 //
994 for (int i=fPntFirst;i<=fPntLast;i++) {
995 double xs = x[i] - xshift;
996 double ys = y[i] - yshift;
997 double w = xs*xs + ys*ys;
998 double scl = 1./(1.+w);
999 int i0 = i-fPntFirst;
1000 wav += parrW[i0] = w*scl;
1001 uav += parrU[i0] = xs*scl;
1002 vav += parrV[i0] = ys*scl;
1003 }
1004 uav /= np; vav /= np; wav /= np;
1005 //
1006 for (int i=fPntFirst;i<=fPntLast;i++) {
1007 //
1008 // point next to closest
1009 int i0 = i-fPntFirst;
1010 if (parrW[i0]<parrW[minRId1-fPntFirst] && i!=minRId) minRId1 = i;
1011 double u = parrU[i0] - uav;
1012 double v = parrV[i0] - vav;
1013 double w = parrW[i0] - wav;
1014 muu += u*u;
1015 muv += u*v;
1016 muw += u*w;
1017 mvv += v*v;
1018 mvw += v*w;
1019 mww += w*w;
1020 }
1021 muu/=np; muv/=np; muw/=np; mvv/=np; mvw/=np; mww/=np;
1022 //
1023 // find eigenvalues:
1024 double trace3 = (muu + mvv + mww)/3.;
1025 double muut = muu-trace3;
1026 double mvvt = mvv-trace3;
1027 double mwwt = mww-trace3;
1028 double q = (muut*(mvvt*mwwt-mvw*mvw) - muv*(muv*mwwt-mvw*muw) + muw*(muv*mvw-mvvt*muw))/2;
1029 double p = (muut*muut+mvvt*mvvt+mwwt*mwwt+2*(muv*muv+muw*muw+mvw*mvw))/6;
1030 double dfpp = p*p*p-q*q;
1031 dfpp = dfpp>0 ? TMath::Sqrt(dfpp)/q : 0;
1032 double ph = TMath::ATan( dfpp )/3.;
1033 if (ph<0) ph += TMath::Pi()/3;
1034 p = p>0 ? TMath::Sqrt(p) : 0;
1035 const double kSqrt3 = 1.73205080;
1036 double snp = TMath::Sin(ph);
1037 double csp = TMath::Cos(ph);
1038 // double eg1 = trace3 + 2*p*csp;
1039 double eg2 = trace3 - p*(csp+kSqrt3*snp); // smallest one
1040 // double eg3 = trace3 - p*(csp-kSqrt3*snp);
1041 // eigenvector for min.eigenvalue
1042 muut = muu-eg2;
1043 mvvt = mvv-eg2;
1044 mwwt = mww-eg2;
1045 double n0 = muv*mvw-muw*mvvt;
1046 double n1 = muv*muw-mvw*muut;
1047 double n2 = muut*mvvt-muv*muv;
1048 // normalize to largest one
1049 double nrm = TMath::Abs(n0);
1050 if (nrm<TMath::Abs(n1)) nrm = TMath::Abs(n1);
1051 if (nrm<TMath::Abs(n2)) nrm = TMath::Abs(n2);
1052 n0/=nrm; n1/=nrm; n2/=nrm;
1053 //
1054 double cpar = -(uav*n0 + vav*n1 + wav*n2);
1055 double xc = -n0/(cpar+n2)/2 + xshift;
1056 double yc = -n1/(cpar+n2)/2 + yshift;
1057 double rad = TMath::Sqrt(n0*n0+n1*n1-4*cpar*(cpar+n2))/2./TMath::Abs(cpar+n2);
1058 //
1059 // printf("Rad: %+e xc: %+e yc: %+e | X0: %+e Y0: %+e | X1: %+e Y1: %+e\n",rad,xc,yc, x[minRId],y[minRId],x[minRId1],y[minRId1]);
1060
1061 // linear circle fit --------------------------------------------------- <<<
1062 //
1063 // decide sign(Q*B) and fill cicrle parameters ------------------------- >>>
1064 int sqb;
1065 if (extQ) {
1066 SetCharge(extQ);
1067 sqb = fBz<0 ? -GetCharge():GetCharge();
1068 }
1069 else {
1070 // determine the charge from the collision type and field sign
1071 // the negative Q*B will have positive Vc x dir product Z component
1072 // with Vc={-xc,-yc} : vector from circle center to the origin
1073 // and V0 - track direction vector (take {0,-1,1} for cosmics)
1074 // If Bz is not provided, assume positive Bz
1075 if ( IsTypeCosmics() ) sqb = xc>0 ? -1:1;
1076 else {
1077 // track direction vector as a - diference between the closest and the next to closest to origin points
1078 double v0X = x[minRId1] - x[minRId];
1079 double v0Y = y[minRId1] - y[minRId];
1080 sqb = (yc*v0X - xc*v0Y)<0 ? -1:1;
1081 }
1082 SetCharge( fBz<0 ? -sqb : sqb);
1083 }
1084 //
1085 Double_t phi = TMath::ATan2(yc,xc);
1086 if (sqb<0) phi += TMath::Pi();
1087 if (phi > TMath::Pi()) phi -= 2.*TMath::Pi();
1088 else if (phi <-TMath::Pi()) phi += 2.*TMath::Pi();
1089 fParams[kPhi0] = phi;
1090 fParams[kR0] = sqb<0 ? -rad:rad;
1091 fParams[kD0] = xc*TMath::Cos(phi) + yc*TMath::Sin(phi) - fParams[kR0];
1092 //
1093 // decide sign(Q*B) and fill cicrle parameters ------------------------- <<<
1094 //
1095 // find z-offset and dip + the parameter t of closest approach to hits - >>>
1096 //
1097 UInt_t hitLrPos=0; // pattern of hit layers at pos
1098 UInt_t hitLrNeg=0; // and negative t's
1099
1100 Double_t ss=0,st=0,sz=0,stt=0,szt=0;
1101 for (int i=fPntFirst;i<=fPntLast;i++) {
1102 //
1103 Double_t ze2 = cov[i*6 + kZZ];
1104 Double_t t = TMath::ATan2(yc-y[i],xc-x[i]) - fParams[kPhi0]; // angle at measured z
1105 if (fParams[kR0]<0) t += TMath::Pi();
1106 if (t > TMath::Pi()) t -= TMath::Pi()*2;
1107 else if (t <-TMath::Pi()) t += TMath::Pi()*2;
1108 if (ze2<fgkAlmostZero) ze2 = 1E-8;
1109 ze2 = 1./ze2;
1110 ss += ze2;
1111 st += t*ze2;
1112 stt+= t*t*ze2;
1113 sz += z[i]*ze2;
1114 szt+= z[i]*t*ze2;
1115 //
1116 fCurT[i] = t; // parameter of the closest approach to the point
1117 // printf("%d %+e %+e %+e %+e\n",i,x[i],y[i],z[i],t);
1118 if (eloss) {
1119 double r = TMath::Sqrt(x[i]*x[i]+y[i]*y[i]);
1120 int lr;
1121 for (lr=kMaxLrITS;lr--;) if ( IsZero(r-fgkRLayITS[ lr ],1.) ) break;
1122 if (lr<kMaxLrITS) {
1123 if (t>0) hitLrPos |= (1<<lr); // set bit of the layer
1124 else hitLrNeg |= (1<<lr); // set bit of the layer
1125 }
1126 }
1127 }
1128 double det = ss*stt - st*st;
1129 if (TMath::Abs(det)<fgkAlmostZero) { // no Z dependence
1130 fParams[kDZ] = sz/ss;
1131 fParams[kDip] = 0;
1132 }
1133 else {
1134 fParams[kDZ] = (sz*stt-st*szt)/det;
1135 fParams[kDip] = -(ss*szt-st*sz)/det/fParams[kR0];
1136 }
1137 //
1138 // find z-offset and dip + the parameter t of closest approach to hits - <<<
1139 //
1140 // fill info needed to account for ELoss ------------------------------- >>>
1141 if (eloss) {
1142 fNElsPnt = fPntLast - fPntFirst + 1;
1143 //
1144 // to account for the energy loss in the passive volumes, calculate the relevant t-parameters
1145 double* tcur = fCurT + fPntFirst;
1146 double* ecur = fElsDR+ fPntFirst;
1147 //
1148 for (int ilp=3;ilp--;) {
1149 int id = fgkPassivLrITS[ilp];
1150 double tp = GetHelixParAtR( fgkRLayITS[ id ] );
1151 if (tp<0) continue; // does not hit this radius
1152 //
1153 tcur[fNElsPnt] = GetSignQB()>0 ? -tp : tp;
1154 ecur[fNElsPnt] = fgRhoLITS[ id ];
1155 fNElsPnt++;
1156 // printf("Passive on lr %d %+e\n",ilp,tcur[fNElsPnt-1]);
1157 //
1158 if (IsTypeCosmics() && !IsZero(tp)) { // 2 crossings for cosmics
1159 tcur[fNElsPnt] = -tcur[fNElsPnt-1];
1160 ecur[fNElsPnt] = ecur[fNElsPnt-1];
1161 fNElsPnt++;
1162 //printf("Passive* on lr %d %+e\n",ilp,-tcur[fNElsPnt-1]);
1163 }
1164 //
1165 }
1166 // check if some active layers did not miss the hit, treat them as passive
1167 for (int ilp=6;ilp--;) {
1168 int id = fgkActiveLrITS[ilp];
1169 double tp = GetHelixParAtR( fgkRLayITS[ id ] );
1170 if (tp<0) continue; // does not hit this radius
1171 //
1172 if ( (GetSignQB()>0||IsTypeCosmics()) && !(hitLrNeg & (1<<id)) ) {
1173 tcur[fNElsPnt] = -tp;
1174 ecur[fNElsPnt] = fgRhoLITS[ id ];
1175 fNElsPnt++;
1176 //printf("Missed on lr %d %+e\n",ilp,-tp);
1177 }
1178 //
1179 if ( (GetSignQB()<0||IsTypeCosmics()) && !(hitLrPos & (1<<id)) ) {
1180 tcur[fNElsPnt] = tp;
1181 ecur[fNElsPnt] = fgRhoLITS[ id ];
1182 fNElsPnt++;
1183 //printf("Missed* on lr %d %e\n",ilp,tp);
1184 }
1185 }
1186 //
1187 TMath::Sort(fNElsPnt,fCurT+fPntFirst,fElsId,kFALSE); // index e-loss points in increasing order
1188 // find the position of smallest positive t-param
1189 for (fFirstPosT=0;fFirstPosT<fNElsPnt;fFirstPosT++) if (fCurT[ fElsId[ fFirstPosT ] ]>0) break;
1190 //
1191 Double_t cdip = 1./TMath::Sqrt(1.+fParams[kDip]*fParams[kDip]);
1192 Double_t ptot = TMath::Abs(fParams[kR0]*fgkCQConv*fBz/cdip); // momentum and energy
1193 Double_t etot = TMath::Sqrt(ptot*ptot + fMass*fMass); // in the point of closest approach to beam
1194 Double_t normS[3];
1195 //
1196 // Positive t-params: along the track direction for negative track, against for positive
1197 Double_t pcur = ptot, ecurr = etot;
1198 for (int ip=fFirstPosT;ip<fNElsPnt;ip++) {
1199 int tID = fElsId[ip];
1200 Double_t t = fCurT[ tID ];
1201 //
1202 if (tID>fPntLast) { // this is not a hit layer but passive layer
1203 double php = TMath::ATan2(yc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t),
1204 xc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t));
1205 normS[0] = -TMath::Cos(php); // normal to the cylinder at intersection point
1206 normS[1] = -TMath::Sin(php);
1207 normS[2] = 0;
1208 }
1209 else GetNormal(normS,fkPoints->GetCov()+tID*6); // vector normal to hit module
1210 fElsDR[tID] = GetDRofELoss(t,cdip,fElsDR[tID],normS,ptot,etot);
1211 }
1212 //
1213 // negaive t-params: against the track direction for negative track, along for positive
1214 pcur = ptot;
1215 ecurr = etot;
1216 for (int ip=fFirstPosT;ip--;) {
1217 int tID = fElsId[ip];
1218 Double_t t = fCurT[ tID ];
1219 //
1220 if (tID>fPntLast) { // this is not a hit layer but passive layer
1221 double php = TMath::ATan2(yc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t),
1222 xc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t));
1223 normS[0] = -TMath::Cos(php); // normal to the cylinder at intersection point
1224 normS[1] = -TMath::Sin(php);
1225 normS[2] = 0;
1226 }
1227 else GetNormal(normS,fkPoints->GetCov()+tID*6); // vector normal to hit module
1228 //
1229 fElsDR[tID] = GetDRofELoss(t,cdip,fElsDR[tID],normS,ptot,etot);
1230 }
1231 }
1232 // fill info needed to account for ELoss ------------------------------- <<<
1233 //
1234 return kTRUE;
1235}
1236
1237/*
6be22b3f 1238//________________________________________________________________________________________________________
1239Bool_t AliITSTPArrayFit::FitHelixCrude(Int_t extQ)
1240{
1241 // crude estimate of helix parameters, w/o errors and Eloss.
1242 // 1st fit the circle (R,xc,yc) by minimizing
1243 // chi2 = sum{ (bx*xi + by*yi + xi^2+yi^2 + rho)^2 } vs bx,by,rho
1244 // with bx = -2*xc, by = -2*yc , rho = xc^2+yc^2 - R2
1245 //
1246 // if charge is not imposed (extQ==0) then it will be determined from the collision type
1247 //
1248 Bool_t eloss = IsELossON();
1249 //
1250 int np = fPntLast - fPntFirst + 1;
1d06ac63 1251 if (np<3) { AliError("At least 3 points are needed for helix fit"); return kFALSE; }
6be22b3f 1252 //
24391cd5 1253 const float *x=fkPoints->GetX(),*y=fkPoints->GetY(),*z=fkPoints->GetZ(),*cov=fkPoints->GetCov();
6be22b3f 1254 //
1255 // linear circle fit --------------------------------------------------- >>>
1256 Double_t sxx=0,sxy=0,syy=0,sx=0,sy=0,rhs0=0,rhs1=0,rhs2=0,minR=1E9;
1257 int minRId = 0;
1258 for (int i=fPntFirst;i<=fPntLast;i++) {
1259 Double_t xx = x[i]*x[i];
1260 Double_t yy = y[i]*y[i];
1261 Double_t xy = x[i]*y[i];
1262 Double_t xxyy = xx + yy;
1263 //
1264 sxx += xx;
1265 sxy += xy;
1266 syy += yy;
1267 sx += x[i];
1268 sy += y[i];
1269 //
1270 rhs0 -= xxyy*x[i];
1271 rhs1 -= xxyy*y[i];
1272 rhs2 -= xxyy;
1273 //
1274 // remember Id of the point closest to origin, to determine the charge
1275 if (xxyy<minR) { minR = xxyy; minRId = i; }
1276 //
1277 if (eloss) { // find layer id
24391cd5 1278 int lrid,volid = fkPoints->GetVolumeID()[i];
1279 if (volid>0) lrid = fgkActiveLrITS[AliGeomManager::VolUIDToLayer(fkPoints->GetVolumeID()[i])-1];
6be22b3f 1280 else { // missing layer info, find from radius
1281 double r = TMath::Sqrt(xxyy);
1282 for (lrid=kMaxLrITS;lrid--;) if ( IsZero(r-fgkRLayITS[ lrid ],1.) ) break;
1283 }
1284 fElsDR[i] = (lrid>=0 && lrid<kMaxLrITS) ? fgRhoLITS[ lrid ] : 0; // eloss for normal track
1285 }
1286 //
1287 }
1288 //
1289 Double_t mn00 = syy*np-sy*sy;
1290 Double_t mn01 = sxy*np-sy*sx;
1291 Double_t mn02 = sxy*sy-syy*sx;
1292 Double_t det = sxx*mn00 - sxy*mn01 + sx*mn02;
1293 if (TMath::Abs(det)<fgkAlmostZero) return kFALSE;
1294 //
1295 Double_t mn11 = sxx*np-sx*sx;
1296 Double_t mn12 = sxx*sy-sxy*sx;
1297 Double_t mn22 = sxx*syy-sxy*sxy;
1298 //
1299 Double_t mi00 = mn00/det;
1300 Double_t mi01 = -mn01/det;
1301 Double_t mi02 = mn02/det;
1302 Double_t mi11 = mn11/det;
1303 Double_t mi12 = -mn12/det;
1304 Double_t mi22 = mn22/det;
1305 //
1306 Double_t xc = -(rhs0*mi00 + rhs1*mi01 + rhs2*mi02)/2;
1307 Double_t yc = -(rhs0*mi01 + rhs1*mi11 + rhs2*mi12)/2;
1308 Double_t rho2 = (rhs0*mi02 + rhs1*mi12 + rhs2*mi22);
ef24eb3b 1309
1310 //
1311 // check
1312 double bx = -2*xc;
1313 double by = -2*yc;
1314 double sm0=0,sm1=0;
1315 for (int i=fPntFirst;i<=fPntLast;i++) {
1316 double dst = bx*x[i]+by*y[i]+x[i]*x[i]+y[i]*y[i]+rho2;
1317 sm0 += dst;
1318 sm1 += dst*dst;
1319 printf("Point %d: %+e %+e %+e\n",i,dst,sm0,sm1);
1320 }
1321
6be22b3f 1322 //
1d06ac63 1323 Double_t rad = xc*xc + yc*yc - rho2;
6be22b3f 1324 rad = (rad>fgkAlmostZero) ? (TMath::Sqrt(rad)):fgkAlmostZero;
1325 //
1326 // printf("Rad: %+e xc: %+e yc: %+e\n",rad,xc,yc);
1d06ac63 1327
6be22b3f 1328 // linear circle fit --------------------------------------------------- <<<
1329 //
1330 // decide sign(Q*B) and fill cicrle parameters ------------------------- >>>
1331 int sqb;
1332 if (extQ) {
1333 SetCharge(extQ);
1334 sqb = fBz<0 ? -GetCharge():GetCharge();
1335 }
1336 else {
1337 // determine the charge from the collision type and field sign
1338 // the negative Q*B will have positive Vc x V0 product Z component
1339 // with Vc={-xc,-yc} : vector from circle center to the origin
1340 // and V0 - track direction vector (take {0,-1,1} for cosmics)
1341 // If Bz is not provided, assume positive Bz
1342 sqb = ( IsTypeCosmics() ? xc:(yc*x[minRId]-xc*y[minRId]) ) > 0 ? -1:1;
1343 SetCharge( fBz<0 ? -sqb : sqb);
1344 }
1345 //
6be22b3f 1346 Double_t phi = TMath::ATan2(yc,xc);
1347 if (sqb<0) phi += TMath::Pi();
1348 if (phi > TMath::Pi()) phi -= 2.*TMath::Pi();
1349 else if (phi <-TMath::Pi()) phi += 2.*TMath::Pi();
1350 fParams[kPhi0] = phi;
1351 fParams[kR0] = sqb<0 ? -rad:rad;
1d06ac63 1352 fParams[kD0] = xc*TMath::Cos(phi) + yc*TMath::Sin(phi) - fParams[kR0];
6be22b3f 1353 //
1354 // decide sign(Q*B) and fill cicrle parameters ------------------------- <<<
1355 //
1356 // find z-offset and dip + the parameter t of closest approach to hits - >>>
1357 //
1358 UInt_t hitLrPos=0; // pattern of hit layers at pos
1359 UInt_t hitLrNeg=0; // and negative t's
1360
1361 Double_t ss=0,st=0,sz=0,stt=0,szt=0;
1362 for (int i=fPntFirst;i<=fPntLast;i++) {
1363 //
1364 Double_t ze2 = cov[i*6 + kZZ];
1365 Double_t t = TMath::ATan2(yc-y[i],xc-x[i]) - fParams[kPhi0]; // angle at measured z
1366 if (fParams[kR0]<0) t += TMath::Pi();
1367 if (t > TMath::Pi()) t -= TMath::Pi()*2;
1368 else if (t <-TMath::Pi()) t += TMath::Pi()*2;
1369 if (ze2<fgkAlmostZero) ze2 = 1E-8;
1370 ze2 = 1./ze2;
1371 ss += ze2;
1372 st += t*ze2;
1373 stt+= t*t*ze2;
1374 sz += z[i]*ze2;
1375 szt+= z[i]*t*ze2;
1376 //
1377 fCurT[i] = t; // parameter of the closest approach to the point
1378 // printf("%d %+e %+e %+e %+e\n",i,x[i],y[i],z[i],t);
1379 if (eloss) {
1380 double r = TMath::Sqrt(x[i]*x[i]+y[i]*y[i]);
1381 int lr;
1382 for (lr=kMaxLrITS;lr--;) if ( IsZero(r-fgkRLayITS[ lr ],1.) ) break;
1383 if (lr<kMaxLrITS) {
1384 if (t>0) hitLrPos |= (1<<lr); // set bit of the layer
1385 else hitLrNeg |= (1<<lr); // set bit of the layer
1386 }
1387 }
1388 }
1389 det = ss*stt - st*st;
1390 if (TMath::Abs(det)<fgkAlmostZero) { // no Z dependence
1391 fParams[kDZ] = sz/ss;
1392 fParams[kDip] = 0;
1393 }
1394 else {
1395 fParams[kDZ] = (sz*stt-st*szt)/det;
1396 fParams[kDip] = -(ss*szt-st*sz)/det/fParams[kR0];
1397 }
1398 //
1399 // find z-offset and dip + the parameter t of closest approach to hits - <<<
1400 //
1401 // fill info needed to account for ELoss ------------------------------- >>>
1402 if (eloss) {
1403 fNElsPnt = fPntLast - fPntFirst + 1;
1404 //
1405 // to account for the energy loss in the passive volumes, calculate the relevant t-parameters
1406 double* tcur = fCurT + fPntFirst;
1407 double* ecur = fElsDR+ fPntFirst;
1408 //
1409 for (int ilp=3;ilp--;) {
1410 int id = fgkPassivLrITS[ilp];
1411 double tp = GetHelixParAtR( fgkRLayITS[ id ] );
1412 if (tp<0) continue; // does not hit this radius
1413 //
1414 tcur[fNElsPnt] = GetSignQB()>0 ? -tp : tp;
1415 ecur[fNElsPnt] = fgRhoLITS[ id ];
1416 fNElsPnt++;
1417 // printf("Passive on lr %d %+e\n",ilp,tcur[fNElsPnt-1]);
1418 //
1419 if (IsTypeCosmics() && !IsZero(tp)) { // 2 crossings for cosmics
1420 tcur[fNElsPnt] = -tcur[fNElsPnt-1];
1421 ecur[fNElsPnt] = ecur[fNElsPnt-1];
1422 fNElsPnt++;
1423 //printf("Passive* on lr %d %+e\n",ilp,-tcur[fNElsPnt-1]);
1424 }
1425 //
1426 }
1427 // check if some active layers did not miss the hit, treat them as passive
1428 for (int ilp=6;ilp--;) {
1429 int id = fgkActiveLrITS[ilp];
1430 double tp = GetHelixParAtR( fgkRLayITS[ id ] );
1431 if (tp<0) continue; // does not hit this radius
1432 //
1433 if ( (GetSignQB()>0||IsTypeCosmics()) && !(hitLrNeg & (1<<id)) ) {
1434 tcur[fNElsPnt] = -tp;
1435 ecur[fNElsPnt] = fgRhoLITS[ id ];
1436 fNElsPnt++;
1437 //printf("Missed on lr %d %+e\n",ilp,-tp);
1438 }
1439 //
1440 if ( (GetSignQB()<0||IsTypeCosmics()) && !(hitLrPos & (1<<id)) ) {
1441 tcur[fNElsPnt] = tp;
1442 ecur[fNElsPnt] = fgRhoLITS[ id ];
1443 fNElsPnt++;
1444 //printf("Missed* on lr %d %e\n",ilp,tp);
1445 }
1446 }
1447 //
1448 TMath::Sort(fNElsPnt,fCurT+fPntFirst,fElsId,kFALSE); // index e-loss points in increasing order
1449 // find the position of smallest positive t-param
1450 for (fFirstPosT=0;fFirstPosT<fNElsPnt;fFirstPosT++) if (fCurT[ fElsId[ fFirstPosT ] ]>0) break;
1451 //
1452 Double_t cdip = 1./TMath::Sqrt(1.+fParams[kDip]*fParams[kDip]);
1453 Double_t ptot = TMath::Abs(fParams[kR0]*fgkCQConv*fBz/cdip); // momentum and energy
1454 Double_t etot = TMath::Sqrt(ptot*ptot + fMass*fMass); // in the point of closest approach to beam
1455 Double_t normS[3];
1456 //
1457 // Positive t-params: along the track direction for negative track, against for positive
1458 Double_t pcur = ptot, ecurr = etot;
1459 for (int ip=fFirstPosT;ip<fNElsPnt;ip++) {
1460 int tID = fElsId[ip];
1461 Double_t t = fCurT[ tID ];
1462 //
1463 if (tID>fPntLast) { // this is not a hit layer but passive layer
1464 double php = TMath::ATan2(yc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t),
1465 xc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t));
1466 normS[0] = -TMath::Cos(php); // normal to the cylinder at intersection point
1467 normS[1] = -TMath::Sin(php);
1468 normS[2] = 0;
1469 }
24391cd5 1470 else GetNormal(normS,fkPoints->GetCov()+tID*6); // vector normal to hit module
6be22b3f 1471 fElsDR[tID] = GetDRofELoss(t,cdip,fElsDR[tID],normS,ptot,etot);
1472 }
1473 //
1474 // negaive t-params: against the track direction for negative track, along for positive
1475 pcur = ptot;
1476 ecurr = etot;
1477 for (int ip=fFirstPosT;ip--;) {
1478 int tID = fElsId[ip];
1479 Double_t t = fCurT[ tID ];
1480 //
1481 if (tID>fPntLast) { // this is not a hit layer but passive layer
1482 double php = TMath::ATan2(yc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t),
1483 xc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t));
1484 normS[0] = -TMath::Cos(php); // normal to the cylinder at intersection point
1485 normS[1] = -TMath::Sin(php);
1486 normS[2] = 0;
1487 }
24391cd5 1488 else GetNormal(normS,fkPoints->GetCov()+tID*6); // vector normal to hit module
6be22b3f 1489 //
1490 fElsDR[tID] = GetDRofELoss(t,cdip,fElsDR[tID],normS,ptot,etot);
1491 }
1492 }
1493 // fill info needed to account for ELoss ------------------------------- <<<
1494 //
1495 return kTRUE;
1496}
ef24eb3b 1497*/
6be22b3f 1498//____________________________________________________
1499Double_t AliITSTPArrayFit::FitHelix(Int_t extQ, Double_t extPT,Double_t extPTerr)
1500{
1501 // fit by helix accounting for the errors of all coordinates (and energy loss if requested)
1502 //
1503 // If extQ is non-0, its sign is imposed as a charge of the track
1504 // If extPT>0 and extPTerr>=0, constrain to measured tr.momentum PT
1505 // with corresponding error (err=0 -> rel.err=1e-6)
1506 //
1507 double chiprev=1e99;
1508 //const Double_t kMaxTEffect = 1E-6;
1509 Double_t dXYZdGlo[3*5],dXYZdLoc[3],xyzRes[3];
1510 //
1511 SetFitDone(kFALSE);
1512 fChi2NDF = -1;
1513 //
1514 if (!FitHelixCrude(extQ)) return -1; // get initial estimate, ignoring the errors
1515 //
1d06ac63 1516 if (TMath::Abs(fParams[kR0])>fMaxRforHelix && extPT<=0) {
1517 fSwitch2Line = kTRUE;
1518 return FitLine();
1519 }
1520 //
1521 //
6be22b3f 1522 if (!IsCovInv()) InvertPointsCovMat(); // prepare inverted errors
1523 if (!fParSol) fParSol = new AliParamSolver(5);
1524 fParSol->SetNGlobal(5);
1525 //
1526 // printf("-1 | %+.2e %+.2e %+.2e %+.2e %+.2e | chi2: %+.4e\n",fParams[0],fParams[1],fParams[2],fParams[3],fParams[4],CalcChi2NDF());
1527 int iter = 0;
1528 fChi2NDF = 1e99;
1529 Bool_t converged = kFALSE;
1530 while(iter<fMaxIter) {
1531 chiprev = fChi2NDF;
1532 fParSol->Clear();
1533 for (int ip=fPntFirst;ip<=fPntLast;ip++) {
1534 //
1535 GetResiduals(xyzRes, ip); // current residuals at point ip
1536 Double_t rrho = fParams[kR0]+fParams[kD0];
1537 Double_t cs0 = TMath::Cos(fParams[kPhi0]);
1538 Double_t sn0 = TMath::Sin(fParams[kPhi0]);
1539 Double_t cst = TMath::Cos(fCurT[ip]+fParams[kPhi0]);
1540 Double_t snt = TMath::Sin(fCurT[ip]+fParams[kPhi0]);
1541 //
1542 int offs = kD0; // dXYZ/dD0
1543 dXYZdGlo[offs + kX] = cs0;
1544 dXYZdGlo[offs + kY] = sn0;
1545 dXYZdGlo[offs + kZ] = 0;
1546 //
1547 offs = kPhi0*3; // dXYZ/dPhi0
1d06ac63 1548 dXYZdGlo[offs + kX] = -rrho*sn0 + fParams[kR0]*snt;
1549 dXYZdGlo[offs + kY] = rrho*cs0 - fParams[kR0]*cst;
6be22b3f 1550 dXYZdGlo[offs + kZ] = 0;
1551 //
1552 offs = kR0*3; // dXYZ/dR0
ef24eb3b 1553 if (TMath::Abs(fParams[kR0])<0.9*fMaxRforHelix) {
1554 dXYZdGlo[offs + kX] = cs0 - cst;
1555 dXYZdGlo[offs + kY] = sn0 - snt;
1556 dXYZdGlo[offs + kZ] = -fParams[kDip]*fCurT[ip];
1557 }
1558 else {
1559 dXYZdGlo[offs + kX] = dXYZdGlo[offs + kY] = dXYZdGlo[offs + kZ] = 0;
1560 fParSol->AddConstraint(kR0,0,1.e2);
1561 }
6be22b3f 1562 //
1563 offs = kDZ*3; // dXYZ/dDZ
1564 dXYZdGlo[offs + kX] = 0;
1565 dXYZdGlo[offs + kY] = 0;
1566 dXYZdGlo[offs + kZ] = 1.;
1567 //
1568 offs = kDip*3; // dXYZ/dDip
1569 dXYZdGlo[offs + kX] = 0;
1570 dXYZdGlo[offs + kY] = 0;
1571 dXYZdGlo[offs + kZ] = -fParams[kR0]*fCurT[ip];
1572 //
ef24eb3b 1573 // /*
6be22b3f 1574 dXYZdLoc[kX] = fParams[kR0]*snt;
1575 dXYZdLoc[kY] = -fParams[kR0]*cst;
1576 dXYZdLoc[kZ] = -fParams[kR0]*fParams[kDip];
ef24eb3b 1577 // */
1578 // dXYZdLoc[0] = dXYZdLoc[1] = dXYZdLoc[2] = 0;
6be22b3f 1579 //
1580 fParSol->AddEquation(dXYZdGlo,dXYZdLoc,xyzRes,GetCovI(ip));
1581 }
1582 //
1583 if (extPT>0) { // add constraint on pt
1584 if (extPTerr<fgkAlmostZero) extPTerr = 1e-6*extPT;
1585 Double_t cf = fBz*GetCharge()*fgkCQConv;
1586 Double_t err2i = extPTerr/cf;
1587 err2i = 1./err2i/err2i;
1588 // printf("Constrain R to %+e\n",extPT/cf);
1589 fParSol->AddConstraint(kR0,-extPT/cf+fParams[kR0],err2i);
1590 }
1591 if (!fParSol->Solve()) { AliError("Failed to fit helix"); return -1; }
1592 Double_t *deltaG = fParSol->GetGlobals();
ef0526f3 1593 // Double_t *deltaT = fParSol->GetLocals();
6be22b3f 1594 for (int ipar=5;ipar--;) fParams[ipar] -= deltaG[ipar];
ef24eb3b 1595 //
1596 if (TMath::Abs(fParams[kR0])>0.9*fMaxRforHelix) fParams[kR0] = TMath::Sign(0.9*fMaxRforHelix,fParams[kR0]);
1597 //
1d06ac63 1598 for (int ip=fPntFirst;ip<=fPntLast;ip++) {
1599 fCurT[ip] = CalcParPCA(ip);
1600 // printf("iter%d : delta%2d : expl: %+e | %+e | R=%+e p0=%+e\n",iter,ip,deltaT[ip-fPntFirst], fCurT[ip],fParams[kR0],fParams[kPhi0]);
1601 // fCurT[ip] -= deltaT[ip-fPntFirst];
1602 }
6be22b3f 1603 iter++;
1604 //
1605 fChi2NDF = CalcChi2NDF();
ef24eb3b 1606 // printf("%2d | %+.2e %+.2e %+.2e %+.2e %+.2e | chi2: %+.4e %+.4e\n",iter,deltaG[0],deltaG[1],deltaG[2],deltaG[3],deltaG[4],fChi2NDF,fChi2NDF-chiprev);
1607 // printf("->> %+.2e %+.2e %+.2e %+.2e %+.2e | Chi2: %+.6e %+.6e\n",fParams[0],fParams[1],fParams[2],fParams[3],fParams[4],fChi2NDF,fChi2NDF-chiprev);
6be22b3f 1608 double difchi2 = chiprev - fChi2NDF;
1609 if ( difchi2<fEps && TMath::Abs(difchi2)<1e-4) {converged = kTRUE; break;}
1610 // if (errT*TMath::Abs(fParams[kR0])<kMaxTEffect && errP<fEps) {converged = kTRUE; break;}
1611 }
1612 //
1613 if (!converged) {
1614 AliDebug(2,Form("Max number of %d iteration reached, Current chi2:%.3e, chi2 change %+.3e",iter,
1615 fChi2NDF,chiprev-fChi2NDF));
1616 for (int ip=fPntFirst;ip<=fPntLast;ip++)
24391cd5 1617 AliDebug(2,Form("P%2d| %+.3e %+.3e %+.3e\n",ip,fkPoints->GetX()[ip],fkPoints->GetY()[ip],fkPoints->GetZ()[ip]));
6be22b3f 1618
1619 }
1620 fIter = iter;
1621 SetCharge( fParams[kR0]>0 ? (fBz<0?-1:1):(fBz>0?-1:1) );
1622 SetFitDone();
1623 // printf("F1>> %+.7e %+.7e %+.7e %+.7e %.7e\n",fParams[0],fParams[1],fParams[2],fParams[3],fParams[4]);
1624 //
1625 return fChi2NDF;
1626}
1627
1628//____________________________________________________
1629Double_t AliITSTPArrayFit::FitLine()
1630{
1631 // fit by helix accounting for the errors of all coordinates (and energy loss if requested)
1632 //
1633 double chiprev=1e99;
1634 // const Double_t kMaxTEffect = 1.e-6;
1635 Double_t dXYZdGlo[3*4],dXYZdLoc[3],xyzRes[3];
1636 SetFitDone(kFALSE);
1637 fChi2NDF = -1;
1638 //
1639 if (fParAxis<0) SetParAxis(ChoseParAxis());
1640 //
24391cd5 1641 const float *xyzp[3]={fkPoints->GetX(),fkPoints->GetY(),fkPoints->GetZ()};
6be22b3f 1642 if (!IsCovInv()) InvertPointsCovMat();
1643 if (!FitLineCrude()) return -1; // get initial estimate, ignoring the errors
1644 //
1645 if (!fParSol) fParSol = new AliParamSolver(5);
1646 fParSol->SetNGlobal(4);
1647 // initial set of parameters
1648 for (int ip=fPntFirst;ip<=fPntLast;ip++) fCurT[ip] = xyzp[fParAxis][ip]; // use measured param-coordinate
1649 //
1650 int iter = 0;
1651 Bool_t converged = kFALSE;
1652 fChi2NDF = 1e99;
1653 while(iter<fMaxIter) {
1654 chiprev = fChi2NDF;
1655 fParSol->Clear();
1656 for (int ip=fPntFirst;ip<=fPntLast;ip++) {
1657 //
1658 int offs;
1659 GetResiduals(xyzRes, ip); // current residuals at point ip
1660 //
1661 offs = kA0*3; // dXYZ/dA0
1662 dXYZdGlo[offs + fkAxID[kX]] = 1;
1663 dXYZdGlo[offs + fkAxID[kY]] = 0;
1664 dXYZdGlo[offs + fParAxis ] = 0;
1665 //
1666 offs = kB0*3; // dXYZ/dB0
1667 dXYZdGlo[offs + fkAxID[kX]] = fCurT[ip];
1668 dXYZdGlo[offs + fkAxID[kY]] = 0;
1669 dXYZdGlo[offs + fParAxis ] = 0;
1670 //
1671 offs = kA1*3; // dXYZ/dA1
1672 dXYZdGlo[offs + fkAxID[kX]] = 0;
1673 dXYZdGlo[offs + fkAxID[kY]] = 1;
1674 dXYZdGlo[offs + fParAxis ] = 0;
1675 //
1676 offs = kB1*3; // dXYZ/dB1
1677 dXYZdGlo[offs + fkAxID[kX]] = 0;
1678 dXYZdGlo[offs + fkAxID[kY]] = fCurT[ip];
1679 dXYZdGlo[offs + fParAxis ] = 0;
1680 //
1681 dXYZdLoc[ fkAxID[kX] ] = fParams[kB0]; // dX/dt
1682 dXYZdLoc[ fkAxID[kY] ] = fParams[kB1]; // dY/dt
1683 dXYZdLoc[ fParAxis ] = 1;
1684 //
1685 fParSol->AddEquation(dXYZdGlo,dXYZdLoc,xyzRes,GetCovI(ip));
1686 }
1687 //
1688 if (!fParSol->Solve()) { AliError("Failed to fit line"); return -1; }
1689 Double_t *deltaG = fParSol->GetGlobals();
1690 Double_t *deltaT = fParSol->GetLocals();
1691 for (int ipar=4;ipar--;) fParams[ipar] -= deltaG[ipar];
1692 for (int ip=fPntFirst;ip<=fPntLast;ip++) fCurT[ip] -= deltaT[ip-fPntFirst];
1693 iter++;
1694 fChi2NDF = CalcChi2NDF();
1695 // printf("%d %+e %+e | %+.2e %+.2e %+.2e %+.2e | chi2: %+.4e %+.4e\n",iter,errP,errT, deltaG[0],deltaG[1],deltaG[2],deltaG[3],fChi2NDF,fChi2NDF-chiprev);
1696 // printf("->> %+.2e %+.2e %+.2e %+.2e %+.2e | Chi2: %+.6e %+.6e\n",fParams[0],fParams[1],fParams[2],fParams[3],fParams[4],fChi2NDF,fChi2NDF-chiprev);
1697 double difchi2 = chiprev - fChi2NDF;
1698 if ( difchi2<fEps && TMath::Abs(difchi2)<1e-4) {converged = kTRUE; break;}
1699 chiprev = fChi2NDF;
1700 // if (errT<kMaxTEffect && errP<fEps) {converged = kTRUE; break;}
1701 }
1702 //
1703 if (!converged) {
1704 AliDebug(2,Form("Max number of %d iteration reached, Current chi2:%.3e, chi2 change %+.3e",iter,
1705 fChi2NDF,chiprev-fChi2NDF));
1706 for (int ip=fPntFirst;ip<=fPntLast;ip++)
24391cd5 1707 AliDebug(2,Form("P%2d| %+.3e %+.3e %+.3e\n",ip,fkPoints->GetX()[ip],fkPoints->GetY()[ip],fkPoints->GetZ()[ip]));
6be22b3f 1708 }
1709 fIter = iter;
1710 SetFitDone();
1711 //printf("F1>> %+.2e %+.2e %+.2e %+.2e\n",fParams[0],fParams[1],fParams[2],fParams[3]);
1712 return fChi2NDF;
1713 //
1714}
1715
1716//____________________________________________________
1717void AliITSTPArrayFit::GetNormal(Double_t *norm, const Float_t *covMat)
1718{
1719 // obtain the lab normal vector to the sensor from the covariance matrix
1720 // in such a way that when the local frame of the sensor coincides with
1721 // the lab frame, the vector {0,1,0} is obtained
1722 Double_t tgxy = TMath::Tan(0.5*TMath::ATan2(2.*covMat[kXY],covMat[kYY]-covMat[kXX]));
1723 Double_t tgyz = TMath::Tan(0.5*TMath::ATan2(2.*covMat[kYZ],covMat[kZZ]-covMat[kYY]));
1724 norm[kY] = 1./TMath::Sqrt(1 + tgxy*tgxy + tgyz*tgyz);
1725 norm[kX] = norm[kY]*tgxy;
1726 norm[kZ] = norm[kY]*tgyz;
1727 //
1728}
1729
1730//____________________________________________________
1731Double_t AliITSTPArrayFit::GetDRofELoss(Double_t t,Double_t cdip,Double_t rhoL,const Double_t *normS,
1732 Double_t &p,Double_t &e) const
1733{
1734 // Calculate energy loss of the particle at given t-param on the layer with rhoL (thickness*density) with
1735 // normal vector normS in the lab. The particle before eloss has energy "e" and momentum "p"
1736 // cdip = cosine of the dip angle = 1/sqrt(1+tgL^2)
1737 // Return the change DR of the radius due to the ELoss
1738 //
1739 // NOTE: with B>0 the negative particles propagate along increasing t-param and positive
1740 // particles - against.
1741 // t-param = 0 corresponds to the point of closest approach of the track to the beam.
1742 // Since the fitted helix parameters of the track are defined in this PCA point, when the correction
1743 // is applied upstream of the PCS, the energy must be increased (DR>0) rather than decreased (DR<0)
1744 //
1745 Double_t dirTr[3];
1746 dirTr[0] = -TMath::Sin(fParams[kPhi0]+t);
1747 dirTr[1] = TMath::Cos(fParams[kPhi0]+t);
1748 dirTr[2] = fParams[kDip];
1749 // cosine of the impact angle
1750 Double_t cosImp = cdip*TMath::Abs(dirTr[0]*normS[0]+dirTr[1]*normS[1]+dirTr[2]*normS[2]);
1751 //
1752 if (cosImp<0.3) cosImp = 0.3; //?
1753 Double_t dE = AliExternalTrackParam::BetheBlochSolid(p/fMass)*rhoL/cosImp;
1754 Double_t dP = e/p*dE;
1755 //
1756 if (t*GetSignQB() < 0) {
1757 dP = -dP;
1758 dE = -dE;
1759 }
1760 //
1761 if (p+dP<0) {
1762 AliInfo(Form("Estimated PLoss %.3f is larger than particle momentum %.3f. Skipping",dP,p));
1763 return 0;
1764 }
1765 //
1766 p += dP;
1767 e += dE;
1768 //
1769 return fCharge*dP*cdip/fBz/fgkCQConv;
1770}
1771
1772//_____________________________________________________________
1773Double_t AliITSTPArrayFit::GetLineOffset(Int_t axis) const
1774{
1775 // return intercept of the parameterization coord = intercept + slope*t for given axis
1776 if (fParAxis<0) return -1E6; // no line fit
1777 if (axis==fParAxis) return 0;
1778 if (fParAxis==kX) return fParams[axis==kY ? kA0 : kA1 ];
1779 if (fParAxis==kY) return fParams[axis==kZ ? kA0 : kA1 ];
1780 return fParams[axis==kX ? kA0 : kA1 ];
1781}
1782
1783//_____________________________________________________________
1784Double_t AliITSTPArrayFit::GetLineSlope(Int_t axis) const
1785{
1786 // return intercept of the parameterization coord = intercept + slope*t for given axis
1787 if (fParAxis<0) return -1E6; // no line fit
1788 if (axis==fParAxis) return 1.;
1789 if (fParAxis==kX) return fParams[axis==kY ? kB0 : kB1 ];
1790 if (fParAxis==kY) return fParams[axis==kZ ? kB0 : kB1 ];
1791 return fParams[axis==kX ? kB0 : kB1 ];
1792}
1793
1794//_____________________________________________________________
1795void AliITSTPArrayFit::Print(Option_t *) const
1796{
24391cd5 1797 // print results of the fit
1798 //
6be22b3f 1799 const char kCxyz[] = "XYZ";
24391cd5 1800 if (!fkPoints) return;
6be22b3f 1801 //
1802 printf("Track of %3d points in Bz=%+.1f |Fit ",fPntLast-fPntFirst+1,fBz);
1803 if ( IsFitDone() ) {
ef24eb3b 1804 if (IsHelix())
6be22b3f 1805 printf("Helix: Chi2: %5.1f | %+.2e %+.2e %+.2e %+.2e %+.2e\n",
1806 fChi2NDF,fParams[kD0],fParams[kPhi0],fParams[kR0],fParams[kDZ],fParams[kDip]);
1807 else
1808 printf("Line%c: Chi2: %5.1f | %+.2e %+.2e %+.2e %+.2e\n",
1809 kCxyz[fParAxis],fChi2NDF,fParams[kA0],fParams[kB0],fParams[kA1],fParams[kB1]);
1810 }
1811 else printf("N/A\n");
1812}
1813
1814
1815
1816
1817//____________________________________________________
1818void AliITSTPArrayFit::BuildMaterialLUT(Int_t ntri)
1819{
1820 // Fill a look-up table with mean material a la AliITSTrackerMI
1821 //
1822 if (!AliGeomManager::GetGeometry()) AliFatal("Geometry is not loaded");
1823 //
1824 // detector layer to check: dX,dZ,Ymin,Ymax
1825 const double kLayr[9][4] = {{0. ,60. , 2.80,3.00}, // beam pipe
1826 {1.28,7.07,-0.20,0.22}, // SPD1
1827 {1.28,7.07,-0.20,0.22}, // SPD2
1828 {0. ,76.0, 10.4,11.8}, // Shield1
1829 {7.02,7.53,-1.00,4.50}, // SDD1
1830 {7.02,7.53,-1.00,4.50}, // SDD2
1831 {0. ,102., 29.0,30.0}, // Shield2
1832 {7.50,4.20,-0.15,4.50}, // SSD1
1833 {7.50,4.20,-0.15,4.50}}; // SSD2
1834 //
1835 //
1836 // build <dens*L> for detectors (track hitting the sensor in normal direction)
1837 double pg1[3],pg2[3],res[7];
1838 //
1839 int sID = 0;
1840 int actLrID = 0;
1841 for (int lr=0;lr<9;lr++) {
1842 //
1843 Bool_t active = kFALSE;
1844 const double* tpars = kLayr[lr];
1845 //
1846 if (IsZero(tpars[0])) { // passive layer
1847 active = kFALSE;
1848 AliInfo(Form("Probing passive layer (total layer #%d)",lr));
1849 }
1850 else {
1851 active = kTRUE;
1852 sID += AliGeomManager::LayerSize(++actLrID);
1853 AliInfo(Form("Probing sensors of active layer #%d (total layers #%d)",actLrID,lr));
1854 }
1855 double shift = TMath::Abs(tpars[2]-tpars[3])*1E-4;
1856 double rhol = 0;
1857 for (int i=ntri;i--;) {
1858 //
1859 if (active) {
1860 int ssID = sID -1 - AliGeomManager::LayerSize(actLrID)*gRandom->Rndm();
1861 pg1[0] = pg2[0] = (gRandom->Rndm()-0.5)*tpars[0] + shift; // local X
1862 pg2[0] -= 2*shift;
1863 pg1[1] = tpars[2];
1864 pg2[1] = tpars[3];
1865 pg1[2] = pg2[2] = (gRandom->Rndm()-0.5)*tpars[1] + shift; // local Z
1866 pg2[2] -= 2*shift;
1867 AliITSgeomTGeo::LocalToGlobal(ssID,pg1,pg1);
1868 AliITSgeomTGeo::LocalToGlobal(ssID,pg2,pg2);
1869 }
1870 else {
1871 double ang = gRandom->Rndm()*TMath::Pi()*2;
1872 pg1[0] = tpars[2]*TMath::Cos(ang)+shift;
1873 pg2[0] = tpars[3]*TMath::Cos(ang)-shift;
1874 pg1[1] = tpars[2]*TMath::Sin(ang);
1875 pg2[1] = tpars[3]*TMath::Sin(ang);
1876 pg1[2] = pg2[2] = (gRandom->Rndm()-0.5)*tpars[1]+shift; // local Z
1877 pg2[2] -= 2*shift;
1878 }
1879
1880 //
1881 AliTracker::MeanMaterialBudget(pg1,pg2,res);
1882 rhol += res[0]*res[4]; // rho*L
1883 }
1884 fgRhoLITS[lr] = rhol/ntri;
1885 AliInfo(Form("Obtained <rho*L> = %e\n",fgRhoLITS[lr]));
1886 }
1887 //
1888 return;
1889}
1890
6be22b3f 1891//____________________________________________________
1892Double_t AliITSTPArrayFit::GetPCA2PlaneInfo(Double_t *xyz, Double_t *dir, Int_t axis, Double_t axval) const
1893{
1894 // calculate the PCA to plane normal ti axis and crossing it at axval
1895 // fill the position and direction cosines at this point
1896 //
1897 double xyzp[3] = {0,0,0}; // create fake point
1898 xyzp[axis] = axval;
1899 double covI[6] = {1e-4,0,0,1e-4,0,1e-4}; // fake cov.matrix loose in all directions
1900 covI[4*axis - axis*(axis+1)/2] = 1e8; // except axis
1901 //
1902 double t = GetPosition(xyz, xyzp, covI); // got pca
1903 //
1904 if (dir) GetDirCos(dir,t);
1905 return t;
1906}
1907
66214d86 1908//____________________________________________________
1909void AliITSTPArrayFit::GetT0Info(Double_t* xyz, Double_t *dir) const
1910{
1911 // get direction cosines for t = 0;
1912 GetPosition(xyz, 0);
1913 if (dir) GetDirCos(dir,0);
1914}
1915
1916//____________________________________________________
1917Bool_t AliITSTPArrayFit::CalcErrorMatrix()
1918{
1919 // compute covariance matrix in lenear approximation of residuals on parameters
1920 static AliSymMatrix cv(5);
1921 static Double_t dres[5][3];
1922 cv.Reset();
ef24eb3b 1923 int npar = IsHelix() ? 5:4;
66214d86 1924 //
1925 for (int ip=fPntFirst;ip<=fPntLast;ip++) {
1926 GetDResDParams(&dres[0][0],ip);
1927 Double_t* covI = GetCovI(ip);
1928 for (int i=npar;i--;)
1929 for (int j=i+1;j--;)
1930 cv(i,j) += dres[i][kX]*(dres[j][kX]*covI[ kXX ] + dres[j][kY]*covI[ kXY ] + dres[j][kZ]*covI[ kXZ ])
1931 + dres[i][kY]*(dres[j][kX]*covI[ kXY ] + dres[j][kY]*covI[ kYY ] + dres[j][kZ]*covI[ kYZ ])
1932 + dres[i][kZ]*(dres[j][kX]*covI[ kXZ ] + dres[j][kY]*covI[ kYZ ] + dres[j][kZ]*covI[ kZZ ]);
1933 }
1934 cv.SetSizeUsed(npar);
1935 if (cv.InvertChol()) {
ef24eb3b 1936 //cv.Print("l");
66214d86 1937 int cnt = 0;
1938 for (int i=npar;i--;) for (int j=i+1;j--;)fParamsCov[cnt++] = cv(i,j);
1939 return kTRUE;
1940 }
1941 //
1942 return kFALSE;
1943}
1d06ac63 1944
1945//____________________________________________________
1946Double_t AliITSTPArrayFit::CalcParPCA(Int_t ipnt) const
1947{
1948 // get parameter for the point with least weighted distance to the point
1949 const double *xyz = GetPoint(ipnt);
1950 const double *covI = GetCovI(ipnt);
ef24eb3b 1951 if (IsHelix()) return GetParPCAHelix(xyz,covI);
1d06ac63 1952 else return GetParPCALine(xyz,covI);
1953}
24391cd5 1954
ef24eb3b 1955//____________________________________________________
1956Double_t AliITSTPArrayFit::GetPt() const
1957{
1958 return IsFieldON()&&IsHelix() ? TMath::Abs(fParams[kR0]*fBz*fgkCQConv) : -1;
1959}
1960
1961//____________________________________________________
1962Double_t AliITSTPArrayFit::GetP() const
1963{
1964 if (!IsFieldON()) return -1;
1965 Double_t cdip = 1./TMath::Sqrt(1.+fParams[kDip]*fParams[kDip]);
1966 return TMath::Abs(fParams[kR0]*fgkCQConv*fBz/cdip); // momentum
1967}
24391cd5 1968