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 // ---------------------------------------------------------------------------------
23 #include <TObjArray.h>
30 #include "AliITSgeom.h"
31 #include "AliITStrackSA.h"
32 #include "AliITSclusterV2.h"
33 #include "AliITStrackV2.h"
35 #include "AliITStrackerANN.h"
40 ClassImp(AliITStrackerANN)
42 //__________________________________________________________________________________
43 AliITStrackerANN::AliITStrackerANN(const AliITSgeom *geom, Int_t msglev)
44 : AliITStrackerV2(geom), fMsgLevel(msglev)
46 /**************************************************************************
48 CONSTRUCTOR (standard-to-use version)
51 1) the ITS geometry used in the generated event
52 2) the flag for log-messages writing
54 The AliITSgeometry is used along the class,
55 in order to translate the local AliITSclusterV2 coordinates
56 into the Global reference frame, which is necessary for the
57 Neural Network algorithm to operate.
58 In case of serialized use, the log messages should be excluded,
59 in order to save real execution time.
62 - all pointer data members are initialized
63 (if possible, otherwise are set to NULL)
64 - all numeric data members are initialized to
65 values which allow the Neural Network to operate
66 even if nothing is externally set.
68 NOTE: it is possible that tracking an event
69 with these default values results in a non-sense.
71 **************************************************************************/
76 fGeom = (AliITSgeom*)geom;
78 // Initialize the array of first module indexes per layer
79 fNLayers = geom->GetNlayers();
80 fFirstModInLayer = new Int_t[fNLayers + 1];
81 for (i = 0; i < fNLayers; i++) {
82 fFirstModInLayer[i] = fGeom->GetModuleIndex(i + 1, 1, 1);
84 fFirstModInLayer[fNLayers] = geom->GetIndexMax();
86 // initialization: no curvature cut steps
90 // initialization: 4 sectors (one for each quadrant)
92 fSectorWidth = fgkHalfPi;
94 // initialization: theta offset of 20 degrees
95 fPolarInterval = 20.0;
97 // initialization: array structure not defined
98 fStructureOK = kFALSE;
100 // initialization: vertex in the origin
105 // initialization: uninitialized point array
108 // initialization: very large (unuseful) cut values
110 for (ilayer = 0; ilayer < 6; ilayer++) {
111 fThetaCut2D[ilayer] = TMath::Pi();
112 fThetaCut3D[ilayer] = TMath::Pi();
113 fHelixMatchCut[ilayer] = 1.0;
116 // initialization: inictial activation range between 0.3 and 0.7
120 // initialization: neural network operation & weights
122 fStabThreshold = 0.001;
123 fGain2CostRatio = 1.0;
127 // initialization: uninitialized array of neurons
130 // initialization: uninitialized array of tracks
134 //__________________________________________________________________________________
135 void AliITStrackerANN::SetCuts
136 (Int_t ncurv, Double_t *curv, Double_t *theta2D, Double_t *theta3D, Double_t *helix)
138 /**************************************************************************
142 Allows for the definition of all kind of geometric cuts
143 which have been studied in for the creation of a neuron
144 from a pair of clusters C1 and C2 on consecutive layers.
145 Neuron will be created only if the pair passes ALL of these cuts.
147 At the moment, we define 4 kinds of geometrical cuts:
148 a) cut on the difference of the polar 'theta' angle;
149 b) cut on the angle between origin->C1 and C1->C2 in space;
150 c) cut on the curvature of the circle passing
151 through C1, C2 and the primary vertex;
152 d) cut on heli-matching of the same three points.
155 1) the number of curvature cut steps
156 2) the array of curvature cuts for each step
157 (its dimension is given by the argument 1)
158 3) array of 5 cut values (one for each consecutive lauer pair)
160 4) array of 5 cut values (one for each consecutive lauer pair)
162 5) array of 5 cut values (one for each consecutive lauer pair)
166 - gets the values for each cut and stores them in data members
167 - in the case of curvature cuts, the cut array
168 (whose size is not fixed) is allocated
170 NOTE: in case that the user wants to set onyl ONE of the 4 cuts array,
171 he can simply pass NULL arguments for other cuts, and (eventually)
172 a ZERO as the first argument (if curvature cuts have not to be set).
174 Anyway, all the cuts have to be set at least once.
176 **************************************************************************/
181 /*** Curvature cut setting ***/
183 // first of all, the curvature cuts are sorted in increasing order
184 // (from the smallest to the largest one)
185 Int_t *ind = new Int_t[ncurv];
186 TMath::Sort(ncurv, curv, ind, kFALSE);
187 // then, the curvature cut array is allocated and filled
188 // (a message with the list of defined cuts can be optionally shown)
189 fCurvCut = new Double_t[ncurv];
190 if (fMsgLevel >= 1) cout << "Number of curvature cuts: " << ncurv << endl;
191 for (i = 0; i < ncurv; i++) {
192 fCurvCut[i] = curv[ind[i]];
193 if (fMsgLevel >= 1) cout << " - " << fCurvCut[i] << endl;
197 /*** 'Fixed' cuts setting ***/
199 // checks what cuts have to be set
200 Bool_t doTheta2D = (theta2D != 0);
201 Bool_t doTheta3D = (theta3D != 0);
202 Bool_t doHelix = (helix != 0);
203 // sets the cuts for all layer pairs
204 for (i = 0; i < fNLayers - 1; i++) {
205 if (doTheta2D) fThetaCut2D[i] = theta2D[i];
206 if (doTheta3D) fThetaCut3D[i] = theta3D[i];
207 if (doHelix) fHelixMatchCut[i] = helix[i];
209 // if required, lists the cuts
210 if (!fMsgLevel < 2) return;
211 cout << "Theta 2D cuts: ";
213 cout << "<not set>" << endl;
217 for (i = 0; i < fNLayers - 1; i++) {
218 cout << "For layers " << i << " --> " << i + 1;
219 cout << " cut = " << fThetaCut2D[i] << endl;
222 cout << "---" << endl;
223 cout << "Theta 3D cuts: ";
225 cout << "<not set>" << endl;
229 for (i = 0; i < fNLayers - 1; i++) {
230 cout << "For layers " << i << " --> " << i + 1;
231 cout << " cut = " << fThetaCut3D[i] << endl;
234 cout << "---" << endl;
235 cout << "Helix-match cuts: ";
237 cout << "<not set>" << endl;
241 for (i = 0; i < fNLayers - 1; i++) {
242 cout << "For layers " << i << " --> " << i + 1;
243 cout << " cut = " << fHelixMatchCut[i] << endl;
246 cout << "---" << endl;
249 //__________________________________________________________________________________
250 Bool_t AliITStrackerANN::GetGlobalXYZ
252 Double_t &x, Double_t &y, Double_t &z, Double_t &ex2, Double_t &ey2, Double_t &ez2)
254 /**************************************************************************
256 LOCAL TO GLOBAL TRANSLATOR
258 Taking information from the ITS geometry stored in the class,
259 gets a stored AliITScluster and calculates it coordinates
260 and errors in the global reference frame.
261 These values are stored in the variables,
262 which are passed by reference.
265 1) reference index for the cluster to use
266 2) (by reference) place to store the global X-coord into
267 3) (by reference) place to store the global Y-coord into
268 4) (by reference) place to store the global Z-coord into
269 5) (by reference) place to store the global X-coord error into
270 6) (by reference) place to store the global Y-coord error into
271 7) (by reference) place to store the global Z-coord error into
274 essentially, determines the ITS module index from the
275 detector index of the AliITSclusterV2 object, and extracts
276 the roto-translation from the ITS geometry, to convert
277 the local module coordinates into the global ones.
280 - kFALSE if the given cluster index points to a non-existing
281 cluster, or if the layer number makes no sense (< 0 or > 6).
282 - otherwise, kTRUE (meaning a successful operation).
284 **************************************************************************/
286 // checks if the layer is correct
287 Int_t ilayer = (refIndex & 0xf0000000) >> 28;
288 if (ilayer < 0 || ilayer >= fNLayers) {
289 Error("GetGlobalXYZ", "Wrong layer number: %d [range: %d - %d]", ilayer, 0, fNLayers);
292 // checks if the referenced cluster exists and corresponds to the passed reference
293 AliITSclusterV2 *refCluster = (AliITSclusterV2*) GetCluster(refIndex);
295 Error("GetGlobalXYZ", "Cluster not found for index %d", refIndex);
299 // determine the detector number
300 Int_t detID = refCluster->GetDetectorIndex() + fFirstModInLayer[ilayer];
302 // get rotation matrix
304 fGeom->GetRotMatrix(detID, rot);
306 // get translation vector
308 fGeom->GetTrans(detID, tx, ty, tz);
310 // determine r and phi for the reference conversion
311 Double_t r = -(Double_t)tx * rot[1] + (Double_t)ty * rot[0];
312 if (ilayer == 0) r = -r;
313 Double_t phi = TMath::ATan2(rot[1],rot[0]);
314 if (ilayer == 0) phi -= fgkPi;
316 // sets values for X, Y, Z in global coordinates and their errors
317 Double_t cosPhi = TMath::Cos(phi);
318 Double_t sinPhi = TMath::Sin(phi);
319 x = r*cosPhi + refCluster->GetY()*sinPhi;
320 y = -r*sinPhi + refCluster->GetY()*cosPhi;
321 z = refCluster->GetZ();
322 ex2 = refCluster->GetSigmaY2()*sinPhi*sinPhi;
323 ey2 = refCluster->GetSigmaY2()*cosPhi*cosPhi;
324 ez2 = refCluster->GetSigmaZ2();
329 //__________________________________________________________________________________
330 AliITStrackerANN::AliITSnode* AliITStrackerANN::AddNode(Int_t refIndex)
332 /**************************************************************************
334 GENERATOR OF NEURAL NODES
336 Fills the array of neural 'nodes', which are the ITS clusters
337 translated in the global reference frame.
338 Given that the global coordinates are used many times, they are
339 stored in a well-defined structure, in the form of an embedded class.
340 Moreover, this class allows a faster navigation among points
341 and neurons, by means of some object arrays, storing only the
342 neurons which start from, or end to, the given node.
343 Finally, each node contains all the other nodes which match it
344 with respect to the fixed walues, in order to perform a faster
345 neuron-creation phase.
348 1) reference index of the correlated AliITSclusterV2 object
351 - allocates the new AliITSnode objects
352 - initializes its object arrays
353 - from the global coordinates, calculates the
354 'phi' and 'theta' coordinates, in order to store it
355 into the correct theta-slot and azimutal sector.
358 - the pointer of the creater AliITSnode object
359 - in case of errors, a waring is given and a NULL is returned
361 **************************************************************************/
363 // create object and set the reference
364 AliITSnode *node = new AliITSnode;
366 Warning("AddNode", "Error occurred when allocating AliITSnode");
369 node->ClusterRef() = refIndex;
371 // calls the conversion function, which makes also some checks
372 // (layer number within range, existence of referenced cluster)
375 node->X(), node->Y(), node->Z(),
376 node->ErrX2(), node->ErrY2(), node->ErrZ2()
379 // initializes the object arrays
380 node->Matches() = new TObjArray;
381 node->InnerOf() = new TObjArray;
382 node->OuterOf() = new TObjArray;
384 // finds azimutal and polar sector (in degrees)
385 Double_t phi = node->GetPhi() * 180.0 / fgkPi;
386 Double_t theta = node->GetTheta() * 180.0 / fgkPi;
387 Int_t isector = (Int_t)(phi / fSectorWidth);
388 Int_t itheta = (Int_t)theta;
389 Int_t ilayer = (refIndex & 0xf0000000) >> 28;
391 // selects the right TObjArray to store object into
392 TObjArray *sector = (TObjArray*)fNodes[ilayer][itheta]->At(isector);
393 sector->AddLast(node);
398 //__________________________________________________________________________________
399 void AliITStrackerANN::CreateArrayStructure(Int_t nsectors)
401 /**************************************************************************
403 ARRAY STRUCTURE CREATOR
405 Creates a structure of nested TObjArray's where the AliITSnode's
407 - the first level is made by 6 arrays (one for each layer)
408 - the second level is made by 180 arrays (one for each int part of theta)
409 - the third level is made by a variable number of arrays
410 (one for each azimutal sector)
413 1) the number of azimutal sectors
416 - calculates the width of each sector, from the argument
417 - allocates and initializes all array levels
418 - sets a flag which tells the user if this NECESSARY operation
419 has been performed (it is needed BEFORE performing tracking)
421 **************************************************************************/
423 // Set the number of sectors and their width.
424 fSectorNum = nsectors;
425 fSectorWidth = 360.0 / (Double_t)fSectorNum;
426 if (fMsgLevel >= 2) {
427 cout << fSectorNum << " sectors --> sector width (degrees) = " << fSectorWidth << endl;
430 // Meaningful indexes
431 Int_t ilayer, isector, itheta;
433 // Mark for the created objects
434 TObjArray *sector = 0;
436 // First index: layer
437 fNodes = new TObjArray**[fNLayers];
438 for (ilayer = 0; ilayer < fNLayers; ilayer++) {
439 fNodes[ilayer] = new TObjArray*[180];
440 for (itheta = 0; itheta < 180; itheta++) fNodes[ilayer][itheta] = 0;
441 for (itheta = 0; itheta < 180; itheta++) {
442 fNodes[ilayer][itheta] = new TObjArray(nsectors);
443 for (isector = 0; isector < nsectors; isector++) {
444 sector = new TObjArray;
446 fNodes[ilayer][itheta]->AddAt(sector, isector);
451 // Sets a checking flag to TRUE.
452 // This flag is checked before filling up the arrays with the points.
453 fStructureOK = kTRUE;
456 //__________________________________________________________________________________
457 Int_t AliITStrackerANN::ArrangePoints(char *exportFile)
459 /**************************************************************************
463 This function assembles the operation from the other above methods,
464 and fills the arrays with the clusters already stored in the
465 layers of the tracker.
466 Then, in order to use this method, the user MUSTs call LoadClusters()
470 1) string for a file name where the global ccordinates
471 of all points can be exported (optional).
472 If this file must not be created, simply pass a NULL argument
475 - for each AliITSclusterV2 in each AliITSlayer, a ne AliITSnode
476 is created and stored in the correct location.
479 - the number of stored points
480 - when errors occur, or no points are found, 0 is returned
482 **************************************************************************/
484 // Check if the array structure has been created
486 Error("ArrangePoints", "Structure NOT defined. Call CreateArrayStructure() first");
490 // meaningful indexes
491 Int_t ientry, ilayer, nentries = 0, index;
494 // if the argument is not NULL, a file is opened
495 fstream file(exportFile, ios::out);
496 if (!exportFile || file.fail()) {
501 // scan all layers for node creation
502 for (ilayer = 0; ilayer < fNLayers; ilayer++) {
503 nPtsLayer = GetNumberOfClustersLayer(ilayer);
504 if (fMsgLevel >= 1) {
505 cout << "Layer " << ilayer << " --> " << nPtsLayer << " clusters" << endl;
507 for (ientry = 0; ientry < nPtsLayer; ientry++) {
508 // calculation of cluster index : (Bit mask LLLLIIIIIIIIIIII)
509 // (L = bits used for layer)
510 // (I = bits used for position in layer)
511 index = ilayer << 28;
513 // add new AliITSnode object
514 AliITSnode *n = AddNode(index);
515 if ( (n != NULL) && exportFile ) {
516 file << index << ' ' << n->X() << ' ' << n->Y() << ' ' << n->Z() << endl;
519 nentries += nPtsLayer;
522 // a conventional final message is put at the end of file
524 file << "-1 0.0 0.0 0.0" << endl;
528 // returns the number of points processed
532 //__________________________________________________________________________________
533 void AliITStrackerANN::StoreOverallMatches()
535 /**************************************************************************
539 Once the nodes have been created, a firs analysis is to check
540 what pairs will satisfy at least the 'fixed' cuts (theta, helix-match)
541 and the most permissive (= larger) curvature cut.
542 All these node pairs are suitable for neuron creation.
543 In fact, when performing a Neural Tracking step, the only further check
544 will be a check against the current curvature step, while the other
546 For thi purpose, each AliITSnode has a data member, named 'fMatches'
547 which contains references to all other AliITSnodes in the successive layer
548 that form, with it, a 'good' pair, with respect to the above cited cuts.
549 Then, in each step for neuron creation, the possible neurons starting from
550 each node will be searched ONLY within the nodes referenced in fMatches.
551 This, of course, speeds up a lot the neuron creation procedure, at the
552 cost of some memory occupation, which results not to be critical.
555 - for each AliITSnode, matches are found according to the criteria
556 expressed above, and stored in the node->fMatches array
558 **************************************************************************/
560 // meaningful counters
561 Int_t ilayer, isector, itheta1, itheta2, check;
562 TObjArray *list1 = 0, *list2 = 0;
563 AliITSnode *node1 = 0, *node2 = 0;
564 Double_t thetaMin, thetaMax;
567 // Scan for each sector
568 for (isector = 0; isector < fSectorNum; isector++) {
569 // sector is chosen once for both lists
570 for (ilayer = 0; ilayer < fNLayers - 1; ilayer++) {
571 for (itheta1 = 0; itheta1 < 180; itheta1++) {
572 list1 = (TObjArray*)fNodes[ilayer][itheta1]->At(isector);
573 TObjArrayIter iter1(list1);
574 while ( (node1 = (AliITSnode*)iter1.Next()) ) {
575 if (node1->IsUsed()) continue;
576 // clear an eventually already present array
577 // node1->Matches()->Clear();
578 // get the global coordinates and defines the theta interval from cut
579 thetaMin = (node1->GetTheta() * 180.0 / fgkPi) - fPolarInterval;
580 thetaMax = (node1->GetTheta() * 180.0 / fgkPi) + fPolarInterval;
581 imin = (Int_t)thetaMin;
582 imax = (Int_t)thetaMax;
583 if (imin < 0) imin = 0;
584 if (imax > 179) imax = 179;
585 // loop on the selected theta slots
586 for (itheta2 = imin; itheta2 <= imax; itheta2++) {
587 list2 = (TObjArray*)fNodes[ilayer + 1][itheta2]->At(isector);
588 TObjArrayIter iter2(list2);
589 while ( (node2 = (AliITSnode*)iter2.Next()) ) {
590 check = PassAllCuts(node1, node2, fCurvNum - 1, fVertexX, fVertexY, fVertexZ);
592 node1->Matches()->AddLast(node2);
594 } // while (node2...)
595 } // for (itheta2...)
596 } // while (node1...)
599 } // for (isector...)
602 //__________________________________________________________________________________
603 Int_t AliITStrackerANN::PassAllCuts
604 (AliITSnode *inner, AliITSnode *outer, Int_t curvStep,
605 Double_t vx, Double_t vy, Double_t vz)
607 // ***********************************************************************************
609 // This check is called in the above method for finding the matches of each node
610 // It check the passed point pair against all the fixed cuts and a specified
611 // curvature cut, among all the ones which have been defined.
612 // The cuts need a vertex-constraint, which is not absolute, but it is passed
616 // 1) the point in the inner layer
617 // 2) the point in the outer layer
618 // 3) curvature step for the curvature cut check (preferably the last)
619 // 4) X of the used vertex
620 // 5) Y of the used vertex
621 // 6) Z of the used vertex
624 // - if necessary, swaps the two points
625 // (the first must be in the innermost of the two layers)
626 // - checks for the theta cuts
627 // - calculates the circle passing through the vertex
628 // and the given points and checks for the curvature cut
629 // - using the radius calculated there, checks for the helix-math cut
632 // 0 - All cuts passed
633 // 1 - theta 2D cut not passed
634 // 2 - theta 3D cut not passed
635 // 3 - curvature calculated but cut not passed
636 // 4 - curvature not calculated (division by zero)
637 // 5 - helix cut not passed
638 // 6 - curvature inxed out of range
640 // ***********************************************************************************
642 // Check for curvature index
643 if (curvStep < 0 || curvStep >= fCurvNum) return 6;
645 // Swap points in order that r1 < r2
646 AliITSnode *temp = 0;
647 if (outer->GetLayer() < inner->GetLayer()) {
653 // All cuts are variable according to the layer of the
654 // innermost point (the other point will surely be
655 // in the further one, because we don't check poin pairs
656 // which are not in adjacent layers)
657 // The reference is given by the innermost point.
658 Int_t layRef = inner->GetLayer();
660 // The calculations in the transverse plane are made in
661 // a shifted reference frame, whose origin corresponds to
662 // the reference point passed in the argument.
663 Double_t xIn = inner->X() - vx;
664 Double_t xOut = outer->X() - vx;
665 Double_t yIn = inner->Y() - vy;
666 Double_t yOut = outer->Y() - vy;
667 Double_t zIn = inner->Z() - vz;
668 Double_t zOut = outer->Z() - vz;
669 Double_t rIn = TMath::Sqrt(xIn*xIn + yIn*yIn);
670 Double_t rOut = TMath::Sqrt(xOut*xOut + yOut*yOut);
672 // Check for theta cut.
673 // There are two different angular cuts:
674 // one w.r. to the angle in the 2-dimensional r-z plane...
676 TVector3 origin2innerRZ(zIn, rIn, 0.0);
677 TVector3 inner2outerRZ(zOut - zIn, rOut - rIn, 0.0);
678 dthetaRZ = origin2innerRZ.Angle(inner2outerRZ) * 180.0 / fgkPi;
679 if (dthetaRZ > fThetaCut2D[layRef]) {
681 // r-z theta cut not passed ---> 1
683 // ...and another w.r. to the angle in the 3-dimensional x-y-z space
685 TVector3 origin2innerXYZ(xIn, yIn, zIn);
686 TVector3 inner2outerXYZ(xOut - xIn, yOut - yIn, zOut - zIn);
687 dthetaXYZ = origin2innerXYZ.Angle(inner2outerXYZ) * 180.0 / fgkPi;
688 if (dthetaXYZ > fThetaCut3D[layRef]) {
690 // x-y-z theta cut not passed ---> 2
693 // Calculation & check of curvature
694 Double_t dx = xIn - xOut;
695 Double_t dy = yIn - yOut;
696 Double_t num = 2.0 * (xIn*yOut - xOut*yIn);
697 Double_t den = rIn*rOut*sqrt(dx*dx + dy*dy);
700 curv = TMath::Abs(num / den);
701 if (curv > fCurvCut[curvStep]) {
703 // curvature too large for cut ---> 3
707 Error("PassAllCuts", "Curvature calculations gives zero denominator");
709 // error: denominator = 0 ---> 4
712 // Calculation & check of helix matching
713 Double_t helMatch = 0.0;
714 Double_t arcIn = 2.0 * rIn * curv;
715 Double_t arcOut = 2.0 * rOut * curv;
716 if (arcIn > -1.0 && arcIn < 1.0)
717 arcIn = TMath::ASin(arcIn);
719 arcIn = ((arcIn > 0.0) ? 0.5 : 1.5) * TMath::Pi();
720 if (arcOut > -1.0 && arcOut < 1.0)
721 arcOut = TMath::ASin(arcOut);
723 arcOut = ((arcOut > 0.0) ? 0.5 : 1.5) * TMath::Pi();
725 arcOut /= 2.0 * curv;
726 if (arcIn == 0.0 || arcOut == 0.0) {
727 Error("PassAllCuts", "Calculation returns zero-length arcs: l1=%f, l2=%f", arcIn, arcOut);
729 // error: circumference arcs seem to equal zero ---> 4
731 helMatch = TMath::Abs(zIn / arcIn - zOut / arcOut);
732 if (helMatch > fHelixMatchCut[layRef]) {
734 // helix match cut not passed ---> 5
737 // ALL CUTS PASSED ---> 0
741 //__________________________________________________________________________________
742 Bool_t AliITStrackerANN::PassCurvCut
743 (AliITSnode *inner, AliITSnode *outer,
745 Double_t vx, Double_t vy, Double_t vz)
747 //***********************************************************************************
749 // This method operates essentially like the above one, but it is used
750 // during a single step of Neural Tracking, where the curvature cut
752 // Then, not necessaryly all the nodes stored in the fMatches array
753 // will be suitable for neuron creation in an intermediate step.
755 // It has the same arguments of the PassAllCuts() method, but
756 // the theta cut is not checked.
757 // Moreover, it has a boolean return value.
759 //***********************************************************************************
761 // Check for curvature index
762 if (curvStep < 0 || curvStep >= fCurvNum) return 6;
764 // Find the reference layer
765 Int_t layIn = inner->GetLayer();
766 Int_t layOut = outer->GetLayer();
767 Int_t layRef = (layIn < layOut) ? layIn : layOut;
769 // The calculations in the transverse plane are made in
770 // a shifted reference frame, whose origin corresponds to
771 // the reference point passed in the argument.
772 Double_t xIn = inner->X() - vx;
773 Double_t xOut = outer->X() - vx;
774 Double_t yIn = inner->Y() - vy;
775 Double_t yOut = outer->Y() - vy;
776 Double_t zIn = inner->Z() - vz;
777 Double_t zOut = outer->Z() - vz;
778 Double_t rIn = TMath::Sqrt(xIn*xIn + yIn*yIn);
779 Double_t rOut = TMath::Sqrt(xOut*xOut + yOut*yOut);
781 // Calculation & check of curvature
782 Double_t dx = xIn - xOut;
783 Double_t dy = yIn - yOut;
784 Double_t num = 2.0 * (xIn*yOut - xOut*yIn);
785 Double_t den = rIn*rOut*sqrt(dx*dx + dy*dy);
789 curv = TMath::Abs(num / den);
790 if (curv > fCurvCut[curvStep]) return kFALSE;
798 curv = TMath::Abs(num / den);
799 if (curv > fCurvCut[curvStep]) {
804 Error("PassAllCuts", "Curvature calculations gives zero denominator");
808 // Calculation & check of helix matching
809 Double_t helMatch = 0.0;
810 Double_t arcIn = 2.0 * rIn * curv;
811 Double_t arcOut = 2.0 * rOut * curv;
812 if (arcIn > -1.0 && arcIn < 1.0)
813 arcIn = TMath::ASin(arcIn);
815 arcIn = ((arcIn > 0.0) ? 0.5 : 1.5) * TMath::Pi();
816 if (arcOut > -1.0 && arcOut < 1.0)
817 arcOut = TMath::ASin(arcOut);
819 arcOut = ((arcOut > 0.0) ? 0.5 : 1.5) * TMath::Pi();
821 arcOut /= 2.0 * curv;
822 if (arcIn == 0.0 || arcOut == 0.0) {
823 Error("PassAllCuts", "Calculation returns zero-length arcs: l1=%f, l2=%f", arcIn, arcOut);
825 // error: circumference arcs seem to equal zero ---> 4
827 helMatch = TMath::Abs(zIn / arcIn - zOut / arcOut);
828 return (helMatch <= fHelixMatchCut[layRef]);
831 //__________________________________________________________________________________
832 void AliITStrackerANN::PrintMatches(Bool_t stop)
834 // Prints the list of points which appear to match
835 // each one of them, according to the preliminary
837 // The arguments states if a pause is required after printing
838 // the matches for each one. In this case, a keypress is needed.
840 TObjArray *sector = 0;
841 Int_t ilayer, isector, itheta, nF;
842 AliITSnode *node1 = 0, *node2 = 0;
843 //AliITSclusterV2 *cluster1 = 0, *cluster2 = 0;
845 for (ilayer = 0; ilayer < 6; ilayer++) {
846 for (isector = 0; isector < fSectorNum; isector++) {
847 for (itheta = 0; itheta < 180; itheta++) {
848 sector = (TObjArray*)fNodes[ilayer][itheta]->At(isector);
849 TObjArrayIter points(sector);
850 while ( (node1 = (AliITSnode*)points.Next()) ) {
851 nF = (Int_t)node1->Matches()->GetEntries();
852 cout << "Node layer: " << node1->GetLayer() << " --> ";
854 cout << "NO Matches!!!" << endl;
857 cout << nF << " Matches" << endl;
858 cout << "Reference cluster: " << hex << node1->ClusterRef() << endl;
859 TObjArrayIter matches(node1->Matches());
860 while ( (node2 = (AliITSnode*)matches.Next()) ) {
861 cout << "Match with " << hex << node2->ClusterRef() << endl;
864 cout << "Press a key" << endl;
873 //__________________________________________________________________________________
874 void AliITStrackerANN::ResetNodes(Int_t isector)
876 /***********************************************************************************
878 NODE NEURON ARRAY CLEANER
880 After a neural tracking step, this method
881 clears the arrays 'fInnerOf' and 'fOuterOf' of each AliITSnode
884 - the sector where the operation is being executed
886 ***********************************************************************************/
888 Int_t ilayer, itheta;
889 TObjArray *sector = 0;
890 AliITSnode *node = 0;
891 for (ilayer = 0; ilayer < fNLayers; ilayer++) {
892 for (itheta = 0; itheta < 180; itheta++) {
893 sector = (TObjArray*)fNodes[ilayer][itheta]->At(isector);
894 TObjArrayIter iter(sector);
896 node = (AliITSnode*)iter.Next();
898 node->InnerOf()->Clear();
899 node->OuterOf()->Clear();
901 delete node->InnerOf();
902 delete node->OuterOf();
903 node->InnerOf() = new TObjArray;
904 node->OuterOf() = new TObjArray;
911 //__________________________________________________________________________________
912 Int_t AliITStrackerANN::CreateNeurons
913 (AliITSnode *node, Int_t curvStep, Double_t vx, Double_t vy, Double_t vz)
915 // This method is used to create alle suitable neurons starting from
916 // a single AliITSnode. Each unit is also stored in the fInnerOf array
917 // of the passed node, and in the fOuterOf array of the other neuron edge.
918 // In the new implementation of the intermediate check steps, a further one
919 // is made, which chechs how well a helix is matched by three points
920 // in three consecutive layers.
921 // Then, a vertex-constrained check is made with vertex located
922 // in a layer L, for points in layers L+1 and L+2.
924 // In order to do this, the creator works recursively, in a tree-visit like operation.
925 // The neurons are effectively created only if the node argument passed is in
926 // the 5th layer (they are created between point of 5th and 6th layer).
927 // If the node is in an inner layer, its coordinates are passet as vertex for a nested
928 // call of the same function in the next two layers.
932 // 2) current curvature cut step
933 // 3) X of referenced temporary (not primary) vertex
934 // 4) Y of referenced temporary (not primary) vertex
935 // 5) Z of referenced temporary (not primary) vertex
938 // - if the layer is the 5th, neurons are created with nodes
939 // in the fMatches array of the passed node
940 // - otherwise, the X, Y, Z of the passed node are given as
941 // vertex and the same procedure is recursively called for all
942 // nodes in the fMatches array of the passed one.
945 // - the total number of neurons created from the passed one
946 // summed with all neurons created from all nodes well matched with it
947 // (assumes a meaning only for nodes in the first layer)
950 Int_t made = 0; // counter
951 Bool_t found = 0; // flag
952 AliITSnode *match = 0; // cursor for a AliITSnode array
953 AliITSneuron *unit = 0; // cursor for a AliITSneuron array
955 // --> Case 0: the passed node has already been used
956 // as member of a track found in a previous step.
957 // In this case, of course, the function exits.
958 if (node->IsUsed()) return 0;
960 // --> Case 1: there are NO well-matched points.
961 // This can happen in all ITS layers, but it happens **for sure**
962 // for a node in the 'last' layer.
963 // Even in this case, the function exits.
964 if (node->Matches()->IsEmpty()) return 0;
966 // --> Case 2: there are well-matched points.
967 // In this case, the function creates a neuron for each
968 // well-matched pair (according to the cuts for the current step)
969 // Moreover, before storing the neuron, a check is necessary
970 // to avoid the duplicate creation of the same neuron twice.
971 // (This could happen if the 3 last arguments of the function
972 // are close enough to cause a good match for the current step
973 // between two points, independently of their difference).
974 // Finally, a node is skipped if it has already been used.
975 // For each matched point for which a neuron is created, the procedure is
976 // recursively called.
977 TObjArrayIter matches(node->Matches());
978 while ( (match = (AliITSnode*)matches.Next()) ) {
979 if (match->IsUsed()) continue;
980 if (!PassCurvCut(node, match, curvStep, vx, vy, vz)) continue;
982 if (!node->InnerOf()->IsEmpty()) {
983 TObjArrayIter presentUnits(node->InnerOf());
984 while ( (unit = (AliITSneuron*)presentUnits.Next()) ) {
985 if (unit->Inner() == node && unit->Outer() == match) {
992 AliITSneuron *unit = new AliITSneuron(node, match, fEdge2, fEdge1);
993 fNeurons->AddLast(unit);
994 node->InnerOf()->AddLast(unit);
995 match->OuterOf()->AddLast(unit);
996 made += CreateNeurons(match, curvStep, node->X(), node->Y(), node->Z());
1000 // Of course, the return value contains the number of neurons
1001 // counting in also the oned created in all levels of recursive calls.
1005 //__________________________________________________________________________________
1006 Int_t AliITStrackerANN::CreateNetwork(Int_t sector, Int_t curvStep)
1008 // This function simply recalls the CreateNeurons() method for each node
1009 // in the first layer, for the current sector.
1010 // This generates the whole network, thanks to the recursive calls.
1013 // 1) current sector
1014 // 2) current curvature step
1017 // - scans the nodes array for all theta's in the current sector
1018 // and layer 0, and calls the CreateNeurons() function.
1020 // removes all eventually present neurons
1021 if (fNeurons) delete fNeurons;
1022 fNeurons = new TObjArray;
1023 fNeurons->SetOwner(kTRUE);
1025 // calls the ResetNodes() function to free the AliITSnode arrays
1026 if (fMsgLevel >= 2) {
1027 cout << "Sector " << sector << " PHI = ";
1028 cout << fSectorWidth * (Double_t)sector << " --> ";
1029 cout << fSectorWidth * (Double_t)(sector + 1) << endl;
1030 cout << "Curvature step " << curvStep << " [cut = " << fCurvCut[curvStep] << "]" << endl;
1034 // meaningful counters
1035 Int_t itheta, neurons = 0;
1036 TObjArray *lstSector = 0;
1039 Double_t vx[6], vy[6], vz[6];
1040 AliITSnode *p[6] = {0, 0, 0, 0, 0, 0};
1041 for (itheta = 0; itheta < 180; itheta++) {
1042 lstSector = (TObjArray*)fNodes[0][itheta]->At(sector);
1043 TObjArrayIter lay0(lstSector);
1044 while ( (p[0] = (AliITSnode*)lay0.Next()) ) {
1045 if (p[0]->IsUsed()) continue;
1049 neurons += CreateNeurons(p[0], curvStep, fVertexX, fVertexY, fVertexZ);
1051 TObjArrayIter lay1(p[0]->Matches());
1052 while ( (p[1] = (AliITSnode*)lay1.Next()) ) {
1053 if (p[1]->IsUsed()) continue;
1054 if (!PassCurvCut(p[0], p[1], curvStep, vx[0], vy[0], vz[0])) continue;
1055 unit = new AliITSneuron;
1056 unit->Inner() = p[0];
1057 unit->Outer() = p[1];
1058 unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1059 unit->fGain = new TObjArray;
1060 fNeurons->AddLast(unit);
1061 p[0]->InnerOf()->AddLast(unit);
1062 p[1]->OuterOf()->AddLast(unit);
1067 TObjArrayIter lay2(p[1]->Matches());
1068 while ( (p[2] = (AliITSnode*)lay2.Next()) ) {
1069 if (p[2]->IsUsed()) continue;
1070 if (!PassCurvCut(p[1], p[2], curvStep, vx[1], vy[1], vz[1])) continue;
1071 unit = new AliITSneuron;
1072 unit->Inner() = p[1];
1073 unit->Outer() = p[2];
1074 unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1075 unit->fGain = new TObjArray;
1076 fNeurons->AddLast(unit);
1077 p[1]->InnerOf()->AddLast(unit);
1078 p[2]->OuterOf()->AddLast(unit);
1083 TObjArrayIter lay3(p[2]->Matches());
1084 while ( (p[3] = (AliITSnode*)lay3.Next()) ) {
1085 if (p[3]->IsUsed()) continue;
1086 if (!PassCurvCut(p[2], p[3], curvStep, vx[2], vy[2], vz[2])) continue;
1087 unit = new AliITSneuron;
1088 unit->Inner() = p[2];
1089 unit->Outer() = p[3];
1090 unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1091 unit->fGain = new TObjArray;
1092 fNeurons->AddLast(unit);
1093 p[2]->InnerOf()->AddLast(unit);
1094 p[3]->OuterOf()->AddLast(unit);
1099 TObjArrayIter lay4(p[3]->Matches());
1100 while ( (p[4] = (AliITSnode*)lay4.Next()) ) {
1101 if (p[4]->IsUsed()) continue;
1102 if (!PassCurvCut(p[3], p[4], curvStep, vx[3], vy[3], vz[3])) continue;
1103 unit = new AliITSneuron;
1104 unit->Inner() = p[3];
1105 unit->Outer() = p[4];
1106 unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1107 unit->fGain = new TObjArray;
1108 fNeurons->AddLast(unit);
1109 p[3]->InnerOf()->AddLast(unit);
1110 p[4]->OuterOf()->AddLast(unit);
1115 TObjArrayIter lay5(p[4]->Matches());
1116 while ( (p[5] = (AliITSnode*)lay5.Next()) ) {
1117 if (p[5]->IsUsed()) continue;
1118 if (!PassCurvCut(p[4], p[5], curvStep, vx[4], vy[4], vz[4])) continue;
1119 unit = new AliITSneuron;
1120 unit->Inner() = p[4];
1121 unit->Outer() = p[5];
1122 unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1123 unit->fGain = new TObjArray;
1124 fNeurons->AddLast(unit);
1125 p[4]->InnerOf()->AddLast(unit);
1126 p[5]->OuterOf()->AddLast(unit);
1135 } // for (itheta...)
1136 // END OF NEW VERSION
1139 for (ilayer = 0; ilayer < 6; ilayer++) {
1140 for (itheta = 0; itheta < 180; itheta++) {
1141 lstSector = (TObjArray*)fNodes[ilayer][itheta]->At(sector_idx);
1142 TObjArrayIter inners(lstSector);
1143 while ( (inner = (AliITSnode*)inners.Next()) ) {
1144 if (inner->GetUser() >= 0) continue;
1145 TObjArrayIter outers(inner->Matches());
1146 while ( (outer = (AliITSnode*)outers.Next()) ) {
1147 if (outer->GetUser() >= 0) continue;
1148 if (!PassCurvCut(inner, outer, curvStep, fVX, fVY, fVZ)) continue;
1149 unit = new AliITSneuron;
1150 unit->Inner() = inner;
1151 unit->Outer() = outer;
1152 unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
1153 unit->fGain = new TObjArray;
1154 fNeurons->AddLast(unit);
1155 inner->InnerOf()->AddLast(unit);
1156 outer->OuterOf()->AddLast(unit);
1160 } // for (itheta...)
1161 } // for (ilayer...)
1164 fNeurons->SetOwner();
1168 //__________________________________________________________________________________
1169 Int_t AliITStrackerANN::LinkNeurons()
1171 /***********************************************************************************
1175 Scans the whole neuron array, in order to find all neuron pairs
1176 which are connected in sequence and share a positive weight.
1177 For each of them, an AliITSlink is created, which stores
1178 the weight value, and will allow for a faster calculation
1179 of the total neural input for each updating cycle.
1181 Every neuron contains an object array which stores all AliITSlink
1182 objects which point to sequenced units, with the respective weights.
1185 - the number of link created for the neural network.
1186 If they are 0, no updating can be done and the step is skipped.
1188 ***********************************************************************************/
1190 // meaningful indexes
1192 Double_t weight = 0.0;
1193 TObjArrayIter neurons(fNeurons), *iter;
1194 AliITSneuron *neuron = 0, *test = 0;
1196 // scan in the neuron array
1198 neuron = (AliITSneuron*)neurons.Next();
1200 // checks for neurons startin from its head ( -> )
1201 iter = (TObjArrayIter*)neuron->Inner()->OuterOf()->MakeIterator();
1203 test = (AliITSneuron*)iter->Next();
1205 weight = Weight(test, neuron);
1206 if (weight > 0.0) neuron->Gain()->AddLast(new AliITSlink(weight, test));
1210 // checks for neurons ending in its tail ( >- )
1211 iter = (TObjArrayIter*)neuron->Outer()->InnerOf()->MakeIterator();
1213 test = (AliITSneuron*)iter->Next();
1215 weight = Weight(neuron, test);
1216 if (weight > 0.0) neuron->Gain()->AddLast(new AliITSlink(weight, test));
1224 //__________________________________________________________________________________
1225 Bool_t AliITStrackerANN::Update()
1227 /***********************************************************************************
1229 Performs a single updating cycle.
1232 - for each neuron, gets the activation with the neuron Activate() method
1233 - checks if stability has been reached (compare mean activation variation
1234 with the stability threshold data member)
1237 - kTRUE means that the neural network has stabilized
1238 - kFALSE means that another updating cycle is needed
1240 ***********************************************************************************/
1242 Double_t actVar = 0.0, totDiff = 0.0;
1243 TObjArrayIter iter(fNeurons);
1246 unit = (AliITSneuron*)iter.Next();
1248 actVar = unit->Activate(fTemperature);
1249 // calculation the relative activation variation
1252 totDiff /= fNeurons->GetSize();
1253 return (totDiff < fStabThreshold);
1256 //__________________________________________________________________________________
1257 void AliITStrackerANN::FollowChains(Int_t sector)
1259 /***********************************************************************************
1263 After that the neural network has stabilized,
1264 the final step is to create polygonal chains
1265 of clusters, one in each layer, which represent
1266 the tracks recognized by the neural algorithm.
1267 This is made by means of a choice of the best
1268 neuron among the ones starting from each point.
1270 Once that such neuron is selected, its inner point
1271 will set the 'fNext' field to its outer point, and
1272 similarly, its outer point will set the 'fPrev' field
1274 This defines a bi-directional sequence.
1276 In this procedure, it can happen that many neurons
1277 which have the head of the arrow in a given node, will
1278 all select as best following the neuron with the largest
1279 activation starting in that point.
1280 This results in MANY nodes which have the same 'fNext'.
1281 But, this field will be set to NULL for all these points,
1282 but the only one which is pointed by the 'fPrev' field
1283 of this shared node.
1285 ***********************************************************************************/
1287 // meaningful counters
1288 Int_t itheta, ilayer;
1289 TObjArray *lstSector = 0;
1290 Double_t test = fActMinimum;
1292 AliITSneuron *n = 0;
1294 // scan the whole collection of nodes
1295 for (ilayer = 0; ilayer < fNLayers; ilayer++) {
1296 for (itheta = 0; itheta < 180; itheta++) {
1297 // get the array of point in a given layer/theta-slot/sector
1298 lstSector = (TObjArray*)fNodes[ilayer][itheta]->At(sector);
1299 TObjArrayIter nodes(lstSector);
1300 while ( (p = (AliITSnode*)nodes.Next()) ) {
1301 // if the point is used, it is skipped
1302 if (p->IsUsed()) continue;
1303 // initially, fNext points to nothing, and
1304 // the comparison value is set to the minimum activation
1305 // which allows to say that a neuron is turned 'on'
1306 // a node from which only 'off' neurons start is probably
1307 // a noise point, which will be excluded from the reconstruction.
1310 TObjArrayIter innerof(p->InnerOf());
1311 while ( (n = (AliITSneuron*)innerof.Next()) ) {
1312 // if the examined neuron has not the largest activation
1313 // it is skipped and removed from array of all neurons
1314 // and of its outer point (its inner is the cursor p)
1315 if (n->Activation() < test) {
1316 p->InnerOf()->Remove(n);
1317 n->Outer()->OuterOf()->Remove(n);
1318 delete fNeurons->Remove(n);
1321 // otherwise, its activation becomes the maximum reference
1322 p->Next() = n->Outer();
1323 // at the exit of the while(), the fNext will point
1324 // to the outer node of the neuron starting in p, whose
1325 // activation is the largest.
1327 // the same procedure is made now for all neurons
1328 // for which p is the outer point
1331 TObjArrayIter outerof(p->OuterOf());
1332 while ( (n = (AliITSneuron*)outerof.Next()) ) {
1333 // if the examined neuron has not the largest activation
1334 // it is skipped and removed from array of all neurons
1335 // and of its inner point (its outer is the cursor p)
1336 if (n->Activation() < test) {
1337 p->OuterOf()->Remove(n);
1338 n->Inner()->InnerOf()->Remove(n);
1339 delete fNeurons->Remove(n);
1342 // otherwise, its activation becomes the maximum reference
1343 p->Prev() = n->Inner();
1344 // at the exit of the while(), the fPrev will point
1345 // to the inner node of the neuron ending in p, whose
1346 // activation is the largest.
1348 } // end while (p ...)
1349 } // end for (itheta ...)
1350 } // end for (ilayer ...)
1352 // now the mismatches are solved
1353 Bool_t matchPrev, matchNext;
1354 for (ilayer = 0; ilayer < fNLayers; ilayer++) {
1355 for (itheta = 0; itheta < 180; itheta++) {
1356 // get the array of point in a given layer/theta-slot/sector
1357 lstSector = (TObjArray*)fNodes[ilayer][itheta]->At(sector);
1358 TObjArrayIter nodes(lstSector);
1359 while ( (p = (AliITSnode*)nodes.Next()) ) {
1360 // now p will point to a fPrev and a fNext node.
1361 // Ideally they are placed this way: fPrev --> P --> fNext
1362 // A mismatch happens if the point addressed as fPrev does NOT
1363 // point to p as its fNext. And the same for the point addressed
1365 // In this case, the fNext and fPrev pointers are set to NULL
1366 // and p is excluded from the reconstruction
1367 matchPrev = matchNext= kFALSE;
1368 if (ilayer > 0 && p->Prev() != NULL)
1369 if (p->Prev()->Next() == p) matchPrev = kTRUE;
1370 if (ilayer < 5 && p->Next() != NULL)
1371 if (p->Next()->Prev() == p) matchNext = kTRUE;
1374 else if (ilayer == 5)
1376 if (!matchNext || !matchPrev) {
1377 p->Prev() = p->Next() = 0;
1379 } // end while (p ...)
1380 } // end for (itheta ...)
1381 } // end for (ilayer ...)
1384 //__________________________________________________________________________________
1385 Int_t AliITStrackerANN::SaveTracks(Int_t sector)
1387 /********************************************************************************
1391 Using the fNext and fPrev pointers, the chain is followed
1392 and the track is fitted and saved.
1393 Of course, the track is followed as a chain with a point
1394 for each layer, then the track following starts always
1395 from the clusters in layer 0.
1397 ***********************************************************************************/
1399 // if not initialized, the tracks TobjArray is initialized
1400 if (!fFoundTracks) fFoundTracks = new TObjArray;
1402 // meaningful counters
1403 Int_t itheta, ilayer, l;
1404 TObjArray *lstSector = 0;
1405 AliITSnode *p = 0, *q = 0, **node = new AliITSnode*[fNLayers];
1406 for (l = 0; l < fNLayers; l++) node[l] = 0;
1409 array = new AliITSnode*[fNLayers + 1];
1410 for (l = 0; l <= fNLayers; l++) array[l] = 0;
1411 array[0] = new AliITSnode();
1412 array[0]->X() = fVertexX;
1413 array[0]->Y() = fVertexY;
1414 array[0]->Z() = fVertexZ;
1415 array[0]->ErrX2() = fVertexErrorX;
1416 array[0]->ErrY2() = fVertexErrorY;
1417 array[0]->ErrZ2() = fVertexErrorZ;
1419 Double_t *param = new Double_t[8];
1421 // scan the whole collection of nodes
1422 for (ilayer = 0; ilayer < 1; ilayer++) {
1423 for (itheta = 0; itheta < 180; itheta++) {
1424 // get the array of point in a given layer/theta-slot/sector
1425 lstSector = (TObjArray*)fNodes[ilayer][itheta]->At(sector);
1426 TObjArrayIter nodes(lstSector);
1427 while ( (p = (AliITSnode*)nodes.Next()) ) {
1428 for (q = p; q; q = q->Next()) {
1432 //if (!RiemannFit(fNLayers, node, param)) continue;
1433 // initialization of Kalman Filter Tracking
1434 AliITSclusterV2 *cluster = (AliITSclusterV2*)GetCluster(node[0]->ClusterRef());
1435 Int_t mod = cluster->GetDetectorIndex();
1436 Int_t lay, lad, det;
1437 fGeom->GetModuleId(mod, lay, lad, det);
1438 Float_t y0 = cluster->GetY();
1439 Float_t z0 = cluster->GetZ();
1440 AliITStrackSA* trac = new AliITStrackSA(lay, lad, det,
1442 param[4], param[7], param[3], 1);
1443 for (l = 0; l < fNLayers; l++) {
1444 cluster = (AliITSclusterV2*)GetCluster(node[l]->ClusterRef());
1445 if (cluster) trac->AddClusterV2(l, (node[l]->ClusterRef() & 0x0fffffff)>>0);
1447 AliITStrackV2* ot = new AliITStrackV2(*trac);
1448 ot->ResetCovariance();
1449 ot->ResetClusters();
1450 if (RefitAt(49.,ot,trac)) { //fit from layer 1 to layer 6
1451 AliITStrackV2 *otrack2 = new AliITStrackV2(*ot);
1452 otrack2->ResetCovariance();
1453 otrack2->ResetClusters();
1454 //fit from layer 6 to layer 1
1455 if (RefitAt(3.7,otrack2,ot)) fFoundTracks->AddLast(otrack2);
1457 // end of Kalman Filter fit
1465 //__________________________________________________________________________________
1466 void AliITStrackerANN::ExportTracks(const char *filename) const
1468 // Exports found tracks into a TTree of AliITStrackV2 objects
1469 TFile *file = new TFile(filename, "RECREATE");
1470 TTree *tree = new TTree("TreeT-ANN", "Tracks found in ITS stand-alone with Neural Tracking");
1471 AliITStrackV2 *track = 0;
1472 tree->Branch("Tracks", &track, "AliITStrackV2");
1473 TObjArrayIter tracks(fFoundTracks);
1474 while ( (track = (AliITStrackV2*)tracks.Next()) ) {
1483 //__________________________________________________________________________________
1484 void AliITStrackerANN::CleanNetwork()
1486 // Removes deactivated units from the network
1488 AliITSneuron *unit = 0;
1489 TObjArrayIter neurons(fNeurons);
1490 while ( (unit = (AliITSneuron*)neurons.Next()) ) {
1491 if (unit->Activation() < fActMinimum) {
1492 unit->Inner()->InnerOf()->Remove(unit);
1493 unit->Outer()->OuterOf()->Remove(unit);
1494 delete fNeurons->Remove(unit);
1500 AliITSneuron *enemy = 0;
1502 while ( (unit = (AliITSneuron*)neurons.Next()) ) {
1503 nIn = (Int_t)unit->Inner()->InnerOf()->GetSize();
1504 nOut = (Int_t)unit->Outer()->OuterOf()->GetSize();
1505 if (nIn < 2 && nOut < 2) continue;
1508 TObjArrayIter competing(unit->Inner()->InnerOf());
1509 while ( (enemy = (AliITSneuron*)competing.Next()) ) {
1510 if (unit->Activation() > enemy->Activation()) {
1511 enemy->Inner()->InnerOf()->Remove(enemy);
1512 enemy->Outer()->OuterOf()->Remove(enemy);
1513 delete fNeurons->Remove(enemy);
1516 unit->Inner()->InnerOf()->Remove(unit);
1517 unit->Outer()->OuterOf()->Remove(unit);
1518 delete fNeurons->Remove(unit);
1523 if (removed) continue;
1526 TObjArrayIter competing(unit->Outer()->OuterOf());
1527 while ( (enemy = (AliITSneuron*)competing.Next()) ) {
1528 if (unit->Activation() > enemy->Activation()) {
1529 enemy->Inner()->InnerOf()->Remove(enemy);
1530 enemy->Outer()->OuterOf()->Remove(enemy);
1531 delete fNeurons->Remove(enemy);
1534 unit->Inner()->InnerOf()->Remove(unit);
1535 unit->Outer()->OuterOf()->Remove(unit);
1536 delete fNeurons->Remove(unit);
1545 //__________________________________________________________________________________
1546 Int_t AliITStrackerANN::StoreTracks()
1548 // Stores the tracks found in a single neural tracking step.
1549 // In order to do this, it sects each neuron which has a point
1550 // in the first layer.
1553 // if not initialized, the tracks TobjArray is initialized
1554 if (!fFoundTracks) fFoundTracks = new TObjArray;
1556 Int_t i, check, stored = 0;
1557 Double_t testAct = 0;
1558 AliITSneuron *unit = 0, *cursor = 0, *fwd = 0;
1559 AliITSnode *node = 0;
1560 TObjArrayIter iter(fNeurons), *fwdIter;
1561 TObjArray *removedUnits = new TObjArray(0);
1562 removedUnits->SetOwner(kFALSE);
1563 AliITStrackANN annTrack(fNLayers);
1566 unit = (AliITSneuron*)iter.Next();
1568 if (unit->Inner()->GetLayer() > 0) continue;
1569 annTrack.SetNode(unit->Inner()->GetLayer(), unit->Inner());
1570 annTrack.SetNode(unit->Outer()->GetLayer(), unit->Outer());
1571 node = unit->Outer();
1572 removedUnits->AddLast(unit);
1574 testAct = fActMinimum;
1575 fwdIter = (TObjArrayIter*)node->InnerOf()->MakeIterator();
1578 cursor = (AliITSneuron*)fwdIter->Next();
1580 if (cursor->Used()) continue;
1581 if (cursor->Activation() >= testAct) {
1582 testAct = cursor->Activation();
1587 removedUnits->AddLast(fwd);
1588 node = fwd->Outer();
1589 annTrack.SetNode(node->GetLayer(), node);
1591 check = annTrack.CheckOccupation();
1595 //if (!RiemannFit(fNLayers, trackitem, param)) continue;
1596 if (!annTrack.RiemannFit()) continue;
1597 // initialization of Kalman Filter Tracking
1598 AliITSclusterV2 *cluster = (AliITSclusterV2*)GetCluster(annTrack[0]->ClusterRef());
1599 Int_t mod = cluster->GetDetectorIndex();
1600 Int_t lay, lad, det;
1601 fGeom->GetModuleId(mod, lay, lad, det);
1602 Float_t y0 = cluster->GetY();
1603 Float_t z0 = cluster->GetZ();
1604 AliITStrackSA* trac = new AliITStrackSA(lay, lad, det, y0, z0,
1605 annTrack.Phi(), annTrack.TanLambda(),
1606 annTrack.Curv(), 1);
1607 for (Int_t l = 0; l < fNLayers; l++) {
1608 if (!annTrack[l]) continue;
1609 cluster = (AliITSclusterV2*)GetCluster(annTrack[l]->ClusterRef());
1610 if (cluster) trac->AddClusterV2(l, (annTrack[l]->ClusterRef() & 0x0fffffff)>>0);
1612 AliITStrackV2* ot = new AliITStrackV2(*trac);
1613 ot->ResetCovariance();
1614 ot->ResetClusters();
1615 if (RefitAt(49.,ot,trac)) { //fit from layer 1 to layer 6
1616 AliITStrackV2 *otrack2 = new AliITStrackV2(*ot);
1617 otrack2->ResetCovariance();
1618 otrack2->ResetClusters();
1619 //fit from layer 6 to layer 1
1620 if (RefitAt(3.7,otrack2,ot)) fFoundTracks->AddLast(otrack2);
1622 // end of Kalman Filter fit
1624 for (i = 0; i < fNLayers; i++) {
1625 //node = (AliITSnode*)removedPoints->At(i);
1629 fwdIter = (TObjArrayIter*)removedUnits->MakeIterator();
1631 cursor = (AliITSneuron*)fwdIter->Next();
1641 Double_t AliITStrackerANN::Weight(AliITSneuron *nAB, AliITSneuron *nBC)
1643 /***********************************************************************************
1644 * Calculation of neural weight.
1645 * The implementation of positive neural weight is set only in the case
1646 * of connected units (e.g.: A->B with B->C).
1647 * Given that B is the **common** point. care should be taken to pass
1648 * as FIRST argument the neuron going "to" B, and
1649 * as SECOND argument the neuron starting "from" B
1650 * anyway, a check is put in order to return 0.0 when arguments are not well given.
1651 ***********************************************************************************/
1653 if (nAB->Outer() != nBC->Inner()) {
1654 if (nBC->Outer() == nAB->Inner()) {
1655 AliITSneuron *temp = nAB;
1659 if (fMsgLevel >= 3) {
1660 Info("Weight", "Switching wrongly ordered arguments.");
1663 Warning("Weight", "Not connected segments. Returning 0.0");
1667 AliITSnode *pA = nAB->Inner();
1668 AliITSnode *pB = nAB->Outer();
1669 AliITSnode *pC = nBC->Outer();
1671 TVector3 vAB(pB->X() - pA->X(), pB->Y() - pA->Y(), pB->Z() - pA->Z());
1672 TVector3 vBC(pC->X() - pB->X(), pC->Y() - pB->Y(), pC->Z() - pB->Z());
1674 Double_t weight = 1.0 - sin(vAB.Angle(vBC));
1675 return fGain2CostRatio * TMath::Power(weight, fExponent);
1680 /******************************************
1681 ******************************************
1682 *** AliITStrackerANN::AliITSnode class ***
1683 ******************************************
1684 ******************************************/
1686 //__________________________________________________________________________________
1687 inline AliITStrackerANN::AliITSnode::AliITSnode()
1688 : fUsed(kFALSE), fClusterRef(-1),
1689 fMatches(NULL), fInnerOf(NULL), fOuterOf(NULL),
1690 fNext(NULL), fPrev(NULL)
1692 // Constructor for the embedded 'AliITSnode' class.
1693 // It initializes all pointer-like objects.
1696 fEX2 = fEY2 = fEZ2 = 0.0;
1699 //__________________________________________________________________________________
1700 AliITStrackerANN::AliITSnode::~AliITSnode()
1702 // Destructor for the embedded 'AliITSnode' class.
1703 // It should clear the object arrays, but it is possible
1704 // that some objects still are useful after the point deletion
1705 // then the arrays are cleared but their objects are owed by
1706 // another TCollection object, and not deleted.
1707 // For safety reasons, all the pointers are set to zero.
1709 fInnerOf->SetOwner(kFALSE);
1713 fOuterOf->SetOwner(kFALSE);
1717 fMatches->SetOwner(kFALSE);
1725 //__________________________________________________________________________________
1726 inline Double_t AliITStrackerANN::AliITSnode::GetPhi() const
1728 // Calculates the 'phi' (azimutal) angle, and returns it
1729 // in the range between 0 and 2Pi radians.
1732 q = TMath::ATan2(fY,fX);
1736 return q + TMath::TwoPi();
1739 //__________________________________________________________________________________
1740 inline Double_t AliITStrackerANN::AliITSnode::GetError(Option_t *option)
1742 // Returns the error or the square error of
1743 // values related to the coordinates in different systems.
1744 // The option argument specifies the coordinate error desired:
1746 // "R2" --> error in transverse radius
1747 // "R3" --> error in spherical radius
1748 // "PHI" --> error in azimuthal angle
1749 // "THETA" --> error in polar angle
1750 // "SQ" --> get the square of error
1752 // In order to get the error on the cartesian coordinates
1753 // reference to the inline ErrX2(), ErrY2() adn ErrZ2() methods.
1755 TString opt(option);
1756 Double_t errorSq = 0.0;
1759 if (opt.Contains("R2")) {
1760 errorSq = fX*fX*fEX2 + fY*fY*fEY2;
1761 errorSq /= GetR2sq();
1763 else if (opt.Contains("R3")) {
1764 errorSq = fX*fX*fEX2 + fY*fY*fEY2 + fZ*fZ*fEZ2;
1765 errorSq /= GetR3sq();
1767 else if (opt.Contains("PHI")) {
1768 errorSq = fY*fY*fEX2;
1769 errorSq += fX*fX*fEY2;
1770 errorSq /= GetR2sq() * GetR2sq();
1772 else if (opt.Contains("THETA")) {
1773 errorSq = fZ*fZ * (fX*fX*fEX2 + fY*fY*fEY2);
1774 errorSq += GetR2sq() * GetR2sq() * fEZ2;
1775 errorSq /= GetR3sq() * GetR3sq() * GetR2() * GetR2();
1778 if (!opt.Contains("SQ"))
1779 return TMath::Sqrt(errorSq);
1786 /********************************************
1787 ********************************************
1788 *** AliITStrackerANN::AliITSneuron class ***
1789 ********************************************
1790 ********************************************/
1792 //__________________________________________________________________________________
1793 AliITStrackerANN::AliITSneuron::AliITSneuron
1794 (AliITSnode *inner, AliITSnode *outer, Double_t minAct, Double_t maxAct)
1795 : fUsed(0), fInner(inner), fOuter(outer)
1797 // Default neuron constructor
1798 fActivation = gRandom->Rndm() * (maxAct-minAct) + minAct;
1799 fGain = new TObjArray;
1802 //__________________________________________________________________________________
1803 inline Double_t AliITStrackerANN::AliITSneuron::Activate(Double_t temperature)
1805 // This computes the new activation of a neuron, and returns
1806 // its activation variation as a consequence of the updating.
1809 // - the 'temperature' parameter for the neural activation logistic function
1812 // - collects the total gain, by summing the products
1813 // of the activation of each sequenced unit by the relative weight.
1814 // - collects the total cost, by summing the activations of
1815 // all competing units
1816 // - passes the sum of gain - cost to the activation function and
1817 // calculates the new activation
1820 // - the difference between the old activation and the new one
1824 Double_t sumGain = 0.0; // total contribution from chained neurons
1825 Double_t sumCost = 0.0; // total contribution from crossing neurons
1826 Double_t input; // total input
1827 Double_t actOld, actNew; // old and new values for the activation
1828 AliITSneuron *linked = 0; // cursor for scanning the neuron arrays (for link check)
1829 AliITSlink *link; // cursor for scanning the synapses arrays (for link check)
1830 TObjArrayIter *iterator = 0; // pointer to the iterator along the neuron arrays
1832 // sum contributions from the correlated units
1833 iterator = (TObjArrayIter*)fGain->MakeIterator();
1835 link = (AliITSlink*)iterator->Next();
1837 sumGain += link->Contribution();
1841 // sum contributions from the competing units:
1842 // the ones which have the same starting point...
1843 iterator = (TObjArrayIter*)fInner->InnerOf()->MakeIterator();
1845 linked = (AliITSneuron*)iterator->Next();
1847 if (linked == this) continue;
1848 sumCost += linked->fActivation;
1851 // ...and the ones which have the same ending point
1852 iterator = (TObjArrayIter*)fOuter->OuterOf()->MakeIterator();
1854 linked = (AliITSneuron*)iterator->Next();
1856 if (linked == this) continue;
1857 sumCost += linked->fActivation;
1860 // calculate the total input as the difference between gain and cost
1861 input = (sumGain - sumCost) / temperature;
1862 actOld = fActivation;
1863 // calculate the final output
1864 #ifdef NEURAL_LINEAR
1865 if (input <= -2.0 * temperature)
1867 else if (input >= 2.0 * temperature)
1870 actNew = input / (4.0 * temperature) + 0.5;
1872 actNew = 1.0 / (1.0 + TMath::Exp(-input));
1874 fActivation = actNew;
1876 // return the activation variation
1877 return TMath::Abs(actNew - actOld);
1882 /******************************************
1883 ******************************************
1884 *** AliITStrackerANN::AliITSlink class ***
1885 ******************************************
1886 ******************************************/
1888 // No methods defined non-inline
1892 /**********************************************
1893 **********************************************
1894 *** AliITStrackerANN::AliITStrackANN class ***
1895 **********************************************
1896 **********************************************/
1898 //__________________________________________________________________________________
1899 AliITStrackerANN::AliITStrackANN::AliITStrackANN(Int_t dim) : fNPoints(dim)
1901 // Default constructor for the AliITStrackANN class
1916 fNode = new AliITSnode*[dim];
1917 for (i = 0; i < dim; i++) fNode[i] = 0;
1921 //__________________________________________________________________________________
1922 Int_t AliITStrackerANN::AliITStrackANN::CheckOccupation() const
1924 // Returns the number of pointers fNode which are not NULL
1927 Int_t count = 0; // counter for how many points are stored in the track
1929 for (i = 0; i < fNPoints; i++) {
1930 if (fNode[i] != NULL) count++;
1936 //__________________________________________________________________________________
1937 Bool_t AliITStrackerANN::AliITStrackANN::RiemannFit()
1939 // Performs the Riemann Sphere fit for the given points to a circle
1940 // and then uses the fit parameters to fit a helix in space.
1943 // - kTRUE if all operations have been performed
1944 // - kFALSE if the numbers risk to lead to an arithmetic violation
1946 Int_t i, j, count, dim = fNPoints;
1948 // First check for all points
1949 count = CheckOccupation();
1950 if (count != fNPoints) {
1951 Error ("AliITStrackANN::RiemannFit", "CheckOccupations returns %d, fNPoints = %d ==> MISMATCH", count, fNPoints);
1957 for (i = 0; i < dim; i++) m1(i,0) = 1.0;
1959 // matrix of Rieman projection coordinates
1960 TMatrixD coords(dim,3);
1961 for (i = 0; i < dim; i++) {
1962 coords(i,0) = fNode[i]->X();
1963 coords(i,1) = fNode[i]->Y();
1964 coords(i,2) = fNode[i]->GetR2sq();
1967 // matrix of weights
1968 Double_t xterm, yterm, ex, ey;
1969 TMatrixD weights(dim,dim);
1970 for (i = 0; i < dim; i++) {
1971 xterm = fNode[i]->X() * fNode[i]->GetPhi() - fNode[i]->Y() / fNode[i]->GetR2();
1972 ex = fNode[i]->ErrX2();
1973 yterm = fNode[i]->Y() * fNode[i]->GetPhi() + fNode[i]->X() / fNode[i]->GetR2();
1974 ey = fNode[i]->ErrY2();
1975 weights(i,i) = fNode[i]->GetR2sq() / (xterm * xterm * ex + yterm * yterm * ey );
1978 // weighted sample mean
1979 Double_t meanX = 0.0, meanY = 0.0, meanW = 0.0, sw = 0.0;
1980 for (i = 0; i < dim; i++) {
1981 meanX += weights(i,i) * coords(i,0);
1982 meanY += weights(i,i) * coords(i,1);
1983 meanW += weights(i,i) * coords(i,2);
1990 // sample covariance matrix
1991 for (i = 0; i < dim; i++) {
1992 coords(i,0) -= meanX;
1993 coords(i,1) -= meanY;
1994 coords(i,2) -= meanW;
1996 TMatrixD coordsT(TMatrixD::kTransposed, coords);
1997 TMatrixD weights4coords(weights, TMatrixD::kMult, coords);
1998 TMatrixD sampleCov(coordsT, TMatrixD::kMult, weights4coords);
1999 for (i = 0; i < 3; i++) {
2000 for (j = i + 1; j < 3; j++) {
2001 sampleCov(i,j) = sampleCov(j,i) = (sampleCov(i,j) + sampleCov(j,i)) * 0.5;
2005 // Eigenvalue problem solving for V matrix
2007 TVectorD eval(3), n(3);
2008 TMatrixD evec = sampleCov.EigenVectors(eval);
2009 if (eval(1) < eval(ileast)) ileast = 1;
2010 if (eval(2) < eval(ileast)) ileast = 2;
2011 n(0) = evec(0, ileast);
2012 n(1) = evec(1, ileast);
2013 n(2) = evec(2, ileast);
2015 // c - known term in the plane intersection with Riemann axes
2016 Double_t c = -(meanX * n(0) + meanY * n(1) + meanW * n(2));
2018 // center and radius of fitted circle
2019 Double_t xc, yc, radius, curv;
2020 xc = -n(0) / (2. * n(2));
2021 yc = -n(1) / (2. * n(2));
2022 radius = (1. - n(2)*n(2) - 4.*c*n(2)) / (4. * n(2) * n(2));
2024 if (radius <= 0.E0) {
2025 Error("RiemannFit", "Radius = %f less than zero!!!", radius);
2028 radius = TMath::Sqrt(radius);
2029 curv = 1.0 / radius;
2031 // evaluating signs for curvature and others
2032 Double_t phi1 = 0.0, phi2, temp1, temp2, phi0, sumdphi = 0.0;
2033 AliITSnode *p = fNode[0];
2035 for (i = 1; i < dim; i++) {
2036 p = (AliITSnode*)fNode[i];
2041 if (temp1 > fgkPi && temp2 < fgkPi)
2043 else if (temp1 < fgkPi && temp2 > fgkPi)
2045 sumdphi += temp2 - temp1;
2048 if (sumdphi < 0.E0) curv = -curv;
2049 Double_t diff, angle = TMath::ATan2(yc, xc);
2051 phi0 = angle + 0.5 * TMath::Pi();
2053 phi0 = angle - 0.5 * TMath::Pi();
2054 diff = angle - phi0;
2056 Double_t dt, temp = TMath::Sqrt(xc*xc + yc*yc) - radius;
2061 //cout << "Dt = " << dt << endl;
2063 Double_t halfC = 0.5 * curv, test;
2064 Double_t *s = new Double_t[dim], *zz = new Double_t[dim], *ws = new Double_t[dim];
2065 for (j = 0; j < 6; j++) {
2069 s[j] = (p->GetR2sq() - dt * dt) / (1. + curv * dt);
2071 if (fabs(s[j]) < 1.E-6) s[j] = 0.;
2073 Error("RiemannFit", "Square root argument error: %17.15g < 0", s[j]);
2077 s[j] = TMath::Sqrt(s[j]);
2078 //cout << "Curv = " << halfC << " --- s[" << j << "] = " << s[j] << endl;
2080 test = TMath::Abs(s[j]);
2083 s[j] = ((s[j] > 0.) ? 0.99999999999 : -0.9999999999);
2085 Error("RiemannFit", "Value too large: %17.15g", s[j]);
2091 s[j] = TMath::ASin(s[j]) / halfC;
2092 ws[j] = 1.0 / (p->ErrZ2());
2095 // second tep final fit
2096 Double_t s2Sum = 0.0, zSum = 0.0, szSum = 0.0, sSum = 0.0, sumw = 0.0;
2097 for (i = 0; i < dim; i++) {
2098 s2Sum += ws[i] * s[i] * s[i];
2099 zSum += ws[i] * zz[i];
2100 sSum += ws[i] * s[i];
2101 szSum += ws[i] * s[i] * zz[i];
2108 temp = s2Sum - sSum*sSum;
2111 dz = (s2Sum*zSum - sSum*szSum) / temp;
2112 tanL = (szSum - sSum*zSum) / temp;