Changes for #95543: request to commit to trunk: TOF code supporting PID for light...
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliPID.cxx
CommitLineData
304864ab 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
a9bbb414 41#include <TClass.h>
42#include <TDatabasePDG.h>
43#include <TPDGCode.h>
304864ab 44
304864ab 45#include "AliLog.h"
35f27279 46#include "AliPDG.h"
a9bbb414 47#include "AliPID.h"
304864ab 48
90e48c0c 49#define M(PID) TDatabasePDG::Instance()->GetParticle(fgkParticleCode[(PID)])->Mass()
304864ab 50
51ClassImp(AliPID)
52
35f27279 53const char* AliPID::fgkParticleName[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
304864ab 54 "electron",
55 "muon",
56 "pion",
57 "kaon",
58 "proton",
59 "photon",
60 "pi0",
61 "neutron",
62 "kaon0",
63 "eleCon",
35f27279 64 "deuteron",
65 "triton",
66 "helium-3",
67 "alpha",
304864ab 68 "unknown"
69};
70
35f27279 71const char* AliPID::fgkParticleShortName[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
91dbd0dc 72 "e",
73 "mu",
74 "pi",
75 "K",
76 "p",
77 "photon",
78 "pi0",
79 "n",
80 "K0",
81 "eleCon",
35f27279 82 "d",
83 "t",
84 "he3",
85 "alpha",
91dbd0dc 86 "unknown"
87};
88
35f27279 89const char* AliPID::fgkParticleLatexName[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
91dbd0dc 90 "e",
91 "#mu",
92 "#pi",
93 "K",
94 "p",
95 "#gamma",
96 "#pi_{0}",
97 "n",
98 "K_{0}",
99 "eleCon",
35f27279 100 "d",
101 "t",
102 "he3",
103 "#alpha",
91dbd0dc 104 "unknown"
105};
106
35f27279 107const Int_t AliPID::fgkParticleCode[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
304864ab 108 ::kElectron,
109 ::kMuonMinus,
110 ::kPiPlus,
111 ::kKPlus,
112 ::kProton,
113 ::kGamma,
114 ::kPi0,
115 ::kNeutron,
116 ::kK0,
0133aa88 117 ::kElectron,
35f27279 118 1000010020,
119 1000010030,
120 1000020030,
121 1000020040,
304864ab 122 0
123};
124
35f27279 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
a9bbb414 127 /*
90e48c0c 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
35f27279 138 M(kDeuteron), // deuteron
139 M(kTriton), // triton
140 M(kHe3), // he3
141 M(kAlpha), // alpha
90e48c0c 142 0.00000 // unknown
a9bbb414 143 */
90e48c0c 144};
67376d1d 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};
304864ab 165
166Double_t AliPID::fgPrior[kSPECIESN] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
167
168
169//_______________________________________________________________________
90e48c0c 170AliPID::AliPID() :
171 TObject(),
172 fCharged(0)
304864ab 173{
a9bbb414 174 //
175 // Default constructor
176 //
177 Init();
178 // set default values (= equal probabilities)
179 for (Int_t i = 0; i < kSPECIESN; i++)
304864ab 180 fProbDensity[i] = 1./kSPECIESN;
304864ab 181}
182
183//_______________________________________________________________________
90e48c0c 184AliPID::AliPID(const Double_t* probDensity, Bool_t charged) :
185 TObject(),
186 fCharged(charged)
304864ab 187{
90e48c0c 188 //
a9bbb414 189 // Standard constructor
90e48c0c 190 //
a9bbb414 191 Init();
192 // set given probability densities
193 for (Int_t i = 0; i < kSPECIES; i++)
304864ab 194 fProbDensity[i] = probDensity[i];
a9bbb414 195
196 for (Int_t i = kSPECIES; i < kSPECIESN; i++)
304864ab 197 fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
304864ab 198}
199
200//_______________________________________________________________________
90e48c0c 201AliPID::AliPID(const Float_t* probDensity, Bool_t charged) :
202 TObject(),
203 fCharged(charged)
304864ab 204{
90e48c0c 205 //
a9bbb414 206 // Standard constructor
90e48c0c 207 //
a9bbb414 208 Init();
209 // set given probability densities
210 for (Int_t i = 0; i < kSPECIES; i++)
304864ab 211 fProbDensity[i] = probDensity[i];
a9bbb414 212
213 for (Int_t i = kSPECIES; i < kSPECIESN; i++)
304864ab 214 fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
304864ab 215}
216
217//_______________________________________________________________________
218AliPID::AliPID(const AliPID& pid) :
219 TObject(pid),
220 fCharged(pid.fCharged)
221{
90e48c0c 222 //
223 // copy constructor
224 //
a9bbb414 225 // We do not call init here, MUST already be done
226 for (Int_t i = 0; i < kSPECIESN; i++)
304864ab 227 fProbDensity[i] = pid.fProbDensity[i];
304864ab 228}
229
3088e2dd 230//_______________________________________________________________________
231void 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
304864ab 243//_______________________________________________________________________
244AliPID& AliPID::operator = (const AliPID& pid)
245{
246// assignment operator
247
049e4804 248 if(this != &pid) {
249 fCharged = pid.fCharged;
250 for (Int_t i = 0; i < kSPECIESN; i++) {
251 fProbDensity[i] = pid.fProbDensity[i];
252 }
304864ab 253 }
254 return *this;
255}
256
a9bbb414 257//_______________________________________________________________________
258void AliPID::Init()
259{
260 //
261 // Initialise the masses
262 //
263 // Initialise only once...
35f27279 264 if(!fgkParticleMass[0]) {
265 AliPDG::AddParticlesToPdgDataBase();
67376d1d 266 for (Int_t i = 0; i < kSPECIESN + kSPECIESLN; i++) {
374ef46b 267 fgkParticleMass[i] = M(i);
67376d1d 268 if (i == kHe3 || i == kAlpha) fgkParticleMassZ[i] = M(i)/2.;
269 else fgkParticleMassZ[i]=M(i);
270 }
35f27279 271 }
a9bbb414 272}
304864ab 273
274//_____________________________________________________________________________
275Double_t AliPID::GetProbability(EParticleType iType,
276 const Double_t* prior) const
277{
a9bbb414 278 //
279 // Get the probability to be a particle of type "iType"
280 // assuming the a priori probabilities "prior"
281 //
304864ab 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//_____________________________________________________________________________
295Double_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//_____________________________________________________________________________
304void 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//_____________________________________________________________________________
326void 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//_____________________________________________________________________________
335AliPID::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//_____________________________________________________________________________
357AliPID::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//_____________________________________________________________________________
367void 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//_____________________________________________________________________________
393void 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
304864ab 407//_____________________________________________________________________________
408AliPID& 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//_____________________________________________________________________________
419AliPID 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}