BmnWeightedLine.h

Go to the documentation of this file.

/*
    BM@N alignment routine
    BM@N experiment at NICA complex, JINR, 2026
 
    Department: Math & Soft Group of HEP lab
    Author:     Claude (Opus 5)
    Supervisor: Igor Polev, polev@jinr.int
 
    BmnWeightedLine header-only implementation
    Weighted least squares straight line in one projection
*/
 
#ifndef BMNWEIGHTEDLINE_H
#define BMNWEIGHTEDLINE_H
 
#include <RtypesCore.h>
#include <TMath.h>
#include <algorithm>
#include <cmath>
 
// v = At0() + Slope() * z. Raw sums lose digits when the points lie far from z = 0 compared with
// their z spread; for the BM@N detectors they keep 12 digits or more.
// A point of zero weight adds nothing, but still counts in n.
class BmnWeightedLine
{
  private:
    Double_t s{0.0}, sz{0.0}, szz{0.0}, sv{0.0}, szv{0.0}, svv{0.0};
    Int_t n{0};
    // det value is heavily reused -- it needs to be cached
    mutable Double_t det{TMath::QuietNaN()};
 
  public:
    void Add(Double_t z, Double_t v, Double_t w = 1.0) noexcept
    {
        s += w;
        sz += w * z;
        szz += w * z * z;
        sv += w * v;
        szv += w * z * v;
        svv += w * v * v;
        n++;
        det = TMath::QuietNaN();
    }
    // the inverse of Add(): the line through the other points
    void Remove(Double_t z, Double_t v, Double_t w = 1.0) noexcept
    {
        s -= w;
        sz -= w * z;
        szz -= w * z * z;
        sv -= w * v;
        szv -= w * z * v;
        svv -= w * v * v;
        n--;
        det = TMath::QuietNaN();
    }
    Int_t N() const noexcept { return n; }
    Double_t Mean() const noexcept { return sv / s; }
    Double_t At(Double_t z) const noexcept { return At0() + Slope() * z; }
    Double_t At0() const noexcept { return (szz * sv - sz * szv) / Det(); }
    Double_t Slope() const noexcept { return (s * szv - sz * sv) / Det(); }
    Double_t Chi2() const noexcept { return svv - At0() * sv - Slope() * szv; }
    Double_t Rms() const noexcept { return std::sqrt(std::max(svv / s - Mean() * Mean(), 0.0)); }
    // Var(At0) + 2 z Cov(At0, Slope) + z^2 Var(Slope), from the weights alone
    Double_t Var(Double_t z) const noexcept { return (szz - 2.0 * z * sz + z * z * s) / Det(); }
    Double_t SlopeVar() const noexcept { return s / Det(); }
    Double_t Cov() const noexcept { return -sz / Det(); }   // Cov(At0, Slope)
    // cached functions used in hot code or heavily reused
    Double_t Det() const noexcept
    {
        if (TMath::IsNaN(det))
            det = s * szz - sz * sz;
        return det;
    }
};
 
#endif   // BMNWEIGHTEDLINE_H