]>
Commit | Line | Data |
---|---|---|
f77f13c8 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, 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 | ||
16 | /* $Id$ */ | |
17 | ||
853a0f19 | 18 | // AliITSNeuralTrack |
19 | // | |
20 | // The format of output data from Neural Tracker | |
21 | // It can export data in the format of AliITSIOTrack | |
22 | // (V1) tracking. | |
23 | // Compatibility adaptation to V2 tracking is on the way. | |
24 | // Author: A. Pulvirenti | |
25 | ||
b9d722bc | 26 | #include <Riostream.h> |
853a0f19 | 27 | //#include <cstdlib> |
28 | //#include <cstring> | |
b9d722bc | 29 | |
853a0f19 | 30 | //#include <TObject.h> |
31 | //#include <TROOT.h> | |
b9d722bc | 32 | #include <TMath.h> |
33 | #include <TString.h> | |
853a0f19 | 34 | //#include <TObjArray.h> |
35 | //#include <TH1.h> | |
b9d722bc | 36 | #include <TMatrixD.h> |
49370114 | 37 | #if ROOT_VERSION_CODE >= 262146 |
38 | #include <TMatrixDEigen.h> | |
39 | #endif | |
b9d722bc | 40 | |
41 | //#include "AliITSVertex.h" | |
42 | #include "AliITSIOTrack.h" | |
43 | #include "AliITSNeuralPoint.h" | |
44 | ||
45 | #include "AliITSNeuralTrack.h" | |
46 | ||
47 | ||
48 | ||
49 | ClassImp(AliITSNeuralTrack) | |
50 | // | |
51 | // | |
52 | // | |
53 | AliITSNeuralTrack::AliITSNeuralTrack() : fMatrix(5,5), fVertex() | |
54 | { | |
853a0f19 | 55 | // Default constructor |
56 | ||
b9d722bc | 57 | Int_t i; |
58 | ||
59 | fMass = 0.1396; // default assumption: pion | |
60 | fField = 2.0; // default assumption: B = 0.4 Tesla | |
61 | ||
62 | fXC = fYC = fR = fC = 0.0; | |
63 | fTanL = fG0 = fDt = fDz = 0.0; | |
64 | fStateR = fStatePhi = fStateZ = fChi2 = fNSteps = 0.0; | |
65 | ||
66 | fLabel = 0; | |
67 | fCount = 0; | |
68 | for (i = 0; i < 6; i++) fPoint[i] = 0; | |
69 | ||
70 | fVertex.X() = 0.0; | |
71 | fVertex.Y() = 0.0; | |
72 | fVertex.Z() = 0.0; | |
73 | fVertex.ErrX() = 0.0; | |
74 | fVertex.ErrY() = 0.0; | |
75 | fVertex.ErrZ() = 0.0; | |
76 | } | |
77 | // | |
78 | // | |
79 | // | |
853a0f19 | 80 | AliITSNeuralTrack::AliITSNeuralTrack(AliITSNeuralTrack &track) |
81 | : TObject((TObject&)track), fMatrix(5,5), fVertex() | |
82 | { | |
83 | // copy constructor | |
84 | ||
85 | Int_t i; | |
86 | ||
87 | fMass = 0.1396; // default assumption: pion | |
88 | fField = 2.0; // default assumption: B = 0.4 Tesla | |
89 | ||
90 | fXC = fYC = fR = fC = 0.0; | |
91 | fTanL = fG0 = fDt = fDz = 0.0; | |
92 | fStateR = fStatePhi = fStateZ = fChi2 = fNSteps = 0.0; | |
93 | ||
94 | fLabel = 0; | |
95 | fCount = 0; | |
96 | for (i = 0; i < 6; i++) fPoint[i] = track.fPoint[i]; | |
97 | ||
98 | fVertex.X() = 0.0; | |
99 | fVertex.Y() = 0.0; | |
100 | fVertex.Z() = 0.0; | |
101 | fVertex.ErrX() = 0.0; | |
102 | fVertex.ErrY() = 0.0; | |
103 | fVertex.ErrZ() = 0.0; | |
104 | } | |
105 | // | |
106 | // | |
107 | // | |
b9d722bc | 108 | AliITSNeuralTrack::~AliITSNeuralTrack() |
109 | { | |
110 | Int_t l; | |
111 | for (l = 0; l < 6; l++) fPoint[l] = 0; | |
112 | } | |
113 | // | |
114 | // | |
115 | // | |
116 | void AliITSNeuralTrack::AssignLabel() | |
853a0f19 | 117 | { |
b9d722bc | 118 | // Assigns a GEANT label to the found track. |
119 | // Every cluster has up to three labels (it can have less). Then each label is | |
120 | // recorded for each point. Then, counts are made to check if some of the labels | |
121 | // appear more than once. Finally, the label which appears most times is assigned | |
122 | // to the track in the field fLabel. | |
123 | // The number of points containing that label is counted in the fCount data-member. | |
853a0f19 | 124 | |
b9d722bc | 125 | Bool_t found; |
126 | Int_t i, j, l, lab, max = 0; | |
127 | ||
128 | // We have up to 6 points for 3 labels each => up to 18 possible different values | |
129 | Int_t idx[18] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; | |
130 | Int_t count[18] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; | |
131 | ||
132 | for (l = 0; l < 6; l++) { | |
133 | if (!fPoint[l]) continue; | |
134 | // Sometimes the same label appears two times in the same recpoint. | |
135 | // With these if statements, such problem is solved by turning | |
136 | // one of them to -1. | |
137 | if (fPoint[l]->GetLabel(1) >= 0 && fPoint[l]->GetLabel(1) == fPoint[l]->GetLabel(0)) | |
138 | fPoint[l]->SetLabel(1, -1); | |
139 | if (fPoint[l]->GetLabel(2) >= 0 && fPoint[l]->GetLabel(2) == fPoint[l]->GetLabel(0)) | |
140 | fPoint[l]->SetLabel(2, -1); | |
141 | if (fPoint[l]->GetLabel(2) >= 0 && fPoint[l]->GetLabel(2) == fPoint[l]->GetLabel(1)) | |
142 | fPoint[l]->SetLabel(2, -1); | |
143 | for (i = 0; i < 3; i++) { | |
144 | lab = fPoint[l]->GetLabel(i); | |
145 | if (lab < 0) continue; | |
146 | found = kFALSE; | |
147 | for (j = 0; j < max; j++) { | |
148 | if (idx[j] == lab) { | |
149 | count[j]++; | |
150 | found = kTRUE; | |
151 | } | |
152 | } | |
153 | if(!found) { | |
154 | max++; | |
155 | idx[max - 1] = lab; | |
156 | count[max - 1] = 1; | |
157 | } | |
158 | } | |
159 | } | |
160 | ||
161 | j = 0, max = count[0]; | |
162 | for (i = 0; i < 18; i++) { | |
163 | if (count[i] > max) { | |
164 | j = i; | |
165 | max = count[i]; | |
166 | } | |
167 | } | |
168 | fLabel = idx[j]; | |
169 | fCount = count[j]; | |
170 | } | |
171 | // | |
172 | // | |
173 | // | |
174 | void AliITSNeuralTrack::CleanSlot(Int_t i, Bool_t del) | |
853a0f19 | 175 | { |
b9d722bc | 176 | // Removes a point from the corresponding layer slot in the found track. |
177 | // If the argument is TRUE, the point object is also deleted from heap. | |
853a0f19 | 178 | |
b9d722bc | 179 | if (i >= 0 && i < 6) { |
180 | if (del) delete fPoint[i]; | |
181 | fPoint[i] = 0; | |
182 | } | |
183 | } | |
184 | // | |
185 | // | |
186 | // | |
187 | void AliITSNeuralTrack::GetModuleData(Int_t layer, Int_t &mod, Int_t &pos) | |
853a0f19 | 188 | { |
b9d722bc | 189 | // Returns the point coordinates according to the TreeR philosophy in galice.root files |
190 | // that consist in the module number (mod) and the position in the TClonesArray of | |
191 | // the points reconstructed in that module for the run being examined. | |
853a0f19 | 192 | |
b9d722bc | 193 | if (layer < 0 || layer > 5) { |
194 | Error("GetModuleData", "Layer out of range: %d", layer); | |
195 | return; | |
196 | } | |
197 | mod = fPoint[layer]->GetModule(); | |
198 | pos = fPoint[layer]->GetIndex(); | |
199 | } | |
200 | // | |
201 | // | |
202 | // | |
203 | void AliITSNeuralTrack::Insert(AliITSNeuralPoint *point) | |
853a0f19 | 204 | { |
b9d722bc | 205 | // A trivial method to insert a point in the tracks; |
206 | // the point is inserted to the slot corresponding to its ITS layer. | |
853a0f19 | 207 | |
b9d722bc | 208 | if (!point) return; |
209 | ||
210 | Int_t layer = point->GetLayer(); | |
211 | if (layer < 0 || layer > 6) { | |
212 | Error("Insert", "Layer index %d out of range", layer); | |
213 | return; | |
214 | } | |
215 | ||
216 | fPoint[layer] = point; | |
217 | } | |
218 | // | |
219 | // | |
220 | // | |
221 | Int_t AliITSNeuralTrack::OccupationMask() | |
853a0f19 | 222 | { |
b9d722bc | 223 | // Returns a byte which maps the occupied slots. |
224 | // Each bit represents a layer going from the less significant on. | |
853a0f19 | 225 | |
b9d722bc | 226 | Int_t i, check, mask = 0; |
227 | for (i = 0; i < 6; i++) { | |
228 | check = 1 << i; | |
229 | if (fPoint[i]) mask |= check; | |
230 | } | |
231 | return mask; | |
232 | } | |
233 | // | |
234 | // | |
235 | // | |
236 | void AliITSNeuralTrack::PrintLabels() | |
853a0f19 | 237 | { |
b9d722bc | 238 | // Prints the results of the AssignLabel() method, together with |
239 | // the GEANT labels assigned to each point, in order to evaluate | |
240 | // how the assigned label is distributed among points. | |
853a0f19 | 241 | |
b9d722bc | 242 | cout << "Assigned label = " << fLabel << " -- counted " << fCount << " times: " << endl << endl; |
243 | for (Int_t i = 0; i < 6; i++) { | |
244 | cout << "Point #" << i + 1 << " --> "; | |
245 | if (fPoint[i]) { | |
246 | cout << "labels = " << fPoint[i]->GetLabel(0) << ", "; | |
247 | cout << fPoint[i]->GetLabel(1) << ", "; | |
248 | cout << fPoint[i]->GetLabel(2) << endl; | |
249 | } | |
250 | else { | |
251 | cout << "not assigned" << endl; | |
252 | } | |
253 | } | |
254 | cout << endl; | |
255 | } | |
256 | // | |
257 | // | |
258 | // | |
259 | Bool_t AliITSNeuralTrack::AddEL(Int_t layer, Double_t sign) | |
260 | { | |
853a0f19 | 261 | // Calculates the correction for energy loss |
262 | ||
b9d722bc | 263 | Double_t width = 0.0; |
264 | switch (layer) { | |
265 | case 0: width = 0.00260 + 0.00283; break; | |
266 | case 1: width = 0.0180; break; | |
267 | case 2: width = 0.0094; break; | |
268 | case 3: width = 0.0095; break; | |
269 | case 4: width = 0.0091; break; | |
270 | case 5: width = 0.0087; break; | |
271 | default: | |
272 | Error("AddEL", "Layer value %d out of range!", layer); | |
273 | return kFALSE; | |
274 | } | |
275 | width *= 1.7; | |
276 | ||
277 | if((layer == 5) && (fStatePhi < 0.174 || fStatePhi > 6.100 || (fStatePhi > 2.960 && fStatePhi < 3.31))) { | |
278 | width += 0.012; | |
279 | } | |
280 | ||
281 | Double_t invSqCosL = 1. + fTanL * fTanL; // = 1 / (cos(lambda)^2) = 1 + tan(lambda)^2 | |
282 | Double_t invCosL = TMath::Sqrt(invSqCosL); // = 1 / cos(lambda) | |
283 | Double_t pt = GetPt(); // = transverse momentum | |
284 | Double_t p2 = pt *pt * invSqCosL; // = square modulus of momentum | |
285 | Double_t energy = TMath::Sqrt(p2 + fMass * fMass); // = energy | |
286 | Double_t beta2 = p2 / (p2 + fMass * fMass); // = (v / c) ^ 2 | |
287 | if (beta2 == 0.0) { | |
288 | printf("Anomaly in AddEL: pt=%8.6f invSqCosL=%8.6f fMass=%8.7f --> beta2 = %8.7f\n", pt, invSqCosL, fMass, beta2); | |
289 | return kFALSE; | |
290 | } | |
291 | ||
292 | Double_t dE = 0.153 / beta2 * (log(5940. * beta2 / (1. - beta2)) - beta2) * width * 21.82 * invCosL; | |
293 | dE = sign * dE * 0.001; | |
294 | ||
295 | energy += dE; | |
296 | p2 = energy * energy - fMass * fMass; | |
297 | pt = TMath::Sqrt(p2) / invCosL; | |
298 | if (fC < 0.) pt = -pt; | |
299 | fC = (0.299792458 * 0.2 * fField) / (pt * 100.); | |
300 | ||
301 | return kTRUE; | |
302 | } | |
303 | // | |
304 | // | |
305 | // | |
306 | Bool_t AliITSNeuralTrack::AddMS(Int_t layer) | |
307 | { | |
853a0f19 | 308 | // Calculates the noise perturbation due to multiple scattering |
309 | ||
b9d722bc | 310 | Double_t width = 0.0; |
311 | switch (layer) { | |
312 | case 0: width = 0.00260 + 0.00283; break; | |
313 | case 1: width = 0.0180; break; | |
314 | case 2: width = 0.0094; break; | |
315 | case 3: width = 0.0095; break; | |
316 | case 4: width = 0.0091; break; | |
317 | case 5: width = 0.0087; break; | |
318 | default: | |
319 | Error("AddEL", "Layer value %d out of range!", layer); | |
320 | return kFALSE; | |
321 | } | |
322 | width *= 1.7; | |
323 | ||
324 | if((layer == 5) && (fStatePhi < 0.174 || fStatePhi > 6.100 || (fStatePhi > 2.960 && fStatePhi < 3.31))) { | |
325 | width += 0.012; | |
326 | } | |
327 | ||
328 | Double_t cosL = TMath::Cos(TMath::ATan(fTanL)); | |
329 | Double_t halfC = fC / 2.; | |
330 | Double_t q20 = 1. / (cosL * cosL); | |
331 | Double_t q30 = fC * fTanL; | |
332 | ||
333 | Double_t q40 = halfC * (fStateR * fStateR - fDt * fDt) / (1. + 2. * halfC * fDt); | |
334 | Double_t dd = fDt + halfC * fDt * fDt - halfC * fStateR * fStateR; | |
335 | Double_t dprova = fStateR * fStateR - dd * dd; | |
336 | Double_t q41 = 0.; | |
337 | if(dprova > 0.) q41 = -1. / cosL * TMath::Sqrt(dprova) / (1. + 2. * halfC *fDt); | |
338 | ||
339 | Double_t p2 = (GetPt()*GetPt()) / (cosL * cosL); | |
340 | Double_t beta2 = p2 / (p2 + fMass * fMass); | |
341 | Double_t theta2 = 14.1 * 14.1 / (beta2 * p2 * 1.e6) * (width / TMath::Abs(cosL)); | |
342 | ||
343 | fMatrix(2,2) += theta2 * (q40 * q40 + q41 * q41); | |
344 | fMatrix(3,2) += theta2 * q20 * q40; | |
345 | fMatrix(2,3) += theta2 * q20 * q40; | |
346 | fMatrix(3,3) += theta2 * q20 * q20; | |
347 | fMatrix(4,2) += theta2 * q30 * q40; | |
348 | fMatrix(2,4) += theta2 * q30 * q40; | |
349 | fMatrix(4,3) += theta2 * q30 * q20; | |
350 | fMatrix(3,4) += theta2 * q30 * q20; | |
351 | fMatrix(4,4) += theta2 * q30 * q30; | |
352 | ||
353 | return kTRUE; | |
354 | } | |
355 | // | |
356 | // | |
357 | // | |
358 | Int_t AliITSNeuralTrack::PropagateTo(Double_t rk) | |
359 | { | |
360 | // Propagation method. | |
361 | // Changes the state vector according to a new radial position | |
362 | // which is specified by the passed 'r' value (in cylindircal coordinates). | |
363 | // The covariance matrix is also propagated (and enlarged) according to | |
364 | // the FCFt technique, where F is the jacobian of the new parameters | |
365 | // w.r.t. their old values. | |
366 | // The option argument forces the method to add also the energy loss | |
367 | // and the multiple scattering effects, which respectively have the effect | |
368 | // of changing the curvature and widening the covariance matrix. | |
369 | ||
b9f05b32 | 370 | if (rk < TMath::Abs(fDt)) { |
b9d722bc | 371 | Error("PropagateTo", Form("Impossible propagation to r (=%17.15g) < Dt (=%17.15g)", rk, fDt)); |
372 | return 0; | |
373 | } | |
374 | ||
375 | Double_t duepi = 2. * TMath::Pi(); | |
376 | Double_t rkm1 = fStateR; | |
377 | Double_t aAk = ArgPhi(rk), aAkm1 = ArgPhi(rkm1); | |
378 | Double_t ak = ArgZ(rk), akm1 = ArgZ(rkm1); | |
379 | ||
380 | fStatePhi += TMath::ASin(aAk) - TMath::ASin(aAkm1); | |
381 | if(fStatePhi > duepi) fStatePhi -= duepi; | |
382 | if(fStatePhi < 0.) fStatePhi += duepi; | |
383 | ||
384 | Double_t halfC = 0.5 * fC; | |
385 | fStateZ += fTanL / halfC * (TMath::ASin(ak)-TMath::ASin(akm1)); | |
386 | ||
387 | Double_t bk = ArgB(rk), bkm1 = ArgB(rkm1); | |
388 | Double_t ck = ArgC(rk), ckm1 = ArgC(rkm1); | |
389 | ||
390 | Double_t f02 = ck / TMath::Sqrt(1. - aAk * aAk) - ckm1 / TMath::Sqrt(1. - aAkm1 * aAkm1); | |
391 | Double_t f04 = bk / TMath::Sqrt(1. - aAk * aAk) - bkm1 / TMath::Sqrt(1. - aAkm1 * aAkm1); | |
392 | Double_t f12 = fTanL * fDt * (1. / rk - 1. / rkm1); | |
393 | Double_t f13 = rk - rkm1; | |
394 | ||
395 | Double_t c00 = fMatrix(0,0); | |
396 | Double_t c10 = fMatrix(1,0); | |
397 | Double_t c11 = fMatrix(1,1); | |
398 | Double_t c20 = fMatrix(2,0); | |
399 | Double_t c21 = fMatrix(2,1); | |
400 | Double_t c22 = fMatrix(2,2); | |
401 | Double_t c30 = fMatrix(3,0); | |
402 | Double_t c31 = fMatrix(3,1); | |
403 | Double_t c32 = fMatrix(3,2); | |
404 | Double_t c33 = fMatrix(3,3); | |
405 | Double_t c40 = fMatrix(4,0); | |
406 | Double_t c41 = fMatrix(4,1); | |
407 | Double_t c42 = fMatrix(4,2); | |
408 | Double_t c43 = fMatrix(4,3); | |
409 | Double_t c44 = fMatrix(4,4); | |
410 | ||
411 | Double_t r10 = c10 + c21*f02 + c41*f04; | |
412 | Double_t r20 = c20 + c22*f02 + c42*f04; | |
413 | Double_t r30 = c30 + c32*f02 + c43*f04; | |
414 | Double_t r40 = c40 + c42*f02 + c44*f04; | |
415 | Double_t r21 = c21 + c22*f12 + c32*f13; | |
416 | Double_t r31 = c31 + c32*f12 + c33*f13; | |
417 | Double_t r41 = c41 + c42*f12 + c43*f13; | |
418 | ||
419 | fMatrix(0,0) = c00 + c20*f02 + c40*f04 + f02*r20 + f04*r40; | |
420 | fMatrix(1,0) = fMatrix(0,1) = r10 + f12*r20 + f13*r30; | |
421 | fMatrix(1,1) = c11 + c21*f12 + c31*f13 + f12*r21 + f13*r31; | |
422 | fMatrix(2,0) = fMatrix(0,2) = r20; | |
423 | fMatrix(2,1) = fMatrix(1,2) = r21; | |
424 | fMatrix(3,0) = fMatrix(0,3) = r30; | |
425 | fMatrix(3,1) = fMatrix(1,3) = r31; | |
426 | fMatrix(4,0) = fMatrix(0,4) = r40; | |
427 | fMatrix(4,1) = fMatrix(1,4) = r41; | |
428 | ||
429 | fStateR = rk; | |
430 | ||
431 | if (rkm1 < fStateR) // going to greater R --> energy LOSS | |
432 | return -1; | |
433 | else // going to smaller R --> energy GAIN | |
434 | return 1; | |
435 | } | |
436 | // | |
437 | // | |
438 | // | |
439 | Bool_t AliITSNeuralTrack::SeedCovariance() | |
440 | { | |
441 | // generate a covariance matrix depending on the results obtained from | |
442 | // the preliminary seeding fit procedure. | |
443 | // It calculates the variances for C, D ans TanL, according to the | |
444 | // differences of the fitted values from the requested ones necessary | |
445 | // to make the curve exactly pass through each point. | |
446 | ||
447 | /* | |
448 | Int_t i, j; | |
449 | AliITSNeuralPoint *p = 0; | |
450 | Double_t r, argPhi, phiC, phiD, argZ, zL; | |
451 | Double_t sumC = 0.0, sumD = 0.0, sumphi = 0., sumz = 0., sumL = 0.; | |
452 | for (i = 0; i < fNum; i++) { | |
453 | p = At(i); | |
454 | if (!p) continue; | |
455 | r = p->GetR2(); | |
456 | // weight and derivatives of phi and zeta w.r.t. various params | |
457 | sumphi += 1./ p->ErrorGetPhi(); | |
458 | argPhi = ArgPhi(r); | |
459 | argZ = ArgZ(r); | |
460 | if (argPhi > 100.0 || argZ > 100.0) { | |
461 | Error("InitCovariance", "Argument error"); | |
462 | return kFALSE; | |
463 | } | |
464 | phiC = DerArgPhiC(r) / TMath::Sqrt(1.0 - argPhi * argPhi); | |
465 | phiD = DerArgPhiD(r) / TMath::Sqrt(1.0 - argPhi * argPhi); | |
466 | if (phiC > 100.0 || phiD > 100.0) { | |
467 | Error("InitCovariance", "Argument error"); | |
468 | return kFALSE; | |
469 | } | |
470 | zL = asin(argZ) / fC; | |
471 | sumL += zL * zL; | |
472 | sumC += phiC * phiC; | |
473 | sumD += phiD * phiD; | |
474 | sumz += 1.0 / (p->fError[2] * p->fError[2]); | |
475 | } | |
476 | ||
477 | for (i = 0; i < 5; i++) for (j = 0; j < 5; j++) fMatrix(i,j) = 0.; | |
478 | fMatrix(0,0) = 1. / sumphi; | |
479 | fMatrix(1,1) = 1. / sumz; | |
480 | fMatrix(2,2) = 1. / sumD; | |
481 | fMatrix(3,3) = 1. / sumL; | |
482 | fMatrix(4,4) = 1. / sumC; | |
483 | fMatrix.Print(); | |
484 | */ | |
485 | ||
486 | AliITSNeuralPoint *p = 0; | |
487 | Double_t delta, cs, sn, r, argz; | |
488 | Double_t diffC, diffD, diffL, calcC, calcD, calcL; | |
489 | ||
490 | Int_t l; | |
491 | for (l = 0; l < 6; l++) { | |
492 | p = fPoint[l]; | |
493 | if (!p) break; | |
494 | sn = TMath::Sin(p->GetPhi() - fG0); | |
495 | cs = TMath::Cos(p->GetPhi() - fG0); | |
496 | r = p->GetR2(); | |
497 | calcC = (fDt/r - sn) / (2.*fDt*sn - r - fDt*fDt/r); | |
498 | argz = ArgZ(r); | |
499 | if (argz > 1000.0) { | |
500 | Error("Covariance", "Value too high"); | |
501 | return kFALSE; | |
502 | } | |
503 | calcL = (p->Z() - fDz) * fC / asin(argz); | |
504 | delta = fR*fR + r*r + 2.0*fR*r*sin(p->GetPhi() - fG0); | |
505 | if (delta < 0.E0) { | |
506 | if (delta >= -0.5) | |
507 | delta = 0.; | |
508 | else { | |
509 | Error("Covariance", Form("Discriminant = %g --- Dt = %g", delta, fDt)); | |
510 | return kFALSE; | |
511 | } | |
512 | } | |
513 | delta = sqrt(delta); | |
514 | if (fC >= 0) | |
515 | calcD = delta - fR; | |
516 | else | |
517 | calcD = fR - delta; | |
518 | diffD = calcD - fDt; | |
519 | diffL = calcL - fTanL; | |
520 | diffC = fC - calcC; | |
521 | fMatrix(0,0) += 100000000.0 * p->GetError("phi") * p->GetError("phi"); | |
522 | fMatrix(1,1) += 10000.0 * p->ErrZ() * p->ErrZ(); | |
523 | fMatrix(2,2) += 100000.0 * diffD * diffD; | |
524 | fMatrix(3,3) += diffL * diffL; | |
525 | fMatrix(4,4) += 100000000.0 * diffC * diffC; | |
526 | } | |
853a0f19 | 527 | Double_t n = 0.; |
528 | for (l = 0; l < 6; l++) if (fPoint[l]) n++; | |
f77f13c8 | 529 | fMatrix *= 1./(n * (n+1)); |
b9d722bc | 530 | return kTRUE; |
531 | } | |
532 | // | |
533 | // | |
534 | // | |
535 | Bool_t AliITSNeuralTrack::Filter(AliITSNeuralPoint *test) | |
536 | { | |
537 | // Makes all calculations which apply the Kalman filter to the | |
538 | // stored guess of the state vector, after propagation to a new layer | |
539 | ||
540 | if (!test) { | |
541 | Error("Filter", "Null pointer passed"); | |
542 | return kFALSE; | |
543 | } | |
544 | ||
545 | Double_t m[2]; | |
546 | Double_t rk, phik, zk; | |
547 | rk = test->GetR2(); | |
548 | phik = test->GetPhi(); | |
549 | zk = test->Z(); | |
550 | m[0]=phik; | |
551 | m[1]=zk; | |
552 | ||
553 | //////////////////////// Evaluation of the error matrix V ///////// | |
554 | Double_t v00 = test->GetError("phi") * rk; | |
555 | Double_t v11 = test->ErrZ(); | |
556 | //////////////////////////////////////////////////////////////////// | |
557 | ||
558 | // Get the covariance matrix | |
559 | Double_t cin00, cin10, cin20, cin30, cin40; | |
560 | Double_t cin11, cin21, cin31, cin41, cin22; | |
561 | Double_t cin32, cin42, cin33, cin43, cin44; | |
562 | cin00 = fMatrix(0,0); | |
563 | cin10 = fMatrix(1,0); | |
564 | cin20 = fMatrix(2,0); | |
565 | cin30 = fMatrix(3,0); | |
566 | cin40 = fMatrix(4,0); | |
567 | cin11 = fMatrix(1,1); | |
568 | cin21 = fMatrix(2,1); | |
569 | cin31 = fMatrix(3,1); | |
570 | cin41 = fMatrix(4,1); | |
571 | cin22 = fMatrix(2,2); | |
572 | cin32 = fMatrix(3,2); | |
573 | cin42 = fMatrix(4,2); | |
574 | cin33 = fMatrix(3,3); | |
575 | cin43 = fMatrix(4,3); | |
576 | cin44 = fMatrix(4,4); | |
577 | ||
578 | // Calculate R matrix | |
579 | Double_t rold00 = cin00 + v00; | |
580 | Double_t rold10 = cin10; | |
581 | Double_t rold11 = cin11 + v11; | |
582 | ||
583 | ////////////////////// R matrix inversion ///////////////////////// | |
584 | Double_t det = rold00*rold11 - rold10*rold10; | |
585 | Double_t r00 = rold11/det; | |
586 | Double_t r10 = -rold10/det; | |
587 | Double_t r11 = rold00/det; | |
588 | //////////////////////////////////////////////////////////////////// | |
589 | ||
590 | // Calculate Kalman matrix | |
591 | Double_t k00 = cin00*r00 + cin10*r10; | |
592 | Double_t k01 = cin00*r10 + cin10*r11; | |
593 | Double_t k10 = cin10*r00 + cin11*r10; | |
594 | Double_t k11 = cin10*r10 + cin11*r11; | |
595 | Double_t k20 = cin20*r00 + cin21*r10; | |
596 | Double_t k21 = cin20*r10 + cin21*r11; | |
597 | Double_t k30 = cin30*r00 + cin31*r10; | |
598 | Double_t k31 = cin30*r10 + cin31*r11; | |
599 | Double_t k40 = cin40*r00 + cin41*r10; | |
600 | Double_t k41 = cin40*r10 + cin41*r11; | |
601 | ||
602 | // Get state vector (will keep the old values for phi and z) | |
603 | Double_t x0, x1, x2, x3, x4, savex0, savex1; | |
604 | x0 = savex0 = fStatePhi; | |
605 | x1 = savex1 = fStateZ; | |
606 | x2 = fDt; | |
607 | x3 = fTanL; | |
608 | x4 = fC; | |
609 | ||
610 | // Update the state vector | |
611 | x0 += k00*(m[0]-savex0) + k01*(m[1]-savex1); | |
612 | x1 += k10*(m[0]-savex0) + k11*(m[1]-savex1); | |
613 | x2 += k20*(m[0]-savex0) + k21*(m[1]-savex1); | |
614 | x3 += k30*(m[0]-savex0) + k31*(m[1]-savex1); | |
615 | x4 += k40*(m[0]-savex0) + k41*(m[1]-savex1); | |
616 | ||
617 | // Update the covariance matrix | |
618 | Double_t cout00, cout10, cout20, cout30, cout40; | |
619 | Double_t cout11, cout21, cout31, cout41, cout22; | |
620 | Double_t cout32, cout42, cout33, cout43, cout44; | |
621 | ||
622 | cout00 = cin00 - k00*cin00 - k01*cin10; | |
623 | cout10 = cin10 - k00*cin10 - k01*cin11; | |
624 | cout20 = cin20 - k00*cin20 - k01*cin21; | |
625 | cout30 = cin30 - k00*cin30 - k01*cin31; | |
626 | cout40 = cin40 - k00*cin40 - k01*cin41; | |
627 | cout11 = cin11 - k10*cin10 - k11*cin11; | |
628 | cout21 = cin21 - k10*cin20 - k11*cin21; | |
629 | cout31 = cin31 - k10*cin30 - k11*cin31; | |
630 | cout41 = cin41 - k10*cin40 - k11*cin41; | |
631 | cout22 = cin22 - k20*cin20 - k21*cin21; | |
632 | cout32 = cin32 - k20*cin30 - k21*cin31; | |
633 | cout42 = cin42 - k20*cin40 - k21*cin41; | |
634 | cout33 = cin33 - k30*cin30 - k31*cin31; | |
635 | cout43 = cin43 - k30*cin40 - k31*cin41; | |
636 | cout44 = cin44 - k40*cin40 - k41*cin41; | |
637 | ||
638 | // Store the new covariance matrix | |
639 | fMatrix(0,0) = cout00; | |
640 | fMatrix(1,0) = fMatrix(0,1) = cout10; | |
641 | fMatrix(2,0) = fMatrix(0,2) = cout20; | |
642 | fMatrix(3,0) = fMatrix(0,3) = cout30; | |
643 | fMatrix(4,0) = fMatrix(0,4) = cout40; | |
644 | fMatrix(1,1) = cout11; | |
645 | fMatrix(2,1) = fMatrix(1,2) = cout21; | |
646 | fMatrix(3,1) = fMatrix(1,3) = cout31; | |
647 | fMatrix(4,1) = fMatrix(1,4) = cout41; | |
648 | fMatrix(2,2) = cout22; | |
649 | fMatrix(3,2) = fMatrix(2,3) = cout32; | |
650 | fMatrix(4,2) = fMatrix(2,4) = cout42; | |
651 | fMatrix(3,3) = cout33; | |
652 | fMatrix(4,3) = fMatrix(3,4) = cout43; | |
653 | fMatrix(4,4) = cout44; | |
654 | ||
655 | // Calculation of the chi2 increment | |
656 | Double_t vmcold00 = v00 - cout00; | |
657 | Double_t vmcold10 = -cout10; | |
658 | Double_t vmcold11 = v11 - cout11; | |
659 | ////////////////////// Matrix vmc inversion /////////////////////// | |
660 | det = vmcold00*vmcold11 - vmcold10*vmcold10; | |
661 | Double_t vmc00=vmcold11/det; | |
662 | Double_t vmc10 = -vmcold10/det; | |
663 | Double_t vmc11 = vmcold00/det; | |
664 | //////////////////////////////////////////////////////////////////// | |
665 | Double_t chi2 = (m[0] - x0)*( vmc00*(m[0] - x0) + 2.*vmc10*(m[1] - x1) ) + (m[1] - x1)*vmc11*(m[1] - x1); | |
666 | fChi2 += chi2; | |
667 | fNSteps++; | |
668 | ||
669 | return kTRUE; | |
670 | } | |
671 | // | |
672 | // | |
673 | // | |
674 | Bool_t AliITSNeuralTrack::KalmanFit() | |
853a0f19 | 675 | { |
b9d722bc | 676 | // Applies the Kalman Filter to improve the track parameters resolution. |
677 | // First, thre point which lies closer to the estimated helix is chosen. | |
678 | // Then, a fit is performed towards the 6th layer | |
679 | // Finally, the track is refitted to the 1st layer | |
853a0f19 | 680 | |
b9d722bc | 681 | Double_t rho; |
682 | Int_t l, layer, sign; | |
683 | ||
684 | fStateR = fPoint[0]->GetR2(); | |
685 | fStatePhi = fPoint[0]->GetPhi(); | |
686 | fStateZ = fPoint[0]->Z(); | |
687 | ||
688 | if (!PropagateTo(3.0)) { | |
689 | Error("KalmanFit", "Unsuccessful initialization"); | |
690 | return kFALSE; | |
691 | } | |
692 | l=0; | |
693 | ||
694 | // Performs a Kalman filter going from the actual state position | |
695 | // towards layer 6 position | |
696 | // Now, the propagation + filtering operations can be performed | |
697 | Double_t argPhi = 0.0, argZ = 0.0; | |
698 | while (l <= 5) { | |
699 | if (!fPoint[l]) { | |
700 | Error("KalmanFit", "Not six points!"); | |
701 | return kFALSE; | |
702 | } | |
703 | layer = fPoint[l]->GetLayer(); | |
704 | rho = fPoint[l]->GetR2(); | |
705 | sign = PropagateTo(rho); | |
706 | if (!sign) return kFALSE; | |
707 | AddEL(layer, -1.0); | |
708 | AddMS(layer); | |
709 | if (!Filter(fPoint[l])) return kFALSE; | |
710 | // these two parameters are update according to the filtered values | |
711 | argPhi = ArgPhi(fStateR); | |
712 | argZ = ArgZ(fStateR); | |
713 | if (argPhi > 1.0 || argPhi < -1.0 || argZ > 1.0 || argZ < -1.0) { | |
714 | Error("Filter", Form("Filtering returns too large values: %g, %g", argPhi, argZ)); | |
715 | return kFALSE; | |
716 | } | |
717 | fG0 = fStatePhi - asin(argPhi); | |
718 | fDz = fStateZ - (2.0 * fTanL / fC) * asin(argZ); | |
719 | l++; | |
720 | } | |
721 | ||
722 | // Now a Kalman filter i performed going from the actual state position | |
723 | // towards layer 1 position and then propagates to vertex | |
724 | if (l >= 5) l = 5; | |
725 | while (l >= 1) { | |
726 | layer = fPoint[l]->GetLayer(); | |
727 | rho = fPoint[l]->GetR2(); | |
728 | AddEL(layer, 1.0); | |
729 | sign = PropagateTo(rho); | |
730 | if (!sign) return kFALSE; | |
731 | AddMS(layer); | |
732 | if (!Filter(fPoint[l])) return kFALSE; | |
733 | // these two parameters are update according to the filtered values | |
734 | argPhi = ArgPhi(fStateR); | |
735 | argZ = ArgZ(fStateR); | |
736 | if (argPhi > 1.0 || argPhi < -1.0 || argZ > 1.0 || argZ < -1.0) { | |
737 | Error("Filter", Form("Filtering returns too large values: %g, %g", argPhi, argZ)); | |
738 | return kFALSE; | |
739 | } | |
740 | fG0 = fStatePhi - asin(argPhi); | |
741 | fDz = fStateZ - (2.0 * fTanL / fC) * asin(argZ); | |
742 | l--; | |
743 | } | |
744 | return kTRUE; | |
745 | } | |
746 | // | |
747 | // | |
748 | // | |
749 | Bool_t AliITSNeuralTrack::RiemannFit() | |
750 | { | |
751 | // Method which executes the circle fit via a Riemann Sphere projection | |
752 | // with the only improvement of a weighted mean, due to different errors | |
753 | // over different point measurements. | |
754 | // As an output, it returns kTRUE or kFALSE respectively if the fit succeeded or not | |
755 | // in fact, if some variables assume strange values, the fit is aborted, | |
756 | // in order to prevent the class from raising a floating point error; | |
757 | ||
758 | Int_t i, j; | |
759 | ||
760 | // M1 - matrix of ones | |
761 | TMatrixD m1(6,1); | |
762 | for (i = 0; i < 6; i++) m1(i,0) = 1.0; | |
763 | ||
764 | // X - matrix of Rieman projection coordinates | |
853a0f19 | 765 | TMatrixD mX(6,3); |
b9d722bc | 766 | for (i = 0; i < 6; i++) { |
853a0f19 | 767 | mX(i,0) = fPoint[i]->X(); |
768 | mX(i,1) = fPoint[i]->Y(); | |
769 | mX(i,2) = fPoint[i]->GetR2sq(); | |
b9d722bc | 770 | } |
771 | ||
772 | // W - matrix of weights | |
773 | Double_t xterm, yterm, ex, ey; | |
853a0f19 | 774 | TMatrixD mW(6,6); |
b9d722bc | 775 | for (i = 0; i < 6; i++) { |
776 | xterm = fPoint[i]->X() * fPoint[i]->GetPhi() - fPoint[i]->Y() / fPoint[i]->GetR2(); | |
777 | ex = fPoint[i]->ErrX(); | |
778 | yterm = fPoint[i]->Y() * fPoint[i]->GetPhi() + fPoint[i]->X() / fPoint[i]->GetR2(); | |
779 | ey = fPoint[i]->ErrY(); | |
853a0f19 | 780 | mW(i,i) = fPoint[i]->GetR2sq() / (xterm * xterm * ex * ex + yterm * yterm * ey * ey); |
b9d722bc | 781 | } |
782 | ||
783 | // Xm - weighted sample mean | |
853a0f19 | 784 | Double_t meanX = 0.0, meanY = 0.0, meanW = 0.0, sw = 0.0; |
b9d722bc | 785 | for (i = 0; i < 6; i++) { |
853a0f19 | 786 | meanX += mW(i,i) * mX(i,0); |
787 | meanY += mW(i,i) * mX(i,1); | |
788 | meanW += mW(i,i) * mX(i,2); | |
789 | sw += mW(i,i); | |
b9d722bc | 790 | } |
853a0f19 | 791 | meanX /= sw; |
792 | meanY /= sw; | |
793 | meanW /= sw; | |
b9d722bc | 794 | |
795 | // V - sample covariance matrix | |
796 | for (i = 0; i < 6; i++) { | |
853a0f19 | 797 | mX(i,0) -= meanX; |
798 | mX(i,1) -= meanY; | |
799 | mX(i,2) -= meanW; | |
b9d722bc | 800 | } |
853a0f19 | 801 | TMatrixD mXt(TMatrixD::kTransposed, mX); |
802 | TMatrixD mWX(mW, TMatrixD::kMult, mX); | |
803 | TMatrixD mV(mXt, TMatrixD::kMult, mWX); | |
b9d722bc | 804 | for (i = 0; i < 3; i++) { |
805 | for (j = i + 1; j < 3; j++) { | |
853a0f19 | 806 | mV(i,j) = mV(j,i) = (mV(i,j) + mV(j,i)) * 0.5; |
b9d722bc | 807 | } |
808 | } | |
809 | ||
810 | // Eigenvalue problem solving for V matrix | |
811 | Int_t ileast = 0; | |
853a0f19 | 812 | TVectorD eval(3), n(3); |
49370114 | 813 | // TMatrixD evec = mV.EigenVectors(eval); |
814 | #if ROOT_VERSION_CODE >= 262146 | |
815 | TMatrixDEigen ei(mV); | |
816 | TMatrixD evec = ei.GetEigenVectors(); | |
817 | eval = ei.GetEigenValues(); | |
818 | #else | |
853a0f19 | 819 | TMatrixD evec = mV.EigenVectors(eval); |
49370114 | 820 | #endif |
821 | ||
853a0f19 | 822 | if (eval(1) < eval(ileast)) ileast = 1; |
823 | if (eval(2) < eval(ileast)) ileast = 2; | |
824 | n(0) = evec(0, ileast); | |
825 | n(1) = evec(1, ileast); | |
826 | n(2) = evec(2, ileast); | |
b9d722bc | 827 | |
828 | // c - known term in the plane intersection with Riemann axes | |
853a0f19 | 829 | Double_t c = -(meanX * n(0) + meanY * n(1) + meanW * n(2)); |
b9d722bc | 830 | |
831 | fXC = -n(0) / (2. * n(2)); | |
832 | fYC = -n(1) / (2. * n(2)); | |
833 | fR = (1. - n(2)*n(2) - 4.*c*n(2)) / (4. * n(2) * n(2)); | |
834 | ||
835 | if (fR <= 0.E0) { | |
836 | Error("RiemannFit", "Radius comed less than zero!!!"); | |
837 | return kFALSE; | |
838 | } | |
839 | fR = TMath::Sqrt(fR); | |
840 | fC = 1.0 / fR; | |
841 | ||
842 | // evaluating signs for curvature and others | |
843 | Double_t phi1 = 0.0, phi2, temp1, temp2, sumdphi = 0.0, ref = TMath::Pi(); | |
844 | AliITSNeuralPoint *p = fPoint[0]; | |
845 | phi1 = p->GetPhi(); | |
846 | for (i = 1; i < 6; i++) { | |
847 | p = (AliITSNeuralPoint*)fPoint[i]; | |
848 | if (!p) break; | |
849 | phi2 = p->GetPhi(); | |
850 | temp1 = phi1; | |
851 | temp2 = phi2; | |
852 | if (temp1 > ref && temp2 < ref) | |
853 | temp2 += 2.0 * ref; | |
854 | else if (temp1 < ref && temp2 > ref) | |
855 | temp1 += 2.0 * ref; | |
856 | sumdphi += temp2 - temp1; | |
857 | phi1 = phi2; | |
858 | } | |
859 | if (sumdphi < 0.E0) fC = -fC; | |
860 | Double_t diff, angle = TMath::ATan2(fYC, fXC); | |
861 | if (fC < 0.E0) | |
862 | fG0 = angle + 0.5 * TMath::Pi(); | |
863 | else | |
864 | fG0 = angle - 0.5 * TMath::Pi(); | |
865 | diff = angle - fG0; | |
866 | ||
853a0f19 | 867 | Double_t d = TMath::Sqrt(fXC*fXC + fYC*fYC) - fR; |
b9d722bc | 868 | if (fC >= 0.E0) |
853a0f19 | 869 | fDt = d; |
b9d722bc | 870 | else |
853a0f19 | 871 | fDt = -d; |
b9d722bc | 872 | |
853a0f19 | 873 | Int_t nn = 6; |
b9d722bc | 874 | Double_t halfC = 0.5 * fC; |
853a0f19 | 875 | Double_t *s = new Double_t[nn], *z = new Double_t[nn], *ws = new Double_t[nn]; |
b9d722bc | 876 | for (j = 0; j < 6; j++) { |
877 | p = fPoint[j]; | |
878 | if (!p) break; | |
879 | s[j] = ArgZ(p->GetR2()); | |
880 | if (s[j] > 100.0) return kFALSE; | |
881 | z[j] = p->Z(); | |
882 | s[j] = asin(s[j]) / halfC; | |
883 | ws[j] = 1.0 / (p->ErrZ()*p->ErrZ()); | |
884 | } | |
885 | ||
886 | // second tep final fit | |
853a0f19 | 887 | Double_t sums2 = 0.0, sumz = 0.0, sumsz = 0.0, sums = 0.0, sumw = 0.0; |
888 | for (i = 0; i < nn; i++) { | |
889 | sums2 += ws[i] * s[i] * s[i]; | |
890 | sumz += ws[i] * z[i]; | |
891 | sums += ws[i] * s[i]; | |
892 | sumsz += ws[i] * s[i] * z[i]; | |
b9d722bc | 893 | sumw += ws[i]; |
894 | } | |
853a0f19 | 895 | sums2 /= sumw; |
896 | sumz /= sumw; | |
897 | sums /= sumw; | |
898 | sumsz /= sumw; | |
899 | d = sums2 - sums*sums; | |
b9d722bc | 900 | |
853a0f19 | 901 | fDz = (sums2*sumz - sums*sumsz) / d; |
902 | fTanL = (sumsz - sums*sumz) / d; | |
b9d722bc | 903 | |
904 | delete [] s; | |
905 | delete [] z; | |
906 | delete [] ws; | |
907 | ||
908 | return kTRUE; | |
909 | } | |
910 | // | |
911 | // | |
912 | // | |
913 | void AliITSNeuralTrack::PrintState(Bool_t matrix) | |
853a0f19 | 914 | { |
b9d722bc | 915 | // Prints the state vector values. |
916 | // The argument switches on or off the printing of the covariance matrix. | |
853a0f19 | 917 | |
b9d722bc | 918 | cout << "\nState vector: " << endl; |
919 | cout << " Rho = " << fStateR << "\n"; | |
920 | cout << " Phi = " << fStatePhi << "\n"; | |
921 | cout << " Z = " << fStateZ << "\n"; | |
922 | cout << " Dt = " << fDt << "\n"; | |
923 | cout << " Dz = " << fDz << "\n"; | |
924 | cout << "TanL = " << fTanL << "\n"; | |
925 | cout << " C = " << fC << "\n"; | |
926 | cout << " G0 = " << fG0 << "\n"; | |
927 | cout << " XC = " << fXC << "\n"; | |
928 | cout << " YC = " << fYC << "\n"; | |
929 | if (matrix) { | |
930 | cout << "\nCovariance Matrix: " << endl; | |
931 | fMatrix.Print(); | |
932 | } | |
933 | cout << "Actual square chi = " << fChi2; | |
934 | } | |
935 | // | |
936 | // | |
937 | // | |
853a0f19 | 938 | Double_t AliITSNeuralTrack::GetDz() const |
b9d722bc | 939 | { |
940 | // Double_t argZ = ArgZ(fStateR); | |
941 | // if (argZ > 9.9) { | |
942 | // Error("GetDz", "Too large value: %g", argZ); | |
943 | // return 0.0; | |
944 | // } | |
945 | // fDz = fStateZ - (2.0 * fTanL / fC) * asin(argZ); | |
946 | return fDz; | |
947 | } | |
948 | // | |
949 | // | |
950 | // | |
853a0f19 | 951 | Double_t AliITSNeuralTrack::GetGamma() const |
b9d722bc | 952 | { |
953 | // these two parameters are update according to the filtered values | |
954 | // Double_t argPhi = ArgPhi(fStateR); | |
955 | // if (argPhi > 9.9) { | |
956 | // Error("Filter", "Too large value: %g", argPhi); | |
957 | // return kFALSE; | |
958 | // } | |
959 | // fG0 = fStatePhi - asin(argPhi); | |
960 | return fG0; | |
961 | } | |
962 | // | |
963 | // | |
964 | // | |
853a0f19 | 965 | Double_t AliITSNeuralTrack::GetPhi(Double_t r) const |
966 | { | |
b9d722bc | 967 | // Gives the value of azymuthal coordinate in the helix |
968 | // as a function of cylindric radius | |
853a0f19 | 969 | |
b9d722bc | 970 | Double_t arg = ArgPhi(r); |
971 | if (arg > 0.9) return 0.0; | |
972 | arg = fG0 + asin(arg); | |
973 | while (arg >= 2.0 * TMath::Pi()) { arg -= 2.0 * TMath::Pi(); } | |
974 | while (arg < 0.0) { arg += 2.0 * TMath::Pi(); } | |
975 | return arg; | |
976 | } | |
977 | // | |
978 | // | |
979 | // | |
853a0f19 | 980 | Double_t AliITSNeuralTrack::GetZ(Double_t r) const |
981 | { | |
b9d722bc | 982 | // gives the value of Z in the helix |
983 | // as a function of cylindric radius | |
853a0f19 | 984 | |
b9d722bc | 985 | Double_t arg = ArgZ(r); |
986 | if (arg > 0.9) return 0.0; | |
987 | return fDz + fTanL * asin(arg) / fC; | |
988 | } | |
989 | // | |
990 | // | |
991 | // | |
992 | Double_t AliITSNeuralTrack::GetdEdX() | |
993 | { | |
853a0f19 | 994 | // total energy loss of the track |
995 | ||
b9d722bc | 996 | Double_t q[4] = {0., 0., 0., 0.}, dedx = 0.0; |
997 | Int_t i = 0, swap = 0; | |
998 | for (i = 2; i < 6; i++) { | |
999 | if (!fPoint[i]) continue; | |
1000 | q[i - 2] = (Double_t)fPoint[i]->GetCharge(); | |
1001 | q[i - 2] /= (1 + fTanL*fTanL); | |
1002 | } | |
1003 | q[0] /= 280.; | |
1004 | q[1] /= 280.; | |
1005 | q[2] /= 38.; | |
1006 | q[3] /= 38.; | |
1007 | do { | |
1008 | swap = 0; | |
1009 | for (i = 0; i < 3; i++) { | |
1010 | if (q[i] <= q[i + 1]) continue; | |
1011 | Double_t tmp = q[i]; | |
1012 | q[i] = q[i + 1]; | |
1013 | q[i+1] = tmp; | |
1014 | swap++; | |
1015 | } | |
1016 | } while(swap); | |
1017 | if(q[0] < 0.) { | |
1018 | q[0] = q[1]; | |
1019 | q[1] = q[2]; | |
1020 | q[2] = q[3]; | |
1021 | q[3] = -1.; | |
1022 | } | |
1023 | dedx = (q[0] + q[1]) / 2.; | |
1024 | return dedx; | |
1025 | } | |
1026 | // | |
1027 | // | |
1028 | // | |
1029 | void AliITSNeuralTrack::SetVertex(Double_t *pos, Double_t *err) | |
1030 | { | |
1031 | // Stores vertex data | |
1032 | ||
1033 | if (!pos || !err) return; | |
1034 | fVertex.ErrX() = err[0]; | |
1035 | fVertex.ErrY() = err[1]; | |
1036 | fVertex.ErrZ() = err[2]; | |
1037 | fVertex.SetLayer(0); | |
1038 | fVertex.SetModule(0); | |
1039 | fVertex.SetIndex(0); | |
1040 | fVertex.SetLabel(0, -1); | |
1041 | fVertex.SetLabel(1, -1); | |
1042 | fVertex.SetLabel(2, -1); | |
1043 | fVertex.SetUser(1); | |
1044 | } | |
1045 | // | |
1046 | // | |
1047 | // | |
1048 | AliITSIOTrack* AliITSNeuralTrack::ExportIOtrack(Int_t min) | |
b9d722bc | 1049 | { |
853a0f19 | 1050 | // Exports an object in the standard format for reconstructed tracks |
1051 | ||
b9d722bc | 1052 | Int_t layer = 0; |
1053 | AliITSIOTrack *track = new AliITSIOTrack; | |
1054 | ||
1055 | // covariance matrix | |
1056 | track->SetCovMatrix(fMatrix(0,0), fMatrix(1,0), fMatrix(1,1), | |
1057 | fMatrix(2,0), fMatrix(2,1), fMatrix(2,2), | |
1058 | fMatrix(3,0), fMatrix(3,1), fMatrix(3,2), | |
1059 | fMatrix(3,3), fMatrix(4,0), fMatrix(4,1), | |
1060 | fMatrix(4,2), fMatrix(4,3), fMatrix(4,4)); | |
1061 | ||
1062 | // labels | |
1063 | track->SetLabel(IsGood(min) ? fLabel : -fLabel); | |
1064 | track->SetTPCLabel(-1); | |
1065 | ||
1066 | // points characteristics | |
1067 | for (layer = 0; layer < 6; layer++) { | |
1068 | if (fPoint[layer]) { | |
1069 | track->SetIdModule(layer, fPoint[layer]->GetModule()); | |
1070 | track->SetIdPoint(layer, fPoint[layer]->GetIndex()); | |
1071 | } | |
1072 | } | |
1073 | ||
1074 | // state vector | |
1075 | track->SetStatePhi(fStatePhi); | |
1076 | track->SetStateZ(fStateZ); | |
1077 | track->SetStateD(fDt); | |
1078 | track->SetStateTgl(fTanL); | |
1079 | track->SetStateC(fC); | |
1080 | track->SetRadius(fStateR); | |
1081 | track->SetCharge((fC > 0.0) ? -1 : 1); | |
1082 | track->SetDz(fDz); | |
1083 | ||
1084 | // track parameters in the closest point | |
1085 | track->SetX(fStateR * cos(fStatePhi)); | |
1086 | track->SetY(fStateR * cos(fStatePhi)); | |
1087 | track->SetZ(fStateZ); | |
1088 | track->SetPx(GetPt() * cos(fG0)); | |
1089 | track->SetPy(GetPt() * sin(fG0)); | |
1090 | track->SetPz(GetPt() * fTanL); | |
1091 | ||
1092 | // PID | |
1093 | track->SetPid(fPDG); | |
1094 | track->SetMass(fMass); | |
1095 | ||
1096 | return track; | |
1097 | } | |
1098 | // | |
1099 | // | |
1100 | //==================================================================================== | |
1101 | //============================ PRIVATE METHODS ============================ | |
1102 | //==================================================================================== | |
1103 | // | |
1104 | // | |
1105 | Double_t AliITSNeuralTrack::ArgPhi(Double_t r) const | |
1106 | { | |
1107 | // calculates the expression ((1/2)Cr + (1 + (1/2)CD) D/r) / (1 + CD) | |
1108 | ||
1109 | Double_t arg, num, den; | |
1110 | num = (0.5 * fC * r) + (1. + (0.5 * fC * fDt)) * (fDt / r); | |
1111 | den = 1. + fC * fDt; | |
1112 | if (den == 0.) { | |
1113 | Error("ArgPhi", "Denominator = 0!"); | |
1114 | return 10.0; | |
1115 | } | |
1116 | arg = num / den; | |
1117 | if (TMath::Abs(arg) < 1.) return arg; | |
1118 | if (TMath::Abs(arg) <= 1.00001) return (arg > 0.) ? 0.99999999999 : -0.9999999999; | |
1119 | Error("ArgPhi", "Value too large: %17.15g", arg); | |
1120 | return 10.0; | |
1121 | } | |
1122 | // | |
1123 | // | |
1124 | // | |
1125 | Double_t AliITSNeuralTrack::ArgZ(Double_t r) const | |
1126 | { | |
1127 | // calculates the expression (1/2)C * sqrt( (r^2 - Dt^2) / (1 + CD) ) | |
1128 | ||
1129 | Double_t arg; | |
1130 | arg = (r * r - fDt * fDt) / (1. + fC * fDt); | |
1131 | if (arg < 0.) { | |
b9f05b32 | 1132 | if (TMath::Abs(arg) < 1.E-6) arg = 0.; |
b9d722bc | 1133 | else { |
1134 | Error("ArgZ", "Square root argument error: %17.15g < 0", arg); | |
1135 | return 10.; | |
1136 | } | |
1137 | } | |
1138 | arg = 0.5 * fC * TMath::Sqrt(arg); | |
1139 | if (TMath::Abs(arg) < 1.) return arg; | |
1140 | if (TMath::Abs(arg) <= 1.00001) return (arg > 0.) ? 0.99999999999 : -0.9999999999; | |
1141 | Error("ArgZ", "Value too large: %17.15g", arg); | |
1142 | return 10.0; | |
1143 | } | |
1144 | // | |
1145 | // | |
1146 | // | |
1147 | Double_t AliITSNeuralTrack::ArgB(Double_t r) const | |
1148 | { | |
853a0f19 | 1149 | // UTILITY FUNCTION |
1150 | ||
b9d722bc | 1151 | Double_t arg; |
1152 | arg = (r*r - fDt*fDt); | |
1153 | arg /= (r*(1.+ fC*fDt)*(1.+ fC*fDt)); | |
1154 | return arg; | |
1155 | } | |
1156 | // | |
1157 | // | |
1158 | // | |
1159 | Double_t AliITSNeuralTrack::ArgC(Double_t r) const | |
1160 | { | |
853a0f19 | 1161 | // UTILITY FUNCTION |
1162 | ||
b9d722bc | 1163 | Double_t arg; |
1164 | arg = (1./r - fC * ArgPhi(r)); | |
1165 | arg /= 1.+ fC*fDt; | |
1166 | return arg; | |
1167 | } |