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