7c5b2613a9ebd4e8a330ba6cf5624c7bb78d6268
[u/mrichter/AliRoot.git] / STEER / STEERBase / 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 /*const*/ Float_t AliPID::fgkParticleMassZ[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
146   0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
147   /*
148   M(kElectron),  // electron
149   M(kMuon), // muon
150   M(kPion),    // pion
151   M(kKaon),     // kaon
152   M(kProton),    // proton
153   M(kPhoton),     // photon
154   M(kPi0),       // pi0
155   M(kNeutron),   // neutron
156   M(kKaon0),        // kaon0
157   M(kEleCon),     // electron conversion
158   M(kDeuteron), // deuteron
159   M(kTriton),   // triton
160   M(kHe3)/2,      // he3
161   M(kAlpha)/2,    // alpha
162   0.00000        // unknown
163   */
164 };
165
166 Double_t AliPID::fgPrior[kSPECIESN] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
167
168
169 //_______________________________________________________________________
170 AliPID::AliPID() :
171   TObject(),
172   fCharged(0)
173 {
174   //
175   // Default constructor
176   //
177   Init();
178   // set default values (= equal probabilities)
179   for (Int_t i = 0; i < kSPECIESN; i++)
180     fProbDensity[i] = 1./kSPECIESN;
181 }
182
183 //_______________________________________________________________________
184 AliPID::AliPID(const Double_t* probDensity, Bool_t charged) : 
185   TObject(),
186   fCharged(charged)
187 {
188   //
189   // Standard constructor
190   //
191   Init();
192   // set given probability densities
193   for (Int_t i = 0; i < kSPECIES; i++) 
194     fProbDensity[i] = probDensity[i];
195
196   for (Int_t i = kSPECIES; i < kSPECIESN; i++) 
197     fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
198 }
199
200 //_______________________________________________________________________
201 AliPID::AliPID(const Float_t* probDensity, Bool_t charged) :
202   TObject(),
203   fCharged(charged)
204 {
205   //
206   // Standard constructor
207   //
208   Init();
209   // set given probability densities
210   for (Int_t i = 0; i < kSPECIES; i++) 
211     fProbDensity[i] = probDensity[i];
212
213   for (Int_t i = kSPECIES; i < kSPECIESN; i++) 
214     fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
215 }
216
217 //_______________________________________________________________________
218 AliPID::AliPID(const AliPID& pid) : 
219   TObject(pid),
220   fCharged(pid.fCharged)
221 {
222   //
223   // copy constructor
224   //
225   // We do not call init here, MUST already be done
226   for (Int_t i = 0; i < kSPECIESN; i++) 
227     fProbDensity[i] = pid.fProbDensity[i];
228 }
229
230 //_______________________________________________________________________
231 void AliPID::SetProbabilities(const Double_t* probDensity, Bool_t charged) 
232 {
233   //
234   // Set the probability densities
235   //
236   for (Int_t i = 0; i < kSPECIES; i++) 
237     fProbDensity[i] = probDensity[i];
238
239   for (Int_t i = kSPECIES; i < kSPECIESN; i++) 
240     fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
241 }
242
243 //_______________________________________________________________________
244 AliPID& AliPID::operator = (const AliPID& pid)
245 {
246 // assignment operator
247
248   if(this != &pid) {
249     fCharged = pid.fCharged;
250     for (Int_t i = 0; i < kSPECIESN; i++) {
251       fProbDensity[i] = pid.fProbDensity[i];
252     }
253   }
254   return *this;
255 }
256
257 //_______________________________________________________________________
258 void AliPID::Init() 
259 {
260   //
261   // Initialise the masses
262   //
263   // Initialise only once... 
264   if(!fgkParticleMass[0]) {
265     AliPDG::AddParticlesToPdgDataBase();
266     for (Int_t i = 0; i < kSPECIESN + kSPECIESLN; i++) {
267       fgkParticleMass[i] = M(i);
268       if (i == kHe3 || i == kAlpha) fgkParticleMassZ[i] = M(i)/2.;
269       else fgkParticleMassZ[i]=M(i);
270     }
271   }
272 }
273
274 //_____________________________________________________________________________
275 Double_t AliPID::GetProbability(EParticleType iType,
276                                 const Double_t* prior) const
277 {
278   //
279   // Get the probability to be a particle of type "iType"
280   // assuming the a priori probabilities "prior"
281   //
282   Double_t sum = 0.;
283   Int_t nSpecies = ((fCharged) ? kSPECIES : kSPECIESN);
284   for (Int_t i = 0; i < nSpecies; i++) {
285     sum += fProbDensity[i] * prior[i];
286   }
287   if (sum <= 0) {
288     AliError("Invalid probability densities or priors");
289     return -1;
290   }
291   return fProbDensity[iType] * prior[iType] / sum;
292 }
293
294 //_____________________________________________________________________________
295 Double_t AliPID::GetProbability(EParticleType iType) const
296 {
297 // get the probability to be a particle of type "iType"
298 // assuming the globaly set a priori probabilities
299
300   return GetProbability(iType, fgPrior);
301 }
302
303 //_____________________________________________________________________________
304 void AliPID::GetProbabilities(Double_t* probabilities,
305                               const Double_t* prior) const
306 {
307 // get the probabilities to be a particle of given type
308 // assuming the a priori probabilities "prior"
309
310   Double_t sum = 0.;
311   Int_t nSpecies = ((fCharged) ? kSPECIES : kSPECIESN);
312   for (Int_t i = 0; i < nSpecies; i++) {
313     sum += fProbDensity[i] * prior[i];
314   }
315   if (sum <= 0) {
316     AliError("Invalid probability densities or priors");
317     for (Int_t i = 0; i < nSpecies; i++) probabilities[i] = -1;
318     return;
319   }
320   for (Int_t i = 0; i < nSpecies; i++) {
321     probabilities[i] = fProbDensity[i] * prior[i] / sum;
322   }
323 }
324
325 //_____________________________________________________________________________
326 void AliPID::GetProbabilities(Double_t* probabilities) const
327 {
328 // get the probabilities to be a particle of given type
329 // assuming the globaly set a priori probabilities
330
331   GetProbabilities(probabilities, fgPrior);
332 }
333
334 //_____________________________________________________________________________
335 AliPID::EParticleType AliPID::GetMostProbable(const Double_t* prior) const
336 {
337 // get the most probable particle id hypothesis
338 // assuming the a priori probabilities "prior"
339
340   Double_t max = 0.;
341   EParticleType id = kPion;
342   Int_t nSpecies = ((fCharged) ? kSPECIES : kSPECIESN);
343   for (Int_t i = 0; i < nSpecies; i++) {
344     Double_t prob = fProbDensity[i] * prior[i];
345     if (prob > max) {
346       max = prob;
347       id = EParticleType(i);
348     }
349   }
350   if (max == 0) {
351     AliError("Invalid probability densities or priors");
352   }
353   return id;
354 }
355
356 //_____________________________________________________________________________
357 AliPID::EParticleType AliPID::GetMostProbable() const
358 {
359 // get the most probable particle id hypothesis
360 // assuming the globaly set a priori probabilities
361
362   return GetMostProbable(fgPrior);
363 }
364
365
366 //_____________________________________________________________________________
367 void AliPID::SetPriors(const Double_t* prior, Bool_t charged)
368 {
369 // use the given priors as global a priori probabilities
370
371   Double_t sum = 0;
372   for (Int_t i = 0; i < kSPECIESN; i++) {
373     if (charged && (i >= kSPECIES)) {
374       fgPrior[i] = 0;      
375     } else {
376       if (prior[i] < 0) {
377         AliWarningClass(Form("negative prior (%g) for %ss. "
378                              "Using 0 instead.", prior[i], 
379                              fgkParticleName[i]));
380         fgPrior[i] = 0;
381       } else {
382         fgPrior[i] = prior[i];
383       }
384     }
385     sum += prior[i];
386   }
387   if (sum == 0) {
388     AliWarningClass("all priors are zero.");
389   }
390 }
391
392 //_____________________________________________________________________________
393 void AliPID::SetPrior(EParticleType iType, Double_t prior)
394 {
395 // use the given prior as global a priori probability for particles
396 // of type "iType"
397
398   if (prior < 0) {
399     AliWarningClass(Form("negative prior (%g) for %ss. Using 0 instead.", 
400                          prior, fgkParticleName[iType]));
401     prior = 0;
402   }
403   fgPrior[iType] = prior;
404 }
405
406
407 //_____________________________________________________________________________
408 AliPID& AliPID::operator *= (const AliPID& pid)
409 {
410 // combine this probability densities with the one of "pid"
411
412   for (Int_t i = 0; i < kSPECIESN; i++) {
413     fProbDensity[i] *= pid.fProbDensity[i];
414   }
415   return *this;
416 }
417
418 //_____________________________________________________________________________
419 AliPID operator * (const AliPID& pid1, const AliPID& pid2)
420 {
421 // combine the two probability densities
422
423   AliPID result;
424   result *= pid1;
425   result *= pid2;
426   return result;
427 }