protect function against missing data
[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 <TGraphErrors.h>
36 #include <TLegend.h>
37 #include <TLinearFitter.h>
38 #include <TMath.h>
39 #include <TMap.h>
40 #include <TObjArray.h>
41 #include <TObject.h>
42 #include <TObjString.h>
43
44 #include <TPad.h>
45 #include <TProfile.h>
46 #include <TProfile2D.h>
47 #include <TROOT.h>
48 #include <TChain.h>
49
50 #include "AliLog.h"
51 #include "AliTRDcluster.h"
52 #include "AliESDHeader.h"
53 #include "AliESDRun.h"
54 #include "AliESDtrack.h"
55 #include "AliExternalTrackParam.h"
56 #include "AliTRDgeometry.h"
57 #include "AliTRDpadPlane.h"
58 #include "AliTRDSimParam.h"
59 #include "AliTRDseedV1.h"
60 #include "AliTRDtrackV1.h"
61 #include "AliTRDtrackerV1.h"
62 #include "AliTRDReconstructor.h"
63 #include "AliTrackReference.h"
64 #include "AliTrackPointArray.h"
65 #include "AliTracker.h"
66 #include "TTreeStream.h"
67
68 #include "info/AliTRDtrackInfo.h"
69 #include "info/AliTRDeventInfo.h"
70 #include "AliTRDcheckDET.h"
71
72 #include <cstdio>
73 #include <iostream>
74
75 ClassImp(AliTRDcheckDET)
76
77 //_______________________________________________________
78 AliTRDcheckDET::AliTRDcheckDET():
79   AliTRDrecoTask()
80   ,fEventInfo(NULL)
81   ,fTriggerNames(NULL)
82   ,fReconstructor(NULL)
83   ,fGeo(NULL)
84   ,fFlags(0)
85 {
86   //
87   // Default constructor
88   //
89   SetNameTitle("checkDET", "Basic TRD data checker");
90 }
91
92 //_______________________________________________________
93 AliTRDcheckDET::AliTRDcheckDET(char* name):
94   AliTRDrecoTask(name, "Basic TRD data checker")
95   ,fEventInfo(NULL)
96   ,fTriggerNames(NULL)
97   ,fReconstructor(NULL)
98   ,fGeo(NULL)
99   ,fFlags(0)
100 {
101   //
102   // Default constructor
103   //
104   DefineInput(2, AliTRDeventInfo::Class());
105
106   fReconstructor = new AliTRDReconstructor;
107   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
108   fGeo = new AliTRDgeometry;
109   InitFunctorList();
110 }
111
112
113 //_______________________________________________________
114 AliTRDcheckDET::~AliTRDcheckDET(){
115   //
116   // Destructor
117   // 
118   if(fTriggerNames) delete fTriggerNames;
119   delete fReconstructor;
120   delete fGeo;
121 }
122
123 //_______________________________________________________
124 void AliTRDcheckDET::ConnectInputData(Option_t *opt){
125   //
126   // Connect the Input data with the task
127   //
128   AliTRDrecoTask::ConnectInputData(opt);
129   fEventInfo = dynamic_cast<AliTRDeventInfo *>(GetInputData(2));
130 }
131
132 //_______________________________________________________
133 void AliTRDcheckDET::UserCreateOutputObjects(){
134   //
135   // Create Output Objects
136   //
137   OpenFile(1,"RECREATE");
138   fContainer = Histos();
139   if(!fTriggerNames) fTriggerNames = new TMap();
140 }
141
142 //_______________________________________________________
143 void AliTRDcheckDET::UserExec(Option_t *opt){
144   //
145   // Execution function
146   // Filling TRD quality histos
147   //
148   if(!HasMCdata() && fEventInfo->GetEventHeader()->GetEventType() != 7) return; // For real data we select only physical events
149   AliTRDrecoTask::UserExec(opt);  
150   Int_t nTracks = 0;            // Count the number of tracks per event
151   Int_t triggermask = fEventInfo->GetEventHeader()->GetTriggerMask();
152   TString triggername =  fEventInfo->GetRunInfo()->GetFiredTriggerClasses(triggermask);
153   AliDebug(6, Form("Trigger cluster: %d, Trigger class: %s\n", triggermask, triggername.Data()));
154   dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger))->Fill(triggermask);
155   for(Int_t iti = 0; iti < fTracks->GetEntriesFast(); iti++){
156     if(!fTracks->UncheckedAt(iti)) continue;
157     AliTRDtrackInfo *fTrackInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iti));
158     if(!fTrackInfo->GetTrack()) continue;
159     nTracks++;
160   }
161
162   if(nTracks){
163     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks))->Fill(triggermask);
164     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtracksEvent))->Fill(nTracks);
165   }
166   if(triggermask <= 20 && !fTriggerNames->FindObject(Form("%d", triggermask))){
167     fTriggerNames->Add(new TObjString(Form("%d", triggermask)), new TObjString(triggername));
168     // also set the label for both histograms
169     TH1 *histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks));
170     histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
171     histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger));
172     histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
173   }
174   PostData(1, fContainer);
175 }
176
177
178 //_______________________________________________________
179 Bool_t AliTRDcheckDET::PostProcess(){
180   //
181   // Do Postprocessing (for the moment set the number of Reference histograms)
182   //
183   
184   TH1 * h = NULL;
185   
186   // Calculate of the trigger clusters purity
187   h = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTrigger"));
188   TH1F *h1 = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTriggerTracks"));
189   h1->Divide(h);
190   Float_t purities[20], val = 0;
191   TString triggernames[20];
192   Int_t nTriggerClasses = 0;
193   for(Int_t ibin = 1; ibin <= h->GetNbinsX(); ibin++){
194     if((val = h1->GetBinContent(ibin))){
195       purities[nTriggerClasses] = val;
196       triggernames[nTriggerClasses] = h1->GetXaxis()->GetBinLabel(ibin);
197       nTriggerClasses++;
198     }
199   }
200   h = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kTriggerPurity));
201   TAxis *ax = h->GetXaxis();
202   for(Int_t itrg = 0; itrg < nTriggerClasses; itrg++){
203     h->Fill(itrg, purities[itrg]);
204     ax->SetBinLabel(itrg+1, triggernames[itrg].Data());
205   }
206   ax->SetRangeUser(-0.5, nTriggerClasses+.5);
207   h->GetYaxis()->SetRangeUser(0,1);
208
209   // track status
210   h=dynamic_cast<TH1F*>(fContainer->At(kTrackStatus));
211   Float_t ok = h->GetBinContent(1);
212   Int_t nerr = h->GetNbinsX();
213   for(Int_t ierr=nerr; ierr--;){
214     h->SetBinContent(ierr+1, 1.e2*h->GetBinContent(ierr+1)/ok);
215   }
216   h->SetBinContent(1, 0.);
217
218   // tracklet status
219   
220   TObjArray *arr = dynamic_cast<TObjArray*>(fContainer->UncheckedAt(kTrackletStatus));
221   for(Int_t ily = AliTRDgeometry::kNlayer; ily--;){
222     h=dynamic_cast<TH1F*>(arr->At(ily));
223     Float_t okB = h->GetBinContent(1);
224     Int_t nerrB = h->GetNbinsX();
225     for(Int_t ierr=nerrB; ierr--;){
226       h->SetBinContent(ierr+1, 1.e2*h->GetBinContent(ierr+1)/okB);
227     }
228     h->SetBinContent(1, 0.);
229   }
230
231   fNRefFigures = 17;
232
233   return kTRUE;
234 }
235
236 //_______________________________________________________
237 Bool_t AliTRDcheckDET::GetRefFigure(Int_t ifig){
238   //
239   // Setting Reference Figures
240   //
241   enum FigureType_t{
242     kFigNclustersTrack,
243     kFigNclustersTracklet,
244     kFigNtrackletsTrack,
245     kFigNTrackletsP,
246     kFigNtrackletsCross,
247     kFigNtrackletsFindable,
248     kFigNtracksEvent,
249     kFigNtracksSector,
250     kFigTrackStatus,
251     kFigTrackletStatus,
252     kFigChi2,
253     kFigPH,
254     kFigChargeCluster,
255     kFigChargeTracklet,
256     kFigNeventsTrigger,
257     kFigNeventsTriggerTracks,
258     kFigTriggerPurity
259   };
260   gPad->SetLogy(0);
261   gPad->SetLogx(0);
262   TH1 *h = NULL; TObjArray *arr=NULL;
263   TLegend *leg = NULL;
264   Bool_t kFIRST(1);
265   switch(ifig){
266   case kFigNclustersTrack:
267     (h=(TH1F*)fContainer->FindObject("hNcls"))->Draw("pl");
268     PutTrendValue("NClustersTrack", h->GetMean());
269     PutTrendValue("NClustersTrackRMS", h->GetRMS());
270     return kTRUE;
271   case kFigNclustersTracklet:
272     (h =(TH1F*)fContainer->FindObject("hNclTls"))->Draw("pc");
273     PutTrendValue("NClustersTracklet", h->GetMean());
274     PutTrendValue("NClustersTrackletRMS", h->GetRMS());
275     return kTRUE;
276   case kFigNtrackletsTrack:
277     h=MakePlotNTracklets();
278     PutTrendValue("NTrackletsTrack", h->GetMean());
279     PutTrendValue("NTrackletsTrackRMS", h->GetRMS());
280     return kTRUE;
281   case kFigNTrackletsP:
282     MakePlotnTrackletsVsP();
283     return kTRUE;
284   case kFigNtrackletsCross:
285     h = (TH1F*)fContainer->FindObject("hNtlsCross");
286     if(!MakeBarPlot(h, kRed)) break;
287     PutTrendValue("NTrackletsCross", h->GetMean());
288     PutTrendValue("NTrackletsCrossRMS", h->GetRMS());
289     return kTRUE;
290   case kFigNtrackletsFindable:
291     h = (TH1F*)fContainer->FindObject("hNtlsFindable");
292     if(!MakeBarPlot(h, kGreen)) break;
293     PutTrendValue("NTrackletsFindable", h->GetMean());
294     PutTrendValue("NTrackletsFindableRMS", h->GetRMS());
295     return kTRUE;
296   case kFigNtracksEvent:
297     (h = (TH1F*)fContainer->FindObject("hNtrks"))->Draw("pl");
298     PutTrendValue("NTracksEvent", h->GetMean());
299     PutTrendValue("NTracksEventRMS", h->GetRMS());
300     return kTRUE;
301   case kFigNtracksSector:
302     h = (TH1F*)fContainer->FindObject("hNtrksSector");
303     if(!MakeBarPlot(h, kGreen)) break;
304     PutTrendValue("NTracksSector", h->Integral()/h->GetNbinsX());
305     return kTRUE;
306   case kFigTrackStatus:
307     if(!(h=(TH1F *)fContainer->FindObject("hTrackStatus"))) break;
308     h->GetXaxis()->SetRangeUser(0.5, -1);
309     h->GetYaxis()->CenterTitle();
310     h->Draw("c");
311     PutTrendValue("TrackStatus", h->Integral());
312     gPad->SetLogy(0);
313     return kTRUE;
314   case kFigTrackletStatus:
315     if(!(arr = dynamic_cast<TObjArray*>(fContainer->At(kTrackletStatus)))) break;
316     leg = new TLegend(.68, .7, .97, .97);
317     leg->SetBorderSize(0);leg->SetFillStyle(0);
318     leg->SetHeader("TRD layer");
319     for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
320       if(!(h=dynamic_cast<TH1F*>(arr->At(ily)))) continue;
321       if(kFIRST){
322         h->Draw("pl");
323         h->GetXaxis()->SetRangeUser(0.5, -1);
324         h->GetYaxis()->CenterTitle();
325         kFIRST = kFALSE;
326       } else h->Draw("samepl");
327       leg->AddEntry(h, Form("ly = %d", ily), "l");
328       PutTrendValue(Form("TrackletStatus%d", ily), h->Integral());
329     }
330     leg->Draw();
331     gPad->SetLogy(0);
332     return kTRUE;
333   case kFigChi2:
334     MakePlotChi2();
335     return kTRUE;
336   case kFigPH:
337     gPad->SetMargin(0.125, 0.015, 0.1, 0.1);
338     MakePlotPulseHeight();
339     gPad->SetLogy(0);
340     return kTRUE;
341   case kFigChargeCluster:
342     (h = (TH1F*)fContainer->FindObject("hQcl"))->Draw("c");
343     gPad->SetLogy(1);
344     PutTrendValue("ChargeCluster", h->GetMaximumBin());
345     PutTrendValue("ChargeClusterRMS", h->GetRMS());
346     return kTRUE;
347   case kFigChargeTracklet:
348     (h=(TH1F*)fContainer->FindObject("hQtrklt"))->Draw("c");
349     PutTrendValue("ChargeTracklet", h->GetMaximumBin());
350     PutTrendValue("ChargeTrackletRMS", h->GetRMS());
351     return kTRUE;
352   case kFigNeventsTrigger:
353     ((TH1F*)fContainer->FindObject("hEventsTrigger"))->Draw("");
354     return kTRUE;
355   case kFigNeventsTriggerTracks:
356     ((TH1F*)fContainer->FindObject("hEventsTriggerTracks"))->Draw("");
357     return kTRUE;
358   case kFigTriggerPurity: 
359     if(!MakeBarPlot((TH1F*)fContainer->FindObject("hTriggerPurity"), kGreen)) break;
360     break;
361   default:
362     break;
363   }
364   AliInfo(Form("Reference plot [%d] missing result", ifig));
365   return kFALSE;
366 }
367
368 //_______________________________________________________
369 TObjArray *AliTRDcheckDET::Histos(){
370   //
371   // Create QA histograms
372   //
373     
374   if(fContainer) return fContainer;
375   
376   fContainer = new TObjArray(20);
377   //fContainer->SetOwner(kTRUE);
378
379   // Register Histograms
380   TH1 * h = NULL;
381   TAxis *ax = NULL;
382   if(!(h = (TH1F *)gROOT->FindObject("hNcls"))){
383     h = new TH1F("hNcls", "N_{clusters} / track", 181, -0.5, 180.5);
384     h->GetXaxis()->SetTitle("N_{clusters}");
385     h->GetYaxis()->SetTitle("Entries");
386   } else h->Reset();
387   fContainer->AddAt(h, kNclustersTrack);
388
389   if(!(h = (TH1F *)gROOT->FindObject("hNclTls"))){
390     h = new TH1F("hNclTls","N_{clusters} / tracklet", 51, -0.5, 50.5);
391     h->GetXaxis()->SetTitle("N_{clusters}");
392     h->GetYaxis()->SetTitle("Entries");
393   } else h->Reset();
394   fContainer->AddAt(h, kNclustersTracklet);
395
396   if(!(h = (TH1F *)gROOT->FindObject("hNtls"))){
397     h = new TH1F("hNtls", "N_{tracklets} / track", AliTRDgeometry::kNlayer, 0.5, 6.5);
398     h->GetXaxis()->SetTitle("N^{tracklet}");
399     h->GetYaxis()->SetTitle("freq. [%]");
400   } else h->Reset();
401   fContainer->AddAt(h, kNtrackletsTrack);
402
403   if(!(h = (TH1F *)gROOT->FindObject("htlsSTA"))){
404     h = new TH1F("hNtlsSTA", "N_{tracklets} / track (Stand Alone)", AliTRDgeometry::kNlayer, 0.5, 6.5);
405     h->GetXaxis()->SetTitle("N^{tracklet}");
406     h->GetYaxis()->SetTitle("freq. [%]");
407   }
408   fContainer->AddAt(h, kNtrackletsSTA);
409
410   // Binning for momentum dependent tracklet Plots
411   const Int_t kNp(30);
412   Float_t P=0.2;
413   Float_t binsP[kNp+1], binsTrklt[AliTRDgeometry::kNlayer+1];
414   for(Int_t i=0;i<kNp+1; i++,P=TMath::Exp(i*i*.001)-1.) binsP[i]=P;
415   for(Int_t il = 0; il <= AliTRDgeometry::kNlayer; il++) binsTrklt[il] = 0.5 + il;
416   if(!(h = (TH1F *)gROOT->FindObject("htlsBAR"))){
417     // Make tracklets for barrel tracks momentum dependent (if we do not exceed min and max values)
418     h = new TH2F("hNtlsBAR", 
419     "N_{tracklets} / track;p [GeV/c];N^{tracklet};freq. [%]", 
420     kNp, binsP, AliTRDgeometry::kNlayer, binsTrklt);
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       return NULL;
723     }
724     AliExternalTrackParam *par(fkTrack->GetTrackIn());
725     if(!par){
726       AliError("Input track params missing");
727       return NULL;
728     }
729     hBarrel->Fill(par->P(), nTracklets);
730   } else {
731     // Stand alone Track: momentum dependence not usefull
732     method = 0;
733     if(!(hSta = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsSTA)))) {
734       AliWarning("Method: StandAlone.  Histogram not processed!");
735       return NULL;
736     }
737     hSta->Fill(nTracklets);
738   }
739
740   if(DebugLevel() > 2){
741     AliTRDseedV1 *tracklet = NULL;
742     Int_t sector = -1;
743     for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
744       if(!(tracklet = fkTrack->GetTracklet(itl)) || !(tracklet->IsOK())) continue;
745       sector = fGeo->GetSector(tracklet->GetDetector());
746       break;
747     }
748     (*DebugStream()) << "NTrackletsTrack"
749       << "Sector="      << sector
750       << "NTracklets="  << nTracklets
751       << "Method="      << method
752       << "p="           << p
753       << "\n";
754   }
755   if(DebugLevel() > 3){
756     if(nTracklets == 1){
757       // If we have one Tracklet, check in which layer this happens
758       Int_t layer = -1;
759       AliTRDseedV1 *tracklet = NULL;
760       for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
761         if((tracklet = fkTrack->GetTracklet(il)) && tracklet->IsOK()){layer =  il; break;}
762       }
763       if(layer >= 0){
764         (*DebugStream()) << "NTrackletsLayer"
765           << "Layer=" << layer
766           << "p=" << p
767           << "\n";
768       }
769     }
770   }
771   return h;
772 }
773
774
775 //_______________________________________________________
776 TH1 *AliTRDcheckDET::PlotNTrackletsRowCross(const AliTRDtrackV1 *track){
777   //
778   // Plot the number of tracklets
779   //
780   if(track) fkTrack = track;
781   if(!fkTrack){
782     AliWarning("No Track defined.");
783     return NULL;
784   }
785   TH1 *h = NULL;
786   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsCross)))){
787     AliWarning("No Histogram defined.");
788     return NULL;
789   }
790
791   Int_t ncross = 0;
792   AliTRDseedV1 *tracklet = NULL;
793   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
794     if(!(tracklet = fkTrack->GetTracklet(il)) || !tracklet->IsOK()) continue;
795
796     if(tracklet->IsRowCross()) ncross++;
797   }
798   h->Fill(ncross);
799   return h;
800 }
801
802 //_______________________________________________________
803 TH1 *AliTRDcheckDET::PlotFindableTracklets(const AliTRDtrackV1 *track){
804   //
805   // Plots the ratio of number of tracklets vs.
806   // number of findable tracklets
807   //
808   // Findable tracklets are defined as track prolongation
809   // to layer i does not hit the dead area +- epsilon
810   //
811   // In order to check whether tracklet hist active area in Layer i, 
812   // the track is refitted and the fitted position + an uncertainty 
813   // range is compared to the chamber border (also with a different
814   // uncertainty)
815   //
816   // For the track fit two cases are distinguished:
817   // If the track is a stand alone track (defined by the status bit 
818   // encoding, then the track is fitted with the tilted Rieman model
819   // Otherwise the track is fitted with the Kalman fitter in two steps:
820   // Since the track parameters are give at the outer point, we first 
821   // fit in direction inwards. Afterwards we fit again in direction outwards
822   // to extrapolate the track to layers which are not reached by the track
823   // For the Kalman model, the radial track points have to be shifted by
824   // a distance epsilon in the direction that we want to fit
825   //
826   const Float_t epsilon = 0.01;   // dead area tolerance
827   const Float_t epsilonR = 1;    // shift in radial direction of the anode wire position (Kalman filter only)
828   const Float_t deltaY = 0.7;    // Tolerance in the track position in y-direction
829   const Float_t deltaZ = 7.0;    // Tolerance in the track position in z-direction (Padlength)
830   Double_t xAnode[AliTRDgeometry::kNlayer] = {300.2, 312.8, 325.4, 338.0, 350.6, 363.2}; // Take the default X0
831  
832   if(track) fkTrack = track;
833   if(!fkTrack){
834     AliWarning("No Track defined.");
835     return NULL;
836   }
837   TH1 *h = NULL;
838   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsFindable)))){
839     AliWarning("No Histogram defined.");
840     return NULL;
841   }
842   Int_t nFound = 0, nFindable = 0;
843   Int_t stack = -1;
844   Double_t ymin = 0., ymax = 0., zmin = 0., zmax = 0.;
845   Double_t y = 0., z = 0.;
846   AliTRDseedV1 *tracklet = NULL;
847   AliTRDpadPlane *pp;  
848   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
849     if((tracklet = fkTrack->GetTracklet(il)) && tracklet->IsOK()){
850       tracklet->SetReconstructor(fReconstructor);
851       nFound++;
852     }
853   }
854   // 2 Different cases:
855   // 1st stand alone: here we cannot propagate, but be can do a Tilted Rieman Fit
856   // 2nd barrel track: here we propagate the track to the layers
857   AliTrackPoint points[6];
858   Float_t xyz[3];
859   memset(xyz, 0, sizeof(Float_t) * 3);
860   if(((fkESD->GetStatus() & AliESDtrack::kTRDout) > 0) && !((fkESD->GetStatus() & AliESDtrack::kTRDin) > 0)){
861     // stand alone track
862     for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
863       xyz[0] = xAnode[il];
864       points[il].SetXYZ(xyz);
865     }
866     AliTRDtrackerV1::FitRiemanTilt(const_cast<AliTRDtrackV1 *>(fkTrack), NULL, kTRUE, 6, points);
867   } else {
868     // barrel track
869     //
870     // 2 Steps:
871     // -> Kalman inwards
872     // -> Kalman outwards
873     AliTRDtrackV1 copyTrack(*fkTrack);  // Do Kalman on a (non-constant) copy of the track
874     AliTrackPoint pointsInward[6], pointsOutward[6];
875     for(Int_t il = AliTRDgeometry::kNlayer; il--;){
876       // In order to avoid complications in the Kalman filter if the track points have the same radial
877       // position like the tracklets, we have to shift the radial postion of the anode wire by epsilon
878       // in the direction we want to go
879       // The track points have to be in reverse order for the Kalman Filter inwards
880       xyz[0] = xAnode[AliTRDgeometry::kNlayer - il - 1] - epsilonR;
881       pointsInward[il].SetXYZ(xyz);
882       xyz[0] = xAnode[il] + epsilonR;
883       pointsOutward[il].SetXYZ(xyz);
884     }
885     /*for(Int_t ipt = 0; ipt < AliTRDgeometry::kNlayer; ipt++)
886       printf("%d. X = %f\n", ipt, points[ipt].GetX());*/
887     // Kalman inwards
888     AliTRDtrackerV1::FitKalman(&copyTrack, NULL, kFALSE, 6, pointsInward);
889     memcpy(points, pointsInward, sizeof(AliTrackPoint) * 6); // Preliminary store the inward results in the Array points
890     // Kalman outwards
891     AliTRDtrackerV1::FitKalman(&copyTrack, NULL, kTRUE, 6, pointsInward);
892     memcpy(points, pointsOutward, sizeof(AliTrackPoint) * AliTRDgeometry::kNlayer);
893   }
894   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
895     y = points[il].GetY();
896     z = points[il].GetZ();
897     if((stack = fGeo->GetStack(z, il)) < 0) continue; // Not findable
898     pp = fGeo->GetPadPlane(il, stack);
899     ymin = pp->GetCol0() + epsilon;
900     ymax = pp->GetColEnd() - epsilon; 
901     zmin = pp->GetRowEnd() + epsilon; 
902     zmax = pp->GetRow0() - epsilon;
903     // ignore y-crossing (material)
904     if((z + deltaZ > zmin && z - deltaZ < zmax) && (y + deltaY > ymin && y - deltaY < ymax)) nFindable++;
905       if(DebugLevel() > 3){
906         Double_t posTracklet[2] = {tracklet ? tracklet->GetYfit(0) : 0, tracklet ? tracklet->GetZfit(0) : 0};
907         Int_t hasTracklet = tracklet ? 1 : 0;
908         (*DebugStream())   << "FindableTracklets"
909           << "layer="     << il
910           << "ytracklet=" << posTracklet[0]
911           << "ytrack="    << y
912           << "ztracklet=" << posTracklet[1]
913           << "ztrack="    << z
914           << "tracklet="  << hasTracklet
915           << "\n";
916       }
917   }
918   
919   h->Fill(nFindable > 0 ? TMath::Min(nFound/static_cast<Double_t>(nFindable), 1.) : 1);
920   AliDebug(2, Form("Findable[Found]: %d[%d|%f]", nFindable, nFound, nFound/static_cast<Float_t>(nFindable > 0 ? nFindable : 1)));
921   return h;
922 }
923
924
925 //_______________________________________________________
926 TH1 *AliTRDcheckDET::PlotChi2(const AliTRDtrackV1 *track){
927   //
928   // Plot the chi2 of the track
929   //
930   if(track) fkTrack = track;
931   if(!fkTrack){
932     AliWarning("No Track defined.");
933     return NULL;
934   }
935   TH1 *h = NULL;
936   if(!(h = dynamic_cast<TH2S*>(fContainer->At(kChi2)))) {
937     AliWarning("No Histogram defined.");
938     return NULL;
939   }
940   Int_t n = fkTrack->GetNumberOfTracklets();
941   if(!n) return NULL;
942
943   h->Fill(n, fkTrack->GetChi2()/n);
944   return h;
945 }
946
947
948 //_______________________________________________________
949 TH1 *AliTRDcheckDET::PlotPHt(const AliTRDtrackV1 *track){
950   //
951   // Plot the average pulse height
952   //
953   if(track) fkTrack = track;
954   if(!fkTrack){
955     AliWarning("No Track defined.");
956     return NULL;
957   }
958   TProfile *h = NULL;
959   if(!(h = dynamic_cast<TProfile *>(((TObjArray*)(fContainer->At(kPH)))->At(0)))){
960     AliWarning("No Histogram defined.");
961     return NULL;
962   }
963   AliTRDseedV1 *tracklet = NULL;
964   AliTRDcluster *c = NULL;
965   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
966     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
967     Int_t crossing = Int_t(tracklet->IsRowCross());
968     Int_t detector = tracklet->GetDetector();
969     tracklet->ResetClusterIter();
970     while((c = tracklet->NextCluster())){
971       if(!IsUsingClustersOutsideChamber() && !c->IsInChamber()) continue;
972       Int_t localtime        = c->GetLocalTimeBin();
973       Double_t absoluteCharge = TMath::Abs(c->GetQ());
974       h->Fill(localtime, absoluteCharge);
975       if(DebugLevel() > 3){
976         Double_t distance[2];
977         GetDistanceToTracklet(distance, tracklet, c);
978         Float_t theta = TMath::ATan(tracklet->GetZref(1));
979         Float_t phi = TMath::ATan(tracklet->GetYref(1));
980         Float_t momentum = 0.;
981         Int_t pdg = 0;
982         Int_t kinkIndex = fkESD ? fkESD->GetKinkIndex() : 0;
983         UShort_t TPCncls = fkESD ? fkESD->GetTPCncls() : 0;
984         if(fkMC){
985           if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
986           pdg = fkMC->GetPDG();
987         }
988         (*DebugStream()) << "PHt"
989           << "Detector="        << detector
990           << "crossing="        << crossing
991           << "Timebin="         << localtime
992           << "Charge="          << absoluteCharge
993           << "momentum="        << momentum
994           << "pdg="                             << pdg
995           << "theta="                   << theta
996           << "phi="                             << phi
997           << "kinkIndex="       << kinkIndex
998           << "TPCncls="         << TPCncls
999           << "dy="        << distance[0]
1000           << "dz="        << distance[1]
1001           << "c.="        << c
1002           << "\n";
1003       }
1004     }
1005   }
1006   return h;
1007 }
1008
1009 //_______________________________________________________
1010 TH1 *AliTRDcheckDET::PlotPHx(const AliTRDtrackV1 *track){
1011   //
1012   // Plots the average pulse height vs the distance from the anode wire
1013   // (plus const anode wire offset)
1014   //
1015   if(track) fkTrack = track;
1016   if(!fkTrack){
1017     AliWarning("No Track defined.");
1018     return NULL;
1019   }
1020   TProfile *h = NULL;
1021   if(!(h = dynamic_cast<TProfile *>(((TObjArray*)(fContainer->At(kPH)))->At(1)))){
1022     AliWarning("No Histogram defined.");
1023     return NULL;
1024   }
1025   Float_t offset = .5*AliTRDgeometry::CamHght();
1026   AliTRDseedV1 *tracklet = NULL;
1027   AliTRDcluster *c = NULL;
1028   Double_t distance = 0;
1029   Double_t x, y;
1030   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
1031     if(!(tracklet = fkTrack->GetTracklet(itl)) || !(tracklet->IsOK())) continue;
1032     tracklet->ResetClusterIter();
1033     while((c = tracklet->NextCluster())){
1034       if(!IsUsingClustersOutsideChamber() && !c->IsInChamber()) continue;
1035       x = c->GetX()-AliTRDcluster::GetXcorr(c->GetLocalTimeBin());
1036       y = c->GetY()-AliTRDcluster::GetYcorr(AliTRDgeometry::GetLayer(c->GetDetector()), c->GetCenter());
1037
1038       distance = tracklet->GetX0() - (c->GetX() + 0.3) + offset;
1039       h->Fill(distance, TMath::Abs(c->GetQ()));
1040     }
1041   }  
1042   return h;
1043 }
1044
1045 //_______________________________________________________
1046 TH1 *AliTRDcheckDET::PlotChargeCluster(const AliTRDtrackV1 *track){
1047   //
1048   // Plot the cluster charge
1049   //
1050   if(track) fkTrack = track;
1051   if(!fkTrack){
1052     AliWarning("No Track defined.");
1053     return NULL;
1054   }
1055   TH1 *h = NULL;
1056   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeCluster)))){
1057     AliWarning("No Histogram defined.");
1058     return NULL;
1059   }
1060   AliTRDseedV1 *tracklet = NULL;
1061   AliTRDcluster *c = NULL;
1062   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
1063     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
1064     for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
1065       if(!(c = tracklet->GetClusters(itime))) continue;
1066       h->Fill(c->GetQ());
1067     }
1068   }
1069   return h;
1070 }
1071
1072 //_______________________________________________________
1073 TH1 *AliTRDcheckDET::PlotChargeTracklet(const AliTRDtrackV1 *track){
1074   //
1075   // Plot the charge deposit per chamber
1076   //
1077   if(track) fkTrack = track;
1078   if(!fkTrack){
1079     AliWarning("No Track defined.");
1080     return NULL;
1081   }
1082   TH1 *h = NULL;
1083   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeTracklet)))){
1084     AliWarning("No Histogram defined.");
1085     return NULL;
1086   }
1087   AliTRDseedV1 *tracklet = NULL;
1088   AliTRDcluster *c = NULL;
1089   Double_t qTot = 0;
1090   Int_t nTracklets =fkTrack->GetNumberOfTracklets();
1091   for(Int_t itl(0); itl < AliTRDgeometry::kNlayer; itl++){
1092     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
1093     qTot = 0.;
1094     for(Int_t ic = AliTRDseedV1::kNclusters; ic--;){
1095       if(!(c = tracklet->GetClusters(ic))) continue;
1096       qTot += TMath::Abs(c->GetQ());
1097     }
1098     h->Fill(qTot);
1099     if(DebugLevel() > 3){
1100       Int_t crossing = (Int_t)tracklet->IsRowCross();
1101       Int_t detector = tracklet->GetDetector();
1102       Float_t theta = TMath::ATan(tracklet->GetZfit(1));
1103       Float_t phi = TMath::ATan(tracklet->GetYfit(1));
1104       Float_t momentum = 0.;
1105       Int_t pdg = 0;
1106       Int_t kinkIndex = fkESD ? fkESD->GetKinkIndex() : 0;
1107       UShort_t nclsTPC = fkESD ? fkESD->GetTPCncls() : 0;
1108       if(fkMC){
1109               if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
1110         pdg = fkMC->GetPDG();
1111       }
1112       (*DebugStream()) << "ChargeTracklet"
1113         << "Detector="  << detector
1114         << "crossing="  << crossing
1115         << "momentum="  << momentum
1116         << "nTracklets="<< nTracklets
1117         << "pdg="                               << pdg
1118         << "theta="                     << theta
1119         << "phi="                               << phi
1120         << "kinkIndex=" << kinkIndex
1121         << "TPCncls="           << nclsTPC
1122         << "QT="        << qTot
1123         << "\n";
1124     }
1125   }
1126   return h;
1127 }
1128
1129 //_______________________________________________________
1130 TH1 *AliTRDcheckDET::PlotNTracksSector(const AliTRDtrackV1 *track){
1131   //
1132   // Plot the number of tracks per Sector
1133   //
1134   if(track) fkTrack = track;
1135   if(!fkTrack){
1136     AliWarning("No Track defined.");
1137     return NULL;
1138   }
1139   TH1 *h = NULL;
1140   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtracksSector)))){
1141     AliWarning("No Histogram defined.");
1142     return NULL;
1143   }
1144
1145   // TODO we should compare with
1146   // sector = Int_t(track->GetAlpha() / AliTRDgeometry::GetAlpha());
1147
1148   AliTRDseedV1 *tracklet = NULL;
1149   Int_t sector = -1;
1150   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
1151     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
1152     sector = static_cast<Int_t>(tracklet->GetDetector()/AliTRDgeometry::kNdets);
1153     break;
1154   }
1155   h->Fill(sector);
1156   return h;
1157 }
1158
1159
1160 //________________________________________________________
1161 void AliTRDcheckDET::SetRecoParam(AliTRDrecoParam *r)
1162 {
1163
1164   fReconstructor->SetRecoParam(r);
1165 }
1166
1167 //________________________________________________________
1168 void AliTRDcheckDET::GetDistanceToTracklet(Double_t *dist, AliTRDseedV1 * const tracklet, AliTRDcluster * const c)
1169 {
1170   Float_t x = c->GetX();
1171   dist[0] = c->GetY() - tracklet->GetYat(x);
1172   dist[1] = c->GetZ() - tracklet->GetZat(x);
1173 }
1174
1175
1176 //_______________________________________________________
1177 TH1* AliTRDcheckDET::MakePlotChi2()
1178 {
1179 // Plot chi2/track normalized to number of degree of freedom 
1180 // (tracklets) and compare with the theoretical distribution.
1181 // 
1182 // Alex Bercuci <A.Bercuci@gsi.de>
1183
1184   TH2S *h2 = (TH2S*)fContainer->At(kChi2);
1185   TF1 f("fChi2", "[0]*pow(x, [1]-1)*exp(-0.5*x)", 0., 50.);
1186   TLegend *leg = new TLegend(.7,.7,.95,.95);
1187   leg->SetBorderSize(1); leg->SetHeader("Tracklets per Track");
1188   TH1D *h1 = NULL;
1189   Bool_t kFIRST = kTRUE;
1190   for(Int_t il=1; il<=h2->GetNbinsX(); il++){
1191     h1 = h2->ProjectionY(Form("pyChi2%d", il), il, il);
1192     if(h1->Integral()<50) continue;
1193     h1->Scale(1./h1->Integral());
1194     h1->SetMarkerStyle(7);h1->SetMarkerColor(il);
1195     h1->SetLineColor(il);h1->SetLineStyle(2);
1196     f.SetParameter(1, .5*il);f.SetLineColor(il);
1197     h1->Fit(&f, "QW+", kFIRST ? "pc": "pcsame");
1198     leg->AddEntry(h1, Form("%d", il), "l");
1199     if(kFIRST){
1200       h1->GetXaxis()->SetRangeUser(0., 25.);
1201     }
1202     kFIRST = kFALSE;
1203   }
1204   leg->Draw();
1205   gPad->SetLogy();
1206   return h1;
1207 }
1208
1209
1210 //________________________________________________________
1211 TH1* AliTRDcheckDET::MakePlotNTracklets(){
1212   //
1213   // Make nice bar plot of the number of tracklets in each method
1214   //
1215   TH2F *tmp = (TH2F *)fContainer->FindObject("hNtlsBAR");
1216   TH1D *hBAR = tmp->ProjectionY();
1217   TH1F *hSTA = (TH1F *)fContainer->FindObject("hNtlsSTA");
1218   TH1F *hCON = (TH1F *)fContainer->FindObject("hNtls");
1219   TLegend *leg = new TLegend(0.13, 0.75, 0.39, 0.89);
1220   leg->SetBorderSize(1);
1221   leg->SetFillColor(0);
1222
1223   Float_t scale = hCON->Integral();
1224   hCON->Scale(100./scale);
1225   hCON->SetFillColor(kRed);hCON->SetLineColor(kRed);
1226   hCON->SetBarWidth(0.2);
1227   hCON->SetBarOffset(0.6);
1228   hCON->SetStats(kFALSE);
1229   hCON->GetYaxis()->SetRangeUser(0.,40.);
1230   hCON->GetYaxis()->SetTitleOffset(1.2);
1231   hCON->Draw("bar1"); leg->AddEntry(hCON, "Total", "f");
1232   hCON->SetMaximum(55.);
1233
1234   hBAR->Scale(100./scale);
1235   hBAR->SetFillColor(kGreen);hBAR->SetLineColor(kGreen);
1236   hBAR->SetBarWidth(0.2);
1237   hBAR->SetBarOffset(0.2);
1238   hBAR->SetTitle("");
1239   hBAR->SetStats(kFALSE);
1240   hBAR->GetYaxis()->SetRangeUser(0.,40.);
1241   hBAR->GetYaxis()->SetTitleOffset(1.2);
1242   hBAR->Draw("bar1same"); leg->AddEntry(hBAR, "Barrel", "f");
1243
1244   hSTA->Scale(100./scale);
1245   hSTA->SetFillColor(kBlue);hSTA->SetLineColor(kBlue);
1246   hSTA->SetBarWidth(0.2);
1247   hSTA->SetBarOffset(0.4);
1248   hSTA->SetTitle("");
1249   hSTA->SetStats(kFALSE);
1250   hSTA->GetYaxis()->SetRangeUser(0.,40.);
1251   hSTA->GetYaxis()->SetTitleOffset(1.2);
1252   hSTA->Draw("bar1same"); leg->AddEntry(hSTA, "Stand Alone", "f");
1253   leg->Draw();
1254   gPad->Update();
1255   return hCON;
1256 }
1257
1258 //________________________________________________________
1259 void AliTRDcheckDET::MakePlotnTrackletsVsP(){
1260   //
1261   // Plot abundance of tracks with number of tracklets as function of momentum
1262   //
1263
1264
1265
1266
1267   Color_t color[AliTRDgeometry::kNlayer] = {kBlue, kOrange, kBlack, kGreen, kCyan, kRed};
1268   TH1 *h(NULL); TGraphErrors *g[AliTRDgeometry::kNlayer];
1269   for(Int_t itl(0); itl<AliTRDgeometry::kNlayer; itl++){
1270     g[itl] = new TGraphErrors();
1271     g[itl]->SetLineColor(color[itl]);
1272     g[itl]->SetMarkerColor(color[itl]);
1273     g[itl]->SetMarkerStyle(20 + itl);
1274   }
1275
1276   TH2 *hBar = (TH2F *)fContainer->FindObject("hNtlsBAR");
1277   TAxis *ax(hBar->GetXaxis());
1278   Int_t np(ax->GetNbins());
1279   for(Int_t ipBin(1); ipBin<np; ipBin++){
1280     h = hBar->ProjectionY("npBin", ipBin, ipBin);
1281     if(!Int_t(h->Integral())) return;
1282     h->Scale(100./h->Integral());
1283     Float_t p(ax->GetBinCenter(ipBin)); 
1284     Float_t dp(ax->GetBinWidth(ipBin)); 
1285     for(Int_t itl(AliTRDgeometry::kNlayer); itl--;){
1286       g[itl]->SetPoint(ipBin-1, p, h->GetBinContent(itl+1));
1287       g[itl]->SetPointError(ipBin-1, dp/3.46, h->GetBinError(itl+1));
1288     }
1289   }
1290
1291   TLegend *leg = new TLegend(0.5, 0.75, 0.9, 0.9);
1292   leg->SetBorderSize(0);
1293   leg->SetHeader("Tracklet/Track");
1294   leg->SetFillStyle(0);
1295   h = hBar->ProjectionX("npxBin"); h->Reset();
1296   h->GetYaxis()->SetRangeUser(0., 90.);
1297   h->GetXaxis()->SetMoreLogLabels();
1298   h->SetYTitle("Prob. [%]");
1299   h->Draw("p");
1300   for(Int_t itl(AliTRDgeometry::kNlayer); itl--;){
1301     g[itl]->Draw("pc");
1302     leg->AddEntry(g[itl], Form("n = %d", itl+1),"pl");
1303   }
1304
1305   leg->Draw();
1306   gPad->SetLogx();
1307 }
1308
1309 //________________________________________________________
1310 TH1* AliTRDcheckDET::MakePlotPulseHeight(){
1311   //
1312   // Create Plot of the Pluse Height Spectrum
1313   //
1314   TH1 *h, *h1, *h2;
1315   TObjArray *arr = (TObjArray*)fContainer->FindObject("<PH>");
1316   h = (TH1F*)arr->At(0);
1317   h->SetMarkerStyle(24);
1318   h->SetMarkerColor(kBlack);
1319   h->SetLineColor(kBlack);
1320   h->Draw("e1");
1321   // Trending for the pulse height: plateau value, slope and timebin of the maximum
1322   TLinearFitter fit(1,"pol1");
1323   Double_t time = 0.;
1324   for(Int_t itime = 10; itime <= 20; itime++){
1325     time = static_cast<Double_t>(itime);
1326     fit.AddPoint(&time, h->GetBinContent(itime + 1), h->GetBinError(itime + 1));
1327   }
1328   fit.Eval();
1329   Double_t plateau = fit.GetParameter(0) + 12 * fit.GetParameter(1);
1330   Double_t slope = fit.GetParameter(1);
1331   PutTrendValue("PHplateau", plateau);
1332   PutTrendValue("PHslope", slope);
1333   PutTrendValue("PHamplificationPeak", static_cast<Double_t>(h->GetMaximumBin()-1));
1334   AliDebug(1, Form("plateau %f, slope %f, MaxTime %f", plateau, slope, static_cast<Double_t>(h->GetMaximumBin()-1)));
1335 //   copy the second histogram in a new one with the same x-dimension as the phs with respect to time
1336   h1 = (TH1F *)arr->At(1);
1337   h2 = new TH1F("hphs1","Average PH", 31, -0.5, 30.5);
1338   for(Int_t ibin = h1->GetXaxis()->GetFirst(); ibin < h1->GetNbinsX(); ibin++) 
1339     h2->SetBinContent(ibin, h1->GetBinContent(ibin));
1340   h2->SetMarkerStyle(22);
1341   h2->SetMarkerColor(kBlue);
1342   h2->SetLineColor(kBlue);
1343   h2->Draw("e1same");
1344   gPad->Update();
1345 //   create axis according to the histogram dimensions of the original second histogram
1346   TGaxis *axis = new TGaxis(gPad->GetUxmin(),
1347                     gPad->GetUymax(),
1348                     gPad->GetUxmax(),
1349                     gPad->GetUymax(),
1350                     -0.08, 4.88, 510,"-L");
1351   axis->SetLineColor(kBlue);
1352   axis->SetLabelColor(kBlue);
1353   axis->SetTextColor(kBlue);
1354   axis->SetTitle("x_{0}-x_{c} [cm]");
1355   axis->Draw();
1356   return h1;
1357 }
1358
1359 //________________________________________________________
1360 Bool_t AliTRDcheckDET::MakeBarPlot(TH1 *histo, Int_t color){
1361   //
1362   // Draw nice bar plots
1363   //
1364   if(!histo->GetEntries()) return kFALSE;
1365   histo->Scale(100./histo->Integral());
1366   histo->SetFillColor(color);
1367   histo->SetBarOffset(.2);
1368   histo->SetBarWidth(.6);
1369   histo->Draw("bar1");
1370   return kTRUE;
1371 }