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