1721 lines
49 KiB
C++
1721 lines
49 KiB
C++
#ifndef BASEINIT
|
||
#define BASEINIT
|
||
|
||
#endif // BASEINIT
|
||
|
||
#include "GammaAnalyAlgLib.h"
|
||
#include <armadillo>
|
||
#include "fitfunc.h"
|
||
#include <QVector>
|
||
|
||
SpecBaseInfo specInfo;
|
||
PAT pat;
|
||
|
||
void findROI(urowvec& l, urowvec& h, rowvec& roiMask,
|
||
vec fLo = vec("2"), uvec idx = uvec(0), vec fHi = vec(0))
|
||
{
|
||
if(fHi.is_empty()) fHi = fLo;
|
||
|
||
colvec C = pat.Peaks.col(i_Centroid);
|
||
colvec FWHM_CH = pat.Peaks.col(i_FWHM_Ch);
|
||
rowvec rg = specInfo.AnalysisRange;
|
||
|
||
int NChan = specInfo.Spectrum.size();
|
||
if(idx.is_empty()) idx = RangeVec(0, C.size()-1);
|
||
|
||
colvec temp_FWHM = IdxMat(FWHM_CH, idx);
|
||
colvec temp_C = IdxMat(C, idx);
|
||
colvec tempL = fLo % temp_FWHM;
|
||
colvec tempH = fHi % temp_FWHM;
|
||
|
||
colvec pkLo = max(rg[0], floor(temp_C - tempL));
|
||
colvec pkHi = min(rg[1], ceil(temp_C + tempH));
|
||
|
||
roiMask = zeros<rowvec>(NChan);
|
||
|
||
for(size_t i=0; i!=pkLo.size(); ++i)
|
||
{
|
||
roiMask(span(pkLo[i]-1, pkHi[i]-1)) = 1;
|
||
}
|
||
|
||
roiMask[0] = 0;
|
||
roiMask[roiMask.size()-1] = 0;
|
||
|
||
rowvec db = diff(roiMask);
|
||
l = find(db == 1) + 1;
|
||
h = find(db == -1);
|
||
}
|
||
|
||
void 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');
|
||
rowvec stripped = specInfo.Stripped;
|
||
step = specInfo.Steps;
|
||
rg = specInfo.AnalysisRange;
|
||
|
||
// same as matlab [pC, pPSS, pFC] = dmps(sn, 'Centroid', 'Sensitivity', 'FWHM_Ch');
|
||
colvec pC = pat.Peaks.col(i_Centroid);
|
||
colvec pPSS = pat.Peaks.col(i_Sensitivity);
|
||
colvec 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;
|
||
|
||
urowvec l;
|
||
urowvec h;
|
||
|
||
findROI(l, h, roiMask, factFWHM); // l,h,roiMask are out params
|
||
|
||
// get borders of ranges containing big peaks
|
||
rowvec isbig = zeros(size(l));
|
||
uvec idx;
|
||
for(size_t i=0; i<l.size(); ++i)
|
||
{
|
||
idx = find(pC > l[i] && pC < h[i]);
|
||
if(any(pPSS(idx) > PSSBreak)) isbig[i] = 1;
|
||
}
|
||
|
||
rowvec Lr = partVec(l, isbig); // row
|
||
rowvec Hr = partVec(h, isbig); // row
|
||
|
||
// remove ROIs that end to near lower cutoff
|
||
idx = find(Hr <= (rg(0) + 2));
|
||
if(any(idx))
|
||
{
|
||
for(size_t i=0; i!=idx.size(); ++i)
|
||
{
|
||
Lr.shed_col(idx(i));
|
||
Hr.shed_col(idx(i));
|
||
}
|
||
}
|
||
|
||
// shift ROI limit if first ROI goes beyond lower cutoff
|
||
if((!Lr.is_empty()) && (Lr[0] <= rg[0]))
|
||
{
|
||
Lr[0] = rg[0] + 1;
|
||
}
|
||
|
||
// add lower cutoff as a ROI
|
||
urowvec V(1);
|
||
V[0] = 1;
|
||
|
||
Lr.insert_cols(0,V);
|
||
V[0] = rg[0];
|
||
Hr.insert_cols(0, V);
|
||
|
||
L = vectorise(Lr);
|
||
H = vectorise(Hr);
|
||
}
|
||
|
||
rowvec restpolyfit(vec& r, vec x, vec y, vec p0, vec w = vec(0u))
|
||
{
|
||
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);
|
||
rowvec p;
|
||
if(fixIdx.size() == 1)
|
||
{
|
||
vec yfix = V.col( fixIdx(0) ) * p0( fixIdx(0) );
|
||
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();
|
||
}
|
||
else {
|
||
qDebug() << "In function restpolyfit, the finite number of p0 must be 1";
|
||
}
|
||
return p;
|
||
}
|
||
|
||
void lNorm(rowvec yy, rowvec mask, double xc, double xh, double yh, double np,
|
||
double& nrm, double& lbda, double& yc, double& dyc, double& xl, double& yh, double& dyh)
|
||
{
|
||
// check argument, get polynomial mask
|
||
rowvec p;
|
||
if(np == 4)
|
||
{
|
||
p.resize(4);
|
||
p << gNaN << gNaN << gNaN << yh;
|
||
}
|
||
else if(np == 2)
|
||
{
|
||
p.resize(2);
|
||
p << gNaN << yh;
|
||
}
|
||
else qDebug() << "invalid np";
|
||
|
||
// get whole fitting interval
|
||
double k = xh - xc;
|
||
xl = qMax(1, xh-2*k);
|
||
rowvec xx = RangeVec2(xl-1,xh-1);
|
||
|
||
// remove flagged points
|
||
uvec idx = find( IdxMat(mask,xx) == true );
|
||
for(int i=idx.size()-1; i>=0; --i)
|
||
{
|
||
xx.shed_col(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;
|
||
|
||
// perform fit
|
||
colvec r;
|
||
p = restpolyfit(r, x, y, 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 yval, dyval;
|
||
rowvec rg(2);
|
||
rg(0) = -k;
|
||
rg(1) = 0;
|
||
yval = Matlab::polyval(p, rg);
|
||
dyval = Matlab::polyval(dp, rg);
|
||
|
||
yc = yval(0);
|
||
yh = yval(1);
|
||
|
||
dyc = dyval(0);
|
||
dyh = dyval(1);
|
||
}
|
||
|
||
void 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 = qMax(1, x - MinDist);
|
||
xn = qMax(1, hi - 1);
|
||
|
||
// binary search
|
||
// get norms for initial new x
|
||
double norm, lambda, x0, y0, dy0;
|
||
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, y0, dy0);
|
||
}
|
||
|
||
hit = (x0 <= xb);
|
||
}
|
||
|
||
void 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) = [];
|
||
uvec idx = find( mask(xx) == true );
|
||
for(int i=idx.size()-1; i>=0; ++i)
|
||
{
|
||
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));
|
||
yv = IdxMat(data, 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, xx);
|
||
xx = xx - x0;
|
||
colvec r;
|
||
rowvec p = restpolyfit(r, xx, yy, p0);
|
||
|
||
// evaluate
|
||
yv = Matlab::polyval(p, xv - x0);
|
||
|
||
mat a, b;
|
||
Matlab::polyder(a, b, 1, p);
|
||
dyv = Matlab::polyval(a, xv - x0);
|
||
}
|
||
|
||
void 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);
|
||
|
||
double x = rg[1];
|
||
double y = gNaN;
|
||
double dy = gNaN;
|
||
bool pingold = true;
|
||
double nP = 0;
|
||
while(x > rg[0])
|
||
{
|
||
// get next controlpoints
|
||
if(pingold) nP = 2;
|
||
else nP = 4;
|
||
|
||
double xn, yn, dyn, yo, dyo;
|
||
bool pingNew = false;
|
||
lSelectNext(x, y, data, xb, roiMask, nP, xn, yn, dyn, pingNew, yo, 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
|
||
lFitRange(x, gNaN, xb, x, data, roiMask, 2, cy, cdy);
|
||
cx << x;
|
||
}
|
||
else if(roiIdx == 1)
|
||
{
|
||
// lowest controlpoint, must be introduced
|
||
lFitRange(xb, gNaN, x, xb, data, roiMask, 2, yn, dyn);
|
||
rowvec tmp;
|
||
tmp << xb;
|
||
cx = join_horiz(tmp, cx);
|
||
tmp << yn;
|
||
cy = join_horiz(tmp, cy);
|
||
tmp << dyn;
|
||
cdy = join_horiz(tmp, cdy);
|
||
}
|
||
else
|
||
{
|
||
if(delta < N1)
|
||
{
|
||
// no control point introduced
|
||
}
|
||
else
|
||
{
|
||
// introducing only one control point
|
||
double xc = round((x + xb)/2);
|
||
lFitRange(xb,gNaN,x, xc, data, roiMask, 2, yn, dyn);
|
||
rowvec tmp;
|
||
tmp << xc;
|
||
cx = join_horiz(tmp, cx);
|
||
tmp << yn;
|
||
cy = join_horiz(tmp, cy);
|
||
tmp << dyn;
|
||
cdy = join_horiz(tmp, cdy);
|
||
}
|
||
}
|
||
}
|
||
else if(delta < N3)
|
||
{
|
||
rowvec tmp; tmp << xb << x;
|
||
lFitRange(x, y, xb, tmp, data, roiMask, 2, yn, dyn);
|
||
cx = join_horiz(tmp, cx); //cx = [xb, x, cx];
|
||
tmp << yn;
|
||
cy = join_horiz(tmp, 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);
|
||
rowvec tmp;
|
||
tmp << xb << xc << x;
|
||
lFitRange(xb, gNaN, x, tmp, data, roiMask, 4, yn, dyn);
|
||
cx = join_horiz(tmp, cx); //cx = [xb,xc,x,cx];
|
||
tmp << yn;
|
||
cy = join_horiz(tmp, 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
|
||
rowvec tmp; tmp << x;
|
||
cx = join_horiz(tmp, cx); //cx = [x, cx];
|
||
tmp << yo;
|
||
cy = join_horiz(tmp, cy); //cy = [yo, cy];
|
||
tmp << dyo;
|
||
cdy = join_horiz(tmp, cdy); //cdy = [dyo, cdy];
|
||
// continue from here
|
||
x = xn;
|
||
y = yn;
|
||
dy = 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
|
||
lFitRange(x, y, xb, xb, data, roiMask, 2, yn, dyn);
|
||
rowvec tmp;
|
||
tmp << xb; cx = join_horiz(tmp, cx); //cx = [xb, cx];
|
||
tmp << yn; cy = join_horiz(tmp, cy); //cy = [yn, cy];
|
||
tmp << dyn; cdy = join_horiz(tmp, cdy); //cdy = [dyn, cdy];
|
||
}
|
||
else
|
||
{
|
||
double xm = round((x+xb)/2);
|
||
rowvec tmp; tmp << xb << xm;
|
||
lFitRange(x, y, xb, tmp, data, roiMask, 4, yn, dyn);
|
||
|
||
cx = join_horiz(tmp, cx); //cx = [xb, xm, cx];
|
||
tmp << yn;
|
||
cy = join_horiz(tmp, 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
|
||
rowvec tmp;
|
||
tmp << xn; cx = join_horiz(tmp, cx); //cx = [xn, cx];
|
||
tmp << yn; cy = join_horiz(tmp, cy); //cy = [yn, cy];
|
||
tmp << dyn; cdy = join_horiz(tmp, cdy); //cdy = [dyn, cdy];
|
||
}
|
||
else {
|
||
double xc = ceil((x+xb)/2);
|
||
rowvec tmp; 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];
|
||
tmp << yn;
|
||
cy = join_horiz(tmp, 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
|
||
rowvec tmp;
|
||
tmp << xn; cx = join_horiz(tmp, cx); //cx = [xn, cx];
|
||
tmp << yn; cy = join_horiz(tmp, cy); //cy = [yn, cy];
|
||
tmp << gNaN; cdy = join_horiz(tmp, cdy); //cdy = [NaN, cdy];
|
||
// continue from here
|
||
x = xn;
|
||
y = yn;
|
||
dy = 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";
|
||
}
|
||
cy = cy + IdxMat(step, cx-1);
|
||
}
|
||
|
||
void tmpbase(rowvec& cx,rowvec& cy,rowvec& cdy)
|
||
{
|
||
|
||
}
|
||
|
||
void baseInit(rowvec& cx, rowvec& cy, rowvec& cdy,
|
||
double PSSBreak = 100, double PSSfwhm = 100)
|
||
{
|
||
//% 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;
|
||
|
||
lGetData(data, step, rg, roiMask, L, H, PSSBreak, PSSfwhm);
|
||
|
||
lInitBase(data, step, rg, roiMask, L, H, NMin, NCritical, NTwo, cx, cy, cdy);
|
||
tmpbase(cx,cy,cdy);
|
||
return;
|
||
}
|
||
|
||
void dmps()
|
||
{
|
||
% some default values
|
||
myNaN=NaN;
|
||
myDummy={'-'};
|
||
|
||
% get PAT tables depending on source
|
||
[s, msg, patEntry, isValidName] = getPATByName(sn, varargin{1});
|
||
patName = patEntry{1};
|
||
pstab = patEntry{2};
|
||
flags = patEntry{3};
|
||
nuclist = patEntry{4};
|
||
calNames = patEntry{5};
|
||
nucIdxStruc = patEntry{6};
|
||
|
||
if isValidName
|
||
varargin(1) = [];
|
||
end
|
||
|
||
% number of PAT entries
|
||
NPeaks=size(pstab,1);
|
||
|
||
%flags which vars are required
|
||
Centroid_n = 1;
|
||
FWHM_Ch_n = 2;
|
||
NetArea_n = 3;
|
||
AreaError_n = 4;
|
||
Step_n = 5;
|
||
Multiplet_n = 6;
|
||
MeanBackCount_n = 7;
|
||
Sensitivity_n = 8;
|
||
Source_n = 9;
|
||
FitFlags_n = 10;
|
||
Spurious_n = 11;
|
||
Reviewed_n = 12;
|
||
Nuclide_n = 13;
|
||
Energy_n = 14;
|
||
FWHM_n = 15;
|
||
Efficiency_n = 16;
|
||
AreaErrorP_n = 17;
|
||
MultipletStr_n = 18;
|
||
LineDeviation_n = 19;
|
||
dE_n = 20;
|
||
Res_n = 21;
|
||
Res_Ch_n = 22;
|
||
SignificancePAT_n = 23;
|
||
SignificanceFlag_n = 24;
|
||
Index_n = 25;
|
||
ManualTail_n = 26;
|
||
Tail_n = 27;
|
||
TailFlag_n = 28;
|
||
TailChanged_n = 29;
|
||
FWHMIsFitted_n = 30;
|
||
CentroidErr_n = 31;
|
||
EnergyPAT_n = 32;
|
||
FWHM_ChPAT_n = 33;
|
||
FWHMPAT_n = 34;
|
||
FWHMErr_n = 35;
|
||
BWWidthPAT_n = 36;
|
||
BWWidth_n = 37;
|
||
StepRatioPAT_n = 38;
|
||
StepErr_n = 39;
|
||
TailErr_n = 40;
|
||
AreaErrorPAT_n = 41;
|
||
EmissionRatePAT_n = 42;
|
||
EmissionRate_n = 43;
|
||
EmissionRateErr_n = 44;
|
||
CCFPAT_n = 45;
|
||
CCF_n = 46;
|
||
CCFErrPAT_n = 47;
|
||
SignificanceErr_n = 48;
|
||
MultipletPAT_n = 49;
|
||
EfficiencyPAT_n = 50;
|
||
EfficiencyErrPAT_n = 51;
|
||
EnergyErr_n = 52;
|
||
BWWidthErr_n = 53;
|
||
ResChanged_n = 54;
|
||
Step_n = 55;
|
||
StepRatio_n = 56;
|
||
StepRatioErr_n = 57;
|
||
TailAlpha_n = 58;
|
||
UpperTail_n = 59;
|
||
UpperTailAlpha_n = 60;
|
||
BWWidthChan_n = 61;
|
||
NetAreaFree_n = 62;
|
||
NetAreaFlag_n = 63;
|
||
LineEnergy_n = 64;
|
||
FWHMFitted_n = 65;
|
||
CountsPerSecond_n = 66;
|
||
appERate_n = 67;
|
||
Yield_n = 68;
|
||
YieldErr_n = 69;
|
||
ReferenceID_n = 70;
|
||
RefPointer_n = 71;
|
||
BranchingRatio_n = 72;
|
||
NuclideHard_n = 73;
|
||
NuclideSoft_n = 74;
|
||
NetAreaCalculated_n = 75;
|
||
LineEnergyErr_n = 76;
|
||
ReferencePN_n = 77;
|
||
LC_n = 78;
|
||
LD_n = 79;
|
||
SignedResidual_n = 80;
|
||
UnsignedResidual_n = 81;
|
||
CCFErr_n = 82;
|
||
SumPeakArea_n = 83;
|
||
SumPeakAreaErrP_n = 84;
|
||
PeakShare_n = 85;
|
||
NuclideInternal_n = 86;
|
||
NACalcMin_n = 87;
|
||
NetAreaMin_n = 88;
|
||
PeakIdentified_n = 89;
|
||
HalfLife_n = 90;
|
||
HalfLifeStr_n = 91;
|
||
Activity_n = 92;
|
||
DecayCorrection_n = 93;
|
||
EfficiencyErr_n = 94;
|
||
BackgroundArea_n = 95;
|
||
Significance_n = 96;
|
||
NetAreaCalcError_n = 97;
|
||
NetAreaCalcUTest_n = 98;
|
||
RefEfficiencyErrPAT_n = 99;
|
||
RefEfficiencyErr_n = 100;
|
||
UpperTailBeta_n = 101;
|
||
RecoilBeta_n = 102;
|
||
RecoilBetaChan_n = 103;
|
||
RecoilDeltaE_n = 104;
|
||
RecoilDeltaChan_n = 105;
|
||
TailAlphaPAT_n = 106;
|
||
UpperTailPAT_n = 107;
|
||
UpperTailAlphaPAT_n = 108;
|
||
EffectiveHalflife_n = 109;
|
||
EffectiveHalflifeStr_n = 110;
|
||
HalfLifeError_n = 111;
|
||
EffectiveHalflifeError_n = 112;
|
||
SigmaCh_n = 113;
|
||
Sigma_n = 114;
|
||
VAR_NUMBER = 114;
|
||
|
||
depend = cell(VAR_NUMBER, 1);
|
||
% dependencies
|
||
depend{Centroid_n} = Centroid_n;
|
||
depend{CentroidErr_n} = CentroidErr_n;
|
||
depend{EnergyPAT_n} = EnergyPAT_n;
|
||
depend{FWHM_ChPAT_n} = FWHM_ChPAT_n;
|
||
depend{FWHM_Ch_n} = [FWHM_Ch_n, FWHM_ChPAT_n, Res_Ch_n, dE_n];
|
||
depend{NetArea_n} = NetArea_n;
|
||
depend{AreaErrorPAT_n} = AreaErrorPAT_n;
|
||
depend{AreaError_n} = [AreaErrorPAT_n, AreaError_n, NetArea_n, LC_n, ...
|
||
BackgroundArea_n];
|
||
depend{StepRatioPAT_n} = StepRatioPAT_n;
|
||
depend{Step_n} = [Step_n, StepRatio_n, NetArea_n];
|
||
depend{StepRatio_n} = [StepRatio_n, StepRatioPAT_n, Energy_n];
|
||
depend{StepRatioErr_n} = StepRatioErr_n;
|
||
depend{StepErr_n} = [StepErr_n, StepRatioErr_n, NetArea_n];
|
||
depend{TailErr_n} = TailErr_n;
|
||
depend{EmissionRatePAT_n} = EmissionRatePAT_n;
|
||
depend{EmissionRateErr_n} = EmissionRateErr_n;
|
||
depend{CCFPAT_n} = CCFPAT_n;
|
||
depend{CCFErrPAT_n} = CCFErrPAT_n;
|
||
depend{Yield_n} = Yield_n;
|
||
depend{YieldErr_n} = YieldErr_n;
|
||
depend{ReferenceID_n} = ReferenceID_n;
|
||
depend{EfficiencyErrPAT_n} = EfficiencyErrPAT_n;
|
||
depend{EfficiencyErr_n} = [EfficiencyErr_n, EfficiencyErrPAT_n, Energy_n];
|
||
depend{RefEfficiencyErr_n} = [RefEfficiencyErr_n, RefEfficiencyErrPAT_n, ...
|
||
EfficiencyErr_n, ReferencePN_n];
|
||
depend{EnergyErr_n} = EnergyErr_n;
|
||
depend{BWWidthErr_n} = BWWidthErr_n;
|
||
depend{MultipletPAT_n} = Multiplet_n;
|
||
depend{Multiplet_n} = [MultipletPAT_n, Multiplet_n, Centroid_n, FWHM_Ch_n, NetArea_n];
|
||
depend{MeanBackCount_n} = MeanBackCount_n;
|
||
depend{Sensitivity_n} = Sensitivity_n;
|
||
depend{SignificancePAT_n} = SignificancePAT_n;
|
||
depend{SignificanceErr_n}= SignificanceErr_n;
|
||
depend{LineEnergy_n} = LineEnergy_n;
|
||
depend{LineEnergyErr_n} = LineEnergyErr_n;
|
||
depend{LC_n} = LC_n;
|
||
depend{LD_n} = LD_n;
|
||
depend{SignedResidual_n} = SignedResidual_n;
|
||
depend{UnsignedResidual_n} = UnsignedResidual_n;
|
||
depend{SumPeakArea_n} = SumPeakArea_n;
|
||
depend{SumPeakAreaErrP_n} = SumPeakAreaErrP_n;
|
||
depend{Source_n} = Source_n;
|
||
depend{FitFlags_n} = FitFlags_n;
|
||
depend{Spurious_n} = Spurious_n;
|
||
depend{Reviewed_n} = Reviewed_n;
|
||
depend{NuclideHard_n} = NuclideHard_n;
|
||
depend{NuclideSoft_n} = NuclideSoft_n;
|
||
depend{Nuclide_n} = Nuclide_n;
|
||
depend{ManualTail_n} = ManualTail_n;
|
||
depend{TailChanged_n} = TailChanged_n;
|
||
depend{ResChanged_n} = ResChanged_n;
|
||
depend{BackgroundArea_n} = BackgroundArea_n;
|
||
depend{RefEfficiencyErrPAT_n} = RefEfficiencyErrPAT_n;
|
||
|
||
depend{Energy_n} = [EnergyPAT_n, Energy_n, Centroid_n];
|
||
depend{FWHM_n} = [FWHM_n, FWHMPAT_n, FWHM_Ch_n, dE_n];
|
||
depend{SigmaCh_n} = [SigmaCh_n, FWHM_Ch_n];
|
||
depend{Sigma_n} = [Sigma_n, FWHM_n];
|
||
depend{EfficiencyPAT_n} = EfficiencyPAT_n;
|
||
depend{Efficiency_n} = [Efficiency_n, EfficiencyPAT_n, Energy_n];
|
||
depend{AreaErrorP_n} = [AreaErrorP_n, NetArea_n, AreaError_n, LC_n];
|
||
depend{MultipletStr_n} = [MultipletStr_n, Multiplet_n];
|
||
depend{LineDeviation_n} = [LineDeviation_n, Energy_n, LineEnergy_n];
|
||
depend{dE_n} = [dE_n, Centroid_n];
|
||
depend{Res_n} = [Res_n, Energy_n];
|
||
depend{Res_Ch_n} = [Res_Ch_n, Res_n, dE_n];
|
||
depend{Significance_n} = [Significance_n, SignificancePAT_n, NetArea_n, ...
|
||
LC_n, Source_n];
|
||
depend{SignificanceFlag_n} = [SignificanceFlag_n, Significance_n];
|
||
depend{Index_n} = Index_n;
|
||
depend{Tail_n} = [Tail_n, ManualTail_n, Energy_n];
|
||
depend{TailFlag_n} = [TailFlag_n, ManualTail_n];
|
||
depend{FWHMIsFitted_n} = [FWHMIsFitted_n, FitFlags_n];
|
||
depend{FWHMFitted_n} = [FWHMFitted_n, FWHMIsFitted_n, FWHM_n];
|
||
depend{FWHMPAT_n} = FWHMPAT_n;
|
||
depend{FWHMErr_n} = FWHMErr_n;
|
||
depend{BWWidthPAT_n} = BWWidthPAT_n;
|
||
depend{BWWidth_n} = [BWWidth_n, BWWidthPAT_n];
|
||
depend{BWWidthChan_n} = [BWWidthChan_n, BWWidth_n, dE_n];
|
||
depend{EmissionRate_n} = [EmissionRate_n, EmissionRatePAT_n, NetArea_n, ...
|
||
Efficiency_n];
|
||
depend{CCF_n} = [CCF_n, CCFPAT_n];
|
||
depend{CCFErr_n} = [CCFErr_n, CCFErrPAT_n, CCF_n];
|
||
depend{TailAlphaPAT_n} = TailAlphaPAT_n;
|
||
depend{TailAlpha_n} = [TailAlpha_n, TailAlphaPAT_n, Energy_n];
|
||
depend{UpperTailPAT_n} = [UpperTailPAT_n];
|
||
depend{UpperTail_n} = [UpperTail_n, UpperTailPAT_n, Energy_n];
|
||
depend{UpperTailAlphaPAT_n} = [UpperTailAlphaPAT_n];
|
||
depend{UpperTailAlpha_n} = [UpperTailAlpha_n, UpperTailAlphaPAT_n, Energy_n];
|
||
depend{NetAreaFlag_n} = [NetAreaFlag_n];
|
||
depend{NetAreaFree_n} = [NetAreaFree_n, Source_n, NetAreaFlag_n];
|
||
depend{CountsPerSecond_n} = [CountsPerSecond_n, NetArea_n];
|
||
depend{appERate_n} = [appERate_n, CountsPerSecond_n, Efficiency_n];
|
||
depend{Nuclide_n} = [Nuclide_n, NuclideHard_n, NuclideSoft_n];
|
||
depend{ReferencePN_n} = [ReferencePN_n, ReferenceID_n, RefPointer_n];
|
||
depend{NuclideInternal_n} = [NuclideInternal_n, NuclideHard_n];
|
||
depend{PeakIdentified_n} = [PeakIdentified_n, NuclideHard_n];
|
||
depend{RefPointer_n} = [RefPointer_n, NuclideInternal_n];
|
||
depend{BranchingRatio_n} = [BranchingRatio_n, NuclideInternal_n, ReferencePN_n];
|
||
depend{NetAreaCalculated_n} = [NetAreaCalculated_n, Yield_n, ReferencePN_n, ...
|
||
BranchingRatio_n, CCF_n, NetArea_n, Efficiency_n];
|
||
depend{NetAreaCalcError_n} = [NetAreaCalcError_n, ReferencePN_n, ...
|
||
YieldErr_n, AreaErrorP_n, CCFErr_n, RefEfficiencyErr_n];
|
||
depend{NetAreaCalcUTest_n} = [NetAreaCalcUTest_n, NetArea_n, AreaError_n, ...
|
||
NetAreaCalculated_n, NACalcMin_n, NetAreaCalcError_n];
|
||
depend{NACalcMin_n} = [NACalcMin_n, Yield_n, ReferencePN_n, ...
|
||
BranchingRatio_n, CCF_n, NetAreaMin_n, Efficiency_n];
|
||
depend{NetAreaMin_n} = [NetAreaMin_n, NetArea_n, LC_n];
|
||
depend{PeakShare_n} = [PeakShare_n, SumPeakArea_n, ...
|
||
NetAreaCalculated_n, NetArea_n];
|
||
depend{HalfLife_n} = [HalfLife_n, NuclideInternal_n];
|
||
depend{HalfLifeError_n} = [HalfLifeError_n, NuclideInternal_n];
|
||
depend{HalfLifeStr_n} = [HalfLifeStr_n, NuclideInternal_n];
|
||
depend{EffectiveHalflife_n} = [EffectiveHalflife_n, NuclideInternal_n];
|
||
depend{EffectiveHalflifeError_n} = [EffectiveHalflifeError_n,NuclideInternal_n];
|
||
depend{EffectiveHalflifeStr_n} = [EffectiveHalflife_n, EffectiveHalflifeStr_n];
|
||
depend{DecayCorrection_n} = [DecayCorrection_n, HalfLife_n];
|
||
depend{Activity_n} = [Activity_n, NetArea_n, CCF_n, Efficiency_n, Yield_n, ...
|
||
DecayCorrection_n];
|
||
depend{UpperTailBeta_n} = [UpperTailBeta_n, UpperTail_n, UpperTailAlpha_n, ...
|
||
SigmaCh_n, dE_n];
|
||
depend{RecoilBeta_n} = RecoilBeta_n;
|
||
depend{RecoilBetaChan_n} = [RecoilBetaChan_n, RecoilBeta_n, dE_n];
|
||
depend{RecoilDeltaE_n} = RecoilDeltaE_n;
|
||
depend{RecoilDeltaChan_n} = [RecoilDeltaChan_n, RecoilDeltaE_n, dE_n];
|
||
|
||
% which variables are needed
|
||
needed = zeros(VAR_NUMBER, 1);
|
||
|
||
% parse descriptors and remember needed variables
|
||
for i = 1 : length(varargin)
|
||
desc = varargin{i};
|
||
switch desc
|
||
case 'Activity'
|
||
needed(depend{Activity_n}) = 1;
|
||
case 'appERate'
|
||
needed(depend{appERate_n}) = 1;
|
||
case 'AreaError'
|
||
needed(depend{AreaError_n}) = 1;
|
||
case 'AreaErrorP'
|
||
needed(depend{AreaErrorP_n}) = 1;
|
||
case 'AreaErrorPAT'
|
||
needed(depend{AreaErrorPAT_n}) = 1;
|
||
case 'BackgroundArea'
|
||
needed(depend{BackgroundArea_n}) = 1;
|
||
case 'BranchingRatio'
|
||
needed(depend{BranchingRatio_n}) = 1;
|
||
case 'BWWidth'
|
||
needed(depend{BWWidth_n}) = 1;
|
||
case 'BWWidthChan'
|
||
needed(depend{BWWidthChan_n}) = 1;
|
||
case 'BWWidthErr'
|
||
needed(depend{BWWidthErr_n}) = 1;
|
||
case 'BWWidthPAT'
|
||
needed(depend{BWWidthPAT_n}) = 1;
|
||
case 'Centroid'
|
||
needed(depend{Centroid_n}) = 1;
|
||
case 'CentroidErr'
|
||
needed(depend{CentroidErr_n}) = 1;
|
||
case 'CCF'
|
||
needed(depend{CCF_n}) = 1;
|
||
case 'CCFPAT'
|
||
needed(depend{CCFPAT_n}) = 1;
|
||
case 'CCFErr'
|
||
needed(depend{CCFErr_n}) = 1;
|
||
case 'CountsPerSecond'
|
||
needed(depend{CountsPerSecond_n}) = 1;
|
||
case 'dE'
|
||
needed(depend{dE_n}) = 1;
|
||
case 'EffectiveHalflife'
|
||
needed(depend{EffectiveHalflife_n}) = 1;
|
||
case 'EffectiveHalflifeError'
|
||
needed(depend{EffectiveHalflifeError_n}) = 1;
|
||
case 'EffectiveHalflifeStr'
|
||
needed(depend{EffectiveHalflifeStr_n}) = 1;
|
||
case 'Efficiency'
|
||
needed(depend{Efficiency_n}) = 1;
|
||
case 'EfficiencyErr'
|
||
needed(depend{EfficiencyErr_n}) = 1;
|
||
case 'EfficiencyPAT'
|
||
needed(depend{EfficiencyPAT_n}) = 1;
|
||
case 'EmissionRate'
|
||
needed(depend{EmissionRate_n}) = 1;
|
||
case 'EmissionRateErr'
|
||
needed(depend{EmissionRateErr_n}) = 1;
|
||
case 'EmissionRatePAT'
|
||
needed(depend{EmissionRatePAT_n}) = 1;
|
||
case 'Energy'
|
||
needed(depend{Energy_n}) = 1;
|
||
case 'EnergyErr'
|
||
needed(depend{EnergyErr_n}) = 1;
|
||
case 'EnergyPAT'
|
||
needed(depend{EnergyPAT_n}) = 1;
|
||
case 'FitFlags'
|
||
needed(depend{FitFlags_n}) = 1;
|
||
case 'FWHM'
|
||
needed(depend{FWHM_n}) = 1;
|
||
case 'FWHM_Ch'
|
||
needed(depend{FWHM_Ch_n}) = 1;
|
||
case 'FWHM_ChPAT'
|
||
needed(depend{FWHM_ChPAT_n}) = 1;
|
||
case 'FWHMErr'
|
||
needed(depend{FWHMErr_n}) = 1;
|
||
case 'FWHMFitted'
|
||
needed(depend{FWHMFitted_n}) = 1;
|
||
case 'FWHMIsFitted'
|
||
needed(depend{FWHMIsFitted_n}) = 1;
|
||
case 'FWHMPAT'
|
||
needed(depend{FWHMPAT_n}) = 1;
|
||
case 'HalfLife'
|
||
needed(depend{HalfLife_n}) = 1;
|
||
case 'HalfLifeError'
|
||
needed(depend{HalfLifeError_n}) = 1;
|
||
case 'HalfLifeStr'
|
||
needed(depend{HalfLifeStr_n}) = 1;
|
||
case 'Index'
|
||
needed(depend{Index_n}) = 1;
|
||
case 'LC'
|
||
needed(depend{LC_n}) = 1;
|
||
case 'LD'
|
||
needed(depend{LD_n}) = 1;
|
||
case 'LineDeviation'
|
||
needed(depend{LineDeviation_n}) = 1;
|
||
case 'LineEnergy'
|
||
needed(depend{LineEnergy_n}) = 1;
|
||
case 'LineEnergyErr'
|
||
needed(depend{LineEnergyErr_n}) = 1;
|
||
case 'ManualTail'
|
||
needed(depend{ManualTail_n}) = 1;
|
||
case 'MeanBackCount'
|
||
needed(depend{MeanBackCount_n}) = 1;
|
||
case 'Multiplet'
|
||
needed(depend{Multiplet_n}) = 1;
|
||
case 'MultipletPAT'
|
||
needed(depend{MultipletPAT_n}) = 1;
|
||
case 'MultipletStr'
|
||
needed(depend{MultipletStr_n}) = 1;
|
||
case 'NACalcMin'
|
||
needed(depend{NACalcMin_n}) = 1;
|
||
case 'NetArea'
|
||
needed(depend{NetArea_n}) = 1;
|
||
case 'NetAreaCalculated'
|
||
needed(depend{NetAreaCalculated_n}) = 1;
|
||
case 'NetAreaCalcError'
|
||
needed(depend{NetAreaCalcError_n}) = 1;
|
||
case 'NetAreaCalcUTest'
|
||
needed(depend{NetAreaCalcUTest_n}) = 1;
|
||
case 'NetAreaFlag'
|
||
needed(depend{NetAreaFlag_n}) = 1;
|
||
case 'NetAreaFree'
|
||
needed(depend{NetAreaFree_n}) = 1;
|
||
case 'Nuclide'
|
||
needed(depend{Nuclide_n}) = 1;
|
||
case 'NuclideHard'
|
||
needed(depend{NuclideHard_n}) = 1;
|
||
case 'NuclideInternal'
|
||
needed(depend{NuclideInternal_n}) = 1;
|
||
case 'NuclideSoft'
|
||
needed(depend{NuclideSoft_n}) = 1;
|
||
case 'PeakIdentified'
|
||
needed(depend{PeakIdentified_n}) = 1;
|
||
case 'PeakShare'
|
||
needed(depend{PeakShare_n}) = 1;
|
||
case 'RecoilBeta'
|
||
needed(depend{RecoilBeta_n}) = 1;
|
||
case 'RecoilBetaChan'
|
||
needed(depend{RecoilBetaChan_n}) = 1;
|
||
case 'RecoilDeltaE'
|
||
needed(depend{RecoilDeltaE_n}) = 1;
|
||
case 'RecoilDeltaChan'
|
||
needed(depend{RecoilDeltaChan_n}) = 1;
|
||
case 'RefEfficiencyErr'
|
||
needed(depend{RefEfficiencyErr_n}) = 1;
|
||
case 'ReferenceID'
|
||
needed(depend{ReferenceID_n}) = 1;
|
||
case 'ReferencePN'
|
||
needed(depend{ReferencePN_n}) = 1;
|
||
case 'RefPointer'
|
||
needed(depend{RefPointer_n}) = 1;
|
||
case 'Res'
|
||
needed(depend{Res_n}) = 1;
|
||
case 'Res_Ch'
|
||
needed(depend{Res_Ch_n}) = 1;
|
||
case 'ResChanged'
|
||
needed(depend{ResChanged_n}) = 1;
|
||
case 'Reviewed'
|
||
needed(depend{Reviewed_n}) = 1;
|
||
case 'Sensitivity'
|
||
needed(depend{Sensitivity_n}) = 1;
|
||
case 'Sigma'
|
||
needed(depend{Sigma_n}) = 1;
|
||
case 'SigmaCh'
|
||
needed(depend{SigmaCh_n}) = 1;
|
||
case 'SignedResidual'
|
||
needed(depend{SignedResidual_n}) = 1;
|
||
case 'Significance'
|
||
needed(depend{Significance_n}) = 1;
|
||
case 'SignificanceErr'
|
||
needed(depend{SignificanceErr_n}) = 1;
|
||
case 'SignificanceFlag'
|
||
needed(depend{SignificanceFlag_n}) = 1;
|
||
case 'Source'
|
||
needed(depend{Source_n}) = 1;
|
||
case 'Spurious'
|
||
needed(depend{Spurious_n}) = 1;
|
||
case 'Step'
|
||
needed(depend{Step_n}) = 1;
|
||
case 'StepErr'
|
||
needed(depend{StepErr_n}) = 1;
|
||
case 'StepRatio'
|
||
needed(depend{StepRatio_n}) = 1;
|
||
case 'StepRatioErr'
|
||
needed(depend{StepRatioErr_n}) = 1;
|
||
case 'StepRatioPAT'
|
||
needed(depend{StepRatioPAT_n}) = 1;
|
||
case 'SumPeakArea'
|
||
needed(depend{SumPeakArea_n}) = 1;
|
||
case 'SumPeakAreaErrP'
|
||
needed(depend{SumPeakAreaErrP_n}) = 1;
|
||
case 'Tail'
|
||
needed(depend{Tail_n}) = 1;
|
||
case 'TailAlpha'
|
||
needed(depend{TailAlpha_n}) = 1;
|
||
case 'TailAlphaPAT'
|
||
needed(depend{TailAlphaPAT_n}) = 1;
|
||
case 'TailChanged'
|
||
needed(depend{TailChanged_n}) = 1;
|
||
case 'TailErr'
|
||
needed(depend{TailErr_n}) = 1;
|
||
case 'TailFlag'
|
||
needed(depend{TailFlag_n}) = 1;
|
||
case 'UnsignedResidual'
|
||
needed(depend{UnsignedResidual_n}) = 1;
|
||
case 'UpperTail'
|
||
needed(depend{UpperTail_n}) = 1;
|
||
case 'UpperTailPAT'
|
||
needed(depend{UpperTailPAT_n}) = 1;
|
||
case 'UpperTailAlpha'
|
||
needed(depend{UpperTailAlpha_n}) = 1;
|
||
case 'UpperTailAlphaPAT'
|
||
needed(depend{UpperTailAlphaPAT_n}) = 1;
|
||
case 'UpperTailBeta'
|
||
needed(depend{UpperTailBeta_n}) = 1;
|
||
case 'Yield'
|
||
needed(depend{Yield_n}) = 1;
|
||
case 'YieldErr'
|
||
needed(depend{YieldErr_n}) = 1;
|
||
otherwise
|
||
disp(desc)
|
||
error('invalid descriptor')
|
||
end
|
||
end
|
||
|
||
% resolve dependencies
|
||
oldNeeded = zeros(VAR_NUMBER, 1);
|
||
% while any new requests added
|
||
while any(needed - oldNeeded)
|
||
oldNeeded = needed;
|
||
for i = 1 : VAR_NUMBER
|
||
% if we need i'th variable
|
||
if needed(i)
|
||
% also request all variables, on which i depends
|
||
needed(depend{i}) = 1;
|
||
end
|
||
end
|
||
end
|
||
|
||
% extract and calculate the actual values
|
||
|
||
if needed(Centroid_n)
|
||
Centroid = pstab(:,1);
|
||
end
|
||
|
||
if needed(CentroidErr_n)
|
||
CentroidErr = pstab(:,2);
|
||
end
|
||
|
||
if needed(EnergyPAT_n)
|
||
EnergyPAT = pstab(:,3);
|
||
end
|
||
|
||
if needed(EnergyErr_n)
|
||
EnergyErr = pstab(:,4);
|
||
end
|
||
|
||
if needed(FWHM_ChPAT_n)
|
||
FWHM_ChPAT = pstab(:,5);
|
||
end
|
||
|
||
if needed(FWHMPAT_n)
|
||
FWHMPAT = pstab(:,6);
|
||
end
|
||
|
||
if needed(FWHMErr_n)
|
||
FWHMErr = pstab(:,7);
|
||
end
|
||
|
||
if needed(StepRatioPAT_n)
|
||
StepRatioPAT = pstab(:,8);
|
||
end
|
||
|
||
if needed(StepRatioErr_n)
|
||
StepRatioErr = pstab(:,9);
|
||
end
|
||
|
||
if needed(ManualTail_n)
|
||
ManualTail = pstab(:,10);
|
||
end
|
||
|
||
if needed(TailErr_n)
|
||
TailErr = pstab(:,11);
|
||
end
|
||
|
||
if needed(NetArea_n)
|
||
NetArea = pstab(:,12);
|
||
end
|
||
|
||
if needed(AreaErrorPAT_n)
|
||
AreaErrorPAT = pstab(:,13);
|
||
end
|
||
|
||
if needed(MeanBackCount_n)
|
||
MeanBackCount = pstab(:,14);
|
||
end
|
||
|
||
if needed(Sensitivity_n)
|
||
Sensitivity = pstab(:,15);
|
||
end
|
||
|
||
if needed(SignificancePAT_n)
|
||
SignificancePAT = pstab(:,16);
|
||
end
|
||
|
||
if needed(SignificanceErr_n)
|
||
SignificanceErr = pstab(:,17);
|
||
end
|
||
|
||
if needed(MultipletPAT_n)
|
||
MultipletPAT = pstab(:,18);
|
||
end
|
||
|
||
if needed(EfficiencyPAT_n)
|
||
EfficiencyPAT = pstab(:,19);
|
||
end
|
||
|
||
if needed(EfficiencyErrPAT_n)
|
||
EfficiencyErrPAT = pstab(:,20);
|
||
end
|
||
|
||
if needed(BWWidthPAT_n)
|
||
BWWidthPAT = pstab(:,21);
|
||
end
|
||
|
||
if needed(BWWidthErr_n)
|
||
BWWidthErr = pstab(:,22);
|
||
end
|
||
|
||
if needed(EmissionRatePAT_n)
|
||
EmissionRatePAT = pstab(:,23);
|
||
end
|
||
|
||
if needed(EmissionRateErr_n)
|
||
EmissionRateErr = pstab(:,24);
|
||
end
|
||
|
||
if needed(CCFPAT_n)
|
||
CCFPAT = pstab(:,25);
|
||
end
|
||
|
||
if needed(CCFErrPAT_n)
|
||
CCFErrPAT = pstab(:,26);
|
||
end
|
||
|
||
if needed(LineEnergy_n)
|
||
LineEnergy = pstab(:,27);
|
||
end
|
||
|
||
if needed(LineEnergyErr_n)
|
||
LineEnergyErr = pstab(:, 28);
|
||
end
|
||
|
||
if needed(Yield_n)
|
||
Yield = pstab(:, 29);
|
||
end
|
||
|
||
if needed(YieldErr_n)
|
||
YieldErr = pstab(:, 30);
|
||
end
|
||
|
||
if needed(ReferenceID_n)
|
||
ReferenceID = pstab(:, 31);
|
||
end
|
||
|
||
if needed(LC_n)
|
||
LC = pstab(:, 34);
|
||
end
|
||
|
||
if needed(LD_n)
|
||
LD = pstab(:, 35);
|
||
end
|
||
|
||
if any(needed([SignedResidual_n, UnsignedResidual_n]))
|
||
|
||
% updated only when requested
|
||
ResidualUpToDate = pstab(:, 38);
|
||
if any(ResidualUpToDate ~= 1)
|
||
[SignedResidual, UnsignedResidual] = patBaseVar(sn, 'all', patName, 1);
|
||
else
|
||
|
||
% variables are up to date
|
||
if needed(SignedResidual_n)
|
||
SignedResidual = pstab(:, 36);
|
||
end
|
||
|
||
if needed(UnsignedResidual_n)
|
||
UnsignedResidual = pstab(:, 37);
|
||
end
|
||
end
|
||
end
|
||
|
||
if needed(SumPeakArea_n)
|
||
SumPeakArea = pstab(:, 39);
|
||
end
|
||
|
||
if needed(SumPeakAreaErrP_n)
|
||
SumPeakAreaErrP = pstab(:, 40);
|
||
end
|
||
|
||
if needed(BackgroundArea_n)
|
||
BackgroundArea = pstab(:, 41);
|
||
end
|
||
|
||
if needed(RefEfficiencyErrPAT_n)
|
||
RefEfficiencyErrPAT = pstab(:, 42);
|
||
end
|
||
|
||
if needed(RecoilBeta_n)
|
||
RecoilBeta = pstab(:, 46);
|
||
end
|
||
|
||
if needed(RecoilDeltaE_n)
|
||
RecoilDeltaE = pstab(:, 47);
|
||
end
|
||
|
||
if needed(Source_n)
|
||
Source = flags(:, 1);
|
||
end
|
||
|
||
if needed(FitFlags_n)
|
||
% old: FitFlags = flags(:, 6 : 17);
|
||
FitFlags = flags(:, [8, 9, 12, 13, 16, 17]);
|
||
end
|
||
if needed(NetAreaFlag_n)
|
||
NetAreaFlag = flags(:, 12);
|
||
end
|
||
|
||
if needed(Spurious_n)
|
||
Spurious = flags(:, 26);
|
||
end
|
||
|
||
if needed(Reviewed_n)
|
||
Reviewed = flags(:, 27);
|
||
end
|
||
|
||
if needed(TailChanged_n)
|
||
TailChanged = flags(:, 28);
|
||
end
|
||
|
||
if needed(ResChanged_n)
|
||
ResChanged = flags(:, 29);
|
||
end
|
||
|
||
if needed(NuclideSoft_n)
|
||
NuclideSoft = nuclist(:, 1);
|
||
for i = 1 : NPeaks
|
||
if ~isstr(NuclideSoft{i})
|
||
NuclideSoft{i} = '';
|
||
end
|
||
end
|
||
end
|
||
|
||
if needed(NuclideHard_n)
|
||
NuclideHard = nuclist(:, 2);
|
||
for i = 1 : NPeaks
|
||
if ~isstr(NuclideHard{i})
|
||
NuclideHard{i} = '';
|
||
end
|
||
end
|
||
end
|
||
|
||
if any(needed([NuclideInternal_n, PeakIdentified_n]))
|
||
[PeakIdentified, NuclideInternal] = libInternalName(NuclideHard);
|
||
end
|
||
|
||
if needed(RefPointer_n)
|
||
[unuc, i, j] = unique(NuclideInternal);
|
||
refptr = repmat(NaN, size(unuc));
|
||
for k = 1 : length(unuc);
|
||
if isempty(unuc{k})
|
||
% skip
|
||
elseif isfield(nucIdxStruc, unuc{k})
|
||
refptr(k) = nucIdxStruc.(unuc{k}).ReferenceID;
|
||
else
|
||
% skip
|
||
end
|
||
end
|
||
RefPointer = refptr(j);
|
||
end
|
||
|
||
if needed(Nuclide_n)
|
||
Nuclide = NuclideHard;
|
||
for i = 1 : NPeaks
|
||
if isempty(NuclideHard{i})
|
||
Nuclide{i} = NuclideSoft{i};
|
||
end
|
||
end
|
||
end
|
||
|
||
if needed(Energy_n)
|
||
Energy = EnergyPAT;
|
||
bool = isnan(Energy);
|
||
if any(bool)
|
||
Energy(bool) = calValuesByName(sn, 'Energy', calNames{1}, Centroid(bool));
|
||
end
|
||
end
|
||
|
||
if needed(dE_n)
|
||
dE = calDerivativeByName(sn, 'Energy', calNames{1}, Centroid);
|
||
end
|
||
|
||
if needed(Res_n)
|
||
Res = calValuesByName(sn, 'Resolution', calNames{2}, Energy);
|
||
end
|
||
|
||
if needed(Res_Ch_n)
|
||
Res_Ch=Res./dE;
|
||
end
|
||
|
||
if needed(FWHM_Ch_n)
|
||
FWHM_Ch = FWHM_ChPAT;
|
||
bool = isnan(FWHM_Ch);
|
||
if any(bool)
|
||
FWHM_Ch(bool) = Res_Ch(bool);
|
||
end
|
||
end
|
||
|
||
if needed(FWHM_n)
|
||
FWHM = FWHMPAT;
|
||
bool = isnan(FWHM);
|
||
if any(bool)
|
||
FWHM(bool) = FWHM_Ch(bool).*dE(bool);
|
||
end
|
||
end
|
||
|
||
if needed(SigmaCh_n)
|
||
SigmaCh = FWHM_Ch/sqrt(8*log(2));
|
||
end
|
||
|
||
if needed(Sigma_n)
|
||
Sigma = FWHM/sqrt(8*log(2));
|
||
end
|
||
|
||
if needed(Efficiency_n)
|
||
Efficiency = EfficiencyPAT;
|
||
bool = isnan(Efficiency);
|
||
if any(bool)
|
||
Efficiency(bool) = calValuesByName(sn, 'Efficiency', calNames{3}, ...
|
||
Energy(bool));
|
||
end
|
||
end
|
||
|
||
if needed(EfficiencyErr_n)
|
||
EfficiencyErr = EfficiencyErrPAT;
|
||
bool = isnan(EfficiencyErr);
|
||
if any(bool)
|
||
[s, msg, p, perr] = calGetPATPara(sn, 'Efficiency', patName);
|
||
if s
|
||
% relative error in percent
|
||
EfficiencyErr(bool) = 100 * calErrorEval(Energy(bool), p, perr);
|
||
end
|
||
end
|
||
end
|
||
|
||
if needed(Tail_n)
|
||
Tail = ManualTail;
|
||
bool = isnan(Tail);
|
||
if any(bool)
|
||
Tail(bool) = calValuesByName(sn, 'Tail', calNames{5}, Energy(bool));
|
||
end
|
||
end
|
||
|
||
if needed(TailAlphaPAT_n)
|
||
TailAlphaPAT = pstab(:, 43);
|
||
end
|
||
|
||
if needed(TailAlpha_n)
|
||
TailAlpha = TailAlphaPAT;
|
||
bool = isnan(TailAlpha);
|
||
if any(bool)
|
||
TailAlpha(bool) = calValuesByName(sn, 'Tail_alpha', calNames{6}, ...
|
||
Energy(bool));
|
||
end
|
||
end
|
||
|
||
if needed(UpperTailPAT_n)
|
||
UpperTailPAT = pstab(:, 44);
|
||
end
|
||
|
||
if needed(UpperTail_n)
|
||
UpperTail = UpperTailPAT;
|
||
bool = isnan(UpperTail);
|
||
if any(bool)
|
||
UpperTail(bool) = calValuesByName(sn, 'Tail_right', calNames{7}, ...
|
||
Energy(bool));
|
||
end
|
||
end
|
||
|
||
if needed(UpperTailAlphaPAT_n)
|
||
UpperTailAlphaPAT = pstab(:, 45);
|
||
end
|
||
|
||
if needed(UpperTailAlpha_n)
|
||
UpperTailAlpha = UpperTailAlphaPAT;
|
||
bool = isnan(UpperTailAlpha);
|
||
if any(bool)
|
||
UpperTailAlpha(bool) = calValuesByName(sn, 'Tail_right_alpha', ...
|
||
calNames{8}, Energy(bool));
|
||
end
|
||
end
|
||
|
||
if needed(UpperTailBeta_n)
|
||
bool = (UpperTail > 0 & UpperTailAlpha == 1);
|
||
SigmaCh = FWHM_Ch/sqrt(8*log(2));
|
||
Beta = 1 ./ (UpperTail(bool) .* SigmaCh(bool));
|
||
UpperTailBeta = repmat(NaN, size(dE));
|
||
UpperTailBeta(bool) = Beta./dE(bool);
|
||
end
|
||
|
||
if needed(StepRatio_n)
|
||
StepRatio = StepRatioPAT;
|
||
bool = isnan(StepRatio);
|
||
if any(bool)
|
||
StepRatio(bool) = calValuesByName(sn, 'Step_ratio', calNames{9}, ...
|
||
Energy(bool));
|
||
end
|
||
end
|
||
|
||
if needed(Step_n)
|
||
Step = NetArea .* StepRatio;
|
||
end
|
||
|
||
if needed(StepErr_n)
|
||
StepErr = NetArea .* StepRatioErr;
|
||
end
|
||
|
||
if needed(Significance_n)
|
||
Significance = lCalcSignificance(SignificancePAT, NetArea, LC);
|
||
% Do not report significance for neutron lumps
|
||
Significance(Source == 'G') = NaN;
|
||
end
|
||
|
||
if needed(AreaError_n)
|
||
AreaError = lCalcAreaError(AreaErrorPAT, NetArea, BackgroundArea, LC);
|
||
end
|
||
|
||
if needed(AreaErrorP_n)
|
||
AreaErrorP = 100 * AreaError ./ max(NetArea, LC);
|
||
end
|
||
|
||
% use original multiplet flag, where they are finite
|
||
if needed(Multiplet_n)
|
||
Multiplet = MultipletPAT;
|
||
bool = isnan(Multiplet);
|
||
if any(bool)
|
||
% where the PAT entry is NaN, set calculated flag; for calculation,
|
||
% all peaks are needed !!
|
||
tmpMultiplet = peakFindMultiplets(Centroid, NetArea, FWHM_Ch, 2);
|
||
Multiplet(bool) = tmpMultiplet(bool);
|
||
end
|
||
end
|
||
|
||
if needed(MultipletStr_n)
|
||
MultipletStr=cell(size(Multiplet));
|
||
MultipletStr(find(Multiplet==0))={'.'};
|
||
MultipletStr(find(Multiplet==1))={'M'};
|
||
MultipletStr(find(Multiplet~=0&Multiplet~=1))={'U'};
|
||
end
|
||
|
||
if needed(LineDeviation_n)
|
||
LineDeviation = LineEnergy-Energy;
|
||
end
|
||
|
||
if needed(SignificanceFlag_n)
|
||
SignificanceFlag=sigflag(Significance);
|
||
end
|
||
|
||
if needed(Index_n)
|
||
Index = [1 : NPeaks]';
|
||
end
|
||
|
||
if needed(TailFlag_n)
|
||
% flag manual tails with 'T', others with '.'
|
||
dot = '.';
|
||
TailFlag = dot(ones(size(ManualTail)));
|
||
TailFlag(isfinite(ManualTail)) = 'T';
|
||
end
|
||
|
||
if needed(FWHMIsFitted_n)
|
||
FWHMFlag = FitFlags(:, 5);
|
||
FWHMIsFitted = (FWHMFlag == 'Q' | FWHMFlag == 'F');
|
||
end
|
||
|
||
if needed(FWHMFitted_n)
|
||
FWHMFitted = FWHM;
|
||
if ~all(FWHMIsFitted)
|
||
FWHMFitted(~FWHMIsFitted) = NaN;
|
||
end
|
||
end
|
||
|
||
if needed(BWWidth_n)
|
||
BWWidth = BWWidthPAT;
|
||
bool = isnan(BWWidth);
|
||
if any(bool)
|
||
BWWidth(bool) = 0;
|
||
end
|
||
end
|
||
|
||
if needed(BWWidthChan_n)
|
||
BWWidthChan = BWWidth ./ dE;
|
||
end
|
||
|
||
if needed(EmissionRate_n)
|
||
EmissionRate = EmissionRatePAT;
|
||
bool = isnan(EmissionRate);
|
||
tLive = dminfo(sn, 'AcquisitionLive');
|
||
if any(bool)
|
||
EmissionRate(bool) = NetArea ./ (tLive*Efficiency);
|
||
end
|
||
end
|
||
|
||
if needed(CCF_n)
|
||
CCF = CCFPAT;
|
||
bool = isnan(CCF);
|
||
if any(bool)
|
||
CCF(bool) = 1;
|
||
end
|
||
end
|
||
|
||
if needed(CCFErr_n)
|
||
CCFErr = CCFErrPAT;
|
||
bool = find(~isfinite(CCFErr));
|
||
if any(bool)
|
||
CCFErr(bool) = 10 * abs(CCF(bool) - 1);
|
||
end
|
||
end
|
||
|
||
if needed(CountsPerSecond_n)
|
||
lt = dminfo(sn, 'AcquisitionLive');
|
||
if lt > 0
|
||
CountsPerSecond = NetArea/lt;
|
||
else
|
||
CountsPerSecond = repmat(NaN, size(NetArea));
|
||
end
|
||
end
|
||
|
||
if needed(appERate_n)
|
||
appERate = CountsPerSecond ./ Efficiency;
|
||
end
|
||
|
||
if needed(NetAreaFree_n);
|
||
% net area not calculated or background
|
||
NetAreaFree = logical(NetAreaFlag ~= 'C' & Source ~= 'B' & Source ~= 'G');
|
||
end
|
||
|
||
if needed(ReferencePN_n)
|
||
ReferencePN = resolveDependency(ReferenceID, RefPointer);
|
||
end
|
||
|
||
if needed(BranchingRatio_n)
|
||
[unuc, i, j] = unique(NuclideInternal);
|
||
bratio = repmat(NaN, size(unuc));
|
||
for k = 1 : length(unuc);
|
||
if isempty(unuc{k})
|
||
% skip
|
||
elseif isfinite(ReferencePN(i(k)))
|
||
onuc = NuclideInternal{ReferencePN(i(k))};
|
||
if strcmp(unuc{k}, onuc)
|
||
bratio(k) = 1;
|
||
else
|
||
[s,t] = libBranchEquilib(sn, unuc{k}, onuc);
|
||
if s
|
||
bratio(k) = t;
|
||
else
|
||
% backwards compability
|
||
if isfield(nucIdxStruc, unuc{k})
|
||
bratio(k) = nucIdxStruc.(unuc{k}).BranchingRatio;
|
||
end
|
||
end
|
||
end
|
||
else
|
||
% skip
|
||
end
|
||
end
|
||
BranchingRatio = bratio(j);
|
||
end
|
||
|
||
if needed(NetAreaCalculated_n)
|
||
NetAreaCalculated = calcNetArea(ReferencePN, NetArea, ...
|
||
BranchingRatio, Yield, Efficiency, CCF);
|
||
end
|
||
|
||
if needed(RefEfficiencyErr_n)
|
||
RefEfficiencyErr = lCalcRefEfficiencyErr(RefEfficiencyErrPAT, ...
|
||
EfficiencyErr, ReferencePN);
|
||
end
|
||
|
||
if needed(NetAreaCalcError_n)
|
||
NetAreaCalcError = lCalcNetAreaCalcError(ReferencePN, AreaErrorP, ...
|
||
YieldErr, RefEfficiencyErr, CCFErr);
|
||
end
|
||
|
||
if needed(NetAreaMin_n)
|
||
NetAreaMin = max(max(NetArea, LC), 1);
|
||
end
|
||
|
||
if needed(NACalcMin_n)
|
||
NACalcMin = calcNetArea(ReferencePN, NetAreaMin, ...
|
||
BranchingRatio, Yield, Efficiency, CCF);
|
||
end
|
||
|
||
if needed(NetAreaCalcUTest_n)
|
||
NetAreaCalcUTest = abs(NetAreaCalculated - NetArea) ./ ...
|
||
sqrt((NACalcMin .* NetAreaCalcError/100) .^ 2 + AreaError.^2);
|
||
end
|
||
|
||
if needed(PeakShare_n)
|
||
myNACalc = NetAreaCalculated;
|
||
bool = isnan(myNACalc);
|
||
if any(bool)
|
||
myNACalc(bool) = SumPeakArea(bool);
|
||
end
|
||
PeakShare = 100 * myNACalc ./ NetArea;
|
||
end
|
||
|
||
if any(needed([ HalfLife_n, HalfLifeStr_n, ...
|
||
EffectiveHalflife_n, EffectiveHalflifeStr_n]))
|
||
[unuc, i, j] = unique(NuclideInternal);
|
||
[hl, hlunc, hle] = libHalfLife('Full', unuc);
|
||
HalfLife = hl(j);
|
||
HalfLifeError = hle(j);
|
||
if needed(HalfLifeStr_n)
|
||
hlstr = repmat({''}, size(unuc));
|
||
for k = 1 : length(hl)
|
||
if isfinite(hl(k))
|
||
hlstr{k} = halfLifeStr(hl(k));
|
||
end
|
||
end
|
||
HalfLifeStr = hlstr(j);
|
||
end
|
||
if any(needed([ EffectiveHalflife_n, EffectiveHalflifeStr_n]))
|
||
lookup = getEquilibData(sn, 'EffectiveHalflifes');
|
||
for k = 1 : length(hl)
|
||
if isfield(lookup, unuc{k})
|
||
hl(k) = lookup.(unuc{k});
|
||
hle(k) = NaN;
|
||
end
|
||
end
|
||
EffectiveHalflife = hl(j);
|
||
EffectiveHalflifeError = hle(j);
|
||
if needed(EffectiveHalflifeStr_n)
|
||
hlstr = repmat({''}, size(unuc));
|
||
for k = 1 : length(hl)
|
||
if isfinite(hl(k))
|
||
hlstr{k} = halfLifeStr(hl(k));
|
||
end
|
||
end
|
||
EffectiveHalflifeStr = hlstr(j);
|
||
end
|
||
end
|
||
end
|
||
|
||
if needed(DecayCorrection_n)
|
||
% decay correction just for acquisition
|
||
t_real = dminfo(sn, 'AcquisitionReal')/(24*3600);
|
||
f_real = log(2)*t_real ./ HalfLife;
|
||
DecayCorrection = (1 - exp(-f_real)) ./ f_real;
|
||
end
|
||
|
||
if needed(Activity_n)
|
||
t_live = dminfo(sn, 'AcquisitionLive');
|
||
Activity = CCF .* NetArea ./ ( ...
|
||
Efficiency .* (Yield/100) .* DecayCorrection * t_live);
|
||
end
|
||
|
||
if needed(RecoilBetaChan_n)
|
||
RecoilBetaChan = RecoilBeta .* dE;
|
||
end
|
||
|
||
if needed(RecoilDeltaChan_n)
|
||
RecoilDeltaChan = RecoilDeltaE ./ dE;
|
||
end
|
||
|
||
for i=1:length(varargin)
|
||
try
|
||
varargout{i}=eval(varargin{i},'error(''foobar'')');
|
||
catch
|
||
logWrite('Descriptor [%s] not properly assigned', varargin{i});
|
||
varargout{i} = repmat(NaN, [NPeaks, 1]);
|
||
end
|
||
end
|
||
}
|