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