Update of combined PID (T.Kuhr)

diff --git a/STEER/AliESDtest.C b/STEER/AliESDtest.C
index 5029952501f5ac60b20c779ef85e5f1b4992bfd5..fc0bf081bc0ae4e993d27607b03c8040bf89ac5a 100644 (file)
--- a/STEER/AliESDtest.C
+++ b/STEER/AliESDtest.C
@@ -76,16 +76,8 @@ Int_t AliESDtest(Int_t nev=1,
    AliTRDtracker trdTracker(trdcf);
 
    //An instance of the TRD PID maker
-   TFile* pidFile = TFile::Open("pid.root");
-   if (!pidFile->IsOpen()) {
-     cerr << "Can't get pid.root !\n";
-     return 7;
-   }
-   AliTRDPartID* trdPID = (AliTRDPartID*) pidFile->Get("AliTRDPartID");
-   if (!trdPID) {
-     cerr << "Can't get PID object !\n";
-     return 8;
-   }
+   AliTRDPartID* trdPID = AliTRDPartID::GetFromFile();
+   if (!trdPID) return 7;
 
    TFile *ef=TFile::Open("AliESDs.root","RECREATE");
    if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
@@ -145,7 +137,6 @@ Int_t AliESDtest(Int_t nev=1,
    }
    timer.Stop(); timer.Print();
 
-   pidFile->Close();
    trdcf->Close();
    delete geom;
    itscf->Close();

diff --git a/STEER/pid.root b/STEER/pid.root
index 6f3383841b59151814968acff4854164b05aa534..5a6c35adc1f35ecaab70996c8240b03305b15ca1 100644 (file)
Binary files a/STEER/pid.root and b/STEER/pid.root differ

diff --git a/TRD/AliTRDPartID.cxx b/TRD/AliTRDPartID.cxx
index 91e3ca334ccf404318c45e054dc4ab27ec12350c..3a095143daae92eedc358edf5831aa682d35b632 100644 (file)
--- a/TRD/AliTRDPartID.cxx
+++ b/TRD/AliTRDPartID.cxx
@@ -15,10 +15,12 @@
 
 #include "AliTRDPartID.h"
 #include "AliESDtrack.h"
-#include "TProfile.h"
-#include "TF1.h"
-
+#include <TProfile.h>
+#include <TF1.h>
+#include <TFile.h>
+#include <TROOT.h>
 #include <TMath.h>
+#include <Riostream.h>
 
 ClassImp(AliTRDPartID)
 
@@ -51,9 +53,8 @@ Bool_t AliTRDPartID::MakePID(AliESDtrack* track)
 {
 // This function calculates the "detector response" PID probabilities 
 
-  static const Double_t masses[]={
-    0.000511, 0.105658, 0.139570, 0.493677, 0.938272
-  };
+  static const Double_t masses[AliESDtrack::kSPECIES] = 
+    {0.000511, 0.105658, 0.139570, 0.493677, 0.938272};
 
   if (((track->GetStatus()&AliESDtrack::kTRDin) == 0) &&
       ((track->GetStatus()&AliESDtrack::kTRDout) == 0)) return kFALSE;
@@ -69,7 +70,9 @@ Bool_t AliTRDPartID::MakePID(AliESDtrack* track)
     if (TMath::Abs(measuredSignal - expectedSignal) > fRange * expectedError) {
       track->SetTRDpid(iSpecies, 0.);
     } else {
-      Double_t p = TMath::Gaus(measuredSignal, expectedSignal, expectedError);
+      Double_t delta = (measuredSignal-expectedSignal) / expectedError;
+      const Double_t kInvSqr2Pi = 0.398942280401432703;
+      Double_t p = kInvSqr2Pi / expectedError * TMath::Exp(-delta*delta / 2.);
       pSum += p;
       track->SetTRDpid(iSpecies, p);
     }
@@ -90,15 +93,14 @@ void AliTRDPartID::FitBetheBloch(TProfile* dEdxVsBetaGamma)
 // determine the parameters of the bethe bloch function
 
   if (fBetheBloch) delete fBetheBloch;
-  fBetheBloch = new TF1("fBetheBlochTRD", fcnBetheBloch, 
-                       0.001, 100000., 5);
+  fBetheBloch = new TF1("fBetheBlochTRD", fcnBetheBloch, 0.001, 100000., 5);
   fBetheBloch->SetParameters(1, 10, 0.00002, 2, 2);
   fBetheBloch->SetParLimits(2, 0., 0.01);
   fBetheBloch->SetParLimits(3, 0., 10.);
   fBetheBloch->SetParLimits(4, 0., 10.);
   fBetheBloch->SetFillStyle(0);
   fBetheBloch->SetLineColor(kRed);
-  dEdxVsBetaGamma->Fit(fBetheBloch, "N");
+  dEdxVsBetaGamma->Fit(fBetheBloch, "NIR", "goff", 0.6, dEdxVsBetaGamma->GetXaxis()->GetXmax());
 }
 
 TF1* AliTRDPartID::CreateBetheBloch(Double_t mass)
@@ -118,6 +120,38 @@ TF1* AliTRDPartID::CreateBetheBloch(Double_t mass)
   return result;
 }
 
+AliTRDPartID* AliTRDPartID::GetFromFile(const char* fileName)
+{
+// read an AliTRDPartID object from a file
+
+  TFile* pidFile = (TFile*) gROOT->GetListOfFiles()->FindObject(fileName);
+  Bool_t fileOpened = kFALSE;
+  if (!pidFile) {
+    pidFile = TFile::Open(fileName);
+    fileOpened = kTRUE;
+  }
+  if (!pidFile->IsOpen()) {
+    cerr << "Can't open " << fileName << " !\n";
+    if (fileOpened) delete pidFile;
+    return NULL;
+  }
+  gROOT->cd();
+
+  AliTRDPartID* trdPID = (AliTRDPartID*) pidFile->Get("AliTRDPartID")->Clone();
+  if (!trdPID) {
+    cerr << "Can't get PID object !\n";
+  } else {
+    trdPID->GetBetheBloch()->SetFunction(fcnBetheBloch);
+  }
+
+  if (fileOpened) {
+    pidFile->Close();
+    delete pidFile;
+  }
+
+  return trdPID;
+}
+
 Double_t AliTRDPartID::fcnBetheBloch(Double_t* xx, Double_t* par)
 {
 // parametrized bethe bloch function (xx[0] = beta*gamma = p/m):

diff --git a/TRD/AliTRDPartID.h b/TRD/AliTRDPartID.h
index 40af6035646f1b611a45293f457dd9150b2657d8..d0e2ac70c5b9a1677e9fa1fd8a8201687cdb2ee2 100644 (file)
--- a/TRD/AliTRDPartID.h
+++ b/TRD/AliTRDPartID.h
@@ -22,6 +22,8 @@ class AliTRDPartID: public TObject {
     inline TF1*     GetBetheBloch() {return fBetheBloch;};
     TF1*            CreateBetheBloch(Double_t mass);
 
+    static AliTRDPartID* GetFromFile(const char* fileName = "pid.root");
+
   private:
     static Double_t fcnBetheBloch(Double_t* xx, Double_t* par);
     static Double_t fcnBetheBlochMass(Double_t* xx, Double_t* par);