init histogram arrays when they read for the Terminate execution
[u/mrichter/AliRoot.git] / STEER / AliTPCpidESD.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 //-----------------------------------------------------------------
17 //           Implementation of the TPC PID class
18 // Very naive one... Should be made better by the detector experts...
19 //      Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
20 // With many additions and modifications suggested by
21 //      Alexander Kalweit, GSI, alexander.philipp.kalweit@cern.ch
22 //      Dariusz Miskowiec, GSI, D.Miskowiec@gsi.de
23 //-----------------------------------------------------------------
24
25 #include "AliTPCpidESD.h"
26 #include "AliESDEvent.h"
27 #include "AliESDtrack.h"
28 #include "AliMathBase.h"
29
30 ClassImp(AliTPCpidESD)
31
32 //_________________________________________________________________________
33 AliTPCpidESD::AliTPCpidESD():
34     fMIP(50.),
35     fRes(0.07),
36     fRange(5.),
37     fKp1(0.76176e-1),
38     fKp2(10.632),
39     fKp3(0.13279e-4),
40     fKp4(1.8631),
41     fKp5(1.9479)
42 {
43   //
44   //  The default constructor
45   //
46 }
47
48 //_________________________________________________________________________
49 AliTPCpidESD::AliTPCpidESD(Double_t *param):
50     fMIP(param[0]),
51     fRes(param[1]),
52     fRange(param[2]),
53     fKp1(0.76176e-1),
54     fKp2(10.632),
55     fKp3(0.13279e-4),
56     fKp4(1.8631),
57     fKp5(1.9479)
58 {
59   //
60   //  The main constructor
61   //
62 }
63
64 Double_t AliTPCpidESD::Bethe(Double_t betaGamma) const {
65   //
66   // This is the Bethe-Bloch function normalised to 1 at the minimum
67   // WARNING
68   // Simulated and reconstructed Bethe-Bloch differs
69   //           Simulated  curve is the dNprim/dx
70   //           Reconstructed is proportianal dNtot/dx
71   // Temporary fix for production -  Simple linear correction function
72   // Future    2 Bethe Bloch formulas needed
73   //           1. for simulation
74   //           2. for reconstructed PID
75   //
76   const Float_t kmeanCorrection =0.1;
77   Double_t bb=AliMathBase::BetheBlochAleph(betaGamma,fKp1,fKp2,fKp3,fKp4,fKp5);
78   Double_t meanCorrection =(1+(bb-1)*kmeanCorrection);
79   bb *= meanCorrection;
80   return bb;
81 }
82
83 //_________________________________________________________________________
84 void AliTPCpidESD::SetBetheBlochParameters(Double_t kp1,
85                              Double_t kp2,
86                              Double_t kp3,
87                              Double_t kp4,
88                              Double_t kp5) {
89   //
90   // Set the parameters of the ALEPH Bethe-Bloch formula
91   //
92   fKp1=kp1;
93   fKp2=kp2;
94   fKp3=kp3;
95   fKp4=kp4;
96   fKp5=kp5;
97 }
98
99 //_________________________________________________________________________
100 Bool_t AliTPCpidESD::ExpectedSigmas(const AliESDtrack *t,
101                                      Double_t s[],
102                                      Int_t n) const {
103   //
104   // Calculate the expected dE/dx resolution as the function of 
105   // the information stored in the track.
106   //
107   // At the moment, this resolution is just proportional to the expected 
108   // signal. This can be improved. By taking into account the number of
109   // assigned clusters, for example.  
110   //
111   Bool_t ok=kFALSE;
112   Double_t signals[AliPID::kSPECIESN]; 
113   if (ExpectedSignals(t,signals,n)) {
114      for (Int_t i=0; i<n; i++) s[i] = fRes*signals[i];
115      ok=kTRUE;
116   }
117   return ok;
118 }
119
120 //_________________________________________________________________________
121 Bool_t AliTPCpidESD::NumberOfSigmas(const AliESDtrack *t,
122                                      Double_t s[],
123                                      Int_t n) const {
124   //
125   // Calculate the deviation of the actual PID signal from the expected
126   // signal, in units of expected sigmas.
127   //
128   Bool_t ok=kFALSE;
129
130   Double_t dedx=t->GetTPCsignal()/fMIP;
131   Double_t sigmas[AliPID::kSPECIESN];
132   if (ExpectedSigmas(t,sigmas,n)) {
133      Double_t signals[AliPID::kSPECIESN];
134      if (ExpectedSignals(t,signals,n)) {
135        for (Int_t i=0; i<n; i++) s[i] = (signals[i] - dedx)/sigmas[i];
136        ok=kTRUE;
137      }
138   }
139   return ok;
140 }
141
142 //_________________________________________________________________________
143 Bool_t AliTPCpidESD::ExpectedSignals(const AliESDtrack *t, 
144                                       Double_t s[],
145                                       Int_t n) const {
146   //
147   // Calculates the expected PID signals as the function of 
148   // the information stored in the track.
149   //  
150   // At the moment, these signals are just the results of calling the 
151   // Bethe-Bloch formula. 
152   // This can be improved. By taking into account the number of
153   // assigned clusters and/or the track dip angle, for example.  
154   //
155   
156   Double_t mom=t->GetP();
157   const AliExternalTrackParam *in=t->GetInnerParam();
158   if (in) mom=in->GetP();
159
160   for (Int_t i=0; i<n; i++) {
161       Double_t mass=AliPID::ParticleMass(i); 
162       s[i]=Bethe(mom/mass);
163   }
164
165   return kTRUE;
166 }
167
168 //_________________________________________________________________________
169 Double_t AliTPCpidESD::GetExpectedSignal(const AliESDtrack *t,
170                                          AliPID::EParticleType n) const {
171   //
172   // Calculates the expected PID signal as the function of 
173   // the information stored in the track, for the specified particle type 
174   //  
175   // At the moment, these signals are just the results of calling the 
176   // Bethe-Bloch formula. 
177   // This can be improved. By taking into account the number of
178   // assigned clusters and/or the track dip angle, for example.  
179   //
180   
181   Double_t mom=t->GetP();
182   const AliExternalTrackParam *in=t->GetInnerParam();
183   if (in) mom=in->GetP();
184
185   Double_t mass=AliPID::ParticleMass(n); 
186   return Bethe(mom/mass);
187 }
188
189 //_________________________________________________________________________
190 Double_t AliTPCpidESD::GetExpectedSigma(const AliESDtrack *t,
191                                          AliPID::EParticleType n) const {
192   //
193   // Calculates the expected sigma of the PID signal as the function of 
194   // the information stored in the track, for the specified particle type 
195   //  
196   //
197   // At the moment, this sigma is just proportional to the expected 
198   // signal. This can be improved. By taking into account the number of
199   // assigned clusters, for example.  
200   //
201   
202   return fRes*GetExpectedSignal(t,n);
203 }
204
205 //_________________________________________________________________________
206 Double_t AliTPCpidESD::GetNumberOfSigmas(const AliESDtrack *t,
207                                          AliPID::EParticleType n) const {
208   // 
209   // Calculate the deviation of the actual PID signal from the expected
210   // signal, in units of expected sigmas, for the specified particle type
211   //
212   
213     Double_t dedx=t->GetTPCsignal()/fMIP;
214     return (dedx - GetExpectedSignal(t,n))/GetExpectedSigma(t,n);
215 }
216
217 //_________________________________________________________________________
218 Int_t AliTPCpidESD::MakePID(AliESDEvent *event)
219 {
220   //
221   //  This function calculates the "detector response" PID probabilities 
222   //
223   Int_t ntrk=event->GetNumberOfTracks();
224   for (Int_t i=0; i<ntrk; i++) {
225     AliESDtrack *t=event->GetTrack(i);
226     if ((t->GetStatus()&AliESDtrack::kTPCin )==0)
227       if ((t->GetStatus()&AliESDtrack::kTPCout)==0) continue;
228     Double_t p[10];
229     Double_t dedx=t->GetTPCsignal()/fMIP;
230     Bool_t mismatch=kTRUE, heavy=kTRUE;
231     for (Int_t j=0; j<AliPID::kSPECIES; j++) {
232       AliPID::EParticleType type=AliPID::EParticleType(j);
233       Double_t bethe=GetExpectedSignal(t,type); 
234       Double_t sigma=GetExpectedSigma(t,type);
235       if (TMath::Abs(dedx-bethe) > fRange*sigma) {
236         p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
237       } else {
238         p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
239         mismatch=kFALSE;
240       }
241
242       // Check for particles heavier than (AliPID::kSPECIES - 1)
243       if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
244
245     }
246
247     if (mismatch)
248        for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1/AliPID::kSPECIES;
249
250     t->SetTPCpid(p);
251
252     if (heavy) t->ResetStatus(AliESDtrack::kTPCpid);
253
254   }
255   return 0;
256 }