ITS/AliITSNeuralTracker.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: Alberto Pulvirenti.                                            *
   5  *                                                                        *
   6  * Permission to use, copy, modify and distribute this software and its   *
   7  * documentation strictly for non-commercial purposes is hereby granted   *
   8  * without fee, provided that the above copyright notice appears in all   *
   9  * copies and that both the copyright notice and this permission notice   *
  10  * appear in the supporting documentation. The authors make no claims     *
  11  * about the suitability of this software for any purpose.                *
  12  * It is provided "as is" without express or implied warranty.            *
  13  *                                                                        *
  14  * AN ITS STAND-ALONE "NEURAL" TRACK FINDER                               *
  15  * ----------------------------------------                               *
  16  * This class implements the Denby-Peterson algorithm for track finding   *
  17  * in the ITS stand-alone, by means of a neural network simulation.       *
  18  * Many parameters have to be set for the neural network to operate       *
  19  * correctly and with a good efficiency.                                  *
  20  * The neural tracker must be feeded with a TTree filled with objects     *
  21  * of the class "AliITSNeuralPoint", into a single branch called       *
  22  * "Points".                                                              *
  23  **************************************************************************/
  24
  25 //#define NEURAL_LINEAR
  26 #include <fstream>
  27 #include <Riostream.h>
  28 #include <stdlib.h>
  29
  30 #include <TROOT.h>
  31 #include <TFile.h>
  32 #include <TTree.h>
  33 #include <TMath.h>
  34 #include <TLine.h>
  35 #include <TMarker.h>
  36 #include <TRandom.h>
  37 #include <TString.h>
  38 #include <TCanvas.h>
  39 #include <TVector3.h>
  40 #include <TParticle.h>
  41 #include <TObjArray.h>
  42 #include <TList.h>
  43
  44 #include "AliITSNeuralPoint.h"
  45 #include "AliITSNeuralTracker.h"
  46
  47 using namespace std;
  48 ClassImp(AliITSNeuralTracker)
  49
  50 //--------------------------------------------------------------------------------------------
  51
  52 AliITSNeuralTracker::AliITSNeuralTracker()
  53 {
  54         // CONSTRUCTOR
  55         //
  56         // Initializes some data-members:
  57         // - all pointers to NULL
  58         // - theta cut to 180 deg. (= no theta cut)
  59         // - initial choice of only 1 azimuthal sector
  60         // - initial values for the neural network parameters.
  61         //
  62         // With these settings the neural tracker can't work
  63         // because it has not any curvature cut set.
  64
  65         fCurvNum = 0;
  66         fCurvCut = 0;
  67
  68         fSectorNum   = 1;
  69         fSectorWidth = TMath::Pi() * 2.0;
  70         fPolarInterval = 45.0;
  71
  72         fStructureOK = kFALSE;
  73
  74         fVX = fVY = fVZ = 0.0;
  75
  76         Int_t ilayer, itheta;
  77         for (ilayer = 0; ilayer < 6; ilayer++) {
  78                 for (itheta = 0; itheta < 180; itheta++) fPoints[ilayer][itheta] = 0;
  79                 fThetaCut2DMin[ilayer] = 0.0;
  80                 fThetaCut2DMax[ilayer] = TMath::Pi();
  81                 fThetaCut3DMin[ilayer] = 0.0;
  82                 fThetaCut3DMax[ilayer] = TMath::Pi();
  83                 fHelixMatchCutMin[ilayer] = 1.0;
  84                 fHelixMatchCutMax[ilayer] = 1.0;
  85         }
  86
  87         fEdge1 = 0.3;
  88         fEdge2=  0.7;
  89
  90         fTemperature = 1.0;
  91         fStabThreshold = 0.001;
  92         fGain2CostRatio = 1.0;
  93         fAlignExponent = 1.0;
  94         fActMinimum = 0.5;
  95
  96         fNeurons = 0;
  97
  98         fChains = new TTree("TreeC", "Sextines of points");
  99         fChains->Branch("l0", &fPoint[0], "l0/I");
 100         fChains->Branch("l1", &fPoint[1], "l1/I");
 101         fChains->Branch("l2", &fPoint[2], "l2/I");
 102         fChains->Branch("l3", &fPoint[3], "l3/I");
 103         fChains->Branch("l4", &fPoint[4], "l4/I");
 104         fChains->Branch("l5", &fPoint[5], "l5/I");
 105 }
 106 //
 107 //--------------------------------------------------------------------------------------------
 108 //
 109 AliITSNeuralTracker::~AliITSNeuralTracker()
 110 {
 111         // DESTRUCTOR
 112         //
 113         // It Destroys all the dynamic arrays and
 114         // clears the TCollections and the points tree
 115
 116         cout << "Starting destructor..." << endl;
 117
 118         delete [] fCurvCut;
 119
 120         Int_t ilayer, itheta;
 121         if (fStructureOK) {
 122                 cout << "Deleting points..." << endl;
 123                 for (ilayer = 0; ilayer < 6; ilayer++) {
 124                         for (itheta = 0; itheta < 180; itheta++) {
 125                                 fPoints[ilayer][itheta]->SetOwner();
 126                                 delete fPoints[ilayer][itheta];
 127                         }
 128                 }
 129                 cout << "Deleting neurons..." << endl;
 130                 fNeurons->SetOwner();
 131                 delete fNeurons;
 132         }
 133
 134         cout << "AliITSNeuralTracker destructed completely!" << endl;
 135 }
 136 //
 137 //--------------------------------------------------------------------------------------------
 138 //
 139 void AliITSNeuralTracker::Display(TCanvas*& canv) const
 140 {
 141         // Displays the neural segments
 142
 143         Double_t x1, y1, x2, y2;
 144         canv->Clear();
 145         TObjArrayIter iter(fNeurons);
 146         for (;;) {
 147                 AliITSneuron *unit = (AliITSneuron*)iter.Next();
 148                 if (!unit) break;
 149                 if (unit->fActivation < fActMinimum) continue;
 150                 x1 = unit->fInner->X();
 151                 x2 = unit->fOuter->X();
 152                 y1 = unit->fInner->Y();
 153                 y2 = unit->fOuter->Y();
 154                 TLine *line = new TLine(x1, y1, x2, y2);
 155                 canv->cd();
 156                 line->Draw();
 157         }
 158         canv->Update();
 159 }
 160 //
 161 //--------------------------------------------------------------------------------------------
 162 //
 163 void AliITSNeuralTracker::SetThetaCuts2D(Double_t *min, Double_t *max)
 164 {
 165         // Cut setter
 166
 167         Int_t i;
 168         Double_t temp;
 169         for (i = 0; i < 5; i++) {
 170                 if (min[i] > max[i]) {
 171                         temp = min[i];
 172                         min[i] = max[i];
 173                         max[i] = temp;
 174                 }
 175                 fThetaCut2DMin[i] = min[i];
 176                 fThetaCut2DMax[i] = max[i];
 177         }
 178         for (i = 0; i < 5; i++) {
 179                 cout << "Theta 2D cut for layer " << i << " = " << fThetaCut2DMin[i] << " to " << fThetaCut2DMax[i] << endl;
 180         }
 181 }
 182 //
 183 //--------------------------------------------------------------------------------------------
 184 //
 185 void AliITSNeuralTracker::SetThetaCuts3D(Double_t *min, Double_t *max)
 186 {
 187         // Cut setter
 188
 189         Int_t i;
 190         Double_t temp;
 191         for (i = 0; i < 5; i++) {
 192                 if (min[i] > max[i]) {
 193                         temp = min[i];
 194                         min[i] = max[i];
 195                         max[i] = temp;
 196                 }
 197                 fThetaCut3DMin[i] = min[i];
 198                 fThetaCut3DMax[i] = max[i];
 199         }
 200         for (i = 0; i < 5; i++) {
 201                 cout << "Theta 3D cut for layer " << i << " = " << fThetaCut3DMin[i] << " to " << fThetaCut3DMax[i] << endl;
 202         }
 203 }
 204 //
 205 //--------------------------------------------------------------------------------------------
 206 //
 207 void AliITSNeuralTracker::SetHelixMatchCuts(Double_t *min, Double_t *max)
 208 {
 209         // Cut setter
 210
 211         Int_t i;
 212         Double_t temp;
 213         for (i = 0; i < 5; i++) {
 214                 if (min[i] > max[i]) {
 215                         temp = min[i];
 216                         min[i] = max[i];
 217                         max[i] = temp;
 218                 }
 219                 fHelixMatchCutMin[i] = min[i];
 220                 fHelixMatchCutMax[i] = max[i];
 221         }
 222         for (i = 0; i < 5; i++) {
 223                 cout << "Helix-match cut for layer " << i << " = " << fHelixMatchCutMin[i] << " to " << fHelixMatchCutMax[i] << endl;
 224         }
 225 }
 226 //
 227 //--------------------------------------------------------------------------------------------
 228 //
 229 void AliITSNeuralTracker::SetCurvatureCuts(Int_t ncuts, Double_t *curv)
 230 {
 231         // CURVATURE CUTS SETTER
 232         //
 233         // Requires an array of double values and its dimension
 234         // After sorting it in increasing order, the array of curvature cuts
 235         // is dinamically allocated, and filled with the sorted cuts array.
 236         // A message is shown which lists all the curvature cuts.
 237
 238         Int_t i, *ind = new Int_t[ncuts];
 239         TMath::Sort(ncuts, curv, ind, kFALSE);
 240         fCurvCut = new Double_t[ncuts];
 241         cout << "\n" << ncuts << " curvature cuts defined" << endl << "-----" << endl;
 242         for (i = 0; i < ncuts; i++) {
 243                 fCurvCut[i] = curv[ind[i]];
 244                 cout << "Cut #" << i + 1 << " : " << fCurvCut[i] << endl;
 245         }
 246         cout << "-----" << endl;
 247         fCurvNum = ncuts;
 248 }
 249 //
 250 //--------------------------------------------------------------------------------------------
 251 //
 252 void AliITSNeuralTracker::CreateArrayStructure(Int_t nsectors)
 253 {
 254         // ARRAY CREATOR
 255         //
 256         // Organizes the array structure to store all points in.
 257         //
 258         // The array is organized into a "multi-level" TCollection:
 259         // - 6 fPoints[] TObjArray containing a TObjArray for each layer
 260         // - each TObject contains a TObjArray for each sector.
 261
 262         // sets the number of sectors and their width.
 263         fSectorNum   = nsectors;
 264         fSectorWidth = TMath::Pi() * 2.0 / (Double_t)fSectorNum;
 265
 266         // creates the TCollection structure
 267         Int_t ilayer, isector, itheta;
 268         TObjArray *sector = 0;
 269         for (ilayer = 0; ilayer < 6; ilayer++) {
 270                 for (itheta = 0; itheta < 180; itheta++) {
 271                         if (fPoints[ilayer][itheta]) {
 272                                 fPoints[ilayer][itheta]->SetOwner();
 273                                 delete fPoints[ilayer][itheta];
 274                         }
 275                         fPoints[ilayer][itheta] = new TObjArray(nsectors);
 276                         for (isector = 0; isector < nsectors; isector++) {
 277                                 sector = new TObjArray;
 278                                 sector->SetOwner();
 279                                 fPoints[ilayer][itheta]->AddAt(sector, isector);
 280                         }
 281                 }
 282         }
 283
 284         // Sets a checking flag to TRUE.
 285         // This flag is checked before filling up the arrays with the points.
 286         fStructureOK = kTRUE;
 287 }
 288 //
 289 //--------------------------------------------------------------------------------------------
 290 //
 291 Int_t AliITSNeuralTracker::ArrangePoints(TTree* ptstree)
 292 {
 293         // POINTS STORAGE INTO ARRAY
 294         //
 295         // Reads points from the tree and creates AliITSNode objects for each one,
 296         // storing them into the array structure defined above.
 297         // Returns the number of points collected (if successful) or 0 (otherwise)
 298
 299         // check: if the points tree is NULL or empty, there is nothing to do...
 300         if ( !ptstree || (ptstree && !(Int_t)ptstree->GetEntries()) ) {
 301                 Error("ArrangePoints", "Points tree is NULL or empty: no points to arrange");
 302                 return 0;
 303         }
 304
 305         if (!fStructureOK) {
 306                 Error("ArrangePoints", "Structure NOT defined. Call CreateArrayStructure() first");
 307                 return 0;
 308         }
 309
 310         Int_t isector, itheta, ientry, ilayer, nentries, pos;
 311         TObjArray *sector = 0;
 312         AliITSNode *created = 0;
 313         AliITSNeuralPoint *cursor = 0;
 314
 315         ptstree->SetBranchAddress("pos", &pos);
 316         ptstree->SetBranchAddress("Points", &cursor);
 317         nentries = (Int_t)ptstree->GetEntries();
 318
 319         for (ientry = 0; ientry < nentries; ientry++) {
 320                 ptstree->GetEntry(ientry);
 321                 // creates the object
 322                 created = new AliITSNode(cursor, kTRUE);
 323                 created->SetUser(-1);
 324                 created->fPosInTree = pos;
 325                 // finds the sector in phi
 326                 isector = created->GetSector(fSectorWidth);
 327                 itheta  = created->GetThetaCell();
 328                 ilayer  = created->GetLayer();
 329                 if (ilayer < 0 || ilayer > 5) {
 330                         Error("ArrangePoints", "Layer value %d not allowed. Aborting...", ilayer);
 331                         return 0;
 332                 }
 333                 // selects the right TObjArray to store object into
 334                 sector = (TObjArray*)fPoints[ilayer][itheta]->At(isector);
 335                 sector->AddLast(created);
 336         }
 337
 338         // returns the number of points processed
 339         return ientry;
 340 }
 341 //
 342 //--------------------------------------------------------------------------------------------
 343 //
 344 void AliITSNeuralTracker::PrintPoints()
 345 {
 346         // creates the TCollection structure
 347         TObjArray *sector = 0;
 348         Int_t ilayer, isector, itheta;
 349         fstream file("test.log", ios::out);
 350         for (ilayer = 0; ilayer < 6; ilayer++) {
 351                 for (isector = 0; isector < fSectorNum; isector++) {
 352                         for (itheta = 0; itheta < 180; itheta++) {
 353                                 sector = (TObjArray*)fPoints[ilayer][itheta]->At(isector);
 354                                 file << ilayer << " " << isector << " " << itheta;
 355                                 file << " " << sector->GetSize() << " points" << endl;
 356                         }
 357                 }
 358         }
 359         file.close();
 360 }
 361 //
 362 //--------------------------------------------------------------------------------------------
 363 //
 364 Bool_t AliITSNeuralTracker::PassCurvCut
 365 (AliITSNode *p1, AliITSNode *p2, Int_t curvindex, Double_t vx, Double_t vy, Double_t vz)
 366 {
 367         // CURVATURE CUT EVALUATOR
 368         //
 369         // Checks the passsed point pair w.r. to the current curvature cut
 370         // Returns the result of the check.
 371
 372         if (curvindex < 0 || curvindex >= fCurvNum) {
 373                 Error("PassCurvCut", "Curv index %d out of range", curvindex);
 374                 return kFALSE;
 375         }
 376
 377         // Find the reference layer
 378         Int_t lay1 = p1->GetLayer();
 379         Int_t lay2 = p2->GetLayer();
 380         Int_t reflayer = (lay1 < lay2) ? lay1 : lay2;
 381
 382         Double_t x1 = p1->X() - vx;
 383         Double_t x2 = p2->X() - vx;
 384         Double_t y1 = p1->Y() - vy;
 385         Double_t y2 = p2->Y() - vy;
 386         Double_t z1 = p1->Z() - vz;
 387         Double_t z2 = p2->Z() - vz;
 388         Double_t r1 = sqrt(x1*x1 + y1*y1);
 389         Double_t r2 = sqrt(x2*x2 + y2*y2);
 390
 391         // calculation of curvature
 392         Double_t dx = p1->X() - p2->X(), dy = p1->Y() - p2->Y();
 393         Double_t num = 2 * (x1*y2 - x2*y1);
 394         Double_t den = r1*r2*sqrt(dx*dx + dy*dy);
 395         Double_t curv = 0.;
 396         /* FOR OLD VERSION
 397         if (den != 0.) {
 398                 curv = fabs(num / den);
 399                 if (curv > fCurvCut[curvindex]) return kFALSE;
 400                 return kTRUE;
 401         }
 402         else
 403                 return kFALSE;
 404         */
 405         // NEW VERSION
 406         if (den != 0.) {
 407                 curv = fabs(num / den);
 408                 if (curv > fCurvCut[curvindex]) return kFALSE;
 409         }
 410         else
 411                 return kFALSE;
 412         // calculation of helix matching
 413         Double_t arc1 = 2.0 * r1 * curv;
 414         Double_t arc2 = 2.0 * r2 * curv;
 415         Double_t helmatch = 0.0;
 416         if (arc1 > -1.0 && arc1 < 1.0) arc1 = asin(arc1);
 417         else arc1 = ((arc1 > 0.0) ? 0.5 : 1.5) * TMath::Pi();
 418         if (arc2 > -1.0 && arc2 < 1.0) arc2 = asin(arc2);
 419         else arc2 = ((arc2 > 0.0) ? 0.5 : 1.5) * TMath::Pi();
 420         arc1 /= 2.0 * curv;
 421         arc2 /= 2.0 * curv;
 422         if (arc1 == 0.0 || arc2 == 0.0) return kFALSE;
 423         helmatch = fabs(z1 / arc1 - z2 / arc2);
 424         return (helmatch >= fHelixMatchCutMin[reflayer] && helmatch <= fHelixMatchCutMax[reflayer]);
 425         // END NEW VERSION
 426 }
 427 //
 428 //--------------------------------------------------------------------------------------------
 429 //
 430 Int_t AliITSNeuralTracker::PassAllCuts
 431 (AliITSNode *p1, AliITSNode *p2, Int_t curvindex, Double_t vx, Double_t vy, Double_t vz)
 432 {
 433         // GLOBAL CUT EVALUATOR
 434         //
 435         // Checks all cuts for the passed point pair.
 436         // Return values:
 437         //
 438         // 0 - All cuts passed
 439         // 1 - theta 2D cut not passed
 440         // 2 - theta 3D cut not passed
 441         // 3 - curvature calculated but cut not passed
 442         // 4 - curvature not calculated (division by zero)
 443         // 5 - helix cut not passed
 444         // 6 - curvature inxed out of range
 445
 446         if (curvindex < 0 || curvindex >= fCurvNum) return 6;
 447
 448         // Find the reference layer
 449         Int_t lay1 = p1->GetLayer();
 450         Int_t lay2 = p2->GetLayer();
 451         Int_t reflayer = (lay1 < lay2) ? lay1 : lay2;
 452
 453         // Swap points in order that r1 < r2
 454         AliITSNode *temp = 0;
 455         if (p2->GetLayer() < p1->GetLayer()) {
 456                 temp = p2;
 457                 p2 = p1;
 458                 p1 = temp;
 459         }
 460
 461         // shift XY coords to the reference to the vertex position,
 462         // for easier calculus of quantities.
 463         Double_t x1 = p1->X() - vx;
 464         Double_t x2 = p2->X() - vx;
 465         Double_t y1 = p1->Y() - vy;
 466         Double_t y2 = p2->Y() - vy;
 467         Double_t z1 = p1->Z() - vz;
 468         Double_t z2 = p2->Z() - vz;
 469         Double_t r1 = sqrt(x1*x1 + y1*y1);
 470         Double_t r2 = sqrt(x2*x2 + y2*y2);
 471
 472         // Check for theta cut
 473         Double_t dtheta, dtheta3;
 474         TVector3 v01(z1, r1, 0.0);
 475         TVector3 v12(z2 - z1, r2 - r1, 0.0);
 476         dtheta = v01.Angle(v12) * 180.0 / TMath::Pi();
 477         TVector3 vv01(x1, y1, z1);
 478         TVector3 vv12(x2 - x1, y2 - y1, z2 - z1);
 479         dtheta3 = vv01.Angle(vv12) * 180.0 / TMath::Pi();
 480         if (dtheta < fThetaCut2DMin[reflayer] || dtheta > fThetaCut2DMax[reflayer]) return 1;
 481         if (dtheta3 < fThetaCut3DMin[reflayer] || dtheta3 > fThetaCut3DMax[reflayer]) return 2;
 482
 483         // calculation of curvature
 484         Double_t dx = x1 - x2, dy = y1 - y2;
 485         Double_t num = 2 * (x1*y2 - x2*y1);
 486         Double_t den = r1*r2*sqrt(dx*dx + dy*dy);
 487         Double_t curv = 0.;
 488         if (den != 0.) {
 489                 curv = fabs(num / den);
 490                 if (curv > fCurvCut[curvindex]) return 3;
 491         }
 492         else
 493                 return 4;
 494
 495         // calculation of helix matching
 496         Double_t arc1 = 2.0 * r1 * curv;
 497         Double_t arc2 = 2.0 * r2 * curv;
 498         Double_t helmatch = 0.0;
 499         if (arc1 > -1.0 && arc1 < 1.0) arc1 = asin(arc1);
 500         else arc1 = ((arc1 > 0.0) ? 0.5 : 1.5) * TMath::Pi();
 501         if (arc2 > -1.0 && arc2 < 1.0) arc2 = asin(arc2);
 502         else arc2 = ((arc2 > 0.0) ? 0.5 : 1.5) * TMath::Pi();
 503         arc1 /= 2.0 * curv;
 504         arc2 /= 2.0 * curv;
 505         if (arc1 == 0.0 || arc2 == 0.0) return kFALSE;
 506         helmatch = fabs(z1 / arc1 - z2 / arc2);
 507         if (helmatch < fHelixMatchCutMin[reflayer] || helmatch > fHelixMatchCutMax[reflayer]) return 5;
 508
 509         return 0;
 510 }
 511 //
 512 //--------------------------------------------------------------------------------------------
 513 //
 514 void AliITSNeuralTracker::StoreAbsoluteMatches()
 515 {
 516         // Stores in the 'fMatches' array of each node all the points in the
 517         // adjacent layers which allow to create neurons accordin to the
 518         // helix and theta cut, and also to the largest curvature cut
 519
 520         Int_t ilayer, isector, itheta1, itheta2, check;
 521         TObjArray *list1 = 0, *list2 = 0;
 522         AliITSNode *node1 = 0, *node2 = 0;
 523         Double_t thetamin, thetamax;
 524         Int_t imin, imax;
 525
 526         for (isector = 0; isector < fSectorNum; isector++) {
 527                 // sector is chosen once for both lists
 528                 for (ilayer = 0; ilayer < 5; ilayer++) {
 529                         for (itheta1 = 0; itheta1 < 180; itheta1++) {
 530                                 list1 = (TObjArray*)fPoints[ilayer][itheta1]->At(isector);
 531                                 TObjArrayIter iter1(list1);
 532                                 while ( (node1 = (AliITSNode*)iter1.Next()) ) {
 533                                         if (node1->GetUser() >= 0) continue;
 534                                         node1->fMatches->Clear();
 535                                         thetamin = node1->ThetaDeg() - fPolarInterval;
 536                                         thetamax = node1->ThetaDeg() + fPolarInterval;
 537                                         imin = (Int_t)thetamin;
 538                                         imax = (Int_t)thetamax;
 539                                         if (imin < 0) imin = 0;
 540                                         if (imax > 179) imax = 179;
 541                                         for (itheta2 = imin; itheta2 <= imax; itheta2++) {
 542                                                 list2 = (TObjArray*)fPoints[ilayer + 1][itheta2]->At(isector);
 543                                                 TObjArrayIter iter2(list2);
 544                                                 while ( (node2 = (AliITSNode*)iter2.Next()) ) {
 545                                                         check = PassAllCuts(node1, node2, fCurvNum - 1, fVX, fVY, fVZ);
 546                                                         switch (check) {
 547                                                                 case 0:
 548                                                                         node1->fMatches->AddLast(node2);
 549                                                                         break;
 550                                                                 case 1:
 551                                                                         //Info("StoreAbsoluteMatches", "Layer %d: THETA 2D cut not passed", ilayer);
 552                                                                         break;
 553                                                                 case 2:
 554                                                                         //Info("StoreAbsoluteMatches", "Layer %d: THETA 3D cut not passed", ilayer);
 555                                                                         break;
 556                                                                 case 3:
 557                                                                         //Info("StoreAbsoluteMatches", "Layer %d: CURVATURE cut not passed", ilayer);
 558                                                                         break;
 559                                                                 case 4:
 560                                                                         //Info("StoreAbsoluteMatches", "Layer %d: Division by ZERO in curvature evaluation", ilayer);
 561                                                                         break;
 562                                                                 case 5:
 563                                                                         //Info("StoreAbsoluteMatches", "Layer %d: HELIX-MATCH cut not passed", ilayer);
 564                                                                         break;
 565                                                                 case 6:
 566                                                                         //Error("PassAllCuts", "Layer %d: Curv index out of range", ilayer);
 567                                                                         break;
 568                                                                 default:
 569                                                                         Warning("StoreAbsoluteMatches", "Layer %d: %d: unrecognized return value", ilayer, check);
 570                                                         }
 571                                                 } // while (node2...)
 572                                         } // for (itheta2...)
 573                                 } // while (node1...)
 574                         } // for (itheta...)
 575                 } // for (ilayer...)
 576         } // for (isector...)
 577 }
 578 //
 579 //--------------------------------------------------------------------------------------------
 580 //
 581 void AliITSNeuralTracker::PrintMatches(Bool_t stop)
 582 {
 583         // Prints the matches for each point
 584
 585         TFile *ft = new TFile("its_findables_v1.root");
 586         TTree *tt = (TTree*)ft->Get("Tracks");
 587         Int_t it, nP, nU, lab, nF = (Int_t)tt->GetEntries();
 588         tt->SetBranchAddress("nhitsP", &nP);
 589         tt->SetBranchAddress("nhitsU", &nU);
 590         tt->SetBranchAddress("label", &lab);
 591         TString strP("|"), strU("|");
 592         for (it = 0; it < nF; it++) {
 593                 tt->GetEntry(it);
 594                 if (nP >= 5) strP.Append(Form("%d|", lab));
 595                 if (nU >= 5) strU.Append(Form("%d|", lab));
 596         }
 597
 598         TObjArray *sector = 0;
 599         Int_t ilayer, isector, itheta, id[3];
 600         AliITSNode *node1 = 0, *node2 = 0;
 601
 602         for (ilayer = 0; ilayer < 6; ilayer++) {
 603                 for (isector = 0; isector < fSectorNum; isector++) {
 604                         for (itheta = 0; itheta < 180; itheta++) {
 605                                 sector = (TObjArray*)fPoints[ilayer][itheta]->At(isector);
 606                                 TObjArrayIter points(sector);
 607                                 while ( (node1 = (AliITSNode*)points.Next()) ) {
 608                                         for (it = 0; it < 3; it++) id[it] = node1->GetLabel(it);
 609                                         nF = (Int_t)node1->fMatches->GetSize();
 610                                         cout << "Node layer: " << node1->GetLayer() << ", labels: ";
 611                                         cout << id[0] << " " << id[1] << " " << id[2] << " --> ";
 612                                         if (!nF) {
 613                                                 cout << "NO MatchES!!!" << endl;
 614                                                 continue;
 615                                         }
 616                                         cout << nF << " Matches" << endl;
 617                                         for (it = 0; it < 3; it++) {
 618                                                 if (strP.Contains(Form("|%d|", id[it])))
 619                                                         cout << "Belongs to findable (originary) track " << id[it] << endl;
 620                                                 if (strU.Contains(Form("|%d|", id[it])))
 621                                                         cout << "Belongs to findable (post-Kalman) track " << id[it] << endl;
 622                                         }
 623                                         TObjArrayIter matches(node1->fMatches);
 624                                         while ( (node2 = (AliITSNode*)matches.Next()) ) {
 625                                                 cout << "Match with " << node2;
 626                                                 Int_t *sh = node1->SharedID(node2);
 627                                                 for (Int_t k = 0; k < 3; k++)
 628                                                         if (sh[k] > -1) cout << " " << sh[k];
 629                                                 cout << endl;
 630                         }
 631                                         if (stop) {
 632                                                 cout << "Press a key" << endl;
 633                                                 cin.get();
 634                                         }
 635                                 }
 636                         }
 637                 }
 638         }
 639 }
 640 //
 641 //--------------------------------------------------------------------------------------------
 642 //
 643 void AliITSNeuralTracker::ResetNodes(Int_t isector)
 644 {
 645         // Clears some utility arrays for each node
 646
 647         Int_t ilayer, itheta;
 648         TObjArray *sector = 0;
 649         AliITSNode *node = 0;
 650         for (ilayer = 0; ilayer < 6; ilayer++) {
 651                 for (itheta = 0; itheta < 180; itheta++) {
 652                         sector = (TObjArray*)fPoints[ilayer][itheta]->At(isector);
 653                         TObjArrayIter iter(sector);
 654                         for (;;) {
 655                                 node = (AliITSNode*)iter.Next();
 656                                 if (!node) break;
 657                                 node->fInnerOf->Clear();
 658                                 node->fOuterOf->Clear();
 659                                 delete node->fInnerOf;
 660                                 delete node->fOuterOf;
 661                                 node->fInnerOf = new TObjArray;
 662                                 node->fOuterOf = new TObjArray;
 663                         }
 664                 }
 665         }
 666 }
 667 //
 668 //--------------------------------------------------------------------------------------------
 669 //
 670 void AliITSNeuralTracker::NeuralTracking(const char* filesave, TCanvas*& display)
 671 {
 672 // This is the public method that operates the tracking.
 673 // It works sector by sector, and at the end  saves the found tracks.
 674 // Other methods are privare because they have no meaning id used alone,
 675 // and sometimes they could get segmentation faults due to uninitialized
 676 // datamembert they require to work on.
 677 // The argument is a file where the final results have to be stored.
 678
 679         Bool_t isStable = kFALSE;
 680         Int_t i, nUnits = 0, nLinks = 0, nTracks = 0, sectTracks = 0, totTracks = 0;
 681
 682         // tracking through sectors
 683         cout << endl;
 684         Int_t sector, curv;
 685         for (sector = 0; sector < fSectorNum; sector++) {
 686                 cout << "\rSector " << sector << ": " << endl;
 687                 sectTracks = 0;
 688                 for(curv = 0; curv < fCurvNum; curv++) {
 689                         cout << "- curvature step " << curv + 1;
 690                         cout << " (" << fCurvCut[curv] << "): ";
 691                         // units creation
 692                         nUnits = CreateNeurons(sector, curv);
 693                         if (!nUnits) {
 694                                 cout << "no neurons --> skipped" << endl;
 695                                 continue;
 696                         }
 697                         cout << nUnits << " units, " << flush;
 698                         // units connection
 699                         nLinks = LinkNeurons();
 700                         if (!nLinks) {
 701                                 cout << "no connections --> skipped" << endl;
 702                                 continue;
 703                         }
 704                         cout << nLinks << " connections, " << flush;
 705                         // neural updating
 706                         for (i = 0;; i++) {
 707                                 isStable = Update();
 708                                 if (display) Display(display);
 709                                 TObjArrayIter iter(fNeurons);
 710                                 for (;;) {
 711                                         AliITSneuron *n = (AliITSneuron*)iter.Next();
 712                                         if (!n) break;
 713                                 }
 714                                 if (isStable) break;
 715                         }
 716                         cout << i << " cycles --> " << flush;
 717                         // tracks saving
 718                         CleanNetwork();
 719                         nTracks = Save(sector);
 720                         cout << nTracks << " tracks" << endl;
 721                         sectTracks += nTracks;
 722                 }
 723                 totTracks += sectTracks;
 724                 //cout << sectTracks << " tracks found (total = " << totTracks << ")      " << flush;
 725         }
 726
 727         cout << endl << totTracks << " tracks found!!!" << endl;
 728         cout << endl << "Saving results in file " << filesave << "..." << flush;
 729         TFile *f = new TFile(filesave, "recreate");
 730         fChains->Write("TreeC");
 731         f->Close();
 732 }
 733 //
 734 //--------------------------------------------------------------------------------------------
 735 //
 736 Int_t AliITSNeuralTracker::CreateNeurons(Int_t sectoridx, Int_t curvidx)
 737 {
 738 // Fills the neuron arrays, following the cut criteria for the selected step
 739 // (secnum = sector to analyze, curvnum = curvature cut step to use)
 740 // It also sets the initial random activation.
 741 // In order to avoid a large number of 'new' operations, all existing neurons
 742 // are reset and initialized with the new values, and are created new unit only if
 743 // it turns out to be necessary
 744 // the 'flag' argument is used to decide if the lower edge in the intevral
 745 // of the accepted curvatures is given by zero (kFALSE) or by the preceding used cut (kTRUE)
 746 // (if this is the first step, anyway, the minimum is always zero)
 747
 748         ResetNodes(sectoridx);
 749
 750         if (fNeurons) delete fNeurons;
 751         fNeurons = new TObjArray;
 752
 753         AliITSneuron *unit = 0;
 754         Int_t itheta, neurons = 0;
 755         TObjArray *lstsector = 0;
 756
 757         // NEW VERSION
 758         Double_t vx[6], vy[6], vz[6];
 759         AliITSNode *p[6] = {0, 0, 0, 0, 0, 0};
 760         for (itheta = 0; itheta < 180; itheta++) {
 761                 lstsector = (TObjArray*)fPoints[0][itheta]->At(sectoridx);
 762                 TObjArrayIter lay0(lstsector);
 763                 while ( (p[0] = (AliITSNode*)lay0.Next()) ) {
 764                         if (p[0]->GetUser() >= 0) continue;
 765                         vx[0] = fVX;
 766                         vy[0] = fVY;
 767                         vz[0] = fVZ;
 768                         TObjArrayIter lay1(p[0]->fMatches);
 769                         while ( (p[1] = (AliITSNode*)lay1.Next()) ) {
 770                                 if (p[1]->GetUser() >= 0) continue;
 771                                 if (!PassCurvCut(p[0], p[1], curvidx, fVX, fVY, fVZ)) continue;
 772                                 unit = new AliITSneuron;
 773                                 unit->fInner = p[0];
 774                                 unit->fOuter = p[1];
 775                                 unit->fActivation = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
 776                                 unit->fGain = new TObjArray;
 777                                 fNeurons->AddLast(unit);
 778                                 p[0]->fInnerOf->AddLast(unit);
 779                                 p[1]->fOuterOf->AddLast(unit);
 780                                 neurons++;
 781                                 vx[1] = p[0]->X();
 782                                 vy[1] = p[0]->Y();
 783                                 vz[1] = p[0]->Z();
 784                                 TObjArrayIter lay2(p[1]->fMatches);
 785                                 while ( (p[2] = (AliITSNode*)lay2.Next()) ) {
 786                                         if (p[2]->GetUser() >= 0) continue;
 787                                         if (!PassCurvCut(p[1], p[2], curvidx, vx[1], vy[1], vz[1])) continue;
 788                                         unit = new AliITSneuron;
 789                                         unit->fInner = p[1];
 790                                         unit->fOuter = p[2];
 791                                         unit->fActivation = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
 792                                         unit->fGain = new TObjArray;
 793                                         fNeurons->AddLast(unit);
 794                                         p[1]->fInnerOf->AddLast(unit);
 795                                         p[2]->fOuterOf->AddLast(unit);
 796                                         neurons++;
 797                                         vx[2] = p[1]->X();
 798                                         vy[2] = p[1]->Y();
 799                                         vz[2] = p[1]->Z();
 800                                         TObjArrayIter lay3(p[2]->fMatches);
 801                                         while ( (p[3] = (AliITSNode*)lay3.Next()) ) {
 802                                                 if (p[3]->GetUser() >= 0) continue;
 803                                                 if (!PassCurvCut(p[2], p[3], curvidx, vx[2], vy[2], vz[2])) continue;
 804                                                 unit = new AliITSneuron;
 805                                                 unit->fInner = p[2];
 806                                                 unit->fOuter = p[3];
 807                                                 unit->fActivation = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
 808                                                 unit->fGain = new TObjArray;
 809                                                 fNeurons->AddLast(unit);
 810                                                 p[2]->fInnerOf->AddLast(unit);
 811                                                 p[3]->fOuterOf->AddLast(unit);
 812                                                 neurons++;
 813                                                 vx[3] = p[2]->X();
 814                                                 vy[3] = p[2]->Y();
 815                                                 vz[3] = p[2]->Z();
 816                                                 TObjArrayIter lay4(p[3]->fMatches);
 817                                                 while ( (p[4] = (AliITSNode*)lay4.Next()) ) {
 818                                                         if (p[4]->GetUser() >= 0) continue;
 819                                                         if (!PassCurvCut(p[3], p[4], curvidx, vx[3], vy[3], vz[3])) continue;
 820                                                         unit = new AliITSneuron;
 821                                                         unit->fInner = p[3];
 822                                                         unit->fOuter = p[4];
 823                                                         unit->fActivation = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
 824                                                         unit->fGain = new TObjArray;
 825                                                         fNeurons->AddLast(unit);
 826                                                         p[3]->fInnerOf->AddLast(unit);
 827                                                         p[4]->fOuterOf->AddLast(unit);
 828                                                         neurons++;
 829                                                         vx[4] = p[3]->X();
 830                                                         vy[4] = p[3]->Y();
 831                                                         vz[4] = p[3]->Z();
 832                                                         TObjArrayIter lay5(p[4]->fMatches);
 833                                                         while ( (p[5] = (AliITSNode*)lay5.Next()) ) {
 834                                                                 if (p[5]->GetUser() >= 0) continue;
 835                                                                 if (!PassCurvCut(p[4], p[5], curvidx, vx[4], vy[4], vz[4])) continue;
 836                                                                 unit = new AliITSneuron;
 837                                                                 unit->fInner = p[4];
 838                                                                 unit->fOuter = p[5];
 839                                                                 unit->fActivation = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
 840                                                                 unit->fGain = new TObjArray;
 841                                                                 fNeurons->AddLast(unit);
 842                                                                 p[4]->fInnerOf->AddLast(unit);
 843                                                                 p[5]->fOuterOf->AddLast(unit);
 844                                                                 neurons++;
 845                                                         } // while (p[5])
 846                                                 } // while (p[4])
 847                                         } // while (p[3])
 848                                 } // while (p[2])
 849                         } // while (p[1])
 850                 } // while (p[0])
 851         } // for (itheta...)
 852         // END OF NEW VERSION
 853
 854         /* OLD VERSION
 855         for (ilayer = 0; ilayer < 6; ilayer++) {
 856                 for (itheta = 0; itheta < 180; itheta++) {
 857                         lstsector = (TObjArray*)fPoints[ilayer][itheta]->At(sectoridx);
 858                         TObjArrayIter inners(lstsector);
 859                         while ( (inner = (AliITSNode*)inners.Next()) ) {
 860                                 if (inner->GetUser() >= 0) continue;
 861                                 TObjArrayIter outers(inner->fMatches);
 862                                 while ( (outer = (AliITSNode*)outers.Next()) ) {
 863                                         if (outer->GetUser() >= 0) continue;
 864                                         if (!PassCurvCut(inner, outer, curvidx, fVX, fVY, fVZ)) continue;
 865                                         unit = new AliITSneuron;
 866                                         unit->fInner = inner;
 867                                         unit->fOuter = outer;
 868                                         unit->fActivation = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
 869                                         unit->fGain = new TObjArray;
 870                                         fNeurons->AddLast(unit);
 871                                         inner->fInnerOf->AddLast(unit);
 872                                         outer->fOuterOf->AddLast(unit);
 873                                         neurons++;
 874                                 } // for (;;)
 875                         } // for (;;)
 876                 } // for (itheta...)
 877         } // for (ilayer...)
 878         */
 879
 880         fNeurons->SetOwner();
 881         return neurons;
 882 }
 883 //
 884 //
 885 //
 886 Int_t AliITSNeuralTracker::LinkNeurons() const
 887 {
 888 // Creates the neural synapses among all neurons
 889 // which share a point.
 890
 891         Int_t total = 0;
 892         TObjArrayIter iter(fNeurons), *testiter;
 893         AliITSneuron *neuron = 0, *test = 0;
 894         for (;;) {
 895                 neuron = (AliITSneuron*)iter.Next();
 896                 if (!neuron) break;
 897                 // gain contributors
 898                 testiter = (TObjArrayIter*)neuron->fInner->fOuterOf->MakeIterator();
 899                 for (;;) {
 900                         test = (AliITSneuron*)testiter->Next();
 901                         if (!test) break;
 902                         neuron->Add2Gain(test, fGain2CostRatio, fAlignExponent);
 903                         total++;
 904                 }
 905                 delete testiter;
 906                 testiter = (TObjArrayIter*)neuron->fOuter->fInnerOf->MakeIterator();
 907                 for (;;) {
 908                         test = (AliITSneuron*)testiter->Next();
 909                         if (!test) break;
 910                         neuron->Add2Gain(test, fGain2CostRatio, fAlignExponent);
 911                         total++;
 912                 }
 913                 delete testiter;
 914         }
 915         return total;
 916 }
 917 //
 918 //
 919 //
 920 Bool_t AliITSNeuralTracker::Update()
 921 {
 922 // A complete updating cycle with the asynchronous method
 923 // when the neural network is stable, kTRUE is returned.
 924
 925         Double_t actVar = 0.0, totDiff = 0.0;
 926         TObjArrayIter iter(fNeurons);
 927         AliITSneuron *unit;
 928         for (;;) {
 929                 unit = (AliITSneuron*)iter.Next();
 930                 if (!unit) break;
 931                 actVar = Activate(unit);
 932                 // calculation the relative activation variation
 933                 totDiff += actVar;
 934         }
 935         totDiff /= fNeurons->GetSize();
 936         return (totDiff < fStabThreshold);
 937 }
 938 //
 939 //
 940 //
 941 void AliITSNeuralTracker::CleanNetwork()
 942 {
 943 // Removes all deactivated neurons, and resolves the competitions
 944 // to the advantage of the most active unit
 945
 946         AliITSneuron *unit = 0;
 947         TObjArrayIter neurons(fNeurons);
 948         while ( (unit = (AliITSneuron*)neurons.Next()) ) {
 949                 if (unit->fActivation < fActMinimum) {
 950                         unit->fInner->fInnerOf->Remove(unit);
 951                         unit->fOuter->fOuterOf->Remove(unit);
 952                         delete fNeurons->Remove(unit);
 953                 }
 954         }
 955         return;
 956         Bool_t removed;
 957         Int_t nIn, nOut;
 958         AliITSneuron *enemy = 0;
 959         neurons.Reset();
 960         while ( (unit = (AliITSneuron*)neurons.Next()) ) {
 961                 nIn = (Int_t)unit->fInner->fInnerOf->GetSize();
 962                 nOut = (Int_t)unit->fOuter->fOuterOf->GetSize();
 963                 if (nIn < 2 && nOut < 2) continue;
 964                 removed = kFALSE;
 965                 if (nIn > 1) {
 966                         TObjArrayIter competing(unit->fInner->fInnerOf);
 967                         while ( (enemy = (AliITSneuron*)competing.Next()) ) {
 968                                 if (unit->fActivation > enemy->fActivation) {
 969                                         enemy->fInner->fInnerOf->Remove(enemy);
 970                                         enemy->fOuter->fOuterOf->Remove(enemy);
 971                                         delete fNeurons->Remove(enemy);
 972                                 }
 973                                 else {
 974                                         unit->fInner->fInnerOf->Remove(unit);
 975                                         unit->fOuter->fOuterOf->Remove(unit);
 976                                         delete fNeurons->Remove(unit);
 977                                         removed = kTRUE;
 978                                         break;
 979                                 }
 980                         }
 981                         if (removed) continue;
 982                 }
 983                 if (nOut > 1) {
 984                         TObjArrayIter competing(unit->fOuter->fOuterOf);
 985                         while ( (enemy = (AliITSneuron*)competing.Next()) ) {
 986                                 if (unit->fActivation > enemy->fActivation) {
 987                                         enemy->fInner->fInnerOf->Remove(enemy);
 988                                         enemy->fOuter->fOuterOf->Remove(enemy);
 989                                         delete fNeurons->Remove(enemy);
 990                                 }
 991                                 else {
 992                                         unit->fInner->fInnerOf->Remove(unit);
 993                                         unit->fOuter->fOuterOf->Remove(unit);
 994                                         delete fNeurons->Remove(unit);
 995                                         removed = kTRUE;
 996                                         break;
 997                                 }
 998                         }
 999                 }
1000         }
1001 }
1002 //
1003 //
1004 //
1005 Bool_t AliITSNeuralTracker::CheckOccupation() const
1006 {
1007 // checks if a track recognized has a point for each layer
1008         Int_t i;
1009         for (i = 0; i < 6; i++) if (fPoint[i] < 0) return kFALSE;
1010         return kTRUE;
1011 }
1012 //
1013 //
1014 //
1015 Int_t AliITSNeuralTracker::Save(Int_t sectorid)
1016 {
1017 // Creates chains of active neurons, in order to
1018 // find the tracks obtained as the neural network output.
1019
1020         // every chain is started from the neurons in the first 2 layers
1021         Int_t i, check, stored = 0;
1022         Int_t a = sectorid; a = 0;
1023         Double_t testact = 0;
1024         AliITSneuron *unit = 0, *cursor = 0, *fwd = 0;
1025         AliITSNode *node = 0;
1026         TObjArrayIter iter(fNeurons), *fwditer;
1027         TObjArray *removedUnits = new TObjArray(0);
1028         TObjArray *removedPoints = new TObjArray(6);
1029         removedUnits->SetOwner(kFALSE);
1030         removedPoints->SetOwner(kFALSE);
1031         for (;;) {
1032                 for (i = 0; i < 6; i++) fPoint[i] = -1;
1033                 unit = (AliITSneuron*)iter.Next();
1034                 if (!unit) break;
1035                 if (unit->fInner->GetLayer() > 0) continue;
1036                 fPoint[unit->fInner->GetLayer()] = unit->fInner->fPosInTree;
1037                 fPoint[unit->fOuter->GetLayer()] = unit->fOuter->fPosInTree;
1038                 node = unit->fOuter;
1039                 removedUnits->AddLast(unit);
1040                 removedPoints->AddAt(unit->fInner, unit->fInner->GetLayer());
1041                 removedPoints->AddAt(unit->fOuter, unit->fOuter->GetLayer());
1042                 while (node) {
1043                         testact = fActMinimum;
1044                         fwditer = (TObjArrayIter*)node->fInnerOf->MakeIterator();
1045                         fwd = 0;
1046                         for (;;) {
1047                                 cursor = (AliITSneuron*)fwditer->Next();
1048                                 if (!cursor) break;
1049                                 if (cursor->fUsed) continue;
1050                                 if (cursor->fActivation >= testact) {
1051                                         testact = cursor->fActivation;
1052                                         fwd = cursor;
1053                                 }
1054                         }
1055                         if (!fwd) break;
1056                         removedUnits->AddLast(fwd);
1057                         node = fwd->fOuter;
1058                         fPoint[node->GetLayer()] = node->fPosInTree;
1059                         removedPoints->AddAt(node, node->GetLayer());
1060                 }
1061                 check = 0;
1062                 for (i = 0; i < 6; i++) if (fPoint[i] >= 0) check++;
1063                 if (check >= 6) {
1064                         stored++;
1065                         fChains->Fill();
1066                         for (i = 0; i < 6; i++) {
1067                                 node = (AliITSNode*)removedPoints->At(i);
1068                                 node->SetUser(1);
1069                         }
1070                         fwditer = (TObjArrayIter*)removedUnits->MakeIterator();
1071                         for (;;) {
1072                                 cursor = (AliITSneuron*)fwditer->Next();
1073                                 if(!cursor) break;
1074                                 cursor->fUsed = 1;
1075                         }
1076                 }
1077         }
1078
1079         return stored;
1080 }
1081 /*
1082 Int_t AliITSNeuralTracker::Save(Int_t sectoridx)
1083 // Creates chains of active neurons, in order to
1084 // find the tracks obtained as the neural network output.
1085 {
1086         // Reads the final map of neurons and removes all
1087         // units with too small activation
1088         cout << "saving..." << flush;
1089
1090         // every chain is started from the neurons in the first 2 layers
1091         //                     00111111  00011111  00101111  00110111
1092         Int_t allowmask[] = {     0x3F,     0x1F,     0x2F,     0x37,
1093                                                                                 0x3B,     0x3D,     0x3E };
1094         //                     00111011  00111101  00111110
1095
1096         Double_t test_fwd = 0., testback = 0.;
1097         Int_t ilayer, itheta;
1098         AliITSNode *node = 0;
1099         AliITSneuron *fwd = 0, *back = 0;
1100         TList *listsector = 0;
1101
1102         cout << "A -" << fActMinimum << "-" << flush;
1103         for (ilayer = 0; ilayer < 6; ilayer++) {
1104                 for (itheta = 0; itheta < 180; itheta++) {
1105                         listsector = (TList*)fPoints[ilayer][itheta]->At(sectoridx);
1106                         TListIter iter(listsector);
1107                         while ( (node = (AliITSNode*)iter.Next()) ) {
1108                                 TListIter fwditer(node->fInnerOf);
1109                                 TListIter backiter(node->fOuterOf);
1110                                 testfwd = testback = fActMinimum;
1111                                 while ( (fwd = (AliITSneuron*)fwditer.Next()) ) {
1112                                         if (fwd->fActivation > testfwd) {
1113                                                 testfwd = fwd->fActivation;
1114                                                 node->fNext = fwd->fOuter;
1115                                         }
1116                                 }
1117                                 while ( (back = (AliITSneuron*)backiter.Next()) ) {
1118                                         if (back->fActivation > testback) {
1119                                                 testback = back->fActivation;
1120                                                 node->fPrev = back->fInner;
1121                                         }
1122                                 }
1123                         }
1124                 }
1125         }
1126
1127         cout << "B" << flush;
1128         for (ilayer = 0; ilayer < 5; ilayer++) {
1129                 for (itheta = 0; itheta < 180; itheta++) {
1130                         listsector = (TList*)fPoints[ilayer][itheta]->At(sectoridx);
1131                         TListIter iter(listsector);
1132                         while ( (node = (AliITSNode*)iter.Next()) ) {
1133                                 if (node->fNext) {
1134                                         if (node->fNext->fPrev != node) node->fNext = 0;
1135                                 }
1136                         }
1137                 }
1138         }
1139
1140         cout << "C" << flush;
1141         Bool_t ok;
1142         Int_t stored = 0, l1, l2, i, mask, im;
1143         AliITSNode *temp = 0;
1144         TList *list[2];
1145         for (itheta = 0; itheta < 180; itheta++) {
1146                 list[0] = (TList*)fPoints[0][itheta]->At(sectoridx);
1147                 list[1] = (TList*)fPoints[1][itheta]->At(sectoridx);
1148                 for (i = 0; i < 2; i++) {
1149                         TListIter iter(list[i]);
1150                         while ( (node = (AliITSNode*)iter.Next()) ) {
1151                                 //cout << endl << node << flush;
1152                                 AliITSneuralChain *chain = new AliITSneuralChain;
1153                                 for (temp = node; temp; temp = temp->fNext) {
1154                                         chain->Insert((AliITSNeuralPoint*)temp);
1155                                 }
1156                                 ok = kFALSE;
1157                                 mask = chain->OccupationMask();
1158                                 for (im = 0; im < 7; im++) ok = ok || (mask == allowmask[im]);
1159                                 if (!ok) {
1160                                         delete chain;
1161                                         continue;
1162                                 }
1163                                 //cout << " lab " << flush;
1164                                 chain->AssignLabel();
1165                                 //cout << " add " << flush;
1166                                 fTracks->AddLast(chain);
1167                                 stored++;
1168                                 //cout << " " << stored << flush;
1169                         }
1170                 }
1171         }
1172
1173         cout << "D" << flush;
1174         return stored;
1175 }
1176 */
1177 //
1178 //
1179 //
1180 Double_t AliITSNeuralTracker::Activate(AliITSneuron* &unit)
1181 {
1182 // Computes the new neural activation and
1183 // returns the variation w.r.t the previous one
1184
1185         if (!unit) return 0.0;
1186         Double_t sumgain = 0.0, sumcost = 0.0, input, actOld, actNew;
1187
1188         // sum gain contributions
1189         TObjArrayIter *iter = (TObjArrayIter*)unit->fGain->MakeIterator();
1190         AliITSlink *link;
1191         for(;;) {
1192                 link = (AliITSlink*)iter->Next();
1193                 if (!link) break;
1194                 sumgain += link->fLinked->fActivation * link->fWeight;
1195         }
1196
1197         // sum cost contributions
1198         TObjArrayIter *testiter = (TObjArrayIter*)unit->fInner->fInnerOf->MakeIterator();
1199         AliITSneuron *linked = 0;
1200         for (;;) {
1201                 linked = (AliITSneuron*)testiter->Next();
1202                 if (!linked) break;
1203                 if (linked == unit) continue;
1204                 sumcost += linked->fActivation;
1205         }
1206         delete testiter;
1207         testiter = (TObjArrayIter*)unit->fOuter->fOuterOf->MakeIterator();
1208         for (;;) {
1209                 linked = (AliITSneuron*)testiter->Next();
1210                 if (!linked) break;
1211                 if (linked == unit) continue;
1212                 sumcost += linked->fActivation;
1213         }
1214
1215         //cout << "gain = " << sumgain << ", cost = " << sumcost << endl;
1216
1217         input = (sumgain - sumcost) / fTemperature;
1218         actOld = unit->fActivation;
1219 #ifdef NEURAL_LINEAR
1220         if (input <= -2.0 * fTemperature)
1221                 actNew = 0.0;
1222         else if (input >= 2.0 * fTemperature)
1223                 actNew = 1.0;
1224         else
1225                 actNew = input / (4.0 * fTemperature) + 0.5;
1226 #else
1227         actNew = 1.0 / (1.0 + TMath::Exp(-input));
1228 #endif
1229         unit->fActivation = actNew;
1230         return TMath::Abs(actNew - actOld);
1231 }
1232 //
1233 //
1234 //
1235 void AliITSNeuralTracker::AliITSneuron::Add2Gain(AliITSneuron *n, Double_t multconst, Double_t exponent)
1236 {
1237 // Adds a neuron to the collection of sequenced ones of another neuron
1238
1239         AliITSlink *link = new AliITSlink;
1240         link->fLinked = n;
1241         Double_t weight = Weight(n);
1242         link->fWeight = multconst * TMath::Power(weight, exponent);
1243         fGain->AddLast(link);
1244 }
1245 //
1246 //
1247 //
1248 Double_t AliITSNeuralTracker::AliITSneuron::Weight(AliITSneuron *n)
1249 {
1250 // computes synaptic weight
1251         TVector3 me(fOuter->X() - fInner->X(), fOuter->Y() - fInner->Y(), fOuter->Z() - fInner->Z());
1252         TVector3 it(n->fOuter->X() - n->fInner->X(), n->fOuter->Y() - n->fInner->Y(), n->fOuter->Z() - n->fInner->Z());
1253
1254         Double_t angle = me.Angle(it);
1255         Double_t weight = 1.0 - sin(angle);
1256         return weight;
1257 }