Updates for better DQM.
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 17 Jul 2012 12:08:27 +0000 (12:08 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 17 Jul 2012 12:08:27 +0000 (12:08 +0000)
Now, we use the reconstructor
to properly pedestal subtract and gain calibrate the
signals.   This allows us - since the signals are now on
the same footing - to fit a function the summed spectra.
The fit function used is a Landau convolved with a Gaussian.
The quality is set based on chi^2/nu of the fits, as well
as the relative errors of the fit parameters.  Both pieces
of information is shown on the plots by adding TLatex
objects to the TH1::fListOfFunctions.   The color of the
function line reflects the quality (green: good, orange:
problem, red: bad).

Also, the average number of read-out errors per sub-detector
FMD1, 2, and 3 are computed and the quality set based on this
number.  The means are shown on the plot together with the
used thresholds.  Again, the color of the text reflects the
quality.

In all cases, the quality of the individual plots are set on
the objects themselves in the member TObject::fUniqueID.  This
can then be retrieved by a Amore agent script to set the quality
of the individual objects in the Amore environment.   This script
is simple and looks like

  // -*- mode: C++ -*-
  #include <map>

  void
  FMDshifterQuality(TObjArray* objs, map<string,int>* map)
  {
    for (Int_t i = 0; i < objs->GetEntriesFast(); i++) {
      TObject* o = objs->At(i);
      if (!o) continue;
      // Printf("Setting quality of %s to %d", o->GetName(), o->GetUniqueID());
      std::string s(o->GetName());
      (*map)[s] = o->GetUniqueID();
    }
  }

This should be added to the AMORE database at P2, and the fmd.configfile
should look like

  # QA species definition. Options: [default] [calib] [low] [high] [cosmic]
  qa_species default calib low high cosmic

  # Low cut for fitting energy loss specrta
   ELossLowCut     0.3
  # Number of RMS to fit beyound maximum
  ELossNRMS       2.
  # Cut on chi^2/nu of fit for bad value
  ELossBadChi2Nu  100.
  # Cut on chi^2/nu of fit for f**ked up value
  ELossFkupChi2Nu 200.
  # Least number of entries before fitting
  ELossMinEntries 10000.
  # Cut on average number of R/O errors per detector
  ROErrorsBad     0.31
  # Cut on average number of R/O errors per detector
  ROErrorsFkup    0.51

  # qaTresholdsDownoad 0

  # _macroDownload 1
  _macroExec1 FMDshifterQuality.C
  #
  # EOF
  #

The AliFMDQAChecker now also properly uses the QA thresholds from
the P2 environment by querying the (specific) CDB storage.  To set
this up, a QAThresholds.configfile should be added for the FMD detector
with the content:

  #
  # Thresholds for the FMD
  #

  # Cut on chi^2/nu of fit for bad value
  ELossBadChi2Nu D 100.
  # Cut on chi^2/nu of fit for f**ked up value
  ELossFkupChi2Nu D 200
  # Cut on parameter relative error
  ELossGoodParError D 0.1
  # Cut on average number of R/O errors per detector
  ROErrorsBad D 0.31
  # Cut on average number of R/O errors per detector
  ROErrorsFkup D 0.51

  #
  # EOF
  #

The revision should eventually be ported to the release, but
I'd like to do some more checks when I'm back from vacation.

Christian

FMD/AliFMDQAChecker.cxx
FMD/AliFMDQAChecker.h
FMD/AliFMDQADataMakerRec.cxx
FMD/AliFMDQADataMakerRec.h
FMD/AliFMDReconstructor.cxx
FMD/AliFMDReconstructor.h

index 4453d11..5a6b7c2 100644 (file)
@@ -27,6 +27,7 @@
 // --- ROOT system ---
 #include <TClass.h>
 #include <TH1F.h> 
+#include <TF1.h> 
 #include <TH2.h> 
 #include <TH1I.h> 
 #include <TIterator.h> 
@@ -37,6 +38,9 @@
 #include <TPaveText.h>
 #include <TStyle.h>
 #include <TLatex.h>
+#include <TFitResult.h>
+#include <TParameter.h>
+#include <TMacro.h>
 
 // --- AliRoot header files ---
 #include "AliLog.h"
 #include "AliQAChecker.h"
 #include "AliFMDQAChecker.h"
 #include "AliRecoParam.h"
+#include <AliCDBManager.h>
+#include <AliCDBEntry.h>
+#include <AliCDBId.h>
+#include <AliQAThresholds.h>
 
 ClassImp(AliFMDQAChecker)
 #if 0
 ; // This is for Emacs! - do not delete
 #endif
 
+namespace {
+  void addFitsMacro(TList* l) {
+    TMacro* m = new TMacro("fits");
+    m->AddLine("void fits() {");
+    m->AddLine("  if (!gPad) { Printf(\"No gPad\"); return; }");
+    m->AddLine("  TList* lp = gPad->GetListOfPrimitives();");
+    m->AddLine("  if (!lp) return;");
+    m->AddLine("  TObject* po = 0;");
+    m->AddLine("  TIter next(lp);");
+    m->AddLine("  while ((po = next())) {");
+    m->AddLine("    if (!po->IsA()->InheritsFrom(TH1::Class())) continue;");
+    m->AddLine("    TH1*   htmp = dynamic_cast<TH1*>(po);");
+    m->AddLine("    TList* lf   = htmp->GetListOfFunctions();");
+    m->AddLine("    TObject* pso = (lf ? lf->FindObject(\"stats\") : 0);");
+    m->AddLine("    if (!pso) continue;");
+    m->AddLine("    TPaveStats* ps = static_cast<TPaveStats*>(pso);");
+    m->AddLine("    ps->SetOptFit(111);");
+    m->AddLine("    UShort_t qual = htmp->GetUniqueID();");
+    m->AddLine("    ps->SetFillColor(qual >= 3 ? kRed-4 : qual >= 2 ? kOrange-4 : qual >= 1 ? kYellow-4 : kGreen-4);");
+    // m->AddLine("    lf->Remove(lf->FindObject(\"fits\"));");
+    // m->AddLine("    ps->Paint();");
+    m->AddLine("    break;");
+    m->AddLine("  }");
+    // m->AddLine("  gPad->Modified(); gPad->Update(); gPad->cd();");
+    m->AddLine("}");
+
+    TObject* old = l->FindObject(m->GetName());
+    if (old) l->Remove(old);
+    l->Add(m);
+  }
+  
+  const Double_t kROErrorsLabelY    = 30.;
+  
+  const Int_t    kConvolutionSteps  = 100;
+  const Double_t kConvolutionNSigma = 5;
+
+  // 
+  // The shift of the most probable value for the ROOT function TMath::Landau 
+  //
+  const Double_t  kMpShift  = -0.22278298;
+  // 
+  // Integration normalisation 
+  //
+  const Double_t  kInvSq2Pi = 1. / TMath::Sqrt(2*TMath::Pi());
+
+  Double_t landau(Double_t x, Double_t delta, Double_t xi)
+  {
+    // 
+    // Calculate the shifted Landau
+    // @f[
+    //    f'_{L}(x;\Delta,\xi) = f_L(x;\Delta+0.22278298\xi)
+    // @f]
+    //
+    // where @f$ f_{L}@f$ is the ROOT implementation of the Landau
+    // distribution (known to have @f$ \Delta_{p}=-0.22278298@f$ for
+    // @f$\Delta=0,\xi=1@f$.
+    // 
+    // Parameters:
+    //    x      Where to evaluate @f$ f'_{L}@f$ 
+    //    delta  Most probable value 
+    //    xi     The 'width' of the distribution 
+    //
+    // Return:
+    //    @f$ f'_{L}(x;\Delta,\xi) @f$
+    //
+    return TMath::Landau(x, delta - xi * kMpShift, xi);
+  }
+  Double_t landauGaus(Double_t x, Double_t delta, Double_t xi,
+                     Double_t sigma, Double_t sigmaN)
+  {
+    // 
+    // Calculate the value of a Landau convolved with a Gaussian 
+    // 
+    // @f[ 
+    // f(x;\Delta,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}}
+    //    \int_{-\infty}^{+\infty} d\Delta' f'_{L}(x;\Delta',\xi)
+    //    \exp{-\frac{(\Delta-\Delta')^2}{2\sigma'^2}}
+    // @f]
+    // 
+    // where @f$ f'_{L}@f$ is the Landau distribution, @f$ \Delta@f$
+    // the energy loss, @f$ \xi@f$ the width of the Landau, and @f$
+    // \sigma'^2=\sigma^2-\sigma_n^2 @f$.  Here, @f$\sigma@f$ is the
+    // variance of the Gaussian, and @f$\sigma_n@f$ is a parameter
+    // modelling noise in the detector.
+    //
+    // Note that this function uses the constants kConvolutionSteps and
+    // kConvolutionNSigma
+    // 
+    // References: 
+    //  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
+    //  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
+    //  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
+    // 
+    // Parameters:
+    //    x         where to evaluate @f$ f@f$
+    //    delta     @f$ \Delta@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
+    //    xi        @f$ \xi@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
+    //    sigma     @f$ \sigma@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
+    //    sigma_n   @f$ \sigma_n@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
+    // 
+    // Return:
+    //    @f$ f@f$ evaluated at @f$ x@f$.  
+    //
+    Double_t deltap = delta - xi * kMpShift;
+    Double_t sigma2 = sigmaN*sigmaN + sigma*sigma;
+    Double_t sigma1 = sigmaN == 0 ? sigma : TMath::Sqrt(sigma2);
+    Double_t xlow   = x - kConvolutionNSigma * sigma1;
+    Double_t xhigh  = x + kConvolutionNSigma * sigma1;
+    Double_t step   = (xhigh - xlow) / kConvolutionSteps;
+    Double_t sum    = 0;
+  
+    for (Int_t i = 0; i <= kConvolutionSteps/2; i++) { 
+      Double_t x1 = xlow  + (i - .5) * step;
+      Double_t x2 = xhigh - (i - .5) * step;
+    
+      sum += TMath::Landau(x1, deltap, xi, kTRUE) * TMath::Gaus(x, x1, sigma1);
+      sum += TMath::Landau(x2, deltap, xi, kTRUE) * TMath::Gaus(x, x2, sigma1);
+    }
+    return step * sum * kInvSq2Pi / sigma1;
+  }
+
+  // 
+  // Utility function to use in TF1 defintition 
+  //
+  Double_t landauGaus1(Double_t* xp, Double_t* pp) 
+  {
+    Double_t x        = xp[0];
+    Double_t constant = pp[0];
+    Double_t delta    = pp[1];
+    Double_t xi       = pp[2];
+    Double_t sigma    = pp[3];
+    Double_t sigmaN   = pp[4];
+
+    return constant * landauGaus(x, delta, xi, sigma, sigmaN);
+  }
+
+  //____________________________________________________________________
+  TF1* makeLandauGaus(const char* , 
+                     Double_t  c=1, 
+                     Double_t  delta=.5,  Double_t xi=0.07, 
+                     Double_t  sigma=.1,  Double_t sigmaN=-1, 
+                     Double_t  xmin=0,    Double_t xmax=15)
+  {
+    // 
+    // Generate a TF1 object of @f$ f_I@f$ 
+    // 
+    // Parameters:
+    //    c        Constant
+    //    delta    @f$ \Delta@f$ 
+    //    xi       @f$ \xi_1@f$               
+    //    sigma    @f$ \sigma_1@f$            
+    //    sigma_n  @f$ \sigma_n@f$            
+    //    xmin     Least value of range
+    //    xmax     Largest value of range
+    // 
+    // Return:
+    //    Newly allocated TF1 object
+    //
+    Int_t npar     = 5;
+    TF1*  func     = new TF1("landauGaus", 
+                            &landauGaus1,xmin,xmax,npar);
+    // func->SetLineStyle(((i-2) % 10)+2); // start at dashed
+    func->SetLineColor(kBlack); 
+    func->SetLineWidth(2);
+    func->SetNpx(500);
+    func->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
+
+    // Set the initial parameters from the seed fit 
+    func->SetParameter(0,      c);       
+    func->SetParameter(1,  delta);   
+    func->SetParameter(2,     xi);      
+    func->SetParameter(3,  sigma);   
+    func->SetParameter(4, sigmaN); 
+
+    func->SetParLimits(1, 0,    xmax);
+    func->SetParLimits(2, 0,    xmax);
+    func->SetParLimits(3, 0.01, 0.4);
+    
+    if (sigmaN < 0) func->FixParameter(4, 0);
+    else            func->SetParLimits(4, 0, xmax);
+
+    return func;
+  }
+}
+  
+//__________________________________________________________________
+AliFMDQAChecker::AliFMDQAChecker() 
+  : AliQACheckerBase("FMD","FMD Quality Assurance Checker") ,
+    fDoScale(false),
+    fDidExternal(false), 
+    fShowFitResults(true),
+    fELossLowCut(0.2), 
+    fELossNRMS(3), 
+    fELossBadChi2Nu(10), 
+    fELossFkupChi2Nu(100), 
+    fELossMinEntries(1000),
+    fELossGoodParError(0.1),
+    fROErrorsBad(0.3), 
+    fROErrorsFkup(0.5)
+{
+}          
+
+//__________________________________________________________________
+void
+AliFMDQAChecker::ProcessExternalParams()
+{
+  ProcessExternalParam("ELossLowCut",          fELossLowCut);
+  ProcessExternalParam("ELossNRMS",            fELossNRMS);
+  ProcessExternalParam("ELossBadChi2Nu",       fELossBadChi2Nu);
+  ProcessExternalParam("ELossFkupChi2Nu",      fELossFkupChi2Nu);
+  ProcessExternalParam("ELossGoodParError",    fELossGoodParError);
+  ProcessExternalParam("ROErrorsBad",          fROErrorsBad);
+  ProcessExternalParam("ROErrorsFkup",         fROErrorsFkup);
+  Double_t tmp = 0;
+  ProcessExternalParam("CommonScale",          tmp);
+  fDoScale = tmp > 0;
+  ProcessExternalParam("ELossMinEntries",      tmp);
+  fELossMinEntries = tmp;
+
+  GetThresholds();
+
+  fDidExternal = true;
+}
+//__________________________________________________________________
+void
+AliFMDQAChecker::ProcessExternalParam(const char* name, Double_t& v)
+{
+  TObject* o = fExternParamList->FindObject(name);
+  if (!o) return; 
+  TParameter<double>* p = static_cast<TParameter<double>*>(o);
+  v = p->GetVal();
+  AliDebugF(3, "External parameter: %-20s=%lf", name, v);
+}
+
+//__________________________________________________________________
+void
+AliFMDQAChecker::GetThresholds()
+{
+  const char*    path   = "GRP/Calib/QAThresholds";
+  AliCDBManager* cdbMan = AliCDBManager::Instance();
+  AliCDBEntry*   cdbEnt = cdbMan->Get(path);
+  if (!cdbEnt) { 
+    AliWarningF("Failed to get CDB entry at %s", path);
+    return;
+  }
+  
+  TObjArray*     cdbObj = static_cast<TObjArray*>(cdbEnt->GetObject());
+  if (!cdbObj) { 
+    AliWarningF("Failed to get CDB object at %s", path);
+    return;
+  }
+  
+  TObject*       fmdObj = cdbObj->FindObject("FMD");
+  if (!fmdObj) { 
+    AliWarningF("Failed to get FMD object at from CDB %s", path);
+    return;
+  }
+  
+  AliQAThresholds* qaThr = static_cast<AliQAThresholds*>(fmdObj);
+  Int_t nThr = qaThr->GetSize();
+  for (Int_t i = 0; i < nThr; i++) { 
+    TObject* thr = qaThr->GetThreshold(i);
+    if (!thr) continue; 
+
+    TParameter<double>* d = dynamic_cast<TParameter<double>*>(thr);
+    if (!d) { 
+      AliWarningF("Parameter %s not of type double", thr->GetName());
+      continue;
+    }
+    Double_t val = d->GetVal();
+    TString  name(thr->GetName());
+    if      (name.EqualTo("ELossBadChi2Nu"))     fELossBadChi2Nu    = val;
+    else if (name.EqualTo("ELossFkupChi2Nu"))    fELossFkupChi2Nu   = val;
+    else if (name.EqualTo("ELossGoodParError"))  fELossGoodParError = val;
+    else if (name.EqualTo("ROErrorsBad"))        fROErrorsBad       = val;
+    else if (name.EqualTo("ROErrorsFkup"))       fROErrorsFkup      = val;    
+    AliDebugF(3, "Threshold %s=%f", name.Data(), val);
+  }
+}
+
+//__________________________________________________________________
+AliQAv1::QABIT_t
+AliFMDQAChecker::Quality2Bit(UShort_t value) const
+{
+  AliQAv1::QABIT_t  ret   = AliQAv1::kINFO; // Assume success 
+  if      (value >= kWhatTheFk) ret = AliQAv1::kFATAL;
+  else if (value >= kBad)       ret = AliQAv1::kERROR;
+  else if (value >= kProblem)   ret = AliQAv1::kWARNING;
+  
+  return ret;
+}
+
+//__________________________________________________________________
+void
+AliFMDQAChecker::SetQA(AliQAv1::ALITASK_t index, Double_t* values) const
+{
+  AliQAv1 * qa = AliQAv1::Instance(index) ;
+
+  for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) {
+    // Check if specie is defined 
+    if (!qa->IsEventSpecieSet(AliRecoParam::ConvertIndex(specie)))
+      continue ;
+    
+    // No checker is implemented, set all QA to Fatal
+    if (!values) { 
+      qa->Set(AliQAv1::kFATAL, specie) ; 
+      continue;
+    }
+    
+    UShort_t          value = values[specie];
+    AliQAv1::QABIT_t  ret   = Quality2Bit(value);
+    
+    qa->Set(ret, AliRecoParam::ConvertIndex(specie));
+    AliDebugF(3, "Quality of %s: %d -> %d", 
+             AliRecoParam::GetEventSpecieName(specie), value, ret);
+  }
+}
+
+//__________________________________________________________________
+UShort_t
+AliFMDQAChecker::BasicCheck(TH1* hist) const
+{
+  if (hist->GetEntries() <= 0) return kOK;
+  return (hist->GetMean() > 0 ? kOK : kProblem);
+}
+
 //__________________________________________________________________
-Double_t 
+UShort_t
 AliFMDQAChecker::CheckOne(AliQAv1::ALITASK_t          what,
                          AliRecoParam::EventSpecie_t specie, 
                          TH1*                        hist) const
@@ -63,32 +397,303 @@ AliFMDQAChecker::CheckOne(AliQAv1::ALITASK_t          what,
   return 0;
 }
 //__________________________________________________________________
-Double_t 
+UShort_t 
 AliFMDQAChecker::CheckESD(AliRecoParam::EventSpecie_t /* specie*/, 
                          TH1*                        hist) const
 {
-  return (hist->GetMean() > 0 ? 1 : 0);
+  return BasicCheck(hist);
+}
+//__________________________________________________________________
+UShort_t
+AliFMDQAChecker::CheckFit(TH1* hist, const TFitResultPtr& res, 
+                         Double_t low, Double_t high, Int_t& color) const
+{
+  color = kGreen+4;
+
+  UShort_t   ret   = 0;
+  Int_t      nPar  = res->NPar();
+  Double_t   dy    = .06;
+  Double_t   x     = .2;
+  Double_t   y     = .9-dy;
+  Double_t   chi2  = res->Chi2();
+  Int_t      nu    = res->Ndf();
+  Double_t   red   = (nu == 0 ? fELossFkupChi2Nu : chi2 / nu);
+  TObjArray* lines = 0;
+  TLatex*    lRed  = 0;
+  TLatex*    ltx   = 0;
+  if (fShowFitResults) { 
+    lines = new TObjArray(nPar+3);
+    lines->SetName("lines");
+    lines->SetOwner(true);
+    
+    ltx = new TLatex(x, y, Form("#chi^{2}/#nu: %7.3f",red));
+    ltx->SetNDC(true);
+    ltx->SetTextColor(color);
+    ltx->SetTextSize(dy-.01);
+    lines->Add(ltx);
+    lRed = ltx;
+    
+    Double_t x1 = .85;
+    Double_t y1 = .5;
+    ltx = new TLatex(x1, y1, Form("[thresholds: %6.2f, %6.2f]", 
+                               fELossBadChi2Nu, fELossFkupChi2Nu));
+    ltx->SetTextColor(kGray+3);
+    ltx->SetTextSize(dy-.01);
+    ltx->SetNDC(true);
+    ltx->SetTextAlign(31);
+    lines->Add(ltx);    
+
+    y1 -= dy;
+    ltx = new TLatex(x1, y1, Form("Fit range: [%6.2f,%6.2f]", low, high));
+    ltx->SetTextColor(kGray+3);
+    ltx->SetTextSize(dy-.01);
+    ltx->SetTextAlign(31);
+    ltx->SetNDC(true);
+    lines->Add(ltx);
+
+    y1 -= dy;
+    ltx = new TLatex(x1, y1, Form("Entries: %d", 
+                                  Int_t(hist->GetEffectiveEntries())));
+    ltx->SetTextColor(kGray+3);
+    ltx->SetTextSize(dy-.01);
+    ltx->SetTextAlign(31);
+    ltx->SetNDC(true);
+    lines->Add(ltx);
+  }
+  
+  if (red > fELossBadChi2Nu) { // || res->Prob() < .01) { 
+    AliWarningF("Fit gave chi^2/nu=%f/%d=%f>%f (%f)", 
+               res->Chi2(), res->Ndf(), red, fELossBadChi2Nu, 
+               fELossFkupChi2Nu);
+    if (lRed) lRed->SetTextColor(kOrange+2);
+    res->Print();
+    ret++;
+    if (red > fELossFkupChi2Nu) { 
+      if (lRed) lRed->SetTextColor(kRed+2);
+      ret++;
+    }
+  }
+  // Now check the relative error on the fit parameters 
+  Int_t parsOk = 0;
+  for (Int_t i = 0; i < nPar; i++) { 
+    if (res->IsParameterFixed(i)) continue; 
+    Double_t pv  = res->Parameter(i);
+    Double_t pe  = res->ParError(i);
+    Double_t rel = (pv == 0 ? 100 : pe / pv);
+    if (lines) {
+      y -= dy;
+      ltx = new TLatex(x, y, Form("#delta%s/%s: %7.3f", 
+                                 res->ParName(i).c_str(),
+                                 res->ParName(i).c_str(),
+                                 /*pv, pe,*/ rel));
+      ltx->SetNDC(true);
+      ltx->SetTextColor(color);
+      ltx->SetTextSize(dy-.01);
+      lines->Add(ltx);
+    }
+    if (i == 3) continue; // Skip sigma 
+    Double_t thr = fELossGoodParError;
+    if (rel < thr) parsOk++;
+    else if (ltx) ltx->SetLineColor(kOrange+2);
+  }
+  if (parsOk > 0) 
+    ret = TMath::Max(ret-(parsOk-1),0);
+  if (ret > 1) color = kRed+2;
+  if (ret > 0) color = kOrange+2;
+
+  TList*   lf  = hist->GetListOfFunctions();
+  TObject* old = lf->FindObject(lines->GetName());
+  if (old) {
+    lf->Remove(old);
+    delete old;
+  }
+  lf->Add(lines);
+  hist->SetStats(false);
+    
+  return ret;
+}
+
+//__________________________________________________________________
+void
+AliFMDQAChecker::AddFitResults(TH1* hist, const TFitResultPtr& res, 
+                              Int_t color, Double_t low, Double_t high) const
+{
+  if (!fShowFitResults) return; 
+
+  Int_t      nPar  = res->NPar();
+  TObjArray* lines = new TObjArray(nPar+1);
+  lines->SetOwner(kTRUE);
+  lines->SetName("fitResults");
+
+  Double_t   dy    = .06;
+  Double_t   x     = .2;
+  Double_t   y     = .9-dy;
+  Double_t   chi2  = res->Chi2();
+  Int_t      nu    = res->Ndf();
+  Double_t   red   = (nu == 0 ? fELossFkupChi2Nu : chi2 / nu);
+  
+  TLatex* ltx = new TLatex(x, y, Form("#chi^{2}/#nu: %7.3f",red));
+  ltx->SetNDC(true);
+  ltx->SetTextColor(color);
+  ltx->SetTextSize(dy-.01);
+  lines->Add(ltx);
+  
+  y -= dy;
+  ltx = new TLatex(x, y, Form("[thresholds: %6.2f, %6.2f]", 
+                             fELossBadChi2Nu, fELossFkupChi2Nu));
+  ltx->SetTextColor(kGray+3);
+  ltx->SetTextSize(dy-.01);
+  ltx->SetNDC(true);
+  lines->Add(ltx);
+
+  y -= dy;
+  ltx = new TLatex(x, y, Form("Fit range: [%6.2f,%6.2f]", low, high));
+  ltx->SetTextColor(kGray+3);
+  ltx->SetTextSize(dy-.01);
+  ltx->SetNDC(true);
+  lines->Add(ltx);
+
+  for (Int_t i = 0; i < nPar; i++) { 
+    if (res->IsParameterFixed(i)) continue; 
+    y -= dy;
+    Double_t pv  = res->Parameter(i);
+    Double_t pe  = res->ParError(i);
+    Double_t rel = (pv == 0 ? 100 : pe / pv);
+    ltx = new TLatex(x, y, Form("#delta%s/%s: %7.3f", 
+                               res->ParName(i).c_str(),
+                               res->ParName(i).c_str(),
+                               /*pv, pe,*/ rel));
+    ltx->SetNDC(true);
+    ltx->SetTextColor(color);
+    ltx->SetTextSize(dy-.01);
+    lines->Add(ltx);
+  }
+  TList*   lf  = hist->GetListOfFunctions();
+  TObject* old = lf->FindObject(lines->GetName());
+  if (old) {
+    lf->Remove(old);
+    delete old;
+  }
+  lf->Add(lines);
+  hist->SetStats(false);
 }
+
 //__________________________________________________________________
-Double_t 
+UShort_t
 AliFMDQAChecker::CheckRaw(AliRecoParam::EventSpecie_t /* specie*/, 
                          TH1*                        hist) const
 {
-  return (hist->GetMean() > 0 ? 1 : 0);
+  Int_t ret = BasicCheck(hist);
+  TString name(hist->GetName());
+  if (name.Contains("readouterrors", TString::kIgnoreCase)) { 
+    // Check the mean number of errors per event 
+    TH2*  roErrors = static_cast<TH2*>(hist);
+    Int_t nY       = roErrors->GetNbinsY();
+
+    TLatex* ltx = new TLatex(.15, .8, Form("Thresholds: %5.2f,%5.2f",
+                                          fROErrorsBad, fROErrorsFkup));
+    ltx->SetName("thresholds");
+    ltx->SetTextColor(kGray+3);
+    ltx->SetNDC();
+
+    TList* ll = hist->GetListOfFunctions(); 
+    TObject* old = ll->FindObject(ltx->GetName());
+    if (old) { 
+      ll->Remove(old);
+      delete old;
+    }
+    ll->Add(ltx);
+    
+    for (Int_t i = 1; i <= 3; i++) {
+      Double_t sum = 0;
+      Int_t    cnt = 0;
+      for (Int_t j = 1; j <= nY; j++) {
+       Int_t n =  roErrors->GetBinContent(i, j);
+       sum     += n * roErrors->GetYaxis()->GetBinCenter(j);
+       cnt     += n;
+      }
+      Double_t mean = sum / cnt;
+
+      ltx = new TLatex(i, kROErrorsLabelY, Form("Mean: %6.3f", mean));
+      ltx->SetName(Form("FMD%d", i));
+      ltx->SetTextAngle(90);
+      ltx->SetTextColor(kGreen+4);
+      old = ll->FindObject(ltx->GetName());
+      if (old) { 
+       ll->Remove(old);
+       delete old;
+      }
+      ll->Add(ltx);
+
+      if (mean > fROErrorsBad) { 
+       AliWarningF("Mean of readout errors for FMD%d = %f > %f (%f)", 
+                   i, mean, fROErrorsBad, fROErrorsFkup);
+       ret++;
+       ltx->SetTextColor(kOrange+2);
+       if (mean > fROErrorsFkup) {
+         ret++;
+         ltx->SetTextColor(kRed+2);
+       }
+      }
+    }
+  }
+  else if (name.Contains("eloss",TString::kIgnoreCase)) { 
+    // Try to fit a function to the histogram 
+    if (hist->GetEntries() < 1000) return ret;
+
+    Double_t xMin  = hist->GetXaxis()->GetXmin();
+    Double_t xMax  = hist->GetXaxis()->GetXmax();
+
+    hist->GetXaxis()->SetRangeUser(fELossLowCut, xMax);
+    Int_t    bMaxY = hist->GetMaximumBin();
+    Double_t xMaxY = hist->GetXaxis()->GetBinCenter(bMaxY);
+    Double_t rms   = hist->GetRMS();
+    Double_t low   = hist->GetXaxis()->GetBinCenter(bMaxY-4);
+    hist->GetXaxis()->SetRangeUser(0.2, xMaxY+(fELossNRMS+1)*rms);
+    rms  = hist->GetRMS();
+    hist->GetXaxis()->SetRange(0,-1);
+    TF1*          func = makeLandauGaus(name);
+    func->SetParameter(1, xMaxY);
+    func->SetLineColor(kGreen+4);
+    // func->SetLineStyle(2);
+    Double_t high = xMax; // xMaxY+fELossNRMS*rms;
+
+    TFitResultPtr res  = hist->Fit(func, "QS", "", low, high);
+
+    Int_t    color = func->GetLineColor();
+    UShort_t qual  = CheckFit(hist, res, low, high, color);
+
+    // Make sure we save the function in the full range of the histogram
+    func = hist->GetFunction("landauGaus");
+    func->SetRange(xMin, xMax);
+    // func->SetParent(hist);
+    func->Save(xMin, xMax, 0, 0, 0, 0);
+    func->SetLineColor(color);
+    
+    if (qual > 0) { 
+      func->SetLineWidth(3);
+      func->SetLineStyle(1);
+      if (qual > 1) 
+       func->SetLineWidth(4);
+    }
+    ret += qual;
+  }
+
+  return ret;
 }
 //__________________________________________________________________
-Double_t 
+UShort_t
 AliFMDQAChecker::CheckSim(AliRecoParam::EventSpecie_t /* specie*/, 
                          TH1*                        hist) const
 {
-  return (hist->GetMean() > 0 ? 1 : 0);
+  return BasicCheck(hist);
 }
 //__________________________________________________________________
-Double_t 
+UShort_t
 AliFMDQAChecker::CheckRec(AliRecoParam::EventSpecie_t /* specie*/, 
                          TH1*                        hist) const
 {
-  return (hist->GetMean() > 0 ? 1 : 0);
+  return BasicCheck(hist);
 }
 
 //__________________________________________________________________
@@ -107,8 +712,17 @@ void AliFMDQAChecker::Check(Double_t*                   rv,
   //         each 'specie'
   //    t    Reconstruction parameters - not used. 
   //
-  
+  // The bounds defined for RV are 
+  // 
+  //   FATAL:      [-1,   0.0]
+  //   ERROR:      (0.0,  0.002]
+  //   WARNING:    (0.002,0.5]
+  //   INFO:       (0.5,  1.0]
+  //
   // Double_t* rv = new Double_t[AliRecoParam::kNSpecies] ; 
+
+  if (!fDidExternal) ProcessExternalParams();
+
   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) {
     // Int_t count   = 0;
     rv[specie]    = 0.; 
@@ -120,12 +734,20 @@ void AliFMDQAChecker::Check(Double_t*                   rv,
     
     TH1* hist  = 0;
     Int_t nHist = list[specie]->GetEntriesFast();
+    UShort_t ret = 0;
     for(Int_t i= 0; i< nHist; i++) {
-      
       if (!(hist = static_cast<TH1*>(list[specie]->At(i)))) continue;
-      
-      rv[specie] += CheckOne(what, AliRecoParam::ConvertIndex(specie), hist);
+      Int_t qual = CheckOne(what, AliRecoParam::ConvertIndex(specie), hist);
+      hist->SetUniqueID(Quality2Bit(qual));
+      ret += qual;
     } // for (int i ...)
+    rv[specie] = ret;
+    // if      (ret > kWhatTheFk) rv[specie] = fLowTestValue[AliQAv1::kFATAL];
+    // else if (ret > kBad)       rv[specie] = fUpTestValue[AliQAv1::kERROR]; 
+    // else if (ret > kProblem)   rv[specie] = fUpTestValue[AliQAv1::kWARNING]; 
+    // else                       rv[specie] = fUpTestValue[AliQAv1::kINFO]; 
+    AliDebugF(3, "Combined sum is %d -> %f", ret, rv[specie]);
+
     // if (count != 0) rv[specie] /= count;
   }
   // return rv;
@@ -219,13 +841,13 @@ AliFMDQAChecker::MakeImage(TObjArray** list,
        FindMinMax(hist, hMin, hMax);
        max = TMath::Max(max, hMax);
        min = TMath::Min(min, hMin);
-       AliInfoF("Min/max of %40s: %f/%f, global -> %f/%f", 
-                hist->GetName(), hMin, hMax, min, max);
+       // AliInfoF("Min/max of %40s: %f/%f, global -> %f/%f", 
+       //          hist->GetName(), hMin, hMax, min, max);
       }
     }
     break ; 
   }
-  AliInfoF("Global min/max=%f/%f", min, max);
+  // AliInfoF("Global min/max=%f/%f", min, max);
   min = TMath::Max(1e-1, min);
   max = TMath::Max(1e5,  max);
 
@@ -274,7 +896,6 @@ AliFMDQAChecker::MakeImage(TObjArray** list,
                         AliQAv1::GetImageFileFormat()));
     fImage[specie]->Print(outName, "ps") ; 
 
-
     // Now set some parameters on the canvas 
     fImage[specie]->Clear(); 
     fImage[specie]->SetTopMargin(0.10);
@@ -327,7 +948,7 @@ AliFMDQAChecker::MakeImage(TObjArray** list,
       logOpts |= CheckForLog(hist->GetZaxis(), pad, 3);
 
       // Figure out special cases 
-      TString opt;
+      TString opt("");
       TString name(hist->GetName());
       if (name.Contains("readouterrors", TString::kIgnoreCase)) {
        pad->SetRightMargin(0.15);
@@ -349,17 +970,19 @@ AliFMDQAChecker::MakeImage(TObjArray** list,
       }
       // Draw (As a copy)
       hist->DrawCopy(opt);
-
+      
       // Special cases 
       if (name.Contains("readouterrors", TString::kIgnoreCase)) {
+#if 0
        for (Int_t kk = 1; kk <= 3; kk++) {
          TH1* proj = static_cast<TH2*>(hist)->ProjectionY("",kk,kk);
-       Double_t m = proj->GetMean(); 
-       TLatex* l = new TLatex(kk, 30, Form("Mean: %f", m));
-       l->SetTextAngle(90);
-       l->SetTextColor(m > 10 ? kRed+1 : m > .7 ? kOrange+2 :kGreen+2);
-       l->Draw();
-      }
+         Double_t m = proj->GetMean(); 
+         TLatex* l = new TLatex(kk, 30, Form("Mean: %f", m));
+         l->SetTextAngle(90);
+         l->SetTextColor(m > 10 ? kRed+1 : m > .7 ? kOrange+2 :kGreen+2);
+         l->Draw();
+       }
+#endif
       }
       else {
        gStyle->SetOptTitle(0);
@@ -390,6 +1013,11 @@ AliFMDQAChecker::MakeImage(TObjArray** list,
     }
     // Print to a post-script file 
     fImage[specie]->Print(outName, "ps");
+#if 0
+    fImage[specie]->Print(Form("%s_%d.png", 
+                              AliRecoParam::GetEventSpecieName(specie), 
+                              AliQAChecker::Instance()->GetRunNumber()));
+#endif
   }
 }
 
index d65e212..7f11b1e 100644 (file)
@@ -9,17 +9,17 @@
 class TFile; 
 class TH1F; 
 class TH1I; 
+class TFitResultPtr;
 
-/** @class AliFMDQAChecker 
-    @brief Quality assurance checker for the FMD */
+/** 
+ * @class AliFMDQAChecker 
+ * @brief Quality assurance checker for the FMD 
+ */
 class AliFMDQAChecker : public AliQACheckerBase 
 {
 public:
   /** Constructor */
-  AliFMDQAChecker() 
-    : AliQACheckerBase("FMD","FMD Quality Assurance Checker") ,
-      fDoScale(false)
-  {}          
+  AliFMDQAChecker();
   /** Destructor */
   virtual ~AliFMDQAChecker() {}
   /** 
@@ -45,33 +45,173 @@ public:
   void  MakeImage(TObjArray** list, 
                  AliQAv1::TASKINDEX_t task, 
                  AliQAv1::MODE_t mode);
+  /** 
+   * Set wether to scale the histograms to a common Y axis or not when
+   * generating plots
+   * 
+   * @param on If true, do scale 
+   */
   void SetDoScale(Bool_t on=true) { fDoScale = on; }
-
 protected:
+  // Return values 
+  enum { 
+    kOK,
+    kProblem, 
+    kBad, 
+    kWhatTheFk
+  };
   /** 
    * Check one histogram 
    * 
-   * @param specie 
-   * @param hist 
+   * @param specie Event specie 
+   * @param hist   Histogram to check 
    * 
-   * @return 
+   * @return 0 if all is good, increasing severity for increasingly
+   * bad data
    */
-  Double_t CheckOne(AliQAv1::ALITASK_t          what,
+  UShort_t CheckOne(AliQAv1::ALITASK_t          what,
                    AliRecoParam::EventSpecie_t specie, 
                    TH1*                        hist) const;
-  Double_t CheckRaw(AliRecoParam::EventSpecie_t specie, 
+  /** 
+   * Check raw input.   If the reconstructor is enabled, then we try to fit 
+   * 
+   * @f[ 
+   *   f(\Delta;\Delta_p,\xi,\sigma) = \int_{-\infty}^{\infty}d\Delta_p' 
+   *     f_{L}(\Delta;\Delta_p',\xi) \times 
+   *     e^{\frac{(\Delta-\Delta_p')^2}{\sigma^2}}
+   * @f]
+   * 
+   * where @f$f_L@f$ is the Landau distribution, to the data.  The
+   * quality is set according to the value of @f$\chi^2/\nu@f$.
+   *
+   * If no reconstructor is set, then simply check that the histogram
+   * isn't empty
+   * 
+   * @param specie Event specie 
+   * @param hist   Histogram to check 
+   * 
+   * @return 0 if all is good, increasing severity for increasingly
+   * bad data
+   */
+  UShort_t CheckRaw(AliRecoParam::EventSpecie_t specie, 
                    TH1*                        hist) const;
-  Double_t CheckSim(AliRecoParam::EventSpecie_t specie, 
+  /** 
+   * Check simulation output.  Does a simple test of whether the
+   * histogram is empty or not.
+   * 
+   * @param specie Event specie 
+   * @param hist   Histogram to check 
+   * 
+   * @return 0 if all is good, increasing severity for increasingly
+   * bad data
+   */
+  UShort_t CheckSim(AliRecoParam::EventSpecie_t specie, 
                    TH1*                        hist) const;
-  Double_t CheckESD(AliRecoParam::EventSpecie_t specie, 
+  /** 
+   * Check ESD.  Does a simple test of whether the histogram is empty
+   * or not.
+   * 
+   * @param specie Event specie 
+   * @param hist   Histogram to check 
+   * 
+   * @return 0 if all is good, increasing severity for increasingly
+   * bad data
+   */
+  UShort_t CheckESD(AliRecoParam::EventSpecie_t specie, 
                    TH1*                        hist) const;
-  Double_t CheckRec(AliRecoParam::EventSpecie_t specie, 
+  /** 
+   * Check reconstruction points.  Does a simple test of whether the
+   * histogram is empty or not.
+   * 
+   * @param specie Event specie 
+   * @param hist   Histogram to check 
+   * 
+   * @return 0 if all is good, increasing severity for increasingly
+   * bad data
+   */
+  UShort_t CheckRec(AliRecoParam::EventSpecie_t specie, 
                    TH1*                        hist) const;
-
-  Bool_t fDoScale;
+  /** 
+   * Set the returned QA from this checker based on the values in the
+   * array @a values.  Note, this by-passes the Low/High setting of
+   * the base class (which are very confusing) 
+   * 
+   * @param index   Task index 
+   * @param values  Array of values 
+   */
+  void SetQA(AliQAv1::ALITASK_t index, Double_t* values) const;
+  /** 
+   * Process external parameters 
+   * 
+   */
+  void ProcessExternalParams();
+  /** 
+   * Process a single external parameter 
+   * 
+   * @param name Name of parameter 
+   * @param v    On return, the value - as a double 
+   */
+  void ProcessExternalParam(const char* name, Double_t& v);
+  /** 
+   * Get the thresholds from OCDB
+   * 
+   */
+  void GetThresholds();
+  /** 
+   * The basic check on a histogram 
+   * 
+   * @param hist Histogram 
+   * 
+   * @return 1 empty - 0 otherwise 
+   */
+  UShort_t BasicCheck(TH1* hist) const;
+  /** 
+   * Translate our internal quality measure to QA framework bit 
+   * 
+   * @param qual Internal quality 
+   * 
+   * @return QA framework quality bit 
+   */
+  AliQAv1::QABIT_t Quality2Bit(UShort_t qual) const;
+  /** 
+   * Add fit results to to plot 
+   * 
+   * @param hist  Histogram 
+   * @param res   Fit result 
+   * @param color Color to use for the text - depends on quality 
+   * @param low   Lower bound on fit range 
+   * @param high  Upper bound on fit range 
+   */
+  void AddFitResults(TH1* hist, const TFitResultPtr& res, Int_t color,
+                    Double_t low, Double_t high) const;
+  UShort_t CheckFit(TH1* hist, const TFitResultPtr& res, 
+                   Double_t low, Double_t high, Int_t& color) const;
+  Bool_t   fDoScale;         // Whether to scale all histograms 
+  Bool_t   fDidExternal;     // Whether we've processed the external params 
+  Bool_t   fShowFitResults;  // Whether to put the fit result on the plots
+  Double_t fELossLowCut;     // Low cut on ELoss fits 
+  Double_t fELossNRMS;       // Number of RMS to fit upward
+  Double_t fELossBadChi2Nu;  // Cut on bad chi2/nu
+  Double_t fELossFkupChi2Nu; // Cut on F**ked up chi2/nu
+  Int_t    fELossMinEntries; // Least number of entries before fitting
+  Double_t fELossGoodParError; // Least relative error
+  Double_t fROErrorsBad;     // Cut on read-out errors 
+  Double_t fROErrorsFkup;    // Cut on read-out errors 
 private:
-  AliFMDQAChecker(const AliFMDQAChecker& qac); // cpy ctor   
-  AliFMDQAChecker &operator=(const AliFMDQAChecker& qac); // assignment operator
+  /** 
+   * Copy constructor - not implemented 
+   * 
+   * @param qac Object to copy from
+   */
+  AliFMDQAChecker(const AliFMDQAChecker& qac); 
+  /** 
+   * assignment operator - not implemented 
+   * 
+   * @param qac Object to assign from 
+   * 
+   * @return Reference to this object 
+   */
+  AliFMDQAChecker& operator=(const AliFMDQAChecker& qac); 
 
   ClassDef(AliFMDQAChecker,0)  // Yves? what to do? 
 };
index 4959bac..c5b1b59 100644 (file)
 #include <TH1F.h> 
 #include <TH1I.h> 
 #include <TH2I.h> 
+#include <TGeoManager.h>
 
 // --- AliRoot header files ---
 #include "AliESDEvent.h"
 #include "AliLog.h"
+#include "AliGeomManager.h"
 #include "AliFMDQADataMakerRec.h"
 #include "AliFMDDigit.h"
 #include "AliFMDRecPoint.h"
@@ -30,6 +32,7 @@
 #include "AliESDFMD.h"
 #include "AliFMDParameters.h"
 #include "AliFMDRawReader.h"
+#include "AliFMDReconstructor.h"
 #include "AliRawReader.h"
 #include "AliFMDAltroMapping.h"
 #include "AliFMDDebug.h"
@@ -54,10 +57,12 @@ ClassImp(AliFMDQADataMakerRec)
 #endif
            
 //_____________________________________________________________________
-AliFMDQADataMakerRec::AliFMDQADataMakerRec() : 
-  AliQADataMakerRec(AliQAv1::GetDetName(AliQAv1::kFMD), 
-                   "FMD Quality Assurance Data Maker"),
-  fRecPointsArray("AliFMDRecPoint", 1000)
+AliFMDQADataMakerRec::AliFMDQADataMakerRec() 
+  : AliQADataMakerRec(AliQAv1::GetDetName(AliQAv1::kFMD), 
+                     "FMD Quality Assurance Data Maker"),
+    fRecPointsArray("AliFMDRecPoint", 1000), 
+    fReconstructor(0),
+    fUseReconstructor(true)
 {
   // ctor
  
@@ -67,7 +72,9 @@ AliFMDQADataMakerRec::AliFMDQADataMakerRec() :
 AliFMDQADataMakerRec::AliFMDQADataMakerRec(const AliFMDQADataMakerRec& qadm) 
   : AliQADataMakerRec(AliQAv1::GetDetName(AliQAv1::kFMD), 
                      "FMD Quality Assurance Data Maker"),
-    fRecPointsArray(qadm.fRecPointsArray)
+    fRecPointsArray(qadm.fRecPointsArray), 
+    fReconstructor(qadm.fReconstructor),
+    fUseReconstructor(qadm.fUseReconstructor)
 {
   // copy ctor 
   // Parameters: 
@@ -87,8 +94,9 @@ AliFMDQADataMakerRec::operator = (const AliFMDQADataMakerRec& qadm )
   // Return:
   //    Reference to this
   //
-  fRecPointsArray = qadm.fRecPointsArray;
-  
+  fRecPointsArray   = qadm.fRecPointsArray;
+  fReconstructor    = qadm.fReconstructor;
+  fUseReconstructor = qadm.fUseReconstructor;
   return *this;
 }
 //_____________________________________________________________________
@@ -159,13 +167,14 @@ TH1* AliFMDQADataMakerRec::MakeELossHist(UShort_t d, Char_t r, Short_t b)
       title.Append(Form("[0x%02x]", b));
     }
   }
-  TH1* hist = new TH1F(name, title,300,0, 15);
-  hist->SetXTitle("#Delta E/#Delta_{mip}");
+  TH1* hist = new TH1F(name, title,600,0, 15);
+  hist->SetXTitle("#Delta/#Delta_{mip}");
   hist->SetYTitle("Events [log]");
   hist->SetFillStyle(3001);
   hist->SetFillColor(color);
   hist->SetLineColor(color);
   hist->SetMarkerColor(color);
+  hist->Sumw2();
   hist->SetDirectory(0);
   // hist->SetStats(0);
 
@@ -177,6 +186,7 @@ TH1* AliFMDQADataMakerRec::MakeELossHist(UShort_t d, Char_t r, Short_t b)
 void AliFMDQADataMakerRec::InitESDs()
 {
   // create Digits histograms in Digits subdir
+  Info("InitESDs", "Initializing ESDs");
   const Bool_t expert   = kTRUE ; 
   const Bool_t image    = kTRUE ; 
   
@@ -189,6 +199,7 @@ void AliFMDQADataMakerRec::InitESDs()
 void AliFMDQADataMakerRec::InitDigits()
 {
   // create Digits histograms in Digits subdir
+  Info("InitDigits", "Initializing Digits");
   const Bool_t expert   = kTRUE ; 
   const Bool_t image    = kTRUE ; 
   
@@ -201,6 +212,7 @@ void AliFMDQADataMakerRec::InitDigits()
 void AliFMDQADataMakerRec::InitRecPoints()
 {
   // create Reconstructed Points histograms in RecPoints subdir
+  Info("InitRecPoints", "Initializing RecPoints");
   const Bool_t expert   = kTRUE ; 
   const Bool_t image    = kTRUE ; 
 
@@ -213,6 +225,7 @@ void AliFMDQADataMakerRec::InitRecPoints()
 void AliFMDQADataMakerRec::InitRaws()
 {
   // create Raws histograms in Raws subdir  
+  Info("InitRaws", "Initializing Raws");
   const Bool_t expert   = kTRUE ; 
   const Bool_t saveCorr = kTRUE ; 
   const Bool_t image    = kTRUE ; 
@@ -227,19 +240,39 @@ void AliFMDQADataMakerRec::InitRaws()
   Add2RawsList(hErrors, 1, !expert, image, !saveCorr);
   //AliInfo(Form("Adding %30s to raw list @ %2d", hErrors->GetName(), 1));
 
+  if (fUseReconstructor && !fReconstructor) {
+    // Int_t oldDbg = AliLog::GetDebugLevel("FMD","");
+    // AliLog::SetModuleDebugLevel("FMD", 5);
+
+    if (!gGeoManager) {
+      Info("InitRaws", "Loading the geometry");
+      AliGeomManager::LoadGeometry();
+    }
+
+    fReconstructor = new AliFMDReconstructor();
+    fReconstructor->SetDiagnose(false);
+    fReconstructor->Init();
+    // AliLog::SetModuleDebugLevel("FMD", oldDbg);
+  }
+
   TH1* hist;
   Int_t idx = 0;
   for(UShort_t d = 1; d<=3; d++) {
     UShort_t nR = (d == 1 ? 1 : 2); 
     for(UShort_t q = 0; q < nR; q++) {
       Char_t r = (q == 1 ? 'O' : 'I');
-      hist = MakeADCHist(d, r, -1);
+      hist = (fUseReconstructor ? 
+             MakeELossHist(d, r, -1) : 
+             MakeADCHist(d, r, -1));
 
       Int_t index1 = GetHalfringIndex(d, r, 0, 1);
       idx          = TMath::Max(index1, idx);
       Add2RawsList(hist, index1, !expert, image, !saveCorr);
       //AliInfo(Form("Adding %30s to raw list @ %2d", hist->GetName(), index1));
       
+      // If we're using the reconstructor, do not make expert histograms 
+      if (fUseReconstructor) continue;
+      
       for(UShort_t b = 0; b <= 1; b++) {
        //Hexadecimal board numbers 0x0, 0x1, 0x10, 0x11;
        UShort_t board = (q == 1 ? 0 : 1) + b*16;
@@ -374,40 +407,65 @@ void AliFMDQADataMakerRec::MakeRaws(AliRawReader* rawReader)
   // Parameters:
   //    rawReader Raw reader
   //
   AliFMDRawReader fmdReader(rawReader,0);
-  
-  if (fDigitsArray) 
-    fDigitsArray->Clear();
-  else 
-    fDigitsArray = new TClonesArray("AliFMDDigit", 1000);
-
+  if (fDigitsArray) fDigitsArray->Clear();
+  else              fDigitsArray = new TClonesArray("AliFMDDigit", 1000);
+    
   TClonesArray* digitsAddress = fDigitsArray;
-  
+    
   rawReader->Reset();
-               
+    
   digitsAddress->Clear();
   fmdReader.ReadAdcs(digitsAddress);
-  for(Int_t i=0;i<digitsAddress->GetEntriesFast();i++) {
-    //Raw ADC counts
-    AliFMDDigit*      digit = static_cast<AliFMDDigit*>(digitsAddress->At(i));
-    UShort_t          det   = digit->Detector();
-    Char_t            ring  = digit->Ring();
-    UShort_t          sec   = digit->Sector();
-    // UShort_t strip = digit->Strip();
-    AliFMDParameters* pars  = AliFMDParameters::Instance();
-    Short_t           board = pars->GetAltroMap()->Sector2Board(ring, sec);
-    
-    Int_t index1 = GetHalfringIndex(det, ring, 0, 1);
-    FillRawsData(index1,digit->Counts());
-    Int_t index2 = GetHalfringIndex(det, ring, board/16,0);
-    FillRawsData(index2,digit->Counts());
-    
-  }
   //
   FillRawsData(1,1, fmdReader.GetNErrors(0));
   FillRawsData(1,2, fmdReader.GetNErrors(1));
   FillRawsData(1,3, fmdReader.GetNErrors(2));
+
+  if (fUseReconstructor) { 
+    AliESDFMD* fmd = fReconstructor->GetESDObject();
+    fmd->Clear();
+    
+    // AliLog::SetModuleDebugLevel("FMD", 15);
+    fReconstructor->ProcessDigits(digitsAddress, fmdReader);
+
+    if (!fmd) AliFatal("No ESD object from reconstructor");
+
+    for(UShort_t det=1;det<=3;det++) {
+      UShort_t nrng = (det == 1 ? 1 : 2);
+      for (UShort_t ir = 0; ir < nrng; ir++) {
+       Char_t   ring = (ir == 0 ? 'I' : 'O');
+       UShort_t nsec = (ir == 0 ? 20  : 40);
+       UShort_t nstr = (ir == 0 ? 512 : 256);
+       for(UShort_t sec =0; sec < nsec;  sec++)  {
+         for(UShort_t strip = 0; strip < nstr; strip++) {
+           Float_t mult = fmd->Multiplicity(det,ring,sec,strip);
+           if(mult == AliESDFMD::kInvalidMult) continue;
+           
+           Int_t index1 = GetHalfringIndex(det, ring, 0, 1);
+           FillRawsData(index1,mult);
+         }
+       }
+      }
+    }
+  }
+  else {
+    for(Int_t i=0;i<digitsAddress->GetEntriesFast();i++) {
+      //Raw ADC counts
+      AliFMDDigit*      digit = static_cast<AliFMDDigit*>(digitsAddress->At(i));
+      UShort_t          det   = digit->Detector();
+      Char_t            ring  = digit->Ring();
+      UShort_t          sec   = digit->Sector();
+      // UShort_t strip = digit->Strip();
+      AliFMDParameters* pars  = AliFMDParameters::Instance();
+      Short_t           board = pars->GetAltroMap()->Sector2Board(ring, sec);
+      
+      Int_t index1 = GetHalfringIndex(det, ring, 0, 1);
+      FillRawsData(index1,digit->Counts());
+      Int_t index2 = GetHalfringIndex(det, ring, board/16,0);
+      FillRawsData(index2,digit->Counts());
+    }
+  }
   //
   IncEvCountCycleRaws();
   IncEvCountTotalRaws();
index cf3ec06..b8813f0 100644 (file)
@@ -10,7 +10,7 @@
 class TH1F; 
 class TH1I; 
 class TList; 
-
+class AliFMDReconstructor;
 
 //_____________________________________________________________________
 // This class implements the AliQADataMakerRec for the FMD. Some
@@ -116,6 +116,8 @@ private:
   Int_t GetHalfringIndex(UShort_t det, Char_t ring, 
                         UShort_t board, UShort_t monitor = 0) const;
   TClonesArray fRecPointsArray; // Rec points
+  AliFMDReconstructor* fReconstructor;
+  Bool_t               fUseReconstructor;
 
   ClassDef(AliFMDQADataMakerRec,0)  // description 
 };
index f9f3f7f..be523f6 100644 (file)
@@ -314,16 +314,9 @@ AliFMDReconstructor::MarkDeadChannels(AliESDFMD* esd) const
 
 //____________________________________________________________________
 void 
-AliFMDReconstructor::Reconstruct(AliRawReader* reader, TTree*) const
+AliFMDReconstructor::Reconstruct(AliFMDRawReader& rawReader) const
 {
-  // Reconstruct directly from raw data (no intermediate output on
-  // digit tree or rec point tree).  
-  // 
-  // Parameters: 
-  //   reader  Raw event reader 
-  //   ctree    Not used - 'cluster tree' to store rec-points in. 
-  AliFMDDebug(1, ("Reconstructing from raw reader"));
-  AliFMDRawReader rawReader(reader, 0);
+  AliFMDDebug(1, ("Reconstructing from FMD raw reader"));
   fBad.Reset(false);
   UShort_t det, sec, str, fac;
   Short_t  adc, oldDet = -1;
@@ -340,7 +333,22 @@ AliFMDReconstructor::Reconstruct(AliRawReader* reader, TTree*) const
     ProcessSignal(det, rng, sec, str, adc);
   }
   UseRecoParam(kFALSE);
-  
+
+}
+
+//____________________________________________________________________
+void 
+AliFMDReconstructor::Reconstruct(AliRawReader* reader, TTree*) const
+{
+  // Reconstruct directly from raw data (no intermediate output on
+  // digit tree or rec point tree).  
+  // 
+  // Parameters: 
+  //   reader  Raw event reader 
+  //   ctree    Not used - 'cluster tree' to store rec-points in. 
+  AliFMDDebug(1, ("Reconstructing from raw reader"));
+  AliFMDRawReader rawReader(reader, 0);
+  Reconstruct(rawReader);
 }
 
 //____________________________________________________________________
@@ -427,6 +435,28 @@ AliFMDReconstructor::Reconstruct(TTree* digitsTree,
 
 //____________________________________________________________________
 void
+AliFMDReconstructor::ProcessDigits(TClonesArray* digits, 
+                                  const AliFMDRawReader& rawRead) const
+{
+  // For each digit, find the pseudo rapdity, azimuthal angle, and
+  // number of corrected ADC counts, and pass it on to the algorithms
+  // used. 
+  // 
+  // Parameters: 
+  //    digits Array of digits
+  // 
+  AliFMDAltroMapping* map = AliFMDParameters::Instance()->GetAltroMap();
+  for (size_t i = 1; i <= 3; i++) { 
+    fZS[i-1]       = rawRead.IsZeroSuppressed(map->Detector2DDL(i));
+    fZSFactor[i-1] = rawRead.NoiseFactor(map->Detector2DDL(i));
+  }
+  UseRecoParam(kTRUE);
+  ProcessDigits(digits);
+  UseRecoParam(kFALSE);
+}
+
+//____________________________________________________________________
+void
 AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
 {
   // For each digit, find the pseudo rapdity, azimuthal angle, and
index 90bfea4..2e94293 100644 (file)
@@ -32,6 +32,7 @@ class TTree;
 class TClonesArray;
 class AliFMDDigit;
 class AliRawReader;
+class AliFMDRawReader;
 class AliESDEvent;
 class AliESDFMD;
 class AliFMDRecoParam;
@@ -103,6 +104,7 @@ public:
    *       TClonesArray of AliFMDDigits
    */
   virtual void   Reconstruct(AliRawReader *, TTree*) const;
+  virtual void Reconstruct(AliFMDRawReader& reader) const;
   /** 
    * Put in the ESD data, the FMD ESD data.  The object created by
    * the Reconstruct member function is copied to the ESD object. 
@@ -119,7 +121,12 @@ public:
    */
   virtual void   FillESD(AliRawReader*, TTree* clusterTree, 
                         AliESDEvent* esd) const;
-
+  /** 
+   * Return the filled FMD ESD object
+   * 
+   * @return FMD ESD object
+   */
+  AliESDFMD* GetESDObject() const { return fESDObj; }
   /** 
    * Create SDigits from raw data
    * 
@@ -161,6 +168,17 @@ public:
    * @param use If true, make the diagnostics file 
    */
   void SetDiagnose(Bool_t use=kTRUE) { fDiagnostics = use; }
+  /** 
+   * Process AliFMDDigit objects in @a digits.  For each digit, find
+   * the psuedo-rapidity @f$ \eta@f$, azimuthal angle @f$ \varphi@f$,
+   * energy deposited @f$ E@f$, and psuedo-inclusive multiplicity @f$
+   * M@f$.
+   * 
+   * @param digits  Array of digits. 
+   * @param rawRead Raw reader used 
+   */
+  virtual void ProcessDigits(TClonesArray* digits,
+                            const AliFMDRawReader& rawRead) const;
     
 protected:
   /** 
@@ -186,6 +204,14 @@ protected:
    */
   virtual void GetVertex(AliESDEvent* esd) const;
   /** 
+   * Set-up reconstructor to use values from reconstruction
+   * parameters, if present, for this event.   If the argument @a set
+   * is @c false, then restore preset values. 
+   * 
+   * @param set 
+   */  
+  virtual void UseRecoParam(Bool_t set=kTRUE) const;
+  /** 
    * Process AliFMDDigit objects in @a digits.  For each digit, find
    * the psuedo-rapidity @f$ \eta@f$, azimuthal angle @f$ \varphi@f$,
    * energy deposited @f$ E@f$, and psuedo-inclusive multiplicity @f$
@@ -193,7 +219,7 @@ protected:
    * 
    * @param digits Array of digits. 
    */
-  virtual void     ProcessDigits(TClonesArray* digits) const;
+  virtual void ProcessDigits(TClonesArray* digits) const;
   /** 
    * Process a single digit 
    * 
@@ -374,14 +400,6 @@ protected:
   void MarkDeadChannels(AliESDFMD* esd) const;
 
   /** 
-   * Set-up reconstructor to use values from reconstruction
-   * parameters, if present, for this event.   If the argument @a set
-   * is @c false, then restore preset values. 
-   * 
-   * @param set 
-   */  
-  virtual void UseRecoParam(Bool_t set=kTRUE) const;
-  /** 
    * Utility member function to get the reconstruction parameters for 
    * this event
    *