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