STEER/STEERBase/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 "AliPDG.h"
  47 #include "AliPID.h"
  48
  49 #define M(PID) TDatabasePDG::Instance()->GetParticle(fgkParticleCode[(PID)])->Mass()
  50
  51 ClassImp(AliPID)
  52
  53 const char* AliPID::fgkParticleName[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
  54   "electron",
  55   "muon",
  56   "pion",
  57   "kaon",
  58   "proton",
  59   "photon",
  60   "pi0",
  61   "neutron",
  62   "kaon0",
  63   "eleCon",
  64   "deuteron",
  65   "triton",
  66   "helium-3",
  67   "alpha",
  68   "unknown"
  69 };
  70
  71 const char* AliPID::fgkParticleShortName[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
  72   "e",
  73   "mu",
  74   "pi",
  75   "K",
  76   "p",
  77   "photon",
  78   "pi0",
  79   "n",
  80   "K0",
  81   "eleCon",
  82   "d",
  83   "t",
  84   "he3",
  85   "alpha",
  86   "unknown"
  87 };
  88
  89 const char* AliPID::fgkParticleLatexName[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
  90   "e",
  91   "#mu",
  92   "#pi",
  93   "K",
  94   "p",
  95   "#gamma",
  96   "#pi_{0}",
  97   "n",
  98   "K_{0}",
  99   "eleCon",
 100   "d",
 101   "t",
 102   "he3",
 103   "#alpha",
 104   "unknown"
 105 };
 106
 107 const Int_t AliPID::fgkParticleCode[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
 108   ::kElectron,
 109   ::kMuonMinus,
 110   ::kPiPlus,
 111   ::kKPlus,
 112   ::kProton,
 113   ::kGamma,
 114   ::kPi0,
 115   ::kNeutron,
 116   ::kK0,
 117   ::kElectron,
 118   1000010020,
 119   1000010030,
 120   1000020030,
 121   1000020040,
 122   0
 123 };
 124
 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
 127   /*
 128   M(kElectron),  // electron
 129   M(kMuon), // muon
 130   M(kPion),    // pion
 131   M(kKaon),     // kaon
 132   M(kProton),    // proton
 133   M(kPhoton),     // photon
 134   M(kPi0),       // pi0
 135   M(kNeutron),   // neutron
 136   M(kKaon0),        // kaon0
 137   M(kEleCon),     // electron conversion
 138   M(kDeuteron), // deuteron
 139   M(kTriton),   // triton
 140   M(kHe3),      // he3
 141   M(kAlpha),    // alpha
 142   0.00000        // unknown
 143   */
 144 };
 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
 147   /*
 148   M(kElectron),  // electron
 149   M(kMuon), // muon
 150   M(kPion),    // pion
 151   M(kKaon),     // kaon
 152   M(kProton),    // proton
 153   M(kPhoton),     // photon
 154   M(kPi0),       // pi0
 155   M(kNeutron),   // neutron
 156   M(kKaon0),        // kaon0
 157   M(kEleCon),     // electron conversion
 158   M(kDeuteron), // deuteron
 159   M(kTriton),   // triton
 160   M(kHe3)/2,      // he3
 161   M(kAlpha)/2,    // alpha
 162   0.00000        // unknown
 163   */
 164 };
 165
 166 Double_t AliPID::fgPrior[kSPECIESN] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 167
 168
 169 //_______________________________________________________________________
 170 AliPID::AliPID() :
 171   TObject(),
 172   fCharged(0)
 173 {
 174   //
 175   // Default constructor
 176   //
 177   Init();
 178   // set default values (= equal probabilities)
 179   for (Int_t i = 0; i < kSPECIESN; i++)
 180     fProbDensity[i] = 1./kSPECIESN;
 181 }
 182
 183 //_______________________________________________________________________
 184 AliPID::AliPID(const Double_t* probDensity, Bool_t charged) :
 185   TObject(),
 186   fCharged(charged)
 187 {
 188   //
 189   // Standard constructor
 190   //
 191   Init();
 192   // set given probability densities
 193   for (Int_t i = 0; i < kSPECIES; i++)
 194     fProbDensity[i] = probDensity[i];
 195
 196   for (Int_t i = kSPECIES; i < kSPECIESN; i++)
 197     fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
 198 }
 199
 200 //_______________________________________________________________________
 201 AliPID::AliPID(const Float_t* probDensity, Bool_t charged) :
 202   TObject(),
 203   fCharged(charged)
 204 {
 205   //
 206   // Standard constructor
 207   //
 208   Init();
 209   // set given probability densities
 210   for (Int_t i = 0; i < kSPECIES; i++)
 211     fProbDensity[i] = probDensity[i];
 212
 213   for (Int_t i = kSPECIES; i < kSPECIESN; i++)
 214     fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
 215 }
 216
 217 //_______________________________________________________________________
 218 AliPID::AliPID(const AliPID& pid) :
 219   TObject(pid),
 220   fCharged(pid.fCharged)
 221 {
 222   //
 223   // copy constructor
 224   //
 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];
 228 }
 229
 230 //_______________________________________________________________________
 231 void AliPID::SetProbabilities(const Double_t* probDensity, Bool_t charged)
 232 {
 233   //
 234   // Set the probability densities
 235   //
 236   for (Int_t i = 0; i < kSPECIES; i++)
 237     fProbDensity[i] = probDensity[i];
 238
 239   for (Int_t i = kSPECIES; i < kSPECIESN; i++)
 240     fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
 241 }
 242
 243 //_______________________________________________________________________
 244 AliPID& AliPID::operator = (const AliPID& pid)
 245 {
 246 // assignment operator
 247
 248   if(this != &pid) {
 249     fCharged = pid.fCharged;
 250     for (Int_t i = 0; i < kSPECIESN; i++) {
 251       fProbDensity[i] = pid.fProbDensity[i];
 252     }
 253   }
 254   return *this;
 255 }
 256
 257 //_______________________________________________________________________
 258 void AliPID::Init()
 259 {
 260   //
 261   // Initialise the masses
 262   //
 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);
 270     }
 271   }
 272 }
 273
 274 //_____________________________________________________________________________
 275 Double_t AliPID::GetProbability(EParticleType iType,
 276                                 const Double_t* prior) const
 277 {
 278   //
 279   // Get the probability to be a particle of type "iType"
 280   // assuming the a priori probabilities "prior"
 281   //
 282   Double_t sum = 0.;
 283   Int_t nSpecies = ((fCharged) ? kSPECIES : kSPECIESN);
 284   for (Int_t i = 0; i < nSpecies; i++) {
 285     sum += fProbDensity[i] * prior[i];
 286   }
 287   if (sum <= 0) {
 288     AliError("Invalid probability densities or priors");
 289     return -1;
 290   }
 291   return fProbDensity[iType] * prior[iType] / sum;
 292 }
 293
 294 //_____________________________________________________________________________
 295 Double_t AliPID::GetProbability(EParticleType iType) const
 296 {
 297 // get the probability to be a particle of type "iType"
 298 // assuming the globaly set a priori probabilities
 299
 300   return GetProbability(iType, fgPrior);
 301 }
 302
 303 //_____________________________________________________________________________
 304 void AliPID::GetProbabilities(Double_t* probabilities,
 305                               const Double_t* prior) const
 306 {
 307 // get the probabilities to be a particle of given type
 308 // assuming the a priori probabilities "prior"
 309
 310   Double_t sum = 0.;
 311   Int_t nSpecies = ((fCharged) ? kSPECIES : kSPECIESN);
 312   for (Int_t i = 0; i < nSpecies; i++) {
 313     sum += fProbDensity[i] * prior[i];
 314   }
 315   if (sum <= 0) {
 316     AliError("Invalid probability densities or priors");
 317     for (Int_t i = 0; i < nSpecies; i++) probabilities[i] = -1;
 318     return;
 319   }
 320   for (Int_t i = 0; i < nSpecies; i++) {
 321     probabilities[i] = fProbDensity[i] * prior[i] / sum;
 322   }
 323 }
 324
 325 //_____________________________________________________________________________
 326 void AliPID::GetProbabilities(Double_t* probabilities) const
 327 {
 328 // get the probabilities to be a particle of given type
 329 // assuming the globaly set a priori probabilities
 330
 331   GetProbabilities(probabilities, fgPrior);
 332 }
 333
 334 //_____________________________________________________________________________
 335 AliPID::EParticleType AliPID::GetMostProbable(const Double_t* prior) const
 336 {
 337 // get the most probable particle id hypothesis
 338 // assuming the a priori probabilities "prior"
 339
 340   Double_t max = 0.;
 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];
 345     if (prob > max) {
 346       max = prob;
 347       id = EParticleType(i);
 348     }
 349   }
 350   if (max == 0) {
 351     AliError("Invalid probability densities or priors");
 352   }
 353   return id;
 354 }
 355
 356 //_____________________________________________________________________________
 357 AliPID::EParticleType AliPID::GetMostProbable() const
 358 {
 359 // get the most probable particle id hypothesis
 360 // assuming the globaly set a priori probabilities
 361
 362   return GetMostProbable(fgPrior);
 363 }
 364
 365
 366 //_____________________________________________________________________________
 367 void AliPID::SetPriors(const Double_t* prior, Bool_t charged)
 368 {
 369 // use the given priors as global a priori probabilities
 370
 371   Double_t sum = 0;
 372   for (Int_t i = 0; i < kSPECIESN; i++) {
 373     if (charged && (i >= kSPECIES)) {
 374       fgPrior[i] = 0;
 375     } else {
 376       if (prior[i] < 0) {
 377         AliWarningClass(Form("negative prior (%g) for %ss. "
 378                              "Using 0 instead.", prior[i],
 379                              fgkParticleName[i]));
 380         fgPrior[i] = 0;
 381       } else {
 382         fgPrior[i] = prior[i];
 383       }
 384     }
 385     sum += prior[i];
 386   }
 387   if (sum == 0) {
 388     AliWarningClass("all priors are zero.");
 389   }
 390 }
 391
 392 //_____________________________________________________________________________
 393 void AliPID::SetPrior(EParticleType iType, Double_t prior)
 394 {
 395 // use the given prior as global a priori probability for particles
 396 // of type "iType"
 397
 398   if (prior < 0) {
 399     AliWarningClass(Form("negative prior (%g) for %ss. Using 0 instead.",
 400                          prior, fgkParticleName[iType]));
 401     prior = 0;
 402   }
 403   fgPrior[iType] = prior;
 404 }
 405
 406
 407 //_____________________________________________________________________________
 408 AliPID& AliPID::operator *= (const AliPID& pid)
 409 {
 410 // combine this probability densities with the one of "pid"
 411
 412   for (Int_t i = 0; i < kSPECIESN; i++) {
 413     fProbDensity[i] *= pid.fProbDensity[i];
 414   }
 415   return *this;
 416 }
 417
 418 //_____________________________________________________________________________
 419 AliPID operator * (const AliPID& pid1, const AliPID& pid2)
 420 {
 421 // combine the two probability densities
 422
 423   AliPID result;
 424   result *= pid1;
 425   result *= pid2;
 426   return result;
 427 }