8f2e69f148d50511bde9d99d0fb23a7102035a13
[u/mrichter/AliRoot.git] / TRD / qaRec / AliTRDcheckDET.cxx
1 #include <TAxis.h>
2 #include <TCanvas.h>
3 #include <TFile.h>
4 #include <TH1F.h>
5 #include <TH1I.h>
6 #include <TF1.h>
7 #include <TGaxis.h>
8 #include <TGraph.h>
9 #include <TLegend.h>
10 #include <TMath.h>
11 #include <TMap.h>
12 #include <TObjArray.h>
13 #include <TObject.h>
14 #include <TObjString.h>
15
16 #include <TPad.h>
17 #include <TProfile.h>
18 #include <TProfile2D.h>
19 #include <TROOT.h>
20
21 #include "AliLog.h"
22 #include "AliTRDcluster.h"
23 #include "AliESDHeader.h"
24 #include "AliESDRun.h"
25 #include "AliESDtrack.h"
26 #include "AliTRDgeometry.h"
27 #include "AliTRDpadPlane.h"
28 #include "AliTRDSimParam.h"
29 #include "AliTRDseedV1.h"
30 #include "AliTRDtrackV1.h"
31 #include "AliTRDtrackerV1.h"
32 #include "AliTRDReconstructor.h"
33 #include "AliTrackReference.h"
34 #include "AliTrackPointArray.h"
35 #include "AliTracker.h"
36 #include "TTreeStream.h"
37
38 #include "info/AliTRDtrackInfo.h"
39 #include "info/AliTRDeventInfo.h"
40 #include "AliTRDcheckDET.h"
41
42 #include <cstdio>
43 #include <iostream>
44
45 ////////////////////////////////////////////////////////////////////////////
46 //                                                                        //
47 //  Reconstruction QA                                                     //
48 //                                                                        //
49 //  Task doing basic checks for tracking and detector performance         //
50 //                                                                        //
51 //  Authors:                                                              //
52 //    Anton Andronic <A.Andronic@gsi.de>                                  //
53 //    Alexandru Bercuci <A.Bercuci@gsi.de>                                //
54 //    Markus Fasel <M.Fasel@gsi.de>                                       //
55 //                                                                        //
56 ////////////////////////////////////////////////////////////////////////////
57
58 //_______________________________________________________
59 AliTRDcheckDET::AliTRDcheckDET():
60   AliTRDrecoTask("checkDET", "Basic TRD data checker")
61   ,fEventInfo(0x0)
62   ,fTriggerNames(0x0)
63   ,fReconstructor(0x0)
64   ,fGeo(0x0)
65 {
66   //
67   // Default constructor
68   //
69   DefineInput(1,AliTRDeventInfo::Class());
70   fReconstructor = new AliTRDReconstructor;
71   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
72   fGeo = new AliTRDgeometry;
73   InitFunctorList();
74 }
75
76 //_______________________________________________________
77 AliTRDcheckDET::~AliTRDcheckDET(){
78   //
79   // Destructor
80   // 
81   if(fTriggerNames) delete fTriggerNames;
82   delete fReconstructor;
83   delete fGeo;
84 }
85
86 //_______________________________________________________
87 void AliTRDcheckDET::ConnectInputData(Option_t *opt){
88   //
89   // Connect the Input data with the task
90   //
91   AliTRDrecoTask::ConnectInputData(opt);
92   fEventInfo = dynamic_cast<AliTRDeventInfo *>(GetInputData(1));
93 }
94
95 //_______________________________________________________
96 void AliTRDcheckDET::CreateOutputObjects(){
97   //
98   // Create Output Objects
99   //
100   OpenFile(0,"RECREATE");
101   fContainer = Histos();
102   if(!fTriggerNames) fTriggerNames = new TMap();
103 }
104
105 //_______________________________________________________
106 void AliTRDcheckDET::Exec(Option_t *opt){
107   //
108   // Execution function
109   // Filling TRD quality histos
110   //
111   if(!HasMCdata() && fEventInfo->GetEventHeader()->GetEventType() != 7) return; // For real data we select only physical events
112   AliTRDrecoTask::Exec(opt);  
113   Int_t nTracks = 0;            // Count the number of tracks per event
114   Int_t triggermask = fEventInfo->GetEventHeader()->GetTriggerMask();
115   TString triggername =  fEventInfo->GetRunInfo()->GetFiredTriggerClasses(triggermask);
116   if(fDebugLevel > 6)printf("Trigger cluster: %d, Trigger class: %s\n", triggermask, triggername.Data());
117   dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger))->Fill(triggermask);
118   for(Int_t iti = 0; iti < fTracks->GetEntriesFast(); iti++){
119     if(!fTracks->UncheckedAt(iti)) continue;
120     AliTRDtrackInfo *fTrackInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iti));
121     if(!fTrackInfo->GetTrack()) continue;
122     nTracks++;
123   }
124   if(nTracks){
125     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks))->Fill(triggermask);
126     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtracksEvent))->Fill(nTracks);
127   }
128   if(triggermask <= 20 && !fTriggerNames->FindObject(Form("%d", triggermask))){
129     fTriggerNames->Add(new TObjString(Form("%d", triggermask)), new TObjString(triggername));
130     // also set the label for both histograms
131     TH1 *histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks));
132     histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
133     histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger));
134     histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
135   }
136   PostData(0, fContainer);
137 }
138
139 //_______________________________________________________
140 void AliTRDcheckDET::Terminate(Option_t *){
141   //
142   // Terminate function
143   //
144 }
145
146 //_______________________________________________________
147 Bool_t AliTRDcheckDET::PostProcess(){
148   //
149   // Do Postprocessing (for the moment set the number of Reference histograms)
150   //
151   
152   TH1 * h = 0x0;
153   
154   // Calculate of the trigger clusters purity
155   h = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTrigger"));
156   TH1F *h1 = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTriggerTracks"));
157   h1->Divide(h);
158   Float_t purities[20], val = 0;
159   TString triggernames[20];
160   Int_t nTriggerClasses = 0;
161   for(Int_t ibin = 1; ibin <= h->GetNbinsX(); ibin++){
162     if((val = h1->GetBinContent(ibin))){
163       purities[nTriggerClasses] = val;
164       triggernames[nTriggerClasses] = h1->GetXaxis()->GetBinLabel(ibin);
165       nTriggerClasses++;
166     }
167   }
168   h = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kTriggerPurity));
169   TAxis *ax = h->GetXaxis();
170   for(Int_t itrg = 0; itrg < nTriggerClasses; itrg++){
171     h->Fill(itrg, purities[itrg]);
172     ax->SetBinLabel(itrg+1, triggernames[itrg].Data());
173   }
174   ax->SetRangeUser(-0.5, nTriggerClasses+.5);
175   h->GetYaxis()->SetRangeUser(0,1);
176
177   fNRefFigures = 14;
178
179   return kTRUE;
180 }
181
182 //_______________________________________________________
183 Bool_t AliTRDcheckDET::GetRefFigure(Int_t ifig){
184   //
185   // Setting Reference Figures
186   //
187   gPad->SetLogy(0);
188   gPad->SetLogx(0);
189   switch(ifig){
190   case kNclustersTrack:
191     ((TH1F*)fContainer->FindObject("hNcls"))->Draw("pl");
192     return kTRUE;
193   case kNclustersTracklet:
194     ((TH1F*)fContainer->FindObject("hNclTls"))->Draw("pc");
195     return kTRUE;
196   case kNtrackletsTrack:
197     MakePlotNTracklets();
198     return kTRUE;
199   case kNtrackletsCross:
200     if(!MakeBarPlot((TH1F*)fContainer->FindObject("hNtlsCross"), kRed)) break;
201     return kTRUE;
202   case kNtrackletsFindable:
203     if(!MakeBarPlot((TH1F*)fContainer->FindObject("hNtlsFindable"), kGreen)) break;
204     return kTRUE;
205   case kNtracksEvent:
206     ((TH1F*)fContainer->FindObject("hNtrks"))->Draw("pl");
207     return kTRUE;
208   case kNtracksSector:
209     if(!MakeBarPlot((TH1F*)fContainer->FindObject("hNtrksSector"), kGreen)) break;
210     return kTRUE;
211   case kTrackStatus:
212     ((TH1I *)fContainer->FindObject("hTrackStatus"))->Draw("c");
213     return kTRUE;
214   case kTrackletStatus:
215     ((TH2S *)fContainer->FindObject("hTrackletStatus"))->Draw("colz");
216     return kTRUE;
217   case kChi2:
218     MakePlotChi2();
219     return kTRUE;
220   case kPH:
221     MakePlotPulseHeight();
222     return kTRUE;
223   case kChargeCluster:
224     ((TH1F*)fContainer->FindObject("hQcl"))->Draw("c");
225     gPad->SetLogy(1);
226     return kTRUE;
227   case kChargeTracklet:
228     ((TH1F*)fContainer->FindObject("hQtrklt"))->Draw("c");
229     return kTRUE;
230   case kNeventsTrigger:
231     ((TH1F*)fContainer->FindObject("hEventsTrigger"))->Draw("");
232     return kTRUE;
233   case kNeventsTriggerTracks:
234     ((TH1F*)fContainer->FindObject("hEventsTriggerTracks"))->Draw("");
235     return kTRUE;
236   case kTriggerPurity: 
237     if(!MakeBarPlot((TH1F*)fContainer->FindObject("hTriggerPurity"), kGreen)) break;
238     break;
239   default:
240     break;
241   }
242   AliInfo(Form("Reference plot [%d] missing result", ifig));
243   return kFALSE;
244 }
245
246 //_______________________________________________________
247 TObjArray *AliTRDcheckDET::Histos(){
248   //
249   // Create QA histograms
250   //
251   if(fContainer) return fContainer;
252   
253   fContainer = new TObjArray(20);
254   //fContainer->SetOwner(kTRUE);
255
256   // Register Histograms
257   TH1 * h = NULL;
258   TAxis *axis = NULL;
259   if(!(h = (TH1F *)gROOT->FindObject("hNcls"))){
260     h = new TH1F("hNcls", "N_{clusters} / track", 181, -0.5, 180.5);
261     h->GetXaxis()->SetTitle("N_{clusters}");
262     h->GetYaxis()->SetTitle("Entries");
263   } else h->Reset();
264   fContainer->AddAt(h, kNclustersTrack);
265
266   if(!(h = (TH1F *)gROOT->FindObject("hNclTls"))){
267     h = new TH1F("hNclTls","N_{clusters} / tracklet", 51, -0.5, 50.5);
268     h->GetXaxis()->SetTitle("N_{clusters}");
269     h->GetYaxis()->SetTitle("Entries");
270   } else h->Reset();
271   fContainer->AddAt(h, kNclustersTracklet);
272
273   if(!(h = (TH1F *)gROOT->FindObject("hNtls"))){
274     h = new TH1F("hNtls", "N_{tracklets} / track", AliTRDgeometry::kNlayer, 0.5, 6.5);
275     h->GetXaxis()->SetTitle("N^{tracklet}");
276     h->GetYaxis()->SetTitle("freq. [%]");
277   } else h->Reset();
278   fContainer->AddAt(h, kNtrackletsTrack);
279
280   if(!(h = (TH1F *)gROOT->FindObject("htlsSTA"))){
281     h = new TH1F("hNtlsSTA", "N_{tracklets} / track (Stand Alone)", AliTRDgeometry::kNlayer, 0.5, 6.5);
282     h->GetXaxis()->SetTitle("N^{tracklet}");
283     h->GetYaxis()->SetTitle("freq. [%]");
284   }
285   fContainer->AddAt(h, kNtrackletsSTA);
286
287   if(!(h = (TH1F *)gROOT->FindObject("htlsBAR"))){
288     h = new TH1F("hNtlsBAR", "N_{tracklets} / track (Barrel)", AliTRDgeometry::kNlayer, 0.5, 6.5);
289     h->GetXaxis()->SetTitle("N^{tracklet}");
290     h->GetYaxis()->SetTitle("freq. [%]");
291   }
292   fContainer->AddAt(h, kNtrackletsBAR);
293
294   // 
295   if(!(h = (TH1F *)gROOT->FindObject("hNtlsCross"))){
296     h = new TH1F("hNtlsCross", "N_{tracklets}^{cross} / track", 7, -0.5, 6.5);
297     h->GetXaxis()->SetTitle("n_{row cross}");
298     h->GetYaxis()->SetTitle("freq. [%]");
299   } else h->Reset();
300   fContainer->AddAt(h, kNtrackletsCross);
301
302   if(!(h = (TH1F *)gROOT->FindObject("hNtlsFindable"))){
303     h = new TH1F("hNtlsFindable", "Found/Findable Tracklets" , 101, -0.005, 1.005);
304     h->GetXaxis()->SetTitle("r [a.u]");
305     h->GetYaxis()->SetTitle("Entries");
306   } else h->Reset();
307   fContainer->AddAt(h, kNtrackletsFindable);
308
309   if(!(h = (TH1F *)gROOT->FindObject("hNtrks"))){
310     h = new TH1F("hNtrks", "N_{tracks} / event", 100, 0, 100);
311     h->GetXaxis()->SetTitle("N_{tracks}");
312     h->GetYaxis()->SetTitle("Entries");
313   } else h->Reset();
314   fContainer->AddAt(h, kNtracksEvent);
315
316   if(!(h = (TH1F *)gROOT->FindObject("hNtrksSector"))){
317     h = new TH1F("hNtrksSector", "N_{tracks} / sector", AliTRDgeometry::kNsector, -0.5, 17.5);
318     h->GetXaxis()->SetTitle("sector");
319     h->GetYaxis()->SetTitle("freq. [%]");
320   } else h->Reset();
321   fContainer->AddAt(h, kNtracksSector);
322
323   if(!(h = (TH1I *)gROOT->FindObject("hTrackStatus"))){
324     h = new TH1I("hTrackStatus", "Track Status", 7,0,7);
325     axis = h->GetXaxis();
326     axis->SetBinLabel(axis->GetFirst() + 0, "OK");
327     axis->SetBinLabel(axis->GetFirst() + 1, "PROL");
328     axis->SetBinLabel(axis->GetFirst() + 2, "PROP");
329     axis->SetBinLabel(axis->GetFirst() + 3, "ADJ");
330     axis->SetBinLabel(axis->GetFirst() + 4, "SNP");
331     axis->SetBinLabel(axis->GetFirst() + 5, "TLIN");
332     axis->SetBinLabel(axis->GetFirst() + 6, "UP");
333   }
334   fContainer->AddAt(h, kTrackStatus);
335
336   if(!(h = (TH2S *)gROOT->FindObject("hTrackletStatus"))){
337     h = new TH2S("hTrackletStatus", "Tracklet status", 6,0,6,8,0,8);
338     axis = h->GetYaxis();
339     axis->SetBinLabel(axis->GetFirst() + 0, "OK");
340     axis->SetBinLabel(axis->GetFirst() + 1, "Geom");
341     axis->SetBinLabel(axis->GetFirst() + 2, "Boun");
342     axis->SetBinLabel(axis->GetFirst() + 3, "NoCl");
343     axis->SetBinLabel(axis->GetFirst() + 4, "NoAttach");
344     axis->SetBinLabel(axis->GetFirst() + 5, "NoClTr");
345     axis->SetBinLabel(axis->GetFirst() + 6, "NoFit");
346     axis->SetBinLabel(axis->GetFirst() + 7, "Chi2");
347   }
348   fContainer->AddAt(h, kTrackletStatus);
349
350   // <PH> histos
351   TObjArray *arr = new TObjArray(2);
352   arr->SetOwner(kTRUE);  arr->SetName("<PH>");
353   fContainer->AddAt(arr, kPH);
354   if(!(h = (TH1F *)gROOT->FindObject("hPHt"))){
355     h = new TProfile("hPHt", "<PH>", 31, -0.5, 30.5);
356     h->GetXaxis()->SetTitle("Time / 100ns");
357     h->GetYaxis()->SetTitle("<PH> [a.u]");
358   } else h->Reset();
359   arr->AddAt(h, 0);
360   if(!(h = (TH1F *)gROOT->FindObject("hPHx")))
361     h = new TProfile("hPHx", "<PH>", 31, -0.08, 4.88);
362   else h->Reset();
363   arr->AddAt(h, 1);
364
365   // Chi2 histos
366   if(!(h = (TH2S*)gROOT->FindObject("hChi2"))){
367     h = new TH2S("hChi2", "#chi^{2} per track", AliTRDgeometry::kNlayer, .5, AliTRDgeometry::kNlayer+.5, 100, 0, 50);
368     h->SetXTitle("ndf");
369     h->SetYTitle("#chi^{2}/ndf");
370     h->SetZTitle("entries");
371   } else h->Reset();
372   fContainer->AddAt(h, kChi2);
373
374   if(!(h = (TH1F *)gROOT->FindObject("hQcl"))){
375     h = new TH1F("hQcl", "Q_{cluster}", 200, 0, 1200);
376     h->GetXaxis()->SetTitle("Q_{cluster} [a.u.]");
377     h->GetYaxis()->SetTitle("Entries");
378   }else h->Reset();
379   fContainer->AddAt(h, kChargeCluster);
380
381   if(!(h = (TH1F *)gROOT->FindObject("hQtrklt"))){
382     h = new TH1F("hQtrklt", "Q_{tracklet}", 6000, 0, 6000);
383     h->GetXaxis()->SetTitle("Q_{tracklet} [a.u.]");
384     h->GetYaxis()->SetTitle("Entries");
385   }else h->Reset();
386   fContainer->AddAt(h, kChargeTracklet);
387
388
389   if(!(h = (TH1F *)gROOT->FindObject("hEventsTrigger")))
390     h = new TH1F("hEventsTrigger", "Trigger Class", 100, 0, 100);
391   else h->Reset();
392   fContainer->AddAt(h, kNeventsTrigger);
393
394   if(!(h = (TH1F *)gROOT->FindObject("hEventsTriggerTracks")))
395     h = new TH1F("hEventsTriggerTracks", "Trigger Class (Tracks)", 100, 0, 100);
396   else h->Reset();
397   fContainer->AddAt(h, kNeventsTriggerTracks);
398
399   if(!(h = (TH1F *)gROOT->FindObject("hTriggerPurity"))){
400     h = new TH1F("hTriggerPurity", "Trigger Purity", 10, -0.5, 9.5);
401     h->GetXaxis()->SetTitle("Trigger Cluster");
402     h->GetYaxis()->SetTitle("freq.");
403   } else h->Reset();
404   fContainer->AddAt(h, kTriggerPurity);
405
406   return fContainer;
407 }
408
409 /*
410 * Plotting Functions
411 */
412
413 //_______________________________________________________
414 TH1 *AliTRDcheckDET::PlotTrackStatus(const AliTRDtrackV1 *track){
415   //
416   // Plot the track status
417   //
418   if(track) fTrack = track;
419   if(!fTrack){
420     AliWarning("No Track defined.");
421     return 0x0;
422   }
423   TH1 *h = 0x0;
424   if(!(h = dynamic_cast<TH1I *>(fContainer->At(kTrackStatus)))){
425     AliWarning("No Histogram defined.");
426     return 0x0;
427   }
428   h->Fill(fTrack->GetStatusTRD());
429   return h;
430 }
431
432 //_______________________________________________________
433 TH1 *AliTRDcheckDET::PlotTrackletStatus(const AliTRDtrackV1 *track){
434   //
435   // Plot the track status
436   //
437   if(track) fTrack = track;
438   if(!fTrack){
439     AliWarning("No Track defined.");
440     return 0x0;
441   }
442   TH1 *h = 0x0;
443   if(!(h = dynamic_cast<TH2S *>(fContainer->At(kTrackletStatus)))){
444     AliWarning("No Histogram defined.");
445     return 0x0;
446   }
447   UChar_t status = 0;
448   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
449     status = fTrack->GetStatusTRD(il);
450     h->Fill(il, status);
451   }
452   return h;
453 }
454
455 //_______________________________________________________
456 TH1 *AliTRDcheckDET::PlotNClustersTracklet(const AliTRDtrackV1 *track){
457   //
458   // Plot the mean number of clusters per tracklet
459   //
460   if(track) fTrack = track;
461   if(!fTrack){
462     AliWarning("No Track defined.");
463     return 0x0;
464   }
465   TH1 *h = 0x0;
466   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclustersTracklet)))){
467     AliWarning("No Histogram defined.");
468     return 0x0;
469   }
470   AliTRDseedV1 *tracklet = 0x0;
471   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
472     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
473     h->Fill(tracklet->GetN2());
474   }
475   return h;
476 }
477
478 //_______________________________________________________
479 TH1 *AliTRDcheckDET::PlotNClustersTrack(const AliTRDtrackV1 *track){
480   //
481   // Plot the number of clusters in one track
482   //
483   if(track) fTrack = track;
484   if(!fTrack){
485     AliWarning("No Track defined.");
486     return 0x0;
487   }
488   TH1 *h = 0x0;
489   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclustersTrack)))){
490     AliWarning("No Histogram defined.");
491     return 0x0;
492   }
493   
494   Int_t nclusters = 0;
495   AliTRDseedV1 *tracklet = 0x0;
496   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
497     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
498     nclusters += tracklet->GetN();
499     if(fDebugLevel > 2){
500       Int_t crossing = Int_t(tracklet->IsRowCross());
501       Int_t detector = tracklet->GetDetector();
502       Float_t theta = TMath::ATan(tracklet->GetZref(1));
503       Float_t phi = TMath::ATan(tracklet->GetYref(1));
504       Float_t momentum = 0.;
505       Int_t pdg = 0;
506       Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
507       UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
508       if(fMC){
509         if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
510         pdg = fMC->GetPDG();
511       }
512       (*fDebugStream) << "NClustersTrack"
513         << "Detector="  << detector
514         << "crossing="  << crossing
515         << "momentum="  << momentum
516         << "pdg="                               << pdg
517         << "theta="                     << theta
518         << "phi="                               << phi
519         << "kinkIndex=" << kinkIndex
520         << "TPCncls="           << TPCncls
521         << "nclusters=" << nclusters
522         << "\n";
523     }
524   }
525   h->Fill(nclusters);
526   return h;
527 }
528
529
530 //_______________________________________________________
531 TH1 *AliTRDcheckDET::PlotNTrackletsTrack(const AliTRDtrackV1 *track){
532   //
533   // Plot the number of tracklets
534   //
535   if(track) fTrack = track;
536   if(!fTrack){
537     AliWarning("No Track defined.");
538     return 0x0;
539   }
540   TH1 *h = 0x0, *hMethod = 0x0;
541   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsTrack)))){
542     AliWarning("No Histogram defined.");
543     return 0x0;
544   }
545   Int_t status = fESD->GetStatus();
546 /*  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);*/
547   if((status & AliESDtrack::kTRDin) != 0){
548     // Full BarrelTrack
549     if(!(hMethod = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsBAR))))
550       AliWarning("Method: Barrel.  Histogram not processed!");
551   } else {
552     // Stand alone Track
553     if(!(hMethod = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsSTA))))
554       AliWarning("Method: StandAlone.  Histogram not processed!");
555   }
556   Int_t nTracklets = fTrack->GetNumberOfTracklets();
557   h->Fill(nTracklets);
558   hMethod->Fill(nTracklets);
559   if(fDebugLevel > 3){
560     if(nTracklets == 1){
561       // If we have one Tracklet, check in which layer this happens
562       Int_t layer = -1;
563       AliTRDseedV1 *tracklet = 0x0;
564       for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
565         if((tracklet = fTrack->GetTracklet(il)) && tracklet->IsOK()){layer =  il; break;}
566       }
567       (*fDebugStream) << "NTrackletsTrack"
568         << "Layer=" << layer
569         << "\n";
570     }
571   }
572   return h;
573 }
574
575
576 //_______________________________________________________
577 TH1 *AliTRDcheckDET::PlotNTrackletsRowCross(const AliTRDtrackV1 *track){
578   //
579   // Plot the number of tracklets
580   //
581   if(track) fTrack = track;
582   if(!fTrack){
583     AliWarning("No Track defined.");
584     return 0x0;
585   }
586   TH1 *h = 0x0;
587   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsCross)))){
588     AliWarning("No Histogram defined.");
589     return 0x0;
590   }
591
592   Int_t ncross = 0;
593   AliTRDseedV1 *tracklet = 0x0;
594   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
595     if(!(tracklet = fTrack->GetTracklet(il)) || !tracklet->IsOK()) continue;
596
597     if(tracklet->IsRowCross()) ncross++;
598   }
599   h->Fill(ncross);
600   return h;
601 }
602
603 //_______________________________________________________
604 TH1 *AliTRDcheckDET::PlotFindableTracklets(const AliTRDtrackV1 *track){
605   //
606   // Plots the ratio of number of tracklets vs.
607   // number of findable tracklets
608   //
609   // Findable tracklets are defined as track prolongation
610   // to layer i does not hit the dead area +- epsilon
611   //
612   // In order to check whether tracklet hist active area in Layer i, 
613   // the track is refitted and the fitted position + an uncertainty 
614   // range is compared to the chamber border (also with a different
615   // uncertainty)
616   //
617   // For the track fit two cases are distinguished:
618   // If the track is a stand alone track (defined by the status bit 
619   // encoding, then the track is fitted with the tilted Rieman model
620   // Otherwise the track is fitted with the Kalman fitter in two steps:
621   // Since the track parameters are give at the outer point, we first 
622   // fit in direction inwards. Afterwards we fit again in direction outwards
623   // to extrapolate the track to layers which are not reached by the track
624   // For the Kalman model, the radial track points have to be shifted by
625   // a distance epsilon in the direction that we want to fit
626   //
627   const Float_t epsilon = 0.01;   // dead area tolerance
628   const Float_t epsilon_R = 1;    // shift in radial direction of the anode wire position (Kalman filter only)
629   const Float_t delta_y = 0.7;    // Tolerance in the track position in y-direction
630   const Float_t delta_z = 7.0;    // Tolerance in the track position in z-direction (Padlength)
631   Double_t x_anode[AliTRDgeometry::kNlayer] = {300.2, 312.8, 325.4, 338.0, 350.6, 363.2}; // Take the default X0
632  
633   if(track) fTrack = track;
634   if(!fTrack){
635     AliWarning("No Track defined.");
636     return 0x0;
637   }
638   TH1 *h = 0x0;
639   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsFindable)))){
640     AliWarning("No Histogram defined.");
641     return 0x0;
642   }
643   Int_t nFound = 0, nFindable = 0;
644   Int_t stack = -1;
645   Double_t ymin = 0., ymax = 0., zmin = 0., zmax = 0.;
646   Double_t y = 0., z = 0.;
647   AliTRDseedV1 *tracklet = 0x0;
648   AliTRDpadPlane *pp;  
649   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
650     if((tracklet = fTrack->GetTracklet(il)) && tracklet->IsOK()){
651       tracklet->SetReconstructor(fReconstructor);
652       nFound++;
653     }
654   }
655   // 2 Different cases:
656   // 1st stand alone: here we cannot propagate, but be can do a Tilted Rieman Fit
657   // 2nd barrel track: here we propagate the track to the layers
658   AliTrackPoint points[6];
659   Float_t xyz[3];
660   memset(xyz, 0, sizeof(Float_t) * 3);
661   if(((fESD->GetStatus() & AliESDtrack::kTRDout) > 0) && !((fESD->GetStatus() & AliESDtrack::kTRDin) > 0)){
662     // stand alone track
663     for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
664       xyz[0] = x_anode[il];
665       points[il].SetXYZ(xyz);
666     }
667     AliTRDtrackerV1::FitRiemanTilt(const_cast<AliTRDtrackV1 *>(fTrack), 0x0, kTRUE, 6, points);
668   } else {
669     // barrel track
670     //
671     // 2 Steps:
672     // -> Kalman inwards
673     // -> Kalman outwards
674     AliTRDtrackV1 copy_track(*fTrack);  // Do Kalman on a (non-constant) copy of the track
675     AliTrackPoint points_inward[6], points_outward[6];
676     for(Int_t il = AliTRDgeometry::kNlayer; il--;){
677       // In order to avoid complications in the Kalman filter if the track points have the same radial
678       // position like the tracklets, we have to shift the radial postion of the anode wire by epsilon
679       // in the direction we want to go
680       // The track points have to be in reverse order for the Kalman Filter inwards
681       xyz[0] = x_anode[AliTRDgeometry::kNlayer - il - 1] - epsilon_R;
682       points_inward[il].SetXYZ(xyz);
683       xyz[0] = x_anode[il] + epsilon_R;
684       points_outward[il].SetXYZ(xyz);
685     }
686     /*for(Int_t ipt = 0; ipt < AliTRDgeometry::kNlayer; ipt++)
687       printf("%d. X = %f\n", ipt, points[ipt].GetX());*/
688     // Kalman inwards
689     AliTRDtrackerV1::FitKalman(&copy_track, 0x0, kFALSE, 6, points_inward);
690     memcpy(points, points_inward, sizeof(AliTrackPoint) * 6); // Preliminary store the inward results in the Array points
691     // Kalman outwards
692     AliTRDtrackerV1::FitKalman(&copy_track, 0x0, kTRUE, 6, points_inward);
693     memcpy(points, points_outward, sizeof(AliTrackPoint) * AliTRDgeometry::kNlayer);
694   }
695   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
696     y = points[il].GetY();
697     z = points[il].GetZ();
698     if((stack = fGeo->GetStack(z, il)) < 0) continue; // Not findable
699     pp = fGeo->GetPadPlane(il, stack);
700     ymin = pp->GetCol0() + epsilon;
701     ymax = pp->GetColEnd() - epsilon; 
702     zmin = pp->GetRowEnd() + epsilon; 
703     zmax = pp->GetRow0() - epsilon;
704     // ignore y-crossing (material)
705     if((z + delta_z > zmin && z - delta_z < zmax) && (y + delta_y > ymin && y - delta_y < ymax)) nFindable++;
706       if(fDebugLevel > 3){
707         Double_t pos_tracklet[2] = {tracklet ? tracklet->GetYfit(0) : 0, tracklet ? tracklet->GetZfit(0) : 0};
708         Int_t hasTracklet = tracklet ? 1 : 0;
709         (*fDebugStream)   << "FindableTracklets"
710           << "layer="     << il
711           << "ytracklet=" << pos_tracklet[0]
712           << "ytrack="    << y
713           << "ztracklet=" << pos_tracklet[1]
714           << "ztrack="    << z
715           << "tracklet="  << hasTracklet
716           << "\n";
717       }
718   }
719   
720   h->Fill(nFindable > 0 ? TMath::Min(nFound/static_cast<Double_t>(nFindable), 1.) : 1);
721   if(fDebugLevel > 2) AliInfo(Form("Findable[Found]: %d[%d|%f]", nFindable, nFound, nFound/static_cast<Float_t>(nFindable > 0 ? nFindable : 1)));
722   return h;
723 }
724
725
726 //_______________________________________________________
727 TH1 *AliTRDcheckDET::PlotChi2(const AliTRDtrackV1 *track){
728   //
729   // Plot the chi2 of the track
730   //
731   if(track) fTrack = track;
732   if(!fTrack){
733     AliWarning("No Track defined.");
734     return 0x0;
735   }
736   TH1 *h = 0x0;
737   if(!(h = dynamic_cast<TH2S*>(fContainer->At(kChi2)))) {
738     AliWarning("No Histogram defined.");
739     return 0x0;
740   }
741   Int_t n = fTrack->GetNumberOfTracklets();
742   if(!n) return 0x0;
743
744   h->Fill(n, fTrack->GetChi2()/n);
745   return h;
746 }
747
748
749 //_______________________________________________________
750 TH1 *AliTRDcheckDET::PlotPHt(const AliTRDtrackV1 *track){
751   //
752   // Plot the average pulse height
753   //
754   if(track) fTrack = track;
755   if(!fTrack){
756     AliWarning("No Track defined.");
757     return 0x0;
758   }
759   TProfile *h = 0x0;
760   if(!(h = dynamic_cast<TProfile *>(((TObjArray*)(fContainer->At(kPH)))->At(0)))){
761     AliWarning("No Histogram defined.");
762     return 0x0;
763   }
764   AliTRDseedV1 *tracklet = 0x0;
765   AliTRDcluster *c = 0x0;
766   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
767     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
768     Int_t crossing = Int_t(tracklet->IsRowCross());
769     Int_t detector = tracklet->GetDetector();
770     tracklet->ResetClusterIter();
771     while((c = tracklet->NextCluster())){
772       if(!c->IsInChamber()) continue;
773       Int_t localtime        = c->GetLocalTimeBin();
774       Double_t absolute_charge = TMath::Abs(c->GetQ());
775       h->Fill(localtime, absolute_charge);
776       if(fDebugLevel > 3){
777         Double_t distance[2];
778         GetDistanceToTracklet(distance, tracklet, c);
779         Float_t theta = TMath::ATan(tracklet->GetZref(1));
780         Float_t phi = TMath::ATan(tracklet->GetYref(1));
781         Float_t momentum = 0.;
782         Int_t pdg = 0;
783         Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
784         UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
785         if(fMC){
786           if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
787           pdg = fMC->GetPDG();
788         }
789         (*fDebugStream) << "PHt"
790           << "Detector="        << detector
791           << "crossing="        << crossing
792           << "Timebin="         << localtime
793           << "Charge="          << absolute_charge
794           << "momentum="        << momentum
795           << "pdg="                             << pdg
796           << "theta="                   << theta
797           << "phi="                             << phi
798           << "kinkIndex="       << kinkIndex
799           << "TPCncls="         << TPCncls
800           << "dy="        << distance[0]
801           << "dz="        << distance[1]
802           << "c.="        << c
803           << "\n";
804       }
805     }
806   }
807   return h;
808 }
809
810 //_______________________________________________________
811 TH1 *AliTRDcheckDET::PlotPHx(const AliTRDtrackV1 *track){
812   //
813   // Plots the average pulse height vs the distance from the anode wire
814   // (plus const anode wire offset)
815   //
816   if(track) fTrack = track;
817   if(!fTrack){
818     AliWarning("No Track defined.");
819     return 0x0;
820   }
821   TProfile *h = 0x0;
822   if(!(h = dynamic_cast<TProfile *>(((TObjArray*)(fContainer->At(kPH)))->At(1)))){
823     AliWarning("No Histogram defined.");
824     return 0x0;
825   }
826   Float_t offset = .5*AliTRDgeometry::CamHght();
827   AliTRDseedV1 *tracklet = 0x0;
828   AliTRDcluster *c = 0x0;
829   Double_t distance = 0;
830   Double_t x, y;
831   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
832     if(!(tracklet = fTrack->GetTracklet(itl)) || !(tracklet->IsOK())) continue;
833     tracklet->ResetClusterIter();
834     while((c = tracklet->NextCluster())){
835       if(!c->IsInChamber()) continue;
836       x = c->GetX()-AliTRDcluster::GetXcorr(c->GetLocalTimeBin());
837       y = c->GetY()-AliTRDcluster::GetYcorr(AliTRDgeometry::GetLayer(c->GetDetector()), c->GetCenter());
838
839       distance = tracklet->GetX0() - (c->GetX() + 0.3) + offset;
840       h->Fill(distance, TMath::Abs(c->GetQ()));
841     }
842   }  
843   return h;
844 }
845
846 //_______________________________________________________
847 TH1 *AliTRDcheckDET::PlotChargeCluster(const AliTRDtrackV1 *track){
848   //
849   // Plot the cluster charge
850   //
851   if(track) fTrack = track;
852   if(!fTrack){
853     AliWarning("No Track defined.");
854     return 0x0;
855   }
856   TH1 *h = 0x0;
857   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeCluster)))){
858     AliWarning("No Histogram defined.");
859     return 0x0;
860   }
861   AliTRDseedV1 *tracklet = 0x0;
862   AliTRDcluster *c = 0x0;
863   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
864     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
865     for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
866       if(!(c = tracklet->GetClusters(itime))) continue;
867       h->Fill(c->GetQ());
868     }
869   }
870   return h;
871 }
872
873 //_______________________________________________________
874 TH1 *AliTRDcheckDET::PlotChargeTracklet(const AliTRDtrackV1 *track){
875   //
876   // Plot the charge deposit per chamber
877   //
878   if(track) fTrack = track;
879   if(!fTrack){
880     AliWarning("No Track defined.");
881     return 0x0;
882   }
883   TH1 *h = 0x0;
884   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeTracklet)))){
885     AliWarning("No Histogram defined.");
886     return 0x0;
887   }
888   AliTRDseedV1 *tracklet = 0x0;
889   AliTRDcluster *c = 0x0;
890   Double_t Qtot = 0;
891   Int_t nTracklets =fTrack->GetNumberOfTracklets();
892   for(Int_t itl = 0x0; itl < AliTRDgeometry::kNlayer; itl++){
893     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
894     Qtot = 0.;
895     for(Int_t ic = AliTRDseedV1::kNclusters; ic--;){
896       if(!(c = tracklet->GetClusters(ic))) continue;
897       Qtot += TMath::Abs(c->GetQ());
898     }
899     h->Fill(Qtot);
900     if(fDebugLevel > 3){
901       Int_t crossing = (Int_t)tracklet->IsRowCross();
902       Int_t detector = tracklet->GetDetector();
903       Float_t theta = TMath::ATan(tracklet->GetZfit(1));
904       Float_t phi = TMath::ATan(tracklet->GetYfit(1));
905       Float_t momentum = 0.;
906       Int_t pdg = 0;
907       Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
908       UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
909       if(fMC){
910               if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
911         pdg = fMC->GetPDG();
912       }
913       (*fDebugStream) << "ChargeTracklet"
914         << "Detector="  << detector
915         << "crossing="  << crossing
916         << "momentum="  << momentum
917         << "nTracklets="<< nTracklets
918         << "pdg="                               << pdg
919         << "theta="                     << theta
920         << "phi="                               << phi
921         << "kinkIndex=" << kinkIndex
922         << "TPCncls="           << TPCncls
923         << "QT="        << Qtot
924         << "\n";
925     }
926   }
927   return h;
928 }
929
930 //_______________________________________________________
931 TH1 *AliTRDcheckDET::PlotNTracksSector(const AliTRDtrackV1 *track){
932   //
933   // Plot the number of tracks per Sector
934   //
935   if(track) fTrack = track;
936   if(!fTrack){
937     AliWarning("No Track defined.");
938     return 0x0;
939   }
940   TH1 *h = 0x0;
941   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtracksSector)))){
942     AliWarning("No Histogram defined.");
943     return 0x0;
944   }
945
946   // TODO we should compare with
947   // sector = Int_t(track->GetAlpha() / AliTRDgeometry::GetAlpha());
948
949   AliTRDseedV1 *tracklet = 0x0;
950   Int_t sector = -1;
951   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
952     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
953     sector = static_cast<Int_t>(tracklet->GetDetector()/AliTRDgeometry::kNdets);
954     break;
955   }
956   h->Fill(sector);
957   return h;
958 }
959
960
961 //________________________________________________________
962 void AliTRDcheckDET::SetRecoParam(AliTRDrecoParam *r)
963 {
964
965   fReconstructor->SetRecoParam(r);
966 }
967
968 //________________________________________________________
969 void AliTRDcheckDET::GetDistanceToTracklet(Double_t *dist, AliTRDseedV1 *tracklet, AliTRDcluster *c)
970 {
971   Float_t x = c->GetX();
972   dist[0] = c->GetY() - tracklet->GetYat(x);
973   dist[1] = c->GetZ() - tracklet->GetZat(x);
974 }
975
976
977 //_______________________________________________________
978 void AliTRDcheckDET::MakePlotChi2()
979 {
980 // Plot chi2/track normalized to number of degree of freedom 
981 // (tracklets) and compare with the theoretical distribution.
982 // 
983 // Alex Bercuci <A.Bercuci@gsi.de>
984
985   TH2S *h2 = (TH2S*)fContainer->At(kChi2);
986   TF1 f("fChi2", "[0]*pow(x, [1]-1)*exp(-0.5*x)", 0., 50.);
987   TLegend *leg = new TLegend(.7,.7,.95,.95);
988   leg->SetBorderSize(1); leg->SetHeader("Tracklets per Track");
989   Bool_t kFIRST = kTRUE;
990   for(Int_t il=1; il<=h2->GetNbinsX(); il++){
991     TH1D *h1 = h2->ProjectionY(Form("pyChi2%d", il), il, il);
992     if(h1->Integral()<50) continue;
993     h1->Scale(1./h1->Integral());
994     h1->SetMarkerStyle(7);h1->SetMarkerColor(il);
995     h1->SetLineColor(il);h1->SetLineStyle(2);
996     f.SetParameter(1, .5*il);f.SetLineColor(il);
997     h1->Fit(&f, "QW+", kFIRST ? "pc": "pcsame");
998     leg->AddEntry(h1, Form("%d", il), "l");
999     if(kFIRST){
1000       h1->GetXaxis()->SetRangeUser(0., 25.);
1001     }
1002     kFIRST = kFALSE;
1003   }
1004   leg->Draw();
1005   gPad->SetLogy();
1006 }
1007
1008
1009 //________________________________________________________
1010 void AliTRDcheckDET::MakePlotNTracklets(){
1011   //
1012   // Make nice bar plot of the number of tracklets in each method
1013   //
1014   TH1F *hBAR = (TH1F *)fContainer->FindObject("hNtlsBAR");
1015   TH1F *hSTA = (TH1F *)fContainer->FindObject("hNtlsSTA");
1016   TH1F *hCON = (TH1F *)fContainer->FindObject("hNtls");
1017   TLegend *leg = new TLegend(0.6, 0.75, 0.89, 0.89);
1018
1019   Float_t scale = hCON->Integral();
1020   hCON->Scale(100./scale);
1021   hCON->SetFillColor(kRed);hCON->SetLineColor(kRed);
1022   hCON->SetBarWidth(0.2);
1023   hCON->SetBarOffset(0.6);
1024   hCON->SetStats(kFALSE);
1025   hCON->GetYaxis()->SetRangeUser(0.,40.);
1026   hCON->GetYaxis()->SetTitleOffset(1.2);
1027   hCON->Draw("bar1"); leg->AddEntry(hCON, "Total", "f");
1028
1029   hBAR->Scale(100./scale);
1030   hBAR->SetFillColor(kGreen);hBAR->SetLineColor(kGreen);
1031   hBAR->SetBarWidth(0.2);
1032   hBAR->SetBarOffset(0.2);
1033   hBAR->SetTitle("");
1034   hBAR->SetStats(kFALSE);
1035   hBAR->GetYaxis()->SetRangeUser(0.,40.);
1036   hBAR->GetYaxis()->SetTitleOffset(1.2);
1037   hBAR->Draw("bar1same"); leg->AddEntry(hBAR, "Barrel", "f");
1038
1039   hSTA->Scale(100./scale);
1040   hSTA->SetFillColor(kBlue);hSTA->SetLineColor(kBlue);
1041   hSTA->SetBarWidth(0.2);
1042   hSTA->SetBarOffset(0.4);
1043   hSTA->SetTitle("");
1044   hSTA->SetStats(kFALSE);
1045   hSTA->GetYaxis()->SetRangeUser(0.,40.);
1046   hSTA->GetYaxis()->SetTitleOffset(1.2);
1047   hSTA->Draw("bar1same"); leg->AddEntry(hSTA, "Stand Alone", "f");
1048   leg->Draw();
1049   gPad->Update();
1050 }
1051
1052 //________________________________________________________
1053 void AliTRDcheckDET::MakePlotPulseHeight(){
1054   //
1055   // Create Plot of the Pluse Height Spectrum
1056   //
1057   TH1 *h, *h1, *h2;
1058   TObjArray *arr = (TObjArray*)fContainer->FindObject("<PH>");
1059   h = (TH1F*)arr->At(0);
1060   h->SetMarkerStyle(24);
1061   h->SetMarkerColor(kBlack);
1062   h->SetLineColor(kBlack);
1063   h->Draw("e1");
1064 //   copy the second histogram in a new one with the same x-dimension as the phs with respect to time
1065   h1 = (TH1F *)arr->At(1);
1066   h2 = new TH1F("hphs1","Average PH", 31, -0.5, 30.5);
1067   for(Int_t ibin = h1->GetXaxis()->GetFirst(); ibin < h1->GetNbinsX(); ibin++) 
1068     h2->SetBinContent(ibin, h1->GetBinContent(ibin));
1069   h2->SetMarkerStyle(22);
1070   h2->SetMarkerColor(kBlue);
1071   h2->SetLineColor(kBlue);
1072   h2->Draw("e1same");
1073   gPad->Update();
1074 //   create axis according to the histogram dimensions of the original second histogram
1075   TGaxis *axis = new TGaxis(gPad->GetUxmin(),
1076                     gPad->GetUymax(),
1077                     gPad->GetUxmax(),
1078                     gPad->GetUymax(),
1079                     -0.08, 4.88, 510,"-L");
1080   axis->SetLineColor(kBlue);
1081   axis->SetLabelColor(kBlue);
1082   axis->SetTextColor(kBlue);
1083   axis->SetTitle("x_{0}-x_{c} [cm]");
1084   axis->Draw();
1085 }
1086
1087 //________________________________________________________
1088 Bool_t AliTRDcheckDET::MakeBarPlot(TH1 *histo, Int_t color){
1089   //
1090   // Draw nice bar plots
1091   //
1092   if(!histo->GetEntries()) return kFALSE;
1093   histo->Scale(100./histo->Integral());
1094   histo->SetFillColor(color);
1095   histo->SetBarOffset(.2);
1096   histo->SetBarWidth(.6);
1097   histo->Draw("bar1");
1098   return kTRUE;
1099 }