1 // ---------------------------------------------------------------------------------
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
16 // ---------------------------------------------------------------------------------
17 // Author: Alberto Pulvirenti (university of Catania)
18 // Email : alberto.pulvirenti@ct.infn.it
19 // ---------------------------------------------------------------------------------
22 #include <TObjArray.h>
28 #if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,2)
29 #include <TMatrixDEigen.h>
32 #include "AliITSgeom.h"
33 #include "AliITStrackSA.h"
34 #include "AliITSRecPoint.h"
35 #include "AliITStrackV2.h"
37 #include "AliITStrackerANN.h"
39 const Double_t AliITStrackerANN::fgkPi = 3.141592653; // pi
40 const Double_t AliITStrackerANN::fgkHalfPi = 1.570796327; // pi / 2
41 const Double_t AliITStrackerANN::fgkTwoPi = 6.283185307; // 2 * pi
43 ClassImp(AliITStrackerANN)
45 //__________________________________________________________________________________
46 AliITStrackerANN::AliITStrackerANN(const AliITSgeom *geom, Int_t msglev)
47 : AliITStrackerV2(geom), fMsgLevel(msglev)
49 /**************************************************************************
51 CONSTRUCTOR (standard-to-use version)
54 1) the ITS geometry used in the generated event
55 2) the flag for log-messages writing
57 The AliITSgeometry is used along the class,
58 in order to translate the local AliITSRecPoint coordinates
59 into the Global reference frame, which is necessary for the
60 Neural Network algorithm to operate.
61 In case of serialized use, the log messages should be excluded,
62 in order to save real execution time.
65 - all pointer data members are initialized
66 (if possible, otherwise are set to NULL)
67 - all numeric data members are initialized to
68 values which allow the Neural Network to operate
69 even if nothing is externally set.
71 NOTE: it is possible that tracking an event
72 with these default values results in a non-sense.
74 **************************************************************************/
79 fGeom = (AliITSgeom*)geom;
81 // Initialize the array of first module indexes per layer
82 fNLayers = geom->GetNlayers();
83 fFirstModInLayer = new Int_t[fNLayers + 1];
84 for (i = 0; i < fNLayers; i++) {
85 fFirstModInLayer[i] = fGeom->GetModuleIndex(i + 1, 1, 1);
87 fFirstModInLayer[fNLayers] = geom->GetIndexMax();
89 // initialization: no curvature cut steps
93 // initialization: 4 sectors (one for each quadrant)
95 fSectorWidth = fgkHalfPi;
97 // initialization: theta offset of 20 degrees
98 fPolarInterval = 20.0;
100 // initialization: array structure not defined
101 fStructureOK = kFALSE;
103 // initialization: vertex in the origin
108 // initialization: uninitialized point array
111 // initialization: very large (unuseful) cut values
113 for (ilayer = 0; ilayer < 6; ilayer++) {
114 fThetaCut2D[ilayer] = TMath::Pi();
115 fThetaCut3D[ilayer] = TMath::Pi();
116 fHelixMatchCut[ilayer] = 1.0;
119 // initialization: inictial activation range between 0.3 and 0.7
123 // initialization: neural network operation & weights
125 fStabThreshold = 0.001;
126 fGain2CostRatio = 1.0;
130 // initialization: uninitialized array of neurons
133 // initialization: uninitialized array of tracks
137 //__________________________________________________________________________________
138 void AliITStrackerANN::SetCuts
139 (Int_t ncurv, Double_t *curv, Double_t *theta2D, Double_t *theta3D, Double_t *helix)
141 /**************************************************************************
145 Allows for the definition of all kind of geometric cuts
146 which have been studied in for the creation of a neuron
147 from a pair of clusters C1 and C2 on consecutive layers.
148 Neuron will be created only if the pair passes ALL of these cuts.
150 At the moment, we define 4 kinds of geometrical cuts:
151 a) cut on the difference of the polar 'theta' angle;
152 b) cut on the angle between origin->C1 and C1->C2 in space;
153 c) cut on the curvature of the circle passing
154 through C1, C2 and the primary vertex;
155 d) cut on heli-matching of the same three points.
158 1) the number of curvature cut steps
159 2) the array of curvature cuts for each step
160 (its dimension is given by the argument 1)
161 3) array of 5 cut values (one for each consecutive lauer pair)
163 4) array of 5 cut values (one for each consecutive lauer pair)
165 5) array of 5 cut values (one for each consecutive lauer pair)
169 - gets the values for each cut and stores them in data members
170 - in the case of curvature cuts, the cut array
171 (whose size is not fixed) is allocated
173 NOTE: in case that the user wants to set onyl ONE of the 4 cuts array,
174 he can simply pass NULL arguments for other cuts, and (eventually)
175 a ZERO as the first argument (if curvature cuts have not to be set).
177 Anyway, all the cuts have to be set at least once.
179 **************************************************************************/
184 /*** Curvature cut setting ***/
186 // first of all, the curvature cuts are sorted in increasing order
187 // (from the smallest to the largest one)
188 Int_t *ind = new Int_t[ncurv];
189 TMath::Sort(ncurv, curv, ind, kFALSE);
190 // then, the curvature cut array is allocated and filled
191 // (a message with the list of defined cuts can be optionally shown)
192 fCurvCut = new Double_t[ncurv];
193 if (fMsgLevel >= 1) cout << "Number of curvature cuts: " << ncurv << endl;
194 for (i = 0; i < ncurv; i++) {
195 fCurvCut[i] = curv[ind[i]];
196 if (fMsgLevel >= 1) cout << " - " << fCurvCut[i] << endl;
200 /*** 'Fixed' cuts setting ***/
202 // checks what cuts have to be set
203 Bool_t doTheta2D = (theta2D != 0);
204 Bool_t doTheta3D = (theta3D != 0);
205 Bool_t doHelix = (helix != 0);
206 // sets the cuts for all layer pairs
207 for (i = 0; i < fNLayers - 1; i++) {
208 if (doTheta2D) fThetaCut2D[i] = theta2D[i];
209 if (doTheta3D) fThetaCut3D[i] = theta3D[i];
210 if (doHelix) fHelixMatchCut[i] = helix[i];
212 // if required, lists the cuts
213 if (!fMsgLevel < 2) return;
214 cout << "Theta 2D cuts: ";
216 cout << "<not set>" << endl;
220 for (i = 0; i < fNLayers - 1; i++) {
221 cout << "For layers " << i << " --> " << i + 1;
222 cout << " cut = " << fThetaCut2D[i] << endl;
225 cout << "---" << endl;
226 cout << "Theta 3D cuts: ";
228 cout << "<not set>" << endl;
232 for (i = 0; i < fNLayers - 1; i++) {
233 cout << "For layers " << i << " --> " << i + 1;
234 cout << " cut = " << fThetaCut3D[i] << endl;
237 cout << "---" << endl;
238 cout << "Helix-match cuts: ";
240 cout << "<not set>" << endl;
244 for (i = 0; i < fNLayers - 1; i++) {
245 cout << "For layers " << i << " --> " << i + 1;
246 cout << " cut = " << fHelixMatchCut[i] << endl;
249 cout << "---" << endl;
252 //__________________________________________________________________________________
253 Bool_t AliITStrackerANN::GetGlobalXYZ
255 Double_t &x, Double_t &y, Double_t &z, Double_t &ex2, Double_t &ey2, Double_t &ez2)
257 /**************************************************************************
259 LOCAL TO GLOBAL TRANSLATOR
261 Taking information from the ITS geometry stored in the class,
262 gets a stored AliITScluster and calculates it coordinates
263 and errors in the global reference frame.
264 These values are stored in the variables,
265 which are passed by reference.
268 1) reference index for the cluster to use
269 2) (by reference) place to store the global X-coord into
270 3) (by reference) place to store the global Y-coord into
271 4) (by reference) place to store the global Z-coord into
272 5) (by reference) place to store the global X-coord error into
273 6) (by reference) place to store the global Y-coord error into
274 7) (by reference) place to store the global Z-coord error into
277 essentially, determines the ITS module index from the
278 detector index of the AliITSRecPoint object, and extracts
279 the roto-translation from the ITS geometry, to convert
280 the local module coordinates into the global ones.
283 - kFALSE if the given cluster index points to a non-existing
284 cluster, or if the layer number makes no sense (< 0 or > 6).
285 - otherwise, kTRUE (meaning a successful operation).
287 **************************************************************************/
289 // checks if the layer is correct
290 Int_t ilayer = (refIndex & 0xf0000000) >> 28;
291 if (ilayer < 0 || ilayer >= fNLayers) {
292 Error("GetGlobalXYZ", "Wrong layer number: %d [range: %d - %d]", ilayer, 0, fNLayers);
295 // checks if the referenced cluster exists and corresponds to the passed reference
296 AliITSRecPoint *refCluster = (AliITSRecPoint*) GetCluster(refIndex);
298 Error("GetGlobalXYZ", "Cluster not found for index %d", refIndex);
302 // determine the detector number
303 Int_t detID = refCluster->GetDetectorIndex() + fFirstModInLayer[ilayer];
305 // get rotation matrix
307 fGeom->GetRotMatrix(detID, rot);
309 // get translation vector
311 fGeom->GetTrans(detID, tx, ty, tz);
313 // determine r and phi for the reference conversion
314 Double_t r = -(Double_t)tx * rot[1] + (Double_t)ty * rot[0];
315 if (ilayer == 0) r = -r;
316 Double_t phi = TMath::ATan2(rot[1],rot[0]);
317 if (ilayer == 0) phi -= fgkPi;
319 // sets values for X, Y, Z in global coordinates and their errors
320 Double_t cosPhi = TMath::Cos(phi);
321 Double_t sinPhi = TMath::Sin(phi);
322 x = r*cosPhi + refCluster->GetY()*sinPhi;
323 y = -r*sinPhi + refCluster->GetY()*cosPhi;
324 z = refCluster->GetZ();
325 ex2 = refCluster->GetSigmaY2()*sinPhi*sinPhi;
326 ey2 = refCluster->GetSigmaY2()*cosPhi*cosPhi;
327 ez2 = refCluster->GetSigmaZ2();
332 //__________________________________________________________________________________
333 AliITStrackerANN::AliITSnode* AliITStrackerANN::AddNode(Int_t refIndex)
335 /**************************************************************************
337 GENERATOR OF NEURAL NODES
339 Fills the array of neural 'nodes', which are the ITS clusters
340 translated in the global reference frame.
341 Given that the global coordinates are used many times, they are
342 stored in a well-defined structure, in the form of an embedded class.
343 Moreover, this class allows a faster navigation among points
344 and neurons, by means of some object arrays, storing only the
345 neurons which start from, or end to, the given node.
346 Finally, each node contains all the other nodes which match it
347 with respect to the fixed walues, in order to perform a faster
348 neuron-creation phase.
351 1) reference index of the correlated AliITSRecPoint object
354 - allocates the new AliITSnode objects
355 - initializes its object arrays
356 - from the global coordinates, calculates the
357 'phi' and 'theta' coordinates, in order to store it
358 into the correct theta-slot and azimutal sector.
361 - the pointer of the creater AliITSnode object
362 - in case of errors, a waring is given and a NULL is returned
364 **************************************************************************/
366 // create object and set the reference
367 AliITSnode *node = new AliITSnode;
369 Warning("AddNode", "Error occurred when allocating AliITSnode");
372 node->ClusterRef() = refIndex;
374 // calls the conversion function, which makes also some checks
375 // (layer number within range, existence of referenced cluster)
378 node->X(), node->Y(), node->Z(),
379 node->ErrX2(), node->ErrY2(), node->ErrZ2()
382 // initializes the object arrays
383 node->Matches() = new TObjArray;
384 node->InnerOf() = new TObjArray;
385 node->OuterOf() = new TObjArray;
387 // finds azimutal and polar sector (in degrees)
388 Double_t phi = node->GetPhi() * 180.0 / fgkPi;
389 Double_t theta = node->GetTheta() * 180.0 / fgkPi;
390 Int_t isector = (Int_t)(phi / fSectorWidth);
391 Int_t itheta = (Int_t)theta;
392 Int_t ilayer = (refIndex & 0xf0000000) >> 28;
394 // selects the right TObjArray to store object into
395 TObjArray *sector = (TObjArray*)fNodes[ilayer][itheta]->At(isector);
396 sector->AddLast(node);
401 //__________________________________________________________________________________
402 void AliITStrackerANN::CreateArrayStructure(Int_t nsectors)
404 /**************************************************************************
406 ARRAY STRUCTURE CREATOR
408 Creates a structure of nested TObjArray's where the AliITSnode's
410 - the first level is made by 6 arrays (one for each layer)
411 - the second level is made by 180 arrays (one for each int part of theta)
412 - the third level is made by a variable number of arrays
413 (one for each azimutal sector)
416 1) the number of azimutal sectors
419 - calculates the width of each sector, from the argument
420 - allocates and initializes all array levels
421 - sets a flag which tells the user if this NECESSARY operation
422 has been performed (it is needed BEFORE performing tracking)
424 **************************************************************************/
426 // Set the number of sectors and their width.
427 fSectorNum = nsectors;
428 fSectorWidth = 360.0 / (Double_t)fSectorNum;
429 if (fMsgLevel >= 2) {
430 cout << fSectorNum << " sectors --> sector width (degrees) = " << fSectorWidth << endl;
433 // Meaningful indexes
434 Int_t ilayer, isector, itheta;
436 // Mark for the created objects
437 TObjArray *sector = 0;
439 // First index: layer
440 fNodes = new TObjArray**[fNLayers];
441 for (ilayer = 0; ilayer < fNLayers; ilayer++) {
442 fNodes[ilayer] = new TObjArray*[180];
443 for (itheta = 0; itheta < 180; itheta++) fNodes[ilayer][itheta] = 0;
444 for (itheta = 0; itheta < 180; itheta++) {
445 fNodes[ilayer][itheta] = new TObjArray(nsectors);
446 for (isector = 0; isector < nsectors; isector++) {
447 sector = new TObjArray;
449 fNodes[ilayer][itheta]->AddAt(sector, isector);
454 // Sets a checking flag to TRUE.
455 // This flag is checked before filling up the arrays with the points.
456 fStructureOK = kTRUE;
459 //__________________________________________________________________________________
460 Int_t AliITStrackerANN::ArrangePoints(char *exportFile)
462 /**************************************************************************
466 This function assembles the operation from the other above methods,
467 and fills the arrays with the clusters already stored in the
468 layers of the tracker.
469 Then, in order to use this method, the user MUSTs call LoadClusters()
473 1) string for a file name where the global ccordinates
474 of all points can be exported (optional).
475 If this file must not be created, simply pass a NULL argument
478 - for each AliITSRecPoint in each AliITSlayer, a ne AliITSnode
479 is created and stored in the correct location.
482 - the number of stored points
483 - when errors occur, or no points are found, 0 is returned
485 **************************************************************************/
487 // Check if the array structure has been created
489 Error("ArrangePoints", "Structure NOT defined. Call CreateArrayStructure() first");
493 // meaningful indexes
494 Int_t ientry, ilayer, nentries = 0, index;
497 // if the argument is not NULL, a file is opened
498 fstream file(exportFile, ios::out);
499 if (!exportFile || file.fail()) {
504 // scan all layers for node creation
505 for (ilayer = 0; ilayer < fNLayers; ilayer++) {
506 nPtsLayer = GetNumberOfClustersLayer(ilayer);
507 if (fMsgLevel >= 1) {
508 cout << "Layer " << ilayer << " --> " << nPtsLayer << " clusters" << endl;
510 for (ientry = 0; ientry < nPtsLayer; ientry++) {
511 // calculation of cluster index : (Bit mask LLLLIIIIIIIIIIII)
512 // (L = bits used for layer)
513 // (I = bits used for position in layer)
514 index = ilayer << 28;
516 // add new AliITSnode object
517 AliITSnode *n = AddNode(index);
518 if ( (n != NULL) && exportFile ) {
519 file << index << ' ' << n->X() << ' ' << n->Y() << ' ' << n->Z() << endl;
522 nentries += nPtsLayer;
525 // a conventional final message is put at the end of file
527 file << "-1 0.0 0.0 0.0" << endl;
531 // returns the number of points processed
535 //__________________________________________________________________________________
536 void AliITStrackerANN::StoreOverallMatches()
538 /**************************************************************************
542 Once the nodes have been created, a firs analysis is to check
543 what pairs will satisfy at least the 'fixed' cuts (theta, helix-match)
544 and the most permissive (= larger) curvature cut.
545 All these node pairs are suitable for neuron creation.
546 In fact, when performing a Neural Tracking step, the only further check
547 will be a check against the current curvature step, while the other
549 For thi purpose, each AliITSnode has a data member, named 'fMatches'
550 which contains references to all other AliITSnodes in the successive layer
551 that form, with it, a 'good' pair, with respect to the above cited cuts.
552 Then, in each step for neuron creation, the possible neurons starting from
553 each node will be searched ONLY within the nodes referenced in fMatches.
554 This, of course, speeds up a lot the neuron creation procedure, at the
555 cost of some memory occupation, which results not to be critical.
558 - for each AliITSnode, matches are found according to the criteria
559 expressed above, and stored in the node->fMatches array
561 **************************************************************************/
563 // meaningful counters
564 Int_t ilayer, isector, itheta1, itheta2, check;
565 TObjArray *list1 = 0, *list2 = 0;
566 AliITSnode *node1 = 0, *node2 = 0;
567 Double_t thetaMin, thetaMax;
570 // Scan for each sector
571 for (isector = 0; isector < fSectorNum; isector++) {
572 // sector is chosen once for both lists
573 for (ilayer = 0; ilayer < fNLayers - 1; ilayer++) {
574 for (itheta1 = 0; itheta1 < 180; itheta1++) {
575 list1 = (TObjArray*)fNodes[ilayer][itheta1]->At(isector);
576 TObjArrayIter iter1(list1);
577 while ( (node1 = (AliITSnode*)iter1.Next()) ) {
578 if (node1->IsUsed()) continue;
579 // clear an eventually already present array
580 // node1->Matches()->Clear();
581 // get the global coordinates and defines the theta interval from cut
582 thetaMin = (node1->GetTheta() * 180.0 / fgkPi) - fPolarInterval;
583 thetaMax = (node1->GetTheta() * 180.0 / fgkPi) + fPolarInterval;
584 imin = (Int_t)thetaMin;
585 imax = (Int_t)thetaMax;
586 if (imin < 0) imin = 0;
587 if (imax > 179) imax = 179;
588 // loop on the selected theta slots
589 for (itheta2 = imin; itheta2 <= imax; itheta2++) {
590 list2 = (TObjArray*)fNodes[ilayer + 1][itheta2]->At(isector);
591 TObjArrayIter iter2(list2);
592 while ( (node2 = (AliITSnode*)iter2.Next()) ) {
593 check = PassAllCuts(node1, node2, fCurvNum - 1, fVertexX, fVertexY, fVertexZ);
595 node1->Matches()->AddLast(node2);
597 } // while (node2...)
598 } // for (itheta2...)
599 } // while (node1...)
602 } // for (isector...)
605 //__________________________________________________________________________________
606 Int_t AliITStrackerANN::PassAllCuts
607 (AliITSnode *inner, AliITSnode *outer, Int_t curvStep,
608 Double_t vx, Double_t vy, Double_t vz)
610 // ***********************************************************************************
612 // This check is called in the above method for finding the matches of each node
613 // It check the passed point pair against all the fixed cuts and a specified
614 // curvature cut, among all the ones which have been defined.
615 // The cuts need a vertex-constraint, which is not absolute, but it is passed
619 // 1) the point in the inner layer
620 // 2) the point in the outer layer
621 // 3) curvature step for the curvature cut check (preferably the last)
622 // 4) X of the used vertex
623 // 5) Y of the used vertex
624 // 6) Z of the used vertex
627 // - if necessary, swaps the two points
628 // (the first must be in the innermost of the two layers)
629 // - checks for the theta cuts
630 // - calculates the circle passing through the vertex
631 // and the given points and checks for the curvature cut
632 // - using the radius calculated there, checks for the helix-math cut
635 // 0 - All cuts passed
636 // 1 - theta 2D cut not passed
637 // 2 - theta 3D cut not passed
638 // 3 - curvature calculated but cut not passed
639 // 4 - curvature not calculated (division by zero)
640 // 5 - helix cut not passed
641 // 6 - curvature inxed out of range
643 // ***********************************************************************************
645 // Check for curvature index
646 if (curvStep < 0 || curvStep >= fCurvNum) return 6;
648 // Swap points in order that r1 < r2
649 AliITSnode *temp = 0;
650 if (outer->GetLayer() < inner->GetLayer()) {
656 // All cuts are variable according to the layer of the
657 // innermost point (the other point will surely be
658 // in the further one, because we don't check poin pairs
659 // which are not in adjacent layers)
660 // The reference is given by the innermost point.
661 Int_t layRef = inner->GetLayer();
663 // The calculations in the transverse plane are made in
664 // a shifted reference frame, whose origin corresponds to
665 // the reference point passed in the argument.
666 Double_t xIn = inner->X() - vx;
667 Double_t xOut = outer->X() - vx;
668 Double_t yIn = inner->Y() - vy;
669 Double_t yOut = outer->Y() - vy;
670 Double_t zIn = inner->Z() - vz;
671 Double_t zOut = outer->Z() - vz;
672 Double_t rIn = TMath::Sqrt(xIn*xIn + yIn*yIn);
673 Double_t rOut = TMath::Sqrt(xOut*xOut + yOut*yOut);
675 // Check for theta cut.
676 // There are two different angular cuts:
677 // one w.r. to the angle in the 2-dimensional r-z plane...
679 TVector3 origin2innerRZ(zIn, rIn, 0.0);
680 TVector3 inner2outerRZ(zOut - zIn, rOut - rIn, 0.0);
681 dthetaRZ = origin2innerRZ.Angle(inner2outerRZ) * 180.0 / fgkPi;
682 if (dthetaRZ > fThetaCut2D[layRef]) {
684 // r-z theta cut not passed ---> 1
686 // ...and another w.r. to the angle in the 3-dimensional x-y-z space
688 TVector3 origin2innerXYZ(xIn, yIn, zIn);
689 TVector3 inner2outerXYZ(xOut - xIn, yOut - yIn, zOut - zIn);
690 dthetaXYZ = origin2innerXYZ.Angle(inner2outerXYZ) * 180.0 / fgkPi;
691 if (dthetaXYZ > fThetaCut3D[layRef]) {
693 // x-y-z theta cut not passed ---> 2
696 // Calculation & check of curvature
697 Double_t dx = xIn - xOut;
698 Double_t dy = yIn - yOut;
699 Double_t num = 2.0 * (xIn*yOut - xOut*yIn);
700 Double_t den = rIn*rOut*sqrt(dx*dx + dy*dy);
703 curv = TMath::Abs(num / den);
704 if (curv > fCurvCut[curvStep]) {
706 // curvature too large for cut ---> 3
710 Error("PassAllCuts", "Curvature calculations gives zero denominator");
712 // error: denominator = 0 ---> 4
715 // Calculation & check of helix matching
716 Double_t helMatch = 0.0;
717 Double_t arcIn = 2.0 * rIn * curv;
718 Double_t arcOut = 2.0 * rOut * curv;
719 if (arcIn > -1.0 && arcIn < 1.0)
720 arcIn = TMath::ASin(arcIn);
722 arcIn = ((arcIn > 0.0) ? 0.5 : 1.5) * TMath::Pi();
723 if (arcOut > -1.0 && arcOut < 1.0)
724 arcOut = TMath::ASin(arcOut);
726 arcOut = ((arcOut > 0.0) ? 0.5 : 1.5) * TMath::Pi();
728 arcOut /= 2.0 * curv;
729 if (arcIn == 0.0 || arcOut == 0.0) {
730 Error("PassAllCuts", "Calculation returns zero-length arcs: l1=%f, l2=%f", arcIn, arcOut);
732 // error: circumference arcs seem to equal zero ---> 4
734 helMatch = TMath::Abs(zIn / arcIn - zOut / arcOut);
735 if (helMatch > fHelixMatchCut[layRef]) {
737 // helix match cut not passed ---> 5
740 // ALL CUTS PASSED ---> 0
744 //__________________________________________________________________________________
745 Bool_t AliITStrackerANN::PassCurvCut
746 (AliITSnode *inner, AliITSnode *outer,
748 Double_t vx, Double_t vy, Double_t vz)
750 //***********************************************************************************
752 // This method operates essentially like the above one, but it is used
753 // during a single step of Neural Tracking, where the curvature cut
755 // Then, not necessaryly all the nodes stored in the fMatches array
756 // will be suitable for neuron creation in an intermediate step.
758 // It has the same arguments of the PassAllCuts() method, but
759 // the theta cut is not checked.
760 // Moreover, it has a boolean return value.
762 //***********************************************************************************
764 // Check for curvature index
765 if (curvStep < 0 || curvStep >= fCurvNum) return 6;
767 // Find the reference layer
768 Int_t layIn = inner->GetLayer();
769 Int_t layOut = outer->GetLayer();
770 Int_t layRef = (layIn < layOut) ? layIn : layOut;
772 // The calculations in the transverse plane are made in
773 // a shifted reference frame, whose origin corresponds to
774 // the reference point passed in the argument.
775 Double_t xIn = inner->X() - vx;
776 Double_t xOut = outer->X() - vx;
777 Double_t yIn = inner->Y() - vy;
778 Double_t yOut = outer->Y() - vy;
779 Double_t zIn = inner->Z() - vz;
780 Double_t zOut = outer->Z() - vz;
781 Double_t rIn = TMath::Sqrt(xIn*xIn + yIn*yIn);
782 Double_t rOut = TMath::Sqrt(xOut*xOut + yOut*yOut);
784 // Calculation & check of curvature
785 Double_t dx = xIn - xOut;
786 Double_t dy = yIn - yOut;
787 Double_t num = 2.0 * (xIn*yOut - xOut*yIn);
788 Double_t den = rIn*rOut*sqrt(dx*dx + dy*dy);
792 curv = TMath::Abs(num / den);
793 if (curv > fCurvCut[curvStep]) return kFALSE;
801 curv = TMath::Abs(num / den);
802 if (curv > fCurvCut[curvStep]) {
807 Error("PassAllCuts", "Curvature calculations gives zero denominator");
811 // Calculation & check of helix matching
812 Double_t helMatch = 0.0;
813 Double_t arcIn = 2.0 * rIn * curv;
814 Double_t arcOut = 2.0 * rOut * curv;
815 if (arcIn > -1.0 && arcIn < 1.0)
816 arcIn = TMath::ASin(arcIn);
818 arcIn = ((arcIn > 0.0) ? 0.5 : 1.5) * TMath::Pi();
819 if (arcOut > -1.0 && arcOut < 1.0)
820 arcOut = TMath::ASin(arcOut);
822 arcOut = ((arcOut > 0.0) ? 0.5 : 1.5) * TMath::Pi();
824 arcOut /= 2.0 * curv;
825 if (arcIn == 0.0 || arcOut == 0.0) {
826 Error("PassAllCuts", "Calculation returns zero-length arcs: l1=%f, l2=%f", arcIn, arcOut);
828 // error: circumference arcs seem to equal zero ---> 4
830 helMatch = TMath::Abs(zIn / arcIn - zOut / arcOut);
831 return (helMatch <= fHelixMatchCut[layRef]);
834 //__________________________________________________________________________________
835 void AliITStrackerANN::PrintMatches(Bool_t stop)
837 // Prints the list of points which appear to match
838 // each one of them, according to the preliminary
840 // The arguments states if a pause is required after printing
841 // the matches for each one. In this case, a keypress is needed.
843 TObjArray *sector = 0;
844 Int_t ilayer, isector, itheta, nF;
845 AliITSnode *node1 = 0, *node2 = 0;
846 //AliITSRecPoint *cluster1 = 0, *cluster2 = 0;
848 for (ilayer = 0; ilayer < 6; ilayer++) {
849 for (isector = 0; isector < fSectorNum; isector++) {
850 for (itheta = 0; itheta < 180; itheta++) {
851 sector = (TObjArray*)fNodes[ilayer][itheta]->At(isector);
852 TObjArrayIter points(sector);
853 while ( (node1 = (AliITSnode*)points.Next()) ) {
854 nF = (Int_t)node1->Matches()->GetEntries();
855 cout << "Node layer: " << node1->GetLayer() << " --> ";
857 cout << "NO Matches!!!" << endl;
860 cout << nF << " Matches" << endl;
861 cout << "Reference cluster: " << hex << node1->ClusterRef() << endl;
862 TObjArrayIter matches(node1->Matches());
863 while ( (node2 = (AliITSnode*)matches.Next()) ) {
864 cout << "Match with " << hex << node2->ClusterRef() << endl;
867 cout << "Press a key" << endl;
876 //__________________________________________________________________________________
877 void AliITStrackerANN::ResetNodes(Int_t isector)
879 /***********************************************************************************
881 NODE NEURON ARRAY CLEANER
883 After a neural tracking step, this method
884 clears the arrays 'fInnerOf' and 'fOuterOf' of each AliITSnode
887 - the sector where the operation is being executed
889 ***********************************************************************************/
891 Int_t ilayer, itheta;
892 TObjArray *sector = 0;
893 AliITSnode *node = 0;
894 for (ilayer = 0; ilayer < fNLayers; ilayer++) {
895 for (itheta = 0; itheta < 180; itheta++) {
896 sector = (TObjArray*)fNodes[ilayer][itheta]->At(isector);
897 TObjArrayIter iter(sector);
899 node = (AliITSnode*)iter.Next();
901 node->InnerOf()->Clear();
902 node->OuterOf()->Clear();
904 delete node->InnerOf();
905 delete node->OuterOf();
906 node->InnerOf() = new TObjArray;
907 node->OuterOf() = new TObjArray;
914 //__________________________________________________________________________________
915 Int_t AliITStrackerANN::CreateNeurons
916 (AliITSnode *node, Int_t curvStep, Double_t vx, Double_t vy, Double_t vz)
918 // This method is used to create alle suitable neurons starting from
919 // a single AliITSnode. Each unit is also stored in the fInnerOf array
920 // of the passed node, and in the fOuterOf array of the other neuron edge.
921 // In the new implementation of the intermediate check steps, a further one
922 // is made, which chechs how well a helix is matched by three points
923 // in three consecutive layers.
924 // Then, a vertex-constrained check is made with vertex located
925 // in a layer L, for points in layers L+1 and L+2.
927 // In order to do this, the creator works recursively, in a tree-visit like operation.
928 // The neurons are effectively created only if the node argument passed is in
929 // the 5th layer (they are created between point of 5th and 6th layer).
930 // If the node is in an inner layer, its coordinates are passet as vertex for a nested
931 // call of the same function in the next two layers.
935 // 2) current curvature cut step
936 // 3) X of referenced temporary (not primary) vertex
937 // 4) Y of referenced temporary (not primary) vertex
938 // 5) Z of referenced temporary (not primary) vertex
941 // - if the layer is the 5th, neurons are created with nodes
942 // in the fMatches array of the passed node
943 // - otherwise, the X, Y, Z of the passed node are given as
944 // vertex and the same procedure is recursively called for all
945 // nodes in the fMatches array of the passed one.
948 // - the total number of neurons created from the passed one
949 // summed with all neurons created from all nodes well matched with it
950 // (assumes a meaning only for nodes in the first layer)
953 Int_t made = 0; // counter
954 Bool_t found = 0; // flag
955 AliITSnode *match = 0; // cursor for a AliITSnode array
956 AliITSneuron *unit = 0; // cursor for a AliITSneuron array
958 // --> Case 0: the passed node has already been used
959 // as member of a track found in a previous step.
960 // In this case, of course, the function exits.
961 if (node->IsUsed()) return 0;
963 // --> Case 1: there are NO well-matched points.
964 // This can happen in all ITS layers, but it happens **for sure**
965 // for a node in the 'last' layer.
966 // Even in this case, the function exits.
967 if (node->Matches()->IsEmpty()) return 0;
969 // --> Case 2: there are well-matched points.
970 // In this case, the function creates a neuron for each
971 // well-matched pair (according to the cuts for the current step)
972 // Moreover, before storing the neuron, a check is necessary
973 // to avoid the duplicate creation of the same neuron twice.
974 // (This could happen if the 3 last arguments of the function
975 // are close enough to cause a good match for the current step
976 // between two points, independently of their difference).
977 // Finally, a node is skipped if it has already been used.
978 // For each matched point for which a neuron is created, the procedure is
979 // recursively called.
980 TObjArrayIter matches(node->Matches());
981 while ( (match = (AliITSnode*)matches.Next()) ) {
982 if (match->IsUsed()) continue;
983 if (!PassCurvCut(node, match, curvStep, vx, vy, vz)) continue;
985 if (!node->InnerOf()->IsEmpty()) {
986 TObjArrayIter presentUnits(node->InnerOf());
987 while ( (unit = (AliITSneuron*)presentUnits.Next()) ) {
988 if (unit->Inner() == node && unit->Outer() == match) {
995 AliITSneuron *unit = new AliITSneuron(node, match, fEdge2, fEdge1);
996 fNeurons->AddLast(unit);
997 node->InnerOf()->AddLast(unit);
998 match->OuterOf()->AddLast(unit);
999 made += CreateNeurons(match, curvStep, node->X(), node->Y(), node->Z());
1003 // Of course, the return value contains the number of neurons
1004 // counting in also the oned created in all levels of recursive calls.
1008 //__________________________________________________________________________________
1009 Int_t AliITStrackerANN::CreateNetwork(Int_t sector, Int_t curvStep)
1011 // This function simply recalls the CreateNeurons() method for each node
1012 // in the first layer, for the current sector.
1013 // This generates the whole network, thanks to the recursive calls.
1016 // 1) current sector
1017 // 2) current curvature step
1020 // - scans the nodes array for all theta's in the current sector
1021 // and layer 0, and calls the CreateNeurons() function.
1023 // removes all eventually present neurons
1024 if (fNeurons) delete fNeurons;
1025 fNeurons = new TObjArray;
1026 fNeurons->SetOwner(kTRUE);
1028 // calls the ResetNodes() function to free the AliITSnode arrays
1029 if (fMsgLevel >= 2) {
1030 cout << "Sector " << sector << " PHI = ";
1031 cout << fSectorWidth * (Double_t)sector << " --> ";
1032 cout << fSectorWidth * (Double_t)(sector + 1) << endl;
1033 cout << "Curvature step " << curvStep << " [cut = " << fCurvCut[curvStep] << "]" << endl;
1037 // meaningful counters
1038 Int_t itheta, neurons = 0;
1039 TObjArray *lstSector = 0;
1042 Double_t vx[6], vy[6], vz[6];
1043 AliITSnode *p[6] = {0, 0, 0, 0, 0, 0};
1044 for (itheta = 0; itheta < 180; itheta++) {
1045 lstSector = (TObjArray*)fNodes[0][itheta]->At(sector);
1046 TObjArrayIter lay0(lstSector);
1047 while ( (p[0] = (AliITSnode*)lay0.Next()) ) {
1048 if (p[0]->IsUsed()) continue;
1052 neurons += CreateNeurons(p[0], curvStep, fVertexX, fVertexY, fVertexZ);
1054 TObjArrayIter lay1(p[0]->Matches());
1055 while ( (p[1] = (AliITSnode*)lay1.Next()) ) {
1056 if (p[1]->IsUsed()) continue;
1057 if (!PassCurvCut(p[0], p[1], curvStep, vx[0], vy[0], vz[0])) continue;
1058 unit = new AliITSneuron;
1059 unit->Inner() = p[0];
1060 unit->Outer() = p[1];
1061 unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1062 unit->fGain = new TObjArray;
1063 fNeurons->AddLast(unit);
1064 p[0]->InnerOf()->AddLast(unit);
1065 p[1]->OuterOf()->AddLast(unit);
1070 TObjArrayIter lay2(p[1]->Matches());
1071 while ( (p[2] = (AliITSnode*)lay2.Next()) ) {
1072 if (p[2]->IsUsed()) continue;
1073 if (!PassCurvCut(p[1], p[2], curvStep, vx[1], vy[1], vz[1])) continue;
1074 unit = new AliITSneuron;
1075 unit->Inner() = p[1];
1076 unit->Outer() = p[2];
1077 unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1078 unit->fGain = new TObjArray;
1079 fNeurons->AddLast(unit);
1080 p[1]->InnerOf()->AddLast(unit);
1081 p[2]->OuterOf()->AddLast(unit);
1086 TObjArrayIter lay3(p[2]->Matches());
1087 while ( (p[3] = (AliITSnode*)lay3.Next()) ) {
1088 if (p[3]->IsUsed()) continue;
1089 if (!PassCurvCut(p[2], p[3], curvStep, vx[2], vy[2], vz[2])) continue;
1090 unit = new AliITSneuron;
1091 unit->Inner() = p[2];
1092 unit->Outer() = p[3];
1093 unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1094 unit->fGain = new TObjArray;
1095 fNeurons->AddLast(unit);
1096 p[2]->InnerOf()->AddLast(unit);
1097 p[3]->OuterOf()->AddLast(unit);
1102 TObjArrayIter lay4(p[3]->Matches());
1103 while ( (p[4] = (AliITSnode*)lay4.Next()) ) {
1104 if (p[4]->IsUsed()) continue;
1105 if (!PassCurvCut(p[3], p[4], curvStep, vx[3], vy[3], vz[3])) continue;
1106 unit = new AliITSneuron;
1107 unit->Inner() = p[3];
1108 unit->Outer() = p[4];
1109 unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1110 unit->fGain = new TObjArray;
1111 fNeurons->AddLast(unit);
1112 p[3]->InnerOf()->AddLast(unit);
1113 p[4]->OuterOf()->AddLast(unit);
1118 TObjArrayIter lay5(p[4]->Matches());
1119 while ( (p[5] = (AliITSnode*)lay5.Next()) ) {
1120 if (p[5]->IsUsed()) continue;
1121 if (!PassCurvCut(p[4], p[5], curvStep, vx[4], vy[4], vz[4])) continue;
1122 unit = new AliITSneuron;
1123 unit->Inner() = p[4];
1124 unit->Outer() = p[5];
1125 unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1126 unit->fGain = new TObjArray;
1127 fNeurons->AddLast(unit);
1128 p[4]->InnerOf()->AddLast(unit);
1129 p[5]->OuterOf()->AddLast(unit);
1138 } // for (itheta...)
1139 // END OF NEW VERSION
1142 for (ilayer = 0; ilayer < 6; ilayer++) {
1143 for (itheta = 0; itheta < 180; itheta++) {
1144 lstSector = (TObjArray*)fNodes[ilayer][itheta]->At(sector_idx);
1145 TObjArrayIter inners(lstSector);
1146 while ( (inner = (AliITSnode*)inners.Next()) ) {
1147 if (inner->GetUser() >= 0) continue;
1148 TObjArrayIter outers(inner->Matches());
1149 while ( (outer = (AliITSnode*)outers.Next()) ) {
1150 if (outer->GetUser() >= 0) continue;
1151 if (!PassCurvCut(inner, outer, curvStep, fVX, fVY, fVZ)) continue;
1152 unit = new AliITSneuron;
1153 unit->Inner() = inner;
1154 unit->Outer() = outer;
1155 unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1156 unit->fGain = new TObjArray;
1157 fNeurons->AddLast(unit);
1158 inner->InnerOf()->AddLast(unit);
1159 outer->OuterOf()->AddLast(unit);
1163 } // for (itheta...)
1164 } // for (ilayer...)
1167 fNeurons->SetOwner();
1171 //__________________________________________________________________________________
1172 Int_t AliITStrackerANN::LinkNeurons()
1174 /***********************************************************************************
1178 Scans the whole neuron array, in order to find all neuron pairs
1179 which are connected in sequence and share a positive weight.
1180 For each of them, an AliITSlink is created, which stores
1181 the weight value, and will allow for a faster calculation
1182 of the total neural input for each updating cycle.
1184 Every neuron contains an object array which stores all AliITSlink
1185 objects which point to sequenced units, with the respective weights.
1188 - the number of link created for the neural network.
1189 If they are 0, no updating can be done and the step is skipped.
1191 ***********************************************************************************/
1193 // meaningful indexes
1195 Double_t weight = 0.0;
1196 TObjArrayIter neurons(fNeurons), *iter;
1197 AliITSneuron *neuron = 0, *test = 0;
1199 // scan in the neuron array
1201 neuron = (AliITSneuron*)neurons.Next();
1203 // checks for neurons startin from its head ( -> )
1204 iter = (TObjArrayIter*)neuron->Inner()->OuterOf()->MakeIterator();
1206 test = (AliITSneuron*)iter->Next();
1208 weight = Weight(test, neuron);
1209 if (weight > 0.0) neuron->Gain()->AddLast(new AliITSlink(weight, test));
1213 // checks for neurons ending in its tail ( >- )
1214 iter = (TObjArrayIter*)neuron->Outer()->InnerOf()->MakeIterator();
1216 test = (AliITSneuron*)iter->Next();
1218 weight = Weight(neuron, test);
1219 if (weight > 0.0) neuron->Gain()->AddLast(new AliITSlink(weight, test));
1227 //__________________________________________________________________________________
1228 Bool_t AliITStrackerANN::Update()
1230 /***********************************************************************************
1232 Performs a single updating cycle.
1235 - for each neuron, gets the activation with the neuron Activate() method
1236 - checks if stability has been reached (compare mean activation variation
1237 with the stability threshold data member)
1240 - kTRUE means that the neural network has stabilized
1241 - kFALSE means that another updating cycle is needed
1243 ***********************************************************************************/
1245 Double_t actVar = 0.0, totDiff = 0.0;
1246 TObjArrayIter iter(fNeurons);
1249 unit = (AliITSneuron*)iter.Next();
1251 actVar = unit->Activate(fTemperature);
1252 // calculation the relative activation variation
1255 totDiff /= fNeurons->GetSize();
1256 return (totDiff < fStabThreshold);
1259 //__________________________________________________________________________________
1260 void AliITStrackerANN::FollowChains(Int_t sector)
1262 /***********************************************************************************
1266 After that the neural network has stabilized,
1267 the final step is to create polygonal chains
1268 of clusters, one in each layer, which represent
1269 the tracks recognized by the neural algorithm.
1270 This is made by means of a choice of the best
1271 neuron among the ones starting from each point.
1273 Once that such neuron is selected, its inner point
1274 will set the 'fNext' field to its outer point, and
1275 similarly, its outer point will set the 'fPrev' field
1277 This defines a bi-directional sequence.
1279 In this procedure, it can happen that many neurons
1280 which have the head of the arrow in a given node, will
1281 all select as best following the neuron with the largest
1282 activation starting in that point.
1283 This results in MANY nodes which have the same 'fNext'.
1284 But, this field will be set to NULL for all these points,
1285 but the only one which is pointed by the 'fPrev' field
1286 of this shared node.
1288 ***********************************************************************************/
1290 // meaningful counters
1291 Int_t itheta, ilayer;
1292 TObjArray *lstSector = 0;
1293 Double_t test = fActMinimum;
1295 AliITSneuron *n = 0;
1297 // scan the whole collection of nodes
1298 for (ilayer = 0; ilayer < fNLayers; ilayer++) {
1299 for (itheta = 0; itheta < 180; itheta++) {
1300 // get the array of point in a given layer/theta-slot/sector
1301 lstSector = (TObjArray*)fNodes[ilayer][itheta]->At(sector);
1302 TObjArrayIter nodes(lstSector);
1303 while ( (p = (AliITSnode*)nodes.Next()) ) {
1304 // if the point is used, it is skipped
1305 if (p->IsUsed()) continue;
1306 // initially, fNext points to nothing, and
1307 // the comparison value is set to the minimum activation
1308 // which allows to say that a neuron is turned 'on'
1309 // a node from which only 'off' neurons start is probably
1310 // a noise point, which will be excluded from the reconstruction.
1313 TObjArrayIter innerof(p->InnerOf());
1314 while ( (n = (AliITSneuron*)innerof.Next()) ) {
1315 // if the examined neuron has not the largest activation
1316 // it is skipped and removed from array of all neurons
1317 // and of its outer point (its inner is the cursor p)
1318 if (n->Activation() < test) {
1319 p->InnerOf()->Remove(n);
1320 n->Outer()->OuterOf()->Remove(n);
1321 delete fNeurons->Remove(n);
1324 // otherwise, its activation becomes the maximum reference
1325 p->Next() = n->Outer();
1326 // at the exit of the while(), the fNext will point
1327 // to the outer node of the neuron starting in p, whose
1328 // activation is the largest.
1330 // the same procedure is made now for all neurons
1331 // for which p is the outer point
1334 TObjArrayIter outerof(p->OuterOf());
1335 while ( (n = (AliITSneuron*)outerof.Next()) ) {
1336 // if the examined neuron has not the largest activation
1337 // it is skipped and removed from array of all neurons
1338 // and of its inner point (its outer is the cursor p)
1339 if (n->Activation() < test) {
1340 p->OuterOf()->Remove(n);
1341 n->Inner()->InnerOf()->Remove(n);
1342 delete fNeurons->Remove(n);
1345 // otherwise, its activation becomes the maximum reference
1346 p->Prev() = n->Inner();
1347 // at the exit of the while(), the fPrev will point
1348 // to the inner node of the neuron ending in p, whose
1349 // activation is the largest.
1351 } // end while (p ...)
1352 } // end for (itheta ...)
1353 } // end for (ilayer ...)
1355 // now the mismatches are solved
1356 Bool_t matchPrev, matchNext;
1357 for (ilayer = 0; ilayer < fNLayers; ilayer++) {
1358 for (itheta = 0; itheta < 180; itheta++) {
1359 // get the array of point in a given layer/theta-slot/sector
1360 lstSector = (TObjArray*)fNodes[ilayer][itheta]->At(sector);
1361 TObjArrayIter nodes(lstSector);
1362 while ( (p = (AliITSnode*)nodes.Next()) ) {
1363 // now p will point to a fPrev and a fNext node.
1364 // Ideally they are placed this way: fPrev --> P --> fNext
1365 // A mismatch happens if the point addressed as fPrev does NOT
1366 // point to p as its fNext. And the same for the point addressed
1368 // In this case, the fNext and fPrev pointers are set to NULL
1369 // and p is excluded from the reconstruction
1370 matchPrev = matchNext= kFALSE;
1371 if (ilayer > 0 && p->Prev() != NULL)
1372 if (p->Prev()->Next() == p) matchPrev = kTRUE;
1373 if (ilayer < 5 && p->Next() != NULL)
1374 if (p->Next()->Prev() == p) matchNext = kTRUE;
1377 else if (ilayer == 5)
1379 if (!matchNext || !matchPrev) {
1380 p->Prev() = p->Next() = 0;
1382 } // end while (p ...)
1383 } // end for (itheta ...)
1384 } // end for (ilayer ...)
1387 //__________________________________________________________________________________
1388 Int_t AliITStrackerANN::SaveTracks(Int_t sector)
1390 /********************************************************************************
1394 Using the fNext and fPrev pointers, the chain is followed
1395 and the track is fitted and saved.
1396 Of course, the track is followed as a chain with a point
1397 for each layer, then the track following starts always
1398 from the clusters in layer 0.
1400 ***********************************************************************************/
1402 // if not initialized, the tracks TobjArray is initialized
1403 if (!fFoundTracks) fFoundTracks = new TObjArray;
1405 // meaningful counters
1406 Int_t itheta, ilayer, l;
1407 TObjArray *lstSector = 0;
1408 AliITSnode *p = 0, *q = 0, **node = new AliITSnode*[fNLayers];
1409 for (l = 0; l < fNLayers; l++) node[l] = 0;
1412 array = new AliITSnode*[fNLayers + 1];
1413 for (l = 0; l <= fNLayers; l++) array[l] = 0;
1414 array[0] = new AliITSnode();
1415 array[0]->X() = fVertexX;
1416 array[0]->Y() = fVertexY;
1417 array[0]->Z() = fVertexZ;
1418 array[0]->ErrX2() = fVertexErrorX;
1419 array[0]->ErrY2() = fVertexErrorY;
1420 array[0]->ErrZ2() = fVertexErrorZ;
1422 Double_t *param = new Double_t[8];
1424 // scan the whole collection of nodes
1425 for (ilayer = 0; ilayer < 1; ilayer++) {
1426 for (itheta = 0; itheta < 180; itheta++) {
1427 // get the array of point in a given layer/theta-slot/sector
1428 lstSector = (TObjArray*)fNodes[ilayer][itheta]->At(sector);
1429 TObjArrayIter nodes(lstSector);
1430 while ( (p = (AliITSnode*)nodes.Next()) ) {
1431 for (q = p; q; q = q->Next()) {
1435 //if (!RiemannFit(fNLayers, node, param)) continue;
1436 // initialization of Kalman Filter Tracking
1437 AliITSRecPoint *cluster = (AliITSRecPoint*)GetCluster(node[0]->ClusterRef());
1438 Int_t mod = cluster->GetDetectorIndex();
1439 Int_t lay, lad, det;
1440 fGeom->GetModuleId(mod, lay, lad, det);
1441 Float_t y0 = cluster->GetY();
1442 Float_t z0 = cluster->GetZ();
1443 AliITStrackSA* trac = new AliITStrackSA(fGeom,lay, lad, det,
1445 param[4], param[7], param[3], 1);
1446 for (l = 0; l < fNLayers; l++) {
1447 cluster = (AliITSRecPoint*)GetCluster(node[l]->ClusterRef());
1448 if (cluster) trac->AddClusterV2(l, (node[l]->ClusterRef() & 0x0fffffff)>>0);
1450 AliITStrackV2* ot = new AliITStrackV2(*trac);
1451 ot->ResetCovariance();
1452 ot->ResetClusters();
1453 if (RefitAt(49.,ot,trac)) { //fit from layer 1 to layer 6
1454 AliITStrackV2 *otrack2 = new AliITStrackV2(*ot);
1455 otrack2->ResetCovariance();
1456 otrack2->ResetClusters();
1457 //fit from layer 6 to layer 1
1458 if (RefitAt(3.7,otrack2,ot)) fFoundTracks->AddLast(otrack2);
1460 // end of Kalman Filter fit
1468 //__________________________________________________________________________________
1469 void AliITStrackerANN::ExportTracks(const char *filename) const
1471 // Exports found tracks into a TTree of AliITStrackV2 objects
1472 TFile *file = new TFile(filename, "RECREATE");
1473 TTree *tree = new TTree("TreeT-ANN", "Tracks found in ITS stand-alone with Neural Tracking");
1474 AliITStrackV2 *track = 0;
1475 tree->Branch("Tracks", &track, "AliITStrackV2");
1476 TObjArrayIter tracks(fFoundTracks);
1477 while ( (track = (AliITStrackV2*)tracks.Next()) ) {
1486 //__________________________________________________________________________________
1487 void AliITStrackerANN::CleanNetwork()
1489 // Removes deactivated units from the network
1491 AliITSneuron *unit = 0;
1492 TObjArrayIter neurons(fNeurons);
1493 while ( (unit = (AliITSneuron*)neurons.Next()) ) {
1494 if (unit->Activation() < fActMinimum) {
1495 unit->Inner()->InnerOf()->Remove(unit);
1496 unit->Outer()->OuterOf()->Remove(unit);
1497 delete fNeurons->Remove(unit);
1503 AliITSneuron *enemy = 0;
1505 while ( (unit = (AliITSneuron*)neurons.Next()) ) {
1506 nIn = (Int_t)unit->Inner()->InnerOf()->GetSize();
1507 nOut = (Int_t)unit->Outer()->OuterOf()->GetSize();
1508 if (nIn < 2 && nOut < 2) continue;
1511 TObjArrayIter competing(unit->Inner()->InnerOf());
1512 while ( (enemy = (AliITSneuron*)competing.Next()) ) {
1513 if (unit->Activation() > enemy->Activation()) {
1514 enemy->Inner()->InnerOf()->Remove(enemy);
1515 enemy->Outer()->OuterOf()->Remove(enemy);
1516 delete fNeurons->Remove(enemy);
1519 unit->Inner()->InnerOf()->Remove(unit);
1520 unit->Outer()->OuterOf()->Remove(unit);
1521 delete fNeurons->Remove(unit);
1526 if (removed) continue;
1529 TObjArrayIter competing(unit->Outer()->OuterOf());
1530 while ( (enemy = (AliITSneuron*)competing.Next()) ) {
1531 if (unit->Activation() > enemy->Activation()) {
1532 enemy->Inner()->InnerOf()->Remove(enemy);
1533 enemy->Outer()->OuterOf()->Remove(enemy);
1534 delete fNeurons->Remove(enemy);
1537 unit->Inner()->InnerOf()->Remove(unit);
1538 unit->Outer()->OuterOf()->Remove(unit);
1539 delete fNeurons->Remove(unit);
1548 //__________________________________________________________________________________
1549 Int_t AliITStrackerANN::StoreTracks()
1551 // Stores the tracks found in a single neural tracking step.
1552 // In order to do this, it sects each neuron which has a point
1553 // in the first layer.
1556 // if not initialized, the tracks TobjArray is initialized
1557 if (!fFoundTracks) fFoundTracks = new TObjArray;
1559 Int_t i, check, stored = 0;
1560 Double_t testAct = 0;
1561 AliITSneuron *unit = 0, *cursor = 0, *fwd = 0;
1562 AliITSnode *node = 0;
1563 TObjArrayIter iter(fNeurons), *fwdIter;
1564 TObjArray *removedUnits = new TObjArray(0);
1565 removedUnits->SetOwner(kFALSE);
1566 AliITStrackANN annTrack(fNLayers);
1569 unit = (AliITSneuron*)iter.Next();
1571 if (unit->Inner()->GetLayer() > 0) continue;
1572 annTrack.SetNode(unit->Inner()->GetLayer(), unit->Inner());
1573 annTrack.SetNode(unit->Outer()->GetLayer(), unit->Outer());
1574 node = unit->Outer();
1575 removedUnits->AddLast(unit);
1577 testAct = fActMinimum;
1578 fwdIter = (TObjArrayIter*)node->InnerOf()->MakeIterator();
1581 cursor = (AliITSneuron*)fwdIter->Next();
1583 if (cursor->Used()) continue;
1584 if (cursor->Activation() >= testAct) {
1585 testAct = cursor->Activation();
1590 removedUnits->AddLast(fwd);
1591 node = fwd->Outer();
1592 annTrack.SetNode(node->GetLayer(), node);
1594 check = annTrack.CheckOccupation();
1598 //if (!RiemannFit(fNLayers, trackitem, param)) continue;
1599 if (!annTrack.RiemannFit()) continue;
1600 // initialization of Kalman Filter Tracking
1601 AliITSRecPoint *cluster = (AliITSRecPoint*)GetCluster(annTrack[0]->ClusterRef());
1602 Int_t mod = cluster->GetDetectorIndex();
1603 Int_t lay, lad, det;
1604 fGeom->GetModuleId(mod, lay, lad, det);
1605 Float_t y0 = cluster->GetY();
1606 Float_t z0 = cluster->GetZ();
1607 AliITStrackSA* trac = new AliITStrackSA(fGeom,lay, lad, det, y0, z0,
1608 annTrack.Phi(), annTrack.TanLambda(),
1609 annTrack.Curv(), 1);
1610 for (Int_t l = 0; l < fNLayers; l++) {
1611 if (!annTrack[l]) continue;
1612 cluster = (AliITSRecPoint*)GetCluster(annTrack[l]->ClusterRef());
1613 if (cluster) trac->AddClusterV2(l, (annTrack[l]->ClusterRef() & 0x0fffffff)>>0);
1615 AliITStrackV2* ot = new AliITStrackV2(*trac);
1616 ot->ResetCovariance();
1617 ot->ResetClusters();
1618 if (RefitAt(49.,ot,trac)) { //fit from layer 1 to layer 6
1619 AliITStrackV2 *otrack2 = new AliITStrackV2(*ot);
1620 otrack2->ResetCovariance();
1621 otrack2->ResetClusters();
1622 //fit from layer 6 to layer 1
1623 if (RefitAt(3.7,otrack2,ot)) fFoundTracks->AddLast(otrack2);
1625 // end of Kalman Filter fit
1627 for (i = 0; i < fNLayers; i++) {
1628 //node = (AliITSnode*)removedPoints->At(i);
1632 fwdIter = (TObjArrayIter*)removedUnits->MakeIterator();
1634 cursor = (AliITSneuron*)fwdIter->Next();
1644 Double_t AliITStrackerANN::Weight(AliITSneuron *nAB, AliITSneuron *nBC)
1646 /***********************************************************************************
1647 * Calculation of neural weight.
1648 * The implementation of positive neural weight is set only in the case
1649 * of connected units (e.g.: A->B with B->C).
1650 * Given that B is the **common** point. care should be taken to pass
1651 * as FIRST argument the neuron going "to" B, and
1652 * as SECOND argument the neuron starting "from" B
1653 * anyway, a check is put in order to return 0.0 when arguments are not well given.
1654 ***********************************************************************************/
1656 if (nAB->Outer() != nBC->Inner()) {
1657 if (nBC->Outer() == nAB->Inner()) {
1658 AliITSneuron *temp = nAB;
1662 if (fMsgLevel >= 3) {
1663 Info("Weight", "Switching wrongly ordered arguments.");
1666 Warning("Weight", "Not connected segments. Returning 0.0");
1670 AliITSnode *pA = nAB->Inner();
1671 AliITSnode *pB = nAB->Outer();
1672 AliITSnode *pC = nBC->Outer();
1674 TVector3 vAB(pB->X() - pA->X(), pB->Y() - pA->Y(), pB->Z() - pA->Z());
1675 TVector3 vBC(pC->X() - pB->X(), pC->Y() - pB->Y(), pC->Z() - pB->Z());
1677 Double_t weight = 1.0 - sin(vAB.Angle(vBC));
1678 return fGain2CostRatio * TMath::Power(weight, fExponent);
1683 /******************************************
1684 ******************************************
1685 *** AliITStrackerANN::AliITSnode class ***
1686 ******************************************
1687 ******************************************/
1689 //__________________________________________________________________________________
1690 AliITStrackerANN::AliITSnode::AliITSnode()
1691 : fUsed(kFALSE), fClusterRef(-1),
1692 fMatches(NULL), fInnerOf(NULL), fOuterOf(NULL),
1693 fNext(NULL), fPrev(NULL)
1695 // Constructor for the embedded 'AliITSnode' class.
1696 // It initializes all pointer-like objects.
1699 fEX2 = fEY2 = fEZ2 = 0.0;
1702 //__________________________________________________________________________________
1703 AliITStrackerANN::AliITSnode::~AliITSnode()
1705 // Destructor for the embedded 'AliITSnode' class.
1706 // It should clear the object arrays, but it is possible
1707 // that some objects still are useful after the point deletion
1708 // then the arrays are cleared but their objects are owed by
1709 // another TCollection object, and not deleted.
1710 // For safety reasons, all the pointers are set to zero.
1712 fInnerOf->SetOwner(kFALSE);
1716 fOuterOf->SetOwner(kFALSE);
1720 fMatches->SetOwner(kFALSE);
1728 //__________________________________________________________________________________
1729 Double_t AliITStrackerANN::AliITSnode::GetPhi() const
1731 // Calculates the 'phi' (azimutal) angle, and returns it
1732 // in the range between 0 and 2Pi radians.
1735 q = TMath::ATan2(fY,fX);
1739 return q + TMath::TwoPi();
1742 //__________________________________________________________________________________
1743 Double_t AliITStrackerANN::AliITSnode::GetError(Option_t *option)
1745 // Returns the error or the square error of
1746 // values related to the coordinates in different systems.
1747 // The option argument specifies the coordinate error desired:
1749 // "R2" --> error in transverse radius
1750 // "R3" --> error in spherical radius
1751 // "PHI" --> error in azimuthal angle
1752 // "THETA" --> error in polar angle
1753 // "SQ" --> get the square of error
1755 // In order to get the error on the cartesian coordinates
1756 // reference to the inline ErrX2(), ErrY2() adn ErrZ2() methods.
1758 TString opt(option);
1759 Double_t errorSq = 0.0;
1762 if (opt.Contains("R2")) {
1763 errorSq = fX*fX*fEX2 + fY*fY*fEY2;
1764 errorSq /= GetR2sq();
1766 else if (opt.Contains("R3")) {
1767 errorSq = fX*fX*fEX2 + fY*fY*fEY2 + fZ*fZ*fEZ2;
1768 errorSq /= GetR3sq();
1770 else if (opt.Contains("PHI")) {
1771 errorSq = fY*fY*fEX2;
1772 errorSq += fX*fX*fEY2;
1773 errorSq /= GetR2sq() * GetR2sq();
1775 else if (opt.Contains("THETA")) {
1776 errorSq = fZ*fZ * (fX*fX*fEX2 + fY*fY*fEY2);
1777 errorSq += GetR2sq() * GetR2sq() * fEZ2;
1778 errorSq /= GetR3sq() * GetR3sq() * GetR2() * GetR2();
1781 if (!opt.Contains("SQ"))
1782 return TMath::Sqrt(errorSq);
1789 /********************************************
1790 ********************************************
1791 *** AliITStrackerANN::AliITSneuron class ***
1792 ********************************************
1793 ********************************************/
1795 //__________________________________________________________________________________
1796 AliITStrackerANN::AliITSneuron::AliITSneuron
1797 (AliITSnode *inner, AliITSnode *outer, Double_t minAct, Double_t maxAct)
1798 : fUsed(0), fInner(inner), fOuter(outer)
1800 // Default neuron constructor
1801 fActivation = gRandom->Rndm() * (maxAct-minAct) + minAct;
1802 fGain = new TObjArray;
1805 //__________________________________________________________________________________
1806 Double_t AliITStrackerANN::AliITSneuron::Activate(Double_t temperature)
1808 // This computes the new activation of a neuron, and returns
1809 // its activation variation as a consequence of the updating.
1812 // - the 'temperature' parameter for the neural activation logistic function
1815 // - collects the total gain, by summing the products
1816 // of the activation of each sequenced unit by the relative weight.
1817 // - collects the total cost, by summing the activations of
1818 // all competing units
1819 // - passes the sum of gain - cost to the activation function and
1820 // calculates the new activation
1823 // - the difference between the old activation and the new one
1827 Double_t sumGain = 0.0; // total contribution from chained neurons
1828 Double_t sumCost = 0.0; // total contribution from crossing neurons
1829 Double_t input; // total input
1830 Double_t actOld, actNew; // old and new values for the activation
1831 AliITSneuron *linked = 0; // cursor for scanning the neuron arrays (for link check)
1832 AliITSlink *link; // cursor for scanning the synapses arrays (for link check)
1833 TObjArrayIter *iterator = 0; // pointer to the iterator along the neuron arrays
1835 // sum contributions from the correlated units
1836 iterator = (TObjArrayIter*)fGain->MakeIterator();
1838 link = (AliITSlink*)iterator->Next();
1840 sumGain += link->Contribution();
1844 // sum contributions from the competing units:
1845 // the ones which have the same starting point...
1846 iterator = (TObjArrayIter*)fInner->InnerOf()->MakeIterator();
1848 linked = (AliITSneuron*)iterator->Next();
1850 if (linked == this) continue;
1851 sumCost += linked->fActivation;
1854 // ...and the ones which have the same ending point
1855 iterator = (TObjArrayIter*)fOuter->OuterOf()->MakeIterator();
1857 linked = (AliITSneuron*)iterator->Next();
1859 if (linked == this) continue;
1860 sumCost += linked->fActivation;
1863 // calculate the total input as the difference between gain and cost
1864 input = (sumGain - sumCost) / temperature;
1865 actOld = fActivation;
1866 // calculate the final output
1867 #ifdef NEURAL_LINEAR
1868 if (input <= -2.0 * temperature)
1870 else if (input >= 2.0 * temperature)
1873 actNew = input / (4.0 * temperature) + 0.5;
1875 actNew = 1.0 / (1.0 + TMath::Exp(-input));
1877 fActivation = actNew;
1879 // return the activation variation
1880 return TMath::Abs(actNew - actOld);
1885 /******************************************
1886 ******************************************
1887 *** AliITStrackerANN::AliITSlink class ***
1888 ******************************************
1889 ******************************************/
1891 // No methods defined non-inline
1895 /**********************************************
1896 **********************************************
1897 *** AliITStrackerANN::AliITStrackANN class ***
1898 **********************************************
1899 **********************************************/
1901 //__________________________________________________________________________________
1902 AliITStrackerANN::AliITStrackANN::AliITStrackANN(Int_t dim) : fNPoints(dim)
1904 // Default constructor for the AliITStrackANN class
1919 fNode = new AliITSnode*[dim];
1920 for (i = 0; i < dim; i++) fNode[i] = 0;
1924 //__________________________________________________________________________________
1925 Int_t AliITStrackerANN::AliITStrackANN::CheckOccupation() const
1927 // Returns the number of pointers fNode which are not NULL
1930 Int_t count = 0; // counter for how many points are stored in the track
1932 for (i = 0; i < fNPoints; i++) {
1933 if (fNode[i] != NULL) count++;
1939 //__________________________________________________________________________________
1940 Bool_t AliITStrackerANN::AliITStrackANN::RiemannFit()
1942 // Performs the Riemann Sphere fit for the given points to a circle
1943 // and then uses the fit parameters to fit a helix in space.
1946 // - kTRUE if all operations have been performed
1947 // - kFALSE if the numbers risk to lead to an arithmetic violation
1949 Int_t i, j, count, dim = fNPoints;
1951 // First check for all points
1952 count = CheckOccupation();
1953 if (count != fNPoints) {
1954 Error ("AliITStrackANN::RiemannFit", "CheckOccupations returns %d, fNPoints = %d ==> MISMATCH", count, fNPoints);
1960 for (i = 0; i < dim; i++) m1(i,0) = 1.0;
1962 // matrix of Rieman projection coordinates
1963 TMatrixD coords(dim,3);
1964 for (i = 0; i < dim; i++) {
1965 coords(i,0) = fNode[i]->X();
1966 coords(i,1) = fNode[i]->Y();
1967 coords(i,2) = fNode[i]->GetR2sq();
1970 // matrix of weights
1971 Double_t xterm, yterm, ex, ey;
1972 TMatrixD weights(dim,dim);
1973 for (i = 0; i < dim; i++) {
1974 xterm = fNode[i]->X() * fNode[i]->GetPhi() - fNode[i]->Y() / fNode[i]->GetR2();
1975 ex = fNode[i]->ErrX2();
1976 yterm = fNode[i]->Y() * fNode[i]->GetPhi() + fNode[i]->X() / fNode[i]->GetR2();
1977 ey = fNode[i]->ErrY2();
1978 weights(i,i) = fNode[i]->GetR2sq() / (xterm * xterm * ex + yterm * yterm * ey );
1981 // weighted sample mean
1982 Double_t meanX = 0.0, meanY = 0.0, meanW = 0.0, sw = 0.0;
1983 for (i = 0; i < dim; i++) {
1984 meanX += weights(i,i) * coords(i,0);
1985 meanY += weights(i,i) * coords(i,1);
1986 meanW += weights(i,i) * coords(i,2);
1993 // sample covariance matrix
1994 for (i = 0; i < dim; i++) {
1995 coords(i,0) -= meanX;
1996 coords(i,1) -= meanY;
1997 coords(i,2) -= meanW;
1999 TMatrixD coordsT(TMatrixD::kTransposed, coords);
2000 TMatrixD weights4coords(weights, TMatrixD::kMult, coords);
2001 TMatrixD sampleCov(coordsT, TMatrixD::kMult, weights4coords);
2002 for (i = 0; i < 3; i++) {
2003 for (j = i + 1; j < 3; j++) {
2004 sampleCov(i,j) = sampleCov(j,i) = (sampleCov(i,j) + sampleCov(j,i)) * 0.5;
2008 // Eigenvalue problem solving for V matrix
2010 TVectorD eval(3), n(3);
2011 #if ROOT_VERSION_CODE < ROOT_VERSION(4,0,2)
2012 TMatrixD evec = sampleCov.EigenVectors(eval);
2014 TMatrixDEigen ei(sampleCov);
2015 TMatrixD evec = ei.GetEigenVectors();
2016 eval = ei.GetEigenValuesRe();
2018 if (eval(1) < eval(ileast)) ileast = 1;
2019 if (eval(2) < eval(ileast)) ileast = 2;
2020 n(0) = evec(0, ileast);
2021 n(1) = evec(1, ileast);
2022 n(2) = evec(2, ileast);
2024 // c - known term in the plane intersection with Riemann axes
2025 Double_t c = -(meanX * n(0) + meanY * n(1) + meanW * n(2));
2027 // center and radius of fitted circle
2028 Double_t xc, yc, radius, curv;
2029 xc = -n(0) / (2. * n(2));
2030 yc = -n(1) / (2. * n(2));
2031 radius = (1. - n(2)*n(2) - 4.*c*n(2)) / (4. * n(2) * n(2));
2033 if (radius <= 0.E0) {
2034 Error("RiemannFit", "Radius = %f less than zero!!!", radius);
2037 radius = TMath::Sqrt(radius);
2038 curv = 1.0 / radius;
2040 // evaluating signs for curvature and others
2041 Double_t phi1 = 0.0, phi2, temp1, temp2, phi0, sumdphi = 0.0;
2042 AliITSnode *p = fNode[0];
2044 for (i = 1; i < dim; i++) {
2045 p = (AliITSnode*)fNode[i];
2050 if (temp1 > fgkPi && temp2 < fgkPi)
2052 else if (temp1 < fgkPi && temp2 > fgkPi)
2054 sumdphi += temp2 - temp1;
2057 if (sumdphi < 0.E0) curv = -curv;
2058 Double_t diff, angle = TMath::ATan2(yc, xc);
2060 phi0 = angle + 0.5 * TMath::Pi();
2062 phi0 = angle - 0.5 * TMath::Pi();
2063 diff = angle - phi0;
2065 Double_t dt, temp = TMath::Sqrt(xc*xc + yc*yc) - radius;
2070 //cout << "Dt = " << dt << endl;
2072 Double_t halfC = 0.5 * curv, test;
2073 Double_t *s = new Double_t[dim], *zz = new Double_t[dim], *ws = new Double_t[dim];
2074 for (j = 0; j < 6; j++) {
2078 s[j] = (p->GetR2sq() - dt * dt) / (1. + curv * dt);
2080 if (TMath::Abs(s[j]) < 1.E-6) s[j] = 0.;
2082 Error("RiemannFit", "Square root argument error: %17.15g < 0", s[j]);
2086 s[j] = TMath::Sqrt(s[j]);
2087 //cout << "Curv = " << halfC << " --- s[" << j << "] = " << s[j] << endl;
2089 test = TMath::Abs(s[j]);
2092 s[j] = ((s[j] > 0.) ? 0.99999999999 : -0.9999999999);
2094 Error("RiemannFit", "Value too large: %17.15g", s[j]);
2100 s[j] = TMath::ASin(s[j]) / halfC;
2101 ws[j] = 1.0 / (p->ErrZ2());
2104 // second tep final fit
2105 Double_t s2Sum = 0.0, zSum = 0.0, szSum = 0.0, sSum = 0.0, sumw = 0.0;
2106 for (i = 0; i < dim; i++) {
2107 s2Sum += ws[i] * s[i] * s[i];
2108 zSum += ws[i] * zz[i];
2109 sSum += ws[i] * s[i];
2110 szSum += ws[i] * s[i] * zz[i];
2117 temp = s2Sum - sSum*sSum;
2120 dz = (s2Sum*zSum - sSum*szSum) / temp;
2121 tanL = (szSum - sSum*zSum) / temp;