#include <TAxis.h>
#include <TBrowser.h>
#include <TH2.h>
#include <THnSparse.h>
#include <TString.h>
#include "AliESDtrack.h"
#include "AliLog.h"
#include "AliPID.h"
#include "AliPIDResponse.h"
#include "AliHFEcollection.h"
#include "AliHFEpidBase.h"
#include "AliHFEpidQAmanager.h"
#include "AliHFEpidTPC.h"
#include "AliHFEpidTRD.h"
#include "AliHFEtrdPIDqaV1.h"
ClassImp(AliHFEtrdPIDqaV1)
AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1():
AliHFEdetPIDqa(),
fHistos(NULL)
{
}
AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1(const Char_t *name):
AliHFEdetPIDqa(name, "QA for TRD"),
fHistos(NULL)
{
}
AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1(const AliHFEtrdPIDqaV1 &o):
AliHFEdetPIDqa(o),
fHistos(NULL)
{
}
AliHFEtrdPIDqaV1 &AliHFEtrdPIDqaV1::operator=(const AliHFEtrdPIDqaV1 &o){
if(this == &o) return *this;
AliHFEdetPIDqa::operator=(o);
fHistos = o.fHistos;
return *this;
}
AliHFEtrdPIDqaV1::~AliHFEtrdPIDqaV1(){
if(fHistos) delete fHistos;
}
Long64_t AliHFEtrdPIDqaV1::Merge(TCollection *coll){
if(!coll) return 0;
if(coll->IsEmpty()) return 1;
TIter it(coll);
AliHFEtrdPIDqaV1 *refQA = NULL;
TObject *o = NULL;
Long64_t count = 0;
TList listHistos;
while((o = it())){
refQA = dynamic_cast<AliHFEtrdPIDqaV1 *>(o);
if(!refQA) continue;
listHistos.Add(refQA->fHistos);
count++;
}
fHistos->Merge(&listHistos);
return count + 1;
}
void AliHFEtrdPIDqaV1::Browse(TBrowser *b){
if(b){
if(fHistos){
b->Add(fHistos, fHistos->GetName());
TString specnames[4] = {"All", "Electrons", "Pions", "Protons"};
Int_t specind[4] = {-1, AliPID::kElectron, AliPID::kPion, AliPID::kProton};
TList *listTM = new TList;
listTM->SetOwner();
TList *listLike = new TList;
listLike->SetOwner();
TList *listCharge = new TList;
listCharge->SetOwner();
TList *listTPCnsigma = new TList;
listTPCnsigma->SetOwner();
TH2 *hptr = NULL;
TList *lctm, *lclike, *lccharge, *lcTPCsig;
for(Int_t icent = -1; icent < 11; icent++){
lctm = new TList;
lctm->SetName(icent < 0 ? "MinBias" : Form("Centrality%d", icent));
lctm->SetOwner();
listTM->Add(lctm);
lclike = new TList;
lclike->SetName(icent < 0 ? "MinBias" : Form("Centrality%d", icent));
lclike->SetOwner();
listLike->Add(lclike);
lccharge = new TList;
lccharge->SetName(icent < 0 ? "MinBias" : Form("Centrality%d", icent));
lccharge->SetOwner();
listCharge->Add(lccharge);
lcTPCsig = new TList;
lcTPCsig->SetName(icent < 0 ? "MinBias" : Form("Centrality%d", icent));
lcTPCsig->SetOwner();
listTPCnsigma->Add(lcTPCsig);
for(Int_t ispec = 0; ispec < 4; ispec++){
for(Int_t istep = 0; istep < 2; istep++){
hptr = MakeTRDspectrumTM(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], icent);
if(hptr){
hptr->SetName(Form("hTRDtm%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
lctm->Add(hptr);
}
hptr = MakeTRDlikelihoodDistribution(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], icent);
hptr->SetName(Form("hTRDlike%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
lclike->Add(hptr);
hptr = MakeTRDchargeDistribution(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], icent);
hptr->SetName(Form("hTRDcharge%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
lccharge->Add(hptr);
hptr = MakeTPCspectrumNsigma(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], icent);
hptr->SetName(Form("hTPCspectrum%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
lcTPCsig->Add(hptr);
}
}
}
b->Add(listTM, "Projections Truncated Mean");
b->Add(listLike, "Projections Likelihood distribution");
b->Add(listCharge, "Projections Tracklet Charge");
b->Add(listTPCnsigma, "Projections TPC spectra");
}
}
}
void AliHFEtrdPIDqaV1::Initialize(){
AliDebug(1, "Initializing PID QA for TRD");
const Int_t kPIDbins = AliPID::kSPECIES + 1;
const Int_t kSteps = 2;
const Double_t kMinPID = -1;
const Double_t kMinP = 0.;
const Double_t kMaxPID = (Double_t)AliPID::kSPECIES;
const Double_t kMaxP = 20.;
const Int_t kCentralityBins = 11;
Int_t kPbins = fQAmanager->HasHighResolutionHistos() ? 1000 : 100;
Int_t tpcSigmaBins = fQAmanager->HasHighResolutionHistos() ? 1400 : 140;
Int_t trdLikelihoodBins = fQAmanager->HasHighResolutionHistos() ? 200 : 100;
fHistos = new AliHFEcollection("trdqahistos", "Collection of TRD QA histograms");
Int_t nBinsTPCSigma[5] = {kPIDbins, kPbins, tpcSigmaBins, kSteps, kCentralityBins};
Double_t minTPCSigma[5] = {kMinPID, kMinP, -12., 0., 0.};
Double_t maxTPCSigma[5] = {kMaxPID, kMaxP, 12., 2., 11.};
fHistos->CreateTHnSparse("hTPCsigma", "TPC sigma; species p [GeV/c]; TPC dEdx - <dE/dx>|_{el} [#sigma]; selection step", 5, nBinsTPCSigma, minTPCSigma, maxTPCSigma);
fHistos->Sumw2("hTPCsigma");
Int_t nBinsTRDlike[5] = {kPIDbins, kPbins, trdLikelihoodBins, kSteps, kCentralityBins};
Double_t minTRDlike[5] = {kMinPID, kMinP, 0., 0., 0.};
Double_t maxTRDlike[5] = {kMaxPID, kMaxP, 1., 2., 11.};
fHistos->CreateTHnSparse("hTRDlikelihood", "TRD Likelihood Distribution; species; p [GeV/c]; TRD electron Likelihood; selection step", 5, nBinsTRDlike, minTRDlike, maxTRDlike);
fHistos->Sumw2("hTRDlikelihood");
const Int_t kTRDchargeBins = 1000;
Int_t nBinsTRDcharge[5] = {kPIDbins, kPbins, kTRDchargeBins, kSteps, kCentralityBins};
Double_t minTRDcharge[5] = {kMinPID, kMinP, 0., 0., 0.};
Double_t maxTRDcharge[5] = {kMaxPID, kMaxP, 100000., 2., 11.};
fHistos->CreateTHnSparse("hTRDcharge", "Total TRD charge; species; p [GeV/c]; TRD charge [a.u.]; selection step", 5, nBinsTRDcharge, minTRDcharge, maxTRDcharge);
fHistos->Sumw2("hTRDcharge");
const Int_t kTRDtmBins = 1000;
Int_t nBinsTRDtm[5] = {kPIDbins, kPbins, kTRDtmBins, kSteps, kCentralityBins};
Double_t minTRDtm[5] = {kMinPID, kMinP, 0., 0., 0.};
Double_t maxTRDtm[5] = {kMaxPID, kMaxP, 20000., 2., 11.};
fHistos->CreateTHnSparse("hTRDtruncatedMean", "TRD truncated Mean; species; p [GeV/c]; TRD signal [a.u.]; selection step", 5, nBinsTRDtm, minTRDtm, maxTRDtm);
fHistos->Sumw2("hTRDtruncatedMean");
fHistos->CreateTH2F("hNtrackletsBefore", "Number of tracklets before PID; species; p (GeV/c), Number of Tracklets", kPbins, kMinP, kMaxP, 7, 0., 7.);
fHistos->Sumw2("hNtrackletsBefore");
fHistos->CreateTH2F("hNtrackletsAfter", "Number of tracklets after PID; species; p (GeV/c), Number of Tracklets", kPbins, kMinP, kMaxP, 7, 0., 7.);
fHistos->Sumw2("hNtrackletsAfter");
}
void AliHFEtrdPIDqaV1::ProcessTrack(const AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step){
AliDebug(1, Form("QA started for TRD PID for step %d", (Int_t)step));
Int_t species = track->GetAbInitioPID();
if(species >= AliPID::kSPECIES) species = -1;
AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
AliHFEpidTRD *trdpid = dynamic_cast<AliHFEpidTRD *>(fQAmanager->GetDetectorPID(AliHFEpid::kTRDpid));
const AliPIDResponse *pidResponse = trdpid ? trdpid->GetPIDResponse() : NULL;
Double_t container[5];
container[0] = species;
container[1] = trdpid ? trdpid->GetP(track->GetRecTrack(), anatype) : 0.;
container[2] = pidResponse ? pidResponse->NumberOfSigmasTPC(track->GetRecTrack(), AliPID::kElectron) : 0.;
container[3] = step;
container[4] = track->GetCentrality();
AliDebug(1, Form("Species %d, p %f, number of sigma TPC %f, step %f, centrality %f", (Int_t)species,(Float_t) container[1],(Float_t) container[2],(Float_t) container[3],(Float_t) container[4]));
fHistos->Fill("hTPCsigma", container);
AliDebug(1, "Filled TPC sigma\n");
container[2] = trdpid ? trdpid->GetElectronLikelihood(static_cast<const AliVTrack*>(track->GetRecTrack()), anatype) : 0;
fHistos->Fill("hTRDlikelihood", container);
AliDebug(1, "Filled likelihood\n");
if(track->IsESDanalysis()){
const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track->GetRecTrack());
if(esdtrack){
container[2] = trdpid ? trdpid->GetTRDSignalV1(esdtrack) : 0;
fHistos->Fill("hTRDtruncatedMean", container);
}
}
for(UInt_t ily = 0; ily < 6; ily++){
container[2] = trdpid ? trdpid->GetChargeLayer(track->GetRecTrack(), ily, anatype) : 0;
if(container[2] < 1e-3) continue;
fHistos->Fill("hTRDcharge", container);
AliDebug(1, "Filled TRD charge\n");
}
Int_t ntracklets = track->GetRecTrack()->GetTRDntrackletsPID();
if(step == AliHFEdetPIDqa::kBeforePID)
fHistos->Fill("hNtrackletsBefore", trdpid ? trdpid->GetP(track->GetRecTrack(), anatype) : 0., ntracklets);
else
fHistos->Fill("hNtrackletsAfter", trdpid ? trdpid->GetP(track->GetRecTrack(), anatype) : 0., ntracklets);
AliDebug(1, "Filled tracklet\n");
}
TH2 *AliHFEtrdPIDqaV1::MakeTPCspectrumNsigma(AliHFEdetPIDqa::EStep_t step, Int_t species, Int_t centralityClass){
THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTPCsigma"));
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 *AliHFEtrdPIDqaV1::MakeTRDspectrumTM(AliHFEdetPIDqa::EStep_t step, Int_t species, Int_t centralityClass){
THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDtruncatedMean"));
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= kFALSE;
}
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("hTMTRD%s%s%s", specID.Data(), stepname.Data(), centname.Data());
TString histtitle = Form("TRD Truncated Mean 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(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
if(hasCentralityInfo)histo->GetAxis(4)->SetRange(0, histo->GetAxis(4)->GetNbins());
return hSpec;
}
TH2 *AliHFEtrdPIDqaV1::MakeTRDlikelihoodDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species, Int_t centralityClass){
THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDlikelihood"));
if(!histo){
AliError("QA histogram monitoring TRD Electron Likelihood 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("hLikeElTRD%s%s%s", specID.Data(), stepname.Data(),centname.Data());
TString histtitle = Form("TRD electron Likelihood 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(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
if(hasCentralityInfo)histo->GetAxis(4)->SetRange(0, histo->GetAxis(4)->GetNbins());
return hSpec;
}
TH2 *AliHFEtrdPIDqaV1::MakeTRDchargeDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species, Int_t centralityClass){
THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDcharge"));
if(!histo){
AliError("QA histogram monitoring TRD total charge 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("hChargeTRD%s%s%s", specID.Data(), stepname.Data(), centname.Data());
TString histtitle = Form("TRD total charge 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(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
if(hasCentralityInfo)histo->GetAxis(4)->SetRange(0, histo->GetAxis(4)->GetNbins());
return hSpec;
}