--- /dev/null
+#include <stdlib.h>
+#include <iostream.h>
+
+#include <TROOT.h>
+#include <TObject.h>
+#include <TMath.h>
+#include <TString.h>
+
+#include "AliITSglobalRecPoint.h"
+
+
+ClassImp(AliITSglobalRecPoint)
+//
+//
+//
+AliITSglobalRecPoint::AliITSglobalRecPoint()
+{
+ // Default constructor:
+ // sets to zero all members
+ fLayer = 0;
+
+ fGX = fGSX = 0.0;
+ fGY = fGSY = 0.0;
+ fGZ = fGSZ = 0.0;
+
+ fR2 = fR3 = fPhi = fTheta = 0.0;
+
+ fLabel[0] = fLabel[1] = fLabel[2] = -1;
+ fKalmanLabel = 0;
+ fSign = 0;
+ fModule = 0;
+ fPosInModule = 0;
+
+ fUsed = 0;
+}
+//
+//
+//
+AliITSglobalRecPoint::AliITSglobalRecPoint(Double_t gx, Double_t gy,
+ Double_t gz, Double_t gsx, Double_t gsy, Double_t gsz, Int_t l) :
+ fGX(gx), fGY(gy), fGZ(gz), fGSX(gsx), fGSY(gsy), fGSZ(gsz), fLayer(l)
+{
+ fR2 = TMath::Sqrt(fGX*fGX + fGY*fGY);
+ fR3 = TMath::Sqrt(fGX*fGX + fGY*fGY + fGZ*fGZ);
+ fPhi = TMath::ATan2(fGY, fGX);
+ if (fPhi < 0.0) fPhi = 2.0 * TMath::Pi() + fPhi;
+ fTheta = TMath::ATan2(fR2, fGZ);
+
+ fLabel[0] = fLabel[1] = fLabel[2] = -1;
+ fUsed = 0;
+ fKalmanLabel = 0;
+ fSign = 0;
+ fModule = 0;
+ fPosInModule = 0;
+}
+//
+//
+/*
+Double_t AliITSglobalRecPoint::GetGlobalSigmaR2()
+{
+ // Square sigma of the 2-D radius:
+ // = X^2 * Sx^2 + Y^2 * Sy^2
+
+ Double_t answer = fGX*fGX*fGSigmaX + fGY*fGY*fGSigmaY;
+ return answer / (fR2 * fR2);
+}
+//
+//
+//
+Double_t AliITSglobalRecPoint::GetGlobalSigmaR3()
+{
+ // Square sigma of the 3-D radius:
+ // = (X^2 * Sx^2 + Y^2 * Sy^2 + Z^2 * Sz^2) / R^2
+ // R in 3-D, r in 2-D
+
+ Double_t answer = fGX*fGX*fGSigmaX + fGY*fGY*fGSigmaY + fGZ*fGZ*fGSigmaZ;
+ return answer / (fR3 * fR3);
+}
+//
+//
+//
+Double_t AliITSglobalRecPoint::GetGlobalSigmaTheta()
+{
+ // Square sigma of theta:
+ // = (Z^2 * (X^2 * Sx^2 + Y^2 * Sy^2) + r^4 * Sz^2) / (R^4 * r^2)
+ // R in 3-D, r in 2-D
+
+ Double_t answer = fGZ*fGZ*(fGX*fGX*fGSigmaX + fGY*fGY*fGSigmaY) + fR2*fR2*fR2*fR2*fGSigmaZ;
+ return answer / (fR3*fR3*fR3*fR3*fR2*fR2);
+}
+//
+//
+//
+Double_t AliITSglobalRecPoint::GetGlobalSigmaPhi()
+{
+ // Square sigma of phi:
+ // = (Y^2 * Sx^2 + X^2 * Sy^2) / r^4
+ // R in 3-D, r in 2-D
+
+ Double_t answer = fGY*fGY*fGSigmaX + fGX*fGX*fGSigmaY;
+ return answer / (fR2 * fR2 * fR2 * fR2);
+}
+//
+//
+//
+Bool_t AliITSglobalRecPoint::SharesID(AliITSRecPoint* pt)
+{
+ // Controls if there is a track index shared by two points
+
+ Bool_t ok = kFALSE;
+ Int_t i, j;
+ for (i = 0; i < 3; i++) {
+ if (fTracks[i] < 0) continue;
+ for (j = 0; j < 3; j++) {
+ if (pt->fTracks[j] < 0) continue;
+ ok = ok || (fTracks[i] == pt->fTracks[j]);
+ }
+ }
+ return ok;
+}
+*/
+//
+//
+Double_t AliITSglobalRecPoint::DPhi(AliITSglobalRecPoint *p)
+{
+ // Absolute value of the difference
+ // between 'phi' coordinates of two points
+ // controlled in order to avoid that, for
+ // reasons of initial values, it come > 180. degrees
+
+ Double_t phi = TMath::Abs(fPhi - p->fPhi);
+ if (phi > TMath::Pi())
+ phi = 2.0 * TMath::Pi() - phi;
+ return phi;
+}
+//
+//
+//
+Double_t AliITSglobalRecPoint::DTheta(AliITSglobalRecPoint *p)
+{
+ // Absolute value of the difference
+ // between 'theta' coordinates of two points
+
+ return TMath::Abs(fTheta - p->fTheta);
+}
+//
+//
+//
+Int_t AliITSglobalRecPoint::Compare(const TObject *O) const
+{
+ // Criterion for sorting:
+ // p1 < p2 if p1.phi < p2.phi;
+
+ // Casting to convert the given argument to rec-point
+ AliITSglobalRecPoint *you = (AliITSglobalRecPoint*)O;
+
+ // Comparation
+ if (fLayer < you->fLayer)
+ return -1;
+ else if (fLayer > you->fLayer)
+ return 1;
+ else {
+ if (fPhi < you->fPhi)
+ return -1;
+ else if (fPhi > you->fPhi)
+ return 1;
+ else
+ return 0;
+ }
+}
+//
+//
+//
--- /dev/null
+#ifndef ALIITSGLOBALRECPOINT_H
+#define ALIITSGLOBALRECPOINT_H
+
+class AliITSglobalRecPoint : public TObject {
+
+ public:
+
+ AliITSglobalRecPoint();
+ AliITSglobalRecPoint(Double_t gx, Double_t gy, Double_t gz, Double_t gsx, Double_t gsy, Double_t gsz, Int_t l);
+
+ virtual ~AliITSglobalRecPoint() { }
+
+ // difference in azymuthal angle
+ Double_t DPhi (AliITSglobalRecPoint *p);
+ // difference in polar angle
+ Double_t DTheta(AliITSglobalRecPoint *p);
+ // checks if the pt contains a given positive label
+ Bool_t HasID (Int_t ID)
+ { if (ID < 0) return kFALSE; else return (fLabel[0]==ID || fLabel[1]==ID || fLabel[2]==ID); }
+ // checks if the pt shares at least a positive label with another one
+ Bool_t SharesID(AliITSglobalRecPoint *p)
+ { return (HasID(p->fLabel[0]) || HasID(p->fLabel[1]) || HasID(p->fLabel[2])); }
+
+ // Parameters for sorting
+ Bool_t IsSortable() const { return kTRUE; }
+ Int_t Compare(const TObject *O) const;
+
+ public:
+
+ Double_t fGX, fGY, fGZ; // (x,y,z) in the global reference
+
+ Double_t fR2; // = sqrt(x^2 + y^2)
+ Double_t fR3; // = sqrt(x^2 + y^2 + z^2)
+ Double_t fPhi; // = atan(y/x)
+ Double_t fTheta; // = acos(z/r3)
+
+ Double_t fGSX; //
+ Double_t fGSY; // sigmas of global coords
+ Double_t fGSZ; //
+
+ Int_t fModule; // ITS module containing the point
+ Int_t fPosInModule; // position in TClonesArray of TreeR
+ Int_t fLayer; // ITS layer containing the point
+ Int_t fLabel[3]; // Generated tracks containing the point
+ Int_t fKalmanLabel; // Kalman recognized track containing the point
+
+ Int_t fSign; // 0 if the point doesn't belong to a Kalman track
+ // + 1 if the point belongs to a Kalman good track
+ // - 1 if the point belongs to a Kalman fake track
+
+ Int_t fUsed; // a flag to avoid point recycling
+
+ ClassDef(AliITSglobalRecPoint, 1) // AliITSglobalRecPoints class
+};
+
+#endif
--- /dev/null
+#include <fstream.h>
+#include <stdlib.h>
+
+#include <TObject.h>
+#include <TROOT.h>
+#include <TMath.h>
+#include <TString.h>
+
+#include "AliITSglobalRecPoint.h"
+#include "AliITSneuralTracker.h"
+
+#include "AliITSneuralTrack.h"
+
+
+
+ClassImp(AliITSneuralTrack)
+//
+//
+//
+AliITSneuralTrack::AliITSneuralTrack()
+{
+ Int_t i;
+ for (i = 0; i < 6; i++) fPoint[i] = 0;
+}
+//
+//
+//
+AliITSneuralTrack::~AliITSneuralTrack()
+{
+ Int_t i;
+ for (i = 0; i < 6; i++) delete fPoint[i];
+}
+//
+//
+//
+Int_t AliITSneuralTrack::CheckMe(Bool_t verbose)
+{
+ Int_t l, stored = 0;
+ TString empty("Not filled slots: ");
+ for (l = 0; l < 6; l++) {
+ if (fPoint[l])
+ stored++;
+ else {
+ empty += l;
+ empty += ' ';
+ }
+ }
+
+ if (stored < 6 && verbose) Warning("", empty);
+ return stored;
+}
+//
+//
+//
+Int_t AliITSneuralTrack::EvaluateTrack(Bool_t verbose, Int_t min, Int_t* &good)
+{
+ Int_t i, j, k = 0, count[18], id[18];
+ for (i = 0; i < 6; i++) {
+ for (j = 0; j < 3; j++) {
+ if (fPoint[i])
+ id[k] = fPoint[i]->fLabel[j];
+ else
+ id[k] = -1;
+ count[k] = 0;
+ k++;
+ }
+ }
+ for (i = 0; i < 18; i++) {
+ for (j = i+1; j < 18; j++) {
+ if (id[i] == id[j] || id[j] < 0) id[j] = -1;
+ }
+ }
+ for (i = 0; i < 18; i++) {
+ if (id[i] < 0) continue;
+ for (j = 0; j < 6; j++) {
+ if (fPoint[j] && fPoint[j]->HasID(id[i])) count[i]++;
+ }
+ }
+ Int_t index[18], best;
+ TMath::Sort(18, count, index);
+ best = id[index[0]];
+ if (verbose) {
+ if (count[index[0]]<min) cout << "\t";
+ cout << best << " " << count[index[0]] << endl;
+ }
+
+ if (good) delete [] good;
+ good = new Int_t[6];
+ for (i = 0; i < 6; i++) {
+ good[i] = 0;
+ if (!fPoint[i])
+ good[i] = -1;
+ else if (fPoint[i] && fPoint[i]->HasID(best))
+ good[i] = 1;
+ else if (fPoint[i]->fLabel[0] < 0 && fPoint[i]->fLabel[1] < 0 && fPoint[i]->fLabel[2] < 0)
+ good[i] = -1;
+ }
+
+ if (count[index[0]] < min) best = -best;
+
+ return best;
+
+
+
+ /*
+ if (good) delete [] good;
+ good = new Int_t[6];
+ Int_t count = CheckMe(verbose);
+
+ if (count < min) return -9999999;
+
+ Int_t i, l, max = 0, best = 0;
+ for (l = 0; l < 6; l++) {
+ for (i = 0; i < 3; i++) {
+ if (fPoint[l].fLabel[i] > max) max = fPoint[l].fLabel[i];
+ }
+ }
+
+ count = max + 1;
+ Int_t *counter = new Int_t[count];
+ for (i = 0; i < count; i++) counter[i] = 0;
+
+ for (l = 0; l < 6; l++) {
+ for (i = 0; i < 3; i++) {
+ if (fPoint[l].fLabel[i] >= 0)
+ counter[fPoint[l].fLabel[i]]++;
+ }
+ }
+
+ for (i = 0; i < count; i++) {
+ if (counter[i] > counter[best]) best = i;
+ }
+
+ for (l = 0; l < 6; l++) {
+ good[l] = 0;
+ if (fPoint[l].fLabel[0] == best || fPoint[l].fLabel[1] == best || fPoint[l].fLabel[2] == best)
+ good[l] = 1;
+ else if (fPoint[l].fLabel[0] < 0 && fPoint[l].fLabel[1] < 0 && fPoint[l].fLabel[2] < 0)
+ good[l] = -1;
+ }
+
+ if (counter[best] < min) best = -best;
+ delete counter;
+ return best;
+ */
+}
+//
+//
+//
+void AliITSneuralTrack::GetCoords(Double_t* &x, Double_t* &y, Double_t* &z)
+{
+ if (x) delete [] x; x = new Double_t[6];
+ if (y) delete [] y; y = new Double_t[6];
+ if (z) delete [] z; z = new Double_t[6];
+
+ Int_t i;
+
+ for (i = 0; i < 6; i++) {
+ x[i] = y[i] = z[i] = 0.0;
+ if (!fPoint[i]) continue;
+ x[i] = fPoint[i]->fGX;
+ y[i] = fPoint[i]->fGY;
+ z[i] = fPoint[i]->fGZ;
+ }
+}
+//
+//
+//
+void AliITSneuralTrack::CopyPoint(AliITSglobalRecPoint *p)
+{
+ Int_t layer = p->fLayer;
+ if (layer < 0 || layer > 6) {
+ Error("", Form("Wrong layer [%d]", layer));
+ return;
+ }
+
+ fPoint[layer] = new AliITSglobalRecPoint;
+ fPoint[layer]->fGX = p->fGX;
+ fPoint[layer]->fGY = p->fGY;
+ fPoint[layer]->fGZ = p->fGZ;
+ fPoint[layer]->fGSX = p->fGSX;
+ fPoint[layer]->fGSY = p->fGSY;
+ fPoint[layer]->fGSZ = p->fGSZ;
+ fPoint[layer]->fLayer = layer;
+ fPoint[layer]->fLabel[0] = p->fLabel[0];
+ fPoint[layer]->fLabel[1] = p->fLabel[1];
+ fPoint[layer]->fLabel[2] = p->fLabel[2];
+}
+//
+//
+//
+void AliITSneuralTrack::Print(Option_t *option, Int_t min)
+{
+ Int_t *vuoto = 0;
+ TString opt(option);
+ opt.ToUpper();
+ Int_t id = EvaluateTrack(0, min, vuoto);
+ if (opt.Contains("A")) {
+ cout << "\nEvaluated ID for this track: ";
+ cout.width(8);
+ if (id >= 0) {
+ cout << id;
+ cout << " [good]";
+ }
+ else {
+ cout << -id;
+ cout << " [fake]";
+ }
+ }
+ cout << endl << endl;
+}
+//
+//
+//
+void AliITSneuralTrack::Kinks(Int_t &pos, Int_t &neg, Int_t &incr, Int_t &decr)
+{
+ Int_t i;
+ Double_t dphi, dphi_old = 0.0;
+ pos = neg = incr = decr = 0;
+ for (i = 1; i < 6; i++) {
+ dphi = fPoint[i]->fPhi - fPoint[i-1]->fPhi;
+ if (dphi > 0.0) pos++; else neg++;
+ if (TMath::Abs(dphi) > dphi_old) incr++; else decr++;
+ dphi_old = TMath::Abs(dphi);
+ }
+}
+//
+//
+//
+Double_t AliITSneuralTrack::FitXY(Double_t VX, Double_t VY)
+{
+ Int_t i;
+ Double_t X, Y, D, R;
+ Double_t rx(0.0), ry(0.0), x2(0.0), y2(0.0), xy(0.0);
+ for (i = 0; i < 6; i++) {
+ X = fPoint[i]->fGX - VX;
+ Y = fPoint[i]->fGY - VY;
+ R = X * X + Y * Y;
+ rx += R * X;
+ ry += R * Y;
+ x2 += X * X;
+ y2 += Y * Y;
+ xy += X * Y;
+ }
+
+ D = 2 * (x2 * y2 - xy * xy);
+ if (D == 0.0)
+ return 1000.0;
+ else {
+ X = (rx * y2 - ry * xy) / D;
+ Y = (ry * x2 - rx * xy) / D;
+ fFitRadius = TMath::Sqrt(X * X + Y * Y);
+ fFitXC = X + VX;
+ fFitYC = Y + VY;
+ }
+
+ fSqChi = 0.0;
+ for (i = 0; i < 6; i++) {
+ X = fPoint[i]->fGX - fFitXC;
+ Y = fPoint[i]->fGY - fFitYC;
+ fSqChi += ((X * X + Y * Y) / (fFitRadius * fFitRadius)) * ((X * X + Y * Y) / (fFitRadius * fFitRadius));
+ }
+ fSqChi /= 6.0;
+ fSqChi = TMath::Sqrt(fSqChi);
+ return fSqChi;
+}
--- /dev/null
+#ifndef ALIITSNEURALTRACK_H
+#define ALIITSNEURALTRACK_H
+
+class AliITSglobalRecPoint;
+
+class AliITSneuralTrack : public TObject {
+
+ friend class AliITSneuralTracker;
+
+public:
+ AliITSneuralTrack();
+ virtual ~AliITSneuralTrack();
+
+ Double_t GetSqChi() {return fSqChi;}
+ Double_t GetR(Double_t &Xc, Double_t &Yc)
+ {Xc=fFitXC; Yc=fFitYC; return fFitRadius;}
+
+ Int_t CheckMe(Bool_t verbose);
+ Int_t EvaluateTrack(Bool_t verbose, Int_t min, Int_t* &good);
+
+ void GetCoords(Double_t* &x, Double_t* &y, Double_t* &z);
+ void CopyPoint(AliITSglobalRecPoint *p);
+ void Print(Option_t *option, Int_t min);
+ void Kinks(Int_t &pos, Int_t &neg, Int_t &decr, Int_t &incr);
+
+ AliITSglobalRecPoint* GetPoint(Int_t i) { return fPoint[i]; }
+
+private:
+
+ Double_t FitXY(Double_t VX, Double_t VY); //!
+
+ Double_t fFitXC;
+ Double_t fFitYC;
+ Double_t fFitRadius;
+ Double_t fFitTanL;
+ Double_t fSqChi;
+
+ AliITSglobalRecPoint *fPoint[6];
+
+ ClassDef(AliITSneuralTrack, 1)
+};
+
+#endif
--- /dev/null
+/**************************************************************************
+ * 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. *
+ * *
+ * track finder for ITS stand-alone with neural network algorithm *
+ * This class defines a Hopfield MFT neural network simulation *
+ * which reads all recpoints from an event and produces a tree with *
+ * the points belonging to recognized tracks *
+ * the TTree obtained as the output file will be saved in a .root file *
+ * and it must be read by groups of six entries at a time *
+ **************************************************************************/
+
+
+#include <fstream.h>
+#include <stdlib.h>
+
+#include <TROOT.h>
+#include <TFile.h>
+#include <TObjArray.h>
+#include <TTree.h>
+#include <TMath.h>
+#include <TRandom.h>
+#include <TVector3.h>
+
+#include "AliRun.h"
+#include "AliITS.h"
+#include "AliITSgeom.h"
+#include "AliITSgeomMatrix.h"
+#include "AliITSRecPoint.h"
+#include "AliITSglobalRecPoint.h"
+#include "AliITSneuralTrack.h"
+#include "AliITSneuralTracker.h"
+
+ClassImp(AliITSneuralTracker)
+
+
+// ====================================================================================================
+// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> AliITSneuralTracker <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+// ====================================================================================================
+
+AliITSneuralTracker::AliITSneuralTracker()
+{
+ //
+ // Default constructor (not for use)
+ //
+ Int_t i = 0;
+ fCurvCut = 0;
+ for (i = 0; i < 6; i++) {
+ fPoints[i] = 0;
+ if (i < 5) fNeurons[i] = 0;
+ }
+ fTracks = 0;
+}
+//---------------------------------------------------------------------------------------------------------
+AliITSneuralTracker::AliITSneuralTracker (Int_t nsecs, Int_t ncurv, Double_t *curv, Double_t theta)
+{
+ // Argumented constructor
+ // ----------------------
+ // gets the number of azymuthal sectors,
+ // the number of curvature cut steps and their cut values,
+ // and then the polar angle cut.
+ // Be careful not to put 'ncurv' greater than the dimension of the
+ // 'curv' array (--> segmentation fault)
+
+ fSectorNum = nsecs;
+ fSectorWidth = 2.0 * TMath::Pi()/(Double_t)nsecs;
+
+ fCurvNum = ncurv;
+ fCurvCut = new Double_t[ncurv];
+ Int_t i;
+ for (i = 0; i < ncurv; i++) fCurvCut[i] = curv[i];
+
+ fThetaCut = theta;
+ fThetaNum = (Int_t)(TMath::Pi() / theta) + 1;
+ if (fThetaNum < 1) fThetaNum = 1;
+
+ cout << "\nClass AliITSneuralTracker invoked with " << fSectorNum << " phi secs.";
+ cout << " and " << fThetaNum << " theta secs.\n" << endl;
+
+ Int_t k, j;
+ for (k = 0; k < 6; k++) {
+ fPoints[k] = new TObjArray**[fSectorNum];
+ for (i = 0; i < fSectorNum; i++) {
+ fPoints[k][i] = new TObjArray*[fThetaNum];
+ for (j = 0; j < fThetaNum; j++) fPoints[k][i][j] = new TObjArray(0);
+ }
+ if (k < 6) fNeurons[k] = new TObjArray(0);
+ }
+
+ fTracks = new TObjArray(1);
+}
+//---------------------------------------------------------------------------------------------------------
+AliITSneuralTracker::~AliITSneuralTracker()
+{
+ // Destructor
+ // ----------
+ // clears the TObjArrays and all heap objects
+
+ Int_t lay, sec, k;
+ delete [] fCurvCut;
+
+ for (lay = 0; lay < 6; lay++) {
+ for (sec = 0; sec < fSectorNum; sec++) {
+ for (k = 0; k < fThetaNum; k++) {
+ fPoints[lay][sec][k]->Delete();
+ }
+ }
+ delete [] fPoints[lay];
+ if (lay < 5) {
+ fNeurons[lay]->Delete();
+ delete fNeurons[lay];
+ }
+ }
+
+ fTracks->Delete();
+}
+//---------------------------------------------------------------------------------------------------------
+Int_t AliITSneuralTracker::ReadFile(const Text_t *fname, Int_t evnum)
+{
+ // File reader
+ // -----------
+ // Reads a ROOT file in order to retrieve recpoints.
+ // the 'evnum' argument can have two meanings:
+ // if it is >= 0 it represents the event number to retrieve the correct gAlice
+ // if it is < 0, instead, this means that the open file contains
+ // a 'TreeP' tree, with all points which remained unused
+ // after the Kalman tracking procedure (produced via the ITSafterKalman.C macro)
+ // the return value is the number of collected points
+ // (then, if errors occur, the method returns 0)
+
+ TFile *file = new TFile(fname, "read");
+ TTree *tree = 0;
+ Int_t i, nentries, total, sector, mesh;
+
+ if (evnum < 0) {
+ tree = (TTree*)file->Get("TreeP0");
+ if (!tree) {
+ Error("", "Specifying evnum < 0, I expect to find a 'TreeP' tree in the file, but there isn't any");
+ return 0;
+ }
+ AliITSglobalRecPoint *p = 0;
+ tree->SetBranchAddress("Points", &p);
+ nentries = (Int_t)tree->GetEntries();
+ total = 0;
+ for (i = 0; i < nentries; i++) {
+ tree->GetEntry(i);
+ AliITSglobalRecPoint *q = new AliITSglobalRecPoint(*p);
+ sector = (Int_t)(q->fPhi / fSectorWidth);
+ mesh = (Int_t)(q->fTheta / fThetaCut);
+ q->fUsed = 0;
+ fPoints[q->fLayer][sector][mesh]->AddLast(q);
+ total++;
+ }
+ return total;
+ }
+
+ if (gAlice) { delete gAlice; gAlice = 0; }
+ gAlice = (AliRun*)file->Get("gAlice");
+ if (!gAlice) {
+ Error("", "gAlice is NULL, maybe wrong filename...");
+ return 0;
+ }
+
+ Int_t nparticles = (Int_t)gAlice->GetEvent(evnum);
+ if (!nparticles) {
+ Error("", "No particles found");
+ return 0;
+ }
+
+ AliITS *its = (AliITS*)gAlice->GetModule("ITS");
+ if (!its) {
+ Error("", "No ITS found");
+ return 0;
+ }
+
+ AliITSgeom *geom = (AliITSgeom*)its->GetITSgeom();
+ if (!geom) {
+ Error("", "AliITSgeom is NULL");
+ return 0;
+ }
+
+ tree = gAlice->TreeR();
+ if (!tree) {
+ Error("", "TreeR() returned NULL");
+ return 0;
+ }
+
+ nentries = (Int_t)tree->GetEntries();
+ if (!nentries) {
+ Error("", "TreeR is empty");
+ return 0;
+ }
+
+ // Reading objects initialization
+ Int_t k, entry, lay, lad, det;
+ Double_t l[3], g[3], ll[3][3], gg[3][3];
+
+ TObjArray *recsArray = its->RecPoints();
+ AliITSRecPoint *recp;
+ AliITSglobalRecPoint *pnt = 0;
+
+ total = 0;
+ for(entry = 0; entry < nentries; entry++) {
+ if (entry > geom->GetLastSSD()) continue;
+ its->ResetRecPoints();
+ tree->GetEvent(entry);
+ TObjArrayIter recs(recsArray);
+ AliITSgeomMatrix *gm = geom->GetGeomMatrix(entry);
+ geom->GetModuleId(entry, lay, lad, det);
+ lay--;
+ while ((recp = (AliITSRecPoint*)recs.Next())) {
+ // conversion to global coordinate system
+ for (i = 0; i < 3; i++) {
+ l[i] = g[i] = 0.0;
+ for (k = 0; k < 3; k++) {
+ ll[i][k] = 0.0;
+ gg[i][k] = 0.0;
+ }
+ }
+ // local to global conversions of coords
+ l[0] = recp->fX;
+ l[1] = 0.0;
+ l[2] = recp->fZ;
+ gm->LtoGPosition(l, g);
+ // local to global conversions of sigmas
+ ll[0][0] = recp->fSigmaX2;
+ ll[2][2] = recp->fSigmaZ2;
+ gm->LtoGPositionError(ll, gg);
+ // instantiation of a new global rec-point
+ pnt = new AliITSglobalRecPoint(g[0], g[1], g[2], gg[0][0], gg[1][1], gg[2][2], lay);
+ // copy of owner track labels
+ for (i = 0; i < 3; i++) pnt->fLabel[i] = recp->fTracks[i];
+ sector = (Int_t)(pnt->fPhi / fSectorWidth);
+ mesh = (Int_t)(pnt->fTheta / fThetaCut);
+ pnt->fUsed = 0;
+ fPoints[lay][sector][mesh]->AddLast(pnt);
+ total++;
+ }
+ }
+
+ return total;
+}
+//---------------------------------------------------------------------------------------------------------
+void AliITSneuralTracker::Go(const char* filesave, Bool_t flag)
+{
+ // Global working method
+ // ---------------------
+ // It's the method which calls all other ones,
+ // in order to analyze the whole ITS points ensemble
+ // It's the only method to use, besides the parameter setters,
+ // so, all other methods are defined as private
+ // (see header file)
+ // the flag meaning is explained in the 'Initialize' method
+
+ Int_t i;
+ for (i = 0; i < fSectorNum; i++) DoSector(i, flag);
+ cout << endl << "Saving results in file " << filesave << "..." << flush;
+ TFile *f = new TFile(filesave, "recreate");
+ fTracks->Write();
+ cout << "done" << endl;
+ f->Close();
+}
+
+//=========================================================================================================
+// * * * P R I V A T E S E C T I O N * * *
+//=========================================================================================================
+
+Int_t AliITSneuralTracker::Initialize(Int_t secnum, Int_t curvnum, Bool_t flag)
+{
+ // Network initializer
+ // -------------------
+ // 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)
+
+ Int_t l, m, k, neurons = 0;
+
+ for (l = 0; l < 5; l++) fNeurons[l]->Delete();
+
+ AliITSneuron *unit = 0;
+ Double_t abscurv, max = fCurvCut[curvnum], min = 0.0;
+ if (flag && curvnum > 0) min = fCurvCut[curvnum - 1];
+ AliITSglobalRecPoint *inner = 0, *outer = 0;
+ for (l = 0; l < 5; l++) {
+ for (m = 0; m < fThetaNum; m++) {
+ TObjArrayIter inners(fPoints[l][secnum][m]);
+ while ((inner = (AliITSglobalRecPoint*)inners.Next())) {
+ if (inner->fUsed > 0) continue; // points can't be recycled
+ for (k = m-1; k <= m+1; k++) {
+ if (k < 0 || k >= fThetaNum) continue; // to avoid seg faults
+ TObjArrayIter outers(fPoints[l+1][secnum][k]);
+ while ((outer = (AliITSglobalRecPoint*)outers.Next())) {
+ if (outer->fUsed > 0) continue; // points can't be recycled
+ if (inner->DTheta(outer) > fThetaCut) continue;
+ unit = new AliITSneuron;
+ unit->fInner = inner;
+ unit->fOuter = outer;
+ CalcParams(unit);
+ abscurv = TMath::Abs(unit->fCurv);
+ if (unit->fDiff > fDiff || abscurv < min || abscurv > max) {
+ delete unit;
+ continue;
+ }
+ unit->fActivation = gRandom->Rndm() * (fMax - fMin) + fMin;
+ fNeurons[l]->AddLast(unit);
+ neurons++;
+ } // end loop on candidates for outer point for neurons
+ }
+ } // end loop on candidates inner points for neuron
+ }
+ } // end loop on layers
+ return neurons;
+}
+//---------------------------------------------------------------------------------------------------------
+Bool_t AliITSneuralTracker::Update()
+{
+ // Updating cycle
+ // --------------
+ // Performs an updating cycle, by summing all the
+ // gain and cost contribution for each neuron and
+ // updating its activation.
+ // An asynchronous method is followed, with the order
+ // according that the neurons have been created
+ // Time wasting is avoided by dividing the neurons
+ // into different arrays, depending on the layer of their inner point.
+ // The return value answers the question: "The network has stabilized?"
+
+ Int_t j, l;
+ Double_t actOld, actNew, argNew, totDiff = 0.0;
+ Double_t sumInh, sumExc, units = 0.0;
+ AliITSneuron *me = 0, *it = 0;
+ for (l = 0; l < 5; l++) {
+ TObjArrayIter meIter(fNeurons[l]);
+ while ((me = (AliITSneuron*)meIter.Next())) {
+ units++;
+ TObjArrayIter inhIter(fNeurons[l]);
+
+ // inhibition (search only for neurons starting in the same layer)
+ sumInh = 0.0;
+ while((it = (AliITSneuron*)inhIter.Next())) {
+ if (it->fOuter == me->fOuter || it->fInner == me->fInner)
+ sumInh += it->fActivation;
+ }
+
+ // excitation (search only for neurons
+ // which start in the previous or next layers)
+ sumExc = 0.0;
+ for (j = l - 1; j <= l + 1; j += 2) {
+ if (j < 0 || j >= 5) continue;
+ TObjArrayIter itIter(fNeurons[j]);
+ while ((it = (AliITSneuron*)itIter.Next())) {
+ if (it->fInner == me->fOuter || it->fOuter == me->fInner && it->fCurv * me->fCurv > 0.0) {
+ sumExc += Weight(me, it) * it->fActivation;
+ }
+ }
+ } // end search for excitories
+
+ actOld = me->fActivation;
+ argNew = fRatio * sumExc - sumInh;
+ actNew = 1.0 / (1.0 + TMath::Exp(-argNew / fTemp));
+ me->fActivation = actNew;
+ // calculation the relative activation variation
+ // (for stabilization check)
+ totDiff += TMath::Abs((actNew - actOld) / actOld);
+ } // end loop over 'me' (updated neuron)
+ } // end loop over layers
+ totDiff /= units;
+
+ return (totDiff < fVar);
+}
+//---------------------------------------------------------------------------------------------------------
+Int_t AliITSneuralTracker::Save()
+{
+ // Tracki saving method
+ // --------------------
+ // Stores the recognized tracks into the TObjArray 'fTracks'
+ // but doesn't save it into a file until the whole work has ended
+
+ Int_t l;
+ Double_t test = 0.5;
+
+ // Before all, for each neuron is found the best active sequences
+ // among the incoming and outgoing other units
+ AliITSneuron *main, *fwd, *back, *start;
+ for (l = 0; l < 5; l++) {
+ TObjArrayIter mainIter(fNeurons[l]);
+ while((main = (AliITSneuron*)mainIter.Next())) {
+ if (l < 4) {
+ TObjArrayIter fwdIter(fNeurons[l+1]);
+ test = 0.5;
+ while ((fwd = (AliITSneuron*)fwdIter.Next())) {
+ if (main->fOuter == fwd->fInner && fwd->fActivation > test) {
+ test = fwd->fActivation;
+ main->fNext = fwd;
+ }
+ };
+ }
+ if (l > 0) {
+ TObjArrayIter backIter(fNeurons[l-1]);
+ test = 0.5;
+ while ((back = (AliITSneuron*)backIter.Next())) {
+ if (main->fInner == back->fOuter && back->fActivation > test) {
+ test = back->fActivation;
+ main->fPrev = back;
+ }
+ };
+ }
+ }
+ }
+
+ // Then, are resolved the mismatches in the chains found above
+ // (when unit->next->prev != unit or unit->prev->next != unit)
+ for (l = 0; l < 5; l++) {
+ TObjArrayIter mainIter(fNeurons[l]);
+ while((main = (AliITSneuron*)mainIter.Next())) {
+ fwd = main->fNext;
+ back = main->fPrev;
+ if (fwd != 0 && fwd->fPrev != main) main->fNext = 0;
+ if (back != 0 && back->fNext != main) main->fPrev = 0;
+ }
+ }
+
+ Int_t pos, neg, incr, decr, ntracks;
+ AliITSneuralTrack *trk = 0;
+
+ TObjArrayIter startIter(fNeurons[0]);
+ for (;;) {
+ start = (AliITSneuron*)startIter.Next();
+ if (!start) break;
+ Int_t pts = 0;
+ // with the chain above defined, a track can be followed
+ // like a linked list structure
+ for (main = start; main; main = main->fNext) pts++;
+ // the count will be > 5 only if the track contains 6 points
+ // (what we want)
+ if (pts < 5) continue;
+ // track storage
+ trk = new AliITSneuralTrack;
+ trk->CopyPoint(start->fInner);
+ for (main = start; main; main = main->fNext) {
+ //main->fInner->fUsed = kTRUE;
+ //main->fOuter->fUsed = kTRUE;
+ trk->CopyPoint(main->fOuter);
+ }
+ trk->Kinks(pos, neg, incr, decr);
+ if (pos != 0 && neg != 0) {
+ cout << "pos, neg kinks = " << pos << " " << neg << endl;
+ continue;
+ }
+ if (incr != 0 && decr != 0) {
+ cout << "increment, decrements in delta phi = " << incr << " " << decr << endl;
+ continue;
+ }
+ for (main = start; main; main = main->fNext) {
+ main->fInner->fUsed++;
+ main->fOuter->fUsed++;
+ }
+ fTracks->AddLast(trk);
+ }
+ ntracks = (Int_t)fTracks->GetEntriesFast();
+ return ntracks;
+}
+//---------------------------------------------------------------------------------------------------------
+void AliITSneuralTracker::DoSector(Int_t sect, Bool_t flag)
+{
+ // Sector recognition
+ // ------------------
+ // This is a private method which works on a single sector
+ // just for an easier readability of the code.
+ // The 'flag' argument is needed to be given to the 'Initialize' method (see it)
+
+ Int_t curvnum, nTracks = 0, nUnits;
+ Bool_t isStable = kFALSE;
+ for (curvnum = 0; curvnum < fCurvNum; curvnum++) {
+ cout << "\rSector " << sect << ":" << ((curvnum+1)*100)/fCurvNum << "%" << flush;
+ nUnits = Initialize(sect, curvnum, flag);
+ if (!nUnits) continue;
+ do {
+ isStable = Update();
+ } while (!isStable);
+ nTracks = Save();
+ }
+ cout << "\rTracks stored after sector " << sect << ": " << nTracks << endl;
+}
+//---------------------------------------------------------------------------------------------------------
+Int_t AliITSneuralTracker::CountExits(AliITSneuron *n)
+{
+ // Counter for neurons which enter in 'n' tail
+
+ Int_t count = 0, l = n->fOuter->fLayer;
+ if (l == 5) return 0;
+
+ AliITSneuron *test = 0;
+ TObjArrayIter iter(fNeurons[l]);
+ while ((test = (AliITSneuron*)iter.Next())) {
+ if (test->fInner == n->fOuter) count++;
+ }
+ return count;
+}
+//---------------------------------------------------------------------------------------------------------
+Int_t AliITSneuralTracker::CountEnters(AliITSneuron *n)
+{
+ // Counter for neurons which exit from 'n' head
+
+ Int_t count = 0, l = n->fInner->fLayer;
+ if (l == 0) return 0;
+
+ AliITSneuron *test = 0;
+ TObjArrayIter iter(fNeurons[l-1]);
+ while ((test = (AliITSneuron*)iter.Next())) {
+ if (test->fOuter == n->fInner) count++;
+ }
+ return count;
+}
+//---------------------------------------------------------------------------------------------------------
+Bool_t AliITSneuralTracker::ResetNeuron(AliITSneuron *n)
+{
+ // A method which resets the neuron's pointers to zero
+
+ if (!n) return kFALSE;
+ n->fNext = 0;
+ n->fPrev = 0;
+ n->fInner = 0;
+ n->fOuter = 0;
+ // the other datamembers don't need to be reinitialized
+ return kTRUE;
+}
+//---------------------------------------------------------------------------------------------------------
+Bool_t AliITSneuralTracker::SetEdge(AliITSneuron *n, AliITSglobalRecPoint *p, const char what)
+{
+ // Sets the indicated edge point to the neuron
+ // (if the arguments are suitable ones...)
+
+ if (!n || !p || (what != 'i' && what != 'I' && what != 'o' && what != 'O')) return kFALSE;
+ if (what == 'i' || what == 'I')
+ n->fInner = p;
+ else
+ n->fOuter = p;
+ return kTRUE;
+}
+//---------------------------------------------------------------------------------------------------------
+Bool_t AliITSneuralTracker::CalcParams(AliITSneuron *n)
+{
+ // This method evaluates the helix parameters from the edged of a neuron
+ // (curvature, phase parameter, tangent of lambda angle)
+ // if the p[arameters assume too strange (and unphysical) values, the
+ // method return kFALSE, else it return kTRUE;
+
+ if (!n) return kFALSE;
+
+ Double_t sign(1.0), det, rc = 0.0;
+ Double_t tanl1, tanl2, l1, l2;
+ // changing the reference frame into the one centered in the vertex
+ Double_t x1 = n->fInner->fGX - fVPos[0];
+ Double_t y1 = n->fInner->fGY - fVPos[1];
+ Double_t z1 = n->fInner->fGZ - fVPos[2];
+ Double_t r1 = x1 * x1 + y1 * y1;
+ Double_t x2 = n->fOuter->fGX - fVPos[0];
+ Double_t y2 = n->fOuter->fGY - fVPos[1];
+ Double_t z2 = n->fOuter->fGZ - fVPos[2];
+ Double_t r2 = x2 * x2 + y2 * y2;
+ // initialization of the XY-plane data (curvature and centre)
+ // (will remain these values if we encounter errors)
+ n->fCurv = n->fCX = n->fCY = n->fTanL = n->fPhase = 0.0;
+ n->fLength = (n->fOuter->fGX - n->fInner->fGX) * (n->fOuter->fGX - n->fInner->fGX);
+ n->fLength += (n->fOuter->fGY - n->fInner->fGY) * (n->fOuter->fGY - n->fInner->fGY);
+ n->fLength += (n->fOuter->fGZ - n->fInner->fGZ) * (n->fOuter->fGZ - n->fInner->fGZ);
+ n->fLength = TMath::Sqrt(n->fLength);
+ // this determinant which is zero when the points are in a line
+ det = 2.0 * (x1 * y2 - x2 * y1);
+ if (det == 0.0)
+ return kFALSE; // it is difficult having perfectly aligned points --- it's more likely to be an error
+ else {
+ n->fCX = (r1*y2 - r2*y1) / det;
+ n->fCY = (r2*x1 - r1*x2) / det;
+ rc = TMath::Sqrt(n->fCX * n->fCX + n->fCY * n->fCY);
+ r1 = TMath::Sqrt(r1);
+ r2 = TMath::Sqrt(r2);
+ sign = (n->fOuter->fPhi >= n->fInner->fPhi) ? 1.0 : -1.0;
+ if (rc > 0.0) {
+ n->fCurv = sign / rc;
+ n->fPhase = TMath::ATan2(-n->fCX, -n->fCY);
+ if (n->fPhase < 0.0) n->fPhase += 2.0 * TMath::Pi(); // angles in 0 --> 2pi
+ if (r1 <= 2.0 * rc && r2 <= 2.0 * rc) {
+ l1 = 2.0 * rc * TMath::ASin(r1 / (2.0 * rc));
+ l2 = 2.0 * rc * TMath::ASin(r2 / (2.0 * rc));
+ tanl1 = z1 / l1;
+ tanl2 = z2 / l2;
+ n->fDiff = TMath::Abs(tanl1 - tanl2);
+ n->fTanL = 0.5 * (tanl1 + tanl2);
+ }
+ else {
+ n->fTanL = 0.0;
+ n->fDiff = 100.0;
+ return kFALSE; // it' even more difficult that the arguments above aren't acceptable for arcsine
+ }
+ }
+ else
+ return kFALSE;
+ }
+ return kTRUE;
+}
+//---------------------------------------------------------------------------------------------------------
+Double_t AliITSneuralTracker::Angle(AliITSneuron *n, AliITSneuron *m)
+{
+ // calculates the angle between two segments
+ // but doesn't check if they have a common point.
+ // The calculation is made by the ratio of their scalar product
+ // to the product of their magnitudes
+
+ Double_t x = (m->fOuter->fGX - m->fInner->fGX) * (n->fOuter->fGX - n->fInner->fGX);
+ Double_t y = (m->fOuter->fGY - m->fInner->fGY) * (n->fOuter->fGY - n->fInner->fGY);
+ Double_t z = (m->fOuter->fGZ - m->fInner->fGZ) * (n->fOuter->fGZ - n->fInner->fGZ);
+
+ Double_t cosine = (x + y + z) / (n->fLength * m->fLength);
+
+ if (cosine > 1.0 || cosine < -1.0) {
+ Warning("AliITSneuronV1::Angle", "Strange value of cosine");
+ return 0.0;
+ }
+
+ return TMath::ACos(cosine);
+}
+//---------------------------------------------------------------------------------------------------------
+Double_t AliITSneuralTracker::Weight(AliITSneuron *n, AliITSneuron *m)
+{
+ // calculation of the excitoy weight
+
+// Double_t v1 = TMath::Abs(fTanL - n->fTanL);
+// Double_t v2 = TMath::Abs(TMath::Abs(fCurv) - TMath::Abs(n->fCurv));
+// Double_t v3 = TMath::Abs(fPhase - n->fPhase);
+// Double_t q = (v1 + v2 + v3) / 3.0;
+// Double_t a = 1.0; // 0.02 / (fDiff + n->fDiff);
+ Double_t b = 1.0 - TMath::Sin(Angle(n, m));
+ return TMath::Power(b, fExp);
+// return TMath::Exp(q / par1);
+
+}
+//---------------------------------------------------------------------------------------------------------
+// ====================================================================================================
+// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> AliITSneuron <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+// ====================================================================================================
+
+AliITSneuron::AliITSneuron()
+{
+ // Default constructor
+ // -------------------
+ // sets all parameters to default (but unuseful)
+ // (in particular all pointers are set to zero)
+ // There is no alternative constructor, because
+ // this class should not be used by itself, but
+ // it is only necessary to allow the neural tracker to work
+ // In fact, all methods which operate on neuron's datamembers
+ // are stored in the AliITSneuralTracker class anyway
+
+ fActivation = 0.0;
+
+ fCurv = 0.0;
+ fCX = 0.0;
+ fCY = 0.0;
+ fTanL = 0.0;
+ fPhase = 0.0;
+ fDiff = 0.0;
+ fLength = 0.0;
+
+ fNext = fPrev = 0;
+ fInner = fOuter = 0;
+}
+//---------------------------------------------------------------------------------------------------------
+AliITSneuron::~AliITSneuron()
+{
+ // Destructor
+ // ----------
+ // does nothing, because there is no need to affect parameters,
+ // while the AliITSglobalRecPoint pointers can't be deleted,
+ // and the AliITSneuron pointers belong to a TObjArray which will
+ // be cleared when necessary...
+}
+
--- /dev/null
+#ifndef AliITSneuralTracker_h
+#define AliITSneuralTracker_h
+
+class TObjArray;
+class AliITSneuron;
+class AliITSneuralTrack;
+class AliITSglobalRecPoint;
+
+///////////////////////////////////////////////////////////////////////
+//
+// 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:
+
+ // constructors & destructor
+ AliITSneuralTracker();
+ AliITSneuralTracker(Int_t nsecs, Int_t nc, Double_t *c, Double_t theta);
+ virtual ~AliITSneuralTracker();
+
+ // file reading and points array population
+ Int_t ReadFile(const Text_t *fname, Int_t evnum);
+
+ // setters for working parameters (NOT cuts)
+ void SetDiff(Double_t a) {fDiff=a;} // helix adaptment parameter
+ void SetExponent(Double_t a) {fExp = a;} // exponent which selects the high excitory values
+ void SetGainToCostRatio(Double_t a) {fRatio = a;} // ratio between the gain and cost contributions
+ void SetTemperature(Double_t a) {fTemp = a;} // temperature parameter in activation function
+ void SetVariationLimit(Double_t a) {fVar = a;} // stability threshold
+ void SetInitLimits(Double_t a, Double_t b) // intervals for initial random activations
+ {fMin = (a<=b) ? a : b; fMax = (a<=b)? b : a;}
+ void SetVertex(Double_t x, Double_t y, Double_t z) // estimated vertex position
+ { fVPos[0] = x; fVPos[1] = y; fVPos[2] = z;}
+
+ // neuron managin methods
+ Bool_t SetEdge(AliITSneuron *n, AliITSglobalRecPoint *p, const char what); // sets the neuron's edges
+ Bool_t CalcParams(AliITSneuron *n); // calculation of helix parameters for the neuron
+
+ // neural tracking (the argument is the ROOT file to store results)
+ void Go(const char* filesave, Bool_t flag = kFALSE);
+
+private:
+
+ // These methods are private to prevent the users to use them in incorrect sequences
+
+ // neuron management methods (not to be used publically)
+ Bool_t ResetNeuron(AliITSneuron *n); // resets all neuron's fields like in the constructor
+ Int_t CountExits (AliITSneuron *n); // counter for # of neurons going into n's tail
+ Int_t CountEnters(AliITSneuron *n); // counter for # of neurons startins from n's head
+ Double_t Angle(AliITSneuron *n, AliITSneuron *m); // angle between two neural segments
+ Double_t Weight(AliITSneuron *n, AliITSneuron *m); // synaptic weight between two neural segments
+
+ // neural network work-flow
+ Int_t Initialize(Int_t snum, Int_t cnum, Bool_t flag = kFALSE);
+ // resets (or creates) neurons array for a new neural tracking
+ Bool_t Update(); // updates the neurons and checks if stabilization has occurred
+ Int_t Save(); // saves the found tracks after stabilization
+ void DoSector(Int_t sect, Bool_t flag); // complete work on a single sector
+
+ Int_t fCurvNum; // # of curvature cut steps
+ Double_t *fCurvCut; //! value of all curvature cuts
+
+ Int_t fSectorNum; // no. of azymuthal sectors
+ Double_t fSectorWidth; // width of an azymuthal sector
+
+ Double_t fVPos[3]; // estimated vertex coords
+
+ Double_t fMin; // min initial random activations
+ Double_t fMax; // max initial random activations
+ Double_t fThetaCut; // polar angle cut
+ Int_t fThetaNum; // size of theta sectionement
+ Double_t fTemp; // logistic function parameter
+ Double_t fVar; // stability threshold (for mean rel. activations)
+ Double_t fRatio; // ratio between inhibitory and excitory contributions
+ Double_t fExp; // alignment weight
+ Double_t fDiff; // max allowed difference between TanL exstimations from the two neuron edges
+
+ TObjArray ***fPoints[6]; //! Collection of recpoints (sectioned in azym. secs)
+ TObjArray *fNeurons[6]; //! Collection of neurons
+ TObjArray *fTracks; //! Collection of tracks
+
+ ClassDef(AliITSneuralTracker, 1)
+};
+
+
+////////////////////////////////////////////////////////////////////////////////
+
+
+class AliITSneuron : public TObject {
+
+ friend class AliITSneuralTracker;
+
+public:
+
+ AliITSneuron();
+ virtual ~AliITSneuron();
+
+ Bool_t ContainsID(Int_t ID)
+ { return (fInner->HasID(ID) && fOuter->HasID(ID)); }
+
+private:
+
+ Double_t fActivation; // Activation value
+
+ Double_t fCurv; // curvature [= 1 / R]
+ Double_t fCX; // curvature centre (X) in changed RF
+ Double_t fCY; // curvature centre (X) in changed RF
+ Double_t fTanL; // tan(dip angle) = C / freq
+ Double_t fPhase; // 'phase' parameter
+ Double_t fDiff; // difference between tan(lambda) estimated by the two different points
+ Double_t fLength; // segment length
+
+ AliITSneuron *fNext; // best outgoing unit
+ AliITSneuron *fPrev; // best incoming unit
+
+ AliITSglobalRecPoint *fInner; // inner point
+ AliITSglobalRecPoint *fOuter; // outer point
+};
+
+
+#endif
#pragma link C++ class AliCascadeVertexer+;
#pragma link C++ class AliITSVertex+;
+
+// Classes for neural tracking
+#pragma link C++ class AliITSglobalRecPoint+;
+#pragma link C++ class AliITSneuron+;
+#pragma link C++ class AliITSneuralTrack+;
+#pragma link C++ class AliITSneuralTracker+;
+
#endif
--- /dev/null
+// Neural performance evaluation
+//
+// before using this macro, you should have make run the macro
+// 'ITSstoreFindableTracks.C' in order to have a file named following
+// the rule to append '_fnd' to the name
+// (e.g. galice.root --> galice_fnd.root)
+// and have in it a tree with some tracks statistics useful for evaluation
+// of performance with the same event without loosing too much time...
+
+void ITSneuralCompleteEval
+(Int_t nsecs, const char *filename, Bool_t low = kFALSE,
+ Bool_t draw = kTRUE, const char *save = 0)
+{
+ Int_t N, M;
+ Float_t *edge = 0;
+ if (!low) {
+ N = 7;
+ edge = new Float_t[8];
+ edge[0] = 0.0;
+ edge[1] = 0.5;
+ edge[2] = 1.0;
+ edge[3] = 1.5;
+ edge[4] = 2.0;
+ edge[5] = 3.0;
+ edge[6] = 4.0;
+ edge[7] = 5.0;
+ }
+ else {
+ N = 5;
+ edge = new Float_t[6];
+ edge[0] = 0.0;
+ edge[1] = 0.2;
+ edge[2] = 0.4;
+ edge[3] = 0.6;
+ edge[4] = 0.8;
+ edge[5] = 1.0;
+ }
+
+ Float_t find[2] = {0.0, 0.0}, find1[2] = {0.0, 0.0};
+ Float_t good[2] = {0.0, 0.0}, good1[2] = {0.0, 0.0};
+ Float_t fake[2] = {0.0, 0.0}, fake1[2] = {0.0, 0.0};
+
+ // histos filled with neural results
+ TH1D *hgood[2], *hfake[2], *hfind[2];
+ TH1D *hgood1[2], *hfake1[2], *hfind1[2];
+ TH1D *hg[2], *hf[2];
+ if (draw) {
+ hgood[0] = new TH1D("hgood0", "Good found tracks", N, edge);
+ hfake[0] = new TH1D("hfake0", "Fake found tracks", N, edge);
+ hfind[0] = new TH1D("hfound0", "Findable tracks", N, edge);
+ hgood[1] = new TH1D("hgood1", "Good found tracks", N, edge);
+ hfake[1] = new TH1D("hfake1", "Fake found tracks", N, edge);
+ hfind[1] = new TH1D("hfound1", "Findable tracks", N, edge);
+
+ hgood[0]->Sumw2();
+ hfake[0]->Sumw2();
+ hfind[0]->Sumw2();
+ hgood[1]->Sumw2();
+ hfake[1]->Sumw2();
+ hfind[1]->Sumw2();
+
+ // histos for evaluating percentual efficiency
+ hg[0] = new TH1D("hg0", "Efficiency (%) for #geq 5 right pts.", N, edge);
+ hf[0] = new TH1D("hf0", "Fake probability (%) for #geq 5 right pts.", N, edge);
+ hg[1] = new TH1D("hg1", "Efficiency (%) for 6 right pts.", N, edge);
+ hf[1] = new TH1D("hf1", "Fake probability (%) for 6 right pts.", N, edge);
+ }
+
+ TFile *ffind;
+ TTree *tree;
+
+ if (low) {
+ ffind = new TFile(Form("%s.root", filename), "READ");
+ tree = (TTree*)ffind->Get("TreeF");
+ }
+ else {
+ ffind = new TFile(Form("%s_fnd.root", filename), "READ");
+ tree = (TTree*)ffind->Get("tree");
+ }
+
+ TFile *ftracks = new TFile(Form("%s_%d.root", filename, nsecs), "READ");
+
+
+ Double_t pt;
+ Int_t i, j, count, prim, *none = 0, div;
+ Int_t entries = tree->GetEntries(), label, min[] = {5,6};
+ tree->SetBranchAddress("pt", &pt);
+ tree->SetBranchAddress("count", &count);
+ tree->SetBranchAddress("prim", &prim);
+
+ AliITSneuralTrack *trk = 0;
+ div = low ? 50 : 500;
+ for(i = 1;;i++) {
+ trk = (AliITSneuralTrack*)ftracks->Get(Form("AliITSneuralTrack;%d", i));
+ if (i%div == 0) cout << "\rEvaluating found track " << i << flush;
+ if (!trk) break;
+ for (j = 0; j < 2; j++) {
+ label = trk->EvaluateTrack(0, min[j], none);
+ tree->GetEntry(abs(label));
+ if (count >= min[j]) {
+ if (label > 0) {
+ if (draw) hgood[j]->Fill(pt);
+ good[j]++;
+ if (pt >= 1.0) good1[j]++;
+ }
+ else {
+ if (draw) hfake[j]->Fill(pt);
+ fake[j]++;
+ if (pt >= 1.0) fake1[j]++;
+ }
+ }
+ }
+ }
+ cout << endl;
+
+ div = low ? 200 : 20000;
+
+ for (i = 0; i < entries; i++) {
+ if (i%div == 0) cout << "\rEvaluating findable track no. " << i << flush;
+ tree->GetEntry(i);
+ for (j = 0; j < 2; j++) {
+ if (count >= min[j]) {
+ find[j]++;
+ if (draw) hfind[j]->Fill(pt);
+ if (pt >= 1.0) find1[j]++;
+ }
+ }
+ }
+ cout << endl;
+ cout << hgood[0]->GetEntries() << " " << hgood[1]->GetEntries() << endl;
+ cout << hfake[0]->GetEntries() << " " << hfake[1]->GetEntries() << endl;
+ cout << hfind[0]->GetEntries() << " " << hfind[1]->GetEntries() << endl << endl;
+
+ if (draw) {
+ TCanvas *canv[2];
+ canv[0] = new TCanvas("c_0", "Tracking efficiency (soft)", 0, 0, 600, 500);
+ canv[1] = new TCanvas("c_1", "Tracking efficiency (hard)", 630, 0, 600, 500);
+
+ TLine *line1 = new TLine(1,100.0,edge[N],100.0); line1->SetLineStyle(4);
+ TLine *line2 = new TLine(1,90,edge[N],90); line2->SetLineStyle(4);
+
+ Bool_t good_drawn;
+ for (i = 0; i < 2; i++) {
+ canv[i]->cd();
+ good_drawn = kFALSE;
+ if (hgood[i]->GetEntries() > 0.0) {
+ good_drawn = kTRUE;
+ hg[i]->Divide(hgood[i], hfind[i], 100.0, 1.0);
+ hg[i]->SetMaximum(120);
+ hg[i]->SetMinimum(0);
+ hg[i]->SetMarkerStyle(21);
+ hg[i]->SetMarkerSize(1);
+ hg[i]->SetStats(kFALSE);
+ hg[i]->GetXaxis()->SetTitle("pt (GeV/c)");
+ hg[i]->Draw("PE1");
+ }
+ if (hfake[i]->GetEntries() > 0.0) {
+ hf[i]->Divide(hfake[i], hfind[i], 100.0, 1.0);
+ hf[i]->SetMaximum(120);
+ hf[i]->SetMinimum(0);
+ hf[i]->SetMarkerStyle(25);
+ hf[i]->SetMarkerSize(1);
+ hf[i]->SetStats(kFALSE);
+ if (good_drawn)
+ hf[i]->Draw("PE1SAME");
+ else
+ hf[i]->Draw("PE1");
+ }
+ line1->Draw("histosame");
+ line2->Draw("histosame");
+ canv[i]->Update();
+ }
+ canv[0]->SaveAs(Form("%s_soft.eps", filename));
+ canv[1]->SaveAs(Form("%s_hard.eps", filename));
+ cout << endl;
+ }
+
+ Float_t sgood[2] = {0.0, 0.0}, sgood1[2] = {0.0, 0.0};
+ Float_t sfake[2] = {0.0, 0.0}, sfake1[2] = {0.0, 0.0};
+ for (i = 0; i < 2; i++) {
+ sgood[i] = error(good[i], find[i]);
+ sgood1[i] = error(good1[i], find1[i]);
+ sfake[i] = error(fake[i], find[i]);
+ sfake1[i] = error(fake1[i], find1[i]);
+
+ good[i] = good[i] * 100.0 / find[i];
+ fake[i] = fake[i] * 100.0 / find[i];
+ good1[i] = good1[i] * 100.0 / find1[i];
+ fake1[i] = fake1[i] * 100.0 / find1[i];
+ }
+
+ if (save) {
+ fstream data(save, ios::app);
+ data.setf(ios::fixed);
+ data.precision(1);
+ data << good1[0] << " " << fake1[0] << " ";
+ data << good1[1] << " " << fake1[1] << endl;
+ data.close();
+ }
+ else {
+ cout.setf(ios::fixed);
+ cout.precision(1);
+ cout << "*****************************************" << endl;
+ cout << "* Tracks with at least 5 correct points *" << endl;
+ cout << "*****************************************" << endl;
+ cout << "(all particles)" << endl;
+ cout << "Efficiency: " << good[0] << " +/- " << sgood[0] << "%" << endl;
+ cout << "Fake prob.: " << fake[0] << " +/- " << sfake[0] << "%" << endl;
+ if (!low) {
+ cout << "(pt >= 1 GeV/c)" << endl;
+ cout << "Efficiency: " << good1[0] << " +/- " << sgood1[0] << "%" << endl;
+ cout << "Fake prob.: " << fake1[0] << " +/- " << sfake1[0] << "%" << endl;
+ }
+ cout << endl;
+ cout << "************************************" << endl;
+ cout << "* Tracks with all 6 correct points *" << endl;
+ cout << "************************************" << endl;
+ cout << "(all particles)" << endl;
+ cout << "Efficiency: " << good[1] << " +/- " << sgood[1] << "%" << endl;
+ cout << "Fake prob.: " << fake[1] << " +/- " << sfake[1] << "%" << endl;
+ if (!low) {
+ cout << "(pt >= 1 GeV/c)" << endl;
+ cout << "Efficiency: " << good1[1] << " +/- " << sgood1[1] << "%" << endl;
+ cout << "Fake prob.: " << fake1[1] << " +/- " << sfake1[1] << "%" << endl;
+ }
+ cout << endl;
+ }
+}
+
+Double_t error(Double_t x, Double_t y) {
+ Double_t radq = x + x * x / y;
+ return 100.0 * sqrt(radq) / y;
+}
--- /dev/null
+// ARGUMENTS:
+// 1. number of azymuthal sectors (it's better not to go under 8 or over 40)
+// 2. the ROOT file to read (WITHOUT exstension)
+// 3. event number
+// 4. if specified a string, a fstream named like the argument is opened and
+// the elapsed CPU time is stored (not useful)
+
+// the macro will save a file named, for example "galice_<nsecs>.root"
+// containing may AliITSneuralTrack objects
+
+void ITSneuralTracking
+(Int_t nsecs = 12, const char* rfile = "galice", Int_t event = 0, const char* save = 0)
+{
+ TStopwatch timer;
+ Double_t CONVERT = TMath::Pi() / 180.0;
+ const char* wfile = Form("%s_%d.root", rfile, nsecs);
+ cout << "Reading file " << rfile << ".root and saving in " << wfile << endl;
+
+// ==================================
+// ==== CURVATURE CUT DEFINITION ====
+// ==================================
+
+ // These values define the curvature cuts for all steps
+ // within a sector.
+ // For a greater clarity, the cuts are indicated in units
+ // of transverse momentum (GeV/c) but these value have no
+ // exact physical meaning, but are useful to understand
+ // well what means a choice in the value of a certain
+ // curvature constraint
+ // NOTE: becareful to make sure that the 'ncuts' variable
+ // have the same value of the dimension of the allocated arrays
+
+ Int_t ncuts;
+ Double_t *p, *cut;
+
+ ncuts = 5;
+ p = new Double_t[5];
+ cut = new Double_t[5];
+ p[0] = 2.0;
+ p[1] = 1.0;
+ p[2] = 0.7;
+ p[3] = 0.5;
+ p[4] = 0.3;
+
+ for (Int_t i = 0; i < ncuts; i++) cut[i] = 0.003 * 0.2 / p[i];
+
+
+// ==========================
+// ==== OTHER PARAMETERS ====
+// ==========================
+
+ Bool_t flag = kFALSE; // for now, don't change this line, please...
+
+ Double_t diff = 0.02; // helicoidal cut
+ Double_t dtheta = 1.0; // delta-theta cut
+ Double_t temp = 1.0; // temperature parameter
+ Double_t var = 0.0001; // stabilization threshold
+
+ Double_t exp = 7.0; // straight-line excitator
+ Double_t gtoc = 3.0; // gain/cost contribution ratio
+
+ Double_t min = 0.4; // minimum in random activation initialization
+ Double_t max = 0.6; // maximum in random activation initialization
+
+
+// =========================
+// ==== NEURAL TRACKING ====
+// =========================
+
+ AliITSneuralTracker *ANN = new AliITSneuralTracker(nsecs, ncuts, cut, CONVERT*dtheta);
+
+ TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(Form("%s.root", rfile));
+ if (!file) file = new TFile(Form("%s.root", rfile),"UPDATE");
+
+ //Double_t Xv = -0.001097;
+ //Double_t Yv = -0.00347647;
+ //Double_t Zv = 0.000631345;
+ //ANN->SetVertex(Xv, Yv, Zv);
+ // You should find the vertex with VertexMacro.C
+ // and then put by hand the found values with
+ // the above method.
+
+ Int_t points = ANN->ReadFile(Form("%s.root", rfile), event);
+
+ ANN->SetTemperature(temp);
+ ANN->SetVariationLimit(var);
+ ANN->SetGainToCostRatio(gtoc);
+ ANN->SetExponent(exp);
+ ANN->SetInitLimits(min, max);
+ ANN->SetDiff(diff);
+
+ cout << points << " points found " << endl << endl;
+
+ TStopwatch timer;
+ timer.Start();
+
+ ANN->Go(wfile, flag);
+
+ timer.Stop();
+ cout << endl;
+ timer.Print();
+
+ if (save) {
+ fstream ksave(save, ios::app);
+ ksave << nsecs << " " << timer->CpuTime() << endl;
+ ksave.close();
+ }
+
+// delete gAlice;
+// gAlice = 0;
+}
--- /dev/null
+void ITSstoreFindableTracks(const char *nfile = "galice", Int_t evnum = 0)
+{
+ TFile *froot = new TFile(Form("%s.root", nfile), "READ");
+
+ gAlice = (AliRun*)froot->Get("gAlice");
+ if (!gAlice) {
+ cout << "gAlice not found in file!!!" << endl;
+ return;
+ }
+
+ AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
+
+ Int_t nparts = gAlice->GetEvent(evnum);
+ cout << "Tracks saved in event " << evnum << ": " << nparts << endl;
+
+ TClonesArray *recPoints = ITS->RecPoints();
+ TTree *treeRec = gAlice->TreeR();
+ Int_t ntracks = gAlice->GetNtrack(); //FCA correction
+ Int_t nmodules = treeRec->GetEntries();
+ Int_t modmin[6];
+ Int_t modmax[6];
+
+ Int_t layer, nlads, ndets;
+ AliITSgeom *gm = ITS->GetITSgeom();
+ for (layer = 0; layer < 6; layer++) {
+ nlads = gm->GetNladders(layer+1);
+ ndets = gm->GetNdetectors(layer+1);
+ modmin[layer] = gm->GetModuleIndex(layer+1, 1, 1);
+ modmax[layer] = gm->GetModuleIndex(layer+1, nlads, ndets);
+ }
+
+ Int_t track, *goodITS[6];
+ for (layer = 0; layer < 6; layer++) {
+ goodITS[layer] = new Int_t[ntracks];
+ for(track = 0; track < ntracks; track++) goodITS[layer][track] = 0;
+ }
+
+ Int_t irec, index, npoints = 0;
+ for (Int_t layer = 1; layer <= 6; layer++) {
+ for (Int_t mod = modmin[layer-1]; mod <= modmax[layer-1]; mod++) {
+ ITS->ResetRecPoints();
+ treeRec->GetEvent(mod);
+ cout << "\rLayer " << layer << ", module: " << mod << flush;
+ TObjArrayIter *iter = (TObjArrayIter*)recPoints->MakeIterator();
+ while ((recp = (AliITSRecPoint*)iter.Next())) {
+ for (index = 0; index < 3; index++) {
+ track = recp->GetLabel(index);
+ if(track < 0) continue;
+ if(track > ntracks) {
+ cerr << " Error on track number \n";
+ continue;
+ }
+ goodITS[layer-1][track] = 1;
+ } //loop over tracks
+ } //loop over points
+ } //loop over modules
+ cout << endl;
+ } //loop over layers
+
+ // histos filled with neural results
+ TFile *fnew = new TFile(Form("%s_fnd.root", nfile), "RECREATE");
+ TParticle *part = 0;
+ TTree *tree = new TTree("tree", "Tree of tracks");
+
+ Int_t count = 0, prim = 0;
+ Double_t pt = 0.0;
+ tree->Branch("count", &count, "count/I");
+ tree->Branch("prim", &prim, "prim/I");
+ tree->Branch("pt", &pt, "pt/D");
+
+ for(track = 0; track < ntracks; track++) {
+ prim = 0;
+ count = 0;
+ part = gAlice->Particle(track);
+ if (!(track % 10000)) cout << "\rStoring track " << track << flush;
+ if (part->GetFirstMother() == -1) prim = 1;
+ for (layer = 0; layer < 6; layer++) count += goodITS[layer][track];
+ pt = part->Pt();
+ tree->Fill();
+ }
+
+ fnew->cd();
+ tree->Write("tree", TObject::kOverwrite);
+ fnew->Close();
+}
AliITSPid.cxx AliITStrackV2Pid.cxx \
AliITSVertex.cxx \
AliV0vertex.cxx AliV0vertexer.cxx \
- AliCascadeVertex.cxx AliCascadeVertexer.cxx
+ AliCascadeVertex.cxx AliCascadeVertexer.cxx \
+ AliITSglobalRecPoint.cxx AliITSneuralTrack.cxx \
+ AliITSneuralTracker.cxx
# AliITSAlignmentTrack.cxx AliITSAlignmentModule.cxx \
# Fortran sources
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff