3 // The format of output data from Neural Tracker
4 // It can export data in the format of AliITSIOTrack
6 // Compatibility adaptation to V2 tracking is on the way.
7 // Author: A. Pulvirenti
13 //#include <TObject.h>
17 //#include <TObjArray.h>
20 #if ROOT_VERSION_CODE >= 262146
21 #include <TMatrixDEigen.h>
24 //#include "AliITSVertex.h"
25 #include "AliITSIOTrack.h"
26 #include "AliITSNeuralPoint.h"
28 #include "AliITSNeuralTrack.h"
32 ClassImp(AliITSNeuralTrack)
36 AliITSNeuralTrack::AliITSNeuralTrack() : fMatrix(5,5), fVertex()
38 // Default constructor
42 fMass = 0.1396; // default assumption: pion
43 fField = 2.0; // default assumption: B = 0.4 Tesla
45 fXC = fYC = fR = fC = 0.0;
46 fTanL = fG0 = fDt = fDz = 0.0;
47 fStateR = fStatePhi = fStateZ = fChi2 = fNSteps = 0.0;
51 for (i = 0; i < 6; i++) fPoint[i] = 0;
63 AliITSNeuralTrack::AliITSNeuralTrack(AliITSNeuralTrack &track)
64 : TObject((TObject&)track), fMatrix(5,5), fVertex()
70 fMass = 0.1396; // default assumption: pion
71 fField = 2.0; // default assumption: B = 0.4 Tesla
73 fXC = fYC = fR = fC = 0.0;
74 fTanL = fG0 = fDt = fDz = 0.0;
75 fStateR = fStatePhi = fStateZ = fChi2 = fNSteps = 0.0;
79 for (i = 0; i < 6; i++) fPoint[i] = track.fPoint[i];
91 AliITSNeuralTrack::~AliITSNeuralTrack()
94 for (l = 0; l < 6; l++) fPoint[l] = 0;
99 void AliITSNeuralTrack::AssignLabel()
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.
109 Int_t i, j, l, lab, max = 0;
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};
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;
130 for (j = 0; j < max; j++) {
144 j = 0, max = count[0];
145 for (i = 0; i < 18; i++) {
146 if (count[i] > max) {
157 void AliITSNeuralTrack::CleanSlot(Int_t i, Bool_t del)
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.
162 if (i >= 0 && i < 6) {
163 if (del) delete fPoint[i];
170 void AliITSNeuralTrack::GetModuleData(Int_t layer, Int_t &mod, Int_t &pos)
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.
176 if (layer < 0 || layer > 5) {
177 Error("GetModuleData", "Layer out of range: %d", layer);
180 mod = fPoint[layer]->GetModule();
181 pos = fPoint[layer]->GetIndex();
186 void AliITSNeuralTrack::Insert(AliITSNeuralPoint *point)
188 // A trivial method to insert a point in the tracks;
189 // the point is inserted to the slot corresponding to its ITS layer.
193 Int_t layer = point->GetLayer();
194 if (layer < 0 || layer > 6) {
195 Error("Insert", "Layer index %d out of range", layer);
199 fPoint[layer] = point;
204 Int_t AliITSNeuralTrack::OccupationMask()
206 // Returns a byte which maps the occupied slots.
207 // Each bit represents a layer going from the less significant on.
209 Int_t i, check, mask = 0;
210 for (i = 0; i < 6; i++) {
212 if (fPoint[i]) mask |= check;
219 void AliITSNeuralTrack::PrintLabels()
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.
225 cout << "Assigned label = " << fLabel << " -- counted " << fCount << " times: " << endl << endl;
226 for (Int_t i = 0; i < 6; i++) {
227 cout << "Point #" << i + 1 << " --> ";
229 cout << "labels = " << fPoint[i]->GetLabel(0) << ", ";
230 cout << fPoint[i]->GetLabel(1) << ", ";
231 cout << fPoint[i]->GetLabel(2) << endl;
234 cout << "not assigned" << endl;
242 Bool_t AliITSNeuralTrack::AddEL(Int_t layer, Double_t sign)
244 // Calculates the correction for energy loss
246 Double_t width = 0.0;
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;
255 Error("AddEL", "Layer value %d out of range!", layer);
260 if((layer == 5) && (fStatePhi < 0.174 || fStatePhi > 6.100 || (fStatePhi > 2.960 && fStatePhi < 3.31))) {
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
271 printf("Anomaly in AddEL: pt=%8.6f invSqCosL=%8.6f fMass=%8.7f --> beta2 = %8.7f\n", pt, invSqCosL, fMass, beta2);
275 Double_t dE = 0.153 / beta2 * (log(5940. * beta2 / (1. - beta2)) - beta2) * width * 21.82 * invCosL;
276 dE = sign * dE * 0.001;
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.);
289 Bool_t AliITSNeuralTrack::AddMS(Int_t layer)
291 // Calculates the noise perturbation due to multiple scattering
293 Double_t width = 0.0;
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;
302 Error("AddEL", "Layer value %d out of range!", layer);
307 if((layer == 5) && (fStatePhi < 0.174 || fStatePhi > 6.100 || (fStatePhi > 2.960 && fStatePhi < 3.31))) {
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;
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;
320 if(dprova > 0.) q41 = -1. / cosL * TMath::Sqrt(dprova) / (1. + 2. * halfC *fDt);
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));
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;
341 Int_t AliITSNeuralTrack::PropagateTo(Double_t rk)
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.
353 if (rk < TMath::Abs(fDt)) {
354 Error("PropagateTo", Form("Impossible propagation to r (=%17.15g) < Dt (=%17.15g)", rk, fDt));
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);
363 fStatePhi += TMath::ASin(aAk) - TMath::ASin(aAkm1);
364 if(fStatePhi > duepi) fStatePhi -= duepi;
365 if(fStatePhi < 0.) fStatePhi += duepi;
367 Double_t halfC = 0.5 * fC;
368 fStateZ += fTanL / halfC * (TMath::ASin(ak)-TMath::ASin(akm1));
370 Double_t bk = ArgB(rk), bkm1 = ArgB(rkm1);
371 Double_t ck = ArgC(rk), ckm1 = ArgC(rkm1);
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;
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);
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;
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;
414 if (rkm1 < fStateR) // going to greater R --> energy LOSS
416 else // going to smaller R --> energy GAIN
422 Bool_t AliITSNeuralTrack::SeedCovariance()
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.
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++) {
439 // weight and derivatives of phi and zeta w.r.t. various params
440 sumphi += 1./ p->ErrorGetPhi();
443 if (argPhi > 100.0 || argZ > 100.0) {
444 Error("InitCovariance", "Argument error");
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");
453 zL = asin(argZ) / fC;
457 sumz += 1.0 / (p->fError[2] * p->fError[2]);
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;
469 AliITSNeuralPoint *p = 0;
470 Double_t delta, cs, sn, r, argz;
471 Double_t diffC, diffD, diffL, calcC, calcD, calcL;
474 for (l = 0; l < 6; l++) {
477 sn = TMath::Sin(p->GetPhi() - fG0);
478 cs = TMath::Cos(p->GetPhi() - fG0);
480 calcC = (fDt/r - sn) / (2.*fDt*sn - r - fDt*fDt/r);
483 Error("Covariance", "Value too high");
486 calcL = (p->Z() - fDz) * fC / asin(argz);
487 delta = fR*fR + r*r + 2.0*fR*r*sin(p->GetPhi() - fG0);
492 Error("Covariance", Form("Discriminant = %g --- Dt = %g", delta, fDt));
502 diffL = calcL - fTanL;
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;
511 for (l = 0; l < 6; l++) if (fPoint[l]) n++;
512 fMatrix *= 1./(n++ * n);
518 Bool_t AliITSNeuralTrack::Filter(AliITSNeuralPoint *test)
520 // Makes all calculations which apply the Kalman filter to the
521 // stored guess of the state vector, after propagation to a new layer
524 Error("Filter", "Null pointer passed");
529 Double_t rk, phik, zk;
531 phik = test->GetPhi();
536 //////////////////////// Evaluation of the error matrix V /////////
537 Double_t v00 = test->GetError("phi") * rk;
538 Double_t v11 = test->ErrZ();
539 ////////////////////////////////////////////////////////////////////
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);
561 // Calculate R matrix
562 Double_t rold00 = cin00 + v00;
563 Double_t rold10 = cin10;
564 Double_t rold11 = cin11 + v11;
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 ////////////////////////////////////////////////////////////////////
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;
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;
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);
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;
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;
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;
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);
657 Bool_t AliITSNeuralTrack::KalmanFit()
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
665 Int_t l, layer, sign;
667 fStateR = fPoint[0]->GetR2();
668 fStatePhi = fPoint[0]->GetPhi();
669 fStateZ = fPoint[0]->Z();
671 if (!PropagateTo(3.0)) {
672 Error("KalmanFit", "Unsuccessful initialization");
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;
683 Error("KalmanFit", "Not six points!");
686 layer = fPoint[l]->GetLayer();
687 rho = fPoint[l]->GetR2();
688 sign = PropagateTo(rho);
689 if (!sign) return kFALSE;
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));
700 fG0 = fStatePhi - asin(argPhi);
701 fDz = fStateZ - (2.0 * fTanL / fC) * asin(argZ);
705 // Now a Kalman filter i performed going from the actual state position
706 // towards layer 1 position and then propagates to vertex
709 layer = fPoint[l]->GetLayer();
710 rho = fPoint[l]->GetR2();
712 sign = PropagateTo(rho);
713 if (!sign) return kFALSE;
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));
723 fG0 = fStatePhi - asin(argPhi);
724 fDz = fStateZ - (2.0 * fTanL / fC) * asin(argZ);
732 Bool_t AliITSNeuralTrack::RiemannFit()
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;
743 // M1 - matrix of ones
745 for (i = 0; i < 6; i++) m1(i,0) = 1.0;
747 // X - matrix of Rieman projection coordinates
749 for (i = 0; i < 6; i++) {
750 mX(i,0) = fPoint[i]->X();
751 mX(i,1) = fPoint[i]->Y();
752 mX(i,2) = fPoint[i]->GetR2sq();
755 // W - matrix of weights
756 Double_t xterm, yterm, ex, ey;
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();
763 mW(i,i) = fPoint[i]->GetR2sq() / (xterm * xterm * ex * ex + yterm * yterm * ey * ey);
766 // Xm - weighted sample mean
767 Double_t meanX = 0.0, meanY = 0.0, meanW = 0.0, sw = 0.0;
768 for (i = 0; i < 6; i++) {
769 meanX += mW(i,i) * mX(i,0);
770 meanY += mW(i,i) * mX(i,1);
771 meanW += mW(i,i) * mX(i,2);
778 // V - sample covariance matrix
779 for (i = 0; i < 6; i++) {
784 TMatrixD mXt(TMatrixD::kTransposed, mX);
785 TMatrixD mWX(mW, TMatrixD::kMult, mX);
786 TMatrixD mV(mXt, TMatrixD::kMult, mWX);
787 for (i = 0; i < 3; i++) {
788 for (j = i + 1; j < 3; j++) {
789 mV(i,j) = mV(j,i) = (mV(i,j) + mV(j,i)) * 0.5;
793 // Eigenvalue problem solving for V matrix
795 TVectorD eval(3), n(3);
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();
802 TMatrixD evec = mV.EigenVectors(eval);
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);
811 // c - known term in the plane intersection with Riemann axes
812 Double_t c = -(meanX * n(0) + meanY * n(1) + meanW * n(2));
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));
819 Error("RiemannFit", "Radius comed less than zero!!!");
822 fR = TMath::Sqrt(fR);
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];
829 for (i = 1; i < 6; i++) {
830 p = (AliITSNeuralPoint*)fPoint[i];
835 if (temp1 > ref && temp2 < ref)
837 else if (temp1 < ref && temp2 > ref)
839 sumdphi += temp2 - temp1;
842 if (sumdphi < 0.E0) fC = -fC;
843 Double_t diff, angle = TMath::ATan2(fYC, fXC);
845 fG0 = angle + 0.5 * TMath::Pi();
847 fG0 = angle - 0.5 * TMath::Pi();
850 Double_t d = TMath::Sqrt(fXC*fXC + fYC*fYC) - fR;
857 Double_t halfC = 0.5 * fC;
858 Double_t *s = new Double_t[nn], *z = new Double_t[nn], *ws = new Double_t[nn];
859 for (j = 0; j < 6; j++) {
862 s[j] = ArgZ(p->GetR2());
863 if (s[j] > 100.0) return kFALSE;
865 s[j] = asin(s[j]) / halfC;
866 ws[j] = 1.0 / (p->ErrZ()*p->ErrZ());
869 // second tep final fit
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];
882 d = sums2 - sums*sums;
884 fDz = (sums2*sumz - sums*sumsz) / d;
885 fTanL = (sumsz - sums*sumz) / d;
896 void AliITSNeuralTrack::PrintState(Bool_t matrix)
898 // Prints the state vector values.
899 // The argument switches on or off the printing of the covariance matrix.
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";
913 cout << "\nCovariance Matrix: " << endl;
916 cout << "Actual square chi = " << fChi2;
921 Double_t AliITSNeuralTrack::GetDz() const
923 // Double_t argZ = ArgZ(fStateR);
925 // Error("GetDz", "Too large value: %g", argZ);
928 // fDz = fStateZ - (2.0 * fTanL / fC) * asin(argZ);
934 Double_t AliITSNeuralTrack::GetGamma() const
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);
942 // fG0 = fStatePhi - asin(argPhi);
948 Double_t AliITSNeuralTrack::GetPhi(Double_t r) const
950 // Gives the value of azymuthal coordinate in the helix
951 // as a function of cylindric radius
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(); }
963 Double_t AliITSNeuralTrack::GetZ(Double_t r) const
965 // gives the value of Z in the helix
966 // as a function of cylindric radius
968 Double_t arg = ArgZ(r);
969 if (arg > 0.9) return 0.0;
970 return fDz + fTanL * asin(arg) / fC;
975 Double_t AliITSNeuralTrack::GetdEdX()
977 // total energy loss of the track
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);
992 for (i = 0; i < 3; i++) {
993 if (q[i] <= q[i + 1]) continue;
1006 dedx = (q[0] + q[1]) / 2.;
1012 void AliITSNeuralTrack::SetVertex(Double_t *pos, Double_t *err)
1014 // Stores vertex data
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);
1031 AliITSIOTrack* AliITSNeuralTrack::ExportIOtrack(Int_t min)
1033 // Exports an object in the standard format for reconstructed tracks
1036 AliITSIOTrack *track = new AliITSIOTrack;
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));
1046 track->SetLabel(IsGood(min) ? fLabel : -fLabel);
1047 track->SetTPCLabel(-1);
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());
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);
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);
1076 track->SetPid(fPDG);
1077 track->SetMass(fMass);
1083 //====================================================================================
1084 //============================ PRIVATE METHODS ============================
1085 //====================================================================================
1088 Double_t AliITSNeuralTrack::ArgPhi(Double_t r) const
1090 // calculates the expression ((1/2)Cr + (1 + (1/2)CD) D/r) / (1 + CD)
1092 Double_t arg, num, den;
1093 num = (0.5 * fC * r) + (1. + (0.5 * fC * fDt)) * (fDt / r);
1094 den = 1. + fC * fDt;
1096 Error("ArgPhi", "Denominator = 0!");
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);
1108 Double_t AliITSNeuralTrack::ArgZ(Double_t r) const
1110 // calculates the expression (1/2)C * sqrt( (r^2 - Dt^2) / (1 + CD) )
1113 arg = (r * r - fDt * fDt) / (1. + fC * fDt);
1115 if (TMath::Abs(arg) < 1.E-6) arg = 0.;
1117 Error("ArgZ", "Square root argument error: %17.15g < 0", arg);
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);
1130 Double_t AliITSNeuralTrack::ArgB(Double_t r) const
1135 arg = (r*r - fDt*fDt);
1136 arg /= (r*(1.+ fC*fDt)*(1.+ fC*fDt));
1142 Double_t AliITSNeuralTrack::ArgC(Double_t r) const
1147 arg = (1./r - fC * ArgPhi(r));