#include <TMath.h>
#include <TObjArray.h>
#include <TPDGCode.h>
#include <TString.h>
#include <TMultiLayerPerceptron.h>
#include <TFile.h>
#include "AliAODMCParticle.h"
#include "AliAODEvent.h"
#include "AliAODTrack.h"
#include "AliESDtrack.h"
#include "AliESDEvent.h"
#include "AliMCEvent.h"
#include "AliMCParticle.h"
#include "AliMultiplicity.h"
#include "AliPID.h"
#include "AliPIDResponse.h"
#include "AliVParticle.h"
#include "AliHFEcollection.h"
#include "AliHFEpidQA.h"
#include "AliHFEV0info.h"
#include "AliHFEV0pid.h"
#include "AliHFEV0pidMC.h"
#include "AliHFEV0cuts.h"
#include "AliHFEtrdPIDqa.h"
ClassImp(AliHFEpidQA)
AliHFEpidQA::AliHFEpidQA():
fEvent(NULL),
fMC(NULL),
fV0pid(NULL),
fV0pidMC(NULL),
fTRDpidQA(NULL),
fOutput(NULL),
fESDpid(NULL),
fNNref(NULL),
fTRDTotalChargeInSlice0(kFALSE),
fIsPbPb(kFALSE),
fIsppMultiBin(kFALSE)
{
for(Int_t mom = 0; mom < 11; mom++){
fNet[mom] = NULL;
}
}
AliHFEpidQA::AliHFEpidQA(const AliHFEpidQA &ref):
TObject(ref),
fEvent(NULL),
fMC(NULL),
fV0pid(NULL),
fV0pidMC(NULL),
fTRDpidQA(NULL),
fOutput(NULL),
fESDpid(NULL),
fNNref(NULL),
fTRDTotalChargeInSlice0(ref.fTRDTotalChargeInSlice0),
fIsPbPb(kFALSE),
fIsppMultiBin(kFALSE)
{
for(Int_t mom = 0; mom < 11; mom++){
fNet[mom] = NULL;
}
ref.Copy(*this);
}
AliHFEpidQA &AliHFEpidQA::operator=(const AliHFEpidQA &ref){
if(this != &ref)
ref.Copy(*this);
return *this;
}
AliHFEpidQA::~AliHFEpidQA(){
}
void AliHFEpidQA::Copy(TObject &o) const {
TObject::Copy(o);
AliHFEpidQA &target = dynamic_cast<AliHFEpidQA &>(o);
target.fMC = fMC;
target.fTRDTotalChargeInSlice0 = fTRDTotalChargeInSlice0;
if(target.fESDpid) delete target.fESDpid;
target.fESDpid = fESDpid;
if(target.fV0pid) delete target.fV0pid;
if(fV0pid)
target.fV0pid = dynamic_cast<AliHFEV0pid *>(fV0pid->Clone());
else
target.fV0pid = NULL;
if(target.fV0pidMC) delete target.fV0pidMC;
if(fV0pidMC)
target.fV0pidMC = dynamic_cast<AliHFEV0pidMC *>(fV0pidMC->Clone());
else
target.fV0pidMC = NULL;
if(target.fTRDpidQA) delete target.fTRDpidQA;
if(fTRDpidQA)
target.fTRDpidQA = dynamic_cast<AliHFEtrdPIDqa *>(fTRDpidQA->Clone());
else
target.fTRDpidQA = NULL;
if(target.fOutput) delete target.fOutput;
if(fOutput)
target.fOutput = dynamic_cast<AliHFEcollection *>(fOutput->Clone());
}
void AliHFEpidQA::Init(){
if(fNNref){
for(Int_t mom = 0; mom < 11; mom++){
fNet[mom] = (TMultiLayerPerceptron*) fNNref->Get(Form("NN_Mom%d", mom));
if(!fNet[mom]){
AliError(Form("No reference network for momentum bin %d!", mom));
}
}
}
fV0pid = new AliHFEV0pid("fV0pid");
if(HasV0pidQA()) fV0pid->InitQA();
fV0pidMC = new AliHFEV0pidMC();
fV0pidMC->Init();
fTRDpidQA = new AliHFEtrdPIDqa;
if(fTRDTotalChargeInSlice0) fTRDpidQA->SetTotalChargeInSlice0();
fTRDpidQA->Init();
fOutput = new AliHFEcollection("pidQA", "PID QA output");
const char *name[4] = {"Electron", "PionK0", "PionL", "Proton"};
const char *title[4] = {"Electron", "K0 Pion", "Lambda Pion", "Proton"};
const char *det[4] = {"ITS", "TPC", "TRD", "TOF"};
const int effs[6] = {70, 75, 80, 85, 90, 95};
for(Int_t i = 0; i < 4; i++){
fOutput->CreateTH2F(Form("purity%s", name[i]), Form("%s Purity", title[i]), 2, -0.5, 1.5, 20, 0.1, 10, 1);
for(Int_t idet = 0; idet < 4; idet++){
if(idet != 2){
fOutput->CreateTH2F(Form("h%s_nSigma_%s", det[idet], name[i]), Form("%s number of sigmas for %ss; p (GeV/c); number of sigmas", det[idet], title[i]), 20, 0.1, 10, 100, -7, 7, 0);
}
fOutput->CreateTH2F(Form("h%s_PID_p_%s", det[idet], name[i]), Form("%s PID for %ss; p (GeV/c); ITS PID", det[idet], title[i]), 100, 0.1, 10, 5, -0.5, 4.5, 0);
fOutput->CreateTH2F(Form("h%s_El_like_%s", det[idet], name[i]), Form("%s Electron Likelihoods for %ss; p (GeV/c); likelihood", det[idet], title[i]), 25, 0.1, 10, 1000, 0., 1., 0);
fOutput->CreateTH2F(Form("h%s_El_like_MC_%s", det[idet], name[i]), Form("%s Electron Likelihoods for MC %ss; p (GeV/c); likelihood", det[idet], title[i]), 25, 0.1, 10, 1000, 0., 1., 0);
}
fOutput->CreateTH2F(Form("hITS_Signal_%s", name[i]), Form("ITS Signal org. for %ss", title[i]), 40, 0.1, 10, 400, 0, 1000, 0);
fOutput->CreateTH2Fvector1(5, Form("hITS_dEdx_%s", name[i]), Form("ITS dEdx for %ss; p (GeV/c); dEdx (a.u.)", title[i]), 40, 0.1, 10, 400, 0, 1000, 0);
fOutput->CreateTH2F(Form("hTPC_dEdx_%s", name[i]), Form("TPC dEdx for %ss; p (GeV/c); dEdx (a.u.)", title[i]), 20, 0.1, 10, 200, 0, 200, 0);
fOutput->CreateTH2F(Form("hTRD_trk_%s", name[i]), Form("%s PID tracklets; p (GeV/c); N TRD tracklets", title[i]), 100, 0.1, 10, 7, -0.5, 6.5, 0);
fOutput->CreateTH2F(Form("hTRD_Nslc_%s", name[i]), Form("%s PID slices > 0; p (GeV/c); N slc > 0", title[i]), 100, 0.1, 10, 9, -0.5, 8.5, 0);
fOutput->CreateTH2F(Form("hTRD_slc_%s", name[i]), Form("%s PID slices > 0 position; p (GeV/c); slice", title[i]), 100, 0.1, 10, 9, -0.5, 8.5, 0);
fOutput->CreateTH2F(Form("hTRD_cls_%s", name[i]), Form("TRD clusters for %s candidates; p (GeV/c); N cls", title[i]), 25, 0.1, 10, 1000, 0, 1000);
fOutput->CreateTH2F(Form("hTRD_dEdx_%s", name[i]), Form("TRD dEdx (trk) for %s candidates; p (GeV/c); tracklet dEdx (a.u.)", title[i]), 25, 0.1, 10, 1000, 0, 100000, 0);
for(Int_t itl = 4; itl < 6; itl++){
for(Int_t ieff = 0; ieff < 6; ieff++){
fOutput->CreateTH2F(Form("hTRD_likeSel_%s_%dtls_%deff", name[i], itl, effs[ieff]), Form(" %s selected as electrons with %d tracklets and %f electron efficiency", title[i], itl, static_cast<Double_t>(effs[ieff])/100.), 44, 0.1, 20., 100, 0., 1.);
fOutput->CreateTH2F(Form("hTRD_likeRej_%s_%dtls_%deff", name[i], itl, effs[ieff]), Form(" %s rejected as electrons with %d tracklets and %f electron efficiency", title[i], itl, static_cast<Double_t>(effs[ieff])/100.), 44, 0.1, 20., 100, 0., 1.);
}
}
fOutput->CreateTH2F(Form("hTOF_beta_%s", name[i]), Form("TOF beta for %s; p (GeV/c); #beta", title[i]), 50, 0.1, 5, 350, 0.4, 1.1, 0);
}
fOutput->CreateTH1F("hITS_dEdx_nSamples", "ITS - number of non 0 dEdx samples", 4, 0.5, 4.5);
fOutput->CreateTH2F("hTPC_PID", "TPC pid all tracks; tpc pid probability; species",100, 0, 1, 5, -0.5, 4.5 );
fOutput->CreateTH2F("hTOF_PID", "TOF pid all tracks; tof pid probability; species",100, 0, 1,5, -0.5, 4.5 );
fOutput->CreateTH2F("hTOF_beta_all", "TOF beta for all nice single tracks; p (GeV/c); #beta", 100, 0.1, 10, 480, 0, 1.2, 0);
fOutput->CreateTH2F("hTOF_qa_TmT0mT", "TOF (t - t0 - t[pion]) qa verus run", 10000, 114000, 124000, 200, -200, 200);
fOutput->CreateTH2F("hTOF_qa_length", "TOF track length verus run", 10000, 114000, 112400, 200, 0, 1000);
fOutput->CreateTH2F("hITS_kFlags", "ITS flags; flag; V0 candidates", 5, 0.5, 5.5, 5, -0.5, 4.5);
fOutput->CreateTH2F("hTPC_kFlags", "TPC flags; flag; V0 candidates", 5, 0.5, 5.5, 5, -0.5, 4.5);
fOutput->CreateTH2F("hTRD_kFlags", "TRD flags; flag; V0 candidates", 5, 0.5, 5.5, 5, -0.5, 4.5);
fOutput->CreateTH2F("hTOF_kFlags", "TOF flags; flag; V0 candidates", 5, 0.5, 5.5, 5, -0.5, 4.5);
Int_t nT0[2] = {10000, 100};
Double_t minT0[2] = {114500, -1000};
Double_t maxT0[2] = {124500, 1000};
fOutput->CreateTHnSparse("hEvent_T0", "T0 as a function of run number; run number; T0 (ns)", 2, nT0, minT0, maxT0);
fOutput->CreateTH2F("h_tender_check_01", "tender -vs- HFEpidQA; HFEpidQA V0 candiadates; tender V0 candidates", 4, -1.5, 2.5, 4, -1.5, 2.5);
fOutput->CreateTH2Fvector1(3, "h_tender_check_02", "tender -vs- HFEpidQA per type; AliHFEpidQA V0 ; tender V0", 4, -1.5, 2.5, 4, -1.5, 2.5);
{
Int_t nbins[6] = {4, 40, 16, 72, 2, 2};
Double_t min[6] = { 0, 0.1, -0.8, 0., 0., 0.};
Double_t max[6] = { 4., 20., 0.8, 2*TMath::Pi(), 2., 2.};
fOutput->CreateTHnSparse("hIllumination", "Illumination", 6, nbins, min, max);
fOutput->BinLogAxis("hIllumination", 1);
}
{
Int_t nbins[6] = { 5, 40, 20, 100, 100, 200};
Double_t min[6] = { -0.5, 0.1, 0., 0., -5., 0.};
Double_t max[6] = { 4.5, 20., 200, 1., 5., 120.};
TString htitle = "TPC info sparse; VO identified species; p (GeV/c); n TPC clusters; TPC N sigma; TPC signal";
fOutput->CreateTHnSparse("hTPCclusters",htitle,6, nbins, min, max);
fOutput->BinLogAxis("hTPCclusters", 1);
}
{
Int_t nbins[7] = {5, 20, 6, 7, 9, 45, 100};
Double_t min[7] = {-0.5, 0.5, -0.5, -0.5, -0.5, -0.5, 0.};
Double_t max[7] = {4.5, 10, 5.5, 6.5, 8.5, 179.5, 1.};
TString htitle = "TRD tracklets; VO identified species; p (GeV/c); tracklet position; No. PID tacklets; No. slices; No. clusters; El. likelihood";
fOutput->CreateTHnSparse("hTRDtracklets",htitle,7, nbins, min, max);
fOutput->BinLogAxis("hTRDtracklets", 1);
}
}
void AliHFEpidQA::Process(){
if(!fV0pid){
AliError("V0pid not available! Forgotten to initialize?");
return;
}
if(!fESDpid){
AliError("fESDpid not initialized, I am leaving this!");
return;
}
if(!fEvent){
AliError("AliVEvent not available, returning");
}
if(fMC) fV0pidMC->SetMCEvent(fMC);
if(fMC) fV0pid->SetMCEvent(fMC);
fV0pid->Process(fEvent);
TObjArray *hfeelectrons = fV0pid->GetListOfElectrons();
TObjArray *hfepionsK0 = fV0pid->GetListOfPionsK0();
TObjArray *hfepionsL = fV0pid->GetListOfPionsL();
TObjArray *hfeprotons = fV0pid->GetListOfProtons();
TObjArray *electrons = MakeTrackList(hfeelectrons);
TObjArray *pionsK0 = MakeTrackList(hfepionsK0);
TObjArray *pionsL = MakeTrackList(hfepionsL);
TObjArray *protons = MakeTrackList(hfeprotons);
TObjArray *cleanElectrons = MakeCleanListElectrons(hfeelectrons);
if(fMC){
fV0pidMC->Process(electrons, AliHFEV0cuts::kRecoElectron);
fV0pidMC->Process(pionsK0, AliHFEV0cuts::kRecoPionK0);
fV0pidMC->Process(pionsL, AliHFEV0cuts::kRecoPionL);
fV0pidMC->Process(protons, AliHFEV0cuts::kRecoProton);
}
AliDebug(2, Form("Number of Electrons : %d", electrons->GetEntries()));
AliDebug(2, Form("Number of K0 Pions : %d", pionsK0->GetEntries()));
AliDebug(2, Form("Number of Lambda Pions : %d", pionsL->GetEntries()));
AliDebug(2, Form("Number of Protons : %d", protons->GetEntries()));
if(fMC){
AliDebug(2, "MC Information available. Doing Purity checks...");
MakePurity(electrons, AliHFEV0cuts::kRecoElectron);
MakePurity(pionsK0, AliHFEV0cuts::kRecoPionK0);
MakePurity(pionsL, AliHFEV0cuts::kRecoPionL);
MakePurity(protons, AliHFEV0cuts::kRecoProton);
}
CheckEvent();
FillIllumination(electrons, AliHFEV0cuts::kRecoElectron);
FillIllumination(pionsK0, AliHFEV0cuts::kRecoPionK0);
FillIllumination(pionsL, AliHFEV0cuts::kRecoPionL);
FillIllumination(protons, AliHFEV0cuts::kRecoProton);
Int_t centrality=-1;
if(fIsPbPb) centrality=GetCentrality(fEvent);
if(fIsppMultiBin) centrality=GetMultiplicityITS(fEvent);
fTRDpidQA->SetCentralityBin(centrality);
fTRDpidQA->ProcessTracks(cleanElectrons, AliPID::kElectron);
fTRDpidQA->ProcessTracks(pionsK0, AliPID::kPion);
fTRDpidQA->ProcessTracks(pionsL, AliPID::kPion);
fTRDpidQA->ProcessTracks(protons, AliPID::kProton);
FillElectronLikelihoods(electrons, AliHFEV0cuts::kRecoElectron);
FillElectronLikelihoods(pionsK0, AliHFEV0cuts::kRecoPionK0);
FillElectronLikelihoods(pionsL, AliHFEV0cuts::kRecoPionL);
FillElectronLikelihoods(protons, AliHFEV0cuts::kRecoProton);
FillPIDresponse(electrons, AliHFEV0cuts::kRecoElectron);
FillPIDresponse(pionsK0, AliHFEV0cuts::kRecoPionK0);
FillPIDresponse(pionsL, AliHFEV0cuts::kRecoPionL);
FillPIDresponse(protons, AliHFEV0cuts::kRecoProton);
CheckTenderV0pid(electrons, AliHFEV0cuts::kRecoElectron);
CheckTenderV0pid(pionsK0, AliHFEV0cuts::kRecoPionK0);
CheckTenderV0pid(pionsL, AliHFEV0cuts::kRecoPionL);
CheckTenderV0pid(protons, AliHFEV0cuts::kRecoProton);
fV0pid->Flush();
delete electrons;
delete pionsL;
delete pionsK0;
delete protons;
delete cleanElectrons;
}
void AliHFEpidQA::FillIllumination(const TObjArray * const tracks, Int_t species){
THnSparseF *hIllumination = dynamic_cast<THnSparseF *>(fOutput->Get("hIllumination"));
if(!hIllumination) return;
Double_t quantities[6]; memset(quantities, 0, sizeof(Double_t) *6);
TIter trackIter(tracks);
quantities[0] = species;
TObject *o = NULL; AliESDtrack *esdtrack = NULL;
while((o = trackIter())){
if(!TString(o->IsA()->GetName()).CompareTo("AliESDtrack")){
esdtrack = new AliESDtrack(*(static_cast<AliESDtrack *>(o)));
if(!esdtrack) continue;
} else if(!TString(o->IsA()->GetName()).CompareTo("AliAODrack")){
esdtrack = new AliESDtrack(static_cast<AliAODTrack *>(o));
if(!esdtrack) continue;
} else {
continue;
}
esdtrack->PropagateTo(300, fEvent->GetMagneticField());
quantities[1] = esdtrack->Pt();
quantities[2] = esdtrack->Eta();
quantities[3] = esdtrack->Phi();
quantities[4] = (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) ? 1 : 0;
quantities[5] = (esdtrack->GetStatus() & AliESDtrack::kTRDout) ? 1. : 0.;
hIllumination->Fill(quantities);
delete esdtrack;
}
}
void AliHFEpidQA::FillTPCinfo(AliESDtrack *const esdtrack, Int_t species){
THnSparseF *hTPCclusters = dynamic_cast<THnSparseF *>(fOutput->Get("hTPCclusters"));
if(!hTPCclusters) return;
Double_t quantities[6]; memset(quantities, 0, sizeof(Double_t) *6);
Double_t pidProbs[5];
const Int_t typePID[5] = {0, 2, 2, 3, 4};
quantities[0] = species;
esdtrack->GetTPCpid(pidProbs);
quantities[1] = esdtrack->P();
quantities[2] = esdtrack->GetTPCNcls();
quantities[3] = pidProbs[0];
quantities[4] = fESDpid->NumberOfSigmasTPC(esdtrack,(AliPID::EParticleType)typePID[species]);
quantities[5] = esdtrack->GetTPCsignal();
hTPCclusters->Fill(quantities);
}
void AliHFEpidQA::TestTRDResponse(const TObjArray * const tracks, Int_t species){
Int_t effInt[6] = {70, 75, 80, 85, 90, 95};
const char *sname[5] = {"Electron", "PionK0", "PionL", "Kaon", "Proton"};
AliESDtrack *track = NULL;
TIter trackIter(tracks);
Float_t momenta[6], p;
Int_t nmomenta;
Double_t probs[5], likeEl;
while((track = static_cast<AliESDtrack *>(trackIter()))){
if(track->GetTRDntrackletsPID() < 4) continue;
memset(momenta, 0, sizeof(Float_t) * 6); nmomenta = 0;
for(Int_t ily = 0; ily < 6; ily++){
if(track->GetTRDmomentum(ily) > 0.01) momenta[nmomenta++] = track->GetTRDmomentum(ily);
}
p = TMath::Mean(nmomenta, momenta);
track->GetTRDpid(probs);
likeEl = probs[0]/(probs[0] + probs[2]);
for(Int_t ieff = 0; ieff < 6; ieff++){
if(fESDpid->IdentifiedAsElectronTRD(track, static_cast<Double_t>(effInt[ieff])/100.)) fOutput->Fill(Form("hTRD_likeSel_%s_%dtls_%deff", sname[species], static_cast<Int_t>(track->GetTRDntrackletsPID()), effInt[ieff]), p, likeEl);
else fOutput->Fill(Form("hTRD_likeRej_%s_%dtls_%deff", sname[species], static_cast<Int_t>(track->GetTRDntrackletsPID()), effInt[ieff]), p, likeEl);
}
}
}
void AliHFEpidQA::MakePurity(const TObjArray *tracks, Int_t species){
if(!fMC) return;
AliDebug(3, Form("Doing Purity checks for species %d", species));
Int_t pdg = GetPDG(species);
TString hname;
switch(species){
case AliHFEV0cuts::kRecoElectron:
hname = "purityElectron";
break;
case AliHFEV0cuts::kRecoPionK0:
hname = "purityPionK0";
break;
case AliHFEV0cuts::kRecoPionL:
hname = "purityPionL";
break;
case AliHFEV0cuts::kRecoProton:
hname = "purityProton";
break;
default:
AliDebug(3, "Species not investigated");
return;
}
AliDebug(3, Form("Number of tracks: %d", tracks->GetEntries()));
TIter trackIter(tracks);
AliVParticle *recTrack = NULL, *mcTrack = NULL;
while((recTrack = dynamic_cast<AliVParticle *>(trackIter()))){
Int_t label = recTrack->GetLabel();
AliDebug(4, Form("MC Label %d", label));
mcTrack =fMC->GetTrack(TMath::Abs(label));
if(!mcTrack){
AliDebug(4, "MC track not available");
continue;
}
Int_t trackPdg = 0;
if(!TString(mcTrack->IsA()->GetName()).CompareTo("AliMCParticle")){
AliMCParticle *mcp = dynamic_cast<AliMCParticle *>(mcTrack);
if(!mcp) continue;
trackPdg = TMath::Abs(mcp->Particle()->GetPdgCode());
} else {
AliAODMCParticle *aodmcp = dynamic_cast<AliAODMCParticle *>(mcTrack);
if(!aodmcp) continue;
trackPdg = TMath::Abs(aodmcp->GetPdgCode());
}
if(trackPdg == pdg)
{
fOutput->Fill(hname, 0., recTrack->Pt());
}
else
fOutput->Fill(hname, 1., recTrack->Pt());
}
}
void AliHFEpidQA::FillElectronLikelihoods(const TObjArray * const particles, Int_t species){
Long_t status = 0;
const TString detname[4] = {"ITS", "TPC", "TRD", "TOF"};
TString specname;
switch(species){
case AliHFEV0cuts::kRecoElectron:
specname = "Electron";
break;
case AliHFEV0cuts::kRecoPionK0:
specname = "PionK0";
break;
case AliHFEV0cuts::kRecoPionL:
specname = "PionL";
break;
case AliHFEV0cuts::kRecoProton:
specname = "Proton";
break;
default:
AliDebug(2, Form("Species %d not investigated", species));
return;
};
AliVParticle *recTrack = NULL;
TIter trackIter(particles);
Double_t quantities[2];
Double_t pidProbs[5];
while((recTrack = dynamic_cast<AliVParticle *>(trackIter()))){
if(!TString(recTrack->IsA()->GetName()).CompareTo("AliESDtrack")){
AliESDtrack *esdTrack = dynamic_cast<AliESDtrack *>(recTrack);
if(!esdTrack) continue;
status = esdTrack->GetStatus();
Double_t pTPC = 0.;
pTPC = esdTrack->GetInnerParam() ? esdTrack->GetInnerParam()->P() : esdTrack->P();
Bool_t mcFound = kFALSE;
if(fMC){
Int_t label = esdTrack->GetLabel();
AliMCParticle *mcTrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(label));
Int_t pdg = GetPDG(species);
Int_t trackPdg = 0;
if(mcTrack){
trackPdg = TMath::Abs(mcTrack->Particle()->GetPdgCode());
}
if(pdg == trackPdg) mcFound = kTRUE;
}
quantities[0] = pTPC;
Bool_t detFlagSet = kFALSE;
for(Int_t idet = 0; idet < 4; idet++){
TString histname, histnameMC;
histname = "h" + detname[idet] + "_El_like_" + specname;
histnameMC = "h" + detname[idet] + "_El_like_MC_" + specname;
switch(idet){
case kITS: esdTrack->GetITSpid(pidProbs);
detFlagSet = status & AliESDtrack::kITSpid;
break;
case kTPC: esdTrack->GetTPCpid(pidProbs);
detFlagSet = status & AliESDtrack::kTPCpid;
break;
case kTRD: esdTrack->GetTRDpid(pidProbs);
detFlagSet = status & AliESDtrack::kTRDpid;
break;
case kTOF: esdTrack->GetTOFpid(pidProbs);
detFlagSet = status & AliESDtrack::kTOFpid;
break;
};
quantities[1] = pidProbs[AliPID::kElectron];
if(kTRD == idet && esdTrack->GetTRDntrackletsPID() != 6) continue;
if(detFlagSet){
fOutput->Fill(histname, quantities[0], quantities[1]);
if(mcFound)
fOutput->Fill(histnameMC, quantities[0], quantities[1]);
}
}
}
else{
}
}
}
void AliHFEpidQA::FillPIDresponse(const TObjArray * const particles, Int_t species){
TString hname, hname2, hname3;
const TString typeName[5] = {"Electron", "PionK0", "PionL", "Kaon", "Proton"};
const Int_t typePID[5] = {0, 2, 2, 3, 4};
Int_t run = fEvent->GetRunNumber();
AliVParticle *recTrack = NULL;
TIter trackIter(particles);
while((recTrack = dynamic_cast<AliVParticle *>(trackIter()))){
if(!TString(recTrack->IsA()->GetName()).CompareTo("AliESDtrack")){
AliESDtrack *esdTrack = dynamic_cast<AliESDtrack *>(recTrack);
if(!esdTrack) continue;
Float_t impactR = -1.;
Float_t impactZ = -1.;
esdTrack->GetImpactParameters(impactR, impactZ);
ULong_t status = 0;
status = esdTrack->GetStatus();
fOutput->Fill("hITS_kFlags", 5., species);
if(status & AliESDtrack::kITSin) fOutput->Fill("hITS_kFlags", 1., species);
if(status & AliESDtrack::kITSout) fOutput->Fill("hITS_kFlags", 2., species);
if(status & AliESDtrack::kITSrefit) fOutput->Fill("hITS_kFlags", 3., species);
if(status & AliESDtrack::kITSpid) fOutput->Fill("hITS_kFlags", 4., species);
fOutput->Fill("hTPC_kFlags", 5., species);
if(status & AliESDtrack::kTPCin) fOutput->Fill("hTPC_kFlags", 1., species);
if(status & AliESDtrack::kTPCout) fOutput->Fill("hTPC_kFlags", 2., species);
if(status & AliESDtrack::kTPCrefit) fOutput->Fill("hTPC_kFlags", 3., species);
if(status & AliESDtrack::kTPCpid) fOutput->Fill("hTPC_kFlags", 4., species);
fOutput->Fill("hTRD_kFlags", 5., species);
if(status & AliESDtrack::kTRDin) fOutput->Fill("hTRD_kFlags", 1., species);
if(status & AliESDtrack::kTRDout) fOutput->Fill("hTRD_kFlags", 2., species);
if(status & AliESDtrack::kTRDrefit) fOutput->Fill("hTRD_kFlags", 3., species);
if(status & AliESDtrack::kTRDpid) fOutput->Fill("hTRD_kFlags", 4., species);
fOutput->Fill("hTOF_kFlags", 5., species);
if(status & AliESDtrack::kTOFin) fOutput->Fill("hTOF_kFlags", 1., species);
if(status & AliESDtrack::kTOFout) fOutput->Fill("hTOF_kFlags", 2., species);
if(status & AliESDtrack::kTOFrefit) fOutput->Fill("hTOF_kFlags", 3., species);
if(status & AliESDtrack::kTOFpid) fOutput->Fill("hTOF_kFlags", 4., species);
if(status & AliESDtrack::kITSpid){
Double_t p = esdTrack->P();
Double_t dEdxSamples[4];
esdTrack->GetITSdEdxSamples(dEdxSamples);
Int_t nSamples = 0;
Double_t dEdxSum = 0.;
hname = "hITS_dEdx_" + typeName[species];
for(Int_t i=0; i<4; ++i){
if(dEdxSamples[i] > 0){
nSamples++;
fOutput->Fill(hname, i+1, p, dEdxSamples[i]);
dEdxSum += dEdxSamples[i];
}
}
if(4 == nSamples)fOutput->Fill(hname, 0, p, dEdxSum);
fOutput->Fill("hITS_dEdx_nSamples", nSamples);
Double_t signal = esdTrack->GetITSsignal();
hname = "hITS_Signal_" + typeName[species];
fOutput->Fill(hname, p, signal);
Double_t nsigma = fESDpid->NumberOfSigmasITS(esdTrack,(AliPID::EParticleType)typePID[species]);
hname = "hITS_nSigma_" + typeName[species];
fOutput->Fill(hname, p, nsigma);
Double_t itsPID[5] = {-1, -1, -1, -1, -1};
esdTrack->GetITSpid(itsPID);
Int_t ix = GetMaxPID(itsPID);
hname = "hITS_PID_p_" + typeName[species];
fOutput->Fill(hname, p, ix);
}
if(status & AliESDtrack::kTPCpid){
FillTPCinfo(esdTrack, species);
Double_t p = esdTrack->GetInnerParam() ? esdTrack->GetInnerParam()->P() : esdTrack->P();
Double_t dEdx = esdTrack->GetTPCsignal();
hname = "hTPC_dEdx_" + typeName[species];
fOutput->Fill(hname, p, dEdx);
Double_t nsigma = fESDpid->NumberOfSigmasTPC(esdTrack,(AliPID::EParticleType)typePID[species]);
hname = "hTPC_nSigma_" + typeName[species];
fOutput->Fill(hname, p, nsigma);
hname = "hTPC_PID_p_" + typeName[species];
Double_t tpcPID[5] = {-1, -1, -1, -1, -1};
esdTrack->GetTPCpid(tpcPID);
Int_t ix = GetMaxPID(tpcPID);
fOutput->Fill(hname, p, ix);
}
if(status & AliESDtrack::kTRDpid){
Double_t p = esdTrack->GetOuterParam() ? esdTrack->GetOuterParam()->P() : esdTrack->P();
Int_t ntrk = esdTrack->GetTRDntrackletsPID();
hname = "hTRD_trk_" + typeName[species];
fOutput->Fill(hname, p, ntrk);
hname = "hTRD_PID_p_" + typeName[species];
Double_t trdPID[5] = {-1., -1., -1., -1., -1};
esdTrack->GetTRDpid(trdPID);
Int_t ix = GetMaxPID(trdPID);
fOutput->Fill(hname, p, ix);
Int_t ncls = esdTrack->GetTRDncls();
hname = "hTRD_cls_" + typeName[species];
fOutput->Fill(hname, p, ncls);
hname = "hTRD_Nslc_" + typeName[species];
hname2 = "hTRD_slc_" + typeName[species];
hname3 = "hTRD_dEdx_" + typeName[species];
Int_t nSlices = esdTrack->GetNumberOfTRDslices();
Double_t sumTot = 0.;
Int_t not0Tot = 0;
for(Int_t l=0; l< 6; ++l){
Double_t trkData[7] = {-1.,-1, -1, -1, -1, -1, -1};
trkData[0] = species;
trkData[1] = p;
trkData[2] = l;
trkData[3] = ntrk;
trkData[5] = ncls;
Double_t sum = 0.;
Int_t not0 = 0;
for(Int_t s=0; s<nSlices; ++s){
Double_t slice = esdTrack->GetTRDslice(l, s);
sum += slice;
if(slice > 0){
not0 += 1;
fOutput->Fill(hname2, p, s);
}
}
trkData[4] = not0;
fOutput->Fill(hname, p, not0);
fOutput->Fill(hname3, p, sum);
if(sum > 0){
sumTot += sum;
not0Tot += 1;
}
if(not0Tot > 0 && fNNref){
Double_t likelihoods[5] = {-1., -1., -1., -1., -1};
TRDlikeTracklet(l, esdTrack, likelihoods);
trkData[6] = likelihoods[0];
}
if(not0Tot) fOutput->Fill("hTRDtracklets", trkData);
}
}
if(status & AliESDtrack::kTOFpid){
Double_t p = esdTrack->GetOuterParam() ? esdTrack->GetOuterParam()->P() : esdTrack->P();
Double_t t0 = fESDpid->GetTOFResponse().GetStartTime(esdTrack->P());
hname = "hTOF_beta_" + typeName[species];
Float_t beta = TOFbeta(esdTrack);
fOutput->Fill(hname, p, beta);
fOutput->Fill("hTOF_beta_all", p, beta);
Double_t nsigma = fESDpid->NumberOfSigmasTOF(esdTrack,(AliPID::EParticleType)typePID[species]);
hname = "hTOF_nSigma_" + typeName[species];
fOutput->Fill(hname, p, nsigma);
hname = "hTOF_PID_p_" + typeName[species];
Double_t tofPID[5] = {-1., -1., -1., -1., -1};
esdTrack->GetTOFpid(tofPID);
Int_t ix = GetMaxPID(tofPID);
fOutput->Fill(hname, p, ix);
Double_t times[AliPID::kSPECIESC];
esdTrack->GetIntegratedTimes(times);
Double_t tItrackL = esdTrack->GetIntegratedLength();
Double_t tTOFsignal = esdTrack->GetTOFsignal();
Double_t dT = tTOFsignal - t0 - times[2];
fOutput->Fill("hTOF_qa_TmT0mT", run*1.0, dT);
fOutput->Fill("hTOF_qa_length", run*1.0, tItrackL);
}
}
else{
continue;
}
}
}
void AliHFEpidQA:: CheckEvent(){
Double_t t0 = fESDpid->GetTOFResponse().GetTimeZero();
Int_t run = fEvent->GetRunNumber();
Double_t data[2] = {run*1.0, t0*1000.};
fOutput->Fill("hEvent_T0", data);
}
TList *AliHFEpidQA::GetOutput(){
return fOutput->GetList();
}
TList *AliHFEpidQA::GetV0pidQA(){
return fV0pid->GetListOfQAhistograms();
}
TList *AliHFEpidQA::GetV0pidMC(){
if(fV0pidMC)
return fV0pidMC->GetListOfQAhistograms();
return NULL;
}
void AliHFEpidQA::RecalculateTRDpid(AliESDtrack * , Double_t * ) const{
}
void AliHFEpidQA::RecalculateTRDpid(AliAODTrack * , Double_t * ) const{
}
Float_t AliHFEpidQA::TOFbeta(const AliESDtrack * const track) const {
Double_t l = track->GetIntegratedLength();
Double_t t = track->GetTOFsignal();
Double_t t0 = fESDpid->GetTOFResponse().GetTimeZero();
if(l < 360. || l > 800.) return 0.;
if(t <= 0.) return 0.;
if(t0 >999990.0) return 0.;
t -= t0;
l *= 0.01;
t *= 1e-12;
Double_t v = l / t;
Float_t beta = v / TMath::C();
return beta;
}
Int_t AliHFEpidQA::GetMaxPID(const Double_t *pidProbs) const {
Int_t ix = -1;
Double_t tmp = 0.2;
for(Int_t i=0; i<5; ++i){
if(pidProbs[i] > tmp){
ix = i;
tmp = pidProbs[i];
}
}
return ix;
}
Int_t AliHFEpidQA::GetPDG(Int_t species){
Int_t pdg = 0;
switch(species){
case AliHFEV0cuts::kRecoElectron:
pdg = TMath::Abs(kElectron);
break;
case AliHFEV0cuts::kRecoPionK0:
pdg = TMath::Abs(kPiPlus);
break;
case AliHFEV0cuts::kRecoPionL:
pdg = TMath::Abs(kPiPlus);
break;
case AliHFEV0cuts::kRecoProton:
pdg = TMath::Abs(kProton);
break;
default:
AliDebug(3, "Species not recognised");
return 0;
}
return pdg;
}
TObjArray * AliHFEpidQA::MakeTrackList(const TObjArray *tracks) const {
TObjArray *output = new TObjArray;
TIter trackInfos(tracks);
AliHFEV0info *trackInfo = NULL;
while((trackInfo = dynamic_cast<AliHFEV0info *>(trackInfos())))
output->Add(trackInfo->GetTrack());
return output;
}
TObjArray * AliHFEpidQA::MakeCleanListElectrons(const TObjArray *electrons) const {
TObjArray *tracks = new TObjArray;
TIter candidates(electrons);
AliESDEvent *esd;
AliHFEV0info *hfetrack;
const Double_t kSigmaTight = 2.;
const Double_t kSigmaLoose = 2.;
const Double_t shift = 0.57;
if((esd = dynamic_cast<AliESDEvent *>(fEvent))){
AliESDtrack *track = NULL, *partnerTrack = NULL;
while((hfetrack = dynamic_cast<AliHFEV0info *>(candidates()))){
track = dynamic_cast<AliESDtrack *>(hfetrack->GetTrack());
if(!track) continue;
partnerTrack = esd->GetTrack(hfetrack->GetPartnerID());
Double_t nSigmaTrack = TMath::Abs(fESDpid->NumberOfSigmasTPC(track, AliPID::kElectron) - shift);
Double_t nSigmaPartner = TMath::Abs(fESDpid->NumberOfSigmasTPC(partnerTrack, AliPID::kElectron) - shift);
if((nSigmaTrack < kSigmaTight && nSigmaPartner < kSigmaLoose) || (nSigmaTrack < kSigmaLoose && nSigmaPartner < kSigmaTight))
tracks->Add(track);
}
}
return tracks;
}
void AliHFEpidQA::CheckTenderV0pid(const TObjArray * const particles, Int_t species){
const Int_t id[5] = {0, 1, 1, -1, 2};
AliVParticle *recTrack = NULL;
TIter trackIter(particles);
while((recTrack = dynamic_cast<AliVParticle *>(trackIter()))){
if(!TString(recTrack->IsA()->GetName()).CompareTo("AliESDtrack")){
AliESDtrack *track = dynamic_cast<AliESDtrack *>(recTrack);
if(!track) continue;
Int_t tPID = GetTenderV0pid(track);
fOutput->Fill("h_tender_check_01", id[species]*1.0, tPID);
fOutput->Fill("h_tender_check_02", id[species], id[species]*1.0, tPID);
}
}
}
Int_t AliHFEpidQA::GetTenderV0pid(AliESDtrack * const track){
Int_t pid = -1;
if(!track){
return pid;
}
Int_t nTimes = 0;
if(track->TestBit(2<<14)){
pid = 0;
nTimes++;
}
if(track->TestBit(2<<15)){
pid = 1;
nTimes++;
}
if(track->TestBit(2<<16)){
pid = 2;
nTimes++;
}
if(nTimes > 1){
AliWarning("V0 track labeled multiple times by the V0 tender");
pid = -1;
}
return pid;
}
Double_t AliHFEpidQA::TRDlikeTracklet(Int_t layer, AliESDtrack * const track, Double_t *likelihood){
const Double_t cScaleGain = 1./ 16000.;
const Float_t pBins[11] ={0.6, 0.8, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0};
if(!track) return kFALSE;
Float_t p = track->GetTRDmomentum(layer);
if(p < 0) return kFALSE;
Int_t mombin = TRDmomBin(p);
if(mombin < 0) return -1.;
Float_t dEdxTRDsum = 0;
Float_t dEdxTRD[8];
Double_t ddEdxTRD[8];
Double_t prob[AliPID::kSPECIES];
for(Int_t is = 0; is < AliPID::kSPECIES; is++){
likelihood[is] = 0.2;
prob[is] = 0.2;
}
Double_t sum = 0.;
for(Int_t islice = 0; islice<8; islice++){
dEdxTRD[islice]=0.;
ddEdxTRD[islice]=0.;
dEdxTRD[islice]=track->GetTRDslice(layer,islice);
dEdxTRDsum += dEdxTRD[islice];
ddEdxTRD[islice]=(Double_t)dEdxTRD[islice]*cScaleGain;
}
for(Int_t is = 0; is < AliPID::kSPECIES; is++){
Double_t probloc1, probloc2;
if(mombin == 0 && mombin < pBins[0]){
prob[is] = fNet[mombin]->Evaluate(is, ddEdxTRD);
}
else if(mombin == 10 && mombin >= pBins[10]){
prob[is] = fNet[mombin]->Evaluate(is, ddEdxTRD);
}
else{
Int_t mombin1 = 0, mombin2 = 0;
if(p < pBins[mombin]) {mombin1 = mombin -1; mombin2 = mombin;}
if(p >= pBins[mombin]) {mombin1 = mombin; mombin2 = mombin+1;}
probloc1 = fNet[mombin1]->Evaluate(is, ddEdxTRD);
probloc2 = fNet[mombin2]->Evaluate(is, ddEdxTRD);
prob[is] = probloc1 + (probloc2-probloc1)*(p-pBins[mombin1])/(pBins[mombin2]-pBins[mombin1]);
}
likelihood[is] = prob[is];
sum += likelihood[is];
}
return sum;
}
Int_t AliHFEpidQA::TRDmomBin(Double_t p) const {
const Float_t pBins[11] ={0.6, 0.8, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0};
Int_t pBin1 = -1;
Int_t pBin2 = -1;
if(p < 0) return -1;
if(p < pBins[0]) return 0;
if(p >= pBins[10]) return 10;
for(Int_t iMomBin = 1; iMomBin < 11; iMomBin++){
if(p < pBins[iMomBin]){
pBin1 = iMomBin - 1;
pBin2 = iMomBin;
}
else
continue;
if(p - pBins[pBin1] >= pBins[pBin2] - p){
return pBin2;
}
else{
return pBin1;
}
}
return -1;
}
Int_t AliHFEpidQA::GetCentrality(AliVEvent* const fInputEvent){
Int_t bin = -1;
if(fIsPbPb)
{
AliCentrality *centrality = fInputEvent->GetCentrality();
Float_t fCentralityFtemp = -1;
fCentralityFtemp = centrality->GetCentralityPercentile("V0M");
Float_t centralityLimits[12] = {0.,5.,10., 20., 30., 40., 50., 60.,70.,80., 90., 100.};
for(Int_t ibin = 0; ibin < 11; ibin++){
if(fCentralityFtemp >= centralityLimits[ibin] && fCentralityFtemp < centralityLimits[ibin+1]){
bin = ibin;
break;
}
}
}
AliDebug(2, Form("Centrality class %d\n", bin));
return bin;
}
Int_t AliHFEpidQA::GetMultiplicityITS(AliVEvent* const fInputEvent) const
{
Int_t nTracklets = 0;
Int_t nAcc = 0;
Double_t etaRange = 1.6;
Int_t bin = -1;
nTracklets = ((AliESDEvent*) fInputEvent)->GetMultiplicity()->GetNumberOfTracklets();
for (Int_t nn = 0; nn < nTracklets; nn++) {
Double_t eta = ((AliESDEvent*)fInputEvent)->GetMultiplicity()->GetEta(nn);
if (TMath::Abs(eta) < etaRange) nAcc++;
}
Int_t itsMultiplicity = nAcc;
Int_t multiplicityLimits[8] = {0, 1, 9, 17, 25, 36, 60, 500};
for(Int_t ibin = 0; ibin < 7; ibin++){
if(itsMultiplicity >= multiplicityLimits[ibin] && itsMultiplicity < multiplicityLimits[ibin + 1]){
bin = ibin;
break;
}
}
return bin;
}