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     fELossMaxEntries(-1),
 258     fELossGoodParError(0.1),
 259     fROErrorsBad(0.3),
 260     fROErrorsFkup(0.5)
 261 {
 262 }
 263
 264 //__________________________________________________________________
 265 void
 266 AliFMDQAChecker::ProcessExternalParams()
 267 {
 268   ProcessExternalParam("ELossLowCut",           fELossLowCut);
 269   ProcessExternalParam("ELossNRMS",             fELossNRMS);
 270   ProcessExternalParam("ELossBadChi2Nu",        fELossBadChi2Nu);
 271   ProcessExternalParam("ELossFkupChi2Nu",       fELossFkupChi2Nu);
 272   ProcessExternalParam("ELossGoodParError",     fELossGoodParError);
 273   ProcessExternalParam("ROErrorsBad",           fROErrorsBad);
 274   ProcessExternalParam("ROErrorsFkup",          fROErrorsFkup);
 275   Double_t tmp = 0;
 276   ProcessExternalParam("CommonScale",           tmp);
 277   fDoScale = tmp > 0;
 278   ProcessExternalParam("ELossMinEntries",       tmp);
 279   fELossMinEntries = tmp;
 280   ProcessExternalParam("ELossMaxEntries",       tmp);
 281   fELossMaxEntries = tmp;
 282
 283   GetThresholds();
 284
 285   fDidExternal = true;
 286 }
 287 //__________________________________________________________________
 288 void
 289 AliFMDQAChecker::ProcessExternalParam(const char* name, Double_t& v)
 290 {
 291   TObject* o = fExternParamList->FindObject(name);
 292   if (!o) return;
 293   TParameter<double>* p = static_cast<TParameter<double>*>(o);
 294   v = p->GetVal();
 295   AliDebugF(3, "External parameter: %-20s=%lf", name, v);
 296 }
 297
 298 //__________________________________________________________________
 299 void
 300 AliFMDQAChecker::GetThresholds()
 301 {
 302   const char*    path   = "GRP/Calib/QAThresholds";
 303   AliCDBManager* cdbMan = AliCDBManager::Instance();
 304   AliCDBEntry*   cdbEnt = cdbMan->Get(path);
 305   if (!cdbEnt) {
 306     AliWarningF("Failed to get CDB entry at %s", path);
 307     return;
 308   }
 309
 310   TObjArray*     cdbObj = static_cast<TObjArray*>(cdbEnt->GetObject());
 311   if (!cdbObj) {
 312     AliWarningF("Failed to get CDB object at %s", path);
 313     return;
 314   }
 315
 316   TObject*       fmdObj = cdbObj->FindObject("FMD");
 317   if (!fmdObj) {
 318     AliWarningF("Failed to get FMD object at from CDB %s", path);
 319     return;
 320   }
 321
 322   AliQAThresholds* qaThr = static_cast<AliQAThresholds*>(fmdObj);
 323   Int_t nThr = qaThr->GetSize();
 324   for (Int_t i = 0; i < nThr; i++) {
 325     TObject* thr = qaThr->GetThreshold(i);
 326     if (!thr) continue;
 327
 328     TParameter<double>* d = dynamic_cast<TParameter<double>*>(thr);
 329     if (!d) {
 330       AliWarningF("Parameter %s not of type double", thr->GetName());
 331       continue;
 332     }
 333     Double_t val = d->GetVal();
 334     TString  name(thr->GetName());
 335     if      (name.EqualTo("ELossBadChi2Nu"))     fELossBadChi2Nu    = val;
 336     else if (name.EqualTo("ELossFkupChi2Nu"))    fELossFkupChi2Nu   = val;
 337     else if (name.EqualTo("ELossGoodParError"))  fELossGoodParError = val;
 338     else if (name.EqualTo("ROErrorsBad"))        fROErrorsBad       = val;
 339     else if (name.EqualTo("ROErrorsFkup"))       fROErrorsFkup      = val;
 340     AliDebugF(3, "Threshold %s=%f", name.Data(), val);
 341   }
 342 }
 343
 344 //__________________________________________________________________
 345 AliQAv1::QABIT_t
 346 AliFMDQAChecker::Quality2Bit(UShort_t value) const
 347 {
 348   AliQAv1::QABIT_t  ret   = AliQAv1::kINFO; // Assume success
 349   if      (value >= kWhatTheFk) ret = AliQAv1::kFATAL;
 350   else if (value >= kBad)       ret = AliQAv1::kERROR;
 351   else if (value >= kProblem)   ret = AliQAv1::kWARNING;
 352
 353   return ret;
 354 }
 355
 356 //__________________________________________________________________
 357 void
 358 AliFMDQAChecker::SetQA(AliQAv1::ALITASK_t index, Double_t* values) const
 359 {
 360   AliQAv1 * qa = AliQAv1::Instance(index) ;
 361
 362   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) {
 363     // Check if specie is defined
 364     if (!qa->IsEventSpecieSet(AliRecoParam::ConvertIndex(specie)))
 365       continue ;
 366
 367     // No checker is implemented, set all QA to Fatal
 368     if (!values) {
 369       qa->Set(AliQAv1::kFATAL, specie) ;
 370       continue;
 371     }
 372
 373     UShort_t          value = values[specie];
 374     AliQAv1::QABIT_t  ret   = Quality2Bit(value);
 375
 376     qa->Set(ret, AliRecoParam::ConvertIndex(specie));
 377     AliDebugF(3, "Quality of %s: %d -> %d",
 378               AliRecoParam::GetEventSpecieName(specie), value, ret);
 379   }
 380 }
 381
 382 //__________________________________________________________________
 383 UShort_t
 384 AliFMDQAChecker::BasicCheck(TH1* hist) const
 385 {
 386   if (hist->GetEntries() <= 0) return kOK;
 387   return (hist->GetMean() > 0 ? kOK : kProblem);
 388 }
 389
 390 //__________________________________________________________________
 391 UShort_t
 392 AliFMDQAChecker::CheckOne(AliQAv1::ALITASK_t          what,
 393                           AliRecoParam::EventSpecie_t specie,
 394                           TH1*                        hist) const
 395 {
 396   if(what == AliQAv1::kESD) return CheckESD(specie, hist);
 397   if(what == AliQAv1::kRAW) return CheckRaw(specie, hist);
 398   if(what == AliQAv1::kSIM) return CheckSim(specie, hist);
 399   if(what == AliQAv1::kREC) return CheckRec(specie, hist);
 400   return 0;
 401 }
 402 //__________________________________________________________________
 403 UShort_t
 404 AliFMDQAChecker::CheckESD(AliRecoParam::EventSpecie_t /* specie*/,
 405                           TH1*                        hist) const
 406 {
 407   return BasicCheck(hist);
 408 }
 409 //__________________________________________________________________
 410 UShort_t
 411 AliFMDQAChecker::CheckFit(TH1* hist, const TFitResultPtr& res,
 412                           Double_t low, Double_t high, Int_t& color) const
 413 {
 414   color = kGreen+4;
 415
 416   UShort_t   ret   = 0;
 417   Int_t      nPar  = res->NPar();
 418   Double_t   dy    = .06;
 419   Double_t   x     = .2;
 420   Double_t   y     = .9-dy;
 421   Double_t   chi2  = res->Chi2();
 422   Int_t      nu    = res->Ndf();
 423   Double_t   red   = (nu == 0 ? fELossFkupChi2Nu : chi2 / nu);
 424   TObjArray* lines = 0;
 425   // TLatex*    lRed  = 0;
 426   TLatex*    ltx   = 0;
 427   Int_t      chi2Check = 0;
 428   Double_t   chi2Lim   = fELossBadChi2Nu;
 429   if (red > fELossBadChi2Nu) { // || res->Prob() < .01) {
 430     // AliWarningF("Fit gave chi^2/nu=%f/%d=%f>%f (%f)",
 431     //             res->Chi2(), res->Ndf(), red, fELossBadChi2Nu,
 432     //             fELossFkupChi2Nu);
 433     // res->Print();
 434     chi2Check++;
 435     if (red > fELossFkupChi2Nu) {
 436       chi2Check++;
 437       chi2Lim = fELossFkupChi2Nu;
 438     }
 439   }
 440   ret += chi2Check;
 441
 442   if (fShowFitResults) {
 443     lines = new TObjArray(nPar+3);
 444     lines->SetName("lines");
 445     lines->SetOwner(true);
 446
 447     ltx = new TLatex(x, y, Form("#chi^{2}/#nu: %7.3f %c %6.2f",
 448                                 red, chi2Check < 1 ? '<' : '>',
 449                                 chi2Lim));
 450     ltx->SetNDC(true);
 451     ltx->SetTextColor(chi2Check < 1 ? color :
 452                       chi2Check < 2 ? kOrange+2 : kRed+2);
 453     // ltx->SetTextColor(color);
 454     ltx->SetTextSize(dy-.01);
 455     lines->Add(ltx);
 456     // lRed = ltx;
 457
 458     Double_t x1 = .85;
 459     Double_t y1 = .5;
 460 #if 0
 461     ltx = new TLatex(x1, y1, Form("[thresholds: %6.2f, %6.2f]",
 462                                 fELossBadChi2Nu, fELossFkupChi2Nu));
 463     ltx->SetTextColor(kGray+3);
 464     ltx->SetTextSize(dy-.01);
 465     ltx->SetNDC(true);
 466     ltx->SetTextAlign(31);
 467     // lines->Add(ltx);
 468 #endif
 469
 470     // y1 -= dy;
 471     ltx = new TLatex(x1, y1, Form("Fit range: [%6.2f,%6.2f]", low, high));
 472     ltx->SetTextColor(kGray+3);
 473     ltx->SetTextSize(dy-.01);
 474     ltx->SetTextAlign(31);
 475     ltx->SetNDC(true);
 476     lines->Add(ltx);
 477
 478     y1 -= dy;
 479     ltx = new TLatex(x1, y1, Form("Entries: %d (%d)",
 480                                   Int_t(hist->GetEffectiveEntries()),
 481                                   fELossMaxEntries));
 482     ltx->SetTextColor(kGray+3);
 483     ltx->SetTextSize(dy-.01);
 484     ltx->SetTextAlign(31);
 485     ltx->SetNDC(true);
 486     lines->Add(ltx);
 487   }
 488
 489   // Now check the relative error on the fit parameters
 490   Int_t parsOk = 0;
 491   for (Int_t i = 0; i < nPar; i++) {
 492     if (res->IsParameterFixed(i)) continue;
 493     Double_t thr = fELossGoodParError;
 494     Double_t pv  = res->Parameter(i);
 495     Double_t pe  = res->ParError(i);
 496     Double_t rel = (pv == 0 ? 100 : pe / pv);
 497     Bool_t   ok  = (i == 3) || (rel < thr);
 498     if (lines) {
 499       y -= dy;
 500       TString txt(Form("#delta%s/%s: %7.3f ",
 501                        res->ParName(i).c_str(),
 502                        res->ParName(i).c_str(),
 503                        /*pv, pe,*/ rel));
 504       if (i != 3) txt.Append(Form("%c %4.2f", ok ? '<' : '>', thr));
 505       else        txt.Append("(ignored)");
 506       ltx = new TLatex(x, y, txt);
 507       ltx->SetNDC(true);
 508       ltx->SetTextSize(dy-.01);
 509       ltx->SetTextColor(ok ? color : kOrange+2);
 510       lines->Add(ltx);
 511     }
 512     if (i == 3) continue; // Skip sigma
 513     if (ok) parsOk++;
 514   }
 515   if (parsOk > 0)
 516     ret = TMath::Max(ret-(parsOk-1),0);
 517   if (ret > 1) color = kRed+2;
 518   if (ret > 0) color = kOrange+2;
 519
 520   TList*   lf  = hist->GetListOfFunctions();
 521   TObject* old = lf->FindObject(lines->GetName());
 522   if (old) {
 523     lf->Remove(old);
 524     delete old;
 525   }
 526   lf->Add(lines);
 527   hist->SetStats(false);
 528
 529   return ret;
 530 }
 531
 532 //__________________________________________________________________
 533 void
 534 AliFMDQAChecker::AddFitResults(TH1* hist, const TFitResultPtr& res,
 535                                Int_t color, Double_t low, Double_t high) const
 536 {
 537   // Obsolete - not used
 538   if (!fShowFitResults) return;
 539
 540   Int_t      nPar  = res->NPar();
 541   TObjArray* lines = new TObjArray(nPar+1);
 542   lines->SetOwner(kTRUE);
 543   lines->SetName("fitResults");
 544
 545   Double_t   dy    = .06;
 546   Double_t   x     = .2;
 547   Double_t   y     = .9-dy;
 548   Double_t   chi2  = res->Chi2();
 549   Int_t      nu    = res->Ndf();
 550   Double_t   red   = (nu == 0 ? fELossFkupChi2Nu : chi2 / nu);
 551
 552   TLatex* ltx = new TLatex(x, y, Form("#chi^{2}/#nu: %7.3f",red));
 553   ltx->SetNDC(true);
 554   ltx->SetTextColor(color);
 555   ltx->SetTextSize(dy-.01);
 556   lines->Add(ltx);
 557
 558   y -= dy;
 559   ltx = new TLatex(x, y, Form("[thresholds: %6.2f, %6.2f]",
 560                               fELossBadChi2Nu, fELossFkupChi2Nu));
 561   ltx->SetTextColor(kGray+3);
 562   ltx->SetTextSize(dy-.01);
 563   ltx->SetNDC(true);
 564   lines->Add(ltx);
 565
 566   y -= dy;
 567   ltx = new TLatex(x, y, Form("Fit range: [%6.2f,%6.2f]", low, high));
 568   ltx->SetTextColor(kGray+3);
 569   ltx->SetTextSize(dy-.01);
 570   ltx->SetNDC(true);
 571   lines->Add(ltx);
 572
 573   for (Int_t i = 0; i < nPar; i++) {
 574     if (res->IsParameterFixed(i)) continue;
 575     y -= dy;
 576     Double_t pv  = res->Parameter(i);
 577     Double_t pe  = res->ParError(i);
 578     Double_t rel = (pv == 0 ? 100 : pe / pv);
 579     ltx = new TLatex(x, y, Form("#delta%s/%s: %7.3f",
 580                                 res->ParName(i).c_str(),
 581                                 res->ParName(i).c_str(),
 582                                 /*pv, pe,*/ rel));
 583     ltx->SetNDC(true);
 584     ltx->SetTextColor(color);
 585     ltx->SetTextSize(dy-.01);
 586     lines->Add(ltx);
 587   }
 588   TList*   lf  = hist->GetListOfFunctions();
 589   TObject* old = lf->FindObject(lines->GetName());
 590   if (old) {
 591     lf->Remove(old);
 592     delete old;
 593   }
 594   lf->Add(lines);
 595   hist->SetStats(false);
 596 }
 597
 598 //__________________________________________________________________
 599 UShort_t
 600 AliFMDQAChecker::CheckRaw(AliRecoParam::EventSpecie_t specie,
 601                           TH1*                        hist) const
 602 {
 603   Int_t ret = BasicCheck(hist);
 604   TString name(hist->GetName());
 605   if (name.Contains("readouterrors", TString::kIgnoreCase)) {
 606     // Check the mean number of errors per event
 607     TH2*  roErrors = static_cast<TH2*>(hist);
 608     Int_t nY       = roErrors->GetNbinsY();
 609
 610     TLatex* ltx = new TLatex(.15, .9, Form("Thresholds: %5.2f,%5.2f",
 611                                            fROErrorsBad, fROErrorsFkup));
 612     ltx->SetName("thresholds");
 613     ltx->SetTextColor(kGray+3);
 614     ltx->SetNDC();
 615
 616     TList* ll = hist->GetListOfFunctions();
 617     TObject* old = ll->FindObject(ltx->GetName());
 618     if (old) {
 619       ll->Remove(old);
 620       delete old;
 621     }
 622     ll->Add(ltx);
 623
 624     for (Int_t i = 1; i <= 3; i++) {
 625       Double_t sum = 0;
 626       Int_t    cnt = 0;
 627       for (Int_t j = 1; j <= nY; j++) {
 628         Int_t n =  roErrors->GetBinContent(i, j);
 629         sum     += n * roErrors->GetYaxis()->GetBinCenter(j);
 630         cnt     += n;
 631       }
 632       Double_t mean = sum / cnt;
 633       Double_t x    = ((i-.5) * (1-0.1-0.1) / 3 + 0.1);
 634
 635       ltx = new TLatex(x, kROErrorsLabelY, Form("Mean: %6.3f", mean));
 636       ltx->SetName(Form("FMD%d", i));
 637       ltx->SetNDC();
 638       ltx->SetTextAngle(90);
 639       ltx->SetTextColor(kGreen+4);
 640       old = ll->FindObject(ltx->GetName());
 641       if (old) {
 642         ll->Remove(old);
 643         delete old;
 644       }
 645       ll->Add(ltx);
 646
 647       if (mean > fROErrorsBad) {
 648         AliWarningF("Mean of readout errors for FMD%d = %f > %f (%f)",
 649                     i, mean, fROErrorsBad, fROErrorsFkup);
 650         ret++;
 651         ltx->SetTextColor(kOrange+2);
 652         if (mean > fROErrorsFkup) {
 653           ret++;
 654           ltx->SetTextColor(kRed+2);
 655         }
 656       }
 657     }
 658   }
 659   else if (name.Contains("eloss",TString::kIgnoreCase)) {
 660     // Try to fit a function to the histogram
 661     if (hist->GetEntries() < 1000) return ret;
 662     if (specie == AliRecoParam::kCosmic ||
 663         specie == AliRecoParam::kCalib) return ret;
 664
 665     Double_t xMin  = hist->GetXaxis()->GetXmin();
 666     Double_t xMax  = hist->GetXaxis()->GetXmax();
 667
 668     hist->GetXaxis()->SetRangeUser(fELossLowCut, xMax);
 669     Int_t    bMaxY = hist->GetMaximumBin();
 670     Double_t xMaxY = hist->GetXaxis()->GetBinCenter(bMaxY);
 671     Double_t rms   = hist->GetRMS();
 672     Double_t low   = hist->GetXaxis()->GetBinCenter(bMaxY-4);
 673     hist->GetXaxis()->SetRangeUser(0.2, xMaxY+(fELossNRMS+1)*rms);
 674     rms  = hist->GetRMS();
 675     hist->GetXaxis()->SetRange(0,-1);
 676     TF1*          func = makeLandauGaus(name);
 677     func->SetParameter(1, xMaxY);
 678     func->SetLineColor(kGreen+4);
 679     // func->SetLineStyle(2);
 680     Double_t high = xMax; // xMaxY+fELossNRMS*rms;
 681
 682     TFitResultPtr res  = hist->Fit(func, "QS", "", low, high);
 683
 684     Int_t    color = func->GetLineColor();
 685     UShort_t qual  = CheckFit(hist, res, low, high, color);
 686
 687     // Make sure we save the function in the full range of the histogram
 688     func = hist->GetFunction("landauGaus");
 689     func->SetRange(xMin, xMax);
 690     // func->SetParent(hist);
 691     func->Save(xMin, xMax, 0, 0, 0, 0);
 692     func->SetLineColor(color);
 693
 694     // Now check if this histogram should be cleared or not
 695     if (fELossMaxEntries > 0 && hist->GetEntries() > fELossMaxEntries)
 696       hist->SetBit(BIT(23));
 697     if (qual > 0) {
 698       func->SetLineWidth(3);
 699       func->SetLineStyle(1);
 700       if (qual > 1)
 701         func->SetLineWidth(4);
 702     }
 703     ret += qual;
 704   }
 705
 706   return ret;
 707 }
 708 //__________________________________________________________________
 709 UShort_t
 710 AliFMDQAChecker::CheckSim(AliRecoParam::EventSpecie_t /* specie*/,
 711                           TH1*                        hist) const
 712 {
 713   return BasicCheck(hist);
 714 }
 715 //__________________________________________________________________
 716 UShort_t
 717 AliFMDQAChecker::CheckRec(AliRecoParam::EventSpecie_t /* specie*/,
 718                           TH1*                        hist) const
 719 {
 720   return BasicCheck(hist);
 721 }
 722
 723 //__________________________________________________________________
 724 void AliFMDQAChecker::Check(Double_t*                   rv,
 725                             AliQAv1::ALITASK_t          what,
 726                             TObjArray**                 list,
 727                             const AliDetectorRecoParam* /*t*/)
 728 {
 729   //
 730   // Member function called to do the actual checking
 731   //
 732   // Parameters:
 733   //    rv   Array of return values.
 734   //    what What to check
 735   //    list Array of arrays of histograms.  There's one arrat for
 736   //         each 'specie'
 737   //    t    Reconstruction parameters - not used.
 738   //
 739   // The bounds defined for RV are
 740   //
 741   //   FATAL:      [-1,   0.0]
 742   //   ERROR:      (0.0,  0.002]
 743   //   WARNING:    (0.002,0.5]
 744   //   INFO:       (0.5,  1.0]
 745   //
 746   // Double_t* rv = new Double_t[AliRecoParam::kNSpecies] ;
 747
 748   if (!fDidExternal) ProcessExternalParams();
 749
 750   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) {
 751     // Int_t count   = 0;
 752     rv[specie]    = 0.;
 753
 754     if (!AliQAv1::Instance()->IsEventSpecieSet(specie) )
 755       continue ;
 756
 757     if(!list[specie]) continue;
 758
 759     TH1* hist  = 0;
 760     Int_t nHist = list[specie]->GetEntriesFast();
 761     UShort_t ret = 0;
 762     for(Int_t i= 0; i< nHist; i++) {
 763       if (!(hist = static_cast<TH1*>(list[specie]->At(i)))) continue;
 764       Int_t qual = CheckOne(what, AliRecoParam::ConvertIndex(specie), hist);
 765       hist->SetUniqueID(Quality2Bit(qual));
 766       ret += qual;
 767     } // for (int i ...)
 768     rv[specie] = ret;
 769     // if      (ret > kWhatTheFk) rv[specie] = fLowTestValue[AliQAv1::kFATAL];
 770     // else if (ret > kBad)       rv[specie] = fUpTestValue[AliQAv1::kERROR];
 771     // else if (ret > kProblem)   rv[specie] = fUpTestValue[AliQAv1::kWARNING];
 772     // else                       rv[specie] = fUpTestValue[AliQAv1::kINFO];
 773     AliDebugF(3, "Combined sum is %d -> %f", ret, rv[specie]);
 774
 775     // if (count != 0) rv[specie] /= count;
 776   }
 777   // return rv;
 778 }
 779
 780 namespace {
 781   Int_t CheckForLog(TAxis*       axis,
 782                     TVirtualPad* pad,
 783                     Int_t        xyz)
 784   {
 785     Int_t ret = 0;
 786     TString t(axis->GetTitle());
 787     if (!t.Contains("[log]", TString::kIgnoreCase)) return 0;
 788     t.ReplaceAll("[log]", "");
 789     switch (xyz) {
 790     case 1: pad->SetLogx(); ret |= 0x1; break;
 791     case 2: pad->SetLogy(); ret |= 0x2; break;
 792     case 3: pad->SetLogz(); ret |= 0x4; break;
 793     }
 794     axis->SetTitle(t);
 795     return ret;
 796   }
 797   void RestoreLog(TAxis* axis, Bool_t log)
 798   {
 799     if (!log) return;
 800     TString t(axis->GetTitle());
 801     t.Append("[log]");
 802     axis->SetTitle(t);
 803   }
 804 }
 805
 806 namespace {
 807   void FindMinMax(TH1* h, Double_t& min, Double_t& max)
 808   {
 809     Double_t tmin = 1e9;
 810     Double_t tmax = 0;
 811     for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
 812       Double_t c = h->GetBinContent(i);
 813       if (c < 1e-8) continue;
 814       tmin = TMath::Min(tmin, c);
 815       tmax = TMath::Max(tmax, c);
 816     }
 817     min = tmin;
 818     max = tmax;
 819   }
 820 }
 821 //____________________________________________________________________________
 822 void
 823 AliFMDQAChecker::MakeImage(TObjArray** list,
 824                            AliQAv1::TASKINDEX_t task,
 825                            AliQAv1::MODE_t mode)
 826 {
 827   // makes the QA image for sim and rec
 828   //
 829   // Parameters:
 830   //    task What to check
 831   //    list Array of arrays of histograms.  There's one array for
 832   //         each 'specie'
 833   //    t    Reconstruction parameters - not used.
 834   //
 835   Int_t    nImages = 0 ;
 836   Double_t max     = 0;
 837   Double_t min     = 10000;
 838
 839   // Loop over all species
 840   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) {
 841     AliRecoParam::EventSpecie_t spe = AliRecoParam::ConvertIndex(specie);
 842     if (!AliQAv1::Instance(AliQAv1::GetDetIndex(GetName()))
 843         ->IsEventSpecieSet(spe))
 844       continue;
 845
 846     // If nothing is defined for this specie, go on.
 847     if(!list[specie] || list[specie]->GetEntriesFast() == 0) continue;
 848
 849     // Loop over the histograms and figure out how many histograms we
 850     // have and the min/max
 851     TH1* hist  = 0;
 852     Int_t nHist = list[specie]->GetEntriesFast();
 853     for(Int_t i= 0; i< nHist; i++) {
 854       hist = static_cast<TH1F*>(list[specie]->At(i));
 855       if (hist && hist->TestBit(AliQAv1::GetImageBit())) {
 856         nImages++;
 857         TString name(hist->GetName());
 858         if (name.Contains("readouterrors", TString::kIgnoreCase)) continue;
 859
 860         // Double_t hMax = hist->GetMaximum();
 861         // hist->GetBinContent(hist->GetMaximumBin());
 862         // Double_t hMin = hist->GetMinimum();
 863         // hist->GetBinContent(hist->GetMinimumBin());
 864         Double_t hMax, hMin;
 865         FindMinMax(hist, hMin, hMax);
 866         max = TMath::Max(max, hMax);
 867         min = TMath::Min(min, hMin);
 868         // AliInfoF("Min/max of %40s: %f/%f, global -> %f/%f",
 869         //          hist->GetName(), hMin, hMax, min, max);
 870       }
 871     }
 872     break ;
 873   }
 874   // AliInfoF("Global min/max=%f/%f", min, max);
 875   min = TMath::Max(1e-1, min);
 876   max = TMath::Max(1e5,  max);
 877
 878   // IF no images, go on.
 879   if (nImages == 0) {
 880     AliDebug(AliQAv1::GetQADebugLevel(),
 881              Form("No histogram will be plotted for %s %s\n", GetName(),
 882                   AliQAv1::GetTaskName(task).Data()));
 883     return;
 884   }
 885
 886   AliDebug(AliQAv1::GetQADebugLevel(),
 887            Form("%d histograms will be plotted for %s %s\n",
 888                 nImages, GetName(), AliQAv1::GetTaskName(task).Data()));
 889   gStyle->SetOptStat(0);
 890
 891   // Again loop over species and draw a canvas
 892   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) {
 893     if (!AliQAv1::Instance(AliQAv1::GetDetIndex(GetName()))
 894         ->IsEventSpecieSet(specie)) continue;
 895
 896     // if Nothing here, go on
 897     if(!list[specie] || list[specie]->GetEntries() <= 0 ||
 898        nImages <= 0) continue;
 899
 900    // Form the title
 901     const Char_t * title = Form("QA_%s_%s_%s", GetName(),
 902                                 AliQAv1::GetTaskName(task).Data(),
 903                                 AliRecoParam::GetEventSpecieName(specie));
 904     if (!fImage[specie]) fImage[specie] = new TCanvas(title, title) ;
 905     fImage[specie]->Clear() ;
 906     fImage[specie]->SetTitle(title) ;
 907     fImage[specie]->cd() ;
 908
 909     // Put something in the canvas - even if empty
 910     TPaveText someText(0.015, 0.015, 0.98, 0.98) ;
 911     someText.AddText(title) ;
 912     someText.SetFillColor(0);
 913     someText.SetFillStyle(0);
 914     someText.SetBorderSize(0);
 915     someText.SetTextColor(kRed+1);
 916     someText.Draw() ;
 917     TString outName(Form("%s%s%d.%s", AliQAv1::GetImageFileName(),
 918                          AliQAv1::GetModeName(mode),
 919                          AliQAChecker::Instance()->GetRunNumber(),
 920                          AliQAv1::GetImageFileFormat()));
 921     fImage[specie]->Print(outName, "ps") ;
 922
 923     // Now set some parameters on the canvas
 924     fImage[specie]->Clear();
 925     fImage[specie]->SetTopMargin(0.10);
 926     fImage[specie]->SetBottomMargin(0.15);
 927     fImage[specie]->SetLeftMargin(0.15);
 928     fImage[specie]->SetRightMargin(0.05);
 929
 930     // Put title on top
 931     const char* topT = Form("Mode: %s, Task: %s, Specie: %s, Run: %d",
 932                             AliQAv1::GetModeName(mode),
 933                             AliQAv1::GetTaskName(task).Data(),
 934                             AliRecoParam::GetEventSpecieName(specie),
 935                             AliQAChecker::Instance()->GetRunNumber());
 936     TLatex* topText = new TLatex(.5, .99, topT);
 937     topText->SetTextAlign(23);
 938     topText->SetTextSize(.038);
 939     topText->SetTextFont(42);
 940     topText->SetTextColor(kBlue+3);
 941     topText->SetNDC();
 942     topText->Draw();
 943
 944     // Divide canvas
 945     Int_t nx = int(nImages + .5) / 2;
 946     Int_t ny = 2;
 947     // if (fDoScale)
 948     fImage[specie]->Divide(nx, ny, 0, 0);
 949     // else fImage[specie]->Divide(nx, ny);
 950
 951
 952     // Loop over histograms
 953     TH1*  hist  = 0;
 954     Int_t nHist = list[specie]->GetEntriesFast();
 955     Int_t j     = 0;
 956     for (Int_t i = 0; i < nHist; i++) {
 957       hist = static_cast<TH1*>(list[specie]->At(i));
 958       if (!hist || !hist->TestBit(AliQAv1::GetImageBit())) continue;
 959
 960       // Go to sub-pad
 961       TVirtualPad* pad = fImage[specie]->cd(++j);
 962       pad->SetRightMargin(0.01);
 963       if (!fDoScale) {
 964         pad->SetLeftMargin(0.10);
 965         pad->SetBottomMargin(0.10);
 966       }
 967
 968       // Check for log scale
 969       Int_t logOpts = 0;
 970       logOpts |= CheckForLog(hist->GetXaxis(), pad, 1);
 971       logOpts |= CheckForLog(hist->GetYaxis(), pad, 2);
 972       logOpts |= CheckForLog(hist->GetZaxis(), pad, 3);
 973
 974       // Figure out special cases
 975       TString opt("");
 976       TString name(hist->GetName());
 977       if (name.Contains("readouterrors", TString::kIgnoreCase)) {
 978         pad->SetRightMargin(0.15);
 979         pad->SetBottomMargin(0.10);
 980         // pad->SetTopMargin(0.02);
 981         opt="COLZ";
 982       }
 983       else {
 984         pad->SetGridx();
 985         pad->SetGridy();
 986         if (fDoScale) {
 987           hist->SetMinimum(min);
 988           hist->SetMaximum(max);
 989         }
 990         else {
 991           hist->SetMinimum();
 992           hist->SetMaximum();
 993         }
 994       }
 995       // Draw (As a copy)
 996       hist->DrawCopy(opt);
 997
 998       // Special cases
 999       if (name.Contains("readouterrors", TString::kIgnoreCase)) {
1000 #if 0
1001         for (Int_t kk = 1; kk <= 3; kk++) {
1002           TH1* proj = static_cast<TH2*>(hist)->ProjectionY("",kk,kk);
1003           Double_t m = proj->GetMean();
1004           TLatex* l = new TLatex(kk, 30, Form("Mean: %f", m));
1005           l->SetTextAngle(90);
1006           l->SetTextColor(m > 10 ? kRed+1 : m > .7 ? kOrange+2 :kGreen+2);
1007           l->Draw();
1008         }
1009 #endif
1010       }
1011       else {
1012         gStyle->SetOptTitle(0);
1013         TPad* insert = new TPad("insert", "Zoom",
1014                                 .4,.4, .99, .95, 0, 0, 0);
1015         insert->SetTopMargin(0.01);
1016         insert->SetRightMargin(0.01);
1017         insert->SetFillColor(0);
1018         insert->SetBorderSize(1);
1019         insert->SetBorderMode(0);
1020         insert->Draw();
1021         insert->cd();
1022         if (logOpts & 0x1) insert->SetLogx();
1023         if (logOpts & 0x2) insert->SetLogy();
1024         if (logOpts & 0x4) insert->SetLogz();
1025         hist->GetXaxis()->SetRange(1, hist->GetNbinsX()/8);
1026         TH1* copy = hist->DrawCopy(opt);
1027         copy->GetXaxis()->SetNdivisions(408, false);
1028         // Restore full range
1029         hist->GetXaxis()->SetRange(0, 0);
1030         gStyle->SetOptTitle(1);
1031       }
1032       pad->cd();
1033       // Possibly restore the log options
1034       RestoreLog(hist->GetXaxis(), logOpts & 0x1);
1035       RestoreLog(hist->GetYaxis(), logOpts & 0x2);
1036       RestoreLog(hist->GetZaxis(), logOpts & 0x4);
1037     }
1038     // Print to a post-script file
1039     fImage[specie]->Print(outName, "ps");
1040 #if 0
1041     fImage[specie]->Print(Form("%s_%d.png",
1042                                AliRecoParam::GetEventSpecieName(specie),
1043                                AliQAChecker::Instance()->GetRunNumber()));
1044 #endif
1045   }
1046 }
1047
1048 //__________________________________________________________________
1049 //
1050 // EOF
1051 //