[u/mrichter/AliRoot.git] / PWG3 / 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 <TClass.h>
#include <TH2.h>
#include <THnSparse.h>
#include <TString.h>

#include "AliAODTrack.h"
#include "AliESDtrack.h"
#include "AliESDpid.h"
#include "AliLog.h"
#include "AliPID.h"

#include "AliHFEcollection.h"
#include "AliHFEpidBase.h"
#include "AliHFEtools.h"
#include "AliHFEtpcPIDqa.h"

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

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

//_________________________________________________________
AliHFEtpcPIDqa::AliHFEtpcPIDqa(const AliHFEtpcPIDqa &o):
    AliHFEdetPIDqa(o)
  , fHistos()
{
  //
  // 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);
}

//_________________________________________________________
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::Initialize(){
  //
  // Define Histograms
  //

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

  // Make common binning
  const Int_t kNdim = 5;
  const Int_t kPIDbins = AliPID::kSPECIES + 1;
  const Int_t kPbins = 1000;
  const Int_t kSteps = 2;
  const Int_t kCentralityBins = 20;
  const Double_t kMinPID = -1;
  const Double_t kMinP = 0.;
  const Double_t kMaxPID = (Double_t)AliPID::kSPECIES;
  const Double_t kMaxP = 20.;
  
  // 1st histogram: TPC dEdx: (species, p, dEdx, step)
  const Int_t kDedxbins = 600;
  Int_t nBinsdEdx[kNdim] = {kPIDbins, kPbins, kDedxbins, kSteps, kCentralityBins};
  Double_t mindEdx[kNdim] =  {kMinPID, kMinP, 0., 0., 0.};
  Double_t maxdEdx[kNdim] =  {kMaxPID, kMaxP, 300, 2., 100.}; 
  fHistos->CreateTHnSparse("tpcDedx", "TPC signal; species; p [GeV/c]; TPC signal [a.u.]; Centrality; Selection Step", kNdim, nBinsdEdx, mindEdx, maxdEdx);
  // 2nd histogram: TPC sigmas: (species, p nsigma, step)
  const Int_t kSigmaBins = 1400;
  Int_t nBinsSigma[kNdim] = {kPIDbins, kPbins, kSigmaBins, kSteps, kCentralityBins};
  Double_t minSigma[kNdim] = {kMinPID, kMinP, -12., 0., 0.};
  Double_t maxSigma[kNdim] = {kMaxPID, kMaxP, 12., 2., 100.};
  fHistos->CreateTHnSparse("tpcnSigma", "TPC signal; species; p [GeV/c]; TPC signal [a.u.]; Centrality; Selection Step", kNdim, nBinsSigma, minSigma, maxSigma);

  // General TPC QA
}

//_________________________________________________________
void AliHFEtpcPIDqa::ProcessTrack(AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step){
  //
  // Fill TPC histograms
  //
  AliDebug(1, Form("QA started for TPC PID for step %d", (Int_t)step));
  Int_t species = track->GetAbInitioPID();
  Float_t centrality = track->GetCentrality();
  const AliVParticle *rectrack = track->GetRecTrack();
  if(species >= AliPID::kSPECIES) species = -1;
  if(!TString(rectrack->IsA()->GetName()).CompareTo("AliESDtrack")) ProcessESDtrack(dynamic_cast<const AliESDtrack *>(rectrack), step, species, centrality);
  else if(!TString(rectrack->IsA()->GetName()).CompareTo("AliAODTrack")) ProcessAODtrack(dynamic_cast<const AliAODTrack *>(rectrack), step, species, centrality);
  else  AliWarning(Form("Object type %s not supported\n", rectrack->IsA()->GetName()));
}

//_________________________________________________________
void AliHFEtpcPIDqa::ProcessESDtrack(const AliESDtrack *track, AliHFEdetPIDqa::EStep_t step, Int_t species, Float_t centrality){
  //
  // Process track as ESD track
  //
  if(!fESDpid){
    AliError("No ESD PID object available");
    return;
  }
  AliDebug(1, Form("Monitoring particle of type %d for step %d", species, step));
  
  Double_t contentSignal[5];
  contentSignal[0] = species;
  contentSignal[1] = track->GetInnerParam() ? track->GetInnerParam()->P() : track->P();
  contentSignal[2] = track->GetTPCsignal();
  contentSignal[3] = step;
  contentSignal[4] = centrality;
  (dynamic_cast<THnSparseF *>(fHistos->Get("tpcDedx")))->Fill(contentSignal);

  contentSignal[2] = fESDpid->NumberOfSigmasTPC(track, AliPID::kElectron); 
  (dynamic_cast<THnSparseF *>(fHistos->Get("tpcnSigma")))->Fill(contentSignal);
}

//_________________________________________________________
void AliHFEtpcPIDqa::ProcessAODtrack(const AliAODTrack * /*track*/, AliHFEdetPIDqa::EStep_t /*step*/, Int_t /*species*/, Float_t /*centrality*/){
  //
  // Process track as AOD track
  //
  AliInfo("Method implemented soon!");
}

//_________________________________________________________
TH2 *AliHFEtpcPIDqa::MakeSpectrumdEdx(AliHFEdetPIDqa::EStep_t istep, Int_t species){
  //
  // Plot the Spectrum
  //
  THnSparseF *hSignal = dynamic_cast<THnSparseF *>(fHistos->Get("tpcDedx"));
  hSignal->GetAxis(3)->SetRange(istep + 1, istep + 1);
  if(species > 0 && species < AliPID::kSPECIES)
    hSignal->GetAxis(0)->SetRange(2 + species, 2 + species);
  TH2 *hTmp = hSignal->Projection(2,1);
  Char_t hname[256], htitle[256];
  sprintf(hname, "hTPCsignal%s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
  sprintf(htitle, "TPC dE/dx Spectrum %s selection", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
  if(species > -1){
    strncat(hname, AliPID::ParticleName(species), strlen(AliPID::ParticleName(species)));
     Char_t speciesname[256];
     sprintf(speciesname, " for %ss", AliPID::ParticleName(species));
     strncat(htitle, speciesname, strlen(speciesname));
  }
  hTmp->SetName(hname);
  hTmp->SetTitle(htitle);
  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());
  return hTmp;
}

//_________________________________________________________
TH2 *AliHFEtpcPIDqa::MakeSpectrumNSigma(AliHFEdetPIDqa::EStep_t istep, Int_t species){
  //
  // Plot the Spectrum
  //
  THnSparseF *hSignal = dynamic_cast<THnSparseF *>(fHistos->Get("tpcnSigma"));
  hSignal->GetAxis(3)->SetRange(istep + 1, istep + 1);
  if(species >= 0 && species < AliPID::kSPECIES)
    hSignal->GetAxis(0)->SetRange(2 + species, 2 + species);
  TH2 *hTmp = hSignal->Projection(2,1);
  Char_t hname[256], htitle[256];
  sprintf(hname, "hTPCsigma%s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
  sprintf(htitle, "TPC dE/dx Spectrum[#sigma] %s selection", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
  if(species > -1){
    strncat(hname, AliPID::ParticleName(species), strlen(AliPID::ParticleName(species)));
     Char_t speciesname[256];
     sprintf(speciesname, " for %ss", AliPID::ParticleName(species));
     strncat(htitle, speciesname, strlen(speciesname));
  }
  hTmp->SetName(hname);
  hTmp->SetTitle(htitle);
  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());
  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	//
	29	#include <TClass.h>
	30	#include <TH2.h>
	31	#include <THnSparse.h>
	32	#include <TString.h>
	33
	34	#include "AliAODTrack.h"
	35	#include "AliESDtrack.h"
	36	#include "AliESDpid.h"
	37	#include "AliLog.h"
	38	#include "AliPID.h"
	39
	40	#include "AliHFEcollection.h"
	41	#include "AliHFEpidBase.h"
	42	#include "AliHFEtools.h"
	43	#include "AliHFEtpcPIDqa.h"
	44
	45	//_________________________________________________________
	46	AliHFEtpcPIDqa::AliHFEtpcPIDqa():
	47	AliHFEdetPIDqa()
	48	, fHistos(NULL)
	49	{
	50	//
	51	// Dummy constructor
	52	//
	53	}
	54
	55	//_________________________________________________________
	56	AliHFEtpcPIDqa::AliHFEtpcPIDqa(const char* name):
	57	AliHFEdetPIDqa(name, "QA for TPC")
	58	, fHistos(NULL)
	59	{
	60	//
	61	// Default constructor
	62	//
	63	}
	64
65	//_________________________________________________________
66	AliHFEtpcPIDqa::AliHFEtpcPIDqa(const AliHFEtpcPIDqa &o):
67	AliHFEdetPIDqa(o)
68	, fHistos()
69	{
70	//
71	// Copy constructor
72	//
73	o.Copy(*this);
74	}
75
76	//_________________________________________________________
77	AliHFEtpcPIDqa &AliHFEtpcPIDqa::operator=(const AliHFEtpcPIDqa &o){
78	//
79	// Do assignment
80	//
81	AliHFEdetPIDqa::operator=(o);
82	if(&o != this) o.Copy(*this);
83	return *this;
84	}
85
86	//_________________________________________________________
87	AliHFEtpcPIDqa::~AliHFEtpcPIDqa(){
88	//
89	// Destructor
90	//
91	if(fHistos) delete fHistos;
92	}
93
94	//_________________________________________________________
95	void AliHFEtpcPIDqa::Copy(TObject &o) const {
96	//
97	// Make copy
98	//
99	AliHFEtpcPIDqa &target = dynamic_cast<AliHFEtpcPIDqa &>(o);
100	if(target.fHistos){
101	delete target.fHistos;
102	target.fHistos = NULL;
103	}
104	if(fHistos) target.fHistos = new AliHFEcollection(*fHistos);
105	}
106
107	//_________________________________________________________
108	Long64_t AliHFEtpcPIDqa::Merge(TCollection *coll){
109	//
110	// Merge with other objects
111	//
112	if(!coll) return 0;
113	if(coll->IsEmpty()) return 1;
114
115	TIter it(coll);
116	AliHFEtpcPIDqa *refQA = NULL;
117	TObject *o = NULL;
118	Long64_t count = 0;
119	TList listHistos;
120	while((o = it())){
121	refQA = dynamic_cast<AliHFEtpcPIDqa *>(o);
122	if(!refQA) continue;
123
124	listHistos.Add(refQA->fHistos);
125	count++;
126	}
127	fHistos->Merge(&listHistos);
128	return count + 1;
129	}
130
131	//_________________________________________________________
132	void AliHFEtpcPIDqa::Initialize(){
133	//
134	// Define Histograms
135	//
136
137	fHistos = new AliHFEcollection("tpcqahistos", "Collection of TPC QA histograms");
138
139	// Make common binning
140	const Int_t kNdim = 5;
141	const Int_t kPIDbins = AliPID::kSPECIES + 1;
142	const Int_t kPbins = 1000;
143	const Int_t kSteps = 2;
144	const Int_t kCentralityBins = 20;
145	const Double_t kMinPID = -1;
146	const Double_t kMinP = 0.;
147	const Double_t kMaxPID = (Double_t)AliPID::kSPECIES;
148	const Double_t kMaxP = 20.;
149
150	// 1st histogram: TPC dEdx: (species, p, dEdx, step)
151	const Int_t kDedxbins = 600;
152	Int_t nBinsdEdx[kNdim] = {kPIDbins, kPbins, kDedxbins, kSteps, kCentralityBins};
153	Double_t mindEdx[kNdim] = {kMinPID, kMinP, 0., 0., 0.};
154	Double_t maxdEdx[kNdim] = {kMaxPID, kMaxP, 300, 2., 100.};
155	fHistos->CreateTHnSparse("tpcDedx", "TPC signal; species; p [GeV/c]; TPC signal [a.u.]; Centrality; Selection Step", kNdim, nBinsdEdx, mindEdx, maxdEdx);
156	// 2nd histogram: TPC sigmas: (species, p nsigma, step)
157	const Int_t kSigmaBins = 1400;
158	Int_t nBinsSigma[kNdim] = {kPIDbins, kPbins, kSigmaBins, kSteps, kCentralityBins};
159	Double_t minSigma[kNdim] = {kMinPID, kMinP, -12., 0., 0.};
160	Double_t maxSigma[kNdim] = {kMaxPID, kMaxP, 12., 2., 100.};
161	fHistos->CreateTHnSparse("tpcnSigma", "TPC signal; species; p [GeV/c]; TPC signal [a.u.]; Centrality; Selection Step", kNdim, nBinsSigma, minSigma, maxSigma);
162
163	// General TPC QA
164	}
165
166	//_________________________________________________________
167	void AliHFEtpcPIDqa::ProcessTrack(AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step){
168	//
169	// Fill TPC histograms
170	//
171	AliDebug(1, Form("QA started for TPC PID for step %d", (Int_t)step));
172	Int_t species = track->GetAbInitioPID();
173	Float_t centrality = track->GetCentrality();
174	const AliVParticle *rectrack = track->GetRecTrack();
175	if(species >= AliPID::kSPECIES) species = -1;
176	if(!TString(rectrack->IsA()->GetName()).CompareTo("AliESDtrack")) ProcessESDtrack(dynamic_cast<const AliESDtrack *>(rectrack), step, species, centrality);
177	else if(!TString(rectrack->IsA()->GetName()).CompareTo("AliAODTrack")) ProcessAODtrack(dynamic_cast<const AliAODTrack *>(rectrack), step, species, centrality);
178	else AliWarning(Form("Object type %s not supported\n", rectrack->IsA()->GetName()));
179	}
180
181	//_________________________________________________________
182	void AliHFEtpcPIDqa::ProcessESDtrack(const AliESDtrack *track, AliHFEdetPIDqa::EStep_t step, Int_t species, Float_t centrality){
183	//
184	// Process track as ESD track
185	//
186	if(!fESDpid){
187	AliError("No ESD PID object available");
188	return;
189	}
190	AliDebug(1, Form("Monitoring particle of type %d for step %d", species, step));
191
192	Double_t contentSignal[5];
193	contentSignal[0] = species;
194	contentSignal[1] = track->GetInnerParam() ? track->GetInnerParam()->P() : track->P();
195	contentSignal[2] = track->GetTPCsignal();
196	contentSignal[3] = step;
197	contentSignal[4] = centrality;
198	(dynamic_cast<THnSparseF *>(fHistos->Get("tpcDedx")))->Fill(contentSignal);
199
200	contentSignal[2] = fESDpid->NumberOfSigmasTPC(track, AliPID::kElectron);
201	(dynamic_cast<THnSparseF *>(fHistos->Get("tpcnSigma")))->Fill(contentSignal);
202	}
203
204	//_________________________________________________________
205	void AliHFEtpcPIDqa::ProcessAODtrack(const AliAODTrack * /track/, AliHFEdetPIDqa::EStep_t /step/, Int_t /species/, Float_t /centrality/){
206	//
207	// Process track as AOD track
208	//
209	AliInfo("Method implemented soon!");
210	}
211
212	//_________________________________________________________
213	TH2 *AliHFEtpcPIDqa::MakeSpectrumdEdx(AliHFEdetPIDqa::EStep_t istep, Int_t species){
214	//
215	// Plot the Spectrum
216	//
217	THnSparseF hSignal = dynamic_cast<THnSparseF >(fHistos->Get("tpcDedx"));
218	hSignal->GetAxis(3)->SetRange(istep + 1, istep + 1);
219	if(species > 0 && species < AliPID::kSPECIES)
220	hSignal->GetAxis(0)->SetRange(2 + species, 2 + species);
221	TH2 *hTmp = hSignal->Projection(2,1);
222	Char_t hname[256], htitle[256];
223	sprintf(hname, "hTPCsignal%s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
224	sprintf(htitle, "TPC dE/dx Spectrum %s selection", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
225	if(species > -1){
226	strncat(hname, AliPID::ParticleName(species), strlen(AliPID::ParticleName(species)));
227	Char_t speciesname[256];
228	sprintf(speciesname, " for %ss", AliPID::ParticleName(species));
229	strncat(htitle, speciesname, strlen(speciesname));
230	}
231	hTmp->SetName(hname);
232	hTmp->SetTitle(htitle);
233	hTmp->SetStats(kFALSE);
234	hTmp->GetXaxis()->SetTitle("p [GeV/c]");
235	hTmp->GetYaxis()->SetTitle("TPC signal [a.u.]");
236	hSignal->GetAxis(3)->SetRange(0, hSignal->GetAxis(3)->GetNbins());
237	hSignal->GetAxis(0)->SetRange(0, hSignal->GetAxis(0)->GetNbins());
238	return hTmp;
239	}
240
241	//_________________________________________________________
242	TH2 *AliHFEtpcPIDqa::MakeSpectrumNSigma(AliHFEdetPIDqa::EStep_t istep, Int_t species){
243	//
244	// Plot the Spectrum
245	//
246	THnSparseF hSignal = dynamic_cast<THnSparseF >(fHistos->Get("tpcnSigma"));
247	hSignal->GetAxis(3)->SetRange(istep + 1, istep + 1);
248	if(species >= 0 && species < AliPID::kSPECIES)
249	hSignal->GetAxis(0)->SetRange(2 + species, 2 + species);
250	TH2 *hTmp = hSignal->Projection(2,1);
251	Char_t hname[256], htitle[256];
252	sprintf(hname, "hTPCsigma%s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
253	sprintf(htitle, "TPC dE/dx Spectrum[#sigma] %s selection", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
254	if(species > -1){
255	strncat(hname, AliPID::ParticleName(species), strlen(AliPID::ParticleName(species)));
256	Char_t speciesname[256];
257	sprintf(speciesname, " for %ss", AliPID::ParticleName(species));
258	strncat(htitle, speciesname, strlen(speciesname));
259	}
260	hTmp->SetName(hname);
261	hTmp->SetTitle(htitle);
262	hTmp->SetStats(kFALSE);
263	hTmp->GetXaxis()->SetTitle("p [GeV/c]");
264	hTmp->GetYaxis()->SetTitle("TPC dE/dx - <dE/dx>\|_{el} [#sigma]");
265	hSignal->GetAxis(3)->SetRange(0, hSignal->GetAxis(3)->GetNbins());
266	hSignal->GetAxis(0)->SetRange(0, hSignal->GetAxis(0)->GetNbins());
267	return hTmp;
268	}
269