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