FMD/AliFMDQAChecker.cxx

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