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 #include <algorithm>
  13 using std::sort;
  14
  15
  16 // Vc library
  17 //#include "Vc/Vc"
  18 //#include "AliITSUTrackerSAauxVc.h" // Structs and other stuff using Vc library
  19 #include "AliLog.h"
  20 #include "AliESDEvent.h"
  21 #include "AliITSUClusterPix.h"
  22 #include "AliITSUTrackerSA.h"
  23 #include "AliITSUReconstructor.h"
  24 #include "AliITSURecoDet.h"
  25 #include "AliESDtrack.h"
  26
  27 #include <Riostream.h>
  28
  29 using std::cout;
  30 using std::endl;
  31 using std::flush;
  32
  33 //#include "AliITSUtrackSA.h"      // Some dedicated SA track class ?
  34
  35 ClassImp(AliITSUTrackerSA)
  36
  37 const Double_t AliITSUTrackerSA::fgkToler =  1e-6;// tolerance for layer on-surface check
  38 const Double_t AliITSUTrackerSA::fgkChi2Cut =  10.f;
  39
  40 //________________________________________________________________________________
  41 AliITSUTrackerSA::AliITSUTrackerSA(AliITSUReconstructor* rec) :
  42 fReconstructor(rec),
  43 fITS(0),
  44 fMatLUT(0),
  45 fUseMatLUT(kFALSE),
  46 fCurrMass(0.14),
  47 //
  48 fClusters(),
  49 fClustersTC(),
  50 fDoublets(),
  51 fIndex(),
  52 fNClusters(),
  53 fNDoublets(),
  54 fPhiCut(0.05),
  55 fRPhiCut(0.03),
  56 fZCut(0.01)
  57 #ifdef __DEBUG__
  58 ,fCv(0x0)
  59 ,fMk(0x0)
  60 ,fLn(0x0)
  61 ,fTx(0x0)
  62 #endif
  63 {
  64   //--------------------------------------------------------------------
  65   // This default constructor needs to be provided
  66   //--------------------------------------------------------------------
  67   for(Int_t i=0;i<7;++i) {
  68     fClusters[i].reserve(5000);
  69   }
  70   if (rec) Init(rec);
  71 }
  72
  73 //________________________________________________________________________________
  74 AliITSUTrackerSA::AliITSUTrackerSA(const AliITSUTrackerSA &t):
  75 AliTracker(t),
  76 fReconstructor(t.fReconstructor),
  77 fITS(t.fITS),
  78 fMatLUT(t.fMatLUT),
  79 fUseMatLUT(t.fUseMatLUT),
  80 fCurrMass(t.fCurrMass),
  81 //
  82 fClusters(),
  83 fClustersTC(),
  84 fIndex(),
  85 fNClusters(),
  86 fNDoublets(),
  87 fPhiCut(),
  88 fRPhiCut(),
  89 fZCut()
  90 #ifdef __DEBUG__
  91 ,fCv(0x0)
  92 ,fMk(0x0)
  93 ,fLn(0x0)
  94 ,fTx(0x0)
  95 #endif
  96 {
  97   //--------------------------------------------------------------------
  98   // The copy constructor is protected
  99   //--------------------------------------------------------------------
 100 }
 101
 102 //________________________________________________________________________________
 103 AliITSUTrackerSA::~AliITSUTrackerSA()
 104 {
 105   // d-tor
 106   delete fMatLUT;
 107 }
 108
 109
 110 //_________________________________________________________________________
 111 void AliITSUTrackerSA::Init(AliITSUReconstructor* rec)
 112 {
 113   // init with external reconstructor
 114   //
 115   fITS = rec->GetITSInterface();
 116   //
 117   // create material lookup table
 118   const int kNTest = 1000;
 119   const double kStepsPerCM=5;
 120   fMatLUT  = new AliITSUMatLUT(fITS->GetRMin(),fITS->GetRMax(),Nint(kStepsPerCM*(fITS->GetRMax()-fITS->GetRMin())));
 121   double zmn = 1e6;
 122   for (int ilr=fITS->GetNLayers();ilr--;) {
 123     AliITSURecoLayer* lr = fITS->GetLayer(ilr);
 124     if (zmn>Abs(lr->GetZMin())) zmn = Abs(lr->GetZMin());
 125     if (zmn>Abs(lr->GetZMax())) zmn = Abs(lr->GetZMax());
 126   }
 127   fMatLUT->FillData(kNTest,-zmn,zmn);
 128   //
 129 }
 130
 131 //________________________________________________________________________________
 132 Int_t AliITSUTrackerSA::Clusters2Tracks(AliESDEvent *event) {
 133   //--------------------------------------------------------------------
 134   // This is the main tracking function
 135   // The clusters must already be loaded
 136   //--------------------------------------------------------------------
 137
 138   // Possibly, create the track "seeds" (combinatorial)
 139
 140   // Possibly, increment the seeds with additional clusters (Kalman)
 141
 142   // Possibly, (re)fit the found tracks
 143
 144   // Three iterations:
 145   // - High momentum first;
 146   // - Low momentum with vertex constraint;
 147   // - Everything else.
 148
 149   MakeDoublets();       // To be checked
 150   //MakeTriplets();       // Are triplets really necessary? MFT does not use them.
 151   CASelection(event);
 152
 153   return 0;
 154 }
 155
 156 //________________________________________________________________________________
 157 Int_t AliITSUTrackerSA::PropagateBack(AliESDEvent * /*event*/) {
 158   //--------------------------------------------------------------------
 159   // Here, we implement the Kalman smoother ?
 160   // The clusters must be already loaded
 161   //--------------------------------------------------------------------
 162
 163   return 0;
 164 }
 165
 166 //________________________________________________________________________________
 167 Int_t AliITSUTrackerSA::RefitInward(AliESDEvent * /*event*/) {
 168   //--------------------------------------------------------------------
 169   // Some final refit, after the outliers get removed by the smoother ?
 170   // The clusters must be loaded
 171   //--------------------------------------------------------------------
 172
 173   return 0;
 174 }
 175
 176 //________________________________________________________________________________
 177 Int_t AliITSUTrackerSA::LoadClusters(TTree *cluTree) {
 178   //--------------------------------------------------------------------
 179   // This function reads the ITSU clusters from the tree,
 180   // sort them, distribute over the internal tracker arrays, etc
 181   //--------------------------------------------------------------------
 182   fITS->LoadClusters(cluTree);
 183   fITS->ProcessClusters();
 184
 185   //
 186   for(int iL=0; iL<7; ++iL) {
 187     fClustersTC[iL]=*fITS->GetLayerActive(iL)->GetClustersAddress();
 188     TClonesArray *clCont=fClustersTC[iL];
 189     fNClusters[iL]=clCont->GetEntriesFast();
 190     Float_t phi[fNClusters[iL]];
 191     fIndex[iL] = new Int_t[fNClusters[iL]];
 192     AliITSURecoLayer* lr = fITS->GetLayerActive(iL) ; // assign the layer which the cluster belongs to
 193     for(int iC=0;iC<fNClusters[iL];++iC) {
 194       const AliITSUClusterPix *cl = (AliITSUClusterPix*)clCont->At(iC);
 195       float pos[3];
 196       cl->GetGlobalXYZ(pos);
 197       phi[iC] = pos[0]==0.f ? TMath::PiOver2() : TMath::ATan2(pos[1]-GetY(),pos[0]-GetX());
 198       AliITSURecoSens* sens = lr->GetSensorFromID(cl->GetVolumeId());
 199       //double x = sens->GetXTF() + clus->GetX();
 200       float angle= sens->GetPhiTF();
 201       #ifdef __DEBUG__
 202       int label=cl->GetLabel(0);
 203       //cout << "Guarda te che label" << label << endl;
 204       fClusters[iL].push_back(itsCluster(pos[0],pos[1],pos[2],cl->GetSigmaY2(),cl->GetSigmaZ2(),cl->GetSigmaYZ(),phi[iC],angle,label));
 205       #else
 206       fClusters[iL].push_back(itsCluster(pos[0],pos[1],pos[2],cl->GetSigmaY2(),cl->GetSigmaZ2(),cl->GetSigmaYZ(),phi[iC],angle));
 207       #endif
 208     }
 209     TMath::Sort(fNClusters[iL],phi,fIndex[iL],kFALSE);
 210   }
 211   #ifdef __DEBUG__
 212   //PrintInfo("clusters");
 213   DrawEvent("clusters");
 214   #endif
 215   return 0;
 216 }
 217
 218 //________________________________________________________________________________
 219 void AliITSUTrackerSA::UnloadClusters() {
 220   //--------------------------------------------------------------------
 221   // This function unloads ITSU clusters from the RAM
 222   //--------------------------------------------------------------------
 223   for(int i=0;i<7;++i) {
 224     fClusters[i].clear();
 225     fNClusters[i]=0;
 226     delete fIndex[i];
 227   }
 228   for(int i=0;i<6;++i) fDoublets[i].clear();
 229 }
 230
 231 //________________________________________________________________________________
 232 AliCluster *AliITSUTrackerSA::GetCluster(Int_t /*index*/) const {
 233   //--------------------------------------------------------------------
 234   //       Return pointer to a given cluster
 235   //--------------------------------------------------------------------
 236   return 0;  // replace with an actual pointer
 237 }
 238
 239 //________________________________________________________________________________
 240 void AliITSUTrackerSA::CASelection(AliESDEvent *event) {
 241   // Here it's implemented the Cellular Automaton routine
 242   // Firstly the level of each doublet is set according to the level of
 243   // the neighbour doublets.
 244   // Doublet are considered to be neighbour if they share one point and the
 245   // phi and theta direction difference of the two is below a cut value.
 246
 247   cout << "Begin of the CA selection" << endl;
 248   for( int iL = 1; iL < 6; ++iL ) {
 249
 250     const itsCluster* clusters1 = &fClusters[iL-1][0];
 251     const itsCluster* clusters2 = &fClusters[iL][0];
 252     const itsCluster* clusters3 = &fClusters[iL+1][0];
 253
 254     const nPlets* doublets1 = &fDoublets[iL-1][0];
 255     nPlets* doublets2 = &fDoublets[iL][0];
 256
 257     for ( int iD2 = 0; iD2 < fNDoublets[iL]; ++iD2 ) {
 258       for ( int iD1 = 0; iD1 < fNDoublets[iL-1]; ++iD1 ) {
 259         const int id1 = doublets1[iD1].id1;
 260         const int id2 = doublets2[iD2].id0;
 261         if ( id1 == id2 ) {
 262           if ( doublets2[iD2].level <= ( doublets1[iD1].level + 1 ) ) {
 263             const int id3 = doublets2[iD2].id1;
 264             const float r3 = Sqrt( clusters3[id3].x * clusters3[id3].x + clusters3[id3].y * clusters3[id3].y );
 265             const float r2 = Sqrt( clusters2[id2].x * clusters2[id2].x + clusters2[id2].y * clusters2[id2].y );
 266             const float extrZ3 = doublets1[iD1].tanLambda * ( r3 - r2 ) + clusters2[id2].z ;
 267             //cout << extrZ3 << " " << clusters3[id3].z << " " << Abs ( extrZ3 - clusters3[id3].z ) << endl;
 268             if ( Abs ( extrZ3 - clusters3[id3].z ) < fZCut ) {
 269               //cout << "OK Z doublets: "<< iL-1 << "," << iD1 << "\t" << iL << "," <<iD2 << endl;
 270               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);
 271               //cout << det << endl;
 272               if ( Abs(det) <= 1e-12 ) {
 273                 // linear extrapolation to next layer
 274                 const float dsq = ( doublets1[iD1].tanPhi * (clusters3[id3].x + clusters2[id2].x) + clusters3[id3].y - clusters2[id2].y ) *
 275                 ( doublets1[iD1].tanPhi * (clusters3[id3].x + clusters2[id2].x) + clusters3[id3].y - clusters2[id2].y ) / (1 + doublets1[iD1].tanPhi * doublets1[iD1].tanPhi );
 276                 if ( dsq < fRPhiCut*fRPhiCut )  {
 277                   doublets2[iD2].level = doublets1[iD1].level+1;
 278                   doublets2[iD2].neighbours.push_back(iD1);
 279                 }
 280               } else {
 281                 const float r1sq = clusters1[id1].x * clusters1[id1].x + clusters1[id1].y * clusters1[id1].y ;
 282                 const float rvsq = GetX() * GetX() + GetY() * GetY();
 283                 const float deta = (rvsq - r1sq) * (clusters2[id2].y - GetY()) - (rvsq - r2*r2) * (clusters1[id1].y - GetY());
 284                 const float detb = - (rvsq - r1sq) * (clusters2[id2].x - GetX()) + (rvsq - r2*r2) * (clusters1[id1].x - GetX()) ;
 285                 const float a = deta/det ;
 286                 const float b = detb/det ;
 287                 const float c = -rvsq - a * GetX() - b * GetY();
 288                 const float rc = Sqrt( a*a/4.f + b*b/4.f - c );
 289                 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) );
 290                 //cout << d << " " << rc << " " << d - rc << endl;
 291                 if ( Abs( d - rc ) < fRPhiCut ) {
 292                   doublets2[iD2].level = doublets1[iD1].level+1;
 293                   doublets2[iD2].neighbours.push_back(iD1);
 294                 }
 295               }
 296             }
 297           }
 298         }
 299       }
 300     }
 301   }
 302   //#ifdef __DEBUG__
 303   //PrintInfo("doublets");
 304   //DrawEvent("doublets+level");
 305   //return;
 306   //#endif
 307   // Hic sunt leones: the following code could be optimised to be iterative. But now I don't have time.
 308   vector<trackC> tracks;
 309   for ( int level = 6; level >= 2 ; --level ) {
 310     cout << "level " << level << endl;
 311     vector<Road> roads;
 312     roads.clear();
 313     roads.reserve(100);
 314     for ( int doubl = 5; doubl >= level-1; --doubl ) {
 315       for ( int iD = 0; iD < fNDoublets[doubl]; ++iD ) {
 316         if ( fDoublets[doubl][iD].level == level ) {
 317           roads.push_back(Road());
 318           roads.back().AddElement(doubl,iD);
 319           cout << "\nseed " << iD << "("<<fDoublets[doubl][iD].id1<<","<<fDoublets[doubl][iD].id0;
 320           cout <<") in array " << doubl << " with level " << level;
 321           for ( unsigned int iN = 0; iN < fDoublets[doubl][iD].neighbours.size(); ++iN ) {
 322             const int currD = doubl - 1 ;
 323             const int neigh = fDoublets[doubl][iD].neighbours[iN];
 324             //if ( level != fDoublets[currD][neigh].level + 1 ) continue;
 325             if ( iN > 0 ) roads.push_back(static_cast<Road>(roads.back()));
 326             CandidatesTreeTraversal(roads,neigh,currD);
 327             cout << endl;// << roads.back() << endl;
 328           }
 329           fDoublets[doubl][iD].level = -1; // mark as used
 330         }
 331         //for ( int j = 0; j < 2*AliITSUAux::kMaxLayers; ++j ) cout << roads.back().fPoints[j] << " ";
 332         //cout << endl;
 333       }
 334     }
 335
 336
 337     Double_t rDest = 0.;//fITS->GetRMax();
 338     DrawRoads(roads);
 339
 340     vector<trackC> candidates;
 341     for ( size_t iR = 0; iR < roads.size(); ++iR ) {
 342       if ( roads[iR].fNElements!=level ) {
 343         cout << "JUMP" << endl;
 344         continue;
 345       }
 346       candidates.push_back(trackC());
 347       candidates.back().fNPoints = level+1;
 348       for ( size_t j = 0; j < 6; ++j ) {
 349         if ( roads[iR].fElements[j] == -1 ) continue;
 350         candidates.back().fPoints[j<<0x1] = fDoublets[j][roads[iR].fElements[j]].id0;
 351         candidates.back().fPoints[(j+1)<<0x1] = fDoublets[j][roads[iR].fElements[j]].id1;
 352         //cout << (j<<0x1) << " " << fDoublets[j][roads[iR].fElements[j]].id0 << "\t" << ((j+1)<<0x1) << " " << fDoublets[j][roads[iR].fElements[j]].id1 << endl;
 353       }
 354       //cout << endl;
 355       cout << "Candidate " << candidates.size() << ", number of points: " << level+1 << endl;
 356       //cout << candidates.back() << endl;
 357       InitTrackParams(candidates.back());
 358       candidates.back().fChi2 = RefitTrack( (AliExternalTrackParam*)&candidates.back(), candidates.back().fPoints, rDest ,-1);
 359       //cout << "Fit cnd: " << cand << " " << candidates[cand] << endl;
 360     }
 361
 362     int index[candidates.size()];
 363     for ( size_t i = 0; i < candidates.size(); ++i ) index[i]=i;
 364     CompDesc comp(&candidates);
 365     sort(index,index+candidates.size(),comp);
 366
 367     for ( size_t cand = 0; cand < candidates.size(); ++cand ) {
 368       const int ii = index[cand];
 369       //#ifdef __DEBUG__
 370       //cout << ii << " " << candidates[ii] << endl;
 371       //#endif
 372       if ( candidates[ii].fChi2 < 0. ) break;
 373       bool goodTrack = true;
 374       for ( unsigned int point = 0; point < 14; ++point ) { //-> here it's necessary to use a temporary array for the used clusters.
 375         if ( candidates[ii].fPoints[point] != -1 ) {
 376           if( !(fClusters[ point/2 ][ candidates[ii].fPoints[point] ].isUsed ) ) {
 377             fClusters[ point/2 ][ candidates[ii].fPoints[point] ].isUsed = true;
 378           } else {
 379             goodTrack = false;
 380           }
 381         }
 382       }
 383       //cout << endl;
 384       if ( goodTrack ) {
 385         tracks.push_back(candidates[ii]);
 386       }
 387     }
 388     //cout << "End of level " << level << endl;*/
 389   }
 390
 391   bool joined[tracks.size()];
 392   MergeTracks(tracks,joined);
 393   for ( unsigned int ii = 0; ii < tracks.size(); ++ii ) {
 394     if ( ! joined[ii] ) {
 395       AliESDtrack outTrack;
 396       outTrack.SetOuterParam((AliExternalTrackParam*)&tracks[ii],AliESDtrack::kITSpureSA);
 397       event->AddTrack(&outTrack);
 398     }
 399   }
 400 }
 401
 402 //________________________________________________________________________________
 403 void AliITSUTrackerSA::MakeDoublets() {
 404   // Make associations between two points on adjacent layers within an azimuthal window.
 405   // Under consideration:
 406   // - track parameter estimation using the primary vertex position
 407   // To do:
 408   // - last iteration
 409
 410   //cout << "Vertex of used by the tracker: " << GetX() << " " << GetY() << " " << GetZ() << endl;
 411
 412   for( int iL = 0 ; iL < 6 ; ++iL ) {
 413     fNDoublets[iL] = 0;
 414     const itsCluster* clusters1 = &fClusters[iL][0];
 415     const itsCluster* clusters2 = &fClusters[iL+1][0];
 416
 417     // 0 - 2Pi junction treatment (part I)
 418     for ( int iCC1 = 0 ; iCC1 < fNClusters[iL] ; ++iCC1 ) {
 419       bool flag = true;
 420       const int iC1 = fIndex[iL][iCC1];
 421       for ( int iCC2 = fNClusters[iL+1]-1; iCC2 >= 0 ; --iCC2 ) {
 422         const int iC2 = fIndex[iL+1][iCC2];
 423         if( (TMath::TwoPi() - (clusters2[iC2].phi-clusters1[iC1].phi) ) < fPhiCut ) {
 424           #ifdef __DEBUG__
 425           fDoublets[iL].push_back(nPlets(iC1,iC2,clusters1[iC1].pid,clusters2[iC2].pid));
 426           #else
 427           fDoublets[iL].push_back(nPlets(iC1,iC2));
 428           #endif
 429           fDoublets[iL][fNDoublets[iL]].tanPhi = (clusters1[iC1].y-clusters2[iC2].y)/(clusters1[iC1].x-clusters2[iC2].x);
 430           float r1  = Sqrt(clusters1[iC1].x * clusters1[iC1].x + clusters1[iC1].y * clusters1[iC1].y);
 431           //cout << clusters2[iC2].x * clusters2[iC2].x + clusters2[iC2].y * clusters2[iC2].y << flush << endl;
 432           float r2  = Sqrt(clusters2[iC2].x * clusters2[iC2].x + clusters2[iC2].y * clusters2[iC2].y);
 433           fDoublets[iL][fNDoublets[iL]].tanLambda = (clusters1[iC1].z-clusters2[iC2].z)/(r1-r2);
 434           ++fNDoublets[iL];
 435           flag = false;
 436         } else break;
 437
 438       }
 439       if (flag) break;
 440     }
 441
 442
 443     // "Central" points
 444     for ( int iCC1 = 0 ; iCC1 < fNClusters[iL] ; ++iCC1 ) {
 445       const int iC1 = fIndex[iL][iCC1];
 446       for ( int iCC2 = 0; iCC2 < fNClusters[iL+1] ; ++iCC2 ) {
 447         const int iC2 = fIndex[iL+1][iCC2];
 448         if( Abs( clusters1[iC1].phi - clusters2[iC2].phi ) < fPhiCut ) {
 449           #ifdef __DEBUG__
 450           fDoublets[iL].push_back(nPlets(iC1,iC2,clusters1[iC1].pid,clusters2[iC2].pid));
 451           #else
 452           fDoublets[iL].push_back(nPlets(iC1,iC2));
 453           #endif
 454           fDoublets[iL][fNDoublets[iL]].tanPhi = (clusters1[iC1].y-clusters2[iC2].y)/(clusters1[iC1].x-clusters2[iC2].x);
 455           float r1  = Sqrt(clusters1[iC1].x * clusters1[iC1].x + clusters1[iC1].y * clusters1[iC1].y);
 456           float r2  = Sqrt(clusters2[iC2].x * clusters2[iC2].x + clusters2[iC2].y * clusters2[iC2].y);
 457           fDoublets[iL][fNDoublets[iL]].tanLambda = (clusters1[iC1].z-clusters2[iC2].z)/(r1-r2);
 458           ++fNDoublets[iL];
 459         } else if( clusters2[iC2].phi - clusters1[iC1].phi > fPhiCut ) break;
 460
 461       }
 462
 463     }
 464
 465     // 0 - 2Pi junction treatment (part II)
 466     for ( int iCC1 = fNClusters[iL]-1; iCC1 > -1 ; --iCC1 ) {
 467       bool flag = true;
 468       const int iC1 = fIndex[iL][iCC1];
 469       for ( int iCC2 = 0; iCC2 < fNClusters[iL+1] ; ++iCC2 ) {
 470         const int iC2 = fIndex[iL+1][iCC2];
 471         if( (TMath::TwoPi() - (clusters1[iC1].phi-clusters2[iC2].phi) ) < fPhiCut ) {
 472           #ifdef __DEBUG__
 473           fDoublets[iL].push_back(nPlets(iC1,iC2,clusters1[iC1].pid,clusters2[iC2].pid));
 474           #else
 475           fDoublets[iL].push_back(nPlets(iC1,iC2));
 476           #endif
 477           fDoublets[iL][fNDoublets[iL]].tanPhi = (clusters1[iC1].y-clusters2[iC2].y)/(clusters1[iC1].x-clusters2[iC2].x);
 478           float r1  = Sqrt(clusters1[iC1].x * clusters1[iC1].x + clusters1[iC1].y * clusters1[iC1].y);
 479           float r2  = Sqrt(clusters2[iC2].x * clusters2[iC2].x + clusters2[iC2].y * clusters2[iC2].y);
 480           fDoublets[iL][fNDoublets[iL]].tanLambda = (clusters1[iC1].z-clusters2[iC2].z)/(r1-r2);
 481           ++fNDoublets[iL];
 482           flag = false;
 483         } else break;
 484
 485       }
 486
 487       if (flag) break;
 488     }
 489
 490   }
 491   // #ifdef __DEBUG__
 492   // PrintInfo("doublets");
 493   // #endif
 494 }
 495 //______________________________________________________________________________
 496 Bool_t AliITSUTrackerSA::InitTrackParams(trackC &track)
 497 {
 498   // Set the initial guess on track kinematics for propagation.
 499   // Assume at least 3 points available
 500   int lrOcc[AliITSUAux::kMaxLayers], nCl=0;
 501   //
 502   // we will need endpoints and middle layer
 503   for (int i=fITS->GetNLayersActive();i--;) {
 504     if (track.fPoints[i<<0x1]>-1) {
 505       lrOcc[nCl++] = i;
 506       track.fInnermostLayer = i;
 507     }
 508   }
 509   track.fNPoints = nCl;
 510   track.fOutermostLayer = track.fInnermostLayer + nCl - 1;
 511   if (nCl<3) {
 512     AliError(Form("Cannot estimate momentum of tracks with %d clusters",nCl));
 513     //cout << track << endl;
 514     return kFALSE;
 515   }
 516   //
 517   int lr0   = lrOcc[0];
 518   int lr1   = lrOcc[nCl/2];
 519   int lr2   = lrOcc[nCl-1];
 520   //
 521   //cout << (lr0<<0x1) << " " <<  (lr1<<0x1)<< " " <<  (lr2<<0x1) << endl;
 522   const itsCluster& cl0 = fClusters[lr0][ track.fPoints[lr0<<0x1] ];
 523   const itsCluster& cl1 = fClusters[lr1][ track.fPoints[lr1<<0x1] ];
 524   const itsCluster& cl2 = fClusters[lr2][ track.fPoints[lr2<<0x1] ];
 525   double cv = Curvature(cl0.x,cl0.y, cl1.x,cl1.y, cl2.x,cl2.y);
 526   double tgl = (cl2.z-cl0.z)/TMath::Sqrt((cl2.x-cl0.x)*(cl2.x-cl0.x)+(cl2.y-cl0.y)*(cl2.y-cl0.y));
 527   //  double phi = TMath::ATan2((cl2.y-cl1.y),(cl2.x-cl1.x));
 528   //
 529   AliITSUClusterPix* clus = (AliITSUClusterPix*)fClustersTC[ lr0 ]->At( track.fPoints[lr0<<0x1] );
 530   AliITSURecoLayer* lr = fITS->GetLayerActive(lr0);
 531   AliITSURecoSens* sens = lr->GetSensorFromID(clus->GetVolumeId());
 532   double x = sens->GetXTF() + clus->GetX();
 533   double alp = sens->GetPhiTF();
 534   //  printf("Alp: %f phi: %f\n",alp,phi);
 535   double par[5] = {clus->GetY(),clus->GetZ(),0,tgl,cv};
 536   double cov[15] = {
 537     5*5,
 538     0, 5*5,
 539     0, 0, 0.7*0.7,
 540     0,0,0,0.7*0.7,
 541     0,0,0,0,10
 542   };
 543   track.Set(x,alp,par,cov);
 544   //cout << track ;
 545   return kTRUE;
 546 }
 547
 548 //______________________________________________________________________________
 549 void AliITSUTrackerSA::CandidatesTreeTraversal(vector<Road> &candidates, const int &iD, const int &doubl) {
 550
 551   if ( doubl < 0     ) {
 552     #ifdef __DEBUG__
 553     //cout << "ERROR IN CandidatesTreeTraversal" << endl;
 554     //cout << endl;
 555     #endif
 556     return;
 557   }
 558   candidates.back().AddElement(doubl,iD);
 559   cout << "\n\tTraversing the tree through " << iD << " in array " << doubl << " at level " << fDoublets[doubl][iD].level;
 560   //cout << doubl*2 << " ";
 561   for ( unsigned int iN = 0; iN < fDoublets[doubl][iD].neighbours.size(); ++iN ) {
 562     const int currD = doubl - 1 ;
 563     const int neigh = fDoublets[doubl][iD].neighbours[iN];
 564     //const int level = fDoublets[doubl][iD].level;
 565     //if ( level != fDoublets[currD][neigh].level + 1 ) continue;
 566     if ( iN > 0 ) {
 567       cout << " -> branching here!";
 568       candidates.push_back(static_cast<Road>(candidates.back()));
 569     }
 570     CandidatesTreeTraversal(candidates,neigh,currD);
 571   }
 572
 573   fDoublets[doubl][iD].level = -1;
 574
 575 }
 576
 577 //______________________________________________________________________________
 578 void AliITSUTrackerSA::MergeTracks(vector<trackC> &tracks, bool flag[]) {
 579   cout << "Merging tracks" << endl;
 580   for ( unsigned int iT = 0; iT < tracks.size(); ++iT ) {
 581     flag[iT]=false;
 582   }
 583
 584   for ( unsigned int iT1 = 0; iT1 < tracks.size(); ++iT1 ) {
 585     if ( tracks[iT1].fNPoints > 4 || flag[iT1] ) continue;
 586     const int inL1 = tracks[iT1].fInnermostLayer;
 587     const int outL1 = tracks[iT1].fOutermostLayer;
 588     for ( unsigned int iT2 = iT1+1; iT2 < tracks.size(); ++iT2 ) {
 589       const int inL2 = tracks[iT2].fInnermostLayer;
 590       const int outL2 = tracks[iT2].fOutermostLayer;
 591       printf("%d: In/Out %d/%d\t%d: In/Out %d/%d",iT1,inL1,outL1,iT2,inL2,outL2);
 592       if ( outL1 < inL2 ) {
 593         cout << " <- check1" << flush;
 594         if ( TransportToLayer((AliExternalTrackParam*)&tracks[iT2],inL2,outL1) ) {
 595           cout << " <- transport";
 596           if ( tracks[iT2].Rotate( fClusters[outL1][ tracks[iT1].fPoints[ outL1<<0x1 ] ].phiM ) ) {
 597             const double chi2 = tracks[iT2].GetPredictedChi2((AliExternalTrackParam*)&tracks[iT1]);
 598             printf(" Merging candidates (%i,%i), chi2 = %f ",iT1,iT2,chi2);
 599             if ( chi2 < fgkChi2Cut ) {
 600               flag[iT2] = true;
 601               for ( int np = tracks[iT2].fInnermostLayer; np < tracks[iT2].fOutermostLayer; ++np ) {
 602                 tracks[iT1].fPoints[np<<0x1] = tracks[iT2].fPoints[np<<0x1];
 603               }
 604               InitTrackParams(tracks[iT1]);
 605               RefitTrack((AliExternalTrackParam*)&tracks[iT1],tracks[iT1].fPoints,0,-1);
 606             }
 607           } else { cout << " <- Failed rotation!!"; }
 608         } else { cout << " <- Failed transport!!"; }
 609       } else if ( inL1 > outL2 ) {
 610         cout << " <- check2" << flush;
 611         if ( TransportToLayer((AliExternalTrackParam*)&tracks[iT1],inL1,outL2) ) {
 612           cout << " <- transport";
 613           if ( tracks[iT1].Rotate( fClusters[outL2][ tracks[iT1].fPoints[ outL2<<0x1 ] ].phiM ) ) {
 614             const double chi2 = tracks[iT2].GetPredictedChi2((AliExternalTrackParam*)&tracks[iT1]);
 615             printf(" Merging candidates (%i,%i), chi2 = %f ",iT1,iT2,chi2);
 616             if ( chi2 < fgkChi2Cut ) {
 617               flag[iT2] = true;
 618               for ( int np = tracks[iT2].fInnermostLayer; np < tracks[iT2].fOutermostLayer; ++np ) {
 619                 tracks[iT1].fPoints[np<<0x1] = tracks[iT2].fPoints[np<<0x1];
 620               }
 621               InitTrackParams(tracks[iT1]);
 622               RefitTrack((AliExternalTrackParam*)&tracks[iT1],tracks[iT1].fPoints,0,-1);
 623             }
 624           } else { cout << " <- Failed rotation!!"; }
 625         } else { cout << " <- Failed transport!!"; }
 626       }
 627       cout << endl;
 628     }
 629   }
 630 }
 631
 632   //______________________________________________________________________________
 633 Double_t AliITSUTrackerSA::RefitTrack(AliExternalTrackParam* trc,
 634 Int_t clInfo[2*AliITSUAux::kMaxLayers],
 635 Double_t rDest, Int_t stopCond)
 636 {
 637   // refit track till radius rDest.
 638   // if stopCond<0 : propagate till last cluster then stop
 639   // if stopCond==0: propagate till last cluster then try to go till limiting rDest, don't mind if fail
 640   // if stopCond>0 : rDest must be reached
 641   //
 642   // The clList should provide the indices of clusters at corresponding layer (as stored in the layer
 643   // TClonesArray, with convention (allowing for up to 2 clusters per layer due to the overlaps):
 644   // if there is a cluster on given layer I, then it should be stored at clInfo[2*I-1]
 645   // if there is an additional cluster on this layer, it goes to clInfo[2*I]
 646   // -1 means no cluster
 647   //
 648   double rCurr = Sqrt(trc->GetX()*trc->GetX() + trc->GetY()*trc->GetY());
 649   int dir,lrStart,lrStop;
 650   //
 651   dir = rCurr<rDest ? 1 : -1;
 652   lrStart = fITS->FindFirstLayerID(rCurr,dir);
 653   lrStop  = fITS->FindLastLayerID(rDest,dir); // lr id before which we have to stop
 654   //cout << "Start/End layers and direction " << lrStart << " " << lrStop << " " <<  dir << endl;
 655   //
 656   if (lrStop<0 || lrStart<0) AliFatal(Form("Failed to find start(%d) or last(%d) layers. "
 657   "Track from %.3f to %.3f",lrStart,lrStop,rCurr,rDest));
 658   //
 659   int nCl = 0;
 660   cout << "pid : ";
 661   for (int i=2*fITS->GetNLayersActive();i--;) {
 662     if (clInfo[i]<0) continue;
 663     cout << fClusters[i/2][clInfo[i]].pid;
 664     ++nCl;
 665   }
 666   cout << endl;
 667   //
 668   //cout << "#ptr: " << nCl << endl;
 669   AliExternalTrackParam tmpTr(*trc);
 670   double chi2 = 0;
 671   int iclLr[2],nclLr;
 672   int nclFit = 0;
 673   //
 674   int lrStop1 = lrStop+dir;
 675   for (int ilr=lrStart;ilr!=lrStop1;ilr+=dir) {
 676     AliITSURecoLayer* lr = fITS->GetLayer(ilr);
 677     if ( dir*(rCurr-lr->GetR(dir))>0) {
 678       cout << ilr << " passed!" << endl;
 679       continue;
 680     } // this layer is already passed
 681     int ilrA2,ilrA = lr->GetActiveID();
 682     // passive layer or active w/o hits will be traversed on the way to next cluster
 683     if (!lr->IsActive() || clInfo[ilrA2=(ilrA<<1)]<0) {
 684       //cout << ilr << " is inactive or without cluster for current candidates" << endl;
 685       continue;
 686     }
 687     //cout << "OK layer " << ilr << endl;
 688     //
 689     // select the order in which possible 2 clusters (in case of the overlap) will be traversed and fitted
 690     nclLr=0;
 691     if (dir>0) { // clusters are stored in increasing radius order
 692       iclLr[nclLr++]=clInfo[ilrA2++];
 693       if (clInfo[ilrA2]>=0) iclLr[nclLr++]=clInfo[ilrA2];
 694     }
 695     else {
 696       if ( clInfo[ilrA2+1]>=0 ) iclLr[nclLr++]=clInfo[ilrA2+1];
 697       iclLr[nclLr++]=clInfo[ilrA2];
 698     }
 699     //
 700     Bool_t transportedToLayer = kFALSE;
 701     for (int icl=0;icl<nclLr;icl++) {
 702       AliITSUClusterPix* clus =  (AliITSUClusterPix*)lr->GetCluster(iclLr[icl]);
 703       AliITSURecoSens* sens = lr->GetSensorFromID(clus->GetVolumeId());
 704       if (!tmpTr.Rotate(sens->GetPhiTF())) { cout << "failed rotation" << endl; return -1; }
 705       //
 706       double xClus = sens->GetXTF()+clus->GetX();
 707       if (!transportedToLayer) {
 708         if (ilr!=lrStart && !TransportToLayerX(&tmpTr,lrStart,ilr,xClus)) {
 709           cout << "failed transport to the entrance" << endl; return -1; } // go to the entrance to the layer
 710           lrStart = ilr;
 711           transportedToLayer = kTRUE;
 712         }
 713         //
 714         if (!PropagateSeed(&tmpTr,xClus,fCurrMass)) { cout << "failed propagation of the seed X:" << xClus << endl; tmpTr.Print(); return -1; }
 715         //
 716         Double_t p[2]={clus->GetY(), clus->GetZ()};
 717         Double_t cov[3]={clus->GetSigmaY2(), clus->GetSigmaYZ(), clus->GetSigmaZ2()};
 718         double chi2cl = tmpTr.GetPredictedChi2(p,cov);
 719         chi2 += chi2cl;
 720         //
 721         if ( !tmpTr.Update(p,cov) ) {
 722           cout << "failed update of the covariance" << endl;
 723           return -1;
 724         }
 725
 726         if (++nclFit==nCl && stopCond<0) {
 727           *trc = tmpTr;
 728           printf("Fit chi2: %f for %d clusters\n",chi2,nclFit);
 729           return chi2; // it was requested to not propagate after last update
 730         }
 731       }
 732       //
 733     }
 734     // All clusters were succesfully fitted. Even if the track does not reach rDest, this is enough to validate it.
 735     // Still, try to go as close as possible to rDest.
 736     //
 737     //  printf("Fit chi2: %f for %d clusters\n",chi2,nclFit);
 738     //
 739     if (lrStart!=lrStop) {
 740       if (!TransportToLayer(&tmpTr,lrStart,lrStop)) return (stopCond>0) ? -chi2 : chi2; // rDest was obligatory
 741       if (!GoToExitFromLayer(&tmpTr,fITS->GetLayer(lrStop),dir)) return (stopCond>0) ? -chi2 : chi2; // rDest was obligatory
 742     }
 743     // go to the destination radius. Note that here we don't select direction to avoid precision problems
 744     if (!tmpTr.GetXatLabR(rDest,rDest,GetBz(),0) || !PropagateSeed(&tmpTr,rDest,fCurrMass, 100, kFALSE)) {
 745       return (stopCond>0) ? -chi2 : chi2; // rDest was obligatory
 746     }
 747     *trc = tmpTr;
 748
 749     return chi2;
 750   }
 751
 752   //______________________________________________________________________________
 753   Bool_t AliITSUTrackerSA::PropagateSeed(AliExternalTrackParam *seed, Double_t xToGo, Double_t mass, Double_t maxStep, Bool_t matCorr)
 754   {
 755     // propagate seed to given x applying material correction if requested
 756     const Double_t kEpsilon = 1e-5;
 757     Double_t xpos     = seed->GetX();
 758     Int_t dir         = (xpos<xToGo) ? 1:-1;
 759     Double_t xyz0[3],xyz1[3];
 760     //
 761     Bool_t updTime = dir>0 && seed->IsStartedTimeIntegral();
 762     if (matCorr || updTime) seed->GetXYZ(xyz1);   //starting global position
 763     while ( (xToGo-xpos)*dir > kEpsilon){
 764       Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
 765       Double_t x    = xpos+step;
 766       Double_t bz=GetBz();   // getting the local Bz
 767       if (!seed->PropagateTo(x,bz)) {
 768         //cout << " failed PropagateTo " << endl;
 769         return kFALSE;
 770       }
 771       double ds = 0;
 772       if (matCorr || updTime) {
 773         xyz0[0]=xyz1[0]; // global pos at the beginning of step
 774         xyz0[1]=xyz1[1];
 775         xyz0[2]=xyz1[2];
 776         seed->GetXYZ(xyz1);    //  // global pos at the end of step
 777         //
 778         if (matCorr) {
 779           Double_t xrho,xx0;
 780           ds = GetMaterialBudget(xyz0,xyz1,xx0,xrho);
 781           if (dir>0) xrho = -xrho; // outward should be negative
 782           if (!seed->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
 783         }
 784         else { // matCorr is not requested but time integral is
 785           double d0 = xyz1[0]-xyz0[0];
 786           double d1 = xyz1[1]-xyz0[1];
 787           double d2 = xyz1[2]-xyz0[2];
 788           ds = TMath::Sqrt(d0*d0+d1*d1+d2*d2);
 789         }
 790       }
 791       if (updTime) seed->AddTimeStep(ds);
 792       //
 793       xpos = seed->GetX();
 794     }
 795     return kTRUE;
 796   }
 797
 798   //_________________________________________________________________________
 799   Bool_t AliITSUTrackerSA::TransportToLayer(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo, Double_t rLim)
 800   {
 801     // transport track from layerFrom to the entrance of layerTo or to rLim (if>0), wathever is closer
 802     //
 803     if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
 804     //
 805     int dir = lTo > lFrom ? 1:-1;
 806     //printf("From %d, to %d, direction %d ",lFrom,lTo,dir);
 807     AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
 808     Bool_t checkFirst = kTRUE;
 809     Bool_t limReached = kFALSE;
 810     while(lFrom!=lTo) {
 811       if (lrFr) {
 812         if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) {
 813           cout << " Failed GoToExitFromLayer ";
 814           return kFALSE; // go till the end of current layer
 815         }
 816         checkFirst = kFALSE;
 817       }
 818       AliITSURecoLayer* lrTo =  fITS->GetLayer( (lFrom+=dir) );
 819       if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
 820       //
 821       // go the entrance of the layer, assuming no materials in between
 822       double xToGo = lrTo->GetR(-dir);
 823       if (rLim>0) {
 824         if (dir>0) {
 825           if (rLim<xToGo) {xToGo = rLim; limReached = kTRUE;}
 826         }
 827         else {
 828           if (rLim>xToGo) {xToGo = rLim; limReached = kTRUE;}
 829         }
 830       }
 831       //    double xts = xToGo;
 832       if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) {
 833         //      printf("FailHere1: %f %f %d\n",xts,xToGo,dir);
 834         //      seed->Print("etp");
 835         cout << " Failed GetXatLabR " << endl;
 836         return kFALSE;
 837       }
 838       if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) {
 839         //printf("FailHere2: %f %f %d\n",xts,xToGo,dir);
 840         //seed->Print("etp");
 841         cout << " Failed PropagateSeed " << endl;
 842         return kFALSE;
 843       }
 844       lrFr = lrTo;
 845       if (limReached) break;
 846     }
 847     return kTRUE;
 848     //
 849   }
 850
 851   //_________________________________________________________________________
 852   Bool_t AliITSUTrackerSA::TransportToLayerX(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo, Double_t xStop)
 853   {
 854     // transport track from layerFrom to the entrance of layerTo but do not pass control parameter X
 855     //
 856     if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
 857     //
 858     int dir = lTo > lFrom ? 1:-1;
 859     AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
 860     Bool_t checkFirst = kTRUE;
 861     while(lFrom!=lTo) {
 862       if (lrFr) {
 863         if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) return kFALSE; // go till the end of current layer
 864         checkFirst = kFALSE;
 865       }
 866       AliITSURecoLayer* lrTo =  fITS->GetLayer( (lFrom+=dir) );
 867       if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
 868       //
 869       // go the entrance of the layer, assuming no materials in between
 870       double xToGo = lrTo->GetR(-dir); // R of the entrance to layer
 871       //
 872       //    double xts = xToGo;
 873       if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) {
 874         //      printf("FailHere1: %f %f %d\n",xts,xToGo,dir);
 875         //      seed->Print("etp");
 876         return kFALSE;
 877       }
 878       if ( (dir>0&&xToGo>xStop) || (dir<0&&xToGo<xStop) ) xToGo = xStop;
 879       //
 880       #ifdef _ITSU_DEBUG_
 881       AliDebug(2,Form("go in dir=%d to R=%.4f(X:%.4f)",dir,lrTo->GetR(-dir), xToGo));
 882       #endif
 883       if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) {
 884         //printf("FailHere2: %f %f %d\n",xts,xToGo,dir);
 885         //seed->Print("etp");
 886         return kFALSE;
 887       }
 888       lrFr = lrTo;
 889     }
 890     return kTRUE;
 891     //
 892   }
 893
 894   //_________________________________________________________________________
 895   Bool_t AliITSUTrackerSA::GoToExitFromLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
 896   {
 897     // go to the exit from lr in direction dir, applying material corrections in steps specific for this layer
 898     // If check is requested, do this only provided the track has not exited the layer already
 899     double xToGo = lr->GetR(dir);
 900     if (check) { // do we need to track till the surface of the current layer ?
 901       double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
 902       if      (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // on the surface or outside of the layer
 903       else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // on the surface or outside of the layer
 904     }
 905     if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
 906     // go via layer to its boundary, applying material correction.
 907     if (!PropagateSeed(seed,xToGo,fCurrMass, lr->GetMaxStep())) return kFALSE;
 908     //
 909     return kTRUE;
 910     //
 911   }
 912
 913   //_________________________________________________________________________
 914   Bool_t AliITSUTrackerSA::GoToEntranceToLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
 915   {
 916     // go to the entrance of lr in direction dir, w/o applying material corrections.
 917     // If check is requested, do this only provided the track did not reach the layer already
 918     double xToGo = lr->GetR(-dir);
 919     if (check) { // do we need to track till the surface of the current layer ?
 920       double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
 921       if      (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // already passed
 922       else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // already passed
 923     }
 924     if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
 925     // go via layer to its boundary, applying material correction.
 926     if (!PropagateSeed(seed,xToGo,fCurrMass, 100, kFALSE)) return kFALSE;
 927     return kTRUE;
 928     //
 929   }
 930
 931   //____________________________________________________
 932   Double_t AliITSUTrackerSA::GetMaterialBudget(const double* pnt0,const double* pnt1, double& x2x0, double& rhol) const
 933   {
 934     double par[7];
 935     if (fUseMatLUT && fMatLUT) {
 936       double d = fMatLUT->GetMatBudget(pnt0,pnt1,par);
 937       x2x0 = par[AliITSUMatLUT::kParX2X0];
 938       rhol = par[AliITSUMatLUT::kParRhoL];
 939       return d;
 940     }
 941     else {
 942       MeanMaterialBudget(pnt0,pnt1,par);
 943       x2x0 = par[1];
 944       rhol = par[0]*par[4];
 945       return par[4];
 946     }
 947   }
 948
 949   //____________________________________________________________________
 950   Double_t AliITSUTrackerSA::Curvature(Double_t x1,Double_t y1,Double_t
 951   x2,Double_t y2,Double_t x3,Double_t y3)
 952   {
 953
 954     //calculates the curvature of track
 955     Double_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
 956     if(den==0) return 0;
 957     Double_t a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
 958     Double_t b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
 959     Double_t c = -x1*x1-y1*y1-a*x1-b*y1;
 960     Double_t xc=-a/2.;
 961
 962     if((a*a+b*b-4*c)<0) return 0;
 963     Double_t rad = TMath::Sqrt(a*a+b*b-4*c)/2.;
 964     if(rad==0) return 0;
 965
 966     if((x1>0 && y1>0 && x1<xc)) rad*=-1;
 967     if((x1<0 && y1>0 && x1<xc)) rad*=-1;
 968     //  if((x1<0 && y1<0 && x1<xc)) rad*=-1;
 969     // if((x1>0 && y1<0 && x1<xc)) rad*=-1;
 970
 971     return 1/rad;
 972
 973   }
 974
 975   #ifdef __DEBUG__
 976   //____________________________________________________
 977   void AliITSUTrackerSA::PrintInfo(TString what) {
 978     //
 979     if( what.Contains("clusters") ) {
 980       cout << "Dumping clusters info" << endl;
 981       for ( int i = 0; i < 7; ++i ) {
 982         cout << "**** Layer " << i << " ****" << endl;
 983         for ( int c = 0; c < fNClusters[i]; ++c ) {
 984           cout << "*** Cluster " << c << " ***" <<endl;
 985           cout << fClusters[i][fIndex[i][c]] << endl;
 986         }
 987       }
 988     }
 989     //
 990     if( what.Contains("doublets") ) {
 991       cout << "Dumping doublets info" << endl;
 992       for ( int i = 0; i < 6; ++i ) {
 993         cout << "**** Doublets array " << i << " ****" << endl;
 994         for ( int c = 0; c < fNDoublets[i]; ++c ) {
 995           cout << "*** Doublet " << c << " ***" <<endl;
 996           cout << fDoublets[i][c] << endl;
 997         }
 998       }
 999     }
1000   }
1001
1002   //____________________________________________________
1003   void AliITSUTrackerSA::DrawEvent(TString what) {
1004     //
1005
1006     const int size = 900;
1007     if (fCv == 0x0 ) {
1008       fCv = new TCanvas("cv_event","cv_event",size,size);
1009       fCv->Range(0,0,size,size);
1010     }
1011     if (fTx == 0x0 ) {
1012       fTx = new TText();
1013       fTx->SetTextSize(0.015);
1014       fTx->SetTextColor(kBlack);
1015       fTx->DrawText(0,0,"Bending plane");
1016     }
1017
1018     const float cX = size/2.f; //pt.GetX();
1019     const float cY = size/2.f; //pt.GetY();
1020     const float sc = 8.5f;
1021     if( what.Contains("clusters") ) {
1022       if(fMk==0x0) fMk=new TMarker();
1023       fMk->SetMarkerStyle(20);
1024       fMk->SetMarkerColor(kRed);
1025       fMk->SetMarkerSize(0.3);
1026       for(int iL=0; iL<7; ++iL) {
1027         for(int iC=0; iC<fNClusters[iL];++iC) {
1028           fMk->DrawMarker(cX+fClusters[iL][iC].x*sc,cY+fClusters[iL][iC].y*sc);
1029           if( what.Contains("clusters+id") )  fTx->DrawText(cX+fClusters[iL][iC].x*sc,cY+fClusters[iL][iC].y*sc,Form("%i",fIndex[iL][iC]));
1030         }
1031       }
1032     }
1033     //
1034     if( what.Contains("doublets") ) {
1035       if(fLn==0x0) fLn=new TLine();
1036       for(int iL=0; iL<6; ++iL) {
1037         for(int iD=0; iD<fNDoublets[iL]; ++iD) {
1038           const int id0 = fDoublets[iL][iD].id0;//fIndex[iL][fDoublets[iL][iD].id0];
1039           const int id1 = fDoublets[iL][iD].id1;//fIndex[iL+1][fDoublets[iL][iD].id1];
1040           fLn->DrawLine(cX+fClusters[iL][id0].x*sc,cY+fClusters[iL][id0].y*sc,cX+fClusters[iL+1][id1].x*sc,cY+fClusters[iL+1][id1].y*sc);
1041           if( what.Contains("doublets+level") ) fTx->DrawText(cX+(fClusters[iL][id0].x+fClusters[iL+1][id1].x)*sc/(2.f),cY+(fClusters[iL][id0].y+fClusters[iL+1][id1].y)*sc/(2.f),Form("%i",fDoublets[iL][iD].level));
1042         }
1043       }
1044     }
1045     //
1046   }
1047
1048 //____________________________________________________
1049 void AliITSUTrackerSA::DrawRoads(vector<Road> &vec) {
1050   const int size = 900;
1051   const float cX = size/2.f;
1052   const float cY = size/2.f;
1053   const float sc = 8.5f;
1054
1055   if (fCv == 0x0 ) {
1056     fCv = new TCanvas("cv_event","cv_event",size,size);
1057     fCv->Range(0,0,size,size);
1058   }
1059   if(fLn==0x0) {
1060     fLn=new TLine();
1061     fLn->SetLineColor(kRed);
1062   }
1063   for ( size_t iV = 0; iV < vec.size(); ++iV ) {
1064    cout << vec[iV] << endl;
1065     for( int iE=0; iE<6; ++iE ) {
1066       if (vec[iV].fElements[iE]!=-1) {
1067         const int iD = vec[iV].fElements[iE];
1068         const int id0 = fDoublets[iE][iD].id0;//fIndex[iL][fDoublets[iL][iD].id0];
1069         const int id1 = fDoublets[iE][iD].id1;//fIndex[iE+1][fDoublets[iE][iE].id1];
1070         fLn->DrawLine(cX+fClusters[iE][id0].x*sc,cY+fClusters[iE][id0].y*sc,cX+fClusters[iE+1][id1].x*sc,cY+fClusters[iE+1][id1].y*sc);
1071       }
1072     }
1073   }
1074 }
1075
1076
1077 #endif