ecc418088e82f8d025577180712795e452b570e4
[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::kSPECIESCN+1] = {
54   "electron",
55   "muon",
56   "pion",
57   "kaon",  
58   "proton",
59   
60   "deuteron",
61   "triton",
62   "helium-3",
63   "alpha",
64   
65   "photon",
66   "pi0",
67   "neutron",
68   "kaon0",
69   "eleCon",
70   
71   "unknown"
72 };
73
74 const char* AliPID::fgkParticleShortName[AliPID::kSPECIESCN+1] = {
75   "e",
76   "mu",
77   "pi",
78   "K",
79   "p",
80
81   "d",
82   "t",
83   "he3",
84   "alpha",
85
86   "photon",
87   "pi0",
88   "n",
89   "K0",
90   "eleCon",
91   
92   "unknown"
93 };
94
95 const char* AliPID::fgkParticleLatexName[AliPID::kSPECIESCN+1] = {
96   "e",
97   "#mu",
98   "#pi",
99   "K",
100   "p",
101
102   "d",
103   "t",
104   "he3",
105   "#alpha",
106
107   "#gamma",
108   "#pi_{0}",
109   "n",
110   "K_{0}",
111   "eleCon",
112
113   "unknown"
114 };
115
116 const Int_t AliPID::fgkParticleCode[AliPID::kSPECIESCN+1] = {
117   ::kElectron, 
118   ::kMuonMinus, 
119   ::kPiPlus, 
120   ::kKPlus, 
121   ::kProton,
122   
123   1000010020,
124   1000010030,
125   1000020030,
126   1000020040,
127
128   ::kGamma,
129   ::kPi0,
130   ::kNeutron,
131   ::kK0,
132   ::kElectron,
133   0
134 };
135
136 /*const*/ Float_t AliPID::fgkParticleMass[AliPID::kSPECIESCN+1] = {
137   0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
138   /*
139   M(kElectron),  // electron
140   M(kMuon), // muon
141   M(kPion),    // pion
142   M(kKaon),     // kaon
143   M(kProton),    // proton
144   M(kPhoton),     // photon
145   M(kPi0),       // pi0
146   M(kNeutron),   // neutron
147   M(kKaon0),        // kaon0
148   M(kEleCon),     // electron conversion
149   M(kDeuteron), // deuteron
150   M(kTriton),   // triton
151   M(kHe3),      // he3
152   M(kAlpha),    // alpha
153   0.00000        // unknown
154   */
155 };
156 /*const*/ Float_t AliPID::fgkParticleMassZ[AliPID::kSPECIESCN+1] = {
157   0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
158   /*
159   M(kElectron),  // electron
160   M(kMuon), // muon
161   M(kPion),    // pion
162   M(kKaon),     // kaon
163   M(kProton),    // proton
164   M(kPhoton),     // photon
165   M(kPi0),       // pi0
166   M(kNeutron),   // neutron
167   M(kKaon0),        // kaon0
168   M(kEleCon),     // electron conversion
169   M(kDeuteron), // deuteron
170   M(kTriton),   // triton
171   M(kHe3)/2,      // he3
172   M(kAlpha)/2,    // alpha
173   0.00000        // unknown
174   */
175 };
176
177 Char_t AliPID::fgkParticleCharge[AliPID::kSPECIESCN+1] = { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 };
178
179 Double_t AliPID::fgPrior[kSPECIESCN] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
180
181
182 //_______________________________________________________________________
183 AliPID::AliPID() :
184   TObject(),
185   fCharged(0)
186 {
187   //
188   // Default constructor
189   //
190   Init();
191   // set default values (= equal probabilities)
192   for (Int_t i = 0; i < kSPECIESCN; i++)
193     fProbDensity[i] = 1./kSPECIESCN;
194 }
195
196 //_______________________________________________________________________
197 AliPID::AliPID(const Double_t* probDensity, Bool_t charged) : 
198   TObject(),
199   fCharged(charged)
200 {
201   //
202   // Standard constructor
203   //
204   Init();
205   // set given probability densities
206   for (Int_t i = 0; i < kSPECIESC; i++)
207     fProbDensity[i] = probDensity[i];
208
209   for (Int_t i = kSPECIESC; i < kSPECIESCN; i++)
210     fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
211 }
212
213 //_______________________________________________________________________
214 AliPID::AliPID(const Float_t* probDensity, Bool_t charged) :
215   TObject(),
216   fCharged(charged)
217 {
218   //
219   // Standard constructor
220   //
221   Init();
222   // set given probability densities
223   for (Int_t i = 0; i < kSPECIESC; i++) 
224     fProbDensity[i] = probDensity[i];
225
226   for (Int_t i = kSPECIESC; i < kSPECIESCN; i++) 
227     fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
228 }
229
230 //_______________________________________________________________________
231 AliPID::AliPID(const AliPID& pid) : 
232   TObject(pid),
233   fCharged(pid.fCharged)
234 {
235   //
236   // copy constructor
237   //
238   // We do not call init here, MUST already be done
239   for (Int_t i = 0; i < kSPECIESCN; i++) 
240     fProbDensity[i] = pid.fProbDensity[i];
241 }
242
243 //_______________________________________________________________________
244 void AliPID::SetProbabilities(const Double_t* probDensity, Bool_t charged) 
245 {
246   //
247   // Set the probability densities
248   //
249   for (Int_t i = 0; i < kSPECIESC; i++) 
250     fProbDensity[i] = probDensity[i];
251
252   for (Int_t i = kSPECIESC; i < kSPECIESCN; i++) 
253     fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
254 }
255
256 //_______________________________________________________________________
257 AliPID& AliPID::operator = (const AliPID& pid)
258 {
259 // assignment operator
260
261   if(this != &pid) {
262     fCharged = pid.fCharged;
263     for (Int_t i = 0; i < kSPECIESCN; i++) {
264       fProbDensity[i] = pid.fProbDensity[i];
265     }
266   }
267   return *this;
268 }
269
270 //_______________________________________________________________________
271 void AliPID::Init() 
272 {
273   //
274   // Initialise the masses, charges
275   //
276   // Initialise only once... 
277   if(!fgkParticleMass[0]) {
278     AliPDG::AddParticlesToPdgDataBase();
279     for (Int_t i = 0; i < kSPECIESC; i++) {
280       fgkParticleMass[i] = M(i);
281       if (i == kHe3 || i == kAlpha) {
282         fgkParticleMassZ[i] = M(i)/2.;
283         fgkParticleCharge[i] = 2;
284       }
285       else {
286         fgkParticleMassZ[i]=M(i);
287         fgkParticleCharge[i]=1;
288       }
289     }
290   }
291 }
292
293 //_____________________________________________________________________________
294 Double_t AliPID::GetProbability(EParticleType iType,
295                                 const Double_t* prior) const
296 {
297   //
298   // Get the probability to be a particle of type "iType"
299   // assuming the a priori probabilities "prior"
300   //
301   Double_t sum = 0.;
302   Int_t nSpecies = ((fCharged) ? kSPECIESC : kSPECIESCN);
303   for (Int_t i = 0; i < nSpecies; i++) {
304     sum += fProbDensity[i] * prior[i];
305   }
306   if (sum <= 0) {
307     AliError("Invalid probability densities or priors");
308     return -1;
309   }
310   return fProbDensity[iType] * prior[iType] / sum;
311 }
312
313 //_____________________________________________________________________________
314 Double_t AliPID::GetProbability(EParticleType iType) const
315 {
316 // get the probability to be a particle of type "iType"
317 // assuming the globaly set a priori probabilities
318
319   return GetProbability(iType, fgPrior);
320 }
321
322 //_____________________________________________________________________________
323 void AliPID::GetProbabilities(Double_t* probabilities,
324                               const Double_t* prior) const
325 {
326 // get the probabilities to be a particle of given type
327 // assuming the a priori probabilities "prior"
328
329   Double_t sum = 0.;
330   Int_t nSpecies = ((fCharged) ? kSPECIESC : kSPECIESCN);
331   for (Int_t i = 0; i < nSpecies; i++) {
332     sum += fProbDensity[i] * prior[i];
333   }
334   if (sum <= 0) {
335     AliError("Invalid probability densities or priors");
336     for (Int_t i = 0; i < nSpecies; i++) probabilities[i] = -1;
337     return;
338   }
339   for (Int_t i = 0; i < nSpecies; i++) {
340     probabilities[i] = fProbDensity[i] * prior[i] / sum;
341   }
342 }
343
344 //_____________________________________________________________________________
345 void AliPID::GetProbabilities(Double_t* probabilities) const
346 {
347 // get the probabilities to be a particle of given type
348 // assuming the globaly set a priori probabilities
349
350   GetProbabilities(probabilities, fgPrior);
351 }
352
353 //_____________________________________________________________________________
354 AliPID::EParticleType AliPID::GetMostProbable(const Double_t* prior) const
355 {
356 // get the most probable particle id hypothesis
357 // assuming the a priori probabilities "prior"
358
359   Double_t max = 0.;
360   EParticleType id = kPion;
361   Int_t nSpecies = ((fCharged) ? kSPECIESC : kSPECIESCN);
362   for (Int_t i = 0; i < nSpecies; i++) {
363     Double_t prob = fProbDensity[i] * prior[i];
364     if (prob > max) {
365       max = prob;
366       id = EParticleType(i);
367     }
368   }
369   if (max == 0) {
370     AliError("Invalid probability densities or priors");
371   }
372   return id;
373 }
374
375 //_____________________________________________________________________________
376 AliPID::EParticleType AliPID::GetMostProbable() const
377 {
378 // get the most probable particle id hypothesis
379 // assuming the globaly set a priori probabilities
380
381   return GetMostProbable(fgPrior);
382 }
383
384
385 //_____________________________________________________________________________
386 void AliPID::SetPriors(const Double_t* prior, Bool_t charged)
387 {
388 // use the given priors as global a priori probabilities
389
390   Double_t sum = 0;
391   for (Int_t i = 0; i < kSPECIESCN; i++) {
392     if (charged && (i >= kSPECIESC)) {
393       fgPrior[i] = 0;      
394     } else {
395       if (prior[i] < 0) {
396         AliWarningClass(Form("negative prior (%g) for %ss. "
397                              "Using 0 instead.", prior[i], 
398                              fgkParticleName[i]));
399         fgPrior[i] = 0;
400       } else {
401         fgPrior[i] = prior[i];
402       }
403     }
404     sum += prior[i];
405   }
406   if (sum == 0) {
407     AliWarningClass("all priors are zero.");
408   }
409 }
410
411 //_____________________________________________________________________________
412 void AliPID::SetPrior(EParticleType iType, Double_t prior)
413 {
414 // use the given prior as global a priori probability for particles
415 // of type "iType"
416
417   if (prior < 0) {
418     AliWarningClass(Form("negative prior (%g) for %ss. Using 0 instead.", 
419                          prior, fgkParticleName[iType]));
420     prior = 0;
421   }
422   fgPrior[iType] = prior;
423 }
424
425
426 //_____________________________________________________________________________
427 AliPID& AliPID::operator *= (const AliPID& pid)
428 {
429 // combine this probability densities with the one of "pid"
430
431   for (Int_t i = 0; i < kSPECIESCN; i++) {
432     fProbDensity[i] *= pid.fProbDensity[i];
433   }
434   return *this;
435 }
436
437 //_____________________________________________________________________________
438 AliPID operator * (const AliPID& pid1, const AliPID& pid2)
439 {
440 // combine the two probability densities
441
442   AliPID result;
443   result *= pid1;
444   result *= pid2;
445   return result;
446 }