ITS/AliITStrackerANN.cxx

   1 // ---------------------------------------------------------------------------------
   2 // AliITStrackerANN
   3 // ---------------------------------------------------------------------------------
   4 // Neural Network-based algorithm for track reconstruction in the ALICE ITS.
   5 // The class contains all that is necessary in order to perform the tracking
   6 // using the ITS clusters in the V2 format.
   7 // The class is organized with some embedded objects which serve for
   8 // data management, and all setters and getters for all parameters.
   9 // The usage of the class from a tracking macro should be based essentially
  10 // in the initialization (constructor) and the call of a method which
  11 // performs in the right sequence all the operations.
  12 // Temporarily (and maybe definitively) all the intermediate steps are
  13 // public functions which could be called independently, but they do not
  14 // produce any result if all the preventively required steps have not been
  15 // performed.
  16 // ---------------------------------------------------------------------------------
  17 // Author: Alberto Pulvirenti (university of Catania)
  18 // Email : alberto.pulvirenti@ct.infn.it
  19 // ---------------------------------------------------------------------------------
  20
  21 #include <Riostream.h>
  22 #include <TMath.h>
  23 #include <TString.h>
  24 #include <TObjArray.h>
  25 #include <TVector3.h>
  26 #include <TFile.h>
  27 #include <TTree.h>
  28 #include <TRandom.h>
  29 #include <TMatrixD.h>
  30 #if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,2)
  31   #include <TMatrixDEigen.h>
  32 #endif
  33
  34 #include "AliITSgeom.h"
  35 #include "AliITStrackSA.h"
  36 #include "AliITSclusterV2.h"
  37 #include "AliITStrackV2.h"
  38
  39 #include "AliITStrackerANN.h"
  40
  41 const Double_t AliITStrackerANN::fgkPi     = 3.141592653; // pi
  42 const Double_t AliITStrackerANN::fgkHalfPi = 1.570796327; // pi / 2
  43 const Double_t AliITStrackerANN::fgkTwoPi  = 6.283185307; // 2 * pi
  44
  45 ClassImp(AliITStrackerANN)
  46
  47 //__________________________________________________________________________________
  48 AliITStrackerANN::AliITStrackerANN(const AliITSgeom *geom, Int_t msglev)
  49 : AliITStrackerV2(geom), fMsgLevel(msglev)
  50 {
  51 /**************************************************************************
  52
  53                 CONSTRUCTOR (standard-to-use version)
  54
  55                 Arguments:
  56                         1) the ITS geometry used in the generated event
  57                         2) the flag for log-messages writing
  58
  59                 The AliITSgeometry is used along the class,
  60                 in order to translate the local AliITSclusterV2 coordinates
  61                 into the Global reference frame, which is necessary for the
  62                 Neural Network algorithm to operate.
  63                 In case of serialized use, the log messages should be excluded,
  64                 in order to save real execution time.
  65
  66                 Operations:
  67                         - all pointer data members are initialized
  68                           (if possible, otherwise are set to NULL)
  69                         - all numeric data members are initialized to
  70                           values which allow the Neural Network to operate
  71                           even if nothing is externally set.
  72
  73                  NOTE: it is possible that tracking an event
  74                        with these default values results in a non-sense.
  75
  76  **************************************************************************/
  77
  78         Int_t i;
  79
  80         // Get ITS geometry
  81         fGeom = (AliITSgeom*)geom;
  82
  83         // Initialize the array of first module indexes per layer
  84         fNLayers = geom->GetNlayers();
  85         fFirstModInLayer = new Int_t[fNLayers + 1];
  86         for (i = 0; i < fNLayers; i++) {
  87                 fFirstModInLayer[i] = fGeom->GetModuleIndex(i + 1, 1, 1);
  88         }
  89         fFirstModInLayer[fNLayers] = geom->GetIndexMax();
  90
  91         // initialization: no curvature cut steps
  92         fCurvNum = 0;
  93         fCurvCut = 0;
  94
  95         // initialization: 4 sectors (one for each quadrant)
  96         fSectorNum   = 4;
  97         fSectorWidth = fgkHalfPi;
  98
  99         // initialization: theta offset of 20 degrees
 100         fPolarInterval = 20.0;
 101
 102         // initialization: array structure not defined
 103         fStructureOK = kFALSE;
 104
 105         // initialization: vertex in the origin
 106         fVertexX = 0.0;
 107         fVertexY = 0.0;
 108         fVertexZ = 0.0;
 109
 110         // initialization: uninitialized point array
 111         fNodes = 0;
 112
 113         // initialization: very large (unuseful) cut values
 114         Int_t ilayer;
 115         for (ilayer = 0; ilayer < 6; ilayer++) {
 116                 fThetaCut2D[ilayer] = TMath::Pi();
 117                 fThetaCut3D[ilayer] = TMath::Pi();
 118                 fHelixMatchCut[ilayer] = 1.0;
 119         }
 120
 121         // initialization: inictial activation range between 0.3 and 0.7
 122         fEdge1 = 0.3;
 123         fEdge2 = 0.7;
 124
 125         // initialization: neural network operation & weights
 126         fTemperature = 1.0;
 127         fStabThreshold = 0.001;
 128         fGain2CostRatio = 1.0;
 129         fExponent = 1.0;
 130         fActMinimum = 0.5;
 131
 132         // initialization: uninitialized array of neurons
 133         fNeurons = 0;
 134
 135         // initialization: uninitialized array of tracks
 136         fFoundTracks = 0;
 137 }
 138
 139 //__________________________________________________________________________________
 140 void AliITStrackerANN::SetCuts
 141 (Int_t ncurv, Double_t *curv, Double_t *theta2D, Double_t *theta3D, Double_t *helix)
 142
 143 /**************************************************************************
 144
 145         CUT SETTER
 146
 147         Allows for the definition of all kind of geometric cuts
 148         which have been studied in for the creation of a neuron
 149         from a pair of clusters C1 and C2 on consecutive layers.
 150         Neuron will be created only if the pair passes ALL of these cuts.
 151
 152         At the moment, we define 4 kinds of geometrical cuts:
 153                 a) cut on the difference of the polar 'theta' angle;
 154                 b) cut on the angle between origin->C1 and C1->C2 in space;
 155                 c) cut on the curvature of the circle passing
 156                    through C1, C2 and the primary vertex;
 157                 d) cut on heli-matching of the same three points.
 158
 159         Arguments:
 160                 1) the number of curvature cut steps
 161                 2) the array of curvature cuts for each step
 162                    (its dimension is given by the argument 1)
 163                 3) array of 5 cut values (one for each consecutive lauer pair)
 164                    related to cut a)
 165                 4) array of 5 cut values (one for each consecutive lauer pair)
 166                    related to cut b)
 167                 5) array of 5 cut values (one for each consecutive lauer pair)
 168                      related to cut c)
 169
 170         Operations:
 171                 - gets the values for each cut and stores them in data members
 172                 - in the case of curvature cuts, the cut array
 173                   (whose size is not fixed) is allocated
 174
 175         NOTE: in case that the user wants to set onyl ONE of the 4 cuts array,
 176               he can simply pass NULL arguments for other cuts, and (eventually)
 177                         a ZERO as the first argument (if curvature cuts have not to be set).
 178
 179         Anyway, all the cuts have to be set at least once.
 180
 181  **************************************************************************/
 182 {
 183         // counter
 184         Int_t i;
 185
 186         /*** Curvature cut setting ***/
 187
 188         // first of all, the curvature cuts are sorted in increasing order
 189         // (from the smallest to the largest one)
 190         Int_t *ind = new Int_t[ncurv];
 191         TMath::Sort(ncurv, curv, ind, kFALSE);
 192         // then, the curvature cut array is allocated and filled
 193         // (a message with the list of defined cuts can be optionally shown)
 194         fCurvCut = new Double_t[ncurv];
 195         if (fMsgLevel >= 1) cout << "Number of curvature cuts: " << ncurv << endl;
 196         for (i = 0; i < ncurv; i++) {
 197                 fCurvCut[i] = curv[ind[i]];
 198                 if (fMsgLevel >= 1) cout << " - " << fCurvCut[i] << endl;
 199         }
 200         fCurvNum = ncurv;
 201
 202         /*** 'Fixed' cuts setting ***/
 203
 204         // checks what cuts have to be set
 205         Bool_t doTheta2D = (theta2D != 0);
 206         Bool_t doTheta3D = (theta3D != 0);
 207         Bool_t doHelix = (helix != 0);
 208         // sets the cuts for all layer pairs
 209         for (i = 0; i < fNLayers - 1; i++) {
 210                 if (doTheta2D) fThetaCut2D[i] = theta2D[i];
 211                 if (doTheta3D) fThetaCut3D[i] = theta3D[i];
 212                 if (doHelix) fHelixMatchCut[i] = helix[i];
 213         }
 214         // if required, lists the cuts
 215         if (!fMsgLevel < 2) return;
 216         cout << "Theta 2D cuts: ";
 217         if (!doTheta2D) {
 218                 cout << "<not set>" << endl;
 219         }
 220         else {
 221                 cout << endl;
 222                 for (i = 0; i < fNLayers - 1; i++) {
 223                         cout << "For layers " << i << " --> " << i + 1;
 224                         cout << " cut = " << fThetaCut2D[i] << endl;
 225                 }
 226         }
 227         cout << "---" << endl;
 228         cout << "Theta 3D cuts: ";
 229         if (!doTheta3D) {
 230                 cout << "<not set>" << endl;
 231         }
 232         else {
 233                 cout << endl;
 234                 for (i = 0; i < fNLayers - 1; i++) {
 235                         cout << "For layers " << i << " --> " << i + 1;
 236                         cout << " cut = " << fThetaCut3D[i] << endl;
 237                 }
 238         }
 239         cout << "---" << endl;
 240         cout << "Helix-match cuts: ";
 241         if (!doHelix) {
 242                 cout << "<not set>" << endl;
 243         }
 244         else {
 245                 cout << endl;
 246                 for (i = 0; i < fNLayers - 1; i++) {
 247                         cout << "For layers " << i << " --> " << i + 1;
 248                         cout << " cut = " << fHelixMatchCut[i] << endl;
 249                 }
 250         }
 251         cout << "---" << endl;
 252 }
 253
 254 //__________________________________________________________________________________
 255 Bool_t AliITStrackerANN::GetGlobalXYZ
 256 (Int_t refIndex,
 257  Double_t &x, Double_t &y, Double_t &z, Double_t &ex2, Double_t &ey2, Double_t &ez2)
 258
 259 /**************************************************************************
 260
 261         LOCAL TO GLOBAL TRANSLATOR
 262
 263         Taking information from the ITS geometry stored in the class,
 264         gets a stored AliITScluster and calculates it coordinates
 265         and errors in the global reference frame.
 266         These values are stored in the variables,
 267         which are passed by reference.
 268
 269         Arguments:
 270                 1) reference index for the cluster to use
 271                 2) (by reference) place to store the global X-coord into
 272                 3) (by reference) place to store the global Y-coord into
 273                 4) (by reference) place to store the global Z-coord into
 274                 5) (by reference) place to store the global X-coord error into
 275                 6) (by reference) place to store the global Y-coord error into
 276                 7) (by reference) place to store the global Z-coord error into
 277
 278         Operations:
 279                 essentially, determines the ITS module index from the
 280                 detector index of the AliITSclusterV2 object, and extracts
 281                 the roto-translation from the ITS geometry, to convert
 282                 the local module coordinates into the global ones.
 283
 284         Return value:
 285                 - kFALSE if the given cluster index points to a non-existing
 286                 cluster, or if the layer number makes no sense (< 0 or > 6).
 287                 - otherwise, kTRUE (meaning a successful operation).
 288
 289  **************************************************************************/
 290 {
 291         // checks if the layer is correct
 292         Int_t ilayer = (refIndex & 0xf0000000) >> 28;
 293         if (ilayer < 0 || ilayer >= fNLayers) {
 294                 Error("GetGlobalXYZ", "Wrong layer number: %d [range: %d - %d]", ilayer, 0, fNLayers);
 295                 return kFALSE;
 296         }
 297         // checks if the referenced cluster exists and corresponds to the passed reference
 298         AliITSclusterV2 *refCluster = (AliITSclusterV2*) GetCluster(refIndex);
 299         if (!refCluster) {
 300                 Error("GetGlobalXYZ", "Cluster not found for index %d", refIndex);
 301                 return kFALSE;
 302         }
 303
 304         // determine the detector number
 305         Int_t detID = refCluster->GetDetectorIndex() + fFirstModInLayer[ilayer];
 306
 307         // get rotation matrix
 308         Double_t rot[9];
 309         fGeom->GetRotMatrix(detID, rot);
 310
 311         // get translation vector
 312         Float_t tx,ty,tz;
 313         fGeom->GetTrans(detID, tx, ty, tz);
 314
 315         // determine r and phi for the reference conversion
 316         Double_t r = -(Double_t)tx * rot[1] + (Double_t)ty * rot[0];
 317         if (ilayer == 0) r = -r;
 318         Double_t phi = TMath::ATan2(rot[1],rot[0]);
 319         if (ilayer == 0) phi -= fgkPi;
 320
 321         // sets values for X, Y, Z in global coordinates and their errors
 322         Double_t cosPhi = TMath::Cos(phi);
 323         Double_t sinPhi = TMath::Sin(phi);
 324         x =  r*cosPhi + refCluster->GetY()*sinPhi;
 325         y = -r*sinPhi + refCluster->GetY()*cosPhi;
 326         z =  refCluster->GetZ();
 327         ex2 = refCluster->GetSigmaY2()*sinPhi*sinPhi;
 328         ey2 = refCluster->GetSigmaY2()*cosPhi*cosPhi;
 329         ez2 = refCluster->GetSigmaZ2();
 330
 331         return kTRUE;
 332 }
 333
 334 //__________________________________________________________________________________
 335 AliITStrackerANN::AliITSnode* AliITStrackerANN::AddNode(Int_t refIndex)
 336
 337 /**************************************************************************
 338
 339         GENERATOR OF NEURAL NODES
 340
 341         Fills the array of neural 'nodes', which are the ITS clusters
 342         translated in the global reference frame.
 343         Given that the global coordinates are used many times, they are
 344         stored in a well-defined structure, in the form of an embedded class.
 345         Moreover, this class allows a faster navigation among points
 346         and neurons, by means of some object arrays, storing only the
 347         neurons which start from, or end to, the given node.
 348         Finally, each node contains all the other nodes which match it
 349         with respect to the fixed walues, in order to perform a faster
 350         neuron-creation phase.
 351
 352         Arguments:
 353                 1) reference index of the correlated AliITSclusterV2 object
 354
 355         Operations:
 356                 - allocates the new AliITSnode objects
 357                 - initializes its object arrays
 358                 - from the global coordinates, calculates the
 359                   'phi' and 'theta' coordinates, in order to store it
 360                   into the correct theta-slot and azimutal sector.
 361
 362         REturn values:
 363                 - the pointer of the creater AliITSnode object
 364                 - in case of errors, a waring is given and a NULL is returned
 365
 366  **************************************************************************/
 367 {
 368         // create object and set the reference
 369         AliITSnode *node = new AliITSnode;
 370         if (!node) {
 371                 Warning("AddNode", "Error occurred when allocating AliITSnode");
 372                 return 0;
 373         }
 374         node->ClusterRef() = refIndex;
 375
 376         // calls the conversion function, which makes also some checks
 377         // (layer number within range, existence of referenced cluster)
 378         if ( !GetGlobalXYZ (
 379                                 refIndex,
 380                                 node->X(), node->Y(), node->Z(),
 381                                 node->ErrX2(), node->ErrY2(), node->ErrZ2()
 382                         ) ) {return 0;}
 383
 384         // initializes the object arrays
 385         node->Matches() = new TObjArray;
 386         node->InnerOf() = new TObjArray;
 387         node->OuterOf() = new TObjArray;
 388
 389         // finds azimutal and polar sector (in degrees)
 390         Double_t phi = node->GetPhi() * 180.0 / fgkPi;
 391         Double_t theta = node->GetTheta() * 180.0 / fgkPi;
 392         Int_t isector = (Int_t)(phi / fSectorWidth);
 393         Int_t itheta = (Int_t)theta;
 394         Int_t ilayer = (refIndex & 0xf0000000) >> 28;
 395
 396         // selects the right TObjArray to store object into
 397         TObjArray *sector = (TObjArray*)fNodes[ilayer][itheta]->At(isector);
 398         sector->AddLast(node);
 399
 400         return node;
 401 }
 402
 403 //__________________________________________________________________________________
 404 void AliITStrackerANN::CreateArrayStructure(Int_t nsectors)
 405
 406 /**************************************************************************
 407
 408         ARRAY STRUCTURE CREATOR
 409
 410         Creates a structure of nested TObjArray's where the AliITSnode's
 411         have to be stored:
 412         - the first level is made by 6 arrays (one for each layer)
 413         - the second level is made by 180 arrays (one for each int part of theta)
 414         - the third level is made by a variable number of arrays
 415           (one for each azimutal sector)
 416
 417         Arguments:
 418                 1) the number of azimutal sectors
 419
 420         Operations:
 421                 - calculates the width of each sector, from the argument
 422                 - allocates and initializes all array levels
 423                 - sets a flag which tells the user if this NECESSARY operation
 424                   has been performed (it is needed BEFORE performing tracking)
 425
 426  **************************************************************************/
 427 {
 428         // Set the number of sectors and their width.
 429         fSectorNum   = nsectors;
 430         fSectorWidth = 360.0 / (Double_t)fSectorNum;
 431         if (fMsgLevel >= 2) {
 432                 cout << fSectorNum << " sectors --> sector width (degrees) = " << fSectorWidth << endl;
 433         }
 434
 435         // Meaningful indexes
 436         Int_t ilayer, isector, itheta;
 437
 438         // Mark for the created objects
 439         TObjArray *sector = 0;
 440
 441         // First index: layer
 442         fNodes = new TObjArray**[fNLayers];
 443         for (ilayer = 0; ilayer < fNLayers; ilayer++) {
 444                 fNodes[ilayer] = new TObjArray*[180];
 445                 for (itheta = 0; itheta < 180; itheta++) fNodes[ilayer][itheta] = 0;
 446                 for (itheta = 0; itheta < 180; itheta++) {
 447                         fNodes[ilayer][itheta] = new TObjArray(nsectors);
 448                         for (isector = 0; isector < nsectors; isector++) {
 449                                 sector = new TObjArray;
 450                                 sector->SetOwner();
 451                                 fNodes[ilayer][itheta]->AddAt(sector, isector);
 452                         }
 453                 }
 454         }
 455
 456         // Sets a checking flag to TRUE.
 457         // This flag is checked before filling up the arrays with the points.
 458         fStructureOK = kTRUE;
 459 }
 460
 461 //__________________________________________________________________________________
 462 Int_t AliITStrackerANN::ArrangePoints(char *exportFile)
 463
 464 /**************************************************************************
 465
 466         POINTS LOCATOR
 467
 468         This function assembles the operation from the other above methods,
 469         and fills the arrays with the clusters already stored in the
 470         layers of the tracker.
 471         Then, in order to use this method, the user MUSTs call LoadClusters()
 472         before.
 473
 474         Arguments:
 475                 1) string for a file name where the global ccordinates
 476                    of all points can be exported (optional).
 477                         If this file must not be created, simply pass a NULL argument
 478
 479         Operations:
 480                 - for each AliITSclusterV2 in each AliITSlayer, a ne AliITSnode
 481                   is created and stored in the correct location.
 482
 483         Return values:
 484                 - the number of stored points
 485                 - when errors occur, or no points are found, 0 is returned
 486
 487  **************************************************************************/
 488 {
 489         // Check if the array structure has been created
 490         if (!fStructureOK) {
 491                 Error("ArrangePoints", "Structure NOT defined. Call CreateArrayStructure() first");
 492                 return 0;
 493         }
 494
 495         // meaningful indexes
 496         Int_t ientry, ilayer, nentries = 0, index;
 497         Int_t nPtsLayer = 0;
 498
 499         // if the argument is not NULL, a file is opened
 500         fstream file(exportFile, ios::out);
 501         if (!exportFile || file.fail()) {
 502                 file.close();
 503                 exportFile = 0;
 504         }
 505
 506         // scan all layers for node creation
 507         for (ilayer = 0; ilayer < fNLayers; ilayer++) {
 508                 nPtsLayer = GetNumberOfClustersLayer(ilayer);
 509                 if (fMsgLevel >= 1) {
 510                         cout << "Layer " << ilayer << " --> " << nPtsLayer << " clusters" << endl;
 511                 }
 512                 for (ientry = 0; ientry < nPtsLayer; ientry++) {
 513                         // calculation of cluster index : (Bit mask LLLLIIIIIIIIIIII)
 514                         // (L = bits used for layer)
 515                         // (I = bits used for position in layer)
 516                         index = ilayer << 28;
 517                         index += ientry;
 518                         // add new AliITSnode object
 519                         AliITSnode *n = AddNode(index);
 520                         if ( (n != NULL) && exportFile ) {
 521                                 file << index << ' ' << n->X() << ' ' << n->Y() << ' ' << n->Z() << endl;
 522                         }
 523                 }
 524                 nentries += nPtsLayer;
 525         }
 526
 527         // a conventional final message is put at the end of file
 528         if (exportFile) {
 529                 file << "-1 0.0 0.0 0.0" << endl;
 530                 file.close();
 531         }
 532
 533         // returns the number of points processed
 534         return nentries;
 535 }
 536
 537 //__________________________________________________________________________________
 538 void AliITStrackerANN::StoreOverallMatches()
 539
 540 /**************************************************************************
 541
 542         NODE-MATCH ANALYSIS
 543
 544         Once the nodes have been created, a firs analysis is to check
 545         what pairs will satisfy at least the 'fixed' cuts (theta, helix-match)
 546         and the most permissive (= larger) curvature cut.
 547         All these node pairs are suitable for neuron creation.
 548         In fact, when performing a Neural Tracking step, the only further check
 549         will be a check against the current curvature step, while the other
 550         are always the same.
 551         For thi purpose, each AliITSnode has a data member, named 'fMatches'
 552         which contains references to all other AliITSnodes in the successive layer
 553         that form, with it, a 'good' pair, with respect to the above cited cuts.
 554         Then, in each step for neuron creation, the possible neurons starting from
 555         each node will be searched ONLY within the nodes referenced in fMatches.
 556         This, of course, speeds up a lot the neuron creation procedure, at the
 557         cost of some memory occupation, which results not to be critical.
 558
 559         Operations:
 560                 - for each AliITSnode, matches are found according to the criteria
 561                   expressed above, and stored in the node->fMatches array
 562
 563  **************************************************************************/
 564 {
 565         // meaningful counters
 566         Int_t ilayer, isector, itheta1, itheta2, check;
 567         TObjArray *list1 = 0, *list2 = 0;
 568         AliITSnode *node1 = 0, *node2 = 0;
 569         Double_t thetaMin, thetaMax;
 570         Int_t imin, imax;
 571
 572         // Scan for each sector
 573         for (isector = 0; isector < fSectorNum; isector++) {
 574                 // sector is chosen once for both lists
 575                 for (ilayer = 0; ilayer < fNLayers - 1; ilayer++) {
 576                         for (itheta1 = 0; itheta1 < 180; itheta1++) {
 577                                 list1 = (TObjArray*)fNodes[ilayer][itheta1]->At(isector);
 578                                 TObjArrayIter iter1(list1);
 579                                 while ( (node1 = (AliITSnode*)iter1.Next()) ) {
 580                                         if (node1->IsUsed()) continue;
 581                                         // clear an eventually already present array
 582                                         // node1->Matches()->Clear();
 583                                         // get the global coordinates and defines the theta interval from cut
 584                                         thetaMin = (node1->GetTheta() * 180.0 / fgkPi) - fPolarInterval;
 585                                         thetaMax = (node1->GetTheta() * 180.0 / fgkPi) + fPolarInterval;
 586                                         imin = (Int_t)thetaMin;
 587                                         imax = (Int_t)thetaMax;
 588                                         if (imin < 0) imin = 0;
 589                                         if (imax > 179) imax = 179;
 590                                         // loop on the selected theta slots
 591                                         for (itheta2 = imin; itheta2 <= imax; itheta2++) {
 592                                                 list2 = (TObjArray*)fNodes[ilayer + 1][itheta2]->At(isector);
 593                                                 TObjArrayIter iter2(list2);
 594                                                 while ( (node2 = (AliITSnode*)iter2.Next()) ) {
 595                                                         check = PassAllCuts(node1, node2, fCurvNum - 1, fVertexX, fVertexY, fVertexZ);
 596                                                         if (check == 0) {
 597                                                                 node1->Matches()->AddLast(node2);
 598                                                         }
 599                                                 } // while (node2...)
 600                                         } // for (itheta2...)
 601                                 } // while (node1...)
 602                         } // for (itheta...)
 603                 } // for (ilayer...)
 604         } // for (isector...)
 605 }
 606
 607 //__________________________________________________________________________________
 608 Int_t AliITStrackerANN::PassAllCuts
 609 (AliITSnode *inner, AliITSnode *outer, Int_t curvStep,
 610  Double_t vx, Double_t vy, Double_t vz)
 611 {
 612 // ***********************************************************************************
 613 //
 614 //      This check is called in the above method for finding the matches of each node
 615 //      It check the passed point pair against all the fixed cuts and a specified
 616 //      curvature cut, among all the ones which have been defined.
 617 //      The cuts need a vertex-constraint, which is not absolute, but it is passed
 618 //      as argument.
 619 //
 620 //      Arguments:
 621 //              1) the point in the inner layer
 622 //              2) the point in the outer layer
 623 //              3) curvature step for the curvature cut check (preferably the last)
 624 //              4) X of the used vertex
 625 //              5) Y of the used vertex
 626 //              6) Z of the used vertex
 627 //
 628 //      Operations:
 629 //              - if necessary, swaps the two points
 630 //                (the first must be in the innermost of the two layers)
 631 //              - checks for the theta cuts
 632 //              - calculates the circle passing through the vertex
 633 //                and the given points and checks for the curvature cut
 634 //              - using the radius calculated there, checks for the helix-math cut
 635 //
 636 //      Return values:
 637 //              0 - All cuts passed
 638 //              1 - theta 2D cut not passed
 639 //              2 - theta 3D cut not passed
 640 //              3 - curvature calculated but cut not passed
 641 //              4 - curvature not calculated (division by zero)
 642 //              5 - helix cut not passed
 643 //              6 - curvature inxed out of range
 644 //
 645 // ***********************************************************************************
 646
 647         // Check for curvature index
 648         if (curvStep < 0 || curvStep >= fCurvNum) return 6;
 649
 650         // Swap points in order that r1 < r2
 651         AliITSnode *temp = 0;
 652         if (outer->GetLayer() < inner->GetLayer()) {
 653                 temp = outer;
 654                 outer = inner;
 655                 inner = temp;
 656         }
 657
 658         // All cuts are variable according to the layer of the
 659         // innermost point (the other point will surely be
 660         // in the further one, because we don't check poin pairs
 661         // which are not in adjacent layers)
 662         // The reference is given by the innermost point.
 663         Int_t layRef = inner->GetLayer();
 664
 665         // The calculations in the transverse plane are made in
 666         // a shifted reference frame, whose origin corresponds to
 667         // the reference point passed in the argument.
 668         Double_t xIn = inner->X() - vx;
 669         Double_t xOut = outer->X() - vx;
 670         Double_t yIn = inner->Y() - vy;
 671         Double_t yOut = outer->Y() - vy;
 672         Double_t zIn = inner->Z() - vz;
 673         Double_t zOut = outer->Z() - vz;
 674         Double_t rIn = TMath::Sqrt(xIn*xIn + yIn*yIn);
 675         Double_t rOut = TMath::Sqrt(xOut*xOut + yOut*yOut);
 676
 677         // Check for theta cut.
 678         // There are two different angular cuts:
 679         // one w.r. to the angle in the 2-dimensional r-z plane...
 680         Double_t dthetaRZ;
 681         TVector3 origin2innerRZ(zIn, rIn, 0.0);
 682         TVector3 inner2outerRZ(zOut - zIn, rOut - rIn, 0.0);
 683         dthetaRZ = origin2innerRZ.Angle(inner2outerRZ) * 180.0 / fgkPi;
 684         if (dthetaRZ > fThetaCut2D[layRef]) {
 685                 return 1;
 686                 // r-z theta cut not passed ---> 1
 687         }
 688         // ...and another w.r. to the angle in the 3-dimensional x-y-z space
 689         Double_t dthetaXYZ;
 690         TVector3 origin2innerXYZ(xIn, yIn, zIn);
 691         TVector3 inner2outerXYZ(xOut - xIn, yOut - yIn, zOut - zIn);
 692         dthetaXYZ = origin2innerXYZ.Angle(inner2outerXYZ) * 180.0 / fgkPi;
 693         if (dthetaXYZ > fThetaCut3D[layRef]) {
 694                 return 2;
 695                 // x-y-z theta cut not passed ---> 2
 696         }
 697
 698         // Calculation & check of curvature
 699         Double_t dx = xIn - xOut;
 700         Double_t dy = yIn - yOut;
 701         Double_t num = 2.0 * (xIn*yOut - xOut*yIn);
 702         Double_t den = rIn*rOut*sqrt(dx*dx + dy*dy);
 703         Double_t curv = 0.;
 704         if (den != 0.) {
 705                 curv = TMath::Abs(num / den);
 706                 if (curv > fCurvCut[curvStep]) {
 707                         return 3;
 708                         // curvature too large for cut ---> 3
 709                 }
 710         }
 711         else {
 712                 Error("PassAllCuts", "Curvature calculations gives zero denominator");
 713                 return 4;
 714                 // error: denominator = 0 ---> 4
 715         }
 716
 717         // Calculation & check of helix matching
 718         Double_t helMatch = 0.0;
 719         Double_t arcIn = 2.0 * rIn * curv;
 720         Double_t arcOut = 2.0 * rOut * curv;
 721         if (arcIn > -1.0 && arcIn < 1.0)
 722                 arcIn = TMath::ASin(arcIn);
 723         else
 724                 arcIn = ((arcIn > 0.0) ? 0.5 : 1.5) * TMath::Pi();
 725         if (arcOut > -1.0 && arcOut < 1.0)
 726                 arcOut = TMath::ASin(arcOut);
 727         else
 728                 arcOut = ((arcOut > 0.0) ? 0.5 : 1.5) * TMath::Pi();
 729         arcIn /= 2.0 * curv;
 730         arcOut /= 2.0 * curv;
 731         if (arcIn == 0.0 || arcOut == 0.0) {
 732                 Error("PassAllCuts", "Calculation returns zero-length arcs: l1=%f, l2=%f", arcIn, arcOut);
 733                 return 4;
 734                 // error: circumference arcs seem to equal zero ---> 4
 735         }
 736         helMatch = TMath::Abs(zIn / arcIn - zOut / arcOut);
 737         if (helMatch > fHelixMatchCut[layRef]) {
 738                 return 5;
 739                 // helix match cut not passed ---> 5
 740         }
 741
 742         // ALL CUTS PASSED ---> 0
 743         return 0;
 744 }
 745
 746 //__________________________________________________________________________________
 747 Bool_t AliITStrackerANN::PassCurvCut
 748 (AliITSnode *inner, AliITSnode *outer,
 749  Int_t curvStep,
 750  Double_t vx, Double_t vy, Double_t vz)
 751 {
 752 //***********************************************************************************
 753 //
 754 //      This method operates essentially like the above one, but it is used
 755 //      during a single step of Neural Tracking, where the curvature cut
 756 //      changes.
 757 //      Then, not necessaryly all the nodes stored in the fMatches array
 758 //      will be suitable for neuron creation in an intermediate step.
 759 //
 760 //      It has the same arguments of the PassAllCuts() method, but
 761 //  the theta cut is not checked.
 762 //  Moreover, it has a boolean return value.
 763 //
 764 //***********************************************************************************
 765
 766         // Check for curvature index
 767         if (curvStep < 0 || curvStep >= fCurvNum) return 6;
 768
 769         // Find the reference layer
 770         Int_t layIn = inner->GetLayer();
 771         Int_t layOut = outer->GetLayer();
 772         Int_t layRef = (layIn < layOut) ? layIn : layOut;
 773
 774         // The calculations in the transverse plane are made in
 775         // a shifted reference frame, whose origin corresponds to
 776         // the reference point passed in the argument.
 777         Double_t xIn = inner->X() - vx;
 778         Double_t xOut = outer->X() - vx;
 779         Double_t yIn = inner->Y() - vy;
 780         Double_t yOut = outer->Y() - vy;
 781         Double_t zIn = inner->Z() - vz;
 782         Double_t zOut = outer->Z() - vz;
 783         Double_t rIn = TMath::Sqrt(xIn*xIn + yIn*yIn);
 784         Double_t rOut = TMath::Sqrt(xOut*xOut + yOut*yOut);
 785
 786         // Calculation & check of curvature
 787         Double_t dx = xIn - xOut;
 788         Double_t dy = yIn - yOut;
 789         Double_t num = 2.0 * (xIn*yOut - xOut*yIn);
 790         Double_t den = rIn*rOut*sqrt(dx*dx + dy*dy);
 791         Double_t curv = 0.;
 792         /* OLD VERSION
 793         if (den != 0.) {
 794                 curv = TMath::Abs(num / den);
 795                 if (curv > fCurvCut[curvStep]) return kFALSE;
 796                 return kTRUE;
 797         }
 798         else
 799                 return kFALSE;
 800         */
 801         // NEW VERSION
 802         if (den != 0.) {
 803                 curv = TMath::Abs(num / den);
 804                 if (curv > fCurvCut[curvStep]) {
 805                         return kFALSE;
 806                 }
 807         }
 808         else {
 809                 Error("PassAllCuts", "Curvature calculations gives zero denominator");
 810                 return kFALSE;
 811         }
 812
 813         // Calculation & check of helix matching
 814         Double_t helMatch = 0.0;
 815         Double_t arcIn = 2.0 * rIn * curv;
 816         Double_t arcOut = 2.0 * rOut * curv;
 817         if (arcIn > -1.0 && arcIn < 1.0)
 818                 arcIn = TMath::ASin(arcIn);
 819         else
 820                 arcIn = ((arcIn > 0.0) ? 0.5 : 1.5) * TMath::Pi();
 821         if (arcOut > -1.0 && arcOut < 1.0)
 822                 arcOut = TMath::ASin(arcOut);
 823         else
 824                 arcOut = ((arcOut > 0.0) ? 0.5 : 1.5) * TMath::Pi();
 825         arcIn /= 2.0 * curv;
 826         arcOut /= 2.0 * curv;
 827         if (arcIn == 0.0 || arcOut == 0.0) {
 828                 Error("PassAllCuts", "Calculation returns zero-length arcs: l1=%f, l2=%f", arcIn, arcOut);
 829                 return 4;
 830                 // error: circumference arcs seem to equal zero ---> 4
 831         }
 832         helMatch = TMath::Abs(zIn / arcIn - zOut / arcOut);
 833         return (helMatch <= fHelixMatchCut[layRef]);
 834 }
 835
 836 //__________________________________________________________________________________
 837 void AliITStrackerANN::PrintMatches(Bool_t stop)
 838 {
 839 //      Prints the list of points which appear to match
 840 //      each one of them, according to the preliminary
 841 //      overall cuts.
 842 //      The arguments states if a pause is required after printing
 843 //      the matches for each one. In this case, a keypress is needed.
 844
 845         TObjArray *sector = 0;
 846         Int_t ilayer, isector, itheta, nF;
 847         AliITSnode *node1 = 0, *node2 = 0;
 848         //AliITSclusterV2 *cluster1 = 0, *cluster2 = 0;
 849
 850         for (ilayer = 0; ilayer < 6; ilayer++) {
 851                 for (isector = 0; isector < fSectorNum; isector++) {
 852                         for (itheta = 0; itheta < 180; itheta++) {
 853                                 sector = (TObjArray*)fNodes[ilayer][itheta]->At(isector);
 854                                 TObjArrayIter points(sector);
 855                                 while ( (node1 = (AliITSnode*)points.Next()) ) {
 856                                         nF = (Int_t)node1->Matches()->GetEntries();
 857                                         cout << "Node layer: " << node1->GetLayer() << " --> ";
 858                                         if (!nF) {
 859                                                 cout << "NO Matches!!!" << endl;
 860                                                 continue;
 861                                         }
 862                                         cout << nF << " Matches" << endl;
 863                                         cout << "Reference cluster: " << hex << node1->ClusterRef() << endl;
 864                                         TObjArrayIter matches(node1->Matches());
 865                                         while ( (node2 = (AliITSnode*)matches.Next()) ) {
 866                                                 cout << "Match with " << hex << node2->ClusterRef() << endl;
 867                                         }
 868                                         if (stop) {
 869                                                 cout << "Press a key" << endl;
 870                                                 cin.get();
 871                                         }
 872                                 }
 873                         }
 874                 }
 875         }
 876 }
 877
 878 //__________________________________________________________________________________
 879 void AliITStrackerANN::ResetNodes(Int_t isector)
 880 {
 881 /***********************************************************************************
 882
 883         NODE NEURON ARRAY CLEANER
 884
 885         After a neural tracking step, this method
 886         clears the arrays 'fInnerOf' and 'fOuterOf' of each AliITSnode
 887
 888         Arguments:
 889                 - the sector where the operation is being executed
 890
 891  ***********************************************************************************/
 892
 893         Int_t ilayer, itheta;
 894         TObjArray *sector = 0;
 895         AliITSnode *node = 0;
 896         for (ilayer = 0; ilayer < fNLayers; ilayer++) {
 897                 for (itheta = 0; itheta < 180; itheta++) {
 898                         sector = (TObjArray*)fNodes[ilayer][itheta]->At(isector);
 899                         TObjArrayIter iter(sector);
 900                         for (;;) {
 901                                 node = (AliITSnode*)iter.Next();
 902                                 if (!node) break;
 903                                 node->InnerOf()->Clear();
 904                                 node->OuterOf()->Clear();
 905                                 /*
 906                                 delete node->InnerOf();
 907                                 delete node->OuterOf();
 908                                 node->InnerOf() = new TObjArray;
 909                                 node->OuterOf() = new TObjArray;
 910                                 */
 911                         }
 912                 }
 913         }
 914 }
 915
 916 //__________________________________________________________________________________
 917 Int_t AliITStrackerANN::CreateNeurons
 918 (AliITSnode *node, Int_t curvStep, Double_t vx, Double_t vy, Double_t vz)
 919 {
 920         // This method is used to create alle suitable neurons starting from
 921         // a single AliITSnode. Each unit is also stored in the fInnerOf array
 922         // of the passed node, and in the fOuterOf array of the other neuron edge.
 923         // In the new implementation of the intermediate check steps, a further one
 924         // is made, which chechs how well a helix is matched by three points
 925         // in three consecutive layers.
 926         // Then, a vertex-constrained check is made with vertex located
 927         // in a layer L, for points in layers L+1 and L+2.
 928         //
 929         // In order to do this, the creator works recursively, in a tree-visit like operation.
 930         // The neurons are effectively created only if the node argument passed is in
 931         // the 5th layer (they are created between point of 5th and 6th layer).
 932         // If the node is in an inner layer, its coordinates are passet as vertex for a nested
 933         // call of the same function in the next two layers.
 934         //
 935         // Arguments:
 936         //      1) reference node
 937         //      2) current curvature cut step
 938         //      3) X of referenced temporary (not primary) vertex
 939         //      4) Y of referenced temporary (not primary) vertex
 940         //      5) Z of referenced temporary (not primary) vertex
 941         //
 942         // Operations:
 943         //      - if the layer is the 5th, neurons are created with nodes
 944         //        in the fMatches array of the passed node
 945         //      - otherwise, the X, Y, Z of the passed node are given as
 946         //        vertex and the same procedure is recursively called for all
 947         //        nodes in the fMatches array of the passed one.
 948         //
 949         // Return values:
 950         //      - the total number of neurons created from the passed one
 951         //        summed with all neurons created from all nodes well matched with it
 952         //        (assumes a meaning only for nodes in the first layer)
 953
 954         // local variables
 955         Int_t made = 0;         // counter
 956         Bool_t found = 0;       // flag
 957         AliITSnode *match = 0;  // cursor for a AliITSnode array
 958         AliITSneuron *unit = 0; // cursor for a AliITSneuron array
 959
 960         // --> Case 0: the passed node has already been used
 961         // as member of a track found in a previous step.
 962         // In this case, of course, the function exits.
 963         if (node->IsUsed()) return 0;
 964
 965         // --> Case 1: there are NO well-matched points.
 966         // This can happen in all ITS layers, but it happens **for sure**
 967         // for a node in the 'last' layer.
 968         // Even in this case, the function exits.
 969         if (node->Matches()->IsEmpty()) return 0;
 970
 971         // --> Case 2: there are well-matched points.
 972         // In this case, the function creates a neuron for each
 973         // well-matched pair (according to the cuts for the current step)
 974         // Moreover, before storing the neuron, a check is necessary
 975         // to avoid the duplicate creation of the same neuron twice.
 976         // (This could happen if the 3 last arguments of the function
 977         // are close enough to cause a good match for the current step
 978         // between two points, independently of their difference).
 979         // Finally, a node is skipped if it has already been used.
 980         // For each matched point for which a neuron is created, the procedure is
 981         // recursively called.
 982         TObjArrayIter matches(node->Matches());
 983         while ( (match = (AliITSnode*)matches.Next()) ) {
 984                 if (match->IsUsed()) continue;
 985                 if (!PassCurvCut(node, match, curvStep, vx, vy, vz)) continue;
 986                 found = kFALSE;
 987                 if (!node->InnerOf()->IsEmpty()) {
 988                         TObjArrayIter presentUnits(node->InnerOf());
 989                         while ( (unit = (AliITSneuron*)presentUnits.Next()) ) {
 990                                 if (unit->Inner() == node && unit->Outer() == match) {
 991                                         found = kTRUE;
 992                                         break;
 993                                 }
 994                         }
 995                 }
 996                 if (found) continue;
 997                 AliITSneuron *unit = new AliITSneuron(node, match, fEdge2, fEdge1);
 998                 fNeurons->AddLast(unit);
 999                 node->InnerOf()->AddLast(unit);
1000                 match->OuterOf()->AddLast(unit);
1001                 made += CreateNeurons(match, curvStep, node->X(), node->Y(), node->Z());
1002                 made++;
1003         }
1004
1005         // Of course, the return value contains the number of neurons
1006         // counting in also the oned created in all levels of recursive calls.
1007         return made;
1008 }
1009
1010 //__________________________________________________________________________________
1011 Int_t AliITStrackerANN::CreateNetwork(Int_t sector, Int_t curvStep)
1012 {
1013         // This function simply recalls the CreateNeurons() method for each node
1014         // in the first layer, for the current sector.
1015         // This generates the whole network, thanks to the recursive calls.
1016         //
1017         // Arguments:
1018         //       1) current sector
1019         //       2) current curvature step
1020         //
1021         // Operations:
1022         //      - scans the nodes array for all theta's in the current sector
1023         //        and layer 0, and calls the CreateNeurons() function.
1024
1025         // removes all eventually present neurons
1026         if (fNeurons) delete fNeurons;
1027         fNeurons = new TObjArray;
1028         fNeurons->SetOwner(kTRUE);
1029
1030         // calls the ResetNodes() function to free the AliITSnode arrays
1031         if (fMsgLevel >= 2) {
1032                 cout << "Sector " << sector << " PHI = ";
1033                 cout << fSectorWidth * (Double_t)sector << " --> ";
1034                 cout << fSectorWidth * (Double_t)(sector + 1) << endl;
1035                 cout << "Curvature step " << curvStep << " [cut = " << fCurvCut[curvStep] << "]" << endl;
1036         }
1037         ResetNodes(sector);
1038
1039         // meaningful counters
1040         Int_t itheta, neurons = 0;
1041         TObjArray *lstSector = 0;
1042
1043         // NEW VERSION
1044         Double_t vx[6], vy[6], vz[6];
1045         AliITSnode *p[6] = {0, 0, 0, 0, 0, 0};
1046         for (itheta = 0; itheta < 180; itheta++) {
1047                 lstSector = (TObjArray*)fNodes[0][itheta]->At(sector);
1048                 TObjArrayIter lay0(lstSector);
1049                 while ( (p[0] = (AliITSnode*)lay0.Next()) ) {
1050                         if (p[0]->IsUsed()) continue;
1051                         vx[0] = fVertexX;
1052                         vy[0] = fVertexY;
1053                         vz[0] = fVertexZ;
1054                         neurons += CreateNeurons(p[0], curvStep, fVertexX, fVertexY, fVertexZ);
1055                         /*
1056                         TObjArrayIter lay1(p[0]->Matches());
1057                         while ( (p[1] = (AliITSnode*)lay1.Next()) ) {
1058                                 if (p[1]->IsUsed()) continue;
1059                                 if (!PassCurvCut(p[0], p[1], curvStep, vx[0], vy[0], vz[0])) continue;
1060                                 unit = new AliITSneuron;
1061                                 unit->Inner() = p[0];
1062                                 unit->Outer() = p[1];
1063                                 unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1064                                 unit->fGain = new TObjArray;
1065                                 fNeurons->AddLast(unit);
1066                                 p[0]->InnerOf()->AddLast(unit);
1067                                 p[1]->OuterOf()->AddLast(unit);
1068                                 neurons++;
1069                                 vx[1] = p[0]->X();
1070                                 vy[1] = p[0]->Y();
1071                                 vz[1] = p[0]->Z();
1072                                 TObjArrayIter lay2(p[1]->Matches());
1073                                 while ( (p[2] = (AliITSnode*)lay2.Next()) ) {
1074                                         if (p[2]->IsUsed()) continue;
1075                                         if (!PassCurvCut(p[1], p[2], curvStep, vx[1], vy[1], vz[1])) continue;
1076                                         unit = new AliITSneuron;
1077                                         unit->Inner() = p[1];
1078                                         unit->Outer() = p[2];
1079                                         unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1080                                         unit->fGain = new TObjArray;
1081                                         fNeurons->AddLast(unit);
1082                                         p[1]->InnerOf()->AddLast(unit);
1083                                         p[2]->OuterOf()->AddLast(unit);
1084                                         neurons++;
1085                                         vx[2] = p[1]->X();
1086                                         vy[2] = p[1]->Y();
1087                                         vz[2] = p[1]->Z();
1088                                         TObjArrayIter lay3(p[2]->Matches());
1089                                         while ( (p[3] = (AliITSnode*)lay3.Next()) ) {
1090                                                 if (p[3]->IsUsed()) continue;
1091                                                 if (!PassCurvCut(p[2], p[3], curvStep, vx[2], vy[2], vz[2])) continue;
1092                                                 unit = new AliITSneuron;
1093                                                 unit->Inner() = p[2];
1094                                                 unit->Outer() = p[3];
1095                                                 unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1096                                                 unit->fGain = new TObjArray;
1097                                                 fNeurons->AddLast(unit);
1098                                                 p[2]->InnerOf()->AddLast(unit);
1099                                                 p[3]->OuterOf()->AddLast(unit);
1100                                                 neurons++;
1101                                                 vx[3] = p[2]->X();
1102                                                 vy[3] = p[2]->Y();
1103                                                 vz[3] = p[2]->Z();
1104                                                 TObjArrayIter lay4(p[3]->Matches());
1105                                                 while ( (p[4] = (AliITSnode*)lay4.Next()) ) {
1106                                                         if (p[4]->IsUsed()) continue;
1107                                                         if (!PassCurvCut(p[3], p[4], curvStep, vx[3], vy[3], vz[3])) continue;
1108                                                         unit = new AliITSneuron;
1109                                                         unit->Inner() = p[3];
1110                                                         unit->Outer() = p[4];
1111                                                         unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1112                                                         unit->fGain = new TObjArray;
1113                                                         fNeurons->AddLast(unit);
1114                                                         p[3]->InnerOf()->AddLast(unit);
1115                                                         p[4]->OuterOf()->AddLast(unit);
1116                                                         neurons++;
1117                                                         vx[4] = p[3]->X();
1118                                                         vy[4] = p[3]->Y();
1119                                                         vz[4] = p[3]->Z();
1120                                                         TObjArrayIter lay5(p[4]->Matches());
1121                                                         while ( (p[5] = (AliITSnode*)lay5.Next()) ) {
1122                                                                 if (p[5]->IsUsed()) continue;
1123                                                                 if (!PassCurvCut(p[4], p[5], curvStep, vx[4], vy[4], vz[4])) continue;
1124                                                                 unit = new AliITSneuron;
1125                                                                 unit->Inner() = p[4];
1126                                                                 unit->Outer() = p[5];
1127                                                                 unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1128                                                                 unit->fGain = new TObjArray;
1129                                                                 fNeurons->AddLast(unit);
1130                                                                 p[4]->InnerOf()->AddLast(unit);
1131                                                                 p[5]->OuterOf()->AddLast(unit);
1132                                                                 neurons++;
1133                                                         } // while (p[5])
1134                                                 } // while (p[4])
1135                                         } // while (p[3])
1136                                 } // while (p[2])
1137                         } // while (p[1])
1138                         */
1139                 } // while (p[0])
1140         } // for (itheta...)
1141         // END OF NEW VERSION
1142
1143         /* OLD VERSION
1144         for (ilayer = 0; ilayer < 6; ilayer++) {
1145                 for (itheta = 0; itheta < 180; itheta++) {
1146                         lstSector = (TObjArray*)fNodes[ilayer][itheta]->At(sector_idx);
1147                         TObjArrayIter inners(lstSector);
1148                         while ( (inner = (AliITSnode*)inners.Next()) ) {
1149                                 if (inner->GetUser() >= 0) continue;
1150                                 TObjArrayIter outers(inner->Matches());
1151                                 while ( (outer = (AliITSnode*)outers.Next()) ) {
1152                                         if (outer->GetUser() >= 0) continue;
1153                                         if (!PassCurvCut(inner, outer, curvStep, fVX, fVY, fVZ)) continue;
1154                                         unit = new AliITSneuron;
1155                                         unit->Inner() = inner;
1156                                         unit->Outer() = outer;
1157                                         unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1158                                         unit->fGain = new TObjArray;
1159                                         fNeurons->AddLast(unit);
1160                                         inner->InnerOf()->AddLast(unit);
1161                                         outer->OuterOf()->AddLast(unit);
1162                                         neurons++;
1163                                 } // for (;;)
1164                         } // for (;;)
1165                 } // for (itheta...)
1166         } // for (ilayer...)
1167         */
1168
1169         fNeurons->SetOwner();
1170         return neurons;
1171  }
1172
1173  //__________________________________________________________________________________
1174 Int_t AliITStrackerANN::LinkNeurons()
1175 {
1176 /***********************************************************************************
1177
1178         SYNAPSIS GENERATOR
1179
1180         Scans the whole neuron array, in order to find all neuron pairs
1181         which are connected in sequence and share a positive weight.
1182         For each of them, an AliITSlink is created, which stores
1183         the weight value, and will allow for a faster calculation
1184         of the total neural input for each updating cycle.
1185
1186         Every neuron contains an object array which stores all AliITSlink
1187         objects which point to sequenced units, with the respective weights.
1188
1189         Return value:
1190                 - the number of link created for the neural network.
1191                   If they are 0, no updating can be done and the step is skipped.
1192
1193  ***********************************************************************************/
1194
1195         // meaningful indexes
1196         Int_t total = 0;
1197         Double_t weight = 0.0;
1198         TObjArrayIter neurons(fNeurons), *iter;
1199         AliITSneuron *neuron = 0, *test = 0;
1200
1201         // scan in the neuron array
1202         for (;;) {
1203                 neuron = (AliITSneuron*)neurons.Next();
1204                 if (!neuron) break;
1205                 // checks for neurons startin from its head ( -> )
1206                 iter = (TObjArrayIter*)neuron->Inner()->OuterOf()->MakeIterator();
1207                 for (;;) {
1208                         test = (AliITSneuron*)iter->Next();
1209                         if (!test) break;
1210                         weight = Weight(test, neuron);
1211                         if (weight > 0.0) neuron->Gain()->AddLast(new AliITSlink(weight, test));
1212                         total++;
1213                 }
1214                 delete iter;
1215                 // checks for neurons ending in its tail ( >- )
1216                 iter = (TObjArrayIter*)neuron->Outer()->InnerOf()->MakeIterator();
1217                 for (;;) {
1218                         test = (AliITSneuron*)iter->Next();
1219                         if (!test) break;
1220                         weight = Weight(neuron, test);
1221                         if (weight > 0.0) neuron->Gain()->AddLast(new AliITSlink(weight, test));
1222                         total++;
1223                 }
1224                 delete iter;
1225         }
1226         return total;
1227 }
1228
1229 //__________________________________________________________________________________
1230 Bool_t AliITStrackerANN::Update()
1231 {
1232 /***********************************************************************************
1233
1234         Performs a single updating cycle.
1235
1236         Operations:
1237                 - for each neuron, gets the activation with the neuron Activate() method
1238                 - checks if stability has been reached (compare mean activation variation
1239                   with the stability threshold data member)
1240
1241         Return values:
1242                 - kTRUE means that the neural network has stabilized
1243                 - kFALSE means that another updating cycle is needed
1244
1245  ***********************************************************************************/
1246
1247         Double_t actVar = 0.0, totDiff = 0.0;
1248         TObjArrayIter iter(fNeurons);
1249         AliITSneuron *unit;
1250         for (;;) {
1251                 unit = (AliITSneuron*)iter.Next();
1252                 if (!unit) break;
1253                 actVar = unit->Activate(fTemperature);
1254                 // calculation the relative activation variation
1255                 totDiff += actVar;
1256         }
1257         totDiff /= fNeurons->GetSize();
1258         return (totDiff < fStabThreshold);
1259 }
1260
1261 //__________________________________________________________________________________
1262 void AliITStrackerANN::FollowChains(Int_t sector)
1263 {
1264 /***********************************************************************************
1265
1266         CHAINS CREATION
1267
1268         After that the neural network has stabilized,
1269         the final step is to create polygonal chains
1270         of clusters, one in each layer, which represent
1271         the tracks recognized by the neural algorithm.
1272         This is made by means of a choice of the best
1273         neuron among the ones starting from each point.
1274
1275         Once that such neuron is selected, its inner point
1276         will set the 'fNext' field to its outer point, and
1277         similarly, its outer point will set the 'fPrev' field
1278         to its inner point.
1279         This defines a bi-directional sequence.
1280
1281         In this procedure, it can happen that many neurons
1282         which have the head of the arrow in a given node, will
1283         all select as best following the neuron with the largest
1284         activation starting in that point.
1285         This results in MANY nodes which have the same 'fNext'.
1286         But, this field will be set to NULL for all these points,
1287         but the only one which is pointed by the 'fPrev' field
1288         of this shared node.
1289
1290  ***********************************************************************************/
1291
1292         // meaningful counters
1293         Int_t itheta, ilayer;
1294         TObjArray *lstSector = 0;
1295         Double_t test = fActMinimum;
1296         AliITSnode *p = 0;
1297         AliITSneuron *n = 0;
1298
1299         // scan the whole collection of nodes
1300         for (ilayer = 0; ilayer < fNLayers; ilayer++) {
1301                 for (itheta = 0; itheta < 180; itheta++) {
1302                         // get the array of point in a given layer/theta-slot/sector
1303                         lstSector = (TObjArray*)fNodes[ilayer][itheta]->At(sector);
1304                         TObjArrayIter nodes(lstSector);
1305                         while ( (p = (AliITSnode*)nodes.Next()) ) {
1306                                 // if the point is used, it is skipped
1307                                 if (p->IsUsed()) continue;
1308                                 // initially, fNext points to nothing, and
1309                                 // the comparison value is set to the minimum activation
1310                                 // which allows to say that a neuron is turned 'on'
1311                                 // a node from which only 'off' neurons start is probably
1312                                 // a noise point, which will be excluded from the reconstruction.
1313                                 test = fActMinimum;
1314                                 p->Next() = 0;
1315                                 TObjArrayIter innerof(p->InnerOf());
1316                                 while ( (n = (AliITSneuron*)innerof.Next()) ) {
1317                                         // if the examined neuron has not the largest activation
1318                                         // it is skipped and removed from array of all neurons
1319                                         // and of its outer point (its inner is the cursor p)
1320                                         if (n->Activation() < test) {
1321                                                 p->InnerOf()->Remove(n);
1322                                                 n->Outer()->OuterOf()->Remove(n);
1323                                                 delete fNeurons->Remove(n);
1324                                                 continue;
1325                                         }
1326                                         // otherwise, its activation becomes the maximum reference
1327                                         p->Next() = n->Outer();
1328                                         // at the exit of the while(), the fNext will point
1329                                         // to the outer node of the neuron starting in p, whose
1330                                         // activation is the largest.
1331                                 }
1332                                 // the same procedure is made now for all neurons
1333                                 // for which p is the outer point
1334                                 test = fActMinimum;
1335                                 p->Prev() = 0;
1336                                 TObjArrayIter outerof(p->OuterOf());
1337                                 while ( (n = (AliITSneuron*)outerof.Next()) ) {
1338                                         // if the examined neuron has not the largest activation
1339                                         // it is skipped and removed from array of all neurons
1340                                         // and of its inner point (its outer is the cursor p)
1341                                         if (n->Activation() < test) {
1342                                                 p->OuterOf()->Remove(n);
1343                                                 n->Inner()->InnerOf()->Remove(n);
1344                                                 delete fNeurons->Remove(n);
1345                                                 continue;
1346                                         }
1347                                         // otherwise, its activation becomes the maximum reference
1348                                         p->Prev() = n->Inner();
1349                                         // at the exit of the while(), the fPrev will point
1350                                         // to the inner node of the neuron ending in p, whose
1351                                         // activation is the largest.
1352                                 }
1353                         } // end while (p ...)
1354                 } // end for (itheta ...)
1355         } // end for (ilayer ...)
1356
1357         // now the mismatches are solved
1358         Bool_t matchPrev, matchNext;
1359         for (ilayer = 0; ilayer < fNLayers; ilayer++) {
1360                 for (itheta = 0; itheta < 180; itheta++) {
1361                         // get the array of point in a given layer/theta-slot/sector
1362                         lstSector = (TObjArray*)fNodes[ilayer][itheta]->At(sector);
1363                         TObjArrayIter nodes(lstSector);
1364                         while ( (p = (AliITSnode*)nodes.Next()) ) {
1365                                 // now p will point to a fPrev and a fNext node.
1366                                 // Ideally they are placed this way: fPrev --> P --> fNext
1367                                 // A mismatch happens if the point addressed as fPrev does NOT
1368                                 // point to p as its fNext. And the same for the point addressed
1369                                 // as fNext.
1370                                 // In this case, the fNext and fPrev pointers are set to NULL
1371                                 // and p is excluded from the reconstruction
1372                                 matchPrev = matchNext= kFALSE;
1373                                 if (ilayer > 0 && p->Prev() != NULL)
1374                                         if (p->Prev()->Next() == p) matchPrev = kTRUE;
1375                                 if (ilayer < 5 && p->Next() != NULL)
1376                                         if (p->Next()->Prev() == p) matchNext = kTRUE;
1377                                 if (ilayer == 0)
1378                                         matchPrev = kTRUE;
1379                                 else if (ilayer == 5)
1380                                         matchNext = kTRUE;
1381                                 if (!matchNext || !matchPrev) {
1382                                         p->Prev() = p->Next() = 0;
1383                                 }
1384                         } // end while (p ...)
1385                 } // end for (itheta ...)
1386         } // end for (ilayer ...)
1387 }
1388
1389 //__________________________________________________________________________________
1390 Int_t AliITStrackerANN::SaveTracks(Int_t sector)
1391 {
1392 /********************************************************************************
1393
1394         TRACK SAVING
1395         ------------
1396         Using the fNext and fPrev pointers, the chain is followed
1397         and the track is fitted and saved.
1398         Of course, the track is followed as a chain with a point
1399         for each layer, then the track following starts always
1400         from the clusters in layer 0.
1401
1402 ***********************************************************************************/
1403
1404         // if not initialized, the tracks TobjArray is initialized
1405         if (!fFoundTracks) fFoundTracks = new TObjArray;
1406
1407         // meaningful counters
1408         Int_t itheta, ilayer, l;
1409         TObjArray *lstSector = 0;
1410         AliITSnode *p = 0, *q = 0, **node = new AliITSnode*[fNLayers];
1411         for (l = 0; l < fNLayers; l++) node[l] = 0;
1412
1413         /*
1414         array = new AliITSnode*[fNLayers + 1];
1415         for (l = 0; l <= fNLayers; l++) array[l] = 0;
1416         array[0] = new AliITSnode();
1417         array[0]->X() = fVertexX;
1418         array[0]->Y() = fVertexY;
1419         array[0]->Z() = fVertexZ;
1420         array[0]->ErrX2() = fVertexErrorX;
1421         array[0]->ErrY2() = fVertexErrorY;
1422         array[0]->ErrZ2() = fVertexErrorZ;
1423         */
1424         Double_t *param = new Double_t[8];
1425
1426         // scan the whole collection of nodes
1427         for (ilayer = 0; ilayer < 1; ilayer++) {
1428                 for (itheta = 0; itheta < 180; itheta++) {
1429                         // get the array of point in a given layer/theta-slot/sector
1430                         lstSector = (TObjArray*)fNodes[ilayer][itheta]->At(sector);
1431                         TObjArrayIter nodes(lstSector);
1432                         while ( (p = (AliITSnode*)nodes.Next()) ) {
1433                                 for (q = p; q; q = q->Next()) {
1434                                         l = q->GetLayer();
1435                                         node[l] = q;
1436                                 }
1437                                 //if (!RiemannFit(fNLayers, node, param)) continue;
1438                                 // initialization of Kalman Filter Tracking
1439                                 AliITSclusterV2 *cluster = (AliITSclusterV2*)GetCluster(node[0]->ClusterRef());
1440                                 Int_t mod = cluster->GetDetectorIndex();
1441                                 Int_t lay, lad, det;
1442                                 fGeom->GetModuleId(mod, lay, lad, det);
1443                                 Float_t y0 = cluster->GetY();
1444                                 Float_t z0 = cluster->GetZ();
1445                                 AliITStrackSA* trac = new AliITStrackSA(lay, lad, det,
1446                                                                         y0, z0,
1447                                                                                                                                          param[4], param[7], param[3], 1);
1448                                 for (l = 0; l < fNLayers; l++) {
1449                                         cluster = (AliITSclusterV2*)GetCluster(node[l]->ClusterRef());
1450                                         if (cluster) trac->AddClusterV2(l, (node[l]->ClusterRef() & 0x0fffffff)>>0);
1451                                 }
1452             AliITStrackV2* ot = new AliITStrackV2(*trac);
1453                                 ot->ResetCovariance();
1454                                 ot->ResetClusters();
1455                                 if (RefitAt(49.,ot,trac)) { //fit from layer 1 to layer 6
1456                                         AliITStrackV2 *otrack2 = new AliITStrackV2(*ot);
1457                                         otrack2->ResetCovariance();
1458                                         otrack2->ResetClusters();
1459                                         //fit from layer 6 to layer 1
1460                                         if (RefitAt(3.7,otrack2,ot)) fFoundTracks->AddLast(otrack2);
1461                                 }
1462                                 // end of Kalman Filter fit
1463                         }
1464                 }
1465         }
1466
1467         return 1;
1468 }
1469
1470 //__________________________________________________________________________________
1471 void AliITStrackerANN::ExportTracks(const char *filename) const
1472 {
1473 // Exports found tracks into a TTree of AliITStrackV2 objects
1474         TFile *file = new TFile(filename, "RECREATE");
1475         TTree *tree = new TTree("TreeT-ANN", "Tracks found in ITS stand-alone with Neural Tracking");
1476         AliITStrackV2 *track = 0;
1477         tree->Branch("Tracks", &track, "AliITStrackV2");
1478         TObjArrayIter tracks(fFoundTracks);
1479         while ( (track = (AliITStrackV2*)tracks.Next()) ) {
1480                 tree->Fill();
1481         }
1482         file->cd();
1483         tree->Write();
1484         file->Close();
1485 }
1486
1487
1488 //__________________________________________________________________________________
1489 void AliITStrackerANN::CleanNetwork()
1490 {
1491         // Removes deactivated units from the network
1492
1493         AliITSneuron *unit = 0;
1494         TObjArrayIter neurons(fNeurons);
1495         while ( (unit = (AliITSneuron*)neurons.Next()) ) {
1496                 if (unit->Activation() < fActMinimum) {
1497                         unit->Inner()->InnerOf()->Remove(unit);
1498                         unit->Outer()->OuterOf()->Remove(unit);
1499                         delete fNeurons->Remove(unit);
1500                 }
1501         }
1502         return;
1503         Bool_t removed;
1504         Int_t nIn, nOut;
1505         AliITSneuron *enemy = 0;
1506         neurons.Reset();
1507         while ( (unit = (AliITSneuron*)neurons.Next()) ) {
1508                 nIn = (Int_t)unit->Inner()->InnerOf()->GetSize();
1509                 nOut = (Int_t)unit->Outer()->OuterOf()->GetSize();
1510                 if (nIn < 2 && nOut < 2) continue;
1511                 removed = kFALSE;
1512                 if (nIn > 1) {
1513                         TObjArrayIter competing(unit->Inner()->InnerOf());
1514                         while ( (enemy = (AliITSneuron*)competing.Next()) ) {
1515                                 if (unit->Activation() > enemy->Activation()) {
1516                                         enemy->Inner()->InnerOf()->Remove(enemy);
1517                                         enemy->Outer()->OuterOf()->Remove(enemy);
1518                                         delete fNeurons->Remove(enemy);
1519                                 }
1520                                 else {
1521                                         unit->Inner()->InnerOf()->Remove(unit);
1522                                         unit->Outer()->OuterOf()->Remove(unit);
1523                                         delete fNeurons->Remove(unit);
1524                                         removed = kTRUE;
1525                                         break;
1526                                 }
1527                         }
1528                         if (removed) continue;
1529                 }
1530                 if (nOut > 1) {
1531                         TObjArrayIter competing(unit->Outer()->OuterOf());
1532                         while ( (enemy = (AliITSneuron*)competing.Next()) ) {
1533                                 if (unit->Activation() > enemy->Activation()) {
1534                                         enemy->Inner()->InnerOf()->Remove(enemy);
1535                                         enemy->Outer()->OuterOf()->Remove(enemy);
1536                                         delete fNeurons->Remove(enemy);
1537                                 }
1538                                 else {
1539                                         unit->Inner()->InnerOf()->Remove(unit);
1540                                         unit->Outer()->OuterOf()->Remove(unit);
1541                                         delete fNeurons->Remove(unit);
1542                                         removed = kTRUE;
1543                                         break;
1544                                 }
1545                         }
1546                 }
1547         }
1548  }
1549
1550 //__________________________________________________________________________________
1551 Int_t AliITStrackerANN::StoreTracks()
1552 {
1553         // Stores the tracks found in a single neural tracking step.
1554         // In order to do this, it sects each neuron which has a point
1555         // in the first layer.
1556         // Then
1557
1558         // if not initialized, the tracks TobjArray is initialized
1559         if (!fFoundTracks) fFoundTracks = new TObjArray;
1560
1561         Int_t i, check, stored = 0;
1562         Double_t testAct = 0;
1563         AliITSneuron *unit = 0, *cursor = 0, *fwd = 0;
1564         AliITSnode *node = 0;
1565         TObjArrayIter iter(fNeurons), *fwdIter;
1566         TObjArray *removedUnits = new TObjArray(0);
1567         removedUnits->SetOwner(kFALSE);
1568         AliITStrackANN annTrack(fNLayers);
1569
1570         for (;;) {
1571                 unit = (AliITSneuron*)iter.Next();
1572                 if (!unit) break;
1573                 if (unit->Inner()->GetLayer() > 0) continue;
1574                 annTrack.SetNode(unit->Inner()->GetLayer(), unit->Inner());
1575                 annTrack.SetNode(unit->Outer()->GetLayer(), unit->Outer());
1576                 node = unit->Outer();
1577                 removedUnits->AddLast(unit);
1578                 while (node) {
1579                         testAct = fActMinimum;
1580                         fwdIter = (TObjArrayIter*)node->InnerOf()->MakeIterator();
1581                         fwd = 0;
1582                         for (;;) {
1583                                 cursor = (AliITSneuron*)fwdIter->Next();
1584                                 if (!cursor) break;
1585                                 if (cursor->Used()) continue;
1586                                 if (cursor->Activation() >= testAct) {
1587                                         testAct = cursor->Activation();
1588                                         fwd = cursor;
1589                                 }
1590                         }
1591                         if (!fwd) break;
1592                         removedUnits->AddLast(fwd);
1593                         node = fwd->Outer();
1594                         annTrack.SetNode(node->GetLayer(), node);
1595                 }
1596                 check = annTrack.CheckOccupation();
1597                 if (check >= 6) {
1598                         stored++;
1599                         // FIT
1600                         //if (!RiemannFit(fNLayers, trackitem, param)) continue;
1601                         if (!annTrack.RiemannFit()) continue;
1602                         // initialization of Kalman Filter Tracking
1603                         AliITSclusterV2 *cluster = (AliITSclusterV2*)GetCluster(annTrack[0]->ClusterRef());
1604                         Int_t mod = cluster->GetDetectorIndex();
1605                         Int_t lay, lad, det;
1606                         fGeom->GetModuleId(mod, lay, lad, det);
1607                         Float_t y0 = cluster->GetY();
1608                         Float_t z0 = cluster->GetZ();
1609                         AliITStrackSA* trac = new AliITStrackSA(lay, lad, det, y0, z0,
1610                                                                                                                                  annTrack.Phi(), annTrack.TanLambda(),
1611                                                                                                                                  annTrack.Curv(), 1);
1612                         for (Int_t l = 0; l < fNLayers; l++) {
1613                                 if (!annTrack[l]) continue;
1614                                 cluster = (AliITSclusterV2*)GetCluster(annTrack[l]->ClusterRef());
1615                                 if (cluster) trac->AddClusterV2(l, (annTrack[l]->ClusterRef() & 0x0fffffff)>>0);
1616                         }
1617                         AliITStrackV2* ot = new AliITStrackV2(*trac);
1618                         ot->ResetCovariance();
1619                         ot->ResetClusters();
1620                         if (RefitAt(49.,ot,trac)) { //fit from layer 1 to layer 6
1621                                 AliITStrackV2 *otrack2 = new AliITStrackV2(*ot);
1622                                 otrack2->ResetCovariance();
1623                                 otrack2->ResetClusters();
1624                                 //fit from layer 6 to layer 1
1625                                 if (RefitAt(3.7,otrack2,ot)) fFoundTracks->AddLast(otrack2);
1626                         }
1627                         // end of Kalman Filter fit
1628                         // END FIT
1629                         for (i = 0; i < fNLayers; i++) {
1630                                 //node = (AliITSnode*)removedPoints->At(i);
1631                                 //node->Use();
1632                                 annTrack[i]->Use();
1633                         }
1634                         fwdIter = (TObjArrayIter*)removedUnits->MakeIterator();
1635                         for (;;) {
1636                                 cursor = (AliITSneuron*)fwdIter->Next();
1637                                 if(!cursor) break;
1638                                 cursor->Used() = 1;
1639                         }
1640                 }
1641         }
1642
1643         return stored;
1644 }
1645
1646 Double_t AliITStrackerANN::Weight(AliITSneuron *nAB, AliITSneuron *nBC)
1647 {
1648 /***********************************************************************************
1649  * Calculation of neural weight.
1650  * The implementation of positive neural weight is set only in the case
1651  * of connected units (e.g.: A->B with B->C).
1652  * Given that B is the **common** point. care should be taken to pass
1653  * as FIRST argument the neuron going "to" B, and
1654  * as SECOND argument the neuron starting "from" B
1655  * anyway, a check is put in order to return 0.0 when arguments are not well given.
1656  ***********************************************************************************/
1657
1658         if (nAB->Outer() != nBC->Inner()) {
1659                 if (nBC->Outer() == nAB->Inner()) {
1660                         AliITSneuron *temp = nAB;
1661                         nAB = nBC;
1662                         nBC = temp;
1663                         temp = 0;
1664                         if (fMsgLevel >= 3) {
1665                                 Info("Weight", "Switching wrongly ordered arguments.");
1666                         }
1667                 }
1668                 Warning("Weight", "Not connected segments. Returning 0.0");
1669                 return 0.0;
1670         }
1671
1672         AliITSnode *pA = nAB->Inner();
1673         AliITSnode *pB = nAB->Outer();
1674         AliITSnode *pC = nBC->Outer();
1675
1676         TVector3 vAB(pB->X() - pA->X(), pB->Y() - pA->Y(), pB->Z() - pA->Z());
1677         TVector3 vBC(pC->X() - pB->X(), pC->Y() - pB->Y(), pC->Z() - pB->Z());
1678
1679         Double_t weight = 1.0 - sin(vAB.Angle(vBC));
1680         return fGain2CostRatio * TMath::Power(weight, fExponent);
1681 }
1682
1683
1684
1685 /******************************************
1686  ******************************************
1687  *** AliITStrackerANN::AliITSnode class ***
1688  ******************************************
1689  ******************************************/
1690
1691 //__________________________________________________________________________________
1692 AliITStrackerANN::AliITSnode::AliITSnode()
1693 : fUsed(kFALSE), fClusterRef(-1),
1694   fMatches(NULL), fInnerOf(NULL), fOuterOf(NULL),
1695   fNext(NULL), fPrev(NULL)
1696 {
1697         // Constructor for the embedded 'AliITSnode' class.
1698         // It initializes all pointer-like objects.
1699
1700         fX = fY = fZ = 0.0;
1701         fEX2 = fEY2 = fEZ2 = 0.0;
1702 }
1703
1704 //__________________________________________________________________________________
1705 AliITStrackerANN::AliITSnode::~AliITSnode()
1706 {
1707         // Destructor for the embedded 'AliITSnode' class.
1708         // It should clear the object arrays, but it is possible
1709         // that some objects still are useful after the point deletion
1710         // then the arrays are cleared but their objects are owed by
1711         // another TCollection object, and not deleted.
1712         // For safety reasons, all the pointers are set to zero.
1713
1714         fInnerOf->SetOwner(kFALSE);
1715         fInnerOf->Clear();
1716         delete fInnerOf;
1717         fInnerOf = 0;
1718         fOuterOf->SetOwner(kFALSE);
1719         fOuterOf->Clear();
1720         delete fOuterOf;
1721         fOuterOf = 0;
1722         fMatches->SetOwner(kFALSE);
1723         fMatches->Clear();
1724         delete fMatches;
1725         fMatches = 0;
1726         fNext = 0;
1727         fPrev = 0;
1728 }
1729
1730 //__________________________________________________________________________________
1731 Double_t AliITStrackerANN::AliITSnode::GetPhi() const
1732 {
1733         // Calculates the 'phi' (azimutal) angle, and returns it
1734         // in the range between 0 and 2Pi radians.
1735
1736         Double_t q;
1737         q = TMath::ATan2(fY,fX);
1738         if (q >= 0.)
1739                 return q;
1740         else
1741                 return q + TMath::TwoPi();
1742 }
1743
1744 //__________________________________________________________________________________
1745 Double_t AliITStrackerANN::AliITSnode::GetError(Option_t *option)
1746 {
1747         // Returns the error or the square error of
1748         // values related to the coordinates in different systems.
1749         // The option argument specifies the coordinate error desired:
1750         //
1751         // "R2"     --> error in transverse radius
1752         // "R3"     --> error in spherical radius
1753         // "PHI"    --> error in azimuthal angle
1754         // "THETA"  --> error in polar angle
1755         // "SQ"     --> get the square of error
1756         //
1757         // In order to get the error on the cartesian coordinates
1758         // reference to the inline ErrX2(), ErrY2() adn ErrZ2() methods.
1759
1760         TString opt(option);
1761         Double_t errorSq = 0.0;
1762         opt.ToUpper();
1763
1764         if (opt.Contains("R2")) {
1765                 errorSq  = fX*fX*fEX2 + fY*fY*fEY2;
1766                 errorSq /= GetR2sq();
1767         }
1768         else if (opt.Contains("R3")) {
1769                 errorSq  = fX*fX*fEX2 + fY*fY*fEY2 + fZ*fZ*fEZ2;
1770                 errorSq /= GetR3sq();
1771         }
1772         else if (opt.Contains("PHI")) {
1773                 errorSq  = fY*fY*fEX2;
1774                 errorSq += fX*fX*fEY2;
1775                 errorSq /= GetR2sq() * GetR2sq();
1776         }
1777         else if (opt.Contains("THETA")) {
1778                 errorSq = fZ*fZ * (fX*fX*fEX2 + fY*fY*fEY2);
1779                 errorSq += GetR2sq() * GetR2sq() * fEZ2;
1780                 errorSq /= GetR3sq() * GetR3sq() * GetR2() * GetR2();
1781         }
1782
1783         if (!opt.Contains("SQ"))
1784                 return TMath::Sqrt(errorSq);
1785         else
1786                 return errorSq;
1787 }
1788
1789
1790
1791 /********************************************
1792  ********************************************
1793  *** AliITStrackerANN::AliITSneuron class ***
1794  ********************************************
1795  ********************************************/
1796
1797 //__________________________________________________________________________________
1798 AliITStrackerANN::AliITSneuron::AliITSneuron
1799 (AliITSnode *inner, AliITSnode *outer, Double_t minAct, Double_t maxAct)
1800         : fUsed(0), fInner(inner), fOuter(outer)
1801 {
1802         // Default neuron constructor
1803         fActivation = gRandom->Rndm() * (maxAct-minAct) + minAct;
1804         fGain = new TObjArray;
1805 }
1806
1807 //__________________________________________________________________________________
1808 Double_t AliITStrackerANN::AliITSneuron::Activate(Double_t temperature)
1809 {
1810         // This computes the new activation of a neuron, and returns
1811         // its activation variation as a consequence of the updating.
1812         //
1813         // Arguments:
1814         // - the 'temperature' parameter for the neural activation logistic function
1815         //
1816         // Operations:
1817         // - collects the total gain, by summing the products
1818         //   of the activation of each sequenced unit by the relative weight.
1819         // - collects the total cost, by summing the activations of
1820         //   all competing units
1821         // - passes the sum of gain - cost to the activation function and
1822         //   calculates the new activation
1823         //
1824         // Return value:
1825         // - the difference between the old activation and the new one
1826         //   (absolute value)
1827
1828         // local variables
1829         Double_t sumGain = 0.0;      // total contribution from chained neurons
1830         Double_t sumCost = 0.0;      // total contribution from crossing neurons
1831         Double_t input;              // total input
1832         Double_t actOld, actNew;     // old and new values for the activation
1833         AliITSneuron *linked = 0;    // cursor for scanning the neuron arrays (for link check)
1834         AliITSlink *link;            // cursor for scanning the synapses arrays (for link check)
1835         TObjArrayIter *iterator = 0; // pointer to the iterator along the neuron arrays
1836
1837         // sum contributions from the correlated units
1838         iterator = (TObjArrayIter*)fGain->MakeIterator();
1839         for(;;) {
1840                 link = (AliITSlink*)iterator->Next();
1841                 if (!link) break;
1842                 sumGain += link->Contribution();
1843         }
1844         delete iterator;
1845
1846         // sum contributions from the competing units:
1847         // the ones which have the same starting point...
1848         iterator = (TObjArrayIter*)fInner->InnerOf()->MakeIterator();
1849         for (;;) {
1850                 linked = (AliITSneuron*)iterator->Next();
1851                 if (!linked) break;
1852                 if (linked == this) continue;
1853                 sumCost += linked->fActivation;
1854         }
1855         delete iterator;
1856         // ...and the ones which have the same ending point
1857         iterator = (TObjArrayIter*)fOuter->OuterOf()->MakeIterator();
1858         for (;;) {
1859                 linked = (AliITSneuron*)iterator->Next();
1860                 if (!linked) break;
1861                 if (linked == this) continue;
1862                 sumCost += linked->fActivation;
1863         }
1864
1865         // calculate the total input as the difference between gain and cost
1866         input = (sumGain - sumCost) / temperature;
1867         actOld = fActivation;
1868         // calculate the final output
1869 #ifdef NEURAL_LINEAR
1870         if (input <= -2.0 * temperature)
1871                 actNew = 0.0;
1872         else if (input >= 2.0 * temperature)
1873                 actNew = 1.0;
1874         else
1875                 actNew = input / (4.0 * temperature) + 0.5;
1876 #else
1877         actNew = 1.0 / (1.0 + TMath::Exp(-input));
1878 #endif
1879         fActivation = actNew;
1880
1881         // return the activation variation
1882         return TMath::Abs(actNew - actOld);
1883 }
1884
1885
1886
1887 /******************************************
1888  ******************************************
1889  *** AliITStrackerANN::AliITSlink class ***
1890  ******************************************
1891  ******************************************/
1892
1893  // No methods defined non-inline
1894
1895
1896
1897  /**********************************************
1898   **********************************************
1899   *** AliITStrackerANN::AliITStrackANN class ***
1900   **********************************************
1901   **********************************************/
1902
1903 //__________________________________________________________________________________
1904 AliITStrackerANN::AliITStrackANN::AliITStrackANN(Int_t dim) : fNPoints(dim)
1905 {
1906         // Default constructor for the AliITStrackANN class
1907
1908         fXCenter = 0.0;
1909         fYCenter = 0.0;
1910         fRadius = 0.0;
1911         fCurv = 0.0;
1912         fDTrans = 0.0;
1913         fDLong = 0.0;
1914         fTanLambda = 0.0;
1915
1916         if (! dim) {
1917                 fNode = 0;
1918         }
1919         else{
1920                 Int_t i = 0;
1921                 fNode = new AliITSnode*[dim];
1922                 for (i = 0; i < dim; i++) fNode[i] = 0;
1923         }
1924 }
1925
1926 //__________________________________________________________________________________
1927 Int_t AliITStrackerANN::AliITStrackANN::CheckOccupation() const
1928 {
1929         // Returns the number of pointers fNode which are not NULL
1930
1931         Int_t i;         // cursor
1932         Int_t count = 0; // counter for how many points are stored in the track
1933
1934         for (i = 0; i < fNPoints; i++) {
1935                 if (fNode[i] != NULL) count++;
1936         }
1937
1938         return count;
1939 }
1940
1941 //__________________________________________________________________________________
1942 Bool_t AliITStrackerANN::AliITStrackANN::RiemannFit()
1943 {
1944         // Performs the Riemann Sphere fit for the given points to a circle
1945         // and then uses the fit parameters to fit a helix in space.
1946         //
1947         // Return values:
1948         // - kTRUE if all operations have been performed
1949         // - kFALSE if the numbers risk to lead to an arithmetic violation
1950
1951         Int_t i, j, count, dim = fNPoints;
1952
1953         // First check for all points
1954         count = CheckOccupation();
1955         if (count != fNPoints) {
1956                 Error ("AliITStrackANN::RiemannFit", "CheckOccupations returns %d, fNPoints = %d ==> MISMATCH", count, fNPoints);
1957                 return kFALSE;
1958         }
1959
1960         // matrix of ones
1961         TMatrixD m1(dim,1);
1962         for (i = 0; i < dim; i++) m1(i,0) = 1.0;
1963
1964         // matrix of Rieman projection coordinates
1965         TMatrixD coords(dim,3);
1966         for (i = 0; i < dim; i++) {
1967                 coords(i,0) = fNode[i]->X();
1968                 coords(i,1) = fNode[i]->Y();
1969                 coords(i,2) = fNode[i]->GetR2sq();
1970         }
1971
1972         // matrix of weights
1973         Double_t xterm, yterm, ex, ey;
1974         TMatrixD weights(dim,dim);
1975         for (i = 0; i < dim; i++) {
1976                 xterm = fNode[i]->X() * fNode[i]->GetPhi() - fNode[i]->Y() / fNode[i]->GetR2();
1977                 ex = fNode[i]->ErrX2();
1978                 yterm = fNode[i]->Y() * fNode[i]->GetPhi() + fNode[i]->X() / fNode[i]->GetR2();
1979                 ey = fNode[i]->ErrY2();
1980                 weights(i,i) = fNode[i]->GetR2sq() / (xterm * xterm * ex + yterm * yterm * ey );
1981         }
1982
1983         // weighted sample mean
1984         Double_t meanX = 0.0, meanY = 0.0, meanW = 0.0, sw = 0.0;
1985         for (i = 0; i < dim; i++) {
1986                 meanX += weights(i,i) * coords(i,0);
1987                 meanY += weights(i,i) * coords(i,1);
1988                 meanW += weights(i,i) * coords(i,2);
1989                 sw += weights(i,i);
1990         }
1991         meanX /= sw;
1992         meanY /= sw;
1993         meanW /= sw;
1994
1995         // sample covariance matrix
1996         for (i = 0; i < dim; i++) {
1997                 coords(i,0) -= meanX;
1998                 coords(i,1) -= meanY;
1999                 coords(i,2) -= meanW;
2000         }
2001         TMatrixD coordsT(TMatrixD::kTransposed, coords);
2002         TMatrixD weights4coords(weights, TMatrixD::kMult, coords);
2003         TMatrixD sampleCov(coordsT, TMatrixD::kMult, weights4coords);
2004         for (i = 0; i < 3; i++) {
2005                 for (j = i + 1; j < 3; j++) {
2006                         sampleCov(i,j)  = sampleCov(j,i)  = (sampleCov(i,j) + sampleCov(j,i)) * 0.5;
2007                 }
2008         }
2009
2010         // Eigenvalue problem solving for V matrix
2011         Int_t ileast = 0;
2012         TVectorD eval(3), n(3);
2013 #if ROOT_VERSION_CODE < ROOT_VERSION(4,0,2)
2014         TMatrixD evec = sampleCov.EigenVectors(eval);
2015 #else
2016         TMatrixDEigen ei(sampleCov);
2017         TMatrixD evec = ei.GetEigenVectors();
2018         eval = ei.GetEigenValuesRe();
2019 #endif
2020         if (eval(1) < eval(ileast)) ileast = 1;
2021         if (eval(2) < eval(ileast)) ileast = 2;
2022         n(0) = evec(0, ileast);
2023         n(1) = evec(1, ileast);
2024         n(2) = evec(2, ileast);
2025
2026         // c - known term in the plane intersection with Riemann axes
2027         Double_t c = -(meanX * n(0) + meanY * n(1) + meanW * n(2));
2028
2029         // center and radius of fitted circle
2030         Double_t xc, yc, radius, curv;
2031         xc = -n(0) / (2. * n(2));
2032         yc = -n(1) / (2. * n(2));
2033         radius = (1. - n(2)*n(2) - 4.*c*n(2)) / (4. * n(2) * n(2));
2034
2035         if (radius <= 0.E0) {
2036                 Error("RiemannFit", "Radius = %f less than zero!!!", radius);
2037                 return kFALSE;
2038         }
2039         radius = TMath::Sqrt(radius);
2040         curv = 1.0 / radius;
2041
2042         // evaluating signs for curvature and others
2043         Double_t phi1 = 0.0, phi2, temp1, temp2, phi0, sumdphi = 0.0;
2044         AliITSnode *p = fNode[0];
2045         phi1 = p->GetPhi();
2046         for (i = 1; i < dim; i++) {
2047                 p = (AliITSnode*)fNode[i];
2048                 if (!p) break;
2049                 phi2 = p->GetPhi();
2050                 temp1 = phi1;
2051                 temp2 = phi2;
2052                 if (temp1 > fgkPi && temp2 < fgkPi)
2053                         temp2 += fgkTwoPi;
2054                 else if (temp1 < fgkPi && temp2 > fgkPi)
2055                         temp1 += fgkTwoPi;
2056                 sumdphi += temp2 - temp1;
2057                 phi1 = phi2;
2058         }
2059         if (sumdphi < 0.E0) curv = -curv;
2060         Double_t diff, angle = TMath::ATan2(yc, xc);
2061         if (curv < 0.E0)
2062                 phi0 = angle + 0.5 * TMath::Pi();
2063         else
2064                 phi0 = angle - 0.5 * TMath::Pi();
2065         diff = angle - phi0;
2066
2067         Double_t dt, temp = TMath::Sqrt(xc*xc + yc*yc) - radius;
2068         if (curv >= 0.E0)
2069                 dt = temp;
2070         else
2071                 dt = -temp;
2072         //cout << "Dt = " << dt << endl;
2073
2074         Double_t halfC = 0.5 * curv, test;
2075         Double_t *s = new Double_t[dim], *zz = new Double_t[dim], *ws = new Double_t[dim];
2076         for (j = 0; j < 6; j++) {
2077                 p = fNode[j];
2078                 if (!p) break;
2079                 //----
2080                 s[j] = (p->GetR2sq() - dt * dt) / (1. + curv * dt);
2081                 if (s[j] < 0.) {
2082                         if (TMath::Abs(s[j]) < 1.E-6) s[j] = 0.;
2083                         else {
2084                                 Error("RiemannFit", "Square root argument error: %17.15g < 0", s[j]);
2085                                 return kFALSE;
2086                         }
2087                 }
2088                 s[j] = TMath::Sqrt(s[j]);
2089                 //cout << "Curv = " << halfC << " --- s[" << j << "] = " << s[j] << endl;
2090                 s[j] *= halfC;
2091                 test = TMath::Abs(s[j]);
2092                 if (test > 1.) {
2093                         if (test <= 1.1)
2094                                 s[j] = ((s[j] > 0.) ? 0.99999999999 : -0.9999999999);
2095                         else {
2096                                 Error("RiemannFit", "Value too large: %17.15g", s[j]);
2097                                 return kFALSE;
2098                         }
2099                 }
2100                 //----
2101                 zz[j] = p->Z();
2102                 s[j] = TMath::ASin(s[j]) / halfC;
2103                 ws[j] = 1.0 / (p->ErrZ2());
2104         }
2105
2106         // second tep final fit
2107         Double_t s2Sum = 0.0, zSum = 0.0, szSum = 0.0, sSum = 0.0, sumw = 0.0;
2108         for (i = 0; i < dim; i++) {
2109                 s2Sum += ws[i] * s[i] * s[i];
2110                 zSum  += ws[i] * zz[i];
2111                 sSum  += ws[i] * s[i];
2112                 szSum += ws[i] * s[i] * zz[i];
2113                 sumw += ws[i];
2114         }
2115         s2Sum /= sumw;
2116         zSum /= sumw;
2117         sSum /= sumw;
2118         szSum /= sumw;
2119         temp = s2Sum - sSum*sSum;
2120
2121         Double_t dz, tanL;
2122         dz = (s2Sum*zSum - sSum*szSum) / temp;
2123         tanL = (szSum - sSum*zSum) / temp;
2124
2125         fXCenter = xc;
2126         fYCenter = yc;
2127         fRadius = radius;
2128         fCurv = curv;
2129         fPhi = phi0;
2130         fDTrans = dt;
2131         fDLong = dz;
2132         fTanLambda = tanL;
2133
2134         delete [] s;
2135         delete [] zz;
2136         delete [] ws;
2137
2138         return kTRUE;
2139 }