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::kSPECIESN+AliPID::kSPECIESLN+1] = {
71 const char* AliPID::fgkParticleShortName[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
89 const char* AliPID::fgkParticleLatexName[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
107 const Int_t AliPID::fgkParticleCode[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
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
128 M(kElectron), // electron
132 M(kProton), // proton
133 M(kPhoton), // photon
135 M(kNeutron), // neutron
137 M(kEleCon), // electron conversion
138 M(kDeuteron), // deuteron
139 M(kTriton), // triton
145 /*const*/ Float_t AliPID::fgkParticleMassZ[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
146 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
148 M(kElectron), // electron
152 M(kProton), // proton
153 M(kPhoton), // photon
155 M(kNeutron), // neutron
157 M(kEleCon), // electron conversion
158 M(kDeuteron), // deuteron
159 M(kTriton), // triton
161 M(kAlpha)/2, // alpha
166 Double_t AliPID::fgPrior[kSPECIESN] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
169 //_______________________________________________________________________
175 // Default constructor
178 // set default values (= equal probabilities)
179 for (Int_t i = 0; i < kSPECIESN; i++)
180 fProbDensity[i] = 1./kSPECIESN;
183 //_______________________________________________________________________
184 AliPID::AliPID(const Double_t* probDensity, Bool_t charged) :
189 // Standard constructor
192 // set given probability densities
193 for (Int_t i = 0; i < kSPECIES; i++)
194 fProbDensity[i] = probDensity[i];
196 for (Int_t i = kSPECIES; i < kSPECIESN; i++)
197 fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
200 //_______________________________________________________________________
201 AliPID::AliPID(const Float_t* probDensity, Bool_t charged) :
206 // Standard constructor
209 // set given probability densities
210 for (Int_t i = 0; i < kSPECIES; i++)
211 fProbDensity[i] = probDensity[i];
213 for (Int_t i = kSPECIES; i < kSPECIESN; i++)
214 fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
217 //_______________________________________________________________________
218 AliPID::AliPID(const AliPID& pid) :
220 fCharged(pid.fCharged)
225 // We do not call init here, MUST already be done
226 for (Int_t i = 0; i < kSPECIESN; i++)
227 fProbDensity[i] = pid.fProbDensity[i];
230 //_______________________________________________________________________
231 void AliPID::SetProbabilities(const Double_t* probDensity, Bool_t charged)
234 // Set the probability densities
236 for (Int_t i = 0; i < kSPECIES; i++)
237 fProbDensity[i] = probDensity[i];
239 for (Int_t i = kSPECIES; i < kSPECIESN; i++)
240 fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
243 //_______________________________________________________________________
244 AliPID& AliPID::operator = (const AliPID& pid)
246 // assignment operator
249 fCharged = pid.fCharged;
250 for (Int_t i = 0; i < kSPECIESN; i++) {
251 fProbDensity[i] = pid.fProbDensity[i];
257 //_______________________________________________________________________
261 // Initialise the masses
263 // Initialise only once...
264 if(!fgkParticleMass[0]) {
265 AliPDG::AddParticlesToPdgDataBase();
266 for (Int_t i = 0; i < kSPECIESN + kSPECIESLN; i++) {
267 fgkParticleMass[i] = M(i);
268 if (i == kHe3 || i == kAlpha) fgkParticleMassZ[i] = M(i)/2.;
269 else fgkParticleMassZ[i]=M(i);
274 //_____________________________________________________________________________
275 Double_t AliPID::GetProbability(EParticleType iType,
276 const Double_t* prior) const
279 // Get the probability to be a particle of type "iType"
280 // assuming the a priori probabilities "prior"
283 Int_t nSpecies = ((fCharged) ? kSPECIES : kSPECIESN);
284 for (Int_t i = 0; i < nSpecies; i++) {
285 sum += fProbDensity[i] * prior[i];
288 AliError("Invalid probability densities or priors");
291 return fProbDensity[iType] * prior[iType] / sum;
294 //_____________________________________________________________________________
295 Double_t AliPID::GetProbability(EParticleType iType) const
297 // get the probability to be a particle of type "iType"
298 // assuming the globaly set a priori probabilities
300 return GetProbability(iType, fgPrior);
303 //_____________________________________________________________________________
304 void AliPID::GetProbabilities(Double_t* probabilities,
305 const Double_t* prior) const
307 // get the probabilities to be a particle of given type
308 // assuming the a priori probabilities "prior"
311 Int_t nSpecies = ((fCharged) ? kSPECIES : kSPECIESN);
312 for (Int_t i = 0; i < nSpecies; i++) {
313 sum += fProbDensity[i] * prior[i];
316 AliError("Invalid probability densities or priors");
317 for (Int_t i = 0; i < nSpecies; i++) probabilities[i] = -1;
320 for (Int_t i = 0; i < nSpecies; i++) {
321 probabilities[i] = fProbDensity[i] * prior[i] / sum;
325 //_____________________________________________________________________________
326 void AliPID::GetProbabilities(Double_t* probabilities) const
328 // get the probabilities to be a particle of given type
329 // assuming the globaly set a priori probabilities
331 GetProbabilities(probabilities, fgPrior);
334 //_____________________________________________________________________________
335 AliPID::EParticleType AliPID::GetMostProbable(const Double_t* prior) const
337 // get the most probable particle id hypothesis
338 // assuming the a priori probabilities "prior"
341 EParticleType id = kPion;
342 Int_t nSpecies = ((fCharged) ? kSPECIES : kSPECIESN);
343 for (Int_t i = 0; i < nSpecies; i++) {
344 Double_t prob = fProbDensity[i] * prior[i];
347 id = EParticleType(i);
351 AliError("Invalid probability densities or priors");
356 //_____________________________________________________________________________
357 AliPID::EParticleType AliPID::GetMostProbable() const
359 // get the most probable particle id hypothesis
360 // assuming the globaly set a priori probabilities
362 return GetMostProbable(fgPrior);
366 //_____________________________________________________________________________
367 void AliPID::SetPriors(const Double_t* prior, Bool_t charged)
369 // use the given priors as global a priori probabilities
372 for (Int_t i = 0; i < kSPECIESN; i++) {
373 if (charged && (i >= kSPECIES)) {
377 AliWarningClass(Form("negative prior (%g) for %ss. "
378 "Using 0 instead.", prior[i],
379 fgkParticleName[i]));
382 fgPrior[i] = prior[i];
388 AliWarningClass("all priors are zero.");
392 //_____________________________________________________________________________
393 void AliPID::SetPrior(EParticleType iType, Double_t prior)
395 // use the given prior as global a priori probability for particles
399 AliWarningClass(Form("negative prior (%g) for %ss. Using 0 instead.",
400 prior, fgkParticleName[iType]));
403 fgPrior[iType] = prior;
407 //_____________________________________________________________________________
408 AliPID& AliPID::operator *= (const AliPID& pid)
410 // combine this probability densities with the one of "pid"
412 for (Int_t i = 0; i < kSPECIESN; i++) {
413 fProbDensity[i] *= pid.fProbDensity[i];
418 //_____________________________________________________________________________
419 AliPID operator * (const AliPID& pid1, const AliPID& pid2)
421 // combine the two probability densities