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