Modifications needed to use PID framework based mass during tracking and
[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
00a38d07 53const char* AliPID::fgkParticleName[AliPID::kSPECIESCN+1] = {
304864ab 54 "electron",
55 "muon",
56 "pion",
00a38d07 57 "kaon",
304864ab 58 "proton",
00a38d07 59
60 "deuteron",
61 "triton",
62 "helium-3",
63 "alpha",
64
304864ab 65 "photon",
66 "pi0",
67 "neutron",
68 "kaon0",
69 "eleCon",
00a38d07 70
304864ab 71 "unknown"
72};
73
00a38d07 74const char* AliPID::fgkParticleShortName[AliPID::kSPECIESCN+1] = {
91dbd0dc 75 "e",
76 "mu",
77 "pi",
78 "K",
79 "p",
00a38d07 80
81 "d",
82 "t",
83 "he3",
84 "alpha",
85
91dbd0dc 86 "photon",
87 "pi0",
88 "n",
89 "K0",
90 "eleCon",
00a38d07 91
91dbd0dc 92 "unknown"
93};
94
00a38d07 95const char* AliPID::fgkParticleLatexName[AliPID::kSPECIESCN+1] = {
91dbd0dc 96 "e",
97 "#mu",
98 "#pi",
99 "K",
100 "p",
00a38d07 101
102 "d",
103 "t",
104 "he3",
105 "#alpha",
106
91dbd0dc 107 "#gamma",
108 "#pi_{0}",
109 "n",
110 "K_{0}",
111 "eleCon",
00a38d07 112
91dbd0dc 113 "unknown"
114};
115
00a38d07 116const Int_t AliPID::fgkParticleCode[AliPID::kSPECIESCN+1] = {
304864ab 117 ::kElectron,
118 ::kMuonMinus,
119 ::kPiPlus,
120 ::kKPlus,
121 ::kProton,
00a38d07 122
123 1000010020,
124 1000010030,
125 1000020030,
126 1000020040,
127
304864ab 128 ::kGamma,
129 ::kPi0,
130 ::kNeutron,
131 ::kK0,
0133aa88 132 ::kElectron,
304864ab 133 0
134};
135
00a38d07 136/*const*/ Float_t AliPID::fgkParticleMass[AliPID::kSPECIESCN+1] = {
35f27279 137 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
a9bbb414 138 /*
90e48c0c 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
35f27279 149 M(kDeuteron), // deuteron
150 M(kTriton), // triton
151 M(kHe3), // he3
152 M(kAlpha), // alpha
90e48c0c 153 0.00000 // unknown
a9bbb414 154 */
90e48c0c 155};
00a38d07 156/*const*/ Float_t AliPID::fgkParticleMassZ[AliPID::kSPECIESCN+1] = {
67376d1d 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};
304864ab 176
1d26da6d 177Char_t AliPID::fgkParticleCharge[AliPID::kSPECIESCN+1] = { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 };
178
00a38d07 179Double_t AliPID::fgPrior[kSPECIESCN] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
304864ab 180
181
182//_______________________________________________________________________
90e48c0c 183AliPID::AliPID() :
184 TObject(),
185 fCharged(0)
304864ab 186{
a9bbb414 187 //
188 // Default constructor
189 //
190 Init();
191 // set default values (= equal probabilities)
00a38d07 192 for (Int_t i = 0; i < kSPECIESCN; i++)
193 fProbDensity[i] = 1./kSPECIESCN;
304864ab 194}
195
196//_______________________________________________________________________
90e48c0c 197AliPID::AliPID(const Double_t* probDensity, Bool_t charged) :
198 TObject(),
199 fCharged(charged)
304864ab 200{
90e48c0c 201 //
a9bbb414 202 // Standard constructor
90e48c0c 203 //
a9bbb414 204 Init();
205 // set given probability densities
00a38d07 206 for (Int_t i = 0; i < kSPECIESC; i++)
304864ab 207 fProbDensity[i] = probDensity[i];
a9bbb414 208
00a38d07 209 for (Int_t i = kSPECIESC; i < kSPECIESCN; i++)
304864ab 210 fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
304864ab 211}
212
213//_______________________________________________________________________
90e48c0c 214AliPID::AliPID(const Float_t* probDensity, Bool_t charged) :
215 TObject(),
216 fCharged(charged)
304864ab 217{
90e48c0c 218 //
a9bbb414 219 // Standard constructor
90e48c0c 220 //
a9bbb414 221 Init();
222 // set given probability densities
00a38d07 223 for (Int_t i = 0; i < kSPECIESC; i++)
304864ab 224 fProbDensity[i] = probDensity[i];
a9bbb414 225
00a38d07 226 for (Int_t i = kSPECIESC; i < kSPECIESCN; i++)
304864ab 227 fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
304864ab 228}
229
230//_______________________________________________________________________
231AliPID::AliPID(const AliPID& pid) :
232 TObject(pid),
233 fCharged(pid.fCharged)
234{
90e48c0c 235 //
236 // copy constructor
237 //
a9bbb414 238 // We do not call init here, MUST already be done
00a38d07 239 for (Int_t i = 0; i < kSPECIESCN; i++)
304864ab 240 fProbDensity[i] = pid.fProbDensity[i];
304864ab 241}
242
3088e2dd 243//_______________________________________________________________________
244void AliPID::SetProbabilities(const Double_t* probDensity, Bool_t charged)
245{
246 //
247 // Set the probability densities
248 //
00a38d07 249 for (Int_t i = 0; i < kSPECIESC; i++)
3088e2dd 250 fProbDensity[i] = probDensity[i];
251
00a38d07 252 for (Int_t i = kSPECIESC; i < kSPECIESCN; i++)
3088e2dd 253 fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
254}
255
304864ab 256//_______________________________________________________________________
257AliPID& AliPID::operator = (const AliPID& pid)
258{
259// assignment operator
260
049e4804 261 if(this != &pid) {
262 fCharged = pid.fCharged;
00a38d07 263 for (Int_t i = 0; i < kSPECIESCN; i++) {
049e4804 264 fProbDensity[i] = pid.fProbDensity[i];
265 }
304864ab 266 }
267 return *this;
268}
269
a9bbb414 270//_______________________________________________________________________
271void AliPID::Init()
272{
273 //
1d26da6d 274 // Initialise the masses, charges
a9bbb414 275 //
276 // Initialise only once...
35f27279 277 if(!fgkParticleMass[0]) {
278 AliPDG::AddParticlesToPdgDataBase();
00a38d07 279 for (Int_t i = 0; i < kSPECIESC; i++) {
374ef46b 280 fgkParticleMass[i] = M(i);
1d26da6d 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 }
67376d1d 289 }
35f27279 290 }
a9bbb414 291}
304864ab 292
293//_____________________________________________________________________________
294Double_t AliPID::GetProbability(EParticleType iType,
295 const Double_t* prior) const
296{
a9bbb414 297 //
298 // Get the probability to be a particle of type "iType"
299 // assuming the a priori probabilities "prior"
300 //
304864ab 301 Double_t sum = 0.;
00a38d07 302 Int_t nSpecies = ((fCharged) ? kSPECIESC : kSPECIESCN);
304864ab 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//_____________________________________________________________________________
314Double_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//_____________________________________________________________________________
323void 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.;
00a38d07 330 Int_t nSpecies = ((fCharged) ? kSPECIESC : kSPECIESCN);
304864ab 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//_____________________________________________________________________________
345void 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//_____________________________________________________________________________
354AliPID::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;
00a38d07 361 Int_t nSpecies = ((fCharged) ? kSPECIESC : kSPECIESCN);
304864ab 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//_____________________________________________________________________________
376AliPID::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//_____________________________________________________________________________
386void 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;
00a38d07 391 for (Int_t i = 0; i < kSPECIESCN; i++) {
392 if (charged && (i >= kSPECIESC)) {
304864ab 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//_____________________________________________________________________________
412void 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
304864ab 426//_____________________________________________________________________________
427AliPID& AliPID::operator *= (const AliPID& pid)
428{
429// combine this probability densities with the one of "pid"
430
00a38d07 431 for (Int_t i = 0; i < kSPECIESCN; i++) {
304864ab 432 fProbDensity[i] *= pid.fProbDensity[i];
433 }
434 return *this;
435}
436
437//_____________________________________________________________________________
438AliPID 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}