From a9638912eac6ea975168979b7034dafa464c9e98 Mon Sep 17 00:00:00 2001
From: "se.gorbunov" <se.gorbunov@gsi.de>
Date: Wed, 2 Aug 2023 15:10:29 +0000
Subject: [PATCH] ca: update of tracker input qa for sts

---
 core/data/CbmPixelHit.h                       |   6 +-
 .../detectors/sts/CbmStsTrackingInterface.cxx |  42 +
 core/detectors/sts/CbmStsTrackingInterface.h  |  36 +-
 core/qa/CbmQaCanvas.h                         |   2 +-
 core/qa/CbmQaIO.cxx                           |  59 +-
 core/qa/CbmQaIO.h                             |  37 +-
 core/qa/CbmQaUtil.cxx                         |  91 +-
 core/qa/CbmQaUtil.h                           |  13 +-
 macro/run/run_qa.C                            |   3 -
 reco/L1/qa/CbmCaInputQaSts.cxx                | 946 +++++++++++-------
 reco/L1/qa/CbmCaInputQaSts.h                  |  66 +-
 11 files changed, 904 insertions(+), 397 deletions(-)

diff --git a/core/data/CbmPixelHit.h b/core/data/CbmPixelHit.h
index 070a61f479..3821064029 100644
--- a/core/data/CbmPixelHit.h
+++ b/core/data/CbmPixelHit.h
@@ -38,7 +38,7 @@ public:
 	 * \param[in] dx X position error of the hit [cm].
 	 * \param[in] dy Y position error of the hit [cm].
 	 * \param[in] dz Z position error of the hit [cm].
-	 * \param[in] dxy XY correlation of the hit.
+	 * \param[in] dxy X-Y covariance of the hit [cm**2].
 	 * \param[in] refId Some reference ID.
          * \param[in] time Hit time [ns].   
          * \param[in] timeError Error of hit time [ns].         
@@ -51,7 +51,7 @@ public:
 	 * \param address Detector unique identifier.
 	 * \param pos Position of the hit as TVector3 [cm].
 	 * \param err Position errors of the hit as TVector3 [cm].
-	 * \param dxy XY correlation of the hit.
+	 * \param dxy X-Y covariance of the hit [cm**2].
 	 * \param refId Some reference ID.
          * \param[in] time Hit time [ns].   
          * \param[in] timeError Error of hit time [ns].         
@@ -110,7 +110,7 @@ public:
 private:
   double fX, fY;    ///< X, Y positions of hit [cm]
   double fDx, fDy;  ///< X, Y errors [cm]
-  double fDxy;      ///< XY correlation
+  double fDxy;      ///< X-Y covariance of the hit [cm**2]
 
   ClassDef(CbmPixelHit, 1);
 };
diff --git a/core/detectors/sts/CbmStsTrackingInterface.cxx b/core/detectors/sts/CbmStsTrackingInterface.cxx
index fdf43d62b0..035a3a2a0e 100644
--- a/core/detectors/sts/CbmStsTrackingInterface.cxx
+++ b/core/detectors/sts/CbmStsTrackingInterface.cxx
@@ -15,6 +15,8 @@
 #include "FairRunAna.h"
 #include <Logger.h>
 
+#include "TGeoPhysicalNode.h"
+
 ClassImp(CbmStsTrackingInterface)
 
   CbmStsTrackingInterface* CbmStsTrackingInterface::fpInstance = nullptr;
@@ -49,6 +51,46 @@ double CbmStsTrackingInterface::GetStripsStereoAngleBack(int stationId) const
   return station->GetSensorRotation() + station->GetSensorStereoAngle(1) * TMath::Pi() / 180.;
 }
 
+//-------------------------------------------------------------------------------------------------------------------------------------
+//
+std::tuple<double, double> CbmStsTrackingInterface::GetStripStereoAnglesSensor(int address) const
+{
+  /// Gets front & back strip angle
+
+  const CbmStsSensor* sensor =
+    dynamic_cast<const CbmStsSensor*>(CbmStsSetup::Instance()->GetElement(address, EStsElementLevel::kStsSensor));
+
+  assert(sensor);
+
+  const TGeoPhysicalNode* node = sensor->GetNode();
+  assert(node);
+  const CbmStsParSensor* param = sensor->GetParams();
+  assert(param);
+
+  double angleF = param->GetPar(8) * TMath::DegToRad();
+  double angleB = param->GetPar(9) * TMath::DegToRad();
+
+  const TGeoHMatrix* matrix = node->GetMatrix();
+
+  if (matrix) {
+    TGeoRotation rot(*matrix);
+    {
+      Double_t local[3] = {TMath::Cos(angleF), TMath::Sin(angleF), 0.};
+      Double_t global[3];
+      rot.LocalToMaster(local, global);
+      angleF = TMath::ATan2(global[1], global[0]);
+    }
+    {
+      Double_t local[3] = {TMath::Cos(angleB), TMath::Sin(angleB), 0.};
+      Double_t global[3];
+      rot.LocalToMaster(local, global);
+      angleB = TMath::ATan2(global[1], global[0]);
+    }
+  }
+
+  return std::tuple(angleF, angleB);
+}
+
 //-------------------------------------------------------------------------------------------------------------------------------------
 //
 InitStatus CbmStsTrackingInterface::Init()
diff --git a/core/detectors/sts/CbmStsTrackingInterface.h b/core/detectors/sts/CbmStsTrackingInterface.h
index 44dddd5519..26694a34a2 100644
--- a/core/detectors/sts/CbmStsTrackingInterface.h
+++ b/core/detectors/sts/CbmStsTrackingInterface.h
@@ -81,6 +81,11 @@ public:
   /// \return Absolute stereo angle for front strips [rad]
   double GetStripsStereoAngleFront(int stationId) const;
 
+  /// Gets front and back strip stereo angle
+  /// \param  address  Unique element address
+  /// \return std::tie(front,back) - Stereo angle for the front and back strips [rad]
+  std::tuple<double, double> GetStripStereoAnglesSensor(int address) const;
+
   /// Gets station thickness along the Z-axis
   /// \param  stationId  Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
   /// \return Station thickness [cm]
@@ -96,20 +101,47 @@ public:
   /// \return Local index of the tracking station
   int GetTrackingStationIndex(int address) const { return CbmStsSetup::Instance()->GetStationNumber(address); }
 
+  /// Gets min size of a station along the X-axis
+  /// \param  stationId  Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
+  /// \return Size of station along the X-axis [cm]
+  double GetXmin(int stationId) const { return GetStsStation(stationId)->GetXmin(); }
+
   /// Gets max size of a station along the X-axis
   /// \param  stationId  Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
   /// \return Size of station along the X-axis [cm]
   double GetXmax(int stationId) const { return GetStsStation(stationId)->GetXmax(); }
 
+  /// Gets min size of a station along the Y-axis
+  /// \param  stationId  Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
+  /// \return Size of station along the Y-axis [cm]
+  double GetYmin(int stationId) const { return GetStsStation(stationId)->GetYmin(); }
+
   /// Gets max size of a station along the Y-axis
   /// \param  stationId  Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
   /// \return Size of station along the Y-axis [cm]
   double GetYmax(int stationId) const { return GetStsStation(stationId)->GetYmax(); }
 
-  /// Gets z component of the station position
+  /// Gets reference Z position of the station: approximately == mean Z of all its component
   /// \param  stationId  Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
   /// \return Z position of station [cm]
-  double GetZ(int stationId) const { return GetStsStation(stationId)->GetZ(); }
+  double GetZref(int stationId) const { return GetStsStation(stationId)->GetZ(); }
+
+  /// Gets reference Z position of the station: approximately == mean Z of all its component
+  /// \param  stationId  Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
+  /// \return Z position of station [cm]
+  double GetZ(int stationId) const { return GetZref(stationId); }
+
+  /// Gets minimal z of the station
+  /// \param  stationId  Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
+  /// \return min z position of station [cm]
+  // TODO: calculate from components
+  double GetZmin(int stationId) const { return GetStsStation(stationId)->GetZ() - 2.; }
+
+  /// Gets maximal z of the station
+  /// \param  stationId  Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
+  /// \return max z position of station [cm]
+  // TODO: calculate from components
+  double GetZmax(int stationId) const { return GetStsStation(stationId)->GetZ() + 2.; }
 
   /// Check if station provides time measurements
   /// \param  stationId  Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
diff --git a/core/qa/CbmQaCanvas.h b/core/qa/CbmQaCanvas.h
index bce7112517..f22cb77fcb 100644
--- a/core/qa/CbmQaCanvas.h
+++ b/core/qa/CbmQaCanvas.h
@@ -52,7 +52,7 @@ public:
   {
     /// the static variable will be initialised at the first call;
     /// deleted at the application end (c++11)
-    static CbmQaCanvas dummy("CbmQaTempCanvas", "CbmQaTempCanvas", 1, 1);
+    static CbmQaCanvas dummy("CbmQaDummyCanvas", "CbmQaDummyCanvas", 1, 1);
     return dummy;
   }
 
diff --git a/core/qa/CbmQaIO.cxx b/core/qa/CbmQaIO.cxx
index e5ba4c24cb..22851af947 100644
--- a/core/qa/CbmQaIO.cxx
+++ b/core/qa/CbmQaIO.cxx
@@ -9,7 +9,11 @@
 
 #include "CbmQaIO.h"
 
+#include "CbmQaCanvas.h"
+#include "CbmQaUtil.h"
+
 #include "TFile.h"
+#include "TPaveStats.h"
 
 // ---------------------------------------------------------------------------------------------------------------------
 //
@@ -25,7 +29,48 @@ CbmQaIO::~CbmQaIO() {}
 
 // ---------------------------------------------------------------------------------------------------------------------
 //
-void CbmQaIO::SetHistoProperties(TH1* pHist) { pHist->SetStats(true); }
+void CbmQaIO::SetTH1Properties(TH1* pHist)
+{
+  // Set default histo properties
+
+  TPaveStats* stats = cbm::qa::util::GetHistStats(pHist);
+  assert(stats);
+  stats->SetY1NDC(0.7);
+  stats->SetOptStat(111110);
+  stats->SetOptFit(100001);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmQaIO::SetTH2Properties(TH2* pHist)
+{
+  // Set default histo properties
+
+  pHist->SetOption("colz");
+  pHist->SetStats(false);
+  /*
+  TPaveStats* stats = cbm::qa::util::GetHistStats(pHist);
+  assert(stats);
+  stats->SetOptStat(10);
+  stats->SetOptFit(0);
+  */
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmQaIO::SetTProfile2DProperties(TProfile2D* pHist)
+{
+  // Set default histo properties
+
+  pHist->SetOption("colz");
+  pHist->SetStats(false);
+  /*
+  TPaveStats* stats = cbm::qa::util::GetHistStats(pHist);
+  assert(stats);
+  stats->SetOptStat(10);
+  stats->SetOptFit(0);
+  */
+}
 
 // ---------------------------------------------------------------------------------------------------------------------
 //
@@ -70,3 +115,15 @@ void CbmQaIO::WriteToFile(TFile* pOutFile) const
     if (pObject) { pObject->Write(); }
   }
 }
+
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmQaIO::MakeQaDirectory(TString sDirectory)
+{
+  // Add parent directory
+  if (fsRootFolderName.Length() != 0) { sDirectory = fsRootFolderName + "/" + sDirectory; }
+
+  // Register the object in the list
+  fpvObjList->push_back(std::make_pair(nullptr, sDirectory));
+}
diff --git a/core/qa/CbmQaIO.h b/core/qa/CbmQaIO.h
index d466756d39..b99f475cd1 100644
--- a/core/qa/CbmQaIO.h
+++ b/core/qa/CbmQaIO.h
@@ -89,6 +89,15 @@ public:
   template<typename T, typename... Args>
   T* MakeQaObject(TString sName, TString sTitle, Args... args);
 
+  /// @brief Creates a directory in the output file.
+  ///        Only needed to place directories in a preffered order when doing Write()
+  ///        Call it prior to MakeQaObject() calls.
+  /// @param dirName    Name of the directory
+  /// @note Tag is defined after the last ';' symbol in the nameBase string
+  /// @example
+  ///    dirName = "/stations/station0"
+  void MakeQaDirectory(TString sName);
+
   /// @brief Creates a ROOT object
   /// @brief Sets config name
   /// @param name  A path to the config
@@ -111,7 +120,15 @@ protected:
 
   /// @brief Applies properties on the histogram created with the MakeQaObject function
   /// @param pHist  Pointer to the histogram
-  virtual void SetHistoProperties(TH1* pHist);
+  virtual void SetTH1Properties(TH1* pHist);
+
+  /// @brief Applies properties on the histogram created with the MakeQaObject function
+  /// @param pHist  Pointer to the histogram
+  virtual void SetTH2Properties(TH2* pHist);
+
+  /// @brief Applies properties on the profile 2D  created with the MakeQaObject function
+  /// @param pHist  Pointer to the profile
+  void SetTProfile2DProperties(TProfile2D* pHist);
 
   /// @brief Applies properties on the canvas created with the MakeQaObject funciton
   /// @param pCanv  Pointer to the canvas
@@ -173,10 +190,20 @@ T* CbmQaIO::ConstructAndRegisterQaObject(TString sName, Args... args)
   // Create a new object
   T* pObj = new T(sName.Data(), args...);
 
-  // Take object ownership from ROOT and apply user-defined properties
-  if constexpr (std::is_base_of_v<TH1, T>) {  // histograms
+  // Take object ownership from ROOT
+  if constexpr (std::is_base_of_v<TH1, T>) {  // all kind of histograms
     pObj->SetDirectory(0);
-    SetHistoProperties(pObj);
+  }
+
+  // apply user-defined properties
+  if constexpr (std::is_same_v<TH1F, T> || std::is_same_v<TH1D, T>) {  // histograms
+    SetTH1Properties(pObj);
+  }
+  else if constexpr (std::is_same_v<TH2F, T> || std::is_same_v<TH2D, T>) {  // histograms
+    SetTH2Properties(pObj);
+  }
+  else if constexpr (std::is_same_v<TProfile2D, T>) {  // profile 2D
+    SetTProfile2DProperties(pObj);
   }
   else if constexpr (std::is_base_of_v<TCanvas, T>) {  // canvases
     auto* pListOfCanvases = gROOT->GetListOfCanvases();
@@ -186,7 +213,7 @@ T* CbmQaIO::ConstructAndRegisterQaObject(TString sName, Args... args)
 
   // Define a "summary" subdirectory of canvases and tables
   if constexpr (std::is_base_of_v<CbmQaTable, T> || std::is_base_of_v<TCanvas, T>) {
-    if (EStoringMode::kSUBDIR == fStoringMode) { sDirectory = TString("summary/") + sDirectory; }
+    if (EStoringMode::kSUBDIR == fStoringMode) { sDirectory = TString("Summary/") + sDirectory; }
   }
 
   // Add parent directory
diff --git a/core/qa/CbmQaUtil.cxx b/core/qa/CbmQaUtil.cxx
index b28128eb07..4f02c68b36 100644
--- a/core/qa/CbmQaUtil.cxx
+++ b/core/qa/CbmQaUtil.cxx
@@ -16,12 +16,68 @@
 #include "TCanvas.h"
 #include "TF1.h"
 #include "TH1.h"
+#include "TPaletteAxis.h"
 #include "TPaveStats.h"
+#include "TStyle.h"
 #include "TVirtualPad.h"
 
 namespace cbm::qa::util
 {
 
+  // ---------------------------------------------------------------------------------------------------------------------
+  //
+  TPaveStats* GetHistStats(TH1* pHist)
+  {
+    //
+    // Returns histogram stats window. Creates is if it doesn't exists.
+    //
+
+    TPaveStats* stats = (TPaveStats*) pHist->FindObject("stats");
+    if (stats) { return stats; }
+
+    TVirtualPad* padsav = gPad;
+    auto styleSave      = gStyle->GetOptStat();
+
+    CbmQaCanvas::GetDummyCanvas().cd();
+    //pHist->SetStats(false);                                  // remove the old stat window
+    pHist->SetStats(true);                          // set the flag to create the stat window during Draw()
+    if (styleSave == 0) { gStyle->SetOptStat(1); }  // set some stat option
+
+    // protection of a warning for automatic axis binning
+    TAxis* x = pHist->GetXaxis();
+    TAxis* y = pHist->GetYaxis();
+    TAxis* z = pHist->GetZaxis();
+    assert(x);
+    assert(y);
+    assert(z);
+    double rx[2] {x->GetXmin(), x->GetXmax()};
+    double ry[2] {y->GetXmin(), y->GetXmax()};
+    double rz[2] {z->GetXmin(), z->GetXmax()};
+
+    x->SetLimits(0., 1.);
+    y->SetLimits(0., 1.);
+    z->SetLimits(0., 1.);
+
+    pHist->Draw("");
+    CbmQaCanvas::GetDummyCanvas().Update();
+    CbmQaCanvas::GetDummyCanvas().Clear();
+
+    x->SetLimits(rx[0], rx[1]);
+    y->SetLimits(ry[0], ry[1]);
+    z->SetLimits(rz[0], rz[1]);
+
+    if (styleSave == 0) { gStyle->SetOptStat(styleSave); }
+    if (padsav) padsav->cd();
+
+    stats = (TPaveStats*) pHist->FindObject("stats");
+
+    // at this point the statistics window should exist under all circumstances
+    assert(stats);
+
+    return stats;
+  }
+
+
   // ---------------------------------------------------------------------------------------------------------------------
   //
   std::tuple<double, double> FitKaniadakisGaussian(TH1* pHist)
@@ -70,21 +126,10 @@ namespace cbm::qa::util
 
     // 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)));
@@ -97,18 +142,32 @@ namespace cbm::qa::util
       f->FixParameter(4, f->GetParameter(4));
     }
 
-    TPaveStats* stats = (TPaveStats*) pHist->FindObject("stats");
+    TPaveStats* stats = GetHistStats(pHist);
     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.);
   }
 
+
+  // ---------------------------------------------------------------------------------------------------------------------
+  //
+  void SetLargeStats(TH1* pHist)
+  {
+    //
+    // Set large stat. window
+    //
+
+    TPaveStats* stats = GetHistStats(pHist);
+    assert(stats);
+    stats->SetX1NDC(0.6);
+    stats->SetY1NDC(0.5);
+    stats->SetOptStat(111110);
+    stats->SetOptFit(100001);
+  }
+
+
 }  // namespace cbm::qa::util
diff --git a/core/qa/CbmQaUtil.h b/core/qa/CbmQaUtil.h
index a6a02b0150..5d0e29709d 100644
--- a/core/qa/CbmQaUtil.h
+++ b/core/qa/CbmQaUtil.h
@@ -10,18 +10,27 @@
 #ifndef CbmQaUtil_h
 #define CbmQaUtil_h 1
 
-#include "TObject.h"
-
 #include <tuple>
 class TH1;
+class TPaveStats;
 
 /// namespace cbm::qa::util contains useful utilities for CBM QA tasks
 namespace cbm::qa::util
 {
 
+  /// @brief Finds/Creates stats window for a histogram
+  /// @param pHist histogram
+  /// @return stats window
+  TPaveStats* GetHistStats(TH1* pHist);
+
   /// @brief Fit a histogram with Kaniadakis Gaussian Distribution
   /// @return  mean and std.dev of the fit
   std::tuple<double, double> FitKaniadakisGaussian(TH1* pHist);
+
+  /// @brief Set large stat. window
+  ///
+  void SetLargeStats(TH1* pHist);
+
 }  // namespace cbm::qa::util
 
 #endif  // CbmQaUtil_h
diff --git a/macro/run/run_qa.C b/macro/run/run_qa.C
index 3d451b40ab..49bc38d931 100644
--- a/macro/run/run_qa.C
+++ b/macro/run/run_qa.C
@@ -24,7 +24,6 @@
 #include "CbmMuchHitFinderQa.h"
 #include "CbmMuchTransportQa.h"
 #include "CbmSetup.h"
-#include "CbmStsFindTracksQa.h"
 
 #include <FairFileSource.h>
 #include <FairMonitor.h>
@@ -237,8 +236,6 @@ void run_qa(TString dataTra = "data/sis100_muon_jpsi_test", TString dataRaw = "d
 
   // ----- STS QA  ---------------------------------
   if (CbmSetup::Instance()->IsActive(ECbmModuleId::kSts)) {
-    //run->AddTask(new CbmStsDigitizeQa()); //opens lots of windows
-    run->AddTask(new CbmStsFindTracksQa());
     auto* pCaInputQaSts = new CbmCaInputQaSts(verbose, bUseMC);
     run->AddTask(pCaInputQaSts);
   }
diff --git a/reco/L1/qa/CbmCaInputQaSts.cxx b/reco/L1/qa/CbmCaInputQaSts.cxx
index 8489caf358..b8dd85b0db 100644
--- a/reco/L1/qa/CbmCaInputQaSts.cxx
+++ b/reco/L1/qa/CbmCaInputQaSts.cxx
@@ -74,6 +74,21 @@ bool CbmCaInputQaSts::Check()
   // ** Basic checks, available both for real and simulated data **
   // **************************************************************
 
+  // Fill efficiency distributions
+
+  for (int iSt = 0; iSt < nSt; ++iSt) {
+    TProfile2D* effxy = fvpe_reco_eff_vs_xy[iSt];
+    for (int i = 1; i < effxy->GetNbinsX() - 1; i++) {
+      for (int j = 1; j < effxy->GetNbinsY() - 1; j++) {
+        int bin = effxy->GetBin(i, j);
+        if (effxy->GetBinEntries(bin) >= 1) {
+          fvph_reco_eff[iSt]->Fill(effxy->GetBinContent(bin));
+          fvph_reco_eff[nSt]->Fill(effxy->GetBinContent(bin));
+        }
+      }
+    }
+  }
+
   // ----- Checks for mismatches in the ordering of the stations
   //
   std::vector<double> vStationPos(nSt, 0.);
@@ -84,8 +99,8 @@ bool CbmCaInputQaSts::Check()
   if (!std::is_sorted(vStationPos.cbegin(), vStationPos.cend(), [](int l, int r) { return l <= r; })) {
     if (fVerbose > 0) {
       LOG(error) << fName << ": stations are ordered improperly along the beam axis:";
-      for (auto zPos : vStationPos) {
-        LOG(error) << "\t- " << zPos;
+      for (auto z : vStationPos) {
+        LOG(error) << "\t- " << z;
       }
     }
     res = false;
@@ -133,12 +148,12 @@ bool CbmCaInputQaSts::Check()
     // Fit efficiency curves
     LOG(info) << "-- Hit efficiency integrated over hit distance from station center";
 
-    auto* pEffTable = MakeQaObject<CbmQaTable>("eff_table", "Efficiency table", nSt, 2);
+    auto* pEffTable = MakeQaObject<CbmQaTable>("vs Station/eff_table", "Efficiency table", nSt, 2);
     pEffTable->SetNamesOfCols({"Station ID", "Efficiency"});
     pEffTable->SetColWidth(20);
 
     for (int iSt = 0; iSt < nSt; ++iSt) {
-      auto eff = fvpe_reco_eff_vs_r[iSt]->GetMean(2);
+      auto eff = fvph_reco_eff[iSt]->GetMean();
       pEffTable->SetCell(iSt, 0, iSt);
       pEffTable->SetCell(iSt, 1, eff);
       res = CheckRange("Hit finder efficiency in station " + std::to_string(iSt), eff, fEffThrsh, 1.000);
@@ -150,17 +165,18 @@ bool CbmCaInputQaSts::Check()
     // Check hits for potential biases from central values
 
     auto* pResidualsTable =
-      MakeQaObject<CbmQaTable>("residuals_mean", "Residual mean values in different stations", nSt, 4);
+      MakeQaObject<CbmQaTable>("vs Station/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);
 
     // 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]);
+
+      cbm::qa::util::SetLargeStats(fvph_res_x[iSt]);
+      cbm::qa::util::SetLargeStats(fvph_res_y[iSt]);
+      cbm::qa::util::SetLargeStats(fvph_res_t[iSt]);
+      cbm::qa::util::SetLargeStats(fvph_res_u[iSt]);
+      cbm::qa::util::SetLargeStats(fvph_res_v[iSt]);
 
       pResidualsTable->SetCell(iSt, 0, iSt);
       pResidualsTable->SetCell(iSt, 1, fvph_res_x[iSt]->GetStdDev());
@@ -171,24 +187,27 @@ bool CbmCaInputQaSts::Check()
     LOG(info) << '\n' << pResidualsTable->ToString(8);
 
     // ----- Checks for pulls
+    //
     // For the hit pull is defined as a ration of the difference between hit and MC-point position or time component
     // values and an error of the hit position (or time) component, calculated in the same z-positions. In ideal case,
     // 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.
 
-    auto* pPullsTable = MakeQaObject<CbmQaTable>("pulls_rms", "Pulls std. dev. values in different stations", nSt, 4);
+    auto* pPullsTable =
+      MakeQaObject<CbmQaTable>("vs Station/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 + 1; ++iSt) {
 
       // 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]);
+
+      cbm::qa::util::SetLargeStats(fvph_pull_x[iSt]);
+      cbm::qa::util::SetLargeStats(fvph_pull_y[iSt]);
+      cbm::qa::util::SetLargeStats(fvph_pull_t[iSt]);
+      cbm::qa::util::SetLargeStats(fvph_pull_u[iSt]);
+      cbm::qa::util::SetLargeStats(fvph_pull_v[iSt]);
 
       // Check the pull quality
       res = CheckRangePull(fvph_pull_x[iSt]) && res;
@@ -204,8 +223,11 @@ bool CbmCaInputQaSts::Check()
     }
 
     for (int idig = 0; idig <= fkMaxDigisInClusterForPulls; idig++) {
-      cbm::qa::util::FitKaniadakisGaussian(fvph_pull_uv_Ndig[idig]);
-      res = CheckRangePull(fvph_pull_uv_Ndig[idig]) && res;
+      cbm::qa::util::SetLargeStats(fvph_pull_u_Ndig[idig]);
+      res = CheckRangePull(fvph_pull_u_Ndig[idig]) && res;
+
+      cbm::qa::util::SetLargeStats(fvph_pull_v_Ndig[idig]);
+      res = CheckRangePull(fvph_pull_v_Ndig[idig]) && res;
     }
 
     LOG(info) << '\n' << pPullsTable->ToString(3);
@@ -219,9 +241,9 @@ bool CbmCaInputQaSts::Check()
 void CbmCaInputQaSts::DeInit()
 {
   // Vectors with pointers to histograms
-  fvph_hit_ypos_vs_xpos.clear();
-  fvph_hit_xpos_vs_zpos.clear();
-  fvph_hit_ypos_vs_zpos.clear();
+  fvph_hit_xy.clear();
+  fvph_hit_zx.clear();
+  fvph_hit_zy.clear();
 
   fvph_hit_station_delta_z.clear();
 
@@ -233,9 +255,9 @@ void CbmCaInputQaSts::DeInit()
 
   fvph_n_points_per_hit.clear();
 
-  fvph_point_ypos_vs_xpos.clear();
-  fvph_point_xpos_vs_zpos.clear();
-  fvph_point_ypos_vs_zpos.clear();
+  fvph_point_xy.clear();
+  fvph_point_zx.clear();
+  fvph_point_zy.clear();
 
   fvph_point_hit_delta_z.clear();
 
@@ -251,7 +273,8 @@ void CbmCaInputQaSts::DeInit()
   fvph_pull_v.clear();
   fvph_pull_t.clear();
 
-  fvph_pull_uv_Ndig.clear();
+  fvph_pull_u_Ndig.clear();
+  fvph_pull_v_Ndig.clear();
 
   fvph_res_x_vs_x.clear();
   fvph_res_y_vs_y.clear();
@@ -266,7 +289,7 @@ void CbmCaInputQaSts::DeInit()
   fvph_pull_t_vs_t.clear();
 
   fvpe_reco_eff_vs_xy.clear();
-  fvpe_reco_eff_vs_r.clear();
+  fvph_reco_eff.clear();
 }
 
 // ---------------------------------------------------------------------------------------------------------------------
@@ -277,15 +300,19 @@ void CbmCaInputQaSts::FillHistograms()
   int nHits     = fpHits->GetEntriesFast();
   int nMCevents = (IsMCUsed()) ? fpMCEventList->GetNofEvents() : -1;
 
-  std::vector<std::vector<int>> vNofHitsPerPoint;  // Number of hits per MC point in different MC events
+  std::vector<std::vector<std::vector<int>>>
+    vNofHitsPerMcTrack;  // Number of hits per MC track per station in different MC events
 
   if (IsMCUsed()) {
-    vNofHitsPerPoint.resize(nMCevents);
+    vNofHitsPerMcTrack.resize(nMCevents);
     for (int iE = 0; iE < nMCevents; ++iE) {
-      int iFile   = fpMCEventList->GetFileIdByIndex(iE);
-      int iEvent  = fpMCEventList->GetEventIdByIndex(iE);
-      int nPoints = fpMCPoints->Size(iFile, iEvent);
-      vNofHitsPerPoint[iE].resize(nPoints, 0);
+      int iFile     = fpMCEventList->GetFileIdByIndex(iE);
+      int iEvent    = fpMCEventList->GetEventIdByIndex(iE);
+      int nMcTracks = fpMCTracks->Size(iFile, iEvent);
+      vNofHitsPerMcTrack[iE].resize(nSt);
+      for (int iSt = 0; iSt < nSt; iSt++) {
+        vNofHitsPerMcTrack[iE][iSt].resize(nMcTracks, 0);
+      }
     }
   }
 
@@ -300,31 +327,45 @@ void CbmCaInputQaSts::FillHistograms()
     int iSt = fpDetInterface->GetTrackingStationIndex(pHit);
     LOG_IF(fatal, iSt < 0 || iSt >= nSt) << fName << ": index of station (" << iSt << ") is out of range "
                                          << "for hit with id = " << iH;
-    double stPhiF = fpDetInterface->GetStripsStereoAngleFront(iSt);  // STS front strips stereo angle
-    double stPhiB = fpDetInterface->GetStripsStereoAngleBack(iSt);   // STS back strips stereo angle
+    double stPhiU;  // STS front strips stereo angle
+    double stPhiV;  // STS back strips stereo angle
+
+    std::tie(stPhiU, stPhiV) = fpDetInterface->GetStripStereoAnglesSensor(pHit->GetAddress());
+
+    double sinU = sin(stPhiU);
+    double cosU = cos(stPhiU);
+    double sinV = sin(stPhiV);
+    double cosV = cos(stPhiV);
+
+    // hit cov matrix
+    double Cxx = pHit->GetDx() * pHit->GetDx();
+    double Cxy = pHit->GetDxy();
+    double Cyy = pHit->GetDy() * pHit->GetDy();
 
     // Hit position
     double xHit = pHit->GetX();
     double yHit = pHit->GetY();
     double zHit = pHit->GetZ();
-    double uHit = xHit * cos(stPhiF) + yHit * sin(stPhiF);
-    double vHit = xHit * cos(stPhiB) + yHit * sin(stPhiB);
+    double uHit = xHit * cosU + yHit * sinU;
+    double vHit = xHit * cosV + yHit * sinV;
 
-    double duHit = pHit->GetDu();
-    double dvHit = pHit->GetDv();
+    // Hit errors
+    //double duHit = pHit->GetDu();
+    //double dvHit = pHit->GetDv();
+    double duHit = sqrt(cosU * cosU * Cxx + 2. * sinU * cosU * Cxy + sinU * sinU * Cyy);
+    double dvHit = sqrt(cosV * cosV * Cxx + 2. * sinV * cosV * Cxy + sinV * sinV * Cyy);
     double dxHit = pHit->GetDx();
     double dyHit = pHit->GetDy();
-    //double dxyHit = pHit->GetDxy();
 
     // Hit time
     double tHit  = pHit->GetTime();
     double dtHit = pHit->GetTimeError();
 
-    fvph_hit_ypos_vs_xpos[iSt]->Fill(xHit, yHit);
-    fvph_hit_xpos_vs_zpos[iSt]->Fill(zHit, xHit);
-    fvph_hit_ypos_vs_zpos[iSt]->Fill(zHit, yHit);
+    fvph_hit_xy[iSt]->Fill(xHit, yHit);
+    fvph_hit_zx[iSt]->Fill(zHit, xHit);
+    fvph_hit_zy[iSt]->Fill(zHit, yHit);
 
-    fvph_hit_station_delta_z[iSt]->Fill(zHit - fpDetInterface->GetZ(iSt));
+    fvph_hit_station_delta_z[iSt]->Fill(zHit - fpDetInterface->GetZref(iSt));
 
     fvph_hit_dx[iSt]->Fill(dxHit);
     fvph_hit_dy[iSt]->Fill(dyHit);
@@ -332,9 +373,12 @@ void CbmCaInputQaSts::FillHistograms()
     fvph_hit_dv[iSt]->Fill(dvHit);
     fvph_hit_dt[iSt]->Fill(dtHit);
 
-    fvph_hit_ypos_vs_xpos[nSt]->Fill(xHit, yHit);
-    fvph_hit_xpos_vs_zpos[nSt]->Fill(zHit, xHit);
-    fvph_hit_ypos_vs_zpos[nSt]->Fill(zHit, yHit);
+    fvph_hit_xy[nSt]->Fill(xHit, yHit);
+    fvph_hit_zx[nSt]->Fill(zHit, xHit);
+    fvph_hit_zy[nSt]->Fill(zHit, yHit);
+
+    fvph_hit_station_delta_z[nSt]->Fill(zHit - fpDetInterface->GetZref(iSt));
+
     fvph_hit_dx[nSt]->Fill(dxHit);
     fvph_hit_dy[nSt]->Fill(dyHit);
     fvph_hit_du[nSt]->Fill(duHit);
@@ -365,16 +409,27 @@ void CbmCaInputQaSts::FillHistograms()
         LOG_IF(fatal, iE < 0 || iE >= nMCevents) << fName << ": id of MC event is out of range (hit id = " << iH
                                                  << ", link id = " << iLink << ", event id = " << iE << ')';
 
-        LOG_IF(fatal, iP >= (int) vNofHitsPerPoint[iE].size())
-          << fName << ": index of MC point is out of range (hit id = " << iH << ", link id = " << iLink
-          << ", event id = " << iE << ", point id = " << iP << ')';
-        vNofHitsPerPoint[iE][iP]++;
+        // matched point
+        const auto* pMCPoint = dynamic_cast<const CbmStsPoint*>(fpMCPoints->Get(link));
+        LOG_IF(fatal, !pMCPoint) << fName << ": MC point object does not exist for hit " << iH;
+
+        int iTr = pMCPoint->GetTrackID();
+
+        LOG_IF(fatal, iTr >= (int) vNofHitsPerMcTrack[iE][iSt].size())
+          << fName << ": index of MC track is out of range (hit id = " << iH << ", link id = " << iLink
+          << ", event id = " << iE << ", track id = " << iTr << ')';
+
+        vNofHitsPerMcTrack[iE][iSt][iTr]++;
       }
 
       fvph_n_points_per_hit[iSt]->Fill(nMCpoints);
       fvph_n_points_per_hit[nSt]->Fill(nMCpoints);
 
-      // fill the following histos only for isolated hits produced by one track
+      //
+      // Fill the following histograms exclusively for isolated hits with a one-to-one correspondence to the mc track
+      // The mc track must satisfy the cuts
+      //
+
       if (pHitMatch->GetNofLinks() != 1) { continue; }
 
       // The best link to in the match (probably, the cut on nMCpoints is meaningless)
@@ -397,6 +452,13 @@ void CbmCaInputQaSts::FillHistograms()
       double t0MC = fpMCEventList->GetEventTime(bestPointLink);
       LOG_IF(fatal, t0MC < 0) << fName << ": MC time zero is lower then 0 ns: " << t0MC;
 
+      // cut on the mc track quality
+      if (!IsTrackSelected(pMCTrack, pMCPoint)) { continue; }
+
+      {  // skip the case when the mc point participates to several hits
+        int iE = fpMCEventList->GetEventIndex(bestTrackLink);
+        if (vNofHitsPerMcTrack[iE][iSt][pMCPoint->GetTrackID()] != 1) { continue; }
+      }
 
       // ----- MC point properties
       //
@@ -407,8 +469,6 @@ void CbmCaInputQaSts::FillHistograms()
       // Entrance position and time
       double xMC = pMCPoint->GetXIn();
       double yMC = pMCPoint->GetYIn();
-      //double uMC = xMC * cos(stPhiF) + yMC * sin(stPhiF);
-      //double vMC = xMC * cos(stPhiB) + yMC * sin(stPhiB);
       double zMC = pMCPoint->GetZIn();
       double tMC = pMCPoint->GetTime() + t0MC;
 
@@ -428,8 +488,8 @@ void CbmCaInputQaSts::FillHistograms()
       double shiftZ = zHit - zMC;  // Difference between hit and point z positions
       double xMCs   = xMC + shiftZ * pxMC / pzMC;
       double yMCs   = yMC + shiftZ * pyMC / pzMC;
-      double uMCs   = xMCs * cos(stPhiF) + yMCs * sin(stPhiF);
-      double vMCs   = xMCs * cos(stPhiB) + yMCs * sin(stPhiB);
+      double uMCs   = xMCs * cosU + yMCs * sinU;
+      double vMCs   = xMCs * cosV + yMCs * sinV;
       double tMCs   = tMC + shiftZ / (pzMC * phys::kSpeedOfLight) * sqrt(mass * mass + pMC * pMC);
 
       // Residuals
@@ -438,17 +498,9 @@ void CbmCaInputQaSts::FillHistograms()
       double uRes = uHit - uMCs;
       double vRes = vHit - vMCs;
       double tRes = tHit - tMCs;
+      double zRes = zHit - 0.5 * (pMCPoint->GetZIn() + pMCPoint->GetZOut());
 
-      // ------ Cuts
-
-      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);
-
-      fvph_point_hit_delta_z[iSt]->Fill(shiftZ);
+      fvph_point_hit_delta_z[iSt]->Fill(zRes);
 
       fvph_res_x[iSt]->Fill(xRes);
       fvph_res_y[iSt]->Fill(yRes);
@@ -474,11 +526,9 @@ void CbmCaInputQaSts::FillHistograms()
       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);
-      fvph_point_ypos_vs_zpos[nSt]->Fill(zMC, yMC);
+      // fill summary histos
 
-      fvph_point_hit_delta_z[nSt]->Fill(shiftZ);
+      fvph_point_hit_delta_z[nSt]->Fill(zRes);
 
       fvph_res_x[nSt]->Fill(xRes);
       fvph_res_y[nSt]->Fill(yRes);
@@ -509,7 +559,7 @@ void CbmCaInputQaSts::FillHistograms()
         assert(pCluster);
         int ndigis = pCluster->GetNofDigis();
         if (ndigis > fkMaxDigisInClusterForPulls) { ndigis = 0; }
-        fvph_pull_uv_Ndig[ndigis]->Fill(uRes / duHit);
+        fvph_pull_u_Ndig[ndigis]->Fill(uRes / duHit);
       }
 
       {
@@ -517,7 +567,7 @@ void CbmCaInputQaSts::FillHistograms()
         assert(pCluster);
         int ndigis = pCluster->GetNofDigis();
         if (ndigis > fkMaxDigisInClusterForPulls) { ndigis = 0; }
-        fvph_pull_uv_Ndig[ndigis]->Fill(vRes / dvHit);
+        fvph_pull_v_Ndig[ndigis]->Fill(vRes / dvHit);
       }
     }
   }  // loop over hits: end
@@ -528,29 +578,63 @@ void CbmCaInputQaSts::FillHistograms()
       int iFile   = fpMCEventList->GetFileIdByIndex(iE);
       int iEvent  = fpMCEventList->GetEventIdByIndex(iE);
       int nPoints = fpMCPoints->Size(iFile, iEvent);
+      int nTracks = fpMCTracks->Size(iFile, iEvent);
+
+      // If efficiency for the track at the station is already evaluated
+      std::vector<std::vector<bool>> vIsTrackProcessed(nSt);
+      for (int iSt = 0; iSt < nSt; iSt++) {
+        vIsTrackProcessed[iSt].resize(nTracks, 0);
+      }
 
       for (int iP = 0; iP < nPoints; ++iP) {
         const auto* pMCPoint = dynamic_cast<const CbmStsPoint*>(fpMCPoints->Get(iFile, iEvent, iP));
         LOG_IF(fatal, !pMCPoint) << fName << ": MC point does not exist for iFile = " << iFile
                                  << ", iEvent = " << iEvent << ", iP = " << iP;
+
         int address = pMCPoint->GetDetectorID();
         int iSt     = fpDetInterface->GetTrackingStationIndex(address);
         LOG_IF(fatal, iSt < 0 || iSt >= nSt)
           << fName << ": MC point for FEI = " << iFile << ", " << iEvent << ", " << iP << " and address " << address
           << " has wrong station index: iSt = " << iSt;
 
+
         double xMC = pMCPoint->GetXIn();
         double yMC = pMCPoint->GetYIn();
-        double rMC = sqrt(xMC * xMC + yMC * yMC);
+        double zMC = pMCPoint->GetZIn();
 
-        // Conditions under which point is accounted as reconstructed: point
-        bool ifPointHasHits = (vNofHitsPerPoint[iE][iP] > 0);
+        {
+          fvph_point_xy[iSt]->Fill(xMC, yMC);
+          fvph_point_zx[iSt]->Fill(zMC, xMC);
+          fvph_point_zy[iSt]->Fill(zMC, yMC);
+
+          fvph_point_xy[nSt]->Fill(xMC, yMC);
+          fvph_point_zx[nSt]->Fill(zMC, xMC);
+          fvph_point_zy[nSt]->Fill(zMC, yMC);
+        }
+
+        int iTr = pMCPoint->GetTrackID();
+
+        LOG_IF(fatal, iTr >= nTracks) << fName << ": index of MC track is out of range (point id = " << iP
+                                      << ", event id = " << iE << ", track id = " << iTr << ')';
+
+        const auto* pMCTrack = dynamic_cast<const CbmMCTrack*>(fpMCTracks->Get(iFile, iEvent, iTr));
+
+        LOG_IF(fatal, pMCTrack == nullptr) << fName << ": null MC track pointer for file id = " << iFile
+                                           << ", event id = " << iEvent << ", track id = " << iTr;
+
+        // check efficiency only once per mc track per station
+        if (vIsTrackProcessed[iSt][iTr]) { continue; }
+
+        vIsTrackProcessed[iSt][iTr] = true;
 
-        fvpe_reco_eff_vs_xy[iSt]->Fill(xMC, yMC, ifPointHasHits);
-        fvpe_reco_eff_vs_xy[nSt]->Fill(xMC, yMC, ifPointHasHits);
+        // cut on the mc track quality
+        if (!IsTrackSelected(pMCTrack, pMCPoint)) { continue; }
 
-        fvpe_reco_eff_vs_r[iSt]->Fill(rMC, ifPointHasHits);
-        fvpe_reco_eff_vs_r[nSt]->Fill(rMC, ifPointHasHits);
+        // Conditions under which point is accounted as reconstructed: point
+        bool ifTrackHasHits = (vNofHitsPerMcTrack[iE][iSt][iTr] > 0);
+
+        fvpe_reco_eff_vs_xy[iSt]->Fill(xMC, yMC, ifTrackHasHits);
+        fvpe_reco_eff_vs_xy[nSt]->Fill(xMC, yMC, ifTrackHasHits);
 
       }  // loop over MC-points: end
     }    // loop over MC-events: end
@@ -608,11 +692,70 @@ InitStatus CbmCaInputQaSts::InitDataBranches()
     LOG_IF(fatal, !fpHitMatches) << "\033[1;31m]" << fName << ": hit match branch is not found for STS\033[0m";
   }
 
+
+  // get hit distribution ranges from the geometry
+
+  int nSt = fpDetInterface->GetNtrackingStations();
+
+  frXmin.resize(nSt + 1, 0.);
+  frXmax.resize(nSt + 1, 0.);
+
+  frYmin.resize(nSt + 1, 0.);
+  frYmax.resize(nSt + 1, 0.);
+
+  frZmin.resize(nSt + 1, 0.);
+  frZmax.resize(nSt + 1, 0.);
+
+  for (int i = nSt; i >= 0; --i) {
+    {  // read the ranges for station i
+      int j = (i == nSt ? 0 : i);
+
+      frXmin[i] = fpDetInterface->GetXmin(j);
+      frXmax[i] = fpDetInterface->GetXmax(j);
+
+      frYmin[i] = fpDetInterface->GetYmin(j);
+      frYmax[i] = fpDetInterface->GetYmax(j);
+
+      frZmin[i] = fpDetInterface->GetZmin(j);
+      frZmax[i] = fpDetInterface->GetZmax(j);
+    }
+
+    // update overall ranges
+
+    frXmin[nSt] = std::min(frXmin[nSt], frXmin[i]);
+    frXmax[nSt] = std::max(frXmax[nSt], frXmax[i]);
+
+    frYmin[nSt] = std::min(frYmin[nSt], frYmin[i]);
+    frYmax[nSt] = std::max(frYmax[nSt], frYmax[i]);
+
+    frZmin[nSt] = std::min(frZmin[nSt], frZmin[i]);
+    frZmax[nSt] = std::max(frZmax[nSt], frZmax[i]);
+  }
+
+  // add 5% margins for a better view
+
+  for (int i = 0; i <= nSt; ++i) {
+    double dx = 0.05 * fabs(frXmax[i] - frXmin[i]);
+    frXmin[i] -= dx;
+    frXmax[i] += dx;
+
+    double dy = 0.05 * fabs(frYmax[i] - frYmin[i]);
+    frYmin[i] -= dy;
+    frYmax[i] += dy;
+
+    double dz = 0.05 * fabs(frZmax[i] - frZmin[i]);
+    frZmin[i] -= dz;
+    frZmax[i] += dz;
+  }
+
+
   return kSUCCESS;
 }
 
 // ---------------------------------------------------------------------------------------------------------------------
 //
+// TODO: rename the method to Finish or Finalise or MakeSummary, since it does more than just canvases
+//
 InitStatus CbmCaInputQaSts::InitCanvases()
 {
   gStyle->SetOptFit(1);
@@ -627,135 +770,98 @@ InitStatus CbmCaInputQaSts::InitCanvases()
   constexpr auto contStyle = 2;  // Line style
   constexpr auto contFill  = 0;  // Fill style
 
-  // *********************************
-  // ** Hit occupancies
-
-  // ** Occupancy in XY plane
-  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);
-    fvph_hit_ypos_vs_xpos[iSt]->DrawCopy("colz", "");
-
-
-    // Reference size of the station saved to the detector interface
-    auto* pInnerCircle = new TEllipse(0., 0., fpDetInterface->GetRmin(iSt));
-    pInnerCircle->SetLineWidth(contWidth);
-    pInnerCircle->SetLineStyle(contStyle);
-    pInnerCircle->SetLineColor(contColor);
-    pInnerCircle->SetFillStyle(contFill);
-    pInnerCircle->Draw("SAME");
-
-    auto* pOuterCircle = new TEllipse(0., 0., fpDetInterface->GetRmax(iSt));
-    pOuterCircle->SetLineWidth(contWidth);
-    pOuterCircle->SetLineStyle(contStyle);
-    pOuterCircle->SetLineColor(contColor);
-    pOuterCircle->SetFillStyle(contFill);
-    pOuterCircle->Draw("SAME");
-
-    // Square boarders of station, which are used for the field approximation
-    double stXmin = +fpDetInterface->GetXmax(iSt);
-    double stXmax = -fpDetInterface->GetXmax(iSt);
-    double stYmin = -fpDetInterface->GetYmax(iSt);
-    double stYmax = +fpDetInterface->GetYmax(iSt);
-
-    auto* pBox = new TBox(stXmin, stYmin, stXmax, stYmax);
-    pBox->SetLineWidth(contWidth);
-    pBox->SetLineStyle(contStyle);
-    pBox->SetLineColor(contColor);
-    pBox->SetFillStyle(contFill);
-    pBox->Draw("SAME");
+  // **********************
+  // ** summary for all stations
+
+  {
+    auto* canv = MakeQaObject<CbmQaCanvas>("occ_hit", "Hit Occupancy", 1600, 800);
+    canv->Divide2D(3);
+    canv->cd(1);
+    fvph_hit_xy[nSt]->DrawCopy("colz", "");
+    canv->cd(2);
+    fvph_hit_zx[nSt]->DrawCopy("colz", "");
+    canv->cd(3);
+    fvph_hit_zy[nSt]->DrawCopy("colz", "");
   }
 
-  // ----- Occupancy projections
-  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);
-    auto* phProj = fvph_hit_ypos_vs_xpos[iSt]->ProjectionY();
-    phProj->DrawCopy();
-  }
+  {
+    auto* canv = MakeQaObject<CbmQaCanvas>("occ_point", "Point Occupancy", 1600, 800);
+    canv->Divide2D(3);
 
-  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);
-    auto* phProj = fvph_hit_ypos_vs_xpos[iSt]->ProjectionX();
-    phProj->DrawCopy();
+    canv->cd(1);
+    fvph_point_xy[nSt]->DrawCopy("colz", "");
+    canv->cd(2);
+    fvph_point_zx[nSt]->DrawCopy("colz", "");
+    canv->cd(3);
+    fvph_point_zy[nSt]->DrawCopy("colz", "");
   }
 
+  {
+    auto* canv = MakeQaObject<CbmQaCanvas>("residual", "Hit Residuals", 1600, 800);
+    canv->Divide2D(6);
+    canv->cd(1);
+    fvph_res_x[nSt]->DrawCopy("colz", "");
+    canv->cd(2);
+    fvph_res_y[nSt]->DrawCopy("colz", "");
+    canv->cd(3);
+    fvph_res_t[nSt]->DrawCopy("colz", "");
+    canv->cd(4);
+    fvph_res_u[nSt]->DrawCopy("colz", "");
+    canv->cd(5);
+    fvph_res_v[nSt]->DrawCopy("colz", "");
+  }
 
-  // ** Occupancy in XZ plane
-  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
-    double stZmin = fpDetInterface->GetZ(iSt) - 0.5 * fpDetInterface->GetThickness(iSt);
-    double stZmax = fpDetInterface->GetZ(iSt) + 0.5 * fpDetInterface->GetThickness(iSt);
-    double stHmin = -fpDetInterface->GetRmax(iSt);
-    double stHmax = +fpDetInterface->GetRmax(iSt);
-
-    auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax);
-    pBox->SetLineWidth(contWidth);
-    pBox->SetLineStyle(contStyle);
-    pBox->SetLineColor(contColor);
-    pBox->SetFillStyle(contFill);
-    pBox->Draw("SAME");
+  {
+    auto* canv = MakeQaObject<CbmQaCanvas>("pull", "Hit Pulls", 1600, 800);
+    canv->Divide2D(6);
+    canv->cd(1);
+    fvph_pull_x[nSt]->DrawCopy("colz", "");
+    canv->cd(2);
+    fvph_pull_y[nSt]->DrawCopy("colz", "");
+    canv->cd(3);
+    fvph_pull_t[nSt]->DrawCopy("colz", "");
+    canv->cd(4);
+    fvph_pull_u[nSt]->DrawCopy("colz", "");
+    canv->cd(5);
+    fvph_pull_v[nSt]->DrawCopy("colz", "");
   }
 
-  // ** Occupancy in YZ plane
-  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
-    double stZmin = fpDetInterface->GetZ(iSt) - 0.5 * fpDetInterface->GetThickness(iSt);
-    double stZmax = fpDetInterface->GetZ(iSt) + 0.5 * fpDetInterface->GetThickness(iSt);
-    double stHmin = -fpDetInterface->GetRmax(iSt);
-    double stHmax = +fpDetInterface->GetRmax(iSt);
-
-    auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax);
-    pBox->SetLineWidth(contWidth);
-    pBox->SetLineStyle(contStyle);
-    pBox->SetLineColor(contColor);
-    pBox->SetFillStyle(contFill);
-    pBox->Draw("SAME");
+  {
+    auto* canv = MakeQaObject<CbmQaCanvas>("eff", "Hit Reconstruction Efficiency", 1600, 800);
+    canv->Divide2D(2);
+    canv->cd(1);
+    fvpe_reco_eff_vs_xy[nSt]->DrawCopy("colz", "");
+    canv->cd(2);
+    fvph_reco_eff[nSt]->DrawCopy("colz", "");
   }
 
-  // ************
-  // ************
+  {
+    auto* canv = MakeQaObject<CbmQaCanvas>("other", "Other histograms", 1600, 800);
+    canv->Divide2D(3);
+    canv->cd(1);
+    fvph_n_points_per_hit[nSt]->DrawCopy("colz", "");
+    canv->cd(2);
+    fvph_hit_station_delta_z[nSt]->DrawCopy("colz", "");
+    canv->cd(3);
+    fvph_point_hit_delta_z[nSt]->DrawCopy("colz", "");
+  }
 
-  if (IsMCUsed()) {
 
-    // **********************
-    // ** Point occupancies
+  // *********************************
+  // ** Hit occupancies
 
-    // ** Occupancy in XY plane
-    auto* pc_point_ypos_vs_xpos = MakeQaObject<CbmQaCanvas>("point_ypos_vs_xpos", "", 1600, 800);
-    pc_point_ypos_vs_xpos->Divide2D(nSt);
+  {  // ** Occupancy in XY plane
+    auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/hit_xy", "", 1600, 800);
+    canv->Divide2D(nSt);
     for (int iSt = 0; iSt < nSt; ++iSt) {
-      pc_point_ypos_vs_xpos->cd(iSt + 1);
-      fvph_point_ypos_vs_xpos[iSt]->DrawCopy("colz", "");
-
-      // Reference size of the station saved to the detector interface
-      auto* pInnerCircle = new TEllipse(0., 0., fpDetInterface->GetRmin(iSt));
-      pInnerCircle->SetLineWidth(contWidth);
-      pInnerCircle->SetLineStyle(contStyle);
-      pInnerCircle->SetLineColor(contColor);
-      pInnerCircle->SetFillStyle(contFill);
-      pInnerCircle->Draw("SAME");
-
-      auto* pOuterCircle = new TEllipse(0., 0., fpDetInterface->GetRmax(iSt));
-      pOuterCircle->SetLineWidth(contWidth);
-      pOuterCircle->SetLineStyle(contStyle);
-      pOuterCircle->SetLineColor(contColor);
-      pOuterCircle->SetFillStyle(contFill);
-      pOuterCircle->Draw("SAME");
-
-      // Square boarders of station, which are used for the field approximation
-      double stXmin = +fpDetInterface->GetXmax(iSt);
-      double stXmax = -fpDetInterface->GetXmax(iSt);
-      double stYmin = -fpDetInterface->GetYmax(iSt);
-      double stYmax = +fpDetInterface->GetYmax(iSt);
+      canv->cd(iSt + 1);
+      fvph_hit_xy[iSt]->DrawCopy("colz", "");
+
+      // Square boarders of the active area of the station
+      double stXmin = fpDetInterface->GetXmin(iSt);
+      double stXmax = fpDetInterface->GetXmax(iSt);
+      double stYmin = fpDetInterface->GetYmin(iSt);
+      double stYmax = fpDetInterface->GetYmax(iSt);
 
       auto* pBox = new TBox(stXmin, stYmin, stXmax, stYmax);
       pBox->SetLineWidth(contWidth);
@@ -764,17 +870,18 @@ InitStatus CbmCaInputQaSts::InitCanvases()
       pBox->SetFillStyle(contFill);
       pBox->Draw("SAME");
     }
+  }
 
-    // ** Occupancy in XZ plane
-    MakeQaObject<CbmQaCanvas>("point_xpos_vs_zpos", "", 600, 400);
-    fvph_point_xpos_vs_zpos[nSt]->SetStats(false);
-    fvph_point_xpos_vs_zpos[nSt]->DrawCopy("colz", "");
+  {  // ** Occupancy in XZ plane
+    auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/hit_zx", "", 1600, 800);
+    canv->cd();
+    fvph_hit_zx[nSt]->DrawCopy("colz", "");
     for (int iSt = 0; iSt < nSt; ++iSt) {
       // Station positions in detector IFS
-      double stZmin = fpDetInterface->GetZ(iSt) - 0.5 * fpDetInterface->GetThickness(iSt);
-      double stZmax = fpDetInterface->GetZ(iSt) + 0.5 * fpDetInterface->GetThickness(iSt);
-      double stHmin = -fpDetInterface->GetRmax(iSt);
-      double stHmax = +fpDetInterface->GetRmax(iSt);
+      double stZmin = fpDetInterface->GetZmin(iSt);
+      double stZmax = fpDetInterface->GetZmax(iSt);
+      double stHmin = fpDetInterface->GetXmin(iSt);
+      double stHmax = fpDetInterface->GetXmax(iSt);
 
       auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax);
       pBox->SetLineWidth(contWidth);
@@ -783,17 +890,18 @@ InitStatus CbmCaInputQaSts::InitCanvases()
       pBox->SetFillStyle(contFill);
       pBox->Draw("SAME");
     }
+  }
 
-    // ** Occupancy in YZ plane
-    MakeQaObject<CbmQaCanvas>("point_ypos_vs_zpos", "", 600, 400);
-    fvph_point_ypos_vs_zpos[nSt]->SetStats(false);
-    fvph_point_ypos_vs_zpos[nSt]->DrawCopy("colz", "");
+  {  // ** Occupancy in YZ plane
+    auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/hit_zy", "", 1600, 800);
+    canv->cd();
+    fvph_hit_zy[nSt]->DrawCopy("colz", "");
     for (int iSt = 0; iSt < nSt; ++iSt) {
       // Station positions in detector IFS
-      double stZmin = fpDetInterface->GetZ(iSt) - 0.5 * fpDetInterface->GetThickness(iSt);
-      double stZmax = fpDetInterface->GetZ(iSt) + 0.5 * fpDetInterface->GetThickness(iSt);
-      double stHmin = -fpDetInterface->GetRmax(iSt);
-      double stHmax = +fpDetInterface->GetRmax(iSt);
+      double stZmin = fpDetInterface->GetZmin(iSt);
+      double stZmax = fpDetInterface->GetZmax(iSt);
+      double stHmin = fpDetInterface->GetYmin(iSt);
+      double stHmax = fpDetInterface->GetYmax(iSt);
 
       auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax);
       pBox->SetLineWidth(contWidth);
@@ -802,48 +910,124 @@ InitStatus CbmCaInputQaSts::InitCanvases()
       pBox->SetFillStyle(contFill);
       pBox->Draw("SAME");
     }
+  }
+
+  // ************
+  // ************
+
+  if (IsMCUsed()) {
+
+    // **********************
+    // ** Point occupancies
+
+    {  // ** Occupancy in XY plane
+      auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/point_xy", "", 1600, 800);
+      canv->Divide2D(nSt);
+      for (int iSt = 0; iSt < nSt; ++iSt) {
+        canv->cd(iSt + 1);
+        fvph_point_xy[iSt]->DrawCopy("colz", "");
+
+        // Square boarders of the active area of the station
+        double stXmin = fpDetInterface->GetXmin(iSt);
+        double stXmax = fpDetInterface->GetXmax(iSt);
+        double stYmin = fpDetInterface->GetYmin(iSt);
+        double stYmax = fpDetInterface->GetYmax(iSt);
+
+        auto* pBox = new TBox(stXmin, stYmin, stXmax, stYmax);
+        pBox->SetLineWidth(contWidth);
+        pBox->SetLineStyle(contStyle);
+        pBox->SetLineColor(contColor);
+        pBox->SetFillStyle(contFill);
+        pBox->Draw("SAME");
+      }
+    }
+
+    {  // ** Occupancy in XZ plane
+      auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/point_zx", "", 1600, 800);
+      canv->cd();
+      fvph_point_zx[nSt]->DrawCopy("colz", "");
+      for (int iSt = 0; iSt < nSt; ++iSt) {
+        // Station positions in detector IFS
+        double stZmin = fpDetInterface->GetZmin(iSt);
+        double stZmax = fpDetInterface->GetZmax(iSt);
+        double stHmin = fpDetInterface->GetXmin(iSt);
+        double stHmax = fpDetInterface->GetXmax(iSt);
+
+        auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax);
+        pBox->SetLineWidth(contWidth);
+        pBox->SetLineStyle(contStyle);
+        pBox->SetLineColor(contColor);
+        pBox->SetFillStyle(contFill);
+        pBox->Draw("SAME");
+      }
+    }
+
+    {  // ** Occupancy in YZ plane
+      auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/point_zy", "", 1600, 800);
+      canv->cd();
+      fvph_point_zy[nSt]->DrawCopy("colz", "");
+      for (int iSt = 0; iSt < nSt; ++iSt) {
+        // Station positions in detector IFS
+        double stZmin = fpDetInterface->GetZmin(iSt);
+        double stZmax = fpDetInterface->GetZmax(iSt);
+        double stHmin = fpDetInterface->GetYmin(iSt);
+        double stHmax = fpDetInterface->GetYmax(iSt);
+
+        auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax);
+        pBox->SetLineWidth(contWidth);
+        pBox->SetLineStyle(contStyle);
+        pBox->SetLineColor(contColor);
+        pBox->SetFillStyle(contFill);
+        pBox->Draw("SAME");
+      }
+    }
 
     // **************
     // ** Residuals
 
-    // x-coordinate
-    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);
-      fvph_res_x[iSt]->DrawCopy("", "");
+    {  // x-coordinate
+      auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/res_x", "Residuals for x coordinate", 1600, 800);
+      canv->Divide2D(nSt);
+      for (int iSt = 0; iSt < nSt; ++iSt) {
+        canv->cd(iSt + 1);
+        fvph_res_x[iSt]->DrawCopy("", "");
+      }
     }
 
-    // y-coordinate
-    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);
-      fvph_res_y[iSt]->DrawCopy("", "");
+    {  // y-coordinate
+      auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/res_y", "Residuals for y coordinate", 1600, 800);
+      canv->Divide2D(nSt);
+      for (int iSt = 0; iSt < nSt; ++iSt) {
+        canv->cd(iSt + 1);
+        fvph_res_y[iSt]->DrawCopy("", "");
+      }
     }
 
-    // x-coordinate
-    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);
-      fvph_res_u[iSt]->DrawCopy("", "");
+    {  // u-coordinate
+      auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/res_u", "Residuals for u coordinate", 1600, 800);
+      canv->Divide2D(nSt);
+      for (int iSt = 0; iSt < nSt; ++iSt) {
+        canv->cd(iSt + 1);
+        fvph_res_u[iSt]->DrawCopy("", "");
+      }
     }
 
-    // x-coordinate
-    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);
-      fvph_res_v[iSt]->DrawCopy("", "");
+    {  // v-coordinate
+      auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/res_v", "Residuals for v coordinate", 1600, 800);
+      canv->Divide2D(nSt);
+      for (int iSt = 0; iSt < nSt; ++iSt) {
+        canv->cd(iSt + 1);
+        fvph_res_v[iSt]->DrawCopy("", "");
+      }
     }
 
-    // x-coordinate
-    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);
-      fvph_res_t[iSt]->DrawCopy("", "");
+    {  // t-coordinate
+      auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/res_t", "Residuals for time", 1600, 800);
+      canv->Divide2D(nSt);
+      for (int iSt = 0; iSt < nSt; ++iSt) {
+        canv->cd(iSt + 1);
+        fvph_res_t[iSt]->DrawCopy("", "");
+      }
     }
 
 
@@ -851,68 +1035,88 @@ InitStatus CbmCaInputQaSts::InitCanvases()
     // ** Pulls
 
     {  // x-coordinate
-      auto* pc_pull_x = MakeQaObject<CbmQaCanvas>("pull_x", "Pulls for x coordinate", 1600, 800);
-      pc_pull_x->Divide2D(nSt);
+      auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/pull_x", "Pulls for x coordinate", 1600, 800);
+      canv->Divide2D(nSt);
       for (int iSt = 0; iSt < nSt; ++iSt) {
-        pc_pull_x->cd(iSt + 1);
+        canv->cd(iSt + 1);
         fvph_pull_x[iSt]->DrawCopy("", "");
       }
     }
 
     {  // y-coordinate
-      auto* pc_pull_y = MakeQaObject<CbmQaCanvas>("pull_y", "Pulls for y coordinate", 1600, 800);
-      pc_pull_y->Divide2D(nSt);
+      auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/pull_y", "Pulls for y coordinate", 1600, 800);
+      canv->Divide2D(nSt);
       for (int iSt = 0; iSt < nSt; ++iSt) {
-        pc_pull_y->cd(iSt + 1);
+        canv->cd(iSt + 1);
         fvph_pull_y[iSt]->DrawCopy("", "");
       }
     }
 
     {  // u-coordinate
-      auto* pc_pull_u = MakeQaObject<CbmQaCanvas>("pull_u", "Pulls for u coordinate", 1600, 800);
-      pc_pull_u->Divide2D(nSt);
+      auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/pull_u", "Pulls for u coordinate", 1600, 800);
+      canv->Divide2D(nSt);
       for (int iSt = 0; iSt < nSt; ++iSt) {
-        pc_pull_u->cd(iSt + 1);
+        canv->cd(iSt + 1);
         fvph_pull_u[iSt]->DrawCopy("", "");
       }
     }
 
     {  // v-coordinate
-      auto* pc_pull_v = MakeQaObject<CbmQaCanvas>("pull_v", "Pulls for v coordinate", 1600, 800);
-      pc_pull_v->Divide2D(nSt);
+      auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/pull_v", "Pulls for v coordinate", 1600, 800);
+      canv->Divide2D(nSt);
       for (int iSt = 0; iSt < nSt; ++iSt) {
-        pc_pull_v->cd(iSt + 1);
+        canv->cd(iSt + 1);
         fvph_pull_v[iSt]->DrawCopy("", "");
       }
     }
 
-    {  // x-coordinate
-      auto* pc_pull_t = MakeQaObject<CbmQaCanvas>("pull_t", "Pulls for time", 1600, 800);
-      pc_pull_t->Divide2D(nSt);
+    {  // t-coordinate
+      auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/pull_t", "Pulls for time", 1600, 800);
+      canv->Divide2D(nSt);
       for (int iSt = 0; iSt < nSt; ++iSt) {
-        pc_pull_t->cd(iSt + 1);
+        canv->cd(iSt + 1);
         fvph_pull_t[iSt]->DrawCopy("", "");
       }
     }
 
+    {  // u-coordinate vs N digis
+      auto* canv = MakeQaObject<CbmQaCanvas>("vs N digi/pull_u", "Pulls for u coordinate vs N digis", 1600, 800);
+      canv->Divide2D(fkMaxDigisInClusterForPulls + 1);
+      for (int i = 0; i <= fkMaxDigisInClusterForPulls; ++i) {
+        canv->cd((i > 0) ? i : fkMaxDigisInClusterForPulls + 1);
+        fvph_pull_u_Ndig[i]->DrawCopy("", "");
+      }
+    }
+
+    {  // v-coordinate vs N digis
+      auto* canv = MakeQaObject<CbmQaCanvas>("vs N digi/pull_v", "Pulls for v coordinate vs N digis", 1600, 800);
+      canv->Divide2D(fkMaxDigisInClusterForPulls + 1);
+      for (int i = 0; i <= fkMaxDigisInClusterForPulls; ++i) {
+        canv->cd((i > 0) ? i : fkMaxDigisInClusterForPulls + 1);
+        fvph_pull_v_Ndig[i]->DrawCopy("", "");
+      }
+    }
 
     // ************************************
     // ** Hit reconstruction efficiencies
 
-    auto* pc_reco_eff_vs_r =
-      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);
-      fvpe_reco_eff_vs_r[iSt]->SetMarkerStyle(20);
-      fvpe_reco_eff_vs_r[iSt]->DrawCopy("AP", "");
+    {
+      auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/reco_eff", "Hit efficiencies in xy bins", 1600, 800);
+      canv->Divide2D(nSt);
+      for (int iSt = 0; iSt < nSt; ++iSt) {
+        canv->cd(iSt + 1);
+        //fvph_reco_eff[iSt]->SetMarkerStyle(20);
+        fvph_reco_eff[iSt]->DrawCopy();
+      }
     }
 
-    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);
-      fvpe_reco_eff_vs_xy[iSt]->DrawCopy("colz", "");
+    {
+      auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/reco_eff_vs_xy", "Hit efficiencies wrt x and y", 1600, 800);
+      canv->Divide2D(nSt);
+      for (int iSt = 0; iSt < nSt; ++iSt) {
+        canv->cd(iSt + 1);
+        fvpe_reco_eff_vs_xy[iSt]->DrawCopy("colz", "");
+      }
     }
   }
 
@@ -926,14 +1130,20 @@ InitStatus CbmCaInputQaSts::InitHistograms()
   int nSt = fpDetInterface->GetNtrackingStations();
 
   //SetStoringMode(EStoringMode::kSAMEDIR);
-  SetStoringMode(EStoringMode::kSAMEDIR);
+  SetStoringMode(EStoringMode::kSUBDIR);
+
+  // ----- Subdirectories
+  MakeQaDirectory("Summary");
+  MakeQaDirectory("Summary/vs Station");
+  MakeQaDirectory("Summary/vs N digi");
+  MakeQaDirectory("All stations");
 
   // ----- Histograms initialization
-  fvph_hit_ypos_vs_xpos.resize(nSt + 1, nullptr);
-  fvph_hit_ypos_vs_zpos.resize(nSt + 1, nullptr);
-  fvph_hit_xpos_vs_zpos.resize(nSt + 1, nullptr);
+  fvph_hit_xy.resize(nSt + 1, nullptr);
+  fvph_hit_zy.resize(nSt + 1, nullptr);
+  fvph_hit_zx.resize(nSt + 1, nullptr);
 
-  fvph_hit_station_delta_z.resize(nSt);
+  fvph_hit_station_delta_z.resize(nSt + 1, nullptr);
 
   fvph_hit_dx.resize(nSt + 1, nullptr);
   fvph_hit_dy.resize(nSt + 1, nullptr);
@@ -949,21 +1159,32 @@ InitStatus CbmCaInputQaSts::InitHistograms()
     TString sN    = "";
     TString sT    = "";
 
+    // place directories in a prefered order
+    MakeQaDirectory(sF + "occup/");
+    if (IsMCUsed()) {
+      MakeQaDirectory(sF + "res/");
+      MakeQaDirectory(sF + "pull/");
+      MakeQaDirectory(sF + "eff/");
+    }
+    MakeQaDirectory(sF + "err/");
+
+    int nBinsZ = (iSt < nSt) ? kNbins : kNbinsZ;
+
     // Hit occupancy
-    sN = (TString) "hit_ypos_vs_xpos" + nsuff;
+    sN = (TString) "hit_xy" + 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]);
+    fvph_hit_xy[iSt] =
+      MakeQaObject<TH2F>(sF + "occup/" + sN, sT, kNbins, frXmin[iSt], frXmax[iSt], kNbins, frYmin[iSt], frYmax[iSt]);
 
-    sN = (TString) "hit_xpos_vs_zpos" + nsuff;
+    sN = (TString) "hit_zx" + 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]);
+    fvph_hit_zx[iSt] =
+      MakeQaObject<TH2F>(sF + "occup/" + sN, sT, nBinsZ, frZmin[iSt], frZmax[iSt], kNbins, frXmin[iSt], frXmax[iSt]);
 
-    sN = (TString) "hit_ypos_vs_zpos" + nsuff;
+    sN = (TString) "hit_zy" + 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]);
+    fvph_hit_zy[iSt] =
+      MakeQaObject<TH2F>(sF + "occup/" + sN, sT, nBinsZ, frZmin[iSt], frZmax[iSt], kNbins, frYmin[iSt], frYmax[iSt]);
 
     // Hit errors
     sN               = (TString) "hit_dx" + nsuff;
@@ -986,11 +1207,9 @@ InitStatus CbmCaInputQaSts::InitHistograms()
     sT               = (TString) "Hit time error" + tsuff + ";dt_{hit} [ns]";
     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] = MakeQaObject<TH1F>(sF + sN, sT, kNbins, -5., 5);
+    fvph_hit_station_delta_z[iSt] = MakeQaObject<TH1F>(sF + sN, sT, kNbins, -5., 5.);
 
   }  // loop over station index: end
 
@@ -998,9 +1217,9 @@ InitStatus CbmCaInputQaSts::InitHistograms()
   if (IsMCUsed()) {
     // Resize histogram vectors
     fvph_n_points_per_hit.resize(nSt + 1, nullptr);
-    fvph_point_ypos_vs_xpos.resize(nSt + 1, nullptr);
-    fvph_point_xpos_vs_zpos.resize(nSt + 1, nullptr);
-    fvph_point_ypos_vs_zpos.resize(nSt + 1, nullptr);
+    fvph_point_xy.resize(nSt + 1, nullptr);
+    fvph_point_zx.resize(nSt + 1, nullptr);
+    fvph_point_zy.resize(nSt + 1, nullptr);
     fvph_point_hit_delta_z.resize(nSt + 1, nullptr);
     fvph_res_x.resize(nSt + 1, nullptr);
     fvph_res_y.resize(nSt + 1, nullptr);
@@ -1011,7 +1230,8 @@ InitStatus CbmCaInputQaSts::InitHistograms()
     fvph_pull_y.resize(nSt + 1, nullptr);
     fvph_pull_u.resize(nSt + 1, nullptr);
 
-    fvph_pull_uv_Ndig.resize(fkMaxDigisInClusterForPulls + 1, nullptr);
+    fvph_pull_u_Ndig.resize(fkMaxDigisInClusterForPulls + 1, nullptr);
+    fvph_pull_v_Ndig.resize(fkMaxDigisInClusterForPulls + 1, nullptr);
 
     fvph_pull_v.resize(nSt + 1, nullptr);
     fvph_pull_t.resize(nSt + 1, nullptr);
@@ -1026,7 +1246,7 @@ InitStatus CbmCaInputQaSts::InitHistograms()
     fvph_pull_v_vs_v.resize(nSt + 1, nullptr);
     fvph_pull_t_vs_t.resize(nSt + 1, nullptr);
     fvpe_reco_eff_vs_xy.resize(nSt + 1, nullptr);
-    fvpe_reco_eff_vs_r.resize(nSt + 1, nullptr);
+    fvph_reco_eff.resize(nSt + 1, nullptr);
 
     for (int iSt = 0; iSt <= nSt; ++iSt) {
       TString sF = "";  // Subfolder
@@ -1041,120 +1261,129 @@ InitStatus CbmCaInputQaSts::InitHistograms()
       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;
+      sN = (TString) "point_xy" + nsuff;
       sT = (TString) "Point occupancy in XY plane" + tsuff + ";x_{MC} [cm];y_{MC} [cm]";
-      fvph_point_ypos_vs_xpos[iSt] =
-        MakeQaObject<TH2F>(sF + "occup/" + sN, sT, kNbins, kRHitX[0], kRHitX[1], kNbins, kRHitY[0], kRHitY[1]);
+      fvph_point_xy[iSt] =
+        MakeQaObject<TH2F>(sF + "occup/" + sN, sT, kNbins, frXmin[iSt], frXmax[iSt], kNbins, frYmin[iSt], frYmax[iSt]);
 
-      sN = (TString) "point_xpos_vs_zpos" + nsuff;
+      sN = (TString) "point_zx" + nsuff;
       sT = (TString) "Point Occupancy in XZ plane" + tsuff + ";z_{MC} [cm];x_{MC} [cm]";
-      fvph_point_xpos_vs_zpos[iSt] =
-        MakeQaObject<TH2F>(sF + "occup/" + sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitX[0], kRHitX[1]);
+      fvph_point_zx[iSt] =
+        MakeQaObject<TH2F>(sF + "occup/" + sN, sT, kNbinsZ, frZmin[iSt], frZmax[iSt], kNbins, frXmin[iSt], frXmax[iSt]);
 
-      sN = (TString) "point_ypos_vs_zpos" + nsuff;
+      sN = (TString) "point_zy" + nsuff;
       sT = (TString) "Point Occupancy in YZ Plane" + tsuff + ";z_{MC} [cm];y_{MC} [cm]";
-      fvph_point_ypos_vs_zpos[iSt] =
-        MakeQaObject<TH2F>(sF + "occup/" + sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitY[0], kRHitY[1]);
+      fvph_point_zy[iSt] =
+        MakeQaObject<TH2F>(sF + "occup/" + sN, sT, kNbinsZ, frZmin[iSt], frZmax[iSt], kNbins, frYmin[iSt], frYmax[iSt]);
 
       // 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] = MakeQaObject<TH1F>(sF + sN, sT, kNbins, -0.04, 0.04);
+      fvph_point_hit_delta_z[iSt] = MakeQaObject<TH1F>(sF + sN, sT, kNbins, -0.05, 0.05);
 
       sN              = (TString) "res_x" + nsuff;
-      sT              = (TString) "Residuals for x-coordinate" + tsuff + ";x_{reco} - x_{MC} [cm]";
+      sT              = (TString) "Residuals for X" + tsuff + ";x_{reco} - x_{MC} [cm]";
       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]";
+      sT              = (TString) "Residuals for Y" + tsuff + ";y_{reco} - y_{MC} [cm]";
       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]";
+      sT              = (TString) "Residuals for Front strip coordinate" + tsuff + ";u_{reco} - u_{MC} [cm]";
       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] = MakeQaObject<TH1F>(sF + "res/" + sN, sT, kNbins, kRResV[1], kRResV[1]);
+      sT              = (TString) "Residuals for Back strip coordinate" + tsuff + ";v_{reco} - v_{MC} [cm]";
+      fvph_res_v[iSt] = MakeQaObject<TH1F>(sF + "res/" + sN, sT, kNbins, kRResV[0], kRResV[1]);
 
       sN              = (TString) "res_t" + nsuff;
-      sT              = (TString) "Residuals for time" + tsuff + ";t_{reco} - t_{MC} [ns]";
+      sT              = (TString) "Residuals for Time" + tsuff + ";t_{reco} - t_{MC} [ns]";
       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}";
+      sT               = (TString) "Pulls for X" + tsuff + ";(x_{reco} - x_{MC}) / #sigma_{x}^{reco}";
       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}";
+      sT               = (TString) "Pulls for Y" + tsuff + ";(y_{reco} - y_{MC}) / #sigma_{y}^{reco}";
       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}";
+      sT = (TString) "Pulls for Front strip coordinate" + tsuff + ";(u_{reco} - u_{MC}) / #sigma_{u}^{reco}";
       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}";
+      sT = (TString) "Pulls for Back strip coordinate" + tsuff + ";(v_{reco} - v_{MC}) / #sigma_{v}^{reco}";
       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}";
+      sT               = (TString) "Pulls for Time" + tsuff + ";(t_{reco} - t_{MC}) / #sigma_{t}^{reco}";
       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] = MakeQaObject<TH2F>(sF + "res/" + sN, sT, kNbins, 0, 0, kNbins, kRResX[0], kRResX[1]);
+      sN = (TString) "res_x_vs_x" + nsuff;
+      sT = (TString) "Residuals for X" + tsuff + ";x_{reco} [cm];x_{reco} - x_{MC} [cm]";
+      fvph_res_x_vs_x[iSt] =
+        MakeQaObject<TH2F>(sF + "res/" + sN, sT, kNbins, frXmin[iSt], frXmax[iSt], 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] = MakeQaObject<TH2F>(sF + "res/" + sN, sT, kNbins, 0, 0, kNbins, kRResY[0], kRResY[1]);
+      sN = (TString) "res_y_vs_y" + nsuff;
+      sT = (TString) "Residuals for Y" + tsuff + ";y_{reco} [cm];y_{reco} - y_{MC} [cm]";
+      fvph_res_y_vs_y[iSt] =
+        MakeQaObject<TH2F>(sF + "res/" + sN, sT, kNbins, frYmin[iSt], frYmax[iSt], 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] = MakeQaObject<TH2F>(sF + "res/" + sN, sT, kNbins, 0, 0, kNbins, kRResU[0], kRResU[1]);
+      sT = (TString) "Residuals for Front strip coordinate" + tsuff + ";u_{reco} [cm];u_{reco} - u_{MC} [cm]";
+      fvph_res_u_vs_u[iSt] =
+        MakeQaObject<TH2F>(sF + "res/" + sN, sT, kNbins, frXmin[iSt], frXmax[iSt], 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] = MakeQaObject<TH2F>(sF + "res/" + sN, sT, kNbins, 0, 0, kNbins, kRResV[0], kRResV[1]);
+      sN = (TString) "res_v_vs_v" + nsuff;
+      sT = (TString) "Residuals for Back strip coordinate" + tsuff + ";v_{reco} [cm];v_{reco} - v_{MC} [cm]";
+      fvph_res_v_vs_v[iSt] =
+        MakeQaObject<TH2F>(sF + "res/" + sN, sT, kNbins, frXmin[iSt], frXmax[iSt], 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]";
+      sT                   = (TString) "Residuals for Time" + tsuff + "t_{reco} [ns];t_{reco} - t_{MC} [ns]";
       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] = MakeQaObject<TH2F>(sF + "pull/" + sN, sT, kNbins, 0, 0, kNbins, kRPull[0], kRPull[1]);
+      sT = (TString) "Pulls for X" + tsuff + "x_{reco} [cm];(x_{reco} - x_{MC}) / #sigma_{x}^{reco}";
+      fvph_pull_x_vs_x[iSt] =
+        MakeQaObject<TH2F>(sF + "pull/" + sN, sT, kNbins, frXmin[iSt], frXmax[iSt], 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] = MakeQaObject<TH2F>(sF + "pull/" + sN, sT, kNbins, 0, 0, kNbins, kRPull[0], kRPull[1]);
+      sT = (TString) "Pulls for Y" + tsuff + ";y_{reco} [cm];(y_{reco} - y_{MC}) / #sigma_{y}^{reco}";
+      fvph_pull_y_vs_y[iSt] =
+        MakeQaObject<TH2F>(sF + "pull/" + sN, sT, kNbins, frYmin[iSt], frYmax[iSt], 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] = MakeQaObject<TH2F>(sF + "pull/" + sN, sT, kNbins, 0, 0, kNbins, kRPull[0], kRPull[1]);
+      sT =
+        (TString) "Pulls for Front strip coordinate" + tsuff + ";u_{reco} [cm];(u_{reco} - u_{MC}) / #sigma_{u}^{reco}";
+      fvph_pull_u_vs_u[iSt] =
+        MakeQaObject<TH2F>(sF + "pull/" + sN, sT, kNbins, frXmin[iSt], frXmax[iSt], 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] = MakeQaObject<TH2F>(sF + "pull/" + sN, sT, kNbins, 0, 0, kNbins, kRPull[0], kRPull[1]);
+        (TString) "Pulls for Back strip coordinate" + tsuff + ";v_{reco} [cm];(v_{reco} - v_{MC}) / #sigma_{v}^{reco}";
+      fvph_pull_v_vs_v[iSt] =
+        MakeQaObject<TH2F>(sF + "pull/" + sN, sT, kNbins, frXmin[iSt], frXmax[iSt], 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}";
+      sT = (TString) "Pulls for Time" + tsuff + ";t_{reco} [ns];(t_{reco} - t_{MC}) / #sigma_{t}^{reco}";
       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] = MakeQaObject<TProfile2D>(sN, sT, 100, -50, 50, 100, -50, 50);
+      sT                       = (TString) "Hit rec. efficiency in XY" + tsuff + ";x_{MC} [cm];y_{MC} [cm];#epsilon";
+      fvpe_reco_eff_vs_xy[iSt] = MakeQaObject<TProfile2D>(sF + "eff/" + sN, sT, kNbins, frXmin[iSt], frXmax[iSt],
+                                                          kNbins, frYmin[iSt], frYmax[iSt]);
 
-      sN                      = (TString) "reco_eff_vs_r" + nsuff;
-      sT                      = (TString) "Hit rec. efficiency" + tsuff + ";r_{MC} [cm];#epsilon";
-      fvpe_reco_eff_vs_r[iSt] = MakeQaObject<TProfile>(sN, sT, 100, 0, 100);
+      sN                 = (TString) "reco_eff" + nsuff;
+      sT                 = (TString) "Hit rec. efficiency in XY bins" + tsuff + ";eff";
+      fvph_reco_eff[iSt] = MakeQaObject<TH1F>(sF + "eff/" + sN, sT, 130, -0.005, 1.30 - 0.005);
     }  // 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 ";
+      TString sN = "All stations/pull/Ndigi/pull_u_";
+      TString sT = "Pulls for Front strip coordinate, hits with ";
       if (idig == 0) {
         sN += Form("%d+_digi", fkMaxDigisInClusterForPulls + 1);
         sT += Form(">=%d digi", fkMaxDigisInClusterForPulls + 1);
@@ -1164,8 +1393,45 @@ InitStatus CbmCaInputQaSts::InitHistograms()
         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]);
+      fvph_pull_u_Ndig[idig] = MakeQaObject<TH1F>(sN, sT, kNbins, kRPull[0], kRPull[1]);
+    }
+
+    for (int idig = 0; idig <= fkMaxDigisInClusterForPulls; idig++) {
+      TString sN = "All stations/pull/Ndigi/pull_v_";
+      TString sT = "Pulls for Back strip 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 += "; (v_{reco} - v_{MC}) / #sigma_{v}^{reco}";
+      fvph_pull_v_Ndig[idig] = MakeQaObject<TH1F>(sN, sT, kNbins, kRPull[0], kRPull[1]);
     }
   }
   return kSUCCESS;
 }
+
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+Bool_t CbmCaInputQaSts::IsTrackSelected(const CbmMCTrack* track, const CbmStsPoint* point) const
+{
+
+  if (kTrackSelectionPrimary && track->GetMotherId() >= 0) { return false; }
+
+  double px = point->GetPxOut();
+  double py = point->GetPyOut();
+  double pz = point->GetPzOut();
+  double p  = sqrt(px * px + py * py + pz * pz);
+
+  if (pz <= 0.) { return false; }
+
+  if (p < kTrackSelectionMomentum) { return false; }
+
+  if (TMath::ATan2(sqrt(px * px + py * py), pz) * TMath::RadToDeg() > kTrackSelectionAngle) { return false; }
+
+  return true;
+}
diff --git a/reco/L1/qa/CbmCaInputQaSts.h b/reco/L1/qa/CbmCaInputQaSts.h
index 7d5fb730d6..c50b2d9ab1 100644
--- a/reco/L1/qa/CbmCaInputQaSts.h
+++ b/reco/L1/qa/CbmCaInputQaSts.h
@@ -21,13 +21,15 @@
 #include <unordered_map>
 #include <vector>
 
+class CbmMatch;
 class CbmMCEventList;
 class CbmMCDataArray;
 class CbmMCDataManager;
-class CbmTimeSlice;
-class CbmMatch;
+class CbmMCTrack;
 class CbmStsHit;
+class CbmStsPoint;
 class CbmStsTrackingInterface;
+class CbmTimeSlice;
 class TClonesArray;
 class TH1F;
 class TH2F;
@@ -42,8 +44,6 @@ public:
   /// \param  isMCUsed  Flag, whether MC information is available for this task
   CbmCaInputQaSts(int verbose, bool isMCUsed);
 
-  ClassDef(CbmCaInputQaSts, 0);
-
 protected:
   // ********************************************
   // ** Virtual method override from CbmQaTask **
@@ -67,13 +67,17 @@ protected:
   /// De-initializes histograms
   void DeInit();
 
-  // Checks ranges for mean and standard deviation
-  // \return  False, if variable exits the range
+  /// 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);
   }
 
+  /// \brief checks if the track at mc point passes the track selection cuts
+  /// \return true when selected
+  Bool_t IsTrackSelected(const CbmMCTrack* track, const CbmStsPoint* point) const;
+
 private:
   // ----- Data branches
   CbmStsTrackingInterface* fpDetInterface = nullptr;  ///< Instance of detector interface
@@ -90,35 +94,46 @@ private:
 
   TClonesArray* fpHitMatches = nullptr;  ///< Array of hit matches
 
+  // ----- Track selection for efficiencies and pulls
+
+  static constexpr bool kTrackSelectionPrimary {true};    // must the track come from the primary vertex
+  static constexpr double kTrackSelectionMomentum {0.1};  // track momentum GeV when it exits the tracking station
+  static constexpr double kTrackSelectionAngle {60};      // track incl. angle [grad] when it exits the tracking station
 
   // ----- Histograms
-  static constexpr double kRHitX[2] = {-50., 50};  ///< Range for hit x coordinate [cm]
-  static constexpr double kRHitY[2] = {-50., 50};  ///< Range for hit y coordinate [cm]
-  static constexpr double kRHitZ[2] = {-20., 60};  ///< Range for hit z coordinate [cm]
 
   static constexpr int kNbins  = 200;  ///< General number of bins
-  static constexpr int kNbinsZ = 800;  ///< Number of bins for z coordinate axis
+  static constexpr int kNbinsZ = 800;  ///< Number of bins for z axis in overall views
 
-  static constexpr double kRHitDx[2] = {-.005, .005};  ///< Range for hit x coordinate [cm]
-  static constexpr double kRHitDy[2] = {-.02, .02};    ///< Range for hit y coordinate [cm]
-  static constexpr double kRHitDu[2] = {-.005, .005};  ///< Range for hit u coordinate [cm]
-  static constexpr double kRHitDv[2] = {-.005, .005};  ///< Range for hit v coordinate [cm]
-  static constexpr double kRHitDt[2] = {-10., 10.};    ///< Range for hit time [ns]
+  static constexpr double kRHitDx[2] = {0., .0050};  ///< Range for hit x coordinate error [cm]
+  static constexpr double kRHitDy[2] = {0., .0200};  ///< Range for hit y coordinate error [cm]
+  static constexpr double kRHitDu[2] = {0., .0050};  ///< Range for hit u coordinate error [cm]
+  static constexpr double kRHitDv[2] = {0., .0050};  ///< Range for hit v coordinate error [cm]
+  static constexpr double kRHitDt[2] = {0., 10.};    ///< Range for hit time error [ns]
 
   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] = {-.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 kRResT[2] = {-25., 25.};  ///< Range for residual of time [ns]
 
   static constexpr double kRPull[2] = {-10., 10.};  ///< Range for pull histograms coordinate
 
+  std::vector<double> frXmin;  ///< Range for hit x coordinate [cm]
+  std::vector<double> frXmax;  ///< Range for hit x coordinate [cm]
+
+  std::vector<double> frYmin;  ///< Range for hit y coordinate [cm]
+  std::vector<double> frYmax;  ///< Range for hit y coordinate [cm]
+
+  std::vector<double> frZmin;  ///< Range for hit z coordinate [cm]
+  std::vector<double> frZmax;  ///< Range for hit z coordinate [cm]
+
   // NOTE: the last element of each vector stands for integral distribution over all stations
 
   // Hit occupancy
-  std::vector<TH2F*> fvph_hit_ypos_vs_xpos;  ///< hit ypos vs xpos in different stations
-  std::vector<TH2F*> fvph_hit_xpos_vs_zpos;  ///< hit xpos vs zpos in different stations
-  std::vector<TH2F*> fvph_hit_ypos_vs_zpos;  ///< hit ypos vs zpos in different stations
+  std::vector<TH2F*> fvph_hit_xy;  ///< hit y vs x in different stations
+  std::vector<TH2F*> fvph_hit_zx;  ///< hit x vs z in different stations
+  std::vector<TH2F*> fvph_hit_zy;  ///< hit y vs z in different stations
 
   std::vector<TH1F*> fvph_hit_station_delta_z;  ///< Difference between station and hit z positions [cm]
 
@@ -132,9 +147,9 @@ private:
   // MC points occupancy
   std::vector<TH1F*> fvph_n_points_per_hit;  ///< number of points per one hit in different stations
 
-  std::vector<TH2F*> fvph_point_ypos_vs_xpos;  ///< point ypos vs xpos in different stations
-  std::vector<TH2F*> fvph_point_xpos_vs_zpos;  ///< point xpos vs zpos in different stations
-  std::vector<TH2F*> fvph_point_ypos_vs_zpos;  ///< point ypos vs zpos in different stations
+  std::vector<TH2F*> fvph_point_xy;  ///< point y vs x in different stations
+  std::vector<TH2F*> fvph_point_zx;  ///< point x vs z in different stations
+  std::vector<TH2F*> fvph_point_zy;  ///< point y vs z in different stations
 
   std::vector<TH1F*> fvph_point_hit_delta_z;  ///< difference between z of hit and MC point in different stations
 
@@ -160,7 +175,8 @@ private:
 
   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<TH1F*> fvph_pull_u_Ndig;  ///< pull for u coordinate, depending on N digis in the cluster
+  std::vector<TH1F*> fvph_pull_v_Ndig;  ///< pull for v coordinate, 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
@@ -170,7 +186,9 @@ private:
 
   // Hit efficiencies
   std::vector<TProfile2D*> fvpe_reco_eff_vs_xy;  ///< Efficiency of hit reconstruction vs. x and y coordinates of MC
-  std::vector<TProfile*> fvpe_reco_eff_vs_r;     ///< Efficiency of hit reconstruction vs. distance from center
+  std::vector<TH1F*> fvph_reco_eff;  ///< Distribution of hit reconstruction efficiency in bins of fvpe_reco_eff_vs_xy
+
+  ClassDef(CbmCaInputQaSts, 0);
 };
 
 #endif  // CbmCaInputQaSts_h
-- 
GitLab