X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITStrackerV2.cxx;h=a1b1d6a2c5099c77776df72b390c4d00b5272709;hb=5a19f42ab781b5b5df08209adb11fe19cf597920;hp=6d6a5335cf1856de8d1b889b1874ad61e5b354ee;hpb=73f446a023af9c141bc518face2989dcf0dfd36e;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/AliITStrackerV2.cxx b/ITS/AliITStrackerV2.cxx index 6d6a5335cf1..a1b1d6a2c50 100644 --- a/ITS/AliITStrackerV2.cxx +++ b/ITS/AliITStrackerV2.cxx @@ -12,10 +12,10 @@ * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ - +/* $Id$ */ //------------------------------------------------------------------------- // Implementation of the ITS tracker class -// It reads AliITSclusterV2 clusters and creates AliITStrackV2 tracks +// It reads AliITSRecPoint clusters and creates AliITStrackV2 tracks // and fills with them the ESD // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch // dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch @@ -23,42 +23,113 @@ #include +#include #include #include #include +#include -#include "AliITSgeom.h" +#include "AliITSgeomTGeo.h" +#include "AliAlignObj.h" +#include "AliITSRecPoint.h" +#include "AliESDEvent.h" +#include "AliESDtrack.h" #include "AliITSRecPoint.h" -#include "AliESD.h" -#include "AliITSclusterV2.h" +#include "AliITSReconstructor.h" #include "AliITStrackerV2.h" ClassImp(AliITStrackerV2) -AliITStrackerV2::AliITSlayer AliITStrackerV2::fgLayers[kMaxLayer]; // ITS layers +AliITStrackerV2::AliITSlayer AliITStrackerV2::fgLayers[AliITSgeomTGeo::kNLayers]; //ITS layers + +AliITStrackerV2::AliITStrackerV2(): + AliTracker(), + fI(AliITSgeomTGeo::GetNLayers()), + fBestTrack(), + fTrackToFollow(), + fPass(0), + fLastLayerToTrackTo(AliITSRecoParam::GetLastLayerToTrackTo()) +{ + //-------------------------------------------------------------------- + //This is the AliITStrackerV2 default constructor + //-------------------------------------------------------------------- + + for (Int_t i=1; iGetXVdef(), + AliITSReconstructor::GetRecoParam()->GetYVdef(), + AliITSReconstructor::GetRecoParam()->GetZVdef()}; + Double_t ers[]={AliITSReconstructor::GetRecoParam()->GetSigmaXVdef(), + AliITSReconstructor::GetRecoParam()->GetSigmaYVdef(), + AliITSReconstructor::GetRecoParam()->GetSigmaZVdef()}; + SetVertex(xyz,ers); + + for (Int_t i=0; iGetXVdef(), + AliITSReconstructor::GetRecoParam()->GetYVdef(), + AliITSReconstructor::GetRecoParam()->GetZVdef()}; + Double_t ers[]={AliITSReconstructor::GetRecoParam()->GetSigmaXVdef(), + AliITSReconstructor::GetRecoParam()->GetSigmaYVdef(), + AliITSReconstructor::GetRecoParam()->GetSigmaZVdef()}; + xyz[0]=t.GetX(); xyz[1]=t.GetY(); xyz[2]=t.GetZ(); + ers[0]=t.GetSigmaX(); ers[1]=t.GetSigmaY(); ers[2]=t.GetSigmaZ(); + SetVertex(xyz,ers); + + for (Int_t i=0; iGetNladders(i); - Int_t ndet=g->GetNdetectors(i); + for (Int_t i=1; iGetTrans(i,1,1,x,y,z); - Double_t r=TMath::Sqrt(x*x + y*y); + Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2]; + AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz); Double_t poff=TMath::ATan2(y,x); Double_t zoff=z; + Double_t r=TMath::Sqrt(x*x + y*y); - g->GetTrans(i,1,2,x,y,z); + AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz); r += TMath::Sqrt(x*x + y*y); - g->GetTrans(i,2,1,x,y,z); + AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz); r += TMath::Sqrt(x*x + y*y); - g->GetTrans(i,2,2,x,y,z); + AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz); r += TMath::Sqrt(x*x + y*y); r*=0.25; @@ -66,14 +137,17 @@ AliITStrackerV2::AliITStrackerV2(const AliITSgeom *geom) : AliTracker() { for (Int_t j=1; jGetTrans(i,j,k,x,y,zshift); - Double_t rot[9]; g->GetRotMatrix(i,j,k,rot); - - Double_t phi=TMath::ATan2(rot[1],rot[0])+TMath::Pi(); - phi+=TMath::Pi()/2; - if (i==1) phi+=TMath::Pi(); - Double_t cp=TMath::Cos(phi), sp=TMath::Sin(phi); - Double_t r=x*cp+y*sp; + TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m); + const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k); + m.Multiply(tm); + Double_t txyz[3]={0.}; + xyz[0]=0.; xyz[1]=0.; xyz[2]=0.; + m.LocalToMaster(txyz,xyz); + r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]); + Double_t phi=TMath::ATan2(xyz[1],xyz[0]); + + if (phi<0) phi+=TMath::TwoPi(); + else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi(); AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1); new(&det) AliITSdetector(r,phi); @@ -82,16 +156,17 @@ AliITStrackerV2::AliITStrackerV2(const AliITSgeom *geom) : AliTracker() { } - fI=kMaxLayer; - - fPass=0; fConstraint[0]=1; fConstraint[1]=0; - Double_t xyz[]={kXV,kYV,kZV}, ers[]={kSigmaXV,kSigmaYV,kSigmaZV}; + Double_t xyz[]={AliITSReconstructor::GetRecoParam()->GetXVdef(), + AliITSReconstructor::GetRecoParam()->GetYVdef(), + AliITSReconstructor::GetRecoParam()->GetZVdef()}; + Double_t ers[]={AliITSReconstructor::GetRecoParam()->GetSigmaXVdef(), + AliITSReconstructor::GetRecoParam()->GetSigmaYVdef(), + AliITSReconstructor::GetRecoParam()->GetSigmaZVdef()}; SetVertex(xyz,ers); - for (Int_t i=0; iGetBranch("Clusters"); + TBranch *branch=cTree->GetBranch("ITSRecPoints"); if (!branch) { Error("LoadClusters"," can't get the branch !\n"); return 1; } - TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy; + TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy; branch->SetAddress(&clusters); Int_t j=0; - Int_t detector=0; - for (Int_t i=0; iGetEvent(j)) continue; Int_t ncl=clusters->GetEntriesFast(); + while (ncl--) { - AliITSclusterV2 *c=(AliITSclusterV2*)clusters->UncheckedAt(ncl); - detector = c->GetDetectorIndex(); - fgLayers[i].InsertCluster(new AliITSclusterV2(*c)); + AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl); + + if (!c->Misalign()) AliWarning("Can't misalign this cluster !"); + + Int_t idx=c->GetDetectorIndex(); + AliITSdetector &det=fgLayers[i].GetDetector(idx); + + Double_t y=r*det.GetPhi()+c->GetY(); + if (y>circ) y-=circ; else if (y<0) y+=circ; + c->SetPhiR(y); + + fgLayers[i].InsertCluster(new AliITSRecPoint(*c)); } clusters->Delete(); - //add dead zone virtual "cluster" - - if (i<1){ - for (Float_t ydead = -2.; ydead < 2. ; ydead+=0.05){ - Int_t lab[4] = {0,0,0,detector}; - Int_t info[3] = {0,0,0}; - Float_t hit[5]={ydead,0,1,0.01,0}; - fgLayers[i].InsertCluster(new AliITSclusterV2(lab, hit, info)); - // - hit[1]=-7.; - fgLayers[i].InsertCluster(new AliITSclusterV2(lab, hit, info)); - } - } - } - // fgLayers[i].ResetRoad(); //road defined by the cluster density } @@ -155,7 +228,7 @@ void AliITStrackerV2::UnloadClusters() { //-------------------------------------------------------------------- //This function unloads ITS clusters //-------------------------------------------------------------------- - for (Int_t i=0; iGetD())>5) { + if (TMath::Abs(t->GetD(GetX(),GetY()))>4) { delete t; continue; } @@ -228,115 +301,54 @@ Int_t AliITStrackerV2::Clusters2Tracks(AliESD *event) { delete t; continue; } - t->SetReconstructed(kFALSE); itsTracks.AddLast(t); } } /* End Read ESD tracks */ itsTracks.Sort(); Int_t nentr=itsTracks.GetEntriesFast(); - fTrackHypothesys.Expand(nentr); + Int_t ntrk=0; for (fPass=0; fPass<2; fPass++) { Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue; for (Int_t i=0; iGetReconstructed()) continue; //this track was already "succesfully" reconstructed - if ( (TMath::Abs(t->GetD(GetX(),GetY())) >2.) && fConstraint[fPass]) continue; - if ( (TMath::Abs(t->GetZat(GetX())-GetZ())>2.) && fConstraint[fPass]) continue; - + if (t==0) continue; //this track has been already tracked Int_t tpcLabel=t->GetLabel(); //save the TPC track label ResetTrackToFollow(*t); ResetBestTrack(); - for (FollowProlongation(); fISetLabel(tpcLabel); - besttrack->CookdEdx(); - besttrack->SetFakeRatio(1.); - CookLabel(besttrack,0.); //For comparison only - // besttrack->UpdateESDtrack(AliESDtrack::kITSin); - - if (besttrack->GetChi2()/besttrack->GetNumberOfClusters()>3.5){ - if ( (TMath::Abs(besttrack->GetD(GetX(),GetY()))>0.4) && fConstraint[fPass]) { - CompressTrackHypothesys(fCurrentEsdTrack,0.0,0); - continue; - } - if ( (TMath::Abs(besttrack->GetZat(GetX()) -GetZ() )>0.4) && fConstraint[fPass]){ - CompressTrackHypothesys(fCurrentEsdTrack,0.0,0); - continue; - } - } - - //delete itsTracks.RemoveAt(i); - t->SetReconstructed(); - ntrk++; - } } - - for (Int_t i=0; iGetLabel(); //save the TPC track label - AliITStrackV2 * besttrack = GetBestHypothesysMIP(fCurrentEsdTrack,t); - if (!besttrack) continue; - - besttrack->SetLabel(tpcLabel); - besttrack->CookdEdx(); - besttrack->SetFakeRatio(1.); - CookLabel(besttrack,0.); //For comparison only - besttrack->UpdateESDtrack(AliESDtrack::kITSin); - } - // itsTracks.Delete(); - // - Int_t entries = fTrackHypothesys.GetEntriesFast(); - for (Int_t ientry=0;ientryDelete(); - delete fTrackHypothesys.RemoveAt(ientry); - } - fTrackHypothesys.Delete(); Info("Clusters2Tracks","Number of prolonged tracks: %d\n",ntrk); return 0; } - -Int_t AliITStrackerV2::PropagateBack(AliESD *event) { +Int_t AliITStrackerV2::PropagateBack(AliESDEvent *event) { //-------------------------------------------------------------------- // This functions propagates reconstructed ITS tracks back // The clusters must be loaded ! @@ -365,13 +377,13 @@ Int_t AliITStrackerV2::PropagateBack(AliESD *event) { // propagete to vertex [SR, GSI 17.02.2003] // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov if (fTrackToFollow.PropagateTo(3.,0.0028,65.19)) { - if (fTrackToFollow.PropagateToVertex()) { + if (fTrackToFollow.PropagateToVertex(event->GetVertex())) { fTrackToFollow.StartTimeIntegral(); } fTrackToFollow.PropagateTo(3.,-0.0028,65.19); } - fTrackToFollow.ResetCovariance(); fTrackToFollow.ResetClusters(); + fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters(); if (RefitAt(49.,&fTrackToFollow,t)) { if (CorrectForDeadZoneMaterial(&fTrackToFollow)!=0) { Warning("PropagateBack", @@ -383,7 +395,7 @@ Int_t AliITStrackerV2::PropagateBack(AliESD *event) { //fTrackToFollow.CookdEdx(); CookLabel(&fTrackToFollow,0.); //For comparison only fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout); - //UseClusters(&fTrackToFollow); + UseClusters(&fTrackToFollow); ntrk++; } delete t; @@ -394,7 +406,7 @@ Int_t AliITStrackerV2::PropagateBack(AliESD *event) { return 0; } -Int_t AliITStrackerV2::RefitInward(AliESD *event) { +Int_t AliITStrackerV2::RefitInward(AliESDEvent *event) { //-------------------------------------------------------------------- // This functions refits ITS tracks using the // "inward propagated" TPC tracks @@ -410,7 +422,7 @@ Int_t AliITStrackerV2::RefitInward(AliESD *event) { if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue; if (esd->GetStatus()&AliESDtrack::kITSrefit) continue; if (esd->GetStatus()&AliESDtrack::kTPCout) - if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue; + if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue; AliITStrackV2 *t=0; try { @@ -431,44 +443,19 @@ Int_t AliITStrackerV2::RefitInward(AliESD *event) { ResetTrackToFollow(*t); fTrackToFollow.ResetClusters(); - if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) - fTrackToFollow.ResetCovariance(); - //Refitting... - if (RefitAt(3.7, &fTrackToFollow, t)) { + if (RefitAt(3.7, &fTrackToFollow, t, kTRUE)) { fTrackToFollow.SetLabel(t->GetLabel()); fTrackToFollow.CookdEdx(); - CookLabel(&fTrackToFollow,0.0); //For comparison only - - if (fTrackToFollow.PropagateTo(3.,0.0028,65.19)) {//The beam pipe - Double_t a=fTrackToFollow.GetAlpha(); - Double_t cs=TMath::Cos(a),sn=TMath::Sin(a); - Double_t xv= GetX()*cs + GetY()*sn; - Double_t yv=-GetX()*sn + GetY()*cs; - - Double_t c=fTrackToFollow.GetC(), snp=fTrackToFollow.GetSnp(); - Double_t x=fTrackToFollow.GetX(), y=fTrackToFollow.GetY(); - Double_t tgfv=-(c*(x-xv)-snp)/(c*(y-yv) + TMath::Sqrt(1.-snp*snp)); - Double_t fv=TMath::ATan(tgfv); - - cs=TMath::Cos(fv); sn=TMath::Sin(fv); - x = xv*cs + yv*sn; - yv=-xv*sn + yv*cs; xv=x; - - if (fTrackToFollow.Propagate(fv+a,xv)) { - fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit); - //UseClusters(&fTrackToFollow); - { - AliITSclusterV2 c; c.SetY(yv); c.SetZ(GetZ()); - c.SetSigmaY2(GetSigmaY()*GetSigmaY()); - c.SetSigmaZ2(GetSigmaZ()*GetSigmaZ()); - Double_t chi2=fTrackToFollow.GetPredictedChi2(&c); - if (chi2RelateToVertex(event->GetVertex(),GetBz(r),maxD); + ntrk++; } } delete t; @@ -512,12 +499,11 @@ void AliITStrackerV2::FollowProlongation() { } //find intersection - Double_t x,y,z; - if (!fTrackToFollow.GetGlobalXYZat(r,x,y,z)) { + Double_t phi,z; + if (!fTrackToFollow.GetPhiZat(r,phi,z)) { //Warning("FollowProlongation","failed to estimate track !\n"); return; } - Double_t phi=TMath::ATan2(y,x); Int_t idet=layer.FindDetectorIndex(phi,z); if (idet<0) { @@ -536,85 +522,57 @@ void AliITStrackerV2::FollowProlongation() { //Select possible prolongations and store the current track estimation track.~AliITStrackV2(); new(&track) AliITStrackV2(fTrackToFollow); - Double_t dz=7*TMath::Sqrt(track.GetSigmaZ2() + kSigmaZ2[i]); - Double_t dy=7*TMath::Sqrt(track.GetSigmaY2() + kSigmaY2[i]); + Double_t dz=7*TMath::Sqrt(track.GetSigmaZ2() + AliITSReconstructor::GetRecoParam()->GetSigmaZ2(i)); + Double_t dy=7*TMath::Sqrt(track.GetSigmaY2() + AliITSReconstructor::GetRecoParam()->GetSigmaY2(i)); Double_t road=layer.GetRoad(); if (dz*dy>road*road) { Double_t dd=TMath::Sqrt(dz*dy), scz=dz/dd, scy=dy/dd; dz=road*scz; dy=road*scy; } - //Double_t dz=4*TMath::Sqrt(track.GetSigmaZ2() + kSigmaZ2[i]); + //Double_t dz=4*TMath::Sqrt(track.GetSigmaZ2() + AliITSReconstructor::GetRecoParam()->GetSigmaZ2(i)); if (dz < 0.5*TMath::Abs(track.GetTgl())) dz=0.5*TMath::Abs(track.GetTgl()); - if (dz > kMaxRoad) { + if (dz > AliITSReconstructor::GetRecoParam()->GetMaxRoad()) { //Warning("FollowProlongation","too broad road in Z !\n"); return; } - // if (TMath::Abs(fTrackToFollow.GetZ()-GetZ()) > r+dz) return; + if (TMath::Abs(fTrackToFollow.GetZ()-GetZ()) > r+dz) return; - //Double_t dy=4*TMath::Sqrt(track.GetSigmaY2() + kSigmaY2[i]); - if (TMath::Abs(track.GetSnp()>kMaxSnp)) { - fI--; - continue; // MI - } + //Double_t dy=4*TMath::Sqrt(track.GetSigmaY2() + AliITSReconstructor::GetRecoParam()->GetSigmaY2(i)); if (dy < 0.5*TMath::Abs(track.GetSnp())) dy=0.5*TMath::Abs(track.GetSnp()); - if (dy > kMaxRoad) { + if (dy > AliITSReconstructor::GetRecoParam()->GetMaxRoad()) { //Warning("FollowProlongation","too broad road in Y !\n"); return; } + fI--; + Double_t zmin=track.GetZ() - dz; Double_t zmax=track.GetZ() + dz; Double_t ymin=track.GetY() + r*phi - dy; Double_t ymax=track.GetY() + r*phi + dy; - layer.SelectClusters(zmin,zmax,ymin,ymax); - fI--; + if (layer.SelectClusters(zmin,zmax,ymin,ymax)==0) + if (fLayersNotToSkip[fI]) return; + + if (!TakeNextProlongation()) + if (fLayersNotToSkip[fI]) return; - //take another prolongation - if (!TakeNextProlongation()){ - //skipped++; - fTrackToFollow.IncrementNSkipped(); - if (fLayersNotToSkip[fI]||fTrackToFollow.GetNSkipped()>1) return; - } - if (fTrackToFollow.GetNUsed()>1) return; - if (fTrackToFollow.GetNUsed()+fTrackToFollow.GetNSkipped()>1) return; - if ( (fI<3) && ( fTrackToFollow.GetChi2()/fTrackToFollow.GetNumberOfClusters()>kChi2PerCluster)) return; } //deal with the best track Int_t ncl=fTrackToFollow.GetNumberOfClusters(); Int_t nclb=fBestTrack.GetNumberOfClusters(); - if (ncl){ - if (ncl<4) return; - if ( (ncl<6) && (fTrackToFollow.GetChi2()/float(ncl))>3) return; - if (fTrackToFollow.GetChi2()/ncl>5.5) return; - fTrackToFollow.CookdEdx(); - if (fTrackToFollow.GetESDtrack()->GetTPCsignal()>80.) - if ((fTrackToFollow.GetdEdx()/fTrackToFollow.GetESDtrack()->GetTPCsignal())<0.35){ - // mismatch in dedx - return; - } - // - fTrackToFollow.SetLabel(fTrackToFollow.GetESDtrack()->GetLabel()); - CookLabel(&fTrackToFollow,0.); // - // - if ( (nclb>3) && ((fTrackToFollow.GetChi2()/ncl)<(3*fBestTrack.GetChi2()/(nclb)))) - AddTrackHypothesys(new AliITStrackV2(fTrackToFollow), fCurrentEsdTrack); - else - if (ncl>3) AddTrackHypothesys(new AliITStrackV2(fTrackToFollow), fCurrentEsdTrack); - if (ncl >= nclb) { - Double_t chi2=fTrackToFollow.GetChi2(); - if (chi2/ncl < kChi2PerCluster) { - if (ncl > nclb ) { - ResetBestTrack(); + if (ncl) + if (ncl >= nclb) { + Double_t chi2=fTrackToFollow.GetChi2(); + if (chi2/ncl < AliITSReconstructor::GetRecoParam()->GetChi2PerCluster()) { + if (ncl > nclb || chi2 < fBestTrack.GetChi2()) { + ResetBestTrack(); } - if ( (ncl == nclb) && chi2 < fBestTrack.GetChi2()) { - ResetBestTrack(); - } - } - } + } } + } Int_t AliITStrackerV2::TakeNextProlongation() { @@ -626,18 +584,20 @@ Int_t AliITStrackerV2::TakeNextProlongation() { AliITSlayer &layer=fgLayers[fI]; ResetTrackToFollow(fTracks[fI]); - Double_t dz=7*TMath::Sqrt(fTrackToFollow.GetSigmaZ2() + kSigmaZ2[fI]); - Double_t dy=7*TMath::Sqrt(fTrackToFollow.GetSigmaY2() + kSigmaY2[fI]); + Double_t dz=7*TMath::Sqrt(fTrackToFollow.GetSigmaZ2() + AliITSReconstructor::GetRecoParam()->GetSigmaZ2(fI)); + Double_t dy=7*TMath::Sqrt(fTrackToFollow.GetSigmaY2() + AliITSReconstructor::GetRecoParam()->GetSigmaY2(fI)); Double_t road=layer.GetRoad(); if (dz*dy>road*road) { Double_t dd=TMath::Sqrt(dz*dy), scz=dz/dd, scy=dy/dd; dz=road*scz; dy=road*scy; } - const AliITSclusterV2 *c=0; Int_t ci=-1; - Double_t chi2=12345.; + const AliITSRecPoint *c=0; Int_t ci=-1; + const AliITSRecPoint *cc=0; Int_t cci=-1; + Double_t chi2=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); while ((c=layer.GetNextCluster(ci))!=0) { Int_t idet=c->GetDetectorIndex(); + if (fTrackToFollow.GetDetectorIndex()!=idet) { const AliITSdetector &det=layer.GetDetector(idet); ResetTrackToFollow(fTracks[fI]); @@ -652,33 +612,29 @@ Int_t AliITStrackerV2::TakeNextProlongation() { if (TMath::Abs(fTrackToFollow.GetZ() - c->GetZ()) > dz) continue; if (TMath::Abs(fTrackToFollow.GetY() - c->GetY()) > dy) continue; - chi2=fTrackToFollow.GetPredictedChi2(c); if (chi2 chi2) continue; + chi2=ch2; + cc=c; cci=ci; + break; } - if (chi2>=kMaxChi2) return 0; - if (!c) return 0; - if (c->IsUsed()&&c->GetNy()<5) { //shared factor - chi2+=1; - chi2*=2*(1./(TMath::Max(c->GetNy(),1))); - } - if (c->GetQ()==0){ //dead zone virtual cluster - chi2*=4.; - chi2+=20; - fTrackToFollow.IncrementNUsed(); - return 1; - } - //if ((fI<2)&&chi2>kMaxChi2In) return 0; + if (!cc) return 0; - if (!fTrackToFollow.Update(c,chi2,(fI<<28)+ci)) { + {// Take into account the mis-alignment + Double_t x = fTrackToFollow.GetX() + cc->GetX(); + if (!fTrackToFollow.PropagateTo(x,0.,0.)) return 0; + } + if (!fTrackToFollow.Update(cc,chi2,(fI<<28)+cci)) { //Warning("TakeNextProlongation","filtering failed !\n"); return 0; } - if (c->IsUsed()&&c->GetNy()<5) fTrackToFollow.IncrementNUsed(); + if (fTrackToFollow.GetNumberOfClusters()>1) - if (TMath::Abs(fTrackToFollow.GetD())>4) return 0; + if (TMath::Abs(fTrackToFollow.GetD(GetX(),GetY()))>4) return 0; fTrackToFollow. - SetSampledEdx(c->GetQ(),fTrackToFollow.GetNumberOfClusters()-1); //b.b. + SetSampledEdx(cc->GetQ(),fI-2); //b.b. { Double_t x0; @@ -690,15 +646,6 @@ Int_t AliITStrackerV2::TakeNextProlongation() { Double_t d=GetEffectiveThickness(0,0); //Think of this !!!! Double_t xyz[]={GetX(),GetY(),GetZ()}; Double_t ers[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()}; - Double_t deltad = TMath::Abs(fTrackToFollow.GetD(GetX(),GetY())); - Double_t deltaz = TMath::Abs(fTrackToFollow.GetZat(GetX())-GetZ()); - - if ( (fI==4) && (deltad>2.0 || deltaz>1.5)) return 0; // don't "improve" secondaries - if ( (fI==3) && (deltad>1.5 || deltaz>0.9)) return 0; // don't "improve" secondaries - if ( (fI==2) && (deltad>0.9 || deltaz>0.6)) return 1; // don't "improve" secondaries - if ( (fI==1) && (deltad>0.3 || deltaz>0.3)) return 1; // don't "improve" secondaries - if ( (fI==0) && (deltad>0.1 || deltaz>0.1)) return 1; // don't "improve" secondaries - fTrackToFollow.Improve(d,xyz,ers); } @@ -706,27 +653,43 @@ Int_t AliITStrackerV2::TakeNextProlongation() { } -AliITStrackerV2::AliITSlayer::AliITSlayer() { +AliITStrackerV2::AliITSlayer::AliITSlayer(): + fR(0.), + fPhiOffset(0.), + fNladders(0), + fZOffset(0.), + fNdetectors(0), + fDetectors(0), + fNsel(0), + fRoad(2*fR*TMath::Sqrt(3.14/1.)) //assuming that there's only one cluster +{ //-------------------------------------------------------------------- //default AliITSlayer constructor //-------------------------------------------------------------------- - fN=0; - fDetectors=0; + + for (Int_t i=0; iGetZ())GetZ())1) fRoad=2*fR*TMath::Sqrt(3.14/n); } -Int_t AliITStrackerV2::AliITSlayer::InsertCluster(AliITSclusterV2 *c) { +Int_t AliITStrackerV2::AliITSlayer::InsertCluster(AliITSRecPoint *c) { //-------------------------------------------------------------------- - //This function adds a cluster to this layer + // This function inserts a cluster to this layer in increasing + // order of the cluster's fZ //-------------------------------------------------------------------- - if (fN==kMaxClusterPerLayer) { - ::Error("InsertCluster","Too many clusters !\n"); - return 1; + Float_t circ=TMath::TwoPi()*fR; + Int_t sec=Int_t(kNsector*c->GetPhiR()/circ); + if (sec>=kNsector) { + ::Error("InsertCluster","Wrong sector !\n"); + return 1; } - - if (fN==0) {fClusters[fN++]=c; return 0;} - Int_t i=FindClusterIndex(c->GetZ()); - memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliITSclusterV2*)); - fClusters[i]=c; fN++; - + Int_t &n=fN[sec]; + if (n>=kMaxClusterPerSector) { + ::Error("InsertCluster","Too many clusters !\n"); + return 1; + } + if (n==0) fClusters[sec*kMaxClusterPerSector]=c; + else { + Int_t i=FindClusterIndex(c->GetZ(),sec); + Int_t k=n-i+sec*kMaxClusterPerSector; + memmove(fClusters+i+1 ,fClusters+i,k*sizeof(AliITSRecPoint*)); + fClusters[i]=c; + } + n++; return 0; } -Int_t AliITStrackerV2::AliITSlayer::FindClusterIndex(Double_t z) const { +Int_t +AliITStrackerV2::AliITSlayer::FindClusterIndex(Float_t z,Int_t s) const { //-------------------------------------------------------------------- - // This function returns the index of the nearest cluster + // For the sector "s", this function returns the index of the first + // with its fZ >= "z". //-------------------------------------------------------------------- - if (fN==0) return 0; - if (z <= fClusters[0]->GetZ()) return 0; - if (z > fClusters[fN-1]->GetZ()) return fN; - Int_t b=0, e=fN-1, m=(b+e)/2; + Int_t nc=fN[s]; + if (nc==0) return kMaxClusterPerSector*s; + + Int_t b=kMaxClusterPerSector*s; + if (z <= fClusters[b]->GetZ()) return b; + + Int_t e=b+nc-1; + if (z > fClusters[e]->GetZ()) return e+1; + + Int_t m=(b+e)/2; for (; b fClusters[m]->GetZ()) b=m+1; else e=m; @@ -798,48 +776,80 @@ Int_t AliITStrackerV2::AliITSlayer::FindClusterIndex(Double_t z) const { return m; } -void AliITStrackerV2::AliITSlayer:: -SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) { +Int_t AliITStrackerV2::AliITSlayer:: +SelectClusters(Float_t zmin,Float_t zmax,Float_t ymin, Float_t ymax) { //-------------------------------------------------------------------- - // This function sets the "window" + // This function selects clusters within the "window" //-------------------------------------------------------------------- - fI=FindClusterIndex(zmin); fZmax=zmax; - Double_t circle=2*TMath::Pi()*fR; - if (ymax>circle) { ymax-=circle; ymin-=circle; } - fYmin=ymin; fYmax=ymax; + Float_t circ=fR*TMath::TwoPi(); + + if (ymin>circ) ymin-=circ; else if (ymin<0) ymin+=circ; + if (ymax>circ) ymax-=circ; else if (ymax<0) ymax+=circ; + + Int_t i1=Int_t(kNsector*ymin/circ); if (i1==kNsector) i1--; + if (fN[i1]!=0) { + Float_t ym = (ymaxIsUsed()) continue; + if (c->GetZ()>zmax) break; + if (c->GetPhiR()<=ymin) continue; + if (c->GetPhiR()>ym) continue; + fIndex[fNsel++]=i; + } + } + + Int_t i2=Int_t(kNsector*ymax/circ); if (i2==kNsector) i2--; + if (i2==i1) return fNsel; + + if (fN[i2]!=0) { + Float_t ym = (ymin>ymax) ? ymin-circ : ymin; + Int_t i=FindClusterIndex(zmin,i2), imax=i2*kMaxClusterPerSector+fN[i2]; + for (; iIsUsed()) continue; + if (c->GetZ()>zmax) break; + if (c->GetPhiR()<=ym) continue; + if (c->GetPhiR()>ymax) continue; + fIndex[fNsel++]=i; + } + } + + return fNsel; } -const AliITSclusterV2 *AliITStrackerV2::AliITSlayer::GetNextCluster(Int_t &ci){ +const AliITSRecPoint *AliITStrackerV2::AliITSlayer::GetNextCluster(Int_t &ci){ //-------------------------------------------------------------------- // This function returns clusters within the "window" //-------------------------------------------------------------------- - const AliITSclusterV2 *cluster=0; - for (Int_t i=fI; iGetZ() > fZmax) break; - // if (c->IsUsed()) continue; - const AliITSdetector &det=GetDetector(c->GetDetectorIndex()); - Double_t y=fR*det.GetPhi() + c->GetY(); - - if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi(); - if (y>1.*fR*TMath::Pi() && fYmaxfYmax) continue; - cluster=c; ci=i; - fI=i+1; - break; + AliITSRecPoint *c=0; + ci=-1; + if (fNsel) { + fNsel--; + ci=fIndex[fNsel]; + c=fClusters[ci]; } + return c; +} - return cluster; +Int_t AliITStrackerV2::AliITSlayer::GetNumberOfClusters() const { + Int_t n=0; + for (Int_t s=0; s= 2*TMath::Pi()) dphi -= 2*TMath::Pi(); Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5); @@ -983,37 +993,17 @@ Double_t AliITStrackerV2::GetEffectiveThickness(Double_t y,Double_t z) const return d/(xn*xn); } -Int_t AliITStrackerV2::AliITSlayer::InRoad() const { - //-------------------------------------------------------------------- - // This function returns number of clusters within the "window" - //-------------------------------------------------------------------- - Int_t ncl=0; - for (Int_t i=fI; iGetZ() > fZmax) break; - if (c->IsUsed()) continue; - const AliITSdetector &det=GetDetector(c->GetDetectorIndex()); - Double_t y=fR*det.GetPhi() + c->GetY(); - - if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi(); - if (y>1.*fR*TMath::Pi() && fYmaxfYmax) continue; - ncl++; - } - return ncl; -} - -Bool_t -AliITStrackerV2::RefitAt(Double_t xx,AliITStrackV2 *t,const AliITStrackV2 *c) { +Bool_t AliITStrackerV2::RefitAt(Double_t xx,AliITStrackV2 *t, + const AliITStrackV2 *c, Bool_t extra) { //-------------------------------------------------------------------- // This function refits the track "t" at the position "x" using // the clusters from "c" + // If "extra"==kTRUE, + // the clusters from overlapped modules get attached to "t" //-------------------------------------------------------------------- - Int_t index[kMaxLayer]; + Int_t index[AliITSgeomTGeo::kNLayers]; Int_t k; - for (k=0; kGetNumberOfClusters(); for (k=0; kGetClusterIndex(k),nl=(idx&0xf0000000)>>28; @@ -1022,10 +1012,10 @@ AliITStrackerV2::RefitAt(Double_t xx,AliITStrackV2 *t,const AliITStrackV2 *c) { Int_t from, to, step; if (xx > t->GetX()) { - from=0; to=kMaxLayer; + from=0; to=AliITSgeomTGeo::GetNLayers(); step=+1; } else { - from=kMaxLayer-1; to=-1; + from=AliITSgeomTGeo::GetNLayers()-1; to=-1; step=-1; } @@ -1052,11 +1042,11 @@ AliITStrackerV2::RefitAt(Double_t xx,AliITStrackV2 *t,const AliITStrackV2 *c) { } // - Double_t x,y,z; - if (!t->GetGlobalXYZat(r,x,y,z)) { + Double_t phi,z; + if (!t->GetPhiZat(r,phi,z)) { return kFALSE; } - Double_t phi=TMath::ATan2(y,x); + Int_t idet=layer.FindDetectorIndex(phi,z); if (idet<0) { return kFALSE; @@ -1068,52 +1058,37 @@ AliITStrackerV2::RefitAt(Double_t xx,AliITStrackV2 *t,const AliITStrackV2 *c) { } t->SetDetectorIndex(idet); - const AliITSclusterV2 *cl=0; - Double_t maxchi2=kMaxChi2; + const AliITSRecPoint *cl=0; + Double_t maxchi2=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); Int_t idx=index[i]; - if (idx>0) { - const AliITSclusterV2 *c=(AliITSclusterV2 *)GetCluster(idx); - if (idet != c->GetDetectorIndex()) { - idet=c->GetDetectorIndex(); - const AliITSdetector &det=layer.GetDetector(idet); - if (!t->Propagate(det.GetPhi(),det.GetR())) { + if (idx>=0) { + const AliITSRecPoint *ccc=(AliITSRecPoint *)GetCluster(idx); + if (idet != ccc->GetDetectorIndex()) { + idet=ccc->GetDetectorIndex(); + const AliITSdetector &det2=layer.GetDetector(idet); + if (!t->Propagate(det2.GetPhi(),det2.GetR())) { return kFALSE; } t->SetDetectorIndex(idet); } - Double_t chi2=t->GetPredictedChi2(c); + Double_t chi2=t->GetPredictedChi2(ccc); if (chi2GetNumberOfClusters()>2) { - Double_t dz=4*TMath::Sqrt(t->GetSigmaZ2()+kSigmaZ2[i]); - Double_t dy=4*TMath::Sqrt(t->GetSigmaY2()+kSigmaY2[i]); - Double_t zmin=t->GetZ() - dz; - Double_t zmax=t->GetZ() + dz; - Double_t ymin=t->GetY() + phi*r - dy; - Double_t ymax=t->GetY() + phi*r + dy; - layer.SelectClusters(zmin,zmax,ymin,ymax); - - const AliITSclusterV2 *c=0; Int_t ci=-1; - while ((c=layer.GetNextCluster(ci))!=0) { - if (idet != c->GetDetectorIndex()) continue; - Double_t chi2=t->GetPredictedChi2(c); - if (chi2GetX()+cl->GetX(); + if (!t->PropagateTo(x,0.,0.)) return kFALSE; if (!t->Update(cl,maxchi2,idx)) { return kFALSE; } - t->SetSampledEdx(cl->GetQ(),t->GetNumberOfClusters()-1); + t->SetSampledEdx(cl->GetQ(),i-2); } { @@ -1122,6 +1097,37 @@ AliITStrackerV2::RefitAt(Double_t xx,AliITStrackV2 *t,const AliITStrackV2 *c) { t->CorrectForMaterial(-step*d,x0); } + if (extra) { //search for extra clusters + AliITStrackV2 tmp(*t); + Double_t dz=4*TMath::Sqrt(tmp.GetSigmaZ2()+AliITSReconstructor::GetRecoParam()->GetSigmaZ2(i)); + if (dz < 0.5*TMath::Abs(tmp.GetTgl())) dz=0.5*TMath::Abs(tmp.GetTgl()); + Double_t dy=4*TMath::Sqrt(t->GetSigmaY2()+AliITSReconstructor::GetRecoParam()->GetSigmaY2(i)); + if (dy < 0.5*TMath::Abs(tmp.GetSnp())) dy=0.5*TMath::Abs(tmp.GetSnp()); + Double_t zmin=t->GetZ() - dz; + Double_t zmax=t->GetZ() + dz; + Double_t ymin=t->GetY() + phi*r - dy; + Double_t ymax=t->GetY() + phi*r + dy; + layer.SelectClusters(zmin,zmax,ymin,ymax); + + const AliITSRecPoint *cx=0; Int_t ci=-1,cci=-1; + maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2(); + Double_t tolerance=0.1; + while ((cx=layer.GetNextCluster(ci))!=0) { + if (idet == cx->GetDetectorIndex()) continue; + + const AliITSdetector &detx=layer.GetDetector(cx->GetDetectorIndex()); + + if (!tmp.Propagate(detx.GetPhi(),detx.GetR())) continue; + + if (TMath::Abs(tmp.GetZ() - cx->GetZ()) > tolerance) continue; + if (TMath::Abs(tmp.GetY() - cx->GetY()) > tolerance) continue; + + Double_t chi2=tmp.GetPredictedChi2(cx); + if (chi2=0) t->SetExtraCluster(i,(i<<28)+cci); + } + // track time update [SR, GSI 17.02.2003] if (t->IsStartedTimeIntegral() && step==1) { Double_t newX, newY, newZ; @@ -1138,566 +1144,23 @@ AliITStrackerV2::RefitAt(Double_t xx,AliITStrackV2 *t,const AliITStrackV2 *c) { return kTRUE; } - -Float_t *AliITStrackerV2::GetWeight(Int_t index) { - //-------------------------------------------------------------------- - // Return pointer to a given cluster - //-------------------------------------------------------------------- - Int_t l=(index & 0xf0000000) >> 28; - Int_t c=(index & 0x0fffffff) >> 00; - return fgLayers[l].GetWeight(c); -} - - - void AliITStrackerV2::UseClusters(const AliKalmanTrack *t, Int_t from) const { //-------------------------------------------------------------------- // This function marks clusters assigned to the track //-------------------------------------------------------------------- AliTracker::UseClusters(t,from); - AliITSclusterV2 *c=(AliITSclusterV2 *)GetCluster(t->GetClusterIndex(0)); - //if (c->GetQ()>2) c->Use(); - if (c->GetSigmaZ2()>0.1) c->Use(); - c=(AliITSclusterV2 *)GetCluster(t->GetClusterIndex(1)); - //if (c->GetQ()>2) c->Use(); - if (c->GetSigmaZ2()>0.1) c->Use(); - -} - - -void AliITStrackerV2::AddTrackHypothesys(AliITStrackV2 * track, Int_t esdindex) -{ - //------------------------------------------------------------------ - // add track to the list of hypothesys - //------------------------------------------------------------------ - - if (esdindex>=fTrackHypothesys.GetEntriesFast()) fTrackHypothesys.Expand(esdindex*2+10); - // - TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex); - if (!array) { - array = new TObjArray(10); - fTrackHypothesys.AddAt(array,esdindex); - } - array->AddLast(track); -} - -void AliITStrackerV2::SortTrackHypothesys(Int_t esdindex, Float_t likelihoodlevel) -{ - //------------------------------------------------------------------- - // compress array of track hypothesys - // keep only maxsize best hypothesys - //------------------------------------------------------------------- - if (esdindex>fTrackHypothesys.GetEntriesFast()) return; - if (! (fTrackHypothesys.At(esdindex)) ) return; - TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex); - Int_t entries = array->GetEntriesFast(); - - Float_t * chi2 = new Float_t[entries]; - Float_t * probability = new Float_t[entries]; - Float_t sumprobabilityall=0; - Int_t * index = new Int_t[entries]; - // - // - for (Int_t itrack=0;itrackGetEntriesFast();itrack++){ - AliITStrackV2 * track = (AliITStrackV2*)array->At(itrack); - // - if (track && track->GetNumberOfClusters()>(track->GetNUsed()+track->GetNSkipped())){ - // - chi2[itrack] = track->GetChi2()/(track->GetNumberOfClusters()-track->GetNUsed()-track->GetNSkipped()); - if (track->GetESDtrack()) - if (track->GetESDtrack()->GetTPCsignal()>80){ - track->CookdEdx(); - if ((track->GetdEdx()/track->GetESDtrack()->GetTPCsignal())<0.4){ - Float_t factor = 2.+10.*(0.6-track->GetdEdx()/track->GetESDtrack()->GetTPCsignal()); - chi2[itrack]*= factor; //mismatch in dEdx - } - } - } - else - chi2[itrack] = 10000000.; - probability[itrack] = 1./(0.3+chi2[itrack]); - sumprobabilityall+=probability[itrack]; - } - - TMath::Sort(entries,chi2,index,kFALSE); - TObjArray * newarray = new TObjArray(); - Float_t sumprobability = 0.; - for (Int_t i=0;iAt(index[i]); - if (!track) break; - if (chi2[index[i]]<30){ - newarray->AddLast(array->RemoveAt(index[i])); - track->SetChi2MIP(0,chi2[index[i]]); // base chi 2 - sumprobability+= probability[index[i]]/sumprobabilityall; - if (sumprobability> likelihoodlevel) break; - } - else{ - delete array->RemoveAt(index[i]); - } - } + Int_t clusterIndex = t->GetClusterIndex(0); + AliITSRecPoint *c= 0x0; - array->Delete(); - delete fTrackHypothesys.RemoveAt(esdindex); - fTrackHypothesys.AddAt(newarray,esdindex); + if (clusterIndex>-1) + c = (AliITSRecPoint *)GetCluster(clusterIndex); + if (c && c->GetSigmaZ2()>0.1) c->UnUse(); - delete [] chi2; - delete [] probability; - delete [] index; + c = 0x0; + clusterIndex = t->GetClusterIndex(1); + if (clusterIndex>-1) + c=(AliITSRecPoint *)GetCluster(clusterIndex); + if (c && c->GetSigmaZ2()>0.1) c->UnUse(); } - - -void AliITStrackerV2::CompressTrackHypothesys(Int_t esdindex, Float_t likelihoodlevel, Int_t maxsize) -{ - // - // - if (esdindex>fTrackHypothesys.GetEntriesFast()) return; - if (! (fTrackHypothesys.At(esdindex)) ) return; - TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex); - Int_t entries = array->GetEntriesFast(); - // - if (likelihoodlevel>0.000001){ - Float_t *probability = new Float_t[entries]; - for (Int_t i=0;iAt(itrack); - probability[itrack]=0; - if (!track) continue; - probability[itrack] = 1./(0.3+track->GetChi2MIP(0)); - sumprobabilityall += probability[itrack]; - // - } - if (sumprobabilityall<=0.000000000000001){ - return; - } - for (Int_t itrack=0;itrackAt(index[i]); - if (!track) continue; - newarray->AddLast(array->RemoveAt(index[i])); - sumprobability+= probability[index[i]]; - if (sumprobability> likelihoodlevel) break; - if (i>maxsize) break; - } - - array->Delete(); - delete fTrackHypothesys.RemoveAt(esdindex); - fTrackHypothesys.AddAt(newarray,esdindex); - // - delete []index; - delete []probability; - } - else{ - array->Delete(); - delete fTrackHypothesys.RemoveAt(esdindex); - } -} - - -AliITStrackV2 * AliITStrackerV2::GetBestHypothesys(Int_t esdindex, AliITStrackV2 * original, Int_t checkmax) -{ - //------------------------------------------------------------- - // try to find best hypothesy - // currently - minimal chi2 of track+backpropagated track+matching to the tpc track - //------------------------------------------------------------- - if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0; - TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex); - if (!array) return 0; - Int_t entries = array->GetEntriesFast(); - if (!entries) return 0; - Float_t minchi2 = 100000; - Int_t maxn = 3; - AliITStrackV2 * besttrack=0; - Int_t accepted =0; - Int_t maxindex=0; - // - Float_t sumz2=0; - Float_t sumy2=0; - Float_t sumall=0; - for (Int_t itrack=0;itrackAt(itrack); - if (!track) continue; - sumall++; - sumz2+=track->GetSigmaZ2(); - sumy2+=track->GetSigmaY2(); - } - sumz2/=sumall; - sumy2/=sumall; - - Float_t dedxmismatch=1; - for (Int_t i=0;iAt(i); - if (!track) continue; - track->SetChi2MIP(1,1000000); - track->SetChi2MIP(2,1000000); - - if ( (track->GetNumberOfClusters()-track->GetNSkipped()-track->GetNUsed())<2) continue; - // - if (track->GetESDtrack()) - if (track->GetESDtrack()->GetTPCsignal()>80){ - track->CookdEdx(); - if ((track->GetdEdx()/track->GetESDtrack()->GetTPCsignal())<0.4){ - //mismatch in dEdx - dedxmismatch= 2.+ 10.*(0.6-track->GetdEdx()/track->GetESDtrack()->GetTPCsignal()); - } - } - - // track->SetLabel(original->GetLabel()); - //CookLabel(track,0.0); - //if (track->GetFakeRatio()>0.01) continue; - // - // - // backtrack - AliITStrackV2 * backtrack = new AliITStrackV2(*track); - backtrack->ResetCovariance(); - backtrack->ResetClusters(); - Double_t x = original->GetX(); - if (!RefitAt(x,backtrack,track)){ - delete backtrack; - delete array->RemoveAt(i); - continue; - } - if ( (backtrack->GetChi2() / float(backtrack->GetNumberOfClusters()-track->GetNSkipped()-track->GetNUsed()-0.5))>6) - { - delete backtrack; - delete array->RemoveAt(i); - continue; - } - Double_t deltac = backtrack->GetC()-original->GetC(); - Double_t deltatgl = backtrack->GetTgl()-original->GetTgl(); - // - Double_t poolc2 = (deltac*deltac)/(original->GetCov44()+backtrack->GetCov44()); - Double_t pooltgl2 = (deltatgl*deltatgl)/(original->GetCov33()+backtrack->GetCov33()); - if ((poolc2+pooltgl2)>32){ //4 sigma - delete backtrack; - delete array->RemoveAt(i); - continue; - } - //Double_t bpoolc = (deltac*deltac)/(original->GetCov44()); - //Double_t bpooltgl = (deltatgl*deltatgl)/(original->GetCov33()); - - // - //forward track - without constraint - AliITStrackV2 * forwardtrack = new AliITStrackV2(*original); - // forwardtrack->ResetCovariance(); - forwardtrack->ResetClusters(); - x = track->GetX(); - if (!RefitAt(x,forwardtrack,track)){ - delete forwardtrack; - delete backtrack; - delete array->RemoveAt(i); - continue; - } - if ( (forwardtrack->GetChi2()/float(forwardtrack->GetNumberOfClusters()-track->GetNSkipped()-track->GetNUsed()))>6) { - delete forwardtrack; - delete backtrack; - delete array->RemoveAt(i); - continue; - } - // - accepted++; - if (accepted>checkmax){ - delete backtrack; - delete forwardtrack; - break; - } - Double_t chi2 = (backtrack->GetChi2()/(backtrack->GetNumberOfClusters()-1-track->GetNSkipped()-track->GetNUsed())+ - forwardtrack->GetChi2()/(forwardtrack->GetNumberOfClusters()-track->GetNSkipped()-track->GetNUsed())); - // bpoolc+bpooltgl; - // chi2 *= (forwardtrack->GetSigmaZ2()/sumz2+forwardtrack->GetSigmaY2()/sumy2); - chi2 *= dedxmismatch; - // - // - track->SetChi2MIP(1,backtrack->GetChi2()/(backtrack->GetNumberOfClusters()-1-track->GetNSkipped()-track->GetNUsed())); - track->SetChi2MIP(2,forwardtrack->GetChi2()/(forwardtrack->GetNumberOfClusters()-track->GetNSkipped()-track->GetNUsed())); - track->SetChi2MIP(3,poolc2+pooltgl2); - // - - if (track->GetNumberOfClusters()>maxn){ - delete besttrack; - besttrack = new AliITStrackV2(*forwardtrack); - maxn = track->GetNumberOfClusters(); - minchi2 = chi2; - delete backtrack; - delete forwardtrack; - continue; - } - // - if (chi2 < minchi2){ - besttrack = new AliITStrackV2(*forwardtrack); - minchi2 = chi2; - } - delete backtrack; - delete forwardtrack; - } - // - // - if (!besttrack || besttrack->GetNumberOfClusters()<4) { - delete besttrack; - return 0; - } - - // - besttrack->SetLabel(original->GetLabel()); - CookLabel(besttrack,0.0); - // - // calculate "weight of the cluster" - // - - { - //sign usage information for clusters - Int_t clusterindex[6][100]; - Double_t clusterweight[6][100]; - for (Int_t ilayer=0;ilayer<6;ilayer++) - for (Int_t icluster=0;icluster<100;icluster++){ - clusterindex[ilayer][icluster] = -1; - clusterweight[ilayer][icluster] = 0; - } - //printf("%d\t%d\n",esdindex, entries); - // - Float_t sumchi2=0; - for (Int_t itrack=0;itrackAt(itrack); - if (!track) continue; - if (track->GetChi2MIP(1)>1000) continue; //not accepted - sumchi2 +=1./(0.3+track->GetChi2MIP(1)+track->GetChi2MIP(2)); - } - for (Int_t itrack=0;itrackAt(itrack); - if (!track) continue; - if (track->GetChi2MIP(1)>1000) continue; //not accepted - for (Int_t icluster=0;iclusterGetNumberOfClusters();icluster++){ - Int_t tindex = track->GetClusterIndex(icluster); - Int_t ilayer = (tindex & 0xf0000000) >> 28; - if (tindex<0) continue; - Int_t cindex =0; - // - for (cindex=0;cindex<100;cindex++){ - if (clusterindex[ilayer][cindex]<0) break; - if (clusterindex[ilayer][cindex]==tindex) break; - } - if (cindex>100) break; - if (clusterindex[ilayer][cindex]!=tindex) clusterindex[ilayer][cindex] = tindex; - clusterweight[ilayer][cindex]+= (1./(0.3+track->GetChi2MIP(1)+track->GetChi2MIP(2)))* (1./sumchi2); - Float_t *weight = GetWeight(tindex); - - if (weight){ - *weight+= (1./(0.3+track->GetChi2MIP(1)+track->GetChi2MIP(2)))* (1./sumchi2); - } - } - } - - if (besttrack->GetChi2()/besttrack->GetNumberOfClusters()>3.5) return besttrack; //don't sign clusters - Int_t current=0; - Double_t deltad = besttrack->GetD(GetX(),GetY()); - Double_t deltaz = besttrack->GetZat(GetX()) - GetZ(); - Double_t deltaprim = TMath::Sqrt(deltad*deltad+deltaz*deltaz); - - for (Int_t icluster=0;iclusterGetNumberOfClusters();icluster++){ - Int_t index = besttrack->GetClusterIndex(icluster); - Int_t ilayer = (index & 0xf0000000) >> 28; - AliITSclusterV2 *c = (AliITSclusterV2*)GetCluster(index); - if (!c) continue; - // - for (Int_t icluster=0;icluster<100;icluster++){ - // sign non "doubt" clusters as used - if (clusterindex[ilayer][icluster]!=index) continue; - // Float_t * weight = GetWeight(index); - //if (weight) if (*weight>1){ - // if (c->IsUsed()) continue; - // c->Use(); - //} - if ( (ilayer*0.2+0.2)GetNy()>4) continue; // don sign cluster - if ( (ilayer>1&&clusterweight[ilayer][icluster]>0.7) || (ilayer<2&&clusterweight[ilayer][icluster]>0.8) ){ - current++; - if (c->IsUsed()) continue; - c->Use(); - } - } - } - } - - // - return besttrack; -} - - -AliITStrackV2 * AliITStrackerV2::GetBestHypothesysMIP(Int_t esdindex, AliITStrackV2 * original) -{ - //------------------------------------------------------------- - // try to find best hypothesy - // currently - minimal chi2 of track+backpropagated track+matching to the tpc track - //------------------------------------------------------------- - if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0; - TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex); - if (!array) return 0; - Int_t entries = array->GetEntriesFast(); - if (!entries) return 0; - AliITStrackV2 * besttrack=0; - // - //sign usage information for clusters - Int_t clusterindex[6][100]; - Double_t clusterweight[6][100]; - for (Int_t ilayer=0;ilayer<6;ilayer++) - for (Int_t icluster=0;icluster<100;icluster++){ - clusterindex[ilayer][icluster] = -1; - clusterweight[ilayer][icluster] = 0; - } - //printf("%d\t%d\n",esdindex, entries); - // - Float_t sumchi2=0; - for (Int_t itrack=0;itrackAt(itrack); - if (!track) continue; - if (track->GetChi2MIP(1)>1000) continue; //not accepted - before - sumchi2 +=1./(0.3+track->GetChi2MIP(1)+track->GetChi2MIP(2)); - } - // - // get cluster weight - for (Int_t itrack=0;itrackAt(itrack); - if (!track) continue; - if (track->GetChi2MIP(1)>1000) continue; // track not accepted in previous iterration - for (Int_t icluster=0;iclusterGetNumberOfClusters();icluster++){ - Int_t tindex = track->GetClusterIndex(icluster); - Int_t ilayer = (tindex & 0xf0000000) >> 28; - if (tindex<0) continue; - Int_t cindex =0; - // - for (cindex=0;cindex<100;cindex++){ - if (clusterindex[ilayer][cindex]<0) break; - if (clusterindex[ilayer][cindex]==tindex) break; - } - if (cindex>100) break; - if (clusterindex[ilayer][cindex]!=tindex) clusterindex[ilayer][cindex] = tindex; - clusterweight[ilayer][cindex]+= (1./(0.3+track->GetChi2MIP(1)+track->GetChi2MIP(2)))* (1./sumchi2); - } - } - // - // get cluster relative sharing - factor - // - // - for (Int_t ilayer=0;ilayer<6;ilayer++) - for (Int_t cindex=0;cindex<100;cindex++){ - if (clusterindex[ilayer][cindex]<0) continue; - Int_t tindex = clusterindex[ilayer][cindex]; - Float_t *weight = GetWeight(tindex); - if (!weight){ - printf("Problem 1\n"); // not existing track - continue; - } - if (*weight<(clusterweight[ilayer][cindex]-0.00001)){ - printf("Problem 2\n"); // not normalized probability - continue; - } - AliITSclusterV2 *c = (AliITSclusterV2*)GetCluster(tindex); - if (c->GetNy()<5){ - clusterweight[ilayer][cindex]/= *weight; - } - else{ - Float_t weight2 = TMath::Max(*weight-0.5*clusterweight[ilayer][cindex],0.0000001); - clusterweight[ilayer][cindex]/= weight2; - if ( clusterweight[ilayer][cindex]>1) clusterweight[ilayer][cindex] =1.; - } - } - // - //take to the account sharing factor - // - Float_t chi2 =10000000; - Float_t sharefactor=0; - Float_t minchi2 = 100000000; - Float_t secchi2 = 100000000; - Float_t norm=0; - for (Int_t itrack=0;itrackAt(itrack); - if (!track) continue; - if (track->GetChi2MIP(1)>1000) continue; //not accepted - before - chi2 = track->GetChi2MIP(1); - Float_t newchi2=0; - sharefactor=0; - norm =0; - // - for (Int_t icluster=0;iclusterGetNumberOfClusters();icluster++){ - Int_t tindex = track->GetClusterIndex(icluster); - Int_t ilayer = (tindex & 0xf0000000) >> 28; - if (tindex<0) continue; - Int_t cindex =0; - Float_t cchi2 = (track->GetDy(ilayer)*track->GetDy(ilayer))/(track->GetSigmaY(ilayer)*track->GetSigmaY(ilayer)) + - (track->GetDz(ilayer)*track->GetDz(ilayer))/(track->GetSigmaZ(ilayer)*track->GetSigmaZ(ilayer)) ; - // - for (cindex=0;cindex<100;cindex++){ - if (clusterindex[ilayer][cindex]<0) break; - if (clusterindex[ilayer][cindex]==tindex) break; - } - if (cindex>100) continue; - if (clusterweight[ilayer][cindex]>0.00001){ - sharefactor+= clusterweight[ilayer][cindex]; - cchi2/=clusterweight[ilayer][cindex]; - norm++; - } - newchi2+=cchi2; - } - newchi2/=(norm-track->GetNSkipped()-track->GetNUsed()); - track->SetChi2MIP(4,newchi2); - //if (norm>0) sharefactor/=norm; - //if (sharefactor<0.5) return 0; - // chi2*=1./(0.5+sharefactor); - if (newchi2GetNumberOfClusters()<4) { - return 0; - } - // - if ((minchi2/secchi2)<0.7){ - // - //increase the weight for clusters if the probability of other hypothesys is small - Double_t deltad = besttrack->GetD(GetX(),GetY()); - Double_t deltaz = besttrack->GetZat(GetX()) - GetZ(); - Double_t deltaprim = TMath::Sqrt(deltad*deltad+deltaz*deltaz); - // - for (Int_t icluster=0;iclusterGetNumberOfClusters();icluster++){ - Int_t index = besttrack->GetClusterIndex(icluster); - Int_t ilayer = (index & 0xf0000000) >> 28; - AliITSclusterV2 *c = (AliITSclusterV2*)GetCluster(index); - if (!c) continue; - if (c->GetNy()>3) continue; - if ( (ilayer*0.2+0.2)SetLabel(original->GetLabel()); - CookLabel(besttrack,0.0); - - // - return besttrack; -} -