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