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