]> 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
changed back to unix style
[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>
47
48ClassImp(AliITSTPArrayFit)
49
50const Int_t AliITSTPArrayFit::fgkAxisID[3][3] = {
51 {AliITSTPArrayFit::kY,AliITSTPArrayFit::kZ,AliITSTPArrayFit::kX},
52 {AliITSTPArrayFit::kZ,AliITSTPArrayFit::kX,AliITSTPArrayFit::kY},
53 {AliITSTPArrayFit::kX,AliITSTPArrayFit::kY,AliITSTPArrayFit::kZ} };
54
55const Int_t AliITSTPArrayFit::fgkAxisCID[3][6] = {
56 {AliITSTPArrayFit::kYY,AliITSTPArrayFit::kYZ,AliITSTPArrayFit::kXY,
57 AliITSTPArrayFit::kZZ,AliITSTPArrayFit::kXZ,AliITSTPArrayFit::kXX},
58 //
59 {AliITSTPArrayFit::kZZ,AliITSTPArrayFit::kXZ,AliITSTPArrayFit::kYZ,
60 AliITSTPArrayFit::kXX,AliITSTPArrayFit::kYX,AliITSTPArrayFit::kYY},
61 //
62 {AliITSTPArrayFit::kXX,AliITSTPArrayFit::kXY,AliITSTPArrayFit::kXZ,
63 AliITSTPArrayFit::kYY,AliITSTPArrayFit::kYZ,AliITSTPArrayFit::kZZ}
64};
65//
66
67const Double_t AliITSTPArrayFit::fgkAlmostZero = 1E-55;
68const Double_t AliITSTPArrayFit::fgkCQConv = 0.299792458e-3;// R = PT/Bz/fgkCQConv with GeV,kGauss,cm
69const Double_t AliITSTPArrayFit::fgkZSpanITS[AliITSTPArrayFit::kMaxLrITS] = {
70 36. ,14.1,14.1, 38., 22.2,29.7, 51. ,43.1,48.9};
71
72const Double_t AliITSTPArrayFit::fgkRLayITS[AliITSTPArrayFit::kMaxLrITS] = {
73 2.94, 3.9,7.6, 11.04, 15.0,23.9, 29.44 ,38.0,43.0};
74
75const Int_t AliITSTPArrayFit::fgkPassivLrITS[3] =
76 {AliITSTPArrayFit::kLrBeamPime,AliITSTPArrayFit::kLrShield1,AliITSTPArrayFit::kLrShield2};
77
78const Int_t AliITSTPArrayFit::fgkActiveLrITS[6] =
79 {AliITSTPArrayFit::kLrSPD1,AliITSTPArrayFit::kLrSPD2,
80 AliITSTPArrayFit::kLrSDD1,AliITSTPArrayFit::kLrSDD2,
81 AliITSTPArrayFit::kLrSSD1,AliITSTPArrayFit::kLrSSD2};
82
83Double_t AliITSTPArrayFit::fgRhoLITS[AliITSTPArrayFit::kMaxLrITS] = {
84 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};
85
86//____________________________________________________
87AliITSTPArrayFit::AliITSTPArrayFit() :
88 fPoints(0),fParSol(0),fBz(0),fCharge(0),fPntFirst(-1),
89 fPntLast(-1),fNPBooked(0),fParAxis(-1),fCovI(0),fChi2NDF(0),
90 fMaxIter(20),fIter(0),fEps(1e-6),fMass(0),fkAxID(0),fkAxCID(0),fCurT(0),
91 fFirstPosT(0),fNElsPnt(0),fElsId(0),fElsDR(0)
92{
93 // default constructor
94 for (int i=kMaxParam;i--;) fParams[i] = 0;
95 for (int i=kMaxParamSq;i--;) fParamsCov[i] = 0;
96 SetMass();
97}
98
99//____________________________________________________
100AliITSTPArrayFit::AliITSTPArrayFit(Int_t np) :
101 fPoints(0),fParSol(0),fBz(0),fCharge(0),fPntFirst(-1),
102 fPntLast(-1),fNPBooked(np),fParAxis(-1),fCovI(0),fChi2NDF(0),
103 fMaxIter(20),fIter(0),fEps(1e-6),fMass(0),fkAxID(0),fkAxCID(0),fCurT(0),
104 fFirstPosT(0),fNElsPnt(0),fElsId(0),fElsDR(0)
105{
106 // constructor with booking of np points
107 for (int i=kMaxParam;i--;) fParams[i] = 0;
108 for (int i=kMaxParamSq;i--;) fParamsCov[i] = 0;
109 InitAux();
110 SetEps();
111 SetMass();
112 SetMaxIterations();
113}
114
115//____________________________________________________
116AliITSTPArrayFit::AliITSTPArrayFit(const AliITSTPArrayFit &src) :
117 TObject(src),fPoints(src.fPoints),fParSol(0),fBz(src.fBz),
118 fCharge(src.fCharge),fPntFirst(src.fPntFirst),fPntLast(src.fPntLast),fNPBooked(src.fNPBooked),
119 fParAxis(src.fParAxis),fCovI(0),fChi2NDF(0),fMaxIter(20),fIter(0),fEps(0),fMass(src.fMass),
120 fkAxID(0),fkAxCID(0),fCurT(0),fFirstPosT(0),fNElsPnt(0),fElsId(0),fElsDR(0)
121{
122 // copy constructor
123 InitAux();
124 memcpy(fCovI,src.fCovI,fNPBooked*6*sizeof(Double_t));
125 for (int i=kMaxParam;i--;) fParams[i] = src.fParams[i];
126 for (int i=kMaxParamSq;i--;) fParamsCov[i] = src.fParamsCov[i];
127 memcpy(fCurT,src.fCurT,fNPBooked*sizeof(Double_t));
128 SetEps(src.fEps);
129 SetMaxIterations(src.fMaxIter);
130 //
131}
132
133//____________________________________________________
134AliITSTPArrayFit &AliITSTPArrayFit::operator =(const AliITSTPArrayFit& src)
135{
136 // assignment operator
137 if (this==&src) return *this;
138 ((TObject*)this)->operator=(src);
139 fPoints = src.fPoints;
140 if (!fParSol) fParSol = new AliParamSolver(*src.fParSol);
141 else *fParSol = *src.fParSol;
142 fBz = src.fBz;
143 fCharge = src.fCharge;
144 fNPBooked = src.fNPBooked;
145 fPntFirst = src.fPntFirst;
146 fPntLast = src.fPntLast;
147 InitAux();
148 memcpy(fCovI,src.fCovI,fNPBooked*6*sizeof(Double_t));
149 for (int i=kMaxParam;i--;) fParams[i] = src.fParams[i];
150 for (int i=kMaxParamSq;i--;) fParamsCov[i] = src.fParamsCov[i];
151 SetParAxis(src.fParAxis);
152 fNElsPnt = src.fNElsPnt;
153 fFirstPosT = src.fFirstPosT;
154 memcpy(fCurT ,src.fCurT ,fNPBooked*sizeof(Double_t));
155 memcpy(fElsId ,src.fElsId ,fNPBooked*sizeof(Int_t));
156 memcpy(fElsDR ,src.fElsDR ,fNPBooked*sizeof(Double_t));
157 memcpy(fCurT ,src.fCurT ,fNPBooked*sizeof(Double_t));
158 SetEps(src.fEps);
159 SetMaxIterations(src.fMaxIter);
160 //
161 return *this;
162 //
163}
164
165//____________________________________________________
166AliITSTPArrayFit::~AliITSTPArrayFit()
167{
168 // destructor
169 delete fParSol;
170 delete[] fCovI;
171 delete[] fCurT;
172 delete[] fElsId;
173 delete[] fElsDR;
174}
175
176//____________________________________________________
177void AliITSTPArrayFit::Reset()
178{
179 // reset to process new track
180 if (fParSol) fParSol->Clear();
181 fPoints=0;
182 fNElsPnt = 0;
183 fFirstPosT = 0;
184 // fBz = 0;
185 fCharge = 0;
186 fIter = 0;
187 fPntFirst=fPntLast=-1;
188 SetParAxis(-1);
189 ResetBit(kFitDoneBit|kCovInvBit);
190}
191
192//____________________________________________________
193void AliITSTPArrayFit::AttachPoints(const AliTrackPointArray* points, Int_t pfirst,Int_t plast)
194{
195 // create from piece of AliTrackPointArray
196 Reset();
197 fPoints = points;
198 int np = points->GetNPoints();
199 if (fNPBooked<np) {
200 fNPBooked = np;
201 InitAux();
202 }
203 fPntFirst = pfirst<0 ? 0 : pfirst;
204 fPntLast = plast<fPntFirst ? np-1 : plast;
205 //
206 for (int i=kMaxParam;i--;) fParams[i] = 0;
207 for (int i=kMaxParamSq;i--;) fParamsCov[i] = 0;
208 //
209 InvertPointsCovMat();
210 //
211}
212
213//____________________________________________________
214Bool_t AliITSTPArrayFit::SetFirstLast(Int_t pfirst,Int_t plast)
215{
216 // set first and last point to fit
217 const AliTrackPointArray* pnts = fPoints;
218 if (!pnts) {AliError("TrackPointArray is not attached yet"); return kFALSE;}
219 AttachPoints(pnts,pfirst,plast);
220 return kTRUE;
221 //
222}
223
224//____________________________________________________
225Bool_t AliITSTPArrayFit::InvertPointsCovMat()
226{
227 // invert the cov.matrices of the points
228 for (int i=fPntFirst;i<=fPntLast;i++) {
229 //
230 const float *cov = fPoints->GetCov() + i*6; // pointer on cov.matrix
231 //
232 Double_t t0 = cov[kYY]*cov[kZZ] - cov[kYZ]*cov[kYZ];
233 Double_t t1 = cov[kXY]*cov[kZZ] - cov[kXZ]*cov[kYZ];
234 Double_t t2 = cov[kXY]*cov[kYZ] - cov[kXZ]*cov[kYY];
235 Double_t det = cov[kXX]*t0 - cov[kXY]*t1 + cov[kXZ]*t2;
abd7ef79 236 AliDebug(2,Form("%+.4e %+.4e %+.4e -> %+.4e",t0,t1,t2,det));
6be22b3f 237 if (TMath::Abs(det)<fgkAlmostZero) {
238 AliInfo(Form("Cov.Matrix for point %d is singular",i));
239 return kFALSE;
240 }
241 //
242 Double_t *covI = GetCovI(i);
243 covI[kXX] = t0/det;
244 covI[kXY] = -t1/det;
245 covI[kXZ] = t2/det;
246 covI[kYY] = (cov[kXX]*cov[kZZ] - cov[kXZ]*cov[kXZ])/det;
247 covI[kYZ] = (cov[kXY]*cov[kXZ] - cov[kXX]*cov[kYZ])/det;
248 covI[kZZ] = (cov[kXX]*cov[kYY] - cov[kXY]*cov[kXY])/det;
249 //
250 }
251 SetCovInv();
252 return kTRUE;
253}
254
255//____________________________________________________
256void AliITSTPArrayFit::InitAux()
257{
258 // init auxiliary space
259 if (fCovI) delete[] fCovI;
260 if (fCurT) delete[] fCurT;
261 //
262 fCovI = new Double_t[6*fNPBooked];
263 fCurT = new Double_t[fNPBooked+kMaxLrITS];
264 fElsId = new Int_t[fNPBooked+kMaxLrITS];
265 fElsDR = new Double_t[fNPBooked+kMaxLrITS];
266 memset(fElsDR,0,(fNPBooked+kMaxLrITS)*sizeof(Double_t));
267 memset(fCovI,0,fNPBooked*6*sizeof(Double_t));
268 //
269}
270
271//____________________________________________________
272Bool_t AliITSTPArrayFit::FitLineCrude()
273{
274 // perform linear fit w/o accounting the errors
275 // fit is done in the parameterization
276 // x = res[0] + res[1]*z
277 // y = res[2] + res[3]*z
278 // where x,y,z are NOT the lab axes but z is the lab axis along which the track
279 // has the largest lever arm and x,y are the remaining 2 axis in
280 // the order of fgkAxisID[z][0], fgkAxisID[z][1]
281 //
282 int np = fPntLast - fPntFirst + 1;
283 if (np<2) {
284 AliError("At least 2 points are needed for straight line fit");
285 return kFALSE;
286 }
287 //
288 if (fParAxis<0) SetParAxis(ChoseParAxis());
289 Double_t sZ=0,sZZ=0,sY=0,sYZ=0,sX=0,sXZ=0,det=0;
290 //
291 const float *coord[3] = {fPoints->GetX(),fPoints->GetY(),fPoints->GetZ()};
292 const Float_t *varZ = coord[ fParAxis ];
293 const Float_t *varX = coord[ fkAxID[kX] ];
294 const Float_t *varY = coord[ fkAxID[kY] ];
295 //
296 for (int i=fPntFirst;i<=fPntLast;i++) {
297 sZ += varZ[i];
298 sZZ += varZ[i]*varZ[i];
299 //
300 sX += varX[i];
301 sXZ += varX[i]*varZ[i];
302 //
303 sY += varY[i];
304 sYZ += varY[i]*varZ[i];
305 }
306 det = sZZ*np-sZ*sZ;
307 if (TMath::Abs(det)<fgkAlmostZero) return kFALSE;
308 fParams[0] = (sX*sZZ-sZ*sXZ)/det;
309 fParams[1] = (sXZ*np-sZ*sX)/det;
310 //
311 fParams[2] = (sY*sZZ-sZ*sYZ)/det;
312 fParams[3] = (sYZ*np-sZ*sY)/det;
313 //
314 return kTRUE;
315 //
316}
317
318//____________________________________________________
319void AliITSTPArrayFit::SetParAxis(Int_t ax)
320{
321 // select the axis which will be used as a parameter for the line: longest baseline
322 if (ax>kZ) {
323 AliInfo(Form("Wrong axis choice: %d",ax));
324 fParAxis = -1;
325 }
326 fParAxis = ax;
327 if (ax>=0) {
328 fkAxID = fgkAxisID[ax];
329 fkAxCID = fgkAxisCID[ax];
330 }
331 else {
332 fkAxID = fkAxCID = 0;
333 }
334 //
335}
336
337//____________________________________________________
338Int_t AliITSTPArrayFit::ChoseParAxis() const
339{
340 // select the variable with largest base as a parameter
341 Double_t cmn[3]={1.e9,1.e9,1.e9},cmx[3]={-1.e9,-1.e9,-1.e9};
342 //
343 const float *coord[3] = {fPoints->GetX(),fPoints->GetY(),fPoints->GetZ()};
344 for (int i=fPntFirst;i<=fPntLast;i++) {
345 for (int j=3;j--;) {
346 Double_t val = coord[j][i];
347 if (cmn[j]>val) cmn[j] = val;
348 if (cmx[j]<val) cmx[j] = val;
349 }
350 }
351 //
352 int axis = kZ;
353 if (cmx[axis]-cmn[axis] < cmx[kX]-cmn[kX]) axis = kX;
354 if (cmx[axis]-cmn[axis] < cmx[kY]-cmn[kY]) axis = kY;
355 return axis;
356 //
357}
358
359//____________________________________________________
360Double_t AliITSTPArrayFit::GetPosition(Double_t *xyzPCA, const Double_t *xyz, const Double_t *covI) const
361{
362 // calculate the position of the track at PCA to xyz
363 Double_t t = GetParPCA(xyz,covI);
364 GetPosition(xyzPCA,t);
365 return t;
366}
367
368//____________________________________________________
369Double_t AliITSTPArrayFit::GetPosition(Double_t *xyzPCA, const AliTrackPoint *pntCovInv) const
370{
371 // calculate the position of the track at PCA to pntCovInv
372 // NOTE: the covariance matrix of the point must be inverted
373 Double_t covI[6],xyz[3] = {pntCovInv->GetX(),pntCovInv->GetY(),pntCovInv->GetZ()};
374 for (int i=6;i--;) covI[i] = pntCovInv->GetCov()[i];
375 Double_t t = GetParPCA(xyz,covI);
376 GetPosition(xyzPCA,t);
377 return t;
378}
379
380//____________________________________________________
381void AliITSTPArrayFit::GetResiduals(Double_t *resPCA, const AliTrackPoint *pntCovInv) const
382{
383 // calculate the residuals of the track at PCA to pntCovInv
384 // NOTE: the covariance matrix of the point must be inverted
385 GetPosition(resPCA,pntCovInv);
386 resPCA[0] -= pntCovInv->GetX();
387 resPCA[1] -= pntCovInv->GetY();
388 resPCA[2] -= pntCovInv->GetZ();
389}
390
391//____________________________________________________
392void AliITSTPArrayFit::GetResiduals(Double_t *resPCA, const Double_t *xyz, const Double_t *covI) const
393{
394 // calculate the residuals of the track at PCA to xyz
395 GetPosition(resPCA,xyz,covI);
396 resPCA[kX] -= xyz[kX];
397 resPCA[kY] -= xyz[kY];
398 resPCA[kZ] -= xyz[kZ];
399}
400
401//____________________________________________________
402Double_t AliITSTPArrayFit::GetParPCALine(const Double_t *xyz, const Double_t *covI) const
403{
404 // get parameter for the point with least weighted distance to the point
405 //
406 Double_t rhs,denom;
407 Double_t dx = fParams[kA0]-xyz[ fkAxID[kX] ];
408 Double_t dy = fParams[kA1]-xyz[ fkAxID[kY] ];
409 Double_t dz = -xyz[ fkAxID[kZ] ];
410 //
411 if (covI) {
412 Double_t tx = fParams[kB0]*covI[ fkAxCID[kXX] ] + fParams[kB1]*covI[ fkAxCID[kXY] ] + covI[ fkAxCID[kXZ] ];
413 Double_t ty = fParams[kB0]*covI[ fkAxCID[kXY] ] + fParams[kB1]*covI[ fkAxCID[kYY] ] + covI[ fkAxCID[kYZ] ];
414 Double_t tz = fParams[kB0]*covI[ fkAxCID[kXZ] ] + fParams[kB1]*covI[ fkAxCID[kYZ] ] + covI[ fkAxCID[kZZ] ];
415 rhs = tx*dx + ty*dy + tz*dz;
416 denom = -(fParams[kB0]*(covI[ fkAxCID[kXZ] ] + tx) + fParams[kB1]*(covI[ fkAxCID[kYZ] ] + ty) + covI[ fkAxCID[kZZ] ]);
417 }
418 else {
419 rhs = fParams[kB0]*dx + fParams[kB1]*dy + dz;
420 denom = -(fParams[kB0]*fParams[kB0] + fParams[kB1]*fParams[kB1] + 1);
421 }
422 //
423 return rhs/denom;
424 //
425}
426
427//____________________________________________________
428void AliITSTPArrayFit::GetDResDPosLine(Double_t *dXYZdP, /*const Double_t *xyz,*/ const Double_t *covI) const
429{
430 // calculate detivative of the PCA residuals vs point position and fill in user provide
431 // array in the format {dXdXp,dY/dXp,dZdXp, ... dXdZp,dYdZp,dZdZp}
432 //
433 Double_t dTdP[3];
434 GetDtDPosLine(dTdP, /*xyz,*/ covI); // derivative of the t-param over point position
435 //
436 for (int i=3;i--;) {
437 int var = fkAxID[i];
438 Double_t *curd = dXYZdP + var*3; // d/dCoord_i
439 curd[ fkAxID[kX] ] = fParams[kB0]*dTdP[var];
440 curd[ fkAxID[kY] ] = fParams[kB1]*dTdP[var];
441 curd[ fkAxID[kZ] ] = dTdP[var];
442 curd[ var ]-= 1.;
443 }
444 //
445}
446
447//____________________________________________________
448void AliITSTPArrayFit::GetDResDParamsLine(Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI) const
449{
450 // calculate detivative of the PCA residuals vs line parameters and fill in user provide
451 // array in the format {dXdP0,dYdP0,dZdP0, ... dXdPn,dYdPn,dZdPn}
452 //
453 Double_t dTdP[4];
454 Double_t t = GetDtDParamsLine(dTdP, xyz, covI); // derivative of the t-param over line params
455 //
456 Double_t *curd = dXYZdP + kA0*3; // d/dA0
457 curd[ fkAxID[kX] ] = fParams[kB0]*dTdP[kA0] + 1.;
458 curd[ fkAxID[kY] ] = fParams[kB1]*dTdP[kA0];
459 curd[ fkAxID[kZ] ] = dTdP[kA0];
460 //
461 curd = dXYZdP + kB0*3; // d/dB0
462 curd[ fkAxID[kX] ] = fParams[kB0]*dTdP[kB0] + t;
463 curd[ fkAxID[kY] ] = fParams[kB1]*dTdP[kB0];
464 curd[ fkAxID[kZ] ] = dTdP[kB0];
465 //
466 curd = dXYZdP + kA1*3; // d/dA1
467 curd[ fkAxID[kX] ] = fParams[kB0]*dTdP[kA1];
468 curd[ fkAxID[kY] ] = fParams[kB1]*dTdP[kA1] + 1.;
469 curd[ fkAxID[kZ] ] = dTdP[kA1];
470 //
471 curd = dXYZdP + kB1*3; // d/dB1
472 curd[ fkAxID[kX] ] = fParams[kB0]*dTdP[kB1];
473 curd[ fkAxID[kY] ] = fParams[kB1]*dTdP[kB1] + t;
474 curd[ fkAxID[kZ] ] = dTdP[kB1];
475 //
476}
477
478//____________________________________________________
479Double_t AliITSTPArrayFit::GetDtDParamsLine(Double_t *dtparam,const Double_t *xyz, const Double_t *covI) const
480{
481 // get t-param detivative over the parameters for the point with least weighted distance to the point
482 //
483 Double_t rhs,denom;
484 Double_t dx = fParams[kA0]-xyz[ fkAxID[kX] ];
485 Double_t dy = fParams[kA1]-xyz[ fkAxID[kY] ];
486 Double_t dz = -xyz[ fkAxID[kZ] ];
487 Double_t rhsDA0,rhsDA1,rhsDB0,rhsDB1,denDB0,denDB1;
488 //
489 if (covI) {
490 Double_t tx = fParams[kB0]*covI[ fkAxCID[kXX] ] + fParams[kB1]*covI[ fkAxCID[kXY] ] + covI[ fkAxCID[kXZ] ];
491 Double_t ty = fParams[kB0]*covI[ fkAxCID[kXY] ] + fParams[kB1]*covI[ fkAxCID[kYY] ] + covI[ fkAxCID[kYZ] ];
492 Double_t tz = fParams[kB0]*covI[ fkAxCID[kXZ] ] + fParams[kB1]*covI[ fkAxCID[kYZ] ] + covI[ fkAxCID[kZZ] ];
493 rhs = tx*dx + ty*dy + tz*dz;
494 denom = -(fParams[kB0]*(covI[ fkAxCID[kXZ] ] + tx) + fParams[kB1]*(covI[ fkAxCID[kYZ] ] + ty) + covI[ fkAxCID[kZZ] ]);
495 //
496 rhsDA0 = tx;
497 rhsDA1 = ty;
498 rhsDB0 = covI[ fkAxCID[kXX] ]*dx + covI[ fkAxCID[kXY] ]*dy + covI[ fkAxCID[kXZ] ]*dz;
499 rhsDB1 = covI[ fkAxCID[kXY] ]*dx + covI[ fkAxCID[kYY] ]*dy + covI[ fkAxCID[kYZ] ]*dz;
500 //
501 denDB0 = -(tx + tx);
502 denDB1 = -(ty + ty);
503 }
504 else {
505 rhs = fParams[kB0]*dx + fParams[kB1]*dy + dz;
506 denom = -(fParams[kB0]*fParams[kB0] + fParams[kB1]*fParams[kB1] + 1);
507 //
508 rhsDA0 = fParams[kB0];
509 rhsDB0 = dx;
510 rhsDA1 = fParams[kB1];
511 rhsDB1 = dy;
512 //
513 denDB0 = -(fParams[kB0]+fParams[kB0]);
514 denDB1 = -(fParams[kB1]+fParams[kB1]);
515 //
516 }
517 //
518 Double_t denom2 = denom*denom;
519 dtparam[kA0] = rhsDA0/denom; // denom does not depend on A0,A1
520 dtparam[kA1] = rhsDA1/denom;
521 dtparam[kB0] = rhsDB0/denom - rhs/denom2 * denDB0;
522 dtparam[kB1] = rhsDB1/denom - rhs/denom2 * denDB1;
523 //
524 return rhs/denom;
525}
526
527//____________________________________________________
528void AliITSTPArrayFit::GetDtDPosLine(Double_t *dtpos,/*const Double_t *xyz,*/ const Double_t *covI) const
529{
530 // get t-param detivative over the parameters for the point with least weighted distance to the point
531 //
532 // Double_t rhs;
533 // Double_t dx = fParams[kA0]-xyz[ fkAxID[kX] ];
534 // Double_t dy = fParams[kA1]-xyz[ fkAxID[kY] ];
535 // Double_t dz = -xyz[ fkAxID[kZ] ];
536 Double_t denom;
537 Double_t rhsDX,rhsDY,rhsDZ;
538 //
539 if (covI) {
540 Double_t tx = fParams[kB0]*covI[ fkAxCID[kXX] ] + fParams[kB1]*covI[ fkAxCID[kXY] ] + covI[ fkAxCID[kXZ] ];
541 Double_t ty = fParams[kB0]*covI[ fkAxCID[kXY] ] + fParams[kB1]*covI[ fkAxCID[kYY] ] + covI[ fkAxCID[kYZ] ];
542 Double_t tz = fParams[kB0]*covI[ fkAxCID[kXZ] ] + fParams[kB1]*covI[ fkAxCID[kYZ] ] + covI[ fkAxCID[kZZ] ];
543 // rhs = tx*dx + ty*dy + tz*dz;
544 denom = -(fParams[kB0]*(covI[ fkAxCID[kXZ] ] + tx) + fParams[kB1]*(covI[ fkAxCID[kYZ] ] + ty) + covI[ fkAxCID[kZZ] ]);
545 //
546 rhsDX = -tx;
547 rhsDY = -ty;
548 rhsDZ = -tz;
549 }
550 else {
551 // rhs = fParams[kB0]*dx + fParams[kB1]*dy + dz;
552 denom = -(fParams[kB0]*fParams[kB0] + fParams[kB1]*fParams[kB1] + 1);
553 //
554 rhsDX = -fParams[kB0];
555 rhsDY = -fParams[kB1];
556 rhsDZ = -1;
557 //
558 }
559 //
560 dtpos[ fkAxID[kX] ] = rhsDX/denom;
561 dtpos[ fkAxID[kY] ] = rhsDY/denom;
562 dtpos[ fkAxID[kZ] ] = rhsDZ/denom;
563 //
564 // return rhs/denom;
565}
566
567//____________________________________________________
568void AliITSTPArrayFit::GetDResDParamsLine(Double_t *dXYZdP, Int_t ipnt) const
569{
570 // calculate detivative of the PCA residuals vs line parameters and fill in user provide
571 // array in the format {dXdP0,dYdP0,dZdP0, ... dXdPn,dYdPn,dZdPn}
572 //
573 if (ipnt<fPntFirst || ipnt>fPntLast) {
574 AliError(Form("Attempt to access the point %d not in the fitted points [%d:%d]",ipnt,fPntFirst,fPntLast));
575 return;
576 }
577 GetDResDParamsLine(dXYZdP, GetPoint(ipnt) , IsCovIgnored() ? 0 : GetCovI(ipnt));
578}
579
580//____________________________________________________
581void AliITSTPArrayFit::GetDResDPosLine(Double_t *dXYZdP, Int_t ipnt) const
582{
583 // calculate detivative of the PCA residuals vs point position and fill in user provide
584 // array in the format {dXdXp,dY/dXp,dZdXp, ... dXdZp,dYdZp,dZdZp}
585 //
586 if (ipnt<fPntFirst || ipnt>fPntLast) {
587 AliError(Form("Attempt to access the point %d not in the fitted points [%d:%d]",ipnt,fPntFirst,fPntLast));
588 return;
589 }
590 GetDResDPosLine(dXYZdP,IsCovIgnored() ? 0 : GetCovI(ipnt));
591}
592
593//____________________________________________________
594void AliITSTPArrayFit::GetDResDParams(Double_t *dXYZdP, Int_t ipnt)
595{
596 // calculate detivative of the PCA residuals vs track parameters and fill in user provide
597 // array in the format {dXdP0,dYdP0,dZdP0, ... dXdPn,dYdPn,dZdPn}
598 //
599 if (ipnt<fPntFirst || ipnt>fPntLast) {
600 AliError(Form("Attempt to access the point %d not in the fitted points [%d:%d]",ipnt,fPntFirst,fPntLast));
601 return;
602 }
603 GetDResDParams(dXYZdP, GetPoint(ipnt) , IsCovIgnored() ? 0 : GetCovI(ipnt));
604}
605
606//____________________________________________________
607void AliITSTPArrayFit::GetDResDPos(Double_t *dXYZdP, Int_t ipnt)
608{
609 // calculate detivative of the PCA residuals vs point position and fill in user provide
610 // array in the format {dXdXp,dY/dXp,dZdXp, ... dXdZp,dYdZp,dZdZp}
611 //
612 if (ipnt<fPntFirst || ipnt>fPntLast) {
613 AliError(Form("Attempt to access the point %d not in the fitted points [%d:%d]",ipnt,fPntFirst,fPntLast));
614 return;
615 }
616 GetDResDPos(dXYZdP, GetPoint(ipnt), IsCovIgnored() ? 0 : GetCovI(ipnt));
617}
618
619//____________________________________________________
620void AliITSTPArrayFit::GetDResDParams(Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI)
621{
622 // get residual detivatives over the track parameters for the point with least weighted distance to the point
623 //
624 if (!IsFieldON()) { // for the straight line calculate analytically
625 GetDResDParamsLine(dXYZdP, xyz, covI);
626 return;
627 }
628 //
629 // calculate derivative numerically
630 const Double_t delta = 0.01;
631 Double_t xyzVar[4][3];
632 //
633 for (int ipar = 5;ipar--;) {
634 double sav = fParams[ipar];
635 fParams[ipar] -= delta;
636 GetPosition(xyzVar[0],xyz,covI);
637 fParams[ipar] += delta/2;
638 GetPosition(xyzVar[1],xyz,covI);
639 fParams[ipar] += delta;
640 GetPosition(xyzVar[2],xyz,covI);
641 fParams[ipar] += delta/2;
642 GetPosition(xyzVar[3],xyz,covI);
643 fParams[ipar] = sav; // restore
644 //
645 double *curd = dXYZdP + 3*ipar;
646 for (int i=3;i--;) curd[i] = (8.*(xyzVar[2][i]-xyzVar[1][i]) - (xyzVar[3][i]-xyzVar[0][i]))/6./delta;
647 }
648 //
649}
650
651
652//____________________________________________________
653void AliITSTPArrayFit::GetDResDPos(Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI)
654{
655 // get residuals detivative over the point position for the point with least weighted distance to the point
656 //
657
658 if (!IsFieldON()) { // for the straight line calculate analytically
659 GetDResDPosLine(dXYZdP, /*xyz,*/ covI);
660 return;
661 }
662 //
663 // calculate derivative numerically
664 const Double_t delta = 0.005;
665 Double_t xyzVar[4][3];
666 Double_t xyzv[3] = {xyz[0],xyz[1],xyz[2]};
667 //
668 for (int ipar = 3;ipar--;) {
669 double sav = xyzv[ipar];
670 xyzv[ipar] -= delta;
671 GetPosition(xyzVar[0],xyzv,covI);
672 xyzv[ipar] += delta/2;
673 GetPosition(xyzVar[1],xyzv,covI);
674 xyzv[ipar] += delta;
675 GetPosition(xyzVar[2],xyzv,covI);
676 xyzv[ipar] += delta/2;
677 GetPosition(xyzVar[3],xyzv,covI);
678 xyzv[ipar] = sav; // restore
679 //
680 double *curd = dXYZdP + 3*ipar;
681 for (int i=3;i--;) curd[i] = (8.*(xyzVar[2][i]-xyzVar[1][i]) - (xyzVar[3][i]-xyzVar[0][i]))/6./delta;
682 curd[ipar] -= 1.;
683 }
684 //
685}
686
687//________________________________________________________________________________________________________
688Double_t AliITSTPArrayFit::GetParPCAHelix(const Double_t* xyz, const Double_t* covI) const
689{
690 // find track parameter t (eq.2) corresponding to point of closest approach to xyz
691 //
692 Double_t phi = GetParPCACircle(xyz[kX],xyz[kY]);
693 Double_t cs = TMath::Cos(fParams[kPhi0]);
694 Double_t sn = TMath::Sin(fParams[kPhi0]);
695 Double_t xc = (fParams[kD0]+fParams[kR0])*cs;
696 Double_t yc = (fParams[kD0]+fParams[kR0])*sn;
697 Double_t dchi2,ddchi2;
698 //
699 Double_t dzD = -fParams[kR0]*fParams[kDip];
700 Double_t dphi = 0;
701 //
702 int it=0;
703 do {
704 cs = TMath::Cos(phi + fParams[kPhi0]);
705 sn = TMath::Sin(phi + fParams[kPhi0]);
706 //
707 Double_t dxD = fParams[kR0]*sn;
708 Double_t dyD = -fParams[kR0]*cs;
709 Double_t dxDD = -dyD;
710 Double_t dyDD = dxD;
711 //
712 Double_t dx = xc - fParams[kR0]*cs - xyz[kX];
713 Double_t dy = yc - fParams[kR0]*sn - xyz[kY];
714 Double_t dz = fParams[kDZ] + dzD*phi- xyz[kZ];
715 //
716 if (covI) {
717 Double_t tx = dx*covI[kXX] + dy*covI[kXY] + dz*covI[kXZ];
718 Double_t ty = dx*covI[kXY] + dy*covI[kYY] + dz*covI[kYZ];
719 Double_t tz = dx*covI[kXZ] + dy*covI[kYZ] + dz*covI[kZZ];
720 //
721 Double_t ttx = dxD*covI[kXX] + dyD*covI[kXY] + dzD*covI[kXZ];
722 Double_t tty = dxD*covI[kXY] + dyD*covI[kYY] + dzD*covI[kYZ];
723 Double_t ttz = dxD*covI[kXZ] + dyD*covI[kYZ] + dzD*covI[kZZ];
724 //
725 // chi2 = dx*tx + dy*ty + dz*tz;
726 dchi2 = dxD*tx + dyD*ty + dzD*tz;
727 ddchi2 = dxDD*tx + dyDD*ty + dxD *ttx + dyD *tty + dzD *ttz;
728 //
729 }
730 else {
731 // chi2 = dx*dx + dy*dy + dz*dz;
732 dchi2 = dxD*dx + dyD*dy + dzD*dz;
733 ddchi2 = dxDD*dx + dyDD*dy + + dxD*dxD + dyD*dyD + dzD*dzD;
734 }
735 //
736 if (TMath::Abs(ddchi2)<fgkAlmostZero || TMath::Abs(dphi=dchi2/ddchi2)<fEps) break;
737 phi -= dphi;
738 } while(++it<fMaxIter);
739 //
740 return phi;
741}
742
743//________________________________________________________________________________________________________
744Double_t AliITSTPArrayFit::GetParPCACircle(Double_t x,Double_t y) const
745{
746 // find track parameter t (eq.2) corresponding to point on the circle with closest approach to x,y
747 //
748 Double_t r = fParams[kD0]+fParams[kR0];
749 Double_t t = TMath::ATan2( r*TMath::Sin(fParams[kPhi0])-y, r*TMath::Cos(fParams[kPhi0])-x ) - fParams[kPhi0];
750 if (fParams[kR0] < 0) t += TMath::Pi();
751 if (t > TMath::Pi()) t -= TMath::Pi()*2;
752 if (t <-TMath::Pi()) t += TMath::Pi()*2;
753 return t;
754}
755
756//________________________________________________________________________________________________________
757Double_t AliITSTPArrayFit::GetHelixParAtR(Double_t r) const
758{
759 // find helix parameter t (eq.2) corresponding to point on the circle of radius t
760 //
761 double gam = 1. - (r-fParams[kD0])*(r+fParams[kD0])/fParams[kR0]/(fParams[kD0]+fParams[kR0])/2.;
762 return (TMath::Abs(gam)>1) ? -1e9 : TMath::ACos(gam);
763}
764
765//________________________________________________________________________________________________________
766Double_t AliITSTPArrayFit::CalcChi2NDF() const
767{
768 // calculate fit chi2/ndf
769 Double_t chi2 = 0;
770 Double_t dr[3]; // residuals
771 //if (!IsFitDone()) return -1;
772 for (int ipnt=fPntFirst;ipnt<=fPntLast;ipnt++) {
773 GetResiduals(dr,ipnt);
774 Double_t* covI = GetCovI(ipnt);
775 chi2 += dr[kX]*(dr[kX]*covI[ kXX ]+dr[kY]*covI[ kXY ]+dr[kZ]*covI[ kXZ ])
776 + dr[kY]*(dr[kX]*covI[ kXY ]+dr[kY]*covI[ kYY ]+dr[kZ]*covI[ kYZ ])
777 + dr[kZ]*(dr[kX]*covI[ kXZ ]+dr[kY]*covI[ kYZ ]+dr[kZ]*covI[ kZZ ]);
778 }
779 int ndf = (fPntLast-fPntFirst+1)*3 - GetNParams();
780 chi2 /= ndf;
781 return chi2;
782}
783
784//________________________________________________________________________________________________________
785void AliITSTPArrayFit::GetResiduals(Double_t *res,Int_t ipnt) const
786{
787 // calculate residuals at point
788 if (ipnt<fPntFirst || ipnt>fPntLast) {
789 AliError(Form("Attempt to access the point %d not in the fitted points [%d:%d]",ipnt,fPntFirst,fPntLast));
790 return;
791 }
792 GetPosition(res,fCurT[ipnt]);
793 res[kX] -= fPoints->GetX()[ipnt];
794 res[kY] -= fPoints->GetY()[ipnt];
795 res[kZ] -= fPoints->GetZ()[ipnt];
796}
797
798//________________________________________________________________________________________________________
799void AliITSTPArrayFit::GetPosition(Double_t *xyz, Double_t t) const
800{
801 // calculate track position for parameter value t
802 if (IsFieldON()) {
803 //
804 Double_t rrho = fParams[kD0]+fParams[kR0];
805 Double_t xc = rrho*TMath::Cos(fParams[kPhi0]);
806 Double_t yc = rrho*TMath::Sin(fParams[kPhi0]);
807 Double_t r = fParams[kR0];
808 Double_t ze = 0;
809 //
810 if (IsELossON()) {
811 if (t>0) {
812 for (int i=fFirstPosT;i<fNElsPnt;i++) { // along the track direction
813 int indE = fElsId[i];
814 if ( t<fCurT[indE] ) break; // does not reach this layer on its way to t
815 xc += fElsDR[indE] * TMath::Cos(fParams[kPhi0] + fCurT[indE]);
816 yc += fElsDR[indE] * TMath::Sin(fParams[kPhi0] + fCurT[indE]);
817 ze += fElsDR[indE] * fCurT[indE];
818 r += fElsDR[indE];
819 //printf("ELoss@ %+.2e r:%+.3e got %+.3e\n",fCurT[indE],r,fElsDR[indE]);
820 }
821 } else {
822 for (int i=fFirstPosT;i--;) { // against the track direction
823 int indE = fElsId[i];
824 if ( t>=fCurT[indE] ) break; // does not reach this layer on its way to t
825 xc += fElsDR[indE] * TMath::Cos(fParams[kPhi0] + fCurT[indE]);
826 yc += fElsDR[indE] * TMath::Sin(fParams[kPhi0] + fCurT[indE]);
827 ze += fElsDR[indE] * fCurT[indE];
828 r += fElsDR[indE];
829 //printf("ELoss@ %+.2e r:%+.3e got %+.3e\n",fCurT[indE],r,fElsDR[indE]);
830 }
831 }
832 }
833 //
834 xyz[kZ] = fParams[kDZ] - fParams[kDip]*(t*r - ze);
835 //
836 t += fParams[kPhi0];
837 xyz[kX] = xc - r*TMath::Cos(t);
838 xyz[kY] = yc - r*TMath::Sin(t);
839 // 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());
840 }
841 else {
842 xyz[ fkAxID[kX] ] = fParams[kA0] + fParams[kB0]*t;
843 xyz[ fkAxID[kY] ] = fParams[kA1] + fParams[kB1]*t;
844 xyz[ fParAxis ] = t;
845 }
846}
847
848//________________________________________________________________________________________________________
849void AliITSTPArrayFit::GetDirCos(Double_t *dircos, Double_t t) const
850{
851 // calculate track direction cosines for parameter value t
852 if (IsFieldON()) {
853 dircos[kZ] = -fParams[kDip];
854 t += fParams[kPhi0];
855 dircos[kX] = TMath::Sin(t);
856 dircos[kY] =-TMath::Cos(t);
857 double gam = TMath::Sign(1/TMath::Sqrt(dircos[kZ]*dircos[kZ]+dircos[kY]*dircos[kY]+dircos[kX]*dircos[kX]),fParams[kR0]);
858 for (int i=3;i--;) dircos[i] *= gam;
859 }
860 else {
861 double gam = 1/TMath::Sqrt( fParams[kB0]*fParams[kB0] + fParams[kB1]*fParams[kB1] + 1.);
862 dircos[ fkAxID[kX] ] = fParams[kB0]*gam;
863 dircos[ fkAxID[kY] ] = fParams[kB1]*gam;
864 dircos[ fParAxis ] = gam;
865 }
866}
867
868//________________________________________________________________________________________________________
869Double_t AliITSTPArrayFit::GetMachinePrec()
870{
871 // estimate machine precision
872 Double_t eps=1.0,a;
873 do { a = 1. + (eps=eps/2.0); } while(a>1.);
874 return TMath::Abs(2.*eps);
875}
876
877//________________________________________________________________________________________________________
878Bool_t AliITSTPArrayFit::FitHelixCrude(Int_t extQ)
879{
880 // crude estimate of helix parameters, w/o errors and Eloss.
881 // 1st fit the circle (R,xc,yc) by minimizing
882 // chi2 = sum{ (bx*xi + by*yi + xi^2+yi^2 + rho)^2 } vs bx,by,rho
883 // with bx = -2*xc, by = -2*yc , rho = xc^2+yc^2 - R2
884 //
885 // if charge is not imposed (extQ==0) then it will be determined from the collision type
886 //
887 Bool_t eloss = IsELossON();
888 //
889 int np = fPntLast - fPntFirst + 1;
890 if (np<2) { AliError("At least 3 points are needed for helix fit"); return kFALSE; }
891 //
892 const float *x=fPoints->GetX(),*y=fPoints->GetY(),*z=fPoints->GetZ(),*cov=fPoints->GetCov();
893 //
894 // linear circle fit --------------------------------------------------- >>>
895 Double_t sxx=0,sxy=0,syy=0,sx=0,sy=0,rhs0=0,rhs1=0,rhs2=0,minR=1E9;
896 int minRId = 0;
897 for (int i=fPntFirst;i<=fPntLast;i++) {
898 Double_t xx = x[i]*x[i];
899 Double_t yy = y[i]*y[i];
900 Double_t xy = x[i]*y[i];
901 Double_t xxyy = xx + yy;
902 //
903 sxx += xx;
904 sxy += xy;
905 syy += yy;
906 sx += x[i];
907 sy += y[i];
908 //
909 rhs0 -= xxyy*x[i];
910 rhs1 -= xxyy*y[i];
911 rhs2 -= xxyy;
912 //
913 // remember Id of the point closest to origin, to determine the charge
914 if (xxyy<minR) { minR = xxyy; minRId = i; }
915 //
916 if (eloss) { // find layer id
917 int lrid,volid = fPoints->GetVolumeID()[i];
918 if (volid>0) lrid = fgkActiveLrITS[AliGeomManager::VolUIDToLayer(fPoints->GetVolumeID()[i])-1];
919 else { // missing layer info, find from radius
920 double r = TMath::Sqrt(xxyy);
921 for (lrid=kMaxLrITS;lrid--;) if ( IsZero(r-fgkRLayITS[ lrid ],1.) ) break;
922 }
923 fElsDR[i] = (lrid>=0 && lrid<kMaxLrITS) ? fgRhoLITS[ lrid ] : 0; // eloss for normal track
924 }
925 //
926 }
927 //
928 Double_t mn00 = syy*np-sy*sy;
929 Double_t mn01 = sxy*np-sy*sx;
930 Double_t mn02 = sxy*sy-syy*sx;
931 Double_t det = sxx*mn00 - sxy*mn01 + sx*mn02;
932 if (TMath::Abs(det)<fgkAlmostZero) return kFALSE;
933 //
934 Double_t mn11 = sxx*np-sx*sx;
935 Double_t mn12 = sxx*sy-sxy*sx;
936 Double_t mn22 = sxx*syy-sxy*sxy;
937 //
938 Double_t mi00 = mn00/det;
939 Double_t mi01 = -mn01/det;
940 Double_t mi02 = mn02/det;
941 Double_t mi11 = mn11/det;
942 Double_t mi12 = -mn12/det;
943 Double_t mi22 = mn22/det;
944 //
945 Double_t xc = -(rhs0*mi00 + rhs1*mi01 + rhs2*mi02)/2;
946 Double_t yc = -(rhs0*mi01 + rhs1*mi11 + rhs2*mi12)/2;
947 Double_t rho2 = (rhs0*mi02 + rhs1*mi12 + rhs2*mi22);
948 //
949 Double_t dcen = xc*xc + yc*yc;
950 Double_t rad = dcen - rho2;
951 rad = (rad>fgkAlmostZero) ? (TMath::Sqrt(rad)):fgkAlmostZero;
952 //
953 // printf("Rad: %+e xc: %+e yc: %+e\n",rad,xc,yc);
954 // linear circle fit --------------------------------------------------- <<<
955 //
956 // decide sign(Q*B) and fill cicrle parameters ------------------------- >>>
957 int sqb;
958 if (extQ) {
959 SetCharge(extQ);
960 sqb = fBz<0 ? -GetCharge():GetCharge();
961 }
962 else {
963 // determine the charge from the collision type and field sign
964 // the negative Q*B will have positive Vc x V0 product Z component
965 // with Vc={-xc,-yc} : vector from circle center to the origin
966 // and V0 - track direction vector (take {0,-1,1} for cosmics)
967 // If Bz is not provided, assume positive Bz
968 sqb = ( IsTypeCosmics() ? xc:(yc*x[minRId]-xc*y[minRId]) ) > 0 ? -1:1;
969 SetCharge( fBz<0 ? -sqb : sqb);
970 }
971 //
972 dcen = TMath::Sqrt(dcen);
973 fParams[kD0] = dcen-rad;
974 Double_t phi = TMath::ATan2(yc,xc);
975 if (sqb<0) phi += TMath::Pi();
976 if (phi > TMath::Pi()) phi -= 2.*TMath::Pi();
977 else if (phi <-TMath::Pi()) phi += 2.*TMath::Pi();
978 fParams[kPhi0] = phi;
979 fParams[kR0] = sqb<0 ? -rad:rad;
980 //
981 // decide sign(Q*B) and fill cicrle parameters ------------------------- <<<
982 //
983 // find z-offset and dip + the parameter t of closest approach to hits - >>>
984 //
985 UInt_t hitLrPos=0; // pattern of hit layers at pos
986 UInt_t hitLrNeg=0; // and negative t's
987
988 Double_t ss=0,st=0,sz=0,stt=0,szt=0;
989 for (int i=fPntFirst;i<=fPntLast;i++) {
990 //
991 Double_t ze2 = cov[i*6 + kZZ];
992 Double_t t = TMath::ATan2(yc-y[i],xc-x[i]) - fParams[kPhi0]; // angle at measured z
993 if (fParams[kR0]<0) t += TMath::Pi();
994 if (t > TMath::Pi()) t -= TMath::Pi()*2;
995 else if (t <-TMath::Pi()) t += TMath::Pi()*2;
996 if (ze2<fgkAlmostZero) ze2 = 1E-8;
997 ze2 = 1./ze2;
998 ss += ze2;
999 st += t*ze2;
1000 stt+= t*t*ze2;
1001 sz += z[i]*ze2;
1002 szt+= z[i]*t*ze2;
1003 //
1004 fCurT[i] = t; // parameter of the closest approach to the point
1005 // printf("%d %+e %+e %+e %+e\n",i,x[i],y[i],z[i],t);
1006 if (eloss) {
1007 double r = TMath::Sqrt(x[i]*x[i]+y[i]*y[i]);
1008 int lr;
1009 for (lr=kMaxLrITS;lr--;) if ( IsZero(r-fgkRLayITS[ lr ],1.) ) break;
1010 if (lr<kMaxLrITS) {
1011 if (t>0) hitLrPos |= (1<<lr); // set bit of the layer
1012 else hitLrNeg |= (1<<lr); // set bit of the layer
1013 }
1014 }
1015 }
1016 det = ss*stt - st*st;
1017 if (TMath::Abs(det)<fgkAlmostZero) { // no Z dependence
1018 fParams[kDZ] = sz/ss;
1019 fParams[kDip] = 0;
1020 }
1021 else {
1022 fParams[kDZ] = (sz*stt-st*szt)/det;
1023 fParams[kDip] = -(ss*szt-st*sz)/det/fParams[kR0];
1024 }
1025 //
1026 // find z-offset and dip + the parameter t of closest approach to hits - <<<
1027 //
1028 // fill info needed to account for ELoss ------------------------------- >>>
1029 if (eloss) {
1030 fNElsPnt = fPntLast - fPntFirst + 1;
1031 //
1032 // to account for the energy loss in the passive volumes, calculate the relevant t-parameters
1033 double* tcur = fCurT + fPntFirst;
1034 double* ecur = fElsDR+ fPntFirst;
1035 //
1036 for (int ilp=3;ilp--;) {
1037 int id = fgkPassivLrITS[ilp];
1038 double tp = GetHelixParAtR( fgkRLayITS[ id ] );
1039 if (tp<0) continue; // does not hit this radius
1040 //
1041 tcur[fNElsPnt] = GetSignQB()>0 ? -tp : tp;
1042 ecur[fNElsPnt] = fgRhoLITS[ id ];
1043 fNElsPnt++;
1044 // printf("Passive on lr %d %+e\n",ilp,tcur[fNElsPnt-1]);
1045 //
1046 if (IsTypeCosmics() && !IsZero(tp)) { // 2 crossings for cosmics
1047 tcur[fNElsPnt] = -tcur[fNElsPnt-1];
1048 ecur[fNElsPnt] = ecur[fNElsPnt-1];
1049 fNElsPnt++;
1050 //printf("Passive* on lr %d %+e\n",ilp,-tcur[fNElsPnt-1]);
1051 }
1052 //
1053 }
1054 // check if some active layers did not miss the hit, treat them as passive
1055 for (int ilp=6;ilp--;) {
1056 int id = fgkActiveLrITS[ilp];
1057 double tp = GetHelixParAtR( fgkRLayITS[ id ] );
1058 if (tp<0) continue; // does not hit this radius
1059 //
1060 if ( (GetSignQB()>0||IsTypeCosmics()) && !(hitLrNeg & (1<<id)) ) {
1061 tcur[fNElsPnt] = -tp;
1062 ecur[fNElsPnt] = fgRhoLITS[ id ];
1063 fNElsPnt++;
1064 //printf("Missed on lr %d %+e\n",ilp,-tp);
1065 }
1066 //
1067 if ( (GetSignQB()<0||IsTypeCosmics()) && !(hitLrPos & (1<<id)) ) {
1068 tcur[fNElsPnt] = tp;
1069 ecur[fNElsPnt] = fgRhoLITS[ id ];
1070 fNElsPnt++;
1071 //printf("Missed* on lr %d %e\n",ilp,tp);
1072 }
1073 }
1074 //
1075 TMath::Sort(fNElsPnt,fCurT+fPntFirst,fElsId,kFALSE); // index e-loss points in increasing order
1076 // find the position of smallest positive t-param
1077 for (fFirstPosT=0;fFirstPosT<fNElsPnt;fFirstPosT++) if (fCurT[ fElsId[ fFirstPosT ] ]>0) break;
1078 //
1079 Double_t cdip = 1./TMath::Sqrt(1.+fParams[kDip]*fParams[kDip]);
1080 Double_t ptot = TMath::Abs(fParams[kR0]*fgkCQConv*fBz/cdip); // momentum and energy
1081 Double_t etot = TMath::Sqrt(ptot*ptot + fMass*fMass); // in the point of closest approach to beam
1082 Double_t normS[3];
1083 //
1084 // Positive t-params: along the track direction for negative track, against for positive
1085 Double_t pcur = ptot, ecurr = etot;
1086 for (int ip=fFirstPosT;ip<fNElsPnt;ip++) {
1087 int tID = fElsId[ip];
1088 Double_t t = fCurT[ tID ];
1089 //
1090 if (tID>fPntLast) { // this is not a hit layer but passive layer
1091 double php = TMath::ATan2(yc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t),
1092 xc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t));
1093 normS[0] = -TMath::Cos(php); // normal to the cylinder at intersection point
1094 normS[1] = -TMath::Sin(php);
1095 normS[2] = 0;
1096 }
1097 else GetNormal(normS,fPoints->GetCov()+tID*6); // vector normal to hit module
1098 fElsDR[tID] = GetDRofELoss(t,cdip,fElsDR[tID],normS,ptot,etot);
1099 }
1100 //
1101 // negaive t-params: against the track direction for negative track, along for positive
1102 pcur = ptot;
1103 ecurr = etot;
1104 for (int ip=fFirstPosT;ip--;) {
1105 int tID = fElsId[ip];
1106 Double_t t = fCurT[ tID ];
1107 //
1108 if (tID>fPntLast) { // this is not a hit layer but passive layer
1109 double php = TMath::ATan2(yc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t),
1110 xc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t));
1111 normS[0] = -TMath::Cos(php); // normal to the cylinder at intersection point
1112 normS[1] = -TMath::Sin(php);
1113 normS[2] = 0;
1114 }
1115 else GetNormal(normS,fPoints->GetCov()+tID*6); // vector normal to hit module
1116 //
1117 fElsDR[tID] = GetDRofELoss(t,cdip,fElsDR[tID],normS,ptot,etot);
1118 }
1119 }
1120 // fill info needed to account for ELoss ------------------------------- <<<
1121 //
1122 return kTRUE;
1123}
1124
1125//____________________________________________________
1126Double_t AliITSTPArrayFit::FitHelix(Int_t extQ, Double_t extPT,Double_t extPTerr)
1127{
1128 // fit by helix accounting for the errors of all coordinates (and energy loss if requested)
1129 //
1130 // If extQ is non-0, its sign is imposed as a charge of the track
1131 // If extPT>0 and extPTerr>=0, constrain to measured tr.momentum PT
1132 // with corresponding error (err=0 -> rel.err=1e-6)
1133 //
1134 double chiprev=1e99;
1135 //const Double_t kMaxTEffect = 1E-6;
1136 Double_t dXYZdGlo[3*5],dXYZdLoc[3],xyzRes[3];
1137 //
1138 SetFitDone(kFALSE);
1139 fChi2NDF = -1;
1140 //
1141 if (!FitHelixCrude(extQ)) return -1; // get initial estimate, ignoring the errors
1142 //
1143 if (!IsCovInv()) InvertPointsCovMat(); // prepare inverted errors
1144 if (!fParSol) fParSol = new AliParamSolver(5);
1145 fParSol->SetNGlobal(5);
1146 //
1147 // printf("-1 | %+.2e %+.2e %+.2e %+.2e %+.2e | chi2: %+.4e\n",fParams[0],fParams[1],fParams[2],fParams[3],fParams[4],CalcChi2NDF());
1148 int iter = 0;
1149 fChi2NDF = 1e99;
1150 Bool_t converged = kFALSE;
1151 while(iter<fMaxIter) {
1152 chiprev = fChi2NDF;
1153 fParSol->Clear();
1154 for (int ip=fPntFirst;ip<=fPntLast;ip++) {
1155 //
1156 GetResiduals(xyzRes, ip); // current residuals at point ip
1157 Double_t rrho = fParams[kR0]+fParams[kD0];
1158 Double_t cs0 = TMath::Cos(fParams[kPhi0]);
1159 Double_t sn0 = TMath::Sin(fParams[kPhi0]);
1160 Double_t cst = TMath::Cos(fCurT[ip]+fParams[kPhi0]);
1161 Double_t snt = TMath::Sin(fCurT[ip]+fParams[kPhi0]);
1162 //
1163 int offs = kD0; // dXYZ/dD0
1164 dXYZdGlo[offs + kX] = cs0;
1165 dXYZdGlo[offs + kY] = sn0;
1166 dXYZdGlo[offs + kZ] = 0;
1167 //
1168 offs = kPhi0*3; // dXYZ/dPhi0
1169 dXYZdGlo[offs + kX] = -rrho*sn0;
1170 dXYZdGlo[offs + kY] = rrho*cs0;
1171 dXYZdGlo[offs + kZ] = 0;
1172 //
1173 offs = kR0*3; // dXYZ/dR0
1174 dXYZdGlo[offs + kX] = cs0 - cst;
1175 dXYZdGlo[offs + kY] = sn0 - snt;
1176 dXYZdGlo[offs + kZ] = -fParams[kDip]*fCurT[ip];
1177 //
1178 offs = kDZ*3; // dXYZ/dDZ
1179 dXYZdGlo[offs + kX] = 0;
1180 dXYZdGlo[offs + kY] = 0;
1181 dXYZdGlo[offs + kZ] = 1.;
1182 //
1183 offs = kDip*3; // dXYZ/dDip
1184 dXYZdGlo[offs + kX] = 0;
1185 dXYZdGlo[offs + kY] = 0;
1186 dXYZdGlo[offs + kZ] = -fParams[kR0]*fCurT[ip];
1187 //
1188 dXYZdLoc[kX] = fParams[kR0]*snt;
1189 dXYZdLoc[kY] = -fParams[kR0]*cst;
1190 dXYZdLoc[kZ] = -fParams[kR0]*fParams[kDip];
1191 //
1192 fParSol->AddEquation(dXYZdGlo,dXYZdLoc,xyzRes,GetCovI(ip));
1193 }
1194 //
1195 if (extPT>0) { // add constraint on pt
1196 if (extPTerr<fgkAlmostZero) extPTerr = 1e-6*extPT;
1197 Double_t cf = fBz*GetCharge()*fgkCQConv;
1198 Double_t err2i = extPTerr/cf;
1199 err2i = 1./err2i/err2i;
1200 // printf("Constrain R to %+e\n",extPT/cf);
1201 fParSol->AddConstraint(kR0,-extPT/cf+fParams[kR0],err2i);
1202 }
1203 if (!fParSol->Solve()) { AliError("Failed to fit helix"); return -1; }
1204 Double_t *deltaG = fParSol->GetGlobals();
1205 Double_t *deltaT = fParSol->GetLocals();
1206 for (int ipar=5;ipar--;) fParams[ipar] -= deltaG[ipar];
1207 for (int ip=fPntFirst;ip<=fPntLast;ip++) fCurT[ip] -= deltaT[ip-fPntFirst];
1208 iter++;
1209 //
1210 fChi2NDF = CalcChi2NDF();
1211 // printf("%d | %+.2e %+.2e %+.2e %+.2e %+.2e | chi2: %+.4e %+.4e\n",iter,deltaG[0],deltaG[1],deltaG[2],deltaG[3],deltaG[4],fChi2NDF,fChi2NDF-chiprev);
1212 // printf("->> %+.2e %+.2e %+.2e %+.2e %+.2e | Chi2: %+.6e %+.6e\n",fParams[0],fParams[1],fParams[2],fParams[3],fParams[4],fChi2NDF,fChi2NDF-chiprev);
1213 double difchi2 = chiprev - fChi2NDF;
1214 if ( difchi2<fEps && TMath::Abs(difchi2)<1e-4) {converged = kTRUE; break;}
1215 // if (errT*TMath::Abs(fParams[kR0])<kMaxTEffect && errP<fEps) {converged = kTRUE; break;}
1216 }
1217 //
1218 if (!converged) {
1219 AliDebug(2,Form("Max number of %d iteration reached, Current chi2:%.3e, chi2 change %+.3e",iter,
1220 fChi2NDF,chiprev-fChi2NDF));
1221 for (int ip=fPntFirst;ip<=fPntLast;ip++)
1222 AliDebug(2,Form("P%2d| %+.3e %+.3e %+.3e\n",ip,fPoints->GetX()[ip],fPoints->GetY()[ip],fPoints->GetZ()[ip]));
1223
1224 }
1225 fIter = iter;
1226 SetCharge( fParams[kR0]>0 ? (fBz<0?-1:1):(fBz>0?-1:1) );
1227 SetFitDone();
1228 // printf("F1>> %+.7e %+.7e %+.7e %+.7e %.7e\n",fParams[0],fParams[1],fParams[2],fParams[3],fParams[4]);
1229 //
1230 return fChi2NDF;
1231}
1232
1233//____________________________________________________
1234Double_t AliITSTPArrayFit::FitLine()
1235{
1236 // fit by helix accounting for the errors of all coordinates (and energy loss if requested)
1237 //
1238 double chiprev=1e99;
1239 // const Double_t kMaxTEffect = 1.e-6;
1240 Double_t dXYZdGlo[3*4],dXYZdLoc[3],xyzRes[3];
1241 SetFitDone(kFALSE);
1242 fChi2NDF = -1;
1243 //
1244 if (fParAxis<0) SetParAxis(ChoseParAxis());
1245 //
1246 const float *xyzp[3]={fPoints->GetX(),fPoints->GetY(),fPoints->GetZ()};
1247 if (!IsCovInv()) InvertPointsCovMat();
1248 if (!FitLineCrude()) return -1; // get initial estimate, ignoring the errors
1249 //
1250 if (!fParSol) fParSol = new AliParamSolver(5);
1251 fParSol->SetNGlobal(4);
1252 // initial set of parameters
1253 for (int ip=fPntFirst;ip<=fPntLast;ip++) fCurT[ip] = xyzp[fParAxis][ip]; // use measured param-coordinate
1254 //
1255 int iter = 0;
1256 Bool_t converged = kFALSE;
1257 fChi2NDF = 1e99;
1258 while(iter<fMaxIter) {
1259 chiprev = fChi2NDF;
1260 fParSol->Clear();
1261 for (int ip=fPntFirst;ip<=fPntLast;ip++) {
1262 //
1263 int offs;
1264 GetResiduals(xyzRes, ip); // current residuals at point ip
1265 //
1266 offs = kA0*3; // dXYZ/dA0
1267 dXYZdGlo[offs + fkAxID[kX]] = 1;
1268 dXYZdGlo[offs + fkAxID[kY]] = 0;
1269 dXYZdGlo[offs + fParAxis ] = 0;
1270 //
1271 offs = kB0*3; // dXYZ/dB0
1272 dXYZdGlo[offs + fkAxID[kX]] = fCurT[ip];
1273 dXYZdGlo[offs + fkAxID[kY]] = 0;
1274 dXYZdGlo[offs + fParAxis ] = 0;
1275 //
1276 offs = kA1*3; // dXYZ/dA1
1277 dXYZdGlo[offs + fkAxID[kX]] = 0;
1278 dXYZdGlo[offs + fkAxID[kY]] = 1;
1279 dXYZdGlo[offs + fParAxis ] = 0;
1280 //
1281 offs = kB1*3; // dXYZ/dB1
1282 dXYZdGlo[offs + fkAxID[kX]] = 0;
1283 dXYZdGlo[offs + fkAxID[kY]] = fCurT[ip];
1284 dXYZdGlo[offs + fParAxis ] = 0;
1285 //
1286 dXYZdLoc[ fkAxID[kX] ] = fParams[kB0]; // dX/dt
1287 dXYZdLoc[ fkAxID[kY] ] = fParams[kB1]; // dY/dt
1288 dXYZdLoc[ fParAxis ] = 1;
1289 //
1290 fParSol->AddEquation(dXYZdGlo,dXYZdLoc,xyzRes,GetCovI(ip));
1291 }
1292 //
1293 if (!fParSol->Solve()) { AliError("Failed to fit line"); return -1; }
1294 Double_t *deltaG = fParSol->GetGlobals();
1295 Double_t *deltaT = fParSol->GetLocals();
1296 for (int ipar=4;ipar--;) fParams[ipar] -= deltaG[ipar];
1297 for (int ip=fPntFirst;ip<=fPntLast;ip++) fCurT[ip] -= deltaT[ip-fPntFirst];
1298 iter++;
1299 fChi2NDF = CalcChi2NDF();
1300 // 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);
1301 // printf("->> %+.2e %+.2e %+.2e %+.2e %+.2e | Chi2: %+.6e %+.6e\n",fParams[0],fParams[1],fParams[2],fParams[3],fParams[4],fChi2NDF,fChi2NDF-chiprev);
1302 double difchi2 = chiprev - fChi2NDF;
1303 if ( difchi2<fEps && TMath::Abs(difchi2)<1e-4) {converged = kTRUE; break;}
1304 chiprev = fChi2NDF;
1305 // if (errT<kMaxTEffect && errP<fEps) {converged = kTRUE; break;}
1306 }
1307 //
1308 if (!converged) {
1309 AliDebug(2,Form("Max number of %d iteration reached, Current chi2:%.3e, chi2 change %+.3e",iter,
1310 fChi2NDF,chiprev-fChi2NDF));
1311 for (int ip=fPntFirst;ip<=fPntLast;ip++)
1312 AliDebug(2,Form("P%2d| %+.3e %+.3e %+.3e\n",ip,fPoints->GetX()[ip],fPoints->GetY()[ip],fPoints->GetZ()[ip]));
1313 }
1314 fIter = iter;
1315 SetFitDone();
1316 //printf("F1>> %+.2e %+.2e %+.2e %+.2e\n",fParams[0],fParams[1],fParams[2],fParams[3]);
1317 return fChi2NDF;
1318 //
1319}
1320
1321//____________________________________________________
1322void AliITSTPArrayFit::GetNormal(Double_t *norm, const Float_t *covMat)
1323{
1324 // obtain the lab normal vector to the sensor from the covariance matrix
1325 // in such a way that when the local frame of the sensor coincides with
1326 // the lab frame, the vector {0,1,0} is obtained
1327 Double_t tgxy = TMath::Tan(0.5*TMath::ATan2(2.*covMat[kXY],covMat[kYY]-covMat[kXX]));
1328 Double_t tgyz = TMath::Tan(0.5*TMath::ATan2(2.*covMat[kYZ],covMat[kZZ]-covMat[kYY]));
1329 norm[kY] = 1./TMath::Sqrt(1 + tgxy*tgxy + tgyz*tgyz);
1330 norm[kX] = norm[kY]*tgxy;
1331 norm[kZ] = norm[kY]*tgyz;
1332 //
1333}
1334
1335//____________________________________________________
1336Double_t AliITSTPArrayFit::GetDRofELoss(Double_t t,Double_t cdip,Double_t rhoL,const Double_t *normS,
1337 Double_t &p,Double_t &e) const
1338{
1339 // Calculate energy loss of the particle at given t-param on the layer with rhoL (thickness*density) with
1340 // normal vector normS in the lab. The particle before eloss has energy "e" and momentum "p"
1341 // cdip = cosine of the dip angle = 1/sqrt(1+tgL^2)
1342 // Return the change DR of the radius due to the ELoss
1343 //
1344 // NOTE: with B>0 the negative particles propagate along increasing t-param and positive
1345 // particles - against.
1346 // t-param = 0 corresponds to the point of closest approach of the track to the beam.
1347 // Since the fitted helix parameters of the track are defined in this PCA point, when the correction
1348 // is applied upstream of the PCS, the energy must be increased (DR>0) rather than decreased (DR<0)
1349 //
1350 Double_t dirTr[3];
1351 dirTr[0] = -TMath::Sin(fParams[kPhi0]+t);
1352 dirTr[1] = TMath::Cos(fParams[kPhi0]+t);
1353 dirTr[2] = fParams[kDip];
1354 // cosine of the impact angle
1355 Double_t cosImp = cdip*TMath::Abs(dirTr[0]*normS[0]+dirTr[1]*normS[1]+dirTr[2]*normS[2]);
1356 //
1357 if (cosImp<0.3) cosImp = 0.3; //?
1358 Double_t dE = AliExternalTrackParam::BetheBlochSolid(p/fMass)*rhoL/cosImp;
1359 Double_t dP = e/p*dE;
1360 //
1361 if (t*GetSignQB() < 0) {
1362 dP = -dP;
1363 dE = -dE;
1364 }
1365 //
1366 if (p+dP<0) {
1367 AliInfo(Form("Estimated PLoss %.3f is larger than particle momentum %.3f. Skipping",dP,p));
1368 return 0;
1369 }
1370 //
1371 p += dP;
1372 e += dE;
1373 //
1374 return fCharge*dP*cdip/fBz/fgkCQConv;
1375}
1376
1377//_____________________________________________________________
1378Double_t AliITSTPArrayFit::GetLineOffset(Int_t axis) const
1379{
1380 // return intercept of the parameterization coord = intercept + slope*t for given axis
1381 if (fParAxis<0) return -1E6; // no line fit
1382 if (axis==fParAxis) return 0;
1383 if (fParAxis==kX) return fParams[axis==kY ? kA0 : kA1 ];
1384 if (fParAxis==kY) return fParams[axis==kZ ? kA0 : kA1 ];
1385 return fParams[axis==kX ? kA0 : kA1 ];
1386}
1387
1388//_____________________________________________________________
1389Double_t AliITSTPArrayFit::GetLineSlope(Int_t axis) const
1390{
1391 // return intercept of the parameterization coord = intercept + slope*t for given axis
1392 if (fParAxis<0) return -1E6; // no line fit
1393 if (axis==fParAxis) return 1.;
1394 if (fParAxis==kX) return fParams[axis==kY ? kB0 : kB1 ];
1395 if (fParAxis==kY) return fParams[axis==kZ ? kB0 : kB1 ];
1396 return fParams[axis==kX ? kB0 : kB1 ];
1397}
1398
1399//_____________________________________________________________
1400void AliITSTPArrayFit::Print(Option_t *) const
1401{
1402 const char kCxyz[] = "XYZ";
1403 if (!fPoints) return;
1404 //
1405 printf("Track of %3d points in Bz=%+.1f |Fit ",fPntLast-fPntFirst+1,fBz);
1406 if ( IsFitDone() ) {
1407 if (IsFieldON())
1408 printf("Helix: Chi2: %5.1f | %+.2e %+.2e %+.2e %+.2e %+.2e\n",
1409 fChi2NDF,fParams[kD0],fParams[kPhi0],fParams[kR0],fParams[kDZ],fParams[kDip]);
1410 else
1411 printf("Line%c: Chi2: %5.1f | %+.2e %+.2e %+.2e %+.2e\n",
1412 kCxyz[fParAxis],fChi2NDF,fParams[kA0],fParams[kB0],fParams[kA1],fParams[kB1]);
1413 }
1414 else printf("N/A\n");
1415}
1416
1417
1418
1419
1420//____________________________________________________
1421void AliITSTPArrayFit::BuildMaterialLUT(Int_t ntri)
1422{
1423 // Fill a look-up table with mean material a la AliITSTrackerMI
1424 //
1425 if (!AliGeomManager::GetGeometry()) AliFatal("Geometry is not loaded");
1426 //
1427 // detector layer to check: dX,dZ,Ymin,Ymax
1428 const double kLayr[9][4] = {{0. ,60. , 2.80,3.00}, // beam pipe
1429 {1.28,7.07,-0.20,0.22}, // SPD1
1430 {1.28,7.07,-0.20,0.22}, // SPD2
1431 {0. ,76.0, 10.4,11.8}, // Shield1
1432 {7.02,7.53,-1.00,4.50}, // SDD1
1433 {7.02,7.53,-1.00,4.50}, // SDD2
1434 {0. ,102., 29.0,30.0}, // Shield2
1435 {7.50,4.20,-0.15,4.50}, // SSD1
1436 {7.50,4.20,-0.15,4.50}}; // SSD2
1437 //
1438 //
1439 // build <dens*L> for detectors (track hitting the sensor in normal direction)
1440 double pg1[3],pg2[3],res[7];
1441 //
1442 int sID = 0;
1443 int actLrID = 0;
1444 for (int lr=0;lr<9;lr++) {
1445 //
1446 Bool_t active = kFALSE;
1447 const double* tpars = kLayr[lr];
1448 //
1449 if (IsZero(tpars[0])) { // passive layer
1450 active = kFALSE;
1451 AliInfo(Form("Probing passive layer (total layer #%d)",lr));
1452 }
1453 else {
1454 active = kTRUE;
1455 sID += AliGeomManager::LayerSize(++actLrID);
1456 AliInfo(Form("Probing sensors of active layer #%d (total layers #%d)",actLrID,lr));
1457 }
1458 double shift = TMath::Abs(tpars[2]-tpars[3])*1E-4;
1459 double rhol = 0;
1460 for (int i=ntri;i--;) {
1461 //
1462 if (active) {
1463 int ssID = sID -1 - AliGeomManager::LayerSize(actLrID)*gRandom->Rndm();
1464 pg1[0] = pg2[0] = (gRandom->Rndm()-0.5)*tpars[0] + shift; // local X
1465 pg2[0] -= 2*shift;
1466 pg1[1] = tpars[2];
1467 pg2[1] = tpars[3];
1468 pg1[2] = pg2[2] = (gRandom->Rndm()-0.5)*tpars[1] + shift; // local Z
1469 pg2[2] -= 2*shift;
1470 AliITSgeomTGeo::LocalToGlobal(ssID,pg1,pg1);
1471 AliITSgeomTGeo::LocalToGlobal(ssID,pg2,pg2);
1472 }
1473 else {
1474 double ang = gRandom->Rndm()*TMath::Pi()*2;
1475 pg1[0] = tpars[2]*TMath::Cos(ang)+shift;
1476 pg2[0] = tpars[3]*TMath::Cos(ang)-shift;
1477 pg1[1] = tpars[2]*TMath::Sin(ang);
1478 pg2[1] = tpars[3]*TMath::Sin(ang);
1479 pg1[2] = pg2[2] = (gRandom->Rndm()-0.5)*tpars[1]+shift; // local Z
1480 pg2[2] -= 2*shift;
1481 }
1482
1483 //
1484 AliTracker::MeanMaterialBudget(pg1,pg2,res);
1485 rhol += res[0]*res[4]; // rho*L
1486 }
1487 fgRhoLITS[lr] = rhol/ntri;
1488 AliInfo(Form("Obtained <rho*L> = %e\n",fgRhoLITS[lr]));
1489 }
1490 //
1491 return;
1492}
1493
1494
1495//____________________________________________________
1496Double_t AliITSTPArrayFit::GetPCA2PlaneInfo(Double_t *xyz, Double_t *dir, Int_t axis, Double_t axval) const
1497{
1498 // calculate the PCA to plane normal ti axis and crossing it at axval
1499 // fill the position and direction cosines at this point
1500 //
1501 double xyzp[3] = {0,0,0}; // create fake point
1502 xyzp[axis] = axval;
1503 double covI[6] = {1e-4,0,0,1e-4,0,1e-4}; // fake cov.matrix loose in all directions
1504 covI[4*axis - axis*(axis+1)/2] = 1e8; // except axis
1505 //
1506 double t = GetPosition(xyz, xyzp, covI); // got pca
1507 //
1508 if (dir) GetDirCos(dir,t);
1509 return t;
1510}
1511