1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////////
20 // particle id probability densities //
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. //
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 //
36 // The implementation of this class is based on the note ... //
37 // by Iouri Belikov and Karel Safarik. //
39 ///////////////////////////////////////////////////////////////////////////////
42 #include <TDatabasePDG.h>
49 #define M(PID) TDatabasePDG::Instance()->GetParticle(fgkParticleCode[(PID)])->Mass()
53 const char* AliPID::fgkParticleName[AliPID::kSPECIESCN+1] = {
74 const char* AliPID::fgkParticleShortName[AliPID::kSPECIESCN+1] = {
95 const char* AliPID::fgkParticleLatexName[AliPID::kSPECIESCN+1] = {
116 const Int_t AliPID::fgkParticleCode[AliPID::kSPECIESCN+1] = {
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
139 M(kElectron), // electron
143 M(kProton), // proton
144 M(kPhoton), // photon
146 M(kNeutron), // neutron
148 M(kEleCon), // electron conversion
149 M(kDeuteron), // deuteron
150 M(kTriton), // triton
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
159 M(kElectron), // electron
163 M(kProton), // proton
164 M(kPhoton), // photon
166 M(kNeutron), // neutron
168 M(kEleCon), // electron conversion
169 M(kDeuteron), // deuteron
170 M(kTriton), // triton
172 M(kAlpha)/2, // alpha
177 Char_t AliPID::fgkParticleCharge[AliPID::kSPECIESCN+1] = { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 };
179 Double_t AliPID::fgPrior[kSPECIESCN] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
182 //_______________________________________________________________________
188 // Default constructor
191 // set default values (= equal probabilities)
192 for (Int_t i = 0; i < kSPECIESCN; i++)
193 fProbDensity[i] = 1./kSPECIESCN;
196 //_______________________________________________________________________
197 AliPID::AliPID(const Double_t* probDensity, Bool_t charged) :
202 // Standard constructor
205 // set given probability densities
206 for (Int_t i = 0; i < kSPECIESC; i++)
207 fProbDensity[i] = probDensity[i];
209 for (Int_t i = kSPECIESC; i < kSPECIESCN; i++)
210 fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
213 //_______________________________________________________________________
214 AliPID::AliPID(const Float_t* probDensity, Bool_t charged) :
219 // Standard constructor
222 // set given probability densities
223 for (Int_t i = 0; i < kSPECIESC; i++)
224 fProbDensity[i] = probDensity[i];
226 for (Int_t i = kSPECIESC; i < kSPECIESCN; i++)
227 fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
230 //_______________________________________________________________________
231 AliPID::AliPID(const AliPID& pid) :
233 fCharged(pid.fCharged)
238 // We do not call init here, MUST already be done
239 for (Int_t i = 0; i < kSPECIESCN; i++)
240 fProbDensity[i] = pid.fProbDensity[i];
243 //_______________________________________________________________________
244 void AliPID::SetProbabilities(const Double_t* probDensity, Bool_t charged)
247 // Set the probability densities
249 for (Int_t i = 0; i < kSPECIESC; i++)
250 fProbDensity[i] = probDensity[i];
252 for (Int_t i = kSPECIESC; i < kSPECIESCN; i++)
253 fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
256 //_______________________________________________________________________
257 AliPID& AliPID::operator = (const AliPID& pid)
259 // assignment operator
262 fCharged = pid.fCharged;
263 for (Int_t i = 0; i < kSPECIESCN; i++) {
264 fProbDensity[i] = pid.fProbDensity[i];
270 //_______________________________________________________________________
274 // Initialise the masses, charges
276 // Initialise only once...
277 if(!fgkParticleMass[0]) {
278 AliPDG::AddParticlesToPdgDataBase();
279 for (Int_t i = 0; i < kSPECIESC; i++) {
280 fgkParticleMass[i] = M(i);
281 if (i == kHe3 || i == kAlpha) {
282 fgkParticleMassZ[i] = M(i)/2.;
283 fgkParticleCharge[i] = 2;
286 fgkParticleMassZ[i]=M(i);
287 fgkParticleCharge[i]=1;
293 //_____________________________________________________________________________
294 Double_t AliPID::GetProbability(EParticleType iType,
295 const Double_t* prior) const
298 // Get the probability to be a particle of type "iType"
299 // assuming the a priori probabilities "prior"
302 Int_t nSpecies = ((fCharged) ? kSPECIESC : kSPECIESCN);
303 for (Int_t i = 0; i < nSpecies; i++) {
304 sum += fProbDensity[i] * prior[i];
307 AliError("Invalid probability densities or priors");
310 return fProbDensity[iType] * prior[iType] / sum;
313 //_____________________________________________________________________________
314 Double_t AliPID::GetProbability(EParticleType iType) const
316 // get the probability to be a particle of type "iType"
317 // assuming the globaly set a priori probabilities
319 return GetProbability(iType, fgPrior);
322 //_____________________________________________________________________________
323 void AliPID::GetProbabilities(Double_t* probabilities,
324 const Double_t* prior) const
326 // get the probabilities to be a particle of given type
327 // assuming the a priori probabilities "prior"
330 Int_t nSpecies = ((fCharged) ? kSPECIESC : kSPECIESCN);
331 for (Int_t i = 0; i < nSpecies; i++) {
332 sum += fProbDensity[i] * prior[i];
335 AliError("Invalid probability densities or priors");
336 for (Int_t i = 0; i < nSpecies; i++) probabilities[i] = -1;
339 for (Int_t i = 0; i < nSpecies; i++) {
340 probabilities[i] = fProbDensity[i] * prior[i] / sum;
344 //_____________________________________________________________________________
345 void AliPID::GetProbabilities(Double_t* probabilities) const
347 // get the probabilities to be a particle of given type
348 // assuming the globaly set a priori probabilities
350 GetProbabilities(probabilities, fgPrior);
353 //_____________________________________________________________________________
354 AliPID::EParticleType AliPID::GetMostProbable(const Double_t* prior) const
356 // get the most probable particle id hypothesis
357 // assuming the a priori probabilities "prior"
360 EParticleType id = kPion;
361 Int_t nSpecies = ((fCharged) ? kSPECIESC : kSPECIESCN);
362 for (Int_t i = 0; i < nSpecies; i++) {
363 Double_t prob = fProbDensity[i] * prior[i];
366 id = EParticleType(i);
370 AliError("Invalid probability densities or priors");
375 //_____________________________________________________________________________
376 AliPID::EParticleType AliPID::GetMostProbable() const
378 // get the most probable particle id hypothesis
379 // assuming the globaly set a priori probabilities
381 return GetMostProbable(fgPrior);
385 //_____________________________________________________________________________
386 void AliPID::SetPriors(const Double_t* prior, Bool_t charged)
388 // use the given priors as global a priori probabilities
391 for (Int_t i = 0; i < kSPECIESCN; i++) {
392 if (charged && (i >= kSPECIESC)) {
396 AliWarningClass(Form("negative prior (%g) for %ss. "
397 "Using 0 instead.", prior[i],
398 fgkParticleName[i]));
401 fgPrior[i] = prior[i];
407 AliWarningClass("all priors are zero.");
411 //_____________________________________________________________________________
412 void AliPID::SetPrior(EParticleType iType, Double_t prior)
414 // use the given prior as global a priori probability for particles
418 AliWarningClass(Form("negative prior (%g) for %ss. Using 0 instead.",
419 prior, fgkParticleName[iType]));
422 fgPrior[iType] = prior;
426 //_____________________________________________________________________________
427 AliPID& AliPID::operator *= (const AliPID& pid)
429 // combine this probability densities with the one of "pid"
431 for (Int_t i = 0; i < kSPECIESCN; i++) {
432 fProbDensity[i] *= pid.fProbDensity[i];
437 //_____________________________________________________________________________
438 AliPID operator * (const AliPID& pid1, const AliPID& pid2)
440 // combine the two probability densities