[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEtpcPIDqa.cxx

/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/
//
// Class AliHFEtpcPIDqa
// Monitoring TPC PID in the HFE PID montioring framework. The following
// quantities are monitored:
//   TPC dE/dx (Number of sigmas)
//   TPC dE/dx (Absolute values)
// (Always as function of momentum, particle species and centrality 
// before and after cut)
// More information about the PID monitoring framework can be found in
// AliHFEpidQAmanager.cxx and AliHFEdetPIDqa.cxx
//
// Author:
//    Markus Fasel <M.Fasel@gsi.de>
//
#include <TBrowser.h>
#include <TClass.h>
#include <TH2.h>
#include <THnSparse.h>
#include <TString.h>

#include "AliAODEvent.h"
#include "AliAODTrack.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 "AliHFEtools.h"
#include "AliHFEtpcPIDqa.h"

ClassImp(AliHFEtpcPIDqa)

//_________________________________________________________
AliHFEtpcPIDqa::AliHFEtpcPIDqa():
    AliHFEdetPIDqa()
  , fHistos(NULL)
  , fBrowseCentrality(-1)
{
  //
  // Dummy constructor
  //
}

//_________________________________________________________
AliHFEtpcPIDqa::AliHFEtpcPIDqa(const char* name):
    AliHFEdetPIDqa(name, "QA for TPC")
  , fHistos(NULL)
  , fBrowseCentrality(-1)
{
  //
  // Default constructor
  //
}

//_________________________________________________________
AliHFEtpcPIDqa::AliHFEtpcPIDqa(const AliHFEtpcPIDqa &o):
    AliHFEdetPIDqa(o)
  , fHistos(NULL)
  , fBrowseCentrality(o.fBrowseCentrality)
{
  //
  // Copy constructor
  //
  o.Copy(*this);
}

//_________________________________________________________
AliHFEtpcPIDqa &AliHFEtpcPIDqa::operator=(const AliHFEtpcPIDqa &o){
  //
  // Do assignment
  //
  AliHFEdetPIDqa::operator=(o);
  if(&o != this) o.Copy(*this);
  return *this;
}

//_________________________________________________________
AliHFEtpcPIDqa::~AliHFEtpcPIDqa(){
  //
  // Destructor
  //
  if(fHistos) delete fHistos;
}

//_________________________________________________________
void AliHFEtpcPIDqa::Copy(TObject &o) const {
  //
  // Make copy
  //
  AliHFEtpcPIDqa &target = dynamic_cast<AliHFEtpcPIDqa &>(o);
  if(target.fHistos){
    delete target.fHistos;
    target.fHistos = NULL;
  }
  if(fHistos) target.fHistos = new AliHFEcollection(*fHistos);
  target.fBrowseCentrality = fBrowseCentrality;
}

//_________________________________________________________
Long64_t AliHFEtpcPIDqa::Merge(TCollection *coll){
  //
  // Merge with other objects
  //
  if(!coll) return 0;
  if(coll->IsEmpty()) return 1;

  TIter it(coll);
  AliHFEtpcPIDqa *refQA = NULL;
  TObject *o = NULL;
  Long64_t count = 0;
  TList listHistos;
  while((o = it())){
    refQA = dynamic_cast<AliHFEtpcPIDqa *>(o);
    if(!refQA) continue;

    listHistos.Add(refQA->fHistos);
    count++; 
  }
  fHistos->Merge(&listHistos);
  return count + 1;
}

//_________________________________________________________
void AliHFEtpcPIDqa::Browse(TBrowser *b){
  //
  // Browse the PID QA
  //
  if(b){
    if(fHistos){
      b->Add(fHistos, fHistos->GetName());

      // Make Projections of the dE/dx Spectra and add them to a new Folder
      TString specnames[AliPID::kSPECIES+1] = {"All", "Electrons", "Muon", "Pions", "Kaon", "Protons"};
      Int_t specind[AliPID::kSPECIES+1] = {-1, AliPID::kElectron, AliPID::kMuon, AliPID::kPion, AliPID::kKaon, AliPID::kProton};
      TList *listdEdx = new TList;
      listdEdx->SetOwner();
      TList *listNsigma = new TList;
      listNsigma->SetOwner();

      TH2 *hptr = NULL; 
      for(Int_t ispec = 0; ispec < AliPID::kSPECIES+1; ispec++){
        for(Int_t istep = 0; istep < 2; istep++){
          hptr = MakeSpectrumdEdx(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], fBrowseCentrality);
          hptr->SetName(Form("hTPCdEdx%s%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After" , fBrowseCentrality == -1 ? "MinBias" : Form("Cent%d", fBrowseCentrality)));
          listdEdx->Add(hptr);
          hptr = MakeSpectrumNSigma(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], fBrowseCentrality);
          hptr->SetName(Form("hTPCnsigma%s%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After", fBrowseCentrality == -1 ? "MinBias" : Form("Cent%d", fBrowseCentrality)));
          listNsigma->Add(hptr);
        }
      }
      
      b->Add(listdEdx, "Projections dE/dx");
      b->Add(listNsigma, "Projections NSigma");
    }
  }
}

//_________________________________________________________
void AliHFEtpcPIDqa::Initialize(){
  //
  // Define Histograms
  //

  fHistos = new AliHFEcollection("tpcqahistos", "Collection of TPC QA histograms");

  // Make common binning
  const Int_t kNdim = 7;
  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.;
  const Double_t kMinEta = -0.9;
  const Double_t kMaxEta = 0.9;

  // Quantities where one can switch between low and high resolution
  Int_t kPbins = (fQAmanager->HasHighResolutionHistos() || fQAmanager->HasMidResolutionHistos()) ? 1000 : 100;
  Int_t kDedxbins = fQAmanager->HasHighResolutionHistos() ? 400 : 200;
  Int_t kSigmaBins = fQAmanager->HasHighResolutionHistos() ? 1400 : 240;
  kSigmaBins = fQAmanager->HasMidResolutionHistos() ? 400 : kSigmaBins;
  Int_t kEtabins = fQAmanager->HasHighResolutionEtaHistos() ? 33 : 18;

  // 1st histogram: TPC dEdx: (species, p, dEdx, step)
  Int_t nBinsdEdx[kNdim] = {kPIDbins, kPbins, kDedxbins, kSteps, kCentralityBins, kEtabins, 33};
  Double_t mindEdx[kNdim] =  {kMinPID, kMinP, 20., 0., 0., kMinEta, 0};
  Double_t maxdEdx[kNdim] =  {kMaxPID, kMaxP, 130, 2., 11., kMaxEta, 2000};
  fHistos->CreateTHnSparse("tpcDedx", "TPC signal; species; p [GeV/c]; TPC signal [a.u.]; Selection Step; Centrality; Eta; pVx contrib", kNdim, nBinsdEdx, mindEdx, maxdEdx);
  
  // 2nd histogram: TPC sigmas: (species, p nsigma, step)
  Int_t nBinsSigma[kNdim] = {kPIDbins, kPbins, kSigmaBins, kSteps, kCentralityBins, kEtabins, 33};
  Double_t minSigma[kNdim] = {kMinPID, kMinP, -12., 0., 0., kMinEta, 0};
  Double_t maxSigma[kNdim] = {kMaxPID, kMaxP, 12., 2., 11., kMaxEta, 2000};
  fHistos->CreateTHnSparse("tpcnSigma", "TPC signal; species; p [GeV/c]; TPC signal [a.u.]; Selection Step; Centrality; Eta; pVx contrib", kNdim, nBinsSigma, minSigma, maxSigma);

  // General TPC QA
}

//_________________________________________________________
void AliHFEtpcPIDqa::ProcessTrack(const AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step){
  //
  // Fill TPC histograms
  //
  AliDebug(1, Form("QA started for TPC PID for step %d", (Int_t)step));
  AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
  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;
  Double_t contentSignal[7];
  contentSignal[0] = species;
  contentSignal[1] = tpcpid ? tpcpid->GetP(track->GetRecTrack(), anatype) : 0.;
  contentSignal[2] = GetTPCsignal(track->GetRecTrack(), anatype);
  contentSignal[3] = step;
  contentSignal[4] = track->GetCentrality();
  contentSignal[5] = GetEta(track->GetRecTrack(), anatype);
  contentSignal[6] = fQAmanager->HasFillMultiplicity() ? track->GetMultiplicity() : 0;
  fHistos->Fill("tpcDedx", contentSignal);

  if(track->HasCorrectedTPCnSigma())
    contentSignal[2] = track->GetCorrectedTPCnSigma();
  else
    contentSignal[2] = pidResponse ? pidResponse->NumberOfSigmasTPC(track->GetRecTrack(), AliPID::kElectron) : -100.;
  
  fHistos->Fill("tpcnSigma", contentSignal);
}

//_________________________________________________________
Double_t AliHFEtpcPIDqa::GetTPCsignal(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anatype){
  //
  // Get TPC signal for ESD and AOD track
  //
  Double_t tpcSignal = 0.;
  if(anatype == AliHFEpidObject::kESDanalysis){
    const AliESDtrack *esdtrack = static_cast<const AliESDtrack *>(track);
    tpcSignal = esdtrack->GetTPCsignal();
  } else {
    const AliAODTrack *aodtrack = static_cast<const AliAODTrack *>(track);
    tpcSignal = aodtrack->GetDetPid() ? aodtrack->GetDetPid()->GetTPCsignal() : 0.;
  }
  return tpcSignal;
}


//_________________________________________________________
Double_t AliHFEtpcPIDqa::GetEta(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anatype){
  //
  // Get TPC signal for ESD and AOD track
  //
  Double_t eta = 0.;
  if(anatype == AliHFEpidObject::kESDanalysis){
    const AliESDtrack *esdtrack = static_cast<const AliESDtrack *>(track);
    eta = esdtrack->Eta();
  } else {
    const AliAODTrack *aodtrack = static_cast<const AliAODTrack *>(track);
    eta = aodtrack->Eta(); 
  }
  return eta;
}

//_________________________________________________________
TH2 *AliHFEtpcPIDqa::MakeSpectrumdEdx(AliHFEdetPIDqa::EStep_t istep, Int_t species, Int_t centralityClass){
  //
  // Plot the Spectrum
  //
  THnSparseF *hSignal = dynamic_cast<THnSparseF *>(fHistos->Get("tpcDedx"));
  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 hname = Form("hTPCsignal%s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after"), 
          htitle = Form("TPC dE/dx Spectrum %s selection", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
  if(species > -1){
    hname += AliPID::ParticleName(species);
    htitle += Form(" for %ss", AliPID::ParticleName(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 {
      // Project centrality classes
      // -1 is Min. Bias
      hSignal->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
      centname = Form("Cent%d", centralityClass);
      centtitle = Form(" Centrality %d", centralityClass);
    }
  } else {
    centname = centtitle = "MinBias";
    hasCentralityInfo= kFALSE;
  }
  hname += centtitle;
  htitle += centtitle;

  TH2 *hTmp = hSignal->Projection(2,1);
  hTmp->SetName(hname.Data());
  hTmp->SetTitle(htitle.Data());
  hTmp->SetStats(kFALSE);
  hTmp->GetXaxis()->SetTitle("p [GeV/c]");
  hTmp->GetYaxis()->SetTitle("TPC signal [a.u.]");
  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;
}

//_________________________________________________________
TH2 *AliHFEtpcPIDqa::MakeSpectrumNSigma(AliHFEdetPIDqa::EStep_t istep, Int_t species, Int_t centralityClass){
  //
  // Plot the Spectrum
  //
  THnSparseF *hSignal = dynamic_cast<THnSparseF *>(fHistos->Get("tpcnSigma"));
  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 hname = Form("hTPCsigma%s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after"), 
          htitle = Form("TPC dE/dx Spectrum[#sigma] %s selection", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
  if(species > -1){
    hname += AliPID::ParticleName(species);
    htitle += Form(" for %ss", AliPID::ParticleName(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 {
      // Project centrality classes
      // -1 is Min. Bias
      hSignal->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
      centname = Form("Cent%d", centralityClass);
      centtitle = Form(" Centrality %d", centralityClass);
    }
  } else {
    centname = centtitle = "MinBias";
    hasCentralityInfo= kFALSE;
  }
  hname += centtitle;
  htitle += centtitle;

  TH2 *hTmp = hSignal->Projection(2,1);
  hTmp->SetName(hname.Data());
  hTmp->SetTitle(htitle.Data());
  hTmp->SetStats(kFALSE);
  hTmp->GetXaxis()->SetTitle("p [GeV/c]");
  hTmp->GetYaxis()->SetTitle("TPC dE/dx - <dE/dx>|_{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;
}
Commit	Line	Data
3a72645a	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15	//
	16	// Class AliHFEtpcPIDqa
	17	// Monitoring TPC PID in the HFE PID montioring framework. The following
	18	// quantities are monitored:
	19	// TPC dE/dx (Number of sigmas)
	20	// TPC dE/dx (Absolute values)
	21	// (Always as function of momentum, particle species and centrality
	22	// before and after cut)
	23	// More information about the PID monitoring framework can be found in
	24	// AliHFEpidQAmanager.cxx and AliHFEdetPIDqa.cxx
	25	//
	26	// Author:
	27	// Markus Fasel <M.Fasel@gsi.de>
	28	//
bf892a6a	29	#include <TBrowser.h>
3a72645a	30	#include <TClass.h>
	31	#include <TH2.h>
	32	#include <THnSparse.h>
	33	#include <TString.h>
	34
bab17338	35	#include "AliAODEvent.h"
3a72645a	36	#include "AliAODTrack.h"
3a72645a	37	#include "AliESDtrack.h"
3a72645a	38	#include "AliLog.h"
3a72645a	39	#include "AliPID.h"
8c1c76e9	40	#include "AliPIDResponse.h"
3a72645a	41
	42	#include "AliHFEcollection.h"
	43	#include "AliHFEpidBase.h"
6555e2ad	44	#include "AliHFEpidQAmanager.h"
6555e2ad	45	#include "AliHFEpidTPC.h"
3a72645a	46	#include "AliHFEtools.h"
	47	#include "AliHFEtpcPIDqa.h"
	48
6555e2ad	49	ClassImp(AliHFEtpcPIDqa)
6555e2ad	50
3a72645a	51	//_________________________________________________________
	52	AliHFEtpcPIDqa::AliHFEtpcPIDqa():
	53	AliHFEdetPIDqa()
	54	, fHistos(NULL)
e156c3bb	55	, fBrowseCentrality(-1)
3a72645a	56	{
	57	//
	58	// Dummy constructor
	59	//
	60	}
	61
	62	//_________________________________________________________
	63	AliHFEtpcPIDqa::AliHFEtpcPIDqa(const char* name):
	64	AliHFEdetPIDqa(name, "QA for TPC")
	65	, fHistos(NULL)
e156c3bb	66	, fBrowseCentrality(-1)
3a72645a	67	{
	68	//
	69	// Default constructor
	70	//
	71	}
	72
	73	//_________________________________________________________
	74	AliHFEtpcPIDqa::AliHFEtpcPIDqa(const AliHFEtpcPIDqa &o):
	75	AliHFEdetPIDqa(o)
71478e06	76	, fHistos(NULL)
e156c3bb	77	, fBrowseCentrality(o.fBrowseCentrality)
3a72645a	78	{
	79	//
	80	// Copy constructor
	81	//
	82	o.Copy(*this);
	83	}
	84
	85	//_________________________________________________________
	86	AliHFEtpcPIDqa &AliHFEtpcPIDqa::operator=(const AliHFEtpcPIDqa &o){
	87	//
	88	// Do assignment
	89	//
	90	AliHFEdetPIDqa::operator=(o);
	91	if(&o != this) o.Copy(*this);
	92	return *this;
	93	}
	94
	95	//_________________________________________________________
	96	AliHFEtpcPIDqa::~AliHFEtpcPIDqa(){
	97	//
	98	// Destructor
	99	//
	100	if(fHistos) delete fHistos;
	101	}
	102
	103	//_________________________________________________________
	104	void AliHFEtpcPIDqa::Copy(TObject &o) const {
	105	//
	106	// Make copy
	107	//
	108	AliHFEtpcPIDqa &target = dynamic_cast<AliHFEtpcPIDqa &>(o);
	109	if(target.fHistos){
	110	delete target.fHistos;
	111	target.fHistos = NULL;
	112	}
	113	if(fHistos) target.fHistos = new AliHFEcollection(*fHistos);
e156c3bb	114	target.fBrowseCentrality = fBrowseCentrality;
3a72645a	115	}
	116
	117	//_________________________________________________________
	118	Long64_t AliHFEtpcPIDqa::Merge(TCollection *coll){
	119	//
	120	// Merge with other objects
	121	//
	122	if(!coll) return 0;
	123	if(coll->IsEmpty()) return 1;
	124
	125	TIter it(coll);
	126	AliHFEtpcPIDqa *refQA = NULL;
	127	TObject *o = NULL;
	128	Long64_t count = 0;
	129	TList listHistos;
	130	while((o = it())){
	131	refQA = dynamic_cast<AliHFEtpcPIDqa *>(o);
	132	if(!refQA) continue;
	133
	134	listHistos.Add(refQA->fHistos);
	135	count++;
	136	}
	137	fHistos->Merge(&listHistos);
	138	return count + 1;
	139	}
	140
bf892a6a	141	//_________________________________________________________
	142	void AliHFEtpcPIDqa::Browse(TBrowser *b){
	143	//
	144	// Browse the PID QA
	145	//
	146	if(b){
	147	if(fHistos){
	148	b->Add(fHistos, fHistos->GetName());
	149
	150	// Make Projections of the dE/dx Spectra and add them to a new Folder
c2690925	151	TString specnames[AliPID::kSPECIES+1] = {"All", "Electrons", "Muon", "Pions", "Kaon", "Protons"};
c2690925	152	Int_t specind[AliPID::kSPECIES+1] = {-1, AliPID::kElectron, AliPID::kMuon, AliPID::kPion, AliPID::kKaon, AliPID::kProton};
bf892a6a	153	TList *listdEdx = new TList;
	154	listdEdx->SetOwner();
	155	TList *listNsigma = new TList;
	156	listNsigma->SetOwner();
	157
	158	TH2 *hptr = NULL;
c2690925	159	for(Int_t ispec = 0; ispec < AliPID::kSPECIES+1; ispec++){
bf892a6a	160	for(Int_t istep = 0; istep < 2; istep++){
e156c3bb	161	hptr = MakeSpectrumdEdx(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], fBrowseCentrality);
e156c3bb	162	hptr->SetName(Form("hTPCdEdx%s%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After" , fBrowseCentrality == -1 ? "MinBias" : Form("Cent%d", fBrowseCentrality)));
bf892a6a	163	listdEdx->Add(hptr);
e156c3bb	164	hptr = MakeSpectrumNSigma(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], fBrowseCentrality);
e156c3bb	165	hptr->SetName(Form("hTPCnsigma%s%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After", fBrowseCentrality == -1 ? "MinBias" : Form("Cent%d", fBrowseCentrality)));
bf892a6a	166	listNsigma->Add(hptr);
	167	}
	168	}
	169
	170	b->Add(listdEdx, "Projections dE/dx");
	171	b->Add(listNsigma, "Projections NSigma");
	172	}
	173	}
	174	}
	175
3a72645a	176	//_________________________________________________________
	177	void AliHFEtpcPIDqa::Initialize(){
	178	//
	179	// Define Histograms
	180	//
	181
	182	fHistos = new AliHFEcollection("tpcqahistos", "Collection of TPC QA histograms");
	183
	184	// Make common binning
bab17338	185	const Int_t kNdim = 7;
3a72645a	186	const Int_t kPIDbins = AliPID::kSPECIES + 1;
3a72645a	187	const Int_t kSteps = 2;
6555e2ad	188	const Int_t kCentralityBins = 11;
3a72645a	189	const Double_t kMinPID = -1;
	190	const Double_t kMinP = 0.;
	191	const Double_t kMaxPID = (Double_t)AliPID::kSPECIES;
	192	const Double_t kMaxP = 20.;
a8d09079	193	const Double_t kMinEta = -0.9;
	194	const Double_t kMaxEta = 0.9;
	195
bf892a6a	196	// Quantities where one can switch between low and high resolution
8b835997	197	Int_t kPbins = (fQAmanager->HasHighResolutionHistos() \|\| fQAmanager->HasMidResolutionHistos()) ? 1000 : 100;
bf892a6a	198	Int_t kDedxbins = fQAmanager->HasHighResolutionHistos() ? 400 : 200;
bf892a6a	199	Int_t kSigmaBins = fQAmanager->HasHighResolutionHistos() ? 1400 : 240;
bab17338	200	kSigmaBins = fQAmanager->HasMidResolutionHistos() ? 400 : kSigmaBins;
ecaf6269	201	Int_t kEtabins = fQAmanager->HasHighResolutionEtaHistos() ? 33 : 18;
bab17338	202
3a72645a	203	// 1st histogram: TPC dEdx: (species, p, dEdx, step)
ecaf6269	204	Int_t nBinsdEdx[kNdim] = {kPIDbins, kPbins, kDedxbins, kSteps, kCentralityBins, kEtabins, 33};
bab17338	205	Double_t mindEdx[kNdim] = {kMinPID, kMinP, 20., 0., 0., kMinEta, 0};
38be5083	206	Double_t maxdEdx[kNdim] = {kMaxPID, kMaxP, 130, 2., 11., kMaxEta, 2000};
bab17338	207	fHistos->CreateTHnSparse("tpcDedx", "TPC signal; species; p [GeV/c]; TPC signal [a.u.]; Selection Step; Centrality; Eta; pVx contrib", kNdim, nBinsdEdx, mindEdx, maxdEdx);
bab17338	208
3a72645a	209	// 2nd histogram: TPC sigmas: (species, p nsigma, step)
ecaf6269	210	Int_t nBinsSigma[kNdim] = {kPIDbins, kPbins, kSigmaBins, kSteps, kCentralityBins, kEtabins, 33};
bab17338	211	Double_t minSigma[kNdim] = {kMinPID, kMinP, -12., 0., 0., kMinEta, 0};
	212	Double_t maxSigma[kNdim] = {kMaxPID, kMaxP, 12., 2., 11., kMaxEta, 2000};
	213	fHistos->CreateTHnSparse("tpcnSigma", "TPC signal; species; p [GeV/c]; TPC signal [a.u.]; Selection Step; Centrality; Eta; pVx contrib", kNdim, nBinsSigma, minSigma, maxSigma);
3a72645a	214
	215	// General TPC QA
	216	}
	217
	218	//_________________________________________________________
6555e2ad	219	void AliHFEtpcPIDqa::ProcessTrack(const AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step){
3a72645a	220	//
	221	// Fill TPC histograms
	222	//
	223	AliDebug(1, Form("QA started for TPC PID for step %d", (Int_t)step));
6555e2ad	224	AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
6555e2ad	225	Int_t species = track->GetAbInitioPID();
3a72645a	226	if(species >= AliPID::kSPECIES) species = -1;
3a72645a	227
6555e2ad	228	AliHFEpidTPC tpcpid = dynamic_cast<AliHFEpidTPC >(fQAmanager->GetDetectorPID(AliHFEpid::kTPCpid));
8c1c76e9	229	const AliPIDResponse *pidResponse = tpcpid ? tpcpid->GetPIDResponse() : NULL;
d3a15def	230	Double_t contentSignal[7];
3a72645a	231	contentSignal[0] = species;
bf892a6a	232	contentSignal[1] = tpcpid ? tpcpid->GetP(track->GetRecTrack(), anatype) : 0.;
6555e2ad	233	contentSignal[2] = GetTPCsignal(track->GetRecTrack(), anatype);
3a72645a	234	contentSignal[3] = step;
38be5083	235	contentSignal[4] = track->GetCentrality();
a8d09079	236	contentSignal[5] = GetEta(track->GetRecTrack(), anatype);
bab17338	237	contentSignal[6] = fQAmanager->HasFillMultiplicity() ? track->GetMultiplicity() : 0;
bf892a6a	238	fHistos->Fill("tpcDedx", contentSignal);
3a72645a	239
38be5083	240	if(track->HasCorrectedTPCnSigma())
	241	contentSignal[2] = track->GetCorrectedTPCnSigma();
	242	else
	243	contentSignal[2] = pidResponse ? pidResponse->NumberOfSigmasTPC(track->GetRecTrack(), AliPID::kElectron) : -100.;
	244
bf892a6a	245	fHistos->Fill("tpcnSigma", contentSignal);
3a72645a	246	}
	247
	248	//_________________________________________________________
6555e2ad	249	Double_t AliHFEtpcPIDqa::GetTPCsignal(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anatype){
3a72645a	250	//
6555e2ad	251	// Get TPC signal for ESD and AOD track
3a72645a	252	//
6555e2ad	253	Double_t tpcSignal = 0.;
6555e2ad	254	if(anatype == AliHFEpidObject::kESDanalysis){
215ffe88	255	const AliESDtrack esdtrack = static_cast<const AliESDtrack >(track);
6555e2ad	256	tpcSignal = esdtrack->GetTPCsignal();
6555e2ad	257	} else {
215ffe88	258	const AliAODTrack aodtrack = static_cast<const AliAODTrack >(track);
6555e2ad	259	tpcSignal = aodtrack->GetDetPid() ? aodtrack->GetDetPid()->GetTPCsignal() : 0.;
	260	}
	261	return tpcSignal;
3a72645a	262	}
3a72645a	263
a8d09079	264
	265	//_________________________________________________________
	266	Double_t AliHFEtpcPIDqa::GetEta(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anatype){
	267	//
	268	// Get TPC signal for ESD and AOD track
	269	//
	270	Double_t eta = 0.;
	271	if(anatype == AliHFEpidObject::kESDanalysis){
	272	const AliESDtrack esdtrack = static_cast<const AliESDtrack >(track);
	273	eta = esdtrack->Eta();
	274	} else {
	275	const AliAODTrack aodtrack = static_cast<const AliAODTrack >(track);
	276	eta = aodtrack->Eta();
	277	}
	278	return eta;
	279	}
	280
3a72645a	281	//_________________________________________________________
e156c3bb	282	TH2 *AliHFEtpcPIDqa::MakeSpectrumdEdx(AliHFEdetPIDqa::EStep_t istep, Int_t species, Int_t centralityClass){
3a72645a	283	//
	284	// Plot the Spectrum
	285	//
	286	THnSparseF hSignal = dynamic_cast<THnSparseF >(fHistos->Get("tpcDedx"));
bf892a6a	287	if(!hSignal) return NULL;
3a72645a	288	hSignal->GetAxis(3)->SetRange(istep + 1, istep + 1);
bf892a6a	289	if(species >= 0 && species < AliPID::kSPECIES)
3a72645a	290	hSignal->GetAxis(0)->SetRange(2 + species, 2 + species);
bf892a6a	291	TString hname = Form("hTPCsignal%s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after"),
bf892a6a	292	htitle = Form("TPC dE/dx Spectrum %s selection", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
3a72645a	293	if(species > -1){
bf892a6a	294	hname += AliPID::ParticleName(species);
bf892a6a	295	htitle += Form(" for %ss", AliPID::ParticleName(species));
3a72645a	296	}
e156c3bb	297	TString centname, centtitle;
	298	Bool_t hasCentralityInfo = kTRUE;
	299	if(centralityClass > -1){
	300	if(hSignal->GetNdimensions() < 5){
	301	AliError("Centrality Information not available");
	302	centname = centtitle = "MinBias";
	303	hasCentralityInfo= kFALSE;
	304	} else {
	305	// Project centrality classes
	306	// -1 is Min. Bias
	307	hSignal->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
	308	centname = Form("Cent%d", centralityClass);
	309	centtitle = Form(" Centrality %d", centralityClass);
	310	}
	311	} else {
	312	centname = centtitle = "MinBias";
	313	hasCentralityInfo= kFALSE;
	314	}
	315	hname += centtitle;
	316	htitle += centtitle;
	317
	318	TH2 *hTmp = hSignal->Projection(2,1);
bf892a6a	319	hTmp->SetName(hname.Data());
bf892a6a	320	hTmp->SetTitle(htitle.Data());
3a72645a	321	hTmp->SetStats(kFALSE);
	322	hTmp->GetXaxis()->SetTitle("p [GeV/c]");
	323	hTmp->GetYaxis()->SetTitle("TPC signal [a.u.]");
	324	hSignal->GetAxis(3)->SetRange(0, hSignal->GetAxis(3)->GetNbins());
	325	hSignal->GetAxis(0)->SetRange(0, hSignal->GetAxis(0)->GetNbins());
e156c3bb	326	if(hasCentralityInfo) hSignal->GetAxis(4)->SetRange(0, hSignal->GetAxis(4)->GetNbins());
3a72645a	327	return hTmp;
	328	}
	329
	330	//_________________________________________________________
e156c3bb	331	TH2 *AliHFEtpcPIDqa::MakeSpectrumNSigma(AliHFEdetPIDqa::EStep_t istep, Int_t species, Int_t centralityClass){
3a72645a	332	//
	333	// Plot the Spectrum
	334	//
	335	THnSparseF hSignal = dynamic_cast<THnSparseF >(fHistos->Get("tpcnSigma"));
bf892a6a	336	if(!hSignal) return NULL;
3a72645a	337	hSignal->GetAxis(3)->SetRange(istep + 1, istep + 1);
	338	if(species >= 0 && species < AliPID::kSPECIES)
	339	hSignal->GetAxis(0)->SetRange(2 + species, 2 + species);
bf892a6a	340	TString hname = Form("hTPCsigma%s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after"),
bf892a6a	341	htitle = Form("TPC dE/dx Spectrum[#sigma] %s selection", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
3a72645a	342	if(species > -1){
bf892a6a	343	hname += AliPID::ParticleName(species);
bf892a6a	344	htitle += Form(" for %ss", AliPID::ParticleName(species));
3a72645a	345	}
e156c3bb	346	TString centname, centtitle;
	347	Bool_t hasCentralityInfo = kTRUE;
	348	if(centralityClass > -1){
	349	if(hSignal->GetNdimensions() < 5){
	350	AliError("Centrality Information not available");
	351	centname = centtitle = "MinBias";
	352	hasCentralityInfo= kFALSE;
	353	} else {
	354	// Project centrality classes
	355	// -1 is Min. Bias
	356	hSignal->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
	357	centname = Form("Cent%d", centralityClass);
	358	centtitle = Form(" Centrality %d", centralityClass);
	359	}
	360	} else {
	361	centname = centtitle = "MinBias";
	362	hasCentralityInfo= kFALSE;
	363	}
	364	hname += centtitle;
	365	htitle += centtitle;
	366
	367	TH2 *hTmp = hSignal->Projection(2,1);
bf892a6a	368	hTmp->SetName(hname.Data());
bf892a6a	369	hTmp->SetTitle(htitle.Data());
3a72645a	370	hTmp->SetStats(kFALSE);
	371	hTmp->GetXaxis()->SetTitle("p [GeV/c]");
	372	hTmp->GetYaxis()->SetTitle("TPC dE/dx - <dE/dx>\|_{el} [#sigma]");
	373	hSignal->GetAxis(3)->SetRange(0, hSignal->GetAxis(3)->GetNbins());
	374	hSignal->GetAxis(0)->SetRange(0, hSignal->GetAxis(0)->GetNbins());
e156c3bb	375	if(hasCentralityInfo) hSignal->GetAxis(4)->SetRange(0, hSignal->GetAxis(4)->GetNbins());
3a72645a	376	return hTmp;
	377	}
	378