TRD/qaRec/AliTRDpidChecker.cxx

   1 #include "TPDGCode.h"
   2 #include "TCanvas.h"
   3 #include "TF1.h"
   4 #include "TH1F.h"
   5 #include "TH1D.h"
   6 #include "TH2F.h"
   7 #include "TProfile.h"
   8 #include "TProfile2D.h"
   9 #include "TGraph.h"
  10 #include "TGraphErrors.h"
  11
  12 #include <TClonesArray.h>
  13 #include <TObjArray.h>
  14 #include <TList.h>
  15
  16 // #include "AliPID.h"
  17 #include "AliESDEvent.h"
  18 #include "AliESDInputHandler.h"
  19 #include "AliTrackReference.h"
  20
  21 #include "AliAnalysisTask.h"
  22 // #include "AliAnalysisManager.h"
  23
  24 #include "AliTRDtrackerV1.h"
  25 #include "AliTRDtrackV1.h"
  26 #include "AliTRDcluster.h"
  27 #include "AliTRDReconstructor.h"
  28 #include "AliCDBManager.h"
  29 // #include "../Cal/AliTRDCalPID.h"
  30 #include "AliTRDpidUtil.h"
  31
  32
  33 #include "AliTRDpidChecker.h"
  34 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
  35
  36 // calculate pion efficiency at 90% electron efficiency for 11 momentum bins
  37 // this task should be used with simulated data only
  38
  39 ClassImp(AliTRDpidChecker)
  40
  41 //________________________________________________________________________
  42 AliTRDpidChecker::AliTRDpidChecker()
  43   :AliTRDrecoTask("PID", "PID Checker")
  44   ,fReconstructor(0x0)
  45 {
  46   //
  47   // Default constructor
  48   //
  49
  50   fReconstructor = new AliTRDReconstructor();
  51   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
  52 }
  53
  54
  55 //________________________________________________________________________
  56 AliTRDpidChecker::~AliTRDpidChecker()
  57 {
  58   if(fReconstructor) delete fReconstructor;
  59 }
  60
  61
  62 //________________________________________________________________________
  63 void AliTRDpidChecker::CreateOutputObjects()
  64 {
  65   // Create histograms
  66   // Called once
  67
  68   OpenFile(0, "RECREATE");
  69   Int_t xBins = AliPID::kSPECIES*AliTRDCalPID::kNMom;
  70   fContainer = new TObjArray();
  71   fContainer -> Expand(kGraphNN + 1);
  72
  73   const Float_t epsilon = 1./(2*(AliTRDpidUtil::kBins-1));     // get nice histos with bin center at 0 and 1
  74
  75   // histos of the electron probability of all 5 particle species and 11 momenta for the 2-dim LQ method
  76   fContainer->AddAt(
  77     new TH2F("PID_LQ", "",
  78       xBins, -0.5, xBins - 0.5,
  79       AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon)
  80   ,kLQlikelihood);
  81
  82
  83   // histos of the electron probability of all 5 particle species and 11 momenta for the neural network method
  84   fContainer->AddAt(
  85     new TH2F("PID_NN", "",
  86       xBins, -0.5, xBins - 0.5,
  87       AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon)
  88   ,kNNlikelihood);
  89
  90   // histos of the dE/dx distribution for all 5 particle species and 11 momenta
  91   fContainer->AddAt(
  92     new TH2F("dEdx", "",
  93       xBins, -0.5, xBins - 0.5,
  94       200, 0, 10000)
  95     ,kdEdx);
  96
  97   // histos of the pulse height distribution for all 5 particle species and 11 momenta
  98   fContainer->AddAt(
  99     new TProfile2D("PH", "",
 100       xBins, -0.5, xBins - 0.5,
 101       AliTRDtrackerV1::GetNTimeBins(), -0.5, AliTRDtrackerV1::GetNTimeBins() - 0.5)
 102     ,kPH);
 103
 104   // histos of the number of clusters distribution for all 5 particle species and 11 momenta
 105   fContainer->AddAt(
 106     new TH2F("NClus", "",
 107       xBins, -0.5, xBins - 0.5,
 108       AliTRDtrackerV1::GetNTimeBins(), -0.5, AliTRDtrackerV1::GetNTimeBins() - 0.5)
 109     ,kNClus);
 110
 111
 112   // momentum distributions - absolute and in momentum bins
 113   fContainer->AddAt(new TH1F("hMom", "momentum distribution", 100, 0., 12.),kMomentum);
 114   fContainer->AddAt(new TH1F("hMomBin", "momentum distribution in momentum bins", AliTRDCalPID::kNMom, 0.5, 11.5),kMomentumBin);
 115
 116
 117   // TGraph of the pion efficiencies
 118
 119   TGraphErrors *gEffisLQ = 0x0;
 120   TGraphErrors *gEffisNN = 0x0;
 121
 122   fContainer->AddAt(gEffisLQ = new TGraphErrors(), kGraphLQ);
 123   gEffisLQ->SetLineColor(kBlue);
 124   gEffisLQ->SetMarkerColor(kBlue);
 125   gEffisLQ->SetMarkerStyle(29);
 126
 127   fContainer -> AddAt(gEffisNN = new TGraphErrors(),kGraphNN);
 128   gEffisNN->SetLineColor(kRed);
 129   gEffisNN->SetMarkerColor(kRed);
 130   gEffisNN->SetMarkerStyle(29);
 131
 132   gEffisLQ -> SetNameTitle("gEffisLQErr", "Efficiencies and Errors of the 2-dim LQ method");
 133   gEffisNN -> SetNameTitle("gEffisNNErr", "Efficiencies and Errors of the NN method");
 134
 135 }
 136
 137 //________________________________________________________________________
 138 void AliTRDpidChecker::Exec(Option_t *)
 139 {
 140   // Main loop
 141   // Called for each event
 142
 143
 144 //   if(!AliTracker::GetFieldMap()){
 145 //     AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., AliMagFMaps::k5kG);
 146 //     AliTracker::SetFieldMap(field, kTRUE);
 147 //   }
 148
 149   TH2F *hPIDLQ;
 150   TH2F *hPIDNN;
 151   TH2F *hdEdx;
 152   TProfile2D *hPH;
 153   TH2F *hNClus;
 154
 155   hPIDLQ = (TH2F*)fContainer->At(kLQlikelihood);
 156   hPIDNN = (TH2F*)fContainer->At(kNNlikelihood);
 157   hdEdx  = (TH2F*)fContainer->At(kdEdx);
 158   hPH    = (TProfile2D*)fContainer->At(kPH);
 159   hNClus  = (TH2F*)fContainer->At(kNClus);
 160
 161   TH1F *hMom    = (TH1F*)fContainer->At(kMomentum);
 162   TH1F *hMomBin = (TH1F*)fContainer->At(kMomentumBin);
 163
 164   Int_t labelsacc[10000];
 165   memset(labelsacc, 0, sizeof(Int_t) * 10000);
 166
 167   Float_t mom;
 168   ULong_t status;
 169   Int_t nTRD = 0;
 170   Float_t *fdEdx;
 171
 172   AliTRDtrackInfo     *track = 0x0;
 173   AliTRDtrackV1    *TRDtrack = 0x0;
 174   AliTrackReference     *ref = 0x0;
 175   AliExternalTrackParam *esd = 0x0;
 176
 177   AliTRDseedV1 *TRDtracklet[AliTRDgeometry::kNlayer];
 178   for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++)
 179     TRDtracklet[iChamb] = 0x0;
 180
 181   AliTRDcluster *TRDcluster = 0x0;
 182
 183   AliTRDpidUtil *util = new AliTRDpidUtil();
 184   for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
 185     track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
 186     if(!track->HasESDtrack()) continue;
 187     status = track->GetStatus();
 188     if(!(status&AliESDtrack::kTPCout)) continue;
 189
 190     if(!(TRDtrack = track->GetTRDtrack())) continue;
 191     //&&(track->GetNumberOfClustersRefit()
 192
 193     // use only tracks that hit 6 chambers
 194     if(!(TRDtrack->GetNumberOfTracklets() == AliTRDgeometry::kNlayer)) continue;
 195
 196     ref = track->GetTrackRef(0);
 197     esd = track->GetOuterParam();
 198     mom = ref ? ref->P(): esd->P();
 199
 200     labelsacc[nTRD] = track->GetLabel();
 201     nTRD++;
 202
 203
 204     // set the 11 momentum bins
 205     Int_t iMomBin = -1;
 206     iMomBin = util -> GetMomentumBin(mom);
 207     if(fDebugLevel>=4) Printf("MomBin[%d] MomTot[%f]", iMomBin, mom);
 208
 209
 210     // fill momentum histo to have the momentum distribution
 211     hMom -> Fill(mom);
 212     hMomBin -> Fill(iMomBin);
 213
 214
 215     // set the reconstructor
 216     TRDtrack -> SetReconstructor(fReconstructor);
 217
 218
 219     // if no monte carlo data available -> use TRDpid
 220     if(!HasMCdata()){
 221       fReconstructor -> SetOption("nn");
 222       TRDtrack -> CookPID();
 223       if(TRDtrack -> GetPID(0) > TRDtrack -> GetPID(1) + TRDtrack -> GetPID(2)  + TRDtrack -> GetPID(3) + TRDtrack -> GetPID(4)){
 224         track -> SetPDG(kElectron);
 225       }
 226       else if(TRDtrack -> GetPID(4) > TRDtrack -> GetPID(2)  && TRDtrack -> GetPID(4) > TRDtrack -> GetPID(3)  && TRDtrack -> GetPID(4) > TRDtrack -> GetPID(1)){
 227         track -> SetPDG(kProton);
 228       }
 229       else if(TRDtrack -> GetPID(3) > TRDtrack -> GetPID(1)  && TRDtrack -> GetPID(3) > TRDtrack -> GetPID(2)){
 230         track -> SetPDG(kKPlus);
 231       }
 232       else if(TRDtrack -> GetPID(1) > TRDtrack -> GetPID(2)){
 233         track -> SetPDG(kMuonPlus);
 234       }
 235       else{
 236         track -> SetPDG(kPiPlus);
 237       }
 238     }
 239
 240
 241     // calculate the probabilities for electron probability using 2-dim LQ, the deposited charge per chamber and the pulse height spectra and fill histograms
 242     fReconstructor -> SetOption("!nn");
 243     TRDtrack -> CookPID();
 244
 245     if(fDebugLevel>=4) Printf("PIDmethod[%d] Slices[%d] PDG[%d] LQLike[%f]", fReconstructor->GetPIDMethod(), fReconstructor->GetNdEdxSlices(), track->GetPDG(), TRDtrack -> GetPID(0));
 246
 247
 248 //     Bool_t bChange = kFALSE;
 249
 250     Float_t SumdEdx[AliTRDgeometry::kNlayer];
 251     Int_t iNClus[AliTRDgeometry::kNlayer];
 252
 253     for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
 254       TRDtracklet[iChamb] = TRDtrack -> GetTracklet(iChamb);
 255
 256 //       if(!(TRDtracklet[iChamb] -> GetNChange() == 0))
 257 //      bChange = 1;
 258
 259       SumdEdx[iChamb] = 0.;
 260       iNClus[iChamb] = 0;
 261       iNClus[iChamb] = TRDtracklet[iChamb] -> GetN();
 262 //       Printf("NClus[%d]", iNClus[iChamb]);
 263       fdEdx = TRDtracklet[iChamb] -> GetdEdx();
 264       SumdEdx[iChamb] += fdEdx[0] + fdEdx[1] + fdEdx[2];
 265     }
 266
 267 //     if(bChange == kTRUE)
 268 //       continue;
 269
 270     switch(track->GetPDG()){
 271     case kElectron:
 272     case kPositron:
 273       hPIDLQ -> Fill(AliPID::kElectron * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
 274       for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
 275         hdEdx -> Fill(AliPID::kElectron * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
 276         hNClus -> Fill(AliPID::kElectron * AliTRDCalPID::kNMom + iMomBin, iNClus[iChamb]);
 277         for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
 278           if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
 279             continue;
 280           hPH -> Fill(AliPID::kElectron * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
 281         }
 282       }
 283       break;
 284     case kMuonPlus:
 285     case kMuonMinus:
 286       hPIDLQ -> Fill(AliPID::kMuon * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
 287       for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
 288         hdEdx -> Fill(AliPID::kMuon * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
 289         hNClus -> Fill(AliPID::kMuon * AliTRDCalPID::kNMom + iMomBin, iNClus[iChamb]);
 290         for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
 291           if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
 292             continue;
 293           hPH -> Fill(AliPID::kMuon * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
 294         }
 295       }
 296       break;
 297     case kPiPlus:
 298     case kPiMinus:
 299       hPIDLQ -> Fill(AliPID::kPion * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
 300       for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
 301         hdEdx -> Fill(AliPID::kPion * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
 302         hNClus -> Fill(AliPID::kPion * AliTRDCalPID::kNMom + iMomBin, iNClus[iChamb]);
 303         for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
 304           if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
 305             continue;
 306           hPH -> Fill(AliPID::kPion * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
 307         }
 308       }
 309       break;
 310     case kKPlus:
 311     case kKMinus:
 312       hPIDLQ -> Fill(AliPID::kKaon * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
 313       for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
 314         hdEdx -> Fill(AliPID::kKaon * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
 315         hNClus -> Fill(AliPID::kKaon * AliTRDCalPID::kNMom + iMomBin, iNClus[iChamb]);
 316         for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
 317           if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
 318             continue;
 319           hPH -> Fill(AliPID::kKaon * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
 320         }
 321       }
 322       break;
 323     case kProton:
 324     case kProtonBar:
 325       hPIDLQ -> Fill(AliPID::kProton * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
 326       for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
 327         hdEdx -> Fill(AliPID::kProton * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
 328         hNClus -> Fill(AliPID::kProton * AliTRDCalPID::kNMom + iMomBin, iNClus[iChamb]);
 329         for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
 330           if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
 331             continue;
 332           hPH -> Fill(AliPID::kProton * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
 333         }
 334       }
 335       break;
 336     }
 337
 338
 339     // calculate the probabilities and fill histograms for electrons using NN
 340     fReconstructor -> SetOption("nn");
 341     TRDtrack->CookPID();
 342
 343
 344     if(fDebugLevel>=4) Printf("PIDmethod[%d] Slices[%d] PDG[%d] NNLike[%f]", fReconstructor->GetPIDMethod(), fReconstructor->GetNdEdxSlices(), track->GetPDG(), TRDtrack -> GetPID(0));
 345
 346
 347     switch(track->GetPDG()){
 348     case kElectron:
 349     case kPositron:
 350       hPIDNN -> Fill(AliPID::kElectron * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
 351       break;
 352     case kMuonPlus:
 353     case kMuonMinus:
 354       hPIDNN -> Fill(AliPID::kMuon * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
 355       break;
 356     case kPiPlus:
 357     case kPiMinus:
 358       hPIDNN -> Fill(AliPID::kPion * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
 359       break;
 360     case kKPlus:
 361     case kKMinus:
 362       hPIDNN -> Fill(AliPID::kKaon * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
 363       break;
 364     case kProton:
 365     case kProtonBar:
 366       hPIDNN -> Fill(AliPID::kProton * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
 367       break;
 368     }
 369
 370   }
 371
 372   util -> Delete();
 373   PostData(0, fContainer);
 374 }
 375
 376
 377 //________________________________________________________
 378 void AliTRDpidChecker::GetRefFigure(Int_t ifig)
 379 {
 380   Bool_t FIRST = kTRUE;
 381   TGraphErrors *g = 0x0;
 382   TH1 *h1 = 0x0;
 383   TH2 *h2 = 0x0;
 384   switch(ifig){
 385   case 0:
 386     g = (TGraphErrors*)fContainer->At(kGraphStart);
 387     g->Draw("apl");
 388     g->GetHistogram()->GetXaxis()->SetTitle("p [GeV/c]");
 389     g->GetHistogram()->GetXaxis()->SetRangeUser(.6, 10.5);
 390     g->GetHistogram()->GetYaxis()->SetTitle("#epsilon_{#pi} [%]");
 391     ((TGraphErrors*)fContainer->At(kGraphStart+1))->Draw("pl");
 392     gPad->SetLogy();
 393     gPad->SetLogx();
 394     gPad->SetGridy();
 395     gPad->SetGridx();
 396     break;
 397   case 1:
 398     // save 2.0 GeV projection as reference
 399     FIRST = kTRUE;
 400     h2 = (TH2F*)(fContainer->At(kdEdx));
 401     for(Int_t is = AliPID::kSPECIES-1; is>=0; is--){
 402       Int_t bin = is*AliTRDCalPID::kNMom+4;
 403       h1 = h2->ProjectionY("px", bin, bin);
 404       if(!h1->GetEntries()) continue;
 405       h1->Scale(1./h1->Integral());
 406       h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
 407       h1->DrawClone(FIRST ? "c" : "samec");
 408       FIRST = kFALSE;
 409     }
 410     gPad->SetLogy();
 411     gPad->SetLogx(0);
 412     gPad->SetGridy();
 413     gPad->SetGridx();
 414     break;
 415   case 2:
 416     // save 2.0 GeV projection as reference
 417     FIRST = kTRUE;
 418     h2 = (TH2F*)(fContainer->At(kPH));
 419     for(Int_t is=0; is<AliPID::kSPECIES; is++){
 420       Int_t bin = is*AliTRDCalPID::kNMom+4;
 421       h1 = h2->ProjectionY("py", bin, bin);
 422       if(!h1->GetEntries()) continue;
 423       h1->SetMarkerStyle(24);
 424       h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
 425       h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
 426       h1->DrawClone(FIRST ? "e2" : "same e2");
 427       FIRST = kFALSE;
 428     }
 429     gPad->SetLogy(0);
 430     gPad->SetLogx(0);
 431     gPad->SetGridy();
 432     gPad->SetGridx();
 433     break;
 434   case 3:
 435     // save 2.0 GeV projection as reference
 436     FIRST = kTRUE;
 437     h2 = (TH2F*)(fContainer->At(kNClus));
 438     for(Int_t is=0; is<AliPID::kSPECIES; is++){
 439       Int_t bin = is*AliTRDCalPID::kNMom+4;
 440       h1 = h2->ProjectionY("py", bin, bin);
 441       if(!h1->GetEntries()) continue;
 442       h1->SetMarkerStyle(24);
 443       h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
 444       h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
 445       h1->DrawClone(FIRST ? "e2" : "same e2");
 446       FIRST = kFALSE;
 447     }
 448     gPad->SetLogy(0);
 449     gPad->SetLogx(0);
 450     gPad->SetGridy();
 451     gPad->SetGridx();
 452     break;
 453   }
 454 }
 455
 456 //________________________________________________________________________
 457 Bool_t AliTRDpidChecker::PostProcess()
 458 {
 459   // Draw result to the screen
 460   // Called once at the end of the query
 461
 462   if (!fContainer) {
 463     Printf("ERROR: list not available");
 464     return kFALSE;
 465   }
 466 //   return kTRUE; // testing protection
 467
 468
 469   // container for the pion efficiencies and the errors
 470   Double_t  PionEffiLQ[AliTRDCalPID::kNMom],
 471             PionEffiErrorLQ[AliTRDCalPID::kNMom],
 472             EleEffiLQ[AliTRDCalPID::kNMom],
 473             ThresholdLQ[AliTRDCalPID::kNMom];
 474
 475   Double_t  PionEffiNN[AliTRDCalPID::kNMom],
 476             PionEffiErrorNN[AliTRDCalPID::kNMom],
 477             EleEffiNN[AliTRDCalPID::kNMom],
 478             ThresholdNN[AliTRDCalPID::kNMom];
 479
 480   Float_t mom = 0.;
 481
 482   TH1D *Histo1=0x0, *Histo2=0x0;
 483
 484   TH2F *hPIDLQ=0x0, *hPIDNN=0x0;
 485   hPIDLQ = (TH2F*)fContainer->At(kLQlikelihood);
 486   hPIDNN = (TH2F*)fContainer->At(kNNlikelihood);
 487
 488   // calculate the pion efficiencies and the errors for 90% electron efficiency (2-dim LQ)
 489   for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
 490
 491     AliTRDpidUtil *util = new AliTRDpidUtil();
 492     mom = AliTRDCalPID::GetMomentum(iMom);
 493
 494     Histo1 = hPIDLQ -> ProjectionY("LQ_ele",AliTRDCalPID::kNMom*AliPID::kElectron+iMom+1,AliTRDCalPID::kNMom*AliPID::kElectron+iMom+1);
 495     Histo2 = hPIDLQ -> ProjectionY("LQ_pio",AliTRDCalPID::kNMom*AliPID::kPion+iMom+1,AliTRDCalPID::kNMom*AliPID::kPion+iMom+1);
 496
 497     util -> CalculatePionEffi(Histo1, Histo2);
 498
 499     PionEffiLQ[iMom] = util -> GetPionEfficiency();
 500     PionEffiErrorLQ[iMom] = util -> GetError();
 501     EleEffiLQ[iMom] = util -> GetCalcElectronEfficiency();
 502     ThresholdLQ[iMom] = util -> GetThreshold();
 503
 504     if(fDebugLevel>=1) Printf("Pion Efficiency for 2-dim LQ is : %f +/- %f\n\n", PionEffiLQ[iMom], PionEffiErrorLQ[iMom]);
 505     util -> Delete();
 506   }
 507
 508
 509
 510   // calculate the pion efficiencies and the errors for 90% electron efficiency (NN)
 511   for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
 512
 513     AliTRDpidUtil *util = new AliTRDpidUtil();
 514     mom = AliTRDCalPID::GetMomentum(iMom);
 515
 516     Histo1 = hPIDNN -> ProjectionY("NN_ele",AliTRDCalPID::kNMom*AliPID::kElectron+iMom+1,AliTRDCalPID::kNMom*AliPID::kElectron+iMom+1);
 517     Histo2 = hPIDNN -> ProjectionY("NN_pio",AliTRDCalPID::kNMom*AliPID::kPion+iMom+1,AliTRDCalPID::kNMom*AliPID::kPion+iMom+1);
 518
 519     util -> CalculatePionEffi(Histo1, Histo2);
 520
 521     PionEffiNN[iMom] = util -> GetPionEfficiency();
 522     PionEffiErrorNN[iMom] = util -> GetError();
 523     EleEffiNN[iMom] = util -> GetCalcElectronEfficiency();
 524     ThresholdNN[iMom] = util -> GetThreshold();
 525
 526     if(fDebugLevel>=1) Printf("Pion Efficiency for NN is : %f +/- %f\n\n", PionEffiNN[iMom], PionEffiErrorNN[iMom]);
 527
 528     util -> Delete();
 529   }
 530
 531
 532   // create TGraph to plot the pion efficiencies
 533   TGraphErrors *gEffisLQ=0x0, *gEffisNN=0x0;
 534   gEffisLQ = (TGraphErrors*)fContainer->At(kGraphLQ);
 535   gEffisNN = (TGraphErrors*)fContainer->At(kGraphNN);
 536
 537
 538   for(Int_t iBin = 0; iBin < AliTRDCalPID::kNMom; iBin++){
 539
 540     Float_t momentum = AliTRDCalPID::GetMomentum(iBin);
 541     gEffisLQ->SetPoint(iBin, momentum, PionEffiLQ[iBin]);
 542     gEffisLQ->SetPointError(iBin, 0., PionEffiErrorLQ[iBin]);
 543
 544     gEffisNN->SetPoint(iBin, momentum, PionEffiNN[iBin]);
 545     gEffisNN->SetPointError(iBin, 0., PionEffiErrorNN[iBin]);
 546   }
 547
 548
 549   fNRefFigures = 4/*1*/;
 550   return kTRUE; // testing protection
 551 }
 552
 553
 554 //________________________________________________________________________
 555 void AliTRDpidChecker::Terminate(Option_t *)
 556 {
 557   // Draw result to the screen
 558   // Called once at the end of the query
 559
 560   fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
 561   if (!fContainer) {
 562     Printf("ERROR: list not available");
 563     return;
 564   }
 565 }
 566
 567