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