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