Minor Update for changes in AliITSv11Geometry class.

[u/mrichter/AliRoot.git] / ITS / AliITStrackerANN.cxx

// ---------------------------------------------------------------------------------
// 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 <Riostream.h>
#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"

const Double_t AliITStrackerANN::fgkPi     = 3.141592653; // pi
const Double_t AliITStrackerANN::fgkHalfPi = 1.570796327; // pi / 2
const Double_t AliITStrackerANN::fgkTwoPi  = 6.283185307; // 2 * pi

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 ***
 ******************************************
 ******************************************/
 
//__________________________________________________________________________________ 
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;
}

//__________________________________________________________________________________ 
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();
}
 
//__________________________________________________________________________________ 
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;
}
 
//__________________________________________________________________________________
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 (TMath::Abs(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;
}