ITS/AliITSNeuralTrack.cxx

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