--- /dev/null
+// ---------------------------------------------------------------------------------
+// AliITStrackerANN
+// ---------------------------------------------------------------------------------
+// Neural Network-based algorithm for track reconstruction in the ALICE ITS.
+// The class contains all that is necessary in order to perform the tracking
+// using the ITS clusters in the V2 format.
+// The class is organized with some embedded objects which serve for
+// data management, and all setters and getters for all parameters.
+// The usage of the class from a tracking macro should be based essentially
+// in the initialization (constructor) and the call of a method which
+// performs in the right sequence all the operations.
+// Temporarily (and maybe definitively) all the intermediate steps are
+// public functions which could be called independently, but they do not
+// produce any result if all the preventively required steps have not been
+// performed.
+// ---------------------------------------------------------------------------------
+// Author: Alberto Pulvirenti (university of Catania)
+// Email : alberto.pulvirenti@ct.infn.it
+// ---------------------------------------------------------------------------------
+
+#include <TMath.h>
+#include <TString.h>
+#include <TObjArray.h>
+#include <TVector3.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TRandom.h>
+#include <TMatrixD.h>
+
+#include "AliITSgeom.h"
+#include "AliITStrackSA.h"
+#include "AliITSclusterV2.h"
+#include "AliITStrackV2.h"
+
+#include "AliITStrackerANN.h"
+
+class iostream;
+using namespace std;
+
+ClassImp(AliITStrackerANN)
+
+//__________________________________________________________________________________
+AliITStrackerANN::AliITStrackerANN(const AliITSgeom *geom, Int_t msglev)
+: AliITStrackerV2(geom), fMsgLevel(msglev)
+{
+/**************************************************************************
+
+ CONSTRUCTOR (standard-to-use version)
+
+ Arguments:
+ 1) the ITS geometry used in the generated event
+ 2) the flag for log-messages writing
+
+ The AliITSgeometry is used along the class,
+ in order to translate the local AliITSclusterV2 coordinates
+ into the Global reference frame, which is necessary for the
+ Neural Network algorithm to operate.
+ In case of serialized use, the log messages should be excluded,
+ in order to save real execution time.
+
+ Operations:
+ - all pointer data members are initialized
+ (if possible, otherwise are set to NULL)
+ - all numeric data members are initialized to
+ values which allow the Neural Network to operate
+ even if nothing is externally set.
+
+ NOTE: it is possible that tracking an event
+ with these default values results in a non-sense.
+
+ **************************************************************************/
+
+ Int_t i;
+
+ // Get ITS geometry
+ fGeom = (AliITSgeom*)geom;
+
+ // Initialize the array of first module indexes per layer
+ fNLayers = geom->GetNlayers();
+ fFirstModInLayer = new Int_t[fNLayers + 1];
+ for (i = 0; i < fNLayers; i++) {
+ fFirstModInLayer[i] = fGeom->GetModuleIndex(i + 1, 1, 1);
+ }
+ fFirstModInLayer[fNLayers] = geom->GetIndexMax();
+
+ // initialization: no curvature cut steps
+ fCurvNum = 0;
+ fCurvCut = 0;
+
+ // initialization: 4 sectors (one for each quadrant)
+ fSectorNum = 4;
+ fSectorWidth = fgkHalfPi;
+
+ // initialization: theta offset of 20 degrees
+ fPolarInterval = 20.0;
+
+ // initialization: array structure not defined
+ fStructureOK = kFALSE;
+
+ // initialization: vertex in the origin
+ fVertexX = 0.0;
+ fVertexY = 0.0;
+ fVertexZ = 0.0;
+
+ // initialization: uninitialized point array
+ fNodes = 0;
+
+ // initialization: very large (unuseful) cut values
+ Int_t ilayer;
+ for (ilayer = 0; ilayer < 6; ilayer++) {
+ fThetaCut2D[ilayer] = TMath::Pi();
+ fThetaCut3D[ilayer] = TMath::Pi();
+ fHelixMatchCut[ilayer] = 1.0;
+ }
+
+ // initialization: inictial activation range between 0.3 and 0.7
+ fEdge1 = 0.3;
+ fEdge2 = 0.7;
+
+ // initialization: neural network operation & weights
+ fTemperature = 1.0;
+ fStabThreshold = 0.001;
+ fGain2CostRatio = 1.0;
+ fExponent = 1.0;
+ fActMinimum = 0.5;
+
+ // initialization: uninitialized array of neurons
+ fNeurons = 0;
+
+ // initialization: uninitialized array of tracks
+ fFoundTracks = 0;
+}
+
+//__________________________________________________________________________________
+void AliITStrackerANN::SetCuts
+(Int_t ncurv, Double_t *curv, Double_t *theta2D, Double_t *theta3D, Double_t *helix)
+
+/**************************************************************************
+
+ CUT SETTER
+
+ Allows for the definition of all kind of geometric cuts
+ which have been studied in for the creation of a neuron
+ from a pair of clusters C1 and C2 on consecutive layers.
+ Neuron will be created only if the pair passes ALL of these cuts.
+
+ At the moment, we define 4 kinds of geometrical cuts:
+ a) cut on the difference of the polar 'theta' angle;
+ b) cut on the angle between origin->C1 and C1->C2 in space;
+ c) cut on the curvature of the circle passing
+ through C1, C2 and the primary vertex;
+ d) cut on heli-matching of the same three points.
+
+ Arguments:
+ 1) the number of curvature cut steps
+ 2) the array of curvature cuts for each step
+ (its dimension is given by the argument 1)
+ 3) array of 5 cut values (one for each consecutive lauer pair)
+ related to cut a)
+ 4) array of 5 cut values (one for each consecutive lauer pair)
+ related to cut b)
+ 5) array of 5 cut values (one for each consecutive lauer pair)
+ related to cut c)
+
+ Operations:
+ - gets the values for each cut and stores them in data members
+ - in the case of curvature cuts, the cut array
+ (whose size is not fixed) is allocated
+
+ NOTE: in case that the user wants to set onyl ONE of the 4 cuts array,
+ he can simply pass NULL arguments for other cuts, and (eventually)
+ a ZERO as the first argument (if curvature cuts have not to be set).
+
+ Anyway, all the cuts have to be set at least once.
+
+ **************************************************************************/
+{
+ // counter
+ Int_t i;
+
+ /*** Curvature cut setting ***/
+
+ // first of all, the curvature cuts are sorted in increasing order
+ // (from the smallest to the largest one)
+ Int_t *ind = new Int_t[ncurv];
+ TMath::Sort(ncurv, curv, ind, kFALSE);
+ // then, the curvature cut array is allocated and filled
+ // (a message with the list of defined cuts can be optionally shown)
+ fCurvCut = new Double_t[ncurv];
+ if (fMsgLevel >= 1) cout << "Number of curvature cuts: " << ncurv << endl;
+ for (i = 0; i < ncurv; i++) {
+ fCurvCut[i] = curv[ind[i]];
+ if (fMsgLevel >= 1) cout << " - " << fCurvCut[i] << endl;
+ }
+ fCurvNum = ncurv;
+
+ /*** 'Fixed' cuts setting ***/
+
+ // checks what cuts have to be set
+ Bool_t doTheta2D = (theta2D != 0);
+ Bool_t doTheta3D = (theta3D != 0);
+ Bool_t doHelix = (helix != 0);
+ // sets the cuts for all layer pairs
+ for (i = 0; i < fNLayers - 1; i++) {
+ if (doTheta2D) fThetaCut2D[i] = theta2D[i];
+ if (doTheta3D) fThetaCut3D[i] = theta3D[i];
+ if (doHelix) fHelixMatchCut[i] = helix[i];
+ }
+ // if required, lists the cuts
+ if (!fMsgLevel < 2) return;
+ cout << "Theta 2D cuts: ";
+ if (!doTheta2D) {
+ cout << "<not set>" << endl;
+ }
+ else {
+ cout << endl;
+ for (i = 0; i < fNLayers - 1; i++) {
+ cout << "For layers " << i << " --> " << i + 1;
+ cout << " cut = " << fThetaCut2D[i] << endl;
+ }
+ }
+ cout << "---" << endl;
+ cout << "Theta 3D cuts: ";
+ if (!doTheta3D) {
+ cout << "<not set>" << endl;
+ }
+ else {
+ cout << endl;
+ for (i = 0; i < fNLayers - 1; i++) {
+ cout << "For layers " << i << " --> " << i + 1;
+ cout << " cut = " << fThetaCut3D[i] << endl;
+ }
+ }
+ cout << "---" << endl;
+ cout << "Helix-match cuts: ";
+ if (!doHelix) {
+ cout << "<not set>" << endl;
+ }
+ else {
+ cout << endl;
+ for (i = 0; i < fNLayers - 1; i++) {
+ cout << "For layers " << i << " --> " << i + 1;
+ cout << " cut = " << fHelixMatchCut[i] << endl;
+ }
+ }
+ cout << "---" << endl;
+}
+
+//__________________________________________________________________________________
+Bool_t AliITStrackerANN::GetGlobalXYZ
+(Int_t refIndex,
+ Double_t &x, Double_t &y, Double_t &z, Double_t &ex2, Double_t &ey2, Double_t &ez2)
+
+/**************************************************************************
+
+ LOCAL TO GLOBAL TRANSLATOR
+
+ Taking information from the ITS geometry stored in the class,
+ gets a stored AliITScluster and calculates it coordinates
+ and errors in the global reference frame.
+ These values are stored in the variables,
+ which are passed by reference.
+
+ Arguments:
+ 1) reference index for the cluster to use
+ 2) (by reference) place to store the global X-coord into
+ 3) (by reference) place to store the global Y-coord into
+ 4) (by reference) place to store the global Z-coord into
+ 5) (by reference) place to store the global X-coord error into
+ 6) (by reference) place to store the global Y-coord error into
+ 7) (by reference) place to store the global Z-coord error into
+
+ Operations:
+ essentially, determines the ITS module index from the
+ detector index of the AliITSclusterV2 object, and extracts
+ the roto-translation from the ITS geometry, to convert
+ the local module coordinates into the global ones.
+
+ Return value:
+ - kFALSE if the given cluster index points to a non-existing
+ cluster, or if the layer number makes no sense (< 0 or > 6).
+ - otherwise, kTRUE (meaning a successful operation).
+
+ **************************************************************************/
+{
+ // checks if the layer is correct
+ Int_t ilayer = (refIndex & 0xf0000000) >> 28;
+ if (ilayer < 0 || ilayer >= fNLayers) {
+ Error("GetGlobalXYZ", "Wrong layer number: %d [range: %d - %d]", ilayer, 0, fNLayers);
+ return kFALSE;
+ }
+ // checks if the referenced cluster exists and corresponds to the passed reference
+ AliITSclusterV2 *refCluster = (AliITSclusterV2*) GetCluster(refIndex);
+ if (!refCluster) {
+ Error("GetGlobalXYZ", "Cluster not found for index %d", refIndex);
+ return kFALSE;
+ }
+
+ // determine the detector number
+ Int_t detID = refCluster->GetDetectorIndex() + fFirstModInLayer[ilayer];
+
+ // get rotation matrix
+ Double_t rot[9];
+ fGeom->GetRotMatrix(detID, rot);
+
+ // get translation vector
+ Float_t tx,ty,tz;
+ fGeom->GetTrans(detID, tx, ty, tz);
+
+ // determine r and phi for the reference conversion
+ Double_t r = -(Double_t)tx * rot[1] + (Double_t)ty * rot[0];
+ if (ilayer == 0) r = -r;
+ Double_t phi = TMath::ATan2(rot[1],rot[0]);
+ if (ilayer == 0) phi -= fgkPi;
+
+ // sets values for X, Y, Z in global coordinates and their errors
+ Double_t cosPhi = TMath::Cos(phi);
+ Double_t sinPhi = TMath::Sin(phi);
+ x = r*cosPhi + refCluster->GetY()*sinPhi;
+ y = -r*sinPhi + refCluster->GetY()*cosPhi;
+ z = refCluster->GetZ();
+ ex2 = refCluster->GetSigmaY2()*sinPhi*sinPhi;
+ ey2 = refCluster->GetSigmaY2()*cosPhi*cosPhi;
+ ez2 = refCluster->GetSigmaZ2();
+
+ return kTRUE;
+}
+
+//__________________________________________________________________________________
+AliITStrackerANN::AliITSnode* AliITStrackerANN::AddNode(Int_t refIndex)
+
+/**************************************************************************
+
+ GENERATOR OF NEURAL NODES
+
+ Fills the array of neural 'nodes', which are the ITS clusters
+ translated in the global reference frame.
+ Given that the global coordinates are used many times, they are
+ stored in a well-defined structure, in the form of an embedded class.
+ Moreover, this class allows a faster navigation among points
+ and neurons, by means of some object arrays, storing only the
+ neurons which start from, or end to, the given node.
+ Finally, each node contains all the other nodes which match it
+ with respect to the fixed walues, in order to perform a faster
+ neuron-creation phase.
+
+ Arguments:
+ 1) reference index of the correlated AliITSclusterV2 object
+
+ Operations:
+ - allocates the new AliITSnode objects
+ - initializes its object arrays
+ - from the global coordinates, calculates the
+ 'phi' and 'theta' coordinates, in order to store it
+ into the correct theta-slot and azimutal sector.
+
+ REturn values:
+ - the pointer of the creater AliITSnode object
+ - in case of errors, a waring is given and a NULL is returned
+
+ **************************************************************************/
+{
+ // create object and set the reference
+ AliITSnode *node = new AliITSnode;
+ if (!node) {
+ Warning("AddNode", "Error occurred when allocating AliITSnode");
+ return 0;
+ }
+ node->ClusterRef() = refIndex;
+
+ // calls the conversion function, which makes also some checks
+ // (layer number within range, existence of referenced cluster)
+ if ( !GetGlobalXYZ (
+ refIndex,
+ node->X(), node->Y(), node->Z(),
+ node->ErrX2(), node->ErrY2(), node->ErrZ2()
+ ) ) {return 0;}
+
+ // initializes the object arrays
+ node->Matches() = new TObjArray;
+ node->InnerOf() = new TObjArray;
+ node->OuterOf() = new TObjArray;
+
+ // finds azimutal and polar sector (in degrees)
+ Double_t phi = node->GetPhi() * 180.0 / fgkPi;
+ Double_t theta = node->GetTheta() * 180.0 / fgkPi;
+ Int_t isector = (Int_t)(phi / fSectorWidth);
+ Int_t itheta = (Int_t)theta;
+ Int_t ilayer = (refIndex & 0xf0000000) >> 28;
+
+ // selects the right TObjArray to store object into
+ TObjArray *sector = (TObjArray*)fNodes[ilayer][itheta]->At(isector);
+ sector->AddLast(node);
+
+ return node;
+}
+
+//__________________________________________________________________________________
+void AliITStrackerANN::CreateArrayStructure(Int_t nsectors)
+
+/**************************************************************************
+
+ ARRAY STRUCTURE CREATOR
+
+ Creates a structure of nested TObjArray's where the AliITSnode's
+ have to be stored:
+ - the first level is made by 6 arrays (one for each layer)
+ - the second level is made by 180 arrays (one for each int part of theta)
+ - the third level is made by a variable number of arrays
+ (one for each azimutal sector)
+
+ Arguments:
+ 1) the number of azimutal sectors
+
+ Operations:
+ - calculates the width of each sector, from the argument
+ - allocates and initializes all array levels
+ - sets a flag which tells the user if this NECESSARY operation
+ has been performed (it is needed BEFORE performing tracking)
+
+ **************************************************************************/
+{
+ // Set the number of sectors and their width.
+ fSectorNum = nsectors;
+ fSectorWidth = 360.0 / (Double_t)fSectorNum;
+ if (fMsgLevel >= 2) {
+ cout << fSectorNum << " sectors --> sector width (degrees) = " << fSectorWidth << endl;
+ }
+
+ // Meaningful indexes
+ Int_t ilayer, isector, itheta;
+
+ // Mark for the created objects
+ TObjArray *sector = 0;
+
+ // First index: layer
+ fNodes = new TObjArray**[fNLayers];
+ for (ilayer = 0; ilayer < fNLayers; ilayer++) {
+ fNodes[ilayer] = new TObjArray*[180];
+ for (itheta = 0; itheta < 180; itheta++) fNodes[ilayer][itheta] = 0;
+ for (itheta = 0; itheta < 180; itheta++) {
+ fNodes[ilayer][itheta] = new TObjArray(nsectors);
+ for (isector = 0; isector < nsectors; isector++) {
+ sector = new TObjArray;
+ sector->SetOwner();
+ fNodes[ilayer][itheta]->AddAt(sector, isector);
+ }
+ }
+ }
+
+ // Sets a checking flag to TRUE.
+ // This flag is checked before filling up the arrays with the points.
+ fStructureOK = kTRUE;
+}
+
+//__________________________________________________________________________________
+Int_t AliITStrackerANN::ArrangePoints(char *exportFile)
+
+/**************************************************************************
+
+ POINTS LOCATOR
+
+ This function assembles the operation from the other above methods,
+ and fills the arrays with the clusters already stored in the
+ layers of the tracker.
+ Then, in order to use this method, the user MUSTs call LoadClusters()
+ before.
+
+ Arguments:
+ 1) string for a file name where the global ccordinates
+ of all points can be exported (optional).
+ If this file must not be created, simply pass a NULL argument
+
+ Operations:
+ - for each AliITSclusterV2 in each AliITSlayer, a ne AliITSnode
+ is created and stored in the correct location.
+
+ Return values:
+ - the number of stored points
+ - when errors occur, or no points are found, 0 is returned
+
+ **************************************************************************/
+{
+ // Check if the array structure has been created
+ if (!fStructureOK) {
+ Error("ArrangePoints", "Structure NOT defined. Call CreateArrayStructure() first");
+ return 0;
+ }
+
+ // meaningful indexes
+ Int_t ientry, ilayer, nentries = 0, index;
+ Int_t nPtsLayer = 0;
+
+ // if the argument is not NULL, a file is opened
+ fstream file(exportFile, ios::out);
+ if (!exportFile || file.fail()) {
+ file.close();
+ exportFile = 0;
+ }
+
+ // scan all layers for node creation
+ for (ilayer = 0; ilayer < fNLayers; ilayer++) {
+ nPtsLayer = GetNumberOfClustersLayer(ilayer);
+ if (fMsgLevel >= 1) {
+ cout << "Layer " << ilayer << " --> " << nPtsLayer << " clusters" << endl;
+ }
+ for (ientry = 0; ientry < nPtsLayer; ientry++) {
+ // calculation of cluster index : (Bit mask LLLLIIIIIIIIIIII)
+ // (L = bits used for layer)
+ // (I = bits used for position in layer)
+ index = ilayer << 28;
+ index += ientry;
+ // add new AliITSnode object
+ AliITSnode *n = AddNode(index);
+ if ( (n != NULL) && exportFile ) {
+ file << index << ' ' << n->X() << ' ' << n->Y() << ' ' << n->Z() << endl;
+ }
+ }
+ nentries += nPtsLayer;
+ }
+
+ // a conventional final message is put at the end of file
+ if (exportFile) {
+ file << "-1 0.0 0.0 0.0" << endl;
+ file.close();
+ }
+
+ // returns the number of points processed
+ return nentries;
+}
+
+//__________________________________________________________________________________
+void AliITStrackerANN::StoreOverallMatches()
+
+/**************************************************************************
+
+ NODE-MATCH ANALYSIS
+
+ Once the nodes have been created, a firs analysis is to check
+ what pairs will satisfy at least the 'fixed' cuts (theta, helix-match)
+ and the most permissive (= larger) curvature cut.
+ All these node pairs are suitable for neuron creation.
+ In fact, when performing a Neural Tracking step, the only further check
+ will be a check against the current curvature step, while the other
+ are always the same.
+ For thi purpose, each AliITSnode has a data member, named 'fMatches'
+ which contains references to all other AliITSnodes in the successive layer
+ that form, with it, a 'good' pair, with respect to the above cited cuts.
+ Then, in each step for neuron creation, the possible neurons starting from
+ each node will be searched ONLY within the nodes referenced in fMatches.
+ This, of course, speeds up a lot the neuron creation procedure, at the
+ cost of some memory occupation, which results not to be critical.
+
+ Operations:
+ - for each AliITSnode, matches are found according to the criteria
+ expressed above, and stored in the node->fMatches array
+
+ **************************************************************************/
+{
+ // meaningful counters
+ Int_t ilayer, isector, itheta1, itheta2, check;
+ TObjArray *list1 = 0, *list2 = 0;
+ AliITSnode *node1 = 0, *node2 = 0;
+ Double_t thetaMin, thetaMax;
+ Int_t imin, imax;
+
+ // Scan for each sector
+ for (isector = 0; isector < fSectorNum; isector++) {
+ // sector is chosen once for both lists
+ for (ilayer = 0; ilayer < fNLayers - 1; ilayer++) {
+ for (itheta1 = 0; itheta1 < 180; itheta1++) {
+ list1 = (TObjArray*)fNodes[ilayer][itheta1]->At(isector);
+ TObjArrayIter iter1(list1);
+ while ( (node1 = (AliITSnode*)iter1.Next()) ) {
+ if (node1->IsUsed()) continue;
+ // clear an eventually already present array
+ // node1->Matches()->Clear();
+ // get the global coordinates and defines the theta interval from cut
+ thetaMin = (node1->GetTheta() * 180.0 / fgkPi) - fPolarInterval;
+ thetaMax = (node1->GetTheta() * 180.0 / fgkPi) + fPolarInterval;
+ imin = (Int_t)thetaMin;
+ imax = (Int_t)thetaMax;
+ if (imin < 0) imin = 0;
+ if (imax > 179) imax = 179;
+ // loop on the selected theta slots
+ for (itheta2 = imin; itheta2 <= imax; itheta2++) {
+ list2 = (TObjArray*)fNodes[ilayer + 1][itheta2]->At(isector);
+ TObjArrayIter iter2(list2);
+ while ( (node2 = (AliITSnode*)iter2.Next()) ) {
+ check = PassAllCuts(node1, node2, fCurvNum - 1, fVertexX, fVertexY, fVertexZ);
+ if (check == 0) {
+ node1->Matches()->AddLast(node2);
+ }
+ } // while (node2...)
+ } // for (itheta2...)
+ } // while (node1...)
+ } // for (itheta...)
+ } // for (ilayer...)
+ } // for (isector...)
+}
+
+//__________________________________________________________________________________
+Int_t AliITStrackerANN::PassAllCuts
+(AliITSnode *inner, AliITSnode *outer, Int_t curvStep,
+ Double_t vx, Double_t vy, Double_t vz)
+{
+// ***********************************************************************************
+//
+// This check is called in the above method for finding the matches of each node
+// It check the passed point pair against all the fixed cuts and a specified
+// curvature cut, among all the ones which have been defined.
+// The cuts need a vertex-constraint, which is not absolute, but it is passed
+// as argument.
+//
+// Arguments:
+// 1) the point in the inner layer
+// 2) the point in the outer layer
+// 3) curvature step for the curvature cut check (preferably the last)
+// 4) X of the used vertex
+// 5) Y of the used vertex
+// 6) Z of the used vertex
+//
+// Operations:
+// - if necessary, swaps the two points
+// (the first must be in the innermost of the two layers)
+// - checks for the theta cuts
+// - calculates the circle passing through the vertex
+// and the given points and checks for the curvature cut
+// - using the radius calculated there, checks for the helix-math cut
+//
+// Return values:
+// 0 - All cuts passed
+// 1 - theta 2D cut not passed
+// 2 - theta 3D cut not passed
+// 3 - curvature calculated but cut not passed
+// 4 - curvature not calculated (division by zero)
+// 5 - helix cut not passed
+// 6 - curvature inxed out of range
+//
+// ***********************************************************************************
+
+ // Check for curvature index
+ if (curvStep < 0 || curvStep >= fCurvNum) return 6;
+
+ // Swap points in order that r1 < r2
+ AliITSnode *temp = 0;
+ if (outer->GetLayer() < inner->GetLayer()) {
+ temp = outer;
+ outer = inner;
+ inner = temp;
+ }
+
+ // All cuts are variable according to the layer of the
+ // innermost point (the other point will surely be
+ // in the further one, because we don't check poin pairs
+ // which are not in adjacent layers)
+ // The reference is given by the innermost point.
+ Int_t layRef = inner->GetLayer();
+
+ // The calculations in the transverse plane are made in
+ // a shifted reference frame, whose origin corresponds to
+ // the reference point passed in the argument.
+ Double_t xIn = inner->X() - vx;
+ Double_t xOut = outer->X() - vx;
+ Double_t yIn = inner->Y() - vy;
+ Double_t yOut = outer->Y() - vy;
+ Double_t zIn = inner->Z() - vz;
+ Double_t zOut = outer->Z() - vz;
+ Double_t rIn = TMath::Sqrt(xIn*xIn + yIn*yIn);
+ Double_t rOut = TMath::Sqrt(xOut*xOut + yOut*yOut);
+
+ // Check for theta cut.
+ // There are two different angular cuts:
+ // one w.r. to the angle in the 2-dimensional r-z plane...
+ Double_t dthetaRZ;
+ TVector3 origin2innerRZ(zIn, rIn, 0.0);
+ TVector3 inner2outerRZ(zOut - zIn, rOut - rIn, 0.0);
+ dthetaRZ = origin2innerRZ.Angle(inner2outerRZ) * 180.0 / fgkPi;
+ if (dthetaRZ > fThetaCut2D[layRef]) {
+ return 1;
+ // r-z theta cut not passed ---> 1
+ }
+ // ...and another w.r. to the angle in the 3-dimensional x-y-z space
+ Double_t dthetaXYZ;
+ TVector3 origin2innerXYZ(xIn, yIn, zIn);
+ TVector3 inner2outerXYZ(xOut - xIn, yOut - yIn, zOut - zIn);
+ dthetaXYZ = origin2innerXYZ.Angle(inner2outerXYZ) * 180.0 / fgkPi;
+ if (dthetaXYZ > fThetaCut3D[layRef]) {
+ return 2;
+ // x-y-z theta cut not passed ---> 2
+ }
+
+ // Calculation & check of curvature
+ Double_t dx = xIn - xOut;
+ Double_t dy = yIn - yOut;
+ Double_t num = 2.0 * (xIn*yOut - xOut*yIn);
+ Double_t den = rIn*rOut*sqrt(dx*dx + dy*dy);
+ Double_t curv = 0.;
+ if (den != 0.) {
+ curv = TMath::Abs(num / den);
+ if (curv > fCurvCut[curvStep]) {
+ return 3;
+ // curvature too large for cut ---> 3
+ }
+ }
+ else {
+ Error("PassAllCuts", "Curvature calculations gives zero denominator");
+ return 4;
+ // error: denominator = 0 ---> 4
+ }
+
+ // Calculation & check of helix matching
+ Double_t helMatch = 0.0;
+ Double_t arcIn = 2.0 * rIn * curv;
+ Double_t arcOut = 2.0 * rOut * curv;
+ if (arcIn > -1.0 && arcIn < 1.0)
+ arcIn = TMath::ASin(arcIn);
+ else
+ arcIn = ((arcIn > 0.0) ? 0.5 : 1.5) * TMath::Pi();
+ if (arcOut > -1.0 && arcOut < 1.0)
+ arcOut = TMath::ASin(arcOut);
+ else
+ arcOut = ((arcOut > 0.0) ? 0.5 : 1.5) * TMath::Pi();
+ arcIn /= 2.0 * curv;
+ arcOut /= 2.0 * curv;
+ if (arcIn == 0.0 || arcOut == 0.0) {
+ Error("PassAllCuts", "Calculation returns zero-length arcs: l1=%f, l2=%f", arcIn, arcOut);
+ return 4;
+ // error: circumference arcs seem to equal zero ---> 4
+ }
+ helMatch = TMath::Abs(zIn / arcIn - zOut / arcOut);
+ if (helMatch > fHelixMatchCut[layRef]) {
+ return 5;
+ // helix match cut not passed ---> 5
+ }
+
+ // ALL CUTS PASSED ---> 0
+ return 0;
+}
+
+//__________________________________________________________________________________
+Bool_t AliITStrackerANN::PassCurvCut
+(AliITSnode *inner, AliITSnode *outer,
+ Int_t curvStep,
+ Double_t vx, Double_t vy, Double_t vz)
+{
+//***********************************************************************************
+//
+// This method operates essentially like the above one, but it is used
+// during a single step of Neural Tracking, where the curvature cut
+// changes.
+// Then, not necessaryly all the nodes stored in the fMatches array
+// will be suitable for neuron creation in an intermediate step.
+//
+// It has the same arguments of the PassAllCuts() method, but
+// the theta cut is not checked.
+// Moreover, it has a boolean return value.
+//
+//***********************************************************************************
+
+ // Check for curvature index
+ if (curvStep < 0 || curvStep >= fCurvNum) return 6;
+
+ // Find the reference layer
+ Int_t layIn = inner->GetLayer();
+ Int_t layOut = outer->GetLayer();
+ Int_t layRef = (layIn < layOut) ? layIn : layOut;
+
+ // The calculations in the transverse plane are made in
+ // a shifted reference frame, whose origin corresponds to
+ // the reference point passed in the argument.
+ Double_t xIn = inner->X() - vx;
+ Double_t xOut = outer->X() - vx;
+ Double_t yIn = inner->Y() - vy;
+ Double_t yOut = outer->Y() - vy;
+ Double_t zIn = inner->Z() - vz;
+ Double_t zOut = outer->Z() - vz;
+ Double_t rIn = TMath::Sqrt(xIn*xIn + yIn*yIn);
+ Double_t rOut = TMath::Sqrt(xOut*xOut + yOut*yOut);
+
+ // Calculation & check of curvature
+ Double_t dx = xIn - xOut;
+ Double_t dy = yIn - yOut;
+ Double_t num = 2.0 * (xIn*yOut - xOut*yIn);
+ Double_t den = rIn*rOut*sqrt(dx*dx + dy*dy);
+ Double_t curv = 0.;
+ /* OLD VERSION
+ if (den != 0.) {
+ curv = TMath::Abs(num / den);
+ if (curv > fCurvCut[curvStep]) return kFALSE;
+ return kTRUE;
+ }
+ else
+ return kFALSE;
+ */
+ // NEW VERSION
+ if (den != 0.) {
+ curv = TMath::Abs(num / den);
+ if (curv > fCurvCut[curvStep]) {
+ return kFALSE;
+ }
+ }
+ else {
+ Error("PassAllCuts", "Curvature calculations gives zero denominator");
+ return kFALSE;
+ }
+
+ // Calculation & check of helix matching
+ Double_t helMatch = 0.0;
+ Double_t arcIn = 2.0 * rIn * curv;
+ Double_t arcOut = 2.0 * rOut * curv;
+ if (arcIn > -1.0 && arcIn < 1.0)
+ arcIn = TMath::ASin(arcIn);
+ else
+ arcIn = ((arcIn > 0.0) ? 0.5 : 1.5) * TMath::Pi();
+ if (arcOut > -1.0 && arcOut < 1.0)
+ arcOut = TMath::ASin(arcOut);
+ else
+ arcOut = ((arcOut > 0.0) ? 0.5 : 1.5) * TMath::Pi();
+ arcIn /= 2.0 * curv;
+ arcOut /= 2.0 * curv;
+ if (arcIn == 0.0 || arcOut == 0.0) {
+ Error("PassAllCuts", "Calculation returns zero-length arcs: l1=%f, l2=%f", arcIn, arcOut);
+ return 4;
+ // error: circumference arcs seem to equal zero ---> 4
+ }
+ helMatch = TMath::Abs(zIn / arcIn - zOut / arcOut);
+ return (helMatch <= fHelixMatchCut[layRef]);
+}
+
+//__________________________________________________________________________________
+void AliITStrackerANN::PrintMatches(Bool_t stop)
+{
+// Prints the list of points which appear to match
+// each one of them, according to the preliminary
+// overall cuts.
+// The arguments states if a pause is required after printing
+// the matches for each one. In this case, a keypress is needed.
+
+ TObjArray *sector = 0;
+ Int_t ilayer, isector, itheta, nF;
+ AliITSnode *node1 = 0, *node2 = 0;
+ //AliITSclusterV2 *cluster1 = 0, *cluster2 = 0;
+
+ for (ilayer = 0; ilayer < 6; ilayer++) {
+ for (isector = 0; isector < fSectorNum; isector++) {
+ for (itheta = 0; itheta < 180; itheta++) {
+ sector = (TObjArray*)fNodes[ilayer][itheta]->At(isector);
+ TObjArrayIter points(sector);
+ while ( (node1 = (AliITSnode*)points.Next()) ) {
+ nF = (Int_t)node1->Matches()->GetEntries();
+ cout << "Node layer: " << node1->GetLayer() << " --> ";
+ if (!nF) {
+ cout << "NO Matches!!!" << endl;
+ continue;
+ }
+ cout << nF << " Matches" << endl;
+ cout << "Reference cluster: " << hex << node1->ClusterRef() << endl;
+ TObjArrayIter matches(node1->Matches());
+ while ( (node2 = (AliITSnode*)matches.Next()) ) {
+ cout << "Match with " << hex << node2->ClusterRef() << endl;
+ }
+ if (stop) {
+ cout << "Press a key" << endl;
+ cin.get();
+ }
+ }
+ }
+ }
+ }
+}
+
+//__________________________________________________________________________________
+void AliITStrackerANN::ResetNodes(Int_t isector)
+{
+/***********************************************************************************
+
+ NODE NEURON ARRAY CLEANER
+
+ After a neural tracking step, this method
+ clears the arrays 'fInnerOf' and 'fOuterOf' of each AliITSnode
+
+ Arguments:
+ - the sector where the operation is being executed
+
+ ***********************************************************************************/
+
+ Int_t ilayer, itheta;
+ TObjArray *sector = 0;
+ AliITSnode *node = 0;
+ for (ilayer = 0; ilayer < fNLayers; ilayer++) {
+ for (itheta = 0; itheta < 180; itheta++) {
+ sector = (TObjArray*)fNodes[ilayer][itheta]->At(isector);
+ TObjArrayIter iter(sector);
+ for (;;) {
+ node = (AliITSnode*)iter.Next();
+ if (!node) break;
+ node->InnerOf()->Clear();
+ node->OuterOf()->Clear();
+ /*
+ delete node->InnerOf();
+ delete node->OuterOf();
+ node->InnerOf() = new TObjArray;
+ node->OuterOf() = new TObjArray;
+ */
+ }
+ }
+ }
+}
+
+//__________________________________________________________________________________
+Int_t AliITStrackerANN::CreateNeurons
+(AliITSnode *node, Int_t curvStep, Double_t vx, Double_t vy, Double_t vz)
+{
+ // This method is used to create alle suitable neurons starting from
+ // a single AliITSnode. Each unit is also stored in the fInnerOf array
+ // of the passed node, and in the fOuterOf array of the other neuron edge.
+ // In the new implementation of the intermediate check steps, a further one
+ // is made, which chechs how well a helix is matched by three points
+ // in three consecutive layers.
+ // Then, a vertex-constrained check is made with vertex located
+ // in a layer L, for points in layers L+1 and L+2.
+ //
+ // In order to do this, the creator works recursively, in a tree-visit like operation.
+ // The neurons are effectively created only if the node argument passed is in
+ // the 5th layer (they are created between point of 5th and 6th layer).
+ // If the node is in an inner layer, its coordinates are passet as vertex for a nested
+ // call of the same function in the next two layers.
+ //
+ // Arguments:
+ // 1) reference node
+ // 2) current curvature cut step
+ // 3) X of referenced temporary (not primary) vertex
+ // 4) Y of referenced temporary (not primary) vertex
+ // 5) Z of referenced temporary (not primary) vertex
+ //
+ // Operations:
+ // - if the layer is the 5th, neurons are created with nodes
+ // in the fMatches array of the passed node
+ // - otherwise, the X, Y, Z of the passed node are given as
+ // vertex and the same procedure is recursively called for all
+ // nodes in the fMatches array of the passed one.
+ //
+ // Return values:
+ // - the total number of neurons created from the passed one
+ // summed with all neurons created from all nodes well matched with it
+ // (assumes a meaning only for nodes in the first layer)
+
+ // local variables
+ Int_t made = 0; // counter
+ Bool_t found = 0; // flag
+ AliITSnode *match = 0; // cursor for a AliITSnode array
+ AliITSneuron *unit = 0; // cursor for a AliITSneuron array
+
+ // --> Case 0: the passed node has already been used
+ // as member of a track found in a previous step.
+ // In this case, of course, the function exits.
+ if (node->IsUsed()) return 0;
+
+ // --> Case 1: there are NO well-matched points.
+ // This can happen in all ITS layers, but it happens **for sure**
+ // for a node in the 'last' layer.
+ // Even in this case, the function exits.
+ if (node->Matches()->IsEmpty()) return 0;
+
+ // --> Case 2: there are well-matched points.
+ // In this case, the function creates a neuron for each
+ // well-matched pair (according to the cuts for the current step)
+ // Moreover, before storing the neuron, a check is necessary
+ // to avoid the duplicate creation of the same neuron twice.
+ // (This could happen if the 3 last arguments of the function
+ // are close enough to cause a good match for the current step
+ // between two points, independently of their difference).
+ // Finally, a node is skipped if it has already been used.
+ // For each matched point for which a neuron is created, the procedure is
+ // recursively called.
+ TObjArrayIter matches(node->Matches());
+ while ( (match = (AliITSnode*)matches.Next()) ) {
+ if (match->IsUsed()) continue;
+ if (!PassCurvCut(node, match, curvStep, vx, vy, vz)) continue;
+ found = kFALSE;
+ if (!node->InnerOf()->IsEmpty()) {
+ TObjArrayIter presentUnits(node->InnerOf());
+ while ( (unit = (AliITSneuron*)presentUnits.Next()) ) {
+ if (unit->Inner() == node && unit->Outer() == match) {
+ found = kTRUE;
+ break;
+ }
+ }
+ }
+ if (found) continue;
+ AliITSneuron *unit = new AliITSneuron(node, match, fEdge2, fEdge1);
+ fNeurons->AddLast(unit);
+ node->InnerOf()->AddLast(unit);
+ match->OuterOf()->AddLast(unit);
+ made += CreateNeurons(match, curvStep, node->X(), node->Y(), node->Z());
+ made++;
+ }
+
+ // Of course, the return value contains the number of neurons
+ // counting in also the oned created in all levels of recursive calls.
+ return made;
+}
+
+//__________________________________________________________________________________
+Int_t AliITStrackerANN::CreateNetwork(Int_t sector, Int_t curvStep)
+{
+ // This function simply recalls the CreateNeurons() method for each node
+ // in the first layer, for the current sector.
+ // This generates the whole network, thanks to the recursive calls.
+ //
+ // Arguments:
+ // 1) current sector
+ // 2) current curvature step
+ //
+ // Operations:
+ // - scans the nodes array for all theta's in the current sector
+ // and layer 0, and calls the CreateNeurons() function.
+
+ // removes all eventually present neurons
+ if (fNeurons) delete fNeurons;
+ fNeurons = new TObjArray;
+ fNeurons->SetOwner(kTRUE);
+
+ // calls the ResetNodes() function to free the AliITSnode arrays
+ if (fMsgLevel >= 2) {
+ cout << "Sector " << sector << " PHI = ";
+ cout << fSectorWidth * (Double_t)sector << " --> ";
+ cout << fSectorWidth * (Double_t)(sector + 1) << endl;
+ cout << "Curvature step " << curvStep << " [cut = " << fCurvCut[curvStep] << "]" << endl;
+ }
+ ResetNodes(sector);
+
+ // meaningful counters
+ Int_t itheta, neurons = 0;
+ TObjArray *lstSector = 0;
+
+ // NEW VERSION
+ Double_t vx[6], vy[6], vz[6];
+ AliITSnode *p[6] = {0, 0, 0, 0, 0, 0};
+ for (itheta = 0; itheta < 180; itheta++) {
+ lstSector = (TObjArray*)fNodes[0][itheta]->At(sector);
+ TObjArrayIter lay0(lstSector);
+ while ( (p[0] = (AliITSnode*)lay0.Next()) ) {
+ if (p[0]->IsUsed()) continue;
+ vx[0] = fVertexX;
+ vy[0] = fVertexY;
+ vz[0] = fVertexZ;
+ neurons += CreateNeurons(p[0], curvStep, fVertexX, fVertexY, fVertexZ);
+ /*
+ TObjArrayIter lay1(p[0]->Matches());
+ while ( (p[1] = (AliITSnode*)lay1.Next()) ) {
+ if (p[1]->IsUsed()) continue;
+ if (!PassCurvCut(p[0], p[1], curvStep, vx[0], vy[0], vz[0])) continue;
+ unit = new AliITSneuron;
+ unit->Inner() = p[0];
+ unit->Outer() = p[1];
+ unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
+ unit->fGain = new TObjArray;
+ fNeurons->AddLast(unit);
+ p[0]->InnerOf()->AddLast(unit);
+ p[1]->OuterOf()->AddLast(unit);
+ neurons++;
+ vx[1] = p[0]->X();
+ vy[1] = p[0]->Y();
+ vz[1] = p[0]->Z();
+ TObjArrayIter lay2(p[1]->Matches());
+ while ( (p[2] = (AliITSnode*)lay2.Next()) ) {
+ if (p[2]->IsUsed()) continue;
+ if (!PassCurvCut(p[1], p[2], curvStep, vx[1], vy[1], vz[1])) continue;
+ unit = new AliITSneuron;
+ unit->Inner() = p[1];
+ unit->Outer() = p[2];
+ unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
+ unit->fGain = new TObjArray;
+ fNeurons->AddLast(unit);
+ p[1]->InnerOf()->AddLast(unit);
+ p[2]->OuterOf()->AddLast(unit);
+ neurons++;
+ vx[2] = p[1]->X();
+ vy[2] = p[1]->Y();
+ vz[2] = p[1]->Z();
+ TObjArrayIter lay3(p[2]->Matches());
+ while ( (p[3] = (AliITSnode*)lay3.Next()) ) {
+ if (p[3]->IsUsed()) continue;
+ if (!PassCurvCut(p[2], p[3], curvStep, vx[2], vy[2], vz[2])) continue;
+ unit = new AliITSneuron;
+ unit->Inner() = p[2];
+ unit->Outer() = p[3];
+ unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
+ unit->fGain = new TObjArray;
+ fNeurons->AddLast(unit);
+ p[2]->InnerOf()->AddLast(unit);
+ p[3]->OuterOf()->AddLast(unit);
+ neurons++;
+ vx[3] = p[2]->X();
+ vy[3] = p[2]->Y();
+ vz[3] = p[2]->Z();
+ TObjArrayIter lay4(p[3]->Matches());
+ while ( (p[4] = (AliITSnode*)lay4.Next()) ) {
+ if (p[4]->IsUsed()) continue;
+ if (!PassCurvCut(p[3], p[4], curvStep, vx[3], vy[3], vz[3])) continue;
+ unit = new AliITSneuron;
+ unit->Inner() = p[3];
+ unit->Outer() = p[4];
+ unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
+ unit->fGain = new TObjArray;
+ fNeurons->AddLast(unit);
+ p[3]->InnerOf()->AddLast(unit);
+ p[4]->OuterOf()->AddLast(unit);
+ neurons++;
+ vx[4] = p[3]->X();
+ vy[4] = p[3]->Y();
+ vz[4] = p[3]->Z();
+ TObjArrayIter lay5(p[4]->Matches());
+ while ( (p[5] = (AliITSnode*)lay5.Next()) ) {
+ if (p[5]->IsUsed()) continue;
+ if (!PassCurvCut(p[4], p[5], curvStep, vx[4], vy[4], vz[4])) continue;
+ unit = new AliITSneuron;
+ unit->Inner() = p[4];
+ unit->Outer() = p[5];
+ unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
+ unit->fGain = new TObjArray;
+ fNeurons->AddLast(unit);
+ p[4]->InnerOf()->AddLast(unit);
+ p[5]->OuterOf()->AddLast(unit);
+ neurons++;
+ } // while (p[5])
+ } // while (p[4])
+ } // while (p[3])
+ } // while (p[2])
+ } // while (p[1])
+ */
+ } // while (p[0])
+ } // for (itheta...)
+ // END OF NEW VERSION
+
+ /* OLD VERSION
+ for (ilayer = 0; ilayer < 6; ilayer++) {
+ for (itheta = 0; itheta < 180; itheta++) {
+ lstSector = (TObjArray*)fNodes[ilayer][itheta]->At(sector_idx);
+ TObjArrayIter inners(lstSector);
+ while ( (inner = (AliITSnode*)inners.Next()) ) {
+ if (inner->GetUser() >= 0) continue;
+ TObjArrayIter outers(inner->Matches());
+ while ( (outer = (AliITSnode*)outers.Next()) ) {
+ if (outer->GetUser() >= 0) continue;
+ if (!PassCurvCut(inner, outer, curvStep, fVX, fVY, fVZ)) continue;
+ unit = new AliITSneuron;
+ unit->Inner() = inner;
+ unit->Outer() = outer;
+ unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
+ unit->fGain = new TObjArray;
+ fNeurons->AddLast(unit);
+ inner->InnerOf()->AddLast(unit);
+ outer->OuterOf()->AddLast(unit);
+ neurons++;
+ } // for (;;)
+ } // for (;;)
+ } // for (itheta...)
+ } // for (ilayer...)
+ */
+
+ fNeurons->SetOwner();
+ return neurons;
+ }
+
+ //__________________________________________________________________________________
+Int_t AliITStrackerANN::LinkNeurons()
+{
+/***********************************************************************************
+
+ SYNAPSIS GENERATOR
+
+ Scans the whole neuron array, in order to find all neuron pairs
+ which are connected in sequence and share a positive weight.
+ For each of them, an AliITSlink is created, which stores
+ the weight value, and will allow for a faster calculation
+ of the total neural input for each updating cycle.
+
+ Every neuron contains an object array which stores all AliITSlink
+ objects which point to sequenced units, with the respective weights.
+
+ Return value:
+ - the number of link created for the neural network.
+ If they are 0, no updating can be done and the step is skipped.
+
+ ***********************************************************************************/
+
+ // meaningful indexes
+ Int_t total = 0;
+ Double_t weight = 0.0;
+ TObjArrayIter neurons(fNeurons), *iter;
+ AliITSneuron *neuron = 0, *test = 0;
+
+ // scan in the neuron array
+ for (;;) {
+ neuron = (AliITSneuron*)neurons.Next();
+ if (!neuron) break;
+ // checks for neurons startin from its head ( -> )
+ iter = (TObjArrayIter*)neuron->Inner()->OuterOf()->MakeIterator();
+ for (;;) {
+ test = (AliITSneuron*)iter->Next();
+ if (!test) break;
+ weight = Weight(test, neuron);
+ if (weight > 0.0) neuron->Gain()->AddLast(new AliITSlink(weight, test));
+ total++;
+ }
+ delete iter;
+ // checks for neurons ending in its tail ( >- )
+ iter = (TObjArrayIter*)neuron->Outer()->InnerOf()->MakeIterator();
+ for (;;) {
+ test = (AliITSneuron*)iter->Next();
+ if (!test) break;
+ weight = Weight(neuron, test);
+ if (weight > 0.0) neuron->Gain()->AddLast(new AliITSlink(weight, test));
+ total++;
+ }
+ delete iter;
+ }
+ return total;
+}
+
+//__________________________________________________________________________________
+Bool_t AliITStrackerANN::Update()
+{
+/***********************************************************************************
+
+ Performs a single updating cycle.
+
+ Operations:
+ - for each neuron, gets the activation with the neuron Activate() method
+ - checks if stability has been reached (compare mean activation variation
+ with the stability threshold data member)
+
+ Return values:
+ - kTRUE means that the neural network has stabilized
+ - kFALSE means that another updating cycle is needed
+
+ ***********************************************************************************/
+
+ Double_t actVar = 0.0, totDiff = 0.0;
+ TObjArrayIter iter(fNeurons);
+ AliITSneuron *unit;
+ for (;;) {
+ unit = (AliITSneuron*)iter.Next();
+ if (!unit) break;
+ actVar = unit->Activate(fTemperature);
+ // calculation the relative activation variation
+ totDiff += actVar;
+ }
+ totDiff /= fNeurons->GetSize();
+ return (totDiff < fStabThreshold);
+}
+
+//__________________________________________________________________________________
+void AliITStrackerANN::FollowChains(Int_t sector)
+{
+/***********************************************************************************
+
+ CHAINS CREATION
+
+ After that the neural network has stabilized,
+ the final step is to create polygonal chains
+ of clusters, one in each layer, which represent
+ the tracks recognized by the neural algorithm.
+ This is made by means of a choice of the best
+ neuron among the ones starting from each point.
+
+ Once that such neuron is selected, its inner point
+ will set the 'fNext' field to its outer point, and
+ similarly, its outer point will set the 'fPrev' field
+ to its inner point.
+ This defines a bi-directional sequence.
+
+ In this procedure, it can happen that many neurons
+ which have the head of the arrow in a given node, will
+ all select as best following the neuron with the largest
+ activation starting in that point.
+ This results in MANY nodes which have the same 'fNext'.
+ But, this field will be set to NULL for all these points,
+ but the only one which is pointed by the 'fPrev' field
+ of this shared node.
+
+ ***********************************************************************************/
+
+ // meaningful counters
+ Int_t itheta, ilayer;
+ TObjArray *lstSector = 0;
+ Double_t test = fActMinimum;
+ AliITSnode *p = 0;
+ AliITSneuron *n = 0;
+
+ // scan the whole collection of nodes
+ for (ilayer = 0; ilayer < fNLayers; ilayer++) {
+ for (itheta = 0; itheta < 180; itheta++) {
+ // get the array of point in a given layer/theta-slot/sector
+ lstSector = (TObjArray*)fNodes[ilayer][itheta]->At(sector);
+ TObjArrayIter nodes(lstSector);
+ while ( (p = (AliITSnode*)nodes.Next()) ) {
+ // if the point is used, it is skipped
+ if (p->IsUsed()) continue;
+ // initially, fNext points to nothing, and
+ // the comparison value is set to the minimum activation
+ // which allows to say that a neuron is turned 'on'
+ // a node from which only 'off' neurons start is probably
+ // a noise point, which will be excluded from the reconstruction.
+ test = fActMinimum;
+ p->Next() = 0;
+ TObjArrayIter innerof(p->InnerOf());
+ while ( (n = (AliITSneuron*)innerof.Next()) ) {
+ // if the examined neuron has not the largest activation
+ // it is skipped and removed from array of all neurons
+ // and of its outer point (its inner is the cursor p)
+ if (n->Activation() < test) {
+ p->InnerOf()->Remove(n);
+ n->Outer()->OuterOf()->Remove(n);
+ delete fNeurons->Remove(n);
+ continue;
+ }
+ // otherwise, its activation becomes the maximum reference
+ p->Next() = n->Outer();
+ // at the exit of the while(), the fNext will point
+ // to the outer node of the neuron starting in p, whose
+ // activation is the largest.
+ }
+ // the same procedure is made now for all neurons
+ // for which p is the outer point
+ test = fActMinimum;
+ p->Prev() = 0;
+ TObjArrayIter outerof(p->OuterOf());
+ while ( (n = (AliITSneuron*)outerof.Next()) ) {
+ // if the examined neuron has not the largest activation
+ // it is skipped and removed from array of all neurons
+ // and of its inner point (its outer is the cursor p)
+ if (n->Activation() < test) {
+ p->OuterOf()->Remove(n);
+ n->Inner()->InnerOf()->Remove(n);
+ delete fNeurons->Remove(n);
+ continue;
+ }
+ // otherwise, its activation becomes the maximum reference
+ p->Prev() = n->Inner();
+ // at the exit of the while(), the fPrev will point
+ // to the inner node of the neuron ending in p, whose
+ // activation is the largest.
+ }
+ } // end while (p ...)
+ } // end for (itheta ...)
+ } // end for (ilayer ...)
+
+ // now the mismatches are solved
+ Bool_t matchPrev, matchNext;
+ for (ilayer = 0; ilayer < fNLayers; ilayer++) {
+ for (itheta = 0; itheta < 180; itheta++) {
+ // get the array of point in a given layer/theta-slot/sector
+ lstSector = (TObjArray*)fNodes[ilayer][itheta]->At(sector);
+ TObjArrayIter nodes(lstSector);
+ while ( (p = (AliITSnode*)nodes.Next()) ) {
+ // now p will point to a fPrev and a fNext node.
+ // Ideally they are placed this way: fPrev --> P --> fNext
+ // A mismatch happens if the point addressed as fPrev does NOT
+ // point to p as its fNext. And the same for the point addressed
+ // as fNext.
+ // In this case, the fNext and fPrev pointers are set to NULL
+ // and p is excluded from the reconstruction
+ matchPrev = matchNext= kFALSE;
+ if (ilayer > 0 && p->Prev() != NULL)
+ if (p->Prev()->Next() == p) matchPrev = kTRUE;
+ if (ilayer < 5 && p->Next() != NULL)
+ if (p->Next()->Prev() == p) matchNext = kTRUE;
+ if (ilayer == 0)
+ matchPrev = kTRUE;
+ else if (ilayer == 5)
+ matchNext = kTRUE;
+ if (!matchNext || !matchPrev) {
+ p->Prev() = p->Next() = 0;
+ }
+ } // end while (p ...)
+ } // end for (itheta ...)
+ } // end for (ilayer ...)
+}
+
+//__________________________________________________________________________________
+Int_t AliITStrackerANN::SaveTracks(Int_t sector)
+{
+/********************************************************************************
+
+ TRACK SAVING
+ ------------
+ Using the fNext and fPrev pointers, the chain is followed
+ and the track is fitted and saved.
+ Of course, the track is followed as a chain with a point
+ for each layer, then the track following starts always
+ from the clusters in layer 0.
+
+***********************************************************************************/
+
+ // if not initialized, the tracks TobjArray is initialized
+ if (!fFoundTracks) fFoundTracks = new TObjArray;
+
+ // meaningful counters
+ Int_t itheta, ilayer, l;
+ TObjArray *lstSector = 0;
+ AliITSnode *p = 0, *q = 0, **node = new AliITSnode*[fNLayers];
+ for (l = 0; l < fNLayers; l++) node[l] = 0;
+
+ /*
+ array = new AliITSnode*[fNLayers + 1];
+ for (l = 0; l <= fNLayers; l++) array[l] = 0;
+ array[0] = new AliITSnode();
+ array[0]->X() = fVertexX;
+ array[0]->Y() = fVertexY;
+ array[0]->Z() = fVertexZ;
+ array[0]->ErrX2() = fVertexErrorX;
+ array[0]->ErrY2() = fVertexErrorY;
+ array[0]->ErrZ2() = fVertexErrorZ;
+ */
+ Double_t *param = new Double_t[8];
+
+ // scan the whole collection of nodes
+ for (ilayer = 0; ilayer < 1; ilayer++) {
+ for (itheta = 0; itheta < 180; itheta++) {
+ // get the array of point in a given layer/theta-slot/sector
+ lstSector = (TObjArray*)fNodes[ilayer][itheta]->At(sector);
+ TObjArrayIter nodes(lstSector);
+ while ( (p = (AliITSnode*)nodes.Next()) ) {
+ for (q = p; q; q = q->Next()) {
+ l = q->GetLayer();
+ node[l] = q;
+ }
+ //if (!RiemannFit(fNLayers, node, param)) continue;
+ // initialization of Kalman Filter Tracking
+ AliITSclusterV2 *cluster = (AliITSclusterV2*)GetCluster(node[0]->ClusterRef());
+ Int_t mod = cluster->GetDetectorIndex();
+ Int_t lay, lad, det;
+ fGeom->GetModuleId(mod, lay, lad, det);
+ Float_t y0 = cluster->GetY();
+ Float_t z0 = cluster->GetZ();
+ AliITStrackSA* trac = new AliITStrackSA(lay, lad, det,
+ y0, z0,
+ param[4], param[7], param[3], 1);
+ for (l = 0; l < fNLayers; l++) {
+ cluster = (AliITSclusterV2*)GetCluster(node[l]->ClusterRef());
+ if (cluster) trac->AddClusterV2(l, (node[l]->ClusterRef() & 0x0fffffff)>>0);
+ }
+ AliITStrackV2* ot = new AliITStrackV2(*trac);
+ ot->ResetCovariance();
+ ot->ResetClusters();
+ if (RefitAt(49.,ot,trac)) { //fit from layer 1 to layer 6
+ AliITStrackV2 *otrack2 = new AliITStrackV2(*ot);
+ otrack2->ResetCovariance();
+ otrack2->ResetClusters();
+ //fit from layer 6 to layer 1
+ if (RefitAt(3.7,otrack2,ot)) fFoundTracks->AddLast(otrack2);
+ }
+ // end of Kalman Filter fit
+ }
+ }
+ }
+
+ return 1;
+}
+
+//__________________________________________________________________________________
+void AliITStrackerANN::ExportTracks(const char *filename) const
+{
+// Exports found tracks into a TTree of AliITStrackV2 objects
+ TFile *file = new TFile(filename, "RECREATE");
+ TTree *tree = new TTree("TreeT-ANN", "Tracks found in ITS stand-alone with Neural Tracking");
+ AliITStrackV2 *track = 0;
+ tree->Branch("Tracks", &track, "AliITStrackV2");
+ TObjArrayIter tracks(fFoundTracks);
+ while ( (track = (AliITStrackV2*)tracks.Next()) ) {
+ tree->Fill();
+ }
+ file->cd();
+ tree->Write();
+ file->Close();
+}
+
+
+//__________________________________________________________________________________
+void AliITStrackerANN::CleanNetwork()
+{
+ // Removes deactivated units from the network
+
+ AliITSneuron *unit = 0;
+ TObjArrayIter neurons(fNeurons);
+ while ( (unit = (AliITSneuron*)neurons.Next()) ) {
+ if (unit->Activation() < fActMinimum) {
+ unit->Inner()->InnerOf()->Remove(unit);
+ unit->Outer()->OuterOf()->Remove(unit);
+ delete fNeurons->Remove(unit);
+ }
+ }
+ return;
+ Bool_t removed;
+ Int_t nIn, nOut;
+ AliITSneuron *enemy = 0;
+ neurons.Reset();
+ while ( (unit = (AliITSneuron*)neurons.Next()) ) {
+ nIn = (Int_t)unit->Inner()->InnerOf()->GetSize();
+ nOut = (Int_t)unit->Outer()->OuterOf()->GetSize();
+ if (nIn < 2 && nOut < 2) continue;
+ removed = kFALSE;
+ if (nIn > 1) {
+ TObjArrayIter competing(unit->Inner()->InnerOf());
+ while ( (enemy = (AliITSneuron*)competing.Next()) ) {
+ if (unit->Activation() > enemy->Activation()) {
+ enemy->Inner()->InnerOf()->Remove(enemy);
+ enemy->Outer()->OuterOf()->Remove(enemy);
+ delete fNeurons->Remove(enemy);
+ }
+ else {
+ unit->Inner()->InnerOf()->Remove(unit);
+ unit->Outer()->OuterOf()->Remove(unit);
+ delete fNeurons->Remove(unit);
+ removed = kTRUE;
+ break;
+ }
+ }
+ if (removed) continue;
+ }
+ if (nOut > 1) {
+ TObjArrayIter competing(unit->Outer()->OuterOf());
+ while ( (enemy = (AliITSneuron*)competing.Next()) ) {
+ if (unit->Activation() > enemy->Activation()) {
+ enemy->Inner()->InnerOf()->Remove(enemy);
+ enemy->Outer()->OuterOf()->Remove(enemy);
+ delete fNeurons->Remove(enemy);
+ }
+ else {
+ unit->Inner()->InnerOf()->Remove(unit);
+ unit->Outer()->OuterOf()->Remove(unit);
+ delete fNeurons->Remove(unit);
+ removed = kTRUE;
+ break;
+ }
+ }
+ }
+ }
+ }
+
+//__________________________________________________________________________________
+Int_t AliITStrackerANN::StoreTracks()
+{
+ // Stores the tracks found in a single neural tracking step.
+ // In order to do this, it sects each neuron which has a point
+ // in the first layer.
+ // Then
+
+ // if not initialized, the tracks TobjArray is initialized
+ if (!fFoundTracks) fFoundTracks = new TObjArray;
+
+ Int_t i, check, stored = 0;
+ Double_t testAct = 0;
+ AliITSneuron *unit = 0, *cursor = 0, *fwd = 0;
+ AliITSnode *node = 0;
+ TObjArrayIter iter(fNeurons), *fwdIter;
+ TObjArray *removedUnits = new TObjArray(0);
+ removedUnits->SetOwner(kFALSE);
+ AliITStrackANN annTrack(fNLayers);
+
+ for (;;) {
+ unit = (AliITSneuron*)iter.Next();
+ if (!unit) break;
+ if (unit->Inner()->GetLayer() > 0) continue;
+ annTrack.SetNode(unit->Inner()->GetLayer(), unit->Inner());
+ annTrack.SetNode(unit->Outer()->GetLayer(), unit->Outer());
+ node = unit->Outer();
+ removedUnits->AddLast(unit);
+ while (node) {
+ testAct = fActMinimum;
+ fwdIter = (TObjArrayIter*)node->InnerOf()->MakeIterator();
+ fwd = 0;
+ for (;;) {
+ cursor = (AliITSneuron*)fwdIter->Next();
+ if (!cursor) break;
+ if (cursor->Used()) continue;
+ if (cursor->Activation() >= testAct) {
+ testAct = cursor->Activation();
+ fwd = cursor;
+ }
+ }
+ if (!fwd) break;
+ removedUnits->AddLast(fwd);
+ node = fwd->Outer();
+ annTrack.SetNode(node->GetLayer(), node);
+ }
+ check = annTrack.CheckOccupation();
+ if (check >= 6) {
+ stored++;
+ // FIT
+ //if (!RiemannFit(fNLayers, trackitem, param)) continue;
+ if (!annTrack.RiemannFit()) continue;
+ // initialization of Kalman Filter Tracking
+ AliITSclusterV2 *cluster = (AliITSclusterV2*)GetCluster(annTrack[0]->ClusterRef());
+ Int_t mod = cluster->GetDetectorIndex();
+ Int_t lay, lad, det;
+ fGeom->GetModuleId(mod, lay, lad, det);
+ Float_t y0 = cluster->GetY();
+ Float_t z0 = cluster->GetZ();
+ AliITStrackSA* trac = new AliITStrackSA(lay, lad, det, y0, z0,
+ annTrack.Phi(), annTrack.TanLambda(),
+ annTrack.Curv(), 1);
+ for (Int_t l = 0; l < fNLayers; l++) {
+ if (!annTrack[l]) continue;
+ cluster = (AliITSclusterV2*)GetCluster(annTrack[l]->ClusterRef());
+ if (cluster) trac->AddClusterV2(l, (annTrack[l]->ClusterRef() & 0x0fffffff)>>0);
+ }
+ AliITStrackV2* ot = new AliITStrackV2(*trac);
+ ot->ResetCovariance();
+ ot->ResetClusters();
+ if (RefitAt(49.,ot,trac)) { //fit from layer 1 to layer 6
+ AliITStrackV2 *otrack2 = new AliITStrackV2(*ot);
+ otrack2->ResetCovariance();
+ otrack2->ResetClusters();
+ //fit from layer 6 to layer 1
+ if (RefitAt(3.7,otrack2,ot)) fFoundTracks->AddLast(otrack2);
+ }
+ // end of Kalman Filter fit
+ // END FIT
+ for (i = 0; i < fNLayers; i++) {
+ //node = (AliITSnode*)removedPoints->At(i);
+ //node->Use();
+ annTrack[i]->Use();
+ }
+ fwdIter = (TObjArrayIter*)removedUnits->MakeIterator();
+ for (;;) {
+ cursor = (AliITSneuron*)fwdIter->Next();
+ if(!cursor) break;
+ cursor->Used() = 1;
+ }
+ }
+ }
+
+ return stored;
+}
+
+Double_t AliITStrackerANN::Weight(AliITSneuron *nAB, AliITSneuron *nBC)
+{
+/***********************************************************************************
+ * Calculation of neural weight.
+ * The implementation of positive neural weight is set only in the case
+ * of connected units (e.g.: A->B with B->C).
+ * Given that B is the **common** point. care should be taken to pass
+ * as FIRST argument the neuron going "to" B, and
+ * as SECOND argument the neuron starting "from" B
+ * anyway, a check is put in order to return 0.0 when arguments are not well given.
+ ***********************************************************************************/
+
+ if (nAB->Outer() != nBC->Inner()) {
+ if (nBC->Outer() == nAB->Inner()) {
+ AliITSneuron *temp = nAB;
+ nAB = nBC;
+ nBC = temp;
+ temp = 0;
+ if (fMsgLevel >= 3) {
+ Info("Weight", "Switching wrongly ordered arguments.");
+ }
+ }
+ Warning("Weight", "Not connected segments. Returning 0.0");
+ return 0.0;
+ }
+
+ AliITSnode *pA = nAB->Inner();
+ AliITSnode *pB = nAB->Outer();
+ AliITSnode *pC = nBC->Outer();
+
+ TVector3 vAB(pB->X() - pA->X(), pB->Y() - pA->Y(), pB->Z() - pA->Z());
+ TVector3 vBC(pC->X() - pB->X(), pC->Y() - pB->Y(), pC->Z() - pB->Z());
+
+ Double_t weight = 1.0 - sin(vAB.Angle(vBC));
+ return fGain2CostRatio * TMath::Power(weight, fExponent);
+}
+
+
+
+/******************************************
+ ******************************************
+ *** AliITStrackerANN::AliITSnode class ***
+ ******************************************
+ ******************************************/
+
+//__________________________________________________________________________________
+inline AliITStrackerANN::AliITSnode::AliITSnode()
+: fUsed(kFALSE), fClusterRef(-1),
+ fMatches(NULL), fInnerOf(NULL), fOuterOf(NULL),
+ fNext(NULL), fPrev(NULL)
+{
+ // Constructor for the embedded 'AliITSnode' class.
+ // It initializes all pointer-like objects.
+
+ fX = fY = fZ = 0.0;
+ fEX2 = fEY2 = fEZ2 = 0.0;
+}
+
+//__________________________________________________________________________________
+AliITStrackerANN::AliITSnode::~AliITSnode()
+{
+ // Destructor for the embedded 'AliITSnode' class.
+ // It should clear the object arrays, but it is possible
+ // that some objects still are useful after the point deletion
+ // then the arrays are cleared but their objects are owed by
+ // another TCollection object, and not deleted.
+ // For safety reasons, all the pointers are set to zero.
+
+ fInnerOf->SetOwner(kFALSE);
+ fInnerOf->Clear();
+ delete fInnerOf;
+ fInnerOf = 0;
+ fOuterOf->SetOwner(kFALSE);
+ fOuterOf->Clear();
+ delete fOuterOf;
+ fOuterOf = 0;
+ fMatches->SetOwner(kFALSE);
+ fMatches->Clear();
+ delete fMatches;
+ fMatches = 0;
+ fNext = 0;
+ fPrev = 0;
+}
+
+//__________________________________________________________________________________
+inline Double_t AliITStrackerANN::AliITSnode::GetPhi() const
+{
+ // Calculates the 'phi' (azimutal) angle, and returns it
+ // in the range between 0 and 2Pi radians.
+
+ Double_t q;
+ q = TMath::ATan2(fY,fX);
+ if (q >= 0.)
+ return q;
+ else
+ return q + TMath::TwoPi();
+}
+
+//__________________________________________________________________________________
+inline Double_t AliITStrackerANN::AliITSnode::GetError(Option_t *option)
+{
+ // Returns the error or the square error of
+ // values related to the coordinates in different systems.
+ // The option argument specifies the coordinate error desired:
+ //
+ // "R2" --> error in transverse radius
+ // "R3" --> error in spherical radius
+ // "PHI" --> error in azimuthal angle
+ // "THETA" --> error in polar angle
+ // "SQ" --> get the square of error
+ //
+ // In order to get the error on the cartesian coordinates
+ // reference to the inline ErrX2(), ErrY2() adn ErrZ2() methods.
+
+ TString opt(option);
+ Double_t errorSq = 0.0;
+ opt.ToUpper();
+
+ if (opt.Contains("R2")) {
+ errorSq = fX*fX*fEX2 + fY*fY*fEY2;
+ errorSq /= GetR2sq();
+ }
+ else if (opt.Contains("R3")) {
+ errorSq = fX*fX*fEX2 + fY*fY*fEY2 + fZ*fZ*fEZ2;
+ errorSq /= GetR3sq();
+ }
+ else if (opt.Contains("PHI")) {
+ errorSq = fY*fY*fEX2;
+ errorSq += fX*fX*fEY2;
+ errorSq /= GetR2sq() * GetR2sq();
+ }
+ else if (opt.Contains("THETA")) {
+ errorSq = fZ*fZ * (fX*fX*fEX2 + fY*fY*fEY2);
+ errorSq += GetR2sq() * GetR2sq() * fEZ2;
+ errorSq /= GetR3sq() * GetR3sq() * GetR2() * GetR2();
+ }
+
+ if (!opt.Contains("SQ"))
+ return TMath::Sqrt(errorSq);
+ else
+ return errorSq;
+}
+
+
+
+/********************************************
+ ********************************************
+ *** AliITStrackerANN::AliITSneuron class ***
+ ********************************************
+ ********************************************/
+
+//__________________________________________________________________________________
+AliITStrackerANN::AliITSneuron::AliITSneuron
+(AliITSnode *inner, AliITSnode *outer, Double_t minAct, Double_t maxAct)
+ : fUsed(0), fInner(inner), fOuter(outer)
+{
+ // Default neuron constructor
+ fActivation = gRandom->Rndm() * (maxAct-minAct) + minAct;
+ fGain = new TObjArray;
+}
+
+//__________________________________________________________________________________
+inline Double_t AliITStrackerANN::AliITSneuron::Activate(Double_t temperature)
+{
+ // This computes the new activation of a neuron, and returns
+ // its activation variation as a consequence of the updating.
+ //
+ // Arguments:
+ // - the 'temperature' parameter for the neural activation logistic function
+ //
+ // Operations:
+ // - collects the total gain, by summing the products
+ // of the activation of each sequenced unit by the relative weight.
+ // - collects the total cost, by summing the activations of
+ // all competing units
+ // - passes the sum of gain - cost to the activation function and
+ // calculates the new activation
+ //
+ // Return value:
+ // - the difference between the old activation and the new one
+ // (absolute value)
+
+ // local variables
+ Double_t sumGain = 0.0; // total contribution from chained neurons
+ Double_t sumCost = 0.0; // total contribution from crossing neurons
+ Double_t input; // total input
+ Double_t actOld, actNew; // old and new values for the activation
+ AliITSneuron *linked = 0; // cursor for scanning the neuron arrays (for link check)
+ AliITSlink *link; // cursor for scanning the synapses arrays (for link check)
+ TObjArrayIter *iterator = 0; // pointer to the iterator along the neuron arrays
+
+ // sum contributions from the correlated units
+ iterator = (TObjArrayIter*)fGain->MakeIterator();
+ for(;;) {
+ link = (AliITSlink*)iterator->Next();
+ if (!link) break;
+ sumGain += link->Contribution();
+ }
+ delete iterator;
+
+ // sum contributions from the competing units:
+ // the ones which have the same starting point...
+ iterator = (TObjArrayIter*)fInner->InnerOf()->MakeIterator();
+ for (;;) {
+ linked = (AliITSneuron*)iterator->Next();
+ if (!linked) break;
+ if (linked == this) continue;
+ sumCost += linked->fActivation;
+ }
+ delete iterator;
+ // ...and the ones which have the same ending point
+ iterator = (TObjArrayIter*)fOuter->OuterOf()->MakeIterator();
+ for (;;) {
+ linked = (AliITSneuron*)iterator->Next();
+ if (!linked) break;
+ if (linked == this) continue;
+ sumCost += linked->fActivation;
+ }
+
+ // calculate the total input as the difference between gain and cost
+ input = (sumGain - sumCost) / temperature;
+ actOld = fActivation;
+ // calculate the final output
+#ifdef NEURAL_LINEAR
+ if (input <= -2.0 * temperature)
+ actNew = 0.0;
+ else if (input >= 2.0 * temperature)
+ actNew = 1.0;
+ else
+ actNew = input / (4.0 * temperature) + 0.5;
+#else
+ actNew = 1.0 / (1.0 + TMath::Exp(-input));
+#endif
+ fActivation = actNew;
+
+ // return the activation variation
+ return TMath::Abs(actNew - actOld);
+}
+
+
+
+/******************************************
+ ******************************************
+ *** AliITStrackerANN::AliITSlink class ***
+ ******************************************
+ ******************************************/
+
+ // No methods defined non-inline
+
+
+
+ /**********************************************
+ **********************************************
+ *** AliITStrackerANN::AliITStrackANN class ***
+ **********************************************
+ **********************************************/
+
+//__________________________________________________________________________________
+AliITStrackerANN::AliITStrackANN::AliITStrackANN(Int_t dim) : fNPoints(dim)
+{
+ // Default constructor for the AliITStrackANN class
+
+ fXCenter = 0.0;
+ fYCenter = 0.0;
+ fRadius = 0.0;
+ fCurv = 0.0;
+ fDTrans = 0.0;
+ fDLong = 0.0;
+ fTanLambda = 0.0;
+
+ if (! dim) {
+ fNode = 0;
+ }
+ else{
+ Int_t i = 0;
+ fNode = new AliITSnode*[dim];
+ for (i = 0; i < dim; i++) fNode[i] = 0;
+ }
+}
+
+//__________________________________________________________________________________
+Int_t AliITStrackerANN::AliITStrackANN::CheckOccupation() const
+{
+ // Returns the number of pointers fNode which are not NULL
+
+ Int_t i; // cursor
+ Int_t count = 0; // counter for how many points are stored in the track
+
+ for (i = 0; i < fNPoints; i++) {
+ if (fNode[i] != NULL) count++;
+ }
+
+ return count;
+}
+
+//__________________________________________________________________________________
+Bool_t AliITStrackerANN::AliITStrackANN::RiemannFit()
+{
+ // Performs the Riemann Sphere fit for the given points to a circle
+ // and then uses the fit parameters to fit a helix in space.
+ //
+ // Return values:
+ // - kTRUE if all operations have been performed
+ // - kFALSE if the numbers risk to lead to an arithmetic violation
+
+ Int_t i, j, count, dim = fNPoints;
+
+ // First check for all points
+ count = CheckOccupation();
+ if (count != fNPoints) {
+ Error ("AliITStrackANN::RiemannFit", "CheckOccupations returns %d, fNPoints = %d ==> MISMATCH", count, fNPoints);
+ return kFALSE;
+ }
+
+ // matrix of ones
+ TMatrixD m1(dim,1);
+ for (i = 0; i < dim; i++) m1(i,0) = 1.0;
+
+ // matrix of Rieman projection coordinates
+ TMatrixD coords(dim,3);
+ for (i = 0; i < dim; i++) {
+ coords(i,0) = fNode[i]->X();
+ coords(i,1) = fNode[i]->Y();
+ coords(i,2) = fNode[i]->GetR2sq();
+ }
+
+ // matrix of weights
+ Double_t xterm, yterm, ex, ey;
+ TMatrixD weights(dim,dim);
+ for (i = 0; i < dim; i++) {
+ xterm = fNode[i]->X() * fNode[i]->GetPhi() - fNode[i]->Y() / fNode[i]->GetR2();
+ ex = fNode[i]->ErrX2();
+ yterm = fNode[i]->Y() * fNode[i]->GetPhi() + fNode[i]->X() / fNode[i]->GetR2();
+ ey = fNode[i]->ErrY2();
+ weights(i,i) = fNode[i]->GetR2sq() / (xterm * xterm * ex + yterm * yterm * ey );
+ }
+
+ // weighted sample mean
+ Double_t meanX = 0.0, meanY = 0.0, meanW = 0.0, sw = 0.0;
+ for (i = 0; i < dim; i++) {
+ meanX += weights(i,i) * coords(i,0);
+ meanY += weights(i,i) * coords(i,1);
+ meanW += weights(i,i) * coords(i,2);
+ sw += weights(i,i);
+ }
+ meanX /= sw;
+ meanY /= sw;
+ meanW /= sw;
+
+ // sample covariance matrix
+ for (i = 0; i < dim; i++) {
+ coords(i,0) -= meanX;
+ coords(i,1) -= meanY;
+ coords(i,2) -= meanW;
+ }
+ TMatrixD coordsT(TMatrixD::kTransposed, coords);
+ TMatrixD weights4coords(weights, TMatrixD::kMult, coords);
+ TMatrixD sampleCov(coordsT, TMatrixD::kMult, weights4coords);
+ for (i = 0; i < 3; i++) {
+ for (j = i + 1; j < 3; j++) {
+ sampleCov(i,j) = sampleCov(j,i) = (sampleCov(i,j) + sampleCov(j,i)) * 0.5;
+ }
+ }
+
+ // Eigenvalue problem solving for V matrix
+ Int_t ileast = 0;
+ TVectorD eval(3), n(3);
+ TMatrixD evec = sampleCov.EigenVectors(eval);
+ if (eval(1) < eval(ileast)) ileast = 1;
+ if (eval(2) < eval(ileast)) ileast = 2;
+ n(0) = evec(0, ileast);
+ n(1) = evec(1, ileast);
+ n(2) = evec(2, ileast);
+
+ // c - known term in the plane intersection with Riemann axes
+ Double_t c = -(meanX * n(0) + meanY * n(1) + meanW * n(2));
+
+ // center and radius of fitted circle
+ Double_t xc, yc, radius, curv;
+ xc = -n(0) / (2. * n(2));
+ yc = -n(1) / (2. * n(2));
+ radius = (1. - n(2)*n(2) - 4.*c*n(2)) / (4. * n(2) * n(2));
+
+ if (radius <= 0.E0) {
+ Error("RiemannFit", "Radius = %f less than zero!!!", radius);
+ return kFALSE;
+ }
+ radius = TMath::Sqrt(radius);
+ curv = 1.0 / radius;
+
+ // evaluating signs for curvature and others
+ Double_t phi1 = 0.0, phi2, temp1, temp2, phi0, sumdphi = 0.0;
+ AliITSnode *p = fNode[0];
+ phi1 = p->GetPhi();
+ for (i = 1; i < dim; i++) {
+ p = (AliITSnode*)fNode[i];
+ if (!p) break;
+ phi2 = p->GetPhi();
+ temp1 = phi1;
+ temp2 = phi2;
+ if (temp1 > fgkPi && temp2 < fgkPi)
+ temp2 += fgkTwoPi;
+ else if (temp1 < fgkPi && temp2 > fgkPi)
+ temp1 += fgkTwoPi;
+ sumdphi += temp2 - temp1;
+ phi1 = phi2;
+ }
+ if (sumdphi < 0.E0) curv = -curv;
+ Double_t diff, angle = TMath::ATan2(yc, xc);
+ if (curv < 0.E0)
+ phi0 = angle + 0.5 * TMath::Pi();
+ else
+ phi0 = angle - 0.5 * TMath::Pi();
+ diff = angle - phi0;
+
+ Double_t dt, temp = TMath::Sqrt(xc*xc + yc*yc) - radius;
+ if (curv >= 0.E0)
+ dt = temp;
+ else
+ dt = -temp;
+ //cout << "Dt = " << dt << endl;
+
+ Double_t halfC = 0.5 * curv, test;
+ Double_t *s = new Double_t[dim], *zz = new Double_t[dim], *ws = new Double_t[dim];
+ for (j = 0; j < 6; j++) {
+ p = fNode[j];
+ if (!p) break;
+ //----
+ s[j] = (p->GetR2sq() - dt * dt) / (1. + curv * dt);
+ if (s[j] < 0.) {
+ if (fabs(s[j]) < 1.E-6) s[j] = 0.;
+ else {
+ Error("RiemannFit", "Square root argument error: %17.15g < 0", s[j]);
+ return kFALSE;
+ }
+ }
+ s[j] = TMath::Sqrt(s[j]);
+ //cout << "Curv = " << halfC << " --- s[" << j << "] = " << s[j] << endl;
+ s[j] *= halfC;
+ test = TMath::Abs(s[j]);
+ if (test > 1.) {
+ if (test <= 1.1)
+ s[j] = ((s[j] > 0.) ? 0.99999999999 : -0.9999999999);
+ else {
+ Error("RiemannFit", "Value too large: %17.15g", s[j]);
+ return kFALSE;
+ }
+ }
+ //----
+ zz[j] = p->Z();
+ s[j] = TMath::ASin(s[j]) / halfC;
+ ws[j] = 1.0 / (p->ErrZ2());
+ }
+
+ // second tep final fit
+ Double_t s2Sum = 0.0, zSum = 0.0, szSum = 0.0, sSum = 0.0, sumw = 0.0;
+ for (i = 0; i < dim; i++) {
+ s2Sum += ws[i] * s[i] * s[i];
+ zSum += ws[i] * zz[i];
+ sSum += ws[i] * s[i];
+ szSum += ws[i] * s[i] * zz[i];
+ sumw += ws[i];
+ }
+ s2Sum /= sumw;
+ zSum /= sumw;
+ sSum /= sumw;
+ szSum /= sumw;
+ temp = s2Sum - sSum*sSum;
+
+ Double_t dz, tanL;
+ dz = (s2Sum*zSum - sSum*szSum) / temp;
+ tanL = (szSum - sSum*zSum) / temp;
+
+ fXCenter = xc;
+ fYCenter = yc;
+ fRadius = radius;
+ fCurv = curv;
+ fPhi = phi0;
+ fDTrans = dt;
+ fDLong = dz;
+ fTanLambda = tanL;
+
+ delete [] s;
+ delete [] zz;
+ delete [] ws;
+
+ return kTRUE;
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff