Corrections needed on Sun and HP
[u/mrichter/AliRoot.git] / ITS / AliITStrackerANN.cxx
1 // ---------------------------------------------------------------------------------
2 // AliITStrackerANN
3 // ---------------------------------------------------------------------------------
4 // Neural Network-based algorithm for track reconstruction in the ALICE ITS.
5 // The class contains all that is necessary in order to perform the tracking
6 // using the ITS clusters in the V2 format.
7 // The class is organized with some embedded objects which serve for 
8 // data management, and all setters and getters for all parameters.
9 // The usage of the class from a tracking macro should be based essentially
10 // in the initialization (constructor) and the call of a method which
11 // performs in the right sequence all the operations.
12 // Temporarily (and maybe definitively) all the intermediate steps are 
13 // public functions which could be called independently, but they do not
14 // produce any result if all the preventively required steps have not been 
15 // performed.
16 // ---------------------------------------------------------------------------------
17 // Author: Alberto Pulvirenti (university of Catania)
18 // Email : alberto.pulvirenti@ct.infn.it
19 // ---------------------------------------------------------------------------------
20
21 #include <Riostream.h>
22 #include <TMath.h>
23 #include <TString.h>
24 #include <TObjArray.h>
25 #include <TVector3.h>
26 #include <TFile.h>
27 #include <TTree.h>
28 #include <TRandom.h>
29 #include <TMatrixD.h>
30
31 #include "AliITSgeom.h"
32 #include "AliITStrackSA.h"
33 #include "AliITSclusterV2.h"
34 #include "AliITStrackV2.h"
35
36 #include "AliITStrackerANN.h"
37
38 const Double_t AliITStrackerANN::fgkPi     = 3.141592653; // pi
39 const Double_t AliITStrackerANN::fgkHalfPi = 1.570796327; // pi / 2
40 const Double_t AliITStrackerANN::fgkTwoPi  = 6.283185307; // 2 * pi
41
42 ClassImp(AliITStrackerANN)
43
44 //__________________________________________________________________________________
45 AliITStrackerANN::AliITStrackerANN(const AliITSgeom *geom, Int_t msglev) 
46 : AliITStrackerV2(geom), fMsgLevel(msglev)
47 {
48 /**************************************************************************
49
50                 CONSTRUCTOR (standard-to-use version)
51                 
52                 Arguments:
53                         1) the ITS geometry used in the generated event 
54                         2) the flag for log-messages writing
55                         
56                 The AliITSgeometry is used along the class, 
57                 in order to translate the local AliITSclusterV2 coordinates
58                 into the Global reference frame, which is necessary for the
59                 Neural Network algorithm to operate.
60                 In case of serialized use, the log messages should be excluded, 
61                 in order to save real execution time.
62                 
63                 Operations:
64                         - all pointer data members are initialized 
65                           (if possible, otherwise are set to NULL)
66                         - all numeric data members are initialized to
67                           values which allow the Neural Network to operate
68                           even if nothing is externally set.
69                           
70                  NOTE: it is possible that tracking an event 
71                        with these default values results in a non-sense.
72
73  **************************************************************************/
74
75         Int_t i;
76         
77         // Get ITS geometry
78         fGeom = (AliITSgeom*)geom;
79         
80         // Initialize the array of first module indexes per layer
81         fNLayers = geom->GetNlayers();
82         fFirstModInLayer = new Int_t[fNLayers + 1];
83         for (i = 0; i < fNLayers; i++) {
84                 fFirstModInLayer[i] = fGeom->GetModuleIndex(i + 1, 1, 1);
85         }
86         fFirstModInLayer[fNLayers] = geom->GetIndexMax();
87         
88         // initialization: no curvature cut steps
89         fCurvNum = 0;
90         fCurvCut = 0;
91
92         // initialization: 4 sectors (one for each quadrant)
93         fSectorNum   = 4;
94         fSectorWidth = fgkHalfPi;
95         
96         // initialization: theta offset of 20 degrees
97         fPolarInterval = 20.0;
98
99         // initialization: array structure not defined
100         fStructureOK = kFALSE;
101
102         // initialization: vertex in the origin
103         fVertexX = 0.0;
104         fVertexY = 0.0;
105         fVertexZ = 0.0;
106         
107         // initialization: uninitialized point array
108         fNodes = 0;
109
110         // initialization: very large (unuseful) cut values
111         Int_t ilayer;
112         for (ilayer = 0; ilayer < 6; ilayer++) {
113                 fThetaCut2D[ilayer] = TMath::Pi();
114                 fThetaCut3D[ilayer] = TMath::Pi();
115                 fHelixMatchCut[ilayer] = 1.0;
116         }
117
118         // initialization: inictial activation range between 0.3 and 0.7
119         fEdge1 = 0.3;
120         fEdge2 = 0.7;
121
122         // initialization: neural network operation & weights
123         fTemperature = 1.0;
124         fStabThreshold = 0.001;
125         fGain2CostRatio = 1.0;
126         fExponent = 1.0;
127         fActMinimum = 0.5;
128
129         // initialization: uninitialized array of neurons
130         fNeurons = 0;
131         
132         // initialization: uninitialized array of tracks
133         fFoundTracks = 0;
134 }
135
136 //__________________________________________________________________________________
137 void AliITStrackerANN::SetCuts
138 (Int_t ncurv, Double_t *curv, Double_t *theta2D, Double_t *theta3D, Double_t *helix)
139
140 /**************************************************************************
141  
142         CUT SETTER
143         
144         Allows for the definition of all kind of geometric cuts
145         which have been studied in for the creation of a neuron
146         from a pair of clusters C1 and C2 on consecutive layers.
147         Neuron will be created only if the pair passes ALL of these cuts.
148         
149         At the moment, we define 4 kinds of geometrical cuts:
150                 a) cut on the difference of the polar 'theta' angle;
151                 b) cut on the angle between origin->C1 and C1->C2 in space;
152                 c) cut on the curvature of the circle passing 
153                    through C1, C2 and the primary vertex;
154                 d) cut on heli-matching of the same three points.
155         
156         Arguments:
157                 1) the number of curvature cut steps
158                 2) the array of curvature cuts for each step 
159                    (its dimension is given by the argument 1)
160                 3) array of 5 cut values (one for each consecutive lauer pair)
161                    related to cut a)
162                 4) array of 5 cut values (one for each consecutive lauer pair)
163                    related to cut b)
164                 5) array of 5 cut values (one for each consecutive lauer pair)
165                      related to cut c)
166                           
167         Operations:
168                 - gets the values for each cut and stores them in data members
169                 - in the case of curvature cuts, the cut array 
170                   (whose size is not fixed) is allocated
171           
172         NOTE: in case that the user wants to set onyl ONE of the 4 cuts array,
173               he can simply pass NULL arguments for other cuts, and (eventually)
174                         a ZERO as the first argument (if curvature cuts have not to be set).
175         
176         Anyway, all the cuts have to be set at least once.
177         
178  **************************************************************************/
179 {
180         // counter
181         Int_t i;
182         
183         /*** Curvature cut setting ***/
184         
185         // first of all, the curvature cuts are sorted in increasing order
186         // (from the smallest to the largest one)
187         Int_t *ind = new Int_t[ncurv];
188         TMath::Sort(ncurv, curv, ind, kFALSE);
189         // then, the curvature cut array is allocated and filled
190         // (a message with the list of defined cuts can be optionally shown)
191         fCurvCut = new Double_t[ncurv];
192         if (fMsgLevel >= 1) cout << "Number of curvature cuts: " << ncurv << endl;
193         for (i = 0; i < ncurv; i++) {
194                 fCurvCut[i] = curv[ind[i]];
195                 if (fMsgLevel >= 1) cout << " - " << fCurvCut[i] << endl;
196         }
197         fCurvNum = ncurv;
198         
199         /*** 'Fixed' cuts setting ***/
200         
201         // checks what cuts have to be set
202         Bool_t doTheta2D = (theta2D != 0);
203         Bool_t doTheta3D = (theta3D != 0);
204         Bool_t doHelix = (helix != 0);
205         // sets the cuts for all layer pairs
206         for (i = 0; i < fNLayers - 1; i++) {
207                 if (doTheta2D) fThetaCut2D[i] = theta2D[i];
208                 if (doTheta3D) fThetaCut3D[i] = theta3D[i];
209                 if (doHelix) fHelixMatchCut[i] = helix[i];
210         }
211         // if required, lists the cuts
212         if (!fMsgLevel < 2) return;
213         cout << "Theta 2D cuts: ";
214         if (!doTheta2D) {
215                 cout << "<not set>" << endl;
216         }
217         else {
218                 cout << endl;
219                 for (i = 0; i < fNLayers - 1; i++) {
220                         cout << "For layers " << i << " --> " << i + 1;
221                         cout << " cut = " << fThetaCut2D[i] << endl;
222                 }
223         }
224         cout << "---" << endl;
225         cout << "Theta 3D cuts: ";
226         if (!doTheta3D) {
227                 cout << "<not set>" << endl;
228         }
229         else {
230                 cout << endl;
231                 for (i = 0; i < fNLayers - 1; i++) {
232                         cout << "For layers " << i << " --> " << i + 1;
233                         cout << " cut = " << fThetaCut3D[i] << endl;
234                 }
235         }
236         cout << "---" << endl;
237         cout << "Helix-match cuts: ";
238         if (!doHelix) {
239                 cout << "<not set>" << endl;
240         }
241         else {
242                 cout << endl;
243                 for (i = 0; i < fNLayers - 1; i++) {
244                         cout << "For layers " << i << " --> " << i + 1;
245                         cout << " cut = " << fHelixMatchCut[i] << endl;
246                 }
247         }
248         cout << "---" << endl;
249 }
250
251 //__________________________________________________________________________________ 
252 Bool_t AliITStrackerANN::GetGlobalXYZ
253 (Int_t refIndex, 
254  Double_t &x, Double_t &y, Double_t &z, Double_t &ex2, Double_t &ey2, Double_t &ez2)
255  
256 /**************************************************************************
257
258         LOCAL TO GLOBAL TRANSLATOR
259
260         Taking information from the ITS geometry stored in the class,
261         gets a stored AliITScluster and calculates it coordinates
262         and errors in the global reference frame.
263         These values are stored in the variables, 
264         which are passed by reference.
265                 
266         Arguments:
267                 1) reference index for the cluster to use
268                 2) (by reference) place to store the global X-coord into
269                 3) (by reference) place to store the global Y-coord into
270                 4) (by reference) place to store the global Z-coord into
271                 5) (by reference) place to store the global X-coord error into
272                 6) (by reference) place to store the global Y-coord error into
273                 7) (by reference) place to store the global Z-coord error into
274                 
275         Operations:
276                 essentially, determines the ITS module index from the 
277                 detector index of the AliITSclusterV2 object, and extracts
278                 the roto-translation from the ITS geometry, to convert
279                 the local module coordinates into the global ones.
280                 
281         Return value:
282                 - kFALSE if the given cluster index points to a non-existing
283                 cluster, or if the layer number makes no sense (< 0 or > 6).
284                 - otherwise, kTRUE (meaning a successful operation).
285
286  **************************************************************************/ 
287 {
288         // checks if the layer is correct
289         Int_t ilayer = (refIndex & 0xf0000000) >> 28;
290         if (ilayer < 0 || ilayer >= fNLayers) {
291                 Error("GetGlobalXYZ", "Wrong layer number: %d [range: %d - %d]", ilayer, 0, fNLayers);
292                 return kFALSE;
293         }
294         // checks if the referenced cluster exists and corresponds to the passed reference
295         AliITSclusterV2 *refCluster = (AliITSclusterV2*) GetCluster(refIndex);
296         if (!refCluster) {
297                 Error("GetGlobalXYZ", "Cluster not found for index %d", refIndex);
298                 return kFALSE;
299         }
300         
301         // determine the detector number
302         Int_t detID = refCluster->GetDetectorIndex() + fFirstModInLayer[ilayer];
303         
304         // get rotation matrix
305         Double_t rot[9];
306         fGeom->GetRotMatrix(detID, rot);
307         
308         // get translation vector
309         Float_t tx,ty,tz;
310         fGeom->GetTrans(detID, tx, ty, tz);
311         
312         // determine r and phi for the reference conversion
313         Double_t r = -(Double_t)tx * rot[1] + (Double_t)ty * rot[0];
314         if (ilayer == 0) r = -r;
315         Double_t phi = TMath::ATan2(rot[1],rot[0]);
316         if (ilayer == 0) phi -= fgkPi;
317         
318         // sets values for X, Y, Z in global coordinates and their errors
319         Double_t cosPhi = TMath::Cos(phi);
320         Double_t sinPhi = TMath::Sin(phi);
321         x =  r*cosPhi + refCluster->GetY()*sinPhi;
322         y = -r*sinPhi + refCluster->GetY()*cosPhi;
323         z =  refCluster->GetZ();
324         ex2 = refCluster->GetSigmaY2()*sinPhi*sinPhi;
325         ey2 = refCluster->GetSigmaY2()*cosPhi*cosPhi;
326         ez2 = refCluster->GetSigmaZ2();
327         
328         return kTRUE;
329 }
330
331 //__________________________________________________________________________________ 
332 AliITStrackerANN::AliITSnode* AliITStrackerANN::AddNode(Int_t refIndex)
333
334 /**************************************************************************
335
336         GENERATOR OF NEURAL NODES
337
338         Fills the array of neural 'nodes', which are the ITS clusters
339         translated in the global reference frame.
340         Given that the global coordinates are used many times, they are 
341         stored in a well-defined structure, in the form of an embedded class.
342         Moreover, this class allows a faster navigation among points
343         and neurons, by means of some object arrays, storing only the 
344         neurons which start from, or end to, the given node.
345         Finally, each node contains all the other nodes which match it
346         with respect to the fixed walues, in order to perform a faster
347         neuron-creation phase.
348         
349         Arguments:
350                 1) reference index of the correlated AliITSclusterV2 object
351                 
352         Operations:
353                 - allocates the new AliITSnode objects
354                 - initializes its object arrays
355                 - from the global coordinates, calculates the
356                   'phi' and 'theta' coordinates, in order to store it
357                   into the correct theta-slot and azimutal sector.
358                   
359         REturn values:
360                 - the pointer of the creater AliITSnode object
361                 - in case of errors, a waring is given and a NULL is returned
362
363  **************************************************************************/
364 {
365         // create object and set the reference
366         AliITSnode *node = new AliITSnode;
367         if (!node) {
368                 Warning("AddNode", "Error occurred when allocating AliITSnode");
369                 return 0;
370         }
371         node->ClusterRef() = refIndex;
372         
373         // calls the conversion function, which makes also some checks 
374         // (layer number within range, existence of referenced cluster)
375         if ( !GetGlobalXYZ (
376                                 refIndex, 
377                                 node->X(), node->Y(), node->Z(), 
378                                 node->ErrX2(), node->ErrY2(), node->ErrZ2()
379                         ) ) {return 0;}
380         
381         // initializes the object arrays
382         node->Matches() = new TObjArray;
383         node->InnerOf() = new TObjArray;
384         node->OuterOf() = new TObjArray;
385         
386         // finds azimutal and polar sector (in degrees)
387         Double_t phi = node->GetPhi() * 180.0 / fgkPi;
388         Double_t theta = node->GetTheta() * 180.0 / fgkPi;
389         Int_t isector = (Int_t)(phi / fSectorWidth);
390         Int_t itheta = (Int_t)theta;
391         Int_t ilayer = (refIndex & 0xf0000000) >> 28;
392         
393         // selects the right TObjArray to store object into
394         TObjArray *sector = (TObjArray*)fNodes[ilayer][itheta]->At(isector);
395         sector->AddLast(node);
396         
397         return node;
398 }
399
400 //__________________________________________________________________________________
401 void AliITStrackerANN::CreateArrayStructure(Int_t nsectors)
402
403 /**************************************************************************
404
405         ARRAY STRUCTURE CREATOR
406         
407         Creates a structure of nested TObjArray's where the AliITSnode's
408         have to be stored:
409         - the first level is made by 6 arrays (one for each layer)
410         - the second level is made by 180 arrays (one for each int part of theta)
411         - the third level is made by a variable number of arrays 
412           (one for each azimutal sector)
413           
414         Arguments:
415                 1) the number of azimutal sectors
416                 
417         Operations:
418                 - calculates the width of each sector, from the argument
419                 - allocates and initializes all array levels
420                 - sets a flag which tells the user if this NECESSARY operation
421                   has been performed (it is needed BEFORE performing tracking)
422
423  **************************************************************************/
424 {
425         // Set the number of sectors and their width.
426         fSectorNum   = nsectors;
427         fSectorWidth = 360.0 / (Double_t)fSectorNum;
428         if (fMsgLevel >= 2) {
429                 cout << fSectorNum << " sectors --> sector width (degrees) = " << fSectorWidth << endl;
430         }
431                 
432         // Meaningful indexes
433         Int_t ilayer, isector, itheta;
434         
435         // Mark for the created objects
436         TObjArray *sector = 0;
437         
438         // First index: layer
439         fNodes = new TObjArray**[fNLayers];
440         for (ilayer = 0; ilayer < fNLayers; ilayer++) {
441                 fNodes[ilayer] = new TObjArray*[180];
442                 for (itheta = 0; itheta < 180; itheta++) fNodes[ilayer][itheta] = 0;
443                 for (itheta = 0; itheta < 180; itheta++) {
444                         fNodes[ilayer][itheta] = new TObjArray(nsectors);
445                         for (isector = 0; isector < nsectors; isector++) {
446                                 sector = new TObjArray;
447                                 sector->SetOwner();
448                                 fNodes[ilayer][itheta]->AddAt(sector, isector);
449                         }
450                 }
451         }
452
453         // Sets a checking flag to TRUE. 
454         // This flag is checked before filling up the arrays with the points.
455         fStructureOK = kTRUE;
456 }
457
458 //__________________________________________________________________________________
459 Int_t AliITStrackerANN::ArrangePoints(char *exportFile)
460
461 /**************************************************************************
462
463         POINTS LOCATOR
464         
465         This function assembles the operation from the other above methods, 
466         and fills the arrays with the clusters already stored in the 
467         layers of the tracker.
468         Then, in order to use this method, the user MUSTs call LoadClusters()
469         before.
470         
471         Arguments: 
472                 1) string for a file name where the global ccordinates
473                    of all points can be exported (optional).
474                         If this file must not be created, simply pass a NULL argument
475                         
476         Operations:
477                 - for each AliITSclusterV2 in each AliITSlayer, a ne AliITSnode
478                   is created and stored in the correct location.
479                   
480         Return values:
481                 - the number of stored points
482                 - when errors occur, or no points are found, 0 is returned
483                 
484  **************************************************************************/
485 {
486         // Check if the array structure has been created
487         if (!fStructureOK) {
488                 Error("ArrangePoints", "Structure NOT defined. Call CreateArrayStructure() first");
489                 return 0;
490         }
491
492         // meaningful indexes
493         Int_t ientry, ilayer, nentries = 0, index;
494         Int_t nPtsLayer = 0;
495         
496         // if the argument is not NULL, a file is opened
497         fstream file(exportFile, ios::out);
498         if (!exportFile || file.fail()) {
499                 file.close();
500                 exportFile = 0;
501         }
502         
503         // scan all layers for node creation
504         for (ilayer = 0; ilayer < fNLayers; ilayer++) {
505                 nPtsLayer = GetNumberOfClustersLayer(ilayer);
506                 if (fMsgLevel >= 1) {
507                         cout << "Layer " << ilayer << " --> " << nPtsLayer << " clusters" << endl;
508                 }
509                 for (ientry = 0; ientry < nPtsLayer; ientry++) {
510                         // calculation of cluster index : (Bit mask LLLLIIIIIIIIIIII)
511                         // (L = bits used for layer)
512                         // (I = bits used for position in layer)
513                         index = ilayer << 28;
514                         index += ientry;
515                         // add new AliITSnode object
516                         AliITSnode *n = AddNode(index);
517                         if ( (n != NULL) && exportFile ) {
518                                 file << index << ' ' << n->X() << ' ' << n->Y() << ' ' << n->Z() << endl;
519                         }
520                 }
521                 nentries += nPtsLayer;
522         }
523         
524         // a conventional final message is put at the end of file
525         if (exportFile) {
526                 file << "-1 0.0 0.0 0.0" << endl;
527                 file.close();
528         }
529
530         // returns the number of points processed
531         return nentries;
532 }
533
534 //__________________________________________________________________________________
535 void AliITStrackerANN::StoreOverallMatches()
536
537 /**************************************************************************
538
539         NODE-MATCH ANALYSIS
540         
541         Once the nodes have been created, a firs analysis is to check
542         what pairs will satisfy at least the 'fixed' cuts (theta, helix-match)
543         and the most permissive (= larger) curvature cut.
544         All these node pairs are suitable for neuron creation. 
545         In fact, when performing a Neural Tracking step, the only further check 
546         will be a check against the current curvature step, while the other 
547         are always the same.
548         For thi purpose, each AliITSnode has a data member, named 'fMatches'
549         which contains references to all other AliITSnodes in the successive layer
550         that form, with it, a 'good' pair, with respect to the above cited cuts.
551         Then, in each step for neuron creation, the possible neurons starting from
552         each node will be searched ONLY within the nodes referenced in fMatches.
553         This, of course, speeds up a lot the neuron creation procedure, at the 
554         cost of some memory occupation, which results not to be critical.
555         
556         Operations:
557                 - for each AliITSnode, matches are found according to the criteria
558                   expressed above, and stored in the node->fMatches array
559
560  **************************************************************************/
561 {
562         // meaningful counters
563         Int_t ilayer, isector, itheta1, itheta2, check;
564         TObjArray *list1 = 0, *list2 = 0;
565         AliITSnode *node1 = 0, *node2 = 0;
566         Double_t thetaMin, thetaMax;
567         Int_t imin, imax;
568
569         // Scan for each sector
570         for (isector = 0; isector < fSectorNum; isector++) {
571                 // sector is chosen once for both lists
572                 for (ilayer = 0; ilayer < fNLayers - 1; ilayer++) {
573                         for (itheta1 = 0; itheta1 < 180; itheta1++) {
574                                 list1 = (TObjArray*)fNodes[ilayer][itheta1]->At(isector);
575                                 TObjArrayIter iter1(list1);
576                                 while ( (node1 = (AliITSnode*)iter1.Next()) ) {
577                                         if (node1->IsUsed()) continue;
578                                         // clear an eventually already present array
579                                         // node1->Matches()->Clear();
580                                         // get the global coordinates and defines the theta interval from cut
581                                         thetaMin = (node1->GetTheta() * 180.0 / fgkPi) - fPolarInterval;
582                                         thetaMax = (node1->GetTheta() * 180.0 / fgkPi) + fPolarInterval;
583                                         imin = (Int_t)thetaMin;
584                                         imax = (Int_t)thetaMax;
585                                         if (imin < 0) imin = 0;
586                                         if (imax > 179) imax = 179;
587                                         // loop on the selected theta slots
588                                         for (itheta2 = imin; itheta2 <= imax; itheta2++) {
589                                                 list2 = (TObjArray*)fNodes[ilayer + 1][itheta2]->At(isector);
590                                                 TObjArrayIter iter2(list2);
591                                                 while ( (node2 = (AliITSnode*)iter2.Next()) ) {
592                                                         check = PassAllCuts(node1, node2, fCurvNum - 1, fVertexX, fVertexY, fVertexZ);
593                                                         if (check == 0) {
594                                                                 node1->Matches()->AddLast(node2);
595                                                         }
596                                                 } // while (node2...)
597                                         } // for (itheta2...)
598                                 } // while (node1...)
599                         } // for (itheta...)
600                 } // for (ilayer...)
601         } // for (isector...)
602 }
603
604 //__________________________________________________________________________________
605 Int_t AliITStrackerANN::PassAllCuts
606 (AliITSnode *inner, AliITSnode *outer, Int_t curvStep,                     
607  Double_t vx, Double_t vy, Double_t vz)  
608 {
609 // ***********************************************************************************
610 //
611 //      This check is called in the above method for finding the matches of each node
612 //      It check the passed point pair against all the fixed cuts and a specified
613 //      curvature cut, among all the ones which have been defined.
614 //      The cuts need a vertex-constraint, which is not absolute, but it is passed
615 //      as argument.
616 //      
617 //      Arguments:
618 //              1) the point in the inner layer
619 //              2) the point in the outer layer
620 //              3) curvature step for the curvature cut check (preferably the last)
621 //              4) X of the used vertex
622 //              5) Y of the used vertex
623 //              6) Z of the used vertex
624 //              
625 //      Operations:
626 //              - if necessary, swaps the two points 
627 //                (the first must be in the innermost of the two layers)
628 //              - checks for the theta cuts
629 //              - calculates the circle passing through the vertex
630 //                and the given points and checks for the curvature cut
631 //              - using the radius calculated there, checks for the helix-math cut
632 //      
633 //      Return values:
634 //              0 - All cuts passed
635 //              1 - theta 2D cut not passed
636 //              2 - theta 3D cut not passed
637 //              3 - curvature calculated but cut not passed
638 //              4 - curvature not calculated (division by zero)
639 //              5 - helix cut not passed
640 //              6 - curvature inxed out of range
641 // 
642 // ***********************************************************************************
643         
644         // Check for curvature index
645         if (curvStep < 0 || curvStep >= fCurvNum) return 6;
646         
647         // Swap points in order that r1 < r2
648         AliITSnode *temp = 0;
649         if (outer->GetLayer() < inner->GetLayer()) {
650                 temp = outer;
651                 outer = inner;
652                 inner = temp;
653         }
654         
655         // All cuts are variable according to the layer of the 
656         // innermost point (the other point will surely be 
657         // in the further one, because we don't check poin pairs
658         // which are not in adjacent layers)
659         // The reference is given by the innermost point.
660         Int_t layRef = inner->GetLayer();
661         
662         // The calculations in the transverse plane are made in 
663         // a shifted reference frame, whose origin corresponds to
664         // the reference point passed in the argument.
665         Double_t xIn = inner->X() - vx;
666         Double_t xOut = outer->X() - vx;
667         Double_t yIn = inner->Y() - vy;
668         Double_t yOut = outer->Y() - vy;
669         Double_t zIn = inner->Z() - vz;
670         Double_t zOut = outer->Z() - vz;
671         Double_t rIn = TMath::Sqrt(xIn*xIn + yIn*yIn);
672         Double_t rOut = TMath::Sqrt(xOut*xOut + yOut*yOut);
673         
674         // Check for theta cut.
675         // There are two different angular cuts:
676         // one w.r. to the angle in the 2-dimensional r-z plane...
677         Double_t dthetaRZ;
678         TVector3 origin2innerRZ(zIn, rIn, 0.0);
679         TVector3 inner2outerRZ(zOut - zIn, rOut - rIn, 0.0);
680         dthetaRZ = origin2innerRZ.Angle(inner2outerRZ) * 180.0 / fgkPi;
681         if (dthetaRZ > fThetaCut2D[layRef]) {
682                 return 1;
683                 // r-z theta cut not passed ---> 1
684         }
685         // ...and another w.r. to the angle in the 3-dimensional x-y-z space
686         Double_t dthetaXYZ;
687         TVector3 origin2innerXYZ(xIn, yIn, zIn);
688         TVector3 inner2outerXYZ(xOut - xIn, yOut - yIn, zOut - zIn);
689         dthetaXYZ = origin2innerXYZ.Angle(inner2outerXYZ) * 180.0 / fgkPi;
690         if (dthetaXYZ > fThetaCut3D[layRef]) {
691                 return 2;
692                 // x-y-z theta cut not passed ---> 2
693         }
694         
695         // Calculation & check of curvature
696         Double_t dx = xIn - xOut;
697         Double_t dy = yIn - yOut;
698         Double_t num = 2.0 * (xIn*yOut - xOut*yIn);
699         Double_t den = rIn*rOut*sqrt(dx*dx + dy*dy);
700         Double_t curv = 0.;
701         if (den != 0.) {
702                 curv = TMath::Abs(num / den);
703                 if (curv > fCurvCut[curvStep]) {
704                         return 3;
705                         // curvature too large for cut ---> 3
706                 }
707         }
708         else {
709                 Error("PassAllCuts", "Curvature calculations gives zero denominator");
710                 return 4;
711                 // error: denominator = 0 ---> 4
712         }
713                 
714         // Calculation & check of helix matching
715         Double_t helMatch = 0.0;
716         Double_t arcIn = 2.0 * rIn * curv;
717         Double_t arcOut = 2.0 * rOut * curv;
718         if (arcIn > -1.0 && arcIn < 1.0) 
719                 arcIn = TMath::ASin(arcIn);
720         else 
721                 arcIn = ((arcIn > 0.0) ? 0.5 : 1.5) * TMath::Pi();
722         if (arcOut > -1.0 && arcOut < 1.0) 
723                 arcOut = TMath::ASin(arcOut);
724         else 
725                 arcOut = ((arcOut > 0.0) ? 0.5 : 1.5) * TMath::Pi();
726         arcIn /= 2.0 * curv;
727         arcOut /= 2.0 * curv;
728         if (arcIn == 0.0 || arcOut == 0.0) {
729                 Error("PassAllCuts", "Calculation returns zero-length arcs: l1=%f, l2=%f", arcIn, arcOut);
730                 return 4;
731                 // error: circumference arcs seem to equal zero ---> 4
732         }
733         helMatch = TMath::Abs(zIn / arcIn - zOut / arcOut);
734         if (helMatch > fHelixMatchCut[layRef]) {
735                 return 5;
736                 // helix match cut not passed ---> 5
737         }
738         
739         // ALL CUTS PASSED ---> 0
740         return 0;
741 }
742
743 //__________________________________________________________________________________
744 Bool_t AliITStrackerANN::PassCurvCut
745 (AliITSnode *inner, AliITSnode *outer, 
746  Int_t curvStep, 
747  Double_t vx, Double_t vy, Double_t vz)
748 {
749 //***********************************************************************************
750 //
751 //      This method operates essentially like the above one, but it is used
752 //      during a single step of Neural Tracking, where the curvature cut
753 //      changes.
754 //      Then, not necessaryly all the nodes stored in the fMatches array
755 //      will be suitable for neuron creation in an intermediate step.
756 //      
757 //      It has the same arguments of the PassAllCuts() method, but 
758 //  the theta cut is not checked.
759 //  Moreover, it has a boolean return value.
760 //      
761 //***********************************************************************************
762         
763         // Check for curvature index
764         if (curvStep < 0 || curvStep >= fCurvNum) return 6;
765         
766         // Find the reference layer
767         Int_t layIn = inner->GetLayer();
768         Int_t layOut = outer->GetLayer();
769         Int_t layRef = (layIn < layOut) ? layIn : layOut;
770         
771         // The calculations in the transverse plane are made in 
772         // a shifted reference frame, whose origin corresponds to
773         // the reference point passed in the argument.
774         Double_t xIn = inner->X() - vx;
775         Double_t xOut = outer->X() - vx;
776         Double_t yIn = inner->Y() - vy;
777         Double_t yOut = outer->Y() - vy;
778         Double_t zIn = inner->Z() - vz;
779         Double_t zOut = outer->Z() - vz;
780         Double_t rIn = TMath::Sqrt(xIn*xIn + yIn*yIn);
781         Double_t rOut = TMath::Sqrt(xOut*xOut + yOut*yOut);
782         
783         // Calculation & check of curvature
784         Double_t dx = xIn - xOut;
785         Double_t dy = yIn - yOut;
786         Double_t num = 2.0 * (xIn*yOut - xOut*yIn);
787         Double_t den = rIn*rOut*sqrt(dx*dx + dy*dy);
788         Double_t curv = 0.;
789         /* OLD VERSION
790         if (den != 0.) {
791                 curv = TMath::Abs(num / den);
792                 if (curv > fCurvCut[curvStep]) return kFALSE;
793                 return kTRUE;
794         }
795         else
796                 return kFALSE;
797         */
798         // NEW VERSION
799         if (den != 0.) {
800                 curv = TMath::Abs(num / den);
801                 if (curv > fCurvCut[curvStep]) {
802                         return kFALSE;
803                 }
804         }
805         else {
806                 Error("PassAllCuts", "Curvature calculations gives zero denominator");
807                 return kFALSE;
808         }
809                 
810         // Calculation & check of helix matching
811         Double_t helMatch = 0.0;
812         Double_t arcIn = 2.0 * rIn * curv;
813         Double_t arcOut = 2.0 * rOut * curv;
814         if (arcIn > -1.0 && arcIn < 1.0) 
815                 arcIn = TMath::ASin(arcIn);
816         else 
817                 arcIn = ((arcIn > 0.0) ? 0.5 : 1.5) * TMath::Pi();
818         if (arcOut > -1.0 && arcOut < 1.0) 
819                 arcOut = TMath::ASin(arcOut);
820         else 
821                 arcOut = ((arcOut > 0.0) ? 0.5 : 1.5) * TMath::Pi();
822         arcIn /= 2.0 * curv;
823         arcOut /= 2.0 * curv;
824         if (arcIn == 0.0 || arcOut == 0.0) {
825                 Error("PassAllCuts", "Calculation returns zero-length arcs: l1=%f, l2=%f", arcIn, arcOut);
826                 return 4;
827                 // error: circumference arcs seem to equal zero ---> 4
828         }
829         helMatch = TMath::Abs(zIn / arcIn - zOut / arcOut);
830         return (helMatch <= fHelixMatchCut[layRef]);
831 }
832
833 //__________________________________________________________________________________
834 void AliITStrackerANN::PrintMatches(Bool_t stop)
835 {
836 //      Prints the list of points which appear to match
837 //      each one of them, according to the preliminary 
838 //      overall cuts.
839 //      The arguments states if a pause is required after printing
840 //      the matches for each one. In this case, a keypress is needed.
841
842         TObjArray *sector = 0;
843         Int_t ilayer, isector, itheta, nF;
844         AliITSnode *node1 = 0, *node2 = 0;
845         //AliITSclusterV2 *cluster1 = 0, *cluster2 = 0;
846
847         for (ilayer = 0; ilayer < 6; ilayer++) {
848                 for (isector = 0; isector < fSectorNum; isector++) {
849                         for (itheta = 0; itheta < 180; itheta++) {
850                                 sector = (TObjArray*)fNodes[ilayer][itheta]->At(isector);
851                                 TObjArrayIter points(sector);
852                                 while ( (node1 = (AliITSnode*)points.Next()) ) {
853                                         nF = (Int_t)node1->Matches()->GetEntries();
854                                         cout << "Node layer: " << node1->GetLayer() << " --> ";
855                                         if (!nF) {
856                                                 cout << "NO Matches!!!" << endl;
857                                                 continue;
858                                         }
859                                         cout << nF << " Matches" << endl;
860                                         cout << "Reference cluster: " << hex << node1->ClusterRef() << endl;
861                                         TObjArrayIter matches(node1->Matches());
862                                         while ( (node2 = (AliITSnode*)matches.Next()) ) {
863                                                 cout << "Match with " << hex << node2->ClusterRef() << endl;
864                                         }
865                                         if (stop) {
866                                                 cout << "Press a key" << endl;
867                                                 cin.get();
868                                         }
869                                 }
870                         }
871                 }
872         }
873 }
874
875 //__________________________________________________________________________________
876 void AliITStrackerANN::ResetNodes(Int_t isector)
877 {
878 /***********************************************************************************
879
880         NODE NEURON ARRAY CLEANER
881         
882         After a neural tracking step, this method
883         clears the arrays 'fInnerOf' and 'fOuterOf' of each AliITSnode
884         
885         Arguments:
886                 - the sector where the operation is being executed
887                 
888  ***********************************************************************************/
889
890         Int_t ilayer, itheta;
891         TObjArray *sector = 0;
892         AliITSnode *node = 0;
893         for (ilayer = 0; ilayer < fNLayers; ilayer++) {
894                 for (itheta = 0; itheta < 180; itheta++) {
895                         sector = (TObjArray*)fNodes[ilayer][itheta]->At(isector);
896                         TObjArrayIter iter(sector);
897                         for (;;) {
898                                 node = (AliITSnode*)iter.Next();
899                                 if (!node) break;
900                                 node->InnerOf()->Clear();
901                                 node->OuterOf()->Clear();
902                                 /*
903                                 delete node->InnerOf();
904                                 delete node->OuterOf();
905                                 node->InnerOf() = new TObjArray;
906                                 node->OuterOf() = new TObjArray;
907                                 */
908                         }
909                 }
910         }
911 }
912
913 //__________________________________________________________________________________
914 Int_t AliITStrackerANN::CreateNeurons
915 (AliITSnode *node, Int_t curvStep, Double_t vx, Double_t vy, Double_t vz)
916 {
917         // This method is used to create alle suitable neurons starting from
918         // a single AliITSnode. Each unit is also stored in the fInnerOf array
919         // of the passed node, and in the fOuterOf array of the other neuron edge.
920         // In the new implementation of the intermediate check steps, a further one
921         // is made, which chechs how well a helix is matched by three points
922         // in three consecutive layers.
923         // Then, a vertex-constrained check is made with vertex located
924         // in a layer L, for points in layers L+1 and L+2.
925         //      
926         // In order to do this, the creator works recursively, in a tree-visit like operation.
927         // The neurons are effectively created only if the node argument passed is in 
928         // the 5th layer (they are created between point of 5th and 6th layer).
929         // If the node is in an inner layer, its coordinates are passet as vertex for a nested
930         // call of the same function in the next two layers.
931         // 
932         // Arguments:
933         //      1) reference node
934         //      2) current curvature cut step
935         //      3) X of referenced temporary (not primary) vertex
936         //      4) Y of referenced temporary (not primary) vertex
937         //      5) Z of referenced temporary (not primary) vertex
938         //      
939         // Operations:
940         //      - if the layer is the 5th, neurons are created with nodes
941         //        in the fMatches array of the passed node
942         //      - otherwise, the X, Y, Z of the passed node are given as 
943         //        vertex and the same procedure is recursively called for all
944         //        nodes in the fMatches array of the passed one.
945         //        
946         // Return values:
947         //      - the total number of neurons created from the passed one
948         //        summed with all neurons created from all nodes well matched with it
949         //        (assumes a meaning only for nodes in the first layer)
950         
951         // local variables
952         Int_t made = 0;         // counter
953         Bool_t found = 0;       // flag
954         AliITSnode *match = 0;  // cursor for a AliITSnode array
955         AliITSneuron *unit = 0; // cursor for a AliITSneuron array
956         
957         // --> Case 0: the passed node has already been used 
958         // as member of a track found in a previous step. 
959         // In this case, of course, the function exits.
960         if (node->IsUsed()) return 0;
961         
962         // --> Case 1: there are NO well-matched points.
963         // This can happen in all ITS layers, but it happens **for sure**
964         // for a node in the 'last' layer.
965         // Even in this case, the function exits.
966         if (node->Matches()->IsEmpty()) return 0;
967         
968         // --> Case 2: there are well-matched points.
969         // In this case, the function creates a neuron for each
970         // well-matched pair (according to the cuts for the current step)
971         // Moreover, before storing the neuron, a check is necessary
972         // to avoid the duplicate creation of the same neuron twice.
973         // (This could happen if the 3 last arguments of the function
974         // are close enough to cause a good match for the current step
975         // between two points, independently of their difference).
976         // Finally, a node is skipped if it has already been used.
977         // For each matched point for which a neuron is created, the procedure is 
978         // recursively called.
979         TObjArrayIter matches(node->Matches());
980         while ( (match = (AliITSnode*)matches.Next()) ) {
981                 if (match->IsUsed()) continue;
982                 if (!PassCurvCut(node, match, curvStep, vx, vy, vz)) continue;
983                 found = kFALSE;
984                 if (!node->InnerOf()->IsEmpty()) {
985                         TObjArrayIter presentUnits(node->InnerOf());
986                         while ( (unit = (AliITSneuron*)presentUnits.Next()) ) {
987                                 if (unit->Inner() == node && unit->Outer() == match) {
988                                         found = kTRUE;
989                                         break;
990                                 }
991                         }
992                 }
993                 if (found) continue;
994                 AliITSneuron *unit = new AliITSneuron(node, match, fEdge2, fEdge1);
995                 fNeurons->AddLast(unit);
996                 node->InnerOf()->AddLast(unit);
997                 match->OuterOf()->AddLast(unit);
998                 made += CreateNeurons(match, curvStep, node->X(), node->Y(), node->Z());
999                 made++;
1000         }
1001         
1002         // Of course, the return value contains the number of neurons
1003         // counting in also the oned created in all levels of recursive calls.
1004         return made;
1005 }
1006  
1007 //__________________________________________________________________________________
1008 Int_t AliITStrackerANN::CreateNetwork(Int_t sector, Int_t curvStep)
1009 {
1010         // This function simply recalls the CreateNeurons() method for each node
1011         // in the first layer, for the current sector. 
1012         // This generates the whole network, thanks to the recursive calls.
1013         //
1014         // Arguments:
1015         //       1) current sector
1016         //       2) current curvature step
1017         //      
1018         // Operations:
1019         //      - scans the nodes array for all theta's in the current sector
1020         //        and layer 0, and calls the CreateNeurons() function.
1021         
1022         // removes all eventually present neurons
1023         if (fNeurons) delete fNeurons;
1024         fNeurons = new TObjArray;
1025         fNeurons->SetOwner(kTRUE);
1026         
1027         // calls the ResetNodes() function to free the AliITSnode arrays
1028         if (fMsgLevel >= 2) {
1029                 cout << "Sector " << sector << " PHI = ";
1030                 cout << fSectorWidth * (Double_t)sector << " --> ";
1031                 cout << fSectorWidth * (Double_t)(sector + 1) << endl;
1032                 cout << "Curvature step " << curvStep << " [cut = " << fCurvCut[curvStep] << "]" << endl;
1033         }
1034         ResetNodes(sector);
1035
1036         // meaningful counters
1037         Int_t itheta, neurons = 0;
1038         TObjArray *lstSector = 0;
1039         
1040         // NEW VERSION
1041         Double_t vx[6], vy[6], vz[6];
1042         AliITSnode *p[6] = {0, 0, 0, 0, 0, 0};
1043         for (itheta = 0; itheta < 180; itheta++) {
1044                 lstSector = (TObjArray*)fNodes[0][itheta]->At(sector);
1045                 TObjArrayIter lay0(lstSector);
1046                 while ( (p[0] = (AliITSnode*)lay0.Next()) ) {
1047                         if (p[0]->IsUsed()) continue;
1048                         vx[0] = fVertexX;
1049                         vy[0] = fVertexY;
1050                         vz[0] = fVertexZ;
1051                         neurons += CreateNeurons(p[0], curvStep, fVertexX, fVertexY, fVertexZ);
1052                         /*
1053                         TObjArrayIter lay1(p[0]->Matches());
1054                         while ( (p[1] = (AliITSnode*)lay1.Next()) ) {
1055                                 if (p[1]->IsUsed()) continue;
1056                                 if (!PassCurvCut(p[0], p[1], curvStep, vx[0], vy[0], vz[0])) continue;
1057                                 unit = new AliITSneuron;
1058                                 unit->Inner() = p[0];
1059                                 unit->Outer() = p[1];
1060                                 unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1061                                 unit->fGain = new TObjArray;
1062                                 fNeurons->AddLast(unit);
1063                                 p[0]->InnerOf()->AddLast(unit);
1064                                 p[1]->OuterOf()->AddLast(unit);
1065                                 neurons++;
1066                                 vx[1] = p[0]->X();
1067                                 vy[1] = p[0]->Y();
1068                                 vz[1] = p[0]->Z();
1069                                 TObjArrayIter lay2(p[1]->Matches());
1070                                 while ( (p[2] = (AliITSnode*)lay2.Next()) ) {
1071                                         if (p[2]->IsUsed()) continue;
1072                                         if (!PassCurvCut(p[1], p[2], curvStep, vx[1], vy[1], vz[1])) continue;
1073                                         unit = new AliITSneuron;
1074                                         unit->Inner() = p[1];
1075                                         unit->Outer() = p[2];
1076                                         unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1077                                         unit->fGain = new TObjArray;
1078                                         fNeurons->AddLast(unit);
1079                                         p[1]->InnerOf()->AddLast(unit);
1080                                         p[2]->OuterOf()->AddLast(unit);
1081                                         neurons++;
1082                                         vx[2] = p[1]->X();
1083                                         vy[2] = p[1]->Y();
1084                                         vz[2] = p[1]->Z();
1085                                         TObjArrayIter lay3(p[2]->Matches());
1086                                         while ( (p[3] = (AliITSnode*)lay3.Next()) ) {
1087                                                 if (p[3]->IsUsed()) continue;
1088                                                 if (!PassCurvCut(p[2], p[3], curvStep, vx[2], vy[2], vz[2])) continue;
1089                                                 unit = new AliITSneuron;
1090                                                 unit->Inner() = p[2];
1091                                                 unit->Outer() = p[3];
1092                                                 unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1093                                                 unit->fGain = new TObjArray;
1094                                                 fNeurons->AddLast(unit);
1095                                                 p[2]->InnerOf()->AddLast(unit);
1096                                                 p[3]->OuterOf()->AddLast(unit);
1097                                                 neurons++;
1098                                                 vx[3] = p[2]->X();
1099                                                 vy[3] = p[2]->Y();
1100                                                 vz[3] = p[2]->Z();
1101                                                 TObjArrayIter lay4(p[3]->Matches());
1102                                                 while ( (p[4] = (AliITSnode*)lay4.Next()) ) {
1103                                                         if (p[4]->IsUsed()) continue;
1104                                                         if (!PassCurvCut(p[3], p[4], curvStep, vx[3], vy[3], vz[3])) continue;
1105                                                         unit = new AliITSneuron;
1106                                                         unit->Inner() = p[3];
1107                                                         unit->Outer() = p[4];
1108                                                         unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1109                                                         unit->fGain = new TObjArray;
1110                                                         fNeurons->AddLast(unit);
1111                                                         p[3]->InnerOf()->AddLast(unit);
1112                                                         p[4]->OuterOf()->AddLast(unit);
1113                                                         neurons++;
1114                                                         vx[4] = p[3]->X();
1115                                                         vy[4] = p[3]->Y();
1116                                                         vz[4] = p[3]->Z();
1117                                                         TObjArrayIter lay5(p[4]->Matches());
1118                                                         while ( (p[5] = (AliITSnode*)lay5.Next()) ) {
1119                                                                 if (p[5]->IsUsed()) continue;
1120                                                                 if (!PassCurvCut(p[4], p[5], curvStep, vx[4], vy[4], vz[4])) continue;
1121                                                                 unit = new AliITSneuron;
1122                                                                 unit->Inner() = p[4];
1123                                                                 unit->Outer() = p[5];
1124                                                                 unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1125                                                                 unit->fGain = new TObjArray;
1126                                                                 fNeurons->AddLast(unit);
1127                                                                 p[4]->InnerOf()->AddLast(unit);
1128                                                                 p[5]->OuterOf()->AddLast(unit);
1129                                                                 neurons++;
1130                                                         } // while (p[5])
1131                                                 } // while (p[4])
1132                                         } // while (p[3])
1133                                 } // while (p[2])
1134                         } // while (p[1])
1135                         */
1136                 } // while (p[0])
1137         } // for (itheta...)
1138         // END OF NEW VERSION
1139
1140         /* OLD VERSION
1141         for (ilayer = 0; ilayer < 6; ilayer++) {
1142                 for (itheta = 0; itheta < 180; itheta++) {
1143                         lstSector = (TObjArray*)fNodes[ilayer][itheta]->At(sector_idx);
1144                         TObjArrayIter inners(lstSector);
1145                         while ( (inner = (AliITSnode*)inners.Next()) ) {
1146                                 if (inner->GetUser() >= 0) continue;
1147                                 TObjArrayIter outers(inner->Matches());
1148                                 while ( (outer = (AliITSnode*)outers.Next()) ) {
1149                                         if (outer->GetUser() >= 0) continue;
1150                                         if (!PassCurvCut(inner, outer, curvStep, fVX, fVY, fVZ)) continue;
1151                                         unit = new AliITSneuron;
1152                                         unit->Inner() = inner;
1153                                         unit->Outer() = outer;
1154                                         unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1155                                         unit->fGain = new TObjArray;
1156                                         fNeurons->AddLast(unit);
1157                                         inner->InnerOf()->AddLast(unit);
1158                                         outer->OuterOf()->AddLast(unit);
1159                                         neurons++;
1160                                 } // for (;;)
1161                         } // for (;;)
1162                 } // for (itheta...)
1163         } // for (ilayer...)
1164         */
1165         
1166         fNeurons->SetOwner();
1167         return neurons;
1168  }
1169  
1170  //__________________________________________________________________________________
1171 Int_t AliITStrackerANN::LinkNeurons()
1172 {
1173 /***********************************************************************************
1174         
1175         SYNAPSIS GENERATOR
1176         
1177         Scans the whole neuron array, in order to find all neuron pairs
1178         which are connected in sequence and share a positive weight.
1179         For each of them, an AliITSlink is created, which stores 
1180         the weight value, and will allow for a faster calculation
1181         of the total neural input for each updating cycle.
1182         
1183         Every neuron contains an object array which stores all AliITSlink
1184         objects which point to sequenced units, with the respective weights.
1185         
1186         Return value:
1187                 - the number of link created for the neural network. 
1188                   If they are 0, no updating can be done and the step is skipped.
1189         
1190  ***********************************************************************************/
1191
1192         // meaningful indexes
1193         Int_t total = 0;
1194         Double_t weight = 0.0;
1195         TObjArrayIter neurons(fNeurons), *iter;
1196         AliITSneuron *neuron = 0, *test = 0;
1197         
1198         // scan in the neuron array
1199         for (;;) {
1200                 neuron = (AliITSneuron*)neurons.Next();
1201                 if (!neuron) break;
1202                 // checks for neurons startin from its head ( -> )
1203                 iter = (TObjArrayIter*)neuron->Inner()->OuterOf()->MakeIterator();
1204                 for (;;) {
1205                         test = (AliITSneuron*)iter->Next();
1206                         if (!test) break;
1207                         weight = Weight(test, neuron);
1208                         if (weight > 0.0) neuron->Gain()->AddLast(new AliITSlink(weight, test));
1209                         total++;
1210                 }
1211                 delete iter;
1212                 // checks for neurons ending in its tail ( >- )
1213                 iter = (TObjArrayIter*)neuron->Outer()->InnerOf()->MakeIterator();
1214                 for (;;) {
1215                         test = (AliITSneuron*)iter->Next();
1216                         if (!test) break;
1217                         weight = Weight(neuron, test);
1218                         if (weight > 0.0) neuron->Gain()->AddLast(new AliITSlink(weight, test));
1219                         total++;
1220                 }
1221                 delete iter;
1222         }
1223         return total;
1224 }
1225  
1226 //__________________________________________________________________________________
1227 Bool_t AliITStrackerANN::Update()
1228 {
1229 /***********************************************************************************
1230         
1231         Performs a single updating cycle.
1232         
1233         Operations:
1234                 - for each neuron, gets the activation with the neuron Activate() method
1235                 - checks if stability has been reached (compare mean activation variation
1236                   with the stability threshold data member)
1237         
1238         Return values:
1239                 - kTRUE means that the neural network has stabilized
1240                 - kFALSE means that another updating cycle is needed
1241                 
1242  ***********************************************************************************/
1243
1244         Double_t actVar = 0.0, totDiff = 0.0;
1245         TObjArrayIter iter(fNeurons);
1246         AliITSneuron *unit;
1247         for (;;) {
1248                 unit = (AliITSneuron*)iter.Next();
1249                 if (!unit) break;
1250                 actVar = unit->Activate(fTemperature);
1251                 // calculation the relative activation variation
1252                 totDiff += actVar;
1253         }
1254         totDiff /= fNeurons->GetSize();
1255         return (totDiff < fStabThreshold);
1256 }
1257
1258 //__________________________________________________________________________________
1259 void AliITStrackerANN::FollowChains(Int_t sector)
1260 {
1261 /***********************************************************************************
1262         
1263         CHAINS CREATION
1264         
1265         After that the neural network has stabilized, 
1266         the final step is to create polygonal chains 
1267         of clusters, one in each layer, which represent 
1268         the tracks recognized by the neural algorithm.
1269         This is made by means of a choice of the best 
1270         neuron among the ones starting from each point.
1271         
1272         Once that such neuron is selected, its inner point
1273         will set the 'fNext' field to its outer point, and 
1274         similarly, its outer point will set the 'fPrev' field
1275         to its inner point. 
1276         This defines a bi-directional sequence. 
1277         
1278         In this procedure, it can happen that many neurons
1279         which have the head of the arrow in a given node, will
1280         all select as best following the neuron with the largest 
1281         activation starting in that point. 
1282         This results in MANY nodes which have the same 'fNext'.
1283         But, this field will be set to NULL for all these points, 
1284         but the only one which is pointed by the 'fPrev' field
1285         of this shared node.
1286  
1287  ***********************************************************************************/
1288
1289         // meaningful counters
1290         Int_t itheta, ilayer;
1291         TObjArray *lstSector = 0;
1292         Double_t test = fActMinimum;
1293         AliITSnode *p = 0;
1294         AliITSneuron *n = 0;
1295         
1296         // scan the whole collection of nodes
1297         for (ilayer = 0; ilayer < fNLayers; ilayer++) {
1298                 for (itheta = 0; itheta < 180; itheta++) {
1299                         // get the array of point in a given layer/theta-slot/sector
1300                         lstSector = (TObjArray*)fNodes[ilayer][itheta]->At(sector);
1301                         TObjArrayIter nodes(lstSector);
1302                         while ( (p = (AliITSnode*)nodes.Next()) ) {
1303                                 // if the point is used, it is skipped
1304                                 if (p->IsUsed()) continue;
1305                                 // initially, fNext points to nothing, and
1306                                 // the comparison value is set to the minimum activation
1307                                 // which allows to say that a neuron is turned 'on'
1308                                 // a node from which only 'off' neurons start is probably
1309                                 // a noise point, which will be excluded from the reconstruction.
1310                                 test = fActMinimum;
1311                                 p->Next() = 0;
1312                                 TObjArrayIter innerof(p->InnerOf());
1313                                 while ( (n = (AliITSneuron*)innerof.Next()) ) {
1314                                         // if the examined neuron has not the largest activation
1315                                         // it is skipped and removed from array of all neurons
1316                                         // and of its outer point (its inner is the cursor p)
1317                                         if (n->Activation() < test) {
1318                                                 p->InnerOf()->Remove(n);
1319                                                 n->Outer()->OuterOf()->Remove(n);
1320                                                 delete fNeurons->Remove(n);
1321                                                 continue;
1322                                         }
1323                                         // otherwise, its activation becomes the maximum reference
1324                                         p->Next() = n->Outer();
1325                                         // at the exit of the while(), the fNext will point
1326                                         // to the outer node of the neuron starting in p, whose
1327                                         // activation is the largest.
1328                                 }
1329                                 // the same procedure is made now for all neurons
1330                                 // for which p is the outer point
1331                                 test = fActMinimum;
1332                                 p->Prev() = 0;
1333                                 TObjArrayIter outerof(p->OuterOf());
1334                                 while ( (n = (AliITSneuron*)outerof.Next()) ) {
1335                                         // if the examined neuron has not the largest activation
1336                                         // it is skipped and removed from array of all neurons
1337                                         // and of its inner point (its outer is the cursor p)
1338                                         if (n->Activation() < test) {
1339                                                 p->OuterOf()->Remove(n);
1340                                                 n->Inner()->InnerOf()->Remove(n);
1341                                                 delete fNeurons->Remove(n);
1342                                                 continue;
1343                                         }
1344                                         // otherwise, its activation becomes the maximum reference
1345                                         p->Prev() = n->Inner();
1346                                         // at the exit of the while(), the fPrev will point
1347                                         // to the inner node of the neuron ending in p, whose
1348                                         // activation is the largest.
1349                                 }
1350                         } // end while (p ...)
1351                 } // end for (itheta ...)
1352         } // end for (ilayer ...)
1353         
1354         // now the mismatches are solved
1355         Bool_t matchPrev, matchNext;
1356         for (ilayer = 0; ilayer < fNLayers; ilayer++) {
1357                 for (itheta = 0; itheta < 180; itheta++) {
1358                         // get the array of point in a given layer/theta-slot/sector
1359                         lstSector = (TObjArray*)fNodes[ilayer][itheta]->At(sector);
1360                         TObjArrayIter nodes(lstSector);
1361                         while ( (p = (AliITSnode*)nodes.Next()) ) {
1362                                 // now p will point to a fPrev and a fNext node.
1363                                 // Ideally they are placed this way: fPrev --> P --> fNext
1364                                 // A mismatch happens if the point addressed as fPrev does NOT
1365                                 // point to p as its fNext. And the same for the point addressed
1366                                 // as fNext.
1367                                 // In this case, the fNext and fPrev pointers are set to NULL
1368                                 // and p is excluded from the reconstruction
1369                                 matchPrev = matchNext= kFALSE;
1370                                 if (ilayer > 0 && p->Prev() != NULL) 
1371                                         if (p->Prev()->Next() == p) matchPrev = kTRUE;
1372                                 if (ilayer < 5 && p->Next() != NULL) 
1373                                         if (p->Next()->Prev() == p) matchNext = kTRUE;
1374                                 if (ilayer == 0) 
1375                                         matchPrev = kTRUE;
1376                                 else if (ilayer == 5)
1377                                         matchNext = kTRUE;
1378                                 if (!matchNext || !matchPrev) {
1379                                         p->Prev() = p->Next() = 0;
1380                                 }
1381                         } // end while (p ...)
1382                 } // end for (itheta ...)
1383         } // end for (ilayer ...)
1384 }
1385
1386 //__________________________________________________________________________________
1387 Int_t AliITStrackerANN::SaveTracks(Int_t sector)
1388 {
1389 /********************************************************************************
1390
1391         TRACK SAVING
1392         ------------    
1393         Using the fNext and fPrev pointers, the chain is followed 
1394         and the track is fitted and saved.
1395         Of course, the track is followed as a chain with a point
1396         for each layer, then the track following starts always
1397         from the clusters in layer 0.
1398         
1399 ***********************************************************************************/
1400  
1401         // if not initialized, the tracks TobjArray is initialized
1402         if (!fFoundTracks) fFoundTracks = new TObjArray;
1403
1404         // meaningful counters
1405         Int_t itheta, ilayer, l;
1406         TObjArray *lstSector = 0;
1407         AliITSnode *p = 0, *q = 0, **node = new AliITSnode*[fNLayers];
1408         for (l = 0; l < fNLayers; l++) node[l] = 0;
1409         
1410         /*
1411         array = new AliITSnode*[fNLayers + 1];
1412         for (l = 0; l <= fNLayers; l++) array[l] = 0;
1413         array[0] = new AliITSnode();
1414         array[0]->X() = fVertexX;
1415         array[0]->Y() = fVertexY;
1416         array[0]->Z() = fVertexZ;
1417         array[0]->ErrX2() = fVertexErrorX;
1418         array[0]->ErrY2() = fVertexErrorY;
1419         array[0]->ErrZ2() = fVertexErrorZ;
1420         */
1421         Double_t *param = new Double_t[8];
1422         
1423         // scan the whole collection of nodes
1424         for (ilayer = 0; ilayer < 1; ilayer++) {
1425                 for (itheta = 0; itheta < 180; itheta++) {
1426                         // get the array of point in a given layer/theta-slot/sector
1427                         lstSector = (TObjArray*)fNodes[ilayer][itheta]->At(sector);
1428                         TObjArrayIter nodes(lstSector);
1429                         while ( (p = (AliITSnode*)nodes.Next()) ) {
1430                                 for (q = p; q; q = q->Next()) {
1431                                         l = q->GetLayer();
1432                                         node[l] = q;
1433                                 }
1434                                 //if (!RiemannFit(fNLayers, node, param)) continue;
1435                                 // initialization of Kalman Filter Tracking
1436                                 AliITSclusterV2 *cluster = (AliITSclusterV2*)GetCluster(node[0]->ClusterRef());
1437                                 Int_t mod = cluster->GetDetectorIndex();
1438                                 Int_t lay, lad, det;
1439                                 fGeom->GetModuleId(mod, lay, lad, det);
1440                                 Float_t y0 = cluster->GetY();
1441                                 Float_t z0 = cluster->GetZ();
1442                                 AliITStrackSA* trac = new AliITStrackSA(lay, lad, det, 
1443                                                                         y0, z0, 
1444                                                                                                                                          param[4], param[7], param[3], 1);
1445                                 for (l = 0; l < fNLayers; l++) {
1446                                         cluster = (AliITSclusterV2*)GetCluster(node[l]->ClusterRef());
1447                                         if (cluster) trac->AddClusterV2(l, (node[l]->ClusterRef() & 0x0fffffff)>>0);
1448                                 }
1449             AliITStrackV2* ot = new AliITStrackV2(*trac);
1450                                 ot->ResetCovariance();
1451                                 ot->ResetClusters();
1452                                 if (RefitAt(49.,ot,trac)) { //fit from layer 1 to layer 6
1453                                         AliITStrackV2 *otrack2 = new AliITStrackV2(*ot);
1454                                         otrack2->ResetCovariance();
1455                                         otrack2->ResetClusters();
1456                                         //fit from layer 6 to layer 1
1457                                         if (RefitAt(3.7,otrack2,ot)) fFoundTracks->AddLast(otrack2);
1458                                 }       
1459                                 // end of Kalman Filter fit
1460                         }
1461                 }
1462         }
1463         
1464         return 1;
1465 }
1466  
1467 //__________________________________________________________________________________
1468 void AliITStrackerANN::ExportTracks(const char *filename) const
1469 {
1470 // Exports found tracks into a TTree of AliITStrackV2 objects
1471         TFile *file = new TFile(filename, "RECREATE");
1472         TTree *tree = new TTree("TreeT-ANN", "Tracks found in ITS stand-alone with Neural Tracking");
1473         AliITStrackV2 *track = 0;
1474         tree->Branch("Tracks", &track, "AliITStrackV2");
1475         TObjArrayIter tracks(fFoundTracks);
1476         while ( (track = (AliITStrackV2*)tracks.Next()) ) {
1477                 tree->Fill();
1478         }
1479         file->cd();
1480         tree->Write();
1481         file->Close();
1482 }
1483
1484
1485 //__________________________________________________________________________________
1486 void AliITStrackerANN::CleanNetwork()
1487 {
1488         // Removes deactivated units from the network
1489
1490         AliITSneuron *unit = 0;
1491         TObjArrayIter neurons(fNeurons);
1492         while ( (unit = (AliITSneuron*)neurons.Next()) ) {
1493                 if (unit->Activation() < fActMinimum) {
1494                         unit->Inner()->InnerOf()->Remove(unit);
1495                         unit->Outer()->OuterOf()->Remove(unit);
1496                         delete fNeurons->Remove(unit);
1497                 }
1498         }
1499         return;
1500         Bool_t removed;
1501         Int_t nIn, nOut;
1502         AliITSneuron *enemy = 0;
1503         neurons.Reset();
1504         while ( (unit = (AliITSneuron*)neurons.Next()) ) {
1505                 nIn = (Int_t)unit->Inner()->InnerOf()->GetSize();
1506                 nOut = (Int_t)unit->Outer()->OuterOf()->GetSize();
1507                 if (nIn < 2 && nOut < 2) continue;
1508                 removed = kFALSE;
1509                 if (nIn > 1) {
1510                         TObjArrayIter competing(unit->Inner()->InnerOf());
1511                         while ( (enemy = (AliITSneuron*)competing.Next()) ) {
1512                                 if (unit->Activation() > enemy->Activation()) {
1513                                         enemy->Inner()->InnerOf()->Remove(enemy);
1514                                         enemy->Outer()->OuterOf()->Remove(enemy);
1515                                         delete fNeurons->Remove(enemy);
1516                                 }
1517                                 else {
1518                                         unit->Inner()->InnerOf()->Remove(unit);
1519                                         unit->Outer()->OuterOf()->Remove(unit);
1520                                         delete fNeurons->Remove(unit);
1521                                         removed = kTRUE;
1522                                         break;
1523                                 }
1524                         }
1525                         if (removed) continue;
1526                 }
1527                 if (nOut > 1) {
1528                         TObjArrayIter competing(unit->Outer()->OuterOf());
1529                         while ( (enemy = (AliITSneuron*)competing.Next()) ) {
1530                                 if (unit->Activation() > enemy->Activation()) {
1531                                         enemy->Inner()->InnerOf()->Remove(enemy);
1532                                         enemy->Outer()->OuterOf()->Remove(enemy);
1533                                         delete fNeurons->Remove(enemy);
1534                                 }
1535                                 else {
1536                                         unit->Inner()->InnerOf()->Remove(unit);
1537                                         unit->Outer()->OuterOf()->Remove(unit);
1538                                         delete fNeurons->Remove(unit);
1539                                         removed = kTRUE;
1540                                         break;
1541                                 }
1542                         }
1543                 }
1544         }
1545  }
1546  
1547 //__________________________________________________________________________________
1548 Int_t AliITStrackerANN::StoreTracks()
1549 {
1550         // Stores the tracks found in a single neural tracking step.
1551         // In order to do this, it sects each neuron which has a point
1552         // in the first layer.
1553         // Then 
1554         
1555         // if not initialized, the tracks TobjArray is initialized
1556         if (!fFoundTracks) fFoundTracks = new TObjArray;
1557         
1558         Int_t i, check, stored = 0;
1559         Double_t testAct = 0;
1560         AliITSneuron *unit = 0, *cursor = 0, *fwd = 0;
1561         AliITSnode *node = 0;
1562         TObjArrayIter iter(fNeurons), *fwdIter;
1563         TObjArray *removedUnits = new TObjArray(0);
1564         removedUnits->SetOwner(kFALSE);
1565         AliITStrackANN annTrack(fNLayers);
1566         
1567         for (;;) {
1568                 unit = (AliITSneuron*)iter.Next();
1569                 if (!unit) break;
1570                 if (unit->Inner()->GetLayer() > 0) continue;
1571                 annTrack.SetNode(unit->Inner()->GetLayer(), unit->Inner());
1572                 annTrack.SetNode(unit->Outer()->GetLayer(), unit->Outer());
1573                 node = unit->Outer();
1574                 removedUnits->AddLast(unit);
1575                 while (node) {
1576                         testAct = fActMinimum;
1577                         fwdIter = (TObjArrayIter*)node->InnerOf()->MakeIterator();
1578                         fwd = 0;
1579                         for (;;) {
1580                                 cursor = (AliITSneuron*)fwdIter->Next();
1581                                 if (!cursor) break;
1582                                 if (cursor->Used()) continue;
1583                                 if (cursor->Activation() >= testAct) {
1584                                         testAct = cursor->Activation();
1585                                         fwd = cursor;
1586                                 }
1587                         }
1588                         if (!fwd) break;
1589                         removedUnits->AddLast(fwd);
1590                         node = fwd->Outer();
1591                         annTrack.SetNode(node->GetLayer(), node);
1592                 }
1593                 check = annTrack.CheckOccupation();
1594                 if (check >= 6) {
1595                         stored++;
1596                         // FIT
1597                         //if (!RiemannFit(fNLayers, trackitem, param)) continue;
1598                         if (!annTrack.RiemannFit()) continue;
1599                         // initialization of Kalman Filter Tracking
1600                         AliITSclusterV2 *cluster = (AliITSclusterV2*)GetCluster(annTrack[0]->ClusterRef());
1601                         Int_t mod = cluster->GetDetectorIndex();
1602                         Int_t lay, lad, det;
1603                         fGeom->GetModuleId(mod, lay, lad, det);
1604                         Float_t y0 = cluster->GetY();
1605                         Float_t z0 = cluster->GetZ();
1606                         AliITStrackSA* trac = new AliITStrackSA(lay, lad, det, y0, z0, 
1607                                                                                                                                  annTrack.Phi(), annTrack.TanLambda(), 
1608                                                                                                                                  annTrack.Curv(), 1);
1609                         for (Int_t l = 0; l < fNLayers; l++) {
1610                                 if (!annTrack[l]) continue;
1611                                 cluster = (AliITSclusterV2*)GetCluster(annTrack[l]->ClusterRef());
1612                                 if (cluster) trac->AddClusterV2(l, (annTrack[l]->ClusterRef() & 0x0fffffff)>>0);
1613                         }
1614                         AliITStrackV2* ot = new AliITStrackV2(*trac);
1615                         ot->ResetCovariance();
1616                         ot->ResetClusters();
1617                         if (RefitAt(49.,ot,trac)) { //fit from layer 1 to layer 6
1618                                 AliITStrackV2 *otrack2 = new AliITStrackV2(*ot);
1619                                 otrack2->ResetCovariance();
1620                                 otrack2->ResetClusters();
1621                                 //fit from layer 6 to layer 1
1622                                 if (RefitAt(3.7,otrack2,ot)) fFoundTracks->AddLast(otrack2);
1623                         }       
1624                         // end of Kalman Filter fit
1625                         // END FIT
1626                         for (i = 0; i < fNLayers; i++) {
1627                                 //node = (AliITSnode*)removedPoints->At(i);
1628                                 //node->Use();
1629                                 annTrack[i]->Use();
1630                         }
1631                         fwdIter = (TObjArrayIter*)removedUnits->MakeIterator();
1632                         for (;;) {
1633                                 cursor = (AliITSneuron*)fwdIter->Next();
1634                                 if(!cursor) break;
1635                                 cursor->Used() = 1;
1636                         }
1637                 }
1638         }
1639
1640         return stored;
1641 }
1642
1643 Double_t AliITStrackerANN::Weight(AliITSneuron *nAB, AliITSneuron *nBC)
1644 {
1645 /***********************************************************************************
1646  * Calculation of neural weight.
1647  * The implementation of positive neural weight is set only in the case
1648  * of connected units (e.g.: A->B with B->C).
1649  * Given that B is the **common** point. care should be taken to pass 
1650  * as FIRST argument the neuron going "to" B, and
1651  * as SECOND argument the neuron starting "from" B
1652  * anyway, a check is put in order to return 0.0 when arguments are not well given.
1653  ***********************************************************************************/
1654  
1655         if (nAB->Outer() != nBC->Inner()) {
1656                 if (nBC->Outer() == nAB->Inner()) {
1657                         AliITSneuron *temp = nAB;
1658                         nAB = nBC;
1659                         nBC = temp;
1660                         temp = 0;
1661                         if (fMsgLevel >= 3) {
1662                                 Info("Weight", "Switching wrongly ordered arguments.");
1663                         }
1664                 }
1665                 Warning("Weight", "Not connected segments. Returning 0.0");
1666                 return 0.0;
1667         }
1668         
1669         AliITSnode *pA = nAB->Inner();
1670         AliITSnode *pB = nAB->Outer();
1671         AliITSnode *pC = nBC->Outer();
1672         
1673         TVector3 vAB(pB->X() - pA->X(), pB->Y() - pA->Y(), pB->Z() - pA->Z());
1674         TVector3 vBC(pC->X() - pB->X(), pC->Y() - pB->Y(), pC->Z() - pB->Z());
1675
1676         Double_t weight = 1.0 - sin(vAB.Angle(vBC));
1677         return fGain2CostRatio * TMath::Power(weight, fExponent);
1678 }
1679
1680
1681
1682 /******************************************
1683  ******************************************
1684  *** AliITStrackerANN::AliITSnode class ***
1685  ******************************************
1686  ******************************************/
1687  
1688 //__________________________________________________________________________________ 
1689 AliITStrackerANN::AliITSnode::AliITSnode()
1690 : fUsed(kFALSE), fClusterRef(-1), 
1691   fMatches(NULL), fInnerOf(NULL), fOuterOf(NULL), 
1692   fNext(NULL), fPrev(NULL)
1693 {
1694         // Constructor for the embedded 'AliITSnode' class.
1695         // It initializes all pointer-like objects.
1696         
1697         fX = fY = fZ = 0.0;
1698         fEX2 = fEY2 = fEZ2 = 0.0;
1699 }
1700
1701 //__________________________________________________________________________________ 
1702 AliITStrackerANN::AliITSnode::~AliITSnode()
1703 {
1704         // Destructor for the embedded 'AliITSnode' class.
1705         // It should clear the object arrays, but it is possible
1706         // that some objects still are useful after the point deletion
1707         // then the arrays are cleared but their objects are owed by
1708         // another TCollection object, and not deleted.
1709         // For safety reasons, all the pointers are set to zero.
1710         
1711         fInnerOf->SetOwner(kFALSE); 
1712         fInnerOf->Clear(); 
1713         delete fInnerOf;
1714         fInnerOf = 0;
1715         fOuterOf->SetOwner(kFALSE); 
1716         fOuterOf->Clear(); 
1717         delete fOuterOf;
1718         fOuterOf = 0;
1719         fMatches->SetOwner(kFALSE); 
1720         fMatches->Clear();
1721         delete fMatches;
1722         fMatches = 0;
1723         fNext = 0;
1724         fPrev = 0;
1725 }
1726
1727 //__________________________________________________________________________________ 
1728 Double_t AliITStrackerANN::AliITSnode::GetPhi() const
1729 {
1730         // Calculates the 'phi' (azimutal) angle, and returns it
1731         // in the range between 0 and 2Pi radians.
1732         
1733         Double_t q;
1734         q = TMath::ATan2(fY,fX); 
1735         if (q >= 0.) 
1736                 return q;
1737         else 
1738                 return q + TMath::TwoPi();
1739 }
1740  
1741 //__________________________________________________________________________________ 
1742 Double_t AliITStrackerANN::AliITSnode::GetError(Option_t *option)
1743 {
1744         // Returns the error or the square error of 
1745         // values related to the coordinates in different systems.
1746         // The option argument specifies the coordinate error desired:
1747         //
1748         // "R2"     --> error in transverse radius
1749         // "R3"     --> error in spherical radius
1750         // "PHI"    --> error in azimuthal angle
1751         // "THETA"  --> error in polar angle
1752         // "SQ"     --> get the square of error
1753         //
1754         // In order to get the error on the cartesian coordinates
1755         // reference to the inline ErrX2(), ErrY2() adn ErrZ2() methods.
1756         
1757         TString opt(option);
1758         Double_t errorSq = 0.0;
1759         opt.ToUpper();
1760
1761         if (opt.Contains("R2")) {
1762                 errorSq  = fX*fX*fEX2 + fY*fY*fEY2;
1763                 errorSq /= GetR2sq();
1764         }
1765         else if (opt.Contains("R3")) {
1766                 errorSq  = fX*fX*fEX2 + fY*fY*fEY2 + fZ*fZ*fEZ2;
1767                 errorSq /= GetR3sq();
1768         }
1769         else if (opt.Contains("PHI")) {
1770                 errorSq  = fY*fY*fEX2;
1771                 errorSq += fX*fX*fEY2;
1772                 errorSq /= GetR2sq() * GetR2sq();
1773         }
1774         else if (opt.Contains("THETA")) {
1775                 errorSq = fZ*fZ * (fX*fX*fEX2 + fY*fY*fEY2);
1776                 errorSq += GetR2sq() * GetR2sq() * fEZ2;
1777                 errorSq /= GetR3sq() * GetR3sq() * GetR2() * GetR2();
1778         }
1779         
1780         if (!opt.Contains("SQ")) 
1781                 return TMath::Sqrt(errorSq);
1782         else 
1783                 return errorSq;
1784 }
1785
1786
1787
1788 /********************************************
1789  ********************************************
1790  *** AliITStrackerANN::AliITSneuron class ***
1791  ********************************************
1792  ********************************************/
1793  
1794 //__________________________________________________________________________________
1795 AliITStrackerANN::AliITSneuron::AliITSneuron
1796 (AliITSnode *inner, AliITSnode *outer, Double_t minAct, Double_t maxAct)
1797         : fUsed(0), fInner(inner), fOuter(outer) 
1798 {
1799         // Default neuron constructor
1800         fActivation = gRandom->Rndm() * (maxAct-minAct) + minAct;
1801         fGain = new TObjArray;
1802 }
1803  
1804 //__________________________________________________________________________________
1805 Double_t AliITStrackerANN::AliITSneuron::Activate(Double_t temperature)
1806 {
1807         // This computes the new activation of a neuron, and returns
1808         // its activation variation as a consequence of the updating.
1809         // 
1810         // Arguments:
1811         // - the 'temperature' parameter for the neural activation logistic function
1812         // 
1813         // Operations:
1814         // - collects the total gain, by summing the products
1815         //   of the activation of each sequenced unit by the relative weight.
1816         // - collects the total cost, by summing the activations of 
1817         //   all competing units
1818         // - passes the sum of gain - cost to the activation function and
1819         //   calculates the new activation
1820         //   
1821         // Return value:
1822         // - the difference between the old activation and the new one
1823         //   (absolute value)
1824
1825         // local variables
1826         Double_t sumGain = 0.0;      // total contribution from chained neurons
1827         Double_t sumCost = 0.0;      // total contribution from crossing neurons
1828         Double_t input;              // total input
1829         Double_t actOld, actNew;     // old and new values for the activation
1830         AliITSneuron *linked = 0;    // cursor for scanning the neuron arrays (for link check)
1831         AliITSlink *link;            // cursor for scanning the synapses arrays (for link check)
1832         TObjArrayIter *iterator = 0; // pointer to the iterator along the neuron arrays
1833
1834         // sum contributions from the correlated units
1835         iterator = (TObjArrayIter*)fGain->MakeIterator();
1836         for(;;) {
1837                 link = (AliITSlink*)iterator->Next();
1838                 if (!link) break;
1839                 sumGain += link->Contribution();
1840         }
1841         delete iterator;
1842
1843         // sum contributions from the competing units:
1844         // the ones which have the same starting point...
1845         iterator = (TObjArrayIter*)fInner->InnerOf()->MakeIterator();
1846         for (;;) {
1847                 linked = (AliITSneuron*)iterator->Next();
1848                 if (!linked) break;
1849                 if (linked == this) continue;
1850                 sumCost += linked->fActivation;
1851         }
1852         delete iterator;
1853         // ...and the ones which have the same ending point
1854         iterator = (TObjArrayIter*)fOuter->OuterOf()->MakeIterator();
1855         for (;;) {
1856                 linked = (AliITSneuron*)iterator->Next();
1857                 if (!linked) break;
1858                 if (linked == this) continue;
1859                 sumCost += linked->fActivation;
1860         }
1861
1862         // calculate the total input as the difference between gain and cost
1863         input = (sumGain - sumCost) / temperature;
1864         actOld = fActivation;
1865         // calculate the final output
1866 #ifdef NEURAL_LINEAR
1867         if (input <= -2.0 * temperature)
1868                 actNew = 0.0;
1869         else if (input >= 2.0 * temperature)
1870                 actNew = 1.0;
1871         else
1872                 actNew = input / (4.0 * temperature) + 0.5;
1873 #else
1874         actNew = 1.0 / (1.0 + TMath::Exp(-input));
1875 #endif
1876         fActivation = actNew;
1877         
1878         // return the activation variation
1879         return TMath::Abs(actNew - actOld);
1880 }
1881
1882
1883
1884 /******************************************
1885  ******************************************
1886  *** AliITStrackerANN::AliITSlink class ***
1887  ******************************************
1888  ******************************************/
1889  
1890  // No methods defined non-inline
1891  
1892  
1893  
1894  /**********************************************
1895   **********************************************
1896   *** AliITStrackerANN::AliITStrackANN class ***
1897   **********************************************
1898   **********************************************/
1899  
1900 //__________________________________________________________________________________
1901 AliITStrackerANN::AliITStrackANN::AliITStrackANN(Int_t dim) : fNPoints(dim)
1902 {
1903         // Default constructor for the AliITStrackANN class
1904         
1905         fXCenter = 0.0;
1906         fYCenter = 0.0;
1907         fRadius = 0.0;
1908         fCurv = 0.0;
1909         fDTrans = 0.0;
1910         fDLong = 0.0;
1911         fTanLambda = 0.0;
1912         
1913         if (! dim) {
1914                 fNode = 0;
1915         }
1916         else{
1917                 Int_t i = 0;
1918                 fNode = new AliITSnode*[dim];
1919                 for (i = 0; i < dim; i++) fNode[i] = 0;
1920         }
1921 }
1922
1923 //__________________________________________________________________________________
1924 Int_t AliITStrackerANN::AliITStrackANN::CheckOccupation() const
1925 {
1926         // Returns the number of pointers fNode which are not NULL      
1927
1928         Int_t i;         // cursor
1929         Int_t count = 0; // counter for how many points are stored in the track
1930         
1931         for (i = 0; i < fNPoints; i++) {
1932                 if (fNode[i] != NULL) count++;
1933         }
1934         
1935         return count;
1936 }
1937
1938 //__________________________________________________________________________________
1939 Bool_t AliITStrackerANN::AliITStrackANN::RiemannFit()
1940
1941         // Performs the Riemann Sphere fit for the given points to a circle
1942         // and then uses the fit parameters to fit a helix in space.
1943         //
1944         // Return values:
1945         // - kTRUE if all operations have been performed
1946         // - kFALSE if the numbers risk to lead to an arithmetic violation
1947
1948         Int_t i, j, count, dim = fNPoints;
1949         
1950         // First check for all points
1951         count = CheckOccupation();
1952         if (count != fNPoints) {
1953                 Error ("AliITStrackANN::RiemannFit", "CheckOccupations returns %d, fNPoints = %d ==> MISMATCH", count, fNPoints);
1954                 return kFALSE;
1955         }
1956
1957         // matrix of ones
1958         TMatrixD m1(dim,1);
1959         for (i = 0; i < dim; i++) m1(i,0) = 1.0;
1960
1961         // matrix of Rieman projection coordinates
1962         TMatrixD coords(dim,3);
1963         for (i = 0; i < dim; i++) {
1964                 coords(i,0) = fNode[i]->X();
1965                 coords(i,1) = fNode[i]->Y();
1966                 coords(i,2) = fNode[i]->GetR2sq();
1967         }
1968
1969         // matrix of weights
1970         Double_t xterm, yterm, ex, ey;
1971         TMatrixD weights(dim,dim);
1972         for (i = 0; i < dim; i++) {
1973                 xterm = fNode[i]->X() * fNode[i]->GetPhi() - fNode[i]->Y() / fNode[i]->GetR2();
1974                 ex = fNode[i]->ErrX2();
1975                 yterm = fNode[i]->Y() * fNode[i]->GetPhi() + fNode[i]->X() / fNode[i]->GetR2();
1976                 ey = fNode[i]->ErrY2();
1977                 weights(i,i) = fNode[i]->GetR2sq() / (xterm * xterm * ex + yterm * yterm * ey );
1978         }
1979
1980         // weighted sample mean
1981         Double_t meanX = 0.0, meanY = 0.0, meanW = 0.0, sw = 0.0;
1982         for (i = 0; i < dim; i++) {
1983                 meanX += weights(i,i) * coords(i,0);
1984                 meanY += weights(i,i) * coords(i,1);
1985                 meanW += weights(i,i) * coords(i,2);
1986                 sw += weights(i,i);
1987         }
1988         meanX /= sw;
1989         meanY /= sw;
1990         meanW /= sw;
1991
1992         // sample covariance matrix
1993         for (i = 0; i < dim; i++) {
1994                 coords(i,0) -= meanX;
1995                 coords(i,1) -= meanY;
1996                 coords(i,2) -= meanW;
1997         }
1998         TMatrixD coordsT(TMatrixD::kTransposed, coords);
1999         TMatrixD weights4coords(weights, TMatrixD::kMult, coords);
2000         TMatrixD sampleCov(coordsT, TMatrixD::kMult, weights4coords);
2001         for (i = 0; i < 3; i++) {
2002                 for (j = i + 1; j < 3; j++) {
2003                         sampleCov(i,j)  = sampleCov(j,i)  = (sampleCov(i,j) + sampleCov(j,i)) * 0.5;
2004                 }
2005         }
2006
2007         // Eigenvalue problem solving for V matrix
2008         Int_t ileast = 0;
2009         TVectorD eval(3), n(3);
2010         TMatrixD evec = sampleCov.EigenVectors(eval);
2011         if (eval(1) < eval(ileast)) ileast = 1;
2012         if (eval(2) < eval(ileast)) ileast = 2;
2013         n(0) = evec(0, ileast);
2014         n(1) = evec(1, ileast);
2015         n(2) = evec(2, ileast);
2016
2017         // c - known term in the plane intersection with Riemann axes
2018         Double_t c = -(meanX * n(0) + meanY * n(1) + meanW * n(2));
2019
2020         // center and radius of fitted circle
2021         Double_t xc, yc, radius, curv;
2022         xc = -n(0) / (2. * n(2));
2023         yc = -n(1) / (2. * n(2));
2024         radius = (1. - n(2)*n(2) - 4.*c*n(2)) / (4. * n(2) * n(2));
2025         
2026         if (radius <= 0.E0) {
2027                 Error("RiemannFit", "Radius = %f less than zero!!!", radius);
2028                 return kFALSE;
2029         }
2030         radius = TMath::Sqrt(radius);
2031         curv = 1.0 / radius;
2032                 
2033         // evaluating signs for curvature and others
2034         Double_t phi1 = 0.0, phi2, temp1, temp2, phi0, sumdphi = 0.0;
2035         AliITSnode *p = fNode[0];
2036         phi1 = p->GetPhi();
2037         for (i = 1; i < dim; i++) {
2038                 p = (AliITSnode*)fNode[i];
2039                 if (!p) break;
2040                 phi2 = p->GetPhi();
2041                 temp1 = phi1;
2042                 temp2 = phi2;
2043                 if (temp1 > fgkPi && temp2 < fgkPi)
2044                         temp2 += fgkTwoPi;
2045                 else if (temp1 < fgkPi && temp2 > fgkPi)
2046                         temp1 += fgkTwoPi;
2047                 sumdphi += temp2 - temp1;
2048                 phi1 = phi2;
2049         }
2050         if (sumdphi < 0.E0) curv = -curv;
2051         Double_t diff, angle = TMath::ATan2(yc, xc);
2052         if (curv < 0.E0)
2053                 phi0 = angle + 0.5 * TMath::Pi();
2054         else
2055                 phi0 = angle - 0.5 * TMath::Pi();
2056         diff = angle - phi0;
2057
2058         Double_t dt, temp = TMath::Sqrt(xc*xc + yc*yc) - radius;
2059         if (curv >= 0.E0)
2060                 dt = temp;
2061         else
2062                 dt = -temp;
2063         //cout << "Dt = " << dt << endl;
2064         
2065         Double_t halfC = 0.5 * curv, test;
2066         Double_t *s = new Double_t[dim], *zz = new Double_t[dim], *ws = new Double_t[dim];
2067         for (j = 0; j < 6; j++) {
2068                 p = fNode[j];
2069                 if (!p) break;
2070                 //----
2071                 s[j] = (p->GetR2sq() - dt * dt) / (1. + curv * dt);
2072                 if (s[j] < 0.) {
2073                         if (TMath::Abs(s[j]) < 1.E-6) s[j] = 0.;
2074                         else {
2075                                 Error("RiemannFit", "Square root argument error: %17.15g < 0", s[j]);
2076                                 return kFALSE;
2077                         }
2078                 }
2079                 s[j] = TMath::Sqrt(s[j]);
2080                 //cout << "Curv = " << halfC << " --- s[" << j << "] = " << s[j] << endl;
2081                 s[j] *= halfC;
2082                 test = TMath::Abs(s[j]);
2083                 if (test > 1.) {
2084                         if (test <= 1.1) 
2085                                 s[j] = ((s[j] > 0.) ? 0.99999999999 : -0.9999999999);
2086                         else {
2087                                 Error("RiemannFit", "Value too large: %17.15g", s[j]);
2088                                 return kFALSE;
2089                         }
2090                 }
2091                 //----
2092                 zz[j] = p->Z();
2093                 s[j] = TMath::ASin(s[j]) / halfC;
2094                 ws[j] = 1.0 / (p->ErrZ2());
2095         }
2096
2097         // second tep final fit
2098         Double_t s2Sum = 0.0, zSum = 0.0, szSum = 0.0, sSum = 0.0, sumw = 0.0;
2099         for (i = 0; i < dim; i++) {
2100                 s2Sum += ws[i] * s[i] * s[i];
2101                 zSum  += ws[i] * zz[i];
2102                 sSum  += ws[i] * s[i];
2103                 szSum += ws[i] * s[i] * zz[i];
2104                 sumw += ws[i];
2105         }
2106         s2Sum /= sumw;
2107         zSum /= sumw;
2108         sSum /= sumw;
2109         szSum /= sumw;
2110         temp = s2Sum - sSum*sSum;
2111
2112         Double_t dz, tanL;
2113         dz = (s2Sum*zSum - sSum*szSum) / temp;
2114         tanL = (szSum - sSum*zSum) / temp;
2115         
2116         fXCenter = xc;
2117         fYCenter = yc;
2118         fRadius = radius;
2119         fCurv = curv;
2120         fPhi = phi0;
2121         fDTrans = dt;
2122         fDLong = dz;
2123         fTanLambda = tanL;
2124
2125         delete [] s;
2126         delete [] zz;
2127         delete [] ws;
2128         
2129         return kTRUE;
2130 }