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