Update of the HFE package
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEtpcPIDqa.cxx
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 <TBrowser.h>
30 #include <TClass.h>
31 #include <TH2.h>
32 #include <THnSparse.h>
33 #include <TString.h>
34
35 #include "AliAODTrack.h"
36 #include "AliESDtrack.h"
37 #include "AliESDpid.h"
38 #include "AliLog.h"
39 #include "AliPID.h"
40
41 #include "AliHFEcollection.h"
42 #include "AliHFEpidBase.h"
43 #include "AliHFEpidQAmanager.h"
44 #include "AliHFEpidTPC.h"
45 #include "AliHFEtools.h"
46 #include "AliHFEtpcPIDqa.h"
47
48 ClassImp(AliHFEtpcPIDqa)
49
50 //_________________________________________________________
51 AliHFEtpcPIDqa::AliHFEtpcPIDqa():
52     AliHFEdetPIDqa()
53   , fHistos(NULL)
54   , fBrowseCentrality(-1)
55 {
56   //
57   // Dummy constructor
58   //
59 }
60
61 //_________________________________________________________
62 AliHFEtpcPIDqa::AliHFEtpcPIDqa(const char* name):
63     AliHFEdetPIDqa(name, "QA for TPC")
64   , fHistos(NULL)
65   , fBrowseCentrality(-1)
66 {
67   //
68   // Default constructor
69   //
70 }
71
72 //_________________________________________________________
73 AliHFEtpcPIDqa::AliHFEtpcPIDqa(const AliHFEtpcPIDqa &o):
74     AliHFEdetPIDqa(o)
75   , fHistos(NULL)
76   , fBrowseCentrality(o.fBrowseCentrality)
77 {
78   //
79   // Copy constructor
80   //
81   o.Copy(*this);
82 }
83
84 //_________________________________________________________
85 AliHFEtpcPIDqa &AliHFEtpcPIDqa::operator=(const AliHFEtpcPIDqa &o){
86   //
87   // Do assignment
88   //
89   AliHFEdetPIDqa::operator=(o);
90   if(&o != this) o.Copy(*this);
91   return *this;
92 }
93
94 //_________________________________________________________
95 AliHFEtpcPIDqa::~AliHFEtpcPIDqa(){
96   //
97   // Destructor
98   //
99   if(fHistos) delete fHistos;
100 }
101
102 //_________________________________________________________
103 void AliHFEtpcPIDqa::Copy(TObject &o) const {
104   //
105   // Make copy
106   //
107   AliHFEtpcPIDqa &target = dynamic_cast<AliHFEtpcPIDqa &>(o);
108   if(target.fHistos){
109     delete target.fHistos;
110     target.fHistos = NULL;
111   }
112   if(fHistos) target.fHistos = new AliHFEcollection(*fHistos);
113   target.fBrowseCentrality = fBrowseCentrality;
114 }
115
116 //_________________________________________________________
117 Long64_t AliHFEtpcPIDqa::Merge(TCollection *coll){
118   //
119   // Merge with other objects
120   //
121   if(!coll) return 0;
122   if(coll->IsEmpty()) return 1;
123
124   TIter it(coll);
125   AliHFEtpcPIDqa *refQA = NULL;
126   TObject *o = NULL;
127   Long64_t count = 0;
128   TList listHistos;
129   while((o = it())){
130     refQA = dynamic_cast<AliHFEtpcPIDqa *>(o);
131     if(!refQA) continue;
132
133     listHistos.Add(refQA->fHistos);
134     count++; 
135   }
136   fHistos->Merge(&listHistos);
137   return count + 1;
138 }
139
140 //_________________________________________________________
141 void AliHFEtpcPIDqa::Browse(TBrowser *b){
142   //
143   // Browse the PID QA
144   //
145   if(b){
146     if(fHistos){
147       b->Add(fHistos, fHistos->GetName());
148
149       // Make Projections of the dE/dx Spectra and add them to a new Folder
150       TString specnames[AliPID::kSPECIES+1] = {"All", "Electrons", "Muon", "Pions", "Kaon", "Protons"};
151       Int_t specind[AliPID::kSPECIES+1] = {-1, AliPID::kElectron, AliPID::kMuon, AliPID::kPion, AliPID::kKaon, AliPID::kProton};
152       TList *listdEdx = new TList;
153       listdEdx->SetOwner();
154       TList *listNsigma = new TList;
155       listNsigma->SetOwner();
156
157       TH2 *hptr = NULL; 
158       for(Int_t ispec = 0; ispec < AliPID::kSPECIES+1; ispec++){
159         for(Int_t istep = 0; istep < 2; istep++){
160           hptr = MakeSpectrumdEdx(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], fBrowseCentrality);
161           hptr->SetName(Form("hTPCdEdx%s%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After" , fBrowseCentrality == -1 ? "MinBias" : Form("Cent%d", fBrowseCentrality)));
162           listdEdx->Add(hptr);
163           hptr = MakeSpectrumNSigma(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], fBrowseCentrality);
164           hptr->SetName(Form("hTPCnsigma%s%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After", fBrowseCentrality == -1 ? "MinBias" : Form("Cent%d", fBrowseCentrality)));
165           listNsigma->Add(hptr);
166         }
167       }
168       
169       b->Add(listdEdx, "Projections dE/dx");
170       b->Add(listNsigma, "Projections NSigma");
171     }
172   }
173 }
174
175 //_________________________________________________________
176 void AliHFEtpcPIDqa::Initialize(){
177   //
178   // Define Histograms
179   //
180
181   fHistos = new AliHFEcollection("tpcqahistos", "Collection of TPC QA histograms");
182
183   // Make common binning
184   const Int_t kNdim = 5;
185   const Int_t kPIDbins = AliPID::kSPECIES + 1;
186   const Int_t kSteps = 2;
187   const Int_t kCentralityBins = 11;
188   const Double_t kMinPID = -1;
189   const Double_t kMinP = 0.;
190   const Double_t kMaxPID = (Double_t)AliPID::kSPECIES;
191   const Double_t kMaxP = 20.;
192   // Quantities where one can switch between low and high resolution
193   Int_t kPbins = fQAmanager->HasHighResolutionHistos() ? 1000 : 100;
194   Int_t kDedxbins = fQAmanager->HasHighResolutionHistos() ? 400 : 200;
195   Int_t kSigmaBins = fQAmanager->HasHighResolutionHistos() ? 1400 : 240;
196  
197   // 1st histogram: TPC dEdx: (species, p, dEdx, step)
198   Int_t nBinsdEdx[kNdim] = {kPIDbins, kPbins, kDedxbins, kSteps, kCentralityBins};
199   Double_t mindEdx[kNdim] =  {kMinPID, kMinP, 0., 0., 0.};
200   Double_t maxdEdx[kNdim] =  {kMaxPID, kMaxP, 200, 2., 11.}; 
201   fHistos->CreateTHnSparse("tpcDedx", "TPC signal; species; p [GeV/c]; TPC signal [a.u.]; Selection Step; Centrality", kNdim, nBinsdEdx, mindEdx, maxdEdx);
202   // 2nd histogram: TPC sigmas: (species, p nsigma, step)
203   Int_t nBinsSigma[kNdim] = {kPIDbins, kPbins, kSigmaBins, kSteps, kCentralityBins};
204   Double_t minSigma[kNdim] = {kMinPID, kMinP, -12., 0., 0.};
205   Double_t maxSigma[kNdim] = {kMaxPID, kMaxP, 12., 2., 11.};
206   fHistos->CreateTHnSparse("tpcnSigma", "TPC signal; species; p [GeV/c]; TPC signal [a.u.]; Selection Step; Centrality", kNdim, nBinsSigma, minSigma, maxSigma);
207
208   // General TPC QA
209 }
210
211 //_________________________________________________________
212 void AliHFEtpcPIDqa::ProcessTrack(const AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step){
213   //
214   // Fill TPC histograms
215   //
216   AliDebug(1, Form("QA started for TPC PID for step %d", (Int_t)step));
217   AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
218   Float_t centrality = track->GetCentrality();
219   Int_t species = track->GetAbInitioPID();
220   if(species >= AliPID::kSPECIES) species = -1;
221
222   AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fQAmanager->GetDetectorPID(AliHFEpid::kTPCpid));
223   Double_t contentSignal[5];
224   contentSignal[0] = species;
225   contentSignal[1] = tpcpid ? tpcpid->GetP(track->GetRecTrack(), anatype) : 0.;
226   contentSignal[2] = GetTPCsignal(track->GetRecTrack(), anatype);
227   contentSignal[3] = step;
228   contentSignal[4] = centrality;
229   fHistos->Fill("tpcDedx", contentSignal);
230
231   contentSignal[2] = tpcpid ? tpcpid->NumberOfSigmas(track->GetRecTrack(), AliPID::kElectron, anatype) : 0.; 
232   fHistos->Fill("tpcnSigma", contentSignal);
233 }
234
235 //_________________________________________________________
236 Double_t AliHFEtpcPIDqa::GetTPCsignal(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anatype){
237   //
238   // Get TPC signal for ESD and AOD track
239   //
240   Double_t tpcSignal = 0.;
241   if(anatype == AliHFEpidObject::kESDanalysis){
242     const AliESDtrack *esdtrack = static_cast<const AliESDtrack *>(track);
243     tpcSignal = esdtrack->GetTPCsignal();
244   } else {
245     const AliAODTrack *aodtrack = static_cast<const AliAODTrack *>(track);
246     tpcSignal = aodtrack->GetDetPid() ? aodtrack->GetDetPid()->GetTPCsignal() : 0.;
247   }
248   return tpcSignal;
249 }
250
251 //_________________________________________________________
252 TH2 *AliHFEtpcPIDqa::MakeSpectrumdEdx(AliHFEdetPIDqa::EStep_t istep, Int_t species, Int_t centralityClass){
253   //
254   // Plot the Spectrum
255   //
256   THnSparseF *hSignal = dynamic_cast<THnSparseF *>(fHistos->Get("tpcDedx"));
257   if(!hSignal) return NULL;
258   hSignal->GetAxis(3)->SetRange(istep + 1, istep + 1);
259   if(species >= 0 && species < AliPID::kSPECIES)
260     hSignal->GetAxis(0)->SetRange(2 + species, 2 + species);
261   TString hname = Form("hTPCsignal%s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after"), 
262           htitle = Form("TPC dE/dx Spectrum %s selection", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
263   if(species > -1){
264     hname += AliPID::ParticleName(species);
265     htitle += Form(" for %ss", AliPID::ParticleName(species));
266   }
267   TString centname, centtitle;
268   Bool_t hasCentralityInfo = kTRUE;
269   if(centralityClass > -1){
270     if(hSignal->GetNdimensions() < 5){
271       AliError("Centrality Information not available");
272       centname = centtitle = "MinBias";
273       hasCentralityInfo= kFALSE;
274     } else {
275       // Project centrality classes
276       // -1 is Min. Bias
277       hSignal->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
278       centname = Form("Cent%d", centralityClass);
279       centtitle = Form(" Centrality %d", centralityClass);
280     }
281   } else {
282     centname = centtitle = "MinBias";
283     hasCentralityInfo= kFALSE;
284   }
285   hname += centtitle;
286   htitle += centtitle;
287
288   TH2 *hTmp = hSignal->Projection(2,1);
289   hTmp->SetName(hname.Data());
290   hTmp->SetTitle(htitle.Data());
291   hTmp->SetStats(kFALSE);
292   hTmp->GetXaxis()->SetTitle("p [GeV/c]");
293   hTmp->GetYaxis()->SetTitle("TPC signal [a.u.]");
294   hSignal->GetAxis(3)->SetRange(0, hSignal->GetAxis(3)->GetNbins());
295   hSignal->GetAxis(0)->SetRange(0, hSignal->GetAxis(0)->GetNbins());
296   if(hasCentralityInfo) hSignal->GetAxis(4)->SetRange(0, hSignal->GetAxis(4)->GetNbins());
297   return hTmp;
298 }
299
300 //_________________________________________________________
301 TH2 *AliHFEtpcPIDqa::MakeSpectrumNSigma(AliHFEdetPIDqa::EStep_t istep, Int_t species, Int_t centralityClass){
302   //
303   // Plot the Spectrum
304   //
305   THnSparseF *hSignal = dynamic_cast<THnSparseF *>(fHistos->Get("tpcnSigma"));
306   if(!hSignal) return NULL;
307   hSignal->GetAxis(3)->SetRange(istep + 1, istep + 1);
308   if(species >= 0 && species < AliPID::kSPECIES)
309     hSignal->GetAxis(0)->SetRange(2 + species, 2 + species);
310   TString hname = Form("hTPCsigma%s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after"), 
311           htitle = Form("TPC dE/dx Spectrum[#sigma] %s selection", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
312   if(species > -1){
313     hname += AliPID::ParticleName(species);
314     htitle += Form(" for %ss", AliPID::ParticleName(species));
315   }
316   TString centname, centtitle;
317   Bool_t hasCentralityInfo = kTRUE;
318   if(centralityClass > -1){
319     if(hSignal->GetNdimensions() < 5){
320       AliError("Centrality Information not available");
321       centname = centtitle = "MinBias";
322       hasCentralityInfo= kFALSE;
323     } else {
324       // Project centrality classes
325       // -1 is Min. Bias
326       hSignal->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
327       centname = Form("Cent%d", centralityClass);
328       centtitle = Form(" Centrality %d", centralityClass);
329     }
330   } else {
331     centname = centtitle = "MinBias";
332     hasCentralityInfo= kFALSE;
333   }
334   hname += centtitle;
335   htitle += centtitle;
336
337   TH2 *hTmp = hSignal->Projection(2,1);
338   hTmp->SetName(hname.Data());
339   hTmp->SetTitle(htitle.Data());
340   hTmp->SetStats(kFALSE);
341   hTmp->GetXaxis()->SetTitle("p [GeV/c]");
342   hTmp->GetYaxis()->SetTitle("TPC dE/dx - <dE/dx>|_{el} [#sigma]");
343   hSignal->GetAxis(3)->SetRange(0, hSignal->GetAxis(3)->GetNbins());
344   hSignal->GetAxis(0)->SetRange(0, hSignal->GetAxis(0)->GetNbins());
345   if(hasCentralityInfo) hSignal->GetAxis(4)->SetRange(0, hSignal->GetAxis(4)->GetNbins());
346   return hTmp;
347 }
348