Implemented additional functionality according to the

https://savannah.cern.ch/bugs/?47411

(J.Belikov)

diff --git a/TPC/AliTPCpidESD.cxx b/TPC/AliTPCpidESD.cxx
index 66dc38e95389582eead12da7904b23bea859ef50..a3dd6d501dcfa2082662d1e179d287b1ba475322 100644 (file)
--- a/TPC/AliTPCpidESD.cxx
+++ b/TPC/AliTPCpidESD.cxx
@@ -17,6 +17,9 @@
 //           Implementation of the TPC PID class
 // Very naive one... Should be made better by the detector experts...
 //      Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
+// With many additions and modifications suggested by
+//      Alexander Kalweit, GSI, alexander.philipp.kalweit@cern.ch
+//      Dariusz Miskowiec, GSI, D.Miskowiec@gsi.de
 //-----------------------------------------------------------------
 
 #include "AliTPCpidESD.h"
@@ -27,24 +30,42 @@
 ClassImp(AliTPCpidESD)
 
 //_________________________________________________________________________
-  AliTPCpidESD::AliTPCpidESD(Double_t *param):
-    fMIP(0.),
-    fRes(0.),
-    fRange(0.)
+AliTPCpidESD::AliTPCpidESD():
+    fMIP(50.),
+    fRes(0.07),
+    fRange(5.),
+    fKp1(0.76176e-1),
+    fKp2(10.632),
+    fKp3(0.13279e-4),
+    fKp4(1.8631),
+    fKp5(1.9479)
+{
+  //
+  //  The default constructor
+  //
+}
+
+//_________________________________________________________________________
+AliTPCpidESD::AliTPCpidESD(Double_t *param):
+    fMIP(param[0]),
+    fRes(param[1]),
+    fRange(param[2]),
+    fKp1(0.76176e-1),
+    fKp2(10.632),
+    fKp3(0.13279e-4),
+    fKp4(1.8631),
+    fKp5(1.9479)
 {
   //
   //  The main constructor
   //
-  fMIP=param[0];
-  fRes=param[1];
-  fRange=param[2];
 }
 
-Double_t AliTPCpidESD::Bethe(Double_t betaGamma) {
+Double_t AliTPCpidESD::Bethe(Double_t betaGamma) const {
   //
   // This is the Bethe-Bloch function normalised to 1 at the minimum
   // WARNING
-  // Simulated and reconstructed Bethe-Bloch differs 
+  // Simulated and reconstructed Bethe-Bloch differs
   //           Simulated  curve is the dNprim/dx
   //           Reconstructed is proportianal dNtot/dx
   // Temporary fix for production -  Simple linear correction function
@@ -53,12 +74,146 @@ Double_t AliTPCpidESD::Bethe(Double_t betaGamma) {
   //           2. for reconstructed PID
   //
   const Float_t kmeanCorrection =0.1;
-  Double_t bb = AliMathBase::BetheBlochAleph(betaGamma);
+  Double_t bb=AliMathBase::BetheBlochAleph(betaGamma,fKp1,fKp2,fKp3,fKp4,fKp5);
   Double_t meanCorrection =(1+(bb-1)*kmeanCorrection);
   bb *= meanCorrection;
   return bb;
 }
 
+//_________________________________________________________________________
+void AliTPCpidESD::SetBetheBlochParameters(Double_t kp1,
+                             Double_t kp2,
+                             Double_t kp3,
+                             Double_t kp4,
+                             Double_t kp5) {
+  //
+  // Set the parameters of the ALEPH Bethe-Bloch formula
+  //
+  fKp1=kp1;
+  fKp2=kp2;
+  fKp3=kp3;
+  fKp4=kp4;
+  fKp5=kp5;
+}
+
+//_________________________________________________________________________
+Bool_t AliTPCpidESD::ExpectedSigmas(const AliESDtrack *t,
+                                     Double_t s[],
+                                     Int_t n) const {
+  //
+  // Calculate the expected dE/dx resolution as the function of 
+  // the information stored in the track.
+  //
+  // At the moment, this resolution is just proportional to the expected 
+  // signal. This can be improved. By taking into account the number of
+  // assigned clusters, for example.  
+  //
+  Bool_t ok=kFALSE;
+  Double_t signals[AliPID::kSPECIESN]; 
+  if (ExpectedSignals(t,signals,n)) {
+     for (Int_t i=0; i<n; i++) s[i] = fRes*signals[i];
+     ok=kTRUE;
+  }
+  return ok;
+}
+
+//_________________________________________________________________________
+Bool_t AliTPCpidESD::NumberOfSigmas(const AliESDtrack *t,
+                                     Double_t s[],
+                                     Int_t n) const {
+  //
+  // Calculate the deviation of the actual PID signal from the expected
+  // signal, in units of expected sigmas.
+  //
+  Bool_t ok=kFALSE;
+
+  Double_t dedx=t->GetTPCsignal()/fMIP;
+  Double_t sigmas[AliPID::kSPECIESN];
+  if (ExpectedSigmas(t,sigmas,n)) {
+     Double_t signals[AliPID::kSPECIESN];
+     if (ExpectedSignals(t,signals,n)) {
+       for (Int_t i=0; i<n; i++) s[i] = (signals[i] - dedx)/sigmas[i];
+       ok=kTRUE;
+     }
+  }
+  return ok;
+}
+
+//_________________________________________________________________________
+Bool_t AliTPCpidESD::ExpectedSignals(const AliESDtrack *t, 
+                                      Double_t s[],
+                                      Int_t n) const {
+  //
+  // Calculates the expected PID signals as the function of 
+  // the information stored in the track.
+  //  
+  // At the moment, these signals are just the results of calling the 
+  // Bethe-Bloch formula. 
+  // This can be improved. By taking into account the number of
+  // assigned clusters and/or the track dip angle, for example.  
+  //
+  
+  Double_t mom=t->GetP();
+  const AliExternalTrackParam *in=t->GetInnerParam();
+  if (in) mom=in->GetP();
+
+  for (Int_t i=0; i<n; i++) {
+      Double_t mass=AliPID::ParticleMass(i); 
+      s[i]=Bethe(mom/mass);
+  }
+
+  return kTRUE;
+}
+
+//_________________________________________________________________________
+Double_t AliTPCpidESD::GetExpectedSignal(const AliESDtrack *t,
+                                         AliPID::EParticleType n) const {
+  //
+  // Calculates the expected PID signal as the function of 
+  // the information stored in the track, for the specified particle type 
+  //  
+  // At the moment, these signals are just the results of calling the 
+  // Bethe-Bloch formula. 
+  // This can be improved. By taking into account the number of
+  // assigned clusters and/or the track dip angle, for example.  
+  //
+  
+  Double_t mom=t->GetP();
+  const AliExternalTrackParam *in=t->GetInnerParam();
+  if (in) mom=in->GetP();
+
+  Double_t mass=AliPID::ParticleMass(n); 
+  return Bethe(mom/mass);
+}
+
+//_________________________________________________________________________
+Double_t AliTPCpidESD::GetExpectedSigma(const AliESDtrack *t,
+                                         AliPID::EParticleType n) const {
+  //
+  // Calculates the expected sigma of the PID signal as the function of 
+  // the information stored in the track, for the specified particle type 
+  //  
+  //
+  // At the moment, this sigma is just proportional to the expected 
+  // signal. This can be improved. By taking into account the number of
+  // assigned clusters, for example.  
+  //
+  
+  return fRes*GetExpectedSignal(t,n);
+}
+
+//_________________________________________________________________________
+Double_t AliTPCpidESD::GetNumberOfSigmas(const AliESDtrack *t,
+                                         AliPID::EParticleType n) const {
+  // 
+  // Calculate the deviation of the actual PID signal from the expected
+  // signal, in units of expected sigmas, for the specified particle type
+  //
+  
+    Double_t dedx=t->GetTPCsignal()/fMIP;
+    return (dedx - GetExpectedSignal(t,n))/GetExpectedSigma(t,n);
+}
+
 //_________________________________________________________________________
 Int_t AliTPCpidESD::MakePID(AliESDEvent *event)
 {
@@ -71,15 +226,12 @@ Int_t AliTPCpidESD::MakePID(AliESDEvent *event)
     if ((t->GetStatus()&AliESDtrack::kTPCin )==0)
       if ((t->GetStatus()&AliESDtrack::kTPCout)==0) continue;
     Double_t p[10];
-    Double_t mom=t->GetP();
-    const AliExternalTrackParam *in=t->GetInnerParam();
-    if (in) mom=in->GetP();
     Double_t dedx=t->GetTPCsignal()/fMIP;
     Bool_t mismatch=kTRUE, heavy=kTRUE;
     for (Int_t j=0; j<AliPID::kSPECIES; j++) {
-      Double_t mass=AliPID::ParticleMass(j);
-      Double_t bethe=Bethe(mom/mass); 
-      Double_t sigma=fRes*bethe;
+      AliPID::EParticleType type=AliPID::EParticleType(j);
+      Double_t bethe=GetExpectedSignal(t,type); 
+      Double_t sigma=GetExpectedSigma(t,type);
       if (TMath::Abs(dedx-bethe) > fRange*sigma) {
        p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
       } else {

diff --git a/TPC/AliTPCpidESD.h b/TPC/AliTPCpidESD.h
index 593e1589f806c5ee6857ab4fe084aee940959dcf..b9ca7e9bdb87bf4918d5de583792b07c77d614ca 100644 (file)
--- a/TPC/AliTPCpidESD.h
+++ b/TPC/AliTPCpidESD.h
@@ -7,23 +7,60 @@
 //                    TPC PID class
 // A very naive design... Should be made better by the detector experts...
 //   Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch 
+// With many additions and modifications suggested by
+//      Alexander Kalweit, GSI, alexander.philipp.kalweit@cern.ch
+//      Dariusz Miskowiec, GSI, D.Miskowiec@gsi.de
 //-------------------------------------------------------
 #include <Rtypes.h>
 
+#include "AliPID.h"
+
 class AliESDEvent;
+class AliESDtrack;
 
 class AliTPCpidESD {
 public:
-  AliTPCpidESD():fMIP(0.),fRes(0.),fRange(0.){}
+  AliTPCpidESD();
   AliTPCpidESD(Double_t *param);
   virtual ~AliTPCpidESD() {}
+  void SetBetheBlochParameters(Double_t kp1,
+                               Double_t kp2,
+                               Double_t kp3,
+                               Double_t kp4,
+                               Double_t kp5
+                               );
   Int_t MakePID(AliESDEvent *event);
-  static Double_t Bethe(Double_t bg);
+  Double_t Bethe(Double_t bg) const;
+
+  Bool_t ExpectedSignals(const AliESDtrack *t, 
+                          Double_t s[], 
+                          Int_t n=AliPID::kSPECIES) const;
+  Bool_t ExpectedSigmas(const AliESDtrack *t, 
+                         Double_t s[],
+                         Int_t n=AliPID::kSPECIES) const;
+  Bool_t NumberOfSigmas(const AliESDtrack *t, 
+                         Double_t s[],
+                         Int_t n=AliPID::kSPECIES) const;
+
+  Double_t GetExpectedSignal(const AliESDtrack *t,
+                     AliPID::EParticleType n=AliPID::kKaon) const;
+  Double_t GetExpectedSigma(const AliESDtrack *t,
+                     AliPID::EParticleType n=AliPID::kKaon) const;
+  Double_t GetNumberOfSigmas(const AliESDtrack *t,
+                     AliPID::EParticleType n=AliPID::kKaon) const;
+
 private:
   Double_t fMIP;          // dEdx for MIP
   Double_t fRes;          // relative dEdx resolution
   Double_t fRange;        // one particle type PID range (in sigmas)
-  ClassDef(AliTPCpidESD,1)   // TPC PID class
+
+  Double_t fKp1;   // Parameters
+  Double_t fKp2;   //    of
+  Double_t fKp3;   // the ALEPH
+  Double_t fKp4;   // Bethe-Bloch
+  Double_t fKp5;   // formula
+
+  ClassDef(AliTPCpidESD,2)   // TPC PID class
 };
 
 #endif