X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TRD%2FAliTRDtrack.cxx;h=7d0bfa1a69587a0cee4e2ca7db3160b9f4f76839;hb=e8974bc9b30bf13ace888b2405dc360a75b5a6bc;hp=be6a363045486f560b2f9e9d33a2b007207757dc;hpb=6c94f3308e8e21efdd662a43a6a60db688334d59;p=u%2Fmrichter%2FAliRoot.git diff --git a/TRD/AliTRDtrack.cxx b/TRD/AliTRDtrack.cxx index be6a3630454..7d0bfa1a695 100644 --- a/TRD/AliTRDtrack.cxx +++ b/TRD/AliTRDtrack.cxx @@ -15,16 +15,15 @@ /* $Id$ */ -#include -#include #include #include "AliTracker.h" -#include "AliESDtrack.h" + #include "AliTRDgeometry.h" #include "AliTRDcluster.h" #include "AliTRDtrack.h" -#include "AliTRDtracklet.h" +#include "AliTRDcalibDB.h" +#include "Cal/AliTRDCalPID.h" ClassImp(AliTRDtrack) @@ -35,258 +34,348 @@ ClassImp(AliTRDtrack) // // /////////////////////////////////////////////////////////////////////////////// -AliTRDtrack::AliTRDtrack(): - AliKalmanTrack(), - fSeedLab(-1), - fdEdx(0), - fdEdxT(0), - fDE(0), - fStopped(kFALSE), - fLhElectron(0), - fNWrong(0), - fNRotate(0), - fNCross(0), - fNExpected(0), - fNLast(0), - fNExpectedLast(0), - fNdedx(0), - fChi2Last(1e10), - fBackupTrack(0x0) +//_____________________________________________________________________________ +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) { - for (Int_t i=0; iGetQ()); - Double_t s = GetSnp(), t=GetTgl(); - if(s*s < 1) q *= TMath::Sqrt((1-s*s)/(1+t*t)); + 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; + } - fdQdl[0] = q; - // initialisation [SR, GSI 18.02.2003] (i startd for 1) - for(UInt_t i=1; iGetX() + ,par->GetAlpha() + ,par->GetParameter() + ,par->GetCovariance()); + + for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) { + fdQdl[i] = 0; + fClusters[i] = 0x0; } - Set(par->GetX(),par->GetAlpha(),par->GetParameter(),par->GetCovariance()); - - for (UInt_t i=0; i110&&fChi2/(Float_t(fN))<3) return 3; //gold - if (fNLast>30&&fChi2Last/(Float_t(fNLast))<3) return 3; //gold - if (fNLast>20&&fChi2Last/(Float_t(fNLast))<2) return 3; //gold - if (fNLast/(fNExpectedLast+3.)>0.8 && fChi2Last/Float_t(fNLast)<5&&fNLast>20) return 2; //silber - if (fNLast>5 &&((fNLast+1.)/(fNExpectedLast+1.))>0.8&&fChi2Last/(fNLast-5.)<6) return 1; - - return status; - } //_____________________________________________________________________________ @@ -334,8 +426,6 @@ Int_t AliTRDtrack::Compare(const TObject *o) const // AliTRDtrack *t = (AliTRDtrack *) o; - // Double_t co=t->GetSigmaY2(); - // Double_t c =GetSigmaY2(); Double_t co = TMath::Abs(t->GetC()); Double_t c = TMath::Abs(GetC()); @@ -346,6 +436,7 @@ Int_t AliTRDtrack::Compare(const TObject *o) const else if (c < co) { return -1; } + return 0; } @@ -357,314 +448,508 @@ void AliTRDtrack::CookdEdx(Double_t low, Double_t up) // Calculates the truncated dE/dx within the "low" and "up" cuts. // - Int_t i = 0; - // Array to sort the dEdx values according to amplitude Float_t sorted[kMAXCLUSTERSPERTRACK]; - - // Number of clusters used for dedx - Int_t nc = fNdedx; - + fdEdx = 0.0; + // Require at least 10 clusters for a dedx measurement - if (nc < 10) { - SetdEdx(0); + if (fNdedx < 10) { return; } - // Lower and upper bound - Int_t nl = Int_t(low * nc); - Int_t nu = Int_t( up * nc); - // Can fdQdl be negative ???? - for (i = 0; i < nc; i++) { + for (Int_t i = 0; i < fNdedx; i++) { sorted[i] = TMath::Abs(fdQdl[i]); } - // Sort the dedx values by amplitude - Int_t *index = new Int_t[nc]; - TMath::Sort(nc,sorted,index,kFALSE); - - // Sum up the truncated charge between nl and nu - Float_t dedx = 0.0; - for (i = nl; i <= nu; i++) { - dedx += sorted[index[i]]; + Int_t *index = new Int_t[fNdedx]; + TMath::Sort(fNdedx, sorted, index, kFALSE); + + // Sum up the truncated charge between lower and upper bounds + Int_t nl = Int_t(low * fNdedx); + Int_t nu = Int_t( up * fNdedx); + for (Int_t i = nl; i <= nu; i++) { + fdEdx += sorted[index[i]]; } - dedx /= (nu - nl + 1.0); - SetdEdx(dedx); + fdEdx /= (nu - nl + 1.0); + + delete[] index; } //_____________________________________________________________________________ -Bool_t AliTRDtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho) +void AliTRDtrack::CookdEdxTimBin(const Int_t/* tid*/) { - // Propagates a track of particle with mass=pm to a reference plane - // defined by x=xk through media of density=rho and radiationLength=x0 - - if (xk == GetX()) return kTRUE; - - Double_t oldX=GetX(), oldY=GetY(), oldZ=GetZ(); - - Double_t bz=GetBz(); - if (!AliExternalTrackParam::PropagateTo(xk,bz)) return kFALSE; - - Double_t x=GetX(), y=GetY(), z=GetZ(); - Double_t d=TMath::Sqrt((x-oldX)*(x-oldX)+(y-oldY)*(y-oldY)+(z-oldZ)*(z-oldZ)); - if (oldX < xk) - if (IsStartedTimeIntegral()) { - Double_t l2=d; - Double_t crv=GetC(); - if (TMath::Abs(l2*crv)>0.0001){ - // make correction for curvature if neccesary - l2 = 0.5*TMath::Sqrt((x-oldX)*(x-oldX) + (y-oldY)*(y-oldY)); - l2 = 2*TMath::ASin(l2*crv)/crv; - l2 = TMath::Sqrt(l2*l2+(z-oldZ)*(z-oldZ)); + // + // Set fdEdxPlane and fTimBinPlane and also get the + // Time bin for Max. Cluster + // + // Authors: + // Prashant Shukla (shukla@physi.uni-heidelberg.de) + // Alexandru Bercuci (A.Bercuci@gsi.de) + // + + // Max charge in chamber + Double_t maxcharge[kNplane]; + // Number of clusters attached to track per chamber and slice + Int_t nCluster[kNplane][kNslice]; + // 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; // 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; } - AddTimeStep(l2); - } - - Double_t ll = (oldX < xk) ? -d : d; - if (!AliExternalTrackParam::CorrectForMaterial(ll*rho/x0,x0,GetMass())) - return kFALSE; - - {//Energy losses************************ - Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()); - Double_t beta2=p2/(p2 + GetMass()*GetMass()); - if ((5940*beta2/(1-beta2+1e-10) - beta2) < 0) return kFALSE; - - Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2+1e-10)) - beta2)*d*rho; - Float_t budget = d*rho; - fBudget[0] += budget; - /* - // 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; - */ } - return kTRUE; + // 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]; + } + } + } + +} //_____________________________________________________________________________ -Bool_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, Int_t index, - Double_t h01) +void AliTRDtrack::CookdEdxNN(Float_t *dedx) { - // Assignes found cluster to the track and updates track information - - Bool_t fNoTilt = kTRUE; - if(TMath::Abs(h01) > 0.003) fNoTilt = kFALSE; - // add angular effect to the error contribution - MI - Float_t tangent2 = GetSnp()*GetSnp(); - if (tangent2 < 0.90000){ - tangent2 = tangent2/(1.-tangent2); + // + // 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; + } } - //Float_t errang = tangent2*0.04; // - Double_t p[2]={c->GetY(), c->GetZ()}; - //Double_t cov[3]={c->GetSigmaY2()+errang, 0., c->GetSigmaZ2()*100.}; - Double_t sy2=c->GetSigmaY2()*4; - Double_t sz2=c->GetSigmaZ2()*4; - Double_t cov[3]={sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2}; + // Start looping over clusters attached to this track + for (Int_t iClus = 0; iClus < GetNumberOfClusters(); iClus++) { - if (!AliExternalTrackParam::Update(p,cov)) return kFALSE; + 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; + } - Int_t n=GetNumberOfClusters(); - fIndex[n]=index; - SetNumberOfClusters(n+1); + 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; + } - SetChi2(GetChi2()+chisq); + slice = tb * kNMLPslice / ntb; + + *(dedx+(kNMLPslice * plane) + slice) += fdQdl[iClus]/kMLPscale; + + } // End of loop over cluster - return kTRUE; -} +} //_____________________________________________________________________________ -Int_t AliTRDtrack::UpdateMI(const AliTRDcluster *c, Double_t chisq, Int_t index, Double_t h01, Int_t /*plane*/) { - // Assignes found cluster to the track and updates track information - Bool_t fNoTilt = kTRUE; - if(TMath::Abs(h01) > 0.003) fNoTilt = kFALSE; - // add angular effect to the error contribution and make correction - MI - // - Double_t tangent2 = GetSnp()*GetSnp(); - if (tangent2 < 0.90000){ - tangent2 = tangent2/(1.-tangent2); - } - Double_t tangent = TMath::Sqrt(tangent2); - if (GetSnp()<0) tangent*=-1; - // Double_t correction = 0*plane; - /* - Double_t errang = tangent2*0.04; // - Double_t errsys =0.025*0.025*20; //systematic error part - - Float_t extend =1; - if (c->GetNPads()==4) extend=2; - */ - //if (c->GetNPads()==5) extend=3; - //if (c->GetNPads()==6) extend=3; - //if (c->GetQ()<15) return 1; - - /* - if (corrector!=0){ - //if (0){ - correction = corrector->GetCorrection(plane,c->GetLocalTimeBin(),tangent); - if (TMath::Abs(correction)>0){ - //if we have info - errang = corrector->GetSigma(plane,c->GetLocalTimeBin(),tangent); - errang *= errang; - errang += tangent2*0.04; - } +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]); + +} + +//_____________________________________________________________________________ +Float_t AliTRDtrack::GetTrackLengthPlane(Int_t plane) const +{ + // + // Calculate the track length of a track segment in a given plane // - //Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12.); - /* - { - Double_t dy=c->GetY() - GetY(), dz=c->GetZ() - GetZ(); - printf("%e %e %e %e\n",dy,dz,padlength/2,h01); - } - */ - Double_t p[2]={c->GetY(), c->GetZ()}; - /* - Double_t cov[3]={(c->GetSigmaY2()+errang+errsys)*extend, 0., - c->GetSigmaZ2()*10000.}; - */ - Double_t sy2=c->GetSigmaY2()*4; - Double_t sz2=c->GetSigmaZ2()*4; - 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; -} -/* + if ((plane < 0) || + (plane >= kNplane)) { + return 0.0; + } + + return (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick()) + / TMath::Sqrt((1.0 - fSnp[plane]*fSnp[plane]) + / (1.0 + fTgl[plane]*fTgl[plane])); + +} + //_____________________________________________________________________________ -Int_t AliTRDtrack::UpdateMI(const AliTRDtracklet &tracklet) +Bool_t AliTRDtrack::CookPID(Int_t &pidQuality) { // - // Assignes found tracklet to the track and updates track information + // 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. // - Double_t r00=(tracklet.GetTrackletSigma2()), r01=0., r11= 10000.; - r00+=fCyy; r01+=fCzy; r11+=fCzz; + // Author + // Alex Bercuci (A.Bercuci@gsi.de) 2nd May 2007 + // + + 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(fPIDmethod == kNN ? AliTRDrecoParam::kNNPID : AliTRDrecoParam::kLQPID); + if (!pd) { + AliError("No access to AliTRDCalPID"); + return kFALSE; + } + + // Calculate the input for the NN if fPIDmethod is kNN + Float_t ldEdxNN[AliTRDCalPID::kNPlane * kNMLPslice], *dedx = 0x0; + if(fPIDmethod == kNN) { + CookdEdxNN(&ldEdxNN[0]); + } + + // Sets the a priori probabilities + for(int ispec=0; ispecGetProbability(iSpecies, fMom[iPlane], dedx, length, iPlane); + } + + } + + if (pidQuality == 0) { + return kTRUE; + } + + // normalize probabilities + Double_t probTotal = 0.0; + for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) { + probTotal += fPID[iSpecies]; + } + + 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; + } + + for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) { + fPID[iSpecies] /= probTotal; + } + + return kTRUE; + +} + +//_____________________________________________________________________________ +Bool_t AliTRDtrack::PropagateTo(Double_t xk, Double_t xx0, Double_t xrho) +{ // - Double_t det=r00*r11 - r01*r01; - Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det; + // Propagates this track to a reference plane defined by "xk" [cm] + // correcting for the mean crossed material. // + // "xx0" - thickness/rad.length [units of the radiation length] + // "xrho" - thickness*density [g/cm^2] + // - Double_t dy=tracklet.GetY() - fY, dz=tracklet.GetZ() - fZ; + if (xk == GetX()) { + return kTRUE; + } - - Double_t s00 = tracklet.GetTrackletSigma2(); // error pad - Double_t s11 = 100000; // error pad-row - Float_t h01 = tracklet.GetTilt(); - // - // r00 = fCyy + 2*fCzy*h01 + fCzz*h01*h01+s00; - r00 = fCyy + fCzz*h01*h01+s00; - // r01 = fCzy + fCzz*h01; - r01 = fCzy ; - r11 = fCzz + s11; - det = r00*r11 - r01*r01; - // inverse matrix - 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=fCey*r00+fCez*r01, k21=fCey*r01+fCez*r11; - Double_t k30=fCty*r00+fCtz*r01, k31=fCty*r01+fCtz*r11; - Double_t k40=fCcy*r00+fCcz*r01, k41=fCcy*r01+fCcz*r11; - - // K matrix -// k00=fCyy*r00+fCzy*(r01+h01*r00),k01=fCyy*r01+fCzy*(r11+h01*r01); -// k10=fCzy*r00+fCzz*(r01+h01*r00),k11=fCzy*r01+fCzz*(r11+h01*r01); -// k20=fCey*r00+fCez*(r01+h01*r00),k21=fCey*r01+fCez*(r11+h01*r01); -// k30=fCty*r00+fCtz*(r01+h01*r00),k31=fCty*r01+fCtz*(r11+h01*r01); -// k40=fCcy*r00+fCcz*(r01+h01*r00),k41=fCcy*r01+fCcz*(r11+h01*r01); - // - //Update measurement - Double_t cur=fC + k40*dy + k41*dz, eta=fE + k20*dy + k21*dz; - // cur=fC + k40*dy + k41*dz; eta=fE + k20*dy + k21*dz; - if (TMath::Abs(cur*fX-eta) >= 0.90000) { - //Int_t n=GetNumberOfClusters(); - // if (n>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; + */ + + } + + 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 oldyy = fCyy, oldzz = fCzz; //, oldee=fCee, oldcc =fCcc; - Double_t oldzy = fCzy, oldey = fCey, oldty=fCty, oldcy =fCcy; - Double_t oldez = fCez, oldtz = fCtz, oldcz=fCcz; - //Double_t oldte = fCte, oldce = fCce; - //Double_t oldct = fCct; - fCyy-=k00*oldyy+k01*oldzy; - fCzy-=k10*oldyy+k11*oldzy; - fCey-=k20*oldyy+k21*oldzy; - fCty-=k30*oldyy+k31*oldzy; - fCcy-=k40*oldyy+k41*oldzy; + 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*/) +{ // - fCzz-=k10*oldzy+k11*oldzz; - fCez-=k20*oldzy+k21*oldzz; - fCtz-=k30*oldzy+k31*oldzz; - fCcz-=k40*oldzy+k41*oldzz; + // Assignes the found cluster to the track and + // updates track information + // : predicted chi2 + // : ?? + // : Tilting factor // - fCee-=k20*oldey+k21*oldez; - fCte-=k30*oldey+k31*oldez; - fCce-=k40*oldey+k41*oldez; + // Difference to Update(AliTRDcluster *c): cluster is added to fClusters // - fCtt-=k30*oldty+k31*oldtz; - fCct-=k40*oldty+k41*oldtz; + + 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; + } + + AliTracker::FillResiduals(this,p,cov,c->GetVolumeId()); + + // Register cluster to track + Int_t n = GetNumberOfClusters(); + fIndex[n] = index; + fClusters[n] = c; + SetNumberOfClusters(n+1); + SetChi2(GetChi2() + chisq); + + return kTRUE; + +} + +//_____________________________________________________________________________ +Bool_t AliTRDtrack::Update(const AliTRDtracklet &t, Double_t chisq, Int_t index) +{ // - fCcc-=k40*oldcy+k41*oldcz; + // Assignes the found tracklet to the track + // and updates track information + // : predicted chi2 + // : ?? // - - //Int_t n=GetNumberOfClusters(); - //fIndex[n]=index; - //SetNumberOfClusters(n+1); - //SetChi2(GetChi2()+chisq); - // cerr<<"in update: fIndex["<GetY(), c->GetZ()}; - Double_t sy2=c->GetSigmaY2()*4; - Double_t sz2=c->GetSigmaZ2()*4; - Double_t cov[3]={sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2}; + 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 }; return AliExternalTrackParam::GetPredictedChi2(p,cov); - /* - Bool_t fNoTilt = kTRUE; - if(TMath::Abs(h01) > 0.003) fNoTilt = kFALSE; - - return (c->GetY() - GetY())*(c->GetY() - GetY())/c->GetSigmaY2(); - */ - - /* - Double_t chi2, dy, r00, r01, r11; - - if(fNoTilt) { - dy=c->GetY() - fY; - r00=c->GetSigmaY2(); - chi2 = (dy*dy)/r00; - } - else { - Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12); - // - r00=c->GetSigmaY2(); r01=0.; r11=c->GetSigmaZ2(); - r00+=fCyy; r01+=fCzy; r11+=fCzz; - - Double_t det=r00*r11 - r01*r01; - if (TMath::Abs(det) < 1.e-10) { - Int_t n=GetNumberOfClusters(); - if (n>4) cerr<GetY() - fY, dz=c->GetZ() - fZ; - Double_t tiltdz = dz; - if (TMath::Abs(tiltdz)>padlength/2.) { - tiltdz = TMath::Sign(padlength/2,dz); - } - // dy=dy+h01*dz; - dy=dy+h01*tiltdz; - - chi2 = (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det; - } - - return chi2; - */ -} +} //_____________________________________________________________________________ void AliTRDtrack::MakeBackupTrack() @@ -742,7 +987,9 @@ void AliTRDtrack::MakeBackupTrack() // Creates a backup track // - if (fBackupTrack) delete fBackupTrack; + if (fBackupTrack) { + delete fBackupTrack; + } fBackupTrack = new AliTRDtrack(*this); } @@ -751,52 +998,70 @@ void AliTRDtrack::MakeBackupTrack() Int_t AliTRDtrack::GetProlongation(Double_t xk, Double_t &y, Double_t &z) { // - // Find prolongation at given x - // return 0 if not exist - - Double_t bz=GetBz(); + // Find a prolongation at given x + // Return 0 if it does not exist + // - if (!AliExternalTrackParam::GetYAt(xk,bz,y)) return 0; - if (!AliExternalTrackParam::GetZAt(xk,bz,z)) return 0; + Double_t bz = GetBz(); + + if (!AliExternalTrackParam::GetYAt(xk,bz,y)) { + return 0; + } + if (!AliExternalTrackParam::GetZAt(xk,bz,z)) { + return 0; + } return 1; + } //_____________________________________________________________________________ -Int_t AliTRDtrack::PropagateToX(Double_t xr, Double_t step) +Int_t AliTRDtrack::PropagateToX(Double_t xr, Double_t step) { // // Propagate track to given x position - // works inside of the 20 degree segmentation (local cooordinate frame for TRD , TPC, TOF) + // Works inside of the 20 degree segmentation + // (local cooordinate frame for TRD , TPC, TOF) // - // material budget from geo manager + // Material budget from geo manager // - Double_t xyz0[3], xyz1[3],y,z; - const Double_t kAlphac = TMath::Pi()/9.; + + Double_t xyz0[3]; + Double_t xyz1[3]; + Double_t y; + Double_t z; + + const Double_t kAlphac = TMath::Pi()/9.0; const Double_t kTalphac = TMath::Tan(kAlphac*0.5); - // critical alpha - cross sector indication - // - Double_t dir = (GetX()>xr) ? -1.:1.; - // direction +- - for (Double_t x=GetX()+dir*step;dir*x 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[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]; - AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param); - // - if (param[0]>0&¶m[1]>0) PropagateTo(x,param[1],param[0]); - if (GetY()>GetX()*kTalphac){ + AliTracker::MeanMaterialBudget(xyz0,xyz1,param); + + if ((param[0] > 0) && + (param[1] > 0)) { + PropagateTo(x,param[1],param[0]*param[4]); + } + + if (GetY() > GetX()*kTalphac) { Rotate(-kAlphac); } - if (GetY()<-GetX()*kTalphac){ - Rotate(kAlphac); + if (GetY() < -GetX()*kTalphac) { + Rotate( kAlphac); } + } - // + PropagateTo(xr); return 0; @@ -807,40 +1072,53 @@ Int_t AliTRDtrack::PropagateToX(Double_t xr, Double_t step) Int_t AliTRDtrack::PropagateToR(Double_t r,Double_t step) { // - // propagate track to the radial position - // rotation always connected to the last track position + // Propagate track to the radial position + // Rotation always connected to the last track position // - Double_t xyz0[3], xyz1[3],y,z; + + Double_t xyz0[3]; + Double_t xyz1[3]; + Double_t y; + Double_t z; + Double_t radius = TMath::Sqrt(GetX()*GetX() + GetY()*GetY()); - Double_t dir = (radius>r) ? -1.:1.; // direction +- - // - for (Double_t x=radius+dir*step;dir*x 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[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]; - AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param); - if (param[1]<=0) param[1] =100000000; - PropagateTo(x,param[1],param[0]); + 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[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]; - AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param); - // - if (param[1]<=0) param[1] =100000000; - PropagateTo(r,param[1],param[0]); + AliTracker::MeanMaterialBudget(xyz0,xyz1,param); + + if (param[1] <= 0) { + param[1] = 100000000; + } + PropagateTo(r,param[1],param[0]*param[4]); return 0; @@ -853,9 +1131,8 @@ Int_t AliTRDtrack::GetSector() const // Return the current sector // - return Int_t(TVector2::Phi_0_2pi(GetAlpha()) - / AliTRDgeometry::GetAlpha()) - % AliTRDgeometry::kNsect; + return Int_t(TVector2::Phi_0_2pi(GetAlpha()) / AliTRDgeometry::GetAlpha()) + % AliTRDgeometry::kNsector; } @@ -868,12 +1145,12 @@ void AliTRDtrack::SetSampledEdx(Float_t q, Int_t i) Double_t s = GetSnp(); Double_t t = GetTgl(); - q *= TMath::Sqrt((1-s*s)/(1+t*t)); + q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t)); fdQdl[i] = q; } - //_____________________________________________________________________________ +//_____________________________________________________________________________ void AliTRDtrack::SetSampledEdx(Float_t q) { // @@ -882,19 +1159,25 @@ void AliTRDtrack::SetSampledEdx(Float_t q) Double_t s = GetSnp(); Double_t t = GetTgl(); - q*= TMath::Sqrt((1-s*s)/(1+t*t)); + q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t)); fdQdl[fNdedx] = q; fNdedx++; } -Double_t AliTRDtrack::GetBz() const { +//_____________________________________________________________________________ +Double_t AliTRDtrack::GetBz() const +{ // - // returns Bz component of the magnetic field (kG) + // Returns Bz component of the magnetic field (kG) // - if (AliTracker::UniformField()) return AliTracker::GetBz(); - Double_t r[3]; GetXYZ(r); - return AliTracker::GetBz(r); -} + if (AliTracker::UniformField()) { + return AliTracker::GetBz(); + } + Double_t r[3]; + GetXYZ(r); + + return AliTracker::GetBz(r); +}