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