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