/************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * 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 phisics analysis of data // Origin: Iouri Belikov, CERN // e-mail: Jouri.Belikov@cern.ch //----------------------------------------------------------------- #include #include #include "AliESDVertex.h" #include "AliESDtrack.h" #include "AliKalmanTrack.h" #include "AliLog.h" #include "AliTrackPointArray.h" #include "TPolyMarker3D.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; if(fTRDnSlices) delete[] fTRDslices; } AliESDtrack &AliESDtrack::operator=(const AliESDtrack &source){ if(&source == this) return *this; AliExternalTrackParam::operator=(source); if(source.fCp){ // we have the trackparam: assign or copy construct if(fCp)*fCp = *source.fCp; else fCp = new AliExternalTrackParam(*source.fCp); } else{ // no track param delete the old one if(fCp)delete fCp; fCp = 0; } if(source.fIp){ // we have the trackparam: assign or copy construct if(fIp)*fIp = *source.fIp; else fIp = new AliExternalTrackParam(*source.fIp); } else{ // no track param delete the old one if(fIp)delete fIp; fIp = 0; } if(source.fTPCInner){ // we have the trackparam: assign or copy construct if(fTPCInner) *fTPCInner = *source.fTPCInner; else fTPCInner = new AliExternalTrackParam(*source.fTPCInner); } else{ // no track param delete the old one if(fTPCInner)delete fTPCInner; fTPCInner = 0; } if(source.fOp){ // we have the trackparam: assign or copy construct if(fOp) *fOp = *source.fOp; else fOp = new AliExternalTrackParam(*source.fOp); } else{ // no track param delete the old one if(fOp)delete fOp; fOp = 0; } // copy also the friend track // use copy constructor if(source.fFriendTrack){ // we have the trackparam: assign or copy construct delete fFriendTrack; fFriendTrack=new AliESDfriendTrack(*source.fFriendTrack); } else{ // no track param delete the old one delete fFriendTrack; fFriendTrack= 0; } fTPCClusterMap = source.fTPCClusterMap; fTPCSharedMap = source.fTPCSharedMap; // the simple stuff fFlags = source.fFlags; fID = source.fID; fLabel = source.fLabel; fITSLabel = source.fITSLabel; for(int i = 0; i< 12;++i){ fITSModule[i] = source.fITSModule[i]; } fTPCLabel = source.fTPCLabel; fTRDLabel = source.fTRDLabel; for(int i = 0; i< 3;++i){ fTOFLabel[i] = source.fTOFLabel[i]; } fTOFCalChannel = source.fTOFCalChannel; fTOFindex = source.fTOFindex; fHMPIDqn = source.fHMPIDqn; fHMPIDcluIdx = source.fHMPIDcluIdx; fEMCALindex = source.fEMCALindex; for(int i = 0; i< 3;++i){ fKinkIndexes[i] = source.fKinkIndexes[i]; fV0Indexes[i] = source.fV0Indexes[i]; } for(int i = 0; i< AliPID::kSPECIES;++i){ fR[i] = source.fR[i]; fITSr[i] = source.fITSr[i]; fTPCr[i] = source.fTPCr[i]; fTRDr[i] = source.fTRDr[i]; fTOFr[i] = source.fTOFr[i]; fHMPIDr[i] = source.fHMPIDr[i]; fTrackTime[i] = source.fTrackTime[i]; } fHMPIDtrkTheta = source.fHMPIDtrkTheta; fHMPIDtrkPhi = source.fHMPIDtrkPhi; fHMPIDsignal = source.fHMPIDsignal; fTrackLength = source. fTrackLength; fD = source.fD; fZ = source.fZ; fCdd = source.fCdd; fCdz = source.fCdz; fCzz = source.fCzz; fCchi2 = source.fCchi2; fITSchi2 = source.fITSchi2; fTPCchi2 = source.fTPCchi2; fTRDchi2 = source.fTRDchi2; fTOFchi2 = source.fTOFchi2; fHMPIDchi2 = source.fHMPIDchi2; fITSsignal = source.fITSsignal; fTPCsignal = source.fTPCsignal; fTPCsignalS = source.fTPCsignalS; for(int i = 0; i< 4;++i){ fTPCPoints[i] = source.fTPCPoints[i]; } fTRDsignal = source.fTRDsignal; for(int i = 0;i < kTRDnPlanes;++i){ fTRDTimBin[i] = source.fTRDTimBin[i]; } if(fTRDnSlices) delete[] fTRDslices; fTRDslices=0; fTRDnSlices=source.fTRDnSlices; if (fTRDnSlices) { fTRDslices=new Double32_t[fTRDnSlices]; for(int j = 0;j < fTRDnSlices;++j) fTRDslices[j] = source.fTRDslices[j]; } fTRDQuality = source.fTRDQuality; fTRDBudget = source.fTRDBudget; fTOFsignal = source.fTOFsignal; fTOFsignalToT = source.fTOFsignalToT; fTOFsignalRaw = source.fTOFsignalRaw; fTOFsignalDz = source.fTOFsignalDz; for(int i = 0;i<10;++i){ fTOFInfo[i] = source.fTOFInfo[i]; } fHMPIDtrkX = source.fHMPIDtrkX; fHMPIDtrkY = source.fHMPIDtrkY; fHMPIDmipX = source.fHMPIDmipX; fHMPIDmipY = source.fHMPIDmipY; fTPCncls = source.fTPCncls; fTPCnclsF = source.fTPCnclsF; fTPCsignalN = source.fTPCsignalN; fITSncls = source.fITSncls; fITSClusterMap = source.fITSClusterMap; fTRDncls = source.fTRDncls; fTRDncls0 = source.fTRDncls0; fTRDpidQuality = source.fTRDpidQuality; return *this; } void AliESDtrack::Copy(TObject &obj) const { // this overwrites the virtual TOBject::Copy() // to allow run time copying without casting // in AliESDEvent if(this==&obj)return; AliESDtrack *robj = dynamic_cast(&obj); if(!robj)return; // not an AliESDtrack *robj = *this; } 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) { // 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(); if (t->IsStartedTimeIntegral()) { SetStatus(kTIME); Double_t times[10];t->GetIntegratedTimes(times); SetIntegratedTimes(times); SetIntegratedLength(t->GetIntegratedLength()); } 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(); fITSsignal=t->GetPIDsignal(); fITSLabel = t->GetLabel(); break; case kTPCin: case kTPCrefit: 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: 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(); {//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(); break; case kTRDout: case kTRDin: case kTRDrefit: index = fFriendTrack->GetTRDindices(); fTRDLabel = t->GetLabel(); fTRDchi2 = t->GetChi2(); fTRDncls = t->GetNumberOfClusters(); 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: AliError("Wrong flag !"); return kFALSE; } return rc; } //_______________________________________________________________________ void AliESDtrack::GetExternalParameters(Double_t &x, Double_t p[5]) const { //--------------------------------------------------------------------- // This function returns external representation of the track parameters //--------------------------------------------------------------------- x=GetX(); for (Int_t i=0; i<5; i++) p[i]=GetParameter()[i]; } //_______________________________________________________________________ void AliESDtrack::GetExternalCovariance(Double_t cov[15]) const { //--------------------------------------------------------------------- // This function returns external representation of the cov. matrix //--------------------------------------------------------------------- for (Int_t i=0; i<15; i++) cov[i]=AliExternalTrackParam::GetCovariance()[i]; } //_______________________________________________________________________ Bool_t AliESDtrack::GetConstrainedExternalParameters (Double_t &alpha, Double_t &x, Double_t p[5]) const { //--------------------------------------------------------------------- // This function returns the constrained external track parameters //--------------------------------------------------------------------- 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; } //_______________________________________________________________________ Bool_t AliESDtrack::GetConstrainedExternalCovariance(Double_t c[15]) const { //--------------------------------------------------------------------- // This function returns the constrained external cov. matrix //--------------------------------------------------------------------- if (!fCp) return kFALSE; for (Int_t i=0; i<15; i++) c[i]=fCp->GetCovariance()[i]; return kTRUE; } Bool_t AliESDtrack::GetInnerExternalParameters (Double_t &alpha, Double_t &x, Double_t p[5]) const { //--------------------------------------------------------------------- // This function returns external representation of the track parameters // at the inner layer of TPC //--------------------------------------------------------------------- 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; } 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; } Bool_t AliESDtrack::GetOuterExternalParameters (Double_t &alpha, Double_t &x, Double_t p[5]) const { //--------------------------------------------------------------------- // This function returns external representation of the track parameters // at the inner layer of TRD //--------------------------------------------------------------------- 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; } 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; } 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) SetPIDValues(fTPCr,p,AliPID::kSPECIES); SetStatus(AliESDtrack::kTPCpid); } //_______________________________________________________________________ void AliESDtrack::GetTPCpid(Double_t *p) const { // Gets 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) SetPIDValues(fTRDr,p,AliPID::kSPECIES); SetStatus(AliESDtrack::kTRDpid); } //_______________________________________________________________________ void AliESDtrack::GetTRDpid(Double_t *p) const { // Gets the probability of each particle type (in TRD) for (Int_t i=0; i=kTRDnPlanes)) { AliError("Wrong TRD plane !"); return; } if ((slice<0) || (slice>=ns)) { AliError("Wrong TRD slice !"); return; } Int_t n=plane*ns + slice; fTRDslices[n]=q; } Double_t AliESDtrack::GetTRDslice(Int_t plane, Int_t slice) const { //Gets the charge from the slice of the plane Int_t ns=GetNumberOfTRDslices(); if (ns==0) { //AliError("No TRD slices allocated for this track !"); return -1.; } if ((plane<0) || (plane>=kTRDnPlanes)) { AliError("Wrong TRD plane !"); return -1.; } if ((slice<-1) || (slice>=ns)) { //AliError("Wrong TRD slice !"); return -1.; } if (slice==-1) { Double_t q=0.; for (Int_t i=0; iGetCovMatrix(covar); Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]}; Double_t c[3]={covar[2],0.,covar[5]}; Double_t chi2=GetPredictedChi2(p,c); if (chi2>77.) return kFALSE; delete fCp; fCp=new AliExternalTrackParam(*this); if (!fCp->Update(p,c)) {delete fCp; fCp=0; return kFALSE;} fCchi2=chi2; 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()) ; } } // // Draw functionality // Origin: Marian Ivanov, Marian.Ivanov@cern.ch // void AliESDtrack::FillPolymarker(TPolyMarker3D *pol, Float_t magF, Float_t minR, Float_t maxR, Float_t stepR){ // // Fill points in the polymarker // TObjArray arrayRef; arrayRef.AddLast(new AliExternalTrackParam(*this)); if (fIp) arrayRef.AddLast(new AliExternalTrackParam(*fIp)); if (fOp) arrayRef.AddLast(new AliExternalTrackParam(*fOp)); // Double_t mpos[3]={0,0,0}; Int_t entries=arrayRef.GetEntries(); for (Int_t i=0;iGetXYZ(pos); mpos[0]+=pos[0]/entries; mpos[1]+=pos[1]/entries; mpos[2]+=pos[2]/entries; } // Rotate to the mean position // Float_t fi= TMath::ATan2(mpos[1],mpos[0]); for (Int_t i=0;iRotate(fi); if (!res) delete arrayRef.RemoveAt(i); } Int_t counter=0; for (Double_t r=minR; rGetXYZAt(r,magF,point)){ Double_t weight = 1./(10.+(r-param->GetX())*(r-param->GetX())); sweight+=weight; mlpos[0]+=point[0]*weight; mlpos[1]+=point[1]*weight; mlpos[2]+=point[2]*weight; } } if (sweight>0){ mlpos[0]/=sweight; mlpos[1]/=sweight; mlpos[2]/=sweight; pol->SetPoint(counter,mlpos[0],mlpos[1], mlpos[2]); printf("xyz\t%f\t%f\t%f\n",mlpos[0], mlpos[1],mlpos[2]); counter++; } } }