#include <TF1.h>
#include <TMath.h>
#include <TH2D.h>
#include "AliTPCdEdxInfo.h"
#include "AliAODPid.h"
#include "AliAODTrack.h"
#include "AliAODMCParticle.h"
#include "AliESDtrack.h"
#include "AliExternalTrackParam.h"
#include "AliLog.h"
#include "AliMCParticle.h"
#include "AliPID.h"
#include "AliPIDResponse.h"
#include "AliHFEpidTPC.h"
#include "AliHFEpidQAmanager.h"
ClassImp(AliHFEpidTPC)
AliHFEpidTPC::AliHFEpidTPC() :
AliHFEpidBase()
, fLineCrossingsEnabled(0)
, fkEtaCorrection(NULL)
, fkCentralityCorrection(NULL)
, fkEtaMeanCorrection(NULL)
, fkEtaWidthCorrection(NULL)
, fkCentralityMeanCorrection(NULL)
, fkCentralityWidthCorrection(NULL)
, fkCentralityEtaCorrectionMeanJpsi(NULL)
, fkCentralityEtaCorrectionWidthJpsi(NULL)
, fHasCutModel(kFALSE)
, fUseOnlyOROC(kFALSE)
, fNsigmaTPC(3)
, fRejectionEnabled(0)
, fUsedEdx(kFALSE)
{
memset(fkUpperSigmaCut, 0, sizeof(const TF1 *) * 12);
memset(fkLowerSigmaCut, 0, sizeof(const TF1 *) * 12);
memset(fRejection, 0, sizeof(Float_t) * 4 * AliPID::kSPECIES);
memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
memset(fPAsigCut, 0, sizeof(Float_t) * 2);
memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2);
}
AliHFEpidTPC::AliHFEpidTPC(const char* name) :
AliHFEpidBase(name)
, fLineCrossingsEnabled(0)
, fkEtaCorrection(NULL)
, fkCentralityCorrection(NULL)
, fkEtaMeanCorrection(NULL)
, fkEtaWidthCorrection(NULL)
, fkCentralityMeanCorrection(NULL)
, fkCentralityWidthCorrection(NULL)
, fkCentralityEtaCorrectionMeanJpsi(NULL)
, fkCentralityEtaCorrectionWidthJpsi(NULL)
, fHasCutModel(kFALSE)
, fUseOnlyOROC(kFALSE)
, fNsigmaTPC(3)
, fRejectionEnabled(0)
, fUsedEdx(kFALSE)
{
memset(fkUpperSigmaCut, 0, sizeof(const TF1 *) * 12);
memset(fkLowerSigmaCut, 0, sizeof(const TF1 *) * 12);
memset(fRejection, 0, sizeof(Float_t) * 4 * AliPID::kSPECIES);
memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
memset(fPAsigCut, 0, sizeof(Float_t) * 2);
memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2);
}
AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) :
AliHFEpidBase("")
, fLineCrossingsEnabled(0)
, fkEtaCorrection(NULL)
, fkCentralityCorrection(NULL)
, fkEtaMeanCorrection(NULL)
, fkEtaWidthCorrection(NULL)
, fkCentralityMeanCorrection(NULL)
, fkCentralityWidthCorrection(NULL)
, fkCentralityEtaCorrectionMeanJpsi(NULL)
, fkCentralityEtaCorrectionWidthJpsi(NULL)
, fHasCutModel(ref.fHasCutModel)
, fUseOnlyOROC(ref.fUseOnlyOROC)
, fNsigmaTPC(2)
, fRejectionEnabled(0)
, fUsedEdx(kFALSE)
{
ref.Copy(*this);
}
AliHFEpidTPC &AliHFEpidTPC::operator=(const AliHFEpidTPC &ref){
if(this != &ref){
ref.Copy(*this);
}
return *this;
}
void AliHFEpidTPC::Copy(TObject &o) const{
AliHFEpidTPC &target = dynamic_cast<AliHFEpidTPC &>(o);
target.fkEtaCorrection = fkEtaCorrection;
target.fkCentralityCorrection = fkCentralityCorrection;
target.fkEtaMeanCorrection = fkEtaMeanCorrection;
target.fkEtaWidthCorrection = fkEtaWidthCorrection;
target.fkCentralityMeanCorrection = fkCentralityMeanCorrection;
target.fkCentralityWidthCorrection = fkCentralityWidthCorrection;
target.fLineCrossingsEnabled = fLineCrossingsEnabled;
target.fHasCutModel = fHasCutModel;
target.fUseOnlyOROC = fUseOnlyOROC;
target.fNsigmaTPC = fNsigmaTPC;
target.fRejectionEnabled = fRejectionEnabled;
memcpy(target.fkUpperSigmaCut, fkUpperSigmaCut, sizeof(const TF1 *) * 12);
memcpy(target.fkLowerSigmaCut, fkLowerSigmaCut, sizeof(const TF1 *) * 12);
memcpy(target.fLineCrossingSigma, fLineCrossingSigma, sizeof(Double_t) * AliPID::kSPECIES);
memcpy(target.fPAsigCut, fPAsigCut, sizeof(Float_t) * 2);
memcpy(target.fNAsigmaTPC, fNAsigmaTPC, sizeof(Float_t) * 2);
AliHFEpidBase::Copy(target);
}
AliHFEpidTPC::~AliHFEpidTPC(){
}
Bool_t AliHFEpidTPC::InitializePID(Int_t ){
return kTRUE;
}
Int_t AliHFEpidTPC::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const
{
if(!fkPIDResponse) return 0;
AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
const AliVTrack *rectrack;
AliESDtrack esdtrack;
AliAODTrack aodtrack;
Double_t correctedTPCnSigma=0.;
Bool_t TPCnSigmaCorrected=kFALSE;
if((fkEtaMeanCorrection&&fkEtaWidthCorrection)||
(fkCentralityMeanCorrection&&fkCentralityWidthCorrection)){
TPCnSigmaCorrected=kTRUE;
correctedTPCnSigma=GetCorrectedTPCnSigma(track->GetRecTrack()->Eta(), track->GetMultiplicity(), fkPIDResponse->NumberOfSigmasTPC(track->GetRecTrack(), AliPID::kElectron));
}
if((fkCentralityEtaCorrectionMeanJpsi)&&
(fkCentralityEtaCorrectionWidthJpsi)){
TPCnSigmaCorrected=kTRUE;
correctedTPCnSigma=GetCorrectedTPCnSigmaJpsi(track->GetRecTrack()->Eta(), track->GetMultiplicity(), fkPIDResponse->NumberOfSigmasTPC(track->GetRecTrack(), AliPID::kElectron));
}
if(fkEtaCorrection || fkCentralityCorrection){
if(track->IsESDanalysis()){
esdtrack.~AliESDtrack();
new(&esdtrack) AliESDtrack(*(static_cast<const AliESDtrack *>(track->GetRecTrack())));
if(track->IsPbPb() && HasCentralityCorrection())
ApplyCentralityCorrection(&esdtrack, track->GetMultiplicity(), anatype);
if(HasEtaCorrection())
ApplyEtaCorrection(&esdtrack, anatype);
rectrack = &esdtrack;
} else {
aodtrack.~AliAODTrack();
new(&aodtrack) AliAODTrack(*(static_cast<const AliAODTrack *>(track->GetRecTrack())));
if(track->IsPbPb() && HasCentralityCorrection())
ApplyCentralityCorrection(&aodtrack, track->GetMultiplicity(), anatype);
if(HasEtaCorrection())
ApplyEtaCorrection(&aodtrack, anatype);
rectrack = &aodtrack;
}
} else {
rectrack = track->GetRecTrack();
}
AliHFEpidObject tpctrack(*track);
tpctrack.SetRecTrack(rectrack);
if(TPCnSigmaCorrected)tpctrack.SetCorrectedTPCnSigma(correctedTPCnSigma);
if(pidqa) pidqa->ProcessTrack(&tpctrack, AliHFEpid::kTPCpid, AliHFEdetPIDqa::kBeforePID);
AliDebug(1, "Doing TPC PID based on n-Sigma cut approach");
Float_t nsigma=correctedTPCnSigma;
if(!TPCnSigmaCorrected)
nsigma = fUsedEdx ? rectrack->GetTPCsignal() : fkPIDResponse->NumberOfSigmasTPC(rectrack, AliPID::kElectron);
AliDebug(1, Form("TPC NSigma: %f", nsigma));
Bool_t isLineCrossing = kFALSE;
for(Int_t ispecies = 0; ispecies < AliPID::kSPECIES; ispecies++){
if(ispecies == AliPID::kElectron) continue;
if(!(fLineCrossingsEnabled & 1 << ispecies)) continue;
if(TMath::Abs(fkPIDResponse->NumberOfSigmasTPC(rectrack, (AliPID::EParticleType)ispecies)) < fLineCrossingSigma[ispecies] && TMath::Abs(nsigma) < fNsigmaTPC){
isLineCrossing = kTRUE;
break;
}
}
if(isLineCrossing) return 0;
if(HasParticleRejection()){
Int_t reject = Reject(rectrack, anatype);
if(reject != 0) return reject;
}
Int_t pdg = 0;
if(fHasCutModel){
pdg = CutSigmaModel(&tpctrack) ? 11 : 0;
} else {
Float_t p = rectrack->P();
if(HasAsymmetricSigmaCut()) {
if(p >= fPAsigCut[0] && p <= fPAsigCut[1]) {
if(nsigma >= fNAsigmaTPC[0] && nsigma <= fNAsigmaTPC[1]) pdg = 11;
else pdg = 0;
}
else pdg = 0;
}
else {
if(TMath::Abs(nsigma) < fNsigmaTPC ) pdg = 11;
}
}
if(pidqa && pdg != 0) pidqa->ProcessTrack(&tpctrack, AliHFEpid::kTPCpid, AliHFEdetPIDqa::kAfterPID);
return pdg;
}
Bool_t AliHFEpidTPC::CutSigmaModel(const AliHFEpidObject * const track) const {
Bool_t isSelected = kTRUE;
AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
Float_t nsigma = fUsedEdx ? track->GetRecTrack()->GetTPCsignal() : fkPIDResponse->NumberOfSigmasTPC(track->GetRecTrack(), AliPID::kElectron);
Double_t p = GetP(track->GetRecTrack(), anatype);
Int_t centrality = track->IsPbPb() ? track->GetCentrality() + 1 : 0;
AliDebug(2, Form("Centrality: %d\n", centrality));
if(centrality > 11) return kFALSE;
const TF1 *cutfunction;
if((cutfunction = fkUpperSigmaCut[centrality]) && nsigma > cutfunction->Eval(p)) isSelected = kFALSE;
if((cutfunction = fkLowerSigmaCut[centrality]) && nsigma < cutfunction->Eval(p)) isSelected = kFALSE;
return isSelected;
}
Int_t AliHFEpidTPC::Reject(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const{
Int_t pdc[AliPID::kSPECIES] = {11,13,211,321,2212};
Double_t p = GetP(track, anaType);
for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
if(!TESTBIT(fRejectionEnabled, ispec)) continue;
if(p < fRejection[4*ispec] || p > fRejection[4*ispec+2]) continue;
Double_t sigma = fkPIDResponse->NumberOfSigmasTPC(track, static_cast<AliPID::EParticleType>(ispec));
if(sigma >= fRejection[4*ispec+1] && sigma <= fRejection[4*ispec+3]) return pdc[ispec] * track->Charge();
}
return 0;
}
void AliHFEpidTPC::ApplyEtaCorrection(AliVTrack *track, AliHFEpidObject::AnalysisType_t anatype) const{
AliDebug(1, Form("Applying correction function %s\n", fkEtaCorrection->GetName()));
Double_t original = track->GetTPCsignal(),
eta = track->Eta();
if(anatype == AliHFEpidObject::kESDanalysis){
AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
if(fkEtaCorrection->Eval(eta)>0.0) esdtrack->SetTPCsignal(original/fkEtaCorrection->Eval(eta), esdtrack->GetTPCsignalSigma(), esdtrack->GetTPCsignalN());
} else {
AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
AliAODPid *pid = aodtrack->GetDetPid();
if(pid && fkEtaCorrection->Eval(eta)>0.0) pid->SetTPCsignal(original/fkEtaCorrection->Eval(eta));
}
}
void AliHFEpidTPC::ApplyCentralityCorrection(AliVTrack *track, Double_t centralityEstimator, AliHFEpidObject::AnalysisType_t anatype) const{
AliDebug(1, Form("Applying correction function %s\n", fkCentralityCorrection->GetName()));
Double_t original = track->GetTPCsignal();
if(anatype == AliHFEpidObject::kESDanalysis){
AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
if(fkCentralityCorrection->Eval(centralityEstimator)>0.0) esdtrack->SetTPCsignal(original/fkCentralityCorrection->Eval(centralityEstimator), esdtrack->GetTPCsignalSigma(), esdtrack->GetTPCsignalN());
} else {
AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
AliAODPid *pid = aodtrack->GetDetPid();
if(pid && fkCentralityCorrection->Eval(centralityEstimator)>0.0) pid->SetTPCsignal(original/fkCentralityCorrection->Eval(centralityEstimator));
}
}
Double_t AliHFEpidTPC::GetCorrectedTPCnSigma(Double_t eta, Double_t centralityEstimator, Double_t tpcNsigma) const{
Double_t corrtpcNsigma = tpcNsigma;
if(fkEtaMeanCorrection&&fkEtaWidthCorrection){
if(TMath::Abs(fkEtaWidthCorrection->Eval(eta))>0.0000001) corrtpcNsigma=(corrtpcNsigma-fkEtaMeanCorrection->Eval(eta))/fkEtaWidthCorrection->Eval(eta);
}
if(fkCentralityMeanCorrection&&fkCentralityWidthCorrection) {
if(TMath::Abs(fkCentralityWidthCorrection->Eval(centralityEstimator))>0.0000001) corrtpcNsigma=(corrtpcNsigma-fkCentralityMeanCorrection->Eval(centralityEstimator))/fkCentralityWidthCorrection->Eval(centralityEstimator);
}
return corrtpcNsigma;
}
Double_t AliHFEpidTPC::GetCorrectedTPCnSigmaJpsi(Double_t eta, Double_t centralityEstimator, Double_t tpcNsigma) const{
Double_t corrtpcNsigma = tpcNsigma;
if(fkCentralityEtaCorrectionMeanJpsi&&fkCentralityEtaCorrectionWidthJpsi){
const TAxis *caxis = fkCentralityEtaCorrectionMeanJpsi->GetXaxis();
const TAxis *eaxis = fkCentralityEtaCorrectionMeanJpsi->GetYaxis();
Int_t cbin = caxis->FindFixBin(centralityEstimator);
Int_t ebin = eaxis->FindFixBin(eta);
Double_t center = fkCentralityEtaCorrectionMeanJpsi->GetBinContent(cbin,ebin);
Double_t width = fkCentralityEtaCorrectionWidthJpsi->GetBinContent(cbin,ebin);
if(TMath::Abs(width)>0.0000001) corrtpcNsigma=(corrtpcNsigma-center)/width;
}
return corrtpcNsigma;
}
void AliHFEpidTPC::UseOROC(AliVTrack *track, AliHFEpidObject::AnalysisType_t anatype) const{
if(anatype == AliHFEpidObject::kESDanalysis){
AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
AliTPCdEdxInfo *dEdxInfo = track->GetTPCdEdxInfo();
Double32_t TPCsignalRegion[4];
Char_t TPCsignalNRegion[3];
Char_t TPCsignalNRowRegion[3];
dEdxInfo->GetTPCSignalRegionInfo(TPCsignalRegion,TPCsignalNRegion,TPCsignalNRowRegion);
esdtrack->SetTPCsignal(TPCsignalRegion[3],esdtrack->GetTPCsignalSigma(),(TPCsignalNRegion[1]+TPCsignalNRegion[2]));
} else {
AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
AliTPCdEdxInfo *dEdxInfo = track->GetTPCdEdxInfo();
Double32_t TPCsignalRegion[4];
Char_t TPCsignalNRegion[3];
Char_t TPCsignalNRowRegion[3];
dEdxInfo->GetTPCSignalRegionInfo(TPCsignalRegion,TPCsignalNRegion,TPCsignalNRowRegion);
AliAODPid *pid = aodtrack->GetDetPid();
if(pid) pid->SetTPCsignal(TPCsignalRegion[3]);
if(pid) pid->SetTPCsignalN((TPCsignalNRegion[1]+TPCsignalNRegion[2]));
}
}
void AliHFEpidTPC::AddTPCdEdxLineCrossing(Int_t species, Double_t sigma){
if(species >= AliPID::kSPECIES){
AliError("Species doesn't exist");
return;
}
fLineCrossingsEnabled |= 1 << species;
fLineCrossingSigma[species] = sigma;
}
Double_t AliHFEpidTPC::GetP(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anatype) const {
Double_t p = -1;
if(anatype == AliHFEpidObject::kESDanalysis){
const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
if(esdtrack) p = esdtrack->GetInnerParam() ? esdtrack->GetInnerParam()->GetP() : esdtrack->P();
}
if(anatype == AliHFEpidObject::kAODanalysis)
{
const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
if(aodtrack) p = aodtrack->GetDetPid() ? aodtrack->GetDetPid()->GetTPCmomentum() : aodtrack->P();
}
return p;
}