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