New version of Neural Tracking which works with AliITSclusterV2 data and returns...
authorbarbera <barbera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 23 Mar 2004 13:52:20 +0000 (13:52 +0000)
committerbarbera <barbera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 23 Mar 2004 13:52:20 +0000 (13:52 +0000)
ITS/AliITSFindTracksANN.C [new file with mode: 0644]
ITS/AliITStrackerANN.cxx [new file with mode: 0644]
ITS/AliITStrackerANN.h [new file with mode: 0644]
ITS/ITSLinkDef.h
ITS/libITS.pkg

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