--- /dev/null
+#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;
+}
--- /dev/null
+#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
--- /dev/null
+#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;
+}
--- /dev/null
+#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
--- /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. *
+ * *
+ * 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;
+}
--- /dev/null
+#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
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff