Update of combined PID (T.Kuhr)
[u/mrichter/AliRoot.git] / TRD / AliTRDPartID.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 #include "AliTRDPartID.h"
17 #include "AliESDtrack.h"
18 #include <TProfile.h>
19 #include <TF1.h>
20 #include <TFile.h>
21 #include <TROOT.h>
22 #include <TMath.h>
23 #include <Riostream.h>
24
25 ClassImp(AliTRDPartID)
26
27
28 AliTRDPartID::AliTRDPartID()
29 {
30 // create a TRD particle identifier
31
32   fBetheBloch = NULL;
33   fRes = 0.2;
34   fRange = 3.;
35 }
36
37 AliTRDPartID::AliTRDPartID(TF1* betheBloch, Double_t res, Double_t range)
38 {
39 // create a TRD particle identifier with custom settings
40
41   fBetheBloch = betheBloch;
42   fRes = res;
43   fRange = range;
44 }
45
46 AliTRDPartID::~AliTRDPartID()
47 {
48   if (fBetheBloch) delete fBetheBloch;
49 }
50
51
52 Bool_t AliTRDPartID::MakePID(AliESDtrack* track)
53 {
54 // This function calculates the "detector response" PID probabilities 
55
56   static const Double_t masses[AliESDtrack::kSPECIES] = 
57     {0.000511, 0.105658, 0.139570, 0.493677, 0.938272};
58
59   if (((track->GetStatus()&AliESDtrack::kTRDin) == 0) &&
60       ((track->GetStatus()&AliESDtrack::kTRDout) == 0)) return kFALSE;
61   Double_t momentum = track->GetP();
62   if (momentum < 0.001) return kFALSE;
63
64   // get the probability densities
65   Double_t pSum = 0;
66   for (Int_t iSpecies = 0; iSpecies < AliESDtrack::kSPECIES; iSpecies++) {
67     Double_t expectedSignal = fBetheBloch->Eval(momentum/masses[iSpecies]);
68     Double_t expectedError = fRes * expectedSignal;
69     Double_t measuredSignal = track->GetTRDsignal();
70     if (TMath::Abs(measuredSignal - expectedSignal) > fRange * expectedError) {
71       track->SetTRDpid(iSpecies, 0.);
72     } else {
73       Double_t delta = (measuredSignal-expectedSignal) / expectedError;
74       const Double_t kInvSqr2Pi = 0.398942280401432703;
75       Double_t p = kInvSqr2Pi / expectedError * TMath::Exp(-delta*delta / 2.);
76       pSum += p;
77       track->SetTRDpid(iSpecies, p);
78     }
79   }
80
81   // "normalize" the probability densities
82   if (pSum <= 0) return kFALSE;
83   for (Int_t iSpecies = 0; iSpecies < AliESDtrack::kSPECIES; iSpecies++) {
84     track->SetTRDpid(iSpecies, track->GetTRDpid(iSpecies) / pSum);
85   }
86
87   return kTRUE;
88 }
89
90
91 void AliTRDPartID::FitBetheBloch(TProfile* dEdxVsBetaGamma)
92 {
93 // determine the parameters of the bethe bloch function
94
95   if (fBetheBloch) delete fBetheBloch;
96   fBetheBloch = new TF1("fBetheBlochTRD", fcnBetheBloch, 0.001, 100000., 5);
97   fBetheBloch->SetParameters(1, 10, 0.00002, 2, 2);
98   fBetheBloch->SetParLimits(2, 0., 0.01);
99   fBetheBloch->SetParLimits(3, 0., 10.);
100   fBetheBloch->SetParLimits(4, 0., 10.);
101   fBetheBloch->SetFillStyle(0);
102   fBetheBloch->SetLineColor(kRed);
103   dEdxVsBetaGamma->Fit(fBetheBloch, "NIR", "goff", 0.6, dEdxVsBetaGamma->GetXaxis()->GetXmax());
104 }
105
106 TF1* AliTRDPartID::CreateBetheBloch(Double_t mass)
107 {
108 // create a function for expected dE/dx vs p
109
110   TF1* result = new TF1("betheBlochMass", fcnBetheBlochMass, 
111                         0.001, 100000., 6);
112   result->SetParameter(0, mass);
113   if (fBetheBloch) {
114     for (Int_t iPar = 0; iPar < 5; iPar++) {
115       result->SetParameter(iPar+1, fBetheBloch->GetParameter(iPar));
116       result->SetParError(iPar+1, fBetheBloch->GetParError(iPar));
117     }
118   }
119   result->SetFillStyle(0);
120   return result;
121 }
122
123 AliTRDPartID* AliTRDPartID::GetFromFile(const char* fileName)
124 {
125 // read an AliTRDPartID object from a file
126
127   TFile* pidFile = (TFile*) gROOT->GetListOfFiles()->FindObject(fileName);
128   Bool_t fileOpened = kFALSE;
129   if (!pidFile) {
130     pidFile = TFile::Open(fileName);
131     fileOpened = kTRUE;
132   }
133   if (!pidFile->IsOpen()) {
134     cerr << "Can't open " << fileName << " !\n";
135     if (fileOpened) delete pidFile;
136     return NULL;
137   }
138   gROOT->cd();
139
140   AliTRDPartID* trdPID = (AliTRDPartID*) pidFile->Get("AliTRDPartID")->Clone();
141   if (!trdPID) {
142     cerr << "Can't get PID object !\n";
143   } else {
144     trdPID->GetBetheBloch()->SetFunction(fcnBetheBloch);
145   }
146
147   if (fileOpened) {
148     pidFile->Close();
149     delete pidFile;
150   }
151
152   return trdPID;
153 }
154
155 Double_t AliTRDPartID::fcnBetheBloch(Double_t* xx, Double_t* par)
156 {
157 // parametrized bethe bloch function (xx[0] = beta*gamma = p/m):
158 //
159 //   p0/beta^p3 * [ p1 - log(p2 + 1/(beta*gamma)^p4) - beta^p3 ]
160 //
161
162   Double_t betaGamma2 = xx[0] * xx[0];
163   Double_t beta2 = betaGamma2 / (1. + betaGamma2);
164   Double_t betaPar3 = TMath::Power(beta2, par[3]/2.);
165   return par[0]/betaPar3 * (par[1] - TMath::Log(TMath::Abs(par[2] + TMath::Power(betaGamma2, -par[4]/2.))) - betaPar3);
166 }
167
168 Double_t AliTRDPartID::fcnBetheBlochMass(Double_t* xx, Double_t* par)
169 {
170 // parametrized bethe bloch function:  xx[0] = p,  p0 = mass
171
172   Double_t betaGamma = xx[0] / par[0];
173   return fcnBetheBloch(&betaGamma, &par[1]);
174 }