From 1306ba553734c54a4beabd3f8e6ca0d888403dbc Mon Sep 17 00:00:00 2001 From: cholm Date: Tue, 17 Jul 2012 12:08:27 +0000 Subject: [PATCH] Updates for better DQM. 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 void FMDshifterQuality(TObjArray* objs, map* 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 | 678 +++++++++++++++++++++++++++++++++-- FMD/AliFMDQAChecker.h | 178 ++++++++- FMD/AliFMDQADataMakerRec.cxx | 128 +++++-- FMD/AliFMDQADataMakerRec.h | 4 +- FMD/AliFMDReconstructor.cxx | 50 ++- FMD/AliFMDReconstructor.h | 38 +- 6 files changed, 976 insertions(+), 100 deletions(-) diff --git a/FMD/AliFMDQAChecker.cxx b/FMD/AliFMDQAChecker.cxx index 4453d11d343..5a6b7c2c258 100644 --- a/FMD/AliFMDQAChecker.cxx +++ b/FMD/AliFMDQAChecker.cxx @@ -27,6 +27,7 @@ // --- ROOT system --- #include #include +#include #include #include #include @@ -37,6 +38,9 @@ #include #include #include +#include +#include +#include // --- AliRoot header files --- #include "AliLog.h" @@ -44,14 +48,344 @@ #include "AliQAChecker.h" #include "AliFMDQAChecker.h" #include "AliRecoParam.h" +#include +#include +#include +#include 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(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(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: + // - Nucl.Instrum.Meth.B1:16 + // - Phys.Rev.A28:615 + // - ROOT implementation + // + // 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* p = static_cast*>(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(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(fmdObj); + Int_t nThr = qaThr->GetSize(); + for (Int_t i = 0; i < nThr; i++) { + TObject* thr = qaThr->GetThreshold(i); + if (!thr) continue; + + TParameter* d = dynamic_cast*>(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(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(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(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 } } diff --git a/FMD/AliFMDQAChecker.h b/FMD/AliFMDQAChecker.h index d65e212b4e0..7f11b1e5709 100644 --- a/FMD/AliFMDQAChecker.h +++ b/FMD/AliFMDQAChecker.h @@ -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? }; diff --git a/FMD/AliFMDQADataMakerRec.cxx b/FMD/AliFMDQADataMakerRec.cxx index 4959baccda2..c5b1b59b440 100644 --- a/FMD/AliFMDQADataMakerRec.cxx +++ b/FMD/AliFMDQADataMakerRec.cxx @@ -19,10 +19,12 @@ #include #include #include +#include // --- 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;iGetEntriesFast();i++) { - //Raw ADC counts - AliFMDDigit* digit = static_cast(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;iGetEntriesFast();i++) { + //Raw ADC counts + AliFMDDigit* digit = static_cast(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(); diff --git a/FMD/AliFMDQADataMakerRec.h b/FMD/AliFMDQADataMakerRec.h index cf3ec06d814..b8813f0b1e4 100644 --- a/FMD/AliFMDQADataMakerRec.h +++ b/FMD/AliFMDQADataMakerRec.h @@ -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 }; diff --git a/FMD/AliFMDReconstructor.cxx b/FMD/AliFMDReconstructor.cxx index f9f3f7f7730..be523f62b5d 100644 --- a/FMD/AliFMDReconstructor.cxx +++ b/FMD/AliFMDReconstructor.cxx @@ -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); } //____________________________________________________________________ @@ -425,6 +433,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 diff --git a/FMD/AliFMDReconstructor.h b/FMD/AliFMDReconstructor.h index 90bfea431c3..2e9429390c6 100644 --- a/FMD/AliFMDReconstructor.h +++ b/FMD/AliFMDReconstructor.h @@ -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: /** @@ -185,6 +203,14 @@ protected: * @param esd ESD structure to get Vz from */ 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$, @@ -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 * @@ -373,14 +399,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 -- 2.43.0