ITS/UPGRADE/AliITSUTrackerSA.cxx

   1 //-------------------------------------------------------------------------
   2 //               Implementation of the ITS tracker class
   3 //    It reads AliITSUClusterPix clusters and and fills the ESD with tracks
   4 //-------------------------------------------------------------------------
   5
   6 #include <TBranch.h>
   7 #include <TMath.h>
   8 using TMath::Abs;
   9 using TMath::Sort;
  10 using TMath::Sqrt;
  11 #include <TTree.h>
  12
  13 // Vc library
  14 //#include "Vc/Vc"
  15 //#include "AliITSUTrackerSAauxVc.h" // Structs and other stuff using Vc library
  16
  17 #include "AliESDEvent.h"
  18 #include "AliITSUClusterPix.h"
  19 #include "AliITSUTrackerSA.h"
  20
  21 //#include "AliITSUtrackSA.h"      // Some dedicated SA track class ?
  22
  23 ClassImp(AliITSUTrackerSA)
  24
  25 //________________________________________________________________________________
  26 AliITSUTrackerSA::AliITSUTrackerSA() : AliTracker(),
  27   fClusters(),
  28   fClustersTC(),
  29   fDoublets(),
  30   fIndex(),
  31   fNClusters(),
  32   fNDoublets(),
  33   fPhiCut(0.05),
  34   fRPhiCut(0.01),
  35   fZCut(0.005)
  36 {
  37   //--------------------------------------------------------------------
  38   // This default constructor needs to be provided
  39   //--------------------------------------------------------------------
  40   for(Int_t i=0;i<7;++i) {
  41     fClusters[i].reserve(5000);
  42   }
  43 }
  44
  45 //________________________________________________________________________________
  46 AliITSUTrackerSA::AliITSUTrackerSA(const AliITSUTrackerSA &t):
  47   AliTracker(t),
  48   fClusters(),
  49   fClustersTC(),
  50   fIndex(),
  51   fNClusters(),
  52   fNDoublets(),
  53   fPhiCut(),
  54   fRPhiCut(),
  55   fZCut()
  56 {
  57   //--------------------------------------------------------------------
  58   // The copy constructor is protected
  59   //--------------------------------------------------------------------
  60 }
  61
  62 //________________________________________________________________________________
  63 Int_t AliITSUTrackerSA::Clusters2Tracks(AliESDEvent */*event*/) {
  64   //--------------------------------------------------------------------
  65   // This is the main tracking function
  66   // The clusters must already be loaded
  67   //--------------------------------------------------------------------
  68
  69   // Possibly, create the track "seeds" (combinatorial)
  70
  71   // Possibly, increment the seeds with additional clusters (Kalman)
  72
  73   // Possibly, (re)fit the found tracks
  74
  75   // Tree iterations:
  76   // - High momentum first;
  77   // - Low momentum with vertex constraint;
  78   // - Everything else.
  79
  80   MakeDoublets();       // To be implemented
  81   //MakeTriplets();       // Are triplets really necessary? MFT does not use them.
  82   CASelection();        // To be implemented
  83   GlobalFit();          // To be implemented
  84   ChiSquareSelection(); // To be implemented
  85
  86   return 0;
  87 }
  88
  89 //________________________________________________________________________________
  90 Int_t AliITSUTrackerSA::PropagateBack(AliESDEvent */*event*/) {
  91   //--------------------------------------------------------------------
  92   // Here, we implement the Kalman smoother ?
  93   // The clusters must be already loaded
  94   //--------------------------------------------------------------------
  95
  96   return 0;
  97 }
  98
  99 //________________________________________________________________________________
 100 Int_t AliITSUTrackerSA::RefitInward(AliESDEvent */*event*/) {
 101   //--------------------------------------------------------------------
 102   // Some final refit, after the outliers get removed by the smoother ?
 103   // The clusters must be loaded
 104   //--------------------------------------------------------------------
 105
 106   return 0;
 107 }
 108
 109 Int_t AliITSUTrackerSA::LoadClusters(TTree *cluTree) {
 110   //--------------------------------------------------------------------
 111   // This function reads the ITSU clusters from the tree,
 112   // sort them, distribute over the internal tracker arrays, etc
 113   //--------------------------------------------------------------------
 114
 115   for(Int_t i=0;i<7;++i) {
 116     TBranch* br = cluTree->GetBranch(Form("ITSRecPoints%d",i));
 117     if (!br) return -1;
 118     br->SetAddress(&fClustersTC[i]);
 119   }
 120   cluTree->GetEntry(0);
 121
 122   for(int iL=0; iL<7; ++iL) {
 123     TClonesArray *clCont=&fClustersTC[iL];
 124     fNClusters[iL]=clCont->GetEntriesFast();
 125     Float_t phi[fNClusters[iL]];
 126     fIndex[iL] = new Int_t[fNClusters[iL]];
 127     for(int iC=0;iC<fNClusters[iL];++iC) {
 128       const AliITSUClusterPix *cl = (AliITSUClusterPix*)clCont->At(iC);
 129       float pos[3];
 130       cl->GetGlobalXYZ(pos);
 131       phi[iC] = pos[0]==0.f ? TMath::PiOver2() : TMath::ATan2(pos[1]-GetY(),pos[0]-GetX());
 132       float angle=0.f; // TO BE UNDERSTOOD: DO I STILL NEED THE STATION ANGLE IF I USE THE GLOBAL COVARIANCE MATRIX?
 133       fClusters[iL].push_back(itsCluster(pos[0],pos[1],pos[2],cl->GetSigmaY2(),cl->GetSigmaZ2(),cl->GetSigmaYZ(),phi[iC],angle));
 134     }
 135     TMath::Sort(fNClusters[iL],phi,fIndex[iL],kFALSE);
 136   }
 137
 138   return 0;
 139 }
 140
 141 //________________________________________________________________________________
 142 void AliITSUTrackerSA::UnloadClusters() {
 143   //--------------------------------------------------------------------
 144   // This function unloads ITSU clusters from the RAM
 145   //--------------------------------------------------------------------
 146   for(int i=0;i<7;++i) {
 147     fClusters[i].clear();
 148     fNClusters[i]=0;
 149     delete fIndex[i];
 150   }
 151   for(int i=0;i<6;++i) fDoublets[i].clear();
 152 }
 153
 154 //________________________________________________________________________________
 155 AliCluster *AliITSUTrackerSA::GetCluster(Int_t /*index*/) const {
 156   //--------------------------------------------------------------------
 157   //       Return pointer to a given cluster
 158   //--------------------------------------------------------------------
 159   return 0;  // replace with an actual pointer
 160 }
 161
 162 //________________________________________________________________________________
 163 void AliITSUTrackerSA::CASelection() {
 164   // Here it's implemented the Cellular Automaton routine
 165   // Firstly the level of each doublet is set according to the level of
 166   // the neighbour doublets.
 167   // Doublet are considered to be neighbour if they share one point and the
 168   // phi and theta direction difference of the two is below a cut value.
 169
 170   for( int iL = 1; iL < 6; ++iL ) {
 171
 172     const itsCluster* clusters1 = &fClusters[iL-1][0];
 173     const itsCluster* clusters2 = &fClusters[iL][0];
 174     const itsCluster* clusters3 = &fClusters[iL+1][0];
 175
 176     const nPlets* doublets1 = &fDoublets[iL-1][0];
 177     nPlets* doublets2 = &fDoublets[iL][0];
 178
 179     for ( int iD2 = 0; iD2 < fNDoublets[iL]; ++iD2 ) {
 180       for ( int iD1 = 0; iD1 < fNDoublets[iL-1]; ++iD1 ) {
 181         const int id1 = doublets1[iD1].id1;
 182         const int id2 = doublets2[iD2].id0;
 183         if ( id1 == id2 ) {
 184           if ( doublets2[iD2].level <= ( doublets1[iD1].level + 1 ) ) {
 185             const int id3 = doublets2[iD2].id1;
 186             const float r3 = Sqrt( clusters3[id3].x * clusters3[id3].x + clusters3[id3].y * clusters3[id3].y );
 187             const float r2 = Sqrt( clusters3[id2].x * clusters2[id2].x + clusters2[id2].y * clusters2[id2].y );
 188             const float extrZ3 = doublets1[iD1].tanLambda * ( r3 - r2 ) + clusters3[id3].z ;
 189             if ( Abs ( extrZ3 - clusters3[id3].z ) < fZCut ) {
 190               const float det = (GetX() - clusters2[id2].x)*(GetY() - clusters1[id1].y) - (GetY() - clusters2[id2].y)*(GetX() - clusters1[id1].x);
 191               if ( Abs(det) <= 1e-12 ) {
 192                 // linear extrapolation to next layer
 193                 const float dsq = ( doublets1[iD1].tanPhi * (clusters3[id3].x + clusters2[id2].x) + clusters3[id3].y - clusters2[id2].y ) *
 194                   ( doublets1[iD1].tanPhi * (clusters3[id3].x + clusters2[id2].x) + clusters3[id3].y - clusters2[id2].y ) / (1 + doublets1[iD1].tanPhi * doublets1[iD1].tanPhi );
 195                 if ( dsq < fRPhiCut*fRPhiCut )  {
 196                   doublets2[iD2].level += doublets1[iD1].level;
 197                   doublets2[iD2].neighbours.push_back(iD1);
 198                 }
 199               } else {
 200                 const float r1sq = clusters1[id1].x * clusters1[id1].x + clusters1[id1].y * clusters1[id1].y ;
 201                 const float rvsq = GetX() * GetX() + GetY() * GetY();
 202                 const float deta = (r2*r2-rvsq) * (GetY() - clusters1[id1].y) - (r1sq-rvsq) * (GetY() - clusters2[id2].y);
 203                 const float detb =  (r2*r2-rvsq) * (GetX() - clusters1[id1].x) - (r1sq-rvsq) * (GetX() - clusters2[id2].x) ;
 204                 const float a = deta/det ;
 205                 const float b = detb/det ;
 206                 const float c = -rvsq - a * GetX() - b * GetY();
 207                 const float rc = Sqrt( a*a/4.f + b*b/4.f - c );
 208                 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) );
 209                 if ( ( d - rc ) < fRPhiCut ) {
 210                   doublets2[iD2].level += doublets1[iD1].level;
 211                   doublets2[iD2].neighbours.push_back(iD1);
 212                 }
 213               }
 214             }
 215           }
 216         }
 217       }
 218     }
 219   }
 220
 221   for ( int level = 6; level >=3 ; --level ) {
 222     //vector<int> points[6];
 223     for ( int doubl = 5; doubl >= 0; --doubl ) {
 224       if( ( doubl + 1 - level ) < 0 ) break;
 225       for ( int iD = 0; iD < fNDoublets[doubl]; ++iD ) {
 226         if ( fDoublets[doubl][iD].level == level ) {
 227
 228         }
 229       }
 230     }
 231   }
 232 }
 233
 234 //________________________________________________________________________________
 235 void AliITSUTrackerSA::MakeDoublets() {
 236   // Make associations between two points on adjacent layers within an azimuthal window.
 237   // Under consideration:
 238   // - track parameter estimation using the primary vertex position
 239   // To do:
 240   // - last iteration
 241
 242   for( int iL = 0 ; iL < 6 ; ++iL ) {
 243     fNDoublets[iL] = 0;
 244     const itsCluster* clusters1 = &fClusters[iL][0];
 245     const itsCluster* clusters2 = &fClusters[iL+1][0];
 246
 247     // 0 - 2Pi junction treatment (part I)
 248     for ( int iC1 = 0 ; iC1 < fNClusters[iL] ; ++iC1 ) {
 249       bool flag = true;
 250       for ( int iC2 = fNClusters[iL]-1; iC2 > -1 ; --iC2 ) {
 251
 252         if( (TMath::TwoPi() - (clusters2[iC2].phi-clusters1[iC1].phi) ) < fPhiCut ) {
 253           fDoublets[iL].push_back(nPlets(iC1,iC2));
 254           fDoublets[iL][fNDoublets[iL]].tanPhi = (clusters1[iC1].y-clusters2[iC2].y)/(clusters1[iC1].x-clusters2[iC2].x);
 255           float r1  = Sqrt(clusters1[iC1].x * clusters1[iC1].x + clusters1[iC1].y * clusters1[iC1].y);
 256           float r2  = Sqrt(clusters2[iC2].x * clusters2[iC2].x + clusters2[iC2].y * clusters2[iC2].y);
 257           fDoublets[iL][fNDoublets[iL]].tanLambda = (clusters1[iC1].z-clusters2[iC2].z)/(r1-r2);
 258           ++fNDoublets[iL];
 259           flag = false;
 260         } else break;
 261
 262       }
 263       if (flag) break;
 264     }
 265
 266
 267     // "Central" points
 268     for ( int iC1 = 0 ; iC1 < fNClusters[iL] ; ++iC1 ) {
 269
 270       for ( int iC2 = 0; iC2 < fNClusters[iL+1] ; ++iC2 ) {
 271
 272         if( Abs( clusters1[iC1].phi - clusters2[iC2].phi ) < fPhiCut ) {
 273           fDoublets[iL].push_back(nPlets(iC1,iC2));
 274           fDoublets[iL][fNDoublets[iL]].tanPhi = (clusters1[iC1].y-clusters2[iC2].y)/(clusters1[iC1].x-clusters2[iC2].x);
 275           float r1  = Sqrt(clusters1[iC1].x * clusters1[iC1].x + clusters1[iC1].y * clusters1[iC1].y);
 276           float r2  = Sqrt(clusters2[iC2].x * clusters2[iC2].x + clusters2[iC2].y * clusters2[iC2].y);
 277           fDoublets[iL][fNDoublets[iL]].tanLambda = (clusters1[iC1].z-clusters2[iC2].z)/(r1-r2);
 278           ++fNDoublets[iL];
 279         } else if( clusters2[iC2].phi - clusters1[iC1].phi > fPhiCut ) break;
 280
 281       }
 282
 283     }
 284
 285     // 0 - 2Pi junction treatment (part II)
 286     for ( int iC1 = fNClusters[iL]-1; iC1 > -1 ; --iC1 ) {
 287       bool flag = true;
 288
 289       for ( int iC2 = 0; iC2 < fNClusters[iL+1] ; ++iC2 ) {
 290
 291         if( (TMath::TwoPi() - (clusters1[iC1].phi-clusters2[iC2].phi) ) < fPhiCut ) {
 292           fDoublets[iL].push_back(nPlets(iC1,iC2));
 293           fDoublets[iL][fNDoublets[iL]].tanPhi = (clusters1[iC1].y-clusters2[iC2].y)/(clusters1[iC1].x-clusters2[iC2].x);
 294           float r1  = Sqrt(clusters1[iC1].x * clusters1[iC1].x + clusters1[iC1].y * clusters1[iC1].y);
 295           float r2  = Sqrt(clusters2[iC2].x * clusters2[iC2].x + clusters2[iC2].y * clusters2[iC2].y);
 296           fDoublets[iL][fNDoublets[iL]].tanLambda = (clusters1[iC1].z-clusters2[iC2].z)/(r1-r2);
 297           ++fNDoublets[iL];
 298           flag = false;
 299         } else break;
 300
 301       }
 302
 303       if (flag) break;
 304     }
 305
 306   }
 307
 308 }