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