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