diff --git a/core/qa/CMakeLists.txt b/core/qa/CMakeLists.txt
index 5c4eeb7907a0658a8c988b3af9be4f7bd9b02956..e62f00327e0fa61ad80426c29374539d5cd82b15 100644
--- a/core/qa/CMakeLists.txt
+++ b/core/qa/CMakeLists.txt
@@ -11,6 +11,7 @@ set(SRCS
   CbmQaTable.cxx
   CbmQaTask.cxx
   CbmQaIO.cxx
+  CbmQaUtil.cxx
   checker/CbmQaCheckerCore.cxx
   checker/CbmQaCheckerFileHandler.cxx
   checker/CbmQaCheckerHist1DHandler.cxx
@@ -52,6 +53,7 @@ Install(FILES
         CbmQaHist.h
         CbmQaEff.h
         CbmQaPie.h
+        CbmQaUtil.h
         CbmQaConstants.h
         CbmQaCmpDrawer.h
         checker/CbmQaCheckerTypedefs.h
diff --git a/core/qa/CbmQaIO.cxx b/core/qa/CbmQaIO.cxx
index e09b028a0a256f886a1f73ff399c6437f5ecc3d2..56873a5d531bca0087f00e9f572d1c2888e472f6 100644
--- a/core/qa/CbmQaIO.cxx
+++ b/core/qa/CbmQaIO.cxx
@@ -13,26 +13,19 @@
 
 // ---------------------------------------------------------------------------------------------------------------------
 //
-CbmQaIO::CbmQaIO(const char* prefixName, std::shared_ptr<ObjList_t> pObjList)
-  : fsPrefix(prefixName)
-  , fpvObjList(pObjList)
+CbmQaIO::CbmQaIO(TString prefixName, std::shared_ptr<ObjList_t> pObjList) : fsPrefix(prefixName), fpvObjList(pObjList)
 {
+  if (!fsPrefix.IsNull()) { fsPrefix += "_"; }
   if (!fpvObjList.get()) { fpvObjList = std::make_shared<ObjList_t>(); }
 }
 
 // ---------------------------------------------------------------------------------------------------------------------
 //
-CbmQaIO::~CbmQaIO()
-{
-}
+CbmQaIO::~CbmQaIO() {}
 
 // ---------------------------------------------------------------------------------------------------------------------
 //
-void CbmQaIO::SetHistoProperties(TH1* pHist)
-{
-  pHist->SetStats(true);
-  pHist->Sumw2();
-}
+void CbmQaIO::SetHistoProperties(TH1* pHist) { pHist->SetStats(true); }
 
 // ---------------------------------------------------------------------------------------------------------------------
 //
@@ -43,6 +36,6 @@ void CbmQaIO::WriteToFile(TFile* pOutFile) const
   for (const auto& [pObject, sPath] : (*fpvObjList)) {
     if (!pOutFile->GetDirectory(sPath)) { pOutFile->mkdir(sPath); }
     pOutFile->cd(sPath);
-    pObject->Write();
+    if (pObject) { pObject->Write(); }
   }
 }
diff --git a/core/qa/CbmQaIO.h b/core/qa/CbmQaIO.h
index f2f235b2b38e32f55910018db3b7fb79849c7a68..199e7603f9eb3a310365e578b6c46092cd64c2bc 100644
--- a/core/qa/CbmQaIO.h
+++ b/core/qa/CbmQaIO.h
@@ -10,6 +10,8 @@
 #ifndef CbmQaIO_h
 #define CbmQaIO_h 1
 
+#include "CbmQaCanvas.h"
+#include "CbmQaEff.h"
 #include "CbmQaTable.h"
 
 #include "Logger.h"
@@ -29,6 +31,7 @@
 
 #include <vector>
 
+
 class TFile;
 
 /// @brief  ROOT object IO interface for QA
@@ -40,14 +43,14 @@ public:
 
   enum class EStoringMode
   {
-    kSAMEDIR,  ///< Objects will be stored to root directory
-    kSUBDIR    ///< Objects will be stored to corresponding subdirectory
+    kSAMEDIR,  ///< Objects of different type will be stored to root directory
+    kSUBDIR    ///< Objects of different type will be stored in different subdirectories like histograms/ canvases/
   };
 
   /// @brief Constructor
   /// @param prefixName    Name of the unique prefix
-  /// @param pParentFolder Pointer to parent folder
-  CbmQaIO(const char* prefixName, std::shared_ptr<ObjList_t> pObjList = nullptr);
+  /// @param pObjList      List of already created objects
+  CbmQaIO(TString prefixName, std::shared_ptr<ObjList_t> pObjList = nullptr);
 
   /// @brief Destructor
   virtual ~CbmQaIO();
@@ -64,42 +67,46 @@ public:
   /// @brief Move assignment operator
   CbmQaIO& operator=(CbmQaIO&&) = delete;
 
+  /// Set storing mode
+  void SetStoringMode(EStoringMode val) { fStoringMode = val; }
+
+  /// @brief Sets a common root path to the objects in the output file
+  /// @param path  A path to the object
+  void SetRootFolderName(TString path) { fsRootFolderName = path; }
+
+  /// @brief Creates a ROOT object
+  /// @param name    Name of the object including the path
+  /// @param args... Arguments passed to the object constructor
+  template<typename T, typename... Args>
+  T* MakeQaObject(TString name, Args... args);
+
   /// @brief Creates, initializes and registers a canvas
   /// @tparam  T     Type of the canvas: TCanvas or CbmQaCanvas (\note T should not be a pointer)
   /// @param   name  Name of canvas
   /// @param   args  The rest of the arguments, which will be passed to the canvas constructor
   template<typename T, typename... Args>
-  T* MakeCanvas(const char* name, Args... args);
+  T* MakeCanvas(TString name, Args... args);
 
   /// @brief Creates, initializes and registers an efficiency
   /// @tparam  T     Type of the canvas: either TProfile, TProfile1D or TEfficiency (\note T should not be a pointer)
   /// @param   name  Name of canvas
   /// @param   args  The rest of the arguments, which will be passed to the efficiency constructor
   template<typename T, typename... Args>
-  T* MakeEfficiency(const char* name, Args... args);
+  T* MakeEfficiency(TString name, Args... args);
 
   /// @brief Creates, initializes and registers a histogram
   /// @tparam  T     Type of the histogram (\note T should not be a pointer)
   /// @param   name  Name of histogram
   /// @param   args  The rest of the arguments, which will be passed to the histogram constructor
   template<typename T, typename... Args>
-  T* MakeHisto(const char* name, Args... args);
+  T* MakeHisto(TString name, Args... args);
 
   /// @brief Creates, initializes and registers a table
   /// @param  name  Name of the table
   /// @param  args  The rest of the arguments, which will be passed to the table constructor
   template<typename... Args>
-  CbmQaTable* MakeTable(const char* name, Args... args);
+  CbmQaTable* MakeTable(TString name, Args... args);
 
-  /// @brief Creates a ROOT object
-  /// @param name    Name of the object
-  /// @param args... Arguments passed to the object constructor
-  template<typename T, typename... Args>
-  T* MakeObject(const char* name, Args... args);
-
-  /// @brief Sets a common root path to the objects in the output file
-  /// @param path  A path to the object
-  void SetRootFolderName(const char* path) { fsRootFolderName = path; }
 
 protected:
   /// @brief Applies properties on the histogram created with MakeHisto funciton
@@ -121,75 +128,81 @@ protected:
 // ---------------------------------------------------------------------------------------------------------------------
 //
 template<typename T, typename... Args>
-T* CbmQaIO::MakeCanvas(const char* nameBase, Args... args)
+T* CbmQaIO::MakeQaObject(TString sObjName, Args... args)
 {
-  TString sFullName = EStoringMode::kSUBDIR == fStoringMode ? Form("canvases/%s", nameBase) : nameBase;
-  auto* pCanv       = MakeObject<T>(sFullName, args...);
-  return pCanv;
+  // Resolve directory name and object name
+  auto iLastSlashPos = static_cast<int>(sObjName.Last('/'));
+  TString sDirName   = sObjName(0, iLastSlashPos);
+  sObjName           = sObjName(iLastSlashPos + 1, sObjName.Length() - iLastSlashPos);
+
+  // Add unique prefix to the object name
+  sObjName = fsPrefix + sObjName;
+
+  // Check, if an object with the same name already exists
+  // TODO: check if this makes sense and doesn't crash the memory
+  //if (gROOT->FindObject(sObjName)) {
+  //  LOG(warn) << fsRootFolderName << ": A previously created object \"" << sObjName << "\" will be deleted";
+  //  auto* pObj = static_cast<T*>(gROOT->FindObject(sObjName));
+  //  delete pObj;
+  //}
+
+  // Create a new object
+  T* pObj = new T(sObjName, args...);
+
+  // specific stuff for different kind of objects
+
+  TH1* pH1 = dynamic_cast<TH1*>(pObj);
+
+  if (pH1) {               // any kind of histogram
+    pH1->SetDirectory(0);  // don't let ROOT own the histogram. The user should take care on deleting it.
+    SetHistoProperties(pH1);
+  }
+
+  TString sPrefix = "";
+  if (EStoringMode::kSUBDIR == fStoringMode) {
+    if (dynamic_cast<CbmQaTable*>(pObj) || dynamic_cast<CbmQaCanvas*>(pObj)) { sPrefix = "summary/"; }
+  }
+
+  // Add parent directory to the folder
+  if (fsRootFolderName.Length() != 0) { sDirName = fsRootFolderName + "/" + sPrefix + sDirName; }
+
+  // Register the object into an object list
+  fpvObjList->push_back(std::make_pair(static_cast<TObject*>(pObj), sDirName));
+  return pObj;
 }
 
+
 // ---------------------------------------------------------------------------------------------------------------------
 //
 template<typename T, typename... Args>
-T* CbmQaIO::MakeEfficiency(const char* nameBase, Args... args)
+T* CbmQaIO::MakeCanvas(TString nameBase, Args... args)
 {
-  TString sFullName = EStoringMode::kSUBDIR == fStoringMode ? Form("efficiencies/%s", nameBase) : nameBase;
-  auto* pEff        = MakeObject<T>(sFullName, args...);
-  return pEff;
+  return MakeQaObject<T>(nameBase, args...);
 }
 
 // ---------------------------------------------------------------------------------------------------------------------
 //
 template<typename T, typename... Args>
-T* CbmQaIO::MakeHisto(const char* nameBase, Args... args)
+T* CbmQaIO::MakeEfficiency(TString nameBase, Args... args)
 {
-  TString sFullName = EStoringMode::kSUBDIR == fStoringMode ? Form("histograms/%s", nameBase) : nameBase;
-  auto* pHist       = MakeObject<T>(sFullName, args...);
-  pHist->SetDirectory(0);
-  SetHistoProperties(pHist);
-  return pHist;
+  return MakeQaObject<T>(nameBase, args...);
 }
 
 // ---------------------------------------------------------------------------------------------------------------------
 //
 template<typename T, typename... Args>
-T* CbmQaIO::MakeObject(const char* name, Args... args)
+T* CbmQaIO::MakeHisto(TString nameBase, Args... args)
 {
-  // Resolve directory name and object name
-  TString sObjName   = name;
-  auto iLastSlashPos = static_cast<int>(sObjName.Last('/'));
-  TString sDirName   = sObjName(0, iLastSlashPos);
-  sObjName           = sObjName(iLastSlashPos + 1, sObjName.Length() - iLastSlashPos);
-
-  // Add unique prefix to the object name
-  sObjName = fsPrefix + "_" + sObjName;
-
-  // Add parent directory to the folder
-  if (fsRootFolderName.Length() != 0) { sDirName = fsRootFolderName + "/" + sDirName; }
-
-  // Check, if an object with the same name already exists
-  if (gROOT->FindObject(sObjName)) {
-    LOG(warn) << fsRootFolderName << ": A previously created object \"" << sObjName << "\" will be deleted";
-    auto* pObj = static_cast<T*>(gROOT->FindObject(sObjName));
-    delete pObj;
-  }
-
-  // Create a new object
-  T* pObj = new T(sObjName, args...);
-
-  // Register the object into an object list
-  fpvObjList->push_back(std::make_pair(static_cast<TObject*>(pObj), sDirName));
-  return pObj;
+  return MakeQaObject<T>(nameBase, args...);
 }
 
 // ---------------------------------------------------------------------------------------------------------------------
 //
 template<typename... Args>
-CbmQaTable* CbmQaIO::MakeTable(const char* nameBase, Args... args)
+CbmQaTable* CbmQaIO::MakeTable(TString nameBase, Args... args)
 {
-  TString sFullName = EStoringMode::kSUBDIR == fStoringMode ? Form("tables/%s", nameBase) : nameBase;
-  auto* pTable      = MakeObject<CbmQaTable>(sFullName, args...);
-  return pTable;
+  return MakeQaObject<CbmQaTable>(nameBase, args...);
 }
 
+
 #endif  // CbmQaIO_h
diff --git a/core/qa/CbmQaTask.cxx b/core/qa/CbmQaTask.cxx
index 53effef64abb1c8aeb85396f887c81f295150efe..2782247f6a260d610c70a5a761590363519cf5e5 100644
--- a/core/qa/CbmQaTask.cxx
+++ b/core/qa/CbmQaTask.cxx
@@ -10,8 +10,6 @@
 
 #include "CbmQaTask.h"
 
-#include "CbmQaCanvas.h"
-
 #include "FairRootFileSink.h"
 #include "FairRootManager.h"
 #include "FairRunAna.h"
@@ -19,6 +17,7 @@
 #include "TAxis.h"
 #include "TCanvas.h"
 #include "TF1.h"
+#include "TString.h"
 
 #include <array>
 
@@ -133,3 +132,25 @@ void CbmQaTask::DeInitBase()
   // De-initialize particular QA-task implementation
   DeInit();
 }
+
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+bool CbmQaTask::CheckRange(TH1* h, double meanMax, double rmsMin, double rmsMax)
+{
+  // Checks ranges for mean and standard deviation
+  // \return  False, if variable exits the range
+
+  double mean    = h->GetMean();
+  double meanErr = h->GetMeanError();
+
+  double rms    = h->GetStdDev();
+  double rmsErr = h->GetStdDevError();
+
+  bool res = true;
+  if (h->GetEntries() >= 10) {  // ask for some statistics, otherwise errors are not properly estimated
+    res = CheckRange(TString("Mean for ") + h->GetTitle(), mean, meanErr, -meanMax, meanMax) && res;
+    res = CheckRange(TString("RMS for ") + h->GetTitle(), rms, rmsErr, rmsMin, rmsMax) && res;
+  }
+  return res;
+}
diff --git a/core/qa/CbmQaTask.h b/core/qa/CbmQaTask.h
index e17a8b9fddd23cdd8b5cc16b11d14ebe06637d33..e79ce19fcaa995e961d96ab6d8eb3e922bc95a4e 100644
--- a/core/qa/CbmQaTask.h
+++ b/core/qa/CbmQaTask.h
@@ -31,6 +31,7 @@
 #include <algorithm>
 #include <regex>
 #include <string_view>
+#include <tuple>
 #include <type_traits>
 
 #include <yaml-cpp/yaml.h>
@@ -81,8 +82,6 @@ public:
   /// parameters: title, number of bins and axis ranges.
   void SetConfigName(const char* filename) { fsConfigName = filename; }
 
-  ClassDef(CbmQaTask, 0);
-
 protected:
   // *****************************************************
   // ** Functions accessible inside the derived classes **
@@ -124,19 +123,34 @@ protected:
   /// \param hi    Upper limit of the variable
   /// \return  False, if variable exits the range
   template<typename T>
-  bool CheckRange(std::string_view name, const T& var, const T& lo, const T& hi) const;
+  bool CheckRange(std::string_view name, T var, T lo, T hi) const;
 
+  /// Checks range of variable
+  /// \param name    Name of the variable
+  /// \param var     Variable to check
+  /// \param varErr  Standard Error of the variable
+  /// \param lo      Lower limit of the variable
+  /// \param hi      Upper limit of the variable
+  /// \return        False, if variable exits the range
+  template<typename T>
+  bool CheckRange(std::string_view name, T var, T varErr, T lo, T hi) const;
+
+  /// Checks ranges for mean and standard deviation
+  /// \return  False, if variable exits the range
+  bool CheckRange(TH1* h, double meanMax, double rmsMin, double rmsMax);
 
 private:
   /// @brief De-initializes this task
   void DeInitBase();
 
-  bool fbUseMC     = false;  ///< Flag, if MC is used
+  bool fbUseMC = false;  ///< Flag, if MC is used
 
   TParameter<int> fNofEvents {"nEvents", 0};  ///< Number of processed events
 
   TString fsConfigName      = "";       ///< Name of YAML configuration file
   YAML::Node* fpCurrentNode = nullptr;  ///< Pointer to current node of the configuration file
+
+  ClassDef(CbmQaTask, 0);
 };
 
 
@@ -147,17 +161,28 @@ private:
 // ---------------------------------------------------------------------------------------------------------------------
 //
 template<typename T>
-bool CbmQaTask::CheckRange(std::string_view name, const T& var, const T& lo, const T& hi) const
+bool CbmQaTask::CheckRange(std::string_view name, T var, T lo, T hi) const
+{
+  return CheckRange<T>(name, var, 0., lo, hi);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T>
+bool CbmQaTask::CheckRange(std::string_view name, T var, T varErr, T lo, T hi) const
 {
-  if (std::clamp(var, lo, hi) == lo) {
-    LOG(error) << fName << ": " << name << " is found to be under the lower limit (" << var << " < " << lo << ')';
-    return false;
+  bool ret = true;
+  if (var + 3.5 * varErr < lo) {
+    LOG(error) << fName << ": " << name << " is found to be under the lower limit (" << var << " +-" << varErr << " < "
+               << lo << ')';
+    ret = false;
   }
-  if (std::clamp(var, lo, hi) == hi) {
-    LOG(error) << fName << ": " << name << " is found to be above the upper limit (" << var << " > " << hi << ')';
-    return false;
+  if (var - 3.5 * varErr > hi) {
+    LOG(error) << fName << ": " << name << " is found to be above the upper limit (" << var << " +-" << varErr << " > "
+               << hi << ')';
+    ret = false;
   }
-  return true;
+  return ret;
 }
 
 #endif  // CbmQaTask_h
diff --git a/core/qa/CbmQaUtil.cxx b/core/qa/CbmQaUtil.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..b28128eb079e0bfc71ae32a036362bdc101dff8a
--- /dev/null
+++ b/core/qa/CbmQaUtil.cxx
@@ -0,0 +1,114 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergey Gorbunov [committer] */
+
+/// \file   CbmQaUtil.cxx
+/// \brief  Useful utilities for CBM QA tasks
+/// \author S.Gorbunov
+/// \data   31.07.2023
+
+
+#include "CbmQaUtil.h"
+
+#include "CbmQaCanvas.h"
+
+#include "TAxis.h"
+#include "TCanvas.h"
+#include "TF1.h"
+#include "TH1.h"
+#include "TPaveStats.h"
+#include "TVirtualPad.h"
+
+namespace cbm::qa::util
+{
+
+  // ---------------------------------------------------------------------------------------------------------------------
+  //
+  std::tuple<double, double> FitKaniadakisGaussian(TH1* pHist)
+  {
+    //
+    // Fit function - Kaniadakis Gaussian Distribution:
+    // https://en.wikipedia.org/wiki/Kaniadakis_Gaussian_distribution
+    //
+    // where exp_k(x) is calculated via arcsinh():
+    // https://en.wikipedia.org/wiki/Kaniadakis_statistics
+    //
+    // parameters:
+    //
+    // [0] - mean
+    // [1] - Std. Dev ( calculated after the fit )
+    // [2] - peak
+    // [3] - hwhm - half width at half maximum, scaled.  ( == standard deviation for a gaussian )
+    // [4] - k - Kaniadakis parameter, k in [0.,1.]. ( == 0.0 for a gaussian. The formula below breaks when k is exactly 0.)
+    //
+    // returns mean and std.dev
+    //
+
+    double xMin = pHist->GetXaxis()->GetXmin();
+    double xMax = pHist->GetXaxis()->GetXmax();
+
+    TF1 fit("FitKaniadakisGaussian", "[2] * TMath::Exp(TMath::ASinH(-0.5*[4]*((x-[0])/[3])**2)/[4]) + 0.0*[1]", xMin,
+            xMax, "NL");
+
+    fit.SetParName(0, "fit_Mean");
+    fit.SetParName(1, "fit_StdDev");
+    fit.SetParName(2, "fit_Peak");
+    fit.SetParName(3, "fit_Hwhm");
+    fit.SetParName(4, "fit_k");
+
+    // set reasonable initial values
+
+    double mean = pHist->GetMean();
+    double peak = pHist->GetBinContent(pHist->GetMaximumBin());
+    double hwhm = pHist->GetStdDev();
+
+    fit.SetParameters(mean, hwhm, peak, hwhm, .001);
+
+    // limit parameter k
+
+    fit.SetParLimits(4, 0.00001, 5.);
+
+    // fit in a quite mode
+
+    TVirtualPad* padsav = gPad;
+    CbmQaCanvas::GetDummyCanvas().cd();
+    pHist->SetStats(0);  // remove the old stat window
+    pHist->SetStats(1);  // set the flag to create the stat window during Draw()
+    //pHist->Draw();
+    //CbmQaCanvas::GetDummyCanvas().Update();
+
+    //pHist->Sumw2();
+    if (pHist->GetEntries() > 0) { pHist->Fit(&fit, "Q"); }
+    pHist->Draw();
+    CbmQaCanvas::GetDummyCanvas().Update();
+
+    TF1* f = pHist->GetFunction("FitKaniadakisGaussian");
+    if (nullptr != f) {
+
+      // calculate the Std Dev
+
+      f->SetParameter(1, sqrt(f->CentralMoment(2, xMin, xMax)));
+      f->SetParError(1, f->GetParError(3));
+
+      // fix some parameters to prevent them from showing up in the stat window
+
+      f->FixParameter(2, f->GetParameter(2));
+      f->FixParameter(3, f->GetParameter(3));
+      f->FixParameter(4, f->GetParameter(4));
+    }
+
+    TPaveStats* stats = (TPaveStats*) pHist->FindObject("stats");
+    assert(stats);
+
+    stats->SetX1NDC(0.7);
+    stats->SetY1NDC(0.5);
+    stats->SetOptStat(111110);
+    //stats->SetOptStat("eMRuo");
+    stats->SetOptFit(100001);
+
+    CbmQaCanvas::GetDummyCanvas().Clear();
+    if (padsav) padsav->cd();
+    return (nullptr != f) ? std::tuple(f->GetParameter(0), f->GetParameter(1)) : std::tuple(0., -1.);
+  }
+
+}  // namespace cbm::qa::util
diff --git a/core/qa/CbmQaUtil.h b/core/qa/CbmQaUtil.h
new file mode 100644
index 0000000000000000000000000000000000000000..a6a02b015079ac23315fca670b48aa6c083d0ba5
--- /dev/null
+++ b/core/qa/CbmQaUtil.h
@@ -0,0 +1,27 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergey Gorbunov [committer] */
+
+/// \file   CbmQaUtil.h
+/// \brief  Useful utilities for CBM QA tasks
+/// \author S.Gorbunov
+/// \data   31.07.2023
+
+#ifndef CbmQaUtil_h
+#define CbmQaUtil_h 1
+
+#include "TObject.h"
+
+#include <tuple>
+class TH1;
+
+/// namespace cbm::qa::util contains useful utilities for CBM QA tasks
+namespace cbm::qa::util
+{
+
+  /// @brief Fit a histogram with Kaniadakis Gaussian Distribution
+  /// @return  mean and std.dev of the fit
+  std::tuple<double, double> FitKaniadakisGaussian(TH1* pHist);
+}  // namespace cbm::qa::util
+
+#endif  // CbmQaUtil_h
diff --git a/macro/run/run_qa.C b/macro/run/run_qa.C
index b259bee06d6b528fd7d3e76e38cc8c876fb7b81e..3d451b40abf30705b42353a99bdf9a984ac2d322 100644
--- a/macro/run/run_qa.C
+++ b/macro/run/run_qa.C
@@ -240,7 +240,6 @@ void run_qa(TString dataTra = "data/sis100_muon_jpsi_test", TString dataRaw = "d
     //run->AddTask(new CbmStsDigitizeQa()); //opens lots of windows
     run->AddTask(new CbmStsFindTracksQa());
     auto* pCaInputQaSts = new CbmCaInputQaSts(verbose, bUseMC);
-    pCaInputQaSts->SetEfficiencyThrsh(0.5, 0, 100);
     run->AddTask(pCaInputQaSts);
   }
   // ------------------------------------------------------------------------
diff --git a/reco/L1/qa/CbmCaInputQaBase.h b/reco/L1/qa/CbmCaInputQaBase.h
index fc00be730eb0e6c572e6ad6b1bdce427741aec9c..0169b0f68cc2275e83b7ed5df4d124456d194d8a 100644
--- a/reco/L1/qa/CbmCaInputQaBase.h
+++ b/reco/L1/qa/CbmCaInputQaBase.h
@@ -24,6 +24,9 @@ public:
   /// Default destructor
   ~CbmCaInputQaBase() = default;
 
+  /// Sets maximum allowed deviation of pull distribution mean from 0
+  /// \param devThrsh  Deviation threshold (- devThrsh, + devThrsh)
+  void SetPullMeanDeviationThrsh(double devThrsh) { fPullMeanThrsh = devThrsh; }
 
   /// Sets maximum allowed deviation of pull distribution width from unity
   /// \param devThrsh  Deviation threshold (1 - devThrsh, 1 + devThrsh)
@@ -54,7 +57,6 @@ public:
   /// \param mom Minimum absolute value of momentum
   void SetMinMomentum(double mom) { fMinMomentum = mom; }
 
-
 protected:
   // *********************
   // ** Data structures **
@@ -68,7 +70,9 @@ protected:
   };
 
   // ----- QA constants (thresholds, ranges etc.)
+  //TODO: remove check of residuals
   double fResMeanThrsh   = .5;   ///< Maximum allowed deviation of residual mean from zero [sigma]
+  double fPullMeanThrsh  = .05;  ///< Maximum allowed deviation of pull mean from zero
   double fPullWidthThrsh = .1;   ///< Maximum allowed deviation of pull width from unity
   double fEffThrsh       = .5;   ///< Threshold for hit efficiency in the selected range
   double fMaxDiffZStHit  = 1.;   ///< Maximum allowed difference between z-position of hit and station [cm]
diff --git a/reco/L1/qa/CbmCaInputQaSts.cxx b/reco/L1/qa/CbmCaInputQaSts.cxx
index 711e69cf033556a2c718a549a9c524ba3544b2cd..67268bbb26ec175d7a6007b5da19dfd24fc8272a 100644
--- a/reco/L1/qa/CbmCaInputQaSts.cxx
+++ b/reco/L1/qa/CbmCaInputQaSts.cxx
@@ -17,6 +17,7 @@
 #include "CbmQaCanvas.h"
 #include "CbmQaEff.h"
 #include "CbmQaTable.h"
+#include "CbmQaUtil.h"
 #include "CbmStsCluster.h"
 #include "CbmStsHit.h"
 #include "CbmStsPoint.h"
@@ -53,6 +54,10 @@ namespace phys = L1Constants::phys;  // from L1Constants.h
 //
 CbmCaInputQaSts::CbmCaInputQaSts(int verbose, bool isMCUsed) : CbmQaTask("CbmCaInputQaSts", "casts", verbose, isMCUsed)
 {
+  // TODO: remove an enlarged threshold when the pulls are fixed
+  SetPullMeanDeviationThrsh(0.1);
+  SetPullWidthDeviationThrsh(2.);
+  SetEfficiencyThrsh(0.5, 0, 100);
 }
 
 // ---------------------------------------------------------------------------------------------------------------------
@@ -127,7 +132,7 @@ bool CbmCaInputQaSts::Check()
     // Fit efficiency curves
     LOG(info) << "-- Hit efficiency integrated over hit distance from station center";
 
-    auto* pEffTable = MakeTable("eff_table", "Efficiency table", nSt, 2);
+    auto* pEffTable = MakeQaObject<CbmQaTable>("eff_table", "Efficiency table", nSt, 2);
     pEffTable->SetNamesOfCols({"Station ID", "Efficiency"});
     pEffTable->SetColWidth(20);
 
@@ -143,45 +148,23 @@ bool CbmCaInputQaSts::Check()
     // ----- Checks for residuals
     // Check hits for potential biases from central values
 
-    // Function to fit a residual distribution, returns a structure
-    auto FitResidualsAndGetMean = [&](TH1* pHist) {
-      // Fit function - Kaniadakis Gaussian Distribution:
-      // [0] - Scale (correlated with [3])
-      // [1] - Mean
-      // [2] - Sigma (biased from standard deviation)
-      // [3] - Kaniadakis parameter (0 - 2)
-      auto fit =
-        TF1("fitres", "[0]*TMath::Exp(TMath::ASinH(-0.5*[3]*((x-[1])/[2])**2)/[3])/([2]*TMath::Sqrt(2*TMath::Pi()))",
-            -10., 10.);
-      fit.SetParameters(100, 0., 1., .3);
-      fit.SetParLimits(3, 0, 2);
-      pHist->Fit("fitres", "Q");
-      // NOTE: Several fit procedures are made to avoid empty fit results
-      std::array<double, 3> result;
-      result[0] = fit.GetParameter(1);
-      double fitK    = fit.GetParameter(3);
-      double fitHWHM = fit.GetParameter(2) * std::sqrt((1 - std::pow(2, -2 * fitK)) / fitK);
-      result[1]      = -fitHWHM * fResMeanThrsh;
-      result[2]      = +fitHWHM * fResMeanThrsh;
-      return result;
-    };
-
-    auto* pResidualsTable = MakeTable("residuals_mean", "Residual mean values in different stations", nSt, 4);
+    auto* pResidualsTable =
+      MakeQaObject<CbmQaTable>("residuals_mean", "Residual mean values in different stations", nSt, 4);
     pResidualsTable->SetNamesOfCols({"Station ID", "Residual(x) [cm]", "Residual(y) [cm]", "Residual(t) [ns]"});
     pResidualsTable->SetColWidth(20);
 
-    // Fill residuals fit results and check boundaries
-    for (int iSt = 0; iSt < nSt; ++iSt) {
-      auto fitX = FitResidualsAndGetMean(fvph_res_x[iSt]);
-      auto fitY = FitResidualsAndGetMean(fvph_res_y[iSt]);
-      auto fitT = FitResidualsAndGetMean(fvph_res_t[iSt]);
-      res = CheckRange("Mean of x residual [cm] in station " + std::to_string(iSt), fitX[0], fitX[1], fitX[2]) && res;
-      res = CheckRange("Mean of y residual [cm] in station " + std::to_string(iSt), fitY[0], fitY[1], fitY[2]) && res;
-      res = CheckRange("Mean of t residual [ns] in station " + std::to_string(iSt), fitT[0], fitT[1], fitT[2]) && res;
+    // Fit residuals
+    for (int iSt = 0; iSt <= nSt; ++iSt) {
+      cbm::qa::util::FitKaniadakisGaussian(fvph_res_x[iSt]);
+      cbm::qa::util::FitKaniadakisGaussian(fvph_res_y[iSt]);
+      cbm::qa::util::FitKaniadakisGaussian(fvph_res_t[iSt]);
+      cbm::qa::util::FitKaniadakisGaussian(fvph_res_u[iSt]);
+      cbm::qa::util::FitKaniadakisGaussian(fvph_res_v[iSt]);
+
       pResidualsTable->SetCell(iSt, 0, iSt);
-      pResidualsTable->SetCell(iSt, 1, fitX[0]);
-      pResidualsTable->SetCell(iSt, 2, fitX[1]);
-      pResidualsTable->SetCell(iSt, 3, fitX[2]);
+      pResidualsTable->SetCell(iSt, 1, fvph_res_x[iSt]->GetStdDev());
+      pResidualsTable->SetCell(iSt, 2, fvph_res_y[iSt]->GetStdDev());
+      pResidualsTable->SetCell(iSt, 3, fvph_res_t[iSt]->GetStdDev());
     }
 
     LOG(info) << '\n' << pResidualsTable->ToString(8);
@@ -192,48 +175,40 @@ bool CbmCaInputQaSts::Check()
     // when the resolutions are well determined for detecting stations, the pull distributions should have a RMS equal
     // to unity. Here we allow a RMS of the pull distributions to be varied in a predefined range. If the RMS runs out
     // this range, QA task fails.
-    // TODO: Add checks for center
-
-    // Fit pull distributions
-
-    auto FitPullsAndGetSigma = [&](TH1* pHist) {
-      // Fit function - Kaniadakis Gaussian Distribution:
-      // [0] - Scale (correlated with [3])
-      // [1] - Mean
-      // [2] - Sigma (biased from standard deviation)
-      // [3] - Kaniadakis parameter (0 - 2)
-      auto fit = TF1("fitpull", "[0] * TMath::Exp(TMath::ASinH(-0.5*[3]*((x-[1])/[2])**2)/[3])", -10., 10.);
-      fit.SetParameters(100, 0., 1., .3);
-      fit.SetParLimits(3, 0, 2);
-      pHist->Fit("fitpull", "Q");
-      double fitK    = fit.GetParameter(3);
-      double fitHWHM = fit.GetParameter(2) * std::sqrt((1 - std::pow(2, -2 * fitK)) / fitK);
-      return fitHWHM;
-    };
-
-    auto* pPullsTable = MakeTable("pulls_rms", "Pulls std. dev. values in different stations", nSt, 4);
-    pPullsTable->SetNamesOfCols({"Station ID", "Pull(x) sigm", "Pull(y) sigm", "Pull(z) sigm"});
+
+    auto* pPullsTable = MakeQaObject<CbmQaTable>("pulls_rms", "Pulls std. dev. values in different stations", nSt, 4);
+    pPullsTable->SetNamesOfCols({"Station ID", "Pull(x) sigm", "Pull(y) sigm", "Pull(t) sigm"});
     pPullsTable->SetColWidth(20);
 
-    for (int iSt = 0; iSt < nSt; ++iSt) {
-      double rmsPullX = FitPullsAndGetSigma(fvph_pull_x[iSt]);
-      double rmsPullY = FitPullsAndGetSigma(fvph_pull_y[iSt]);
-      double rmsPullT = FitPullsAndGetSigma(fvph_pull_t[iSt]);
+    for (int iSt = 0; iSt < nSt + 1; ++iSt) {
 
-      double rmsPullHi = 1. + fPullWidthThrsh;
-      double rmsPullLo = 1. - fPullWidthThrsh;
+      // Fit pull distributions for nicer representation. Fit results are not used in further checks.
+      cbm::qa::util::FitKaniadakisGaussian(fvph_pull_x[iSt]);
+      cbm::qa::util::FitKaniadakisGaussian(fvph_pull_y[iSt]);
+      cbm::qa::util::FitKaniadakisGaussian(fvph_pull_t[iSt]);
+      cbm::qa::util::FitKaniadakisGaussian(fvph_pull_u[iSt]);
+      cbm::qa::util::FitKaniadakisGaussian(fvph_pull_v[iSt]);
+
+      // Check the pull quality
+      res = CheckRangePull(fvph_pull_x[iSt]) && res;
+      res = CheckRangePull(fvph_pull_y[iSt]) && res;
+      res = CheckRangePull(fvph_pull_t[iSt]) && res;
+      res = CheckRangePull(fvph_pull_u[iSt]) && res;
+      res = CheckRangePull(fvph_pull_v[iSt]) && res;
 
-      res = CheckRange("Rms for x pull in station " + std::to_string(iSt), rmsPullX, rmsPullLo, rmsPullHi) && res;
-      res = CheckRange("Rms for y pull in station " + std::to_string(iSt), rmsPullY, rmsPullLo, rmsPullHi) && res;
-      res = CheckRange("Rms for t pull in station " + std::to_string(iSt), rmsPullT, rmsPullLo, rmsPullHi) && res;
       pPullsTable->SetCell(iSt, 0, iSt);
-      pPullsTable->SetCell(iSt, 1, rmsPullX);
-      pPullsTable->SetCell(iSt, 2, rmsPullY);
-      pPullsTable->SetCell(iSt, 3, rmsPullT);
+      pPullsTable->SetCell(iSt, 1, fvph_pull_x[iSt]->GetStdDev());
+      pPullsTable->SetCell(iSt, 2, fvph_pull_y[iSt]->GetStdDev());
+      pPullsTable->SetCell(iSt, 3, fvph_pull_t[iSt]->GetStdDev());
+    }
+
+    for (int idig = 0; idig <= fkMaxDigisInClusterForPulls; idig++) {
+      cbm::qa::util::FitKaniadakisGaussian(fvph_pull_uv_Ndig[idig]);
+      res = CheckRangePull(fvph_pull_uv_Ndig[idig]) && res;
     }
 
     LOG(info) << '\n' << pPullsTable->ToString(3);
-  }
+  }  // McUsed
 
   return res;
 }
@@ -275,6 +250,8 @@ void CbmCaInputQaSts::DeInit()
   fvph_pull_v.clear();
   fvph_pull_t.clear();
 
+  fvph_pull_uv_Ndig.clear();
+
   fvph_res_x_vs_x.clear();
   fvph_res_y_vs_y.clear();
   fvph_res_u_vs_u.clear();
@@ -329,7 +306,7 @@ void CbmCaInputQaSts::FillHistograms()
     double xHit = pHit->GetX();
     double yHit = pHit->GetY();
     double zHit = pHit->GetZ();
-    double uHit = xHit * cos(stPhiF) + yHit * sin(stPhiB);
+    double uHit = xHit * cos(stPhiF) + yHit * sin(stPhiF);
     double vHit = xHit * cos(stPhiB) + yHit * sin(stPhiB);
 
     double duHit = pHit->GetDu();
@@ -363,7 +340,6 @@ void CbmCaInputQaSts::FillHistograms()
     fvph_hit_dv[nSt]->Fill(dvHit);
     fvph_hit_dt[nSt]->Fill(dtHit);
 
-
     // **********************
     // ** MC information QA
 
@@ -397,7 +373,8 @@ void CbmCaInputQaSts::FillHistograms()
       fvph_n_points_per_hit[iSt]->Fill(nMCpoints);
       fvph_n_points_per_hit[nSt]->Fill(nMCpoints);
 
-      if (nMCpoints != 1) { continue; }  // ??
+      // fill the following histos only for isolated hits produced by one track
+      if (pHitMatch->GetNofLinks() != 1) { continue; }
 
       // The best link to in the match (probably, the cut on nMCpoints is meaningless)
       assert(pHitMatch->GetNofLinks() > 0);  // Should be always true due to the cut above
@@ -463,10 +440,9 @@ void CbmCaInputQaSts::FillHistograms()
 
       // ------ Cuts
 
-      if (std::fabs(pMCPoint->GetPzOut()) < 0.01) { continue; }  // CUT ON MINIMUM MOMENTUM
+      if (pMCPoint->GetPzOut() < 0.01) { continue; }  // CUT ON MINIMUM MOMENTUM
       //if (pMCo < cuts::kMinP) { continue; }  // Momentum threshold
 
-
       fvph_point_ypos_vs_xpos[iSt]->Fill(xMC, yMC);
       fvph_point_xpos_vs_zpos[iSt]->Fill(zMC, xMC);
       fvph_point_ypos_vs_zpos[iSt]->Fill(zMC, yMC);
@@ -491,11 +467,11 @@ void CbmCaInputQaSts::FillHistograms()
       fvph_res_v_vs_v[iSt]->Fill(vHit, vRes);
       fvph_res_t_vs_t[iSt]->Fill(tHit, tRes);
 
-      fvph_pull_x_vs_x[iSt]->Fill(xHit, xRes / dxHit);
-      fvph_pull_y_vs_y[iSt]->Fill(yHit, yRes / dyHit);
-      fvph_pull_u_vs_u[iSt]->Fill(uHit, uRes / duHit);
-      fvph_pull_v_vs_v[iSt]->Fill(vHit, vRes / dvHit);
-      fvph_pull_t_vs_t[iSt]->Fill(tHit, tRes / dtHit);
+      fvph_pull_x_vs_x[iSt]->Fill(xMCs, xRes / dxHit);
+      fvph_pull_y_vs_y[iSt]->Fill(yMCs, yRes / dyHit);
+      fvph_pull_u_vs_u[iSt]->Fill(uMCs, uRes / duHit);
+      fvph_pull_v_vs_v[iSt]->Fill(vMCs, vRes / dvHit);
+      fvph_pull_t_vs_t[iSt]->Fill(tMCs, tRes / dtHit);
 
       fvph_point_ypos_vs_xpos[nSt]->Fill(xMC, yMC);
       fvph_point_xpos_vs_zpos[nSt]->Fill(zMC, xMC);
@@ -515,17 +491,33 @@ void CbmCaInputQaSts::FillHistograms()
       fvph_pull_v[nSt]->Fill(vRes / dvHit);
       fvph_pull_t[nSt]->Fill(tRes / dtHit);
 
-      fvph_res_x_vs_x[nSt]->Fill(xHit, xRes);
-      fvph_res_y_vs_y[nSt]->Fill(yHit, yRes);
-      fvph_res_u_vs_u[nSt]->Fill(uHit, uRes);
-      fvph_res_v_vs_v[nSt]->Fill(vHit, vRes);
-      fvph_res_t_vs_t[nSt]->Fill(tHit, tRes);
-
-      fvph_pull_x_vs_x[nSt]->Fill(xHit, xRes / dxHit);
-      fvph_pull_y_vs_y[nSt]->Fill(yHit, yRes / dyHit);
-      fvph_pull_u_vs_u[nSt]->Fill(uHit, uRes / duHit);
-      fvph_pull_v_vs_v[nSt]->Fill(vHit, vRes / dvHit);
-      fvph_pull_t_vs_t[nSt]->Fill(tHit, tRes / dtHit);
+      fvph_res_x_vs_x[nSt]->Fill(xMCs, xRes);
+      fvph_res_y_vs_y[nSt]->Fill(yMCs, yRes);
+      fvph_res_u_vs_u[nSt]->Fill(uMCs, uRes);
+      fvph_res_v_vs_v[nSt]->Fill(vMCs, vRes);
+      fvph_res_t_vs_t[nSt]->Fill(tMCs, tRes);
+
+      fvph_pull_x_vs_x[nSt]->Fill(xMCs, xRes / dxHit);
+      fvph_pull_y_vs_y[nSt]->Fill(yMCs, yRes / dyHit);
+      fvph_pull_u_vs_u[nSt]->Fill(uMCs, uRes / duHit);
+      fvph_pull_v_vs_v[nSt]->Fill(vMCs, vRes / dvHit);
+      fvph_pull_t_vs_t[nSt]->Fill(tMCs, tRes / dtHit);
+
+      {
+        const auto* pCluster = dynamic_cast<const CbmStsCluster*>(fpClusters->At(pHit->GetFrontClusterId()));
+        assert(pCluster);
+        int ndigis = pCluster->GetNofDigis();
+        if (ndigis > fkMaxDigisInClusterForPulls) { ndigis = 0; }
+        fvph_pull_uv_Ndig[ndigis]->Fill(uRes / duHit);
+      }
+
+      {
+        const auto* pCluster = dynamic_cast<const CbmStsCluster*>(fpClusters->At(pHit->GetBackClusterId()));
+        assert(pCluster);
+        int ndigis = pCluster->GetNofDigis();
+        if (ndigis > fkMaxDigisInClusterForPulls) { ndigis = 0; }
+        fvph_pull_uv_Ndig[ndigis]->Fill(vRes / dvHit);
+      }
     }
   }  // loop over hits: end
 
@@ -638,7 +630,7 @@ InitStatus CbmCaInputQaSts::InitCanvases()
   // ** Hit occupancies
 
   // ** Occupancy in XY plane
-  auto* pc_hit_ypos_vs_xpos = MakeCanvas<CbmQaCanvas>("hit_ypos_vs_xpos", "", 1600, 800);
+  auto* pc_hit_ypos_vs_xpos = MakeQaObject<CbmQaCanvas>("hit_ypos_vs_xpos", "", 1600, 800);
   pc_hit_ypos_vs_xpos->Divide2D(nSt);
   for (int iSt = 0; iSt < nSt; ++iSt) {
     pc_hit_ypos_vs_xpos->cd(iSt + 1);
@@ -675,7 +667,7 @@ InitStatus CbmCaInputQaSts::InitCanvases()
   }
 
   // ----- Occupancy projections
-  auto* pc_hit_xpos = MakeCanvas<CbmQaCanvas>("hit_xpos", "", 1400, 800);
+  auto* pc_hit_xpos = MakeQaObject<CbmQaCanvas>("hit_xpos", "", 1400, 800);
   pc_hit_xpos->Divide2D(nSt);
   for (int iSt = 0; iSt < nSt; ++iSt) {
     pc_hit_xpos->cd(iSt + 1);
@@ -683,7 +675,7 @@ InitStatus CbmCaInputQaSts::InitCanvases()
     phProj->DrawCopy();
   }
 
-  auto* pc_hit_ypos = MakeCanvas<CbmQaCanvas>("hit_ypos", "", 1400, 800);
+  auto* pc_hit_ypos = MakeQaObject<CbmQaCanvas>("hit_ypos", "", 1400, 800);
   pc_hit_ypos->Divide2D(nSt);
   for (int iSt = 0; iSt < nSt; ++iSt) {
     pc_hit_ypos->cd(iSt + 1);
@@ -693,7 +685,7 @@ InitStatus CbmCaInputQaSts::InitCanvases()
 
 
   // ** Occupancy in XZ plane
-  MakeCanvas<CbmQaCanvas>("hit_xpos_vs_zpos", "", 600, 400);
+  MakeQaObject<CbmQaCanvas>("hit_xpos_vs_zpos", "", 600, 400);
   fvph_hit_xpos_vs_zpos[nSt]->DrawCopy("colz", "");
   for (int iSt = 0; iSt < nSt; ++iSt) {
     // Station positions in detector IFS
@@ -711,7 +703,7 @@ InitStatus CbmCaInputQaSts::InitCanvases()
   }
 
   // ** Occupancy in YZ plane
-  MakeCanvas<CbmQaCanvas>("hit_ypos_vs_zpos", "", 600, 400);
+  MakeQaObject<CbmQaCanvas>("hit_ypos_vs_zpos", "", 600, 400);
   fvph_hit_ypos_vs_zpos[nSt]->DrawCopy("colz", "");
   for (int iSt = 0; iSt < nSt; ++iSt) {
     // Station positions in detector IFS
@@ -737,7 +729,7 @@ InitStatus CbmCaInputQaSts::InitCanvases()
     // ** Point occupancies
 
     // ** Occupancy in XY plane
-    auto* pc_point_ypos_vs_xpos = MakeCanvas<CbmQaCanvas>("point_ypos_vs_xpos", "", 1600, 800);
+    auto* pc_point_ypos_vs_xpos = MakeQaObject<CbmQaCanvas>("point_ypos_vs_xpos", "", 1600, 800);
     pc_point_ypos_vs_xpos->Divide2D(nSt);
     for (int iSt = 0; iSt < nSt; ++iSt) {
       pc_point_ypos_vs_xpos->cd(iSt + 1);
@@ -773,7 +765,7 @@ InitStatus CbmCaInputQaSts::InitCanvases()
     }
 
     // ** Occupancy in XZ plane
-    MakeCanvas<CbmQaCanvas>("point_xpos_vs_zpos", "", 600, 400);
+    MakeQaObject<CbmQaCanvas>("point_xpos_vs_zpos", "", 600, 400);
     fvph_point_xpos_vs_zpos[nSt]->SetStats(false);
     fvph_point_xpos_vs_zpos[nSt]->DrawCopy("colz", "");
     for (int iSt = 0; iSt < nSt; ++iSt) {
@@ -792,7 +784,7 @@ InitStatus CbmCaInputQaSts::InitCanvases()
     }
 
     // ** Occupancy in YZ plane
-    MakeCanvas<CbmQaCanvas>("point_ypos_vs_zpos", "", 600, 400);
+    MakeQaObject<CbmQaCanvas>("point_ypos_vs_zpos", "", 600, 400);
     fvph_point_ypos_vs_zpos[nSt]->SetStats(false);
     fvph_point_ypos_vs_zpos[nSt]->DrawCopy("colz", "");
     for (int iSt = 0; iSt < nSt; ++iSt) {
@@ -814,7 +806,7 @@ InitStatus CbmCaInputQaSts::InitCanvases()
     // ** Residuals
 
     // x-coordinate
-    auto* pc_res_x = MakeCanvas<CbmQaCanvas>("res_x", "Residuals for x coordinate", 1600, 800);
+    auto* pc_res_x = MakeQaObject<CbmQaCanvas>("res_x", "Residuals for x coordinate", 1600, 800);
     pc_res_x->Divide2D(nSt);
     for (int iSt = 0; iSt < nSt; ++iSt) {
       pc_res_x->cd(iSt + 1);
@@ -822,7 +814,7 @@ InitStatus CbmCaInputQaSts::InitCanvases()
     }
 
     // y-coordinate
-    auto* pc_res_y = MakeCanvas<CbmQaCanvas>("res_y", "Residuals for y coordinate", 1600, 800);
+    auto* pc_res_y = MakeQaObject<CbmQaCanvas>("res_y", "Residuals for y coordinate", 1600, 800);
     pc_res_y->Divide2D(nSt);
     for (int iSt = 0; iSt < nSt; ++iSt) {
       pc_res_y->cd(iSt + 1);
@@ -830,7 +822,7 @@ InitStatus CbmCaInputQaSts::InitCanvases()
     }
 
     // x-coordinate
-    auto* pc_res_u = MakeCanvas<CbmQaCanvas>("res_u", "Residuals for u coordinate", 1600, 800);
+    auto* pc_res_u = MakeQaObject<CbmQaCanvas>("res_u", "Residuals for u coordinate", 1600, 800);
     pc_res_u->Divide2D(nSt);
     for (int iSt = 0; iSt < nSt; ++iSt) {
       pc_res_u->cd(iSt + 1);
@@ -838,7 +830,7 @@ InitStatus CbmCaInputQaSts::InitCanvases()
     }
 
     // x-coordinate
-    auto* pc_res_v = MakeCanvas<CbmQaCanvas>("res_v", "Residuals for v coordinate", 1600, 800);
+    auto* pc_res_v = MakeQaObject<CbmQaCanvas>("res_v", "Residuals for v coordinate", 1600, 800);
     pc_res_v->Divide2D(nSt);
     for (int iSt = 0; iSt < nSt; ++iSt) {
       pc_res_v->cd(iSt + 1);
@@ -846,7 +838,7 @@ InitStatus CbmCaInputQaSts::InitCanvases()
     }
 
     // x-coordinate
-    auto* pc_res_t = MakeCanvas<CbmQaCanvas>("res_t", "Residuals for time", 1600, 800);
+    auto* pc_res_t = MakeQaObject<CbmQaCanvas>("res_t", "Residuals for time", 1600, 800);
     pc_res_t->Divide2D(nSt);
     for (int iSt = 0; iSt < nSt; ++iSt) {
       pc_res_t->cd(iSt + 1);
@@ -857,44 +849,49 @@ InitStatus CbmCaInputQaSts::InitCanvases()
     // **********
     // ** Pulls
 
-    // x-coordinate
-    auto* pc_pull_x = MakeCanvas<CbmQaCanvas>("pull_x", "Pulls for x coordinate", 1600, 800);
-    pc_pull_x->Divide2D(nSt);
-    for (int iSt = 0; iSt < nSt; ++iSt) {
-      pc_pull_x->cd(iSt + 1);
-      fvph_pull_x[iSt]->DrawCopy("", "");
+    {  // x-coordinate
+      auto* pc_pull_x = MakeQaObject<CbmQaCanvas>("pull_x", "Pulls for x coordinate", 1600, 800);
+      pc_pull_x->Divide2D(nSt);
+      for (int iSt = 0; iSt < nSt; ++iSt) {
+        pc_pull_x->cd(iSt + 1);
+        fvph_pull_x[iSt]->DrawCopy("", "");
+      }
     }
 
-    // y-coordinate
-    auto* pc_pull_y = MakeCanvas<CbmQaCanvas>("pull_y", "Pulls for y coordinate", 1600, 800);
-    pc_pull_y->Divide2D(nSt);
-    for (int iSt = 0; iSt < nSt; ++iSt) {
-      pc_pull_y->cd(iSt + 1);
-      fvph_pull_y[iSt]->DrawCopy("", "");
+    {  // y-coordinate
+      auto* pc_pull_y = MakeQaObject<CbmQaCanvas>("pull_y", "Pulls for y coordinate", 1600, 800);
+      pc_pull_y->Divide2D(nSt);
+      for (int iSt = 0; iSt < nSt; ++iSt) {
+        pc_pull_y->cd(iSt + 1);
+        fvph_pull_y[iSt]->DrawCopy("", "");
+      }
     }
 
-    // x-coordinate
-    auto* pc_pull_u = MakeCanvas<CbmQaCanvas>("pull_u", "Pulls for u coordinate", 1600, 800);
-    pc_pull_u->Divide2D(nSt);
-    for (int iSt = 0; iSt < nSt; ++iSt) {
-      pc_pull_u->cd(iSt + 1);
-      fvph_pull_u[iSt]->DrawCopy("", "");
+    {  // u-coordinate
+      auto* pc_pull_u = MakeQaObject<CbmQaCanvas>("pull_u", "Pulls for u coordinate", 1600, 800);
+      pc_pull_u->Divide2D(nSt);
+      for (int iSt = 0; iSt < nSt; ++iSt) {
+        pc_pull_u->cd(iSt + 1);
+        fvph_pull_u[iSt]->DrawCopy("", "");
+      }
     }
 
-    // x-coordinate
-    auto* pc_pull_v = MakeCanvas<CbmQaCanvas>("pull_v", "Pulls for v coordinate", 1600, 800);
-    pc_pull_v->Divide2D(nSt);
-    for (int iSt = 0; iSt < nSt; ++iSt) {
-      pc_pull_v->cd(iSt + 1);
-      fvph_pull_v[iSt]->DrawCopy("", "");
+    {  // v-coordinate
+      auto* pc_pull_v = MakeQaObject<CbmQaCanvas>("pull_v", "Pulls for v coordinate", 1600, 800);
+      pc_pull_v->Divide2D(nSt);
+      for (int iSt = 0; iSt < nSt; ++iSt) {
+        pc_pull_v->cd(iSt + 1);
+        fvph_pull_v[iSt]->DrawCopy("", "");
+      }
     }
 
-    // x-coordinate
-    auto* pc_pull_t = MakeCanvas<CbmQaCanvas>("pull_t", "Pulls for time", 1600, 800);
-    pc_pull_t->Divide2D(nSt);
-    for (int iSt = 0; iSt < nSt; ++iSt) {
-      pc_pull_t->cd(iSt + 1);
-      fvph_pull_t[iSt]->DrawCopy("", "");
+    {  // x-coordinate
+      auto* pc_pull_t = MakeQaObject<CbmQaCanvas>("pull_t", "Pulls for time", 1600, 800);
+      pc_pull_t->Divide2D(nSt);
+      for (int iSt = 0; iSt < nSt; ++iSt) {
+        pc_pull_t->cd(iSt + 1);
+        fvph_pull_t[iSt]->DrawCopy("", "");
+      }
     }
 
 
@@ -902,7 +899,7 @@ InitStatus CbmCaInputQaSts::InitCanvases()
     // ** Hit reconstruction efficiencies
 
     auto* pc_reco_eff_vs_r =
-      MakeCanvas<CbmQaCanvas>("reco_eff_vs_r", "Hit efficiencies wrt distance from center", 1600, 800);
+      MakeQaObject<CbmQaCanvas>("reco_eff_vs_r", "Hit efficiencies wrt distance from center", 1600, 800);
     pc_reco_eff_vs_r->Divide2D(nSt);
     for (int iSt = 0; iSt < nSt; ++iSt) {
       pc_reco_eff_vs_r->cd(iSt + 1);
@@ -910,7 +907,7 @@ InitStatus CbmCaInputQaSts::InitCanvases()
       fvpe_reco_eff_vs_r[iSt]->DrawCopy("AP", "");
     }
 
-    auto* pc_reco_eff_vs_xy = MakeCanvas<CbmQaCanvas>("reco_eff_vs_xy", "Hit efficiencies wrt x and y", 1600, 800);
+    auto* pc_reco_eff_vs_xy = MakeQaObject<CbmQaCanvas>("reco_eff_vs_xy", "Hit efficiencies wrt x and y", 1600, 800);
     pc_reco_eff_vs_xy->Divide2D(nSt);
     for (int iSt = 0; iSt < nSt; ++iSt) {
       pc_reco_eff_vs_xy->cd(iSt + 1);
@@ -927,6 +924,9 @@ InitStatus CbmCaInputQaSts::InitHistograms()
 {
   int nSt = fpDetInterface->GetNtrackingStations();
 
+  //SetStoringMode(EStoringMode::kSAMEDIR);
+  SetStoringMode(EStoringMode::kSAMEDIR);
+
   // ----- Histograms initialization
   fvph_hit_ypos_vs_xpos.resize(nSt + 1, nullptr);
   fvph_hit_ypos_vs_zpos.resize(nSt + 1, nullptr);
@@ -941,50 +941,55 @@ InitStatus CbmCaInputQaSts::InitHistograms()
   fvph_hit_du.resize(nSt + 1, nullptr);
 
   for (int iSt = 0; iSt <= nSt; ++iSt) {
+    TString sF = "";  // Subfolder
+    sF += (iSt == nSt) ? "All stations/" : Form("Station %d/", iSt);
     TString nsuff = (iSt == nSt) ? "" : Form("_st%d", iSt);               // Histogram name suffix
     TString tsuff = (iSt == nSt) ? "" : Form(" in STS station %d", iSt);  // Histogram title suffix
     TString sN    = "";
     TString sT    = "";
 
     // Hit occupancy
-    sN                         = (TString) "hit_ypos_vs_xpos" + nsuff;
-    sT                         = (TString) "Hit occupancy in xy-Plane" + tsuff + ";x_{hit} [cm];y_{hit} [cm]";
-    fvph_hit_ypos_vs_xpos[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, kRHitX[0], kRHitX[1], kNbins, kRHitY[0], kRHitY[1]);
+    sN = (TString) "hit_ypos_vs_xpos" + nsuff;
+    sT = (TString) "Hit occupancy in xy-Plane" + tsuff + ";x_{hit} [cm];y_{hit} [cm]";
+    fvph_hit_ypos_vs_xpos[iSt] =
+      MakeQaObject<TH2F>(sF + "occup/" + sN, sT, kNbins, kRHitX[0], kRHitX[1], kNbins, kRHitY[0], kRHitY[1]);
 
-    sN                         = (TString) "hit_xpos_vs_zpos" + nsuff;
-    sT                         = (TString) "Hit occupancy in xz-Plane" + tsuff + ";z_{hit} [cm];x_{hit} [cm]";
-    fvph_hit_xpos_vs_zpos[iSt] = MakeHisto<TH2F>(sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitX[0], kRHitX[1]);
+    sN = (TString) "hit_xpos_vs_zpos" + nsuff;
+    sT = (TString) "Hit occupancy in xz-Plane" + tsuff + ";z_{hit} [cm];x_{hit} [cm]";
+    fvph_hit_xpos_vs_zpos[iSt] =
+      MakeQaObject<TH2F>(sF + "occup/" + sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitX[0], kRHitX[1]);
 
-    sN                         = (TString) "hit_ypos_vs_zpos" + nsuff;
-    sT                         = (TString) "Hit occupancy in yz-plane" + tsuff + ";z_{hit} [cm];y_{hit} [cm]";
-    fvph_hit_ypos_vs_zpos[iSt] = MakeHisto<TH2F>(sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitY[0], kRHitY[1]);
+    sN = (TString) "hit_ypos_vs_zpos" + nsuff;
+    sT = (TString) "Hit occupancy in yz-plane" + tsuff + ";z_{hit} [cm];y_{hit} [cm]";
+    fvph_hit_ypos_vs_zpos[iSt] =
+      MakeQaObject<TH2F>(sF + "occup/" + sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitY[0], kRHitY[1]);
 
     // Hit errors
     sN               = (TString) "hit_dx" + nsuff;
     sT               = (TString) "Hit position error along x-axis" + tsuff + ";dx_{hit} [cm]";
-    fvph_hit_dx[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRHitDx[0], kRHitDx[1]);
+    fvph_hit_dx[iSt] = MakeQaObject<TH1F>(sF + "err/" + sN, sT, kNbins, kRHitDx[0], kRHitDx[1]);
 
     sN               = (TString) "hit_dy" + nsuff;
     sT               = (TString) "Hit position error along y-axis" + tsuff + ";dy_{hit} [cm]";
-    fvph_hit_dy[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRHitDy[0], kRHitDy[1]);
+    fvph_hit_dy[iSt] = MakeQaObject<TH1F>(sF + "err/" + sN, sT, kNbins, kRHitDy[0], kRHitDy[1]);
 
     sN               = (TString) "hit_du" + nsuff;
     sT               = (TString) "Hit position error across front strips" + tsuff + ";du_{hit} [cm]";
-    fvph_hit_du[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRHitDu[0], kRHitDu[1]);
+    fvph_hit_du[iSt] = MakeQaObject<TH1F>(sF + "err/" + sN, sT, kNbins, kRHitDu[0], kRHitDu[1]);
 
     sN               = (TString) "hit_dv" + nsuff;
     sT               = (TString) "Hit position error across back strips" + tsuff + ";dv_{hit} [cm]";
-    fvph_hit_dv[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRHitDv[0], kRHitDv[1]);
+    fvph_hit_dv[iSt] = MakeQaObject<TH1F>(sF + "err/" + sN, sT, kNbins, kRHitDv[0], kRHitDv[1]);
 
     sN               = (TString) "hit_dt" + nsuff;
     sT               = (TString) "Hit time error" + tsuff + ";dt_{hit} [ns]";
-    fvph_hit_dt[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRHitDt[0], kRHitDt[1]);
+    fvph_hit_dt[iSt] = MakeQaObject<TH1F>(sF + "err/" + sN, sT, kNbins, kRHitDt[0], kRHitDt[1]);
 
     if (iSt == nSt) { continue; }  // Following histograms are defined only for separate station
 
     sN = (TString) "hit_station_delta_z" + nsuff;
     sT = (TString) "Different between hit and station z-positions" + tsuff + ";z_{hit} - z_{st} [cm]";
-    fvph_hit_station_delta_z[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, -5., 5);
+    fvph_hit_station_delta_z[iSt] = MakeQaObject<TH1F>(sF + sN, sT, kNbins, -5., 5);
 
   }  // loop over station index: end
 
@@ -1004,6 +1009,9 @@ InitStatus CbmCaInputQaSts::InitHistograms()
     fvph_pull_x.resize(nSt + 1, nullptr);
     fvph_pull_y.resize(nSt + 1, nullptr);
     fvph_pull_u.resize(nSt + 1, nullptr);
+
+    fvph_pull_uv_Ndig.resize(fkMaxDigisInClusterForPulls + 1, nullptr);
+
     fvph_pull_v.resize(nSt + 1, nullptr);
     fvph_pull_t.resize(nSt + 1, nullptr);
     fvph_res_x_vs_x.resize(nSt + 1, nullptr);
@@ -1020,6 +1028,8 @@ InitStatus CbmCaInputQaSts::InitHistograms()
     fvpe_reco_eff_vs_r.resize(nSt + 1, nullptr);
 
     for (int iSt = 0; iSt <= nSt; ++iSt) {
+      TString sF = "";  // Subfolder
+      sF += (iSt == nSt) ? "All stations/" : Form("Station %d/", iSt);
       TString nsuff = (iSt == nSt) ? "" : Form("_st%d", iSt);               // Histogram name suffix
       TString tsuff = (iSt == nSt) ? "" : Form(" in STS station %d", iSt);  // Histogram title suffix
       TString sN    = "";
@@ -1027,121 +1037,134 @@ InitStatus CbmCaInputQaSts::InitHistograms()
 
       sN                         = (TString) "n_points_per_hit" + nsuff;
       sT                         = (TString) "Number of points per hit" + tsuff + ";N_{point}/hit";
-      fvph_n_points_per_hit[iSt] = MakeHisto<TH1F>(sN, sT, 10, -0.5, 9.5);
+      fvph_n_points_per_hit[iSt] = MakeQaObject<TH1F>(sF + sN, sT, 10, -0.5, 9.5);
 
       // Point occupancy
       sN = (TString) "point_ypos_vs_xpos" + nsuff;
       sT = (TString) "Point occupancy in XY plane" + tsuff + ";x_{MC} [cm];y_{MC} [cm]";
       fvph_point_ypos_vs_xpos[iSt] =
-        MakeHisto<TH2F>(sN, sT, kNbins, kRHitX[0], kRHitX[1], kNbins, kRHitY[0], kRHitY[1]);
+        MakeQaObject<TH2F>(sF + "occup/" + sN, sT, kNbins, kRHitX[0], kRHitX[1], kNbins, kRHitY[0], kRHitY[1]);
 
       sN = (TString) "point_xpos_vs_zpos" + nsuff;
       sT = (TString) "Point Occupancy in XZ plane" + tsuff + ";z_{MC} [cm];x_{MC} [cm]";
       fvph_point_xpos_vs_zpos[iSt] =
-        MakeHisto<TH2F>(sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitX[0], kRHitX[1]);
+        MakeQaObject<TH2F>(sF + "occup/" + sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitX[0], kRHitX[1]);
 
       sN = (TString) "point_ypos_vs_zpos" + nsuff;
       sT = (TString) "Point Occupancy in YZ Plane" + tsuff + ";z_{MC} [cm];y_{MC} [cm]";
       fvph_point_ypos_vs_zpos[iSt] =
-        MakeHisto<TH2F>(sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitY[0], kRHitY[1]);
+        MakeQaObject<TH2F>(sF + "occup/" + sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitY[0], kRHitY[1]);
 
       // Difference between MC input z and hit z coordinates
       sN = (TString) "point_hit_delta_z" + nsuff;
       sT = (TString) "Distance between STS point and hit along z axis" + tsuff + ";z_{reco} - z_{MC} [cm]";
-      fvph_point_hit_delta_z[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, -0.04, 0.04);
+      fvph_point_hit_delta_z[iSt] = MakeQaObject<TH1F>(sF + sN, sT, kNbins, -0.04, 0.04);
 
       sN              = (TString) "res_x" + nsuff;
       sT              = (TString) "Residuals for x-coordinate" + tsuff + ";x_{reco} - x_{MC} [cm]";
-      fvph_res_x[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRResX[0], kRResX[1]);
+      fvph_res_x[iSt] = MakeQaObject<TH1F>(sF + "res/" + sN, sT, kNbins, kRResX[0], kRResX[1]);
 
       sN              = (TString) "res_y" + nsuff;
       sT              = (TString) "Residuals for y-coordinate" + tsuff + ";y_{reco} - y_{MC} [cm]";
-      fvph_res_y[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRResY[0], kRResY[1]);
+      fvph_res_y[iSt] = MakeQaObject<TH1F>(sF + "res/" + sN, sT, kNbins, kRResY[0], kRResY[1]);
 
       sN              = (TString) "res_u" + nsuff;
       sT              = (TString) "Residuals for front strips coordinate" + tsuff + ";u_{reco} - u_{MC} [cm]";
-      fvph_res_u[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRResU[0], kRResU[1]);
+      fvph_res_u[iSt] = MakeQaObject<TH1F>(sF + "res/" + sN, sT, kNbins, kRResU[0], kRResU[1]);
 
       sN              = (TString) "res_v" + nsuff;
       sT              = (TString) "Residuals for back strips coordinate" + tsuff + ";v_{reco} - v_{MC} [cm]";
-      fvph_res_v[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRResV[1], kRResV[1]);
+      fvph_res_v[iSt] = MakeQaObject<TH1F>(sF + "res/" + sN, sT, kNbins, kRResV[1], kRResV[1]);
 
       sN              = (TString) "res_t" + nsuff;
       sT              = (TString) "Residuals for time" + tsuff + ";t_{reco} - t_{MC} [ns]";
-      fvph_res_t[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRResT[0], kRResT[1]);
+      fvph_res_t[iSt] = MakeQaObject<TH1F>(sF + "res/" + sN, sT, kNbins, kRResT[0], kRResT[1]);
 
       sN               = (TString) "pull_x" + nsuff;
       sT               = (TString) "Pulls for x-coordinate" + tsuff + ";(x_{reco} - x_{MC}) / #sigma_{x}^{reco}";
-      fvph_pull_x[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRPullX[0], kRPullX[1]);
+      fvph_pull_x[iSt] = MakeQaObject<TH1F>(sF + "pull/" + sN, sT, kNbins, kRPull[0], kRPull[1]);
 
       sN               = (TString) "pull_y" + nsuff;
       sT               = (TString) "Pulls for y-coordinate" + tsuff + ";(y_{reco} - y_{MC}) / #sigma_{y}^{reco}";
-      fvph_pull_y[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRPullY[0], kRPullY[1]);
+      fvph_pull_y[iSt] = MakeQaObject<TH1F>(sF + "pull/" + sN, sT, kNbins, kRPull[0], kRPull[1]);
 
       sN = (TString) "pull_u" + nsuff;
       sT = (TString) "Pulls for front strips coordinate" + tsuff + ";(u_{reco} - u_{MC}) / #sigma_{u}^{reco}";
-      fvph_pull_u[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRPullU[0], kRPullU[1]);
+      fvph_pull_u[iSt] = MakeQaObject<TH1F>(sF + "pull/" + sN, sT, kNbins, kRPull[0], kRPull[1]);
 
       sN = (TString) "pull_v" + nsuff;
       sT = (TString) "Pulls for back strips coordinate" + tsuff + ";(v_{reco} - v_{MC}) / #sigma_{v}^{reco}";
-      fvph_pull_v[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRPullV[0], kRPullV[1]);
+      fvph_pull_v[iSt] = MakeQaObject<TH1F>(sF + "pull/" + sN, sT, kNbins, kRPull[0], kRPull[1]);
 
       sN               = (TString) "pull_t" + nsuff;
       sT               = (TString) "Pulls for time" + tsuff + ";(t_{reco} - t_{MC}) / #sigma_{t}^{reco}";
-      fvph_pull_t[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRPullT[0], kRPullT[1]);
+      fvph_pull_t[iSt] = MakeQaObject<TH1F>(sF + "pull/" + sN, sT, kNbins, kRPull[0], kRPull[1]);
 
       sN                   = (TString) "res_x_vs_x" + nsuff;
       sT                   = (TString) "Residuals for x-coordinate" + tsuff + ";x_{reco} [cm];x_{reco} - x_{MC} [cm]";
-      fvph_res_x_vs_x[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRResX[0], kRResX[1]);
+      fvph_res_x_vs_x[iSt] = MakeQaObject<TH2F>(sF + "res/" + sN, sT, kNbins, 0, 0, kNbins, kRResX[0], kRResX[1]);
 
       sN                   = (TString) "res_y_vs_y" + nsuff;
       sT                   = (TString) "Residuals for y-coordinate" + tsuff + ";y_{reco} [cm];y_{reco} - y_{MC} [cm]";
-      fvph_res_y_vs_y[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRResY[0], kRResY[1]);
+      fvph_res_y_vs_y[iSt] = MakeQaObject<TH2F>(sF + "res/" + sN, sT, kNbins, 0, 0, kNbins, kRResY[0], kRResY[1]);
 
       sN = (TString) "res_u_vs_u" + nsuff;
       sT = (TString) "Residuals for front strips coordinate" + tsuff + ";u_{reco} [cm];u_{reco} - u_{MC} [cm]";
-      fvph_res_u_vs_u[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRResU[0], kRResU[1]);
+      fvph_res_u_vs_u[iSt] = MakeQaObject<TH2F>(sF + "res/" + sN, sT, kNbins, 0, 0, kNbins, kRResU[0], kRResU[1]);
 
       sN = (TString) "res_v_vs_u" + nsuff;
       sT = (TString) "Residuals for back strips coordinate" + tsuff + ";v_{reco} [cm];v_{reco} - v_{MC} [cm]";
-      fvph_res_v_vs_v[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRResV[0], kRResV[1]);
+      fvph_res_v_vs_v[iSt] = MakeQaObject<TH2F>(sF + "res/" + sN, sT, kNbins, 0, 0, kNbins, kRResV[0], kRResV[1]);
 
       sN                   = (TString) "res_t_vs_t" + nsuff;
       sT                   = (TString) "Residuals for time" + tsuff + "t_{reco} [ns];t_{reco} - t_{MC} [ns]";
-      fvph_res_t_vs_t[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRResT[0], kRResT[1]);
+      fvph_res_t_vs_t[iSt] = MakeQaObject<TH2F>(sF + "res/" + sN, sT, kNbins, 0, 0, kNbins, kRResT[0], kRResT[1]);
 
       sN = (TString) "pull_x_vs_x" + nsuff;
       sT = (TString) "Pulls for x-coordinate" + tsuff + "x_{reco} [cm];(x_{reco} - x_{MC}) / #sigma_{x}^{reco}";
-      fvph_pull_x_vs_x[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRPullX[0], kRPullX[1]);
+      fvph_pull_x_vs_x[iSt] = MakeQaObject<TH2F>(sF + "pull/" + sN, sT, kNbins, 0, 0, kNbins, kRPull[0], kRPull[1]);
 
       sN = (TString) "pull_y_vs_y" + nsuff;
       sT = (TString) "Pulls for y-coordinate" + tsuff + ";y_{reco} [cm];(y_{reco} - y_{MC}) / #sigma_{y}^{reco}";
-      fvph_pull_y_vs_y[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRPullY[0], kRPullY[1]);
+      fvph_pull_y_vs_y[iSt] = MakeQaObject<TH2F>(sF + "pull/" + sN, sT, kNbins, 0, 0, kNbins, kRPull[0], kRPull[1]);
 
       sN = (TString) "pull_u_vs_u" + nsuff;
       sT = (TString) "Pulls for front strips coordinate" + tsuff
            + ";u_{reco} [cm];(u_{reco} - u_{MC}) / #sigma_{u}^{reco}";
-      fvph_pull_u_vs_u[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRPullU[0], kRPullU[1]);
+      fvph_pull_u_vs_u[iSt] = MakeQaObject<TH2F>(sF + "pull/" + sN, sT, kNbins, 0, 0, kNbins, kRPull[0], kRPull[1]);
 
       sN = (TString) "pull_v_vs_v" + nsuff;
       sT =
         (TString) "Pulls for back strips coordinate" + tsuff + ";v_{reco} [cm];(v_{reco} - v_{MC}) / #sigma_{v}^{reco}";
-      fvph_pull_v_vs_v[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRPullV[0], kRPullV[1]);
+      fvph_pull_v_vs_v[iSt] = MakeQaObject<TH2F>(sF + "pull/" + sN, sT, kNbins, 0, 0, kNbins, kRPull[0], kRPull[1]);
 
       sN = (TString) "pull_t_vs_t" + nsuff;
       sT = (TString) "Pulls for time" + tsuff + ";t_{reco} [ns];(t_{reco} - t_{MC}) / #sigma_{t}^{reco}";
-      fvph_pull_t_vs_t[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRPullT[0], kRPullT[1]);
-
+      fvph_pull_t_vs_t[iSt] = MakeQaObject<TH2F>(sF + "pull/" + sN, sT, kNbins, 0, 0, kNbins, kRPull[0], kRPull[1]);
 
       sN                       = (TString) "reco_eff_vs_xy" + nsuff;
       sT                       = (TString) "Hit rec. efficiency" + tsuff + ";x_{MC} [cm];y_{MC} [cm];#epsilon";
-      fvpe_reco_eff_vs_xy[iSt] = MakeEfficiency<CbmQaEff>(sN, sT, 100, -50, 50, 100, -50, 50);
+      fvpe_reco_eff_vs_xy[iSt] = MakeQaObject<CbmQaEff>(sN, sT, 100, -50, 50, 100, -50, 50);
 
       sN                      = (TString) "reco_eff_vs_r" + nsuff;
       sT                      = (TString) "Hit rec. efficiency" + tsuff + ";r_{MC} [cm];#epsilon";
-      fvpe_reco_eff_vs_r[iSt] = MakeEfficiency<CbmQaEff>(sN, sT, 100, 0, 100);
+      fvpe_reco_eff_vs_r[iSt] = MakeQaObject<CbmQaEff>(sN, sT, 100, 0, 100);
+    }  // stations
+
+    for (int idig = 0; idig <= fkMaxDigisInClusterForPulls; idig++) {
+      TString sN = "All stations/pull/pull_uv_";
+      TString sT = "Pulls for front & back strips coordinate, hits with ";
+      if (idig == 0) {
+        sN += Form("%d+_digi", fkMaxDigisInClusterForPulls + 1);
+        sT += Form(">=%d digi", fkMaxDigisInClusterForPulls + 1);
+      }
+      else {
+        sN += Form("%d_digi", idig);
+        sT += Form("%d digi", idig);
+      }
+      sT += "; (u_{reco} - u_{MC}) / #sigma_{u}^{reco}";
+      fvph_pull_uv_Ndig[idig] = MakeQaObject<TH1F>(sN, sT, kNbins, kRPull[0], kRPull[1]);
     }
   }
   return kSUCCESS;
 }
-
diff --git a/reco/L1/qa/CbmCaInputQaSts.h b/reco/L1/qa/CbmCaInputQaSts.h
index 2a11a0c10a4620a597d2f9c8ee8f6e4e57aaaffd..d8c670aa3570d6065ff3b544a914cc25b06d1477 100644
--- a/reco/L1/qa/CbmCaInputQaSts.h
+++ b/reco/L1/qa/CbmCaInputQaSts.h
@@ -66,6 +66,13 @@ protected:
   /// De-initializes histograms
   void DeInit();
 
+  // Checks ranges for mean and standard deviation
+  // \return  False, if variable exits the range
+  bool CheckRangePull(TH1* h)
+  {
+    return CbmQaTask::CheckRange(h, fPullMeanThrsh, 1. - fPullWidthThrsh, 1. + fPullWidthThrsh);
+  }
+
 private:
   // ----- Data branches
   CbmStsTrackingInterface* fpDetInterface = nullptr;  ///< Instance of detector interface
@@ -99,15 +106,11 @@ private:
 
   static constexpr double kRResX[2] = {-.02, .02};  ///< Range for residual of x coordinate [cm]
   static constexpr double kRResY[2] = {-.1, .1};    ///< Range for residual of y coordinate [cm]
-  static constexpr double kRResU[2] = {-5., 5.};    ///< Range for residual of v coordinate [cm]
-  static constexpr double kRResV[2] = {-.4, .4};    ///< Range for residual of y coordinate [cm]
+  static constexpr double kRResU[2] = {-.02, .02};  ///< Range for residual of u coordinate [cm]
+  static constexpr double kRResV[2] = {-.02, .02};  ///< Range for residual of v coordinate [cm]
   static constexpr double kRResT[2] = {-15., 15.};  ///< Range for residual of time [ns]
 
-  static constexpr double kRPullX[2] = {-10., 10.};    ///< Range for pull of x coordinate
-  static constexpr double kRPullY[2] = {-10., 10.};    ///< Range for pull of y coordinate
-  static constexpr double kRPullU[2] = {-400., 400.};  ///< Range for pull of v coordinate
-  static constexpr double kRPullV[2] = {-10., 10.};    ///< Range for pull of y coordinate
-  static constexpr double kRPullT[2] = {-5., 5.};      ///< Range for pull of time
+  static constexpr double kRPull[2] = {-10., 10.};  ///< Range for pull histograms coordinate
 
   // NOTE: the last element of each vector stands for integral distribution over all stations
 
@@ -154,6 +157,10 @@ private:
   std::vector<TH1F*> fvph_pull_v;  ///< pull for v coordinate in different stations
   std::vector<TH1F*> fvph_pull_t;  ///< pull for t coordinate in different stations
 
+  const int fkMaxDigisInClusterForPulls {5};  ///< max digis in cluster for separate histogramming of puls
+
+  std::vector<TH1F*> fvph_pull_uv_Ndig;  ///< pull for u & v coordinates, depending on N digis in the cluster
+
   std::vector<TH2F*> fvph_pull_x_vs_x;  ///< pull for x coordinate in different stations
   std::vector<TH2F*> fvph_pull_y_vs_y;  ///< pull for y coordinate in different stations
   std::vector<TH2F*> fvph_pull_u_vs_u;  ///< pull for u coordinate in different stations