STEER/AliPID.cxx

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