xiaoguangbin/AnalysisSystemForRadionuclide_qt_alg
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
main
AnalysisSystemForRadionucli.../GammaAnalyALG.h
wanglong fe493fb636 Initial commit on main branch
2024-06-04 15:25:02 +08:00

1074 lines
49 KiB
C++
Raw Permalink Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

#ifndef GAMMAANALYALG_H
#define GAMMAANALYALG_H
#include <jni.h>
#include "gamma_alg_global.h"
#pragma warning(disable:4190)
namespace AlgFunc
{
/** 加载解析原始谱文件
* @param phd：用于存储谱文件的各项信息
* @param fileContents: 为空时读取phd中filepath指定文件的内容进行分析，否则直接分析fileContents的内容
* @return bool 是否解析成功
*/
extern "C" bool LoadSpectrum(PHDFile* phd, QString fileContents = QString());
/** 读用户核素库
* @param sample_type：样品类型
* @param filename：为空则读默认用户核素库文件
* @return QStringList：用户核素列表
*/
extern "C" QStringList UserNuclide(QString sample_type, QString filename = "");
/** 核素识别
* @param phd：谱信息结构指针
* @param mapLines：用户核素库中核素的射线信息
* @param xmlPath:xml文件路径
* @return void
*/
extern "C" void NuclidesIdent(PHDFile* phd, QMap<QString, NuclideLines> &mapLines, QString xmlPath);
/** 交互分析界面中给某峰添加识别的核素时计算该核素的活度浓度等信息
* @param phd：谱信息结构指针
* @param lines：该核素射线信息
* @param nucName：核素名称
* @param vPeakIdx：与该核素相关的峰索引集合
* @return void
*/
extern "C" void CalcNuclideMDA(PHDFile* phd, NuclideLines &lines, const QString &nucName, std::vector<int> &vPeakIdx);
/** 执行质量检查
* @param phd：谱信息结构指针
* @return void
*/
extern "C" void RunQC(PHDFile *phd, QString strFPath);
/** 效率刻度数据和拟合方程更改后重新计算各峰的效率
* @param phd：谱信息结构指针
* @return void
*/
extern "C" void UpdateEfficiency(PHDFile *phd);
/** 根据选择的中心道，插入峰
* @param phd：谱信息结构指针
* @param spec：原始谱计数
* @param iC：新插入的峰的中心道
* @param minNA：最小净面积
* @return vLeft：新插入的峰的左边界
* @return vRight：新插入的峰的右边界
* @return QVector<int>：新插入的峰在峰列表中的索引
*/
//stdvec &vLeft, stdvec &vRight, PHDFile *phd, arma::rowvec spec, arma::vec iC, double minNA = 1.0
extern "C" QVector<int> insertpeak(stdvec &vLeft, stdvec &vRight, PHDFile *phd, arma::rowvec spec, arma::vec iC, double minNA = 1.0);
/** 峰拟合函数
* @param phd：谱信息结构指针
* @param Af：净面积需重新拟合的峰的索引
* @param Cf：中心道需重新拟合的峰的索引
* @param Ff：半高宽需重新拟合的峰的索引
* @return void
*/
extern "C" void fitPeakFull(PHDFile *phd, arma::uvec Af, arma::uvec Cf, arma::uvec Ff);
/** 能量转换为道
* @param e：能量值集合
* @param p：能量刻度系数
* @return stdvec：返回能量值对应的道的集合
*/
extern "C" stdvec energyToChannel(arma::rowvec e, arma::vec p);
/** 根据新的控制点信息重新计算基线数据
* @param bc：基线控制点信息
* @return void
*/
extern "C" void UpdateBaseControl(BaseControls& bc);
/** 峰曲线绘制相关函数，功能是取比1道更密集的点来计算各点上的峰计数以绘制峰曲线
* @param peak：峰信息
* @param regBase：该峰区域的基线数据
* @param u：道集合，道间距一般小于1，道间距越小，绘制出的曲线越光滑，但绘制效率越低
* @return stdvec：返回与道集合一一对应的峰计数
*/
extern "C" stdvec interp1(const PeakInfo &peak, arma::vec regBase, arma::vec u);
/** 拟合函数外部接口，根据道计算能量或根据能量计算分辨率或根据能量计算效率
* @param Chan：拟合函数的x值
* @param p：刻度系数
* @return stdvec：返回拟合函数的y值
*/
extern "C" stdvec calValuesOut(arma::rowvec Chan, arma::vec p);
/** 求导拟合函数
* @param Chan：道
* @param p：能量刻度系数
* @return double：dE
*/
extern "C" double calDerivaOut(double Chan, arma::vec p);
/** 求刻度数据的拟合系数
* @param cal：刻度类型
* @param funcId：拟合方程代号
* @param x：刻度数据x值
* @param y：刻度数据y值
* @param err：刻度数据不确定度
*/
extern "C" ParameterInfo calFitPara(CalibrationType cal, int funcId, stdvec x, stdvec y, stdvec err);
/** 在整个峰列表范围内计算所有峰的左右边界
* @param vPeaks：峰信息列表
* @param n：谱总道数
* @param c：峰列表中所有道的集合
* @param fwhmch：峰列表中所有半高宽的集合
* @param tlo：峰列表中所有tail的集合
* @param thi：峰列表中所有upperTail的集合
* @return void
*/
extern "C" void ComputePeakRange(std::vector<PeakInfo> &vPeaks, int n, arma::vec c, arma::vec fwhmch, arma::vec tlo, arma::vec thi);
/** 将Lc、Scac数据按事先约定好的格式写到文件中
* @param vData：数据集合
* @param type：数据标识
* @param filename：文件路径
* @return bool：是否写成功
*/
extern "C" bool WriteLcScac(const stdvec &vData, QString type, const QString &filename);
/** 将基线控制信息按事先约定好的格式写到文件中
* @param baseCtrl：基线控制信息
* @param filename：文件路径
* @return bool：是否写成功
*/
extern "C" bool WriteBaseInfo(BaseControls &baseCtrl, const QString &filename);
/** 根据方程代号获取方程描述信息（公式字符串）
* @param funcId：方程代号
* @return QString：方程描述
*/
extern "C" QString EquationDescription(int funcId);
/** 根据方程代号获取方程名称
* @param funcId：方程代号
* @return QString：方程名称
*/
extern "C" QString EquationName(int funcId);
/** 计算核素活度浓度
* @param phd：谱信息结构指针
* @param mainPeakIdx：用于计算核素活度浓度的峰的索引
* @param lambda：log(2.0)除以核素半衰期
* @param keyLineYield：主射线的分支比
* @param CCF：校正因子
* @return nucActMda：返回核素的活度浓度等计算结果
* @return bool：若activity大于mda，则返回true，否则返回false
*/
extern "C" bool CalculateMDCs(PHDFile *phd, NuclideActMda &nucActMda, int mainPeakIdx, double lambda, double keyLineYield, double CCF = 1.0);
/** 质量检查时计算 Ba140 和 Xe133 的 MDC 值
* @param phd：谱信息结构指针
* @param vMdcInfo：有三个元素，依次是“能量、分支比、半衰期”
* @param CCF：校正因子
* @return double：mdc的值
*/
double CalculateMDC(PHDFile *phd, stdvec vMdcInfo, double CCF = 1.0);
/** 从应用程序执行目录下读 SystemManager.xml 配置文件中的质量检查标准
* @param ener_Be7：Be7的能量
* @param vMdcInfo：有三个元素，依次是“能量、分支比、半衰期”
* @param system_type：谱类型：P 或 G
* @return qcItems：用于返回QC标准
* @return bool：信息读取是否成功
*/
bool ReadQCLimit(QMap<QString, QcCheckItem>& qcItems, double &ener_Be7, stdvec &vMdcInfo, QString system_type, QString strFPath);
/** 过滤核素能量小于 lowE 的射线
* @param lines：核素射线信息
* @param lowE：最小射线能量限制
* @return void
*/
void FilterNuclideLine(NuclideLines &lines, double lowE);
/** 插入新峰时计算新插入峰的左右边界
* @param pnr：新插入的峰在峰列表中的索引
* @param C：峰列表中所有峰的中心道的集合
* @param F_Ch：峰列表中所有峰的半高宽的集合
* @param rg_low：分析范围的左边界（单位：道）
* @param rg_high：分析范围的右边界（单位：道）
* @return mat：矩阵共两列，第一列为新插入峰的左边界，第二列为新插入峰的右边界
*/
arma::mat findPeakRange(arma::uvec pnr, arma::vec C, arma::vec F_Ch, double rg_low, double rg_high);
arma::rowvec findVPeakRange(int pNr, arma::vec C, arma::vec FC, arma::rowvec rg, double LD = 2.0, double RD = 2.5);
}
class GammaAnalyALG : public QObject
{
Q_OBJECT
public:
GammaAnalyALG();
GammaAnalyALG(SpecSetup setup);
GammaAnalyALG(PHDFile *phd);
~GammaAnalyALG();
/** 谱算法分析入口
* @param phd：谱信息结构指针
* @param mapLines：用户核素库中核素的射线信息
* @return bool：分析是否成功
*/
bool AnalyseSpectrum(PHDFile *phd, QMap<QString, NuclideLines> mapLines);
// 谱分析配置信息
void SetSpecSetup(SpecSetup setup);
// 设置刻度数据
void SetCalData(CalibrationType cal, stdvec x, stdvec y, stdvec err, QString calName = CalPHD); // 注：当为能量刻度时，x为道，y为能量
void SetCalData(QStringList names, G_EnergyBlock EnerKD, G_ResolutionBlock ResoKD, G_EfficiencyBlock EffiKD, TotaleffBlock TotEKD);
stdvec2 GetCalData(CalibrationType cal, QString calName = CalPHD);
// 设置/获取刻度拟合参数
void SetCalPara(CalibrationType cal, ParameterInfo para, QString calName = CalPHD);
void SetCalPara(QStringList paraNames, ParameterInfo EnerPara, ParameterInfo ResoPara, ParameterInfo EffiPara, ParameterInfo TotEPara);
ParameterInfo GetCalPara(CalibrationType cal, QString calName = CalPHD);
// 谱基本信息
void SetBaseInfo(const QVector<long> &spec, int begin_channel);
void SetBaseInfo(stdvec spec);
void SetBaseInfo(stdvec2 info);
stdvec2 getBaseInfo();
stdvec getBaseInfo(SpecBaseIdx idx);
// PAT分析结果存储表
void SetPAT(const std::vector<PeakInfo> &vPeaks);
std::vector<PeakInfo> getPAT();
stdvec getPatInfo(PeakIndex idx);
void setFilePath(QString strPath);
void setJniBl(JNIEnv* pEnv, jobject jobj, QString strUid);
signals:
void process_state(int);
private:
void Init();
bool Process(bool bUpdate = false, QChar dataType = QChar(), arma::vec certEne = arma::vec());
void callGammaProcess(int npro);
arma::field<arma::vec> FitCalPara(CalibrationType cal, arma::mat calData);
/* =====================================================寻峰 20======================================================== */
/* 寻峰入口函数
* @syntax [idx, C, NA, CNL, CNR, PSSfound] = peakSearch(sn, PSSs, type, PSSb)
* @param PSSs: peak search sensitivity treshold, default is setup variable PSS_low
* @param type: search full or stripped spectrum: "Full" or "Residual", default is "Full"
* @param ELim: peak search limits, default [ECutAnalysis_Low ECutAnalysis_High]
* @return idx: index of found peaks in PAT */
arma::uvec PeakSearch(double PSSs = qQNaN(), PeakSearchType PStype = search_full, arma::rowvec ELim = arma::rowvec());
arma::uvec PeakSearch(arma::vec &C, arma::vec &NA, arma::vec &CNL, arma::vec &CNR, arma::vec &PSSfound,
double PSSs = qQNaN(), PeakSearchType PStype = search_full, arma::rowvec ELim = arma::rowvec());
/* Returns the y values of a calibration that is valid for a named PAT of a given spectrum.
* @syntax [y, s] = calValues(sn, cal, pat, x)
* @param sn spectrum running number
* @param calibration type string ('energy', 'resolution', ...)
* @param pat PAT name string
* @param x values of independent variable at which calibration shall be evaluated
* @return y values of the dependent (calibration) variable at x
* @return s status 0: failure, 1: ok */
public:
arma::mat calValues(CalibrationType cal, arma::mat Chan);
private:
/* Returns the slope of a calibration curve of a given named PAT of a spectrum.
* function [dy, s] = calDerivative(sn, cal, pat, x) */
arma::mat calDerivative(CalibrationType cal, arma::mat Chan);
/* @syntax [s, msg, p, perr, dname, truenm] = calGetPATPara(sn, cal, pat)
* @param sn spectrum running number
* @param cal calibration type string
* @param pat pat name (compare getPATByName)
* @return s status 0: failed, 1: ok
* @return msg error message
* @return p calibration parameter
* @return perr parameter error estimate
* @return dname data name
* @return truenm calibration parameter name */
bool calGetPATPara(arma::colvec &p, CalibrationType cal);
/* @syntax [s, msg, out] = getCalParaGlob(sn, cal)
%% @syntax [s, msg, out] = getCalParaGlob(sn, cal, desc)
%% @param sn spectrum running number
%% @param cal calibration type string
%% @param desc descriptor string, optional<BR>
%% 'all' - default, return complete list<BR>
%% 'names' - return only column containing parameter names<BR>
%% 'datanames' - return only column containing names of related data points */
bool getCalParaGlob(arma::colvec &out, CalibrationType cal);
/* @syntax [s, msg, p, perr] = calDefaultPara(cal)
* @param cal calibration type string
* @return s status 0: unknown calibration, 1: ok
* @return msg error message
* @return p default calibration parameter vector
* @return perr default calibration parameter error vector */
bool calDefaultPara(arma::colvec &p, CalibrationType cal);
/* @syntax [s, msg, p, perr, dname, truenm] = getCalibrationPara(sn, cal, name)
* @param sn spectrum running number
* @param cal string specifying for which calibration the parameters shall be
* retrieved. One of: <BR> 'energy', 'resolution',
* 'efficiency', 'tot_efficiency', 'tail'
* @param name name of the calibration parameter to be returned
* @return s status 0: failure, 1: parameter ok
* @return p calibration parameter
* @return perr parameter error estimates
* @return dname name of related data points
* @return truenm true name of calibration parameter */
bool getCalibrationPara(arma::colvec &p, CalibrationType cal);
/* get analysis range
* @syntax rg = lGetAnalysisRange(ECutLow, ECutHigh, spec, E, res_ch)
* @param ECutLow lower energy cutoff
* @param ECutHigh upper energy cutoff
* @param spec counts in channels 1...NChan
* @param E energy of channels 1...NChan */
arma::rowvec lGetAnalysisRange(double ECutLow, double ECutHigh, arma::mat spec, arma::mat E, arma::mat res_ch);
/* @syntax pss = pksens(spec, res_c, LROI, RROI)
* @param spec counts in channels 1 to length(spec)
* @param res_c resolution in channel units in each ROI
* @param LROI left ROI channel
* @param RROI right ROI channel
* @return pss minimum peak search sensitivities in each ROI */
arma::colvec pksens(arma::rowvec spec, arma::mat res_c, arma::mat LROI, arma::mat RROI);
/* sdfilt - generalized second difference (peak search sensitivity) filter
* @syntax [f,k]=sdfilt(RC,CO)
* @syntax [f,k]=sdfilt(RC)
* @param RC resolution (full width at half maximum) in channel of peaks that shall be detected
* @param CO filter cutoff percentage (defaults to 1)
* @return f filter
* @return k index of filter mid point */
int sdfilt(arma::rowvec &f, double RC, double CO = 1);
/* @syntax [base, roiLo, roiHi] = specCutROIs(spec, C, FC)
* @param spec spectrum counts, starting with channel 1
* @param C peak centroids, in channels
* @param FC peak FWHMs, in channels
* @return base spectrum with interpolation in Regions of interests
* @return roiLo lower borders of regions of interest
* @return roiHi upper borders of regions of interest */
void specCutROIs(arma::rowvec &base, arma::rowvec &roiLo, arma::rowvec &roiHi, arma::rowvec spec, arma::vec C, arma::vec FC);
/* patAddPeaks(sn, DESC1, VAR1,...) creates a temporary pat for the sn'th
* spectrum, adding the peaks with the properties given by descriptor -
* value pairs. The new peak numbers are being returned
* Results are not valid until a call to patCommit and may be revoked by patRollback
*
* Valid descriptors are:
* DESC REQUIRED
* 'Source' yes
* 'Centroid' yes
* 'FWHM_Ch' no
* 'NetArea' yes
* 'AreaError' no
* 'Step' no
* 'Sensitivity' yes
* 'CentroidMethod' yes
* 'NetAreaMethod' yes
* 'FWHMMethod' yes* *if FWHM_Ch given, else 'C' will be applied
* 'NuclideSoft' no
* @syntax [nidx, oidx] = patAddPeaks(sn, D_1, V_1, D_2, V_2, ...)
* @param sn spectrum running number
* @param D_i descriptor of the i'th peak variable
* @param V_i value of the variable described by D_i
* @return nidx index of inserted peaks in temporary PAT
* @return oidx index of old peaks in new pat. oidx(pn) is the peak number
* of the peak whose initial peak number was pn */
bool patAddPeaks(arma::uvec Idx, arma::mat Datas, arma::mat& peaks, arma::uvec &nidx, arma::uvec &oidx);
/* function [sqi, uqi] = patBaseVar(sn, idx, pat, updr)
* @param sn spectrum running number for which update shall be performed
* @param idx index of peaks to be updated. May also be string 'all' (default)
* @param pat PAT to be modified.
* @param updr flag 1: update residual dependent parameters as well, 0: just others (default)
* @return sqi mean signed residual, empty if updr is 0
* @return uqi mean unsigned residual, empty if updr is 0 */
void patBaseVar(arma::colvec& sqi, arma::colvec& uqi, arma::uvec idx = arma::uvec(), int updr = 0);
/* write baseline dependent variables to PAT
*
* @syntax lSetPat(sn, idx, patin, patout, mbc, cba, lc, ld, sqi, uqi)
* @param sn spectrum running number
* @param idx index of modified peaks
* @param patin pat to be modified
* @param patout PAT where result goes to
* @param mbc mean back counts
* @param cba corrected background area
* @param lc currie's LCs
* @param ld currie's LDs
* @param sqi signed quality index
* @param uqi unsigned quality index */
void lSetPat(arma::uvec idx, arma::vec mbc, arma::vec cba, arma::vec lc,
arma::vec ld, arma::vec sqi = arma::vec(), arma::vec uqi = arma::vec());
/* tmpbase - create temporary baseline
*
* tmpbase creates the global variable 'Ac_baseline_tmp'
* that overrides the baseline. Baseline dependend PAT variables are also
* updated in a temporary PAT
*
* @syntax tmpbase(sn,cx,cy,cdy)
* @syntax tmpbase(sn,cx,cy,cdy,ar)
*
* @param sn spectrum running number for which the temporary baseline shall be created
* @param cx baseline controlpoint channel values
* @param cy baseline controlpoint height (count) values
* @param cdy baseline control point slope values, NaN for not-a-knot
* @param ar analysis range (lower and upper cutoff, optional) default read */
void tmpbase(arma::mat cx, arma::mat cy, arma::mat cdy, arma::mat ar = arma::mat(0u,0u));
/* =====================================================基线拟合 7======================================================== */
/* @syntax [cx, cy, cdy] = baseInit(sn, PSSBreak, PSSfwhm)
* @param sn spectrum number
* @param PSSBreak peak search sensitivity at which CPs are introduced
* by force left and right of ROI
* @param PSSfwhm peak search sensitivity at which the ROI width is +-5FWHM
* @return cx control point x values
* @return cy control point y values
* @return cdy control point slopes, NaN where continuos 2nd derivative */
void baseInit(double PSSBreak = 100, double PSSfwhm = 100);
/* @syntax [roiLo, roiHi, roiMask] = findROI(sn, fLo, idx, pat, fHi)
* @param sn spectrum number
* @param fLo multiplyer for FWHM on left side, optional, default 2
* @param idx index of peaks, optional, default 'all'
* @param pat name of pat to be used, optional, default 'Current'
* @param fHi multiplyer for FWHM on right side, optional, default fLo
* @return roiLo lower borders of regions of interest
* @return roiHi upper borders of regions of interest
* @return roiMask boolean mask of length of spectrum; zero outside ROIs one inside ROIs */
void findROI(arma::urowvec &l, arma::urowvec &h, arma::rowvec &roiMask,
arma::vec fLo = arma::vec("2"), arma::uvec idx = arma::uvec(), arma::vec fHi = arma::vec());
/* get data for baseline initialization
* @syntax [data, step, rg, roiMask, L, H] = lGetData(sn, PSSBreak, PSSfwhm)
* @param sn spectrum number
* @param PSSBreak sensitivity of big peaks, beneath which CPs are introduced
* @param PSSfwhm sensitivity of peaks for which ROI width is +-5FWHM
* @return data fitting data, stripped spectrum minus step
* @return step cumulative step height
* @return rg analysis range
* @return roiMask boolean flag of invalid data points (ROIs)
* @return L left borders of Big ROIs
* @return H right borders of Big ROIs */
void lGetData(arma::rowvec &data, arma::rowvec &step, arma::rowvec &rg, arma::rowvec &roiMask,
arma::uvec &L, arma::uvec &H, double PSSBreak, double PSSfwhm);
/* restricted polynomial fit
*
* fits the free coefficients of a polynomial to data x,y
* The data points may be weighted (big weight for reliable data points)
*
* @syntax [p, r] = restpolyfit(x, y, p0, w)
* @param x data points of ordinate
* @param y data points at abscissa
* @param p0 polynomial coefficients, NaN entries for coefficients to be fitted
* @param w weights, optional, default ones(size(x)) */
arma::mat restpolyfit(arma::vec &r, arma::vec x, arma::vec y, arma::vec p0, arma::vec w = arma::vec());
/* get norm of least-squares fit polynomial to data
*
* Fits a first or third order polynomial to the used data points in an
* interval; the function value at the upper border can be fixed.
*
* The returned characteristics of the fit are:
* nrm - the square sum of the residual
* lbda - the mean of the square root of the fitted data
*
* @syntax [nrm,lbda,yc,dyc,xl,yh,dyh] = lNorm(yy, mask, xc, xh, yh, np)
* @param yy whole data
* @param mask boolean mask of data to be ignored
* @param xc midpoint of fitting interval
* @param xh upper end of fitting interval
* @param yh value at upper end, NaN if free parameter
* @param np number of fitting parameters; valid values are 2 (linear) or
* 4 (cubic)
* @return n square sum of residual
* @return m mean square root of data
* @return yc value of fit at center
* @return dyc derivative of fit at center
* @return xl lower interval border
* @return yh value of fit at upper border
* @return dyh derivative of fit at upper border */
void lNorm(arma::rowvec yy, arma::rowvec mask, double xc, double xh, double yh, double np,
double &nrm, double &lbda, double &yc, double &dyc, double &xl, double &yho, double &dyh);
/* selection of next control point
*
* [xn, yn, dyn, hit, y1, dy1] = lSelectNext(x, y, data, xb, imask, np)
* @param x starting point x value
* @param y starting point y value, NaN for calculated
* @param data data to be fitted
* @param xb ultimate x endpoint
* @param imask mask of data to be ignored
* @return xn next control point x value
* @return yn next control point y value
* @return dyn next control point slope
* @return hit flag if fitting interval includes xb
* @return y1 current control point y value
* @return dy1 current control point slope */
void lSelectNext(double x, double y, arma::rowvec data, double xb, arma::rowvec imask,
double np, double &xn, double &yn, double &dyn, bool &hit, double &y1, double &dy1);
/* fit polynomial in predefined range, return value and derivative at
* evaluation points.
*
* @syntax [yv, dyv] = lFitRange(x0, y0, x1, xv, data, mask, npara)
* @param x0 interval end at which values are prescibed
* @param y0 value at x0
* @param dy0 derivative at x0
* @param x1 other interval end
* @param xv points inside interval where fit shall be evaluated
* @param data data to be fitted
* @param mask mask of data to be ignored
* @param npara number of fitting parameters (order of polynomial + 1)
* valid values are 2 (linear) and 4 (cubic)
* @return yv values of fit at xv
* @return dyv values of derivative at xv */
void lFitRange(double x0, double y0, double x1, arma::rowvec xv, arma::rowvec data, arma::rowvec mask,
int npara, arma::rowvec &yv, arma::rowvec &dyv);
void lInitBase(const arma::rowvec &data, const arma::rowvec &step, const arma::rowvec &rg, const arma::rowvec &roiMask,
const arma::uvec &L, const arma::uvec &H, double N1, double N2, double N3, arma::rowvec &cx,
arma::rowvec &cy, arma::rowvec &cdy);
/* =====================================================基线优化 4======================================================== */
/*baseImprove - improve baseline
*
* improves baseline by fitting control points to the residual. Additional
* peaks in the residual are also searched and substracted from the
* spectrum before fitting the baseline. These peaks are not written to PAT.
* The improved baseline is written to temporary baseline.
*
* @syntax baseImprove(sn, pss_r)
* @param sn spectrum running number
* @param pss_r peak search sensitivity for residual peak search */
void baseImprove(double pss_r);
/* fitPeakFull - fit peaks above a full baseline
*
* fit peak parameters to match spectrum counts by a spline baseline plus
* gaussian peaks, possibly distorted by tailing
* Initial values for fitting are taken from PAT, results are written to
* temporary PAT, but also are returned
*
* @syntax [s, Ao, Co, Fo, X2] = fitPeakFull(sn, Af, Cf, Ff, SRf)
* @syntax [s, Ao, Co, Fo, X2] = fitPeakFull(sn, Af, Cf, Ff, SRf, optargs)
* @param sn running number of spectrum to be fitted
* @param Af index of peaks for which Area shall be free
* @param Cf index of peaks for which Centroid shall be free
* @param Ff index of peaks for which FWHM shall be free
* @param SRf index of peaks to which a step ratio shall be fitted
* @param optargs the required arguments may be followed by a list of optional
* arguments. Each set of optional arguments is started with an option
* string and followed by the arguments the option takes. Option strings and
* their arguments are:
* 'BaseFree' one argument: Cidx - also fit baseline control points indexed by
* Cidx
* 'Verbose' no arguments, show waitbar while fitting
* @return s status of success
* @return Ao fitted net areas
* @return Cf fitted centroids
* @return Ff fitted FHWM
* @return X2 mean square residual over analysis range */
void fitPeakFull(arma::uvec Af, arma::uvec Cf, arma::uvec Ff, arma::uvec SRf,
arma::uvec CPf = arma::uvec(), bool bVerbose = false);
/* @syntax patSetPeaks(sn, idx, D_1, V_1, D_2, V_2, ...)
* @param sn spectrum running number
* @param idx index of peaks for which the specified variables shall be set,
* may be string 'all'
* @param D_i descriptor of a variable that shall be changed, see abstract
* @param V_i new value of the variable described by D_i */
void patSetPeaks(arma::uvec idx, QStringList Descs, arma::field<arma::vec> Values, bool bAll = false);
/* fitBase - fit baseline to residual
* fit the Y values of baseline control points of a certain spectrum
* to the smoothened stripped spectrum using weights to punish peak regions.
* Results are written to temporary baseline and need to be accepted.
* The fitting may be restricted to a set of indexed control points
*
* @syntax fitBase(sn)
* @syntax fitBase(sn, cidx)
* @syntax fitBase(sn, cidx, sidx)
* @param sn running number of spectrum to be fitted
* @param cidx fit only the control points with index cidx
* @param sidx fit slope for control points with index sidx */
void FitBase(arma::uvec cidx = arma::uvec(), arma::uvec sidx = arma::uvec());
/* =====================================================谱基线结构 21======================================================== */
/* @param items: one or more DESCRIPTOR QString
* DESCRIPTOR VARIABLE
* "Spectrum" channel contents
* "SpectrumLong" channel contents, padded with NaNs to length of longest spectrum
* "BaseLine" fitted baseline
* "AnalysisRange" Analysis Range (first and last valid channel)
* "AnalysisRangeSave" Save Analysis range, i.e. if empty will be [1,NCh]
* "XControl" baseline control points x-data
* "YControl" baseline control points y-data
* "YSlope" slope at breakpoints, NaN at normal control points
* "Residual" residual (Spectrum-BaseLine-Peaks)
* "Stripped" stripped spectrum (gaussian peaks substracted)
* "Steps" non-spline baseline components
* "ROISmooth" stripped spectrum minus Steps, with smoothened ROIs
* "SmoothStripped" smoothened stripped spectrum
* "Approximation" spectrum approximation (baseline + peaks)
* "ChiSquare" mean weighted residual and number of channels with zero counts ( [ChiSquare, nzero] )
* "ArtXControl" artificial baseline control points x-data
* "ArtYControl" artificial baseline control points y-data
* "ArtYSlope" slope of artificial baseline control points
* "ArtificialBase" Artificial baseline (if not empty)
* "BaseFittingWeights" Weight function used in baseline fitting
* "PSS" value of peak search sensitivity in each channel
* @return f: the results of DESCRIPTORs in DESCRIPTOR order */
arma::field<arma::rowvec> dmspec(QStringList items);
/* @syntax [d, nzero] = msrd(Y, F)
* @param Y data values
* @param F approximation values
* @return d mean relative residual
* @return nzero number of values in Y that were excluded from calculation */
void msrd(arma::mat &d, int &nzero, arma::mat Y, arma::mat F);
/* smoothen data using savitzky-golay filter proportional to fwhm
*
* data is smoothend piecewise with a savitzky-golay filter with framewidth
* that is about 3*FWHM. Filter changes wherever the framewidth jumps to the
* next odd number.
*
* The data is assumed to start with channel 1.
*
* @syntax s = lSmoothen(sn, y)
* @param sn spectrum running number
* @param y data to be smoothened
* @return s smoothened data */
arma::mat lSmoothen(arma::rowvec y);
/* savitzgy-golay numerical approximation of derivative
*
* In the Savitzky-Golay method, the derivative of equally-spaced data is
* approximated by the derivative of the least-squares fitted polynomial
* through the data points, evaluated at the center of the fitting interval.
*
* Implemented orders are 0, 2, 3
* Note that sgderive(y, np, 0) yields np-point smoothened data.
*
* @syntax s = sgderive(y, np, order)
* @param y data
* @param np frame width
* @param order degree of derivation */
arma::mat sgderive(arma::mat y, int np, int order);
/* savitzky golay filter
*
* Cubic (quartic) savitzky golay filter with arbitrary frameWidth
*
* @syntax [f, k] = sgfilter(frameWidth, deriveOrder)
* @param frameWidth filter length (odd)
* @param deriveOrder order of derivation provided by the filter
* @return f filter (row vector)
* @return k filter half width (frameWidth-1)/2 */
void sgfilter(arma::rowvec &f, int &k, int frameWidth, int deriveOrder);
/* get smooth data inside regions of interest by using data from outside
*
* DOC
*
* The regions of interest are defined as the range from
* ROI_LOW = max(1, floor( c - f*fwhm )) to
* ROI_HIGH = min(n, ceil( c + f*fwhm ))
*
*
* @syntax s = lSmoothROIs(y, c, w)
* @param y data to be smoothened
* @param c peak centroids
* @param w peak widths
* @return s input data with smoothened ROIs */
arma::mat lSmoothROIs(arma::rowvec y, arma::mat c, arma::mat w);
/* =====================================================PAT结构 3======================================================== */
/* evaluate named calibration of a spectrum
*
* Returns the y values of a named calibration of a certain type.
*
* @syntax [y, s] = calValuesByName(sn, cal, name, x)
* @param sn spectrum running number
* @param cal calibration type string
* @param name name of the parameter in the list of parameters of that
* calibration and spectrum
* @param x values of independent variable at which calibration shall be
* evaluated
* @return values of the dependent (calibration) variable at x
* @return s status 0: failure, 1: ok */
arma::mat calValuesByName(CalibrationType cal, arma::mat x);
/* slope of a named calibration parameter
*
* Returns the slope of a calibration described by a named calibration parameter
*
* @syntax [dy, s] = calDerivativeByName(sn, cal, name, x)
* @param sn spectrum running number
* @param cal calibration type string
* @param name name of the calibration parameter (compare getCalibrationPara)
* @param x values at which slope shall be returned
* @return dy slope at x
* @return s status 0: some problem occurred, 1: ok */
arma::mat calDerivativeByName(CalibrationType cal, arma::mat x);
/* peak analysis table data managment
* dmps provides access to variables stored in the Peak Analysis Table (PAT) or
* some basic peak variables that can be derived from them */
void dmps();
arma::vec dmps(PeakIndex idx);
arma::field<arma::vec> dmps(QStringList items);
arma::field<arma::vec> dmps(QList<PeakIndex> list_idx);
/* function calUpdate(sn, E1, R, E2)
* calUpdate - update internal calibration parameters
*
* calUpdate supervises the calibration update.
* It gets the library lines for peaks to be used for updates, calls the
* calibration peak search and the calibration update procedures
*
* Calibration lines are read from the file indicated by the setup variable
* 'CalibrationLinesFile'. For C-type spectra, lines in the certificate block
* from message are used instead.
*
* @syntax calUpdate(sn, E1, R, E2)
* @param sn running number of spectrum
* @param E1 boolean: perform first energy calibration update or not
* @param R boolean: perform resolution calibration update or nor
* @param E2 boolean: perform second energy calibration update or not */
public:
void calUpdate(QChar dataType, arma::vec certEne, bool E1, bool R, bool E2);
private:
/* calPeakSearch - search calibration peaks
* calPeakSearch selects the biggest singlet peaks from the internal PAT with
* energies near library lines.
* If the internal PAT is empty, a peak search is performed first.
*
* @syntax [IDX, EF] = calPeakSearch(SN, ELib)
* @param SN running number of examined spectrum
* @param ELib all relevant library energy lines
* @return IDX index of found peaks in internal PAT
* @return EF library energies of the found peaks */
arma::uvec calPeakSearch(arma::vec &EF, arma::vec ELib);
/* peakFindMultiplets(C, NA, FC, RW) finds the multiplets in the list of
* peaks with centroids C, net areas NA and channel resolutions FC.
* RW gives the ROI width (curr: C +- RW*FC)
*
* [Mult, CL, CR] = peakFindMultiplets(C, NA, FC, RW) also gives left and
* right ROI channels
*
* @syntax [Mult, CL, CR] = peakFindMultiplets(C, NA, FC, RW)
* @param C peak centroids
* @param NA peak net areas
* @param FC peak fwhms in channel units
* @param RW ROI widths in fwhm units
* @return Mult multiplet indicators (0: singlet, 1: multiplet)
* @return CL left ROI channels
* @return CR right ROI channels */
arma::vec peakFindMultiplets(arma::vec C, arma::vec NA, arma::vec FC, double RW);
/*calEnergyUpdate1 - perform first energy calibration update
*
* calEnergyUpdate1 fits a curve to centroids from the PAT versus library
* energies. If that fitting was successful, an update of the internal
* energy calibration coefficients and the calibration data points is performed
* and the flag EneCalUpdate1 is raised.
*
* @syntax calEnergyUpdate1(SN, IDX, ELib)
* @param SN running number of spectrum, for which an update shall be performed
* @param IDX index of peaks used for fitting
* @param ELib library energies of those peaks */
void calEnergyUpdate1(arma::uvec idx, arma::vec ELib);
/* calResUpdate - perform resolution calibration update
*
* calResUpdate tries to improve the resolution calibration for a given
* spectrum. First a good estimate of the FWHM is determined for the data peaks
* using quickfit with free Net Area and FWHM. Then a curve is fitted to the
* peak energies versus FWHM. If the fitting was successful, the internal
* resolution calibration parameters and data points are updated.
*
* @syntax calResUpdate(SN, IDX)
* @param SN running number of spectrum to be updated
* @param IDX index of peaks in internal PAT to be used for fitting */
void calResUpdate(arma::uvec idx);
/* calEnergyUpdate2 - perform second energy calibration update
*
* calEnergyUpdate2 is intended to improve the energy calibration after the
* resolution of certain peaks has been improved. It first improves the
* centroid of used peaks by quickfit with free centroid and net area. Then
* a curve is fitted to the resulting centroids versus library energies.
* If that fitting was successful, the energy calibration parameters are
* updated and the flag EneCalUpdate2 is set
*
* @syntax calEnergyUpdate2(SN, IDX, ELib)
* @param SN running number of spectrum to be updated
* @param IDX index of peaks in current PAT to be used for fitting
* @param ELib library energies of those peaks */
void calEnergyUpdate2(arma::uvec idx, arma::vec ELib);
/* check validity of calibration data points
*
* @syntax [s, msg, data] = lCalDataCheck(data)
* @syntax [s, msg, data] = lCalDataCheck(x, y, err)
* @param data calibration data point matrix
* @param x data points x data
* @param y data points y data
* @param err data points error estimates
* @return s status 0: invalid, 1:valid
* @return msg error message
* @return data calibration data point matrix */
bool lCalDataCheck(arma::mat &data, arma::mat x, arma::mat y, arma::mat err = arma::mat(0u, 0u));
/* check parameter error estimates
*
* @syntax [s, msg, perr] = lCalCheckPerr(p, perr)
* @param p calibration parameter vector
* @param perr parameter error estimates
* @return s status, 0: invalid, 1: valid
* @return msg error message
* @return perr transformed parameter error estimates */
bool lCalCheckPerr(arma::mat p, arma::mat& perr);
/* add new calibration data to list
*
* adds a set of calibration data points to the list of data sets.
*
* @syntax s = setCalibrationData(sn, cal, src, data)
* @syntax s = setCalibrationData(sn, cal, src, x, y, err)
* @param sn spectrum running number
* @param cal string specifying for which calibration the parameters shall be
* set. One of: <BR> 'energy', 'resolution',
* 'efficiency', 'tot_efficiency', 'tail', 'tail_alpha', 'tail_right',
* 'tail_right_alpha', 'step_ratio'
* @param src source of the calibration data, must not contain spaces
* @param data nx3 matrix containing the independent variable as first row,
* the dependent variable as second row and the error estimate as third row
* @param x independent variable data
* @param y dependent variable data
* @param err error estimate
* @return s status 0: parameter not set, 1: parameter set
* @return msg error message
* @return name name of the calibration data set */
//bool setCalibrationData(QString src, arma::mat x, arma::mat y, arma::mat err = arma::mat(0u, 0u));
/* add new calibration parameter to list
*
* Adds a new calibration parameter to the list of calibration parameters
* of type cal for the snth spectrum.
*
* The special parameter [1] (interpolation) will be converted to the
* parameter [1, x1, y1, x2, y2, ...] .
*
* @syntax [s, name] = setCalibrationPara(sn, cal, src, p, perr, data)
* @param sn spectrum running number
* @param cal string specifying for which calibration the parameters shall be
* set. One of: <BR> 'energy', 'resolution',
* 'efficiency', 'tot_efficiency', 'tail', 'tail_alpha', 'tail_right',
* 'tail_right_alpha', 'step_ratio'
* @param src source of the calibration parameters
* @param p calibration parameter vector
* @param perr vector of error estimates, may be empty
* @param data string, name of the calibration data points used; or Nx3 matrix
* of calibration data points [x,y,err]
* @return s status 0: parameter not set, 1: parameter set
* @return name name of the parameter in the list */
//bool setCalibrationPara(QString src, arma::mat p, arma::mat perr);
/* sets current calibration to named parameter
*
* Links the current (true current and temporary) PAT calibration of type cal
* to a named calibration parameter.
*
* @syntax [s, msg] = setCurrentCalibration(sn, cal, name)
* @param sn spectrum running number
* @param cal calibration type string
* @param name name of calibration parameter
* @return s status 0: failed, 1: ok
* @return msg error string */
//bool setCurrentCalibration(CalibrationType cal, QString name);
/* fitPeakQuick - fit peaks using linear baseline
*
* fits a linear baseline plus gaussian peaks to spectrum data. The fitting
* interval is defined as the range of peaks that are used for fitting (any
* peak parameter is free) plus a factor times the fwhm of the outermost
* peaks.
* Results are written to temporary PAT.
*
* @syntax [s, Ao, Co, Fo, X2] = fitPeakQuick(sn, Af, Cf, Ff)
* @syntax [s, Ao, Co, Fo, X2] = fitPeakQuick(sn, Af, Cf, Ff, vb)
* @syntax [s, Ao, Co, Fo, X2] = fitPeakQuick(sn, Af, Cf, Ff, vb, FitWidth)
* @param sn spectrum running number
* @param Af index of peaks for which net area shall be fitting parameter
* @param Cf index of peaks for which centroid shall be fitting parameter
* @param Ff index of peaks for which fwhm shall be fitting parameter
* @param vb optional parameter, if called with this parameter, fitting is
* called verbose (i.e. show waitbar)
* @return s status of success (0 if failed)
* @return Ao fitted net areas (only free parameters)
* @return Co fitted centroids (only free parameters)
* @return Fo fitted fwhm (only free parameters)
* @return X2 mean square relative residual over the fitting interval */
bool fitPeakQuick(arma::vec &Ao, arma::vec &Co, arma::vec &Fo, double &X2,
arma::uvec Af, arma::uvec Cf, arma::uvec Ff,
QString vb = QString(""), double FitWidth = 2.5);
public:
arma::rowvec GetFwhmcAll();
stdvec calculateSCAC(arma::rowvec &spec, arma::rowvec &baseline, arma::rowvec &FwhmC);
stdvec calculateLC(arma::rowvec &baseline, arma::rowvec &FwhmC, double RiskLevelK = 4.264890);
private:
/*function [r, l, h] = lGetPlotRanges(n, c, fwhm, tlo, thi, bhi, rdelta)
get plotting ranges for peaks */
arma::umat lGetPlotRanges(int n, const arma::vec &c, const arma::vec &fwhmch, const arma::vec &tlo, const arma::vec &thi);
/* =================================================== 峰插入、删除、修改操作 =========================================== */
/*insertpeak - insert peak in temporary PAT
*
* insertpeak inserts peaks with centroids iC to the PAT of the sn'th specrum.
* Net areas are estimated by summation method, fwhms are taken from calibration.
* Results are written to temporary PAT and need to be committed or rolled back.
*
* @syntax [idx, gidx, oidx]=insertpeak(sn,iC)
* @syntax [idx, gidx, oidx]=insertpeak(sn,iC,src)
* @syntax [idx, gidx, oidx]=insertpeak(sn,iC,src, minNA)
* @param sn spectrum running number
* @param iC centroids of peaks to be inserted
* @param src source string, optional, default 'Manual', valid values:
* 'Manual' manual insertion
* 'H0Test' inserted to test library line
* 'MultFit' inserted by multiplet fitting routine
* 'SumPeakModel' inserted by summation peak model
* @param minNA minimum net area assigned, default 1
* @return idx index of inserted peaks in PAT
* @return gidx index of successfully inserted peaks in iC
* @return oidx index of old peaks. oidx(pn) is the new peak number of the peak
* with the initial peak number pn.*/
arma::uvec insertpeak(arma::vec iC, PeakInsertMethod src = Ins_Manual, double minNA = 1.0);
/*findPeakRange - find channel range of a multiplet
* rg = findPeakRange(sn, pnr) finds the range of the multiplet the pnr'th
* peak of spectrum sn belogs to
*
* @syntax [rg, iL, iH] = findPeakRange(sn, pnr)
* @syntax [rg, iL, iH] = findPeakRange(sn, pnr, wL, wR)
* @syntax [rg, iL, iH] = findPeakRange(sn, pnr, pat)
* @syntax [rg, iL, iH] = findPeakRange(sn, pnr, pat, wL, wR)
* @param sn spectrum running number
* @param pnr list of peak numbers, with length n
* @param pat pat name, optional, default 'Current'
* @param wL left ROI width in units of FWHM (optional)
* @param wR right ROI width in units of FWHM (optional)
* @return rg nx2 matrix of multiplet ranges; rg(1, i) is the left border
* channel of the multiplet that peak number pnr(i) belongs to, rg(2, i) the
* right border channel
* @return iL nx1 vector, index of lowest peak in multiplet
* @return iH nx1 vector, index of highest peak in multiplet*/
arma::mat findPeakRange(arma::uvec pnr, arma::vec C, arma::vec F_Ch, double rg_low, double rg_high);
/* function [l, h, i, j] = findVPeakRange(pNr, C, FC, rg, LD, RD)
* findVPeakRange - find channel range of a multiplet
*
* finds the low border l and the high border h (channels) of the multiplet,
* to which the peak with index pNr belongs from centroid and fwhm information.
*
* @syntax [l, h, i, j] = findVPeakRange(pNr, C, FC, rg)
* @param pNr peak number
* @param C peak centroids
* @param FC peak fwhms in channel units
* @param rg analysis range [low, high]
* @param LD left width of ROI in units of FWHM
* @param RD right width of ROI in units of FWHM
* @return l leftmost channel of multiplet containing the peak
* @return h rightmost channel of multiplet
* @return i index of leftmost peak in multiplet
* @return j index of rightmost peak in multiplet */
arma::rowvec findVPeakRange(int pNr, arma::vec C, arma::vec FC, arma::rowvec rg, double LD = 2.0, double RD = 2.5);
/* peakfitdlg - dialog for aatami peak fitting
*
* creates a peak fitting dialog for a given energy range and returns the
* the selected exit pushbutton name
* but=peakfitdlg('create',ELim) opens a dialog with a list of all peaks
* in the range ELim(1) to ELim(2) and allows fitting of these peaks;
* 'but' specifies, which branch has been executed:
* 'Fit' peaks fitted using spline baseline
* 'QuickFit' peaks fitted using linear baseline
* 'FitBase' baseline fitted to stripped spectrum
* 'Cancel' no actions are performed
*
* The function also handles the uicontrol callbacks of that GUI.
*
* @syntax but = peakfitdlg(action)
* @param action pushbutton name for which callback shall be executed
* @syntax but = peakfitdlg('create', ELim)
* @param ELim energy limits for which dialog shall be opened */
void peakfitdlg(QString action, arma::mat ELim);
public:
SpecSetup specSetup;
int m_nChans;
PAT pat;
arma::field<arma::rowvec> baseInfo;
// 存储 PHD 文件刻度数据点
QMap<QString, arma::mat> Cal_Ener_Data;
QMap<QString, arma::mat> Cal_Reso_Data;
arma::mat Cal_Effi_Data;
arma::mat Cal_Tot_Effi_Data;
private:
// 存储 PHD 文件刻度拟合参数
QMap<QString, arma::field<arma::vec> > Acal_energy_para;
QMap<QString, arma::field<arma::vec> > Acal_resolution_para;
arma::field<arma::vec> Acal_efficiency_para;
arma::field<arma::vec> Acal_tot_efficiency_para;
arma::vec Acal_tail_para;
arma::vec Acal_tail_alpha_para;
arma::vec Acal_tail_right_para;
arma::vec Acal_tail_right_alpha_para;
arma::vec Acal_step_ratio_para;
QString m_curEner; // 当前使用的能量刻度名称
QString m_curReso; // 当前使用的分辨率刻度名称
QString m_strFilePath;
QString m_strUserId;
// 存储刻度更新时的刻度数据点和对应的拟合参数
arma::field<arma::mat> Cal_update_data;
arma::field<arma::vec> Cal_update_para;
JNIEnv* m_pJniEnv;
jobject m_pSendObj;
//QMap2Data m_CalData; // 混合存储所有刻度数据点
//QMap2Para m_CalPara; // 混合存储所有刻度拟合参数
//QMap<CalibrationType, QString> m_mapCurCal; // 当前刻度
};
#endif // GAMMAANALYALG_H
Reference in New Issue View Git Blame Copy Permalink