ITS/AliITSneuralTracker.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: Alberto Pulvirenti.                                            *
   5  *                                                                        *
   6  * Permission to use, copy, modify and distribute this software and its   *
   7  * documentation strictly for non-commercial purposes is hereby granted   *
   8  * without fee, provided that the above copyright notice appears in all   *
   9  * copies and that both the copyright notice and this permission notice   *
  10  * appear in the supporting documentation. The authors make no claims     *
  11  * about the suitability of this software for any purpose. It is          *
  12  * provided "as is" without express or implied warranty.                  *
  13  *                                                                        *
  14  * track finder for ITS stand-alone with neural network algorithm         *
  15  * This class defines a Hopfield MFT neural network simulation            *
  16  * which reads all recpoints from an event and produces a tree with       *
  17  * the points belonging to recognized tracks                              *
  18  * the TTree obtained as the output file will be saved in a .root file    *
  19  * and it must be read by groups of six entries at a time                 *
  20  **************************************************************************/
  21
  22 #include <Riostream.h>
  23 #include <Riostream.h>
  24 #include <stdlib.h>
  25
  26 #include <TROOT.h>
  27 #include <TFile.h>
  28 #include <TObjArray.h>
  29 #include <TTree.h>
  30 #include <TMath.h>
  31 #include <TRandom.h>
  32 #include <TVector3.h>
  33
  34 #include "AliRun.h"
  35 #include "AliITS.h"
  36 #include "AliITSgeom.h"
  37 #include "AliITSgeomMatrix.h"
  38 #include "AliITSRecPoint.h"
  39 #include "AliITSglobalRecPoint.h"
  40 #include "AliITSneuralTrack.h"
  41 #include "AliITSneuralTracker.h"
  42
  43 ClassImp(AliITSneuralTracker)
  44
  45
  46 // ====================================================================================================
  47 // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> AliITSneuralTracker <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  48 // ====================================================================================================
  49
  50 AliITSneuralTracker::AliITSneuralTracker()
  51 {
  52         //
  53         // Default constructor (not for use)
  54         //
  55         Int_t i = 0;
  56         fCurvCut = 0;
  57         for (i = 0; i < 6; i++) {
  58                 fPoints[i] = 0;
  59                 if (i < 5) fNeurons[i] = 0;
  60         }
  61         fTracks = 0;
  62 }
  63 //---------------------------------------------------------------------------------------------------------
  64 AliITSneuralTracker::AliITSneuralTracker (Int_t nsecs, Int_t ncurv, Double_t *curv, Double_t theta)
  65 {
  66         // Argumented constructor
  67         // ----------------------
  68         // gets the number of azymuthal sectors,
  69         // the number of curvature cut steps and their cut values,
  70         // and then the polar angle cut.
  71         // Be careful not to put 'ncurv' greater than the dimension of the
  72         // 'curv' array (--> segmentation fault)
  73
  74         fSectorNum   = nsecs;
  75         fSectorWidth = 2.0 * TMath::Pi()/(Double_t)nsecs;
  76
  77         fCurvNum = ncurv;
  78         fCurvCut = new Double_t[ncurv];
  79         Int_t i;
  80         for (i = 0; i < ncurv; i++) fCurvCut[i] = curv[i];
  81
  82         fThetaCut = theta;
  83         fThetaNum = (Int_t)(TMath::Pi() / theta) + 1;
  84         if (fThetaNum < 1) fThetaNum = 1;
  85
  86         cout << "\nClass AliITSneuralTracker invoked with " << fSectorNum << " phi secs.";
  87         cout << " and " << fThetaNum << " theta secs.\n" << endl;
  88
  89         Int_t k, j;
  90         for (k = 0; k < 6; k++) {
  91                 fPoints[k] = new TObjArray**[fSectorNum];
  92                 for (i = 0; i < fSectorNum; i++) {
  93                         fPoints[k][i] = new TObjArray*[fThetaNum];
  94                         for (j = 0; j < fThetaNum; j++) fPoints[k][i][j] = new TObjArray(0);
  95                 }
  96                 if (k < 6) fNeurons[k] = new TObjArray(0);
  97         }
  98
  99         fTracks = new TObjArray(1);
 100 }
 101 //---------------------------------------------------------------------------------------------------------
 102 AliITSneuralTracker::~AliITSneuralTracker()
 103 {
 104         // Destructor
 105         // ----------
 106         // clears the TObjArrays and all heap objects
 107
 108         Int_t lay, sec, k;
 109         delete [] fCurvCut;
 110
 111         for (lay = 0; lay < 6; lay++) {
 112                 for (sec = 0; sec < fSectorNum; sec++) {
 113                         for (k = 0; k < fThetaNum; k++) {
 114                                 fPoints[lay][sec][k]->Delete();
 115                         }
 116                 }
 117                 delete [] fPoints[lay];
 118                 if (lay < 5) {
 119                         fNeurons[lay]->Delete();
 120                         delete fNeurons[lay];
 121                 }
 122         }
 123
 124         fTracks->Delete();
 125 }
 126 //---------------------------------------------------------------------------------------------------------
 127 Int_t AliITSneuralTracker::ReadFile(const Text_t *fname, Int_t evnum)
 128 {
 129         // File reader
 130         // -----------
 131         // Reads a ROOT file in order to retrieve recpoints.
 132         // the 'evnum' argument can have two meanings:
 133         // if it is >= 0 it represents the event number to retrieve the correct gAlice
 134         // if it is < 0, instead, this means that the open file contains
 135         // a 'TreeP' tree, with all points which remained unused
 136         // after the Kalman tracking procedure (produced via the ITSafterKalman.C macro)
 137         // the return value is the number of collected points
 138         // (then, if errors occur, the method returns 0)
 139
 140         TFile *file = new TFile(fname, "read");
 141         TTree *tree = 0;
 142         Int_t i, nentries, total, sector, mesh;
 143
 144         if (evnum < 0) {
 145                 tree = (TTree*)file->Get("TreeP0");
 146                 if (!tree) {
 147                         Error("", "Specifying evnum < 0, I expect to find a 'TreeP' tree in the file, but there isn't any");
 148                         return 0;
 149                 }
 150                 AliITSglobalRecPoint *p = 0;
 151                 tree->SetBranchAddress("Points", &p);
 152                 nentries = (Int_t)tree->GetEntries();
 153                 total = 0;
 154                 for (i = 0; i < nentries; i++) {
 155                         tree->GetEntry(i);
 156                         AliITSglobalRecPoint *q = new AliITSglobalRecPoint(*p);
 157                         sector = (Int_t)(q->fPhi / fSectorWidth);
 158                         mesh = (Int_t)(q->fTheta / fThetaCut);
 159                         q->fUsed = 0;
 160                         fPoints[q->fLayer][sector][mesh]->AddLast(q);
 161                         total++;
 162                 }
 163                 return total;
 164         }
 165
 166         if (gAlice) { delete gAlice; gAlice = 0; }
 167         gAlice = (AliRun*)file->Get("gAlice");
 168         if (!gAlice) {
 169                 Error("", "gAlice is NULL, maybe wrong filename...");
 170                 return 0;
 171         }
 172
 173         Int_t nparticles = (Int_t)gAlice->GetEvent(evnum);
 174         if (!nparticles) {
 175                 Error("", "No particles found");
 176                 return 0;
 177         }
 178
 179         AliITS *its = (AliITS*)gAlice->GetModule("ITS");
 180         if (!its) {
 181                 Error("", "No ITS found");
 182                 return 0;
 183         }
 184
 185         AliITSgeom *geom = (AliITSgeom*)its->GetITSgeom();
 186         if (!geom) {
 187                 Error("", "AliITSgeom is NULL");
 188                 return 0;
 189         }
 190
 191         tree = gAlice->TreeR();
 192         if (!tree) {
 193                 Error("", "TreeR() returned NULL");
 194                 return 0;
 195         }
 196
 197         nentries = (Int_t)tree->GetEntries();
 198         if (!nentries) {
 199                 Error("", "TreeR is empty");
 200                 return 0;
 201         }
 202
 203         // Reading objects initialization
 204         Int_t k, entry, lay, lad, det;
 205         Double_t l[3], g[3], ll[3][3], gg[3][3];
 206
 207         TObjArray *recsArray = its->RecPoints();
 208         AliITSRecPoint *recp;
 209         AliITSglobalRecPoint *pnt = 0;
 210
 211         total = 0;
 212         for(entry = 0; entry < nentries; entry++) {
 213                 if (entry > geom->GetLastSSD()) continue;
 214                 its->ResetRecPoints();
 215                 tree->GetEvent(entry);
 216                 TObjArrayIter recs(recsArray);
 217                 AliITSgeomMatrix *gm = geom->GetGeomMatrix(entry);
 218                 geom->GetModuleId(entry, lay, lad, det);
 219                 lay--;
 220                 while ((recp = (AliITSRecPoint*)recs.Next())) {
 221                         // conversion to global coordinate system
 222                         for (i = 0; i < 3; i++) {
 223                                 l[i] = g[i] = 0.0;
 224                                 for (k = 0; k < 3; k++) {
 225                                         ll[i][k] = 0.0;
 226                                         gg[i][k] = 0.0;
 227                                 }
 228                         }
 229                         // local to global conversions of coords
 230                         l[0] = recp->fX;
 231                         l[1] = 0.0;
 232                         l[2] = recp->fZ;
 233                         gm->LtoGPosition(l, g);
 234                         // local to global conversions of sigmas
 235                         ll[0][0] = recp->fSigmaX2;
 236                         ll[2][2] = recp->fSigmaZ2;
 237                         gm->LtoGPositionError(ll, gg);
 238                         // instantiation of a new global rec-point
 239                         pnt = new AliITSglobalRecPoint(g[0], g[1], g[2], gg[0][0], gg[1][1], gg[2][2], lay);
 240                         // copy of owner track labels
 241                         for (i = 0; i < 3; i++) pnt->fLabel[i] = recp->fTracks[i];
 242                         sector = (Int_t)(pnt->fPhi / fSectorWidth);
 243                         mesh = (Int_t)(pnt->fTheta / fThetaCut);
 244                         pnt->fUsed = 0;
 245                         fPoints[lay][sector][mesh]->AddLast(pnt);
 246                         total++;
 247                 }
 248         }
 249
 250         return total;
 251 }
 252 //---------------------------------------------------------------------------------------------------------
 253 void AliITSneuralTracker::Go(const char* filesave, Bool_t flag)
 254 {
 255         // Global working method
 256         // ---------------------
 257         // It's the method which calls all other ones,
 258         // in order to analyze the whole ITS points ensemble
 259         // It's the only method to use, besides the parameter setters,
 260         // so, all other methods are defined as private
 261         // (see header file)
 262         // the flag meaning is explained in the 'Initialize' method
 263
 264         Int_t i;
 265         for (i = 0; i < fSectorNum; i++) DoSector(i, flag);
 266         cout << endl << "Saving results in file " << filesave << "..." << flush;
 267         TFile *f = new TFile(filesave, "recreate");
 268         fTracks->Write();
 269         cout << "done" << endl;
 270         f->Close();
 271 }
 272
 273 //=========================================================================================================
 274 // * * *   P R I V A T E   S E C T I O N   * * *
 275 //=========================================================================================================
 276
 277 Int_t AliITSneuralTracker::Initialize(Int_t secnum, Int_t curvnum, Bool_t flag)
 278 {
 279         // Network initializer
 280         // -------------------
 281         // Fills the neuron arrays, following the cut criteria for the selected step
 282         // (secnum = sector to analyze, curvnum = curvature cut step to use)
 283         // It also sets the initial random activation.
 284         // In order to avoid a large number of 'new' operations, all existing neurons
 285         // are reset and initialized with the new values, and are created new unit only if
 286         // it turns out to be necessary
 287         // the 'flag' argument is used to decide if the lower edge in the intevral
 288         // of the accepted curvatures is given by zero (kFALSE) or by the preceding used cut (kTRUE)
 289         // (if this is the first step, anyway, the minimum is always zero)
 290
 291         Int_t l, m, k, neurons = 0;
 292
 293         for (l = 0; l < 5; l++) fNeurons[l]->Delete();
 294
 295         AliITSneuron *unit = 0;
 296         Double_t abscurv, max = fCurvCut[curvnum], min = 0.0;
 297         if (flag && curvnum > 0) min = fCurvCut[curvnum - 1];
 298         AliITSglobalRecPoint *inner = 0, *outer = 0;
 299         for (l = 0; l < 5; l++) {
 300                 for (m = 0; m < fThetaNum; m++) {
 301                         TObjArrayIter inners(fPoints[l][secnum][m]);
 302                         while ((inner = (AliITSglobalRecPoint*)inners.Next())) {
 303                                 if (inner->fUsed > 0) continue; // points can't be recycled
 304                                 for (k = m-1; k <= m+1; k++) {
 305                                         if (k < 0 || k >= fThetaNum) continue; // to avoid seg faults
 306                                         TObjArrayIter outers(fPoints[l+1][secnum][k]);
 307                                         while ((outer = (AliITSglobalRecPoint*)outers.Next())) {
 308                                                 if (outer->fUsed > 0) continue;  // points can't be recycled
 309                                                 if (inner->DTheta(outer) > fThetaCut) continue;
 310                                                 unit = new AliITSneuron;
 311                                                 unit->fInner = inner;
 312                                                 unit->fOuter = outer;
 313                                                 CalcParams(unit);
 314                                                 abscurv = TMath::Abs(unit->fCurv);
 315                                                 if (unit->fDiff > fDiff || abscurv < min || abscurv > max) {
 316                                                         delete unit;
 317                                                         continue;
 318                                                 }
 319                                                 unit->fActivation = gRandom->Rndm() * (fMax - fMin) + fMin;
 320                                                 fNeurons[l]->AddLast(unit);
 321                                                 neurons++;
 322                                         } // end loop on candidates for outer point for neurons
 323                                 }
 324                         } // end loop on candidates inner points for neuron
 325                 }
 326         } // end loop on layers
 327         return neurons;
 328 }
 329 //---------------------------------------------------------------------------------------------------------
 330 Bool_t AliITSneuralTracker::Update()
 331 {
 332         // Updating cycle
 333         // --------------
 334         // Performs an updating cycle, by summing all the
 335         // gain and cost contribution for each neuron and
 336         // updating its activation.
 337         // An asynchronous method is followed, with the order
 338         // according that the neurons have been created
 339         // Time wasting is avoided by dividing the neurons
 340         // into different arrays, depending on the layer of their inner point.
 341         // The return value answers the question: "The network has stabilized?"
 342
 343         Int_t j, l;
 344         Double_t actOld, actNew, argNew, totDiff = 0.0;
 345         Double_t sumInh, sumExc, units = 0.0;
 346         AliITSneuron *me = 0, *it = 0;
 347         for (l = 0; l < 5; l++) {
 348                 TObjArrayIter meIter(fNeurons[l]);
 349                 while ((me = (AliITSneuron*)meIter.Next())) {
 350                         units++;
 351                         TObjArrayIter inhIter(fNeurons[l]);
 352
 353                         // inhibition (search only for neurons starting in the same layer)
 354                         sumInh = 0.0;
 355                         while((it = (AliITSneuron*)inhIter.Next())) {
 356                                 if (it->fOuter == me->fOuter || it->fInner == me->fInner)
 357                                         sumInh += it->fActivation;
 358                         }
 359
 360                         // excitation (search only for neurons
 361                         // which start in the previous or next layers)
 362                         sumExc = 0.0;
 363                         for (j = l - 1; j <= l + 1; j += 2) {
 364                                 if (j < 0 || j >= 5) continue;
 365                                 TObjArrayIter itIter(fNeurons[j]);
 366                                 while ((it = (AliITSneuron*)itIter.Next())) {
 367                                         if (it->fInner == me->fOuter || it->fOuter == me->fInner && it->fCurv * me->fCurv > 0.0) {
 368                                                 sumExc += Weight(me, it) * it->fActivation;
 369                                         }
 370                                 }
 371                         } // end search for excitories
 372
 373                         actOld = me->fActivation;
 374                         argNew = fRatio * sumExc - sumInh;
 375                         actNew = 1.0 / (1.0 + TMath::Exp(-argNew / fTemp));
 376                         me->fActivation = actNew;
 377                         // calculation the relative activation variation
 378                         // (for stabilization check)
 379                         totDiff += TMath::Abs((actNew - actOld) / actOld);
 380                 } // end loop over 'me' (updated neuron)
 381         } // end loop over layers
 382         totDiff /= units;
 383
 384         return (totDiff < fVar);
 385 }
 386 //---------------------------------------------------------------------------------------------------------
 387 Int_t AliITSneuralTracker::Save()
 388 {
 389         // Tracki saving method
 390         // --------------------
 391         // Stores the recognized tracks into the TObjArray 'fTracks'
 392         // but doesn't save it into a file until the whole work has ended
 393
 394         Int_t l;
 395         Double_t test = 0.5;
 396
 397         // Before all, for each neuron is found the best active sequences
 398         // among the incoming and outgoing other units
 399         AliITSneuron *main, *fwd, *back, *start;
 400         for (l = 0; l < 5; l++) {
 401                 TObjArrayIter mainIter(fNeurons[l]);
 402                 while((main = (AliITSneuron*)mainIter.Next())) {
 403                         if (l < 4) {
 404                                 TObjArrayIter fwdIter(fNeurons[l+1]);
 405                                 test = 0.5;
 406                                 while ((fwd = (AliITSneuron*)fwdIter.Next())) {
 407                                         if (main->fOuter == fwd->fInner && fwd->fActivation > test) {
 408                                                 test = fwd->fActivation;
 409                                                 main->fNext = fwd;
 410                                         }
 411                                 };
 412                         }
 413                         if (l > 0) {
 414                                 TObjArrayIter backIter(fNeurons[l-1]);
 415                                 test = 0.5;
 416                                 while ((back = (AliITSneuron*)backIter.Next())) {
 417                                         if (main->fInner == back->fOuter && back->fActivation > test) {
 418                                                 test = back->fActivation;
 419                                                 main->fPrev = back;
 420                                         }
 421                                 };
 422                         }
 423                 }
 424         }
 425
 426         // Then, are resolved the mismatches in the chains found above
 427         // (when unit->next->prev != unit or unit->prev->next != unit)
 428         for (l = 0; l < 5; l++) {
 429                 TObjArrayIter mainIter(fNeurons[l]);
 430                 while((main = (AliITSneuron*)mainIter.Next())) {
 431                         fwd  = main->fNext;
 432                         back = main->fPrev;
 433                         if (fwd  != 0 && fwd->fPrev  != main) main->fNext = 0;
 434                         if (back != 0 && back->fNext != main) main->fPrev = 0;
 435                 }
 436         }
 437
 438         Int_t pos, neg, incr, decr, ntracks;
 439         AliITSneuralTrack *trk = 0;
 440
 441         TObjArrayIter startIter(fNeurons[0]);
 442         for (;;) {
 443                 start = (AliITSneuron*)startIter.Next();
 444                 if (!start) break;
 445                 Int_t pts = 0;
 446                 // with the chain above defined, a track can be followed
 447                 // like a linked list structure
 448                 for (main = start; main; main = main->fNext) pts++;
 449                 // the count will be  > 5 only if the track contains 6 points
 450                 // (what we want)
 451                 if (pts < 5) continue;
 452                 // track storage
 453                 trk = new AliITSneuralTrack;
 454                 trk->CopyPoint(start->fInner);
 455                 for (main = start; main; main = main->fNext) {
 456                         //main->fInner->fUsed = kTRUE;
 457                         //main->fOuter->fUsed = kTRUE;
 458                         trk->CopyPoint(main->fOuter);
 459                 }
 460                 trk->Kinks(pos, neg, incr, decr);
 461                 if (pos != 0 && neg != 0) {
 462                         cout << "pos, neg kinks = " << pos << " " << neg << endl;
 463                         continue;
 464                 }
 465                 if (incr != 0 && decr != 0) {
 466                         cout << "increment, decrements in delta phi = " << incr << " " << decr << endl;
 467                         continue;
 468                 }
 469                 for (main = start; main; main = main->fNext) {
 470                         main->fInner->fUsed++;
 471                         main->fOuter->fUsed++;
 472                 }
 473                 fTracks->AddLast(trk);
 474         }
 475         ntracks = (Int_t)fTracks->GetEntriesFast();
 476         return ntracks;
 477 }
 478 //---------------------------------------------------------------------------------------------------------
 479 void AliITSneuralTracker::DoSector(Int_t sect, Bool_t flag)
 480 {
 481         // Sector recognition
 482         // ------------------
 483         // This is a private method which works on a single sector
 484         // just for an easier readability of the code.
 485         // The 'flag' argument is needed to be given to the 'Initialize' method (see it)
 486
 487         Int_t curvnum, nTracks = 0, nUnits;
 488         Bool_t isStable = kFALSE;
 489         for (curvnum = 0; curvnum < fCurvNum; curvnum++) {
 490                 cout << "\rSector " << sect << ":" << ((curvnum+1)*100)/fCurvNum << "%" << flush;
 491                 nUnits = Initialize(sect, curvnum, flag);
 492                 if (!nUnits) continue;
 493                 do {
 494                         isStable = Update();
 495                 } while (!isStable);
 496                 nTracks = Save();
 497         }
 498         cout << "\rTracks stored after sector " << sect << ": " << nTracks << endl;
 499 }
 500 //---------------------------------------------------------------------------------------------------------
 501 Int_t AliITSneuralTracker::CountExits(AliITSneuron *n)
 502 {
 503         // Counter for neurons which enter in 'n' tail
 504
 505         Int_t count = 0, l = n->fOuter->fLayer;
 506         if (l == 5) return 0;
 507
 508         AliITSneuron *test = 0;
 509         TObjArrayIter iter(fNeurons[l]);
 510         while ((test = (AliITSneuron*)iter.Next())) {
 511                 if (test->fInner == n->fOuter) count++;
 512         }
 513         return count;
 514 }
 515 //---------------------------------------------------------------------------------------------------------
 516 Int_t AliITSneuralTracker::CountEnters(AliITSneuron *n)
 517 {
 518         // Counter for neurons which exit from 'n' head
 519
 520         Int_t count = 0, l = n->fInner->fLayer;
 521         if (l == 0) return 0;
 522
 523         AliITSneuron *test = 0;
 524         TObjArrayIter iter(fNeurons[l-1]);
 525         while ((test = (AliITSneuron*)iter.Next())) {
 526                 if (test->fOuter == n->fInner) count++;
 527         }
 528         return count;
 529 }
 530 //---------------------------------------------------------------------------------------------------------
 531 Bool_t AliITSneuralTracker::ResetNeuron(AliITSneuron *n)
 532 {
 533         // A method which resets the neuron's pointers to zero
 534
 535         if (!n) return kFALSE;
 536         n->fNext = 0;
 537         n->fPrev = 0;
 538         n->fInner = 0;
 539         n->fOuter = 0;
 540         // the other datamembers don't need to be reinitialized
 541         return kTRUE;
 542 }
 543 //---------------------------------------------------------------------------------------------------------
 544 Bool_t AliITSneuralTracker::SetEdge(AliITSneuron *n, AliITSglobalRecPoint *p, const char what)
 545 {
 546         // Sets the indicated edge point to the neuron
 547         // (if the arguments are suitable ones...)
 548
 549         if (!n || !p || (what != 'i' && what != 'I' && what != 'o' && what != 'O')) return kFALSE;
 550         if (what == 'i' || what == 'I')
 551                 n->fInner = p;
 552         else
 553                 n->fOuter = p;
 554         return kTRUE;
 555 }
 556 //---------------------------------------------------------------------------------------------------------
 557 Bool_t AliITSneuralTracker::CalcParams(AliITSneuron *n)
 558 {
 559         // This method evaluates the helix parameters from the edged of a neuron
 560         // (curvature, phase parameter, tangent of lambda angle)
 561         // if the p[arameters assume too strange (and unphysical) values, the
 562         // method return kFALSE, else it return kTRUE;
 563
 564         if (!n) return kFALSE;
 565
 566         Double_t sign(1.0), det, rc = 0.0;
 567         Double_t tanl1, tanl2, l1, l2;
 568         // changing the reference frame into the one centered in the vertex
 569         Double_t x1 = n->fInner->fGX - fVPos[0];
 570         Double_t y1 = n->fInner->fGY - fVPos[1];
 571         Double_t z1 = n->fInner->fGZ - fVPos[2];
 572         Double_t r1 = x1 * x1 + y1 * y1;
 573         Double_t x2 = n->fOuter->fGX - fVPos[0];
 574         Double_t y2 = n->fOuter->fGY - fVPos[1];
 575         Double_t z2 = n->fOuter->fGZ - fVPos[2];
 576         Double_t r2 = x2 * x2 + y2 * y2;
 577         // initialization of the XY-plane data (curvature and centre)
 578         // (will remain these values if we encounter errors)
 579         n->fCurv = n->fCX = n->fCY = n->fTanL = n->fPhase = 0.0;
 580         n->fLength = (n->fOuter->fGX - n->fInner->fGX) * (n->fOuter->fGX - n->fInner->fGX);
 581         n->fLength += (n->fOuter->fGY - n->fInner->fGY) * (n->fOuter->fGY - n->fInner->fGY);
 582         n->fLength += (n->fOuter->fGZ - n->fInner->fGZ) * (n->fOuter->fGZ - n->fInner->fGZ);
 583         n->fLength = TMath::Sqrt(n->fLength);
 584         // this determinant which is zero when the points are in a line
 585         det = 2.0 * (x1 * y2 - x2 * y1);
 586         if (det == 0.0)
 587                 return kFALSE; // it is difficult having perfectly aligned points --- it's more likely to be an error
 588         else {
 589                 n->fCX = (r1*y2 - r2*y1) / det;
 590                 n->fCY = (r2*x1 - r1*x2) / det;
 591                 rc = TMath::Sqrt(n->fCX * n->fCX + n->fCY * n->fCY);
 592                 r1 = TMath::Sqrt(r1);
 593                 r2 = TMath::Sqrt(r2);
 594                 sign = (n->fOuter->fPhi >= n->fInner->fPhi) ? 1.0 : -1.0;
 595                 if (rc > 0.0) {
 596                         n->fCurv = sign / rc;
 597                         n->fPhase = TMath::ATan2(-n->fCX, -n->fCY);
 598                         if (n->fPhase < 0.0) n->fPhase += 2.0 * TMath::Pi(); // angles in 0 --> 2pi
 599                         if (r1 <= 2.0 * rc && r2 <= 2.0 * rc) {
 600                                 l1 = 2.0 * rc * TMath::ASin(r1 / (2.0 * rc));
 601                                 l2 = 2.0 * rc * TMath::ASin(r2 / (2.0 * rc));
 602                                 tanl1 = z1 / l1;
 603                                 tanl2 = z2 / l2;
 604                                 n->fDiff = TMath::Abs(tanl1 - tanl2);
 605                                 n->fTanL = 0.5 * (tanl1 + tanl2);
 606                         }
 607                         else {
 608                                 n->fTanL = 0.0;
 609                                 n->fDiff = 100.0;
 610                                 return kFALSE; // it' even more difficult that the arguments above aren't acceptable for arcsine
 611                         }
 612                 }
 613                 else
 614                         return kFALSE;
 615         }
 616         return kTRUE;
 617 }
 618 //---------------------------------------------------------------------------------------------------------
 619 Double_t AliITSneuralTracker::Angle(AliITSneuron *n, AliITSneuron *m)
 620 {
 621         // calculates the angle between two segments
 622         // but doesn't check if they have a common point.
 623         // The calculation is made by the ratio of their scalar product
 624         // to the product of their magnitudes
 625
 626         Double_t x = (m->fOuter->fGX - m->fInner->fGX) * (n->fOuter->fGX - n->fInner->fGX);
 627         Double_t y = (m->fOuter->fGY - m->fInner->fGY) * (n->fOuter->fGY - n->fInner->fGY);
 628         Double_t z = (m->fOuter->fGZ - m->fInner->fGZ) * (n->fOuter->fGZ - n->fInner->fGZ);
 629
 630         Double_t cosine = (x + y + z) / (n->fLength * m->fLength);
 631
 632         if (cosine > 1.0 || cosine < -1.0) {
 633                 Warning("AliITSneuronV1::Angle", "Strange value of cosine");
 634                 return 0.0;
 635         }
 636
 637         return TMath::ACos(cosine);
 638 }
 639 //---------------------------------------------------------------------------------------------------------
 640 Double_t AliITSneuralTracker::Weight(AliITSneuron *n, AliITSneuron *m)
 641 {
 642         // calculation of the excitoy weight
 643
 644 //      Double_t v1 = TMath::Abs(fTanL - n->fTanL);
 645 //      Double_t v2 = TMath::Abs(TMath::Abs(fCurv) - TMath::Abs(n->fCurv));
 646 //      Double_t v3 = TMath::Abs(fPhase - n->fPhase);
 647 //      Double_t q = (v1 + v2 + v3) / 3.0;
 648 //      Double_t a = 1.0; // 0.02 / (fDiff + n->fDiff);
 649         Double_t b = 1.0 - TMath::Sin(Angle(n, m));
 650         return TMath::Power(b, fExp);
 651 //      return TMath::Exp(q / par1);
 652
 653 }
 654 //---------------------------------------------------------------------------------------------------------
 655 // ====================================================================================================
 656 // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> AliITSneuron <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 657 // ====================================================================================================
 658 ClassImp(AliITSneuron)
 659 AliITSneuron::AliITSneuron()
 660 {
 661         // Default constructor
 662         // -------------------
 663         // sets all parameters to default (but unuseful)
 664         // (in particular all pointers are set to zero)
 665         // There is no alternative constructor, because
 666         // this class should not be used by itself, but
 667         // it is only necessary to allow the neural tracker to work
 668         // In fact, all methods which operate on neuron's datamembers
 669         // are stored in the AliITSneuralTracker class anyway
 670
 671         fActivation = 0.0;
 672
 673         fCurv = 0.0;
 674         fCX = 0.0;
 675         fCY = 0.0;
 676         fTanL = 0.0;
 677         fPhase = 0.0;
 678         fDiff = 0.0;
 679         fLength = 0.0;
 680
 681         fNext = fPrev = 0;
 682         fInner = fOuter = 0;
 683 }
 684 //---------------------------------------------------------------------------------------------------------
 685 AliITSneuron::~AliITSneuron()
 686 {
 687         // Destructor
 688         // ----------
 689         // does nothing, because there is no need to affect parameters,
 690         // while the AliITSglobalRecPoint pointers can't be deleted,
 691         // and the AliITSneuron pointers belong to a TObjArray which will
 692         // be cleared when necessary...
 693 }
 694