Method to get the total number of ITS modules (Jouri)
[u/mrichter/AliRoot.git] / ITS / AliITStrackerMI.cxx
CommitLineData
e43c066c 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
bf6adc12 16/* $Id$ */
17
e43c066c 18//-------------------------------------------------------------------------
19// Implementation of the ITS tracker class
00a7cc50 20// It reads AliITSRecPoint clusters and creates AliITStrackMI tracks
e43c066c 21// and fills with them the ESD
15dd636f 22// Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
e43c066c 23// dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
f43e383b 24//
e43c066c 25//-------------------------------------------------------------------------
bf6adc12 26
27#include <TMatrixD.h>
e43c066c 28#include <TTree.h>
bf6adc12 29#include <TTreeStream.h>
30#include <TTree.h>
31
e43c066c 32#include "AliESD.h"
6c94f330 33#include "AliV0.h"
bf6adc12 34#include "AliHelix.h"
00a7cc50 35#include "AliITSRecPoint.h"
e341247d 36#include "AliITSgeomTGeo.h"
e43c066c 37#include "AliITStrackerMI.h"
df29e9a4 38#include "AliTrackPointArray.h"
39#include "AliAlignObj.h"
e43c066c 40
41ClassImp(AliITStrackerMI)
e43c066c 42
43
44
45AliITStrackerMI::AliITSlayer AliITStrackerMI::fgLayers[kMaxLayer]; // ITS layers
8221b41b 46AliITStrackerMI::AliITStrackerMI():AliTracker(),
47fI(0),
48fBestTrack(),
49fTrackToFollow(),
50fTrackHypothesys(),
51fBestHypothesys(),
52fOriginal(),
53fCurrentEsdTrack(),
54fPass(0),
55fAfterV0(kFALSE),
56fLastLayerToTrackTo(0),
57fCoeficients(0),
58fEsd(0),
59fDebugStreamer(0){
60 //Default constructor
61}
62
e43c066c 63
e341247d 64AliITStrackerMI::AliITStrackerMI(const AliITSgeomTGeo *geom) : AliTracker(),
8221b41b 65fI(kMaxLayer),
66fBestTrack(),
67fTrackToFollow(),
68fTrackHypothesys(),
69fBestHypothesys(),
70fOriginal(),
71fCurrentEsdTrack(),
72fPass(0),
73fAfterV0(kFALSE),
74fLastLayerToTrackTo(kLastLayerToTrackTo),
75fCoeficients(0),
76fEsd(0),
77fDebugStreamer(0){
e43c066c 78 //--------------------------------------------------------------------
79 //This is the AliITStrackerMI constructor
80 //--------------------------------------------------------------------
e341247d 81 if (geom) {
82 AliWarning("\"geom\" is actually a dummy argument !");
83 }
84
628e7bb0 85 fCoeficients = 0;
86 fAfterV0 = kFALSE;
e341247d 87
88 for (Int_t i=1; i<kMaxLayer+1; i++) {
89 Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
90 Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
91
92 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
93 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
e43c066c 94 Double_t poff=TMath::ATan2(y,x);
95 Double_t zoff=z;
e341247d 96 Double_t r=TMath::Sqrt(x*x + y*y);
e43c066c 97
e341247d 98 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
e43c066c 99 r += TMath::Sqrt(x*x + y*y);
e341247d 100 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
e43c066c 101 r += TMath::Sqrt(x*x + y*y);
e341247d 102 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
e43c066c 103 r += TMath::Sqrt(x*x + y*y);
104 r*=0.25;
105
106 new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
107
108 for (Int_t j=1; j<nlad+1; j++) {
109 for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
e341247d 110 TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
111 Double_t *xyz=m.GetTranslation(), x=xyz[0], y=xyz[1];
112 Double_t *rot=m.GetRotationMatrix();
e43c066c 113
114 Double_t phi=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
115 phi+=TMath::Pi()/2;
116 if (i==1) phi+=TMath::Pi();
e341247d 117
118 if (phi<0) phi+=TMath::TwoPi();
119 else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
120
e43c066c 121 Double_t cp=TMath::Cos(phi), sp=TMath::Sin(phi);
122 Double_t r=x*cp+y*sp;
123
124 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
125 new(&det) AliITSdetector(r,phi);
126 }
127 }
128
129 }
130
131 fI=kMaxLayer;
132
133 fPass=0;
134 fConstraint[0]=1; fConstraint[1]=0;
135
136 Double_t xyz[]={kXV,kYV,kZV}, ers[]={kSigmaXV,kSigmaYV,kSigmaZV};
137 SetVertex(xyz,ers);
138
139 for (Int_t i=0; i<kMaxLayer; i++) fLayersNotToSkip[i]=kLayersNotToSkip[i];
140 fLastLayerToTrackTo=kLastLayerToTrackTo;
628e7bb0 141 for (Int_t i=0;i<100000;i++){
142 fBestTrackIndex[i]=0;
143 }
81e97e0d 144 //
145 fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
146
628e7bb0 147}
e43c066c 148
8221b41b 149AliITStrackerMI::AliITStrackerMI(const AliITStrackerMI &tracker):AliTracker(tracker),
150fI(tracker.fI),
151fBestTrack(tracker.fBestTrack),
152fTrackToFollow(tracker.fTrackToFollow),
153fTrackHypothesys(tracker.fTrackHypothesys),
154fBestHypothesys(tracker.fBestHypothesys),
155fOriginal(tracker.fOriginal),
156fCurrentEsdTrack(tracker.fCurrentEsdTrack),
157fPass(tracker.fPass),
158fAfterV0(tracker.fAfterV0),
159fLastLayerToTrackTo(tracker.fLastLayerToTrackTo),
160fCoeficients(tracker.fCoeficients),
161fEsd(tracker.fEsd),
162fDebugStreamer(tracker.fDebugStreamer){
163 //Copy constructor
164}
165
166AliITStrackerMI & AliITStrackerMI::operator=(const AliITStrackerMI &tracker){
167 //Assignment operator
168 this->~AliITStrackerMI();
169 new(this) AliITStrackerMI(tracker);
170 return *this;
171}
172
173
628e7bb0 174AliITStrackerMI::~AliITStrackerMI()
175{
176 //
177 //destructor
178 //
179 if (fCoeficients) delete []fCoeficients;
81e97e0d 180 if (fDebugStreamer) {
181 //fDebugStreamer->Close();
182 delete fDebugStreamer;
183 }
e43c066c 184}
185
186void AliITStrackerMI::SetLayersNotToSkip(Int_t *l) {
187 //--------------------------------------------------------------------
188 //This function set masks of the layers which must be not skipped
189 //--------------------------------------------------------------------
190 for (Int_t i=0; i<kMaxLayer; i++) fLayersNotToSkip[i]=l[i];
191}
192
193Int_t AliITStrackerMI::LoadClusters(TTree *cTree) {
194 //--------------------------------------------------------------------
195 //This function loads ITS clusters
196 //--------------------------------------------------------------------
00a7cc50 197 TBranch *branch=cTree->GetBranch("ITSRecPoints");
e43c066c 198 if (!branch) {
199 Error("LoadClusters"," can't get the branch !\n");
200 return 1;
201 }
202
00a7cc50 203 TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy;
e43c066c 204 branch->SetAddress(&clusters);
205
206 Int_t j=0;
207 Int_t detector=0;
208 for (Int_t i=0; i<kMaxLayer; i++) {
209 Int_t ndet=fgLayers[i].GetNdetectors();
210 Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
211 for (; j<jmax; j++) {
212 if (!cTree->GetEvent(j)) continue;
213 Int_t ncl=clusters->GetEntriesFast();
214 SignDeltas(clusters,GetZ());
215 while (ncl--) {
00a7cc50 216 AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
a504d56f 217
e43c066c 218 detector = c->GetDetectorIndex();
a504d56f 219 AliITSdetector &det=fgLayers[i].GetDetector(detector);
220
e341247d 221 //Shift the cluster to the misaligned position (temporary solution)
222 Double_t x=det.GetR();
223 if (i==0) x-=0.0075;
224 else if (i==1) x+=0.0075;
225 c->SetX(x);
a504d56f 226
00a7cc50 227 fgLayers[i].InsertCluster(new AliITSRecPoint(*c));
e43c066c 228 }
229 clusters->Delete();
230 //add dead zone virtual "cluster"
231 if (i<2){
232 for (Float_t ydead = 0; ydead < 1.31 ; ydead+=(i+1.)*0.018){
233 Int_t lab[4] = {0,0,0,detector};
234 Int_t info[3] = {0,0,0};
235 Float_t hit[5]={0,0,0.004/12.,0.001/12.,0};
236 if (i==0) hit[0] =ydead-0.4;
237 if (i==1) hit[0]=ydead-3.75;
238 hit[1] =-0.04;
239 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<2.)
00a7cc50 240 fgLayers[i].InsertCluster(new AliITSRecPoint(lab, hit, info));
e43c066c 241 hit[1]=-7.05;
242 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<2.)
00a7cc50 243 fgLayers[i].InsertCluster(new AliITSRecPoint(lab, hit, info));
e43c066c 244 hit[1]=-7.15;
245 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<2.)
00a7cc50 246 fgLayers[i].InsertCluster(new AliITSRecPoint(lab, hit, info));
e43c066c 247 hit[1] =0.06;
248 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<2.)
00a7cc50 249 fgLayers[i].InsertCluster(new AliITSRecPoint(lab, hit, info));
e43c066c 250 hit[1]=7.05;
251 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<2.)
00a7cc50 252 fgLayers[i].InsertCluster(new AliITSRecPoint(lab, hit, info));
e43c066c 253 hit[1]=7.25;
254 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<2.)
00a7cc50 255 fgLayers[i].InsertCluster(new AliITSRecPoint(lab, hit, info));
e43c066c 256 }
257 }
258
259 }
260 //
261 fgLayers[i].ResetRoad(); //road defined by the cluster density
262 fgLayers[i].SortClusters();
263 }
264
265 return 0;
266}
267
268void AliITStrackerMI::UnloadClusters() {
269 //--------------------------------------------------------------------
270 //This function unloads ITS clusters
271 //--------------------------------------------------------------------
272 for (Int_t i=0; i<kMaxLayer; i++) fgLayers[i].ResetClusters();
273}
274
15dd636f 275static Int_t CorrectForDeadZoneMaterial(AliITStrackMI *t) {
e43c066c 276 //--------------------------------------------------------------------
277 // Correction for the material between the TPC and the ITS
278 // (should it belong to the TPC code ?)
279 //--------------------------------------------------------------------
280 Double_t riw=80., diw=0.0053, x0iw=30; // TPC inner wall ?
281 Double_t rcd=61., dcd=0.0053, x0cd=30; // TPC "central drum" ?
282 Double_t yr=12.8, dr=0.03; // rods ?
283 Double_t zm=0.2, dm=0.40; // membrane
284 //Double_t rr=52., dr=0.19, x0r=24., yyr=7.77; //rails
285 Double_t rs=50., ds=0.001; // something belonging to the ITS (screen ?)
286
287 if (t->GetX() > riw) {
288 if (!t->PropagateTo(riw,diw,x0iw)) return 1;
289 if (TMath::Abs(t->GetY())>yr) t->CorrectForMaterial(dr);
290 if (TMath::Abs(t->GetZ())<zm) t->CorrectForMaterial(dm);
291 if (!t->PropagateTo(rcd,dcd,x0cd)) return 1;
292 //Double_t x,y,z; t->GetGlobalXYZat(rr,x,y,z);
293 //if (TMath::Abs(y)<yyr) t->PropagateTo(rr,dr,x0r);
294 if (!t->PropagateTo(rs,ds)) return 1;
295 } else if (t->GetX() < rs) {
296 if (!t->PropagateTo(rs,-ds)) return 1;
297 //Double_t x,y,z; t->GetGlobalXYZat(rr,x,y,z);
298 //if (TMath::Abs(y)<yyr) t->PropagateTo(rr,-dr,x0r);
299 if (!t->PropagateTo(rcd,-dcd,x0cd)) return 1;
300 if (!t->PropagateTo(riw+0.001,-diw,x0iw)) return 1;
301 } else {
302 ::Error("CorrectForDeadZoneMaterial","track is already in the dead zone !");
303 return 1;
304 }
305
306 return 0;
307}
308
309Int_t AliITStrackerMI::Clusters2Tracks(AliESD *event) {
310 //--------------------------------------------------------------------
311 // This functions reconstructs ITS tracks
312 // The clusters must be already loaded !
313 //--------------------------------------------------------------------
314 TObjArray itsTracks(15000);
628e7bb0 315 fOriginal.Clear();
81e97e0d 316 fEsd = event; // store pointer to the esd
e43c066c 317 {/* Read ESD tracks */
318 Int_t nentr=event->GetNumberOfTracks();
319 Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
320 while (nentr--) {
321 AliESDtrack *esd=event->GetTrack(nentr);
322
323 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
324 if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
325 if (esd->GetStatus()&AliESDtrack::kITSin) continue;
51ad6848 326 if (esd->GetKinkIndex(0)>0) continue; //kink daughter
15dd636f 327 AliITStrackMI *t=0;
e43c066c 328 try {
15dd636f 329 t=new AliITStrackMI(*esd);
e43c066c 330 } catch (const Char_t *msg) {
628e7bb0 331 //Warning("Clusters2Tracks",msg);
e43c066c 332 delete t;
333 continue;
334 }
791f9a2a 335 //t->fD[0] = t->GetD(GetX(),GetY());
336 //t->fD[1] = t->GetZat(GetX())-GetZ();
b9671574 337 t->GetDZ(GetX(),GetY(),GetZ(),t->GetDP()); //I.B.
338 Double_t vdist = TMath::Sqrt(t->GetD(0)*t->GetD(0)+t->GetD(1)*t->GetD(1));
628e7bb0 339 if (t->GetMass()<0.13) t->SetMass(0.13957); // MI look to the esd - mass hypothesys !!!!!!!!!!!
e43c066c 340 // write expected q
b9671574 341 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
e43c066c 342
b9671574 343 if (esd->GetV0Index(0)>0 && t->GetD(0)<30){
81e97e0d 344 //track - can be V0 according to TPC
e43c066c 345 }
81e97e0d 346 else{
b9671574 347 if (TMath::Abs(t->GetD(0))>10) {
81e97e0d 348 delete t;
349 continue;
350 }
351
352 if (TMath::Abs(vdist)>20) {
353 delete t;
354 continue;
355 }
356 if (TMath::Abs(1/t->Get1Pt())<0.120) {
357 delete t;
358 continue;
359 }
360
361 if (CorrectForDeadZoneMaterial(t)!=0) {
362 //Warning("Clusters2Tracks",
363 // "failed to correct for the material in the dead zone !\n");
364 delete t;
365 continue;
366 }
e43c066c 367 }
b9671574 368 t->SetReconstructed(kFALSE);
e43c066c 369 itsTracks.AddLast(t);
628e7bb0 370 fOriginal.AddLast(t);
e43c066c 371 }
372 } /* End Read ESD tracks */
373
374 itsTracks.Sort();
628e7bb0 375 fOriginal.Sort();
e43c066c 376 Int_t nentr=itsTracks.GetEntriesFast();
377 fTrackHypothesys.Expand(nentr);
81e97e0d 378 fBestHypothesys.Expand(nentr);
e43c066c 379 MakeCoeficients(nentr);
380 Int_t ntrk=0;
381 for (fPass=0; fPass<2; fPass++) {
382 Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
383 for (Int_t i=0; i<nentr; i++) {
4c53440c 384// cerr<<fPass<<" "<<i<<'\r';
e43c066c 385 fCurrentEsdTrack = i;
15dd636f 386 AliITStrackMI *t=(AliITStrackMI*)itsTracks.UncheckedAt(i);
e43c066c 387 if (t==0) continue; //this track has been already tracked
b9671574 388 if (t->GetReconstructed()&&(t->GetNUsed()<1.5)) continue; //this track was already "succesfully" reconstructed
791f9a2a 389 //if ( (TMath::Abs(t->GetD(GetX(),GetY())) >3.) && fConstraint[fPass]) continue;
390 //if ( (TMath::Abs(t->GetZat(GetX())-GetZ())>3.) && fConstraint[fPass]) continue;
391 Float_t dz[2]; t->GetDZ(GetX(),GetY(),GetZ(),dz); //I.B.
392 if ( (TMath::Abs(dz[0])>3.) && fConstraint[fPass]) continue;
393 if ( (TMath::Abs(dz[1])>3.) && fConstraint[fPass]) continue;
e43c066c 394
395 Int_t tpcLabel=t->GetLabel(); //save the TPC track label
396 fI = 6;
397 ResetTrackToFollow(*t);
398 ResetBestTrack();
81e97e0d 399 FollowProlongationTree(t,i,fConstraint[fPass]);
e43c066c 400
401 SortTrackHypothesys(fCurrentEsdTrack,20,0); //MI change
402 //
15dd636f 403 AliITStrackMI * besttrack = GetBestHypothesys(fCurrentEsdTrack,t,15);
e43c066c 404 if (!besttrack) continue;
405 besttrack->SetLabel(tpcLabel);
406 // besttrack->CookdEdx();
407 CookdEdx(besttrack);
b9671574 408 besttrack->SetFakeRatio(1.);
e43c066c 409 CookLabel(besttrack,0.); //For comparison only
e43c066c 410 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
628e7bb0 411
412 /*
413 if ( besttrack->GetNumberOfClusters()<6 && fConstraint[fPass]) {
e43c066c 414 continue;
415 }
416 if (besttrack->fChi2MIP[0]+besttrack->fNUsed>3.5) continue;
417 if ( (TMath::Abs(besttrack->fD[0]*besttrack->fD[0]+besttrack->fD[1]*besttrack->fD[1])>0.1) && fConstraint[fPass]) continue;
418 //delete itsTracks.RemoveAt(i);
628e7bb0 419 */
420 if (fConstraint[fPass]&&(!besttrack->IsGoldPrimary())) continue; //to be tracked also without vertex constrain
421
422
b9671574 423 t->SetReconstructed(kTRUE);
e43c066c 424 ntrk++;
425 }
426 GetBestHypothesysMIP(itsTracks);
427 }
428
429 //GetBestHypothesysMIP(itsTracks);
81e97e0d 430 UpdateTPCV0(event);
431 FindV02(event);
628e7bb0 432 fAfterV0 = kTRUE;
433 //GetBestHypothesysMIP(itsTracks);
434 //
e43c066c 435 itsTracks.Delete();
436 //
437 Int_t entries = fTrackHypothesys.GetEntriesFast();
438 for (Int_t ientry=0;ientry<entries;ientry++){
439 TObjArray * array =(TObjArray*)fTrackHypothesys.UncheckedAt(ientry);
440 if (array) array->Delete();
441 delete fTrackHypothesys.RemoveAt(ientry);
442 }
443
444 fTrackHypothesys.Delete();
81e97e0d 445 fBestHypothesys.Delete();
628e7bb0 446 fOriginal.Clear();
447 delete []fCoeficients;
448 fCoeficients=0;
e43c066c 449 Info("Clusters2Tracks","Number of prolonged tracks: %d\n",ntrk);
628e7bb0 450
e43c066c 451 return 0;
452}
453
454
e43c066c 455Int_t AliITStrackerMI::PropagateBack(AliESD *event) {
456 //--------------------------------------------------------------------
457 // This functions propagates reconstructed ITS tracks back
458 // The clusters must be loaded !
459 //--------------------------------------------------------------------
460 Int_t nentr=event->GetNumberOfTracks();
461 Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
462
463 Int_t ntrk=0;
464 for (Int_t i=0; i<nentr; i++) {
465 AliESDtrack *esd=event->GetTrack(i);
466
467 if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
468 if (esd->GetStatus()&AliESDtrack::kITSout) continue;
469
15dd636f 470 AliITStrackMI *t=0;
e43c066c 471 try {
15dd636f 472 t=new AliITStrackMI(*esd);
e43c066c 473 } catch (const Char_t *msg) {
628e7bb0 474 //Warning("PropagateBack",msg);
e43c066c 475 delete t;
476 continue;
477 }
b9671574 478 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
e43c066c 479
480 ResetTrackToFollow(*t);
481
482 // propagete to vertex [SR, GSI 17.02.2003]
483 // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
484 if (fTrackToFollow.PropagateTo(3.,0.0028,65.19)) {
6c94f330 485 if (fTrackToFollow.PropagateToVertex(event->GetVertex())) {
e43c066c 486 fTrackToFollow.StartTimeIntegral();
487 }
488 fTrackToFollow.PropagateTo(3.,-0.0028,65.19);
489 }
490
6c94f330 491 fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters();
e43c066c 492 if (RefitAt(49.,&fTrackToFollow,t)) {
493 if (CorrectForDeadZoneMaterial(&fTrackToFollow)!=0) {
628e7bb0 494 //Warning("PropagateBack",
495 // "failed to correct for the material in the dead zone !\n");
e43c066c 496 delete t;
497 continue;
498 }
499 fTrackToFollow.SetLabel(t->GetLabel());
500 //fTrackToFollow.CookdEdx();
501 CookLabel(&fTrackToFollow,0.); //For comparison only
502 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
503 //UseClusters(&fTrackToFollow);
504 ntrk++;
505 }
506 delete t;
507 }
508
509 Info("PropagateBack","Number of back propagated ITS tracks: %d\n",ntrk);
510
511 return 0;
512}
513
514Int_t AliITStrackerMI::RefitInward(AliESD *event) {
515 //--------------------------------------------------------------------
516 // This functions refits ITS tracks using the
517 // "inward propagated" TPC tracks
518 // The clusters must be loaded !
519 //--------------------------------------------------------------------
67c1e979 520 RefitV02(event);
e43c066c 521 Int_t nentr=event->GetNumberOfTracks();
522 Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
523
524 Int_t ntrk=0;
525 for (Int_t i=0; i<nentr; i++) {
526 AliESDtrack *esd=event->GetTrack(i);
527
528 if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
529 if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
530 if (esd->GetStatus()&AliESDtrack::kTPCout)
531 if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
532
15dd636f 533 AliITStrackMI *t=0;
e43c066c 534 try {
15dd636f 535 t=new AliITStrackMI(*esd);
e43c066c 536 } catch (const Char_t *msg) {
628e7bb0 537 //Warning("RefitInward",msg);
e43c066c 538 delete t;
539 continue;
540 }
b9671574 541 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
e43c066c 542 if (CorrectForDeadZoneMaterial(t)!=0) {
628e7bb0 543 //Warning("RefitInward",
544 // "failed to correct for the material in the dead zone !\n");
e43c066c 545 delete t;
546 continue;
547 }
548
549 ResetTrackToFollow(*t);
550 fTrackToFollow.ResetClusters();
551
552 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0)
6c94f330 553 fTrackToFollow.ResetCovariance(10.);
e43c066c 554
555 //Refitting...
ef7253ac 556 if (RefitAt(3.7, &fTrackToFollow, t,kTRUE)) {
e43c066c 557 fTrackToFollow.SetLabel(t->GetLabel());
558 // fTrackToFollow.CookdEdx();
559 CookdEdx(&fTrackToFollow);
560
561 CookLabel(&fTrackToFollow,0.0); //For comparison only
562
563 if (fTrackToFollow.PropagateTo(3.,0.0028,65.19)) {//The beam pipe
49d13e89 564 AliESDtrack *esdTrack =fTrackToFollow.GetESDtrack();
565 esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit);
566 Float_t r[3]={0.,0.,0.};
567 Double_t maxD=3.;
568 esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
569 ntrk++;
e43c066c 570 }
571 }
572 delete t;
573 }
574
575 Info("RefitInward","Number of refitted tracks: %d\n",ntrk);
576
577 return 0;
578}
579
580AliCluster *AliITStrackerMI::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;
586 return fgLayers[l].GetCluster(c);
587}
588
4e65ff9e 589#include "TGeoManager.h"
590#include "TGeoMatrix.h"
591#include "TGeoPhysicalNode.h"
592
df29e9a4 593Bool_t AliITStrackerMI::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
594 //
595 // Get track space point with index i
596 //
78b03929 597 Int_t l=(index & 0xf0000000) >> 28;
598 Int_t c=(index & 0x0fffffff) >> 00;
00a7cc50 599 AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
df29e9a4 600 Int_t idet = cl->GetDetectorIndex();
4e65ff9e 601
602 const char* name = AliAlignObj::SymName((AliAlignObj::ELayerID)
603 (l+AliAlignObj::kFirstLayer), idet);
604 TGeoPNEntry *mapPN = gGeoManager->GetAlignableEntry(name);
605
606 if (!mapPN) return kFALSE;
607 TGeoPhysicalNode *node = mapPN->GetPhysicalNode();
608 if (!node) {
609 gGeoManager->MakeAlignablePN(name);
610 node = mapPN->GetPhysicalNode();
611 }
612 if (!node) return kFALSE;
613 TGeoHMatrix* matrix = node->GetMatrix();
614 if (!matrix) return kFALSE;
615
616 //
617 // Calculate the global coordinates
618 //
619 Double_t localCoord[] = {cl->GetDetLocalX(), 0, cl->GetDetLocalZ()};
620 // LG AliAlignObj makes life simple but keep in mind that alignable
621 // LG volume doesn't mean sensitive volume. There might be a shift
622 // LG between the 2, has it is here for the SPD :
623 if (l<2) localCoord[1] = 0.01;
624 // LG !!! Check for this when the new geometry comes !!!
625 Double_t globalCoord[3] = {0};
626 matrix->LocalToMaster(localCoord, globalCoord);
627 Float_t xyz[3]= {globalCoord[0], globalCoord[1], globalCoord[2]};
628
629 //
630 // Calculate the cov matrix
631 //
632 TGeoRotation rotMatrix(*matrix);
633 TGeoRotation rotMatrixTr(rotMatrix.Inverse());
634 Double_t sigmaArray[] = {cl->GetSigmaY2(),0,0, 0,0,0, 0,0,cl->GetSigmaZ2()};
635 TGeoRotation sigmMatrix;
636 sigmMatrix.SetMatrix( sigmaArray );
637 sigmMatrix.MultiplyBy(&rotMatrixTr, kFALSE);
638 sigmMatrix.MultiplyBy(&rotMatrix, kTRUE);
639 const Double_t *globalSigma = sigmMatrix.GetRotationMatrix();
640 Float_t cov[6]= { globalSigma[0], globalSigma[1], globalSigma[2],
641 globalSigma[4], globalSigma[5],
642 globalSigma[8] };
643
28469fac 644 p.SetXYZ(xyz[0],xyz[1],xyz[2],cov);
df29e9a4 645 AliAlignObj::ELayerID iLayer = AliAlignObj::kInvalidLayer;
646 switch (l) {
647 case 0:
648 iLayer = AliAlignObj::kSPD1;
649 break;
650 case 1:
651 iLayer = AliAlignObj::kSPD2;
652 break;
653 case 2:
654 iLayer = AliAlignObj::kSDD1;
655 break;
656 case 3:
657 iLayer = AliAlignObj::kSDD2;
658 break;
659 case 4:
660 iLayer = AliAlignObj::kSSD1;
661 break;
662 case 5:
663 iLayer = AliAlignObj::kSSD2;
664 break;
665 default:
666 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
667 break;
668 };
669 UShort_t volid = AliAlignObj::LayerToVolUID(iLayer,idet);
670 p.SetVolumeID((UShort_t)volid);
671 return kTRUE;
672}
e43c066c 673
81e97e0d 674void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain)
e43c066c 675{
676 //--------------------------------------------------------------------
677 // Follow prolongation tree
678 //--------------------------------------------------------------------
81e97e0d 679 //
b9671574 680 AliESDtrack * esd = otrack->GetESDtrack();
81e97e0d 681 if (esd->GetV0Index(0)>0){
682 //
683 // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
684 // mapping of esd track is different as its track in Containers
685 // Need something more stable
686 // Indexes are set back againg to the ESD track indexes in UpdateTPCV0
687 for (Int_t i=0;i<3;i++){
688 Int_t index = esd->GetV0Index(i);
689 if (index==0) break;
d6a49f20 690 AliESDv0 * vertex = fEsd->GetV0(index);
81e97e0d 691 if (vertex->GetStatus()<0) continue; // rejected V0
692 //
b75d63a7 693 if (esd->GetSign()>0) {
694 vertex->SetIndex(0,esdindex);
695 }
696 else{
697 vertex->SetIndex(1,esdindex);
698 }
81e97e0d 699 }
700 }
701 TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
702 if (!bestarray){
703 bestarray = new TObjArray(5);
704 fBestHypothesys.AddAt(bestarray,esdindex);
705 }
e43c066c 706
81e97e0d 707 //
e43c066c 708 //setup tree of the prolongations
709 //
15dd636f 710 static AliITStrackMI tracks[7][100];
711 AliITStrackMI *currenttrack;
712 static AliITStrackMI currenttrack1;
713 static AliITStrackMI currenttrack2;
714 static AliITStrackMI backuptrack;
e43c066c 715 Int_t ntracks[7];
716 Int_t nindexes[7][100];
717 Float_t normalizedchi2[100];
718 for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0;
b9671574 719 otrack->SetNSkipped(0);
15dd636f 720 new (&(tracks[6][0])) AliITStrackMI(*otrack);
e43c066c 721 ntracks[6]=1;
81e97e0d 722 for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
e43c066c 723 //
724 //
725 // follow prolongations
726 for (Int_t ilayer=5;ilayer>=0;ilayer--){
727 //
728 AliITSlayer &layer=fgLayers[ilayer];
729 Double_t r=layer.GetR();
730 ntracks[ilayer]=0;
731 //
732 //
733 Int_t nskipped=0;
734 Float_t nused =0;
735 for (Int_t itrack =0;itrack<ntracks[ilayer+1];itrack++){
736 //set current track
737 if (ntracks[ilayer]>=100) break;
b9671574 738 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0) nskipped++;
739 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2.) nused++;
e43c066c 740 if (ntracks[ilayer]>15+ilayer){
b9671574 741 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0 && nskipped>4+ilayer) continue;
742 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2. && nused>3) continue;
e43c066c 743 }
744
15dd636f 745 new(&currenttrack1) AliITStrackMI(tracks[ilayer+1][nindexes[ilayer+1][itrack]]);
e43c066c 746 if (ilayer==3 || ilayer==1) {
747 Double_t rs=0.5*(fgLayers[ilayer+1].GetR() + r);
748 Double_t d=0.0034, x0=38.6;
749 if (ilayer==1) {rs=9.; d=0.0097; x0=42;}
750 if (!currenttrack1.PropagateTo(rs,d,x0)) {
751 continue;
752 }
753 }
754 //
755 //find intersection with layer
756 Double_t x,y,z;
757 if (!currenttrack1.GetGlobalXYZat(r,x,y,z)) {
758 continue;
759 }
760 Double_t phi=TMath::ATan2(y,x);
761 Int_t idet=layer.FindDetectorIndex(phi,z);
762 if (idet<0) {
763 continue;
764 }
765 //propagate to the intersection
766 const AliITSdetector &det=layer.GetDetector(idet);
767 phi=det.GetPhi();
15dd636f 768 new(&currenttrack2) AliITStrackMI(currenttrack1);
e43c066c 769 if (!currenttrack1.Propagate(phi,det.GetR())) {
770 continue;
771 }
772 currenttrack2.Propagate(phi,det.GetR()); //
773 currenttrack1.SetDetectorIndex(idet);
774 currenttrack2.SetDetectorIndex(idet);
775
776 //
777 //
778 Double_t dz=7.5*TMath::Sqrt(currenttrack1.GetSigmaZ2() + 16.*kSigmaZ2[ilayer]);
779 Double_t dy=7.5*TMath::Sqrt(currenttrack1.GetSigmaY2() + 16.*kSigmaY2[ilayer]);
780 //
781 Bool_t isBoundary=kFALSE;
782 if (currenttrack1.GetY()-dy< det.GetYmin()+0.2) isBoundary = kTRUE;
783 if (currenttrack1.GetY()+dy> det.GetYmax()-0.2) isBoundary = kTRUE;
784 if (currenttrack1.GetZ()-dz< det.GetZmin()+0.2) isBoundary = kTRUE;
785 if (currenttrack1.GetZ()+dz> det.GetZmax()-0.2) isBoundary = kTRUE;
786
787 if (isBoundary){ // track at boundary between detectors
788 Float_t maxtgl = TMath::Abs(currenttrack1.GetTgl());
789 if (maxtgl>1) maxtgl=1;
790 dz = TMath::Sqrt(dz*dz+0.25*maxtgl*maxtgl);
791 //
792 Float_t maxsnp = TMath::Abs(currenttrack1.GetSnp());
793 if (maxsnp>0.95) continue;
794 //if (maxsnp>0.5) maxsnp=0.5;
795 dy=TMath::Sqrt(dy*dy+0.25*maxsnp*maxsnp);
796 }
797
798 Double_t zmin=currenttrack1.GetZ() - dz;
799 Double_t zmax=currenttrack1.GetZ() + dz;
800 Double_t ymin=currenttrack1.GetY() + r*phi - dy;
801 Double_t ymax=currenttrack1.GetY() + r*phi + dy;
802 layer.SelectClusters(zmin,zmax,ymin,ymax);
803 //
804 // loop over all possible prolongations
805 //
806 Double_t msz=1./((currenttrack1.GetSigmaZ2() + 16.*kSigmaZ2[ilayer]));
807 Double_t msy=1./((currenttrack1.GetSigmaY2() + 16.*kSigmaY2[ilayer]));
81e97e0d 808 if (constrain){
e43c066c 809 msy/=60; msz/=60.;
810 }
811 else{
812 msy/=50; msz/=50.;
813 }
814 //
00a7cc50 815 const AliITSRecPoint *c=0; Int_t ci=-1;
e43c066c 816 Double_t chi2=12345.;
817 Int_t deadzone=0;
818 currenttrack = &currenttrack1;
819 while ((c=layer.GetNextCluster(ci))!=0) {
820 if (ntracks[ilayer]>95) break; //space for skipped clusters
821 Bool_t change =kFALSE;
822 if (c->GetQ()==0 && (deadzone==1)) continue;
823 Int_t idet=c->GetDetectorIndex();
824 if (currenttrack->GetDetectorIndex()!=idet) {
825 const AliITSdetector &det=layer.GetDetector(idet);
826 Double_t y,z;
827 if (!currenttrack2.GetProlongationFast(det.GetPhi(),det.GetR(),y,z)) continue;
828 Float_t pz = (z - c->GetZ()) , py=(y - c->GetY());
829 if (pz*pz*msz+py*py*msy>1.) continue;
830 //
15dd636f 831 new (&backuptrack) AliITStrackMI(currenttrack2);
e43c066c 832 change = kTRUE;
833 currenttrack =&currenttrack2;
834 if (!currenttrack->Propagate(det.GetPhi(),det.GetR())) {
15dd636f 835 new (currenttrack) AliITStrackMI(backuptrack);
e43c066c 836 change = kFALSE;
837 continue;
838 }
839 currenttrack->SetDetectorIndex(idet);
840 }
841 else{
842 Float_t pz = (currenttrack->GetZ() - c->GetZ()) , py=(currenttrack->GetY() - c->GetY());
843 if (pz*pz*msz+py*py*msy>1.) continue;
844 }
845
846 chi2=GetPredictedChi2MI(currenttrack,c,ilayer);
847 if (chi2<kMaxChi2s[ilayer]){
848 if (c->GetQ()==0) deadzone=1; // take dead zone only once
849 if (ntracks[ilayer]>=100) continue;
15dd636f 850 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
b9671574 851 updatetrack->SetClIndex(ilayer,0);
e43c066c 852 if (change){
15dd636f 853 new (&currenttrack2) AliITStrackMI(backuptrack);
e43c066c 854 }
855 if (c->GetQ()!=0){
856 if (!UpdateMI(updatetrack,c,chi2,(ilayer<<28)+ci)) continue;
857 updatetrack->SetSampledEdx(c->GetQ(),updatetrack->GetNumberOfClusters()-1); //b.b.
858 }
859 else {
b9671574 860 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
861 updatetrack->SetDeadZoneProbability(GetDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
e43c066c 862 }
863 if (c->IsUsed()){
b9671574 864 updatetrack->IncrementNUsed();
e43c066c 865 }
866 Double_t x0;
867 Double_t d=layer.GetThickness(updatetrack->GetY(),updatetrack->GetZ(),x0);
868 updatetrack->CorrectForMaterial(d,x0);
81e97e0d 869 if (constrain) {
b9671574 870 updatetrack->SetConstrain(constrain);
e43c066c 871 fI = ilayer;
872 Double_t d=GetEffectiveThickness(0,0); //Think of this !!!!
873 Double_t xyz[]={GetX(),GetY(),GetZ()};
874 Double_t ptfactor = 1;
875 Double_t ers[]={GetSigmaX()*ptfactor,GetSigmaY()*ptfactor,GetSigmaZ()};
876 Bool_t isPrim = kTRUE;
877 if (ilayer<4){
791f9a2a 878 //updatetrack->fD[0] = updatetrack->GetD(GetX(),GetY());
879 //updatetrack->fD[1] = updatetrack->GetZat(GetX())-GetZ();
b9671574 880 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
881 if ( TMath::Abs(updatetrack->GetD(0)/(1.+ilayer))>0.4 || TMath::Abs(updatetrack->GetD(1)/(1.+ilayer))>0.4) isPrim=kFALSE;
e43c066c 882 }
883 if (isPrim) updatetrack->Improve(d,xyz,ers);
884 } //apply vertex constrain
885 ntracks[ilayer]++;
886 } // create new hypothesy
887 } // loop over possible cluster prolongation
81e97e0d 888 // if (constrain&&itrack<2&&currenttrack1.fNSkipped==0 && deadzone==0){
b9671574 889 if (constrain&&itrack<2&&currenttrack1.GetNSkipped()==0 && deadzone==0&&ntracks[ilayer]<100){
15dd636f 890 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
b9671574 891 vtrack->SetClIndex(ilayer,0);
e43c066c 892 fI = ilayer;
893 Double_t d=GetEffectiveThickness(0,0); //Think of this !!!!
894 Double_t xyz[]={GetX(),GetY(),GetZ()};
895 Double_t ers[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
896 vtrack->Improve(d,xyz,ers);
b9671574 897 vtrack->IncrementNSkipped();
e43c066c 898 ntracks[ilayer]++;
899 }
628e7bb0 900
6c94f330 901 if (constrain&&itrack<1&&TMath::Abs(currenttrack1.GetTgl())>1.1){ //big theta -- for low mult. runs
628e7bb0 902 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
b9671574 903 vtrack->SetClIndex(ilayer,0);
628e7bb0 904 fI = ilayer;
905 Double_t d=GetEffectiveThickness(0,0); //Think of this !!!!
906 Double_t xyz[]={GetX(),GetY(),GetZ()};
907 Double_t ers[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
908 vtrack->Improve(d,xyz,ers);
b9671574 909 vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
628e7bb0 910 ntracks[ilayer]++;
911 }
912
e43c066c 913
914 } //loop over track candidates
915 //
916 //
917 Int_t accepted=0;
918
919 Int_t golds=0;
920 for (Int_t itrack=0;itrack<ntracks[ilayer];itrack++){
921 normalizedchi2[itrack] = NormalizedChi2(&tracks[ilayer][itrack],ilayer);
922 if ( normalizedchi2[itrack]<3+0.5*ilayer) golds++;
923 if (ilayer>4) accepted++;
924 else{
81e97e0d 925 if ( constrain && normalizedchi2[itrack]<kMaxNormChi2C[ilayer]+1) accepted++;
926 if (!constrain && normalizedchi2[itrack]<kMaxNormChi2NonC[ilayer]+1) accepted++;
e43c066c 927 }
928 }
929 TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE);
930 ntracks[ilayer] = TMath::Min(accepted,7+2*ilayer);
931 if (ntracks[ilayer]<golds+2+ilayer) ntracks[ilayer]=TMath::Min(golds+2+ilayer,accepted);
932 if (ntracks[ilayer]>90) ntracks[ilayer]=90;
933 } //loop over layers
934 //printf("%d\t%d\t%d\t%d\t%d\t%d\n",ntracks[0],ntracks[1],ntracks[2],ntracks[3],ntracks[4],ntracks[5]);
81e97e0d 935 Int_t max = constrain? 20: 5;
e43c066c 936
628e7bb0 937 for (Int_t i=0;i<TMath::Min(max,ntracks[0]);i++) {
15dd636f 938 AliITStrackMI & track= tracks[0][nindexes[0][i]];
628e7bb0 939 if (track.GetNumberOfClusters()<2) continue;
b9671574 940 if (!constrain&&track.GetNormChi2(0)>7.)continue;
15dd636f 941 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
e43c066c 942 }
628e7bb0 943 for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
15dd636f 944 AliITStrackMI & track= tracks[1][nindexes[1][i]];
e43c066c 945 if (track.GetNumberOfClusters()<4) continue;
b9671574 946 if (!constrain&&track.GetNormChi2(1)>7.)continue;
947 if (constrain) track.IncrementNSkipped();
81e97e0d 948 if (!constrain) {
b9671574 949 track.SetD(0,track.GetD(GetX(),GetY()));
950 track.SetNSkipped(track.GetNSkipped()+4./(4.+8.*TMath::Abs(track.GetD(0))));
951 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
952 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
e43c066c 953 }
954 }
15dd636f 955 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
e43c066c 956 }
957 //}
958
81e97e0d 959 if (!constrain){
628e7bb0 960 for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
15dd636f 961 AliITStrackMI & track= tracks[2][nindexes[2][i]];
628e7bb0 962 if (track.GetNumberOfClusters()<3) continue;
b9671574 963 if (!constrain&&track.GetNormChi2(2)>7.)continue;
964 if (constrain) track.SetNSkipped(track.GetNSkipped()+2);
81e97e0d 965 if (!constrain){
b9671574 966 track.SetD(0,track.GetD(GetX(),GetY()));
967 track.SetNSkipped(track.GetNSkipped()+7./(7.+8.*TMath::Abs(track.GetD(0))));
968 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
969 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
e43c066c 970 }
971 }
15dd636f 972 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
e43c066c 973 }
974 }
81e97e0d 975
976 if (!constrain){
977 //
978 // register best tracks - important for V0 finder
979 //
980 for (Int_t ilayer=0;ilayer<5;ilayer++){
981 if (ntracks[ilayer]==0) continue;
982 AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
983 if (track.GetNumberOfClusters()<1) continue;
984 CookLabel(&track,0);
985 bestarray->AddAt(new AliITStrackMI(track),ilayer);
986 }
987 }
988 //
989 // update TPC V0 information
990 //
b9671574 991 if (otrack->GetESDtrack()->GetV0Index(0)>0){
81e97e0d 992 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
993 for (Int_t i=0;i<3;i++){
b9671574 994 Int_t index = otrack->GetESDtrack()->GetV0Index(i);
81e97e0d 995 if (index==0) break;
d6a49f20 996 AliV0 * vertex = (AliV0*)fEsd->GetV0(index);
81e97e0d 997 if (vertex->GetStatus()<0) continue; // rejected V0
998 //
6c94f330 999 if (otrack->GetSign()>0) {
81e97e0d 1000 vertex->SetIndex(0,esdindex);
1001 }
1002 else{
1003 vertex->SetIndex(1,esdindex);
1004 }
1005 //find nearest layer with track info
b75d63a7 1006 Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
1007 Int_t nearestold = GetNearestLayer(xrp); //I.B.
81e97e0d 1008 Int_t nearest = nearestold;
1009 for (Int_t ilayer =nearest;ilayer<8;ilayer++){
1010 if (ntracks[nearest]==0){
1011 nearest = ilayer;
1012 }
1013 }
1014 //
1015 AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
1016 if (nearestold<5&&nearest<5){
b9671574 1017 Bool_t accept = track.GetNormChi2(nearest)<10;
81e97e0d 1018 if (accept){
6c94f330 1019 if (track.GetSign()>0) {
b75d63a7 1020 vertex->SetParamP(track);
81e97e0d 1021 vertex->Update(fprimvertex);
1022 // vertex->SetIndex(0,track.fESDtrack->GetID());
1023 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1024 }else{
b75d63a7 1025 vertex->SetParamN(track);
81e97e0d 1026 vertex->Update(fprimvertex);
1027 //vertex->SetIndex(1,track.fESDtrack->GetID());
1028 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1029 }
1030 vertex->SetStatus(vertex->GetStatus()+1);
1031 }else{
1032 // vertex->SetStatus(-2); // reject V0 - not enough clusters
1033 }
1034 }
1035 // if (nearestold>3){
1036// Int_t indexlayer = (ntracks[0]>0)? 0:1;
1037// if (ntracks[indexlayer]>0){
1038// AliITStrackMI & track= tracks[indexlayer][nindexes[indexlayer][0]];
1039// if (track.GetNumberOfClusters()>4&&track.fNormChi2[indexlayer]<4){
1040// vertex->SetStatus(-1); // reject V0 - clusters before
1041// }
1042// }
1043// }
1044 }
1045 }
e43c066c 1046}
1047
1048
1049AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const
1050{
1051 //--------------------------------------------------------------------
1052 //
1053 //
1054 return fgLayers[layer];
1055}
8221b41b 1056
1057AliITStrackerMI::AliITSlayer::AliITSlayer():
1058fR(0),
1059fPhiOffset(0),
1060fNladders(0),
1061fZOffset(0),
1062fNdetectors(0),
1063fDetectors(0),
1064fN(0),
1065fDy5(0),
1066fDy10(0),
1067fDy20(0),
1068fClustersCs(0),
1069fClusterIndexCs(0),
1070fYcs(0),
1071fZcs(0),
1072fNcs(0),
1073fCurrentSlice(-1),
1074fZmax(0),
1075fYmin(0),
1076fYmax(0),
1077fI(0),
1078fImax(0),
1079fSkip(0),
1080fAccepted(0),
1081fRoad(0){
e43c066c 1082 //--------------------------------------------------------------------
1083 //default AliITSlayer constructor
1084 //--------------------------------------------------------------------
e43c066c 1085 for (Int_t i=0; i<kMaxClusterPerLayer;i++) {
1086 fClusterWeight[i]=0;
1087 fClusterTracks[0][i]=-1;
1088 fClusterTracks[1][i]=-1;
1089 fClusterTracks[2][i]=-1;
1090 fClusterTracks[3][i]=-1;
1091 }
1092}
1093
1094AliITStrackerMI::AliITSlayer::
8221b41b 1095AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
1096fR(r),
1097fPhiOffset(p),
1098fNladders(nl),
1099fZOffset(z),
1100fNdetectors(nd),
1101fDetectors(0),
1102fN(0),
1103fDy5(0),
1104fDy10(0),
1105fDy20(0),
1106fClustersCs(0),
1107fClusterIndexCs(0),
1108fYcs(0),
1109fZcs(0),
1110fNcs(0),
1111fCurrentSlice(-1),
1112fZmax(0),
1113fYmin(0),
1114fYmax(0),
1115fI(0),
1116fImax(0),
1117fSkip(0),
1118fAccepted(0),
1119fRoad(0) {
e43c066c 1120 //--------------------------------------------------------------------
1121 //main AliITSlayer constructor
1122 //--------------------------------------------------------------------
e43c066c 1123 fDetectors=new AliITSdetector[fNladders*fNdetectors];
e43c066c 1124 fRoad=2*fR*TMath::Sqrt(3.14/1.);//assuming that there's only one cluster
1125}
1126
8221b41b 1127AliITStrackerMI::AliITSlayer::AliITSlayer(const AliITSlayer& layer):
1128fR(layer.fR),
1129fPhiOffset(layer.fPhiOffset),
1130fNladders(layer.fNladders),
1131fZOffset(layer.fZOffset),
1132fNdetectors(layer.fNdetectors),
1133fDetectors(layer.fDetectors),
1134fN(layer.fN),
1135fDy5(layer.fDy5),
1136fDy10(layer.fDy10),
1137fDy20(layer.fDy20),
1138fClustersCs(layer.fClustersCs),
1139fClusterIndexCs(layer.fClusterIndexCs),
1140fYcs(layer.fYcs),
1141fZcs(layer.fZcs),
1142fNcs(layer.fNcs),
1143fCurrentSlice(layer.fCurrentSlice),
1144fZmax(layer.fZmax),
1145fYmin(layer.fYmin),
1146fYmax(layer.fYmax),
1147fI(layer.fI),
1148fImax(layer.fImax),
1149fSkip(layer.fSkip),
1150fAccepted(layer.fAccepted),
1151fRoad(layer.fRoad){
1152 //Copy constructor
1153}
1154
1155
e43c066c 1156AliITStrackerMI::AliITSlayer::~AliITSlayer() {
1157 //--------------------------------------------------------------------
1158 // AliITSlayer destructor
1159 //--------------------------------------------------------------------
1160 delete[] fDetectors;
1161 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1162 for (Int_t i=0; i<kMaxClusterPerLayer;i++) {
1163 fClusterWeight[i]=0;
1164 fClusterTracks[0][i]=-1;
1165 fClusterTracks[1][i]=-1;
1166 fClusterTracks[2][i]=-1;
1167 fClusterTracks[3][i]=-1;
1168 }
1169}
1170
1171void AliITStrackerMI::AliITSlayer::ResetClusters() {
1172 //--------------------------------------------------------------------
1173 // This function removes loaded clusters
1174 //--------------------------------------------------------------------
1175 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1176 for (Int_t i=0; i<kMaxClusterPerLayer;i++){
1177 fClusterWeight[i]=0;
1178 fClusterTracks[0][i]=-1;
1179 fClusterTracks[1][i]=-1;
1180 fClusterTracks[2][i]=-1;
1181 fClusterTracks[3][i]=-1;
1182 }
1183
1184 fN=0;
1185 fI=0;
1186}
1187
1188void AliITStrackerMI::AliITSlayer::ResetWeights() {
1189 //--------------------------------------------------------------------
1190 // This function reset weights of the clusters
1191 //--------------------------------------------------------------------
1192 for (Int_t i=0; i<kMaxClusterPerLayer;i++) {
1193 fClusterWeight[i]=0;
1194 fClusterTracks[0][i]=-1;
1195 fClusterTracks[1][i]=-1;
1196 fClusterTracks[2][i]=-1;
1197 fClusterTracks[3][i]=-1;
1198 }
1199 for (Int_t i=0; i<fN;i++) {
00a7cc50 1200 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster(i);
e43c066c 1201 if (cl&&cl->IsUsed()) cl->Use();
1202 }
1203
1204}
1205
1206void AliITStrackerMI::AliITSlayer::ResetRoad() {
1207 //--------------------------------------------------------------------
1208 // This function calculates the road defined by the cluster density
1209 //--------------------------------------------------------------------
1210 Int_t n=0;
1211 for (Int_t i=0; i<fN; i++) {
1212 if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
1213 }
1214 //if (n>1) fRoad=2*fR*TMath::Sqrt(3.14/n);
1215 if (n>1) fRoad=2*fR*TMath::Sqrt(3.14/n);
1216}
1217
628e7bb0 1218
00a7cc50 1219Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSRecPoint *c) {
e43c066c 1220 //--------------------------------------------------------------------
1221 //This function adds a cluster to this layer
1222 //--------------------------------------------------------------------
1223 if (fN==kMaxClusterPerLayer) {
1224 ::Error("InsertCluster","Too many clusters !\n");
1225 return 1;
1226 }
1227 fCurrentSlice=-1;
1d4090b7 1228 fClusters[fN]=c;
1229 fN++;
e43c066c 1230 AliITSdetector &det=GetDetector(c->GetDetectorIndex());
e43c066c 1231 if (c->GetY()<det.GetYmin()) det.SetYmin(c->GetY());
1232 if (c->GetY()>det.GetYmax()) det.SetYmax(c->GetY());
1233 if (c->GetZ()<det.GetZmin()) det.SetZmin(c->GetZ());
1234 if (c->GetZ()>det.GetZmax()) det.SetZmax(c->GetZ());
1235
1236 return 0;
1237}
1238
1239void AliITStrackerMI::AliITSlayer::SortClusters()
1240{
1241 //
1242 //sort clusters
1243 //
00a7cc50 1244 AliITSRecPoint **clusters = new AliITSRecPoint*[fN];
1d4090b7 1245 Float_t *z = new Float_t[fN];
1246 Int_t * index = new Int_t[fN];
1247 //
1248 for (Int_t i=0;i<fN;i++){
1249 z[i] = fClusters[i]->GetZ();
1250 }
1251 TMath::Sort(fN,z,index,kFALSE);
1252 for (Int_t i=0;i<fN;i++){
1253 clusters[i] = fClusters[index[i]];
1254 }
1255 //
1256 for (Int_t i=0;i<fN;i++){
1257 fClusters[i] = clusters[i];
1258 fZ[i] = fClusters[i]->GetZ();
1259 AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());
1260 Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
1261 if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
1262 fY[i] = y;
1263 }
1264 delete[] index;
1265 delete[] z;
1266 delete[] clusters;
1267 //
1268
e43c066c 1269 fYB[0]=10000000;
1270 fYB[1]=-10000000;
1271 for (Int_t i=0;i<fN;i++){
1272 if (fY[i]<fYB[0]) fYB[0]=fY[i];
1273 if (fY[i]>fYB[1]) fYB[1]=fY[i];
1274 fClusterIndex[i] = i;
1275 }
1276 //
1277 // fill slices
1278 fDy5 = (fYB[1]-fYB[0])/5.;
1279 fDy10 = (fYB[1]-fYB[0])/10.;
1280 fDy20 = (fYB[1]-fYB[0])/20.;
1281 for (Int_t i=0;i<6;i++) fN5[i] =0;
1282 for (Int_t i=0;i<11;i++) fN10[i]=0;
1283 for (Int_t i=0;i<21;i++) fN20[i]=0;
1284 //
1285 for (Int_t i=0;i<6;i++) {fBy5[i][0] = fYB[0]+(i-0.75)*fDy5; fBy5[i][1] = fYB[0]+(i+0.75)*fDy5;}
1286 for (Int_t i=0;i<11;i++) {fBy10[i][0] = fYB[0]+(i-0.75)*fDy10; fBy10[i][1] = fYB[0]+(i+0.75)*fDy10;}
1287 for (Int_t i=0;i<21;i++) {fBy20[i][0] = fYB[0]+(i-0.75)*fDy20; fBy20[i][1] = fYB[0]+(i+0.75)*fDy20;}
1288 //
1289 //
1d4090b7 1290 for (Int_t i=0;i<fN;i++)
1291 for (Int_t irot=-1;irot<=1;irot++){
1292 Float_t curY = fY[i]+irot*TMath::TwoPi()*fR;
1293 // slice 5
1294 for (Int_t slice=0; slice<6;slice++){
1295 if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<kMaxClusterPerLayer5){
1296 fClusters5[slice][fN5[slice]] = fClusters[i];
1297 fY5[slice][fN5[slice]] = curY;
1298 fZ5[slice][fN5[slice]] = fZ[i];
1299 fClusterIndex5[slice][fN5[slice]]=i;
1300 fN5[slice]++;
1301 }
e43c066c 1302 }
1d4090b7 1303 // slice 10
1304 for (Int_t slice=0; slice<11;slice++){
1305 if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<kMaxClusterPerLayer10){
1306 fClusters10[slice][fN10[slice]] = fClusters[i];
1307 fY10[slice][fN10[slice]] = curY;
1308 fZ10[slice][fN10[slice]] = fZ[i];
1309 fClusterIndex10[slice][fN10[slice]]=i;
1310 fN10[slice]++;
1311 }
e43c066c 1312 }
1d4090b7 1313 // slice 20
1314 for (Int_t slice=0; slice<21;slice++){
1315 if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<kMaxClusterPerLayer20){
1316 fClusters20[slice][fN20[slice]] = fClusters[i];
1317 fY20[slice][fN20[slice]] = curY;
1318 fZ20[slice][fN20[slice]] = fZ[i];
1319 fClusterIndex20[slice][fN20[slice]]=i;
1320 fN20[slice]++;
1321 }
1322 }
e43c066c 1323 }
1d4090b7 1324
1325 //
1326 // consistency check
1327 //
1328 for (Int_t i=0;i<fN-1;i++){
1329 if (fZ[i]>fZ[i+1]){
1330 printf("Bugg\n");
e43c066c 1331 }
1332 }
1d4090b7 1333 //
1334 for (Int_t slice=0;slice<21;slice++)
1335 for (Int_t i=0;i<fN20[slice]-1;i++){
1336 if (fZ20[slice][i]>fZ20[slice][i+1]){
1337 printf("Bugg\n");
1338 }
1339 }
1340
1341
e43c066c 1342}
1343
1344
e43c066c 1345Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
1346 //--------------------------------------------------------------------
1347 // This function returns the index of the nearest cluster
1348 //--------------------------------------------------------------------
1349 Int_t ncl=0;
1350 const Float_t *zcl;
1351 if (fCurrentSlice<0) {
1352 ncl = fN;
1353 zcl = fZ;
1354 }
1355 else{
1356 ncl = fNcs;
1357 zcl = fZcs;;
1358 }
1359
1360 if (ncl==0) return 0;
1361 Int_t b=0, e=ncl-1, m=(b+e)/2;
1362 for (; b<e; m=(b+e)/2) {
1363 // if (z > fClusters[m]->GetZ()) b=m+1;
1364 if (z > zcl[m]) b=m+1;
1365 else e=m;
1366 }
1367 return m;
1368}
e43c066c 1369
1370
1371void AliITStrackerMI::AliITSlayer::
1372SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
1373 //--------------------------------------------------------------------
1374 // This function sets the "window"
1375 //--------------------------------------------------------------------
1376
628e7bb0 1377 Double_t circle=2*TMath::Pi()*fR;
e43c066c 1378 fYmin = ymin; fYmax =ymax;
1379 Float_t ymiddle = (fYmin+fYmax)*0.5;
1380 if (ymiddle<fYB[0]) {fYmin+=circle; fYmax+=circle;ymiddle+=circle;}
1381 else{
1382 if (ymiddle>fYB[1]) {fYmin-=circle; fYmax-=circle;ymiddle-=circle;}
1383 }
1384 //
1385 fCurrentSlice =-1;
1386 // defualt take all
1387 fClustersCs = fClusters;
1388 fClusterIndexCs = fClusterIndex;
1389 fYcs = fY;
1390 fZcs = fZ;
1391 fNcs = fN;
628e7bb0 1392 //
e43c066c 1393 //is in 20 slice?
1394 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
1395 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
1396 if (slice<0) slice=0;
1397 if (slice>20) slice=20;
1398 Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
1399 if (isOK) {
1400 fCurrentSlice=slice;
1401 fClustersCs = fClusters20[fCurrentSlice];
1402 fClusterIndexCs = fClusterIndex20[fCurrentSlice];
1403 fYcs = fY20[fCurrentSlice];
1404 fZcs = fZ20[fCurrentSlice];
1405 fNcs = fN20[fCurrentSlice];
1406 }
1407 }
1408 //
1409 //is in 10 slice?
1410 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
1411 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
1412 if (slice<0) slice=0;
1413 if (slice>10) slice=10;
1414 Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
1415 if (isOK) {
1416 fCurrentSlice=slice;
1417 fClustersCs = fClusters10[fCurrentSlice];
1418 fClusterIndexCs = fClusterIndex10[fCurrentSlice];
1419 fYcs = fY10[fCurrentSlice];
1420 fZcs = fZ10[fCurrentSlice];
1421 fNcs = fN10[fCurrentSlice];
1422 }
1423 }
1424 //
1425 //is in 5 slice?
1426 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
1427 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
1428 if (slice<0) slice=0;
1429 if (slice>5) slice=5;
1430 Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
1431 if ( isOK){
1432 fCurrentSlice=slice;
1433 fClustersCs = fClusters5[fCurrentSlice];
1434 fClusterIndexCs = fClusterIndex5[fCurrentSlice];
1435 fYcs = fY5[fCurrentSlice];
1436 fZcs = fZ5[fCurrentSlice];
1437 fNcs = fN5[fCurrentSlice];
1438 }
1439 }
628e7bb0 1440 //
e43c066c 1441 fI=FindClusterIndex(zmin); fZmax=zmax;
1442 fImax = TMath::Min(FindClusterIndex(zmax)+1,fNcs);
1443 fSkip = 0;
1444 fAccepted =0;
1445}
1446
e43c066c 1447
628e7bb0 1448
1449
1450Int_t AliITStrackerMI::AliITSlayer::
1451FindDetectorIndex(Double_t phi, Double_t z) const {
1452 //--------------------------------------------------------------------
1453 //This function finds the detector crossed by the track
1454 //--------------------------------------------------------------------
1455 Double_t dphi=-(phi-fPhiOffset);
1456 if (dphi < 0) dphi += 2*TMath::Pi();
1457 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
1458 Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
1459 if (np>=fNladders) np-=fNladders;
1460 if (np<0) np+=fNladders;
1461
1462 Double_t dz=fZOffset-z;
1463 Int_t nz=Int_t(dz*(fNdetectors-1)*0.5/fZOffset+0.5);
1464 if (nz>=fNdetectors) return -1;
1465 if (nz<0) return -1;
1466
1467 return np*fNdetectors + nz;
1468}
1469
1470
00a7cc50 1471const AliITSRecPoint *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci){
e43c066c 1472 //--------------------------------------------------------------------
1473 // This function returns clusters within the "window"
1474 //--------------------------------------------------------------------
1475
1476 if (fCurrentSlice<0){
1477 Double_t rpi2 = 2.*fR*TMath::Pi();
1478 for (Int_t i=fI; i<fImax; i++) {
1479 Double_t y = fY[i];
1480 if (fYmax<y) y -= rpi2;
1d4090b7 1481 if (fYmin>y) y += rpi2;
e43c066c 1482 if (y<fYmin) continue;
1483 if (y>fYmax) continue;
1484 if (fClusters[i]->GetQ()==0&&fSkip==2) continue;
1485 ci=i;
1486 fI=i+1;
1487 return fClusters[i];
1488 }
1489 }
1490 else{
1491 for (Int_t i=fI; i<fImax; i++) {
1492 if (fYcs[i]<fYmin) continue;
1493 if (fYcs[i]>fYmax) continue;
1494 if (fClustersCs[i]->GetQ()==0&&fSkip==2) continue;
1495 ci=fClusterIndexCs[i];
1496 fI=i+1;
1497 return fClustersCs[i];
1498 }
1499 }
1500 return 0;
1501}
1502
1503
1504
e43c066c 1505Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
1506const {
1507 //--------------------------------------------------------------------
1508 //This function returns the layer thickness at this point (units X0)
1509 //--------------------------------------------------------------------
1510 Double_t d=0.0085;
1511 x0=21.82;
1512 if (43<fR&&fR<45) { //SSD2
1513 Double_t dd=0.0034;
1514 d=dd;
1515 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1516 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
1517 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
1518 for (Int_t i=0; i<12; i++) {
1519 if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
1520 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1521 d+=0.0034;
1522 break;
1523 }
1524 if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
1525 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1526 d+=0.0034;
1527 break;
1528 }
1529 if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1530 if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1531 }
1532 } else
1533 if (37<fR&&fR<41) { //SSD1
1534 Double_t dd=0.0034;
1535 d=dd;
1536 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1537 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
1538 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
1539 for (Int_t i=0; i<11; i++) {
1540 if (TMath::Abs(z-3.9*i)<0.15) {
1541 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1542 d+=dd;
1543 break;
1544 }
1545 if (TMath::Abs(z+3.9*i)<0.15) {
1546 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1547 d+=dd;
1548 break;
1549 }
1550 if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1551 if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1552 }
1553 } else
1554 if (13<fR&&fR<26) { //SDD
1555 Double_t dd=0.0033;
1556 d=dd;
1557 if (TMath::Abs(y-0.00)>3.30) d+=dd;
1558
1559 if (TMath::Abs(y-1.80)<0.55) {
1560 d+=0.016;
1561 for (Int_t j=0; j<20; j++) {
1562 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1563 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1564 }
1565 }
1566 if (TMath::Abs(y+1.80)<0.55) {
1567 d+=0.016;
1568 for (Int_t j=0; j<20; j++) {
1569 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1570 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1571 }
1572 }
1573
1574 for (Int_t i=0; i<4; i++) {
1575 if (TMath::Abs(z-7.3*i)<0.60) {
1576 d+=dd;
1577 if (TMath::Abs(y-0.00)>3.30) d+=dd;
1578 break;
1579 }
1580 if (TMath::Abs(z+7.3*i)<0.60) {
1581 d+=dd;
1582 if (TMath::Abs(y-0.00)>3.30) d+=dd;
1583 break;
1584 }
1585 }
1586 } else
1587 if (6<fR&&fR<8) { //SPD2
1588 Double_t dd=0.0063; x0=21.5;
1589 d=dd;
1590 if (TMath::Abs(y-3.08)>0.5) d+=dd;
1591 //if (TMath::Abs(y-3.08)>0.45) d+=dd;
1592 if (TMath::Abs(y-3.03)<0.10) {d+=0.014;}
1593 } else
1594 if (3<fR&&fR<5) { //SPD1
1595 Double_t dd=0.0063; x0=21.5;
1596 d=dd;
1597 if (TMath::Abs(y+0.21)>0.6) d+=dd;
1598 //if (TMath::Abs(y+0.21)>0.45) d+=dd;
1599 if (TMath::Abs(y+0.10)<0.10) {d+=0.014;}
1600 }
1601
1602 return d;
1603}
1604
1605Double_t AliITStrackerMI::GetEffectiveThickness(Double_t y,Double_t z) const
1606{
1607 //--------------------------------------------------------------------
1608 //Returns the thickness between the current layer and the vertex (units X0)
1609 //--------------------------------------------------------------------
1610 Double_t d=0.0028*3*3; //beam pipe
1611 Double_t x0=0;
1612
1613 Double_t xn=fgLayers[fI].GetR();
1614 for (Int_t i=0; i<fI; i++) {
1615 Double_t xi=fgLayers[i].GetR();
1616 d+=fgLayers[i].GetThickness(y,z,x0)*xi*xi;
1617 }
1618
1619 if (fI>1) {
1620 Double_t xi=9.;
1621 d+=0.0097*xi*xi;
1622 }
1623
1624 if (fI>3) {
1625 Double_t xi=0.5*(fgLayers[3].GetR()+fgLayers[4].GetR());
1626 d+=0.0034*xi*xi;
1627 }
1628
1629 return d/(xn*xn);
1630}
1631
1632Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
1633 //--------------------------------------------------------------------
1634 // This function returns number of clusters within the "window"
1635 //--------------------------------------------------------------------
1636 Int_t ncl=0;
1637 for (Int_t i=fI; i<fN; i++) {
00a7cc50 1638 const AliITSRecPoint *c=fClusters[i];
e43c066c 1639 if (c->GetZ() > fZmax) break;
1640 if (c->IsUsed()) continue;
1641 const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
1642 Double_t y=fR*det.GetPhi() + c->GetY();
1643
1644 if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
1645 if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
1646
1647 if (y<fYmin) continue;
1648 if (y>fYmax) continue;
1649 ncl++;
1650 }
1651 return ncl;
1652}
1653
ef7253ac 1654Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *t,
1655 const AliITStrackMI *c, Bool_t extra) {
e43c066c 1656 //--------------------------------------------------------------------
1657 // This function refits the track "t" at the position "x" using
1658 // the clusters from "c"
ef7253ac 1659 // If "extra"==kTRUE,
1660 // the clusters from overlapped modules get attached to "t"
e43c066c 1661 //--------------------------------------------------------------------
1662 Int_t index[kMaxLayer];
1663 Int_t k;
1664 for (k=0; k<kMaxLayer; k++) index[k]=-1;
1665 Int_t nc=c->GetNumberOfClusters();
1666 for (k=0; k<nc; k++) {
1667 Int_t idx=c->GetClusterIndex(k),nl=(idx&0xf0000000)>>28;
1668 index[nl]=idx;
1669 }
1670
1671 Int_t from, to, step;
1672 if (xx > t->GetX()) {
1673 from=0; to=kMaxLayer;
1674 step=+1;
1675 } else {
1676 from=kMaxLayer-1; to=-1;
1677 step=-1;
1678 }
1679
1680 for (Int_t i=from; i != to; i += step) {
1681 AliITSlayer &layer=fgLayers[i];
1682 Double_t r=layer.GetR();
1683
1684 {
1685 Double_t hI=i-0.5*step;
1686 if (TMath::Abs(hI-1.5)<0.01 || TMath::Abs(hI-3.5)<0.01) {
1687 Double_t rs=0.5*(fgLayers[i-step].GetR() + r);
1688 Double_t d=0.0034, x0=38.6;
1689 if (TMath::Abs(hI-1.5)<0.01) {rs=9.; d=0.0097; x0=42;}
1690 if (!t->PropagateTo(rs,-step*d,x0)) {
1691 return kFALSE;
1692 }
1693 }
1694 }
1695
1696 // remember old position [SR, GSI 18.02.2003]
1697 Double_t oldX=0., oldY=0., oldZ=0.;
1698 if (t->IsStartedTimeIntegral() && step==1) {
1699 t->GetGlobalXYZat(t->GetX(),oldX,oldY,oldZ);
1700 }
1701 //
1702
1703 Double_t x,y,z;
1704 if (!t->GetGlobalXYZat(r,x,y,z)) {
1705 return kFALSE;
1706 }
1707 Double_t phi=TMath::ATan2(y,x);
1708 Int_t idet=layer.FindDetectorIndex(phi,z);
1709 if (idet<0) {
1710 return kFALSE;
1711 }
1712 const AliITSdetector &det=layer.GetDetector(idet);
1713 phi=det.GetPhi();
1714 if (!t->Propagate(phi,det.GetR())) {
1715 return kFALSE;
1716 }
1717 t->SetDetectorIndex(idet);
1718
00a7cc50 1719 const AliITSRecPoint *cl=0;
e43c066c 1720 Double_t maxchi2=1000.*kMaxChi2;
1721
1722 Int_t idx=index[i];
1723 if (idx>0) {
00a7cc50 1724 const AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(idx);
67c1e979 1725 if (c){
1726 if (idet != c->GetDetectorIndex()) {
1727 idet=c->GetDetectorIndex();
1728 const AliITSdetector &det=layer.GetDetector(idet);
1729 if (!t->Propagate(det.GetPhi(),det.GetR())) {
1730 return kFALSE;
1731 }
1732 t->SetDetectorIndex(idet);
1733 }
1734 //Double_t chi2=t->GetPredictedChi2(c);
1735 Int_t layer = (idx & 0xf0000000) >> 28;;
1736 Double_t chi2=GetPredictedChi2MI(t,c,layer);
1737 if (chi2<maxchi2) {
1738 cl=c;
1739 maxchi2=chi2;
1740 } else {
1741 return kFALSE;
1742 }
1743 }
1744 }
67c1e979 1745
67c1e979 1746 if (cl) {
1747 //if (!t->Update(cl,maxchi2,idx)) {
1748 if (!UpdateMI(t,cl,maxchi2,idx)) {
1749 return kFALSE;
1750 }
1751 t->SetSampledEdx(cl->GetQ(),t->GetNumberOfClusters()-1);
1752 }
1753
1754 {
1755 Double_t x0;
1756 Double_t d=layer.GetThickness(t->GetY(),t->GetZ(),x0);
1757 t->CorrectForMaterial(-step*d,x0);
1758 }
1759
ef7253ac 1760 if (extra) { //search for extra clusters
1761 AliITStrackV2 tmp(*t);
1762 Double_t dz=4*TMath::Sqrt(tmp.GetSigmaZ2()+kSigmaZ2[i]);
1763 if (dz < 0.5*TMath::Abs(tmp.GetTgl())) dz=0.5*TMath::Abs(tmp.GetTgl());
1764 Double_t dy=4*TMath::Sqrt(t->GetSigmaY2()+kSigmaY2[i]);
1765 if (dy < 0.5*TMath::Abs(tmp.GetSnp())) dy=0.5*TMath::Abs(tmp.GetSnp());
1766 Double_t zmin=t->GetZ() - dz;
1767 Double_t zmax=t->GetZ() + dz;
1768 Double_t ymin=t->GetY() + phi*r - dy;
1769 Double_t ymax=t->GetY() + phi*r + dy;
1770 layer.SelectClusters(zmin,zmax,ymin,ymax);
1771
1772 const AliITSRecPoint *c=0; Int_t ci=-1,cci=-1;
1773 Double_t maxchi2=1000.*kMaxChi2, tolerance=0.1;
1774 while ((c=layer.GetNextCluster(ci))!=0) {
1775 if (idet == c->GetDetectorIndex()) continue;
1776
1777 const AliITSdetector &det=layer.GetDetector(c->GetDetectorIndex());
1778
1779 if (!tmp.Propagate(det.GetPhi(),det.GetR())) continue;
1780
1781 if (TMath::Abs(tmp.GetZ() - c->GetZ()) > tolerance) continue;
1782 if (TMath::Abs(tmp.GetY() - c->GetY()) > tolerance) continue;
1783
1784 Double_t chi2=tmp.GetPredictedChi2(c);
1785 if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
1786 }
1787 if (cci>=0) t->SetExtraCluster(i,(i<<28)+cci);
1788 }
1789
67c1e979 1790 // track time update [SR, GSI 17.02.2003]
1791 if (t->IsStartedTimeIntegral() && step==1) {
1792 Double_t newX, newY, newZ;
1793 t->GetGlobalXYZat(t->GetX(),newX,newY,newZ);
1794 Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) +
1795 (oldZ-newZ)*(oldZ-newZ);
1796 t->AddTimeStep(TMath::Sqrt(dL2));
1797 }
1798 //
1799
1800 }
1801
1802 if (!t->PropagateTo(xx,0.,0.)) return kFALSE;
1803 return kTRUE;
1804}
1805
67c1e979 1806Bool_t
1807AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *t,const Int_t *clindex) {
1808 //--------------------------------------------------------------------
1809 // This function refits the track "t" at the position "x" using
1810 // the clusters from array
1811 //--------------------------------------------------------------------
1812 Int_t index[kMaxLayer];
1813 Int_t k;
1814 for (k=0; k<kMaxLayer; k++) index[k]=-1;
1815 //
1816 for (k=0; k<kMaxLayer; k++) {
1817 index[k]=clindex[k];
1818 }
1819
1820 Int_t from, to, step;
1821 if (xx > t->GetX()) {
1822 from=0; to=kMaxLayer;
1823 step=+1;
1824 } else {
1825 from=kMaxLayer-1; to=-1;
1826 step=-1;
1827 }
1828
1829 for (Int_t i=from; i != to; i += step) {
1830 AliITSlayer &layer=fgLayers[i];
1831 Double_t r=layer.GetR();
1832 if (step<0 && xx>r) break; //
1833 {
1834 Double_t hI=i-0.5*step;
1835 if (TMath::Abs(hI-1.5)<0.01 || TMath::Abs(hI-3.5)<0.01) {
1836 Double_t rs=0.5*(fgLayers[i-step].GetR() + r);
1837 Double_t d=0.0034, x0=38.6;
1838 if (TMath::Abs(hI-1.5)<0.01) {rs=9.; d=0.0097; x0=42;}
1839 if (!t->PropagateTo(rs,-step*d,x0)) {
1840 return kFALSE;
e43c066c 1841 }
67c1e979 1842 }
1843 }
1844
1845 // remember old position [SR, GSI 18.02.2003]
1846 Double_t oldX=0., oldY=0., oldZ=0.;
1847 if (t->IsStartedTimeIntegral() && step==1) {
1848 t->GetGlobalXYZat(t->GetX(),oldX,oldY,oldZ);
1849 }
1850 //
1851
1852 Double_t x,y,z;
1853 if (!t->GetGlobalXYZat(r,x,y,z)) {
1854 return kFALSE;
1855 }
1856 Double_t phi=TMath::ATan2(y,x);
1857 Int_t idet=layer.FindDetectorIndex(phi,z);
1858 if (idet<0) {
1859 return kFALSE;
1860 }
1861 const AliITSdetector &det=layer.GetDetector(idet);
1862 phi=det.GetPhi();
1863 if (!t->Propagate(phi,det.GetR())) {
1864 return kFALSE;
1865 }
1866 t->SetDetectorIndex(idet);
1867
00a7cc50 1868 const AliITSRecPoint *cl=0;
67c1e979 1869 Double_t maxchi2=1000.*kMaxChi2;
1870
1871 Int_t idx=index[i];
1872 if (idx>0) {
00a7cc50 1873 const AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(idx);
67c1e979 1874 if (c){
1875 if (idet != c->GetDetectorIndex()) {
1876 idet=c->GetDetectorIndex();
1877 const AliITSdetector &det=layer.GetDetector(idet);
1878 if (!t->Propagate(det.GetPhi(),det.GetR())) {
1879 return kFALSE;
1880 }
1881 t->SetDetectorIndex(idet);
1882 }
1883 //Double_t chi2=t->GetPredictedChi2(c);
1884 Int_t layer = (idx & 0xf0000000) >> 28;;
1885 Double_t chi2=GetPredictedChi2MI(t,c,layer);
1886 if (chi2<maxchi2) {
1887 cl=c;
1888 maxchi2=chi2;
1889 } else {
1890 return kFALSE;
1891 }
e43c066c 1892 }
1893 }
1894 /*
1895 if (cl==0)
1896 if (t->GetNumberOfClusters()>2) {
1897 Double_t dz=4*TMath::Sqrt(t->GetSigmaZ2()+kSigmaZ2[i]);
1898 Double_t dy=4*TMath::Sqrt(t->GetSigmaY2()+kSigmaY2[i]);
1899 Double_t zmin=t->GetZ() - dz;
1900 Double_t zmax=t->GetZ() + dz;
1901 Double_t ymin=t->GetY() + phi*r - dy;
1902 Double_t ymax=t->GetY() + phi*r + dy;
1903 layer.SelectClusters(zmin,zmax,ymin,ymax);
1904
00a7cc50 1905 const AliITSRecPoint *c=0; Int_t ci=-1;
e43c066c 1906 while ((c=layer.GetNextCluster(ci))!=0) {
1907 if (idet != c->GetDetectorIndex()) continue;
1908 Double_t chi2=t->GetPredictedChi2(c);
1909 if (chi2<maxchi2) { cl=c; maxchi2=chi2; idx=ci; }
1910 }
1911 }
1912 */
1913 if (cl) {
1914 //if (!t->Update(cl,maxchi2,idx)) {
1915 if (!UpdateMI(t,cl,maxchi2,idx)) {
1916 return kFALSE;
1917 }
1918 t->SetSampledEdx(cl->GetQ(),t->GetNumberOfClusters()-1);
1919 }
1920
1921 {
1922 Double_t x0;
1923 Double_t d=layer.GetThickness(t->GetY(),t->GetZ(),x0);
1924 t->CorrectForMaterial(-step*d,x0);
1925 }
1926
1927 // track time update [SR, GSI 17.02.2003]
1928 if (t->IsStartedTimeIntegral() && step==1) {
1929 Double_t newX, newY, newZ;
1930 t->GetGlobalXYZat(t->GetX(),newX,newY,newZ);
1931 Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) +
1932 (oldZ-newZ)*(oldZ-newZ);
1933 t->AddTimeStep(TMath::Sqrt(dL2));
1934 }
1935 //
1936
1937 }
1938
1939 if (!t->PropagateTo(xx,0.,0.)) return kFALSE;
1940 return kTRUE;
1941}
1942
15dd636f 1943Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
e43c066c 1944{
1945 //
1946 // calculate normalized chi2
1947 // return NormalizedChi2(track,0);
1948 Float_t chi2 = 0;
1949 Float_t sum=0;
1950 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
1951 // track->fdEdxMismatch=0;
1952 Float_t dedxmismatch =0;
1953 Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
1954 if (mode<100){
1955 for (Int_t i = 0;i<6;i++){
b9671574 1956 if (track->GetClIndex(i)>0){
e43c066c 1957 Float_t cerry, cerrz;
1958 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
1959 else
b9671574 1960 { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
e43c066c 1961 cerry*=cerry;
1962 cerrz*=cerrz;
b9671574 1963 Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
e43c066c 1964 if (i>1){
b9671574 1965 Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
e43c066c 1966 if (ratio<0.5) {
1967 cchi2+=(0.5-ratio)*10.;
1968 //track->fdEdxMismatch+=(0.5-ratio)*10.;
1969 dedxmismatch+=(0.5-ratio)*10.;
1970 }
1971 }
1972 if (i<2 ||i>3){
b9671574 1973 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
e43c066c 1974 Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
1975 if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
1976 if (i<2) chi2+=2*cl->GetDeltaProbability();
1977 }
1978 chi2+=cchi2;
1979 sum++;
1980 }
1981 }
b9671574 1982 if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
1983 track->SetdEdxMismatch(dedxmismatch);
e43c066c 1984 }
1985 }
1986 else{
1987 for (Int_t i = 0;i<4;i++){
b9671574 1988 if (track->GetClIndex(i)>0){
e43c066c 1989 Float_t cerry, cerrz;
1990 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
b9671574 1991 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
e43c066c 1992 cerry*=cerry;
1993 cerrz*=cerrz;
b9671574 1994 chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
1995 chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
e43c066c 1996 sum++;
1997 }
1998 }
1999 for (Int_t i = 4;i<6;i++){
b9671574 2000 if (track->GetClIndex(i)>0){
e43c066c 2001 Float_t cerry, cerrz;
2002 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
b9671574 2003 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
e43c066c 2004 cerry*=cerry;
2005 cerrz*=cerrz;
2006 Float_t cerryb, cerrzb;
2007 if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
b9671574 2008 else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
e43c066c 2009 cerryb*=cerryb;
2010 cerrzb*=cerrzb;
b9671574 2011 chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2012 chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
e43c066c 2013 sum++;
2014 }
2015 }
2016 }
b9671574 2017 if (track->GetESDtrack()->GetTPCsignal()>85){
2018 Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
e43c066c 2019 if (ratio<0.5) {
2020 chi2+=(0.5-ratio)*5.;
2021 }
2022 if (ratio>2){
2023 chi2+=(ratio-2.0)*3;
2024 }
2025 }
2026 //
b9671574 2027 Double_t match = TMath::Sqrt(track->GetChi22());
2028 if (track->GetConstrain()) match/=track->GetNumberOfClusters();
2029 if (!track->GetConstrain()) match/=track->GetNumberOfClusters()-2.;
e43c066c 2030 if (match<0) match=0;
b9671574 2031 Float_t deadzonefactor = (track->GetNDeadZone()>0) ? 3*(1.1-track->GetDeadZoneProbability()):0.;
2032 Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2033 (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2034 1./(1.+track->GetNSkipped()));
e43c066c 2035
2036 return normchi2;
2037}
2038
2039
15dd636f 2040Double_t AliITStrackerMI::GetMatchingChi2(AliITStrackMI * track1, AliITStrackMI * track2)
e43c066c 2041{
2042 //
2043 // return matching chi2 between two tracks
15dd636f 2044 AliITStrackMI track3(*track2);
e43c066c 2045 track3.Propagate(track1->GetAlpha(),track1->GetX());
2046 TMatrixD vec(5,1);
6c94f330 2047 vec(0,0)=track1->GetY() - track3.GetY();
2048 vec(1,0)=track1->GetZ() - track3.GetZ();
2049 vec(2,0)=track1->GetSnp() - track3.GetSnp();
2050 vec(3,0)=track1->GetTgl() - track3.GetTgl();
2051 vec(4,0)=track1->Get1Pt() - track3.Get1Pt();
e43c066c 2052 //
2053 TMatrixD cov(5,5);
6c94f330 2054 cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2055 cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2056 cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2057 cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2058 cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
e43c066c 2059
6c94f330 2060 cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2061 cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2062 cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2063 cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
e43c066c 2064 //
6c94f330 2065 cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2066 cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2067 cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
e43c066c 2068 //
6c94f330 2069 cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2070 cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
e43c066c 2071 //
6c94f330 2072 cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
e43c066c 2073
2074 cov.Invert();
2075 TMatrixD vec2(cov,TMatrixD::kMult,vec);
2076 TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2077 return chi2(0,0);
2078}
2079
2080Double_t AliITStrackerMI::GetDeadZoneProbability(Double_t zpos, Double_t zerr)
2081{
2082 //
2083 // return probability that given point - characterized by z position and error is in dead zone
2084 //
2085 Double_t probability =0;
2086 Double_t absz = TMath::Abs(zpos);
2087 Double_t nearestz = (absz<2)? 0.:7.1;
2088 if (TMath::Abs(absz-nearestz)>0.25+3*zerr) return 0;
2089 Double_t zmin=0, zmax=0;
2090 if (zpos<-6.){
2091 zmin = -7.25; zmax = -6.95;
2092 }
2093 if (zpos>6){
2094 zmin = 7.0; zmax =7.3;
2095 }
2096 if (absz<2){
2097 zmin = -0.75; zmax = 1.5;
2098 }
2099 probability = (TMath::Erf((zpos-zmin)/zerr) - TMath::Erf((zpos-zmax)/zerr))*0.5;
2100 return probability;
2101}
2102
2103
15dd636f 2104Double_t AliITStrackerMI::GetTruncatedChi2(AliITStrackMI * track, Float_t fac)
e43c066c 2105{
2106 //
2107 // calculate normalized chi2
2108 Float_t chi2[6];
2109 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2110 Float_t ncl = 0;
2111 for (Int_t i = 0;i<6;i++){
b9671574 2112 if (TMath::Abs(track->GetDy(i))>0){
2113 chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2114 chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
e43c066c 2115 ncl++;
2116 }
2117 else{chi2[i]=10000;}
2118 }
2119 Int_t index[6];
2120 TMath::Sort(6,chi2,index,kFALSE);
2121 Float_t max = float(ncl)*fac-1.;
2122 Float_t sumchi=0, sumweight=0;
2123 for (Int_t i=0;i<max+1;i++){
2124 Float_t weight = (i<max)?1.:(max+1.-i);
2125 sumchi+=weight*chi2[index[i]];
2126 sumweight+=weight;
2127 }
2128 Double_t normchi2 = sumchi/sumweight;
2129 return normchi2;
2130}
2131
2132
15dd636f 2133Double_t AliITStrackerMI::GetInterpolatedChi2(AliITStrackMI * forwardtrack, AliITStrackMI * backtrack)
e43c066c 2134{
2135 //
2136 // calculate normalized chi2
2137 // if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2138 Int_t npoints = 0;
2139 Double_t res =0;
2140 for (Int_t i=0;i<6;i++){
b9671574 2141 if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2142 Double_t sy1 = forwardtrack->GetSigmaY(i);
2143 Double_t sz1 = forwardtrack->GetSigmaZ(i);
2144 Double_t sy2 = backtrack->GetSigmaY(i);
2145 Double_t sz2 = backtrack->GetSigmaZ(i);
e43c066c 2146 if (i<2){ sy2=1000.;sz2=1000;}
2147 //
b9671574 2148 Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2149 Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
e43c066c 2150 //
2151 Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2152 Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2153 //
2154 res+= nz0*nz0+ny0*ny0;
2155 npoints++;
2156 }
2157 if (npoints>1) return
2158 TMath::Max(TMath::Abs(0.3*forwardtrack->Get1Pt())-0.5,0.)+
2159 //2*forwardtrack->fNUsed+
b9671574 2160 res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2161 1./(1.+forwardtrack->GetNSkipped()));
e43c066c 2162 return 1000;
2163}
2164
2165
2166
2167
2168
2169Float_t *AliITStrackerMI::GetWeight(Int_t index) {
2170 //--------------------------------------------------------------------
2171 // Return pointer to a given cluster
2172 //--------------------------------------------------------------------
2173 Int_t l=(index & 0xf0000000) >> 28;
2174 Int_t c=(index & 0x0fffffff) >> 00;
2175 return fgLayers[l].GetWeight(c);
2176}
2177
15dd636f 2178void AliITStrackerMI::RegisterClusterTracks(AliITStrackMI* track,Int_t id)
e43c066c 2179{
2180 //---------------------------------------------
2181 // register track to the list
628e7bb0 2182 //
b9671574 2183 if (track->GetESDtrack()->GetKinkIndex(0)!=0) return; //don't register kink tracks
628e7bb0 2184 //
2185 //
e43c066c 2186 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2187 Int_t index = track->GetClusterIndex(icluster);
2188 Int_t l=(index & 0xf0000000) >> 28;
2189 Int_t c=(index & 0x0fffffff) >> 00;
b9671574 2190 if (c>fgLayers[l].GetNumberOfClusters()) continue;
e43c066c 2191 for (Int_t itrack=0;itrack<4;itrack++){
b9671574 2192 if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2193 fgLayers[l].SetClusterTracks(itrack,c,id);
e43c066c 2194 break;
2195 }
2196 }
2197 }
2198}
15dd636f 2199void AliITStrackerMI::UnRegisterClusterTracks(AliITStrackMI* track, Int_t id)
e43c066c 2200{
2201 //---------------------------------------------
2202 // unregister track from the list
2203 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2204 Int_t index = track->GetClusterIndex(icluster);
2205 Int_t l=(index & 0xf0000000) >> 28;
2206 Int_t c=(index & 0x0fffffff) >> 00;
b9671574 2207 if (c>fgLayers[l].GetNumberOfClusters()) continue;
e43c066c 2208 for (Int_t itrack=0;itrack<4;itrack++){
b9671574 2209 if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2210 fgLayers[l].SetClusterTracks(itrack,c,-1);
e43c066c 2211 }
2212 }
2213 }
2214}
00a7cc50 2215Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
e43c066c 2216{
2217 //-------------------------------------------------------------
2218 //get number of shared clusters
2219 //-------------------------------------------------------------
2220 Float_t shared=0;
2221 for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2222 // mean number of clusters
2223 Float_t *ny = GetNy(id), *nz = GetNz(id);
2224
2225
2226 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2227 Int_t index = track->GetClusterIndex(icluster);
2228 Int_t l=(index & 0xf0000000) >> 28;
2229 Int_t c=(index & 0x0fffffff) >> 00;
b9671574 2230 if (c>fgLayers[l].GetNumberOfClusters()) continue;
e43c066c 2231 if (ny[l]==0){
2232 printf("problem\n");
2233 }
00a7cc50 2234 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
e43c066c 2235 Float_t weight=1;
2236 //
2237 Float_t deltan = 0;
2238 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
b9671574 2239 if (l>2&&track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
e43c066c 2240 if (l<2 || l>3){
2241 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2242 }
2243 else{
2244 deltan = (cl->GetNz()-nz[l]);
2245 }
2246 if (deltan>2.0) continue; // extended - highly probable shared cluster
2247 weight = 2./TMath::Max(3.+deltan,2.);
2248 //
2249 for (Int_t itrack=0;itrack<4;itrack++){
b9671574 2250 if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
e43c066c 2251 list[l]=index;
00a7cc50 2252 clist[l] = (AliITSRecPoint*)GetCluster(index);
e43c066c 2253 shared+=weight;
2254 break;
2255 }
2256 }
2257 }
b9671574 2258 track->SetNUsed(shared);
e43c066c 2259 return shared;
2260}
2261
15dd636f 2262Int_t AliITStrackerMI::GetOverlapTrack(AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
e43c066c 2263{
2264 //
2265 // find first shared track
2266 //
2267 // mean number of clusters
2268 Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
2269 //
2270 for (Int_t i=0;i<6;i++) overlist[i]=-1;
2271 Int_t sharedtrack=100000;
2272 Int_t tracks[24],trackindex=0;
2273 for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
2274 //
2275 for (Int_t icluster=0;icluster<6;icluster++){
2276 if (clusterlist[icluster]<0) continue;
2277 Int_t index = clusterlist[icluster];
2278 Int_t l=(index & 0xf0000000) >> 28;
2279 Int_t c=(index & 0x0fffffff) >> 00;
2280 if (ny[l]==0){
2281 printf("problem\n");
2282 }
b9671574 2283 if (c>fgLayers[l].GetNumberOfClusters()) continue;
e43c066c 2284 //if (l>3) continue;
00a7cc50 2285 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
e43c066c 2286 //
2287 Float_t deltan = 0;
2288 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
b9671574 2289 if (l>2&&track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
e43c066c 2290 if (l<2 || l>3){
2291 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2292 }
2293 else{
2294 deltan = (cl->GetNz()-nz[l]);
2295 }
2296 if (deltan>2.0) continue; // extended - highly probable shared cluster
2297 //
2298 for (Int_t itrack=3;itrack>=0;itrack--){
b9671574 2299 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2300 if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
2301 tracks[trackindex] = fgLayers[l].GetClusterTracks(itrack,c);
e43c066c 2302 trackindex++;
2303 }
2304 }
2305 }
2306 if (trackindex==0) return -1;
2307 if (trackindex==1){
2308 sharedtrack = tracks[0];
2309 }else{
2310 if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
2311 else{
2312 //
2313 Int_t track[24], cluster[24];
2314 for (Int_t i=0;i<trackindex;i++){ track[i]=-1; cluster[i]=0;}
2315 Int_t index =0;
2316 //
2317 for (Int_t i=0;i<trackindex;i++){
2318 if (tracks[i]<0) continue;
2319 track[index] = tracks[i];
2320 cluster[index]++;
2321 for (Int_t j=i+1;j<trackindex;j++){
2322 if (tracks[j]<0) continue;
2323 if (tracks[j]==tracks[i]){
2324 cluster[index]++;
2325 tracks[j]=-1;
2326 }
2327 }
2328 index++;
2329 }
2330 Int_t max=0;
2331 for (Int_t i=0;i<index;i++){
2332 if (cluster[index]>max) {
2333 sharedtrack=track[index];
2334 max=cluster[index];
2335 }
2336 }
2337 }
2338 }
2339
2340 if (sharedtrack>=100000) return -1;
2341 //
2342 // make list of overlaps
2343 shared =0;
2344 for (Int_t icluster=0;icluster<6;icluster++){
2345 if (clusterlist[icluster]<0) continue;
2346 Int_t index = clusterlist[icluster];
2347 Int_t l=(index & 0xf0000000) >> 28;
2348 Int_t c=(index & 0x0fffffff) >> 00;
b9671574 2349 if (c>fgLayers[l].GetNumberOfClusters()) continue;
00a7cc50 2350 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
e43c066c 2351 if (l==0 || l==1){
2352 if (cl->GetNy()>2) continue;
2353 if (cl->GetNz()>2) continue;
2354 }
2355 if (l==4 || l==5){
2356 if (cl->GetNy()>3) continue;
2357 if (cl->GetNz()>3) continue;
2358 }
2359 //
2360 for (Int_t itrack=3;itrack>=0;itrack--){
b9671574 2361 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2362 if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
e43c066c 2363 overlist[l]=index;
2364 shared++;
2365 }
2366 }
2367 }
2368 return sharedtrack;
2369}
2370
2371
15dd636f 2372AliITStrackMI * AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1){
e43c066c 2373 //
2374 // try to find track hypothesys without conflicts
2375 // with minimal chi2;
2376 TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
2377 Int_t entries1 = arr1->GetEntriesFast();
2378 TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
15dd636f 2379 if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
e43c066c 2380 Int_t entries2 = arr2->GetEntriesFast();
15dd636f 2381 if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
e43c066c 2382 //
15dd636f 2383 AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
2384 AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
e43c066c 2385 if (TMath::Abs(1./track10->Get1Pt())>0.5+TMath::Abs(1/track20->Get1Pt())) return track10;
2386
2387 for (Int_t itrack=0;itrack<entries1;itrack++){
15dd636f 2388 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
e43c066c 2389 UnRegisterClusterTracks(track,trackID1);
2390 }
2391 //
2392 for (Int_t itrack=0;itrack<entries2;itrack++){
15dd636f 2393 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
e43c066c 2394 UnRegisterClusterTracks(track,trackID2);
2395 }
2396 Int_t index1=0;
2397 Int_t index2=0;
2398 Float_t maxconflicts=6;
2399 Double_t maxchi2 =1000.;
2400 //
2401 // get weight of hypothesys - using best hypothesy
2402 Double_t w1,w2;
2403
2404 Int_t list1[6],list2[6];
00a7cc50 2405 AliITSRecPoint *clist1[6], *clist2[6] ;
e43c066c 2406 RegisterClusterTracks(track10,trackID1);
2407 RegisterClusterTracks(track20,trackID2);
2408 Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
2409 Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
2410 UnRegisterClusterTracks(track10,trackID1);
2411 UnRegisterClusterTracks(track20,trackID2);
2412 //
2413 // normalized chi2
2414 Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
2415 Float_t nerry[6],nerrz[6];
2416 Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
2417 Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
2418 for (Int_t i=0;i<6;i++){
2419 if ( (erry1[i]>0) && (erry2[i]>0)) {
2420 nerry[i] = TMath::Min(erry1[i],erry2[i]);
2421 nerrz[i] = TMath::Min(errz1[i],errz2[i]);
2422 }else{
2423 nerry[i] = TMath::Max(erry1[i],erry2[i]);
2424 nerrz[i] = TMath::Max(errz1[i],errz2[i]);
2425 }
b9671574 2426 if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
2427 chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
2428 chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
e43c066c 2429 ncl1++;
2430 }
b9671574 2431 if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
2432 chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
2433 chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
e43c066c 2434 ncl2++;
2435 }
2436 }
2437 chi21/=ncl1;
2438 chi22/=ncl2;
2439 //
2440 //
b9671574 2441 Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
2442 Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
e43c066c 2443 Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
2444 Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
2445 //
2446 w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
2447 +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
2448 +1.*TMath::Abs(1./track10->Get1Pt())/(TMath::Abs(1./track10->Get1Pt())+TMath::Abs(1./track20->Get1Pt()))
2449 );
2450 w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
2451 s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
2452 +1.*TMath::Abs(1./track20->Get1Pt())/(TMath::Abs(1./track10->Get1Pt())+TMath::Abs(1./track20->Get1Pt()))
2453 );
2454
2455 Double_t sumw = w1+w2;
2456 w1/=sumw;
2457 w2/=sumw;
2458 if (w1<w2*0.5) {
2459 w1 = (d2+0.5)/(d1+d2+1.);
2460 w2 = (d1+0.5)/(d1+d2+1.);
2461 }
2462 // Float_t maxmax = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
2463 //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
2464 //
2465 // get pair of "best" hypothesys
2466 //
2467 Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
2468 Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
2469
2470 for (Int_t itrack1=0;itrack1<entries1;itrack1++){
15dd636f 2471 AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
e43c066c 2472 //if (track1->fFakeRatio>0) continue;
2473 RegisterClusterTracks(track1,trackID1);
2474 for (Int_t itrack2=0;itrack2<entries2;itrack2++){
15dd636f 2475 AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
e43c066c 2476
2477 // Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
2478 //if (track2->fFakeRatio>0) continue;
2479 Float_t nskipped=0;
2480 RegisterClusterTracks(track2,trackID2);
2481 Int_t list1[6],list2[6];
00a7cc50 2482 AliITSRecPoint *clist1[6], *clist2[6] ;
e43c066c 2483 Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
2484 Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
2485 UnRegisterClusterTracks(track2,trackID2);
2486 //
b9671574 2487 if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
2488 if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
e43c066c 2489 if (nskipped>0.5) continue;
2490 //
2491 //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
2492 if (conflict1+1<cconflict1) continue;
2493 if (conflict2+1<cconflict2) continue;
2494 Float_t conflict=0;
2495 Float_t sumchi2=0;
2496 Float_t sum=0;
2497 for (Int_t i=0;i<6;i++){
2498 //
2499 Float_t c1 =0.; // conflict coeficients
2500 Float_t c2 =0.;
2501 if (clist1[i]&&clist2[i]){
2502 Float_t deltan = 0;
2503 if (i<2 || i>3){
2504 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
2505 }
2506 else{
2507 deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
2508 }
2509 c1 = 2./TMath::Max(3.+deltan,2.);
2510 c2 = 2./TMath::Max(3.+deltan,2.);
2511 }
2512 else{
2513 if (clist1[i]){
2514 Float_t deltan = 0;
2515 if (i<2 || i>3){
2516 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
2517 }
2518 else{
2519 deltan = (clist1[i]->GetNz()-nz1[i]);
2520 }
2521 c1 = 2./TMath::Max(3.+deltan,2.);
2522 c2 = 0;
2523 }
2524
2525 if (clist2[i]){
2526 Float_t deltan = 0;
2527 if (i<2 || i>3){
2528 deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
2529 }
2530 else{
2531 deltan = (clist2[i]->GetNz()-nz2[i]);
2532 }
2533 c2 = 2./TMath::Max(3.+deltan,2.);
2534 c1 = 0;
2535 }
2536 }
2537 //
2538 Double_t chi21=0,chi22=0;
b9671574 2539 if (TMath::Abs(track1->GetDy(i))>0.) {
2540 chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
2541 (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
e43c066c 2542 //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
b9671574 2543 // (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
e43c066c 2544 }else{
b9671574 2545 if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
e43c066c 2546 }
2547 //
b9671574 2548 if (TMath::Abs(track2->GetDy(i))>0.) {
2549 chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
2550 (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
e43c066c 2551 //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
2552 // (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
2553 }
2554 else{
b9671574 2555 if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
e43c066c 2556 }
2557 sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
2558 if (chi21>0) sum+=w1;
2559 if (chi22>0) sum+=w2;
2560 conflict+=(c1+c2);
2561 }
b9671574 2562 Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
2563 if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
e43c066c 2564 Double_t normchi2 = 2*conflict+sumchi2/sum;
2565 if ( normchi2 <maxchi2 ){
2566 index1 = itrack1;
2567 index2 = itrack2;
2568 maxconflicts = conflict;
2569 maxchi2 = normchi2;
2570 }
2571 }
2572 UnRegisterClusterTracks(track1,trackID1);
2573 }
2574 //
2575 // if (maxconflicts<4 && maxchi2<th0){
2576 if (maxchi2<th0*2.){
b9671574 2577 Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
15dd636f 2578 AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
b9671574 2579 track1->SetChi2MIP(5,maxconflicts);
2580 track1->SetChi2MIP(6,maxchi2);
2581 track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
e43c066c 2582 // track1->UpdateESDtrack(AliESDtrack::kITSin);
b9671574 2583 track1->SetChi2MIP(8,index1);
e43c066c 2584 fBestTrackIndex[trackID1] =index1;
2585 UpdateESDtrack(track1, AliESDtrack::kITSin);
2586 }
b9671574 2587 else if (track10->GetChi2MIP(0)<th1){
2588 track10->SetChi2MIP(5,maxconflicts);
2589 track10->SetChi2MIP(6,maxchi2);
e43c066c 2590 // track10->UpdateESDtrack(AliESDtrack::kITSin);
2591 UpdateESDtrack(track10,AliESDtrack::kITSin);
2592 }
2593
2594 for (Int_t itrack=0;itrack<entries1;itrack++){
15dd636f 2595 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
e43c066c 2596 UnRegisterClusterTracks(track,trackID1);
2597 }
2598 //
2599 for (Int_t itrack=0;itrack<entries2;itrack++){
15dd636f 2600 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
e43c066c 2601 UnRegisterClusterTracks(track,trackID2);
2602 }
2603
b9671574 2604 if (track10->GetConstrain()&&track10->GetChi2MIP(0)<kMaxChi2PerCluster[0]&&track10->GetChi2MIP(1)<kMaxChi2PerCluster[1]
2605 &&track10->GetChi2MIP(2)<kMaxChi2PerCluster[2]&&track10->GetChi2MIP(3)<kMaxChi2PerCluster[3]){
e43c066c 2606 // if (track10->fChi2MIP[0]<kMaxChi2PerCluster[0]&&track10->fChi2MIP[1]<kMaxChi2PerCluster[1]
2607 // &&track10->fChi2MIP[2]<kMaxChi2PerCluster[2]&&track10->fChi2MIP[3]<kMaxChi2PerCluster[3]){
2608 RegisterClusterTracks(track10,trackID1);
2609 }
b9671574 2610 if (track20->GetConstrain()&&track20->GetChi2MIP(0)<kMaxChi2PerCluster[0]&&track20->GetChi2MIP(1)<kMaxChi2PerCluster[1]
2611 &&track20->GetChi2MIP(2)<kMaxChi2PerCluster[2]&&track20->GetChi2MIP(3)<kMaxChi2PerCluster[3]){
e43c066c 2612 //if (track20->fChi2MIP[0]<kMaxChi2PerCluster[0]&&track20->fChi2MIP[1]<kMaxChi2PerCluster[1]
2613 // &&track20->fChi2MIP[2]<kMaxChi2PerCluster[2]&&track20->fChi2MIP[3]<kMaxChi2PerCluster[3]){
2614 RegisterClusterTracks(track20,trackID2);
2615 }
2616 return track10;
2617
2618}
2619
2620void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
2621 //--------------------------------------------------------------------
2622 // This function marks clusters assigned to the track
2623 //--------------------------------------------------------------------
2624 AliTracker::UseClusters(t,from);
2625
00a7cc50 2626 AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
e43c066c 2627 //if (c->GetQ()>2) c->Use();
2628 if (c->GetSigmaZ2()>0.1) c->Use();
00a7cc50 2629 c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
e43c066c 2630 //if (c->GetQ()>2) c->Use();
2631 if (c->GetSigmaZ2()>0.1) c->Use();
2632
2633}
2634
2635
15dd636f 2636void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
e43c066c 2637{
2638 //------------------------------------------------------------------
2639 // add track to the list of hypothesys
2640 //------------------------------------------------------------------
2641
2642 if (esdindex>=fTrackHypothesys.GetEntriesFast()) fTrackHypothesys.Expand(esdindex*2+10);
2643 //
2644 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
2645 if (!array) {
2646 array = new TObjArray(10);
2647 fTrackHypothesys.AddAt(array,esdindex);
2648 }
2649 array->AddLast(track);
2650}
2651
2652void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
2653{
2654 //-------------------------------------------------------------------
2655 // compress array of track hypothesys
2656 // keep only maxsize best hypothesys
2657 //-------------------------------------------------------------------
2658 if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
2659 if (! (fTrackHypothesys.At(esdindex)) ) return;
2660 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
2661 Int_t entries = array->GetEntriesFast();
2662 //
2663 //- find preliminary besttrack as a reference
2664 Float_t minchi2=10000;
2665 Int_t maxn=0;
15dd636f 2666 AliITStrackMI * besttrack=0;
e43c066c 2667 for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
15dd636f 2668 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
e43c066c 2669 if (!track) continue;
2670 Float_t chi2 = NormalizedChi2(track,0);
2671 //
b9671574 2672 Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
e43c066c 2673 track->SetLabel(tpcLabel);
2674 CookdEdx(track);
b9671574 2675 track->SetFakeRatio(1.);
e43c066c 2676 CookLabel(track,0.); //For comparison only
2677 //
2678 //if (chi2<kMaxChi2PerCluster[0]&&track->fFakeRatio==0){
2679 if (chi2<kMaxChi2PerCluster[0]){
2680 if (track->GetNumberOfClusters()<maxn) continue;
2681 maxn = track->GetNumberOfClusters();
2682 if (chi2<minchi2){
2683 minchi2=chi2;
2684 besttrack=track;
2685 }
2686 }
2687 else{
b9671574 2688 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
628e7bb0 2689 delete array->RemoveAt(itrack);
2690 }
e43c066c 2691 }
2692 }
2693 if (!besttrack) return;
2694 //
2695 //
2696 //take errors of best track as a reference
2697 Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
2698 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
2699 for (Int_t i=0;i<6;i++) {
b9671574 2700 if (besttrack->GetClIndex(i)>0){
2701 erry[i] = besttrack->GetSigmaY(i); erry[i+6] = besttrack->GetSigmaY(i+6);
2702 errz[i] = besttrack->GetSigmaZ(i); errz[i+6] = besttrack->GetSigmaZ(i+6);
2703 ny[i] = besttrack->GetNy(i);
2704 nz[i] = besttrack->GetNz(i);
e43c066c 2705 }
2706 }
2707 //
2708 // calculate normalized chi2
2709 //
2710 Float_t * chi2 = new Float_t[entries];
2711 Int_t * index = new Int_t[entries];
2712 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
2713 for (Int_t itrack=0;itrack<entries;itrack++){
15dd636f 2714 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
e43c066c 2715 if (track){
b9671574 2716 track->SetChi2MIP(0,GetNormalizedChi2(track, mode));
2717 if (track->GetChi2MIP(0)<kMaxChi2PerCluster[0])
2718 chi2[itrack] = track->GetChi2MIP(0);
628e7bb0 2719 else{
b9671574 2720 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
628e7bb0 2721 delete array->RemoveAt(itrack);
2722 }
2723 }
e43c066c 2724 }
2725 }
2726 //
2727 TMath::Sort(entries,chi2,index,kFALSE);
15dd636f 2728 besttrack = (AliITStrackMI*)array->At(index[0]);
b9671574 2729 if (besttrack&&besttrack->GetChi2MIP(0)<kMaxChi2PerCluster[0]){
e43c066c 2730 for (Int_t i=0;i<6;i++){
b9671574 2731 if (besttrack->GetClIndex(i)>0){
2732 erry[i] = besttrack->GetSigmaY(i); erry[i+6] = besttrack->GetSigmaY(i+6);
2733 errz[i] = besttrack->GetSigmaZ(i); erry[i+6] = besttrack->GetSigmaY(i+6);
2734 ny[i] = besttrack->GetNy(i);
2735 nz[i] = besttrack->GetNz(i);
e43c066c 2736 }
2737 }
2738 }
2739 //
2740 // calculate one more time with updated normalized errors
2741 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
2742 for (Int_t itrack=0;itrack<entries;itrack++){
15dd636f 2743 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
e43c066c 2744 if (track){
b9671574 2745 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
2746 if (track->GetChi2MIP(0)<kMaxChi2PerCluster[0])
2747 chi2[itrack] = track->GetChi2MIP(0)-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
628e7bb0 2748 else
2749 {
b9671574 2750 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
628e7bb0 2751 delete array->RemoveAt(itrack);
2752 }
2753 }
e43c066c 2754 }
2755 }
2756 entries = array->GetEntriesFast();
2757 //
628e7bb0 2758 //
e43c066c 2759 if (entries>0){
2760 TObjArray * newarray = new TObjArray();
2761 TMath::Sort(entries,chi2,index,kFALSE);
15dd636f 2762 besttrack = (AliITStrackMI*)array->At(index[0]);
e43c066c 2763 if (besttrack){
2764 //
2765 for (Int_t i=0;i<6;i++){
b9671574 2766 if (besttrack->GetNz(i)>0&&besttrack->GetNy(i)>0){
2767 erry[i] = besttrack->GetSigmaY(i); erry[i+6] = besttrack->GetSigmaY(i+6);
2768 errz[i] = besttrack->GetSigmaZ(i); errz[i+6] = besttrack->GetSigmaZ(i+6);
2769 ny[i] = besttrack->GetNy(i);
2770 nz[i] = besttrack->GetNz(i);
e43c066c 2771 }
2772 }
b9671574 2773 besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
2774 Float_t minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(kMaxChi2PerCluster[0]));
e43c066c 2775 Float_t minn = besttrack->GetNumberOfClusters()-3;
2776 Int_t accepted=0;
2777 for (Int_t i=0;i<entries;i++){
15dd636f 2778 AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
e43c066c 2779 if (!track) continue;
2780 if (accepted>maxcut) break;
b9671574 2781 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
2782 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
2783 if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
628e7bb0 2784 delete array->RemoveAt(index[i]);
2785 continue;
2786 }
e43c066c 2787 }
b9671574 2788 Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
2789 if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
628e7bb0 2790 if (!shortbest) accepted++;
e43c066c 2791 //
2792 newarray->AddLast(array->RemoveAt(index[i]));
2793 for (Int_t i=0;i<6;i++){
2794 if (nz[i]==0){
b9671574 2795 erry[i] = track->GetSigmaY(i); erry[i+6] = track->GetSigmaY(i+6);
2796 errz[i] = track->GetSigmaZ(i); errz[i] = track->GetSigmaZ(i+6);
2797 ny[i] = track->GetNy(i);
2798 nz[i] = track->GetNz(i);
e43c066c 2799 }
2800 }
2801 }
2802 else{
2803 delete array->RemoveAt(index[i]);
2804 }
2805 }
2806 array->Delete();
2807 delete fTrackHypothesys.RemoveAt(esdindex);
2808 fTrackHypothesys.AddAt(newarray,esdindex);
2809 }
2810 else{
2811 array->Delete();
2812 delete fTrackHypothesys.RemoveAt(esdindex);
2813 }
2814 }
2815 delete [] chi2;
2816 delete [] index;
2817}
2818
2819
2820
15dd636f 2821AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
e43c066c 2822{
2823 //-------------------------------------------------------------
2824 // try to find best hypothesy
2825 // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
2826 //-------------------------------------------------------------
2827 if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
2828 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
2829 if (!array) return 0;
2830 Int_t entries = array->GetEntriesFast();
2831 if (!entries) return 0;
2832 Float_t minchi2 = 100000;
15dd636f 2833 AliITStrackMI * besttrack=0;
e43c066c 2834 //
15dd636f 2835 AliITStrackMI * backtrack = new AliITStrackMI(*original);
2836 AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
628e7bb0 2837 Double_t xyzv[]={GetX(),GetY(),GetZ()};
2838 Double_t ersv[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
e43c066c 2839 //
2840 for (Int_t i=0;i<entries;i++){
15dd636f 2841 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
e43c066c 2842 if (!track) continue;
628e7bb0 2843 Float_t sigmarfi,sigmaz;
2844 GetDCASigma(track,sigmarfi,sigmaz);
b9671574 2845 track->SetDnorm(0,sigmarfi);
2846 track->SetDnorm(1,sigmaz);
628e7bb0 2847 //
b9671574 2848 track->SetChi2MIP(1,1000000);
2849 track->SetChi2MIP(2,1000000);
2850 track->SetChi2MIP(3,1000000);
e43c066c 2851 //
2852 // backtrack
15dd636f 2853 backtrack = new(backtrack) AliITStrackMI(*track);
b9671574 2854 if (track->GetConstrain()){
628e7bb0 2855 if (!backtrack->PropagateTo(3.,0.0028,65.19)) continue;
2856 if (!backtrack->Improve(0,xyzv,ersv)) continue;
6c94f330 2857 //if (!backtrack->PropagateTo(2.,0.0028,0)) continue; // This
2858 //if (!backtrack->Improve(0,xyzv,ersv)) continue; // is
2859 //if (!backtrack->PropagateTo(1.,0.0028,0)) continue; // an over-kill
2860 //if (!backtrack->Improve(0,xyzv,ersv)) continue; // (I.B.)
2861 //if (!backtrack->PropagateToVertex()) continue; //
2862 backtrack->ResetCovariance(10.);
2863 //if (!backtrack->Improve(0,xyzv,ersv)) continue;
628e7bb0 2864 }else{
6c94f330 2865 backtrack->ResetCovariance(10.);
628e7bb0 2866 }
e43c066c 2867 backtrack->ResetClusters();
628e7bb0 2868
e43c066c 2869 Double_t x = original->GetX();
2870 if (!RefitAt(x,backtrack,track)) continue;
628e7bb0 2871 //
b9671574 2872 track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
e43c066c 2873 //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
b9671574 2874 if (track->GetChi2MIP(1)>kMaxChi2PerCluster[1]*6.) continue;
2875 track->SetChi22(GetMatchingChi2(backtrack,original));
e43c066c 2876
b9671574 2877 if ((track->GetConstrain()) && track->GetChi22()>90.) continue;
2878 if ((!track->GetConstrain()) && track->GetChi22()>30.) continue;
2879 if ( track->GetChi22()/track->GetNumberOfClusters()>11.) continue;
e43c066c 2880
2881
b9671574 2882 if (!(track->GetConstrain())&&track->GetChi2MIP(1)>kMaxChi2PerCluster[1]) continue;
e43c066c 2883 Bool_t isOK=kTRUE;
e43c066c 2884 if(!isOK) continue;
2885 //
2886 //forward track - without constraint
15dd636f 2887 forwardtrack = new(forwardtrack) AliITStrackMI(*original);
e43c066c 2888 forwardtrack->ResetClusters();
2889 x = track->GetX();
628e7bb0 2890 RefitAt(x,forwardtrack,track);
b9671574 2891 track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
2892 if (track->GetChi2MIP(2)>kMaxChi2PerCluster[2]*6.0) continue;
2893 if (!(track->GetConstrain())&&track->GetChi2MIP(2)>kMaxChi2PerCluster[2]) continue;
e43c066c 2894
791f9a2a 2895 //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
2896 //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
b9671574 2897 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP()); //I.B.
2898 forwardtrack->SetD(0,track->GetD(0));
2899 forwardtrack->SetD(1,track->GetD(1));
e43c066c 2900 {
2901 Int_t list[6];
00a7cc50 2902 AliITSRecPoint* clist[6];
b9671574 2903 track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
2904 if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
e43c066c 2905 }
2906
b9671574 2907 track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
2908 if ( (track->GetChi2MIP(3)>6.*kMaxChi2PerCluster[3])) continue;
2909 if ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*kMaxChi2PerCluster[3])) {
2910 track->SetChi2MIP(3,1000);
e43c066c 2911 continue;
2912 }
b9671574 2913 Double_t chi2 = track->GetChi2MIP(0)+track->GetNUsed();
e43c066c 2914 //
2915 for (Int_t ichi=0;ichi<5;ichi++){
b9671574 2916 forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
e43c066c 2917 }
2918 if (chi2 < minchi2){
15dd636f 2919 //besttrack = new AliITStrackMI(*forwardtrack);
99c2ee26 2920 besttrack = track;
e43c066c 2921 besttrack->SetLabel(track->GetLabel());
b9671574 2922 besttrack->SetFakeRatio(track->GetFakeRatio());
e43c066c 2923 minchi2 = chi2;
791f9a2a 2924 //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
2925 //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
b9671574 2926 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP()); //I.B.
e43c066c 2927 }
2928 }
2929 delete backtrack;
2930 delete forwardtrack;
2931 Int_t accepted=0;
2932 for (Int_t i=0;i<entries;i++){
15dd636f 2933 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
e43c066c 2934 if (!track) continue;
628e7bb0 2935
b9671574 2936 if (accepted>checkmax || track->GetChi2MIP(3)>kMaxChi2PerCluster[3]*6. ||
e43c066c 2937 (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)||
b9671574 2938 track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.){
2939 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
628e7bb0 2940 delete array->RemoveAt(i);
2941 continue;
2942 }
e43c066c 2943 }
2944 else{
2945 accepted++;
2946 }
2947 }
2948 //
2949 array->Compress();
2950 SortTrackHypothesys(esdindex,checkmax,1);
2951 array = (TObjArray*) fTrackHypothesys.At(esdindex);
afdf9b9c 2952 if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
15dd636f 2953 besttrack = (AliITStrackMI*)array->At(0);
e43c066c 2954 if (!besttrack) return 0;
b9671574 2955 besttrack->SetChi2MIP(8,0);
e43c066c 2956 fBestTrackIndex[esdindex]=0;
2957 entries = array->GetEntriesFast();
15dd636f 2958 AliITStrackMI *longtrack =0;
e43c066c 2959 minchi2 =1000;
b9671574 2960 Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
e43c066c 2961 for (Int_t itrack=entries-1;itrack>0;itrack--){
15dd636f 2962 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
b9671574 2963 if (!track->GetConstrain()) continue;
2964 if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
2965 if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
2966 if (track->GetChi2MIP(0)>4.) continue;
2967 minn = track->GetNumberOfClusters()+track->GetNDeadZone();
e43c066c 2968 longtrack =track;
2969 }
2970 //if (longtrack) besttrack=longtrack;
2971
2972 Int_t list[6];
00a7cc50 2973 AliITSRecPoint * clist[6];
e43c066c 2974 Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
b9671574 2975 if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<kMaxChi2PerCluster[0]&&besttrack->GetChi2MIP(1)<kMaxChi2PerCluster[1]
2976 &&besttrack->GetChi2MIP(2)<kMaxChi2PerCluster[2]&&besttrack->GetChi2MIP(3)<kMaxChi2PerCluster[3]){
e43c066c 2977 RegisterClusterTracks(besttrack,esdindex);
2978 }
2979 //
2980 //
2981 if (shared>0.0){
2982 Int_t nshared;
2983 Int_t overlist[6];
2984 Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
2985 if (sharedtrack>=0){
2986 //
2987 besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5);
2988 if (besttrack){
2989 shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
2990 }
2991 else return 0;
2992 }
2993 }
2994
2995 if (shared>2.5) return 0;
2996 if (shared>1.0) return besttrack;
2997 //
628e7bb0 2998 // Don't sign clusters if not gold track
2999 //
3000 if (!besttrack->IsGoldPrimary()) return besttrack;
b9671574 3001 if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack; //track belong to kink
628e7bb0 3002 //
e43c066c 3003 if (fConstraint[fPass]){
3004 //
3005 // sign clusters
3006 //
3007 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3008 for (Int_t i=0;i<6;i++){
b9671574 3009 Int_t index = besttrack->GetClIndex(i);
e43c066c 3010 if (index<=0) continue;
3011 Int_t ilayer = (index & 0xf0000000) >> 28;
b9671574 3012 if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
00a7cc50 3013 AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
e43c066c 3014 if (!c) continue;
3015 if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3016 if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3017 if ( c->GetNz()> nz[i]+0.7) continue; //shared track
b9671574 3018 if ( ilayer>2&& besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
e43c066c 3019 //if ( c->GetNy()> ny[i]+0.7) continue; //shared track
3020
3021 Bool_t cansign = kTRUE;
3022 for (Int_t itrack=0;itrack<entries; itrack++){
15dd636f 3023 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
e43c066c 3024 if (!track) continue;
b9671574 3025 if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3026 if ( (track->GetClIndex(ilayer)>0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
e43c066c 3027 cansign = kFALSE;
3028 break;
3029 }
3030 }
3031 if (cansign){
b9671574 3032 if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3033 if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
e43c066c 3034 if (!c->IsUsed()) c->Use();
3035 }
3036 }
3037 }
3038 return besttrack;
3039}
3040
3041
3042
3043void AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3044{
3045 //
3046 // get "best" hypothesys
628e7bb0 3047 //
e43c066c 3048
3049 Int_t nentries = itsTracks.GetEntriesFast();
3050 for (Int_t i=0;i<nentries;i++){
15dd636f 3051 AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
e43c066c 3052 if (!track) continue;
3053 TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3054 if (!array) continue;
3055 if (array->GetEntriesFast()<=0) continue;
3056 //
15dd636f 3057 AliITStrackMI* longtrack=0;
e43c066c 3058 Float_t minn=0;
3059 Float_t maxchi2=1000;
3060 for (Int_t j=0;j<array->GetEntriesFast();j++){
15dd636f 3061 AliITStrackMI* track = (AliITStrackMI*)array->At(j);
e43c066c 3062 if (!track) continue;
b9671574 3063 if (track->GetGoldV0()) {
628e7bb0 3064 longtrack = track; //gold V0 track taken
3065 break;
3066 }
b9671574 3067 if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3068 Float_t chi2 = track->GetChi2MIP(0);
628e7bb0 3069 if (fAfterV0){
b9671574 3070 if (!track->GetGoldV0()&&track->GetConstrain()==kFALSE) chi2+=5;
628e7bb0 3071 }
b9671574 3072 if (track->GetNumberOfClusters()+track->GetNDeadZone()>minn) maxchi2 = track->GetChi2MIP(0);
628e7bb0 3073 //
3074 if (chi2 > maxchi2) continue;
b9671574 3075 minn= track->GetNumberOfClusters()+track->GetNDeadZone();
628e7bb0 3076 maxchi2 = chi2;
e43c066c 3077 longtrack=track;
e43c066c 3078 }
628e7bb0 3079 //
3080 //
3081 //
15dd636f 3082 AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
e43c066c 3083 if (!longtrack) {longtrack = besttrack;}
3084 else besttrack= longtrack;
628e7bb0 3085 //
e43c066c 3086 if (besttrack){
3087 Int_t list[6];
00a7cc50 3088 AliITSRecPoint * clist[6];
e43c066c 3089 Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3090 //
b9671574 3091 track->SetNUsed(shared);
3092 track->SetNSkipped(besttrack->GetNSkipped());
3093 track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
e43c066c 3094 if (shared>0){
3095 Int_t nshared;
3096 Int_t overlist[6];
3097 //
3098 Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3099 //if (sharedtrack==-1) sharedtrack=0;
3100 if (sharedtrack>=0){
3101 besttrack = GetBest2Tracks(i,sharedtrack,10,5.5);
3102 }
3103 }
628e7bb0 3104 if (besttrack&&fAfterV0){
3105 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3106 }
e43c066c 3107 if (besttrack&&fConstraint[fPass])
3108 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3109 //if (besttrack&&besttrack->fConstrain)
3110 // UpdateESDtrack(besttrack,AliESDtrack::kITSin);
b9671574 3111 if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5){
3112 if ( (TMath::Abs(besttrack->GetD(0))>0.1) && fConstraint[fPass]) {
3113 track->SetReconstructed(kFALSE);
e43c066c 3114 }
b9671574 3115 if ( (TMath::Abs(besttrack->GetD(1))>0.1) && fConstraint[fPass]){
3116 track->SetReconstructed(kFALSE);
e43c066c 3117 }
3118 }
3119
3120 }
3121 }
3122}
3123
3124
15dd636f 3125void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
e43c066c 3126 //--------------------------------------------------------------------
3127 //This function "cooks" a track label. If label<0, this track is fake.
3128 //--------------------------------------------------------------------
3129 Int_t tpcLabel=-1;
3130
b9671574 3131 if ( track->GetESDtrack()) tpcLabel = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
e43c066c 3132
b9671574 3133 track->SetChi2MIP(9,0);
e43c066c 3134 Int_t nwrong=0;
3135 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3136 Int_t cindex = track->GetClusterIndex(i);
3137 Int_t l=(cindex & 0xf0000000) >> 28;
00a7cc50 3138 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
e43c066c 3139 Int_t isWrong=1;
3140 for (Int_t ind=0;ind<3;ind++){
3141 if (tpcLabel>0)
3142 if (cl->GetLabel(ind)==tpcLabel) isWrong=0;
3143 }
b9671574 3144 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
e43c066c 3145 nwrong+=isWrong;
3146 }
b9671574 3147 track->SetFakeRatio(double(nwrong)/double(track->GetNumberOfClusters()));
e43c066c 3148 if (tpcLabel>0){
b9671574 3149 if (track->GetFakeRatio()>wrong) track->SetLabel(-tpcLabel);
e43c066c 3150 else
b9671574 3151 track->SetLabel(tpcLabel);
e43c066c 3152 }
3153
3154}
3155
3156
3157
15dd636f 3158void AliITStrackerMI::CookdEdx(AliITStrackMI* track)
e43c066c 3159{
3160 //
3161 //
3162 // Int_t list[6];
00a7cc50 3163 //AliITSRecPoint * clist[6];
e43c066c 3164 // Int_t shared = GetNumberOfSharedClusters(track,index,list,clist);
3165 Float_t dedx[4];
3166 Int_t accepted=0;
b9671574 3167 track->SetChi2MIP(9,0);
e43c066c 3168 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3169 Int_t cindex = track->GetClusterIndex(i);
3170 Int_t l=(cindex & 0xf0000000) >> 28;
00a7cc50 3171 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
b9671574 3172 Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
e43c066c 3173 Int_t isWrong=1;
3174 for (Int_t ind=0;ind<3;ind++){
3175 if (cl->GetLabel(ind)==lab) isWrong=0;
3176 }
b9671574 3177 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
e43c066c 3178 if (l<2) continue;
3179 //if (l>3 && (cl->GetNy()>4) || (cl->GetNz()>4)) continue; //shared track
3180 //if (l>3&& !(cl->GetType()==1||cl->GetType()==10)) continue;
3181 //if (l<4&& !(cl->GetType()==1)) continue;
b9671574 3182 dedx[accepted]= track->GetSampledEdx(i);
e43c066c 3183 //dedx[accepted]= track->fNormQ[l];
3184 accepted++;
3185 }
3186 if (accepted<1) {
3187 track->SetdEdx(0);
3188 return;
3189 }
3190 Int_t indexes[4];
3191 TMath::Sort(accepted,dedx,indexes,kFALSE);
3192 Double_t low=0.;
3193 Double_t up=0.51;
3194 Double_t nl=low*accepted, nu =up*accepted;
3195 Float_t sumamp = 0;
3196 Float_t sumweight =0;
3197 for (Int_t i=0; i<accepted; i++) {
3198 Float_t weight =1;
3199 if (i<nl+0.1) weight = TMath::Max(1.-(nl-i),0.);
3200 if (i>nu-1) weight = TMath::Max(nu-i,0.);
3201 sumamp+= dedx[indexes[i]]*weight;
3202 sumweight+=weight;
3203 }
3204 track->SetdEdx(sumamp/sumweight);
3205}
3206
3207
3208void AliITStrackerMI::MakeCoeficients(Int_t ntracks){
3209 //
3210 //
628e7bb0 3211 if (fCoeficients) delete []fCoeficients;
e43c066c 3212 fCoeficients = new Float_t[ntracks*48];
3213 for (Int_t i=0;i<ntracks*48;i++) fCoeficients[i]=-1.;
3214}
3215
3216
00a7cc50 3217Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
e43c066c 3218{
3219 //
3220 //
3221 //
3222 Float_t erry,errz;
3223 Float_t theta = track->GetTgl();
3224 Float_t phi = track->GetSnp();
3225 phi = TMath::Sqrt(phi*phi/(1.-phi*phi));
b9671574 3226 GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz);
e43c066c 3227 Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz);
3228 Float_t ny,nz;
3229 GetNTeor(layer,cluster, theta,phi,ny,nz);
3230 Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
3231 if (delta>1){
3232 chi2+=0.5*TMath::Min(delta/2,2.);
3233 chi2+=2.*cluster->GetDeltaProbability();
3234 }
3235 //
b9671574 3236 track->SetNy(layer,ny);
3237 track->SetNz(layer,nz);
3238 track->SetSigmaY(layer,erry);
3239 track->SetSigmaZ(layer, errz);
e43c066c 3240 //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
b9671574 3241 track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/(1.- track->GetSnp()*track->GetSnp())));
e43c066c 3242 return chi2;
3243
3244}
3245
00a7cc50 3246Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
e43c066c 3247{
3248 //
3249 //
3250 //
3251 Int_t layer = (index & 0xf0000000) >> 28;
b9671574 3252 track->SetClIndex(layer, index);
3253 if ( (layer>1) &&track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
3254 chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
3255 track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
e43c066c 3256 }
a504d56f 3257
b067e048 3258 if (cl->GetQ()>0) // ingore the "virtual" clusters
3259 if (!track->PropagateTo(cl->GetX(),0.,0.)) return 0; // Alignment
a504d56f 3260
b9671574 3261 return track->UpdateMI(cl->GetY(),cl->GetZ(),track->GetSigmaY(layer),track->GetSigmaZ(layer),chi2,index);
e43c066c 3262}
3263
00a7cc50 3264void AliITStrackerMI::GetNTeor(Int_t layer, const AliITSRecPoint* /*cl*/, Float_t theta, Float_t phi, Float_t &ny, Float_t &nz)
e43c066c 3265{
3266 //
3267 //get "mean shape"
3268 //
3269 if (layer==0){
3270 ny = 1.+TMath::Abs(phi)*3.2;
3271 nz = 1.+TMath::Abs(theta)*0.34;
3272 return;
3273 }
3274 if (layer==1){
3275 ny = 1.+TMath::Abs(phi)*3.2;
3276 nz = 1.+TMath::Abs(theta)*0.28;
3277 return;
3278 }
3279
3280 if (layer>3){
3281 ny = 2.02+TMath::Abs(phi)*1.95;
3282 nz = 2.02+TMath::Abs(phi)*2.35;
3283 return;
3284 }
4c53440c 3285 ny = 6.6-2.7*TMath::Abs(phi);
e43c066c 3286 nz = 2.8-3.11*TMath::Abs(phi)+0.45*TMath::Abs(theta);
3287}
3288
3289
3290