mirror of
https://github.com/slsdetectorgroup/aare.git
synced 2026-06-04 18:08:43 +02:00
Merge branch 'main' into dev/strixels/remap_simple
This commit is contained in:
@@ -0,0 +1,154 @@
|
||||
// SPDX-License-Identifier: MPL-2.0
|
||||
#pragma once
|
||||
#include "aare/Models.hpp"
|
||||
#include "aare/NDView.hpp"
|
||||
#include <Minuit2/FCNGradientBase.h>
|
||||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <stdexcept>
|
||||
#include <vector>
|
||||
|
||||
namespace aare {
|
||||
|
||||
namespace func {
|
||||
|
||||
/**
|
||||
* @brief Generic chi-squared FCN with analytic gradient for a 1D model.
|
||||
*
|
||||
* @tparam Model A model struct that satisfies:
|
||||
* - static constexpr std::size_t npar;
|
||||
* - static double eval(double x, const std::vector<double>& par);
|
||||
* - static void eval_and_grad(double x, const std::vector<double>& par,
|
||||
* double& f, std::array<double, npar>& g);
|
||||
* - static bool is_valid(const std::vector<double>& par);
|
||||
*
|
||||
* Gradient:
|
||||
* d(chi2)/dp_k = -2 * sum_i w_i * (y_i - f_i) * df_i/dp_k
|
||||
*
|
||||
* where w_i = 1/sigma_i^2 (weighted) or 1 (unweighted).
|
||||
*
|
||||
* By providing analytic gradients we avoid 2*npar extra function evaluations
|
||||
* per Minuit step that would otherwise be spent on finite differences.
|
||||
*
|
||||
* @throws std::invalid_argument if par.size() != Model::npar.
|
||||
*
|
||||
* Invalid model parameters do not throw; they return a large penalty
|
||||
* (and a zero gradient fallback) so the minimizer can remain in control.
|
||||
*/
|
||||
template <class Model>
|
||||
class Chi2Model1DGrad : public ROOT::Minuit2::FCNGradientBase {
|
||||
public:
|
||||
Chi2Model1DGrad(NDView<double, 1> x, NDView<double, 1> y)
|
||||
: x_(x), y_(y), s_(), weighted_(false) {}
|
||||
|
||||
Chi2Model1DGrad(NDView<double, 1> x, NDView<double, 1> y,
|
||||
NDView<double, 1> y_err)
|
||||
: x_(x), y_(y), s_(y_err), weighted_(true) {}
|
||||
|
||||
~Chi2Model1DGrad() override = default;
|
||||
|
||||
double operator()(const std::vector<double> &par) const override {
|
||||
if (par.size() != Model::npar) {
|
||||
throw std::invalid_argument(
|
||||
"Chi2Model1DGrad: wrong parameter vector size.");
|
||||
}
|
||||
|
||||
if (!Model::is_valid(par))
|
||||
return 1e20;
|
||||
|
||||
double chi2 = 0.0;
|
||||
|
||||
if (weighted_) {
|
||||
for (ssize_t i = 0; i < x_.size(); ++i) {
|
||||
const double si = s_[i];
|
||||
if (si == 0.0)
|
||||
continue;
|
||||
|
||||
const double f_i = Model::eval(x_[i], par);
|
||||
const double r_i = y_[i] - f_i;
|
||||
chi2 += (r_i * r_i) / (si * si);
|
||||
}
|
||||
} else {
|
||||
for (ssize_t i = 0; i < x_.size(); ++i) {
|
||||
const double f_i = Model::eval(x_[i], par);
|
||||
const double r_i = y_[i] - f_i;
|
||||
chi2 += r_i * r_i;
|
||||
}
|
||||
}
|
||||
|
||||
return chi2;
|
||||
}
|
||||
|
||||
std::vector<double>
|
||||
Gradient(const std::vector<double> &par) const override {
|
||||
if (par.size() != Model::npar) {
|
||||
throw std::invalid_argument(
|
||||
"Chi2Model1DGrad: wrong parameter vector size.");
|
||||
}
|
||||
|
||||
std::vector<double> grad(Model::npar, 0.0);
|
||||
if (!Model::is_valid(par))
|
||||
return grad;
|
||||
|
||||
std::array<double, Model::npar> df{};
|
||||
double f_i = 0.0;
|
||||
|
||||
if (weighted_) {
|
||||
for (ssize_t i = 0; i < x_.size(); ++i) {
|
||||
const double si = s_[i];
|
||||
if (si == 0.0)
|
||||
continue;
|
||||
|
||||
Model::eval_and_grad(x_[i], par, f_i, df);
|
||||
|
||||
const double r_i = y_[i] - f_i;
|
||||
const double c = -2.0 * r_i / (si * si);
|
||||
|
||||
for (std::size_t k = 0; k < Model::npar; ++k) {
|
||||
grad[k] += c * df[k];
|
||||
}
|
||||
}
|
||||
} else {
|
||||
for (ssize_t i = 0; i < x_.size(); ++i) {
|
||||
Model::eval_and_grad(x_[i], par, f_i, df);
|
||||
|
||||
const double r_i = y_[i] - f_i;
|
||||
const double c = -2.0 * r_i;
|
||||
|
||||
for (std::size_t k = 0; k < Model::npar; ++k) {
|
||||
grad[k] += c * df[k];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return grad;
|
||||
}
|
||||
|
||||
/** @brief Error definition: 1.0 for chi-squared (delta_chi2 = 1 ->
|
||||
* 1-sigma). */
|
||||
double Up() const override { return 1.0; }
|
||||
|
||||
private:
|
||||
NDView<double, 1> x_;
|
||||
NDView<double, 1> y_;
|
||||
NDView<double, 1> s_;
|
||||
bool weighted_;
|
||||
};
|
||||
|
||||
// ── Convenient aliases ──────────────────────────────────────────────
|
||||
|
||||
using Chi2Gaussian = Chi2Model1DGrad<aare::model::Gaussian>;
|
||||
using Chi2GaussianErfcPlateau =
|
||||
Chi2Model1DGrad<aare::model::GaussianErfcPlateau>;
|
||||
using Chi2GaussianChargeSharing =
|
||||
Chi2Model1DGrad<aare::model::GaussianChargeSharing>;
|
||||
using Chi2GaussianChargeSharingKb =
|
||||
Chi2Model1DGrad<aare::model::GaussianChargeSharingKb>;
|
||||
using Chi2RisingScurve = Chi2Model1DGrad<aare::model::RisingScurve>;
|
||||
using Chi2FallingScurve = Chi2Model1DGrad<aare::model::FallingScurve>;
|
||||
using Chi2Pol1 = Chi2Model1DGrad<aare::model::Pol1>;
|
||||
using Chi2Pol2 = Chi2Model1DGrad<aare::model::Pol2>;
|
||||
|
||||
} // namespace func
|
||||
|
||||
} // namespace aare
|
||||
@@ -17,7 +17,7 @@ template <class ItemType> class CircularFifo {
|
||||
aare::ProducerConsumerQueue<ItemType> filled_slots;
|
||||
|
||||
public:
|
||||
CircularFifo() : CircularFifo(100){};
|
||||
CircularFifo() : CircularFifo(100) {};
|
||||
CircularFifo(uint32_t size)
|
||||
: fifo_size(size), free_slots(size + 1), filled_slots(size + 1) {
|
||||
|
||||
|
||||
@@ -383,7 +383,7 @@ ClusterFile<ClusterType, Enable>::read_frame_without_cut() {
|
||||
"Unexpected error (not feof or ferror) when reading frame number");
|
||||
}
|
||||
|
||||
int32_t n_clusters; // Saved as 32bit integer in the cluster file
|
||||
uint32_t n_clusters;
|
||||
if (fread(&n_clusters, sizeof(n_clusters), 1, fp) != 1) {
|
||||
throw std::runtime_error(LOCATION +
|
||||
"Could not read number of clusters");
|
||||
|
||||
@@ -23,7 +23,7 @@ template <typename ClusterType = Cluster<int32_t, 3, 3>,
|
||||
typename = std::enable_if_t<no_2x2_cluster<ClusterType>::value>>
|
||||
class ClusterFinder {
|
||||
Shape<2> m_image_size;
|
||||
const PEDESTAL_TYPE m_nSigma;
|
||||
PEDESTAL_TYPE m_nSigma;
|
||||
const PEDESTAL_TYPE c2;
|
||||
const PEDESTAL_TYPE c3;
|
||||
Pedestal<PEDESTAL_TYPE> m_pedestal;
|
||||
@@ -53,6 +53,10 @@ class ClusterFinder {
|
||||
<< ", nSigma: " << nSigma << ", capacity: " << capacity;
|
||||
}
|
||||
|
||||
void set_nSigma(PEDESTAL_TYPE nSigma) { m_nSigma = nSigma; }
|
||||
|
||||
PEDESTAL_TYPE get_nSigma() const { return m_nSigma; }
|
||||
|
||||
void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) {
|
||||
m_pedestal.push(frame);
|
||||
}
|
||||
|
||||
@@ -284,6 +284,23 @@ class ClusterFinderMT {
|
||||
// auto rc = m_input_queue.write(std::move(frame));
|
||||
// fmt::print("pushed frame {}\n", rc);
|
||||
// }
|
||||
|
||||
/**
|
||||
* @brief Set the nSigma value for all the cluster finders.
|
||||
* @param nSigma number of sigma above the pedestal to consider a photon
|
||||
* during cluster finding.
|
||||
*/
|
||||
void set_nSigma(const PEDESTAL_TYPE nSigma) {
|
||||
// Wait for all queues to be empty before changing the sigma
|
||||
for (auto &q : m_input_queues) {
|
||||
while (!q->isEmpty()) {
|
||||
std::this_thread::sleep_for(m_default_wait);
|
||||
}
|
||||
}
|
||||
for (auto &cf : m_cluster_finders) {
|
||||
cf->set_nSigma(nSigma);
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
} // namespace aare
|
||||
@@ -59,6 +59,27 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
|
||||
other.m_data.clear();
|
||||
}
|
||||
|
||||
/**
|
||||
* @brief Create a copy of the clustervector by filtering clusters in the
|
||||
* ClusterVector using a boolean mask.
|
||||
* @param mask boolean 1d mask
|
||||
* @return ClusterVector containing only the clusters where the mask is
|
||||
* true
|
||||
*/
|
||||
ClusterVector operator()(NDView<bool, 1> mask) {
|
||||
if (static_cast<size_t>(mask.size()) != m_data.size()) {
|
||||
throw std::runtime_error(
|
||||
LOCATION + "Mask size does not match number of clusters");
|
||||
}
|
||||
ClusterVector result(capacity(), frame_number());
|
||||
for (size_t i = 0; i < m_data.size(); ++i) {
|
||||
if (mask(i)) {
|
||||
result.push_back(m_data[i]);
|
||||
}
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
// Move assignment operator
|
||||
ClusterVector &operator=(ClusterVector &&other) noexcept {
|
||||
if (this != &other) {
|
||||
|
||||
+229
-14
@@ -5,7 +5,17 @@
|
||||
#include <fmt/core.h>
|
||||
#include <vector>
|
||||
|
||||
#include "aare/Chi2.hpp"
|
||||
#include "aare/FitModel.hpp"
|
||||
#include "aare/NDArray.hpp"
|
||||
#include "aare/utils/par.hpp"
|
||||
#include "aare/utils/task.hpp"
|
||||
|
||||
#include "Minuit2/FunctionMinimum.h"
|
||||
#include "Minuit2/MnHesse.h"
|
||||
#include "Minuit2/MnMigrad.h"
|
||||
#include "Minuit2/MnPrint.h"
|
||||
#include "Minuit2/MnUserParameters.h"
|
||||
|
||||
namespace aare {
|
||||
|
||||
@@ -24,20 +34,6 @@ NDArray<double, 1> scurve2(NDView<double, 1> x, NDView<double, 1> par);
|
||||
|
||||
} // namespace func
|
||||
|
||||
/**
|
||||
* @brief Estimate the initial parameters for a Gaussian fit
|
||||
*/
|
||||
std::array<double, 3> gaus_init_par(const NDView<double, 1> x,
|
||||
const NDView<double, 1> y);
|
||||
|
||||
std::array<double, 2> pol1_init_par(const NDView<double, 1> x,
|
||||
const NDView<double, 1> y);
|
||||
|
||||
std::array<double, 6> scurve_init_par(const NDView<double, 1> x,
|
||||
const NDView<double, 1> y);
|
||||
std::array<double, 6> scurve2_init_par(const NDView<double, 1> x,
|
||||
const NDView<double, 1> y);
|
||||
|
||||
static constexpr int DEFAULT_NUM_THREADS = 4;
|
||||
|
||||
/**
|
||||
@@ -118,4 +114,223 @@ void fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
|
||||
NDView<double, 3> y_err, NDView<double, 3> par_out,
|
||||
NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
|
||||
int n_threads);
|
||||
|
||||
// Minuit2 fit_pixel / fit_3d object based API
|
||||
|
||||
// _____________________________________________________________________
|
||||
//
|
||||
// fit_pixel — single-pixel minimisation
|
||||
// _____________________________________________________________________
|
||||
|
||||
/**
|
||||
* @brief Fit a single pixel's data using Minuit2.
|
||||
*
|
||||
* The caller provides a thread-local clone of MnUserParameters so that
|
||||
* no heap allocation happens here (only SetValue/SetError stores).
|
||||
*
|
||||
* User-precedence rules:
|
||||
* - Fixed parameters: untouched (value and fixed flag preserved from clone).
|
||||
* - User-set start: value preserved, step size auto-filled.
|
||||
* - Neither: both value and step size auto-filled from data.
|
||||
*
|
||||
* @tparam Model Model struct (Gaussian, RisingScurve, …).
|
||||
* @tparam FCN Chi2 functor type (Chi2Model1D or Chi2Model1DGrad
|
||||
* instantiation).
|
||||
*
|
||||
* @param model The FitModel configuration (read-only).
|
||||
* @param upar_local Thread-local clone of model.upar(). Modified in place.
|
||||
* @param x Scan points (shared across all pixels).
|
||||
* @param y Measured values for this pixel.
|
||||
* @param y_err Per-point uncertainties (empty view -> unweighted fit).
|
||||
*
|
||||
* @return NDArray<double,1> of size:
|
||||
* - compute_errors: [p0..pN, err0..errN, chi2] -> 2*npar + 1
|
||||
* - otherwise: [p0..pN, chi2] -> npar + 1
|
||||
*/
|
||||
template <typename Model, typename FCN>
|
||||
NDArray<double, 1> fit_pixel(const FitModel<Model> &model,
|
||||
ROOT::Minuit2::MnUserParameters &upar_local,
|
||||
NDView<double, 1> x, NDView<double, 1> y,
|
||||
NDView<double, 1> y_err) {
|
||||
|
||||
constexpr std::size_t npar = Model::npar;
|
||||
const bool want_errors = model.compute_errors();
|
||||
const ssize_t result_size = want_errors ? (2 * npar + 1) : (npar + 1);
|
||||
|
||||
// ──── automatic parameter estimation ─────────────
|
||||
auto start = Model::estimate_par(x, y);
|
||||
|
||||
// dead / degenerate pixel guard
|
||||
if (!Model::is_valid(std::vector<double>(start.begin(), start.end()))) {
|
||||
return NDArray<double, 1>({result_size}, 0.0);
|
||||
}
|
||||
|
||||
// ──── data-range statistics for step sizes ─────────────
|
||||
double x_range, y_range, slope_scale;
|
||||
model::compute_ranges(x, y, x_range, y_range, slope_scale);
|
||||
|
||||
std::array<double, npar> steps{};
|
||||
Model::compute_steps(start, x_range, y_range, slope_scale, steps);
|
||||
|
||||
// ── apply auto-estimates respecting user precedence ─────────────
|
||||
for (std::size_t i = 0; i < npar; ++i) {
|
||||
// fixed: do not touch at all
|
||||
if (model.is_user_fixed(i)) {
|
||||
continue;
|
||||
}
|
||||
|
||||
if (!model.is_user_start(i)) {
|
||||
upar_local.SetValue(i, start[i]);
|
||||
}
|
||||
|
||||
upar_local.SetError(i, steps[i]);
|
||||
}
|
||||
|
||||
// ──── build functor ────────
|
||||
auto chi2 = (y_err.size() > 0) ? FCN(x, y, y_err) : FCN(x, y);
|
||||
|
||||
// ──── run minimizer ────────
|
||||
ROOT::Minuit2::MnMigrad migrad(chi2, upar_local, model.strategy());
|
||||
ROOT::Minuit2::FunctionMinimum min =
|
||||
migrad(model.max_calls(), model.tolerance());
|
||||
|
||||
if (!min.IsValid())
|
||||
return NDArray<double, 1>({result_size}, 0.0);
|
||||
|
||||
// ──── pack results ────────
|
||||
if (want_errors) {
|
||||
ROOT::Minuit2::MnHesse hesse;
|
||||
hesse(chi2, min);
|
||||
|
||||
const auto &values = min.UserState().Params();
|
||||
const auto &errors = min.UserState().Errors();
|
||||
|
||||
NDArray<double, 1> result({result_size});
|
||||
for (std::size_t k = 0; k < npar; ++k) {
|
||||
result[k] = values[k];
|
||||
result[npar + k] = errors[k];
|
||||
}
|
||||
result[2 * npar] = min.Fval();
|
||||
return result;
|
||||
}
|
||||
|
||||
const auto &values = min.UserState().Params();
|
||||
NDArray<double, 1> result({result_size});
|
||||
for (std::size_t k = 0; k < npar; ++k)
|
||||
result[k] = values[k];
|
||||
result[npar] = min.Fval();
|
||||
return result;
|
||||
}
|
||||
|
||||
// ── self-contained for 1D / standalone use ─────────
|
||||
template <typename Model, typename FCN>
|
||||
NDArray<double, 1> fit_pixel(const FitModel<Model> &model, NDView<double, 1> x,
|
||||
NDView<double, 1> y, NDView<double, 1> y_err) {
|
||||
auto upar_local = model.upar();
|
||||
return fit_pixel<Model, FCN>(model, upar_local, x, y, y_err);
|
||||
}
|
||||
|
||||
// Overload: uncertainties not provided
|
||||
template <typename Model, typename FCN>
|
||||
NDArray<double, 1> fit_pixel(const FitModel<Model> &model, NDView<double, 1> x,
|
||||
NDView<double, 1> y) {
|
||||
auto upar_local = model.upar();
|
||||
return fit_pixel<Model, FCN>(model, upar_local, x, y, NDView<double, 1>{});
|
||||
}
|
||||
|
||||
// _____________________________________________________________________
|
||||
//
|
||||
// fit_3d — row-parallel fitting over (rows, cols) pixel grid
|
||||
// _____________________________________________________________________
|
||||
/**
|
||||
* @brief Fit all pixels in a 3D data cube (rows x cols x n_scan).
|
||||
*
|
||||
* @tparam Model Model struct.
|
||||
* @tparam FCN Chi2 functor type.
|
||||
*
|
||||
* @param model Fit configuration shared by all pixels.
|
||||
* @param x Scan points, shape `(n_scan)`.
|
||||
* @param y Measured values, shape `(rows, cols, n_scan)`.
|
||||
* @param y_err Uncertainties, same shape as y, or empty for unweighted
|
||||
* fits.
|
||||
* @param par_out Output parameters, shape `(rows, cols, npar)`.
|
||||
* @param err_out Output parameter errors, shape `(rows, cols, npar)`, if
|
||||
* used.
|
||||
* @param chi2_out Output chi-squared / objective values, shape `(rows,
|
||||
* cols)`.
|
||||
* @param n_threads Number of threads used to split rows.
|
||||
*
|
||||
*/
|
||||
template <typename Model, typename FCN>
|
||||
void fit_3d(
|
||||
const FitModel<Model> &model, NDView<double, 1> x, // (n_scan)
|
||||
NDView<double, 3> y, // (rows, cols, n_scan)
|
||||
NDView<double, 3> y_err, // (rows, cols, n_scan) or empty for unweighted fit
|
||||
NDView<double, 3> par_out, NDView<double, 3> err_out,
|
||||
NDView<double, 2> chi2_out, int n_threads) {
|
||||
const std::size_t npar = Model::npar;
|
||||
|
||||
// ──── checks ───────
|
||||
if (x.size() != y.shape(2))
|
||||
throw std::runtime_error("fit_3d: x.size() must match y.shape(2).");
|
||||
|
||||
if (par_out.shape(0) != y.shape(0) || par_out.shape(1) != y.shape(1) ||
|
||||
par_out.shape(2) != npar)
|
||||
throw std::runtime_error("par_out must have shape [rows, cols, npar].");
|
||||
|
||||
if (chi2_out.shape(0) != y.shape(0) || chi2_out.shape(1) != y.shape(1))
|
||||
throw std::runtime_error("chi2_out must have shape [rows, cols].");
|
||||
|
||||
const bool has_errors = (y_err.size() > 0);
|
||||
const bool want_par_errors = (err_out.size() > 0) && model.compute_errors();
|
||||
|
||||
if (has_errors) {
|
||||
if (y.shape(0) != y_err.shape(0) || y.shape(1) != y_err.shape(1) ||
|
||||
y.shape(2) != y_err.shape(2))
|
||||
throw std::runtime_error(
|
||||
"fit_3d: y and y_err must have identical shape.");
|
||||
|
||||
if (err_out.shape(0) != y.shape(0) || err_out.shape(1) != y.shape(1) ||
|
||||
err_out.shape(2) != npar)
|
||||
throw std::runtime_error(
|
||||
"err_out must have shape [rows, cols, npar].");
|
||||
}
|
||||
|
||||
// ──── parallel dispatch ───────
|
||||
auto process = [&](ssize_t first_row, ssize_t last_row) {
|
||||
// one clone per thread
|
||||
auto upar_local = model.upar();
|
||||
|
||||
for (ssize_t row = first_row; row < last_row; row++) {
|
||||
for (ssize_t col = 0; col < y.shape(1); col++) {
|
||||
|
||||
NDView<double, 1> values(&y(row, col, 0), {y.shape(2)});
|
||||
NDView<double, 1> errors =
|
||||
has_errors ? NDView<double, 1>(&y_err(row, col, 0),
|
||||
{y_err.shape(2)})
|
||||
: NDView<double, 1>{};
|
||||
|
||||
auto res =
|
||||
fit_pixel<Model, FCN>(model, upar_local, x, values, errors);
|
||||
|
||||
for (std::size_t k = 0; k < npar; ++k) {
|
||||
par_out(row, col, k) = res(k);
|
||||
}
|
||||
|
||||
if (want_par_errors) {
|
||||
for (std::size_t k = 0; k < npar; ++k) {
|
||||
err_out(row, col, k) = res(npar + k);
|
||||
}
|
||||
chi2_out(row, col) = res(2 * npar);
|
||||
} else {
|
||||
chi2_out(row, col) = res(npar);
|
||||
}
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
auto tasks = split_task(0, static_cast<int>(y.shape(0)), n_threads);
|
||||
RunInParallel(process, tasks);
|
||||
}
|
||||
|
||||
} // namespace aare
|
||||
@@ -0,0 +1,137 @@
|
||||
// SPDX-License-Identifier: MPL-2.0
|
||||
#pragma once
|
||||
|
||||
#include "aare/Models.hpp"
|
||||
#include <type_traits>
|
||||
|
||||
#include "Minuit2/MnStrategy.h"
|
||||
#include "Minuit2/MnUserParameters.h"
|
||||
|
||||
namespace aare {
|
||||
|
||||
template <typename Model> class FitModel {
|
||||
ROOT::Minuit2::MnUserParameters upar_;
|
||||
ROOT::Minuit2::MnStrategy strategy_;
|
||||
unsigned int max_calls_;
|
||||
double tolerance_;
|
||||
bool compute_errors_;
|
||||
|
||||
std::array<bool, Model::npar> user_fixed_{};
|
||||
std::array<bool, Model::npar> user_start_{};
|
||||
|
||||
/** @brief Safely resolve a parameter name to its index. */
|
||||
unsigned int checked_index(const std::string &name) const {
|
||||
for (std::size_t i = 0; i < npar; ++i) {
|
||||
if (upar_.Name(i) == name)
|
||||
return static_cast<unsigned int>(i);
|
||||
}
|
||||
throw std::runtime_error("FitModel: unknown parameter name '" + name +
|
||||
"'");
|
||||
}
|
||||
|
||||
public:
|
||||
static constexpr std::size_t npar = Model::npar;
|
||||
|
||||
/**
|
||||
* @brief Construct a fit model with sensible defaults.
|
||||
*
|
||||
* @param strategy Minuit2 strategy level (0 = fast/gradient, 1 =
|
||||
* default).
|
||||
* @param max_calls Maximum FCN calls per pixel minimisation.
|
||||
* @param tolerance Minuit2 EDM tolerance.
|
||||
* @param compute_errors If true, run MnHesse after minimisation.
|
||||
*/
|
||||
FitModel(unsigned int strategy = 0, unsigned int max_calls = 100,
|
||||
double tolerance = 0.5, bool compute_errors = false)
|
||||
: strategy_(strategy), max_calls_(max_calls), tolerance_(tolerance),
|
||||
compute_errors_(compute_errors) {
|
||||
for (std::size_t i = 0; i < npar; ++i) {
|
||||
const auto pi = Model::param_info[i];
|
||||
const bool has_lo = std::isfinite(pi.default_lo);
|
||||
const bool has_hi = std::isfinite(pi.default_hi);
|
||||
|
||||
// Add parameters and valid bounds
|
||||
if (has_lo && has_hi) {
|
||||
upar_.Add(pi.name, 0.0, 1.0, pi.default_lo, pi.default_hi);
|
||||
} else if (has_lo) {
|
||||
upar_.Add(pi.name, 0.0, 1.0, pi.default_lo, 1e6);
|
||||
} else {
|
||||
upar_.Add(pi.name, 0.0, 1.0);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/** @brief Set lower and upper bounds for parameter idx.*/
|
||||
void SetParLimits(unsigned int idx, double lo, double hi) {
|
||||
upar_.SetLimits(idx, lo, hi);
|
||||
}
|
||||
|
||||
/**
|
||||
* @brief Fix parameter idx at value val.
|
||||
*
|
||||
* Excluded from minimisation. Automatic estimates will not touch it.
|
||||
*/
|
||||
void FixParameter(unsigned int idx, double val) {
|
||||
SetParameter(idx, val);
|
||||
upar_.Fix(idx);
|
||||
user_fixed_[idx] = true;
|
||||
}
|
||||
|
||||
/** @brief Release a previously fixed parameter, re-enabling auto estimates.
|
||||
*/
|
||||
void ReleaseParameter(unsigned int idx) {
|
||||
upar_.Release(idx);
|
||||
user_fixed_[idx] = false;
|
||||
}
|
||||
|
||||
void ReleaseParameter(const std::string &name) {
|
||||
ReleaseParameter(checked_index(name));
|
||||
}
|
||||
|
||||
/** @brief Set an explicit starting value for parameter idx.*/
|
||||
void SetParameter(unsigned int idx, double val) {
|
||||
upar_.SetValue(idx, val);
|
||||
user_start_[idx] = true;
|
||||
}
|
||||
|
||||
void SetParameter(const std::string &name, double val) {
|
||||
// go through index to maintain user_start_ bookkeeping
|
||||
SetParameter(checked_index(name), val);
|
||||
}
|
||||
|
||||
void FixParameter(const std::string &name, double val) {
|
||||
// go through index to maintain user_fixed_ bookkeeping
|
||||
FixParameter(checked_index(name), val);
|
||||
}
|
||||
|
||||
void SetParLimits(const std::string &name, double lo, double hi) {
|
||||
SetParLimits(checked_index(name), lo, hi);
|
||||
}
|
||||
|
||||
std::string GetParName(unsigned int idx) const {
|
||||
return upar_.GetName(idx);
|
||||
}
|
||||
|
||||
std::vector<std::string> GetParNames() const {
|
||||
std::vector<std::string> names;
|
||||
for (std::size_t i = 0; i < npar; ++i)
|
||||
names.push_back(GetParName(i));
|
||||
return names;
|
||||
}
|
||||
static constexpr std::size_t GetNpar() noexcept { return npar; }
|
||||
|
||||
void SetMaxCalls(unsigned int n) { max_calls_ = n; }
|
||||
void SetTolerance(double t) { tolerance_ = t; }
|
||||
void SetComputeErrors(bool b) { compute_errors_ = b; }
|
||||
|
||||
// accessors
|
||||
const ROOT::Minuit2::MnUserParameters &upar() const { return upar_; }
|
||||
const ROOT::Minuit2::MnStrategy &strategy() const { return strategy_; }
|
||||
unsigned int max_calls() const { return max_calls_; }
|
||||
double tolerance() const { return tolerance_; }
|
||||
bool compute_errors() const { return compute_errors_; }
|
||||
bool is_user_fixed(unsigned int idx) const { return user_fixed_[idx]; }
|
||||
bool is_user_start(unsigned int idx) const { return user_start_[idx]; }
|
||||
};
|
||||
|
||||
} // namespace aare
|
||||
@@ -93,6 +93,26 @@ class Interpolator {
|
||||
std::vector<Photon>
|
||||
interpolate(const ClusterVector<ClusterType> &clusters) const;
|
||||
|
||||
/**
|
||||
* @brief interpolates the cluster centers for all clusters to a better
|
||||
* precision
|
||||
* @param clusters clusters of photon hits to interpolate
|
||||
* @param etas precomputed eta values for each cluster (must be in the same
|
||||
* order as the clusters)
|
||||
* @return interpolated photons (photon positions are given as double but
|
||||
* following row column format e.g. x=0, y=0 means top row and first column
|
||||
* of frame) (An interpolated photon position of (1.5, 2.5) corresponds to
|
||||
* an estimated photon hit at the pixel center of pixel (1,2))
|
||||
*/
|
||||
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
|
||||
typename CoordType = uint16_t,
|
||||
typename Enable = std::enable_if_t<is_cluster_v<
|
||||
Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>>>
|
||||
std::vector<Photon> interpolate(
|
||||
const ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>
|
||||
&clusters,
|
||||
const std::vector<Eta2<T>> &etas) const;
|
||||
|
||||
/**
|
||||
* @brief transforms the eta values to uniform coordinates based on the CDF
|
||||
* ieta_x and ieta_y
|
||||
@@ -171,11 +191,14 @@ Coordinate2D Interpolator::transform_eta_values(const Eta2<T> &eta) const {
|
||||
|
||||
if (static_cast<ssize_t>(ix) >= m_etabinsx.size() - 1 ||
|
||||
static_cast<ssize_t>(iy) >= m_etabinsy.size() - 1 ||
|
||||
static_cast<ssize_t>(ie) >= m_energy_bins.size() - 1)
|
||||
throw std::runtime_error(
|
||||
fmt::format("Eta values out of bounds of eta distribution: eta.x = "
|
||||
"{:.4f}, eta.y = {:.4f}, energy = {:.4f}",
|
||||
eta.x, eta.y, eta.sum));
|
||||
static_cast<ssize_t>(ie) >= m_energy_bins.size() - 1) {
|
||||
throw std::runtime_error(fmt::format(
|
||||
"Eta values (eta.x = {:4f}, eta.y = {:4f}, energy = {}) out of "
|
||||
"bounds of eta distribution with largest values: (eta.x = {:4f}, "
|
||||
"eta.y = {:4f}, energy = {:4f})",
|
||||
eta.x, eta.y, eta.sum, *(m_etabinsx.end() - 1),
|
||||
*(m_etabinsy.end() - 1), *(m_energy_bins.end() - 1)));
|
||||
}
|
||||
|
||||
// TODO: bilinear interpolation only works if all bins have a size > 1 -
|
||||
// otherwise bilinear interpolation with zero values which skew the
|
||||
@@ -203,16 +226,16 @@ Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) const {
|
||||
photon.y = cluster.y;
|
||||
photon.energy = static_cast<decltype(photon.energy)>(eta.sum);
|
||||
|
||||
Coordinate2D uniform_coordinates{};
|
||||
|
||||
try {
|
||||
// check if eta values are within bounds
|
||||
transform_eta_values(eta);
|
||||
uniform_coordinates = transform_eta_values(eta);
|
||||
} catch (const std::runtime_error &e) {
|
||||
throw std::runtime_error(
|
||||
fmt::format("{} for cluster: {}", e.what(), cluster_index));
|
||||
}
|
||||
|
||||
auto uniform_coordinates = transform_eta_values(eta);
|
||||
|
||||
if (EtaFunction == &calculate_eta2<typename ClusterType::value_type,
|
||||
ClusterType::cluster_size_x,
|
||||
ClusterType::cluster_size_y,
|
||||
@@ -262,4 +285,48 @@ Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) const {
|
||||
return photons;
|
||||
}
|
||||
|
||||
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
|
||||
typename CoordType, typename Enable>
|
||||
std::vector<Photon> Interpolator::interpolate(
|
||||
const ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>
|
||||
&clusters,
|
||||
const std::vector<Eta2<T>> &etas) const {
|
||||
|
||||
if (clusters.size() != etas.size()) {
|
||||
throw std::runtime_error(
|
||||
fmt::format("Size of clusters and precomputed etas must be the "
|
||||
"same, but got {} clusters and {} etas",
|
||||
clusters.size(), etas.size()));
|
||||
}
|
||||
|
||||
std::vector<Photon> photons;
|
||||
photons.reserve(clusters.size());
|
||||
|
||||
for (size_t i = 0; i < clusters.size(); ++i) {
|
||||
const auto &cluster = clusters[i];
|
||||
const auto &eta = etas[i];
|
||||
|
||||
Photon photon;
|
||||
photon.x = cluster.x;
|
||||
photon.y = cluster.y;
|
||||
photon.energy = static_cast<decltype(photon.energy)>(eta.sum);
|
||||
|
||||
Coordinate2D uniform_coordinates{};
|
||||
try {
|
||||
// check if eta values are within bounds
|
||||
uniform_coordinates = transform_eta_values(eta);
|
||||
} catch (const std::runtime_error &e) {
|
||||
throw std::runtime_error(
|
||||
fmt::format("{} for cluster: {}", e.what(), i));
|
||||
}
|
||||
|
||||
photon.x += uniform_coordinates.x;
|
||||
photon.y += uniform_coordinates.y;
|
||||
|
||||
photons.push_back(photon);
|
||||
}
|
||||
|
||||
return photons;
|
||||
}
|
||||
|
||||
} // namespace aare
|
||||
File diff suppressed because it is too large
Load Diff
@@ -36,7 +36,7 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
|
||||
* @brief Default constructor. Constructs an empty NDArray.
|
||||
*
|
||||
*/
|
||||
NDArray() : shape_(), strides_(c_strides<Ndim>(shape_)), data_(nullptr){};
|
||||
NDArray() : shape_(), strides_(c_strides<Ndim>(shape_)), data_(nullptr) {};
|
||||
|
||||
/**
|
||||
* @brief Construct a new NDArray object with a given shape.
|
||||
|
||||
@@ -63,7 +63,7 @@ template <class T> struct ProducerConsumerQueue {
|
||||
return *this;
|
||||
}
|
||||
|
||||
ProducerConsumerQueue() : ProducerConsumerQueue(2){};
|
||||
ProducerConsumerQueue() : ProducerConsumerQueue(2) {};
|
||||
// size must be >= 2.
|
||||
//
|
||||
// Also, note that the number of usable slots in the queue at any
|
||||
|
||||
@@ -45,7 +45,11 @@ template <typename T> class VarClusterFinder {
|
||||
{-1, -1, 0, 1}}; // col ### 8-neighbour by scaning from top to bottom
|
||||
const std::array<int, 8> di_{{0, 0, -1, 1, -1, 1, -1, 1}}; // row
|
||||
const std::array<int, 8> dj_{{-1, 1, 0, 0, 1, -1, -1, 1}}; // col
|
||||
std::map<int, int> child; // heirachy: key: child; val: parent
|
||||
std::map<int, int> child; // heirachy: key: child; val: parent
|
||||
int numberOfNeighbours = 8; // 4 or 8
|
||||
bool empty_surroundingPixels =
|
||||
true; // whether to set peripheral pixels to 0, to avoid potential
|
||||
// influence for pedestal updating
|
||||
std::unordered_map<int, Hit> h_size;
|
||||
std::vector<Hit> hits;
|
||||
// std::vector<std::vector<int16_t>> row
|
||||
@@ -67,6 +71,16 @@ template <typename T> class VarClusterFinder {
|
||||
void set_peripheralThresholdFactor(int factor) {
|
||||
peripheralThresholdFactor_ = factor;
|
||||
}
|
||||
void set_numberOfNeighbours(int num) {
|
||||
if (num == 4 || num == 8)
|
||||
numberOfNeighbours = num;
|
||||
else
|
||||
throw std::runtime_error(
|
||||
LOCATION + "number of neighbours should be either 4 or 8");
|
||||
}
|
||||
void set_empty_surroundingPixels(bool empty) {
|
||||
empty_surroundingPixels = empty;
|
||||
}
|
||||
void find_clusters(NDView<T, 2> img);
|
||||
void find_clusters_X(NDView<T, 2> img);
|
||||
void rec_FillHit(int clusterIndex, int i, int j);
|
||||
@@ -201,7 +215,7 @@ void VarClusterFinder<T>::rec_FillHit(int clusterIndex, int i, int j) {
|
||||
}
|
||||
original_(i, j) = 0;
|
||||
|
||||
for (int k = 0; k < 8; ++k) { // 8 for 8-neighbour
|
||||
for (int k = 0; k < numberOfNeighbours; ++k) { //
|
||||
const auto row = i + di_[k];
|
||||
const auto col = j + dj_[k];
|
||||
if (row >= 0 && col >= 0 && row < shape_[0] && col < shape_[1]) {
|
||||
@@ -218,9 +232,11 @@ void VarClusterFinder<T>::rec_FillHit(int clusterIndex, int i, int j) {
|
||||
// h_size[clusterIndex].enes[h_size[clusterIndex].size] =
|
||||
// original_(row, col);
|
||||
// }// ? weather to include peripheral pixels
|
||||
original_(row, col) =
|
||||
0; // remove peripheral pixels, to avoid potential influence
|
||||
// for pedestal updating
|
||||
if (empty_surroundingPixels) {
|
||||
original_(row, col) =
|
||||
0; // remove peripheral pixels, to avoid potential
|
||||
// influence for pedestal updating
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@@ -67,7 +67,7 @@ class Logger {
|
||||
|
||||
public:
|
||||
Logger() = default;
|
||||
explicit Logger(TLogLevel level) : m_level(level){};
|
||||
explicit Logger(TLogLevel level) : m_level(level) {};
|
||||
~Logger() {
|
||||
// output in the destructor to allow for << syntax
|
||||
os << RESET << '\n';
|
||||
|
||||
@@ -1,11 +1,11 @@
|
||||
// SPDX-License-Identifier: MPL-2.0
|
||||
#pragma once
|
||||
#include "aare/NDView.hpp"
|
||||
#include "aare/utils/task.hpp"
|
||||
#include <thread>
|
||||
#include <utility>
|
||||
#include <vector>
|
||||
|
||||
#include "aare/utils/task.hpp"
|
||||
|
||||
namespace aare {
|
||||
|
||||
template <typename F>
|
||||
|
||||
Reference in New Issue
Block a user