TRD/qaRec/AliTRDcheckPID.cxx

   1 #include "TAxis.h"
   2 #include "TROOT.h"
   3 #include "TPDGCode.h"
   4 #include "TCanvas.h"
   5 #include "TF1.h"
   6 #include "TH1F.h"
   7 #include "TH1D.h"
   8 #include "TH2F.h"
   9 #include "TProfile.h"
  10 #include "TProfile2D.h"
  11 #include "TGraph.h"
  12 #include "TGraphErrors.h"
  13 #include "TLegend.h"
  14
  15 #include <TClonesArray.h>
  16 #include <TObjArray.h>
  17 #include <TList.h>
  18
  19 #include "AliESDEvent.h"
  20 #include "AliESDInputHandler.h"
  21 #include "AliTrackReference.h"
  22
  23 #include "AliAnalysisTask.h"
  24
  25 #include "AliTRDtrackerV1.h"
  26 #include "AliTRDtrackV1.h"
  27 #include "AliTRDcluster.h"
  28 #include "AliTRDReconstructor.h"
  29 #include "AliCDBManager.h"
  30 #include "AliTRDpidUtil.h"
  31
  32 #include "AliTRDcheckPID.h"
  33 #include "info/AliTRDtrackInfo.h"
  34
  35 // calculate pion efficiency at 90% electron efficiency for 11 momentum bins
  36 // this task should be used with simulated data only
  37
  38 ClassImp(AliTRDcheckPID)
  39
  40 //________________________________________________________________________
  41 AliTRDcheckPID::AliTRDcheckPID()
  42   :AliTRDrecoTask("checkPID", "TRD PID checker")
  43   ,fReconstructor(0x0)
  44   ,fUtil(0x0)
  45   ,fGraph(0x0)
  46   ,fEfficiency(0x0)
  47   ,fMomentumAxis(0x0)
  48   ,fMinNTracklets(AliTRDgeometry::kNlayer)
  49   ,fMaxNTracklets(AliTRDgeometry::kNlayer)
  50  {
  51   //
  52   // Default constructor
  53   //
  54
  55   fReconstructor = new AliTRDReconstructor();
  56   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
  57
  58   // Initialize momentum axis with default values
  59   Double_t defaultMomenta[AliTRDCalPID::kNMom+1];
  60   for(Int_t imom = 0; imom < AliTRDCalPID::kNMom+1; imom++)
  61     defaultMomenta[imom] = AliTRDCalPID::GetMomentumBinning(imom);
  62   SetMomentumBinning(AliTRDCalPID::kNMom, defaultMomenta);
  63
  64   fUtil = new AliTRDpidUtil();
  65   InitFunctorList();
  66 }
  67
  68
  69 //________________________________________________________________________
  70 AliTRDcheckPID::~AliTRDcheckPID()
  71 {
  72  if(fGraph){fGraph->Delete(); delete fGraph;}
  73  if(fReconstructor) delete fReconstructor;
  74  if(fUtil) delete fUtil;
  75 }
  76
  77
  78 //________________________________________________________________________
  79 void AliTRDcheckPID::CreateOutputObjects()
  80 {
  81   // Create histograms
  82   // Called once
  83
  84   OpenFile(0, "RECREATE");
  85   fContainer = Histos();
  86 }
  87
  88
  89 //_______________________________________________________
  90 TObjArray * AliTRDcheckPID::Histos(){
  91
  92   //
  93   // Create QA histograms
  94   //
  95   if(fContainer) return fContainer;
  96
  97   Int_t xBins = AliPID::kSPECIES*fMomentumAxis->GetNbins();
  98   fContainer = new TObjArray(); fContainer->Expand(kNPlots);
  99
 100   const Float_t epsilon = 1./(2*(AliTRDpidUtil::kBins-1));     // get nice histos with bin center at 0 and 1
 101
 102   // histos of the electron probability of all 5 particle species and 11 momenta for the 2-dim LQ method
 103   fEfficiency = new TObjArray(3);
 104   fEfficiency->SetOwner(); fEfficiency->SetName("Efficiency");
 105   fContainer->AddAt(fEfficiency, kEfficiency);
 106
 107   TH1 *h = 0x0;
 108   if(!(h = (TH2F*)gROOT->FindObject("PID_LQ"))){
 109     h = new TH2F("PID_LQ", "", xBins, -0.5, xBins - 0.5,
 110       AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon);
 111   } else h->Reset();
 112   fEfficiency->AddAt(h, AliTRDpidUtil::kLQ);
 113
 114   // histos of the electron probability of all 5 particle species and 11 momenta for the neural network method
 115   if(!(h = (TH2F*)gROOT->FindObject("PID_NN"))){
 116     h = new TH2F("PID_NN", "",
 117       xBins, -0.5, xBins - 0.5,
 118       AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon);
 119   } else h->Reset();
 120   fEfficiency->AddAt(h, AliTRDpidUtil::kNN);
 121
 122   // histos of the electron probability of all 5 particle species and 11 momenta for the ESD output
 123   if(!(h = (TH2F*)gROOT->FindObject("PID_ESD"))){
 124     h = new TH2F("PID_ESD", "",
 125       xBins, -0.5, xBins - 0.5,
 126       AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon);
 127   } else h->Reset();
 128   fEfficiency->AddAt(h, AliTRDpidUtil::kESD);
 129
 130   // histos of the dE/dx distribution for all 5 particle species and 11 momenta
 131   if(!(h = (TH2F*)gROOT->FindObject("dEdx"))){
 132     h = new TH2F("dEdx", "",
 133       xBins, -0.5, xBins - 0.5,
 134       200, 0, 10000);
 135   } else h->Reset();
 136   fContainer->AddAt(h, kdEdx);
 137
 138   // histos of the dE/dx slices for all 5 particle species and 11 momenta
 139   if(!(h = (TH2F*)gROOT->FindObject("dEdxSlice"))){
 140     h = new TH2F("dEdxSlice", "",
 141       xBins*AliTRDpidUtil::kLQslices, -0.5, xBins*AliTRDpidUtil::kLQslices - 0.5,
 142       200, 0, 5000);
 143   } else h->Reset();
 144   fContainer->AddAt(h, kdEdxSlice);
 145
 146   // histos of the pulse height distribution for all 5 particle species and 11 momenta
 147   if(!(h = (TH2F*)gROOT->FindObject("PH"))){
 148     h = new TProfile2D("PH", "",
 149       xBins, -0.5, xBins - 0.5,
 150       AliTRDtrackerV1::GetNTimeBins(), -0.5, AliTRDtrackerV1::GetNTimeBins() - 0.5);
 151   } else h->Reset();
 152   fContainer->AddAt(h, kPH);
 153
 154   // histos of the number of clusters distribution for all 5 particle species and 11 momenta
 155   if(!(h = (TH2F*)gROOT->FindObject("NClus"))){
 156     h = new TH2F("NClus", "",
 157       xBins, -0.5, xBins - 0.5,
 158       AliTRDtrackerV1::GetNTimeBins(), -0.5, AliTRDtrackerV1::GetNTimeBins() - 0.5);
 159   } else h->Reset();
 160   fContainer->AddAt(h, kNClus);
 161
 162
 163   // momentum distributions - absolute and in momentum bins
 164   if(!(h = (TH1F*)gROOT->FindObject("hMom"))){
 165     h = new TH1F("hMom", "momentum distribution", fMomentumAxis->GetNbins(), fMomentumAxis->GetXmin(), fMomentumAxis->GetXmax());
 166   } else h->Reset();
 167   fContainer->AddAt(h, kMomentum);
 168
 169   if(!(h = (TH1F*)gROOT->FindObject("hMomBin"))){
 170     h = new TH1F("hMomBin", "momentum distribution in momentum bins", fMomentumAxis->GetNbins(), fMomentumAxis->GetXmin(), fMomentumAxis->GetXmax());
 171   } else h->Reset();
 172   fContainer->AddAt(h, kMomentumBin);
 173
 174   // Number of tracklets per track
 175   if(!(h = (TH2F*)gROOT->FindObject("nTracklets"))){
 176     h = new TH2F("nTracklets", "",
 177       xBins, -0.5, xBins - 0.5,
 178       AliTRDgeometry::kNlayer, 0.5, AliTRDgeometry::kNlayer+.5);
 179   } else h->Reset();
 180   fContainer->AddAt(h, kNTracklets);
 181
 182   return fContainer;
 183 }
 184
 185
 186 //________________________________________________________________________
 187 Bool_t AliTRDcheckPID::CheckTrackQuality(const AliTRDtrackV1* track)
 188 {
 189   //
 190   // Check if the track is ok for PID
 191   //
 192
 193   Int_t ntracklets = track->GetNumberOfTracklets();
 194   if(ntracklets >= fMinNTracklets && ntracklets <= fMaxNTracklets) return 1;
 195 //   if(!fESD)
 196 //     return 0;
 197
 198   return 0;
 199 }
 200
 201 //________________________________________________________________________
 202 Int_t AliTRDcheckPID::CalcPDG(AliTRDtrackV1* track)
 203 {
 204
 205   track -> SetReconstructor(fReconstructor);
 206
 207   fReconstructor -> SetOption("nn");
 208   track -> CookPID();
 209
 210   if(track -> GetPID(AliPID::kElectron) > track -> GetPID(AliPID::kMuon) + track -> GetPID(AliPID::kPion)  + track -> GetPID(AliPID::kKaon) + track -> GetPID(AliPID::kProton)){
 211     return kElectron;
 212   }
 213   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)){
 214     return kProton;
 215   }
 216   else if(track -> GetPID(AliPID::kKaon) > track -> GetPID(AliPID::kMuon)  && track -> GetPID(AliPID::kKaon) > track -> GetPID(AliPID::kPion)){
 217     return kKPlus;
 218   }
 219   else if(track -> GetPID(AliPID::kMuon) > track -> GetPID(AliPID::kPion)){
 220     return kMuonPlus;
 221   }
 222   else{
 223     return kPiPlus;
 224   }
 225 }
 226
 227
 228 //_______________________________________________________
 229 TH1 *AliTRDcheckPID::PlotLQ(const AliTRDtrackV1 *track)
 230 {
 231   //
 232   // Plot the probabilities for electrons using 2-dim LQ
 233   //
 234
 235   if(!fESD){
 236     AliWarning("No ESD info available.");
 237     return 0x0;
 238   }
 239
 240   if(track) fTrack = track;
 241   if(!fTrack){
 242     AliWarning("No Track defined.");
 243     return 0x0;
 244   }
 245
 246   ULong_t status;
 247   status = fESD -> GetStatus();
 248   if(!(status&AliESDtrack::kTRDin)) return 0x0;
 249
 250   if(!CheckTrackQuality(fTrack)) return 0x0;
 251
 252   if(!(fEfficiency = dynamic_cast<TObjArray *>(fContainer->At(kEfficiency)))){
 253     AliWarning("No Histogram defined.");
 254     return 0x0;
 255   }
 256   TH2F *hPIDLQ = 0x0;
 257   if(!(hPIDLQ = dynamic_cast<TH2F *>(fEfficiency->At(AliTRDpidUtil::kLQ)))){
 258     AliWarning("No Histogram defined.");
 259     return 0x0;
 260   }
 261
 262   AliTRDtrackV1 cTrack(*fTrack);
 263   cTrack.SetReconstructor(fReconstructor);
 264
 265   Int_t pdg = 0;
 266   Float_t momentum = 0.;
 267
 268   if(fMC){
 269     if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
 270     pdg = fMC->GetPDG();
 271   } else{
 272     //AliWarning("No MC info available!");
 273     momentum = cTrack.GetMomentum(0);
 274     pdg = CalcPDG(&cTrack);
 275   }
 276   if(!IsInRange(momentum)) return 0x0;
 277
 278   fReconstructor -> SetOption("!nn");
 279   cTrack.CookPID();
 280   if(cTrack.GetNumberOfTrackletsPID() < fMinNTracklets) return 0x0;
 281
 282   Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
 283   hPIDLQ -> Fill(FindBin(species, momentum), cTrack.GetPID(AliPID::kElectron));
 284   return hPIDLQ;
 285 }
 286
 287
 288 //_______________________________________________________
 289 TH1 *AliTRDcheckPID::PlotNN(const AliTRDtrackV1 *track)
 290 {
 291   //
 292   // Plot the probabilities for electrons using 2-dim LQ
 293   //
 294
 295   if(!fESD){
 296     AliWarning("No ESD info available.");
 297     return 0x0;
 298   }
 299
 300   if(track) fTrack = track;
 301   if(!fTrack){
 302     AliWarning("No Track defined.");
 303     return 0x0;
 304   }
 305
 306   ULong_t status;
 307   status = fESD -> GetStatus();
 308   if(!(status&AliESDtrack::kTRDin)) return 0x0;
 309
 310   if(!CheckTrackQuality(fTrack)) return 0x0;
 311
 312   if(!(fEfficiency = dynamic_cast<TObjArray *>(fContainer->At(kEfficiency)))){
 313     AliWarning("No Histogram defined.");
 314     return 0x0;
 315   }
 316   TH2F *hPIDNN;
 317   if(!(hPIDNN = dynamic_cast<TH2F *>(fEfficiency->At(AliTRDpidUtil::kNN)))){
 318     AliWarning("No Histogram defined.");
 319     return 0x0;
 320   }
 321
 322   AliTRDtrackV1 cTrack(*fTrack);
 323   cTrack.SetReconstructor(fReconstructor);
 324
 325   Int_t pdg = 0;
 326   Float_t momentum = 0.;
 327   if(fMC){
 328     if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
 329     pdg = fMC->GetPDG();
 330   } else {
 331     //AliWarning("No MC info available!");
 332     momentum = cTrack.GetMomentum(0);
 333     pdg = CalcPDG(&cTrack);
 334   }
 335   if(!IsInRange(momentum)) return 0x0;
 336
 337   fReconstructor -> SetOption("nn");
 338   cTrack.CookPID();
 339   if(cTrack.GetNumberOfTrackletsPID() < fMinNTracklets) return 0x0;
 340
 341   Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
 342   hPIDNN -> Fill(FindBin(species, momentum), cTrack.GetPID(AliPID::kElectron));
 343   return hPIDNN;
 344 }
 345
 346
 347 //_______________________________________________________
 348 TH1 *AliTRDcheckPID::PlotESD(const AliTRDtrackV1 *track)
 349 {
 350   //
 351   // Plot the probabilities for electrons using 2-dim LQ
 352   //
 353
 354   if(!fESD){
 355     AliWarning("No ESD info available.");
 356     return 0x0;
 357   }
 358
 359   if(track) fTrack = track;
 360   if(!fTrack){
 361     AliWarning("No Track defined.");
 362     return 0x0;
 363   }
 364
 365   ULong_t status;
 366   status = fESD -> GetStatus();
 367   if(!(status&AliESDtrack::kTRDin)) return 0x0;
 368
 369   if(!CheckTrackQuality(fTrack)) return 0x0;
 370   if(fTrack->GetNumberOfTrackletsPID() < fMinNTracklets) return 0x0;
 371
 372   if(!(fEfficiency = dynamic_cast<TObjArray *>(fContainer->At(kEfficiency)))){
 373     AliWarning("No Histogram defined.");
 374     return 0x0;
 375   }
 376   TH2F *hPIDESD = 0x0;
 377   if(!(hPIDESD = dynamic_cast<TH2F *>(fEfficiency->At(AliTRDpidUtil::kESD)))){
 378     AliWarning("No Histogram defined.");
 379     return 0x0;
 380   }
 381
 382
 383   Int_t pdg = 0;
 384   Float_t momentum = 0.;
 385   if(fMC){
 386     if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
 387     pdg = fMC->GetPDG();
 388   } else {
 389     //AliWarning("No MC info available!");
 390     AliTRDtrackV1 cTrack(*fTrack);
 391     cTrack.SetReconstructor(fReconstructor);
 392     momentum = cTrack.GetMomentum(0);
 393     pdg = CalcPDG(&cTrack);
 394   }
 395   if(!IsInRange(momentum)) return 0x0;
 396
 397 //   Double32_t pidESD[AliPID::kSPECIES];
 398   const Double32_t *pidESD = fESD->GetResponseIter();
 399   Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
 400   hPIDESD -> Fill(FindBin(species, momentum), pidESD[0]);
 401   return hPIDESD;
 402 }
 403
 404
 405 //_______________________________________________________
 406 TH1 *AliTRDcheckPID::PlotdEdx(const AliTRDtrackV1 *track)
 407 {
 408   //
 409   // Plot the probabilities for electrons using 2-dim LQ
 410   //
 411
 412   if(track) fTrack = track;
 413   if(!fTrack){
 414     AliWarning("No Track defined.");
 415     return 0x0;
 416   }
 417
 418   if(!CheckTrackQuality(fTrack)) return 0x0;
 419
 420   TH2F *hdEdx;
 421   if(!(hdEdx = dynamic_cast<TH2F *>(fContainer->At(kdEdx)))){
 422     AliWarning("No Histogram defined.");
 423     return 0x0;
 424   }
 425
 426   AliTRDtrackV1 cTrack(*fTrack);
 427   cTrack.SetReconstructor(fReconstructor);
 428   Int_t pdg = 0;
 429   Float_t momentum = 0.;
 430   if(fMC){
 431     if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
 432     pdg = fMC->GetPDG();
 433   } else {
 434     //AliWarning("No MC info available!");
 435     momentum = cTrack.GetMomentum(0);
 436     pdg = CalcPDG(&cTrack);
 437   }
 438   Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
 439   if(!IsInRange(momentum)) return 0x0;
 440
 441   fReconstructor -> SetOption("!nn");
 442   Float_t SumdEdx = 0;
 443   Int_t iBin = FindBin(species, momentum);
 444   AliTRDseedV1 *tracklet = 0x0;
 445   for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
 446     SumdEdx = 0;
 447     tracklet = cTrack.GetTracklet(iChamb);
 448     if(!tracklet) continue;
 449     tracklet -> CookdEdx(AliTRDpidUtil::kLQslices);
 450     for(Int_t i = 3; i--;) SumdEdx += tracklet->GetdEdx()[i];
 451     hdEdx -> Fill(iBin, SumdEdx);
 452   }
 453   return hdEdx;
 454 }
 455
 456
 457 //_______________________________________________________
 458 TH1 *AliTRDcheckPID::PlotdEdxSlice(const AliTRDtrackV1 *track)
 459 {
 460   //
 461   // Plot the probabilities for electrons using 2-dim LQ
 462   //
 463
 464   if(track) fTrack = track;
 465   if(!fTrack){
 466     AliWarning("No Track defined.");
 467     return 0x0;
 468   }
 469
 470   if(!CheckTrackQuality(fTrack)) return 0x0;
 471
 472   TH2F *hdEdxSlice;
 473   if(!(hdEdxSlice = dynamic_cast<TH2F *>(fContainer->At(kdEdxSlice)))){
 474     AliWarning("No Histogram defined.");
 475     return 0x0;
 476   }
 477
 478   AliTRDtrackV1 cTrack(*fTrack);
 479   cTrack.SetReconstructor(fReconstructor);
 480   Int_t pdg = 0;
 481   Float_t momentum = 0.;
 482   if(fMC){
 483     if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
 484     pdg = fMC->GetPDG();
 485   } else {
 486     //AliWarning("No MC info available!");
 487     momentum = cTrack.GetMomentum(0);
 488     pdg = CalcPDG(&cTrack);
 489   }
 490   if(!IsInRange(momentum)) return 0x0;
 491
 492   fReconstructor -> SetOption("!nn");
 493   Int_t iMomBin = fMomentumAxis->FindBin(momentum);
 494   Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
 495   Float_t *fdEdx;
 496   AliTRDseedV1 *tracklet = 0x0;
 497   for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
 498     tracklet = cTrack.GetTracklet(iChamb);
 499     if(!tracklet) continue;
 500     tracklet -> CookdEdx(AliTRDpidUtil::kLQslices);
 501     fdEdx = tracklet->GetdEdx();
 502     for(Int_t iSlice = 0; iSlice < AliTRDpidUtil::kLQslices; iSlice++){
 503       hdEdxSlice -> Fill(species * fMomentumAxis->GetNbins() * AliTRDpidUtil::kLQslices + (iMomBin-1) * AliTRDpidUtil::kLQslices + iSlice, fdEdx[iSlice]);
 504     }
 505   }
 506
 507   return hdEdxSlice;
 508 }
 509
 510
 511 //_______________________________________________________
 512 TH1 *AliTRDcheckPID::PlotPH(const AliTRDtrackV1 *track)
 513 {
 514   //
 515   // Plot the probabilities for electrons using 2-dim LQ
 516   //
 517
 518   if(track) fTrack = track;
 519   if(!fTrack){
 520     AliWarning("No Track defined.");
 521     return 0x0;
 522   }
 523
 524   if(!CheckTrackQuality(fTrack)) return 0x0;
 525
 526   TProfile2D *hPH;
 527   if(!(hPH = dynamic_cast<TProfile2D *>(fContainer->At(kPH)))){
 528     AliWarning("No Histogram defined.");
 529     return 0x0;
 530   }
 531
 532   Int_t pdg = 0;
 533   Float_t momentum = 0.;
 534   if(fMC){
 535     if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
 536     pdg = fMC->GetPDG();
 537   } else {
 538     //AliWarning("No MC info available!");
 539     AliTRDtrackV1 cTrack(*fTrack);
 540     cTrack.SetReconstructor(fReconstructor);
 541     momentum = cTrack.GetMomentum(0);
 542     pdg = CalcPDG(&cTrack);
 543   }
 544   if(!IsInRange(momentum)) return 0x0;;
 545
 546   AliTRDseedV1 *tracklet = 0x0;
 547   AliTRDcluster *TRDcluster = 0x0;
 548   Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
 549   Int_t iBin = FindBin(species, momentum);
 550   for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
 551     tracklet = fTrack->GetTracklet(iChamb);
 552     if(!tracklet) continue;
 553     for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
 554       if(!(TRDcluster = tracklet->GetClusters(iClus))) continue;
 555       hPH -> Fill(iBin, TRDcluster->GetLocalTimeBin(), tracklet->GetdQdl(iClus));
 556     }
 557   }
 558   return hPH;
 559 }
 560
 561
 562 //_______________________________________________________
 563 TH1 *AliTRDcheckPID::PlotNClus(const AliTRDtrackV1 *track)
 564 {
 565   //
 566   // Plot the probabilities for electrons using 2-dim LQ
 567   //
 568
 569   if(track) fTrack = track;
 570   if(!fTrack){
 571     AliWarning("No Track defined.");
 572     return 0x0;
 573   }
 574
 575   if(!CheckTrackQuality(fTrack)) return 0x0;
 576
 577   TH2F *hNClus;
 578   if(!(hNClus = dynamic_cast<TH2F *>(fContainer->At(kNClus)))){
 579     AliWarning("No Histogram defined.");
 580     return 0x0;
 581   }
 582
 583
 584   Int_t pdg = 0;
 585   Float_t momentum = 0.;
 586   if(fMC){
 587     if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
 588     pdg = fMC->GetPDG();
 589   } else {
 590     //AliWarning("No MC info available!");
 591     AliTRDtrackV1 cTrack(*fTrack);
 592     cTrack.SetReconstructor(fReconstructor);
 593     momentum = cTrack.GetMomentum(0);
 594     pdg = CalcPDG(&cTrack);
 595   }
 596   if(!IsInRange(momentum)) return 0x0;
 597
 598   Int_t species = AliTRDpidUtil::AliTRDpidUtil::Pdg2Pid(pdg);
 599   Int_t iBin = FindBin(species, momentum);
 600   AliTRDseedV1 *tracklet = 0x0;
 601   for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
 602     tracklet = fTrack->GetTracklet(iChamb);
 603     if(!tracklet) continue;
 604     hNClus -> Fill(iBin, tracklet->GetN());
 605   }
 606
 607   return hNClus;
 608 }
 609
 610 //_______________________________________________________
 611 TH1 *AliTRDcheckPID::PlotNTracklets(const AliTRDtrackV1 *track)
 612 {
 613   //
 614   // Plot the probabilities for electrons using 2-dim LQ
 615   //
 616
 617   if(track) fTrack = track;
 618   if(!fTrack){
 619     AliWarning("No Track defined.");
 620     return 0x0;
 621   }
 622
 623   TH2F *hTracklets;
 624   if(!(hTracklets = dynamic_cast<TH2F *>(fContainer->At(kNTracklets)))){
 625     AliWarning("No Histogram defined.");
 626     return 0x0;
 627   }
 628
 629   AliTRDtrackV1 cTrack(*fTrack);
 630   cTrack.SetReconstructor(fReconstructor);
 631   Int_t pdg = 0;
 632   Float_t momentum = 0.;
 633   if(fMC){
 634     if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
 635     pdg = fMC->GetPDG();
 636   } else {
 637     //AliWarning("No MC info available!");
 638     momentum = cTrack.GetMomentum(0);
 639     pdg = CalcPDG(&cTrack);
 640   }
 641   Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
 642   if(!IsInRange(momentum)) return 0x0;
 643
 644   Int_t iBin = FindBin(species, momentum);
 645   hTracklets -> Fill(iBin, cTrack.GetNumberOfTracklets());
 646   return hTracklets;
 647 }
 648
 649 //_______________________________________________________
 650 TH1 *AliTRDcheckPID::PlotMom(const AliTRDtrackV1 *track)
 651 {
 652   //
 653   // Plot the probabilities for electrons using 2-dim LQ
 654   //
 655
 656   if(track) fTrack = track;
 657   if(!fTrack){
 658     AliWarning("No Track defined.");
 659     return 0x0;
 660   }
 661
 662   if(!CheckTrackQuality(fTrack)) return 0x0;
 663
 664   TH1F *hMom;
 665   if(!(hMom = dynamic_cast<TH1F *>(fContainer->At(kMomentum)))){
 666     AliWarning("No Histogram defined.");
 667     return 0x0;
 668   }
 669
 670
 671   Int_t pdg = 0;
 672   Float_t momentum = 0.;
 673   if(fMC){
 674     if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
 675     pdg = fMC->GetPDG();
 676   } else {
 677     //AliWarning("No MC info available!");
 678     AliTRDtrackV1 cTrack(*fTrack);
 679     cTrack.SetReconstructor(fReconstructor);
 680     momentum = cTrack.GetMomentum(0);
 681     pdg = CalcPDG(&cTrack);
 682   }
 683   if(IsInRange(momentum)) hMom -> Fill(momentum);
 684   return hMom;
 685 }
 686
 687
 688 //_______________________________________________________
 689 TH1 *AliTRDcheckPID::PlotMomBin(const AliTRDtrackV1 *track)
 690 {
 691   //
 692   // Plot the probabilities for electrons using 2-dim LQ
 693   //
 694
 695   if(track) fTrack = track;
 696   if(!fTrack){
 697     AliWarning("No Track defined.");
 698     return 0x0;
 699   }
 700
 701   if(!CheckTrackQuality(fTrack)) return 0x0;
 702
 703   TH1F *hMomBin;
 704   if(!(hMomBin = dynamic_cast<TH1F *>(fContainer->At(kMomentumBin)))){
 705     AliWarning("No Histogram defined.");
 706     return 0x0;
 707   }
 708
 709
 710   Int_t pdg = 0;
 711   Float_t momentum = 0.;
 712
 713   if(fMC){
 714     if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
 715     pdg = fMC->GetPDG();
 716   } else {
 717     //AliWarning("No MC info available!");
 718     AliTRDtrackV1 cTrack(*fTrack);
 719     cTrack.SetReconstructor(fReconstructor);
 720     momentum = cTrack.GetMomentum(0);
 721   }
 722   if(IsInRange(momentum)) hMomBin -> Fill(fMomentumAxis->FindBin(momentum));
 723   return hMomBin;
 724 }
 725
 726
 727 //________________________________________________________
 728 Bool_t AliTRDcheckPID::GetRefFigure(Int_t ifig)
 729 {
 730   Bool_t FIRST = kTRUE;
 731   TGraphErrors *g = 0x0;
 732   TAxis *ax = 0x0;
 733   TH1 *h1 = 0x0, *h=0x0;
 734   TH2 *h2 = 0x0;
 735   TList *content = 0x0;
 736   switch(ifig){
 737   case kEfficiency:{
 738     gPad->Divide(2, 1, 0.,0.);
 739     TList *l=gPad->GetListOfPrimitives();
 740     TVirtualPad *pad = ((TVirtualPad*)l->At(0));pad->cd();
 741     pad->SetMargin(0.07, 0.07, 0.1, 0.001);
 742
 743     TLegend *legEff = new TLegend(.7, .7, .98, .98);
 744     legEff->SetBorderSize(1);
 745     content = (TList *)fGraph->FindObject("Efficiencies");
 746     if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kLQ))) break;
 747     if(!g->GetN()) break;
 748     legEff->SetHeader("PID Method");
 749     g->Draw("apl");
 750     ax = g->GetHistogram()->GetXaxis();
 751     ax->SetTitle("p (GeV/c)");
 752     ax->SetRangeUser(.5, 11.);
 753     ax->SetMoreLogLabels();
 754     ax = g->GetHistogram()->GetYaxis();
 755     ax->SetTitle("#epsilon_{#pi}");
 756     ax->SetRangeUser(1.e-3, 1.e-1);
 757     legEff->AddEntry(g, "2D LQ", "pl");
 758     if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kNN))) break;
 759     g->Draw("pl");
 760     legEff->AddEntry(g, "NN", "pl");
 761     if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kESD))) break;
 762     g->Draw("p");
 763     legEff->AddEntry(g, "ESD", "pl");
 764     legEff->Draw();
 765     gPad->SetLogy();
 766     gPad->SetLogx();
 767     gPad->SetGridy();
 768     gPad->SetGridx();
 769
 770     pad = ((TVirtualPad*)l->At(1));pad->cd();
 771     pad->SetMargin(0.07, 0.07, 0.1, 0.001);
 772     content = (TList *)fGraph->FindObject("Thresholds");
 773     if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kLQ))) break;
 774     if(!g->GetN()) break;
 775     g->Draw("apl");
 776     ax = g->GetHistogram()->GetXaxis();
 777     ax->SetTitle("p (GeV/c)");
 778     ax->SetRangeUser(.5, 11.5);
 779     ax->SetMoreLogLabels();
 780     ax = g->GetHistogram()->GetYaxis();
 781     ax->SetTitle("threshold (%)");
 782     ax->SetRangeUser(5.e-2, 1.);
 783     if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kNN))) break;
 784     g->Draw("pl");
 785     if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kESD))) break;
 786     g->Draw("p");
 787     gPad->SetLogx();
 788     gPad->SetGridy();
 789     gPad->SetGridx();
 790     return kTRUE;
 791   }
 792   case kdEdx:{
 793     // save 2.0 GeV projection as reference
 794     TLegend *legdEdx = new TLegend(.7, .7, .98, .98);
 795     legdEdx->SetBorderSize(1);
 796     FIRST = kTRUE;
 797     if(!(h2 = (TH2F*)(fContainer->At(kdEdx)))) break;
 798     legdEdx->SetHeader("Particle Species");
 799     for(Int_t is = AliPID::kSPECIES-1; is>=0; is--){
 800       Int_t bin = FindBin(is, 2.);
 801       h1 = h2->ProjectionY("px", bin, bin);
 802       if(!h1->GetEntries()) continue;
 803       h1->Scale(1./h1->Integral());
 804       h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
 805       if(FIRST){
 806         h1->GetXaxis()->SetTitle("dE/dx (a.u.)");
 807         h1->GetYaxis()->SetTitle("<Entries>");
 808       }
 809       h = (TH1F*)h1->DrawClone(FIRST ? "c" : "samec");
 810       legdEdx->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "l");
 811       FIRST = kFALSE;
 812     }
 813     if(FIRST) break;
 814     legdEdx->Draw();
 815     gPad->SetLogy();
 816     gPad->SetLogx(0);
 817     gPad->SetGridy();
 818     gPad->SetGridx();
 819     return kTRUE;
 820   }
 821   case kdEdxSlice:
 822     break;
 823   case kPH:{
 824     // save 2.0 GeV projection as reference
 825     TLegend *legPH = new TLegend(.4, .7, .68, .98);
 826     legPH->SetBorderSize(1);
 827     FIRST = kTRUE;
 828     if(!(h2 = (TH2F*)(fContainer->At(kPH)))) break;;
 829     legPH->SetHeader("Particle Species");
 830     for(Int_t is=0; is<AliPID::kSPECIES; is++){
 831       Int_t bin = FindBin(is, 2.);
 832       h1 = h2->ProjectionY("py", bin, bin);
 833       if(!h1->GetEntries()) continue;
 834       h1->SetMarkerStyle(24);
 835       h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
 836       h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
 837       if(FIRST){
 838         h1->GetXaxis()->SetTitle("tb(1/100 ns^{-1})");
 839         h1->GetYaxis()->SetTitle("<PH> (a.u.)");
 840       }
 841       h = (TH1F*)h1->DrawClone(FIRST ? "c" : "samec");
 842       legPH->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "pl");
 843       FIRST = kFALSE;
 844     }
 845     if(FIRST) break;
 846     legPH->Draw();
 847     gPad->SetLogy(0);
 848     gPad->SetLogx(0);
 849     gPad->SetGridy();
 850     gPad->SetGridx();
 851     return kTRUE;
 852   }
 853   case kNClus:{
 854     // save 2.0 GeV projection as reference
 855     TLegend *legNClus = new TLegend(.4, .7, .68, .98);
 856     legNClus->SetBorderSize(1);
 857     FIRST = kTRUE;
 858     if(!(h2 = (TH2F*)(fContainer->At(kNClus)))) break;
 859     legNClus->SetHeader("Particle Species");
 860     for(Int_t is=0; is<AliPID::kSPECIES; is++){
 861       Int_t bin = FindBin(is, 2.);
 862       h1 = h2->ProjectionY("py", bin, bin);
 863       if(!h1->GetEntries()) continue;
 864       h1->Scale(1./h1->Integral());
 865       //h1->SetMarkerStyle(24);
 866       //h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
 867       h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
 868       if(FIRST) h1->GetXaxis()->SetTitle("N^{cl}/tracklet");
 869       if(FIRST) h1->GetYaxis()->SetTitle("<Entries>");
 870       h = (TH1F*)h1->DrawClone(FIRST ? "c" : "samec");
 871       legNClus->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "l");
 872       FIRST = kFALSE;
 873     }
 874     if(FIRST) break;
 875     legNClus->Draw();
 876     gPad->SetLogy();
 877     gPad->SetLogx(0);
 878     gPad->SetGridy();
 879     gPad->SetGridx();
 880     return kTRUE;
 881   }
 882   case kMomentum:
 883   case kMomentumBin:
 884     break;
 885   case kNTracklets:{
 886     TLegend *legNClus = new TLegend(.4, .7, .68, .98);
 887     legNClus->SetBorderSize(1);
 888     FIRST = kTRUE;
 889     if(!(h2 = (TH2F*)(fContainer->At(kNTracklets)))) break;
 890     legNClus->SetHeader("Particle Species");
 891     for(Int_t is=0; is<AliPID::kSPECIES; is++){
 892       Int_t bin = FindBin(is, 2.);
 893       h1 = h2->ProjectionY("py", bin, bin);
 894       if(!h1->GetEntries()) continue;
 895       h1->Scale(1./h1->Integral());
 896       //h1->SetMarkerStyle(24);
 897       //h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
 898       h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
 899       if(FIRST) h1->GetXaxis()->SetTitle("N^{tl}/track");
 900       if(FIRST) h1->GetYaxis()->SetTitle("<Entries>");
 901       h = (TH1F*)h1->DrawClone(FIRST ? "c" : "samec");
 902       legNClus->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "l");
 903       FIRST = kFALSE;
 904     }
 905     if(FIRST) break;
 906     legNClus->Draw();
 907     gPad->SetLogy();
 908     gPad->SetLogx(0);
 909     gPad->SetGridy();
 910     gPad->SetGridx();
 911     return kTRUE;
 912   }
 913   }
 914   AliInfo(Form("Reference plot [%d] missing result", ifig));
 915   return kFALSE;
 916 }
 917
 918 //________________________________________________________________________
 919 Bool_t AliTRDcheckPID::PostProcess()
 920 {
 921   // Draw result to the screen
 922   // Called once at the end of the query
 923
 924   if (!fContainer) {
 925     Printf("ERROR: list not available");
 926     return kFALSE;
 927   }
 928   if(!(fEfficiency = dynamic_cast<TObjArray *>(fContainer->At(kEfficiency)))){
 929     AliError("Efficiency container missing.");
 930     return 0x0;
 931   }
 932   if(!fGraph){
 933     fGraph = new TObjArray(6);
 934     fGraph->SetOwner();
 935     EvaluatePionEfficiency(fEfficiency, fGraph, 0.9);
 936   }
 937   fNRefFigures = 9;
 938   return kTRUE;
 939 }
 940
 941 //________________________________________________________________________
 942 void AliTRDcheckPID::EvaluatePionEfficiency(TObjArray *histoContainer, TObjArray *results, Float_t electron_efficiency){
 943   fUtil->SetElectronEfficiency(electron_efficiency);
 944
 945   Color_t colors[3] = {kBlue, kGreen+2, kRed};
 946   Int_t markerStyle[3] = {7, 7, 24};
 947   TString MethodName[3] = {"LQ", "NN", "ESD"};
 948   // efficiency graphs
 949   TGraphErrors *g, *gPtrEff[3], *gPtrThres[3];
 950   TObjArray *eff = new TObjArray(3); eff->SetOwner(); eff->SetName("Efficiencies");
 951   results->AddAt(eff, 0);
 952   for(Int_t iMethod = 0; iMethod < 3; iMethod++){
 953     eff->AddAt(g = gPtrEff[iMethod] = new TGraphErrors(), iMethod);
 954     g->SetName(Form("efficiency_%s", MethodName[iMethod].Data()));
 955     g->SetLineColor(colors[iMethod]);
 956     g->SetMarkerColor(colors[iMethod]);
 957     g->SetMarkerStyle(markerStyle[iMethod]);
 958   }
 959
 960   // Threshold graphs
 961   TObjArray *thres = new TObjArray(3); thres->SetOwner(); thres->SetName("Thresholds");
 962   results->AddAt(thres, 1);
 963   for(Int_t iMethod = 0; iMethod < 3; iMethod++){
 964     thres->AddAt(g = gPtrThres[iMethod] = new TGraphErrors(), iMethod);
 965     g->SetName(Form("threshold_%s", MethodName[iMethod].Data()));
 966     g->SetLineColor(colors[iMethod]);
 967     g->SetMarkerColor(colors[iMethod]);
 968     g->SetMarkerStyle(markerStyle[iMethod]);
 969   }
 970
 971   Float_t mom = 0.;
 972   TH1D *Histo1=0x0, *Histo2=0x0;
 973
 974   TH2F *hPtr[3];
 975   hPtr[0] = (TH2F*)histoContainer->At(AliTRDpidUtil::kLQ);
 976   hPtr[1] = (TH2F*)histoContainer->At(AliTRDpidUtil::kNN);
 977   hPtr[2] = (TH2F*)histoContainer->At(AliTRDpidUtil::kESD);
 978
 979   for(Int_t iMom = 0; iMom < fMomentumAxis->GetNbins(); iMom++){
 980     mom = fMomentumAxis->GetBinCenter(iMom+1);
 981
 982     Int_t binEl = fMomentumAxis->GetNbins() * AliPID::kElectron + iMom + 1,
 983           binPi = fMomentumAxis->GetNbins() * AliPID::kPion + iMom + 1;
 984     for(Int_t iMethod = 0; iMethod < 3; iMethod++){
 985       // Calculate the Pion Efficiency at 90% electron efficiency for each Method
 986       Histo1 = hPtr[iMethod] -> ProjectionY(Form("%s_ele", MethodName[iMethod].Data()), binEl, binEl);
 987       Histo2 = hPtr[iMethod] -> ProjectionY(Form("%s_pio", MethodName[iMethod].Data()), binPi, binPi);
 988
 989       if(!fUtil->CalculatePionEffi(Histo1, Histo2)) continue;
 990
 991       gPtrEff[iMethod]->SetPoint(iMom, mom, fUtil->GetPionEfficiency());
 992       gPtrEff[iMethod]->SetPointError(iMom, 0., fUtil->GetError());
 993       gPtrThres[iMethod]->SetPoint(iMom, mom, fUtil->GetThreshold());
 994       gPtrThres[iMethod]->SetPointError(iMom, 0., 0.);
 995
 996       if(fDebugLevel>=2) Printf(Form("Pion Efficiency for 2-dim LQ is : %f +/- %f\n\n", MethodName[iMethod].Data()), fUtil->GetPionEfficiency(), fUtil->GetError());
 997     }
 998   }
 999 }
1000
1001 //________________________________________________________________________
1002 void AliTRDcheckPID::Terminate(Option_t *)
1003 {
1004   // Draw result to the screen
1005   // Called once at the end of the query
1006
1007   fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
1008   if (!fContainer) {
1009     Printf("ERROR: list not available");
1010     return;
1011   }
1012 }
1013
1014