Moving the classes that belong to the following libraries: STEERBase, ESD, CDB, AOD...
[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};
304864ab 145
146Double_t AliPID::fgPrior[kSPECIESN] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
147
148
149//_______________________________________________________________________
90e48c0c 150AliPID::AliPID() :
151 TObject(),
152 fCharged(0)
304864ab 153{
a9bbb414 154 //
155 // Default constructor
156 //
157 Init();
158 // set default values (= equal probabilities)
159 for (Int_t i = 0; i < kSPECIESN; i++)
304864ab 160 fProbDensity[i] = 1./kSPECIESN;
304864ab 161}
162
163//_______________________________________________________________________
90e48c0c 164AliPID::AliPID(const Double_t* probDensity, Bool_t charged) :
165 TObject(),
166 fCharged(charged)
304864ab 167{
90e48c0c 168 //
a9bbb414 169 // Standard constructor
90e48c0c 170 //
a9bbb414 171 Init();
172 // set given probability densities
173 for (Int_t i = 0; i < kSPECIES; i++)
304864ab 174 fProbDensity[i] = probDensity[i];
a9bbb414 175
176 for (Int_t i = kSPECIES; i < kSPECIESN; i++)
304864ab 177 fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
304864ab 178}
179
180//_______________________________________________________________________
90e48c0c 181AliPID::AliPID(const Float_t* probDensity, Bool_t charged) :
182 TObject(),
183 fCharged(charged)
304864ab 184{
90e48c0c 185 //
a9bbb414 186 // Standard constructor
90e48c0c 187 //
a9bbb414 188 Init();
189 // set given probability densities
190 for (Int_t i = 0; i < kSPECIES; i++)
304864ab 191 fProbDensity[i] = probDensity[i];
a9bbb414 192
193 for (Int_t i = kSPECIES; i < kSPECIESN; i++)
304864ab 194 fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
304864ab 195}
196
197//_______________________________________________________________________
198AliPID::AliPID(const AliPID& pid) :
199 TObject(pid),
200 fCharged(pid.fCharged)
201{
90e48c0c 202 //
203 // copy constructor
204 //
a9bbb414 205 // We do not call init here, MUST already be done
206 for (Int_t i = 0; i < kSPECIESN; i++)
304864ab 207 fProbDensity[i] = pid.fProbDensity[i];
304864ab 208}
209
210//_______________________________________________________________________
3088e2dd 211void 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//_______________________________________________________________________
304864ab 224AliPID& 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
a9bbb414 235//_______________________________________________________________________
236void AliPID::Init()
237{
238 //
239 // Initialise the masses
240 //
241 // Initialise only once...
35f27279 242 if(!fgkParticleMass[0]) {
243 AliPDG::AddParticlesToPdgDataBase();
244 for (Int_t i = 0; i < kSPECIESN + kSPECIESLN; i++)
374ef46b 245 fgkParticleMass[i] = M(i);
35f27279 246 }
a9bbb414 247}
304864ab 248
249//_____________________________________________________________________________
250Double_t AliPID::GetProbability(EParticleType iType,
251 const Double_t* prior) const
252{
a9bbb414 253 //
254 // Get the probability to be a particle of type "iType"
255 // assuming the a priori probabilities "prior"
256 //
304864ab 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//_____________________________________________________________________________
270Double_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//_____________________________________________________________________________
279void 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//_____________________________________________________________________________
301void 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//_____________________________________________________________________________
310AliPID::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//_____________________________________________________________________________
332AliPID::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//_____________________________________________________________________________
342void 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//_____________________________________________________________________________
368void 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//_____________________________________________________________________________
304864ab 383AliPID& 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//_____________________________________________________________________________
394AliPID 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}