ITS/AliITSNeuralTrack.cxx

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