X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=STEER%2FAliESDtrack.cxx;h=b048d0dc4ec305b172365f2895dc40a5bd70dac6;hp=6fbe66e8a6ac5844bdc1cf9efa14b4dfb50ae566;hb=70d0b23e438246ab2b407e74b03afa4cd2f8113a;hpb=672b5f431b0c30fd97dc3e080fa06accab03bc9a diff --git a/STEER/AliESDtrack.cxx b/STEER/AliESDtrack.cxx index 6fbe66e8a6a..b048d0dc4ec 100644 --- a/STEER/AliESDtrack.cxx +++ b/STEER/AliESDtrack.cxx @@ -12,83 +12,582 @@ * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ - //----------------------------------------------------------------- // Implementation of the ESD track class // ESD = Event Summary Data -// This is the class to deal with during the phisical analysis of data +// This is the class to deal with during the phisics analysis of data // Origin: Iouri Belikov, CERN // e-mail: Jouri.Belikov@cern.ch //----------------------------------------------------------------- -#include "TMath.h" +#include +#include +#include "AliESDVertex.h" #include "AliESDtrack.h" #include "AliKalmanTrack.h" +#include "AliLog.h" +#include "AliTrackPointArray.h" ClassImp(AliESDtrack) +void SetPIDValues(Double_t * dest, const Double_t * src, Int_t n) { + // This function copies "n" PID weights from "scr" to "dest" + // and normalizes their sum to 1 thus producing conditional probabilities. + // The negative weights are set to 0. + // In case all the weights are non-positive they are replaced by + // uniform probabilities + + if (n<=0) return; + + Float_t uniform = 1./(Float_t)n; + + Float_t sum = 0; + for (Int_t i=0; i=0) { + sum+=src[i]; + dest[i] = src[i]; + } + else { + dest[i] = 0; + } + + if(sum>0) + for (Int_t i=0; iPhi()*180./TMath::Pi(); + if (alpha<0) alpha+= 360.; + if (alpha>360) alpha -= 360.; + + Int_t sector = (Int_t)(alpha/20.); + alpha = 10. + 20.*sector; + alpha /= 180; + alpha *= TMath::Pi(); + + // Covariance matrix: no errors, the parameters are exact + for (i=0; i<15; i++) covar[i]=0.; + + // Get the vertex of origin and the momentum + TVector3 ver(part->Vx(),part->Vy(),part->Vz()); + TVector3 mom(part->Px(),part->Py(),part->Pz()); + + // Rotate to the local coordinate system (TPC sector) + ver.RotateZ(-alpha); + mom.RotateZ(-alpha); + + // X of the referense plane + xref = ver.X(); + + Int_t pdgCode = part->GetPdgCode(); + + Double_t charge = + TDatabasePDG::Instance()->GetParticle(pdgCode)->Charge(); + + param[0] = ver.Y(); + param[1] = ver.Z(); + param[2] = TMath::Sin(mom.Phi()); + param[3] = mom.Pz()/mom.Pt(); + param[4] = TMath::Sign(1/mom.Pt(),charge); + + // Set AliExternalTrackParam + Set(xref, alpha, param, covar); + + // Set the PID + Int_t indexPID = 99; + + switch (TMath::Abs(pdgCode)) { + + case 11: // electron + indexPID = 0; + break; + + case 13: // muon + indexPID = 1; + break; + + case 211: // pion + indexPID = 2; + break; + + case 321: // kaon + indexPID = 3; + break; + + case 2212: // proton + indexPID = 4; + break; + + default: + break; + } + + // If the particle is not e,mu,pi,K or p the PID probabilities are set to 0 + if (indexPID < AliPID::kSPECIES) { + fR[indexPID]=1.; + fITSr[indexPID]=1.; + fTPCr[indexPID]=1.; + fTRDr[indexPID]=1.; + fTOFr[indexPID]=1.; + fHMPIDr[indexPID]=1.; + + } + // AliESD track label + SetLabel(part->GetUniqueID()); + +} + +//_______________________________________________________________________ +AliESDtrack::~AliESDtrack(){ + // + // This is destructor according Coding Conventrions + // + //printf("Delete track\n"); + delete fIp; + delete fTPCInner; + delete fOp; + delete fCp; + delete fFriendTrack; +} + +void AliESDtrack::AddCalibObject(TObject * object){ + // + // add calib object to the list + // + if (!fFriendTrack) fFriendTrack = new AliESDfriendTrack; + fFriendTrack->AddCalibObject(object); +} + +TObject * AliESDtrack::GetCalibObject(Int_t index){ + // + // return calib objct at given position + // + if (!fFriendTrack) return 0; + return fFriendTrack->GetCalibObject(index); +} + + +//_______________________________________________________________________ +void AliESDtrack::MakeMiniESDtrack(){ + // Resets everything except + // fFlags: Reconstruction status flags + // fLabel: Track label + // fID: Unique ID of the track + // fD: Impact parameter in XY-plane + // fZ: Impact parameter in Z + // fR[AliPID::kSPECIES]: combined "detector response probability" + // Running track parameters in the base class (AliExternalTrackParam) + + fTrackLength = 0; + + for (Int_t i=0;imax) {k=i; max=fR[i];} } - if (k==0) return 0.00051; - if (k==1) return 0.10566; - if (k==2||k==-1) return 0.13957; - if (k==3) return 0.49368; - if (k==4) return 0.93827; - Warning("GetMass()","Undefined mass !"); - return 0.13957; + if (k==0) { // dE/dx "crossing points" in the TPC + Double_t p=GetP(); + if ((p>0.38)&&(p<0.48)) + if (fR[0]0.75)&&(p<0.85)) + if (fR[0]GetLabel(); @@ -98,61 +597,124 @@ Bool_t AliESDtrack::UpdateTrackParams(AliKalmanTrack *t, ULong_t flags) { SetIntegratedLength(t->GetIntegratedLength()); } - fRalpha=t->GetAlpha(); - t->GetExternalParameters(fRx,fRp); - t->GetExternalCovariance(fRc); - + Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance()); + switch (flags) { case kITSin: case kITSout: case kITSrefit: + fITSClusterMap=0; fITSncls=t->GetNumberOfClusters(); + index=fFriendTrack->GetITSindices(); + for (Int_t i=0;iGetClusterIndex(i); + if (i> 28; + SETBIT(fITSClusterMap,l); + } + } fITSchi2=t->GetChi2(); - for (Int_t i=0;iGetClusterIndex(i); fITSsignal=t->GetPIDsignal(); + fITSLabel = t->GetLabel(); break; case kTPCin: case kTPCrefit: - fIalpha=fRalpha; - fIx=fRx; - { - Int_t i; - for (i=0; i<5; i++) fIp[i]=fRp[i]; - for (i=0; i<15;i++) fIc[i]=fRc[i]; - } + fTPCLabel = t->GetLabel(); + if (flags==kTPCin) fTPCInner=new AliExternalTrackParam(*t); + if (!fIp) fIp=new AliExternalTrackParam(*t); + else + fIp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance()); case kTPCout: - fTPCncls=t->GetNumberOfClusters(); + index=fFriendTrack->GetTPCindices(); + if (flags & kTPCout){ + if (!fOp) fOp=new AliExternalTrackParam(*t); + else + fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance()); + } + fTPCncls=t->GetNumberOfClusters(); fTPCchi2=t->GetChi2(); - for (Int_t i=0;iGetClusterIndex(i); + + {//prevrow must be declared in separate namespace, otherwise compiler cries: + //"jump to case label crosses initialization of `Int_t prevrow'" + Int_t prevrow = -1; + // for (Int_t i=0;iGetClusterIndex(i); + Int_t idx = index[i]; + + if (idx<0) continue; + + // Piotr's Cluster Map for HBT + // ### please change accordingly if cluster array is changing + // to "New TPC Tracking" style (with gaps in array) + Int_t sect = (idx&0xff000000)>>24; + Int_t row = (idx&0x00ff0000)>>16; + if (sect > 18) row +=63; //if it is outer sector, add number of inner sectors + + fTPCClusterMap.SetBitNumber(row,kTRUE); + + //Fill the gap between previous row and this row with 0 bits + //In case ### pleas change it as well - just set bit 0 in case there + //is no associated clusters for current "i" + if (prevrow < 0) + { + prevrow = row;//if previous bit was not assigned yet == this is the first one + } + else + { //we don't know the order (inner to outer or reverse) + //just to be save in case it is going to change + Int_t n = 0, m = 0; + if (prevrow < row) + { + n = prevrow; + m = row; + } + else + { + n = row; + m = prevrow; + } + + for (Int_t j = n+1; j < m; j++) + { + fTPCClusterMap.SetBitNumber(j,kFALSE); + } + prevrow = row; + } + // End Of Piotr's Cluster Map for HBT + } + } fTPCsignal=t->GetPIDsignal(); - {Double_t mass=t->GetMass(); // preliminary mass setting - if (mass>0.5) fR[4]=1.; // used by - else if (mass<0.4) fR[2]=1.; // the ITS reconstruction - else fR[3]=1.;} // break; - case kTRDout: - { //requested by the PHOS ("temporary solution") - Double_t r=474.; - t->PropagateTo(r,30.,0.); - fOalpha=fRalpha; - fOx=fRx; - Int_t i; - for (i=0; i<5; i++) fOp[i]=fRp[i]; - for (i=0; i<15;i++) fOc[i]=fRc[i]; - } - case kTRDin: case kTRDrefit: - fTRDncls=t->GetNumberOfClusters(); + case kTRDout: case kTRDin: case kTRDrefit: + index=fFriendTrack->GetTRDindices(); + fTRDLabel = t->GetLabel(); fTRDchi2=t->GetChi2(); - for (Int_t i=0;iGetClusterIndex(i); + fTRDncls=6;//t->GetNumberOfTracklets(); //t->GetNumberOfClusters(); + //for (Int_t i=0;iGetClusterIndex(i); + for (Int_t i=0;i<6;i++) index[i]=t->GetTrackletIndex(i); + fTRDsignal=t->GetPIDsignal(); break; - + case kTRDbackup: + if (!fOp) fOp=new AliExternalTrackParam(*t); + else + fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance()); + fTRDncls0 = t->GetNumberOfClusters(); + break; + case kTOFin: + break; + case kTOFout: + break; + case kTRDStop: + break; default: - Error("UpdateTrackParams()","Wrong flag !\n"); + AliError("Wrong flag !"); return kFALSE; } - return kTRUE; + return rc; } //_______________________________________________________________________ @@ -160,161 +722,318 @@ void AliESDtrack::GetExternalParameters(Double_t &x, Double_t p[5]) const { //--------------------------------------------------------------------- // This function returns external representation of the track parameters //--------------------------------------------------------------------- - x=fRx; - for (Int_t i=0; i<5; i++) p[i]=fRp[i]; + x=GetX(); + for (Int_t i=0; i<5; i++) p[i]=GetParameter()[i]; } -Double_t AliESDtrack::GetP() const { +//_______________________________________________________________________ +void AliESDtrack::GetExternalCovariance(Double_t cov[15]) const { //--------------------------------------------------------------------- - // This function returns the track momentum + // This function returns external representation of the cov. matrix //--------------------------------------------------------------------- - Double_t lam=TMath::ATan(fRp[3]); - Double_t pt=1./TMath::Abs(fRp[4]); - return pt/TMath::Cos(lam); + for (Int_t i=0; i<15; i++) cov[i]=AliExternalTrackParam::GetCovariance()[i]; } -void AliESDtrack::GetPxPyPz(Double_t *p) const { +//_______________________________________________________________________ +Bool_t AliESDtrack::GetConstrainedExternalParameters + (Double_t &alpha, Double_t &x, Double_t p[5]) const { //--------------------------------------------------------------------- - // This function returns the global track momentum components + // This function returns the constrained external track parameters //--------------------------------------------------------------------- - Double_t phi=TMath::ASin(fRp[2]) + fRalpha; - Double_t pt=1./TMath::Abs(fRp[4]); - p[0]=pt*TMath::Cos(phi); p[1]=pt*TMath::Sin(phi); p[2]=pt*fRp[3]; + if (!fCp) return kFALSE; + alpha=fCp->GetAlpha(); + x=fCp->GetX(); + for (Int_t i=0; i<5; i++) p[i]=fCp->GetParameter()[i]; + return kTRUE; } -void AliESDtrack::GetXYZ(Double_t *xyz) const { +//_______________________________________________________________________ +Bool_t +AliESDtrack::GetConstrainedExternalCovariance(Double_t c[15]) const { //--------------------------------------------------------------------- - // This function returns the global track position + // This function returns the constrained external cov. matrix //--------------------------------------------------------------------- - Double_t phi=TMath::ATan2(fRp[0],fRx) + fRalpha; - Double_t r=TMath::Sqrt(fRx*fRx + fRp[0]*fRp[0]); - xyz[0]=r*TMath::Cos(phi); xyz[1]=r*TMath::Sin(phi); xyz[2]=fRp[1]; + if (!fCp) return kFALSE; + for (Int_t i=0; i<15; i++) c[i]=fCp->GetCovariance()[i]; + return kTRUE; } -void AliESDtrack::GetInnerPxPyPz(Double_t *p) const { +Bool_t +AliESDtrack::GetInnerExternalParameters + (Double_t &alpha, Double_t &x, Double_t p[5]) const { //--------------------------------------------------------------------- - // This function returns the global track momentum components - // af the entrance of the TPC + // This function returns external representation of the track parameters + // at the inner layer of TPC //--------------------------------------------------------------------- - Double_t phi=TMath::ASin(fIp[2]) + fIalpha; - Double_t pt=1./TMath::Abs(fIp[4]); - p[0]=pt*TMath::Cos(phi); p[1]=pt*TMath::Sin(phi); p[2]=pt*fIp[3]; + if (!fIp) return kFALSE; + alpha=fIp->GetAlpha(); + x=fIp->GetX(); + for (Int_t i=0; i<5; i++) p[i]=fIp->GetParameter()[i]; + return kTRUE; } -void AliESDtrack::GetInnerXYZ(Double_t *xyz) const { - //--------------------------------------------------------------------- - // This function returns the global track position - // af the entrance of the TPC - //--------------------------------------------------------------------- - Double_t phi=TMath::ATan2(fIp[0],fIx) + fIalpha; - Double_t r=TMath::Sqrt(fIx*fIx + fIp[0]*fIp[0]); - xyz[0]=r*TMath::Cos(phi); xyz[1]=r*TMath::Sin(phi); xyz[2]=fIp[1]; +Bool_t +AliESDtrack::GetInnerExternalCovariance(Double_t cov[15]) const { + //--------------------------------------------------------------------- + // This function returns external representation of the cov. matrix + // at the inner layer of TPC + //--------------------------------------------------------------------- + if (!fIp) return kFALSE; + for (Int_t i=0; i<15; i++) cov[i]=fIp->GetCovariance()[i]; + return kTRUE; } -void AliESDtrack::GetOuterPxPyPz(Double_t *p) const { +Bool_t +AliESDtrack::GetOuterExternalParameters + (Double_t &alpha, Double_t &x, Double_t p[5]) const { //--------------------------------------------------------------------- - // This function returns the global track momentum components - // af the radius of the PHOS + // This function returns external representation of the track parameters + // at the inner layer of TRD //--------------------------------------------------------------------- - Double_t phi=TMath::ASin(fOp[2]) + fOalpha; - Double_t pt=1./TMath::Abs(fOp[4]); - p[0]=pt*TMath::Cos(phi); p[1]=pt*TMath::Sin(phi); p[2]=pt*fOp[3]; + if (!fOp) return kFALSE; + alpha=fOp->GetAlpha(); + x=fOp->GetX(); + for (Int_t i=0; i<5; i++) p[i]=fOp->GetParameter()[i]; + return kTRUE; } -void AliESDtrack::GetOuterXYZ(Double_t *xyz) const { - //--------------------------------------------------------------------- - // This function returns the global track position - // af the radius of the PHOS - //--------------------------------------------------------------------- - Double_t phi=TMath::ATan2(fOp[0],fOx) + fOalpha; - Double_t r=TMath::Sqrt(fOx*fOx + fOp[0]*fOp[0]); - xyz[0]=r*TMath::Cos(phi); xyz[1]=r*TMath::Sin(phi); xyz[2]=fOp[1]; +Bool_t +AliESDtrack::GetOuterExternalCovariance(Double_t cov[15]) const { + //--------------------------------------------------------------------- + // This function returns external representation of the cov. matrix + // at the inner layer of TRD + //--------------------------------------------------------------------- + if (!fOp) return kFALSE; + for (Int_t i=0; i<15; i++) cov[i]=fOp->GetCovariance()[i]; + return kTRUE; } -//_______________________________________________________________________ -void AliESDtrack::GetExternalCovariance(Double_t c[15]) const { - //--------------------------------------------------------------------- - // This function returns external representation of the cov. matrix - //--------------------------------------------------------------------- - for (Int_t i=0; i<15; i++) c[i]=fRc[i]; +Int_t AliESDtrack::GetNcls(Int_t idet) const +{ + // Get number of clusters by subdetector index + // + Int_t ncls = 0; + switch(idet){ + case 0: + ncls = fITSncls; + break; + case 1: + ncls = fTPCncls; + break; + case 2: + ncls = fTRDncls; + break; + case 3: + if (fTOFindex != 0) + ncls = 1; + break; + default: + break; + } + return ncls; +} + +Int_t AliESDtrack::GetClusters(Int_t idet, Int_t *idx) const +{ + // Get cluster index array by subdetector index + // + Int_t ncls = 0; + switch(idet){ + case 0: + ncls = GetITSclusters(idx); + break; + case 1: + ncls = GetTPCclusters(idx); + break; + case 2: + ncls = GetTRDclusters(idx); + break; + case 3: + if (fTOFindex != 0) { + idx[0] = GetTOFcluster(); + ncls = 1; + } + break; + default: + break; + } + return ncls; } //_______________________________________________________________________ void AliESDtrack::GetIntegratedTimes(Double_t *times) const { // Returns the array with integrated times for each particle hypothesis - for (Int_t i=0; iGetITSindices(); + for (Int_t i=0; iGetTPCindices(); + for (Int_t i=0; iGetTPCindices(); + for (Int_t i=row0;i<=row1;i++){ + Int_t idx = index[i]; + if (idx!=-1) good++; // track outside of dead zone + if (idx>0) found++; + } + Float_t density=0.5; + if (good>(row1-row0)*0.5) density = Float_t(found)/Float_t(good); + return density; +} + //_______________________________________________________________________ void AliESDtrack::SetTPCpid(const Double_t *p) { // Sets values for the probability of each particle type (in TPC) - for (Int_t i=0; iGetTRDindices(); + for (Int_t i=0; iGetTRDindices(); + for (Int_t i=0; i<6/*AliESDfriendTrack::kMaxTRDcluster*/; i++) idx[i]=index[i]; + } return fTRDncls; } //_______________________________________________________________________ void AliESDtrack::SetTRDpid(const Double_t *p) { // Sets values for the probability of each particle type (in TRD) - for (Int_t i=0; iGetXv()*cs + vtx->GetYv()*sn; + Double_t yv=-vtx->GetXv()*sn + vtx->GetYv()*cs, zv=vtx->GetZv(); + x-=xv; y-=yv; + + //Estimate the impact parameter neglecting the track curvature + Double_t d=TMath::Abs(x*snp - y*TMath::Sqrt(1.- snp*snp)); + if (d > maxd) return kFALSE; + + //Propagate to the DCA + Double_t crv=kB2C*b*GetParameter()[4]; + if (TMath::Abs(b) < kAlmost0Field) crv=0.; + + Double_t tgfv=-(crv*x - snp)/(crv*y + TMath::Sqrt(1.-snp*snp)); + sn=tgfv/TMath::Sqrt(1.+ tgfv*tgfv); + if (TMath::Abs(tgfv)>0.) cs = sn/tgfv; + else cs=1.; + + x = xv*cs + yv*sn; + yv=-xv*sn + yv*cs; xv=x; + + if (!Propagate(alpha+TMath::ASin(sn),xv,b)) return kFALSE; + + fD = GetParameter()[0] - yv; + fZ = GetParameter()[1] - zv; + + Double_t cov[6]; vtx->GetCovMatrix(cov); + + //***** Improvements by A.Dainese + alpha=GetAlpha(); sn=TMath::Sin(alpha); cs=TMath::Cos(alpha); + Double_t s2ylocvtx = cov[0]*sn*sn + cov[2]*cs*cs - 2.*cov[1]*cs*sn; + fCdd = GetCovariance()[0] + s2ylocvtx; // neglecting correlations + fCdz = GetCovariance()[1]; // between (x,y) and z + fCzz = GetCovariance()[2] + cov[5]; // in vertex's covariance matrix + //***** + + {//Try to constrain + Double_t p[2]={yv,zv}, c[3]={cov[2],0.,cov[5]}; + Double_t chi2=GetPredictedChi2(p,c); + + if (chi2>77.) return kFALSE; + + AliExternalTrackParam tmp(*this); + if (!tmp.Update(p,c)) return kFALSE; + + fCchi2=chi2; + if (!fCp) fCp=new AliExternalTrackParam(); + new (fCp) AliExternalTrackParam(tmp); + } + + return kTRUE; } +//_______________________________________________________________________ +void AliESDtrack::Print(Option_t *) const { + // Prints info on the track + + printf("ESD track info\n") ; + Double_t p[AliPID::kSPECIESN] ; + Int_t index = 0 ; + if( IsOn(kITSpid) ){ + printf("From ITS: ") ; + GetITSpid(p) ; + for(index = 0 ; index < AliPID::kSPECIES; index++) + printf("%f, ", p[index]) ; + printf("\n signal = %f\n", GetITSsignal()) ; + } + if( IsOn(kTPCpid) ){ + printf("From TPC: ") ; + GetTPCpid(p) ; + for(index = 0 ; index < AliPID::kSPECIES; index++) + printf("%f, ", p[index]) ; + printf("\n signal = %f\n", GetTPCsignal()) ; + } + if( IsOn(kTRDpid) ){ + printf("From TRD: ") ; + GetTRDpid(p) ; + for(index = 0 ; index < AliPID::kSPECIES; index++) + printf("%f, ", p[index]) ; + printf("\n signal = %f\n", GetTRDsignal()) ; + } + if( IsOn(kTOFpid) ){ + printf("From TOF: ") ; + GetTOFpid(p) ; + for(index = 0 ; index < AliPID::kSPECIES; index++) + printf("%f, ", p[index]) ; + printf("\n signal = %f\n", GetTOFsignal()) ; + } + if( IsOn(kHMPIDpid) ){ + printf("From HMPID: ") ; + GetHMPIDpid(p) ; + for(index = 0 ; index < AliPID::kSPECIES; index++) + printf("%f, ", p[index]) ; + printf("\n signal = %f\n", GetHMPIDsignal()) ; + } +} +