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