1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 // Several studies done on clean samples of electrons, pions and kaons
19 // Compatible with both ESDs and AODs
22 // Matus Kalisky <matus.kalisky@cern.ch>
23 // Markus Heide <mheide@uni-muenster.de>
24 // Markus Fasel <M.Fasel@gsi.de>
29 #include <TObjArray.h>
32 #include <TMultiLayerPerceptron.h>
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 "AliMultiplicity.h"
44 #include "AliPIDResponse.h"
45 #include "AliVParticle.h"
48 #include "AliHFEcollection.h"
49 #include "AliHFEpidQA.h"
50 #include "AliHFEV0info.h"
51 #include "AliHFEV0pid.h"
52 #include "AliHFEV0pidMC.h"
53 #include "AliHFEV0cuts.h"
54 #include "AliHFEtrdPIDqa.h"
58 //__________________________________________
59 AliHFEpidQA::AliHFEpidQA():
68 fTRDTotalChargeInSlice0(kFALSE),
73 // Default constructor
75 for(Int_t mom = 0; mom < 11; mom++){
80 //__________________________________________
81 AliHFEpidQA::AliHFEpidQA(const AliHFEpidQA &ref):
91 fTRDTotalChargeInSlice0(ref.fTRDTotalChargeInSlice0),
98 for(Int_t mom = 0; mom < 11; mom++){
104 //__________________________________________
105 AliHFEpidQA &AliHFEpidQA::operator=(const AliHFEpidQA &ref){
107 // Assignment operator
114 //__________________________________________
115 AliHFEpidQA::~AliHFEpidQA(){
120 // the pointers bellow are not dfeleted to prevend double-deleting of some of the content
121 // these pointers are defined only once during the program call and should not cause a problem,
122 // but cleaner solution is necessary.
123 //if(fV0pid) delete fV0pid;
124 //if(fV0pidMC) delete fV0pidMC;
125 //if(fOutput) delete fOutput;
127 // if(fTRDpidResponse) delete fTRDpidResponse;
130 //__________________________________________
131 void AliHFEpidQA::Copy(TObject &o) const {
138 AliHFEpidQA &target = dynamic_cast<AliHFEpidQA &>(o);
140 target.fTRDTotalChargeInSlice0 = fTRDTotalChargeInSlice0;
142 if(target.fESDpid) delete target.fESDpid;
143 target.fESDpid = fESDpid;
144 if(target.fV0pid) delete target.fV0pid;
146 target.fV0pid = dynamic_cast<AliHFEV0pid *>(fV0pid->Clone());
148 target.fV0pid = NULL;
149 if(target.fV0pidMC) delete target.fV0pidMC;
151 target.fV0pidMC = dynamic_cast<AliHFEV0pidMC *>(fV0pidMC->Clone());
153 target.fV0pidMC = NULL;
154 if(target.fTRDpidQA) delete target.fTRDpidQA;
156 target.fTRDpidQA = dynamic_cast<AliHFEtrdPIDqa *>(fTRDpidQA->Clone());
158 target.fTRDpidQA = NULL;
159 if(target.fOutput) delete target.fOutput;
161 target.fOutput = dynamic_cast<AliHFEcollection *>(fOutput->Clone());
164 //__________________________________________
165 void AliHFEpidQA::Init(){
167 // Prepare task output
172 for(Int_t mom = 0; mom < 11; mom++){
173 fNet[mom] = (TMultiLayerPerceptron*) fNNref->Get(Form("NN_Mom%d", mom));
175 AliError(Form("No reference network for momentum bin %d!", mom));
180 fV0pid = new AliHFEV0pid("fV0pid");
181 if(HasV0pidQA()) fV0pid->InitQA();
182 fV0pidMC = new AliHFEV0pidMC();
185 fTRDpidQA = new AliHFEtrdPIDqa;
186 if(fTRDTotalChargeInSlice0) fTRDpidQA->SetTotalChargeInSlice0();
189 fOutput = new AliHFEcollection("pidQA", "PID QA output");
191 const char *name[4] = {"Electron", "PionK0", "PionL", "Proton"};
192 const char *title[4] = {"Electron", "K0 Pion", "Lambda Pion", "Proton"};
193 const char *det[4] = {"ITS", "TPC", "TRD", "TOF"};
194 const int effs[6] = {70, 75, 80, 85, 90, 95};
195 for(Int_t i = 0; i < 4; i++){
196 fOutput->CreateTH2F(Form("purity%s", name[i]), Form("%s Purity", title[i]), 2, -0.5, 1.5, 20, 0.1, 10, 1);
198 for(Int_t idet = 0; idet < 4; idet++){
199 // create all the histograms which all the detectors have in common
200 if(idet != 2){ // no nSigma histogram for TRD
201 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);
203 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);
204 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);
205 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);
211 fOutput->CreateTH2F(Form("hITS_Signal_%s", name[i]), Form("ITS Signal org. for %ss", title[i]), 40, 0.1, 10, 400, 0, 1000, 0);
212 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);
217 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);
222 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);
223 // number of the non 0 slices per tracklet
224 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);
225 // location of the slices > 0 - where are the emtpy slices located ?
226 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);
227 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);
228 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);
230 for(Int_t itl = 4; itl < 6; itl++){
231 for(Int_t ieff = 0; ieff < 6; ieff++){
232 fOutput->CreateTH2F(Form("hTRD_likeSel_%s_%dtls_%deff", name[i], itl, effs[ieff]), Form(" %s selected as electrons with %d tracklets and %f electron efficiency", title[i], itl, static_cast<Double_t>(effs[ieff])/100.), 44, 0.1, 20., 100, 0., 1.);
233 fOutput->CreateTH2F(Form("hTRD_likeRej_%s_%dtls_%deff", name[i], itl, effs[ieff]), Form(" %s rejected as electrons with %d tracklets and %f electron efficiency", title[i], itl, static_cast<Double_t>(effs[ieff])/100.), 44, 0.1, 20., 100, 0., 1.);
239 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);
240 }//.. loop over identified particle species
243 fOutput->CreateTH1F("hITS_dEdx_nSamples", "ITS - number of non 0 dEdx samples", 4, 0.5, 4.5);
244 fOutput->CreateTH2F("hTPC_PID", "TPC pid all tracks; tpc pid probability; species",100, 0, 1, 5, -0.5, 4.5 );
245 fOutput->CreateTH2F("hTOF_PID", "TOF pid all tracks; tof pid probability; species",100, 0, 1,5, -0.5, 4.5 );
246 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);
247 fOutput->CreateTH2F("hTOF_qa_TmT0mT", "TOF (t - t0 - t[pion]) qa verus run", 10000, 114000, 124000, 200, -200, 200);
248 fOutput->CreateTH2F("hTOF_qa_length", "TOF track length verus run", 10000, 114000, 112400, 200, 0, 1000);
253 fOutput->CreateTH2F("hITS_kFlags", "ITS flags; flag; V0 candidates", 5, 0.5, 5.5, 5, -0.5, 4.5);
254 fOutput->CreateTH2F("hTPC_kFlags", "TPC flags; flag; V0 candidates", 5, 0.5, 5.5, 5, -0.5, 4.5);
255 fOutput->CreateTH2F("hTRD_kFlags", "TRD flags; flag; V0 candidates", 5, 0.5, 5.5, 5, -0.5, 4.5);
256 fOutput->CreateTH2F("hTOF_kFlags", "TOF flags; flag; V0 candidates", 5, 0.5, 5.5, 5, -0.5, 4.5);
260 // event based histograms
262 Int_t nT0[2] = {10000, 100};
263 Double_t minT0[2] = {114500, -1000};
264 Double_t maxT0[2] = {124500, 1000};
265 fOutput->CreateTHnSparse("hEvent_T0", "T0 as a function of run number; run number; T0 (ns)", 2, nT0, minT0, maxT0);
268 // test the tender V0 supply
270 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);
271 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);
278 // Create Illumination Plot
279 // bins: species, pt, eta, phi, TPC status, TRD status
281 Int_t nbins[6] = {4, 40, 16, 72, 2, 2};
282 Double_t min[6] = { 0, 0.1, -0.8, 0., 0., 0.};
283 Double_t max[6] = { 4., 20., 0.8, 2*TMath::Pi(), 2., 2.};
284 fOutput->CreateTHnSparse("hIllumination", "Illumination", 6, nbins, min, max);
285 fOutput->BinLogAxis("hIllumination", 1);
289 // bins: species, pt, n TPC clusters., TPC electron likelihood, TPC n sigmas, TPC signal
291 Int_t nbins[6] = { 5, 40, 20, 100, 100, 200};
292 Double_t min[6] = { -0.5, 0.1, 0., 0., -5., 0.};
293 Double_t max[6] = { 4.5, 20., 200, 1., 5., 120.};
294 TString htitle = "TPC info sparse; VO identified species; p (GeV/c); n TPC clusters; TPC N sigma; TPC signal";
295 fOutput->CreateTHnSparse("hTPCclusters",htitle,6, nbins, min, max);
296 fOutput->BinLogAxis("hTPCclusters", 1);
298 // TRD THnSparse - entries per tracklet
299 // species, p, tracklet position, number of PID tracklets, number of slices (non 0), number of clusters, electron likelihood,
302 Int_t nbins[7] = {5, 20, 6, 7, 9, 45, 100};
303 Double_t min[7] = {-0.5, 0.5, -0.5, -0.5, -0.5, -0.5, 0.};
304 Double_t max[7] = {4.5, 10, 5.5, 6.5, 8.5, 179.5, 1.};
305 TString htitle = "TRD tracklets; VO identified species; p (GeV/c); tracklet position; No. PID tacklets; No. slices; No. clusters; El. likelihood";
306 fOutput->CreateTHnSparse("hTRDtracklets",htitle,7, nbins, min, max);
307 fOutput->BinLogAxis("hTRDtracklets", 1);
312 //__________________________________________
313 void AliHFEpidQA::Process(){
319 AliError("V0pid not available! Forgotten to initialize?");
323 AliError("fESDpid not initialized, I am leaving this!");
327 // to be udpated to AOD save mdoe
329 AliError("AliVEvent not available, returning");
333 if(fMC) fV0pidMC->SetMCEvent(fMC);
334 if(fMC) fV0pid->SetMCEvent(fMC);
336 fV0pid->Process(fEvent);
337 TObjArray *hfeelectrons = fV0pid->GetListOfElectrons();
338 TObjArray *hfepionsK0 = fV0pid->GetListOfPionsK0();
339 TObjArray *hfepionsL = fV0pid->GetListOfPionsL();
340 TObjArray *hfeprotons = fV0pid->GetListOfProtons();
342 // Get Track list for normal purpose
343 TObjArray *electrons = MakeTrackList(hfeelectrons);
344 TObjArray *pionsK0 = MakeTrackList(hfepionsK0);
345 TObjArray *pionsL = MakeTrackList(hfepionsL);
346 TObjArray *protons = MakeTrackList(hfeprotons);
347 TObjArray *cleanElectrons = MakeCleanListElectrons(hfeelectrons);
350 fV0pidMC->Process(electrons, AliHFEV0cuts::kRecoElectron);
351 fV0pidMC->Process(pionsK0, AliHFEV0cuts::kRecoPionK0);
352 fV0pidMC->Process(pionsL, AliHFEV0cuts::kRecoPionL);
353 fV0pidMC->Process(protons, AliHFEV0cuts::kRecoProton);
356 AliDebug(2, Form("Number of Electrons : %d", electrons->GetEntries()));
357 AliDebug(2, Form("Number of K0 Pions : %d", pionsK0->GetEntries()));
358 AliDebug(2, Form("Number of Lambda Pions : %d", pionsL->GetEntries()));
359 AliDebug(2, Form("Number of Protons : %d", protons->GetEntries()));
361 AliDebug(2, "MC Information available. Doing Purity checks...");
362 // Calculate the purity of the clean samples using MC
363 MakePurity(electrons, AliHFEV0cuts::kRecoElectron);
364 MakePurity(pionsK0, AliHFEV0cuts::kRecoPionK0);
365 MakePurity(pionsL, AliHFEV0cuts::kRecoPionL);
366 MakePurity(protons, AliHFEV0cuts::kRecoProton);
369 // some event wise checks
372 // Make Illumination Plot
373 FillIllumination(electrons, AliHFEV0cuts::kRecoElectron);
374 FillIllumination(pionsK0, AliHFEV0cuts::kRecoPionK0);
375 FillIllumination(pionsL, AliHFEV0cuts::kRecoPionL);
376 FillIllumination(protons, AliHFEV0cuts::kRecoProton);
378 // Now we can do studies on the PID itself
379 // For TRD use the TRD PID QA object
381 if(fIsPbPb) centrality=GetCentrality(fEvent);
382 if(fIsppMultiBin) centrality=GetMultiplicityITS(fEvent);
383 fTRDpidQA->SetCentralityBin(centrality);
384 fTRDpidQA->ProcessTracks(cleanElectrons, AliPID::kElectron);
385 fTRDpidQA->ProcessTracks(pionsK0, AliPID::kPion);
386 fTRDpidQA->ProcessTracks(pionsL, AliPID::kPion);
387 fTRDpidQA->ProcessTracks(protons, AliPID::kProton);
389 // Monitor TRD PID Response
390 /*TestTRDResponse(cleanElectrons, AliHFEV0cuts::kRecoElectron);
391 TestTRDResponse(pionsK0, AliHFEV0cuts::kRecoPionK0);
392 TestTRDResponse(pionsL, AliHFEV0cuts::kRecoPionL);
393 TestTRDResponse(protons, AliHFEV0cuts::kRecoProton);*/
395 FillElectronLikelihoods(electrons, AliHFEV0cuts::kRecoElectron);
396 FillElectronLikelihoods(pionsK0, AliHFEV0cuts::kRecoPionK0);
397 FillElectronLikelihoods(pionsL, AliHFEV0cuts::kRecoPionL);
398 FillElectronLikelihoods(protons, AliHFEV0cuts::kRecoProton);
400 FillPIDresponse(electrons, AliHFEV0cuts::kRecoElectron);
401 FillPIDresponse(pionsK0, AliHFEV0cuts::kRecoPionK0);
402 FillPIDresponse(pionsL, AliHFEV0cuts::kRecoPionL);
403 FillPIDresponse(protons, AliHFEV0cuts::kRecoProton);
405 // check the tender V0s
406 CheckTenderV0pid(electrons, AliHFEV0cuts::kRecoElectron);
407 CheckTenderV0pid(pionsK0, AliHFEV0cuts::kRecoPionK0);
408 CheckTenderV0pid(pionsL, AliHFEV0cuts::kRecoPionL);
409 CheckTenderV0pid(protons, AliHFEV0cuts::kRecoProton);
411 // Analysis done, flush the containers
418 delete cleanElectrons;
421 //__________________________________________
422 void AliHFEpidQA::FillIllumination(const TObjArray * const tracks, Int_t species){
424 // Fill Illumination Plot
426 THnSparseF *hIllumination = dynamic_cast<THnSparseF *>(fOutput->Get("hIllumination"));
427 if(!hIllumination) return;
429 Double_t quantities[6]; memset(quantities, 0, sizeof(Double_t) *6);
430 TIter trackIter(tracks);
432 quantities[0] = species;
433 TObject *o = NULL; AliESDtrack *esdtrack = NULL;
434 while((o = trackIter())){
435 if(!TString(o->IsA()->GetName()).CompareTo("AliESDtrack")){
436 // work on local copy in order to not spoil others
437 esdtrack = new AliESDtrack(*(static_cast<AliESDtrack *>(o)));
438 if(!esdtrack) continue;
439 } else if(!TString(o->IsA()->GetName()).CompareTo("AliAODrack")){
440 // Bad hack: Fill ESD track with AOD information
441 esdtrack = new AliESDtrack(static_cast<AliAODTrack *>(o));
442 if(!esdtrack) continue;
448 // Propagate to the entrance of the TRD
449 esdtrack->PropagateTo(300, fEvent->GetMagneticField());
450 quantities[1] = esdtrack->Pt();
451 quantities[2] = esdtrack->Eta();
452 quantities[3] = esdtrack->Phi();
453 quantities[4] = (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) ? 1 : 0;
454 quantities[5] = (esdtrack->GetStatus() & AliESDtrack::kTRDout) ? 1. : 0.;
455 hIllumination->Fill(quantities);
460 //__________________________________________
461 void AliHFEpidQA::FillTPCinfo(AliESDtrack *const esdtrack, Int_t species){
463 // Fill TPC Cluster Plots
465 THnSparseF *hTPCclusters = dynamic_cast<THnSparseF *>(fOutput->Get("hTPCclusters"));
466 if(!hTPCclusters) return;
468 Double_t quantities[6]; memset(quantities, 0, sizeof(Double_t) *6);
470 Double_t pidProbs[5];
471 const Int_t typePID[5] = {0, 2, 2, 3, 4};
474 quantities[0] = species;
477 esdtrack->GetTPCpid(pidProbs);
479 quantities[1] = esdtrack->P();
480 quantities[2] = esdtrack->GetTPCNcls();
481 quantities[3] = pidProbs[0];
482 quantities[4] = fESDpid->NumberOfSigmasTPC(esdtrack,(AliPID::EParticleType)typePID[species]);
483 quantities[5] = esdtrack->GetTPCsignal();
484 hTPCclusters->Fill(quantities);
488 //__________________________________________
489 void AliHFEpidQA::TestTRDResponse(const TObjArray * const tracks, Int_t species){
491 // Test PID Response function of the TRD
493 Int_t effInt[6] = {70, 75, 80, 85, 90, 95};
494 const char *sname[5] = {"Electron", "PionK0", "PionL", "Kaon", "Proton"};
495 AliESDtrack *track = NULL;
496 TIter trackIter(tracks);
497 Float_t momenta[6], p;
499 Double_t probs[5], likeEl;
500 while((track = static_cast<AliESDtrack *>(trackIter()))){
501 if(track->GetTRDntrackletsPID() < 4) continue;
503 // calculate momentum at TRD level
504 memset(momenta, 0, sizeof(Float_t) * 6); nmomenta = 0;
505 for(Int_t ily = 0; ily < 6; ily++){
506 if(track->GetTRDmomentum(ily) > 0.01) momenta[nmomenta++] = track->GetTRDmomentum(ily);
508 p = TMath::Mean(nmomenta, momenta);
510 // Get Electron likelihood
511 track->GetTRDpid(probs);
512 likeEl = probs[0]/(probs[0] + probs[2]);
514 for(Int_t ieff = 0; ieff < 6; ieff++){
515 if(fESDpid->IdentifiedAsElectronTRD(track, static_cast<Double_t>(effInt[ieff])/100.)) fOutput->Fill(Form("hTRD_likeSel_%s_%dtls_%deff", sname[species], static_cast<Int_t>(track->GetTRDntrackletsPID()), effInt[ieff]), p, likeEl);
516 else fOutput->Fill(Form("hTRD_likeRej_%s_%dtls_%deff", sname[species], static_cast<Int_t>(track->GetTRDntrackletsPID()), effInt[ieff]), p, likeEl);
521 //__________________________________________
522 void AliHFEpidQA::MakePurity(const TObjArray *tracks, Int_t species){
524 // Fill the QA histos for a given species
527 AliDebug(3, Form("Doing Purity checks for species %d", species));
528 Int_t pdg = GetPDG(species);
532 case AliHFEV0cuts::kRecoElectron:
533 hname = "purityElectron";
535 case AliHFEV0cuts::kRecoPionK0:
536 hname = "purityPionK0";
538 case AliHFEV0cuts::kRecoPionL:
539 hname = "purityPionL";
541 case AliHFEV0cuts::kRecoProton:
542 hname = "purityProton";
544 default: // non investigated species
545 AliDebug(3, "Species not investigated");
549 AliDebug(3, Form("Number of tracks: %d", tracks->GetEntries()));
550 TIter trackIter(tracks);
551 AliVParticle *recTrack = NULL, *mcTrack = NULL;
552 while((recTrack = dynamic_cast<AliVParticle *>(trackIter()))){
553 Int_t label = recTrack->GetLabel();
554 AliDebug(4, Form("MC Label %d", label));
555 mcTrack =fMC->GetTrack(TMath::Abs(label));
557 AliDebug(4, "MC track not available");
558 continue; // we don't know
563 if(!TString(mcTrack->IsA()->GetName()).CompareTo("AliMCParticle")){
565 AliMCParticle *mcp = dynamic_cast<AliMCParticle *>(mcTrack);
567 trackPdg = TMath::Abs(mcp->Particle()->GetPdgCode());
570 AliAODMCParticle *aodmcp = dynamic_cast<AliAODMCParticle *>(mcTrack);
571 if(!aodmcp) continue;
572 trackPdg = TMath::Abs(aodmcp->GetPdgCode());
574 if(trackPdg == pdg) // Correct identification
576 fOutput->Fill(hname, 0., recTrack->Pt());
578 else // Wrong identification
579 fOutput->Fill(hname, 1., recTrack->Pt());
583 //__________________________________________
584 void AliHFEpidQA::FillElectronLikelihoods(const TObjArray * const particles, Int_t species){
586 // Fill electron Likelihoods for the ITS, TPC and TOF
587 // Required for the calculation of the electron efficiency,
588 // pion and proton efficiency and the thresholds
591 const TString detname[4] = {"ITS", "TPC", "TRD", "TOF"};
595 case AliHFEV0cuts::kRecoElectron:
596 specname = "Electron";
598 case AliHFEV0cuts::kRecoPionK0:
601 case AliHFEV0cuts::kRecoPionL:
604 case AliHFEV0cuts::kRecoProton:
608 AliDebug(2, Form("Species %d not investigated", species));
611 AliVParticle *recTrack = NULL;
612 // mcTrack =fMC->GetTrack(TMath::Abs(label));
614 // AliDebug(4, "MC track not available");
615 // continue; // we don't know
618 TIter trackIter(particles);
620 Double_t quantities[2];
621 Double_t pidProbs[5];
623 while((recTrack = dynamic_cast<AliVParticle *>(trackIter()))){
624 if(!TString(recTrack->IsA()->GetName()).CompareTo("AliESDtrack")){
626 AliESDtrack *esdTrack = dynamic_cast<AliESDtrack *>(recTrack);
627 if(!esdTrack) continue;
628 status = esdTrack->GetStatus();
630 //TPC momentum and likelihoods
632 pTPC = esdTrack->GetInnerParam() ? esdTrack->GetInnerParam()->P() : esdTrack->P();
633 Bool_t mcFound = kFALSE;
635 Int_t label = esdTrack->GetLabel();
636 AliMCParticle *mcTrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(label));
637 Int_t pdg = GetPDG(species);
640 trackPdg = TMath::Abs(mcTrack->Particle()->GetPdgCode());
642 if(pdg == trackPdg) mcFound = kTRUE;
644 quantities[0] = pTPC;
645 Bool_t detFlagSet = kFALSE;
646 for(Int_t idet = 0; idet < 4; idet++){
647 TString histname, histnameMC;
648 histname = "h" + detname[idet] + "_El_like_" + specname;
649 histnameMC = "h" + detname[idet] + "_El_like_MC_" + specname;
652 case kITS: esdTrack->GetITSpid(pidProbs);
653 detFlagSet = status & AliESDtrack::kITSpid;
655 case kTPC: esdTrack->GetTPCpid(pidProbs);
656 detFlagSet = status & AliESDtrack::kTPCpid;
658 case kTRD: esdTrack->GetTRDpid(pidProbs);
659 detFlagSet = status & AliESDtrack::kTRDpid;
661 case kTOF: esdTrack->GetTOFpid(pidProbs);
662 detFlagSet = status & AliESDtrack::kTOFpid;
665 quantities[1] = pidProbs[AliPID::kElectron];
666 // in case of TRD require 6 PID tracklets
667 if(kTRD == idet && esdTrack->GetTRDntrackletsPID() != 6) continue;
669 fOutput->Fill(histname, quantities[0], quantities[1]);
671 fOutput->Fill(histnameMC, quantities[0], quantities[1]);
680 //__________________________________________
681 void AliHFEpidQA::FillPIDresponse(const TObjArray * const particles, Int_t species){
683 // Fill the PID response of different detectors to V0 daughter particles
685 TString hname, hname2, hname3;
687 const TString typeName[5] = {"Electron", "PionK0", "PionL", "Kaon", "Proton"};
688 const Int_t typePID[5] = {0, 2, 2, 3, 4};
692 // 0) species, 1) momentum, 2) DCA xy, 3) DCA z
694 // 5) TPC Ncls 6) TPC signal 7) TPC nSigma,
695 // 8) TRD Ntrk, 9) TRD Ncls, 10) TRD dEdx,
700 Int_t run = fEvent->GetRunNumber();
702 AliVParticle *recTrack = NULL;
703 TIter trackIter(particles);
704 while((recTrack = dynamic_cast<AliVParticle *>(trackIter()))){
705 //for(Int_t i=0; i<12; ++i) data[i] = -99.;
707 if(!TString(recTrack->IsA()->GetName()).CompareTo("AliESDtrack")){
709 AliESDtrack *esdTrack = dynamic_cast<AliESDtrack *>(recTrack);
710 if(!esdTrack) continue;
712 // for the PID THnSparse
714 //data[1] = esdTrack->P();
715 Float_t impactR = -1.;
716 Float_t impactZ = -1.;
717 esdTrack->GetImpactParameters(impactR, impactZ);
720 //data[11] = 0; // initialize the TOF pid cut on elecgrons to false
721 // use ONLY tracks with PID flag TRUE
723 status = esdTrack->GetStatus();
729 fOutput->Fill("hITS_kFlags", 5., species);
730 if(status & AliESDtrack::kITSin) fOutput->Fill("hITS_kFlags", 1., species);
731 if(status & AliESDtrack::kITSout) fOutput->Fill("hITS_kFlags", 2., species);
732 if(status & AliESDtrack::kITSrefit) fOutput->Fill("hITS_kFlags", 3., species);
733 if(status & AliESDtrack::kITSpid) fOutput->Fill("hITS_kFlags", 4., species);
735 fOutput->Fill("hTPC_kFlags", 5., species);
736 if(status & AliESDtrack::kTPCin) fOutput->Fill("hTPC_kFlags", 1., species);
737 if(status & AliESDtrack::kTPCout) fOutput->Fill("hTPC_kFlags", 2., species);
738 if(status & AliESDtrack::kTPCrefit) fOutput->Fill("hTPC_kFlags", 3., species);
739 if(status & AliESDtrack::kTPCpid) fOutput->Fill("hTPC_kFlags", 4., species);
741 fOutput->Fill("hTRD_kFlags", 5., species);
742 if(status & AliESDtrack::kTRDin) fOutput->Fill("hTRD_kFlags", 1., species);
743 if(status & AliESDtrack::kTRDout) fOutput->Fill("hTRD_kFlags", 2., species);
744 if(status & AliESDtrack::kTRDrefit) fOutput->Fill("hTRD_kFlags", 3., species);
745 if(status & AliESDtrack::kTRDpid) fOutput->Fill("hTRD_kFlags", 4., species);
747 fOutput->Fill("hTOF_kFlags", 5., species);
748 if(status & AliESDtrack::kTOFin) fOutput->Fill("hTOF_kFlags", 1., species);
749 if(status & AliESDtrack::kTOFout) fOutput->Fill("hTOF_kFlags", 2., species);
750 if(status & AliESDtrack::kTOFrefit) fOutput->Fill("hTOF_kFlags", 3., species);
751 if(status & AliESDtrack::kTOFpid) fOutput->Fill("hTOF_kFlags", 4., species);
757 if(status & AliESDtrack::kITSpid){
758 Double_t p = esdTrack->P();
761 //Double_t itsSignal = esdTrack->GetITSsignal();
764 Double_t dEdxSamples[4];
765 esdTrack->GetITSdEdxSamples(dEdxSamples);
767 Double_t dEdxSum = 0.;
768 hname = "hITS_dEdx_" + typeName[species];
769 for(Int_t i=0; i<4; ++i){
770 if(dEdxSamples[i] > 0){
772 fOutput->Fill(hname, i+1, p, dEdxSamples[i]);
773 dEdxSum += dEdxSamples[i];
776 if(4 == nSamples)fOutput->Fill(hname, 0, p, dEdxSum);
777 fOutput->Fill("hITS_dEdx_nSamples", nSamples);
779 Double_t signal = esdTrack->GetITSsignal();
780 hname = "hITS_Signal_" + typeName[species];
781 fOutput->Fill(hname, p, signal);
784 // ITS number of signas
785 Double_t nsigma = fESDpid->NumberOfSigmasITS(esdTrack,(AliPID::EParticleType)typePID[species]);
786 hname = "hITS_nSigma_" + typeName[species];
787 fOutput->Fill(hname, p, nsigma);
789 Double_t itsPID[5] = {-1, -1, -1, -1, -1};
790 esdTrack->GetITSpid(itsPID);
791 Int_t ix = GetMaxPID(itsPID);
792 hname = "hITS_PID_p_" + typeName[species];
793 fOutput->Fill(hname, p, ix);
799 if(status & AliESDtrack::kTPCpid){
800 // Make TPC clusters Plot
801 //data[5] = esdTrack->GetTPCNcls();
802 FillTPCinfo(esdTrack, species);
804 Double_t p = esdTrack->GetInnerParam() ? esdTrack->GetInnerParam()->P() : esdTrack->P();
806 Double_t dEdx = esdTrack->GetTPCsignal();
807 hname = "hTPC_dEdx_" + typeName[species];
808 fOutput->Fill(hname, p, dEdx);
811 //TPC number of sigmas
812 Double_t nsigma = fESDpid->NumberOfSigmasTPC(esdTrack,(AliPID::EParticleType)typePID[species]);
813 hname = "hTPC_nSigma_" + typeName[species];
814 fOutput->Fill(hname, p, nsigma);
818 hname = "hTPC_PID_p_" + typeName[species];
819 Double_t tpcPID[5] = {-1, -1, -1, -1, -1};
820 esdTrack->GetTPCpid(tpcPID);
821 Int_t ix = GetMaxPID(tpcPID);
822 fOutput->Fill(hname, p, ix);
829 if(status & AliESDtrack::kTRDpid){
830 Double_t p = esdTrack->GetOuterParam() ? esdTrack->GetOuterParam()->P() : esdTrack->P();
832 // TRD number of tracklets
833 Int_t ntrk = esdTrack->GetTRDntrackletsPID();
834 hname = "hTRD_trk_" + typeName[species];
835 fOutput->Fill(hname, p, ntrk);
839 hname = "hTRD_PID_p_" + typeName[species];
840 Double_t trdPID[5] = {-1., -1., -1., -1., -1};
841 esdTrack->GetTRDpid(trdPID);
842 Int_t ix = GetMaxPID(trdPID);
843 fOutput->Fill(hname, p, ix);
845 Int_t ncls = esdTrack->GetTRDncls();
846 hname = "hTRD_cls_" + typeName[species];
847 fOutput->Fill(hname, p, ncls);
850 // TRD - per tracklet - dEdx per, likelihood
851 hname = "hTRD_Nslc_" + typeName[species];
852 hname2 = "hTRD_slc_" + typeName[species];
853 hname3 = "hTRD_dEdx_" + typeName[species];
854 Int_t nSlices = esdTrack->GetNumberOfTRDslices();
855 Double_t sumTot = 0.;
857 for(Int_t l=0; l< 6; ++l){
858 Double_t trkData[7] = {-1.,-1, -1, -1, -1, -1, -1};
859 trkData[0] = species;
866 for(Int_t s=0; s<nSlices; ++s){
867 Double_t slice = esdTrack->GetTRDslice(l, s);
871 fOutput->Fill(hname2, p, s);
876 fOutput->Fill(hname, p, not0);
877 fOutput->Fill(hname3, p, sum);
882 // lkelihoods per layer
883 if(not0Tot > 0 && fNNref){
884 Double_t likelihoods[5] = {-1., -1., -1., -1., -1};
885 TRDlikeTracklet(l, esdTrack, likelihoods);
886 trkData[6] = likelihoods[0];
887 //printf(" -D: species: %i, P; %f : %f, s: %f\n", species, p, likelihoods[0], s);
889 if(not0Tot) fOutput->Fill("hTRDtracklets", trkData);
891 // average dEx per number of tracklets
893 // data[10] = sumTot / not0Tot;
900 if(status & AliESDtrack::kTOFpid){
901 Double_t p = esdTrack->GetOuterParam() ? esdTrack->GetOuterParam()->P() : esdTrack->P();
902 Double_t t0 = fESDpid->GetTOFResponse().GetStartTime(esdTrack->P());
905 hname = "hTOF_beta_" + typeName[species];
906 Float_t beta = TOFbeta(esdTrack);
907 fOutput->Fill(hname, p, beta);
908 fOutput->Fill("hTOF_beta_all", p, beta);
910 Double_t nsigma = fESDpid->NumberOfSigmasTOF(esdTrack,(AliPID::EParticleType)typePID[species]);
911 hname = "hTOF_nSigma_" + typeName[species];
912 fOutput->Fill(hname, p, nsigma);
913 //if(beta > 0.97 && beta < 1.03){
918 hname = "hTOF_PID_p_" + typeName[species];
919 Double_t tofPID[5] = {-1., -1., -1., -1., -1};
920 esdTrack->GetTOFpid(tofPID);
921 Int_t ix = GetMaxPID(tofPID);
922 fOutput->Fill(hname, p, ix);
925 // - distribution of (time - t0 - pion_time)
927 esdTrack->GetIntegratedTimes(times);
928 Double_t tItrackL = esdTrack->GetIntegratedLength();
929 Double_t tTOFsignal = esdTrack->GetTOFsignal();
930 Double_t dT = tTOFsignal - t0 - times[2];
931 fOutput->Fill("hTOF_qa_TmT0mT", run*1.0, dT);
932 fOutput->Fill("hTOF_qa_length", run*1.0, tItrackL);
936 // temporary - the PIDsparse needs rebuilding
937 //fOutput->Fill("PIDsparse", data);
939 // AOD - comming soon
943 }// .. tracks in TObjArray
947 //__________________________________________
948 void AliHFEpidQA:: CheckEvent(){
950 // check some event variables
953 // check the T0 as a function of run number (less than one bin per number
954 Double_t t0 = fESDpid->GetTOFResponse().GetTimeZero();
955 Int_t run = fEvent->GetRunNumber();
956 Double_t data[2] = {run*1.0, t0*1000.};
957 fOutput->Fill("hEvent_T0", data);
961 //__________________________________________
962 TList *AliHFEpidQA::GetOutput(){
964 // Getter for Output histograms
966 return fOutput->GetList();
969 //__________________________________________
970 TList *AliHFEpidQA::GetV0pidQA(){
972 // Getter for V0 PID QA histograms
974 return fV0pid->GetListOfQAhistograms();
977 //__________________________________________
978 TList *AliHFEpidQA::GetV0pidMC(){
980 // Getter for V0 PID QA histograms
983 return fV0pidMC->GetListOfQAhistograms();
987 //__________________________________________
988 void AliHFEpidQA::RecalculateTRDpid(AliESDtrack * /*track*/, Double_t * /*pidProbs*/) const{
989 // fTRDpidResponse->MakePID(track);
990 // track->GetTRDpid(pidProbs);
993 //__________________________________________
994 void AliHFEpidQA::RecalculateTRDpid(AliAODTrack * /*track*/, Double_t * /*pidProbs*/) const{
995 // fTRDpidResponse->MakePID(track, pidProbs);
997 //___________________________________________________________________
998 Float_t AliHFEpidQA::TOFbeta(const AliESDtrack * const track) const {
999 // computes the TOF beta
1000 Double_t l = track->GetIntegratedLength(); // cm
1001 Double_t t = track->GetTOFsignal();
1002 Double_t t0 = fESDpid->GetTOFResponse().GetTimeZero(); // ps
1004 //printf("-D: l: %f, t: %f, t0: %f\n", l, t, t0);
1006 if(l < 360. || l > 800.) return 0.;
1007 if(t <= 0.) return 0.;
1008 if(t0 >999990.0) return 0.;
1011 t -= t0; // subtract the T0
1013 l *= 0.01; // cm ->m
1014 t *= 1e-12; //ps -> s
1018 Float_t beta = v / TMath::C();
1022 //____________________________________________
1023 Int_t AliHFEpidQA::GetMaxPID(const Double_t *pidProbs) const {
1025 // return the index of maximal PID probability
1029 for(Int_t i=0; i<5; ++i){
1030 if(pidProbs[i] > tmp){
1037 //_____________________________________________
1038 Int_t AliHFEpidQA::GetPDG(Int_t species){
1040 // return the PDG particle code
1046 case AliHFEV0cuts::kRecoElectron:
1047 pdg = TMath::Abs(kElectron);
1049 case AliHFEV0cuts::kRecoPionK0:
1050 pdg = TMath::Abs(kPiPlus);
1052 case AliHFEV0cuts::kRecoPionL:
1053 pdg = TMath::Abs(kPiPlus);
1055 case AliHFEV0cuts::kRecoProton:
1056 pdg = TMath::Abs(kProton);
1058 default: // non investigated species
1059 AliDebug(3, "Species not recognised");
1067 //_____________________________________________
1068 TObjArray * AliHFEpidQA::MakeTrackList(const TObjArray *tracks) const {
1070 // convert list of AliHFEV0Info into a list of AliVParticle
1072 TObjArray *output = new TObjArray;
1073 TIter trackInfos(tracks);
1074 AliHFEV0info *trackInfo = NULL;
1075 while((trackInfo = dynamic_cast<AliHFEV0info *>(trackInfos())))
1076 output->Add(trackInfo->GetTrack());
1081 //_____________________________________________
1082 TObjArray * AliHFEpidQA::MakeCleanListElectrons(const TObjArray *electrons) const {
1084 // Cleanup electron sample using TPC PID
1085 // PID requirement will allways be implemented to the pair
1088 TObjArray *tracks = new TObjArray;
1089 TIter candidates(electrons);
1092 AliHFEV0info *hfetrack;
1093 // const Double_t kSigmaTight = 1.;
1094 // const Double_t kSigmaLoose = 4.;
1095 const Double_t kSigmaTight = 2.;
1096 const Double_t kSigmaLoose = 2.;
1097 const Double_t shift = 0.57;
1098 if((esd = dynamic_cast<AliESDEvent *>(fEvent))){
1099 AliESDtrack *track = NULL, *partnerTrack = NULL;
1100 while((hfetrack = dynamic_cast<AliHFEV0info *>(candidates()))){
1101 track = dynamic_cast<AliESDtrack *>(hfetrack->GetTrack());
1102 if(!track) continue;
1103 partnerTrack = esd->GetTrack(hfetrack->GetPartnerID());
1104 Double_t nSigmaTrack = TMath::Abs(fESDpid->NumberOfSigmasTPC(track, AliPID::kElectron) - shift);
1105 Double_t nSigmaPartner = TMath::Abs(fESDpid->NumberOfSigmasTPC(partnerTrack, AliPID::kElectron) - shift);
1106 if((nSigmaTrack < kSigmaTight && nSigmaPartner < kSigmaLoose) || (nSigmaTrack < kSigmaLoose && nSigmaPartner < kSigmaTight))
1111 aod = dynamic_cast<AliAODEvent *>(fEvent);
1112 if(!aod) return NULL;
1113 //AliAODTrack *track = NULL, *partnerTrack = NULL;
1114 while((hfetrack = dynamic_cast<AliHFEV0info *>(candidates()))){
1115 if(!hfetrack) continue;
1116 //track = dynamic_cast<AliAODTrack *>(hfetrack->GetTrack());
1117 //partnerTrack = aod->GetTrack(hfetrack->GetPartnerID());
1118 // will be coming soon
1123 //___________________________________________________________
1124 void AliHFEpidQA::CheckTenderV0pid(const TObjArray * const particles, Int_t species){
1127 // retrieve the status bits from the TObject used to flag
1128 // the V0 daughter tracks and return the found PID
1129 // 0 - electron, 1 - pion, 2 - proton
1132 const Int_t id[5] = {0, 1, 1, -1, 2}; // convert AliHFEpid to simple 0, 1, 2
1134 AliVParticle *recTrack = NULL;
1135 TIter trackIter(particles);
1136 while((recTrack = dynamic_cast<AliVParticle *>(trackIter()))){
1137 if(!TString(recTrack->IsA()->GetName()).CompareTo("AliESDtrack")){
1139 AliESDtrack *track = dynamic_cast<AliESDtrack *>(recTrack);
1140 if(!track) continue;
1141 Int_t tPID = GetTenderV0pid(track);
1142 fOutput->Fill("h_tender_check_01", id[species]*1.0, tPID);
1143 fOutput->Fill("h_tender_check_02", id[species], id[species]*1.0, tPID);
1147 }//.. iterate of tracks
1149 //___________________________________________________________
1150 Int_t AliHFEpidQA::GetTenderV0pid(AliESDtrack * const track){
1152 // retrieve the PID nformation stored in the status flags by the train tender
1163 if(track->TestBit(2<<14)){
1167 if(track->TestBit(2<<15)){
1171 if(track->TestBit(2<<16)){
1177 AliWarning("V0 track labeled multiple times by the V0 tender");
1184 //___________________________________________________________
1185 Double_t AliHFEpidQA::TRDlikeTracklet(Int_t layer, AliESDtrack * const track, Double_t *likelihood){
1187 // compute the TRD electron likelihoods for 1 tracklet
1188 // based on teh AliTRDpidRecalculator in train/until/tender
1189 // returns sum of the likelihoods (which should be 1)
1192 const Double_t cScaleGain = 1./ 16000.;
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
1195 if(!track) return kFALSE;
1196 Float_t p = track->GetTRDmomentum(layer); // momentum for a tracklet in the ESDtrack
1197 if(p < 0) return kFALSE;
1199 Int_t mombin = TRDmomBin(p); // momentum bin
1200 if(mombin < 0) return -1.;
1201 Float_t dEdxTRDsum = 0; // dEdxTRDsum for checking if tracklet is available
1202 Float_t dEdxTRD[8]; // dEdx for a tracklet in the ESD slices
1203 Double_t ddEdxTRD[8]; // dEdx as Double_t for TMultiLayerPerceptron::Evaluate()
1205 Double_t prob[AliPID::kSPECIES]; // probabilities for all species in all layers
1207 for(Int_t is = 0; is < AliPID::kSPECIES; is++){
1208 likelihood[is] = 0.2; // init probabilities
1214 for(Int_t islice = 0; islice<8; islice++){
1215 dEdxTRD[islice]=0.; // init dE/dx
1216 ddEdxTRD[islice]=0.; // init dE/dx
1217 dEdxTRD[islice]=track->GetTRDslice(layer,islice); // get the dE/dx values
1218 dEdxTRDsum += dEdxTRD[islice];
1219 ddEdxTRD[islice]=(Double_t)dEdxTRD[islice]*cScaleGain; // rescale dE/dx
1222 for(Int_t is = 0; is < AliPID::kSPECIES; is++){
1223 Double_t probloc1, probloc2;
1224 if(mombin == 0 && mombin < pBins[0]){ // calculate PID for p > 0.6 GeV/c
1225 prob[is] = fNet[mombin]->Evaluate(is, ddEdxTRD);
1227 else if(mombin == 10 && mombin >= pBins[10]){ // calculate PID for p >= 10.0 GeV/c
1228 prob[is] = fNet[mombin]->Evaluate(is, ddEdxTRD);
1230 else{ // standard calculation
1231 Int_t mombin1 = 0, mombin2 = 0; // lower and upper momentum bin
1232 if(p < pBins[mombin]) {mombin1 = mombin -1; mombin2 = mombin;}
1233 if(p >= pBins[mombin]) {mombin1 = mombin; mombin2 = mombin+1;}
1234 probloc1 = fNet[mombin1]->Evaluate(is, ddEdxTRD);
1235 probloc2 = fNet[mombin2]->Evaluate(is, ddEdxTRD);
1236 // weighting of the probability with regard to the track momentum
1237 prob[is] = probloc1 + (probloc2-probloc1)*(p-pBins[mombin1])/(pBins[mombin2]-pBins[mombin1]);
1239 likelihood[is] = prob[is];
1240 sum += likelihood[is];
1245 //__________________________________________________________________________
1246 Int_t AliHFEpidQA::TRDmomBin(Double_t p) const {
1248 // compute the momentum bin position
1251 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
1253 Int_t pBin1 = -1; // check bin1
1254 Int_t pBin2 = -1; // check bin2
1256 if(p < 0) return -1; // return -1 if momentum < 0
1257 if(p < pBins[0]) return 0; // smallest momentum bin
1258 if(p >= pBins[10]) return 10; // largest momentum bin
1261 // calculate momentum bin for non extremal momenta
1262 for(Int_t iMomBin = 1; iMomBin < 11; iMomBin++){
1263 if(p < pBins[iMomBin]){
1264 pBin1 = iMomBin - 1;
1270 if(p - pBins[pBin1] >= pBins[pBin2] - p){
1283 //___________________________________________________
1284 Int_t AliHFEpidQA::GetCentrality(AliVEvent* const fInputEvent){
1286 // Recover the centrality of the event from ESD or AOD
1293 AliCentrality *centrality = fInputEvent->GetCentrality();
1294 Float_t fCentralityFtemp = -1;
1295 fCentralityFtemp = centrality->GetCentralityPercentile("V0M");
1296 Float_t centralityLimits[12] = {0.,5.,10., 20., 30., 40., 50., 60.,70.,80., 90., 100.};
1297 for(Int_t ibin = 0; ibin < 11; ibin++){
1298 if(fCentralityFtemp >= centralityLimits[ibin] && fCentralityFtemp < centralityLimits[ibin+1]){
1304 AliDebug(2, Form("Centrality class %d\n", bin));
1311 //___________________________________________________
1312 Int_t AliHFEpidQA::GetMultiplicityITS(AliVEvent* const fInputEvent) const
1315 // Definition of the Multiplicity according to the JPSI group (F. Kramer)
1317 Int_t nTracklets = 0;
1319 Double_t etaRange = 1.6;
1322 nTracklets = ((AliESDEvent*) fInputEvent)->GetMultiplicity()->GetNumberOfTracklets();
1323 for (Int_t nn = 0; nn < nTracklets; nn++) {
1324 Double_t eta = ((AliESDEvent*)fInputEvent)->GetMultiplicity()->GetEta(nn);
1325 if (TMath::Abs(eta) < etaRange) nAcc++;
1328 Int_t itsMultiplicity = nAcc;
1330 Int_t multiplicityLimits[8] = {0, 1, 9, 17, 25, 36, 60, 500};
1331 for(Int_t ibin = 0; ibin < 7; ibin++){
1332 if(itsMultiplicity >= multiplicityLimits[ibin] && itsMultiplicity < multiplicityLimits[ibin + 1]){