diff --git a/CMakeLists.txt b/CMakeLists.txt
index 817373e6531ce5836736db48a81913243859728c..53c90ccff1faffb5408e25a19eb74108e4a80b20 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -35,8 +35,8 @@ list(APPEND CMAKE_PREFIX_PATH $ENV{QN_PATH})
 set(EXTERNAL_DIR ${CMAKE_BINARY_DIR}/external)
 set(EXTERNAL_INSTALL_DIR ${CMAKE_INSTALL_PREFIX}/external)
 
-include(cmake_modules/cuts.cmake)
-include(cmake_modules/Flow.cmake)
+include(cmake_modules/Cuts.cmake)
+include(cmake_modules/QnTools.cmake)
 
 if (ROOT_FOUND)
     message(STATUS "Using ROOT: ${ROOT_VERSION} <${ROOT_CONFIG}>")
diff --git a/cmake_modules/cuts.cmake b/cmake_modules/Cuts.cmake
similarity index 100%
rename from cmake_modules/cuts.cmake
rename to cmake_modules/Cuts.cmake
diff --git a/cmake_modules/Flow.cmake b/cmake_modules/QnTools.cmake
similarity index 83%
rename from cmake_modules/Flow.cmake
rename to cmake_modules/QnTools.cmake
index f116628289ab580e3d16e5a10f2db5ddac767cfb..6fac32b7daa2fb6f85a1727995d03264cd825150 100644
--- a/cmake_modules/Flow.cmake
+++ b/cmake_modules/QnTools.cmake
@@ -1,13 +1,9 @@
-#
-#
-#
-
 set(Flow_INSTALL_DIR ${EXTERNAL_INSTALL_DIR})
 set(Flow_INCLUDE_DIR ${Flow_INSTALL_DIR}/include)
 set(Flow_LIBRARY_DIR ${Flow_INSTALL_DIR}/lib)
 
 ExternalProject_Add(Flow_Ext
-        GIT_REPOSITORY  "https://github.com/viktorklochkov/flow.git"
+        GIT_REPOSITORY  "https://github.com/viktorklochkov/QnTools.git"
         GIT_TAG         "master"
         UPDATE_DISCONNECTED ${UPDATE_DISCONNECTED}
         SOURCE_DIR      "${EXTERNAL_DIR}/Flow_src"
@@ -19,8 +15,6 @@ ExternalProject_Add(Flow_Ext
             "-DCMAKE_PREFIX_PATH=${CMAKE_PREFIX_PATH}"
 )
 
-
 list(APPEND PROJECT_DEPENDENCIES Flow_Ext)
 list(APPEND PROJECT_LINK_DIRECTORIES ${Flow_LIBRARY_DIR})
-list(APPEND PROJECT_INCLUDE_DIRECTORIES ${Flow_INCLUDE_DIR})
-
+list(APPEND PROJECT_INCLUDE_DIRECTORIES ${Flow_INCLUDE_DIR})
\ No newline at end of file
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 4e36d448324cba36c3b22ce01edcebd71029297e..22d58d67b952c847d1527fe14e68e69c63e942fa 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,8 +1,9 @@
 add_subdirectory(base)
 add_subdirectory(correct)
+add_subdirectory(correlate)
 add_subdirectory(tasks)
-#add_subdirectory(analyze)
-#add_subdirectory(tools)
+add_subdirectory(analyze)
+add_subdirectory(tools)
 
 
 
diff --git a/src/analyze/FlowCoefficients.h b/src/analyze/FlowCoefficients.h
index e7a18ddcb57dc7957e313823a1eb842a7811eeeb..3391042a4e8f294a51255b78cc577a461cc669ca 100644
--- a/src/analyze/FlowCoefficients.h
+++ b/src/analyze/FlowCoefficients.h
@@ -18,16 +18,14 @@ static Result Combine(std::vector<Result> resolutionMHParticles, const std::vect
 
 static Result FlowV1(std::vector<Result> args, const std::vector<std::string> &) {
   assert(args.size() == 2);
-
-  auto uQ = args[0];
-  auto R = args[1];
+  const auto& uQ = args[0];
+  const auto& R = args[1];
   return uQ / R;
 }
 
 static Result FlowV1MC(std::vector<Result> args, const std::vector<std::string> &) {
   assert(args.size() == 1);
-
-  auto nom1 = args[0];
+  const auto& nom1 = args[0];
   return nom1;
 }
 
@@ -43,7 +41,6 @@ static Result FlowV2Opt1(std::vector<Result> args, const std::vector<std::string
   auto R1 = args[1];
   auto R2 = args[2];
   return uQQ * 4 / (R1 * R2);
-
 }
 
 static Result FlowV2Opt2(std::vector<Result> args, const std::vector<std::string> &argNames) {
diff --git a/src/analyze/Resolution.h b/src/analyze/Resolution.h
index eafede168a46b79c3e8c5616518d723fab1f59c2..1cbf37339c6169de6e04b9286063ba0d13b3e212 100644
--- a/src/analyze/Resolution.h
+++ b/src/analyze/Resolution.h
@@ -3,108 +3,110 @@
 
 #include "Utils.h"
 
-static Result Resolution3Sub(std::vector<Result> args, const std::vector<std::string> &) {
+[[maybe_unused]] static Result Resolution3Sub(std::vector<Result> args, const std::vector<std::string> &) {
   assert(args.size() == 3);
 
-  auto &nom1 = args[0];
-  auto &nom2 = args[1];
-  auto &denom = args[2];
+  const auto& nom1 = args[0];
+  const auto& nom2 = args[1];
+  const auto& denom = args[2];
   return Qn::Sqrt(nom1 * nom2 / denom);
 }
 
-static Result ResolutionMC(std::vector<Result> args, const std::vector<std::string> &) {
+[[maybe_unused]] static Result ResolutionMC(std::vector<Result> args, const std::vector<std::string> &) {
   assert(args.size() == 1);
-
-  auto nom1 = args[0];
+  const auto& nom1= args[0];
   return nom1;
 }
 
-struct ResolutionTrack {
-  float yMin{0};
-  float yMax{0};
+ class ResolutionTrack {
+ public:
 
   ResolutionTrack () = delete;
-  ResolutionTrack (float _yMin, float _yMax) :
-    yMin(_yMin), yMax(_yMax) {};
+  ResolutionTrack (float y_min, float y_max) :
+    y_min_(y_min), y_max_(y_max) {};
 
-  Result operator()(std::vector<Result> args, const std::vector<std::string> &argNames) const {
+  Result operator()(std::vector<Result> args, const std::vector<std::string> &) const {
     assert(args.size() == 3);
 
     // multiplication of the two observables is undefined
     // -> convert to REFERENCE
-    Result nom1 = args[0];
+    auto& nom1 = args[0];
     nom1 = SetReference(nom1);
-    nom1 = RebinRapidity(nom1, {yMin, yMax}).Projection({"Centrality"});
-    Result nom2 = args[1];
+    nom1 = RebinRapidity(nom1, {y_min_, y_max_}).Projection({"Centrality"});
+    auto& nom2 = args[1];
     nom2 = SetReference(nom2);
-    nom2 = RebinRapidity(nom2, {yMin, yMax}).Projection({"Centrality"});
-    Result denom = args[2];
+    nom2 = RebinRapidity(nom2, {y_min_, y_max_}).Projection({"Centrality"});
+    const auto& denom = args[2];
     return Qn::Sqrt(nom1 * nom2 / denom);
   }
-};
+ private:
+   float y_min_{0};
+   float y_max_{0};
+ };
 
-struct Resolution4S {
-  float yMin{0};
-  float yMax{0};
+class Resolution4S {
+ public:
 
   Resolution4S () = delete;
-  Resolution4S (float _yMin, float _yMax) :
-    yMin(_yMin), yMax(_yMax) {};
+  Resolution4S (float y_min, float y_max) :
+    y_min_(y_min), y_max_(y_max) {};
 
   Result operator()(std::vector<Result> args, const std::vector<std::string> &) {
 
     if (args.size() == 3) {
-      Result nom = args[0];
-      Result RT = args[1];
-      Result denom = args[2];
+      auto& nom = args[0];
+      auto& RT = args[1];
+      auto& denom = args[2];
       denom = SetReference(denom);
-      denom = RebinRapidity(denom, {yMin, yMax}).Projection({"Centrality"});
+      denom = RebinRapidity(denom, {y_min_, y_max_}).Projection({"Centrality"});
       nom = nom * (-1); //TODO
       return nom * RT / denom;
     } else if (args.size() == 2) {
-      Result nom = args[0];
-      Result RT = args[1];
+      auto& nom = args[0];
+      auto& RT = args[1];
       nom = SetReference(nom);
-      nom = RebinRapidity(nom, {yMin, yMax}).Projection({"Centrality"});
+      nom = RebinRapidity(nom, {y_min_, y_max_}).Projection({"Centrality"});
       nom = nom * (-1); //TODO
       return nom / RT;
     }
-
     assert(0);
-    return {};
   }
+ private:
+  float y_min_{0};
+  float y_max_{0};
 };
 
-struct ResolutionMH {
-  float yMin{0};
-  float yMax{0};
-
+class ResolutionMH {
+ public:
   ResolutionMH () = delete;
-  ResolutionMH (float _yMin, float _yMax) :
-    yMin(_yMin), yMax(_yMax) {};
+  ResolutionMH (float y_min, float y_max) :
+    y_min_(y_min), y_max_(y_max) {};
 
   Result operator()(std::vector<Result> args, const std::vector<std::string> &argNames) const {
     assert(args.size() == 3);
 
-    auto nom1 = args[0];
-    auto nom2 = args[1];
-    auto denom = args[2];
+    auto& nom1 = args[0];
+    auto& nom2 = args[1];
+    auto& denom = args[2];
 
     // flip sign
-    if (argNames[0].find("X2YY") < argNames[0].size()) {
+    if (argNames[0].find("Q2x_Q1y_Q1y") < argNames[0].size()) {
       nom1 = nom1 * (-1);
     }
 
     // flip sign
-    if (argNames[2].find("X2YY") < argNames[2].size()) {
+    if (argNames[2].find("Q2x_Q1y_Q1y") < argNames[2].size()) {
       denom = denom * (-1);
     }
 
-    nom1 = RebinRapidity(nom1, {yMin, yMax}).Projection({"Centrality"});
-    denom = RebinRapidity(denom, {yMin, yMax}).Projection({"Centrality"});
+    nom1 = RebinRapidity(nom1, {y_min_, y_max_}).Projection({"Centrality"});
+    denom = RebinRapidity(denom, {y_min_, y_max_}).Projection({"Centrality"});
 
     return Qn::Sqrt(nom1 * nom2 / denom);
   }
+ private:
+  float y_min_{0};
+  float y_max_{0};
 };
 
 #endif //FLOW_PROCESSING__RESOLUTION_H_
diff --git a/src/analyze/Utils.h b/src/analyze/Utils.h
index 33f6d655374f0220258a5c3585bd185726fcdea5..8718cce7bad8076e814f440dcd5613421a0d386d 100644
--- a/src/analyze/Utils.h
+++ b/src/analyze/Utils.h
@@ -8,10 +8,12 @@
 #include <functional>
 
 #include <Rtypes.h>
-#include <DataContainer.h>
 #include <TDirectory.h>
 #include <TKey.h>
 #include <TPaveText.h>
+#include <TFile.h>
+
+#include <QnTools/DataContainer.hpp>
 
 #define FMT_HEADER_ONLY
 #include "fmt/core.h"
@@ -22,7 +24,7 @@ using Result = Qn::DataContainerStats;
 using ResultPtr = std::shared_ptr<Result>;
 
 static TPaveText *MakePaveText(const std::vector<std::string> &lines, std::vector<Double_t> position) {
-  auto pave = new TPaveText;
+  auto *pave = new TPaveText;
   pave->SetX1NDC(position.at(0));
   pave->SetY1NDC(position.at(1));
   pave->SetX2NDC(position.at(2));
@@ -54,7 +56,7 @@ struct Selection {
 
     std::string axis_name{};
     std::vector<std::string> axis_names{};
-    for (const Qn::Axis &ax : r.GetAxes()) {
+    for (const Qn::AxisD &ax : r.GetAxes()) {
       if (ax.Name().find(axis_) < ax.Name().size()) {
         axis_name = ax.Name();
       }
@@ -65,13 +67,13 @@ struct Selection {
     if (axis_name.empty()) {
       throw std::runtime_error("No axis found");
     }
-    Qn::Axis newAxis(axis_name, 1, min_, max_);
+    Qn::AxisD newAxis(axis_name, 1, min_, max_);
     return r.Rebin(newAxis).Projection(axis_names);
   }
 };
 
 static Result Select(const Result &r, std::string axis, float min, float max) {
-  for (const Qn::Axis &ax : r.GetAxes()) {
+  for (const Qn::AxisD &ax : r.GetAxes()) {
     if (ax.Name().find(axis) < ax.Name().size()) {
       axis = ax.Name();
       break;
@@ -80,15 +82,15 @@ static Result Select(const Result &r, std::string axis, float min, float max) {
   if (axis.empty()) {
     throw std::runtime_error("No axis found");
   }
-  Qn::Axis newAxis(axis, 1, min, max);
+  Qn::AxisD newAxis(axis, 1, min, max);
   return r.Select(newAxis);
 }
 
 
-static Result RebinRapidity(const Result &r, const std::vector<float> &binEdges) {
+static Result RebinRapidity(const Result &r, const std::vector<double> &binEdges) {
   std::string axRapidityName;
   // find rapidity axis
-  for (const Qn::Axis &ax : r.GetAxes()) {
+  for (const Qn::AxisD &ax : r.GetAxes()) {
     if (ax.Name().find("rapidity") < ax.Name().size()) {
       axRapidityName = ax.Name();
       break;
@@ -99,7 +101,7 @@ static Result RebinRapidity(const Result &r, const std::vector<float> &binEdges)
     throw std::runtime_error("No rapidity axis found");
   }
 
-  Qn::Axis newRapidityAxis(axRapidityName, binEdges);
+  Qn::AxisD newRapidityAxis(axRapidityName, binEdges);
 
   return r.Rebin(newRapidityAxis);
 }
@@ -108,7 +110,7 @@ static Result SetReference(const Result &r) {
   Result rNew{r};
 
   for (Qn::Stats &stats : rNew) {
-    stats.SetStatus(Qn::Stats::Status::REFERENCE);
+    stats.SetWeights(Qn::Stats::Weights::REFERENCE);
   }
 
   return rNew;
@@ -127,9 +129,9 @@ static bool GetResource(TDirectory *d, const std::string &name, std::shared_ptr<
 }
 
 /**
- * 
- * @param mg 
- * @param dX 
+ *
+ * @param mg
+ * @param dX
  */
 static void ShiftGraphsX(TMultiGraph &mg, double dX = 0.0) {
 
@@ -216,9 +218,9 @@ struct ProfileExporter {
     return *this;
   }
 
-  std::shared_ptr<Qn::Axis> rebinAxisPtr{};
-  ProfileExporter &Rebin(const Qn::Axis &_axis) {
-    rebinAxisPtr = std::make_shared<Qn::Axis>(_axis);
+  std::shared_ptr<Qn::AxisD> rebinAxisPtr{};
+  ProfileExporter &Rebin(const Qn::AxisD &_axis) {
+    rebinAxisPtr = std::make_shared<Qn::AxisD>(_axis);
     return *this;
   }
   ProfileExporter &Rebin() {
@@ -290,7 +292,7 @@ struct ProfileExporter {
     if (result.GetAxes().size() == 2) {
       std::string axisname{"Centrality"};
       auto profileMultigraph = new TMultiGraph();
-      Qn::Axis axis;
+      Qn::AxisD axis;
       try { axis = result.GetAxis(axisname); } // TODO
       catch (std::exception &) {
         // TODO return?
diff --git a/src/base/GlobalConfig.cpp b/src/base/GlobalConfig.cpp
index 1c52069a41b2c196846e51acf284f8f39001bb78..2520f3a0185382e1dd8b6c93d78326cf8f782358 100644
--- a/src/base/GlobalConfig.cpp
+++ b/src/base/GlobalConfig.cpp
@@ -22,6 +22,10 @@ void Qn::GlobalConfig::Print() const {
   for(const auto& det : channel_qvectors_){
     det.Print();
   }
+  if(is_simulation_){
+    std::cout << " ***** Psi_RP ***** :" << std::endl;
+    psi_qvector_.Print();
+  }
 }
 
 void Qn::GlobalConfig::DieWith(const std::string& message) {
diff --git a/src/base/GlobalConfig.h b/src/base/GlobalConfig.h
index ccab2a1dc7e98dfb9ca40ed94b590fb8b11153a7..49dc4736d01ac27ad36b75f8f652ed2911f72081 100644
--- a/src/base/GlobalConfig.h
+++ b/src/base/GlobalConfig.h
@@ -9,17 +9,16 @@
 #include "TString.h"
 #include "TError.h"
 
-#include "QvectorTracksConfig.h"
-#include "QvectorChannelsConfig.h"
-#include "Axis.h"
-
-//#include "BranchReader.h"
+#include "QnTools/Axis.hpp"
 
 #include "AnalysisTree/Cuts.h"
 #include "AnalysisTree/EventHeader.h"
 #include "AnalysisTree/BranchConfig.h"
 #include "AnalysisTree/Variable.h"
 
+#include "QvectorTracksConfig.h"
+#include "QvectorChannelsConfig.h"
+
 #define gConfig Qn::GlobalConfig::GetInstance()
 
 namespace Qn{
@@ -72,27 +71,47 @@ class GlobalConfig : public TObject {
   void SetEventCuts(AnalysisTree::Cuts* cuts) { event_cuts_ = cuts; }
   AnalysisTree::Cuts* GetEventCuts() { return event_cuts_; }
 
-  void AddCorrectionAxis(const Qn::Axis<double>& axis) {
+  void AddCorrectionAxis(const Qn::AxisD& axis) {
     correction_axes_.emplace_back(axis);
   }
 
-  const std::vector <Qn::Axis<double>>& GetCorrectionAxes() const { return correction_axes_; }
+  const std::vector <Qn::AxisD>& GetCorrectionAxes() const { return correction_axes_; }
   const QvectorChannelsConfig& GetPsiQvector() const {return psi_qvector_;}
+
+  const std::vector<std::string>& GetCorrelationNames() const { return correlation_names_; }
+
+  std::string ConstructCorrelationName(const std::vector<std::string>& q_vectors, const std::string& type) const {
+    std::string name;
+    for(const auto& qn : q_vectors){
+      name += qn + "_";
+    }
+    name += type;
+    return name;
+  }
+  void SetNSamples(int n_samples) {
+    n_samples_ = n_samples;
+  }
+
  protected:
 
   void DieWith(const std::string &message);
 
+  // Correction parameters
   AnalysisTree::Cuts* event_cuts_{nullptr};  //!
   std::vector <QvectorTracksConfig> track_qvectors_{};
   std::vector <QvectorChannelsConfig> channel_qvectors_{};
   QvectorChannelsConfig psi_qvector_;
+//  std::vector<QvectorConfig*> u_vectors_{};
+//  std::vector<QvectorConfig*> q_vectors_;
   std::vector <AnalysisTree::Variable> event_vars_;
-  std::vector <Qn::Axis<double>> correction_axes_{};
+  std::vector <Qn::AxisD> correction_axes_{};
 
-  // sampling
+  // Correlation parameters
   std::string sampling_mode_{""};
   int n_samples_{50};
+  std::vector<std::string> correlation_names_{};
 
+  // Global parameters
   bool is_simulation_{false};
 
  private:
diff --git a/src/base/QvectorConfig.h b/src/base/QvectorConfig.h
index 01297aa22b054fee3220aae1e75060f7b3787f81..7a316aaf16dce07686883f18863f845ea9d93855 100644
--- a/src/base/QvectorConfig.h
+++ b/src/base/QvectorConfig.h
@@ -72,7 +72,7 @@ public:
 
 protected:
 
-  std::string name_{""};          ///<  Name of the Q-vector
+  std::string name_{""};    ///<  Name of the Q-vector
 
   std::string branch_name_{""};
   int branch_id_{-999};
diff --git a/src/base/QvectorTracksConfig.h b/src/base/QvectorTracksConfig.h
index 485fe5e48af7166414dd4ee614688fe8df311b0b..9856435497de988f2768fcb7049d86794a746c5e 100644
--- a/src/base/QvectorTracksConfig.h
+++ b/src/base/QvectorTracksConfig.h
@@ -9,7 +9,7 @@
 
 #include "TTreeReader.h"
 
-#include "Axis.h"
+#include "QnTools/Axis.hpp"
 
 #include "QvectorConfig.h"
 
@@ -22,13 +22,13 @@ public:
   QvectorTracksConfig() = default;
   QvectorTracksConfig(const QvectorTracksConfig& ) = default;
   QvectorTracksConfig(QvectorTracksConfig&& ) = default;
-  QvectorTracksConfig(const std::string& name, const std::string& branch, const std::vector<Qn::Axis<double>>& axes) :
+  QvectorTracksConfig(const std::string& name, const std::string& branch, const std::vector<Qn::AxisD>& axes) :
   QvectorConfig((const std::string&&) name, branch),
   axes_(axes) {}
 
-  void SetAxes(const std::vector<Qn::Axis<double>>& axes) { axes_ = axes; }
-  const std::vector<Qn::Axis<double>>& GetAxes() const { return axes_; }
-  std::vector<Qn::Axis<double>>& GetAxes() { return axes_; }
+  void SetAxes(const std::vector<Qn::AxisD>& axes) { axes_ = axes; }
+  const std::vector<Qn::AxisD>& GetAxes() const { return axes_; }
+  std::vector<Qn::AxisD>& GetAxes() { return axes_; }
 
   void AddCut(const std::vector<std::string>& names, const std::function<bool(double)>& lambda) {
     cuts_.emplace_back(names, lambda);
@@ -41,7 +41,7 @@ public:
 
  protected:
 
-  std::vector<Qn::Axis<double>> axes_{};
+  std::vector<Qn::AxisD> axes_{};
   std::vector<std::pair<std::vector<std::string>, std::function<bool(double)>>> cuts_{};  //!
   bool is_sim_{false};
 };
diff --git a/src/correct/AnalysisTreeVarManager.cpp b/src/correct/AnalysisTreeVarManager.cpp
index 4b0c57e210d7bc8ce6328d0c96f1b9041ef7690c..13e622c4ef0d32871fe93a344f8196793c3aa5ba 100644
--- a/src/correct/AnalysisTreeVarManager.cpp
+++ b/src/correct/AnalysisTreeVarManager.cpp
@@ -1,7 +1,8 @@
 #include <TMath.h>
-#include <TFile.h>
+
+#include "base/GlobalConfig.h"
+
 #include "AnalysisTreeVarManager.h"
-#include "src/base/GlobalConfig.h"
 
 void AnalysisTreeVarManager::FillChannelDetector(const AnalysisTree::ModuleDetector &det, 
                                                  const AnalysisTree::DataHeader &info,
@@ -49,7 +50,6 @@ void AnalysisTreeVarManager::InitializeCorrectionManager(Qn::CorrectionManager &
       var.Init(*config);
       manager.AddVariable(var.GetName(), var.GetId(), var.GetSize());
       ivar += var.GetSize();
-//      var.Print();
     }
   }
 
diff --git a/src/correct/AnalysisTreeVarManager.h b/src/correct/AnalysisTreeVarManager.h
index 8767b80bdb27bebed336f1f95cb926ebd6458119..62ec1c64615f7e3874330becaa0a18f3354934e6 100644
--- a/src/correct/AnalysisTreeVarManager.h
+++ b/src/correct/AnalysisTreeVarManager.h
@@ -3,7 +3,7 @@
 
 #include <list>
 
-#include "CorrectionManager.h"
+#include "QnTools/CorrectionManager.hpp"
 
 #include "AnalysisTree/Configuration.h"
 #include "AnalysisTree/EventHeader.h"
@@ -12,8 +12,8 @@
 #include "AnalysisTree/DataHeader.h"
 #include "AnalysisTree/Variable.h"
 
-#include "src/base/QvectorTracksConfig.h"
-#include "src/base/QvectorChannelsConfig.h"
+#include "base/QvectorTracksConfig.h"
+#include "base/QvectorChannelsConfig.h"
 
 #define VAR AnalysisTreeVarManager
 #define gVarManager VAR::GetInstance()
diff --git a/src/correct/BranchReader.h b/src/correct/BranchReader.h
index 9ae6d2d3acc9601da68211d0b89de69f041a35a4..464364eeebfa24e9f2f2212ed9f37172721cda2c 100644
--- a/src/correct/BranchReader.h
+++ b/src/correct/BranchReader.h
@@ -12,7 +12,7 @@ class BranchReader {
  public:
 
   BranchReader() = default;
-  BranchReader(const std::string& name, AnalysisTree::Cuts* cuts= nullptr) : name_(name), quality_cuts_(cuts) {};
+  explicit BranchReader(const std::string& name, AnalysisTree::Cuts* cuts= nullptr) : name_(name), quality_cuts_(cuts) {};
 
   TTreeReaderValue<AnalysisTree::Particles>* GetReader() const {
     return reader_;
diff --git a/src/correct/CMakeLists.txt b/src/correct/CMakeLists.txt
index 6df83629beb6b5c39780c553710a90dc139fd076..eaa76034b0232499b5cdbd874a84b113a3580f39 100644
--- a/src/correct/CMakeLists.txt
+++ b/src/correct/CMakeLists.txt
@@ -9,7 +9,7 @@ string(REPLACE ".cpp" ".h" HEADERS "${SOURCES}")
 add_library(FlowCorrect SHARED ${SOURCES} G__FlowCorrect.cxx)
 link_directories(${PROJECT_LINK_DIRECTORIES})
 add_dependencies(FlowCorrect ${PROJECT_DEPENDENCIES} FlowBase)
-include_directories(${PROJECT_INCLUDE_DIRECTORIES} ${CMAKE_CURRENT_SOURCE_DIR})
+include_directories(${PROJECT_INCLUDE_DIRECTORIES} ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_SOURCE_DIR}/src)
 
 ROOT_GENERATE_DICTIONARY(G__FlowCorrect ${HEADERS} LINKDEF FlowCorrectLinkDef.h)
 target_link_libraries(FlowCorrect FlowBase AnalysisTreeBase AnalysisTreeInfra ${ROOT_LIBRARIES} )
diff --git a/src/correct/CorrectionTask.cpp b/src/correct/CorrectionTask.cpp
index 2e3486dd9c3ad72a30868a9e7a0a3175541740ef..c7e31dec35cb2267fe106e9f0d5227a21dee220d 100644
--- a/src/correct/CorrectionTask.cpp
+++ b/src/correct/CorrectionTask.cpp
@@ -5,7 +5,7 @@
 
 #include "AnalysisTree/TreeReader.h"
 
-#include "src/base/GlobalConfig.h"
+#include "base/GlobalConfig.h"
 #include "Utils.h"
 
 namespace Qn {
@@ -34,9 +34,6 @@ void CorrectionTask::Initialize()
   psd_ = new TTreeReaderValue<AnalysisTree::ModuleDetector>(tree_reader_, channels_name_.c_str());
   ana_event_ = new TTreeReaderValue<AnalysisTree::Container>(tree_reader_, ana_event_name_.c_str());
 
-//  config_->Print();
-//  static_info_->Print();
-
   out_file_->cd();
   out_tree_ = new TTree("tree", "tree");
   manager_.SetCalibrationInputFileName(in_calibration_file_name_);
@@ -77,7 +74,7 @@ void CorrectionTask::Initialize()
 // Add Q-vectors of tracking detectors
   for(auto& qvec : gConfig->GetQvectorsConfig()) {
     const string& name = qvec.GetName();
-    manager_.AddDetector(name, DetectorType::TRACK, qvec.GetPhiVar().GetName(), qvec.GetWeightVar().GetName(), qvec.GetAxes(), {1,2});
+    manager_.AddDetector(name, DetectorType::TRACK, qvec.GetPhiVar().GetName(), qvec.GetWeightVar().GetName(), qvec.GetAxes(), {1,2}, QVector::Normalization::M);
     Utils::SetCorrectionSteps(manager_, qvec);
     manager_.AddCutOnDetector(name, {(qvec.GetBranchName()+"_Filled").c_str()}, [](const double is){ return is>0.; }, "filled" );
     for(const auto& cut : qvec.GetCuts()) {
@@ -89,13 +86,13 @@ void CorrectionTask::Initialize()
 // Add Q-vectors of channelized detectors
   for(const auto& qvec : gConfig->GetChannelConfig()) {
     const string name = qvec.GetName();
-    manager_.AddDetector(name, DetectorType::CHANNEL, qvec.GetPhiVar().GetName(), qvec.GetWeightVar().GetName(), {}, {1}, QVector::Normalization::MAGNITUDE);
+    manager_.AddDetector(name, DetectorType::CHANNEL, qvec.GetPhiVar().GetName(), qvec.GetWeightVar().GetName(), {}, {1});
     Utils::SetCorrectionSteps(manager_, qvec);
   }
 
   if(gConfig->IsSimulation()){
     const auto & qvec = gConfig->GetPsiQvector();
-    manager_.AddDetector(qvec.GetName(), DetectorType::CHANNEL, qvec.GetPhiVar().GetName(), qvec.GetWeightVar().GetName(), {}, {1,2}, QVector::Normalization::MAGNITUDE);
+    manager_.AddDetector(qvec.GetName(), DetectorType::CHANNEL, qvec.GetPhiVar().GetName(), qvec.GetWeightVar().GetName(), {}, {1,2});
     Utils::SetCorrectionSteps(manager_, qvec);
   }
 
@@ -109,7 +106,6 @@ void CorrectionTask::Initialize()
     gConfig->GetEventCuts()->Init(*config_);
 
   gConfig->Print();
-
 }
 
 void CorrectionTask::Run() {
@@ -193,7 +189,7 @@ void CorrectionTask::AddQAHisto()   //TODO make it automatic
   manager_.AddHisto1D("psd3", {"cbm44_3_phi", 500, -4, 4}, "Ones");
 
   if(is_sim_) {
-//    manager_.AddHisto1D("psi", {"SimEventHeader_psi_RP", 100, - 4, 4}, "Ones");
+    manager_.AddHisto1D("psi", {"SimEventHeader_psi_RP", 100, 0, 7}, "Ones");
     manager_.AddEventHisto1D({"b", 100, 0, 20});
   }
   manager_.AddEventHisto1D({"Centrality", 20, 0, 100});
@@ -212,7 +208,6 @@ void CorrectionTask::Finalize() {
   gConfig->Write("Config");
 
   out_file_->Close();
-
 }
 
 
diff --git a/src/correct/CorrectionTask.h b/src/correct/CorrectionTask.h
index f05bc3e1ac1901cc27f5a89fad78229e6056355d..f0ca884eaa598bd97ec0aa6d026464c4503942a8 100644
--- a/src/correct/CorrectionTask.h
+++ b/src/correct/CorrectionTask.h
@@ -11,12 +11,7 @@
 #include "TTree.h"
 #include "TTreeReader.h"
 
-#include "CorrectionManager.h"
-
-#include "src/base/QvectorConfig.h"
-
-#include "AnalysisTreeVarManager.h"
-#include "BranchReader.h"
+#include "QnTools/CorrectionManager.hpp"
 
 #include "AnalysisTree/Detector.h"
 #include "AnalysisTree/Track.h"
@@ -25,6 +20,12 @@
 #include "AnalysisTree/Variable.h"
 #include "AnalysisTree/DataHeader.h"
 
+#include "base/QvectorConfig.h"
+
+#include "AnalysisTreeVarManager.h"
+#include "BranchReader.h"
+
+
 #define VAR AnalysisTreeVarManager
 
 namespace Qn {
diff --git a/src/correct/Utils.h b/src/correct/Utils.h
index f4cebbbfc8f4770f04d5b0b6bcf662b2cdaa79b4..a3a61bcd8bb9c1c4c3b8994b3d85595127e2a937 100644
--- a/src/correct/Utils.h
+++ b/src/correct/Utils.h
@@ -1,9 +1,9 @@
 #ifndef DATATREEFLOW_UTILS_H
 #define DATATREEFLOW_UTILS_H
 
-#include "CorrectionManager.h"
+#include "QnTools/CorrectionManager.hpp"
 
-#include "src/base/QvectorConfig.h"
+#include "base/QvectorConfig.h"
 
 namespace Utils {
 
diff --git a/src/correlate/CMakeLists.txt b/src/correlate/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..033a8e4afb71b0af3ae643e5b2c20f2744417716
--- /dev/null
+++ b/src/correlate/CMakeLists.txt
@@ -0,0 +1,41 @@
+set(SOURCES CorrelationTask.cpp)
+
+set(HEADERS
+    CorrelationTask.h
+    Utils.h
+    Functions.h
+        )
+
+add_library(FlowCorrelate SHARED ${SOURCES} G__FlowCorrelate.cxx)
+link_directories(${PROJECT_LINK_DIRECTORIES})
+add_dependencies(FlowCorrelate ${PROJECT_DEPENDENCIES})
+include_directories(${PROJECT_INCLUDE_DIRECTORIES} ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_SOURCE_DIR}/src)
+
+ROOT_GENERATE_DICTIONARY(G__FlowCorrelate ${HEADERS} LINKDEF FlowCorrelateLinkDef.h)
+target_link_libraries(FlowCorrelate ${ROOT_LIBRARIES})
+
+install(TARGETS FlowCorrelate EXPORT FlowTargets
+        LIBRARY DESTINATION lib
+        ARCHIVE DESTINATION lib
+        RUNTIME DESTINATION bin
+        INCLUDES DESTINATION include
+        )
+
+install(
+        FILES
+        ${HEADERS}
+        DESTINATION
+        include
+        COMPONENT
+        Devel
+)
+
+set(PCM_FILE_NAME libFlowCorrelate)
+
+install(
+        FILES
+        "${CMAKE_CURRENT_BINARY_DIR}/${PCM_FILE_NAME}_rdict.pcm"
+        DESTINATION
+        lib
+        OPTIONAL
+)
\ No newline at end of file
diff --git a/src/correlate/CorrelationTask.cpp b/src/correlate/CorrelationTask.cpp
index 6cde5b2aef23aaf40ccc6c1ad6b2a2ea83b803d5..cb0a37033cffc4564f6340558cbd1d3f37ce9a45 100644
--- a/src/correlate/CorrelationTask.cpp
+++ b/src/correlate/CorrelationTask.cpp
@@ -1,13 +1,11 @@
 #include <TFileCollection.h>
-#include "TCanvas.h"
 #include "TTreeReaderValue.h"
 
-#include "src/base/GlobalConfig.h"
+#include <QnTools/QnDataFrame.hpp>
 
 #include "CorrelationTask.h"
 #include "Utils.h"
 
-typedef std::map<std::string, ROOT::RDF::RResultPtr<Qn::DataContainerStats>> ResultMap;
 constexpr auto obs = Qn::Stats::Weights::OBSERVABLE;
 constexpr auto ref = Qn::Stats::Weights::REFERENCE;
 
@@ -17,7 +15,7 @@ CorrelationTask::CorrelationTask(const std::string& file, const std::string& tre
   in_tree_(in_file_->Get<TTree>(treename.c_str())),
   reader_(TTreeReader(in_tree_)),
   df_(ROOT::RDataFrame(*in_tree_)),
-  df_samples_( df_.Define("Samples", Qn::Correlation::ReSampler(config_->GetNSamples()), {}) ),
+  dfs_(Qn::Correlation::Resample(df_, config_->GetNSamples())),
   event_(config_->GetCorrectionAxes().at(0)) // NOTE
 {}
 
@@ -34,28 +32,28 @@ void CorrelationTask::FillCorrelations() {
       AddQ1Q1Correlations(qn, config_->GetPsiQvector());
     }
   }
-  // Q1_channel * Q1_track
-  for (const auto& qn_ch : config_->GetChannelConfig()) {
-    for (const auto& qn_tr : config_->GetQvectorsConfig()) {
-      AddQ1Q1Correlations(qn_tr, qn_ch);
-    }
-  }
-  // Q1_channel * Q1_channel
-  for (const auto& qn_ch1 : config_->GetChannelConfig()) {
-    for (const auto& qn_ch2 : config_->GetChannelConfig()) {
-      if( qn_ch1.GetName() != qn_ch2.GetName() )
-      AddQ1Q1Correlations(qn_ch1, qn_ch2);
-    }
-  }
-  // Q2_track * Q1_channel * Q1_channel
-  for (const auto& qn_tr : config_->GetQvectorsConfig()) {
-    for (const auto& qn_ch1 : config_->GetChannelConfig()) {
-      for (const auto& qn_ch2 : config_->GetChannelConfig()) {
-        if(qn_ch1.GetName() != qn_ch2.GetName())
-          AddQ2Q1Q1Correlations(qn_tr, qn_ch1, qn_ch2);
-      }
-    }
-  }
+//  // Q1_channel * Q1_track
+//  for (const auto& qn_ch : config_->GetChannelConfig()) {
+//    for (const auto& qn_tr : config_->GetQvectorsConfig()) {
+//      AddQ1Q1Correlations(qn_tr, qn_ch);
+//    }
+//  }
+//  // Q1_channel * Q1_channel
+//  for (const auto& qn_ch1 : config_->GetChannelConfig()) {
+//    for (const auto& qn_ch2 : config_->GetChannelConfig()) {
+//      if( qn_ch1.GetName() != qn_ch2.GetName() )
+//      AddQ1Q1Correlations(qn_ch1, qn_ch2);
+//    }
+//  }
+//  // Q2_track * Q1_channel * Q1_channel
+//  for (const auto& qn_tr : config_->GetQvectorsConfig()) {
+//    for (const auto& qn_ch1 : config_->GetChannelConfig()) {
+//      for (const auto& qn_ch2 : config_->GetChannelConfig()) {
+//        if(qn_ch1.GetName() != qn_ch2.GetName())
+//          AddQ2Q1Q1Correlations(qn_tr, qn_ch1, qn_ch2);
+//      }
+//    }
+//  }
 }
 
 void CorrelationTask::Run() {
@@ -65,8 +63,13 @@ void CorrelationTask::Run() {
 
   auto* out_file = new TFile("correlationout.root", "RECREATE");
   out_file->cd();
+
+  std::vector<Qn::DataContainerStats> stats;
   for (auto &correlation : correlations_) {
-    correlation.second.GetValue().Write(correlation.first.data());
+    stats.push_back(correlation.second.GetValue().GetDataContainer());
+  }
+  for (auto &correlation : correlations_) {
+    correlation.second->Write();
   }
   out_file->Close();
 
@@ -80,10 +83,17 @@ void CorrelationTask::AddQ1Q1Correlations(const Qn::QvectorConfig& a, const Qn::
     if(verbosity_ >= 1){
       std::cout << "Adding correlation: " << a_name << " " << b_name << " " << cor.first << std::endl;
     }
-    auto q1q1 = Qn::Correlation::MakeCorrelation(cor.first, cor.second, Qn::Correlation::MakeAxes(event_))
-      .SetInputNames(a_name, b_name).SetWeights(obs, ref).BookMe(df_samples_, reader_, config_->GetNSamples());
+    const auto name = a_name + "_" + b_name + "_" + cor.first;
+    auto q11 = Qn::EventAverage(Qn::Correlation::Correlation(
+      name,
+      cor.second,
+      {a_name, b_name},
+      {obs, ref},
+      Qn::EventAxes(event_),
+      config_->GetNSamples()))
+      .BookMe(dfs_);
 
-    correlations_.emplace(a_name + "_" + b_name + "_" + cor.first, q1q1);
+    correlations_.emplace(name, q11);
   }
 }
 
@@ -95,9 +105,17 @@ void CorrelationTask::AddQ2Q1Q1Correlations(const Qn::QvectorTracksConfig& t, co
     if(verbosity_ >= 1){
       std::cout << "Adding correlation: " << t_name << " " << a_name << " " << b_name << " " << cor.first << std::endl;
     }
-    auto q2q1q1 = Qn::Correlation::MakeCorrelation(cor.first, cor.second, Qn::Correlation::MakeAxes(event_))
-      .SetInputNames(t_name, a_name, b_name).SetWeights(ref, ref, ref).BookMe(df_samples_, reader_, config_->GetNSamples());
 
-    correlations_.emplace(t.GetName() + "_" + a_name + "_" + b_name + "_" + cor.first, q2q1q1);
+    const auto name = t_name + "_" + a_name + "_" + b_name + "_" + cor.first;
+    auto q2q1q1 = Qn::EventAverage(Qn::Correlation::Correlation(
+      name,
+      cor.second,
+      {t_name, a_name, b_name},
+      {obs, ref, ref},
+      Qn::EventAxes(event_),
+      config_->GetNSamples()))
+      .BookMe(dfs_);
+
+    correlations_.emplace(name, q2q1q1);
   }
 }
diff --git a/src/correlate/CorrelationTask.h b/src/correlate/CorrelationTask.h
index b3ef6e03eb5c9a8ff649ea71026bfe7da5c1cd79..2c387b0a161e0c030fd2948f0c62e26637624e30 100644
--- a/src/correlate/CorrelationTask.h
+++ b/src/correlate/CorrelationTask.h
@@ -4,12 +4,16 @@
 #include "TChain.h"
 #include "TTreeReader.h"
 
-#include <DFCorrelation.h>
+#include "QnTools/ReSampleHelper.hpp"
+#include "QnTools/DataContainer.hpp"
+#include "QnTools/Axis.hpp"
+#include "QnTools/QnDataFrame.hpp"
 
-#include "src/base/QvectorConfig.h"
-#include "src/base/GlobalConfig.h"
-#include "Axis.h"
+#include "ROOT/RDataFrame.hxx"
+#include "ROOT/RDataSource.hxx"
 
+#include "base/QvectorConfig.h"
+#include "base/GlobalConfig.h"
 
 class CorrelationTask {
  public:
@@ -28,13 +32,13 @@ class CorrelationTask {
   TTree* in_tree_{nullptr};
   TTreeReader reader_;
   ROOT::RDataFrame df_;
-  ROOT::RDF::RInterface<ROOT::Detail::RDF::RLoopManager> df_samples_;
+  ROOT::RDF::RInterface<ROOT::Detail::RDF::RLoopManager> dfs_;
   Qn::AxisD event_;
 
-  std::map<std::string, ROOT::RDF::RResultPtr<Qn::DataContainerStats>> correlations_{};
+  std::map<std::string, ROOT::RDF::RResultPtr<Qn::Correlation::CorrelationActionBase>> correlations_{};
 
   int verbosity_{2};
-  bool non_zero_only_{false};
+  bool non_zero_only_{true};
 };
 
 #endif //CORRELATION_TASK_H
diff --git a/src/correlate/FlowCorrelateLinkDef.h b/src/correlate/FlowCorrelateLinkDef.h
new file mode 100644
index 0000000000000000000000000000000000000000..08feb58f3949d6b9b69e9a9263bdb4a529fb019c
--- /dev/null
+++ b/src/correlate/FlowCorrelateLinkDef.h
@@ -0,0 +1,6 @@
+#ifndef FLOW_SRC_CORRELATE_FLOWCORRELATELINKDEF_H_
+#define FLOW_SRC_CORRELATE_FLOWCORRELATELINKDEF_H_
+
+#pragma link C++ namespace Qn;
+
+#endif //FLOW_SRC_CORRELATE_FLOWCORRELATELINKDEF_H_
diff --git a/src/correlate/Functions.h b/src/correlate/Functions.h
index ad3e4eac241fae902a9dea6737db4266fbaf95de..70b8b719800a0f5f09ede69d1083b0163f1345a4 100644
--- a/src/correlate/Functions.h
+++ b/src/correlate/Functions.h
@@ -1,104 +1,103 @@
 #ifndef FLOW_SRC_CORRELATE_FUNCTIONS_H_
 #define FLOW_SRC_CORRELATE_FUNCTIONS_H_
 
-#include <QVector.h>
-#include <Correlation.h>
+#include <QnTools/QVector.hpp>
 
 namespace Qn{
 
 // 2 particle correlations
 
-[[maybe_unused]] auto XY = [](const Qn::QVector &a, const Qn::QVector &b) {
+[[maybe_unused]] static auto XY = [](const Qn::QVector &a, const Qn::QVector &b) {
   return 2 * a.x(1) * b.y(1);
 };
-[[maybe_unused]] auto YX = [](const Qn::QVector &a, const Qn::QVector &b) {
+[[maybe_unused]] static auto YX = [](const Qn::QVector &a, const Qn::QVector &b) {
   return 2 * a.y(1) * b.x(1);
 };
-[[maybe_unused]] auto XX = [](const Qn::QVector &a, const Qn::QVector &b) {
+[[maybe_unused]] static auto XX = [](const Qn::QVector &a, const Qn::QVector &b) {
   return 2 * a.x(1) * b.x(1);
 };
-[[maybe_unused]] auto YY = [](const Qn::QVector &a, const Qn::QVector &b) {
+[[maybe_unused]] static auto YY = [](const Qn::QVector &a, const Qn::QVector &b) {
   return 2 * a.y(1) * b.y(1);
 };
-[[maybe_unused]] auto scalar = [](const Qn::QVector &a, const Qn::QVector &b) {
+[[maybe_unused]] static auto scalar = [](const Qn::QVector &a, const Qn::QVector &b) {
   return Qn::ScalarProduct(a, b, 2);
 };
 // 3 particle correlations
 
-[[maybe_unused]] auto Y2XY = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c) {
+[[maybe_unused]] static auto Y2XY = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c) {
   return a.y(2) * b.x(1) * c.y(1);
 };
-[[maybe_unused]] auto Y2YX = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c) {
+[[maybe_unused]] static auto Y2YX = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c) {
   return a.y(2) * b.y(1) * c.x(1);
 };
-[[maybe_unused]] auto X2XX = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c) {
+[[maybe_unused]] static auto X2XX = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c) {
   return a.x(2) * b.x(1) * c.x(1);
 };
-[[maybe_unused]] auto X2YY = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c) {
+[[maybe_unused]] static auto X2YY = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c) {
   return a.x(2) * b.y(1) * c.y(1);
 };
 
-[[maybe_unused]] auto X2XY = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c) {
+[[maybe_unused]] static auto X2XY = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c) {
   return a.x(2) * b.x(1) * c.y(1);
 };
-[[maybe_unused]] auto X2YX = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c) {
+[[maybe_unused]] static auto X2YX = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c) {
   return a.x(2) * b.y(1) * c.x(1);
 };
-[[maybe_unused]] auto Y2XX = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c) {
+[[maybe_unused]] static auto Y2XX = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c) {
   return a.y(2) * b.x(1) * c.x(1);
 };
-[[maybe_unused]] auto Y2YY = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c) {
+[[maybe_unused]] static auto Y2YY = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c) {
   return a.y(2) * b.y(1) * c.y(1);
 };
 
 // 4 particle correlations
 
-[[maybe_unused]] auto Y3YYY = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
+[[maybe_unused]] static auto Y3YYY = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
   return a.y(3) * b.y(1) * c.y(1) * d.y(1);
 };
-[[maybe_unused]] auto Y3XYY = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
+[[maybe_unused]] static auto Y3XYY = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
   return a.y(3) * b.x(1) * c.y(1) * d.y(1);
 };
-[[maybe_unused]] auto Y3YXY = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
+[[maybe_unused]] static auto Y3YXY = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
   return a.y(3) * b.y(1) * c.x(1) * d.y(1);
 };
-[[maybe_unused]] auto Y3YYX = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
+[[maybe_unused]] static auto Y3YYX = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
   return a.y(3) * b.y(1) * c.y(1) * d.x(1);
 };
-[[maybe_unused]] auto Y3YXX = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
+[[maybe_unused]] static auto Y3YXX = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
   return a.y(3) * b.y(1) * c.x(1) * d.x(1);
 };
-[[maybe_unused]] auto Y3XYX = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
+[[maybe_unused]] static auto Y3XYX = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
   return a.y(3) * b.x(1) * c.y(1) * d.x(1);
 };
-[[maybe_unused]] auto Y3XXY = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
+[[maybe_unused]] static auto Y3XXY = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
   return a.y(3) * b.x(1) * c.x(1) * d.y(1);
 };
-[[maybe_unused]] auto Y3XXX = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
+[[maybe_unused]] static auto Y3XXX = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
   return a.y(3) * b.x(1) * c.x(1) * d.x(1);
 };
-[[maybe_unused]] auto X3YYY = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
+[[maybe_unused]] static auto X3YYY = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
   return a.x(3) * b.y(1) * c.y(1) * d.y(1);
 };
-[[maybe_unused]] auto X3XYY = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
+[[maybe_unused]] static auto X3XYY = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
   return a.x(3) * b.x(1) * c.y(1) * d.y(1);
 };
-[[maybe_unused]] auto X3YXY = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
+[[maybe_unused]] static auto X3YXY = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
   return a.x(3) * b.y(1) * c.x(1) * d.y(1);
 };
-[[maybe_unused]] auto X3YYX = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
+[[maybe_unused]] static auto X3YYX = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
   return a.x(3) * b.y(1) * c.y(1) * d.x(1);
 };
-[[maybe_unused]] auto X3YXX = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
+[[maybe_unused]] static auto X3YXX = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
   return a.x(3) * b.y(1) * c.x(1) * d.x(1);
 };
-[[maybe_unused]] auto X3XYX = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
+[[maybe_unused]] static auto X3XYX = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
   return a.x(3) * b.x(1) * c.y(1) * d.x(1);
 };
-[[maybe_unused]] auto X3XXY = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
+[[maybe_unused]] static auto X3XXY = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
   return a.x(3) * b.x(1) * c.x(1) * d.y(1);
 };
-[[maybe_unused]] auto X3XXX = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
+[[maybe_unused]] static auto X3XXX = [](const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d) {
   return a.x(3) * b.x(1) * c.x(1) * d.x(1);
 };
 
diff --git a/src/correlate/Utils.h b/src/correlate/Utils.h
index 5bb93ee578d58fc2d40604969155f16cf9061a75..19b2b9780a8dce87c00f4a9bb4b95c6b0d0b6c65 100644
--- a/src/correlate/Utils.h
+++ b/src/correlate/Utils.h
@@ -4,11 +4,9 @@
 #include <vector>
 #include <map>
 
-#include <QVector.h>
-#include <Correlation.h>
-#include <ReSampler.h>
+#include "QnTools/ReSampleHelper.hpp"
 
-#include "src/base/QvectorConfig.h"
+#include "base/QvectorConfig.h"
 
 #include "Functions.h"
 
@@ -18,26 +16,26 @@ using function2_t = const std::function<double(const Qn::QVector &a, const Qn::Q
 using function3_t = const std::function<double(const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c)>;
 using function4_t = const std::function<double(const Qn::QVector &a, const Qn::QVector &b, const Qn::QVector &c, const Qn::QVector &d)>;
 
-[[maybe_unused]] std::vector<std::pair<std::string, function2_t>> Q1Q1(bool non_zero_only) {
+[[maybe_unused]] static std::vector<std::pair<std::string, function2_t>> Q1Q1(bool non_zero_only) {
   if(non_zero_only){
-    return {{"XX", XX}, {"YY", YY}};
+    return {{"Q1x_Q1x", XX}, {"Q1y_Q1y", YY}};
   }
   else{
-    return {{"XX", XX}, {"YY", YY}, {"XY", XY}, {"YX", YX}};
+    return {{"Q1x_Q1x", XX}, {"Q1y_Q1y", YY}, {"Q1x_Q1y", XY}, {"Q1y_Q1x", YX}};
   }
 }
 
-[[maybe_unused]] std::vector<std::pair<std::string, function3_t>> Q2Q1Q1(bool non_zero_only) {
+[[maybe_unused]] static std::vector<std::pair<std::string, function3_t>> Q2Q1Q1(bool non_zero_only) {
   if(non_zero_only){
-    return {{"Y2XY", Y2XY}, {"Y2YX", Y2YX}, {"X2XX", X2XX}, {"X2YY", X2YY}};
+    return {{"Q2y_Q1x_Q1y", Y2XY}, {"Q2y_Q1y_Q1x", Y2YX}, {"Q2x_Q1x_Q1x", X2XX}, {"Q2x_Q1y_Q1y", X2YY}};
   }
   else{
-    return{ {"Y2XY", Y2XY}, {"Y2YX", Y2YX}, {"X2XX", X2XX}, {"X2YY", X2YY},
-            {"X2XY", X2XY}, {"X2YX", X2YX}, {"Y2XX", Y2XX}, {"Y2YY", Y2YY}};
+    return{ {"Q2y_Q1x_Q1y", Y2XY}, {"Q2y_Q1y_Q1x", Y2YX}, {"Q2x_Q1x_Q1x", X2XX}, {"Q2x_Q1y_Q1y", X2YY},
+            {"Q2x_Q1x_Q1y", X2XY}, {"Q2x_Q1y_Q1x", X2YX}, {"Q2y_Q1x_Q1x", Y2XX}, {"Q2y_Q1y_Q1y", Y2YY}};
   }
 }
 
-[[maybe_unused]] std::vector<std::pair<std::string, function3_t>> Q3Q1Q1Q1 = {};  // TODO
+[[maybe_unused]] static std::vector<std::pair<std::string, function3_t>> Q3Q1Q1Q1 = {};  // TODO
 
 }
 
diff --git a/src/tasks/Analyze.cpp b/src/tasks/Analyze.cpp
index f87f676852112aca3738d335163d05208e7a0f04..8cd141f31842cb4f438e225708b6b4ef7cfa6045 100644
--- a/src/tasks/Analyze.cpp
+++ b/src/tasks/Analyze.cpp
@@ -1,10 +1,10 @@
 #include <string>
 #include <vector>
 
-#include "src/analyze/ResourceManager.h"
-#include "src/analyze/Resolution.h"
-#include "src/analyze/FlowCoefficients.h"
-#include "src/analyze/Utils.h"
+#include "analyze/ResourceManager.h"
+#include "analyze/Resolution.h"
+#include "analyze/FlowCoefficients.h"
+#include "analyze/Utils.h"
 
 void DefineResolution3S(ResourceManager& rm);
 void DefineResolution4S(ResourceManager& rm);
diff --git a/src/tasks/CMakeLists.txt b/src/tasks/CMakeLists.txt
index 6e1d25d00d01efd5a4a842fb160f9e9bbcc73289..05657afed527720706801c9d4abd557b074a2808 100644
--- a/src/tasks/CMakeLists.txt
+++ b/src/tasks/CMakeLists.txt
@@ -1,21 +1,21 @@
-include_directories(${PROJECT_INCLUDE_DIRECTORIES} ${CMAKE_CURRENT_SOURCE_DIR})
+include_directories(${PROJECT_INCLUDE_DIRECTORIES} ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_SOURCE_DIR}/src)
 link_directories(${PROJECT_LINK_DIRECTORIES})
 
 add_executable(correct Correct.cpp)
 add_dependencies(correct ${PROJECT_DEPENDENCIES})
-target_link_libraries(correct FlowCorrect FlowBase Correction Base AnalysisTreeBase AnalysisTreeInfra ${ROOT_LIBRARIES})
+target_link_libraries(correct FlowCorrect FlowBase  QnToolsCorrection QnToolsBase AnalysisTreeBase AnalysisTreeInfra ${ROOT_LIBRARIES})
 
-add_executable(correlate Correlate.cpp ../correlate/CorrelationTask.cpp ../correlate/Utils.h ../correlate/Functions.h)
+add_executable(correlate Correlate.cpp)
 add_dependencies(correlate ${PROJECT_DEPENDENCIES})
-target_link_libraries(correlate FlowBase Correction Base Correlation ${ROOT_LIBRARIES} ROOTVecOps)
+target_link_libraries(correlate FlowBase FlowCorrelate QnToolsBase)
 
 #add_executable(analyze Analyze.cpp)
 #add_dependencies(analyze FlowAna ${PROJECT_DEPENDENCIES})
-#target_link_libraries(analyze Base Correlation FlowAna ${ROOT_LIBRARIES})
+#target_link_libraries(analyze QnToolsBase FlowAna ${ROOT_LIBRARIES})
 #
 #add_executable(profiles Profiles.cpp)
 #add_dependencies(profiles FlowAna ${PROJECT_DEPENDENCIES})
-#target_link_libraries(profiles Base Correlation FlowAna ${ROOT_LIBRARIES})
+#target_link_libraries(profiles QnToolsBase FlowAna ${ROOT_LIBRARIES})
 
 install(TARGETS correct correlate DESTINATION bin )
 #install(TARGETS correct correlate analyze profiles DESTINATION bin )
diff --git a/src/tasks/Correct.cpp b/src/tasks/Correct.cpp
index dce40040c60892478be6fe3b37e38c1698782017..52c7830cb1b4c6faa237ee1b578fff3fb7af1c26 100644
--- a/src/tasks/Correct.cpp
+++ b/src/tasks/Correct.cpp
@@ -2,15 +2,16 @@
 #include <chrono>
 
 #include <TSystem.h>
-#include "src/correct/CorrectionTask.h"
-#include "src/base/GlobalConfig.h"
-#include "Axis.h"
-
-#include "src/base/QvectorChannelsConfig.h"
-#include "src/base/QvectorTracksConfig.h"
 
+#include "QnTools/Axis.hpp"
 #include "CbmCuts.h"
 
+#include "base/QvectorChannelsConfig.h"
+#include "base/QvectorTracksConfig.h"
+#include "base/GlobalConfig.h"
+#include "correct/CorrectionTask.h"
+
+
 int main(int argc, char **argv) {
   using namespace std;
 //  ROOT::EnableImplicitMT(8);
@@ -44,6 +45,7 @@ int main(int argc, char **argv) {
 
 //  gConfig->SetEventCuts(AnalysisTree::GetCbmEventCuts(rec_event));
   gConfig->AddEventVar({"Centrality", ana_event, "tracks_centrality"});
+  gConfig->AddEventVar({"Psi_RP", sim_event, "psi_RP"});
   gConfig->AddCorrectionAxis( {"Centrality", {0., 5, 10, 20, 30, 50, 70, 100.}} );
 
   gConfig->AddChannelQvector({psd_modules, psd_config, 1});
@@ -83,8 +85,8 @@ int main(int argc, char **argv) {
  // ***********************************************
   // first filelist should contain DataHeader
   Qn::CorrectionTask task( {fileList, friend_filelist}, {"aTree", "cTree"}, calibFileName);
-  task.AddBranch({rec_particles});
-  task.AddBranch({sim_particles, AnalysisTree::GetCbmMcTracksCuts(sim_particles)});
+  task.AddBranch(BranchReader(rec_particles));
+  task.AddBranch(BranchReader(sim_particles, AnalysisTree::GetCbmMcTracksCuts(sim_particles)) );
 
   task.SetRecEventName(rec_event);
   task.SetSimEventName(sim_event);
diff --git a/src/tasks/Correlate.cpp b/src/tasks/Correlate.cpp
index b17cee5d22c3ca8a58f8089294611a9b9c738076..6d3379127d886b6b5344b81001ea920ec253fac8 100644
--- a/src/tasks/Correlate.cpp
+++ b/src/tasks/Correlate.cpp
@@ -2,7 +2,7 @@
 #include <chrono>
 #include <string>
 
-#include "src/correlate/CorrelationTask.h"
+#include "correlate/CorrelationTask.h"
 
 int main(int argc, char **argv) {
   using namespace std;
@@ -16,14 +16,7 @@ int main(int argc, char **argv) {
   const std::string& file{argv[1]};
   CorrelationTask st(file, "tree");
 
-//  st.AddEventAxis( {"Centrality", {0., 5, 10, 20, 30, 50, 70, 100.}} );
-//
-//  st.SetDetectorsChannels({"psd1", "psd2", "psd3"});
-//  st.SetDetectorsTracks({"pion_neg", "pion", "proton", "mc_pion_neg", "mc_pion", "mc_proton"});
-////  st.SetDetectorsMc();
-//
-//  st.SetIsSimulation(true);
-//  st.SetNSamples(50);
+//  st.GetConfig()->SetNSamples();
 
   std::cout << "go" << std::endl;
   auto start = std::chrono::system_clock::now();
diff --git a/src/tasks/Profiles.cpp b/src/tasks/Profiles.cpp
index b13b52d10c952219fa0f31db6355112f6e32dd3c..3a0864ca50583f3a6cda6e80c0d75bfca046fa02 100644
--- a/src/tasks/Profiles.cpp
+++ b/src/tasks/Profiles.cpp
@@ -1,5 +1,5 @@
 #include <string>
-#include <src/analyze/ResourceManager.h>
+#include <analyze/ResourceManager.h>
 int main(int argc, char* argv[]) {
 
   std::string inputfile = "/home/vklochkov/Data/cbm/oct19_fr_18.2.1_fs_jun19p1/dcmqgsm_smm_pluto/auau/12agev/mbias/psd44_hole20_pipe0/TGeant4/flow_psdmc/analysis.root";
diff --git a/src/tools/CMakeLists.txt b/src/tools/CMakeLists.txt
index 197ab00a314b374dd8448f481efe06c7b783a871..40f989c5c9f85c5a5643c00894db3f8a8eef780b 100644
--- a/src/tools/CMakeLists.txt
+++ b/src/tools/CMakeLists.txt
@@ -1,6 +1,8 @@
 include_directories(${PROJECT_INCLUDE_DIRECTORIES} ${CMAKE_CURRENT_SOURCE_DIR})
 link_directories(${PROJECT_LINK_DIRECTORIES})
 
-add_executable(QVectorQA QVectorQA.cxx)
-add_dependencies(QVectorQA ${PROJECT_DEPENDENCIES})
-target_link_libraries(QVectorQA QnCorrections Base ${ROOT_LIBRARIES} )
+add_executable(qvector_qa QVectorQA.cxx)
+add_dependencies(qvector_qa ${PROJECT_DEPENDENCIES})
+target_link_libraries(qvector_qa QnToolsBase ${ROOT_LIBRARIES} )
+
+install(TARGETS qvector_qa DESTINATION bin )
diff --git a/src/tools/QVectorQA.cxx b/src/tools/QVectorQA.cxx
index 66c76dd722dc6d0076f21a99923362826250901b..14fe32fb91d225da7c860e827b8f04be3fda3b0e 100644
--- a/src/tools/QVectorQA.cxx
+++ b/src/tools/QVectorQA.cxx
@@ -4,9 +4,11 @@
 
 #include <iostream>
 
-#include <DataContainer.h>
-#include <QVector.h>
+#include <QnTools/DataContainer.hpp>
+#include <QnTools/QVector.hpp>
+
 #include <TTree.h>
+#include <TFile.h>
 #include <TTreeReaderValue.h>
 #include <TTreeReader.h>
 #include <TDirectory.h>
@@ -38,7 +40,7 @@ T *createOrGet(TDirectory *tgt, const std::string &namecycle, T &&proto) {
     return objPtr;
   }
 
-  auto pwd = gDirectory;
+  auto *pwd = gDirectory;
   tgt->cd();
   objPtr = new T(proto);
   objPtr->SetName(namecycle.c_str());
@@ -87,10 +89,12 @@ int main(int argc, char **argv) {
 
   TObjArray *branchList = tree->GetListOfBranches();
 
-  for (auto objPtr : *branchList) {
-    auto branchPtr = dynamic_cast<TBranch *>(objPtr);
+  for (auto *objPtr : *branchList) {
+    auto *branchPtr = dynamic_cast<TBranch *>(objPtr);
+
+    std::cout << std::string(branchPtr->GetClassName()) << std::endl;
 
-    if (std::string(branchPtr->GetClassName()) == "Qn::DataContainer<Qn::QVector>") {
+    if (std::string(branchPtr->GetClassName()) == "Qn::DataContainer<Qn::QVector,Qn::Axis<double> >") {
       qVecQAEntries.emplace_back(
           std::string(branchPtr->GetName()), QVecQAEntry::Reader_t(treeReader, branchPtr->GetName()));
 
@@ -100,29 +104,30 @@ int main(int argc, char **argv) {
   Long64_t iEvent{0};
   for (auto &qVecQAEntry : qVecQAEntries) {
     cout << qVecQAEntry.name << endl;
-    auto qVecQAEntryTopDir = mkdirOrGet(qaFile.get(), qVecQAEntry.name);
+    auto *qVecQAEntryTopDir = mkdirOrGet(qaFile.get(), qVecQAEntry.name);
 
     iEvent = 0;
     treeReader.Restart();
     while (treeReader.Next()) {
 
-      auto &qVecContainer = qVecQAEntry.getVal();
+      std::cout << "Event # " << iEvent+1 << "\r" << std::flush;
+
+      const auto &qVecContainer = qVecQAEntry.getVal();
 
       unsigned int iContainerBin = 0;
       for (const auto &qVecPtr : qVecContainer) {
-        auto nHarmonics = qVecPtr.q_.size();
-
-        auto binDir = mkdirOrGet(qVecQAEntryTopDir, Form("bin_%d", iContainerBin));
+        auto nHarmonics = qVecPtr.GetNoOfHarmonics();
+        auto *binDir = mkdirOrGet(qVecQAEntryTopDir, Form("bin_%d", iContainerBin));
         for (decltype(nHarmonics) iHarmonic = 1; iHarmonic < nHarmonics; ++iHarmonic) {
           double qX = qVecPtr.x(iHarmonic);
           double qY = qVecPtr.y(iHarmonic);
-          double weight = qVecPtr.sum_weights_;
+          double weight = qVecPtr.sumweights();
 
-          auto h2QyVsQx =
-              createOrGet(binDir, Form("QyVsQx_H%lu", iHarmonic), TH2D("test", "", 400, -2, 2, 400, -2, 2));
+          auto *h2QyVsQx =
+              createOrGet(binDir, Form("QyVsQx_H%u", iHarmonic), TH2D("test", "", 400, -2, 2, 400, -2, 2));
           h2QyVsQx->Fill(qX, qY);
-          auto h2QyVsQx_weighted =
-              createOrGet(binDir, Form("QyVsQxWeighted_H%lu", iHarmonic), TH2D("test", "", 400, -2, 2, 400, -2, 2));
+          auto *h2QyVsQx_weighted =
+              createOrGet(binDir, Form("QyVsQxWeighted_H%u", iHarmonic), TH2D("test", "", 400, -2, 2, 400, -2, 2));
           h2QyVsQx_weighted->Fill(qX, qY, weight);
         }