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