#include <TClass.h>
#include <TDatabasePDG.h>
#include <TPDGCode.h>
#include "AliLog.h"
#include "AliPDG.h"
#include "AliPID.h"
#define M(PID) TDatabasePDG::Instance()->GetParticle(fgkParticleCode[(PID)])->Mass()
ClassImp(AliPID)
const char* AliPID::fgkParticleName[AliPID::kSPECIESCN+1] = {
"electron",
"muon",
"pion",
"kaon",
"proton",
"deuteron",
"triton",
"helium-3",
"alpha",
"photon",
"pi0",
"neutron",
"kaon0",
"eleCon",
"unknown"
};
const char* AliPID::fgkParticleShortName[AliPID::kSPECIESCN+1] = {
"e",
"mu",
"pi",
"K",
"p",
"d",
"t",
"he3",
"alpha",
"photon",
"pi0",
"n",
"K0",
"eleCon",
"unknown"
};
const char* AliPID::fgkParticleLatexName[AliPID::kSPECIESCN+1] = {
"e",
"#mu",
"#pi",
"K",
"p",
"d",
"t",
"he3",
"#alpha",
"#gamma",
"#pi_{0}",
"n",
"K_{0}",
"eleCon",
"unknown"
};
const Int_t AliPID::fgkParticleCode[AliPID::kSPECIESCN+1] = {
::kElectron,
::kMuonMinus,
::kPiPlus,
::kKPlus,
::kProton,
1000010020,
1000010030,
1000020030,
1000020040,
::kGamma,
::kPi0,
::kNeutron,
::kK0,
::kElectron,
0
};
Float_t AliPID::fgkParticleMass[AliPID::kSPECIESCN+1] = {
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
};
Float_t AliPID::fgkParticleMassZ[AliPID::kSPECIESCN+1] = {
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
};
Char_t AliPID::fgkParticleCharge[AliPID::kSPECIESCN+1] = { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 };
Double_t AliPID::fgPrior[kSPECIESCN] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
AliPID::AliPID() :
TObject(),
fCharged(0)
{
Init();
for (Int_t i = 0; i < kSPECIESCN; i++)
fProbDensity[i] = 1./kSPECIESCN;
}
AliPID::AliPID(const Double_t* probDensity, Bool_t charged) :
TObject(),
fCharged(charged)
{
Init();
for (Int_t i = 0; i < kSPECIESC; i++)
fProbDensity[i] = probDensity[i];
for (Int_t i = kSPECIESC; i < kSPECIESCN; i++)
fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
}
AliPID::AliPID(const Float_t* probDensity, Bool_t charged) :
TObject(),
fCharged(charged)
{
Init();
for (Int_t i = 0; i < kSPECIESC; i++)
fProbDensity[i] = probDensity[i];
for (Int_t i = kSPECIESC; i < kSPECIESCN; i++)
fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
}
AliPID::AliPID(const AliPID& pid) :
TObject(pid),
fCharged(pid.fCharged)
{
for (Int_t i = 0; i < kSPECIESCN; i++)
fProbDensity[i] = pid.fProbDensity[i];
}
void AliPID::SetProbabilities(const Double_t* probDensity, Bool_t charged)
{
for (Int_t i = 0; i < kSPECIESC; i++)
fProbDensity[i] = probDensity[i];
for (Int_t i = kSPECIESC; i < kSPECIESCN; i++)
fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
}
AliPID& AliPID::operator = (const AliPID& pid)
{
if(this != &pid) {
fCharged = pid.fCharged;
for (Int_t i = 0; i < kSPECIESCN; i++) {
fProbDensity[i] = pid.fProbDensity[i];
}
}
return *this;
}
void AliPID::Init()
{
if(!fgkParticleMass[0]) {
AliPDG::AddParticlesToPdgDataBase();
for (Int_t i = 0; i < kSPECIESC; i++) {
fgkParticleMass[i] = M(i);
if (i == kHe3 || i == kAlpha) {
fgkParticleMassZ[i] = M(i)/2.;
fgkParticleCharge[i] = 2;
}
else {
fgkParticleMassZ[i]=M(i);
fgkParticleCharge[i]=1;
}
}
}
}
Double_t AliPID::GetProbability(EParticleType iType,
const Double_t* prior) const
{
Double_t sum = 0.;
Int_t nSpecies = ((fCharged) ? kSPECIESC : kSPECIESCN);
for (Int_t i = 0; i < nSpecies; i++) {
sum += fProbDensity[i] * prior[i];
}
if (sum <= 0) {
AliError("Invalid probability densities or priors");
return -1;
}
return fProbDensity[iType] * prior[iType] / sum;
}
Double_t AliPID::GetProbability(EParticleType iType) const
{
return GetProbability(iType, fgPrior);
}
void AliPID::GetProbabilities(Double_t* probabilities,
const Double_t* prior) const
{
Double_t sum = 0.;
Int_t nSpecies = ((fCharged) ? kSPECIESC : kSPECIESCN);
for (Int_t i = 0; i < nSpecies; i++) {
sum += fProbDensity[i] * prior[i];
}
if (sum <= 0) {
AliError("Invalid probability densities or priors");
for (Int_t i = 0; i < nSpecies; i++) probabilities[i] = -1;
return;
}
for (Int_t i = 0; i < nSpecies; i++) {
probabilities[i] = fProbDensity[i] * prior[i] / sum;
}
}
void AliPID::GetProbabilities(Double_t* probabilities) const
{
GetProbabilities(probabilities, fgPrior);
}
AliPID::EParticleType AliPID::GetMostProbable(const Double_t* prior) const
{
Double_t max = 0.;
EParticleType id = kPion;
Int_t nSpecies = ((fCharged) ? kSPECIESC : kSPECIESCN);
for (Int_t i = 0; i < nSpecies; i++) {
Double_t prob = fProbDensity[i] * prior[i];
if (prob > max) {
max = prob;
id = EParticleType(i);
}
}
if (max == 0) {
AliError("Invalid probability densities or priors");
}
return id;
}
AliPID::EParticleType AliPID::GetMostProbable() const
{
return GetMostProbable(fgPrior);
}
void AliPID::SetPriors(const Double_t* prior, Bool_t charged)
{
Double_t sum = 0;
for (Int_t i = 0; i < kSPECIESCN; i++) {
if (charged && (i >= kSPECIESC)) {
fgPrior[i] = 0;
} else {
if (prior[i] < 0) {
AliWarningClass(Form("negative prior (%g) for %ss. "
"Using 0 instead.", prior[i],
fgkParticleName[i]));
fgPrior[i] = 0;
} else {
fgPrior[i] = prior[i];
}
}
sum += prior[i];
}
if (sum == 0) {
AliWarningClass("all priors are zero.");
}
}
void AliPID::SetPrior(EParticleType iType, Double_t prior)
{
if (prior < 0) {
AliWarningClass(Form("negative prior (%g) for %ss. Using 0 instead.",
prior, fgkParticleName[iType]));
prior = 0;
}
fgPrior[iType] = prior;
}
AliPID& AliPID::operator *= (const AliPID& pid)
{
for (Int_t i = 0; i < kSPECIESCN; i++) {
fProbDensity[i] *= pid.fProbDensity[i];
}
return *this;
}
AliPID operator * (const AliPID& pid1, const AliPID& pid2)
{
AliPID result;
result *= pid1;
result *= pid2;
return result;
}