Number of sigma pedestal cut increased to 4
[u/mrichter/AliRoot.git] / STEER / AliPID.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 /* $Id$ */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 //                                                                           //
20 // particle id probability densities                                         //
21 //                                                                           //
22 // The AliPID class stores the probability densities for the different       //
23 // particle type hypotheses electron, muon, pion, kaon, proton, photon,      //
24 // pi0, neutron, K0 and electron conversion. These probability densities     //
25 // are determined from the detector response functions.                      //
26 // The * and *= operators are overloaded for AliPID to combine the PIDs      //
27 // from different detectors.                                                 //
28 //                                                                           //
29 // The Bayesian probability to be a particle of a given type can be          //
30 // calculated from the probability densities, if the a priori probabilities  //
31 // (or abundences, concentrations) of particle species are known. These      //
32 // priors can be given as argument to the GetProbability or GetMostProbable  //
33 // method or they can be set globally by calling the static method           //
34 // SetPriors().                                                              //
35 //                                                                           //
36 // The implementation of this class is based on the note ...                 //
37 // by Iouri Belikov and Karel Safarik.                                       //
38 //                                                                           //
39 ///////////////////////////////////////////////////////////////////////////////
40
41 #include <TClass.h>
42 #include <TDatabasePDG.h>
43 #include <TPDGCode.h>
44
45 #include "AliLog.h"
46 #include "AliPDG.h"
47 #include "AliPID.h"
48
49 #define M(PID) TDatabasePDG::Instance()->GetParticle(fgkParticleCode[(PID)])->Mass()
50
51 ClassImp(AliPID)
52
53 const char* AliPID::fgkParticleName[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
54   "electron",
55   "muon",
56   "pion",
57   "kaon",
58   "proton",
59   "photon",
60   "pi0",
61   "neutron",
62   "kaon0",
63   "eleCon",
64   "deuteron",
65   "triton",
66   "helium-3",
67   "alpha",
68   "unknown"
69 };
70
71 const char* AliPID::fgkParticleShortName[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
72   "e",
73   "mu",
74   "pi",
75   "K",
76   "p",
77   "photon",
78   "pi0",
79   "n",
80   "K0",
81   "eleCon",
82   "d",
83   "t",
84   "he3",
85   "alpha",
86   "unknown"
87 };
88
89 const char* AliPID::fgkParticleLatexName[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
90   "e",
91   "#mu",
92   "#pi",
93   "K",
94   "p",
95   "#gamma",
96   "#pi_{0}",
97   "n",
98   "K_{0}",
99   "eleCon",
100   "d",
101   "t",
102   "he3",
103   "#alpha",
104   "unknown"
105 };
106
107 const Int_t AliPID::fgkParticleCode[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
108   ::kElectron, 
109   ::kMuonMinus, 
110   ::kPiPlus, 
111   ::kKPlus, 
112   ::kProton,
113   ::kGamma,
114   ::kPi0,
115   ::kNeutron,
116   ::kK0,
117   ::kElectron,
118   1000010020,
119   1000010030,
120   1000020030,
121   1000020040,
122   0
123 };
124
125 /*const*/ Float_t AliPID::fgkParticleMass[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
126   0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
127   /*
128   M(kElectron),  // electron
129   M(kMuon), // muon
130   M(kPion),    // pion
131   M(kKaon),     // kaon
132   M(kProton),    // proton
133   M(kPhoton),     // photon
134   M(kPi0),       // pi0
135   M(kNeutron),   // neutron
136   M(kKaon0),        // kaon0
137   M(kEleCon),     // electron conversion
138   M(kDeuteron), // deuteron
139   M(kTriton),   // triton
140   M(kHe3),      // he3
141   M(kAlpha),    // alpha
142   0.00000        // unknown
143   */
144 };
145
146 Double_t AliPID::fgPrior[kSPECIESN] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
147
148
149 //_______________________________________________________________________
150 AliPID::AliPID() :
151   TObject(),
152   fCharged(0)
153 {
154   //
155   // Default constructor
156   //
157   Init();
158   // set default values (= equal probabilities)
159   for (Int_t i = 0; i < kSPECIESN; i++)
160     fProbDensity[i] = 1./kSPECIESN;
161 }
162
163 //_______________________________________________________________________
164 AliPID::AliPID(const Double_t* probDensity, Bool_t charged) : 
165   TObject(),
166   fCharged(charged)
167 {
168   //
169   // Standard constructor
170   //
171   Init();
172   // set given probability densities
173   for (Int_t i = 0; i < kSPECIES; i++) 
174     fProbDensity[i] = probDensity[i];
175
176   for (Int_t i = kSPECIES; i < kSPECIESN; i++) 
177     fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
178 }
179
180 //_______________________________________________________________________
181 AliPID::AliPID(const Float_t* probDensity, Bool_t charged) :
182   TObject(),
183   fCharged(charged)
184 {
185   //
186   // Standard constructor
187   //
188   Init();
189   // set given probability densities
190   for (Int_t i = 0; i < kSPECIES; i++) 
191     fProbDensity[i] = probDensity[i];
192
193   for (Int_t i = kSPECIES; i < kSPECIESN; i++) 
194     fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
195 }
196
197 //_______________________________________________________________________
198 AliPID::AliPID(const AliPID& pid) : 
199   TObject(pid),
200   fCharged(pid.fCharged)
201 {
202   //
203   // copy constructor
204   //
205   // We do not call init here, MUST already be done
206   for (Int_t i = 0; i < kSPECIESN; i++) 
207     fProbDensity[i] = pid.fProbDensity[i];
208 }
209
210 //_______________________________________________________________________
211 void AliPID::SetProbabilities(const Double_t* probDensity, Bool_t charged) 
212 {
213   //
214   // Set the probability densities
215   //
216   for (Int_t i = 0; i < kSPECIES; i++) 
217     fProbDensity[i] = probDensity[i];
218
219   for (Int_t i = kSPECIES; i < kSPECIESN; i++) 
220     fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
221 }
222
223 //_______________________________________________________________________
224 AliPID& AliPID::operator = (const AliPID& pid)
225 {
226 // assignment operator
227
228   fCharged = pid.fCharged;
229   for (Int_t i = 0; i < kSPECIESN; i++) {
230     fProbDensity[i] = pid.fProbDensity[i];
231   }
232   return *this;
233 }
234
235 //_______________________________________________________________________
236 void AliPID::Init() 
237 {
238   //
239   // Initialise the masses
240   //
241   // Initialise only once... 
242   if(!fgkParticleMass[0]) {
243     AliPDG::AddParticlesToPdgDataBase();
244     for (Int_t i = 0; i < kSPECIESN + kSPECIESLN; i++) 
245       fgkParticleMass[i] = M(i);
246   }
247 }
248
249 //_____________________________________________________________________________
250 Double_t AliPID::GetProbability(EParticleType iType,
251                                 const Double_t* prior) const
252 {
253   //
254   // Get the probability to be a particle of type "iType"
255   // assuming the a priori probabilities "prior"
256   //
257   Double_t sum = 0.;
258   Int_t nSpecies = ((fCharged) ? kSPECIES : kSPECIESN);
259   for (Int_t i = 0; i < nSpecies; i++) {
260     sum += fProbDensity[i] * prior[i];
261   }
262   if (sum <= 0) {
263     AliError("Invalid probability densities or priors");
264     return -1;
265   }
266   return fProbDensity[iType] * prior[iType] / sum;
267 }
268
269 //_____________________________________________________________________________
270 Double_t AliPID::GetProbability(EParticleType iType) const
271 {
272 // get the probability to be a particle of type "iType"
273 // assuming the globaly set a priori probabilities
274
275   return GetProbability(iType, fgPrior);
276 }
277
278 //_____________________________________________________________________________
279 void AliPID::GetProbabilities(Double_t* probabilities,
280                               const Double_t* prior) const
281 {
282 // get the probabilities to be a particle of given type
283 // assuming the a priori probabilities "prior"
284
285   Double_t sum = 0.;
286   Int_t nSpecies = ((fCharged) ? kSPECIES : kSPECIESN);
287   for (Int_t i = 0; i < nSpecies; i++) {
288     sum += fProbDensity[i] * prior[i];
289   }
290   if (sum <= 0) {
291     AliError("Invalid probability densities or priors");
292     for (Int_t i = 0; i < nSpecies; i++) probabilities[i] = -1;
293     return;
294   }
295   for (Int_t i = 0; i < nSpecies; i++) {
296     probabilities[i] = fProbDensity[i] * prior[i] / sum;
297   }
298 }
299
300 //_____________________________________________________________________________
301 void AliPID::GetProbabilities(Double_t* probabilities) const
302 {
303 // get the probabilities to be a particle of given type
304 // assuming the globaly set a priori probabilities
305
306   GetProbabilities(probabilities, fgPrior);
307 }
308
309 //_____________________________________________________________________________
310 AliPID::EParticleType AliPID::GetMostProbable(const Double_t* prior) const
311 {
312 // get the most probable particle id hypothesis
313 // assuming the a priori probabilities "prior"
314
315   Double_t max = 0.;
316   EParticleType id = kPion;
317   Int_t nSpecies = ((fCharged) ? kSPECIES : kSPECIESN);
318   for (Int_t i = 0; i < nSpecies; i++) {
319     Double_t prob = fProbDensity[i] * prior[i];
320     if (prob > max) {
321       max = prob;
322       id = EParticleType(i);
323     }
324   }
325   if (max == 0) {
326     AliError("Invalid probability densities or priors");
327   }
328   return id;
329 }
330
331 //_____________________________________________________________________________
332 AliPID::EParticleType AliPID::GetMostProbable() const
333 {
334 // get the most probable particle id hypothesis
335 // assuming the globaly set a priori probabilities
336
337   return GetMostProbable(fgPrior);
338 }
339
340
341 //_____________________________________________________________________________
342 void AliPID::SetPriors(const Double_t* prior, Bool_t charged)
343 {
344 // use the given priors as global a priori probabilities
345
346   Double_t sum = 0;
347   for (Int_t i = 0; i < kSPECIESN; i++) {
348     if (charged && (i >= kSPECIES)) {
349       fgPrior[i] = 0;      
350     } else {
351       if (prior[i] < 0) {
352         AliWarningClass(Form("negative prior (%g) for %ss. "
353                              "Using 0 instead.", prior[i], 
354                              fgkParticleName[i]));
355         fgPrior[i] = 0;
356       } else {
357         fgPrior[i] = prior[i];
358       }
359     }
360     sum += prior[i];
361   }
362   if (sum == 0) {
363     AliWarningClass("all priors are zero.");
364   }
365 }
366
367 //_____________________________________________________________________________
368 void AliPID::SetPrior(EParticleType iType, Double_t prior)
369 {
370 // use the given prior as global a priori probability for particles
371 // of type "iType"
372
373   if (prior < 0) {
374     AliWarningClass(Form("negative prior (%g) for %ss. Using 0 instead.", 
375                          prior, fgkParticleName[iType]));
376     prior = 0;
377   }
378   fgPrior[iType] = prior;
379 }
380
381
382 //_____________________________________________________________________________
383 AliPID& AliPID::operator *= (const AliPID& pid)
384 {
385 // combine this probability densities with the one of "pid"
386
387   for (Int_t i = 0; i < kSPECIESN; i++) {
388     fProbDensity[i] *= pid.fProbDensity[i];
389   }
390   return *this;
391 }
392
393 //_____________________________________________________________________________
394 AliPID operator * (const AliPID& pid1, const AliPID& pid2)
395 {
396 // combine the two probability densities
397
398   AliPID result;
399   result *= pid1;
400   result *= pid2;
401   return result;
402 }