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 Char_t AliPID::fgkParticleCharge[AliPID::kSPECIESCN+1] = { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 };
 178
 179 Double_t AliPID::fgPrior[kSPECIESCN] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 180
 181
 182 //_______________________________________________________________________
 183 AliPID::AliPID() :
 184   TObject(),
 185   fCharged(0)
 186 {
 187   //
 188   // Default constructor
 189   //
 190   Init();
 191   // set default values (= equal probabilities)
 192   for (Int_t i = 0; i < kSPECIESCN; i++)
 193     fProbDensity[i] = 1./kSPECIESCN;
 194 }
 195
 196 //_______________________________________________________________________
 197 AliPID::AliPID(const Double_t* probDensity, Bool_t charged) :
 198   TObject(),
 199   fCharged(charged)
 200 {
 201   //
 202   // Standard constructor
 203   //
 204   Init();
 205   // set given probability densities
 206   for (Int_t i = 0; i < kSPECIESC; i++)
 207     fProbDensity[i] = probDensity[i];
 208
 209   for (Int_t i = kSPECIESC; i < kSPECIESCN; i++)
 210     fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
 211 }
 212
 213 //_______________________________________________________________________
 214 AliPID::AliPID(const Float_t* probDensity, Bool_t charged) :
 215   TObject(),
 216   fCharged(charged)
 217 {
 218   //
 219   // Standard constructor
 220   //
 221   Init();
 222   // set given probability densities
 223   for (Int_t i = 0; i < kSPECIESC; i++)
 224     fProbDensity[i] = probDensity[i];
 225
 226   for (Int_t i = kSPECIESC; i < kSPECIESCN; i++)
 227     fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
 228 }
 229
 230 //_______________________________________________________________________
 231 AliPID::AliPID(const AliPID& pid) :
 232   TObject(pid),
 233   fCharged(pid.fCharged)
 234 {
 235   //
 236   // copy constructor
 237   //
 238   // We do not call init here, MUST already be done
 239   for (Int_t i = 0; i < kSPECIESCN; i++)
 240     fProbDensity[i] = pid.fProbDensity[i];
 241 }
 242
 243 //_______________________________________________________________________
 244 void AliPID::SetProbabilities(const Double_t* probDensity, Bool_t charged)
 245 {
 246   //
 247   // Set the probability densities
 248   //
 249   for (Int_t i = 0; i < kSPECIESC; i++)
 250     fProbDensity[i] = probDensity[i];
 251
 252   for (Int_t i = kSPECIESC; i < kSPECIESCN; i++)
 253     fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
 254 }
 255
 256 //_______________________________________________________________________
 257 AliPID& AliPID::operator = (const AliPID& pid)
 258 {
 259 // assignment operator
 260
 261   if(this != &pid) {
 262     fCharged = pid.fCharged;
 263     for (Int_t i = 0; i < kSPECIESCN; i++) {
 264       fProbDensity[i] = pid.fProbDensity[i];
 265     }
 266   }
 267   return *this;
 268 }
 269
 270 //_______________________________________________________________________
 271 void AliPID::Init()
 272 {
 273   //
 274   // Initialise the masses, charges
 275   //
 276   // Initialise only once...
 277   if(!fgkParticleMass[0]) {
 278     AliPDG::AddParticlesToPdgDataBase();
 279     for (Int_t i = 0; i < kSPECIESC; i++) {
 280       fgkParticleMass[i] = M(i);
 281       if (i == kHe3 || i == kAlpha) {
 282         fgkParticleMassZ[i] = M(i)/2.;
 283         fgkParticleCharge[i] = 2;
 284       }
 285       else {
 286         fgkParticleMassZ[i]=M(i);
 287         fgkParticleCharge[i]=1;
 288       }
 289     }
 290   }
 291 }
 292
 293 //_____________________________________________________________________________
 294 Double_t AliPID::GetProbability(EParticleType iType,
 295                                 const Double_t* prior) const
 296 {
 297   //
 298   // Get the probability to be a particle of type "iType"
 299   // assuming the a priori probabilities "prior"
 300   //
 301   Double_t sum = 0.;
 302   Int_t nSpecies = ((fCharged) ? kSPECIESC : kSPECIESCN);
 303   for (Int_t i = 0; i < nSpecies; i++) {
 304     sum += fProbDensity[i] * prior[i];
 305   }
 306   if (sum <= 0) {
 307     AliError("Invalid probability densities or priors");
 308     return -1;
 309   }
 310   return fProbDensity[iType] * prior[iType] / sum;
 311 }
 312
 313 //_____________________________________________________________________________
 314 Double_t AliPID::GetProbability(EParticleType iType) const
 315 {
 316 // get the probability to be a particle of type "iType"
 317 // assuming the globaly set a priori probabilities
 318
 319   return GetProbability(iType, fgPrior);
 320 }
 321
 322 //_____________________________________________________________________________
 323 void AliPID::GetProbabilities(Double_t* probabilities,
 324                               const Double_t* prior) const
 325 {
 326 // get the probabilities to be a particle of given type
 327 // assuming the a priori probabilities "prior"
 328
 329   Double_t sum = 0.;
 330   Int_t nSpecies = ((fCharged) ? kSPECIESC : kSPECIESCN);
 331   for (Int_t i = 0; i < nSpecies; i++) {
 332     sum += fProbDensity[i] * prior[i];
 333   }
 334   if (sum <= 0) {
 335     AliError("Invalid probability densities or priors");
 336     for (Int_t i = 0; i < nSpecies; i++) probabilities[i] = -1;
 337     return;
 338   }
 339   for (Int_t i = 0; i < nSpecies; i++) {
 340     probabilities[i] = fProbDensity[i] * prior[i] / sum;
 341   }
 342 }
 343
 344 //_____________________________________________________________________________
 345 void AliPID::GetProbabilities(Double_t* probabilities) const
 346 {
 347 // get the probabilities to be a particle of given type
 348 // assuming the globaly set a priori probabilities
 349
 350   GetProbabilities(probabilities, fgPrior);
 351 }
 352
 353 //_____________________________________________________________________________
 354 AliPID::EParticleType AliPID::GetMostProbable(const Double_t* prior) const
 355 {
 356 // get the most probable particle id hypothesis
 357 // assuming the a priori probabilities "prior"
 358
 359   Double_t max = 0.;
 360   EParticleType id = kPion;
 361   Int_t nSpecies = ((fCharged) ? kSPECIESC : kSPECIESCN);
 362   for (Int_t i = 0; i < nSpecies; i++) {
 363     Double_t prob = fProbDensity[i] * prior[i];
 364     if (prob > max) {
 365       max = prob;
 366       id = EParticleType(i);
 367     }
 368   }
 369   if (max == 0) {
 370     AliError("Invalid probability densities or priors");
 371   }
 372   return id;
 373 }
 374
 375 //_____________________________________________________________________________
 376 AliPID::EParticleType AliPID::GetMostProbable() const
 377 {
 378 // get the most probable particle id hypothesis
 379 // assuming the globaly set a priori probabilities
 380
 381   return GetMostProbable(fgPrior);
 382 }
 383
 384
 385 //_____________________________________________________________________________
 386 void AliPID::SetPriors(const Double_t* prior, Bool_t charged)
 387 {
 388 // use the given priors as global a priori probabilities
 389
 390   Double_t sum = 0;
 391   for (Int_t i = 0; i < kSPECIESCN; i++) {
 392     if (charged && (i >= kSPECIESC)) {
 393       fgPrior[i] = 0;
 394     } else {
 395       if (prior[i] < 0) {
 396         AliWarningClass(Form("negative prior (%g) for %ss. "
 397                              "Using 0 instead.", prior[i],
 398                              fgkParticleName[i]));
 399         fgPrior[i] = 0;
 400       } else {
 401         fgPrior[i] = prior[i];
 402       }
 403     }
 404     sum += prior[i];
 405   }
 406   if (sum == 0) {
 407     AliWarningClass("all priors are zero.");
 408   }
 409 }
 410
 411 //_____________________________________________________________________________
 412 void AliPID::SetPrior(EParticleType iType, Double_t prior)
 413 {
 414 // use the given prior as global a priori probability for particles
 415 // of type "iType"
 416
 417   if (prior < 0) {
 418     AliWarningClass(Form("negative prior (%g) for %ss. Using 0 instead.",
 419                          prior, fgkParticleName[iType]));
 420     prior = 0;
 421   }
 422   fgPrior[iType] = prior;
 423 }
 424
 425
 426 //_____________________________________________________________________________
 427 AliPID& AliPID::operator *= (const AliPID& pid)
 428 {
 429 // combine this probability densities with the one of "pid"
 430
 431   for (Int_t i = 0; i < kSPECIESCN; i++) {
 432     fProbDensity[i] *= pid.fProbDensity[i];
 433   }
 434   return *this;
 435 }
 436
 437 //_____________________________________________________________________________
 438 AliPID operator * (const AliPID& pid1, const AliPID& pid2)
 439 {
 440 // combine the two probability densities
 441
 442   AliPID result;
 443   result *= pid1;
 444   result *= pid2;
 445   return result;
 446 }