[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEitsPIDqa.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 AliHFEitsPIDqa
// Monitoring ITS PID in the HFE PID montioring framework. The following
// quantities are monitored:
//   ITS dE/dx (Number of sigmas)
// (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:
//    Yvonne Pachmayer <pachmay@physi.uni-heidelberg.de>
//
#include <TBrowser.h>
#include <TClass.h>
#include <TH2.h>
#include <THnSparse.h>
#include <TString.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 "AliHFEpidITS.h"
#include "AliHFEtools.h"
#include "AliHFEitsPIDqa.h"

ClassImp(AliHFEitsPIDqa)

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

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

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

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

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

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

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

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

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

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

      // Make Projections of the nsigma 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 *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 = MakeSpectrumNSigma(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], fBrowseCentrality);
          hptr->SetName(Form("hITSnsigma%s%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After", fBrowseCentrality == -1 ? "MinBias" : Form("Cent%d", fBrowseCentrality)));
          listNsigma->Add(hptr);
        }
      }
      
      b->Add(listNsigma, "Projections NSigma");
    }
  }
}

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

  fHistos = new AliHFEcollection("itsqahistos", "Collection of ITS QA histograms");

  // Make common binning
  const Int_t kNdim = 5;
  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 = 2.;

  // Quantities where one can switch between low and high resolution
  Int_t kPbins = fQAmanager->HasHighResolutionHistos() ? 100 : 10;
  Int_t kSigmaBins = fQAmanager->HasHighResolutionHistos() ? 1400 : 240;
 
  // 1st histogram: ITS sigmas: (species, p, nsigma, step, centrality)
  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., 11.};
  fHistos->CreateTHnSparse("itsnSigma", "ITS signal; species; p [GeV/c]; ITS signal [a.u.]; Selection Step; Centrality; Eta", kNdim, nBinsSigma, minSigma, maxSigma);

  // General ITS QA
}

//_________________________________________________________
void AliHFEitsPIDqa::ProcessTrack(const AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step){
  //
  // Fill ITS histograms
  //
  AliDebug(1, Form("QA started for ITS PID for step %d", (Int_t)step));
  //AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
  Float_t centrality = track->GetCentrality();
  Int_t species = track->GetAbInitioPID();
  if(species >= AliPID::kSPECIES) species = -1;

  AliHFEpidITS *itspid = dynamic_cast<AliHFEpidITS *>(fQAmanager->GetDetectorPID(AliHFEpid::kITSpid));
  Double_t contentSignal[5];
  contentSignal[0] = species;
  contentSignal[1] = track->GetRecTrack()->P();
  contentSignal[2] = itspid ? itspid->GetITSNsigmaCorrected(track->GetRecTrack()) : -100.; 
  contentSignal[3] = step;
  contentSignal[4] = centrality;
 
  fHistos->Fill("itsnSigma", contentSignal);
}

//_________________________________________________________
TH2 *AliHFEitsPIDqa::MakeSpectrumNSigma(AliHFEdetPIDqa::EStep_t istep, Int_t species, Int_t centralityClass){
  //
  // Plot the Spectrum
  //
  THnSparseF *hSignal = dynamic_cast<THnSparseF *>(fHistos->Get("itsnSigma"));
  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("hITSsigma%s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after"), 
          htitle = Form("ITS 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("ITS 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
4437a0d2	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 AliHFEitsPIDqa
	17	// Monitoring ITS PID in the HFE PID montioring framework. The following
	18	// quantities are monitored:
	19	// ITS dE/dx (Number of sigmas)
	20	// (Always as function of momentum, particle species and centrality
	21	// before and after cut)
	22	// More information about the PID monitoring framework can be found in
	23	// AliHFEpidQAmanager.cxx and AliHFEdetPIDqa.cxx
	24	//
	25	// Author:
	26	// Yvonne Pachmayer <pachmay@physi.uni-heidelberg.de>
	27	//
	28	#include <TBrowser.h>
	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 "AliLog.h"
	37	#include "AliPID.h"
	38	#include "AliPIDResponse.h"
	39
	40	#include "AliHFEcollection.h"
	41	#include "AliHFEpidBase.h"
	42	#include "AliHFEpidQAmanager.h"
	43	#include "AliHFEpidITS.h"
	44	#include "AliHFEtools.h"
	45	#include "AliHFEitsPIDqa.h"
	46
	47	ClassImp(AliHFEitsPIDqa)
	48
	49	//_________________________________________________________
	50	AliHFEitsPIDqa::AliHFEitsPIDqa():
	51	AliHFEdetPIDqa()
	52	, fHistos(NULL)
	53	, fBrowseCentrality(-1)
	54	{
	55	//
	56	// Dummy constructor
	57	//
	58	}
	59
	60	//_________________________________________________________
	61	AliHFEitsPIDqa::AliHFEitsPIDqa(const char* name):
	62	AliHFEdetPIDqa(name, "QA for ITS")
	63	, fHistos(NULL)
	64	, fBrowseCentrality(-1)
65	{
66	//
67	// Default constructor
68	//
69	}
70
71	//_________________________________________________________
72	AliHFEitsPIDqa::AliHFEitsPIDqa(const AliHFEitsPIDqa &o):
73	AliHFEdetPIDqa(o)
74	, fHistos(NULL)
75	, fBrowseCentrality(o.fBrowseCentrality)
76	{
77	//
78	// Copy constructor
79	//
80	o.Copy(*this);
81	}
82
83	//_________________________________________________________
84	AliHFEitsPIDqa &AliHFEitsPIDqa::operator=(const AliHFEitsPIDqa &o){
85	//
86	// Do assignment
87	//
88	AliHFEdetPIDqa::operator=(o);
89	if(&o != this) o.Copy(*this);
90	return *this;
91	}
92
93	//_________________________________________________________
94	AliHFEitsPIDqa::~AliHFEitsPIDqa(){
95	//
96	// Destructor
97	//
98	if(fHistos) delete fHistos;
99	}
100
101	//_________________________________________________________
102	void AliHFEitsPIDqa::Copy(TObject &o) const {
103	//
104	// Make copy
105	//
106	AliHFEitsPIDqa &target = dynamic_cast<AliHFEitsPIDqa &>(o);
107	if(target.fHistos){
108	delete target.fHistos;
109	target.fHistos = NULL;
110	}
111	if(fHistos) target.fHistos = new AliHFEcollection(*fHistos);
112	target.fBrowseCentrality = fBrowseCentrality;
113	}
114
115	//_________________________________________________________
116	Long64_t AliHFEitsPIDqa::Merge(TCollection *coll){
117	//
118	// Merge with other objects
119	//
120	if(!coll) return 0;
121	if(coll->IsEmpty()) return 1;
122
123	TIter it(coll);
124	AliHFEitsPIDqa *refQA = NULL;
125	TObject *o = NULL;
126	Long64_t count = 0;
127	TList listHistos;
128	while((o = it())){
129	refQA = dynamic_cast<AliHFEitsPIDqa *>(o);
130	if(!refQA) continue;
131
132	listHistos.Add(refQA->fHistos);
133	count++;
134	}
135	fHistos->Merge(&listHistos);
136	return count + 1;
137	}
138
139	//_________________________________________________________
140	void AliHFEitsPIDqa::Browse(TBrowser *b){
141	//
142	// Browse the PID QA
143	//
144	if(b){
145	if(fHistos){
146	b->Add(fHistos, fHistos->GetName());
147
148	// Make Projections of the nsigma Spectra and add them to a new Folder
149	TString specnames[AliPID::kSPECIES+1] = {"All", "Electrons", "Muon", "Pions", "Kaon", "Protons"};
150	Int_t specind[AliPID::kSPECIES+1] = {-1, AliPID::kElectron, AliPID::kMuon, AliPID::kPion, AliPID::kKaon, AliPID::kProton};
151	TList *listNsigma = new TList;
152	listNsigma->SetOwner();
153
154	TH2 *hptr = NULL;
155	for(Int_t ispec = 0; ispec < AliPID::kSPECIES+1; ispec++){
156	for(Int_t istep = 0; istep < 2; istep++){
157	hptr = MakeSpectrumNSigma(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], fBrowseCentrality);
158	hptr->SetName(Form("hITSnsigma%s%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After", fBrowseCentrality == -1 ? "MinBias" : Form("Cent%d", fBrowseCentrality)));
159	listNsigma->Add(hptr);
160	}
161	}
162
163	b->Add(listNsigma, "Projections NSigma");
164	}
165	}
166	}
167
168	//_________________________________________________________
169	void AliHFEitsPIDqa::Initialize(){
170	//
171	// Define Histograms
172	//
173
174	fHistos = new AliHFEcollection("itsqahistos", "Collection of ITS QA histograms");
175
176	// Make common binning
177	const Int_t kNdim = 5;
178	const Int_t kPIDbins = AliPID::kSPECIES + 1;
179	const Int_t kSteps = 2;
180	const Int_t kCentralityBins = 11;
181	const Double_t kMinPID = -1;
182	const Double_t kMinP = 0.;
183	const Double_t kMaxPID = (Double_t)AliPID::kSPECIES;
184	const Double_t kMaxP = 2.;
185
186	// Quantities where one can switch between low and high resolution
187	Int_t kPbins = fQAmanager->HasHighResolutionHistos() ? 100 : 10;
188	Int_t kSigmaBins = fQAmanager->HasHighResolutionHistos() ? 1400 : 240;
189
190	// 1st histogram: ITS sigmas: (species, p, nsigma, step, centrality)
191	Int_t nBinsSigma[kNdim] = {kPIDbins, kPbins, kSigmaBins, kSteps, kCentralityBins};
192	Double_t minSigma[kNdim] = {kMinPID, kMinP, -12., 0., 0.};
193	Double_t maxSigma[kNdim] = {kMaxPID, kMaxP, 12., 2., 11.};
194	fHistos->CreateTHnSparse("itsnSigma", "ITS signal; species; p [GeV/c]; ITS signal [a.u.]; Selection Step; Centrality; Eta", kNdim, nBinsSigma, minSigma, maxSigma);
195
196	// General ITS QA
197	}
198
199	//_________________________________________________________
200	void AliHFEitsPIDqa::ProcessTrack(const AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step){
201	//
202	// Fill ITS histograms
203	//
204	AliDebug(1, Form("QA started for ITS PID for step %d", (Int_t)step));
f06b970d	205	//AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
4437a0d2	206	Float_t centrality = track->GetCentrality();
	207	Int_t species = track->GetAbInitioPID();
	208	if(species >= AliPID::kSPECIES) species = -1;
	209
	210	AliHFEpidITS itspid = dynamic_cast<AliHFEpidITS >(fQAmanager->GetDetectorPID(AliHFEpid::kITSpid));
4437a0d2	211	Double_t contentSignal[5];
	212	contentSignal[0] = species;
	213	contentSignal[1] = track->GetRecTrack()->P();
f06b970d	214	contentSignal[2] = itspid ? itspid->GetITSNsigmaCorrected(track->GetRecTrack()) : -100.;
4437a0d2	215	contentSignal[3] = step;
	216	contentSignal[4] = centrality;
	217
	218	fHistos->Fill("itsnSigma", contentSignal);
	219	}
	220
	221	//_________________________________________________________
	222	TH2 *AliHFEitsPIDqa::MakeSpectrumNSigma(AliHFEdetPIDqa::EStep_t istep, Int_t species, Int_t centralityClass){
	223	//
	224	// Plot the Spectrum
	225	//
	226	THnSparseF hSignal = dynamic_cast<THnSparseF >(fHistos->Get("itsnSigma"));
	227	if(!hSignal) return NULL;
	228	hSignal->GetAxis(3)->SetRange(istep + 1, istep + 1);
	229	if(species >= 0 && species < AliPID::kSPECIES)
	230	hSignal->GetAxis(0)->SetRange(2 + species, 2 + species);
	231	TString hname = Form("hITSsigma%s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after"),
	232	htitle = Form("ITS dE/dx Spectrum[#sigma] %s selection", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
	233	if(species > -1){
	234	hname += AliPID::ParticleName(species);
	235	htitle += Form(" for %ss", AliPID::ParticleName(species));
	236	}
	237	TString centname, centtitle;
	238	Bool_t hasCentralityInfo = kTRUE;
	239	if(centralityClass > -1){
	240	if(hSignal->GetNdimensions() < 5){
	241	AliError("Centrality Information not available");
	242	centname = centtitle = "MinBias";
	243	hasCentralityInfo= kFALSE;
	244	} else {
	245	// Project centrality classes
	246	// -1 is Min. Bias
	247	hSignal->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
	248	centname = Form("Cent%d", centralityClass);
	249	centtitle = Form(" Centrality %d", centralityClass);
	250	}
	251	} else {
	252	centname = centtitle = "MinBias";
	253	hasCentralityInfo= kFALSE;
	254	}
	255	hname += centtitle;
	256	htitle += centtitle;
	257
	258	TH2 *hTmp = hSignal->Projection(2,1);
	259	hTmp->SetName(hname.Data());
	260	hTmp->SetTitle(htitle.Data());
	261	hTmp->SetStats(kFALSE);
	262	hTmp->GetXaxis()->SetTitle("p [GeV/c]");
	263	hTmp->GetYaxis()->SetTitle("ITS dE/dx - <dE/dx>\|_{el} [#sigma]");
	264	hSignal->GetAxis(3)->SetRange(0, hSignal->GetAxis(3)->GetNbins());
	265	hSignal->GetAxis(0)->SetRange(0, hSignal->GetAxis(0)->GetNbins());
	266	if(hasCentralityInfo) hSignal->GetAxis(4)->SetRange(0, hSignal->GetAxis(4)->GetNbins());
	267	return hTmp;
	268	}
	269