#include <TClass.h>
#include <TH2.h>
#include <THnSparse.h>
#include <TString.h>
#include <TMath.h>
#include "AliESDInputHandler.h"
#include "AliMagF.h"
#include "AliPIDResponse.h"
#include "AliLog.h"
#include "AliPID.h"
#include "AliVParticle.h"
#include "AliHFEcollection.h"
#include "AliHFEpid.h"
#include "AliHFEpidBase.h"
#include "AliHFEpidQAmanager.h"
#include "AliHFEpidTPC.h"
#include "AliHFEpidEMCAL.h"
#include "AliHFEtools.h"
#include "AliHFEemcalPIDqa.h"
#include "AliTrackerBase.h"
ClassImp(AliHFEemcalPIDqa)
AliHFEemcalPIDqa::AliHFEemcalPIDqa():
AliHFEdetPIDqa()
, fHistos(NULL)
{
}
AliHFEemcalPIDqa::AliHFEemcalPIDqa(const char* name):
AliHFEdetPIDqa(name, "QA for EMCAL")
, fHistos(NULL)
{
}
AliHFEemcalPIDqa::AliHFEemcalPIDqa(const AliHFEemcalPIDqa &o):
AliHFEdetPIDqa(o)
, fHistos()
{
o.Copy(*this);
}
AliHFEemcalPIDqa &AliHFEemcalPIDqa::operator=(const AliHFEemcalPIDqa &o){
AliHFEdetPIDqa::operator=(o);
if(&o != this) o.Copy(*this);
return *this;
}
AliHFEemcalPIDqa::~AliHFEemcalPIDqa(){
if(fHistos) delete fHistos;
}
void AliHFEemcalPIDqa::Copy(TObject &o) const {
AliHFEemcalPIDqa &target = dynamic_cast<AliHFEemcalPIDqa &>(o);
if(target.fHistos){
delete target.fHistos;
target.fHistos = NULL;
}
if(fHistos) target.fHistos = new AliHFEcollection(*fHistos);
}
Long64_t AliHFEemcalPIDqa::Merge(TCollection *coll){
if(!coll) return 0;
if(coll->IsEmpty()) return 1;
TIter it(coll);
AliHFEemcalPIDqa *refQA = NULL;
TObject *o = NULL;
Long64_t count = 0;
TList listHistos;
while((o = it())){
refQA = dynamic_cast<AliHFEemcalPIDqa *>(o);
if(!refQA) continue;
listHistos.Add(refQA->fHistos);
count++;
}
fHistos->Merge(&listHistos);
return count + 1;
}
void AliHFEemcalPIDqa::Initialize(){
fHistos = new AliHFEcollection("emcalqahistos", "Collection of EMCAL QA histograms");
const Int_t kCentralityBins = 11;
const Double_t kMinP = 0.;
const Double_t kMaxP = 50.;
const Double_t kTPCSigMim = 40.;
const Double_t kTPCSigMax = 140.;
Int_t nBins[16] = {AliPID::kSPECIES + 1, 500, 500, 400, 300, 100, 600, 300, 100, 40, 300, 300, 300, kCentralityBins, 6, 2};
Double_t min[16] = {-1, kMinP, kMinP, kTPCSigMim, 0., -1.0, -8.0, 0, 0, 0, 0.0, 0.0, 0.0, 0, -3.0, 0.};
Double_t max[16] = {AliPID::kSPECIES, kMaxP, kMaxP, kTPCSigMax, 6.0, 1.0, 4.0, 3.0, 0.1, 40, 3.0, 3.0, 3.0, 11., 3.0, 2.};
fHistos->CreateTHnSparse("EMCAL_TPCdedx", "EMCAL signal; species; p [GeV/c]; pT [GeV/c] ; TPC signal [a.u.]; phi ; eta ; nSig ; E/p ; Rmatch ; Ncell ; M02 ; M20 ; Disp ; Centrality;charge ;PID Step; ", 16, nBins, min, max);
}
void AliHFEemcalPIDqa::ProcessTrack(const AliHFEpidObject *track,AliHFEdetPIDqa::EStep_t step){
Float_t centrality = track->GetCentrality();
const AliESDtrack *esdtrack = static_cast<const AliESDtrack *>(track->GetRecTrack());
Int_t species = track->GetAbInitioPID();
if(species >= AliPID::kSPECIES) species = -1;
AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fQAmanager->GetDetectorPID(AliHFEpid::kTPCpid));
const AliPIDResponse *pidResponse = tpcpid ? tpcpid->GetPIDResponse() : NULL;
Float_t nSigmatpc = pidResponse ? pidResponse->NumberOfSigmasTPC(static_cast<const AliVTrack *>(track->GetRecTrack()), AliPID::kElectron) : 0.;
AliDebug(2, Form("nSigmatpc = %f\n",nSigmatpc));
double charge = esdtrack->Charge();
Double_t showerEMC[4];
TVector3 emcsignal = MomentumEnergyMatchV2(esdtrack,showerEMC);
AliDebug(2, Form("shower shape ; %f %f %f %f \n",showerEMC[0],showerEMC[1],showerEMC[2],showerEMC[3]));
Double_t contentSignal[16];
contentSignal[0] = species;
contentSignal[1] = track->GetRecTrack()->P();
contentSignal[2] = track->GetRecTrack()->Pt();
contentSignal[3] = esdtrack->GetTPCsignal();
double phi = track->GetRecTrack()->Phi();
double eta = track->GetRecTrack()->Eta();
contentSignal[4] = phi;
contentSignal[5] = eta;
contentSignal[6] = nSigmatpc;
contentSignal[7] = emcsignal(0);
contentSignal[8] = emcsignal(2);
contentSignal[9] = showerEMC[0];
contentSignal[10] = showerEMC[1];
contentSignal[11] = showerEMC[2];
contentSignal[12] = showerEMC[3];
contentSignal[13] = centrality;
contentSignal[14] = charge;
contentSignal[15] = step == AliHFEdetPIDqa::kBeforePID ? 0. : 1.;
fHistos->Fill("EMCAL_TPCdedx", contentSignal);
}
TH1 *AliHFEemcalPIDqa::GetHistogram(const char *name){
if(!fHistos) return NULL;
TObject *histo = fHistos->Get(name);
if(!histo->InheritsFrom("TH1")) return NULL;
return dynamic_cast<TH1 *>(histo);
}
TVector3 AliHFEemcalPIDqa::MomentumEnergyMatchV2(const AliESDtrack *esdtrack, Double_t *shower) const
{
TVector3 refVec(-9999,-9999,-9999);
Double_t matchclsE = 9999.9;
for(int i=0; i<4; i++)shower[i] = -9999;
const AliESDEvent *evt = esdtrack->GetESDEvent();
Int_t icl = esdtrack->GetEMCALcluster();
AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
if(cluster==NULL) return refVec;
if(!cluster->IsEMCAL()) return refVec;
double delphi = cluster->GetTrackDx();
double deleta = cluster->GetTrackDz();
double rmatch1 = sqrt(pow(delphi,2)+pow(deleta,2));
matchclsE = cluster->E();
double feop = matchclsE/esdtrack->P();
shower[0] = cluster->GetNCells();
shower[1] = cluster->GetM02();
shower[2] = cluster->GetM20();
shower[3] = cluster->GetDispersion();
TVector3 emcsignal(feop,0.0,rmatch1);
return emcsignal;
}