#include <TClass.h>
#include <TArrayD.h>
#include <TH2.h>
#include <THnSparse.h>
#include <TString.h>
#include "AliLog.h"
#include "AliPID.h"
#include "AliPIDResponse.h"
#include "AliVParticle.h"
#include "AliHFEcollection.h"
#include "AliHFEpid.h"
#include "AliHFEpidBase.h"
#include "AliHFEpidQAmanager.h"
#include "AliHFEpidTPC.h"
#include "AliHFEpidTRD.h"
#include "AliHFEpidTOF.h"
#include "AliHFEtools.h"
#include "AliHFEtofPIDqa.h"
ClassImp(AliHFEpidTOF)
AliHFEtofPIDqa::AliHFEtofPIDqa():
AliHFEdetPIDqa()
, fHistos(NULL)
{
}
AliHFEtofPIDqa::AliHFEtofPIDqa(const char* name):
AliHFEdetPIDqa(name, "QA for TOF")
, fHistos(NULL)
{
}
AliHFEtofPIDqa::AliHFEtofPIDqa(const AliHFEtofPIDqa &o):
AliHFEdetPIDqa(o)
, fHistos()
{
o.Copy(*this);
}
AliHFEtofPIDqa &AliHFEtofPIDqa::operator=(const AliHFEtofPIDqa &o){
AliHFEdetPIDqa::operator=(o);
if(&o != this) o.Copy(*this);
return *this;
}
AliHFEtofPIDqa::~AliHFEtofPIDqa(){
if(fHistos) delete fHistos;
}
void AliHFEtofPIDqa::Copy(TObject &o) const {
AliHFEtofPIDqa &target = dynamic_cast<AliHFEtofPIDqa &>(o);
if(target.fHistos){
delete target.fHistos;
target.fHistos = NULL;
}
if(fHistos) target.fHistos = new AliHFEcollection(*fHistos);
}
Long64_t AliHFEtofPIDqa::Merge(TCollection *coll){
if(!coll) return 0;
if(coll->IsEmpty()) return 1;
TIter it(coll);
AliHFEtofPIDqa *refQA = NULL;
TObject *o = NULL;
Long64_t count = 0;
TList listHistos;
while((o = it())){
refQA = dynamic_cast<AliHFEtofPIDqa *>(o);
if(!refQA) continue;
listHistos.Add(refQA->fHistos);
count++;
}
fHistos->Merge(&listHistos);
return count + 1;
}
void AliHFEtofPIDqa::Initialize(){
fHistos = new AliHFEcollection("tofqahistos", "Collection of TOF QA histograms");
const Int_t kPIDbins = AliPID::kSPECIES + 1;
const Int_t kSteps = 2;
const Int_t kCentralityBins = 11;
const Double_t kMinPID = -1;
const Double_t kMinP = 0.;
const Double_t kMaxPID = (Double_t)AliPID::kSPECIES;
const Double_t kMaxP = 20.;
Int_t kPbins = fQAmanager->HasHighResolutionHistos() ? 1000 : 100;
Int_t kSigmaBins = fQAmanager->HasHighResolutionHistos() ? 1400 : 140;
Int_t nBinsSigma[5] = {kPIDbins, kPbins, kSigmaBins, kSteps, kCentralityBins};
Double_t minSigma[5] = {kMinPID, kMinP, -12., 0, 0};
Double_t maxSigma[5] = {kMaxPID, kMaxP, 12., 2., 11.};
fHistos->CreateTHnSparse("tofnSigma", "TOF signal; species; p [GeV/c]; TOF signal [a.u.]; Selection Step; Centrality", 5, nBinsSigma, minSigma, maxSigma);
fHistos->CreateTHnSparse("tofMonitorTPC", "TPC signal; species; p [GeV/c]; TPC signal [a.u.]; Selection Step; Centrality", 5, nBinsSigma, minSigma, maxSigma);
AliHFEpidTOF *tofpid = dynamic_cast<AliHFEpidTOF *>(fQAmanager->GetDetectorPID(AliHFEpid::kTOFpid));
if(!tofpid) AliDebug(1, "TOF PID not available");
if(tofpid && tofpid->IsGenerateTOFmismatch()){
AliDebug(1, "Prepare histogram for TOF mismatch tracks");
fHistos->CreateTHnSparse("tofMismatch", "TOF signal; species; p [GeV/c]; TOF signal [a.u.]; Selection Step; Centrality", 5, nBinsSigma, minSigma, maxSigma);
}
}
void AliHFEtofPIDqa::ProcessTrack(const AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step){
Float_t centrality = track->GetCentrality();
Int_t species = track->GetAbInitioPID();
if(species >= AliPID::kSPECIES) species = -1;
AliDebug(1, Form("Monitoring particle of type %d for step %d", species, step));
AliHFEpidTOF *tofpid= dynamic_cast<AliHFEpidTOF *>(fQAmanager->GetDetectorPID(AliHFEpid::kTOFpid));
const AliPIDResponse *pidResponse = tofpid ? tofpid->GetPIDResponse() : NULL;
Double_t contentSignal[5];
contentSignal[0] = species;
contentSignal[1] = track->GetRecTrack()->P();
contentSignal[2] = pidResponse ? pidResponse->NumberOfSigmasTOF(track->GetRecTrack(), AliPID::kElectron): 0.;
contentSignal[3] = step;
contentSignal[4] = centrality;
fHistos->Fill("tofnSigma", contentSignal);
contentSignal[2] = pidResponse ? pidResponse->NumberOfSigmasTPC(track->GetRecTrack(), AliPID::kElectron) : 0.;
fHistos->Fill("tofMonitorTPC", contentSignal);
if(tofpid && tofpid->IsGenerateTOFmismatch()){
AliDebug(1,"Filling TOF mismatch histos");
TArrayD nsigmismatch(tofpid->GetNmismatchTracks());
tofpid->GenerateTOFmismatch(track->GetRecTrack(),tofpid->GetNmismatchTracks(), nsigmismatch);
for(int itrk = 0; itrk < tofpid->GetNmismatchTracks(); itrk++){
contentSignal[2] = nsigmismatch[itrk];
fHistos->Fill("tofMismatch",contentSignal);
}
}
}
TH2 *AliHFEtofPIDqa::MakeTPCspectrumNsigma(AliHFEdetPIDqa::EStep_t step, Int_t species, Int_t centralityClass){
THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("tofMonitorTPC"));
if(!histo){
AliError("QA histogram monitoring TPC nSigma not available");
return NULL;
}
if(species > -1 && species < AliPID::kSPECIES){
histo->GetAxis(0)->SetRange(species + 2, species + 2);
}
TString centname, centtitle;
Bool_t hasCentralityInfo = kTRUE;
if(centralityClass > -1){
if(histo->GetNdimensions() < 5){
AliError("Centrality Information not available");
centname = centtitle = "MinBias";
hasCentralityInfo = kFALSE;
} else {
histo->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
centname = Form("Cent%d", centralityClass);
centtitle = Form("Centrality %d", centralityClass);
}
} else {
histo->GetAxis(4)->SetRange(1, histo->GetAxis(4)->GetNbins()-1);
centname = centtitle = "MinBias";
hasCentralityInfo = kTRUE;
}
histo->GetAxis(3)->SetRange(step + 1, step + 1);
TH2 *hSpec = histo->Projection(2, 1);
TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
TString histname = Form("hSigmaTPC%s%s%s", specID.Data(), stepname.Data(), centname.Data());
TString histtitle = Form("TPC Sigma for %s %s PID %s", speciesname.Data(), stepname.Data(), centtitle.Data());
hSpec->SetName(histname.Data());
hSpec->SetTitle(histtitle.Data());
histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
histo->GetAxis(3)->SetRange(0, histo->GetAxis(3)->GetNbins());
if(hasCentralityInfo)histo->GetAxis(4)->SetRange(0, histo->GetAxis(4)->GetNbins());
return hSpec;
}
TH2 *AliHFEtofPIDqa::MakeSpectrumNSigma(AliHFEdetPIDqa::EStep_t istep, Int_t species, Int_t centralityClass){
THnSparseF *hSignal = dynamic_cast<THnSparseF *>(fHistos->Get("tofnSigma"));
if(!hSignal) return NULL;
hSignal->GetAxis(3)->SetRange(istep + 1, istep + 1);
if(species >= 0 && species < AliPID::kSPECIES)
hSignal->GetAxis(0)->SetRange(2 + species, 2 + species);
TString centname, centtitle;
Bool_t hasCentralityInfo = kTRUE;
if(centralityClass > -1){
if(hSignal->GetNdimensions() < 5){
AliError("Centrality Information not available");
centname = centtitle = "MinBias";
hasCentralityInfo = kFALSE;
} else {
AliInfo(Form("Centrality Class: %d", centralityClass));
hSignal->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
centname = Form("Cent%d", centralityClass);
centtitle = Form("Centrality %d", centralityClass);
}
} else {
hSignal->GetAxis(4)->SetRange(1, hSignal->GetAxis(4)->GetNbins()-1);
centname = centtitle = "MinBias";
hasCentralityInfo = kTRUE;
}
TH2 *hTmp = hSignal->Projection(2,1);
TString hname = Form("hTOFsigma%s%s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after", centname.Data()),
htitle = Form("Time-of-flight spectrum[#sigma] %s selection %s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after", centtitle.Data());
if(species > -1){
hname += AliPID::ParticleName(species);
htitle += Form(" for %ss", AliPID::ParticleName(species));
}
hTmp->SetName(hname.Data());
hTmp->SetTitle(htitle.Data());
hTmp->SetStats(kFALSE);
hTmp->GetXaxis()->SetTitle("p [GeV/c]");
hTmp->GetYaxis()->SetTitle("TOF time|_{el} - expected time|_{el} [#sigma]");
hSignal->GetAxis(3)->SetRange(0, hSignal->GetAxis(3)->GetNbins());
hSignal->GetAxis(0)->SetRange(0, hSignal->GetAxis(0)->GetNbins());
if(hasCentralityInfo)hSignal->GetAxis(4)->SetRange(0, hSignal->GetAxis(4)->GetNbins());
return hTmp;
}
TH1 *AliHFEtofPIDqa::GetHistogram(const char *name){
if(!fHistos) return NULL;
TObject *histo = fHistos->Get(name);
if(!histo->InheritsFrom("TH1")) return NULL;
return dynamic_cast<TH1 *>(histo);
}