X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITStrackV2.cxx;h=f62e7c63a9aa8b3b062357fc71e6d010b77a1b94;hb=0e73ed90793f328b8f0fbeebb2e2d0994ccf4ddc;hp=ec1fb078b40b47a1d224519c5550f8d9012c1356;hpb=6c4ef2ed40aa03c28c60811a781da1bb8c545974;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/AliITStrackV2.cxx b/ITS/AliITStrackV2.cxx index ec1fb078b40..f62e7c63a9a 100644 --- a/ITS/AliITStrackV2.cxx +++ b/ITS/AliITStrackV2.cxx @@ -13,6 +13,8 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ +/* $Id$ */ + /////////////////////////////////////////////////////////////////////////// // Implementation of the ITS track class // @@ -22,27 +24,32 @@ #include #include "AliCluster.h" -#include "AliTracker.h" #include "AliESDtrack.h" +#include "AliESDVertex.h" +#include "AliITSReconstructor.h" #include "AliITStrackV2.h" +#include "AliTracker.h" + +const Int_t AliITStrackV2::fgkWARN = 5; ClassImp(AliITStrackV2) -const Int_t kWARN=5; //____________________________________________________________________________ AliITStrackV2::AliITStrackV2() : AliKalmanTrack(), + fCheckInvariant(kTRUE), fdEdx(0), fESDtrack(0) { - for(Int_t i=0; i<2*kMaxLayer; i++) fIndex[i]=-1; - for(Int_t i=0; i<4; i++) fdEdxSample[i]=0; + for(Int_t i=0; i<2*AliITSgeomTGeo::kNLayers; i++) {fIndex[i]=-1; fModule[i]=-1;} + for(Int_t i=0; i<4; i++) fdEdxSample[i]=0; } //____________________________________________________________________________ AliITStrackV2::AliITStrackV2(AliESDtrack& t,Bool_t c) throw (const Char_t *) : AliKalmanTrack(), + fCheckInvariant(kTRUE), fdEdx(t.GetITSsignal()), fESDtrack(&t) { @@ -72,13 +79,35 @@ AliITStrackV2::AliITStrackV2(AliESDtrack& t,Bool_t c) throw (const Char_t *) : for(Int_t i=0; i<4; i++) fdEdxSample[i]=0; } +//____________________________________________________________________________ +void AliITStrackV2::ResetClusters() { + //------------------------------------------------------------------ + // Reset the array of attached clusters. + //------------------------------------------------------------------ + for (Int_t i=0; i<2*AliITSgeomTGeo::kNLayers; i++) fIndex[i]=-1; + SetChi2(0.); + SetNumberOfClusters(0); +} + +//____________________________________________________________________________ void AliITStrackV2::UpdateESDtrack(ULong_t flags) const { + // Update track params fESDtrack->UpdateTrackParams(this,flags); + // copy the module indices + for(Int_t i=0;i<12;i++) { + // printf(" %d\n",GetModuleIndex(i)); + fESDtrack->SetITSModuleIndex(i,GetModuleIndex(i)); + } + // copy the 4 dedx samples + Double_t sdedx[4]={0.,0.,0.,0.}; + for(Int_t i=0; i<4; i++) sdedx[i]=fdEdxSample[i]; + fESDtrack->SetITSdEdxSamples(sdedx); } //____________________________________________________________________________ AliITStrackV2::AliITStrackV2(const AliITStrackV2& t) : AliKalmanTrack(t), + fCheckInvariant(t.fCheckInvariant), fdEdx(t.fdEdx), fESDtrack(t.fESDtrack) { @@ -87,7 +116,10 @@ AliITStrackV2::AliITStrackV2(const AliITStrackV2& t) : //------------------------------------------------------------------ Int_t i; for (i=0; i<4; i++) fdEdxSample[i]=t.fdEdxSample[i]; - for (i=0; i<2*kMaxLayer; i++) fIndex[i]=t.fIndex[i]; + for (i=0; i<2*AliITSgeomTGeo::GetNLayers(); i++) { + fIndex[i]=t.fIndex[i]; + fModule[i]=t.fModule[i]; + } } //_____________________________________________________________________________ @@ -96,8 +128,8 @@ Int_t AliITStrackV2::Compare(const TObject *o) const { // This function compares tracks according to the their curvature //----------------------------------------------------------------- AliITStrackV2 *t=(AliITStrackV2*)o; - //Double_t co=TMath::Abs(t->Get1Pt()); - //Double_t c =TMath::Abs(Get1Pt()); + //Double_t co=OneOverPt(); + //Double_t c =OneOverPt(); Double_t co=t->GetSigmaY2()*t->GetSigmaZ2(); Double_t c =GetSigmaY2()*GetSigmaZ2(); if (c>co) return 1; @@ -113,19 +145,22 @@ AliITStrackV2::PropagateToVertex(const AliESDVertex *v,Double_t d,Double_t x0) //This function propagates a track to the minimal distance from the origin //------------------------------------------------------------------ Double_t bz=GetBz(); - if (PropagateToDCA(v,bz,kVeryBig)) - if (AliExternalTrackParam::CorrectForMaterial(d,x0,GetMass())) return kTRUE; + if (PropagateToDCA(v,bz,kVeryBig)) { + Double_t xOverX0,xTimesRho; + xOverX0 = d; xTimesRho = d*x0; + if (CorrectForMeanMaterial(xOverX0,xTimesRho,kTRUE)) return kTRUE; + } return kFALSE; } //____________________________________________________________________________ Bool_t AliITStrackV2:: -GetGlobalXYZat(Double_t xk, Double_t &x, Double_t &y, Double_t &z) const { +GetGlobalXYZat(Double_t xloc, Double_t &x, Double_t &y, Double_t &z) const { //------------------------------------------------------------------ //This function returns a track position in the global system //------------------------------------------------------------------ Double_t r[3]; - Bool_t rc=GetXYZAt(xk, AliTracker::GetBz(), r); + Bool_t rc=GetXYZAt(xloc, GetBz(), r); x=r[0]; y=r[1]; z=r[2]; return rc; } @@ -145,13 +180,18 @@ Bool_t AliITStrackV2::PropagateTo(Double_t xk, Double_t d, Double_t x0) { //------------------------------------------------------------------ //This function propagates a track //------------------------------------------------------------------ + Double_t oldX=GetX(), oldY=GetY(), oldZ=GetZ(); - Double_t bz=GetBz(); - if (!AliExternalTrackParam::PropagateTo(xk,bz)) return kFALSE; - if (!AliExternalTrackParam::CorrectForMaterial(d,x0,GetMass())) return kFALSE; - - Double_t x=GetX(), y=GetZ(), z=GetZ(); + //Double_t bz=GetBz(); + //if (!AliExternalTrackParam::PropagateTo(xk,bz)) return kFALSE; + Double_t b[3]; GetBxByBz(b); + if (!AliExternalTrackParam::PropagateToBxByBz(xk,b)) return kFALSE; + Double_t xOverX0,xTimesRho; + xOverX0 = d; xTimesRho = d*x0; + if (!CorrectForMeanMaterial(xOverX0,xTimesRho,kTRUE)) return kFALSE; + + Double_t x=GetX(), y=GetY(), z=GetZ(); if (IsStartedTimeIntegral() && x>oldX) { Double_t l2 = (x-oldX)*(x-oldX) + (y-oldY)*(y-oldY) + (z-oldZ)*(z-oldZ); AddTimeStep(TMath::Sqrt(l2)); @@ -160,6 +200,53 @@ Bool_t AliITStrackV2::PropagateTo(Double_t xk, Double_t d, Double_t x0) { return kTRUE; } +//____________________________________________________________________________ +Bool_t AliITStrackV2::PropagateToTGeo(Double_t xToGo, Int_t nstep, Double_t &xOverX0, Double_t &xTimesRho, Bool_t addTime) { + //------------------------------------------------------------------- + // Propagates the track to a reference plane x=xToGo in n steps. + // These n steps are only used to take into account the curvature. + // The material is calculated with TGeo. (L.Gaudichet) + //------------------------------------------------------------------- + + Double_t startx = GetX(), starty = GetY(), startz = GetZ(); + Double_t sign = (startxstartx) { + Double_t l2 = ( (GetX()-startx)*(GetX()-startx) + + (GetY()-starty)*(GetY()-starty) + + (GetZ()-startz)*(GetZ()-startz) ); + AddTimeStep(TMath::Sqrt(l2)); + } + + return kTRUE; +} + //____________________________________________________________________________ Bool_t AliITStrackV2::Update(const AliCluster* c, Double_t chi2, Int_t index) { @@ -167,18 +254,25 @@ Bool_t AliITStrackV2::Update(const AliCluster* c, Double_t chi2, Int_t index) //This function updates track parameters //------------------------------------------------------------------ Double_t p[2]={c->GetY(), c->GetZ()}; - Double_t cov[3]={c->GetSigmaY2(), 0., c->GetSigmaZ2()}; + Double_t cov[3]={c->GetSigmaY2(), c->GetSigmaYZ(), c->GetSigmaZ2()}; if (!AliExternalTrackParam::Update(p,cov)) return kFALSE; + Int_t n=GetNumberOfClusters(); if (!Invariant()) { - AliWarning("Wrong invariant !"); + if (n>fgkWARN) AliWarning("Wrong invariant !"); return kFALSE; } if (chi2<0) return kTRUE; - Int_t n=GetNumberOfClusters(); + // fill residuals for ITS+TPC tracks + if (fESDtrack) { + if (fESDtrack->GetStatus()&AliESDtrack::kTPCin) { + AliTracker::FillResiduals(this,p,cov,c->GetVolumeId()); + } + } + fIndex[n]=index; SetNumberOfClusters(n+1); SetChi2(GetChi2()+chi2); @@ -190,36 +284,51 @@ Bool_t AliITStrackV2::Invariant() const { //------------------------------------------------------------------ // This function is for debugging purpose only //------------------------------------------------------------------ + if(!fCheckInvariant) return kTRUE; + Int_t n=GetNumberOfClusters(); + // take into account the misalignment error + Float_t maxMisalErrY2=0,maxMisalErrZ2=0; + for (Int_t lay=0; layGetClusterMisalErrorY(lay,GetBz())); + maxMisalErrZ2 = TMath::Max(maxMisalErrZ2,AliITSReconstructor::GetRecoParam()->GetClusterMisalErrorZ(lay,GetBz())); + } + maxMisalErrY2 *= maxMisalErrY2; + maxMisalErrZ2 *= maxMisalErrZ2; + // this is because when we reset before refitting, we multiply the + // matrix by 10 + maxMisalErrY2 *= 10.; + maxMisalErrZ2 *= 10.; + Double_t sP2=GetParameter()[2]; if (TMath::Abs(sP2) >= kAlmost1){ - if (n>kWARN) Warning("Invariant","fP2=%f\n",sP2); + if (n>fgkWARN) Warning("Invariant","fP2=%f\n",sP2); return kFALSE; } Double_t sC00=GetCovariance()[0]; - if (sC00<=0 || sC00>9.) { - if (n>kWARN) Warning("Invariant","fC00=%f\n",sC00); + if (sC00<=0 || sC00>(9.+maxMisalErrY2)) { + if (n>fgkWARN) Warning("Invariant","fC00=%f\n",sC00); return kFALSE; } Double_t sC11=GetCovariance()[2]; - if (sC11<=0 || sC11>9.) { - if (n>kWARN) Warning("Invariant","fC11=%f\n",sC11); + if (sC11<=0 || sC11>(9.+maxMisalErrZ2)) { + if (n>fgkWARN) Warning("Invariant","fC11=%f\n",sC11); return kFALSE; } Double_t sC22=GetCovariance()[5]; if (sC22<=0 || sC22>1.) { - if (n>kWARN) Warning("Invariant","fC22=%f\n",sC22); + if (n>fgkWARN) Warning("Invariant","fC22=%f\n",sC22); return kFALSE; } Double_t sC33=GetCovariance()[9]; if (sC33<=0 || sC33>1.) { - if (n>kWARN) Warning("Invariant","fC33=%f\n",sC33); + if (n>fgkWARN) Warning("Invariant","fC33=%f\n",sC33); return kFALSE; } Double_t sC44=GetCovariance()[14]; if (sC44<=0 /*|| sC44>6e-5*/) { - if (n>kWARN) Warning("Invariant","fC44=%f\n",sC44); + if (n>fgkWARN) Warning("Invariant","fC44=%f\n",sC44); return kFALSE; } @@ -231,12 +340,51 @@ Bool_t AliITStrackV2::Propagate(Double_t alp,Double_t xk) { //------------------------------------------------------------------ //This function propagates a track //------------------------------------------------------------------ - Double_t bz=GetBz(); - if (!AliExternalTrackParam::Propagate(alp,xk,bz)) return kFALSE; + //Double_t bz=GetBz(); + //if (!AliExternalTrackParam::Propagate(alp,xk,bz)) return kFALSE; + Double_t b[3]; GetBxByBz(b); + if (!AliExternalTrackParam::PropagateBxByBz(alp,xk,b)) return kFALSE; if (!Invariant()) { - AliWarning("Wrong invariant !"); - return kFALSE; + Int_t n=GetNumberOfClusters(); + if (n>fgkWARN) AliWarning("Wrong invariant !"); + return kFALSE; + } + + return kTRUE; +} + +Bool_t AliITStrackV2::MeanBudgetToPrimVertex(Double_t xyz[3], Double_t step, Double_t &d) const { + + //------------------------------------------------------------------- + // Get the mean material budget between the actual point and the + // primary vertex. (L.Gaudichet) + //------------------------------------------------------------------- + + Double_t cs=TMath::Cos(GetAlpha()), sn=TMath::Sin(GetAlpha()); + Double_t vertexX = xyz[0]*cs + xyz[1]*sn; + + Int_t nstep = Int_t((GetX()-vertexX)/step); + if (nstep<1) nstep = 1; + step = (GetX()-vertexX)/nstep; + + // Double_t mparam[7], densMean=0, radLength=0, length=0; + Double_t mparam[7]; + Double_t p1[3], p2[3], x = GetX(), bz = GetBz(); + GetXYZ(p1); + + d=0.; + + for (Int_t i=0; i900000) return kFALSE; + d += mparam[1]; + + p1[0] = p2[0]; + p1[1] = p2[1]; + p1[2] = p2[2]; } return kTRUE; @@ -246,14 +394,40 @@ Bool_t AliITStrackV2::Improve(Double_t x0,Double_t xyz[3],Double_t ers[3]) { //------------------------------------------------------------------ //This function improves angular track parameters //------------------------------------------------------------------ + //Store the initail track parameters + + return kTRUE; //PH temporary switched off + + Double_t x = GetX(); + Double_t alpha = GetAlpha(); + Double_t par[] = {GetY(),GetZ(),GetSnp(),GetTgl(),GetSigned1Pt()}; + Double_t cov[] = { + GetSigmaY2(), + GetSigmaZY(), + GetSigmaZ2(), + GetSigmaSnpY(), + GetSigmaSnpZ(), + GetSigmaSnp2(), + GetSigmaTglY(), + GetSigmaTglZ(), + GetSigmaTglSnp(), + GetSigmaTgl2(), + GetSigma1PtY(), + GetSigma1PtZ(), + GetSigma1PtSnp(), + GetSigma1PtTgl(), + GetSigma1Pt2() + }; + + Double_t cs=TMath::Cos(GetAlpha()), sn=TMath::Sin(GetAlpha()); //Double_t xv = xyz[0]*cs + xyz[1]*sn; // vertex Double_t yv =-xyz[0]*sn + xyz[1]*cs; // in the Double_t zv = xyz[2]; // local frame - Double_t dy = Par(0) - yv, dz = Par(1) - zv; + Double_t dy = par[0] - yv, dz = par[1] - zv; Double_t r2=GetX()*GetX() + dy*dy; - Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()); + Double_t p2=(1.+ GetTgl()*GetTgl())/(GetSigned1Pt()*GetSigned1Pt()); Double_t beta2=p2/(p2 + GetMass()*GetMass()); x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp())); Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0; @@ -264,29 +438,34 @@ Bool_t AliITStrackV2::Improve(Double_t x0,Double_t xyz[3],Double_t ers[3]) { Double_t dummy = 4/r2 - GetC()*GetC(); if (dummy < 0) return kFALSE; Double_t parp = 0.5*(GetC()*GetX() + dy*TMath::Sqrt(dummy)); - Double_t sigma2p = theta2*(1.- GetSnp()*GetSnp())*(1. + GetTgl()*GetTgl()); - sigma2p += Cov(0)/r2*(1.- dy*dy/r2)*(1.- dy*dy/r2); + Double_t sigma2p = theta2*(1.-GetSnp())*(1.+GetSnp())*(1. + GetTgl()*GetTgl()); + Double_t ovSqr2 = 1./TMath::Sqrt(r2); + Double_t tfact = ovSqr2*(1.-dy*ovSqr2)*(1.+dy*ovSqr2); + sigma2p += cov[0]*tfact*tfact; sigma2p += ers[1]*ers[1]/r2; - sigma2p += 0.25*Cov(14)*cnv*cnv*GetX()*GetX(); - Double_t eps2p=sigma2p/(Cov(2) + sigma2p); - Par(0) += Cov(3)/(Cov(5) + sigma2p)*(parp - GetSnp()); - Par(2) = eps2p*GetSnp() + (1 - eps2p)*parp; - Cov(5) *= eps2p; - Cov(3) *= eps2p; + sigma2p += 0.25*cov[14]*cnv*cnv*GetX()*GetX(); + Double_t eps2p=sigma2p/(cov[5] + sigma2p); + par[0] += cov[3]/(cov[5] + sigma2p)*(parp - GetSnp()); + par[2] = eps2p*GetSnp() + (1 - eps2p)*parp; + cov[5] *= eps2p; + cov[3] *= eps2p; } { Double_t parl=0.5*GetC()*dz/TMath::ASin(0.5*GetC()*TMath::Sqrt(r2)); Double_t sigma2l=theta2; - sigma2l += Cov(2)/r2 + Cov(0)*dy*dy*dz*dz/(r2*r2*r2); + sigma2l += cov[2]/r2 + cov[0]*dy*dy*dz*dz/(r2*r2*r2); sigma2l += ers[2]*ers[2]/r2; - Double_t eps2l = sigma2l/(Cov(9) + sigma2l); - Par(1) += Cov(7 )/(Cov(9) + sigma2l)*(parl - Par(3)); - Par(4) += Cov(13)/(Cov(9) + sigma2l)*(parl - Par(3)); - Par(3) = eps2l*Par(3) + (1-eps2l)*parl; - Cov(9) *= eps2l; - Cov(13)*= (eps2l/cnv/cnv); - Cov(7) *= eps2l; + Double_t eps2l = sigma2l/(cov[9] + sigma2l); + par[1] += cov[7 ]/(cov[9] + sigma2l)*(parl - par[3]); + par[4] += cov[13]/(cov[9] + sigma2l)*(parl - par[3]); + par[3] = eps2l*par[3] + (1-eps2l)*parl; + cov[9] *= eps2l; + cov[13]*= eps2l; + cov[7] *= eps2l; } + + Set(x,alpha,par,cov); + if (!Invariant()) return kFALSE; return kTRUE; @@ -296,44 +475,106 @@ void AliITStrackV2::CookdEdx(Double_t low, Double_t up) { //----------------------------------------------------------------- // This function calculates dE/dX within the "low" and "up" cuts. // Origin: Boris Batyunya, JINR, Boris.Batiounia@cern.ch + // Updated: F. Prino 8-June-2009 //----------------------------------------------------------------- - // The clusters order is: SSD-2, SSD-1, SDD-2, SDD-1, SPD-2, SPD-1 + // The cluster order is: SDD-1, SDD-2, SSD-1, SSD-2 - Int_t i; Int_t nc=0; - for (i=0; i>28; - if (idx>1) nc++; // Take only SSD and SDD + Float_t dedx[4]; + for (Int_t il=0; il<4; il++) { // count good (>0) dE/dx values + if(fdEdxSample[il]>0.){ + dedx[nc]= fdEdxSample[il]; + nc++; + } + } + if(nc<1){ + SetdEdx(0.); + return; } - Int_t swap;//stupid sorting + Int_t swap; // sort in ascending order do { swap=0; - for (i=0; i0) dedx /= (nu-nl); - SetdEdx(dedx); + Double_t nl=low*nc, nu =up*nc; + Float_t sumamp = 0; + Float_t sumweight =0; + for (Int_t i=0; inu-1) weight = TMath::Max(nu-i,0.); + sumamp+= dedx[i]*weight; + sumweight+=weight; + } + SetdEdx(sumamp/sumweight); } -Double_t AliITStrackV2::GetBz() const { - // - // returns Bz component of the magnetic field (kG) - // - if (AliTracker::UniformField()) return AliTracker::GetBz(); - Double_t r[3]; GetXYZ(r); - return AliTracker::GetBz(r); +//____________________________________________________________________________ +Bool_t AliITStrackV2:: +GetPhiZat(Double_t r, Double_t &phi, Double_t &z) const { + //------------------------------------------------------------------ + // This function returns the global cylindrical (phi,z) of the track + // position estimated at the radius r. + // The track curvature is neglected. + //------------------------------------------------------------------ + Double_t d=GetD(0.,0.); + if (TMath::Abs(d) > r) { + if (r>1e-1) return kFALSE; + r = TMath::Abs(d); + } + + Double_t rcurr=TMath::Sqrt(GetX()*GetX() + GetY()*GetY()); + if (TMath::Abs(d) > rcurr) return kFALSE; + Double_t globXYZcurr[3]; GetXYZ(globXYZcurr); + Double_t phicurr=TMath::ATan2(globXYZcurr[1],globXYZcurr[0]); + + if (GetX()>=0.) { + phi=phicurr+TMath::ASin(d/r)-TMath::ASin(d/rcurr); + } else { + phi=phicurr+TMath::ASin(d/r)+TMath::ASin(d/rcurr)-TMath::Pi(); + } + + // return a phi in [0,2pi[ + if (phi<0.) phi+=2.*TMath::Pi(); + else if (phi>=2.*TMath::Pi()) phi-=2.*TMath::Pi(); + z=GetZ()+GetTgl()*(TMath::Sqrt((r-d)*(r+d))-TMath::Sqrt((rcurr-d)*(rcurr+d))); + return kTRUE; } +//____________________________________________________________________________ +Bool_t AliITStrackV2:: +GetLocalXat(Double_t r,Double_t &xloc) const { + //------------------------------------------------------------------ + // This function returns the local x of the track + // position estimated at the radius r. + // The track curvature is neglected. + //------------------------------------------------------------------ + Double_t d=GetD(0.,0.); + if (TMath::Abs(d) > r) { + if (r>1e-1) return kFALSE; + r = TMath::Abs(d); + } + + Double_t rcurr=TMath::Sqrt(GetX()*GetX() + GetY()*GetY()); + Double_t globXYZcurr[3]; GetXYZ(globXYZcurr); + Double_t phicurr=TMath::ATan2(globXYZcurr[1],globXYZcurr[0]); + Double_t phi; + if (GetX()>=0.) { + phi=phicurr+TMath::ASin(d/r)-TMath::ASin(d/rcurr); + } else { + phi=phicurr+TMath::ASin(d/r)+TMath::ASin(d/rcurr)-TMath::Pi(); + } + + xloc=r*(TMath::Cos(phi)*TMath::Cos(GetAlpha()) + +TMath::Sin(phi)*TMath::Sin(GetAlpha())); + return kTRUE; +}