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