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