Updates for better DQM.
[u/mrichter/AliRoot.git] / FMD / AliFMDQAChecker.cxx
1 /**************************************************************************
2  * Copyright(c) 2004, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15 //__________________________________________________________________
16 //
17 // Yves?
18 // What 
19 // is 
20 // this 
21 // class 
22 // supposed 
23 // to
24 // do?
25 //__________________________________________________________________
26 //
27 // --- ROOT system ---
28 #include <TClass.h>
29 #include <TH1F.h> 
30 #include <TF1.h> 
31 #include <TH2.h> 
32 #include <TH1I.h> 
33 #include <TIterator.h> 
34 #include <TKey.h> 
35 #include <TFile.h> 
36 #include <iostream>
37 #include <TCanvas.h>
38 #include <TPaveText.h>
39 #include <TStyle.h>
40 #include <TLatex.h>
41 #include <TFitResult.h>
42 #include <TParameter.h>
43 #include <TMacro.h>
44
45 // --- AliRoot header files ---
46 #include "AliLog.h"
47 #include "AliQAv1.h"
48 #include "AliQAChecker.h"
49 #include "AliFMDQAChecker.h"
50 #include "AliRecoParam.h"
51 #include <AliCDBManager.h>
52 #include <AliCDBEntry.h>
53 #include <AliCDBId.h>
54 #include <AliQAThresholds.h>
55
56 ClassImp(AliFMDQAChecker)
57 #if 0
58 ; // This is for Emacs! - do not delete
59 #endif
60
61 namespace {
62   void addFitsMacro(TList* l) {
63     TMacro* m = new TMacro("fits");
64     m->AddLine("void fits() {");
65     m->AddLine("  if (!gPad) { Printf(\"No gPad\"); return; }");
66     m->AddLine("  TList* lp = gPad->GetListOfPrimitives();");
67     m->AddLine("  if (!lp) return;");
68     m->AddLine("  TObject* po = 0;");
69     m->AddLine("  TIter next(lp);");
70     m->AddLine("  while ((po = next())) {");
71     m->AddLine("    if (!po->IsA()->InheritsFrom(TH1::Class())) continue;");
72     m->AddLine("    TH1*   htmp = dynamic_cast<TH1*>(po);");
73     m->AddLine("    TList* lf   = htmp->GetListOfFunctions();");
74     m->AddLine("    TObject* pso = (lf ? lf->FindObject(\"stats\") : 0);");
75     m->AddLine("    if (!pso) continue;");
76     m->AddLine("    TPaveStats* ps = static_cast<TPaveStats*>(pso);");
77     m->AddLine("    ps->SetOptFit(111);");
78     m->AddLine("    UShort_t qual = htmp->GetUniqueID();");
79     m->AddLine("    ps->SetFillColor(qual >= 3 ? kRed-4 : qual >= 2 ? kOrange-4 : qual >= 1 ? kYellow-4 : kGreen-4);");
80     // m->AddLine("    lf->Remove(lf->FindObject(\"fits\"));");
81     // m->AddLine("    ps->Paint();");
82     m->AddLine("    break;");
83     m->AddLine("  }");
84     // m->AddLine("  gPad->Modified(); gPad->Update(); gPad->cd();");
85     m->AddLine("}");
86
87     TObject* old = l->FindObject(m->GetName());
88     if (old) l->Remove(old);
89     l->Add(m);
90   }
91   
92   const Double_t kROErrorsLabelY    = 30.;
93   
94   const Int_t    kConvolutionSteps  = 100;
95   const Double_t kConvolutionNSigma = 5;
96
97   // 
98   // The shift of the most probable value for the ROOT function TMath::Landau 
99   //
100   const Double_t  kMpShift  = -0.22278298;
101   // 
102   // Integration normalisation 
103   //
104   const Double_t  kInvSq2Pi = 1. / TMath::Sqrt(2*TMath::Pi());
105
106   Double_t landau(Double_t x, Double_t delta, Double_t xi)
107   {
108     // 
109     // Calculate the shifted Landau
110     // @f[
111     //    f'_{L}(x;\Delta,\xi) = f_L(x;\Delta+0.22278298\xi)
112     // @f]
113     //
114     // where @f$ f_{L}@f$ is the ROOT implementation of the Landau
115     // distribution (known to have @f$ \Delta_{p}=-0.22278298@f$ for
116     // @f$\Delta=0,\xi=1@f$.
117     // 
118     // Parameters:
119     //    x      Where to evaluate @f$ f'_{L}@f$ 
120     //    delta  Most probable value 
121     //    xi     The 'width' of the distribution 
122     //
123     // Return:
124     //    @f$ f'_{L}(x;\Delta,\xi) @f$
125     //
126     return TMath::Landau(x, delta - xi * kMpShift, xi);
127   }
128   Double_t landauGaus(Double_t x, Double_t delta, Double_t xi,
129                       Double_t sigma, Double_t sigmaN)
130   {
131     // 
132     // Calculate the value of a Landau convolved with a Gaussian 
133     // 
134     // @f[ 
135     // f(x;\Delta,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}}
136     //    \int_{-\infty}^{+\infty} d\Delta' f'_{L}(x;\Delta',\xi)
137     //    \exp{-\frac{(\Delta-\Delta')^2}{2\sigma'^2}}
138     // @f]
139     // 
140     // where @f$ f'_{L}@f$ is the Landau distribution, @f$ \Delta@f$
141     // the energy loss, @f$ \xi@f$ the width of the Landau, and @f$
142     // \sigma'^2=\sigma^2-\sigma_n^2 @f$.  Here, @f$\sigma@f$ is the
143     // variance of the Gaussian, and @f$\sigma_n@f$ is a parameter
144     // modelling noise in the detector.
145     //
146     // Note that this function uses the constants kConvolutionSteps and
147     // kConvolutionNSigma
148     // 
149     // References: 
150     //  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
151     //  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
152     //  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
153     // 
154     // Parameters:
155     //    x         where to evaluate @f$ f@f$
156     //    delta     @f$ \Delta@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
157     //    xi        @f$ \xi@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
158     //    sigma     @f$ \sigma@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
159     //    sigma_n   @f$ \sigma_n@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
160     // 
161     // Return:
162     //    @f$ f@f$ evaluated at @f$ x@f$.  
163     //
164     Double_t deltap = delta - xi * kMpShift;
165     Double_t sigma2 = sigmaN*sigmaN + sigma*sigma;
166     Double_t sigma1 = sigmaN == 0 ? sigma : TMath::Sqrt(sigma2);
167     Double_t xlow   = x - kConvolutionNSigma * sigma1;
168     Double_t xhigh  = x + kConvolutionNSigma * sigma1;
169     Double_t step   = (xhigh - xlow) / kConvolutionSteps;
170     Double_t sum    = 0;
171   
172     for (Int_t i = 0; i <= kConvolutionSteps/2; i++) { 
173       Double_t x1 = xlow  + (i - .5) * step;
174       Double_t x2 = xhigh - (i - .5) * step;
175     
176       sum += TMath::Landau(x1, deltap, xi, kTRUE) * TMath::Gaus(x, x1, sigma1);
177       sum += TMath::Landau(x2, deltap, xi, kTRUE) * TMath::Gaus(x, x2, sigma1);
178     }
179     return step * sum * kInvSq2Pi / sigma1;
180   }
181
182   // 
183   // Utility function to use in TF1 defintition 
184   //
185   Double_t landauGaus1(Double_t* xp, Double_t* pp) 
186   {
187     Double_t x        = xp[0];
188     Double_t constant = pp[0];
189     Double_t delta    = pp[1];
190     Double_t xi       = pp[2];
191     Double_t sigma    = pp[3];
192     Double_t sigmaN   = pp[4];
193
194     return constant * landauGaus(x, delta, xi, sigma, sigmaN);
195   }
196
197   //____________________________________________________________________
198   TF1* makeLandauGaus(const char* , 
199                       Double_t  c=1, 
200                       Double_t  delta=.5,  Double_t xi=0.07, 
201                       Double_t  sigma=.1,  Double_t sigmaN=-1, 
202                       Double_t  xmin=0,    Double_t xmax=15)
203   {
204     // 
205     // Generate a TF1 object of @f$ f_I@f$ 
206     // 
207     // Parameters:
208     //    c        Constant
209     //    delta    @f$ \Delta@f$ 
210     //    xi       @f$ \xi_1@f$        
211     //    sigma    @f$ \sigma_1@f$             
212     //    sigma_n  @f$ \sigma_n@f$             
213     //    xmin     Least value of range
214     //    xmax     Largest value of range
215     // 
216     // Return:
217     //    Newly allocated TF1 object
218     //
219     Int_t npar     = 5;
220     TF1*  func     = new TF1("landauGaus", 
221                              &landauGaus1,xmin,xmax,npar);
222     // func->SetLineStyle(((i-2) % 10)+2); // start at dashed
223     func->SetLineColor(kBlack); 
224     func->SetLineWidth(2);
225     func->SetNpx(500);
226     func->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
227
228     // Set the initial parameters from the seed fit 
229     func->SetParameter(0,      c);       
230     func->SetParameter(1,  delta);   
231     func->SetParameter(2,     xi);      
232     func->SetParameter(3,  sigma);   
233     func->SetParameter(4, sigmaN); 
234
235     func->SetParLimits(1, 0,    xmax);
236     func->SetParLimits(2, 0,    xmax);
237     func->SetParLimits(3, 0.01, 0.4);
238     
239     if (sigmaN < 0) func->FixParameter(4, 0);
240     else            func->SetParLimits(4, 0, xmax);
241
242     return func;
243   }
244 }
245   
246 //__________________________________________________________________
247 AliFMDQAChecker::AliFMDQAChecker() 
248   : AliQACheckerBase("FMD","FMD Quality Assurance Checker") ,
249     fDoScale(false),
250     fDidExternal(false), 
251     fShowFitResults(true),
252     fELossLowCut(0.2), 
253     fELossNRMS(3), 
254     fELossBadChi2Nu(10), 
255     fELossFkupChi2Nu(100), 
256     fELossMinEntries(1000),
257     fELossGoodParError(0.1),
258     fROErrorsBad(0.3), 
259     fROErrorsFkup(0.5)
260 {
261 }          
262
263 //__________________________________________________________________
264 void
265 AliFMDQAChecker::ProcessExternalParams()
266 {
267   ProcessExternalParam("ELossLowCut",           fELossLowCut);
268   ProcessExternalParam("ELossNRMS",             fELossNRMS);
269   ProcessExternalParam("ELossBadChi2Nu",        fELossBadChi2Nu);
270   ProcessExternalParam("ELossFkupChi2Nu",       fELossFkupChi2Nu);
271   ProcessExternalParam("ELossGoodParError",     fELossGoodParError);
272   ProcessExternalParam("ROErrorsBad",           fROErrorsBad);
273   ProcessExternalParam("ROErrorsFkup",          fROErrorsFkup);
274   Double_t tmp = 0;
275   ProcessExternalParam("CommonScale",           tmp);
276   fDoScale = tmp > 0;
277   ProcessExternalParam("ELossMinEntries",       tmp);
278   fELossMinEntries = tmp;
279
280   GetThresholds();
281
282   fDidExternal = true;
283 }
284 //__________________________________________________________________
285 void
286 AliFMDQAChecker::ProcessExternalParam(const char* name, Double_t& v)
287 {
288   TObject* o = fExternParamList->FindObject(name);
289   if (!o) return; 
290   TParameter<double>* p = static_cast<TParameter<double>*>(o);
291   v = p->GetVal();
292   AliDebugF(3, "External parameter: %-20s=%lf", name, v);
293 }
294
295 //__________________________________________________________________
296 void
297 AliFMDQAChecker::GetThresholds()
298 {
299   const char*    path   = "GRP/Calib/QAThresholds";
300   AliCDBManager* cdbMan = AliCDBManager::Instance();
301   AliCDBEntry*   cdbEnt = cdbMan->Get(path);
302   if (!cdbEnt) { 
303     AliWarningF("Failed to get CDB entry at %s", path);
304     return;
305   }
306   
307   TObjArray*     cdbObj = static_cast<TObjArray*>(cdbEnt->GetObject());
308   if (!cdbObj) { 
309     AliWarningF("Failed to get CDB object at %s", path);
310     return;
311   }
312   
313   TObject*       fmdObj = cdbObj->FindObject("FMD");
314   if (!fmdObj) { 
315     AliWarningF("Failed to get FMD object at from CDB %s", path);
316     return;
317   }
318   
319   AliQAThresholds* qaThr = static_cast<AliQAThresholds*>(fmdObj);
320   Int_t nThr = qaThr->GetSize();
321   for (Int_t i = 0; i < nThr; i++) { 
322     TObject* thr = qaThr->GetThreshold(i);
323     if (!thr) continue; 
324
325     TParameter<double>* d = dynamic_cast<TParameter<double>*>(thr);
326     if (!d) { 
327       AliWarningF("Parameter %s not of type double", thr->GetName());
328       continue;
329     }
330     Double_t val = d->GetVal();
331     TString  name(thr->GetName());
332     if      (name.EqualTo("ELossBadChi2Nu"))     fELossBadChi2Nu    = val;
333     else if (name.EqualTo("ELossFkupChi2Nu"))    fELossFkupChi2Nu   = val;
334     else if (name.EqualTo("ELossGoodParError"))  fELossGoodParError = val;
335     else if (name.EqualTo("ROErrorsBad"))        fROErrorsBad       = val;
336     else if (name.EqualTo("ROErrorsFkup"))       fROErrorsFkup      = val;    
337     AliDebugF(3, "Threshold %s=%f", name.Data(), val);
338   }
339 }
340
341 //__________________________________________________________________
342 AliQAv1::QABIT_t
343 AliFMDQAChecker::Quality2Bit(UShort_t value) const
344 {
345   AliQAv1::QABIT_t  ret   = AliQAv1::kINFO; // Assume success 
346   if      (value >= kWhatTheFk) ret = AliQAv1::kFATAL;
347   else if (value >= kBad)       ret = AliQAv1::kERROR;
348   else if (value >= kProblem)   ret = AliQAv1::kWARNING;
349   
350   return ret;
351 }
352
353 //__________________________________________________________________
354 void
355 AliFMDQAChecker::SetQA(AliQAv1::ALITASK_t index, Double_t* values) const
356 {
357   AliQAv1 * qa = AliQAv1::Instance(index) ;
358
359   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) {
360     // Check if specie is defined 
361     if (!qa->IsEventSpecieSet(AliRecoParam::ConvertIndex(specie)))
362       continue ;
363     
364     // No checker is implemented, set all QA to Fatal
365     if (!values) { 
366       qa->Set(AliQAv1::kFATAL, specie) ; 
367       continue;
368     }
369     
370     UShort_t          value = values[specie];
371     AliQAv1::QABIT_t  ret   = Quality2Bit(value);
372     
373     qa->Set(ret, AliRecoParam::ConvertIndex(specie));
374     AliDebugF(3, "Quality of %s: %d -> %d", 
375               AliRecoParam::GetEventSpecieName(specie), value, ret);
376   }
377 }
378
379 //__________________________________________________________________
380 UShort_t
381 AliFMDQAChecker::BasicCheck(TH1* hist) const
382 {
383   if (hist->GetEntries() <= 0) return kOK;
384   return (hist->GetMean() > 0 ? kOK : kProblem);
385 }
386
387 //__________________________________________________________________
388 UShort_t
389 AliFMDQAChecker::CheckOne(AliQAv1::ALITASK_t          what,
390                           AliRecoParam::EventSpecie_t specie, 
391                           TH1*                        hist) const
392 {
393   if(what == AliQAv1::kESD) return CheckESD(specie, hist);
394   if(what == AliQAv1::kRAW) return CheckRaw(specie, hist);
395   if(what == AliQAv1::kSIM) return CheckSim(specie, hist);
396   if(what == AliQAv1::kREC) return CheckRec(specie, hist);
397   return 0;
398 }
399 //__________________________________________________________________
400 UShort_t 
401 AliFMDQAChecker::CheckESD(AliRecoParam::EventSpecie_t /* specie*/, 
402                           TH1*                        hist) const
403 {
404   return BasicCheck(hist);
405 }
406 //__________________________________________________________________
407 UShort_t
408 AliFMDQAChecker::CheckFit(TH1* hist, const TFitResultPtr& res, 
409                           Double_t low, Double_t high, Int_t& color) const
410 {
411   color = kGreen+4;
412
413   UShort_t   ret   = 0;
414   Int_t      nPar  = res->NPar();
415   Double_t   dy    = .06;
416   Double_t   x     = .2;
417   Double_t   y     = .9-dy;
418   Double_t   chi2  = res->Chi2();
419   Int_t      nu    = res->Ndf();
420   Double_t   red   = (nu == 0 ? fELossFkupChi2Nu : chi2 / nu);
421   TObjArray* lines = 0;
422   TLatex*    lRed  = 0;
423   TLatex*    ltx   = 0;
424   if (fShowFitResults) { 
425     lines = new TObjArray(nPar+3);
426     lines->SetName("lines");
427     lines->SetOwner(true);
428     
429     ltx = new TLatex(x, y, Form("#chi^{2}/#nu: %7.3f",red));
430     ltx->SetNDC(true);
431     ltx->SetTextColor(color);
432     ltx->SetTextSize(dy-.01);
433     lines->Add(ltx);
434     lRed = ltx;
435     
436     Double_t x1 = .85;
437     Double_t y1 = .5;
438     ltx = new TLatex(x1, y1, Form("[thresholds: %6.2f, %6.2f]", 
439                                 fELossBadChi2Nu, fELossFkupChi2Nu));
440     ltx->SetTextColor(kGray+3);
441     ltx->SetTextSize(dy-.01);
442     ltx->SetNDC(true);
443     ltx->SetTextAlign(31);
444     lines->Add(ltx);    
445
446     y1 -= dy;
447     ltx = new TLatex(x1, y1, Form("Fit range: [%6.2f,%6.2f]", low, high));
448     ltx->SetTextColor(kGray+3);
449     ltx->SetTextSize(dy-.01);
450     ltx->SetTextAlign(31);
451     ltx->SetNDC(true);
452     lines->Add(ltx);
453
454     y1 -= dy;
455     ltx = new TLatex(x1, y1, Form("Entries: %d", 
456                                    Int_t(hist->GetEffectiveEntries())));
457     ltx->SetTextColor(kGray+3);
458     ltx->SetTextSize(dy-.01);
459     ltx->SetTextAlign(31);
460     ltx->SetNDC(true);
461     lines->Add(ltx);
462   }
463   
464   if (red > fELossBadChi2Nu) { // || res->Prob() < .01) { 
465     AliWarningF("Fit gave chi^2/nu=%f/%d=%f>%f (%f)", 
466                 res->Chi2(), res->Ndf(), red, fELossBadChi2Nu, 
467                 fELossFkupChi2Nu);
468     if (lRed) lRed->SetTextColor(kOrange+2);
469     res->Print();
470     ret++;
471     if (red > fELossFkupChi2Nu) { 
472       if (lRed) lRed->SetTextColor(kRed+2);
473       ret++;
474     }
475   }
476   // Now check the relative error on the fit parameters 
477   Int_t parsOk = 0;
478   for (Int_t i = 0; i < nPar; i++) { 
479     if (res->IsParameterFixed(i)) continue; 
480     Double_t pv  = res->Parameter(i);
481     Double_t pe  = res->ParError(i);
482     Double_t rel = (pv == 0 ? 100 : pe / pv);
483     if (lines) {
484       y -= dy;
485       ltx = new TLatex(x, y, Form("#delta%s/%s: %7.3f", 
486                                   res->ParName(i).c_str(),
487                                   res->ParName(i).c_str(),
488                                   /*pv, pe,*/ rel));
489       ltx->SetNDC(true);
490       ltx->SetTextColor(color);
491       ltx->SetTextSize(dy-.01);
492       lines->Add(ltx);
493     }
494     if (i == 3) continue; // Skip sigma 
495     Double_t thr = fELossGoodParError;
496     if (rel < thr) parsOk++;
497     else if (ltx) ltx->SetLineColor(kOrange+2);
498   }
499   if (parsOk > 0) 
500     ret = TMath::Max(ret-(parsOk-1),0);
501   if (ret > 1) color = kRed+2;
502   if (ret > 0) color = kOrange+2;
503
504   TList*   lf  = hist->GetListOfFunctions();
505   TObject* old = lf->FindObject(lines->GetName());
506   if (old) {
507     lf->Remove(old);
508     delete old;
509   }
510   lf->Add(lines);
511   hist->SetStats(false);
512     
513   return ret;
514 }
515
516 //__________________________________________________________________
517 void
518 AliFMDQAChecker::AddFitResults(TH1* hist, const TFitResultPtr& res, 
519                                Int_t color, Double_t low, Double_t high) const
520 {
521   if (!fShowFitResults) return; 
522
523   Int_t      nPar  = res->NPar();
524   TObjArray* lines = new TObjArray(nPar+1);
525   lines->SetOwner(kTRUE);
526   lines->SetName("fitResults");
527
528   Double_t   dy    = .06;
529   Double_t   x     = .2;
530   Double_t   y     = .9-dy;
531   Double_t   chi2  = res->Chi2();
532   Int_t      nu    = res->Ndf();
533   Double_t   red   = (nu == 0 ? fELossFkupChi2Nu : chi2 / nu);
534   
535   TLatex* ltx = new TLatex(x, y, Form("#chi^{2}/#nu: %7.3f",red));
536   ltx->SetNDC(true);
537   ltx->SetTextColor(color);
538   ltx->SetTextSize(dy-.01);
539   lines->Add(ltx);
540   
541   y -= dy;
542   ltx = new TLatex(x, y, Form("[thresholds: %6.2f, %6.2f]", 
543                               fELossBadChi2Nu, fELossFkupChi2Nu));
544   ltx->SetTextColor(kGray+3);
545   ltx->SetTextSize(dy-.01);
546   ltx->SetNDC(true);
547   lines->Add(ltx);
548
549   y -= dy;
550   ltx = new TLatex(x, y, Form("Fit range: [%6.2f,%6.2f]", low, high));
551   ltx->SetTextColor(kGray+3);
552   ltx->SetTextSize(dy-.01);
553   ltx->SetNDC(true);
554   lines->Add(ltx);
555
556   for (Int_t i = 0; i < nPar; i++) { 
557     if (res->IsParameterFixed(i)) continue; 
558     y -= dy;
559     Double_t pv  = res->Parameter(i);
560     Double_t pe  = res->ParError(i);
561     Double_t rel = (pv == 0 ? 100 : pe / pv);
562     ltx = new TLatex(x, y, Form("#delta%s/%s: %7.3f", 
563                                 res->ParName(i).c_str(),
564                                 res->ParName(i).c_str(),
565                                 /*pv, pe,*/ rel));
566     ltx->SetNDC(true);
567     ltx->SetTextColor(color);
568     ltx->SetTextSize(dy-.01);
569     lines->Add(ltx);
570   }
571   TList*   lf  = hist->GetListOfFunctions();
572   TObject* old = lf->FindObject(lines->GetName());
573   if (old) {
574     lf->Remove(old);
575     delete old;
576   }
577   lf->Add(lines);
578   hist->SetStats(false);
579 }
580
581 //__________________________________________________________________
582 UShort_t
583 AliFMDQAChecker::CheckRaw(AliRecoParam::EventSpecie_t /* specie*/, 
584                           TH1*                        hist) const
585 {
586   Int_t ret = BasicCheck(hist);
587   TString name(hist->GetName());
588   if (name.Contains("readouterrors", TString::kIgnoreCase)) { 
589     // Check the mean number of errors per event 
590     TH2*  roErrors = static_cast<TH2*>(hist);
591     Int_t nY       = roErrors->GetNbinsY();
592
593     TLatex* ltx = new TLatex(.15, .8, Form("Thresholds: %5.2f,%5.2f",
594                                            fROErrorsBad, fROErrorsFkup));
595     ltx->SetName("thresholds");
596     ltx->SetTextColor(kGray+3);
597     ltx->SetNDC();
598
599     TList* ll = hist->GetListOfFunctions(); 
600     TObject* old = ll->FindObject(ltx->GetName());
601     if (old) { 
602       ll->Remove(old);
603       delete old;
604     }
605     ll->Add(ltx);
606     
607     for (Int_t i = 1; i <= 3; i++) {
608       Double_t sum = 0;
609       Int_t    cnt = 0;
610       for (Int_t j = 1; j <= nY; j++) {
611         Int_t n =  roErrors->GetBinContent(i, j);
612         sum     += n * roErrors->GetYaxis()->GetBinCenter(j);
613         cnt     += n;
614       }
615       Double_t mean = sum / cnt;
616
617       ltx = new TLatex(i, kROErrorsLabelY, Form("Mean: %6.3f", mean));
618       ltx->SetName(Form("FMD%d", i));
619       ltx->SetTextAngle(90);
620       ltx->SetTextColor(kGreen+4);
621       old = ll->FindObject(ltx->GetName());
622       if (old) { 
623         ll->Remove(old);
624         delete old;
625       }
626       ll->Add(ltx);
627
628       if (mean > fROErrorsBad) { 
629         AliWarningF("Mean of readout errors for FMD%d = %f > %f (%f)", 
630                     i, mean, fROErrorsBad, fROErrorsFkup);
631         ret++;
632         ltx->SetTextColor(kOrange+2);
633         if (mean > fROErrorsFkup) {
634           ret++;
635           ltx->SetTextColor(kRed+2);
636         }
637       }
638     }
639   }
640   else if (name.Contains("eloss",TString::kIgnoreCase)) { 
641     // Try to fit a function to the histogram 
642     if (hist->GetEntries() < 1000) return ret;
643
644     Double_t xMin  = hist->GetXaxis()->GetXmin();
645     Double_t xMax  = hist->GetXaxis()->GetXmax();
646
647     hist->GetXaxis()->SetRangeUser(fELossLowCut, xMax);
648     Int_t    bMaxY = hist->GetMaximumBin();
649     Double_t xMaxY = hist->GetXaxis()->GetBinCenter(bMaxY);
650     Double_t rms   = hist->GetRMS();
651     Double_t low   = hist->GetXaxis()->GetBinCenter(bMaxY-4);
652     hist->GetXaxis()->SetRangeUser(0.2, xMaxY+(fELossNRMS+1)*rms);
653     rms  = hist->GetRMS();
654     hist->GetXaxis()->SetRange(0,-1);
655     TF1*          func = makeLandauGaus(name);
656     func->SetParameter(1, xMaxY);
657     func->SetLineColor(kGreen+4);
658     // func->SetLineStyle(2);
659     Double_t high = xMax; // xMaxY+fELossNRMS*rms;
660
661     TFitResultPtr res  = hist->Fit(func, "QS", "", low, high);
662
663     Int_t    color = func->GetLineColor();
664     UShort_t qual  = CheckFit(hist, res, low, high, color);
665
666     // Make sure we save the function in the full range of the histogram
667     func = hist->GetFunction("landauGaus");
668     func->SetRange(xMin, xMax);
669     // func->SetParent(hist);
670     func->Save(xMin, xMax, 0, 0, 0, 0);
671     func->SetLineColor(color);
672     
673     if (qual > 0) { 
674       func->SetLineWidth(3);
675       func->SetLineStyle(1);
676       if (qual > 1) 
677         func->SetLineWidth(4);
678     }
679     ret += qual;
680   }
681
682   return ret;
683 }
684 //__________________________________________________________________
685 UShort_t
686 AliFMDQAChecker::CheckSim(AliRecoParam::EventSpecie_t /* specie*/, 
687                           TH1*                        hist) const
688 {
689   return BasicCheck(hist);
690 }
691 //__________________________________________________________________
692 UShort_t
693 AliFMDQAChecker::CheckRec(AliRecoParam::EventSpecie_t /* specie*/, 
694                           TH1*                        hist) const
695 {
696   return BasicCheck(hist);
697 }
698
699 //__________________________________________________________________
700 void AliFMDQAChecker::Check(Double_t*                   rv, 
701                             AliQAv1::ALITASK_t          what, 
702                             TObjArray**                 list, 
703                             const AliDetectorRecoParam* /*t*/) 
704 {
705   // 
706   // Member function called to do the actual checking
707   //
708   // Parameters: 
709   //    rv   Array of return values. 
710   //    what What to check 
711   //    list Array of arrays of histograms.  There's one arrat for
712   //         each 'specie'
713   //    t    Reconstruction parameters - not used. 
714   //
715   // The bounds defined for RV are 
716   // 
717   //   FATAL:      [-1,   0.0]
718   //   ERROR:      (0.0,  0.002]
719   //   WARNING:    (0.002,0.5]
720   //   INFO:       (0.5,  1.0]
721   //
722   // Double_t* rv = new Double_t[AliRecoParam::kNSpecies] ; 
723
724   if (!fDidExternal) ProcessExternalParams();
725
726   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) {
727     // Int_t count   = 0;
728     rv[specie]    = 0.; 
729
730     if (!AliQAv1::Instance()->IsEventSpecieSet(specie) ) 
731       continue ;
732     
733     if(!list[specie]) continue;
734     
735     TH1* hist  = 0;
736     Int_t nHist = list[specie]->GetEntriesFast();
737     UShort_t ret = 0;
738     for(Int_t i= 0; i< nHist; i++) {
739       if (!(hist = static_cast<TH1*>(list[specie]->At(i)))) continue;
740       Int_t qual = CheckOne(what, AliRecoParam::ConvertIndex(specie), hist);
741       hist->SetUniqueID(Quality2Bit(qual));
742       ret += qual;
743     } // for (int i ...)
744     rv[specie] = ret;
745     // if      (ret > kWhatTheFk) rv[specie] = fLowTestValue[AliQAv1::kFATAL];
746     // else if (ret > kBad)       rv[specie] = fUpTestValue[AliQAv1::kERROR]; 
747     // else if (ret > kProblem)   rv[specie] = fUpTestValue[AliQAv1::kWARNING]; 
748     // else                       rv[specie] = fUpTestValue[AliQAv1::kINFO]; 
749     AliDebugF(3, "Combined sum is %d -> %f", ret, rv[specie]);
750
751     // if (count != 0) rv[specie] /= count;
752   }
753   // return rv;
754 }
755
756 namespace {
757   Int_t CheckForLog(TAxis*       axis,
758                     TVirtualPad* pad, 
759                     Int_t        xyz)
760   {
761     Int_t ret = 0;
762     TString t(axis->GetTitle());
763     if (!t.Contains("[log]", TString::kIgnoreCase)) return 0;
764     t.ReplaceAll("[log]", "");
765     switch (xyz) { 
766     case 1: pad->SetLogx(); ret |= 0x1; break;
767     case 2: pad->SetLogy(); ret |= 0x2; break;
768     case 3: pad->SetLogz(); ret |= 0x4; break;
769     }
770     axis->SetTitle(t);
771     return ret;
772   }
773   void RestoreLog(TAxis* axis, Bool_t log) 
774   {
775     if (!log) return;
776     TString t(axis->GetTitle());
777     t.Append("[log]");
778     axis->SetTitle(t);
779   }
780 }
781
782 namespace {
783   void FindMinMax(TH1* h, Double_t& min, Double_t& max)
784   {
785     Double_t tmin = 1e9;
786     Double_t tmax = 0;
787     for (Int_t i = 1; i <= h->GetNbinsX(); i++) { 
788       Double_t c = h->GetBinContent(i);
789       if (c < 1e-8) continue;
790       tmin = TMath::Min(tmin, c);
791       tmax = TMath::Max(tmax, c);
792     }
793     min = tmin;
794     max = tmax;
795   }
796 }
797 //____________________________________________________________________________ 
798 void 
799 AliFMDQAChecker::MakeImage(TObjArray** list, 
800                            AliQAv1::TASKINDEX_t task, 
801                            AliQAv1::MODE_t mode) 
802 {
803   // makes the QA image for sim and rec
804   // 
805   // Parameters: 
806   //    task What to check 
807   //    list Array of arrays of histograms.  There's one array for
808   //         each 'specie'
809   //    t    Reconstruction parameters - not used. 
810   // 
811   Int_t    nImages = 0 ;
812   Double_t max     = 0;
813   Double_t min     = 10000;
814
815   // Loop over all species 
816   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) {
817     AliRecoParam::EventSpecie_t spe = AliRecoParam::ConvertIndex(specie);
818     if (!AliQAv1::Instance(AliQAv1::GetDetIndex(GetName()))
819         ->IsEventSpecieSet(spe)) 
820       continue;
821                                                                         
822     // If nothing is defined for this specie, go on. 
823     if(!list[specie] || list[specie]->GetEntriesFast() == 0) continue;
824
825     // Loop over the histograms and figure out how many histograms we
826     // have and the min/max 
827     TH1* hist  = 0;
828     Int_t nHist = list[specie]->GetEntriesFast();
829     for(Int_t i= 0; i< nHist; i++) {
830       hist = static_cast<TH1F*>(list[specie]->At(i));
831       if (hist && hist->TestBit(AliQAv1::GetImageBit())) {
832         nImages++; 
833         TString name(hist->GetName());
834         if (name.Contains("readouterrors", TString::kIgnoreCase)) continue;
835
836         // Double_t hMax = hist->GetMaximum(); 
837         // hist->GetBinContent(hist->GetMaximumBin());
838         // Double_t hMin = hist->GetMinimum();
839         // hist->GetBinContent(hist->GetMinimumBin());
840         Double_t hMax, hMin;
841         FindMinMax(hist, hMin, hMax);
842         max = TMath::Max(max, hMax);
843         min = TMath::Min(min, hMin);
844         // AliInfoF("Min/max of %40s: %f/%f, global -> %f/%f", 
845         //          hist->GetName(), hMin, hMax, min, max);
846       }
847     }
848     break ; 
849   }
850   // AliInfoF("Global min/max=%f/%f", min, max);
851   min = TMath::Max(1e-1, min);
852   max = TMath::Max(1e5,  max);
853
854   // IF no images, go on. 
855   if (nImages == 0) {
856     AliDebug(AliQAv1::GetQADebugLevel(), 
857              Form("No histogram will be plotted for %s %s\n", GetName(), 
858                   AliQAv1::GetTaskName(task).Data()));
859     return;
860   }
861
862   AliDebug(AliQAv1::GetQADebugLevel(), 
863            Form("%d histograms will be plotted for %s %s\n", 
864                 nImages, GetName(), AliQAv1::GetTaskName(task).Data()));  
865   gStyle->SetOptStat(0);
866   
867   // Again loop over species and draw a canvas 
868   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) {
869     if (!AliQAv1::Instance(AliQAv1::GetDetIndex(GetName()))
870         ->IsEventSpecieSet(specie)) continue;
871
872     // if Nothing here, go on
873     if(!list[specie] || list[specie]->GetEntries() <= 0 || 
874        nImages <= 0) continue;
875
876    // Form the title 
877     const Char_t * title = Form("QA_%s_%s_%s", GetName(), 
878                                 AliQAv1::GetTaskName(task).Data(), 
879                                 AliRecoParam::GetEventSpecieName(specie)); 
880     if (!fImage[specie]) fImage[specie] = new TCanvas(title, title) ;
881     fImage[specie]->Clear() ; 
882     fImage[specie]->SetTitle(title) ; 
883     fImage[specie]->cd() ; 
884
885     // Put something in the canvas - even if empty 
886     TPaveText someText(0.015, 0.015, 0.98, 0.98) ;
887     someText.AddText(title) ;
888     someText.SetFillColor(0);
889     someText.SetFillStyle(0);
890     someText.SetBorderSize(0);
891     someText.SetTextColor(kRed+1);
892     someText.Draw() ; 
893     TString outName(Form("%s%s%d.%s", AliQAv1::GetImageFileName(), 
894                          AliQAv1::GetModeName(mode), 
895                          AliQAChecker::Instance()->GetRunNumber(), 
896                          AliQAv1::GetImageFileFormat()));
897     fImage[specie]->Print(outName, "ps") ; 
898
899     // Now set some parameters on the canvas 
900     fImage[specie]->Clear(); 
901     fImage[specie]->SetTopMargin(0.10);
902     fImage[specie]->SetBottomMargin(0.15);
903     fImage[specie]->SetLeftMargin(0.15);
904     fImage[specie]->SetRightMargin(0.05);
905     
906     // Put title on top 
907     const char* topT = Form("Mode: %s, Task: %s, Specie: %s, Run: %d",
908                             AliQAv1::GetModeName(mode), 
909                             AliQAv1::GetTaskName(task).Data(), 
910                             AliRecoParam::GetEventSpecieName(specie),
911                             AliQAChecker::Instance()->GetRunNumber());
912     TLatex* topText = new TLatex(.5, .99, topT);
913     topText->SetTextAlign(23);
914     topText->SetTextSize(.038);
915     topText->SetTextFont(42);
916     topText->SetTextColor(kBlue+3);
917     topText->SetNDC();
918     topText->Draw();
919                                  
920     // Divide canvas 
921     Int_t nx = int(nImages + .5) / 2;
922     Int_t ny = 2;
923     // if (fDoScale) 
924     fImage[specie]->Divide(nx, ny, 0, 0);
925     // else fImage[specie]->Divide(nx, ny);
926     
927     
928     // Loop over histograms 
929     TH1*  hist  = 0;
930     Int_t nHist = list[specie]->GetEntriesFast();
931     Int_t j     = 0;
932     for (Int_t i = 0; i < nHist; i++) { 
933       hist = static_cast<TH1*>(list[specie]->At(i));
934       if (!hist || !hist->TestBit(AliQAv1::GetImageBit())) continue;
935
936       // Go to sub-pad 
937       TVirtualPad* pad = fImage[specie]->cd(++j);
938       pad->SetRightMargin(0.01);
939       if (!fDoScale) { 
940         pad->SetLeftMargin(0.10);
941         pad->SetBottomMargin(0.10);
942       }
943
944       // Check for log scale 
945       Int_t logOpts = 0;
946       logOpts |= CheckForLog(hist->GetXaxis(), pad, 1);
947       logOpts |= CheckForLog(hist->GetYaxis(), pad, 2);
948       logOpts |= CheckForLog(hist->GetZaxis(), pad, 3);
949
950       // Figure out special cases 
951       TString opt("");
952       TString name(hist->GetName());
953       if (name.Contains("readouterrors", TString::kIgnoreCase)) {
954         pad->SetRightMargin(0.15);
955         pad->SetBottomMargin(0.10);
956         // pad->SetTopMargin(0.02);
957         opt="COLZ";
958       }
959       else {
960         pad->SetGridx();
961         pad->SetGridy();
962         if (fDoScale) { 
963           hist->SetMinimum(min);
964           hist->SetMaximum(max);
965         }
966         else { 
967           hist->SetMinimum();
968           hist->SetMaximum();
969         }
970       }
971       // Draw (As a copy)
972       hist->DrawCopy(opt);
973       
974       // Special cases 
975       if (name.Contains("readouterrors", TString::kIgnoreCase)) {
976 #if 0
977         for (Int_t kk = 1; kk <= 3; kk++) {
978           TH1* proj = static_cast<TH2*>(hist)->ProjectionY("",kk,kk);
979           Double_t m = proj->GetMean(); 
980           TLatex* l = new TLatex(kk, 30, Form("Mean: %f", m));
981           l->SetTextAngle(90);
982           l->SetTextColor(m > 10 ? kRed+1 : m > .7 ? kOrange+2 :kGreen+2);
983           l->Draw();
984         }
985 #endif
986       }
987       else {
988         gStyle->SetOptTitle(0);
989         TPad* insert = new TPad("insert", "Zoom", 
990                                 .4,.4, .99, .95, 0, 0, 0);
991         insert->SetTopMargin(0.01);
992         insert->SetRightMargin(0.01);
993         insert->SetFillColor(0);
994         insert->SetBorderSize(1);
995         insert->SetBorderMode(0);
996         insert->Draw();
997         insert->cd();
998         if (logOpts & 0x1) insert->SetLogx();
999         if (logOpts & 0x2) insert->SetLogy();
1000         if (logOpts & 0x4) insert->SetLogz();
1001         hist->GetXaxis()->SetRange(1, hist->GetNbinsX()/8);
1002         TH1* copy = hist->DrawCopy(opt);
1003         copy->GetXaxis()->SetNdivisions(408, false);
1004         // Restore full range 
1005         hist->GetXaxis()->SetRange(0, 0);
1006         gStyle->SetOptTitle(1);
1007       }
1008       pad->cd();
1009       // Possibly restore the log options 
1010       RestoreLog(hist->GetXaxis(), logOpts & 0x1);
1011       RestoreLog(hist->GetYaxis(), logOpts & 0x2);
1012       RestoreLog(hist->GetZaxis(), logOpts & 0x4);
1013     }
1014     // Print to a post-script file 
1015     fImage[specie]->Print(outName, "ps");
1016 #if 0
1017     fImage[specie]->Print(Form("%s_%d.png", 
1018                                AliRecoParam::GetEventSpecieName(specie), 
1019                                AliQAChecker::Instance()->GetRunNumber()));
1020 #endif
1021   }
1022 }
1023
1024 //__________________________________________________________________
1025 //
1026 // EOF
1027 //