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