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