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() :
  54 fXC(0.0),
  55 fYC(0.0),
  56 fR(0.0),
  57 fC(0.0),
  58 fTanL(0.0),
  59 fG0(0.0),
  60 fDt(0.0),
  61 fDz(0.0),
  62 fStateR(0.0),
  63 fStatePhi(0.0),
  64 fStateZ(0.0),
  65 fMatrix(5,5),
  66 fChi2(0.0),
  67 fNSteps(0.0),
  68 fMass(0.1396),// default assumption: pion
  69 fField(2.0),// default assumption: B = 0.4 Tesla
  70 fPDG(0),
  71 fLabel(0),
  72 fCount(0),
  73 fVertex(){
  74         // Default constructor
  75
  76         Int_t i;
  77
  78         for (i = 0; i < 6; i++) fPoint[i] = 0;
  79
  80         fVertex.X() = 0.0;
  81         fVertex.Y() = 0.0;
  82         fVertex.Z() = 0.0;
  83         fVertex.ErrX() = 0.0;
  84         fVertex.ErrY() = 0.0;
  85         fVertex.ErrZ() = 0.0;
  86 }
  87 //
  88 //
  89 //
  90 AliITSNeuralTrack::AliITSNeuralTrack(const AliITSNeuralTrack &track)
  91 : TObject((TObject&)track),
  92 fXC(0.0),
  93 fYC(0.0),
  94 fR(0.0),
  95 fC(0.0),
  96 fTanL(0.0),
  97 fG0(0.0),
  98 fDt(0.0),
  99 fDz(0.0),
 100 fStateR(0.0),
 101 fStatePhi(0.0),
 102 fStateZ(0.0),
 103 fMatrix(5,5),
 104 fChi2(0.0),
 105 fNSteps(0.0),
 106 fMass(0.1396),// default assumption: pion
 107 fField(2.0),// default assumption: B = 0.4 Tesla
 108 fPDG(0),
 109 fLabel(0),
 110 fCount(0),
 111 fVertex(){
 112 // copy constructor
 113
 114         Int_t i;
 115
 116         fMass  = 0.1396;   // default assumption: pion
 117         fField = 2.0;      // default assumption: B = 0.4 Tesla
 118
 119         fXC = fYC = fR = fC = 0.0;
 120         fTanL = fG0 = fDt = fDz = 0.0;
 121         fStateR = fStatePhi = fStateZ = fChi2 = fNSteps = 0.0;
 122
 123         fLabel = 0;
 124         fCount = 0;
 125         for (i = 0; i < 6; i++) fPoint[i] = track.fPoint[i];
 126
 127         fVertex.X() = 0.0;
 128         fVertex.Y() = 0.0;
 129         fVertex.Z() = 0.0;
 130         fVertex.ErrX() = 0.0;
 131         fVertex.ErrY() = 0.0;
 132         fVertex.ErrZ() = 0.0;
 133 }
 134 //
 135 //
 136 //
 137
 138 AliITSNeuralTrack& AliITSNeuralTrack::operator=(const AliITSNeuralTrack& track){
 139   //assignment operator
 140   this->~AliITSNeuralTrack();
 141   new(this) AliITSNeuralTrack(track);
 142   return *this;
 143 }
 144 AliITSNeuralTrack::~AliITSNeuralTrack()
 145 {
 146         Int_t l;
 147         for (l = 0; l < 6; l++) fPoint[l] = 0;
 148 }
 149 //
 150 //
 151 //
 152 void AliITSNeuralTrack::AssignLabel()
 153 {
 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.
 160
 161         Bool_t found;
 162         Int_t i, j, l, lab, max = 0;
 163
 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};
 167
 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;
 182                         found = kFALSE;
 183                         for (j = 0; j < max; j++) {
 184                                 if (idx[j] == lab) {
 185                                         count[j]++;
 186                                         found = kTRUE;
 187                                 }
 188                         }
 189                         if(!found) {
 190                                 max++;
 191                                 idx[max - 1] = lab;
 192                                 count[max - 1] = 1;
 193                         }
 194                 }
 195         }
 196
 197         j = 0, max = count[0];
 198         for (i = 0; i < 18; i++) {
 199                 if (count[i] > max) {
 200                         j = i;
 201                         max = count[i];
 202                 }
 203         }
 204         fLabel = idx[j];
 205         fCount = count[j];
 206 }
 207 //
 208 //
 209 //
 210 void AliITSNeuralTrack::CleanSlot(Int_t i, Bool_t del)
 211 {
 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.
 214
 215         if (i >= 0 && i < 6) {
 216                 if (del) delete fPoint[i];
 217                 fPoint[i] = 0;
 218         }
 219 }
 220 //
 221 //
 222 //
 223 void AliITSNeuralTrack::GetModuleData(Int_t layer, Int_t &mod, Int_t &pos)
 224 {
 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.
 228
 229         if (layer < 0 || layer > 5) {
 230                 Error("GetModuleData", "Layer out of range: %d", layer);
 231                 return;
 232         }
 233         mod = fPoint[layer]->GetModule();
 234         pos = fPoint[layer]->GetIndex();
 235 }
 236 //
 237 //
 238 //
 239 void AliITSNeuralTrack::Insert(AliITSNeuralPoint *point)
 240 {
 241 // A trivial method to insert a point in the tracks;
 242 // the point is inserted to the slot corresponding to its ITS layer.
 243
 244         if (!point) return;
 245
 246         Int_t layer = point->GetLayer();
 247         if (layer < 0 || layer > 6) {
 248                 Error("Insert", "Layer index %d out of range", layer);
 249                 return;
 250         }
 251
 252         fPoint[layer] = point;
 253 }
 254 //
 255 //
 256 //
 257 Int_t AliITSNeuralTrack::OccupationMask() const
 258 {
 259 // Returns a byte which maps the occupied slots.
 260 // Each bit represents a layer going from the less significant on.
 261
 262         Int_t i, check, mask = 0;
 263         for (i = 0; i < 6; i++) {
 264                 check = 1 << i;
 265                 if (fPoint[i]) mask |= check;
 266         }
 267         return mask;
 268 }
 269 //
 270 //
 271 //
 272 void AliITSNeuralTrack::PrintLabels()
 273 {
 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.
 277
 278         cout << "Assigned label = " << fLabel << " -- counted " << fCount << " times: " << endl << endl;
 279         for (Int_t i = 0; i < 6; i++) {
 280                 cout << "Point #" << i + 1 << " --> ";
 281                 if (fPoint[i]) {
 282                         cout << "labels = " << fPoint[i]->GetLabel(0) << ", ";
 283                         cout << fPoint[i]->GetLabel(1) << ", ";
 284                         cout << fPoint[i]->GetLabel(2) << endl;
 285                 }
 286                 else {
 287                         cout << "not assigned" << endl;
 288                 }
 289         }
 290         cout << endl;
 291 }
 292 //
 293 //
 294 //
 295 Bool_t AliITSNeuralTrack::AddEL(Int_t layer, Double_t sign)
 296 {
 297 // Calculates the correction for energy loss
 298
 299         Double_t width = 0.0;
 300         switch (layer) {
 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;
 307                 default:
 308                         Error("AddEL", "Layer value %d out of range!", layer);
 309                         return kFALSE;
 310         }
 311         width *= 1.7;
 312
 313         if((layer == 5) && (fStatePhi < 0.174 || fStatePhi > 6.100 || (fStatePhi > 2.960 && fStatePhi < 3.31))) {
 314                 width += 0.012;
 315         }
 316
 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
 323         if (beta2 == 0.0) {
 324                 printf("Anomaly in AddEL: pt=%8.6f invSqCosL=%8.6f fMass=%8.7f --> beta2 = %8.7f\n", pt, invSqCosL, fMass, beta2);
 325                 return kFALSE;
 326         }
 327
 328         Double_t dE = 0.153 / beta2 * (log(5940. * beta2 / (1. - beta2)) - beta2) * width * 21.82 * invCosL;
 329         dE = sign * dE * 0.001;
 330
 331         energy += dE;
 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.);
 336
 337         return kTRUE;
 338 }
 339 //
 340 //
 341 //
 342 Bool_t AliITSNeuralTrack::AddMS(Int_t layer)
 343 {
 344 // Calculates the noise perturbation due to multiple scattering
 345
 346         Double_t width = 0.0;
 347         switch (layer) {
 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;
 354                 default:
 355                         Error("AddEL", "Layer value %d out of range!", layer);
 356                         return kFALSE;
 357         }
 358         width *= 1.7;
 359
 360         if((layer == 5) && (fStatePhi < 0.174 || fStatePhi > 6.100 || (fStatePhi > 2.960 && fStatePhi < 3.31))) {
 361                 width += 0.012;
 362         }
 363
 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;
 368
 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;
 372         Double_t q41 = 0.;
 373         if(dprova > 0.) q41 = -1. / cosL * TMath::Sqrt(dprova) / (1. + 2. * halfC *fDt);
 374
 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));
 378
 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;
 388
 389         return kTRUE;
 390 }
 391 //
 392 //
 393 //
 394 Int_t AliITSNeuralTrack::PropagateTo(Double_t rk)
 395 {
 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.
 405
 406         if (rk < TMath::Abs(fDt)) {
 407                 Error("PropagateTo", Form("Impossible propagation to r (=%17.15g) < Dt (=%17.15g)", rk, fDt));
 408                 return 0;
 409         }
 410
 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);
 415
 416         fStatePhi += TMath::ASin(aAk) - TMath::ASin(aAkm1);
 417         if(fStatePhi > duepi) fStatePhi -= duepi;
 418         if(fStatePhi < 0.) fStatePhi += duepi;
 419
 420         Double_t halfC = 0.5 * fC;
 421         fStateZ += fTanL / halfC * (TMath::ASin(ak)-TMath::ASin(akm1));
 422
 423         Double_t bk = ArgB(rk), bkm1 = ArgB(rkm1);
 424         Double_t ck = ArgC(rk), ckm1 = ArgC(rkm1);
 425
 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;
 430
 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);
 446
 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;
 454
 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;
 464
 465         fStateR = rk;
 466
 467         if (rkm1 < fStateR)  // going to greater R --> energy LOSS
 468                 return -1;
 469         else                 // going to smaller R --> energy GAIN
 470                 return 1;
 471 }
 472 //
 473 //
 474 //
 475 Bool_t AliITSNeuralTrack::SeedCovariance()
 476 {
 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.
 482
 483         /*
 484         Int_t i, j;
 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++) {
 489                 p = At(i);
 490                 if (!p) continue;
 491                 r = p->GetR2();
 492                 // weight and derivatives of phi and zeta w.r.t. various params
 493                 sumphi += 1./ p->ErrorGetPhi();
 494                 argPhi = ArgPhi(r);
 495                 argZ = ArgZ(r);
 496                 if (argPhi > 100.0 || argZ > 100.0) {
 497                         Error("InitCovariance", "Argument error");
 498                         return kFALSE;
 499                 }
 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");
 504                         return kFALSE;
 505                 }
 506                 zL = asin(argZ) / fC;
 507                 sumL += zL * zL;
 508                 sumC += phiC * phiC;
 509                 sumD += phiD * phiD;
 510                 sumz += 1.0 / (p->fError[2] * p->fError[2]);
 511         }
 512
 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;
 519         fMatrix.Print();
 520         */
 521
 522         AliITSNeuralPoint *p = 0;
 523         Double_t delta, cs, sn, r, argz;
 524         Double_t diffC, diffD, diffL, calcC, calcD, calcL;
 525
 526         Int_t l;
 527         for (l = 0; l < 6; l++) {
 528                 p = fPoint[l];
 529                 if (!p) break;
 530                 sn = TMath::Sin(p->GetPhi() - fG0);
 531                 cs = TMath::Cos(p->GetPhi() - fG0);
 532                 r  = p->GetR2();
 533                 calcC = (fDt/r - sn) / (2.*fDt*sn - r - fDt*fDt/r);
 534                 argz = ArgZ(r);
 535                 if (argz > 1000.0) {
 536                         Error("Covariance", "Value too high");
 537                         return kFALSE;
 538                 }
 539                 calcL = (p->Z() - fDz) * fC / asin(argz);
 540                 delta = fR*fR + r*r + 2.0*fR*r*sin(p->GetPhi() - fG0);
 541                 if (delta < 0.E0) {
 542                         if (delta >= -0.5)
 543                                 delta = 0.;
 544                         else {
 545                                 Error("Covariance", Form("Discriminant = %g --- Dt = %g", delta, fDt));
 546                                 return kFALSE;
 547                         }
 548                 }
 549                 delta = sqrt(delta);
 550                 if (fC >= 0)
 551                         calcD = delta - fR;
 552                 else
 553                         calcD = fR - delta;
 554                 diffD = calcD - fDt;
 555                 diffL = calcL - fTanL;
 556                 diffC = fC - calcC;
 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;
 562         }
 563         Double_t n = 0.;
 564         for (l = 0; l < 6; l++) if (fPoint[l]) n++;
 565         fMatrix *= 1./(n * (n+1));
 566         return kTRUE;
 567 }
 568 //
 569 //
 570 //
 571 Bool_t AliITSNeuralTrack::Filter(AliITSNeuralPoint *test)
 572 {
 573         // Makes all calculations which apply the Kalman filter to the
 574         // stored guess of the state vector, after propagation to a new layer
 575
 576         if (!test) {
 577                 Error("Filter", "Null pointer passed");
 578                 return kFALSE;
 579         }
 580
 581         Double_t m[2];
 582         Double_t rk, phik, zk;
 583         rk = test->GetR2();
 584         phik = test->GetPhi();
 585         zk = test->Z();
 586         m[0]=phik;
 587         m[1]=zk;
 588
 589         //////////////////////// Evaluation of the error matrix V  /////////
 590         Double_t v00 = test->GetError("phi") * rk;
 591         Double_t v11 = test->ErrZ();
 592         ////////////////////////////////////////////////////////////////////
 593
 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);
 613
 614         // Calculate R matrix
 615         Double_t rold00 = cin00 + v00;
 616         Double_t rold10 = cin10;
 617         Double_t rold11 = cin11 + v11;
 618
 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         ////////////////////////////////////////////////////////////////////
 625
 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;
 637
 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;
 642         x2 = fDt;
 643         x3 = fTanL;
 644         x4 = fC;
 645
 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);
 652
 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;
 657
 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;
 673
 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;
 690
 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);
 702         fChi2 += chi2;
 703         fNSteps++;
 704
 705         return kTRUE;
 706 }
 707 //
 708 //
 709 //
 710 Bool_t AliITSNeuralTrack::KalmanFit()
 711 {
 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
 716
 717         Double_t rho;
 718         Int_t l, layer, sign;
 719
 720         fStateR = fPoint[0]->GetR2();
 721         fStatePhi = fPoint[0]->GetPhi();
 722         fStateZ = fPoint[0]->Z();
 723
 724         if (!PropagateTo(3.0)) {
 725                 Error("KalmanFit", "Unsuccessful initialization");
 726                 return kFALSE;
 727         }
 728         l=0;
 729
 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;
 734         while (l <= 5) {
 735                 if (!fPoint[l]) {
 736                         Error("KalmanFit", "Not six points!");
 737                         return kFALSE;
 738                 }
 739                 layer = fPoint[l]->GetLayer();
 740                 rho = fPoint[l]->GetR2();
 741                 sign = PropagateTo(rho);
 742                 if (!sign) return kFALSE;
 743                 AddEL(layer, -1.0);
 744                 AddMS(layer);
 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));
 751                         return kFALSE;
 752                 }
 753                 fG0 = fStatePhi - asin(argPhi);
 754                 fDz = fStateZ - (2.0 * fTanL / fC) * asin(argZ);
 755                 l++;
 756         }
 757
 758         // Now a Kalman filter i performed going from the actual state position
 759         // towards layer 1 position and then propagates to vertex
 760         if (l >= 5) l = 5;
 761         while (l >= 1) {
 762                 layer = fPoint[l]->GetLayer();
 763                 rho = fPoint[l]->GetR2();
 764                 AddEL(layer, 1.0);
 765                 sign = PropagateTo(rho);
 766                 if (!sign) return kFALSE;
 767                 AddMS(layer);
 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));
 774                         return kFALSE;
 775                 }
 776                 fG0 = fStatePhi - asin(argPhi);
 777                 fDz = fStateZ - (2.0 * fTanL / fC) * asin(argZ);
 778                 l--;
 779         }
 780         return kTRUE;
 781 }
 782 //
 783 //
 784 //
 785 Bool_t AliITSNeuralTrack::RiemannFit()
 786 {
 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;
 793
 794         Int_t i, j;
 795
 796         // M1 - matrix of ones
 797         TMatrixD m1(6,1);
 798         for (i = 0; i < 6; i++) m1(i,0) = 1.0;
 799
 800         // X - matrix of Rieman projection coordinates
 801         TMatrixD mX(6,3);
 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();
 806         }
 807
 808         // W - matrix of weights
 809         Double_t xterm, yterm, ex, ey;
 810         TMatrixD mW(6,6);
 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);
 817         }
 818
 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);
 825                 sw += mW(i,i);
 826         }
 827         meanX /= sw;
 828         meanY /= sw;
 829         meanW /= sw;
 830
 831         // V - sample covariance matrix
 832         for (i = 0; i < 6; i++) {
 833                 mX(i,0) -= meanX;
 834                 mX(i,1) -= meanY;
 835                 mX(i,2) -= meanW;
 836         }
 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;
 843                 }
 844         }
 845
 846         // Eigenvalue problem solving for V matrix
 847         Int_t ileast = 0;
 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();
 854 #else
 855         TMatrixD evec = mV.EigenVectors(eval);
 856 #endif
 857
 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);
 863
 864         // c - known term in the plane intersection with Riemann axes
 865         Double_t c = -(meanX * n(0) + meanY * n(1) + meanW * n(2));
 866
 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));
 870
 871         if (fR <= 0.E0) {
 872                 Error("RiemannFit", "Radius comed less than zero!!!");
 873                 return kFALSE;
 874         }
 875         fR = TMath::Sqrt(fR);
 876         fC = 1.0 / fR;
 877
 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];
 881         phi1 = p->GetPhi();
 882         for (i = 1; i < 6; i++) {
 883                 p = (AliITSNeuralPoint*)fPoint[i];
 884                 if (!p) break;
 885                 phi2 = p->GetPhi();
 886                 temp1 = phi1;
 887                 temp2 = phi2;
 888                 if (temp1 > ref && temp2 < ref)
 889                         temp2 += 2.0 * ref;
 890                 else if (temp1 < ref && temp2 > ref)
 891                         temp1 += 2.0 * ref;
 892                 sumdphi += temp2 - temp1;
 893                 phi1 = phi2;
 894         }
 895         if (sumdphi < 0.E0) fC = -fC;
 896         Double_t diff, angle = TMath::ATan2(fYC, fXC);
 897         if (fC < 0.E0)
 898                 fG0 = angle + 0.5 * TMath::Pi();
 899         else
 900                 fG0 = angle - 0.5 * TMath::Pi();
 901         diff = angle - fG0;
 902
 903         Double_t d = TMath::Sqrt(fXC*fXC + fYC*fYC) - fR;
 904         if (fC >= 0.E0)
 905                 fDt = d;
 906         else
 907                 fDt = -d;
 908
 909         Int_t nn = 6;
 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++) {
 913                 p = fPoint[j];
 914                 if (!p) break;
 915                 s[j] = ArgZ(p->GetR2());
 916                 if (s[j] > 100.0) return kFALSE;
 917                 z[j] = p->Z();
 918                 s[j] = asin(s[j]) / halfC;
 919                 ws[j] = 1.0 / (p->ErrZ()*p->ErrZ());
 920         }
 921
 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];
 929                 sumw += ws[i];
 930         }
 931         sums2 /= sumw;
 932         sumz /= sumw;
 933         sums /= sumw;
 934         sumsz /= sumw;
 935         d = sums2 - sums*sums;
 936
 937         fDz = (sums2*sumz - sums*sumsz) / d;
 938         fTanL = (sumsz - sums*sumz) / d;
 939
 940         delete [] s;
 941         delete [] z;
 942         delete [] ws;
 943
 944         return kTRUE;
 945 }
 946 //
 947 //
 948 //
 949 void AliITSNeuralTrack::PrintState(Bool_t matrix)
 950 {
 951 // Prints the state vector values.
 952 // The argument switches on or off the printing of the covariance matrix.
 953
 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";
 965         if (matrix) {
 966                 cout << "\nCovariance Matrix: " << endl;
 967                 fMatrix.Print();
 968         }
 969         cout << "Actual square chi = " << fChi2;
 970 }
 971 //
 972 //
 973 //
 974 Double_t AliITSNeuralTrack::GetDz() const
 975 {
 976 //      Double_t argZ = ArgZ(fStateR);
 977 //      if (argZ > 9.9) {
 978 //              Error("GetDz", "Too large value: %g", argZ);
 979 //              return 0.0;
 980 //      }
 981 //      fDz = fStateZ - (2.0 * fTanL / fC) * asin(argZ);
 982         return fDz;
 983 }
 984 //
 985 //
 986 //
 987 Double_t AliITSNeuralTrack::GetGamma() const
 988 {
 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);
 993 //              return kFALSE;
 994 //      }
 995 //      fG0 = fStatePhi - asin(argPhi);
 996         return fG0;
 997 }
 998 //
 999 //
1000 //
1001 Double_t AliITSNeuralTrack::GetPhi(Double_t r) const
1002 {
1003 // Gives the value of azymuthal coordinate in the helix
1004 // as a function of cylindric radius
1005
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(); }
1011         return arg;
1012 }
1013 //
1014 //
1015 //
1016 Double_t AliITSNeuralTrack::GetZ(Double_t r) const
1017 {
1018 // gives the value of Z in the helix
1019 // as a function of cylindric radius
1020
1021         Double_t arg = ArgZ(r);
1022         if (arg > 0.9) return 0.0;
1023         return fDz + fTanL * asin(arg) / fC;
1024 }
1025 //
1026 //
1027 //
1028 Double_t AliITSNeuralTrack::GetdEdX()
1029 {
1030 // total energy loss of the track
1031
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);
1038         }
1039         q[0] /= 280.;
1040         q[1] /= 280.;
1041         q[2] /= 38.;
1042         q[3] /= 38.;
1043         do {
1044                 swap = 0;
1045                 for (i = 0; i < 3; i++) {
1046                         if (q[i] <= q[i + 1]) continue;
1047                         Double_t tmp = q[i];
1048                         q[i] = q[i + 1];
1049                         q[i+1] = tmp;
1050                         swap++;
1051                 }
1052         } while(swap);
1053         if(q[0] < 0.) {
1054                 q[0] = q[1];
1055                 q[1] = q[2];
1056                 q[2] = q[3];
1057                 q[3] = -1.;
1058         }
1059         dedx = (q[0] + q[1]) / 2.;
1060         return dedx;
1061 }
1062 //
1063 //
1064 //
1065 void AliITSNeuralTrack::SetVertex(Double_t *pos, Double_t *err)
1066 {
1067         // Stores vertex data
1068
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);
1079         fVertex.SetUser(1);
1080 }
1081 //
1082 //
1083 //
1084 AliITSIOTrack* AliITSNeuralTrack::ExportIOtrack(Int_t min)
1085 {
1086 // Exports an object in the standard format for reconstructed tracks
1087
1088         Int_t layer = 0;
1089         AliITSIOTrack *track = new AliITSIOTrack;
1090
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));
1097
1098         // labels
1099         track->SetLabel(IsGood(min) ? fLabel : -fLabel);
1100         track->SetTPCLabel(-1);
1101
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());
1107                 }
1108         }
1109
1110         // state vector
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);
1118         track->SetDz(fDz);
1119
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);
1127
1128         // PID
1129         track->SetPid(fPDG);
1130         track->SetMass(fMass);
1131
1132         return track;
1133 }
1134 //
1135 //
1136 //====================================================================================
1137 //============================       PRIVATE METHODS      ============================
1138 //====================================================================================
1139 //
1140 //
1141 Double_t AliITSNeuralTrack::ArgPhi(Double_t r) const
1142 {
1143         // calculates the expression ((1/2)Cr + (1 + (1/2)CD) D/r) / (1 + CD)
1144
1145         Double_t arg, num, den;
1146         num = (0.5 * fC * r) + (1. + (0.5 * fC * fDt)) * (fDt / r);
1147         den = 1. + fC * fDt;
1148         if (den == 0.) {
1149                 Error("ArgPhi", "Denominator = 0!");
1150                 return 10.0;
1151         }
1152         arg = num / den;
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);
1156         return 10.0;
1157 }
1158 //
1159 //
1160 //
1161 Double_t AliITSNeuralTrack::ArgZ(Double_t r) const
1162 {
1163         // calculates the expression (1/2)C * sqrt( (r^2 - Dt^2) / (1 + CD) )
1164
1165         Double_t arg;
1166         arg = (r * r - fDt * fDt) / (1. + fC * fDt);
1167         if (arg < 0.) {
1168                 if (TMath::Abs(arg) < 1.E-6) arg = 0.;
1169                 else {
1170                         Error("ArgZ", "Square root argument error: %17.15g < 0", arg);
1171                         return 10.;
1172                 }
1173         }
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);
1178         return 10.0;
1179 }
1180 //
1181 //
1182 //
1183 Double_t AliITSNeuralTrack::ArgB(Double_t r) const
1184 {
1185 // UTILITY FUNCTION
1186
1187         Double_t arg;
1188         arg = (r*r - fDt*fDt);
1189         arg /= (r*(1.+ fC*fDt)*(1.+ fC*fDt));
1190         return arg;
1191 }
1192 //
1193 //
1194 //
1195 Double_t AliITSNeuralTrack::ArgC(Double_t r) const
1196 {
1197 // UTILITY FUNCTION
1198
1199         Double_t arg;
1200         arg = (1./r - fC * ArgPhi(r));
1201         arg /= 1.+ fC*fDt;
1202         return arg;
1203 }