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