#include <TAxis.h>
#include <TClass.h>
#include <TF1.h>
#include <TIterator.h>
#include <TList.h>
#include <TMath.h>
#include <TObjArray.h>
#include <TObjString.h>
#include <TString.h>
#include "AliLog.h"
#include "AliPID.h"
#include "AliVParticle.h"
#include "AliHFEcontainer.h"
#include "AliHFEpid.h"
#include "AliHFEpidQAmanager.h"
#include "AliHFEpidITS.h"
#include "AliHFEpidTPC.h"
#include "AliHFEpidTRD.h"
#include "AliHFEpidTOF.h"
#include "AliHFEpidEMCAL.h"
#include "AliHFEpidMC.h"
#include "AliHFEpidBayes.h"
#include "AliHFEvarManager.h"
ClassImp(AliHFEpid)
const Char_t* AliHFEpid::fgkDetectorName[AliHFEpid::kNdetectorPID + 1] = {
"MCPID",
"BAYESPID",
"ITSPID",
"TPCPID",
"TRDPID",
"TOFPID",
"EMCALPID",
"UndefinedPID"
};
AliHFEpid::AliHFEpid():
TNamed(),
fEnabledDetectors(0),
fNPIDdetectors(0),
fVarManager(NULL),
fCommonObjects(NULL)
{
memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
memset(fDetectorOrder, kUndefined, sizeof(UInt_t) * kNdetectorPID);
memset(fSortedOrder, 0, sizeof(UInt_t) * kNdetectorPID);
}
AliHFEpid::AliHFEpid(const Char_t *name):
TNamed(name, ""),
fEnabledDetectors(0),
fNPIDdetectors(0),
fVarManager(NULL),
fCommonObjects(NULL)
{
memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
memset(fDetectorOrder, kUndefined, sizeof(UInt_t) * kNdetectorPID);
memset(fSortedOrder, 0, sizeof(UInt_t) * kNdetectorPID);
fDetectorPID[kMCpid] = new AliHFEpidMC("MCPID");
fDetectorPID[kBAYESpid] = new AliHFEpidBayes("BAYESPID");
fDetectorPID[kTPCpid] = new AliHFEpidTPC("TPCPID");
fDetectorPID[kTRDpid] = new AliHFEpidTRD("TRDPID");
fDetectorPID[kITSpid] = new AliHFEpidITS("ITSPID");
fDetectorPID[kTOFpid] = new AliHFEpidTOF("TOFPID");
fDetectorPID[kEMCALpid] = new AliHFEpidEMCAL("EMCALPID");
}
AliHFEpid::AliHFEpid(const AliHFEpid &c):
TNamed(c),
fEnabledDetectors(c.fEnabledDetectors),
fNPIDdetectors(c.fNPIDdetectors),
fVarManager(c.fVarManager),
fCommonObjects(NULL)
{
memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
c.Copy(*this);
}
AliHFEpid& AliHFEpid::operator=(const AliHFEpid &c){
if(&c != this) c.Copy(*this);
return *this;
}
AliHFEpid::~AliHFEpid(){
for(Int_t idet = 0; idet < kNdetectorPID; idet++)
if(fDetectorPID[idet]) delete fDetectorPID[idet];
ClearCommonObjects();
}
void AliHFEpid::Copy(TObject &o) const{
TNamed::Copy(o);
AliHFEpid &target = dynamic_cast<AliHFEpid &>(o);
target.ClearCommonObjects();
target.fEnabledDetectors = fEnabledDetectors;
target.fNPIDdetectors = fNPIDdetectors;
target.fVarManager = fVarManager;
for(Int_t idet = 0; idet < kNdetectorPID; idet++){
if(target.fDetectorPID[idet])
delete target.fDetectorPID[idet];
if(fDetectorPID[idet])
target.fDetectorPID[idet] = dynamic_cast<AliHFEpidBase *>(fDetectorPID[idet]->Clone());
}
memcpy(target.fDetectorOrder, fDetectorOrder, sizeof(UInt_t) * kNdetectorPID);
memcpy(target.fSortedOrder, fSortedOrder, sizeof(UInt_t) * kNdetectorPID);
}
void AliHFEpid::AddCommonObject(TObject * const o){
if(!fCommonObjects) fCommonObjects = new TObjArray;
fCommonObjects->Add(o);
}
void AliHFEpid::ClearCommonObjects(){
if(fCommonObjects){
fCommonObjects->Delete();
delete fCommonObjects;
fCommonObjects = NULL;
}
}
void AliHFEpid::SetDetectorsForAnalysis(TString detectors){
TObjArray *detarray = detectors.Tokenize(",");
TObjString *detString(NULL);
int ndet(0);
for(int idet = 0; idet < detarray->GetEntries(); idet++){
detString = dynamic_cast<TObjString *>(detarray->At(idet));
if(detString) AddDetector(detString->String(), ndet++);
}
AliDebug(1, Form("%d detectors used in Analysis", ndet));
}
void AliHFEpid::AddDetector(TString detector, UInt_t position){
UInt_t detectorID = kUndefined;
detector.ToUpper();
if(!detector.CompareTo("MC")) detectorID = kMCpid;
else if(!detector.CompareTo("BAYES")) detectorID = kBAYESpid;
else if(!detector.CompareTo("TPC")) detectorID = kTPCpid;
else if(!detector.CompareTo("TRD")) detectorID = kTRDpid;
else if(!detector.CompareTo("TOF")) detectorID = kTOFpid;
else if(!detector.CompareTo("ITS")) detectorID = kITSpid;
else if(!detector.CompareTo("EMCAL")) detectorID = kEMCALpid;
else AliError("Detector not available");
if(detectorID == kUndefined) return;
if(IsDetectorOn(detectorID)) return;
SwitchOnDetector(detectorID);
fDetectorOrder[detectorID] = position;
fNPIDdetectors++;
}
Bool_t AliHFEpid::InitializePID(Int_t run){
if(!TestBit(kDetectorsSorted)) SortDetectors();
Bool_t status = kTRUE;
for(Int_t idet = 0; idet < kNdetectorPID; idet++){
if(!IsDetectorOn(idet)) continue;
if(fDetectorPID[idet]){
status &= fDetectorPID[idet]->InitializePID(run);
if(HasMCData() && status) fDetectorPID[idet]->SetHasMCData();
}
}
SetBit(kIsInit);
PrintStatus();
return status;
}
Bool_t AliHFEpid::IsSelected(const AliHFEpidObject * const track, AliHFEcontainer *cont, const Char_t *contname, AliHFEpidQAmanager *pidqa){
Bool_t isSelected = kTRUE;
AliDebug(1, Form("Particle used for PID, QA available: %s", pidqa ? "Yes" : "No"));
for(UInt_t idet = 0; idet < fNPIDdetectors; idet++){
AliDebug(2, Form("Using Detector %s\n", SortedDetectorName(idet)));
if(TMath::Abs(fDetectorPID[fSortedOrder[idet]]->IsSelected(track, pidqa)) != 11){
isSelected = kFALSE;
break;
}
AliDebug(2, "Particlae selected by detector");
if(fVarManager && cont){
TString reccontname = contname; reccontname += "Reco";
AliDebug(2, Form("Filling container %s", reccontname.Data()));
if(fVarManager->IsSignalTrack())
fVarManager->FillContainerStepname(cont, reccontname.Data(), SortedDetectorName(idet));
if(HasMCData()){
TString mccontname = contname; mccontname += "MC";
AliDebug(2, Form("MC Information available, Filling container %s", mccontname.Data()));
if(fVarManager->IsSignalTrack()) {
fVarManager->FillContainerStepname(cont, mccontname.Data(), SortedDetectorName(idet), kTRUE);
if(cont->GetCorrelationMatrix("correlationstepafterTOF")){
TString tstept("TOFPID");
if(!tstept.CompareTo(SortedDetectorName(idet))) {
fVarManager->FillCorrelationMatrix(cont->GetCorrelationMatrix("correlationstepafterTOF"));
}
}
}
}
}
}
return isSelected;
}
void AliHFEpid::SortDetectors(){
if(TestBit(kDetectorsSorted)) return;
TMath::Sort(static_cast<UInt_t>(kNdetectorPID), fDetectorOrder, fSortedOrder, kFALSE);
SetBit(kDetectorsSorted);
}
void AliHFEpid::SetPIDResponse(const AliPIDResponse * const pid){
for(Int_t idet = 0; idet < kNdetectorPID; idet++){
if(fDetectorPID[idet]) fDetectorPID[idet]->SetPIDResponse(pid);
}
}
const AliPIDResponse *AliHFEpid::GetPIDResponse() const {
const AliPIDResponse *response = NULL;
for(Int_t idet = 0; idet < kNdetectorPID; idet++){
if(fDetectorPID[idet]){
response = fDetectorPID[idet]->GetPIDResponse();
break;
}
}
return response;
}
void AliHFEpid::ConfigureTPCasymmetric(Double_t pmin, Double_t pmax, Double_t sigmamin, Double_t sigmamax){
AliHFEpidTPC *pid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
if(pid){
pid->SetTPCnSigma(3);
pid->SetAsymmetricTPCsigmaCut(pmin, pmax, sigmamin, sigmamax);
}
}
void AliHFEpid::ConfigureTPCrejectionSimple(){
AliHFEpidTPC *pid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
if(pid){
pid->SetTPCnSigma(3);
pid->SetRejectParticle(AliPID::kPion, 0., -100., 10., 1.);
}
}
void AliHFEpid::ConfigureTOF(Float_t TOFCut){
AliHFEpidTOF *tofpid = dynamic_cast<AliHFEpidTOF *>(fDetectorPID[kTOFpid]);
if(tofpid) tofpid->SetTOFnSigma(TOFCut);
}
void AliHFEpid::ConfigureTPCcentralityCut(Int_t centralityBin, const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut){
ConfigureTPCcut(centralityBin, lowerCutParam, params, upperTPCCut);
}
void AliHFEpid::ConfigureTPCdefaultCut(const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut){
ConfigureTPCcut(-1, lowerCutParam, params, upperTPCCut);
}
void AliHFEpid::ConfigureTPCcut(Int_t centralityBin, const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut){
if(HasMCData()) AliInfo("Configuring TPC for MC\n");
AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
TF1 *upperCut = new TF1(Form("upperCut%s", centralityBin < 0 ? "Default" : Form("Bin%d", centralityBin)), "[0]", 0, 20);
TF1 *lowerCut = new TF1(Form("lowerCut%s", centralityBin < 0 ? "Default" : Form("Bin%d", centralityBin)), lowerCutParam == NULL ? "[0] * TMath::Exp([1]*x) + [2]": lowerCutParam, 0, 20);
upperCut->SetParameter(0, upperTPCCut);
if(params){
for(Int_t ipar = 0; ipar < lowerCut->GetNpar(); ipar++)
{
lowerCut->SetParameter(ipar, params[ipar]);
}
} else {
if(HasMCData()) lowerCut->SetParameter(0, -2.5);
else lowerCut->SetParameter(0, -4.03);
lowerCut->SetParameter(1, -0.22);
if(HasMCData()) lowerCut->SetParameter(2, -2.2);
else lowerCut->SetParameter(2, 0.92);
}
if(tpcpid){
tpcpid->SetTPCnSigma(2);
if(centralityBin < 0){
tpcpid->SetUpperSigmaCutDefault(upperCut);
tpcpid->SetLowerSigmaCutDefault(lowerCut);
} else {
tpcpid->SetUpperSigmaCutCentrality(upperCut, centralityBin);
tpcpid->SetLowerSigmaCutCentrality(lowerCut, centralityBin);
}
}
AddCommonObject(upperCut);
AddCommonObject(lowerCut);
}
void AliHFEpid::ConfigureBayesDetectorMask(Int_t detmask){
if(HasMCData()) AliInfo("Configuring Bayes for MC\n");
AliHFEpidBayes *bayespid = dynamic_cast<AliHFEpidBayes *>(fDetectorPID[kBAYESpid]);
if(bayespid)
{
bayespid->SetBayesDetectorMask(detmask);
printf("detector mask in pid class %i \n",detmask);
}
}
void AliHFEpid::ConfigureBayesPIDThreshold(Float_t pidthres){
if(HasMCData()) AliInfo("Configuring Bayes for MC\n");
AliHFEpidBayes *bayespid = dynamic_cast<AliHFEpidBayes *>(fDetectorPID[kBAYESpid]);
if(bayespid)
{
bayespid->SetBayesPIDThreshold(pidthres);
printf("combined pidthreshold %f \n",pidthres);
}
}
void AliHFEpid::PrintStatus() const {
printf("\n%s: Printing configuration\n", GetName());
printf("===============================================\n");
printf("PID Detectors: \n");
Int_t npid = 0;
TString detectors[kNdetectorPID] = {"MC", "BAYES", "ITS", "TPC", "TRD", "TOF", "EMCAL"};
for(Int_t idet = 0; idet < kNdetectorPID; idet++){
if(IsDetectorOn(idet)){
printf("\t%s\n", detectors[idet].Data());
npid++;
}
}
if(!npid) printf("\tNone\n");
printf("\n");
}