]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/UPGRADE/AliITStrackerUpgrade.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
moving convcorr files to gammaconv dir
[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITStrackerUpgrade.cxx
CommitLineData
1d9af2d5 1/**************************************************************************
2 * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
18////////////////////////////////////////////////////
19// ITS Upgrade Stand alone tracker class //
20// Authors: A.Mastroserio C.Terrevoli //
21// e-mail:annalisa.mastroserio@cern.ch //
22// cristina.terrevoli@ba.infn.it //
23////////////////////////////////////////////////////
24
25#include <stdlib.h>
26
27#include <TArrayI.h>
28#include <TBranch.h>
29#include <TObjArray.h>
30#include <TTree.h>
31#include "AliRunLoader.h"
32#include "AliITSLoader.h"
33#include "AliESDEvent.h"
34#include "AliESDVertex.h"
35#include "AliESDtrack.h"
36#include "AliITSVertexer.h"
37#include "AliITSclusterTable.h"
0ac80088 38#include "AliITSRecPointU.h"
0b255a22 39#include "AliITStrackU.h"
1d9af2d5 40#include "AliITStrackerMI.h"
41#include "AliITSlayerUpgrade.h"
42#include "AliITSsegmentationUpgrade.h"
43#include "AliITStrackerUpgrade.h"
44#include "AliITSReconstructor.h"
45#include "AliLog.h"
46#include "AliRun.h"
47
48ClassImp(AliITStrackerUpgrade)
49
50//____________________________________________________________________________
556b741c 51 AliITStrackerUpgrade::AliITStrackerUpgrade():
52 AliITStrackerMI(),
53 fNLayers(), fPhiEstimate(0), fITSStandAlone(0), fLambdac(0),
54 fPhic(0), fCoef1(0), fCoef2(0), fCoef3(0), fNloop(0),
55 fPhiWin(0), fLambdaWin(0), fVert(0), fVertexer(0),
56 fListOfTracks(0), fListOfUTracks(0), fITSclusters(0),
57 fInwardFlag(0), fOuterStartLayer(0), fInnerStartLayer(0),
58 fMinNPoints(0), fMinQ(0), fTrackToFollow(), fLayers(0), fSegmentation(0x0),
59 fCluLayer(0), fCluCoord(0)
4211bbc5 60{
61 // Default constructor
62 Init();
665c0bc7 63
4211bbc5 64}
65//____________________________________________________________________________
556b741c 66AliITStrackerUpgrade::AliITStrackerUpgrade(Int_t nLay):
67 AliITStrackerMI(),
68 fNLayers(nLay), fPhiEstimate(0), fITSStandAlone(0), fLambdac(0),
69 fPhic(0), fCoef1(0), fCoef2(0), fCoef3(0), fNloop(0),
70 fPhiWin(0), fLambdaWin(0), fVert(0), fVertexer(0),
71 fListOfTracks(0), fListOfUTracks(0), fITSclusters(0),
72 fInwardFlag(0), fOuterStartLayer(0), fInnerStartLayer(nLay),
73 fMinNPoints(0), fMinQ(0), fTrackToFollow(), fLayers(0), fSegmentation(0x0),
74 fCluLayer(0), fCluCoord(0)
4211bbc5 75{
76 //
77 // constructor
78 //
79 Init();
1d9af2d5 80
4211bbc5 81}
665c0bc7 82
1d9af2d5 83//________________________________________________________________________
556b741c 84AliITStrackerUpgrade::AliITStrackerUpgrade(const AliITStrackerUpgrade& tracker):
85 AliITStrackerMI(),
86 fNLayers(tracker.fNLayers), fPhiEstimate(tracker.fPhiEstimate),
87 fITSStandAlone(tracker.fITSStandAlone), fLambdac(tracker.fLambdac),
88 fPhic(tracker.fPhic), fCoef1(tracker.fCoef1), fCoef2(tracker.fCoef2),
89 fCoef3(tracker.fCoef3), fNloop(tracker.fNloop), fPhiWin(tracker.fPhiWin),
90 fLambdaWin(tracker.fLambdaWin), fVert(tracker.fVert), fVertexer(tracker.fVertexer),
91 fListOfTracks(tracker.fListOfTracks), fListOfUTracks(tracker.fListOfUTracks),
92 fITSclusters(tracker.fITSclusters), fInwardFlag(tracker.fInwardFlag),
93 fOuterStartLayer(tracker.fOuterStartLayer), fInnerStartLayer(tracker.fInnerStartLayer),
94 fMinNPoints(tracker.fMinNPoints), fMinQ(tracker.fMinQ), fTrackToFollow(tracker.fTrackToFollow), fLayers(tracker.fLayers),
95 fSegmentation(tracker.fSegmentation), fCluLayer(tracker.fCluLayer), fCluCoord(tracker.fCluCoord)
4211bbc5 96{
97 //
98 // Copy constructor
99 for(Int_t i=0;i<2;i++){
100 fPoint1[i]=tracker.fPoint1[i];
101 fPoint2[i]=tracker.fPoint2[i];
102 fPoint3[i]=tracker.fPoint3[i];
103 fPointc[i]=tracker.fPointc[i];
104 }
105 if(tracker.fVertexer && tracker.fVert){
106 fVert = new AliESDVertex(*tracker.fVert);
107 }
108 else {
109 fVert = tracker.fVert;
110 }
111 for(Int_t i=0;i<fNLayers;i++){
112 fCluLayer[i] = tracker.fCluLayer[i];
113 fCluCoord[i] = tracker.fCluCoord[i];
114 }
1d9af2d5 115}
116//______________________________________________________________________
117AliITStrackerUpgrade& AliITStrackerUpgrade::operator=(const AliITStrackerUpgrade& source){
118 // Assignment operator.
119 this->~AliITStrackerUpgrade();
120 new(this) AliITStrackerUpgrade(source);
121 return *this;
122
123}
124
125//____________________________________________________________________________
126AliITStrackerUpgrade::~AliITStrackerUpgrade(){
127 // destructor
128 // if fVertexer is not null, the AliESDVertex obj. is owned by this class
129 // and is deleted here
130 if(fVertexer){
131 if(fVert)delete fVert;
132 }
133 fVert = 0;
134 fVertexer = 0;
135
136 if(fPhiWin)delete []fPhiWin;
137 if(fLambdaWin)delete []fLambdaWin;
138 fListOfTracks->Delete();
139 delete fListOfTracks;
0b255a22 140 fListOfUTracks->Delete();
141 delete fListOfUTracks;
1d9af2d5 142 if(fCluLayer){
665c0bc7 143 for(Int_t i=0;i<fNLayers;i++){
1d9af2d5 144 if(fCluLayer[i]){
145 fCluLayer[i]->Delete();
146 delete fCluLayer[i];
147 }
148 }
149 delete [] fCluLayer;
150 }
151 if(fCluCoord){
665c0bc7 152 for(Int_t i=0;i<fNLayers;i++){
1d9af2d5 153 if(fCluCoord[i]){
154 fCluCoord[i]->Delete();
155 delete fCluCoord[i];
156 }
157 }
158 delete [] fCluCoord;
159 }
160
f059e2be 161 if(fSegmentation) delete fSegmentation;
1d9af2d5 162}
163
164//____________________________________________________________________________
165Int_t AliITStrackerUpgrade::Clusters2Tracks(AliESDEvent *event){
166 // This method is used to find and fit the tracks. By default the corresponding
167 // method in the parent class is invoked. In this way a combined tracking
168 // TPC+ITS is performed. If the flag fITSStandAlone is true, the tracking
169 // is done in the ITS only. In the standard reconstruction chain this option
170 // can be set via AliReconstruction::SetOption("ITS","onlyITS")
171 Int_t rc=0;
172 if(!fITSStandAlone){
9b615954 173 rc=AliITStrackerMI::Clusters2Tracks(event); // should not be used !!
174 // Note standalone tracking is implememted in "AliITStrackerU"
1d9af2d5 175 }
176 else {
177 AliDebug(1,"Stand Alone flag set: doing tracking in ITS alone\n");
178 }
179 if(!rc){
180 rc=FindTracks(event,kFALSE);
181 if(AliITSReconstructor::GetRecoParam()->GetSAUseAllClusters()==kTRUE) {
9b615954 182 rc=FindTracks(event,kTRUE); // eventually a 2nd time ??????
1d9af2d5 183 }
184 }
185 return rc;
186}
187//_________________________________________________-
188void AliITStrackerUpgrade::Init(){
189 // Reset all data members
665c0bc7 190 const Int_t nL = fNLayers;
1d9af2d5 191 fPhiEstimate=0;
3e1215c5 192 for(Int_t i=0;i<2;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
1d9af2d5 193 fLambdac=0;
194 fPhic=0;
195 fCoef1=0;
196 fCoef2=0;
197 fCoef3=0;
198 fPointc[0]=0;
199 fPointc[1]=0;
200 fVert = 0;
201 fVertexer = 0;
202 Int_t nLoops=AliITSReconstructor::GetRecoParam()->GetNLoopsSA();
203 if(nLoops==33){
204 SetFixedWindowSizes();
205 }else{
206 Double_t phimin=AliITSReconstructor::GetRecoParam()->GetMinPhiSA();
207 Double_t phimax=AliITSReconstructor::GetRecoParam()->GetMaxPhiSA();
208 Double_t lambmin=AliITSReconstructor::GetRecoParam()->GetMinLambdaSA();
209 Double_t lambmax=AliITSReconstructor::GetRecoParam()->GetMaxLambdaSA();
210 SetCalculatedWindowSizes(nLoops,phimin,phimax,lambmin,lambmax);
211 }
212 fMinQ=AliITSReconstructor::GetRecoParam()->GetSAMinClusterCharge();
213 fITSclusters = 0;
214 SetOuterStartLayer(1);
215 SetSAFlag(kFALSE);
4211bbc5 216 fListOfTracks=new TClonesArray("AliITStrackU",100);
0b255a22 217 fListOfUTracks=new TClonesArray("AliITStrackU",100);
665c0bc7 218 fLayers=new AliITSlayerUpgrade*[nL];//to be fixed
1d9af2d5 219 fCluLayer = 0;
220 fCluCoord = 0;
221 fMinNPoints = 3;
665c0bc7 222 for(Int_t layer=0; layer<fNLayers; layer++){
1d9af2d5 223 Double_t p=0.;
224 Double_t zC= 0.;
225 fLayers[layer] = new AliITSlayerUpgrade(p,zC);
226 }
227
db688ae2 228 fSegmentation = new AliITSsegmentationUpgrade();
1d9af2d5 229
230}
231//_______________________________________________________________________
232Int_t AliITStrackerUpgrade::LoadClusters(TTree *clusTree){
f059e2be 233 //
234 // Load clusters for tracking
235 //
8d499867 236
0ac80088 237 TClonesArray statITSCluster("AliITSRecPointU");
8d499867 238 TClonesArray *ITSCluster = &statITSCluster;
239
f059e2be 240 TBranch* itsClusterBranch=clusTree->GetBranch("ITSRecPoints");
241 if (!itsClusterBranch){
242 AliError("can't get the branch with the ITS clusters ! \n");
243 return 1;
244 }
245 itsClusterBranch->SetAddress(&ITSCluster);
246 clusTree->GetEvent(0);
247 Int_t nCluster = ITSCluster->GetEntriesFast();
248 for(Int_t i=0; i<nCluster; i++){
0ac80088 249 AliITSRecPointU *recp = (AliITSRecPointU*)ITSCluster->UncheckedAt(i);
250 fLayers[recp->GetLayer()]->InsertCluster(new AliITSRecPointU(*recp));
f059e2be 251 }//loop clusters
8d499867 252
1d9af2d5 253 SetClusterTree(clusTree);
254 return 0;
255}
256
257//_______________________________________________________
258void AliITStrackerUpgrade::ResetForFinding(){
259 // Reset data members used in all loops during track finding
260 fPhiEstimate=0;
3e1215c5 261 for(Int_t i=0;i<2;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
1d9af2d5 262 fLambdac=0;
263 fPhic=0;
264 fCoef1=0;
265 fCoef2=0;
266 fCoef3=0;
267 fPointc[0]=0;
268 fPointc[1]=0;
269 fListOfTracks->Clear();
0b255a22 270 fListOfUTracks->Clear();
1d9af2d5 271}
272
1d9af2d5 273//______________________________________________________________________
274Int_t AliITStrackerUpgrade::FindTracks(AliESDEvent* event,Bool_t useAllClusters){
275
276 // Track finder using the ESD object
277 AliDebug(2,Form(" field is %f",event->GetMagneticField()));
9b615954 278
1d9af2d5 279 if(!fITSclusters){
280 Fatal("FindTracks","ITS cluster tree is not accessed - Abort!!!\n Please use method SetClusterTree to pass the pointer to the tree\n");
281 return -1;
282 }
283
9b615954 284 // Reads event and mark clusters of tracks already found, with flag kITSin
1d9af2d5 285 Int_t nentr=event->GetNumberOfTracks();
286 if(!useAllClusters) {
287 while (nentr--) {
288 AliESDtrack *track=event->GetTrack(nentr);
289 if ((track->GetStatus()&AliESDtrack::kITSin) == AliESDtrack::kITSin){
290 Int_t idx[12];
291 Int_t ncl = track->GetITSclusters(idx);
292 for(Int_t k=0;k<ncl;k++){
0ac80088 293 AliITSRecPointU* cll = (AliITSRecPointU*)GetCluster(idx[k]);
1d9af2d5 294 cll->SetBit(kSAflag);
295 }
296 }
297 }
298 }else{
299 while (nentr--) {
300 AliESDtrack *track=event->GetTrack(nentr);
301 if ((track->GetStatus()&AliESDtrack::kITSin) == AliESDtrack::kITSin){
302 Int_t idx[12];
303 Int_t ncl = track->GetITSclusters(idx);
304 for(Int_t k=0;k<ncl;k++){
0ac80088 305 AliITSRecPointU* cll = (AliITSRecPointU*)GetCluster(idx[k]);
1d9af2d5 306 cll->ResetBit(kSAflag);
307 }
308 }
309 }
310 }
311 //Get primary vertex
312 Double_t primaryVertex[3];
313 event->GetVertex()->GetXYZ(primaryVertex);
4211bbc5 314 Int_t nclusters[fNLayers];
665c0bc7 315 for(Int_t i=0; i<fNLayers; i++){
0b255a22 316 nclusters[i]=0;
317 }
318
4211bbc5 319 Int_t dmar[fNLayers];
665c0bc7 320 for(Int_t i=0; i<fNLayers; i++){
0b255a22 321 dmar[i]=0;
322 }
323
1d9af2d5 324 if (fCluLayer == 0) {
665c0bc7 325 fCluLayer = new TClonesArray*[fNLayers];
326 fCluCoord = new TClonesArray*[fNLayers];
327 for(Int_t i=0;i<fNLayers;i++) {
1d9af2d5 328 fCluLayer[i]=0;
329 fCluCoord[i]=0;
330 }
331 }
665c0bc7 332 for(Int_t i=0;i<fNLayers;i++){
1d9af2d5 333 if (!ForceSkippingOfLayer(i)) {
334 for(Int_t cli=0;cli<fLayers[i]->GetNumberOfClusters();cli++){
0ac80088 335 AliITSRecPointU* cls = (AliITSRecPointU*)fLayers[i]->GetCluster(cli);
1d9af2d5 336 if(cls->TestBit(kSAflag)==kTRUE) continue; //clusters used by TPC prol.
337 if(cls->GetQ()==0) continue; //fake clusters dead zones
338 if(i>1 && cls->GetQ()<=fMinQ) continue; // cut on SDD and SSD cluster charge
339 nclusters[i]++;
340 }
341 }
342 dmar[i]=0;
343 if(!fCluLayer[i]){
0ac80088 344 fCluLayer[i] = new TClonesArray("AliITSRecPointU",nclusters[i]);
1d9af2d5 345 }else{
346 fCluLayer[i]->Delete();
347 fCluLayer[i]->Expand(nclusters[i]);
348 }
349 if(!fCluCoord[i]){
350 fCluCoord[i] = new TClonesArray("AliITSclusterTable",nclusters[i]);
351 }else{
352 fCluCoord[i]->Delete();
353 fCluCoord[i]->Expand(nclusters[i]);
354 }
355 }
665c0bc7 356 for(Int_t ilay=0;ilay<fNLayers;ilay++){
1d9af2d5 357 TClonesArray &clulay = *fCluLayer[ilay];
358 TClonesArray &clucoo = *fCluCoord[ilay];
359 if (!ForceSkippingOfLayer(ilay)){
f059e2be 360 AliDebug(2,Form("number of clusters in layer %i : %i",ilay,fLayers[ilay]->GetNumberOfClusters()));
1d9af2d5 361 for(Int_t cli=0;cli<fLayers[ilay]->GetNumberOfClusters();cli++){
0ac80088 362 AliITSRecPointU* cls = (AliITSRecPointU*)fLayers[ilay]->GetCluster(cli);
1d9af2d5 363 if(cls->TestBit(kSAflag)==kTRUE) continue;
364 if(cls->GetQ()==0) continue;
365 Double_t phi=0;Double_t lambda=0;
366 Double_t x=0;Double_t y=0;Double_t z=0;
367 Float_t sx=0;Float_t sy=0;Float_t sz=0;
368 GetCoorAngles(cls,phi,lambda,x,y,z,primaryVertex);
369 GetCoorErrors(cls,sx,sy,sz);
370 new (clulay[dmar[ilay]]) AliITSRecPoint(*cls);
371 new (clucoo[dmar[ilay]]) AliITSclusterTable(x,y,z,sx,sy,sz,phi,lambda,cli);
372 dmar[ilay]++;
373 }
db688ae2 374 } else AliDebug(2,Form("Force skipping layer %i",ilay));
1d9af2d5 375 }
376
377
378 // track counter
379 Int_t ntrack=0;
380
4211bbc5 381 static Int_t nClusLay[fgMaxNLayer];//counter for clusters on each layer
1d9af2d5 382 Int_t startLayForSeed=0;
383 Int_t lastLayForSeed=fOuterStartLayer;
384 Int_t nSeedSteps=lastLayForSeed-startLayForSeed;
385 Int_t seedStep=1;
386 if(fInwardFlag){
665c0bc7 387 startLayForSeed=fNLayers-1;
1d9af2d5 388 lastLayForSeed=fInnerStartLayer;
389 nSeedSteps=startLayForSeed-lastLayForSeed;
390 seedStep=-1;
391 }
392 // loop on minimum number of points
665c0bc7 393 for(Int_t iMinNPoints=fNLayers; iMinNPoints>=fMinNPoints; iMinNPoints--) {
1d9af2d5 394 // loop on starting layer for track finding
395 for(Int_t iSeedLay=0; iSeedLay<=nSeedSteps; iSeedLay++) {
396 Int_t theLay=startLayForSeed+iSeedLay*seedStep;
397 if(ForceSkippingOfLayer(theLay)) continue;
398 Int_t minNPoints=iMinNPoints-theLay;
665c0bc7 399 if(fInwardFlag) minNPoints=iMinNPoints-(fNLayers-1-theLay);
400 for(Int_t i=theLay+1;i<fNLayers;i++)
1d9af2d5 401 if(ForceSkippingOfLayer(i))
402 minNPoints--;
403 if(minNPoints<fMinNPoints) continue;
404
405 // loop on phi and lambda window size
406 for(Int_t nloop=0;nloop<fNloop;nloop++){
eaaed734 407 Int_t nclTheLay=fCluLayer[theLay]->GetEntriesFast();
1d9af2d5 408 while(nclTheLay--){
409 ResetForFinding();
410 Bool_t useRP=SetFirstPoint(theLay,nclTheLay,primaryVertex);
411 if(!useRP) continue;
665c0bc7 412 AliITStrackU trs(fNLayers);
1d9af2d5 413
414 Int_t pflag=0;
415 Int_t kk;
665c0bc7 416 for(kk=0;kk<fNLayers;kk++) nClusLay[kk] = 0;
1d9af2d5 417 kk=0;
418 nClusLay[kk] = SearchClusters(theLay,fPhiWin[nloop],fLambdaWin[nloop],
419 &trs,primaryVertex[2],pflag);
420 Int_t nextLay=theLay+seedStep;
dd290e92 421 if(nextLay<0){
422 AliWarning("Negative layer id. Please check!");
423 continue;
424 }
1d9af2d5 425 Bool_t goon=kTRUE;
426 while(goon){
427
428 kk++;
429 nClusLay[kk] = SearchClusters(nextLay,fPhiWin[nloop],fLambdaWin[nloop],
430 &trs,primaryVertex[2],pflag);
431 if(nClusLay[kk]!=0){
432
433 pflag=1;
434 if(kk==1) {
435 fPoint3[0]=fPointc[0];
436 fPoint3[1]=fPointc[1];
437 } else {
438 UpdatePoints();
439 }
440 }
441 nextLay+=seedStep;
665c0bc7 442 if(nextLay<0 || nextLay==fNLayers) goon=kFALSE;
1d9af2d5 443 }
444
445
446 Int_t layOK=0;
447 if(!fInwardFlag){
665c0bc7 448 for(Int_t nnp=0;nnp<fNLayers-theLay;nnp++){
1d9af2d5 449 if(nClusLay[nnp]!=0) layOK+=1;
450 }
451 }else{
452 for(Int_t nnp=theLay; nnp>=0; nnp--){
453 if(nClusLay[nnp]!=0) layOK+=1;
454 }
455 }
456 if(layOK>=minNPoints){
9b615954 457 // AliDebug(2,Form("---NPOINTS: %d; MAP: %d %d %d %d %d %d %d\n",layOK,nClusLay[0],nClusLay[1],nClusLay[2],nClusLay[3],nClusLay[4],nClusLay[5],nClusLay[6]));
1d9af2d5 458 AliITStrackV2* tr2 = 0;
459 Bool_t onePoint = kFALSE;
4211bbc5 460 tr2 = FitTrack(&trs,primaryVertex,onePoint);
461 if(!tr2){
462 continue;
463 }
9b615954 464 // AliDebug(2,Form("---NPOINTS fit: %d\n",tr2->GetNumberOfClusters()));
1d9af2d5 465
4211bbc5 466 StoreTrack(tr2,event,useAllClusters);
467 ntrack++;
468
469 }
1d9af2d5 470
471 }//end loop on clusters of theLay
472 } //end loop on window sizes
473 } //end loop on theLay
474 }//end loop on min points
475 if(!useAllClusters) AliInfo(Form("Number of found tracks: %d",event->GetNumberOfTracks()));
476
477 ResetForFinding();
478 return 0;
479
480}
481
482//________________________________________________________________________
483
0b255a22 484AliITStrackV2* AliITStrackerUpgrade::FitTrack(AliITStrackU* tr,Double_t *primaryVertex,Bool_t onePoint) {
1d9af2d5 485
0b255a22 486 const Int_t kMaxClu=AliITStrackU::kMaxNumberOfClusters;
1d9af2d5 487
4211bbc5 488 static Int_t clind[fgMaxNLayer][kMaxClu];
489 static Int_t clmark[fgMaxNLayer][kMaxClu];
490 static Int_t end[fgMaxNLayer];
0ac80088 491 static AliITSRecPointU *listlayer[fgMaxNLayer][kMaxClu];
556b741c 492 static Int_t indices[fgMaxNLayer];
493
4211bbc5 494 for(Int_t k=0;k<fgMaxNLayer; k++)end[k]=0;
1d9af2d5 495
665c0bc7 496 for(Int_t i=0;i<fNLayers;i++) {
1d9af2d5 497 for(Int_t j=0;j<kMaxClu; j++){
498 clind[i][j]=0;
499 clmark[i][j]=0;
500 listlayer[i][j]=0;
501 }
502 }
665c0bc7 503 Int_t inx[fNLayers];
504 for (Int_t k=0; k<fNLayers; k++) inx[k]=-1;
0b255a22 505 Int_t nclusters = tr->GetNumberOfClustersU();
1d9af2d5 506 for(Int_t ncl=0;ncl<nclusters;ncl++){
0b255a22 507 Int_t index = tr->GetClusterIndexU(ncl);
1d9af2d5 508 Int_t lay = (index & 0xf0000000) >> 28;
7dc6449a 509 //Int_t cli = index&0x0fffffff;
9e396337 510 //AliITSRecPoint* cl = (AliITSRecPoint*)fLayers[lay]->GetCluster(cli);
0ac80088 511 AliITSRecPointU* cl = (AliITSRecPointU*)GetCluster(index);
1d9af2d5 512 if(cl->TestBit(kSAflag)==kTRUE) cl->ResetBit(kSAflag);
513 inx[lay]=index;
514
515 Int_t nInLay=end[lay];
516 listlayer[lay][nInLay]=cl;
517 clind[lay][nInLay]=index;
518 end[lay]++;
519 }
520
665c0bc7 521 for(Int_t nlay=0;nlay<fNLayers;nlay++){
1d9af2d5 522 for(Int_t ncl=0;ncl<tr->GetNumberOfMarked(nlay);ncl++){
523 Int_t mark = tr->GetClusterMark(nlay,ncl);
524 clmark[nlay][ncl]=mark;
525 }
526 }
527
528
529 Int_t firstLay=-1,secondLay=-1;
665c0bc7 530 for(Int_t i=0;i<fNLayers;i++) {
1d9af2d5 531 if(end[i]==0) {
532 end[i]=1;
533 }else{
534 if(firstLay==-1) {
535 firstLay=i;
536 } else if(secondLay==-1) {
537 secondLay=i;
538 }
539 }
540 }
541
4211bbc5 542 //---adding initialization for non existent layers
543 for(Int_t i=fNLayers;i<fgMaxNLayer;i++) {
544 if(end[i]==0) {
545 end[i]=1;
546 }
547 }
548
549
550
1d9af2d5 551 if(firstLay==-1 || (secondLay==-1 && !onePoint)) return 0;
4211bbc5 552
0b255a22 553 TClonesArray &arrSA= *fListOfUTracks;
1d9af2d5 554 Int_t nFoundTracks=0;
555
4211bbc5 556
1d9af2d5 557 for(Int_t l0=0;l0<end[0];l0++){ //loop on layer 1
556b741c 558 indices[0]=l0;
1d9af2d5 559 for(Int_t l1=0;l1<end[1];l1++){ //loop on layer 2
556b741c 560 indices[1]=l1;
1d9af2d5 561 for(Int_t l2=0;l2<end[2];l2++){ //loop on layer 3
556b741c 562 indices[2]=l2;
1d9af2d5 563 for(Int_t l3=0;l3<end[3];l3++){ //loop on layer 4
556b741c 564 indices[3]=l3;
1d9af2d5 565 for(Int_t l4=0;l4<end[4];l4++){ //loop on layer 5
556b741c 566 indices[4]=l4;
1d9af2d5 567 for(Int_t l5=0;l5<end[5];l5++){ //loop on layer 6
556b741c 568 indices[5]=l5;
4211bbc5 569 for(Int_t l6=0;l6<end[6];l6++){//loop on layer 7
556b741c 570 indices[6]=l6;
4211bbc5 571 for(Int_t l7=0;l7<end[7];l7++){ //loop on layer 8
556b741c 572 indices[7]=l7;
573
574 Int_t iFirstLay=indices[firstLay];
575 Int_t mrk1=clmark[firstLay][iFirstLay];
576
4211bbc5 577 Double_t x1,y1,z1,sx1,sy1,sz1;
578 Double_t x2,y2,z2,sx2,sy2,sz2;
556b741c 579 AliITSRecPointU* p1=(AliITSRecPointU*)listlayer[firstLay][iFirstLay];
4211bbc5 580 Int_t layer = 0;
581 Double_t radius = 0.;
582 Double_t yclu1 = p1->GetY();
583 Double_t zclu1 = p1->GetZ();
584 layer=p1->GetLayer();
585 radius = fSegmentation->GetRadius(layer);
586 Double_t cv=0.,tgl2=0.,phi2=0.;
556b741c 587 AliITSclusterTable* arr1 = (AliITSclusterTable*)GetClusterCoord(firstLay,mrk1);
4211bbc5 588 x1 = arr1->GetX();
589 y1 = arr1->GetY();
590 z1 = arr1->GetZ();
591 sx1 = arr1->GetSx();
592 sy1 = arr1->GetSy();
593 sz1 = arr1->GetSz();
594
595 if(secondLay>0) {
556b741c 596 Int_t iSecondLay=indices[secondLay];
597 Int_t mrk2=clmark[secondLay][iSecondLay];
598 AliITSclusterTable* arr2 = (AliITSclusterTable*)GetClusterCoord(secondLay,mrk2);
4211bbc5 599 x2 = arr2->GetX();
600 y2 = arr2->GetY();
601 z2 = arr2->GetZ();
602 sx2 = arr2->GetSx();
603 sy2 = arr2->GetSy();
604 sz2 = arr2->GetSz();
605 cv = Curvature(primaryVertex[0],primaryVertex[1],x1,y1,x2,y2);
606 tgl2 = (z2-z1)/TMath::Sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
607 phi2 = TMath::ATan2((y2-y1),(x2-x1));
608 }
556b741c 609
4211bbc5 610 Double_t x,y,z = 0.;
611 Double_t xz[2];
556b741c 612 xz[0]= p1->GetDetLocalX();
613 xz[1]= p1->GetDetLocalZ();
614 fSegmentation->DetToGlobal(layer,p1->GetModule(),xz[0], xz[1],x,y,z);
4211bbc5 615
616
0ac80088 617 Double_t phiclrad;
4211bbc5 618 phiclrad=TMath::ATan2(y,x);
619 if (phiclrad<0) phiclrad+=TMath::TwoPi();
620 else if (phiclrad>=TMath::TwoPi()) phiclrad-=TMath::TwoPi();
621
556b741c 622 Double_t alpha = fSegmentation->GetAlpha(fSegmentation->GetModule(phiclrad));
317075a3 623
624 Int_t clIndice[fNLayers];
625 for (Int_t k=0; k<fNLayers; k++) clIndice[k]=-1;
0ac80088 626
4211bbc5 627 AliITStrackU trac(alpha,radius,yclu1,zclu1,phi2,tgl2,cv,1,fNLayers);
317075a3 628 for(Int_t iLay=fNLayers-1; iLay>=0; iLay--){
556b741c 629 Int_t iInLay=indices[iLay];
630 AliITSRecPointU* cl=(AliITSRecPointU*)listlayer[iLay][iInLay];
631 if(cl!=0){
632 trac.AddClusterV2(iLay,(clind[iLay][iInLay] & 0x0fffffff)>>0);
633 trac.AddClusterMark(iLay,clmark[iLay][iInLay]);
317075a3 634 clIndice[iLay] = clind[iLay][iInLay];
556b741c 635 }
4211bbc5 636 }
556b741c 637
4211bbc5 638 //fit with Kalman filter using RefitAtBase()
639 AliITStrackU ot(trac, kFALSE);
640 ot.ResetCovariance(10.);
641 ot.ResetClusters();
1d9af2d5 642
4211bbc5 643 // Propagate inside the innermost layer with a cluster
644 if(ot.Propagate(ot.GetX()-0.1*ot.GetX())) {
317075a3 645 if(RefitAtBase(AliITSRecoParam::GetrInsideITSscreen(),&ot,clIndice)){ //fit from layer 1 to layer N
4211bbc5 646 AliITStrackU otrack2(ot,kFALSE);
647 otrack2.ResetCovariance(10.);
648 otrack2.ResetClusters();
4211bbc5 649 //fit from last layer to layer 1
317075a3 650 if(RefitAtBase(fSegmentation->GetRadius(0)-0.1,&otrack2,clIndice)) {
4211bbc5 651 new(arrSA[nFoundTracks]) AliITStrackU(otrack2);
652 ++nFoundTracks;
653 }
1d9af2d5 654
4211bbc5 655 }
656 }
657
4211bbc5 658 }//end loop layer 8
659 }//end loop layer 7
1d9af2d5 660 }//end loop layer 6
661 }//end loop layer 5
662 }//end loop layer 4
663 }//end loop layer 3
664 }//end loop layer 2
665 }//end loop layer 1
666
667
668
669
4211bbc5 670 if(fListOfUTracks->GetEntries()==0) return 0;
1d9af2d5 671
672 Int_t lowchi2 = FindTrackLowChiSquare();
4211bbc5 673 AliITStrackU* otrack = (AliITStrackU*)fListOfUTracks->At(lowchi2);
674
1d9af2d5 675 if(otrack==0) {
676 return 0;
677 }
0a64515a 678
9b615954 679 // CookLabel(otrack,0.); //MI change - to see fake ratio
680 // printf(" label from CookLabel: %d \n ",otrack->GetLabel());
681
0a64515a 682 Int_t label=FindLabel(otrack);
683 otrack->SetLabel(label);
9b615954 684 // printf(" label from FindLabel: %d \n ",otrack->GetLabel());
0a64515a 685
665c0bc7 686 Int_t indexc[fNLayers];
687 for(Int_t i=0;i<fNLayers;i++) indexc[i]=0;
1d9af2d5 688 for(Int_t nind=0;nind<otrack->GetNumberOfClusters();nind++){
689 indexc[nind] = otrack->GetClusterIndex(nind);
690 }
9b615954 691
692
1d9af2d5 693 //remove clusters of found track
665c0bc7 694 for(Int_t nlay=0;nlay<fNLayers;nlay++){
4211bbc5 695 for(Int_t cln=0;cln<otrack->GetNumberOfMarked(nlay);cln++){
696 Int_t index = otrack->GetClusterMark(nlay,cln);
1d9af2d5 697 fCluLayer[nlay]->RemoveAt(index);
698 RemoveClusterCoord(nlay,index);
699 fCluLayer[nlay]->Compress();
700 }
701 }
9b615954 702
1d9af2d5 703 return otrack;
704
705}
706
707//_______________________________________________________
708void AliITStrackerUpgrade::StoreTrack(AliITStrackV2 *t,AliESDEvent *event, Bool_t pureSA) /*const*/
709{
710 //
711 // Add new track to the ESD
712 //
713 AliESDtrack outtrack;
714 outtrack.UpdateTrackParams(t,AliESDtrack::kITSin);
715 if(pureSA) outtrack.SetStatus(AliESDtrack::kITSpureSA);
716 outtrack.SetStatus(AliESDtrack::kITSout);
717 outtrack.SetStatus(AliESDtrack::kITSrefit);
718 Double_t sdedx[4]={0.,0.,0.,0.};
719 for(Int_t i=0; i<4; i++) sdedx[i]=t->GetSampledEdx(i);
720 outtrack.SetITSdEdxSamples(sdedx);
721 event->AddTrack(&outtrack);
722 return;
723}
724
725
726//_______________________________________________________
0b255a22 727Int_t AliITStrackerUpgrade::SearchClusters(Int_t layer,Double_t phiwindow,Double_t lambdawindow, AliITStrackU* trs,Double_t /*zvertex*/,Int_t pflag){
1d9af2d5 728 //function used to to find the clusters associated to the track
729
8d499867 730 AliDebug(2,"Starting...");
731 if(ForceSkippingOfLayer(layer)) {
f059e2be 732 AliDebug(2,Form("Forcing skipping of layer %i. Exiting",layer));
733 return 0;
8d499867 734 }
1d9af2d5 735
736 Int_t nc=0;
db688ae2 737 Double_t r=fSegmentation->GetRadius(layer);
1d9af2d5 738 if(pflag==1){
739 Float_t cx1,cx2,cy1,cy2;
740 FindEquation(fPoint1[0],fPoint1[1],fPoint2[0],fPoint2[1],fPoint3[0],fPoint3[1],fCoef1,fCoef2,fCoef3);
741 if (FindIntersection(fCoef1,fCoef2,fCoef3,-r*r,cx1,cy1,cx2,cy2)==0)
742 return 0;
743 Double_t fi1=TMath::ATan2(cy1-fPoint1[1],cx1-fPoint1[0]);
744 Double_t fi2=TMath::ATan2(cy2-fPoint1[1],cx2-fPoint1[0]);
745 fPhiEstimate=ChoosePoint(fi1,fi2,fPhic);
746 }
747
eaaed734 748
749 Double_t phiExpect=fPhiEstimate;
750 Double_t lamExpect=fLambdac;
751
752 Int_t ncl = fCluLayer[layer]->GetEntriesFast();
1d9af2d5 753 for (Int_t index=0; index<ncl; index++) {
eaaed734 754
1d9af2d5 755 AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(layer,index);
1d9af2d5 756
757 Double_t lambda = arr->GetLambda();
eaaed734 758 if (TMath::Abs(lambda-lamExpect)>lambdawindow) continue;
759
760 Double_t phi = arr->GetPhi();
761 Double_t deltaPhi = phi-phiExpect;
762 if(deltaPhi>TMath::Pi()) deltaPhi-=2*TMath::Pi();
763 else if(deltaPhi<-TMath::Pi()) deltaPhi+=2*TMath::Pi();
764 if (TMath::Abs(deltaPhi)>phiwindow) continue;
1d9af2d5 765
0b255a22 766 if(trs->GetNumberOfClustersU()==trs->GetMaxNumberOfClusters()) return 0;
1d9af2d5 767 if(trs->GetNumberOfMarked(layer)==trs->GetMaxNMarkedPerLayer()) return 0;
768 Int_t orind = arr->GetOrInd();
0b255a22 769 trs->AddClusterU(layer,orind);
1d9af2d5 770 trs->AddClusterMark(layer,index);
1d9af2d5 771 nc++;
772 fLambdac=lambda;
773 fPhiEstimate=phi;
774
775 fPointc[0]=arr->GetX();
776 fPointc[1]=arr->GetY();
eaaed734 777
1d9af2d5 778 }
779 return nc;
eaaed734 780
1d9af2d5 781}
782
783//________________________________________________________________
784Bool_t AliITStrackerUpgrade::SetFirstPoint(Int_t lay, Int_t clu, Double_t* primaryVertex){
785 // Sets the first point (seed) for tracking
786
0ac80088 787 AliITSRecPointU* cl = (AliITSRecPointU*)fCluLayer[lay]->At(clu);
1d9af2d5 788 if(!cl) return kFALSE;
789 if (cl->GetQ()<=0) return kFALSE;
790 if(lay>1 && cl->GetQ()<=fMinQ) return kFALSE;
791 AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(lay,clu);
792 fPhic = arr->GetPhi();
793 fLambdac = arr->GetLambda();
794 fPhiEstimate = fPhic;
795 fPoint1[0]=primaryVertex[0];
796 fPoint1[1]=primaryVertex[1];
797 fPoint2[0]=arr->GetX();
798 fPoint2[1]=arr->GetY();
799 return kTRUE;
800}
801
802//________________________________________________________________
803void AliITStrackerUpgrade::UpdatePoints(){
804 //update of points for the estimation of the curvature
805
806 fPoint2[0]=fPoint3[0];
807 fPoint2[1]=fPoint3[1];
808 fPoint3[0]=fPointc[0];
809 fPoint3[1]=fPointc[1];
1d9af2d5 810}
1d9af2d5 811//___________________________________________________________________
812Int_t AliITStrackerUpgrade::FindEquation(Float_t x1, Float_t y1, Float_t x2, Float_t y2, Float_t x3, Float_t y3,Float_t& a, Float_t& b, Float_t& c){
813
814 //given (x,y) of three recpoints (in global coordinates)
815 //returns the parameters a,b,c of circonference x*x + y*y +a*x + b*y +c
816
817 Float_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
818 if(den==0) return 0;
819 a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
820 b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
821 c = -x1*x1-y1*y1-a*x1-b*y1;
822 return 1;
823}
824//__________________________________________________________________________
825Int_t AliITStrackerUpgrade::FindIntersection(Float_t a1, Float_t b1, Float_t c1, Float_t c2,Float_t& x1,Float_t& y1, Float_t& x2, Float_t& y2){
826
827 //Finds the intersection between the circonference of the track and the circonference centered in (0,0) represented by one layer
828 //c2 is -rlayer*rlayer
829
830 if(a1==0) return 0;
831 Double_t m = c2-c1;
832 Double_t aA = (b1*b1)/(a1*a1)+1;
833 Double_t bB = (-2*m*b1/(a1*a1));
834 Double_t cC = c2+(m*m)/(a1*a1);
835 Double_t dD = bB*bB-4*aA*cC;
836 if(dD<0) return 0;
837
838 y1 = (-bB+TMath::Sqrt(dD))/(2*aA);
839 y2 = (-bB-TMath::Sqrt(dD))/(2*aA);
840 x1 = (c2-c1-b1*y1)/a1;
841 x2 = (c2-c1-b1*y2)/a1;
842
843 return 1;
844}
845//____________________________________________________________________
846Double_t AliITStrackerUpgrade::Curvature(Double_t x1,Double_t y1,Double_t
847 x2,Double_t y2,Double_t x3,Double_t y3){
848
849 //calculates the curvature of track
850 Double_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
851 if(den==0) return 0;
852 Double_t a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
853 Double_t b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
854 Double_t c = -x1*x1-y1*y1-a*x1-b*y1;
855 Double_t xc=-a/2.;
856
857 if((a*a+b*b-4*c)<0) return 0;
858 Double_t rad = TMath::Sqrt(a*a+b*b-4*c)/2.;
859 if(rad==0) return 0;
860
861 if((x1>0 && y1>0 && x1<xc)) rad*=-1;
862 if((x1<0 && y1>0 && x1<xc)) rad*=-1;
863 // if((x1<0 && y1<0 && x1<xc)) rad*=-1;
864 // if((x1>0 && y1<0 && x1<xc)) rad*=-1;
865
866 return 1/rad;
867
868}
1d9af2d5 869//____________________________________________________________________
870Double_t AliITStrackerUpgrade::ChoosePoint(Double_t p1, Double_t p2, Double_t pp){
871
872 //Returns the point closest to pp
873
874 Double_t diff1 = p1-pp;
875 Double_t diff2 = p2-pp;
876
877 if(TMath::Abs(diff1)<TMath::Abs(diff2)) fPhiEstimate=p1;
878 else fPhiEstimate=p2;
879 return fPhiEstimate;
880
881}
1d9af2d5 882//_________________________________________________________________
883Int_t AliITStrackerUpgrade::FindTrackLowChiSquare() const {
884 // returns track with lowest chi square
4211bbc5 885 Int_t dim=fListOfUTracks->GetEntries();
1d9af2d5 886 if(dim<=1) return 0;
4211bbc5 887 AliITStrackV2* trk = (AliITStrackV2*)fListOfUTracks->At(0);
1d9af2d5 888 Double_t minChi2=trk->GetChi2();
889 Int_t index=0;
890 for(Int_t i=1;i<dim;i++){
4211bbc5 891 trk = (AliITStrackV2*)fListOfUTracks->At(i);
1d9af2d5 892 Double_t chi2=trk->GetChi2();
893 if(chi2<minChi2){
894 minChi2=chi2;
895 index=i;
896 }
897 }
898 return index;
899}
900
901//__________________________________________________________
0a64515a 902Int_t AliITStrackerUpgrade::FindLabel(AliITStrackV2* track) const {
903 //
9b615954 904 // max number of labels in AliITSRecPointU is 24;
905 enum {kMaxLab=24};
1d9af2d5 906
9b615954 907 Int_t labl[fNLayers][kMaxLab];
908 Int_t cnts[fNLayers][kMaxLab];
665c0bc7 909 for(Int_t j=0;j<fNLayers;j++){
9b615954 910 for(Int_t k=0;k<kMaxLab;k++){
0a64515a 911 labl[j][k]=-2;
912 cnts[j][k]=1;
913 }
914 }
915 Int_t iNotLabel=0;
916 for(Int_t i=0;i<track->GetNumberOfClusters(); i++) {
917 Int_t indexc = track->GetClusterIndex(i);
0ac80088 918 AliITSRecPointU* cl = (AliITSRecPointU*)GetCluster(indexc);
9b615954 919 // AliDebug(2,Form("Layer %d; cluster index %i; MClabels=(%d,%d,%d) ",
920 // cl->GetLayer(),indexc,cl->GetLabel(0),cl->GetLabel(1),cl->GetLabel(2)));
0a64515a 921 Int_t iLayer=cl->GetLayer();
9b615954 922 for(Int_t k=0;k<kMaxLab;k++){
923 labl[iLayer][k]=cl->GetTrackID(k);
0a64515a 924 if(labl[iLayer][k]<0) iNotLabel++;
925 }
926 }
9b615954 927 if(iNotLabel==kMaxLab*track->GetNumberOfClusters()) return -2; // No MC label found at all
0a64515a 928
665c0bc7 929 for(Int_t j1=0;j1<fNLayers; j1++) {
0a64515a 930 for(Int_t j2=0; j2<j1; j2++){
9b615954 931 for(Int_t k1=0; k1<kMaxLab; k1++){
932 for(Int_t k2=0; k2<kMaxLab; k2++){
0a64515a 933 if(labl[j1][k1]>=0 && labl[j1][k1]==labl[j2][k2] && cnts[j2][k2]>0){
934 cnts[j2][k2]++;
935 cnts[j1][k1]=0;
936 }
937 }
938 }
939 }
940 }
1d9af2d5 941
0a64515a 942
943 Int_t cntMax=0;
944 Int_t label=-1;
665c0bc7 945 for(Int_t j=0;j<fNLayers;j++){
9b615954 946 for(Int_t k=0;k<kMaxLab;k++){
0a64515a 947 if(cnts[j][k]>cntMax && labl[j][k]>=0){
948 cntMax=cnts[j][k];
949 label=labl[j][k];
950 }
951 }
952 }
953
954 Int_t lflag=0;
9b615954 955 for(Int_t i=0;i<fNLayers;i++) {
956 Bool_t labelOkInLayer = kFALSE;
957 for(Int_t k=0;k<kMaxLab;k++){
958 if (labl[i][k]==label) {
959 labelOkInLayer=kTRUE;
960 break;
961 }
962 }
963 if(labelOkInLayer) lflag++;
964 }
965 if(lflag<track->GetNumberOfClusters()) {
966 if (label==0) {
967 label = -1; // to avoid label 0 which can not be identified as as fake with -0
968 } else {
969 label = -label;
970 }
971 }
972 if (label<0) {
973 for(Int_t i=0;i<track->GetNumberOfClusters(); i++) {
974 Int_t indexc = track->GetClusterIndex(i);
975 AliITSRecPointU* cl = (AliITSRecPointU*)GetCluster(indexc);
976 AliDebug(2,Form("Layer %d; cluster index %i; MClabels=(%d,%d,%d...) ",
977 cl->GetLayer(),indexc,cl->GetTrackID(0),cl->GetTrackID(1),cl->GetTrackID(2)));
978 }
979 }
980
981
0a64515a 982 return label;
1d9af2d5 983}
9b615954 984
1d9af2d5 985//_____________________________________________________________________________
1d9af2d5 986void AliITStrackerUpgrade::SetCalculatedWindowSizes(Int_t n, Float_t phimin, Float_t phimax, Float_t lambdamin, Float_t lambdamax){
987 // Set sizes of the phi and lambda windows used for track finding
988 fNloop = n;
989 if(fPhiWin) delete [] fPhiWin;
990 if(fLambdaWin) delete [] fLambdaWin;
991 fPhiWin = new Double_t[fNloop];
992 fLambdaWin = new Double_t[fNloop];
993 Float_t stepPhi=(phimax-phimin)/(Float_t)(fNloop-1);
994 Float_t stepLambda=(lambdamax-lambdamin)/(Float_t)(fNloop-1);
995 for(Int_t k=0;k<fNloop;k++){
996 Float_t phi=phimin+k*stepPhi;
997 Float_t lam=lambdamin+k*stepLambda;
998 fPhiWin[k]=phi;
999 fLambdaWin[k]=lam;
1000 }
1001}
1002//_____________________________________________________________________________
1003void AliITStrackerUpgrade::SetFixedWindowSizes(Int_t n, Double_t *phi, Double_t *lam){
1004 // Set sizes of the phi and lambda windows used for track finding
1005 fNloop = n;
1006 if(phi){ // user defined values
1007 fPhiWin = new Double_t[fNloop];
1008 fLambdaWin = new Double_t[fNloop];
1009 for(Int_t k=0;k<fNloop;k++){
1010 fPhiWin[k]=phi[k];
1011 fLambdaWin[k]=lam[k];
1012 }
1013 }
1014 else { // default values
1015
1016 Double_t phid[33] = {0.002,0.003,0.004,0.0045,0.0047,
1017 0.005,0.0053,0.0055,
1018 0.006,0.0063,0.0065,0.007,0.0073,0.0075,0.0077,
1019 0.008,0.0083,0.0085,0.0087,0.009,0.0095,0.0097,
1020 0.01,0.0105,0.011,0.0115,0.012,0.0125,0.013,0.0135,0.0140,0.0145};
1021 Double_t lambdad[33] = {0.003,0.004,0.005,0.005,0.005,
1022 0.005,0.005,0.006,
1023 0.006,0.006,0.006,0.007,0.007,0.007,0.007,
1024 0.007,0.007,0.007,0.007,0.007,0.007,0.007,
1025 0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008};
1026
1027 if(fNloop!=33){
1028 fNloop = 33;
1029 }
1030
1031
1032 fPhiWin = new Double_t[fNloop];
1033 fLambdaWin = new Double_t[fNloop];
1034
1035 Double_t factor=AliITSReconstructor::GetRecoParam()->GetFactorSAWindowSizes();
1036
1037 for(Int_t k=0;k<fNloop;k++){
1038 fPhiWin[k]=phid[k]*factor;
1039 fLambdaWin[k]=lambdad[k]*factor;
1040 }
1041
1042 }
1043
1044}
1045//_______________________________________________________________________
0ac80088 1046void AliITStrackerUpgrade::GetCoorAngles(AliITSRecPointU* cl,Double_t &phi,Double_t &lambda, Double_t &x, Double_t &y,Double_t &z,Double_t* vertex){
1d9af2d5 1047 //Returns values of phi (azimuthal) and lambda angles for a given cluster
1d9af2d5 1048
1049 Double_t xz[2];
1050 xz[0]= cl->GetDetLocalX();
1051 xz[1]= cl->GetDetLocalZ() ;
1d9af2d5 1052 Int_t ilayer;
1053 ilayer = cl->GetLayer();
556b741c 1054 fSegmentation->DetToGlobal(ilayer,cl->GetModule(),xz[0], xz[1],x,y,z);
1d9af2d5 1055
1056 if(x!=0 && y!=0)
1057 phi=TMath::ATan2(y-vertex[1],x-vertex[0]);
1058 lambda=TMath::ATan2(z-vertex[2],TMath::Sqrt((x-vertex[0])*(x-vertex[0])+(y-vertex[1])*(y-vertex[1])));
1059}
1060
1061//________________________________________________________________________
0ac80088 1062void AliITStrackerUpgrade::GetCoorErrors(AliITSRecPointU* cl,Float_t &sx,Float_t &sy, Float_t &sz){
1d9af2d5 1063
1064 //returns sigmax, y, z of cluster in global coordinates
1065 /*
1066 Double_t rot[9]; fGeom->GetRotMatrix(module,rot);
1067 Int_t lay,lad,det;
1068 AliITSgeomTGeo::GetModuleId(module,lay,lad,det);
1069
1070 Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1071 Double_t phi=TMath::Pi()/2+alpha;
1072 if (lay==1) phi+=TMath::Pi();
1073
1074 Float_t cp=TMath::Cos(phi), sp=TMath::Sin(phi);
1075 */
1076 Float_t covm[6];
1077 cl->GetGlobalCov(covm);
1078 sx=TMath::Sqrt(covm[0]);
1079 sy=TMath::Sqrt(covm[3]);
1080 sz=TMath::Sqrt(covm[5]);
1081
1082
1083
1084 // Float_t covm[6];
1085 // cl->GetGlobalCov(covm);
1086 // sx=12*1e-04;//TMath::Sqrt(covm[0]);
1087 // sy=12*1e-04;//TMath::Sqrt(covm[3]);
1088 // sz=40*1e-04;//TMath::Sqrt(covm[5]);
1089 /*
1090 sx = TMath::Sqrt(sp*sp*cl->GetSigmaY2());
1091 sy = TMath::Sqrt(cp*cp*cl->GetSigmaY2());
1092 sz = TMath::Sqrt(cl->GetSigmaZ2());
1093 */
1094}
1095//_____________________________________________________________
1096void AliITStrackerUpgrade::UnloadClusters() {
1097 //--------------------------------------------------------------------
1098 //This function unloads ITS clusters
1099 //--------------------------------------------------------------------
665c0bc7 1100 for (Int_t i=0; i<fNLayers; i++) fLayers[i]->ResetClusters();
1d9af2d5 1101}
1102//______________________________________________________________________
4211bbc5 1103Bool_t AliITStrackerUpgrade::RefitAtBase(Double_t xx,AliITStrackU *track,
1d9af2d5 1104 const Int_t *clusters)
1105{
1106 //--------------------------------------------------------------------
1107 // Simplified version for ITS upgrade studies -- does not use module info
1108 //--------------------------------------------------------------------
665c0bc7 1109 Int_t index[fNLayers];
1d9af2d5 1110 Int_t k;
665c0bc7 1111 for (k=0; k<fNLayers; k++) index[k]=-1;
1d9af2d5 1112 //
665c0bc7 1113 for (k=0; k<fNLayers; k++) {
1d9af2d5 1114 index[k]=clusters[k];
1115 }
1d9af2d5 1116
1117 ULong_t trStatus=0;
1118 if(track->GetESDtrack()) trStatus=track->GetStatus();
1119 Int_t innermostlayer=0;
1120 if(trStatus&AliESDtrack::kTPCin) {
665c0bc7 1121 innermostlayer=fNLayers-1;
1d9af2d5 1122 Double_t drphi = TMath::Abs(track->GetD(0.,0.));
665c0bc7 1123 for(innermostlayer=0; innermostlayer<fNLayers; innermostlayer++) {
db688ae2 1124 if( (drphi < (fSegmentation->GetRadius(innermostlayer)+1.)) ||
1d9af2d5 1125 index[innermostlayer] >= 0 ) break;
1126 }
1127 AliDebug(2,Form(" drphi %f innermost %d",drphi,innermostlayer));
1128 }
1129
1130 Int_t from, to, step;
1131 if (xx > track->GetX()) {
665c0bc7 1132 from=innermostlayer; to=fNLayers;
1d9af2d5 1133 step=+1;
1134 } else {
665c0bc7 1135 from=fNLayers-1; to=innermostlayer-1;
1d9af2d5 1136 step=-1;
1137 }
1138 TString dir = (step>0 ? "outward" : "inward");
1139
1140 for (Int_t ilayer = from; ilayer != to; ilayer += step) {
1141 Double_t r=0.;
db688ae2 1142 r=fSegmentation->GetRadius(ilayer);
1d9af2d5 1143
1144 if (step<0 && xx>r){
1145 break;
1146 }
1147
1d9af2d5 1148
1149 Double_t oldGlobXYZ[3];
1150
1151 if (!track->GetXYZ(oldGlobXYZ)) {
1152 return kFALSE;
1153 }
1154 // continue if we are already beyond this layer
1155 Double_t oldGlobR = TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
1156 if(step>0 && oldGlobR > r){
1157 continue; // going outward
1158 }
1159 if(step<0 && oldGlobR < r) {
1160 continue; // going inward
1161 }
1162 Double_t phi,z;
1163 if (!track->GetPhiZat(r,phi,z)){
1164 return kFALSE;
556b741c 1165 } // only for ITS-SA tracks refit
1d9af2d5 1166 track->SetCheckInvariant(kFALSE);
1167
1168 if (phi<0) phi+=TMath::TwoPi();//from 0 to 360 (rad)
1169 else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();//
0ac80088 1170
9b615954 1171 Int_t trkModule = fSegmentation->GetModule(phi);
1172
1173 Double_t alpha = fSegmentation->GetAlpha(trkModule);
1d9af2d5 1174 if (!track->Propagate(alpha,r)) {
1d9af2d5 1175 return kFALSE;
1176 }
1177
0ac80088 1178 const AliITSRecPointU *clAcc=0;
1d9af2d5 1179 Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
1180
1181 Int_t idx=index[ilayer];
1182 if (idx>=0) { // cluster in this layer
1183 Int_t cli = idx&0x0fffffff;
9b615954 1184 AliITSRecPointU *cl=(AliITSRecPointU *)fLayers[ilayer]->GetCluster(cli);
1d9af2d5 1185
1186 if (cl) {
1187 Int_t cllayer = (idx & 0xf0000000) >> 28;;
9b615954 1188 // edge effect correction
9b615954 1189 if(trkModule!=cl->GetModule()){
1190 if (!track->Propagate(fSegmentation->GetAlpha(cl->GetModule()),fSegmentation->GetRadius(ilayer))) return kFALSE;
1191 }
1d9af2d5 1192 Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
1193 if (chi2<maxchi2) {
1194 clAcc=cl;
1195 maxchi2=chi2;
1196 } else {
f059e2be 1197 return kFALSE;
1d9af2d5 1198 }
1199 }
1200 }
1201
1202 if (clAcc) {
1203 if (!UpdateMI(track,clAcc,maxchi2,idx)){
1204 return kFALSE;
1205 }
1206 track->SetSampledEdx(clAcc->GetQ(),ilayer-2);
1207 }
1208
1209
1210 // Correct for material of the current layer
1211 // cross material
1212 // add time if going outward
1213 if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)){
1214 return kFALSE;
1215 }
1216 track->SetCheckInvariant(kFALSE);
1217 } // end loop on layers
1218
1219 if (!track->PropagateTo(xx,0.,0.)){
1220 return kFALSE;
8d499867 1221 }
556b741c 1222
1d9af2d5 1223 return kTRUE;
1224}
1225
1226//_____________________________________________________________________
0ac80088 1227Int_t AliITStrackerUpgrade::UpdateMI(AliITStrackU* track, const AliITSRecPointU* cl,Double_t chi2,Int_t index) const
1d9af2d5 1228{
1229 //
1230 // Update ITS track
1231 //
1232 Int_t layer = (index & 0xf0000000) >> 28;
1233 track->SetClIndex(layer, index);
1234 if (TMath::Abs(cl->GetQ())<1.e-13) return 0; // ingore the "virtual" clusters
1235
1236
1237 // Take into account the mis-alignment (bring track to cluster plane)
1238 Double_t xTrOrig=track->GetX();
1239 Double_t trxyz[3]; track->GetXYZ(trxyz);
1240 Float_t xclu1 = cl->GetX();//upgrade clusters tracking coordinate
1241 //Float_t yclu1 = cl->GetY();
1242 //Float_t zclu1 = cl->GetZ();
1243 if (!track->Propagate(xTrOrig+xclu1)){
1244 return 0;
1245 }
1246 AliCluster c(*cl);
1247 c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
1248 c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
1249 c.SetSigmaYZ(track->GetSigmaYZ(layer));
1250
1251 // Bring the track back to detector plane in ideal geometry
1252 Int_t updated = track->UpdateMI(&c,chi2,index);
1253 if (!track->Propagate(xTrOrig)) {
1254 return 0;
1255 }
1256 if(!updated){
1257 AliDebug(2,"update failed");
1258 }
1259 return updated;
1260}
1261//____________________________________________________________________-
4211bbc5 1262Int_t AliITStrackerUpgrade::CorrectForLayerMaterial(AliITStrackU *t,
1d9af2d5 1263 Int_t layerindex,
1264 Double_t oldGlobXYZ[3],
1265 TString direction) {
1266 //-------------------------------------------------------------------
1267 // Propagate beyond layer and correct for material
1268 // (material budget in different ways according to fUseTGeo value)
1269 // Add time if going outward (PropagateTo or PropagateToTGeo)
1270 //-------------------------------------------------------------------
1271
1272 // Define budget mode:
1273 // 0: material from AliITSRecoParam (hard coded)
1274 // 1: material from TGeo in stepsof X cm (on the fly)
1275 // X = AliITSRecoParam::GetStepSizeTGeo()
1276 // 2: material from lut
1277 // 3: material from TGeo in one step (same for all hypotheses)
1278 Int_t mode;
1279 switch(fUseTGeo) {
1280 case 0:
1281 mode=0;
1282 break;
1283 case 1:
1284 mode=1;
1285 break;
1286 case 2:
1287 mode=2;
1288 break;
1289 case 3:
1290 if(fTrackingPhase.Contains("Clusters2Tracks"))
1291 { mode=3; } else { mode=1; }
1292 break;
1293 case 4:
1294 if(fTrackingPhase.Contains("Clusters2Tracks"))
1295 { mode=3; } else { mode=2; }
1296 break;
1297 default:
1298 mode=0;
1299 break;
1300 }
1301 if(fTrackingPhase.Contains("Default")) {
1302 mode=0;
1303 }
1304 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
317075a3 1305 //Double_t r = fSegmentation->GetRadius(layerindex);
1306 //Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
1307 Double_t deltar=0.05;
1d9af2d5 1308 Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
1309 Double_t xToGo;
1310 if (!t->GetLocalXat(rToGo,xToGo)) {
1311 return 0;
1312 }
665c0bc7 1313 Int_t index=fNLayers*fCurrentEsdTrack+layerindex;
1d9af2d5 1314
1d9af2d5 1315 Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
1316 Int_t nsteps=1;
1317 // back before material (no correction)
1318 Double_t rOld,xOld;
1319 rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
1320 if (!t->GetLocalXat(rOld,xOld)) return 0;
1321 if (!t->Propagate(xOld)) return 0;
1322
1323 switch(mode) {
1324 case 1:
1325 x0=21.82;
db688ae2 1326 xOverX0 = fSegmentation->GetThickness(layerindex)/x0;
1d9af2d5 1327
1328 lengthTimesMeanDensity = xOverX0*x0;
1329 lengthTimesMeanDensity *= dir;
1330 // Bring the track beyond the material
1331 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)){
1332 return 0;
1333 }
1334 break;
1335 case 0:
1336 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
1337 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
1338 break;
1339 case 2:
1340 if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
1341 xOverX0 = fxOverX0Layer[layerindex];
1342 lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
1343 lengthTimesMeanDensity *= dir;
1344 // Bring the track beyond the material
1345 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1346 break;
1347 case 3:
665c0bc7 1348 if(!fxOverX0LayerTrks || index<0 || index>=fNLayers*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
1d9af2d5 1349 if(fxOverX0LayerTrks[index]<0) {
1350 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
1351 if (!t->PropagateToTGeo(xToGo,nsteps,xOverX0,lengthTimesMeanDensity)) return 0;
1352 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
1353 ((1.-t->GetSnp())*(1.+t->GetSnp())));
1354 fxOverX0LayerTrks[index] = TMath::Abs(xOverX0)/angle;
1355 fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
1356 return 1;
1357 }
1358 xOverX0 = fxOverX0LayerTrks[index];
1359 lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
1360 lengthTimesMeanDensity *= dir;
1361 // Bring the track beyond the material
1362 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)){
1363 return 0;
1364 }
1365 break;
1366 }
1367
1368
1369 return 1;
1370}
1371
1372
1373//_____________________________________________________________________________
0ac80088 1374Double_t AliITStrackerUpgrade::GetPredictedChi2MI(AliITStrackU* track, const AliITSRecPointU *cluster,Int_t layer)
1d9af2d5 1375{
1376 //
1377 // Compute predicted chi2
1378 //
1379 Float_t erry,errz,covyz;
1380 Float_t theta = track->GetTgl();
1381 Float_t phi = track->GetSnp();
1382 phi = TMath::Abs(phi)*TMath::Sqrt(1./((1.-phi)*(1.+phi)));
1383 GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz,covyz);
9b615954 1384 AliDebug(3,Form(" chi2: tr-cl %f %f tr X %f cl X %f",track->GetY()-cluster->GetY(),track->GetZ()-cluster->GetZ(),track->GetX(),cluster->GetX()));
1d9af2d5 1385 // Take into account the mis-alignment (bring track to cluster plane)
1386 Double_t xTrOrig=track->GetX();
1387 if (!track->Propagate(xTrOrig+cluster->GetX())) return 1000.;
9b615954 1388 AliDebug(3,Form(" chi2: tr-cl %f %f tr X %f cl X %f",track->GetY()-cluster->GetY(),track->GetZ()-cluster->GetZ(),track->GetX(),cluster->GetX()));
1d9af2d5 1389 Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz,covyz);
1390 // Bring the track back to detector plane in ideal geometry
1391 // [mis-alignment will be accounted for in UpdateMI()]
1392 if (!track->Propagate(xTrOrig)) return 1000.;
1393 Float_t ny=0.,nz=0.;
1394 //GetNTeor(layer,cluster,theta,phi,ny,nz);
1395 Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
1396 if (delta>1){
1397 chi2+=0.5*TMath::Min(delta/2,2.);
1398 chi2+=2.*cluster->GetDeltaProbability();
1399 }
1400 //
1401 track->SetNy(layer,ny);
1402 track->SetNz(layer,nz);
1403 track->SetSigmaY(layer,erry);
1404 track->SetSigmaZ(layer, errz);
1405 track->SetSigmaYZ(layer,covyz);
4211bbc5 1406
1d9af2d5 1407 track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/((1.-track->GetSnp())*(1.+track->GetSnp()))));
1408 return chi2;
1409}
1410//________________________________________________________________
1411Int_t AliITStrackerUpgrade::GetError(Int_t /*layer*/,
0ac80088 1412 const AliITSRecPointU *cl,
1d9af2d5 1413 Float_t /*tgl*/,Float_t /*tgphitr*/,Float_t /*expQ*/,
1414 Float_t &erry,Float_t &errz,Float_t &covyz,
1415 Bool_t /*addMisalErr*/)
1416{
1417 //
1418 // Calculate cluster position error for the upgrade
1419 //
1420 Int_t retval=0;
1421 covyz=0.;
1422 retval = GetErrorOrigRecPoint(cl,erry,errz,covyz);
1423
1424
1425
1426 return retval;
1427
1428}
1429//____________________________________________________________________-
1430
0ac80088 1431Int_t AliITStrackerUpgrade::GetErrorOrigRecPoint(const AliITSRecPointU *cl,
1d9af2d5 1432 Float_t &erry,Float_t &errz,Float_t &covyz)
1433{
1434 //
1435 // Calculate cluster position error (just take error from AliITSRecPoint)
1436 //
1437 erry = TMath::Sqrt(cl->GetSigmaY2());
1438 errz = TMath::Sqrt(cl->GetSigmaZ2());
1439 covyz = cl->GetSigmaYZ();
1440 return 1;
1441}
1442//__________________________
9b615954 1443//void AliITStrackerUpgrade::GetNTeor(Int_t layer,const AliITSRecPointU* /*cl*/,
1444/* Float_t tgl,Float_t tgphitr,
1d9af2d5 1445 Float_t &ny,Float_t &nz)
1446{
1447 //
1448 // Get "mean shape" (original parametrization from AliITStrackerMI)
1449 //
1450 tgl = TMath::Abs(tgl);
1451 tgphitr = TMath::Abs(tgphitr);
1452
1453 // SPD
1454 if (layer==0) {
1455 ny = 1.+tgphitr*3.2;
1456 nz = 1.+tgl*0.34;
1457 return;
1458 }
1459 if (layer==1) {
1460 ny = 1.+tgphitr*3.2;
1461 nz = 1.+tgl*0.28;
1462 return;
1463 }
1464 // SSD
1465 if (layer==4 || layer==5) {
1466 ny = 2.02+tgphitr*1.95;
1467 nz = 2.02+tgphitr*2.35;
1468 return;
1469 }
1470 // SDD
1471 ny = 6.6-2.7*tgphitr;
1472 nz = 2.8-3.11*tgphitr+0.45*tgl;
1473 return;
9b615954 1474 }*/
1475
1d9af2d5 1476//_________________________________________________________________
1477Int_t AliITStrackerUpgrade::PropagateBack(AliESDEvent *event) {
1478 //--------------------------------------------------------------------
1479 // This functions propagates reconstructed ITS tracks back
1480 // The clusters must be loaded !
1481 //--------------------------------------------------------------------
1482 fTrackingPhase="PropagateBack";
1483 Int_t nentr=event->GetNumberOfTracks();
1484 Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
1485 Int_t ntrk=0;
1486 for (Int_t i=0; i<nentr; i++) {
1487 AliESDtrack *esd=event->GetTrack(i);
1488
1489 if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
1490 if (esd->GetStatus()&AliESDtrack::kITSout) continue;
1491
f059e2be 1492 /*
1493 AliITStrackMI *t=0;
1494 try {
1d9af2d5 1495 t=new AliITStrackMI(*esd);
f059e2be 1496 } catch (const Char_t *msg) {
1d9af2d5 1497 //Warning("PropagateBack",msg);
1498 delete t;
1499 continue;
f059e2be 1500 }
1501 */
4211bbc5 1502 AliITStrackU *t = new AliITStrackU(*esd);
1d9af2d5 1503 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
1504
1505 ResetTrackToFollow(*t);
1506
1507 /*
1508 // propagate to vertex [SR, GSI 17.02.2003]
1509 // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
1510 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
1511 if (fTrackToFollow.PropagateToVertex(event->GetVertex()))
1512 fTrackToFollow.StartTimeIntegral();
1513 // from vertex to outside pipe
1514 CorrectForPipeMaterial(&fTrackToFollow,"outward");
1515 }*/
1516 // Start time integral and add distance from current position to vertex
1517 Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()};
1518 fTrackToFollow.GetXYZ(xyzTrk);
1519 Double_t dst2 = 0.;
1520 for (Int_t icoord=0; icoord<3; icoord++) {
1521 Double_t di = xyzTrk[icoord] - xyzVtx[icoord];
1522 dst2 += di*di;
1523 }
1524 fTrackToFollow.StartTimeIntegral();
1525 fTrackToFollow.AddTimeStep(TMath::Sqrt(dst2));
1526
1527 fTrackToFollow.ResetCovariance(10.);
1528 fTrackToFollow.ResetClusters();
1529 //
665c0bc7 1530 Int_t inx[fNLayers];
1531 for (Int_t k=0; k<fNLayers; k++) inx[k]=-1;
1d9af2d5 1532 Int_t nclusters = t->GetNumberOfClusters();
1533 for(Int_t ncl=0;ncl<nclusters;ncl++){
1534 Int_t index = t-> GetClIndex(ncl);
1535 Int_t lay = (index & 0xf0000000) >> 28;
1536 inx[lay]=index;
1537 }
1538 //
665c0bc7 1539 if (RefitAtBase(fSegmentation->GetRadius(fNLayers-1),&fTrackToFollow,inx)) {
1d9af2d5 1540 //fTrackToFollow.SetLabel(t->GetLabel());//
1541 //fTrackToFollow.CookdEdx();
1542 //CookLabel(&fTrackToFollow,0.); //For comparison only
1543 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
db688ae2 1544 //UseClusters(&fTrackToFollow);
1d9af2d5 1545 ntrk++;
1546 }
1547 delete t;
1548 }
1549
1d9af2d5 1550 fTrackingPhase="Default";
1551
1552 return 0;
1553}
1554//_________________________________________________________________________
1555AliCluster *AliITStrackerUpgrade::GetCluster(Int_t index) const {
1556 //--------------------------------------------------------------------
1557 // Return pointer to a given cluster
1558 //--------------------------------------------------------------------
1559 Int_t l=(index & 0xf0000000) >> 28;
0a64515a 1560 Int_t c=(index & 0x0fffffff) >> 0;
9b615954 1561
1562 AliCluster *cl = fLayers[l]->GetCluster(c);
1563
1564 // AliDebug(2,Form("index %i; cluster index %i; layer %i; Pointer to cluster:%p", index,c,l,cl));
1565
1566 return cl;
1d9af2d5 1567}
1568//______________________________________________________________________________
4211bbc5 1569Int_t AliITStrackerUpgrade::CorrectForPipeMaterial(AliITStrackU *t, TString direction) {
f059e2be 1570 //-------------------------------------------------------------------
1571 // Propagate beyond beam pipe and correct for material
1572 // (material budget in different ways according to fUseTGeo value)
1573 // Add time if going outward (PropagateTo or PropagateToTGeo)
1574 //-------------------------------------------------------------------
1575
1576 // Define budget mode:
1577 // 0: material from AliITSRecoParam (hard coded)
1578 // 1: material from TGeo in one step (on the fly)
1579 // 2: material from lut
1580 // 3: material from TGeo in one step (same for all hypotheses)
1581 Int_t mode;
1582 switch(fUseTGeo) {
1583 case 0:
1584 mode=0;
1585 break;
1586 case 1:
1587 mode=1;
1588 break;
1589 case 2:
1590 mode=2;
1591 break;
1592 case 3:
1593 if(fTrackingPhase.Contains("Clusters2Tracks"))
1594 { mode=3; } else { mode=1; }
1595 break;
1596 case 4:
1597 if(fTrackingPhase.Contains("Clusters2Tracks"))
1598 { mode=3; } else { mode=2; }
1599 break;
1600 default:
1601 mode=0;
1602 break;
1603 }
1604 if(fTrackingPhase.Contains("Default")) mode=0;
1605 Int_t index=fCurrentEsdTrack;
1606
1607 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
1608 Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
1609 Double_t xToGo;
1610 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
1611
1612 Double_t xOverX0,x0,lengthTimesMeanDensity;
1613
1614 switch(mode) {
1615 case 0:
1616 xOverX0 = AliITSRecoParam::GetdPipe();
1617 x0 = AliITSRecoParam::GetX0Be();
1618 lengthTimesMeanDensity = xOverX0*x0;
1619 lengthTimesMeanDensity *= dir;
1620 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1621 break;
1622 case 1:
1623 if (!t->PropagateToTGeo(xToGo,1)) return 0;
1624 break;
1625 case 2:
1626 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
1627 xOverX0 = fxOverX0Pipe;
1628 lengthTimesMeanDensity = fxTimesRhoPipe;
1629 lengthTimesMeanDensity *= dir;
1630 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1631 break;
1632 case 3:
1633 if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
1634 if(fxOverX0PipeTrks[index]<0) {
1635 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
1636 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
1637 ((1.-t->GetSnp())*(1.+t->GetSnp())));
1638 fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
1639 fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
1640 return 1;
1641 }
1642 xOverX0 = fxOverX0PipeTrks[index];
1643 lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
1644 lengthTimesMeanDensity *= dir;
1645 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1646 break;
1647 }
1648
1649 return 1;
1d9af2d5 1650}
1651//__________________________________________________________________________
1652Int_t AliITStrackerUpgrade::RefitInward(AliESDEvent *event) {
1653 //--------------------------------------------------------------------
1654 // This functions refits ITS tracks using the
1655 // "inward propagated" TPC tracks
1656 // The clusters must be loaded !
1657 //--------------------------------------------------------------------
1658 fTrackingPhase="RefitInward";
1659 Int_t nentr=event->GetNumberOfTracks();
1660 for (Int_t i=0; i<nentr; i++) {
1661 AliESDtrack *esd=event->GetTrack(i);
1662
1663 // if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
1664 // if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
1665 //if (esd->GetStatus()&AliESDtrack::kTPCout)
1666 // if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
4211bbc5 1667 AliITStrackU *t=0;
1d9af2d5 1668 try {
4211bbc5 1669 t=new AliITStrackU(*esd);
1d9af2d5 1670 } catch (const Char_t *msg) {
1671 //Warning("RefitInward",msg);
1672 delete t;
1673 continue;
1674 }
1675
1676 ResetTrackToFollow(*t);
1677 // fTrackToFollow.ResetClusters();
1678 // if ((esd->GetStatus()&AliESDtrack::kTPCin)==0)
1679 // fTrackToFollow.ResetCovariance(10.);
1680 //Refitting...
1681 //The beam pipe
f059e2be 1682 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
1683 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
1684 AliESDtrack *esdTrack =fTrackToFollow.GetESDtrack();
1685 Double_t r[3]={0.,0.,0.};
1686 Double_t maxD=3.;
1687 esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
1688 }
1d9af2d5 1689 delete t;
1690 }
1691
1692 fTrackingPhase="Default";
1693 return 0;
1694}
1695
1696
1697