- Reshuffling of the particle codes in AliPID. Now the light nuclei are between the
[u/mrichter/AliRoot.git] / STEER / ESD / AliESDtrack.cxx
index f6c75e5..68a77bd 100644 (file)
 #include "TPolyMarker3D.h"
 #include "AliTrackerBase.h"
 #include "AliTPCdEdxInfo.h"
+#include "AliDetectorPID.h"
 
 ClassImp(AliESDtrack)
 
@@ -200,6 +201,7 @@ AliESDtrack::AliESDtrack() :
   fGlobalChi2(0),
   fITSsignal(0),
   fTPCsignal(0),
+  fTPCsignalTuned(0),
   fTPCsignalS(0),
   fTPCdEdxInfo(0),
   fTRDsignal(0),
@@ -229,13 +231,18 @@ AliESDtrack::AliESDtrack() :
   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)
+  fCacheChi2TPCConstrainedVsGlobalVertex(0),
+  fDetectorPID(0x0),
+  fTrackPhiOnEMCal(-999),
+  fTrackEtaOnEMCal(-999)
 {
   //
   // The default ESD constructor 
@@ -310,6 +317,7 @@ AliESDtrack::AliESDtrack(const AliESDtrack& track):
   fGlobalChi2(track.fGlobalChi2),
   fITSsignal(track.fITSsignal),
   fTPCsignal(track.fTPCsignal),
+  fTPCsignalTuned(track.fTPCsignalTuned),
   fTPCsignalS(track.fTPCsignalS),
   fTPCdEdxInfo(0),
   fTRDsignal(track.fTRDsignal),
@@ -339,13 +347,18 @@ AliESDtrack::AliESDtrack(const AliESDtrack& track):
   fTRDncls(track.fTRDncls),
   fTRDncls0(track.fTRDncls0),
   fTRDntracklets(track.fTRDntracklets),
+  fTRDNchamberdEdx(track.fTRDNchamberdEdx),
+  fTRDNclusterdEdx(track.fTRDNclusterdEdx),
   fTRDnSlices(track.fTRDnSlices),
   fTRDslices(0x0),
   fVertexID(track.fVertexID),
   fESDEvent(track.fESDEvent),
   fCacheNCrossedRows(track.fCacheNCrossedRows),
   fCacheChi2TPCConstrainedVsGlobal(track.fCacheChi2TPCConstrainedVsGlobal),
-  fCacheChi2TPCConstrainedVsGlobalVertex(track.fCacheChi2TPCConstrainedVsGlobalVertex)
+  fCacheChi2TPCConstrainedVsGlobalVertex(track.fCacheChi2TPCConstrainedVsGlobalVertex),
+  fDetectorPID(0x0),
+  fTrackPhiOnEMCal(track.fTrackPhiOnEMCal),
+  fTrackEtaOnEMCal(track.fTrackEtaOnEMCal)
 {
   //
   //copy constructor
@@ -371,6 +384,7 @@ AliESDtrack::AliESDtrack(const AliESDtrack& track):
     for (Int_t i=0; i<fTRDnSlices; i++) fTRDslices[i]=track.fTRDslices[i];
   }
 
+  if (track.fDetectorPID) fDetectorPID = new AliDetectorPID(*track.fDetectorPID);
 
   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTRDr[i]=track.fTRDr[i]; 
   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTOFr[i]=track.fTOFr[i];
@@ -432,6 +446,7 @@ AliESDtrack::AliESDtrack(const AliVTrack *track) :
   fGlobalChi2(0),
   fITSsignal(0),
   fTPCsignal(0),
+  fTPCsignalTuned(0),
   fTPCsignalS(0),
   fTPCdEdxInfo(0),
   fTRDsignal(0),
@@ -461,13 +476,18 @@ AliESDtrack::AliESDtrack(const AliVTrack *track) :
   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)
+  fCacheChi2TPCConstrainedVsGlobalVertex(0),
+  fDetectorPID(0x0),
+  fTrackPhiOnEMCal(-999),
+  fTrackEtaOnEMCal(-999)
 {
   //
   // ESD track from AliVTrack.
@@ -574,6 +594,7 @@ AliESDtrack::AliESDtrack(TParticle * part) :
   fGlobalChi2(0),
   fITSsignal(0),
   fTPCsignal(0),
+  fTPCsignalTuned(0),
   fTPCsignalS(0),
   fTPCdEdxInfo(0),
   fTRDsignal(0),
@@ -603,13 +624,18 @@ AliESDtrack::AliESDtrack(TParticle * part) :
   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)  
+  fCacheChi2TPCConstrainedVsGlobalVertex(0),
+  fDetectorPID(0x0),
+  fTrackPhiOnEMCal(-999),
+  fTrackEtaOnEMCal(-999)
 {
   //
   // ESD track from TParticle
@@ -739,11 +765,17 @@ AliESDtrack::~AliESDtrack(){
   delete fOp;
   delete fHMPIDp;
   delete fCp; 
-  if (fFriendTrack) delete fFriendTrack;
-  if (fTPCdEdxInfo) delete fTPCdEdxInfo;
-  fFriendTrack=NULL;
+  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){
@@ -760,7 +792,7 @@ AliESDtrack &AliESDtrack::operator=(const AliESDtrack &source){
   }
   else{
     // no track param delete the old one
-    if(fCp)delete fCp;
+    delete fCp;
     fCp = 0;
   }
 
@@ -771,7 +803,7 @@ AliESDtrack &AliESDtrack::operator=(const AliESDtrack &source){
   }
   else{
     // no track param delete the old one
-    if(fIp)delete fIp;
+    delete fIp;
     fIp = 0;
   }
 
@@ -783,7 +815,7 @@ AliESDtrack &AliESDtrack::operator=(const AliESDtrack &source){
   }
   else{
     // no track param delete the old one
-    if(fTPCInner)delete fTPCInner;
+    delete fTPCInner;
     fTPCInner = 0;
   }
 
@@ -799,7 +831,7 @@ AliESDtrack &AliESDtrack::operator=(const AliESDtrack &source){
   }
   else{
     // no track param delete the old one
-    if(fOp)delete fOp;
+    delete fOp;
     fOp = 0;
   }
 
@@ -811,7 +843,7 @@ AliESDtrack &AliESDtrack::operator=(const AliESDtrack &source){
   }
   else{
     // no track param delete the old one
-    if(fHMPIDp)delete fHMPIDp;
+    delete fHMPIDp;
     fHMPIDp = 0;
   }
 
@@ -895,11 +927,14 @@ AliESDtrack &AliESDtrack::operator=(const AliESDtrack &source){
   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];   
@@ -946,6 +981,18 @@ AliESDtrack &AliESDtrack::operator=(const AliESDtrack &source){
   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;
+
   return *this;
 }
 
@@ -1034,6 +1081,7 @@ Bool_t AliESDtrack::FillTPCOnlyTrack(AliESDtrack &track){
   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];
 
@@ -1113,6 +1161,7 @@ void AliESDtrack::MakeMiniESDtrack(){
   fTPCClusterMap = 0;  
   fTPCSharedMap = 0;  
   fTPCsignal= 0;      
+  fTPCsignalTuned= 0;
   fTPCsignalS= 0;      
   fTPCsignalN= 0;      
   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTPCr[i]=0; 
@@ -1126,6 +1175,9 @@ void AliESDtrack::MakeMiniESDtrack(){
   fTRDncls = 0;       
   fTRDncls0 = 0;       
   fTRDsignal = 0;      
+  fTRDNchamberdEdx = 0;
+  fTRDNclusterdEdx = 0;
+
   for (Int_t i=0;i<kTRDnPlanes;i++) {
     fTRDTimBin[i]  = 0;
   }
@@ -1177,24 +1229,32 @@ void AliESDtrack::MakeMiniESDtrack(){
 } 
 
 //_______________________________________________________________________
-Int_t AliESDtrack::GetPID() const 
+Int_t AliESDtrack::GetPID(Bool_t tpcOnly) const 
 {
   // Returns the particle most probable id
   Int_t i;
-  for (i=0; i<AliPID::kSPECIES-1; i++) if (fR[i] != fR[i+1]) break;
-  //
-  if (i == AliPID::kSPECIES-1) return AliPID::kPion;  // If all the probabilities are equal, return the pion mass
+  const Double32_t *prob = 0;
+  if (tpcOnly) { // check if TPCpid is valid
+    prob = fTPCr;
+    for (i=0; i<AliPID::kSPECIES-1; i++) if (prob[i] != prob[i+1]) break;
+    if (i == AliPID::kSPECIES-1) prob = 0; // not valid, try with combined pid
+  }
+  if (!prob) { // either requested TPCpid is not valid or comb.pid is requested 
+    prob = fR;
+    for (i=0; i<AliPID::kSPECIES-1; i++) if (prob[i] != prob[i+1]) break;
+    if (i == AliPID::kSPECIES-1) return AliPID::kPion;  // If all the probabilities are equal, return the pion mass
+  }
   //
   Float_t max=0.;
   Int_t k=-1;
-  for (i=0; i<AliPID::kSPECIES; i++) if (fR[i]>max) {k=i; max=fR[i];}
+  for (i=0; i<AliPID::kSPECIES; i++) if (prob[i]>max) {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 (fR[0]<fR[3]*10.) return AliPID::kKaon;
+      if (prob[0]<prob[3]*10.) return AliPID::kKaon;
     if ((p>0.75)&&(p<0.85))
-      if (fR[0]<fR[4]*10.) return AliPID::kProton;
+      if (prob[0]<prob[4]*10.) return AliPID::kProton;
     return AliPID::kElectron;
   }
   if (k==1) return AliPID::kMuon; 
@@ -1206,7 +1266,7 @@ Int_t AliESDtrack::GetPID() const
 }
 
 //_______________________________________________________________________
-Int_t AliESDtrack::GetTOFBunchCrossing(Double_t b) const 
+Int_t AliESDtrack::GetTOFBunchCrossing(Double_t b, Bool_t pidTPConly) const 
 {
   // Returns the number of bunch crossings after trigger (assuming 25ns spacing)
   const double kSpacing = 25e3; // min interbanch spacing
@@ -1216,7 +1276,7 @@ Int_t AliESDtrack::GetTOFBunchCrossing(Double_t b) const
   //
   double tdif = fTOFsignal;
   if (IsOn(kTIME)) { // integrated time info is there
-    int pid = GetPID();
+    int pid = GetPID(pidTPConly);
     tdif -= fTrackTime[pid];
   }
   else { // assume integrated time info from TOF radius and momentum
@@ -1224,7 +1284,7 @@ Int_t AliESDtrack::GetTOFBunchCrossing(Double_t b) const
     const double kCSpeed = 3.e-2; // cm/ps
     double p = GetP();
     if (p<0.01) return bcid;
-    double m = GetMass();
+    double m = GetMass(pidTPConly);
     double curv = GetC(b);
     double path = TMath::Abs(curv)>kAlmost0 ? // account for curvature
       2./curv*TMath::ASin(kRTOF*curv/2.)*TMath::Sqrt(1.+GetTgl()*GetTgl()) : kRTOF;
@@ -1384,7 +1444,8 @@ Bool_t AliESDtrack::UpdateTrackParams(const AliKalmanTrack *t, ULong_t flags){
       delete [] indexTRD;
     }    
     
-    fTRDsignal=t->GetPIDsignal();
+    //commented out by Xianguo
+    //fTRDsignal=t->GetPIDsignal();
     }
     break;
   case kTRDbackup:
@@ -1755,31 +1816,33 @@ Float_t AliESDtrack::GetTPCCrossedRows() const
 }
 
 //_______________________________________________________________________
-Float_t AliESDtrack::GetTPCClusterInfo(Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1) const
+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.
   //
 
-  if (type==2) return fTPCClusterMap.CountBits();
   
   Int_t found=0;
   Int_t findable=0;
   Int_t last=-nNeighbours;
+  const TBits & clusterMap = (bitType%2==0) ? fTPCClusterMap : fTPCFitMap;
   
-  Int_t upperBound=fTPCClusterMap.GetNbits();
+  Int_t upperBound=clusterMap.GetNbits();
   if (upperBound>row1) upperBound=row1;
   for (Int_t i=row0; i<upperBound; ++i){
     //look to current row
-    if (fTPCClusterMap[i]) {
+    if (clusterMap[i]) {
       last=i;
       ++found;
       ++findable;
@@ -1792,12 +1855,62 @@ Float_t AliESDtrack::GetTPCClusterInfo(Int_t nNeighbours/*=3*/, Int_t type/*=0*/
     }
     //look to nNeighbours after
     for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
-      if (fTPCClusterMap[j]){
+      if (clusterMap[j]){
         ++findable;
         break;
       }
     }
   }
+  if (type==2) return found;
+  if (type==1) return findable;
+  
+  if (type==0){
+    Float_t fraction=0;
+    if (findable>0) 
+      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<upperBound; ++i){
+    Bool_t isUp=kFALSE;
+    Bool_t isDown=kFALSE;
+    for (Int_t idelta=1; idelta<=nNeighbours; idelta++){
+      if (i-idelta>=0 && clusterMap[i-idelta]) isDown=kTRUE;
+      if (i+idelta<upperBound && clusterMap[i+idelta]) isUp=kTRUE;
+    }
+    if (isUp&&isDown){
+      ++findable;
+      if (clusterMap[i]) ++found;
+    }
+  }
+  if (type==2) return found;
   if (type==1) return findable;
   
   if (type==0){
@@ -1811,6 +1924,9 @@ Float_t AliESDtrack::GetTPCClusterInfo(Int_t nNeighbours/*=3*/, Int_t type/*=0*/
   return 0;  // undefined type - default value
 }
 
+
+
+
 //_______________________________________________________________________
 Double_t AliESDtrack::GetTPCdensity(Int_t row0, Int_t row1) const{
   //
@@ -1828,7 +1944,7 @@ Double_t AliESDtrack::GetTPCdensity(Int_t row0, Int_t row1) const{
     if (idx>0)    found++;
   }
   Float_t density=0.5;
-  if (good>(row1-row0)*0.5) density = Float_t(found)/Float_t(good);
+  if (good>TMath::Max((row1-row0)*0.5,0.0)) density = Float_t(found)/Float_t(good);
   return density;
 }
 
@@ -2309,7 +2425,7 @@ void AliESDtrack::Print(Option_t *) const {
   // Prints info on the track
   AliExternalTrackParam::Print();
   printf("ESD track info\n") ; 
-  Double_t p[AliPID::kSPECIESN] ; 
+  Double_t p[AliPID::kSPECIES] ;
   Int_t index = 0 ; 
   if( IsOn(kITSpid) ){
     printf("From ITS: ") ; 
@@ -2331,6 +2447,8 @@ void AliESDtrack::Print(Option_t *) const {
     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: ") ; 
@@ -2525,3 +2643,13 @@ Double_t AliESDtrack::GetChi2TPCConstrainedVsGlobal(const AliESDVertex* vtx) con
   fCacheChi2TPCConstrainedVsGlobal = chi2(0,0);
   return fCacheChi2TPCConstrainedVsGlobal;
 }
+
+void AliESDtrack::SetDetectorPID(const AliDetectorPID *pid)
+{
+  //
+  // Set the detector PID
+  //
+  if (fDetectorPID) delete fDetectorPID;
+  fDetectorPID=pid;
+  
+}