Line data Source code
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 176 : 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 10 : TObject(),
185 10 : fCharged(0)
186 50 : {
187 : //
188 : // Default constructor
189 : //
190 10 : Init();
191 : // set default values (= equal probabilities)
192 300 : for (Int_t i = 0; i < kSPECIESCN; i++)
193 140 : fProbDensity[i] = 1./kSPECIESCN;
194 20 : }
195 :
196 : //_______________________________________________________________________
197 : AliPID::AliPID(const Double_t* probDensity, Bool_t charged) :
198 0 : TObject(),
199 0 : fCharged(charged)
200 0 : {
201 : //
202 : // Standard constructor
203 : //
204 0 : Init();
205 : // set given probability densities
206 0 : for (Int_t i = 0; i < kSPECIESC; i++)
207 0 : fProbDensity[i] = probDensity[i];
208 :
209 0 : for (Int_t i = kSPECIESC; i < kSPECIESCN; i++)
210 0 : fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
211 0 : }
212 :
213 : //_______________________________________________________________________
214 : AliPID::AliPID(const Float_t* probDensity, Bool_t charged) :
215 0 : TObject(),
216 0 : fCharged(charged)
217 0 : {
218 : //
219 : // Standard constructor
220 : //
221 0 : Init();
222 : // set given probability densities
223 0 : for (Int_t i = 0; i < kSPECIESC; i++)
224 0 : fProbDensity[i] = probDensity[i];
225 :
226 0 : for (Int_t i = kSPECIESC; i < kSPECIESCN; i++)
227 0 : fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
228 0 : }
229 :
230 : //_______________________________________________________________________
231 : AliPID::AliPID(const AliPID& pid) :
232 0 : TObject(pid),
233 0 : fCharged(pid.fCharged)
234 0 : {
235 : //
236 : // copy constructor
237 : //
238 : // We do not call init here, MUST already be done
239 0 : for (Int_t i = 0; i < kSPECIESCN; i++)
240 0 : fProbDensity[i] = pid.fProbDensity[i];
241 0 : }
242 :
243 : //_______________________________________________________________________
244 : void AliPID::SetProbabilities(const Double_t* probDensity, Bool_t charged)
245 : {
246 : //
247 : // Set the probability densities
248 : //
249 0 : for (Int_t i = 0; i < kSPECIESC; i++)
250 0 : fProbDensity[i] = probDensity[i];
251 :
252 0 : for (Int_t i = kSPECIESC; i < kSPECIESCN; i++)
253 0 : fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
254 0 : }
255 :
256 : //_______________________________________________________________________
257 : AliPID& AliPID::operator = (const AliPID& pid)
258 : {
259 : // assignment operator
260 :
261 0 : if(this != &pid) {
262 0 : fCharged = pid.fCharged;
263 0 : for (Int_t i = 0; i < kSPECIESCN; i++) {
264 0 : fProbDensity[i] = pid.fProbDensity[i];
265 : }
266 0 : }
267 0 : return *this;
268 : }
269 :
270 : //_______________________________________________________________________
271 : void AliPID::Init()
272 : {
273 : //
274 : // Initialise the masses, charges
275 : //
276 : // Initialise only once...
277 24 : if(!fgkParticleMass[0]) {
278 4 : AliPDG::AddParticlesToPdgDataBase();
279 80 : for (Int_t i = 0; i < kSPECIESC; i++) {
280 36 : fgkParticleMass[i] = M(i);
281 72 : if (i == kHe3 || i == kAlpha) {
282 44 : fgkParticleMassZ[i] = M(i)/2.;
283 8 : fgkParticleCharge[i] = 2;
284 8 : }
285 : else {
286 28 : fgkParticleMassZ[i]=M(i);
287 28 : fgkParticleCharge[i]=1;
288 : }
289 : }
290 4 : }
291 12 : }
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 0 : Int_t nSpecies = ((fCharged) ? kSPECIESC : kSPECIESCN);
303 0 : for (Int_t i = 0; i < nSpecies; i++) {
304 0 : sum += fProbDensity[i] * prior[i];
305 : }
306 0 : if (sum <= 0) {
307 0 : AliError("Invalid probability densities or priors");
308 0 : return -1;
309 : }
310 0 : return fProbDensity[iType] * prior[iType] / sum;
311 0 : }
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 0 : 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 0 : Int_t nSpecies = ((fCharged) ? kSPECIESC : kSPECIESCN);
331 0 : for (Int_t i = 0; i < nSpecies; i++) {
332 0 : sum += fProbDensity[i] * prior[i];
333 : }
334 0 : if (sum <= 0) {
335 0 : AliError("Invalid probability densities or priors");
336 0 : for (Int_t i = 0; i < nSpecies; i++) probabilities[i] = -1;
337 0 : return;
338 : }
339 0 : for (Int_t i = 0; i < nSpecies; i++) {
340 0 : probabilities[i] = fProbDensity[i] * prior[i] / sum;
341 : }
342 0 : }
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 0 : GetProbabilities(probabilities, fgPrior);
351 0 : }
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 0 : Int_t nSpecies = ((fCharged) ? kSPECIESC : kSPECIESCN);
362 0 : for (Int_t i = 0; i < nSpecies; i++) {
363 0 : Double_t prob = fProbDensity[i] * prior[i];
364 0 : if (prob > max) {
365 : max = prob;
366 : id = EParticleType(i);
367 0 : }
368 : }
369 0 : if (max == 0) {
370 0 : AliError("Invalid probability densities or priors");
371 0 : }
372 0 : 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 0 : 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 0 : for (Int_t i = 0; i < kSPECIESCN; i++) {
392 0 : if (charged && (i >= kSPECIESC)) {
393 0 : fgPrior[i] = 0;
394 0 : } else {
395 0 : if (prior[i] < 0) {
396 0 : AliWarningClass(Form("negative prior (%g) for %ss. "
397 : "Using 0 instead.", prior[i],
398 : fgkParticleName[i]));
399 0 : fgPrior[i] = 0;
400 0 : } else {
401 0 : fgPrior[i] = prior[i];
402 : }
403 : }
404 0 : sum += prior[i];
405 : }
406 0 : if (sum == 0) {
407 0 : AliWarningClass("all priors are zero.");
408 0 : }
409 0 : }
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 0 : if (prior < 0) {
418 0 : AliWarningClass(Form("negative prior (%g) for %ss. Using 0 instead.",
419 : prior, fgkParticleName[iType]));
420 : prior = 0;
421 0 : }
422 0 : fgPrior[iType] = prior;
423 0 : }
424 :
425 :
426 : //_____________________________________________________________________________
427 : AliPID& AliPID::operator *= (const AliPID& pid)
428 : {
429 : // combine this probability densities with the one of "pid"
430 :
431 0 : for (Int_t i = 0; i < kSPECIESCN; i++) {
432 0 : fProbDensity[i] *= pid.fProbDensity[i];
433 : }
434 0 : return *this;
435 : }
436 :
437 : //_____________________________________________________________________________
438 : AliPID operator * (const AliPID& pid1, const AliPID& pid2)
439 : {
440 : // combine the two probability densities
441 :
442 0 : AliPID result;
443 0 : result *= pid1;
444 0 : result *= pid2;
445 : return result;
446 0 : }
|
