- fix bug https://savannah.cern.ch/bugs/?68562
[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::UserCreateOutputObjects(){
125   //
126   // Create Output Objects
127   //
128   if(!HasFunctorList()) InitFunctorList();
129   fContainer = Histos();
130   if(!fTriggerNames) fTriggerNames = new TMap();
131 }
132
133 //_______________________________________________________
134 void AliTRDcheckDET::UserExec(Option_t *opt){
135   //
136   // Execution function
137   // Filling TRD quality histos
138   //
139
140   fEventInfo = dynamic_cast<AliTRDeventInfo *>(GetInputData(2));
141   AliDebug(2, Form("EventInfo[%p] Header[%p]", (void*)fEventInfo, (void*)(fEventInfo?fEventInfo->GetEventHeader():NULL)));
142
143   AliTRDrecoTask::UserExec(opt);  
144   Int_t nTracks = 0;            // Count the number of tracks per event
145   Int_t triggermask = fEventInfo->GetEventHeader()->GetTriggerMask();
146   TString triggername =  fEventInfo->GetRunInfo()->GetFiredTriggerClasses(triggermask);
147   AliDebug(6, Form("Trigger cluster: %d, Trigger class: %s\n", triggermask, triggername.Data()));
148   dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger))->Fill(triggermask);
149   for(Int_t iti = 0; iti < fTracks->GetEntriesFast(); iti++){
150     if(!fTracks->UncheckedAt(iti)) continue;
151     AliTRDtrackInfo *fTrackInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iti));
152     if(!fTrackInfo->GetTrack()) continue;
153     nTracks++;
154   }
155
156   if(nTracks){
157     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks))->Fill(triggermask);
158     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtracksEvent))->Fill(nTracks);
159   }
160   if(triggermask <= 20 && !fTriggerNames->FindObject(Form("%d", triggermask))){
161     fTriggerNames->Add(new TObjString(Form("%d", triggermask)), new TObjString(triggername));
162     // also set the label for both histograms
163     TH1 *histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks));
164     histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
165     histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger));
166     histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
167   }
168   PostData(1, fContainer);
169 }
170
171
172 //_______________________________________________________
173 Bool_t AliTRDcheckDET::PostProcess(){
174   //
175   // Do Postprocessing (for the moment set the number of Reference histograms)
176   //
177   
178   TH1 * h = NULL;
179   
180   // Calculate of the trigger clusters purity
181   h = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTrigger"));
182   TH1F *h1 = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTriggerTracks"));
183   h1->Divide(h);
184   Float_t purities[20], val = 0;
185   TString triggernames[20];
186   Int_t nTriggerClasses = 0;
187   for(Int_t ibin = 1; ibin <= h->GetNbinsX(); ibin++){
188     if((val = h1->GetBinContent(ibin))){
189       purities[nTriggerClasses] = val;
190       triggernames[nTriggerClasses] = h1->GetXaxis()->GetBinLabel(ibin);
191       nTriggerClasses++;
192     }
193   }
194   h = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kTriggerPurity));
195   TAxis *ax = h->GetXaxis();
196   for(Int_t itrg = 0; itrg < nTriggerClasses; itrg++){
197     h->Fill(itrg, purities[itrg]);
198     ax->SetBinLabel(itrg+1, triggernames[itrg].Data());
199   }
200   ax->SetRangeUser(-0.5, nTriggerClasses+.5);
201   h->GetYaxis()->SetRangeUser(0,1);
202
203   // track status
204   h=dynamic_cast<TH1F*>(fContainer->At(kTrackStatus));
205   Float_t ok = h->GetBinContent(1);
206   Int_t nerr = h->GetNbinsX();
207   for(Int_t ierr=nerr; ierr--;){
208     h->SetBinContent(ierr+1, 1.e2*h->GetBinContent(ierr+1)/ok);
209   }
210   h->SetBinContent(1, 0.);
211
212   // tracklet status
213   
214   TObjArray *arr = dynamic_cast<TObjArray*>(fContainer->UncheckedAt(kTrackletStatus));
215   for(Int_t ily = AliTRDgeometry::kNlayer; ily--;){
216     h=dynamic_cast<TH1F*>(arr->At(ily));
217     Float_t okB = h->GetBinContent(1);
218     Int_t nerrB = h->GetNbinsX();
219     for(Int_t ierr=nerrB; ierr--;){
220       h->SetBinContent(ierr+1, 1.e2*h->GetBinContent(ierr+1)/okB);
221     }
222     h->SetBinContent(1, 0.);
223   }
224
225   fNRefFigures = 17;
226
227   return kTRUE;
228 }
229
230 //_______________________________________________________
231 Bool_t AliTRDcheckDET::GetRefFigure(Int_t ifig){
232   //
233   // Setting Reference Figures
234   //
235   enum FigureType_t{
236     kFigNclustersTrack,
237     kFigNclustersTracklet,
238     kFigNtrackletsTrack,
239     kFigNTrackletsP,
240     kFigNtrackletsCross,
241     kFigNtrackletsFindable,
242     kFigNtracksEvent,
243     kFigNtracksSector,
244     kFigTrackStatus,
245     kFigTrackletStatus,
246     kFigChi2,
247     kFigPH,
248     kFigChargeCluster,
249     kFigChargeTracklet,
250     kFigNeventsTrigger,
251     kFigNeventsTriggerTracks,
252     kFigTriggerPurity
253   };
254   gPad->SetLogy(0);
255   gPad->SetLogx(0);
256   TH1 *h = NULL; TObjArray *arr=NULL;
257   TLegend *leg = NULL;
258   Bool_t kFIRST(1);
259   switch(ifig){
260   case kFigNclustersTrack:
261     (h=(TH1F*)fContainer->FindObject("hNcls"))->Draw("pl");
262     PutTrendValue("NClustersTrack", h->GetMean());
263     PutTrendValue("NClustersTrackRMS", h->GetRMS());
264     return kTRUE;
265   case kFigNclustersTracklet:
266     (h =(TH1F*)fContainer->FindObject("hNclTls"))->Draw("pc");
267     PutTrendValue("NClustersTracklet", h->GetMean());
268     PutTrendValue("NClustersTrackletRMS", h->GetRMS());
269     return kTRUE;
270   case kFigNtrackletsTrack:
271     h=MakePlotNTracklets();
272     PutTrendValue("NTrackletsTrack", h->GetMean());
273     PutTrendValue("NTrackletsTrackRMS", h->GetRMS());
274     return kTRUE;
275   case kFigNTrackletsP:
276     MakePlotnTrackletsVsP();
277     return kTRUE;
278   case kFigNtrackletsCross:
279     h = (TH1F*)fContainer->FindObject("hNtlsCross");
280     if(!MakeBarPlot(h, kRed)) break;
281     PutTrendValue("NTrackletsCross", h->GetMean());
282     PutTrendValue("NTrackletsCrossRMS", h->GetRMS());
283     return kTRUE;
284   case kFigNtrackletsFindable:
285     h = (TH1F*)fContainer->FindObject("hNtlsFindable");
286     if(!MakeBarPlot(h, kGreen)) break;
287     PutTrendValue("NTrackletsFindable", h->GetMean());
288     PutTrendValue("NTrackletsFindableRMS", h->GetRMS());
289     return kTRUE;
290   case kFigNtracksEvent:
291     (h = (TH1F*)fContainer->FindObject("hNtrks"))->Draw("pl");
292     PutTrendValue("NTracksEvent", h->GetMean());
293     PutTrendValue("NTracksEventRMS", h->GetRMS());
294     return kTRUE;
295   case kFigNtracksSector:
296     h = (TH1F*)fContainer->FindObject("hNtrksSector");
297     if(!MakeBarPlot(h, kGreen)) break;
298     PutTrendValue("NTracksSector", h->Integral()/h->GetNbinsX());
299     return kTRUE;
300   case kFigTrackStatus:
301     if(!(h=(TH1F *)fContainer->FindObject("hTrackStatus"))) break;
302     h->GetXaxis()->SetRangeUser(0.5, -1);
303     h->GetYaxis()->CenterTitle();
304     h->Draw("c");
305     PutTrendValue("TrackStatus", h->Integral());
306     gPad->SetLogy(0);
307     return kTRUE;
308   case kFigTrackletStatus:
309     if(!(arr = dynamic_cast<TObjArray*>(fContainer->At(kTrackletStatus)))) break;
310     leg = new TLegend(.68, .7, .97, .97);
311     leg->SetBorderSize(0);leg->SetFillStyle(0);
312     leg->SetHeader("TRD layer");
313     for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
314       if(!(h=dynamic_cast<TH1F*>(arr->At(ily)))) continue;
315       if(kFIRST){
316         h->Draw("pl");
317         h->GetXaxis()->SetRangeUser(0.5, -1);
318         h->GetYaxis()->CenterTitle();
319         kFIRST = kFALSE;
320       } else h->Draw("samepl");
321       leg->AddEntry(h, Form("ly = %d", ily), "l");
322       PutTrendValue(Form("TrackletStatus%d", ily), h->Integral());
323     }
324     leg->Draw();
325     gPad->SetLogy(0);
326     return kTRUE;
327   case kFigChi2:
328     MakePlotChi2();
329     return kTRUE;
330   case kFigPH:
331     gPad->SetMargin(0.125, 0.015, 0.1, 0.1);
332     MakePlotPulseHeight();
333     gPad->SetLogy(0);
334     return kTRUE;
335   case kFigChargeCluster:
336     (h = (TH1F*)fContainer->FindObject("hQcl"))->Draw("c");
337     gPad->SetLogy(1);
338     PutTrendValue("ChargeCluster", h->GetMaximumBin());
339     PutTrendValue("ChargeClusterRMS", h->GetRMS());
340     return kTRUE;
341   case kFigChargeTracklet:
342     (h=(TH1F*)fContainer->FindObject("hQtrklt"))->Draw("c");
343     PutTrendValue("ChargeTracklet", h->GetMaximumBin());
344     PutTrendValue("ChargeTrackletRMS", h->GetRMS());
345     return kTRUE;
346   case kFigNeventsTrigger:
347     ((TH1F*)fContainer->FindObject("hEventsTrigger"))->Draw("");
348     return kTRUE;
349   case kFigNeventsTriggerTracks:
350     ((TH1F*)fContainer->FindObject("hEventsTriggerTracks"))->Draw("");
351     return kTRUE;
352   case kFigTriggerPurity: 
353     if(!MakeBarPlot((TH1F*)fContainer->FindObject("hTriggerPurity"), kGreen)) break;
354     break;
355   default:
356     break;
357   }
358   AliInfo(Form("Reference plot [%d] missing result", ifig));
359   return kFALSE;
360 }
361
362 //_______________________________________________________
363 TObjArray *AliTRDcheckDET::Histos(){
364   //
365   // Create QA histograms
366   //
367     
368   if(fContainer) return fContainer;
369   
370   fContainer = new TObjArray(20);
371   //fContainer->SetOwner(kTRUE);
372
373   // Register Histograms
374   TH1 * h = NULL;
375   TAxis *ax = NULL;
376   if(!(h = (TH1F *)gROOT->FindObject("hNcls"))){
377     h = new TH1F("hNcls", "N_{clusters} / track", 181, -0.5, 180.5);
378     h->GetXaxis()->SetTitle("N_{clusters}");
379     h->GetYaxis()->SetTitle("Entries");
380   } else h->Reset();
381   fContainer->AddAt(h, kNclustersTrack);
382
383   if(!(h = (TH1F *)gROOT->FindObject("hNclTls"))){
384     h = new TH1F("hNclTls","N_{clusters} / tracklet", 51, -0.5, 50.5);
385     h->GetXaxis()->SetTitle("N_{clusters}");
386     h->GetYaxis()->SetTitle("Entries");
387   } else h->Reset();
388   fContainer->AddAt(h, kNclustersTracklet);
389
390   if(!(h = (TH1F *)gROOT->FindObject("hNtls"))){
391     h = new TH1F("hNtls", "N_{tracklets} / track", AliTRDgeometry::kNlayer, 0.5, 6.5);
392     h->GetXaxis()->SetTitle("N^{tracklet}");
393     h->GetYaxis()->SetTitle("freq. [%]");
394   } else h->Reset();
395   fContainer->AddAt(h, kNtrackletsTrack);
396
397   if(!(h = (TH1F *)gROOT->FindObject("htlsSTA"))){
398     h = new TH1F("hNtlsSTA", "N_{tracklets} / track (Stand Alone)", AliTRDgeometry::kNlayer, 0.5, 6.5);
399     h->GetXaxis()->SetTitle("N^{tracklet}");
400     h->GetYaxis()->SetTitle("freq. [%]");
401   }
402   fContainer->AddAt(h, kNtrackletsSTA);
403
404   // Binning for momentum dependent tracklet Plots
405   const Int_t kNp(30);
406   Float_t P=0.2;
407   Float_t binsP[kNp+1], binsTrklt[AliTRDgeometry::kNlayer+1];
408   for(Int_t i=0;i<kNp+1; i++,P+=(TMath::Exp(i*i*.001)-1.)) binsP[i]=P;
409   for(Int_t il = 0; il <= AliTRDgeometry::kNlayer; il++) binsTrklt[il] = 0.5 + il;
410   if(!(h = (TH1F *)gROOT->FindObject("htlsBAR"))){
411     // Make tracklets for barrel tracks momentum dependent (if we do not exceed min and max values)
412     h = new TH2F("hNtlsBAR", 
413     "N_{tracklets} / track;p [GeV/c];N^{tracklet};freq. [%]", 
414     kNp, binsP, AliTRDgeometry::kNlayer, binsTrklt);
415   }
416   fContainer->AddAt(h, kNtrackletsBAR);
417
418   // 
419   if(!(h = (TH1F *)gROOT->FindObject("hNtlsCross"))){
420     h = new TH1F("hNtlsCross", "N_{tracklets}^{cross} / track", 7, -0.5, 6.5);
421     h->GetXaxis()->SetTitle("n_{row cross}");
422     h->GetYaxis()->SetTitle("freq. [%]");
423   } else h->Reset();
424   fContainer->AddAt(h, kNtrackletsCross);
425
426   if(!(h = (TH1F *)gROOT->FindObject("hNtlsFindable"))){
427     h = new TH1F("hNtlsFindable", "Found/Findable Tracklets" , 101, -0.005, 1.005);
428     h->GetXaxis()->SetTitle("r [a.u]");
429     h->GetYaxis()->SetTitle("Entries");
430   } else h->Reset();
431   fContainer->AddAt(h, kNtrackletsFindable);
432
433   if(!(h = (TH1F *)gROOT->FindObject("hNtrks"))){
434     h = new TH1F("hNtrks", "N_{tracks} / event", 100, 0, 100);
435     h->GetXaxis()->SetTitle("N_{tracks}");
436     h->GetYaxis()->SetTitle("Entries");
437   } else h->Reset();
438   fContainer->AddAt(h, kNtracksEvent);
439
440   if(!(h = (TH1F *)gROOT->FindObject("hNtrksSector"))){
441     h = new TH1F("hNtrksSector", "N_{tracks} / sector", AliTRDgeometry::kNsector, -0.5, 17.5);
442     h->GetXaxis()->SetTitle("sector");
443     h->GetYaxis()->SetTitle("freq. [%]");
444   } else h->Reset();
445   fContainer->AddAt(h, kNtracksSector);
446
447   if(!(h = (TH1F*)gROOT->FindObject("hTrackStatus"))){
448     const Int_t nerr = 7;
449     h = new TH1F("hTrackStatus", "Track Status", nerr, -0.5, nerr-0.5);
450     const Char_t *label[nerr] = {"OK", "PROL", "PROP", "AJST", "SNP", "TINI", "UPDT"};
451     ax = h->GetXaxis();
452     for(Int_t ierr=nerr; ierr--;) ax->SetBinLabel(ierr+1, label[ierr]);
453     h->SetYTitle("Relative Error to Good [%]");
454   }
455   fContainer->AddAt(h, kTrackStatus);
456
457   TObjArray *arr = new TObjArray(AliTRDgeometry::kNlayer);
458   arr->SetOwner(kTRUE);  arr->SetName("TrackletStatus");
459   fContainer->AddAt(arr, kTrackletStatus);
460   for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
461     if(!(h = (TH1F *)gROOT->FindObject(Form("hTrackletStatus%d", ily)))){
462       const Int_t nerr = 8;
463       h = new TH1F(Form("hTrackletStatus%d", ily), "Tracklet status", nerr, -0.5, nerr-0.5);
464       h->SetLineColor(ily+1);
465       const Char_t *label[nerr] = {"OK", "Geom", "Bound", "NoCl", "NoAttach", "NoClTr", "NoFit", "Chi2"};
466       ax = h->GetXaxis();
467       for(Int_t ierr=nerr; ierr--;) ax->SetBinLabel(ierr+1, label[ierr]);
468       h->SetYTitle("Relative Error to Good [%]");
469     } else h->Reset();
470     arr->AddAt(h, ily);
471   }
472
473   // <PH> histos
474   arr = new TObjArray(3);
475   arr->SetOwner(kTRUE);  arr->SetName("<PH>");
476   fContainer->AddAt(arr, kPH);
477   if(!(h = (TH1F *)gROOT->FindObject("hPHt"))){
478     h = new TProfile("hPHt", "<PH>", 31, -0.5, 30.5);
479     h->GetXaxis()->SetTitle("Time / 100ns");
480     h->GetYaxis()->SetTitle("<PH> [a.u]");
481   } else h->Reset();
482   arr->AddAt(h, 0);
483   if(!(h = (TH1F *)gROOT->FindObject("hPHx")))
484     h = new TProfile("hPHx", "<PH>", 31, -0.08, 4.88);
485   else h->Reset();
486   arr->AddAt(h, 1);
487   if(!(h = (TH2F *)gROOT->FindObject("hPH2D"))){
488     h = new TH2F("hPH2D", "Charge Distribution / time", 31, -0.5, 30.5, 100, 0, 1024);
489     h->GetXaxis()->SetTitle("Time / 100ns");
490     h->GetYaxis()->SetTitle("Charge / a.u.");
491   } else h->Reset();
492   arr->AddAt(h, 2);
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     AliDebug(4, "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     AliDebug(4, "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     AliDebug(4, "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     AliDebug(4, "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   if(!par) return NULL;
658   Double_t momentumRec = par->P();
659   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
660     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
661     Int_t n(tracklet->GetN());
662     nclusters += n;
663     if(DebugLevel() > 2){
664       Int_t crossing = Int_t(tracklet->IsRowCross());
665       Int_t detector = tracklet->GetDetector();
666       Float_t theta = TMath::ATan(tracklet->GetZref(1));
667       Float_t phi = TMath::ATan(tracklet->GetYref(1));
668       Float_t momentumMC = 0.;
669       Int_t pdg = 0;
670       Int_t kinkIndex = fkESD ? fkESD->GetKinkIndex() : 0;
671       UShort_t nclsTPC = fkESD ? fkESD->GetTPCncls() : 0;
672       if(fkMC){
673         if(fkMC->GetTrackRef()) momentumMC = fkMC->GetTrackRef()->P();
674         pdg = fkMC->GetPDG();
675       }
676       (*DebugStream()) << "NClustersTrack"
677         << "Detector="  << detector
678         << "crossing="  << crossing
679         << "momentumMC="<< momentumMC
680         << "momentumRec=" << momentumRec
681         << "pdg="                               << pdg
682         << "theta="                     << theta
683         << "phi="                               << phi
684         << "kinkIndex=" << kinkIndex
685         << "TPCncls="           << nclsTPC
686         << "TRDncls="   << n
687         << "\n";
688     }
689   }
690   h->Fill(nclusters);
691   return h;
692 }
693
694
695 //_______________________________________________________
696 TH1 *AliTRDcheckDET::PlotNTrackletsTrack(const AliTRDtrackV1 *track){
697   //
698   // Plot the number of tracklets
699   //
700   if(track) fkTrack = track;
701   if(!fkTrack){
702     AliDebug(4, "No Track defined.");
703     return NULL;
704   }
705   TH1 *h = NULL, *hSta = NULL; TH2 *hBarrel = NULL;
706   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsTrack)))){
707     AliWarning("No Histogram defined.");
708     return NULL;
709   }
710   Int_t nTracklets = fkTrack->GetNumberOfTracklets();
711   h->Fill(nTracklets);
712   if(!fkESD) return h;
713   Int_t status = fkESD->GetStatus();
714
715 /*  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);*/
716   Double_t p = 0.;
717   Int_t method = -1;    // to distinguish between stand alone and full barrel tracks in the debugging
718   if((status & AliESDtrack::kTRDin) != 0){
719     method = 1;
720     // Full Barrel Track: Save momentum dependence
721     if(!(hBarrel = dynamic_cast<TH2F *>(fContainer->At(kNtrackletsBAR)))){
722       AliWarning("Method: Barrel.  Histogram not processed!");
723       return NULL;
724     }
725     AliExternalTrackParam *par(fkTrack->GetTrackIn());
726     if(!par){
727       AliError("Input track params missing");
728       return NULL;
729     }
730     hBarrel->Fill(par->P(), nTracklets);
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       return NULL;
737     }
738     hSta->Fill(nTracklets);
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     AliDebug(4, "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     AliDebug(4, "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     AliDebug(4, "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     AliDebug(4, "No Track defined.");
957     return NULL;
958   }
959   TProfile *h = NULL; TH2F *phs2D = NULL;
960   if(!(h = dynamic_cast<TProfile *>(((TObjArray*)(fContainer->At(kPH)))->At(0)))){
961     AliWarning("No Histogram defined.");
962     return NULL;
963   }
964   if(!(phs2D = dynamic_cast<TH2F *>(((TObjArray*)(fContainer->At(kPH)))->At(2)))){
965     AliWarning("2D Pulse Height histogram not defined. Histogramm cannot be filled");
966   }
967   AliTRDseedV1 *tracklet = NULL;
968   AliTRDcluster *c = NULL;
969   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
970     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
971     Int_t crossing = Int_t(tracklet->IsRowCross());
972     Int_t detector = tracklet->GetDetector();
973     tracklet->ResetClusterIter();
974     while((c = tracklet->NextCluster())){
975       if(!IsUsingClustersOutsideChamber() && !c->IsInChamber()) continue;
976       Int_t localtime        = c->GetLocalTimeBin();
977       Double_t absoluteCharge = TMath::Abs(c->GetQ());
978       h->Fill(localtime, absoluteCharge);
979       phs2D->Fill(localtime, absoluteCharge); 
980       if(DebugLevel() > 3){
981         Double_t distance[2];
982         GetDistanceToTracklet(distance, tracklet, c);
983         Float_t theta = TMath::ATan(tracklet->GetZref(1));
984         Float_t phi = TMath::ATan(tracklet->GetYref(1));
985         Float_t momentum = 0.;
986         Int_t pdg = 0;
987         Int_t kinkIndex = fkESD ? fkESD->GetKinkIndex() : 0;
988         UShort_t TPCncls = fkESD ? fkESD->GetTPCncls() : 0;
989         if(fkMC){
990           if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
991           pdg = fkMC->GetPDG();
992         }
993         (*DebugStream()) << "PHt"
994           << "Detector="        << detector
995           << "crossing="        << crossing
996           << "Timebin="         << localtime
997           << "Charge="          << absoluteCharge
998           << "momentum="        << momentum
999           << "pdg="                             << pdg
1000           << "theta="                   << theta
1001           << "phi="                             << phi
1002           << "kinkIndex="       << kinkIndex
1003           << "TPCncls="         << TPCncls
1004           << "dy="        << distance[0]
1005           << "dz="        << distance[1]
1006           << "c.="        << c
1007           << "\n";
1008       }
1009     }
1010   }
1011   return h;
1012 }
1013
1014 //_______________________________________________________
1015 TH1 *AliTRDcheckDET::PlotPHx(const AliTRDtrackV1 *track){
1016   //
1017   // Plots the average pulse height vs the distance from the anode wire
1018   // (plus const anode wire offset)
1019   //
1020   if(track) fkTrack = track;
1021   if(!fkTrack){
1022      AliDebug(4, "No Track defined.");
1023      return NULL;
1024   }
1025   TProfile *h = NULL;
1026   if(!(h = dynamic_cast<TProfile *>(((TObjArray*)(fContainer->At(kPH)))->At(1)))){
1027     AliWarning("No Histogram defined.");
1028     return NULL;
1029   }
1030   AliTRDseedV1 *tracklet(NULL);
1031   AliTRDcluster *c(NULL);
1032   Double_t xd(0.), dqdl(0.);
1033   TVectorD vq(AliTRDseedV1::kNtb), vxd(AliTRDseedV1::kNtb), vdqdl(AliTRDseedV1::kNtb);
1034   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
1035     if(!(tracklet = fkTrack->GetTracklet(itl)) || !(tracklet->IsOK())) continue;
1036     Int_t det(tracklet->GetDetector());
1037     Bool_t rc(tracklet->IsRowCross());
1038     for(Int_t ic(0); ic<AliTRDseedV1::kNtb; ic++){
1039       Bool_t kFIRST(kFALSE);
1040       if(!(c = tracklet->GetClusters(ic))){
1041          if(!(c = tracklet->GetClusters(AliTRDseedV1::kNtb+ic))) continue;
1042       } else kFIRST=kTRUE;
1043       if(!IsUsingClustersOutsideChamber() && !c->IsInChamber()) continue;
1044       xd = tracklet->GetX0() - c->GetX(); vxd[ic] = xd;
1045       dqdl=tracklet->GetdQdl(ic); vdqdl[ic] = dqdl;
1046       vq[ic]=c->GetQ();
1047       if(kFIRST && (c = tracklet->GetClusters(AliTRDseedV1::kNtb+ic))) vq[ic]+=c->GetQ();
1048       h->Fill(xd, dqdl);
1049     }
1050     if(DebugLevel() > 3){
1051       (*DebugStream()) << "PHx"
1052         << "det="  << det
1053         << "rc="   << rc
1054         << "xd="   << &vxd
1055         << "q="    << &vq
1056         << "dqdl=" << &vdqdl
1057         << "\n";
1058     }
1059   }  
1060   return h;
1061 }
1062
1063 //_______________________________________________________
1064 TH1 *AliTRDcheckDET::PlotChargeCluster(const AliTRDtrackV1 *track){
1065   //
1066   // Plot the cluster charge
1067   //
1068   if(track) fkTrack = track;
1069   if(!fkTrack){
1070     AliDebug(4, "No Track defined.");
1071     return NULL;
1072   }
1073   TH1 *h = NULL;
1074   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeCluster)))){
1075     AliWarning("No Histogram defined.");
1076     return NULL;
1077   }
1078   AliTRDseedV1 *tracklet = NULL;
1079   AliTRDcluster *c = NULL;
1080   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
1081     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
1082     for(Int_t ic(0); ic < AliTRDseedV1::kNtb; ic++){
1083       Bool_t kFIRST(kFALSE);
1084       if(!(c = tracklet->GetClusters(ic))) {
1085         if(!(c = tracklet->GetClusters(AliTRDseedV1::kNtb+ic))) continue;
1086       } else kFIRST = kTRUE;
1087       Float_t q(c->GetQ());
1088       if(kFIRST && (c = tracklet->GetClusters(AliTRDseedV1::kNtb+ic))) q+=c->GetQ();
1089       h->Fill(q);
1090     }
1091   }
1092   return h;
1093 }
1094
1095 //_______________________________________________________
1096 TH1 *AliTRDcheckDET::PlotChargeTracklet(const AliTRDtrackV1 *track){
1097   //
1098   // Plot the charge deposit per chamber
1099   //
1100   if(track) fkTrack = track;
1101   if(!fkTrack){
1102     AliDebug(4, "No Track defined.");
1103     return NULL;
1104   }
1105   TH1 *h = NULL;
1106   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeTracklet)))){
1107     AliWarning("No Histogram defined.");
1108     return NULL;
1109   }
1110   AliTRDseedV1 *tracklet = NULL;
1111   AliTRDcluster *c = NULL;
1112   Double_t qTot = 0;
1113   Int_t nTracklets =fkTrack->GetNumberOfTracklets();
1114   for(Int_t itl(0); itl < AliTRDgeometry::kNlayer; itl++){
1115     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
1116     qTot = 0.;
1117     for(Int_t ic = AliTRDseedV1::kNclusters; ic--;){
1118       if(!(c = tracklet->GetClusters(ic))) continue;
1119       qTot += TMath::Abs(c->GetQ());
1120     }
1121     h->Fill(qTot);
1122     if(DebugLevel() > 3){
1123       Int_t crossing = (Int_t)tracklet->IsRowCross();
1124       Int_t detector = tracklet->GetDetector();
1125       Float_t theta = TMath::ATan(tracklet->GetZfit(1));
1126       Float_t phi = TMath::ATan(tracklet->GetYfit(1));
1127       Float_t momentum = 0.;
1128       Int_t pdg = 0;
1129       Int_t kinkIndex = fkESD ? fkESD->GetKinkIndex() : 0;
1130       UShort_t nclsTPC = fkESD ? fkESD->GetTPCncls() : 0;
1131       if(fkMC){
1132               if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
1133         pdg = fkMC->GetPDG();
1134       }
1135       (*DebugStream()) << "ChargeTracklet"
1136         << "Detector="  << detector
1137         << "crossing="  << crossing
1138         << "momentum="  << momentum
1139         << "nTracklets="<< nTracklets
1140         << "pdg="                               << pdg
1141         << "theta="                     << theta
1142         << "phi="                               << phi
1143         << "kinkIndex=" << kinkIndex
1144         << "TPCncls="           << nclsTPC
1145         << "QT="        << qTot
1146         << "\n";
1147     }
1148   }
1149   return h;
1150 }
1151
1152 //_______________________________________________________
1153 TH1 *AliTRDcheckDET::PlotNTracksSector(const AliTRDtrackV1 *track){
1154   //
1155   // Plot the number of tracks per Sector
1156   //
1157   if(track) fkTrack = track;
1158   if(!fkTrack){
1159     AliDebug(4, "No Track defined.");
1160     return NULL;
1161   }
1162   TH1 *h = NULL;
1163   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtracksSector)))){
1164     AliWarning("No Histogram defined.");
1165     return NULL;
1166   }
1167
1168   // TODO we should compare with
1169   // sector = Int_t(track->GetAlpha() / AliTRDgeometry::GetAlpha());
1170
1171   AliTRDseedV1 *tracklet = NULL;
1172   Int_t sector = -1;
1173   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
1174     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
1175     sector = static_cast<Int_t>(tracklet->GetDetector()/AliTRDgeometry::kNdets);
1176     break;
1177   }
1178   h->Fill(sector);
1179   return h;
1180 }
1181
1182
1183 //________________________________________________________
1184 void AliTRDcheckDET::SetRecoParam(AliTRDrecoParam *r)
1185 {
1186
1187   fReconstructor->SetRecoParam(r);
1188 }
1189
1190 //________________________________________________________
1191 void AliTRDcheckDET::GetDistanceToTracklet(Double_t *dist, AliTRDseedV1 * const tracklet, AliTRDcluster * const c)
1192 {
1193   Float_t x = c->GetX();
1194   dist[0] = c->GetY() - tracklet->GetYat(x);
1195   dist[1] = c->GetZ() - tracklet->GetZat(x);
1196 }
1197
1198
1199 //_______________________________________________________
1200 TH1* AliTRDcheckDET::MakePlotChi2()
1201 {
1202 // Plot chi2/track normalized to number of degree of freedom 
1203 // (tracklets) and compare with the theoretical distribution.
1204 // 
1205 // Alex Bercuci <A.Bercuci@gsi.de>
1206
1207   return NULL;
1208
1209   TH2S *h2 = (TH2S*)fContainer->At(kChi2);
1210   TF1 f("fChi2", "[0]*pow(x, [1]-1)*exp(-0.5*x)", 0., 50.);
1211   f.SetParLimits(1,1, 1e100);
1212   TLegend *leg = new TLegend(.7,.7,.95,.95);
1213   leg->SetBorderSize(1); leg->SetHeader("Tracklets per Track");
1214   TH1D *h1 = NULL;
1215   Bool_t kFIRST = kTRUE;
1216   for(Int_t il=1; il<=h2->GetNbinsX(); il++){
1217     h1 = h2->ProjectionY(Form("pyChi2%d", il), il, il);
1218     if(h1->Integral()<50) continue;
1219     h1->Scale(1./h1->Integral());
1220     h1->SetMarkerStyle(7);h1->SetMarkerColor(il);
1221     h1->SetLineColor(il);h1->SetLineStyle(2);
1222     f.SetParameter(1, .5*il);f.SetLineColor(il);
1223     h1->Fit(&f, "QW+", kFIRST ? "pc": "pcsame");
1224     leg->AddEntry(h1, Form("%d", il), "l");
1225     if(kFIRST){
1226       h1->GetXaxis()->SetRangeUser(0., 25.);
1227     }
1228     kFIRST = kFALSE;
1229   }
1230   leg->Draw();
1231   gPad->SetLogy();
1232   return h1;
1233 }
1234
1235
1236 //________________________________________________________
1237 TH1* AliTRDcheckDET::MakePlotNTracklets(){
1238   //
1239   // Make nice bar plot of the number of tracklets in each method
1240   //
1241   TH2F *tmp = (TH2F *)fContainer->FindObject("hNtlsBAR");
1242   TH1D *hBAR = tmp->ProjectionY();
1243   TH1F *hSTA = (TH1F *)fContainer->FindObject("hNtlsSTA");
1244   TH1F *hCON = (TH1F *)fContainer->FindObject("hNtls");
1245   TLegend *leg = new TLegend(0.13, 0.75, 0.39, 0.89);
1246   leg->SetBorderSize(1);
1247   leg->SetFillColor(0);
1248
1249   Float_t scale = hCON->Integral();
1250   hCON->Scale(100./scale);
1251   hCON->SetFillColor(kRed);hCON->SetLineColor(kRed);
1252   hCON->SetBarWidth(0.2);
1253   hCON->SetBarOffset(0.6);
1254   hCON->SetStats(kFALSE);
1255   hCON->GetYaxis()->SetRangeUser(0.,40.);
1256   hCON->GetYaxis()->SetTitleOffset(1.2);
1257   hCON->Draw("bar1"); leg->AddEntry(hCON, "Total", "f");
1258   hCON->SetMaximum(55.);
1259
1260   hBAR->Scale(100./scale);
1261   hBAR->SetFillColor(kGreen);hBAR->SetLineColor(kGreen);
1262   hBAR->SetBarWidth(0.2);
1263   hBAR->SetBarOffset(0.2);
1264   hBAR->SetTitle("");
1265   hBAR->SetStats(kFALSE);
1266   hBAR->GetYaxis()->SetRangeUser(0.,40.);
1267   hBAR->GetYaxis()->SetTitleOffset(1.2);
1268   hBAR->Draw("bar1same"); leg->AddEntry(hBAR, "Barrel", "f");
1269
1270   hSTA->Scale(100./scale);
1271   hSTA->SetFillColor(kBlue);hSTA->SetLineColor(kBlue);
1272   hSTA->SetBarWidth(0.2);
1273   hSTA->SetBarOffset(0.4);
1274   hSTA->SetTitle("");
1275   hSTA->SetStats(kFALSE);
1276   hSTA->GetYaxis()->SetRangeUser(0.,40.);
1277   hSTA->GetYaxis()->SetTitleOffset(1.2);
1278   hSTA->Draw("bar1same"); leg->AddEntry(hSTA, "Stand Alone", "f");
1279   leg->Draw();
1280   gPad->Update();
1281   return hCON;
1282 }
1283
1284 //________________________________________________________
1285 void AliTRDcheckDET::MakePlotnTrackletsVsP(){
1286   //
1287   // Plot abundance of tracks with number of tracklets as function of momentum
1288   //
1289
1290
1291
1292
1293   Color_t color[AliTRDgeometry::kNlayer] = {kBlue, kOrange, kBlack, kGreen, kCyan, kRed};
1294   TH1 *h(NULL); TGraphErrors *g[AliTRDgeometry::kNlayer];
1295   for(Int_t itl(0); itl<AliTRDgeometry::kNlayer; itl++){
1296     g[itl] = new TGraphErrors();
1297     g[itl]->SetLineColor(color[itl]);
1298     g[itl]->SetMarkerColor(color[itl]);
1299     g[itl]->SetMarkerStyle(20 + itl);
1300   }
1301
1302   TH2 *hBar = (TH2F *)fContainer->FindObject("hNtlsBAR");
1303   TAxis *ax(hBar->GetXaxis());
1304   Int_t np(ax->GetNbins());
1305   for(Int_t ipBin(1); ipBin<np; ipBin++){
1306     h = hBar->ProjectionY("npBin", ipBin, ipBin);
1307     if(!Int_t(h->Integral())) continue;
1308     h->Scale(100./h->Integral());
1309     Float_t p(ax->GetBinCenter(ipBin)); 
1310     Float_t dp(ax->GetBinWidth(ipBin)); 
1311     Int_t ip(g[0]->GetN());
1312     for(Int_t itl(AliTRDgeometry::kNlayer); itl--;){
1313       g[itl]->SetPoint(ip, p, h->GetBinContent(itl+1));
1314       g[itl]->SetPointError(ip, dp/3.46, h->GetBinError(itl+1));
1315     }
1316   }
1317
1318   TLegend *leg = new TLegend(0.76, 0.6, 1., 0.9);
1319   leg->SetBorderSize(0);
1320   leg->SetHeader("Tracklet/Track");
1321   leg->SetFillStyle(0);
1322   h = hBar->ProjectionX("npxBin"); h->Reset();
1323   h->SetTitle("");
1324   h->GetYaxis()->SetRangeUser(1., 99.);
1325   h->GetYaxis()->SetMoreLogLabels();
1326   h->GetYaxis()->CenterTitle();
1327   h->GetYaxis()->SetTitleOffset(1.2);
1328   h->SetYTitle("Prob. [%]");
1329   h->GetXaxis()->SetRangeUser(0.4, 12.);
1330   h->GetXaxis()->SetMoreLogLabels();
1331   h->GetXaxis()->CenterTitle();
1332   h->Draw("p");
1333   for(Int_t itl(AliTRDgeometry::kNlayer); itl--;){
1334     g[itl]->Draw("pc");
1335     leg->AddEntry(g[itl], Form("n = %d", itl+1),"pl");
1336   }
1337
1338   leg->Draw();
1339   gPad->SetLogx();gPad->SetLogy();
1340 }
1341
1342 //________________________________________________________
1343 Bool_t AliTRDcheckDET::MakePlotPulseHeight(){
1344   //
1345   // Create Plot of the Pluse Height Spectrum
1346   //
1347   TCanvas *output = gPad->GetCanvas();
1348   output->Divide(2);
1349   output->cd(1);
1350   TH1 *h, *h1, *h2;
1351   TObjArray *arr = (TObjArray*)fContainer->FindObject("<PH>");
1352   h = (TH1F*)arr->At(0);
1353   h->SetMarkerStyle(24);
1354   h->SetMarkerColor(kBlack);
1355   h->SetLineColor(kBlack);
1356   h->Draw("e1");
1357   // Trending for the pulse height: plateau value, slope and timebin of the maximum
1358   TLinearFitter fit(1,"pol1");
1359   Double_t time = 0.;
1360   for(Int_t itime = 10; itime <= 20; itime++){
1361     time = static_cast<Double_t>(itime);
1362     fit.AddPoint(&time, h->GetBinContent(itime + 1), h->GetBinError(itime + 1));
1363   }
1364   fit.Eval();
1365   Double_t plateau = fit.GetParameter(0) + 12 * fit.GetParameter(1);
1366   Double_t slope = fit.GetParameter(1);
1367   PutTrendValue("PHplateau", plateau);
1368   PutTrendValue("PHslope", slope);
1369   PutTrendValue("PHamplificationPeak", static_cast<Double_t>(h->GetMaximumBin()-1));
1370   AliDebug(1, Form("plateau %f, slope %f, MaxTime %f", plateau, slope, static_cast<Double_t>(h->GetMaximumBin()-1)));
1371 //   copy the second histogram in a new one with the same x-dimension as the phs with respect to time
1372   h1 = (TH1F *)arr->At(1);
1373   h2 = new TH1F("hphs1","Average PH", 31, -0.5, 30.5);
1374   for(Int_t ibin = h1->GetXaxis()->GetFirst(); ibin < h1->GetNbinsX(); ibin++) 
1375     h2->SetBinContent(ibin, h1->GetBinContent(ibin));
1376   h2->SetMarkerStyle(22);
1377   h2->SetMarkerColor(kBlue);
1378   h2->SetLineColor(kBlue);
1379   h2->Draw("e1same");
1380   gPad->Update();
1381 //   create axis according to the histogram dimensions of the original second histogram
1382   TGaxis *axis = new TGaxis(gPad->GetUxmin(),
1383                     gPad->GetUymax(),
1384                     gPad->GetUxmax(),
1385                     gPad->GetUymax(),
1386                     -0.08, 4.88, 510,"-L");
1387   axis->SetLineColor(kBlue);
1388   axis->SetLabelColor(kBlue);
1389   axis->SetTextColor(kBlue);
1390   axis->SetTitle("x_{0}-x_{c} [cm]");
1391   axis->Draw();
1392
1393   output->cd(2);
1394   TH2 *ph2d = (TH2F *)arr->At(2);
1395   ph2d->SetStats(kFALSE);
1396   ph2d->Draw("colz");
1397   return kTRUE;
1398 }
1399
1400 //________________________________________________________
1401 Bool_t AliTRDcheckDET::MakeBarPlot(TH1 *histo, Int_t color){
1402   //
1403   // Draw nice bar plots
1404   //
1405   if(!histo->GetEntries()) return kFALSE;
1406   histo->Scale(100./histo->Integral());
1407   histo->SetFillColor(color);
1408   histo->SetBarOffset(.2);
1409   histo->SetBarWidth(.6);
1410   histo->Draw("bar1");
1411   return kTRUE;
1412 }