Added methods for finding a given module in the DDL map (F. Prino)
[u/mrichter/AliRoot.git] / STEER / 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"
a9bbb414 46#include "AliPID.h"
304864ab 47
90e48c0c 48#define M(PID) TDatabasePDG::Instance()->GetParticle(fgkParticleCode[(PID)])->Mass()
304864ab 49
50ClassImp(AliPID)
51
304864ab 52const char* AliPID::fgkParticleName[AliPID::kSPECIESN+1] = {
53 "electron",
54 "muon",
55 "pion",
56 "kaon",
57 "proton",
58 "photon",
59 "pi0",
60 "neutron",
61 "kaon0",
62 "eleCon",
63 "unknown"
64};
65
66const Int_t AliPID::fgkParticleCode[AliPID::kSPECIESN+1] = {
67 ::kElectron,
68 ::kMuonMinus,
69 ::kPiPlus,
70 ::kKPlus,
71 ::kProton,
72 ::kGamma,
73 ::kPi0,
74 ::kNeutron,
75 ::kK0,
76 ::kGamma,
77 0
78};
79
a9bbb414 80/*const*/ Float_t AliPID::fgkParticleMass[AliPID::kSPECIESN+1] = {
81 0,0,0,0,0,0,0,0,0,0,0
82 /*
90e48c0c 83 M(kElectron), // electron
84 M(kMuon), // muon
85 M(kPion), // pion
86 M(kKaon), // kaon
87 M(kProton), // proton
88 M(kPhoton), // photon
89 M(kPi0), // pi0
90 M(kNeutron), // neutron
91 M(kKaon0), // kaon0
92 M(kEleCon), // electron conversion
93 0.00000 // unknown
a9bbb414 94 */
90e48c0c 95};
304864ab 96
97Double_t AliPID::fgPrior[kSPECIESN] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
98
99
100//_______________________________________________________________________
90e48c0c 101AliPID::AliPID() :
102 TObject(),
103 fCharged(0)
304864ab 104{
a9bbb414 105 //
106 // Default constructor
107 //
108 Init();
109 // set default values (= equal probabilities)
110 for (Int_t i = 0; i < kSPECIESN; i++)
304864ab 111 fProbDensity[i] = 1./kSPECIESN;
304864ab 112}
113
114//_______________________________________________________________________
90e48c0c 115AliPID::AliPID(const Double_t* probDensity, Bool_t charged) :
116 TObject(),
117 fCharged(charged)
304864ab 118{
90e48c0c 119 //
a9bbb414 120 // Standard constructor
90e48c0c 121 //
a9bbb414 122 Init();
123 // set given probability densities
124 for (Int_t i = 0; i < kSPECIES; i++)
304864ab 125 fProbDensity[i] = probDensity[i];
a9bbb414 126
127 for (Int_t i = kSPECIES; i < kSPECIESN; i++)
304864ab 128 fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
304864ab 129}
130
131//_______________________________________________________________________
90e48c0c 132AliPID::AliPID(const Float_t* probDensity, Bool_t charged) :
133 TObject(),
134 fCharged(charged)
304864ab 135{
90e48c0c 136 //
a9bbb414 137 // Standard constructor
90e48c0c 138 //
a9bbb414 139 Init();
140 // set given probability densities
141 for (Int_t i = 0; i < kSPECIES; i++)
304864ab 142 fProbDensity[i] = probDensity[i];
a9bbb414 143
144 for (Int_t i = kSPECIES; i < kSPECIESN; i++)
304864ab 145 fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
304864ab 146}
147
148//_______________________________________________________________________
149AliPID::AliPID(const AliPID& pid) :
150 TObject(pid),
151 fCharged(pid.fCharged)
152{
90e48c0c 153 //
154 // copy constructor
155 //
a9bbb414 156 // We do not call init here, MUST already be done
157 for (Int_t i = 0; i < kSPECIESN; i++)
304864ab 158 fProbDensity[i] = pid.fProbDensity[i];
304864ab 159}
160
161//_______________________________________________________________________
3088e2dd 162void AliPID::SetProbabilities(const Double_t* probDensity, Bool_t charged)
163{
164 //
165 // Set the probability densities
166 //
167 for (Int_t i = 0; i < kSPECIES; i++)
168 fProbDensity[i] = probDensity[i];
169
170 for (Int_t i = kSPECIES; i < kSPECIESN; i++)
171 fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
172}
173
174//_______________________________________________________________________
304864ab 175AliPID& AliPID::operator = (const AliPID& pid)
176{
177// assignment operator
178
179 fCharged = pid.fCharged;
180 for (Int_t i = 0; i < kSPECIESN; i++) {
181 fProbDensity[i] = pid.fProbDensity[i];
182 }
183 return *this;
184}
185
a9bbb414 186//_______________________________________________________________________
187void AliPID::Init()
188{
189 //
190 // Initialise the masses
191 //
192 // Initialise only once...
193 if(!fgkParticleMass[0])
194 for (Int_t i = 0; i < kSPECIESN; i++)
374ef46b 195 fgkParticleMass[i] = M(i);
a9bbb414 196}
304864ab 197
198//_____________________________________________________________________________
199Double_t AliPID::GetProbability(EParticleType iType,
200 const Double_t* prior) const
201{
a9bbb414 202 //
203 // Get the probability to be a particle of type "iType"
204 // assuming the a priori probabilities "prior"
205 //
304864ab 206 Double_t sum = 0.;
207 Int_t nSpecies = ((fCharged) ? kSPECIES : kSPECIESN);
208 for (Int_t i = 0; i < nSpecies; i++) {
209 sum += fProbDensity[i] * prior[i];
210 }
211 if (sum <= 0) {
212 AliError("Invalid probability densities or priors");
213 return -1;
214 }
215 return fProbDensity[iType] * prior[iType] / sum;
216}
217
218//_____________________________________________________________________________
219Double_t AliPID::GetProbability(EParticleType iType) const
220{
221// get the probability to be a particle of type "iType"
222// assuming the globaly set a priori probabilities
223
224 return GetProbability(iType, fgPrior);
225}
226
227//_____________________________________________________________________________
228void AliPID::GetProbabilities(Double_t* probabilities,
229 const Double_t* prior) const
230{
231// get the probabilities to be a particle of given type
232// assuming the a priori probabilities "prior"
233
234 Double_t sum = 0.;
235 Int_t nSpecies = ((fCharged) ? kSPECIES : kSPECIESN);
236 for (Int_t i = 0; i < nSpecies; i++) {
237 sum += fProbDensity[i] * prior[i];
238 }
239 if (sum <= 0) {
240 AliError("Invalid probability densities or priors");
241 for (Int_t i = 0; i < nSpecies; i++) probabilities[i] = -1;
242 return;
243 }
244 for (Int_t i = 0; i < nSpecies; i++) {
245 probabilities[i] = fProbDensity[i] * prior[i] / sum;
246 }
247}
248
249//_____________________________________________________________________________
250void AliPID::GetProbabilities(Double_t* probabilities) const
251{
252// get the probabilities to be a particle of given type
253// assuming the globaly set a priori probabilities
254
255 GetProbabilities(probabilities, fgPrior);
256}
257
258//_____________________________________________________________________________
259AliPID::EParticleType AliPID::GetMostProbable(const Double_t* prior) const
260{
261// get the most probable particle id hypothesis
262// assuming the a priori probabilities "prior"
263
264 Double_t max = 0.;
265 EParticleType id = kPion;
266 Int_t nSpecies = ((fCharged) ? kSPECIES : kSPECIESN);
267 for (Int_t i = 0; i < nSpecies; i++) {
268 Double_t prob = fProbDensity[i] * prior[i];
269 if (prob > max) {
270 max = prob;
271 id = EParticleType(i);
272 }
273 }
274 if (max == 0) {
275 AliError("Invalid probability densities or priors");
276 }
277 return id;
278}
279
280//_____________________________________________________________________________
281AliPID::EParticleType AliPID::GetMostProbable() const
282{
283// get the most probable particle id hypothesis
284// assuming the globaly set a priori probabilities
285
286 return GetMostProbable(fgPrior);
287}
288
289
290//_____________________________________________________________________________
291void AliPID::SetPriors(const Double_t* prior, Bool_t charged)
292{
293// use the given priors as global a priori probabilities
294
295 Double_t sum = 0;
296 for (Int_t i = 0; i < kSPECIESN; i++) {
297 if (charged && (i >= kSPECIES)) {
298 fgPrior[i] = 0;
299 } else {
300 if (prior[i] < 0) {
301 AliWarningClass(Form("negative prior (%g) for %ss. "
302 "Using 0 instead.", prior[i],
303 fgkParticleName[i]));
304 fgPrior[i] = 0;
305 } else {
306 fgPrior[i] = prior[i];
307 }
308 }
309 sum += prior[i];
310 }
311 if (sum == 0) {
312 AliWarningClass("all priors are zero.");
313 }
314}
315
316//_____________________________________________________________________________
317void AliPID::SetPrior(EParticleType iType, Double_t prior)
318{
319// use the given prior as global a priori probability for particles
320// of type "iType"
321
322 if (prior < 0) {
323 AliWarningClass(Form("negative prior (%g) for %ss. Using 0 instead.",
324 prior, fgkParticleName[iType]));
325 prior = 0;
326 }
327 fgPrior[iType] = prior;
328}
329
330
331//_____________________________________________________________________________
304864ab 332AliPID& AliPID::operator *= (const AliPID& pid)
333{
334// combine this probability densities with the one of "pid"
335
336 for (Int_t i = 0; i < kSPECIESN; i++) {
337 fProbDensity[i] *= pid.fProbDensity[i];
338 }
339 return *this;
340}
341
342//_____________________________________________________________________________
343AliPID operator * (const AliPID& pid1, const AliPID& pid2)
344{
345// combine the two probability densities
346
347 AliPID result;
348 result *= pid1;
349 result *= pid2;
350 return result;
351}