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