X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TRD%2FAliTRDtrack.cxx;h=4df4749a96c677a1fffd8b2b275a94049c2e6bf5;hb=d161aeeb2d9dd6fb9acb53467e8c62a2a8b96cfe;hp=c9f11d1b665f343574c9965f1aaeed0426bb66ae;hpb=a819a5f7bddaac1e6625eb9531e9e0734bba8c83;p=u%2Fmrichter%2FAliRoot.git diff --git a/TRD/AliTRDtrack.cxx b/TRD/AliTRDtrack.cxx index c9f11d1b665..4df4749a96c 100644 --- a/TRD/AliTRDtrack.cxx +++ b/TRD/AliTRDtrack.cxx @@ -13,463 +13,1171 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -Revision 1.4 2000/12/08 16:07:02 cblume -Update of the tracking by Sergei +/* $Id$ */ -Revision 1.3 2000/10/15 23:40:01 cblume -Remove AliTRDconst +#include -Revision 1.2 2000/10/06 16:49:46 cblume -Made Getters const +#include "AliTracker.h" -Revision 1.1.2.1 2000/09/22 14:47:52 cblume -Add the tracking code - -*/ - -#include - -#include - -#include "AliTRD.h" #include "AliTRDgeometry.h" #include "AliTRDcluster.h" #include "AliTRDtrack.h" +#include "AliTRDcalibDB.h" +#include "Cal/AliTRDCalPID.h" ClassImp(AliTRDtrack) +/////////////////////////////////////////////////////////////////////////////// +// // +// Represents a reconstructed TRD track // +// Local TRD Kalman track // +// // +/////////////////////////////////////////////////////////////////////////////// + +//_____________________________________________________________________________ +AliTRDtrack::AliTRDtrack() + :AliKalmanTrack() + ,fSeedLab(-1) + ,fdEdx(0) + ,fDE(0) + ,fPIDquality(0) + ,fClusterOwner(kFALSE) + ,fPIDmethod(kLQ) + ,fStopped(kFALSE) + ,fLhElectron(0) + ,fNWrong(0) + ,fNRotate(0) + ,fNCross(0) + ,fNExpected(0) + ,fNLast(0) + ,fNExpectedLast(0) + ,fNdedx(0) + ,fChi2Last(1e10) + ,fBackupTrack(0x0) +{ + // + // AliTRDtrack default constructor + // + + for (Int_t i = 0; i < kNplane; i++) { + for (Int_t j = 0; j < kNslice; j++) { + fdEdxPlane[i][j] = 0.0; + } + fTimBinPlane[i] = -1; + fMom[i] = -1.; + fSnp[i] = 0.; + fTgl[i] = 0.; + fTrackletIndex[i] = -1; + } + for(int ispec=0; ispecGetQ()); + Double_t s = GetSnp(); + Double_t t = GetTgl(); + if (s*s < 1) { + q *= TMath::Sqrt((1-s*s)/(1+t*t)); + } + + fdQdl[0] = q; + for (UInt_t i = 1; i < kMAXCLUSTERSPERTRACK; i++) { + fdQdl[i] = 0; + fIndex[i] = 0; + fIndexBackup[i] = 0; + fClusters[i] = 0x0; + } + + for (Int_t i = 0; i < 3;i++) { + fBudget[i] = 0; + } -AliTRDtrack::AliTRDtrack(const AliTRDcluster *c, UInt_t index, -const Double_t xx[5], const Double_t cc[15], Double_t xref, Double_t alpha) { - //----------------------------------------------------------------- - // This is the main track constructor. - //----------------------------------------------------------------- - fLab=-1; - fChi2=0.; - fdEdx=0.; - - fAlpha=alpha; - fX=xref; - - fY=xx[0]; fZ=xx[1]; fC=xx[2]; fE=xx[3]; fT=xx[4]; - - fCyy=cc[0]; - fCzy=cc[1]; fCzz=cc[2]; - fCcy=cc[3]; fCcz=cc[4]; fCcc=cc[5]; - fCey=cc[6]; fCez=cc[7]; fCec=cc[8]; fCee=cc[9]; - fCty=cc[10]; fCtz=cc[11]; fCtc=cc[12]; fCte=cc[13]; fCtt=cc[14]; - - fN=0; - fIndex[fN]=index; - - Float_t q = c->GetQ(); - Double_t s = fX*fC - fE, t=fT; - q *= TMath::Sqrt((1-s*s)/(1+t*t)); - fdQdl[fN++] = q; } //_____________________________________________________________________________ -AliTRDtrack::AliTRDtrack(const AliTRDtrack& t) { +AliTRDtrack::AliTRDtrack(const AliTRDtrack &t/*, const Bool_t owner*/) + :AliKalmanTrack(t) + ,fSeedLab(t.GetSeedLabel()) + ,fdEdx(t.fdEdx) + ,fDE(t.fDE) + ,fPIDquality(t.fPIDquality) + ,fClusterOwner(kTRUE) + ,fPIDmethod(t.fPIDmethod) + ,fStopped(t.fStopped) + ,fLhElectron(0) + ,fNWrong(t.fNWrong) + ,fNRotate(t.fNRotate) + ,fNCross(t.fNCross) + ,fNExpected(t.fNExpected) + ,fNLast(t.fNLast) + ,fNExpectedLast(t.fNExpectedLast) + ,fNdedx(t.fNdedx) + ,fChi2Last(t.fChi2Last) + ,fBackupTrack(0x0) +{ // // Copy constructor. // - fLab=t.fLab; + for (Int_t i = 0; i < kNplane; i++) { + for (Int_t j = 0; j < kNslice; j++) { + fdEdxPlane[i][j] = t.fdEdxPlane[i][j]; + } + fTimBinPlane[i] = t.fTimBinPlane[i]; + fTracklets[i] = t.fTracklets[i]; + fMom[i] = t.fMom[i]; + fSnp[i] = t.fSnp[i]; + fTgl[i] = t.fTgl[i]; + fTrackletIndex[i] = t.fTrackletIndex[i]; + } + + Int_t n = t.GetNumberOfClusters(); + SetNumberOfClusters(n); - fChi2=t.fChi2; - fdEdx=t.fdEdx; + for (Int_t i = 0; i < n; i++) { + fIndex[i] = t.fIndex[i]; + fIndexBackup[i] = t.fIndex[i]; + fdQdl[i] = t.fdQdl[i]; + if (fClusterOwner && t.fClusters[i]) { + fClusters[i] = new AliTRDcluster(*(t.fClusters[i])); + } + else { + fClusters[i] = t.fClusters[i]; + } + } - fAlpha=t.fAlpha; - fX=t.fX; + for (UInt_t i = n; i < kMAXCLUSTERSPERTRACK; i++) { + fdQdl[i] = 0; + fIndex[i] = 0; + fIndexBackup[i] = 0; + fClusters[i] = 0x0; + } - fY=t.fY; fZ=t.fZ; fC=t.fC; fE=t.fE; fT=t.fT; + for (Int_t i = 0; i < 3;i++) { + fBudget[i] = t.fBudget[i]; + } - fCyy=t.fCyy; - fCzy=t.fCzy; fCzz=t.fCzz; - fCcy=t.fCcy; fCcz=t.fCcz; fCcc=t.fCcc; - fCey=t.fCey; fCez=t.fCez; fCec=t.fCec; fCee=t.fCee; - fCty=t.fCty; fCtz=t.fCtz; fCtc=t.fCtc; fCte=t.fCte; fCtt=t.fCtt; + for(Int_t ispec = 0; ispecGetX() + ,par->GetAlpha() + ,par->GetParameter() + ,par->GetCovariance()); + + for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) { + fdQdl[i] = 0; + fClusters[i] = 0x0; + } + + for (Int_t i = 0; i < 3; i++) { + fBudget[i] = 0; + } + for(int ispec=0; ispecGetSigmaY2(); - // Double_t c =GetSigmaY2(); + Double_t co = TMath::Abs(t->GetC()); + Double_t c = TMath::Abs(GetC()); - Double_t co=TMath::Abs(t->GetC()); - Double_t c =TMath::Abs(GetC()); + if (c > co) { + return 1; + } + else if (c < co) { + return -1; + } - if (c>co) return 1; - else if (cGetNumberOfTimeBins(); + Int_t plane; // Plane of current cluster + Int_t tb; // Time bin of current cluster + Int_t slice; // Current slice + AliTRDcluster *cluster = 0x0; // Pointer to current cluster + + // Reset class and local counters/variables + for (Int_t iPlane = 0; iPlane < kNplane; iPlane++) { + fTimBinPlane[iPlane] = -1; + maxcharge[iPlane] = 0.0; + for (Int_t iSlice = 0; iSlice < kNslice; iSlice++) { + fdEdxPlane[iPlane][iSlice] = 0.0; + nCluster[iPlane][iSlice] = 0; + } + } + + // Start looping over clusters attached to this track + for (Int_t iClus = 0; iClus < GetNumberOfClusters(); iClus++) { + + cluster = fClusters[iClus]; //(AliTRDcluster*)tracker->GetCluster(fIndex[iClus]); + if (!cluster) continue; + + // Read info from current cluster + plane = AliTRDgeometry::GetLayer(cluster->GetDetector()); + if (plane < 0 || plane >= kNplane) { + AliError(Form("Wrong plane %d", plane)); + continue; + } + + tb = cluster->GetLocalTimeBin(); + if ((tb < 0) || (tb >= ntb)) { + AliWarning(Form("time bin < 0 or > %d in cluster %d", ntb, iClus)); + AliInfo(Form("dQ/dl %f", fdQdl[iClus])); + continue; + } + + slice = tb * kNslice / ntb; + fdEdxPlane[plane][slice] += fdQdl[iClus]; + if (fdQdl[iClus] > maxcharge[plane]) { + maxcharge[plane] = fdQdl[iClus]; + fTimBinPlane[plane] = tb; + } + + nCluster[plane][slice]++; + + } // End of loop over cluster + + // Normalize fdEdxPlane to number of clusters and set track segments + for (Int_t iPlane = 0; iPlane < kNplane; iPlane++) { + for (Int_t iSlice = 0; iSlice < kNslice; iSlice++) { + if (nCluster[iPlane][iSlice]) { + fdEdxPlane[iPlane][iSlice] /= nCluster[iPlane][iSlice]; + } + } + } + +} //_____________________________________________________________________________ -Int_t AliTRDtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm) +void AliTRDtrack::CookdEdxNN(Float_t *dedx) { - // Propagates a track of particle with mass=pm to a reference plane - // defined by x=xk through media of density=rho and radiationLength=x0 + // + // This function calcuates the input for the neural networks + // which are used for the PID. + // + + //number of time bins in chamber + Int_t ntb = AliTRDcalibDB::Instance()->GetNumberOfTimeBins(); + + Int_t plane; // plane of current cluster + Int_t tb; // time bin of current cluster + Int_t slice; // curent slice + AliTRDcluster *cluster = 0x0; // pointer to current cluster + const Int_t kMLPscale = 16000; // scaling of the MLP input to be smaller than 1 + + // Reset class and local contors/variables + for (Int_t iPlane = 0; iPlane < kNplane; iPlane++){ + for (Int_t iSlice = 0; iSlice < kNMLPslice; iSlice++) { + *(dedx + (kNMLPslice * iPlane) + iSlice) = 0.0; + } + } - if (TMath::Abs(fC*xk - fE) >= 0.99999) { - if (fN>4) cerr<GetCluster(fIndex[iClus]); + if (!cluster) { + continue; + } + + // Read info from current cluster + plane = AliTRDgeometry::GetLayer(cluster->GetDetector()); + if (plane < 0 || plane >= kNplane) { + AliError(Form("Wrong plane %d",plane)); + continue; + } + + tb = cluster->GetLocalTimeBin(); + if (tb == 0 || tb >= ntb) { + AliWarning(Form("time bin 0 or > %d in cluster %d",ntb,iClus)); + AliInfo(Form("dQ/dl %f",fdQdl[iClus])); + continue; + } + + slice = tb * kNMLPslice / ntb; + + *(dedx+(kNMLPslice * plane) + slice) += fdQdl[iClus]/kMLPscale; + + } // End of loop over cluster + +} + +//_____________________________________________________________________________ +void AliTRDtrack::SetTrackSegmentDirMom(const Int_t plane) +{ + // + // Set the momenta for a track segment in a given plane + // + + if ((plane < 0) || + (plane>= kNplane)) { + AliError(Form("Trying to access out of range plane (%d)", plane)); + return; } + + fSnp[plane] = GetSnp(); + fTgl[plane] = GetTgl(); + Double_t p[3]; + GetPxPyPz(p); + fMom[plane] = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]); - Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1, y1=fY, z1=fZ; - Double_t c1=fC*x1 - fE, r1=sqrt(1.- c1*c1); - Double_t c2=fC*x2 - fE, r2=sqrt(1.- c2*c2); +} - fY += dx*(c1+c2)/(r1+r2); - fZ += dx*(c1+c2)/(c1*r2 + c2*r1)*fT; +//_____________________________________________________________________________ +Float_t AliTRDtrack::GetTrackLengthPlane(Int_t plane) const +{ + // + // Calculate the track length of a track segment in a given plane + // - //f = F - 1 - Double_t rr=r1+r2, cc=c1+c2, xx=x1+x2; - Double_t f02= dx*(rr*xx + cc*(c1*x1/r1+c2*x2/r2))/(rr*rr); - Double_t f03=-dx*(2*rr + cc*(c1/r1 + c2/r2))/(rr*rr); - Double_t cr=c1*r2+c2*r1; - Double_t f12= dx*fT*(cr*xx-cc*(r1*x2-c2*c1*x1/r1+r2*x1-c1*c2*x2/r2))/(cr*cr); - Double_t f13=-dx*fT*(2*cr + cc*(c2*c1/r1-r1 + c1*c2/r2-r2))/(cr*cr); - Double_t f14= dx*cc/cr; + if ((plane < 0) || + (plane >= kNplane)) { + return 0.0; + } - //b = C*ft - Double_t b00=f02*fCcy + f03*fCey, b01=f12*fCcy + f13*fCey + f14*fCty; - Double_t b10=f02*fCcz + f03*fCez, b11=f12*fCcz + f13*fCez + f14*fCtz; - Double_t b20=f02*fCcc + f03*fCec, b21=f12*fCcc + f13*fCec + f14*fCtc; - Double_t b30=f02*fCec + f03*fCee, b31=f12*fCec + f13*fCee + f14*fCte; - Double_t b40=f02*fCtc + f03*fCte, b41=f12*fCtc + f13*fCte + f14*fCtt; + return (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick()) + / TMath::Sqrt((1.0 - fSnp[plane]*fSnp[plane]) + / (1.0 + fTgl[plane]*fTgl[plane])); - //a = f*b = f*C*ft - Double_t a00=f02*b20+f03*b30,a01=f02*b21+f03*b31,a11=f12*b21+f13*b31+f14*b41; +} - //F*C*Ft = C + (a + b + bt) - fCyy += a00 + 2*b00; - fCzy += a01 + b01 + b10; - fCcy += b20; - fCey += b30; - fCty += b40; - fCzz += a11 + 2*b11; - fCcz += b21; - fCez += b31; - fCtz += b41; +//_____________________________________________________________________________ +Bool_t AliTRDtrack::CookPID(Int_t &pidQuality) +{ + // + // This function calculates the PID probabilities based on TRD signals + // + // The method produces probabilities based on the charge + // and the position of the maximum time bin in each layer. + // The dE/dx information can be used as global charge or 2 to 3 + // slices. Check AliTRDCalPID and AliTRDCalPIDRefMaker for the actual + // implementation. + // + // Author + // Alex Bercuci (A.Bercuci@gsi.de) 2nd May 2007 + // - fX=x2; + AliTRDcalibDB *calibration = AliTRDcalibDB::Instance(); + if (!calibration) { + AliError("No access to calibration data"); + return kFALSE; + } + + // Retrieve the CDB container class with the probability distributions + const AliTRDCalPID *pd = calibration->GetPIDObject(AliTRDReconstructor::kLQPID); + if (!pd) { + AliError("No access to AliTRDCalPID"); + return kFALSE; + } + // Calculate the input for the NN if fPIDmethod is kNN + Float_t ldEdxNN[AliTRDgeometry::kNlayer * kNMLPslice], *dedx = 0x0; + if(fPIDmethod == kNN) { + CookdEdxNN(&ldEdxNN[0]); + } - //Multiple scattering ****************** + // Sets the a priori probabilities + for(int ispec=0; ispecGetProbability(iSpecies, fMom[iPlane], dedx, length, iPlane); + } - Double_t theta2=14.1*14.1/(beta2*p2*1e6)*d/x0*rho; - fCcc += xz*xz*theta2; - fCec += xz*ez*xy*theta2; - fCtc += xz*zz1*theta2; - fCee += (2*ey*ez*ez*fE+1-ey*ey+ez*ez+fE*fE*ez*ez)*theta2; - fCte += ez*zz1*xy*theta2; - fCtt += zz1*zz1*theta2; + } + if (pidQuality == 0) { + return kTRUE; + } - //Energy losses************************ + // normalize probabilities + Double_t probTotal = 0.0; + for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) { + probTotal += fPID[iSpecies]; + } - Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d*rho; - if (x1 < x2) dE=-dE; - fC*=(1.- sqrt(p2+pm*pm)/p2*dE); - //fE*=(1.- sqrt(p2+pm*pm)/p2*dE); + if (probTotal <= 0.0) { + AliWarning("The total probability over all species <= 0."); + AliWarning("This may be caused by some error in the reference histograms."); + return kFALSE; + } - return 1; + for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) { + fPID[iSpecies] /= probTotal; + } -} + return kTRUE; +} //_____________________________________________________________________________ -void AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, UInt_t index) +Bool_t AliTRDtrack::PropagateTo(Double_t xk, Double_t xx0, Double_t xrho) { - // Assignes found cluster to the track and updates track information - - Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2(); - r00+=fCyy; r01+=fCzy; r11+=fCzz; - Double_t det=r00*r11 - r01*r01; - Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det; - - Double_t k00=fCyy*r00+fCzy*r01, k01=fCyy*r01+fCzy*r11; - Double_t k10=fCzy*r00+fCzz*r01, k11=fCzy*r01+fCzz*r11; - Double_t k20=fCcy*r00+fCcz*r01, k21=fCcy*r01+fCcz*r11; - Double_t k30=fCey*r00+fCez*r01, k31=fCey*r01+fCez*r11; - Double_t k40=fCty*r00+fCtz*r01, k41=fCty*r01+fCtz*r11; - - Double_t dy=c->GetY() - fY, dz=c->GetZ() - fZ; - Double_t cur=fC + k20*dy + k21*dz, eta=fE + k30*dy + k31*dz; - if (TMath::Abs(cur*fX-eta) >= 0.99999) { - if (fN>4) cerr< 0.0001) { + // Make correction for curvature if neccesary + l2 = 0.5 * TMath::Sqrt((x-oldX)*(x-oldX) + + (y-oldY)*(y-oldY)); + l2 = 2.0 * TMath::ASin(l2 * crv) / crv; + l2 = TMath::Sqrt(l2*l2 + (z-oldZ)*(z-oldZ)); + } + AddTimeStep(l2); + } + } + + if (!AliExternalTrackParam::CorrectForMeanMaterial(xx0,xrho,GetMass())) { + return kFALSE; + } + + { + + // Energy losses + Double_t p2 = (1.0 + GetTgl()*GetTgl()) / (GetSigned1Pt()*GetSigned1Pt()); + Double_t beta2 = p2 / (p2 + GetMass()*GetMass()); + if ((beta2 < 1.0e-10) || + ((5940.0 * beta2/(1.0 - beta2 + 1.0e-10) - beta2) < 0.0)) { + return kFALSE; + } + + Double_t dE = 0.153e-3 / beta2 + * (TMath::Log(5940.0 * beta2/(1.0 - beta2 + 1.0e-10)) - beta2) + * xrho; + fBudget[0] += xrho; + + /* + // Suspicious part - think about it ? + Double_t kinE = TMath::Sqrt(p2); + if (dE > 0.8*kinE) dE = 0.8 * kinE; // + if (dE < 0) dE = 0.0; // Not valid region for Bethe bloch + */ + + fDE += dE; + + /* + // Suspicious ! I.B. + Double_t sigmade = 0.07 * TMath::Sqrt(TMath::Abs(dE)); // Energy loss fluctuation + Double_t sigmac2 = sigmade*sigmade*fC*fC*(p2+GetMass()*GetMass())/(p2*p2); + fCcc += sigmac2; + fCee += fX*fX * sigmac2; + */ + } - fY += k00*dy + k01*dz; - fZ += k10*dy + k11*dz; - fC = cur; - fE = eta; - fT += k40*dy + k41*dz; + return kTRUE; + +} + +//_____________________________________________________________________________ +Bool_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq + , Int_t index, Double_t h01) +{ + // + // Assignes the found cluster to the track and updates + // track information. + // : predicted chi2 + // : ?? + // : Tilting factor + // + + Double_t p[2] = { c->GetY() + , c->GetZ() }; + Double_t sy2 = c->GetSigmaY2() * 4.0; + Double_t sz2 = c->GetSigmaZ2() * 4.0; + Double_t cov[3] = { sy2 + h01*h01*sz2 + , h01 * (sy2-sz2) + , sz2 + h01*h01*sy2 }; + + if (!AliExternalTrackParam::Update(p,cov)) { + return kFALSE; + } + + Int_t n = GetNumberOfClusters(); + fIndex[n] = index; + SetNumberOfClusters(n+1); + + SetChi2(GetChi2()+chisq); + + return kTRUE; + +} + +//_____________________________________________________________________________ +Int_t AliTRDtrack::UpdateMI(AliTRDcluster *c, Double_t chisq, Int_t index + , Double_t h01, Int_t /*plane*/, Int_t /*tid*/) +{ + // + // Assignes the found cluster to the track and + // updates track information + // : predicted chi2 + // : ?? + // : Tilting factor + // + // Difference to Update(AliTRDcluster *c): cluster is added to fClusters + // - Double_t c01=fCzy, c02=fCcy, c03=fCey, c04=fCty; - Double_t c12=fCcz, c13=fCez, c14=fCtz; + Double_t p[2] = { c->GetY() + , c->GetZ() }; + Double_t sy2 = c->GetSigmaY2() * 4.0; + Double_t sz2 = c->GetSigmaZ2() * 4.0; + Double_t cov[3] = { sy2 + h01*h01*sz2 + , h01 * (sy2-sz2) + , sz2 + h01*h01*sy2 }; - fCyy-=k00*fCyy+k01*fCzy; fCzy-=k00*c01+k01*fCzz; - fCcy-=k00*c02+k01*c12; fCey-=k00*c03+k01*c13; - fCty-=k00*c04+k01*c14; + if (!AliExternalTrackParam::Update(p,cov)) { + return kFALSE; + } + + AliTracker::FillResiduals(this,p,cov,c->GetVolumeId()); - fCzz-=k10*c01+k11*fCzz; - fCcz-=k10*c02+k11*c12; fCez-=k10*c03+k11*c13; - fCtz-=k10*c04+k11*c14; + // Register cluster to track + Int_t n = GetNumberOfClusters(); + fIndex[n] = index; + fClusters[n] = c; + SetNumberOfClusters(n+1); + SetChi2(GetChi2() + chisq); - fCcc-=k20*c02+k21*c12; fCec-=k20*c03+k21*c13; - fCtc-=k20*c04+k21*c14; + return kTRUE; - fCee-=k30*c03+k31*c13; - fCte-=k30*c04+k31*c14; +} - fCtt-=k40*c04+k41*c14; +//_____________________________________________________________________________ +Bool_t AliTRDtrack::Update(const AliTRDtracklet &t, Double_t chisq, Int_t index) +{ + // + // Assignes the found tracklet to the track + // and updates track information + // : predicted chi2 + // : ?? + // - fIndex[fN]=index; + Double_t h01 = t.GetTilt(); // Is this correct???? + Double_t p[2] = { t.GetY(), t.GetZ() }; + Double_t sy2 = t.GetTrackletSigma2(); // Error pad-column + Double_t sz2 = 100000.0; // Error pad-row (????) + Double_t cov[3] = { sy2 + h01*h01*sz2 // Does this have any sense???? + , h01 * (sy2 - sz2) + , sz2 + h01*h01*sy2 }; + if (!AliExternalTrackParam::Update(p,cov)) { + return kFALSE; + } - Float_t q = c->GetQ(); - Double_t s = fX*fC - fE, t=fT; - q *= TMath::Sqrt((1-s*s)/(1+t*t)); - fdQdl[fN++] = q; + Int_t n = GetNumberOfClusters(); + fIndex[n] = index; + SetNumberOfClusters(n+1); + SetChi2(GetChi2()+chisq); - fChi2 += chisq; + return kTRUE; - // cerr<<"in update: fIndex["<= 0.99999) { - if (fN>4) cerr<GetY() + , c->GetZ() }; + Double_t sy2 = c->GetSigmaY2() * 4.0; + Double_t sz2 = c->GetSigmaZ2() * 4.0; + Double_t cov[3] = { sy2 + h01*h01*sz2 + , h01*(sy2-sz2) + , sz2 + h01*h01*sy2 }; + + return AliExternalTrackParam::GetPredictedChi2(p,cov); + +} + +//_____________________________________________________________________________ +void AliTRDtrack::MakeBackupTrack() +{ + // + // Creates a backup track + // + + if (fBackupTrack) { + delete fBackupTrack; } + fBackupTrack = new AliTRDtrack(*this); + +} + +//_____________________________________________________________________________ +Int_t AliTRDtrack::GetProlongation(Double_t xk, Double_t &y, Double_t &z) +{ + // + // Find a prolongation at given x + // Return 0 if it does not exist + // + + Double_t bz = GetBz(); - Double_t y0=fY + sqrt(1.- r2*r2)/fC; - if ((fY-y0)*fC >= 0.) { - if (fN>4) cerr< xr) ? -1.0 : 1.0; + + // Direction +- + for (Double_t x = GetX()+dir*step; dir*x < dir*xr; x += dir*step) { + + GetXYZ(xyz0); + GetProlongation(x,y,z); + xyz1[0] = x * TMath::Cos(GetAlpha()) + y * TMath::Sin(GetAlpha()); + xyz1[1] = x * TMath::Sin(GetAlpha()) - y * TMath::Cos(GetAlpha()); + xyz1[2] = z; + Double_t param[7]; + AliTracker::MeanMaterialBudget(xyz0,xyz1,param); + + if ((param[0] > 0) && + (param[1] > 0)) { + PropagateTo(x,param[1],param[0]*param[4]); + } - // *** Double_t dy2=fCyy; - - //F*C*Ft = C + (a + b + bt) - fCyy += a00 + 2*b00; - fCzy += b10; - fCcy += b20; - fCey += a03+b30+b03; - fCty += b40; - fCez += b13; - fCec += b23; - fCee += a33 + 2*b33; - fCte += b43; + if (GetY() > GetX()*kTalphac) { + Rotate(-kAlphac); + } + if (GetY() < -GetX()*kTalphac) { + Rotate( kAlphac); + } + + } - // *** fCyy+=dy2*sa*sa*r1*r1/(1.- r1*r1); - // *** fCzz+=d2y*sa*sa*fT*fT/(1.- r1*r1); + PropagateTo(xr); - return 1; -} + return 0; +} //_____________________________________________________________________________ -Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c) const +Int_t AliTRDtrack::PropagateToR(Double_t r,Double_t step) { - /* - Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2(); - r00+=fCyy; r01+=fCzy; r11+=fCzz; + // + // Propagate track to the radial position + // Rotation always connected to the last track position + // - Double_t det=r00*r11 - r01*r01; - if (TMath::Abs(det) < 1.e-10) { - if (fN>4) cerr< r) ? -1.0 : 1.0; + + for (Double_t x = radius+dir*step; dir*x < dir*r; x += dir*step) { + + GetXYZ(xyz0); + Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]); + Rotate(alpha,kTRUE); + GetXYZ(xyz0); + GetProlongation(x,y,z); + xyz1[0] = x * TMath::Cos(alpha) + y * TMath::Sin(alpha); + xyz1[1] = x * TMath::Sin(alpha) - y * TMath::Cos(alpha); + xyz1[2] = z; + Double_t param[7]; + AliTracker::MeanMaterialBudget(xyz0,xyz1,param); + if (param[1] <= 0) { + param[1] = 100000000; + } + PropagateTo(x,param[1],param[0]*param[4]); + + } + + GetXYZ(xyz0); + Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]); + Rotate(alpha,kTRUE); + GetXYZ(xyz0); + GetProlongation(r,y,z); + xyz1[0] = r * TMath::Cos(alpha) + y * TMath::Sin(alpha); + xyz1[1] = r * TMath::Sin(alpha) - y * TMath::Cos(alpha); + xyz1[2] = z; + Double_t param[7]; + AliTracker::MeanMaterialBudget(xyz0,xyz1,param); + + if (param[1] <= 0) { + param[1] = 100000000; } - Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01; + PropagateTo(r,param[1],param[0]*param[4]); + + return 0; + +} - Double_t dy=c->GetY() - fY, dz=c->GetZ() - fZ; +//_____________________________________________________________________________ +Int_t AliTRDtrack::GetSector() const +{ + // + // Return the current sector + // - return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det; - */ + return Int_t(TVector2::Phi_0_2pi(GetAlpha()) / AliTRDgeometry::GetAlpha()) + % AliTRDgeometry::kNsector; - Double_t dy=c->GetY() - fY; - Double_t r00=c->GetSigmaY2(); +} - return (dy*dy)/r00; +//_____________________________________________________________________________ +void AliTRDtrack::SetSampledEdx(Float_t q, Int_t i) +{ + // + // The sampled energy loss + // -} + Double_t s = GetSnp(); + Double_t t = GetTgl(); + q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t)); + fdQdl[i] = q; +} -//_________________________________________________________________________ -void AliTRDtrack::GetPxPyPz(Double_t& px, Double_t& py, Double_t& pz) const +//_____________________________________________________________________________ +void AliTRDtrack::SetSampledEdx(Float_t q) { - // Returns reconstructed track momentum in the global system. + // + // The sampled energy loss + // - Double_t pt=TMath::Abs(GetPt()); // GeV/c - Double_t r=fC*fX-fE; - Double_t y0=fY + sqrt(1.- r*r)/fC; - px=-pt*(fY-y0)*fC; //cos(phi); - py=-pt*(fE-fX*fC); //sin(phi); - pz=pt*fT; - Double_t tmp=px*TMath::Cos(fAlpha) - py*TMath::Sin(fAlpha); - py=px*TMath::Sin(fAlpha) + py*TMath::Cos(fAlpha); - px=tmp; + Double_t s = GetSnp(); + Double_t t = GetTgl(); + q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t)); + fdQdl[fNdedx] = q; + fNdedx++; -} +} -//____________________________________________________________________________ -void AliTRDtrack::Streamer(TBuffer &R__b) +//_____________________________________________________________________________ +Double_t AliTRDtrack::GetBz() const { - if (R__b.IsReading()) { - Version_t R__v = R__b.ReadVersion(); if (R__v) { } - TObject::Streamer(R__b); - R__b >> fLab; - R__b >> fChi2; - R__b >> fdEdx; - R__b >> fAlpha; - R__b >> fX; - R__b >> fY; - R__b >> fZ; - R__b >> fC; - R__b >> fE; - R__b >> fT; - R__b >> fCyy; - R__b >> fCzy; - R__b >> fCzz; - R__b >> fCcy; - R__b >> fCcz; - R__b >> fCcc; - R__b >> fCey; - R__b >> fCez; - R__b >> fCec; - R__b >> fCee; - R__b >> fCty; - R__b >> fCtz; - R__b >> fCtc; - R__b >> fCte; - R__b >> fCtt; - R__b >> fN; - for (Int_t i=0; i> fIndex[i]; - for (Int_t i=0; i> fdQdl[i]; - } else { - R__b.WriteVersion(AliTRDtrack::IsA()); - TObject::Streamer(R__b); - R__b << fLab; - R__b << fChi2; - R__b << fdEdx; - R__b << fAlpha; - R__b << fX; - R__b << fY; - R__b << fZ; - R__b << fC; - R__b << fE; - R__b << fT; - R__b << fCyy; - R__b << fCzy; - R__b << fCzz; - R__b << fCcy; - R__b << fCcz; - R__b << fCcc; - R__b << fCey; - R__b << fCez; - R__b << fCec; - R__b << fCee; - R__b << fCty; - R__b << fCtz; - R__b << fCtc; - R__b << fCte; - R__b << fCtt; - R__b << fN; - for (Int_t i=0; i