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