From 2ad70269a07e90b3c1ce52ef481e26028b00acce Mon Sep 17 00:00:00 2001 From: YanzhaoW Date: Sat, 15 Feb 2025 23:24:59 +0100 Subject: [PATCH 1/6] add data processor for mille --- neuland/calibration/cal_to_hit/CMakeLists.txt | 2 + .../cal_to_hit/R3BNeulandCalToHitParTask.cxx | 9 ++- .../engine/R3BNeulandCosmicEngine.h | 2 +- .../engine/R3BNeulandLSQREngineAdaptor.cxx | 9 ++- .../engine/R3BNeulandLSQREngineAdaptor.h | 2 +- .../R3BNeulandMilleCalDataProcessor.cxx | 75 +++++++++++++++++++ .../engine/R3BNeulandMilleCalDataProcessor.h | 50 +++++++++++++ .../cal_to_hit/engine/R3BNeulandMillepede.cxx | 19 +++-- .../cal_to_hit/engine/R3BNeulandMillepede.h | 3 +- .../engine/R3BNeulandPredecessor.cxx | 11 ++- .../cal_to_hit/engine/R3BNeulandPredecessor.h | 2 +- 11 files changed, 161 insertions(+), 23 deletions(-) create mode 100644 neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx create mode 100644 neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.h diff --git a/neuland/calibration/cal_to_hit/CMakeLists.txt b/neuland/calibration/cal_to_hit/CMakeLists.txt index 4b7222190..82b7bd1a2 100644 --- a/neuland/calibration/cal_to_hit/CMakeLists.txt +++ b/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 diff --git a/neuland/calibration/cal_to_hit/R3BNeulandCalToHitParTask.cxx b/neuland/calibration/cal_to_hit/R3BNeulandCalToHitParTask.cxx index 0b634e586..081cf0634 100644 --- a/neuland/calibration/cal_to_hit/R3BNeulandCalToHitParTask.cxx +++ b/neuland/calibration/cal_to_hit/R3BNeulandCalToHitParTask.cxx @@ -70,10 +70,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) { diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandCosmicEngine.h b/neuland/calibration/cal_to_hit/engine/R3BNeulandCosmicEngine.h index c9aec6ab5..05f27ff8d 100644 --- a/neuland/calibration/cal_to_hit/engine/R3BNeulandCosmicEngine.h +++ b/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& /*signals*/) -> bool { return true; } - virtual void AddSignal(const BarCalData& signal) = 0; + virtual void AddSignals(const std::vector& signal) = 0; virtual void Calibrate(Cal2HitPar& hit_par) = 0; virtual void SetMinStat(int min) {} virtual void BeginOfEvent(unsigned int event_num) {}; diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandLSQREngineAdaptor.cxx b/neuland/calibration/cal_to_hit/engine/R3BNeulandLSQREngineAdaptor.cxx index 0b3306587..6b87541ee 100644 --- a/neuland/calibration/cal_to_hit/engine/R3BNeulandLSQREngineAdaptor.cxx +++ b/neuland/calibration/cal_to_hit/engine/R3BNeulandLSQREngineAdaptor.cxx @@ -74,10 +74,13 @@ namespace R3B::Neuland::Calibration hit_cal_engine_.Init(&hit_par_temp); } - void LSQREngineAdaptor::AddSignal(const BarCalData& signal) + void LSQREngineAdaptor::AddSignals(const std::vector& 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) diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandLSQREngineAdaptor.h b/neuland/calibration/cal_to_hit/engine/R3BNeulandLSQREngineAdaptor.h index 71f24e0ca..aa3580036 100644 --- a/neuland/calibration/cal_to_hit/engine/R3BNeulandLSQREngineAdaptor.h +++ b/neuland/calibration/cal_to_hit/engine/R3BNeulandLSQREngineAdaptor.h @@ -31,7 +31,7 @@ namespace R3B::Neuland::Calibration CosmicTracker cosmic_tracker_; void Init() override; - void AddSignal(const BarCalData& signal) override; + void AddSignals(const std::vector& signals) override; void Calibrate(Cal2HitPar& hit_par) override; void EndOfEvent(unsigned int event_num = 0) override; void EventReset() override diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx b/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx new file mode 100644 index 000000000..267012486 --- /dev/null +++ b/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx @@ -0,0 +1,75 @@ +#include "R3BNeulandMilleCalDataProcessor.h" +#include +#include +#include + +namespace R3B::Neuland::Calibration +{ + MilleDataProcessor::MilleDataProcessor(int num_of_modules) { init_data_registers(num_of_modules); } + + 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{}).first; + data_iter->second.reserve(BarsPerPlane); + } + } + + void MilleDataProcessor::reset() + { + for (auto& [plane_id, bar_data] : data_regsiters_) + { + bar_data.clear(); + } + fit_result_.x_z = FitPar{}; + fit_result_.y_z = FitPar{}; + } + + auto MilleDataProcessor::operator()(const std::vector& signals) -> const auto& + { + // 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(signal.module_num); + const auto plane_id = ModuleID2PlaneID(bar_num - 1); + data_regsiters_[plane_id].emplace_back(signal); + } + } + + remove_isolated_bar_signal(); + + return *this; + } + + 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(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()); + } + } + + void MilleDataProcessor::fit_planes() {} +} // namespace R3B::Neuland::Calibration diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.h b/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.h new file mode 100644 index 000000000..8ac2e3990 --- /dev/null +++ b/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.h @@ -0,0 +1,50 @@ +#pragma once + +#include + +namespace R3B::Neuland::Calibration +{ + struct MilleCalData + { + public: + MilleCalData() = default; + explicit MilleCalData(const BarCalData& bar_cal_data) + : module_num{ bar_cal_data.module_num } + , left{ bar_cal_data.left.front() } + , right{ bar_cal_data.right.front() } + { + } + unsigned int module_num = 0; // 1 based bar num + CalDataSignal left; + CalDataSignal right; + }; + + class MilleDataProcessor + { + public: + struct FitPar + { + double offset = 0.; + double slope = 0.; + }; + struct FitResult + { + FitPar x_z; + FitPar y_z; + }; + + explicit MilleDataProcessor(int num_of_modules); + auto operator()(const std::vector& signals) -> const auto&; + [[nodiscard]] auto get_data() const -> const auto& { return data_regsiters_; } + [[nodiscard]] auto get_fit() const -> const auto& { return fit_result_; } + void reset(); + + private: + std::unordered_map> data_regsiters_; + FitResult fit_result_; + + void init_data_registers(int num_of_modules); + void remove_isolated_bar_signal(); + void fit_planes(); + }; +} // namespace R3B::Neuland::Calibration diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandMillepede.cxx b/neuland/calibration/cal_to_hit/engine/R3BNeulandMillepede.cxx index 72d3c7eda..37989653f 100644 --- a/neuland/calibration/cal_to_hit/engine/R3BNeulandMillepede.cxx +++ b/neuland/calibration/cal_to_hit/engine/R3BNeulandMillepede.cxx @@ -290,17 +290,20 @@ namespace R3B::Neuland::Calibration "Writting Mille data to binary file with meas = {} and z = {}", input_data_buffer_.measurement, pos_z)); } - void MillepedeEngine::AddSignal(const BarCalData& signal) + void MillepedeEngine::AddSignals(const std::vector& signals) { - // all bar signal must have one signal on both sides - if (signal.left.size() != 1 or signal.right.size() != 1) + for (const auto& signal : signals) { - return; - } + // all bar signal must have one signal on both sides + if (signal.left.size() != 1 or signal.right.size() != 1) + { + return; + } - add_signal_t_sum(signal); - add_signal_t_diff(signal); - add_spacial_local_constraint(static_cast(signal.module_num)); + add_signal_t_sum(signal); + add_signal_t_diff(signal); + add_spacial_local_constraint(static_cast(signal.module_num)); + } } void MillepedeEngine::Calibrate(Cal2HitPar& hit_par) diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandMillepede.h b/neuland/calibration/cal_to_hit/engine/R3BNeulandMillepede.h index 72120f57a..d029cf515 100644 --- a/neuland/calibration/cal_to_hit/engine/R3BNeulandMillepede.h +++ b/neuland/calibration/cal_to_hit/engine/R3BNeulandMillepede.h @@ -49,6 +49,7 @@ namespace R3B::Neuland::Calibration std::string input_data_filename_ = "neuland_cosmic_mille.bin"; std::string pede_steer_filename_ = "neuland_steer.txt"; std::string parameter_filename_ = "neuland_pars.txt"; + Mille binary_data_writer_{ input_data_filename_ }; Millepede::ResultReader par_result_; Millepede::Launcher pede_launcher_; @@ -62,7 +63,7 @@ namespace R3B::Neuland::Calibration Cal2HitPar* cal_to_hit_par_ = nullptr; void Init() override; - void AddSignal(const BarCalData& signal) override; + void AddSignals(const std::vector& signals) override; void Calibrate(Cal2HitPar& hit_par) override; void EndOfEvent(unsigned int event_num = 0) override; void EventReset() override; diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandPredecessor.cxx b/neuland/calibration/cal_to_hit/engine/R3BNeulandPredecessor.cxx index 965b01692..386dca783 100644 --- a/neuland/calibration/cal_to_hit/engine/R3BNeulandPredecessor.cxx +++ b/neuland/calibration/cal_to_hit/engine/R3BNeulandPredecessor.cxx @@ -91,14 +91,17 @@ namespace R3B::Neuland::Calibration void Predecessor::Init() { cal_to_hit_par_ = GetTask()->GetCal2HitPar(); } - void Predecessor::AddSignal(const BarCalData& signal) + void Predecessor::AddSignals(const std::vector& signals) { // all bar signal must have one signal on both sides - if (signal.left.size() != 1 or signal.right.size() != 1) + for (const auto& signal : signals) { - return; + if (signal.left.size() != 1 or signal.right.size() != 1) + { + continue; + } + fill_hist(signal); } - fill_hist(signal); } void Predecessor::fill_hist(const BarCalData& signal) diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandPredecessor.h b/neuland/calibration/cal_to_hit/engine/R3BNeulandPredecessor.h index 711a960fd..bf839c72d 100644 --- a/neuland/calibration/cal_to_hit/engine/R3BNeulandPredecessor.h +++ b/neuland/calibration/cal_to_hit/engine/R3BNeulandPredecessor.h @@ -19,7 +19,7 @@ namespace R3B::Neuland::Calibration // private virtual functions: void Init() override; - void AddSignal(const BarCalData& signal) override; + void AddSignals(const std::vector& signals) override; void Calibrate(Cal2HitPar& hit_par) override; void EndOfEvent(unsigned int event_num = 0) override {} void EventReset() override {} From 27bdfff2497ea436714689740ef30486735a96bb Mon Sep 17 00:00:00 2001 From: YanzhaoW Date: Sun, 16 Feb 2025 01:05:42 +0100 Subject: [PATCH 2/6] complete mille data processor --- .../R3BNeulandMilleCalDataProcessor.cxx | 70 ++++++++++++++++++- .../engine/R3BNeulandMilleCalDataProcessor.h | 26 +++++++ 2 files changed, 94 insertions(+), 2 deletions(-) diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx b/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx index 267012486..c2ff0a969 100644 --- a/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx +++ b/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx @@ -1,11 +1,18 @@ #include "R3BNeulandMilleCalDataProcessor.h" +#include +#include #include #include #include namespace R3B::Neuland::Calibration { - MilleDataProcessor::MilleDataProcessor(int num_of_modules) { init_data_registers(num_of_modules); } + MilleDataProcessor::MilleDataProcessor(int num_of_modules) + { + init_data_registers(num_of_modules); + fitter.SetFunction( + ROOT::Math::WrappedMultiTF1{ fit_function, static_cast(fit_function.GetNdim()) }, false); + } void MilleDataProcessor::init_data_registers(int num_of_modules) { @@ -25,6 +32,8 @@ namespace R3B::Neuland::Calibration } fit_result_.x_z = FitPar{}; fit_result_.y_z = FitPar{}; + x_z_vals.clear(); + y_z_vals.clear(); } auto MilleDataProcessor::operator()(const std::vector& signals) -> const auto& @@ -41,6 +50,7 @@ namespace R3B::Neuland::Calibration } remove_isolated_bar_signal(); + fit_planes(); return *this; } @@ -71,5 +81,61 @@ namespace R3B::Neuland::Calibration } } - void MilleDataProcessor::fit_planes() {} + void MilleDataProcessor::fit_planes() + { + fill_fit_data(); + + fit_plane_data(); + } + + void MilleDataProcessor::fill_fit_data() + { + for (auto& [plane_id, bar_data] : data_regsiters_) + { + 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(signal.module_num)); }) / + static_cast(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(0.); + fit_data.vals.push_back(displacement); + } + } + + void MilleDataProcessor::fit_plane_data() + { + const auto x_z_data = ROOT::Fit::BinData{ static_cast(x_z_vals.size()), + x_z_vals.z_vals.data(), + x_z_vals.vals.data(), + x_z_vals.z_errs.data(), + x_z_vals.errs.data() }; + const auto y_z_data = ROOT::Fit::BinData{ static_cast(y_z_vals.size()), + y_z_vals.z_vals.data(), + y_z_vals.vals.data(), + y_z_vals.z_errs.data(), + y_z_vals.errs.data() }; + fitter.Fit(x_z_data); + fit_result_.x_z.slope = fitter.Result().Parameter(0); + fit_result_.x_z.offset = fitter.Result().Parameter(1); + + fitter.Fit(y_z_data); + fit_result_.y_z.slope = fitter.Result().Parameter(0); + fit_result_.y_z.offset = fitter.Result().Parameter(1); + } + + auto MilleDataProcessor::calculate_residual(double z_val, double val, int module_num) const -> double + { + 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; + } } // namespace R3B::Neuland::Calibration diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.h b/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.h index 8ac2e3990..dbbaae7bf 100644 --- a/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.h +++ b/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.h @@ -1,6 +1,8 @@ #pragma once +#include #include +#include namespace R3B::Neuland::Calibration { @@ -32,19 +34,43 @@ namespace R3B::Neuland::Calibration FitPar x_z; FitPar y_z; }; + struct FitData + { + std::vector z_vals; + std::vector z_errs; + std::vector vals; + std::vector errs; + + void clear() + { + z_vals.clear(); + z_errs.clear(); + vals.clear(); + errs.clear(); + } + [[nodiscard]] auto size() const { return z_vals.size(); } + }; explicit MilleDataProcessor(int num_of_modules); auto operator()(const std::vector& signals) -> const auto&; [[nodiscard]] auto get_data() const -> const auto& { return data_regsiters_; } [[nodiscard]] auto get_fit() const -> const auto& { return fit_result_; } + [[nodiscard]] auto calculate_residual(double z_val, double val, int module_num) const -> double; void reset(); private: std::unordered_map> data_regsiters_; FitResult fit_result_; + FitData x_z_vals; + FitData y_z_vals; + + ROOT::Fit::Fitter fitter; + TF1 fit_function{ "mille_fitting", "[0] * x + [1]" }; void init_data_registers(int num_of_modules); void remove_isolated_bar_signal(); void fit_planes(); + void fill_fit_data(); + void fit_plane_data(); }; } // namespace R3B::Neuland::Calibration From 12f271fca1c6a9edf5057c0ee4a33d267cdf7c58 Mon Sep 17 00:00:00 2001 From: YanzhaoW Date: Fri, 21 Feb 2025 20:24:37 +0100 Subject: [PATCH 3/6] t_offset and effective light work --- .../cal_to_hit/R3BNeulandCalToHitPar.h | 4 + .../cal_to_hit/R3BNeulandCalToHitParTask.cxx | 6 +- .../cal_to_hit/R3BNeulandCalToHitParTask.h | 2 +- .../R3BNeulandMilleCalDataProcessor.cxx | 98 ++++--- .../engine/R3BNeulandMilleCalDataProcessor.h | 35 ++- .../cal_to_hit/engine/R3BNeulandMillepede.cxx | 240 ++++++++++++------ .../cal_to_hit/engine/R3BNeulandMillepede.h | 19 +- neuland/calibration/mille/Mille.cxx | 1 + neuland/calibration/mille/Mille.h | 1 + neuland/executables/neuland_map_to_cal.cxx | 2 +- neuland/executables/readme.md | 2 +- neuland/shared/R3BNeulandCommon.h | 6 +- neuland/util/neuland_read_hit_par.py | 7 + 13 files changed, 303 insertions(+), 120 deletions(-) diff --git a/neuland/calibration/cal_to_hit/R3BNeulandCalToHitPar.h b/neuland/calibration/cal_to_hit/R3BNeulandCalToHitPar.h index 99fb61489..0bceac509 100644 --- a/neuland/calibration/cal_to_hit/R3BNeulandCalToHitPar.h +++ b/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& { return module_pars_; diff --git a/neuland/calibration/cal_to_hit/R3BNeulandCalToHitParTask.cxx b/neuland/calibration/cal_to_hit/R3BNeulandCalToHitParTask.cxx index 081cf0634..4bb72ddc7 100644 --- a/neuland/calibration/cal_to_hit/R3BNeulandCalToHitParTask.cxx +++ b/neuland/calibration/cal_to_hit/R3BNeulandCalToHitParTask.cxx @@ -29,7 +29,6 @@ namespace R3B::Neuland : CalibrationTask(name, iVerbose) , cal_data_{ cal_data_name } , base_par_{ InputPar(base_par_name) } - , hit_par_{ OutputPar(hit_par_name) } // NOLINTEND { switch (method) @@ -37,14 +36,17 @@ namespace R3B::Neuland case Cal2HitParMethod::LSQT: R3BLOG(info, "Cal2HitPar method: LSQT."); engine_ = std::make_unique(); + hit_par_ = OutputPar(hit_par_name); break; - case Cal2HitParMethod::millipede: + case Cal2HitParMethod::millepede: R3BLOG(info, "Cal2HitPar method: Millepede."); engine_ = std::make_unique(); + hit_par_ = InputPar(hit_par_name); break; case Cal2HitParMethod::predecessor: R3BLOG(info, "Cal2HitPar method: predecessor."); engine_ = std::make_unique(); + hit_par_ = OutputPar(hit_par_name); break; } } diff --git a/neuland/calibration/cal_to_hit/R3BNeulandCalToHitParTask.h b/neuland/calibration/cal_to_hit/R3BNeulandCalToHitParTask.h index cd5a2151a..efc3a27f3 100644 --- a/neuland/calibration/cal_to_hit/R3BNeulandCalToHitParTask.h +++ b/neuland/calibration/cal_to_hit/R3BNeulandCalToHitParTask.h @@ -28,7 +28,7 @@ namespace R3B::Neuland { LSQT, predecessor, - millipede + millepede }; class Cal2HitParTask : public CalibrationTask diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx b/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx index c2ff0a969..3c70aa9c7 100644 --- a/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx +++ b/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx @@ -3,6 +3,7 @@ #include #include #include +#include #include namespace R3B::Neuland::Calibration @@ -10,8 +11,9 @@ 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(fit_function.GetNdim()) }, false); + fitter_.SetFunction( + ROOT::Math::WrappedMultiTF1{ fit_function_, static_cast(fit_function_.GetNdim()) }, true); + fitter_.Config().SetMinimizer("Linear"); } void MilleDataProcessor::init_data_registers(int num_of_modules) @@ -24,6 +26,12 @@ namespace R3B::Neuland::Calibration } } + 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_) @@ -32,11 +40,12 @@ namespace R3B::Neuland::Calibration } fit_result_.x_z = FitPar{}; fit_result_.y_z = FitPar{}; - x_z_vals.clear(); - y_z_vals.clear(); + x_z_vals_.clear(); + y_z_vals_.clear(); + reset_fitpars(); } - auto MilleDataProcessor::operator()(const std::vector& signals) -> const auto& + auto MilleDataProcessor::filter(const std::vector& signals) -> bool { // fill only the bar_cal_data with only one pmt signal on both sides for (const auto& signal : signals) @@ -50,9 +59,8 @@ namespace R3B::Neuland::Calibration } remove_isolated_bar_signal(); - fit_planes(); - return *this; + return fit_planes(); } void MilleDataProcessor::remove_isolated_bar_signal() @@ -81,19 +89,23 @@ namespace R3B::Neuland::Calibration } } - void MilleDataProcessor::fit_planes() + auto MilleDataProcessor::fit_planes() -> bool { fill_fit_data(); - fit_plane_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; + auto& fit_data = is_plane_horizontal ? y_z_vals_ : x_z_vals_; const auto z_val = PlaneID2ZPos(plane_id); const auto displacement = @@ -105,37 +117,63 @@ namespace R3B::Neuland::Calibration static_cast(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(0.); + fit_data.errs.push_back(BarSize_XY / 2.); fit_data.vals.push_back(displacement); } } - void MilleDataProcessor::fit_plane_data() + auto MilleDataProcessor::fit_plane_data() -> bool { - const auto x_z_data = ROOT::Fit::BinData{ static_cast(x_z_vals.size()), - x_z_vals.z_vals.data(), - x_z_vals.vals.data(), - x_z_vals.z_errs.data(), - x_z_vals.errs.data() }; - const auto y_z_data = ROOT::Fit::BinData{ static_cast(y_z_vals.size()), - y_z_vals.z_vals.data(), - y_z_vals.vals.data(), - y_z_vals.z_errs.data(), - y_z_vals.errs.data() }; - fitter.Fit(x_z_data); - fit_result_.x_z.slope = fitter.Result().Parameter(0); - fit_result_.x_z.offset = fitter.Result().Parameter(1); - - fitter.Fit(y_z_data); - fit_result_.y_z.slope = fitter.Result().Parameter(0); - fit_result_.y_z.offset = fitter.Result().Parameter(1); + 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 z_val, double val, int module_num) const -> double + 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(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 diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.h b/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.h index dbbaae7bf..13f638a98 100644 --- a/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.h +++ b/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.h @@ -2,6 +2,7 @@ #include #include +#include #include namespace R3B::Neuland::Calibration @@ -52,25 +53,43 @@ namespace R3B::Neuland::Calibration }; explicit MilleDataProcessor(int num_of_modules); - auto operator()(const std::vector& signals) -> const auto&; + auto filter(const std::vector& signals) -> bool; [[nodiscard]] auto get_data() const -> const auto& { return data_regsiters_; } [[nodiscard]] auto get_fit() const -> const auto& { return fit_result_; } - [[nodiscard]] auto calculate_residual(double z_val, double val, int module_num) const -> double; + [[nodiscard]] auto calculate_residual(double val, int module_num) const -> double; void reset(); + void set_p_value_cut(double val) { p_value_cut_ = val; } + private: + double p_value_cut_ = DEFAULT_CALIBRATION_P_VALUE_CUT; std::unordered_map> data_regsiters_; FitResult fit_result_; - FitData x_z_vals; - FitData y_z_vals; + FitData x_z_vals_; + FitData y_z_vals_; - ROOT::Fit::Fitter fitter; - TF1 fit_function{ "mille_fitting", "[0] * x + [1]" }; + ROOT::Fit::Fitter fitter_; + TF1 fit_function_{ "mille_fitting", "[0] * x + [1]" }; void init_data_registers(int num_of_modules); void remove_isolated_bar_signal(); - void fit_planes(); void fill_fit_data(); - void fit_plane_data(); + void reset_fitpars(); + auto fit_planes() -> bool; + auto fit_plane_data() -> bool; + auto linear_fit(const FitData& data, FitPar& fit_par) -> bool; }; } // namespace R3B::Neuland::Calibration + +template <> +class fmt::formatter +{ + public: + static constexpr auto parse(format_parse_context& ctx) { return ctx.end(); } + template + constexpr auto format(const R3B::Neuland::Calibration::MilleCalData& signal, FmtContent& ctn) const + { + return format_to( + ctn.out(), "ModuleNum: {}, left bar: {}, right bar: {}", signal.module_num, signal.left, signal.right); + } +}; diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandMillepede.cxx b/neuland/calibration/cal_to_hit/engine/R3BNeulandMillepede.cxx index 37989653f..cb87e0904 100644 --- a/neuland/calibration/cal_to_hit/engine/R3BNeulandMillepede.cxx +++ b/neuland/calibration/cal_to_hit/engine/R3BNeulandMillepede.cxx @@ -30,9 +30,10 @@ namespace rng = ranges; constexpr auto DEFAULT_RES_FILENAME = "millepede.res"; -constexpr auto SCALE_FACTOR = 100.F; +constexpr auto SCALE_FACTOR = 10.F; constexpr auto REFERENCE_BAR_NUM = 25; constexpr auto MILLE_BUFFER_SIZE = std::size_t{ 100000 }; +constexpr auto DEFAULT_T_ERROR = 2; // ns namespace { @@ -47,6 +48,15 @@ namespace } } } + + void change_time_offset(R3B::Neuland::Cal2HitPar& cal_to_hit_par) + { + auto& module_pars = cal_to_hit_par.GetListOfModuleParRef(); + for (auto& [module_num, module_par] : module_pars) + { + module_par.t_diff = module_par.t_diff * module_par.effective_speed; + } + } } // namespace namespace R3B::Neuland::Calibration @@ -60,7 +70,7 @@ namespace R3B::Neuland::Calibration pede_launcher_.set_parameter_filename(parameter_filename_); if (const auto* r3b_dir = std::getenv("R3BROOTPATH"); r3b_dir != nullptr) { - pede_launcher_.set_binary_dir(r3b_dir); + pede_launcher_.set_binary_dir(fmt::format("{}/bin", r3b_dir)); } else { @@ -68,6 +78,8 @@ namespace R3B::Neuland::Calibration "Environment variable R3BROOTPATH is not defined! Did you forget to source the \"config.sh\" file?"); } binary_data_writer_.set_buffer_size(MILLE_BUFFER_SIZE); + data_preprocessor_ = std::make_unique(GetModuleSize()); + data_preprocessor_->set_p_value_cut(p_value_cut_); init_steer_writer(); init_parameter(); @@ -83,21 +95,21 @@ namespace R3B::Neuland::Calibration switch (factor) { - case 0: - res.second = GlobalLabel::tsync; - break; - case 1: - res.second = GlobalLabel::offset_effective_c; - break; - case 2: - res.second = GlobalLabel::effective_c; - break; // case 0: - // res.second = GlobalLabel::offset_effective_c; + // res.second = GlobalLabel::tsync; // break; // case 1: + // res.second = GlobalLabel::offset_effective_c; + // break; + // case 2: // res.second = GlobalLabel::effective_c; // break; + case 0: + res.second = GlobalLabel::offset_effective_c; + break; + case 1: + res.second = GlobalLabel::effective_c; + break; default: throw R3B::logic_error(fmt::format("An error occured with unrecognized global par id: {}", par_num)); } @@ -110,24 +122,25 @@ namespace R3B::Neuland::Calibration const auto num_of_module = GetModuleSize(); switch (label) { - case GlobalLabel::tsync: - return module_num; - case GlobalLabel::offset_effective_c: - return module_num + num_of_module; - case GlobalLabel::effective_c: - return module_num + (2 * num_of_module); - // case GlobalLabel::offset_effective_c: + // case GlobalLabel::tsync: // return module_num; - // case GlobalLabel::effective_c: + // case GlobalLabel::offset_effective_c: // return module_num + num_of_module; + // case GlobalLabel::effective_c: + // return module_num + (2 * num_of_module); + case GlobalLabel::offset_effective_c: + return module_num; + case GlobalLabel::effective_c: + return module_num + num_of_module; default: - throw std::runtime_error("An error occured with unrecognized global tag"); + throw R3B::logic_error("An error occured with unrecognized global tag"); } } void MillepedeEngine::fill_module_parameters(const Millepede::ResultReader& result, Neuland::Cal2HitPar& cal_to_hit_par) { + change_time_offset(cal_to_hit_par); const auto& pars = result.get_pars(); for (const auto& [par_id, par] : pars) { @@ -143,12 +156,12 @@ namespace R3B::Neuland::Calibration break; case GlobalLabel::offset_effective_c: // The value here is the product of tDiff and effectiveSped. Real tDiff will be calculated later - par_ref.t_diff.value = par.value * SCALE_FACTOR; - par_ref.t_diff.error = par.error * SCALE_FACTOR; + par_ref.t_diff.value += par.value * SCALE_FACTOR; + par_ref.t_diff.error += par.error * SCALE_FACTOR; break; case GlobalLabel::effective_c: - par_ref.effective_speed.value = par.value; - par_ref.effective_speed.error = par.error; + par_ref.effective_speed.value += par.value; + par_ref.effective_speed.error += par.error; break; default: throw std::runtime_error("An error occured with unrecognized global tag"); @@ -213,7 +226,7 @@ namespace R3B::Neuland::Calibration return true; } - void MillepedeEngine::add_signal_t_sum(const BarCalData& signal) + void MillepedeEngine::add_signal_t_sum(const MilleCalData& signal) { buffer_clear(); const auto module_num = static_cast(signal.module_num); @@ -221,8 +234,8 @@ namespace R3B::Neuland::Calibration auto init_effective_c = cal_to_hit_par_->GetModuleParAt(module_num).effective_speed.value; - const auto& left_signal = signal.left.front(); - const auto& right_signal = signal.right.front(); + const auto& left_signal = signal.left; + const auto& right_signal = signal.right; const auto t_sum = (left_signal.leading_time - left_signal.trigger_time) + (right_signal.leading_time - right_signal.trigger_time) - average_t_sum_.value_or(0.F); @@ -239,50 +252,77 @@ namespace R3B::Neuland::Calibration write_to_buffer(); } - void MillepedeEngine::add_spacial_local_constraint(int module_num) + void MillepedeEngine::add_spacial_local_constraint(int plane_id, const std::vector& plane_signals) { buffer_clear(); - const auto plane_id = ModuleID2PlaneID(module_num - 1); const auto pos_z = PlaneID2ZPos(plane_id); const auto is_horizontal = IsPlaneIDHorizontal(plane_id); - const auto pos_bar_vert_disp = GetBarVerticalDisplacement(module_num); - // const auto local_derivs = is_horizontal ? std::array{ 0.F, pos_z / SCALE_FACTOR, 0.F, 1.F } - // : std::array{ pos_z / SCALE_FACTOR, 0.F, 1.F, 0.F }; - const auto local_derivs = is_horizontal ? std::array{ 0.F, pos_z / SCALE_FACTOR, 0.F, 0.F, 1.F, 0.F } - : std::array{ pos_z / SCALE_FACTOR, 0.F, 0.F, 1.F, 0.F, 0.F }; + // const auto pos_bar_vert_disp = GetBarVerticalDisplacement(module_num); + const auto pos_bar_vert_disp = std::accumulate(plane_signals.begin(), + plane_signals.end(), + 0., + [](double sum, const auto& signal) { + return sum + GetBarVerticalDisplacement(signal.module_num); + }) / + static_cast(plane_signals.size()); + const auto local_derivs = is_horizontal ? std::array{ 0.F, pos_z / SCALE_FACTOR, 0.F, 1.F } + : std::array{ pos_z / SCALE_FACTOR, 0.F, 1.F, 0.F }; + // const auto local_derivs = is_horizontal ? std::array{ 0.F, pos_z / SCALE_FACTOR, 0.F, 0.F, 1.F, 0.F } + // : std::array{ pos_z / SCALE_FACTOR, 0.F, 0.F, 1.F, 0.F, 0.F }; input_data_buffer_.measurement = static_cast(pos_bar_vert_disp / SCALE_FACTOR); input_data_buffer_.sigma = static_cast(BarSize_XY / SQRT_12 / SCALE_FACTOR * error_scale_factor_); + // if (not is_horizontal) + // { + // fmt::println("c_value_1.append({})\na_value_1.append({})\nb_value_1.append({})", + // pos_bar_vert_disp / SCALE_FACTOR, + // pos_z / SCALE_FACTOR, + // 1.); + // } std::copy(local_derivs.begin(), local_derivs.end(), std::back_inserter(input_data_buffer_.locals)); write_to_buffer(); } - void MillepedeEngine::add_signal_t_diff(const BarCalData& signal) + void MillepedeEngine::add_signal_t_diff(const MilleCalData& signal) { buffer_clear(); const auto module_num = static_cast(signal.module_num); const auto plane_id = ModuleID2PlaneID(static_cast(module_num) - 1); const auto is_horizontal = IsPlaneIDHorizontal(plane_id); - auto init_effective_c = cal_to_hit_par_->GetModuleParAt(module_num).effective_speed.value; + const auto& module_par = cal_to_hit_par_->GetModuleParAt(module_num); + auto init_effective_c = module_par.effective_speed.value; + auto init_t_offset = module_par.t_diff.value; const auto pos_z = static_cast(PlaneID2ZPos(plane_id)); - const auto& left_signal = signal.left.front(); - const auto& right_signal = signal.right.front(); + const auto& left_signal = signal.left; + const auto& right_signal = signal.right; const auto t_diff = (right_signal.leading_time - right_signal.trigger_time) - (left_signal.leading_time - left_signal.trigger_time); - input_data_buffer_.measurement = 0.F; + const auto t_error = t_diff.error == 0 ? DEFAULT_T_ERROR : t_diff.error; + input_data_buffer_.measurement = static_cast((init_effective_c * init_t_offset / 2. / SCALE_FACTOR) - + (init_effective_c * t_diff.value / SCALE_FACTOR / 2.)); input_data_buffer_.sigma = - static_cast(t_diff.error / SCALE_FACTOR / 2. * std::abs(init_effective_c) * error_scale_factor_); - // const auto local_derivs = is_horizontal ? std::array{ pos_z / SCALE_FACTOR, 0.F, 1.F, 0.F } - // : std::array{ 0.F, pos_z / SCALE_FACTOR, 0.F, 1.F }; - const auto local_derivs = is_horizontal ? std::array{ pos_z / SCALE_FACTOR, 0.F, 0.F, 1.F, 0.F, 0.F } - : std::array{ 0.F, pos_z / SCALE_FACTOR, 0.F, 0.F, 1.F, 0.F }; + static_cast(t_error / SCALE_FACTOR / 2. * std::abs(init_effective_c) * error_scale_factor_); + const auto local_derivs = is_horizontal ? std::array{ pos_z / SCALE_FACTOR, 0.F, 1.F, 0.F } + : std::array{ 0.F, pos_z / SCALE_FACTOR, 0.F, 1.F }; + // const auto local_derivs = is_horizontal ? std::array{ pos_z / SCALE_FACTOR, 0.F, 0.F, 1.F, 0.F, 0.F } + // : std::array{ 0.F, pos_z / SCALE_FACTOR, 0.F, 0.F, 1.F, 0.F }; std::copy(local_derivs.begin(), local_derivs.end(), std::back_inserter(input_data_buffer_.locals)); - input_data_buffer_.globals.emplace_back(get_global_label_id(module_num, GlobalLabel::offset_effective_c), 0.5F); + // fmt::println("Adding global: {}", get_global_label_id(module_num, GlobalLabel::offset_effective_c)); + input_data_buffer_.globals.emplace_back(get_global_label_id(module_num, GlobalLabel::offset_effective_c), + -0.5F); input_data_buffer_.globals.emplace_back(get_global_label_id(module_num, GlobalLabel::effective_c), static_cast(t_diff.value / SCALE_FACTOR / 2.)); + + // if (is_horizontal) + // { + // const auto c_val = (0.5 * (module_par.t_diff * module_par.effective_speed).value / SCALE_FACTOR) - + // (module_par.effective_speed.value * t_diff.value / SCALE_FACTOR / 2.); + // fmt::println( + // "c_value_2.append({})\na_value_2.append({})\nb_value_2.append({})", c_val, pos_z / SCALE_FACTOR, 1.); + // } write_to_buffer(); R3BLOG( debug, @@ -290,26 +330,71 @@ namespace R3B::Neuland::Calibration "Writting Mille data to binary file with meas = {} and z = {}", input_data_buffer_.measurement, pos_z)); } - void MillepedeEngine::AddSignals(const std::vector& signals) + auto MillepedeEngine::select_t_diff_signal(const std::vector& plane_data) { - for (const auto& signal : signals) + if (plane_data.size() < 1) + { + return plane_data.end(); + } + auto calculate_residual = [this](const MilleCalData& signal) -> double + { + const auto t_diff = (signal.right.leading_time - signal.right.trigger_time) - + (signal.left.leading_time - signal.left.trigger_time); + const auto& module_par = cal_to_hit_par_->GetModuleParAt(signal.module_num); + const auto position = (-t_diff + module_par.t_diff) / 2 * module_par.effective_speed; + const auto res = + data_preprocessor_->calculate_residual(position.value, static_cast(signal.module_num)); + return res; + }; + auto iter = std::min_element(plane_data.begin(), + plane_data.end(), + [calculate_residual](const auto& first, const auto& second) + { return calculate_residual(first) < calculate_residual(second); }); + if (iter != plane_data.end()) { - // all bar signal must have one signal on both sides - if (signal.left.size() != 1 or signal.right.size() != 1) + auto residual = calculate_residual(*iter); + hist_t_offset_residual_->Fill(residual); + if (residual > t_diff_residual_cut_) { - return; + return plane_data.end(); } + // fmt::println("selected signal residual: {}", residual); + } + return iter; + } - add_signal_t_sum(signal); - add_signal_t_diff(signal); - add_spacial_local_constraint(static_cast(signal.module_num)); + void MillepedeEngine::AddSignals(const std::vector& signals) + { + + // fmt::println("==================new event===============\n"); + if (not data_preprocessor_->filter(signals)) + { + return; + } + const auto& processed_data = data_preprocessor_->get_data(); + for (const auto& [plane_id, plane_signals] : + processed_data | + rng::views::filter([](const auto& planeid_signals) { return not planeid_signals.second.empty(); })) + { + // fmt::println("\n--------------------\n"); + auto iter = select_t_diff_signal(plane_signals); + if (iter == plane_signals.end()) + { + continue; + } + // fmt::println("selected signal: {}", *iter); + // for (const auto& signal : plane_signals) + // { + // // add_signal_t_sum(signal); + add_signal_t_diff(*iter); + add_spacial_local_constraint(plane_id, plane_signals); + // } } } void MillepedeEngine::Calibrate(Cal2HitPar& hit_par) { - hit_par.Reset(); - + binary_data_writer_.close(); R3BLOG(info, "Launching pede algorithm.."); pede_launcher_.launch(); pede_launcher_.end(); @@ -337,6 +422,7 @@ namespace R3B::Neuland::Calibration { // TODO: could be an empty event binary_data_writer_.end(); + data_preprocessor_->reset(); } void MillepedeEngine::EventReset() {} @@ -353,6 +439,10 @@ namespace R3B::Neuland::Calibration graph_effective_c_ = histograms.add_graph("effective_c", std::make_unique(module_size)); graph_effective_c_->SetTitle("Effective c vs BarNum"); + + static constexpr auto RESIDUAL_BIN_NUM = 500; + hist_t_offset_residual_ = histograms.add_hist( + "t_diff_residual", "Residual values of the positios calculated from t_dff", RESIDUAL_BIN_NUM, 0., 1000.); } void MillepedeEngine::buffer_clear() @@ -380,14 +470,14 @@ namespace R3B::Neuland::Calibration throw R3B::runtime_error("Pointer to cal_to_hit_par is nullptr!"); } - auto& module_pars = cal_to_hit_par_->GetListOfModuleParRef(); - module_pars.clear(); + // auto& module_pars = cal_to_hit_par_->GetListOfModuleParRef(); + // module_pars.clear(); - for (unsigned int module_num{ 1 }; module_num <= num_of_modules; ++module_num) - { - auto module_par_iter = module_pars.try_emplace(module_num).first; - module_par_iter->second.effective_speed.value = init_effective_c_; - } + // for (unsigned int module_num{ 1 }; module_num <= num_of_modules; ++module_num) + // { + // auto module_par_iter = module_pars.try_emplace(module_num).first; + // module_par_iter->second.effective_speed.value = init_effective_c_; + // } } void MillepedeEngine::init_steer_writer() @@ -400,14 +490,22 @@ namespace R3B::Neuland::Calibration steer_writer.add_other_options(std::vector{ "hugecut", "50000" }); steer_writer.add_other_options(std::vector{ "outlierdownweighting", "4" }); - const auto module_size = GetModuleSize(); - for (int module_num{ 1 }; module_num <= module_size; ++module_num) - { - steer_writer.add_parameter_default(get_global_label_id(module_num, GlobalLabel::effective_c), - std::make_pair(init_effective_c_, 0.F)); - } - steer_writer.add_parameter_default(get_global_label_id(REFERENCE_BAR_NUM, GlobalLabel::tsync), - std::make_pair(0.F, -1.F)); + // const auto module_size = GetModuleSize(); + // for (int module_num{ 1 }; module_num <= module_size; ++module_num) + // { + // const auto& module_par = cal_to_hit_par_->GetModuleParAt(module_num); + // steer_writer.add_parameter_default( + // get_global_label_id(module_num, GlobalLabel::effective_c), + // // std::make_pair(module_par.effective_speed.value, module_par.effective_speed.error)); + // std::make_pair(DEFAULT_EFFECTIVE_C, module_par.effective_speed.error)); + + // const auto offset_effective_c = module_par.t_diff * module_par.effective_speed; + // steer_writer.add_parameter_default( + // get_global_label_id(module_num, GlobalLabel::offset_effective_c), + // std::make_pair(offset_effective_c.value / SCALE_FACTOR, offset_effective_c.error / SCALE_FACTOR)); + // } + // steer_writer.add_parameter_default(get_global_label_id(REFERENCE_BAR_NUM, GlobalLabel::tsync), + // std::make_pair(0.F, -1.F)); steer_writer.write(); } } // namespace R3B::Neuland::Calibration diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandMillepede.h b/neuland/calibration/cal_to_hit/engine/R3BNeulandMillepede.h index d029cf515..ea5891d17 100644 --- a/neuland/calibration/cal_to_hit/engine/R3BNeulandMillepede.h +++ b/neuland/calibration/cal_to_hit/engine/R3BNeulandMillepede.h @@ -13,6 +13,7 @@ #pragma once #include "R3BNeulandCosmicEngine.h" +#include "R3BNeulandMilleCalDataProcessor.h" #include #include #include @@ -24,7 +25,7 @@ class TGraphErrors; namespace R3B::Neuland::Calibration { - enum class GlobalLabel + enum class GlobalLabel : int8_t { tsync, // tsync offset_effective_c, // offset times effective_C @@ -36,6 +37,8 @@ namespace R3B::Neuland::Calibration public: MillepedeEngine() = default; void enable_rank_check(bool rank_check = true) { has_rank_check_ = rank_check; } + void set_t_diff_residual_cut(double val) { t_diff_residual_cut_ = val; } + void set_p_value_cut(double val) { p_value_cut_ = val; } private: bool has_rank_check_ = false; @@ -45,6 +48,9 @@ namespace R3B::Neuland::Calibration // float smallest_time_sum_ = 0.; std::optional average_t_sum_; float init_effective_c_ = DEFAULT_EFFECTIVE_C; + double t_diff_residual_cut_ = DEFAULT_T_DIFF_RESIDUAL_CUT; + double p_value_cut_ = DEFAULT_CALIBRATION_P_VALUE_CUT; + MilleDataPoint input_data_buffer_; std::string input_data_filename_ = "neuland_cosmic_mille.bin"; std::string pede_steer_filename_ = "neuland_steer.txt"; @@ -58,10 +64,13 @@ namespace R3B::Neuland::Calibration TGraphErrors* graph_time_offset_ = nullptr; TGraphErrors* graph_time_sync_ = nullptr; TGraphErrors* graph_effective_c_ = nullptr; + TH1D* hist_t_offset_residual_ = nullptr; // parameter: Cal2HitPar* cal_to_hit_par_ = nullptr; + std::unique_ptr data_preprocessor_; + void Init() override; void AddSignals(const std::vector& signals) override; void Calibrate(Cal2HitPar& hit_par) override; @@ -79,9 +88,9 @@ namespace R3B::Neuland::Calibration void buffer_clear(); void write_to_buffer(); - void add_signal_t_sum(const BarCalData& signal); - void add_signal_t_diff(const BarCalData& signal); - void add_spacial_local_constraint(int module_num); + void add_signal_t_sum(const MilleCalData& signal); + void add_signal_t_diff(const MilleCalData& signal); + void add_spacial_local_constraint(int plane_id, const std::vector& plane_signals); auto set_minimum_values(const std::vector& signals) -> bool; inline auto get_global_label_id(int module_num, GlobalLabel label) -> int; inline auto to_module_num_label(int par_num) -> std::pair; @@ -90,6 +99,8 @@ namespace R3B::Neuland::Calibration void init_parameter(); void init_steer_writer(); + + auto select_t_diff_signal(const std::vector& plane_data); }; } // namespace R3B::Neuland::Calibration diff --git a/neuland/calibration/mille/Mille.cxx b/neuland/calibration/mille/Mille.cxx index 09fd96900..10aea3f50 100644 --- a/neuland/calibration/mille/Mille.cxx +++ b/neuland/calibration/mille/Mille.cxx @@ -132,6 +132,7 @@ namespace R3B output_file_ << fmt::format("{}\n", fmt::join(buffer_.get_values(), " ")); } + void Mille::close() { output_file_.close(); } void Mille::check_buffer_size(std::size_t nLocal, std::size_t nGlobal) { if (buffer_.get_current_size() >= max_buffer_size_) diff --git a/neuland/calibration/mille/Mille.h b/neuland/calibration/mille/Mille.h index 4132afb29..8b436ffc9 100644 --- a/neuland/calibration/mille/Mille.h +++ b/neuland/calibration/mille/Mille.h @@ -47,6 +47,7 @@ namespace R3B has_special_done_ = false; } void end(); + void close(); private: bool has_special_done_ = false; // if true, special(..) already called for this record diff --git a/neuland/executables/neuland_map_to_cal.cxx b/neuland/executables/neuland_map_to_cal.cxx index 618102c54..95c05ba52 100644 --- a/neuland/executables/neuland_map_to_cal.cxx +++ b/neuland/executables/neuland_map_to_cal.cxx @@ -110,7 +110,7 @@ auto main(int argc, char** argv) -> int } auto cal2hit_method = - (enable_mille.value()) ? R3B::Neuland::Cal2HitParMethod::millipede : R3B::Neuland::Cal2HitParMethod::LSQT; + (enable_mille.value()) ? R3B::Neuland::Cal2HitParMethod::millepede : R3B::Neuland::Cal2HitParMethod::LSQT; auto cal2hitParTask = std::make_unique(cal2hit_method); auto* cal2hitParTaskPtr = cal2hitParTask.get(); cal2hitParTaskPtr->SetTrigger(R3B::Neuland::CalTrigger::all); diff --git a/neuland/executables/readme.md b/neuland/executables/readme.md index 8ba7e58fb..121de8eb2 100644 --- a/neuland/executables/readme.md +++ b/neuland/executables/readme.md @@ -158,7 +158,7 @@ neuland ana [-h] [options] - `NeulandNeutronsRValue`. See `R3BNeulandNeutronsRValue`. - `NeulandCal2HitParTask`: See `R3B::Neuland::Cal2HitParTask`. - `min-stat`: Minimal number of hits from the events that are used for the calibration. - - `method`: Method of the calibration. Available options: `LSQT`, `predecessor` and `millepede`. See `R3B::neuland::Cal2HitParMethod`. + - `method`: Method of the calibration. Available options: `LSQT`, `predecessor` and `millepede`. See `R3B::Neuland::Cal2HitParMethod`. All tasks listed above have three common options: "enable", "read" and "write". The "enable" option specifies whether the task is added or not. If some tasks should be added, simply change its value to `true`. diff --git a/neuland/shared/R3BNeulandCommon.h b/neuland/shared/R3BNeulandCommon.h index fda3a01d4..a1cd0f5f4 100644 --- a/neuland/shared/R3BNeulandCommon.h +++ b/neuland/shared/R3BNeulandCommon.h @@ -35,8 +35,10 @@ namespace R3B::Neuland return val * val; } - // Miscellaneous defaults: - constexpr auto DEFAULT_EFFECTIVE_C = 8.; // cm/ns + // Millepede calibration defaults: + constexpr auto DEFAULT_EFFECTIVE_C = 8.; // cm/ns + constexpr auto DEFAULT_CALIBRATION_P_VALUE_CUT = 1e-10; // any smaller values will be discarded + constexpr auto DEFAULT_T_DIFF_RESIDUAL_CUT = 400; // any larger values will be discarded // Initialize variables from Birk' s Law diff --git a/neuland/util/neuland_read_hit_par.py b/neuland/util/neuland_read_hit_par.py index 3b9dff965..7d4bf505e 100644 --- a/neuland/util/neuland_read_hit_par.py +++ b/neuland/util/neuland_read_hit_par.py @@ -2,6 +2,8 @@ Checking the values of cal_to_hit parameters in the root file. The corresponding pandas dataframe can be retrieved. """ +import warnings + import numpy as np import pandas as pd import ROOT @@ -101,6 +103,11 @@ def _fill_data(self, hit_par): num_of_modules = hit_par.GetNumOfModules() for bar_id in range(num_of_modules): self._par_dict["bar_id"][bar_id] = bar_id + if not hit_par.HasModuleParAt(bar_id + 1): + warnings.warn( + f"Parameter with module num {bar_id + 1} doesn't exist in the file!" + ) + continue module_par = hit_par.GetModuleParAt(bar_id + 1) for par_name in self._par_name_list: dicts = dir(getattr(ROOT.R3B.Neuland.HitModulePar(), par_name)) From 421fba19a07e025d385642b8e265809d4935f1fc Mon Sep 17 00:00:00 2001 From: YanzhaoW Date: Tue, 25 Feb 2025 14:42:42 +0100 Subject: [PATCH 4/6] add clangd config and fix include warnings --- .clangd | 2 + .../engine/R3BNeulandLSQREngineAdaptor.cxx | 12 ++++-- .../engine/R3BNeulandLSQREngineAdaptor.h | 4 ++ .../R3BNeulandMilleCalDataProcessor.cxx | 7 ++++ .../engine/R3BNeulandMilleCalDataProcessor.h | 3 ++ .../cal_to_hit/engine/R3BNeulandMillepede.cxx | 40 ++++++++++++------- .../cal_to_hit/engine/R3BNeulandMillepede.h | 12 ++++++ .../engine/R3BNeulandPredecessor.cxx | 11 +++++ .../cal_to_hit/engine/R3BNeulandPredecessor.h | 8 ++++ 9 files changed, 82 insertions(+), 17 deletions(-) create mode 100644 .clangd diff --git a/.clangd b/.clangd new file mode 100644 index 000000000..c37e381f6 --- /dev/null +++ b/.clangd @@ -0,0 +1,2 @@ +Diagnostics: + MissingIncludes: Strict diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandLSQREngineAdaptor.cxx b/neuland/calibration/cal_to_hit/engine/R3BNeulandLSQREngineAdaptor.cxx index 6b87541ee..7adc77e08 100644 --- a/neuland/calibration/cal_to_hit/engine/R3BNeulandLSQREngineAdaptor.cxx +++ b/neuland/calibration/cal_to_hit/engine/R3BNeulandLSQREngineAdaptor.cxx @@ -12,13 +12,19 @@ ******************************************************************************/ #include "R3BNeulandLSQREngineAdaptor.h" -#include "R3BDataMonitor.h" -#include "R3BNeulandHitPar.h" + #include +#include #include +#include +#include +#include +#include +#include #include #include -#include +#include +#include namespace { diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandLSQREngineAdaptor.h b/neuland/calibration/cal_to_hit/engine/R3BNeulandLSQREngineAdaptor.h index aa3580036..118de9084 100644 --- a/neuland/calibration/cal_to_hit/engine/R3BNeulandLSQREngineAdaptor.h +++ b/neuland/calibration/cal_to_hit/engine/R3BNeulandLSQREngineAdaptor.h @@ -14,9 +14,13 @@ #pragma once #include "R3BNeulandCosmicEngine.h" + +#include +#include #include #include #include +#include namespace R3B::Neuland::Calibration { diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx b/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx index 3c70aa9c7..1f33320a9 100644 --- a/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx +++ b/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx @@ -1,10 +1,17 @@ #include "R3BNeulandMilleCalDataProcessor.h" +#include "R3BLogger.h" +#include "R3BNeulandCalData2.h" #include #include #include +#include #include +#include #include +#include #include +#include +#include namespace R3B::Neuland::Calibration { diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.h b/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.h index 13f638a98..5ecc7be02 100644 --- a/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.h +++ b/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.h @@ -4,6 +4,9 @@ #include #include #include +#include +#include +#include namespace R3B::Neuland::Calibration { diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandMillepede.cxx b/neuland/calibration/cal_to_hit/engine/R3BNeulandMillepede.cxx index cb87e0904..8b1c1837d 100644 --- a/neuland/calibration/cal_to_hit/engine/R3BNeulandMillepede.cxx +++ b/neuland/calibration/cal_to_hit/engine/R3BNeulandMillepede.cxx @@ -12,26 +12,46 @@ ******************************************************************************/ #include "R3BNeulandMillepede.h" +#include "ParResultReader.h" +#include "R3BDataMonitor.h" +#include "R3BLogger.h" +#include "R3BNeulandCalData2.h" +#include "R3BNeulandCalToHitPar.h" +#include "R3BNeulandMilleCalDataProcessor.h" #include #include #include #include #include +#include +#include +#include +#include +#include +#include +#include #include +#include #include #include #include #include +#include +#include +#include +#include +#include #include +#include #include namespace rng = ranges; constexpr auto DEFAULT_RES_FILENAME = "millepede.res"; constexpr auto SCALE_FACTOR = 10.F; -constexpr auto REFERENCE_BAR_NUM = 25; +// constexpr auto REFERENCE_BAR_NUM = 25; constexpr auto MILLE_BUFFER_SIZE = std::size_t{ 100000 }; constexpr auto DEFAULT_T_ERROR = 2; // ns @@ -332,7 +352,7 @@ namespace R3B::Neuland::Calibration auto MillepedeEngine::select_t_diff_signal(const std::vector& plane_data) { - if (plane_data.size() < 1) + if (plane_data.empty()) { return plane_data.end(); } @@ -463,21 +483,10 @@ namespace R3B::Neuland::Calibration void MillepedeEngine::init_parameter() { - auto num_of_modules = GetModuleSize(); - if (cal_to_hit_par_ == nullptr) { throw R3B::runtime_error("Pointer to cal_to_hit_par is nullptr!"); } - - // auto& module_pars = cal_to_hit_par_->GetListOfModuleParRef(); - // module_pars.clear(); - - // for (unsigned int module_num{ 1 }; module_num <= num_of_modules; ++module_num) - // { - // auto module_par_iter = module_pars.try_emplace(module_num).first; - // module_par_iter->second.effective_speed.value = init_effective_c_; - // } } void MillepedeEngine::init_steer_writer() @@ -486,7 +495,10 @@ namespace R3B::Neuland::Calibration steer_writer.set_filepath(pede_steer_filename_); steer_writer.set_parameter_file(parameter_filename_); steer_writer.set_data_filepath(input_data_filename_); - steer_writer.add_method(SteerWriter::Method::inversion, std::make_pair(3.F, 0.001F)); + static constexpr auto NUMBER_OF_ITERARTION = 3.F; + static constexpr auto CONVERGENCE_RECOGNITION = 0.001F; + steer_writer.add_method(SteerWriter::Method::inversion, + std::make_pair(NUMBER_OF_ITERARTION, CONVERGENCE_RECOGNITION)); steer_writer.add_other_options(std::vector{ "hugecut", "50000" }); steer_writer.add_other_options(std::vector{ "outlierdownweighting", "4" }); diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandMillepede.h b/neuland/calibration/cal_to_hit/engine/R3BNeulandMillepede.h index ea5891d17..2601d4378 100644 --- a/neuland/calibration/cal_to_hit/engine/R3BNeulandMillepede.h +++ b/neuland/calibration/cal_to_hit/engine/R3BNeulandMillepede.h @@ -12,13 +12,25 @@ ******************************************************************************/ #pragma once +#include "MilleEntry.h" +#include "R3BDataMonitor.h" +#include "R3BLogger.h" +#include "R3BNeulandCalData2.h" +#include "R3BNeulandCalToHitPar.h" #include "R3BNeulandCosmicEngine.h" #include "R3BNeulandMilleCalDataProcessor.h" #include #include #include #include +#include +#include +#include +#include #include +#include +#include +#include // #include class TGraphErrors; diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandPredecessor.cxx b/neuland/calibration/cal_to_hit/engine/R3BNeulandPredecessor.cxx index 386dca783..05484d38e 100644 --- a/neuland/calibration/cal_to_hit/engine/R3BNeulandPredecessor.cxx +++ b/neuland/calibration/cal_to_hit/engine/R3BNeulandPredecessor.cxx @@ -1,8 +1,19 @@ #include "R3BNeulandPredecessor.h" +#include "R3BDataMonitor.h" +#include "R3BLogger.h" +#include "R3BNeulandCalData2.h" +#include "R3BNeulandCalToHitPar.h" +#include "R3BNeulandCommon.h" +#include "R3BValueError.h" #include #include #include #include +#include +#include +#include +#include +#include namespace R3B::Neuland::Calibration { diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandPredecessor.h b/neuland/calibration/cal_to_hit/engine/R3BNeulandPredecessor.h index bf839c72d..1d86ac8e6 100644 --- a/neuland/calibration/cal_to_hit/engine/R3BNeulandPredecessor.h +++ b/neuland/calibration/cal_to_hit/engine/R3BNeulandPredecessor.h @@ -1,7 +1,15 @@ #pragma once +#include "R3BDataMonitor.h" +#include "R3BLogger.h" +#include "R3BNeulandCalData2.h" +#include "R3BNeulandCalToHitPar.h" #include "R3BNeulandCosmicEngine.h" +#include +#include +#include + namespace R3B::Neuland::Calibration { class Predecessor : public CosmicEngineInterface From 1369865ba907cb2e25f223e961877fe67805fab1 Mon Sep 17 00:00:00 2001 From: YanzhaoW Date: Tue, 25 Feb 2025 14:55:36 +0100 Subject: [PATCH 5/6] fix further header issues --- .clangd | 2 ++ .../cal_to_hit/R3BNeulandCalToHitParTask.cxx | 10 ++++++++++ .../engine/R3BNeulandMilleCalDataProcessor.cxx | 2 -- 3 files changed, 12 insertions(+), 2 deletions(-) diff --git a/.clangd b/.clangd index c37e381f6..fcf7dc0a6 100644 --- a/.clangd +++ b/.clangd @@ -1,2 +1,4 @@ Diagnostics: + UnusedIncludes: Strict MissingIncludes: Strict + AnalyzeAngledIncludes: Yes diff --git a/neuland/calibration/cal_to_hit/R3BNeulandCalToHitParTask.cxx b/neuland/calibration/cal_to_hit/R3BNeulandCalToHitParTask.cxx index 4bb72ddc7..f260be014 100644 --- a/neuland/calibration/cal_to_hit/R3BNeulandCalToHitParTask.cxx +++ b/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 +#include #include #include #include #include +#include +#include namespace R3B::Neuland { diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx b/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx index 1f33320a9..5b20ff760 100644 --- a/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx +++ b/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx @@ -7,9 +7,7 @@ #include #include #include -#include #include -#include #include #include From 47404b6468fa4fe96191e937492b0053f2081570 Mon Sep 17 00:00:00 2001 From: YanzhaoW Date: Tue, 25 Feb 2025 15:09:34 +0100 Subject: [PATCH 6/6] fix fmt ranges issue --- .../cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx | 3 +++ 1 file changed, 3 insertions(+) diff --git a/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx b/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx index 5b20ff760..fb0adf00d 100644 --- a/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx +++ b/neuland/calibration/cal_to_hit/engine/R3BNeulandMilleCalDataProcessor.cxx @@ -11,6 +11,9 @@ #include #include +// NOLINTNEXTLINE +#include + namespace R3B::Neuland::Calibration { MilleDataProcessor::MilleDataProcessor(int num_of_modules)