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