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