#include <TMatrixF.h>
#include "TFormula.h"
#include "TBenchmark.h"
#include "TPrincipal.h"
#include "TFile.h"
#include "TSystem.h"
#include "AliPHOS.h"
#include "AliPHOSPIDv1.h"
#include "AliESDEvent.h"
#include "AliESDVertex.h"
#include "AliPHOSTrackSegment.h"
#include "AliPHOSEmcRecPoint.h"
#include "AliPHOSRecParticle.h"
ClassImp( AliPHOSPIDv1)
AliPHOSPIDv1::AliPHOSPIDv1() :
AliPHOSPID(),
fBayesian(kFALSE),
fDefaultInit(kFALSE),
fWrite(kFALSE),
fFileNamePrincipalPhoton(),
fFileNamePrincipalPi0(),
fFileNameParameters(),
fPrincipalPhoton(0),
fPrincipalPi0(0),
fX(0),
fPPhoton(0),
fPPi0(0),
fParameters(0),
fVtx(0.,0.,0.),
fTFphoton(0),
fTFpiong(0),
fTFkaong(0),
fTFkaonl(0),
fTFhhadrong(0),
fTFhhadronl(0),
fDFmuon(0),
fERecWeight(0),
fChargedNeutralThreshold(0.),
fTOFEnThreshold(0),
fDispEnThreshold(0),
fDispMultThreshold(0)
{
InitParameters() ;
fDefaultInit = kTRUE ;
}
AliPHOSPIDv1::AliPHOSPIDv1(const AliPHOSPIDv1 & pid ) :
AliPHOSPID(pid),
fBayesian(kFALSE),
fDefaultInit(kFALSE),
fWrite(kFALSE),
fFileNamePrincipalPhoton(),
fFileNamePrincipalPi0(),
fFileNameParameters(),
fPrincipalPhoton(0),
fPrincipalPi0(0),
fX(0),
fPPhoton(0),
fPPi0(0),
fParameters(0),
fVtx(0.,0.,0.),
fTFphoton(0),
fTFpiong(0),
fTFkaong(0),
fTFkaonl(0),
fTFhhadrong(0),
fTFhhadronl(0),
fDFmuon(0),
fERecWeight(0),
fChargedNeutralThreshold(0.),
fTOFEnThreshold(0),
fDispEnThreshold(0),
fDispMultThreshold(0)
{
InitParameters() ;
}
AliPHOSPIDv1::AliPHOSPIDv1(AliPHOSGeometry *geom):
AliPHOSPID(geom),
fBayesian(kFALSE),
fDefaultInit(kFALSE),
fWrite(kFALSE),
fFileNamePrincipalPhoton(),
fFileNamePrincipalPi0(),
fFileNameParameters(),
fPrincipalPhoton(0),
fPrincipalPi0(0),
fX(0),
fPPhoton(0),
fPPi0(0),
fParameters(0),
fVtx(0.,0.,0.),
fTFphoton(0),
fTFpiong(0),
fTFkaong(0),
fTFkaonl(0),
fTFhhadrong(0),
fTFhhadronl(0),
fDFmuon(0),
fERecWeight(0),
fChargedNeutralThreshold(0.),
fTOFEnThreshold(0),
fDispEnThreshold(0),
fDispMultThreshold(0)
{
InitParameters() ;
fDefaultInit = kFALSE ;
}
AliPHOSPIDv1::~AliPHOSPIDv1()
{
fPrincipalPhoton = 0;
fPrincipalPi0 = 0;
delete [] fX ;
delete [] fPPhoton ;
delete [] fPPi0 ;
delete fParameters;
delete fTFphoton;
delete fTFpiong;
delete fTFkaong;
delete fTFkaonl;
delete fTFhhadrong;
delete fTFhhadronl;
delete fDFmuon;
}
void AliPHOSPIDv1::InitParameters()
{
fWrite = kTRUE ;
fBayesian = kTRUE ;
SetParameters() ;
fTphoton[0] = 7.83E8 ;
fTphoton[1] = 1.55E-8 ;
fTphoton[2] = 5.09E-10 ;
fTFphoton = new TFormula("ToF response to photons" , "gaus") ;
fTFphoton->SetParameters( fTphoton[0], fTphoton[1], fTphoton[2]) ;
fTpiong[0] = 6.73E8 ;
fTpiong[1] = 1.58E-8 ;
fTpiong[2] = 5.87E-10 ;
fTFpiong = new TFormula("ToF response to pions" , "gaus") ;
fTFpiong->SetParameters( fTpiong[0], fTpiong[1], fTpiong[2]) ;
fTkaong[0] = 3.93E8 ;
fTkaong[1] = 1.64E-8 ;
fTkaong[2] = 6.07E-10 ;
fTFkaong = new TFormula("ToF response to kaon" , "gaus") ;
fTFkaong->SetParameters( fTkaong[0], fTkaong[1], fTkaong[2]) ;
fTkaonl[0] = 2.0E9 ;
fTkaonl[1] = 1.68E-8 ;
fTkaonl[2] = 4.10E-10 ;
fTFkaonl = new TFormula("ToF response to kaon" , "landau") ;
fTFkaonl->SetParameters( fTkaonl[0], fTkaonl[1], fTkaonl[2]) ;
fThhadrong[0] = 2.02E8 ;
fThhadrong[1] = 1.73E-8 ;
fThhadrong[2] = 9.52E-10 ;
fTFhhadrong = new TFormula("ToF response to heavy hadrons" , "gaus") ;
fTFhhadrong->SetParameters( fThhadrong[0], fThhadrong[1], fThhadrong[2]) ;
fThhadronl[0] = 1.10E9 ;
fThhadronl[1] = 1.74E-8 ;
fThhadronl[2] = 1.00E-9 ;
fTFhhadronl = new TFormula("ToF response to heavy hadrons" , "landau") ;
fTFhhadronl->SetParameters( fThhadronl[0], fThhadronl[1], fThhadronl[2]) ;
fDphoton[0] = 1.5 ; fDphoton[1] = 0.49 ; fDphoton[2] =-1.7E-2 ;
fDphoton[3] = 1.5 ; fDphoton[4] = 4.0E-2 ; fDphoton[5] = 0.21 ;
fDphoton[6] = 4.8E-2 ; fDphoton[7] =-0.12 ; fDphoton[8] = 0.27 ;
fDphoton[9] = 16.;
fDpi0[0] = 0.25 ; fDpi0[1] = 3.3E-2 ; fDpi0[2] =-1.0e-5 ;
fDpi0[3] = 1.50 ; fDpi0[4] = 398. ; fDpi0[5] = 12. ;
fDpi0[6] =-7.0E-2 ; fDpi0[7] =-524. ; fDpi0[8] = 22. ;
fDpi0[9] = 110.;
fDhadron[0] = 6.5 ; fDhadron[1] =-5.3 ; fDhadron[2] = 1.5 ;
fDhadron[3] = 3.8 ; fDhadron[4] = 0.23 ; fDhadron[5] =-1.2 ;
fDhadron[6] = 0.88 ; fDhadron[7] = 9.3E-2 ; fDhadron[8] =-0.51 ;
fDhadron[9] = 2.;
fDmuon[0] = 0.0631 ;
fDmuon[1] = 1.4 ;
fDmuon[2] = 0.0557 ;
fDFmuon = new TFormula("Shower shape response to muons" , "landau") ;
fDFmuon->SetParameters( fDmuon[0], fDmuon[1], fDmuon[2]) ;
fXelectron[0] =-1.6E-2 ; fXelectron[1] = 0.77 ; fXelectron[2] =-0.15 ;
fXelectron[3] = 0.35 ; fXelectron[4] = 0.25 ; fXelectron[5] = 4.12 ;
fXelectron[6] = 0.30 ; fXelectron[7] = 0.11 ; fXelectron[8] = 0.16 ;
fXelectron[9] = 3.;
fXcharged[0] = 0.14 ; fXcharged[1] =-3.0E-2 ; fXcharged[2] = 0 ;
fXcharged[3] = 1.4 ; fXcharged[4] =-9.3E-2 ; fXcharged[5] = 1.4 ;
fXcharged[6] = 5.7 ; fXcharged[7] = 0.27 ; fXcharged[8] =-1.8 ;
fXcharged[9] = 1.2;
fZelectron[0] = 0.49 ; fZelectron[1] = 0.53 ; fZelectron[2] =-9.8E-2 ;
fZelectron[3] = 2.8E-2 ; fZelectron[4] = 5.0E-2 ; fZelectron[5] =-8.2E-2 ;
fZelectron[6] = 0.25 ; fZelectron[7] =-1.7E-2 ; fZelectron[8] = 0.17 ;
fZelectron[9] = 3.;
fZcharged[0] = 0.46 ; fZcharged[1] =-0.65 ; fZcharged[2] = 0.52 ;
fZcharged[3] = 1.1E-2 ; fZcharged[4] = 0. ; fZcharged[5] = 0. ;
fZcharged[6] = 0.60 ; fZcharged[7] =-8.2E-2 ; fZcharged[8] = 0.45 ;
fZcharged[9] = 1.2;
fChargedNeutralThreshold = 1e-5;
fTOFEnThreshold = 2;
fDispEnThreshold = 0.5;
fDispMultThreshold = 3;
fERecWeightPar[0] = 0.32 ;
fERecWeightPar[1] = 3.8 ;
fERecWeightPar[2] = 5.4E-3 ;
fERecWeightPar[3] = 5.6E-2 ;
fERecWeight = new TFormula("Weight for hadrons" , "[0]*exp(-x*[1])+[2]*exp(-x*[3])") ;
fERecWeight ->SetParameters(fERecWeightPar[0],fERecWeightPar[1] ,fERecWeightPar[2] ,fERecWeightPar[3]) ;
for (Int_t i =0; i< AliPID::kSPECIESCN ; i++)
fInitPID[i] = 1.;
}
void AliPHOSPIDv1::TrackSegments2RecParticles(Option_t *option)
{
if(strstr(option,"tim"))
gBenchmark->Start("PHOSPID");
if(strstr(option,"print")) {
Print() ;
return ;
}
if(fTrackSegments &&
fTrackSegments->GetEntriesFast()) {
GetVertex() ;
MakeRecParticles() ;
if(fBayesian)
MakePID() ;
if(strstr(option,"deb"))
PrintRecParticles(option) ;
}
if(strstr(option,"deb"))
PrintRecParticles(option);
if(strstr(option,"tim")){
gBenchmark->Stop("PHOSPID");
AliInfo(Form("took %f seconds for PID",
gBenchmark->GetCpuTime("PHOSPID")));
}
}
Double_t AliPHOSPIDv1::GausF(Double_t x, Double_t y, Double_t * par)
{
if (x > par[9]) x = par[9];
Double_t cnt = par[0] + par[1] * x + par[2] * x * x ;
Double_t mean = par[4] / (x*x) + par[5] / x + par[3] ;
Double_t sigma = par[7] / (x*x) + par[8] / x + par[6] ;
if(TMath::Abs(sigma) > 1.e-10){
return cnt*TMath::Gaus(y,mean,sigma);
}
else
return 0.;
}
Double_t AliPHOSPIDv1::GausPol2(Double_t x, Double_t y, Double_t * par)
{
Double_t cnt = par[0] + par[1] * x + par[2] * x * x ;
Double_t mean = par[3] + par[4] * x + par[5] * x * x ;
Double_t sigma = par[6] + par[7] * x + par[8] * x * x ;
if(TMath::Abs(sigma) > 1.e-10){
return cnt*TMath::Gaus(y,mean,sigma);
}
else
return 0.;
}
const TString AliPHOSPIDv1::GetFileNamePrincipal(TString particle) const
{
particle.ToLower();
TString name;
if (particle=="photon")
name = fFileNamePrincipalPhoton ;
else if (particle=="pi0" )
name = fFileNamePrincipalPi0 ;
else
AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n",
particle.Data()));
return name;
}
Float_t AliPHOSPIDv1::GetParameterCalibration(Int_t i) const
{
Float_t param = 0.;
if (i>2 || i<0) {
AliError(Form("Invalid parameter number: %d",i));
} else
param = (*fParameters)(0,i);
return param;
}
Float_t AliPHOSPIDv1::GetParameterCpv2Emc(Int_t i, TString axis) const
{
Float_t param = 0.;
if(i>2 || i<0) {
AliError(Form("Invalid parameter number: %d",i));
} else {
axis.ToLower();
if (axis == "x")
param = (*fParameters)(1,i);
else if (axis == "z")
param = (*fParameters)(2,i);
else {
AliError(Form("Invalid axis name: %s",axis.Data()));
}
}
return param;
}
Float_t AliPHOSPIDv1::GetCpv2EmcDistanceCut(TString axis, Float_t e) const
{
axis.ToLower();
Float_t p[]={0.,0.,0.};
for (Int_t i=0; i<3; i++) p[i] = GetParameterCpv2Emc(i,axis);
Float_t sig = p[0] + TMath::Exp(p[1] - p[2]*e);
return sig;
}
Float_t AliPHOSPIDv1::GetEllipseParameter(TString particle, TString param, Float_t e) const
{
particle.ToLower();
param. ToLower();
Float_t p[4]={0.,0.,0.,0.};
Float_t value = 0.0;
for (Int_t i=0; i<4; i++) p[i] = GetParameterToCalculateEllipse(particle,param,i);
if (particle == "photon") {
if (param.Contains("a")) e = TMath::Min((Double_t)e,70.);
else if (param.Contains("b")) e = TMath::Min((Double_t)e,70.);
else if (param.Contains("x0")) e = TMath::Max((Double_t)e,1.1);
}
if (particle == "photon")
value = p[0]/TMath::Sqrt(e) + p[1]*e + p[2]*e*e + p[3];
else if (particle == "pi0")
value = p[0] + p[1]*e + p[2]*e*e;
return value;
}
Float_t AliPHOSPIDv1::GetParameterPhotonBoundary (Int_t i) const
{
Float_t param = 0;
if (i>3 || i<0) {
AliError(Form("Wrong parameter number: %d\n",i));
} else
param = (*fParameters)(14,i) ;
return param;
}
Float_t AliPHOSPIDv1::GetParameterPi0Boundary (Int_t i) const
{
Float_t param = 0;
if (i>2 || i<0) {
AliError(Form("Wrong parameter number: %d\n",i));
} else
param = (*fParameters)(15,i) ;
return param;
}
Float_t AliPHOSPIDv1::GetParameterTimeGate(Int_t i) const
{
Float_t param = 0.;
if(i>2 || i<0) {
AliError(Form("Invalid Efficiency-Purity choice %d",i));
} else
param = (*fParameters)(3,i) ;
return param;
}
Float_t AliPHOSPIDv1::GetParameterToCalculateEllipse(TString particle, TString param, Int_t i) const
{
particle.ToLower();
param. ToLower();
Int_t offset = -1;
if (particle == "photon")
offset=0;
else if (particle == "pi0")
offset=5;
else
AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n",
particle.Data()));
Int_t p= -1;
Float_t par = 0;
if (param.Contains("a")) p=4+offset;
else if(param.Contains("b")) p=5+offset;
else if(param.Contains("c")) p=6+offset;
else if(param.Contains("x0"))p=7+offset;
else if(param.Contains("y0"))p=8+offset;
if (i>4 || i<0) {
AliError(Form("No parameter with index %d", i)) ;
} else if (p==-1) {
AliError(Form("No parameter with name %s", param.Data() )) ;
} else
par = (*fParameters)(p,i) ;
return par;
}
Int_t AliPHOSPIDv1::GetCPVBit(AliPHOSTrackSegment * ts, Int_t effPur, Float_t e) const
{
if(effPur>2 || effPur<0)
AliError(Form("Invalid Efficiency-Purity choice %d",effPur));
Float_t sigX = GetCpv2EmcDistanceCut("X",e);
Float_t sigZ = GetCpv2EmcDistanceCut("Z",e);
Float_t deltaX = TMath::Abs(ts->GetCpvDistance("X"));
Float_t deltaZ = TMath::Abs(ts->GetCpvDistance("Z"));
if((deltaX>sigX*(effPur+1)) && (deltaZ>sigZ*(effPur+1)))
return 1;
else
return 0;
}
Int_t AliPHOSPIDv1::GetPrincipalBit(TString particle, const Double_t* p, Int_t effPur, Float_t e)const
{
particle.ToLower();
Int_t prinbit = 0 ;
Float_t a = GetEllipseParameter(particle,"a" , e);
Float_t b = GetEllipseParameter(particle,"b" , e);
Float_t c = GetEllipseParameter(particle,"c" , e);
Float_t x0 = GetEllipseParameter(particle,"x0", e);
Float_t y0 = GetEllipseParameter(particle,"y0", e);
Float_t r = TMath::Power((p[0] - x0)/a,2) +
TMath::Power((p[1] - y0)/b,2) +
c*(p[0] - x0)*(p[1] - y0)/(a*b) ;
if((effPur==2) && (r<1./2.)) prinbit= 1;
if((effPur==1) && (r<2. )) prinbit= 1;
if((effPur==0) && (r<9./2.)) prinbit= 1;
if(r<0)
AliError("Negative square?") ;
return prinbit;
}
Int_t AliPHOSPIDv1::GetHardPhotonBit(AliPHOSEmcRecPoint * emc) const
{
Float_t e = emc->GetEnergy();
if (e < 30.0) return 0;
Float_t m2x = emc->GetM2x();
Float_t m2xBoundary = GetParameterPhotonBoundary(0) *
TMath::Exp(-TMath::Power(e-GetParameterPhotonBoundary(1),2)/2.0/
TMath::Power(GetParameterPhotonBoundary(2),2)) +
GetParameterPhotonBoundary(3);
AliDebug(1, Form("E=%f, m2x=%f, boundary=%f", e,m2x,m2xBoundary));
if (m2x < m2xBoundary)
return 1;
else
return 0;
}
Int_t AliPHOSPIDv1::GetHardPi0Bit(AliPHOSEmcRecPoint * emc) const
{
Float_t e = emc->GetEnergy();
if (e < 30.0) return 0;
Float_t m2x = emc->GetM2x();
Float_t m2xBoundary = GetParameterPi0Boundary(0) +
e * GetParameterPi0Boundary(1);
AliDebug(1,Form("E=%f, m2x=%f, boundary=%f",e,m2x,m2xBoundary));
if (m2x > m2xBoundary)
return 1;
else
return 0;
}
TVector3 AliPHOSPIDv1::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSCpvRecPoint * )const
{
TVector3 local ;
emc->GetLocalPosition(local) ;
AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ;
Float_t para = 0.925 ;
Float_t parb = 6.52 ;
TVector3 vtxOld(0.,0.,0.) ;
TVector3 vInc ;
Float_t x=local.X() ;
Float_t z=local.Z() ;
phosgeom->GetIncidentVector(vtxOld,emc->GetPHOSMod(),x,z,vInc) ;
Float_t depthxOld = 0.;
Float_t depthzOld = 0.;
Float_t energy = emc->GetEnergy() ;
if (energy > 0 && vInc.Y()!=0.) {
depthxOld = ( para * TMath::Log(energy) + parb ) * vInc.X()/TMath::Abs(vInc.Y()) ;
depthzOld = ( para * TMath::Log(energy) + parb ) * vInc.Z()/TMath::Abs(vInc.Y()) ;
}
else{
AliError("Cluster with zero energy \n");
}
phosgeom->GetIncidentVector(fVtx,emc->GetPHOSMod(),x,z,vInc) ;
Float_t depthx = 0.;
Float_t depthz = 0.;
if (energy > 0 && vInc.Y()!=0.) {
depthx = ( para * TMath::Log(energy) + parb ) * vInc.X()/TMath::Abs(vInc.Y()) ;
depthz = ( para * TMath::Log(energy) + parb ) * vInc.Z()/TMath::Abs(vInc.Y()) ;
}
Double_t xd=x+(depthxOld-depthx) ;
Double_t zd=z+(depthzOld-depthz) ;
TVector3 dir(0,0,0) ;
phosgeom->Local2Global(emc->GetPHOSMod(),xd,zd,dir) ;
dir-=fVtx ;
dir.SetMag(1.) ;
return dir ;
}
Double_t AliPHOSPIDv1::LandauF(Double_t x, Double_t y, Double_t * par)
{
if (x > par[9]) x = par[9];
Double_t cnt = par[0] + par[1] * x + par[2] * x * x ;
Double_t mean = par[4] / (x*x) + par[5] / x + par[3] ;
Double_t sigma = par[7] / (x*x) + par[8] / x + par[6] ;
if(TMath::Abs(sigma) > 1.e-10){
return cnt*TMath::Landau(y,mean,sigma);
}
else
return 0.;
}
Double_t AliPHOSPIDv1::LandauPol2(Double_t x, Double_t y, Double_t * par)
{
Double_t cnt = par[2] * (x*x) + par[1] * x + par[0] ;
Double_t mean = par[5] * (x*x) + par[4] * x + par[3] ;
Double_t sigma = par[8] * (x*x) + par[7] * x + par[6] ;
if(TMath::Abs(sigma) > 1.e-10){
return cnt*TMath::Landau(y,mean,sigma);
}
else
return 0.;
}
void AliPHOSPIDv1::MakePID()
{
const Int_t kSPECIES = AliPID::kSPECIESCN ;
Int_t nparticles = fRecParticles->GetEntriesFast() ;
if ( !fEMCRecPoints || !fCPVRecPoints || !fTrackSegments ) {
AliFatal("RecPoints or TrackSegments not found !") ;
}
TIter next(fTrackSegments) ;
AliPHOSTrackSegment * ts ;
Int_t index = 0 ;
Double_t * stof[kSPECIES] ;
Double_t * sdp [kSPECIES] ;
Double_t * scpv[kSPECIES] ;
Double_t * sw [kSPECIES] ;
for (Int_t i =0; i< kSPECIES; i++){
stof[i] = new Double_t[nparticles] ;
sdp [i] = new Double_t[nparticles] ;
scpv[i] = new Double_t[nparticles] ;
sw [i] = new Double_t[nparticles] ;
}
while ( (ts = (AliPHOSTrackSegment *)next()) ) {
AliPHOSEmcRecPoint * emc = 0 ;
if(ts->GetEmcIndex()>=0)
emc = (AliPHOSEmcRecPoint *) fEMCRecPoints->At(ts->GetEmcIndex()) ;
if (!emc) {
AliFatal(Form("-> emc(%d)", ts->GetEmcIndex())) ;
}
Float_t en = emc->GetEnergy();
Double_t time = emc->GetTime() ;
for(Int_t iii=0; iii<kSPECIES; iii++)stof[iii][index]=0. ;
stof[AliPID::kPhoton][index] = 1.;
stof[AliPID::kElectron][index] = 1.;
stof[AliPID::kEleCon][index] = 1.;
stof[AliPID::kPion][index] = 1./3.;
stof[AliPID::kKaon][index] = 1./3.;
stof[AliPID::kProton][index] = 1./3.;
stof[AliPID::kNeutron][index] = 1./2.;
stof[AliPID::kKaon0][index] = 1./2.;
stof[AliPID::kMuon][index] = 1.;
if(en < fTOFEnThreshold) {
Double_t pTofPion = fTFpiong ->Eval(time) ;
Double_t pTofKaon = 0;
if(time < fTkaonl[1])
pTofKaon = fTFkaong ->Eval(time) ;
else
pTofKaon = fTFkaonl ->Eval(time) ;
Double_t pTofNucleon = 0;
if(time < fThhadronl[1])
pTofNucleon = fTFhhadrong ->Eval(time) ;
else
pTofNucleon = fTFhhadronl ->Eval(time) ;
Double_t pTofNeHadron = (pTofKaon + pTofNucleon)/2. ;
Double_t pTofChHadron = (pTofPion + pTofKaon + pTofNucleon)/3. ;
stof[AliPID::kPhoton][index] = fTFphoton ->Eval(time) ;
stof[AliPID::kEleCon][index] = stof[AliPID::kPhoton][index] ;
stof[AliPID::kMuon][index] = stof[AliPID::kPhoton][index] ;
stof[AliPID::kElectron][index] = pTofPion ;
stof[AliPID::kPion][index] = pTofChHadron ;
stof[AliPID::kKaon][index] = pTofChHadron ;
stof[AliPID::kProton][index] = pTofChHadron ;
stof[AliPID::kKaon0][index] = pTofNeHadron ;
stof[AliPID::kNeutron][index] = pTofNeHadron ;
}
Float_t dispersion = emc->GetDispersion();
for(Int_t iii=0; iii<kSPECIES; iii++)sdp[iii][index]=0. ;
sdp[AliPID::kPhoton][index] = 1. ;
sdp[AliPID::kElectron][index] = 1. ;
sdp[AliPID::kPion][index] = 1. ;
sdp[AliPID::kKaon][index] = 1. ;
sdp[AliPID::kProton][index] = 1. ;
sdp[AliPID::kNeutron][index] = 1. ;
sdp[AliPID::kEleCon][index] = 1. ;
sdp[AliPID::kKaon0][index] = 1. ;
sdp[AliPID::kMuon][index] = 1. ;
if(en > fDispEnThreshold && emc->GetMultiplicity() > fDispMultThreshold){
sdp[AliPID::kPhoton][index] = GausF(en , dispersion, fDphoton) ;
sdp[AliPID::kElectron][index] = sdp[AliPID::kPhoton][index] ;
sdp[AliPID::kPion][index] = LandauF(en , dispersion, fDhadron ) ;
sdp[AliPID::kKaon][index] = sdp[AliPID::kPion][index] ;
sdp[AliPID::kProton][index] = sdp[AliPID::kPion][index] ;
sdp[AliPID::kNeutron][index] = sdp[AliPID::kPion][index] ;
sdp[AliPID::kEleCon][index] = sdp[AliPID::kPhoton][index];
sdp[AliPID::kKaon0][index] = sdp[AliPID::kPion][index] ;
sdp[AliPID::kMuon][index] = fDFmuon ->Eval(dispersion) ;
}
Float_t x = TMath::Abs(ts->GetCpvDistance("X")) ;
Float_t z = ts->GetCpvDistance("Z") ;
Double_t pcpv = 0 ;
Double_t pcpvneutral = 0. ;
Double_t elprobx = GausF(en , x, fXelectron) ;
Double_t elprobz = GausF(en , z, fZelectron) ;
Double_t chprobx = GausF(en , x, fXcharged) ;
Double_t chprobz = GausF(en , z, fZcharged) ;
Double_t pcpvelectron = elprobx * elprobz;
Double_t pcpvcharged = chprobx * chprobz;
if(pcpvelectron >= pcpvcharged)
pcpv = pcpvelectron ;
else
pcpv = pcpvcharged ;
if(pcpv < fChargedNeutralThreshold)
{
pcpvneutral = 1. ;
pcpvcharged = 0. ;
pcpvelectron = 0. ;
}
for(Int_t iii=0; iii<kSPECIES; iii++)scpv[iii][index]=0. ;
scpv[AliPID::kPion][index] = pcpvcharged ;
scpv[AliPID::kKaon][index] = pcpvcharged ;
scpv[AliPID::kProton][index] = pcpvcharged ;
scpv[AliPID::kMuon][index] = pcpvelectron ;
scpv[AliPID::kElectron][index] = pcpvelectron ;
scpv[AliPID::kEleCon][index] = pcpvelectron ;
scpv[AliPID::kPhoton][index] = pcpvneutral ;
scpv[AliPID::kNeutron][index] = pcpvneutral ;
scpv[AliPID::kKaon0][index] = pcpvneutral ;
stof[AliPID::kPi0][index] = 0. ;
scpv[AliPID::kPi0][index] = 0. ;
sdp [AliPID::kPi0][index] = 0. ;
if(en > 30.){
stof[AliPID::kPi0][index] = stof[AliPID::kPhoton][index];
scpv[AliPID::kPi0][index] = pcpvneutral ;
if(emc->GetMultiplicity() > fDispMultThreshold)
sdp [AliPID::kPi0][index] = GausF(en , dispersion, fDpi0) ;
}
if(en > 0.5){
scpv[AliPID::kMuon][index] = 0 ;
stof[AliPID::kMuon][index] = 0 ;
sdp [AliPID::kMuon][index] = 0 ;
}
for(Int_t iii=0; iii<kSPECIES; iii++)sw[iii][index]=1. ;
Float_t weight = fERecWeight ->Eval(en) ;
sw[AliPID::kPhoton][index] = 1. ;
sw[AliPID::kElectron][index] = 1. ;
sw[AliPID::kPion][index] = weight ;
sw[AliPID::kKaon][index] = weight ;
sw[AliPID::kProton][index] = weight ;
sw[AliPID::kNeutron][index] = weight ;
sw[AliPID::kEleCon][index] = 1. ;
sw[AliPID::kKaon0][index] = weight ;
sw[AliPID::kMuon][index] = weight ;
sw[AliPID::kPi0][index] = 1. ;
index++;
}
for(index = 0 ; index < nparticles ; index ++) {
AliPHOSRecParticle * recpar = static_cast<AliPHOSRecParticle *>(fRecParticles->At(index));
if(recpar->IsEleCon()){
fInitPID[AliPID::kEleCon] = 1. ;
fInitPID[AliPID::kPhoton] = 0. ;
fInitPID[AliPID::kElectron] = 0. ;
}
else{
fInitPID[AliPID::kEleCon] = 0. ;
fInitPID[AliPID::kPhoton] = 1. ;
fInitPID[AliPID::kElectron] = 1. ;
}
Int_t jndex ;
Double_t wn = 0.0 ;
for (jndex = 0 ; jndex < kSPECIES ; jndex++)
wn += stof[jndex][index] * sdp[jndex][index] * scpv[jndex][index] *
sw[jndex][index] * fInitPID[jndex] ;
if (TMath::Abs(wn)>0)
for (jndex = 0 ; jndex < kSPECIES ; jndex++) {
recpar->SetPID(jndex, stof[jndex][index] * sdp[jndex][index] *
sw[jndex][index] * scpv[jndex][index] *
fInitPID[jndex] / wn) ;
}
}
for (Int_t i =0; i< kSPECIES; i++){
delete [] stof[i];
delete [] sdp [i];
delete [] scpv[i];
delete [] sw [i];
}
}
void AliPHOSPIDv1::MakeRecParticles()
{
if ( !fEMCRecPoints || !fCPVRecPoints || !fTrackSegments ) {
AliFatal("RecPoints or TrackSegments not found !") ;
}
fRecParticles->Clear();
TIter next(fTrackSegments) ;
AliPHOSTrackSegment * ts ;
Int_t index = 0 ;
AliPHOSRecParticle * rp ;
while ( (ts = (AliPHOSTrackSegment *)next()) ) {
new( (*fRecParticles)[index] ) AliPHOSRecParticle() ;
rp = (AliPHOSRecParticle *)fRecParticles->At(index) ;
rp->SetTrackSegment(index) ;
rp->SetIndexInList(index) ;
AliPHOSEmcRecPoint * emc = 0 ;
if(ts->GetEmcIndex()>=0)
emc = (AliPHOSEmcRecPoint *) fEMCRecPoints->At(ts->GetEmcIndex()) ;
AliPHOSCpvRecPoint * cpv = 0 ;
if(ts->GetCpvIndex()>=0)
cpv = (AliPHOSCpvRecPoint *) fCPVRecPoints->At(ts->GetCpvIndex()) ;
Int_t track = 0 ;
track = ts->GetTrackIndex() ;
if (!emc) {
AliFatal(Form("-> emc(%d)", ts->GetEmcIndex())) ;
}
Float_t e = emc->GetEnergy() ;
Float_t lambda[2]={0.,0.} ;
emc->GetElipsAxis(lambda) ;
if((lambda[0]>0.01) && (lambda[1]>0.01)){
Float_t spher = 0. ;
Float_t emaxdtotal = 0. ;
if((lambda[0]+lambda[1])!=0)
spher=TMath::Abs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]);
emaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy();
fX[0] = lambda[0] ;
fX[1] = lambda[1] ;
fX[2] = emc->GetDispersion() ;
fX[3] = spher ;
fX[4] = emc->GetMultiplicity() ;
fX[5] = emaxdtotal ;
fX[6] = emc->GetCoreEnergy() ;
fPrincipalPhoton->X2P(fX,fPPhoton);
fPrincipalPi0 ->X2P(fX,fPPi0);
}
else{
fPPhoton[0]=-100.0;
fPPhoton[1]=-100.0;
fPPi0[0] =-100.0;
fPPi0[1] =-100.0;
}
Float_t time = emc->GetTime() ;
rp->SetTof(time) ;
for(Int_t effPur = 0; effPur < 3 ; effPur++){
if(GetCPVBit(ts, effPur,e) == 1 ){
rp->SetPIDBit(effPur) ;
}
if(time< (*fParameters)(3,effPur))
rp->SetPIDBit(effPur+3) ;
if(GetPrincipalBit("photon",fPPhoton,effPur,e) == 1)
rp->SetPIDBit(effPur+6) ;
if(GetPrincipalBit("pi0" ,fPPi0 ,effPur,e) == 1)
rp->SetPIDBit(effPur+9) ;
}
if(GetHardPhotonBit(emc))
rp->SetPIDBit(12) ;
if(GetHardPi0Bit (emc))
rp->SetPIDBit(13) ;
if(track >= 0)
rp->SetPIDBit(14) ;
TVector3 dir = GetMomentumDirection(emc,cpv) ;
dir.SetMag(e) ;
rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),e) ;
rp->SetCalcMass(0);
rp->Name();
rp->SetProductionVertex(fVtx.X(),fVtx.Y(),fVtx.Z(),0);
rp->SetFirstMother(-1);
rp->SetLastMother(-1);
rp->SetFirstDaughter(-1);
rp->SetLastDaughter(-1);
rp->SetPolarisation(0,0,0);
AliPHOSTrackSegment * ts1 = static_cast<AliPHOSTrackSegment *>(fTrackSegments->At(rp->GetPHOSTSIndex()));
AliPHOSEmcRecPoint * erp = static_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(ts1->GetEmcIndex()));
TVector3 pos ;
fGeom->GetGlobalPHOS(erp, pos) ;
rp->SetPos(pos);
index++ ;
}
}
void AliPHOSPIDv1::Print(const Option_t *) const
{
AliInfo("=============== AliPHOSPIDv1 ================") ;
printf("Making PID\n") ;
printf(" Pricipal analysis file from 0.5 to 100 %s\n", fFileNamePrincipalPhoton.Data() ) ;
printf(" Name of parameters file %s\n", fFileNameParameters.Data() ) ;
printf(" Matrix of Parameters: 14x4\n") ;
printf(" Energy Calibration 1x3 [3 parametres to calibrate energy: A + B* E + C * E^2]\n") ;
printf(" RCPV 2x3 rows x and z, columns function cut parameters\n") ;
printf(" TOF 1x3 [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]\n") ;
printf(" PCA 5x4 [5 ellipse parametres and 4 parametres to calculate them: A/Sqrt(E) + B* E + C * E^2 + D]\n") ;
printf(" Pi0 PCA 5x3 [5 ellipse parametres and 3 parametres to calculate them: A + B* E + C * E^2]\n") ;
fParameters->Print() ;
}
void AliPHOSPIDv1::PrintRecParticles(Option_t * option)
{
TString message ;
message = " found " ;
message += fRecParticles->GetEntriesFast();
message += " RecParticles\n" ;
if(strstr(option,"all")) {
message += "\n PARTICLE Index \n" ;
Int_t index ;
for (index = 0 ; index < fRecParticles->GetEntries() ; index++) {
AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) fRecParticles->At(index) ;
message += "\n" ;
message += rp->Name().Data() ;
message += " " ;
message += rp->GetIndexInList() ;
message += " " ;
message += rp->GetType() ;
}
}
AliInfo(message.Data() ) ;
}
void AliPHOSPIDv1::SetParameters()
{
fX = new double[7];
fPPhoton = new double[7];
fPPi0 = new double[7];
fFileNamePrincipalPhoton = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root" ;
TFile f( fFileNamePrincipalPhoton.Data(), "read" ) ;
fPrincipalPhoton = dynamic_cast<TPrincipal*> (f.Get("principal")) ;
f.Close() ;
fFileNamePrincipalPi0 = "$ALICE_ROOT/PHOS/PCA_pi0_40-120.root" ;
TFile fPi0( fFileNamePrincipalPi0.Data(), "read" ) ;
fPrincipalPi0 = dynamic_cast<TPrincipal*> (fPi0.Get("principal")) ;
fPi0.Close() ;
fFileNameParameters = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters.dat");
fParameters = new TMatrixF(16,4) ;
const Int_t kMaxLeng=255;
char string[kMaxLeng];
FILE *fd = fopen(fFileNameParameters.Data(),"r");
if (!fd)
AliFatal(Form("File %s with a PID parameters cannot be opened\n",
fFileNameParameters.Data()));
Int_t i=0;
while (fgets(string,kMaxLeng,fd) != NULL) {
if (string[0] == '\n' ) continue;
if (string[0] == '!' ) continue;
sscanf(string, "%f %f %f %f",
&(*fParameters)(i,0), &(*fParameters)(i,1),
&(*fParameters)(i,2), &(*fParameters)(i,3));
i++;
AliDebug(1, Form("Line %d: %s",i,string));
}
fclose(fd);
}
void AliPHOSPIDv1::SetParameterCalibration(Int_t i,Float_t param)
{
if(i>2 || i<0) {
AliError(Form("Invalid parameter number: %d",i));
} else
(*fParameters)(0,i) = param ;
}
void AliPHOSPIDv1::SetParameterCpv2Emc(Int_t i, TString axis, Float_t cut)
{
if(i>2 || i<0) {
AliError(Form("Invalid parameter number: %d",i));
} else {
axis.ToLower();
if (axis == "x") (*fParameters)(1,i) = cut;
else if (axis == "z") (*fParameters)(2,i) = cut;
else {
AliError(Form("Invalid axis name: %s",axis.Data()));
}
}
}
void AliPHOSPIDv1::SetParameterPhotonBoundary(Int_t i,Float_t param)
{
if(i>4 || i<0) {
AliError(Form("Invalid parameter number: %d",i));
} else
(*fParameters)(14,i) = param ;
}
void AliPHOSPIDv1::SetParameterPi0Boundary(Int_t i,Float_t param)
{
if(i>1 || i<0) {
AliError(Form("Invalid parameter number: %d",i));
} else
(*fParameters)(15,i) = param ;
}
void AliPHOSPIDv1::SetParameterTimeGate(Int_t i, Float_t gate)
{
if (i>2 || i<0) {
AliError(Form("Invalid Efficiency-Purity choice %d",i));
} else
(*fParameters)(3,i)= gate ;
}
void AliPHOSPIDv1::SetParameterToCalculateEllipse(TString particle, TString param, Int_t i, Float_t par)
{
particle.ToLower();
param. ToLower();
Int_t p= -1;
Int_t offset=0;
if (particle == "photon") offset=0;
else if (particle == "pi0") offset=5;
else
AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n",
particle.Data()));
if (param.Contains("a")) p=4+offset;
else if(param.Contains("b")) p=5+offset;
else if(param.Contains("c")) p=6+offset;
else if(param.Contains("x0"))p=7+offset;
else if(param.Contains("y0"))p=8+offset;
if((i>4)||(i<0)) {
AliError(Form("No parameter with index %d", i)) ;
} else if(p==-1) {
AliError(Form("No parameter with name %s", param.Data() )) ;
} else
(*fParameters)(p,i) = par ;
}
void AliPHOSPIDv1::GetVertex(void)
{
if(fESD){
const AliESDVertex *esdVtx = fESD->GetVertex() ;
if(esdVtx && esdVtx->GetChi2()!=0.){
fVtx.SetXYZ(esdVtx->GetX(),esdVtx->GetY(),esdVtx->GetZ()) ;
return ;
}
}
AliWarning("Can not read vertex from data, use fixed \n") ;
fVtx.SetXYZ(0.,0.,0.) ;
}
void AliPHOSPIDv1::SetInitPID(const Double_t *p) {
for (Int_t i=0; i<AliPID::kSPECIESCN; i++) fInitPID[i] = p[i];
}
void AliPHOSPIDv1::GetInitPID(Double_t *p) const {
for (Int_t i=0; i<AliPID::kSPECIESCN; i++) p[i] = fInitPID[i];
}