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