ITS/AliITSNeuralTrack.cxx

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