/************************************************************************** * 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 // // // // What do you need to know before starting analysis // (by Marian Ivanov: marian.ivanov@cern.ch) // // // AliESDtrack: // 1. What is the AliESDtrack // 2. What informations do we store // 3. How to use the information for analysis // // // 1.AliESDtrack is the container of the information about the track/particle // reconstructed during Barrel Tracking. // The track information is propagated from one tracking detector to // other using the functionality of AliESDtrack - Current parameters. // // No global fit model is used. // Barrel tracking use Kalman filtering technique, it gives optimal local // track parameters at given point under certian assumptions. // // Kalman filter take into account additional effect which are // difficult to handle using global fit. // Effects: // a.) Multiple scattering // b.) Energy loss // c.) Non homogenous magnetic field // // In general case, following barrel detectors are contributing to // the Kalman track information: // a. TPC // b. ITS // c. TRD // // In general 3 reconstruction itteration are performed: // 1. Find tracks - sequence TPC->ITS // 2. PropagateBack - sequence ITS->TPC->TRD -> Outer PID detectors // 3. Refit invward - sequence TRD->TPC->ITS // The current tracks are updated after each detector (see bellow). // In specical cases a track sanpshots are stored. // // // For some type of analysis (+visualization) track local parameters at // different position are neccesary. A snapshots during the track // propagation are created. // (See AliExternalTrackParam class for desctiption of variables and // functionality) // Snapshots: // a. Current parameters - class itself (AliExternalTrackParam) // Contributors: general case TRD->TPC->ITS // Preferable usage: Decission - primary or secondary track // NOTICE - By default the track parameters are stored at the DCA point // to the primary vertex. optimal for primary tracks, // far from optimal for secondary tracks. // b. Constrained parameters - Kalman information updated with // the Primary vertex information // Contributors: general case TRD->TPC->ITS // Preferable usage: Use only for tracks selected as primary // NOTICE - not real constrain - taken as additional measurement // with corresponding error // Function: // const AliExternalTrackParam *GetConstrainedParam() const {return fCp;} // c. Inner parameters - Track parameters at inner wall of the TPC // Contributors: general case TRD->TPC // function: // const AliExternalTrackParam *GetInnerParam() const { return fIp;} // // d. TPCinnerparam - contributors - TPC only // Contributors: TPC // Preferable usage: Requested for HBT study // (smaller correlations as using also ITS information) // NOTICE - the track parameters are propagated to the DCA to // to primary vertex // Optimal for primary, far from optimal for secondary tracks // Function: // const AliExternalTrackParam *GetTPCInnerParam() const {return fTPCInner;} // // e. Outer parameters - // Contributors- general case - ITS-> TPC -> TRD // The last point - Outer parameters radius is determined // e.a) Local inclination angle bigger than threshold - // Low momenta tracks // e.a) Catastrofic energy losss in material // e.b) Not further improvement (no space points) // Usage: // a.) Tracking: Starting parameter for Refit inward // b.) Visualization // c.) QA // NOTICE: Should be not used for the physic analysis // Function: // const AliExternalTrackParam *GetOuterParam() const { return fOp;} // //----------------------------------------------------------------- #include #include #include #include #include "AliESDVertex.h" #include "AliESDtrack.h" #include "AliESDEvent.h" #include "AliKalmanTrack.h" #include "AliVTrack.h" #include "AliLog.h" #include "AliTrackPointArray.h" #include "TPolyMarker3D.h" #include "AliTrackerBase.h" #include "AliTPCdEdxInfo.h" #include "AliDetectorPID.h" #include "TTreeStream.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; iInheritsFrom("AliExternalTrackParam")) { AliError("This is not a copy constructor. Use AliESDtrack(const AliESDtrack &) !"); AliWarning("Calling the default constructor..."); AliESDtrack(); return; } // Reset all the arrays Int_t i; for (i=kNITSchi2Std;i--;) fITSchi2Std[i] = 0; for (i=0; iGetID()); // Set ITS cluster map fITSClusterMap=track->GetITSClusterMap(); fITSSharedMap=0; fITSncls=0; for(i=0; i<6; i++) { if(HasPointOnITSLayer(i)) fITSncls++; } // Set TPC ncls fTPCncls=track->GetTPCNcls(); fTPCnclsF=track->GetTPCNclsF(); // TPC cluster maps const TBits* bmap = track->GetTPCClusterMapPtr(); if (bmap) SetTPCClusterMap(*bmap); bmap = GetTPCFitMapPtr(); if (bmap) SetTPCFitMap(*bmap); bmap = GetTPCSharedMapPtr(); if (bmap) SetTPCSharedMap(*bmap); // // Set the combined PID const Double_t *pid = track->PID(); if(pid) for (i=0; iGetEMCALcluster()); // AliESD track label // // PID info fITSsignal = track->GetITSsignal(); double itsdEdx[4]; track->GetITSdEdxSamples(itsdEdx); SetITSdEdxSamples(itsdEdx); // SetTPCsignal(track->GetTPCsignal(),fTPCsignalS,track->GetTPCsignalN()); // No signalS in AODPi AliTPCdEdxInfo * dEdxInfo = track->GetTPCdEdxInfo(); if (dEdxInfo) SetTPCdEdxInfo(new AliTPCdEdxInfo(*dEdxInfo)); // SetTRDsignal(track->GetTRDsignal()); int ntrdsl = track->GetNumberOfTRDslices(); if (ntrdsl>0) { SetNumberOfTRDslices((ntrdsl+2)*kTRDnPlanes); for (int ipl=kTRDnPlanes;ipl--;){ for (int isl=ntrdsl;isl--;) SetTRDslice(track->GetTRDslice(ipl,isl),ipl,isl); Double_t sp, p = track->GetTRDmomentum(ipl, &sp); SetTRDmomentum(p, ipl, &sp); } } // fTRDncls = track->GetTRDncls(); fTRDntracklets &= 0xff & track->GetTRDntrackletsPID(); fTRDchi2 = track->GetTRDchi2(); // SetTOFsignal(track->GetTOFsignal()); Double_t expt[AliPID::kSPECIES]; track->GetIntegratedTimes(expt); SetIntegratedTimes(expt); // SetTrackPhiEtaPtOnEMCal(track->GetTrackPhiOnEMCal(),track->GetTrackEtaOnEMCal(),track->GetTrackPtOnEMCal()); // SetLabel(track->GetLabel()); // Set the status SetStatus(track->GetStatus()); } //_______________________________________________________________________ AliESDtrack::AliESDtrack(TParticle * part) : AliExternalTrackParam(), fCp(0), fIp(0), fTPCInner(0), fOp(0), fHMPIDp(0), fFriendTrack(0), fTPCFitMap(159),//number of padrows fTPCClusterMap(159),//number of padrows fTPCSharedMap(159),//number of padrows fFlags(0), fID(0), fLabel(0), fITSLabel(0), fTPCLabel(0), fTRDLabel(0), fTOFCalChannel(-1), fTOFindex(-1), fHMPIDqn(0), fHMPIDcluIdx(-1), fCaloIndex(kEMCALNoMatch), fHMPIDtrkTheta(0), fHMPIDtrkPhi(0), fHMPIDsignal(0), fTrackLength(0), fdTPC(0),fzTPC(0), fCddTPC(0),fCdzTPC(0),fCzzTPC(0), fCchi2TPC(0), fD(0),fZ(0), fCdd(0),fCdz(0),fCzz(0), fCchi2(0), fITSchi2(0), fTPCchi2(0), fTPCchi2Iter1(0), fTRDchi2(0), fTOFchi2(0), fHMPIDchi2(0), fGlobalChi2(0), fITSsignal(0), fTPCsignal(0), fTPCsignalTuned(0), fTPCsignalS(0), fTPCdEdxInfo(0), fTRDsignal(0), fTRDQuality(0), fTRDBudget(0), fTOFsignal(99999), fTOFsignalTuned(99999), fTOFsignalToT(99999), fTOFsignalRaw(99999), fTOFsignalDz(999), fTOFsignalDx(999), fTOFdeltaBC(999), fTOFl0l1(999), fCaloDx(0), fCaloDz(0), fHMPIDtrkX(0), fHMPIDtrkY(0), fHMPIDmipX(0), fHMPIDmipY(0), fTPCncls(0), fTPCnclsF(0), fTPCsignalN(0), fTPCnclsIter1(0), fTPCnclsFIter1(0), fITSncls(0), fITSClusterMap(0), fITSSharedMap(0), fTRDncls(0), fTRDncls0(0), fTRDntracklets(0), fTRDNchamberdEdx(0), fTRDNclusterdEdx(0), fTRDnSlices(0), fTRDslices(0x0), fVertexID(-2), // -2 means an orphan track fESDEvent(0), fCacheNCrossedRows(-10), fCacheChi2TPCConstrainedVsGlobal(-10), fCacheChi2TPCConstrainedVsGlobalVertex(0), fDetectorPID(0x0), fTrackPhiOnEMCal(-999), fTrackEtaOnEMCal(-999), fTrackPtOnEMCal(-999) { // // ESD track from TParticle // // Reset all the arrays Int_t i; for (i=kNITSchi2Std;i--;) fITSchi2Std[i] = 0; for (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 fHMPIDp; delete fCp; delete fFriendTrack; delete fTPCdEdxInfo; if(fTRDnSlices) delete[] fTRDslices; //Reset cached values - needed for TClonesArray in AliESDInputHandler fCacheNCrossedRows = -10.; fCacheChi2TPCConstrainedVsGlobal = -10.; if(fCacheChi2TPCConstrainedVsGlobalVertex) fCacheChi2TPCConstrainedVsGlobalVertex = 0; delete fDetectorPID; } 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 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 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 delete fTPCInner; fTPCInner = 0; } if(source.fTPCdEdxInfo) { if(fTPCdEdxInfo) *fTPCdEdxInfo = *source.fTPCdEdxInfo; fTPCdEdxInfo = new AliTPCdEdxInfo(*source.fTPCdEdxInfo); } 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 delete fOp; fOp = 0; } if(source.fHMPIDp){ // we have the trackparam: assign or copy construct if(fHMPIDp) *fHMPIDp = *source.fHMPIDp; else fHMPIDp = new AliExternalTrackParam(*source.fHMPIDp); } else{ // no track param delete the old one delete fHMPIDp; fHMPIDp = 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; } fTPCFitMap = source.fTPCFitMap; 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; fCaloIndex = source.fCaloIndex; for (int i=kNITSchi2Std;i--;) fITSchi2Std[i] = source.fITSchi2Std[i]; 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; fdTPC = source.fdTPC; fzTPC = source.fzTPC; fCddTPC = source.fCddTPC; fCdzTPC = source.fCdzTPC; fCzzTPC = source.fCzzTPC; fCchi2TPC = source.fCchi2TPC; fD = source.fD; fZ = source.fZ; fCdd = source.fCdd; fCdz = source.fCdz; fCzz = source.fCzz; fCchi2 = source.fCchi2; fITSchi2 = source.fITSchi2; fTPCchi2 = source.fTPCchi2; fTPCchi2Iter1 = source.fTPCchi2Iter1; fTRDchi2 = source.fTRDchi2; fTOFchi2 = source.fTOFchi2; fHMPIDchi2 = source.fHMPIDchi2; fGlobalChi2 = source.fGlobalChi2; fITSsignal = source.fITSsignal; for (Int_t i=0;i<4;i++) {fITSdEdxSamples[i]=source.fITSdEdxSamples[i];} fTPCsignal = source.fTPCsignal; fTPCsignalTuned = source.fTPCsignalTuned; fTPCsignalS = source.fTPCsignalS; for(int i = 0; i< 4;++i){ fTPCPoints[i] = source.fTPCPoints[i]; } fTRDsignal = source.fTRDsignal; fTRDNchamberdEdx = source.fTRDNchamberdEdx; fTRDNclusterdEdx = source.fTRDNclusterdEdx; 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; fTOFsignalTuned = source.fTOFsignalTuned; fTOFsignalToT = source.fTOFsignalToT; fTOFsignalRaw = source.fTOFsignalRaw; fTOFsignalDz = source.fTOFsignalDz; fTOFsignalDx = source.fTOFsignalDx; fTOFdeltaBC = source.fTOFdeltaBC; fTOFl0l1 = source.fTOFl0l1; 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; fTPCnclsIter1 = source.fTPCnclsIter1; fTPCnclsFIter1 = source.fTPCnclsFIter1; fITSncls = source.fITSncls; fITSClusterMap = source.fITSClusterMap; fITSSharedMap = source.fITSSharedMap; fTRDncls = source.fTRDncls; fTRDncls0 = source.fTRDncls0; fTRDntracklets = source.fTRDntracklets; fVertexID = source.fVertexID; fCacheNCrossedRows = source.fCacheNCrossedRows; fCacheChi2TPCConstrainedVsGlobal = source.fCacheChi2TPCConstrainedVsGlobal; fCacheChi2TPCConstrainedVsGlobalVertex = source.fCacheChi2TPCConstrainedVsGlobalVertex; delete fDetectorPID; fDetectorPID=0x0; if (source.fDetectorPID) fDetectorPID = new AliDetectorPID(*source.fDetectorPID); fTrackPhiOnEMCal= source.fTrackPhiOnEMCal; fTrackEtaOnEMCal= source.fTrackEtaOnEMCal; fTrackPtOnEMCal= source.fTrackPtOnEMCal; 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; if (!fFriendTrack) return; fFriendTrack->AddCalibObject(object); } TObject * AliESDtrack::GetCalibObject(Int_t index){ // // return calib objct at given position // if (!fFriendTrack) return 0; return fFriendTrack->GetCalibObject(index); } Bool_t AliESDtrack::FillTPCOnlyTrack(AliESDtrack &track){ // Fills the information of the TPC-only first reconstruction pass // into the passed ESDtrack object. For consistency fTPCInner is also filled // again // For data produced before r26675 // RelateToVertexTPC was not properly called during reco // so you'll have to call it again, before FillTPCOnlyTrack // Float_t p[2],cov[3]; // track->GetImpactParametersTPC(p,cov); // if(p[0]==0&&p[1]==0) // <- Default values // track->RelateToVertexTPC(esd->GetPrimaryVertexTPC(),esd->GetMagneticField(),kVeryBig); if(!fTPCInner)return kFALSE; // fill the TPC track params to the global track parameters track.Set(fTPCInner->GetX(),fTPCInner->GetAlpha(),fTPCInner->GetParameter(),fTPCInner->GetCovariance()); track.fD = fdTPC; track.fZ = fzTPC; track.fCdd = fCddTPC; track.fCdz = fCdzTPC; track.fCzz = fCzzTPC; // copy the inner params if(track.fIp) *track.fIp = *fIp; else track.fIp = new AliExternalTrackParam(*fIp); // copy the TPCinner parameters if(track.fTPCInner) *track.fTPCInner = *fTPCInner; else track.fTPCInner = new AliExternalTrackParam(*fTPCInner); track.fdTPC = fdTPC; track.fzTPC = fzTPC; track.fCddTPC = fCddTPC; track.fCdzTPC = fCdzTPC; track.fCzzTPC = fCzzTPC; track.fCchi2TPC = fCchi2TPC; // copy all other TPC specific parameters // replace label by TPC label track.fLabel = fTPCLabel; track.fTPCLabel = fTPCLabel; track.fTPCchi2 = fTPCchi2; track.fTPCchi2Iter1 = fTPCchi2Iter1; track.fTPCsignal = fTPCsignal; track.fTPCsignalTuned = fTPCsignalTuned; track.fTPCsignalS = fTPCsignalS; for(int i = 0;i<4;++i)track.fTPCPoints[i] = fTPCPoints[i]; track.fTPCncls = fTPCncls; track.fTPCnclsF = fTPCnclsF; track.fTPCsignalN = fTPCsignalN; track.fTPCnclsIter1 = fTPCnclsIter1; track.fTPCnclsFIter1 = fTPCnclsFIter1; // PID for(int i=0;imax) {k=i; max=prob[i];} // if (k==0) { // dE/dx "crossing points" in the TPC Double_t p=GetP(); if ((p>0.38)&&(p<0.48)) if (prob[0]0.75)&&(p<0.85)) if (prob[0]kAlmost0 ? // account for curvature 2./curv*TMath::ASin(kRTOF*curv/2.)*TMath::Sqrt(1.+GetTgl()*GetTgl()) : kRTOF; tdif -= path/kCSpeed*TMath::Sqrt(1.+m*m/(p*p)); } bcid = TMath::Nint((tdif - kShift)/kSpacing); return bcid; } //______________________________________________________________________________ Double_t AliESDtrack::M() const { // Returns the assumed mass // (the pion mass, if the particle can't be identified properly). static Bool_t printerr=kTRUE; if (printerr) { AliWarning("WARNING !!! ... THIS WILL BE PRINTED JUST ONCE !!!"); printerr = kFALSE; AliWarning("This is the ESD mass. Use it with care !"); } return GetMass(); } //______________________________________________________________________________ Double_t AliESDtrack::E() const { // Returns the energy of the particle given its assumed mass. // Assumes the pion mass if the particle can't be identified properly. Double_t m = M(); Double_t p = P(); return TMath::Sqrt(p*p + m*m); } //______________________________________________________________________________ Double_t AliESDtrack::Y() const { // Returns the rapidity of a particle given its assumed mass. // Assumes the pion mass if the particle can't be identified properly. Double_t e = E(); Double_t pz = Pz(); if (e != TMath::Abs(pz)) { // energy was not equal to pz return 0.5*TMath::Log((e+pz)/(e-pz)); } else { // energy was equal to pz return -999.; } } //_______________________________________________________________________ Bool_t AliESDtrack::UpdateTrackParams(const AliKalmanTrack *t, ULong_t flags){ // // This function updates track's running parameters // Bool_t rc=kTRUE; SetStatus(flags); fLabel=t->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()); if (fFriendTrack) { if (flags==kITSout) fFriendTrack->SetITSOut(*t); if (flags==kTPCout) fFriendTrack->SetTPCOut(*t); if (flags==kTRDrefit) fFriendTrack->SetTRDIn(*t); } switch (flags) { case kITSin: fITSchi2Std[0] = t->GetChi2(); // case kITSout: fITSchi2Std[1] = t->GetChi2(); case kITSrefit: { fITSchi2Std[2] = t->GetChi2(); fITSClusterMap=0; fITSncls=t->GetNumberOfClusters(); if (fFriendTrack) { Int_t* indexITS = new Int_t[AliESDfriendTrack::kMaxITScluster]; for (Int_t i=0;iGetClusterIndex(i); if (i> 28; SETBIT(fITSClusterMap,l); } } fFriendTrack->SetITSIndices(indexITS,AliESDfriendTrack::kMaxITScluster); delete [] indexITS; } fITSchi2=t->GetChi2(); fITSsignal=t->GetPIDsignal(); fITSLabel = t->GetLabel(); // keep in fOp the parameters outside ITS for ITS stand-alone tracks if (flags==kITSout) { if (!fOp) fOp=new AliExternalTrackParam(*t); else fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance()); } } break; case kTPCin: case kTPCrefit: { fTPCLabel = t->GetLabel(); if (flags==kTPCin) { fTPCInner=new AliExternalTrackParam(*t); fTPCnclsIter1=t->GetNumberOfClusters(); fTPCchi2Iter1=t->GetChi2(); } if (!fIp) fIp=new AliExternalTrackParam(*t); else fIp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance()); } // Intentionally no break statement; need to set general TPC variables as well case kTPCout: { 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(); if (fFriendTrack) { // Copy cluster indices Int_t* indexTPC = new Int_t[AliESDfriendTrack::kMaxTPCcluster]; for (Int_t i=0;iGetClusterIndex(i); fFriendTrack->SetTPCIndices(indexTPC,AliESDfriendTrack::kMaxTPCcluster); delete [] indexTPC; } fTPCsignal=t->GetPIDsignal(); } break; case kTRDin: case kTRDrefit: break; case kTRDout: { fTRDLabel = t->GetLabel(); fTRDchi2 = t->GetChi2(); fTRDncls = t->GetNumberOfClusters(); if (fFriendTrack) { Int_t* indexTRD = new Int_t[AliESDfriendTrack::kMaxTRDcluster]; for (Int_t i=0;iGetTrackletIndex(i); fFriendTrack->SetTRDIndices(indexTRD,AliESDfriendTrack::kMaxTRDcluster); delete [] indexTRD; } //commented out by Xianguo //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; case kHMPIDout: if (!fHMPIDp) fHMPIDp=new AliExternalTrackParam(*t); else fHMPIDp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance()); 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; } void AliESDtrack::SetOuterParam(const AliExternalTrackParam *p, ULong_t flags) { // // This is a direct setter for the outer track parameters // SetStatus(flags); if (fOp) delete fOp; fOp=new AliExternalTrackParam(*p); } void AliESDtrack::SetOuterHmpParam(const AliExternalTrackParam *p, ULong_t flags) { // // This is a direct setter for the outer track parameters // SetStatus(flags); if (fHMPIDp) delete fHMPIDp; fHMPIDp=new AliExternalTrackParam(*p); } 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::GetOuterHmpExternalParameters (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 (!fHMPIDp) return kFALSE; alpha=fHMPIDp->GetAlpha(); x=fHMPIDp->GetX(); for (Int_t i=0; i<5; i++) p[i]=fHMPIDp->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; } Bool_t AliESDtrack::GetOuterHmpExternalCovariance(Double_t cov[15]) const { //--------------------------------------------------------------------- // This function returns external representation of the cov. matrix // at the inner layer of TRD //--------------------------------------------------------------------- if (!fHMPIDp) return kFALSE; for (Int_t i=0; i<15; i++) cov[i]=fHMPIDp->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 != -1) ncls = 1; break; case 4: //PHOS break; case 5: //HMPID if ((fHMPIDcluIdx >= 0) && (fHMPIDcluIdx < 7000000)) { if ((fHMPIDcluIdx%1000000 != 9999) && (fHMPIDcluIdx%1000000 != 99999)) { 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 != -1) { idx[0] = fTOFindex; ncls = 1; } break; case 4: //PHOS break; case 5: if ((fHMPIDcluIdx >= 0) && (fHMPIDcluIdx < 7000000)) { if ((fHMPIDcluIdx%1000000 != 9999) && (fHMPIDcluIdx%1000000 != 99999)) { idx[0] = GetHMPIDcluIdx(); ncls = 1; } } break; case 6: //EMCAL 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; i=fITSncls) && (i<6) ) idx[i]=-1; else { if (index) { idx[i]=index[i]; } else idx[i]= -2; } } } return fITSncls; } //_______________________________________________________________________ Bool_t AliESDtrack::GetITSModuleIndexInfo(Int_t ilayer,Int_t &idet,Int_t &status, Float_t &xloc,Float_t &zloc) const { //---------------------------------------------------------------------- // This function encodes in the module number also the status of cluster association // "status" can have the following values: // 1 "found" (cluster is associated), // 2 "dead" (module is dead from OCDB), // 3 "skipped" (module or layer forced to be skipped), // 4 "outinz" (track out of z acceptance), // 5 "nocls" (no clusters in the road), // 6 "norefit" (cluster rejected during refit), // 7 "deadzspd" (holes in z in SPD) // Also given are the coordinates of the crossing point of track and module // (in the local module ref. system) // WARNING: THIS METHOD HAS TO BE SYNCHRONIZED WITH AliITStrackV2::GetModuleIndexInfo()! //---------------------------------------------------------------------- if(fITSModule[ilayer]==-1) { idet = -1; status=0; xloc=-99.; zloc=-99.; return kFALSE; } Int_t module = fITSModule[ilayer]; idet = Int_t(module/1000000); module -= idet*1000000; status = Int_t(module/100000); module -= status*100000; Int_t signs = Int_t(module/10000); module-=signs*10000; Int_t xInt = Int_t(module/100); module -= xInt*100; Int_t zInt = module; if(signs==1) { xInt*=1; zInt*=1; } if(signs==2) { xInt*=1; zInt*=-1; } if(signs==3) { xInt*=-1; zInt*=1; } if(signs==4) { xInt*=-1; zInt*=-1; } xloc = 0.1*(Float_t)xInt; zloc = 0.1*(Float_t)zInt; if(status==4) idet = -1; return kTRUE; } //_______________________________________________________________________ UShort_t AliESDtrack::GetTPCclusters(Int_t *idx) const { //--------------------------------------------------------------------- // This function returns indices of the assgined ITS clusters //--------------------------------------------------------------------- if (idx && fFriendTrack) { Int_t *index=fFriendTrack->GetTPCindices(); if (index){ for (Int_t i=0; i -1) return fCacheNCrossedRows; fCacheNCrossedRows = GetTPCClusterInfo(2, 1); return fCacheNCrossedRows; } //_______________________________________________________________________ Float_t AliESDtrack::GetTPCClusterInfo(Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1, Int_t bitType ) const { // // TPC cluster information // type 0: get fraction of found/findable clusters with neighbourhood definition // 1: findable clusters with neighbourhood definition // 2: found clusters // bitType: // 0 - all cluster used // 1 - clusters used for the kalman update // definition of findable clusters: // a cluster is defined as findable if there is another cluster // within +- nNeighbours pad rows. The idea is to overcome threshold // effects with a very simple algorithm. // Int_t found=0; Int_t findable=0; Int_t last=-nNeighbours; const TBits & clusterMap = (bitType%2==0) ? fTPCClusterMap : fTPCFitMap; Int_t upperBound=clusterMap.GetNbits(); if (upperBound>row1) upperBound=row1; for (Int_t i=row0; i0) fraction=(Float_t)found/(Float_t)findable; else fraction=0; return fraction; } return 0; // undefined type - default value } //_______________________________________________________________________ Float_t AliESDtrack::GetTPCClusterDensity(Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1, Int_t bitType ) const { // // TPC cluster density - only rows where signal before and after given row are used // - slower function // type 0: get fraction of found/findable clusters with neighbourhood definition // 1: findable clusters with neighbourhood definition // 2: found clusters // bitType: // 0 - all cluster used // 1 - clusters used for the kalman update // definition of findable clusters: // a cluster is defined as findable if there is another cluster // within +- nNeighbours pad rows. The idea is to overcome threshold // effects with a very simple algorithm. // Int_t found=0; Int_t findable=0; // Int_t last=-nNeighbours; const TBits & clusterMap = (bitType%2==0) ? fTPCClusterMap : fTPCFitMap; Int_t upperBound=clusterMap.GetNbits(); if (upperBound>row1) upperBound=row1; for (Int_t i=row0; i=0 && clusterMap[i-idelta]) isDown=kTRUE; if (i+idelta0) fraction=(Float_t)found/(Float_t)findable; else fraction=0; return fraction; } return 0; // undefined type - default value } //_______________________________________________________________________ Double_t AliESDtrack::GetTPCdensity(Int_t row0, Int_t row1) const{ // // GetDensity of the clusters on given region between row0 and row1 // Dead zone effect takin into acoount // if (!fFriendTrack) return 0.0; Int_t good = 0; Int_t found = 0; // Int_t *index=fFriendTrack->GetTPCindices(); 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>TMath::Max((row1-row0)*0.5,0.0)) 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(); if (index) { for (Int_t i=0; i= AliESDtrack::kTRDnPlanes // 2. The idx array store not only the index but also the layer of the tracklet. // Therefore tracks with TRD gaps contain default values for indices [-1] if (!fFriendTrack) return 0; if (!idx) return GetTRDntracklets(); Int_t *index=fFriendTrack->GetTRDindices(); Int_t n = 0; for (Int_t i=0; i=0) n++; idx[i]=index[i]; } else idx[i] = -2; } return n; } //_______________________________________________________________________ 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)) { AliWarning(Form("Request for TRD plane[%d] outside range.", plane)); return -1.; } Int_t idx = fTRDnSlices-(kTRDnPlanes<<1)+plane; // Protection for backward compatibility if(idx<(GetNumberOfTRDslices()*kTRDnPlanes)) return -1.; if(sp) (*sp) = fTRDslices[idx+kTRDnPlanes]; return fTRDslices[idx]; } //____________________________________________________ Double_t AliESDtrack::GetTRDslice(Int_t plane, Int_t slice) const { //Gets the charge from the slice of the plane if(!fTRDslices) { //AliError("No TRD slices allocated for this track !"); return -1.; } if ((plane<0) || (plane>=kTRDnPlanes)) { AliError("Info for TRD plane not available !"); return -1.; } Int_t ns=GetNumberOfTRDslices(); if ((slice<-1) || (slice>=ns)) { //AliError("Wrong TRD slice !"); return -1.; } if(slice>=0) return fTRDslices[plane*ns + slice]; // return average of the dEdx measurements Double_t q=0.; Double32_t *s = &fTRDslices[plane*ns]; for (Int_t i=0; i0.) q+=(*s); return q/ns; } //____________________________________________________ void AliESDtrack::SetNumberOfTRDslices(Int_t n) { //Sets the number of slices used for PID if (fTRDnSlices) return; fTRDnSlices=n; fTRDslices=new Double32_t[fTRDnSlices]; // set-up correctly the allocated memory memset(fTRDslices, 0, n*sizeof(Double32_t)); for (Int_t i=GetNumberOfTRDslices(); i--;) fTRDslices[i]=-1.; } //____________________________________________________ void AliESDtrack::SetTRDslice(Double_t q, Int_t plane, Int_t slice) { //Sets the charge q in the slice of the plane if(!fTRDslices) { AliError("No TRD slices allocated for this track !"); return; } if ((plane<0) || (plane>=kTRDnPlanes)) { AliError("Info for TRD plane not allocated !"); return; } Int_t ns=GetNumberOfTRDslices(); if ((slice<0) || (slice>=ns)) { AliError(Form("Wrong TRD slice %d/%d, NSlices=%d",plane,slice,ns)); return; } Int_t n=plane*ns + slice; fTRDslices[n]=q; } //____________________________________________________ void AliESDtrack::SetTRDmomentum(Double_t p, Int_t plane, Double_t *sp) { if(!fTRDslices) { AliError("No TRD slices allocated for this track !"); return; } if ((plane<0) || (plane>=kTRDnPlanes)) { AliError("Info for TRD plane not allocated !"); return; } Int_t idx = fTRDnSlices-(kTRDnPlanes<<1)+plane; // Protection for backward compatibility if(idxPropagateToDCA(vtx, b, maxd, dz, cov)) return kFALSE; fdTPC = dz[0]; fzTPC = dz[1]; fCddTPC = cov[0]; fCdzTPC = cov[1]; fCzzTPC = cov[2]; Double_t covar[6]; vtx->GetCovMatrix(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>kVeryBig) return kFALSE; fCchi2TPC=chi2; if (!cParam) return kTRUE; *cParam = *fTPCInner; if (!cParam->Update(p,c)) return kFALSE; return kTRUE; } //_______________________________________________________________________ Bool_t AliESDtrack::RelateToVertexTPCBxByBz(const AliESDVertex *vtx, Double_t b[3], Double_t maxd, AliExternalTrackParam *cParam) { // // Try to relate the TPC-only track parameters to the vertex "vtx", // if the (rough) transverse impact parameter is not bigger then "maxd". // // All three components of the magnetic field ,"b[3]" (kG), // are taken into account. // // a) The TPC-only paramters are extapolated to the DCA to the vertex. // b) The impact parameters and their covariance matrix are calculated. // c) An attempt to constrain the TPC-only params to the vertex is done. // The constrained params are returned via "cParam". // // In the case of success, the returned value is kTRUE // otherwise, it's kFALSE) // if (!fTPCInner) return kFALSE; if (!vtx) return kFALSE; Double_t dz[2],cov[3]; if (!fTPCInner->PropagateToDCABxByBz(vtx, b, maxd, dz, cov)) return kFALSE; fdTPC = dz[0]; fzTPC = dz[1]; fCddTPC = cov[0]; fCdzTPC = cov[1]; fCzzTPC = cov[2]; Double_t covar[6]; vtx->GetCovMatrix(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>kVeryBig) return kFALSE; fCchi2TPC=chi2; if (!cParam) return kTRUE; *cParam = *fTPCInner; if (!cParam->Update(p,c)) return kFALSE; return kTRUE; } //_______________________________________________________________________ Bool_t AliESDtrack::RelateToVertex(const AliESDVertex *vtx, Double_t b, Double_t maxd, AliExternalTrackParam *cParam) { // // Try to relate this track to the vertex "vtx", // if the (rough) transverse impact parameter is not bigger then "maxd". // Magnetic field is "b" (kG). // // a) The track gets extapolated to the DCA to the vertex. // b) The impact parameters and their covariance matrix are calculated. // c) An attempt to constrain this track to the vertex is done. // The constrained params are returned via "cParam". // // In the case of success, the returned value is kTRUE // (otherwise, it's kFALSE) // if (!vtx) return kFALSE; Double_t dz[2],cov[3]; if (!PropagateToDCA(vtx, b, maxd, dz, cov)) return kFALSE; fD = dz[0]; fZ = dz[1]; fCdd = cov[0]; fCdz = cov[1]; fCzz = cov[2]; Double_t covar[6]; vtx->GetCovMatrix(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>kVeryBig) return kFALSE; fCchi2=chi2; //--- Could now these lines be removed ? --- delete fCp; fCp=new AliExternalTrackParam(*this); if (!fCp->Update(p,c)) {delete fCp; fCp=0; return kFALSE;} //---------------------------------------- fVertexID = vtx->GetID(); if (!cParam) return kTRUE; *cParam = *this; if (!cParam->Update(p,c)) return kFALSE; return kTRUE; } //_______________________________________________________________________ Bool_t AliESDtrack::RelateToVertexBxByBz(const AliESDVertex *vtx, Double_t b[3], Double_t maxd, AliExternalTrackParam *cParam) { // // Try to relate this track to the vertex "vtx", // if the (rough) transverse impact parameter is not bigger then "maxd". // Magnetic field is "b" (kG). // // a) The track gets extapolated to the DCA to the vertex. // b) The impact parameters and their covariance matrix are calculated. // c) An attempt to constrain this track to the vertex is done. // The constrained params are returned via "cParam". // // In the case of success, the returned value is kTRUE // (otherwise, it's kFALSE) // if (!vtx) return kFALSE; Double_t dz[2],cov[3]; if (!PropagateToDCABxByBz(vtx, b, maxd, dz, cov)) return kFALSE; fD = dz[0]; fZ = dz[1]; fCdd = cov[0]; fCdz = cov[1]; fCzz = cov[2]; Double_t covar[6]; vtx->GetCovMatrix(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>kVeryBig) return kFALSE; fCchi2=chi2; //--- Could now these lines be removed ? --- delete fCp; fCp=new AliExternalTrackParam(*this); if (!fCp->Update(p,c)) {delete fCp; fCp=0; return kFALSE;} //---------------------------------------- fVertexID = vtx->GetID(); if (!cParam) return kTRUE; *cParam = *this; if (!cParam->Update(p,c)) return kFALSE; return kTRUE; } //_______________________________________________________________________ void AliESDtrack::Print(Option_t *) const { // Prints info on the track AliExternalTrackParam::Print(); printf("ESD track info\n") ; Double_t p[AliPID::kSPECIES] ; 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()) ; printf("\n NchamberdEdx = %d\n", GetTRDNchamberdEdx()) ; printf("\n NclusterdEdx = %d\n", GetTRDNclusterdEdx()) ; } 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)); if (fHMPIDp) arrayRef.AddLast(new AliExternalTrackParam(*fHMPIDp)); // 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++; } } } //_______________________________________________________________________ void AliESDtrack::SetITSdEdxSamples(const Double_t s[4]) { // // Store the dE/dx samples measured by the two SSD and two SDD layers. // These samples are corrected for the track segment length. // for (Int_t i=0; i<4; i++) fITSdEdxSamples[i]=s[i]; } //_______________________________________________________________________ void AliESDtrack::GetITSdEdxSamples(Double_t s[4]) const { // // Get the dE/dx samples measured by the two SSD and two SDD layers. // These samples are corrected for the track segment length. // for (Int_t i=0; i<4; i++) s[i]=fITSdEdxSamples[i]; } UShort_t AliESDtrack::GetTPCnclsS(Int_t i0,Int_t i1) const{ // // get number of shared TPC clusters // return fTPCSharedMap.CountBits(i0)-fTPCSharedMap.CountBits(i1); } UShort_t AliESDtrack::GetTPCncls(Int_t i0,Int_t i1) const{ // // get number of TPC clusters // return fTPCClusterMap.CountBits(i0)-fTPCClusterMap.CountBits(i1); } //____________________________________________________________________ Double_t AliESDtrack::GetChi2TPCConstrainedVsGlobal(const AliESDVertex* vtx) const { // Calculates the chi2 between the TPC track (TPCinner) constrained to the primary vertex and the global track // // Returns -1 in case the calculation failed // // Value is cached as a non-persistent member. // // Code adapted from original code by GSI group (Jacek, Marian, Michael) // cache, ignoring that a different vertex might be passed if (fCacheChi2TPCConstrainedVsGlobalVertex == vtx) return fCacheChi2TPCConstrainedVsGlobal; fCacheChi2TPCConstrainedVsGlobal = -1; fCacheChi2TPCConstrainedVsGlobalVertex = vtx; Double_t x[3]; GetXYZ(x); Double_t b[3]; AliTrackerBase::GetBxByBz(x,b); if (!fTPCInner) { AliWarning("Could not get TPC Inner Param."); return fCacheChi2TPCConstrainedVsGlobal; } // clone for constraining AliExternalTrackParam* tpcInnerC = new AliExternalTrackParam(*fTPCInner); if (!tpcInnerC) { AliWarning("Clone of TPCInnerParam failed."); return fCacheChi2TPCConstrainedVsGlobal; } // transform to the track reference frame Bool_t isOK = tpcInnerC->Rotate(GetAlpha()); isOK &= tpcInnerC->PropagateTo(GetX(), b[2]); if (!isOK) { delete tpcInnerC; tpcInnerC = 0; AliWarning("Rotation/Propagation of track failed.") ; return fCacheChi2TPCConstrainedVsGlobal; } // constrain TPCinner isOK = tpcInnerC->ConstrainToVertex(vtx, b); // transform to the track reference frame isOK &= tpcInnerC->Rotate(GetAlpha()); isOK &= tpcInnerC->PropagateTo(GetX(), b[2]); if (!isOK) { AliWarning("ConstrainTPCInner failed.") ; delete tpcInnerC; tpcInnerC = 0; return fCacheChi2TPCConstrainedVsGlobal; } // calculate chi2 between vi and vj vectors // with covi and covj covariance matrices // chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj) TMatrixD deltaT(5,1); TMatrixD delta(1,5); TMatrixD covarM(5,5); for (Int_t ipar=0; ipar<5; ipar++) { deltaT(ipar,0) = tpcInnerC->GetParameter()[ipar] - GetParameter()[ipar]; delta(0,ipar) = tpcInnerC->GetParameter()[ipar] - GetParameter()[ipar]; for (Int_t jpar=0; jpar<5; jpar++) { Int_t index = GetIndex(ipar,jpar); covarM(ipar,jpar) = GetCovariance()[index]+tpcInnerC->GetCovariance()[index]; } } // chi2 distance TPC constrained and TPC+ITS TMatrixD covarMInv = covarM.Invert(); TMatrixD mat2 = covarMInv*deltaT; TMatrixD chi2 = delta*mat2; delete tpcInnerC; tpcInnerC = 0; fCacheChi2TPCConstrainedVsGlobal = chi2(0,0); return fCacheChi2TPCConstrainedVsGlobal; } void AliESDtrack::SetDetectorPID(const AliDetectorPID *pid) { // // Set the detector PID // if (fDetectorPID) delete fDetectorPID; fDetectorPID=pid; } Double_t AliESDtrack::GetLengthInActiveZone( Int_t mode, Double_t deltaY, Double_t deltaZ, Double_t bz, Double_t exbPhi , TTreeSRedirector * pcstream){ // // Input parameters: // mode - type of external track parameters // deltaY - user defined "dead region" in cm // deltaZ - user defined "active region" in cm (250 cm drift lenght - 14 cm L1 delay // bz - magnetic field // exbPhi - optional rotation due to the ExB effect // return value: // the length of the track in cm in "active volume" of the TPC // if (mode==0) return GetLengthInActiveZone(this, deltaY,deltaZ,bz, exbPhi,pcstream); if (mode==1) return GetLengthInActiveZone(fIp, deltaY,deltaZ,bz, exbPhi,pcstream); if (mode==2) return GetLengthInActiveZone(fOp, deltaY,deltaZ,bz, exbPhi,pcstream); return 0; } Double_t AliESDtrack::GetLengthInActiveZone( AliExternalTrackParam *paramT, Double_t deltaY, Double_t deltaZ, Double_t bz, Double_t exbPhi , TTreeSRedirector * pcstream){ // // Numerical code to calculate the length of the track in active region of the TPC // ( can be speed up if somebody wants to invest time - analysical version shoult be possible) // // Input parameters: // paramT - external track parameters // deltaY - user defined "dead region" in cm // deltaZ - user defined "active region" in cm (250 cm drift lenght - 14 cm L1 delay // bz - magnetic field // exbPhi - optional rotation due to the ExB effect // return value: // the length of the track in cm in "active volume" of the TPC // const Double_t rIn=85; const Double_t rOut=245; Double_t xyz[3], pxyz[3]; if (paramT->GetXYZAt(rIn,bz,xyz)){ paramT->GetPxPyPzAt(rIn,bz,pxyz); }else{ paramT->GetXYZ(xyz); paramT->GetPxPyPz(pxyz); } // Double_t dca = -paramT->GetD(0,0,bz); // get impact parameter distance to point (0,0) Double_t radius= TMath::Abs(1/paramT->GetC(bz)); // Double_t sign = paramT->GetSign(); Double_t R0 = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); // radius at current point Double_t phiR0 = TMath::ATan2(xyz[1],xyz[0]); // angle of given point Double_t dPhiR0= -TMath::ASin((dca*dca-2*dca*radius*sign+R0*R0)/(2*R0*(dca-radius*sign))); Double_t phi0 = phiR0-(dPhiR0); // global phi offset to be added // // AliExternalTrackParam paramR=(*paramT); Double_t length=0; for (Double_t R=rIn; R<=rOut; R++){ Double_t sinPhi=(dca*dca-2*dca*radius*sign+R*R)/(2*R*(dca-radius*sign)); if (TMath::Abs(sinPhi)>=1) continue; Double_t dphi = -TMath::ASin(sinPhi); Double_t phi = phi0+dphi; // global phi Int_t sector = TMath::Nint(9*phi/(TMath::Pi())); Double_t dPhiEdge = phi-(sector*TMath::Pi()/9)+exbPhi; // distance to sector boundary in rphi Double_t dX = R*TMath::Cos(phi)-xyz[0]; Double_t dY = R*TMath::Sin(phi)-xyz[1]; Double_t deltaPhi = 2*TMath::ASin(0.5*TMath::Sqrt(dX*dX+dY*dY)/radius); Double_t z = xyz[2]+deltaPhi*radius*paramT->GetTgl(); if (TMath::Abs(dPhiEdge*R)>deltaY && TMath::Abs(z)