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