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