plot fixes
[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDcheckDET.cxx
1 ////////////////////////////////////////////////////////////////////////////
2 //                                                                        //
3 //                                                                        //
4 //  Basic checks for tracking and detector performance                    //
5 //  
6 //     For the moment (15.10.2009) the following checks are implemented    //
7 //       - Number of clusters/track
8 //       - Number of clusters/tracklet
9 //       - Number of tracklets/track from different sources (Barrel, StandAlone)
10 //       - Number of findable tracklets
11 //       - Number of tracks per event and TRD sector
12 //       - <PH>
13 //       - Chi2 distribution for tracks
14 //       - Charge distribution per cluster and tracklet
15 //       - Number of events and tracks per trigger source 
16 //       - Trigger purity
17 //       - Track and Tracklet propagation status
18 //
19 //  Authors:                                                              //
20 //    Anton Andronic <A.Andronic@gsi.de>                                  //
21 //    Alexandru Bercuci <A.Bercuci@gsi.de>                                //
22 //    Markus Fasel <M.Fasel@gsi.de>                                       //
23 //                                                                        //
24 ////////////////////////////////////////////////////////////////////////////
25
26 #include <TAxis.h>
27 #include <TCanvas.h>
28 #include <TFile.h>
29 #include <TH1F.h>
30 #include <TH1I.h>
31 #include <TH2F.h>
32 #include <TF1.h>
33 #include <TGaxis.h>
34 #include <TGraph.h>
35 #include <TLegend.h>
36 #include <TLinearFitter.h>
37 #include <TMath.h>
38 #include <TMap.h>
39 #include <TObjArray.h>
40 #include <TObject.h>
41 #include <TObjString.h>
42
43 #include <TPad.h>
44 #include <TProfile.h>
45 #include <TProfile2D.h>
46 #include <TROOT.h>
47 #include <TChain.h>
48
49 #include "AliLog.h"
50 #include "AliTRDcluster.h"
51 #include "AliESDHeader.h"
52 #include "AliESDRun.h"
53 #include "AliESDtrack.h"
54 #include "AliExternalTrackParam.h"
55 #include "AliTRDgeometry.h"
56 #include "AliTRDpadPlane.h"
57 #include "AliTRDSimParam.h"
58 #include "AliTRDseedV1.h"
59 #include "AliTRDtrackV1.h"
60 #include "AliTRDtrackerV1.h"
61 #include "AliTRDReconstructor.h"
62 #include "AliTrackReference.h"
63 #include "AliTrackPointArray.h"
64 #include "AliTracker.h"
65 #include "TTreeStream.h"
66
67 #include "info/AliTRDtrackInfo.h"
68 #include "info/AliTRDeventInfo.h"
69 #include "AliTRDcheckDET.h"
70
71 #include <cstdio>
72 #include <iostream>
73
74 ClassImp(AliTRDcheckDET)
75
76 //_______________________________________________________
77 AliTRDcheckDET::AliTRDcheckDET():
78   AliTRDrecoTask()
79   ,fEventInfo(NULL)
80   ,fTriggerNames(NULL)
81   ,fReconstructor(NULL)
82   ,fGeo(NULL)
83   ,fFlags(0)
84 {
85   //
86   // Default constructor
87   //
88   SetNameTitle("checkDET", "Basic TRD data checker");
89 }
90
91 //_______________________________________________________
92 AliTRDcheckDET::AliTRDcheckDET(char* name):
93   AliTRDrecoTask(name, "Basic TRD data checker")
94   ,fEventInfo(NULL)
95   ,fTriggerNames(NULL)
96   ,fReconstructor(NULL)
97   ,fGeo(NULL)
98   ,fFlags(0)
99 {
100   //
101   // Default constructor
102   //
103   DefineInput(2, AliTRDeventInfo::Class());
104
105   fReconstructor = new AliTRDReconstructor;
106   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
107   fGeo = new AliTRDgeometry;
108   InitFunctorList();
109 }
110
111
112 //_______________________________________________________
113 AliTRDcheckDET::~AliTRDcheckDET(){
114   //
115   // Destructor
116   // 
117   if(fTriggerNames) delete fTriggerNames;
118   delete fReconstructor;
119   delete fGeo;
120 }
121
122 //_______________________________________________________
123 void AliTRDcheckDET::ConnectInputData(Option_t *opt){
124   //
125   // Connect the Input data with the task
126   //
127   AliTRDrecoTask::ConnectInputData(opt);
128   fEventInfo = dynamic_cast<AliTRDeventInfo *>(GetInputData(2));
129 }
130
131 //_______________________________________________________
132 void AliTRDcheckDET::UserCreateOutputObjects(){
133   //
134   // Create Output Objects
135   //
136   OpenFile(1,"RECREATE");
137   fContainer = Histos();
138   if(!fTriggerNames) fTriggerNames = new TMap();
139 }
140
141 //_______________________________________________________
142 void AliTRDcheckDET::UserExec(Option_t *opt){
143   //
144   // Execution function
145   // Filling TRD quality histos
146   //
147   if(!HasMCdata() && fEventInfo->GetEventHeader()->GetEventType() != 7) return; // For real data we select only physical events
148   AliTRDrecoTask::UserExec(opt);  
149   Int_t nTracks = 0;            // Count the number of tracks per event
150   Int_t triggermask = fEventInfo->GetEventHeader()->GetTriggerMask();
151   TString triggername =  fEventInfo->GetRunInfo()->GetFiredTriggerClasses(triggermask);
152   AliDebug(6, Form("Trigger cluster: %d, Trigger class: %s\n", triggermask, triggername.Data()));
153   dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger))->Fill(triggermask);
154   for(Int_t iti = 0; iti < fTracks->GetEntriesFast(); iti++){
155     if(!fTracks->UncheckedAt(iti)) continue;
156     AliTRDtrackInfo *fTrackInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iti));
157     if(!fTrackInfo->GetTrack()) continue;
158     nTracks++;
159   }
160
161   if(nTracks){
162     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks))->Fill(triggermask);
163     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtracksEvent))->Fill(nTracks);
164   }
165   if(triggermask <= 20 && !fTriggerNames->FindObject(Form("%d", triggermask))){
166     fTriggerNames->Add(new TObjString(Form("%d", triggermask)), new TObjString(triggername));
167     // also set the label for both histograms
168     TH1 *histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks));
169     histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
170     histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger));
171     histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
172   }
173   PostData(1, fContainer);
174 }
175
176
177 //_______________________________________________________
178 Bool_t AliTRDcheckDET::PostProcess(){
179   //
180   // Do Postprocessing (for the moment set the number of Reference histograms)
181   //
182   
183   TH1 * h = NULL;
184   
185   // Calculate of the trigger clusters purity
186   h = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTrigger"));
187   TH1F *h1 = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTriggerTracks"));
188   h1->Divide(h);
189   Float_t purities[20], val = 0;
190   TString triggernames[20];
191   Int_t nTriggerClasses = 0;
192   for(Int_t ibin = 1; ibin <= h->GetNbinsX(); ibin++){
193     if((val = h1->GetBinContent(ibin))){
194       purities[nTriggerClasses] = val;
195       triggernames[nTriggerClasses] = h1->GetXaxis()->GetBinLabel(ibin);
196       nTriggerClasses++;
197     }
198   }
199   h = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kTriggerPurity));
200   TAxis *ax = h->GetXaxis();
201   for(Int_t itrg = 0; itrg < nTriggerClasses; itrg++){
202     h->Fill(itrg, purities[itrg]);
203     ax->SetBinLabel(itrg+1, triggernames[itrg].Data());
204   }
205   ax->SetRangeUser(-0.5, nTriggerClasses+.5);
206   h->GetYaxis()->SetRangeUser(0,1);
207
208   // track status
209   h=dynamic_cast<TH1F*>(fContainer->At(kTrackStatus));
210   Float_t ok = h->GetBinContent(1);
211   Int_t nerr = h->GetNbinsX();
212   for(Int_t ierr=nerr; ierr--;){
213     h->SetBinContent(ierr+1, 1.e2*h->GetBinContent(ierr+1)/ok);
214   }
215   h->SetBinContent(1, 0.);
216
217   // tracklet status
218   
219   TObjArray *arr = dynamic_cast<TObjArray*>(fContainer->UncheckedAt(kTrackletStatus));
220   for(Int_t ily = AliTRDgeometry::kNlayer; ily--;){
221     h=dynamic_cast<TH1F*>(arr->At(ily));
222     Float_t okB = h->GetBinContent(1);
223     Int_t nerrB = h->GetNbinsX();
224     for(Int_t ierr=nerrB; ierr--;){
225       h->SetBinContent(ierr+1, 1.e2*h->GetBinContent(ierr+1)/okB);
226     }
227     h->SetBinContent(1, 0.);
228   }
229
230   fNRefFigures = 17;
231
232   return kTRUE;
233 }
234
235 //_______________________________________________________
236 Bool_t AliTRDcheckDET::GetRefFigure(Int_t ifig){
237   //
238   // Setting Reference Figures
239   //
240   enum FigureType_t{
241     kFigNclustersTrack,
242     kFigNclustersTracklet,
243     kFigNtrackletsTrack,
244     kFigNTrackletsP,
245     kFigNtrackletsCross,
246     kFigNtrackletsFindable,
247     kFigNtracksEvent,
248     kFigNtracksSector,
249     kFigTrackStatus,
250     kFigTrackletStatus,
251     kFigChi2,
252     kFigPH,
253     kFigChargeCluster,
254     kFigChargeTracklet,
255     kFigNeventsTrigger,
256     kFigNeventsTriggerTracks,
257     kFigTriggerPurity
258   };
259   gPad->SetLogy(0);
260   gPad->SetLogx(0);
261   TH1 *h = NULL; TObjArray *arr=NULL;
262   TLegend *leg = NULL;
263   Bool_t kFIRST(1);
264   switch(ifig){
265   case kFigNclustersTrack:
266     (h=(TH1F*)fContainer->FindObject("hNcls"))->Draw("pl");
267     PutTrendValue("NClustersTrack", h->GetMean());
268     PutTrendValue("NClustersTrackRMS", h->GetRMS());
269     return kTRUE;
270   case kFigNclustersTracklet:
271     (h =(TH1F*)fContainer->FindObject("hNclTls"))->Draw("pc");
272     PutTrendValue("NClustersTracklet", h->GetMean());
273     PutTrendValue("NClustersTrackletRMS", h->GetRMS());
274     return kTRUE;
275   case kFigNtrackletsTrack:
276     h=MakePlotNTracklets();
277     PutTrendValue("NTrackletsTrack", h->GetMean());
278     PutTrendValue("NTrackletsTrackRMS", h->GetRMS());
279     return kTRUE;
280   case kFigNTrackletsP:
281     MakePlotnTrackletsVsP();
282     return kTRUE;
283   case kFigNtrackletsCross:
284     h = (TH1F*)fContainer->FindObject("hNtlsCross");
285     if(!MakeBarPlot(h, kRed)) break;
286     PutTrendValue("NTrackletsCross", h->GetMean());
287     PutTrendValue("NTrackletsCrossRMS", h->GetRMS());
288     return kTRUE;
289   case kFigNtrackletsFindable:
290     h = (TH1F*)fContainer->FindObject("hNtlsFindable");
291     if(!MakeBarPlot(h, kGreen)) break;
292     PutTrendValue("NTrackletsFindable", h->GetMean());
293     PutTrendValue("NTrackletsFindableRMS", h->GetRMS());
294     return kTRUE;
295   case kFigNtracksEvent:
296     (h = (TH1F*)fContainer->FindObject("hNtrks"))->Draw("pl");
297     PutTrendValue("NTracksEvent", h->GetMean());
298     PutTrendValue("NTracksEventRMS", h->GetRMS());
299     return kTRUE;
300   case kFigNtracksSector:
301     h = (TH1F*)fContainer->FindObject("hNtrksSector");
302     if(!MakeBarPlot(h, kGreen)) break;
303     PutTrendValue("NTracksSector", h->Integral()/h->GetNbinsX());
304     return kTRUE;
305   case kFigTrackStatus:
306     if(!(h=(TH1F *)fContainer->FindObject("hTrackStatus"))) break;
307     h->GetXaxis()->SetRangeUser(0.5, -1);
308     h->GetYaxis()->CenterTitle();
309     h->Draw("c");
310     PutTrendValue("TrackStatus", h->Integral());
311     gPad->SetLogy(0);
312     return kTRUE;
313   case kFigTrackletStatus:
314     if(!(arr = dynamic_cast<TObjArray*>(fContainer->At(kTrackletStatus)))) break;
315     leg = new TLegend(.68, .7, .98, .98);
316     leg->SetBorderSize(1);leg->SetFillColor(0);
317     leg->SetHeader("TRD layer");
318     for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
319       if(!(h=dynamic_cast<TH1F*>(arr->At(ily)))) continue;
320       if(kFIRST){
321         h->Draw("c");
322         h->GetXaxis()->SetRangeUser(0.5, -1);
323         h->GetYaxis()->CenterTitle();
324         kFIRST = kFALSE;
325       } else h->Draw("samec");
326       leg->AddEntry(h, Form("%d", ily), "l");
327       PutTrendValue(Form("TrackletStatus%d", ily), h->Integral());
328     }
329     leg->Draw();
330     gPad->SetLogy(0);
331     return kTRUE;
332   case kFigChi2:
333     MakePlotChi2();
334     return kTRUE;
335   case kFigPH:
336     gPad->SetMargin(0.125, 0.015, 0.1, 0.1);
337     MakePlotPulseHeight();
338     gPad->SetLogy(0);
339     return kTRUE;
340   case kFigChargeCluster:
341     (h = (TH1F*)fContainer->FindObject("hQcl"))->Draw("c");
342     gPad->SetLogy(1);
343     PutTrendValue("ChargeCluster", h->GetMaximumBin());
344     PutTrendValue("ChargeClusterRMS", h->GetRMS());
345     return kTRUE;
346   case kFigChargeTracklet:
347     (h=(TH1F*)fContainer->FindObject("hQtrklt"))->Draw("c");
348     PutTrendValue("ChargeTracklet", h->GetMaximumBin());
349     PutTrendValue("ChargeTrackletRMS", h->GetRMS());
350     return kTRUE;
351   case kFigNeventsTrigger:
352     ((TH1F*)fContainer->FindObject("hEventsTrigger"))->Draw("");
353     return kTRUE;
354   case kFigNeventsTriggerTracks:
355     ((TH1F*)fContainer->FindObject("hEventsTriggerTracks"))->Draw("");
356     return kTRUE;
357   case kFigTriggerPurity: 
358     if(!MakeBarPlot((TH1F*)fContainer->FindObject("hTriggerPurity"), kGreen)) break;
359     break;
360   default:
361     break;
362   }
363   AliInfo(Form("Reference plot [%d] missing result", ifig));
364   return kFALSE;
365 }
366
367 //_______________________________________________________
368 TObjArray *AliTRDcheckDET::Histos(){
369   //
370   // Create QA histograms
371   //
372     
373   if(fContainer) return fContainer;
374   
375   fContainer = new TObjArray(20);
376   //fContainer->SetOwner(kTRUE);
377
378   // Register Histograms
379   TH1 * h = NULL;
380   TAxis *ax = NULL;
381   if(!(h = (TH1F *)gROOT->FindObject("hNcls"))){
382     h = new TH1F("hNcls", "N_{clusters} / track", 181, -0.5, 180.5);
383     h->GetXaxis()->SetTitle("N_{clusters}");
384     h->GetYaxis()->SetTitle("Entries");
385   } else h->Reset();
386   fContainer->AddAt(h, kNclustersTrack);
387
388   if(!(h = (TH1F *)gROOT->FindObject("hNclTls"))){
389     h = new TH1F("hNclTls","N_{clusters} / tracklet", 51, -0.5, 50.5);
390     h->GetXaxis()->SetTitle("N_{clusters}");
391     h->GetYaxis()->SetTitle("Entries");
392   } else h->Reset();
393   fContainer->AddAt(h, kNclustersTracklet);
394
395   if(!(h = (TH1F *)gROOT->FindObject("hNtls"))){
396     h = new TH1F("hNtls", "N_{tracklets} / track", AliTRDgeometry::kNlayer, 0.5, 6.5);
397     h->GetXaxis()->SetTitle("N^{tracklet}");
398     h->GetYaxis()->SetTitle("freq. [%]");
399   } else h->Reset();
400   fContainer->AddAt(h, kNtrackletsTrack);
401
402   if(!(h = (TH1F *)gROOT->FindObject("htlsSTA"))){
403     h = new TH1F("hNtlsSTA", "N_{tracklets} / track (Stand Alone)", AliTRDgeometry::kNlayer, 0.5, 6.5);
404     h->GetXaxis()->SetTitle("N^{tracklet}");
405     h->GetYaxis()->SetTitle("freq. [%]");
406   }
407   fContainer->AddAt(h, kNtrackletsSTA);
408
409   // Binning for momentum dependent tracklet Plots
410   const Int_t kNPbins = 11;
411   Double_t binTracklets[AliTRDgeometry::kNlayer + 1];
412   for(Int_t il = 0; il <= AliTRDgeometry::kNlayer; il++) binTracklets[il] = 0.5 + il;
413   Double_t binMomentum[kNPbins + 1] = {0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.5, 2., 3., 4., 5., 10};
414
415   if(!(h = (TH1F *)gROOT->FindObject("htlsBAR"))){
416     // Make tracklets for barrel tracks momentum dependent (if we do not exceed min and max values)
417     h = new TH2F("hNtlsBAR", "N_{tracklets} / track (Barrel)", kNPbins, binMomentum, AliTRDgeometry::kNlayer, binTracklets);
418     h->GetXaxis()->SetTitle("p / GeV/c");
419     h->GetYaxis()->SetTitle("N^{tracklet}");
420     h->GetZaxis()->SetTitle("freq. [%]");
421   }
422   fContainer->AddAt(h, kNtrackletsBAR);
423
424   // 
425   if(!(h = (TH1F *)gROOT->FindObject("hNtlsCross"))){
426     h = new TH1F("hNtlsCross", "N_{tracklets}^{cross} / track", 7, -0.5, 6.5);
427     h->GetXaxis()->SetTitle("n_{row cross}");
428     h->GetYaxis()->SetTitle("freq. [%]");
429   } else h->Reset();
430   fContainer->AddAt(h, kNtrackletsCross);
431
432   if(!(h = (TH1F *)gROOT->FindObject("hNtlsFindable"))){
433     h = new TH1F("hNtlsFindable", "Found/Findable Tracklets" , 101, -0.005, 1.005);
434     h->GetXaxis()->SetTitle("r [a.u]");
435     h->GetYaxis()->SetTitle("Entries");
436   } else h->Reset();
437   fContainer->AddAt(h, kNtrackletsFindable);
438
439   if(!(h = (TH1F *)gROOT->FindObject("hNtrks"))){
440     h = new TH1F("hNtrks", "N_{tracks} / event", 100, 0, 100);
441     h->GetXaxis()->SetTitle("N_{tracks}");
442     h->GetYaxis()->SetTitle("Entries");
443   } else h->Reset();
444   fContainer->AddAt(h, kNtracksEvent);
445
446   if(!(h = (TH1F *)gROOT->FindObject("hNtrksSector"))){
447     h = new TH1F("hNtrksSector", "N_{tracks} / sector", AliTRDgeometry::kNsector, -0.5, 17.5);
448     h->GetXaxis()->SetTitle("sector");
449     h->GetYaxis()->SetTitle("freq. [%]");
450   } else h->Reset();
451   fContainer->AddAt(h, kNtracksSector);
452
453   if(!(h = (TH1F*)gROOT->FindObject("hTrackStatus"))){
454     const Int_t nerr = 7;
455     h = new TH1F("hTrackStatus", "Track Status", nerr, -0.5, nerr-0.5);
456     const Char_t *label[nerr] = {"OK", "PROL", "PROP", "AJST", "SNP", "TINI", "UPDT"};
457     ax = h->GetXaxis();
458     for(Int_t ierr=nerr; ierr--;) ax->SetBinLabel(ierr+1, label[ierr]);
459     h->SetYTitle("Relative Error to Good [%]");
460   }
461   fContainer->AddAt(h, kTrackStatus);
462
463   TObjArray *arr = new TObjArray(AliTRDgeometry::kNlayer);
464   arr->SetOwner(kTRUE);  arr->SetName("TrackletStatus");
465   fContainer->AddAt(arr, kTrackletStatus);
466   for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
467     if(!(h = (TH1F *)gROOT->FindObject(Form("hTrackletStatus%d", ily)))){
468       const Int_t nerr = 8;
469       h = new TH1F(Form("hTrackletStatus%d", ily), "Tracklet status", nerr, -0.5, nerr-0.5);
470       h->SetLineColor(ily+1);
471       const Char_t *label[nerr] = {"OK", "Geom", "Bound", "NoCl", "NoAttach", "NoClTr", "NoFit", "Chi2"};
472       ax = h->GetXaxis();
473       for(Int_t ierr=nerr; ierr--;) ax->SetBinLabel(ierr+1, label[ierr]);
474       h->SetYTitle("Relative Error to Good [%]");
475     } else h->Reset();
476     arr->AddAt(h, ily);
477   }
478
479   // <PH> histos
480   arr = new TObjArray(2);
481   arr->SetOwner(kTRUE);  arr->SetName("<PH>");
482   fContainer->AddAt(arr, kPH);
483   if(!(h = (TH1F *)gROOT->FindObject("hPHt"))){
484     h = new TProfile("hPHt", "<PH>", 31, -0.5, 30.5);
485     h->GetXaxis()->SetTitle("Time / 100ns");
486     h->GetYaxis()->SetTitle("<PH> [a.u]");
487   } else h->Reset();
488   arr->AddAt(h, 0);
489   if(!(h = (TH1F *)gROOT->FindObject("hPHx")))
490     h = new TProfile("hPHx", "<PH>", 31, -0.08, 4.88);
491   else h->Reset();
492   arr->AddAt(h, 1);
493
494   // Chi2 histos
495   if(!(h = (TH2S*)gROOT->FindObject("hChi2"))){
496     h = new TH2S("hChi2", "#chi^{2} per track", AliTRDgeometry::kNlayer, .5, AliTRDgeometry::kNlayer+.5, 100, 0, 50);
497     h->SetXTitle("ndf");
498     h->SetYTitle("#chi^{2}/ndf");
499     h->SetZTitle("entries");
500   } else h->Reset();
501   fContainer->AddAt(h, kChi2);
502
503   if(!(h = (TH1F *)gROOT->FindObject("hQcl"))){
504     h = new TH1F("hQcl", "Q_{cluster}", 200, 0, 1200);
505     h->GetXaxis()->SetTitle("Q_{cluster} [a.u.]");
506     h->GetYaxis()->SetTitle("Entries");
507   }else h->Reset();
508   fContainer->AddAt(h, kChargeCluster);
509
510   if(!(h = (TH1F *)gROOT->FindObject("hQtrklt"))){
511     h = new TH1F("hQtrklt", "Q_{tracklet}", 6000, 0, 6000);
512     h->GetXaxis()->SetTitle("Q_{tracklet} [a.u.]");
513     h->GetYaxis()->SetTitle("Entries");
514   }else h->Reset();
515   fContainer->AddAt(h, kChargeTracklet);
516
517
518   if(!(h = (TH1F *)gROOT->FindObject("hEventsTrigger")))
519     h = new TH1F("hEventsTrigger", "Trigger Class", 100, 0, 100);
520   else h->Reset();
521   printf("Histos Adding \n");
522   
523   fContainer->AddAt(h, kNeventsTrigger);
524
525   if(!(h = (TH1F *)gROOT->FindObject("hEventsTriggerTracks")))
526     h = new TH1F("hEventsTriggerTracks", "Trigger Class (Tracks)", 100, 0, 100);
527   else h->Reset();
528   fContainer->AddAt(h, kNeventsTriggerTracks);
529
530   if(!(h = (TH1F *)gROOT->FindObject("hTriggerPurity"))){
531     h = new TH1F("hTriggerPurity", "Trigger Purity", 10, -0.5, 9.5);
532     h->GetXaxis()->SetTitle("Trigger Cluster");
533     h->GetYaxis()->SetTitle("freq.");
534   } else h->Reset();
535   fContainer->AddAt(h, kTriggerPurity);
536
537   return fContainer;
538 }
539
540 /*
541 * Plotting Functions
542 */
543
544 //_______________________________________________________
545 TH1 *AliTRDcheckDET::PlotTrackStatus(const AliTRDtrackV1 *track)
546 {
547 //
548 // Plot the track propagation status. The following errors are defined (see AliTRDtrackV1::ETRDtrackError)
549 //   PROL - track prolongation failure
550 //   PROP - track propagation failure
551 //   AJST - crossing sectors failure
552 //   SNP  - too large bending
553 //   TINI - tracklet initialization failure
554 //   UPDT - track position/covariance update failure 
555 //
556 // Performance plot looks as below:
557 //Begin_Html
558 //<img src="TRD/trackStatus.gif">
559 //End_Html
560 //
561   if(track) fkTrack = track;
562   if(!fkTrack){
563     AliWarning("No Track defined.");
564     return NULL;
565   }
566   TH1 *h = NULL;
567   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kTrackStatus)))){
568     AliWarning("No Histogram defined.");
569     return NULL;
570   }
571   h->Fill(fkTrack->GetStatusTRD());
572   return h;
573 }
574
575 //_______________________________________________________
576 TH1 *AliTRDcheckDET::PlotTrackletStatus(const AliTRDtrackV1 *track)
577 {
578 //
579 // Plot the tracklet propagation status. The following errors are defined for tracklet (see AliTRDtrackV1::ETRDlayerError)
580 //   Geom   - 
581 //   Bound  - tracklet too close to chamber walls
582 //   NoCl   - no clusters in the track roads
583 //   NoAttach - fail to attach clusters
584 //   NoClTr - fail to use clusters for fit
585 //   NoFit  - tracklet fit failled
586 //   Chi2   - chi2 tracklet-track over threshold
587 //
588 // Performance plot looks as below:
589 //Begin_Html
590 //<img src="TRD/trackletStatus.gif">
591 //End_Html
592 //
593   if(track) fkTrack = track;
594   if(!fkTrack){
595     AliWarning("No Track defined.");
596     return NULL;
597   }
598   TObjArray *arr =NULL;
599   if(!(arr = dynamic_cast<TObjArray*>(fContainer->At(kTrackletStatus)))){
600     AliWarning("Histograms not defined.");
601     return NULL;
602   }
603
604   TH1 *h = NULL;
605   for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
606     if(!(h = dynamic_cast<TH1F*>(arr->At(ily)))){
607       AliWarning(Form("Missing histo for layer %d.", ily));
608       continue;
609     }
610     h->Fill(fkTrack->GetStatusTRD(ily));
611   }
612   return h;
613 }
614
615 //_______________________________________________________
616 TH1 *AliTRDcheckDET::PlotNClustersTracklet(const AliTRDtrackV1 *track){
617   //
618   // Plot the mean number of clusters per tracklet
619   //
620   if(track) fkTrack = track;
621   if(!fkTrack){
622     AliWarning("No Track defined.");
623     return NULL;
624   }
625   TH1 *h = NULL;
626   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclustersTracklet)))){
627     AliWarning("No Histogram defined.");
628     return NULL;
629   }
630   AliTRDseedV1 *tracklet = NULL;
631   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
632     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
633     h->Fill(tracklet->GetN2());
634   }
635   return h;
636 }
637
638 //_______________________________________________________
639 TH1 *AliTRDcheckDET::PlotNClustersTrack(const AliTRDtrackV1 *track){
640   //
641   // Plot the number of clusters in one track
642   //
643   if(track) fkTrack = track;
644   if(!fkTrack){
645     AliWarning("No Track defined.");
646     return NULL;
647   }
648   TH1 *h = NULL;
649   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclustersTrack)))){
650     AliWarning("No Histogram defined.");
651     return NULL;
652   }
653   
654   Int_t nclusters = 0;
655   AliTRDseedV1 *tracklet = NULL;
656   AliExternalTrackParam *par = fkTrack->GetTrackOut() ? fkTrack->GetTrackOut() : fkTrack->GetTrackIn();
657   Double_t momentumRec = par->P();
658   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
659     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
660     Int_t n(tracklet->GetN());
661     nclusters += n;
662     if(DebugLevel() > 2){
663       Int_t crossing = Int_t(tracklet->IsRowCross());
664       Int_t detector = tracklet->GetDetector();
665       Float_t theta = TMath::ATan(tracklet->GetZref(1));
666       Float_t phi = TMath::ATan(tracklet->GetYref(1));
667       Float_t momentumMC = 0.;
668       Int_t pdg = 0;
669       Int_t kinkIndex = fkESD ? fkESD->GetKinkIndex() : 0;
670       UShort_t nclsTPC = fkESD ? fkESD->GetTPCncls() : 0;
671       if(fkMC){
672         if(fkMC->GetTrackRef()) momentumMC = fkMC->GetTrackRef()->P();
673         pdg = fkMC->GetPDG();
674       }
675       (*DebugStream()) << "NClustersTrack"
676         << "Detector="  << detector
677         << "crossing="  << crossing
678         << "momentumMC="<< momentumMC
679         << "momentumRec=" << momentumRec
680         << "pdg="                               << pdg
681         << "theta="                     << theta
682         << "phi="                               << phi
683         << "kinkIndex=" << kinkIndex
684         << "TPCncls="           << nclsTPC
685         << "TRDncls="   << n
686         << "\n";
687     }
688   }
689   h->Fill(nclusters);
690   return h;
691 }
692
693
694 //_______________________________________________________
695 TH1 *AliTRDcheckDET::PlotNTrackletsTrack(const AliTRDtrackV1 *track){
696   //
697   // Plot the number of tracklets
698   //
699   if(track) fkTrack = track;
700   if(!fkTrack){
701     AliWarning("No Track defined.");
702     return NULL;
703   }
704   TH1 *h = NULL, *hSta = NULL; TH2 *hBarrel = NULL;
705   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsTrack)))){
706     AliWarning("No Histogram defined.");
707     return NULL;
708   }
709   Int_t nTracklets = fkTrack->GetNumberOfTracklets();
710   h->Fill(nTracklets);
711   if(!fkESD) return h;
712   Int_t status = fkESD->GetStatus();
713
714 /*  printf("in/out/refit/pid: TRD[%d|%d|%d|%d]\n", status &AliESDtrack::kTRDin ? 1 : 0, status &AliESDtrack::kTRDout ? 1 : 0, status &AliESDtrack::kTRDrefit ? 1 : 0, status &AliESDtrack::kTRDpid ? 1 : 0);*/
715   Double_t p = 0.;
716   Int_t method = -1;    // to distinguish between stand alone and full barrel tracks in the debugging
717   if((status & AliESDtrack::kTRDin) != 0){
718     method = 1;
719     // Full Barrel Track: Save momentum dependence
720     if(!(hBarrel = dynamic_cast<TH2F *>(fContainer->At(kNtrackletsBAR)))){
721       AliWarning("Method: Barrel.  Histogram not processed!");
722     } else {
723       AliExternalTrackParam *par = fkTrack->GetTrackOut() ? fkTrack->GetTrackOut() : fkTrack->GetTrackIn();
724       if(!par){
725        AliError("Outer track params missing");
726       } else {
727         p = par->P();
728       }
729       hBarrel->Fill(p, nTracklets);
730     }
731   } else {
732     // Stand alone Track: momentum dependence not usefull
733     method = 0;
734     if(!(hSta = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsSTA)))) {
735       AliWarning("Method: StandAlone.  Histogram not processed!");
736     } else {
737       hSta->Fill(nTracklets);
738     }
739   }
740
741   if(DebugLevel() > 2){
742     AliTRDseedV1 *tracklet = NULL;
743     Int_t sector = -1;
744     for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
745       if(!(tracklet = fkTrack->GetTracklet(itl)) || !(tracklet->IsOK())) continue;
746       sector = fGeo->GetSector(tracklet->GetDetector());
747       break;
748     }
749     (*DebugStream()) << "NTrackletsTrack"
750       << "Sector="      << sector
751       << "NTracklets="  << nTracklets
752       << "Method="      << method
753       << "p="           << p
754       << "\n";
755   }
756   if(DebugLevel() > 3){
757     if(nTracklets == 1){
758       // If we have one Tracklet, check in which layer this happens
759       Int_t layer = -1;
760       AliTRDseedV1 *tracklet = NULL;
761       for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
762         if((tracklet = fkTrack->GetTracklet(il)) && tracklet->IsOK()){layer =  il; break;}
763       }
764       if(layer >= 0){
765         (*DebugStream()) << "NTrackletsLayer"
766           << "Layer=" << layer
767           << "p=" << p
768           << "\n";
769       }
770     }
771   }
772   return h;
773 }
774
775
776 //_______________________________________________________
777 TH1 *AliTRDcheckDET::PlotNTrackletsRowCross(const AliTRDtrackV1 *track){
778   //
779   // Plot the number of tracklets
780   //
781   if(track) fkTrack = track;
782   if(!fkTrack){
783     AliWarning("No Track defined.");
784     return NULL;
785   }
786   TH1 *h = NULL;
787   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsCross)))){
788     AliWarning("No Histogram defined.");
789     return NULL;
790   }
791
792   Int_t ncross = 0;
793   AliTRDseedV1 *tracklet = NULL;
794   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
795     if(!(tracklet = fkTrack->GetTracklet(il)) || !tracklet->IsOK()) continue;
796
797     if(tracklet->IsRowCross()) ncross++;
798   }
799   h->Fill(ncross);
800   return h;
801 }
802
803 //_______________________________________________________
804 TH1 *AliTRDcheckDET::PlotFindableTracklets(const AliTRDtrackV1 *track){
805   //
806   // Plots the ratio of number of tracklets vs.
807   // number of findable tracklets
808   //
809   // Findable tracklets are defined as track prolongation
810   // to layer i does not hit the dead area +- epsilon
811   //
812   // In order to check whether tracklet hist active area in Layer i, 
813   // the track is refitted and the fitted position + an uncertainty 
814   // range is compared to the chamber border (also with a different
815   // uncertainty)
816   //
817   // For the track fit two cases are distinguished:
818   // If the track is a stand alone track (defined by the status bit 
819   // encoding, then the track is fitted with the tilted Rieman model
820   // Otherwise the track is fitted with the Kalman fitter in two steps:
821   // Since the track parameters are give at the outer point, we first 
822   // fit in direction inwards. Afterwards we fit again in direction outwards
823   // to extrapolate the track to layers which are not reached by the track
824   // For the Kalman model, the radial track points have to be shifted by
825   // a distance epsilon in the direction that we want to fit
826   //
827   const Float_t epsilon = 0.01;   // dead area tolerance
828   const Float_t epsilonR = 1;    // shift in radial direction of the anode wire position (Kalman filter only)
829   const Float_t deltaY = 0.7;    // Tolerance in the track position in y-direction
830   const Float_t deltaZ = 7.0;    // Tolerance in the track position in z-direction (Padlength)
831   Double_t xAnode[AliTRDgeometry::kNlayer] = {300.2, 312.8, 325.4, 338.0, 350.6, 363.2}; // Take the default X0
832  
833   if(track) fkTrack = track;
834   if(!fkTrack){
835     AliWarning("No Track defined.");
836     return NULL;
837   }
838   TH1 *h = NULL;
839   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsFindable)))){
840     AliWarning("No Histogram defined.");
841     return NULL;
842   }
843   Int_t nFound = 0, nFindable = 0;
844   Int_t stack = -1;
845   Double_t ymin = 0., ymax = 0., zmin = 0., zmax = 0.;
846   Double_t y = 0., z = 0.;
847   AliTRDseedV1 *tracklet = NULL;
848   AliTRDpadPlane *pp;  
849   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
850     if((tracklet = fkTrack->GetTracklet(il)) && tracklet->IsOK()){
851       tracklet->SetReconstructor(fReconstructor);
852       nFound++;
853     }
854   }
855   // 2 Different cases:
856   // 1st stand alone: here we cannot propagate, but be can do a Tilted Rieman Fit
857   // 2nd barrel track: here we propagate the track to the layers
858   AliTrackPoint points[6];
859   Float_t xyz[3];
860   memset(xyz, 0, sizeof(Float_t) * 3);
861   if(((fkESD->GetStatus() & AliESDtrack::kTRDout) > 0) && !((fkESD->GetStatus() & AliESDtrack::kTRDin) > 0)){
862     // stand alone track
863     for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
864       xyz[0] = xAnode[il];
865       points[il].SetXYZ(xyz);
866     }
867     AliTRDtrackerV1::FitRiemanTilt(const_cast<AliTRDtrackV1 *>(fkTrack), NULL, kTRUE, 6, points);
868   } else {
869     // barrel track
870     //
871     // 2 Steps:
872     // -> Kalman inwards
873     // -> Kalman outwards
874     AliTRDtrackV1 copyTrack(*fkTrack);  // Do Kalman on a (non-constant) copy of the track
875     AliTrackPoint pointsInward[6], pointsOutward[6];
876     for(Int_t il = AliTRDgeometry::kNlayer; il--;){
877       // In order to avoid complications in the Kalman filter if the track points have the same radial
878       // position like the tracklets, we have to shift the radial postion of the anode wire by epsilon
879       // in the direction we want to go
880       // The track points have to be in reverse order for the Kalman Filter inwards
881       xyz[0] = xAnode[AliTRDgeometry::kNlayer - il - 1] - epsilonR;
882       pointsInward[il].SetXYZ(xyz);
883       xyz[0] = xAnode[il] + epsilonR;
884       pointsOutward[il].SetXYZ(xyz);
885     }
886     /*for(Int_t ipt = 0; ipt < AliTRDgeometry::kNlayer; ipt++)
887       printf("%d. X = %f\n", ipt, points[ipt].GetX());*/
888     // Kalman inwards
889     AliTRDtrackerV1::FitKalman(&copyTrack, NULL, kFALSE, 6, pointsInward);
890     memcpy(points, pointsInward, sizeof(AliTrackPoint) * 6); // Preliminary store the inward results in the Array points
891     // Kalman outwards
892     AliTRDtrackerV1::FitKalman(&copyTrack, NULL, kTRUE, 6, pointsInward);
893     memcpy(points, pointsOutward, sizeof(AliTrackPoint) * AliTRDgeometry::kNlayer);
894   }
895   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
896     y = points[il].GetY();
897     z = points[il].GetZ();
898     if((stack = fGeo->GetStack(z, il)) < 0) continue; // Not findable
899     pp = fGeo->GetPadPlane(il, stack);
900     ymin = pp->GetCol0() + epsilon;
901     ymax = pp->GetColEnd() - epsilon; 
902     zmin = pp->GetRowEnd() + epsilon; 
903     zmax = pp->GetRow0() - epsilon;
904     // ignore y-crossing (material)
905     if((z + deltaZ > zmin && z - deltaZ < zmax) && (y + deltaY > ymin && y - deltaY < ymax)) nFindable++;
906       if(DebugLevel() > 3){
907         Double_t posTracklet[2] = {tracklet ? tracklet->GetYfit(0) : 0, tracklet ? tracklet->GetZfit(0) : 0};
908         Int_t hasTracklet = tracklet ? 1 : 0;
909         (*DebugStream())   << "FindableTracklets"
910           << "layer="     << il
911           << "ytracklet=" << posTracklet[0]
912           << "ytrack="    << y
913           << "ztracklet=" << posTracklet[1]
914           << "ztrack="    << z
915           << "tracklet="  << hasTracklet
916           << "\n";
917       }
918   }
919   
920   h->Fill(nFindable > 0 ? TMath::Min(nFound/static_cast<Double_t>(nFindable), 1.) : 1);
921   AliDebug(2, Form("Findable[Found]: %d[%d|%f]", nFindable, nFound, nFound/static_cast<Float_t>(nFindable > 0 ? nFindable : 1)));
922   return h;
923 }
924
925
926 //_______________________________________________________
927 TH1 *AliTRDcheckDET::PlotChi2(const AliTRDtrackV1 *track){
928   //
929   // Plot the chi2 of the track
930   //
931   if(track) fkTrack = track;
932   if(!fkTrack){
933     AliWarning("No Track defined.");
934     return NULL;
935   }
936   TH1 *h = NULL;
937   if(!(h = dynamic_cast<TH2S*>(fContainer->At(kChi2)))) {
938     AliWarning("No Histogram defined.");
939     return NULL;
940   }
941   Int_t n = fkTrack->GetNumberOfTracklets();
942   if(!n) return NULL;
943
944   h->Fill(n, fkTrack->GetChi2()/n);
945   return h;
946 }
947
948
949 //_______________________________________________________
950 TH1 *AliTRDcheckDET::PlotPHt(const AliTRDtrackV1 *track){
951   //
952   // Plot the average pulse height
953   //
954   if(track) fkTrack = track;
955   if(!fkTrack){
956     AliWarning("No Track defined.");
957     return NULL;
958   }
959   TProfile *h = NULL;
960   if(!(h = dynamic_cast<TProfile *>(((TObjArray*)(fContainer->At(kPH)))->At(0)))){
961     AliWarning("No Histogram defined.");
962     return NULL;
963   }
964   AliTRDseedV1 *tracklet = NULL;
965   AliTRDcluster *c = NULL;
966   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
967     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
968     Int_t crossing = Int_t(tracklet->IsRowCross());
969     Int_t detector = tracklet->GetDetector();
970     tracklet->ResetClusterIter();
971     while((c = tracklet->NextCluster())){
972       if(!IsUsingClustersOutsideChamber() && !c->IsInChamber()) continue;
973       Int_t localtime        = c->GetLocalTimeBin();
974       Double_t absoluteCharge = TMath::Abs(c->GetQ());
975       h->Fill(localtime, absoluteCharge);
976       if(DebugLevel() > 3){
977         Double_t distance[2];
978         GetDistanceToTracklet(distance, tracklet, c);
979         Float_t theta = TMath::ATan(tracklet->GetZref(1));
980         Float_t phi = TMath::ATan(tracklet->GetYref(1));
981         Float_t momentum = 0.;
982         Int_t pdg = 0;
983         Int_t kinkIndex = fkESD ? fkESD->GetKinkIndex() : 0;
984         UShort_t TPCncls = fkESD ? fkESD->GetTPCncls() : 0;
985         if(fkMC){
986           if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
987           pdg = fkMC->GetPDG();
988         }
989         (*DebugStream()) << "PHt"
990           << "Detector="        << detector
991           << "crossing="        << crossing
992           << "Timebin="         << localtime
993           << "Charge="          << absoluteCharge
994           << "momentum="        << momentum
995           << "pdg="                             << pdg
996           << "theta="                   << theta
997           << "phi="                             << phi
998           << "kinkIndex="       << kinkIndex
999           << "TPCncls="         << TPCncls
1000           << "dy="        << distance[0]
1001           << "dz="        << distance[1]
1002           << "c.="        << c
1003           << "\n";
1004       }
1005     }
1006   }
1007   return h;
1008 }
1009
1010 //_______________________________________________________
1011 TH1 *AliTRDcheckDET::PlotPHx(const AliTRDtrackV1 *track){
1012   //
1013   // Plots the average pulse height vs the distance from the anode wire
1014   // (plus const anode wire offset)
1015   //
1016   if(track) fkTrack = track;
1017   if(!fkTrack){
1018     AliWarning("No Track defined.");
1019     return NULL;
1020   }
1021   TProfile *h = NULL;
1022   if(!(h = dynamic_cast<TProfile *>(((TObjArray*)(fContainer->At(kPH)))->At(1)))){
1023     AliWarning("No Histogram defined.");
1024     return NULL;
1025   }
1026   Float_t offset = .5*AliTRDgeometry::CamHght();
1027   AliTRDseedV1 *tracklet = NULL;
1028   AliTRDcluster *c = NULL;
1029   Double_t distance = 0;
1030   Double_t x, y;
1031   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
1032     if(!(tracklet = fkTrack->GetTracklet(itl)) || !(tracklet->IsOK())) continue;
1033     tracklet->ResetClusterIter();
1034     while((c = tracklet->NextCluster())){
1035       if(!IsUsingClustersOutsideChamber() && !c->IsInChamber()) continue;
1036       x = c->GetX()-AliTRDcluster::GetXcorr(c->GetLocalTimeBin());
1037       y = c->GetY()-AliTRDcluster::GetYcorr(AliTRDgeometry::GetLayer(c->GetDetector()), c->GetCenter());
1038
1039       distance = tracklet->GetX0() - (c->GetX() + 0.3) + offset;
1040       h->Fill(distance, TMath::Abs(c->GetQ()));
1041     }
1042   }  
1043   return h;
1044 }
1045
1046 //_______________________________________________________
1047 TH1 *AliTRDcheckDET::PlotChargeCluster(const AliTRDtrackV1 *track){
1048   //
1049   // Plot the cluster charge
1050   //
1051   if(track) fkTrack = track;
1052   if(!fkTrack){
1053     AliWarning("No Track defined.");
1054     return NULL;
1055   }
1056   TH1 *h = NULL;
1057   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeCluster)))){
1058     AliWarning("No Histogram defined.");
1059     return NULL;
1060   }
1061   AliTRDseedV1 *tracklet = NULL;
1062   AliTRDcluster *c = NULL;
1063   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
1064     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
1065     for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
1066       if(!(c = tracklet->GetClusters(itime))) continue;
1067       h->Fill(c->GetQ());
1068     }
1069   }
1070   return h;
1071 }
1072
1073 //_______________________________________________________
1074 TH1 *AliTRDcheckDET::PlotChargeTracklet(const AliTRDtrackV1 *track){
1075   //
1076   // Plot the charge deposit per chamber
1077   //
1078   if(track) fkTrack = track;
1079   if(!fkTrack){
1080     AliWarning("No Track defined.");
1081     return NULL;
1082   }
1083   TH1 *h = NULL;
1084   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeTracklet)))){
1085     AliWarning("No Histogram defined.");
1086     return NULL;
1087   }
1088   AliTRDseedV1 *tracklet = NULL;
1089   AliTRDcluster *c = NULL;
1090   Double_t qTot = 0;
1091   Int_t nTracklets =fkTrack->GetNumberOfTracklets();
1092   for(Int_t itl(0); itl < AliTRDgeometry::kNlayer; itl++){
1093     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
1094     qTot = 0.;
1095     for(Int_t ic = AliTRDseedV1::kNclusters; ic--;){
1096       if(!(c = tracklet->GetClusters(ic))) continue;
1097       qTot += TMath::Abs(c->GetQ());
1098     }
1099     h->Fill(qTot);
1100     if(DebugLevel() > 3){
1101       Int_t crossing = (Int_t)tracklet->IsRowCross();
1102       Int_t detector = tracklet->GetDetector();
1103       Float_t theta = TMath::ATan(tracklet->GetZfit(1));
1104       Float_t phi = TMath::ATan(tracklet->GetYfit(1));
1105       Float_t momentum = 0.;
1106       Int_t pdg = 0;
1107       Int_t kinkIndex = fkESD ? fkESD->GetKinkIndex() : 0;
1108       UShort_t nclsTPC = fkESD ? fkESD->GetTPCncls() : 0;
1109       if(fkMC){
1110               if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
1111         pdg = fkMC->GetPDG();
1112       }
1113       (*DebugStream()) << "ChargeTracklet"
1114         << "Detector="  << detector
1115         << "crossing="  << crossing
1116         << "momentum="  << momentum
1117         << "nTracklets="<< nTracklets
1118         << "pdg="                               << pdg
1119         << "theta="                     << theta
1120         << "phi="                               << phi
1121         << "kinkIndex=" << kinkIndex
1122         << "TPCncls="           << nclsTPC
1123         << "QT="        << qTot
1124         << "\n";
1125     }
1126   }
1127   return h;
1128 }
1129
1130 //_______________________________________________________
1131 TH1 *AliTRDcheckDET::PlotNTracksSector(const AliTRDtrackV1 *track){
1132   //
1133   // Plot the number of tracks per Sector
1134   //
1135   if(track) fkTrack = track;
1136   if(!fkTrack){
1137     AliWarning("No Track defined.");
1138     return NULL;
1139   }
1140   TH1 *h = NULL;
1141   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtracksSector)))){
1142     AliWarning("No Histogram defined.");
1143     return NULL;
1144   }
1145
1146   // TODO we should compare with
1147   // sector = Int_t(track->GetAlpha() / AliTRDgeometry::GetAlpha());
1148
1149   AliTRDseedV1 *tracklet = NULL;
1150   Int_t sector = -1;
1151   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
1152     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
1153     sector = static_cast<Int_t>(tracklet->GetDetector()/AliTRDgeometry::kNdets);
1154     break;
1155   }
1156   h->Fill(sector);
1157   return h;
1158 }
1159
1160
1161 //________________________________________________________
1162 void AliTRDcheckDET::SetRecoParam(AliTRDrecoParam *r)
1163 {
1164
1165   fReconstructor->SetRecoParam(r);
1166 }
1167
1168 //________________________________________________________
1169 void AliTRDcheckDET::GetDistanceToTracklet(Double_t *dist, AliTRDseedV1 * const tracklet, AliTRDcluster * const c)
1170 {
1171   Float_t x = c->GetX();
1172   dist[0] = c->GetY() - tracklet->GetYat(x);
1173   dist[1] = c->GetZ() - tracklet->GetZat(x);
1174 }
1175
1176
1177 //_______________________________________________________
1178 TH1* AliTRDcheckDET::MakePlotChi2()
1179 {
1180 // Plot chi2/track normalized to number of degree of freedom 
1181 // (tracklets) and compare with the theoretical distribution.
1182 // 
1183 // Alex Bercuci <A.Bercuci@gsi.de>
1184
1185   TH2S *h2 = (TH2S*)fContainer->At(kChi2);
1186   TF1 f("fChi2", "[0]*pow(x, [1]-1)*exp(-0.5*x)", 0., 50.);
1187   TLegend *leg = new TLegend(.7,.7,.95,.95);
1188   leg->SetBorderSize(1); leg->SetHeader("Tracklets per Track");
1189   TH1D *h1 = NULL;
1190   Bool_t kFIRST = kTRUE;
1191   for(Int_t il=1; il<=h2->GetNbinsX(); il++){
1192     h1 = h2->ProjectionY(Form("pyChi2%d", il), il, il);
1193     if(h1->Integral()<50) continue;
1194     h1->Scale(1./h1->Integral());
1195     h1->SetMarkerStyle(7);h1->SetMarkerColor(il);
1196     h1->SetLineColor(il);h1->SetLineStyle(2);
1197     f.SetParameter(1, .5*il);f.SetLineColor(il);
1198     h1->Fit(&f, "QW+", kFIRST ? "pc": "pcsame");
1199     leg->AddEntry(h1, Form("%d", il), "l");
1200     if(kFIRST){
1201       h1->GetXaxis()->SetRangeUser(0., 25.);
1202     }
1203     kFIRST = kFALSE;
1204   }
1205   leg->Draw();
1206   gPad->SetLogy();
1207   return h1;
1208 }
1209
1210
1211 //________________________________________________________
1212 TH1* AliTRDcheckDET::MakePlotNTracklets(){
1213   //
1214   // Make nice bar plot of the number of tracklets in each method
1215   //
1216   TH2F *tmp = (TH2F *)fContainer->FindObject("hNtlsBAR");
1217   TH1D *hBAR = tmp->ProjectionY();
1218   TH1F *hSTA = (TH1F *)fContainer->FindObject("hNtlsSTA");
1219   TH1F *hCON = (TH1F *)fContainer->FindObject("hNtls");
1220   TLegend *leg = new TLegend(0.13, 0.75, 0.39, 0.89);
1221   leg->SetBorderSize(1);
1222   leg->SetFillColor(0);
1223
1224   Float_t scale = hCON->Integral();
1225   hCON->Scale(100./scale);
1226   hCON->SetFillColor(kRed);hCON->SetLineColor(kRed);
1227   hCON->SetBarWidth(0.2);
1228   hCON->SetBarOffset(0.6);
1229   hCON->SetStats(kFALSE);
1230   hCON->GetYaxis()->SetRangeUser(0.,40.);
1231   hCON->GetYaxis()->SetTitleOffset(1.2);
1232   hCON->Draw("bar1"); leg->AddEntry(hCON, "Total", "f");
1233   hCON->SetMaximum(55.);
1234
1235   hBAR->Scale(100./scale);
1236   hBAR->SetFillColor(kGreen);hBAR->SetLineColor(kGreen);
1237   hBAR->SetBarWidth(0.2);
1238   hBAR->SetBarOffset(0.2);
1239   hBAR->SetTitle("");
1240   hBAR->SetStats(kFALSE);
1241   hBAR->GetYaxis()->SetRangeUser(0.,40.);
1242   hBAR->GetYaxis()->SetTitleOffset(1.2);
1243   hBAR->Draw("bar1same"); leg->AddEntry(hBAR, "Barrel", "f");
1244
1245   hSTA->Scale(100./scale);
1246   hSTA->SetFillColor(kBlue);hSTA->SetLineColor(kBlue);
1247   hSTA->SetBarWidth(0.2);
1248   hSTA->SetBarOffset(0.4);
1249   hSTA->SetTitle("");
1250   hSTA->SetStats(kFALSE);
1251   hSTA->GetYaxis()->SetRangeUser(0.,40.);
1252   hSTA->GetYaxis()->SetTitleOffset(1.2);
1253   hSTA->Draw("bar1same"); leg->AddEntry(hSTA, "Stand Alone", "f");
1254   leg->Draw();
1255   gPad->Update();
1256   return hCON;
1257 }
1258
1259 //________________________________________________________
1260 void AliTRDcheckDET::MakePlotnTrackletsVsP(){
1261   //
1262   // Plot abundance of tracks with number of tracklets as function of momentum
1263   //
1264   TH1 *hLayer[6];
1265   TH2 *hBar = (TH2F *)fContainer->FindObject("hNtlsBAR");
1266   TLegend *leg = new TLegend(0.13, 0.75, 0.2, 0.9);
1267   leg->SetBorderSize(0);
1268   leg->SetHeader("Tracklet/Track");
1269   leg->SetFillStyle(0);
1270   Color_t color[6] = {kBlue, kOrange, kBlack, kGreen, kCyan, kRed};
1271   Bool_t first = kTRUE;
1272   for(Int_t itl = 1; itl <= 6; itl++){
1273     hLayer[itl-1]= hBar->ProjectionX(Form("ntl%d", itl), itl, itl);
1274     hLayer[itl-1]->Scale(100./hLayer[itl-1]->Integral());
1275     hLayer[itl-1]->SetLineColor(color[itl-1]);
1276     hLayer[itl-1]->GetYaxis()->SetRangeUser(0., 60.);
1277     hLayer[itl-1]->GetXaxis()->SetMoreLogLabels();
1278     hLayer[itl-1]->SetYTitle("Prob. [%]");
1279     if(first){
1280       hLayer[itl-1]->Draw("c");
1281       first=kFALSE;
1282     } else
1283       hLayer[itl-1]->Draw("csame");
1284     leg->AddEntry(hLayer[itl-1], Form("n = %d", itl),"l");
1285   }
1286   leg->Draw();
1287   gPad->Update();
1288 }
1289
1290 //________________________________________________________
1291 TH1* AliTRDcheckDET::MakePlotPulseHeight(){
1292   //
1293   // Create Plot of the Pluse Height Spectrum
1294   //
1295   TH1 *h, *h1, *h2;
1296   TObjArray *arr = (TObjArray*)fContainer->FindObject("<PH>");
1297   h = (TH1F*)arr->At(0);
1298   h->SetMarkerStyle(24);
1299   h->SetMarkerColor(kBlack);
1300   h->SetLineColor(kBlack);
1301   h->Draw("e1");
1302   // Trending for the pulse height: plateau value, slope and timebin of the maximum
1303   TLinearFitter fit(1,"pol1");
1304   Double_t time = 0.;
1305   for(Int_t itime = 10; itime <= 20; itime++){
1306     time = static_cast<Double_t>(itime);
1307     fit.AddPoint(&time, h->GetBinContent(itime + 1), h->GetBinError(itime + 1));
1308   }
1309   fit.Eval();
1310   Double_t plateau = fit.GetParameter(0) + 12 * fit.GetParameter(1);
1311   Double_t slope = fit.GetParameter(1);
1312   PutTrendValue("PHplateau", plateau);
1313   PutTrendValue("PHslope", slope);
1314   PutTrendValue("PHamplificationPeak", static_cast<Double_t>(h->GetMaximumBin()-1));
1315   AliDebug(1, Form("plateau %f, slope %f, MaxTime %f", plateau, slope, static_cast<Double_t>(h->GetMaximumBin()-1)));
1316 //   copy the second histogram in a new one with the same x-dimension as the phs with respect to time
1317   h1 = (TH1F *)arr->At(1);
1318   h2 = new TH1F("hphs1","Average PH", 31, -0.5, 30.5);
1319   for(Int_t ibin = h1->GetXaxis()->GetFirst(); ibin < h1->GetNbinsX(); ibin++) 
1320     h2->SetBinContent(ibin, h1->GetBinContent(ibin));
1321   h2->SetMarkerStyle(22);
1322   h2->SetMarkerColor(kBlue);
1323   h2->SetLineColor(kBlue);
1324   h2->Draw("e1same");
1325   gPad->Update();
1326 //   create axis according to the histogram dimensions of the original second histogram
1327   TGaxis *axis = new TGaxis(gPad->GetUxmin(),
1328                     gPad->GetUymax(),
1329                     gPad->GetUxmax(),
1330                     gPad->GetUymax(),
1331                     -0.08, 4.88, 510,"-L");
1332   axis->SetLineColor(kBlue);
1333   axis->SetLabelColor(kBlue);
1334   axis->SetTextColor(kBlue);
1335   axis->SetTitle("x_{0}-x_{c} [cm]");
1336   axis->Draw();
1337   return h1;
1338 }
1339
1340 //________________________________________________________
1341 Bool_t AliTRDcheckDET::MakeBarPlot(TH1 *histo, Int_t color){
1342   //
1343   // Draw nice bar plots
1344   //
1345   if(!histo->GetEntries()) return kFALSE;
1346   histo->Scale(100./histo->Integral());
1347   histo->SetFillColor(color);
1348   histo->SetBarOffset(.2);
1349   histo->SetBarWidth(.6);
1350   histo->Draw("bar1");
1351   return kTRUE;
1352 }