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