CA tracker - updates
[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITSUTrackerSA.cxx
CommitLineData
df380774 1//--------------------------------------------------------------------------------
2cbb5f31 2// Implementation of the ITS tracker class
3// It reads AliITSUClusterPix clusters and and fills the ESD with tracks
df380774 4//
5// The algorithm implemented here takes inspiration from UniCA code of FIAS
6// group.
7//--------------------------------------------------------------------------------
2cbb5f31 8
9#include <TBranch.h>
10#include <TMath.h>
11using TMath::Abs;
12using TMath::Sort;
13using TMath::Sqrt;
14#include <TTree.h>
a0c47fdb 15#include <algorithm>
16using std::sort;
2cbb5f31 17
18// Vc library
19//#include "Vc/Vc"
57a336dd 20//#include "AliITSUTrackerSAauxVc.h" // Structs and other stuff using Vc library
a0c47fdb 21#include "AliLog.h"
2cbb5f31 22#include "AliESDEvent.h"
23#include "AliITSUClusterPix.h"
24#include "AliITSUTrackerSA.h"
a0c47fdb 25#include "AliITSUReconstructor.h"
26#include "AliITSURecoDet.h"
27#include "AliESDtrack.h"
28
29#include <Riostream.h>
30
d44ab054 31using std::cout;
32using std::endl;
33using std::flush;
2cbb5f31 34
57a336dd 35//#include "AliITSUtrackSA.h" // Some dedicated SA track class ?
2cbb5f31 36
37ClassImp(AliITSUTrackerSA)
38
a0c47fdb 39const Double_t AliITSUTrackerSA::fgkToler = 1e-6;// tolerance for layer on-surface check
df380774 40const Double_t AliITSUTrackerSA::fgkChi2Cut = 600.f;
a0c47fdb 41
2cbb5f31 42//________________________________________________________________________________
57a336dd 43AliITSUTrackerSA::AliITSUTrackerSA(AliITSUReconstructor* rec) :
44fReconstructor(rec),
45fITS(0),
46fMatLUT(0),
47fUseMatLUT(kFALSE),
48fCurrMass(0.14),
49//
57a336dd 50fClustersTC(),
df380774 51fChi2Cut( fgkChi2Cut ),
52fPhiCut( 1 ),
53fRPhiCut( 1 ),
54fZCut( 1 )
2cbb5f31 55{
56 //--------------------------------------------------------------------
57 // This default constructor needs to be provided
58 //--------------------------------------------------------------------
a0c47fdb 59 if (rec) Init(rec);
2cbb5f31 60}
61
62//________________________________________________________________________________
df380774 63AliITSUTrackerSA::AliITSUTrackerSA(const AliITSUTrackerSA &t): AliTracker(t),
64 fReconstructor(t.fReconstructor),
65 fITS(t.fITS),
66 fMatLUT(t.fMatLUT),
67 fUseMatLUT(t.fUseMatLUT),
68 fCurrMass(t.fCurrMass),
69 //
70 fClustersTC(),
71 fChi2Cut(fgkChi2Cut),
72 fPhiCut(),
73 fRPhiCut(),
74 fZCut()
2cbb5f31 75{
76 //--------------------------------------------------------------------
77 // The copy constructor is protected
78 //--------------------------------------------------------------------
79}
80
81//________________________________________________________________________________
57a336dd 82AliITSUTrackerSA::~AliITSUTrackerSA()
a0c47fdb 83{
84 // d-tor
85 delete fMatLUT;
86}
87
88
89//_________________________________________________________________________
90void AliITSUTrackerSA::Init(AliITSUReconstructor* rec)
91{
92 // init with external reconstructor
93 //
94 fITS = rec->GetITSInterface();
95 //
96 // create material lookup table
97 const int kNTest = 1000;
98 const double kStepsPerCM=5;
99 fMatLUT = new AliITSUMatLUT(fITS->GetRMin(),fITS->GetRMax(),Nint(kStepsPerCM*(fITS->GetRMax()-fITS->GetRMin())));
100 double zmn = 1e6;
101 for (int ilr=fITS->GetNLayers();ilr--;) {
102 AliITSURecoLayer* lr = fITS->GetLayer(ilr);
103 if (zmn>Abs(lr->GetZMin())) zmn = Abs(lr->GetZMin());
104 if (zmn>Abs(lr->GetZMax())) zmn = Abs(lr->GetZMax());
105 }
106 fMatLUT->FillData(kNTest,-zmn,zmn);
107 //
108}
109
110//________________________________________________________________________________
111Int_t AliITSUTrackerSA::Clusters2Tracks(AliESDEvent *event) {
2cbb5f31 112 //--------------------------------------------------------------------
113 // This is the main tracking function
114 // The clusters must already be loaded
115 //--------------------------------------------------------------------
116
117 // Possibly, create the track "seeds" (combinatorial)
118
119 // Possibly, increment the seeds with additional clusters (Kalman)
120
57a336dd 121 // Possibly, (re)fit the found tracks
2cbb5f31 122
57a336dd 123 // Three iterations:
2cbb5f31 124 // - High momentum first;
57a336dd 125 // - Low momentum with vertex constraint;
126 // - Everything else.
2cbb5f31 127
df380774 128 CellsCreation(0);
129 CellularAutomaton(event);
130 // VertexFinding();
131 // CellsCreation(1);
132 // CellularAutomaton(event);
2cbb5f31 133
134 return 0;
135}
136
137//________________________________________________________________________________
df380774 138Int_t AliITSUTrackerSA::PropagateBack(AliESDEvent * event) {
2cbb5f31 139 //--------------------------------------------------------------------
140 // Here, we implement the Kalman smoother ?
141 // The clusters must be already loaded
142 //--------------------------------------------------------------------
df380774 143 Int_t n=event->GetNumberOfTracks();
144 Int_t ntrk=0;
145 Int_t ngood=0;
146 for (Int_t i=0; i<n; i++) {
147 AliESDtrack *esdTrack=event->GetTrack(i);
148
149 if ((esdTrack->GetStatus()&AliESDtrack::kITSin)==0) continue;
150
151 AliITSUTrackCooked track(*esdTrack);
152
153 track.ResetCovariance(10.);
154
155 int points[2*AliITSUAux::kMaxLayers];
156 for (UInt_t k=0; k<2*AliITSUAux::kMaxLayers; k++)
157 points[k]=-1;
158 Int_t nc=track.GetNumberOfClusters();
159 for (Int_t k=0; k<nc; k++) {
160 const int layer = (track.GetClusterIndex(k)&0xf0000000)>>28;
161 const int idx = (track.GetClusterIndex(k)&0x0fffffff);
162 points[layer<<1]=idx;
163 }
164
165 if (RefitTrack(&track,points,40,1)>=0) {
166
167 CookLabel(&track, 0.); //For comparison only
168 Int_t label=track.GetLabel();
169 if (label>0) ngood++;
170
171 esdTrack->UpdateTrackParams(&track,AliESDtrack::kITSout);
172 ntrk++;
173 }
174 }
175
176 Info("PropagateBack","Back propagated tracks: %d",ntrk);
177 if (ntrk)
178 Info("PropagateBack","Good tracks/back propagated: %f",Float_t(ngood)/ntrk);
2cbb5f31 179
180 return 0;
181}
182
183//________________________________________________________________________________
df380774 184Int_t AliITSUTrackerSA::RefitInward(AliESDEvent * event) {
2cbb5f31 185 //--------------------------------------------------------------------
57a336dd 186 // Some final refit, after the outliers get removed by the smoother ?
2cbb5f31 187 // The clusters must be loaded
188 //--------------------------------------------------------------------
df380774 189 Int_t n=event->GetNumberOfTracks();
190 Int_t ntrk=0;
191 Int_t ngood=0;
192 for (Int_t i=0; i<n; i++) {
193 AliESDtrack *esdTrack=event->GetTrack(i);
194
195 if ((esdTrack->GetStatus()&AliESDtrack::kITSout)==0) continue;
196
197 AliITSUTrackCooked track(*esdTrack);
198
199 track.ResetCovariance(10.);
200
201 int points[2*AliITSUAux::kMaxLayers];
202 for (UInt_t k=0; k<2*AliITSUAux::kMaxLayers; k++)
203 points[k]=-1;
204 Int_t nc=track.GetNumberOfClusters();
205 for (Int_t k=0; k<nc; k++) {
206 const int layer = (track.GetClusterIndex(k)&0xf0000000)>>28;
207 const int idx = (track.GetClusterIndex(k)&0x0fffffff);
208 points[layer<<1]=idx;
209 }
210
211 if (RefitTrack(&track,points,1.8,1)>=0) { //2.1,1)>=0) {
212
213 //if (!track.PropagateTo(1.8, 2.27e-3, 35.28*1.848)) continue;
214 CookLabel(&track, 0.); //For comparison only
215 Int_t label=track.GetLabel();
216 if (label>0) ngood++;
217
218 //cout << esdTrack->GetStatus() << " ";
219 esdTrack->UpdateTrackParams(&track,AliESDtrack::kITSrefit);
220 //cout << esdTrack->GetStatus() << endl;
221 ntrk++;
222 }
223 }
2cbb5f31 224
df380774 225 Info("RefitInward","Refitted tracks: %d",ntrk);
226 if (ntrk)
227 Info("RefitInward","Good tracks/refitted: %f",Float_t(ngood)/ntrk);
228
2cbb5f31 229 return 0;
230}
231
57a336dd 232//________________________________________________________________________________
2cbb5f31 233Int_t AliITSUTrackerSA::LoadClusters(TTree *cluTree) {
234 //--------------------------------------------------------------------
235 // This function reads the ITSU clusters from the tree,
236 // sort them, distribute over the internal tracker arrays, etc
237 //--------------------------------------------------------------------
a0c47fdb 238 fITS->LoadClusters(cluTree);
239 fITS->ProcessClusters();
240 //
2cbb5f31 241 for(int iL=0; iL<7; ++iL) {
a0c47fdb 242 fClustersTC[iL]=*fITS->GetLayerActive(iL)->GetClustersAddress();
57a336dd 243 AliITSURecoLayer* lr = fITS->GetLayerActive(iL) ; // assign the layer which the cluster belongs to
df380774 244 for(int iC=0;iC<fClustersTC[iL]->GetEntriesFast();++iC) {
245 const AliITSUClusterPix *cl = (AliITSUClusterPix*)fClustersTC[iL]->At(iC);
2cbb5f31 246 float pos[3];
247 cl->GetGlobalXYZ(pos);
df380774 248 float phi = TMath::PiOver2();
249 if(Abs(pos[0])>1e-9) {
250 phi=TMath::ATan2(pos[1]-GetY(),pos[0]-GetX());
251 if(phi<0.f) phi+=TMath::TwoPi();
252 } else if(pos[1]<0.f) phi *= 3.f ;
57a336dd 253 AliITSURecoSens* sens = lr->GetSensorFromID(cl->GetVolumeId());
df380774 254 const float alpha = sens->GetPhiTF();
255 const float cov[3]={cl->GetSigmaZ2(),cl->GetSigmaYZ(),cl->GetSigmaY2()};
256
257 fLayer[iL].AddPoint(pos,cov,phi,alpha);
258 for ( int i=0 ; i<3; ++i ) {
259 fLayer[iL].Points.back().Label[i] = (cl->GetLabel(i)<0) ? -1 : cl->GetLabel(i);
260 }
2cbb5f31 261 }
df380774 262
263 fLayer[iL].Sort();
264
2cbb5f31 265 }
2cbb5f31 266 return 0;
267}
268
269//________________________________________________________________________________
270void AliITSUTrackerSA::UnloadClusters() {
271 //--------------------------------------------------------------------
272 // This function unloads ITSU clusters from the RAM
273 //--------------------------------------------------------------------
df380774 274 for(int i=0;i<7;++i)
275 fLayer[i].Clear();
276 for(int i=0;i<5;++i)
277 fCells[i].clear();
2cbb5f31 278
2cbb5f31 279}
280
281//________________________________________________________________________________
df380774 282void AliITSUTrackerSA::CellularAutomaton(AliESDEvent *event) {
283
2cbb5f31 284 // Here it's implemented the Cellular Automaton routine
57a336dd 285 // Firstly the level of each doublet is set according to the level of
286 // the neighbour doublets.
287 // Doublet are considered to be neighbour if they share one point and the
2cbb5f31 288 // phi and theta direction difference of the two is below a cut value.
df380774 289 Int_t ntrk=0,ngood=0;
290
291 for( int iL = 1; iL < 5; ++iL ) {
292 for( size_t iC1 = 0; iC1 < fCells[iL-1].size(); ++iC1 ) {
293 for( size_t iC2 = 0; iC2 < fCells[iL].size(); ++iC2 ) {
294 if( fCells[iL-1][iC1].Points[1]==fCells[iL][iC2].Points[0] &&
295 fCells[iL-1][iC1].Points[2]==fCells[iL][iC2].Points[1] &&
296 fCells[iL-1][iC1].Level >= fCells[iL][iC2].Level - 1 ) {
297 // The implementation of the curvature based matching has to be studied.
298 fCells[iL][iC2].Level = fCells[iL-1][iC1].Level+1;
299 fCells[iL][iC2].Neighbours.push_back(iC1);
57a336dd 300 }
2cbb5f31 301 }
302 }
303 }
df380774 304
305 for (int level = 5; level > 0; --level ) {
306 vector<Road> roads;
307 roads.reserve(100); // should reserve() be based on number of clusters on outermost layer?
308 for (int iCL=4; iCL >= level-1; --iCL ) {
309 for (size_t iCell = 0; iCell < fCells[iCL].size(); ++iCell) {
310 if ( fCells[iCL][iCell].Level != level ) continue;
311 roads.push_back( Road(iCL,iCell) );
312 for( size_t iN=0;iN<fCells[iCL][iCell].Neighbours.size(); ++iN ) {
313 const int currD = iCL - 1;
314 const int neigh = fCells[iCL][iCell].Neighbours[iN];
315 if( iN > 0 ) roads.push_back(roads.back());
316 CandidatesTreeTraversal(roads,neigh,currD);
57a336dd 317 }
df380774 318 fCells[iCL][iCell].Level = -1;
a0c47fdb 319 }
320 }
321
df380774 322 // for(size_t iR=0; iR<roads.size(); ++iR) {
323 // cout << "ROAD " << iR << " | ";
324 // for(int i=0;i<5;++i) {
325 // if(roads[iR][i]<0) continue;
326 // else {
327 // if(roads[iR].Label==-1){
328 // roads[iR].Label = fCells[i][roads[iR][i]].Label;
329 // if(roads[iR].Label==-1) roads[iR].Label--;
330 // }
331 // if (fCells[i][roads[iR][i]].Label!=roads[iR].Label&&roads[iR].Label>-1) {
332 // roads[iR].Label = -1;
333 // if(fCells[i][roads[iR][i]].Label==-1) roads[iR].Label--;
334 // }
335
336 // cout << fCells[i][roads[iR][i]].Label << " ";
337 // }
338 // }
339 // cout << " | " << roads[iR].Label << " | " << roads[iR].N << endl;
340 // }
341 vector<AliITSUTrackCooked> candidates;
342 candidates.reserve(roads.size());
343
344 for (size_t iR=0; iR<roads.size(); ++iR) {
345 if(roads[iR].N != level) {
57a336dd 346 continue;
a0c47fdb 347 }
df380774 348
349 int points[2*AliITSUAux::kMaxLayers];
350 for(unsigned int i=0;i<2*AliITSUAux::kMaxLayers;++i) points[i] = -1;
351 for(int i=0;i<5;++i) {
352 if(roads[iR].Elements[i]<0) continue;
353 points[( i )<<1]=fLayer[ i ](fCells[i][roads[iR].Elements[i]].Points[0]);
354 points[(i+1)<<1]=fLayer[i+1](fCells[i][roads[iR].Elements[i]].Points[1]);
355 points[(i+2)<<1]=fLayer[i+2](fCells[i][roads[iR].Elements[i]].Points[2]);
356 }
357
358 candidates.push_back(AliITSUTrackCooked());
359
360 InitTrackParams(candidates.back(),points);
361 const double chi2 = RefitTrack( (AliExternalTrackParam*)&candidates.back(), points, 0. ,-1);
362
363 if ( chi2 < 0. ) {
364 // cout << "FAIL: " << chi2 << endl;
365 // for(unsigned int i=0;i<2*AliITSUAux::kMaxLayers;++i)
366 // cout << points[i] << " ";
367 // cout << endl;
368 candidates.back().SetChi2( 1e27 );
369 } else candidates.back().SetChi2( chi2 );
370 candidates.back().SetLabel(roads[iR].Label);
a0c47fdb 371 }
df380774 372
373 vector<int> index;
374 index.reserve(candidates.size());
375 for ( size_t i = 0; i < candidates.size(); ++i ) index.push_back(i);
376 Comparison<AliITSUTrackCooked> comp(&candidates);
377 sort(index.begin(),index.end(),comp);
378
57a336dd 379 for ( size_t cand = 0; cand < candidates.size(); ++cand ) {
a0c47fdb 380 const int ii = index[cand];
df380774 381
382 if ( candidates[ii].GetChi2() < 0. ) continue;
383
384 // cout << candidates[ii].GetChi2() << " " << candidates[ii].GetNumberOfClusters() << " | " << candidates[ii].GetLabel() << " | ";
385 // for(int i=0;i<candidates[ii].GetNumberOfClusters();++i) {
386 // cout<< GetCluster(candidates[ii].GetClusterIndex(i))->GetLabel(0) << " ";
387 // }
388 // cout << endl;
389
390 if( candidates[ii].GetChi2()/candidates[ii].GetNumberOfClusters() > fgkChi2Cut ) {
391 break;
392 }
a0c47fdb 393 bool goodTrack = true;
df380774 394 for ( int point = 0; point < candidates[ii].GetNumberOfClusters(); ++point ) {
395 int layer = (candidates[ii].GetClusterIndex(point)&0xf0000000)>>28;
396 int ind = (candidates[ii].GetClusterIndex(point)&0x0fffffff);
397
398 if( (fLayer[ layer ].Points[ ind ].Used ) ) {
399 goodTrack = false;
57a336dd 400 }
df380774 401
2cbb5f31 402 }
df380774 403 if(!goodTrack) {
404 continue;
405 }
406 for ( int point = 0; point < candidates[ii].GetNumberOfClusters(); ++point ) {
407 int layer = (candidates[ii].GetClusterIndex(point)&0xf0000000)>>28;
408 int ind = (candidates[ii].GetClusterIndex(point)&0x0fffffff);
409 fLayer[ layer ].Points[ ind ].Used = true;
a0c47fdb 410 }
a0c47fdb 411
57a336dd 412 AliESDtrack outTrack;
df380774 413 CookLabel((AliKalmanTrack*)&candidates[ii],0.f);
414 ntrk++;
415 if(candidates[ii].GetChi2()>0) ngood++;
416
417 // cout << candidates[ii].GetChi2() << " " << candidates[ii].GetNumberOfClusters() << " | " << candidates[ii].GetLabel() << " | ";
418 // for(int i=0;i<candidates[ii].GetNumberOfClusters();++i) {
419 // cout<< GetCluster(candidates[ii].GetClusterIndex(i))->GetLabel(0) << " ";
420 // }
421 // cout << endl;
422
423 outTrack.UpdateTrackParams((AliKalmanTrack*)&candidates[ii],AliESDtrack::kITSin);
424 outTrack.SetLabel(candidates[ii].GetLabel());
57a336dd 425 event->AddTrack(&outTrack);
426 }
427 }
df380774 428 Info("Clusters2Tracks","Reconstructed tracks: %d",ntrk);
429 if (ntrk)
430 Info("Clusters2Tracks","Good tracks/reconstructed: %f",Float_t(ngood)/ntrk);
2cbb5f31 431}
432
df380774 433
2cbb5f31 434//________________________________________________________________________________
df380774 435void AliITSUTrackerSA::CellsCreation(const int &cutLevel) {
57a336dd 436 // Make associations between two points on adjacent layers within an azimuthal window.
2cbb5f31 437 // Under consideration:
438 // - track parameter estimation using the primary vertex position
439 // To do:
440 // - last iteration
441
df380774 442 float phiCut = 7.f;
443 if( cutLevel==0 ) phiCut = fPhiCut;
a0c47fdb 444
df380774 445 // Doublets creation
446 vector<Cell> doublets[6];
2cbb5f31 447 for( int iL = 0 ; iL < 6 ; ++iL ) {
df380774 448 for ( int iC1 = 0 ; iC1 < fLayer[iL].N ; ++iC1 ) {
449 for ( int iC2 = 0; iC2 < fLayer[iL+1].N ; ++iC2 ) {
450 const float dPhi = Abs( fLayer[iL][iC1].Phi - fLayer[iL+1][iC2].Phi );
451 if( dPhi < phiCut || Abs( dPhi - TMath::TwoPi() ) < phiCut) {
452 doublets[iL].push_back(Cell(iC1,iC2));
453 if(Abs(fLayer[iL][iC1].XYZ[0]-fLayer[iL+1][iC2].XYZ[0])<1e-32) {
454 doublets[iL].back().Param[0] = 1e32;
455 } else {
456 doublets[iL].back().Param[0] = (fLayer[iL][iC1].XYZ[1]-fLayer[iL+1][iC2].XYZ[1])/(fLayer[iL][iC1].XYZ[0]-fLayer[iL+1][iC2].XYZ[0]);
457 }
458 const float r1 = Sqrt(fLayer[iL][iC1].XYZ[0] * fLayer[iL][iC1].XYZ[0] + fLayer[iL][iC1].XYZ[1] * fLayer[iL][iC1].XYZ[1]);
459 const float r2 = Sqrt(fLayer[iL+1][iC2].XYZ[0] * fLayer[iL+1][iC2].XYZ[0] + fLayer[iL+1][iC2].XYZ[1] * fLayer[iL+1][iC2].XYZ[1]);
460 doublets[iL].back().Param[1] = (fLayer[iL][iC1].XYZ[2]-fLayer[iL+1][iC2].XYZ[2])/(r1-r2);
461 doublets[iL].back().Label=-1;
462 for(int i=0;i<3;++i) {
463 for(int j=0;j<3;++j) {
464 if(fLayer[iL][iC1].Label[i]>-1&&fLayer[iL][iC1].Label[i]==fLayer[iL+1][iC2].Label[j])
465 doublets[iL].back().Label = fLayer[iL][iC1].Label[i];
466 }
467 }
468 } else if( fLayer[iL+1][iC2].Phi - fLayer[iL][iC1].Phi > phiCut ) break;
2cbb5f31 469 }
470
471 }
df380774 472 }
2cbb5f31 473
df380774 474 // Triplets creation
475 for( int iL = 5; iL > 0; --iL ) {
476 fCells[iL-1].clear();
477 for ( size_t iD2 = 0; iD2 < doublets[iL].size(); ++iD2 ) {
478 for ( size_t iD1 = 0; iD1 < doublets[iL-1].size(); ++iD1 ) {
479 const int id1 = doublets[iL-1][iD1].Points[1];
480 const int id2 = doublets[iL][iD2].Points[0];
481 if ( id1 == id2 ) {
482 const int id3 = doublets[iL][iD2].Points[1];
483 const float r3 = Sqrt( fLayer[iL+1][id3].XYZ[0] * fLayer[iL+1][id3].XYZ[0] + fLayer[iL+1][id3].XYZ[1] * fLayer[iL+1][id3].XYZ[1] );
484 const float r2 = Sqrt( fLayer[iL][id2].XYZ[0] * fLayer[iL][id2].XYZ[0] + fLayer[iL][id2].XYZ[1] * fLayer[iL][id2].XYZ[1] );
485 const float extrZ3 = doublets[iL-1][iD1].Param[1] * ( r3 - r2 ) + fLayer[iL][id2].XYZ[2] ;
486 const int iii = doublets[iL-1][iD1].Points[0];
487 if ( Abs ( extrZ3 - fLayer[iL+1][id3].XYZ[2] ) < fZCut ) {
488 fCells[iL-1].push_back(Cell(doublets[iL-1][iD1].Points[0],id2,id3));
489 fCells[iL-1].back().Param[0] = Curvature(fLayer[iL+1][id3].XYZ[0],fLayer[iL+1][id3].XYZ[1],fLayer[iL][id2].XYZ[0],fLayer[iL][id2].XYZ[1],fLayer[iL-1][iii].XYZ[0],fLayer[iL-1][iii].XYZ[1]);
490 fCells[iL-1].back().Param[1] = doublets[iL][iD2].Param[1];
491 if(doublets[iL-1][iD1].Label==doublets[iL][iD2].Label&&doublets[iL][iD2].Label!=-1)
492 fCells[iL-1].back().Label=doublets[iL][iD2].Label;
493 else
494 fCells[iL-1].back().Label=-1;
495 }
496 }
2cbb5f31 497 }
2cbb5f31 498 }
2cbb5f31 499 }
df380774 500
a0c47fdb 501}
df380774 502
a0c47fdb 503//______________________________________________________________________________
df380774 504Bool_t AliITSUTrackerSA::InitTrackParams(AliITSUTrackCooked &track, int points[])
a0c47fdb 505{
506 // Set the initial guess on track kinematics for propagation.
507 // Assume at least 3 points available
508 int lrOcc[AliITSUAux::kMaxLayers], nCl=0;
509 //
510 // we will need endpoints and middle layer
df380774 511 for (int i=fITS->GetNLayersActive()-1; i>=0; i--) {
512 if (points[i<<0x1]>-1) {
57a336dd 513 lrOcc[nCl++] = i;
df380774 514 track.SetClusterIndex(i,points[i<<0x1]);
57a336dd 515 }
516 }
df380774 517
a0c47fdb 518 if (nCl<3) {
519 AliError(Form("Cannot estimate momentum of tracks with %d clusters",nCl));
a0c47fdb 520 return kFALSE;
521 }
522 //
df380774 523 const int lr0 = lrOcc[0];
524 const int lr1 = lrOcc[nCl/2];
525 const int lr2 = lrOcc[nCl-1];
a0c47fdb 526 //
df380774 527 const SpacePoint& cl0 = fLayer[lr0].Points[ points[lr0<<1] ];
528 const SpacePoint& cl1 = fLayer[lr1].Points[ points[lr1<<1] ];
529 const SpacePoint& cl2 = fLayer[lr2].Points[ points[lr2<<1] ];
530 double cv = Curvature(cl0.XYZ[0],cl0.XYZ[1], cl1.XYZ[0],cl1.XYZ[1], cl2.XYZ[0],cl2.XYZ[1]);
531
532 double tgl = (cl2.XYZ[2]-cl0.XYZ[2])/TMath::Sqrt((cl2.XYZ[0]-cl0.XYZ[0])*(cl2.XYZ[0]-cl0.XYZ[0])+(cl2.XYZ[1]-cl0.XYZ[1])*(cl2.XYZ[1]-cl0.XYZ[1]));
a0c47fdb 533 //
df380774 534 AliITSUClusterPix* clus = (AliITSUClusterPix*)fClustersTC[ lr0 ]->At( points[lr0<<1] );
a0c47fdb 535 AliITSURecoLayer* lr = fITS->GetLayerActive(lr0);
536 AliITSURecoSens* sens = lr->GetSensorFromID(clus->GetVolumeId());
537 double x = sens->GetXTF() + clus->GetX();
538 double alp = sens->GetPhiTF();
539 // printf("Alp: %f phi: %f\n",alp,phi);
540 double par[5] = {clus->GetY(),clus->GetZ(),0,tgl,cv};
541 double cov[15] = {
542 5*5,
df380774 543 0, 5*5,
544 0, 0 , 0.7*0.7,
545 0, 0, 0, 0.7*0.7,
546 0, 0, 0, 0, 10
a0c47fdb 547 };
548 track.Set(x,alp,par,cov);
df380774 549
a0c47fdb 550 return kTRUE;
551}
552
553//______________________________________________________________________________
57a336dd 554void AliITSUTrackerSA::CandidatesTreeTraversal(vector<Road> &candidates, const int &iD, const int &doubl) {
a0c47fdb 555
df380774 556 if ( doubl < 0 ) return;
557
57a336dd 558 candidates.back().AddElement(doubl,iD);
df380774 559 const int currentN = candidates.back().N;
560 for ( size_t iN = 0; iN < fCells[doubl][iD].Neighbours.size(); ++iN ) {
57a336dd 561 const int currD = doubl - 1 ;
df380774 562 const int neigh = fCells[doubl][iD].Neighbours[iN];
563
57a336dd 564 if ( iN > 0 ) {
57a336dd 565 candidates.push_back(static_cast<Road>(candidates.back()));
df380774 566 candidates.back().N = currentN;
a0c47fdb 567 }
df380774 568
57a336dd 569 CandidatesTreeTraversal(candidates,neigh,currD);
a0c47fdb 570 }
57a336dd 571
df380774 572 fCells[doubl][iD].Level = -1;
57a336dd 573
a0c47fdb 574}
575
576//______________________________________________________________________________
df380774 577Double_t AliITSUTrackerSA::RefitTrack(AliExternalTrackParam* trc,
578 Int_t clInfo[2*AliITSUAux::kMaxLayers],
579 Double_t rDest, Int_t stopCond)
a0c47fdb 580{
df380774 581 // refit track till radius rDest.
a0c47fdb 582 // if stopCond<0 : propagate till last cluster then stop
583 // if stopCond==0: propagate till last cluster then try to go till limiting rDest, don't mind if fail
584 // if stopCond>0 : rDest must be reached
585 //
586 // The clList should provide the indices of clusters at corresponding layer (as stored in the layer
587 // TClonesArray, with convention (allowing for up to 2 clusters per layer due to the overlaps):
588 // if there is a cluster on given layer I, then it should be stored at clInfo[2*I-1]
589 // if there is an additional cluster on this layer, it goes to clInfo[2*I]
590 // -1 means no cluster
591 //
592 double rCurr = Sqrt(trc->GetX()*trc->GetX() + trc->GetY()*trc->GetY());
593 int dir,lrStart,lrStop;
594 //
595 dir = rCurr<rDest ? 1 : -1;
596 lrStart = fITS->FindFirstLayerID(rCurr,dir);
597 lrStop = fITS->FindLastLayerID(rDest,dir); // lr id before which we have to stop
a0c47fdb 598 //
599 if (lrStop<0 || lrStart<0) AliFatal(Form("Failed to find start(%d) or last(%d) layers. "
df380774 600 "Track from %.3f to %.3f",lrStart,lrStop,rCurr,rDest));
a0c47fdb 601 //
602 int nCl = 0;
df380774 603 for (int i=2*fITS->GetNLayersActive();i--;) {if (clInfo[i]<0) continue; nCl++;}
a0c47fdb 604 //
605 AliExternalTrackParam tmpTr(*trc);
606 double chi2 = 0;
607 int iclLr[2],nclLr;
608 int nclFit = 0;
609 //
610 int lrStop1 = lrStop+dir;
611 for (int ilr=lrStart;ilr!=lrStop1;ilr+=dir) {
612 AliITSURecoLayer* lr = fITS->GetLayer(ilr);
df380774 613 if ( dir*(rCurr-lr->GetR(dir))>0) continue; // this layer is already passed
a0c47fdb 614 int ilrA2,ilrA = lr->GetActiveID();
615 // passive layer or active w/o hits will be traversed on the way to next cluster
df380774 616 if (!lr->IsActive() || clInfo[ilrA2=(ilrA<<1)]<0) continue;
a0c47fdb 617 //
618 // select the order in which possible 2 clusters (in case of the overlap) will be traversed and fitted
619 nclLr=0;
620 if (dir>0) { // clusters are stored in increasing radius order
621 iclLr[nclLr++]=clInfo[ilrA2++];
622 if (clInfo[ilrA2]>=0) iclLr[nclLr++]=clInfo[ilrA2];
623 }
624 else {
625 if ( clInfo[ilrA2+1]>=0 ) iclLr[nclLr++]=clInfo[ilrA2+1];
626 iclLr[nclLr++]=clInfo[ilrA2];
627 }
628 //
629 Bool_t transportedToLayer = kFALSE;
630 for (int icl=0;icl<nclLr;icl++) {
631 AliITSUClusterPix* clus = (AliITSUClusterPix*)lr->GetCluster(iclLr[icl]);
632 AliITSURecoSens* sens = lr->GetSensorFromID(clus->GetVolumeId());
df380774 633 if (!tmpTr.Rotate(sens->GetPhiTF())) return -1;
a0c47fdb 634 //
635 double xClus = sens->GetXTF()+clus->GetX();
636 if (!transportedToLayer) {
df380774 637 if (ilr!=lrStart && !TransportToLayerX(&tmpTr,lrStart,ilr,xClus)) return -1; // go to the entrance to the layer
638 lrStart = ilr;
639 transportedToLayer = kTRUE;
a0c47fdb 640 }
641 //
df380774 642 if (!PropagateSeed(&tmpTr,xClus,fCurrMass)) return -1;
57a336dd 643 //
df380774 644 Double_t p[2]={clus->GetY(), clus->GetZ()};
645 Double_t cov[3]={clus->GetSigmaY2(), clus->GetSigmaYZ(), clus->GetSigmaZ2()};
646 double chi2cl = tmpTr.GetPredictedChi2(p,cov);
647 chi2 += chi2cl;
a0c47fdb 648 //
df380774 649 if ( !tmpTr.Update(p,cov) ) return -1;
650 if (++nclFit==nCl && stopCond<0) {
651 *trc = tmpTr;
652 return chi2; // it was requested to not propagate after last update
a0c47fdb 653 }
654 }
a0c47fdb 655 //
a0c47fdb 656 }
df380774 657 // All clusters were succesfully fitted. Even if the track does not reach rDest, this is enough to validate it.
658 // Still, try to go as close as possible to rDest.
659 //
660 if (lrStart!=lrStop) {
661 if (!TransportToLayer(&tmpTr,lrStart,lrStop)) return (stopCond>0) ? -chi2 : chi2; // rDest was obligatory
662 if (!GoToExitFromLayer(&tmpTr,fITS->GetLayer(lrStop),dir)) return (stopCond>0) ? -chi2 : chi2; // rDest was obligatory
663 }
664 // go to the destination radius. Note that here we don't select direction to avoid precision problems
665 if (!tmpTr.GetXatLabR(rDest,rDest,GetBz(),0) || !PropagateSeed(&tmpTr,rDest,fCurrMass, 100, kFALSE)) {
666 return (stopCond>0) ? -chi2 : chi2; // rDest was obligatory
667 }
668 *trc = tmpTr;
a0c47fdb 669
df380774 670 return chi2;
671}
672
673//______________________________________________________________________________
674Bool_t AliITSUTrackerSA::PropagateSeed(AliExternalTrackParam *seed, Double_t xToGo, Double_t mass, Double_t maxStep, Bool_t matCorr)
675{
676 // propagate seed to given x applying material correction if requested
677 const Double_t kEpsilon = 1e-5;
678 Double_t xpos = seed->GetX();
679 Int_t dir = (xpos<xToGo) ? 1:-1;
680 Double_t xyz0[3],xyz1[3];
681 //
682 Bool_t updTime = dir>0 && seed->IsStartedTimeIntegral();
683 if (matCorr || updTime) seed->GetXYZ(xyz1); //starting global position
684 while ( (xToGo-xpos)*dir > kEpsilon){
685 Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
686 Double_t x = xpos+step;
687 Double_t bz=GetBz(); // getting the local Bz
688 if (!seed->PropagateTo(x,bz)) return kFALSE;
689 double ds = 0;
690 if (matCorr || updTime) {
691 xyz0[0]=xyz1[0]; // global pos at the beginning of step
692 xyz0[1]=xyz1[1];
693 xyz0[2]=xyz1[2];
694 seed->GetXYZ(xyz1); // // global pos at the end of step
57a336dd 695 //
df380774 696 if (matCorr) {
697 Double_t xrho,xx0;
698 ds = GetMaterialBudget(xyz0,xyz1,xx0,xrho);
699 if (dir>0) xrho = -xrho; // outward should be negative
700 if (!seed->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
a0c47fdb 701 }
df380774 702 else { // matCorr is not requested but time integral is
703 double d0 = xyz1[0]-xyz0[0];
704 double d1 = xyz1[1]-xyz0[1];
705 double d2 = xyz1[2]-xyz0[2];
706 ds = TMath::Sqrt(d0*d0+d1*d1+d2*d2);
57a336dd 707 }
a0c47fdb 708 }
df380774 709 if (updTime) seed->AddTimeStep(ds);
57a336dd 710 //
df380774 711 xpos = seed->GetX();
a0c47fdb 712 }
df380774 713 return kTRUE;
714}
a0c47fdb 715
df380774 716//_________________________________________________________________________
717Bool_t AliITSUTrackerSA::TransportToLayer(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo, Double_t rLim)
718{
719 // transport track from layerFrom to the entrance of layerTo or to rLim (if>0), wathever is closer
720 //
721 if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
722 //
723 int dir = lTo > lFrom ? 1:-1;
724 AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
725 Bool_t checkFirst = kTRUE;
726 Bool_t limReached = kFALSE;
727 while(lFrom!=lTo) {
728 if (lrFr) {
729 if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) return kFALSE; // go till the end of current layer
730 checkFirst = kFALSE;
731 }
732 AliITSURecoLayer* lrTo = fITS->GetLayer( (lFrom+=dir) );
733 if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
a0c47fdb 734 //
df380774 735 // go the entrance of the layer, assuming no materials in between
736 double xToGo = lrTo->GetR(-dir);
737 if (rLim>0) {
738 if (dir>0) {
739 if (rLim<xToGo) {xToGo = rLim; limReached = kTRUE;}
740 }
741 else {
742 if (rLim>xToGo) {xToGo = rLim; limReached = kTRUE;}
743 }
a0c47fdb 744 }
df380774 745 // double xts = xToGo;
746 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) {
747 // printf("FailHere1: %f %f %d\n",xts,xToGo,dir);
748 // seed->Print("etp");
749 return kFALSE;
57a336dd 750 }
df380774 751 if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) {
752 //printf("FailHere2: %f %f %d\n",xts,xToGo,dir);
753 //seed->Print("etp");
754 return kFALSE;
57a336dd 755 }
df380774 756 lrFr = lrTo;
757 if (limReached) break;
a0c47fdb 758 }
df380774 759 return kTRUE;
760 //
761}
a0c47fdb 762
df380774 763//_________________________________________________________________________
764Bool_t AliITSUTrackerSA::TransportToLayerX(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo, Double_t xStop)
765{
766 // transport track from layerFrom to the entrance of layerTo but do not pass control parameter X
767 //
768 if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
769 //
770 int dir = lTo > lFrom ? 1:-1;
771 AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
772 Bool_t checkFirst = kTRUE;
773 while(lFrom!=lTo) {
774 if (lrFr) {
775 if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) return kFALSE; // go till the end of current layer
776 checkFirst = kFALSE;
777 }
778 AliITSURecoLayer* lrTo = fITS->GetLayer( (lFrom+=dir) );
779 if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
57a336dd 780 //
df380774 781 // go the entrance of the layer, assuming no materials in between
782 double xToGo = lrTo->GetR(-dir); // R of the entrance to layer
783 //
784 // double xts = xToGo;
785 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) {
786 // printf("FailHere1: %f %f %d\n",xts,xToGo,dir);
787 // seed->Print("etp");
788 return kFALSE;
57a336dd 789 }
df380774 790 if ( (dir>0&&xToGo>xStop) || (dir<0&&xToGo<xStop) ) xToGo = xStop;
57a336dd 791 //
df380774 792#ifdef _ITSU_DEBUG_
793 AliDebug(2,Form("go in dir=%d to R=%.4f(X:%.4f)",dir,lrTo->GetR(-dir), xToGo));
794#endif
795 if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) {
796 //printf("FailHere2: %f %f %d\n",xts,xToGo,dir);
797 //seed->Print("etp");
798 return kFALSE;
57a336dd 799 }
df380774 800 lrFr = lrTo;
a0c47fdb 801 }
df380774 802 return kTRUE;
803 //
804}
a0c47fdb 805
df380774 806//_________________________________________________________________________
807Bool_t AliITSUTrackerSA::GoToExitFromLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
808{
809 // go to the exit from lr in direction dir, applying material corrections in steps specific for this layer
810 // If check is requested, do this only provided the track has not exited the layer already
811 double xToGo = lr->GetR(dir);
812 if (check) { // do we need to track till the surface of the current layer ?
813 double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
814 if (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // on the surface or outside of the layer
815 else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // on the surface or outside of the layer
816 }
817 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
818 // go via layer to its boundary, applying material correction.
819 if (!PropagateSeed(seed,xToGo,fCurrMass, lr->GetMaxStep())) return kFALSE;
820 //
821 return kTRUE;
822 //
823}
a0c47fdb 824
df380774 825//_________________________________________________________________________
826Bool_t AliITSUTrackerSA::GoToEntranceToLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
827{
828 // go to the entrance of lr in direction dir, w/o applying material corrections.
829 // If check is requested, do this only provided the track did not reach the layer already
830 double xToGo = lr->GetR(-dir);
831 if (check) { // do we need to track till the surface of the current layer ?
832 double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
833 if (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // already passed
834 else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // already passed
57a336dd 835 }
df380774 836 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
837 // go via layer to its boundary, applying material correction.
838 if (!PropagateSeed(seed,xToGo,fCurrMass, 100, kFALSE)) return kFALSE;
839 return kTRUE;
840 //
841}
57a336dd 842
843//____________________________________________________
df380774 844Double_t AliITSUTrackerSA::GetMaterialBudget(const double* pnt0,const double* pnt1, double& x2x0, double& rhol) const
845{
846 double par[7];
847 if (fUseMatLUT && fMatLUT) {
848 double d = fMatLUT->GetMatBudget(pnt0,pnt1,par);
849 x2x0 = par[AliITSUMatLUT::kParX2X0];
850 rhol = par[AliITSUMatLUT::kParRhoL];
851 return d;
a0c47fdb 852 }
df380774 853 else {
854 MeanMaterialBudget(pnt0,pnt1,par);
855 x2x0 = par[1];
856 rhol = par[0]*par[4];
857 return par[4];
57a336dd 858 }
a0c47fdb 859}
860
df380774 861//____________________________________________________________________
862Double_t AliITSUTrackerSA::Curvature(Double_t x1,Double_t y1,Double_t x2,Double_t y2,Double_t x3,Double_t y3) {
863
864 //calculates the curvature of track
865 Double_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
866 if(den*den<1e-32) return 0.;
867 Double_t a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
868 if((y2-y1)*(y2-y1)<1e-32) return 0.;
869 Double_t b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
870 Double_t c = -x1*x1-y1*y1-a*x1-b*y1;
871 Double_t xc= -a/2.;
872
873 if((a*a+b*b-4*c)<0) return 0.;
874 Double_t rad = TMath::Sqrt(a*a+b*b-4*c)/2.;
875 if(rad*rad<1e-32) return 1e16;
876
877 if((x1>0 && y1>0 && x1<xc)) rad*=-1;
878 if((x1<0 && y1>0 && x1<xc)) rad*=-1;
879 // if((x1<0 && y1<0 && x1<xc)) rad*=-1;
880 // if((x1>0 && y1<0 && x1<xc)) rad*=-1;
881
882 return 1/rad;
883
884}
57a336dd 885