1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //-------------------------------------------------------------------------
17 // Implementation of the ITS tracker class
18 // It reads AliITSclusterV2 clusters and creates AliITStrackV2 tracks
19 // and fills with them the ESD
20 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
21 // dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
22 //-------------------------------------------------------------------------
30 #include "AliITSgeom.h"
31 #include "AliITSRecPoint.h"
33 #include "AliITSclusterV2.h"
34 #include "AliITStrackerV2.h"
36 ClassImp(AliITStrackerV2)
38 AliITStrackerV2::AliITSlayer AliITStrackerV2::fgLayers[kMaxLayer]; // ITS layers
40 AliITStrackerV2::AliITStrackerV2(const AliITSgeom *geom) : AliTracker() {
41 //--------------------------------------------------------------------
42 //This is the AliITStrackerV2 constructor
43 //--------------------------------------------------------------------
44 AliITSgeom *g=(AliITSgeom*)geom;
48 for (i=1; i<kMaxLayer+1; i++) {
49 Int_t nlad=g->GetNladders(i);
50 Int_t ndet=g->GetNdetectors(i);
52 g->GetTrans(i,1,1,x,y,z);
53 Double_t r=TMath::Sqrt(x*x + y*y);
54 Double_t poff=TMath::ATan2(y,x);
57 g->GetTrans(i,1,2,x,y,z);
58 r += TMath::Sqrt(x*x + y*y);
59 g->GetTrans(i,2,1,x,y,z);
60 r += TMath::Sqrt(x*x + y*y);
61 g->GetTrans(i,2,2,x,y,z);
62 r += TMath::Sqrt(x*x + y*y);
65 new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
67 for (Int_t j=1; j<nlad+1; j++) {
68 for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
69 Float_t x,y,zshift; g->GetTrans(i,j,k,x,y,zshift);
70 Double_t rot[9]; g->GetRotMatrix(i,j,k,rot);
72 Double_t phi=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
74 if (i==1) phi+=TMath::Pi();
75 Double_t cp=TMath::Cos(phi), sp=TMath::Sin(phi);
78 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
79 new(&det) AliITSdetector(r,phi);
88 fConstraint[0]=1; fConstraint[1]=0;
90 Double_t xyz[]={kXV,kYV,kZV}, ers[]={kSigmaXV,kSigmaYV,kSigmaZV};
93 for (Int_t i=0; i<kMaxLayer; i++) fLayersNotToSkip[i]=kLayersNotToSkip[i];
94 fLastLayerToTrackTo=kLastLayerToTrackTo;
98 void AliITStrackerV2::SetLayersNotToSkip(Int_t *l) {
99 //--------------------------------------------------------------------
100 //This function set masks of the layers which must be not skipped
101 //--------------------------------------------------------------------
102 for (Int_t i=0; i<kMaxLayer; i++) fLayersNotToSkip[i]=l[i];
105 Int_t AliITStrackerV2::LoadClusters(TTree *cTree) {
106 //--------------------------------------------------------------------
107 //This function loads ITS clusters
108 //--------------------------------------------------------------------
109 TBranch *branch=cTree->GetBranch("Clusters");
111 Error("LoadClusters"," can't get the branch !\n");
115 TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy;
116 branch->SetAddress(&clusters);
120 for (Int_t i=0; i<kMaxLayer; i++) {
121 Int_t ndet=fgLayers[i].GetNdetectors();
122 Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
123 for (; j<jmax; j++) {
124 if (!cTree->GetEvent(j)) continue;
125 Int_t ncl=clusters->GetEntriesFast();
127 AliITSclusterV2 *c=(AliITSclusterV2*)clusters->UncheckedAt(ncl);
128 detector = c->GetDetectorIndex();
129 fgLayers[i].InsertCluster(new AliITSclusterV2(*c));
132 //add dead zone virtual "cluster"
135 for (Float_t ydead = -2.; ydead < 2. ; ydead+=0.05){
136 Int_t lab[4] = {0,0,0,detector};
137 Int_t info[3] = {0,0,0};
138 Float_t hit[5]={ydead,0,1,0.01,0};
139 fgLayers[i].InsertCluster(new AliITSclusterV2(lab, hit, info));
142 fgLayers[i].InsertCluster(new AliITSclusterV2(lab, hit, info));
148 fgLayers[i].ResetRoad(); //road defined by the cluster density
154 void AliITStrackerV2::UnloadClusters() {
155 //--------------------------------------------------------------------
156 //This function unloads ITS clusters
157 //--------------------------------------------------------------------
158 for (Int_t i=0; i<kMaxLayer; i++) fgLayers[i].ResetClusters();
161 static Int_t CorrectForDeadZoneMaterial(AliITStrackV2 *t) {
162 //--------------------------------------------------------------------
163 // Correction for the material between the TPC and the ITS
164 // (should it belong to the TPC code ?)
165 //--------------------------------------------------------------------
166 Double_t riw=80., diw=0.0053, x0iw=30; // TPC inner wall ?
167 Double_t rcd=61., dcd=0.0053, x0cd=30; // TPC "central drum" ?
168 Double_t yr=12.8, dr=0.03; // rods ?
169 Double_t zm=0.2, dm=0.40; // membrane
170 //Double_t rr=52., dr=0.19, x0r=24., yyr=7.77; //rails
171 Double_t rs=50., ds=0.001; // something belonging to the ITS (screen ?)
173 if (t->GetX() > riw) {
174 if (!t->PropagateTo(riw,diw,x0iw)) return 1;
175 if (TMath::Abs(t->GetY())>yr) t->CorrectForMaterial(dr);
176 if (TMath::Abs(t->GetZ())<zm) t->CorrectForMaterial(dm);
177 if (!t->PropagateTo(rcd,dcd,x0cd)) return 1;
178 //Double_t x,y,z; t->GetGlobalXYZat(rr,x,y,z);
179 //if (TMath::Abs(y)<yyr) t->PropagateTo(rr,dr,x0r);
180 if (!t->PropagateTo(rs,ds)) return 1;
181 } else if (t->GetX() < rs) {
182 if (!t->PropagateTo(rs,-ds)) return 1;
183 //Double_t x,y,z; t->GetGlobalXYZat(rr,x,y,z);
184 //if (TMath::Abs(y)<yyr) t->PropagateTo(rr,-dr,x0r);
185 if (!t->PropagateTo(rcd,-dcd,x0cd)) return 1;
186 if (!t->PropagateTo(riw+0.001,-diw,x0iw)) return 1;
188 ::Error("CorrectForDeadZoneMaterial","track is already in the dead zone !");
195 Int_t AliITStrackerV2::Clusters2Tracks(AliESD *event) {
196 //--------------------------------------------------------------------
197 // This functions reconstructs ITS tracks
198 // The clusters must be already loaded !
199 //--------------------------------------------------------------------
200 TObjArray itsTracks(15000);
202 {/* Read ESD tracks */
203 Int_t nentr=event->GetNumberOfTracks();
204 Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
206 AliESDtrack *esd=event->GetTrack(nentr);
208 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
209 if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
210 if (esd->GetStatus()&AliESDtrack::kITSin) continue;
214 t=new AliITStrackV2(*esd);
215 } catch (const Char_t *msg) {
216 Warning("Clusters2Tracks",msg);
220 if (TMath::Abs(t->GetD())>5) {
225 if (CorrectForDeadZoneMaterial(t)!=0) {
226 Warning("Clusters2Tracks",
227 "failed to correct for the material in the dead zone !\n");
231 t->fReconstructed = kFALSE;
232 itsTracks.AddLast(t);
234 } /* End Read ESD tracks */
237 Int_t nentr=itsTracks.GetEntriesFast();
238 fTrackHypothesys.Expand(nentr);
240 for (fPass=0; fPass<2; fPass++) {
241 Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
242 for (Int_t i=0; i<nentr; i++) {
243 fCurrentEsdTrack = i;
244 AliITStrackV2 *t=(AliITStrackV2*)itsTracks.UncheckedAt(i);
245 if (t==0) continue; //this track has been already tracked
246 if (t->fReconstructed) continue; //this track was already "succesfully" reconstructed
247 if ( (TMath::Abs(t->GetD(GetX(),GetY())) >2.) && fConstraint[fPass]) continue;
248 if ( (TMath::Abs(t->GetZat(GetX())-GetZ())>2.) && fConstraint[fPass]) continue;
250 Int_t tpcLabel=t->GetLabel(); //save the TPC track label
252 ResetTrackToFollow(*t);
255 for (FollowProlongation(); fI<kMaxLayer; fI++) {
256 while (TakeNextProlongation()) {
257 FollowProlongation();
261 SortTrackHypothesys(fCurrentEsdTrack,0.9); //MI change
262 CompressTrackHypothesys(fCurrentEsdTrack,0.90,2); //MI change
265 if (fBestTrack.GetNumberOfClusters() == 0) continue;
267 if (fConstraint[fPass]) {
268 ResetTrackToFollow(*t);
269 if (!RefitAt(3.7, &fTrackToFollow, &fBestTrack)) continue;
273 fBestTrack.SetLabel(tpcLabel);
274 fBestTrack.CookdEdx();
275 CookLabel(&fBestTrack,0.); //For comparison only
276 fBestTrack.UpdateESDtrack(AliESDtrack::kITSin);
277 //UseClusters(&fBestTrack);
278 delete itsTracks.RemoveAt(i);
281 AliITStrackV2 * besttrack = GetBestHypothesys(fCurrentEsdTrack,t,5);
282 if (!besttrack) continue;
283 besttrack->SetLabel(tpcLabel);
284 besttrack->CookdEdx();
285 besttrack->fFakeRatio=1.;
286 CookLabel(besttrack,0.); //For comparison only
287 // besttrack->UpdateESDtrack(AliESDtrack::kITSin);
289 if (besttrack->GetChi2()/besttrack->GetNumberOfClusters()>3.5){
290 if ( (TMath::Abs(besttrack->GetD(GetX(),GetY()))>0.4) && fConstraint[fPass]) {
291 CompressTrackHypothesys(fCurrentEsdTrack,0.0,0);
294 if ( (TMath::Abs(besttrack->GetZat(GetX()) -GetZ() )>0.4) && fConstraint[fPass]){
295 CompressTrackHypothesys(fCurrentEsdTrack,0.0,0);
300 //delete itsTracks.RemoveAt(i);
301 t->fReconstructed = kTRUE;
307 for (Int_t i=0; i<nentr; i++) {
308 fCurrentEsdTrack = i;
309 AliITStrackV2 *t=(AliITStrackV2*)itsTracks.UncheckedAt(i);
311 Int_t tpcLabel=t->GetLabel(); //save the TPC track label
312 AliITStrackV2 * besttrack = GetBestHypothesysMIP(fCurrentEsdTrack,t);
313 if (!besttrack) continue;
315 besttrack->SetLabel(tpcLabel);
316 besttrack->CookdEdx();
317 besttrack->fFakeRatio=1.;
318 CookLabel(besttrack,0.); //For comparison only
319 besttrack->UpdateESDtrack(AliESDtrack::kITSin);
325 Int_t entries = fTrackHypothesys.GetEntriesFast();
326 for (Int_t ientry=0;ientry<entries;ientry++){
327 TObjArray * array =(TObjArray*)fTrackHypothesys.UncheckedAt(ientry);
328 if (array) array->Delete();
329 delete fTrackHypothesys.RemoveAt(ientry);
332 fTrackHypothesys.Delete();
333 Info("Clusters2Tracks","Number of prolonged tracks: %d\n",ntrk);
339 Int_t AliITStrackerV2::PropagateBack(AliESD *event) {
340 //--------------------------------------------------------------------
341 // This functions propagates reconstructed ITS tracks back
342 // The clusters must be loaded !
343 //--------------------------------------------------------------------
344 Int_t nentr=event->GetNumberOfTracks();
345 Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
348 for (Int_t i=0; i<nentr; i++) {
349 AliESDtrack *esd=event->GetTrack(i);
351 if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
352 if (esd->GetStatus()&AliESDtrack::kITSout) continue;
356 t=new AliITStrackV2(*esd);
357 } catch (const Char_t *msg) {
358 Warning("PropagateBack",msg);
363 ResetTrackToFollow(*t);
365 // propagete to vertex [SR, GSI 17.02.2003]
366 // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
367 if (fTrackToFollow.PropagateTo(3.,0.0028,65.19)) {
368 if (fTrackToFollow.PropagateToVertex()) {
369 fTrackToFollow.StartTimeIntegral();
371 fTrackToFollow.PropagateTo(3.,-0.0028,65.19);
374 fTrackToFollow.ResetCovariance(); fTrackToFollow.ResetClusters();
375 if (RefitAt(49.,&fTrackToFollow,t)) {
376 if (CorrectForDeadZoneMaterial(&fTrackToFollow)!=0) {
377 Warning("PropagateBack",
378 "failed to correct for the material in the dead zone !\n");
382 fTrackToFollow.SetLabel(t->GetLabel());
383 //fTrackToFollow.CookdEdx();
384 CookLabel(&fTrackToFollow,0.); //For comparison only
385 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
386 //UseClusters(&fTrackToFollow);
392 Info("PropagateBack","Number of back propagated ITS tracks: %d\n",ntrk);
397 Int_t AliITStrackerV2::RefitInward(AliESD *event) {
398 //--------------------------------------------------------------------
399 // This functions refits ITS tracks using the
400 // "inward propagated" TPC tracks
401 // The clusters must be loaded !
402 //--------------------------------------------------------------------
403 Int_t nentr=event->GetNumberOfTracks();
404 Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
407 for (Int_t i=0; i<nentr; i++) {
408 AliESDtrack *esd=event->GetTrack(i);
410 if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
411 if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
412 if (esd->GetStatus()&AliESDtrack::kTPCout)
413 if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
417 t=new AliITStrackV2(*esd);
418 } catch (const Char_t *msg) {
419 Warning("RefitInward",msg);
424 if (CorrectForDeadZoneMaterial(t)!=0) {
425 Warning("RefitInward",
426 "failed to correct for the material in the dead zone !\n");
431 ResetTrackToFollow(*t);
432 fTrackToFollow.ResetClusters();
434 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0)
435 fTrackToFollow.ResetCovariance();
438 if (RefitAt(3.7, &fTrackToFollow, t)) {
439 fTrackToFollow.SetLabel(t->GetLabel());
440 fTrackToFollow.CookdEdx();
441 CookLabel(&fTrackToFollow,0.0); //For comparison only
443 if (fTrackToFollow.PropagateTo(3.,0.0028,65.19)) {//The beam pipe
444 Double_t a=fTrackToFollow.GetAlpha();
445 Double_t cs=TMath::Cos(a),sn=TMath::Sin(a);
446 Double_t xv= GetX()*cs + GetY()*sn;
447 Double_t yv=-GetX()*sn + GetY()*cs;
449 Double_t c=fTrackToFollow.GetC(), snp=fTrackToFollow.GetSnp();
450 Double_t x=fTrackToFollow.GetX(), y=fTrackToFollow.GetY();
451 Double_t tgfv=-(c*(x-xv)-snp)/(c*(y-yv) + TMath::Sqrt(1.-snp*snp));
452 Double_t fv=TMath::ATan(tgfv);
454 cs=TMath::Cos(fv); sn=TMath::Sin(fv);
456 yv=-xv*sn + yv*cs; xv=x;
458 if (fTrackToFollow.Propagate(fv+a,xv)) {
459 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
460 //UseClusters(&fTrackToFollow);
462 AliITSclusterV2 c; c.SetY(yv); c.SetZ(GetZ());
463 c.SetSigmaY2(GetSigmaY()*GetSigmaY());
464 c.SetSigmaZ2(GetSigmaZ()*GetSigmaZ());
465 Double_t chi2=fTrackToFollow.GetPredictedChi2(&c);
467 if (fTrackToFollow.Update(&c,-chi2,0))
468 fTrackToFollow.SetConstrainedESDtrack(chi2);
477 Info("RefitInward","Number of refitted tracks: %d\n",ntrk);
482 AliCluster *AliITStrackerV2::GetCluster(Int_t index) const {
483 //--------------------------------------------------------------------
484 // Return pointer to a given cluster
485 //--------------------------------------------------------------------
486 Int_t l=(index & 0xf0000000) >> 28;
487 Int_t c=(index & 0x0fffffff) >> 00;
488 return fgLayers[l].GetCluster(c);
492 void AliITStrackerV2::FollowProlongation() {
493 //--------------------------------------------------------------------
494 //This function finds a track prolongation
495 //--------------------------------------------------------------------
496 while (fI>fLastLayerToTrackTo) {
499 AliITSlayer &layer=fgLayers[i];
500 AliITStrackV2 &track=fTracks[i];
502 Double_t r=layer.GetR();
505 Double_t rs=0.5*(fgLayers[i+1].GetR() + r);
506 Double_t d=0.0034, x0=38.6;
507 if (i==1) {rs=9.; d=0.0097; x0=42;}
508 if (!fTrackToFollow.PropagateTo(rs,d,x0)) {
509 //Warning("FollowProlongation","propagation failed !\n");
516 if (!fTrackToFollow.GetGlobalXYZat(r,x,y,z)) {
517 //Warning("FollowProlongation","failed to estimate track !\n");
520 Double_t phi=TMath::ATan2(y,x);
522 Int_t idet=layer.FindDetectorIndex(phi,z);
524 //Warning("FollowProlongation","failed to find a detector !\n");
528 //propagate to the intersection
529 const AliITSdetector &det=layer.GetDetector(idet);
531 if (!fTrackToFollow.Propagate(phi,det.GetR())) {
532 //Warning("FollowProlongation","propagation failed !\n");
535 fTrackToFollow.SetDetectorIndex(idet);
537 //Select possible prolongations and store the current track estimation
538 track.~AliITStrackV2(); new(&track) AliITStrackV2(fTrackToFollow);
539 Double_t dz=7*TMath::Sqrt(track.GetSigmaZ2() + kSigmaZ2[i]);
540 Double_t dy=7*TMath::Sqrt(track.GetSigmaY2() + kSigmaY2[i]);
541 Double_t road=layer.GetRoad();
542 if (dz*dy>road*road) {
543 Double_t dd=TMath::Sqrt(dz*dy), scz=dz/dd, scy=dy/dd;
544 dz=road*scz; dy=road*scy;
547 //Double_t dz=4*TMath::Sqrt(track.GetSigmaZ2() + kSigmaZ2[i]);
548 if (dz < 0.5*TMath::Abs(track.GetTgl())) dz=0.5*TMath::Abs(track.GetTgl());
550 //Warning("FollowProlongation","too broad road in Z !\n");
554 // if (TMath::Abs(fTrackToFollow.GetZ()-GetZ()) > r+dz) return;
556 //Double_t dy=4*TMath::Sqrt(track.GetSigmaY2() + kSigmaY2[i]);
557 if (TMath::Abs(track.GetSnp()>kMaxSnp)) {
561 if (dy < 0.5*TMath::Abs(track.GetSnp())) dy=0.5*TMath::Abs(track.GetSnp());
563 //Warning("FollowProlongation","too broad road in Y !\n");
567 Double_t zmin=track.GetZ() - dz;
568 Double_t zmax=track.GetZ() + dz;
569 Double_t ymin=track.GetY() + r*phi - dy;
570 Double_t ymax=track.GetY() + r*phi + dy;
571 layer.SelectClusters(zmin,zmax,ymin,ymax);
574 //take another prolongation
575 if (!TakeNextProlongation()){
577 fTrackToFollow.fNSkipped++;
578 if (fLayersNotToSkip[fI]||fTrackToFollow.fNSkipped>1) return;
580 if (fTrackToFollow.fNUsed>1) return;
581 if (fTrackToFollow.fNUsed+fTrackToFollow.fNSkipped>1) return;
582 if ( (fI<3) && ( fTrackToFollow.GetChi2()/fTrackToFollow.GetNumberOfClusters()>kChi2PerCluster)) return;
585 //deal with the best track
586 Int_t ncl=fTrackToFollow.GetNumberOfClusters();
587 Int_t nclb=fBestTrack.GetNumberOfClusters();
590 if ( (ncl<6) && (fTrackToFollow.GetChi2()/float(ncl))>3) return;
591 if (fTrackToFollow.GetChi2()/ncl>5.5) return;
592 fTrackToFollow.CookdEdx();
593 if (fTrackToFollow.fESDtrack->GetTPCsignal()>80.)
594 if ((fTrackToFollow.GetdEdx()/fTrackToFollow.fESDtrack->GetTPCsignal())<0.35){
599 fTrackToFollow.SetLabel(fTrackToFollow.fESDtrack->GetLabel());
600 CookLabel(&fTrackToFollow,0.); //
602 if ( (nclb>3) && ((fTrackToFollow.GetChi2()/ncl)<(3*fBestTrack.GetChi2()/(nclb))))
603 AddTrackHypothesys(new AliITStrackV2(fTrackToFollow), fCurrentEsdTrack);
605 if (ncl>3) AddTrackHypothesys(new AliITStrackV2(fTrackToFollow), fCurrentEsdTrack);
607 Double_t chi2=fTrackToFollow.GetChi2();
608 if (chi2/ncl < kChi2PerCluster) {
612 if ( (ncl == nclb) && chi2 < fBestTrack.GetChi2()) {
620 Int_t AliITStrackerV2::TakeNextProlongation() {
621 //--------------------------------------------------------------------
622 // This function takes another track prolongation
624 // dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
625 //--------------------------------------------------------------------
626 AliITSlayer &layer=fgLayers[fI];
627 ResetTrackToFollow(fTracks[fI]);
629 Double_t dz=7*TMath::Sqrt(fTrackToFollow.GetSigmaZ2() + kSigmaZ2[fI]);
630 Double_t dy=7*TMath::Sqrt(fTrackToFollow.GetSigmaY2() + kSigmaY2[fI]);
631 Double_t road=layer.GetRoad();
632 if (dz*dy>road*road) {
633 Double_t dd=TMath::Sqrt(dz*dy), scz=dz/dd, scy=dy/dd;
634 dz=road*scz; dy=road*scy;
637 const AliITSclusterV2 *c=0; Int_t ci=-1;
638 Double_t chi2=12345.;
639 while ((c=layer.GetNextCluster(ci))!=0) {
640 Int_t idet=c->GetDetectorIndex();
641 if (fTrackToFollow.GetDetectorIndex()!=idet) {
642 const AliITSdetector &det=layer.GetDetector(idet);
643 ResetTrackToFollow(fTracks[fI]);
644 if (!fTrackToFollow.Propagate(det.GetPhi(),det.GetR())) {
645 //Warning("TakeNextProlongation","propagation failed !\n");
648 fTrackToFollow.SetDetectorIndex(idet);
649 if (TMath::Abs(fTrackToFollow.GetZ()-GetZ())>layer.GetR()+dz) continue;
652 if (TMath::Abs(fTrackToFollow.GetZ() - c->GetZ()) > dz) continue;
653 if (TMath::Abs(fTrackToFollow.GetY() - c->GetY()) > dy) continue;
655 chi2=fTrackToFollow.GetPredictedChi2(c); if (chi2<kMaxChi2) break;
658 if (chi2>=kMaxChi2) return 0;
660 if (c->IsUsed()&&c->GetNy()<5) { //shared factor
662 chi2*=2*(1./(TMath::Max(c->GetNy(),1)));
664 if (c->GetQ()==0){ //dead zone virtual cluster
667 fTrackToFollow.fNUsed++;
670 //if ((fI<2)&&chi2>kMaxChi2In) return 0;
672 if (!fTrackToFollow.Update(c,chi2,(fI<<28)+ci)) {
673 //Warning("TakeNextProlongation","filtering failed !\n");
676 if (c->IsUsed()&&c->GetNy()<5) fTrackToFollow.fNUsed++;
677 if (fTrackToFollow.GetNumberOfClusters()>1)
678 if (TMath::Abs(fTrackToFollow.GetD())>4) return 0;
681 SetSampledEdx(c->GetQ(),fTrackToFollow.GetNumberOfClusters()-1); //b.b.
685 Double_t d=layer.GetThickness(fTrackToFollow.GetY(),fTrackToFollow.GetZ(),x0);
686 fTrackToFollow.CorrectForMaterial(d,x0);
689 if (fConstraint[fPass]) {
690 Double_t d=GetEffectiveThickness(0,0); //Think of this !!!!
691 Double_t xyz[]={GetX(),GetY(),GetZ()};
692 Double_t ers[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
693 Double_t deltad = TMath::Abs(fTrackToFollow.GetD(GetX(),GetY()));
694 Double_t deltaz = TMath::Abs(fTrackToFollow.GetZat(GetX())-GetZ());
696 if ( (fI==4) && (deltad>2.0 || deltaz>1.5)) return 0; // don't "improve" secondaries
697 if ( (fI==3) && (deltad>1.5 || deltaz>0.9)) return 0; // don't "improve" secondaries
698 if ( (fI==2) && (deltad>0.9 || deltaz>0.6)) return 1; // don't "improve" secondaries
699 if ( (fI==1) && (deltad>0.3 || deltaz>0.3)) return 1; // don't "improve" secondaries
700 if ( (fI==0) && (deltad>0.1 || deltaz>0.1)) return 1; // don't "improve" secondaries
702 fTrackToFollow.Improve(d,xyz,ers);
709 AliITStrackerV2::AliITSlayer::AliITSlayer() {
710 //--------------------------------------------------------------------
711 //default AliITSlayer constructor
712 //--------------------------------------------------------------------
717 AliITStrackerV2::AliITSlayer::
718 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd) {
719 //--------------------------------------------------------------------
720 //main AliITSlayer constructor
721 //--------------------------------------------------------------------
722 fR=r; fPhiOffset=p; fZOffset=z;
723 fNladders=nl; fNdetectors=nd;
724 fDetectors=new AliITSdetector[fNladders*fNdetectors];
729 fRoad=2*fR*TMath::Sqrt(3.14/1.);//assuming that there's only one cluster
732 AliITStrackerV2::AliITSlayer::~AliITSlayer() {
733 //--------------------------------------------------------------------
734 // AliITSlayer destructor
735 //--------------------------------------------------------------------
737 for (Int_t i=0; i<fN; i++) delete fClusters[i];
738 for (Int_t i=0; i<kMaxClusterPerLayer;i++) fClusterWeight[i]=0;
741 void AliITStrackerV2::AliITSlayer::ResetClusters() {
742 //--------------------------------------------------------------------
743 // This function removes loaded clusters
744 //--------------------------------------------------------------------
745 for (Int_t i=0; i<fN; i++) delete fClusters[i];
746 for (Int_t i=0; i<kMaxClusterPerLayer;i++) fClusterWeight[i]=0;
751 void AliITStrackerV2::AliITSlayer::ResetWeights() {
752 //--------------------------------------------------------------------
753 // This function reset weights of the clusters
754 //--------------------------------------------------------------------
755 for (Int_t i=0; i<kMaxClusterPerLayer;i++) fClusterWeight[i]=0;
758 void AliITStrackerV2::AliITSlayer::ResetRoad() {
759 //--------------------------------------------------------------------
760 // This function calculates the road defined by the cluster density
761 //--------------------------------------------------------------------
763 for (Int_t i=0; i<fN; i++) {
764 if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
766 if (n>1) fRoad=2*fR*TMath::Sqrt(3.14/n);
769 Int_t AliITStrackerV2::AliITSlayer::InsertCluster(AliITSclusterV2 *c) {
770 //--------------------------------------------------------------------
771 //This function adds a cluster to this layer
772 //--------------------------------------------------------------------
773 if (fN==kMaxClusterPerLayer) {
774 ::Error("InsertCluster","Too many clusters !\n");
778 if (fN==0) {fClusters[fN++]=c; return 0;}
779 Int_t i=FindClusterIndex(c->GetZ());
780 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliITSclusterV2*));
781 fClusters[i]=c; fN++;
786 Int_t AliITStrackerV2::AliITSlayer::FindClusterIndex(Double_t z) const {
787 //--------------------------------------------------------------------
788 // This function returns the index of the nearest cluster
789 //--------------------------------------------------------------------
791 if (z <= fClusters[0]->GetZ()) return 0;
792 if (z > fClusters[fN-1]->GetZ()) return fN;
793 Int_t b=0, e=fN-1, m=(b+e)/2;
794 for (; b<e; m=(b+e)/2) {
795 if (z > fClusters[m]->GetZ()) b=m+1;
801 void AliITStrackerV2::AliITSlayer::
802 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
803 //--------------------------------------------------------------------
804 // This function sets the "window"
805 //--------------------------------------------------------------------
806 fI=FindClusterIndex(zmin); fZmax=zmax;
807 Double_t circle=2*TMath::Pi()*fR;
808 if (ymax>circle) { ymax-=circle; ymin-=circle; }
809 fYmin=ymin; fYmax=ymax;
812 const AliITSclusterV2 *AliITStrackerV2::AliITSlayer::GetNextCluster(Int_t &ci){
813 //--------------------------------------------------------------------
814 // This function returns clusters within the "window"
815 //--------------------------------------------------------------------
816 const AliITSclusterV2 *cluster=0;
817 for (Int_t i=fI; i<fN; i++) {
818 const AliITSclusterV2 *c=fClusters[i];
819 if (c->GetZ() > fZmax) break;
820 // if (c->IsUsed()) continue;
821 const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
822 Double_t y=fR*det.GetPhi() + c->GetY();
824 if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
825 if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
827 if (y<fYmin) continue;
828 if (y>fYmax) continue;
837 Int_t AliITStrackerV2::AliITSlayer::
838 FindDetectorIndex(Double_t phi, Double_t z) const {
839 //--------------------------------------------------------------------
840 //This function finds the detector crossed by the track
841 //--------------------------------------------------------------------
842 Double_t dphi=-(phi-fPhiOffset);
843 if (dphi < 0) dphi += 2*TMath::Pi();
844 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
845 Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
846 if (np>=fNladders) np-=fNladders;
847 if (np<0) np+=fNladders;
849 Double_t dz=fZOffset-z;
850 Int_t nz=Int_t(dz*(fNdetectors-1)*0.5/fZOffset+0.5);
851 if (nz>=fNdetectors) return -1;
854 return np*fNdetectors + nz;
858 AliITStrackerV2::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
860 //--------------------------------------------------------------------
861 //This function returns the layer thickness at this point (units X0)
862 //--------------------------------------------------------------------
866 if (43<fR&&fR<45) { //SSD2
869 if (TMath::Abs(y-0.00)>3.40) d+=dd;
870 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
871 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
872 for (Int_t i=0; i<12; i++) {
873 if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
874 if (TMath::Abs(y-0.00)>3.40) d+=dd;
878 if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
879 if (TMath::Abs(y-0.00)>3.40) d+=dd;
883 if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
884 if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
887 if (37<fR&&fR<41) { //SSD1
890 if (TMath::Abs(y-0.00)>3.40) d+=dd;
891 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
892 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
893 for (Int_t i=0; i<11; i++) {
894 if (TMath::Abs(z-3.9*i)<0.15) {
895 if (TMath::Abs(y-0.00)>3.40) d+=dd;
899 if (TMath::Abs(z+3.9*i)<0.15) {
900 if (TMath::Abs(y-0.00)>3.40) d+=dd;
904 if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
905 if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
908 if (13<fR&&fR<26) { //SDD
911 if (TMath::Abs(y-0.00)>3.30) d+=dd;
913 if (TMath::Abs(y-1.80)<0.55) {
915 for (Int_t j=0; j<20; j++) {
916 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
917 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
920 if (TMath::Abs(y+1.80)<0.55) {
922 for (Int_t j=0; j<20; j++) {
923 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
924 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
928 for (Int_t i=0; i<4; i++) {
929 if (TMath::Abs(z-7.3*i)<0.60) {
931 if (TMath::Abs(y-0.00)>3.30) d+=dd;
934 if (TMath::Abs(z+7.3*i)<0.60) {
936 if (TMath::Abs(y-0.00)>3.30) d+=dd;
941 if (6<fR&&fR<8) { //SPD2
942 Double_t dd=0.0063; x0=21.5;
944 if (TMath::Abs(y-3.08)>0.5) d+=dd;
945 //if (TMath::Abs(y-3.08)>0.45) d+=dd;
946 if (TMath::Abs(y-3.03)<0.10) {d+=0.014;}
948 if (3<fR&&fR<5) { //SPD1
949 Double_t dd=0.0063; x0=21.5;
951 if (TMath::Abs(y+0.21)>0.6) d+=dd;
952 //if (TMath::Abs(y+0.21)>0.45) d+=dd;
953 if (TMath::Abs(y+0.10)<0.10) {d+=0.014;}
959 Double_t AliITStrackerV2::GetEffectiveThickness(Double_t y,Double_t z) const
961 //--------------------------------------------------------------------
962 //Returns the thickness between the current layer and the vertex (units X0)
963 //--------------------------------------------------------------------
964 Double_t d=0.0028*3*3; //beam pipe
967 Double_t xn=fgLayers[fI].GetR();
968 for (Int_t i=0; i<fI; i++) {
969 Double_t xi=fgLayers[i].GetR();
970 d+=fgLayers[i].GetThickness(y,z,x0)*xi*xi;
979 Double_t xi=0.5*(fgLayers[3].GetR()+fgLayers[4].GetR());
986 Int_t AliITStrackerV2::AliITSlayer::InRoad() const {
987 //--------------------------------------------------------------------
988 // This function returns number of clusters within the "window"
989 //--------------------------------------------------------------------
991 for (Int_t i=fI; i<fN; i++) {
992 const AliITSclusterV2 *c=fClusters[i];
993 if (c->GetZ() > fZmax) break;
994 if (c->IsUsed()) continue;
995 const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
996 Double_t y=fR*det.GetPhi() + c->GetY();
998 if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
999 if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
1001 if (y<fYmin) continue;
1002 if (y>fYmax) continue;
1009 AliITStrackerV2::RefitAt(Double_t xx,AliITStrackV2 *t,const AliITStrackV2 *c) {
1010 //--------------------------------------------------------------------
1011 // This function refits the track "t" at the position "x" using
1012 // the clusters from "c"
1013 //--------------------------------------------------------------------
1014 Int_t index[kMaxLayer];
1016 for (k=0; k<kMaxLayer; k++) index[k]=-1;
1017 Int_t nc=c->GetNumberOfClusters();
1018 for (k=0; k<nc; k++) {
1019 Int_t idx=c->GetClusterIndex(k),nl=(idx&0xf0000000)>>28;
1023 Int_t from, to, step;
1024 if (xx > t->GetX()) {
1025 from=0; to=kMaxLayer;
1028 from=kMaxLayer-1; to=-1;
1032 for (Int_t i=from; i != to; i += step) {
1033 AliITSlayer &layer=fgLayers[i];
1034 Double_t r=layer.GetR();
1037 Double_t hI=i-0.5*step;
1038 if (TMath::Abs(hI-1.5)<0.01 || TMath::Abs(hI-3.5)<0.01) {
1039 Double_t rs=0.5*(fgLayers[i-step].GetR() + r);
1040 Double_t d=0.0034, x0=38.6;
1041 if (TMath::Abs(hI-1.5)<0.01) {rs=9.; d=0.0097; x0=42;}
1042 if (!t->PropagateTo(rs,-step*d,x0)) {
1048 // remember old position [SR, GSI 18.02.2003]
1049 Double_t oldX=0., oldY=0., oldZ=0.;
1050 if (t->IsStartedTimeIntegral() && step==1) {
1051 t->GetGlobalXYZat(t->GetX(),oldX,oldY,oldZ);
1056 if (!t->GetGlobalXYZat(r,x,y,z)) {
1059 Double_t phi=TMath::ATan2(y,x);
1060 Int_t idet=layer.FindDetectorIndex(phi,z);
1064 const AliITSdetector &det=layer.GetDetector(idet);
1066 if (!t->Propagate(phi,det.GetR())) {
1069 t->SetDetectorIndex(idet);
1071 const AliITSclusterV2 *cl=0;
1072 Double_t maxchi2=kMaxChi2;
1076 const AliITSclusterV2 *c=(AliITSclusterV2 *)GetCluster(idx);
1077 if (idet != c->GetDetectorIndex()) {
1078 idet=c->GetDetectorIndex();
1079 const AliITSdetector &det=layer.GetDetector(idet);
1080 if (!t->Propagate(det.GetPhi(),det.GetR())) {
1083 t->SetDetectorIndex(idet);
1085 Double_t chi2=t->GetPredictedChi2(c);
1095 if (t->GetNumberOfClusters()>2) {
1096 Double_t dz=4*TMath::Sqrt(t->GetSigmaZ2()+kSigmaZ2[i]);
1097 Double_t dy=4*TMath::Sqrt(t->GetSigmaY2()+kSigmaY2[i]);
1098 Double_t zmin=t->GetZ() - dz;
1099 Double_t zmax=t->GetZ() + dz;
1100 Double_t ymin=t->GetY() + phi*r - dy;
1101 Double_t ymax=t->GetY() + phi*r + dy;
1102 layer.SelectClusters(zmin,zmax,ymin,ymax);
1104 const AliITSclusterV2 *c=0; Int_t ci=-1;
1105 while ((c=layer.GetNextCluster(ci))!=0) {
1106 if (idet != c->GetDetectorIndex()) continue;
1107 Double_t chi2=t->GetPredictedChi2(c);
1108 if (chi2<maxchi2) { cl=c; maxchi2=chi2; idx=ci; }
1113 if (!t->Update(cl,maxchi2,idx)) {
1116 t->SetSampledEdx(cl->GetQ(),t->GetNumberOfClusters()-1);
1121 Double_t d=layer.GetThickness(t->GetY(),t->GetZ(),x0);
1122 t->CorrectForMaterial(-step*d,x0);
1125 // track time update [SR, GSI 17.02.2003]
1126 if (t->IsStartedTimeIntegral() && step==1) {
1127 Double_t newX, newY, newZ;
1128 t->GetGlobalXYZat(t->GetX(),newX,newY,newZ);
1129 Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) +
1130 (oldZ-newZ)*(oldZ-newZ);
1131 t->AddTimeStep(TMath::Sqrt(dL2));
1137 if (!t->PropagateTo(xx,0.,0.)) return kFALSE;
1142 Float_t *AliITStrackerV2::GetWeight(Int_t index) {
1143 //--------------------------------------------------------------------
1144 // Return pointer to a given cluster
1145 //--------------------------------------------------------------------
1146 Int_t l=(index & 0xf0000000) >> 28;
1147 Int_t c=(index & 0x0fffffff) >> 00;
1148 return fgLayers[l].GetWeight(c);
1153 void AliITStrackerV2::UseClusters(const AliKalmanTrack *t, Int_t from) const {
1154 //--------------------------------------------------------------------
1155 // This function marks clusters assigned to the track
1156 //--------------------------------------------------------------------
1157 AliTracker::UseClusters(t,from);
1159 AliITSclusterV2 *c=(AliITSclusterV2 *)GetCluster(t->GetClusterIndex(0));
1160 //if (c->GetQ()>2) c->Use();
1161 if (c->GetSigmaZ2()>0.1) c->Use();
1162 c=(AliITSclusterV2 *)GetCluster(t->GetClusterIndex(1));
1163 //if (c->GetQ()>2) c->Use();
1164 if (c->GetSigmaZ2()>0.1) c->Use();
1169 void AliITStrackerV2::AddTrackHypothesys(AliITStrackV2 * track, Int_t esdindex)
1171 //------------------------------------------------------------------
1172 // add track to the list of hypothesys
1173 //------------------------------------------------------------------
1175 if (esdindex>=fTrackHypothesys.GetEntriesFast()) fTrackHypothesys.Expand(esdindex*2+10);
1177 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
1179 array = new TObjArray(10);
1180 fTrackHypothesys.AddAt(array,esdindex);
1182 array->AddLast(track);
1185 void AliITStrackerV2::SortTrackHypothesys(Int_t esdindex, Float_t likelihoodlevel)
1187 //-------------------------------------------------------------------
1188 // compress array of track hypothesys
1189 // keep only maxsize best hypothesys
1190 //-------------------------------------------------------------------
1191 if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
1192 if (! (fTrackHypothesys.At(esdindex)) ) return;
1193 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
1194 Int_t entries = array->GetEntriesFast();
1196 Float_t * chi2 = new Float_t[entries];
1197 Float_t * probability = new Float_t[entries];
1198 Float_t sumprobabilityall=0;
1199 Int_t * index = new Int_t[entries];
1202 for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
1203 AliITStrackV2 * track = (AliITStrackV2*)array->At(itrack);
1205 if (track && track->GetNumberOfClusters()>(track->fNUsed+track->fNSkipped)){
1207 chi2[itrack] = track->GetChi2()/(track->GetNumberOfClusters()-track->fNUsed-track->fNSkipped);
1208 if (track->fESDtrack)
1209 if (track->fESDtrack->GetTPCsignal()>80){
1211 if ((track->GetdEdx()/track->fESDtrack->GetTPCsignal())<0.4){
1212 Float_t factor = 2.+10.*(0.6-track->GetdEdx()/track->fESDtrack->GetTPCsignal());
1213 chi2[itrack]*= factor; //mismatch in dEdx
1218 chi2[itrack] = 10000000.;
1219 probability[itrack] = 1./(0.3+chi2[itrack]);
1220 sumprobabilityall+=probability[itrack];
1223 TMath::Sort(entries,chi2,index,kFALSE);
1224 TObjArray * newarray = new TObjArray();
1225 Float_t sumprobability = 0.;
1226 for (Int_t i=0;i<entries;i++){
1227 AliITStrackV2 * track = (AliITStrackV2*)array->At(index[i]);
1229 if (chi2[index[i]]<30){
1230 newarray->AddLast(array->RemoveAt(index[i]));
1231 track->fChi2MIP[0] = chi2[index[i]]; // base chi 2
1232 sumprobability+= probability[index[i]]/sumprobabilityall;
1233 if (sumprobability> likelihoodlevel) break;
1236 delete array->RemoveAt(index[i]);
1241 delete fTrackHypothesys.RemoveAt(esdindex);
1242 fTrackHypothesys.AddAt(newarray,esdindex);
1250 void AliITStrackerV2::CompressTrackHypothesys(Int_t esdindex, Float_t likelihoodlevel, Int_t maxsize)
1254 if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
1255 if (! (fTrackHypothesys.At(esdindex)) ) return;
1256 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
1257 Int_t entries = array->GetEntriesFast();
1259 if (likelihoodlevel>0.000001){
1260 Float_t *probability = new Float_t[entries];
1261 for (Int_t i=0;i<entries;i++) probability[i]=0.;
1262 Float_t sumprobabilityall=0;
1263 Int_t *index = new Int_t[entries];
1265 for (Int_t itrack=0;itrack<entries;itrack++){
1266 AliITStrackV2 * track = (AliITStrackV2*)array->At(itrack);
1267 probability[itrack]=0;
1268 if (!track) continue;
1269 probability[itrack] = 1./(0.3+track->fChi2MIP[0]);
1270 sumprobabilityall += probability[itrack];
1273 if (sumprobabilityall<=0.000000000000001){
1276 for (Int_t itrack=0;itrack<entries;itrack++){
1277 probability[itrack]/=sumprobabilityall;
1280 TMath::Sort(entries, probability, index, kTRUE);
1282 TObjArray * newarray = new TObjArray();
1283 Float_t sumprobability = 0.;
1284 for (Int_t i=0;i<entries;i++){
1285 AliITStrackV2 * track = (AliITStrackV2*)array->At(index[i]);
1286 if (!track) continue;
1287 newarray->AddLast(array->RemoveAt(index[i]));
1288 sumprobability+= probability[index[i]];
1289 if (sumprobability> likelihoodlevel) break;
1290 if (i>maxsize) break;
1294 delete fTrackHypothesys.RemoveAt(esdindex);
1295 fTrackHypothesys.AddAt(newarray,esdindex);
1298 delete []probability;
1302 delete fTrackHypothesys.RemoveAt(esdindex);
1307 AliITStrackV2 * AliITStrackerV2::GetBestHypothesys(Int_t esdindex, AliITStrackV2 * original, Int_t checkmax)
1309 //-------------------------------------------------------------
1310 // try to find best hypothesy
1311 // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
1312 //-------------------------------------------------------------
1313 if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
1314 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
1315 if (!array) return 0;
1316 Int_t entries = array->GetEntriesFast();
1317 if (!entries) return 0;
1318 Float_t minchi2 = 100000;
1320 AliITStrackV2 * besttrack=0;
1327 for (Int_t itrack=0;itrack<entries; itrack++){
1328 AliITStrackV2 * track = (AliITStrackV2*)array->At(itrack);
1329 if (!track) continue;
1331 sumz2+=track->GetSigmaZ2();
1332 sumy2+=track->GetSigmaY2();
1337 Float_t dedxmismatch=1;
1338 for (Int_t i=0;i<entries;i++){
1340 AliITStrackV2 * track = (AliITStrackV2*)array->At(i);
1341 if (!track) continue;
1342 track->fChi2MIP[1] = 1000000;
1343 track->fChi2MIP[2] = 1000000;
1345 if ( (track->GetNumberOfClusters()-track->fNSkipped-track->fNUsed)<2) continue;
1347 if (track->fESDtrack)
1348 if (track->fESDtrack->GetTPCsignal()>80){
1350 if ((track->GetdEdx()/track->fESDtrack->GetTPCsignal())<0.4){
1352 dedxmismatch= 2.+ 10.*(0.6-track->GetdEdx()/track->fESDtrack->GetTPCsignal());
1356 // track->SetLabel(original->GetLabel());
1357 //CookLabel(track,0.0);
1358 //if (track->GetFakeRatio()>0.01) continue;
1362 AliITStrackV2 * backtrack = new AliITStrackV2(*track);
1363 backtrack->ResetCovariance();
1364 backtrack->ResetClusters();
1365 Double_t x = original->GetX();
1366 if (!RefitAt(x,backtrack,track)){
1368 delete array->RemoveAt(i);
1371 if ( (backtrack->GetChi2() / float(backtrack->GetNumberOfClusters()-track->fNSkipped-track->fNUsed-0.5))>6)
1374 delete array->RemoveAt(i);
1377 Double_t deltac = backtrack->GetC()-original->GetC();
1378 Double_t deltatgl = backtrack->GetTgl()-original->GetTgl();
1380 Double_t poolc2 = (deltac*deltac)/(original->fC44+backtrack->fC44);
1381 Double_t pooltgl2 = (deltatgl*deltatgl)/(original->fC33+backtrack->fC33);
1382 if ((poolc2+pooltgl2)>32){ //4 sigma
1384 delete array->RemoveAt(i);
1387 //Double_t bpoolc = (deltac*deltac)/(original->fC44);
1388 //Double_t bpooltgl = (deltatgl*deltatgl)/(original->fC33);
1391 //forward track - without constraint
1392 AliITStrackV2 * forwardtrack = new AliITStrackV2(*original);
1393 // forwardtrack->ResetCovariance();
1394 forwardtrack->ResetClusters();
1396 if (!RefitAt(x,forwardtrack,track)){
1397 delete forwardtrack;
1399 delete array->RemoveAt(i);
1402 if ( (forwardtrack->GetChi2()/float(forwardtrack->GetNumberOfClusters()-track->fNSkipped-track->fNUsed))>6)
1404 delete forwardtrack;
1405 delete array->RemoveAt(i);
1410 if (accepted>checkmax){
1412 delete forwardtrack;
1415 Double_t chi2 = (backtrack->GetChi2()/(backtrack->GetNumberOfClusters()-1-track->fNSkipped-track->fNUsed)+
1416 forwardtrack->GetChi2()/(forwardtrack->GetNumberOfClusters()-track->fNSkipped-track->fNUsed));
1418 // chi2 *= (forwardtrack->GetSigmaZ2()/sumz2+forwardtrack->GetSigmaY2()/sumy2);
1419 chi2 *= dedxmismatch;
1422 track->fChi2MIP[1] = backtrack->GetChi2()/(backtrack->GetNumberOfClusters()-1-track->fNSkipped-track->fNUsed);
1423 track->fChi2MIP[2] = forwardtrack->GetChi2()/(forwardtrack->GetNumberOfClusters()-track->fNSkipped-track->fNUsed);
1424 track->fChi2MIP[3] = poolc2+pooltgl2;
1427 if (track->GetNumberOfClusters()>maxn){
1428 besttrack = new AliITStrackV2(*forwardtrack);
1429 maxn = track->GetNumberOfClusters();
1432 delete forwardtrack;
1436 if (chi2 < minchi2){
1437 besttrack = new AliITStrackV2(*forwardtrack);
1441 delete forwardtrack;
1445 if (!besttrack || besttrack->GetNumberOfClusters()<4) {
1450 besttrack->SetLabel(original->GetLabel());
1451 CookLabel(besttrack,0.0);
1453 // calculate "weight of the cluster"
1457 //sign usage information for clusters
1458 Int_t clusterindex[6][100];
1459 Double_t clusterweight[6][100];
1460 for (Int_t ilayer=0;ilayer<6;ilayer++)
1461 for (Int_t icluster=0;icluster<100;icluster++){
1462 clusterindex[ilayer][icluster] = -1;
1463 clusterweight[ilayer][icluster] = 0;
1465 //printf("%d\t%d\n",esdindex, entries);
1468 for (Int_t itrack=0;itrack<entries; itrack++){
1469 AliITStrackV2 * track = (AliITStrackV2*)array->At(itrack);
1470 if (!track) continue;
1471 if (track->fChi2MIP[1]>1000) continue; //not accepted
1472 sumchi2 +=1./(0.3+track->fChi2MIP[1]+track->fChi2MIP[2]);
1474 for (Int_t itrack=0;itrack<entries;itrack++){
1475 AliITStrackV2 * track = (AliITStrackV2*)array->At(itrack);
1476 if (!track) continue;
1477 if (track->fChi2MIP[1]>1000) continue; //not accepted
1478 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
1479 Int_t tindex = track->GetClusterIndex(icluster);
1480 Int_t ilayer = (tindex & 0xf0000000) >> 28;
1481 if (tindex<0) continue;
1484 for (cindex=0;cindex<100;cindex++){
1485 if (clusterindex[ilayer][cindex]<0) break;
1486 if (clusterindex[ilayer][cindex]==tindex) break;
1488 if (cindex>100) break;
1489 if (clusterindex[ilayer][cindex]!=tindex) clusterindex[ilayer][cindex] = tindex;
1490 clusterweight[ilayer][cindex]+= (1./(0.3+track->fChi2MIP[1]+track->fChi2MIP[2]))* (1./sumchi2);
1491 Float_t *weight = GetWeight(tindex);
1494 *weight+= (1./(0.3+track->fChi2MIP[1]+track->fChi2MIP[2]))* (1./sumchi2);
1499 if (besttrack->GetChi2()/besttrack->GetNumberOfClusters()>3.5) return besttrack; //don't sign clusters
1501 Double_t deltad = besttrack->GetD(GetX(),GetY());
1502 Double_t deltaz = besttrack->GetZat(GetX()) - GetZ();
1503 Double_t deltaprim = TMath::Sqrt(deltad*deltad+deltaz*deltaz);
1505 for (Int_t icluster=0;icluster<besttrack->GetNumberOfClusters();icluster++){
1506 Int_t index = besttrack->GetClusterIndex(icluster);
1507 Int_t ilayer = (index & 0xf0000000) >> 28;
1508 AliITSclusterV2 *c = (AliITSclusterV2*)GetCluster(index);
1511 for (Int_t icluster=0;icluster<100;icluster++){
1512 // sign non "doubt" clusters as used
1513 if (clusterindex[ilayer][icluster]!=index) continue;
1514 // Float_t * weight = GetWeight(index);
1515 //if (weight) if (*weight>1){
1516 // if (c->IsUsed()) continue;
1519 if ( (ilayer*0.2+0.2)<deltaprim) continue; // secondaries
1520 if (c->GetNy()>4) continue; // don sign cluster
1521 if ( (ilayer>1&&clusterweight[ilayer][icluster]>0.7) || (ilayer<2&&clusterweight[ilayer][icluster]>0.8) ){
1523 if (c->IsUsed()) continue;
1535 AliITStrackV2 * AliITStrackerV2::GetBestHypothesysMIP(Int_t esdindex, AliITStrackV2 * original)
1537 //-------------------------------------------------------------
1538 // try to find best hypothesy
1539 // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
1540 //-------------------------------------------------------------
1541 if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
1542 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
1543 if (!array) return 0;
1544 Int_t entries = array->GetEntriesFast();
1545 if (!entries) return 0;
1546 AliITStrackV2 * besttrack=0;
1548 //sign usage information for clusters
1549 Int_t clusterindex[6][100];
1550 Double_t clusterweight[6][100];
1551 for (Int_t ilayer=0;ilayer<6;ilayer++)
1552 for (Int_t icluster=0;icluster<100;icluster++){
1553 clusterindex[ilayer][icluster] = -1;
1554 clusterweight[ilayer][icluster] = 0;
1556 //printf("%d\t%d\n",esdindex, entries);
1559 for (Int_t itrack=0;itrack<entries; itrack++){
1560 AliITStrackV2 * track = (AliITStrackV2*)array->At(itrack);
1561 if (!track) continue;
1562 if (track->fChi2MIP[1]>1000) continue; //not accepted - before
1563 sumchi2 +=1./(0.3+track->fChi2MIP[1]+track->fChi2MIP[2]);
1566 // get cluster weight
1567 for (Int_t itrack=0;itrack<entries;itrack++){
1568 AliITStrackV2 * track = (AliITStrackV2*)array->At(itrack);
1569 if (!track) continue;
1570 if (track->fChi2MIP[1]>1000) continue; // track not accepted in previous iterration
1571 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
1572 Int_t tindex = track->GetClusterIndex(icluster);
1573 Int_t ilayer = (tindex & 0xf0000000) >> 28;
1574 if (tindex<0) continue;
1577 for (cindex=0;cindex<100;cindex++){
1578 if (clusterindex[ilayer][cindex]<0) break;
1579 if (clusterindex[ilayer][cindex]==tindex) break;
1581 if (cindex>100) break;
1582 if (clusterindex[ilayer][cindex]!=tindex) clusterindex[ilayer][cindex] = tindex;
1583 clusterweight[ilayer][cindex]+= (1./(0.3+track->fChi2MIP[1]+track->fChi2MIP[2]))* (1./sumchi2);
1587 // get cluster relative sharing - factor
1590 for (Int_t ilayer=0;ilayer<6;ilayer++)
1591 for (Int_t cindex=0;cindex<100;cindex++){
1592 if (clusterindex[ilayer][cindex]<0) continue;
1593 Int_t tindex = clusterindex[ilayer][cindex];
1594 Float_t *weight = GetWeight(tindex);
1596 printf("Problem 1\n"); // not existing track
1599 if (*weight<(clusterweight[ilayer][cindex]-0.00001)){
1600 printf("Problem 2\n"); // not normalized probability
1603 AliITSclusterV2 *c = (AliITSclusterV2*)GetCluster(tindex);
1605 clusterweight[ilayer][cindex]/= *weight;
1608 Float_t weight2 = TMath::Max(*weight-0.5*clusterweight[ilayer][cindex],0.0000001);
1609 clusterweight[ilayer][cindex]/= weight2;
1610 if ( clusterweight[ilayer][cindex]>1) clusterweight[ilayer][cindex] =1.;
1614 //take to the account sharing factor
1616 Float_t chi2 =10000000;
1617 Float_t sharefactor=0;
1618 Float_t minchi2 = 100000000;
1619 Float_t secchi2 = 100000000;
1621 for (Int_t itrack=0;itrack<entries; itrack++){
1622 AliITStrackV2 * track = (AliITStrackV2*)array->At(itrack);
1623 if (!track) continue;
1624 if (track->fChi2MIP[1]>1000) continue; //not accepted - before
1625 chi2 = track->fChi2MIP[1];
1630 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
1631 Int_t tindex = track->GetClusterIndex(icluster);
1632 Int_t ilayer = (tindex & 0xf0000000) >> 28;
1633 if (tindex<0) continue;
1635 Float_t cchi2 = (track->fDy[ilayer]*track->fDy[ilayer])/(track->fSigmaY[ilayer]*track->fSigmaY[ilayer]) +
1636 (track->fDz[ilayer]*track->fDz[ilayer])/(track->fSigmaZ[ilayer]*track->fSigmaZ[ilayer]) ;
1638 for (cindex=0;cindex<100;cindex++){
1639 if (clusterindex[ilayer][cindex]<0) break;
1640 if (clusterindex[ilayer][cindex]==tindex) break;
1642 if (cindex>100) continue;
1643 if (clusterweight[ilayer][cindex]>0.00001){
1644 sharefactor+= clusterweight[ilayer][cindex];
1645 cchi2/=clusterweight[ilayer][cindex];
1650 newchi2/=(norm-track->fNSkipped-track->fNUsed);
1651 track->fChi2MIP[4] = newchi2;
1652 //if (norm>0) sharefactor/=norm;
1653 //if (sharefactor<0.5) return 0;
1654 // chi2*=1./(0.5+sharefactor);
1655 if (newchi2<minchi2){
1660 if (newchi2<secchi2){
1669 if (!besttrack || besttrack->GetNumberOfClusters()<4) {
1673 if ((minchi2/secchi2)<0.7){
1675 //increase the weight for clusters if the probability of other hypothesys is small
1676 Double_t deltad = besttrack->GetD(GetX(),GetY());
1677 Double_t deltaz = besttrack->GetZat(GetX()) - GetZ();
1678 Double_t deltaprim = TMath::Sqrt(deltad*deltad+deltaz*deltaz);
1680 for (Int_t icluster=0;icluster<besttrack->GetNumberOfClusters();icluster++){
1681 Int_t index = besttrack->GetClusterIndex(icluster);
1682 Int_t ilayer = (index & 0xf0000000) >> 28;
1683 AliITSclusterV2 *c = (AliITSclusterV2*)GetCluster(index);
1685 if (c->GetNy()>3) continue;
1686 if ( (ilayer*0.2+0.2)<deltaprim) continue; // secondaries
1687 Float_t * weight = GetWeight(index);
1694 besttrack->SetLabel(original->GetLabel());
1695 CookLabel(besttrack,0.0);