#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
}