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