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