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