Improved cuts for the reconstruction of V0s (M.Ivanov)
[u/mrichter/AliRoot.git] / ITS / AliITStrackerV2.cxx
CommitLineData
006b5f7f 1/**************************************************************************
2 * Copyright(c) 1998-1999, 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//-------------------------------------------------------------------------
17// Implementation of the ITS tracker class
18856a77 18// It reads AliITSclusterV2 clusters and creates AliITStrackV2 tracks
19// and fills with them the ESD
006b5f7f 20// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
a9a2d814 21// dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
006b5f7f 22//-------------------------------------------------------------------------
1f9a65c4 23
9c369bfe 24#include <new>
25
006b5f7f 26#include <TFile.h>
27#include <TTree.h>
28#include <TRandom.h>
006b5f7f 29
30#include "AliITSgeom.h"
31#include "AliITSRecPoint.h"
517b130f 32#include "AliTPCtrack.h"
c630aafd 33#include "AliESD.h"
006b5f7f 34#include "AliITSclusterV2.h"
35#include "AliITStrackerV2.h"
36
8676d691 37ClassImp(AliITStrackerV2)
006b5f7f 38
18856a77 39AliITStrackerV2::AliITSlayer AliITStrackerV2::fgLayers[kMaxLayer]; // ITS layers
006b5f7f 40
c630aafd 41AliITStrackerV2::AliITStrackerV2(const AliITSgeom *geom) : AliTracker() {
006b5f7f 42 //--------------------------------------------------------------------
61ab8ea8 43 //This is the AliITStrackerV2 constructor
006b5f7f 44 //--------------------------------------------------------------------
006b5f7f 45 AliITSgeom *g=(AliITSgeom*)geom;
46
47 Float_t x,y,z;
48 Int_t i;
c630aafd 49 for (i=1; i<kMaxLayer+1; i++) {
006b5f7f 50 Int_t nlad=g->GetNladders(i);
51 Int_t ndet=g->GetNdetectors(i);
52
53 g->GetTrans(i,1,1,x,y,z);
54 Double_t r=TMath::Sqrt(x*x + y*y);
55 Double_t poff=TMath::ATan2(y,x);
56 Double_t zoff=z;
57
58 g->GetTrans(i,1,2,x,y,z);
59 r += TMath::Sqrt(x*x + y*y);
60 g->GetTrans(i,2,1,x,y,z);
61 r += TMath::Sqrt(x*x + y*y);
62 g->GetTrans(i,2,2,x,y,z);
63 r += TMath::Sqrt(x*x + y*y);
64 r*=0.25;
65
18856a77 66 new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
006b5f7f 67
68 for (Int_t j=1; j<nlad+1; j++) {
69 for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
70 Float_t x,y,zshift; g->GetTrans(i,j,k,x,y,zshift);
71 Double_t rot[9]; g->GetRotMatrix(i,j,k,rot);
72
f363d377 73 Double_t phi=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
74 phi+=TMath::Pi()/2;
75 if (i==1) phi+=TMath::Pi();
76 Double_t cp=TMath::Cos(phi), sp=TMath::Sin(phi);
77 Double_t r=x*cp+y*sp;
006b5f7f 78
18856a77 79 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
006b5f7f 80 new(&det) AliITSdetector(r,phi);
81 }
82 }
c630aafd 83
006b5f7f 84 }
85
61ab8ea8 86 fI=kMaxLayer;
743a19f2 87
61ab8ea8 88 fPass=0;
89 fConstraint[0]=1; fConstraint[1]=0;
8676d691 90
91 Double_t xyz[]={kXV,kYV,kZV}, ers[]={kSigmaXV,kSigmaYV,kSigmaZV};
92 SetVertex(xyz,ers);
c0dd5278 93
94 for (Int_t i=0; i<kMaxLayer; i++) fLayersNotToSkip[i]=kLayersNotToSkip[i];
95 fLastLayerToTrackTo=kLastLayerToTrackTo;
96
97}
98
99void AliITStrackerV2::SetLayersNotToSkip(Int_t *l) {
100 //--------------------------------------------------------------------
101 //This function set masks of the layers which must be not skipped
102 //--------------------------------------------------------------------
103 for (Int_t i=0; i<kMaxLayer; i++) fLayersNotToSkip[i]=l[i];
61ab8ea8 104}
006b5f7f 105
c630aafd 106Int_t AliITStrackerV2::LoadClusters(TTree *cTree) {
61ab8ea8 107 //--------------------------------------------------------------------
108 //This function loads ITS clusters
109 //--------------------------------------------------------------------
61ab8ea8 110 TBranch *branch=cTree->GetBranch("Clusters");
111 if (!branch) {
c630aafd 112 Error("LoadClusters"," can't get the branch !\n");
8676d691 113 return 1;
61ab8ea8 114 }
006b5f7f 115
61ab8ea8 116 TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy;
117 branch->SetAddress(&clusters);
006b5f7f 118
61ab8ea8 119 Int_t j=0;
babd135a 120 Int_t detector=0;
61ab8ea8 121 for (Int_t i=0; i<kMaxLayer; i++) {
18856a77 122 Int_t ndet=fgLayers[i].GetNdetectors();
123 Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
61ab8ea8 124 for (; j<jmax; j++) {
125 if (!cTree->GetEvent(j)) continue;
126 Int_t ncl=clusters->GetEntriesFast();
127 while (ncl--) {
128 AliITSclusterV2 *c=(AliITSclusterV2*)clusters->UncheckedAt(ncl);
babd135a 129 detector = c->GetDetectorIndex();
18856a77 130 fgLayers[i].InsertCluster(new AliITSclusterV2(*c));
61ab8ea8 131 }
132 clusters->Delete();
babd135a 133 //add dead zone virtual "cluster"
134
135 if (i<1){
136 for (Float_t ydead = -2.; ydead < 2. ; ydead+=0.05){
137 Int_t lab[4] = {0,0,0,detector};
138 Int_t info[3] = {0,0,0};
139 Float_t hit[5]={ydead,0,1,0.01,0};
140 fgLayers[i].InsertCluster(new AliITSclusterV2(lab, hit, info));
141 //
142 hit[1]=-7.;
143 fgLayers[i].InsertCluster(new AliITSclusterV2(lab, hit, info));
144 }
145 }
146
61ab8ea8 147 }
babd135a 148 //
18856a77 149 fgLayers[i].ResetRoad(); //road defined by the cluster density
006b5f7f 150 }
c630aafd 151
8676d691 152 return 0;
61ab8ea8 153}
006b5f7f 154
61ab8ea8 155void AliITStrackerV2::UnloadClusters() {
156 //--------------------------------------------------------------------
157 //This function unloads ITS clusters
158 //--------------------------------------------------------------------
18856a77 159 for (Int_t i=0; i<kMaxLayer; i++) fgLayers[i].ResetClusters();
006b5f7f 160}
161
880f41b9 162static Int_t CorrectForDeadZoneMaterial(AliITStrackV2 *t) {
163 //--------------------------------------------------------------------
164 // Correction for the material between the TPC and the ITS
165 // (should it belong to the TPC code ?)
166 //--------------------------------------------------------------------
167 Double_t riw=80., diw=0.0053, x0iw=30; // TPC inner wall ?
168 Double_t rcd=61., dcd=0.0053, x0cd=30; // TPC "central drum" ?
169 Double_t yr=12.8, dr=0.03; // rods ?
170 Double_t zm=0.2, dm=0.40; // membrane
171 //Double_t rr=52., dr=0.19, x0r=24., yyr=7.77; //rails
172 Double_t rs=50., ds=0.001; // something belonging to the ITS (screen ?)
173
174 if (t->GetX() > riw) {
175 if (!t->PropagateTo(riw,diw,x0iw)) return 1;
176 if (TMath::Abs(t->GetY())>yr) t->CorrectForMaterial(dr);
177 if (TMath::Abs(t->GetZ())<zm) t->CorrectForMaterial(dm);
178 if (!t->PropagateTo(rcd,dcd,x0cd)) return 1;
179 //Double_t x,y,z; t->GetGlobalXYZat(rr,x,y,z);
180 //if (TMath::Abs(y)<yyr) t->PropagateTo(rr,dr,x0r);
181 if (!t->PropagateTo(rs,ds)) return 1;
182 } else if (t->GetX() < rs) {
183 if (!t->PropagateTo(rs,-ds)) return 1;
184 //Double_t x,y,z; t->GetGlobalXYZat(rr,x,y,z);
185 //if (TMath::Abs(y)<yyr) t->PropagateTo(rr,-dr,x0r);
186 if (!t->PropagateTo(rcd,-dcd,x0cd)) return 1;
ba1a3e1a 187 if (!t->PropagateTo(riw+0.001,-diw,x0iw)) return 1;
880f41b9 188 } else {
189 ::Error("CorrectForDeadZoneMaterial","track is already in the dead zone !");
190 return 1;
191 }
192
193 return 0;
194}
195
c630aafd 196Int_t AliITStrackerV2::Clusters2Tracks(AliESD *event) {
006b5f7f 197 //--------------------------------------------------------------------
c630aafd 198 // This functions reconstructs ITS tracks
199 // The clusters must be already loaded !
006b5f7f 200 //--------------------------------------------------------------------
c630aafd 201 TObjArray itsTracks(15000);
006b5f7f 202
c630aafd 203 {/* Read ESD tracks */
204 Int_t nentr=event->GetNumberOfTracks();
205 Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
206 while (nentr--) {
207 AliESDtrack *esd=event->GetTrack(nentr);
208
ba1a3e1a 209 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
210 if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
211 if (esd->GetStatus()&AliESDtrack::kITSin) continue;
c630aafd 212
213 AliITStrackV2 *t=0;
214 try {
215 t=new AliITStrackV2(*esd);
216 } catch (const Char_t *msg) {
217 Warning("Clusters2Tracks",msg);
218 delete t;
219 continue;
220 }
babd135a 221 if (TMath::Abs(t->GetD())>5) {
ebd8369c 222 delete t;
223 continue;
224 }
006b5f7f 225
c630aafd 226 if (CorrectForDeadZoneMaterial(t)!=0) {
227 Warning("Clusters2Tracks",
228 "failed to correct for the material in the dead zone !\n");
229 delete t;
230 continue;
231 }
babd135a 232 t->fReconstructed = kFALSE;
c630aafd 233 itsTracks.AddLast(t);
234 }
235 } /* End Read ESD tracks */
236
237 itsTracks.Sort();
238 Int_t nentr=itsTracks.GetEntriesFast();
de5d3595 239 fTrackHypothesys.Expand(nentr);
c630aafd 240 Int_t ntrk=0;
241 for (fPass=0; fPass<2; fPass++) {
242 Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
243 for (Int_t i=0; i<nentr; i++) {
babd135a 244 cerr<<fPass<<" "<<i<<'\r';
de5d3595 245 fCurrentEsdTrack = i;
c630aafd 246 AliITStrackV2 *t=(AliITStrackV2*)itsTracks.UncheckedAt(i);
babd135a 247 if (t==0) continue; //this track has been already tracked
248 if (t->fReconstructed) continue; //this track was already "succesfully" reconstructed
249 if ( (TMath::Abs(t->GetD(GetX(),GetY())) >2.) && fConstraint[fPass]) continue;
250 if ( (TMath::Abs(t->GetZat(GetX())-GetZ())>2.) && fConstraint[fPass]) continue;
251
c630aafd 252 Int_t tpcLabel=t->GetLabel(); //save the TPC track label
253
254 ResetTrackToFollow(*t);
255 ResetBestTrack();
256
257 for (FollowProlongation(); fI<kMaxLayer; fI++) {
de5d3595 258 while (TakeNextProlongation()) {
259 FollowProlongation();
260 }
c630aafd 261 }
262
babd135a 263 SortTrackHypothesys(fCurrentEsdTrack,0.9); //MI change
264 CompressTrackHypothesys(fCurrentEsdTrack,0.90,2); //MI change
c630aafd 265
de5d3595 266 /*
267 if (fBestTrack.GetNumberOfClusters() == 0) continue;
268
c630aafd 269 if (fConstraint[fPass]) {
de5d3595 270 ResetTrackToFollow(*t);
271 if (!RefitAt(3.7, &fTrackToFollow, &fBestTrack)) continue;
272 ResetBestTrack();
273 }
274
c630aafd 275 fBestTrack.SetLabel(tpcLabel);
276 fBestTrack.CookdEdx();
277 CookLabel(&fBestTrack,0.); //For comparison only
278 fBestTrack.UpdateESDtrack(AliESDtrack::kITSin);
de5d3595 279 //UseClusters(&fBestTrack);
c630aafd 280 delete itsTracks.RemoveAt(i);
281 ntrk++;
de5d3595 282 */
babd135a 283 AliITStrackV2 * besttrack = GetBestHypothesys(fCurrentEsdTrack,t,5);
de5d3595 284 if (!besttrack) continue;
285 besttrack->SetLabel(tpcLabel);
286 besttrack->CookdEdx();
babd135a 287 besttrack->fFakeRatio=1.;
de5d3595 288 CookLabel(besttrack,0.); //For comparison only
babd135a 289 // besttrack->UpdateESDtrack(AliESDtrack::kITSin);
290
291 if (besttrack->GetChi2()/besttrack->GetNumberOfClusters()>3.5){
292 if ( (TMath::Abs(besttrack->GetD(GetX(),GetY()))>0.4) && fConstraint[fPass]) {
293 CompressTrackHypothesys(fCurrentEsdTrack,0.0,0);
294 continue;
295 }
296 if ( (TMath::Abs(besttrack->GetZat(GetX()) -GetZ() )>0.4) && fConstraint[fPass]){
297 CompressTrackHypothesys(fCurrentEsdTrack,0.0,0);
298 continue;
299 }
300 }
301
302 //delete itsTracks.RemoveAt(i);
303 t->fReconstructed = kTRUE;
de5d3595 304 ntrk++;
305
c630aafd 306 }
307 }
babd135a 308
309 for (Int_t i=0; i<nentr; i++) {
310 fCurrentEsdTrack = i;
311 AliITStrackV2 *t=(AliITStrackV2*)itsTracks.UncheckedAt(i);
312 if (t==0) continue;
313 Int_t tpcLabel=t->GetLabel(); //save the TPC track label
314 cerr<<i<<'\r';
315 AliITStrackV2 * besttrack = GetBestHypothesysMIP(fCurrentEsdTrack,t);
316 if (!besttrack) continue;
317
318 besttrack->SetLabel(tpcLabel);
319 besttrack->CookdEdx();
320 besttrack->fFakeRatio=1.;
321 CookLabel(besttrack,0.); //For comparison only
322 besttrack->UpdateESDtrack(AliESDtrack::kITSin);
323 }
324 //
325
c630aafd 326 itsTracks.Delete();
de5d3595 327 //
328 Int_t entries = fTrackHypothesys.GetEntriesFast();
329 for (Int_t ientry=0;ientry<entries;ientry++){
330 TObjArray * array =(TObjArray*)fTrackHypothesys.UncheckedAt(ientry);
2257f27e 331 if (array) array->Delete();
332 delete fTrackHypothesys.RemoveAt(ientry);
de5d3595 333 }
2257f27e 334
de5d3595 335 fTrackHypothesys.Delete();
c630aafd 336 Info("Clusters2Tracks","Number of prolonged tracks: %d\n",ntrk);
337
338 return 0;
339}
340
341Int_t AliITStrackerV2::Clusters2Tracks(TTree *tpcTree, TTree *itsTree) {
342 //--------------------------------------------------------------------
343 // This functions reconstructs ITS tracks
344 // The clusters must be already loaded !
345 //--------------------------------------------------------------------
a7554aa9 346 Int_t nentr=0; TObjArray itsTracks(15000);
7f6ddf58 347
df02fd67 348 Warning("Clusters2Tracks(TTree *, TTree *)",
349 "Will be removed soon ! Use Clusters2Tracks(AliESD *) instead.");
350
7f6ddf58 351 {/* Read TPC tracks */
7f6ddf58 352 AliTPCtrack *itrack=new AliTPCtrack;
c630aafd 353 TBranch *branch=tpcTree->GetBranch("tracks");
354 if (!branch) {
355 Error("Clusters2Tracks","Can't get the branch !");
356 return 1;
357 }
7f6ddf58 358 tpcTree->SetBranchAddress("tracks",&itrack);
359 nentr=(Int_t)tpcTree->GetEntries();
c630aafd 360
880f41b9 361 Info("Clusters2Tracks","Number of TPC tracks: %d\n",nentr);
362
7f6ddf58 363 for (Int_t i=0; i<nentr; i++) {
364 tpcTree->GetEvent(i);
4ab260d6 365 AliITStrackV2 *t=0;
366 try {
367 t=new AliITStrackV2(*itrack);
368 } catch (const Char_t *msg) {
8676d691 369 Warning("Clusters2Tracks",msg);
4ab260d6 370 delete t;
371 continue;
372 }
ebd8369c 373 if (TMath::Abs(t->GetD())>4) {
374 delete t;
375 continue;
376 }
a9a2d814 377
880f41b9 378 if (CorrectForDeadZoneMaterial(t)!=0) {
379 Warning("Clusters2Tracks",
380 "failed to correct for the material in the dead zone !\n");
ebd8369c 381 delete t;
880f41b9 382 continue;
383 }
a9a2d814 384
385 itsTracks.AddLast(t);
7f6ddf58 386 }
7f6ddf58 387 delete itrack;
006b5f7f 388 }
7f6ddf58 389 itsTracks.Sort();
880f41b9 390 nentr=itsTracks.GetEntriesFast();
006b5f7f 391
7f6ddf58 392
006b5f7f 393 AliITStrackV2 *otrack=&fBestTrack;
c630aafd 394 TBranch *branch=itsTree->GetBranch("tracks");
395 if (!branch) itsTree->Branch("tracks","AliITStrackV2",&otrack,32000,3);
396 else branch->SetAddress(&otrack);
006b5f7f 397
7f6ddf58 398 for (fPass=0; fPass<2; fPass++) {
399 Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
400 for (Int_t i=0; i<nentr; i++) {
7f6ddf58 401 AliITStrackV2 *t=(AliITStrackV2*)itsTracks.UncheckedAt(i);
402 if (t==0) continue; //this track has been already tracked
403 Int_t tpcLabel=t->GetLabel(); //save the TPC track label
a9a2d814 404
4ab260d6 405 ResetTrackToFollow(*t);
7f6ddf58 406 ResetBestTrack();
006b5f7f 407
c630aafd 408 for (FollowProlongation(); fI<kMaxLayer; fI++) {
409 while (TakeNextProlongation()) FollowProlongation();
410 }
006b5f7f 411
c0dd5278 412 if (fBestTrack.GetNumberOfClusters() == 0) continue;
8676d691 413
414 if (fConstraint[fPass]) {
c630aafd 415 ResetTrackToFollow(*t);
416 if (!RefitAt(3.7, &fTrackToFollow, &fBestTrack)) continue;
417 ResetBestTrack();
8676d691 418 }
4ab260d6 419
8676d691 420 fBestTrack.SetLabel(tpcLabel);
421 fBestTrack.CookdEdx();
422 CookLabel(&fBestTrack,0.); //For comparison only
c630aafd 423 itsTree->Fill();
de5d3595 424 //UseClusters(&fBestTrack);
8676d691 425 delete itsTracks.RemoveAt(i);
7f6ddf58 426 }
006b5f7f 427 }
006b5f7f 428
c630aafd 429 nentr=(Int_t)itsTree->GetEntries();
880f41b9 430 Info("Clusters2Tracks","Number of prolonged tracks: %d\n",nentr);
8676d691 431
7f6ddf58 432 itsTracks.Delete();
61ab8ea8 433
c630aafd 434 return 0;
435}
436
437Int_t AliITStrackerV2::PropagateBack(AliESD *event) {
438 //--------------------------------------------------------------------
439 // This functions propagates reconstructed ITS tracks back
440 // The clusters must be loaded !
441 //--------------------------------------------------------------------
442 Int_t nentr=event->GetNumberOfTracks();
443 Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
444
445 Int_t ntrk=0;
446 for (Int_t i=0; i<nentr; i++) {
447 AliESDtrack *esd=event->GetTrack(i);
448
ba1a3e1a 449 if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
450 if (esd->GetStatus()&AliESDtrack::kITSout) continue;
c630aafd 451
452 AliITStrackV2 *t=0;
453 try {
454 t=new AliITStrackV2(*esd);
455 } catch (const Char_t *msg) {
456 Warning("PropagateBack",msg);
457 delete t;
458 continue;
459 }
460
461 ResetTrackToFollow(*t);
462
463 // propagete to vertex [SR, GSI 17.02.2003]
a9459f9c 464 // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
465 if (fTrackToFollow.PropagateTo(3.,0.0028,65.19)) {
466 if (fTrackToFollow.PropagateToVertex()) {
467 fTrackToFollow.StartTimeIntegral();
468 }
469 fTrackToFollow.PropagateTo(3.,-0.0028,65.19);
470 }
c630aafd 471
472 fTrackToFollow.ResetCovariance(); fTrackToFollow.ResetClusters();
473 if (RefitAt(49.,&fTrackToFollow,t)) {
474 if (CorrectForDeadZoneMaterial(&fTrackToFollow)!=0) {
475 Warning("PropagateBack",
476 "failed to correct for the material in the dead zone !\n");
477 delete t;
478 continue;
479 }
480 fTrackToFollow.SetLabel(t->GetLabel());
ba1a3e1a 481 //fTrackToFollow.CookdEdx();
c630aafd 482 CookLabel(&fTrackToFollow,0.); //For comparison only
483 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
de5d3595 484 //UseClusters(&fTrackToFollow);
c630aafd 485 ntrk++;
486 }
487 delete t;
488 }
489
490 Info("PropagateBack","Number of back propagated ITS tracks: %d\n",ntrk);
61ab8ea8 491
83d73500 492 return 0;
493}
494
c630aafd 495Int_t AliITStrackerV2::RefitInward(AliESD *event) {
83d73500 496 //--------------------------------------------------------------------
497 // This functions refits ITS tracks using the
498 // "inward propagated" TPC tracks
c630aafd 499 // The clusters must be loaded !
83d73500 500 //--------------------------------------------------------------------
c630aafd 501 Int_t nentr=event->GetNumberOfTracks();
502 Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
503
504 Int_t ntrk=0;
505 for (Int_t i=0; i<nentr; i++) {
506 AliESDtrack *esd=event->GetTrack(i);
507
ba1a3e1a 508 if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
509 if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
510 if (esd->GetStatus()&AliESDtrack::kTPCout)
2257f27e 511 if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
c630aafd 512
513 AliITStrackV2 *t=0;
514 try {
515 t=new AliITStrackV2(*esd);
516 } catch (const Char_t *msg) {
517 Warning("RefitInward",msg);
518 delete t;
519 continue;
520 }
521
522 if (CorrectForDeadZoneMaterial(t)!=0) {
523 Warning("RefitInward",
524 "failed to correct for the material in the dead zone !\n");
525 delete t;
526 continue;
527 }
528
529 ResetTrackToFollow(*t);
530 fTrackToFollow.ResetClusters();
531
2257f27e 532 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0)
533 fTrackToFollow.ResetCovariance();
534
c630aafd 535 //Refitting...
536 if (RefitAt(3.7, &fTrackToFollow, t)) {
537 fTrackToFollow.SetLabel(t->GetLabel());
538 fTrackToFollow.CookdEdx();
de5d3595 539 CookLabel(&fTrackToFollow,0.0); //For comparison only
18856a77 540
67c3dcbe 541 if (fTrackToFollow.PropagateTo(3.,0.0028,65.19)) {//The beam pipe
542 Double_t a=fTrackToFollow.GetAlpha();
543 Double_t cs=TMath::Cos(a),sn=TMath::Sin(a);
544 Double_t xv= GetX()*cs + GetY()*sn;
545 Double_t yv=-GetX()*sn + GetY()*cs;
546
547 Double_t c=fTrackToFollow.GetC(), snp=fTrackToFollow.GetSnp();
548 Double_t x=fTrackToFollow.GetX(), y=fTrackToFollow.GetY();
549 Double_t tgfv=-(c*(x-xv)-snp)/(c*(y-yv) + TMath::Sqrt(1.-snp*snp));
550 Double_t fv=TMath::ATan(tgfv);
551
552 cs=TMath::Cos(fv); sn=TMath::Sin(fv);
553 x = xv*cs + yv*sn;
554 yv=-xv*sn + yv*cs; xv=x;
555
556 if (fTrackToFollow.Propagate(fv+a,xv)) {
df02fd67 557 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
de5d3595 558 //UseClusters(&fTrackToFollow);
67c3dcbe 559 {
560 AliITSclusterV2 c; c.SetY(yv); c.SetZ(GetZ());
561 c.SetSigmaY2(GetSigmaY()*GetSigmaY());
562 c.SetSigmaZ2(GetSigmaZ()*GetSigmaZ());
563 Double_t chi2=fTrackToFollow.GetPredictedChi2(&c);
564 if (chi2<kMaxChi2)
565 if (fTrackToFollow.Update(&c,-chi2,0))
566 fTrackToFollow.SetConstrainedESDtrack(chi2);
567 }
df02fd67 568 ntrk++;
569 }
67c3dcbe 570 }
c630aafd 571 }
572 delete t;
573 }
574
575 Info("RefitInward","Number of refitted tracks: %d\n",ntrk);
83d73500 576
c630aafd 577 return 0;
578}
83d73500 579
006b5f7f 580AliCluster *AliITStrackerV2::GetCluster(Int_t index) const {
581 //--------------------------------------------------------------------
582 // Return pointer to a given cluster
583 //--------------------------------------------------------------------
584 Int_t l=(index & 0xf0000000) >> 28;
585 Int_t c=(index & 0x0fffffff) >> 00;
18856a77 586 return fgLayers[l].GetCluster(c);
006b5f7f 587}
588
589
590void AliITStrackerV2::FollowProlongation() {
591 //--------------------------------------------------------------------
592 //This function finds a track prolongation
593 //--------------------------------------------------------------------
c0dd5278 594 while (fI>fLastLayerToTrackTo) {
a9a2d814 595 Int_t i=fI-1;
8676d691 596
18856a77 597 AliITSlayer &layer=fgLayers[i];
a9a2d814 598 AliITStrackV2 &track=fTracks[i];
006b5f7f 599
14825d5a 600 Double_t r=layer.GetR();
61ab8ea8 601
a9a2d814 602 if (i==3 || i==1) {
18856a77 603 Double_t rs=0.5*(fgLayers[i+1].GetR() + r);
61ab8ea8 604 Double_t d=0.0034, x0=38.6;
605 if (i==1) {rs=9.; d=0.0097; x0=42;}
606 if (!fTrackToFollow.PropagateTo(rs,d,x0)) {
8676d691 607 //Warning("FollowProlongation","propagation failed !\n");
c0dd5278 608 return;
a9a2d814 609 }
006b5f7f 610 }
611
612 //find intersection
613 Double_t x,y,z;
14825d5a 614 if (!fTrackToFollow.GetGlobalXYZat(r,x,y,z)) {
8676d691 615 //Warning("FollowProlongation","failed to estimate track !\n");
c0dd5278 616 return;
006b5f7f 617 }
618 Double_t phi=TMath::ATan2(y,x);
f363d377 619
006b5f7f 620 Int_t idet=layer.FindDetectorIndex(phi,z);
621 if (idet<0) {
8676d691 622 //Warning("FollowProlongation","failed to find a detector !\n");
c0dd5278 623 return;
006b5f7f 624 }
625
626 //propagate to the intersection
627 const AliITSdetector &det=layer.GetDetector(idet);
14825d5a 628 phi=det.GetPhi();
a9a2d814 629 if (!fTrackToFollow.Propagate(phi,det.GetR())) {
8676d691 630 //Warning("FollowProlongation","propagation failed !\n");
c0dd5278 631 return;
006b5f7f 632 }
633 fTrackToFollow.SetDetectorIndex(idet);
634
635 //Select possible prolongations and store the current track estimation
636 track.~AliITStrackV2(); new(&track) AliITStrackV2(fTrackToFollow);
8676d691 637 Double_t dz=7*TMath::Sqrt(track.GetSigmaZ2() + kSigmaZ2[i]);
638 Double_t dy=7*TMath::Sqrt(track.GetSigmaY2() + kSigmaY2[i]);
639 Double_t road=layer.GetRoad();
640 if (dz*dy>road*road) {
641 Double_t dd=TMath::Sqrt(dz*dy), scz=dz/dd, scy=dy/dd;
642 dz=road*scz; dy=road*scy;
643 }
b87bd7cd 644
645 //Double_t dz=4*TMath::Sqrt(track.GetSigmaZ2() + kSigmaZ2[i]);
14825d5a 646 if (dz < 0.5*TMath::Abs(track.GetTgl())) dz=0.5*TMath::Abs(track.GetTgl());
7f6ddf58 647 if (dz > kMaxRoad) {
8676d691 648 //Warning("FollowProlongation","too broad road in Z !\n");
c0dd5278 649 return;
006b5f7f 650 }
a9a2d814 651
babd135a 652 // if (TMath::Abs(fTrackToFollow.GetZ()-GetZ()) > r+dz) return;
a9a2d814 653
b87bd7cd 654 //Double_t dy=4*TMath::Sqrt(track.GetSigmaY2() + kSigmaY2[i]);
babd135a 655 if (TMath::Abs(track.GetSnp()>kMaxSnp)) {
656 fI--;
657 continue; // MI
658 }
14825d5a 659 if (dy < 0.5*TMath::Abs(track.GetSnp())) dy=0.5*TMath::Abs(track.GetSnp());
7f6ddf58 660 if (dy > kMaxRoad) {
8676d691 661 //Warning("FollowProlongation","too broad road in Y !\n");
c0dd5278 662 return;
006b5f7f 663 }
a9a2d814 664
7f6ddf58 665 Double_t zmin=track.GetZ() - dz;
006b5f7f 666 Double_t zmax=track.GetZ() + dz;
14825d5a 667 Double_t ymin=track.GetY() + r*phi - dy;
668 Double_t ymax=track.GetY() + r*phi + dy;
a9a2d814 669 layer.SelectClusters(zmin,zmax,ymin,ymax);
670 fI--;
006b5f7f 671
006b5f7f 672 //take another prolongation
de5d3595 673 if (!TakeNextProlongation()){
babd135a 674 //skipped++;
675 fTrackToFollow.fNSkipped++;
676 if (fLayersNotToSkip[fI]||fTrackToFollow.fNSkipped>1) return;
de5d3595 677 }
babd135a 678 if (fTrackToFollow.fNUsed>1) return;
679 if (fTrackToFollow.fNUsed+fTrackToFollow.fNSkipped>1) return;
680 if ( (fI<3) && ( fTrackToFollow.GetChi2()/fTrackToFollow.GetNumberOfClusters()>kChi2PerCluster)) return;
006b5f7f 681 }
682
683 //deal with the best track
684 Int_t ncl=fTrackToFollow.GetNumberOfClusters();
685 Int_t nclb=fBestTrack.GetNumberOfClusters();
de5d3595 686 if (ncl){
babd135a 687 if (ncl<4) return;
688 if ( (ncl<6) && (fTrackToFollow.GetChi2()/float(ncl))>3) return;
689 if (fTrackToFollow.GetChi2()/ncl>5.5) return;
690 fTrackToFollow.CookdEdx();
691 if (fTrackToFollow.fESDtrack->GetTPCsignal()>80.)
692 if ((fTrackToFollow.GetdEdx()/fTrackToFollow.fESDtrack->GetTPCsignal())<0.35){
693 // mismatch in dedx
694 return;
695 }
696 //
697 fTrackToFollow.SetLabel(fTrackToFollow.fESDtrack->GetLabel());
698 CookLabel(&fTrackToFollow,0.); //
699 //
700 if ( (nclb>3) && ((fTrackToFollow.GetChi2()/ncl)<(3*fBestTrack.GetChi2()/(nclb))))
de5d3595 701 AddTrackHypothesys(new AliITStrackV2(fTrackToFollow), fCurrentEsdTrack);
babd135a 702 else
703 if (ncl>3) AddTrackHypothesys(new AliITStrackV2(fTrackToFollow), fCurrentEsdTrack);
de5d3595 704 if (ncl >= nclb) {
705 Double_t chi2=fTrackToFollow.GetChi2();
706 if (chi2/ncl < kChi2PerCluster) {
707 if (ncl > nclb ) {
708 ResetBestTrack();
006b5f7f 709 }
de5d3595 710 if ( (ncl == nclb) && chi2 < fBestTrack.GetChi2()) {
711 ResetBestTrack();
712 }
713 }
714 }
006b5f7f 715 }
006b5f7f 716}
717
006b5f7f 718Int_t AliITStrackerV2::TakeNextProlongation() {
719 //--------------------------------------------------------------------
a9a2d814 720 // This function takes another track prolongation
721 //
722 // dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
006b5f7f 723 //--------------------------------------------------------------------
18856a77 724 AliITSlayer &layer=fgLayers[fI];
4ab260d6 725 ResetTrackToFollow(fTracks[fI]);
006b5f7f 726
b87bd7cd 727 Double_t dz=7*TMath::Sqrt(fTrackToFollow.GetSigmaZ2() + kSigmaZ2[fI]);
728 Double_t dy=7*TMath::Sqrt(fTrackToFollow.GetSigmaY2() + kSigmaY2[fI]);
8676d691 729 Double_t road=layer.GetRoad();
730 if (dz*dy>road*road) {
731 Double_t dd=TMath::Sqrt(dz*dy), scz=dz/dd, scy=dy/dd;
732 dz=road*scz; dy=road*scy;
733 }
006b5f7f 734
735 const AliITSclusterV2 *c=0; Int_t ci=-1;
736 Double_t chi2=12345.;
737 while ((c=layer.GetNextCluster(ci))!=0) {
006b5f7f 738 Int_t idet=c->GetDetectorIndex();
4ab260d6 739 if (fTrackToFollow.GetDetectorIndex()!=idet) {
006b5f7f 740 const AliITSdetector &det=layer.GetDetector(idet);
4ab260d6 741 ResetTrackToFollow(fTracks[fI]);
742 if (!fTrackToFollow.Propagate(det.GetPhi(),det.GetR())) {
8676d691 743 //Warning("TakeNextProlongation","propagation failed !\n");
006b5f7f 744 continue;
745 }
4ab260d6 746 fTrackToFollow.SetDetectorIndex(idet);
747 if (TMath::Abs(fTrackToFollow.GetZ()-GetZ())>layer.GetR()+dz) continue;
006b5f7f 748 }
749
4ab260d6 750 if (TMath::Abs(fTrackToFollow.GetZ() - c->GetZ()) > dz) continue;
751 if (TMath::Abs(fTrackToFollow.GetY() - c->GetY()) > dy) continue;
006b5f7f 752
4ab260d6 753 chi2=fTrackToFollow.GetPredictedChi2(c); if (chi2<kMaxChi2) break;
006b5f7f 754 }
755
006b5f7f 756 if (chi2>=kMaxChi2) return 0;
757 if (!c) return 0;
babd135a 758 if (c->IsUsed()&&c->GetNy()<5) { //shared factor
759 chi2+=1;
760 chi2*=2*(1./(TMath::Max(c->GetNy(),1)));
761 }
762 if (c->GetQ()==0){ //dead zone virtual cluster
763 chi2*=4.;
764 chi2+=20;
765 fTrackToFollow.fNUsed++;
766 return 1;
767 }
768 //if ((fI<2)&&chi2>kMaxChi2In) return 0;
006b5f7f 769
006b5f7f 770 if (!fTrackToFollow.Update(c,chi2,(fI<<28)+ci)) {
8676d691 771 //Warning("TakeNextProlongation","filtering failed !\n");
006b5f7f 772 return 0;
773 }
babd135a 774 if (c->IsUsed()&&c->GetNy()<5) fTrackToFollow.fNUsed++;
4ab260d6 775 if (fTrackToFollow.GetNumberOfClusters()>1)
776 if (TMath::Abs(fTrackToFollow.GetD())>4) return 0;
777
a9a2d814 778 fTrackToFollow.
779 SetSampledEdx(c->GetQ(),fTrackToFollow.GetNumberOfClusters()-1); //b.b.
780
781 {
880f41b9 782 Double_t x0;
61ab8ea8 783 Double_t d=layer.GetThickness(fTrackToFollow.GetY(),fTrackToFollow.GetZ(),x0);
880f41b9 784 fTrackToFollow.CorrectForMaterial(d,x0);
a9a2d814 785 }
4ab260d6 786
a9a2d814 787 if (fConstraint[fPass]) {
4ab260d6 788 Double_t d=GetEffectiveThickness(0,0); //Think of this !!!!
8676d691 789 Double_t xyz[]={GetX(),GetY(),GetZ()};
790 Double_t ers[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
babd135a 791 Double_t deltad = TMath::Abs(fTrackToFollow.GetD(GetX(),GetY()));
792 Double_t deltaz = TMath::Abs(fTrackToFollow.GetZat(GetX())-GetZ());
793
794 if ( (fI==4) && (deltad>2.0 || deltaz>1.5)) return 0; // don't "improve" secondaries
795 if ( (fI==3) && (deltad>1.5 || deltaz>0.9)) return 0; // don't "improve" secondaries
796 if ( (fI==2) && (deltad>0.9 || deltaz>0.6)) return 1; // don't "improve" secondaries
797 if ( (fI==1) && (deltad>0.3 || deltaz>0.3)) return 1; // don't "improve" secondaries
798 if ( (fI==0) && (deltad>0.1 || deltaz>0.1)) return 1; // don't "improve" secondaries
799
8676d691 800 fTrackToFollow.Improve(d,xyz,ers);
a9a2d814 801 }
006b5f7f 802
006b5f7f 803 return 1;
804}
805
806
006b5f7f 807AliITStrackerV2::AliITSlayer::AliITSlayer() {
808 //--------------------------------------------------------------------
809 //default AliITSlayer constructor
810 //--------------------------------------------------------------------
811 fN=0;
812 fDetectors=0;
813}
814
815AliITStrackerV2::AliITSlayer::
816AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd) {
817 //--------------------------------------------------------------------
818 //main AliITSlayer constructor
819 //--------------------------------------------------------------------
820 fR=r; fPhiOffset=p; fZOffset=z;
821 fNladders=nl; fNdetectors=nd;
822 fDetectors=new AliITSdetector[fNladders*fNdetectors];
823
824 fN=0;
825 fI=0;
b87bd7cd 826
827 fRoad=2*fR*TMath::Sqrt(3.14/1.);//assuming that there's only one cluster
006b5f7f 828}
829
830AliITStrackerV2::AliITSlayer::~AliITSlayer() {
831 //--------------------------------------------------------------------
832 // AliITSlayer destructor
833 //--------------------------------------------------------------------
834 delete[] fDetectors;
835 for (Int_t i=0; i<fN; i++) delete fClusters[i];
babd135a 836 for (Int_t i=0; i<kMaxClusterPerLayer;i++) fClusterWeight[i]=0;
006b5f7f 837}
838
61ab8ea8 839void AliITStrackerV2::AliITSlayer::ResetClusters() {
840 //--------------------------------------------------------------------
841 // This function removes loaded clusters
842 //--------------------------------------------------------------------
843 for (Int_t i=0; i<fN; i++) delete fClusters[i];
babd135a 844 for (Int_t i=0; i<kMaxClusterPerLayer;i++) fClusterWeight[i]=0;
61ab8ea8 845 fN=0;
846 fI=0;
847}
848
babd135a 849void AliITStrackerV2::AliITSlayer::ResetWeights() {
850 //--------------------------------------------------------------------
851 // This function reset weights of the clusters
852 //--------------------------------------------------------------------
853 for (Int_t i=0; i<kMaxClusterPerLayer;i++) fClusterWeight[i]=0;
854}
855
b87bd7cd 856void AliITStrackerV2::AliITSlayer::ResetRoad() {
857 //--------------------------------------------------------------------
858 // This function calculates the road defined by the cluster density
859 //--------------------------------------------------------------------
860 Int_t n=0;
861 for (Int_t i=0; i<fN; i++) {
862 if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
863 }
864 if (n>1) fRoad=2*fR*TMath::Sqrt(3.14/n);
865}
866
006b5f7f 867Int_t AliITStrackerV2::AliITSlayer::InsertCluster(AliITSclusterV2 *c) {
868 //--------------------------------------------------------------------
869 //This function adds a cluster to this layer
870 //--------------------------------------------------------------------
871 if (fN==kMaxClusterPerLayer) {
8676d691 872 ::Error("InsertCluster","Too many clusters !\n");
873 return 1;
006b5f7f 874 }
875
876 if (fN==0) {fClusters[fN++]=c; return 0;}
877 Int_t i=FindClusterIndex(c->GetZ());
878 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliITSclusterV2*));
879 fClusters[i]=c; fN++;
880
881 return 0;
882}
883
884Int_t AliITStrackerV2::AliITSlayer::FindClusterIndex(Double_t z) const {
885 //--------------------------------------------------------------------
886 // This function returns the index of the nearest cluster
887 //--------------------------------------------------------------------
888 if (fN==0) return 0;
889 if (z <= fClusters[0]->GetZ()) return 0;
890 if (z > fClusters[fN-1]->GetZ()) return fN;
891 Int_t b=0, e=fN-1, m=(b+e)/2;
892 for (; b<e; m=(b+e)/2) {
893 if (z > fClusters[m]->GetZ()) b=m+1;
894 else e=m;
895 }
896 return m;
897}
898
899void AliITStrackerV2::AliITSlayer::
900SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
901 //--------------------------------------------------------------------
902 // This function sets the "window"
903 //--------------------------------------------------------------------
904 fI=FindClusterIndex(zmin); fZmax=zmax;
14825d5a 905 Double_t circle=2*TMath::Pi()*fR;
906 if (ymax>circle) { ymax-=circle; ymin-=circle; }
006b5f7f 907 fYmin=ymin; fYmax=ymax;
908}
909
910const AliITSclusterV2 *AliITStrackerV2::AliITSlayer::GetNextCluster(Int_t &ci){
911 //--------------------------------------------------------------------
912 // This function returns clusters within the "window"
913 //--------------------------------------------------------------------
914 const AliITSclusterV2 *cluster=0;
915 for (Int_t i=fI; i<fN; i++) {
916 const AliITSclusterV2 *c=fClusters[i];
917 if (c->GetZ() > fZmax) break;
babd135a 918 // if (c->IsUsed()) continue;
006b5f7f 919 const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
920 Double_t y=fR*det.GetPhi() + c->GetY();
921
922 if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
923 if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
924
925 if (y<fYmin) continue;
926 if (y>fYmax) continue;
927 cluster=c; ci=i;
928 fI=i+1;
929 break;
930 }
7f6ddf58 931
006b5f7f 932 return cluster;
933}
934
935Int_t AliITStrackerV2::AliITSlayer::
936FindDetectorIndex(Double_t phi, Double_t z) const {
937 //--------------------------------------------------------------------
938 //This function finds the detector crossed by the track
939 //--------------------------------------------------------------------
f363d377 940 Double_t dphi=-(phi-fPhiOffset);
006b5f7f 941 if (dphi < 0) dphi += 2*TMath::Pi();
942 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
943 Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
3e3e8427 944 if (np>=fNladders) np-=fNladders;
945 if (np<0) np+=fNladders;
006b5f7f 946
947 Double_t dz=fZOffset-z;
948 Int_t nz=Int_t(dz*(fNdetectors-1)*0.5/fZOffset+0.5);
3e3e8427 949 if (nz>=fNdetectors) return -1;
950 if (nz<0) return -1;
006b5f7f 951
006b5f7f 952 return np*fNdetectors + nz;
953}
954
955Double_t
61ab8ea8 956AliITStrackerV2::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
957const {
a9a2d814 958 //--------------------------------------------------------------------
959 //This function returns the layer thickness at this point (units X0)
960 //--------------------------------------------------------------------
961 Double_t d=0.0085;
61ab8ea8 962 x0=21.82;
a9a2d814 963
964 if (43<fR&&fR<45) { //SSD2
61ab8ea8 965 Double_t dd=0.0034;
966 d=dd;
967 if (TMath::Abs(y-0.00)>3.40) d+=dd;
968 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
969 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
a9a2d814 970 for (Int_t i=0; i<12; i++) {
61ab8ea8 971 if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
972 if (TMath::Abs(y-0.00)>3.40) d+=dd;
973 d+=0.0034;
974 break;
975 }
976 if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
977 if (TMath::Abs(y-0.00)>3.40) d+=dd;
978 d+=0.0034;
979 break;
980 }
981 if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
982 if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
a9a2d814 983 }
984 } else
985 if (37<fR&&fR<41) { //SSD1
61ab8ea8 986 Double_t dd=0.0034;
987 d=dd;
988 if (TMath::Abs(y-0.00)>3.40) d+=dd;
989 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
990 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
a9a2d814 991 for (Int_t i=0; i<11; i++) {
61ab8ea8 992 if (TMath::Abs(z-3.9*i)<0.15) {
993 if (TMath::Abs(y-0.00)>3.40) d+=dd;
994 d+=dd;
995 break;
996 }
997 if (TMath::Abs(z+3.9*i)<0.15) {
998 if (TMath::Abs(y-0.00)>3.40) d+=dd;
999 d+=dd;
1000 break;
1001 }
1002 if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1003 if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
a9a2d814 1004 }
1005 } else
1006 if (13<fR&&fR<26) { //SDD
61ab8ea8 1007 Double_t dd=0.0033;
1008 d=dd;
1009 if (TMath::Abs(y-0.00)>3.30) d+=dd;
1010
1011 if (TMath::Abs(y-1.80)<0.55) {
1012 d+=0.016;
1013 for (Int_t j=0; j<20; j++) {
1014 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1015 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1016 }
1017 }
1018 if (TMath::Abs(y+1.80)<0.55) {
1019 d+=0.016;
1020 for (Int_t j=0; j<20; j++) {
1021 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1022 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1023 }
1024 }
1025
1026 for (Int_t i=0; i<4; i++) {
1027 if (TMath::Abs(z-7.3*i)<0.60) {
1028 d+=dd;
1029 if (TMath::Abs(y-0.00)>3.30) d+=dd;
1030 break;
1031 }
1032 if (TMath::Abs(z+7.3*i)<0.60) {
1033 d+=dd;
1034 if (TMath::Abs(y-0.00)>3.30) d+=dd;
1035 break;
1036 }
a9a2d814 1037 }
1038 } else
1039 if (6<fR&&fR<8) { //SPD2
61ab8ea8 1040 Double_t dd=0.0063; x0=21.5;
1041 d=dd;
1042 if (TMath::Abs(y-3.08)>0.5) d+=dd;
1043 //if (TMath::Abs(y-3.08)>0.45) d+=dd;
1044 if (TMath::Abs(y-3.03)<0.10) {d+=0.014;}
a9a2d814 1045 } else
1046 if (3<fR&&fR<5) { //SPD1
61ab8ea8 1047 Double_t dd=0.0063; x0=21.5;
1048 d=dd;
1049 if (TMath::Abs(y+0.21)>0.6) d+=dd;
1050 //if (TMath::Abs(y+0.21)>0.45) d+=dd;
1051 if (TMath::Abs(y+0.10)<0.10) {d+=0.014;}
a9a2d814 1052 }
1053
a9a2d814 1054 return d;
1055}
006b5f7f 1056
a9a2d814 1057Double_t AliITStrackerV2::GetEffectiveThickness(Double_t y,Double_t z) const
006b5f7f 1058{
1059 //--------------------------------------------------------------------
a9a2d814 1060 //Returns the thickness between the current layer and the vertex (units X0)
006b5f7f 1061 //--------------------------------------------------------------------
61ab8ea8 1062 Double_t d=0.0028*3*3; //beam pipe
1063 Double_t x0=0;
006b5f7f 1064
18856a77 1065 Double_t xn=fgLayers[fI].GetR();
006b5f7f 1066 for (Int_t i=0; i<fI; i++) {
18856a77 1067 Double_t xi=fgLayers[i].GetR();
1068 d+=fgLayers[i].GetThickness(y,z,x0)*xi*xi;
006b5f7f 1069 }
1070
1071 if (fI>1) {
61ab8ea8 1072 Double_t xi=9.;
1073 d+=0.0097*xi*xi;
006b5f7f 1074 }
1075
1076 if (fI>3) {
18856a77 1077 Double_t xi=0.5*(fgLayers[3].GetR()+fgLayers[4].GetR());
61ab8ea8 1078 d+=0.0034*xi*xi;
006b5f7f 1079 }
61ab8ea8 1080
006b5f7f 1081 return d/(xn*xn);
1082}
1083
006b5f7f 1084Int_t AliITStrackerV2::AliITSlayer::InRoad() const {
1085 //--------------------------------------------------------------------
1086 // This function returns number of clusters within the "window"
1087 //--------------------------------------------------------------------
1088 Int_t ncl=0;
1089 for (Int_t i=fI; i<fN; i++) {
1090 const AliITSclusterV2 *c=fClusters[i];
1091 if (c->GetZ() > fZmax) break;
a9a2d814 1092 if (c->IsUsed()) continue;
006b5f7f 1093 const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
1094 Double_t y=fR*det.GetPhi() + c->GetY();
1095
1096 if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
1097 if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
1098
1099 if (y<fYmin) continue;
1100 if (y>fYmax) continue;
1101 ncl++;
1102 }
1103 return ncl;
1104}
1105
880f41b9 1106Bool_t
c630aafd 1107AliITStrackerV2::RefitAt(Double_t xx,AliITStrackV2 *t,const AliITStrackV2 *c) {
4ab260d6 1108 //--------------------------------------------------------------------
c630aafd 1109 // This function refits the track "t" at the position "x" using
1110 // the clusters from "c"
4ab260d6 1111 //--------------------------------------------------------------------
880f41b9 1112 Int_t index[kMaxLayer];
1113 Int_t k;
1114 for (k=0; k<kMaxLayer; k++) index[k]=-1;
c630aafd 1115 Int_t nc=c->GetNumberOfClusters();
880f41b9 1116 for (k=0; k<nc; k++) {
c630aafd 1117 Int_t idx=c->GetClusterIndex(k),nl=(idx&0xf0000000)>>28;
880f41b9 1118 index[nl]=idx;
1119 }
880f41b9 1120
4ab260d6 1121 Int_t from, to, step;
ee34184f 1122 if (xx > t->GetX()) {
4ab260d6 1123 from=0; to=kMaxLayer;
1124 step=+1;
1125 } else {
1126 from=kMaxLayer-1; to=-1;
1127 step=-1;
1128 }
1129
1130 for (Int_t i=from; i != to; i += step) {
18856a77 1131 AliITSlayer &layer=fgLayers[i];
4ab260d6 1132 Double_t r=layer.GetR();
1133
1134 {
1135 Double_t hI=i-0.5*step;
880f41b9 1136 if (TMath::Abs(hI-1.5)<0.01 || TMath::Abs(hI-3.5)<0.01) {
18856a77 1137 Double_t rs=0.5*(fgLayers[i-step].GetR() + r);
61ab8ea8 1138 Double_t d=0.0034, x0=38.6;
880f41b9 1139 if (TMath::Abs(hI-1.5)<0.01) {rs=9.; d=0.0097; x0=42;}
1140 if (!t->PropagateTo(rs,-step*d,x0)) {
4ab260d6 1141 return kFALSE;
1142 }
1143 }
1144 }
1145
880f41b9 1146 // remember old position [SR, GSI 18.02.2003]
1147 Double_t oldX=0., oldY=0., oldZ=0.;
1148 if (t->IsStartedTimeIntegral() && step==1) {
1149 t->GetGlobalXYZat(t->GetX(),oldX,oldY,oldZ);
1150 }
1151 //
1152
4ab260d6 1153 Double_t x,y,z;
1154 if (!t->GetGlobalXYZat(r,x,y,z)) {
1155 return kFALSE;
1156 }
1157 Double_t phi=TMath::ATan2(y,x);
1158 Int_t idet=layer.FindDetectorIndex(phi,z);
1159 if (idet<0) {
1160 return kFALSE;
1161 }
1162 const AliITSdetector &det=layer.GetDetector(idet);
1163 phi=det.GetPhi();
1164 if (!t->Propagate(phi,det.GetR())) {
1165 return kFALSE;
1166 }
1167 t->SetDetectorIndex(idet);
1168
1169 const AliITSclusterV2 *cl=0;
1170 Double_t maxchi2=kMaxChi2;
1171
1172 Int_t idx=index[i];
1173 if (idx>0) {
1174 const AliITSclusterV2 *c=(AliITSclusterV2 *)GetCluster(idx);
1175 if (idet != c->GetDetectorIndex()) {
1176 idet=c->GetDetectorIndex();
1177 const AliITSdetector &det=layer.GetDetector(idet);
1178 if (!t->Propagate(det.GetPhi(),det.GetR())) {
1179 return kFALSE;
1180 }
1181 t->SetDetectorIndex(idet);
1182 }
1183 Double_t chi2=t->GetPredictedChi2(c);
c630aafd 1184 if (chi2<maxchi2) {
1185 cl=c;
1186 maxchi2=chi2;
1187 } else {
1188 return kFALSE;
1189 }
4ab260d6 1190 }
4ab260d6 1191 /*
1192 if (cl==0)
1193 if (t->GetNumberOfClusters()>2) {
1194 Double_t dz=4*TMath::Sqrt(t->GetSigmaZ2()+kSigmaZ2[i]);
1195 Double_t dy=4*TMath::Sqrt(t->GetSigmaY2()+kSigmaY2[i]);
1196 Double_t zmin=t->GetZ() - dz;
1197 Double_t zmax=t->GetZ() + dz;
1198 Double_t ymin=t->GetY() + phi*r - dy;
1199 Double_t ymax=t->GetY() + phi*r + dy;
1200 layer.SelectClusters(zmin,zmax,ymin,ymax);
1201
1202 const AliITSclusterV2 *c=0; Int_t ci=-1;
1203 while ((c=layer.GetNextCluster(ci))!=0) {
1204 if (idet != c->GetDetectorIndex()) continue;
1205 Double_t chi2=t->GetPredictedChi2(c);
1206 if (chi2<maxchi2) { cl=c; maxchi2=chi2; idx=ci; }
1207 }
1208 }
1209 */
4ab260d6 1210 if (cl) {
1211 if (!t->Update(cl,maxchi2,idx)) {
1212 return kFALSE;
1213 }
880f41b9 1214 t->SetSampledEdx(cl->GetQ(),t->GetNumberOfClusters()-1);
4ab260d6 1215 }
1216
1217 {
61ab8ea8 1218 Double_t x0;
1219 Double_t d=layer.GetThickness(t->GetY(),t->GetZ(),x0);
1220 t->CorrectForMaterial(-step*d,x0);
4ab260d6 1221 }
880f41b9 1222
1223 // track time update [SR, GSI 17.02.2003]
1224 if (t->IsStartedTimeIntegral() && step==1) {
1225 Double_t newX, newY, newZ;
1226 t->GetGlobalXYZat(t->GetX(),newX,newY,newZ);
1227 Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) +
1228 (oldZ-newZ)*(oldZ-newZ);
1229 t->AddTimeStep(TMath::Sqrt(dL2));
1230 }
1231 //
4ab260d6 1232
1233 }
1234
ee34184f 1235 if (!t->PropagateTo(xx,0.,0.)) return kFALSE;
4ab260d6 1236 return kTRUE;
1237}
1238
babd135a 1239
1240Float_t *AliITStrackerV2::GetWeight(Int_t index) {
1241 //--------------------------------------------------------------------
1242 // Return pointer to a given cluster
1243 //--------------------------------------------------------------------
1244 Int_t l=(index & 0xf0000000) >> 28;
1245 Int_t c=(index & 0x0fffffff) >> 00;
1246 return fgLayers[l].GetWeight(c);
1247}
1248
1249
1250
b87bd7cd 1251void AliITStrackerV2::UseClusters(const AliKalmanTrack *t, Int_t from) const {
1252 //--------------------------------------------------------------------
1253 // This function marks clusters assigned to the track
1254 //--------------------------------------------------------------------
1255 AliTracker::UseClusters(t,from);
880f41b9 1256
b87bd7cd 1257 AliITSclusterV2 *c=(AliITSclusterV2 *)GetCluster(t->GetClusterIndex(0));
1258 //if (c->GetQ()>2) c->Use();
1259 if (c->GetSigmaZ2()>0.1) c->Use();
1260 c=(AliITSclusterV2 *)GetCluster(t->GetClusterIndex(1));
1261 //if (c->GetQ()>2) c->Use();
1262 if (c->GetSigmaZ2()>0.1) c->Use();
880f41b9 1263
b87bd7cd 1264}
de5d3595 1265
1266
1267void AliITStrackerV2::AddTrackHypothesys(AliITStrackV2 * track, Int_t esdindex)
1268{
1269 //------------------------------------------------------------------
1270 // add track to the list of hypothesys
1271 //------------------------------------------------------------------
1272
1273 if (esdindex>=fTrackHypothesys.GetEntriesFast()) fTrackHypothesys.Expand(esdindex*2+10);
1274 //
2257f27e 1275 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
de5d3595 1276 if (!array) {
1277 array = new TObjArray(10);
1278 fTrackHypothesys.AddAt(array,esdindex);
1279 }
1280 array->AddLast(track);
1281}
1282
babd135a 1283void AliITStrackerV2::SortTrackHypothesys(Int_t esdindex, Float_t likelihoodlevel)
de5d3595 1284{
1285 //-------------------------------------------------------------------
1286 // compress array of track hypothesys
1287 // keep only maxsize best hypothesys
1288 //-------------------------------------------------------------------
2257f27e 1289 if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
1290 if (! (fTrackHypothesys.At(esdindex)) ) return;
1291 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
babd135a 1292 Int_t entries = array->GetEntriesFast();
1293
1294 Float_t * chi2 = new Float_t[entries];
1295 Float_t * probability = new Float_t[entries];
1296 Float_t sumprobabilityall=0;
1297 Int_t * index = new Int_t[entries];
1298 //
de5d3595 1299 //
babd135a 1300 for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
1301 AliITStrackV2 * track = (AliITStrackV2*)array->At(itrack);
1302 //
1303 if (track && track->GetNumberOfClusters()>(track->fNUsed+track->fNSkipped)){
1304 //
1305 chi2[itrack] = track->GetChi2()/(track->GetNumberOfClusters()-track->fNUsed-track->fNSkipped);
1306 if (track->fESDtrack)
1307 if (track->fESDtrack->GetTPCsignal()>80){
1308 track->CookdEdx();
1309 if ((track->GetdEdx()/track->fESDtrack->GetTPCsignal())<0.4){
1310 Float_t factor = 2.+10.*(0.6-track->GetdEdx()/track->fESDtrack->GetTPCsignal());
1311 chi2[itrack]*= factor; //mismatch in dEdx
1312 }
1313 }
1314 }
2257f27e 1315 else
babd135a 1316 chi2[itrack] = 10000000.;
1317 probability[itrack] = 1./(0.3+chi2[itrack]);
1318 sumprobabilityall+=probability[itrack];
de5d3595 1319 }
babd135a 1320
1321 TMath::Sort(entries,chi2,index,kFALSE);
de5d3595 1322 TObjArray * newarray = new TObjArray();
babd135a 1323 Float_t sumprobability = 0.;
1324 for (Int_t i=0;i<entries;i++){
1325 AliITStrackV2 * track = (AliITStrackV2*)array->At(index[i]);
1326 if (!track) break;
1327 if (chi2[index[i]]<30){
1328 newarray->AddLast(array->RemoveAt(index[i]));
1329 track->fChi2MIP[0] = chi2[index[i]]; // base chi 2
1330 sumprobability+= probability[index[i]]/sumprobabilityall;
1331 if (sumprobability> likelihoodlevel) break;
1332 }
1333 else{
1334 delete array->RemoveAt(index[i]);
1335 }
de5d3595 1336 }
babd135a 1337
de5d3595 1338 array->Delete();
2257f27e 1339 delete fTrackHypothesys.RemoveAt(esdindex);
de5d3595 1340 fTrackHypothesys.AddAt(newarray,esdindex);
1341
2257f27e 1342 delete [] chi2;
1343 delete [] index;
1344
de5d3595 1345}
1346
1347
babd135a 1348void AliITStrackerV2::CompressTrackHypothesys(Int_t esdindex, Float_t likelihoodlevel, Int_t maxsize)
1349{
1350 //
1351 //
1352 if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
1353 if (! (fTrackHypothesys.At(esdindex)) ) return;
1354 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
1355 Int_t entries = array->GetEntriesFast();
1356 //
1357 if (likelihoodlevel>0.000001){
1358 Float_t *probability = new Float_t[entries];
1359 for (Int_t i=0;i<entries;i++) probability[i]=0.;
1360 Float_t sumprobabilityall=0;
1361 Int_t *index = new Int_t[entries];
1362
1363 for (Int_t itrack=0;itrack<entries;itrack++){
1364 AliITStrackV2 * track = (AliITStrackV2*)array->At(itrack);
1365 probability[itrack]=0;
1366 if (!track) continue;
1367 probability[itrack] = 1./(0.3+track->fChi2MIP[0]);
1368 sumprobabilityall += probability[itrack];
1369 //
1370 }
1371 if (sumprobabilityall<=0.000000000000001){
1372 return;
1373 }
1374 for (Int_t itrack=0;itrack<entries;itrack++){
1375 probability[itrack]/=sumprobabilityall;
1376 }
1377
1378 TMath::Sort(entries, probability, index, kTRUE);
1379 //
1380 TObjArray * newarray = new TObjArray();
1381 Float_t sumprobability = 0.;
1382 for (Int_t i=0;i<entries;i++){
1383 AliITStrackV2 * track = (AliITStrackV2*)array->At(index[i]);
1384 if (!track) continue;
1385 newarray->AddLast(array->RemoveAt(index[i]));
1386 sumprobability+= probability[index[i]];
1387 if (sumprobability> likelihoodlevel) break;
1388 if (i>maxsize) break;
1389 }
1390
1391 array->Delete();
1392 delete fTrackHypothesys.RemoveAt(esdindex);
1393 fTrackHypothesys.AddAt(newarray,esdindex);
1394 //
1395 delete []index;
1396 delete []probability;
1397 }
1398 else{
1399 array->Delete();
1400 delete fTrackHypothesys.RemoveAt(esdindex);
1401 }
1402}
1403
de5d3595 1404
babd135a 1405AliITStrackV2 * AliITStrackerV2::GetBestHypothesys(Int_t esdindex, AliITStrackV2 * original, Int_t checkmax)
de5d3595 1406{
1407 //-------------------------------------------------------------
1408 // try to find best hypothesy
1409 // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
1410 //-------------------------------------------------------------
babd135a 1411 if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
de5d3595 1412 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
1413 if (!array) return 0;
1414 Int_t entries = array->GetEntriesFast();
1415 if (!entries) return 0;
1416 Float_t minchi2 = 100000;
babd135a 1417 Int_t maxn = 3;
de5d3595 1418 AliITStrackV2 * besttrack=0;
babd135a 1419 Int_t accepted =0;
1420 Int_t maxindex=0;
1421 //
1422 Float_t sumz2=0;
1423 Float_t sumy2=0;
1424 Float_t sumall=0;
1425 for (Int_t itrack=0;itrack<entries; itrack++){
1426 AliITStrackV2 * track = (AliITStrackV2*)array->At(itrack);
1427 if (!track) continue;
1428 sumall++;
1429 sumz2+=track->GetSigmaZ2();
1430 sumy2+=track->GetSigmaY2();
1431 }
1432 sumz2/=sumall;
1433 sumy2/=sumall;
1434
1435 Float_t dedxmismatch=1;
1436 for (Int_t i=0;i<entries;i++){
1437 maxindex = i;
1438 AliITStrackV2 * track = (AliITStrackV2*)array->At(i);
1439 if (!track) continue;
1440 track->fChi2MIP[1] = 1000000;
1441 track->fChi2MIP[2] = 1000000;
1442
1443 if ( (track->GetNumberOfClusters()-track->fNSkipped-track->fNUsed)<2) continue;
1444 //
1445 if (track->fESDtrack)
1446 if (track->fESDtrack->GetTPCsignal()>80){
1447 track->CookdEdx();
1448 if ((track->GetdEdx()/track->fESDtrack->GetTPCsignal())<0.4){
1449 //mismatch in dEdx
1450 dedxmismatch= 2.+ 10.*(0.6-track->GetdEdx()/track->fESDtrack->GetTPCsignal());
1451 }
1452 }
1453
1454 // track->SetLabel(original->GetLabel());
1455 //CookLabel(track,0.0);
1456 //if (track->GetFakeRatio()>0.01) continue;
1457 //
1458 //
1459 // backtrack
1460 AliITStrackV2 * backtrack = new AliITStrackV2(*track);
1461 backtrack->ResetCovariance();
1462 backtrack->ResetClusters();
1463 Double_t x = original->GetX();
1464 if (!RefitAt(x,backtrack,track)){
1465 delete backtrack;
1466 delete array->RemoveAt(i);
1467 continue;
1468 }
1469 if ( (backtrack->GetChi2() / float(backtrack->GetNumberOfClusters()-track->fNSkipped-track->fNUsed-0.5))>6)
1470 {
1471 delete backtrack;
1472 delete array->RemoveAt(i);
1473 continue;
1474 }
1475 Double_t deltac = backtrack->GetC()-original->GetC();
1476 Double_t deltatgl = backtrack->GetTgl()-original->GetTgl();
1477 //
1478 Double_t poolc2 = (deltac*deltac)/(original->fC44+backtrack->fC44);
1479 Double_t pooltgl2 = (deltatgl*deltatgl)/(original->fC33+backtrack->fC33);
1480 if ((poolc2+pooltgl2)>32){ //4 sigma
1481 delete backtrack;
1482 delete array->RemoveAt(i);
1483 continue;
1484 }
1485 //Double_t bpoolc = (deltac*deltac)/(original->fC44);
1486 //Double_t bpooltgl = (deltatgl*deltatgl)/(original->fC33);
1487
1488 //
1489 //forward track - without constraint
1490 AliITStrackV2 * forwardtrack = new AliITStrackV2(*original);
1491 // forwardtrack->ResetCovariance();
1492 forwardtrack->ResetClusters();
1493 x = track->GetX();
1494 if (!RefitAt(x,forwardtrack,track)){
1495 delete forwardtrack;
1496 delete backtrack;
1497 delete array->RemoveAt(i);
1498 continue;
1499 }
1500 if ( (forwardtrack->GetChi2()/float(forwardtrack->GetNumberOfClusters()-track->fNSkipped-track->fNUsed))>6)
1501 {
1502 delete forwardtrack;
1503 delete array->RemoveAt(i);
1504 continue;
1505 }
1506 //
1507 accepted++;
1508 if (accepted>checkmax){
1509 delete backtrack;
1510 delete forwardtrack;
1511 break;
1512 }
1513 Double_t chi2 = (backtrack->GetChi2()/(backtrack->GetNumberOfClusters()-1-track->fNSkipped-track->fNUsed)+
1514 forwardtrack->GetChi2()/(forwardtrack->GetNumberOfClusters()-track->fNSkipped-track->fNUsed));
1515 // bpoolc+bpooltgl;
1516 // chi2 *= (forwardtrack->GetSigmaZ2()/sumz2+forwardtrack->GetSigmaY2()/sumy2);
1517 chi2 *= dedxmismatch;
1518 //
1519 //
1520 track->fChi2MIP[1] = backtrack->GetChi2()/(backtrack->GetNumberOfClusters()-1-track->fNSkipped-track->fNUsed);
1521 track->fChi2MIP[2] = forwardtrack->GetChi2()/(forwardtrack->GetNumberOfClusters()-track->fNSkipped-track->fNUsed);
1522 track->fChi2MIP[3] = poolc2+pooltgl2;
1523 //
1524
1525 if (track->GetNumberOfClusters()>maxn){
1526 besttrack = new AliITStrackV2(*forwardtrack);
1527 maxn = track->GetNumberOfClusters();
1528 minchi2 = chi2;
1529 delete backtrack;
1530 delete forwardtrack;
1531 continue;
1532 }
1533 //
1534 if (chi2 < minchi2){
1535 besttrack = new AliITStrackV2(*forwardtrack);
1536 minchi2 = chi2;
1537 }
1538 delete backtrack;
1539 delete forwardtrack;
1540 }
1541 //
1542 //
1543 if (!besttrack || besttrack->GetNumberOfClusters()<4) {
1544 return 0;
1545 }
1546
1547 //
1548 besttrack->SetLabel(original->GetLabel());
1549 CookLabel(besttrack,0.0);
de5d3595 1550 //
1551 // calculate "weight of the cluster"
1552 //
babd135a 1553
1554 {
1555 //sign usage information for clusters
1556 Int_t clusterindex[6][100];
1557 Double_t clusterweight[6][100];
1558 for (Int_t ilayer=0;ilayer<6;ilayer++)
1559 for (Int_t icluster=0;icluster<100;icluster++){
1560 clusterindex[ilayer][icluster] = -1;
1561 clusterweight[ilayer][icluster] = 0;
1562 }
1563 //printf("%d\t%d\n",esdindex, entries);
1564 //
1565 Float_t sumchi2=0;
1566 for (Int_t itrack=0;itrack<entries; itrack++){
1567 AliITStrackV2 * track = (AliITStrackV2*)array->At(itrack);
1568 if (!track) continue;
1569 if (track->fChi2MIP[1]>1000) continue; //not accepted
1570 sumchi2 +=1./(0.3+track->fChi2MIP[1]+track->fChi2MIP[2]);
1571 }
1572 for (Int_t itrack=0;itrack<entries;itrack++){
1573 AliITStrackV2 * track = (AliITStrackV2*)array->At(itrack);
1574 if (!track) continue;
1575 if (track->fChi2MIP[1]>1000) continue; //not accepted
1576 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
1577 Int_t tindex = track->GetClusterIndex(icluster);
1578 Int_t ilayer = (tindex & 0xf0000000) >> 28;
1579 if (tindex<0) continue;
1580 Int_t cindex =0;
1581 //
1582 for (cindex=0;cindex<100;cindex++){
1583 if (clusterindex[ilayer][cindex]<0) break;
1584 if (clusterindex[ilayer][cindex]==tindex) break;
1585 }
1586 if (cindex>100) break;
1587 if (clusterindex[ilayer][cindex]!=tindex) clusterindex[ilayer][cindex] = tindex;
1588 clusterweight[ilayer][cindex]+= (1./(0.3+track->fChi2MIP[1]+track->fChi2MIP[2]))* (1./sumchi2);
1589 Float_t *weight = GetWeight(tindex);
1590
1591 if (weight){
1592 *weight+= (1./(0.3+track->fChi2MIP[1]+track->fChi2MIP[2]))* (1./sumchi2);
1593 }
1594 }
1595 }
1596
1597 if (besttrack->GetChi2()/besttrack->GetNumberOfClusters()>3.5) return besttrack; //don't sign clusters
1598 Int_t current=0;
1599 Double_t deltad = besttrack->GetD(GetX(),GetY());
1600 Double_t deltaz = besttrack->GetZat(GetX()) - GetZ();
1601 Double_t deltaprim = TMath::Sqrt(deltad*deltad+deltaz*deltaz);
1602
1603 for (Int_t icluster=0;icluster<besttrack->GetNumberOfClusters();icluster++){
1604 Int_t index = besttrack->GetClusterIndex(icluster);
1605 Int_t ilayer = (index & 0xf0000000) >> 28;
1606 AliITSclusterV2 *c = (AliITSclusterV2*)GetCluster(index);
1607 if (!c) continue;
1608 //
1609 for (Int_t icluster=0;icluster<100;icluster++){
1610 // sign non "doubt" clusters as used
1611 if (clusterindex[ilayer][icluster]!=index) continue;
1612 // Float_t * weight = GetWeight(index);
1613 //if (weight) if (*weight>1){
1614 // if (c->IsUsed()) continue;
1615 // c->Use();
1616 //}
1617 if ( (ilayer*0.2+0.2)<deltaprim) continue; // secondaries
1618 if (c->GetNy()>4) continue; // don sign cluster
1619 if ( (ilayer>1&&clusterweight[ilayer][icluster]>0.7) || (ilayer<2&&clusterweight[ilayer][icluster]>0.8) ){
1620 current++;
1621 if (c->IsUsed()) continue;
1622 c->Use();
1623 }
1624 }
1625 }
1626 }
1627
1628 //
1629 return besttrack;
1630}
1631
1632
1633AliITStrackV2 * AliITStrackerV2::GetBestHypothesysMIP(Int_t esdindex, AliITStrackV2 * original)
1634{
1635 //-------------------------------------------------------------
1636 // try to find best hypothesy
1637 // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
1638 //-------------------------------------------------------------
1639 if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
1640 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
1641 if (!array) return 0;
1642 Int_t entries = array->GetEntriesFast();
1643 if (!entries) return 0;
1644 AliITStrackV2 * besttrack=0;
1645 //
1646 //sign usage information for clusters
de5d3595 1647 Int_t clusterindex[6][100];
1648 Double_t clusterweight[6][100];
1649 for (Int_t ilayer=0;ilayer<6;ilayer++)
1650 for (Int_t icluster=0;icluster<100;icluster++){
1651 clusterindex[ilayer][icluster] = -1;
1652 clusterweight[ilayer][icluster] = 0;
1653 }
1654 //printf("%d\t%d\n",esdindex, entries);
1655 //
1656 Float_t sumchi2=0;
babd135a 1657 for (Int_t itrack=0;itrack<entries; itrack++){
de5d3595 1658 AliITStrackV2 * track = (AliITStrackV2*)array->At(itrack);
1659 if (!track) continue;
babd135a 1660 if (track->fChi2MIP[1]>1000) continue; //not accepted - before
1661 sumchi2 +=1./(0.3+track->fChi2MIP[1]+track->fChi2MIP[2]);
de5d3595 1662 }
babd135a 1663 //
1664 // get cluster weight
de5d3595 1665 for (Int_t itrack=0;itrack<entries;itrack++){
1666 AliITStrackV2 * track = (AliITStrackV2*)array->At(itrack);
1667 if (!track) continue;
babd135a 1668 if (track->fChi2MIP[1]>1000) continue; // track not accepted in previous iterration
de5d3595 1669 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
1670 Int_t tindex = track->GetClusterIndex(icluster);
1671 Int_t ilayer = (tindex & 0xf0000000) >> 28;
1672 if (tindex<0) continue;
1673 Int_t cindex =0;
1674 //
1675 for (cindex=0;cindex<100;cindex++){
1676 if (clusterindex[ilayer][cindex]<0) break;
1677 if (clusterindex[ilayer][cindex]==tindex) break;
1678 }
1679 if (cindex>100) break;
1680 if (clusterindex[ilayer][cindex]!=tindex) clusterindex[ilayer][cindex] = tindex;
babd135a 1681 clusterweight[ilayer][cindex]+= (1./(0.3+track->fChi2MIP[1]+track->fChi2MIP[2]))* (1./sumchi2);
de5d3595 1682 }
1683 }
babd135a 1684 //
1685 // get cluster relative sharing - factor
de5d3595 1686 //
babd135a 1687 //
1688 for (Int_t ilayer=0;ilayer<6;ilayer++)
1689 for (Int_t cindex=0;cindex<100;cindex++){
1690 if (clusterindex[ilayer][cindex]<0) continue;
1691 Int_t tindex = clusterindex[ilayer][cindex];
1692 Float_t *weight = GetWeight(tindex);
1693 if (!weight){
1694 printf("Problem 1\n"); // not existing track
1695 continue;
1696 }
1697 if (*weight<(clusterweight[ilayer][cindex]-0.00001)){
1698 printf("Problem 2\n"); // not normalized probability
1699 continue;
1700 }
1701 AliITSclusterV2 *c = (AliITSclusterV2*)GetCluster(tindex);
1702 if (c->GetNy()<5){
1703 clusterweight[ilayer][cindex]/= *weight;
1704 }
1705 else{
1706 Float_t weight2 = TMath::Max(*weight-0.5*clusterweight[ilayer][cindex],0.0000001);
1707 clusterweight[ilayer][cindex]/= weight2;
1708 if ( clusterweight[ilayer][cindex]>1) clusterweight[ilayer][cindex] =1.;
1709 }
de5d3595 1710 }
babd135a 1711 //
1712 //take to the account sharing factor
1713 //
1714 Float_t chi2 =10000000;
1715 Float_t sharefactor=0;
1716 Float_t minchi2 = 100000000;
1717 Float_t secchi2 = 100000000;
1718 Float_t norm=0;
1719 for (Int_t itrack=0;itrack<entries; itrack++){
1720 AliITStrackV2 * track = (AliITStrackV2*)array->At(itrack);
1721 if (!track) continue;
1722 if (track->fChi2MIP[1]>1000) continue; //not accepted - before
1723 chi2 = track->fChi2MIP[1];
1724 Float_t newchi2=0;
1725 sharefactor=0;
1726 norm =0;
1727 //
1728 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
1729 Int_t tindex = track->GetClusterIndex(icluster);
1730 Int_t ilayer = (tindex & 0xf0000000) >> 28;
1731 if (tindex<0) continue;
1732 Int_t cindex =0;
1733 Float_t cchi2 = (track->fDy[ilayer]*track->fDy[ilayer])/(track->fSigmaY[ilayer]*track->fSigmaY[ilayer]) +
1734 (track->fDz[ilayer]*track->fDz[ilayer])/(track->fSigmaZ[ilayer]*track->fSigmaZ[ilayer]) ;
1735 //
1736 for (cindex=0;cindex<100;cindex++){
1737 if (clusterindex[ilayer][cindex]<0) break;
1738 if (clusterindex[ilayer][cindex]==tindex) break;
1739 }
1740 if (cindex>100) continue;
1741 if (clusterweight[ilayer][cindex]>0.00001){
1742 sharefactor+= clusterweight[ilayer][cindex];
1743 cchi2/=clusterweight[ilayer][cindex];
1744 norm++;
1745 }
1746 newchi2+=cchi2;
de5d3595 1747 }
babd135a 1748 newchi2/=(norm-track->fNSkipped-track->fNUsed);
1749 track->fChi2MIP[4] = newchi2;
1750 //if (norm>0) sharefactor/=norm;
1751 //if (sharefactor<0.5) return 0;
1752 // chi2*=1./(0.5+sharefactor);
1753 if (newchi2<minchi2){
de5d3595 1754 besttrack = track;
babd135a 1755 minchi2 = newchi2;
1756 }
1757 else{
1758 if (newchi2<secchi2){
1759 secchi2 = newchi2;
1760 }
1761 }
1762 //
de5d3595 1763 }
babd135a 1764
1765 //
de5d3595 1766 //
babd135a 1767 if (!besttrack || besttrack->GetNumberOfClusters()<4) {
1768 return 0;
1769 }
de5d3595 1770 //
babd135a 1771 if ((minchi2/secchi2)<0.7){
1772 //
1773 //increase the weight for clusters if the probability of other hypothesys is small
1774 Double_t deltad = besttrack->GetD(GetX(),GetY());
1775 Double_t deltaz = besttrack->GetZat(GetX()) - GetZ();
1776 Double_t deltaprim = TMath::Sqrt(deltad*deltad+deltaz*deltaz);
1777 //
1778 for (Int_t icluster=0;icluster<besttrack->GetNumberOfClusters();icluster++){
1779 Int_t index = besttrack->GetClusterIndex(icluster);
1780 Int_t ilayer = (index & 0xf0000000) >> 28;
1781 AliITSclusterV2 *c = (AliITSclusterV2*)GetCluster(index);
1782 if (!c) continue;
1783 if (c->GetNy()>3) continue;
1784 if ( (ilayer*0.2+0.2)<deltaprim) continue; // secondaries
1785 Float_t * weight = GetWeight(index);
1786 *weight*=2.;
1787 *weight+=1.;
1788 }
1789 }
1790
1791
1792 besttrack->SetLabel(original->GetLabel());
1793 CookLabel(besttrack,0.0);
de5d3595 1794
de5d3595 1795 //
1796 return besttrack;
1797}
babd135a 1798