1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
20 // The format of output data from Neural Tracker
21 // It can export data in the format of AliITSIOTrack
23 // Compatibility adaptation to V2 tracking is on the way.
24 // Author: A. Pulvirenti
26 #include <Riostream.h>
30 //#include <TObject.h>
34 //#include <TObjArray.h>
37 #if ROOT_VERSION_CODE >= 262146
38 #include <TMatrixDEigen.h>
41 //#include "AliITSVertex.h"
42 #include "AliITSIOTrack.h"
43 #include "AliITSNeuralPoint.h"
45 #include "AliITSNeuralTrack.h"
49 ClassImp(AliITSNeuralTrack)
53 AliITSNeuralTrack::AliITSNeuralTrack() :
68 fMass(0.1396),// default assumption: pion
69 fField(2.0),// default assumption: B = 0.4 Tesla
74 // Default constructor
78 for (i = 0; i < 6; i++) fPoint[i] = 0;
90 AliITSNeuralTrack::AliITSNeuralTrack(const AliITSNeuralTrack &track)
91 : TObject((TObject&)track),
106 fMass(0.1396),// default assumption: pion
107 fField(2.0),// default assumption: B = 0.4 Tesla
116 fMass = 0.1396; // default assumption: pion
117 fField = 2.0; // default assumption: B = 0.4 Tesla
119 fXC = fYC = fR = fC = 0.0;
120 fTanL = fG0 = fDt = fDz = 0.0;
121 fStateR = fStatePhi = fStateZ = fChi2 = fNSteps = 0.0;
125 for (i = 0; i < 6; i++) fPoint[i] = track.fPoint[i];
130 fVertex.ErrX() = 0.0;
131 fVertex.ErrY() = 0.0;
132 fVertex.ErrZ() = 0.0;
138 AliITSNeuralTrack& AliITSNeuralTrack::operator=(const AliITSNeuralTrack& track){
139 //assignment operator
140 this->~AliITSNeuralTrack();
141 new(this) AliITSNeuralTrack(track);
144 AliITSNeuralTrack::~AliITSNeuralTrack()
147 for (l = 0; l < 6; l++) fPoint[l] = 0;
152 void AliITSNeuralTrack::AssignLabel()
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.
162 Int_t i, j, l, lab, max = 0;
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};
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;
183 for (j = 0; j < max; j++) {
197 j = 0, max = count[0];
198 for (i = 0; i < 18; i++) {
199 if (count[i] > max) {
210 void AliITSNeuralTrack::CleanSlot(Int_t i, Bool_t del)
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.
215 if (i >= 0 && i < 6) {
216 if (del) delete fPoint[i];
223 void AliITSNeuralTrack::GetModuleData(Int_t layer, Int_t &mod, Int_t &pos)
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.
229 if (layer < 0 || layer > 5) {
230 Error("GetModuleData", "Layer out of range: %d", layer);
233 mod = fPoint[layer]->GetModule();
234 pos = fPoint[layer]->GetIndex();
239 void AliITSNeuralTrack::Insert(AliITSNeuralPoint *point)
241 // A trivial method to insert a point in the tracks;
242 // the point is inserted to the slot corresponding to its ITS layer.
246 Int_t layer = point->GetLayer();
247 if (layer < 0 || layer > 6) {
248 Error("Insert", "Layer index %d out of range", layer);
252 fPoint[layer] = point;
257 Int_t AliITSNeuralTrack::OccupationMask() const
259 // Returns a byte which maps the occupied slots.
260 // Each bit represents a layer going from the less significant on.
262 Int_t i, check, mask = 0;
263 for (i = 0; i < 6; i++) {
265 if (fPoint[i]) mask |= check;
272 void AliITSNeuralTrack::PrintLabels()
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.
278 cout << "Assigned label = " << fLabel << " -- counted " << fCount << " times: " << endl << endl;
279 for (Int_t i = 0; i < 6; i++) {
280 cout << "Point #" << i + 1 << " --> ";
282 cout << "labels = " << fPoint[i]->GetLabel(0) << ", ";
283 cout << fPoint[i]->GetLabel(1) << ", ";
284 cout << fPoint[i]->GetLabel(2) << endl;
287 cout << "not assigned" << endl;
295 Bool_t AliITSNeuralTrack::AddEL(Int_t layer, Double_t sign)
297 // Calculates the correction for energy loss
299 Double_t width = 0.0;
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;
308 Error("AddEL", "Layer value %d out of range!", layer);
313 if((layer == 5) && (fStatePhi < 0.174 || fStatePhi > 6.100 || (fStatePhi > 2.960 && fStatePhi < 3.31))) {
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
324 printf("Anomaly in AddEL: pt=%8.6f invSqCosL=%8.6f fMass=%8.7f --> beta2 = %8.7f\n", pt, invSqCosL, fMass, beta2);
328 Double_t dE = 0.153 / beta2 * (log(5940. * beta2 / (1. - beta2)) - beta2) * width * 21.82 * invCosL;
329 dE = sign * dE * 0.001;
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.);
342 Bool_t AliITSNeuralTrack::AddMS(Int_t layer)
344 // Calculates the noise perturbation due to multiple scattering
346 Double_t width = 0.0;
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;
355 Error("AddEL", "Layer value %d out of range!", layer);
360 if((layer == 5) && (fStatePhi < 0.174 || fStatePhi > 6.100 || (fStatePhi > 2.960 && fStatePhi < 3.31))) {
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;
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;
373 if(dprova > 0.) q41 = -1. / cosL * TMath::Sqrt(dprova) / (1. + 2. * halfC *fDt);
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));
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;
394 Int_t AliITSNeuralTrack::PropagateTo(Double_t rk)
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.
406 if (rk < TMath::Abs(fDt)) {
407 Error("PropagateTo", Form("Impossible propagation to r (=%17.15g) < Dt (=%17.15g)", rk, fDt));
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);
416 fStatePhi += TMath::ASin(aAk) - TMath::ASin(aAkm1);
417 if(fStatePhi > duepi) fStatePhi -= duepi;
418 if(fStatePhi < 0.) fStatePhi += duepi;
420 Double_t halfC = 0.5 * fC;
421 fStateZ += fTanL / halfC * (TMath::ASin(ak)-TMath::ASin(akm1));
423 Double_t bk = ArgB(rk), bkm1 = ArgB(rkm1);
424 Double_t ck = ArgC(rk), ckm1 = ArgC(rkm1);
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;
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);
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;
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;
467 if (rkm1 < fStateR) // going to greater R --> energy LOSS
469 else // going to smaller R --> energy GAIN
475 Bool_t AliITSNeuralTrack::SeedCovariance()
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.
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++) {
492 // weight and derivatives of phi and zeta w.r.t. various params
493 sumphi += 1./ p->ErrorGetPhi();
496 if (argPhi > 100.0 || argZ > 100.0) {
497 Error("InitCovariance", "Argument error");
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");
506 zL = asin(argZ) / fC;
510 sumz += 1.0 / (p->fError[2] * p->fError[2]);
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;
522 AliITSNeuralPoint *p = 0;
523 Double_t delta, cs, sn, r, argz;
524 Double_t diffC, diffD, diffL, calcC, calcD, calcL;
527 for (l = 0; l < 6; l++) {
530 sn = TMath::Sin(p->GetPhi() - fG0);
531 cs = TMath::Cos(p->GetPhi() - fG0);
533 calcC = (fDt/r - sn) / (2.*fDt*sn - r - fDt*fDt/r);
536 Error("Covariance", "Value too high");
539 calcL = (p->Z() - fDz) * fC / asin(argz);
540 delta = fR*fR + r*r + 2.0*fR*r*sin(p->GetPhi() - fG0);
545 Error("Covariance", Form("Discriminant = %g --- Dt = %g", delta, fDt));
555 diffL = calcL - fTanL;
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;
564 for (l = 0; l < 6; l++) if (fPoint[l]) n++;
565 fMatrix *= 1./(n * (n+1));
571 Bool_t AliITSNeuralTrack::Filter(AliITSNeuralPoint *test)
573 // Makes all calculations which apply the Kalman filter to the
574 // stored guess of the state vector, after propagation to a new layer
577 Error("Filter", "Null pointer passed");
582 Double_t rk, phik, zk;
584 phik = test->GetPhi();
589 //////////////////////// Evaluation of the error matrix V /////////
590 Double_t v00 = test->GetError("phi") * rk;
591 Double_t v11 = test->ErrZ();
592 ////////////////////////////////////////////////////////////////////
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);
614 // Calculate R matrix
615 Double_t rold00 = cin00 + v00;
616 Double_t rold10 = cin10;
617 Double_t rold11 = cin11 + v11;
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 ////////////////////////////////////////////////////////////////////
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;
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;
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);
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;
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;
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;
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);
710 Bool_t AliITSNeuralTrack::KalmanFit()
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
718 Int_t l, layer, sign;
720 fStateR = fPoint[0]->GetR2();
721 fStatePhi = fPoint[0]->GetPhi();
722 fStateZ = fPoint[0]->Z();
724 if (!PropagateTo(3.0)) {
725 Error("KalmanFit", "Unsuccessful initialization");
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;
736 Error("KalmanFit", "Not six points!");
739 layer = fPoint[l]->GetLayer();
740 rho = fPoint[l]->GetR2();
741 sign = PropagateTo(rho);
742 if (!sign) return kFALSE;
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));
753 fG0 = fStatePhi - asin(argPhi);
754 fDz = fStateZ - (2.0 * fTanL / fC) * asin(argZ);
758 // Now a Kalman filter i performed going from the actual state position
759 // towards layer 1 position and then propagates to vertex
762 layer = fPoint[l]->GetLayer();
763 rho = fPoint[l]->GetR2();
765 sign = PropagateTo(rho);
766 if (!sign) return kFALSE;
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));
776 fG0 = fStatePhi - asin(argPhi);
777 fDz = fStateZ - (2.0 * fTanL / fC) * asin(argZ);
785 Bool_t AliITSNeuralTrack::RiemannFit()
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;
796 // M1 - matrix of ones
798 for (i = 0; i < 6; i++) m1(i,0) = 1.0;
800 // X - matrix of Rieman projection coordinates
802 for (i = 0; i < 6; i++) {
803 mX(i,0) = fPoint[i]->X();
804 mX(i,1) = fPoint[i]->Y();
805 mX(i,2) = fPoint[i]->GetR2sq();
808 // W - matrix of weights
809 Double_t xterm, yterm, ex, ey;
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();
816 mW(i,i) = fPoint[i]->GetR2sq() / (xterm * xterm * ex * ex + yterm * yterm * ey * ey);
819 // Xm - weighted sample mean
820 Double_t meanX = 0.0, meanY = 0.0, meanW = 0.0, sw = 0.0;
821 for (i = 0; i < 6; i++) {
822 meanX += mW(i,i) * mX(i,0);
823 meanY += mW(i,i) * mX(i,1);
824 meanW += mW(i,i) * mX(i,2);
831 // V - sample covariance matrix
832 for (i = 0; i < 6; i++) {
837 TMatrixD mXt(TMatrixD::kTransposed, mX);
838 TMatrixD mWX(mW, TMatrixD::kMult, mX);
839 TMatrixD mV(mXt, TMatrixD::kMult, mWX);
840 for (i = 0; i < 3; i++) {
841 for (j = i + 1; j < 3; j++) {
842 mV(i,j) = mV(j,i) = (mV(i,j) + mV(j,i)) * 0.5;
846 // Eigenvalue problem solving for V matrix
848 TVectorD eval(3), n(3);
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();
855 TMatrixD evec = mV.EigenVectors(eval);
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);
864 // c - known term in the plane intersection with Riemann axes
865 Double_t c = -(meanX * n(0) + meanY * n(1) + meanW * n(2));
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));
872 Error("RiemannFit", "Radius comed less than zero!!!");
875 fR = TMath::Sqrt(fR);
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];
882 for (i = 1; i < 6; i++) {
883 p = (AliITSNeuralPoint*)fPoint[i];
888 if (temp1 > ref && temp2 < ref)
890 else if (temp1 < ref && temp2 > ref)
892 sumdphi += temp2 - temp1;
895 if (sumdphi < 0.E0) fC = -fC;
896 Double_t diff, angle = TMath::ATan2(fYC, fXC);
898 fG0 = angle + 0.5 * TMath::Pi();
900 fG0 = angle - 0.5 * TMath::Pi();
903 Double_t d = TMath::Sqrt(fXC*fXC + fYC*fYC) - fR;
910 Double_t halfC = 0.5 * fC;
911 Double_t *s = new Double_t[nn], *z = new Double_t[nn], *ws = new Double_t[nn];
912 for (j = 0; j < 6; j++) {
915 s[j] = ArgZ(p->GetR2());
916 if (s[j] > 100.0) return kFALSE;
918 s[j] = asin(s[j]) / halfC;
919 ws[j] = 1.0 / (p->ErrZ()*p->ErrZ());
922 // second tep final fit
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];
935 d = sums2 - sums*sums;
937 fDz = (sums2*sumz - sums*sumsz) / d;
938 fTanL = (sumsz - sums*sumz) / d;
949 void AliITSNeuralTrack::PrintState(Bool_t matrix)
951 // Prints the state vector values.
952 // The argument switches on or off the printing of the covariance matrix.
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";
966 cout << "\nCovariance Matrix: " << endl;
969 cout << "Actual square chi = " << fChi2;
974 Double_t AliITSNeuralTrack::GetDz() const
976 // Double_t argZ = ArgZ(fStateR);
978 // Error("GetDz", "Too large value: %g", argZ);
981 // fDz = fStateZ - (2.0 * fTanL / fC) * asin(argZ);
987 Double_t AliITSNeuralTrack::GetGamma() const
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);
995 // fG0 = fStatePhi - asin(argPhi);
1001 Double_t AliITSNeuralTrack::GetPhi(Double_t r) const
1003 // Gives the value of azymuthal coordinate in the helix
1004 // as a function of cylindric radius
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(); }
1016 Double_t AliITSNeuralTrack::GetZ(Double_t r) const
1018 // gives the value of Z in the helix
1019 // as a function of cylindric radius
1021 Double_t arg = ArgZ(r);
1022 if (arg > 0.9) return 0.0;
1023 return fDz + fTanL * asin(arg) / fC;
1028 Double_t AliITSNeuralTrack::GetdEdX()
1030 // total energy loss of the track
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);
1045 for (i = 0; i < 3; i++) {
1046 if (q[i] <= q[i + 1]) continue;
1047 Double_t tmp = q[i];
1059 dedx = (q[0] + q[1]) / 2.;
1065 void AliITSNeuralTrack::SetVertex(Double_t *pos, Double_t *err)
1067 // Stores vertex data
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);
1084 AliITSIOTrack* AliITSNeuralTrack::ExportIOtrack(Int_t min)
1086 // Exports an object in the standard format for reconstructed tracks
1089 AliITSIOTrack *track = new AliITSIOTrack;
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));
1099 track->SetLabel(IsGood(min) ? fLabel : -fLabel);
1100 track->SetTPCLabel(-1);
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());
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);
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);
1129 track->SetPid(fPDG);
1130 track->SetMass(fMass);
1136 //====================================================================================
1137 //============================ PRIVATE METHODS ============================
1138 //====================================================================================
1141 Double_t AliITSNeuralTrack::ArgPhi(Double_t r) const
1143 // calculates the expression ((1/2)Cr + (1 + (1/2)CD) D/r) / (1 + CD)
1145 Double_t arg, num, den;
1146 num = (0.5 * fC * r) + (1. + (0.5 * fC * fDt)) * (fDt / r);
1147 den = 1. + fC * fDt;
1149 Error("ArgPhi", "Denominator = 0!");
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);
1161 Double_t AliITSNeuralTrack::ArgZ(Double_t r) const
1163 // calculates the expression (1/2)C * sqrt( (r^2 - Dt^2) / (1 + CD) )
1166 arg = (r * r - fDt * fDt) / (1. + fC * fDt);
1168 if (TMath::Abs(arg) < 1.E-6) arg = 0.;
1170 Error("ArgZ", "Square root argument error: %17.15g < 0", arg);
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);
1183 Double_t AliITSNeuralTrack::ArgB(Double_t r) const
1188 arg = (r*r - fDt*fDt);
1189 arg /= (r*(1.+ fC*fDt)*(1.+ fC*fDt));
1195 Double_t AliITSNeuralTrack::ArgC(Double_t r) const
1200 arg = (1./r - fC * ArgPhi(r));