Changes needed on Sun with Root v4-03-04
[u/mrichter/AliRoot.git] / 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