Cleanup the code. Fix memory leak. Now inherit from AliAnalysisTaskSE (Antoine, Phili...
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEtrdPIDqaV1.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 /* $Id$ */
17
18 //
19 // Class AliHFEtrdPIDqaV1
20 // Monitoring TRD PID in the HFE PID montioring framework. The following
21 // quantities are monitored:
22 //   TRD electron likelihood
23 //   TRD dE/dx (Absolute values)
24 //   TPC dE/dx (Number of sigmas, control histogram)
25 // (Always as function of momentum, particle species and centrality 
26 // before and after cut)
27 // More information about the PID monitoring framework can be found in
28 // AliHFEpidQAmanager.cxx and AliHFEdetPIDqa.cxx
29 //
30 // Author:
31 //    Markus Fasel <M.Fasel@gsi.de>
32 //
33
34 #include <TAxis.h>
35 #include <TBrowser.h>
36 #include <TH2.h>
37 #include <THnSparse.h>
38 #include <TString.h>
39
40 #include "AliESDtrack.h"
41 #include "AliLog.h"
42 #include "AliPID.h"
43
44 #include "AliHFEcollection.h"
45 #include "AliHFEpidBase.h"
46 #include "AliHFEpidQAmanager.h"
47 #include "AliHFEpidTPC.h"
48 #include "AliHFEpidTRD.h"
49 #include "AliHFEtrdPIDqaV1.h"
50
51 ClassImp(AliHFEtrdPIDqaV1)
52
53 //____________________________________________________________
54 AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1():
55     AliHFEdetPIDqa(),
56     fHistos(NULL)
57 {
58   //
59   // Dummy constructor
60   //
61 }
62
63 //____________________________________________________________
64 AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1(const Char_t *name):
65     AliHFEdetPIDqa(name, "QA for TRD"),
66     fHistos(NULL)
67 {
68   //
69   // Default constructor
70   //
71 }
72
73 //____________________________________________________________
74 AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1(const AliHFEtrdPIDqaV1 &o):
75     AliHFEdetPIDqa(o),
76     fHistos(NULL)
77 {
78   //
79   // Copy constructor
80   //
81 }
82
83 //____________________________________________________________
84 AliHFEtrdPIDqaV1 &AliHFEtrdPIDqaV1::operator=(const AliHFEtrdPIDqaV1 &o){
85   //
86   // Make assignment
87   //
88   AliHFEdetPIDqa::operator=(o);
89   fHistos = o.fHistos;
90   
91   return *this;
92 }
93
94 //_________________________________________________________
95 Long64_t AliHFEtrdPIDqaV1::Merge(TCollection *coll){
96   //
97   // Merge with other objects
98   //
99   if(!coll) return 0;
100   if(coll->IsEmpty()) return 1;
101
102   TIter it(coll);
103   AliHFEtrdPIDqaV1 *refQA = NULL;
104   TObject *o = NULL;
105   Long64_t count = 0;
106   TList listHistos;
107   while((o = it())){
108     refQA = dynamic_cast<AliHFEtrdPIDqaV1 *>(o);
109     if(!refQA) continue;
110
111     listHistos.Add(refQA->fHistos);
112     count++; 
113   }
114   fHistos->Merge(&listHistos);
115   return count + 1;
116 }
117
118 //_________________________________________________________
119 void AliHFEtrdPIDqaV1::Browse(TBrowser *b){
120   //
121   // Browse the PID QA
122   //
123   if(b){
124     if(fHistos){
125       b->Add(fHistos, fHistos->GetName());
126
127       // Make Projections of the dE/dx Spectra and add them to a new Folder
128       TString specnames[4] = {"All", "Electrons", "Pions", "Protons"};
129       Int_t specind[4] = {-1, AliPID::kElectron, AliPID::kPion, AliPID::kProton};
130       TList *listTM = new TList;
131       listTM->SetOwner();
132       TList *listLike = new TList;
133       listLike->SetOwner();
134       TList *listCharge = new TList;
135       listCharge->SetOwner();
136       TList *listTPCnsigma = new TList;
137       listTPCnsigma->SetOwner();
138
139       TH2 *hptr = NULL; 
140       for(Int_t ispec = 0; ispec < 4; ispec++){
141         for(Int_t istep = 0; istep < 2; istep++){
142           hptr = MakeTRDspectrumTM(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
143           if(hptr){
144             hptr->SetName(Form("hTRDtm%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
145             listTM->Add(hptr);
146           }
147           hptr = MakeTRDlikelihoodDistribution(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
148           hptr->SetName(Form("hTRDlike%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
149           listLike->Add(hptr);
150           hptr = MakeTRDchargeDistribution(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
151           hptr->SetName(Form("hTRDcharge%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
152           listCharge->Add(hptr);
153           hptr = MakeTPCspectrumNsigma(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
154           hptr->SetName(Form("hTPCspectrum%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
155           listTPCnsigma->Add(hptr);
156         }
157       }
158       
159       b->Add(listTM, "Projections Truncated Mean");
160       b->Add(listLike, "Projections Likelihood distribution");
161       b->Add(listCharge, "Projections Tracklet Charge");
162       b->Add(listTPCnsigma, "Projections TPC spectra");
163     }
164   }
165 }
166
167 //____________________________________________________________
168 void AliHFEtrdPIDqaV1::Initialize(){
169   //
170   // Initialize QA histos for TRD PID
171   //
172
173   AliDebug(1, "Initializing PID QA for TRD");
174   // Make common binning
175   const Int_t kPIDbins = AliPID::kSPECIES + 1;
176   const Int_t kPbins = 100;
177   const Int_t kSteps = 2;
178   const Double_t kMinPID = -1;
179   const Double_t kMinP = 0.;
180   const Double_t kMaxPID = (Double_t)AliPID::kSPECIES;
181   const Double_t kMaxP = 20.;
182
183   fHistos = new AliHFEcollection("trdqahistos", "Collection of TRD QA histograms");
184   
185   // Create Control Histogram monitoring the TPC sigma between the selection steps
186   const Int_t kTPCSigmaBins = 140;
187   Int_t nBinsTPCSigma[4] = {kPIDbins, kPbins, kTPCSigmaBins, kSteps};
188   Double_t minTPCSigma[4] = {kMinPID, kMinP, -12., 0};
189   Double_t maxTPCSigma[4] = {kMaxPID, kMaxP, 12., 2.};
190   fHistos->CreateTHnSparse("hTPCsigma", "TPC sigma; species p [GeV/c]; TPC dEdx - <dE/dx>|_{el} [#sigma]; selection step", 4, nBinsTPCSigma, minTPCSigma, maxTPCSigma);
191   fHistos->Sumw2("hTPCsigma");
192   // Create Monitoring histogram for the Likelihood distribution
193   const Int_t kTRDLikelihoodBins = 100;
194   Int_t nBinsTRDlike[4] = {kPIDbins, kPbins, kTRDLikelihoodBins, kSteps};
195   Double_t minTRDlike[4] = {kMinPID, kMinP, 0., 0.};
196   Double_t maxTRDlike[4] = {kMaxPID, kMaxP, 1., 2.};
197   fHistos->CreateTHnSparse("hTRDlikelihood", "TRD Likelihood Distribution; species; p [GeV/c]; TRD electron Likelihood; selection step", 4, nBinsTRDlike, minTRDlike, maxTRDlike);
198   fHistos->Sumw2("hTRDlikelihood");
199   // Create Monitoring histogram for the TRD total charge
200   const Int_t kTRDchargeBins = 1000;
201   Int_t nBinsTRDcharge[4] = {kPIDbins, kPbins, kTRDchargeBins, kSteps};
202   Double_t minTRDcharge[4] = {kMinPID, kMinP, 0., 0.};
203   Double_t maxTRDcharge[4] = {kMaxPID, kMaxP, 100000., 2.};
204   fHistos->CreateTHnSparse("hTRDcharge", "Total TRD charge; species; p [GeV/c]; TRD charge [a.u.]; selection step", 4, nBinsTRDcharge, minTRDcharge, maxTRDcharge);
205   fHistos->Sumw2("hTRDcharge");
206   // Monitoring of the TRD truncated mean according to version 1
207   const Int_t kTRDtmBins = 1000;
208   Int_t nBinsTRDtm[4] = {kPIDbins, kPbins, kTRDtmBins, kSteps};
209   Double_t minTRDtm[4] = {kMinPID, kMinP, 0., 0.};
210   Double_t maxTRDtm[4] = {kMaxPID, kMaxP, 20000., 2.};
211   fHistos->CreateTHnSparse("hTRDtruncatedMean", "TRD truncated Mean; species; p [GeV/c]; TRD signal [a.u.]; selection step", 4, nBinsTRDtm, minTRDtm, maxTRDtm);
212   fHistos->Sumw2("hTRDtruncatedMean");
213 }
214
215 //____________________________________________________________
216 void AliHFEtrdPIDqaV1::ProcessTrack(const AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step){
217   //
218   // Process the track, fill the containers 
219   //
220   AliDebug(1, Form("QA started for TRD PID for step %d", (Int_t)step));
221   Int_t species = track->GetAbInitioPID();
222   if(species >= AliPID::kSPECIES) species = -1;
223   AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
224
225   AliHFEpidTRD *trdpid = dynamic_cast<AliHFEpidTRD *>(fQAmanager->GetDetectorPID(AliHFEpid::kTRDpid));
226   AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fQAmanager->GetDetectorPID(AliHFEpid::kTPCpid));
227  
228   Double_t container[4];
229   container[0] = species;
230   container[1] = trdpid ? trdpid->GetP(track->GetRecTrack(), anatype) : 0.;
231   container[2] = tpcpid ? tpcpid->NumberOfSigmas(track->GetRecTrack(), AliPID::kElectron, anatype) : 0.;
232   container[3] = step;
233   fHistos->Fill("hTPCsigma", container);
234
235   container[2] = trdpid->GetElectronLikelihood(track->GetRecTrack(), anatype);
236   fHistos->Fill("hTRDlikelihood", container);
237
238   if(track->IsESDanalysis()){
239     container[2] = trdpid->GetTRDSignalV1(dynamic_cast<const AliESDtrack *>(track->GetRecTrack()));
240     fHistos->Fill("hTRDtruncatedMean", container);
241   }
242   for(UInt_t ily = 0; ily < 6; ily++){
243     container[2] = trdpid->GetChargeLayer(track->GetRecTrack(), ily, anatype);
244     fHistos->Fill("hTRDcharge", container);
245   }
246 }
247
248 //_________________________________________________________
249 TH2 *AliHFEtrdPIDqaV1::MakeTPCspectrumNsigma(AliHFEdetPIDqa::EStep_t step, Int_t species){
250   //
251   // Get the TPC control histogram for the TRD selection step (either before or after PID)
252   //
253   THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTPCsigma"));
254   if(!histo){
255     AliError("QA histogram monitoring TPC nSigma not available");
256     return NULL;
257   }
258   if(species > -1 && species < AliPID::kSPECIES){
259     // cut on species (if available)
260     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
261   }
262   histo->GetAxis(3)->SetRange(step + 1, step + 1); 
263
264   TH2 *hSpec = histo->Projection(2, 1);
265   // construct title and name
266   TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
267   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
268   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
269   TString histname = Form("hSigmaTPC%s%s", specID.Data(), stepname.Data());
270   TString histtitle = Form("TPC Sigma for %s %s PID", speciesname.Data(), stepname.Data());
271   hSpec->SetName(histname.Data());
272   hSpec->SetTitle(histtitle.Data());
273
274   // Unset range on the original histogram
275   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
276   histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
277   return hSpec; 
278 }
279
280 //_________________________________________________________
281 TH2 *AliHFEtrdPIDqaV1::MakeTRDspectrumTM(AliHFEdetPIDqa::EStep_t step, Int_t species){
282   //
283   // Get the TPC control histogram for the TRD selection step (either before or after PID)
284   //
285   THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDtruncatedMean"));
286   if(!histo){
287     AliError("QA histogram monitoring TPC nSigma not available");
288     return NULL;
289   }
290   if(species > -1 && species < AliPID::kSPECIES){
291     // cut on species (if available)
292     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
293   }
294   histo->GetAxis(3)->SetRange(step + 1, step + 1); 
295
296   TH2 *hSpec = histo->Projection(2, 1);
297   // construct title and name
298   TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
299   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
300   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
301   TString histname = Form("hTMTRD%s%s", specID.Data(), stepname.Data());
302   TString histtitle = Form("TRD Truncated Mean for %s %s PID", speciesname.Data(), stepname.Data());
303   hSpec->SetName(histname.Data());
304   hSpec->SetTitle(histtitle.Data());
305
306   // Unset range on the original histogram
307   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
308   histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
309   return hSpec; 
310 }
311
312 //_________________________________________________________
313 TH2 *AliHFEtrdPIDqaV1::MakeTRDlikelihoodDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species){
314   //
315   // Make Histogram for TRD Likelihood distribution
316   //
317   THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDlikelihood"));
318   if(!histo){
319     AliError("QA histogram monitoring TRD Electron Likelihood not available");
320     return NULL;
321   }
322   if(species > -1 && species < AliPID::kSPECIES){
323     // cut on species (if available)
324     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
325   }
326   histo->GetAxis(3)->SetRangeUser(step + 1, step + 1);
327
328   TH2 *hSpec = histo->Projection(2, 1);
329   // construct title and name
330   TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
331   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
332   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
333   TString histname = Form("hLikeElTRD%s%s", specID.Data(), stepname.Data());
334   TString histtitle = Form("TRD electron Likelihood for %s %s PID", speciesname.Data(), stepname.Data());
335   hSpec->SetName(histname.Data());
336   hSpec->SetTitle(histtitle.Data());
337
338   // Unset range on the original histogram
339   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
340   histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
341   return hSpec; 
342 }
343
344 //_________________________________________________________
345 TH2 *AliHFEtrdPIDqaV1::MakeTRDchargeDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species){
346   //
347   // Make Histogram for TRD Likelihood distribution
348   //
349   THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDcharge"));
350   if(!histo){
351     AliError("QA histogram monitoring TRD total charge not available");
352     return NULL;
353   }
354   if(species > -1 && species < AliPID::kSPECIES){
355     // cut on species (if available)
356     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
357   }
358   histo->GetAxis(3)->SetRange(step + 1, step + 1);
359
360   TH2 *hSpec = histo->Projection(2, 1);
361   // construct title and name
362   TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
363   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
364   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
365   TString histname = Form("hChargeTRD%s%s", specID.Data(), stepname.Data());
366   TString histtitle = Form("TRD total charge for %s %s PID", speciesname.Data(), stepname.Data());
367   hSpec->SetName(histname.Data());
368   hSpec->SetTitle(histtitle.Data());
369
370   // Unset range on the original histogram
371   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
372   histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
373   return hSpec; 
374 }
375