PWG3/hfe/AliHFEpidQA.cxx

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