]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/hfe/AliHFEpidQA.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Add fast merging option (Diego)
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEpidQA.cxx
CommitLineData
70da6c5a 1/**************************************************************************
faee3b18 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 **************************************************************************/
70da6c5a 15//
16// Class for PID QA
17// Several studies done on clean samples of electrons, pions and kaons
18// coming from V0 PID
19// Compatible with both ESDs and AODs
20//
21// Autors:
22// Matus Kalisky <matus.kalisky@cern.ch>
23// Markus Heide <mheide@uni-muenster.de>
24// Markus Fasel <M.Fasel@gsi.de>
25//
faee3b18 26
27
70da6c5a 28#include <TMath.h>
29#include <TObjArray.h>
70da6c5a 30#include <TPDGCode.h>
31#include <TString.h>
faee3b18 32#include <TMultiLayerPerceptron.h>
33#include <TFile.h>
70da6c5a 34
35#include "AliAODMCParticle.h"
faee3b18 36#include "AliAODEvent.h"
70da6c5a 37#include "AliAODTrack.h"
38#include "AliESDtrack.h"
faee3b18 39#include "AliESDEvent.h"
70da6c5a 40#include "AliMCEvent.h"
41#include "AliMCParticle.h"
42#include "AliPID.h"
8c1c76e9 43#include "AliPIDResponse.h"
70da6c5a 44#include "AliVParticle.h"
70da6c5a 45
46
47#include "AliHFEcollection.h"
48#include "AliHFEpidQA.h"
faee3b18 49#include "AliHFEV0info.h"
70da6c5a 50#include "AliHFEV0pid.h"
51#include "AliHFEV0pidMC.h"
3a72645a 52#include "AliHFEV0cuts.h"
faee3b18 53#include "AliHFEtrdPIDqa.h"
54
70da6c5a 55
56ClassImp(AliHFEpidQA)
57
faee3b18 58 //__________________________________________
59 AliHFEpidQA::AliHFEpidQA():
60 fEvent(NULL),
61 fMC(NULL),
62 fV0pid(NULL),
63 fV0pidMC(NULL),
64 fTRDpidQA(NULL),
65 fOutput(NULL),
66 fESDpid(NULL),
e156c3bb 67 fNNref(NULL),
68 fTRDTotalChargeInSlice0(kFALSE)
70da6c5a 69{
70 //
71 // Default constructor
72 //
faee3b18 73 for(Int_t mom = 0; mom < 11; mom++){
74 fNet[mom] = NULL;
75 }
76}
77
78//__________________________________________
79AliHFEpidQA::AliHFEpidQA(const AliHFEpidQA &ref):
80 TObject(ref),
81 fEvent(NULL),
82 fMC(NULL),
83 fV0pid(NULL),
84 fV0pidMC(NULL),
85 fTRDpidQA(NULL),
86 fOutput(NULL),
87 fESDpid(NULL),
e156c3bb 88 fNNref(NULL),
89 fTRDTotalChargeInSlice0(ref.fTRDTotalChargeInSlice0)
faee3b18 90{
91 //
92 // Copy constructor
93 //
94 for(Int_t mom = 0; mom < 11; mom++){
95 fNet[mom] = NULL;
96 }
97 ref.Copy(*this);
98}
99
100//__________________________________________
101AliHFEpidQA &AliHFEpidQA::operator=(const AliHFEpidQA &ref){
102 //
103 // Assignment operator
104 //
105 if(this != &ref)
106 ref.Copy(*this);
107 return *this;
70da6c5a 108}
109
110//__________________________________________
111AliHFEpidQA::~AliHFEpidQA(){
112 //
113 // Destructor
114 //
3e03bd67 115
116 // the pointers bellow are not dfeleted to prevend double-deleting of some of the content
117 // these pointers are defined only once during the program call and should not cause a problem,
118 // but cleaner solution is necessary.
119 //if(fV0pid) delete fV0pid;
120 //if(fV0pidMC) delete fV0pidMC;
121 //if(fOutput) delete fOutput;
122
faee3b18 123 // if(fTRDpidResponse) delete fTRDpidResponse;
124}
125
126//__________________________________________
127void AliHFEpidQA::Copy(TObject &o) const {
128 //
129 // Copy function
130 //
131
132 TObject::Copy(o);
133
134 AliHFEpidQA &target = dynamic_cast<AliHFEpidQA &>(o);
135 target.fMC = fMC;
e156c3bb 136 target.fTRDTotalChargeInSlice0 = fTRDTotalChargeInSlice0;
faee3b18 137
138 if(target.fESDpid) delete target.fESDpid;
8c1c76e9 139 target.fESDpid = fESDpid;
faee3b18 140 if(target.fV0pid) delete target.fV0pid;
141 if(fV0pid)
142 target.fV0pid = dynamic_cast<AliHFEV0pid *>(fV0pid->Clone());
143 else
144 target.fV0pid = NULL;
145 if(target.fV0pidMC) delete target.fV0pidMC;
146 if(fV0pidMC)
147 target.fV0pidMC = dynamic_cast<AliHFEV0pidMC *>(fV0pidMC->Clone());
148 else
149 target.fV0pidMC = NULL;
150 if(target.fTRDpidQA) delete target.fTRDpidQA;
151 if(fTRDpidQA)
152 target.fTRDpidQA = dynamic_cast<AliHFEtrdPIDqa *>(fTRDpidQA->Clone());
153 else
154 target.fTRDpidQA = NULL;
155 if(target.fOutput) delete target.fOutput;
156 if(fOutput)
157 target.fOutput = dynamic_cast<AliHFEcollection *>(fOutput->Clone());
70da6c5a 158}
159
160//__________________________________________
161void AliHFEpidQA::Init(){
162 //
163 // Prepare task output
164 //
165
3a72645a 166 // load networks
faee3b18 167 if(fNNref){
3a72645a 168 for(Int_t mom = 0; mom < 11; mom++){
faee3b18 169 fNet[mom] = (TMultiLayerPerceptron*) fNNref->Get(Form("NN_Mom%d", mom));
170 if(!fNet[mom]){
171 AliError(Form("No reference network for momentum bin %d!", mom));
172 }
173 }
174 }
175
63bffcf1 176 fV0pid = new AliHFEV0pid("fV0pid");
70da6c5a 177 if(HasV0pidQA()) fV0pid->InitQA();
178 fV0pidMC = new AliHFEV0pidMC();
179 fV0pidMC->Init();
180
faee3b18 181 fTRDpidQA = new AliHFEtrdPIDqa;
e156c3bb 182 if(fTRDTotalChargeInSlice0) fTRDpidQA->SetTotalChargeInSlice0();
faee3b18 183 fTRDpidQA->Init();
184
70da6c5a 185 fOutput = new AliHFEcollection("pidQA", "PID QA output");
186
faee3b18 187 const char *name[4] = {"Electron", "PionK0", "PionL", "Proton"};
188 const char *title[4] = {"Electron", "K0 Pion", "Lambda Pion", "Proton"};
189 const char *det[4] = {"ITS", "TPC", "TRD", "TOF"};
8c1c76e9 190 const int effs[6] = {70, 75, 80, 85, 90, 95};
faee3b18 191 for(Int_t i = 0; i < 4; i++){
192 fOutput->CreateTH2F(Form("purity%s", name[i]), Form("%s Purity", title[i]), 2, -0.5, 1.5, 20, 0.1, 10, 1);
193
194 for(Int_t idet = 0; idet < 4; idet++){
195 // create all the histograms which all the detectors have in common
196 if(idet != 2){ // no nSigma histogram for TRD
197 fOutput->CreateTH2F(Form("h%s_nSigma_%s", det[idet], name[i]), Form("%s number of sigmas for %ss; p (GeV/c); number of sigmas", det[idet], title[i]), 20, 0.1, 10, 100, -7, 7, 0);
198 }
199 fOutput->CreateTH2F(Form("h%s_PID_p_%s", det[idet], name[i]), Form("%s PID for %ss; p (GeV/c); ITS PID", det[idet], title[i]), 100, 0.1, 10, 5, -0.5, 4.5, 0);
200 fOutput->CreateTH2F(Form("h%s_El_like_%s", det[idet], name[i]), Form("%s Electron Likelihoods for %ss; p (GeV/c); likelihood", det[idet], title[i]), 25, 0.1, 10, 1000, 0., 1., 0);
201 fOutput->CreateTH2F(Form("h%s_El_like_MC_%s", det[idet], name[i]), Form("%s Electron Likelihoods for MC %ss; p (GeV/c); likelihood", det[idet], title[i]), 25, 0.1, 10, 1000, 0., 1., 0);
202 }
70da6c5a 203
faee3b18 204 //
205 // ITS pid response
206 //
207 fOutput->CreateTH2F(Form("hITS_Signal_%s", name[i]), Form("ITS Signal org. for %ss", title[i]), 40, 0.1, 10, 400, 0, 1000, 0);
208 fOutput->CreateTH2Fvector1(5, Form("hITS_dEdx_%s", name[i]), Form("ITS dEdx for %ss; p (GeV/c); dEdx (a.u.)", title[i]), 40, 0.1, 10, 400, 0, 1000, 0);
209
210 //
211 // TPC pid response
212 //
213 fOutput->CreateTH2F(Form("hTPC_dEdx_%s", name[i]), Form("TPC dEdx for %ss; p (GeV/c); dEdx (a.u.)", title[i]), 20, 0.1, 10, 200, 0, 200, 0);
214
215 //
216 // TRD pid response
217 //
218 fOutput->CreateTH2F(Form("hTRD_trk_%s", name[i]), Form("%s PID tracklets; p (GeV/c); N TRD tracklets", title[i]), 100, 0.1, 10, 7, -0.5, 6.5, 0);
219 // number of the non 0 slices per tracklet
220 fOutput->CreateTH2F(Form("hTRD_Nslc_%s", name[i]), Form("%s PID slices > 0; p (GeV/c); N slc > 0", title[i]), 100, 0.1, 10, 9, -0.5, 8.5, 0);
221 // location of the slices > 0 - where are the emtpy slices located ?
222 fOutput->CreateTH2F(Form("hTRD_slc_%s", name[i]), Form("%s PID slices > 0 position; p (GeV/c); slice", title[i]), 100, 0.1, 10, 9, -0.5, 8.5, 0);
223 fOutput->CreateTH2F(Form("hTRD_cls_%s", name[i]), Form("TRD clusters for %s candidates; p (GeV/c); N cls", title[i]), 25, 0.1, 10, 1000, 0, 1000);
224 fOutput->CreateTH2F(Form("hTRD_dEdx_%s", name[i]), Form("TRD dEdx (trk) for %s candidates; p (GeV/c); tracklet dEdx (a.u.)", title[i]), 25, 0.1, 10, 1000, 0, 100000, 0);
8c1c76e9 225
226 for(Int_t itl = 4; itl < 6; itl++){
227 for(Int_t ieff = 0; ieff < 6; ieff++){
228 fOutput->CreateTH2F(Form("hTRD_likeSel_%s_%dtls_%deff", name[i], itl, effs[ieff]), Form(" %s selected as electrons with %d tracklets and %f electron efficiency", title[i], itl, static_cast<Double_t>(effs[ieff])/100.), 44, 0.1, 20., 100, 0., 1.);
229 fOutput->CreateTH2F(Form("hTRD_likeRej_%s_%dtls_%deff", name[i], itl, effs[ieff]), Form(" %s rejected as electrons with %d tracklets and %f electron efficiency", title[i], itl, static_cast<Double_t>(effs[ieff])/100.), 44, 0.1, 20., 100, 0., 1.);
230 }
231 }
faee3b18 232 //
233 // TOF pid response
234 //
235 fOutput->CreateTH2F(Form("hTOF_beta_%s", name[i]), Form("TOF beta for %s; p (GeV/c); #beta", title[i]), 50, 0.1, 5, 350, 0.4, 1.1, 0);
236 }//.. loop over identified particle species
237
238 // Global histograms
239 fOutput->CreateTH1F("hITS_dEdx_nSamples", "ITS - number of non 0 dEdx samples", 4, 0.5, 4.5);
70da6c5a 240 fOutput->CreateTH2F("hTPC_PID", "TPC pid all tracks; tpc pid probability; species",100, 0, 1, 5, -0.5, 4.5 );
70da6c5a 241 fOutput->CreateTH2F("hTOF_PID", "TOF pid all tracks; tof pid probability; species",100, 0, 1,5, -0.5, 4.5 );
faee3b18 242 fOutput->CreateTH2F("hTOF_beta_all", "TOF beta for all nice single tracks; p (GeV/c); #beta", 100, 0.1, 10, 480, 0, 1.2, 0);
243 fOutput->CreateTH2F("hTOF_qa_TmT0mT", "TOF (t - t0 - t[pion]) qa verus run", 10000, 114000, 124000, 200, -200, 200);
244 fOutput->CreateTH2F("hTOF_qa_length", "TOF track length verus run", 10000, 114000, 112400, 200, 0, 1000);
245
246 //
247 // debug histograms
248 //
249 fOutput->CreateTH2F("hITS_kFlags", "ITS flags; flag; V0 candidates", 5, 0.5, 5.5, 5, -0.5, 4.5);
250 fOutput->CreateTH2F("hTPC_kFlags", "TPC flags; flag; V0 candidates", 5, 0.5, 5.5, 5, -0.5, 4.5);
251 fOutput->CreateTH2F("hTRD_kFlags", "TRD flags; flag; V0 candidates", 5, 0.5, 5.5, 5, -0.5, 4.5);
252 fOutput->CreateTH2F("hTOF_kFlags", "TOF flags; flag; V0 candidates", 5, 0.5, 5.5, 5, -0.5, 4.5);
70da6c5a 253
254
faee3b18 255 //
256 // event based histograms
257 //
e3fc062d 258 Int_t nT0[2] = {10000, 100};
259 Double_t minT0[2] = {114500, -1000};
260 Double_t maxT0[2] = {124500, 1000};
261 fOutput->CreateTHnSparse("hEvent_T0", "T0 as a function of run number; run number; T0 (ns)", 2, nT0, minT0, maxT0);
faee3b18 262
263 //
264 // test the tender V0 supply
265 //
266 fOutput->CreateTH2F("h_tender_check_01", "tender -vs- HFEpidQA; HFEpidQA V0 candiadates; tender V0 candidates", 4, -1.5, 2.5, 4, -1.5, 2.5);
267 fOutput->CreateTH2Fvector1(3, "h_tender_check_02", "tender -vs- HFEpidQA per type; AliHFEpidQA V0 ; tender V0", 4, -1.5, 2.5, 4, -1.5, 2.5);
70da6c5a 268
70da6c5a 269
70da6c5a 270 //
faee3b18 271 // THnSpasre objects
70da6c5a 272 //
273
faee3b18 274 // Create Illumination Plot
275 // bins: species, pt, eta, phi, TPC status, TRD status
276 {
277 Int_t nbins[6] = {4, 40, 16, 72, 2, 2};
278 Double_t min[6] = { 0, 0.1, -0.8, 0., 0., 0.};
279 Double_t max[6] = { 4., 20., 0.8, 2*TMath::Pi(), 2., 2.};
280 fOutput->CreateTHnSparse("hIllumination", "Illumination", 6, nbins, min, max);
281 fOutput->BinLogAxis("hIllumination", 1);
282 }
283
284 // TPC THnSparse
285 // bins: species, pt, n TPC clusters., TPC electron likelihood, TPC n sigmas, TPC signal
286 {
287 Int_t nbins[6] = { 5, 40, 20, 100, 100, 200};
288 Double_t min[6] = { -0.5, 0.1, 0., 0., -5., 0.};
289 Double_t max[6] = { 4.5, 20., 200, 1., 5., 120.};
290 TString htitle = "TPC info sparse; VO identified species; p (GeV/c); n TPC clusters; TPC N sigma; TPC signal";
291 fOutput->CreateTHnSparse("hTPCclusters",htitle,6, nbins, min, max);
292 fOutput->BinLogAxis("hTPCclusters", 1);
293 }
294 // TRD THnSparse - entries per tracklet
295 // species, p, tracklet position, number of PID tracklets, number of slices (non 0), number of clusters, electron likelihood,
3a72645a 296
faee3b18 297 {
298 Int_t nbins[7] = {5, 20, 6, 7, 9, 45, 100};
299 Double_t min[7] = {-0.5, 0.5, -0.5, -0.5, -0.5, -0.5, 0.};
300 Double_t max[7] = {4.5, 10, 5.5, 6.5, 8.5, 179.5, 1.};
301 TString htitle = "TRD tracklets; VO identified species; p (GeV/c); tracklet position; No. PID tacklets; No. slices; No. clusters; El. likelihood";
302 fOutput->CreateTHnSparse("hTRDtracklets",htitle,7, nbins, min, max);
303 fOutput->BinLogAxis("hTRDtracklets", 1);
70da6c5a 304 }
faee3b18 305
306
307}
308//__________________________________________
309void AliHFEpidQA::Process(){
310 //
311 // Run PID QA
312 //
70da6c5a 313
314 if(!fV0pid){
315 AliError("V0pid not available! Forgotten to initialize?");
316 return;
317 }
faee3b18 318 if(!fESDpid){
319 AliError("fESDpid not initialized, I am leaving this!");
320 return;
321 }
322
323 // to be udpated to AOD save mdoe
324 if(!fEvent){
325 AliError("AliVEvent not available, returning");
326 }
70da6c5a 327
3a72645a 328
70da6c5a 329 if(fMC) fV0pidMC->SetMCEvent(fMC);
3a72645a 330 if(fMC) fV0pid->SetMCEvent(fMC);
70da6c5a 331
faee3b18 332 fV0pid->Process(fEvent);
333 TObjArray *hfeelectrons = fV0pid->GetListOfElectrons();
334 TObjArray *hfepionsK0 = fV0pid->GetListOfPionsK0();
335 TObjArray *hfepionsL = fV0pid->GetListOfPionsL();
336 TObjArray *hfeprotons = fV0pid->GetListOfProtons();
337
338 // Get Track list for normal purpose
339 TObjArray *electrons = MakeTrackList(hfeelectrons);
340 TObjArray *pionsK0 = MakeTrackList(hfepionsK0);
341 TObjArray *pionsL = MakeTrackList(hfepionsL);
342 TObjArray *protons = MakeTrackList(hfeprotons);
343 TObjArray *cleanElectrons = MakeCleanListElectrons(hfeelectrons);
70da6c5a 344
345 if(fMC){
3a72645a 346 fV0pidMC->Process(electrons, AliHFEV0cuts::kRecoElectron);
347 fV0pidMC->Process(pionsK0, AliHFEV0cuts::kRecoPionK0);
348 fV0pidMC->Process(pionsL, AliHFEV0cuts::kRecoPionL);
349 fV0pidMC->Process(protons, AliHFEV0cuts::kRecoProton);
70da6c5a 350 }
351
352 AliDebug(2, Form("Number of Electrons : %d", electrons->GetEntries()));
353 AliDebug(2, Form("Number of K0 Pions : %d", pionsK0->GetEntries()));
354 AliDebug(2, Form("Number of Lambda Pions : %d", pionsL->GetEntries()));
355 AliDebug(2, Form("Number of Protons : %d", protons->GetEntries()));
356 if(fMC){
357 AliDebug(2, "MC Information available. Doing Purity checks...");
358 // Calculate the purity of the clean samples using MC
3a72645a 359 MakePurity(electrons, AliHFEV0cuts::kRecoElectron);
360 MakePurity(pionsK0, AliHFEV0cuts::kRecoPionK0);
361 MakePurity(pionsL, AliHFEV0cuts::kRecoPionL);
362 MakePurity(protons, AliHFEV0cuts::kRecoProton);
70da6c5a 363 }
faee3b18 364
365 // some event wise checks
366 CheckEvent();
367
368 // Make Illumination Plot
3a72645a 369 FillIllumination(electrons, AliHFEV0cuts::kRecoElectron);
370 FillIllumination(pionsK0, AliHFEV0cuts::kRecoPionK0);
371 FillIllumination(pionsL, AliHFEV0cuts::kRecoPionL);
372 FillIllumination(protons, AliHFEV0cuts::kRecoProton);
faee3b18 373
70da6c5a 374 // Now we can do studies on the PID itself
faee3b18 375 // For TRD use the TRD PID QA object
376 fTRDpidQA->ProcessTracks(cleanElectrons, AliPID::kElectron);
377 fTRDpidQA->ProcessTracks(pionsK0, AliPID::kPion);
378 fTRDpidQA->ProcessTracks(pionsL, AliPID::kPion);
379 fTRDpidQA->ProcessTracks(protons, AliPID::kProton);
380
8c1c76e9 381 // Monitor TRD PID Response
382 /*TestTRDResponse(cleanElectrons, AliHFEV0cuts::kRecoElectron);
383 TestTRDResponse(pionsK0, AliHFEV0cuts::kRecoPionK0);
384 TestTRDResponse(pionsL, AliHFEV0cuts::kRecoPionL);
385 TestTRDResponse(protons, AliHFEV0cuts::kRecoProton);*/
386
3a72645a 387 FillElectronLikelihoods(electrons, AliHFEV0cuts::kRecoElectron);
388 FillElectronLikelihoods(pionsK0, AliHFEV0cuts::kRecoPionK0);
389 FillElectronLikelihoods(pionsL, AliHFEV0cuts::kRecoPionL);
390 FillElectronLikelihoods(protons, AliHFEV0cuts::kRecoProton);
70da6c5a 391
3a72645a 392 FillPIDresponse(electrons, AliHFEV0cuts::kRecoElectron);
393 FillPIDresponse(pionsK0, AliHFEV0cuts::kRecoPionK0);
394 FillPIDresponse(pionsL, AliHFEV0cuts::kRecoPionL);
395 FillPIDresponse(protons, AliHFEV0cuts::kRecoProton);
70da6c5a 396
faee3b18 397 // check the tender V0s
3a72645a 398 CheckTenderV0pid(electrons, AliHFEV0cuts::kRecoElectron);
399 CheckTenderV0pid(pionsK0, AliHFEV0cuts::kRecoPionK0);
400 CheckTenderV0pid(pionsL, AliHFEV0cuts::kRecoPionL);
401 CheckTenderV0pid(protons, AliHFEV0cuts::kRecoProton);
faee3b18 402
70da6c5a 403 // Analysis done, flush the containers
404 fV0pid->Flush();
faee3b18 405
406 delete electrons;
407 delete pionsL;
408 delete pionsK0;
409 delete protons;
410 delete cleanElectrons;
411}
412
413//__________________________________________
3a72645a 414void AliHFEpidQA::FillIllumination(const TObjArray * const tracks, Int_t species){
faee3b18 415 //
416 // Fill Illumination Plot
417 //
418 THnSparseF *hIllumination = dynamic_cast<THnSparseF *>(fOutput->Get("hIllumination"));
419 if(!hIllumination) return;
420
421 Double_t quantities[6]; memset(quantities, 0, sizeof(Double_t) *6);
422 TIter trackIter(tracks);
423
424 quantities[0] = species;
425 TObject *o = NULL; AliESDtrack *esdtrack = NULL;
426 while((o = trackIter())){
427 if(!TString(o->IsA()->GetName()).CompareTo("AliESDtrack")){
428 // work on local copy in order to not spoil others
e3ae862b 429 esdtrack = new AliESDtrack(*(static_cast<AliESDtrack *>(o)));
bf892a6a 430 if(!esdtrack) continue;
faee3b18 431 } else if(!TString(o->IsA()->GetName()).CompareTo("AliAODrack")){
432 // Bad hack: Fill ESD track with AOD information
215ffe88 433 esdtrack = new AliESDtrack(static_cast<AliAODTrack *>(o));
bf892a6a 434 if(!esdtrack) continue;
faee3b18 435 } else {
436 // Non usable
437 continue;
438 }
439
440 // Propagate to the entrance of the TRD
441 esdtrack->PropagateTo(300, fEvent->GetMagneticField());
442 quantities[1] = esdtrack->Pt();
443 quantities[2] = esdtrack->Eta();
444 quantities[3] = esdtrack->Phi();
445 quantities[4] = (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) ? 1 : 0;
446 quantities[5] = (esdtrack->GetStatus() & AliESDtrack::kTRDout) ? 1. : 0.;
447 hIllumination->Fill(quantities);
448
449 delete esdtrack;
450 }
451}
452//__________________________________________
453void AliHFEpidQA::FillTPCinfo(AliESDtrack *const esdtrack, Int_t species){
454 //
455 // Fill TPC Cluster Plots
456 //
457 THnSparseF *hTPCclusters = dynamic_cast<THnSparseF *>(fOutput->Get("hTPCclusters"));
458 if(!hTPCclusters) return;
459
460 Double_t quantities[6]; memset(quantities, 0, sizeof(Double_t) *6);
461
462 Double_t pidProbs[5];
463 const Int_t typePID[5] = {0, 2, 2, 3, 4};
464
465
466 quantities[0] = species;
467
468
469 esdtrack->GetTPCpid(pidProbs);
470
471 quantities[1] = esdtrack->P();
472 quantities[2] = esdtrack->GetTPCNcls();
473 quantities[3] = pidProbs[0];
474 quantities[4] = fESDpid->NumberOfSigmasTPC(esdtrack,(AliPID::EParticleType)typePID[species]);
475 quantities[5] = esdtrack->GetTPCsignal();
476 hTPCclusters->Fill(quantities);
477
70da6c5a 478}
479
8c1c76e9 480//__________________________________________
481void AliHFEpidQA::TestTRDResponse(const TObjArray * const tracks, Int_t species){
482 //
483 // Test PID Response function of the TRD
484 //
485 Int_t effInt[6] = {70, 75, 80, 85, 90, 95};
486 const char *sname[5] = {"Electron", "PionK0", "PionL", "Kaon", "Proton"};
487 AliESDtrack *track = NULL;
488 TIter trackIter(tracks);
489 Float_t momenta[6], p;
490 Int_t nmomenta;
491 Double_t probs[5], likeEl;
492 while((track = static_cast<AliESDtrack *>(trackIter()))){
493 if(track->GetTRDntrackletsPID() < 4) continue;
494
495 // calculate momentum at TRD level
496 memset(momenta, 0, sizeof(Float_t) * 6); nmomenta = 0;
497 for(Int_t ily = 0; ily < 6; ily++){
498 if(track->GetTRDmomentum(ily) > 0.01) momenta[nmomenta++] = track->GetTRDmomentum(ily);
499 }
500 p = TMath::Mean(nmomenta, momenta);
501
502 // Get Electron likelihood
503 track->GetTRDpid(probs);
504 likeEl = probs[0]/(probs[0] + probs[2]);
505
506 for(Int_t ieff = 0; ieff < 6; ieff++){
507 if(fESDpid->IdentifiedAsElectronTRD(track, static_cast<Double_t>(effInt[ieff])/100.)) fOutput->Fill(Form("hTRD_likeSel_%s_%dtls_%deff", sname[species], static_cast<Int_t>(track->GetTRDntrackletsPID()), effInt[ieff]), p, likeEl);
508 else fOutput->Fill(Form("hTRD_likeRej_%s_%dtls_%deff", sname[species], static_cast<Int_t>(track->GetTRDntrackletsPID()), effInt[ieff]), p, likeEl);
509 }
510 }
511}
512
70da6c5a 513//__________________________________________
3a72645a 514void AliHFEpidQA::MakePurity(const TObjArray *tracks, Int_t species){
70da6c5a 515 //
516 // Fill the QA histos for a given species
517 //
518 if(!fMC) return;
519 AliDebug(3, Form("Doing Purity checks for species %d", species));
faee3b18 520 Int_t pdg = GetPDG(species);
e3ae862b 521 TString hname;
faee3b18 522
70da6c5a 523 switch(species){
3a72645a 524 case AliHFEV0cuts::kRecoElectron:
e3ae862b 525 hname = "purityElectron";
faee3b18 526 break;
3a72645a 527 case AliHFEV0cuts::kRecoPionK0:
e3ae862b 528 hname = "purityPionK0";
faee3b18 529 break;
3a72645a 530 case AliHFEV0cuts::kRecoPionL:
e3ae862b 531 hname = "purityPionL";
faee3b18 532 break;
3a72645a 533 case AliHFEV0cuts::kRecoProton:
e3ae862b 534 hname = "purityProton";
faee3b18 535 break;
536 default: // non investigated species
537 AliDebug(3, "Species not investigated");
538 return;
539 }
540
70da6c5a 541 AliDebug(3, Form("Number of tracks: %d", tracks->GetEntries()));
faee3b18 542 TIter trackIter(tracks);
70da6c5a 543 AliVParticle *recTrack = NULL, *mcTrack = NULL;
faee3b18 544 while((recTrack = dynamic_cast<AliVParticle *>(trackIter()))){
70da6c5a 545 Int_t label = recTrack->GetLabel();
546 AliDebug(4, Form("MC Label %d", label));
547 mcTrack =fMC->GetTrack(TMath::Abs(label));
548 if(!mcTrack){
549 AliDebug(4, "MC track not available");
550 continue; // we don't know
551 }
552
553 // Get the pdg code
554 Int_t trackPdg = 0;
555 if(!TString(mcTrack->IsA()->GetName()).CompareTo("AliMCParticle")){
556 // case ESD
557 AliMCParticle *mcp = dynamic_cast<AliMCParticle *>(mcTrack);
bf892a6a 558 if(!mcp) continue;
70da6c5a 559 trackPdg = TMath::Abs(mcp->Particle()->GetPdgCode());
560 } else {
561 // case AOD
562 AliAODMCParticle *aodmcp = dynamic_cast<AliAODMCParticle *>(mcTrack);
bf892a6a 563 if(!aodmcp) continue;
70da6c5a 564 trackPdg = TMath::Abs(aodmcp->GetPdgCode());
565 }
566 if(trackPdg == pdg) // Correct identification
faee3b18 567 {
568 fOutput->Fill(hname, 0., recTrack->Pt());
569 }
70da6c5a 570 else // Wrong identification
571 fOutput->Fill(hname, 1., recTrack->Pt());
572 }
70da6c5a 573}
574
575//__________________________________________
3a72645a 576void AliHFEpidQA::FillElectronLikelihoods(const TObjArray * const particles, Int_t species){
70da6c5a 577 //
faee3b18 578 // Fill electron Likelihoods for the ITS, TPC and TOF
70da6c5a 579 // Required for the calculation of the electron efficiency,
580 // pion and proton efficiency and the thresholds
581 //
faee3b18 582 Long_t status = 0;
e3ae862b 583 const TString detname[4] = {"ITS", "TPC", "TRD", "TOF"};
584 TString specname;
faee3b18 585
70da6c5a 586 switch(species){
3a72645a 587 case AliHFEV0cuts::kRecoElectron:
e3ae862b 588 specname = "Electron";
faee3b18 589 break;
3a72645a 590 case AliHFEV0cuts::kRecoPionK0:
e3ae862b 591 specname = "PionK0";
faee3b18 592 break;
3a72645a 593 case AliHFEV0cuts::kRecoPionL:
e3ae862b 594 specname = "PionL";
faee3b18 595 break;
3a72645a 596 case AliHFEV0cuts::kRecoProton:
e3ae862b 597 specname = "Proton";
faee3b18 598 break;
599 default:
600 AliDebug(2, Form("Species %d not investigated", species));
601 return;
70da6c5a 602 };
603 AliVParticle *recTrack = NULL;
faee3b18 604 // mcTrack =fMC->GetTrack(TMath::Abs(label));
605 // if(!mcTrack){
606 // AliDebug(4, "MC track not available");
607 // continue; // we don't know
608 // }
609
610 TIter trackIter(particles);
611
612 Double_t quantities[2];
613 Double_t pidProbs[5];
614
615 while((recTrack = dynamic_cast<AliVParticle *>(trackIter()))){
70da6c5a 616 if(!TString(recTrack->IsA()->GetName()).CompareTo("AliESDtrack")){
617 // case ESD
618 AliESDtrack *esdTrack = dynamic_cast<AliESDtrack *>(recTrack);
bf892a6a 619 if(!esdTrack) continue;
faee3b18 620 status = esdTrack->GetStatus();
621
622 //TPC momentum and likelihoods
623 Double_t pTPC = 0.;
624 pTPC = esdTrack->GetInnerParam() ? esdTrack->GetInnerParam()->P() : esdTrack->P();
625 Bool_t mcFound = kFALSE;
626 if(fMC){
627 Int_t label = esdTrack->GetLabel();
628 AliMCParticle *mcTrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(label));
629 Int_t pdg = GetPDG(species);
630 Int_t trackPdg = 0;
631 if(mcTrack){
632 trackPdg = TMath::Abs(mcTrack->Particle()->GetPdgCode());
633 }
634 if(pdg == trackPdg) mcFound = kTRUE;
70da6c5a 635 }
faee3b18 636 quantities[0] = pTPC;
637 Bool_t detFlagSet = kFALSE;
638 for(Int_t idet = 0; idet < 4; idet++){
e3ae862b 639 TString histname, histnameMC;
e97c2edf 640 histname = "h" + detname[idet] + "_El_like_" + specname;
641 histnameMC = "h" + detname[idet] + "_El_like_MC_" + specname;
e3ae862b 642
faee3b18 643 switch(idet){
644 case kITS: esdTrack->GetITSpid(pidProbs);
645 detFlagSet = status & AliESDtrack::kITSpid;
646 break;
647 case kTPC: esdTrack->GetTPCpid(pidProbs);
648 detFlagSet = status & AliESDtrack::kTPCpid;
649 break;
650 case kTRD: esdTrack->GetTRDpid(pidProbs);
651 detFlagSet = status & AliESDtrack::kTRDpid;
652 break;
653 case kTOF: esdTrack->GetTOFpid(pidProbs);
654 detFlagSet = status & AliESDtrack::kTOFpid;
655 break;
656 };
657 quantities[1] = pidProbs[AliPID::kElectron];
658 // in case of TRD require 6 PID tracklets
659 if(kTRD == idet && esdTrack->GetTRDntrackletsPID() != 6) continue;
660 if(detFlagSet){
661 fOutput->Fill(histname, quantities[0], quantities[1]);
662 if(mcFound)
663 fOutput->Fill(histnameMC, quantities[0], quantities[1]);
664 }
70da6c5a 665 }
faee3b18 666 }//.. ESD
667 else{
668 //AOD
669 }//.. aod
670 }//.. while tracks
70da6c5a 671}
672//__________________________________________
3a72645a 673void AliHFEpidQA::FillPIDresponse(const TObjArray * const particles, Int_t species){
70da6c5a 674 //
675 // Fill the PID response of different detectors to V0 daughter particles
676 //
e3ae862b 677 TString hname, hname2, hname3;
faee3b18 678
e3ae862b 679 const TString typeName[5] = {"Electron", "PionK0", "PionL", "Kaon", "Proton"};
70da6c5a 680 const Int_t typePID[5] = {0, 2, 2, 3, 4};
faee3b18 681
682 // PID THnSparse
683 // axes:
684 // 0) species, 1) momentum, 2) DCA xy, 3) DCA z
685 // 4) ITS signal
686 // 5) TPC Ncls 6) TPC signal 7) TPC nSigma,
687 // 8) TRD Ntrk, 9) TRD Ncls, 10) TRD dEdx,
688
8c1c76e9 689 //Double_t data[12];
bf892a6a 690
faee3b18 691
692 Int_t run = fEvent->GetRunNumber();
693
70da6c5a 694 AliVParticle *recTrack = NULL;
faee3b18 695 TIter trackIter(particles);
696 while((recTrack = dynamic_cast<AliVParticle *>(trackIter()))){
8c1c76e9 697 //for(Int_t i=0; i<12; ++i) data[i] = -99.;
70da6c5a 698 // ESD
699 if(!TString(recTrack->IsA()->GetName()).CompareTo("AliESDtrack")){
700 // case ESD
701 AliESDtrack *esdTrack = dynamic_cast<AliESDtrack *>(recTrack);
faee3b18 702 if(!esdTrack) continue;
70da6c5a 703
faee3b18 704 // for the PID THnSparse
8c1c76e9 705 //data[0] = species;
706 //data[1] = esdTrack->P();
faee3b18 707 Float_t impactR = -1.;
708 Float_t impactZ = -1.;
709 esdTrack->GetImpactParameters(impactR, impactZ);
8c1c76e9 710 //data[2] = impactR;
711 //data[3] = impactZ;
712 //data[11] = 0; // initialize the TOF pid cut on elecgrons to false
faee3b18 713 // use ONLY tracks with PID flag TRUE
714 ULong_t status = 0;
715 status = esdTrack->GetStatus();
716
717 //
718 // DEBUG
719 //
720
721 fOutput->Fill("hITS_kFlags", 5., species);
722 if(status & AliESDtrack::kITSin) fOutput->Fill("hITS_kFlags", 1., species);
723 if(status & AliESDtrack::kITSout) fOutput->Fill("hITS_kFlags", 2., species);
724 if(status & AliESDtrack::kITSrefit) fOutput->Fill("hITS_kFlags", 3., species);
725 if(status & AliESDtrack::kITSpid) fOutput->Fill("hITS_kFlags", 4., species);
726
727 fOutput->Fill("hTPC_kFlags", 5., species);
728 if(status & AliESDtrack::kTPCin) fOutput->Fill("hTPC_kFlags", 1., species);
729 if(status & AliESDtrack::kTPCout) fOutput->Fill("hTPC_kFlags", 2., species);
730 if(status & AliESDtrack::kTPCrefit) fOutput->Fill("hTPC_kFlags", 3., species);
731 if(status & AliESDtrack::kTPCpid) fOutput->Fill("hTPC_kFlags", 4., species);
732
733 fOutput->Fill("hTRD_kFlags", 5., species);
734 if(status & AliESDtrack::kTRDin) fOutput->Fill("hTRD_kFlags", 1., species);
735 if(status & AliESDtrack::kTRDout) fOutput->Fill("hTRD_kFlags", 2., species);
736 if(status & AliESDtrack::kTRDrefit) fOutput->Fill("hTRD_kFlags", 3., species);
737 if(status & AliESDtrack::kTRDpid) fOutput->Fill("hTRD_kFlags", 4., species);
738
739 fOutput->Fill("hTOF_kFlags", 5., species);
740 if(status & AliESDtrack::kTOFin) fOutput->Fill("hTOF_kFlags", 1., species);
741 if(status & AliESDtrack::kTOFout) fOutput->Fill("hTOF_kFlags", 2., species);
742 if(status & AliESDtrack::kTOFrefit) fOutput->Fill("hTOF_kFlags", 3., species);
743 if(status & AliESDtrack::kTOFpid) fOutput->Fill("hTOF_kFlags", 4., species);
744
745
746 //
747 // ITS -
748 //
749 if(status & AliESDtrack::kITSpid){
750 Double_t p = esdTrack->P();
751
752 // ITS signal
753 //Double_t itsSignal = esdTrack->GetITSsignal();
754
755 // ITS dEdx
756 Double_t dEdxSamples[4];
757 esdTrack->GetITSdEdxSamples(dEdxSamples);
758 Int_t nSamples = 0;
759 Double_t dEdxSum = 0.;
e3ae862b 760 hname = "hITS_dEdx_" + typeName[species];
faee3b18 761 for(Int_t i=0; i<4; ++i){
762 if(dEdxSamples[i] > 0){
763 nSamples++;
764 fOutput->Fill(hname, i+1, p, dEdxSamples[i]);
765 dEdxSum += dEdxSamples[i];
766 }
70da6c5a 767 }
faee3b18 768 if(4 == nSamples)fOutput->Fill(hname, 0, p, dEdxSum);
769 fOutput->Fill("hITS_dEdx_nSamples", nSamples);
770
771 Double_t signal = esdTrack->GetITSsignal();
e3ae862b 772 hname = "hITS_Signal_" + typeName[species];
faee3b18 773 fOutput->Fill(hname, p, signal);
8c1c76e9 774 //data[4] = signal;
faee3b18 775
776 // ITS number of signas
777 Double_t nsigma = fESDpid->NumberOfSigmasITS(esdTrack,(AliPID::EParticleType)typePID[species]);
e3ae862b 778 hname = "hITS_nSigma_" + typeName[species];
faee3b18 779 fOutput->Fill(hname, p, nsigma);
faee3b18 780 // ITS PID response
781 Double_t itsPID[5] = {-1, -1, -1, -1, -1};
782 esdTrack->GetITSpid(itsPID);
783 Int_t ix = GetMaxPID(itsPID);
e97c2edf 784 hname = "hITS_PID_p_" + typeName[species];
91c7e1ec 785 fOutput->Fill(hname, p, ix);
faee3b18 786 }//.. kITSpid
787
788 //
789 // TPC
790 //
791 if(status & AliESDtrack::kTPCpid){
792 // Make TPC clusters Plot
8c1c76e9 793 //data[5] = esdTrack->GetTPCNcls();
faee3b18 794 FillTPCinfo(esdTrack, species);
70da6c5a 795
faee3b18 796 Double_t p = esdTrack->GetInnerParam() ? esdTrack->GetInnerParam()->P() : esdTrack->P();
797 // TPC dEdx
798 Double_t dEdx = esdTrack->GetTPCsignal();
e3ae862b 799 hname = "hTPC_dEdx_" + typeName[species];
faee3b18 800 fOutput->Fill(hname, p, dEdx);
8c1c76e9 801 //data[6] = dEdx;
faee3b18 802
803 //TPC number of sigmas
804 Double_t nsigma = fESDpid->NumberOfSigmasTPC(esdTrack,(AliPID::EParticleType)typePID[species]);
e3ae862b 805 hname = "hTPC_nSigma_" + typeName[species];
faee3b18 806 fOutput->Fill(hname, p, nsigma);
8c1c76e9 807 //data[7] = nsigma;
faee3b18 808
809 // TPC PID response
e3ae862b 810 hname = "hTPC_PID_p_" + typeName[species];
faee3b18 811 Double_t tpcPID[5] = {-1, -1, -1, -1, -1};
812 esdTrack->GetTPCpid(tpcPID);
813 Int_t ix = GetMaxPID(tpcPID);
814 fOutput->Fill(hname, p, ix);
815 }//.. kTPCpid
816
817 //
818 // TRD
819 //
820
821 if(status & AliESDtrack::kTRDpid){
822 Double_t p = esdTrack->GetOuterParam() ? esdTrack->GetOuterParam()->P() : esdTrack->P();
823
824 // TRD number of tracklets
825 Int_t ntrk = esdTrack->GetTRDntrackletsPID();
e3ae862b 826 hname = "hTRD_trk_" + typeName[species];
faee3b18 827 fOutput->Fill(hname, p, ntrk);
8c1c76e9 828 //data[8] = ntrk;
faee3b18 829
830 // TRD PID response
e3ae862b 831 hname = "hTRD_PID_p_" + typeName[species];
faee3b18 832 Double_t trdPID[5] = {-1., -1., -1., -1., -1};
833 esdTrack->GetTRDpid(trdPID);
834 Int_t ix = GetMaxPID(trdPID);
91c7e1ec 835 fOutput->Fill(hname, p, ix);
faee3b18 836 // TRD n clusters
837 Int_t ncls = esdTrack->GetTRDncls();
e3ae862b 838 hname = "hTRD_cls_" + typeName[species];
faee3b18 839 fOutput->Fill(hname, p, ncls);
8c1c76e9 840 //data[9] = ncls;
faee3b18 841
842 // TRD - per tracklet - dEdx per, likelihood
e3ae862b 843 hname = "hTRD_Nslc_" + typeName[species];
844 hname2 = "hTRD_slc_" + typeName[species];
845 hname3 = "hTRD_dEdx_" + typeName[species];
faee3b18 846 Int_t nSlices = esdTrack->GetNumberOfTRDslices();
847 Double_t sumTot = 0.;
848 Int_t not0Tot = 0;
849 for(Int_t l=0; l< 6; ++l){
850 Double_t trkData[7] = {-1.,-1, -1, -1, -1, -1, -1};
851 trkData[0] = species;
852 trkData[1] = p;
853 trkData[2] = l;
854 trkData[3] = ntrk;
855 trkData[5] = ncls;
856 Double_t sum = 0.;
857 Int_t not0 = 0;
858 for(Int_t s=0; s<nSlices; ++s){
859 Double_t slice = esdTrack->GetTRDslice(l, s);
860 sum += slice;
861 if(slice > 0){
862 not0 += 1;
863 fOutput->Fill(hname2, p, s);
864 }
865 }//..slices
866
867 trkData[4] = not0;
868 fOutput->Fill(hname, p, not0);
869 fOutput->Fill(hname3, p, sum);
870 if(sum > 0){
871 sumTot += sum;
872 not0Tot += 1;
873 }
874 // lkelihoods per layer
875 if(not0Tot > 0 && fNNref){
876 Double_t likelihoods[5] = {-1., -1., -1., -1., -1};
877 TRDlikeTracklet(l, esdTrack, likelihoods);
878 trkData[6] = likelihoods[0];
879 //printf(" -D: species: %i, P; %f : %f, s: %f\n", species, p, likelihoods[0], s);
880 }
881 if(not0Tot) fOutput->Fill("hTRDtracklets", trkData);
882 }//..layers
883 // average dEx per number of tracklets
8c1c76e9 884 //if(0 < not0Tot)
885 // data[10] = sumTot / not0Tot;
faee3b18 886 }//.. kTRDpid
70da6c5a 887
faee3b18 888
889 //
890 // TOF
891 //
892 if(status & AliESDtrack::kTOFpid){
893 Double_t p = esdTrack->GetOuterParam() ? esdTrack->GetOuterParam()->P() : esdTrack->P();
bf892a6a 894 Double_t t0 = fESDpid->GetTOFResponse().GetStartTime(esdTrack->P());
faee3b18 895
896 //TOF beta
e3ae862b 897 hname = "hTOF_beta_" + typeName[species];
faee3b18 898 Float_t beta = TOFbeta(esdTrack);
899 fOutput->Fill(hname, p, beta);
900 fOutput->Fill("hTOF_beta_all", p, beta);
901 // TOF nSigma
8c1c76e9 902 Double_t nsigma = fESDpid->NumberOfSigmasTOF(esdTrack,(AliPID::EParticleType)typePID[species]);
e3ae862b 903 hname = "hTOF_nSigma_" + typeName[species];
faee3b18 904 fOutput->Fill(hname, p, nsigma);
8c1c76e9 905 //if(beta > 0.97 && beta < 1.03){
906 // data[11] = 1;
907 //}
faee3b18 908
909 // TOF PID response
e3ae862b 910 hname = "hTOF_PID_p_" + typeName[species];
faee3b18 911 Double_t tofPID[5] = {-1., -1., -1., -1., -1};
912 esdTrack->GetTOFpid(tofPID);
913 Int_t ix = GetMaxPID(tofPID);
914 fOutput->Fill(hname, p, ix);
915
916 // time of flight QA
917 // - distribution of (time - t0 - pion_time)
918 Double_t times[5];
919 esdTrack->GetIntegratedTimes(times);
920 Double_t tItrackL = esdTrack->GetIntegratedLength();
921 Double_t tTOFsignal = esdTrack->GetTOFsignal();
922 Double_t dT = tTOFsignal - t0 - times[2];
923 fOutput->Fill("hTOF_qa_TmT0mT", run*1.0, dT);
924 fOutput->Fill("hTOF_qa_length", run*1.0, tItrackL);
70da6c5a 925
faee3b18 926
927 }//.. kTOFpid
928 // temporary - the PIDsparse needs rebuilding
929 //fOutput->Fill("PIDsparse", data);
70da6c5a 930 }
931 // AOD - comming soon
932 else{
933 continue;
934 }
935 }// .. tracks in TObjArray
936
70da6c5a 937
938}
faee3b18 939//__________________________________________
940void AliHFEpidQA:: CheckEvent(){
941 //
942 // check some event variables
943 //
70da6c5a 944
faee3b18 945 // check the T0 as a function of run number (less than one bin per number
946 Double_t t0 = fESDpid->GetTOFResponse().GetTimeZero();
947 Int_t run = fEvent->GetRunNumber();
948 Double_t data[2] = {run*1.0, t0*1000.};
949 fOutput->Fill("hEvent_T0", data);
950
951
952}
70da6c5a 953//__________________________________________
954TList *AliHFEpidQA::GetOutput(){
955 //
956 // Getter for Output histograms
957 //
958 return fOutput->GetList();
959}
960
961//__________________________________________
962TList *AliHFEpidQA::GetV0pidQA(){
963 //
964 // Getter for V0 PID QA histograms
965 //
966 return fV0pid->GetListOfQAhistograms();
967}
968
969//__________________________________________
970TList *AliHFEpidQA::GetV0pidMC(){
971 //
972 // Getter for V0 PID QA histograms
973 //
974 if(fV0pidMC)
975 return fV0pidMC->GetListOfQAhistograms();
976 return NULL;
977}
978
979//__________________________________________
980void AliHFEpidQA::RecalculateTRDpid(AliESDtrack * /*track*/, Double_t * /*pidProbs*/) const{
faee3b18 981 // fTRDpidResponse->MakePID(track);
982 // track->GetTRDpid(pidProbs);
70da6c5a 983}
984
985//__________________________________________
986void AliHFEpidQA::RecalculateTRDpid(AliAODTrack * /*track*/, Double_t * /*pidProbs*/) const{
faee3b18 987 // fTRDpidResponse->MakePID(track, pidProbs);
70da6c5a 988}
989//___________________________________________________________________
e156c3bb 990Float_t AliHFEpidQA::TOFbeta(const AliESDtrack * const track) const {
70da6c5a 991 // computes the TOF beta
992 Double_t l = track->GetIntegratedLength(); // cm
993 Double_t t = track->GetTOFsignal();
faee3b18 994 Double_t t0 = fESDpid->GetTOFResponse().GetTimeZero(); // ps
995
996 //printf("-D: l: %f, t: %f, t0: %f\n", l, t, t0);
997
70da6c5a 998 if(l < 360. || l > 800.) return 0.;
999 if(t <= 0.) return 0.;
faee3b18 1000 if(t0 >999990.0) return 0.;
1001
70da6c5a 1002
1003 t -= t0; // subtract the T0
1004
1005 l *= 0.01; // cm ->m
1006 t *= 1e-12; //ps -> s
1007
1008
1009 Double_t v = l / t;
1010 Float_t beta = v / TMath::C();
1011
1012 return beta;
1013}
faee3b18 1014//____________________________________________
3a72645a 1015Int_t AliHFEpidQA::GetMaxPID(const Double_t *pidProbs) const {
faee3b18 1016 //
1017 // return the index of maximal PID probability
1018 //
1019 Int_t ix = -1;
1020 Double_t tmp = 0.2;
1021 for(Int_t i=0; i<5; ++i){
1022 if(pidProbs[i] > tmp){
1023 ix = i;
1024 tmp = pidProbs[i];
1025 }
1026 }
1027 return ix;
1028}
1029//_____________________________________________
1030Int_t AliHFEpidQA::GetPDG(Int_t species){
1031 //
1032 // return the PDG particle code
1033 //
1034
1035 Int_t pdg = 0;
1036
1037 switch(species){
3a72645a 1038 case AliHFEV0cuts::kRecoElectron:
faee3b18 1039 pdg = TMath::Abs(kElectron);
1040 break;
3a72645a 1041 case AliHFEV0cuts::kRecoPionK0:
faee3b18 1042 pdg = TMath::Abs(kPiPlus);
1043 break;
3a72645a 1044 case AliHFEV0cuts::kRecoPionL:
faee3b18 1045 pdg = TMath::Abs(kPiPlus);
1046 break;
3a72645a 1047 case AliHFEV0cuts::kRecoProton:
faee3b18 1048 pdg = TMath::Abs(kProton);
1049 break;
1050 default: // non investigated species
1051 AliDebug(3, "Species not recognised");
1052 return 0;
1053 }
1054
1055 return pdg;
1056
1057}
1058
1059//_____________________________________________
3a72645a 1060TObjArray * AliHFEpidQA::MakeTrackList(const TObjArray *tracks) const {
faee3b18 1061 //
1062 // convert list of AliHFEV0Info into a list of AliVParticle
1063 //
1064 TObjArray *output = new TObjArray;
1065 TIter trackInfos(tracks);
1066 AliHFEV0info *trackInfo = NULL;
1067 while((trackInfo = dynamic_cast<AliHFEV0info *>(trackInfos())))
1068 output->Add(trackInfo->GetTrack());
1069
1070 return output;
1071}
1072
1073//_____________________________________________
3a72645a 1074TObjArray * AliHFEpidQA::MakeCleanListElectrons(const TObjArray *electrons) const {
faee3b18 1075 //
1076 // Cleanup electron sample using TPC PID
1077 // PID requirement will allways be implemented to the pair
1078 // Strategy
1079 //
1080 TObjArray *tracks = new TObjArray;
1081 TIter candidates(electrons);
fd282aa0 1082 AliESDEvent *esd;
1083 //AliAODEvent *aod;
faee3b18 1084 AliHFEV0info *hfetrack;
e3fc062d 1085 // const Double_t kSigmaTight = 1.;
1086 // const Double_t kSigmaLoose = 4.;
3a72645a 1087 const Double_t kSigmaTight = 2.;
1088 const Double_t kSigmaLoose = 2.;
1089 const Double_t shift = 0.57;
faee3b18 1090 if((esd = dynamic_cast<AliESDEvent *>(fEvent))){
1091 AliESDtrack *track = NULL, *partnerTrack = NULL;
1092 while((hfetrack = dynamic_cast<AliHFEV0info *>(candidates()))){
1093 track = dynamic_cast<AliESDtrack *>(hfetrack->GetTrack());
bf892a6a 1094 if(!track) continue;
faee3b18 1095 partnerTrack = esd->GetTrack(hfetrack->GetPartnerID());
e3fc062d 1096 Double_t nSigmaTrack = TMath::Abs(fESDpid->NumberOfSigmasTPC(track, AliPID::kElectron) - shift);
1097 Double_t nSigmaPartner = TMath::Abs(fESDpid->NumberOfSigmasTPC(partnerTrack, AliPID::kElectron) - shift);
1098 if((nSigmaTrack < kSigmaTight && nSigmaPartner < kSigmaLoose) || (nSigmaTrack < kSigmaLoose && nSigmaPartner < kSigmaTight))
faee3b18 1099 tracks->Add(track);
1100 }
3e03bd67 1101 }
1102 /*else {
faee3b18 1103 aod = dynamic_cast<AliAODEvent *>(fEvent);
bf892a6a 1104 if(!aod) return NULL;
1105 //AliAODTrack *track = NULL, *partnerTrack = NULL;
faee3b18 1106 while((hfetrack = dynamic_cast<AliHFEV0info *>(candidates()))){
bf892a6a 1107 if(!hfetrack) continue;
1108 //track = dynamic_cast<AliAODTrack *>(hfetrack->GetTrack());
1109 //partnerTrack = aod->GetTrack(hfetrack->GetPartnerID());
faee3b18 1110 // will be coming soon
1111 }
3e03bd67 1112 }*/
faee3b18 1113 return tracks;
1114}
1115//___________________________________________________________
3a72645a 1116void AliHFEpidQA::CheckTenderV0pid(const TObjArray * const particles, Int_t species){
faee3b18 1117
1118 //
1119 // retrieve the status bits from the TObject used to flag
1120 // the V0 daughter tracks and return the found PID
1121 // 0 - electron, 1 - pion, 2 - proton
1122 //
1123
1124 const Int_t id[5] = {0, 1, 1, -1, 2}; // convert AliHFEpid to simple 0, 1, 2
1125
1126 AliVParticle *recTrack = NULL;
1127 TIter trackIter(particles);
1128 while((recTrack = dynamic_cast<AliVParticle *>(trackIter()))){
1129 if(!TString(recTrack->IsA()->GetName()).CompareTo("AliESDtrack")){
1130 // case ESD
1131 AliESDtrack *track = dynamic_cast<AliESDtrack *>(recTrack);
1132 if(!track) continue;
1133 Int_t tPID = GetTenderV0pid(track);
1134 fOutput->Fill("h_tender_check_01", id[species]*1.0, tPID);
1135 fOutput->Fill("h_tender_check_02", id[species], id[species]*1.0, tPID);
1136
1137 } //.. case ESD
1138
1139 }//.. iterate of tracks
1140}
1141//___________________________________________________________
1142Int_t AliHFEpidQA::GetTenderV0pid(AliESDtrack * const track){
1143 //
1144 // retrieve the PID nformation stored in the status flags by the train tender
1145 //
1146
1147
1148 Int_t pid = -1;
1149 if(!track){
1150 return pid;
1151 }
1152
1153 Int_t nTimes = 0;
1154
1155 if(track->TestBit(2<<14)){
1156 pid = 0;
1157 nTimes++;
1158 }
1159 if(track->TestBit(2<<15)){
1160 pid = 1;
1161 nTimes++;
1162 }
1163 if(track->TestBit(2<<16)){
1164 pid = 2;
1165 nTimes++;
1166 }
1167
1168 if(nTimes > 1){
1169 AliWarning("V0 track labeled multiple times by the V0 tender");
1170 pid = -1;
1171 }
1172
1173 return pid;
1174
1175}
1176//___________________________________________________________
1177Double_t AliHFEpidQA::TRDlikeTracklet(Int_t layer, AliESDtrack * const track, Double_t *likelihood){
1178 //
1179 // compute the TRD electron likelihoods for 1 tracklet
1180 // based on teh AliTRDpidRecalculator in train/until/tender
1181 // returns sum of the likelihoods (which should be 1)
1182 //
1183
1184 const Double_t cScaleGain = 1./ 16000.;
1185 const Float_t pBins[11] ={0.6, 0.8, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0}; // momentum bins
1186
1187 if(!track) return kFALSE;
1188 Float_t p = track->GetTRDmomentum(layer); // momentum for a tracklet in the ESDtrack
1189 if(p < 0) return kFALSE;
1190
1191 Int_t mombin = TRDmomBin(p); // momentum bin
bf892a6a 1192 if(mombin < 0) return -1.;
faee3b18 1193 Float_t dEdxTRDsum = 0; // dEdxTRDsum for checking if tracklet is available
1194 Float_t dEdxTRD[8]; // dEdx for a tracklet in the ESD slices
1195 Double_t ddEdxTRD[8]; // dEdx as Double_t for TMultiLayerPerceptron::Evaluate()
1196
1197 Double_t prob[AliPID::kSPECIES]; // probabilities for all species in all layers
1198
1199 for(Int_t is = 0; is < AliPID::kSPECIES; is++){
1200 likelihood[is] = 0.2; // init probabilities
1201 prob[is] = 0.2;
1202 }
1203
1204 Double_t sum = 0.;
1205
1206 for(Int_t islice = 0; islice<8; islice++){
1207 dEdxTRD[islice]=0.; // init dE/dx
1208 ddEdxTRD[islice]=0.; // init dE/dx
1209 dEdxTRD[islice]=track->GetTRDslice(layer,islice); // get the dE/dx values
1210 dEdxTRDsum += dEdxTRD[islice];
1211 ddEdxTRD[islice]=(Double_t)dEdxTRD[islice]*cScaleGain; // rescale dE/dx
1212
1213 }
1214 for(Int_t is = 0; is < AliPID::kSPECIES; is++){
1215 Double_t probloc1, probloc2;
1216 if(mombin == 0 && mombin < pBins[0]){ // calculate PID for p > 0.6 GeV/c
1217 prob[is] = fNet[mombin]->Evaluate(is, ddEdxTRD);
1218 }
1219 else if(mombin == 10 && mombin >= pBins[10]){ // calculate PID for p >= 10.0 GeV/c
1220 prob[is] = fNet[mombin]->Evaluate(is, ddEdxTRD);
1221 }
1222 else{ // standard calculation
1223 Int_t mombin1 = 0, mombin2 = 0; // lower and upper momentum bin
1224 if(p < pBins[mombin]) {mombin1 = mombin -1; mombin2 = mombin;}
1225 if(p >= pBins[mombin]) {mombin1 = mombin; mombin2 = mombin+1;}
1226 probloc1 = fNet[mombin1]->Evaluate(is, ddEdxTRD);
1227 probloc2 = fNet[mombin2]->Evaluate(is, ddEdxTRD);
1228 // weighting of the probability with regard to the track momentum
1229 prob[is] = probloc1 + (probloc2-probloc1)*(p-pBins[mombin1])/(pBins[mombin2]-pBins[mombin1]);
1230 }
1231 likelihood[is] = prob[is];
1232 sum += likelihood[is];
1233 }
1234
1235 return sum;
1236}
1237//__________________________________________________________________________
6555e2ad 1238Int_t AliHFEpidQA::TRDmomBin(Double_t p) const {
faee3b18 1239 //
1240 // compute the momentum bin position
1241 //
1242
1243 const Float_t pBins[11] ={0.6, 0.8, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0}; // momentum bins
1244
1245 Int_t pBin1 = -1; // check bin1
1246 Int_t pBin2 = -1; // check bin2
1247
1248 if(p < 0) return -1; // return -1 if momentum < 0
1249 if(p < pBins[0]) return 0; // smallest momentum bin
1250 if(p >= pBins[10]) return 10; // largest momentum bin
1251
1252
1253 // calculate momentum bin for non extremal momenta
1254 for(Int_t iMomBin = 1; iMomBin < 11; iMomBin++){
1255 if(p < pBins[iMomBin]){
1256 pBin1 = iMomBin - 1;
1257 pBin2 = iMomBin;
1258 }
1259 else
1260 continue;
1261
1262 if(p - pBins[pBin1] >= pBins[pBin2] - p){
1263 return pBin2;
1264 }
1265 else{
1266 return pBin1;
1267 }
1268 }
1269
1270 return -1;
1271
1272
1273}
1274//__________________________________________________________________________
1275
1276