12#ifndef ALICEO2_MATHUTILS_MATHBASE_H_
13#define ALICEO2_MATHUTILS_MATHBASE_H_
26#include "TLinearFitter.h"
31#include "HFitInterface.h"
32#include "TFitResultPtr.h"
33#include "TFitResult.h"
34#include "Fit/Fitter.h"
35#include "Fit/BinData.h"
36#include "Math/WrappedMultiTF1.h"
37#include <Math/SMatrix.h>
38#include <Math/SVector.h>
60TFitResultPtr
fit(
const size_t nBins,
const T* arr,
const T xMin,
const T xMax,
TF1&
func, std::string_view option =
"")
63 ROOT::Fit::FitOptionsMake(ROOT::Fit::EFitObjectType::kHistogram, option.data(), fitOption);
65 ROOT::Fit::DataRange
range(xMin, xMax);
66 ROOT::Fit::DataOptions opt;
67 ROOT::Fit::BinData fitdata(opt,
range);
68 fitdata.Initialize(nBins, 1);
71 std::shared_ptr<TFitResult> tfr(
new TFitResult());
74 ROOT::Fit::Fitter fitter(tfr);
77 const double binWidth = double(xMax - xMin) / double(nBins);
79 for (Int_t ibin = 0; ibin < nBins; ibin++) {
80 const double x = double(xMin) + double(ibin + 0.5) * binWidth;
81 const double y = double(arr[ibin]);
82 const double ey = std::sqrt(
y);
83 fitdata.Add(
x,
y, ey);
86 const int special =
func.GetNumber();
87 const int npar =
func.GetNpar();
88 bool linear =
func.IsLinear();
89 if (special == 299 + npar) {
93 if (fitOption.Bound || fitOption.Like || fitOption.Errors || fitOption.Gradient || fitOption.More || fitOption.User || fitOption.Integral || fitOption.Minuit) {
97 if (special != 0 && !fitOption.Bound && !linear) {
99 ROOT::Fit::InitGaus(fitdata, &
func);
100 }
else if (special == 400) {
101 ROOT::Fit::InitGaus(fitdata, &
func);
102 }
else if (special == 200) {
103 ROOT::Fit::InitExpo(fitdata, &
func);
107 if ((linear || fitOption.Gradient)) {
108 fitter.SetFunction(ROOT::Math::WrappedMultiTF1(
func));
110 fitter.SetFunction(
static_cast<const ROOT::Math::IParamMultiFunction&
>(ROOT::Math::WrappedMultiTF1(
func)));
114 const bool fitok = fitter.Fit(fitdata, fitOption.ExecPolicy);
116 LOGP(warning,
"bad fit");
119 return TFitResultPtr(tfr);
133bool medmadGaus(
size_t nBins,
const T* arr,
const T xMin,
const T xMax, std::array<double, 3>&
param)
135 int bStart = 0, bEnd = -1;
136 double sum = 0, binW = double(xMax - xMin) / nBins, medVal = xMin;
137 for (
int i = 0;
i < (
int)nBins;
i++) {
151 double cum = 0, thresh = 0.5 *
sum, frac0 = 0;
152 int bid = bStart, prevbid = bid;
157 frac0 = 1. + (thresh - cum) /
float(arr[bid]);
158 medVal = xMin + binW * (bid + frac0);
159 int bdiff = bid - prevbid - 1;
161 medVal -= bdiff * binW * 0.5;
171 double edgeL = frac0 + bid, edgeR = edgeL, dist = 0., wL = 0, wR = 0;
174 int bL = edgeL, bR = edgeR;
175 if (edgeL > bStart) {
182 wR = 1. + bR - edgeR;
187 auto wdt = std::min(wL, wR);
189 wdt = std::max(wL, wR);
195 dist += wdt * (cum - thresh) / amp * 0.5;
203 constexpr double SQRT2PI = 2.5066283;
205 param[2] = dist * binW * 1.4826;
232Double_t
fitGaus(
const size_t nBins,
const T* arr,
const T xMin,
const T xMax, std::vector<T>&
param)
234 static TLinearFitter fitter(3,
"pol2");
235 static TMatrixD mat(3, 3);
236 static Double_t kTol = mat.GetTol();
237 fitter.StoreData(kFALSE);
238 fitter.ClearPoints();
243 T rms = TMath::RMS(nBins, arr);
244 T
max = TMath::MaxElement(nBins, arr);
245 T binWidth = (xMax - xMin) / T(nBins);
260 for (
size_t i = 0;
i < nBins;
i++) {
282 for (
size_t ibin = 0; ibin < nBins; ibin++) {
285 Double_t xcenter = xMin + (ibin + 0.5) * binWidth;
286 Double_t error = 1. / TMath::Sqrt(entriesI);
288 fitter.AddPoint(&xcenter,
val, error);
291 A(npoints, 1) = xcenter;
292 A(npoints, 2) = xcenter * xcenter;
294 meanCOG += xcenter * entriesI;
295 rms2COG += xcenter * entriesI * xcenter;
316 fitter.GetParameters(par);
317 fitter.GetCovarianceMatrix(mat);
318 chi2 = fitter.GetChisquare() / Double_t(npoints);
320 if (TMath::Abs(par[1]) < kTol) {
323 if (TMath::Abs(par[2]) < kTol) {
326 param[1] = T(par[1] / (-2. * par[2]));
327 param[2] = T(1. / TMath::Sqrt(TMath::Abs(-2. * par[2])));
328 Double_t lnparam0 = par[0] + par[1] *
param[1] + par[2] *
param[1] *
param[1];
329 if (lnparam0 > 307) {
332 param[0] = TMath::Exp(lnparam0);
343 param[2] = TMath::Sqrt(TMath::Abs(meanCOG * meanCOG - rms2COG));
350 param[2] = binWidth / TMath::Sqrt(12);
360double fitGaus(
size_t nBins,
const T* arr,
const T xMin,
const T xMax, std::array<double, 3>&
param,
362 int minVal = 2,
bool applyMAD =
true)
364 double binW = double(xMax - xMin) / nBins,
s0 = 0,
s1 = 0, s2 = 0, s3 = 0, s4 = 0, sy0 = 0, sy1 = 0, sy2 = 0, syy = 0;
366 int bStart = 0, bEnd = (
int)nBins;
367 const float nSigmaMAD = 2.;
369 std::array<double, 3> madPar;
370 if (!
medmadGaus(nBins, arr, xMin, xMax, madPar)) {
373 bStart = std::max(bStart,
int((madPar[1] - nSigmaMAD * madPar[2] - xMin) / binW));
374 bEnd = std::min(bEnd, 1 +
int((madPar[1] + nSigmaMAD * madPar[2] - xMin) / binW));
376 float x = xMin + (bStart - 0.5) * binW;
377 for (
int i = bStart;
i < bEnd;
i++) {
381 throw std::runtime_error(
"Log-normal fit is possible only with non-negative data");
386 double y = std::log(
v), err2i =
v, err2iX = err2i, err2iY = err2i *
y;
394 sy1 += (err2iY *=
x);
395 sy2 += (err2iY *=
x);
401 auto recover = [&
param, binW, np,
s0,
s1, s2, sy0]() {
404 param[2] = np == 1 ? binW / std::sqrt(12) : std::sqrt(std::abs(
param[1] *
param[1] - s2 /
s0));
414 m33(1, 1) = m33(2, 0) = s2;
418 auto m33i = m33.Inverse(
res);
421 LOG(error) << np <<
" points collected, matrix inversion failed " << m33;
430 double chi2 =
v(0) *
v(0) *
s0 +
v(1) *
v(1) * s2 +
v(2) *
v(2) * s4 + syy +
431 2. * (
v(0) *
v(1) *
s1 +
v(0) *
v(2) * s2 +
v(1) *
v(2) * s3 -
v(0) * sy0 -
v(1) * sy1 -
v(2) * sy2);
433 param[2] = 1. / std::sqrt(-2. *
v(2));
435 if (std::isnan(
param[0]) || std::isnan(
param[1]) || std::isnan(
param[2])) {
442 j33(0, 0) =
param[0];
444 j33(0, 2) = j33(0, 1) *
param[1];
445 j33(1, 1) = -0.5 /
v(2);
446 j33(1, 2) = -
param[1] /
v(2);
447 j33(2, 2) =
param[2] * j33(1, 1);
448 *covMat = ROOT::Math::Similarity(j33, m33i);
450 return np > 3 ? chi2 / (np - 3.) : 0.;
478 double binWidth = (xMax - xMin) / (
double)nBins;
484 for (
size_t ibin = 0; ibin < nBins; ++ibin) {
485 double entriesI = (double)arr[ibin];
486 double xcenter = xMin + (ibin + 0.5) * binWidth;
488 mean += xcenter * entriesI;
489 rms2 += xcenter * entriesI * xcenter;
504 data.mStdDev = binWidth / std::sqrt(12.);
506 data.mStdDev = std::sqrt(std::abs(rms2 - mean * mean));
519template <
typename T,
typename R =
double>
525 auto n =
v.size() / 2;
526 nth_element(
v.begin(),
v.begin() +
n,
v.end());
528 if (!(
v.size() & 1)) {
529 auto max_it = max_element(
v.begin(),
v.begin() +
n);
530 med =
R{(*max_it + med) / 2.0};
543 LOG(error) <<
"Vector with values must have same size as vector for indices";
546 std::iota(
index.begin(),
index.end(),
static_cast<size_t>(0));
547 std::sort(
index.begin(),
index.end(), [&](
size_t a,
size_t b) { return values[a] < values[b]; });
572 int nPoints =
data.size();
573 std::vector<float>
w(2 * nPoints);
574 int nKeep = nPoints * fracKeep;
575 if (nKeep > nPoints) {
586 for (
int i = 0;
i < nPoints;
i++) {
591 w[
i + nPoints] = sum2;
593 double maxRMS = sum2 + 1e6;
595 int limI = nPoints - nKeep + 1;
596 for (
int i = 0;
i < limI;
i++) {
597 const int limJ =
i + nKeep - 1;
598 sum1 =
static_cast<double>(
w[limJ]) -
static_cast<double>(
i ?
w[
i - 1] : 0.);
599 sum2 =
static_cast<double>(
w[nPoints + limJ]) -
static_cast<double>(
i ?
w[nPoints +
i - 1] : 0.);
600 const double mean = sum1 / nKeep;
601 const double rms2 = sum2 / nKeep - mean * mean;
613 LOG(error) <<
"Rounding error: RMS = " <<
params[2] <<
" < 0";
630 LOG(error) <<
"Reordering not possible if number of elements in index container different from the data container";
633 std::vector<T> tmp(
data);
634 for (
size_t i = 0;
i <
data.size(); ++
i) {
649bool LTMUnbinnedSig(
const std::vector<T>&
data, std::vector<size_t>&
index, std::array<float, 7>&
params,
float fracKeepMin,
float sigTgt,
bool sorted =
false)
651 int nPoints =
data.size();
652 std::vector<double> wx(nPoints);
653 std::vector<double> wx2(nPoints);
665 for (
int i = 0;
i < nPoints;
i++) {
672 int keepMax = nPoints;
673 int keepMin = fracKeepMin * nPoints;
674 if (keepMin > keepMax) {
677 float sigTgt2 = sigTgt * sigTgt;
680 double maxRMS = wx2.back() + 1e6;
681 int keepN = (keepMax + keepMin) / 2;
686 int limI = nPoints - keepN + 1;
687 for (
int i = 0;
i < limI; ++
i) {
688 const int limJ =
i + keepN - 1;
689 sum1 = wx[limJ] - (
i ? wx[
i - 1] : 0.);
690 sum2 = wx2[limJ] - (
i ? wx2[
i - 1] : 0.);
691 const double mean = sum1 / keepN;
692 const double rms2 = sum2 / keepN - mean * mean;
702 if (maxRMS < sigTgt2) {
707 if (keepMin >= keepMax - 1) {
726 int i,
j, mid,
ir = np - 1, l = 0;
730 if (
ir == l + 1 && arr[
ir] < arr[l]) {
731 std::swap(arr[l], arr[
ir]);
735 int mid = (l +
ir) >> 1,
i = l + 1;
736 std::swap(arr[mid], arr[
i]);
737 if (arr[
i] > arr[
ir]) {
738 std::swap(arr[
i], arr[
ir]);
740 if (arr[l] > arr[
ir]) {
741 std::swap(arr[l], arr[
ir]);
743 if (arr[
i] > arr[l]) {
744 std::swap(arr[
i], arr[l]);
751 }
while (arr[
i] <
a);
754 }
while (arr[
j] >
a);
758 std::swap(arr[
i], arr[
j]);
784 for (
int i = np;
i--;) {
795template <
typename DataTimeType,
typename DataTime>
796std::optional<std::pair<size_t, size_t>>
findClosestIndices(
const std::vector<DataTimeType>& timestamps, DataTime timestamp)
798 if (timestamps.empty()) {
799 LOGP(warning,
"Timestamp vector is empty!");
803 if (timestamp <= timestamps.front()) {
804 return std::pair{0, 0};
805 }
else if (timestamp >= timestamps.back()) {
806 return std::pair{timestamps.size() - 1, timestamps.size() - 1};
809 const auto it = std::lower_bound(timestamps.begin(), timestamps.end(), timestamp);
810 const size_t idx = std::distance(timestamps.begin(), it);
811 const auto prevTimestamp = timestamps[idx - 1];
812 const auto nextTimestamp = timestamps[idx];
813 return std::pair{(idx - 1), idx};
845template <
typename DataTimeType,
typename DataType,
typename DataTime>
849 const size_t vecSize =
times.size();
852 if (!std::is_sorted(timeData.begin(), timeData.end())) {
853 LOGP(error,
"Input data is NOT sorted!");
857 if (timeData.empty()) {
858 LOGP(error,
"Input data is empty!");
863 const size_t timeDataSize = timeData.size();
865 LOGP(error,
"Input data has different sizes {}!={}", timeDataSize,
dataSize);
869 auto myThread = [&](
int iThread) {
872 for (
size_t i = iThread;
i < vecSize;
i += mNthreads) {
873 const double timeI =
times[
i];
876 const double timeStampLower = timeI - deltaMax;
877 const auto lower = std::lower_bound(timeData.begin(), timeData.end(), timeStampLower);
878 size_t idxStart = std::distance(timeData.begin(), lower);
881 const double timeStampUpper = timeI + deltaMax;
882 const auto upper = std::lower_bound(timeData.begin(), timeData.end(), timeStampUpper);
883 size_t idxEnd = std::distance(timeData.begin(), upper);
887 auto [idxLeft, idxRight] = *idxClosest;
888 const auto closestL = std::abs(timeData[idxLeft] - timeI);
889 const auto closestR = std::abs(timeData[idxRight] - timeI);
890 stats.closestDistanceL[
i] = closestL;
891 stats.closestDistanceR[
i] = closestR;
894 const size_t reqSize = idxEnd - idxStart;
895 if (reqSize < minPoints) {
897 idxStart = (idxRight > nClosestPoints) ? (idxRight - nClosestPoints) : 0;
898 idxEnd = std::min(
data.size(), idxRight + nClosestPoints);
899 constexpr float epsilon = 1e-6f;
900 double weightedSum = 0.0;
901 double weightTotal = 0.0;
902 for (
size_t j = idxStart;
j < idxEnd; ++
j) {
903 const double dist = std::abs(timeI - timeData[
j]);
904 const double weight = 1.0 / (dist + epsilon);
908 stats.median[
i] = (weightTotal > 0.) ? (weightedSum / weightTotal) : 0.0f;
911 stats.nPoints[
i] = reqSize;
913 if (idxStart >=
data.size()) {
914 stats.median[
i] =
data.back();
919 stats.median[
i] =
data[idxStart];
925 if (reqSize > window.capacity()) {
926 window.reserve(
static_cast<size_t>(reqSize * 1.5));
928 window.insert(window.end(),
data.begin() + idxStart,
data.begin() + idxEnd);
929 const size_t middle = window.size() / 2;
930 std::nth_element(window.begin(), window.begin() + middle, window.end());
931 stats.median[
i] = (window.size() % 2 == 1) ? window[middle] : ((window[middle - 1] + window[middle]) / 2.0);
934 const float mean = std::accumulate(window.begin(), window.end(), 0.0f) / window.size();
935 std::transform(window.begin(), window.end(), window.begin(), [mean](
const float val) { return val - mean; });
936 const float sqsum = std::inner_product(window.begin(), window.end(), window.begin(), 0.0f);
937 const float stdev = std::sqrt(sqsum / window.size());
938 stats.std[
i] = stdev;
944 std::vector<std::thread> threads(mNthreads);
945 for (
int i = 0;
i < mNthreads;
i++) {
946 threads[
i] = std::thread(myThread,
i);
949 for (
auto& th : threads) {
std::vector< unsigned long > times
float sum(float s, o2::dcs::DataPointValue v)
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat s1
GLuint GLuint GLfloat weight
GLboolean GLboolean GLboolean b
GLenum const GLfloat * params
GLenum GLsizei GLsizei GLint * values
GLfloat GLfloat GLfloat GLfloat v3
GLboolean GLboolean GLboolean GLboolean a
GLubyte GLubyte GLubyte GLubyte w
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat s0
bool medmadGaus(size_t nBins, const T *arr, const T xMin, const T xMax, std::array< double, 3 > ¶m)
void Reorder(std::vector< T > &data, const std::vector< size_t > &index)
TFitResultPtr fit(const size_t nBins, const T *arr, const T xMin, const T xMax, TF1 &func, std::string_view option="")
bool LTMUnbinnedSig(const std::vector< T > &data, std::vector< size_t > &index, std::array< float, 7 > ¶ms, float fracKeepMin, float sigTgt, bool sorted=false)
void SortData(std::vector< T > const &values, std::vector< size_t > &index)
RollingStats getRollingStatistics(const DataTimeType &timeData, const DataType &data, const DataTime ×, const double deltaMax, const int mNthreads, const size_t minPoints=4, const size_t nClosestPoints=4)
calculates the rolling statistics of the input data
StatisticsData getStatisticsData(const T *arr, const size_t nBins, const double xMin, const double xMax)
T MAD2Sigma(int np, T *y)
bool LTMUnbinned(const std::vector< T > &data, std::vector< size_t > &index, std::array< float, 7 > ¶ms, float fracKeep)
T selKthMin(int k, int np, T *arr)
std::optional< std::pair< size_t, size_t > > findClosestIndices(const std::vector< DataTimeType > ×tamps, DataTime timestamp)
R median(std::vector< T > v)
Double_t fitGaus(const size_t nBins, const T *arr, const T xMin, const T xMax, std::vector< T > ¶m)
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
Defining DataPointCompositeObject explicitly as copiable.
std::vector< float > closestDistanceL
distance of closest point to the left
std::vector< float > std
std of rolling data
std::vector< int > nPoints
number of points used for the calculation
std::vector< float > closestDistanceR
distance of closest point to the right
RollingStats(const int nValues)
std::vector< float > median
median of rolling data
ClassDefNV(RollingStats, 1)
double mCOG
calculated centre of gravity
double mStdDev
standard deviation
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"
o2::InteractionRecord ir(0, 0)