update of the tracker
authormasera <massimo.masera@cern.ch>
Thu, 8 May 2014 08:14:49 +0000 (10:14 +0200)
committermasera <massimo.masera@cern.ch>
Thu, 8 May 2014 08:14:49 +0000 (10:14 +0200)
ITS/UPGRADE/AliITSUTrackerSA.cxx [changed mode: 0755->0644]
ITS/UPGRADE/AliITSUTrackerSA.h [changed mode: 0755->0644]
ITS/UPGRADE/AliITSUTrackerSAaux.h [changed mode: 0755->0644]

old mode 100755 (executable)
new mode 100644 (file)
index 883f282..adb62ab
@@ -9,30 +9,47 @@ using TMath::Abs;
 using TMath::Sort;
 using TMath::Sqrt;
 #include <TTree.h>
+#include <algorithm>
+using std::sort;
 
 // Vc library
 //#include "Vc/Vc"
 //#include "AliITSUTrackerSAauxVc.h" // Structs and other stuff using Vc library  
-
+#include "AliLog.h"
 #include "AliESDEvent.h"
 #include "AliITSUClusterPix.h"
 #include "AliITSUTrackerSA.h"
+#include "AliITSUReconstructor.h"
+#include "AliITSURecoDet.h"
+#include "AliESDtrack.h"
+
+#include <Riostream.h>
+
+
 
 //#include "AliITSUtrackSA.h"      // Some dedicated SA track class ?  
 
 ClassImp(AliITSUTrackerSA)
 
+const Double_t AliITSUTrackerSA::fgkToler =  1e-6;// tolerance for layer on-surface check
+
 //________________________________________________________________________________
-AliITSUTrackerSA::AliITSUTrackerSA() : AliTracker(), 
-  fClusters(),
+AliITSUTrackerSA::AliITSUTrackerSA(AliITSUReconstructor* rec) :  
+  fReconstructor(rec)
+  ,fITS(0)
+  ,fMatLUT(0)
+  ,fUseMatLUT(kFALSE)
+  ,fCurrMass(0.14)
+  //
+  ,fClusters(),
   fClustersTC(),
   fDoublets(),
   fIndex(),
   fNClusters(),
   fNDoublets(),
   fPhiCut(0.05),
-  fRPhiCut(0.01),
-  fZCut(0.005)
+  fRPhiCut(0.03),
+  fZCut(0.01)
 {
   //--------------------------------------------------------------------
   // This default constructor needs to be provided
@@ -40,11 +57,18 @@ AliITSUTrackerSA::AliITSUTrackerSA() : AliTracker(),
   for(Int_t i=0;i<7;++i) {
     fClusters[i].reserve(5000);
   }
+  if (rec) Init(rec);
 }
 
 //________________________________________________________________________________
 AliITSUTrackerSA::AliITSUTrackerSA(const AliITSUTrackerSA &t): 
   AliTracker(t),
+  fReconstructor(t.fReconstructor),
+  fITS(t.fITS),
+  fMatLUT(t.fMatLUT),
+  fUseMatLUT(t.fUseMatLUT),
+  fCurrMass(t.fCurrMass),
+  //
   fClusters(),
   fClustersTC(),
   fIndex(),
@@ -60,7 +84,36 @@ AliITSUTrackerSA::AliITSUTrackerSA(const AliITSUTrackerSA &t):
 }
 
 //________________________________________________________________________________
-Int_t AliITSUTrackerSA::Clusters2Tracks(AliESDEvent */*event*/) {
+AliITSUTrackerSA::~AliITSUTrackerSA() 
+{
+  // d-tor
+  delete fMatLUT;
+}
+
+
+//_________________________________________________________________________
+void AliITSUTrackerSA::Init(AliITSUReconstructor* rec)
+{
+  // init with external reconstructor
+  //
+  fITS = rec->GetITSInterface();
+  //
+  // create material lookup table
+  const int kNTest = 1000;
+  const double kStepsPerCM=5;
+  fMatLUT  = new AliITSUMatLUT(fITS->GetRMin(),fITS->GetRMax(),Nint(kStepsPerCM*(fITS->GetRMax()-fITS->GetRMin())));
+  double zmn = 1e6;
+  for (int ilr=fITS->GetNLayers();ilr--;) {
+    AliITSURecoLayer* lr = fITS->GetLayer(ilr);
+    if (zmn>Abs(lr->GetZMin())) zmn = Abs(lr->GetZMin());
+    if (zmn>Abs(lr->GetZMax())) zmn = Abs(lr->GetZMax());
+  }
+  fMatLUT->FillData(kNTest,-zmn,zmn);
+  //
+}
+
+//________________________________________________________________________________
+Int_t AliITSUTrackerSA::Clusters2Tracks(AliESDEvent *event) {
   //--------------------------------------------------------------------
   // This is the main tracking function
   // The clusters must already be loaded
@@ -72,16 +125,16 @@ Int_t AliITSUTrackerSA::Clusters2Tracks(AliESDEvent */*event*/) {
 
   // Possibly, (re)fit the found tracks 
 
-  // Tree iterations: 
+  // Three iterations: 
   // - High momentum first;
   // - Low momentum with vertex constraint; 
   // - Everything else. 
 
   MakeDoublets();       // To be implemented
   //MakeTriplets();       // Are triplets really necessary? MFT does not use them. 
-  CASelection();        // To be implemented
-  GlobalFit();          // To be implemented
-  ChiSquareSelection(); // To be implemented
+  CASelection(event);       
+  //GlobalFit();          // To be implemented
+  //ChiSquareSelection(); // To be implemented
 
   return 0;
 }
@@ -111,16 +164,12 @@ Int_t AliITSUTrackerSA::LoadClusters(TTree *cluTree) {
   // This function reads the ITSU clusters from the tree,
   // sort them, distribute over the internal tracker arrays, etc
   //--------------------------------------------------------------------
-  
-  for(Int_t i=0;i<7;++i) {
-    TBranch* br = cluTree->GetBranch(Form("ITSRecPoints%d",i));
-    if (!br) return -1;
-    br->SetAddress(&fClustersTC[i]);
-  }    
-  cluTree->GetEntry(0);
-  
+  fITS->LoadClusters(cluTree);
+  fITS->ProcessClusters();
+  //
   for(int iL=0; iL<7; ++iL) {
-    TClonesArray *clCont=&fClustersTC[iL];
+    fClustersTC[iL]=*fITS->GetLayerActive(iL)->GetClustersAddress();
+    TClonesArray *clCont=fClustersTC[iL];
     fNClusters[iL]=clCont->GetEntriesFast();
     Float_t phi[fNClusters[iL]];
     fIndex[iL] = new Int_t[fNClusters[iL]];
@@ -134,7 +183,9 @@ Int_t AliITSUTrackerSA::LoadClusters(TTree *cluTree) {
     }
     TMath::Sort(fNClusters[iL],phi,fIndex[iL],kFALSE);
   }
-  
+// #ifdef __DEBUG__
+//   PrintInfo("clusters");
+// #endif
   return 0;
 }
 
@@ -160,13 +211,14 @@ AliCluster *AliITSUTrackerSA::GetCluster(Int_t /*index*/) const {
 }
 
 //________________________________________________________________________________
-void AliITSUTrackerSA::CASelection() {
+void AliITSUTrackerSA::CASelection(AliESDEvent *event) {
   // Here it's implemented the Cellular Automaton routine
   // Firstly the level of each doublet is set according to the level of 
   // the neighbour doublets. 
   // Doublet are considered to be neighbour if they share one point and the 
   // phi and theta direction difference of the two is below a cut value.
 
+  cout << "Begin of the CA selection" << endl;
   for( int iL = 1; iL < 6; ++iL ) {
     
     const itsCluster* clusters1 = &fClusters[iL-1][0];
@@ -178,16 +230,19 @@ void AliITSUTrackerSA::CASelection() {
 
     for ( int iD2 = 0; iD2 < fNDoublets[iL]; ++iD2 ) {
       for ( int iD1 = 0; iD1 < fNDoublets[iL-1]; ++iD1 ) {
-       const int id1 = doublets1[iD1].id1;
+       const int id1 = doublets1[iD1].id0;
        const int id2 = doublets2[iD2].id0;
        if ( id1 == id2 ) {
          if ( doublets2[iD2].level <= ( doublets1[iD1].level + 1 ) ) {
            const int id3 = doublets2[iD2].id1;
            const float r3 = Sqrt( clusters3[id3].x * clusters3[id3].x + clusters3[id3].y * clusters3[id3].y );
-           const float r2 = Sqrt( clusters3[id2].x * clusters2[id2].x + clusters2[id2].y * clusters2[id2].y );
-           const float extrZ3 = doublets1[iD1].tanLambda * ( r3 - r2 ) + clusters3[id3].z ;
+           const float r2 = Sqrt( clusters2[id2].x * clusters2[id2].x + clusters2[id2].y * clusters2[id2].y );
+           const float extrZ3 = doublets1[iD1].tanLambda * ( r3 - r2 ) + clusters2[id2].z ;
+           //cout << extrZ3 << " " << clusters3[id3].z << " " << Abs ( extrZ3 - clusters3[id3].z ) << endl;
            if ( Abs ( extrZ3 - clusters3[id3].z ) < fZCut ) {
-             const float det = (GetX() - clusters2[id2].x)*(GetY() - clusters1[id1].y) - (GetY() - clusters2[id2].y)*(GetX() - clusters1[id1].x);
+             //cout << "OK Z doublets: "<< iL-1 << "," << iD1 << "\t" << iL << "," <<iD2 << endl;
+             const float det = (clusters1[id1].x - GetX())*(clusters2[id2].y - GetY()) - (clusters2[id2].x-GetX() )*(clusters1[id1].y - GetY()); // (GetX() - clusters2[id2].x)*(clusters1[id1].y - clusters2[id2].y) - (GetY() - clusters2[id2].y)*(clusters1[id1].x - clusters2[id2].x);
+             //cout << det << endl;
              if ( Abs(det) <= 1e-12 ) {
                // linear extrapolation to next layer
                const float dsq = ( doublets1[iD1].tanPhi * (clusters3[id3].x + clusters2[id2].x) + clusters3[id3].y - clusters2[id2].y ) * 
@@ -199,14 +254,15 @@ void AliITSUTrackerSA::CASelection() {
              } else {
                const float r1sq = clusters1[id1].x * clusters1[id1].x + clusters1[id1].y * clusters1[id1].y ;
                const float rvsq = GetX() * GetX() + GetY() * GetY();
-               const float deta = (r2*r2-rvsq) * (GetY() - clusters1[id1].y) - (r1sq-rvsq) * (GetY() - clusters2[id2].y);
-               const float detb =  (r2*r2-rvsq) * (GetX() - clusters1[id1].x) - (r1sq-rvsq) * (GetX() - clusters2[id2].x) ;
+               const float deta = (rvsq - r1sq) * (clusters2[id2].y - GetY()) - (rvsq - r2*r2) * (clusters1[id1].y - GetY());
+               const float detb = - (rvsq - r1sq) * (clusters2[id2].x - GetX()) + (rvsq - r2*r2) * (clusters1[id1].x - GetX()) ;
                const float a = deta/det ;
                const float b = detb/det ;
                const float c = -rvsq - a * GetX() - b * GetY();
                const float rc = Sqrt( a*a/4.f + b*b/4.f - c );
                const float d = Sqrt( (a/2.f + clusters3[id3].x) * (a/2.f + clusters3[id3].x) + (b/2.f + clusters3[id3].y) * (b/2.f + clusters3[id3].y) );
-               if ( ( d - rc ) < fRPhiCut ) {
+               //cout << d << " " << rc << " " << d - rc << endl;
+               if ( Abs( d - rc ) < fRPhiCut ) {
                  doublets2[iD2].level += doublets1[iD1].level;
                  doublets2[iD2].neighbours.push_back(iD1);
                }
@@ -218,17 +274,96 @@ void AliITSUTrackerSA::CASelection() {
     }
   }
 
-  for ( int level = 6; level >=3 ; --level ) {
-    //vector<int> points[6];
-    for ( int doubl = 5; doubl >= 0; --doubl ) {
-      if( ( doubl + 1 - level ) < 0 ) break;
+  /*#ifdef __DEBUG__
+  PrintInfo("doublets");
+  #endif*/
+
+  // Hic sunt leones: the following code could be optimised to be iterative. But now I don't have time.
+  for ( int level = 6; level >= 2 ; --level ) {
+    vector<trackC> candidates;
+    int nCand=0;
+    const int limit = level > 5 ? 5 : level;
+    for ( int doubl = 5; doubl >= limit; --doubl ) {
       for ( int iD = 0; iD < fNDoublets[doubl]; ++iD ) {
        if ( fDoublets[doubl][iD].level == level ) {
-         
+         cout << "Pushing new candidate" << endl;
+         candidates.push_back(trackC());
+         candidates[nCand].fPoints[doubl*2+2] = fDoublets[doubl][iD].id1;
+         candidates[nCand].fPoints[doubl*2] = fDoublets[doubl][iD].id0;
+         CandidatesTreeTraversal(candidates, iD, doubl, nCand);
+       }
+      }
+    }
+
+    nCand = candidates.size();
+    Double_t rDest = fITS->GetRMax();
+    cout << "Number of candidates at level " << level << ": " << candidates.size() << endl;
+    int index[nCand];
+    for( int i = 0; i < nCand; ++i ) index[i] = i;
+    for ( int cand = 0; cand < nCand; ++cand ) {
+      int count = 0;
+      for ( int j = 0; j < 14; ++j ) {
+       if( candidates[cand].fPoints[j] > -1 ) ++count; 
+      }
+      cout << "Candidates " << cand << ", number of points: " << count << endl;
+      //      candidates[cand].ResetCovariance(50.);
+      InitTrackParams(candidates[cand]);
+      cout << "Fit cnd: " << cand << " " << candidates[cand] << endl;
+      candidates[cand].fChi2 = RefitTrack( (AliExternalTrackParam*)&candidates[cand], candidates[cand].fPoints, rDest ,-1);
+    }
+    CompDesc comp(&candidates);
+    sort(index,index+nCand,comp);
+    for ( int cand = 0; cand < nCand; ++cand ) {
+      const int ii = index[cand];
+#ifdef __DEBUG__
+      //cout << ii << " " << candidates[ii] << endl;
+#endif
+      if ( candidates[ii].fChi2 < 0. ) break;
+      bool goodTrack = true;
+      for ( unsigned int point = 0; point < 14; ++point ) {
+       if ( candidates[ii].fPoints[point] != -1 ) {
+         if( !(fClusters[ point/2 ][ candidates[ii].fPoints[point] ].isUsed ) ) {
+           fClusters[ point/2 ][ candidates[ii].fPoints[point] ].isUsed = true;
+         } else {
+           goodTrack = false;
+         }
        }
       }
+      if ( goodTrack ) {
+       AliESDtrack outTrack;
+       outTrack.SetOuterParam((AliExternalTrackParam*)&candidates[ii],AliESDtrack::kITSpureSA);
+       event->AddTrack(&outTrack);
+      }
     }
+
   }
+
+  // A first try... 
+  // for ( int level = 6; level >=3 ; --level ) {
+  //   vector<int> points[7];
+  //   for ( int doubl = 5; doubl >= 2; --doubl ) {
+  //     if( ( doubl + 1 - level ) < 0 ) break;
+  //     for ( int iD = 0; iD < fNDoublets[doubl]; ++iD ) {
+  //   if ( fDoublets[doubl][iD].level == level ) {
+  //     points[doubl+1].push_back( fDoublets[doubl][iD].id1 );
+  //     points[doubl].push_back( fDoublets[doubl][iD].id0 );
+     
+  //     vector<int> currentVector = (fDoublets[doubl][iD].neighbours);
+  //     for( int iL = 1; iL <= doubl; ++iL ) {
+  //       const nPlets * doublets = &fDoublets[doubl-iL][0];
+  //       vector<int> tmp;
+  //       for ( unsigned int iV = 0 ; iV < currentVector.size() ; ++iV ) {
+  //         points[doubl-iL].push_back( doublets[ currentVector[ iV ] ].id0 ) ;
+  //         for ( unsigned int iN = 0 ; iN < doublets[ currentVector[ iV ] ].neighbours.size(); ++iN ) {
+  //           tmp.push_back( doublets[ currentVector[ iV ] ].neighbours[ iN ] );
+  //         }
+  //       }
+  //       currentVector.swap( tmp );
+  //     }
+  //   }
+  //     }
+  //   }
+  // }
 }
 
 //________________________________________________________________________________
@@ -239,20 +374,24 @@ void AliITSUTrackerSA::MakeDoublets() {
   // To do:
   // - last iteration
 
+  cout << "Vertex of used by the tracker: " << GetX() << " " << GetY() << " " << GetZ() << endl;
+
   for( int iL = 0 ; iL < 6 ; ++iL ) {
     fNDoublets[iL] = 0; 
     const itsCluster* clusters1 = &fClusters[iL][0];
     const itsCluster* clusters2 = &fClusters[iL+1][0];
 
     // 0 - 2Pi junction treatment (part I)
-    for ( int iC1 = 0 ; iC1 < fNClusters[iL] ; ++iC1 ) {
+    for ( int iCC1 = 0 ; iCC1 < fNClusters[iL] ; ++iCC1 ) {
       bool flag = true;
-      for ( int iC2 = fNClusters[iL]-1; iC2 > -1 ; --iC2 ) {
-       
+      const int iC1 = fIndex[iL][iCC1];
+      for ( int iCC2 = fNClusters[iL+1]-1; iCC2 > -1 ; --iCC2 ) {
+       const int iC2 = fIndex[iL+1][iCC2];
        if( (TMath::TwoPi() - (clusters2[iC2].phi-clusters1[iC1].phi) ) < fPhiCut ) {
          fDoublets[iL].push_back(nPlets(iC1,iC2));
          fDoublets[iL][fNDoublets[iL]].tanPhi = (clusters1[iC1].y-clusters2[iC2].y)/(clusters1[iC1].x-clusters2[iC2].x);
          float r1  = Sqrt(clusters1[iC1].x * clusters1[iC1].x + clusters1[iC1].y * clusters1[iC1].y);
+         cout << clusters2[iC2].x * clusters2[iC2].x + clusters2[iC2].y * clusters2[iC2].y << flush << endl;
          float r2  = Sqrt(clusters2[iC2].x * clusters2[iC2].x + clusters2[iC2].y * clusters2[iC2].y);
          fDoublets[iL][fNDoublets[iL]].tanLambda = (clusters1[iC1].z-clusters2[iC2].z)/(r1-r2);
          ++fNDoublets[iL];
@@ -265,10 +404,10 @@ void AliITSUTrackerSA::MakeDoublets() {
 
     
     // "Central" points 
-    for ( int iC1 = 0 ; iC1 < fNClusters[iL] ; ++iC1 ) {
-
-      for ( int iC2 = 0; iC2 < fNClusters[iL+1] ; ++iC2 ) {
-       
+    for ( int iCC1 = 0 ; iCC1 < fNClusters[iL] ; ++iCC1 ) {
+      const int iC1 = fIndex[iL][iCC1];
+      for ( int iCC2 = 0; iCC2 < fNClusters[iL+1] ; ++iCC2 ) {
+       const int iC2 = fIndex[iL+1][iCC2];
        if( Abs( clusters1[iC1].phi - clusters2[iC2].phi ) < fPhiCut ) {
          fDoublets[iL].push_back(nPlets(iC1,iC2));
          fDoublets[iL][fNDoublets[iL]].tanPhi = (clusters1[iC1].y-clusters2[iC2].y)/(clusters1[iC1].x-clusters2[iC2].x);
@@ -283,11 +422,11 @@ void AliITSUTrackerSA::MakeDoublets() {
     }
 
     // 0 - 2Pi junction treatment (part II)
-    for ( int iC1 = fNClusters[iL]-1; iC1 > -1 ; --iC1 ) {
+    for ( int iCC1 = fNClusters[iL]-1; iCC1 > -1 ; --iCC1 ) {
       bool flag = true;
-
-      for ( int iC2 = 0; iC2 < fNClusters[iL+1] ; ++iC2 ) {
-
+      const int iC1 = fIndex[iL][iCC1];
+      for ( int iCC2 = 0; iCC2 < fNClusters[iL+1] ; ++iCC2 ) {
+       const int iC2 = fIndex[iL+1][iCC2];
        if( (TMath::TwoPi() - (clusters1[iC1].phi-clusters2[iC2].phi) ) < fPhiCut ) { 
          fDoublets[iL].push_back(nPlets(iC1,iC2));
          fDoublets[iL][fNDoublets[iL]].tanPhi = (clusters1[iC1].y-clusters2[iC2].y)/(clusters1[iC1].x-clusters2[iC2].x);
@@ -304,5 +443,467 @@ void AliITSUTrackerSA::MakeDoublets() {
     }
   
   }
+  // #ifdef __DEBUG__
+  // PrintInfo("doublets");
+  // #endif
+}
+//______________________________________________________________________________
+Bool_t AliITSUTrackerSA::InitTrackParams(trackC &track)
+{
+  // Set the initial guess on track kinematics for propagation.
+  // Assume at least 3 points available
+  int lrOcc[AliITSUAux::kMaxLayers], nCl=0;
+  //
+  // we will need endpoints and middle layer
+  for (int i=fITS->GetNLayersActive();i--;) if (track.fPoints[i<<0x1]>-1) lrOcc[nCl++] = i;
+  if (nCl<3) {
+    AliError(Form("Cannot estimate momentum of tracks with %d clusters",nCl));
+    cout << track << endl;
+    return kFALSE;
+  }
+  //
+  int lr0   = lrOcc[0];
+  int lr1   = lrOcc[nCl/2];
+  int lr2   = lrOcc[nCl-1];
+  //
+  const itsCluster& cl0 = fClusters[lr0][ track.fPoints[lr0<<0x1] ];
+  const itsCluster& cl1 = fClusters[lr1][ track.fPoints[lr1<<0x1] ];
+  const itsCluster& cl2 = fClusters[lr2][ track.fPoints[lr2<<0x1] ];
+  double cv = Curvature(cl0.x,cl0.y, cl1.x,cl1.y, cl2.x,cl2.y);
+  double tgl = (cl2.z-cl0.z)/TMath::Sqrt((cl2.x-cl0.x)*(cl2.x-cl0.x)+(cl2.y-cl0.y)*(cl2.y-cl0.y));
+  //  double phi = TMath::ATan2((cl2.y-cl1.y),(cl2.x-cl1.x));
+  //
+  AliITSUClusterPix* clus = (AliITSUClusterPix*)fClustersTC[ lr0 ]->At( track.fPoints[lr0<<0x1] );
+  AliITSURecoLayer* lr = fITS->GetLayerActive(lr0);
+  AliITSURecoSens* sens = lr->GetSensorFromID(clus->GetVolumeId());
+  double x = sens->GetXTF() + clus->GetX();
+  double alp = sens->GetPhiTF();
+  //  printf("Alp: %f phi: %f\n",alp,phi);
+  double par[5] = {clus->GetY(),clus->GetZ(),0,tgl,cv};
+  double cov[15] = {
+    5*5,
+    0, 5*5,
+    0, 0, 0.7*0.7,
+    0,0,0,0.7*0.7,
+    0,0,0,0,10
+  };
+  track.Set(x,alp,par,cov);
+  return kTRUE;
+}
+
+//______________________________________________________________________________
+void AliITSUTrackerSA::CandidatesTreeTraversal(vector<trackC> &track,  int &iD, int &doubl, int &nCand) {
+  
+  if ( doubl < 1 ) {
+#ifdef __DEBUG__
+    cout << "ERROR IN CandidatesTreeTraversal" << endl;
+#endif
+    return;
+  }
+
+  int currD = doubl-1;
+  track[nCand].fPoints[ 2 * currD ] = fDoublets[ currD ] [ fDoublets[ doubl ][ iD ].neighbours[0] ].id0;
+
+
+  if ( fDoublets[ currD ] [ fDoublets[ doubl ][ iD ].neighbours[0] ].level == 1 ) {
+
+    cout << "Setting track " << nCand << " information using its first cluster: " << endl;
+    cout << "Layer " << currD << ", id " << track[nCand].fPoints[ 2 * currD ] <<endl;
+    AliITSUClusterPix* clus = (AliITSUClusterPix*)fClustersTC[ currD ]->At( track[nCand].fPoints[ 2 * currD ] ); // assign your 1st cluster (from which the fit is starting)
+    AliITSURecoLayer* lr = fITS->GetLayerActive(currD) ; // assign the layer which the cluster belongs to
+    
+    AliITSURecoSens* sens = lr->GetSensorFromID(clus->GetVolumeId());
+    double x = sens->GetXTF() + clus->GetX();
+    double alp = sens->GetPhiTF();
+    double par[5] = {clus->GetY(),clus->GetZ(),0,0,0.01};
+    double cov[15] = {
+      5*5,
+      0, 5*5,
+      0, 0, 0.7*0.7,
+      0,0,0,0.7*0.7,
+      0,0,0,0,10
+    };
+    track[nCand].Set(x,alp,par,cov);
+    //    cout << "after set " << track[nCand] << endl;
+    return;
+  } else {
+    CandidatesTreeTraversal(track,fDoublets[doubl][iD].neighbours[0],currD,nCand);
+  }
+  
+  for( unsigned int iD2 = 1 ; iD2 < fDoublets[ doubl ][ iD ].neighbours.size() ; ++iD2 ) {
+    track.push_back(track[nCand++]);
+    track[nCand].fPoints[ 2 * currD ] = fDoublets[ currD ] [ fDoublets[ doubl ][ iD ].neighbours[ iD2 ] ].id0;
+    if ( fDoublets[ currD ] [ fDoublets[ doubl ][ iD ].neighbours[iD2] ].level == 1 ) {
+      AliITSUClusterPix* clus = (AliITSUClusterPix*)fClustersTC[ currD ]->At( track[nCand].fPoints[ 2 * currD ] ); // assign your 1st cluster (from which the fit is starting)
+      AliITSURecoLayer* lr = fITS->GetLayerActive(currD) ; // assign the layer which the cluster belongs to
+      
+      AliITSURecoSens* sens = lr->GetSensorFromID(clus->GetVolumeId());
+      double x = sens->GetXTF() + clus->GetX();
+      double alp = sens->GetPhiTF();
+      double par[5] = {clus->GetY(),clus->GetZ(),0,0,0};
+      double cov[15] = {
+       5*5,
+       0, 5*5,
+       0, 0, 0.7*0.7,
+       0,0,0,0.7*0.7,
+       0,0,0,0,10
+      };
+      track[nCand].Set(x,alp,par,cov);
+      return;
+    } else {
+      CandidatesTreeTraversal(track,fDoublets[doubl][iD].neighbours[iD2],currD,nCand);
+    }
+  }
+}
+
+//______________________________________________________________________________
+Double_t AliITSUTrackerSA::RefitTrack(AliExternalTrackParam* trc, 
+                                     Int_t clInfo[2*AliITSUAux::kMaxLayers],
+                                     Double_t rDest, Int_t stopCond)
+{
+  // refit track till radius rDest. 
+  // if stopCond<0 : propagate till last cluster then stop
+  // if stopCond==0: propagate till last cluster then try to go till limiting rDest, don't mind if fail
+  // if stopCond>0 : rDest must be reached
+  //
+  // The clList should provide the indices of clusters at corresponding layer (as stored in the layer
+  // TClonesArray, with convention (allowing for up to 2 clusters per layer due to the overlaps):
+  // if there is a cluster on given layer I, then it should be stored at clInfo[2*I-1]
+  // if there is an additional cluster on this layer, it goes to clInfo[2*I]
+  // -1 means no cluster
+  //
+  double rCurr = Sqrt(trc->GetX()*trc->GetX() + trc->GetY()*trc->GetY());
+  int dir,lrStart,lrStop;
+  //
+  dir = rCurr<rDest ? 1 : -1;
+  lrStart = fITS->FindFirstLayerID(rCurr,dir);
+  lrStop  = fITS->FindLastLayerID(rDest,dir); // lr id before which we have to stop
+  cout << "Start/End layers and direction " << lrStart << " " << lrStop << " " <<  dir << endl;
+  //
+  if (lrStop<0 || lrStart<0) AliFatal(Form("Failed to find start(%d) or last(%d) layers. "
+                                          "Track from %.3f to %.3f",lrStart,lrStop,rCurr,rDest));
+  //
+  int nCl = 0;
+  for (int i=2*fITS->GetNLayersActive();i--;) {if (clInfo[i]<0) continue; nCl++;}
+  //
+  AliExternalTrackParam tmpTr(*trc);
+  double chi2 = 0;
+  int iclLr[2],nclLr;
+  int nclFit = 0;
+  //
+  int lrStop1 = lrStop+dir;
+  for (int ilr=lrStart;ilr!=lrStop1;ilr+=dir) {
+    AliITSURecoLayer* lr = fITS->GetLayer(ilr);
+    if ( dir*(rCurr-lr->GetR(dir))>0) {
+      cout << ilr << " passed!" << endl;
+      continue;
+    } // this layer is already passed
+    int ilrA2,ilrA = lr->GetActiveID();
+    // passive layer or active w/o hits will be traversed on the way to next cluster
+    if (!lr->IsActive() || clInfo[ilrA2=(ilrA<<1)]<0) {
+      cout << ilr << " is inactive or without cluster for current candidates" << endl;
+      continue;
+    } 
+    cout << "OK layer " << ilr << endl;
+    //
+    // select the order in which possible 2 clusters (in case of the overlap) will be traversed and fitted
+    nclLr=0;
+    if (dir>0) { // clusters are stored in increasing radius order
+      iclLr[nclLr++]=clInfo[ilrA2++];
+      if (clInfo[ilrA2]>=0) iclLr[nclLr++]=clInfo[ilrA2];
+    }
+    else {
+      if ( clInfo[ilrA2+1]>=0 ) iclLr[nclLr++]=clInfo[ilrA2+1];
+      iclLr[nclLr++]=clInfo[ilrA2];
+    }
+    //
+    Bool_t transportedToLayer = kFALSE;
+    for (int icl=0;icl<nclLr;icl++) {
+      AliITSUClusterPix* clus =  (AliITSUClusterPix*)lr->GetCluster(iclLr[icl]);
+      AliITSURecoSens* sens = lr->GetSensorFromID(clus->GetVolumeId());
+      if (!tmpTr.Rotate(sens->GetPhiTF())) { cout << "failed rotation" << endl; return -1; }
+      //
+      double xClus = sens->GetXTF()+clus->GetX();
+      if (!transportedToLayer) {
+       if (ilr!=lrStart && !TransportToLayerX(&tmpTr,lrStart,ilr,xClus)) {
+         cout << "failed transport to the entrance" << endl; return -1; } // go to the entrance to the layer
+       lrStart = ilr;
+       transportedToLayer = kTRUE;
+      }
+      //
+      if (!PropagateSeed(&tmpTr,xClus,fCurrMass)) { cout << "failed propagation of the seed X:" << xClus << endl; tmpTr.Print(); return -1; }
+      //
+      Double_t p[2]={clus->GetY(), clus->GetZ()};
+      Double_t cov[3]={clus->GetSigmaY2(), clus->GetSigmaYZ(), clus->GetSigmaZ2()};
+      double chi2cl = tmpTr.GetPredictedChi2(p,cov);
+      chi2 += chi2cl;
+      //
+      if ( !tmpTr.Update(p,cov) ) { cout << "failed update of the covariance" << endl; return -1; }
+      if (++nclFit==nCl && stopCond<0) {
+       *trc = tmpTr;
+       //      printf("Fit chi2: %f for %d clusters\n",chi2,nclFit);
+       return chi2; // it was requested to not propagate after last update
+      }
+    }
+    //
+  }
+  // All clusters were succesfully fitted. Even if the track does not reach rDest, this is enough to validate it.
+  // Still, try to go as close as possible to rDest.
+  //
+  //  printf("Fit chi2: %f for %d clusters\n",chi2,nclFit);
+  //
+  if (lrStart!=lrStop) {
+    if (!TransportToLayer(&tmpTr,lrStart,lrStop)) return (stopCond>0) ? -chi2 : chi2; // rDest was obligatory
+    if (!GoToExitFromLayer(&tmpTr,fITS->GetLayer(lrStop),dir)) return (stopCond>0) ? -chi2 : chi2; // rDest was obligatory
+  }
+  // go to the destination radius. Note that here we don't select direction to avoid precision problems
+  if (!tmpTr.GetXatLabR(rDest,rDest,GetBz(),0) || !PropagateSeed(&tmpTr,rDest,fCurrMass, 100, kFALSE)) {
+    return (stopCond>0) ? -chi2 : chi2; // rDest was obligatory
+  }
+  *trc = tmpTr;
+
+  return chi2;
+}
 
+//______________________________________________________________________________
+Bool_t AliITSUTrackerSA::PropagateSeed(AliExternalTrackParam *seed, Double_t xToGo, Double_t mass, Double_t maxStep, Bool_t matCorr) 
+{
+  // propagate seed to given x applying material correction if requested
+  const Double_t kEpsilon = 1e-5;
+  Double_t xpos     = seed->GetX();
+  Int_t dir         = (xpos<xToGo) ? 1:-1;
+  Double_t xyz0[3],xyz1[3];
+  //
+  Bool_t updTime = dir>0 && seed->IsStartedTimeIntegral();
+  if (matCorr || updTime) seed->GetXYZ(xyz1);   //starting global position
+  while ( (xToGo-xpos)*dir > kEpsilon){
+    Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
+    Double_t x    = xpos+step;
+    Double_t bz=GetBz();   // getting the local Bz
+    if (!seed->PropagateTo(x,bz)) return kFALSE;
+    double ds = 0;
+    if (matCorr || updTime) {
+      xyz0[0]=xyz1[0]; // global pos at the beginning of step
+      xyz0[1]=xyz1[1];
+      xyz0[2]=xyz1[2];
+      seed->GetXYZ(xyz1);    //  // global pos at the end of step
+      //
+      if (matCorr) {
+       Double_t xrho,xx0;
+       ds = GetMaterialBudget(xyz0,xyz1,xx0,xrho);     
+       if (dir>0) xrho = -xrho; // outward should be negative
+       if (!seed->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
+      }
+      else { // matCorr is not requested but time integral is
+       double d0 = xyz1[0]-xyz0[0];
+       double d1 = xyz1[1]-xyz0[1];
+       double d2 = xyz1[2]-xyz0[2];    
+       ds = TMath::Sqrt(d0*d0+d1*d1+d2*d2);
+      }
+    }
+    if (updTime) seed->AddTimeStep(ds);
+    //
+    xpos = seed->GetX();
+  }
+  return kTRUE;
+}
+
+//_________________________________________________________________________
+Bool_t AliITSUTrackerSA::TransportToLayer(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo, Double_t rLim)
+{
+  // transport track from layerFrom to the entrance of layerTo or to rLim (if>0), wathever is closer
+  //  
+  if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
+  //
+  int dir = lTo > lFrom ? 1:-1;
+  AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
+  Bool_t checkFirst = kTRUE;
+  Bool_t limReached = kFALSE;
+  while(lFrom!=lTo) {
+    if (lrFr) {
+      if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) return kFALSE; // go till the end of current layer
+      checkFirst = kFALSE;
+    }
+    AliITSURecoLayer* lrTo =  fITS->GetLayer( (lFrom+=dir) );
+    if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
+    //
+    // go the entrance of the layer, assuming no materials in between
+    double xToGo = lrTo->GetR(-dir);
+    if (rLim>0) {
+      if (dir>0) {
+       if (rLim<xToGo) {xToGo = rLim; limReached = kTRUE;}
+      }
+      else {
+       if (rLim>xToGo) {xToGo = rLim; limReached = kTRUE;}
+      }
+    }
+    //    double xts = xToGo;
+    if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) {
+      //      printf("FailHere1: %f %f %d\n",xts,xToGo,dir);
+      //      seed->Print("etp");
+      return kFALSE;
+    }
+    if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) {
+      //printf("FailHere2: %f %f %d\n",xts,xToGo,dir);
+      //seed->Print("etp");
+      return kFALSE;
+    }
+    lrFr = lrTo;
+    if (limReached) break;
+  }
+  return kTRUE;
+  //
+}
+
+//_________________________________________________________________________
+Bool_t AliITSUTrackerSA::TransportToLayerX(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo, Double_t xStop)
+{
+  // transport track from layerFrom to the entrance of layerTo but do not pass control parameter X
+  //  
+  if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
+  //
+  int dir = lTo > lFrom ? 1:-1;
+  AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
+  Bool_t checkFirst = kTRUE;
+  while(lFrom!=lTo) {
+    if (lrFr) {
+      if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) return kFALSE; // go till the end of current layer
+      checkFirst = kFALSE;
+    }
+    AliITSURecoLayer* lrTo =  fITS->GetLayer( (lFrom+=dir) );
+    if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
+    //
+    // go the entrance of the layer, assuming no materials in between
+    double xToGo = lrTo->GetR(-dir); // R of the entrance to layer
+    //
+    //    double xts = xToGo;
+    if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) {
+      //      printf("FailHere1: %f %f %d\n",xts,xToGo,dir);
+      //      seed->Print("etp");
+      return kFALSE;
+    }
+    if ( (dir>0&&xToGo>xStop) || (dir<0&&xToGo<xStop) ) xToGo = xStop;
+    //
+#ifdef _ITSU_DEBUG_
+    AliDebug(2,Form("go in dir=%d to R=%.4f(X:%.4f)",dir,lrTo->GetR(-dir), xToGo));
+#endif
+    if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) {
+      //printf("FailHere2: %f %f %d\n",xts,xToGo,dir);
+      //seed->Print("etp");
+      return kFALSE;
+    }
+    lrFr = lrTo;
+  }
+  return kTRUE;
+  //
 }
+
+//_________________________________________________________________________
+Bool_t AliITSUTrackerSA::GoToExitFromLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
+{
+  // go to the exit from lr in direction dir, applying material corrections in steps specific for this layer
+  // If check is requested, do this only provided the track has not exited the layer already
+  double xToGo = lr->GetR(dir);
+  if (check) { // do we need to track till the surface of the current layer ?
+    double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
+    if      (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // on the surface or outside of the layer
+    else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // on the surface or outside of the layer
+  }
+  if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
+  // go via layer to its boundary, applying material correction.
+  if (!PropagateSeed(seed,xToGo,fCurrMass, lr->GetMaxStep())) return kFALSE;
+  //
+  return kTRUE;
+  //
+}
+
+//_________________________________________________________________________
+Bool_t AliITSUTrackerSA::GoToEntranceToLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
+{
+  // go to the entrance of lr in direction dir, w/o applying material corrections.
+  // If check is requested, do this only provided the track did not reach the layer already
+  double xToGo = lr->GetR(-dir);
+  if (check) { // do we need to track till the surface of the current layer ?
+    double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
+    if      (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // already passed
+    else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // already passed
+  }
+  if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
+  // go via layer to its boundary, applying material correction.
+  if (!PropagateSeed(seed,xToGo,fCurrMass, 100, kFALSE)) return kFALSE;
+  return kTRUE;
+  //
+}
+
+//____________________________________________________
+Double_t AliITSUTrackerSA::GetMaterialBudget(const double* pnt0,const double* pnt1, double& x2x0, double& rhol) const
+{
+  double par[7];
+  if (fUseMatLUT && fMatLUT) {
+    double d = fMatLUT->GetMatBudget(pnt0,pnt1,par);
+    x2x0 = par[AliITSUMatLUT::kParX2X0];
+    rhol = par[AliITSUMatLUT::kParRhoL];    
+    return d;
+  }
+  else {
+    MeanMaterialBudget(pnt0,pnt1,par);
+    x2x0 = par[1];
+    rhol = par[0]*par[4];    
+    return par[4];
+  }
+}
+
+//____________________________________________________________________
+Double_t AliITSUTrackerSA::Curvature(Double_t x1,Double_t y1,Double_t 
+x2,Double_t y2,Double_t x3,Double_t y3)
+{
+
+  //calculates the curvature of track  
+  Double_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
+  if(den==0) return 0;
+  Double_t a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
+  Double_t b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
+  Double_t c = -x1*x1-y1*y1-a*x1-b*y1;
+  Double_t xc=-a/2.;
+
+  if((a*a+b*b-4*c)<0) return 0;
+  Double_t rad = TMath::Sqrt(a*a+b*b-4*c)/2.;
+  if(rad==0) return 0;
+  
+  if((x1>0 && y1>0 && x1<xc)) rad*=-1;
+  if((x1<0 && y1>0 && x1<xc)) rad*=-1;
+  //  if((x1<0 && y1<0 && x1<xc)) rad*=-1;
+  // if((x1>0 && y1<0 && x1<xc)) rad*=-1;
+  
+  return 1/rad;
+}
+
+#ifdef __DEBUG__ 
+//____________________________________________________
+void AliITSUTrackerSA::PrintInfo(TString what) {
+  //
+  if( what.Contains("clusters") ) {
+    cout << "Dumping clusters info" << endl;
+    for ( int i = 0; i < 7; ++i ) {
+      cout << "**** Layer " << i << " ****" << endl;
+      for ( int c = 0; c < fNClusters[i]; ++c ) {
+       cout << "*** Cluster " << c << " ***" <<endl;
+       cout << fClusters[i][c] << endl;
+      }
+    }
+  }
+  //
+  if( what.Contains("doublets") ) {
+    cout << "Dumping doublets info" << endl;
+    for ( int i = 0; i < 6; ++i ) {
+      cout << "**** Doublets array " << i << " ****" << endl;
+      for ( int c = 0; c < fNDoublets[i]; ++c ) {
+       cout << "*** Doublet " << c << " ***" <<endl;
+       cout << fDoublets[i][c] << endl;
+      }
+    }
+  }
+}
+
+#endif
old mode 100755 (executable)
new mode 100644 (file)
index 48834b1..94c9d24
@@ -1,6 +1,11 @@
 #ifndef ALIITSUTRACKERSA_H
 #define ALIITSUTRACKERSA_H
 
+#define __DEBUG__
+#ifdef __DEBUG__ 
+#include <TString.h>
+#endif
+
 //-------------------------------------------------------------------------
 //                   The stand-alone ITSU tracker
 //     It reads AliITSUClusterPix clusters and writes the tracks to the ESD
 #include <TClonesArray.h>
 #include <vector>
 #include "AliITSUTrackerSAaux.h"   // Structs and other stuff 
+#include "AliITSUMatLUT.h"
+#include "AliITSUAux.h"
+#include "AliExternalTrackParam.h"
+class AliITSUReconstructor;
+class AliITSURecoDet;
+class AliITSURecoLayer;
 
 class TTree;
 class AliCluster;
@@ -21,8 +32,8 @@ using std::vector;
 //-------------------------------------------------------------------------
 class AliITSUTrackerSA : public AliTracker {
 public:
-  AliITSUTrackerSA();
-  virtual ~AliITSUTrackerSA() {} ;
+  AliITSUTrackerSA(AliITSUReconstructor* rec=0);
+  virtual ~AliITSUTrackerSA();
 
   // These functions must be implemented 
   Int_t Clusters2Tracks(AliESDEvent *event);
@@ -33,26 +44,49 @@ public:
   AliCluster *GetCluster(Int_t index) const;
 
   // Possibly, other public functions
-
-
+  void     Init(AliITSUReconstructor* rec);
+  Double_t RefitTrack(AliExternalTrackParam* trc, Int_t clInfo[2*AliITSUAux::kMaxLayers], Double_t rDest, Int_t stopCond);  
+  Bool_t   PropagateSeed(AliExternalTrackParam *seed, Double_t xToGo, Double_t mass, Double_t maxStep=1.0, Bool_t matCorr=kTRUE);
+  Double_t GetMaterialBudget(const double* pnt0, const double* pnt1, double& x2x0, double& rhol) const;
+  Bool_t   GoToEntranceToLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check=kFALSE);
+  Bool_t   GoToExitFromLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check=kTRUE);
+  Bool_t   TransportToLayerX(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo, Double_t xStop);  
+  Bool_t   TransportToLayer(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo, Double_t rLim=-1);
+  //
 protected:
   AliITSUTrackerSA(const AliITSUTrackerSA&);
 
   void MakeDoublets();
-  void MakeTriplets();
-  void CASelection();
+  //  void MakeTriplets();
+  void CandidatesTreeTraversal( vector<trackC> &vec, int &iD, int &doubl, int &nCand);
+  Bool_t InitTrackParams(trackC &track);
+  void CASelection(AliESDEvent *ev);
   void GlobalFit();
   void ChiSquareSelection();
   // Other protected functions
-  // (Sorting, labeling, calculations of "roads", etc)
+  // (Sorting, labeling, calculations of "roads", etc)  
+  static Double_t Curvature(Double_t x1,Double_t y1,Double_t x2,Double_t y2,Double_t x3,Double_t y3);
+
+#ifdef __DEBUG__ 
+  void PrintInfo(TString opt);
+#endif
 
 private:
   AliITSUTrackerSA &operator=(const AliITSUTrackerSA &tr);
 
   // Data members
+
+  // classes for interfacing the geometry, materials etc.
+  AliITSUReconstructor*           fReconstructor;  // ITS global reconstructor 
+  AliITSURecoDet*                 fITS;            // interface to ITS, borrowed from reconstructor
+  AliITSUMatLUT*                  fMatLUT;         // material lookup table
+  Bool_t                          fUseMatLUT;      //! use material lookup table rather than TGeo
+  Double_t                        fCurrMass;       // assumption about particle mass
+
+
   // Internal tracker arrays, layers, modules, etc
   vector<itsCluster> fClusters[7];
-  TClonesArray fClustersTC[7];
+  TClonesArray *fClustersTC[7];
   vector<nPlets> fDoublets[6];
   Int_t *fIndex[7];
   Int_t fNClusters[7];
@@ -61,6 +95,9 @@ private:
   Float_t fRPhiCut;
   Float_t fZCut;
 
+  //
+  static const Double_t           fgkToler;        // tracking tolerance
+  //
   ClassDef(AliITSUTrackerSA,1)   //ITSU stand-alone tracker
 };
 
old mode 100755 (executable)
new mode 100644 (file)
index 9f1bb2a..d4c74b9
@@ -1,18 +1,16 @@
 #ifndef ALIITSUTRACKERSAAUX_H
 #define ALIITSUTRACKERSAAUX_H
 
+#ifdef __DEBUG__ 
+#include <iostream>
+using std::ostream;
+using std::endl;
+#endif
+
 #include <vector>
 using std::vector;
-
-/* struct CompareAsc { // Adapted from TMath ROOT code */
-/* CompareAsc(Float_t *d,Float_t offset) : fData(d) {} */
-  
-/*   bool operator()(int i1, int i2) { */
-/*     return *(fData + i1) < *(fData + i2); */
-/*   } */
-  
-/*   Float_t fData; */
-/* }; */
+#include "AliExternalTrackParam.h"
+#include "AliITSUAux.h"
 
 struct itsCluster {
 itsCluster():isUsed(false),x(0.f),y(0.f),z(0.f),varx(0.f),covxy(0.f),vary(0.f),phi(0.f),phiM(0.f) {}
@@ -22,15 +20,76 @@ itsCluster(const float &X,const float &Y, const float &Z, const float &varX, con
   float x,y,z;            // Global coordinates
   float varx,covxy,vary;  // Local covariance matrix
   float phi,phiM;         // phi of the cluster and phi angle of the module containing the cluster
+  
+#ifdef __DEBUG__
+  friend ostream& operator<<(ostream& out, const itsCluster& cl) {
+    out << "pos = (" << cl.x << ", " << cl.y << ", "<< cl.z <<")"<<" phi="<<cl.phi <<endl;
+    return out;
+  }
+#endif
 };
 
 struct nPlets {
-nPlets() : id0(-1),id1(-1),id2(-1),level(0),tanPhi(),tanLambda(),neighbours() {}
-nPlets(int arg0,int arg1) : id0(arg0),id1(arg1),id2(-1),level(0),tanPhi(),tanLambda(),neighbours()  {}
+nPlets() : id0(-1),id1(-1),id2(-1),level(1),tanPhi(),tanLambda(),neighbours() {}
+nPlets(int arg0,int arg1) : id0(arg0),id1(arg1),id2(-1),level(1),tanPhi(),tanLambda(),neighbours()  {}
   int id0,id1,id2;
   int level;
   float tanPhi,tanLambda;
   vector<int> neighbours;
+#ifdef __DEBUG__
+  friend ostream& operator<<(ostream& out, const nPlets& cl) {
+    out << "id = (" << cl.id0 << ", " << cl.id1 << ", "<< cl.id2 <<")"<< endl;
+    out << "tanPhi="<< cl.tanPhi <<" tanLambda="<<cl.tanLambda << " level=" << cl.level <<endl;
+    out << "neighbours= ";
+    for( unsigned int i = 0; i< cl.neighbours.size(); ++i ) out << cl.neighbours[i];
+    out << endl;
+    return out;
+  }
+#endif
 };
 
+class trackC : public AliExternalTrackParam {
+ public : 
+ trackC() : AliExternalTrackParam(), fChi2( 0. ), fPoints() {
+    for ( unsigned int i = 0; i < 2*AliITSUAux::kMaxLayers; ++i ) fPoints[i]=-1;
+  }
+#ifdef __DEBUG__
+  friend ostream& operator<<(ostream& out, const trackC& cl) {
+    out << "points = (";
+    for( unsigned int i = 0; i < 2*AliITSUAux::kMaxLayers; ++i ) out << cl.fPoints[i] << " ";
+    out << "), chi2: " << cl.fChi2 <<endl;
+    const double* par = cl.GetParameter();
+    const double* cov = cl.GetCovariance();    
+    out << "X: " << cl.GetX() << " Alpha: " << cl.GetAlpha() << endl;
+    out << "Param: \n";
+    for (int i=0;i<5;i++) out << par[i] << " "; out << endl;
+    out << "Covar: \n";
+    int cnt = 0;
+    for (int i=0;i<5;i++) {for (int j=i+1;j--;) out << cov[cnt++] << " "; out << endl;}
+    return out;
+  }
+#endif
+  
+  /* bool operator < ( const trackC &t ) { return fChi2 < t.fChi2; } */
+  /* bool operator <= ( const trackC &t ) { return fChi2 <= t.fChi2; } */
+  /* bool operator > ( const trackC &t ) { return fChi2 > t.fChi2; } */
+  /* bool operator >= ( const trackC &t ) { return fChi2 >= t.fChi2; } */
+  /* bool operator == ( const trackC &t ) { return fChi2 == t.fChi2; } */
+  /* bool operator != ( const trackC &t ) { return fChi2 != t.fChi2; } */
+  Double_t fChi2;
+  Int_t fPoints[2*AliITSUAux::kMaxLayers]; 
+  
+};
+
+
+struct CompDesc { //Adapted from TMath ROOT code 
+CompDesc(vector<trackC> *d) : fData(d) {} 
+  
+  bool operator()(int i1, int i2) { 
+    return fData->at(i1).fChi2 > fData->at(i2).fChi2; 
+  } 
+  
+  vector<trackC> *fData; 
+}; 
+
 #endif