AnalysisSystemForRadionucli.../GammaAnalyALG.cpp

6511 lines
210 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.

#include "GammaAnalyALG.h"
#include <QtMath>
#include <QFile>
#include <QCoreApplication>
#include <QTextStream>
#include <QStringList>
#include "fitfunc.h"
#include "genenalfunc.h"
#include "matlab_func.h"
#include "IndependentAlg.h"
#include <QDateTime>
#include <QDomDocument>
#include <QJsonDocument>
#include <QJsonObject>
#include <QJsonArray>
//#define gNaN qSNaN()
using namespace arma;
void PrintMat22(QString name, mat& dataMat);
void PrintColvec(QString name, colvec& dataMat);
void PrintMat2fvec(QString name, field<vec>& vecc);
namespace AlgFunc
{
void jsonArrToStdvec(QJsonObject jObj, QString strKey, stdvec& val)
{
QJsonArray usedArr = jObj.value(strKey).toArray();
val.clear();
for (int i = 0; i < usedArr.size(); i++)
{
val.push_back(usedArr.at(i).toDouble());
}
}
void formJson(QString strJson, PHDFile *phd, rowvec& spec, vec& iC)
{
QJsonParseError err_rpt;
QJsonDocument root_Doc = QJsonDocument::fromJson(strJson.toStdString().c_str(), &err_rpt);//字符串格式化为JSON
if (err_rpt.error != QJsonParseError::NoError)
{
return;
}
QJsonObject rootObj = root_Doc.object();
QJsonArray arr = rootObj.value("vPeak").toArray();
phd->vPeak.clear();
for (int i = 0; i < arr.size(); i++)
{
QJsonObject obj = arr[i].toObject();
PeakInfo peak;
peak.peakCentroid = obj.value("peakCentroid").toDouble();
peak.fwhmc = obj.value("fwhmc").toDouble();
peak.tail = obj.value("tail").toDouble();
peak.tailAlpha = obj.value("tailAlpha").toDouble();
peak.upperTail = obj.value("upperTail").toDouble();
peak.upperTailAlpha = obj.value("upperTailAlpha").toDouble();
peak.area = obj.value("area").toDouble();
peak.stepRatio = obj.value("stepRatio").toDouble();
phd->vPeak.push_back(peak);
}
jsonArrToStdvec(rootObj, "usedEnerPara", phd->usedEnerPara.p);
jsonArrToStdvec(rootObj, "usedResoPara", phd->usedResoPara.p);
jsonArrToStdvec(rootObj, "XCtrl", phd->baseCtrls.XCtrl);
jsonArrToStdvec(rootObj, "YCtrl", phd->baseCtrls.YCtrl);
jsonArrToStdvec(rootObj, "YSlope", phd->baseCtrls.YSlope);
for (int i = 0; i < phd->baseCtrls.YSlope.size(); i++)
{
if (phd->baseCtrls.YSlope[i] == 0.0)
phd->baseCtrls.YSlope[i] = qQNaN();
}
jsonArrToStdvec(rootObj, "para_tail", phd->para_tail.p);
jsonArrToStdvec(rootObj, "para_tailAlpha", phd->para_tailAlpha.p);
jsonArrToStdvec(rootObj, "para_tailRight", phd->para_tailRight.p);
jsonArrToStdvec(rootObj, "para_tailRightAlpha", phd->para_tailRightAlpha.p);
jsonArrToStdvec(rootObj, "para_stepRatio", phd->para_stepRatio.p);
jsonArrToStdvec(rootObj, "usedEffiPara", phd->usedEffiPara.p);
phd->baseCtrls.rg_low = rootObj.value("rg_low").toInt();
phd->baseCtrls.rg_high = rootObj.value("rg_high").toInt();
phd->Spec.num_g_channel = rootObj.value("num_g_channel").toInt();
phd->Spec.begin_channel = rootObj.value("begin_channel").toInt();
//spec = v;// .clear();
stdvec v;
QJsonArray usedArr = rootObj.value("spec").toArray();
for (int i = 0; i < usedArr.size(); i++)
{
v.push_back(usedArr.at(i).toDouble());
}
spec = v;
stdvec v1;
QJsonArray iCArr = rootObj.value("iC").toArray();
for (int i = 0; i < iCArr.size(); i++)
{
v1.push_back(iCArr.at(i).toDouble());
}
iC = v1;
}
void formJson2(QString strJson, PHDFile *phd, rowvec& spec, vec& iC)
{
QJsonParseError err_rpt;
QJsonDocument root_Doc = QJsonDocument::fromJson(strJson.toStdString().c_str(), &err_rpt);//字符串格式化为JSON
if (err_rpt.error != QJsonParseError::NoError)
{
return;
}
QJsonObject rootObj = root_Doc.object();
phd->vPeak.clear();
QJsonArray arr1 = rootObj.value("peakCentroid").toArray();
QJsonArray arr2 = rootObj.value("fwhmc").toArray();
QJsonArray arr3 = rootObj.value("tail").toArray();
QJsonArray arr4 = rootObj.value("tailAlpha").toArray();
QJsonArray arr5 = rootObj.value("upperTail").toArray();
QJsonArray arr6 = rootObj.value("upperTailAlpha").toArray();
QJsonArray arr7 = rootObj.value("area").toArray();
QJsonArray arr8 = rootObj.value("stepRatio").toArray();
for (int i = 0; i < arr1.size(); i++)
{
PeakInfo peak;
peak.peakCentroid = arr1.at(i).toDouble();
peak.fwhmc = arr2.at(i).toDouble();
peak.tail = arr3.at(i).toDouble();
peak.tailAlpha = arr4.at(i).toDouble();
peak.upperTail = arr5.at(i).toDouble();
peak.upperTailAlpha = arr6.at(i).toDouble();
peak.area = arr7.at(i).toDouble();
peak.stepRatio = arr8.at(i).toDouble();
phd->vPeak.push_back(peak);
}
jsonArrToStdvec(rootObj, "usedEnerPara", phd->usedEnerPara.p);
jsonArrToStdvec(rootObj, "usedResoPara", phd->usedResoPara.p);
jsonArrToStdvec(rootObj, "XCtrl", phd->baseCtrls.XCtrl);
jsonArrToStdvec(rootObj, "YCtrl", phd->baseCtrls.YCtrl);
jsonArrToStdvec(rootObj, "YSlope", phd->baseCtrls.YSlope);
jsonArrToStdvec(rootObj, "para_tail", phd->para_tail.p);
jsonArrToStdvec(rootObj, "para_tailAlpha", phd->para_tailAlpha.p);
jsonArrToStdvec(rootObj, "para_tailRight", phd->para_tailRight.p);
jsonArrToStdvec(rootObj, "para_tailRightAlpha", phd->para_tailRightAlpha.p);
jsonArrToStdvec(rootObj, "para_stepRatio", phd->para_stepRatio.p);
jsonArrToStdvec(rootObj, "usedEffiPara", phd->usedEffiPara.p);
phd->baseCtrls.rg_low = rootObj.value("rg_low").toInt();
phd->baseCtrls.rg_high = rootObj.value("rg_high").toInt();
phd->Spec.num_g_channel = rootObj.value("num_g_channel").toInt();
phd->Spec.begin_channel = rootObj.value("begin_channel").toInt();
//spec = v;// .clear();
stdvec v;
QJsonArray usedArr = rootObj.value("vCount").toArray();
for (int i = 0; i < usedArr.size(); i++)
{
v.push_back(usedArr.at(i).toDouble());
}
spec = v;
stdvec v1;
QJsonArray iCArr = rootObj.value("iC").toArray();
for (int i = 0; i < iCArr.size(); i++)
{
v1.push_back(iCArr.at(i).toDouble());
}
iC = v1;
}
QJsonArray toJsonArr(stdvec val)
{
QJsonArray arr;
for (int i = 0; i < val.size(); i++)
{
arr.push_back(val[i]);
}
return arr;
}
QString toJsonObj(PHDFile *phd, rowvec spec, vec iC)
{
QJsonObject retobj;
QJsonArray arr;
foreach(const PeakInfo &peak, phd->vPeak)
{
QJsonObject obj;
obj.insert("peakCentroid", peak.peakCentroid);
obj.insert("fwhmc", peak.fwhmc);
obj.insert("tail", peak.tail);
obj.insert("tailAlpha", peak.tailAlpha);
obj.insert("upperTail", peak.upperTail);
obj.insert("upperTailAlpha", peak.upperTailAlpha);
obj.insert("area", peak.area);
obj.insert("stepRatio", peak.stepRatio);
arr.push_back(obj);
}
retobj.insert("vPeak", arr);
QJsonArray usedEnerPara = toJsonArr(phd->usedEnerPara.p);
QJsonArray usedResoPara = toJsonArr(phd->usedResoPara.p);
retobj.insert("usedEnerPara", usedEnerPara);
retobj.insert("usedResoPara", usedResoPara);
QJsonArray XCtrl = toJsonArr(phd->baseCtrls.XCtrl);
QJsonArray YCtrl = toJsonArr(phd->baseCtrls.YCtrl);
QJsonArray YSlope = toJsonArr(phd->baseCtrls.YSlope);
retobj.insert("XCtrl", XCtrl);
retobj.insert("YCtrl", YCtrl);
retobj.insert("YSlope", YSlope);
QJsonArray para_tail = toJsonArr(phd->para_tail.p);
QJsonArray para_tailAlpha = toJsonArr(phd->para_tailAlpha.p);
QJsonArray para_tailRight = toJsonArr(phd->para_tailRight.p);
QJsonArray para_tailRightAlpha = toJsonArr(phd->para_tailRightAlpha.p);
QJsonArray para_stepRatio = toJsonArr(phd->para_stepRatio.p);
QJsonArray usedEffiPara = toJsonArr(phd->usedEffiPara.p);
retobj.insert("para_tail", para_tail);
retobj.insert("para_tailAlpha", para_tailAlpha);
retobj.insert("para_tailRight", para_tailRight);
retobj.insert("para_tailRightAlpha", para_tailRightAlpha);
retobj.insert("para_stepRatio", para_stepRatio);
retobj.insert("usedEffiPara", usedEffiPara);
retobj.insert("rg_low", phd->baseCtrls.rg_low);
retobj.insert("rg_high", phd->baseCtrls.rg_high);
retobj.insert("num_g_channel", (int)phd->Spec.num_g_channel);
retobj.insert("begin_channel", (int)phd->Spec.begin_channel);
QJsonArray specArr;
for (int i = 0; i < spec.size(); i++)
{
specArr.push_back(spec[i]);
}
retobj.insert("spec", specArr);
QJsonArray iCArr;
for (int i = 0; i < iC.size(); i++)
{
iCArr.push_back(iC[i]);
}
retobj.insert("iC", iCArr);
QJsonDocument doc;
doc.setObject(retobj);
return QString(doc.toJson());
}
//
QVector<int> insertpeak(stdvec &vLeft, stdvec &vRight, PHDFile *phd, rowvec spec, vec iC, double minNA)
{
// QString str = toJsonObj(phd, spec, iC);
// formJson(str, phd, spec, iC);
//
// QFile file("E:/AWork/QINGH/insertVal.txt");
// file.open(QIODevice::ReadOnly);
// QByteArray bytearray = file.readAll();
// formJson(QString(bytearray), phd, spec, iC);
//
// rowvec dddy(phd->baseCtrls.YSlope);
// mat mt(dddy);
// uvec idxxxx = find(mt < 1000000.0);
// QString SSS = "";
// for (size_t i = 0; i < idxxxx.size(); i++)
// {
// if (SSS != "")
// SSS += ",";
// SSS += QString("%1").arg(idxxxx(i));
// }
uvec idx;
uword peakNum = phd->vPeak.size();
vec pC(peakNum), pFC(peakNum), pT(peakNum), pUT(peakNum), pNA(peakNum), pSt(peakNum), pTa(peakNum), pUTa(peakNum), pBWW(peakNum), pRBC(peakNum), pRDC(peakNum);
uword t_i = 0;
foreach (const PeakInfo &peak, phd->vPeak)
{
pC(t_i) = peak.peakCentroid;
pFC(t_i) = peak.fwhmc;
pT(t_i) = peak.tail;
pTa(t_i) = peak.tailAlpha;
pUT(t_i) = peak.upperTail;
pUTa(t_i) = peak.upperTailAlpha;
pNA(t_i) = peak.area;
pSt(t_i) = peak.area * peak.stepRatio;
++t_i;
}
pBWW.fill(gNaN);
pRBC.fill(gNaN);
pRDC.fill(gNaN);
// get resolution in channels
vec iE = Independ::calFcnEval(iC, phd->usedEnerPara.p);
vec idE = Independ::calDerivEval(iC, phd->usedEnerPara.p);
vec iR = Independ::calFcnEval(iE, phd->usedResoPara.p);
vec iRC = iR / idE;
//iE.print("iE = ");
//idE.print("idE = ");
//iR.print("iR = ");
//iRC.print("iRC = ");
// residual above 0
int nCounts = phd->Spec.num_g_channel;
if (phd->Spec.begin_channel == 0)
spec(nCounts - 1) = gNaN;
rowvec BaseChan = cumsum(ones<rowvec>(nCounts));
PiecePoly bpp = Independ::makeBasePP(phd->baseCtrls.XCtrl, phd->baseCtrls.YCtrl, phd->baseCtrls.YSlope, pC, pFC, pSt);
rowvec approx = Independ::peakBaseSpline(BaseChan, bpp, pNA, pC, pFC, pSt, pT, pTa, pUT, pUTa, pBWW, pRBC, pRDC);
rowvec resid = spec - approx;
uvec idx_resid = find(resid<0);
resid(idx_resid) = zeros<rowvec>(idx_resid.size());
// get net area estimate
vec iNA(iC.size());
for(uword i=0; i<iC.size(); ++i) //for i=1:length(iC);
{
// ch=[floor(iC(i)-iRC(i)): ceil(iC(i)+iRC(i))];
urowvec ch = RangeVec(floor(iC(i)-iRC(i))-1, ceil(iC(i)+iRC(i))-1);
//iNA(i) = eval('max(minNA,sum(resid(ch)))','minNA');
iNA(i) = max(minNA, sum(resid(ch)));
}
vec vT = Independ::calFcnEval(iE, phd->para_tail.p);
vec vTa = Independ::calFcnEval(iE, phd->para_tailAlpha.p);
vec vUT = Independ::calFcnEval(iE, phd->para_tailRight.p);
vec vUTa = Independ::calFcnEval(iE, phd->para_tailRightAlpha.p);
vec vSr = Independ::calFcnEval(iE, phd->para_stepRatio.p);
vec vEffi = Independ::calFcnEval(iE, phd->usedEffiPara.p);
idx.set_size(iC.n_elem);
for(uword i=0; i<iC.size(); ++i)
{
int j = 0;
while (j < peakNum && iC(i) > phd->vPeak[j].peakCentroid) ++j;
PeakInfo peak;
peak.peakCentroid = iC(i);
peak.energy = iE(i);
peak.area = iNA(i);
peak.sensitivity = gNaN;
peak.fwhm = iR(i);
peak.fwhmc = iRC(i);
peak.stepRatio = vSr(i);
peak.tail = vT(i);
peak.tailAlpha = vTa(i);
peak.upperTail = vUT(i);
peak.upperTailAlpha = vUTa(i);
peak.efficiency = vEffi(i);
peak.BWWidthChan = 0;
peak.recoilBetaChan = gNaN;
peak.recoilDeltaChan = gNaN;
phd->vPeak.insert(phd->vPeak.begin() + j, peak);
idx(i) = j;
pC.insert_rows(j, 1);
pFC.insert_rows(j, 1);
pC(j) = iC(i);
pFC(j) = iRC(i);
}
mat vvRg = findPeakRange(idx, pC, pFC, phd->baseCtrls.rg_low, phd->baseCtrls.rg_high);
vLeft = conv_to<stdvec>::from(vvRg.col(0));
vRight = conv_to<stdvec>::from(vvRg.col(1));
QVector<int> vRetIdx;
for(uword i=0; i<idx.n_elem; ++i)
{
vRetIdx.push_back(idx(i));
}
return vRetIdx;
}
mat findPeakRange(uvec pnr, vec C, vec F_Ch, double rg_low, double rg_high)
{
// get all peak centroids and FWHM in channel units
//[C, F_Ch] = dmps(sn, pat, 'Centroid', 'FWHM_Ch');
rowvec ar; //ar = dmspec(sn, 'AnalysisRange');
ar << rg_low << rg_high;
uword n = pnr.size();
mat rg = zeros(n, 2);
vec iL = zeros<vec>(n);
vec iH = zeros<vec>(n);
// left and right borders of multiplet
for(uword i=0; i<n; ++i) //for i = 1 : n
{
//[rg(i, 1), rg(i, 2), iL(i), iH(i)] = findVPeakRange(pnr(i), C, F_Ch, ar, varargin{:});
rowvec temp = findVPeakRange(pnr(i), C, F_Ch, ar);
rg(i, 0) = temp(0);
rg(i, 1) = temp(1);
iL(i) = temp(2);
iH(i) = temp(3);
}
return rg;
}
rowvec findVPeakRange(int pNr, vec C, vec FC, rowvec rg, double LD, double RD)
{
int NP = C.size();
// find leftmost peak in multiplet, and minimum left border
int i = pNr;
double l = C(i) - LD * FC(i);
while (i > 0 && l < C(i-1) + RD * FC(i-1))
{
i = i - 1;
l = min(l, C(i) - LD * FC(i));
}
// find rightmost peak in multiplet, and maximum right border
int j = pNr;
double h = C(j) + RD * FC(j);
while (j < NP-1 && h > C(j+1) - LD * FC(j+1))
{
j = j + 1;
h = max(h, C(j) + RD * FC(j));
}
// left and right borders of multiplet
l = max(rg(0), floor(l));
h = min(rg(1), ceil(h));
rowvec results;
results << l << h << i << j;
return results;
}
void fitPeakFull(PHDFile *phd, uvec Af, uvec Cf, uvec Ff)
{
int nCount = phd->Spec.num_g_channel;
rowvec spec(nCount);
for(int i=0; i<nCount; ++i)
{
spec(i) = phd->Spec.counts[i];
}
if(phd->Spec.begin_channel == 0)
{
spec.insert_cols(nCount, 1);
spec(nCount) = gNaN;
spec.shed_col(0);
}
// get index to all peaks where any parameter is free
// get common min and max
uvec allidx = join_vert( join_vert(Af, Cf), Ff );
if(allidx.is_empty()) return;
// get unique allidx
ivec Bool = zeros<ivec>( max(allidx) + 1u);
Bool(allidx).fill(1);
allidx = find(Bool);
int ip = 0, np = phd->vPeak.size();
vec C(np), FC(np), NA(np), Sr(np), T(np), Ta(np), UT(np), UTa(np), BW(np), RB(np), RD(np);
foreach (const PeakInfo &peak, phd->vPeak)
{
C(ip) = peak.peakCentroid;
FC(ip) = peak.fwhmc;
NA(ip) = peak.area;
Sr(ip) = peak.stepRatio;
T(ip) = peak.tail;
Ta(ip) = peak.tailAlpha;
UT(ip) = peak.upperTail;
UTa(ip) = peak.upperTailAlpha;
BW(ip) = peak.BWWidthChan;
RB(ip) = peak.recoilBetaChan;
RD(ip) = peak.recoilDeltaChan;
++ip;
}
// PrintColvec("na.txt", NA);
rowvec cx = phd->baseCtrls.XCtrl;
rowvec cy = phd->baseCtrls.YCtrl;
rowvec cdy = phd->baseCtrls.YSlope;
// get fitting region
rowvec x = RangeVec2(phd->baseCtrls.rg_low, phd->baseCtrls.rg_high); // x = [rg(1) : rg(2)];
rowvec t_x = x - 1;
rowvec y = IdxMat(spec, t_x);
/* % fit parameters
[s, CXo, cy, cdy, Ao, Co, Fo, Sro, To, Tao, UTo, UTao, BWo, RBo, RDo, X2] = ...
fitFunction(x, y, 'wrapStepRatio', ...
cx, [], cy, CPf, cdy, [], NA, Af, C, Cf, FC, Ff, Sr, SRf, ...
T, [], Ta, [], UT, [], UTa, [], BW, [], RB, [], RD, [], ...
'Restrict', addArg{:});*/
double X2;
uvec ep;
field<vec> Ps;
field<uvec> free_idx;
Ps << vectorise(cx) << vectorise(cy) << vectorise(cdy) << NA << C << FC << Sr << T << Ta << UT << UTa << BW << RB << RD;
// PrintMat2fvec("Ps1.txt", Ps);
free_idx << ep << ep << ep << Af << Cf << Ff << ep << ep << ep << ep << ep << ep << ep << ep;
QString addArg = "Peaks above spline baseline";
bool s = Independ::fitFunction(Ps, X2, free_idx, vectorise(x), vectorise(y), "wrapStepRatio", QStringList("Restrict"), addArg);
// PrintMat2fvec("Ps1.txt", Ps);
/*mat& CXo = Ps(0);
//if(Ps(1).is_colvec()) cy = Ps(1).t();
//else cy = Ps(1);
//if(Ps(2).is_colvec()) cdy = Ps(2).t();
//else cdy = Ps(2);*/
mat& Ao = Ps(3);
// PrintMat22("area.txt",Ao);
mat& Co = Ps(4);
mat& Fo = Ps(5);
mat& Sro = Ps(6);
mat& To = Ps(7);
mat& Tao = Ps(8);
mat& UTo = Ps(9);
mat& UTao = Ps(10);
mat& BWo = Ps(11);
mat& RBo = Ps(12);
mat& RDo = Ps(13);
// only write to temp pat if fitting converged
if(s)
{
double k = phd->setting.k_back;
double k_alpha = phd->setting.k_alpha;
double k_beta = phd->setting.k_beta;
// get baseline and peak parameters
rowvec BaseChan = RangeVec2(1, nCount);
PiecePoly bpp = Independ::makeBasePP(vectorise(cx), vectorise(cy), vectorise(cdy), Co, Fo, Ao%Sro);
rowvec BL = Matlab::ppval(BaseChan, bpp)+Independ::peakSteps(BaseChan, Co, Fo, Ao%Sro);
phd->vBase = conv_to<stdvec>::from(BL);
vec Sc = Fo/sqrt(8*log(2));
/* // modify all peaks, if idx is string
if isstr(idx)
idx = [1 : length(Ct)];
end*/
if(Co.is_empty()) return; // no peaks in PAT or called needlessly
if(allidx.is_empty())
{
allidx = vectorise( RangeVec(0, Co.size()-1) );
}
vec MBC, LC, LD, BC;
// restrict to indexed peaks
Co = Co(allidx); //Ct = Ct(idx);
Fo = Fo(allidx); //Fc = Fc(idx);
Ao = Ao(allidx); //NA = NA(idx);
// PrintMat22("ao_allidx.txt",Ao);
Sro = Sro(allidx);
To = To(allidx);
Tao = Tao(allidx);
UTo = UTo(allidx);
UTao = UTao(allidx);
BWo = BWo(allidx);
RBo = RBo(allidx);
RDo = RDo(allidx);
Sc = Sc(allidx); //Sc = Sc(idx);
vec W = 2 * k * Fo; //W = 2 * k * Fc;
// total background counts from Ct +- k*Fc
BC = Independ::abkcnt( max(0, BL), Co, Fo, k); //BC = abkcnt(max(BL, 0), Ct, Fc, k);
// Mean Back Counts
MBC = BC / W;
double c = 1.1;
double f = 1.9562;
vec DA = c*f*sqrt(Sc % MBC); //DA = c*f*sqrt(Sc .* MBC);
// Currie's LC
LC = k_alpha * DA;
// Currie's LD
// old formula: LD = k_alpha^2 + 2* LC;
double kb2 = pow(k_beta, 2);
// assymmetric formula, ref TBD
LD = LC + (kb2/2) * ( 1 + sqrt(1 + (4/kb2)*( LC+pow(DA,2) ) ) );
mat dE = Independ::calDerivEval(Co, phd->usedEnerPara.p);
vec vE = Independ::calFcnEval(Co, phd->usedEnerPara.p);
for(uword i=0; i<allidx.size(); ++i)
{
phd->vPeak[ allidx(i) ].peakCentroid = Co(i);
phd->vPeak[ allidx(i) ].energy = vE(i);
phd->vPeak[ allidx(i) ].area = Ao(i);
phd->vPeak[ allidx(i) ].areaErr = sqrt( max( Ao(i), LC(i) ) + BC(i) );
phd->vPeak[ allidx(i) ].fwhm = Fo(i) * dE(i);
phd->vPeak[ allidx(i) ].fwhmc = Fo(i);
phd->vPeak[ allidx(i) ].significance = Ao(i) / LC(i);
phd->vPeak[ allidx(i) ].meanBackCount = MBC(i);
phd->vPeak[ allidx(i) ].backgroundArea = BC(i);
phd->vPeak[ allidx(i) ].lc = LC(i);
phd->vPeak[ allidx(i) ].ld = LD(i);
phd->vPeak[ allidx(i) ].tail = To(i);
phd->vPeak[ allidx(i) ].tailAlpha = Tao(i);
phd->vPeak[ allidx(i) ].upperTail = UTo(i);
phd->vPeak[ allidx(i) ].upperTailAlpha = UTao(i);
phd->vPeak[ allidx(i) ].stepRatio = Sro(i);
phd->vPeak[ allidx(i) ].BWWidthChan = BWo(i);
phd->vPeak[ allidx(i) ].recoilBetaChan = RBo(i);
phd->vPeak[ allidx(i) ].recoilDeltaChan = RDo(i);
}
}
}
void UpdateBaseControl(BaseControls& bc)
{
rowvec cx(bc.XCtrl);
rowvec cy(bc.YCtrl);
rowvec cdy(bc.YSlope);
rowvec base = rowvec(bc.Baseline);
rowvec sc = rowvec(bc.StepCounts);
rowvec tmp_cx = cx - 1;
rowvec cyr = cy - IdxMat(sc, tmp_cx); //cyr = cy - sc(cx);
rowvec chans = RangeVec2(bc.rg_low, bc.rg_high);//chans = [rg(1):rg(2)];
base.cols(bc.rg_low-1, bc.rg_high-1) = sc.cols(bc.rg_low-1, bc.rg_high-1) + Independ::wrapBaseSpline(chans, cx, cyr, cdy);
bc.Baseline = conv_to<stdvec>::from(base);
//hb = plot(e, bc, 'm-');
//hc = plot(ce, cy, 'ms');
//cpd.SlopeHandle = plot([ceLo; ceHi], [cy-2*cdy;cy+2*cdy], 'g+-');
//uvec idx = find_finite(cdy); //idx = find(isfinite(cdy));
//vec tmp1 = cy(idx), tmp2 = 2 * cdy(idx);
//bc.SlopeHandle.clear();
//bc.SlopeHandle = conv_to<stdvec>::from(join_vert(tmp1-tmp2,tmp1+tmp2));
}
stdvec calValuesOut(rowvec Chan, vec p)
{
return conv_to<stdvec>::from(Independ::calFcnEval(Chan, p));
}
double calDerivaOut(double Chan, vec p)
{
if(p.size() < 2) return 0;
vec x; x << Chan;
vec dy = Independ::calDerivEval(x, p);
return as_scalar(dy);
}
ParameterInfo calFitPara(CalibrationType cal, int funcId, stdvec x, stdvec y, stdvec err)
{
ParameterInfo param;
int s = x.size();
if(s != y.size())
{
qDebug() << "The size of param x and y must be same in function calFitPara";
return param;
}
if(err.size() != s) err.assign(s, gNaN);
vec p, perr;
if(funcId == 1)
{
Independ::lMakeInterp(p, perr, x, y, err);
}
else {
vec p0;
p0 << funcId;
p0 = join_vert(p0, Independ::calDefIniPara(funcId, x, y));
//p0.print("p0 = ");
Independ::calFitPara(p, perr, cal, x, y, err, p0);
}
if(Independ::calParaCheck(p))
{
param.p = conv_to<stdvec>::from(p);
param.perr = conv_to<stdvec>::from(perr);
}
return param;
}
stdvec interp1(const PeakInfo &peak, vec regBase, vec u)
{
vec x = vectorise(RangeVec2(peak.left, peak.right));
PiecePoly pp;
vec na, ct, fc, st, tail, ta, ut, uta, bw, rb, rd;
na << peak.area;
ct << peak.peakCentroid;
fc << peak.fwhmc;
st << 0;
tail << peak.tail;
ta << peak.tailAlpha;
ut << peak.upperTail;
uta << peak.upperTailAlpha;
bw << peak.BWWidthChan;
rb << peak.recoilBetaChan;
rd << peak.recoilDeltaChan;
vec y = regBase + vectorise(Independ::peakBaseSpline(x, pp, na, ct, fc, st, tail, ta, ut, uta, bw, rb, rd));
stdvec v;
//bool ix = 1; // Is x given as the first argument?
// The v5 option, '*method', asserts that x is equally spaced.
bool eqsp = false;
QString extrapval = "extrap";
// check for NaN's
if (!x.is_finite())
{
qDebug() << "MATLAB:interp1:NaNinX','NaN is not an appropriate value for X.";
}
vec ua = u;
// NANS are allowed as a value for F(X), since a function may be undefined
// for a given value.
if (!y.is_finite())
{
qDebug() << "MATLAB:interp1:NaNinY, NaN found in Y, interpolation at undefined values will result in undefined values.";
}
// Check dimensions. Work with column vectors.
//SizeCube siz_y = size(y); // obtaining the size of y in a vector
uword m = y.n_rows;
//uword n = y.n_cols;
vec h;
if (x.size() != m)
{
qDebug() << "MATLAB:interp1:YInvalidNumRows','Y must have length(X) rows.";
}
if (!eqsp)
{
h = diff(x);
eqsp = (norm(diff(h), "inf") <= eps(norm(x, "inf")));
if (!x.is_finite())
{
eqsp = 0; // if an INF in x, x is not equally spaced
}
}
if (eqsp)
{
h << (x(m-1)-x(0)) / (m-1);
}
if (m < 2)
{
if (ua.is_empty())
{
return v;
} else {
qDebug() << "MATLAB:interp1:NotEnoughPts','There should be at least two data points.";
}
}
//SizeCube siz = size(ua);
vec p;
// Interpolate
rowvec t_y = y.t();
//v = reshape(pchip(x,y,u),prod(siz_y(2:end)),prod(siz)).';
v = conv_to<stdvec>::from(reshape(Independ::pchip(x, t_y, u), size(u)));
return v;
}
bool CalculateMDCs(PHDFile *phd, NuclideActMda &nucActMda, int mainPeakIdx, double lambda, double keyLineYield, double CCF)
{
//double Ld = 2 * Fwhmc * (Lcc - baseline) / 0.8591;
//double MDA = Ld * CCF * DCF / (Yi * Effki * Tl);
//double MDC = MDA / Svol;
// 计算衰变校正因子——DCF
QDateTime collectStart = QDateTime::fromString(phd->collect.collection_start_date + " " + phd->collect.collection_start_time, QString(DATATIME_FORMAT));
QDateTime collectStop = QDateTime::fromString(phd->collect.collection_stop_date + " " + phd->collect.collection_stop_time, QString(DATATIME_FORMAT));
QDateTime acqStart = QDateTime::fromString(phd->acq.acquisition_start_date + " " + phd->acq.acquisition_start_time, QString(DATATIME_FORMAT));
double Ts = collectStart.secsTo(collectStop); // 采样时间
double Td = collectStop.secsTo(acqStart); // 衰变时间
double Ta = phd->acq.acquisition_real_time; // 能谱获取实时间
double Tl = phd->acq.acquisition_live_time; // 能谱获取活时间
double Svol = phd->collect.air_volume; // 样品采样体积
double DCF1, DCF2, DCF3;
if ( Ts == 0 ) DCF1 = 1;
else DCF1 = lambda * Ts / (1-exp(-lambda*Ts));
if ( Td == 0 ) DCF2 = 1;
else DCF2 = exp(lambda*Td);
if ( Ta == 0 ) DCF3 = 1;
else DCF3 = lambda * Ta / (1-exp(-lambda*Ta));
double DCF_act = exp(lambda * phd->usedSetting.refTime_act.secsTo(acqStart));
double DCF_conc = exp(lambda * phd->usedSetting.refTime_conc.secsTo(collectStart));
//double DCF = DCF1 * DCF2 * DCF3;
// calculate line activities
const PeakInfo &peak = phd->vPeak[mainPeakIdx];
double netKeyPeakArea = peak.area;
double linePeakAreaErr = peak.areaErr;
double detectorEfficiency = peak.efficiency;
double detectorEfficiencyError = DETECTOR_EFFICIENCY_ERROR_RATIO * peak.efficiency;
double fwhmc = peak.fwhmc;
int index = (int)(peak.peakCentroid + 0.5) - 1;
double lcc = phd->vLc[index];
double baseline = phd->vBase[index];
double activity, mda;
if ( keyLineYield == 0 )
{
activity = 0;
mda = 0;
}
else
{
activity = ( netKeyPeakArea * CCF * DCF3*DCF_act ) / ( keyLineYield * detectorEfficiency * Tl ) ;//Bq,参考时间为referenceTimeActivity
mda = (2 * fwhmc * (lcc - baseline) / 0.8591 ) * CCF * DCF3*DCF_act / ( keyLineYield * detectorEfficiency * Tl ) ;//Bq,参考时间为referenceTimeActivity
}
double concentration = ( netKeyPeakArea * CCF*DCF1*DCF2*DCF3*DCF_conc ) / ( keyLineYield * detectorEfficiency * Tl ) *1e6/ Svol;//mBq/m3,参考时间为referenceTimeConc
double mdc = (2 * fwhmc * (lcc - baseline) / 0.8591 ) * CCF*DCF1*DCF2*DCF3*DCF_conc / ( keyLineYield * detectorEfficiency * Tl ) *1e6 / Svol;
double activityErr = sqrt( pow(linePeakAreaErr/netKeyPeakArea, 2) + pow(detectorEfficiencyError/detectorEfficiency, 2) ) * activity;
//double concentrationErr = activityErr / Svol;
nucActMda.activity = activity;
nucActMda.act_err = activityErr;
nucActMda.mda = mda;
nucActMda.mdc = mdc;
nucActMda.efficiency = peak.efficiency;
nucActMda.effi_err = detectorEfficiencyError;
nucActMda.concentration = concentration;
return (activity > mda);
}
double CalculateMDC(PHDFile *phd, stdvec vMdcInfo, double CCF)
{
if(vMdcInfo.size() < 3 || vMdcInfo[2] == 0) return 0;
QDateTime collectStart = QDateTime::fromString(phd->collect.collection_start_date + " " + phd->collect.collection_start_time, QString(DATATIME_FORMAT));
QDateTime collectStop = QDateTime::fromString(phd->collect.collection_stop_date + " " + phd->collect.collection_stop_time, QString(DATATIME_FORMAT));
QDateTime acqStart = QDateTime::fromString(phd->acq.acquisition_start_date + " " + phd->acq.acquisition_start_time, QString(DATATIME_FORMAT));
double Ts = collectStart.secsTo(collectStop); // 采样时间
double Td = collectStop.secsTo(acqStart); // 衰变时间
double Ta = phd->acq.acquisition_real_time; // 能谱获取实时间
double Tl = phd->acq.acquisition_live_time; // 能谱获取活时间
double Svol = phd->collect.air_volume; // 样品采样体积
double DCF1, DCF2, DCF3;
double lambda = log(2.0) / (vMdcInfo[2] * 86400);
if ( Ts == 0 ) DCF1 = 1;
else DCF1 = lambda * Ts / (1-exp(-lambda*Ts));
if ( Td == 0 ) DCF2 = 1;
else DCF2 = exp(lambda*Td);
if ( Ta == 0 ) DCF3 = 1;
else DCF3 = lambda * Ta / (1-exp(-lambda*Ta));
//double DCF = DCF1 * DCF2 * DCF3;
//double DCF_act = exp(lambda * phd->setting.refTime_act.secsTo(acqStart));
double DCF_conc = exp(lambda * phd->usedSetting.refTime_conc.secsTo(collectStart));
// calculate line activities
//PeakInfo &peak = phd->vPeak[mainPeakIdx];
vec energy; energy << vMdcInfo[0];
stdvec channel = energyToChannel(energy, phd->usedEnerPara.p);
vec dE = Independ::calDerivEval(channel, phd->usedEnerPara.p);
double fwhmc = as_scalar(Independ::calFcnEval(energy, phd->usedResoPara.p) / dE);
double effi = as_scalar( Independ::calFcnEval(energy, phd->usedEffiPara.p) );
size_t index = 0;
for(size_t i=1; i<phd->vEnergy.size(); ++i)
{
if(phd->vEnergy[i] >= vMdcInfo[0])
{
index = i;
if(phd->vEnergy[i] - vMdcInfo[0] > vMdcInfo[0] - phd->vEnergy[i-1]) index = i-1;
break;
}
}
double lcc = phd->vLc[index];
double baseline = phd->vBase[index];
/*double mda = 0;
if ( vMdcInfo[1] != 0 )
{
mda = (2 * fwhmc * (lcc - baseline) / 0.8591 ) * CCF * DCF / ( vMdcInfo[1] * effi * Tl ) * 1e06;
mda = (2 * fwhmc * (lcc - baseline) / 0.8591 ) * CCF * DCF3*DCF_act / ( vMdcInfo[1] * effi * Tl ) ;
}
double concentration = activity / Svol;
double mdc = mda / Svol;*/
double mdc = (2 * fwhmc * (lcc - baseline) / 0.8591 ) * CCF*DCF1*DCF2*DCF3*DCF_conc / ( vMdcInfo[1] * effi * Tl ) *1e6 / Svol;
return mdc;
}
bool ReadQCLimit(QMap<QString, QcCheckItem>& qcItems, double &ener_Be7, stdvec &vMdcInfo, QString system_type, QString strFPath)
{
qcItems.clear();
vMdcInfo.clear();
//QFile file(qApp->applicationDirPath() + "/SystemManager.xml");
QFile file(strFPath + "/SystemManager.xml");
if(!file.open(QIODevice::ReadOnly))
{
return false;
}
QDomDocument doc;
if(!doc.setContent(&file))
{
file.close();
return false;
}
file.close();
QDomNode node = doc.documentElement().firstChild();
while(!node.isNull())
{
if(node.nodeName() == "QCFlags-P" && system_type.toUpper() == "P")
{
QDomNode child = node.firstChild();
while(!child.isNull())
{
QString name = child.nodeName();
if(name == QC_COL_TIME || name == QC_ACQ_TIME || name == QC_DECAY_TIME || name == QC_SAMP_VOL
|| name == QC_AIR_FLOW || name == QC_Be7_FWHM || name == QC_Ba140_MDC)
{
qcItems[name].standard = child.toElement().attribute("green");
}
else if(name == "Be7")
{
ener_Be7 = child.toElement().attribute("energy").toDouble();
}
else if(name == "Ba140")
{
vMdcInfo.push_back(child.toElement().attribute("energy").toDouble());
vMdcInfo.push_back(child.toElement().attribute("yield").toDouble());
vMdcInfo.push_back(child.toElement().attribute("halflife").toDouble());
}
child = child.nextSibling();
}
break;
}
else if(node.nodeName() == "QCFlags-G" && system_type.toUpper() == "G")
{
QDomNode child = node.firstChild();
while(!child.isNull())
{
QString name = child.nodeName();
if(name == QC_COL_TIME || name == QC_ACQ_TIME || name == QC_DECAY_TIME || name == QC_SAMP_VOL
|| name == QC_AIR_FLOW || name == QC_Xe133_MDC)
{
qcItems[name].standard = child.toElement().attribute("green");
}
else if(name == "Xe133")
{
vMdcInfo.push_back(child.toElement().attribute("energy").toDouble());
vMdcInfo.push_back(child.toElement().attribute("yield").toDouble());
vMdcInfo.push_back(child.toElement().attribute("halflife").toDouble());
}
child = child.nextSibling();
}
break;
}
node = node.nextSibling();
}
return true;
}
void RunQC(PHDFile *phd, QString strFPath)
{
QDateTime start = QDateTime::fromString(phd->collect.collection_start_date + " " + phd->collect.collection_start_time, QString(DATATIME_FORMAT));
QDateTime end = QDateTime::fromString(phd->collect.collection_stop_date + " " + phd->collect.collection_stop_time, QString(DATATIME_FORMAT));
QDateTime acq = QDateTime::fromString(phd->acq.acquisition_start_date + " " + phd->acq.acquisition_start_time, QString(DATATIME_FORMAT));
double collect_hour = start.secsTo(end) / 3600.0;
double acq_hour = phd->acq.acquisition_real_time / 3600.0;
double Decay_hour = end.secsTo(acq) / 3600.0;
double ener_Be7;
stdvec vMdcInfo;
QMap<QString, QcCheckItem> qcItems;
if(!ReadQCLimit(qcItems, ener_Be7, vMdcInfo, phd->header.system_type, strFPath))
{
qDebug() << "Read QC Flags from SystemManager.xml Failed!";
}
qcItems[QC_COL_TIME].value = collect_hour;
qcItems[QC_ACQ_TIME].value = acq_hour;
qcItems[QC_DECAY_TIME].value = Decay_hour;
qcItems[QC_SAMP_VOL].value = phd->collect.air_volume;
qcItems[QC_AIR_FLOW].value = phd->collect.air_volume / collect_hour;
if(phd->isValid && phd->vBase.size() == phd->Spec.num_g_channel)
{
if(phd->header.system_type == "P")
{
vec energy; energy << ener_Be7;
vec fwhm = Independ::calFcnEval(energy, phd->usedResoPara.p);
qcItems[QC_Be7_FWHM].value = fwhm(0);
qcItems[QC_Ba140_MDC].value = CalculateMDC(phd, vMdcInfo, 1.0);
}
else qcItems[QC_Xe133_MDC].value = CalculateMDC(phd, vMdcInfo, 1.0);
}
for(QMap<QString, QcCheckItem>::iterator iter = qcItems.begin(); iter != qcItems.end(); ++iter)
{
if(iter.value().standard.isEmpty()) continue;
QStringList lists = iter.value().standard.split(',', QString::SkipEmptyParts);
bool bSatisfy = true;
foreach(QString str, lists)
{
if(str.contains('-')) continue;
else if(str.contains('('))
{
if(iter.value().value <= str.remove('(').toDouble()) { bSatisfy = false; break; }
}
else if(str.contains(')'))
{
if(iter.value().value >= str.remove(')').toDouble()) { bSatisfy = false; break; }
}
else if(str.contains('['))
{
if(iter.value().value < str.remove('[').toDouble()) { bSatisfy = false; break; }
}
else if(str.contains(']'))
{
if(iter.value().value > str.remove(']').toDouble()) { bSatisfy = false; break; }
}
}
iter.value().bPass = bSatisfy;
}
phd->QcItems = qcItems;
}
void FilterNuclideLine(NuclideLines &lines, double lowE)
{
stdvec vE = lines.vEnergy;
for(size_t i=0, j=0; i<vE.size(); ++i)
{
if(vE[i] < lowE)
{
lines.fullNames.removeAt(j);
lines.vEnergy.erase(lines.vEnergy.begin()+j);
lines.vYield.erase(lines.vYield.begin()+j);
lines.vUncertE.erase(lines.vUncertE.begin()+j);
lines.vUncertY.erase(lines.vUncertY.begin()+j);
if(i == lines.key_flag) lines.key_flag = -1;
else if(i < lines.key_flag) --lines.key_flag;
}
else ++j;
}
if(lines.key_flag < 0 && lines.vEnergy.size() > 0)
{
stdvec &vY = lines.vYield;
lines.maxYeildIdx = 0;
double maxYield = vY[0];
for(size_t ii=1; ii<vY.size(); ++ii)
{
if(vY[ii] > maxYield)
{
maxYield = vY[ii];
lines.maxYeildIdx = ii;
}
}
}
else lines.maxYeildIdx = lines.key_flag;
}
void ReadSpecialNuclides(QMap<QString, double> &mapHalflife, QStringList &vNuclides,QString xmlPath = "")
{
QString fileName = xmlPath + "/nuclide_ActMdc.txt";
QFile t_file(fileName);
if(!t_file.open(QIODevice::ReadOnly))
{
return;
}
QTextStream in(&t_file);
while(!in.atEnd())
{
QString line = in.readLine();
if(line.contains("#MDA"))
{
line = in.readLine();
while(!line.isEmpty() && !line.contains('#'))
{
QStringList strList = line.split(QRegExp("\\s"), QString::SkipEmptyParts);
if(strList.size() == 3)
{
mapHalflife[ strList[0] ] = strList[2].toDouble() * 86400;
}
line = in.readLine();
}
}
if(line.contains("#Identify"))
{
line = in.readLine();
while(!line.isEmpty() && !line.contains('#'))
{
QStringList strList = line.split(QRegExp("\\s"), QString::SkipEmptyParts);
vNuclides.append(strList);
line = in.readLine();
}
}
}
}
void NuclidesIdent(PHDFile* phd, QMap<QString, NuclideLines> &mapLines, QString xmlPath)
{
// 当重新分析时先清除上一次的分析结果
phd->mapNucActMda.clear();
for(size_t i=0; i<phd->vPeak.size(); ++i) phd->vPeak[i].nuclides.clear();
// 获取需特殊处理的核素
QMap<QString, double> mapHalflife; // 用其他核素半衰期计算活度/浓度的核素
QStringList vNuclides; // 只识别不计算活度/浓度的核素
ReadSpecialNuclides(mapHalflife, vNuclides, xmlPath);
double energyWidth = phd->usedSetting.EnergyTolerance;
int peakNum = phd->vPeak.size();
for (QMap<QString, NuclideLines>::iterator iter = mapLines.begin(); iter != mapLines.end(); ++iter)
{
if(iter.value().halflife <= 0) continue;
FilterNuclideLine(iter.value(), phd->usedSetting.ECutAnalysis_Low); // 过滤核素能量小于ECutAnalysis_Low的射线
const stdvec &vEnergy = iter.value().vEnergy; // 该核素的所有γ射线能量
const stdvec &vYield = iter.value().vYield;
stdvec vEffi = calValuesOut(vEnergy, phd->usedEffiPara.p); // 该核素所有γ射线能量处的探测效率
std::vector<int> vFit(vEnergy.size(), -1); // γ射线能量与峰中心道能量匹配标识
double sum_total = 0; // 该核素所有γ射线能量处效率乘以分支比的和
double sum_found = 0; // 所有匹配的γ射线能量处效率乘以分支比的和
int mainPeakIdx = -1; // 记录核素主γ峰的索引下标
for (size_t i=0, j=0; i<vEnergy.size(); ++i)
{
for(; j<peakNum; ++j)
{
if(phd->vPeak[j].energy >= 510 && phd->vPeak[j].energy <= 512)
{
continue; // 峰中心道能量为511的峰不进行核素识别
}
if(vEnergy[i] < phd->vPeak[j].energy - energyWidth) break;
else if(vEnergy[i] <= phd->vPeak[j].energy + energyWidth)
{
sum_found += vEffi[i] * vYield[i];
vFit[i] = j;
if(iter.value().maxYeildIdx == i) mainPeakIdx = j;
break;
}
}
sum_total += vEffi[i] * vYield[i];
}
// 核素匹配到峰
if(sum_total > 0)
{
// 如果该核素属特殊核素,则用“特殊核素配置文件”中指明的其他核素的半衰期
double halflife = (mapHalflife.find(iter.key()) == mapHalflife.end() ? iter.value().halflife : mapHalflife[iter.key()]);
double lambda = log(2.0) / halflife;
double rate = sum_found / sum_total;
if(rate > NUCL_IDENT_VALUE2)
{
// 获取用于计算Activity、MDC的主γ峰和最大分支比
double maxFoundYield = vYield[iter.value().maxYeildIdx];
if(mainPeakIdx < 0)
{
maxFoundYield = 0;
for(size_t ii=0; ii<vFit.size(); ++ii)
{
if(vFit[ii] >= 0 && vYield[ii] > maxFoundYield)
{
mainPeakIdx = vFit[ii];
maxFoundYield = vYield[ii];
}
}
if(mainPeakIdx < 0) continue;
}
NuclideActMda ActMda;
bool bActBigger = CalculateMDCs(phd, ActMda, mainPeakIdx, lambda, maxFoundYield, 1);
if(rate > NUCL_IDENT_VALUE1 || bActBigger)
{
ActMda.halflife = halflife;
ActMda.key_flag = -1;
if(!vNuclides.contains(iter.key())) // 需要计算活度浓度的核素
{
ActMda.bCalculateMDA = true;
}
else {
ActMda.activity = 0;
ActMda.act_err = 0;
ActMda.efficiency = 0;
ActMda.effi_err = 0;
ActMda.mda = 0;
ActMda.mdc = 0;
ActMda.concentration = 0;
}
int fitLineNum = 0, peakIdx = -1;
for(size_t ii=0; ii<vFit.size(); ++ii)
{
if(vFit[ii] >= 0)
{
// 向峰信息表中添加识别到的核素
if(vFit[ii] != peakIdx)
{
phd->vPeak[ vFit[ii] ].nuclides.push_back( iter.key() );
peakIdx = vFit[ii];
}
// 添加匹配的γ射线的信息
if(vFit[ii] == mainPeakIdx) ActMda.calculateIdx = fitLineNum;
if(iter.value().maxYeildIdx == ii && iter.value().key_flag >= 0) ActMda.key_flag = fitLineNum;
ActMda.vPeakIdx.push_back(peakIdx+1);
ActMda.fullNames.push_back(iter.value().fullNames[ii]);
ActMda.vEnergy.push_back(vEnergy[ii]);
ActMda.vUncertE.push_back(iter.value().vUncertE[ii]);
ActMda.vYield.push_back(vYield[ii]);
ActMda.vUncertY.push_back(iter.value().vUncertY[ii]);
++fitLineNum;
}
}
phd->mapNucActMda[iter.key()] = ActMda;
}
}
}
}
}
void CalcNuclideMDA(PHDFile* phd, NuclideLines &lines, const QString &nucName, std::vector<int> &vPeakIdx)
{
if(lines.halflife <= 0) return;
// 过滤核素能量小于ECutAnalysis_Low的射线
FilterNuclideLine(lines, phd->usedSetting.ECutAnalysis_Low);
// 获取需特殊处理的核素
QMap<QString, double> mapHalflife; // 用其他核素半衰期计算活度/浓度的核素
QStringList vNuclides; // 只识别不计算活度/浓度的核素
ReadSpecialNuclides(mapHalflife, vNuclides);
double energyWidth = phd->usedSetting.EnergyTolerance;
const stdvec &vEnergy = lines.vEnergy; // 该核素的所有γ射线能量
double maxYield = 0;
int mainPeakIdx = -1; // 记录核素主γ峰的索引下标
NuclideActMda ActMda;
ActMda.halflife = (mapHalflife.find(nucName) == mapHalflife.end() ? lines.halflife : mapHalflife[nucName]);
for (size_t i=0, j=0; i<vEnergy.size(); ++i)
{
for(; j<vPeakIdx.size(); ++j)
{
double energy = phd->vPeak.at(vPeakIdx[j]).energy;
if(vEnergy[i] < energy - energyWidth) break;
else if(vEnergy[i] <= energy + energyWidth)
{
ActMda.vEnergy.push_back(vEnergy[i]);
ActMda.vUncertE.push_back(lines.vUncertE[i]);
ActMda.vYield.push_back(lines.vYield[i]);
ActMda.vUncertY.push_back(lines.vUncertY[i]);
ActMda.fullNames.push_back(lines.fullNames[i]);
ActMda.vPeakIdx.push_back(vPeakIdx[j]+1);
if(lines.key_flag == i) ActMda.key_flag = ActMda.vEnergy.size()-1;
break;
}
}
}
for(size_t i=0; i<ActMda.vYield.size(); ++i)
{
if(ActMda.vYield[i] > maxYield)
{
maxYield = ActMda.vYield[i];
mainPeakIdx = ActMda.vPeakIdx[i]-1;
ActMda.calculateIdx = i;
}
}
if(mainPeakIdx < 0) return;
// 如果该核素属特殊核素,则用“特殊核素配置文件”中指明的其他核素的半衰期
double halflife = (mapHalflife.find(nucName) == mapHalflife.end() ? lines.halflife : mapHalflife[nucName]);
double lambda = log(2.0) / halflife;
CalculateMDCs(phd, ActMda, mainPeakIdx, lambda, maxYield, 1);
ActMda.bCalculateMDA = true;
phd->mapNucActMda[nucName] = ActMda;
}
bool LoadSpectrum(PHDFile *phd, QString fileContents)
{
std::string s = arma_version::as_string();
bool bRet = true;
RadionuclideMessage sample_data;
if(fileContents.isNull())
{
bRet = sample_data.AnalysePHD_File(phd->filepath);
}
else
{
bRet = sample_data.AnalysePHD_Msg(fileContents);
}
if(!bRet) return bRet;
phd->msgInfo = sample_data.GetMessageInfo();
QString block_name = QLatin1String("#Header");
QVariant header = sample_data.GetBlockData(block_name);
if (!header.isNull())
{
phd->header = header.value<HeaderBlock>();
}
block_name = QLatin1String("#Comment");
QVariant comment = sample_data.GetBlockData(block_name);
if(!comment.isNull())
{
phd->oriTotalCmt = comment.value<CommentBlock>();
//phd->totalCmt = phd->oriTotalCmt;
}
block_name = QLatin1String("#Collection");
QVariant collection = sample_data.GetBlockData(block_name);
if (!collection.isNull())
{
phd->collect = collection.value<CollectionBlock>();
if(!phd->collect.collection_start_time.contains('.'))
{
phd->collect.collection_start_time.append(".0");
}
if(!phd->collect.collection_stop_time.contains('.'))
{
phd->collect.collection_stop_time.append(".0");
}
}
else {
phd->collect.air_volume = 0.0;
}
block_name = QLatin1String("#Acquisition");
QVariant acquisition = sample_data.GetBlockData(block_name);
if (!acquisition.isNull())
{
phd->acq = acquisition.value<AcquisitionBlock>();
if(!phd->acq.acquisition_start_time.contains('.'))
{
phd->acq.acquisition_start_time.append(".0");
}
}
else {
phd->acq.acquisition_live_time = 0.0;
phd->acq.acquisition_real_time = 0.0;
}
block_name = QLatin1String("#Processing");
QVariant process = sample_data.GetBlockData(block_name);
if (!process.isNull())
{
phd->process = process.value<ProcessingBlock>();
}
else {
phd->process.sample_volume_of_Xe = 0.0;
phd->process.Xe_collection_yield = 0.0;
phd->process.uncertainty_1 = 0.0;
phd->process.uncertainty_2 = 0.0;
}
block_name = QLatin1String("#Sample");
QVariant sample = sample_data.GetBlockData(block_name);
if(!sample.isNull())
{
phd->sampleBlock = sample.value<SampleBlock>();
}
else {
phd->sampleBlock.dimension_1 = 0.0;
phd->sampleBlock.dimension_2 = 0.0;
}
block_name = QLatin1String("#Calibration");
QVariant calibra = sample_data.GetBlockData(block_name);
if (!calibra.isNull())
{
phd->calibration = calibra.value<CalibrationBlock>();
}
block_name = QLatin1String("#Certificate");
QVariant certify = sample_data.GetBlockData(block_name);
if(!certify.isNull())
{
phd->certificate = certify.value<CertificateBlock>();
}
block_name = QLatin1String("#g_Spectrum");
QVariant var_spec = sample_data.GetBlockData(block_name);
if(!var_spec.isNull())
{
phd->Spec = var_spec.value<G_SpectrumBlock>();
int i=0;
for(; i<phd->Spec.num_g_channel; ++i)
{
if(phd->Spec.counts[i] > 0) break;
}
if(i == phd->Spec.num_g_channel) phd->isValid = false;
}
block_name = QLatin1String("#g_Energy");
QVariant g_ener = sample_data.GetBlockData(block_name);
if(!g_ener.isNull())
{
phd->mapEnerKD[CalPHD] = g_ener.value<G_EnergyBlock>();
}
block_name = QLatin1String("#g_Resolution");
QVariant g_reso = sample_data.GetBlockData(block_name);
if(!g_reso.isNull())
{
phd->mapResoKD[CalPHD] = g_reso.value<G_ResolutionBlock>();
}
block_name = QLatin1String("#g_Efficiency");
QVariant g_effi = sample_data.GetBlockData(block_name);
if(!g_effi.isNull())
{
if(phd->header.system_type == 'C')
{
phd->nMapEffiKD[CalPHD] = g_effi.value<n_G_EfficiencyBlock>();
}
else {
phd->mapEffiKD[CalPHD] = g_effi.value<G_EfficiencyBlock>();
}
}
block_name = QLatin1String("#TotalEff");
QVariant g_tote = sample_data.GetBlockData(block_name);
if(!g_tote.isNull())
{
phd->mapTotEKD[CalPHD] = g_tote.value<TotaleffBlock>();
}
block_name = QLatin1String("#b_self_Attenuation");
QVariant bsel = sample_data.GetBlockData(block_name);
if(!bsel.isNull())
{
phd->bSelfAttenuation = bsel.value<BSelfAttenuationBlock>();
}
block_name = QLatin1String("#b-gEfficiency");
QVariant bgeff = sample_data.GetBlockData(block_name);
if(!bgeff.isNull() && phd->header.system_type != 'C')
{
phd->mapbgEffiKD[CalPHD] = bgeff.value<BG_EfficiencyBlock>();
}
block_name = QLatin1String("#b_Efficiency");
QVariant beff = sample_data.GetBlockData(block_name);
if(!beff.isNull() && phd->header.system_type == 'C')
{
phd->nMapbgEffiKD[CalPHD] = beff.value<NBG_EfficiencyBlock>();
}
// 初始化默认分析设置
if(phd->header.system_type == "P")
{
phd->setting.ECutAnalysis_Low = 35.0;
phd->setting.bUpdateCal = true;
}
phd->setting.refTime_conc = QDateTime::fromString(phd->collect.collection_start_date + " " + phd->collect.collection_start_time, QString(DATATIME_FORMAT));
phd->setting.refTime_act = QDateTime::fromString(phd->acq.acquisition_start_date + " " + phd->acq.acquisition_start_time, QString(DATATIME_FORMAT));
phd->usedSetting = phd->setting;
phd->bAnalyed = false;
phd->analy_start_time = QDateTime::currentDateTimeUtc().toString(DATATIME_FORMAT_SPACE_SECONDS);
qDebug() << QString("Load %1:").arg(phd->filepath) << LINE_END
<< QString("\tChannel number: %1").arg(phd->Spec.num_g_channel) << LINE_END
<< QString("\tStartChannel: %1").arg(phd->Spec.begin_channel) << LINE_END
<< QString("\tEnergy Calibration number: %1").arg(phd->mapEnerKD[CalPHD].record_count) << LINE_END
<< QString("\tResolution Calibration number: %1").arg(phd->mapResoKD[CalPHD].record_count) << LINE_END
<< QString("\tEfficiency Calibration number: %1").arg(phd->mapEffiKD[CalPHD].record_count) << LINE_END
<< QString("\tTotal Efficiency Calibration number: %1").arg(phd->mapTotEKD[CalPHD].record_count) << LINE_END;
return bRet;
}
QStringList UserNuclide(QString sample_type, QString filename)
{
if(filename.isEmpty())
{
if(sample_type.toUpper() == "P")
{
filename = qApp->applicationDirPath() + "/setup/P_default.nuclide";
}
else if(sample_type.toUpper() == "G")
{
filename = qApp->applicationDirPath() + "/setup/G_default.nuclide";
}
}
QStringList nuclideList;
QFile file(filename);
if(file.open(QIODevice::ReadOnly))
{
QTextStream in(&file);
QString nuclide = in.readLine();
while(!nuclide.isNull())
{
nuclideList.push_back(nuclide);
nuclide = in.readLine();
}
file.close();
}
else
{
//QString logName = GammaLogInfo::GetSysLogName(Module_G_Interactive);
//MyLog4qt::MakeNewFileLogger(logName);
//MyLog4qt::SetLoggerText(logName, GammaLogInfo::AddTime("Open default user nuclide library failed When Analyse ") + m_phd->filepath);
}
return nuclideList;
}
void ComputePeakRange(std::vector<PeakInfo> &vPeaks, int n, vec c, vec fwhmch, vec tlo, vec thi)
{
vec fL = 3 + tlo;
vec fH = 3 + thi;
vec l = max(1, floor(c - fL % fwhmch));
vec h = min(n, ceil(c + fH % fwhmch));
ivec Bool = zeros<ivec>(n+1);
for(uword i=0; i<l.size(); ++i)
{
Bool.rows(l(i)-1, h(i)-1) = ones<ivec>(h(i)-l(i)+1);
}
Bool(0) = 0;
Bool(n) = 0;
ivec d = diff(Bool);
uvec rg_low = find(d == 1);
uvec rg_high = find(d == -1);
// 求出每个峰的左右边界及重峰信息
uvec idx;
int multiPeaks = 0, L, R, size_idx;
for(uword i=0; i<rg_low.n_rows; ++i)
{
L = rg_low(i) + 1; // Ranges返回的是数组下标但道从1开始
R = rg_high(i) + 1;
idx = find(c >= L && c <= R);
size_idx = idx.size();
if(size_idx > 1)
{
++multiPeaks;
for(uword j=0; j<size_idx; ++j)
{
vPeaks[idx(j)].left = L;
vPeaks[idx(j)].right = R;
vPeaks[idx(j)].multiIndex = multiPeaks;
//vPeaks[idx(j)].no_peaks = j+1;
//vPeaks[idx(j)].num_peaks = size_idx;
}
} else {
vPeaks[idx(0)].left = L;
vPeaks[idx(0)].right = R;
//vPeaks[idx(0)].no_peaks = 0;
//vPeaks[idx(0)].num_peaks = 0;
vPeaks[idx(0)].multiIndex = 0;
}
}
}
stdvec energyToChannel(rowvec e, vec p)
{
// maximal iterations in newton search
int maxiter = 100;
// relative error tolerance
double tol = 1e-6;
// get linear approximation
rowvec x_test("0, 8000"); // x_test = [0, 8000];
rowvec y_test = Independ::calFcnEval(x_test, p);
double d = y_test(0);
rowvec k = diff(y_test) / diff(x_test);
// y_test.print("y_test = ");
// k.print("k = ");
// find bad energies
uvec k2 = find_nonfinite(e);
k2.insert_rows(k2.n_elem, find(e<1));
// calculate with 1 and assing NaN afterwards
e(k2).fill(1);
// first estimate for channel
rowvec c = (e - d) / k;
rowvec eApp = Independ::calFcnEval(c, p);
uvec idx = find(abs(eApp - e) > tol*e);
//idx.print("idx = ");
int iter = 1;
while (iter < maxiter && !idx.is_empty())
{
iter = iter + 1;
// newton step
rowvec dE = Independ::calDerivEval(c(idx), p);
// c.print("c = ");
// e.print("e = ");
// eApp.print("eApp = ");
// dE.print("dE = ");
c(idx) = c(idx) + (e(idx) - eApp(idx)) / dE;
// check converged
eApp(idx) = Independ::calFcnEval(c(idx), p);
uvec subidx = find( abs(eApp(idx) - e(idx)) > tol*e(idx) );
idx = idx(subidx);
}
if (!idx.is_empty())
{
qDebug() << QString("Could not find Channel for %1 energies").arg(idx.n_elem);
c(idx).fill(gNaN);
}
// bad Energies
c(k2).fill(gNaN);
return conv_to<stdvec>::from(c);
}
void UpdateEfficiency(PHDFile *phd)
{
size_t peakNum = phd->vPeak.size();
stdvec vEner(peakNum);
for (size_t i=0; i<peakNum; ++i)
{
vEner[i] = phd->vPeak[i].energy;
}
stdvec vEffi = calValuesOut(vEner, phd->usedEffiPara.p);
for (size_t i=0; i<peakNum; ++i)
{
phd->vPeak[i].efficiency = vEffi[i];
}
}
bool WriteBaseInfo(BaseControls &baseCtrl, const QString &filename)
{
QFile file(filename);
if(!file.open(QIODevice::WriteOnly | QIODevice::Text)) return false;
int numPerLine = 5;
QTextStream out(&file);
out.setFieldWidth(15);
out.setFieldAlignment(QTextStream::AlignLeft);
out.setPadChar(' ');
out << "#AnalyseRange" << LINE_END;
out << baseCtrl.rg_low << baseCtrl.rg_high << LINE_END;
size_t i, nCP = baseCtrl.XCtrl.size(), nGroupCP = nCP / numPerLine * numPerLine;
out << "#XCtrl" << LINE_END;
stdvec &cx = baseCtrl.XCtrl;
for(i=0; i<nGroupCP; i+=numPerLine)
{
out << cx[i] << cx[i+1] << cx[i+2] << cx[i+3] << cx[i+4] << LINE_END;
}
if(i < nCP)
{
for(; i<nCP; ++i) out << cx[i];
out << LINE_END;
}
out << "#YCtrl" << LINE_END;
stdvec &cy = baseCtrl.YCtrl;
for(i=0; i<nGroupCP; i+=numPerLine)
{
out << cy[i] << cy[i+1] << cy[i+2] << cy[i+3] << cy[i+4] << LINE_END;
}
if(i < nCP)
{
for(; i<nCP; ++i) out << cy[i];
out << LINE_END;
}
out << "#YSlope" << LINE_END;
stdvec &cdy = baseCtrl.YSlope;
for(i=0; i<nGroupCP; i+=numPerLine)
{
out << cdy[i] << cdy[i+1] << cdy[i+2] << cdy[i+3] << cdy[i+4] << LINE_END;
}
if(i < nCP)
{
for(; i<nCP; ++i) out << cdy[i];
out << LINE_END;
}
size_t nBL = baseCtrl.Baseline.size(), nGroupBL = nBL / numPerLine * numPerLine;
out << "#Baseline" << LINE_END;
out << nBL << LINE_END;
stdvec &bl = baseCtrl.Baseline;
for(i=0; i<nGroupBL; i+=numPerLine)
{
out << bl[i] << bl[i+1] << bl[i+2] << bl[i+3] << bl[i+4] << LINE_END;
}
if(i < nBL)
{
for(; i<nBL; ++i) out << bl[i];
out << LINE_END;
}
out << "#StepCounts" << LINE_END;
out << nBL << LINE_END;
stdvec &sc = baseCtrl.StepCounts;
for(i=0; i<nGroupBL; i+=numPerLine)
{
out << sc[i] << sc[i+1] << sc[i+2] << sc[i+3] << sc[i+4] << LINE_END;
}
if(i < nBL)
{
for(; i<nBL; ++i) out << sc[i];
out << LINE_END;
}
//Add By XYL 20170927
file.close();
return true;
}
bool WriteLcScac(const stdvec &vData, QString type, const QString &filename)
{
QFile file(filename);
if(!file.open(QIODevice::WriteOnly | QIODevice::Text)) return false;
int numPerLine = 5;
QTextStream out(&file);
out.setFieldWidth(15);
out.setFieldAlignment(QTextStream::AlignLeft);
out.setPadChar(' ');
out << QString("#%1").arg(type) << LINE_END;
out << vData.size() << LINE_END;
size_t i, n = vData.size(), nGroup = n / numPerLine * numPerLine;
for(i=0; i<nGroup; i+=numPerLine)
{
out << vData[i] << vData[i+1] << vData[i+2] << vData[i+3] << vData[i+4] << LINE_END;
}
if(i < n)
{
for(; i<n; ++i) out << vData[i];
out << LINE_END;
}
//Add By XYL 20170927
file.close();
return true;
}
QString EquationDescription(int funcId)
{
QString desc;
switch (funcId)
{
case 1: desc = "y=yi+(y(i+1)-yi)*(x-xi)/(x(i+1)-xi) for xi<=x<x(i+1)"; break;
case 2: desc = "y=a0+a1*x+a2*x^2+a3*x^3"; break;
case 3: desc = "y=a0+a1*x^(1/2)+a2*x+a3*x^(3/2)"; break;
case 4: desc = "y=sqrt(a0+a1*x+a2*x^2+a3*x^3)"; break;
case 5: desc = "y=A*exp(-(E1/x)^k)*(1-exp(-(E2/x)^n))"; break;
case 6: desc = "log(y)=a0+a1*log(x)+a2*log(x)^2+a3*log(x)^3"; break;
case 7: desc = "log(y)=a0*x+a1+a2/x+a3/(x^2)"; break;
case 8: desc = "log(y)=a0+a1*log(c/x)+a2*log(c/x)^2+a3*log(c/x)^3+a4*log(c/x)^4"; break;
case 9: desc = "y=1/(a*x^(-k)+b*x^(-n))"; break;
case 93: desc = "y=S*exp(-(E1/x)^k)*(1-exp(-(2*E3/(x-E3))^n"; break;
case 94: desc = "y=S*exp(-(E1/x)^k)*(1-exp(-b*(1/(x-E2))^m))"; break;
case 95: desc = "y=S*exp(-(E1/x)^k)*(1-exp(-b*(1/(x-E2))^m))*(1-exp(-(2*E3/(E-E3))^n))"; break;
case 96: desc = "y=A+E1*x^(-n)"; break;
case 97: desc = "y=a0+a1*exp(-lam1*x)"; break;
case 98: desc = "y=c"; break;
case 99: desc = "y=t_0+k_t*x if x>ECut"; break;
default: desc = ""; break;
}
return desc;
}
QString EquationName(int funcId)
{
QString name;
switch (funcId)
{
case 1: name = "Interpolation"; break;
case 2: name = "Polynomial"; break;
case 3: name = "Square root polynomial"; break;
case 4: name = "Square root of polynomial"; break;
case 5: name = "HT Efficiency"; break;
case 6: name = "Polynomial in log(y) against log(x)"; break;
case 7: name = "Polynomial in log(y)"; break;
case 8: name = "Polynomial in log(y) against log(1/x)"; break;
case 9: name = "Inverse exponential"; break;
case 93: name = "HAE Efficiency (1-3)"; break;
case 94: name = "HAE Efficiency (1-2)"; break;
case 95: name = "HAE Efficiency (1-2-3)"; break;
case 96: name = "Inverse power"; break;
case 97: name = "Exponential sum"; break;
case 98: name = "Constant"; break;
case 99: name = "Linear with cut-off"; break;
default: name = ""; break;
}
return name;
}
}
void PrintMat2fvec(QString name, field<vec>& vecc)
{
QFile file("C:/Result_Debug3/" + name);
if (file.open(QIODevice::WriteOnly))
{
QTextStream out(&file);
for (uword row = 0; row < vecc.size(); ++row)
{
for (uword col = 0; col < vecc(row).size(); ++col)
{
out << QString::number(vecc(row)(col), 'f', 20) << "\t";
}
out << "\r\n";
}
file.close();
}
else
{
}
}
void PrintMat22(QString name, mat& dataMat)
{
QFile file("C:/Result_Debug3/" + name);
if (file.open(QIODevice::WriteOnly))
{
QTextStream out(&file);
if (dataMat.is_vec())
{
for (int i = 0; i < dataMat.size(); ++i)
{
out << QString::number(dataMat(i), 'g', 10) << "\r\n";
}
}
else {
for (uword row = 0; row < dataMat.n_rows; ++row)
{
for (uword col = 0; col < dataMat.n_cols; ++col)
{
out << QString::number(dataMat(row, col), 'g', 10) << "\t";
}
out << "\r\n";
}
}
file.close();
}
}
void PrintColvec(QString name, colvec& dataMat)
{
QFile file("C:/Result_Debug3/" + name);
if (file.open(QIODevice::WriteOnly))
{
QTextStream out(&file);
for (int i = 0; i < dataMat.size(); ++i)
{
out << QString::number(dataMat(i), 'g', 10) << "\r\n";
}
file.close();
}
}
void PrintFieldVec(QString name, field<vec>& vecc)
{
QFile file(QCoreApplication::applicationDirPath() + "/Result_Debug3/" + name);
if (file.open(QIODevice::WriteOnly))
{
QTextStream out(&file);
for (uword row = 0; row < vecc.size(); ++row)
{
for (uword col = 0; col < vecc(row).size(); ++col)
{
out << QString::number(vecc(row)(col), 'g', 20) << "\t";
}
out << "\r\n";
}
file.close();
}
}
/* ***********************************************************导出的类成员函数实现************************************************************ */
GammaAnalyALG::GammaAnalyALG()
{
Init();
}
GammaAnalyALG::GammaAnalyALG(SpecSetup setup)
{
Init();
specSetup = setup;
}
GammaAnalyALG::GammaAnalyALG(PHDFile *phd)
{
specSetup = phd->setting;
SetBaseInfo(phd->Spec.counts, phd->Spec.begin_channel);
QStringList names;
names << phd->newEner << phd->newReso;
SetCalPara(names, phd->mapEnerPara[names[0]], phd->mapResoPara[names[1]], phd->mapEffiPara[phd->newEffi], phd->mapTotEPara[phd->newTotE]);
SetPAT(phd->vPeak);
baseInfo(s_AnalysisRange) << phd->baseCtrls.rg_low << phd->baseCtrls.rg_high;
baseInfo(s_XControl) = phd->baseCtrls.XCtrl;
baseInfo(s_YControl) = phd->baseCtrls.YCtrl;
baseInfo(s_YSlope) = phd->baseCtrls.YSlope;
baseInfo(s_BaseLine) = phd->baseCtrls.Baseline;
baseInfo(s_Steps) = phd->baseCtrls.StepCounts;
}
GammaAnalyALG::~GammaAnalyALG()
{
}
void GammaAnalyALG::Init()
{
setJniBl(NULL, NULL, "");
m_strFilePath = "";
m_curEner = CalPHD;
m_curReso = CalPHD;
baseInfo.set_size(MAX_INDEX_SPEC);
field<vec> f(2);
f(0) << 2 << 0 << 0.3 << 0;
f(1) = gNaN * ones<vec>(f(0).size()-1u);
Acal_energy_para[CalPHD] = f;
f(0) << 4 << 1 << 0.003 << 0;
f(1) = gNaN * ones<vec>(f(0).size()-1u);
Acal_resolution_para[CalPHD] = f;
f(0) << 94 << 0.1 << 30 << 3 << 50 << 70 << 0.8;
f(1) = gNaN * ones<vec>(f(0).size()-1u);
Acal_efficiency_para = f;
Acal_tot_efficiency_para = f;
Acal_tail_para << 99 << 0 << 0 << 0;
Acal_tail_alpha_para << 98 << 1;
Acal_tail_right_para << 98 << 0.35;
Acal_tail_right_alpha_para << 98 << 1;
Acal_step_ratio_para << 98 << 0;
}
bool GammaAnalyALG::Process(bool bUpdate, QChar dataType, vec certEne)
{
if( !any(baseInfo(s_Spectrum) > 0) )
{
qDebug() << "The counts of spectrum are all zero!";
return false;
}
callGammaProcess(0);
stdvec vy;
vy = getBaseInfo(s_YControl);
/* **********************************刻度更新*********************************************** */
if(bUpdate)
{
calUpdate(dataType, certEne, true, true, true);
callGammaProcess(1);
}
vy = getBaseInfo(s_YControl);
/* **********************************寻峰*********************************************** */
PeakSearch();
callGammaProcess(2);
////("____ps_peaks.txt", pat.Peaks);
vy = getBaseInfo(s_YControl);
/* **********************************基线优化*********************************************** */
baseImprove(specSetup.BaseImprovePSS);
callGammaProcess(3);
////("____bi_peaks.txt", pat.Peaks);
vy = getBaseInfo(s_YControl);
/* **********************************净面积拟合*********************************************** */
vec t_idx = dmps(QStringList("NetAreaFree"))(0);
//t_idx.print("t_idx");
uvec idx = find(t_idx), ep;
vy = getBaseInfo(s_YControl);
//idx.print("idx = ");
fitPeakFull(idx, ep, ep, ep);
vy = getBaseInfo(s_YControl);
callGammaProcess(4);
baseInfo(s_Energy) = calValues(Cal_Energy, RangeVec2(1, m_nChans));
baseInfo(s_FwhmcAll) = GetFwhmcAll();
baseInfo(s_Lc) = calculateLC(baseInfo(s_BaseLine), baseInfo(s_FwhmcAll), specSetup.RiskLevelK);
baseInfo(s_Scac) = calculateSCAC(baseInfo(s_Spectrum), baseInfo(s_BaseLine), baseInfo(s_FwhmcAll));
vy = getBaseInfo(s_YControl);
/*//("____na_peaks.txt", pat.Peaks);
////("____na_bl.txt", baseInfo(s_BaseLine));
////("____na_cx.txt", baseInfo(s_XControl));
////("____na_cy.txt", baseInfo(s_YControl));
////("____na_cdy.txt", baseInfo(s_YSlope));
////("____na_steps.txt", baseInfo(s_Steps));*/
return true;
}
void GammaAnalyALG::callGammaProcess(int npro)
{
if (m_pJniEnv == NULL || m_pSendObj == NULL)
return;
jclass jcls = m_pJniEnv->GetObjectClass(m_pSendObj);
if (jcls == NULL)
{
return;
}
jmethodID jmth = m_pJniEnv->GetMethodID(jcls, "gammaProcess", "(Ljava/lang/String;Ljava/lang/String;)V");
if (jmth == NULL)
{
return;
}
jobject obj = m_pJniEnv->AllocObject(jcls);
if (obj == NULL)
{
return;
}
jstring juser = m_pJniEnv->NewStringUTF(m_strUserId.toStdString().c_str());
char szbuf[32] = { 0 };
itoa(npro, szbuf, 10);
jstring jpro = m_pJniEnv->NewStringUTF(szbuf);
m_pJniEnv->CallVoidMethod(obj, jmth, juser, jpro);
}
bool GammaAnalyALG::AnalyseSpectrum(PHDFile *phd, QMap<QString, NuclideLines> mapLines)
{
// 分析处理调用
//GammaAnalyALG alg(phd->setting);
specSetup = phd->setting;
SetBaseInfo(phd->Spec.counts, phd->Spec.begin_channel);
QStringList names;
names << phd->newEner << phd->newReso;
bool bUpdate = false;
if(phd->bAnalyed)
{
SetCalPara(names, phd->mapEnerPara[ names[0] ], phd->mapResoPara[ names[1] ], phd->mapEffiPara[phd->newEffi], phd->mapTotEPara[phd->newTotE]);
}
else {
SetCalData(names, phd->mapEnerKD[ names[0] ], phd->mapResoKD[ names[1] ], phd->mapEffiKD[phd->newEffi], phd->mapTotEKD[phd->newTotE]);
if(phd->setting.bUpdateCal) bUpdate = true;
}
if(!Process(bUpdate, phd->msgInfo.data_type.at(0), phd->certificate.g_energy.toStdVector())) return false;
// 获取分析结果
if(!phd->bAnalyed)
{
phd->bAnalyed = true;
//?
if(!phd->mapEnerKD.isEmpty()) phd->mapEnerPara[phd->newEner] = GetCalPara(Cal_Energy, phd->newEner);
if(!phd->mapResoKD.isEmpty()) phd->mapResoPara[phd->newReso] = GetCalPara(Cal_Resolution, phd->newReso);
if(!phd->mapEffiKD.isEmpty()) phd->mapEffiPara[phd->newEffi] = GetCalPara(Cal_Efficiency);
if(!phd->mapTotEKD.isEmpty()) phd->mapTotEPara[phd->newTotE] = GetCalPara(Cal_Tot_efficiency);
phd->para_stepRatio = GetCalPara(Cal_Step_ratio);
phd->para_tail = GetCalPara(Cal_Tail);
phd->para_tailAlpha = GetCalPara(Cal_Tail_alpha);
phd->para_tailRight = GetCalPara(Cal_Tail_right);
phd->para_tailRightAlpha = GetCalPara(Cal_Tail_right_alpha);
stdvec2 datas = GetCalData(Cal_Energy, CalUpdate1);
if(datas.size() == 3)
{
G_EnergyBlock g_ener;
g_ener.centroid_channel = QVector<double>::fromStdVector(datas[0]);
g_ener.g_energy = QVector<double>::fromStdVector(datas[1]);
g_ener.uncertainty = QVector<double>::fromStdVector(datas[2]);
g_ener.record_count = g_ener.centroid_channel.size();
phd->newEner = CalUpdate1;
phd->mapEnerKD[CalUpdate1] = g_ener;
phd->mapEnerPara[CalUpdate1] = GetCalPara(Cal_Energy, CalUpdate1);
}
datas = GetCalData(Cal_Energy, CalUpdate2);
if(datas.size() == 3)
{
G_EnergyBlock g_ener;
g_ener.centroid_channel = QVector<double>::fromStdVector(datas[0]);
g_ener.g_energy = QVector<double>::fromStdVector(datas[1]);
g_ener.uncertainty = QVector<double>::fromStdVector(datas[2]);
g_ener.record_count = g_ener.centroid_channel.size();
phd->newEner = CalUpdate2;
phd->mapEnerKD[CalUpdate2] = g_ener;
phd->mapEnerPara[CalUpdate2] = GetCalPara(Cal_Energy, CalUpdate2);
}
datas = GetCalData(Cal_Resolution, CalResUpdate);
if(datas.size() == 3)
{
G_ResolutionBlock g_reso;
g_reso.g_energy = QVector<double>::fromStdVector(datas[0]);
g_reso.FWHM = QVector<double>::fromStdVector(datas[1]);
g_reso.uncertainty = QVector<double>::fromStdVector(datas[2]);
g_reso.record_count = g_reso.FWHM.size();
phd->newReso = CalResUpdate;
phd->mapResoKD[CalResUpdate] = g_reso;
phd->mapResoPara[CalResUpdate] = GetCalPara(Cal_Resolution, CalResUpdate);
}
}
phd->vEnergy = getBaseInfo(s_Energy);
phd->vBase = getBaseInfo(s_BaseLine);
phd->vLc = getBaseInfo(s_Lc);
phd->vScac = getBaseInfo(s_Scac);
phd->vPeak = getPAT();
stdvec AnalyRg = getBaseInfo(s_AnalysisRange);
if(AnalyRg.size() > 1)
{
phd->baseCtrls.rg_low = AnalyRg[0];
phd->baseCtrls.rg_high = AnalyRg[1];
}
else {
phd->baseCtrls.rg_low = 1;
phd->baseCtrls.rg_high = phd->Spec.num_g_channel;
}
phd->baseCtrls.Baseline = phd->vBase;
phd->baseCtrls.XCtrl = getBaseInfo(s_XControl);
phd->baseCtrls.YCtrl = getBaseInfo(s_YControl);
phd->baseCtrls.YSlope = getBaseInfo(s_YSlope);
phd->baseCtrls.StepCounts = getBaseInfo(s_Steps);
BaseCtrlStack bcStack;
bcStack.cx = phd->baseCtrls.XCtrl;
bcStack.cy = phd->baseCtrls.YCtrl;
bcStack.cdy = phd->baseCtrls.YSlope;
phd->baseCtrls.BaseStack.clear();
phd->baseCtrls.BaseStack.push_back(bcStack);
// 保存当前分析所用的参数信息
phd->usedSetting = phd->setting;
if(!phd->mapEnerKD.isEmpty())
{
phd->usedEner = phd->newEner;
phd->usedEnerKD = phd->mapEnerKD[phd->newEner];
phd->usedEnerPara = phd->mapEnerPara[phd->newEner];
}
if(!phd->mapResoKD.isEmpty())
{
phd->usedReso = phd->newReso;
phd->usedResoKD = phd->mapResoKD[phd->newReso];
phd->usedResoPara = phd->mapResoPara[phd->newReso];
}
if(!phd->mapEffiKD.isEmpty())
{
phd->usedEffi = phd->newEffi;
phd->usedEffiKD = phd->mapEffiKD[phd->newEffi];
phd->usedEffiPara = phd->mapEffiPara[phd->newEffi];
}
if(!phd->mapTotEKD.isEmpty())
{
phd->usedTotE = phd->newTotE;
phd->usedTotEKD = phd->mapTotEKD[phd->newTotE];
phd->usedTotEPara = phd->mapTotEPara[phd->newTotE];
}
callGammaProcess(5);
AlgFunc::NuclidesIdent(phd, mapLines, this->m_strFilePath);
AlgFunc::RunQC(phd, m_strFilePath);
callGammaProcess(6);
return true;
}
void GammaAnalyALG::SetSpecSetup(SpecSetup setup)
{
specSetup = setup;
}
void GammaAnalyALG::SetBaseInfo(stdvec spec)
{
m_nChans = spec.size();
baseInfo(s_Spectrum) = spec;
baseInfo(s_Approximation) = zeros<rowvec>(m_nChans);
baseInfo(s_BaseLine) = zeros<rowvec>(m_nChans);
}
void GammaAnalyALG::SetBaseInfo(const QVector<long> &spec, int begin_channel)
{
m_nChans = spec.size();
baseInfo(s_Spectrum) = zeros<rowvec>(m_nChans);
baseInfo(s_Approximation) = zeros<rowvec>(m_nChans);
baseInfo(s_BaseLine) = zeros<rowvec>(m_nChans);
for(int i=0; i<m_nChans; ++i)
{
baseInfo(s_Spectrum)[i] = spec[i];
}
if(begin_channel == 0)
{
baseInfo(s_Spectrum).insert_cols(m_nChans, 1);
baseInfo(s_Spectrum)(m_nChans) = gNaN;
baseInfo(s_Spectrum).shed_col(0);
}
}
void GammaAnalyALG::SetBaseInfo(stdvec2 info)
{
if(info.size() != MAX_INDEX_SPEC)
{
qDebug() << QString("The size of BaseInfo must be %1.").arg(MAX_INDEX_SPEC);
return;
}
if(baseInfo.size() < MAX_INDEX_SPEC) baseInfo.set_size(MAX_INDEX_SPEC);
for(int i=0; i<MAX_INDEX_SPEC; ++i)
{
baseInfo(i) = rowvec( info[i] );
}
m_nChans = baseInfo(s_Spectrum).size();
/*baseInfo(s_BaseLine) = info.BaseLine;
baseInfo(s_AnalysisRange) = info.AnalysisRange;
baseInfo(s_XControl) = info.XControl;
baseInfo(s_YControl) = info.YControl;
baseInfo(s_YSlope) = info.YSlope;
baseInfo(s_Residual) = info.Residual;
baseInfo(s_Stripped) = info.Stripped;
baseInfo(s_Steps) = info.Steps;
baseInfo(s_ROISmooth) = info.ROISmooth;
baseInfo(s_SmoothStripped) = info.SmoothStripped;
baseInfo(s_Approximation) = info.Approx;
baseInfo(s_ChiSquare) = info.ChiSquare;
baseInfo(s_ArtXControl) = info.ArtXControl;
baseInfo(s_ArtYControl) = info.ArtYControl;
baseInfo(s_ArtYSlope) = info.ArtYSlope;
baseInfo(s_ArtificialBase) = info.ArtificialBase;
baseInfo(s_BaseFittingWeights) = info.BaseFittingWeights;
baseInfo(s_PSS) = info.PSS;*/
}
stdvec2 GammaAnalyALG::getBaseInfo()
{
stdvec2 info;
for(uword i=0; i<baseInfo.size(); ++i)
{
info.push_back( conv_to< stdvec >::from( baseInfo(i) ) );
}
return info;
/*AlgBaseInfo info;
info.BaseLine = conv_to< vector<double> >::from( baseInfo(s_BaseLine) );
info.AnalysisRange = conv_to< vector<double> >::from( baseInfo(s_AnalysisRange) );
info.XControl = conv_to< vector<double> >::from( baseInfo(s_XControl) );
info.YControl = conv_to< vector<double> >::from( baseInfo(s_YControl) );
info.YSlope = conv_to< vector<double> >::from( baseInfo(s_YSlope) );
info.Residual = conv_to< vector<double> >::from( baseInfo(s_Residual) );
info.Stripped = conv_to< vector<double> >::from( baseInfo(s_Stripped) );
info.Steps = conv_to< vector<double> >::from( baseInfo(s_Steps) );
info.ROISmooth = conv_to< vector<double> >::from( baseInfo(s_ROISmooth) );
info.SmoothStripped = conv_to< vector<double> >::from( baseInfo(s_SmoothStripped) );
info.Approx = conv_to< vector<double> >::from( baseInfo(s_Approximation) );
info.ChiSquare = conv_to< vector<double> >::from( baseInfo(s_ChiSquare) );
info.ArtXControl = conv_to< vector<double> >::from( baseInfo(s_ArtXControl) );
info.ArtYControl = conv_to< vector<double> >::from( baseInfo(s_ArtYControl) );
info.ArtYSlope = conv_to< vector<double> >::from( baseInfo(s_ArtYSlope) );
info.ArtificialBase = conv_to< vector<double> >::from( baseInfo(s_ArtificialBase) );
info.BaseFittingWeights = conv_to< vector<double> >::from( baseInfo(s_BaseFittingWeights) );
info.PSS = conv_to< vector<double> >::from( baseInfo(s_PSS) );
return info;*/
}
stdvec GammaAnalyALG::getBaseInfo(SpecBaseIdx idx)
{
return conv_to<stdvec>::from(baseInfo(idx));
}
void GammaAnalyALG::SetPAT(const std::vector<PeakInfo> &vPeaks)
{
pat.Peaks.set_size(vPeaks.size(), MAX_INDEX_PEAK);
pat.Peaks.fill(gNaN);
uword i_row = 0;
foreach(PeakInfo peak, vPeaks)
{
pat.Peaks(i_row, i_Centroid) = peak.peakCentroid;
pat.Peaks(i_row, i_Energy) = peak.energy;
pat.Peaks(i_row, i_FWHM) = peak.fwhm;
pat.Peaks(i_row, i_FWHM_Ch) = peak.fwhmc;
pat.Peaks(i_row, i_NetArea) = peak.area;
pat.Peaks(i_row, i_AreaError) = peak.areaErr;
pat.Peaks(i_row, i_Efficiency) = peak.efficiency;
pat.Peaks(i_row, i_Significance) = peak.significance;
pat.Peaks(i_row, i_Sensitivity) = peak.sensitivity;
pat.Peaks(i_row, i_StepRatio) = peak.stepRatio;
pat.Peaks(i_row, i_Tail) = peak.tail;
pat.Peaks(i_row, i_TailAlpha) = peak.tailAlpha;
pat.Peaks(i_row, i_UpperTail) = peak.upperTail;
pat.Peaks(i_row, i_UpperTailAlpha) = peak.upperTailAlpha;
pat.Peaks(i_row, i_LC) = peak.lc;
pat.Peaks(i_row, i_LD) = peak.ld;
pat.Peaks(i_row, i_MeanBackCount) = peak.meanBackCount;
pat.Peaks(i_row, i_BackgroundArea) = peak.backgroundArea;
++i_row;
}
}
std::vector<PeakInfo> GammaAnalyALG::getPAT()
{
// 向 PAT 表中存在但尚未写入值的项写入值(这里只写入需要的)
dmps();
// pat表中没有的项只能通过字符串获取
QStringList items;
items << "BWWidthChan" << "RecoilBetaChan" << "RecoilDeltaChan";
field<vec> f = dmps(items);
std::vector<PeakInfo> vPI;
for(uword row=0; row<pat.Peaks.n_rows; ++row)
{
PeakInfo pi;
pi.index = row + 1;
pi.peakCentroid = pat.Peaks(row, i_Centroid);
pi.energy = pat.Peaks(row, i_Energy);
pi.fwhm = pat.Peaks(row, i_FWHM);
pi.fwhmc = pat.Peaks(row, i_FWHM_Ch);
pi.area = pat.Peaks(row, i_NetArea);
pi.areaErr = pat.Peaks(row, i_AreaError);
pi.significance = pat.Peaks(row, i_Significance);
pi.sensitivity = pat.Peaks(row, i_Sensitivity);
pi.multiIndex = pat.Peaks(row, i_Multiplet);
pi.efficiency = pat.Peaks(row, i_Efficiency);
pi.stepRatio = pat.Peaks(row, i_StepRatio);
pi.tail = pat.Peaks(row, i_Tail);
pi.tailAlpha = pat.Peaks(row, i_TailAlpha);
pi.upperTail = pat.Peaks(row, i_UpperTail);
pi.upperTailAlpha = pat.Peaks(row, i_UpperTailAlpha);
pi.lc = pat.Peaks(row, i_LC);
pi.ld = pat.Peaks(row, i_LD);
pi.meanBackCount = pat.Peaks(row, i_MeanBackCount);
pi.backgroundArea = pat.Peaks(row, i_BackgroundArea);
pi.BWWidthChan = f(0)(row);
pi.recoilBetaChan = f(1)(row);
pi.recoilDeltaChan = f(2)(row);
vPI.push_back(pi);
}
// 求出每个峰的左右边界及重峰信息
uvec idx;
umat Ranges = lGetPlotRanges(m_nChans, pat.Peaks.col(i_Centroid), pat.Peaks.col(i_FWHM_Ch), pat.Peaks.col(i_Tail), pat.Peaks.col(i_UpperTail));
int multiPeaks = 0, L, R, size_idx;
for(uword i=0; i<Ranges.n_rows; ++i)
{
L = Ranges(i, 0) + 1; // Ranges返回的是数组下标但道从1开始
R = Ranges(i, 1) + 1;
idx = find(pat.Peaks.col(i_Centroid) >= L && pat.Peaks.col(i_Centroid) <= R);
size_idx = idx.size();
if(size_idx > 1)
{
++multiPeaks;
for(uword j=0; j<size_idx; ++j)
{
vPI[idx(j)].left = L;
vPI[idx(j)].right = R;
vPI[idx(j)].multiIndex = multiPeaks;
//vPI[idx(j)].no_peaks = j+1;
//vPI[idx(j)].num_peaks = size_idx;
}
} else {
vPI[idx(0)].left = L;
vPI[idx(0)].right = R;
vPI[idx(0)].multiIndex = 0;
//vPI[idx(0)].num_peaks = 0;
}
}
return vPI;
}
stdvec GammaAnalyALG::getPatInfo(PeakIndex idx)
{
if( !pat.Peaks.col(idx).is_finite() )
{
return conv_to<stdvec>::from( dmps(idx) );
}
else return conv_to<stdvec>::from( pat.Peaks.col(idx) );
}
void GammaAnalyALG::setFilePath(QString strPath)
{
m_strFilePath = strPath;
}
void GammaAnalyALG::setJniBl(JNIEnv* pEnv, jobject jobj, QString strUid)
{
m_pJniEnv = pEnv;
m_pSendObj = jobj;
m_strUserId = strUid;
}
void GammaAnalyALG::SetCalData(CalibrationType cal, stdvec x, stdvec y, stdvec err, QString calName)
{
int s = x.size();
if(s != y.size())
{
qDebug() << "The size of x must be same as the size of y in function SetCalData.";
return;
}
if(err.size() != s) { err.assign(s, gNaN); }
mat calDatas(s, 3);
calDatas.col(0) = vec(x);
calDatas.col(1) = vec(y);
calDatas.col(2) = vec(err);
switch (cal)
{
case Cal_Energy: Cal_Ener_Data[calName] = calDatas; break;
case Cal_Resolution: Cal_Reso_Data[calName] = calDatas; break;
case Cal_Efficiency: Cal_Effi_Data = calDatas; break;
case Cal_Tot_efficiency: Cal_Tot_Effi_Data = calDatas; break;
default: break;
}
}
void GammaAnalyALG::SetCalData(QStringList names, G_EnergyBlock EnerKD, G_ResolutionBlock ResoKD, G_EfficiencyBlock EffiKD, TotaleffBlock TotEKD)
{
if(names.size() < 2)
{
qDebug() << "The parameter names must has at least 2 element in GammaAnalyALG::SetCalData.";
return;
}
m_curEner = names[0];
m_curReso = names[1];
if(EnerKD.record_count > 0)
{
mat calDatas(EnerKD.record_count, 3);
calDatas.col(0) = vec( EnerKD.centroid_channel.toStdVector() );
calDatas.col(1) = vec( EnerKD.g_energy.toStdVector() );
calDatas.col(2) = vec( EnerKD.uncertainty.toStdVector() );
Cal_Ener_Data[m_curEner] = calDatas;
Acal_energy_para[m_curEner] = FitCalPara(Cal_Energy, calDatas);
//Acal_energy_para[m_curEner].print("Acal_energy_para = ");
}
if(ResoKD.record_count > 0)
{
mat calDatas(ResoKD.record_count, 3);
calDatas.col(0) = vec( ResoKD.g_energy.toStdVector() );
calDatas.col(1) = vec( ResoKD.FWHM.toStdVector() );
calDatas.col(2) = vec( ResoKD.uncertainty.toStdVector() );
Cal_Reso_Data[m_curReso] = calDatas;
Acal_resolution_para[m_curReso] = FitCalPara(Cal_Resolution, calDatas);
PrintFieldVec("Acal_resolution_para[m_curReso].txt", Acal_resolution_para[m_curReso]);
//Acal_resolution_para[m_curReso].print("Acal_resolution_para = ");
}
if(EffiKD.record_count > 0)
{
mat calDatas(EffiKD.record_count, 3);
calDatas.col(0) = vec( EffiKD.g_energy.toStdVector() );
calDatas.col(1) = vec( EffiKD.efficiency.toStdVector() );
calDatas.col(2) = vec( EffiKD.uncertainty.toStdVector() );
Cal_Effi_Data = calDatas;
Acal_efficiency_para = FitCalPara(Cal_Efficiency, Cal_Effi_Data);
//Acal_efficiency_para.print("Acal_efficiency_para = ");
}
if(TotEKD.record_count > 0)
{
mat calDatas(TotEKD.record_count, 3);
calDatas.col(0) = vec( TotEKD.g_energy.toStdVector() );
calDatas.col(1) = vec( TotEKD.total_efficiency.toStdVector() );
calDatas.col(2) = vec( TotEKD.uncertainty.toStdVector() );
Cal_Tot_Effi_Data = calDatas;
Acal_tot_efficiency_para = FitCalPara(Cal_Tot_efficiency, Cal_Tot_Effi_Data);
//Acal_tot_efficiency_para.print("Acal_tot_efficiency_para = ");
}
/*QString name;
setCalibrationData(name, Cal_Energy, "PHD", tmpEner);
setCalibrationData(name, Cal_Resolution, "PHD", QVec2ToMat(ResoKD));
setCalibrationData(name, Cal_Efficiency, "PHD", QVec2ToMat(EffiKD));
setCalibrationData(name, Cal_Tot_efficiency, "PHD", QVec2ToMat(TotEKD)); */
}
stdvec2 GammaAnalyALG::GetCalData(CalibrationType cal, QString calName)
{
stdvec2 calData;
mat matData;
switch (cal)
{
case Cal_Energy: matData = Cal_Ener_Data[calName]; break;
case Cal_Resolution: matData = Cal_Reso_Data[calName]; break;
case Cal_Efficiency: matData = Cal_Effi_Data; break;
case Cal_Tot_efficiency: matData = Cal_Tot_Effi_Data; break;
}
if(matData.n_cols == 3)
{
calData.push_back( conv_to<stdvec>::from(matData.col(0u)) );
calData.push_back( conv_to<stdvec>::from(matData.col(1u)) );
calData.push_back( conv_to<stdvec>::from(matData.col(2u)) );
}
return calData;
}
void GammaAnalyALG::SetCalPara(CalibrationType cal, ParameterInfo para, QString calName)
{
field<vec> fPara(2);
fPara(0) = para.p;
if(para.perr.empty()) { fPara(1) = gNaN * ones<vec>(fPara(0).size()-1); }
else { fPara(1) = para.perr; }
switch (cal) {
case Cal_Energy: Acal_energy_para[calName] = fPara; break;
case Cal_Resolution: Acal_resolution_para[calName] = fPara; break;
case Cal_Efficiency: Acal_efficiency_para = fPara; break;
case Cal_Tot_efficiency: Acal_tot_efficiency_para = fPara; break;
case Cal_Tail: Acal_tail_para = vec(para.p); break;
case Cal_Tail_alpha: Acal_tail_alpha_para = vec(para.p); break;
case Cal_Tail_right: Acal_tail_right_para = vec(para.p); break;
case Cal_Tail_right_alpha: Acal_tail_right_alpha_para = vec(para.p); break;
case Cal_Step_ratio: Acal_step_ratio_para = vec(para.p);
default: break;
}
}
void GammaAnalyALG::SetCalPara(QStringList paraNames, ParameterInfo EnerPara, ParameterInfo ResoPara, ParameterInfo EffiPara, ParameterInfo TotEPara)
{
if(paraNames.size() < 2)
{
qDebug() << "The parameter paraNames must has at least 2 element in GammaAnalyALG::SetCalPara.";
return;
}
m_curEner = paraNames[0];
m_curReso = paraNames[1];
field<vec> f(2);
f(0) = EnerPara.p;
f(1) = EnerPara.perr;
Acal_energy_para[m_curEner] = f;
f(0) = ResoPara.p;
f(1) = ResoPara.perr;
Acal_resolution_para[m_curReso] = f;
f(0) = EffiPara.p;
f(1) = EffiPara.perr;
Acal_efficiency_para = f;
f(0) = TotEPara.p;
f(1) = TotEPara.perr;
Acal_tot_efficiency_para = f;
}
ParameterInfo GammaAnalyALG::GetCalPara(CalibrationType cal, QString calName)
{
ParameterInfo param;
field<vec> matPara;
switch(cal)
{
case Cal_Energy: matPara = Acal_energy_para[calName]; break;
case Cal_Resolution: matPara = Acal_resolution_para[calName]; break;
case Cal_Efficiency: matPara = Acal_efficiency_para; break;
case Cal_Tot_efficiency: matPara = Acal_tot_efficiency_para; break;
case Cal_Step_ratio: matPara << Acal_step_ratio_para; break;
case Cal_Tail: matPara << Acal_tail_para; break;
case Cal_Tail_alpha: matPara << Acal_tail_alpha_para; break;
case Cal_Tail_right: matPara << Acal_tail_right_para; break;
case Cal_Tail_right_alpha: matPara << Acal_tail_right_alpha_para; break;
}
if(matPara.size() > 0)
{
param.p = conv_to<stdvec>::from(matPara(0));
if(matPara.size() == 2) param.perr = conv_to<stdvec>::from(matPara(1));
}
return param;
}
field<vec> GammaAnalyALG::FitCalPara(CalibrationType cal, mat calData)
{
field<vec> resField(2);
if(calData.n_cols < 2) {
qDebug() << "Calibration data need at least two columns!";
return resField;
}
else if(calData.n_cols < 3)
{
calData.insert_cols(1, 2);
calData.col(2).fill(gNaN);
}
vec p, perr;
if(!Independ::calFitPara(p, perr, cal, calData.col(0u), calData.col(1u), calData.col(2u)))
{
bool b = Independ::lMakeInterp(p, perr, calData.col(0u), calData.col(1u), calData.col(2u));
if(!b || p.size()<2 || !Independ::calParaCheck(p))
{
p.clear();
perr.clear();
switch (cal)
{
case Cal_Energy:
{
p << 2 << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0;
if (m_nChans < 3000) p(2) = 1;
else if (m_nChans < 5000) p(2) = 0.5;
else if (m_nChans < 10000) p(2) = 0.33;
else p(2) = 0.2;
break;
}
case Cal_Resolution: p << 4 << 1 << 0.003 << 0; break;
case Cal_Efficiency:
case Cal_Tot_efficiency: p << 5 << 0.3 << 80 << 100 << 2.95 << 0.8; break;
default: qDebug() << "invalid calibration string in calFitPara";
}
}
}
resField(0) = p;
resField(1) = perr;
return resField;
}
/* =====================================================寻峰模块======================================================== */
uvec GammaAnalyALG::PeakSearch(double PSSs, PeakSearchType PStype, rowvec ELim)
{
vec C, NA, CNL, CNR, PSSfound;
return PeakSearch(C, NA, CNL, CNR, PSSfound, PSSs, PStype, ELim);
}
uvec GammaAnalyALG::PeakSearch(vec &C, vec &NA, vec &CNL, vec &CNR,
vec &PSSfound, double PSSs, PeakSearchType PStype, rowvec ELim)
{
uvec idx;
C.clear();
NA.clear();
CNL.clear();
CNR.clear();
PSSfound.clear();
// lower cutoff energy
double ECutLow = specSetup.ECutAnalysis_Low;
double ECutHigh = qInf();//specSetup.ECutAnalysis_High;
double deltaE = specSetup.EnergyTolerance;
// least net area estimate, lower NA estimates will be increased
double NAMin = 1; // min net area
//double fROI = 2;
double NWidth = 3;
// get default variables for parameters
if( qIsNaN(PSSs) ) PSSs = specSetup.PSS_low;
//if(PSSb < 0.000001) PSSb = PSSs;
// set parameter strings and flags for full or residual peak search
bool fullSearch;
char Src;
colvec COld; // 原来的Centroid
QString type, approxType;
if(PStype == search_residual)
{
fullSearch = false;
type = "Stripped";
approxType = "BaseLine";
Src = 'R';
}
else {
fullSearch = true;
// full peak search: get old centroids, don't insert old peaks again
if(pat.Peaks.size() > 0) COld = pat.Peaks.col(0); // COld = dmps("Centroid");
type = "Spectrum";
approxType = "Approximation";
Src = 'M';
}
//get spectrum counts, analysis range and baseline control points
//[spec, analysisRg, cx, cy, approx] = dmspec(sn, type, 'AnalysisRange', 'XControl', 'YControl', approxType);
QStringList Opts;
Opts << type << "AnalysisRange" << "XControl" << approxType;
field<rowvec> t_f = dmspec(Opts);
rowvec& spec = t_f(0);
rowvec analysisRg = t_f(1);
rowvec& cx = t_f(2);
//cx.print("cx");
//rowvec cy = t_f(3);
rowvec& approx = t_f(3);
//("peakSearch_spec.txt", spec);
//("peakSearch_approx.txt", approx);
// flag if we need to initialize
bool doBaseInit;
rowvec missedCounts;
if(cx.is_empty()) doBaseInit = true;
else {
doBaseInit = false;
missedCounts = spec - approx;
NWidth = 0;
}
uword NChan = spec.size();
rowvec Chan = RangeVec2(1, NChan);
//get energy, resolution and minimum distance in channel units for every channel
rowvec E = calValues(Cal_Energy, Chan);
//("peakSearch_E.txt", E);
rowvec dE = calDerivative(Cal_Energy, Chan);
//("peakSearch_dE.txt", dE);
rowvec res = calValues(Cal_Resolution, E);
//("peakSearch_res.txt", res);
if(res.size() != dE.size())
{
qDebug() << "The size of Energy Calibration must be Same as The size of Resolution Calibration!";
return idx;
}
rowvec res_ch = res / dE;
rowvec deltaC = deltaE / dE;
////("peakSearch_res_ch.txt", res_ch);
////("peakSearch_deltaC.txt", deltaC);
// transform ELim to allowed range (channels)
rowvec searchRg;
if(!ELim.is_empty())
{
searchRg = lGetAnalysisRange(ELim(0), ELim(1), spec, E, res_ch);
// make ELim the analysis range if not defined
if(analysisRg.is_empty()) analysisRg = searchRg;
}
// no analysis range defined; get it
if(analysisRg.is_empty())
{
analysisRg = lGetAnalysisRange(ECutLow, ECutHigh, spec, E, res_ch);
}
if(searchRg.is_empty()) searchRg = analysisRg;
else {
if(searchRg(0) < analysisRg(0)) searchRg(0) = analysisRg(0);
if(searchRg(1) > analysisRg(1)) searchRg(1) = analysisRg(1);
}
int ChanLow = searchRg(0);
int ChanHigh = searchRg(1);
// qDebug() << "ChanLow: " << ChanLow;
// qDebug() << "ChanHigh: " << ChanHigh;
// channels where peak search shall be performed
rowvec SearchChan = RangeVec2(ChanLow, ChanHigh);
//SearchChan.print("SearchChan");
// get peak search sensitivity in every search channel (sign so pss > 0 at peaks)
mat pss = - pksens(spec, res_ch.cols(ChanLow-1, ChanHigh-1), SearchChan, SearchChan);
pss.reshape(1, pss.size());
//("peakSearch_pss.txt", pss);
// look for regions of interest (ROIs)
int NPeaks = 0;
int i = 1; // i = 1;
int NSearch = SearchChan.size();
int ChanOffset = ChanLow - 1;
int l, L, r, R;
//pss.print("pss");
while(i < NSearch)
{
if(pss(i-1) > PSSs)
{
// we found a left ROI channel
l = i;
L = l + ChanOffset;
// find the right ROI channel
while(i < NSearch && pss(i-1) > PSSs)
{
i = i + 1;
}
r = i-1;
R = r + ChanOffset;
rowvec ROIidx = RangeVec2(l, r);
rowvec ROI = ROIidx + ChanOffset;
// Centroid ~ Center of Gravity
// COG = sum( ROI.*pss(ROIidx) ) / sum(pss(ROIidx));
rowvec temp_pss = pss.cols(l-1, r-1);
double COG = sum(ROI % temp_pss) / sum(temp_pss);
double maxPSS = temp_pss.max();
// accept, if residual peak search, or not too near existing peak
//... Centroid deviation max delta at channel floor(Centroid)
if( !fullSearch || all( abs(COld - COG) > deltaC( qFloor(COG) ) ) )
{
// add peak information to list
CNL.resize(NPeaks+1);
CNR.resize(NPeaks+1);
C.resize(NPeaks+1);
PSSfound.resize(NPeaks+1);
CNL(NPeaks) = L;
CNR(NPeaks) = R;
C(NPeaks) = COG;
PSSfound(NPeaks) = maxPSS;
NPeaks = NPeaks + 1;
}
}
i = i + 1;
}
/* ---------------------------------打印输出-------------------------------- */
//("peakSearch_CNL.txt", CNL);
//("peakSearch_CNR.txt", CNR);
//("peakSearch_C.txt", C);
//("peakSearch_PSSfound.txt", PSSfound);
NA = zeros<vec>(C.size());
// Get Net Area:
if(doBaseInit)
{
// get peak fwhm in channel units
vec nE = calValues(Cal_Energy, C);
vec ndE = calDerivative(Cal_Energy, C);
vec nFE = calValues(Cal_Resolution, nE);
vec FC = nFE / ndE; // mat FC = nFE ./ ndE;
/* ---------------------------------打印输出-------------------------------- */
//("peakSearch_FC.txt", FC);
// strip peaks from spectrum
rowvec baseCounts, roiLo, roiHi;
specCutROIs(baseCounts, roiLo, roiHi, spec, C, FC);
uword j = 0;
for(uword i = 0, j = 0; i<roiLo.size(); ++i)
{
while(j < C.size() && C(j) < roiHi(i))
{
// WORK: sumIdx may be wider, but should not exceed halfpoint to next peak
double yplus = sum( spec.cols(CNL(j)-1, CNR(j)-1) ); // double yplus = sum(spec(sumIdx));
double ybase = sum( baseCounts.cols(CNL(j)-1, CNR(j)-1) ); // double ybase = sum(baseCounts(sumIdx));
NA(j) = max(NAMin, yplus - ybase); // NA(j) = max(NAMin, yplus - ybase);
++j; //j = j+1;
}
}
} else {
for(uword i=0; i<C.size(); ++i)
{
NA(i) = max(NAMin, sum(missedCounts.cols(CNL(i)-1, CNR(i)-1))); // NA(i) = max(NAMin, temp_sum);
}
}
/* ---------------------------------打印输出-------------------------------- */
//("peakSearch_NA.txt", NA);
if(!C.is_empty())
{
//idx = patAddPeaks(sn, 'Source', Src, 'Centroid', C, 'NetArea', NA, 'Sensitivity', PSSfound, 'CentroidMethod', 'M', 'NetAreaMethod', 's');
uvec PropIdx;
PropIdx << i_Centroid << i_NetArea << i_Sensitivity;
mat Datas(C.size(), 3);
Datas.col(0) = C;
Datas.col(1) = NA;
Datas.col(2) = PSSfound;
uvec oldIdx;
patAddPeaks(PropIdx, Datas, pat.Peaks, idx, oldIdx);
}
if(doBaseInit)
{
rowvec tmp; tmp << 0 << 0;
tmpbase(analysisRg, tmp, tmp, analysisRg);
baseInit();
}
return idx;
}
mat GammaAnalyALG::calValues(CalibrationType cal, mat Chan)
{
mat y;
colvec p;
if(calGetPATPara(p, cal)) // get calibration parameter for named PAT
{
y = Independ::calFcnEval(Chan, p); // evaluate calibration function
}
return y;
}
mat GammaAnalyALG::calDerivative(CalibrationType cal, mat Chan)
{
mat dy = zeros(Chan.n_rows, Chan.n_cols);
colvec p;
if(calGetPATPara(p, cal)) // get calibration parameter for named PAT
{
dy = Independ::calDerivEval(Chan, p); // evaluate calibration function
}
return dy;
}
bool GammaAnalyALG::calGetPATPara(colvec& p, CalibrationType cal)
{
bool bRet;
if(cal == Cal_Default) // return default parameter, if name is invalid
{
bRet = calDefaultPara(p, cal);
}
else { bRet = getCalibrationPara(p, cal); }
return bRet;
}
bool GammaAnalyALG::getCalParaGlob(colvec &out, CalibrationType cal)
{
bool bRet = true;
// get global variable depending on calibration type
QString curCal = CalPHD;
switch (cal)
{
case Cal_Energy:
if(Acal_energy_para[m_curEner].size() > 0) out = Acal_energy_para[m_curEner](0);
break;
case Cal_Resolution:
if(Acal_resolution_para[m_curReso].size() > 0) out = Acal_resolution_para[m_curReso](0);
break;
case Cal_Efficiency:
if(Acal_efficiency_para.size() > 0) out = Acal_efficiency_para(0);
break;
case Cal_Tot_efficiency:
if(Acal_tot_efficiency_para.size() > 0) out = Acal_tot_efficiency_para(0);
break;
case Cal_Tail: out = Acal_tail_para; break;
case Cal_Tail_alpha: out = Acal_tail_alpha_para; break;
case Cal_Tail_right: out = Acal_tail_right_para; break;
case Cal_Tail_right_alpha: out = Acal_tail_right_alpha_para; break;
case Cal_Step_ratio: out = Acal_step_ratio_para; break;
default:
qDebug() << "invalid calibration string " << cal;
bRet = false;
}
// get list for given spectrum number
if (out.n_elem < 2)
{
bRet = false;
}
return bRet;
}
bool GammaAnalyALG::calDefaultPara(colvec& p, CalibrationType cal)
{
bool bRet = true;
// assign parameter p depending on calibration string
switch(cal)
{
case Cal_Energy:
p << 2 << 0 << 0.3 << 0; // linear
break;
case Cal_Resolution:
p << 4 << 1 << 0.003 << 0; // square root of polynomial
break;
case Cal_Efficiency:
case Cal_Tot_efficiency:
p << 94 << 0.1 << 30 << 3 << 50 << 70 << 0.8; // HAE Efficiency (1-2)
break;
case Cal_Tail:
p << 99 << 0 << 0 << 0; // cutoff linear, 0
break;
case Cal_Tail_alpha:
p << 98 << 1; // constant 1
break;
case Cal_Tail_right:
p << 98 << 0.35; // constant 0.35
break;
case Cal_Tail_right_alpha:
p << 98 << 1; // constant 1
break;
case Cal_Step_ratio:
p << 98 << 0; // constant 0
break;
case Cal_HalflifeAnalysis:
p << 97 << 0 << 1 << log(2); // exponential sum, no constant part, halflife 1 unit
break;
default:
bRet = false;
qDebug() << "invalid calibration string";
}
/* error estimates not reasonable for default calibration
perr = NaN * ones(1, length(p) - 1); */
return bRet;
}
bool GammaAnalyALG::getCalibrationPara(colvec& p, CalibrationType cal)
{
bool bRet = true;
// get all calibration parameters plus names for this spectrum
// [s, msg, calentries] = getCalParaGlob(sn, cal);
bRet = getCalParaGlob(p, cal);
if(!bRet) bRet = calDefaultPara(p, cal);
return bRet;
}
rowvec GammaAnalyALG::lGetAnalysisRange(double ECutLow, double ECutHigh, mat spec, mat E, mat res_ch)
{
// upper cutoff channels
int CutHighChans = 3;
// additional cutoff in resolution units
int CutC = 2;
int NChan = spec.n_elem;
// lower cutoff: at least ECutLow, but skip zero's
int ChanLow = 1;
while(ChanLow < NChan && spec(ChanLow-1) == 0)
{
ChanLow = ChanLow + 1;
}
// get additional distance from first non-zero channel
int t1 = ceil( CutC * res_ch(ChanLow-1) );
if(t1 > 0) ChanLow = ChanLow + t1;
// skip more channels if ChanLow is below ECutLow
if(ChanLow < NChan)
{
while( E(ChanLow-1) < ECutLow)
{
ChanLow = ChanLow + 1;
if(E.size()<=ChanLow)
{
throw std::string(" Mat::operator(): index out of bounds");
}
}
}
// ignore certain amount of channels in the end (may be checksum)
int ChanHigh = NChan - CutHighChans;
// skip zeros
while( ChanHigh > ChanLow && spec(ChanHigh-1) == 0)
{
ChanHigh = ChanHigh - 1;
}
// additional distance for peak search filter
int t2 = floor( CutC * res_ch(ChanHigh-1) );
if(t2 > 0) ChanHigh = ChanHigh - t2;
while(ChanHigh > ChanLow && E(ChanHigh-1) > ECutHigh)
{
ChanHigh = ChanHigh - 1;
}
if(ChanHigh <= ChanLow)
{
qDebug() << "Left Channel bigger than right channel, no peak search!";
ChanHigh = min(NChan, ChanLow + 1);
ChanLow = max(1, ChanHigh - 1);
}
rowvec rg;
rg << ChanLow << ChanHigh;
return rg;
}
colvec GammaAnalyALG::pksens(rowvec spec, mat res_c, mat LROI, mat RROI)
{
// Cut off filter at 1percent of total area
// int FilterCutoff = 1; // FilterCutoff=1;
// check input parameters
uword nROI = res_c.size(); // nROI = length(res_c);
if(LROI.size() != nROI || RROI.size() != nROI)
{
qDebug() << "2nd to 4th argument must be double vectors of same length";
}
res_c.reshape(nROI, 1);
LROI.reshape(nROI, 1);
RROI.reshape(nROI, 1);
uword NSpec = spec.size(); // NSpec=length(spec);
// myNaN=NaN;
colvec pss(nROI); // pss=myNaN(ones(size(res_c)));
pss.fill(gNaN);
// go through all ROIs
for(uword i=0; i<nROI; ++i) // for i=1:nROI
{
// get filter for that ROI
rowvec filt;
uword FilterMean = sdfilt(filt, res_c(i)); // [filt,FilterMean]=sdfilt(res_c(i));
// find the least sensitivity value in the ROI
double sens=qInf(); // sens=Inf;
for(uword ch=LROI(i); ch<=RROI(i); ++ch)
{
// find the channels that meet the filter, if centered in channel ch
urowvec idx = RangeVec(1, filt.size()); // idx=1:length(filt);
uvec iidx = find((ch+idx-FilterMean>0) && (ch+idx-FilterMean<=NSpec));
idx = IdxMat(idx, iidx); // idx=idx(iidx);
// apply the filter
urowvec idx_temp = idx - 1u;
urowvec idx_temp_2 = ch+idx_temp-FilterMean;
rowvec cc = IdxMat(filt, idx_temp) % IdxMat(spec,idx_temp_2); // cc=filt(idx).*spec(ch+idx-FilterMean);
double dd = sum(cc);
double sd = sqrt(sum(cc % IdxMat(filt,idx_temp)));
if(sd>0) // idx not empty
{
double ss = dd / sd;
// lower sensitivity has been found
if(ss < sens)
{
sens = ss;
}
}
}
// record pss for this ROI
pss(i) = sens; // pss(i)=sens;
}
return pss;
}
int GammaAnalyALG::sdfilt(rowvec &f, double RC, double CO)
{
// if resolution value bad, return zero filter
// if isempty(RC) | ~isfinite(RC) | RC < eps
uword k = 1; // index of filter mid point
if(qIsInf(RC) || RC < eps(1.0))
{
f << 0; // f = 0;
// k = 1;
return k;
}
double minPrec = 0.001;
// catch bad CO values
if(qIsInf(CO) || CO <= 0) // if isempty(CO) | ~isfinite(CO) | CO <= 0
{
CO = 1;
}
else if(CO < minPrec)
{
// limit precision (and thereby filter width)
CO = minPrec;
}
// transform resolution to sigma
double f2s = 4 * log(2); // f2s=4*log(2);
double c2 = pow(RC, 2) / f2s; // c2=RC^2/f2s;
// build one side of filter: tf
rowvec tf; tf << -100; // tf=-100;
// k=1;
double k2 = 1; // k2=1;
double v = 100*(k2-c2)*exp(-k2/(2*c2))/c2; // v=100*(k2-c2)*exp(-k2/(2*c2))/c2;
// build untill filter decreases and value is less than cutoff
while(v > min(tf(k-1u), CO)) // while v>min(tf(k),CO)
{
tf.resize(k+1u);
tf(k) = v; // tf(k+1)=v;
++k; // k=k+1;
k2 = pow(k, 2); // k2=k^2;
v = 100*(k2-c2)*exp(-k2/(2*c2))/c2;
}
// make it a 0-sum filter by applying deficiency to entry beside mid-point
double sum_tf;
if(k < 3) sum_tf = 0;
else sum_tf = 50-sum(tf.cols(2, k-1));
if(k < 2) tf.resize(2);
tf(1) = sum_tf; // tf(2)=50-sum(tf(3:k));
// create the symmetric filter arround mid-point
f = zeros<rowvec>(2*k-1); // f=zeros(1,2*k-1);
for(uword i=0; i<k; ++i)
{
f(i+k-1u) = tf(i); // f(k:2*k-1)=tf;
f(i) = tf(k-1u-i); // f(1:k-1)=tf(k:-1:2);
}
return k;
}
void GammaAnalyALG::specCutROIs(rowvec &base, rowvec &roiLo, rowvec &roiHi, rowvec spec, vec C, vec FC)
{
// number of channels for interpolation
int NWidth = 3; // NWidth = 3;
// width of roi in fwhm units
double fROI = 2.5; // fROI = 2.5;
uword len = spec.size();
base = spec; // base = spec;
// define ROIs
vec pkLo = max(1, floor(C - fROI*FC) - NWidth); // pkLo = max(1, floor(C - fROI*FC) - NWidth);
vec pkHi = min(len, ceil(C + fROI*FC) + NWidth); // pkHi = min(length(spec), ceil(C + fROI*FC) + NWidth);
rowvec roiBool = zeros<rowvec>(len); // roiBool = zeros(size(spec));
for(uword i=0; i<C.size(); ++i) // for i = 1 : length(C)
{
// roiBool(pkLo(i) : pkHi(i)) = 1;
for(int j=pkLo(i)-1; j<pkHi(i); ++j)
{
if(j < 0) continue;
if(j >= len) roiBool.resize(j+1);
roiBool(j) = 1;
}
}
// roiBool([1,end+1]) = 0;
roiBool.resize(roiBool.size()+1);
roiBool(0) = 0;
roiBool(roiBool.size()-1) = 0;
rowvec diff_roiBool = diff(roiBool);
uvec temp_1 = find(diff_roiBool == 1);
uvec temp_2 = find(diff_roiBool == -1);
roiLo.set_size(temp_1.size());
roiHi.set_size(temp_2.size());
for(int i=0; i<temp_1.size(); ++i)
{
roiLo(i) = temp_1(i);
}
for(int i=0; i<temp_2.size(); ++i)
{
roiHi(i) = temp_2(i);
}
for(uword i=0; i<roiLo.size() && i<roiHi.size(); ++i) // for i = 1 : length(roiLo)
{
double yLo = sum(spec.cols(roiLo(i), roiLo(i)+NWidth-1))/NWidth; // yLo = sum(spec(roiLo(i) : roiLo(i)+NWidth-1))/NWidth;
double yHi = sum(spec.cols(roiHi(i)-NWidth+1, roiHi(i)))/NWidth; // yHi = sum(spec(roiHi(i)-NWidth+1 : roiHi(i)))/NWidth;
rowvec localSum = cumsum(spec.cols(roiLo(i), roiHi(i))); // localSum = cumsum(spec(roiLo(i) : roiHi(i)));
// base(roiLo(i) : roiHi(i)) = yLo + (yHi - yLo) * localSum / localSum(end);
rowvec temp = yLo + (yHi - yLo) * localSum / localSum(localSum.size()-1);
for(int j=roiLo(i); j<=roiHi(i); ++j)
{
base(j) = temp(j-roiLo(i));
}
}
}
bool GammaAnalyALG::patAddPeaks(uvec Idx, mat Datas, mat &peaks, uvec &nidx, uvec &oidx)
{
// 如果索引的个数与数据的个数不对应则返回
if(Datas.n_cols != Idx.size())
{
qDebug() << "The size of Datas and Idx must be same!";
return false;
}
/* 创建一个新的存储峰信息的矩阵,然后把传入的数据 Datas 插入到矩阵对应列(由 Idx 指定)中
* nold 峰列表中峰的个数, nnew 是新插入峰的个数 */
int nold = peaks.n_rows;
int nnew = Datas.n_rows;
mat newPeaks(nnew, MAX_INDEX_PEAK);
newPeaks.fill(gNaN);
for(uword i=0; i<Idx.size(); ++i)
{
if(Idx(i) > MAX_INDEX_PEAK) continue;
newPeaks.col(Idx(i)) = Datas.col(i);
}
if(!newPeaks.col(i_Centroid).is_finite())
{
qDebug() << "Centroid required";
return false;
}
if(!newPeaks.col(i_FWHM_Ch).is_finite())
{
mat mt = newPeaks.col(i_Centroid);
// PrintMat22("newPeaks.col(i_Centroid).txt", mt);
colvec E = calValues(Cal_Energy, newPeaks.col(i_Centroid));
colvec dE = calDerivative(Cal_Energy, newPeaks.col(i_Centroid));
colvec F = calValues(Cal_Resolution, E);
newPeaks.col(i_FWHM_Ch) = F / dE; // FC = F ./ dE;
// PrintColvec("E.txt", E);
// PrintColvec("dE.txt", dE);
// PrintColvec("F.txt", F);
vec pFC = newPeaks.col(i_FWHM_Ch);
// PrintMat22("newPeaks.col(i_FWHM_Ch).txt", pFC);
}
if(!newPeaks.col(i_NetArea).is_finite())
{
qDebug() << "NetArea required";
return false;
}
if(!newPeaks.col(i_Sensitivity).is_finite())
{
qDebug() << "Sensitivity unknown, setting to NaN";
}
peaks = join_vert(peaks, newPeaks);
uvec index = sortrows(peaks, peaks); // [patjoint, index] = sortrows(patjoint);
//[sindex, order] = sort(index);
mat sindex;
uvec order = sortrows(sindex, UMatToMat(index));
nidx = order.rows(nold, order.size()-1);// nidx = order(nold + 1 : end);
if(nold > 0) oidx = order.rows(0, nold-1); // oidx = order(1 : nold);
colvec sqi, uqi;
patBaseVar(sqi, uqi, nidx);
return true;
}
/* =====================================================基线初始化======================================================= */
void GammaAnalyALG::patBaseVar(colvec &sqi, colvec &uqi, uvec idx, int updr)
{
// patin and patout always assigned, even when nargin == 4
/*if nargin < 3
if nargin < 2
idx = 'all';
end
patin = 'Current';
patout = 'Temporary';
else
patin = pat;
patout = pat;
end*/
// setup parameters: back-counts width in fwhm-units, sensitivity factor
//[k, k_alpha, k_beta] = getSpecSetup(sn, 'k_back', 'k_alpha', 'k_beta');
double k = specSetup.k_back;
double k_alpha = specSetup.k_alpha;
double k_beta = specSetup.k_beta;
// get baseline and peak parameters
field<rowvec> t_f1;
if(updr)
{
//[BL, spec, resid] = dmspec(sn, 'BaseLine', 'Spectrum', 'Residual');
QStringList specItems;
specItems << "BaseLine" << "Residual";
t_f1 = dmspec(specItems);
/* -------------------------------打印输出 ---------------------------------- */
//("patBaseVar_BaseLine.txt", t_f1(0));
//("patBaseVar_Residual.txt", t_f1(1));
}
else {
t_f1 = dmspec(QStringList("BaseLine")); //BL = dmspec(sn, 'BaseLine');
/* -------------------------------打印输出 ---------------------------------- */
//("patBaseVar_BaseLine.txt", t_f1(0));
}
//[Ct, Fc, NA, Sc] = dmps(sn, patin, 'Centroid', 'FWHM_Ch', 'NetArea', 'SigmaCh');
QStringList dmpsItems;
dmpsItems << "Centroid" << "FWHM_Ch" << "NetArea" << "SigmaCh";
field<vec> t_f2 = dmps(dmpsItems);
/* // modify all peaks, if idx is string
if isstr(idx)
idx = [1 : length(Ct)];
end*/
if(t_f2(0).is_empty()) return; // no peaks in PAT or called needlessly
if(idx.is_empty())
{
idx = vectorise( RangeVec(0, t_f2(0).size()-1) );
}
vec MBC, LC, LD, BC;
if(t_f1(0).is_empty()) //if isempty(BL)
{
// no baseline, baseline parameters are undefined
MBC << gNaN;
//S = gNaN;
//Err = gNaN;
LC << gNaN;
LD << gNaN;
sqi << gNaN;
uqi << gNaN;
}
else {
// restrict to indexed peaks
t_f2(0) = t_f2(0)(idx); //Ct = Ct(idx);
t_f2(1) = t_f2(1)(idx); //Fc = Fc(idx);
t_f2(2) = t_f2(2)(idx); //NA = NA(idx);
t_f2(3) = t_f2(3)(idx); //Sc = Sc(idx);
vec W = 2 * k * t_f2(1); //W = 2 * k * Fc;
/* -------------------------------打印输出 ---------------------------------- */
// total background counts from Ct +- k*Fc
rowvec tmp1 = max(0, t_f1(0));
BC = Independ::abkcnt( tmp1, t_f2(0), t_f2(1), k); //BC = abkcnt(max(BL, 0), Ct, Fc, k);
/* -------------------------------打印输出 ---------------------------------- */
// Mean Back Counts
MBC = BC / W;
/* -------------------------------打印输出 ---------------------------------- */
/* standard deviation of A is:
* D(A) = c * f * sqrt(w*b) with
* w - sigma of gaussian peak (channels)
* b - mean baseline/channel under peak
* f - algorithm-dependent factor if baseline known
* c - baseline estimation correction factor;
* formerly: D(A)=sqrt(m*2*kb*fc*b), thus
* f = sqrt(m*2*kb*sqrt(8*log(2))),
* yields: f=1.9562 for m=0.65, kb=1.25 */
double c = 1.1;
double f = 1.9562;
vec DA = c*f*sqrt(t_f2(3) % MBC); //DA = c*f*sqrt(Sc .* MBC);
/* -------------------------------打印输出 ---------------------------------- */
//("patBaseVar_DA.txt", DA);
// Currie's LC
LC = k_alpha * DA;
// Currie's LD
// old formula: LD = k_alpha^2 + 2* LC;
double kb2 = pow(k_beta, 2);
// assymmetric formula, ref TBD
LD = LC + (kb2/2) * ( 1 + sqrt(1 + (4/kb2)*( LC+pow(DA,2) ) ) );
/* -------------------------------打印输出 ---------------------------------- */
//("patBaseVar_LD.txt", LD);
if(updr)
{
// total counts
vec STot = Independ::abkcnt(baseInfo(s_Spectrum), t_f2(0), t_f2(1), k); //STot = abkcnt(spec, Ct, Fc, k);
// signed residual
vec signedResidual = Independ::abkcnt(t_f1(1), t_f2(0), t_f2(1), k);
// unsigned residual
vec unsignedResidual = Independ::abkcnt(abs( t_f1(1) ), t_f2(0), t_f2(1), k);
/* -------------------------------打印输出 ---------------------------------- */
//("patBaseVar_STot.txt", STot);
//("patBaseVar_signedResidual.txt", signedResidual);
//("patBaseVar_unsignedResidual.txt", unsignedResidual);
vec sqrtCT = sqrt(STot);
sqi = signedResidual / sqrtCT;
uqi = unsignedResidual / sqrtCT;
}
}
// write results to temporary PAT
if(updr) lSetPat(idx, MBC, BC, LC, LD, sqi, uqi);
else lSetPat(idx, MBC, BC, LC, LD);
}
void GammaAnalyALG::lSetPat(uvec idx, vec mbc, vec cba, vec lc, vec ld, vec sqi, vec uqi)
{
/*% get PAT table
[s, msg, patEntry] = getPATByName(sn, patin);
if s
pd = patEntry{2};
% write PAT entries, compare helpApat_list */
//MatAssign(pat.Peaks.col(i_MeanBackCount), idx, mbc); //pd(idx, 14) = mbc;
//MatAssign(pat.Peaks.col(i_LC), idx, lc); //pd(idx, 34) = lc;
//MatAssign(pat.Peaks.col(i_LD), idx, ld); //pd(idx, 35) = ld;
//MatAssign(pat.Peaks.col(i_BackgroundArea), idx, cba); //pd(idx, 41) = cba;
for(uword i=0; i<idx.size(); ++i)
{
pat.Peaks.col(i_MeanBackCount)(idx(i)) = mbc(i);
pat.Peaks.col(i_LC)(idx(i)) = lc(i);
pat.Peaks.col(i_LD)(idx(i)) = ld(i);
pat.Peaks.col(i_BackgroundArea)(idx(i)) = cba(i);
pat.Peaks.col(i_Significance)(idx(i)) = gNaN;
}
// either update residuals, or raise flag that they are not up to date
if(!sqi.is_empty() && !uqi.is_empty()) //if nargin == 10
{
//MatAssign(pat.Peaks.col(i_SignedResidual), idx, sqi); //pd(idx, 36) = sqi;
//MatAssign(pat.Peaks.col(i_UnsignedResidual), idx, uqi); //pd(idx, 37) = uqi;
//MatAssign(pat.Peaks.col(i_ResidualUpToDate), idx, 1); //pd(idx, 38) = 1;
for(uword i=0; i<idx.size(); ++i)
{
pat.Peaks.col(i_SignedResidual)(idx(i)) = sqi(i);
pat.Peaks.col(i_UnsignedResidual)(idx(i)) = uqi(i);
pat.Peaks.col(i_ResidualUpToDate)(idx(i)) = 1;
}
}
else {
//MatAssign(pat.Peaks.col(i_ResidualUpToDate), idx, 0); //pd(idx, 38) = 0;
for(uword i=0; i<idx.size(); ++i)
{
pat.Peaks.col(i_ResidualUpToDate)(idx(i)) = 0;
}
}
// significance overwritten
//MatAssign(pat.Peaks.col(i_Significance), idx, gNaN); //pd(idx, 16) = NaN;
/*patEntry{1} = patout;
patEntry{2} = pd;
[s, msg] = setPATEntry(sn, patEntry);
if ~s
logWrite('Failed writing Baseline PAT variables: %s', msg);
end
else
logWrite('Problem while updating Baseline PAT variables: %s', msg);
end*/
}
void GammaAnalyALG::baseInit(double PSSBreak, double PSSfwhm)
{
rowvec cx, cy, cdy;
//% least number of data points to be fitted for 1 CP
double NMin = 5;
//% absolute minimum distance of CPs
double NCritical = 10;
//% minimum number of data points between ROIs to use 3 CPs
double NTwo = 50; //% should come from data condition
rowvec data;
rowvec step;
rowvec rg;
rowvec roiMask;
uvec L;
uvec H;
stdvec vy;
vy = getBaseInfo(s_YControl);
lGetData(data, step, rg, roiMask, L, H, PSSBreak, PSSfwhm);
// PrintMat22("baseInit_data.txt", data);
// PrintMat22("baseInit_step.txt", step);
// PrintMat22("baseInit_rg.txt", rg);
// PrintMat22("baseInit_roiMask.txt", roiMask);
vy = getBaseInfo(s_YControl);
lInitBase(data, step, rg, roiMask, L, H, NMin, NCritical, NTwo, cx, cy, cdy);
// PrintMat22("baseInit_cx.txt", cx);
// PrintMat22("baseInit_cy.txt", cy);
// PrintMat22("baseInit_cdy.txt", cdy);
vy = getBaseInfo(s_YControl);
tmpbase(cx,cy,cdy);
vy = getBaseInfo(s_YControl);
return;
}
void GammaAnalyALG::tmpbase(mat cx, mat cy, mat cdy, mat ar)
{
//error(nargchk(4,5,nargin));
/*if(ar.is_empty()) //if nargin < 5
{
if(!cx.is_empty()) // isempty(cx)
{
ar = baseInfo(s_AnalysisRange);
}
}*/
if(!ar.is_empty())
{
if(ar.is_colvec()) baseInfo(s_AnalysisRange) = ar.t();
baseInfo(s_AnalysisRange) = ar;
}
// create baseline structure
/*global Ac_baseline_tmp
Ac_baseline_tmp=struct('AnalysisRange', {ar}, ...
'XControl', {cx(:)'}, ...
'YControl', {cy(:)'}, ...
'YSlope', {cdy(:)'}, ...
'SpectrumNumber', {sn});*/
if(cx.is_colvec()) baseInfo(s_XControl) = cx.t();
else baseInfo(s_XControl) = cx;
if(cy.is_colvec()) baseInfo(s_YControl) = cy.t();
else baseInfo(s_YControl) = cy;
if(cdy.is_colvec()) baseInfo(s_YSlope) = cdy.t();
else baseInfo(s_YSlope) = cdy;
// update baseline dependent variables
vec sqi, uqi;
patBaseVar(sqi, uqi);
}
void GammaAnalyALG::findROI(urowvec& l, urowvec& h, rowvec& roiMask,
vec fLo, uvec idx, vec fHi)
{
if(fHi.is_empty()) fHi = fLo;
QStringList Opts; Opts << "Centroid" << "FWHM_Ch";
field<vec> t_dmps = dmps(Opts);
vec& C = t_dmps(0);
vec& FWHM_CH = t_dmps(1);
rowvec rg = baseInfo(s_AnalysisRange);
if(idx.is_empty()) idx = RangeVec(0, C.size()-1).t();
vec temp_FWHM = FWHM_CH(idx);
vec temp_C = C(idx);
vec tempL = fLo % temp_FWHM;
vec tempH = fHi % temp_FWHM;
vec pkLo = max(rg(0), floor(temp_C - tempL));
vec pkHi = min(rg(1), ceil(temp_C + tempH));
//("findROI_pkLo.txt", pkLo);
//("findROI_pkHi.txt", pkHi);
roiMask = zeros<rowvec>(m_nChans);
for(size_t i=0; i!=pkLo.size(); ++i)
{
roiMask(span(pkLo(i)-1, pkHi(i)-1)) = ones<rowvec>(pkHi(i)-pkLo(i)+1);
}
roiMask[0] = 0;
roiMask[roiMask.size()-1] = 0;
rowvec db = diff(roiMask);
l = find(db == 1).t() + 1;
h = find(db == -1).t();
}
void GammaAnalyALG::lGetData(rowvec& data, rowvec& step, rowvec& rg, rowvec& roiMask, uvec& L, uvec& H,
double PSSBreak, double PSSfwhm)
{
double maxFWHMfact = 5;
// same as matlab [stripped, step, rg] = dmspec(sn, 'Stripped', 'Steps', 'AnalysisRange');
QStringList Opts;
Opts << "Stripped" << "Steps" << "AnalysisRange";
field<rowvec> t_f = dmspec(Opts);
rowvec& stripped = t_f(0);
step = t_f(1);
// PrintMat22("stripped.txt", stripped);
// PrintMat22("step.txt", step);
rg = t_f(2);
if(rg.is_empty())
{
rg = RangeVec2(1, m_nChans);
}
// same as matlab [pC, pPSS, pFC] = dmps(sn, 'Centroid', 'Sensitivity', 'FWHM_Ch');
Opts.clear();
Opts << "Centroid" << "Sensitivity";
field<vec> f_dmps = dmps(Opts);
vec& pC = f_dmps(0);
vec& pPSS = f_dmps(1);
//vec pFC = pat.Peaks.col(i_FWHM_Ch);
data = stripped - step;
// get peak ranges
colvec factFWHM = 1 + (maxFWHMfact -1)*pPSS/PSSfwhm;
factFWHM = min(maxFWHMfact,factFWHM); // if factFWHM's members bigger then maxFWHMfact set them equal maxFWHMfact
factFWHM = NaNto1(factFWHM); //factFWHM(isnan(factFWHM)) = 1;
//("lGetData_factFWHM.txt", factFWHM);
urowvec l;
urowvec h;
findROI(l, h, roiMask, factFWHM); // l,h,roiMask are out params
l = l + 1u;
h = h + 1u;
//("lGetData_roiMask.txt", roiMask);
// get borders of ranges containing big peaks
urowvec isbig = zeros<urowvec>(l.size()); // isbig = false(size(l));
uvec idx;
for(size_t i=0; i<l.size(); ++i) // for i = 1 : length(l)
{
idx = find(pC > l(i) && pC < h(i)); // idx = find(pC > l(i) & pC < h(i));
if(any(pPSS(idx) > PSSBreak)) isbig(i) = 1; //if any(pPSS(idx) > PSSBreak), isbig(i) = true; end
}
idx = find(isbig > 0);
L = l(idx); // L = l(isbig);
H = h(idx); // H = h(isbig);
// remove ROIs that end to near lower cutoff
idx = find(H <= (rg(0) + 2)); // idx = find(H <= (rg(1)+2));
if(!idx.is_empty()) // if any(idx), L(idx)=[]; H(idx)=[]; end
{
for(int i=idx.size()-1; i>=0; --i)
{
L.shed_col(idx(i));
H.shed_col(idx(i));
}
}
// shift ROI limit if first ROI goes beyond lower cutoff
if((!L.is_empty()) && (L[0] <= rg[0]))
{
L[0] = rg[0] + 1;
}
// add lower cutoff as a ROI
uvec tmp;
tmp << 1u;
L = join_vert(tmp, L); // L = [1; L(:)];
tmp << rg(0);
H = join_vert(tmp, H); // H = [rg(1); H(:)];
}
mat GammaAnalyALG::restpolyfit(vec& r, vec x, vec y, vec p0, vec w)
{
if( x.size() != y.size() )
{
qDebug() << "X and Y vectors must be the same size."; // messagebox
}
if(w.is_empty())
{
w = ones<vec>(x.size());
}
else {
y = y % w;
}
int n = p0.size()-1;
// Construct 'weighted' Vandermonde matrix.
mat V = zeros<mat>(w.size(),n+1);
V.col(n) = w;
for(int i=n-1; i>=0; --i)
{
V.col(i) = x % V.col(i+1);
}
// get free and fixed component
uvec freeIdx = find_nonfinite(p0);
uvec fixIdx = find_finite(p0);
mat p;
vec yfix = MatCols(V, fixIdx) * p0(fixIdx);//yfix = V(:, fixIdx)*p0(fixIdx);
vec yfree = y - yfix;
mat Q, R;
qr_econ(Q, R, MatCols(V, freeIdx));
p = solve(R, trans(Q) * yfree); // 待测试P应为行向量
colvec yhat = yfix + MatCols(V, freeIdx) * p;
r = y - yhat;
p0(freeIdx) = p;
p = p0.t();
return p;
}
void GammaAnalyALG::lNorm(rowvec yy, rowvec mask, double xc, double xh, double yh, double np,
double& nrm, double& lbda, double& yc, double& dyc, double& xl, double& yho, double& dyh)
{
// check argument, get polynomial mask
mat p;
if(np == 4)
{
p << gNaN << gNaN << gNaN << yh;
}
else if(np == 2)
{
p << gNaN << yh;
}
else qDebug() << "invalid np";
// get whole fitting interval
double k = xh - xc;
xl = max(1, xh-2*k); //xl = max(1, xh-2*k);
rowvec xx = RangeVec2(xl-1,xh-1); //xx = [xl : xh];
// remove flagged points
uvec idx = find( IdxMat(mask,xx) == true );
for(int i=idx.size()-1; i>=0; --i)
{
xx.shed_col(idx(i));
}
if(xx.size()<5)
{
nrm = 0;
lbda = 1;
yc = yy(xc);
dyc = 0;
yh = yy(xh);
dyh = 0;
return;
}
// limit to actual points
rowvec y = IdxMat(yy,xx);
rowvec x = xx - xh + 1;
// perform fit
colvec r;
p = restpolyfit(r, vectorise(x), vectorise(y), vectorise(p));
// polyder.m
rowvec dp, t_null;
Matlab::polyder(dp, t_null, 1, p);
nrm = accu(r % r);
lbda = max(1, mean(sqrt(abs(y))));
rowvec rg; rg << -k << 0;
rowvec yval = Matlab::polyval(p, rg);
rowvec dyval = Matlab::polyval(dp, rg);
yc = yval(0);
yho = yval(1);
dyc = dyval(0);
dyh = dyval(1);
}
void GammaAnalyALG::lSelectNext(double x, double y, rowvec data, double xb, rowvec imask, double np,
double& xn, double& yn, double& dyn, bool& hit, double& y1, double& dy1)
{
// selection algorithm parameters:
// lambda multiplier for norm threshold
double NMax = 4000;
// minimum interval length
double MinDist = 10;
// get next step
double hi = max(1, x - MinDist);
xn = max(1, hi - 1);
// binary search
// get norms for initial new x
double norm, lambda, x0, y0, dy0;
//[norm, lambda, yn, dyn, x0, y0, dy0] = lNorm(data, imask, xn, x, y, np);
lNorm(data, imask, xn, x, y, np, norm, lambda, yn, dyn, x0, y0, dy0);
while(hi != xn)
{
if(norm < NMax * lambda)
{
double tmp = xn;
xn = max(1, 3*xn - 2*hi);
hi = tmp;
}
else
{
xn = hi - floor((hi-xn)/2);
}
lNorm(data, imask, xn, x, y, np, norm, lambda, yn, dyn, x0, y1, dy1);
}
if(xn < xb)
{
xn = xb;
lNorm(data, imask, xn, x, y, np, norm, lambda, yn, dyn, x0, y1, dy1);
}
hit = (x0 <= xb);
}
void GammaAnalyALG::lFitRange(double x0, double y0, double x1, rowvec xv, rowvec data, rowvec mask, int npara,
rowvec& yv, rowvec& dyv)
{
rowvec xx = RangeVec2(min(x0, x1), max(x0,x1));
//xx(mask(xx)==true) = [];
rowvec t_xx = xx - 1;
uvec idx = find( IdxMat(mask,t_xx) == true );
for(int i=idx.size()-1; i>=0; --i)
{
if(idx(i) < xx.size())
{
xx.shed_col(idx(i));
t_xx.shed_col(idx(i));
}
}
if(xx.size()<5)
{
qDebug() << QString("Warning: only %1 points between %2 and %3")
.arg(xx.size()).arg(min(x0,x1)).arg(max(x0,x1));
rowvec t_xv = xv - 1;
yv = IdxMat(data, t_xv);
dyv = zeros( size(xv) );
return;
}
// fit linear or cubic polynomial
rowvec p0;
switch(npara)
{
case 2: p0 << gNaN << y0; break;
case 4: p0 << gNaN << gNaN << gNaN << y0; break;
default: qDebug() << "Bug: invalid npara value"; break;
}
// fit polynomial
rowvec yy = IdxMat(data, t_xx);
xx = xx - x0;
colvec r;
mat p = restpolyfit(r, vectorise(xx), vectorise(yy), vectorise(p0));
// evaluate
yv = Matlab::polyval(p, xv - x0);
mat a, b;
Matlab::polyder(a, b, 1, p);
dyv = Matlab::polyval(a, xv - x0);
}
void GammaAnalyALG::lInitBase(const rowvec& data, const rowvec& step, const rowvec& rg,
const rowvec& roiMask, const uvec& L, const uvec& H,
double N1, double N2, double N3,
rowvec& cx, rowvec& cy, rowvec& cdy)
{
size_t roiIdx = L.size();
double xb = H[roiIdx-1];
cx.set_size(0);
cy.set_size(0);
cdy.set_size(0);
rowvec yn, dyn, yo, dyo, tmp;
int nP = 0;
bool pingold = true, pingNew = false;
double x = rg[1], xn, y = gNaN, dy = gNaN;
double t_yn, t_dyn, t_yo, t_dyo;
while(x > rg[0])
{
// get next controlpoints
if(pingold) nP = 2;
else nP = 4;
//[xn, yn, dyn, pingNew, yo, dyo] = lSelectNext(x,y,data,xb,roiMask,nP);
lSelectNext(x, y, data, xb, roiMask, nP, xn, t_yn, t_dyn, pingNew, t_yo, t_dyo);
yn << t_yn; dyn << t_dyn; yo << t_yo; dyo = t_dyo;
if(pingold)
{
// starting left of a ROI
if(pingNew)
{
// delta = x - xb + 1;
double delta = accu( roiMask.cols(xb,x) == false );
if(delta < N2)
{
if(cx.is_empty())
{
// highest controlpoint, must be introduced
//[cy, cdy] = lFitRange(x,NaN,xb, x, data, roiMask, 2);
tmp << x;
lFitRange(x, gNaN, xb, tmp, data, roiMask, 2, cy, cdy);
cx << x;
}
else if(roiIdx == 1)
{
// lowest controlpoint, must be introduced
//[yn, dyn] = lFitRange(xb,NaN,x, xb, data, roiMask, 2);
tmp << xb;
lFitRange(xb, gNaN, x, tmp, data, roiMask, 2, yn, dyn);
cx = join_horiz(tmp, cx); //cx = [xb, cx];
cy = join_horiz(yn, cy); //cy = [yn, cy];
cdy = join_horiz(dyn, cdy); //cdy = [dyn, cdy];
}
else
{
if(delta < N1)
{
// no control point introduced
}
else
{
// introducing only one control point
double xc = round((x + xb)/2);
tmp << xc;
lFitRange(xb,gNaN,x, tmp, data, roiMask, 2, yn, dyn);
cx = join_horiz(tmp, cx); //cx = [xc, cx];
cy = join_horiz(yn, cy); //cy = [yn, cy];
cdy = join_horiz(dyn, cdy); //cdy = [dyn, cdy];
}
}
}
else if(delta < N3)
{
//[yn, dyn] = lFitRange(x, y, xb, [xb, x], data, roiMask, 2);
tmp << xb << x;
lFitRange(x, y, xb, tmp, data, roiMask, 2, yn, dyn);
cx = join_horiz(tmp, cx); //cx = [xb, x, cx];
cy = join_horiz(yn, cy); //cy = [yn, cy];
tmp << gNaN << gNaN;
cdy = join_horiz(tmp, cdy); //cdy = [NaN, NaN, cdy];
}
else
{
// introducing three control points
double xc = round((x + xb)/2);
//[yn, dyn] = lFitRange(xb,NaN,x,[xb,xc,x],data,roiMask,4);
tmp << xb << xc << x;
lFitRange(xb, gNaN, x, tmp, data, roiMask, 4, yn, dyn);
cx = join_horiz(tmp, cx); //cx = [xb,xc,x,cx];
cy = join_horiz(yn, cy); //cy = [yn, cy];
tmp << dyn(0) << gNaN << gNaN;
cdy = join_horiz(tmp, cdy); //cdy = [dyn(1), NaN, NaN, cdy];
}
if(roiIdx == 1)
{
break;
}
x = L(roiIdx -1);
y = gNaN;
dy = gNaN;
roiIdx = roiIdx -1;
xb = H(roiIdx -1);
}
else {
// normal algorithm left of ROI
tmp << x;
cx = join_horiz(tmp, cx); //cx = [x, cx];
cy = join_horiz(yo, cy); //cy = [yo, cy];
cdy = join_horiz(dyo, cdy); //cdy = [dyo, cdy];
// continue from here
x = xn;
y = t_yn;
dy = t_dyn;
}
}
else
{
// starting from normal CP
if(pingNew)
{
// hitting a ROI
if(roiIdx == 1)
{
// hitting lower analysis range
if(xn == xb)
{
// next control point at end
//[yn, dyn] = lFitRange(x, y, xb, xb, data,roiMask, 2);
tmp << xb;
lFitRange(x, y, xb, tmp, data, roiMask, 2, yn, dyn);
cx = join_horiz(tmp, cx); //cx = [xb, cx];
cy = join_horiz(yn, cy); //cy = [yn, cy];
cdy = join_horiz(dyn, cdy); //cdy = [dyn, cdy];
}
else
{
double xm = round((x+xb)/2);
tmp << xb << xm;
lFitRange(x, y, xb, tmp, data, roiMask, 4, yn, dyn);
cx = join_horiz(tmp, cx); //cx = [xb, xm, cx];
cy = join_horiz(yn, cy); //cy = [yn, cy];
tmp << dyn(0) << gNaN;
cdy = join_horiz(tmp, cdy); //cdy = [dyn(1), NaN, cdy];
}
break;
}
else if(xn == xb)
{
// hitting ROI completely
tmp << xn; cx = join_horiz(tmp, cx); //cx = [xn, cx];
cy = join_horiz(yn, cy); //cy = [yn, cy];
cdy = join_horiz(dyn, cdy); //cdy = [dyn, cdy];
}
else {
double xc = ceil((x+xb)/2);
tmp << xb << xc;
lFitRange(x, y, xb, tmp, data, roiMask, 4, yn, dyn);
//add righthand CP, with fixed dy;
cx = join_horiz(tmp, cx); //cx = [xb, xc, cx];
cy = join_horiz(yn, cy); //cy = [yn, cy];
tmp << dyn(0) << gNaN;
cdy = join_horiz(tmp, cdy); //cdy = [dyn(1), NaN, cdy];
}
// continue on lefthand, starting with free slope
x = L(roiIdx-1);
y = gNaN;
dy = gNaN;
roiIdx = roiIdx - 1;
xb = H(roiIdx-1);
}
else
{
// normal procedure
tmp << xn; cx = join_horiz(tmp, cx); //cx = [xn, cx];
cy = join_horiz(yn, cy); //cy = [yn, cy];
tmp << gNaN; cdy = join_horiz(tmp, cdy); //cdy = [NaN, cdy];
// continue from here
x = xn;
y = t_yn;
dy = t_dyn;
}
}
pingold = pingNew;
}
if( qIsNaN( cdy(0) ) )
{
cdy(0) = 0;
}
if( qIsNaN( cdy( cdy.size()-1 ) ) )
{
cdy( cdy.size()-1 ) = 0;
}
if(cx.size() == 1)
{
cx = rg;
double tmp = cy(0);
cy << tmp << tmp;
cdy << 0 << 0;
qDebug() << "Whole baseline was fitted at once";
}
rowvec t_cx = cx - 1;
cy = cy + IdxMat(step, t_cx);
}
/* =====================================================基线优化模块======================================================== */
void GammaAnalyALG::baseImprove(double pss_r)
{
// accept old temporary results to make sure they are not lost
//patCommit(sn);
vec sqi, uqi;
patBaseVar(sqi, uqi); //tbaseac(sn);
// search peaks in residual, written to temporary PAT
uvec idx = PeakSearch(pss_r, search_residual);
//idx.print("idx = ");
// fit net areas of new peaks
if(!idx.is_empty())
{
uvec ep;
fitPeakFull(idx, ep, ep, ep);
}
// fit the baseline to the residual
FitBase();
// undo insertion of residual peaks
//patRollback;
int size_idx = idx.size();
for(int i=size_idx-1; i>=0; --i)
{
pat.Peaks.shed_row(idx(i));
}
}
void GammaAnalyALG::fitPeakFull(uvec Af, uvec Cf, uvec Ff, uvec SRf, uvec CPf, bool bVerbose)
{
QString addArg; //addArg = {};
if(bVerbose) addArg = "Peaks above spline baseline";
// get index to all peaks where any parameter is free
// get common min and max
uvec allidx = join_vert( join_vert(Af, Cf), join_vert(Ff, SRf) );
if(allidx.is_empty()) return;
//uword LIdx = min(allidx);
uword RIdx = max(allidx);
// get unique allidx
vec Bool = zeros<vec>(RIdx+1u);
MatAssign(Bool, allidx, 1);
allidx = find(Bool);
/*% get peak and baseline parameters, and spectrum data
[C, FC, NA, Sr, T, Ta, UT, UTa, BW, RB, RD] = dmps(sn, ...
'Centroid', 'FWHM_Ch', 'NetArea', 'StepRatio', 'Tail', 'TailAlpha', ...
'UpperTail', 'UpperTailAlpha', 'BWWidthChan', ...
'RecoilBetaChan', 'RecoilDeltaChan');*/
QStringList Opt;
Opt << "Centroid" << "FWHM_Ch" << "NetArea" << "StepRatio" << "Tail" << "TailAlpha"
<< "UpperTail" << "UpperTailAlpha" << "BWWidthChan" << "RecoilBetaChan" << "RecoilDeltaChan";
field<vec> t_p = dmps(Opt);
vec& C = t_p(0);
vec& FC = t_p(1);
vec& NA = t_p(2);
vec& Sr = t_p(3);
vec& T = t_p(4);
vec& Ta = t_p(5);
vec& UT = t_p(6);
vec& UTa = t_p(7);
vec& BW = t_p(8);
vec& RB = t_p(9);
vec& RD = t_p(10);
//[S, rg, cx, cy, cdy] = dmspec(sn, 'ShortSpectrum', 'AnalysisRange', 'XControl', 'YControl', 'YSlope');
Opt.clear();
Opt << "ShortSpectrum" << "AnalysisRange" << "XControl" << "YControl" << "YSlope";
field<rowvec> t_s = dmspec(Opt);
rowvec& S = t_s(0);
rowvec& rg = t_s(1);
rowvec& cx = t_s(2);
rowvec cy = t_s(3);
rowvec cdy = t_s(4);
// get fitting region
rowvec x;
if( rg.size() != 2 ) //if length(rg) ~= 2
{
qDebug() << "AnalysisRange not an interval, fitting peaks might result in error";
x = RangeVec2(1, S.size()); // x = [1 : length(S)];
}
else
{
x = RangeVec2(rg(0), rg(1)); // x = [rg(1) : rg(2)];
}
rowvec t_x = x - 1;
rowvec y = IdxMat(S, t_x);
/* % fit parameters
[s, CXo, cy, cdy, Ao, Co, Fo, Sro, To, Tao, UTo, UTao, BWo, RBo, RDo, X2] = ...
fitFunction(x, y, 'wrapStepRatio', ...
cx, [], cy, CPf, cdy, [], NA, Af, C, Cf, FC, Ff, Sr, SRf, ...
T, [], Ta, [], UT, [], UTa, [], BW, [], RB, [], RD, [], ...
'Restrict', addArg{:});*/
double X2;
uvec ep;
field<vec> Ps;
field<uvec> free_idx;
Ps << vectorise(cx) << vectorise(cy) << vectorise(cdy) << NA << C << FC << Sr << T << Ta << UT << UTa << BW << RB << RD;
free_idx << ep << CPf << ep << Af << Cf << Ff << SRf << ep << ep << ep << ep << ep << ep << ep;
bool s = Independ::fitFunction(Ps, X2, free_idx, vectorise(x), vectorise(y), "wrapStepRatio", QStringList("Restrict"), addArg);
//mat& CXo = Ps(0);
cy = Ps(1).t();
cdy = Ps(2).t();
mat& Ao = Ps(3);
mat& Co = Ps(4);
mat& Fo = Ps(5);
mat& Sro = Ps(6);
//mat& To = Ps(7);
//mat& Tao = Ps(8);
//mat& UTo = Ps(9);
//mat& UTao = Ps(10);
//mat& BWo = Ps(11);
//mat& RBo = Ps(12);
//mat& RDo = Ps(13);
// only write to temp pat if fitting converged
if(s)
{
tmpbase(cx, cy, cdy); //tmpbase(sn, cx, cy, cdy);
//("cy.txt", cy);
//("cdy.txt", cdy);
/*patSetPeaks(allidx, 'FittingMethod', 'F', 'NetArea', Ao(allidx), ...
'Centroid', Co(allidx), 'FWHM_Ch', Fo(allidx), ...
'NetAreaFree', Af, 'FWHMFree', Ff, 'CentroidFree', Cf);*/
QStringList Descs;
Descs << "NetArea" << "Centroid" << "FWHM_Ch";// << "NetAreaFree" << "FWHMFree" << "CentroidFree";
field<vec> Vals;
Vals << vectorise(Ao(allidx)) << vectorise(Co(allidx)) << vectorise(Fo(allidx));
//<< UMatToMat(Af) << UMatToMat(Ff) << UMatToMat(Cf);
patSetPeaks(allidx, Descs, Vals);
Vals << vectorise(Sro(SRf));
patSetPeaks(SRf, QStringList("StepRatio"), Vals);
vec sqi, uqi;
patBaseVar(sqi, uqi, allidx);
}
// return only the values of free parameters, if any are requested
//if nargin > 1
//Ao = Ao(Af);
//Co = Co(Cf);
//Fo = Fo(Ff);
//Sro = Sro(SRf);
//end
}
void GammaAnalyALG::patSetPeaks(uvec idx, QStringList Descs, field<vec> Values, bool bAll)
{
// get spectrum number and peak indizes, expand idx 'all' to indizes for all peaks
if(idx.is_empty() && bAll) //if isstr(idx) & strcmp(idx, 'all')
{
idx = RangeVec(0, pat.Peaks.n_rows); //idx = dmps(sn, 'Index');
}
// list of possible descriptors and their number in the list
// list must be unique, variable numbers must match
/*QString Descts[40] = {"FittingMethod", "Centroid", "FWHM_Ch", "NetArea", "AreaError", "StepRatio", "Sensitivity",
"CentroidMethod", "NetAreaMethod", "FWHMMethod", "NuclideSoft", "CentroidFree", "NetAreaFree",
"FWHMFree", "Spurious", "Reviewed", "MeanBackCount", "Significance", "ManualTail", "TailChanged",
"ResChanged", "BWWidth", "LineEnergy", "LineEnergyErr", "Yield", "YieldErr", "ReferenceID",
"SumPeakArea", "SumPeakAreaErrP", "NuclideHard", "Nuclide", "BWWidthErr", "CCF", "CCFErr",
"RefEfficiencyErr", "TailAlpha", "UpperTail", "UpperTailAlpha", "RecoilBeta", "RecoilDeltaE"};*/
//int nold = pat.Peaks.n_rows;
// parse descriptor-value pairs, write values to numbered positions in Value cell
uword i=0;
foreach(QString item, Descs) //for i= 3 : 2 : nargin-1
{
if(item == "Centroid")
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_Centroid)(idx(ii)) = Values(i)(ii);
}
}
else if(item == "FWHM_Ch")
{
vec FC = Values(i);
if(any(FC < 0))
{
FC = abs(FC);
qDebug() << "Negative FWHM forbidden, using absolute value";
}
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_FWHM_Ch)(idx(ii)) = FC(ii);
}
}
else if(item == "StepRatio")
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_StepRatio)(idx(ii)) = Values(i)(ii);
}
}
else if(item == "NetArea")
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_NetArea)(idx(ii)) = Values(i)(ii);
}
/* reset area error, if not specified
if isempty(Value{AreaError_n})
pat(idx, 13) = NaN;
end*/
}
else if(item == "AreaError")
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_AreaError)(idx(ii)) = Values(i)(ii);
}
}
else if(item == "MeanBackCount")
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_MeanBackCount)(idx(ii)) = Values(i)(ii);
}
}
else if(item == "Sensitivity")
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_Sensitivity)(idx(ii)) = Values(i)(ii);
}
}
else if(item == "Significance")
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_Significance)(idx(ii)) = Values(i)(ii);
}
}
else if(item == "BWWidth")
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_BWWidth)(idx(ii)) = Values(i)(ii);
}
}
else if(item == "BWWidthErr")
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_BWWidthErr)(idx(ii)) = Values(i)(ii);
}
}
else if(item == "CCF")
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_CCF)(idx(ii)) = Values(i)(ii);
}
}
else if(item == "CCFErr")
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_CCFErr)(idx(ii)) = Values(i)(ii);
}
}
else if(item == "LineEnergy")
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_LineEnergy)(idx(ii)) = Values(i)(ii);
}
}
else if(item == "LineEnergyErr")
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_LineEnergyErr)(idx(ii)) = Values(i)(ii);
}
}
else if(item == "Yield")
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_Yield)(idx(ii)) = Values(i)(ii);
}
}
else if(item == "YieldErr")
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_YieldErr)(idx(ii)) = Values(i)(ii);
}
}
else if(item == "ReferenceID")
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_ReferenceID)(idx(ii)) = Values(i)(ii);
}
}
else if(item == "SumPeakArea")
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_SumPeakArea)(idx(ii)) = Values(i)(ii);
}
}
else if(item == "SumPeakAreaErrP")
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_SumPeakAreaErr)(idx(ii)) = Values(i)(ii);
}
}
else if(item == "RefEfficiencyErr")
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_RefEfficiencyErr)(idx(ii)) = Values(i)(ii);
}
}
else if(item == "ManualTail")
{
/*% if tail parameter changes, raise a flag
Tl = Value{ManualTail_n};
if ~isempty(Tl)
oldTl = pat(idx, 10);
pat(idx, 10) = Tl;
subIdx = find(Tl ~= oldTl);
if any(subIdx)
flags(idx(subIdx), 28) = 'C';
end
end*/
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_Tail)(idx(ii)) = Values(i)(ii);
}
}
else if(item == "TailAlpha")
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_TailAlpha)(idx(ii)) = Values(i)(ii);
}
}
else if(item == "UpperTail")
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_UpperTail)(idx(ii)) = Values(i)(ii);
}
}
else if(item == "UpperTailAlpha")
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_UpperTailAlpha)(idx(ii)) = Values(i)(ii);
}
}
else if(item == "RecoilBeta")
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_RecoilBeta)(idx(ii)) = Values(i)(ii);
}
}
else if(item == "RecoilDeltaE")
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_RecoilDeltaE)(idx(ii)) = Values(i)(ii);
}
}
++i;
}
// reset area error, if not specified
if( Descs.contains("NetArea") && !Descs.contains("AreaError") )//if isempty(Value{AreaError_n})
{
for(uword ii=0; ii<idx.size(); ++ii)
{
pat.Peaks.col(i_AreaError)(idx(ii)) = gNaN;
}
}
}
void GammaAnalyALG::FitBase(uvec cidx, uvec sidx)
{
// get spectrum
/* [R, baseWeights, Step, rg, cx, cy, cdy] = dmspec(sn, ...
'ROISmooth', 'BaseFittingWeights', 'Steps', 'AnalysisRange', 'XControl', 'YControl', 'YSlope');*/
QStringList descs;
descs << "ROISmooth" << "BaseFittingWeights" << "Steps" << "AnalysisRange" << "XControl" << "YControl" << "YSlope";
field<rowvec> t_f = dmspec(descs);
rowvec &R = t_f(0);
rowvec &baseWeights = t_f(1);
rowvec &Step = t_f(2);
rowvec &rg = t_f(3);
rowvec cx = t_f(4);
rowvec cy = t_f(5);
rowvec cdy = t_f(6);
// fit all control points, if no index given
if(cidx.is_empty())
{
cidx = RangeVec(0, cy.size()-1).t();
}
/*if(sidx.is_empty())
{
sidx = [];
}*/
// get interesting channels
rowvec x;
if(rg.size() != 2) //if length(rg) ~= 2
{
qDebug() << "AnalysisRange not an interval";
x = RangeVec2(1, R.size()); //x = [1 : length(R)];
} else {
x = RangeVec2(rg(0), rg(1)); //x = [rg(1) : rg(2)];
}
rowvec t_x = x - 1;
rowvec y = IdxMat(R,t_x);
rowvec w = IdxMat(baseWeights,t_x);
uvec idx = find_nonfinite( cdy(sidx) );
if( any(idx) )
{
MatAssign(cdy, sidx(idx), 0);
}
// get stepfree control points
rowvec t_cx = cx - 1;
cy = cy - IdxMat(Step, t_cx);
// fit control points to residual
/*[s, cx, cy, cdy] = fitFunction(x, y, 'wrapBaseSpline', ...
cx, [], cy, cidx, cdy, sidx, ...
'Weights', w, 'Waitbar', 'Baseline');*/
mat t_cidx = UMatToMat(cidx);
mat t_sidx = UMatToMat(sidx);
double X2;
field<vec> Ps;
field<uvec> free_idx;
Ps << vectorise(cx) << vectorise(cy) << vectorise(cdy);
free_idx << uvec() << cidx << sidx;
QStringList Opts; Opts << "Weights" << "Waitbar";
Independ::fitFunction(Ps, X2, free_idx, vectorise(x), vectorise(y), "wrapBaseSpline", Opts, "Baseline", w);
cx = Ps(0).t();
cy = Ps(1).t();
cdy = Ps(2).t();
// add step back on
t_cx = cx-1;
cy = vectorise(cy).t() + IdxMat(Step,t_cx);
//("FitBase_cx.txt", cx);
//("FitBase_cy.txt", cy);
//("FitBase_cdy.txt", cdy);
// create temporary baseline
tmpbase(cx, cy, cdy);
}
/* =====================================================谱基线结构======================================================== */
void GammaAnalyALG::msrd(mat &d, int &nzero, mat Y, mat F)
{
// convert to column vectors
vec yt = vectorise(Y);
vec ft = vectorise(F);
// check parameter size
int ntotal = yt.size();
if (ft.size() != ntotal)
{
qDebug() << "Sizes of parameters must match";
}
// find "non-zero" data points
uvec idx = find(abs(yt) > eps(1.0));
if(idx.is_empty()) // isempty(idx)
{
d << gNaN;
nzero = ntotal;
}
else {
int n = idx.size();
vec y = yt(idx);
vec r = ft(idx) - y;
// calculate the msrd
d = ( (r.t() / abs(y).t() ) * r) / n;
nzero = ntotal - n;
}
}
mat GammaAnalyALG::lSmoothen(rowvec y)
{
// filter half width in units of FWHM
double f = 1.5;
rowvec s = y;
// get resolution in channel units
int nchan = y.size();
rowvec c = RangeVec2(1, nchan); //c = cumsum(ones(size(y)));
rowvec e = calValues(Cal_Energy, c);
rowvec de = calDerivative(Cal_Energy, c);
rowvec r = calValues(Cal_Resolution, e);
// filter half width: factor times FWHM/channels
rowvec n = max(2, ceil(f * r / de));
// find positions where half width changes
uvec crossover = find(diff(n) != 0) + 1u; //crossover = find(diff(n) ~= 0);
crossover.resize(crossover.size() + 1u);
crossover( crossover.size() ) = nchan; //crossover(end+1) = nchan;
// filter piecewise
rowvec smooth;
int chlow = 1, chhigh, hw, dlow, dhigh, width;
for(int i=0; i<crossover.size(); ++i) //for i = 1 : length(crossover)
{
// next upper border
chhigh = crossover(i);
hw = n(chhigh-1);
// use true data instead of transient at borders
dlow = min(hw, chlow-1);
dhigh = min(hw, nchan - chhigh);
width = 2*hw+1;
if( (chhigh+dhigh - (chlow-dlow)) < width )
{
if(nchan < width)
{
// don't smoothen, but warn, should hardly ever happen
// smooth = y(chlow-dlow : chhigh+dhigh);
smooth = y.cols(chlow-dlow-1, chhigh+dhigh-1);
qDebug() << QString("Can't smoothen: spectrum length: %1, frame width: %2").arg(nchan).arg(width);
}
else // increase data width
{
while( (chhigh+dhigh - (chlow-dlow)) < width )
{
if(dhigh < nchan - chhigh) dhigh = dhigh+1;
if(dlow < chlow - 1) dlow = dlow+1;
}
// smoothen with wider data
smooth = sgderive(y.cols(chlow-dlow-1, chhigh+dhigh-1), width, 0);
}
}
else // smoothen
{
smooth = sgderive(y.cols(chlow-dlow-1, chhigh+dhigh-1), width, 0);
}
// take component where filter width constant
// s(chlow:chhigh) = smooth(1+dlow : end-dhigh);
s.cols(chlow-1,chhigh-1) = smooth.cols(dlow, smooth.size()-dhigh-1);
// next region
chlow = chhigh;
}
return s;
}
mat GammaAnalyALG::sgderive(mat y, int np, int order)
{
mat s = zeros(size(y));
y = vectorise(y); //y = y(:);
// calculate filter
rowvec f;
int k;
sgfilter(f, k, np, order); //[f, k] = sgfilter(np, order);
//("sgderive_f.txt", f);
// continue data symetrically
// dd = [y(k+1:-1:2); y; y(end-1:-1:end-k-1)];
int size_y = y.size();
vec dd = join_vert( join_vert( y(RangeVec(k, 1, -1)), y ),
y( RangeVec(size_y-2, size_y-k-2, -1) ) );
// fold
urowvec idx = RangeVec(0, 2*k); //idx = [0:2*k];
for(uword i=0; i<y.size(); ++i) //for i = 1:length(y)
{
s(i) = as_scalar(f * dd(i + idx)); //s(i) = f*(dd(i+idx));
}
return s;
}
void GammaAnalyALG::sgfilter(rowvec &f, int &k, int frameWidth, int deriveOrder)
{
k = (frameWidth-1) / 2;
rowvec F = RangeVec2(-k, k); //F = [-k : k];
double S0; //S0 = lQS(k, 0);
double S2 = accu( pow( RangeVec2(-k, k), 2) ); //S2 = lQS(k, 2);
double S4 = accu( pow( RangeVec2(-k, k), 4) ); //S4 = lQS(k, 4);
double S6;
double N;
switch(deriveOrder)
{
case 0: // f = ( S2*F2 - S4*F0 ) / ( S2^2 - S4*S0 )
S0 = accu( pow( RangeVec2(-k, k), 0) );
f = S2 * pow(F,2) - S4 * pow(F,0);
N = -( S4*S0 - pow(S2,2) );
f = f/N;
break;
case 2: // 2*f = ( S0*F2 - S2*F0 ) / ( S4*S0 - S2^2 )
S0 = accu( pow( RangeVec2(-k, k), 0) );
f = S0 * pow(F,2) - S2 * pow(F,0);
N = ( S4*S0 - pow(S2,2) )/2;
f = f/N;
break;
case 3: // 6*f = ( S2*F3 - S4*F1 ) / ( S6*S2 - S4^2 )
S6 = accu( pow( RangeVec2(-k, k), 6) ); //S6 = lQS(k, 6);
f = S2 * pow(F,3) - S4 * F;
N = ( S6*S2 - pow(S4,2) )/6;
f = f/N;
break;
default: qDebug() << "derivation Order not implemented yet";
}
}
mat GammaAnalyALG::lSmoothROIs(rowvec y, mat c, mat w)
{
// roi width in FWHM units
int f_roi = 2;
// smoothing filter width (9-point smoothing)
int NSmooth = 19;
//NMean = 11;
int n = y.size();
irowvec bool_smooth = zeros<irowvec>(n); //bool_smooth = false(size(y));
mat low_roi = max(1, floor( c - f_roi * w ) );
mat high_roi = min(n, ceil( c + f_roi * w ) );
//("lSmoothROIs_low_roi.txt", low_roi);
//("lSmoothROIs_high_roi.txt", high_roi);
for(int i=0; i<low_roi.size(); ++i)
{
//bool_smooth(low_roi(i):high_roi(i)) = true;
bool_smooth.cols(low_roi(i)-1, high_roi(i)-1)
= ones<irowvec>(high_roi(i)-low_roi(i)+1);
}
bool_smooth(0) = false; // bool_smooth([1,end]) = false;
bool_smooth(bool_smooth.size()-1) = false;
irowvec d_smooth = diff(bool_smooth);
uvec low_smooth = find(d_smooth == 1);
uvec high_smooth = find(d_smooth == -1);
//ss = sgderive(y(~bool_smooth), NMean, 0);
rowvec smoothSpec = sgderive(y, NSmooth, 0);
//("lSmoothROIs_smoothSpec.txt", smoothSpec);
rowvec s = y;
int skip = 0;
for(int i=0; i<low_smooth.size(); ++i) // i = 1 : length(low_smooth)
{
int dx = high_smooth(i) - low_smooth(i) + 1;
int l = max(1, low_smooth(i)-skip);
skip = skip + dx;
rowvec localSmoothData = smoothSpec.cols(low_smooth(i), high_smooth(i));
s.cols(low_smooth(i), high_smooth(i)) = localSmoothData;// - mean(localSmoothData) + ss(l);
}
return s;
}
field<rowvec> GammaAnalyALG::dmspec(QStringList items)
{
field<rowvec> Results(items.size());
int i=0;
rowvec BaseChan = cumsum(ones<rowvec>(m_nChans));
/*[pC, pNA, pFC, pSt, pT, pTa, pUT, pUTa, pBWW, pRBC, pRDC] = dmps(sn, ...
'Centroid', 'NetArea', 'FWHM_Ch', 'Step', ...
'Tail', 'TailAlpha', 'UpperTail', 'UpperTailAlpha', 'BWWidthChan', ...
'RecoilBetaChan', 'RecoilDeltaChan');*/
vec pC, pNA, pFC, pSt, pT, pTa, pUT, pUTa, pBWW, pRBC, pRDC, z, o;
QStringList Opts;
if(items.contains("BaseLine") || items.contains("Approximation") || items.contains("Stripped")
|| items.contains("Steps") || items.contains("ROISmooth") || items.contains("BaseFittingWeights"))
{
Opts << "Centroid" << "NetArea" << "FWHM_Ch" << "Step" << "Tail" << "TailAlpha" << "UpperTail"
<< "UpperTailAlpha" << "BWWidthChan" << "RecoilBetaChan" << "RecoilDeltaChan";
}
if(!Opts.isEmpty())
{
field<vec> f = dmps(Opts);
pC = f(0); pNA = f(1); pFC = f(2); pSt = f(3); pT = f(4); pTa = f(5);
pUT = f(6); pUTa = f(7); pBWW = f(8); pRBC = f(9); pRDC = f(10);
z = zeros(size(pC));
o = ones(size(pC));
}
if(items.contains("Approximation") || items.contains("Residual") || items.contains("ChiSquare"))
{
PiecePoly bpp = Independ::makeBasePP(vectorise(baseInfo(s_XControl)), vectorise(baseInfo(s_YControl)),
vectorise(baseInfo(s_YSlope)), pC, pFC, pSt);
baseInfo(s_Approximation) = Independ::peakBaseSpline(BaseChan, bpp, pNA, pC, pFC, pSt, pT, pTa, pUT, pUTa, pBWW, pRBC, pRDC);
}
if(items.contains("Stripped") || items.contains("SmoothStripped") || items.contains("ROISmooth"))
{
if(!pC.is_empty())
{
PiecePoly bpp;
rowvec pk = Independ::peakBaseSpline(BaseChan, bpp, pNA, pC, pFC, z, pT, pTa, pUT, pUTa, pBWW, pRBC, pRDC);
// PrintMat22("BaseChan.txt", BaseChan);
// PrintMat22("pNA.txt", pNA);
// PrintMat22("pC.txt", pC);
// PrintMat22("pFC.txt", pFC);
// PrintMat22("z.txt", z);
// PrintMat22("pT.txt", pT);
// PrintMat22("pTa.txt", pTa);
// PrintMat22("pUT.txt", pUT);
// PrintMat22("pUTa.txt", pUTa);
// PrintMat22("pBWW.txt", pBWW);
// PrintMat22("pRBC.txt", pRBC);
// PrintMat22("pRDC.txt", pRDC);
// PrintMat22("baseInfo(s_Spectrum).txt", baseInfo(s_Spectrum));
// PrintMat22("pk.txt", pk);
baseInfo(s_Stripped) = baseInfo(s_Spectrum) - pk;
}
else baseInfo(s_Stripped) = baseInfo(s_Spectrum);
}
if(items.contains("Steps") || items.contains("ROISmooth"))
{
baseInfo(s_Steps) = Independ::peakSteps(BaseChan, pC, pFC, pSt);
}
foreach (QString item, items)
{
if(item == "Spectrum" || item == "LongSpectrum" || item == "ShortSpectrum")
{
Results(i) = baseInfo(s_Spectrum);
}
else if(item == "XControl")
{
Results(i) = baseInfo(s_XControl);
}
else if(item == "YControl")
{
Results(i) = baseInfo(s_YControl);
}
else if(item == "YSlope")
{
Results(i) = baseInfo(s_YSlope);
}
else if(item == "AnalysisRange")
{
Results(i) = baseInfo(s_AnalysisRange);
}
else if(item == "AnalysisRangeSave")
{
if(baseInfo(s_AnalysisRange).is_empty())
{
rowvec r; r << 1 << m_nChans;
Results(i) = r;
}
else {
Results(i) = baseInfo(s_AnalysisRange);
}
}
else if(item == "BaseLine")
{
if(pat.Peaks.n_cols == MAX_INDEX_PEAK)
{
PiecePoly bpp = Independ::makeBasePP(vectorise(baseInfo(s_XControl)), vectorise(baseInfo(s_YControl)),
vectorise(baseInfo(s_YSlope)), pC, pFC, pSt);
Results(i) = Matlab::ppval(BaseChan, bpp)+Independ::peakSteps(BaseChan, pC, pFC, pSt);
baseInfo(s_BaseLine) = Results(i);
}
}
else if(item == "BaseFittingWeights")
{
//BaseFittingWeights = lGetWeights(BaseChan);
// get peak data
//[C, NA, FC, St] = dmps(sn, 'Centroid', 'NetArea', 'FWHM_Ch', 'Step');
// get fitting weights: punish regions near peaks
mat w = ones(size(BaseChan));
for(int i=0; i<pC.size(); ++i) // for i = 1 : length(C)
{
if(pNA(i) > 1) // if NA(i) > 1
{
// w = w + sqrt(NA(i)) ./ ((max(1, abs((x-C(i))/FC(i)))).^3 );
w = w + sqrt(pNA(i)) / pow( max(1, abs((BaseChan-pC(i))/pFC(i))), 3);
}
}
Results(i) = 1 / (abs(w) + 1);
baseInfo(s_BaseFittingWeights) = Results(i);
}
else if(item == "Approximation")
{
Results(i) = baseInfo(s_Approximation);
}
else if(item == "Residual")
{
Results(i) = baseInfo(s_Spectrum) - baseInfo(s_Approximation);
baseInfo(s_Residual) = Results(i);
}
else if(item == "Stripped")
{
Results(i) = baseInfo(s_Stripped);
}
else if(item == "PSS")
{
rowvec e = calValues(Cal_Energy, BaseChan);
rowvec dE = calDerivative(Cal_Energy, BaseChan);
rowvec r = calValues(Cal_Resolution, e);
Results(i) = -pksens(baseInfo(s_Spectrum), r/dE, BaseChan, BaseChan);
baseInfo(s_PSS) = Results(i);
}
else if(item == "SmoothStripped")
{
Results(i) = lSmoothen(baseInfo(s_Stripped));
baseInfo(s_SmoothStripped) = Results(i);
}
else if(item == "ChiSquare")
{
if(baseInfo(s_AnalysisRange).is_empty()) // isempty(AnalysisRange)
{
Results(i) << gNaN << gNaN; // ChiSquare = [NaN, NaN];
}
else {
//ch = AnalysisRange(1) : AnalysisRange(2);
rowvec F = baseInfo(s_Approximation).cols(baseInfo(s_AnalysisRange)(0), baseInfo(s_AnalysisRange)(1));
rowvec S = baseInfo(s_Spectrum).cols(baseInfo(s_AnalysisRange)(0), baseInfo(s_AnalysisRange)(1));
mat MSRD;
int nzero;
msrd(MSRD, nzero, S, F); // [MSRD, nzero] = msrd(S, F);
// return ChiSquare and number of channels with zero counts
rowvec tmp; tmp << nzero;
Results(i) = join_horiz(MSRD, tmp); //ChiSquare=[ MSRD, nzero];
}
baseInfo(s_ChiSquare) = Results(i);
}
/*else if(item == "ArtXControl")
{
Results(i) = baseInfo(s_ArtXControl);
}
else if(item == "ArtYControl")
{
Results(i) = baseInfo(s_ArtYControl);
}
else if(item == "ArtYSlope")
{
Results(i) = baseInfo(s_ArtYSlope);
}
else if(item == "ArtificialBase")
{
if(baseInfo(s_ArtXControl).is_empty())
{
Results(i) = zeros(size(baseInfo(s_Spectrum)));
}
else //if ~isempty(artbasestruc)
{
mat empty = mat(0, 0);
PiecePoly ArtBaselinePP = makeBasePP(baseInfo(s_ArtXControl), baseInfo(s_ArtYControl),
baseInfo(s_ArtYSlope), empty, empty, empty);
Results(i) = peakBaseSpline(BaseChan, ArtBaselinePP, empty, empty, empty, empty,
empty, empty, empty, empty, empty, empty, empty);
}
}*/
else if(item == "Steps")
{
Results(i) = baseInfo(s_Steps);
}
/*else if(item == "SteplessBase")
{
Results(i) = baseInfo(s_Stripped) - baseInfo(s_Steps);
}*/
else if(item == "ROISmooth")
{
Results(i) = lSmoothROIs(baseInfo(s_Stripped) - baseInfo(s_Steps), pC, pFC);
baseInfo(s_ROISmooth) = Results(i);
}
++i;
}
/*for(int i=0; i<Results.size(); ++i)
{
QString name = QString("dmspec_%1.txt").arg(items.at(i));
//(name, Results(i));
}*/
return Results;
}
/* ====================================================PAT结构======================================================== */
mat GammaAnalyALG::calValuesByName(CalibrationType cal, mat x)
{
mat y = zeros(size(x));
// get named calibration parameter
vec p;
if(getCalibrationPara(p, cal)) //[s, msg, p] = getCalibrationPara(p, cal, name);
{
// evaluate calibration function
y = Independ::calFcnEval(x, p); //[y, s] = calFcnEval(x, p);
}
return y;
}
mat GammaAnalyALG::calDerivativeByName(CalibrationType cal, mat x)
{
mat dy = zeros(size(x));
// get named calibration parameter
colvec p;
if(getCalibrationPara(p, cal)) //[s, msg, p] = getCalibrationPara(sn, cal, name);
{
// evaluate slope
dy = Independ::calDerivEval(x, p); //[dy, s] = calDerivEval(x, p);
}
return dy;
}
void GammaAnalyALG::dmps()
{
QStringList items;
items << "Energy" << "FWHM" << "Multiplet" << "AreaError" << "Significance" << "Tail"
<< "TailAlpha" << "UpperTail" << "UpperTailAlpha" << "StepRatio" << "Efficiency";
field<vec> results = dmps(items);
//results.print("results");
pat.Peaks.col(i_Energy) = results(0);
pat.Peaks.col(i_FWHM) = results(1);
pat.Peaks.col(i_Multiplet) = results(2);
pat.Peaks.col(i_AreaError) = results(3);
pat.Peaks.col(i_Significance) = results(4);
pat.Peaks.col(i_Tail) = results(5);
pat.Peaks.col(i_TailAlpha) = results(6);
pat.Peaks.col(i_UpperTail) = results(7);
pat.Peaks.col(i_UpperTailAlpha) = results(8);
pat.Peaks.col(i_StepRatio) = results(9);
pat.Peaks.col(i_Efficiency) = results(10);
}
vec GammaAnalyALG::dmps(PeakIndex idx)
{
vec vRet;
vec Centroid = pat.Peaks.col(i_Centroid);
if(idx == i_Energy)
{
vRet = pat.Peaks.col(i_Energy); //Energy = EnergyPAT;
uvec Bool = find_nonfinite(vRet);
if( !Bool.is_empty() )
{
//Energy(bool) = calValuesByName(sn, 'Energy', calNames{1}, Centroid(bool));
vRet(Bool) = calValuesByName(Cal_Energy, Centroid(Bool));
}
}
else if(idx == i_Multiplet)
{
vec Multiplet = pat.Peaks.col(i_Multiplet);
uvec Bool = find_nonfinite(Multiplet);
if ( !Bool.is_empty() )
{
vec FWHM_Ch = pat.Peaks.col(i_FWHM_Ch); //FWHM_Ch = FWHM_ChPAT;
uvec Bool2 = find_nonfinite(FWHM_Ch); //bool = isnan(FWHM_Ch);
if( !Bool2.is_empty() ) //if any(bool)
{
vec Energy = pat.Peaks.col(i_Energy); //Energy = EnergyPAT;
uvec Bool3 = find_nonfinite(Energy);
if( !Bool3.is_empty() )
{
Energy(Bool3) = calValuesByName(Cal_Energy, Centroid(Bool3));
}
vec Res = calValuesByName(Cal_Resolution, Energy);
vec dE = calDerivativeByName(Cal_Energy, Centroid);
vec Res_Ch = Res / dE; //Res_Ch=Res./dE;
FWHM_Ch(Bool2) = Res_Ch(Bool2); //FWHM_Ch(bool) = Res_Ch(bool);
}
// where the PAT entry is NaN, set calculated flag; for calculation,
// all peaks are needed !!
vec tmpMultiplet = peakFindMultiplets(Centroid, pat.Peaks.col(i_NetArea), FWHM_Ch, 2);
Multiplet(Bool) = tmpMultiplet(Bool);
}
vRet = Multiplet;
//vRet.print("vRet = ");
}
else if(idx == i_FWHM)
{
vec FWHM = pat.Peaks.col(i_FWHM);
uvec Bool = find_nonfinite(FWHM);
if(!Bool.is_empty())
{
vec FWHM_Ch = pat.Peaks.col(i_FWHM_Ch); //FWHM_Ch = FWHM_ChPAT;
vec dE = calDerivativeByName(Cal_Energy, Centroid);
uvec Bool2 = find_nonfinite(FWHM_Ch); //bool = isnan(FWHM_Ch);
if( !Bool2.is_empty() ) //if any(bool)
{
vec Energy = pat.Peaks.col(i_Energy); //Energy = EnergyPAT;
uvec Bool3 = find_nonfinite(Energy);
if( !Bool3.is_empty() )
{
Energy(Bool3) = calValuesByName(Cal_Energy, Centroid(Bool3));
}
vec Res = calValuesByName(Cal_Resolution, Energy);
vec Res_Ch = Res / dE; //Res_Ch=Res./dE;
FWHM_Ch(Bool2) = Res_Ch(Bool2); //FWHM_Ch(bool) = Res_Ch(bool);
}
FWHM(Bool) = FWHM_Ch(Bool) % dE(Bool);
}
vRet = FWHM;
}
else if(idx == i_AreaError)
{
//AreaError = lCalcAreaError(AreaErrorPAT, NetArea, BackgroundArea, LC);
vec dNA = pat.Peaks.col(i_AreaError);
uvec Bool = find_nonfinite(dNA);
if(!Bool.is_empty())
{
vec NA = pat.Peaks.col(i_NetArea);
vec LC = pat.Peaks.col(i_LC);
vec BA = pat.Peaks.col(i_BackgroundArea);
dNA(Bool) = sqrt( max( NA(Bool), LC(Bool) ) + BA(Bool) );
}
vRet = dNA;
}
else if(idx == i_Significance)
{
//Significance = lCalcSignificance(SignificancePAT, NetArea, LC);
vec signif = pat.Peaks.col(i_Significance);
uvec Bool = find_nonfinite(signif);
if(!Bool.is_empty())
{
vec NA = pat.Peaks.col(i_NetArea);
vec LC = pat.Peaks.col(i_LC);
signif(Bool) = NA(Bool) / LC(Bool);
}
// Do not report significance for neutron lumps
//Significance(Source == 'G') = NaN;
vRet = signif;
}
else if(idx == i_Efficiency)
{
vec vRet = pat.Peaks.col(i_Efficiency);
uvec Bool = find_nonfinite(vRet);
if(!Bool.is_empty())
{
vec Energy = pat.Peaks.col(i_Energy); //Energy = EnergyPAT;
uvec Bool2 = find_nonfinite(Energy);
if( !Bool2.is_empty() )
{
//Energy(bool) = calValuesByName(sn, 'Energy', calNames{1}, Centroid(bool));
Energy(Bool2) = calValuesByName(Cal_Energy, Centroid(Bool2));
}
vRet(Bool) = calValuesByName(Cal_Efficiency, Energy(Bool));
}
}
return vRet;
}
field<vec> GammaAnalyALG::dmps(QStringList items)
{
field<vec> Results(items.size());
vec FWHM_Ch, StepRatio, Res, dE, Res_Ch, Energy, Centroid = pat.Peaks.col(i_Centroid);
bool bFwhmCh = false, bdE = false;
/* 需要 FWHM_Ch 的项FWHM、SigmaCh、UpperTailBeta、Multiplet */
if(items.contains("FWHM_Ch") || items.contains("FWHM") || items.contains("SigmaCh") || items.contains("Multiplet"))
{
bFwhmCh = true;
}
/* 需要 dE 的项Res_Ch、FWHM、UpperTailBeta、BWWidthChan、RecoilBetaChan、RecoilDeltaChan
* 1、Res_ChFWHM_Ch
* 2、FWHM_ChFWHM、SigmaCh、UpperTailBeta、Multiplet*/
if(bFwhmCh || items.contains("BWWidthChan") || items.contains("RecoilBetaChan") || items.contains("RecoilDeltaChan"))
{
bdE = true;
}
/* 需要 Energy 的项Res、Efficiency、EfficiencyErr、Tail、TailAlpha、UpperTail、UpperTailAlpha、StepRatio、LineDeviation
* 1、ResRes_Ch
* 2、Res_ChFWHM_Ch */
if(bdE || items.contains("Energy") || items.contains("Res") || items.contains("Tail")
|| items.contains("TailAlpha") || items.contains("UpperTail") || items.contains("UpperTailAlpha")
|| items.contains("StepRatio") || items.contains("Efficiency"))
{
Energy = pat.Peaks.col(i_Energy); //Energy = EnergyPAT;
uvec Bool = find_nonfinite(Energy);
if( !Bool.is_empty() )
{
//Energy(bool) = calValuesByName(sn, 'Energy', calNames{1}, Centroid(bool));
Energy(Bool) = calValuesByName(Cal_Energy, Centroid(Bool));
}
}
if(bdE)
{
//dE = calDerivativeByName(sn, 'Energy', calNames{1}, Centroid);
dE = calDerivativeByName(Cal_Energy, Centroid);
}
if(bFwhmCh)
{
FWHM_Ch = pat.Peaks.col(i_FWHM_Ch); //FWHM_Ch = FWHM_ChPAT;
// PrintMat22("FWHM_CH.txt", FWHM_Ch);
uvec Bool = find_nonfinite(FWHM_Ch); //bool = isnan(FWHM_Ch);
if( !Bool.is_empty() ) //if any(bool)
{
//Res = calValuesByName(sn, 'Resolution', calNames{2}, Energy);
Res = calValuesByName(Cal_Resolution, Energy);
Res_Ch = Res / dE; //Res_Ch=Res./dE;
FWHM_Ch(Bool) = Res_Ch(Bool); //FWHM_Ch(bool) = Res_Ch(bool);
}
}
if(items.contains("StepRatio") || items.contains("Step"))
{
StepRatio = pat.Peaks.col(i_StepRatio); //StepRatio = StepRatioPAT;
uvec Bool = find_nonfinite(StepRatio); //bool = isnan(StepRatio);
if( !Bool.is_empty() ) //if any(bool)
{
StepRatio(Bool) = calValuesByName(Cal_Step_ratio, Energy(Bool));
}
}
int i=0;
foreach (QString item, items)
{
if(item == "Centroid")
{
Results(i) = Centroid;
}
else if(item == "Energy")
{
Results(i) = Energy;
}
else if(item == "Res")
{
if(Res.is_empty()) Results(i) = calValuesByName(Cal_Resolution, Energy);
else Results(i) = Res;
}
else if(item == "FWHM")
{
vec FWHM = pat.Peaks.col(i_FWHM);
uvec Bool = find_nonfinite(FWHM);
if(!Bool.is_empty())
{
FWHM(Bool) = FWHM_Ch(Bool) % dE(Bool);
}
Results(i) = FWHM;
}
else if(item == "FWHM_Ch")
{
Results(i) = FWHM_Ch;
}
else if(item == "NetArea")
{
Results(i) = pat.Peaks.col(i_NetArea);
}
else if(item == "AreaError")
{
//AreaError = lCalcAreaError(AreaErrorPAT, NetArea, BackgroundArea, LC);
vec dNA = pat.Peaks.col(i_AreaError);
uvec Bool = find_nonfinite(dNA);
if(!Bool.is_empty())
{
vec NA = pat.Peaks.col(i_NetArea);
vec LC = pat.Peaks.col(i_LC);
vec BA = pat.Peaks.col(i_BackgroundArea);
dNA(Bool) = sqrt( max( NA(Bool), LC(Bool) ) + BA(Bool) );
}
Results(i) = dNA;
}
else if(item == "Significance")
{
//Significance = lCalcSignificance(SignificancePAT, NetArea, LC);
vec signif = pat.Peaks.col(i_Significance);
uvec Bool = find_nonfinite(signif);
if(!Bool.is_empty())
{
vec NA = pat.Peaks.col(i_NetArea);
vec LC = pat.Peaks.col(i_LC);
signif(Bool) = NA(Bool) / LC(Bool);
}
// Do not report significance for neutron lumps
//Significance(Source == 'G') = NaN;
Results(i) = signif;
}
else if(item == "Sensitivity")
{
Results(i) = pat.Peaks.col(i_Sensitivity);
}
else if(item == "Multiplet")
{
vec Multiplet = pat.Peaks.col(i_Multiplet);
uvec Bool = find_nonfinite(Multiplet);
if ( !Bool.is_empty() )
{
// where the PAT entry is NaN, set calculated flag; for calculation,
// all peaks are needed !!
vec tmpMultiplet = peakFindMultiplets(Centroid, pat.Peaks.col(i_NetArea), FWHM_Ch, 2);
Multiplet(Bool) = tmpMultiplet(Bool);
}
Results(i) = Multiplet;
}
else if(item == "Index")
{
Results(i) = RangeVec2(1, pat.Peaks.n_rows).t(); //Index = [1 : NPeaks]';
}
else if(item == "NetAreaFree")
{
Results(i) = ones<vec>(Centroid.size());
}
else if(item == "SigmaCh")
{
Results(i) = FWHM_Ch/sqrt(8*log(2)); //SigmaCh = FWHM_Ch/sqrt(8*log(2));
}
else if(item == "Step")
{
Results(i) = pat.Peaks.col(i_NetArea) % StepRatio; //Step = NetArea .* StepRatio;
}
else if(item == "StepRatio")
{
Results(i) = StepRatio;
}
else if(item == "Tail")
{
Results(i) = pat.Peaks.col(i_Tail); //Tail = ManualTail;
uvec Bool = find_nonfinite(Results(i)); //bool = isnan(Tail);
if( !Bool.is_empty() )
{
Results(i)(Bool) = calValuesByName(Cal_Tail, Energy(Bool));
}
}
else if(item == "TailAlpha")
{
Results(i) = pat.Peaks.col(i_TailAlpha); //TailAlpha = TailAlphaPAT;
uvec Bool = find_nonfinite(Results(i));
if( !Bool.is_empty() )
{
Results(i)(Bool) = calValuesByName(Cal_Tail_alpha, Energy(Bool));
}
}
else if(item == "UpperTail")
{
Results(i) = pat.Peaks.col(i_UpperTail); //TailAlpha = TailAlphaPAT;
uvec Bool = find_nonfinite(Results(i));
if( !Bool.is_empty() )
{
Results(i)(Bool) = calValuesByName(Cal_Tail_right, Energy(Bool));
}
}
else if(item == "UpperTailAlpha")
{
Results(i) = pat.Peaks.col(i_UpperTailAlpha); //TailAlpha = TailAlphaPAT;
uvec Bool = find_nonfinite(Results(i));
if( !Bool.is_empty() )
{
Results(i)(Bool) = calValuesByName(Cal_Tail_right_alpha, Energy(Bool));
}
}
else if(item == "BWWidthChan")
{
vec BWWidth = pat.Peaks.col(i_BWWidth);
uvec Bool = find_nonfinite(BWWidth);
if( !Bool.is_empty() )
{
MatAssign(BWWidth, Bool, 0);
}
Results(i) = BWWidth / dE;
}
else if(item == "RecoilBetaChan")
{
Results(i) = pat.Peaks.col(i_RecoilBeta) % dE; //RecoilBetaChan = RecoilBeta .* dE;
}
else if(item == "RecoilDeltaChan")
{
Results(i) = pat.Peaks.col(i_RecoilDeltaE) / dE; //RecoilDeltaChan = RecoilDeltaE ./ dE;
}
else if(item == "Efficiency")
{
vec effi = pat.Peaks.col(i_Efficiency);
uvec Bool = find_nonfinite(effi);
if(!Bool.is_empty())
{
effi(Bool) = calValuesByName(Cal_Efficiency, Energy(Bool));
}
Results(i) = effi;
}
else if(item == "EfficiencyErr")
{
vec eff_err = pat.Peaks.col(i_EfficiencyErr);
uvec Bool = find_nonfinite(eff_err);
if(!Bool.is_empty())
{
//[s, msg, p, perr] = calGetPATPara(sn, 'Efficiency', patName);
vec p;
if( calGetPATPara(p, Cal_Efficiency) )
{
// relative error in percent
//eff_err(Bool) = 100 * calErrorEval(Energy(Bool), p, perr);
}
}
}
++i;
}
return Results;
}
field<vec> GammaAnalyALG::dmps(QList<PeakIndex> list_idx)
{
//pat.Peaks.print("132465");
field<vec> Results(list_idx.size());
vec FWHM_Ch, Res, dE, Res_Ch, Energy, Centroid = pat.Peaks.col(i_Centroid);
bool bFwhmCh = false, bdE = false;
/* 需要 FWHM_Ch 的项FWHM、SigmaCh、UpperTailBeta、Multiplet */
if(list_idx.contains(i_FWHM_Ch) || list_idx.contains(i_FWHM) || list_idx.contains(i_Multiplet))
{
bFwhmCh = true;
bdE = true;
}
/* 需要 Energy 的项Res、Efficiency、EfficiencyErr、Tail、TailAlpha、UpperTail、UpperTailAlpha、StepRatio、LineDeviation
* 1、ResRes_Ch 2、Res_ChFWHM_Ch */
if(bdE || list_idx.contains(i_Energy))
{
Energy = pat.Peaks.col(i_Energy); //Energy = EnergyPAT;
uvec Bool = find_nonfinite(Energy);
if( !Bool.is_empty() )
{
//Energy(bool) = calValuesByName(sn, 'Energy', calNames{1}, Centroid(bool));
Energy(Bool) = calValuesByName(Cal_Energy, Centroid(Bool));
}
//("dmps_Energy.txt", Energy);
}
if(bdE)
{
//dE = calDerivativeByName(sn, 'Energy', calNames{1}, Centroid);
dE = calDerivativeByName(Cal_Energy, Centroid);
//("dmps_dE.txt", dE);
}
if(bFwhmCh)
{
FWHM_Ch = pat.Peaks.col(i_FWHM_Ch); //FWHM_Ch = FWHM_ChPAT;
uvec Bool = find_nonfinite(FWHM_Ch); //bool = isnan(FWHM_Ch);
if( !Bool.is_empty() ) //if any(bool)
{
//Res = calValuesByName(sn, 'Resolution', calNames{2}, Energy);
Res = calValuesByName(Cal_Resolution, Energy);
Res_Ch = Res / dE; //Res_Ch=Res./dE;
FWHM_Ch(Bool) = Res_Ch(Bool); //FWHM_Ch(bool) = Res_Ch(bool);
}
}
int i=0;
foreach (PeakIndex idx, list_idx)
{
if(idx == i_Centroid)
{
Results(i) = Centroid;
}
else if(idx == i_Energy)
{
Results(i) = Energy;
}
else if(idx == i_Multiplet)
{
vec Multiplet = pat.Peaks.col(i_Multiplet);
uvec Bool = find_nonfinite(Multiplet);
if ( !Bool.is_empty() )
{
// where the PAT entry is NaN, set calculated flag; for calculation,
// all peaks are needed !!
vec tmpMultiplet = peakFindMultiplets(Centroid, pat.Peaks.col(i_NetArea), FWHM_Ch, 2);
Multiplet(Bool) = tmpMultiplet(Bool);
}
Results(i) = Multiplet;
}
else if(idx == i_FWHM)
{
vec FWHM = pat.Peaks.col(i_FWHM);
uvec Bool = find_nonfinite(FWHM);
if(!Bool.is_empty())
{
FWHM(Bool) = FWHM_Ch(Bool) % dE(Bool);
}
Results(i) = FWHM;
}
else if(idx == i_NetArea)
{
Results(i) = pat.Peaks.col(i_NetArea);
}
else if(idx == i_AreaError)
{
//AreaError = lCalcAreaError(AreaErrorPAT, NetArea, BackgroundArea, LC);
vec dNA = pat.Peaks.col(i_AreaError);
uvec Bool = find_nonfinite(dNA);
if(!Bool.is_empty())
{
vec NA = pat.Peaks.col(i_NetArea);
vec LC = pat.Peaks.col(i_LC);
vec BA = pat.Peaks.col(i_BackgroundArea);
dNA(Bool) = sqrt( max( NA(Bool), LC(Bool) ) + BA(Bool) );
}
Results(i) = dNA;
}
else if(idx == i_Significance)
{
//Significance = lCalcSignificance(SignificancePAT, NetArea, LC);
vec signif = pat.Peaks.col(i_Significance);
uvec Bool = find_nonfinite(signif);
if(!Bool.is_empty())
{
vec NA = pat.Peaks.col(i_NetArea);
vec LC = pat.Peaks.col(i_LC);
signif(Bool) = NA(Bool) / LC(Bool);
}
// Do not report significance for neutron lumps
//Significance(Source == 'G') = NaN;
Results(i) = signif;
}
else if(idx == i_Sensitivity)
{
Results(i) = pat.Peaks.col(i_Sensitivity);
}
++i;
}
return Results;
}
vec GammaAnalyALG::peakFindMultiplets(vec C, vec NA, vec FC, double RW)
{
uword NPeaks = C.size();
vec Mult = zeros<vec>(NPeaks);
vec CL = floor(C - RW * FC);
vec CR = ceil(C + RW * FC);
for (uword i=0; i<NPeaks; ++i)
{
uvec idx = find(CL <= CR(i) && CR >= CL(i));
if (idx.size() > 1)
{
int L = min(CL(idx));
int R = max(CR(idx));
CL(idx) = ones<vec>(size(idx)) * L; //CL(idx) = L;
CR(idx) = ones<vec>(size(idx)) * R; //CR(idx) = R;
Mult(idx) = ones<vec>(size(idx)); //Mult(idx) = 1;
}
}
return Mult;
}
/* ====================================================刻度更新======================================================== */
void GammaAnalyALG::calUpdate(QChar dataType, vec certEne, bool E1, bool R, bool E2)
{
// if true, peaks found in peak search will be kept
//write_peaks_to_pat = getSpecSetup(sn, 'KeepCalPeakSearchPeaks');
bool write_peaks_to_pat = specSetup.KeepCalPeakSearchPeaks;
// for calibration spectra, use energies from certificate instead of setup file
//[dataType, certEne] = dminfo(sn, 'DataType', 'CertificateGEnergies');
vec ELibEne, ELibRes;
if(dataType == 'C' && !certEne.is_empty())
{
ELibEne = certEne;
ELibRes = certEne;
}
else // getting library energy lines
{
//QString filename = qApp->applicationDirPath() + "/setup/ene_cal_update.txt";
QString filename = m_strFilePath + "/ene_cal_update.txt";
QFile file(filename);
if(file.open(QIODevice::ReadOnly))
{
QTextStream in(&file);
mat tab;
rowvec t_r;
QString line = in.readLine().trimmed();
uword n_row = 0;
while(!in.atEnd())
{
if(!line.isEmpty() && line[0] != '%')
{
QStringList list = line.split(" ", QString::SkipEmptyParts);
if(list.size() >= 5)
{
t_r << list[0].toDouble() << list[1].toDouble() << list[2].toDouble() << list[3].toDouble() << list[4].toDouble();
tab.resize(++n_row, 5u);
tab.row(n_row-1u) = t_r;
}
}
line = in.readLine().trimmed();
}
file.close();
// lines for energy calibration: flag ~= 0
uvec eneIdx = find(tab.col(0));
ELibEne = tab(eneIdx, uvec("1"));
// lines for resolution calibration: flag == 2
uvec resIdx = find(tab.col(0) == 2);
ELibRes = tab(resIdx, uvec("1"));
}
if(ELibEne.is_empty() || ELibRes.is_empty())
{
return;
}
}
////("calUpdate_ELibEne.txt", ELibEne);
////("calUpdate_ELibRes.txt", ELibRes);
if(E1)
{
// first energy calibration update
vec ELibEne1;
uvec ene1PIdx = calPeakSearch(ELibEne1, ELibEne);
calEnergyUpdate1(ene1PIdx, ELibEne1);
}
if(R)
{
// update resolution calibration
vec ELibRes1;
uvec resPIdx = calPeakSearch(ELibRes1, ELibRes);
calResUpdate(resPIdx);
}
if(E2)
{
// second energy calibration update
vec ELibEne2;
uvec ene2PIdx = calPeakSearch(ELibEne2, ELibEne);
calEnergyUpdate2(ene2PIdx, ELibEne2);
}
// keep changes (quickfit) in PAT, if requested
if(write_peaks_to_pat)
{
//patCommit(sn);
//tbaseac(sn);
colvec sqi, uqi;
patBaseVar(sqi, uqi, uvec(), 1);
}
else
{
//patRollback;
pat.Peaks.clear();
pat.Peaks.set_size(0, MAX_INDEX_PEAK);
//tbaseac('reject');
baseInfo(s_XControl).clear();
baseInfo(s_YControl).clear();
baseInfo(s_YSlope).clear();
baseInfo(s_AnalysisRange).clear();
}
}
uvec GammaAnalyALG::calPeakSearch(vec &EF, vec ELib)
{
// max deviation (oldCalibEnergy(PATCentroid) - LibraryLine) < EDevMaxGain*Res
double EDevMaxGain = 0.7;
// first try to find NPeaksDesired peaks, possibly at lower PSS
double NPeaksDesired = 6;
// if not enough peaks, restart at high PSS, and search NPeaksMin peaks
double NPeaksMin = 3;
// Peak search sensitivity
double PSSs = 4;
// get PSS tresholds: high use level, minimum use level
double PSSh = specSetup.CalibrationPSS_high;
double PSSl = specSetup.CalibrationPSS_low;
// get peaks from PAT or do peak search
//[allIdx, C, E, NA, pss, Mult, res] = dmps(sn, 'Index', 'Centroid', 'Energy', 'NetArea', 'Sensitivity', 'Multiplet', 'Res');
QStringList Opts;
Opts << "Index" << "Centroid" << "Energy" << "NetArea" << "Sensitivity" << "Multiplet" << "Res";
field<vec> t_f = dmps(Opts);
uvec allIdx;
vec E = t_f(2), res = t_f(6), pss = t_f(4), Mult = t_f(5);
if (t_f(0).is_empty())
{
vec C, NA, CNL, CNR;
allIdx = PeakSearch(C, NA, CNL, CNR, pss, PSSs, search_full);
E = calValues(Cal_Energy, C);
vec dE = calDerivative(Cal_Energy, C);
res = calValues(Cal_Resolution, E);
vec FC = res / dE;
Mult = peakFindMultiplets(C, NA, FC, 2.25);
}
//("calPeakSearch_E.txt", E);
//("calPeakSearch_res.txt", res);
//("calPeakSearch_pss.txt", pss);
//("calPeakSearch_Mult.txt", Mult);
uvec Bool = zeros<uvec>(E.size());
EF = zeros<vec>(E.size());
uvec idx;
// find out for every singlet peak, if it is close to a library peak
for(int i=0; i<E.size(); ++i) //for i = 1 : length(E)
{
if (!Mult(i))
{
vec EDev = abs(ELib - E(i));
double EDevMin = EDev.min();
if (EDevMin <= EDevMaxGain * res(i))
{
// peak is close, therefor mark it
idx = find(EDev == EDevMin);
Bool(i) = 1u;
EF(i) = ELib(idx(0));
}
}
}
// define steps to reduce PSS level
double PSS = PSSh;
double step = (PSSh - PSSl) / 10;
if (step <= eps(1.0))
{
step = 1.0;
}
// first look, if we can find at least NPeaksDesired peaks at a PSS level
// between PSSh and PSSl
while (PSS >= PSSl)
{
idx = find(Bool && pss > PSS);
if (idx.size() >= NPeaksDesired) break;
PSS = PSS - step;
}
// if failed, loof if we can find at least NPeaksMin peaks at a PSS level
// between PSSh and PSSl
if (idx.size() < NPeaksDesired)
{
PSS = PSSh;
while (PSS >= PSSl)
{
idx = find(Bool && pss > PSS);
if (idx.size() >= NPeaksMin) break;
PSS = PSS - step;
}
}
// reduce vector from all peaks to usable peaks (other entries are 0)
EF = EF(idx);
//EF.print("EF = ");
//idx.print("idx = ");
return idx;
}
void GammaAnalyALG::calEnergyUpdate1(uvec idx, vec ELib)
{
// get pat centroids
vec C = pat.Peaks.col(i_Centroid);
C = C(idx);
// fit parameter
vec p, perr;
bool s = Independ::calFitPara(p, perr, Cal_Energy, C, ELib);
//p.print("p = ");
//perr.print("perr = ");
if (s) // success, set parameter, data points and flag
{
//s = setCalibrationData("CalUpdate1", C, ELib);
//s = setCalibrationPara("CalUpdate1", p, ep);
mat data;
if(lCalDataCheck(data, C, ELib) && Independ::calParaCheck(p))
{
vec F0 = Independ::calFcnEval(C, Acal_energy_para[m_curEner](0));
double sum0 = sum(pow(F0 - ELib, 2)) / C.size();
vec F1 = Independ::calFcnEval(C, p);
double sum1 = sum(pow(F1 - ELib, 2)) / C.size();
if(sum1 < sum0)
{
m_curEner = CalUpdate1;
Cal_Ener_Data[m_curEner] = data;
field<vec> f_p(2);
f_p(0) = p; f_p(1) = perr;
Acal_energy_para[m_curEner] = f_p;
}
}
else { qDebug() << QString("First energy calibration update failed!"); }
}
else { qDebug() << QString("First energy calibration update failed!"); }
}
void GammaAnalyALG::calResUpdate(uvec idx)
{
// Peaks required for fitting
int NPeaksMin = 4;
uword NPeaks = idx.size();
//[NA, C, FC] = dmps(sn, 'NetArea', 'Centroid', 'FWHM_Ch');
QStringList Opts; Opts << "NetArea" << "Centroid" << "FWHM_Ch";
field<vec> t_f = dmps(Opts);
vec NA = t_f(0)(idx);
vec C = t_f(1)(idx);
vec FC = t_f(2)(idx);
//NA.print("NA = ");
//C.print("C = ");
// FC.print("FC = ");
// fit peaks
uvec Bool = zeros<uvec>(NPeaks);
uvec tmp;
vec ret1, ret2, ret3;
double X2;
for(uword i=0; i<NPeaks; ++i)
{
uword j = idx(i);
// quick fixed centroid fit, no multiplets!!
//[bool(i), NA(i), garbage, FC(i), X2] = fitPeakQuick(sn, j, [], j);
tmp << j;
Bool(i) = fitPeakQuick(ret1, ret2, ret3, X2, tmp, uvec(), tmp);
NA(i) = as_scalar(ret1);
FC(i) = as_scalar(ret3);
}
//NA.print("NA = ");
//FC.print("FC = ");
// get peak energies
vec E = calValues(Cal_Energy, C);
vec dE = calDerivative(Cal_Energy, C);
//E.print("E = ");
//dE.print("dE = ");
// fwhm in energy units
vec FE = FC % dE;
//FE.print("FE = ");
idx = find(Bool);
//idx.print("idx = ");
// save data points in any case
//[stat_sv, msg, name] = setCalibrationData(sn, 'Resolution', 'CalUpdate', E(idx), FE(idx), []);
/*bool stat_sv = setCalibrationData("CalUpdate", E(idx), FE(idx));
if (!stat_sv)
{
qDebug() << "Can't set resolution calibration update data in calResUpdate";
}*/
vec p, perr;
bool s;
if (idx.size() > NPeaksMin)
{
s = Independ::calFitPara(p, perr, Cal_Resolution, E(idx), FE(idx)); // fit new resolution
//p.print("p = ");
//perr.print(" perr = ");
} else {
qDebug() << "Not enough peaks for fitting";
s = false;
}
// write results to global
mat data;
if (s && lCalDataCheck(data, E(idx), FE(idx)) && Independ::calParaCheck(p) && lCalCheckPerr(p, perr))
{
//s = setCalibrationPara("CalUpdate", p, perr);
//setCurrentCalibration(Cal_Resolution, pname);
m_curReso = CalResUpdate;
Cal_Reso_Data[m_curReso] = data;
field<vec> f_p(2);
f_p(0) = p; f_p(1) = perr;
Acal_resolution_para[m_curReso] = f_p;
}
else { qDebug() << "Resolution Calibration Update Failed"; }
}
void GammaAnalyALG::calEnergyUpdate2(uvec idx, vec ELib)
{
uword NPeaks = idx.size();
vec C = zeros<vec>(NPeaks);
// fit peaks
vec ret1, ret2, ret3;
uvec tmp;
double X2;
bool s;
for(uword i=0; i<NPeaks; ++i) //for i = 1 : NPeaks
{
uword j = idx(i);
// quick fix width fit
//[s, NA, C(i), FC, X2] = fitPeakQuick(sn, j, j, []);
tmp << j;
s = fitPeakQuick(ret1, ret2, ret3, X2, tmp, tmp, uvec());
C(i) = as_scalar(ret2);
}
//C.print("C = ");
// fit new centroids to library energies
vec p, perr;
s = Independ::calFitPara(p, perr, Cal_Energy, C, ELib);
//p.print("p = ");
//perr.print("perr = ");
// write results to global
if (s)
{
//s = setCalibrationData("CalUpdate2", C, ELib);
//s = setCalibrationPara("CalUpdate2", p, perr);
mat data;
if(lCalDataCheck(data, C, ELib) && Independ::calParaCheck(p) && lCalCheckPerr(p, perr))
{
vec F02 = Independ::calFcnEval(C, Acal_energy_para[m_curEner](0));
double sum02 = sum(pow(F02 - ELib, 2)) / C.size();
vec F2 = Independ::calFcnEval(C, p);
double sum2 = sum(pow(F2 - ELib, 2)) / C.size();
if(sum2 < sum02)
{
m_curEner = CalUpdate2;
Cal_Ener_Data[m_curEner] = data;
field<vec> f_p(2);
f_p(0) = p; f_p(1) = perr;
Acal_energy_para[m_curEner] = f_p;
}
}
else { qDebug() << QString("Calibration update failed in calEnergyUpdate2"); }
}
else // failed
{
qDebug() << "Calibration update failed in calEnergyUpdate2";
}
}
bool GammaAnalyALG::lCalDataCheck(mat &data, mat x, mat y, mat err)
{
bool bRet = false;
uword l = x.size();
if(y.size() != l)
{
return bRet;
}
else {
x = vectorise(x);
y = vectorise(y);
if(err.is_empty())
{
err = gNaN * ones(l, 1);
}
else { err = vectorise(err); }
}
data = join_horiz(x, y);
data.insert_cols(2, err);
if(!data.col(0).is_finite())
{
qDebug() << "Independent variable must be finite";
}
else if (!data.col(1).is_finite())
{
qDebug() << "Dependent variable must be finite";
}
else { bRet = true; }
return bRet;
}
bool GammaAnalyALG::lCalCheckPerr(mat p, mat &perr)
{
bool s = true;
if (perr.is_empty()) perr = gNaN * ones(p.size()-1u, 1u);
else {
if (perr.size() != p.size() - 1u)
{
qDebug() << "error estimate of wrong size";
s = false;
}
else if (any(vectorise(perr) < 0))
{
qDebug() << "negative error estimate";
s = false;
}
}
return s;
}
/*bool GammaAnalyALG::setCalibrationData(QString src, mat x, mat y, mat err)
{
bool bRet = true;
// Check if new data is valid, reshape if necessary
mat data;
bRet = lCalDataCheck(data, x, y, err);
if (!bRet)
{
qDebug() << QString("Invalid calibration data from source %1").arg(src);
return bRet;
}
if(Cal_update_data.size() != 3u) Cal_update_data.set_size(3u);
if(src == "CalUpdate1")
{
Cal_update_data(0u) = data;
}
else if(src == "CalUpdate")
{
Cal_update_data(1u) = data;
}
else if(src == "CalUpdate2")
{
Cal_update_data(2u) = data;
}
else bRet = false;
return bRet;
}
bool GammaAnalyALG::setCalibrationPara(QString src, mat p, mat perr)
{
bool bRet = true;
if(p.is_empty()) return false;
//[s, msg] = lCalParaCheck(cal, p);
bRet = Independ::calParaCheck(p);
if (!bRet) return bRet;
//[s, msg, perr] = lCalCheckPerr(p, perr);
bRet = lCalCheckPerr(p, perr);
if(!bRet) return bRet;
if(Cal_update_para.size() != 3u) Cal_update_para.set_size(3u);
if(src == "CalUpdate1")
{
Cal_update_para(0u) = vectorise(p);
}
else if(src == "CalUpdate")
{
Cal_update_para(1u) = vectorise(p);
}
else if(src == "CalUpdate2")
{
Cal_update_para(2u) = vectorise(p);
}
else bRet = false;
return bRet;
}*/
/*bool GammaAnalyALG::setCurrentCalibration(CalibrationType cal, QString name)
{
// check only one spectrum is affected
bool s = true;
// get index in list of calibration names
vec sqi, uqi;
field<vec> f; f << vec("1");
int n;
switch (cal)
{
case Cal_Energy: n = 1; break;
case Cal_Resolution:
patSetPeaks(uvec(0u), QStringList("ResChanged"), f, true);
n = 2;
break;
case Cal_Efficiency: n = 3; break;
case Cal_Tot_efficiency: n = 4; break;
case Cal_Tail:
patSetPeaks(uvec(0u), QStringList("TailChanged"), f, true);
n = 5;
break;
case Cal_Tail_alpha: n = 6; break;
case Cal_Tail_right: n = 7; break;
case Cal_Tail_right_alpha: n = 8; break;
case Cal_Step_ratio:
patBaseVar(sqi, uqi);
n = 9;
break;
default: s = false; qDebug() << "invalid calibration string in setCurrentCalibration";
}
pat.CalibName[n-1] = name;
}*/
bool GammaAnalyALG::fitPeakQuick(vec &Ao, vec &Co, vec &Fo, double &X2, uvec Af, uvec Cf, uvec Ff, QString vb, double FitWidth)
{
/*% if called with parameter vb, define additional parameters for fitFunction
if strcmpi(vb, 'verbose')
addArg = {'Waitbar', 'peaks above linear baseline'};
else
addArg = {};
end */
bool s = true;
QStringList addArg("");
// get index of leftmost and rightmost peak
//allidx = [Af; Cf; Ff];
uvec allidx = join_vert(Af, Cf);
allidx = join_vert(allidx, Ff);
if (allidx.is_empty())
{
s = false;
Ao.clear();
Co.clear();
Fo.clear();
X2 = gNaN;
return s;
}
uword LIdx = min(allidx);
uword RIdx = max(allidx);
// get unique allidx
uvec Bool = zeros<uvec>(RIdx + 1u);
Bool(allidx) = ones<uvec>(allidx.n_elem);
allidx = find(Bool);
// get peak parameters
//[C, FC, NA] = dmps(sn, 'Centroid', 'FWHM_Ch', 'NetArea');
QStringList Opts; Opts << "Centroid" << "FWHM_Ch" << "NetArea";
field<vec> t_f = dmps(Opts);
vec C = t_f(0);
vec FC = t_f(1);
vec NA = t_f(2);
// interval for fitting: left centroid - width left FWHM to right centroid + ...
int CL = floor(C(LIdx) - FitWidth * FC(LIdx));
int CR = ceil(C(RIdx) + FitWidth * FC(RIdx));
// spectrum and fitting data points (x,y)
rowvec S = baseInfo(s_Spectrum); //S = dmspec(sn, 'Spectrum');
rowvec x = RangeVec2(CL, CR);
rowvec y = S.cols(CL-1, CR-1);
// initial baseline y = k*x + d through extremal points (spectrum data)
double k = (S(CR) - S(CL)) / (CR -CL);
double d = S(CL) - k * CL;
// put seed to too small net areas
double minNA = 0.1;
if (any(NA(Af) < minNA))
{
uvec zeroIdx = find(NA(Af) < minNA);
MatAssign(NA, Af(zeroIdx), minNA);
}
// fit parameters
//[s, BL, Ao, Co, Fo, X2] = fitFunction(x, y, 'peakBaseLin', [k d], [1 2], NA, Af, C, Cf, FC, Ff, addArg{:});
field<vec> para;
field<uvec> free_idx;
vec tmp; tmp << k << d;
para << tmp << NA << C << FC;
free_idx << uvec("0 1") << Af << Cf << Ff;
s = Independ::fitFunction(para, X2, free_idx, vectorise(x), vectorise(y), "peakBaseLin", addArg);
if(para.size() > 3)
{
Ao = para(1);
Co = para(2);
Fo = para(3);
}
// write results to temporary PAT
/*patSetPeaks(allidx, 'FittingMethod', 'Q', 'NetArea', Ao(allidx), ...
'Centroid', Co(allidx), 'FWHM_Ch', Fo(allidx), 'NetAreaFree', Af, ...
'FWHMFree', Ff, 'CentroidFree', Cf);*/
QStringList items;
items << "NetArea" << "Centroid" << "FWHM_Ch";// << "NetAreaFree" << "FWHMFree" << "CentroidFree";
field<vec> vals;
//allidx.print("allidx = ");
//Ao.print("Ao = ");
//Co.print("Co = ");
//Fo.print("Fo = ");
vals << Ao(allidx) << Co(allidx) << Fo(allidx);
patSetPeaks(allidx, items, vals);
vec sqi, uqi;
patBaseVar(sqi, uqi, allidx);
// return only fitted parameters
Ao = Ao(Af);
Co = Co(Cf);
Fo = Fo(Ff);
return s;
}
rowvec GammaAnalyALG::GetFwhmcAll()
{
rowvec Chan = RangeVec2(1, m_nChans);
rowvec E = calValues(Cal_Energy, Chan);
rowvec dE = calDerivative(Cal_Energy, Chan);
rowvec res = calValues(Cal_Resolution, E);
rowvec FwhmC = res / dE;
return FwhmC;
}
stdvec GammaAnalyALG::calculateSCAC(rowvec &spec, rowvec &baseline, rowvec &FwhmC)
{
rowvec scac;
int num = spec.size();
if(num > 0 && baseline.size() == num && FwhmC.size() == num)
{
int iSum1,is;
double s, ds, Sum1, Var1, Var2;
scac = baseline;
/* changed by AWST - the ported routine was insp38, while the correct routine is insp70 */
for (int i = 1; i < num-1; ++i)
{
s = (1.25 * FwhmC(i) - 1) / 2 + 0.0000001;
is = qFloor(s);
ds = s - is;
if ( i > is && ( i + is < num -1 ) && s > 0 )
{
Sum1 = 0;
for (iSum1=i-is; iSum1<=i+is; ++iSum1)
{
Sum1 += spec(iSum1);
}
Var1 = ds * spec(i - is - 1);
Var2 = ds * spec(i + is + 1);
scac(i) = (Sum1 + Var1 + Var2) / (1.25 * FwhmC(i));
}
}
}
return conv_to<stdvec>::from(scac);
}
stdvec GammaAnalyALG::calculateLC(rowvec &baseline, rowvec &FwhmC, double RiskLevelK)
{
rowvec lc;
int num = baseline.size();
if(num > 0 && FwhmC.size() == num)
{
int is;
double s;
/* Initilize the the SCAC array with the Spectrum Array */
lc = baseline;
for (int i=1; i<num-1; ++i)
{
s = (1.25 * FwhmC(i) - 1) / 2 + 0.00000001;
is = qFloor(s);
if (i-is>0 && i+is<num-1)
{
if ( baseline(i) > 0 && FwhmC(i) > 0 ) lc(i) = baseline(i) + RiskLevelK * sqrt(baseline(i) / (1.25 * FwhmC(i)));
else lc(i) = 0;
}
}
}
return conv_to<stdvec>::from(lc);
}
umat GammaAnalyALG::lGetPlotRanges(int n, const vec &c, const vec &fwhmch, const vec &tlo, const vec &thi/*, vec bhi, vec rdelta*/)
{
vec fL = 3 + tlo;
vec fH = 3 + thi;
/*uvec idx = find_finite(bhi);
if (!idx.is_empty())
{
fH(idx) = fH(idx) + 1.5 / (fwhm(idx) % bhi(idx));
fL(idx) = fL(idx) + rdelta(idx) / fwhm(idx);
}*/
vec l = max(1, floor(c - fL % fwhmch));
vec h = min(n, ceil(c + fH % fwhmch));
ivec Bool = zeros<ivec>(n+1);
//l.print("l = ");
//h.print("h = ");
for(uword i=0; i<l.size(); ++i) //for i = 1 : length(l)
{
Bool.rows(l(i)-1, h(i)-1) = ones<ivec>(h(i)-l(i)+1);
}
//Bool([1, end+1]) = false;
Bool(0) = 0;
Bool(n) = 0;
ivec d = diff(Bool);
umat r = join_horiz(find(d == 1), find(d == -1));
return r;
}
uvec GammaAnalyALG::insertpeak(vec iC, PeakInsertMethod src, double minNA)
{
//uword peakNum = phd->vPeak.size();
/*vec pC(peakNum), pFC(peakNum), pT(peakNum), pUT(peakNum), pNA(peakNum), pSt(peakNum), pTa(peakNum), pUTa(peakNum);
vec pBWW(peakNum), pRBC(peakNum), pRDC(peakNum);
pBWW.fill(gNaN);
pRBC.fill(gNaN);
pRDC.fill(gNaN);
uword t_i = 0;
foreach (PeakInfo peak, phd->vPeak)
{
pC(t_i) = peak.peakCentroid;
pFC(t_i) = peak.fwhmc;
pT(t_i) = peak.tail;
pTa(t_i) = peak.tailAlpha;
pUT(t_i) = peak.upperTail;
pUTa(t_i) = peak.upperTailAlpha;
pNA(t_i) = peak.area;
pSt(t_i) = peak.area * peak.stepRatio;
++t_i;
}*/
uvec idx;
rowvec rg = dmspec(QStringList("AnalysisRange"))(0);
//accept only peaks in range
//uvec gidx = find(iC > phd->baseCtrls.rg_low + 1 && iC < phd->baseCtrls.rg_high - 1);
uvec gidx = find(iC > rg(0)+1 && iC < rg(1)-1);
iC = iC(gidx);
// are any peaks to be inserted?
if(iC.is_empty())
{
return idx;
}
/*char src_fl, cent_fl;
switch(src) // peak source flag and centroid source flag
{
case Ins_Manual: src_fl = 'I'; cent_fl = 'U'; break;
case Ins_H0Test:
case Ins_MultFit: src_fl = 'L'; cent_fl = 'L'; break;
case Ins_SumPeakModel: src_fl = 'S'; cent_fl = 'L'; break;
default: qDebug() << "Unkonwn source type"; break;
}*/
// get resolution in channels
/*vec iE = Independ::calFcnEval(Cal_Energy, iC, phd->mapEnerPara[phd->usedEner].p);
vec idE;
Independ::calDerivEval(idE, iC, phd->mapEnerPara[phd->usedEner].p);
vec iR = Independ::calFcnEval(Cal_Resolution, iE, phd->mapResoPara[phd->usedReso].p);*/
vec iE = calValues(Cal_Energy, iC);
vec idE = calDerivative(Cal_Energy, iC);
vec iR = calValues(Cal_Resolution, iE);
vec iRC = iR / idE;
// residual above 0
/*rowvec BaseChan = cumsum(ones<rowvec>(phd->Spec.num_g_channel));
PiecePoly bpp = Independ::makeBasePP(phd->baseCtrls.XCtrl, phd->baseCtrls.YCtrl, phd->baseCtrls.YSlope, pC, pFC, pSt);
rowvec approx = Independ::peakBaseSpline(BaseChan, bpp, pNA, pC, pFC, pSt, pT, pTa, pUT, pUTa, pBWW, pRBC, pRDC);
rowvec resid = spec - approx;*/
rowvec resid = dmspec(QStringList("Residual"))(0);
uvec idx_resid = find(resid<0);
resid(idx_resid) = zeros<rowvec>(idx_resid.size());
// get net area estimate
vec iNA(iC.size());
for(int i=0; i<iC.size(); ++i) //for i=1:length(iC);
{
// ch=[floor(iC(i)-iRC(i)): ceil(iC(i)+iRC(i))];
urowvec ch = RangeVec(floor(iC(i)-iRC(i))-1, ceil(iC(i)+iRC(i))-1);
//iNA(i) = eval('max(minNA,sum(resid(ch)))','minNA');
iNA(i) = max(minNA, sum(resid(ch)));
}
/*idx.set_size(iC.n_elem);
for(int i=0; i<iC.size(); ++i)
{
int j = 0;
while (j < peakNum && iC(i) > phd->vPeak[j].peakCentroid) ++j;
PeakInfo peak;
peak.peakCentroid = iC(i);
peak.area = iNA(i);
peak.sensitivity = gNaN;
phd->vPeak.insert(phd->vPeak.begin() + j, peak);
idx(i) = j;
}*/
// write peaks to temporary pat
/*[idx, oidx] = patAddPeaks(sn, 'Centroid', iC, 'NetArea', iNA, ...
'Source', src_fl, 'CentroidMethod', cent_fl, 'NetAreaMethod', 'S', ...
'Sensitivity', NaN);*/
uvec PropIdx;
PropIdx << i_Centroid << i_NetArea << i_Sensitivity;
mat Datas(iC.size(), 3);
Datas.col(0) = iC;
Datas.col(1) = iNA;
Datas.col(2) = gNaN * ones<vec>(iC.size());
uvec oldIdx;
patAddPeaks(PropIdx, Datas, pat.Peaks, idx, oldIdx);
return idx;
// write peaks to temporary pat
/*[idx, oidx] = patAddPeaks(sn, 'Centroid', iC, 'NetArea', iNA, ...
'Source', src_fl, 'CentroidMethod', cent_fl, 'NetAreaMethod', 'S', ...
'Sensitivity', NaN);*/
}
mat GammaAnalyALG::findPeakRange(uvec pnr, vec C, vec F_Ch, double rg_low, double rg_high)
{
// get all peak centroids and FWHM in channel units
//[C, F_Ch] = dmps(sn, pat, 'Centroid', 'FWHM_Ch');
rowvec ar; //ar = dmspec(sn, 'AnalysisRange');
ar << rg_low << rg_high;
uword n = pnr.size();
mat rg = zeros(n, 2);
vec iL = zeros<vec>(n);
vec iH = zeros<vec>(n);
// left and right borders of multiplet
for(int i=0; i<n; ++i) //for i = 1 : n
{
//[rg(i, 1), rg(i, 2), iL(i), iH(i)] = findVPeakRange(pnr(i), C, F_Ch, ar, varargin{:});
rowvec temp = findVPeakRange(pnr(i), C, F_Ch, ar);
rg(i, 0) = temp(0);
rg(i, 1) = temp(1);
iL(i) = temp(2);
iH(i) = temp(3);
}
return rg;
}
rowvec GammaAnalyALG::findVPeakRange(int pNr, vec C, vec FC, rowvec rg, double LD, double RD)
{
uword NP = C.size();
// find leftmost peak in multiplet, and minimum left border
int i = pNr;
double l = C(i) - LD * FC(i);
while (i > 0 && l < C(i-1) + RD * FC(i-1))
{
i = i - 1;
l = min(l, C(i) - LD * FC(i));
}
// find rightmost peak in multiplet, and maximum right border
int j = pNr;
double h = C(j) + RD * FC(j);
while (j < NP-1 && h > C(j+1) - LD * FC(j+1))
{
j = j + 1;
h = max(h, C(j) + RD * FC(j));
}
// left and right borders of multiplet
l = max(rg(0), floor(l));
h = min(rg(1), ceil(h));
rowvec results;
results << l << h << i << j;
return results;
}
/*void GammaAnalyALG::peakfitdlg(QString action, mat ELim)
{
// write manual changes to PAT
lUseChanges;
// spectrum number and index to regarded peaks
idx = ud.Index;
// unpack fix/free selections
if iscell(NAbool), NAbool=[NAbool{:}]; end
if iscell(Cbool), Cbool=[Cbool{:}]; end
if iscell(Fbool), Fbool=[Fbool{:}]; end
iNA = idx(find(~NAbool)); //net areas not fixed
iCT =idx(find(~Cbool)); //centroid not fixed
iFC =idx(find(~Fbool)); //fwhm not fixed
if(action == "Fit")
{
fitPeakFull(sn, iNA, iCT, iFC, [], 'Verbose');
}
else if(action == "QuickFit")
{
fitPeakQuick(sn, iNA, iCT, iFC, 'Verbose');
}
lPromptAccept( action);
}*/