]>
Commit | Line | Data |
---|---|---|
6be22b3f | 1 | /************************************************************************** |
2 | * Copyright(c) 2009-2011, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | /* $Id$ */ | |
16 | /////////////////////////////////////////////////////////////////////////////////////////////// | |
17 | // // | |
18 | // The line is defined by equations (1) // | |
19 | // a0*z+a1*x-a0*a1=0 and // | |
20 | // b0*z+b1*y-b0*b1=0 // | |
21 | // where x,y,z are NOT the lab axes but z is the lab axis along which the track // | |
22 | // has the largest lever arm and x,y are the remaining 2 axis in // | |
23 | // the order of fgkAxisID[z][0], fgkAxisID[z][1] // | |
24 | // The parameters are fParams[kA0,kB0,kA1,kB1] and the axis chosen as the independent // | |
25 | // var. is fParAxis (i.e. if fParAxis==kZ, then a0=ax,b0=bx, a1=ay,b1=by) // | |
26 | // // | |
27 | // // | |
28 | // The helix is defined by the equations (2) // | |
29 | // X(t) = (dr+R)*cos(phi0) - (R+sum{dRi})*cos(t+phi0) + sum{dRi*cos(phi0+ti)} // | |
30 | // Y(t) = (dr+R)*sin(phi0) - (R+sum{dRi})*sin(t+phi0) + sum{dRi*sin(phi0+ti)} // | |
31 | // Z(t) = dz - (R+sum{dRi})*t*tg(dip) + sum{dRi*ti}*tg(dip) // | |
32 | // where dRi is the change of the radius due to the ELoss at parameter ti // | |
33 | // // | |
34 | // Author: ruben.shahoyan@cern.ch // | |
35 | // // | |
36 | /////////////////////////////////////////////////////////////////////////////////////////////// | |
37 | ||
38 | #include "AliITSTPArrayFit.h" | |
39 | #include "AliExternalTrackParam.h" | |
40 | #include "AliSymMatrix.h" | |
41 | #include "AliLog.h" | |
42 | #include "AliParamSolver.h" | |
43 | #include "AliGeomManager.h" | |
44 | #include "AliITSgeomTGeo.h" | |
45 | #include "AliTracker.h" | |
46 | #include <TRandom.h> | |
ef24eb3b | 47 | #include <TArrayD.h> |
6be22b3f | 48 | |
49 | 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 | AliITSTPArrayFit::AliITSTPArrayFit() : | |
24391cd5 | 89 | fkPoints(0),fParSol(0),fBz(0),fCharge(0),fPntFirst(-1), |
6be22b3f | 90 | fPntLast(-1),fNPBooked(0),fParAxis(-1),fCovI(0),fChi2NDF(0), |
ef24eb3b | 91 | fMaxIter(20),fIter(0),fEps(1e-6),fMass(0),fSwitch2Line(kFALSE),fMaxRforHelix(6.e5), |
1d06ac63 | 92 | fkAxID(0),fkAxCID(0),fCurT(0), |
6be22b3f | 93 | fFirstPosT(0),fNElsPnt(0),fElsId(0),fElsDR(0) |
94 | { | |
95 | // default constructor | |
96 | for (int i=kMaxParam;i--;) fParams[i] = 0; | |
97 | for (int i=kMaxParamSq;i--;) fParamsCov[i] = 0; | |
98 | SetMass(); | |
99 | } | |
100 | ||
101 | //____________________________________________________ | |
102 | AliITSTPArrayFit::AliITSTPArrayFit(Int_t np) : | |
24391cd5 | 103 | fkPoints(0),fParSol(0),fBz(0),fCharge(0),fPntFirst(-1), |
6be22b3f | 104 | fPntLast(-1),fNPBooked(np),fParAxis(-1),fCovI(0),fChi2NDF(0), |
ef24eb3b | 105 | fMaxIter(20),fIter(0),fEps(1e-6),fMass(0),fSwitch2Line(kFALSE),fMaxRforHelix(6.e5), |
1d06ac63 | 106 | fkAxID(0),fkAxCID(0),fCurT(0),fFirstPosT(0),fNElsPnt(0),fElsId(0),fElsDR(0) |
6be22b3f | 107 | { |
108 | // constructor with booking of np points | |
109 | for (int i=kMaxParam;i--;) fParams[i] = 0; | |
110 | for (int i=kMaxParamSq;i--;) fParamsCov[i] = 0; | |
111 | InitAux(); | |
112 | SetEps(); | |
113 | SetMass(); | |
114 | SetMaxIterations(); | |
115 | } | |
116 | ||
117 | //____________________________________________________ | |
118 | AliITSTPArrayFit::AliITSTPArrayFit(const AliITSTPArrayFit &src) : | |
24391cd5 | 119 | TObject(src),fkPoints(src.fkPoints),fParSol(0),fBz(src.fBz), |
6be22b3f | 120 | fCharge(src.fCharge),fPntFirst(src.fPntFirst),fPntLast(src.fPntLast),fNPBooked(src.fNPBooked), |
121 | fParAxis(src.fParAxis),fCovI(0),fChi2NDF(0),fMaxIter(20),fIter(0),fEps(0),fMass(src.fMass), | |
1d06ac63 | 122 | fSwitch2Line(src.fSwitch2Line),fMaxRforHelix(src.fMaxRforHelix),fkAxID(0),fkAxCID(0),fCurT(0), |
123 | fFirstPosT(0),fNElsPnt(0),fElsId(0),fElsDR(0) | |
6be22b3f | 124 | { |
125 | // copy constructor | |
126 | InitAux(); | |
8102b2c9 | 127 | memcpy(fCovI,src.fCovI,fNPBooked*kNCov*sizeof(Double_t)); |
6be22b3f | 128 | for (int i=kMaxParam;i--;) fParams[i] = src.fParams[i]; |
129 | for (int i=kMaxParamSq;i--;) fParamsCov[i] = src.fParamsCov[i]; | |
130 | memcpy(fCurT,src.fCurT,fNPBooked*sizeof(Double_t)); | |
131 | SetEps(src.fEps); | |
132 | SetMaxIterations(src.fMaxIter); | |
133 | // | |
134 | } | |
135 | ||
136 | //____________________________________________________ | |
137 | AliITSTPArrayFit &AliITSTPArrayFit::operator =(const AliITSTPArrayFit& src) | |
138 | { | |
139 | // assignment operator | |
140 | if (this==&src) return *this; | |
141 | ((TObject*)this)->operator=(src); | |
24391cd5 | 142 | fkPoints = src.fkPoints; |
6be22b3f | 143 | if (!fParSol) fParSol = new AliParamSolver(*src.fParSol); |
144 | else *fParSol = *src.fParSol; | |
145 | fBz = src.fBz; | |
146 | fCharge = src.fCharge; | |
147 | fNPBooked = src.fNPBooked; | |
148 | fPntFirst = src.fPntFirst; | |
149 | fPntLast = src.fPntLast; | |
150 | InitAux(); | |
8102b2c9 | 151 | memcpy(fCovI,src.fCovI,fNPBooked*kNCov*sizeof(Double_t)); |
6be22b3f | 152 | for (int i=kMaxParam;i--;) fParams[i] = src.fParams[i]; |
153 | for (int i=kMaxParamSq;i--;) fParamsCov[i] = src.fParamsCov[i]; | |
154 | SetParAxis(src.fParAxis); | |
155 | fNElsPnt = src.fNElsPnt; | |
156 | fFirstPosT = src.fFirstPosT; | |
157 | memcpy(fCurT ,src.fCurT ,fNPBooked*sizeof(Double_t)); | |
158 | memcpy(fElsId ,src.fElsId ,fNPBooked*sizeof(Int_t)); | |
159 | memcpy(fElsDR ,src.fElsDR ,fNPBooked*sizeof(Double_t)); | |
160 | memcpy(fCurT ,src.fCurT ,fNPBooked*sizeof(Double_t)); | |
161 | SetEps(src.fEps); | |
162 | SetMaxIterations(src.fMaxIter); | |
163 | // | |
164 | return *this; | |
165 | // | |
166 | } | |
167 | ||
168 | //____________________________________________________ | |
169 | AliITSTPArrayFit::~AliITSTPArrayFit() | |
170 | { | |
171 | // destructor | |
172 | delete fParSol; | |
173 | delete[] fCovI; | |
174 | delete[] fCurT; | |
175 | delete[] fElsId; | |
176 | delete[] fElsDR; | |
177 | } | |
178 | ||
179 | //____________________________________________________ | |
180 | void AliITSTPArrayFit::Reset() | |
181 | { | |
182 | // reset to process new track | |
183 | if (fParSol) fParSol->Clear(); | |
24391cd5 | 184 | fkPoints=0; |
6be22b3f | 185 | fNElsPnt = 0; |
186 | fFirstPosT = 0; | |
187 | // fBz = 0; | |
188 | fCharge = 0; | |
189 | fIter = 0; | |
190 | fPntFirst=fPntLast=-1; | |
191 | SetParAxis(-1); | |
1d06ac63 | 192 | fSwitch2Line = kFALSE; |
6be22b3f | 193 | ResetBit(kFitDoneBit|kCovInvBit); |
194 | } | |
195 | ||
196 | //____________________________________________________ | |
197 | void AliITSTPArrayFit::AttachPoints(const AliTrackPointArray* points, Int_t pfirst,Int_t plast) | |
198 | { | |
199 | // create from piece of AliTrackPointArray | |
200 | Reset(); | |
24391cd5 | 201 | fkPoints = points; |
6be22b3f | 202 | int np = points->GetNPoints(); |
203 | if (fNPBooked<np) { | |
204 | fNPBooked = np; | |
205 | InitAux(); | |
206 | } | |
207 | fPntFirst = pfirst<0 ? 0 : pfirst; | |
208 | fPntLast = plast<fPntFirst ? np-1 : plast; | |
209 | // | |
210 | for (int i=kMaxParam;i--;) fParams[i] = 0; | |
211 | for (int i=kMaxParamSq;i--;) fParamsCov[i] = 0; | |
212 | // | |
213 | InvertPointsCovMat(); | |
8102b2c9 | 214 | ResetCovIScale(); |
6be22b3f | 215 | // |
216 | } | |
217 | ||
218 | //____________________________________________________ | |
219 | Bool_t AliITSTPArrayFit::SetFirstLast(Int_t pfirst,Int_t plast) | |
220 | { | |
221 | // set first and last point to fit | |
24391cd5 | 222 | const AliTrackPointArray* pnts = fkPoints; |
6be22b3f | 223 | if (!pnts) {AliError("TrackPointArray is not attached yet"); return kFALSE;} |
224 | AttachPoints(pnts,pfirst,plast); | |
225 | return kTRUE; | |
226 | // | |
227 | } | |
228 | ||
229 | //____________________________________________________ | |
230 | Bool_t AliITSTPArrayFit::InvertPointsCovMat() | |
231 | { | |
232 | // invert the cov.matrices of the points | |
233 | for (int i=fPntFirst;i<=fPntLast;i++) { | |
234 | // | |
24391cd5 | 235 | float *cov = (float*)fkPoints->GetCov() + i*6; // pointer on cov.matrix |
6be22b3f | 236 | // |
237 | Double_t t0 = cov[kYY]*cov[kZZ] - cov[kYZ]*cov[kYZ]; | |
238 | Double_t t1 = cov[kXY]*cov[kZZ] - cov[kXZ]*cov[kYZ]; | |
239 | Double_t t2 = cov[kXY]*cov[kYZ] - cov[kXZ]*cov[kYY]; | |
240 | Double_t det = cov[kXX]*t0 - cov[kXY]*t1 + cov[kXZ]*t2; | |
66214d86 | 241 | if (IsZero(det,1e-18)) { // one of errors is 0, fix this |
242 | double norm[3]; | |
243 | TGeoHMatrix hcov; | |
244 | TGeoRotation hrot,hrotI; | |
245 | GetNormal(norm,cov); | |
246 | double phi = TMath::ATan2(norm[1],norm[0]); | |
247 | hrot.SetAngles(-phi*TMath::RadToDeg(),0.,0.); | |
248 | hrotI.SetAngles(phi*TMath::RadToDeg(),0.,0.); | |
249 | // | |
250 | Double_t hcovel[9]; | |
251 | hcovel[0] = cov[kXX]; | |
252 | hcovel[1] = cov[kXY]; | |
253 | hcovel[2] = cov[kXZ]; | |
254 | hcovel[3] = cov[kXY]; | |
255 | hcovel[4] = cov[kYY]; | |
256 | hcovel[5] = cov[kYZ]; | |
257 | hcovel[6] = cov[kXZ]; | |
258 | hcovel[7] = cov[kYZ]; | |
259 | hcovel[8] = cov[kZZ]; | |
260 | hcov.SetRotation(hcovel); | |
261 | // | |
262 | 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 | hcov.Multiply(&hrotI); | |
266 | 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 | if (hcovscl[0]<1e-8) hcovscl[0] = 1e-8; | |
269 | if (hcovscl[4]<1e-8) hcovscl[4] = 1e-8; | |
270 | 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 | hcov.Multiply(&hrot); | |
273 | 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 | cov[kXX] = hcovscl[0]; | |
276 | cov[kXY] = hcovscl[1]; | |
277 | cov[kXZ] = hcovscl[2]; | |
278 | cov[kYY] = hcovscl[4]; | |
279 | cov[kYZ] = hcovscl[5]; | |
280 | cov[kZZ] = hcovscl[8]; | |
281 | } | |
282 | t0 = cov[kYY]*cov[kZZ] - cov[kYZ]*cov[kYZ]; | |
283 | t1 = cov[kXY]*cov[kZZ] - cov[kXZ]*cov[kYZ]; | |
284 | t2 = cov[kXY]*cov[kYZ] - cov[kXZ]*cov[kYY]; | |
285 | det = cov[kXX]*t0 - cov[kXY]*t1 + cov[kXZ]*t2; | |
286 | // | |
abd7ef79 | 287 | AliDebug(2,Form("%+.4e %+.4e %+.4e -> %+.4e",t0,t1,t2,det)); |
66214d86 | 288 | if (IsZero(det,fgkAlmostZero)) { |
6be22b3f | 289 | AliInfo(Form("Cov.Matrix for point %d is singular",i)); |
290 | return kFALSE; | |
291 | } | |
292 | // | |
293 | Double_t *covI = GetCovI(i); | |
294 | covI[kXX] = t0/det; | |
295 | covI[kXY] = -t1/det; | |
296 | covI[kXZ] = t2/det; | |
297 | covI[kYY] = (cov[kXX]*cov[kZZ] - cov[kXZ]*cov[kXZ])/det; | |
298 | covI[kYZ] = (cov[kXY]*cov[kXZ] - cov[kXX]*cov[kYZ])/det; | |
299 | covI[kZZ] = (cov[kXX]*cov[kYY] - cov[kXY]*cov[kXY])/det; | |
300 | // | |
301 | } | |
302 | SetCovInv(); | |
303 | return kTRUE; | |
304 | } | |
305 | ||
306 | //____________________________________________________ | |
307 | void AliITSTPArrayFit::InitAux() | |
308 | { | |
309 | // init auxiliary space | |
310 | if (fCovI) delete[] fCovI; | |
311 | if (fCurT) delete[] fCurT; | |
312 | // | |
8102b2c9 | 313 | fCovI = new Double_t[kNCov*fNPBooked]; |
6be22b3f | 314 | fCurT = new Double_t[fNPBooked+kMaxLrITS]; |
315 | fElsId = new Int_t[fNPBooked+kMaxLrITS]; | |
316 | fElsDR = new Double_t[fNPBooked+kMaxLrITS]; | |
317 | memset(fElsDR,0,(fNPBooked+kMaxLrITS)*sizeof(Double_t)); | |
8102b2c9 | 318 | memset(fCovI,0,fNPBooked*kNCov*sizeof(Double_t)); |
319 | ResetCovIScale(); | |
6be22b3f | 320 | // |
321 | } | |
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 | int np = fPntLast - fPntFirst + 1; | |
335 | if (np<2) { | |
336 | AliError("At least 2 points are needed for straight line fit"); | |
337 | return kFALSE; | |
338 | } | |
339 | // | |
340 | if (fParAxis<0) SetParAxis(ChoseParAxis()); | |
341 | Double_t sZ=0,sZZ=0,sY=0,sYZ=0,sX=0,sXZ=0,det=0; | |
342 | // | |
24391cd5 | 343 | const float *coord[3] = {fkPoints->GetX(),fkPoints->GetY(),fkPoints->GetZ()}; |
6be22b3f | 344 | const Float_t *varZ = coord[ fParAxis ]; |
345 | const Float_t *varX = coord[ fkAxID[kX] ]; | |
346 | const Float_t *varY = coord[ fkAxID[kY] ]; | |
347 | // | |
348 | for (int i=fPntFirst;i<=fPntLast;i++) { | |
349 | sZ += varZ[i]; | |
350 | sZZ += varZ[i]*varZ[i]; | |
351 | // | |
352 | sX += varX[i]; | |
353 | sXZ += varX[i]*varZ[i]; | |
354 | // | |
355 | sY += varY[i]; | |
356 | sYZ += varY[i]*varZ[i]; | |
357 | } | |
358 | det = sZZ*np-sZ*sZ; | |
359 | if (TMath::Abs(det)<fgkAlmostZero) return kFALSE; | |
360 | fParams[0] = (sX*sZZ-sZ*sXZ)/det; | |
361 | fParams[1] = (sXZ*np-sZ*sX)/det; | |
362 | // | |
363 | fParams[2] = (sY*sZZ-sZ*sYZ)/det; | |
364 | fParams[3] = (sYZ*np-sZ*sY)/det; | |
365 | // | |
366 | return kTRUE; | |
367 | // | |
368 | } | |
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 | if (ax>kZ) { | |
375 | AliInfo(Form("Wrong axis choice: %d",ax)); | |
376 | fParAxis = -1; | |
377 | } | |
378 | fParAxis = ax; | |
379 | if (ax>=0) { | |
380 | fkAxID = fgkAxisID[ax]; | |
381 | fkAxCID = fgkAxisCID[ax]; | |
382 | } | |
383 | else { | |
384 | fkAxID = fkAxCID = 0; | |
385 | } | |
386 | // | |
387 | } | |
388 | ||
389 | //____________________________________________________ | |
390 | Int_t AliITSTPArrayFit::ChoseParAxis() const | |
391 | { | |
392 | // select the variable with largest base as a parameter | |
393 | Double_t cmn[3]={1.e9,1.e9,1.e9},cmx[3]={-1.e9,-1.e9,-1.e9}; | |
394 | // | |
24391cd5 | 395 | const float *coord[3] = {fkPoints->GetX(),fkPoints->GetY(),fkPoints->GetZ()}; |
6be22b3f | 396 | for (int i=fPntFirst;i<=fPntLast;i++) { |
397 | for (int j=3;j--;) { | |
398 | Double_t val = coord[j][i]; | |
399 | if (cmn[j]>val) cmn[j] = val; | |
400 | if (cmx[j]<val) cmx[j] = val; | |
401 | } | |
402 | } | |
403 | // | |
404 | int axis = kZ; | |
405 | if (cmx[axis]-cmn[axis] < cmx[kX]-cmn[kX]) axis = kX; | |
406 | if (cmx[axis]-cmn[axis] < cmx[kY]-cmn[kY]) axis = kY; | |
407 | return axis; | |
408 | // | |
409 | } | |
410 | ||
411 | //____________________________________________________ | |
8102b2c9 | 412 | Double_t AliITSTPArrayFit::GetPosition(Double_t *xyzPCA, const Double_t *xyz, const Double_t *covI, Double_t sclCovI) const |
6be22b3f | 413 | { |
414 | // calculate the position of the track at PCA to xyz | |
8102b2c9 | 415 | Double_t t = GetParPCA(xyz,covI,sclCovI); |
6be22b3f | 416 | GetPosition(xyzPCA,t); |
417 | return t; | |
418 | } | |
419 | ||
420 | //____________________________________________________ | |
ef24eb3b | 421 | Double_t AliITSTPArrayFit::GetPosition(Double_t *xyzPCA, const AliTrackPoint *pntCovInv, Bool_t useErr) const |
6be22b3f | 422 | { |
423 | // calculate the position of the track at PCA to pntCovInv | |
424 | // NOTE: the covariance matrix of the point must be inverted | |
ef24eb3b | 425 | double t; |
426 | double xyz[3] = {pntCovInv->GetX(),pntCovInv->GetY(),pntCovInv->GetZ()}; | |
427 | if (useErr) { | |
428 | Double_t covI[6];; | |
429 | for (int i=6;i--;) covI[i] = pntCovInv->GetCov()[i]; | |
430 | t = GetParPCA(xyz,covI); | |
431 | } | |
432 | else t = GetParPCA(xyz); | |
6be22b3f | 433 | GetPosition(xyzPCA,t); |
434 | return t; | |
435 | } | |
436 | ||
437 | //____________________________________________________ | |
ef24eb3b | 438 | void AliITSTPArrayFit::GetResiduals(Double_t *resPCA, const AliTrackPoint *pntCovInv, Bool_t useErr) const |
6be22b3f | 439 | { |
440 | // calculate the residuals of the track at PCA to pntCovInv | |
441 | // NOTE: the covariance matrix of the point must be inverted | |
ef24eb3b | 442 | GetPosition(resPCA,pntCovInv,useErr); |
6be22b3f | 443 | resPCA[0] -= pntCovInv->GetX(); |
444 | resPCA[1] -= pntCovInv->GetY(); | |
445 | resPCA[2] -= pntCovInv->GetZ(); | |
446 | } | |
447 | ||
448 | //____________________________________________________ | |
8102b2c9 | 449 | void AliITSTPArrayFit::GetResiduals(Double_t *resPCA, const Double_t *xyz, const Double_t *covI, Double_t sclCovI) const |
6be22b3f | 450 | { |
451 | // calculate the residuals of the track at PCA to xyz | |
8102b2c9 | 452 | GetPosition(resPCA,xyz,covI,sclCovI); |
6be22b3f | 453 | resPCA[kX] -= xyz[kX]; |
454 | resPCA[kY] -= xyz[kY]; | |
455 | resPCA[kZ] -= xyz[kZ]; | |
456 | } | |
457 | ||
458 | //____________________________________________________ | |
8102b2c9 | 459 | Double_t AliITSTPArrayFit::GetParPCALine(const Double_t *xyz, const Double_t *covI/*, Double_t sclCovI*/) const |
6be22b3f | 460 | { |
461 | // get parameter for the point with least weighted distance to the point | |
462 | // | |
463 | Double_t rhs,denom; | |
464 | Double_t dx = fParams[kA0]-xyz[ fkAxID[kX] ]; | |
465 | Double_t dy = fParams[kA1]-xyz[ fkAxID[kY] ]; | |
466 | Double_t dz = -xyz[ fkAxID[kZ] ]; | |
467 | // | |
468 | if (covI) { | |
469 | Double_t tx = fParams[kB0]*covI[ fkAxCID[kXX] ] + fParams[kB1]*covI[ fkAxCID[kXY] ] + covI[ fkAxCID[kXZ] ]; | |
470 | Double_t ty = fParams[kB0]*covI[ fkAxCID[kXY] ] + fParams[kB1]*covI[ fkAxCID[kYY] ] + covI[ fkAxCID[kYZ] ]; | |
471 | Double_t tz = fParams[kB0]*covI[ fkAxCID[kXZ] ] + fParams[kB1]*covI[ fkAxCID[kYZ] ] + covI[ fkAxCID[kZZ] ]; | |
472 | rhs = tx*dx + ty*dy + tz*dz; | |
473 | denom = -(fParams[kB0]*(covI[ fkAxCID[kXZ] ] + tx) + fParams[kB1]*(covI[ fkAxCID[kYZ] ] + ty) + covI[ fkAxCID[kZZ] ]); | |
474 | } | |
475 | else { | |
476 | rhs = fParams[kB0]*dx + fParams[kB1]*dy + dz; | |
477 | denom = -(fParams[kB0]*fParams[kB0] + fParams[kB1]*fParams[kB1] + 1); | |
478 | } | |
479 | // | |
480 | return rhs/denom; | |
481 | // | |
482 | } | |
483 | ||
484 | //____________________________________________________ | |
8102b2c9 | 485 | void AliITSTPArrayFit::GetDResDPosLine(Double_t *dXYZdP, /*const Double_t *xyz,*/ const Double_t *covI/*,Double_t sclCovI*/) const |
6be22b3f | 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 | Double_t dTdP[3]; | |
8102b2c9 | 491 | GetDtDPosLine(dTdP, /*xyz,*/ covI/*,sclCovI*/); // derivative of the t-param over point position |
6be22b3f | 492 | // |
493 | for (int i=3;i--;) { | |
494 | int var = fkAxID[i]; | |
495 | Double_t *curd = dXYZdP + var*3; // d/dCoord_i | |
496 | curd[ fkAxID[kX] ] = fParams[kB0]*dTdP[var]; | |
497 | curd[ fkAxID[kY] ] = fParams[kB1]*dTdP[var]; | |
498 | curd[ fkAxID[kZ] ] = dTdP[var]; | |
499 | curd[ var ]-= 1.; | |
500 | } | |
501 | // | |
502 | } | |
503 | ||
504 | //____________________________________________________ | |
8102b2c9 | 505 | void AliITSTPArrayFit::GetDResDParamsLine(Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI/*,Double_t sclCovI*/) const |
6be22b3f | 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 | Double_t dTdP[4]; | |
8102b2c9 | 511 | Double_t t = GetDtDParamsLine(dTdP, xyz, covI /*,sclCovI*/); // derivative of the t-param over line params |
6be22b3f | 512 | // |
513 | Double_t *curd = dXYZdP + kA0*3; // d/dA0 | |
514 | curd[ fkAxID[kX] ] = fParams[kB0]*dTdP[kA0] + 1.; | |
515 | curd[ fkAxID[kY] ] = fParams[kB1]*dTdP[kA0]; | |
516 | curd[ fkAxID[kZ] ] = dTdP[kA0]; | |
517 | // | |
518 | curd = dXYZdP + kB0*3; // d/dB0 | |
519 | curd[ fkAxID[kX] ] = fParams[kB0]*dTdP[kB0] + t; | |
520 | curd[ fkAxID[kY] ] = fParams[kB1]*dTdP[kB0]; | |
521 | curd[ fkAxID[kZ] ] = dTdP[kB0]; | |
522 | // | |
523 | curd = dXYZdP + kA1*3; // d/dA1 | |
524 | curd[ fkAxID[kX] ] = fParams[kB0]*dTdP[kA1]; | |
525 | curd[ fkAxID[kY] ] = fParams[kB1]*dTdP[kA1] + 1.; | |
526 | curd[ fkAxID[kZ] ] = dTdP[kA1]; | |
527 | // | |
528 | curd = dXYZdP + kB1*3; // d/dB1 | |
529 | curd[ fkAxID[kX] ] = fParams[kB0]*dTdP[kB1]; | |
530 | curd[ fkAxID[kY] ] = fParams[kB1]*dTdP[kB1] + t; | |
531 | curd[ fkAxID[kZ] ] = dTdP[kB1]; | |
532 | // | |
533 | } | |
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 | Double_t dx = fParams[kA0]-xyz[ fkAxID[kX] ]; | |
542 | Double_t dy = fParams[kA1]-xyz[ fkAxID[kY] ]; | |
543 | Double_t dz = -xyz[ fkAxID[kZ] ]; | |
544 | Double_t rhsDA0,rhsDA1,rhsDB0,rhsDB1,denDB0,denDB1; | |
545 | // | |
546 | if (covI) { | |
547 | Double_t tx = fParams[kB0]*covI[ fkAxCID[kXX] ] + fParams[kB1]*covI[ fkAxCID[kXY] ] + covI[ fkAxCID[kXZ] ]; | |
548 | Double_t ty = fParams[kB0]*covI[ fkAxCID[kXY] ] + fParams[kB1]*covI[ fkAxCID[kYY] ] + covI[ fkAxCID[kYZ] ]; | |
549 | Double_t tz = fParams[kB0]*covI[ fkAxCID[kXZ] ] + fParams[kB1]*covI[ fkAxCID[kYZ] ] + covI[ fkAxCID[kZZ] ]; | |
550 | rhs = tx*dx + ty*dy + tz*dz; | |
551 | denom = -(fParams[kB0]*(covI[ fkAxCID[kXZ] ] + tx) + fParams[kB1]*(covI[ fkAxCID[kYZ] ] + ty) + covI[ fkAxCID[kZZ] ]); | |
552 | // | |
553 | rhsDA0 = tx; | |
554 | rhsDA1 = ty; | |
555 | rhsDB0 = covI[ fkAxCID[kXX] ]*dx + covI[ fkAxCID[kXY] ]*dy + covI[ fkAxCID[kXZ] ]*dz; | |
556 | rhsDB1 = covI[ fkAxCID[kXY] ]*dx + covI[ fkAxCID[kYY] ]*dy + covI[ fkAxCID[kYZ] ]*dz; | |
557 | // | |
558 | denDB0 = -(tx + tx); | |
559 | denDB1 = -(ty + ty); | |
560 | } | |
561 | else { | |
562 | rhs = fParams[kB0]*dx + fParams[kB1]*dy + dz; | |
563 | 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 | denDB0 = -(fParams[kB0]+fParams[kB0]); | |
571 | denDB1 = -(fParams[kB1]+fParams[kB1]); | |
572 | // | |
573 | } | |
574 | // | |
575 | Double_t denom2 = denom*denom; | |
576 | dtparam[kA0] = rhsDA0/denom; // denom does not depend on A0,A1 | |
577 | dtparam[kA1] = rhsDA1/denom; | |
578 | dtparam[kB0] = rhsDB0/denom - rhs/denom2 * denDB0; | |
579 | dtparam[kB1] = rhsDB1/denom - rhs/denom2 * denDB1; | |
580 | // | |
581 | 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 | if (covI) { | |
597 | Double_t tx = fParams[kB0]*covI[ fkAxCID[kXX] ] + fParams[kB1]*covI[ fkAxCID[kXY] ] + covI[ fkAxCID[kXZ] ]; | |
598 | Double_t ty = fParams[kB0]*covI[ fkAxCID[kXY] ] + fParams[kB1]*covI[ fkAxCID[kYY] ] + covI[ fkAxCID[kYZ] ]; | |
599 | 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 | denom = -(fParams[kB0]*(covI[ fkAxCID[kXZ] ] + tx) + fParams[kB1]*(covI[ fkAxCID[kYZ] ] + ty) + covI[ fkAxCID[kZZ] ]); | |
602 | // | |
603 | rhsDX = -tx; | |
604 | rhsDY = -ty; | |
605 | rhsDZ = -tz; | |
606 | } | |
607 | else { | |
608 | // rhs = fParams[kB0]*dx + fParams[kB1]*dy + dz; | |
609 | denom = -(fParams[kB0]*fParams[kB0] + fParams[kB1]*fParams[kB1] + 1); | |
610 | // | |
611 | rhsDX = -fParams[kB0]; | |
612 | rhsDY = -fParams[kB1]; | |
613 | rhsDZ = -1; | |
614 | // | |
615 | } | |
616 | // | |
617 | dtpos[ fkAxID[kX] ] = rhsDX/denom; | |
618 | dtpos[ fkAxID[kY] ] = rhsDY/denom; | |
619 | dtpos[ fkAxID[kZ] ] = rhsDZ/denom; | |
620 | // | |
621 | // return rhs/denom; | |
622 | } | |
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 | if (ipnt<fPntFirst || ipnt>fPntLast) { | |
631 | AliError(Form("Attempt to access the point %d not in the fitted points [%d:%d]",ipnt,fPntFirst,fPntLast)); | |
632 | return; | |
633 | } | |
634 | GetDResDParamsLine(dXYZdP, GetPoint(ipnt) , IsCovIgnored() ? 0 : GetCovI(ipnt)); | |
635 | } | |
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 | if (ipnt<fPntFirst || ipnt>fPntLast) { | |
644 | AliError(Form("Attempt to access the point %d not in the fitted points [%d:%d]",ipnt,fPntFirst,fPntLast)); | |
645 | return; | |
646 | } | |
8102b2c9 | 647 | GetDResDPosLine(dXYZdP,IsCovIgnored() ? 0 : GetCovI(ipnt)/*,GetCovIScale(ipnt)*/); |
6be22b3f | 648 | } |
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 | if (ipnt<fPntFirst || ipnt>fPntLast) { | |
657 | AliError(Form("Attempt to access the point %d not in the fitted points [%d:%d]",ipnt,fPntFirst,fPntLast)); | |
658 | return; | |
659 | } | |
8102b2c9 | 660 | GetDResDParams(dXYZdP, GetPoint(ipnt) , IsCovIgnored() ? 0 : GetCovI(ipnt),GetCovIScale(ipnt)); |
6be22b3f | 661 | } |
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 | if (ipnt<fPntFirst || ipnt>fPntLast) { | |
670 | AliError(Form("Attempt to access the point %d not in the fitted points [%d:%d]",ipnt,fPntFirst,fPntLast)); | |
671 | return; | |
672 | } | |
8102b2c9 | 673 | GetDResDPos(dXYZdP, GetPoint(ipnt), IsCovIgnored() ? 0 : GetCovI(ipnt), GetCovIScale(ipnt)); |
6be22b3f | 674 | } |
675 | ||
676 | //____________________________________________________ | |
8102b2c9 | 677 | void AliITSTPArrayFit::GetDResDParams(Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI, Double_t sclCovI) |
6be22b3f | 678 | { |
679 | // get residual detivatives over the track parameters for the point with least weighted distance to the point | |
680 | // | |
ef24eb3b | 681 | if (!IsHelix()) { // for the straight line calculate analytically |
8102b2c9 | 682 | GetDResDParamsLine(dXYZdP, xyz, covI /*,sclCovI*/); |
6be22b3f | 683 | return; |
684 | } | |
685 | // | |
686 | // calculate derivative numerically | |
687 | const Double_t delta = 0.01; | |
688 | Double_t xyzVar[4][3]; | |
689 | // | |
690 | for (int ipar = 5;ipar--;) { | |
691 | double sav = fParams[ipar]; | |
692 | fParams[ipar] -= delta; | |
8102b2c9 | 693 | GetPosition(xyzVar[0],xyz,covI,sclCovI); |
6be22b3f | 694 | fParams[ipar] += delta/2; |
8102b2c9 | 695 | GetPosition(xyzVar[1],xyz,covI,sclCovI); |
6be22b3f | 696 | fParams[ipar] += delta; |
8102b2c9 | 697 | GetPosition(xyzVar[2],xyz,covI,sclCovI); |
6be22b3f | 698 | fParams[ipar] += delta/2; |
8102b2c9 | 699 | GetPosition(xyzVar[3],xyz,covI,sclCovI); |
6be22b3f | 700 | fParams[ipar] = sav; // restore |
701 | // | |
702 | double *curd = dXYZdP + 3*ipar; | |
703 | 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 | } | |
707 | ||
708 | ||
709 | //____________________________________________________ | |
8102b2c9 | 710 | void AliITSTPArrayFit::GetDResDPos(Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI,Double_t sclCovI) |
6be22b3f | 711 | { |
712 | // get residuals detivative over the point position for the point with least weighted distance to the point | |
713 | // | |
714 | ||
ef24eb3b | 715 | if (!IsHelix()) { // for the straight line calculate analytically |
8102b2c9 | 716 | GetDResDPosLine(dXYZdP, /*xyz,*/ covI /*,sclCovI*/); |
6be22b3f | 717 | return; |
718 | } | |
719 | // | |
720 | // calculate derivative numerically | |
721 | const Double_t delta = 0.005; | |
722 | Double_t xyzVar[4][3]; | |
723 | Double_t xyzv[3] = {xyz[0],xyz[1],xyz[2]}; | |
724 | // | |
725 | for (int ipar = 3;ipar--;) { | |
726 | double sav = xyzv[ipar]; | |
727 | xyzv[ipar] -= delta; | |
8102b2c9 | 728 | GetPosition(xyzVar[0],xyzv,covI,sclCovI); |
6be22b3f | 729 | xyzv[ipar] += delta/2; |
8102b2c9 | 730 | GetPosition(xyzVar[1],xyzv,covI,sclCovI); |
6be22b3f | 731 | xyzv[ipar] += delta; |
8102b2c9 | 732 | GetPosition(xyzVar[2],xyzv,covI,sclCovI); |
6be22b3f | 733 | xyzv[ipar] += delta/2; |
8102b2c9 | 734 | GetPosition(xyzVar[3],xyzv,covI,sclCovI); |
6be22b3f | 735 | xyzv[ipar] = sav; // restore |
736 | // | |
737 | double *curd = dXYZdP + 3*ipar; | |
738 | for (int i=3;i--;) curd[i] = (8.*(xyzVar[2][i]-xyzVar[1][i]) - (xyzVar[3][i]-xyzVar[0][i]))/6./delta; | |
739 | curd[ipar] -= 1.; | |
740 | } | |
741 | // | |
742 | } | |
743 | ||
744 | //________________________________________________________________________________________________________ | |
8102b2c9 | 745 | Double_t AliITSTPArrayFit::GetParPCAHelix(const Double_t* xyz, const Double_t* covI,Double_t sclCovI) const |
6be22b3f | 746 | { |
747 | // find track parameter t (eq.2) corresponding to point of closest approach to xyz | |
748 | // | |
749 | Double_t phi = GetParPCACircle(xyz[kX],xyz[kY]); | |
750 | Double_t cs = TMath::Cos(fParams[kPhi0]); | |
751 | Double_t sn = TMath::Sin(fParams[kPhi0]); | |
752 | Double_t xc = (fParams[kD0]+fParams[kR0])*cs; | |
753 | Double_t yc = (fParams[kD0]+fParams[kR0])*sn; | |
754 | Double_t dchi2,ddchi2; | |
755 | // | |
756 | Double_t dzD = -fParams[kR0]*fParams[kDip]; | |
757 | Double_t dphi = 0; | |
758 | // | |
ef24eb3b | 759 | double rEps = 1e-5/TMath::Abs(fParams[kR0]); // dphi corresponding to 0.1 micron |
760 | if (rEps>fEps) rEps = fEps; | |
761 | // | |
6be22b3f | 762 | int it=0; |
763 | do { | |
764 | cs = TMath::Cos(phi + fParams[kPhi0]); | |
765 | sn = TMath::Sin(phi + fParams[kPhi0]); | |
766 | // | |
767 | Double_t dxD = fParams[kR0]*sn; | |
768 | Double_t dyD = -fParams[kR0]*cs; | |
769 | Double_t dxDD = -dyD; | |
770 | Double_t dyDD = dxD; | |
771 | // | |
772 | Double_t dx = xc - fParams[kR0]*cs - xyz[kX]; | |
773 | Double_t dy = yc - fParams[kR0]*sn - xyz[kY]; | |
774 | Double_t dz = fParams[kDZ] + dzD*phi- xyz[kZ]; | |
775 | // | |
776 | if (covI) { | |
777 | Double_t tx = dx*covI[kXX] + dy*covI[kXY] + dz*covI[kXZ]; | |
778 | Double_t ty = dx*covI[kXY] + dy*covI[kYY] + dz*covI[kYZ]; | |
779 | Double_t tz = dx*covI[kXZ] + dy*covI[kYZ] + dz*covI[kZZ]; | |
780 | // | |
781 | Double_t ttx = dxD*covI[kXX] + dyD*covI[kXY] + dzD*covI[kXZ]; | |
782 | Double_t tty = dxD*covI[kXY] + dyD*covI[kYY] + dzD*covI[kYZ]; | |
783 | Double_t ttz = dxD*covI[kXZ] + dyD*covI[kYZ] + dzD*covI[kZZ]; | |
784 | // | |
785 | // chi2 = dx*tx + dy*ty + dz*tz; | |
786 | dchi2 = dxD*tx + dyD*ty + dzD*tz; | |
787 | ddchi2 = dxDD*tx + dyDD*ty + dxD *ttx + dyD *tty + dzD *ttz; | |
788 | // | |
8102b2c9 | 789 | if (sclCovI>0) {dchi2 *= sclCovI; ddchi2 *= sclCovI;} |
6be22b3f | 790 | } |
791 | else { | |
792 | // chi2 = dx*dx + dy*dy + dz*dz; | |
793 | dchi2 = dxD*dx + dyD*dy + dzD*dz; | |
794 | ddchi2 = dxDD*dx + dyDD*dy + + dxD*dxD + dyD*dyD + dzD*dzD; | |
795 | } | |
796 | // | |
ef24eb3b | 797 | if (TMath::Abs(ddchi2)<fgkAlmostZero || TMath::Abs(dphi=dchi2/ddchi2)<rEps) break; |
6be22b3f | 798 | phi -= dphi; |
799 | } while(++it<fMaxIter); | |
ef24eb3b | 800 | |
6be22b3f | 801 | // |
802 | 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 | Double_t r = fParams[kD0]+fParams[kR0]; | |
811 | Double_t t = TMath::ATan2( r*TMath::Sin(fParams[kPhi0])-y, r*TMath::Cos(fParams[kPhi0])-x ) - fParams[kPhi0]; | |
812 | if (fParams[kR0] < 0) t += TMath::Pi(); | |
813 | if (t > TMath::Pi()) t -= TMath::Pi()*2; | |
814 | if (t <-TMath::Pi()) t += TMath::Pi()*2; | |
815 | 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 | double gam = 1. - (r-fParams[kD0])*(r+fParams[kD0])/fParams[kR0]/(fParams[kD0]+fParams[kR0])/2.; | |
824 | 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 | Double_t dr[3]; // residuals | |
833 | //if (!IsFitDone()) return -1; | |
834 | for (int ipnt=fPntFirst;ipnt<=fPntLast;ipnt++) { | |
835 | GetResiduals(dr,ipnt); | |
836 | Double_t* covI = GetCovI(ipnt); | |
8102b2c9 | 837 | Double_t chi2p = dr[kX]*(dr[kX]*covI[ kXX ]+dr[kY]*covI[ kXY ]+dr[kZ]*covI[ kXZ ]) |
838 | + dr[kY]*(dr[kX]*covI[ kXY ]+dr[kY]*covI[ kYY ]+dr[kZ]*covI[ kYZ ]) | |
839 | + dr[kZ]*(dr[kX]*covI[ kXZ ]+dr[kY]*covI[ kYZ ]+dr[kZ]*covI[ kZZ ]); | |
840 | if (covI[kScl]>0) chi2p *= covI[kScl]; // rescaling was requested for this point's errors | |
841 | chi2 += chi2p; | |
6be22b3f | 842 | } |
843 | int ndf = (fPntLast-fPntFirst+1)*3 - GetNParams(); | |
844 | chi2 /= ndf; | |
845 | return chi2; | |
846 | } | |
847 | ||
848 | //________________________________________________________________________________________________________ | |
849 | void AliITSTPArrayFit::GetResiduals(Double_t *res,Int_t ipnt) const | |
850 | { | |
851 | // calculate residuals at point | |
852 | if (ipnt<fPntFirst || ipnt>fPntLast) { | |
853 | AliError(Form("Attempt to access the point %d not in the fitted points [%d:%d]",ipnt,fPntFirst,fPntLast)); | |
854 | return; | |
855 | } | |
856 | GetPosition(res,fCurT[ipnt]); | |
24391cd5 | 857 | res[kX] -= fkPoints->GetX()[ipnt]; |
858 | res[kY] -= fkPoints->GetY()[ipnt]; | |
859 | res[kZ] -= fkPoints->GetZ()[ipnt]; | |
6be22b3f | 860 | } |
861 | ||
862 | //________________________________________________________________________________________________________ | |
863 | void AliITSTPArrayFit::GetPosition(Double_t *xyz, Double_t t) const | |
864 | { | |
865 | // calculate track position for parameter value t | |
ef24eb3b | 866 | if (IsHelix()) { |
6be22b3f | 867 | // |
868 | Double_t rrho = fParams[kD0]+fParams[kR0]; | |
869 | Double_t xc = rrho*TMath::Cos(fParams[kPhi0]); | |
870 | Double_t yc = rrho*TMath::Sin(fParams[kPhi0]); | |
871 | Double_t r = fParams[kR0]; | |
872 | Double_t ze = 0; | |
873 | // | |
874 | if (IsELossON()) { | |
875 | if (t>0) { | |
876 | for (int i=fFirstPosT;i<fNElsPnt;i++) { // along the track direction | |
877 | int indE = fElsId[i]; | |
878 | if ( t<fCurT[indE] ) break; // does not reach this layer on its way to t | |
879 | xc += fElsDR[indE] * TMath::Cos(fParams[kPhi0] + fCurT[indE]); | |
880 | yc += fElsDR[indE] * TMath::Sin(fParams[kPhi0] + fCurT[indE]); | |
881 | ze += fElsDR[indE] * fCurT[indE]; | |
882 | r += fElsDR[indE]; | |
883 | //printf("ELoss@ %+.2e r:%+.3e got %+.3e\n",fCurT[indE],r,fElsDR[indE]); | |
884 | } | |
885 | } else { | |
886 | for (int i=fFirstPosT;i--;) { // against the track direction | |
887 | int indE = fElsId[i]; | |
888 | if ( t>=fCurT[indE] ) break; // does not reach this layer on its way to t | |
889 | xc += fElsDR[indE] * TMath::Cos(fParams[kPhi0] + fCurT[indE]); | |
890 | yc += fElsDR[indE] * TMath::Sin(fParams[kPhi0] + fCurT[indE]); | |
891 | ze += fElsDR[indE] * fCurT[indE]; | |
892 | r += fElsDR[indE]; | |
893 | //printf("ELoss@ %+.2e r:%+.3e got %+.3e\n",fCurT[indE],r,fElsDR[indE]); | |
894 | } | |
895 | } | |
896 | } | |
897 | // | |
898 | xyz[kZ] = fParams[kDZ] - fParams[kDip]*(t*r - ze); | |
899 | // | |
900 | t += fParams[kPhi0]; | |
901 | xyz[kX] = xc - r*TMath::Cos(t); | |
902 | 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 | } | |
905 | else { | |
906 | xyz[ fkAxID[kX] ] = fParams[kA0] + fParams[kB0]*t; | |
907 | xyz[ fkAxID[kY] ] = fParams[kA1] + fParams[kB1]*t; | |
908 | xyz[ fParAxis ] = t; | |
909 | } | |
910 | } | |
911 | ||
912 | //________________________________________________________________________________________________________ | |
913 | void AliITSTPArrayFit::GetDirCos(Double_t *dircos, Double_t t) const | |
914 | { | |
915 | // calculate track direction cosines for parameter value t | |
ef24eb3b | 916 | if (IsHelix()) { |
6be22b3f | 917 | dircos[kZ] = -fParams[kDip]; |
918 | t += fParams[kPhi0]; | |
919 | dircos[kX] = TMath::Sin(t); | |
920 | dircos[kY] =-TMath::Cos(t); | |
921 | double gam = TMath::Sign(1/TMath::Sqrt(dircos[kZ]*dircos[kZ]+dircos[kY]*dircos[kY]+dircos[kX]*dircos[kX]),fParams[kR0]); | |
922 | for (int i=3;i--;) dircos[i] *= gam; | |
66214d86 | 923 | if (GetSignQB()>0) for (int i=3;i--;) dircos[i] = -dircos[i]; // positive tracks move along decreasing t |
6be22b3f | 924 | } |
925 | else { | |
926 | double gam = 1/TMath::Sqrt( fParams[kB0]*fParams[kB0] + fParams[kB1]*fParams[kB1] + 1.); | |
927 | dircos[ fkAxID[kX] ] = fParams[kB0]*gam; | |
928 | dircos[ fkAxID[kY] ] = fParams[kB1]*gam; | |
929 | dircos[ fParAxis ] = gam; | |
66214d86 | 930 | // decide direction |
931 | if (IsTypeCollision()) { | |
932 | static double xyzF[3],xyzL[3]; | |
933 | GetPosition(xyzF,fPntFirst); | |
934 | GetPosition(xyzL,fPntLast); | |
935 | double dif = fCurT[fPntLast] - fCurT[fPntFirst]; | |
936 | double dr = (xyzL[kX]-xyzF[kX])*(xyzL[kX]+xyzF[kX]) + (xyzL[kY]-xyzF[kY])*(xyzL[kY]+xyzF[kY]); | |
937 | if (dr*dif<0) for (int i=3;i--;) dircos[i] = -dircos[i]; // with increasing t the tracks comes closer to origin | |
938 | } | |
939 | else if (dircos[kY]>0) for (int i=3;i--;) dircos[i] = -dircos[i]; // cosmic tracks have negative angle to Y axis | |
6be22b3f | 940 | } |
66214d86 | 941 | // |
6be22b3f | 942 | } |
943 | ||
944 | //________________________________________________________________________________________________________ | |
945 | Double_t AliITSTPArrayFit::GetMachinePrec() | |
946 | { | |
947 | // estimate machine precision | |
948 | Double_t eps=1.0,a; | |
949 | do { a = 1. + (eps=eps/2.0); } while(a>1.); | |
950 | return TMath::Abs(2.*eps); | |
951 | } | |
952 | ||
ef24eb3b | 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 | static TArrayD arrU,arrV,arrW; | |
962 | double *parrW,*parrU,*parrV; | |
963 | Bool_t eloss = IsELossON(); | |
964 | // | |
965 | int np = fPntLast - fPntFirst + 1; | |
966 | if (np<3) { AliError("At least 3 points are needed for helix fit"); return kFALSE; } | |
967 | // | |
968 | const float *x=fkPoints->GetX(),*y=fkPoints->GetY(),*z=fkPoints->GetZ(),*cov=fkPoints->GetCov(); | |
969 | // | |
970 | if (fPntLast>arrU.GetSize()) { | |
971 | arrU.Set(2*fPntLast); | |
972 | arrV.Set(2*fPntLast); | |
973 | arrW.Set(2*fPntLast); | |
974 | } | |
975 | parrU = arrU.GetArray(); | |
976 | parrV = arrV.GetArray(); | |
977 | 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 | int minRId = fPntFirst; | |
981 | // | |
982 | // get points span | |
983 | double xmn=1e9,xmx=-1e9, ymn=1e9,ymx=-1e9; | |
984 | for (int i=fPntFirst;i<=fPntLast;i++) { | |
985 | parrW[i] = x[i]*x[i]+y[i]*y[i]; | |
986 | if (parrW[i]<parrW[minRId]) minRId = i; // point closest to origin | |
987 | if (xmn>x[i]) xmn = x[i]; | |
988 | if (xmx<x[i]) xmx = x[i]; | |
989 | if (ymn>y[i]) ymn = y[i]; | |
990 | if (ymx<y[i]) ymx = y[i]; | |
991 | } | |
992 | int minRId1 = minRId>fPntFirst ? fPntFirst:fPntFirst+1; | |
993 | for (int i=fPntFirst;i<=fPntLast;i++) if (parrW[i]<parrW[minRId1] && i!=minRId) minRId1 = i; | |
994 | // | |
995 | double xshift = (xmx+xmn)/2 + 10*(ymx-ymn); // shift origin to have uniform weights | |
996 | 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 | for (int i=fPntFirst;i<=fPntLast;i++) { | |
1000 | double xs = x[i] - xshift; | |
1001 | double ys = y[i] - yshift; | |
1002 | double w = xs*xs + ys*ys; | |
1003 | double scl = 1./(1.+w); | |
1004 | int i0 = i-fPntFirst; | |
1005 | wav += parrW[i0] = w*scl; | |
1006 | uav += parrU[i0] = xs*scl; | |
1007 | vav += parrV[i0] = ys*scl; | |
1008 | } | |
1009 | uav /= np; vav /= np; wav /= np; | |
1010 | // | |
1011 | for (int i=fPntFirst;i<=fPntLast;i++) { | |
1012 | // | |
1013 | // point next to closest | |
1014 | int i0 = i-fPntFirst; | |
1015 | if (parrW[i0]<parrW[minRId1-fPntFirst] && i!=minRId) minRId1 = i; | |
1016 | double u = parrU[i0] - uav; | |
1017 | double v = parrV[i0] - vav; | |
1018 | double w = parrW[i0] - wav; | |
1019 | muu += u*u; | |
1020 | muv += u*v; | |
1021 | muw += u*w; | |
1022 | mvv += v*v; | |
1023 | mvw += v*w; | |
1024 | mww += w*w; | |
1025 | } | |
1026 | muu/=np; muv/=np; muw/=np; mvv/=np; mvw/=np; mww/=np; | |
1027 | // | |
1028 | // find eigenvalues: | |
1029 | double trace3 = (muu + mvv + mww)/3.; | |
1030 | double muut = muu-trace3; | |
1031 | double mvvt = mvv-trace3; | |
1032 | double mwwt = mww-trace3; | |
1033 | double q = (muut*(mvvt*mwwt-mvw*mvw) - muv*(muv*mwwt-mvw*muw) + muw*(muv*mvw-mvvt*muw))/2; | |
1034 | double p = (muut*muut+mvvt*mvvt+mwwt*mwwt+2*(muv*muv+muw*muw+mvw*mvw))/6; | |
1035 | double dfpp = p*p*p-q*q; | |
1036 | dfpp = dfpp>0 ? TMath::Sqrt(dfpp)/q : 0; | |
1037 | double ph = TMath::ATan( dfpp )/3.; | |
1038 | if (ph<0) ph += TMath::Pi()/3; | |
1039 | p = p>0 ? TMath::Sqrt(p) : 0; | |
1040 | const double kSqrt3 = 1.73205080; | |
1041 | double snp = TMath::Sin(ph); | |
1042 | double csp = TMath::Cos(ph); | |
1043 | // double eg1 = trace3 + 2*p*csp; | |
1044 | double eg2 = trace3 - p*(csp+kSqrt3*snp); // smallest one | |
1045 | // double eg3 = trace3 - p*(csp-kSqrt3*snp); | |
1046 | // eigenvector for min.eigenvalue | |
1047 | muut = muu-eg2; | |
1048 | mvvt = mvv-eg2; | |
1049 | mwwt = mww-eg2; | |
1050 | double n0 = muv*mvw-muw*mvvt; | |
1051 | double n1 = muv*muw-mvw*muut; | |
1052 | double n2 = muut*mvvt-muv*muv; | |
1053 | // normalize to largest one | |
1054 | double nrm = TMath::Abs(n0); | |
1055 | if (nrm<TMath::Abs(n1)) nrm = TMath::Abs(n1); | |
1056 | if (nrm<TMath::Abs(n2)) nrm = TMath::Abs(n2); | |
1057 | n0/=nrm; n1/=nrm; n2/=nrm; | |
1058 | // | |
1059 | double cpar = -(uav*n0 + vav*n1 + wav*n2); | |
1060 | double xc = -n0/(cpar+n2)/2 + xshift; | |
1061 | double yc = -n1/(cpar+n2)/2 + yshift; | |
1062 | 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 | if (extQ) { | |
1071 | SetCharge(extQ); | |
1072 | sqb = fBz<0 ? -GetCharge():GetCharge(); | |
1073 | } | |
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 | 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 | double v0X = x[minRId1] - x[minRId]; | |
1084 | double v0Y = y[minRId1] - y[minRId]; | |
8102b2c9 | 1085 | sqb = (yc*v0X - xc*v0Y)>0 ? -1:1; |
ef24eb3b | 1086 | } |
1087 | SetCharge( fBz<0 ? -sqb : sqb); | |
1088 | } | |
1089 | // | |
1090 | Double_t phi = TMath::ATan2(yc,xc); | |
1091 | if (sqb<0) phi += TMath::Pi(); | |
1092 | if (phi > TMath::Pi()) phi -= 2.*TMath::Pi(); | |
1093 | else if (phi <-TMath::Pi()) phi += 2.*TMath::Pi(); | |
1094 | fParams[kPhi0] = phi; | |
1095 | fParams[kR0] = sqb<0 ? -rad:rad; | |
1096 | 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 | for (int i=fPntFirst;i<=fPntLast;i++) { | |
1107 | // | |
1108 | Double_t ze2 = cov[i*6 + kZZ]; | |
1109 | Double_t t = TMath::ATan2(yc-y[i],xc-x[i]) - fParams[kPhi0]; // angle at measured z | |
1110 | if (fParams[kR0]<0) t += TMath::Pi(); | |
1111 | if (t > TMath::Pi()) t -= TMath::Pi()*2; | |
1112 | else if (t <-TMath::Pi()) t += TMath::Pi()*2; | |
1113 | if (ze2<fgkAlmostZero) ze2 = 1E-8; | |
1114 | ze2 = 1./ze2; | |
1115 | ss += ze2; | |
1116 | st += t*ze2; | |
1117 | stt+= t*t*ze2; | |
1118 | sz += z[i]*ze2; | |
1119 | szt+= z[i]*t*ze2; | |
1120 | // | |
1121 | 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 | if (eloss) { | |
1124 | double r = TMath::Sqrt(x[i]*x[i]+y[i]*y[i]); | |
1125 | int lr; | |
1126 | for (lr=kMaxLrITS;lr--;) if ( IsZero(r-fgkRLayITS[ lr ],1.) ) break; | |
1127 | if (lr<kMaxLrITS) { | |
1128 | if (t>0) hitLrPos |= (1<<lr); // set bit of the layer | |
1129 | else hitLrNeg |= (1<<lr); // set bit of the layer | |
1130 | } | |
1131 | } | |
1132 | } | |
1133 | double det = ss*stt - st*st; | |
1134 | if (TMath::Abs(det)<fgkAlmostZero) { // no Z dependence | |
1135 | fParams[kDZ] = sz/ss; | |
1136 | fParams[kDip] = 0; | |
1137 | } | |
1138 | else { | |
1139 | fParams[kDZ] = (sz*stt-st*szt)/det; | |
1140 | 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 | if (eloss) { | |
1147 | fNElsPnt = fPntLast - fPntFirst + 1; | |
1148 | // | |
1149 | // to account for the energy loss in the passive volumes, calculate the relevant t-parameters | |
1150 | double* tcur = fCurT + fPntFirst; | |
1151 | double* ecur = fElsDR+ fPntFirst; | |
1152 | // | |
1153 | for (int ilp=3;ilp--;) { | |
1154 | int id = fgkPassivLrITS[ilp]; | |
1155 | double tp = GetHelixParAtR( fgkRLayITS[ id ] ); | |
1156 | if (tp<0) continue; // does not hit this radius | |
1157 | // | |
1158 | tcur[fNElsPnt] = GetSignQB()>0 ? -tp : tp; | |
1159 | ecur[fNElsPnt] = fgRhoLITS[ id ]; | |
1160 | fNElsPnt++; | |
1161 | // printf("Passive on lr %d %+e\n",ilp,tcur[fNElsPnt-1]); | |
1162 | // | |
1163 | if (IsTypeCosmics() && !IsZero(tp)) { // 2 crossings for cosmics | |
1164 | tcur[fNElsPnt] = -tcur[fNElsPnt-1]; | |
1165 | ecur[fNElsPnt] = ecur[fNElsPnt-1]; | |
1166 | fNElsPnt++; | |
1167 | //printf("Passive* on lr %d %+e\n",ilp,-tcur[fNElsPnt-1]); | |
1168 | } | |
1169 | // | |
1170 | } | |
1171 | // check if some active layers did not miss the hit, treat them as passive | |
1172 | for (int ilp=6;ilp--;) { | |
1173 | int id = fgkActiveLrITS[ilp]; | |
1174 | double tp = GetHelixParAtR( fgkRLayITS[ id ] ); | |
1175 | if (tp<0) continue; // does not hit this radius | |
1176 | // | |
1177 | if ( (GetSignQB()>0||IsTypeCosmics()) && !(hitLrNeg & (1<<id)) ) { | |
1178 | tcur[fNElsPnt] = -tp; | |
1179 | ecur[fNElsPnt] = fgRhoLITS[ id ]; | |
1180 | fNElsPnt++; | |
1181 | //printf("Missed on lr %d %+e\n",ilp,-tp); | |
1182 | } | |
1183 | // | |
1184 | if ( (GetSignQB()<0||IsTypeCosmics()) && !(hitLrPos & (1<<id)) ) { | |
1185 | tcur[fNElsPnt] = tp; | |
1186 | ecur[fNElsPnt] = fgRhoLITS[ id ]; | |
1187 | fNElsPnt++; | |
1188 | //printf("Missed* on lr %d %e\n",ilp,tp); | |
1189 | } | |
1190 | } | |
1191 | // | |
1192 | TMath::Sort(fNElsPnt,fCurT+fPntFirst,fElsId,kFALSE); // index e-loss points in increasing order | |
1193 | // find the position of smallest positive t-param | |
1194 | for (fFirstPosT=0;fFirstPosT<fNElsPnt;fFirstPosT++) if (fCurT[ fElsId[ fFirstPosT ] ]>0) break; | |
1195 | // | |
1196 | Double_t cdip = 1./TMath::Sqrt(1.+fParams[kDip]*fParams[kDip]); | |
1197 | Double_t ptot = TMath::Abs(fParams[kR0]*fgkCQConv*fBz/cdip); // momentum and energy | |
1198 | Double_t etot = TMath::Sqrt(ptot*ptot + fMass*fMass); // in the point of closest approach to beam | |
1199 | 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 | for (int ip=fFirstPosT;ip<fNElsPnt;ip++) { | |
1204 | int tID = fElsId[ip]; | |
1205 | Double_t t = fCurT[ tID ]; | |
1206 | // | |
1207 | if (tID>fPntLast) { // this is not a hit layer but passive layer | |
1208 | double php = TMath::ATan2(yc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t), | |
1209 | xc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t)); | |
1210 | normS[0] = -TMath::Cos(php); // normal to the cylinder at intersection point | |
1211 | normS[1] = -TMath::Sin(php); | |
1212 | normS[2] = 0; | |
1213 | } | |
1214 | else GetNormal(normS,fkPoints->GetCov()+tID*6); // vector normal to hit module | |
1215 | 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 | for (int ip=fFirstPosT;ip--;) { | |
1222 | int tID = fElsId[ip]; | |
1223 | Double_t t = fCurT[ tID ]; | |
1224 | // | |
1225 | if (tID>fPntLast) { // this is not a hit layer but passive layer | |
1226 | double php = TMath::ATan2(yc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t), | |
1227 | xc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t)); | |
1228 | normS[0] = -TMath::Cos(php); // normal to the cylinder at intersection point | |
1229 | normS[1] = -TMath::Sin(php); | |
1230 | normS[2] = 0; | |
1231 | } | |
1232 | else GetNormal(normS,fkPoints->GetCov()+tID*6); // vector normal to hit module | |
1233 | // | |
1234 | fElsDR[tID] = GetDRofELoss(t,cdip,fElsDR[tID],normS,ptot,etot); | |
1235 | } | |
1236 | } | |
1237 | // fill info needed to account for ELoss ------------------------------- <<< | |
1238 | // | |
1239 | return kTRUE; | |
1240 | } | |
1241 | ||
1242 | /* | |
6be22b3f | 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; | |
1d06ac63 | 1256 | if (np<3) { AliError("At least 3 points are needed for helix fit"); return kFALSE; } |
6be22b3f | 1257 | // |
24391cd5 | 1258 | const float *x=fkPoints->GetX(),*y=fkPoints->GetY(),*z=fkPoints->GetZ(),*cov=fkPoints->GetCov(); |
6be22b3f | 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 | |
24391cd5 | 1283 | int lrid,volid = fkPoints->GetVolumeID()[i]; |
1284 | if (volid>0) lrid = fgkActiveLrITS[AliGeomManager::VolUIDToLayer(fkPoints->GetVolumeID()[i])-1]; | |
6be22b3f | 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); | |
ef24eb3b | 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 | ||
6be22b3f | 1327 | // |
1d06ac63 | 1328 | Double_t rad = xc*xc + yc*yc - rho2; |
6be22b3f | 1329 | rad = (rad>fgkAlmostZero) ? (TMath::Sqrt(rad)):fgkAlmostZero; |
1330 | // | |
1331 | // printf("Rad: %+e xc: %+e yc: %+e\n",rad,xc,yc); | |
1d06ac63 | 1332 | |
6be22b3f | 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 | // | |
6be22b3f | 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; | |
1d06ac63 | 1357 | fParams[kD0] = xc*TMath::Cos(phi) + yc*TMath::Sin(phi) - fParams[kR0]; |
6be22b3f | 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 | } | |
24391cd5 | 1475 | else GetNormal(normS,fkPoints->GetCov()+tID*6); // vector normal to hit module |
6be22b3f | 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 | } | |
24391cd5 | 1493 | else GetNormal(normS,fkPoints->GetCov()+tID*6); // vector normal to hit module |
6be22b3f | 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 | } | |
ef24eb3b | 1502 | */ |
6be22b3f | 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 | Double_t dXYZdGlo[3*5],dXYZdLoc[3],xyzRes[3]; | |
1515 | // | |
1516 | SetFitDone(kFALSE); | |
1517 | fChi2NDF = -1; | |
1518 | // | |
1519 | if (!FitHelixCrude(extQ)) return -1; // get initial estimate, ignoring the errors | |
1520 | // | |
1d06ac63 | 1521 | if (TMath::Abs(fParams[kR0])>fMaxRforHelix && extPT<=0) { |
1522 | fSwitch2Line = kTRUE; | |
1523 | return FitLine(); | |
1524 | } | |
1525 | // | |
1526 | // | |
6be22b3f | 1527 | if (!IsCovInv()) InvertPointsCovMat(); // prepare inverted errors |
1528 | if (!fParSol) fParSol = new AliParamSolver(5); | |
1529 | 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 | fChi2NDF = 1e99; | |
1534 | Bool_t converged = kFALSE; | |
1535 | while(iter<fMaxIter) { | |
1536 | chiprev = fChi2NDF; | |
1537 | fParSol->Clear(); | |
1538 | for (int ip=fPntFirst;ip<=fPntLast;ip++) { | |
1539 | // | |
1540 | GetResiduals(xyzRes, ip); // current residuals at point ip | |
1541 | Double_t rrho = fParams[kR0]+fParams[kD0]; | |
1542 | Double_t cs0 = TMath::Cos(fParams[kPhi0]); | |
1543 | Double_t sn0 = TMath::Sin(fParams[kPhi0]); | |
1544 | Double_t cst = TMath::Cos(fCurT[ip]+fParams[kPhi0]); | |
1545 | Double_t snt = TMath::Sin(fCurT[ip]+fParams[kPhi0]); | |
1546 | // | |
1547 | int offs = kD0; // dXYZ/dD0 | |
1548 | dXYZdGlo[offs + kX] = cs0; | |
1549 | dXYZdGlo[offs + kY] = sn0; | |
1550 | dXYZdGlo[offs + kZ] = 0; | |
1551 | // | |
1552 | offs = kPhi0*3; // dXYZ/dPhi0 | |
1d06ac63 | 1553 | dXYZdGlo[offs + kX] = -rrho*sn0 + fParams[kR0]*snt; |
1554 | dXYZdGlo[offs + kY] = rrho*cs0 - fParams[kR0]*cst; | |
6be22b3f | 1555 | dXYZdGlo[offs + kZ] = 0; |
1556 | // | |
1557 | offs = kR0*3; // dXYZ/dR0 | |
ef24eb3b | 1558 | if (TMath::Abs(fParams[kR0])<0.9*fMaxRforHelix) { |
1559 | dXYZdGlo[offs + kX] = cs0 - cst; | |
1560 | dXYZdGlo[offs + kY] = sn0 - snt; | |
1561 | dXYZdGlo[offs + kZ] = -fParams[kDip]*fCurT[ip]; | |
1562 | } | |
1563 | else { | |
1564 | dXYZdGlo[offs + kX] = dXYZdGlo[offs + kY] = dXYZdGlo[offs + kZ] = 0; | |
1565 | fParSol->AddConstraint(kR0,0,1.e2); | |
1566 | } | |
6be22b3f | 1567 | // |
1568 | offs = kDZ*3; // dXYZ/dDZ | |
1569 | dXYZdGlo[offs + kX] = 0; | |
1570 | dXYZdGlo[offs + kY] = 0; | |
1571 | dXYZdGlo[offs + kZ] = 1.; | |
1572 | // | |
1573 | offs = kDip*3; // dXYZ/dDip | |
1574 | dXYZdGlo[offs + kX] = 0; | |
1575 | dXYZdGlo[offs + kY] = 0; | |
1576 | dXYZdGlo[offs + kZ] = -fParams[kR0]*fCurT[ip]; | |
1577 | // | |
ef24eb3b | 1578 | // /* |
6be22b3f | 1579 | dXYZdLoc[kX] = fParams[kR0]*snt; |
1580 | dXYZdLoc[kY] = -fParams[kR0]*cst; | |
1581 | dXYZdLoc[kZ] = -fParams[kR0]*fParams[kDip]; | |
ef24eb3b | 1582 | // */ |
1583 | // dXYZdLoc[0] = dXYZdLoc[1] = dXYZdLoc[2] = 0; | |
6be22b3f | 1584 | // |
8102b2c9 | 1585 | fParSol->AddEquation(dXYZdGlo,dXYZdLoc,xyzRes,GetCovI(ip),GetCovIScale(ip)); |
6be22b3f | 1586 | } |
1587 | // | |
1588 | if (extPT>0) { // add constraint on pt | |
1589 | if (extPTerr<fgkAlmostZero) extPTerr = 1e-6*extPT; | |
1590 | Double_t cf = fBz*GetCharge()*fgkCQConv; | |
1591 | Double_t err2i = extPTerr/cf; | |
1592 | err2i = 1./err2i/err2i; | |
1593 | // printf("Constrain R to %+e\n",extPT/cf); | |
1594 | fParSol->AddConstraint(kR0,-extPT/cf+fParams[kR0],err2i); | |
1595 | } | |
1596 | if (!fParSol->Solve()) { AliError("Failed to fit helix"); return -1; } | |
1597 | Double_t *deltaG = fParSol->GetGlobals(); | |
ef0526f3 | 1598 | // Double_t *deltaT = fParSol->GetLocals(); |
6be22b3f | 1599 | for (int ipar=5;ipar--;) fParams[ipar] -= deltaG[ipar]; |
ef24eb3b | 1600 | // |
1601 | if (TMath::Abs(fParams[kR0])>0.9*fMaxRforHelix) fParams[kR0] = TMath::Sign(0.9*fMaxRforHelix,fParams[kR0]); | |
1602 | // | |
1d06ac63 | 1603 | for (int ip=fPntFirst;ip<=fPntLast;ip++) { |
1604 | 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 | } | |
6be22b3f | 1608 | iter++; |
1609 | // | |
1610 | fChi2NDF = CalcChi2NDF(); | |
ef24eb3b | 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); | |
6be22b3f | 1613 | double difchi2 = chiprev - fChi2NDF; |
1614 | 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 | } | |
1617 | // | |
1618 | if (!converged) { | |
1619 | AliDebug(2,Form("Max number of %d iteration reached, Current chi2:%.3e, chi2 change %+.3e",iter, | |
1620 | fChi2NDF,chiprev-fChi2NDF)); | |
1621 | for (int ip=fPntFirst;ip<=fPntLast;ip++) | |
24391cd5 | 1622 | AliDebug(2,Form("P%2d| %+.3e %+.3e %+.3e\n",ip,fkPoints->GetX()[ip],fkPoints->GetY()[ip],fkPoints->GetZ()[ip])); |
6be22b3f | 1623 | |
1624 | } | |
1625 | fIter = iter; | |
1626 | SetCharge( fParams[kR0]>0 ? (fBz<0?-1:1):(fBz>0?-1:1) ); | |
1627 | SetFitDone(); | |
1628 | // printf("F1>> %+.7e %+.7e %+.7e %+.7e %.7e\n",fParams[0],fParams[1],fParams[2],fParams[3],fParams[4]); | |
1629 | // | |
1630 | return fChi2NDF; | |
1631 | } | |
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 | Double_t dXYZdGlo[3*4],dXYZdLoc[3],xyzRes[3]; | |
1641 | SetFitDone(kFALSE); | |
1642 | fChi2NDF = -1; | |
1643 | // | |
1644 | if (fParAxis<0) SetParAxis(ChoseParAxis()); | |
1645 | // | |
24391cd5 | 1646 | const float *xyzp[3]={fkPoints->GetX(),fkPoints->GetY(),fkPoints->GetZ()}; |
6be22b3f | 1647 | if (!IsCovInv()) InvertPointsCovMat(); |
1648 | if (!FitLineCrude()) return -1; // get initial estimate, ignoring the errors | |
1649 | // | |
1650 | if (!fParSol) fParSol = new AliParamSolver(5); | |
1651 | fParSol->SetNGlobal(4); | |
1652 | // initial set of parameters | |
1653 | 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 | fChi2NDF = 1e99; | |
1658 | while(iter<fMaxIter) { | |
1659 | chiprev = fChi2NDF; | |
1660 | fParSol->Clear(); | |
1661 | for (int ip=fPntFirst;ip<=fPntLast;ip++) { | |
1662 | // | |
1663 | int offs; | |
1664 | GetResiduals(xyzRes, ip); // current residuals at point ip | |
1665 | // | |
1666 | offs = kA0*3; // dXYZ/dA0 | |
1667 | dXYZdGlo[offs + fkAxID[kX]] = 1; | |
1668 | dXYZdGlo[offs + fkAxID[kY]] = 0; | |
1669 | dXYZdGlo[offs + fParAxis ] = 0; | |
1670 | // | |
1671 | offs = kB0*3; // dXYZ/dB0 | |
1672 | dXYZdGlo[offs + fkAxID[kX]] = fCurT[ip]; | |
1673 | dXYZdGlo[offs + fkAxID[kY]] = 0; | |
1674 | dXYZdGlo[offs + fParAxis ] = 0; | |
1675 | // | |
1676 | offs = kA1*3; // dXYZ/dA1 | |
1677 | dXYZdGlo[offs + fkAxID[kX]] = 0; | |
1678 | dXYZdGlo[offs + fkAxID[kY]] = 1; | |
1679 | dXYZdGlo[offs + fParAxis ] = 0; | |
1680 | // | |
1681 | offs = kB1*3; // dXYZ/dB1 | |
1682 | dXYZdGlo[offs + fkAxID[kX]] = 0; | |
1683 | dXYZdGlo[offs + fkAxID[kY]] = fCurT[ip]; | |
1684 | dXYZdGlo[offs + fParAxis ] = 0; | |
1685 | // | |
1686 | dXYZdLoc[ fkAxID[kX] ] = fParams[kB0]; // dX/dt | |
1687 | dXYZdLoc[ fkAxID[kY] ] = fParams[kB1]; // dY/dt | |
1688 | dXYZdLoc[ fParAxis ] = 1; | |
1689 | // | |
8102b2c9 | 1690 | fParSol->AddEquation(dXYZdGlo,dXYZdLoc,xyzRes,GetCovI(ip),GetCovIScale(ip)); |
6be22b3f | 1691 | } |
1692 | // | |
1693 | if (!fParSol->Solve()) { AliError("Failed to fit line"); return -1; } | |
1694 | Double_t *deltaG = fParSol->GetGlobals(); | |
1695 | Double_t *deltaT = fParSol->GetLocals(); | |
1696 | for (int ipar=4;ipar--;) fParams[ipar] -= deltaG[ipar]; | |
1697 | for (int ip=fPntFirst;ip<=fPntLast;ip++) fCurT[ip] -= deltaT[ip-fPntFirst]; | |
1698 | iter++; | |
1699 | 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 | double difchi2 = chiprev - fChi2NDF; | |
1703 | if ( difchi2<fEps && TMath::Abs(difchi2)<1e-4) {converged = kTRUE; break;} | |
1704 | chiprev = fChi2NDF; | |
1705 | // if (errT<kMaxTEffect && errP<fEps) {converged = kTRUE; break;} | |
1706 | } | |
1707 | // | |
1708 | if (!converged) { | |
1709 | AliDebug(2,Form("Max number of %d iteration reached, Current chi2:%.3e, chi2 change %+.3e",iter, | |
1710 | fChi2NDF,chiprev-fChi2NDF)); | |
1711 | for (int ip=fPntFirst;ip<=fPntLast;ip++) | |
24391cd5 | 1712 | AliDebug(2,Form("P%2d| %+.3e %+.3e %+.3e\n",ip,fkPoints->GetX()[ip],fkPoints->GetY()[ip],fkPoints->GetZ()[ip])); |
6be22b3f | 1713 | } |
1714 | fIter = iter; | |
1715 | SetFitDone(); | |
1716 | //printf("F1>> %+.2e %+.2e %+.2e %+.2e\n",fParams[0],fParams[1],fParams[2],fParams[3]); | |
1717 | return fChi2NDF; | |
1718 | // | |
1719 | } | |
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 | Double_t tgxy = TMath::Tan(0.5*TMath::ATan2(2.*covMat[kXY],covMat[kYY]-covMat[kXX])); | |
1728 | Double_t tgyz = TMath::Tan(0.5*TMath::ATan2(2.*covMat[kYZ],covMat[kZZ]-covMat[kYY])); | |
1729 | norm[kY] = 1./TMath::Sqrt(1 + tgxy*tgxy + tgyz*tgyz); | |
1730 | norm[kX] = norm[kY]*tgxy; | |
1731 | norm[kZ] = norm[kY]*tgyz; | |
1732 | // | |
1733 | } | |
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 | dirTr[0] = -TMath::Sin(fParams[kPhi0]+t); | |
1752 | dirTr[1] = TMath::Cos(fParams[kPhi0]+t); | |
1753 | dirTr[2] = fParams[kDip]; | |
1754 | // cosine of the impact angle | |
1755 | Double_t cosImp = cdip*TMath::Abs(dirTr[0]*normS[0]+dirTr[1]*normS[1]+dirTr[2]*normS[2]); | |
1756 | // | |
1757 | if (cosImp<0.3) cosImp = 0.3; //? | |
1758 | Double_t dE = AliExternalTrackParam::BetheBlochSolid(p/fMass)*rhoL/cosImp; | |
1759 | Double_t dP = e/p*dE; | |
1760 | // | |
1761 | if (t*GetSignQB() < 0) { | |
1762 | dP = -dP; | |
1763 | dE = -dE; | |
1764 | } | |
1765 | // | |
1766 | if (p+dP<0) { | |
1767 | AliInfo(Form("Estimated PLoss %.3f is larger than particle momentum %.3f. Skipping",dP,p)); | |
1768 | return 0; | |
1769 | } | |
1770 | // | |
1771 | p += dP; | |
1772 | e += dE; | |
1773 | // | |
1774 | return fCharge*dP*cdip/fBz/fgkCQConv; | |
1775 | } | |
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 | if (fParAxis<0) return -1E6; // no line fit | |
1782 | if (axis==fParAxis) return 0; | |
1783 | if (fParAxis==kX) return fParams[axis==kY ? kA0 : kA1 ]; | |
1784 | if (fParAxis==kY) return fParams[axis==kZ ? kA0 : kA1 ]; | |
1785 | return fParams[axis==kX ? kA0 : kA1 ]; | |
1786 | } | |
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 | if (fParAxis<0) return -1E6; // no line fit | |
1793 | if (axis==fParAxis) return 1.; | |
1794 | if (fParAxis==kX) return fParams[axis==kY ? kB0 : kB1 ]; | |
1795 | if (fParAxis==kY) return fParams[axis==kZ ? kB0 : kB1 ]; | |
1796 | return fParams[axis==kX ? kB0 : kB1 ]; | |
1797 | } | |
1798 | ||
1799 | //_____________________________________________________________ | |
1800 | void AliITSTPArrayFit::Print(Option_t *) const | |
1801 | { | |
24391cd5 | 1802 | // print results of the fit |
1803 | // | |
6be22b3f | 1804 | const char kCxyz[] = "XYZ"; |
24391cd5 | 1805 | if (!fkPoints) return; |
6be22b3f | 1806 | // |
1807 | printf("Track of %3d points in Bz=%+.1f |Fit ",fPntLast-fPntFirst+1,fBz); | |
1808 | if ( IsFitDone() ) { | |
ef24eb3b | 1809 | if (IsHelix()) |
6be22b3f | 1810 | printf("Helix: Chi2: %5.1f | %+.2e %+.2e %+.2e %+.2e %+.2e\n", |
1811 | fChi2NDF,fParams[kD0],fParams[kPhi0],fParams[kR0],fParams[kDZ],fParams[kDip]); | |
1812 | else | |
1813 | printf("Line%c: Chi2: %5.1f | %+.2e %+.2e %+.2e %+.2e\n", | |
1814 | kCxyz[fParAxis],fChi2NDF,fParams[kA0],fParams[kB0],fParams[kA1],fParams[kB1]); | |
1815 | } | |
1816 | else printf("N/A\n"); | |
1817 | } | |
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 | 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 | double pg1[3],pg2[3],res[7]; | |
1843 | // | |
1844 | int sID = 0; | |
1845 | int actLrID = 0; | |
1846 | for (int lr=0;lr<9;lr++) { | |
1847 | // | |
1848 | Bool_t active = kFALSE; | |
1849 | const double* tpars = kLayr[lr]; | |
1850 | // | |
1851 | if (IsZero(tpars[0])) { // passive layer | |
1852 | active = kFALSE; | |
1853 | AliInfo(Form("Probing passive layer (total layer #%d)",lr)); | |
1854 | } | |
1855 | else { | |
1856 | active = kTRUE; | |
1857 | sID += AliGeomManager::LayerSize(++actLrID); | |
1858 | AliInfo(Form("Probing sensors of active layer #%d (total layers #%d)",actLrID,lr)); | |
1859 | } | |
1860 | double shift = TMath::Abs(tpars[2]-tpars[3])*1E-4; | |
1861 | double rhol = 0; | |
1862 | for (int i=ntri;i--;) { | |
1863 | // | |
1864 | if (active) { | |
1865 | int ssID = sID -1 - AliGeomManager::LayerSize(actLrID)*gRandom->Rndm(); | |
1866 | pg1[0] = pg2[0] = (gRandom->Rndm()-0.5)*tpars[0] + shift; // local X | |
1867 | pg2[0] -= 2*shift; | |
1868 | pg1[1] = tpars[2]; | |
1869 | pg2[1] = tpars[3]; | |
1870 | pg1[2] = pg2[2] = (gRandom->Rndm()-0.5)*tpars[1] + shift; // local Z | |
1871 | pg2[2] -= 2*shift; | |
1872 | AliITSgeomTGeo::LocalToGlobal(ssID,pg1,pg1); | |
1873 | AliITSgeomTGeo::LocalToGlobal(ssID,pg2,pg2); | |
1874 | } | |
1875 | else { | |
1876 | double ang = gRandom->Rndm()*TMath::Pi()*2; | |
1877 | pg1[0] = tpars[2]*TMath::Cos(ang)+shift; | |
1878 | pg2[0] = tpars[3]*TMath::Cos(ang)-shift; | |
1879 | pg1[1] = tpars[2]*TMath::Sin(ang); | |
1880 | pg2[1] = tpars[3]*TMath::Sin(ang); | |
1881 | pg1[2] = pg2[2] = (gRandom->Rndm()-0.5)*tpars[1]+shift; // local Z | |
1882 | pg2[2] -= 2*shift; | |
1883 | } | |
1884 | ||
1885 | // | |
1886 | AliTracker::MeanMaterialBudget(pg1,pg2,res); | |
1887 | rhol += res[0]*res[4]; // rho*L | |
1888 | } | |
1889 | fgRhoLITS[lr] = rhol/ntri; | |
1890 | AliInfo(Form("Obtained <rho*L> = %e\n",fgRhoLITS[lr])); | |
1891 | } | |
1892 | // | |
1893 | return; | |
1894 | } | |
1895 | ||
6be22b3f | 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 | double xyzp[3] = {0,0,0}; // create fake point | |
1903 | xyzp[axis] = axval; | |
1904 | double covI[6] = {1e-4,0,0,1e-4,0,1e-4}; // fake cov.matrix loose in all directions | |
1905 | covI[4*axis - axis*(axis+1)/2] = 1e8; // except axis | |
1906 | // | |
1907 | double t = GetPosition(xyz, xyzp, covI); // got pca | |
1908 | // | |
1909 | if (dir) GetDirCos(dir,t); | |
1910 | return t; | |
1911 | } | |
1912 | ||
66214d86 | 1913 | //____________________________________________________ |
1914 | void AliITSTPArrayFit::GetT0Info(Double_t* xyz, Double_t *dir) const | |
1915 | { | |
1916 | // get direction cosines for t = 0; | |
1917 | GetPosition(xyz, 0); | |
1918 | if (dir) GetDirCos(dir,0); | |
1919 | } | |
1920 | ||
1921 | //____________________________________________________ | |
1922 | Bool_t AliITSTPArrayFit::CalcErrorMatrix() | |
1923 | { | |
1924 | // compute covariance matrix in lenear approximation of residuals on parameters | |
1925 | static AliSymMatrix cv(5); | |
1926 | static Double_t dres[5][3]; | |
1927 | cv.Reset(); | |
ef24eb3b | 1928 | int npar = IsHelix() ? 5:4; |
66214d86 | 1929 | // |
1930 | for (int ip=fPntFirst;ip<=fPntLast;ip++) { | |
1931 | GetDResDParams(&dres[0][0],ip); | |
1932 | Double_t* covI = GetCovI(ip); | |
1933 | for (int i=npar;i--;) | |
8102b2c9 | 1934 | for (int j=i+1;j--;) { |
1935 | double cvadd = dres[i][kX]*(dres[j][kX]*covI[ kXX ] + dres[j][kY]*covI[ kXY ] + dres[j][kZ]*covI[ kXZ ]) | |
1936 | + dres[i][kY]*(dres[j][kX]*covI[ kXY ] + dres[j][kY]*covI[ kYY ] + dres[j][kZ]*covI[ kYZ ]) | |
1937 | + dres[i][kZ]*(dres[j][kX]*covI[ kXZ ] + dres[j][kY]*covI[ kYZ ] + dres[j][kZ]*covI[ kZZ ]); | |
1938 | if (covI[kScl]>0) cvadd *= covI[kScl]; | |
1939 | cv(i,j) += cvadd; | |
1940 | } | |
66214d86 | 1941 | } |
1942 | cv.SetSizeUsed(npar); | |
1943 | if (cv.InvertChol()) { | |
ef24eb3b | 1944 | //cv.Print("l"); |
66214d86 | 1945 | int cnt = 0; |
1946 | for (int i=npar;i--;) for (int j=i+1;j--;)fParamsCov[cnt++] = cv(i,j); | |
1947 | return kTRUE; | |
1948 | } | |
1949 | // | |
1950 | return kFALSE; | |
1951 | } | |
1d06ac63 | 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 | const double *xyz = GetPoint(ipnt); | |
1958 | const double *covI = GetCovI(ipnt); | |
8102b2c9 | 1959 | if (IsHelix()) return GetParPCAHelix(xyz,covI,covI[kScl]); |
1960 | else return GetParPCALine(xyz,covI/*,covI[kScl]*/); | |
1d06ac63 | 1961 | } |
24391cd5 | 1962 | |
ef24eb3b | 1963 | //____________________________________________________ |
1964 | Double_t AliITSTPArrayFit::GetPt() const | |
1965 | { | |
1966 | return IsFieldON()&&IsHelix() ? TMath::Abs(fParams[kR0]*fBz*fgkCQConv) : -1; | |
1967 | } | |
1968 | ||
1969 | //____________________________________________________ | |
1970 | Double_t AliITSTPArrayFit::GetP() const | |
1971 | { | |
1972 | if (!IsFieldON()) return -1; | |
1973 | Double_t cdip = 1./TMath::Sqrt(1.+fParams[kDip]*fParams[kDip]); | |
1974 | return TMath::Abs(fParams[kR0]*fgkCQConv*fBz/cdip); // momentum | |
1975 | } | |
24391cd5 | 1976 |