X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITStrackerV2.cxx;h=a1b1d6a2c5099c77776df72b390c4d00b5272709;hb=50f91c8bdc57f4d3b826ff69d617e2672177ed2a;hp=3b447b20266df4c2b8ce6d731f50e4a9204d2211;hpb=743a19f2b0103cdccad62e316effbd84ccb23b55;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/AliITStrackerV2.cxx b/ITS/AliITStrackerV2.cxx index 3b447b20266..a1b1d6a2c50 100644 --- a/ITS/AliITStrackerV2.cxx +++ b/ITS/AliITStrackerV2.cxx @@ -12,447 +12,456 @@ * 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 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 //------------------------------------------------------------------------- + +#include + +#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 "AliTPCtrack.h" -#include "AliITSclusterV2.h" +#include "AliITSReconstructor.h" #include "AliITStrackerV2.h" -//#define DEBUG +ClassImp(AliITStrackerV2) + +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); -AliITStrackerV2::AliITSlayer AliITStrackerV2::fLayers[kMaxLayer]; // ITS layers + for (Int_t i=0; iGetNladders(i); - Int_t ndet=g->GetNdetectors(i); + fConstraint[0]=t.fConstraint[0]; fConstraint[1]=t.fConstraint[1]; - g->GetTrans(i,1,1,x,y,z); - Double_t r=TMath::Sqrt(x*x + y*y); + Double_t xyz[]={AliITSReconstructor::GetRecoParam()->GetXVdef(), + 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; iGetTrans(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; - new (fLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet); + new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet); for (Int_t j=1; jGetTrans(i,j,k,x,y,zshift); - Double_t rot[9]; g->GetRotMatrix(i,j,k,rot); - - Double_t r =-x*rot[1] + y*rot[0]; if (i==1) r=-r; - Double_t phi=TMath::ATan2(rot[1],rot[0]); if (i==1) phi-=3.1415927; - phi+=0.5*TMath::Pi(); - AliITSdetector &det=fLayers[i-1].GetDetector((j-1)*ndet + k-1); - -if (phi<0) phi += 2*TMath::Pi(); - + 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); } } } - try { - //Read clusters - //MI change - char cname[100]; - sprintf(cname,"TreeC_ITS_%d",eventn); - TTree *cTree=(TTree*)gDirectory->Get(cname); - // TTree *cTree=(TTree*)gDirectory->Get("cTree"); - - if (!cTree) throw - ("AliITStrackerV2::AliITStrackerV2 can't get cTree !\n"); - - TBranch *branch=cTree->GetBranch("Clusters"); - if (!branch) throw - ("AliITStrackerV2::AliITStrackerV2 can't get Clusters branch !\n"); - - TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy; - branch->SetAddress(&clusters); - - Int_t nentr=(Int_t)cTree->GetEntries(); - for (i=0; iGetEvent(i)) continue; - //Int_t lay,lad,det; g->GetModuleId(i-1,lay,lad,det); - Int_t lay,lad,det; g->GetModuleId(i,lay,lad,det); - Int_t ncl=clusters->GetEntriesFast(); - while (ncl--) { - AliITSclusterV2 *c=(AliITSclusterV2*)clusters->UncheckedAt(ncl); - -#ifdef DEBUG -if (c->GetLabel(2)!=LAB) -if (c->GetLabel(1)!=LAB) -if (c->GetLabel(0)!=LAB) continue; -cout<GetY()<<' '<GetZ()<GetDetectorIndex()); - Double_t r=det.GetR(); r=TMath::Sqrt(r*r + c->GetY()*c->GetY()); - if (r > TMath::Abs(c->GetZ())-0.5) - fLayers[lay-1].InsertCluster(new AliITSclusterV2(*c)); - } - clusters->Delete(); - } - } - catch (const Char_t *msg) { - cerr<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; iIsOpen()) { - cerr<<"AliITStrackerV2::Clusters2Tracks(): "; - cerr<<"file with TPC tracks is not open !\n"; - return 1; +Int_t AliITStrackerV2::LoadClusters(TTree *cTree) { + //-------------------------------------------------------------------- + //This function loads ITS clusters + //-------------------------------------------------------------------- + TBranch *branch=cTree->GetBranch("ITSRecPoints"); + if (!branch) { + Error("LoadClusters"," can't get the branch !\n"); + return 1; } - if (!out->IsOpen()) { - cerr<<"AliITStrackerV2::Clusters2Tracks(): "; - cerr<<"file for ITS tracks is not open !\n"; - return 2; - } + TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy; + branch->SetAddress(&clusters); - in->cd(); - // - char tname[100]; - sprintf(tname,"TreeT_TPC_%d",fEventN); - TTree *tpcTree=(TTree*)in->Get(tname); - //TTree *tpcTree=(TTree*)in->Get("TPCf"); - - if (!tpcTree) { - cerr<<"AliITStrackerV2::Clusters2Tracks() "; - cerr<<"can't get a tree with TPC tracks !\n"; - return 3; - } + Int_t j=0; + for (Int_t i=0; iSetBranchAddress("tracks",&itrack); + Double_t r=fgLayers[i].GetR(); + Double_t circ=TMath::TwoPi()*r; - out->cd(); - TTree itsTree("ITSf","Tree with ITS tracks"); - AliITStrackV2 *otrack=&fBestTrack; - itsTree.Branch("tracks","AliITStrackV2",&otrack,32000,0); + for (; jGetEvent(j)) continue; + Int_t ncl=clusters->GetEntriesFast(); + + while (ncl--) { + AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl); - Int_t ntrk=0; + if (!c->Misalign()) AliWarning("Can't misalign this cluster !"); - Int_t nentr=(Int_t)tpcTree->GetEntries(); - for (Int_t i=0; iGetDetectorIndex(); + 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); - if (!tpcTree->GetEvent(i)) continue; - Int_t tpcLabel=itrack->GetLabel(); //save the TPC track label + fgLayers[i].InsertCluster(new AliITSRecPoint(*c)); + } + clusters->Delete(); + } + fgLayers[i].ResetRoad(); //road defined by the cluster density + } -#ifdef DEBUG -lbl=tpcLabel; -if (TMath::Abs(tpcLabel)!=LAB) continue; -cout<Gaus(0.,kSigmaYV); fZV=gRandom->Gaus(0.,kSigmaZV); - - Double_t r2=fTrackToFollow.GetX()*fTrackToFollow.GetX(); - Double_t x0=0.082/21.82*2.33*(45*45+40*40)/r2+2.000/41.28*0.03*75*75/r2; - fTrackToFollow.Improve(x0,fYV,fZV); - - Double_t xk=80.; - fTrackToFollow.PropagateTo(xk,0.,0.); //Ne if it's still there - - xk-=0.005; - fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar - xk-=0.02; - fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45); //Kevlar - xk-=2.0; - fTrackToFollow.PropagateTo(xk, 2.0/41.28*0.029, 2.0*0.029);//Nomex - xk-=0.02; - fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45); //Kevlar - xk-=0.005; - fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar - - xk=61.; - fTrackToFollow.PropagateTo(xk,0.,0.); //C02 - - xk -=0.005; - fTrackToFollow.PropagateTo(xk, 0.005/24.01*2.7, 0.005*2.7); //Al - xk -=0.005; - fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar - xk -=0.02; - fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45); //Kevlar - xk -=0.5; - fTrackToFollow.PropagateTo(xk, 0.5/41.28*.029, 0.5*0.029); //Nomex - xk -=0.02; - fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45); //Kevlar - xk -=0.005; - fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar - xk -=0.005; - fTrackToFollow.PropagateTo(xk, 0.005/24.01*2.7, 0.005*2.7); //Al - - - for (FollowProlongation(); fIGetX() > riw) { + if (!t->PropagateTo(riw,diw,x0iw)) return 1; + if (TMath::Abs(t->GetY())>yr) t->CorrectForMaterial(dr); + if (TMath::Abs(t->GetZ())CorrectForMaterial(dm); + if (!t->PropagateTo(rcd,dcd,x0cd)) return 1; + //Double_t x,y,z; t->GetGlobalXYZat(rr,x,y,z); + //if (TMath::Abs(y)PropagateTo(rr,dr,x0r); + if (!t->PropagateTo(rs,ds)) return 1; + } else if (t->GetX() < rs) { + if (!t->PropagateTo(rs,-ds)) return 1; + //Double_t x,y,z; t->GetGlobalXYZat(rr,x,y,z); + //if (TMath::Abs(y)PropagateTo(rr,-dr,x0r); + if (!t->PropagateTo(rcd,-dcd,x0cd)) return 1; + if (!t->PropagateTo(riw+0.001,-diw,x0iw)) return 1; + } else { + ::Error("CorrectForDeadZoneMaterial","track is already in the dead zone !"); + return 1; + } + + return 0; +} + +Int_t AliITStrackerV2::Clusters2Tracks(AliESDEvent *event) { + //-------------------------------------------------------------------- + // This functions reconstructs ITS tracks + // The clusters must be already loaded ! + //-------------------------------------------------------------------- + TObjArray itsTracks(15000); + + {/* Read ESD tracks */ + Int_t nentr=event->GetNumberOfTracks(); + Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr); + while (nentr--) { + AliESDtrack *esd=event->GetTrack(nentr); + + if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue; + if (esd->GetStatus()&AliESDtrack::kTPCout) continue; + if (esd->GetStatus()&AliESDtrack::kITSin) continue; + + AliITStrackV2 *t=0; + try { + t=new AliITStrackV2(*esd); + } catch (const Char_t *msg) { + Warning("Clusters2Tracks",msg); + delete t; + continue; + } + if (TMath::Abs(t->GetD(GetX(),GetY()))>4) { + delete t; + continue; + } + + if (CorrectForDeadZoneMaterial(t)!=0) { + Warning("Clusters2Tracks", + "failed to correct for the material in the dead zone !\n"); + delete t; + continue; + } + itsTracks.AddLast(t); } + } /* End Read ESD tracks */ + + itsTracks.Sort(); + Int_t nentr=itsTracks.GetEntriesFast(); + + Int_t ntrk=0; + for (fPass=0; fPass<2; fPass++) { + Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue; + for (Int_t i=0; iGetLabel(); //save the TPC track label + + ResetTrackToFollow(*t); + ResetBestTrack(); + + for (FollowProlongation(); fI= kMaxLayer-kLayersToSkip) { - cerr << ++ntrk << " \r"; fBestTrack.SetLabel(tpcLabel); + fBestTrack.CookdEdx(); + CookLabel(&fBestTrack,0.); //For comparison only + fBestTrack.UpdateESDtrack(AliESDtrack::kITSin); UseClusters(&fBestTrack); - itsTree.Fill(); - } - + delete itsTracks.RemoveAt(i); + ntrk++; + } } - sprintf(tname,"TreeT_ITS_%d",fEventN); - itsTree.Write(tname); - savedir->cd(); - cerr<<"Number of TPC tracks: "<GetNumberOfTracks(); + Info("PropagateBack", "Number of ESD tracks: %d\n", nentr); - if (!in->IsOpen()) { - cerr<<"AliITStrackerV2::PropagateBack(): "; - cerr<<"file with ITS tracks is not open !\n"; - return 1; - } + Int_t ntrk=0; + for (Int_t i=0; iGetTrack(i); + + if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue; + if (esd->GetStatus()&AliESDtrack::kITSout) continue; + + AliITStrackV2 *t=0; + try { + t=new AliITStrackV2(*esd); + } catch (const Char_t *msg) { + Warning("PropagateBack",msg); + delete t; + continue; + } - if (!out->IsOpen()) { - cerr<<"AliITStrackerV2::PropagateBack(): "; - cerr<<"file for back propagated ITS tracks is not open !\n"; - return 2; - } + ResetTrackToFollow(*t); - in->cd(); - TTree *itsTree=(TTree*)in->Get("ITSf"); - if (!itsTree) { - cerr<<"AliITStrackerV2::PropagateBack() "; - cerr<<"can't get a tree with ITS tracks !\n"; - return 3; + // 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(event->GetVertex())) { + fTrackToFollow.StartTimeIntegral(); + } + fTrackToFollow.PropagateTo(3.,-0.0028,65.19); + } + + fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters(); + if (RefitAt(49.,&fTrackToFollow,t)) { + if (CorrectForDeadZoneMaterial(&fTrackToFollow)!=0) { + Warning("PropagateBack", + "failed to correct for the material in the dead zone !\n"); + delete t; + continue; + } + fTrackToFollow.SetLabel(t->GetLabel()); + //fTrackToFollow.CookdEdx(); + CookLabel(&fTrackToFollow,0.); //For comparison only + fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout); + UseClusters(&fTrackToFollow); + ntrk++; + } + delete t; } - AliITStrackV2 *itrack=new AliITStrackV2; - itsTree->SetBranchAddress("tracks",&itrack); - out->cd(); - TTree backTree("ITSb","Tree with back propagated ITS tracks"); - AliTPCtrack *otrack=0; - backTree.Branch("tracks","AliTPCtrack",&otrack,32000,0); + Info("PropagateBack","Number of back propagated ITS tracks: %d\n",ntrk); - Int_t ntrk=0; + return 0; +} - Int_t nentr=(Int_t)itsTree->GetEntries(); +Int_t AliITStrackerV2::RefitInward(AliESDEvent *event) { + //-------------------------------------------------------------------- + // This functions refits ITS tracks using the + // "inward propagated" TPC tracks + // The clusters must be loaded ! + //-------------------------------------------------------------------- + Int_t nentr=event->GetNumberOfTracks(); + Info("RefitInward", "Number of ESD tracks: %d\n", nentr); + + Int_t ntrk=0; for (Int_t i=0; iGetEvent(i); - ResetTrackToFollow(*itrack); - fTrackToFollow.ResetCovariance(); fTrackToFollow.ResetClusters(); - Int_t itsLabel=fTrackToFollow.GetLabel(); //save the ITS track label + AliESDtrack *esd=event->GetTrack(i); -#ifdef DEBUG -if (TMath::Abs(itsLabel)!=LAB) continue; -cout<GetStatus()&AliESDtrack::kITSout) == 0) continue; + if (esd->GetStatus()&AliESDtrack::kITSrefit) continue; + if (esd->GetStatus()&AliESDtrack::kTPCout) + if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue; + AliITStrackV2 *t=0; try { - Int_t nc=itrack->GetNumberOfClusters(); -#ifdef DEBUG -for (Int_t k=0; kGetClusterIndex(k); - AliITSclusterV2 *c=(AliITSclusterV2*)GetCluster(index); - cout<GetLabel(0)<<' '<GetLabel(1)<<' '<GetLabel(2)<GetClusterIndex(nc); l=(idx&0xf0000000)>>28; - c=(AliITSclusterV2*)GetCluster(idx); - } - for (fI=0; fIGetDetectorIndex()) { - idet=c->GetDetectorIndex(); - const AliITSdetector &det=layer.GetDetector(idet); - r=det.GetR(); phi=det.GetPhi(); - if (!fTrackToFollow.Propagate(phi,r,0.,0.)) throw ""; - } - Double_t chi2=fTrackToFollow.GetPredictedChi2(c); - if (chi2GetClusterIndex(nc); l=(idx&0xf0000000)>>28; - c=(AliITSclusterV2*)GetCluster(idx); - } - } - - if (fTrackToFollow.GetNumberOfClusters()>2) { - Double_t dz=3*TMath::Sqrt(fTrackToFollow.GetSigmaZ2()+kSigmaZ2[fI]); - Double_t dy=3*TMath::Sqrt(fTrackToFollow.GetSigmaY2()+kSigmaY2[fI]); - Double_t zmin=fTrackToFollow.GetZ() - dz; - Double_t zmax=fTrackToFollow.GetZ() + dz; - Double_t ymin=fTrackToFollow.GetY() + phi*r - dy; - Double_t ymax=fTrackToFollow.GetY() + phi*r + dy; - layer.SelectClusters(zmin,zmax,ymin,ymax); - - const AliITSclusterV2 *cc=0; Int_t ci; - while ((cc=layer.GetNextCluster(ci))!=0) { - if (idet != cc->GetDetectorIndex()) continue; - Double_t chi2=fTrackToFollow.GetPredictedChi2(cc); - if (chi2GetLabel()); + fTrackToFollow.CookdEdx(); + CookLabel(&fTrackToFollow,0.); //For comparison only + + if (fTrackToFollow.PropagateTo(3.,0.0028,65.19)) {//The beam pipe + fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit); + AliESDtrack *esdTrack =fTrackToFollow.GetESDtrack(); + Double_t r[3]={0.,0.,0.}; + Double_t maxD=3.; + esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD); + ntrk++; + } } + delete t; } - backTree.Write(); - savedir->cd(); - cerr<<"Number of ITS tracks: "<> 28; Int_t c=(index & 0x0fffffff) >> 00; - return fLayers[l].GetCluster(c); + return fgLayers[l].GetCluster(c); } @@ -471,71 +480,83 @@ void AliITStrackerV2::FollowProlongation() { //-------------------------------------------------------------------- //This function finds a track prolongation //-------------------------------------------------------------------- - Int_t tryAgain=kLayersToSkip; + while (fI>fLastLayerToTrackTo) { + Int_t i=fI-1; - while (fI) { - fI--; -#ifdef DEBUG -cout<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() + AliITSReconstructor::GetRecoParam()->GetSigmaZ2(i)); if (dz < 0.5*TMath::Abs(track.GetTgl())) dz=0.5*TMath::Abs(track.GetTgl()); - if (dz > kMaxRoad/4) { - //cerr<<"AliITStrackerV2::FollowProlongation: too broad road in Z !\n"; - break; + if (dz > AliITSReconstructor::GetRecoParam()->GetMaxRoad()) { + //Warning("FollowProlongation","too broad road in Z !\n"); + return; } - Double_t dy=4*TMath::Sqrt(track.GetSigmaY2() + kSigmaY2[fI]); + + if (TMath::Abs(fTrackToFollow.GetZ()-GetZ()) > r+dz) return; + + //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/4) { - //cerr<<"AliITStrackerV2::FollowProlongation: too broad road in Y !\n"; - break; + if (dy > AliITSReconstructor::GetRecoParam()->GetMaxRoad()) { + //Warning("FollowProlongation","too broad road in Y !\n"); + return; } - Double_t zmin=track.GetZ() - dz; + 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); + if (layer.SelectClusters(zmin,zmax,ymin,ymax)==0) + if (fLayersNotToSkip[fI]) return; - //take another prolongation - if (!TakeNextProlongation()) if (!tryAgain--) break; - tryAgain=kLayersToSkip; + if (!TakeNextProlongation()) + if (fLayersNotToSkip[fI]) return; } @@ -545,112 +566,130 @@ cout<= nclb) { Double_t chi2=fTrackToFollow.GetChi2(); - if (chi2/ncl < kChi2PerCluster) { + if (chi2/ncl < AliITSReconstructor::GetRecoParam()->GetChi2PerCluster()) { if (ncl > nclb || chi2 < fBestTrack.GetChi2()) { ResetBestTrack(); } } } - if (fI) fI++; } - Int_t AliITStrackerV2::TakeNextProlongation() { //-------------------------------------------------------------------- - //This function takes another track prolongation + // This function takes another track prolongation + // + // dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch //-------------------------------------------------------------------- - //Double_t m[20]; - Double_t d=GetEffectiveThickness(0,0); //Think of this !!!! - - AliITSlayer &layer=fLayers[fI]; - AliITStrackV2 &t=fTracks[fI]; - - Double_t dz=4*TMath::Sqrt(t.GetSigmaZ2() + kSigmaZ2[fI]); - Double_t dy=4*TMath::Sqrt(t.GetSigmaY2() + kSigmaY2[fI]); + AliITSlayer &layer=fgLayers[fI]; + ResetTrackToFollow(fTracks[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) { - -#ifdef DEBUG -//if (fI==0) -//if (c->GetLabel(0)!=TMath::Abs(lbl)) continue; -#endif - Int_t idet=c->GetDetectorIndex(); - if (t.GetDetectorIndex()!=idet) { + if (fTrackToFollow.GetDetectorIndex()!=idet) { const AliITSdetector &det=layer.GetDetector(idet); - if (!t.Propagate(det.GetPhi(),det.GetR(),0.,0.)) { - cerr<<"AliITStrackerV2::TakeNextProlongation: propagation failed !\n"; + ResetTrackToFollow(fTracks[fI]); + if (!fTrackToFollow.Propagate(det.GetPhi(),det.GetR())) { + //Warning("TakeNextProlongation","propagation failed !\n"); continue; } - t.SetDetectorIndex(idet); - -#ifdef DEBUG -cout<layer.GetR()+dz) continue; } - if (TMath::Abs(t.GetZ() - c->GetZ()) > dz) continue; - if (TMath::Abs(t.GetY() - c->GetY()) > dy) continue; - - //m[0]=fYV; m[1]=fZV; - //chi2=t.GetPredictedChi2(c,m,d/21.82*2.33); - chi2=t.GetPredictedChi2(c); + if (TMath::Abs(fTrackToFollow.GetZ() - c->GetZ()) > dz) continue; + if (TMath::Abs(fTrackToFollow.GetY() - c->GetY()) > dy) continue; - if (chi2 chi2) continue; + chi2=ch2; + cc=c; cci=ci; + break; } -#ifdef DEBUG -cout<1) + if (TMath::Abs(fTrackToFollow.GetD(GetX(),GetY()))>4) return 0; - //if (!fTrackToFollow.Update(m,chi2,(fI<<28)+ci)) { - if (!fTrackToFollow.Update(c,chi2,(fI<<28)+ci)) { - cerr<<"AliITStrackerV2::TakeNextProlongation: filtering failed !\n"; - return 0; + fTrackToFollow. + SetSampledEdx(cc->GetQ(),fI-2); //b.b. + + { + Double_t x0; + Double_t d=layer.GetThickness(fTrackToFollow.GetY(),fTrackToFollow.GetZ(),x0); + fTrackToFollow.CorrectForMaterial(d,x0); } - fTrackToFollow.Improve(d/21.82*2.33,fYV,fZV); -#ifdef DEBUG -cout<<"accepted lab="<GetLabel(0)<<' ' - <GetZ()); - memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliITSclusterV2*)); - fClusters[i]=c; fN++; +void AliITStrackerV2::AliITSlayer::ResetRoad() { + //-------------------------------------------------------------------- + // This function calculates the road defined by the cluster density + //-------------------------------------------------------------------- + Int_t n=0; + for (Int_t s=0; sGetZ())1) fRoad=2*fR*TMath::Sqrt(3.14/n); +} +Int_t AliITStrackerV2::AliITSlayer::InsertCluster(AliITSRecPoint *c) { + //-------------------------------------------------------------------- + // This function inserts a cluster to this layer in increasing + // order of the cluster's fZ + //-------------------------------------------------------------------- + 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; + } + 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; @@ -693,126 +776,391 @@ 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 cluster; + return c; } -Int_t AliITStrackerV2::AliITSlayer:: -FindDetectorIndex(Double_t phi, Double_t z) const { +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); - - Double_t dz=fZOffset-z; - Int_t nz=Int_t(dz*(fNdetectors-1)*0.5/fZOffset+0.5); - if (np>=fNladders) np-=fNladders; if (np<0) np+=fNladders; -#ifdef DEBUG -cout<=fNdetectors) return -1; + if (nz<0) return -1; return np*fNdetectors + nz; } Double_t -AliITStrackerV2::AliITSlayer::GetThickness(Double_t phi, Double_t z) const { +AliITStrackerV2::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0) +const { //-------------------------------------------------------------------- - //This function returns the thickness of this layer + //This function returns the layer thickness at this point (units X0) //-------------------------------------------------------------------- - //-pi3.40) d+=dd; + if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);} + if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);} + for (Int_t i=0; i<12; i++) { + if (TMath::Abs(z-3.9*(i+0.5))<0.15) { + if (TMath::Abs(y-0.00)>3.40) d+=dd; + d+=0.0034; + break; + } + if (TMath::Abs(z+3.9*(i+0.5))<0.15) { + if (TMath::Abs(y-0.00)>3.40) d+=dd; + d+=0.0034; + break; + } + if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;} + if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;} + } + } else + if (373.40) d+=dd; + if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);} + if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);} + for (Int_t i=0; i<11; i++) { + if (TMath::Abs(z-3.9*i)<0.15) { + if (TMath::Abs(y-0.00)>3.40) d+=dd; + d+=dd; + break; + } + if (TMath::Abs(z+3.9*i)<0.15) { + if (TMath::Abs(y-0.00)>3.40) d+=dd; + d+=dd; + break; + } + if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;} + if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;} + } + } else + if (133.30) d+=dd; + + if (TMath::Abs(y-1.80)<0.55) { + d+=0.016; + for (Int_t j=0; j<20; j++) { + if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;} + if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;} + } + } + if (TMath::Abs(y+1.80)<0.55) { + d+=0.016; + for (Int_t j=0; j<20; j++) { + if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;} + if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;} + } + } + for (Int_t i=0; i<4; i++) { + if (TMath::Abs(z-7.3*i)<0.60) { + d+=dd; + if (TMath::Abs(y-0.00)>3.30) d+=dd; + break; + } + if (TMath::Abs(z+7.3*i)<0.60) { + d+=dd; + if (TMath::Abs(y-0.00)>3.30) d+=dd; + break; + } + } + } else + if (60.5) d+=dd; + //if (TMath::Abs(y-3.08)>0.45) d+=dd; + if (TMath::Abs(y-3.03)<0.10) {d+=0.014;} + } else + if (30.6) d+=dd; + //if (TMath::Abs(y+0.21)>0.45) d+=dd; + if (TMath::Abs(y+0.10)<0.10) {d+=0.014;} + } -Double_t AliITStrackerV2::GetEffectiveThickness(Double_t phi,Double_t z) const + return d; +} + +Double_t AliITStrackerV2::GetEffectiveThickness(Double_t y,Double_t z) const { //-------------------------------------------------------------------- - //Returns the thickness between the current layer and the vertex + //Returns the thickness between the current layer and the vertex (units X0) //-------------------------------------------------------------------- - //-pi1) { - Double_t xi=0.5*(fLayers[1].GetR()+fLayers[2].GetR()); - d+=0.034*xi*xi; + Double_t xi=9.; + d+=0.0097*xi*xi; } if (fI>3) { - Double_t xi=0.5*(fLayers[3].GetR()+fLayers[4].GetR()); - d+=0.039*xi*xi; + Double_t xi=0.5*(fgLayers[3].GetR()+fgLayers[4].GetR()); + d+=0.0034*xi*xi; } + return d/(xn*xn); } - - -Int_t AliITStrackerV2::AliITSlayer::InRoad() const { +Bool_t AliITStrackerV2::RefitAt(Double_t xx,AliITStrackV2 *t, + const AliITStrackV2 *c, Bool_t extra) { //-------------------------------------------------------------------- - // This function returns number of clusters within the "window" + // 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 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(); + Int_t index[AliITSgeomTGeo::kNLayers]; + Int_t k; + for (k=0; kGetNumberOfClusters(); + for (k=0; kGetClusterIndex(k),nl=(idx&0xf0000000)>>28; + index[nl]=idx; + } + + Int_t from, to, step; + if (xx > t->GetX()) { + from=0; to=AliITSgeomTGeo::GetNLayers(); + step=+1; + } else { + from=AliITSgeomTGeo::GetNLayers()-1; to=-1; + step=-1; + } + + for (Int_t i=from; i != to; i += step) { + AliITSlayer &layer=fgLayers[i]; + Double_t r=layer.GetR(); + + { + Double_t hI=i-0.5*step; + if (TMath::Abs(hI-1.5)<0.01 || TMath::Abs(hI-3.5)<0.01) { + Double_t rs=0.5*(fgLayers[i-step].GetR() + r); + Double_t d=0.0034, x0=38.6; + if (TMath::Abs(hI-1.5)<0.01) {rs=9.; d=0.0097; x0=42;} + if (!t->PropagateTo(rs,-step*d,x0)) { + return kFALSE; + } + } + } - if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi(); - if (y>1.*fR*TMath::Pi() && fYmaxIsStartedTimeIntegral() && step==1) { + t->GetGlobalXYZat(t->GetX(),oldX,oldY,oldZ); + } + // + + Double_t phi,z; + if (!t->GetPhiZat(r,phi,z)) { + return kFALSE; + } + + Int_t idet=layer.FindDetectorIndex(phi,z); + if (idet<0) { + return kFALSE; + } + const AliITSdetector &det=layer.GetDetector(idet); + phi=det.GetPhi(); + if (!t->Propagate(phi,det.GetR())) { + return kFALSE; + } + t->SetDetectorIndex(idet); + + const AliITSRecPoint *cl=0; + Double_t maxchi2=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); + + Int_t idx=index[i]; + 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(ccc); + if (chi2GetX()+cl->GetX(); + if (!t->PropagateTo(x,0.,0.)) return kFALSE; + if (!t->Update(cl,maxchi2,idx)) { + return kFALSE; + } + t->SetSampledEdx(cl->GetQ(),i-2); + } + + { + Double_t x0; + Double_t d=layer.GetThickness(t->GetY(),t->GetZ(),x0); + 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; + t->GetGlobalXYZat(t->GetX(),newX,newY,newZ); + Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) + + (oldZ-newZ)*(oldZ-newZ); + t->AddTimeStep(TMath::Sqrt(dL2)); + } + // - if (yfYmax) continue; - ncl++; } - return ncl; + + if (!t->PropagateTo(xx,0.,0.)) return kFALSE; + return kTRUE; } +void AliITStrackerV2::UseClusters(const AliKalmanTrack *t, Int_t from) const { + //-------------------------------------------------------------------- + // This function marks clusters assigned to the track + //-------------------------------------------------------------------- + AliTracker::UseClusters(t,from); + + Int_t clusterIndex = t->GetClusterIndex(0); + AliITSRecPoint *c= 0x0; + + if (clusterIndex>-1) + c = (AliITSRecPoint *)GetCluster(clusterIndex); + if (c && c->GetSigmaZ2()>0.1) c->UnUse(); + c = 0x0; + clusterIndex = t->GetClusterIndex(1); + if (clusterIndex>-1) + c=(AliITSRecPoint *)GetCluster(clusterIndex); + if (c && c->GetSigmaZ2()>0.1) c->UnUse(); +}