Major update of the HFE package (comments inside the code
[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 // Class AliHFEtrdPIDqaV1
17 // Monitoring TRD PID in the HFE PID montioring framework. The following
18 // quantities are monitored:
19 //   TRD electron likelihood
20 //   TRD dE/dx (Absolute values)
21 //   TPC dE/dx (Number of sigmas, control histogram)
22 // (Always as function of momentum, particle species and centrality 
23 // before and after cut)
24 // More information about the PID monitoring framework can be found in
25 // AliHFEpidQAmanager.cxx and AliHFEdetPIDqa.cxx
26 //
27 // Author:
28 //    Markus Fasel <M.Fasel@gsi.de>
29 //
30
31 #include <TAxis.h>
32 #include <TH2.h>
33 #include <THnSparse.h>
34 #include <TString.h>
35
36 #include "AliAODTrack.h"
37 #include "AliESDtrack.h"
38 #include "AliESDpid.h"
39 #include "AliExternalTrackParam.h"
40 #include "AliPID.h"
41
42 #include "AliHFEcollection.h"
43 #include "AliHFEpidBase.h"
44 #include "AliHFEtrdPIDqaV1.h"
45
46 ClassImp(AliHFEtrdPIDqaV1)
47
48 //____________________________________________________________
49 AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1():
50     AliHFEdetPIDqa(),
51     fHistos(NULL)
52 {
53   //
54   // Dummy constructor
55   //
56 }
57
58 //____________________________________________________________
59 AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1(const Char_t *name):
60     AliHFEdetPIDqa(name, "QA for TRD"),
61     fHistos(NULL)
62 {
63   //
64   // Default constructor
65   //
66 }
67
68 //____________________________________________________________
69 AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1(const AliHFEtrdPIDqaV1 &o):
70     AliHFEdetPIDqa(o),
71     fHistos(NULL)
72 {
73   //
74   // Copy constructor
75   //
76 }
77
78 //____________________________________________________________
79 AliHFEtrdPIDqaV1 &AliHFEtrdPIDqaV1::operator=(const AliHFEtrdPIDqaV1 &o){
80   //
81   // Make assignment
82   //
83   AliHFEdetPIDqa::operator=(o);
84   fHistos = o.fHistos;
85   
86   return *this;
87 }
88
89 //_________________________________________________________
90 Long64_t AliHFEtrdPIDqaV1::Merge(TCollection *coll){
91   //
92   // Merge with other objects
93   //
94   if(!coll) return 0;
95   if(coll->IsEmpty()) return 1;
96
97   TIter it(coll);
98   AliHFEtrdPIDqaV1 *refQA = NULL;
99   TObject *o = NULL;
100   Long64_t count = 0;
101   TList listHistos;
102   while((o = it())){
103     refQA = dynamic_cast<AliHFEtrdPIDqaV1 *>(o);
104     if(!refQA) continue;
105
106     listHistos.Add(refQA->fHistos);
107     count++; 
108   }
109   fHistos->Merge(&listHistos);
110   return count + 1;
111 }
112
113 //____________________________________________________________
114 void AliHFEtrdPIDqaV1::Initialize(){
115   //
116   // Initialize QA histos for TRD PID
117   //
118
119   AliDebug(1, "Initializing PID QA for TRD");
120   // Make common binning
121   const Int_t kPIDbins = AliPID::kSPECIES + 1;
122   const Int_t kPbins = 100;
123   const Int_t kSteps = 2;
124   const Double_t kMinPID = -1;
125   const Double_t kMinP = 0.;
126   const Double_t kMaxPID = (Double_t)AliPID::kSPECIES;
127   const Double_t kMaxP = 20.;
128
129   fHistos = new AliHFEcollection("trdqahistos", "Collection of TRD QA histograms");
130   
131   // Create Control Histogram monitoring the TPC sigma between the selection steps
132   const Int_t kTPCSigmaBins = 140;
133   Int_t nBinsTPCSigma[4] = {kPIDbins, kPbins, kTPCSigmaBins, kSteps};
134   Double_t minTPCSigma[4] = {kMinPID, kMinP, -12., 0};
135   Double_t maxTPCSigma[4] = {kMaxPID, kMaxP, 12., 2.};
136   fHistos->CreateTHnSparse("hTPCsigma", "TPC sigma; species p [GeV/c]; TPC dEdx - <dE/dx>|_{el} [#sigma]; selection step", 4, nBinsTPCSigma, minTPCSigma, maxTPCSigma);
137   // Create Monitoring histogram for the Likelihood distribution
138   const Int_t kTRDLikelihoodBins = 100;
139   Int_t nBinsTRDlike[4] = {kPIDbins, kPbins, kTRDLikelihoodBins, kSteps};
140   Double_t minTRDlike[4] = {kMinPID, kMinP, 0., 0.};
141   Double_t maxTRDlike[4] = {kMaxPID, kMaxP, 1., 2.};
142   fHistos->CreateTHnSparse("hTRDlikelihood", "TRD Likelihood Distribution; species; p [GeV/c]; TRD electron Likelihood; selection step", 4, nBinsTRDlike, minTRDlike, maxTRDlike);
143   // Create Monitoring histogram for the TRD total charge
144   const Int_t kTRDchargeBins = 1000;
145   Int_t nBinsTRDcharge[4] = {kPIDbins, kPbins, kTRDchargeBins, kSteps};
146   Double_t minTRDcharge[4] = {kMinPID, kMinP, 0., 0.};
147   Double_t maxTRDcharge[4] = {kMaxPID, kMaxP, 100000., 2.};
148   fHistos->CreateTHnSparse("hTRDcharge", "Total TRD charge; species; p [GeV/c]; TRD charge [a.u.]; selection step", 4, nBinsTRDcharge, minTRDcharge, maxTRDcharge);
149 }
150
151 //____________________________________________________________
152 void AliHFEtrdPIDqaV1::ProcessTrack(AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step){
153   //
154   // Process the track, fill the containers 
155   //
156   AliDebug(1, Form("QA started for TRD PID for step %d", (Int_t)step));
157   Int_t species = track->GetAbInitioPID();
158   const AliVParticle *rectrack = track->GetRecTrack();
159   if(species >= AliPID::kSPECIES) species = -1;
160   if(!TString(rectrack->IsA()->GetName()).CompareTo("AliESDtrack")) ProcessESDtrack(dynamic_cast<const AliESDtrack *>(rectrack), step, species);
161   else if(!TString(rectrack->IsA()->GetName()).CompareTo("AliAODTrack")) ProcessAODtrack(dynamic_cast<const AliAODTrack *>(rectrack), step, species);
162   else  AliWarning(Form("Object type %s not supported\n", rectrack->IsA()->GetName()));
163 }
164
165 //____________________________________________________________
166 void AliHFEtrdPIDqaV1::ProcessESDtrack(const AliESDtrack *track, AliHFEdetPIDqa::EStep_t step, Int_t species){
167   //
168   // Fill QA histograms
169   //
170   Double_t container[4];
171   container[0] = species;
172   container[1] = track->GetOuterParam() ? track->GetOuterParam()->P() : track->P();
173   container[2] = fESDpid->NumberOfSigmasTPC(track, AliPID::kElectron);
174   container[3] = step;
175   fHistos->Fill("hTPCsigma", container);
176
177   Double_t pid[5];
178   track->GetTRDpid(pid);
179   container[2] = pid[AliPID::kElectron];
180   fHistos->Fill("hTRDlikelihood", container);
181
182   Double_t charge = 0.;
183   for(Int_t ily = 0; ily < 6; ily++){
184     charge = 0.;
185     for(Int_t isl = 0; isl < track->GetNumberOfTRDslices(); isl++){
186       charge += track->GetTRDslice(ily, isl);
187     }
188     fHistos->Fill("hTRDcharge", container);
189   }
190 }
191
192 //_________________________________________________________
193 void AliHFEtrdPIDqaV1::ProcessAODtrack(const AliAODTrack * /*track*/, AliHFEdetPIDqa::EStep_t /*step*/, Int_t /*species*/){
194   //
195   // To be implemented
196   //
197 }
198
199 //_________________________________________________________
200 TH2 *AliHFEtrdPIDqaV1::MakeTPCspectrumNsigma(AliHFEdetPIDqa::EStep_t step, Int_t species){
201   //
202   // Get the TPC control histogram for the TRD selection step (either before or after PID)
203   //
204   printf("histos :%p\n", fHistos);
205   THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTPCsigma"));
206   if(!histo){
207     AliError("QA histogram monitoring TPC nSigma not available");
208     return NULL;
209   }
210   if(species > -1 && species < AliPID::kSPECIES){
211     // cut on species (if available)
212     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
213   }
214   histo->GetAxis(3)->SetRangeUser(species + 1, species + 1);
215
216   TH2 *hSpec = histo->Projection(2, 1);
217   // construct title and name
218   TString stepname = step ? "before" : "after";
219   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
220   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
221   TString histname = Form("hSigmaTPC%s%s", specID.Data(), stepname.Data());
222   TString histtitle = Form("TPC Sigma for %s %s PID", speciesname.Data(), stepname.Data());
223   hSpec->SetName(histname.Data());
224   hSpec->SetTitle(histtitle.Data());
225
226   // Unset range on the original histogram
227   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
228   histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
229   return hSpec; 
230 }
231
232 //_________________________________________________________
233 TH2 *AliHFEtrdPIDqaV1::MakeTRDlikelihoodDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species){
234   //
235   // Make Histogram for TRD Likelihood distribution
236   //
237   THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDlikelihood"));
238   if(!histo){
239     AliError("QA histogram monitoring TRD Electron Likelihood not available");
240     return NULL;
241   }
242   if(species > -1 && species < AliPID::kSPECIES){
243     // cut on species (if available)
244     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
245   }
246   histo->GetAxis(3)->SetRangeUser(species + 1, species + 1);
247
248   TH2 *hSpec = histo->Projection(2, 1);
249   // construct title and name
250   TString stepname = step ? "before" : "after";
251   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
252   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
253   TString histname = Form("hLikeElTRD%s%s", specID.Data(), stepname.Data());
254   TString histtitle = Form("TRD electron Likelihood for %s %s PID", speciesname.Data(), stepname.Data());
255   hSpec->SetName(histname.Data());
256   hSpec->SetTitle(histtitle.Data());
257
258   // Unset range on the original histogram
259   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
260   histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
261   return hSpec; 
262 }
263
264 //_________________________________________________________
265 TH2 *AliHFEtrdPIDqaV1::MakeTRDchargeDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species){
266   //
267   // Make Histogram for TRD Likelihood distribution
268   //
269   THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDcharge"));
270   if(!histo){
271     AliError("QA histogram monitoring TRD total charge not available");
272     return NULL;
273   }
274   if(species > -1 && species < AliPID::kSPECIES){
275     // cut on species (if available)
276     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
277   }
278   histo->GetAxis(3)->SetRangeUser(species + 1, species + 1);
279
280   TH2 *hSpec = histo->Projection(2, 1);
281   // construct title and name
282   TString stepname = step ? "before" : "after";
283   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
284   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
285   TString histname = Form("hChargeTRD%s%s", specID.Data(), stepname.Data());
286   TString histtitle = Form("TRD total charge for %s %s PID", speciesname.Data(), stepname.Data());
287   hSpec->SetName(histname.Data());
288   hSpec->SetTitle(histtitle.Data());
289
290   // Unset range on the original histogram
291   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
292   histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
293   return hSpec; 
294 }
295