PWG1/TRD/AliTRDcheckPID.cxx

   1 //////////////////////////////////////////////////////
   2 //
   3 // PID performance checker of the TRD
   4 //
   5 // Performs checks of ESD information for TRD-PID and recalculate PID based on OCDB information
   6 // Also provides performance plots on detector based on the PID information - for the moment only
   7 // MC source is used but V0 is also possible. Here is a list of detector performance checks
   8 //   - Integrated dE/dx per chamber
   9 //   - <PH> as function of time bin and local radial position
  10 //   - number of clusters/tracklet
  11 //   - number of tracklets/track
  12 //
  13 // Author : Alex Wilk <wilka@uni-muenster.de>
  14 //          Alex Bercuci <A.Bercuci@gsi.de>
  15 //          Markus Fasel <M.Fasel@gsi.de>
  16 //
  17 ///////////////////////////////////////////////////////
  18
  19 #include "TAxis.h"
  20 #include "TROOT.h"
  21 #include "TPDGCode.h"
  22 #include "TCanvas.h"
  23 #include "TF1.h"
  24 #include "TH1F.h"
  25 #include "TH1D.h"
  26 #include "TH2F.h"
  27 #include "TProfile.h"
  28 #include "TProfile2D.h"
  29 #include "TGraph.h"
  30 #include "TGraphErrors.h"
  31 #include "TLegend.h"
  32
  33 #include <TClonesArray.h>
  34 #include <TObjArray.h>
  35 #include <TList.h>
  36
  37 #include "AliESDEvent.h"
  38 #include "AliESDInputHandler.h"
  39 #include "AliTrackReference.h"
  40
  41 #include "AliAnalysisTask.h"
  42
  43 #include "AliTRDtrackerV1.h"
  44 #include "AliTRDtrackV1.h"
  45 #include "AliTRDcluster.h"
  46 #include "AliTRDReconstructor.h"
  47 #include "AliCDBManager.h"
  48 #include "AliTRDpidUtil.h"
  49
  50 #include "Cal/AliTRDCalPID.h"
  51 #include "Cal/AliTRDCalPIDNN.h"
  52 #include "AliTRDcheckPID.h"
  53 #include "AliTRDinfoGen.h"
  54 #include "AliAnalysisManager.h"
  55 #include "info/AliTRDtrackInfo.h"
  56 #include "info/AliTRDpidInfo.h"
  57 #include "info/AliTRDv0Info.h"
  58
  59 Char_t const * AliTRDcheckPID::fgMethod[3] = {"LQ", "NN", "ESD"};
  60
  61 ClassImp(AliTRDcheckPID)
  62
  63 //________________________________________________________________________
  64 AliTRDcheckPID::AliTRDcheckPID()
  65   :AliTRDrecoTask()
  66   ,fUtil(NULL)
  67   ,fGraph(NULL)
  68   ,fPID(NULL)
  69   ,fV0s(NULL)
  70   ,fMomentumAxis(NULL)
  71   ,fMinNTracklets(AliTRDgeometry::kNlayer)
  72   ,fMaxNTracklets(AliTRDgeometry::kNlayer)
  73  {
  74   //
  75   // Default constructor
  76   //
  77   SetNameTitle("TRDcheckPID", "Check TRD PID");
  78   LocalInit();
  79 }
  80
  81 //________________________________________________________________________
  82 AliTRDcheckPID::AliTRDcheckPID(char* name )
  83   :AliTRDrecoTask(name, "Check TRD PID")
  84   ,fUtil(NULL)
  85   ,fGraph(NULL)
  86   ,fPID(NULL)
  87   ,fV0s(NULL)
  88   ,fMomentumAxis(NULL)
  89   ,fMinNTracklets(AliTRDgeometry::kNlayer)
  90   ,fMaxNTracklets(AliTRDgeometry::kNlayer)
  91  {
  92   //
  93   // Default constructor
  94   //
  95
  96   LocalInit();
  97   InitFunctorList();
  98
  99   DefineInput(3, TObjArray::Class());  // v0 list
 100   DefineOutput(2, TObjArray::Class()); // pid info list
 101 }
 102
 103
 104 //________________________________________________________________________
 105 void AliTRDcheckPID::LocalInit()
 106 {
 107 // Initialize working data
 108
 109   // Initialize momentum axis with default values
 110   Double_t defaultMomenta[AliTRDCalPID::kNMom+1];
 111   for(Int_t imom = 0; imom < AliTRDCalPID::kNMom+1; imom++)
 112     defaultMomenta[imom] = AliTRDCalPID::GetMomentumBinning(imom);
 113   SetMomentumBinning(AliTRDCalPID::kNMom, defaultMomenta);
 114
 115   memset(fEfficiency, 0, AliPID::kSPECIES*sizeof(TObjArray*));
 116
 117   fUtil = new AliTRDpidUtil();
 118 }
 119
 120 //________________________________________________________________________
 121 AliTRDcheckPID::~AliTRDcheckPID()
 122 {
 123   if (AliAnalysisManager::GetAnalysisManager()->IsProofMode()) return;
 124   if(fPID){fPID->Delete(); delete fPID;}
 125   if(fGraph){fGraph->Delete(); delete fGraph;}
 126   if(fUtil) delete fUtil;
 127 }
 128
 129
 130 //________________________________________________________________________
 131 void AliTRDcheckPID::UserCreateOutputObjects()
 132 {
 133   // Create histograms
 134   // Called once
 135
 136   AliTRDrecoTask::UserCreateOutputObjects();
 137   fPID = new TObjArray();
 138   fPID->SetOwner(kTRUE);
 139   PostData(2, fPID);
 140 }
 141
 142 //________________________________________________________
 143 void AliTRDcheckPID::UserExec(Option_t *opt)
 144 {
 145   //
 146   // Execution part
 147   //
 148
 149   fV0s = dynamic_cast<TObjArray *>(GetInputData(3));
 150   fPID->Delete();
 151
 152   AliTRDrecoTask::UserExec(opt);
 153 }
 154
 155
 156 //_______________________________________________________
 157 TObjArray * AliTRDcheckPID::Histos(){
 158
 159   //
 160   // Create QA histograms
 161   //
 162   if(fContainer) return fContainer;
 163
 164   Int_t xBins = AliPID::kSPECIES*fMomentumAxis->GetNbins();
 165   fContainer = new TObjArray(); fContainer->Expand(kNPlots); fContainer->SetOwner(kTRUE);
 166
 167   const Float_t epsilon = 1./(2*(AliTRDpidUtil::kBins-1));     // get nice histos with bin center at 0 and 1
 168   TH1 *h = NULL;
 169
 170   const Int_t kNmethodsPID=Int_t(sizeof(fgMethod)/sizeof(Char_t*));
 171   // histos with posterior probabilities for all particle species
 172   for(Int_t is=AliPID::kSPECIES; is--;){
 173     fEfficiency[is] = new TObjArray(kNmethodsPID);
 174     fEfficiency[is]->SetOwner();
 175     fEfficiency[is]->SetName(Form("Eff_%s", AliPID::ParticleShortName(is)));
 176     fContainer->AddAt(fEfficiency[is], is?kEfficiencyMu+is-1:kEfficiency);
 177     for(Int_t im=kNmethodsPID; im--;){
 178       if(!(h = (TH2F*)gROOT->FindObject(Form("PID_%s_%s", fgMethod[im], AliPID::ParticleShortName(is))))) {
 179         h = new TH2F(Form("PID_%s_%s", fgMethod[im], AliPID::ParticleShortName(is)), "", xBins, -0.5, xBins - 0.5,
 180           AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon);
 181       } else h->Reset();
 182       fEfficiency[is]->AddAt(h, im);
 183     }
 184   }
 185
 186   // histos of the dE/dx distribution for all 5 particle species and 11 momenta
 187   if(!(h = (TH2F*)gROOT->FindObject("dEdx"))){
 188     h = new TH2F("dEdx", "",
 189       xBins, -0.5, xBins - 0.5,
 190       200, 0, 15);
 191 //       200, 0, 10000);
 192   } else h->Reset();
 193   fContainer->AddAt(h, kdEdx);
 194
 195   // histos of the dE/dx slices for all 5 particle species and 11 momenta
 196   if(!(h = (TH2F*)gROOT->FindObject("dEdxSlice"))){
 197     h = new TH2F("dEdxSlice", "",
 198       xBins*AliTRDpidUtil::kLQslices, -0.5, xBins*AliTRDpidUtil::kLQslices - 0.5,
 199       200, 0, 5000);
 200   } else h->Reset();
 201   fContainer->AddAt(h, kdEdxSlice);
 202
 203   // histos of the pulse height distribution for all 5 particle species and 11 momenta
 204   TObjArray *fPH = new TObjArray(2);
 205   fPH->SetOwner(); fPH->SetName("PH");
 206   fContainer->AddAt(fPH, kPH);
 207   if(!(h = (TProfile2D*)gROOT->FindObject("PHT"))){
 208     h = new TProfile2D("PHT", "<PH>(tb);p*species;tb [100*ns];entries",
 209       xBins, -0.5, xBins - 0.5,
 210       AliTRDseedV1::kNtb, -0.5, AliTRDseedV1::kNtb - 0.5);
 211   } else h->Reset();
 212   fPH->AddAt(h, 0);
 213   if(!(h = (TProfile2D*)gROOT->FindObject("PHX"))){
 214     h = new TProfile2D("PHX", "<PH>(x);p*species;x_{drift} [cm];entries",
 215       xBins, -0.5, xBins - 0.5,
 216       40, 0., 4.5);
 217   } else h->Reset();
 218   fPH->AddAt(h, 1);
 219
 220   // histos of the number of clusters distribution for all 5 particle species and 11 momenta
 221   if(!(h = (TH2F*)gROOT->FindObject("NClus"))){
 222     h = new TH2F("NClus", "",
 223       xBins, -0.5, xBins - 0.5,
 224       50, -0.5, 49.5);
 225   } else h->Reset();
 226   fContainer->AddAt(h, kNClus);
 227
 228
 229   // momentum distributions - absolute and in momentum bins
 230   if(!(h = (TH1F*)gROOT->FindObject("hMom"))){
 231     h = new TH1F("hMom", "momentum distribution", fMomentumAxis->GetNbins(), fMomentumAxis->GetXmin(), fMomentumAxis->GetXmax());
 232   } else h->Reset();
 233   fContainer->AddAt(h, kMomentum);
 234
 235   if(!(h = (TH1F*)gROOT->FindObject("hMomBin"))){
 236     h = new TH1F("hMomBin", "momentum distribution in momentum bins", fMomentumAxis->GetNbins(), fMomentumAxis->GetXmin(), fMomentumAxis->GetXmax());
 237   } else h->Reset();
 238   fContainer->AddAt(h, kMomentumBin);
 239
 240   // Number of tracklets per track
 241   if(!(h = (TH2F*)gROOT->FindObject("nTracklets"))){
 242     h = new TH2F("nTracklets", "",
 243       xBins, -0.5, xBins - 0.5,
 244       AliTRDgeometry::kNlayer, 0.5, AliTRDgeometry::kNlayer+.5);
 245   } else h->Reset();
 246   fContainer->AddAt(h, kNTracklets);
 247
 248   // V0 performance
 249   if(!(h = (TH1F*)gROOT->FindObject("nV0"))){
 250     h = new TH1F("nV0", "V0s/track;v0/track;entries",
 251       6, -0.5, 5.5);
 252   } else h->Reset();
 253   fContainer->AddAt(h, kV0);
 254
 255   // dQ/dl for 1D-Likelihood
 256   if(!(h = (TH1F *)gROOT->FindObject("dQdl"))){
 257     h = new TH2F("dQdl", "dQ/dl per layer;p*species;dQ/dl [a.u.]", xBins, -0.5, xBins - 0.5, 800, 0., 40000.);
 258   } else h->Reset();
 259   fContainer->AddAt(h, kdQdl);
 260
 261   return fContainer;
 262 }
 263
 264
 265 //________________________________________________________________________
 266 Bool_t AliTRDcheckPID::CheckTrackQuality(const AliTRDtrackV1* track) const
 267 {
 268   //
 269   // Check if the track is ok for PID
 270   //
 271
 272   Int_t ntracklets = track->GetNumberOfTracklets();
 273   if(ntracklets >= fMinNTracklets && ntracklets <= fMaxNTracklets) return 1;
 274 //   if(!fESD)
 275 //     return 0;
 276
 277   return 0;
 278 }
 279
 280 //________________________________________________________________________
 281 Int_t AliTRDcheckPID::CalcPDG(AliTRDtrackV1* track)
 282 {
 283 // Documentation to come
 284
 285  /* track -> SetReconstructor(AliTRDinfoGen::Reconstructor());
 286   (const_cast<AliTRDrecoParam*>(AliTRDinfoGen::Reconstructor()->GetRecoParam()))->SetPIDNeuralNetwork();
 287   track -> CookPID();*/
 288
 289   if(track -> GetPID(AliPID::kElectron) > track -> GetPID(AliPID::kMuon) + track -> GetPID(AliPID::kPion)  + track -> GetPID(AliPID::kKaon) + track -> GetPID(AliPID::kProton)){
 290     return kElectron;
 291   }
 292   else if(track -> GetPID(kProton) > track -> GetPID(AliPID::kPion)  && track -> GetPID(AliPID::kProton) > track -> GetPID(AliPID::kKaon)  && track -> GetPID(AliPID::kProton) > track -> GetPID(AliPID::kMuon)){
 293     return kProton;
 294   }
 295   else if(track -> GetPID(AliPID::kKaon) > track -> GetPID(AliPID::kMuon)  && track -> GetPID(AliPID::kKaon) > track -> GetPID(AliPID::kPion)){
 296     return kKPlus;
 297   }
 298   else if(track -> GetPID(AliPID::kMuon) > track -> GetPID(AliPID::kPion)){
 299     return kMuonPlus;
 300   }
 301   else{
 302     return kPiPlus;
 303   }
 304 }
 305
 306
 307 //_______________________________________________________
 308 TH1 *AliTRDcheckPID::PlotLQ(const AliTRDtrackV1 *track)
 309 {
 310   //
 311   // Plot the probabilities for electrons using 2-dim LQ
 312   //
 313
 314   if(!fkESD){
 315     AliDebug(2, "No ESD info available.");
 316     return NULL;
 317   }
 318
 319   if(track) fkTrack = track;
 320   if(!fkTrack){
 321     AliDebug(2, "No Track defined.");
 322     return NULL;
 323   }
 324
 325   ULong_t status;
 326   status = fkESD -> GetStatus();
 327   if(!(status&AliESDtrack::kTRDin)) return NULL;
 328
 329   if(!CheckTrackQuality(fkTrack)) return NULL;
 330
 331   AliTRDtrackV1 cTrack(*fkTrack);
 332   cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
 333
 334   Int_t pdg = 0;
 335   Float_t momentum = 0.;
 336
 337   if(fkMC){
 338     if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
 339     pdg = fkMC->GetPDG();
 340   } else{
 341     //AliWarning("No MC info available!");
 342     momentum = cTrack.GetMomentum(0);
 343     pdg = CalcPDG(&cTrack);
 344   }
 345   if(!IsInRange(momentum)) return NULL;
 346
 347   (const_cast<AliTRDrecoParam*>(AliTRDinfoGen::Reconstructor()->GetRecoParam()))->SetPIDNeuralNetwork(kFALSE);
 348   cTrack.CookPID();
 349   if(cTrack.GetNumberOfTrackletsPID() < fMinNTracklets) return NULL;
 350   Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
 351
 352   TH2F *hPIDLQ(NULL);
 353   TObjArray *eff(NULL);
 354   for(Int_t is=AliPID::kSPECIES; is--;){
 355     if(!(eff = dynamic_cast<TObjArray *>(fContainer->At(is?kEfficiencyMu+is-1:kEfficiency)))){
 356       AliWarning("No Histogram List defined.");
 357       return NULL;
 358     }
 359     if(!(hPIDLQ = dynamic_cast<TH2F *>(eff->At(AliTRDpidUtil::kLQ)))){
 360       AliWarning("No Histogram defined.");
 361       return NULL;
 362     }
 363
 364     hPIDLQ -> Fill(FindBin(species, momentum), cTrack.GetPID(is));
 365   }
 366   return hPIDLQ;
 367 }
 368
 369
 370
 371
 372 //_______________________________________________________
 373 TH1 *AliTRDcheckPID::PlotNN(const AliTRDtrackV1 *track)
 374 {
 375   //
 376   // Plot the probabilities for electrons using 2-dim LQ
 377   //
 378
 379   if(!fkESD){
 380     AliDebug(2, "No ESD info available.");
 381     return NULL;
 382   }
 383
 384   if(track) fkTrack = track;
 385   if(!fkTrack){
 386     AliDebug(2, "No Track defined.");
 387     return NULL;
 388   }
 389
 390   ULong_t status;
 391   status = fkESD -> GetStatus();
 392   if(!(status&AliESDtrack::kTRDin)) return NULL;
 393
 394   if(!CheckTrackQuality(fkTrack)) return NULL;
 395
 396   AliTRDtrackV1 cTrack(*fkTrack);
 397   cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
 398
 399   Int_t pdg = 0;
 400   Float_t momentum = 0.;
 401   if(fkMC){
 402     if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
 403     pdg = fkMC->GetPDG();
 404   } else {
 405     //AliWarning("No MC info available!");
 406     momentum = cTrack.GetMomentum(0);
 407     pdg = CalcPDG(&cTrack);
 408   }
 409   if(!IsInRange(momentum)) return NULL;
 410
 411   (const_cast<AliTRDrecoParam*>(AliTRDinfoGen::Reconstructor()->GetRecoParam()))->SetPIDNeuralNetwork();
 412   cTrack.CookPID();
 413   if(cTrack.GetNumberOfTrackletsPID() < fMinNTracklets) return NULL;
 414
 415   Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
 416
 417
 418   TH2F *hPIDNN(NULL);
 419   TObjArray *eff(NULL);
 420   for(Int_t is=AliPID::kSPECIES; is--;){
 421     if(!(eff = dynamic_cast<TObjArray *>(fContainer->At(is?kEfficiencyMu+is-1:kEfficiency)))){
 422       AliWarning("No Histogram List defined.");
 423       return NULL;
 424     }
 425     if(!(hPIDNN = dynamic_cast<TH2F *>(eff->At(AliTRDpidUtil::kNN)))){
 426       AliWarning("No Histogram defined.");
 427       return NULL;
 428     }
 429
 430     hPIDNN->Fill(FindBin(species, momentum), cTrack.GetPID(is));
 431   }
 432   return hPIDNN;
 433 }
 434
 435
 436
 437 //_______________________________________________________
 438 TH1 *AliTRDcheckPID::PlotESD(const AliTRDtrackV1 *track)
 439 {
 440   //
 441   // Plot the probabilities for electrons using 2-dim LQ
 442   //
 443
 444   if(!fkESD){
 445     AliDebug(2, "No ESD info available.");
 446     return NULL;
 447   }
 448
 449   if(track) fkTrack = track;
 450   if(!fkTrack){
 451     AliDebug(2, "No Track defined.");
 452     return NULL;
 453   }
 454
 455   ULong_t status;
 456   status = fkESD -> GetStatus();
 457   if(!(status&AliESDtrack::kTRDin)) return NULL;
 458
 459   if(!CheckTrackQuality(fkTrack)) return NULL;
 460   if(fkTrack->GetNumberOfTrackletsPID() < fMinNTracklets) return NULL;
 461
 462   Int_t pdg = 0;
 463   Float_t momentum = 0.;
 464   if(fkMC){
 465     if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
 466     pdg = fkMC->GetPDG();
 467   } else {
 468     //AliWarning("No MC info available!");
 469     AliTRDtrackV1 cTrack(*fkTrack);
 470     cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
 471     momentum = cTrack.GetMomentum(0);
 472     pdg = CalcPDG(&cTrack);
 473   }
 474   if(!IsInRange(momentum)) return NULL;
 475
 476   const Double32_t *pidESD = fkESD->GetResponseIter();
 477   Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
 478
 479   TH2F *hPID(NULL);
 480   TObjArray *eff(NULL);
 481   for(Int_t is=AliPID::kSPECIES; is--;){
 482     if(!(eff = dynamic_cast<TObjArray *>(fContainer->At(is?kEfficiencyMu+is-1:kEfficiency)))){
 483       AliWarning("No Histogram List defined.");
 484       return NULL;
 485     }
 486     if(!(hPID = dynamic_cast<TH2F *>(eff->At(AliTRDpidUtil::kESD)))){
 487       AliWarning("No Histogram defined.");
 488       return NULL;
 489     }
 490
 491     hPID->Fill(FindBin(species, momentum), pidESD[is]);
 492   }
 493   return hPID;
 494 }
 495
 496
 497
 498 //_______________________________________________________
 499 TH1 *AliTRDcheckPID::PlotdQdl(const AliTRDtrackV1 *track){
 500   //
 501   // Plot the total charge for the 1D Likelihood method
 502   //
 503   if(track) fkTrack = track;
 504   if(!fkTrack){
 505     AliDebug(2, "No Track defined");
 506     return NULL;
 507   }
 508   TH2 *hdQdl = dynamic_cast<TH2F *>(fContainer->At(kdQdl));
 509   if(!hdQdl){
 510     AliWarning("No Histogram defined");
 511     return NULL;
 512   }
 513
 514   if(!CheckTrackQuality(fkTrack)) return NULL;
 515
 516   Int_t pdg = 0;
 517   Float_t momentum = 0.;
 518   AliTRDtrackV1 cTrack(*fkTrack);
 519   if(fkMC){
 520     if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
 521     pdg = fkMC->GetPDG();
 522   } else {
 523     //AliWarning("No MC info available!");
 524     momentum = cTrack.GetMomentum(0);
 525     pdg = CalcPDG(&cTrack);
 526   }
 527   if(!IsInRange(momentum)) return NULL;
 528
 529   // Init exchange container
 530   Int_t s(AliTRDpidUtil::Pdg2Pid(pdg));
 531   Int_t ibin = FindBin(s, momentum);
 532
 533   AliTRDseedV1 *tracklet = NULL;
 534   for(Int_t iseed = 0; iseed < 6; iseed++){
 535     if(!((tracklet = fkTrack->GetTracklet(iseed)) && tracklet->IsOK())) continue;
 536     hdQdl->Fill(ibin, tracklet->GetdQdl());
 537   }
 538   return hdQdl;
 539 }
 540
 541 //_______________________________________________________
 542 TH1 *AliTRDcheckPID::PlotdEdx(const AliTRDtrackV1 *track)
 543 {
 544   //
 545   // Plot the probabilities for electrons using 2-dim LQ
 546   //
 547
 548   if(track) fkTrack = track;
 549   if(!fkTrack){
 550     AliDebug(2, "No Track defined.");
 551     return NULL;
 552   }
 553
 554   if(!CheckTrackQuality(fkTrack)) return NULL;
 555
 556   TH2F *hdEdx;
 557   if(!(hdEdx = dynamic_cast<TH2F *>(fContainer->At(kdEdx)))){
 558     AliWarning("No Histogram defined.");
 559     return NULL;
 560   }
 561
 562   AliTRDtrackV1 cTrack(*fkTrack);
 563   cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
 564   Int_t pdg = 0;
 565   Float_t momentum = 0.;
 566   if(fkMC){
 567     if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
 568     pdg = fkMC->GetPDG();
 569   } else {
 570     //AliWarning("No MC info available!");
 571     momentum = cTrack.GetMomentum(0);
 572     pdg = CalcPDG(&cTrack);
 573   }
 574   if(!IsInRange(momentum)) return NULL;
 575
 576   // Init exchange container
 577   Int_t s(AliTRDpidUtil::Pdg2Pid(pdg));
 578   AliTRDpidInfo *pid = new AliTRDpidInfo(s);
 579
 580   (const_cast<AliTRDrecoParam*>(AliTRDinfoGen::Reconstructor()->GetRecoParam()))->SetPIDNeuralNetwork(kTRUE);
 581
 582   Float_t sumdEdx(0.);
 583   Int_t iBin = FindBin(s, momentum);
 584   AliTRDseedV1 *tracklet = NULL;
 585   for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++){
 586     tracklet = cTrack.GetTracklet(ily);
 587     if(!tracklet) continue;
 588     tracklet -> CookdEdx(AliTRDpidUtil::kNNslices);
 589
 590     // fill exchange container
 591     pid->PushBack(tracklet->GetPlane(),
 592                   AliTRDpidUtil::GetMomentumBin(tracklet->GetMomentum()), tracklet->GetdEdx());
 593
 594     sumdEdx = 0.;
 595     for(Int_t i = AliTRDpidUtil::kNNslices; i--;) sumdEdx += tracklet->GetdEdx()[i];
 596     sumdEdx /= AliTRDCalPIDNN::kMLPscale;
 597     hdEdx -> Fill(iBin, sumdEdx);
 598   }
 599   fPID->Add(pid);
 600
 601   return hdEdx;
 602 }
 603
 604
 605 //_______________________________________________________
 606 TH1 *AliTRDcheckPID::PlotdEdxSlice(const AliTRDtrackV1 *track)
 607 {
 608   //
 609   // Plot the probabilities for electrons using 2-dim LQ
 610   //
 611
 612   if(track) fkTrack = track;
 613   if(!fkTrack){
 614     AliDebug(2, "No Track defined.");
 615     return NULL;
 616   }
 617
 618   if(!CheckTrackQuality(fkTrack)) return NULL;
 619
 620   TH2F *hdEdxSlice;
 621   if(!(hdEdxSlice = dynamic_cast<TH2F *>(fContainer->At(kdEdxSlice)))){
 622     AliWarning("No Histogram defined.");
 623     return NULL;
 624   }
 625
 626   AliTRDtrackV1 cTrack(*fkTrack);
 627   cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
 628   Int_t pdg = 0;
 629   Float_t momentum = 0.;
 630   if(fkMC){
 631     if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
 632     pdg = fkMC->GetPDG();
 633   } else {
 634     //AliWarning("No MC info available!");
 635     momentum = cTrack.GetMomentum(0);
 636     pdg = CalcPDG(&cTrack);
 637   }
 638   if(!IsInRange(momentum)) return NULL;
 639
 640   (const_cast<AliTRDrecoParam*>(AliTRDinfoGen::Reconstructor()->GetRecoParam()))->SetPIDNeuralNetwork(kFALSE);
 641   Int_t iMomBin = fMomentumAxis->FindBin(momentum);
 642   Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
 643   Float_t *fdEdx;
 644   AliTRDseedV1 *tracklet = NULL;
 645   for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
 646     tracklet = cTrack.GetTracklet(iChamb);
 647     if(!tracklet) continue;
 648     tracklet -> CookdEdx(AliTRDpidUtil::kLQslices);
 649     fdEdx = const_cast<Float_t *>(tracklet->GetdEdx());
 650     for(Int_t iSlice = 0; iSlice < AliTRDpidUtil::kLQslices; iSlice++){
 651       hdEdxSlice -> Fill(species * fMomentumAxis->GetNbins() * AliTRDpidUtil::kLQslices + (iMomBin-1) * AliTRDpidUtil::kLQslices + iSlice, fdEdx[iSlice]);
 652     }
 653   }
 654
 655   return hdEdxSlice;
 656 }
 657
 658
 659 //_______________________________________________________
 660 TH1 *AliTRDcheckPID::PlotPH(const AliTRDtrackV1 *track)
 661 {
 662   //
 663   // Plot the probabilities for electrons using 2-dim LQ
 664   //
 665
 666   if(track) fkTrack = track;
 667   if(!fkTrack){
 668     AliDebug(2, "No Track defined.");
 669     return NULL;
 670   }
 671
 672   if(!CheckTrackQuality(fkTrack)) return NULL;
 673
 674   TObjArray *arr = NULL;
 675   TProfile2D *hPHX, *hPHT;
 676   if(!(arr = dynamic_cast<TObjArray *>(fContainer->At(kPH)))){
 677     AliWarning("No Histogram defined.");
 678     return NULL;
 679   }
 680   hPHT = (TProfile2D*)arr->At(0);
 681   hPHX = (TProfile2D*)arr->At(1);
 682
 683   Int_t pdg = 0;
 684   Float_t momentum = 0.;
 685   if(fkMC){
 686     if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
 687     pdg = fkMC->GetPDG();
 688   } else {
 689     //AliWarning("No MC info available!");
 690     AliTRDtrackV1 cTrack(*fkTrack);
 691     cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
 692     momentum = cTrack.GetMomentum(0);
 693     pdg = CalcPDG(&cTrack);
 694   }
 695   if(!IsInRange(momentum)) return NULL;;
 696
 697   AliTRDseedV1 *tracklet = NULL;
 698   AliTRDcluster *cluster = NULL;
 699   Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
 700   Int_t iBin = FindBin(species, momentum);
 701   for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
 702     tracklet = fkTrack->GetTracklet(iChamb);
 703     if(!tracklet) continue;
 704     Float_t x0 = tracklet->GetX0();
 705     for(Int_t ic = 0; ic < AliTRDseedV1::kNclusters; ic++){
 706       if(!(cluster = tracklet->GetClusters(ic))) continue;
 707       hPHT -> Fill(iBin, cluster->GetLocalTimeBin(), TMath::Abs(cluster->GetQ()));
 708       if(ic<AliTRDseedV1::kNtb) hPHX -> Fill(iBin, x0 - cluster->GetX(), tracklet->GetdQdl(ic));
 709     }
 710   }
 711   return hPHT;
 712 }
 713
 714
 715 //_______________________________________________________
 716 TH1 *AliTRDcheckPID::PlotNClus(const AliTRDtrackV1 *track)
 717 {
 718   //
 719   // Plot the probabilities for electrons using 2-dim LQ
 720   //
 721
 722   if(track) fkTrack = track;
 723   if(!fkTrack){
 724     AliDebug(2, "No Track defined.");
 725     return NULL;
 726   }
 727
 728   if(!CheckTrackQuality(fkTrack)) return NULL;
 729
 730   TH2F *hNClus;
 731   if(!(hNClus = dynamic_cast<TH2F *>(fContainer->At(kNClus)))){
 732     AliWarning("No Histogram defined.");
 733     return NULL;
 734   }
 735
 736
 737   Int_t pdg = 0;
 738   Float_t momentum = 0.;
 739   if(fkMC){
 740     if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
 741     pdg = fkMC->GetPDG();
 742   } else {
 743     //AliWarning("No MC info available!");
 744     AliTRDtrackV1 cTrack(*fkTrack);
 745     cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
 746     momentum = cTrack.GetMomentum(0);
 747     pdg = CalcPDG(&cTrack);
 748   }
 749   if(!IsInRange(momentum)) return NULL;
 750
 751   Int_t species = AliTRDpidUtil::AliTRDpidUtil::Pdg2Pid(pdg);
 752   Int_t iBin = FindBin(species, momentum);
 753   AliTRDseedV1 *tracklet = NULL;
 754   for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
 755     tracklet = fkTrack->GetTracklet(iChamb);
 756     if(!tracklet) continue;
 757     hNClus -> Fill(iBin, tracklet->GetN());
 758   }
 759
 760   return hNClus;
 761 }
 762
 763 //_______________________________________________________
 764 TH1 *AliTRDcheckPID::PlotNTracklets(const AliTRDtrackV1 *track)
 765 {
 766   //
 767   // Plot the probabilities for electrons using 2-dim LQ
 768   //
 769
 770   if(track) fkTrack = track;
 771   if(!fkTrack){
 772     AliDebug(2, "No Track defined.");
 773     return NULL;
 774   }
 775
 776   TH2F *hTracklets;
 777   if(!(hTracklets = dynamic_cast<TH2F *>(fContainer->At(kNTracklets)))){
 778     AliWarning("No Histogram defined.");
 779     return NULL;
 780   }
 781
 782   AliTRDtrackV1 cTrack(*fkTrack);
 783   cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
 784   Int_t pdg = 0;
 785   Float_t momentum = 0.;
 786   if(fkMC){
 787     if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
 788     pdg = fkMC->GetPDG();
 789   } else {
 790     //AliWarning("No MC info available!");
 791     momentum = cTrack.GetMomentum();
 792     pdg = CalcPDG(&cTrack);
 793   }
 794   Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
 795   if(!IsInRange(momentum)) return NULL;
 796
 797   Int_t iBin = FindBin(species, momentum);
 798   hTracklets -> Fill(iBin, cTrack.GetNumberOfTracklets());
 799   return hTracklets;
 800 }
 801
 802 //_______________________________________________________
 803 TH1 *AliTRDcheckPID::PlotMom(const AliTRDtrackV1 *track)
 804 {
 805   //
 806   // Plot the probabilities for electrons using 2-dim LQ
 807   //
 808
 809   if(track) fkTrack = track;
 810   if(!fkTrack){
 811     AliDebug(2, "No Track defined.");
 812     return NULL;
 813   }
 814
 815   if(!CheckTrackQuality(fkTrack)) return NULL;
 816
 817   TH1F *hMom;
 818   if(!(hMom = dynamic_cast<TH1F *>(fContainer->At(kMomentum)))){
 819     AliWarning("No Histogram defined.");
 820     return NULL;
 821   }
 822
 823
 824   Int_t pdg = 0;
 825   Float_t momentum = 0.;
 826   if(fkMC){
 827     if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
 828     pdg = fkMC->GetPDG();
 829   } else {
 830     //AliWarning("No MC info available!");
 831     AliTRDtrackV1 cTrack(*fkTrack);
 832     cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
 833     momentum = cTrack.GetMomentum(0);
 834     pdg = CalcPDG(&cTrack);
 835   }
 836   if(IsInRange(momentum)) hMom -> Fill(momentum);
 837   return hMom;
 838 }
 839
 840
 841 //_______________________________________________________
 842 TH1 *AliTRDcheckPID::PlotMomBin(const AliTRDtrackV1 *track)
 843 {
 844   //
 845   // Plot the probabilities for electrons using 2-dim LQ
 846   //
 847
 848   if(track) fkTrack = track;
 849   if(!fkTrack){
 850     AliDebug(2, "No Track defined.");
 851     return NULL;
 852   }
 853
 854   if(!CheckTrackQuality(fkTrack)) return NULL;
 855
 856   TH1F *hMomBin;
 857   if(!(hMomBin = dynamic_cast<TH1F *>(fContainer->At(kMomentumBin)))){
 858     AliWarning("No Histogram defined.");
 859     return NULL;
 860   }
 861
 862
 863   Int_t pdg = 0;
 864   Float_t momentum = 0.;
 865
 866   if(fkMC){
 867     if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
 868     pdg = fkMC->GetPDG();
 869   } else {
 870     //AliWarning("No MC info available!");
 871     AliTRDtrackV1 cTrack(*fkTrack);
 872     cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
 873     momentum = cTrack.GetMomentum(0);
 874   }
 875   if(IsInRange(momentum)) hMomBin -> Fill(fMomentumAxis->FindBin(momentum));
 876   return hMomBin;
 877 }
 878
 879 //_______________________________________________________
 880 TH1 *AliTRDcheckPID::PlotV0(const AliTRDtrackV1 *track)
 881 {
 882   //
 883   // Plot the V0 performance against MC
 884   //
 885
 886   if(track) fkTrack = track;
 887   if(!fkTrack){
 888     AliDebug(2, "No Track defined.");
 889     return NULL;
 890   }
 891   if(!fkESD->HasV0()) return NULL;
 892   if(!HasMCdata()){
 893     AliDebug(1, "No MC defined.");
 894     return NULL;
 895   }
 896   if(!fContainer){
 897     AliWarning("No output container defined.");
 898     return NULL;
 899   }
 900   AliDebug(2, Form("TRACK[%d] species[%s][%d]\n", fkESD->GetId(), fkMC->GetPID()>=0?AliPID::ParticleShortName(fkMC->GetPID()):"none", fkMC->GetPDG()));
 901
 902   TH1 *h(NULL);
 903   if(!(h = dynamic_cast<TH1F*>(fContainer->At(kV0)))) return NULL;
 904   Int_t sgn(0), n(0); AliTRDv0Info *v0(NULL);
 905   for(Int_t iv0(fV0s->GetEntriesFast()); iv0--;){
 906     if(!(v0=(AliTRDv0Info*)fV0s->At(iv0))) continue;
 907     if(!(sgn = v0->HasTrack(fkESD->GetId()))) continue;
 908     //for(Int_t is=AliPID::kSPECIES; is--;) v0->GetPID(is, track);
 909     //v0->Print();
 910     n++;
 911     //break;
 912   }
 913   h->Fill(n);
 914   return h;
 915 }
 916
 917 //________________________________________________________
 918 Bool_t AliTRDcheckPID::GetRefFigure(Int_t ifig)
 919 {
 920 // Steering function to retrieve performance plots
 921
 922   Bool_t kFIRST = kTRUE;
 923   TGraphErrors *g = NULL;
 924   TAxis *ax = NULL;
 925   TObjArray *arr = NULL;
 926   TH1 *h1 = NULL, *h=NULL;
 927   TH2 *h2 = NULL;
 928   TList *content = NULL;
 929   switch(ifig){
 930   case kEfficiency:{
 931     gPad->Divide(2, 1, 1.e-5, 1.e-5);
 932     TList *l=gPad->GetListOfPrimitives();
 933     TVirtualPad *pad = ((TVirtualPad*)l->At(0));pad->cd();
 934     pad->SetMargin(0.1, 0.01, 0.1, 0.01);
 935
 936     TLegend *legEff = new TLegend(.64, .84, .98, .98);
 937     legEff->SetBorderSize(1);legEff->SetTextSize(0.03255879);
 938     legEff->SetFillColor(0);
 939     h=new TH1S("hEff", "", 1, .5, 11.);
 940     h->SetLineColor(1);h->SetLineWidth(1);
 941     ax = h->GetXaxis();
 942     ax->SetTitle("p [GeV/c]");
 943     ax->SetRangeUser(.5, 11.);
 944     ax->SetMoreLogLabels();
 945     ax = h->GetYaxis();
 946     ax->SetTitle("#epsilon_{#pi} [%]");
 947     ax->CenterTitle();
 948     ax->SetRangeUser(1.e-2, 10.);
 949     h->Draw();
 950     content = (TList *)fGraph->FindObject(Form("Eff_%s", AliTRDCalPID::GetPartName(AliPID::kPion)));
 951     if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kLQ))) break;
 952     if(!g->GetN()) break;
 953     legEff->SetHeader("PID Method [PION]");
 954     g->Draw("pc"); legEff->AddEntry(g, "LQ 2D", "pl");
 955     if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kNN))) break;
 956     g->Draw("pc"); legEff->AddEntry(g, "NN", "pl");
 957     if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kESD))) break;
 958     g->Draw("p"); legEff->AddEntry(g, "ESD", "pl");
 959     legEff->Draw();
 960     gPad->SetLogy();
 961     gPad->SetLogx();
 962     gPad->SetGridy();
 963     gPad->SetGridx();
 964
 965
 966     pad = ((TVirtualPad*)l->At(1));pad->cd();
 967     pad->SetMargin(0.1, 0.01, 0.1, 0.01);
 968     h=new TH1S("hThr", "", 1, .5, 11.);
 969     h->SetLineColor(1);h->SetLineWidth(1);
 970     ax = h->GetXaxis();
 971     ax->SetTitle("p [GeV/c]");
 972     ax->SetMoreLogLabels();
 973     ax = h->GetYaxis();
 974     ax->SetTitle("Threshold [%]");
 975     ax->SetRangeUser(5.e-2, 1.);
 976     h->Draw();
 977     content = (TList *)fGraph->FindObject("Thres");
 978     if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kLQ))) break;
 979     if(!g->GetN()) break;
 980     g->Draw("pc");
 981     if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kNN))) break;
 982     g->Draw("pc");
 983     if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kESD))) break;
 984     g->Draw("p");
 985     gPad->SetLogx();
 986     gPad->SetGridy();
 987     gPad->SetGridx();
 988     return kTRUE;
 989   }
 990   case kEfficiencyKa:{
 991     gPad->Divide(2, 1, 1.e-5, 1.e-5);
 992     TList *l=gPad->GetListOfPrimitives();
 993     TVirtualPad *pad = ((TVirtualPad*)l->At(0));pad->cd();
 994     pad->SetMargin(0.1, 0.01, 0.1, 0.01);
 995
 996     TLegend *legEff = new TLegend(.64, .84, .98, .98);
 997     legEff->SetBorderSize(1);legEff->SetTextSize(0.03255879);
 998     legEff->SetFillColor(0);
 999     h = (TH1S*)gROOT->FindObject("hEff");
1000     h=(TH1S*)h->Clone("hEff_K");
1001     h->SetYTitle("#epsilon_{K} [%]");
1002     h->GetYaxis()->SetRangeUser(1.e-2, 1.e2);
1003     h->Draw();
1004     content = (TList *)fGraph->FindObject(Form("Eff_%s", AliTRDCalPID::GetPartName(AliPID::kKaon)));
1005     if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kLQ))) break;
1006     if(!g->GetN()) break;
1007     legEff->SetHeader("PID Method [KAON]");
1008     g->Draw("pc"); legEff->AddEntry(g, "LQ 2D", "pl");
1009     if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kNN))) break;
1010     g->Draw("pc"); legEff->AddEntry(g, "NN", "pl");
1011     if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kESD))) break;
1012     g->Draw("p"); legEff->AddEntry(g, "ESD", "pl");
1013     legEff->Draw();
1014     gPad->SetLogy();
1015     gPad->SetLogx();
1016     gPad->SetGridy();
1017     gPad->SetGridx();
1018
1019     TLegend *legEff2 = new TLegend(.64, .84, .98, .98);
1020     legEff2->SetBorderSize(1);legEff2->SetTextSize(0.03255879);
1021     legEff2->SetFillColor(0);
1022     pad = ((TVirtualPad*)l->At(1));pad->cd();
1023     pad->SetMargin(0.1, 0.01, 0.1, 0.01);
1024     h=(TH1S*)h->Clone("hEff_p");
1025     h->SetYTitle("#epsilon_{p} [%]");
1026     h->GetYaxis()->SetRangeUser(1.e-2, 1.e2);
1027     h->Draw();
1028     content = (TList *)fGraph->FindObject(Form("Eff_%s", AliTRDCalPID::GetPartName(AliPID::kProton)));
1029     if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kLQ))) break;
1030     if(!g->GetN()) break;
1031     legEff2->SetHeader("PID Method [PROTON]");
1032     g->Draw("pc"); legEff2->AddEntry(g, "LQ 2D", "pl");
1033     if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kNN))) break;
1034     g->Draw("pc"); legEff2->AddEntry(g, "NN", "pl");
1035     if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kESD))) break;
1036     g->Draw("p"); legEff2->AddEntry(g, "ESD", "pl");
1037     legEff2->Draw();
1038     gPad->SetLogy();
1039     gPad->SetLogx();
1040     gPad->SetGridy();
1041     gPad->SetGridx();
1042     return kTRUE;
1043   }
1044   case kdEdx:{
1045     // save 2.0 GeV projection as reference
1046     TLegend *legdEdx = new TLegend(.7, .7, .98, .98);
1047     legdEdx->SetBorderSize(1);
1048     kFIRST = kTRUE;
1049     if(!(h2 = (TH2F*)(fContainer->At(kdEdx)))) break;
1050     legdEdx->SetHeader("Particle Species");
1051     gPad->SetMargin(0.1, 0.01, 0.1, 0.01);
1052     for(Int_t is = AliPID::kSPECIES-1; is>=0; is--){
1053       Int_t bin = FindBin(is, 2.);
1054       h1 = h2->ProjectionY(Form("px%d", is), bin, bin);
1055       if(!h1->GetEntries()) continue;
1056       h1->Scale(1./h1->Integral());
1057       h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
1058       if(kFIRST){
1059         h1->GetXaxis()->SetTitle("dE/dx (a.u.)");
1060         h1->GetYaxis()->SetTitle("<Entries>");
1061       }
1062       h = (TH1F*)h1->DrawClone(kFIRST ? "c" : "samec");
1063       legdEdx->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "l");
1064       kFIRST = kFALSE;
1065     }
1066     if(kFIRST) break;
1067     legdEdx->Draw();
1068     gPad->SetLogy();
1069     gPad->SetLogx(0);
1070     gPad->SetGridy();
1071     gPad->SetGridx();
1072     return kTRUE;
1073   }
1074   case kdEdxSlice:
1075     break;
1076   case kPH:{
1077     gPad->Divide(2, 1, 1.e-5, 1.e-5);
1078     TList *l=gPad->GetListOfPrimitives();
1079
1080     // save 2.0 GeV projection as reference
1081     TLegend *legPH = new TLegend(.4, .7, .68, .98);
1082     legPH->SetBorderSize(1);legPH->SetFillColor(0);
1083     legPH->SetHeader("Particle Species");
1084     if(!(arr = (TObjArray*)(fContainer->At(kPH)))) break;
1085     if(!(h2 = (TProfile2D*)(arr->At(0)))) break;
1086
1087     TVirtualPad *pad = ((TVirtualPad*)l->At(0));pad->cd();
1088     pad->SetMargin(0.1, 0.01, 0.1, 0.01);
1089     kFIRST = kTRUE;
1090     for(Int_t is=0; is<AliPID::kSPECIES; is++){
1091       Int_t bin = FindBin(is, 2.);
1092       h1 = h2->ProjectionY(Form("pyt%d", is), bin, bin);
1093       if(!h1->GetEntries()) continue;
1094       h1->SetMarkerStyle(24);
1095       h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
1096       h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
1097       if(kFIRST){
1098         h1->GetXaxis()->SetTitle("t_{drift} [100*ns]");
1099         h1->GetYaxis()->SetTitle("<dQ/dt> [a.u.]");
1100       }
1101       h = (TH1F*)h1->DrawClone(kFIRST ? "c" : "samec");
1102       legPH->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "pl");
1103       kFIRST = kFALSE;
1104     }
1105
1106     pad = ((TVirtualPad*)l->At(1));pad->cd();
1107     pad->SetMargin(0.1, 0.01, 0.1, 0.01);
1108     if(!(h2 = (TProfile2D*)(arr->At(1)))) break;
1109     kFIRST = kTRUE;
1110     for(Int_t is=0; is<AliPID::kSPECIES; is++){
1111       Int_t bin = FindBin(is, 2.);
1112       h1 = h2->ProjectionY(Form("pyx%d", is), bin, bin);
1113       if(!h1->GetEntries()) continue;
1114       h1->SetMarkerStyle(24);
1115       h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
1116       h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
1117       if(kFIRST){
1118         h1->GetXaxis()->SetTitle("x_{drift} [cm]");
1119         h1->GetYaxis()->SetTitle("<dQ/dl> [a.u./cm]");
1120       }
1121       h1->DrawClone(kFIRST ? "c" : "samec");
1122       kFIRST = kFALSE;
1123     }
1124
1125     if(kFIRST) break;
1126     legPH->Draw();
1127     gPad->SetLogy(0);
1128     gPad->SetLogx(0);
1129     gPad->SetGridy();
1130     gPad->SetGridx();
1131     return kTRUE;
1132   }
1133   case kNClus:{
1134     // save 2.0 GeV projection as reference
1135     TLegend *legNClus = new TLegend(.13, .7, .4, .98);
1136     legNClus->SetBorderSize(1);
1137     legNClus->SetFillColor(0);
1138
1139     kFIRST = kTRUE;
1140     if(!(h2 = (TH2F*)(fContainer->At(kNClus)))) break;
1141     legNClus->SetHeader("Particle Species");
1142     for(Int_t is=0; is<AliPID::kSPECIES; is++){
1143       Int_t bin = FindBin(is, 2.);
1144       h1 = h2->ProjectionY(Form("pyNClus%d", is), bin, bin);
1145       if(!h1->GetEntries()) continue;
1146       h1->Scale(100./h1->Integral());
1147       //h1->SetMarkerStyle(24);
1148       //h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
1149       h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
1150       if(kFIRST){
1151         h1->GetXaxis()->SetTitle("N^{cl}/tracklet");
1152         h1->GetYaxis()->SetTitle("Prob. [%]");
1153         h = (TH1F*)h1->DrawClone("c");
1154         h->SetMaximum(20.);
1155         h->GetXaxis()->SetRangeUser(0., 35.);
1156         kFIRST = kFALSE;
1157       } else h = (TH1F*)h1->DrawClone("samec");
1158
1159       legNClus->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "l");
1160     }
1161     if(kFIRST) break;
1162     legNClus->Draw();
1163     gPad->SetLogy(0);
1164     gPad->SetLogx(0);
1165     gPad->SetGridy();
1166     gPad->SetGridx();
1167     return kTRUE;
1168   }
1169   case kMomentum:
1170   case kMomentumBin:
1171     break;
1172   case kNTracklets:{
1173     TLegend *legNClus = new TLegend(.4, .7, .68, .98);
1174     legNClus->SetBorderSize(1);
1175     kFIRST = kTRUE;
1176     if(!(h2 = (TH2F*)(fContainer->At(kNTracklets)))) break;
1177     legNClus->SetHeader("Particle Species");
1178     for(Int_t is=0; is<AliPID::kSPECIES; is++){
1179       Int_t bin = FindBin(is, 2.);
1180       h1 = h2->ProjectionY(Form("pyNTracklets%d", is), bin, bin);
1181       if(!h1->GetEntries()) continue;
1182       h1->Scale(100./h1->Integral());
1183       //h1->SetMarkerStyle(24);
1184       //h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
1185       h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
1186       if(kFIRST){
1187         h1->GetXaxis()->SetTitle("N^{trklt}/track");
1188         h1->GetXaxis()->SetRangeUser(1.,6.);
1189         h1->GetYaxis()->SetTitle("Prob. [%]");
1190       }
1191       h = (TH1F*)h1->DrawClone(kFIRST ? "c" : "samec");
1192       legNClus->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "l");
1193       kFIRST = kFALSE;
1194     }
1195     if(kFIRST) break;
1196     legNClus->Draw();
1197     gPad->SetLogy(0);
1198     gPad->SetLogx(0);
1199     gPad->SetGridy();
1200     gPad->SetGridx();
1201     return kTRUE;
1202   }
1203   }
1204   AliInfo(Form("Reference plot [%d] missing result", ifig));
1205   return kFALSE;
1206 }
1207
1208 //________________________________________________________________________
1209 Bool_t AliTRDcheckPID::PostProcess()
1210 {
1211   // Draw result to the screen
1212   // Called once at the end of the query
1213
1214   if (!fContainer) {
1215     Printf("ERROR: list not available");
1216     return kFALSE;
1217   }
1218
1219   TObjArray *eff(NULL);
1220   if(!fGraph){
1221     fGraph = new TObjArray(2*AliPID::kSPECIES);
1222     fGraph->SetOwner();
1223
1224     if(!(eff = dynamic_cast<TObjArray *>(fContainer->At(kEfficiency)))){
1225       AliError("Efficiency container for Electrons missing.");
1226       return kFALSE;
1227     }
1228     EvaluateEfficiency(eff, fGraph, AliPID::kPion, 0.9);
1229     EvaluateEfficiency(eff, fGraph, AliPID::kKaon, 0.9);
1230     EvaluateEfficiency(eff, fGraph, AliPID::kProton, 0.9);
1231   }
1232   fNRefFigures = 12;
1233   return kTRUE;
1234 }
1235
1236 //________________________________________________________________________
1237 void AliTRDcheckPID::EvaluateEfficiency(const TObjArray * const histoContainer, TObjArray *results, Int_t species, Float_t electronEfficiency){
1238 // Process PID information for pion efficiency
1239
1240   fUtil->SetElectronEfficiency(electronEfficiency);
1241
1242
1243   const Int_t kNmethodsPID=Int_t(sizeof(fgMethod)/sizeof(Char_t*));
1244   Color_t colors[kNmethodsPID] = {kBlue, kGreen+2, kRed};
1245   Int_t markerStyle[kNmethodsPID] = {7, 7, 24};
1246   // efficiency graphs
1247   TGraphErrors *g(NULL);
1248   TObjArray *eff = new TObjArray(kNmethodsPID); eff->SetOwner(); eff->SetName(Form("Eff_%s", AliTRDCalPID::GetPartName(species)));
1249   results->AddAt(eff, species);
1250   for(Int_t iMethod = 0; iMethod < kNmethodsPID; iMethod++){
1251     eff->AddAt(g = new TGraphErrors(), iMethod);
1252     g->SetName(Form("%s", fgMethod[iMethod]));
1253     g->SetLineColor(colors[iMethod]);
1254     g->SetMarkerColor(colors[iMethod]);
1255     g->SetMarkerStyle(markerStyle[iMethod]);
1256   }
1257
1258   // Threshold graphs if not already
1259   TObjArray *thres(NULL);
1260   if(!(results->At(AliPID::kSPECIES))){
1261     thres = new TObjArray(kNmethodsPID); thres->SetOwner();
1262     thres->SetName("Thres");
1263     results->AddAt(thres, AliPID::kSPECIES);
1264     for(Int_t iMethod = 0; iMethod < kNmethodsPID; iMethod++){
1265       thres->AddAt(g = new TGraphErrors(), iMethod);
1266       g->SetName(Form("%s", fgMethod[iMethod]));
1267       g->SetLineColor(colors[iMethod]);
1268       g->SetMarkerColor(colors[iMethod]);
1269       g->SetMarkerStyle(markerStyle[iMethod]);
1270     }
1271   }
1272
1273   TH2F *hPtr[kNmethodsPID]={
1274     (TH2F*)histoContainer->At(AliTRDpidUtil::kLQ),
1275     (TH2F*)histoContainer->At(AliTRDpidUtil::kNN),
1276     (TH2F*)histoContainer->At(AliTRDpidUtil::kESD)
1277   };
1278   for(Int_t iMom = 0; iMom < fMomentumAxis->GetNbins(); iMom++){
1279     Float_t mom(fMomentumAxis->GetBinCenter(iMom+1));
1280
1281     Int_t binEl(fMomentumAxis->GetNbins() * AliPID::kElectron + iMom + 1),
1282                 binXX(fMomentumAxis->GetNbins() * species + iMom + 1);
1283
1284     for(Int_t iMethod = 0; iMethod < kNmethodsPID; iMethod++){
1285       // Calculate the Species Efficiency at electronEfficiency% electron efficiency for each Method
1286
1287       TH1D *histo1 = hPtr[iMethod] -> ProjectionY(Form("%s_el", fgMethod[iMethod]), binEl, binEl);
1288       TH1D *histo2 = hPtr[iMethod] -> ProjectionY(Form("%s_%s", fgMethod[iMethod], AliTRDCalPID::GetPartName(species)), binXX, binXX);
1289
1290       if(!fUtil->CalculatePionEffi(histo1, histo2)) continue;
1291
1292       g=(TGraphErrors*)eff->At(iMethod);
1293       g->SetPoint(iMom, mom, 1.e2*fUtil->GetPionEfficiency());
1294       g->SetPointError(iMom, 0., 1.e2*fUtil->GetError());
1295       AliDebug(2, Form("%s Efficiency for %s is : %f +/- %f", AliTRDCalPID::GetPartName(species), fgMethod[iMethod], fUtil->GetPionEfficiency(), fUtil->GetError()));
1296
1297       if(!thres) continue;
1298       g=(TGraphErrors*)thres->At(iMethod);
1299       g->SetPoint(iMom, mom, fUtil->GetThreshold());
1300       g->SetPointError(iMom, 0., 0.);
1301     }
1302   }
1303 }