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 // ---------------------------------------------------------------------------------
21 #include <Riostream.h>
24 #include <TObjArray.h>
30 #if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,2)
31 #include <TMatrixDEigen.h>
34 #include "AliITSgeom.h"
35 #include "AliITStrackSA.h"
36 #include "AliITSclusterV2.h"
37 #include "AliITStrackV2.h"
39 #include "AliITStrackerANN.h"
41 const Double_t AliITStrackerANN::fgkPi = 3.141592653; // pi
42 const Double_t AliITStrackerANN::fgkHalfPi = 1.570796327; // pi / 2
43 const Double_t AliITStrackerANN::fgkTwoPi = 6.283185307; // 2 * pi
45 ClassImp(AliITStrackerANN)
47 //__________________________________________________________________________________
48 AliITStrackerANN::AliITStrackerANN(const AliITSgeom *geom, Int_t msglev)
49 : AliITStrackerV2(geom), fMsgLevel(msglev)
51 /**************************************************************************
53 CONSTRUCTOR (standard-to-use version)
56 1) the ITS geometry used in the generated event
57 2) the flag for log-messages writing
59 The AliITSgeometry is used along the class,
60 in order to translate the local AliITSclusterV2 coordinates
61 into the Global reference frame, which is necessary for the
62 Neural Network algorithm to operate.
63 In case of serialized use, the log messages should be excluded,
64 in order to save real execution time.
67 - all pointer data members are initialized
68 (if possible, otherwise are set to NULL)
69 - all numeric data members are initialized to
70 values which allow the Neural Network to operate
71 even if nothing is externally set.
73 NOTE: it is possible that tracking an event
74 with these default values results in a non-sense.
76 **************************************************************************/
81 fGeom = (AliITSgeom*)geom;
83 // Initialize the array of first module indexes per layer
84 fNLayers = geom->GetNlayers();
85 fFirstModInLayer = new Int_t[fNLayers + 1];
86 for (i = 0; i < fNLayers; i++) {
87 fFirstModInLayer[i] = fGeom->GetModuleIndex(i + 1, 1, 1);
89 fFirstModInLayer[fNLayers] = geom->GetIndexMax();
91 // initialization: no curvature cut steps
95 // initialization: 4 sectors (one for each quadrant)
97 fSectorWidth = fgkHalfPi;
99 // initialization: theta offset of 20 degrees
100 fPolarInterval = 20.0;
102 // initialization: array structure not defined
103 fStructureOK = kFALSE;
105 // initialization: vertex in the origin
110 // initialization: uninitialized point array
113 // initialization: very large (unuseful) cut values
115 for (ilayer = 0; ilayer < 6; ilayer++) {
116 fThetaCut2D[ilayer] = TMath::Pi();
117 fThetaCut3D[ilayer] = TMath::Pi();
118 fHelixMatchCut[ilayer] = 1.0;
121 // initialization: inictial activation range between 0.3 and 0.7
125 // initialization: neural network operation & weights
127 fStabThreshold = 0.001;
128 fGain2CostRatio = 1.0;
132 // initialization: uninitialized array of neurons
135 // initialization: uninitialized array of tracks
139 //__________________________________________________________________________________
140 void AliITStrackerANN::SetCuts
141 (Int_t ncurv, Double_t *curv, Double_t *theta2D, Double_t *theta3D, Double_t *helix)
143 /**************************************************************************
147 Allows for the definition of all kind of geometric cuts
148 which have been studied in for the creation of a neuron
149 from a pair of clusters C1 and C2 on consecutive layers.
150 Neuron will be created only if the pair passes ALL of these cuts.
152 At the moment, we define 4 kinds of geometrical cuts:
153 a) cut on the difference of the polar 'theta' angle;
154 b) cut on the angle between origin->C1 and C1->C2 in space;
155 c) cut on the curvature of the circle passing
156 through C1, C2 and the primary vertex;
157 d) cut on heli-matching of the same three points.
160 1) the number of curvature cut steps
161 2) the array of curvature cuts for each step
162 (its dimension is given by the argument 1)
163 3) array of 5 cut values (one for each consecutive lauer pair)
165 4) array of 5 cut values (one for each consecutive lauer pair)
167 5) array of 5 cut values (one for each consecutive lauer pair)
171 - gets the values for each cut and stores them in data members
172 - in the case of curvature cuts, the cut array
173 (whose size is not fixed) is allocated
175 NOTE: in case that the user wants to set onyl ONE of the 4 cuts array,
176 he can simply pass NULL arguments for other cuts, and (eventually)
177 a ZERO as the first argument (if curvature cuts have not to be set).
179 Anyway, all the cuts have to be set at least once.
181 **************************************************************************/
186 /*** Curvature cut setting ***/
188 // first of all, the curvature cuts are sorted in increasing order
189 // (from the smallest to the largest one)
190 Int_t *ind = new Int_t[ncurv];
191 TMath::Sort(ncurv, curv, ind, kFALSE);
192 // then, the curvature cut array is allocated and filled
193 // (a message with the list of defined cuts can be optionally shown)
194 fCurvCut = new Double_t[ncurv];
195 if (fMsgLevel >= 1) cout << "Number of curvature cuts: " << ncurv << endl;
196 for (i = 0; i < ncurv; i++) {
197 fCurvCut[i] = curv[ind[i]];
198 if (fMsgLevel >= 1) cout << " - " << fCurvCut[i] << endl;
202 /*** 'Fixed' cuts setting ***/
204 // checks what cuts have to be set
205 Bool_t doTheta2D = (theta2D != 0);
206 Bool_t doTheta3D = (theta3D != 0);
207 Bool_t doHelix = (helix != 0);
208 // sets the cuts for all layer pairs
209 for (i = 0; i < fNLayers - 1; i++) {
210 if (doTheta2D) fThetaCut2D[i] = theta2D[i];
211 if (doTheta3D) fThetaCut3D[i] = theta3D[i];
212 if (doHelix) fHelixMatchCut[i] = helix[i];
214 // if required, lists the cuts
215 if (!fMsgLevel < 2) return;
216 cout << "Theta 2D cuts: ";
218 cout << "<not set>" << endl;
222 for (i = 0; i < fNLayers - 1; i++) {
223 cout << "For layers " << i << " --> " << i + 1;
224 cout << " cut = " << fThetaCut2D[i] << endl;
227 cout << "---" << endl;
228 cout << "Theta 3D cuts: ";
230 cout << "<not set>" << endl;
234 for (i = 0; i < fNLayers - 1; i++) {
235 cout << "For layers " << i << " --> " << i + 1;
236 cout << " cut = " << fThetaCut3D[i] << endl;
239 cout << "---" << endl;
240 cout << "Helix-match cuts: ";
242 cout << "<not set>" << endl;
246 for (i = 0; i < fNLayers - 1; i++) {
247 cout << "For layers " << i << " --> " << i + 1;
248 cout << " cut = " << fHelixMatchCut[i] << endl;
251 cout << "---" << endl;
254 //__________________________________________________________________________________
255 Bool_t AliITStrackerANN::GetGlobalXYZ
257 Double_t &x, Double_t &y, Double_t &z, Double_t &ex2, Double_t &ey2, Double_t &ez2)
259 /**************************************************************************
261 LOCAL TO GLOBAL TRANSLATOR
263 Taking information from the ITS geometry stored in the class,
264 gets a stored AliITScluster and calculates it coordinates
265 and errors in the global reference frame.
266 These values are stored in the variables,
267 which are passed by reference.
270 1) reference index for the cluster to use
271 2) (by reference) place to store the global X-coord into
272 3) (by reference) place to store the global Y-coord into
273 4) (by reference) place to store the global Z-coord into
274 5) (by reference) place to store the global X-coord error into
275 6) (by reference) place to store the global Y-coord error into
276 7) (by reference) place to store the global Z-coord error into
279 essentially, determines the ITS module index from the
280 detector index of the AliITSclusterV2 object, and extracts
281 the roto-translation from the ITS geometry, to convert
282 the local module coordinates into the global ones.
285 - kFALSE if the given cluster index points to a non-existing
286 cluster, or if the layer number makes no sense (< 0 or > 6).
287 - otherwise, kTRUE (meaning a successful operation).
289 **************************************************************************/
291 // checks if the layer is correct
292 Int_t ilayer = (refIndex & 0xf0000000) >> 28;
293 if (ilayer < 0 || ilayer >= fNLayers) {
294 Error("GetGlobalXYZ", "Wrong layer number: %d [range: %d - %d]", ilayer, 0, fNLayers);
297 // checks if the referenced cluster exists and corresponds to the passed reference
298 AliITSclusterV2 *refCluster = (AliITSclusterV2*) GetCluster(refIndex);
300 Error("GetGlobalXYZ", "Cluster not found for index %d", refIndex);
304 // determine the detector number
305 Int_t detID = refCluster->GetDetectorIndex() + fFirstModInLayer[ilayer];
307 // get rotation matrix
309 fGeom->GetRotMatrix(detID, rot);
311 // get translation vector
313 fGeom->GetTrans(detID, tx, ty, tz);
315 // determine r and phi for the reference conversion
316 Double_t r = -(Double_t)tx * rot[1] + (Double_t)ty * rot[0];
317 if (ilayer == 0) r = -r;
318 Double_t phi = TMath::ATan2(rot[1],rot[0]);
319 if (ilayer == 0) phi -= fgkPi;
321 // sets values for X, Y, Z in global coordinates and their errors
322 Double_t cosPhi = TMath::Cos(phi);
323 Double_t sinPhi = TMath::Sin(phi);
324 x = r*cosPhi + refCluster->GetY()*sinPhi;
325 y = -r*sinPhi + refCluster->GetY()*cosPhi;
326 z = refCluster->GetZ();
327 ex2 = refCluster->GetSigmaY2()*sinPhi*sinPhi;
328 ey2 = refCluster->GetSigmaY2()*cosPhi*cosPhi;
329 ez2 = refCluster->GetSigmaZ2();
334 //__________________________________________________________________________________
335 AliITStrackerANN::AliITSnode* AliITStrackerANN::AddNode(Int_t refIndex)
337 /**************************************************************************
339 GENERATOR OF NEURAL NODES
341 Fills the array of neural 'nodes', which are the ITS clusters
342 translated in the global reference frame.
343 Given that the global coordinates are used many times, they are
344 stored in a well-defined structure, in the form of an embedded class.
345 Moreover, this class allows a faster navigation among points
346 and neurons, by means of some object arrays, storing only the
347 neurons which start from, or end to, the given node.
348 Finally, each node contains all the other nodes which match it
349 with respect to the fixed walues, in order to perform a faster
350 neuron-creation phase.
353 1) reference index of the correlated AliITSclusterV2 object
356 - allocates the new AliITSnode objects
357 - initializes its object arrays
358 - from the global coordinates, calculates the
359 'phi' and 'theta' coordinates, in order to store it
360 into the correct theta-slot and azimutal sector.
363 - the pointer of the creater AliITSnode object
364 - in case of errors, a waring is given and a NULL is returned
366 **************************************************************************/
368 // create object and set the reference
369 AliITSnode *node = new AliITSnode;
371 Warning("AddNode", "Error occurred when allocating AliITSnode");
374 node->ClusterRef() = refIndex;
376 // calls the conversion function, which makes also some checks
377 // (layer number within range, existence of referenced cluster)
380 node->X(), node->Y(), node->Z(),
381 node->ErrX2(), node->ErrY2(), node->ErrZ2()
384 // initializes the object arrays
385 node->Matches() = new TObjArray;
386 node->InnerOf() = new TObjArray;
387 node->OuterOf() = new TObjArray;
389 // finds azimutal and polar sector (in degrees)
390 Double_t phi = node->GetPhi() * 180.0 / fgkPi;
391 Double_t theta = node->GetTheta() * 180.0 / fgkPi;
392 Int_t isector = (Int_t)(phi / fSectorWidth);
393 Int_t itheta = (Int_t)theta;
394 Int_t ilayer = (refIndex & 0xf0000000) >> 28;
396 // selects the right TObjArray to store object into
397 TObjArray *sector = (TObjArray*)fNodes[ilayer][itheta]->At(isector);
398 sector->AddLast(node);
403 //__________________________________________________________________________________
404 void AliITStrackerANN::CreateArrayStructure(Int_t nsectors)
406 /**************************************************************************
408 ARRAY STRUCTURE CREATOR
410 Creates a structure of nested TObjArray's where the AliITSnode's
412 - the first level is made by 6 arrays (one for each layer)
413 - the second level is made by 180 arrays (one for each int part of theta)
414 - the third level is made by a variable number of arrays
415 (one for each azimutal sector)
418 1) the number of azimutal sectors
421 - calculates the width of each sector, from the argument
422 - allocates and initializes all array levels
423 - sets a flag which tells the user if this NECESSARY operation
424 has been performed (it is needed BEFORE performing tracking)
426 **************************************************************************/
428 // Set the number of sectors and their width.
429 fSectorNum = nsectors;
430 fSectorWidth = 360.0 / (Double_t)fSectorNum;
431 if (fMsgLevel >= 2) {
432 cout << fSectorNum << " sectors --> sector width (degrees) = " << fSectorWidth << endl;
435 // Meaningful indexes
436 Int_t ilayer, isector, itheta;
438 // Mark for the created objects
439 TObjArray *sector = 0;
441 // First index: layer
442 fNodes = new TObjArray**[fNLayers];
443 for (ilayer = 0; ilayer < fNLayers; ilayer++) {
444 fNodes[ilayer] = new TObjArray*[180];
445 for (itheta = 0; itheta < 180; itheta++) fNodes[ilayer][itheta] = 0;
446 for (itheta = 0; itheta < 180; itheta++) {
447 fNodes[ilayer][itheta] = new TObjArray(nsectors);
448 for (isector = 0; isector < nsectors; isector++) {
449 sector = new TObjArray;
451 fNodes[ilayer][itheta]->AddAt(sector, isector);
456 // Sets a checking flag to TRUE.
457 // This flag is checked before filling up the arrays with the points.
458 fStructureOK = kTRUE;
461 //__________________________________________________________________________________
462 Int_t AliITStrackerANN::ArrangePoints(char *exportFile)
464 /**************************************************************************
468 This function assembles the operation from the other above methods,
469 and fills the arrays with the clusters already stored in the
470 layers of the tracker.
471 Then, in order to use this method, the user MUSTs call LoadClusters()
475 1) string for a file name where the global ccordinates
476 of all points can be exported (optional).
477 If this file must not be created, simply pass a NULL argument
480 - for each AliITSclusterV2 in each AliITSlayer, a ne AliITSnode
481 is created and stored in the correct location.
484 - the number of stored points
485 - when errors occur, or no points are found, 0 is returned
487 **************************************************************************/
489 // Check if the array structure has been created
491 Error("ArrangePoints", "Structure NOT defined. Call CreateArrayStructure() first");
495 // meaningful indexes
496 Int_t ientry, ilayer, nentries = 0, index;
499 // if the argument is not NULL, a file is opened
500 fstream file(exportFile, ios::out);
501 if (!exportFile || file.fail()) {
506 // scan all layers for node creation
507 for (ilayer = 0; ilayer < fNLayers; ilayer++) {
508 nPtsLayer = GetNumberOfClustersLayer(ilayer);
509 if (fMsgLevel >= 1) {
510 cout << "Layer " << ilayer << " --> " << nPtsLayer << " clusters" << endl;
512 for (ientry = 0; ientry < nPtsLayer; ientry++) {
513 // calculation of cluster index : (Bit mask LLLLIIIIIIIIIIII)
514 // (L = bits used for layer)
515 // (I = bits used for position in layer)
516 index = ilayer << 28;
518 // add new AliITSnode object
519 AliITSnode *n = AddNode(index);
520 if ( (n != NULL) && exportFile ) {
521 file << index << ' ' << n->X() << ' ' << n->Y() << ' ' << n->Z() << endl;
524 nentries += nPtsLayer;
527 // a conventional final message is put at the end of file
529 file << "-1 0.0 0.0 0.0" << endl;
533 // returns the number of points processed
537 //__________________________________________________________________________________
538 void AliITStrackerANN::StoreOverallMatches()
540 /**************************************************************************
544 Once the nodes have been created, a firs analysis is to check
545 what pairs will satisfy at least the 'fixed' cuts (theta, helix-match)
546 and the most permissive (= larger) curvature cut.
547 All these node pairs are suitable for neuron creation.
548 In fact, when performing a Neural Tracking step, the only further check
549 will be a check against the current curvature step, while the other
551 For thi purpose, each AliITSnode has a data member, named 'fMatches'
552 which contains references to all other AliITSnodes in the successive layer
553 that form, with it, a 'good' pair, with respect to the above cited cuts.
554 Then, in each step for neuron creation, the possible neurons starting from
555 each node will be searched ONLY within the nodes referenced in fMatches.
556 This, of course, speeds up a lot the neuron creation procedure, at the
557 cost of some memory occupation, which results not to be critical.
560 - for each AliITSnode, matches are found according to the criteria
561 expressed above, and stored in the node->fMatches array
563 **************************************************************************/
565 // meaningful counters
566 Int_t ilayer, isector, itheta1, itheta2, check;
567 TObjArray *list1 = 0, *list2 = 0;
568 AliITSnode *node1 = 0, *node2 = 0;
569 Double_t thetaMin, thetaMax;
572 // Scan for each sector
573 for (isector = 0; isector < fSectorNum; isector++) {
574 // sector is chosen once for both lists
575 for (ilayer = 0; ilayer < fNLayers - 1; ilayer++) {
576 for (itheta1 = 0; itheta1 < 180; itheta1++) {
577 list1 = (TObjArray*)fNodes[ilayer][itheta1]->At(isector);
578 TObjArrayIter iter1(list1);
579 while ( (node1 = (AliITSnode*)iter1.Next()) ) {
580 if (node1->IsUsed()) continue;
581 // clear an eventually already present array
582 // node1->Matches()->Clear();
583 // get the global coordinates and defines the theta interval from cut
584 thetaMin = (node1->GetTheta() * 180.0 / fgkPi) - fPolarInterval;
585 thetaMax = (node1->GetTheta() * 180.0 / fgkPi) + fPolarInterval;
586 imin = (Int_t)thetaMin;
587 imax = (Int_t)thetaMax;
588 if (imin < 0) imin = 0;
589 if (imax > 179) imax = 179;
590 // loop on the selected theta slots
591 for (itheta2 = imin; itheta2 <= imax; itheta2++) {
592 list2 = (TObjArray*)fNodes[ilayer + 1][itheta2]->At(isector);
593 TObjArrayIter iter2(list2);
594 while ( (node2 = (AliITSnode*)iter2.Next()) ) {
595 check = PassAllCuts(node1, node2, fCurvNum - 1, fVertexX, fVertexY, fVertexZ);
597 node1->Matches()->AddLast(node2);
599 } // while (node2...)
600 } // for (itheta2...)
601 } // while (node1...)
604 } // for (isector...)
607 //__________________________________________________________________________________
608 Int_t AliITStrackerANN::PassAllCuts
609 (AliITSnode *inner, AliITSnode *outer, Int_t curvStep,
610 Double_t vx, Double_t vy, Double_t vz)
612 // ***********************************************************************************
614 // This check is called in the above method for finding the matches of each node
615 // It check the passed point pair against all the fixed cuts and a specified
616 // curvature cut, among all the ones which have been defined.
617 // The cuts need a vertex-constraint, which is not absolute, but it is passed
621 // 1) the point in the inner layer
622 // 2) the point in the outer layer
623 // 3) curvature step for the curvature cut check (preferably the last)
624 // 4) X of the used vertex
625 // 5) Y of the used vertex
626 // 6) Z of the used vertex
629 // - if necessary, swaps the two points
630 // (the first must be in the innermost of the two layers)
631 // - checks for the theta cuts
632 // - calculates the circle passing through the vertex
633 // and the given points and checks for the curvature cut
634 // - using the radius calculated there, checks for the helix-math cut
637 // 0 - All cuts passed
638 // 1 - theta 2D cut not passed
639 // 2 - theta 3D cut not passed
640 // 3 - curvature calculated but cut not passed
641 // 4 - curvature not calculated (division by zero)
642 // 5 - helix cut not passed
643 // 6 - curvature inxed out of range
645 // ***********************************************************************************
647 // Check for curvature index
648 if (curvStep < 0 || curvStep >= fCurvNum) return 6;
650 // Swap points in order that r1 < r2
651 AliITSnode *temp = 0;
652 if (outer->GetLayer() < inner->GetLayer()) {
658 // All cuts are variable according to the layer of the
659 // innermost point (the other point will surely be
660 // in the further one, because we don't check poin pairs
661 // which are not in adjacent layers)
662 // The reference is given by the innermost point.
663 Int_t layRef = inner->GetLayer();
665 // The calculations in the transverse plane are made in
666 // a shifted reference frame, whose origin corresponds to
667 // the reference point passed in the argument.
668 Double_t xIn = inner->X() - vx;
669 Double_t xOut = outer->X() - vx;
670 Double_t yIn = inner->Y() - vy;
671 Double_t yOut = outer->Y() - vy;
672 Double_t zIn = inner->Z() - vz;
673 Double_t zOut = outer->Z() - vz;
674 Double_t rIn = TMath::Sqrt(xIn*xIn + yIn*yIn);
675 Double_t rOut = TMath::Sqrt(xOut*xOut + yOut*yOut);
677 // Check for theta cut.
678 // There are two different angular cuts:
679 // one w.r. to the angle in the 2-dimensional r-z plane...
681 TVector3 origin2innerRZ(zIn, rIn, 0.0);
682 TVector3 inner2outerRZ(zOut - zIn, rOut - rIn, 0.0);
683 dthetaRZ = origin2innerRZ.Angle(inner2outerRZ) * 180.0 / fgkPi;
684 if (dthetaRZ > fThetaCut2D[layRef]) {
686 // r-z theta cut not passed ---> 1
688 // ...and another w.r. to the angle in the 3-dimensional x-y-z space
690 TVector3 origin2innerXYZ(xIn, yIn, zIn);
691 TVector3 inner2outerXYZ(xOut - xIn, yOut - yIn, zOut - zIn);
692 dthetaXYZ = origin2innerXYZ.Angle(inner2outerXYZ) * 180.0 / fgkPi;
693 if (dthetaXYZ > fThetaCut3D[layRef]) {
695 // x-y-z theta cut not passed ---> 2
698 // Calculation & check of curvature
699 Double_t dx = xIn - xOut;
700 Double_t dy = yIn - yOut;
701 Double_t num = 2.0 * (xIn*yOut - xOut*yIn);
702 Double_t den = rIn*rOut*sqrt(dx*dx + dy*dy);
705 curv = TMath::Abs(num / den);
706 if (curv > fCurvCut[curvStep]) {
708 // curvature too large for cut ---> 3
712 Error("PassAllCuts", "Curvature calculations gives zero denominator");
714 // error: denominator = 0 ---> 4
717 // Calculation & check of helix matching
718 Double_t helMatch = 0.0;
719 Double_t arcIn = 2.0 * rIn * curv;
720 Double_t arcOut = 2.0 * rOut * curv;
721 if (arcIn > -1.0 && arcIn < 1.0)
722 arcIn = TMath::ASin(arcIn);
724 arcIn = ((arcIn > 0.0) ? 0.5 : 1.5) * TMath::Pi();
725 if (arcOut > -1.0 && arcOut < 1.0)
726 arcOut = TMath::ASin(arcOut);
728 arcOut = ((arcOut > 0.0) ? 0.5 : 1.5) * TMath::Pi();
730 arcOut /= 2.0 * curv;
731 if (arcIn == 0.0 || arcOut == 0.0) {
732 Error("PassAllCuts", "Calculation returns zero-length arcs: l1=%f, l2=%f", arcIn, arcOut);
734 // error: circumference arcs seem to equal zero ---> 4
736 helMatch = TMath::Abs(zIn / arcIn - zOut / arcOut);
737 if (helMatch > fHelixMatchCut[layRef]) {
739 // helix match cut not passed ---> 5
742 // ALL CUTS PASSED ---> 0
746 //__________________________________________________________________________________
747 Bool_t AliITStrackerANN::PassCurvCut
748 (AliITSnode *inner, AliITSnode *outer,
750 Double_t vx, Double_t vy, Double_t vz)
752 //***********************************************************************************
754 // This method operates essentially like the above one, but it is used
755 // during a single step of Neural Tracking, where the curvature cut
757 // Then, not necessaryly all the nodes stored in the fMatches array
758 // will be suitable for neuron creation in an intermediate step.
760 // It has the same arguments of the PassAllCuts() method, but
761 // the theta cut is not checked.
762 // Moreover, it has a boolean return value.
764 //***********************************************************************************
766 // Check for curvature index
767 if (curvStep < 0 || curvStep >= fCurvNum) return 6;
769 // Find the reference layer
770 Int_t layIn = inner->GetLayer();
771 Int_t layOut = outer->GetLayer();
772 Int_t layRef = (layIn < layOut) ? layIn : layOut;
774 // The calculations in the transverse plane are made in
775 // a shifted reference frame, whose origin corresponds to
776 // the reference point passed in the argument.
777 Double_t xIn = inner->X() - vx;
778 Double_t xOut = outer->X() - vx;
779 Double_t yIn = inner->Y() - vy;
780 Double_t yOut = outer->Y() - vy;
781 Double_t zIn = inner->Z() - vz;
782 Double_t zOut = outer->Z() - vz;
783 Double_t rIn = TMath::Sqrt(xIn*xIn + yIn*yIn);
784 Double_t rOut = TMath::Sqrt(xOut*xOut + yOut*yOut);
786 // Calculation & check of curvature
787 Double_t dx = xIn - xOut;
788 Double_t dy = yIn - yOut;
789 Double_t num = 2.0 * (xIn*yOut - xOut*yIn);
790 Double_t den = rIn*rOut*sqrt(dx*dx + dy*dy);
794 curv = TMath::Abs(num / den);
795 if (curv > fCurvCut[curvStep]) return kFALSE;
803 curv = TMath::Abs(num / den);
804 if (curv > fCurvCut[curvStep]) {
809 Error("PassAllCuts", "Curvature calculations gives zero denominator");
813 // Calculation & check of helix matching
814 Double_t helMatch = 0.0;
815 Double_t arcIn = 2.0 * rIn * curv;
816 Double_t arcOut = 2.0 * rOut * curv;
817 if (arcIn > -1.0 && arcIn < 1.0)
818 arcIn = TMath::ASin(arcIn);
820 arcIn = ((arcIn > 0.0) ? 0.5 : 1.5) * TMath::Pi();
821 if (arcOut > -1.0 && arcOut < 1.0)
822 arcOut = TMath::ASin(arcOut);
824 arcOut = ((arcOut > 0.0) ? 0.5 : 1.5) * TMath::Pi();
826 arcOut /= 2.0 * curv;
827 if (arcIn == 0.0 || arcOut == 0.0) {
828 Error("PassAllCuts", "Calculation returns zero-length arcs: l1=%f, l2=%f", arcIn, arcOut);
830 // error: circumference arcs seem to equal zero ---> 4
832 helMatch = TMath::Abs(zIn / arcIn - zOut / arcOut);
833 return (helMatch <= fHelixMatchCut[layRef]);
836 //__________________________________________________________________________________
837 void AliITStrackerANN::PrintMatches(Bool_t stop)
839 // Prints the list of points which appear to match
840 // each one of them, according to the preliminary
842 // The arguments states if a pause is required after printing
843 // the matches for each one. In this case, a keypress is needed.
845 TObjArray *sector = 0;
846 Int_t ilayer, isector, itheta, nF;
847 AliITSnode *node1 = 0, *node2 = 0;
848 //AliITSclusterV2 *cluster1 = 0, *cluster2 = 0;
850 for (ilayer = 0; ilayer < 6; ilayer++) {
851 for (isector = 0; isector < fSectorNum; isector++) {
852 for (itheta = 0; itheta < 180; itheta++) {
853 sector = (TObjArray*)fNodes[ilayer][itheta]->At(isector);
854 TObjArrayIter points(sector);
855 while ( (node1 = (AliITSnode*)points.Next()) ) {
856 nF = (Int_t)node1->Matches()->GetEntries();
857 cout << "Node layer: " << node1->GetLayer() << " --> ";
859 cout << "NO Matches!!!" << endl;
862 cout << nF << " Matches" << endl;
863 cout << "Reference cluster: " << hex << node1->ClusterRef() << endl;
864 TObjArrayIter matches(node1->Matches());
865 while ( (node2 = (AliITSnode*)matches.Next()) ) {
866 cout << "Match with " << hex << node2->ClusterRef() << endl;
869 cout << "Press a key" << endl;
878 //__________________________________________________________________________________
879 void AliITStrackerANN::ResetNodes(Int_t isector)
881 /***********************************************************************************
883 NODE NEURON ARRAY CLEANER
885 After a neural tracking step, this method
886 clears the arrays 'fInnerOf' and 'fOuterOf' of each AliITSnode
889 - the sector where the operation is being executed
891 ***********************************************************************************/
893 Int_t ilayer, itheta;
894 TObjArray *sector = 0;
895 AliITSnode *node = 0;
896 for (ilayer = 0; ilayer < fNLayers; ilayer++) {
897 for (itheta = 0; itheta < 180; itheta++) {
898 sector = (TObjArray*)fNodes[ilayer][itheta]->At(isector);
899 TObjArrayIter iter(sector);
901 node = (AliITSnode*)iter.Next();
903 node->InnerOf()->Clear();
904 node->OuterOf()->Clear();
906 delete node->InnerOf();
907 delete node->OuterOf();
908 node->InnerOf() = new TObjArray;
909 node->OuterOf() = new TObjArray;
916 //__________________________________________________________________________________
917 Int_t AliITStrackerANN::CreateNeurons
918 (AliITSnode *node, Int_t curvStep, Double_t vx, Double_t vy, Double_t vz)
920 // This method is used to create alle suitable neurons starting from
921 // a single AliITSnode. Each unit is also stored in the fInnerOf array
922 // of the passed node, and in the fOuterOf array of the other neuron edge.
923 // In the new implementation of the intermediate check steps, a further one
924 // is made, which chechs how well a helix is matched by three points
925 // in three consecutive layers.
926 // Then, a vertex-constrained check is made with vertex located
927 // in a layer L, for points in layers L+1 and L+2.
929 // In order to do this, the creator works recursively, in a tree-visit like operation.
930 // The neurons are effectively created only if the node argument passed is in
931 // the 5th layer (they are created between point of 5th and 6th layer).
932 // If the node is in an inner layer, its coordinates are passet as vertex for a nested
933 // call of the same function in the next two layers.
937 // 2) current curvature cut step
938 // 3) X of referenced temporary (not primary) vertex
939 // 4) Y of referenced temporary (not primary) vertex
940 // 5) Z of referenced temporary (not primary) vertex
943 // - if the layer is the 5th, neurons are created with nodes
944 // in the fMatches array of the passed node
945 // - otherwise, the X, Y, Z of the passed node are given as
946 // vertex and the same procedure is recursively called for all
947 // nodes in the fMatches array of the passed one.
950 // - the total number of neurons created from the passed one
951 // summed with all neurons created from all nodes well matched with it
952 // (assumes a meaning only for nodes in the first layer)
955 Int_t made = 0; // counter
956 Bool_t found = 0; // flag
957 AliITSnode *match = 0; // cursor for a AliITSnode array
958 AliITSneuron *unit = 0; // cursor for a AliITSneuron array
960 // --> Case 0: the passed node has already been used
961 // as member of a track found in a previous step.
962 // In this case, of course, the function exits.
963 if (node->IsUsed()) return 0;
965 // --> Case 1: there are NO well-matched points.
966 // This can happen in all ITS layers, but it happens **for sure**
967 // for a node in the 'last' layer.
968 // Even in this case, the function exits.
969 if (node->Matches()->IsEmpty()) return 0;
971 // --> Case 2: there are well-matched points.
972 // In this case, the function creates a neuron for each
973 // well-matched pair (according to the cuts for the current step)
974 // Moreover, before storing the neuron, a check is necessary
975 // to avoid the duplicate creation of the same neuron twice.
976 // (This could happen if the 3 last arguments of the function
977 // are close enough to cause a good match for the current step
978 // between two points, independently of their difference).
979 // Finally, a node is skipped if it has already been used.
980 // For each matched point for which a neuron is created, the procedure is
981 // recursively called.
982 TObjArrayIter matches(node->Matches());
983 while ( (match = (AliITSnode*)matches.Next()) ) {
984 if (match->IsUsed()) continue;
985 if (!PassCurvCut(node, match, curvStep, vx, vy, vz)) continue;
987 if (!node->InnerOf()->IsEmpty()) {
988 TObjArrayIter presentUnits(node->InnerOf());
989 while ( (unit = (AliITSneuron*)presentUnits.Next()) ) {
990 if (unit->Inner() == node && unit->Outer() == match) {
997 AliITSneuron *unit = new AliITSneuron(node, match, fEdge2, fEdge1);
998 fNeurons->AddLast(unit);
999 node->InnerOf()->AddLast(unit);
1000 match->OuterOf()->AddLast(unit);
1001 made += CreateNeurons(match, curvStep, node->X(), node->Y(), node->Z());
1005 // Of course, the return value contains the number of neurons
1006 // counting in also the oned created in all levels of recursive calls.
1010 //__________________________________________________________________________________
1011 Int_t AliITStrackerANN::CreateNetwork(Int_t sector, Int_t curvStep)
1013 // This function simply recalls the CreateNeurons() method for each node
1014 // in the first layer, for the current sector.
1015 // This generates the whole network, thanks to the recursive calls.
1018 // 1) current sector
1019 // 2) current curvature step
1022 // - scans the nodes array for all theta's in the current sector
1023 // and layer 0, and calls the CreateNeurons() function.
1025 // removes all eventually present neurons
1026 if (fNeurons) delete fNeurons;
1027 fNeurons = new TObjArray;
1028 fNeurons->SetOwner(kTRUE);
1030 // calls the ResetNodes() function to free the AliITSnode arrays
1031 if (fMsgLevel >= 2) {
1032 cout << "Sector " << sector << " PHI = ";
1033 cout << fSectorWidth * (Double_t)sector << " --> ";
1034 cout << fSectorWidth * (Double_t)(sector + 1) << endl;
1035 cout << "Curvature step " << curvStep << " [cut = " << fCurvCut[curvStep] << "]" << endl;
1039 // meaningful counters
1040 Int_t itheta, neurons = 0;
1041 TObjArray *lstSector = 0;
1044 Double_t vx[6], vy[6], vz[6];
1045 AliITSnode *p[6] = {0, 0, 0, 0, 0, 0};
1046 for (itheta = 0; itheta < 180; itheta++) {
1047 lstSector = (TObjArray*)fNodes[0][itheta]->At(sector);
1048 TObjArrayIter lay0(lstSector);
1049 while ( (p[0] = (AliITSnode*)lay0.Next()) ) {
1050 if (p[0]->IsUsed()) continue;
1054 neurons += CreateNeurons(p[0], curvStep, fVertexX, fVertexY, fVertexZ);
1056 TObjArrayIter lay1(p[0]->Matches());
1057 while ( (p[1] = (AliITSnode*)lay1.Next()) ) {
1058 if (p[1]->IsUsed()) continue;
1059 if (!PassCurvCut(p[0], p[1], curvStep, vx[0], vy[0], vz[0])) continue;
1060 unit = new AliITSneuron;
1061 unit->Inner() = p[0];
1062 unit->Outer() = p[1];
1063 unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1064 unit->fGain = new TObjArray;
1065 fNeurons->AddLast(unit);
1066 p[0]->InnerOf()->AddLast(unit);
1067 p[1]->OuterOf()->AddLast(unit);
1072 TObjArrayIter lay2(p[1]->Matches());
1073 while ( (p[2] = (AliITSnode*)lay2.Next()) ) {
1074 if (p[2]->IsUsed()) continue;
1075 if (!PassCurvCut(p[1], p[2], curvStep, vx[1], vy[1], vz[1])) continue;
1076 unit = new AliITSneuron;
1077 unit->Inner() = p[1];
1078 unit->Outer() = p[2];
1079 unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1080 unit->fGain = new TObjArray;
1081 fNeurons->AddLast(unit);
1082 p[1]->InnerOf()->AddLast(unit);
1083 p[2]->OuterOf()->AddLast(unit);
1088 TObjArrayIter lay3(p[2]->Matches());
1089 while ( (p[3] = (AliITSnode*)lay3.Next()) ) {
1090 if (p[3]->IsUsed()) continue;
1091 if (!PassCurvCut(p[2], p[3], curvStep, vx[2], vy[2], vz[2])) continue;
1092 unit = new AliITSneuron;
1093 unit->Inner() = p[2];
1094 unit->Outer() = p[3];
1095 unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1096 unit->fGain = new TObjArray;
1097 fNeurons->AddLast(unit);
1098 p[2]->InnerOf()->AddLast(unit);
1099 p[3]->OuterOf()->AddLast(unit);
1104 TObjArrayIter lay4(p[3]->Matches());
1105 while ( (p[4] = (AliITSnode*)lay4.Next()) ) {
1106 if (p[4]->IsUsed()) continue;
1107 if (!PassCurvCut(p[3], p[4], curvStep, vx[3], vy[3], vz[3])) continue;
1108 unit = new AliITSneuron;
1109 unit->Inner() = p[3];
1110 unit->Outer() = p[4];
1111 unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1112 unit->fGain = new TObjArray;
1113 fNeurons->AddLast(unit);
1114 p[3]->InnerOf()->AddLast(unit);
1115 p[4]->OuterOf()->AddLast(unit);
1120 TObjArrayIter lay5(p[4]->Matches());
1121 while ( (p[5] = (AliITSnode*)lay5.Next()) ) {
1122 if (p[5]->IsUsed()) continue;
1123 if (!PassCurvCut(p[4], p[5], curvStep, vx[4], vy[4], vz[4])) continue;
1124 unit = new AliITSneuron;
1125 unit->Inner() = p[4];
1126 unit->Outer() = p[5];
1127 unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1128 unit->fGain = new TObjArray;
1129 fNeurons->AddLast(unit);
1130 p[4]->InnerOf()->AddLast(unit);
1131 p[5]->OuterOf()->AddLast(unit);
1140 } // for (itheta...)
1141 // END OF NEW VERSION
1144 for (ilayer = 0; ilayer < 6; ilayer++) {
1145 for (itheta = 0; itheta < 180; itheta++) {
1146 lstSector = (TObjArray*)fNodes[ilayer][itheta]->At(sector_idx);
1147 TObjArrayIter inners(lstSector);
1148 while ( (inner = (AliITSnode*)inners.Next()) ) {
1149 if (inner->GetUser() >= 0) continue;
1150 TObjArrayIter outers(inner->Matches());
1151 while ( (outer = (AliITSnode*)outers.Next()) ) {
1152 if (outer->GetUser() >= 0) continue;
1153 if (!PassCurvCut(inner, outer, curvStep, fVX, fVY, fVZ)) continue;
1154 unit = new AliITSneuron;
1155 unit->Inner() = inner;
1156 unit->Outer() = outer;
1157 unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1158 unit->fGain = new TObjArray;
1159 fNeurons->AddLast(unit);
1160 inner->InnerOf()->AddLast(unit);
1161 outer->OuterOf()->AddLast(unit);
1165 } // for (itheta...)
1166 } // for (ilayer...)
1169 fNeurons->SetOwner();
1173 //__________________________________________________________________________________
1174 Int_t AliITStrackerANN::LinkNeurons()
1176 /***********************************************************************************
1180 Scans the whole neuron array, in order to find all neuron pairs
1181 which are connected in sequence and share a positive weight.
1182 For each of them, an AliITSlink is created, which stores
1183 the weight value, and will allow for a faster calculation
1184 of the total neural input for each updating cycle.
1186 Every neuron contains an object array which stores all AliITSlink
1187 objects which point to sequenced units, with the respective weights.
1190 - the number of link created for the neural network.
1191 If they are 0, no updating can be done and the step is skipped.
1193 ***********************************************************************************/
1195 // meaningful indexes
1197 Double_t weight = 0.0;
1198 TObjArrayIter neurons(fNeurons), *iter;
1199 AliITSneuron *neuron = 0, *test = 0;
1201 // scan in the neuron array
1203 neuron = (AliITSneuron*)neurons.Next();
1205 // checks for neurons startin from its head ( -> )
1206 iter = (TObjArrayIter*)neuron->Inner()->OuterOf()->MakeIterator();
1208 test = (AliITSneuron*)iter->Next();
1210 weight = Weight(test, neuron);
1211 if (weight > 0.0) neuron->Gain()->AddLast(new AliITSlink(weight, test));
1215 // checks for neurons ending in its tail ( >- )
1216 iter = (TObjArrayIter*)neuron->Outer()->InnerOf()->MakeIterator();
1218 test = (AliITSneuron*)iter->Next();
1220 weight = Weight(neuron, test);
1221 if (weight > 0.0) neuron->Gain()->AddLast(new AliITSlink(weight, test));
1229 //__________________________________________________________________________________
1230 Bool_t AliITStrackerANN::Update()
1232 /***********************************************************************************
1234 Performs a single updating cycle.
1237 - for each neuron, gets the activation with the neuron Activate() method
1238 - checks if stability has been reached (compare mean activation variation
1239 with the stability threshold data member)
1242 - kTRUE means that the neural network has stabilized
1243 - kFALSE means that another updating cycle is needed
1245 ***********************************************************************************/
1247 Double_t actVar = 0.0, totDiff = 0.0;
1248 TObjArrayIter iter(fNeurons);
1251 unit = (AliITSneuron*)iter.Next();
1253 actVar = unit->Activate(fTemperature);
1254 // calculation the relative activation variation
1257 totDiff /= fNeurons->GetSize();
1258 return (totDiff < fStabThreshold);
1261 //__________________________________________________________________________________
1262 void AliITStrackerANN::FollowChains(Int_t sector)
1264 /***********************************************************************************
1268 After that the neural network has stabilized,
1269 the final step is to create polygonal chains
1270 of clusters, one in each layer, which represent
1271 the tracks recognized by the neural algorithm.
1272 This is made by means of a choice of the best
1273 neuron among the ones starting from each point.
1275 Once that such neuron is selected, its inner point
1276 will set the 'fNext' field to its outer point, and
1277 similarly, its outer point will set the 'fPrev' field
1279 This defines a bi-directional sequence.
1281 In this procedure, it can happen that many neurons
1282 which have the head of the arrow in a given node, will
1283 all select as best following the neuron with the largest
1284 activation starting in that point.
1285 This results in MANY nodes which have the same 'fNext'.
1286 But, this field will be set to NULL for all these points,
1287 but the only one which is pointed by the 'fPrev' field
1288 of this shared node.
1290 ***********************************************************************************/
1292 // meaningful counters
1293 Int_t itheta, ilayer;
1294 TObjArray *lstSector = 0;
1295 Double_t test = fActMinimum;
1297 AliITSneuron *n = 0;
1299 // scan the whole collection of nodes
1300 for (ilayer = 0; ilayer < fNLayers; ilayer++) {
1301 for (itheta = 0; itheta < 180; itheta++) {
1302 // get the array of point in a given layer/theta-slot/sector
1303 lstSector = (TObjArray*)fNodes[ilayer][itheta]->At(sector);
1304 TObjArrayIter nodes(lstSector);
1305 while ( (p = (AliITSnode*)nodes.Next()) ) {
1306 // if the point is used, it is skipped
1307 if (p->IsUsed()) continue;
1308 // initially, fNext points to nothing, and
1309 // the comparison value is set to the minimum activation
1310 // which allows to say that a neuron is turned 'on'
1311 // a node from which only 'off' neurons start is probably
1312 // a noise point, which will be excluded from the reconstruction.
1315 TObjArrayIter innerof(p->InnerOf());
1316 while ( (n = (AliITSneuron*)innerof.Next()) ) {
1317 // if the examined neuron has not the largest activation
1318 // it is skipped and removed from array of all neurons
1319 // and of its outer point (its inner is the cursor p)
1320 if (n->Activation() < test) {
1321 p->InnerOf()->Remove(n);
1322 n->Outer()->OuterOf()->Remove(n);
1323 delete fNeurons->Remove(n);
1326 // otherwise, its activation becomes the maximum reference
1327 p->Next() = n->Outer();
1328 // at the exit of the while(), the fNext will point
1329 // to the outer node of the neuron starting in p, whose
1330 // activation is the largest.
1332 // the same procedure is made now for all neurons
1333 // for which p is the outer point
1336 TObjArrayIter outerof(p->OuterOf());
1337 while ( (n = (AliITSneuron*)outerof.Next()) ) {
1338 // if the examined neuron has not the largest activation
1339 // it is skipped and removed from array of all neurons
1340 // and of its inner point (its outer is the cursor p)
1341 if (n->Activation() < test) {
1342 p->OuterOf()->Remove(n);
1343 n->Inner()->InnerOf()->Remove(n);
1344 delete fNeurons->Remove(n);
1347 // otherwise, its activation becomes the maximum reference
1348 p->Prev() = n->Inner();
1349 // at the exit of the while(), the fPrev will point
1350 // to the inner node of the neuron ending in p, whose
1351 // activation is the largest.
1353 } // end while (p ...)
1354 } // end for (itheta ...)
1355 } // end for (ilayer ...)
1357 // now the mismatches are solved
1358 Bool_t matchPrev, matchNext;
1359 for (ilayer = 0; ilayer < fNLayers; ilayer++) {
1360 for (itheta = 0; itheta < 180; itheta++) {
1361 // get the array of point in a given layer/theta-slot/sector
1362 lstSector = (TObjArray*)fNodes[ilayer][itheta]->At(sector);
1363 TObjArrayIter nodes(lstSector);
1364 while ( (p = (AliITSnode*)nodes.Next()) ) {
1365 // now p will point to a fPrev and a fNext node.
1366 // Ideally they are placed this way: fPrev --> P --> fNext
1367 // A mismatch happens if the point addressed as fPrev does NOT
1368 // point to p as its fNext. And the same for the point addressed
1370 // In this case, the fNext and fPrev pointers are set to NULL
1371 // and p is excluded from the reconstruction
1372 matchPrev = matchNext= kFALSE;
1373 if (ilayer > 0 && p->Prev() != NULL)
1374 if (p->Prev()->Next() == p) matchPrev = kTRUE;
1375 if (ilayer < 5 && p->Next() != NULL)
1376 if (p->Next()->Prev() == p) matchNext = kTRUE;
1379 else if (ilayer == 5)
1381 if (!matchNext || !matchPrev) {
1382 p->Prev() = p->Next() = 0;
1384 } // end while (p ...)
1385 } // end for (itheta ...)
1386 } // end for (ilayer ...)
1389 //__________________________________________________________________________________
1390 Int_t AliITStrackerANN::SaveTracks(Int_t sector)
1392 /********************************************************************************
1396 Using the fNext and fPrev pointers, the chain is followed
1397 and the track is fitted and saved.
1398 Of course, the track is followed as a chain with a point
1399 for each layer, then the track following starts always
1400 from the clusters in layer 0.
1402 ***********************************************************************************/
1404 // if not initialized, the tracks TobjArray is initialized
1405 if (!fFoundTracks) fFoundTracks = new TObjArray;
1407 // meaningful counters
1408 Int_t itheta, ilayer, l;
1409 TObjArray *lstSector = 0;
1410 AliITSnode *p = 0, *q = 0, **node = new AliITSnode*[fNLayers];
1411 for (l = 0; l < fNLayers; l++) node[l] = 0;
1414 array = new AliITSnode*[fNLayers + 1];
1415 for (l = 0; l <= fNLayers; l++) array[l] = 0;
1416 array[0] = new AliITSnode();
1417 array[0]->X() = fVertexX;
1418 array[0]->Y() = fVertexY;
1419 array[0]->Z() = fVertexZ;
1420 array[0]->ErrX2() = fVertexErrorX;
1421 array[0]->ErrY2() = fVertexErrorY;
1422 array[0]->ErrZ2() = fVertexErrorZ;
1424 Double_t *param = new Double_t[8];
1426 // scan the whole collection of nodes
1427 for (ilayer = 0; ilayer < 1; ilayer++) {
1428 for (itheta = 0; itheta < 180; itheta++) {
1429 // get the array of point in a given layer/theta-slot/sector
1430 lstSector = (TObjArray*)fNodes[ilayer][itheta]->At(sector);
1431 TObjArrayIter nodes(lstSector);
1432 while ( (p = (AliITSnode*)nodes.Next()) ) {
1433 for (q = p; q; q = q->Next()) {
1437 //if (!RiemannFit(fNLayers, node, param)) continue;
1438 // initialization of Kalman Filter Tracking
1439 AliITSclusterV2 *cluster = (AliITSclusterV2*)GetCluster(node[0]->ClusterRef());
1440 Int_t mod = cluster->GetDetectorIndex();
1441 Int_t lay, lad, det;
1442 fGeom->GetModuleId(mod, lay, lad, det);
1443 Float_t y0 = cluster->GetY();
1444 Float_t z0 = cluster->GetZ();
1445 AliITStrackSA* trac = new AliITStrackSA(lay, lad, det,
1447 param[4], param[7], param[3], 1);
1448 for (l = 0; l < fNLayers; l++) {
1449 cluster = (AliITSclusterV2*)GetCluster(node[l]->ClusterRef());
1450 if (cluster) trac->AddClusterV2(l, (node[l]->ClusterRef() & 0x0fffffff)>>0);
1452 AliITStrackV2* ot = new AliITStrackV2(*trac);
1453 ot->ResetCovariance();
1454 ot->ResetClusters();
1455 if (RefitAt(49.,ot,trac)) { //fit from layer 1 to layer 6
1456 AliITStrackV2 *otrack2 = new AliITStrackV2(*ot);
1457 otrack2->ResetCovariance();
1458 otrack2->ResetClusters();
1459 //fit from layer 6 to layer 1
1460 if (RefitAt(3.7,otrack2,ot)) fFoundTracks->AddLast(otrack2);
1462 // end of Kalman Filter fit
1470 //__________________________________________________________________________________
1471 void AliITStrackerANN::ExportTracks(const char *filename) const
1473 // Exports found tracks into a TTree of AliITStrackV2 objects
1474 TFile *file = new TFile(filename, "RECREATE");
1475 TTree *tree = new TTree("TreeT-ANN", "Tracks found in ITS stand-alone with Neural Tracking");
1476 AliITStrackV2 *track = 0;
1477 tree->Branch("Tracks", &track, "AliITStrackV2");
1478 TObjArrayIter tracks(fFoundTracks);
1479 while ( (track = (AliITStrackV2*)tracks.Next()) ) {
1488 //__________________________________________________________________________________
1489 void AliITStrackerANN::CleanNetwork()
1491 // Removes deactivated units from the network
1493 AliITSneuron *unit = 0;
1494 TObjArrayIter neurons(fNeurons);
1495 while ( (unit = (AliITSneuron*)neurons.Next()) ) {
1496 if (unit->Activation() < fActMinimum) {
1497 unit->Inner()->InnerOf()->Remove(unit);
1498 unit->Outer()->OuterOf()->Remove(unit);
1499 delete fNeurons->Remove(unit);
1505 AliITSneuron *enemy = 0;
1507 while ( (unit = (AliITSneuron*)neurons.Next()) ) {
1508 nIn = (Int_t)unit->Inner()->InnerOf()->GetSize();
1509 nOut = (Int_t)unit->Outer()->OuterOf()->GetSize();
1510 if (nIn < 2 && nOut < 2) continue;
1513 TObjArrayIter competing(unit->Inner()->InnerOf());
1514 while ( (enemy = (AliITSneuron*)competing.Next()) ) {
1515 if (unit->Activation() > enemy->Activation()) {
1516 enemy->Inner()->InnerOf()->Remove(enemy);
1517 enemy->Outer()->OuterOf()->Remove(enemy);
1518 delete fNeurons->Remove(enemy);
1521 unit->Inner()->InnerOf()->Remove(unit);
1522 unit->Outer()->OuterOf()->Remove(unit);
1523 delete fNeurons->Remove(unit);
1528 if (removed) continue;
1531 TObjArrayIter competing(unit->Outer()->OuterOf());
1532 while ( (enemy = (AliITSneuron*)competing.Next()) ) {
1533 if (unit->Activation() > enemy->Activation()) {
1534 enemy->Inner()->InnerOf()->Remove(enemy);
1535 enemy->Outer()->OuterOf()->Remove(enemy);
1536 delete fNeurons->Remove(enemy);
1539 unit->Inner()->InnerOf()->Remove(unit);
1540 unit->Outer()->OuterOf()->Remove(unit);
1541 delete fNeurons->Remove(unit);
1550 //__________________________________________________________________________________
1551 Int_t AliITStrackerANN::StoreTracks()
1553 // Stores the tracks found in a single neural tracking step.
1554 // In order to do this, it sects each neuron which has a point
1555 // in the first layer.
1558 // if not initialized, the tracks TobjArray is initialized
1559 if (!fFoundTracks) fFoundTracks = new TObjArray;
1561 Int_t i, check, stored = 0;
1562 Double_t testAct = 0;
1563 AliITSneuron *unit = 0, *cursor = 0, *fwd = 0;
1564 AliITSnode *node = 0;
1565 TObjArrayIter iter(fNeurons), *fwdIter;
1566 TObjArray *removedUnits = new TObjArray(0);
1567 removedUnits->SetOwner(kFALSE);
1568 AliITStrackANN annTrack(fNLayers);
1571 unit = (AliITSneuron*)iter.Next();
1573 if (unit->Inner()->GetLayer() > 0) continue;
1574 annTrack.SetNode(unit->Inner()->GetLayer(), unit->Inner());
1575 annTrack.SetNode(unit->Outer()->GetLayer(), unit->Outer());
1576 node = unit->Outer();
1577 removedUnits->AddLast(unit);
1579 testAct = fActMinimum;
1580 fwdIter = (TObjArrayIter*)node->InnerOf()->MakeIterator();
1583 cursor = (AliITSneuron*)fwdIter->Next();
1585 if (cursor->Used()) continue;
1586 if (cursor->Activation() >= testAct) {
1587 testAct = cursor->Activation();
1592 removedUnits->AddLast(fwd);
1593 node = fwd->Outer();
1594 annTrack.SetNode(node->GetLayer(), node);
1596 check = annTrack.CheckOccupation();
1600 //if (!RiemannFit(fNLayers, trackitem, param)) continue;
1601 if (!annTrack.RiemannFit()) continue;
1602 // initialization of Kalman Filter Tracking
1603 AliITSclusterV2 *cluster = (AliITSclusterV2*)GetCluster(annTrack[0]->ClusterRef());
1604 Int_t mod = cluster->GetDetectorIndex();
1605 Int_t lay, lad, det;
1606 fGeom->GetModuleId(mod, lay, lad, det);
1607 Float_t y0 = cluster->GetY();
1608 Float_t z0 = cluster->GetZ();
1609 AliITStrackSA* trac = new AliITStrackSA(lay, lad, det, y0, z0,
1610 annTrack.Phi(), annTrack.TanLambda(),
1611 annTrack.Curv(), 1);
1612 for (Int_t l = 0; l < fNLayers; l++) {
1613 if (!annTrack[l]) continue;
1614 cluster = (AliITSclusterV2*)GetCluster(annTrack[l]->ClusterRef());
1615 if (cluster) trac->AddClusterV2(l, (annTrack[l]->ClusterRef() & 0x0fffffff)>>0);
1617 AliITStrackV2* ot = new AliITStrackV2(*trac);
1618 ot->ResetCovariance();
1619 ot->ResetClusters();
1620 if (RefitAt(49.,ot,trac)) { //fit from layer 1 to layer 6
1621 AliITStrackV2 *otrack2 = new AliITStrackV2(*ot);
1622 otrack2->ResetCovariance();
1623 otrack2->ResetClusters();
1624 //fit from layer 6 to layer 1
1625 if (RefitAt(3.7,otrack2,ot)) fFoundTracks->AddLast(otrack2);
1627 // end of Kalman Filter fit
1629 for (i = 0; i < fNLayers; i++) {
1630 //node = (AliITSnode*)removedPoints->At(i);
1634 fwdIter = (TObjArrayIter*)removedUnits->MakeIterator();
1636 cursor = (AliITSneuron*)fwdIter->Next();
1646 Double_t AliITStrackerANN::Weight(AliITSneuron *nAB, AliITSneuron *nBC)
1648 /***********************************************************************************
1649 * Calculation of neural weight.
1650 * The implementation of positive neural weight is set only in the case
1651 * of connected units (e.g.: A->B with B->C).
1652 * Given that B is the **common** point. care should be taken to pass
1653 * as FIRST argument the neuron going "to" B, and
1654 * as SECOND argument the neuron starting "from" B
1655 * anyway, a check is put in order to return 0.0 when arguments are not well given.
1656 ***********************************************************************************/
1658 if (nAB->Outer() != nBC->Inner()) {
1659 if (nBC->Outer() == nAB->Inner()) {
1660 AliITSneuron *temp = nAB;
1664 if (fMsgLevel >= 3) {
1665 Info("Weight", "Switching wrongly ordered arguments.");
1668 Warning("Weight", "Not connected segments. Returning 0.0");
1672 AliITSnode *pA = nAB->Inner();
1673 AliITSnode *pB = nAB->Outer();
1674 AliITSnode *pC = nBC->Outer();
1676 TVector3 vAB(pB->X() - pA->X(), pB->Y() - pA->Y(), pB->Z() - pA->Z());
1677 TVector3 vBC(pC->X() - pB->X(), pC->Y() - pB->Y(), pC->Z() - pB->Z());
1679 Double_t weight = 1.0 - sin(vAB.Angle(vBC));
1680 return fGain2CostRatio * TMath::Power(weight, fExponent);
1685 /******************************************
1686 ******************************************
1687 *** AliITStrackerANN::AliITSnode class ***
1688 ******************************************
1689 ******************************************/
1691 //__________________________________________________________________________________
1692 AliITStrackerANN::AliITSnode::AliITSnode()
1693 : fUsed(kFALSE), fClusterRef(-1),
1694 fMatches(NULL), fInnerOf(NULL), fOuterOf(NULL),
1695 fNext(NULL), fPrev(NULL)
1697 // Constructor for the embedded 'AliITSnode' class.
1698 // It initializes all pointer-like objects.
1701 fEX2 = fEY2 = fEZ2 = 0.0;
1704 //__________________________________________________________________________________
1705 AliITStrackerANN::AliITSnode::~AliITSnode()
1707 // Destructor for the embedded 'AliITSnode' class.
1708 // It should clear the object arrays, but it is possible
1709 // that some objects still are useful after the point deletion
1710 // then the arrays are cleared but their objects are owed by
1711 // another TCollection object, and not deleted.
1712 // For safety reasons, all the pointers are set to zero.
1714 fInnerOf->SetOwner(kFALSE);
1718 fOuterOf->SetOwner(kFALSE);
1722 fMatches->SetOwner(kFALSE);
1730 //__________________________________________________________________________________
1731 Double_t AliITStrackerANN::AliITSnode::GetPhi() const
1733 // Calculates the 'phi' (azimutal) angle, and returns it
1734 // in the range between 0 and 2Pi radians.
1737 q = TMath::ATan2(fY,fX);
1741 return q + TMath::TwoPi();
1744 //__________________________________________________________________________________
1745 Double_t AliITStrackerANN::AliITSnode::GetError(Option_t *option)
1747 // Returns the error or the square error of
1748 // values related to the coordinates in different systems.
1749 // The option argument specifies the coordinate error desired:
1751 // "R2" --> error in transverse radius
1752 // "R3" --> error in spherical radius
1753 // "PHI" --> error in azimuthal angle
1754 // "THETA" --> error in polar angle
1755 // "SQ" --> get the square of error
1757 // In order to get the error on the cartesian coordinates
1758 // reference to the inline ErrX2(), ErrY2() adn ErrZ2() methods.
1760 TString opt(option);
1761 Double_t errorSq = 0.0;
1764 if (opt.Contains("R2")) {
1765 errorSq = fX*fX*fEX2 + fY*fY*fEY2;
1766 errorSq /= GetR2sq();
1768 else if (opt.Contains("R3")) {
1769 errorSq = fX*fX*fEX2 + fY*fY*fEY2 + fZ*fZ*fEZ2;
1770 errorSq /= GetR3sq();
1772 else if (opt.Contains("PHI")) {
1773 errorSq = fY*fY*fEX2;
1774 errorSq += fX*fX*fEY2;
1775 errorSq /= GetR2sq() * GetR2sq();
1777 else if (opt.Contains("THETA")) {
1778 errorSq = fZ*fZ * (fX*fX*fEX2 + fY*fY*fEY2);
1779 errorSq += GetR2sq() * GetR2sq() * fEZ2;
1780 errorSq /= GetR3sq() * GetR3sq() * GetR2() * GetR2();
1783 if (!opt.Contains("SQ"))
1784 return TMath::Sqrt(errorSq);
1791 /********************************************
1792 ********************************************
1793 *** AliITStrackerANN::AliITSneuron class ***
1794 ********************************************
1795 ********************************************/
1797 //__________________________________________________________________________________
1798 AliITStrackerANN::AliITSneuron::AliITSneuron
1799 (AliITSnode *inner, AliITSnode *outer, Double_t minAct, Double_t maxAct)
1800 : fUsed(0), fInner(inner), fOuter(outer)
1802 // Default neuron constructor
1803 fActivation = gRandom->Rndm() * (maxAct-minAct) + minAct;
1804 fGain = new TObjArray;
1807 //__________________________________________________________________________________
1808 Double_t AliITStrackerANN::AliITSneuron::Activate(Double_t temperature)
1810 // This computes the new activation of a neuron, and returns
1811 // its activation variation as a consequence of the updating.
1814 // - the 'temperature' parameter for the neural activation logistic function
1817 // - collects the total gain, by summing the products
1818 // of the activation of each sequenced unit by the relative weight.
1819 // - collects the total cost, by summing the activations of
1820 // all competing units
1821 // - passes the sum of gain - cost to the activation function and
1822 // calculates the new activation
1825 // - the difference between the old activation and the new one
1829 Double_t sumGain = 0.0; // total contribution from chained neurons
1830 Double_t sumCost = 0.0; // total contribution from crossing neurons
1831 Double_t input; // total input
1832 Double_t actOld, actNew; // old and new values for the activation
1833 AliITSneuron *linked = 0; // cursor for scanning the neuron arrays (for link check)
1834 AliITSlink *link; // cursor for scanning the synapses arrays (for link check)
1835 TObjArrayIter *iterator = 0; // pointer to the iterator along the neuron arrays
1837 // sum contributions from the correlated units
1838 iterator = (TObjArrayIter*)fGain->MakeIterator();
1840 link = (AliITSlink*)iterator->Next();
1842 sumGain += link->Contribution();
1846 // sum contributions from the competing units:
1847 // the ones which have the same starting point...
1848 iterator = (TObjArrayIter*)fInner->InnerOf()->MakeIterator();
1850 linked = (AliITSneuron*)iterator->Next();
1852 if (linked == this) continue;
1853 sumCost += linked->fActivation;
1856 // ...and the ones which have the same ending point
1857 iterator = (TObjArrayIter*)fOuter->OuterOf()->MakeIterator();
1859 linked = (AliITSneuron*)iterator->Next();
1861 if (linked == this) continue;
1862 sumCost += linked->fActivation;
1865 // calculate the total input as the difference between gain and cost
1866 input = (sumGain - sumCost) / temperature;
1867 actOld = fActivation;
1868 // calculate the final output
1869 #ifdef NEURAL_LINEAR
1870 if (input <= -2.0 * temperature)
1872 else if (input >= 2.0 * temperature)
1875 actNew = input / (4.0 * temperature) + 0.5;
1877 actNew = 1.0 / (1.0 + TMath::Exp(-input));
1879 fActivation = actNew;
1881 // return the activation variation
1882 return TMath::Abs(actNew - actOld);
1887 /******************************************
1888 ******************************************
1889 *** AliITStrackerANN::AliITSlink class ***
1890 ******************************************
1891 ******************************************/
1893 // No methods defined non-inline
1897 /**********************************************
1898 **********************************************
1899 *** AliITStrackerANN::AliITStrackANN class ***
1900 **********************************************
1901 **********************************************/
1903 //__________________________________________________________________________________
1904 AliITStrackerANN::AliITStrackANN::AliITStrackANN(Int_t dim) : fNPoints(dim)
1906 // Default constructor for the AliITStrackANN class
1921 fNode = new AliITSnode*[dim];
1922 for (i = 0; i < dim; i++) fNode[i] = 0;
1926 //__________________________________________________________________________________
1927 Int_t AliITStrackerANN::AliITStrackANN::CheckOccupation() const
1929 // Returns the number of pointers fNode which are not NULL
1932 Int_t count = 0; // counter for how many points are stored in the track
1934 for (i = 0; i < fNPoints; i++) {
1935 if (fNode[i] != NULL) count++;
1941 //__________________________________________________________________________________
1942 Bool_t AliITStrackerANN::AliITStrackANN::RiemannFit()
1944 // Performs the Riemann Sphere fit for the given points to a circle
1945 // and then uses the fit parameters to fit a helix in space.
1948 // - kTRUE if all operations have been performed
1949 // - kFALSE if the numbers risk to lead to an arithmetic violation
1951 Int_t i, j, count, dim = fNPoints;
1953 // First check for all points
1954 count = CheckOccupation();
1955 if (count != fNPoints) {
1956 Error ("AliITStrackANN::RiemannFit", "CheckOccupations returns %d, fNPoints = %d ==> MISMATCH", count, fNPoints);
1962 for (i = 0; i < dim; i++) m1(i,0) = 1.0;
1964 // matrix of Rieman projection coordinates
1965 TMatrixD coords(dim,3);
1966 for (i = 0; i < dim; i++) {
1967 coords(i,0) = fNode[i]->X();
1968 coords(i,1) = fNode[i]->Y();
1969 coords(i,2) = fNode[i]->GetR2sq();
1972 // matrix of weights
1973 Double_t xterm, yterm, ex, ey;
1974 TMatrixD weights(dim,dim);
1975 for (i = 0; i < dim; i++) {
1976 xterm = fNode[i]->X() * fNode[i]->GetPhi() - fNode[i]->Y() / fNode[i]->GetR2();
1977 ex = fNode[i]->ErrX2();
1978 yterm = fNode[i]->Y() * fNode[i]->GetPhi() + fNode[i]->X() / fNode[i]->GetR2();
1979 ey = fNode[i]->ErrY2();
1980 weights(i,i) = fNode[i]->GetR2sq() / (xterm * xterm * ex + yterm * yterm * ey );
1983 // weighted sample mean
1984 Double_t meanX = 0.0, meanY = 0.0, meanW = 0.0, sw = 0.0;
1985 for (i = 0; i < dim; i++) {
1986 meanX += weights(i,i) * coords(i,0);
1987 meanY += weights(i,i) * coords(i,1);
1988 meanW += weights(i,i) * coords(i,2);
1995 // sample covariance matrix
1996 for (i = 0; i < dim; i++) {
1997 coords(i,0) -= meanX;
1998 coords(i,1) -= meanY;
1999 coords(i,2) -= meanW;
2001 TMatrixD coordsT(TMatrixD::kTransposed, coords);
2002 TMatrixD weights4coords(weights, TMatrixD::kMult, coords);
2003 TMatrixD sampleCov(coordsT, TMatrixD::kMult, weights4coords);
2004 for (i = 0; i < 3; i++) {
2005 for (j = i + 1; j < 3; j++) {
2006 sampleCov(i,j) = sampleCov(j,i) = (sampleCov(i,j) + sampleCov(j,i)) * 0.5;
2010 // Eigenvalue problem solving for V matrix
2012 TVectorD eval(3), n(3);
2013 #if ROOT_VERSION_CODE < ROOT_VERSION(4,0,2)
2014 TMatrixD evec = sampleCov.EigenVectors(eval);
2016 TMatrixDEigen ei(sampleCov);
2017 TMatrixD evec = ei.GetEigenVectors();
2018 eval = ei.GetEigenValuesRe();
2020 if (eval(1) < eval(ileast)) ileast = 1;
2021 if (eval(2) < eval(ileast)) ileast = 2;
2022 n(0) = evec(0, ileast);
2023 n(1) = evec(1, ileast);
2024 n(2) = evec(2, ileast);
2026 // c - known term in the plane intersection with Riemann axes
2027 Double_t c = -(meanX * n(0) + meanY * n(1) + meanW * n(2));
2029 // center and radius of fitted circle
2030 Double_t xc, yc, radius, curv;
2031 xc = -n(0) / (2. * n(2));
2032 yc = -n(1) / (2. * n(2));
2033 radius = (1. - n(2)*n(2) - 4.*c*n(2)) / (4. * n(2) * n(2));
2035 if (radius <= 0.E0) {
2036 Error("RiemannFit", "Radius = %f less than zero!!!", radius);
2039 radius = TMath::Sqrt(radius);
2040 curv = 1.0 / radius;
2042 // evaluating signs for curvature and others
2043 Double_t phi1 = 0.0, phi2, temp1, temp2, phi0, sumdphi = 0.0;
2044 AliITSnode *p = fNode[0];
2046 for (i = 1; i < dim; i++) {
2047 p = (AliITSnode*)fNode[i];
2052 if (temp1 > fgkPi && temp2 < fgkPi)
2054 else if (temp1 < fgkPi && temp2 > fgkPi)
2056 sumdphi += temp2 - temp1;
2059 if (sumdphi < 0.E0) curv = -curv;
2060 Double_t diff, angle = TMath::ATan2(yc, xc);
2062 phi0 = angle + 0.5 * TMath::Pi();
2064 phi0 = angle - 0.5 * TMath::Pi();
2065 diff = angle - phi0;
2067 Double_t dt, temp = TMath::Sqrt(xc*xc + yc*yc) - radius;
2072 //cout << "Dt = " << dt << endl;
2074 Double_t halfC = 0.5 * curv, test;
2075 Double_t *s = new Double_t[dim], *zz = new Double_t[dim], *ws = new Double_t[dim];
2076 for (j = 0; j < 6; j++) {
2080 s[j] = (p->GetR2sq() - dt * dt) / (1. + curv * dt);
2082 if (TMath::Abs(s[j]) < 1.E-6) s[j] = 0.;
2084 Error("RiemannFit", "Square root argument error: %17.15g < 0", s[j]);
2088 s[j] = TMath::Sqrt(s[j]);
2089 //cout << "Curv = " << halfC << " --- s[" << j << "] = " << s[j] << endl;
2091 test = TMath::Abs(s[j]);
2094 s[j] = ((s[j] > 0.) ? 0.99999999999 : -0.9999999999);
2096 Error("RiemannFit", "Value too large: %17.15g", s[j]);
2102 s[j] = TMath::ASin(s[j]) / halfC;
2103 ws[j] = 1.0 / (p->ErrZ2());
2106 // second tep final fit
2107 Double_t s2Sum = 0.0, zSum = 0.0, szSum = 0.0, sSum = 0.0, sumw = 0.0;
2108 for (i = 0; i < dim; i++) {
2109 s2Sum += ws[i] * s[i] * s[i];
2110 zSum += ws[i] * zz[i];
2111 sSum += ws[i] * s[i];
2112 szSum += ws[i] * s[i] * zz[i];
2119 temp = s2Sum - sSum*sSum;
2122 dz = (s2Sum*zSum - sSum*szSum) / temp;
2123 tanL = (szSum - sSum*zSum) / temp;