Neural ITS standalone tracking revisited
authorbarbera <barbera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 8 Mar 2004 16:35:43 +0000 (16:35 +0000)
committerbarbera <barbera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 8 Mar 2004 16:35:43 +0000 (16:35 +0000)
ITS/AliITSNeuralPoint.cxx [new file with mode: 0644]
ITS/AliITSNeuralPoint.h [new file with mode: 0644]
ITS/AliITSNeuralTrack.cxx [new file with mode: 0644]
ITS/AliITSNeuralTrack.h [new file with mode: 0644]
ITS/AliITSNeuralTracker.cxx [new file with mode: 0644]
ITS/AliITSNeuralTracker.h [new file with mode: 0644]

diff --git a/ITS/AliITSNeuralPoint.cxx b/ITS/AliITSNeuralPoint.cxx
new file mode 100644 (file)
index 0000000..299b868
--- /dev/null
@@ -0,0 +1,226 @@
+#include <stdlib.h>
+#include <Riostream.h>
+
+#include <TString.h>
+
+#include "AliITSRecPoint.h"
+#include "AliITSclusterV2.h"
+#include "AliITSgeom.h"
+#include "AliITSgeomMatrix.h"
+
+#include "AliITSNeuralPoint.h"
+
+
+ClassImp(AliITSNeuralPoint)
+//
+//------------------------------------------------------------------------------------------------------
+//
+AliITSNeuralPoint::AliITSNeuralPoint()
+{
+       // Default constructor.
+       // Defines the point as a noise point in the origin.
+       
+       fX = fY = fZ = 0.;
+       fEX = fEY = fEZ = 0.;
+       fLayer = 0;
+       fLabel[0] = fLabel[1] = fLabel[2] = -1;
+       fModule = 0;
+       fIndex = 0;
+       fUser = 0;
+}
+//
+//------------------------------------------------------------------------------------------------------
+//
+AliITSNeuralPoint::AliITSNeuralPoint(AliITSNeuralPoint *p) :
+fX(p->fX), fY(p->fY), fZ(p->fZ), fEX(p->fEX), fEY(p->fEY), fEZ(p->fEZ)
+{
+       // Modified copy constructor.
+       // Accepts a pointer to a like object and copies its datamembers.
+       
+       fLayer = p->fLayer;
+       for (Int_t i = 0; i < 3; i++) fLabel[i] = p->fLabel[i];
+       fModule = p->fModule;
+       fIndex = p->fIndex;
+       fUser = p->fUser;
+       fCharge = p->fCharge;
+}
+//
+//------------------------------------------------------------------------------------------------------
+//
+AliITSNeuralPoint::AliITSNeuralPoint(AliITSRecPoint *rp, AliITSgeomMatrix *gm)
+{
+       // Conversion constructor.
+       // Accepts a AliITSRecPoint and a AliITSgeomMatrix,
+       // and converts the local coord of the AliITSRecPoint object into global
+       
+       Int_t i, k;
+       Double_t locPos[3], globPos[3], locErr[3][3], globErr[3][3];
+       for (i = 0; i < 3; i++) {
+               globPos[i] = 0.0;
+               for (k = 0; k < 3; k++) {
+                       locErr[i][k] = 0.0;
+                       globErr[i][k] = 0.0;
+               }
+       }
+       
+       // local to global conversions of coords
+       locPos[0] = rp->fX;
+       locPos[1] = 0.0;
+       locPos[2] = rp->fZ;
+       gm->LtoGPosition(locPos, globPos);
+       fX = globPos[0];
+       fY = globPos[1];
+       fZ = globPos[2];
+
+       // local to global conversions of sigmas
+       locErr[0][0] = rp->fSigmaX2;
+       locErr[2][2] = rp->fSigmaZ2;
+       gm->LtoGPositionError(locErr, globErr);
+       for (i = 0; i < 3; i++) fLabel[i] = rp->fTracks[i];
+       fEX = TMath::Sqrt(globErr[0][0]);
+       fEY = TMath::Sqrt(globErr[1][1]);
+       fEZ = TMath::Sqrt(globErr[2][2]);
+
+       // copy of other data-members
+       fCharge = rp->fQ;
+       fLayer = 0;
+       fIndex = 0;
+       fModule = 0;
+       fUser = 0;
+}
+//
+//-------------------------------------------------------------------------------------------------
+//
+AliITSNeuralPoint::AliITSNeuralPoint
+(AliITSclusterV2 *rp, AliITSgeom *geom, Short_t module, Short_t index)
+{
+       Int_t mod = (Int_t)module, lay, lad, det;
+       fModule = module;
+       fIndex = index;
+       geom->GetModuleId(mod, lay, lad, det);
+       fLayer = (Short_t)lay;
+       
+       Double_t rot[9];  
+       Float_t tx, ty, tz;
+       geom->GetRotMatrix(fModule, rot);
+       geom->GetTrans(fLayer, lad, det, tx, ty, tz);
+       
+       Double_t r, phi, cosPhi, sinPhi;
+       r = -tx*rot[1] + ty*rot[0];
+       if (lay == 1) r = -r;
+       phi = TMath::ATan2(rot[1], rot[0]);
+       if (lay==1) phi -= 3.1415927;
+       cosPhi = TMath::Cos(phi);
+       sinPhi = TMath::Sin(phi);
+       fX =  r*cosPhi + rp->GetY()*sinPhi;
+       fY = -r*sinPhi + rp->GetY()*cosPhi;
+       fZ = rp->GetZ();
+       fEX = TMath::Sqrt(rp->GetSigmaY2())*sinPhi;
+       fEY = TMath::Sqrt(rp->GetSigmaY2())*cosPhi;
+       fEZ = TMath::Sqrt(rp->GetSigmaZ2());
+       fLayer--;
+}
+//
+//-------------------------------------------------------------------------------------------------
+//
+Double_t AliITSNeuralPoint::GetPhi() const
+// Returns the azimuthal coordinate in the range 0-2pi
+{
+       Double_t q;
+       q = TMath::ATan2(fY,fX); 
+       if (q >= 0.) 
+               return q;
+       else 
+               return q + 2.0*TMath::Pi();
+}
+//
+//------------------------------------------------------------------------------------------------------
+//
+Double_t AliITSNeuralPoint::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 ErrX(), ErrY() adn ErrZ() methods.
+{
+       TString opt(option);
+       Double_t errorSq = 0.0;
+       opt.ToUpper();
+
+       if (opt.Contains("R2")) {
+               errorSq  = fX*fX*fEX*fEX + fY*fY*fEY*fEY;
+               errorSq /= GetR2sq();
+       }
+       else if (opt.Contains("R3")) {
+               errorSq  = fX*fX*fEX*fEX + fY*fY*fEY*fEY + fZ*fZ*fEZ*fEZ;
+               errorSq /= GetR3sq();
+       }
+       else if (opt.Contains("PHI")) {
+               errorSq  = fY*fY*fEX*fEX;
+               errorSq += fX*fX*fEY*fEY;
+               errorSq /= GetR2sq() * GetR2sq();
+       }
+       else if (opt.Contains("THETA")) {
+               errorSq = fZ*fZ * (fX*fX*fEX*fEX + fY*fY*fEY*fEY);
+               errorSq += GetR2sq() * GetR2sq() * fEZ*fEZ;
+               errorSq /= GetR3sq() * GetR3sq() * GetR2() * GetR2();
+       }
+       
+       if (!opt.Contains("SQ")) 
+               return TMath::Sqrt(errorSq);
+       else 
+               return errorSq;
+}
+//
+//------------------------------------------------------------------------------------------------------
+//
+Bool_t AliITSNeuralPoint::HasID(Int_t ID)
+// Checks if the recpoint belongs to the GEANT track
+// whose label is specified in the argument
+{
+       if (ID<0) 
+               return kFALSE; 
+       else 
+               return (fLabel[0]==ID || fLabel[1]==ID || fLabel[2]==ID);
+}
+//
+//------------------------------------------------------------------------------------------------------
+//
+Int_t* AliITSNeuralPoint::SharedID(AliITSNeuralPoint *p)
+// Checks if there is a GEANT track owning both
+// <this> and the recpoint in the argument
+// The return value is an array of 4 integers.
+// The firs integer returns the count of matches between labels of 
+// <this> and labels of the argument (0 to 3)
+// The other three return the matched labels.
+// If a NULL pointer is passed, the array will be returned as:
+// {0, -1, -1, -1}
+{
+       Int_t i, *shared = new Int_t[4];
+       for (i = 0; i < 4; i++) shared[i] = -1;
+       shared[0] = 0;
+       if (!p) return shared;
+       for (i = 0; i < 3; i++) {
+               if (HasID(p->fLabel[i])) shared[i + 1] = p->fLabel[i];
+               shared[0]++;
+       }
+       return shared;
+}
+//
+//
+//
+void AliITSNeuralPoint::ConfMap(Double_t vx, Double_t vy)
+{
+       Double_t dx = fX - vx;
+       Double_t dy = vy - fY;
+       Double_t r2 = dx*dx + dy*dy;
+       fConfX = dx / r2;
+       fConfY = dy / r2;
+}
diff --git a/ITS/AliITSNeuralPoint.h b/ITS/AliITSNeuralPoint.h
new file mode 100644 (file)
index 0000000..631746d
--- /dev/null
@@ -0,0 +1,85 @@
+#ifndef ALIITSNEURALPOINT_H
+#define ALIITSNEURALPOINT_H
+
+#include <TMath.h>
+
+class AliITSgeom;
+class AliITSgeomMatrix;
+class AliITSRecPoint;
+class AliITSclusterV2;
+
+class AliITSNeuralPoint : public TObject {
+
+public:
+
+       AliITSNeuralPoint();
+       AliITSNeuralPoint(AliITSNeuralPoint *p);
+       AliITSNeuralPoint(AliITSRecPoint *rp, AliITSgeomMatrix *gm);
+       AliITSNeuralPoint(AliITSclusterV2 *rp, AliITSgeom *geom, Short_t module, Short_t index);
+
+       virtual ~ AliITSNeuralPoint() { }
+
+       Double_t& X()    {return fX;}   // reference to X coord
+       Double_t& Y()    {return fY;}   // reference to Y coord
+       Double_t& Z()    {return fZ;}   // reference to Z coord
+       Double_t& ErrX() {return fEX;}  // reference to X error
+       Double_t& ErrY() {return fEY;}  // reference to Y error
+       Double_t& ErrZ() {return fEZ;}  // reference to Z error
+               
+       Double_t  GetR2()    const  {return TMath::Sqrt(GetR2sq());} // xy radius
+       Double_t  GetR3()    const  {return TMath::Sqrt(GetR3sq());} // 3D radius
+       Double_t  GetR2sq()  const  {return fX*fX+fY*fY;}            // xy rad. square
+       Double_t  GetR3sq()  const  {return GetR2sq()+fZ*fZ;}        // 3D rad. square
+       Double_t  GetPhi()   const;
+       Double_t  GetTheta() const  {return TMath::ATan2(GetR2(),fZ);} // polar angle
+       Double_t  GetConfX() const  {return fConfX;}
+       Double_t  GetConfY() const  {return fConfY;}
+       Double_t  GetError(Option_t *opt);
+       void      ConfMap(Double_t vx, Double_t vy);
+
+       Double_t  GetCharge()       const {return fCharge;}        // ADC signal
+       Short_t   GetIndex()        const {return fIndex;}         // Reference in TreeR
+       Long_t    GetLabel(Int_t i) const {return fLabel[Chk(i)];} // GEANT owner particle
+       Short_t   GetLayer()        const {return fLayer;}         // ITS layer
+       Short_t   GetModule()       const {return fModule;}        // ITS module 
+       Short_t   GetUser()         const {return fUser;}          // Found track owner
+                               
+       void      SetCharge(Double_t val)       {fCharge = val;}
+       void      SetIndex(Short_t val)         {fIndex = val;}
+       void      SetLabel(Int_t i, Long_t val) {fLabel[Chk(i)] = val;}
+       void      SetLayer(Short_t val)         {fLayer = val;}
+       void      SetModule(Short_t val)        {fModule = val;}
+       void      SetUser(Short_t val)          {fUser = val;}
+       
+       Bool_t    HasID (Int_t ID);
+       Int_t*    SharedID(AliITSNeuralPoint *p);
+
+protected:
+       
+       Int_t     Chk(Int_t i) const {if(i<0)i=0;if(i>=3)i=3;return i;}
+
+       Double_t  fX;   //
+       Double_t  fY;   // position
+       Double_t  fZ;   //
+       
+       Double_t  fConfX; // conformal mapping X
+       Double_t  fConfY; // conformal mapping Y
+               
+       Double_t  fEX;  //
+       Double_t  fEY;  // position error
+       Double_t  fEZ;  //
+
+       Double_t  fCharge;   // total charge signal in cluster
+
+       Short_t   fModule;   // ITS module containing the point (0 - 2197)
+       Short_t   fIndex;    // index as TClonesArray entry in TreeR (usually not > 600)
+       Short_t   fLayer;    // ITS layer containing the point
+       Short_t   fUser;     // owner recognized track or flag to evidence using
+       Short_t   fZSort;    // order as a function of local Z
+
+       Int_t     fLabel[3]; // GEANT labels of the owner tracks
+
+       ClassDef(AliITSNeuralPoint, 1) // AliITSNeuralPoints class
+};
+
+#endif
diff --git a/ITS/AliITSNeuralTrack.cxx b/ITS/AliITSNeuralTrack.cxx
new file mode 100644 (file)
index 0000000..8b67365
--- /dev/null
@@ -0,0 +1,1080 @@
+#include <Riostream.h>
+#include <cstdlib>
+#include <cstring>
+
+#include <TObject.h>
+#include <TROOT.h>
+#include <TMath.h>
+#include <TString.h>
+#include <TObjArray.h>
+#include <TH1.h>
+#include <TMatrixD.h>
+
+//#include "AliITSVertex.h"
+#include "AliITSIOTrack.h"
+#include "AliITSNeuralPoint.h"
+
+#include "AliITSNeuralTrack.h"
+
+
+
+ClassImp(AliITSNeuralTrack)
+//
+//
+//
+AliITSNeuralTrack::AliITSNeuralTrack() : fMatrix(5,5), fVertex()
+{
+       Int_t i;
+
+       fMass  = 0.1396;   // default assumption: pion
+       fField = 2.0;      // default assumption: B = 0.4 Tesla
+
+       fXC = fYC = fR = fC = 0.0;
+       fTanL = fG0 = fDt = fDz = 0.0;
+       fStateR = fStatePhi = fStateZ = fChi2 = fNSteps = 0.0;
+
+       fLabel = 0;
+       fCount = 0;
+       for (i = 0; i < 6; i++) fPoint[i] = 0;
+
+       fVertex.X() = 0.0;
+       fVertex.Y() = 0.0;
+       fVertex.Z() = 0.0;
+       fVertex.ErrX() = 0.0;
+       fVertex.ErrY() = 0.0;
+       fVertex.ErrZ() = 0.0;
+}
+//
+//
+//
+AliITSNeuralTrack::~AliITSNeuralTrack()
+{
+       Int_t l;
+       for (l = 0; l < 6; l++) fPoint[l] = 0;
+}
+//
+//
+//
+void AliITSNeuralTrack::AssignLabel()
+// Assigns a GEANT label to the found track. 
+// Every cluster has up to three labels (it can have less). Then each label is 
+// recorded for each point. Then, counts are made to check if some of the labels
+// appear more than once. Finally, the label which appears most times is assigned
+// to the track in the field fLabel.
+// The number of points containing that label is counted in the fCount data-member. 
+{      
+       Bool_t found;
+       Int_t i, j, l, lab, max = 0;
+       
+       // We have up to 6 points for 3 labels each => up to 18 possible different values
+       Int_t idx[18] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
+       Int_t count[18] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
+       
+       for (l = 0; l < 6; l++) {
+               if (!fPoint[l]) continue;
+               // Sometimes the same label appears two times in the same recpoint.
+               // With these if statements, such problem is solved by turning
+               // one of them to -1.
+               if (fPoint[l]->GetLabel(1) >= 0 && fPoint[l]->GetLabel(1) == fPoint[l]->GetLabel(0))
+                       fPoint[l]->SetLabel(1, -1);
+               if (fPoint[l]->GetLabel(2) >= 0 && fPoint[l]->GetLabel(2) == fPoint[l]->GetLabel(0))
+                       fPoint[l]->SetLabel(2, -1);
+               if (fPoint[l]->GetLabel(2) >= 0 && fPoint[l]->GetLabel(2) == fPoint[l]->GetLabel(1))
+                       fPoint[l]->SetLabel(2, -1);
+               for (i = 0; i < 3; i++) {
+                       lab = fPoint[l]->GetLabel(i);
+                       if (lab < 0) continue;
+                       found = kFALSE;
+                       for (j = 0; j < max; j++) {
+                               if (idx[j] == lab) {
+                                       count[j]++;
+                                       found = kTRUE;
+                               }
+                       }
+                       if(!found) {
+                               max++;
+                               idx[max - 1] = lab;
+                               count[max - 1] = 1;
+                       }
+               }
+       }
+       
+       j = 0, max = count[0];
+       for (i = 0; i < 18; i++) {
+               if (count[i] > max) {
+                       j = i;
+                       max = count[i];
+               }
+       }
+       fLabel = idx[j];
+       fCount = count[j];
+}
+//
+//
+//
+void AliITSNeuralTrack::CleanSlot(Int_t i, Bool_t del)
+// Removes a point from the corresponding layer slot in the found track.
+// If the argument is TRUE, the point object is also deleted from heap.
+{
+       if (i >= 0 && i < 6) {
+               if (del) delete fPoint[i];
+               fPoint[i] = 0;
+       }
+}
+//
+//
+//
+void AliITSNeuralTrack::GetModuleData(Int_t layer, Int_t &mod, Int_t &pos)
+// Returns the point coordinates according to the TreeR philosophy in galice.root files
+// that consist in the module number (mod) and the position in the TClonesArray of
+// the points reconstructed in that module for the run being examined.
+{
+       if (layer < 0 || layer > 5) {
+               Error("GetModuleData", "Layer out of range: %d", layer);
+               return;
+       }
+       mod = fPoint[layer]->GetModule();
+       pos = fPoint[layer]->GetIndex();
+}
+//
+//
+//
+void AliITSNeuralTrack::Insert(AliITSNeuralPoint *point)
+// A trivial method to insert a point in the tracks;
+// the point is inserted to the slot corresponding to its ITS layer.
+{
+       if (!point) return;
+       
+       Int_t layer = point->GetLayer();
+       if (layer < 0 || layer > 6) {
+               Error("Insert", "Layer index %d out of range", layer);
+               return;
+       } 
+       
+       fPoint[layer] = point;
+}
+//
+//
+//
+Int_t AliITSNeuralTrack::OccupationMask()
+// Returns a byte which maps the occupied slots. 
+// Each bit represents a layer going from the less significant on.
+{
+       Int_t i, check, mask = 0;
+       for (i = 0; i < 6; i++) {
+               check = 1 << i;
+               if (fPoint[i]) mask |= check;
+       }
+       return mask;
+}
+//
+//
+//
+void AliITSNeuralTrack::PrintLabels()
+// Prints the results of the AssignLabel() method, together with
+// the GEANT labels assigned to each point, in order to evaluate 
+// how the assigned label is distributed among points.
+{
+       cout << "Assigned label = " << fLabel << " -- counted " << fCount << " times: " << endl << endl;
+       for (Int_t i = 0; i < 6; i++) {
+               cout << "Point #" << i + 1 << " --> ";
+               if (fPoint[i]) {
+                       cout << "labels = " << fPoint[i]->GetLabel(0) << ", ";
+                       cout << fPoint[i]->GetLabel(1) << ", ";
+                       cout << fPoint[i]->GetLabel(2) << endl;
+               }
+               else {
+                       cout << "not assigned" << endl;
+               }
+       }
+       cout << endl;
+}
+//
+//
+//
+Bool_t AliITSNeuralTrack::AddEL(Int_t layer, Double_t sign)
+{
+       Double_t width = 0.0;
+       switch (layer) {
+               case 0: width = 0.00260 + 0.00283; break;
+               case 1: width = 0.0180; break;
+               case 2: width = 0.0094; break;
+               case 3: width = 0.0095; break;
+               case 4: width = 0.0091; break;
+               case 5: width = 0.0087; break;
+               default:
+                       Error("AddEL", "Layer value %d out of range!", layer);
+                       return kFALSE;
+       }
+       width *= 1.7;
+
+       if((layer == 5) && (fStatePhi < 0.174 || fStatePhi > 6.100 || (fStatePhi > 2.960 && fStatePhi < 3.31))) {
+               width += 0.012;
+       }
+
+       Double_t invSqCosL = 1. + fTanL * fTanL;            // = 1 / (cos(lambda)^2) = 1 + tan(lambda)^2
+       Double_t invCosL   = TMath::Sqrt(invSqCosL);        // = 1  / cos(lambda)
+       Double_t pt = GetPt();                              // = transverse momentum
+       Double_t p2 = pt *pt * invSqCosL;                   // = square modulus of momentum
+       Double_t energy = TMath::Sqrt(p2 + fMass * fMass);  // = energy
+       Double_t beta2 = p2 / (p2 + fMass * fMass);         // = (v / c) ^ 2
+       if (beta2 == 0.0) {
+               printf("Anomaly in AddEL: pt=%8.6f invSqCosL=%8.6f fMass=%8.7f --> beta2 = %8.7f\n", pt, invSqCosL, fMass, beta2);
+               return kFALSE;
+       }
+
+       Double_t dE = 0.153 / beta2 * (log(5940. * beta2 / (1. - beta2)) - beta2) * width * 21.82 * invCosL;
+       dE = sign * dE * 0.001;
+
+       energy += dE;
+       p2 = energy * energy - fMass * fMass;
+       pt = TMath::Sqrt(p2) / invCosL;
+       if (fC < 0.) pt = -pt;
+       fC = (0.299792458 * 0.2 * fField) / (pt * 100.);
+
+       return kTRUE;
+}
+//
+//
+//
+Bool_t AliITSNeuralTrack::AddMS(Int_t layer)
+{
+       Double_t width = 0.0;
+       switch (layer) {
+               case 0: width = 0.00260 + 0.00283; break;
+               case 1: width = 0.0180; break;
+               case 2: width = 0.0094; break;
+               case 3: width = 0.0095; break;
+               case 4: width = 0.0091; break;
+               case 5: width = 0.0087; break;
+               default:
+                       Error("AddEL", "Layer value %d out of range!", layer);
+                       return kFALSE;
+       }
+       width *= 1.7;
+
+       if((layer == 5) && (fStatePhi < 0.174 || fStatePhi > 6.100 || (fStatePhi > 2.960 && fStatePhi < 3.31))) {
+               width += 0.012;
+       }
+
+       Double_t cosL = TMath::Cos(TMath::ATan(fTanL));
+       Double_t halfC = fC / 2.;
+       Double_t q20 = 1. / (cosL * cosL);
+       Double_t q30 = fC * fTanL;
+
+       Double_t q40 = halfC * (fStateR * fStateR - fDt * fDt) / (1. + 2. * halfC * fDt);
+       Double_t dd  = fDt + halfC * fDt * fDt - halfC * fStateR * fStateR;
+       Double_t dprova = fStateR * fStateR - dd * dd;
+       Double_t q41 = 0.;
+       if(dprova > 0.) q41 = -1. / cosL * TMath::Sqrt(dprova) / (1. + 2. * halfC *fDt);
+
+       Double_t p2 = (GetPt()*GetPt()) / (cosL * cosL);
+       Double_t beta2 = p2 / (p2 + fMass * fMass);
+       Double_t theta2 = 14.1 * 14.1 / (beta2 * p2 * 1.e6) * (width / TMath::Abs(cosL));
+
+       fMatrix(2,2) += theta2 * (q40 * q40 + q41 * q41);
+       fMatrix(3,2) += theta2 * q20 * q40;
+       fMatrix(2,3) += theta2 * q20 * q40;
+       fMatrix(3,3) += theta2 * q20 * q20;
+       fMatrix(4,2) += theta2 * q30 * q40;
+       fMatrix(2,4) += theta2 * q30 * q40;
+       fMatrix(4,3) += theta2 * q30 * q20;
+       fMatrix(3,4) += theta2 * q30 * q20;
+       fMatrix(4,4) += theta2 * q30 * q30;
+       
+       return kTRUE;
+}
+//
+//
+//
+Int_t AliITSNeuralTrack::PropagateTo(Double_t rk)
+{
+       // Propagation method.
+       // Changes the state vector according to a new radial position
+       // which is specified by the passed 'r' value (in cylindircal coordinates).
+       // The covariance matrix is also propagated (and enlarged) according to
+       // the FCFt technique, where F is the jacobian of the new parameters
+       // w.r.t. their old values.
+       // The option argument forces the method to add also the energy loss
+       // and the multiple scattering effects, which respectively have the effect
+       // of changing the curvature and widening the covariance matrix.
+
+       if (rk < fabs(fDt)) {
+               Error("PropagateTo", Form("Impossible propagation to r (=%17.15g) < Dt (=%17.15g)", rk, fDt));
+               return 0;
+       }
+       
+       Double_t duepi = 2. * TMath::Pi();
+       Double_t rkm1 = fStateR;
+       Double_t aAk = ArgPhi(rk), aAkm1 = ArgPhi(rkm1);
+       Double_t ak = ArgZ(rk), akm1 = ArgZ(rkm1);
+
+       fStatePhi += TMath::ASin(aAk) - TMath::ASin(aAkm1);
+       if(fStatePhi > duepi) fStatePhi -= duepi;
+       if(fStatePhi < 0.) fStatePhi += duepi;
+
+       Double_t halfC = 0.5 * fC;
+       fStateZ += fTanL / halfC * (TMath::ASin(ak)-TMath::ASin(akm1));
+       
+       Double_t bk = ArgB(rk), bkm1 = ArgB(rkm1);
+       Double_t ck = ArgC(rk), ckm1 = ArgC(rkm1);
+       
+       Double_t f02 = ck / TMath::Sqrt(1. - aAk * aAk) - ckm1 / TMath::Sqrt(1. - aAkm1 * aAkm1);
+       Double_t f04 = bk / TMath::Sqrt(1. - aAk * aAk) - bkm1 / TMath::Sqrt(1. - aAkm1 * aAkm1);
+       Double_t f12 = fTanL * fDt * (1. / rk - 1. / rkm1);
+       Double_t f13 = rk - rkm1;
+       
+       Double_t c00 = fMatrix(0,0);
+       Double_t c10 = fMatrix(1,0);
+       Double_t c11 = fMatrix(1,1);
+       Double_t c20 = fMatrix(2,0);
+       Double_t c21 = fMatrix(2,1);
+       Double_t c22 = fMatrix(2,2);
+       Double_t c30 = fMatrix(3,0);
+       Double_t c31 = fMatrix(3,1);
+       Double_t c32 = fMatrix(3,2);
+       Double_t c33 = fMatrix(3,3);
+       Double_t c40 = fMatrix(4,0);
+       Double_t c41 = fMatrix(4,1);
+       Double_t c42 = fMatrix(4,2);
+       Double_t c43 = fMatrix(4,3);
+       Double_t c44 = fMatrix(4,4);
+
+       Double_t r10 = c10 + c21*f02 + c41*f04;
+       Double_t r20 = c20 + c22*f02 + c42*f04;
+       Double_t r30 = c30 + c32*f02 + c43*f04;
+       Double_t r40 = c40 + c42*f02 + c44*f04;
+       Double_t r21 = c21 + c22*f12 + c32*f13;
+       Double_t r31 = c31 + c32*f12 + c33*f13;
+       Double_t r41 = c41 + c42*f12 + c43*f13;
+
+       fMatrix(0,0) = c00 + c20*f02 + c40*f04 + f02*r20 + f04*r40;
+       fMatrix(1,0) = fMatrix(0,1) = r10 + f12*r20 + f13*r30;
+       fMatrix(1,1) = c11 + c21*f12 + c31*f13 + f12*r21 + f13*r31;
+       fMatrix(2,0) = fMatrix(0,2) = r20;
+       fMatrix(2,1) = fMatrix(1,2) = r21;
+       fMatrix(3,0) = fMatrix(0,3) = r30;
+       fMatrix(3,1) = fMatrix(1,3) = r31;
+       fMatrix(4,0) = fMatrix(0,4) = r40;
+       fMatrix(4,1) = fMatrix(1,4) = r41;
+
+       fStateR = rk;
+
+       if (rkm1 < fStateR)  // going to greater R --> energy LOSS
+               return -1;
+       else                 // going to smaller R --> energy GAIN
+               return 1;
+}
+//
+//
+//
+Bool_t AliITSNeuralTrack::SeedCovariance()
+{
+       // generate a covariance matrix depending on the results obtained from
+       // the preliminary seeding fit procedure.
+       // It calculates the variances for C, D ans TanL, according to the
+       // differences of the fitted values from the requested ones necessary
+       // to make the curve exactly pass through each point.
+
+       /*
+       Int_t i, j;
+       AliITSNeuralPoint *p = 0;
+       Double_t r, argPhi, phiC, phiD, argZ, zL;
+       Double_t sumC = 0.0, sumD = 0.0, sumphi = 0., sumz = 0., sumL = 0.;
+       for (i = 0; i < fNum; i++) {
+               p = At(i);
+               if (!p) continue;
+               r = p->GetR2();
+               // weight and derivatives of phi and zeta w.r.t. various params
+               sumphi += 1./ p->ErrorGetPhi();
+               argPhi = ArgPhi(r);
+               argZ = ArgZ(r);
+               if (argPhi > 100.0 || argZ > 100.0) {
+                       Error("InitCovariance", "Argument error");
+                       return kFALSE;
+               }
+               phiC = DerArgPhiC(r) / TMath::Sqrt(1.0 - argPhi * argPhi);
+               phiD = DerArgPhiD(r) / TMath::Sqrt(1.0 - argPhi * argPhi);
+               if (phiC > 100.0 || phiD > 100.0) {
+                       Error("InitCovariance", "Argument error");
+                       return kFALSE;
+               }
+               zL = asin(argZ) / fC;
+               sumL += zL * zL;
+               sumC += phiC * phiC;
+               sumD += phiD * phiD;
+               sumz += 1.0 / (p->fError[2] * p->fError[2]);
+       }
+
+       for (i = 0; i < 5; i++) for (j = 0; j < 5; j++) fMatrix(i,j) = 0.;
+       fMatrix(0,0) = 1. / sumphi;
+       fMatrix(1,1) = 1. / sumz;
+       fMatrix(2,2) = 1. / sumD;
+       fMatrix(3,3) = 1. / sumL;
+       fMatrix(4,4) = 1. / sumC;
+       fMatrix.Print();
+       */
+       
+       AliITSNeuralPoint *p = 0;
+       Double_t delta, cs, sn, r, argz;
+       Double_t diffC, diffD, diffL, calcC, calcD, calcL;
+
+       Int_t l;
+       for (l = 0; l < 6; l++) {
+               p = fPoint[l];
+               if (!p) break;
+               sn = TMath::Sin(p->GetPhi() - fG0);
+               cs = TMath::Cos(p->GetPhi() - fG0);
+               r  = p->GetR2();
+               calcC = (fDt/r - sn) / (2.*fDt*sn - r - fDt*fDt/r);
+               argz = ArgZ(r);
+               if (argz > 1000.0) {
+                       Error("Covariance", "Value too high");
+                       return kFALSE;
+               }
+               calcL = (p->Z() - fDz) * fC / asin(argz);
+               delta = fR*fR + r*r + 2.0*fR*r*sin(p->GetPhi() - fG0);
+               if (delta < 0.E0) {
+                       if (delta >= -0.5)
+                               delta = 0.;
+                       else {
+                               Error("Covariance", Form("Discriminant = %g --- Dt = %g", delta, fDt));
+                               return kFALSE;
+                       }
+               }
+               delta = sqrt(delta);
+               if (fC >= 0)
+                       calcD = delta - fR;
+               else
+                       calcD = fR - delta;
+               diffD = calcD - fDt;
+               diffL = calcL - fTanL;
+               diffC = fC - calcC;
+               fMatrix(0,0) += 100000000.0 * p->GetError("phi") * p->GetError("phi");
+               fMatrix(1,1) += 10000.0 * p->ErrZ() * p->ErrZ();
+               fMatrix(2,2) += 100000.0 * diffD * diffD;
+               fMatrix(3,3) += diffL * diffL;
+               fMatrix(4,4) += 100000000.0 * diffC * diffC;
+       }
+       Double_t N = 0.;
+       for (l = 0; l < 6; l++) if (fPoint[l]) N++;
+       fMatrix *= 1./(N++ * N);
+       return kTRUE;
+}
+//
+//
+//
+Bool_t AliITSNeuralTrack::Filter(AliITSNeuralPoint *test)
+{
+       // Makes all calculations which apply the Kalman filter to the
+       // stored guess of the state vector, after propagation to a new layer
+
+       if (!test) {
+               Error("Filter", "Null pointer passed");
+               return kFALSE;
+       }
+       
+       Double_t m[2];
+       Double_t rk, phik, zk;
+       rk = test->GetR2();
+       phik = test->GetPhi();
+       zk = test->Z();
+       m[0]=phik;
+       m[1]=zk;
+       
+       //////////////////////// Evaluation of the error matrix V  /////////
+       Double_t v00 = test->GetError("phi") * rk;
+       Double_t v11 = test->ErrZ();
+       ////////////////////////////////////////////////////////////////////  
+       
+       // Get the covariance matrix
+       Double_t cin00, cin10, cin20, cin30, cin40;
+       Double_t cin11, cin21, cin31, cin41, cin22;
+       Double_t cin32, cin42, cin33, cin43, cin44;
+       cin00 = fMatrix(0,0);
+       cin10 = fMatrix(1,0);
+       cin20 = fMatrix(2,0);
+       cin30 = fMatrix(3,0);
+       cin40 = fMatrix(4,0);
+       cin11 = fMatrix(1,1);
+       cin21 = fMatrix(2,1);
+       cin31 = fMatrix(3,1);
+       cin41 = fMatrix(4,1);
+       cin22 = fMatrix(2,2);
+       cin32 = fMatrix(3,2);
+       cin42 = fMatrix(4,2);
+       cin33 = fMatrix(3,3);
+       cin43 = fMatrix(4,3);
+       cin44 = fMatrix(4,4);
+       
+       // Calculate R matrix
+       Double_t rold00 = cin00 + v00;
+       Double_t rold10 = cin10;
+       Double_t rold11 = cin11 + v11;
+       
+       ////////////////////// R matrix inversion  /////////////////////////
+       Double_t det = rold00*rold11 - rold10*rold10;
+       Double_t r00 = rold11/det;
+       Double_t r10 = -rold10/det;
+       Double_t r11 = rold00/det;
+       ////////////////////////////////////////////////////////////////////
+       
+       // Calculate Kalman matrix
+       Double_t k00 = cin00*r00 + cin10*r10;
+       Double_t k01 = cin00*r10 + cin10*r11;
+       Double_t k10 = cin10*r00 + cin11*r10;  
+       Double_t k11 = cin10*r10 + cin11*r11;
+       Double_t k20 = cin20*r00 + cin21*r10;  
+       Double_t k21 = cin20*r10 + cin21*r11;  
+       Double_t k30 = cin30*r00 + cin31*r10;  
+       Double_t k31 = cin30*r10 + cin31*r11;  
+       Double_t k40 = cin40*r00 + cin41*r10;
+       Double_t k41 = cin40*r10 + cin41*r11;
+       
+       // Get state vector (will keep the old values for phi and z)
+       Double_t x0, x1, x2, x3, x4, savex0, savex1;
+       x0 = savex0 = fStatePhi;
+       x1 = savex1 = fStateZ;
+       x2 = fDt;
+       x3 = fTanL;
+       x4 = fC;
+
+       // Update the state vector
+       x0 += k00*(m[0]-savex0) + k01*(m[1]-savex1);
+       x1 += k10*(m[0]-savex0) + k11*(m[1]-savex1);
+       x2 += k20*(m[0]-savex0) + k21*(m[1]-savex1);
+       x3 += k30*(m[0]-savex0) + k31*(m[1]-savex1);
+       x4 += k40*(m[0]-savex0) + k41*(m[1]-savex1);
+       
+       // Update the covariance matrix
+       Double_t cout00, cout10, cout20, cout30, cout40;
+       Double_t cout11, cout21, cout31, cout41, cout22;
+       Double_t cout32, cout42, cout33, cout43, cout44;
+       
+       cout00 = cin00 - k00*cin00 - k01*cin10;
+       cout10 = cin10 - k00*cin10 - k01*cin11;
+       cout20 = cin20 - k00*cin20 - k01*cin21;
+       cout30 = cin30 - k00*cin30 - k01*cin31;
+       cout40 = cin40 - k00*cin40 - k01*cin41;
+       cout11 = cin11 - k10*cin10 - k11*cin11;
+       cout21 = cin21 - k10*cin20 - k11*cin21;
+       cout31 = cin31 - k10*cin30 - k11*cin31;
+       cout41 = cin41 - k10*cin40 - k11*cin41;
+       cout22 = cin22 - k20*cin20 - k21*cin21;
+       cout32 = cin32 - k20*cin30 - k21*cin31;
+       cout42 = cin42 - k20*cin40 - k21*cin41;
+       cout33 = cin33 - k30*cin30 - k31*cin31;
+       cout43 = cin43 - k30*cin40 - k31*cin41;
+       cout44 = cin44 - k40*cin40 - k41*cin41;
+       
+       // Store the new covariance matrix
+       fMatrix(0,0) = cout00;
+       fMatrix(1,0) = fMatrix(0,1) = cout10;
+       fMatrix(2,0) = fMatrix(0,2) = cout20;
+       fMatrix(3,0) = fMatrix(0,3) = cout30;
+       fMatrix(4,0) = fMatrix(0,4) = cout40;
+       fMatrix(1,1) = cout11;
+       fMatrix(2,1) = fMatrix(1,2) = cout21;
+       fMatrix(3,1) = fMatrix(1,3) = cout31;
+       fMatrix(4,1) = fMatrix(1,4) = cout41;
+       fMatrix(2,2) = cout22;
+       fMatrix(3,2) = fMatrix(2,3) = cout32;
+       fMatrix(4,2) = fMatrix(2,4) = cout42;
+       fMatrix(3,3) = cout33;
+       fMatrix(4,3) = fMatrix(3,4) = cout43;
+       fMatrix(4,4) = cout44;
+       
+       // Calculation of the chi2 increment
+       Double_t vmcold00 = v00 - cout00;
+       Double_t vmcold10 = -cout10;
+       Double_t vmcold11 = v11 - cout11;
+       ////////////////////// Matrix vmc inversion  ///////////////////////
+       det = vmcold00*vmcold11 - vmcold10*vmcold10;
+       Double_t vmc00=vmcold11/det;
+       Double_t vmc10 = -vmcold10/det;
+       Double_t vmc11 = vmcold00/det;
+       ////////////////////////////////////////////////////////////////////
+       Double_t chi2 = (m[0] - x0)*( vmc00*(m[0] - x0) + 2.*vmc10*(m[1] - x1) ) + (m[1] - x1)*vmc11*(m[1] - x1);
+       fChi2 += chi2;
+       fNSteps++;
+       
+       return kTRUE;
+}
+//
+//
+//
+Bool_t AliITSNeuralTrack::KalmanFit()
+// Applies the Kalman Filter to improve the track parameters resolution.
+// First, thre point which lies closer to the estimated helix is chosen.
+// Then, a fit is performed towards the 6th layer
+// Finally, the track is refitted to the 1st layer
+{
+       Double_t rho;
+       Int_t l, layer, sign;
+       
+       fStateR = fPoint[0]->GetR2();
+       fStatePhi = fPoint[0]->GetPhi();
+       fStateZ = fPoint[0]->Z();
+       
+       if (!PropagateTo(3.0)) {
+               Error("KalmanFit", "Unsuccessful initialization");
+               return kFALSE;
+       }
+       l=0;
+
+       // Performs a Kalman filter going from the actual state position
+       // towards layer 6 position
+       // Now, the propagation + filtering operations can be performed
+       Double_t argPhi = 0.0, argZ = 0.0;
+       while (l <= 5) {
+               if (!fPoint[l]) {
+                       Error("KalmanFit", "Not six points!");
+                       return kFALSE;
+               }
+               layer = fPoint[l]->GetLayer();
+               rho = fPoint[l]->GetR2();
+               sign = PropagateTo(rho);
+               if (!sign) return kFALSE;
+               AddEL(layer, -1.0);
+               AddMS(layer);
+               if (!Filter(fPoint[l])) return kFALSE;
+               // these two parameters are update according to the filtered values
+               argPhi = ArgPhi(fStateR);
+               argZ = ArgZ(fStateR);
+               if (argPhi > 1.0 || argPhi < -1.0 || argZ > 1.0 || argZ < -1.0) {
+                       Error("Filter", Form("Filtering returns too large values: %g, %g", argPhi, argZ));
+                       return kFALSE;
+               }
+               fG0 = fStatePhi - asin(argPhi);
+               fDz = fStateZ - (2.0 * fTanL / fC) * asin(argZ);
+               l++;
+       }
+
+       // Now a Kalman filter i performed going from the actual state position
+       // towards layer 1 position and then propagates to vertex
+       if (l >= 5) l = 5;
+       while (l >= 1) {
+               layer = fPoint[l]->GetLayer();
+               rho = fPoint[l]->GetR2();
+               AddEL(layer, 1.0);
+               sign = PropagateTo(rho);
+               if (!sign) return kFALSE;
+               AddMS(layer);
+               if (!Filter(fPoint[l])) return kFALSE;
+               // these two parameters are update according to the filtered values
+               argPhi = ArgPhi(fStateR);
+               argZ = ArgZ(fStateR);
+               if (argPhi > 1.0 || argPhi < -1.0 || argZ > 1.0 || argZ < -1.0) {
+                       Error("Filter", Form("Filtering returns too large values: %g, %g", argPhi, argZ));
+                       return kFALSE;
+               }
+               fG0 = fStatePhi - asin(argPhi);
+               fDz = fStateZ - (2.0 * fTanL / fC) * asin(argZ);
+               l--;
+       }
+       return kTRUE;
+}
+//
+//
+//
+Bool_t AliITSNeuralTrack::RiemannFit()
+{
+       // Method which executes the circle fit via a Riemann Sphere projection
+       // with the only improvement of a weighted mean, due to different errors
+       // over different point measurements.
+       // As an output, it returns kTRUE or kFALSE respectively if the fit succeeded or not
+       // in fact, if some variables assume strange values, the fit is aborted,
+       // in order to prevent the class from raising a floating point error;
+
+       Int_t i, j;
+
+       // M1 - matrix of ones
+       TMatrixD m1(6,1);
+       for (i = 0; i < 6; i++) m1(i,0) = 1.0;
+
+       // X - matrix of Rieman projection coordinates
+       TMatrixD X(6,3);
+       for (i = 0; i < 6; i++) {
+               X(i,0) = fPoint[i]->X();
+               X(i,1) = fPoint[i]->Y();
+               X(i,2) = fPoint[i]->GetR2sq();
+       }
+
+       // W - matrix of weights
+       Double_t xterm, yterm, ex, ey;
+       TMatrixD W(6,6);
+       for (i = 0; i < 6; i++) {
+               xterm = fPoint[i]->X() * fPoint[i]->GetPhi() - fPoint[i]->Y() / fPoint[i]->GetR2();
+               ex = fPoint[i]->ErrX();
+               yterm = fPoint[i]->Y() * fPoint[i]->GetPhi() + fPoint[i]->X() / fPoint[i]->GetR2();
+               ey = fPoint[i]->ErrY();
+               W(i,i) = fPoint[i]->GetR2sq() / (xterm * xterm * ex * ex + yterm * yterm * ey * ey);
+       }
+
+       // Xm - weighted sample mean
+       Double_t Xm = 0.0, Ym = 0.0, Wm = 0.0, sw = 0.0;
+       for (i = 0; i < 6; i++) {
+               Xm += W(i,i) * X(i,0);
+               Ym += W(i,i) * X(i,1);
+               Wm += W(i,i) * X(i,2);
+               sw += W(i,i);
+       }
+       Xm /= sw;
+       Ym /= sw;
+       Wm /= sw;
+
+       // V - sample covariance matrix
+       for (i = 0; i < 6; i++) {
+               X(i,0) -= Xm;
+               X(i,1) -= Ym;
+               X(i,2) -= Wm;
+       }
+       TMatrixD Xt(TMatrixD::kTransposed, X);
+       TMatrixD WX(W, TMatrixD::kMult, X);
+       TMatrixD V(Xt, TMatrixD::kMult, WX);
+       for (i = 0; i < 3; i++) {
+               for (j = i + 1; j < 3; j++) {
+                       V(i,j)  = V(j,i)  = (V(i,j) + V(j,i)) * 0.5;
+               }
+       }
+
+       // Eigenvalue problem solving for V matrix
+       Int_t ileast = 0;
+       TVectorD Eval(3), n(3);
+       TMatrixD Evec = V.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 = -(Xm * n(0) + Ym * n(1) + Wm * n(2));
+
+       fXC = -n(0) / (2. * n(2));
+       fYC = -n(1) / (2. * n(2));
+       fR  = (1. - n(2)*n(2) - 4.*c*n(2)) / (4. * n(2) * n(2));
+
+       if (fR <= 0.E0) {
+               Error("RiemannFit", "Radius comed less than zero!!!");
+               return kFALSE;
+       }
+       fR = TMath::Sqrt(fR);
+       fC = 1.0 / fR;
+
+       // evaluating signs for curvature and others
+       Double_t phi1 = 0.0, phi2, temp1, temp2, sumdphi = 0.0, ref = TMath::Pi();
+       AliITSNeuralPoint *p = fPoint[0];
+       phi1 = p->GetPhi();
+       for (i = 1; i < 6; i++) {
+               p = (AliITSNeuralPoint*)fPoint[i];
+               if (!p) break;
+               phi2 = p->GetPhi();
+               temp1 = phi1;
+               temp2 = phi2;
+               if (temp1 > ref && temp2 < ref)
+                       temp2 += 2.0 * ref;
+               else if (temp1 < ref && temp2 > ref)
+                       temp1 += 2.0 * ref;
+               sumdphi += temp2 - temp1;
+               phi1 = phi2;
+       }
+       if (sumdphi < 0.E0) fC = -fC;
+       Double_t diff, angle = TMath::ATan2(fYC, fXC);
+       if (fC < 0.E0)
+               fG0 = angle + 0.5 * TMath::Pi();
+       else
+               fG0 = angle - 0.5 * TMath::Pi();
+       diff = angle - fG0;
+
+       Double_t D = TMath::Sqrt(fXC*fXC + fYC*fYC) - fR;
+       if (fC >= 0.E0)
+               fDt = D;
+       else
+               fDt = -D;
+
+       Int_t N = 6;
+       Double_t halfC = 0.5 * fC;
+       Double_t *s = new Double_t[N], *z = new Double_t[N], *ws = new Double_t[N];
+       for (j = 0; j < 6; j++) {
+               p = fPoint[j];
+               if (!p) break;
+               s[j] = ArgZ(p->GetR2());
+               if (s[j] > 100.0) return kFALSE;
+               z[j] = p->Z();
+               s[j] = asin(s[j]) / halfC;
+               ws[j] = 1.0 / (p->ErrZ()*p->ErrZ());
+       }
+
+       // second tep final fit
+       Double_t Ss2 = 0.0, Sz = 0.0, Ssz = 0.0, Ss = 0.0, sumw = 0.0;
+       for (i = 0; i < N; i++) {
+               Ss2 += ws[i] * s[i] * s[i];
+               Sz  += ws[i] * z[i];
+               Ss  += ws[i] * s[i];
+               Ssz += ws[i] * s[i] * z[i];
+               sumw += ws[i];
+       }
+       Ss2 /= sumw;
+       Sz /= sumw;
+       Ss /= sumw;
+       Ssz /= sumw;
+       D = Ss2 - Ss*Ss;
+
+       fDz = (Ss2*Sz - Ss*Ssz) / D;
+       fTanL = (Ssz - Ss*Sz) / D;
+
+       delete [] s;
+       delete [] z;
+       delete [] ws;
+
+       return kTRUE;
+}
+//
+//
+//
+void AliITSNeuralTrack::PrintState(Bool_t matrix)
+// Prints the state vector values.
+// The argument switches on or off the printing of the covariance matrix.
+{
+       cout << "\nState vector: " << endl;
+       cout << " Rho = " << fStateR << "\n";
+       cout << " Phi = " << fStatePhi << "\n";
+       cout << "   Z = " << fStateZ << "\n";
+       cout << "  Dt = " << fDt << "\n";
+       cout << "  Dz = " << fDz << "\n";
+       cout << "TanL = " << fTanL << "\n";
+       cout << "   C = " << fC << "\n";
+       cout << "  G0 = " << fG0 << "\n";
+       cout << "  XC = " << fXC << "\n";
+       cout << "  YC = " << fYC << "\n";
+       if (matrix) {
+               cout << "\nCovariance Matrix: " << endl;
+               fMatrix.Print();
+       }
+       cout << "Actual square chi = " << fChi2;
+}
+//
+//
+//
+Double_t AliITSNeuralTrack::GetDz()
+{
+//     Double_t argZ = ArgZ(fStateR);
+//     if (argZ > 9.9) {
+//             Error("GetDz", "Too large value: %g", argZ);
+//             return 0.0;
+//     }
+//     fDz = fStateZ - (2.0 * fTanL / fC) * asin(argZ);
+       return fDz;
+}
+//
+//
+//
+Double_t AliITSNeuralTrack::GetGamma()
+{
+// these two parameters are update according to the filtered values
+//     Double_t argPhi = ArgPhi(fStateR);
+//     if (argPhi > 9.9) {
+//             Error("Filter", "Too large value: %g", argPhi);
+//             return kFALSE;
+//     }
+//     fG0 = fStatePhi - asin(argPhi);
+       return fG0;
+}
+//
+//
+//
+Double_t AliITSNeuralTrack::GetPhi(Double_t r)
+// Gives the value of azymuthal coordinate in the helix
+// as a function of cylindric radius
+{
+       Double_t arg = ArgPhi(r);
+       if (arg > 0.9) return 0.0;
+       arg = fG0 + asin(arg);
+       while (arg >= 2.0 * TMath::Pi()) { arg -= 2.0 * TMath::Pi(); }
+       while (arg < 0.0) { arg += 2.0 * TMath::Pi(); }
+       return arg;
+}
+//
+//
+//
+Double_t AliITSNeuralTrack::GetZ(Double_t r)
+// gives the value of Z in the helix
+// as a function of cylindric radius
+{
+       Double_t arg = ArgZ(r);
+       if (arg > 0.9) return 0.0;
+       return fDz + fTanL * asin(arg) / fC;
+}
+//
+//
+//
+Double_t AliITSNeuralTrack::GetdEdX()
+{
+       Double_t q[4] = {0., 0., 0., 0.}, dedx = 0.0;
+        Int_t i = 0, swap = 0;
+        for (i = 2; i < 6; i++) {
+                if (!fPoint[i]) continue;
+                q[i - 2] = (Double_t)fPoint[i]->GetCharge();
+               q[i - 2] /= (1 + fTanL*fTanL);
+        }
+       q[0] /= 280.;
+       q[1] /= 280.;
+       q[2] /= 38.;
+       q[3] /= 38.;
+       do {
+               swap = 0;
+               for (i = 0; i < 3; i++) {
+                       if (q[i] <= q[i + 1]) continue;
+                       Double_t tmp = q[i];
+                       q[i] = q[i + 1]; 
+                       q[i+1] = tmp;
+                       swap++;
+               }
+       } while(swap); 
+       if(q[0] < 0.) {
+               q[0] = q[1];
+               q[1] = q[2];
+               q[2] = q[3];
+               q[3] = -1.;
+       } 
+       dedx = (q[0] + q[1]) / 2.;
+        return dedx;
+}
+//
+//
+//
+void AliITSNeuralTrack::SetVertex(Double_t *pos, Double_t *err)
+{
+       // Stores vertex data
+
+       if (!pos || !err) return;
+       fVertex.ErrX() = err[0];
+       fVertex.ErrY() = err[1];
+       fVertex.ErrZ() = err[2];
+       fVertex.SetLayer(0);
+       fVertex.SetModule(0);
+       fVertex.SetIndex(0);
+       fVertex.SetLabel(0, -1);
+       fVertex.SetLabel(1, -1);
+       fVertex.SetLabel(2, -1);
+       fVertex.SetUser(1);
+}
+//
+//
+//
+AliITSIOTrack* AliITSNeuralTrack::ExportIOtrack(Int_t min)
+// Exports an object in the standard format for reconstructed tracks
+{
+       Int_t layer = 0;
+       AliITSIOTrack *track = new AliITSIOTrack;
+
+       // covariance matrix
+       track->SetCovMatrix(fMatrix(0,0), fMatrix(1,0), fMatrix(1,1),
+                           fMatrix(2,0), fMatrix(2,1), fMatrix(2,2),
+                           fMatrix(3,0), fMatrix(3,1), fMatrix(3,2),
+                           fMatrix(3,3), fMatrix(4,0), fMatrix(4,1),
+                           fMatrix(4,2), fMatrix(4,3), fMatrix(4,4));
+       
+       // labels
+       track->SetLabel(IsGood(min) ? fLabel : -fLabel);
+       track->SetTPCLabel(-1);
+
+       // points characteristics
+       for (layer = 0; layer < 6; layer++) {
+               if (fPoint[layer]) {
+                       track->SetIdModule(layer, fPoint[layer]->GetModule());
+                       track->SetIdPoint(layer, fPoint[layer]->GetIndex());
+               }
+       }
+       
+       // state vector
+       track->SetStatePhi(fStatePhi);
+       track->SetStateZ(fStateZ);
+       track->SetStateD(fDt);
+       track->SetStateTgl(fTanL);
+       track->SetStateC(fC);
+       track->SetRadius(fStateR);
+       track->SetCharge((fC > 0.0) ? -1 : 1);
+       track->SetDz(fDz);
+
+       // track parameters in the closest point
+       track->SetX(fStateR * cos(fStatePhi));
+       track->SetY(fStateR * cos(fStatePhi));
+       track->SetZ(fStateZ);
+       track->SetPx(GetPt() * cos(fG0));
+       track->SetPy(GetPt() * sin(fG0));
+       track->SetPz(GetPt() * fTanL);
+
+       // PID
+       track->SetPid(fPDG);
+       track->SetMass(fMass);
+
+       return track;
+}
+//
+//
+//====================================================================================
+//============================       PRIVATE METHODS      ============================
+//====================================================================================
+//
+//
+Double_t AliITSNeuralTrack::ArgPhi(Double_t r) const
+{
+       // calculates the expression ((1/2)Cr + (1 + (1/2)CD) D/r) / (1 + CD)
+
+       Double_t arg, num, den;
+       num = (0.5 * fC * r) + (1. + (0.5 * fC * fDt)) * (fDt / r);
+       den = 1. + fC * fDt;
+       if (den == 0.) {
+               Error("ArgPhi", "Denominator = 0!");
+               return 10.0;
+       }
+       arg = num / den;
+       if (TMath::Abs(arg) < 1.) return arg;
+       if (TMath::Abs(arg) <= 1.00001) return (arg > 0.) ? 0.99999999999 : -0.9999999999;
+       Error("ArgPhi", "Value too large: %17.15g", arg);
+       return 10.0;
+}
+//
+//
+//
+Double_t AliITSNeuralTrack::ArgZ(Double_t r) const
+{
+       // calculates the expression (1/2)C * sqrt( (r^2 - Dt^2) / (1 + CD) )
+
+       Double_t arg;
+       arg = (r * r - fDt * fDt) / (1. + fC * fDt);
+       if (arg < 0.) {
+               if (fabs(arg) < 1.E-6) arg = 0.;
+               else {
+                       Error("ArgZ", "Square root argument error: %17.15g < 0", arg);
+                       return 10.;
+               }
+       }
+       arg = 0.5 * fC * TMath::Sqrt(arg);
+       if (TMath::Abs(arg) < 1.) return arg;
+       if (TMath::Abs(arg) <= 1.00001) return (arg > 0.) ? 0.99999999999 : -0.9999999999;
+       Error("ArgZ", "Value too large: %17.15g", arg);
+       return 10.0;
+}
+//
+//
+//
+Double_t AliITSNeuralTrack::ArgB(Double_t r) const 
+{
+       Double_t arg;
+       arg = (r*r - fDt*fDt);
+       arg /= (r*(1.+ fC*fDt)*(1.+ fC*fDt));
+       return arg;
+}
+//
+//
+//
+Double_t AliITSNeuralTrack::ArgC(Double_t r) const 
+{
+       Double_t arg;
+       arg = (1./r - fC * ArgPhi(r));
+       arg /= 1.+ fC*fDt;
+       return arg;
+}
diff --git a/ITS/AliITSNeuralTrack.h b/ITS/AliITSNeuralTrack.h
new file mode 100644 (file)
index 0000000..c9c6e78
--- /dev/null
@@ -0,0 +1,126 @@
+#ifndef ALIITSNEURALTRACK_H
+#define ALIITSNEURALTRACK_H
+
+#include <TMatrixD.h>
+
+class TObjArray;
+class AliITSNeuralPoint;
+//class AliITSVertex;
+class AliITSIOTrack;
+
+class AliITSNeuralTrack : public TObject {
+
+public:
+                AliITSNeuralTrack();
+       virtual ~AliITSNeuralTrack();
+               
+       // Points insertion and goodness evaluation
+
+       void     AssignLabel();
+       void     CleanSlot(Int_t i, Bool_t del = kFALSE);
+       void     CleanAllSlots(Bool_t del = kFALSE)        {Int_t i; for(i=0;i<6;i++) CleanSlot(i,del);}
+       void     GetModuleData(Int_t i, Int_t &mod, Int_t &pos);
+       void     Insert(AliITSNeuralPoint *point);
+       Bool_t   IsGood(Int_t min)                         {return (fCount >= min);}
+       Int_t    OccupationMask();
+       void     PrintLabels();
+               
+       // Fit procedures
+               
+       Bool_t   AddEL(Int_t layer, Double_t sign);
+       Bool_t   AddMS(Int_t layer);
+       void     ForceSign(Double_t sign)            { fC *= sign; }  // externally imposed trach charge
+       void     ResetChi2()                         { fChi2 = fNSteps = 0.0; }
+       Bool_t   SeedCovariance();
+       Int_t    PropagateTo(Double_t r);
+       Bool_t   Filter(AliITSNeuralPoint *test);
+       Bool_t   KalmanFit();
+       Bool_t   RiemannFit();
+       void     PrintState(Bool_t matrix);
+
+       // Getters
+
+       Int_t    GetLabel()        {return fLabel;}
+       Int_t    GetCount()        {return fCount;}
+       Double_t GetDt()           {return fDt;}
+       Double_t GetDz();          
+       Double_t GetC()            {return fC;}
+       Double_t GetR()            {return fR;}
+       Double_t GetXC()           {return fXC;}
+       Double_t GetYC()           {return fYC;}
+       Double_t GetTanL()         {return fTanL;}
+       Double_t GetGamma();       
+       Double_t GetChi2()         {return fChi2;}
+       Double_t GetStateR()       {return fStateR;}
+       Double_t GetStatePhi()     {return fStatePhi;}
+       Double_t GetStateZ()       {return fStateZ;}
+       Double_t GetCovElement(Int_t i, Int_t j) {return fMatrix(i,j);}
+       Double_t GetPhi(Double_t r);        // phi = gamma0 + asin(argphi(rho))
+       Double_t GetZ(Double_t r);          // z = dz + (tanl / C) * asin(argz(rho))
+       
+
+       Double_t GetP()       {return GetPt() * (1.0 + fTanL * fTanL);}
+       Double_t GetPt()      {return 0.299792658 * 0.2 * fField * fabs(1./fC/100.);}
+       Double_t GetPz()      {return GetPt() * fTanL;}
+       Double_t GetE()       {return sqrt(fMass*fMass + GetPt()*GetPt());}
+       Double_t GetLambda()  {return atan(fTanL);}
+       Int_t    GetPDGcode() {return fPDG;}
+       Double_t GetdEdX();
+
+       // Setters
+
+       void     SetFieldFactor(Double_t fact) {fField=fact;}
+       void     SetMass(Double_t mass)        {fMass=mass;}
+       void     SetPDGcode(Int_t code)        {fPDG=code;}
+       void     SetVertex(Double_t *pos, Double_t *err);
+       /*
+       void     SetRho(Double_t a)  {fStateR=a;}
+       void     SetPhi(Double_t a)  {fStatePhi=a;}
+       void     SetZ(Double_t a)    {fStateZ=a;}
+       void     SetDt(Double_t a)   {fDt=a;}
+       void     SetTanL(Double_t a) {fTanL=a;}
+       void     SetC(Double_t a)    {fC=a;}
+       void     SetChi2(Double_t a) {fChi2=a;}
+       void     SetGamma(Double_t a){fG0=a;}
+       void     SetDz(Double_t a)   {fDz=a;}
+       void     SetCovElement(Int_t i, Int_t j, Double_t a) {fMatrix(i,j)=a; if(i!=j) fMatrix(j,i)=a;}
+       */
+       AliITSIOTrack* ExportIOtrack(Int_t min);
+
+private:
+       
+       Double_t ArgPhi(Double_t r) const;
+       Double_t ArgZ  (Double_t r) const;
+       Double_t ArgB  (Double_t r) const;
+       Double_t ArgC  (Double_t r) const;
+
+       Double_t fXC;       // X ofcurvature center
+       Double_t fYC;       // Y of curvature center
+       Double_t fR;        // curvature radius
+       Double_t fC;        // semi-curvature of the projected circle (signed)
+       Double_t fTanL;     // tangent of dip angle
+       Double_t fG0;       // phase coefficient
+       Double_t fDt;       // transverse impact parameter
+       Double_t fDz;       // longitudinal impact parameter
+
+       Double_t fStateR;   //
+       Double_t fStatePhi; // state vector coordinates
+       Double_t fStateZ;   //
+       TMatrixD fMatrix;   // covariance matrix
+       Double_t fChi2;     // square chi (calculated by Kalman filter)
+       Double_t fNSteps;   // number of Kalman steps
+
+       Double_t fMass;     // the particle mass
+       Double_t fField;    // B field = 0.2 * fField (Tesla)
+               
+       Int_t    fPDG;      // PDG code of the recognized particle
+       Int_t    fLabel;    // the GEANT label most appearing among track recpoints
+       Int_t    fCount;    // number of counts of above label
+               
+       AliITSNeuralPoint  fVertex;   // vertex position data
+       AliITSNeuralPoint *fPoint[6]; // track points
+
+       ClassDef(AliITSNeuralTrack, 1)
+};
+
+#endif
diff --git a/ITS/AliITSNeuralTracker.cxx b/ITS/AliITSNeuralTracker.cxx
new file mode 100644 (file)
index 0000000..3962318
--- /dev/null
@@ -0,0 +1,1233 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: Alberto Pulvirenti.                                            *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose.                *
+ * It is provided "as is" without express or implied warranty.            *
+ *                                                                        *
+ * AN ITS STAND-ALONE "NEURAL" TRACK FINDER                               *
+ * ----------------------------------------                               *
+ * This class implements the Denby-Peterson algorithm for track finding   *
+ * in the ITS stand-alone, by means of a neural network simulation.       *
+ * Many parameters have to be set for the neural network to operate       *
+ * correctly and with a good efficiency.                                  *
+ * The neural tracker must be feeded with a TTree filled with objects     *
+ * of the class "AliITSNeuralPoint", into a single branch called       *
+ * "Points".                                                              *
+ **************************************************************************/
+
+//#define NEURAL_LINEAR
+#include <fstream>
+#include <Riostream.h>
+#include <stdlib.h>
+
+#include <TROOT.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TMath.h>
+#include <TLine.h>
+#include <TMarker.h>
+#include <TRandom.h>
+#include <TString.h>
+#include <TCanvas.h>
+#include <TVector3.h>
+#include <TParticle.h>
+#include <TObjArray.h>
+#include <TList.h>
+
+#include "AliITSNeuralPoint.h"
+#include "AliITSNeuralTracker.h"
+
+using namespace std;
+ClassImp(AliITSNeuralTracker)
+
+//--------------------------------------------------------------------------------------------
+
+AliITSNeuralTracker::AliITSNeuralTracker()
+{
+       // CONSTRUCTOR
+       //
+       // Initializes some data-members:
+       // - all pointers to NULL
+       // - theta cut to 180 deg. (= no theta cut)
+       // - initial choice of only 1 azimuthal sector
+       // - initial values for the neural network parameters.
+       //
+       // With these settings the neural tracker can't work
+       // because it has not any curvature cut set.
+       
+       fCurvNum = 0;
+       fCurvCut = 0;
+
+       fSectorNum   = 1;
+       fSectorWidth = TMath::Pi() * 2.0;
+       fPolarInterval = 45.0;
+
+       fStructureOK = kFALSE;
+
+       fVX = fVY = fVZ = 0.0;
+
+       Int_t ilayer, itheta;
+       for (ilayer = 0; ilayer < 6; ilayer++) {
+               for (itheta = 0; itheta < 180; itheta++) fPoints[ilayer][itheta] = 0;
+               fThetaCut2DMin[ilayer] = 0.0;
+               fThetaCut2DMax[ilayer] = TMath::Pi();
+               fThetaCut3DMin[ilayer] = 0.0;
+               fThetaCut3DMax[ilayer] = TMath::Pi();
+               fHelixMatchCutMin[ilayer] = 1.0;
+               fHelixMatchCutMax[ilayer] = 1.0;
+       }
+
+       fEdge1 = 0.3;
+       fEdge2=  0.7;
+
+       fTemperature = 1.0;
+       fStabThreshold = 0.001;
+       fGain2CostRatio = 1.0;
+       fAlignExponent = 1.0;
+       fActMinimum = 0.5;
+
+       fNeurons = 0;
+
+       fChains = new TTree("TreeC", "Sextines of points");
+       fChains->Branch("l0", &fPoint[0], "l0/I");
+       fChains->Branch("l1", &fPoint[1], "l1/I");
+       fChains->Branch("l2", &fPoint[2], "l2/I");
+       fChains->Branch("l3", &fPoint[3], "l3/I");
+       fChains->Branch("l4", &fPoint[4], "l4/I");
+       fChains->Branch("l5", &fPoint[5], "l5/I");
+}
+//
+//--------------------------------------------------------------------------------------------
+//
+AliITSNeuralTracker::~AliITSNeuralTracker()
+{
+       // DESTRUCTOR
+       //
+       // It Destroys all the dynamic arrays and
+       // clears the TCollections and the points tree
+       
+       cout << "Starting destructor..." << endl;
+       
+       delete [] fCurvCut;
+       
+       Int_t ilayer, itheta;
+       if (fStructureOK) {
+               cout << "Deleting points..." << endl;
+               for (ilayer = 0; ilayer < 6; ilayer++) {
+                       for (itheta = 0; itheta < 180; itheta++) {
+                               fPoints[ilayer][itheta]->SetOwner();
+                               delete fPoints[ilayer][itheta];
+                       }
+               }
+               cout << "Deleting neurons..." << endl;
+               fNeurons->SetOwner();
+               delete fNeurons;
+       }
+       
+       cout << "AliITSNeuralTracker destructed completely!" << endl;
+}
+//
+//--------------------------------------------------------------------------------------------
+//
+void AliITSNeuralTracker::Display(TCanvas*& canv)
+{
+       Double_t x1, y1, x2, y2;
+       canv->Clear();
+       TObjArrayIter iter(fNeurons);
+       for (;;) {
+               AliITSneuron *unit = (AliITSneuron*)iter.Next();
+               if (!unit) break;
+               if (unit->fActivation < fActMinimum) continue;
+               x1 = unit->fInner->X();
+               x2 = unit->fOuter->X();
+               y1 = unit->fInner->Y();
+               y2 = unit->fOuter->Y();
+               TLine *line = new TLine(x1, y1, x2, y2);
+               canv->cd();
+               line->Draw();
+       }
+       canv->Update();
+}
+//
+//--------------------------------------------------------------------------------------------
+//
+void AliITSNeuralTracker::SetThetaCuts2D(Double_t *min, Double_t *max)
+{
+       Int_t i;
+       Double_t temp;
+       for (i = 0; i < 5; i++) {
+               if (min[i] > max[i]) {
+                       temp = min[i];
+                       min[i] = max[i];
+                       max[i] = temp;
+               }
+               fThetaCut2DMin[i] = min[i];
+               fThetaCut2DMax[i] = max[i];
+       }
+       for (i = 0; i < 5; i++) {
+               cout << "Theta 2D cut for layer " << i << " = " << fThetaCut2DMin[i] << " to " << fThetaCut2DMax[i] << endl;
+       }
+}
+//
+//--------------------------------------------------------------------------------------------
+//
+void AliITSNeuralTracker::SetThetaCuts3D(Double_t *min, Double_t *max)
+{
+       Int_t i;
+       Double_t temp;
+       for (i = 0; i < 5; i++) {
+               if (min[i] > max[i]) {
+                       temp = min[i];
+                       min[i] = max[i];
+                       max[i] = temp;
+               }
+               fThetaCut3DMin[i] = min[i];
+               fThetaCut3DMax[i] = max[i];
+       }
+       for (i = 0; i < 5; i++) {
+               cout << "Theta 3D cut for layer " << i << " = " << fThetaCut3DMin[i] << " to " << fThetaCut3DMax[i] << endl;
+       }
+}
+//
+//--------------------------------------------------------------------------------------------
+//
+void AliITSNeuralTracker::SetHelixMatchCuts(Double_t *min, Double_t *max)
+{
+       Int_t i;
+       Double_t temp;
+       for (i = 0; i < 5; i++) {
+               if (min[i] > max[i]) {
+                       temp = min[i];
+                       min[i] = max[i];
+                       max[i] = temp;
+               }
+               fHelixMatchCutMin[i] = min[i];
+               fHelixMatchCutMax[i] = max[i];
+       }
+       for (i = 0; i < 5; i++) {
+               cout << "Helix-match cut for layer " << i << " = " << fHelixMatchCutMin[i] << " to " << fHelixMatchCutMax[i] << endl;
+       }
+}
+//
+//--------------------------------------------------------------------------------------------
+//
+void AliITSNeuralTracker::SetCurvatureCuts(Int_t ncuts, Double_t *curv)
+{
+       // CURVATURE CUTS SETTER
+       //
+       // Requires an array of double values and its dimension
+       // After sorting it in increasing order, the array of curvature cuts
+       // is dinamically allocated, and filled with the sorted cuts array.
+       // A message is shown which lists all the curvature cuts.
+
+       Int_t i, *ind = new Int_t[ncuts];
+       TMath::Sort(ncuts, curv, ind, kFALSE);
+       fCurvCut = new Double_t[ncuts];
+       cout << "\n" << ncuts << " curvature cuts defined" << endl << "-----" << endl;
+       for (i = 0; i < ncuts; i++) {
+               fCurvCut[i] = curv[ind[i]];
+               cout << "Cut #" << i + 1 << " : " << fCurvCut[i] << endl;
+       }
+       cout << "-----" << endl;
+       fCurvNum = ncuts;
+}
+//
+//--------------------------------------------------------------------------------------------
+//
+void AliITSNeuralTracker::CreateArrayStructure(Int_t nsectors)
+{
+       // ARRAY CREATOR
+       //
+       // Organizes the array structure to store all points in.
+       //
+       // The array is organized into a "multi-level" TCollection:
+       // - 6 fPoints[] TObjArray containing a TObjArray for each layer
+       // - each TObject contains a TObjArray for each sector.
+
+       // sets the number of sectors and their width.
+       fSectorNum   = nsectors;
+       fSectorWidth = TMath::Pi() * 2.0 / (Double_t)fSectorNum;
+
+       // creates the TCollection structure
+       Int_t ilayer, isector, itheta;
+       TObjArray *sector = 0;
+       for (ilayer = 0; ilayer < 6; ilayer++) {
+               for (itheta = 0; itheta < 180; itheta++) {
+                       if (fPoints[ilayer][itheta]) {
+                               fPoints[ilayer][itheta]->SetOwner();
+                               delete fPoints[ilayer][itheta];
+                       }
+                       fPoints[ilayer][itheta] = new TObjArray(nsectors);
+                       for (isector = 0; isector < nsectors; isector++) {
+                               sector = new TObjArray;
+                               sector->SetOwner();
+                               fPoints[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 AliITSNeuralTracker::ArrangePoints(TTree* pts_tree)
+{
+       // POINTS STORAGE INTO ARRAY
+       //
+       // Reads points from the tree and creates AliITSNode objects for each one,
+       // storing them into the array structure defined above.
+       // Returns the number of points collected (if successful) or 0 (otherwise)
+
+       // check: if the points tree is NULL or empty, there is nothing to do...
+       if ( !pts_tree || (pts_tree && !(Int_t)pts_tree->GetEntries()) ) {
+               Error("ArrangePoints", "Points tree is NULL or empty: no points to arrange");
+               return 0;
+       }
+
+       if (!fStructureOK) {
+               Error("ArrangePoints", "Structure NOT defined. Call CreateArrayStructure() first");
+               return 0;
+       }
+
+       Int_t isector, itheta, ientry, ilayer, nentries, pos;
+       TObjArray *sector = 0;
+       AliITSNode *created = 0;
+       AliITSNeuralPoint *cursor = 0;
+
+       pts_tree->SetBranchAddress("pos", &pos);
+       pts_tree->SetBranchAddress("Points", &cursor);
+       nentries = (Int_t)pts_tree->GetEntries();
+
+       for (ientry = 0; ientry < nentries; ientry++) {
+               pts_tree->GetEntry(ientry);
+               // creates the object
+               created = new AliITSNode(cursor, kTRUE);
+               created->SetUser(-1);
+               created->fPosInTree = pos;
+               // finds the sector in phi
+               isector = created->GetSector(fSectorWidth);
+               itheta  = created->GetThetaCell();
+               ilayer  = created->GetLayer();
+               if (ilayer < 0 || ilayer > 5) {
+                       Error("ArrangePoints", "Layer value %d not allowed. Aborting...", ilayer);
+                       return 0;
+               }
+               // selects the right TObjArray to store object into
+               sector = (TObjArray*)fPoints[ilayer][itheta]->At(isector);
+               sector->AddLast(created);
+       }
+
+       // returns the number of points processed
+       return ientry;
+}
+//
+//--------------------------------------------------------------------------------------------
+//
+void AliITSNeuralTracker::PrintPoints()
+{
+       // creates the TCollection structure
+       TObjArray *sector = 0;
+       Int_t ilayer, isector, itheta;
+       fstream file("test.log", ios::out);
+       for (ilayer = 0; ilayer < 6; ilayer++) {
+               for (isector = 0; isector < fSectorNum; isector++) {
+                       for (itheta = 0; itheta < 180; itheta++) {
+                               sector = (TObjArray*)fPoints[ilayer][itheta]->At(isector);
+                               file << ilayer << " " << isector << " " << itheta;
+                               file << " " << sector->GetSize() << " points" << endl;
+                       }
+               }
+       }
+       file.close();
+}
+//
+//--------------------------------------------------------------------------------------------
+//
+Bool_t AliITSNeuralTracker::PassCurvCut
+(AliITSNode *p1, AliITSNode *p2, Int_t curv_index, Double_t vx, Double_t vy, Double_t vz)
+{
+       // CURVATURE CUT EVALUATOR
+       //
+       // Checks the passsed point pair w.r. to the current curvature cut
+       // Returns the result of the check.
+
+       if (curv_index < 0 || curv_index >= fCurvNum) {
+               Error("PassCurvCut", "Curv index %d out of range", curv_index);
+               return kFALSE;
+       }
+       
+       // Find the reference layer
+       Int_t lay1 = p1->GetLayer();
+       Int_t lay2 = p2->GetLayer();
+       Int_t ref_layer = (lay1 < lay2) ? lay1 : lay2;
+
+       Double_t x1 = p1->X() - vx;
+       Double_t x2 = p2->X() - vx;
+       Double_t y1 = p1->Y() - vy;
+       Double_t y2 = p2->Y() - vy;
+       Double_t z1 = p1->Z() - vz;
+       Double_t z2 = p2->Z() - vz;
+       Double_t r1 = sqrt(x1*x1 + y1*y1);
+       Double_t r2 = sqrt(x2*x2 + y2*y2);
+
+       // calculation of curvature
+       Double_t dx = p1->X() - p2->X(), dy = p1->Y() - p2->Y();
+       Double_t num = 2 * (x1*y2 - x2*y1);
+       Double_t den = r1*r2*sqrt(dx*dx + dy*dy);
+       Double_t curv = 0.;
+       /* FOR OLD VERSION
+       if (den != 0.) {
+               curv = fabs(num / den);
+               if (curv > fCurvCut[curv_index]) return kFALSE;
+               return kTRUE;
+       }
+       else
+               return kFALSE;
+       */
+       // NEW VERSION
+       if (den != 0.) {
+               curv = fabs(num / den);
+               if (curv > fCurvCut[curv_index]) return kFALSE;
+       }
+       else
+               return kFALSE;
+       // calculation of helix matching
+       Double_t arc1 = 2.0 * r1 * curv;
+       Double_t arc2 = 2.0 * r2 * curv;
+       Double_t hel_match = 0.0;
+       if (arc1 > -1.0 && arc1 < 1.0) arc1 = asin(arc1);
+       else arc1 = ((arc1 > 0.0) ? 0.5 : 1.5) * TMath::Pi();
+       if (arc2 > -1.0 && arc2 < 1.0) arc2 = asin(arc2);
+       else arc2 = ((arc2 > 0.0) ? 0.5 : 1.5) * TMath::Pi();
+       arc1 /= 2.0 * curv;
+       arc2 /= 2.0 * curv;
+       if (arc1 == 0.0 || arc2 == 0.0) return kFALSE;
+       hel_match = fabs(z1 / arc1 - z2 / arc2);
+       return (hel_match >= fHelixMatchCutMin[ref_layer] && hel_match <= fHelixMatchCutMax[ref_layer]);
+       // END NEW VERSION
+}
+//
+//--------------------------------------------------------------------------------------------
+//
+Int_t AliITSNeuralTracker::PassAllCuts
+(AliITSNode *p1, AliITSNode *p2, Int_t curv_index, Double_t vx, Double_t vy, Double_t vz)
+{
+       // GLOBAL CUT EVALUATOR
+       //
+       // Checks all cuts for the passed point pair.
+       // 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
+
+       if (curv_index < 0 || curv_index >= fCurvNum) return 6;
+
+       // Find the reference layer
+       Int_t lay1 = p1->GetLayer();
+       Int_t lay2 = p2->GetLayer();
+       Int_t ref_layer = (lay1 < lay2) ? lay1 : lay2;
+       
+       // Swap points in order that r1 < r2
+       AliITSNode *temp = 0;
+       if (p2->GetLayer() < p1->GetLayer()) {
+               temp = p2;
+               p2 = p1;
+               p1 = temp;
+       }
+
+       // shift XY coords to the reference to the vertex position,
+       // for easier calculus of quantities.
+       Double_t x1 = p1->X() - vx;
+       Double_t x2 = p2->X() - vx;
+       Double_t y1 = p1->Y() - vy;
+       Double_t y2 = p2->Y() - vy;
+       Double_t z1 = p1->Z() - vz;
+       Double_t z2 = p2->Z() - vz;
+       Double_t r1 = sqrt(x1*x1 + y1*y1);
+       Double_t r2 = sqrt(x2*x2 + y2*y2);
+       
+       // Check for theta cut
+       Double_t dtheta, dtheta3;
+       TVector3 v01(z1, r1, 0.0);
+       TVector3 v12(z2 - z1, r2 - r1, 0.0);
+       dtheta = v01.Angle(v12) * 180.0 / TMath::Pi();
+       TVector3 V01(x1, y1, z1);
+       TVector3 V12(x2 - x1, y2 - y1, z2 - z1);
+       dtheta3 = V01.Angle(V12) * 180.0 / TMath::Pi();
+       if (dtheta < fThetaCut2DMin[ref_layer] || dtheta > fThetaCut2DMax[ref_layer]) return 1;
+       if (dtheta3 < fThetaCut3DMin[ref_layer] || dtheta3 > fThetaCut3DMax[ref_layer]) return 2;
+
+       // calculation of curvature
+       Double_t dx = x1 - x2, dy = y1 - y2;
+       Double_t num = 2 * (x1*y2 - x2*y1);
+       Double_t den = r1*r2*sqrt(dx*dx + dy*dy);
+       Double_t curv = 0.;
+       if (den != 0.) {
+               curv = fabs(num / den);
+               if (curv > fCurvCut[curv_index]) return 3;
+       }
+       else
+               return 4;
+
+       // calculation of helix matching
+       Double_t arc1 = 2.0 * r1 * curv;
+       Double_t arc2 = 2.0 * r2 * curv;
+       Double_t hel_match = 0.0;
+       if (arc1 > -1.0 && arc1 < 1.0) arc1 = asin(arc1);
+       else arc1 = ((arc1 > 0.0) ? 0.5 : 1.5) * TMath::Pi();
+       if (arc2 > -1.0 && arc2 < 1.0) arc2 = asin(arc2);
+       else arc2 = ((arc2 > 0.0) ? 0.5 : 1.5) * TMath::Pi();
+       arc1 /= 2.0 * curv;
+       arc2 /= 2.0 * curv;
+       if (arc1 == 0.0 || arc2 == 0.0) return kFALSE;
+       hel_match = fabs(z1 / arc1 - z2 / arc2);
+       if (hel_match < fHelixMatchCutMin[ref_layer] || hel_match > fHelixMatchCutMax[ref_layer]) return 5;
+       
+       return 0;
+}
+//
+//--------------------------------------------------------------------------------------------
+//
+void AliITSNeuralTracker::StoreAbsoluteMatches()
+{
+       // Stores in the 'fMatches' array of each node all the points in the
+       // adjacent layers which allow to create neurons accordin to the
+       // helix and theta cut, and also to the largest curvature cut
+
+       Int_t ilayer, isector, itheta1, itheta2, check;
+       TObjArray *list1 = 0, *list2 = 0;
+       AliITSNode *node1 = 0, *node2 = 0;
+       Double_t theta_min, theta_max;
+       Int_t imin, imax;
+
+       for (isector = 0; isector < fSectorNum; isector++) {
+               // sector is chosen once for both lists
+               for (ilayer = 0; ilayer < 5; ilayer++) {
+                       for (itheta1 = 0; itheta1 < 180; itheta1++) {
+                               list1 = (TObjArray*)fPoints[ilayer][itheta1]->At(isector);
+                               TObjArrayIter iter1(list1);
+                               while ( (node1 = (AliITSNode*)iter1.Next()) ) {
+                                       if (node1->GetUser() >= 0) continue;
+                                       node1->fMatches->Clear();
+                                       theta_min = node1->ThetaDeg() - fPolarInterval;
+                                       theta_max = node1->ThetaDeg() + fPolarInterval;
+                                       imin = (Int_t)theta_min;
+                                       imax = (Int_t)theta_max;
+                                       if (imin < 0) imin = 0;
+                                       if (imax > 179) imax = 179;
+                                       for (itheta2 = imin; itheta2 <= imax; itheta2++) {
+                                               list2 = (TObjArray*)fPoints[ilayer + 1][itheta2]->At(isector);
+                                               TObjArrayIter iter2(list2);
+                                               while ( (node2 = (AliITSNode*)iter2.Next()) ) {
+                                                       check = PassAllCuts(node1, node2, fCurvNum - 1, fVX, fVY, fVZ);
+                                                       switch (check) {
+                                                               case 0:
+                                                                       node1->fMatches->AddLast(node2);
+                                                                       break;
+                                                               case 1:
+                                                                       //Info("StoreAbsoluteMatches", "Layer %d: THETA 2D cut not passed", ilayer);
+                                                                       break;
+                                                               case 2:
+                                                                       //Info("StoreAbsoluteMatches", "Layer %d: THETA 3D cut not passed", ilayer);
+                                                                       break;
+                                                               case 3:
+                                                                       //Info("StoreAbsoluteMatches", "Layer %d: CURVATURE cut not passed", ilayer);
+                                                                       break;
+                                                               case 4:
+                                                                       //Info("StoreAbsoluteMatches", "Layer %d: Division by ZERO in curvature evaluation", ilayer);
+                                                                       break;
+                                                               case 5:
+                                                                       //Info("StoreAbsoluteMatches", "Layer %d: HELIX-MATCH cut not passed", ilayer);
+                                                                       break;
+                                                               case 6:
+                                                                       //Error("PassAllCuts", "Layer %d: Curv index out of range", ilayer);
+                                                                       break;
+                                                               default:
+                                                                       Warning("StoreAbsoluteMatches", "Layer %d: %d: unrecognized return value", ilayer, check);
+                                                       }
+                                               } // while (node2...)
+                                       } // for (itheta2...)
+                               } // while (node1...)
+                       } // for (itheta...)
+               } // for (ilayer...)
+       } // for (isector...)
+}
+//
+//--------------------------------------------------------------------------------------------
+//
+void AliITSNeuralTracker::PrintMatches(Bool_t stop)
+{
+       TFile *ft = new TFile("its_findables_v1.root");
+       TTree *tt = (TTree*)ft->Get("Tracks");
+       Int_t it, nP, nU, lab, nF = (Int_t)tt->GetEntries();
+       tt->SetBranchAddress("nhitsP", &nP);
+       tt->SetBranchAddress("nhitsU", &nU);
+       tt->SetBranchAddress("label", &lab);
+       TString strP("|"), strU("|");
+       for (it = 0; it < nF; it++) {
+               tt->GetEntry(it);
+               if (nP >= 5) strP.Append(Form("%d|", lab));
+               if (nU >= 5) strU.Append(Form("%d|", lab));
+       }
+
+       TObjArray *sector = 0;
+       Int_t ilayer, isector, itheta, id[3];
+       AliITSNode *node1 = 0, *node2 = 0;
+
+       for (ilayer = 0; ilayer < 6; ilayer++) {
+               for (isector = 0; isector < fSectorNum; isector++) {
+                       for (itheta = 0; itheta < 180; itheta++) {
+                               sector = (TObjArray*)fPoints[ilayer][itheta]->At(isector);
+                               TObjArrayIter points(sector);
+                               while ( (node1 = (AliITSNode*)points.Next()) ) {
+                                       for (it = 0; it < 3; it++) id[it] = node1->GetLabel(it);
+                                       nF = (Int_t)node1->fMatches->GetSize();
+                                       cout << "Node layer: " << node1->GetLayer() << ", labels: ";
+                                       cout << id[0] << " " << id[1] << " " << id[2] << " --> ";
+                                       if (!nF) {
+                                               cout << "NO MatchES!!!" << endl;
+                                               continue;
+                                       }
+                                       cout << nF << " Matches" << endl;
+                                       for (it = 0; it < 3; it++) {
+                                               if (strP.Contains(Form("|%d|", id[it])))
+                                                       cout << "Belongs to findable (originary) track " << id[it] << endl;
+                                               if (strU.Contains(Form("|%d|", id[it])))
+                                                       cout << "Belongs to findable (post-Kalman) track " << id[it] << endl;
+                                       }
+                                       TObjArrayIter Matches(node1->fMatches);
+                                       while ( (node2 = (AliITSNode*)Matches.Next()) ) {
+                                               cout << "Match with " << node2;
+                                               Int_t *sh = node1->SharedID(node2);
+                                               for (Int_t k = 0; k < 3; k++)
+                                                       if (sh[k] > -1) cout << " " << sh[k];
+                                               cout << endl;
+                       }
+                                       if (stop) {
+                                               cout << "Press a key" << endl;
+                                               cin.get();
+                                       }
+                               }
+                       }
+               }
+       }
+}
+//
+//--------------------------------------------------------------------------------------------
+//
+void AliITSNeuralTracker::ResetNodes(Int_t isector)
+{
+       Int_t ilayer, itheta;
+       TObjArray *sector = 0;
+       AliITSNode *node = 0;
+       for (ilayer = 0; ilayer < 6; ilayer++) {
+               for (itheta = 0; itheta < 180; itheta++) {
+                       sector = (TObjArray*)fPoints[ilayer][itheta]->At(isector);
+                       TObjArrayIter iter(sector);
+                       for (;;) {
+                               node = (AliITSNode*)iter.Next();
+                               if (!node) break;
+                               node->fInnerOf->Clear();
+                               node->fOuterOf->Clear();
+                               delete node->fInnerOf;
+                               delete node->fOuterOf;
+                               node->fInnerOf = new TObjArray;
+                               node->fOuterOf = new TObjArray;
+                       }
+               }
+       }
+}
+//
+//--------------------------------------------------------------------------------------------
+//
+void AliITSNeuralTracker::NeuralTracking(const char* filesave, TCanvas*& display)
+// This is the public method that operates the tracking.
+// It works sector by sector, and at the end  saves the found tracks.
+// Other methods are privare because they have no meaning id used alone,
+// and sometimes they could get segmentation faults due to uninitialized
+// datamembert they require to work on.
+// The argument is a file where the final results have to be stored.
+{
+       Bool_t isStable = kFALSE;
+       Int_t i, nUnits = 0, nLinks = 0, nTracks = 0, sectTracks = 0, totTracks = 0;
+
+       // tracking through sectors
+       cout << endl;
+       Int_t sector, curv;
+       for (sector = 0; sector < fSectorNum; sector++) {
+               cout << "\rSector " << sector << ": " << endl;
+               sectTracks = 0;
+               for(curv = 0; curv < fCurvNum; curv++) {
+                       cout << "- curvature step " << curv + 1;
+                       cout << " (" << fCurvCut[curv] << "): ";
+                       // units creation
+                       nUnits = CreateNeurons(sector, curv);
+                       if (!nUnits) {
+                               cout << "no neurons --> skipped" << endl;
+                               continue;
+                       }
+                       cout << nUnits << " units, " << flush;
+                       // units connection
+                       nLinks = LinkNeurons();
+                       if (!nLinks) {
+                               cout << "no connections --> skipped" << endl;
+                               continue;
+                       }
+                       cout << nLinks << " connections, " << flush;
+                       // neural updating
+                       for (i = 0;; i++) {
+                               isStable = Update();
+                               if (display) Display(display);
+                               TObjArrayIter iter(fNeurons);
+                               for (;;) {
+                                       AliITSneuron *n = (AliITSneuron*)iter.Next();
+                                       if (!n) break;
+                               }
+                               if (isStable) break;
+                       }
+                       cout << i << " cycles --> " << flush;
+                       // tracks saving
+                       CleanNetwork();
+                       nTracks = Save(sector);
+                       cout << nTracks << " tracks" << endl;
+                       sectTracks += nTracks;
+               }
+               totTracks += sectTracks;
+               //cout << sectTracks << " tracks found (total = " << totTracks << ")      " << flush;
+       }
+
+       cout << endl << totTracks << " tracks found!!!" << endl;
+       cout << endl << "Saving results in file " << filesave << "..." << flush;
+       TFile *f = new TFile(filesave, "recreate");
+       fChains->Write("TreeC");
+       f->Close();
+}
+//
+//--------------------------------------------------------------------------------------------
+//
+Int_t AliITSNeuralTracker::CreateNeurons(Int_t sector_idx, Int_t curv_idx)
+// Fills the neuron arrays, following the cut criteria for the selected step
+// (secnum = sector to analyze, curvnum = curvature cut step to use)
+// It also sets the initial random activation.
+// In order to avoid a large number of 'new' operations, all existing neurons
+// are reset and initialized with the new values, and are created new unit only if
+// it turns out to be necessary
+// the 'flag' argument is used to decide if the lower edge in the intevral
+// of the accepted curvatures is given by zero (kFALSE) or by the preceding used cut (kTRUE)
+// (if this is the first step, anyway, the minimum is always zero)
+{
+       ResetNodes(sector_idx);
+
+       if (fNeurons) delete fNeurons;
+       fNeurons = new TObjArray;
+
+       AliITSneuron *unit = 0;
+       Int_t itheta, neurons = 0;
+       TObjArray *lst_sector = 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++) {
+               lst_sector = (TObjArray*)fPoints[0][itheta]->At(sector_idx);
+               TObjArrayIter lay0(lst_sector);
+               while ( (p[0] = (AliITSNode*)lay0.Next()) ) {
+                       if (p[0]->GetUser() >= 0) continue;
+                       vx[0] = fVX;
+                       vy[0] = fVY;
+                       vz[0] = fVZ;
+                       TObjArrayIter lay1(p[0]->fMatches);
+                       while ( (p[1] = (AliITSNode*)lay1.Next()) ) {
+                               if (p[1]->GetUser() >= 0) continue;
+                               if (!PassCurvCut(p[0], p[1], curv_idx, fVX, fVY, fVZ)) continue;
+                               unit = new AliITSneuron;
+                               unit->fInner = p[0];
+                               unit->fOuter = p[1];
+                               unit->fActivation = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
+                               unit->fGain = new TObjArray;
+                               fNeurons->AddLast(unit);
+                               p[0]->fInnerOf->AddLast(unit);
+                               p[1]->fOuterOf->AddLast(unit);
+                               neurons++;
+                               vx[1] = p[0]->X();
+                               vy[1] = p[0]->Y();
+                               vz[1] = p[0]->Z();
+                               TObjArrayIter lay2(p[1]->fMatches);
+                               while ( (p[2] = (AliITSNode*)lay2.Next()) ) {
+                                       if (p[2]->GetUser() >= 0) continue;
+                                       if (!PassCurvCut(p[1], p[2], curv_idx, vx[1], vy[1], vz[1])) continue;
+                                       unit = new AliITSneuron;
+                                       unit->fInner = p[1];
+                                       unit->fOuter = p[2];
+                                       unit->fActivation = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
+                                       unit->fGain = new TObjArray;
+                                       fNeurons->AddLast(unit);
+                                       p[1]->fInnerOf->AddLast(unit);
+                                       p[2]->fOuterOf->AddLast(unit);
+                                       neurons++;
+                                       vx[2] = p[1]->X();
+                                       vy[2] = p[1]->Y();
+                                       vz[2] = p[1]->Z();
+                                       TObjArrayIter lay3(p[2]->fMatches);
+                                       while ( (p[3] = (AliITSNode*)lay3.Next()) ) {
+                                               if (p[3]->GetUser() >= 0) continue;
+                                               if (!PassCurvCut(p[2], p[3], curv_idx, vx[2], vy[2], vz[2])) continue;
+                                               unit = new AliITSneuron;
+                                               unit->fInner = p[2];
+                                               unit->fOuter = p[3];
+                                               unit->fActivation = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
+                                               unit->fGain = new TObjArray;
+                                               fNeurons->AddLast(unit);
+                                               p[2]->fInnerOf->AddLast(unit);
+                                               p[3]->fOuterOf->AddLast(unit);
+                                               neurons++;
+                                               vx[3] = p[2]->X();
+                                               vy[3] = p[2]->Y();
+                                               vz[3] = p[2]->Z();
+                                               TObjArrayIter lay4(p[3]->fMatches);
+                                               while ( (p[4] = (AliITSNode*)lay4.Next()) ) {
+                                                       if (p[4]->GetUser() >= 0) continue;
+                                                       if (!PassCurvCut(p[3], p[4], curv_idx, vx[3], vy[3], vz[3])) continue;
+                                                       unit = new AliITSneuron;
+                                                       unit->fInner = p[3];
+                                                       unit->fOuter = p[4];
+                                                       unit->fActivation = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
+                                                       unit->fGain = new TObjArray;
+                                                       fNeurons->AddLast(unit);
+                                                       p[3]->fInnerOf->AddLast(unit);
+                                                       p[4]->fOuterOf->AddLast(unit);
+                                                       neurons++;
+                                                       vx[4] = p[3]->X();
+                                                       vy[4] = p[3]->Y();
+                                                       vz[4] = p[3]->Z();
+                                                       TObjArrayIter lay5(p[4]->fMatches);
+                                                       while ( (p[5] = (AliITSNode*)lay5.Next()) ) {
+                                                               if (p[5]->GetUser() >= 0) continue;
+                                                               if (!PassCurvCut(p[4], p[5], curv_idx, vx[4], vy[4], vz[4])) continue;
+                                                               unit = new AliITSneuron;
+                                                               unit->fInner = p[4];
+                                                               unit->fOuter = p[5];
+                                                               unit->fActivation = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
+                                                               unit->fGain = new TObjArray;
+                                                               fNeurons->AddLast(unit);
+                                                               p[4]->fInnerOf->AddLast(unit);
+                                                               p[5]->fOuterOf->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++) {
+                       lst_sector = (TObjArray*)fPoints[ilayer][itheta]->At(sector_idx);
+                       TObjArrayIter inners(lst_sector);
+                       while ( (inner = (AliITSNode*)inners.Next()) ) {
+                               if (inner->GetUser() >= 0) continue;
+                               TObjArrayIter outers(inner->fMatches);
+                               while ( (outer = (AliITSNode*)outers.Next()) ) {
+                                       if (outer->GetUser() >= 0) continue;
+                                       if (!PassCurvCut(inner, outer, curv_idx, fVX, fVY, fVZ)) continue;
+                                       unit = new AliITSneuron;
+                                       unit->fInner = inner;
+                                       unit->fOuter = outer;
+                                       unit->fActivation = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
+                                       unit->fGain = new TObjArray;
+                                       fNeurons->AddLast(unit);
+                                       inner->fInnerOf->AddLast(unit);
+                                       outer->fOuterOf->AddLast(unit);
+                                       neurons++;
+                               } // for (;;)
+                       } // for (;;)
+               } // for (itheta...)
+       } // for (ilayer...)
+       */
+       
+       fNeurons->SetOwner();
+       return neurons;
+}
+//
+//
+//
+Int_t AliITSNeuralTracker::LinkNeurons()
+// Creates the neural synapses among all neurons
+// which share a point.
+{
+       Int_t total = 0;
+       TObjArrayIter iter(fNeurons), *test_iter;
+       AliITSneuron *neuron = 0, *test = 0;
+       for (;;) {
+               neuron = (AliITSneuron*)iter.Next();
+               if (!neuron) break;
+               // gain contributors
+               test_iter = (TObjArrayIter*)neuron->fInner->fOuterOf->MakeIterator();
+               for (;;) {
+                       test = (AliITSneuron*)test_iter->Next();
+                       if (!test) break;
+                       neuron->Add2Gain(test, fGain2CostRatio, fAlignExponent);
+                       total++;
+               }
+               delete test_iter;
+               test_iter = (TObjArrayIter*)neuron->fOuter->fInnerOf->MakeIterator();
+               for (;;) {
+                       test = (AliITSneuron*)test_iter->Next();
+                       if (!test) break;
+                       neuron->Add2Gain(test, fGain2CostRatio, fAlignExponent);
+                       total++;
+               }
+               delete test_iter;
+       }
+       return total;
+}
+//
+//
+//
+Bool_t AliITSNeuralTracker::Update()
+// A complete updating cycle with the asynchronous method
+// when the neural network is stable, kTRUE is returned.
+{
+       Double_t actVar = 0.0, totDiff = 0.0;
+       TObjArrayIter iter(fNeurons);
+       AliITSneuron *unit;
+       for (;;) {
+               unit = (AliITSneuron*)iter.Next();
+               if (!unit) break;
+               actVar = Activate(unit);
+               // calculation the relative activation variation
+               totDiff += actVar;
+       }
+       totDiff /= fNeurons->GetSize();
+       return (totDiff < fStabThreshold);
+}
+//
+//
+//
+void AliITSNeuralTracker::CleanNetwork()
+// Removes all deactivated neurons, and resolves the competitions
+// to the advantage of the most active unit
+{
+       AliITSneuron *unit = 0;
+       TObjArrayIter neurons(fNeurons);
+       while ( (unit = (AliITSneuron*)neurons.Next()) ) {
+               if (unit->fActivation < fActMinimum) {
+                       unit->fInner->fInnerOf->Remove(unit);
+                       unit->fOuter->fOuterOf->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->fInner->fInnerOf->GetSize();
+               nOut = (Int_t)unit->fOuter->fOuterOf->GetSize();
+               if (nIn < 2 && nOut < 2) continue;
+               removed = kFALSE;
+               if (nIn > 1) {
+                       TObjArrayIter competing(unit->fInner->fInnerOf);
+                       while ( (enemy = (AliITSneuron*)competing.Next()) ) {
+                               if (unit->fActivation > enemy->fActivation) {
+                                       enemy->fInner->fInnerOf->Remove(enemy);
+                                       enemy->fOuter->fOuterOf->Remove(enemy);
+                                       delete fNeurons->Remove(enemy);
+                               }
+                               else {
+                                       unit->fInner->fInnerOf->Remove(unit);
+                                       unit->fOuter->fOuterOf->Remove(unit);
+                                       delete fNeurons->Remove(unit);
+                                       removed = kTRUE;
+                                       break;
+                               }
+                       }
+                       if (removed) continue;
+               }
+               if (nOut > 1) {
+                       TObjArrayIter competing(unit->fOuter->fOuterOf);
+                       while ( (enemy = (AliITSneuron*)competing.Next()) ) {
+                               if (unit->fActivation > enemy->fActivation) {
+                                       enemy->fInner->fInnerOf->Remove(enemy);
+                                       enemy->fOuter->fOuterOf->Remove(enemy);
+                                       delete fNeurons->Remove(enemy);
+                               }
+                               else {
+                                       unit->fInner->fInnerOf->Remove(unit);
+                                       unit->fOuter->fOuterOf->Remove(unit);
+                                       delete fNeurons->Remove(unit);
+                                       removed = kTRUE;
+                                       break;
+                               }
+                       }
+               }
+       }
+}
+//
+//
+//
+Bool_t AliITSNeuralTracker::CheckOccupation()
+{
+       Int_t i;
+       for (i = 0; i < 6; i++) if (fPoint[i] < 0) return kFALSE;
+       return kTRUE;
+}
+//
+//
+//
+Int_t AliITSNeuralTracker::Save(Int_t sector_id)
+// Creates chains of active neurons, in order to
+// find the tracks obtained as the neural network output.
+{
+       // every chain is started from the neurons in the first 2 layers
+       Int_t i, check, stored = 0;
+       Double_t test_act = 0;
+       AliITSneuron *unit = 0, *cursor = 0, *fwd = 0;
+       AliITSNode *node = 0;
+       TObjArrayIter iter(fNeurons), *fwd_iter;
+       TObjArray *removedUnits = new TObjArray(0);
+       TObjArray *removedPoints = new TObjArray(6);
+       removedUnits->SetOwner(kFALSE);
+       removedPoints->SetOwner(kFALSE);
+       for (;;) {
+               for (i = 0; i < 6; i++) fPoint[i] = -1;
+               unit = (AliITSneuron*)iter.Next();
+               if (!unit) break;
+               if (unit->fInner->GetLayer() > 0) continue;
+               fPoint[unit->fInner->GetLayer()] = unit->fInner->fPosInTree;
+               fPoint[unit->fOuter->GetLayer()] = unit->fOuter->fPosInTree;
+               node = unit->fOuter;
+               removedUnits->AddLast(unit);
+               removedPoints->AddAt(unit->fInner, unit->fInner->GetLayer());
+               removedPoints->AddAt(unit->fOuter, unit->fOuter->GetLayer());
+               while (node) {
+                       test_act = fActMinimum;
+                       fwd_iter = (TObjArrayIter*)node->fInnerOf->MakeIterator();
+                       fwd = 0;
+                       for (;;) {
+                               cursor = (AliITSneuron*)fwd_iter->Next();
+                               if (!cursor) break;
+                               if (cursor->fUsed) continue;
+                               if (cursor->fActivation >= test_act) {
+                                       test_act = cursor->fActivation;
+                                       fwd = cursor;
+                               }
+                       }
+                       if (!fwd) break;
+                       removedUnits->AddLast(fwd);
+                       node = fwd->fOuter;
+                       fPoint[node->GetLayer()] = node->fPosInTree;
+                       removedPoints->AddAt(node, node->GetLayer());
+               }
+               check = 0;
+               for (i = 0; i < 6; i++) if (fPoint[i] >= 0) check++;
+               if (check >= 6) {
+                       stored++;
+                       fChains->Fill();
+                       for (i = 0; i < 6; i++) {
+                               node = (AliITSNode*)removedPoints->At(i);
+                               node->SetUser(1);
+                       }
+                       fwd_iter = (TObjArrayIter*)removedUnits->MakeIterator();
+                       for (;;) {
+                               cursor = (AliITSneuron*)fwd_iter->Next();
+                               if(!cursor) break;
+                               cursor->fUsed = 1;
+                       }
+               }
+       }
+
+       return stored;
+}
+/*
+Int_t AliITSNeuralTracker::Save(Int_t sector_idx)
+// Creates chains of active neurons, in order to
+// find the tracks obtained as the neural network output.
+{
+       // Reads the final map of neurons and removes all
+       // units with too small activation
+       cout << "saving..." << flush;
+       
+       // every chain is started from the neurons in the first 2 layers
+       //                     00111111  00011111  00101111  00110111
+       Int_t allow_mask[] = {     0x3F,     0x1F,     0x2F,     0x37,
+                                                                               0x3B,     0x3D,     0x3E };
+       //                     00111011  00111101  00111110
+
+       Double_t test_fwd = 0., test_back = 0.;
+       Int_t ilayer, itheta;
+       AliITSNode *node = 0;
+       AliITSneuron *fwd = 0, *back = 0;
+       TList *list_sector = 0;
+
+       cout << "A -" << fActMinimum << "-" << flush;
+       for (ilayer = 0; ilayer < 6; ilayer++) {
+               for (itheta = 0; itheta < 180; itheta++) {
+                       list_sector = (TList*)fPoints[ilayer][itheta]->At(sector_idx);
+                       TListIter iter(list_sector);
+                       while ( (node = (AliITSNode*)iter.Next()) ) {
+                               TListIter fwd_iter(node->fInnerOf);
+                               TListIter back_iter(node->fOuterOf);
+                               test_fwd = test_back = fActMinimum;
+                               while ( (fwd = (AliITSneuron*)fwd_iter.Next()) ) {
+                                       if (fwd->fActivation > test_fwd) {
+                                               test_fwd = fwd->fActivation;
+                                               node->fNext = fwd->fOuter;
+                                       }
+                               }
+                               while ( (back = (AliITSneuron*)back_iter.Next()) ) {
+                                       if (back->fActivation > test_back) {
+                                               test_back = back->fActivation;
+                                               node->fPrev = back->fInner;
+                                       }
+                               }
+                       }
+               }
+       }
+
+       cout << "B" << flush;
+       for (ilayer = 0; ilayer < 5; ilayer++) {
+               for (itheta = 0; itheta < 180; itheta++) {
+                       list_sector = (TList*)fPoints[ilayer][itheta]->At(sector_idx);
+                       TListIter iter(list_sector);
+                       while ( (node = (AliITSNode*)iter.Next()) ) {
+                               if (node->fNext) {
+                                       if (node->fNext->fPrev != node) node->fNext = 0;
+                               }
+                       }
+               }
+       }
+
+       cout << "C" << flush;
+       Bool_t ok;
+       Int_t stored = 0, l1, l2, i, mask, im;
+       AliITSNode *temp = 0;
+       TList *list[2];
+       for (itheta = 0; itheta < 180; itheta++) {
+               list[0] = (TList*)fPoints[0][itheta]->At(sector_idx);
+               list[1] = (TList*)fPoints[1][itheta]->At(sector_idx);
+               for (i = 0; i < 2; i++) {
+                       TListIter iter(list[i]);
+                       while ( (node = (AliITSNode*)iter.Next()) ) {
+                               //cout << endl << node << flush;
+                               AliITSneuralChain *chain = new AliITSneuralChain;
+                               for (temp = node; temp; temp = temp->fNext) {
+                                       chain->Insert((AliITSNeuralPoint*)temp);
+                               }
+                               ok = kFALSE;
+                               mask = chain->OccupationMask();
+                               for (im = 0; im < 7; im++) ok = ok || (mask == allow_mask[im]);
+                               if (!ok) {
+                                       delete chain;
+                                       continue;
+                               }
+                               //cout << " lab " << flush;
+                               chain->AssignLabel();
+                               //cout << " add " << flush;
+                               fTracks->AddLast(chain);
+                               stored++;
+                               //cout << " " << stored << flush;
+                       }
+               }
+       }
+
+       cout << "D" << flush;
+       return stored;
+}
+*/
+//
+//
+//
+Double_t AliITSNeuralTracker::Activate(AliITSneuron* &unit)
+// Computes the new neural activation and
+// returns the variation w.r.t the previous one
+{
+       if (!unit) return 0.0;
+       Double_t sum_gain = 0.0, sum_cost = 0.0, input, actOld, actNew;
+
+       // sum gain contributions
+       TObjArrayIter *iter = (TObjArrayIter*)unit->fGain->MakeIterator();
+       AliITSlink *link;
+       for(;;) {
+               link = (AliITSlink*)iter->Next();
+               if (!link) break;
+               sum_gain += link->fLinked->fActivation * link->fWeight;
+       }
+
+       // sum cost contributions
+       TObjArrayIter *test_iter = (TObjArrayIter*)unit->fInner->fInnerOf->MakeIterator();
+       AliITSneuron *linked = 0;
+       for (;;) {
+               linked = (AliITSneuron*)test_iter->Next();
+               if (!linked) break;
+               if (linked == unit) continue;
+               sum_cost += linked->fActivation;
+       }
+       delete test_iter;
+       test_iter = (TObjArrayIter*)unit->fOuter->fOuterOf->MakeIterator();
+       for (;;) {
+               linked = (AliITSneuron*)test_iter->Next();
+               if (!linked) break;
+               if (linked == unit) continue;
+               sum_cost += linked->fActivation;
+       }
+
+       //cout << "gain = " << sum_gain << ", cost = " << sum_cost << endl;
+
+       input = (sum_gain - sum_cost) / fTemperature;
+       actOld = unit->fActivation;
+#ifdef NEURAL_LINEAR
+       if (input <= -2.0 * fTemperature)
+               actNew = 0.0;
+       else if (input >= 2.0 * fTemperature)
+               actNew = 1.0;
+       else
+               actNew = input / (4.0 * fTemperature) + 0.5;
+#else
+       actNew = 1.0 / (1.0 + TMath::Exp(-input));
+#endif
+       unit->fActivation = actNew;
+       return TMath::Abs(actNew - actOld);
+}
+//
+//
+//
+void AliITSNeuralTracker::AliITSneuron::Add2Gain(AliITSneuron *n, Double_t mult_const, Double_t exponent)
+{
+       AliITSlink *link = new AliITSlink;
+       link->fLinked = n;
+       Double_t weight = Weight(n);
+       link->fWeight = mult_const * TMath::Power(weight, exponent);
+       fGain->AddLast(link);
+}
+//
+//
+//
+Double_t AliITSNeuralTracker::AliITSneuron::Weight(AliITSneuron *n)
+{
+       TVector3 me(fOuter->X() - fInner->X(), fOuter->Y() - fInner->Y(), fOuter->Z() - fInner->Z());
+       TVector3 it(n->fOuter->X() - n->fInner->X(), n->fOuter->Y() - n->fInner->Y(), n->fOuter->Z() - n->fInner->Z());
+
+       Double_t angle = me.Angle(it);
+       Double_t weight = 1.0 - sin(angle);
+       return weight;
+}
diff --git a/ITS/AliITSNeuralTracker.h b/ITS/AliITSNeuralTracker.h
new file mode 100644 (file)
index 0000000..1d5aa96
--- /dev/null
@@ -0,0 +1,199 @@
+#ifndef ALIITSNNEURALTRACKER_H
+#define ALIITSNNEURALTRACKER_H
+
+class TObjArray;
+class TCanvas;
+
+///////////////////////////////////////////////////////////////////////
+//
+// AliITSneuralTracker:
+//
+// neural network MFT algorithm
+// for track finding in ITS stand alone
+// according to the Denby-Peterson model with adaptments to the
+// ALICE multiplicity
+//
+///////////////////////////////////////////////////////////////////////
+
+class AliITSNeuralTracker : public TObject {
+
+public:
+
+                AliITSNeuralTracker();
+       virtual ~AliITSNeuralTracker();
+
+       // ******************************************************************************
+       // * Embedded utility class --> >>> NODE <<<
+       // ******************************************************************************
+       // * This class inherits from AliITSNeuralPoint and adds some
+       // * utility pointers for quick path-finding among neurons.
+       // ******************************************************************************
+       class AliITSNode : public AliITSNeuralPoint {
+       public:
+               AliITSNode() 
+                       {fInnerOf = fOuterOf = fMatches = 0; fNext = fPrev = 0;}
+                       
+               AliITSNode(AliITSNeuralPoint *p, Bool_t init = kTRUE) // declared inline
+                       : AliITSNeuralPoint(p)
+                       {
+                               fInnerOf = fOuterOf = fMatches = 0;
+                               fNext = fPrev = 0;
+                               if (init) {
+                                       fInnerOf = new TObjArray;
+                                       fOuterOf = new TObjArray;
+                                       fMatches = new TObjArray;
+                               }
+                       }
+                       
+               AliITSNode(AliITSRecPoint *p, AliITSgeomMatrix *gm)
+                       : AliITSNeuralPoint(p,gm) 
+                       {fInnerOf = fOuterOf = fMatches = 0; fNext = fPrev = 0;}
+
+               virtual  ~AliITSNode() 
+                       {fInnerOf = fOuterOf = fMatches = 0; fNext = fPrev = 0;}
+
+               Double_t  ThetaDeg()                    {return GetTheta()*180.0/TMath::Pi();}
+
+               Int_t     GetSector(Double_t sec_width) {return (Int_t)(GetPhi()/sec_width);}
+               Int_t     GetThetaCell()                {return (Int_t)(ThetaDeg());}
+               
+               Int_t        fPosInTree;
+
+               TObjArray   *fInnerOf; //!
+               TObjArray   *fOuterOf; //! 
+               TObjArray   *fMatches; //!
+
+               AliITSNode  *fNext; //!
+               AliITSNode  *fPrev; //!
+       };
+       // ******************************************************************************
+
+
+
+       // ******************************************************************************
+       // * Embedded utility class --> >>> NEURON <<<
+       // ******************************************************************************
+       // * Simple class implementing the neural unit.
+       // * Contains the activation and some other pointers
+       // * for purposes similar to the ones in AliITSnode.
+       // ******************************************************************************
+       class AliITSneuron : public TObject {
+       public:
+                           AliITSneuron():fUsed(0),fActivation(0.),fInner(0),fOuter(0),fGain(0) { }
+               virtual    ~AliITSneuron() {fInner=fOuter=0;fGain=0;}
+
+               Double_t    Weight(AliITSneuron *n);
+               void        Add2Gain(AliITSneuron *n, Double_t mult_const, Double_t exponent);
+
+               Int_t             fUsed;        //  utility flag
+               Double_t          fActivation;  //  Activation value
+               AliITSNode       *fInner;       //! inner point
+               AliITSNode       *fOuter;       //! outer point
+               TObjArray        *fGain;        //! list of sequenced units
+       };
+       // ******************************************************************************
+
+
+
+       // ******************************************************************************
+       // * Embedded utility class --> >>> CONNECTION <<<
+       // ******************************************************************************
+       // * Used to implement the neural weighted connection
+       // * in such a way to speed up the retrieval of the
+       // * links among neuron, for a fast update procedure.
+       // ******************************************************************************
+       class AliITSlink : public TObject {
+       public:
+                        AliITSlink() : fWeight(0.), fLinked(0) { }
+               virtual ~AliITSlink()   {fLinked = 0;}
+
+               Double_t      fWeight;  //  Weight value
+               AliITSneuron *fLinked;  //! the connected neuron
+       };
+       // ******************************************************************************
+
+
+       // Cut related setters
+
+       void     SetHelixMatchCuts(Double_t *min, Double_t *max);
+       void     SetThetaCuts2D(Double_t *min, Double_t *max);
+       void     SetThetaCuts3D(Double_t *min, Double_t *max);
+       void     SetCurvatureCuts(Int_t n, Double_t *cuts);
+       void     SetVertex(Double_t x, Double_t y, Double_t z)  {fVX=x; fVY=y; fVZ=z;}
+       void     SetPolarInterval(Double_t dtheta) {fPolarInterval=dtheta;}
+
+       // Neural work-flow related setters
+
+       void     SetActThreshold(Double_t val)            {fActMinimum = val;}
+       void     SetWeightExponent(Double_t a)            {fAlignExponent = 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;}
+
+       // Points array arrangement & control
+
+       void     CreateArrayStructure(Int_t nsecs);
+       Int_t    ArrangePoints(TTree *pts_tree);
+       void     StoreAbsoluteMatches();
+       Bool_t   PassCurvCut(AliITSNode *p1, AliITSNode *p2, Int_t curv_idx, Double_t vx, Double_t vy, Double_t vz);
+       Int_t    PassAllCuts(AliITSNode *p1, AliITSNode *p2, Int_t curv_idx, Double_t vx, Double_t vy, Double_t vz);
+       void     PrintPoints();
+       void     PrintMatches(Bool_t stop = kTRUE);
+
+       // Neural tracker work-flow
+
+       void     NeuralTracking(const char* filesave, TCanvas*& display);
+       void     Display(TCanvas*& canvas);
+       void     ResetNodes(Int_t isector);
+       Int_t    CreateNeurons(Int_t sector, Int_t curv);  // create neurons
+       Int_t    LinkNeurons();                // create neural connections
+       Double_t Activate(AliITSneuron* &n);   // calculates the new neural activation
+       Bool_t   Update();                     // an updating cycle
+       void     CleanNetwork();               // removes deactivated units and resolves competitions
+       Int_t    Save(Int_t sector_idx);       // save found tracks for # sector
+       TTree*   GetChains()                   {return fChains;}
+       void     WriteChains()                 {fChains->Write();}
+
+private:
+
+       Bool_t       CheckOccupation(); 
+
+       Int_t        fSectorNum;            //  number of azymuthal sectors
+       Double_t     fSectorWidth;          //  width of an azymuthal sector (in RADIANS) [used internally]
+       Double_t     fPolarInterval;        //  width of a polar sector (in DEGREES)
+
+       Double_t     fThetaCut2DMin[5];     //  lower edge of theta cut range (in DEGREES)
+       Double_t     fThetaCut2DMax[5];     //  upper edge of theta cut range (in DEGREES)
+       Double_t     fThetaCut3DMin[5];     //  lower edge of theta cut range (in DEGREES)
+       Double_t     fThetaCut3DMax[5];     //  upper edge of theta cut range (in DEGREES)
+       Double_t     fHelixMatchCutMin[5];  //  lower edge of helix matching cut range
+       Double_t     fHelixMatchCutMax[5];  //  lower edge of helix matching cut range
+       Int_t        fCurvNum;              //  # of curvature cut steps
+       Double_t    *fCurvCut;              //! value of all curvature cuts
+
+       Bool_t       fStructureOK;          // flag to avoid MANY segfault errors
+
+       Double_t     fVX, fVY, fVZ;         //  estimated vertex coords (for helix matching cut)
+
+       Double_t     fActMinimum;           //  minimum activation to turn 'on' the unit at the end
+       Double_t     fEdge1, fEdge2;        //  initialization interval for activations
+
+       Double_t     fTemperature;          //  logistic function temperature parameter
+       Double_t     fStabThreshold;        //  stability threshold
+       Double_t     fGain2CostRatio;       //  ratio between inhibitory and excitory contributions
+       Double_t     fAlignExponent;        //  alignment-dependent weight term
+
+       Int_t        fPoint[6];             //  Track point in layers
+       TTree       *fChains;               //! Output tree
+
+       TObjArray   *fPoints[6][180];       //! recpoints arranged into sectors for processing
+       TObjArray   *fNeurons;              //! neurons
+
+       ClassDef(AliITSNeuralTracker, 1)
+};
+
+
+////////////////////////////////////////////////////////////////////////////////
+
+#endif