Edwin update mille

.clangd

@@ -0,0 +1,4 @@
+Diagnostics:
+  UnusedIncludes: Strict
+  MissingIncludes: Strict
+  AnalyzeAngledIncludes: Yes

neuland/calibration/cal_to_hit/CMakeLists.txt

@@ -19,6 +19,7 @@ set(LINKDEF NeulandCalToHitLinkDef.h)
 set(SRCS
     # cmake-format: sortable
     engine/R3BNeulandLSQREngineAdaptor.cxx
+    engine/R3BNeulandMilleCalDataProcessor.cxx
     engine/R3BNeulandMillepede.cxx
     engine/R3BNeulandPredecessor.cxx
     R3BNeulandCalToHitPar.cxx
@@ -29,6 +30,7 @@ set(HEADERS
     # cmake-format: sortable
     engine/R3BNeulandCosmicEngine.h
     engine/R3BNeulandLSQREngineAdaptor.h
+    engine/R3BNeulandMilleCalDataProcessor.h
     engine/R3BNeulandMillepede.h
     engine/R3BNeulandPredecessor.h
     R3BNeulandCalToHitPar.h

neuland/calibration/cal_to_hit/R3BNeulandCalToHitPar.h

@@ -77,6 +77,10 @@ namespace R3B::Neuland
         {
             return module_pars_.at(module_num);
         }
+        auto HasModuleParAt(int module_num) const -> bool
+        {
+            return module_pars_.find(module_num) != module_pars_.end();
+        }
         auto GetModulePars() const -> const std::unordered_map<unsigned int, ::R3B::Neuland::HitModulePar>&
         {
             return module_pars_;

neuland/calibration/cal_to_hit/R3BNeulandCalToHitParTask.cxx

@@ -12,10 +12,20 @@
  ******************************************************************************/
 
 #include "R3BNeulandCalToHitParTask.h"
+#include "R3BDataMonitor.h"
+#include "R3BException.h"
+#include "R3BNeulandBasePar.h"
+#include "R3BNeulandCalToHitPar.h"
+#include "R3BNeulandCalibrationTask.h"
+#include "R3BNeulandCommon.h"
+#include <FairRootManager.h>
+#include <FairRuntimeDb.h>
 #include <R3BLogger.h>
 #include <R3BNeulandLSQREngineAdaptor.h>
 #include <R3BNeulandMillepede.h>
 #include <R3BNeulandPredecessor.h>
+#include <memory>
+#include <string_view>
 
 namespace R3B::Neuland
 {
@@ -29,22 +39,24 @@ namespace R3B::Neuland
         : CalibrationTask(name, iVerbose)
         , cal_data_{ cal_data_name }
         , base_par_{ InputPar<CalibrationBasePar>(base_par_name) }
-        , hit_par_{ OutputPar<Cal2HitPar>(hit_par_name) }
     // NOLINTEND
     {
         switch (method)
         {
             case Cal2HitParMethod::LSQT:
                 R3BLOG(info, "Cal2HitPar method: LSQT.");
                 engine_ = std::make_unique<Calibration::LSQREngineAdaptor>();
+                hit_par_ = OutputPar<Cal2HitPar>(hit_par_name);
                 break;
-            case Cal2HitParMethod::millipede:
+            case Cal2HitParMethod::millepede:
                 R3BLOG(info, "Cal2HitPar method: Millepede.");
                 engine_ = std::make_unique<Calibration::MillepedeEngine>();
+                hit_par_ = InputPar<Cal2HitPar>(hit_par_name);
                 break;
             case Cal2HitParMethod::predecessor:
                 R3BLOG(info, "Cal2HitPar method: predecessor.");
                 engine_ = std::make_unique<Calibration::Predecessor>();
+                hit_par_ = OutputPar<Cal2HitPar>(hit_par_name);
                 break;
         }
     }
@@ -70,10 +82,11 @@ namespace R3B::Neuland
     void Cal2HitParTask::TriggeredExec()
     {
         engine_->EventReset();
-        for (const auto& bar_signal : cal_data_)
-        {
-            engine_->AddSignal(bar_signal);
-        }
+        engine_->AddSignals(cal_data_.get());
+        // for (const auto& bar_signal : cal_data_)
+        // {
+        //     engine_->AddSignal(bar_signal);
+        // }
         auto* eventHeader = GetEventHeader();
         if (eventHeader != nullptr)
         {

neuland/calibration/cal_to_hit/R3BNeulandCalToHitParTask.h

@@ -28,7 +28,7 @@ namespace R3B::Neuland
     {
         LSQT,
         predecessor,
-        millipede
+        millepede
     };
 
     class Cal2HitParTask : public CalibrationTask

neuland/calibration/cal_to_hit/engine/R3BNeulandCosmicEngine.h

@@ -43,7 +43,7 @@ namespace R3B::Neuland::Calibration
 
         virtual void Init() {}
         virtual auto SignalFilter(const std::vector<BarCalData>& /*signals*/) -> bool { return true; }
-        virtual void AddSignal(const BarCalData& signal) = 0;
+        virtual void AddSignals(const std::vector<BarCalData>& signal) = 0;
         virtual void Calibrate(Cal2HitPar& hit_par) = 0;
         virtual void SetMinStat(int min) {}
         virtual void BeginOfEvent(unsigned int event_num) {};

neuland/calibration/cal_to_hit/engine/R3BNeulandLSQREngineAdaptor.cxx

@@ -12,13 +12,19 @@
  ******************************************************************************/
 
 #include "R3BNeulandLSQREngineAdaptor.h"
-#include "R3BDataMonitor.h"
-#include "R3BNeulandHitPar.h"
+
 #include <FairRuntimeDb.h>
+#include <R3BException.h>
 #include <R3BLogger.h>
+#include <R3BNeulandCalData2.h>
+#include <R3BNeulandCalToHitPar.h>
+#include <R3BNeulandCommon.h>
+#include <R3BNeulandHitModulePar.h>
+#include <R3BNeulandHitPar.h>
 #include <R3BNeulandParDirCreator.h>
 #include <R3BShared.h>
-#include <algorithm>
+#include <cmath>
+#include <vector>
 
 namespace
 {
@@ -74,10 +80,13 @@ namespace R3B::Neuland::Calibration
         hit_cal_engine_.Init(&hit_par_temp);
     }
 
-    void LSQREngineAdaptor::AddSignal(const BarCalData& signal)
+    void LSQREngineAdaptor::AddSignals(const std::vector<BarCalData>& signals)
     {
-        add_bar_signal(signal, Side::left);
-        add_bar_signal(signal, Side::right);
+        for (const auto& signal : signals)
+        {
+            add_bar_signal(signal, Side::left);
+            add_bar_signal(signal, Side::right);
+        }
     }
 
     void LSQREngineAdaptor::add_bar_signal(const BarCalData& barSignal, Side side)

neuland/calibration/cal_to_hit/engine/R3BNeulandLSQREngineAdaptor.h

@@ -14,9 +14,13 @@
 #pragma once
 
 #include "R3BNeulandCosmicEngine.h"
+
+#include <R3BNeulandCalData2.h>
+#include <R3BNeulandCalToHitPar.h>
 #include <R3BNeulandCosmicTracker.h>
 #include <R3BNeulandHitCalibrationEngine.h>
 #include <R3BShared.h>
+#include <vector>
 
 namespace R3B::Neuland::Calibration
 {
@@ -31,7 +35,7 @@ namespace R3B::Neuland::Calibration
         CosmicTracker cosmic_tracker_;
 
         void Init() override;
-        void AddSignal(const BarCalData& signal) override;
+        void AddSignals(const std::vector<BarCalData>& signals) override;
         void Calibrate(Cal2HitPar& hit_par) override;
         void EndOfEvent(unsigned int event_num = 0) override;
         void EventReset() override

neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx

@@ -0,0 +1,187 @@
+#include "R3BNeulandMilleCalDataProcessor.h"
+#include "R3BLogger.h"
+#include "R3BNeulandCalData2.h"
+#include <Fit/BinData.h>
+#include <Math/WrappedMultiTF1.h>
+#include <R3BNeulandCommon.h>
+#include <TError.h>
+#include <algorithm>
+#include <fmt/core.h>
+#include <numeric>
+#include <range/v3/algorithm/find_if.hpp>
+#include <vector>
+
+// NOLINTNEXTLINE
+#include <fmt/ranges.h>
+
+namespace R3B::Neuland::Calibration
+{
+    MilleDataProcessor::MilleDataProcessor(int num_of_modules)
+    {
+        init_data_registers(num_of_modules);
+        fitter_.SetFunction(
+            ROOT::Math::WrappedMultiTF1{ fit_function_, static_cast<unsigned int>(fit_function_.GetNdim()) }, true);
+        fitter_.Config().SetMinimizer("Linear");
+    }
+
+    void MilleDataProcessor::init_data_registers(int num_of_modules)
+    {
+        const auto num_of_planes = num_of_modules / BarsPerPlane;
+        for (int plane_id{}; plane_id < num_of_planes; ++plane_id)
+        {
+            auto data_iter = data_regsiters_.insert_or_assign(plane_id, std::vector<MilleCalData>{}).first;
+            data_iter->second.reserve(BarsPerPlane);
+        }
+    }
+
+    void MilleDataProcessor::reset_fitpars()
+    {
+        fitter_.Config().ParSettings(0).SetValue(0.);
+        fitter_.Config().ParSettings(1).SetValue(0.);
+    }
+
+    void MilleDataProcessor::reset()
+    {
+        for (auto& [plane_id, bar_data] : data_regsiters_)
+        {
+            bar_data.clear();
+        }
+        fit_result_.x_z = FitPar{};
+        fit_result_.y_z = FitPar{};
+        x_z_vals_.clear();
+        y_z_vals_.clear();
+        reset_fitpars();
+    }
+
+    auto MilleDataProcessor::filter(const std::vector<BarCalData>& signals) -> bool
+    {
+        // fill only the bar_cal_data with only one pmt signal on both sides
+        for (const auto& signal : signals)
+        {
+            if (signal.left.size() == 1 && signal.right.size() == 1)
+            {
+                const auto bar_num = static_cast<int>(signal.module_num);
+                const auto plane_id = ModuleID2PlaneID(bar_num - 1);
+                data_regsiters_[plane_id].emplace_back(signal);
+            }
+        }
+
+        remove_isolated_bar_signal();
+
+        return fit_planes();
+    }
+
+    void MilleDataProcessor::remove_isolated_bar_signal()
+    {
+        namespace rng = ranges;
+        for (auto& [plane_id, bar_data] : data_regsiters_)
+        {
+            if (bar_data.size() < 2)
+            {
+                continue;
+            }
+
+            const auto& bar_signals = bar_data;
+            auto check_if_isolated = [&bar_signals](const auto& signal) -> bool
+            {
+                const auto module_num = static_cast<int>(signal.module_num);
+                return (rng::find_if(bar_signals,
+                                     [module_num](const auto& bar_signal)
+                                     { return bar_signal.module_num == module_num - 1; }) == bar_signals.end()) and
+                       (rng::find_if(bar_signals,
+                                     [module_num](const auto& bar_signal)
+                                     { return bar_signal.module_num == module_num + 1; }) == bar_signals.end());
+            };
+
+            bar_data.erase(std::remove_if(bar_data.begin(), bar_data.end(), check_if_isolated), bar_data.end());
+        }
+    }
+
+    auto MilleDataProcessor::fit_planes() -> bool
+    {
+        fill_fit_data();
+
+        return fit_plane_data();
+    }
+
+    void MilleDataProcessor::fill_fit_data()
+    {
+        for (auto& [plane_id, bar_data] : data_regsiters_)
+        {
+            if (bar_data.empty())
+            {
+                continue;
+            }
+            const auto is_plane_horizontal = IsPlaneIDHorizontal(plane_id);
+            auto& fit_data = is_plane_horizontal ? y_z_vals_ : x_z_vals_;
+            const auto z_val = PlaneID2ZPos(plane_id);
+
+            const auto displacement =
+                std::accumulate(bar_data.begin(),
+                                bar_data.end(),
+                                0.,
+                                [](double sum, const MilleCalData& signal)
+                                { return sum + GetBarVerticalDisplacement(static_cast<int>(signal.module_num)); }) /
+                static_cast<double>(bar_data.size());
+            fit_data.z_vals.push_back(z_val);
+            fit_data.z_errs.push_back(BarSize_Z / 2.);
+            fit_data.errs.push_back(BarSize_XY / 2.);
+            fit_data.vals.push_back(displacement);
+        }
+    }
+
+    auto MilleDataProcessor::fit_plane_data() -> bool
+    {
+        return linear_fit(x_z_vals_, fit_result_.x_z) and linear_fit(y_z_vals_, fit_result_.y_z);
+    }
+
+    auto MilleDataProcessor::calculate_residual(double val, int module_num) const -> double
+    {
+        const auto z_val = ModuleNum2ZPos(module_num);
+        const auto is_plane_horizontal = IsPlaneIDHorizontal(ModuleID2PlaneID(module_num - 1));
+        const auto& fit_result = is_plane_horizontal ? fit_result_.x_z : fit_result_.y_z;
+        const auto diff = val - (fit_result.slope * z_val) - fit_result.offset;
+        return diff * diff;
+    }
+
+    auto MilleDataProcessor::linear_fit(const FitData& data, FitPar& fit_par) -> bool
+    {
+
+        const auto bin_data = ROOT::Fit::BinData{ static_cast<unsigned int>(data.size()),
+                                                  data.z_vals.data(),
+                                                  data.vals.data(),
+                                                  data.z_errs.data(),
+                                                  data.errs.data() };
+        reset_fitpars();
+
+        // disable annoying root printouts
+        auto old_var = gErrorIgnoreLevel;
+        gErrorIgnoreLevel = kFatal;
+        auto res = fitter_.Fit(bin_data);
+        gErrorIgnoreLevel = old_var;
+
+        if (not res)
+        {
+            R3BLOG(debug, "Linear fitting on x_z data failed");
+            return false;
+        }
+        fit_par.slope = fitter_.Result().Parameter(0);
+        fit_par.offset = fitter_.Result().Parameter(1);
+        if (fitter_.Result().Prob() < p_value_cut_)
+        {
+            R3BLOG(debug,
+                   fmt::format("p-value ({}) is too small from the fit.\n\
+                       x = np.array({}) \n\
+                       x_err = np.array({}) \n\
+                       y = np.array({}) \n\
+                       y_err = np.array({})",
+                               fitter_.Result().Prob(),
+                               data.z_vals,
+                               data.z_errs,
+                               data.vals,
+                               data.errs));
+            return false;
+        }
+        return true;
+    }
+} // namespace R3B::Neuland::Calibration