#include "GammaAnalyALG.h"

#include <QtMath>
#include <QFile>
#include <QCoreApplication>
#include <QTextStream>
#include <QStringList>
#include "fitfunc.h"
#include "genenalfunc.h"
#include "matlab_func.h"
#include "IndependentAlg.h"
#include <QDateTime>
#include <QDomDocument>

#include <QJsonDocument>
#include <QJsonObject>
#include <QJsonArray>
//#define gNaN qSNaN()

using namespace arma;
void PrintMat22(QString name, mat& dataMat);
void PrintColvec(QString name, colvec& dataMat);
void PrintMat2fvec(QString name, field<vec>& vecc);
namespace AlgFunc
{
	void jsonArrToStdvec(QJsonObject jObj, QString strKey, stdvec& val)
	{
		QJsonArray usedArr = jObj.value(strKey).toArray();
		val.clear();
		for (int i = 0; i < usedArr.size(); i++)
		{
			val.push_back(usedArr.at(i).toDouble());
		}
	}
	void formJson(QString strJson, PHDFile *phd, rowvec& spec, vec& iC)
	{
		QJsonParseError err_rpt;
		QJsonDocument root_Doc = QJsonDocument::fromJson(strJson.toStdString().c_str(), &err_rpt);//字符串格式化为JSON
		if (err_rpt.error != QJsonParseError::NoError)
		{
			return;
		}
		QJsonObject rootObj = root_Doc.object();
		QJsonArray arr = rootObj.value("vPeak").toArray();
		phd->vPeak.clear();
		for (int i = 0; i < arr.size(); i++)
		{
			QJsonObject obj = arr[i].toObject();

			PeakInfo peak;
			peak.peakCentroid = obj.value("peakCentroid").toDouble();
			peak.fwhmc = obj.value("fwhmc").toDouble();
			peak.tail = obj.value("tail").toDouble();
			peak.tailAlpha = obj.value("tailAlpha").toDouble();
			peak.upperTail = obj.value("upperTail").toDouble();
			peak.upperTailAlpha = obj.value("upperTailAlpha").toDouble();
			peak.area = obj.value("area").toDouble();
			peak.stepRatio = obj.value("stepRatio").toDouble();

			phd->vPeak.push_back(peak);
		}

		jsonArrToStdvec(rootObj, "usedEnerPara", phd->usedEnerPara.p);
		jsonArrToStdvec(rootObj, "usedResoPara", phd->usedResoPara.p);

		jsonArrToStdvec(rootObj, "XCtrl", phd->baseCtrls.XCtrl);
		jsonArrToStdvec(rootObj, "YCtrl", phd->baseCtrls.YCtrl);
		jsonArrToStdvec(rootObj, "YSlope", phd->baseCtrls.YSlope);
		for (int i = 0; i < phd->baseCtrls.YSlope.size(); i++)
		{
			if (phd->baseCtrls.YSlope[i] == 0.0)
				phd->baseCtrls.YSlope[i] = qQNaN();
		}

		jsonArrToStdvec(rootObj, "para_tail", phd->para_tail.p);
		jsonArrToStdvec(rootObj, "para_tailAlpha", phd->para_tailAlpha.p);
		jsonArrToStdvec(rootObj, "para_tailRight", phd->para_tailRight.p);
		jsonArrToStdvec(rootObj, "para_tailRightAlpha", phd->para_tailRightAlpha.p);
		jsonArrToStdvec(rootObj, "para_stepRatio", phd->para_stepRatio.p);
		jsonArrToStdvec(rootObj, "usedEffiPara", phd->usedEffiPara.p);

		phd->baseCtrls.rg_low = rootObj.value("rg_low").toInt();
		phd->baseCtrls.rg_high = rootObj.value("rg_high").toInt();

		phd->Spec.num_g_channel = rootObj.value("num_g_channel").toInt();
		phd->Spec.begin_channel = rootObj.value("begin_channel").toInt();

		
		//spec = v;// .clear();
		stdvec v;
		QJsonArray usedArr = rootObj.value("spec").toArray();
		for (int i = 0; i < usedArr.size(); i++)
		{
			v.push_back(usedArr.at(i).toDouble());
		}
		spec = v;

		stdvec v1;
		QJsonArray iCArr = rootObj.value("iC").toArray();
		for (int i = 0; i < iCArr.size(); i++)
		{
			v1.push_back(iCArr.at(i).toDouble());
		}
		iC = v1;

	}

	void formJson2(QString strJson, PHDFile *phd, rowvec& spec, vec& iC)
	{
		QJsonParseError err_rpt;
		QJsonDocument root_Doc = QJsonDocument::fromJson(strJson.toStdString().c_str(), &err_rpt);//字符串格式化为JSON
		if (err_rpt.error != QJsonParseError::NoError)
		{
			return;
		}
		QJsonObject rootObj = root_Doc.object();
		phd->vPeak.clear();

		QJsonArray arr1 = rootObj.value("peakCentroid").toArray();
		QJsonArray arr2 = rootObj.value("fwhmc").toArray();
		QJsonArray arr3 = rootObj.value("tail").toArray();
		QJsonArray arr4 = rootObj.value("tailAlpha").toArray();
		QJsonArray arr5 = rootObj.value("upperTail").toArray();
		QJsonArray arr6 = rootObj.value("upperTailAlpha").toArray();
		QJsonArray arr7 = rootObj.value("area").toArray();
		QJsonArray arr8 = rootObj.value("stepRatio").toArray();

		for (int i = 0; i < arr1.size(); i++)
		{
			PeakInfo peak;
			peak.peakCentroid		= arr1.at(i).toDouble();
			peak.fwhmc				= arr2.at(i).toDouble();
			peak.tail				= arr3.at(i).toDouble();
			peak.tailAlpha			= arr4.at(i).toDouble();
			peak.upperTail			= arr5.at(i).toDouble();
			peak.upperTailAlpha		= arr6.at(i).toDouble();
			peak.area				= arr7.at(i).toDouble();
			peak.stepRatio			= arr8.at(i).toDouble();

			phd->vPeak.push_back(peak);
		}

		jsonArrToStdvec(rootObj, "usedEnerPara", phd->usedEnerPara.p);
		jsonArrToStdvec(rootObj, "usedResoPara", phd->usedResoPara.p);

		jsonArrToStdvec(rootObj, "XCtrl", phd->baseCtrls.XCtrl);
		jsonArrToStdvec(rootObj, "YCtrl", phd->baseCtrls.YCtrl);
		jsonArrToStdvec(rootObj, "YSlope", phd->baseCtrls.YSlope);

		jsonArrToStdvec(rootObj, "para_tail", phd->para_tail.p);
		jsonArrToStdvec(rootObj, "para_tailAlpha", phd->para_tailAlpha.p);
		jsonArrToStdvec(rootObj, "para_tailRight", phd->para_tailRight.p);
		jsonArrToStdvec(rootObj, "para_tailRightAlpha", phd->para_tailRightAlpha.p);
		jsonArrToStdvec(rootObj, "para_stepRatio", phd->para_stepRatio.p);
		jsonArrToStdvec(rootObj, "usedEffiPara", phd->usedEffiPara.p);

		phd->baseCtrls.rg_low = rootObj.value("rg_low").toInt();
		phd->baseCtrls.rg_high = rootObj.value("rg_high").toInt();

		phd->Spec.num_g_channel = rootObj.value("num_g_channel").toInt();
		phd->Spec.begin_channel = rootObj.value("begin_channel").toInt();


		//spec = v;// .clear();
		stdvec v;
		QJsonArray usedArr = rootObj.value("vCount").toArray();
		for (int i = 0; i < usedArr.size(); i++)
		{
			v.push_back(usedArr.at(i).toDouble());
		}
		spec = v;

		stdvec v1;
		QJsonArray iCArr = rootObj.value("iC").toArray();
		for (int i = 0; i < iCArr.size(); i++)
		{
			v1.push_back(iCArr.at(i).toDouble());
		}
		iC = v1;

	}


	QJsonArray toJsonArr(stdvec val)
	{
		QJsonArray arr;
		for (int i = 0; i < val.size(); i++)
		{
			arr.push_back(val[i]);
		}
		return arr;
	}
QString toJsonObj(PHDFile *phd, rowvec spec, vec iC)
	{
		QJsonObject retobj;

		QJsonArray arr;
		foreach(const PeakInfo &peak, phd->vPeak)
		{
			QJsonObject obj;
			obj.insert("peakCentroid", peak.peakCentroid);
			obj.insert("fwhmc", peak.fwhmc);
			obj.insert("tail", peak.tail);
			obj.insert("tailAlpha", peak.tailAlpha);
			obj.insert("upperTail", peak.upperTail);
			obj.insert("upperTailAlpha", peak.upperTailAlpha);
			obj.insert("area", peak.area);
			obj.insert("stepRatio", peak.stepRatio);
			arr.push_back(obj);
		}
		retobj.insert("vPeak", arr);

		QJsonArray usedEnerPara = toJsonArr(phd->usedEnerPara.p);
		QJsonArray usedResoPara = toJsonArr(phd->usedResoPara.p);
		retobj.insert("usedEnerPara", usedEnerPara);
		retobj.insert("usedResoPara", usedResoPara);

		QJsonArray XCtrl = toJsonArr(phd->baseCtrls.XCtrl);
		QJsonArray YCtrl = toJsonArr(phd->baseCtrls.YCtrl);
		QJsonArray YSlope = toJsonArr(phd->baseCtrls.YSlope);
		retobj.insert("XCtrl", XCtrl);
		retobj.insert("YCtrl", YCtrl);
		retobj.insert("YSlope", YSlope);

		QJsonArray para_tail = toJsonArr(phd->para_tail.p);
		QJsonArray para_tailAlpha = toJsonArr(phd->para_tailAlpha.p);
		QJsonArray para_tailRight = toJsonArr(phd->para_tailRight.p);
		QJsonArray para_tailRightAlpha = toJsonArr(phd->para_tailRightAlpha.p);
		QJsonArray para_stepRatio = toJsonArr(phd->para_stepRatio.p);
		QJsonArray usedEffiPara = toJsonArr(phd->usedEffiPara.p);

		retobj.insert("para_tail", para_tail);
		retobj.insert("para_tailAlpha", para_tailAlpha);
		retobj.insert("para_tailRight", para_tailRight);
		retobj.insert("para_tailRightAlpha", para_tailRightAlpha);
		retobj.insert("para_stepRatio", para_stepRatio);
		retobj.insert("usedEffiPara", usedEffiPara);

		retobj.insert("rg_low", phd->baseCtrls.rg_low);
		retobj.insert("rg_high", phd->baseCtrls.rg_high);

        retobj.insert("num_g_channel", (int)phd->Spec.num_g_channel);
        retobj.insert("begin_channel", (int)phd->Spec.begin_channel);

		QJsonArray specArr;
		for (int i = 0; i < spec.size(); i++)
		{
			specArr.push_back(spec[i]);
		}
		retobj.insert("spec", specArr);

		QJsonArray iCArr;
		for (int i = 0; i < iC.size(); i++)
		{
			iCArr.push_back(iC[i]);
		}
		retobj.insert("iC", iCArr);

		QJsonDocument doc;
		doc.setObject(retobj);

		return QString(doc.toJson());
	}
//
QVector<int> insertpeak(stdvec &vLeft, stdvec &vRight, PHDFile *phd, rowvec spec, vec iC, double minNA)
{
// 	QString str = toJsonObj(phd, spec, iC);
// 	formJson(str, phd, spec, iC);
// 
// 	QFile file("E:/AWork/QINGH/insertVal.txt");
// 	file.open(QIODevice::ReadOnly);
// 	QByteArray bytearray = file.readAll();
// 	formJson(QString(bytearray), phd, spec, iC);
// 
// 	rowvec dddy(phd->baseCtrls.YSlope);
// 	mat mt(dddy);
// 	uvec idxxxx = find(mt < 1000000.0);
// 	QString SSS = "";
// 	for (size_t i = 0; i < idxxxx.size(); i++)
// 	{
// 		if (SSS != "")
// 			SSS += ",";
// 		SSS += QString("%1").arg(idxxxx(i));
// 	}

    uvec idx;
    uword peakNum = phd->vPeak.size();
    vec pC(peakNum), pFC(peakNum), pT(peakNum), pUT(peakNum), pNA(peakNum), pSt(peakNum), pTa(peakNum), pUTa(peakNum), pBWW(peakNum), pRBC(peakNum), pRDC(peakNum);
    uword t_i = 0;
    foreach (const PeakInfo &peak, phd->vPeak)
    {
        pC(t_i) = peak.peakCentroid;
        pFC(t_i) = peak.fwhmc;
        pT(t_i) = peak.tail;
        pTa(t_i) = peak.tailAlpha;
        pUT(t_i) = peak.upperTail;
        pUTa(t_i) = peak.upperTailAlpha;
        pNA(t_i) = peak.area;
        pSt(t_i) = peak.area * peak.stepRatio;
        ++t_i;
    }
    pBWW.fill(gNaN);
    pRBC.fill(gNaN);
    pRDC.fill(gNaN);

    // get resolution in channels
    vec iE = Independ::calFcnEval(iC, phd->usedEnerPara.p);
    vec idE = Independ::calDerivEval(iC, phd->usedEnerPara.p);
    vec iR = Independ::calFcnEval(iE, phd->usedResoPara.p);
    vec iRC = iR / idE;

    //iE.print("iE = ");
    //idE.print("idE = ");
    //iR.print("iR = ");
    //iRC.print("iRC = ");
    // residual above 0
    int nCounts = phd->Spec.num_g_channel;
	if (phd->Spec.begin_channel == 0)
		spec(nCounts - 1) = gNaN;

    rowvec BaseChan = cumsum(ones<rowvec>(nCounts));
    PiecePoly bpp = Independ::makeBasePP(phd->baseCtrls.XCtrl, phd->baseCtrls.YCtrl, phd->baseCtrls.YSlope, pC, pFC, pSt);
    rowvec approx = Independ::peakBaseSpline(BaseChan, bpp, pNA, pC, pFC, pSt, pT, pTa, pUT, pUTa, pBWW, pRBC, pRDC);
    rowvec resid = spec - approx;

    uvec idx_resid = find(resid<0);
    resid(idx_resid) = zeros<rowvec>(idx_resid.size());

    // get net area estimate
    vec iNA(iC.size());
    for(uword i=0; i<iC.size(); ++i) //for i=1:length(iC);
    {
        // ch=[floor(iC(i)-iRC(i)): ceil(iC(i)+iRC(i))];
        urowvec ch = RangeVec(floor(iC(i)-iRC(i))-1, ceil(iC(i)+iRC(i))-1);
        //iNA(i) = eval('max(minNA,sum(resid(ch)))','minNA');
        iNA(i) = max(minNA, sum(resid(ch)));
    }

    vec vT = Independ::calFcnEval(iE, phd->para_tail.p);
    vec vTa = Independ::calFcnEval(iE, phd->para_tailAlpha.p);
    vec vUT = Independ::calFcnEval(iE, phd->para_tailRight.p);
    vec vUTa = Independ::calFcnEval(iE, phd->para_tailRightAlpha.p);
    vec vSr = Independ::calFcnEval(iE, phd->para_stepRatio.p);
    vec vEffi = Independ::calFcnEval(iE, phd->usedEffiPara.p);

    idx.set_size(iC.n_elem);
    for(uword i=0; i<iC.size(); ++i)
    {
        int j = 0;
        while (j < peakNum && iC(i) > phd->vPeak[j].peakCentroid) ++j;

        PeakInfo peak;
        peak.peakCentroid = iC(i);
        peak.energy = iE(i);
        peak.area = iNA(i);
        peak.sensitivity = gNaN;
        peak.fwhm = iR(i);
        peak.fwhmc = iRC(i);
        peak.stepRatio = vSr(i);
        peak.tail = vT(i);
        peak.tailAlpha = vTa(i);
        peak.upperTail = vUT(i);
        peak.upperTailAlpha = vUTa(i);
        peak.efficiency = vEffi(i);
        peak.BWWidthChan = 0;
        peak.recoilBetaChan = gNaN;
        peak.recoilDeltaChan = gNaN;

        phd->vPeak.insert(phd->vPeak.begin() + j, peak);

        idx(i) = j;
        pC.insert_rows(j, 1);
        pFC.insert_rows(j, 1);
        pC(j) = iC(i);
        pFC(j) = iRC(i);
    }

    mat vvRg = findPeakRange(idx, pC, pFC, phd->baseCtrls.rg_low, phd->baseCtrls.rg_high);
    vLeft = conv_to<stdvec>::from(vvRg.col(0));
    vRight = conv_to<stdvec>::from(vvRg.col(1));

    QVector<int> vRetIdx;
    for(uword i=0; i<idx.n_elem; ++i)
    {
        vRetIdx.push_back(idx(i));
    }
    return vRetIdx;
}

mat findPeakRange(uvec pnr, vec C, vec F_Ch, double rg_low, double rg_high)
{
    // get all peak centroids and FWHM in channel units
    //[C, F_Ch] = dmps(sn, pat, 'Centroid', 'FWHM_Ch');
    rowvec ar; //ar = dmspec(sn, 'AnalysisRange');
    ar << rg_low << rg_high;

    uword n = pnr.size();
    mat rg = zeros(n, 2);
    vec iL = zeros<vec>(n);
    vec iH = zeros<vec>(n);
    // left and right borders of multiplet
    for(uword i=0; i<n; ++i) //for i = 1 : n
    {
        //[rg(i, 1), rg(i, 2), iL(i), iH(i)] = findVPeakRange(pnr(i), C, F_Ch, ar, varargin{:});
        rowvec temp = findVPeakRange(pnr(i), C, F_Ch, ar);
        rg(i, 0) = temp(0);
        rg(i, 1) = temp(1);
        iL(i) = temp(2);
        iH(i) = temp(3);
    }
    return rg;
}

rowvec findVPeakRange(int pNr, vec C, vec FC, rowvec rg, double LD, double RD)
{
    int NP = C.size();

    // find leftmost peak in multiplet, and minimum left border
    int i = pNr;
    double l = C(i) - LD * FC(i);
    while (i > 0 && l < C(i-1) + RD * FC(i-1))
    {
        i = i - 1;
        l = min(l, C(i) - LD * FC(i));
    }

    // find rightmost peak in multiplet, and maximum right border
    int j = pNr;
    double h = C(j) + RD * FC(j);
    while (j < NP-1 && h > C(j+1) - LD * FC(j+1))
    {
        j = j + 1;
        h = max(h, C(j) + RD * FC(j));
    }

    // left and right borders of multiplet
    l = max(rg(0), floor(l));
    h = min(rg(1), ceil(h));

    rowvec results;
    results << l << h << i << j;
    return results;
}



void fitPeakFull(PHDFile *phd, uvec Af, uvec Cf, uvec Ff)
{
    int nCount = phd->Spec.num_g_channel;
    rowvec spec(nCount);
    for(int i=0; i<nCount; ++i)
    {
        spec(i) = phd->Spec.counts[i];
    }
    if(phd->Spec.begin_channel == 0)
    {
        spec.insert_cols(nCount, 1);
        spec(nCount) = gNaN;
        spec.shed_col(0);
    }

    // get index to all peaks where any parameter is free
    // get common min and max
    uvec allidx = join_vert( join_vert(Af, Cf), Ff );
    if(allidx.is_empty()) return;

    // get unique allidx
    ivec Bool = zeros<ivec>( max(allidx) + 1u);
    Bool(allidx).fill(1);
    allidx = find(Bool);

    int ip = 0, np = phd->vPeak.size();
    vec C(np), FC(np), NA(np), Sr(np), T(np), Ta(np), UT(np), UTa(np), BW(np), RB(np), RD(np);
    foreach (const PeakInfo &peak, phd->vPeak)
    {
        C(ip) = peak.peakCentroid;
        FC(ip) = peak.fwhmc;
        NA(ip) = peak.area;
        Sr(ip) = peak.stepRatio;
        T(ip) = peak.tail;
        Ta(ip) = peak.tailAlpha;
        UT(ip) = peak.upperTail;
        UTa(ip) = peak.upperTailAlpha;
        BW(ip) = peak.BWWidthChan;
        RB(ip) = peak.recoilBetaChan;
        RD(ip) = peak.recoilDeltaChan;
        ++ip;
    }

//    PrintColvec("na.txt", NA);

    rowvec cx = phd->baseCtrls.XCtrl;
    rowvec cy = phd->baseCtrls.YCtrl;
    rowvec cdy = phd->baseCtrls.YSlope;

    // get fitting region
    rowvec x = RangeVec2(phd->baseCtrls.rg_low, phd->baseCtrls.rg_high); // x = [rg(1) : rg(2)];

    rowvec t_x = x - 1;
    rowvec y = IdxMat(spec, t_x);

    /* % fit parameters
    [s, CXo, cy, cdy, Ao, Co, Fo, Sro, To, Tao, UTo, UTao, BWo, RBo, RDo, X2] = ...
       fitFunction(x, y, 'wrapStepRatio', ...
       cx, [], cy, CPf, cdy, [], NA, Af, C, Cf, FC, Ff, Sr, SRf, ...
       T, [], Ta, [], UT, [], UTa, [], BW, [], RB, [], RD, [], ...
       'Restrict', addArg{:});*/
    double X2;
    uvec ep;
    field<vec> Ps;
    field<uvec> free_idx;
    Ps << vectorise(cx) << vectorise(cy) << vectorise(cdy) << NA << C << FC << Sr << T << Ta << UT << UTa << BW << RB << RD;
//    PrintMat2fvec("Ps1.txt", Ps);
    free_idx << ep << ep << ep << Af << Cf << Ff << ep << ep << ep << ep << ep << ep << ep << ep;
    QString addArg = "Peaks above spline baseline";
    bool s = Independ::fitFunction(Ps, X2, free_idx, vectorise(x), vectorise(y), "wrapStepRatio", QStringList("Restrict"), addArg);
//    PrintMat2fvec("Ps1.txt", Ps);
    /*mat& CXo = Ps(0);
    //if(Ps(1).is_colvec()) cy = Ps(1).t();
    //else cy = Ps(1);
    //if(Ps(2).is_colvec()) cdy = Ps(2).t();
    //else cdy = Ps(2);*/
    mat& Ao = Ps(3);
//    PrintMat22("area.txt",Ao);
    mat& Co = Ps(4);
    mat& Fo = Ps(5);
    mat& Sro = Ps(6);
    mat& To = Ps(7);
    mat& Tao = Ps(8);
    mat& UTo = Ps(9);
    mat& UTao = Ps(10);
    mat& BWo = Ps(11);
    mat& RBo = Ps(12);
    mat& RDo = Ps(13);

    // only write to temp pat if fitting converged
    if(s)
    {
        double k = phd->setting.k_back;
        double k_alpha = phd->setting.k_alpha;
        double k_beta = phd->setting.k_beta;

        // get baseline and peak parameters
        rowvec BaseChan = RangeVec2(1, nCount);
        PiecePoly bpp = Independ::makeBasePP(vectorise(cx), vectorise(cy), vectorise(cdy), Co, Fo, Ao%Sro);
        rowvec BL = Matlab::ppval(BaseChan, bpp)+Independ::peakSteps(BaseChan, Co, Fo, Ao%Sro);
        phd->vBase = conv_to<stdvec>::from(BL);

        vec Sc = Fo/sqrt(8*log(2));

        /* // modify all peaks, if idx is string
        if isstr(idx)
           idx = [1 : length(Ct)];
        end*/
        if(Co.is_empty()) return; // no peaks in PAT or called needlessly
        if(allidx.is_empty())
        {
            allidx = vectorise( RangeVec(0, Co.size()-1) );
        }

        vec MBC, LC, LD, BC;
        // restrict to indexed peaks
        Co = Co(allidx); //Ct = Ct(idx);
        Fo = Fo(allidx); //Fc = Fc(idx);
        Ao = Ao(allidx); //NA = NA(idx);
//        PrintMat22("ao_allidx.txt",Ao);
        Sro = Sro(allidx);
        To = To(allidx);
        Tao = Tao(allidx);
        UTo = UTo(allidx);
        UTao = UTao(allidx);
        BWo = BWo(allidx);
        RBo = RBo(allidx);
        RDo = RDo(allidx);
        Sc = Sc(allidx); //Sc = Sc(idx);

        vec W = 2 * k * Fo; //W = 2 * k * Fc;

        // total background counts from Ct +- k*Fc
        BC = Independ::abkcnt( max(0, BL), Co, Fo, k); //BC = abkcnt(max(BL, 0), Ct, Fc, k);

        // Mean Back Counts
        MBC = BC / W;

        double c = 1.1;
        double f = 1.9562;
        vec DA = c*f*sqrt(Sc % MBC); //DA = c*f*sqrt(Sc .* MBC);

        // Currie's LC
        LC = k_alpha * DA;

        // Currie's LD
        // old formula: LD = k_alpha^2 + 2* LC;
        double kb2 = pow(k_beta, 2);
        // assymmetric formula, ref TBD
        LD = LC + (kb2/2) * ( 1 + sqrt(1 + (4/kb2)*( LC+pow(DA,2) ) ) );

        mat dE = Independ::calDerivEval(Co, phd->usedEnerPara.p);
        vec vE = Independ::calFcnEval(Co, phd->usedEnerPara.p);

        for(uword i=0; i<allidx.size(); ++i)
        {
            phd->vPeak[ allidx(i) ].peakCentroid = Co(i);
            phd->vPeak[ allidx(i) ].energy = vE(i);
            phd->vPeak[ allidx(i) ].area = Ao(i);
            phd->vPeak[ allidx(i) ].areaErr = sqrt( max( Ao(i), LC(i) ) + BC(i) );
            phd->vPeak[ allidx(i) ].fwhm = Fo(i) * dE(i);
            phd->vPeak[ allidx(i) ].fwhmc = Fo(i);
            phd->vPeak[ allidx(i) ].significance = Ao(i) / LC(i);
            phd->vPeak[ allidx(i) ].meanBackCount = MBC(i);
            phd->vPeak[ allidx(i) ].backgroundArea = BC(i);
            phd->vPeak[ allidx(i) ].lc = LC(i);
            phd->vPeak[ allidx(i) ].ld = LD(i);
            phd->vPeak[ allidx(i) ].tail = To(i);
            phd->vPeak[ allidx(i) ].tailAlpha = Tao(i);
            phd->vPeak[ allidx(i) ].upperTail = UTo(i);
            phd->vPeak[ allidx(i) ].upperTailAlpha = UTao(i);
            phd->vPeak[ allidx(i) ].stepRatio = Sro(i);
            phd->vPeak[ allidx(i) ].BWWidthChan = BWo(i);
            phd->vPeak[ allidx(i) ].recoilBetaChan = RBo(i);
            phd->vPeak[ allidx(i) ].recoilDeltaChan = RDo(i);
        }
    }
}

void UpdateBaseControl(BaseControls& bc)
{
    rowvec cx(bc.XCtrl);
    rowvec cy(bc.YCtrl);
    rowvec cdy(bc.YSlope);

    rowvec base = rowvec(bc.Baseline);
    rowvec sc = rowvec(bc.StepCounts);
    rowvec tmp_cx = cx - 1;
    rowvec cyr = cy - IdxMat(sc, tmp_cx); //cyr = cy - sc(cx);
    rowvec chans = RangeVec2(bc.rg_low, bc.rg_high);//chans = [rg(1):rg(2)];
    base.cols(bc.rg_low-1, bc.rg_high-1) = sc.cols(bc.rg_low-1, bc.rg_high-1) + Independ::wrapBaseSpline(chans, cx, cyr, cdy);

    bc.Baseline = conv_to<stdvec>::from(base);

    //hb = plot(e, bc, 'm-');
    //hc = plot(ce, cy, 'ms');
    //cpd.SlopeHandle = plot([ceLo; ceHi], [cy-2*cdy;cy+2*cdy], 'g+-');
    //uvec idx = find_finite(cdy); //idx = find(isfinite(cdy));
    //vec tmp1 = cy(idx), tmp2 = 2 * cdy(idx);
    //bc.SlopeHandle.clear();
    //bc.SlopeHandle = conv_to<stdvec>::from(join_vert(tmp1-tmp2,tmp1+tmp2));
}

stdvec calValuesOut(rowvec Chan, vec p)
{
    return conv_to<stdvec>::from(Independ::calFcnEval(Chan, p));
}

double calDerivaOut(double Chan, vec p)
{
    if(p.size() < 2) return 0;
    vec x; x << Chan;
    vec dy = Independ::calDerivEval(x, p);
    return as_scalar(dy);
}

ParameterInfo calFitPara(CalibrationType cal, int funcId, stdvec x, stdvec y, stdvec err)
{
    ParameterInfo param;
    int s = x.size();
    if(s != y.size())
    {
        qDebug() << "The size of param x and y must be same in function calFitPara";
        return param;
    }
    if(err.size() != s) err.assign(s, gNaN);

    vec p, perr;
    if(funcId == 1)
    {
        Independ::lMakeInterp(p, perr, x, y, err);
    }
    else {
        vec p0;
        p0 << funcId;
        p0 = join_vert(p0, Independ::calDefIniPara(funcId, x, y));
        //p0.print("p0 = ");
        Independ::calFitPara(p, perr, cal, x, y, err, p0);
    }

    if(Independ::calParaCheck(p))
    {
        param.p = conv_to<stdvec>::from(p);
        param.perr = conv_to<stdvec>::from(perr);
    }
    return param;
}

stdvec interp1(const PeakInfo &peak, vec regBase, vec u)
{
    vec x = vectorise(RangeVec2(peak.left, peak.right));

    PiecePoly pp;
    vec na, ct, fc, st, tail, ta, ut, uta, bw, rb, rd;
    na << peak.area;
    ct << peak.peakCentroid;
    fc << peak.fwhmc;
    st << 0;
    tail << peak.tail;
    ta << peak.tailAlpha;
    ut << peak.upperTail;
    uta << peak.upperTailAlpha;
    bw << peak.BWWidthChan;
    rb << peak.recoilBetaChan;
    rd << peak.recoilDeltaChan;
    vec y = regBase + vectorise(Independ::peakBaseSpline(x, pp, na, ct, fc, st, tail, ta, ut, uta, bw, rb, rd));

    stdvec v;
    //bool ix = 1; // Is x given as the first argument?

    // The v5 option, '*method', asserts that x is equally spaced.
    bool eqsp = false;

    QString extrapval = "extrap";

    // check for NaN's
    if (!x.is_finite())
    {
        qDebug() << "MATLAB:interp1:NaNinX','NaN is not an appropriate value for X.";
    }

    vec ua = u;

    // NANS are allowed as a value for F(X), since a function may be undefined
    // for a given value.
    if (!y.is_finite())
    {
        qDebug() << "MATLAB:interp1:NaNinY, NaN found in Y, interpolation at undefined values will result in undefined values.";
    }

    // Check dimensions.  Work with column vectors.
    //SizeCube siz_y = size(y); // obtaining the size of y in a vector
    uword m = y.n_rows;
    //uword n = y.n_cols;
    vec h;
    if (x.size() != m)
    {
        qDebug() << "MATLAB:interp1:YInvalidNumRows','Y must have length(X) rows.";
    }
    if (!eqsp)
    {
        h = diff(x);
        eqsp = (norm(diff(h), "inf") <= eps(norm(x, "inf")));
        if (!x.is_finite())
        {
            eqsp = 0; // if an INF in x, x is not equally spaced
        }
    }
    if (eqsp)
    {
        h << (x(m-1)-x(0)) / (m-1);
    }

    if (m < 2)
    {
        if (ua.is_empty())
        {
            return v;
        } else {
            qDebug() << "MATLAB:interp1:NotEnoughPts','There should be at least two data points.";
        }
    }

    //SizeCube siz = size(ua);
    vec p;

    // Interpolate
    rowvec t_y = y.t();
    //v = reshape(pchip(x,y,u),prod(siz_y(2:end)),prod(siz)).';
    v = conv_to<stdvec>::from(reshape(Independ::pchip(x, t_y, u), size(u)));

    return v;
}

bool CalculateMDCs(PHDFile *phd, NuclideActMda &nucActMda, int mainPeakIdx, double lambda, double keyLineYield, double CCF)
{
    //double Ld = 2 * Fwhmc * (Lcc - baseline) / 0.8591;
    //double MDA = Ld * CCF * DCF / (Yi * Effki * Tl);
    //double MDC = MDA / Svol;

    // 计算衰变校正因子——DCF
    QDateTime collectStart = QDateTime::fromString(phd->collect.collection_start_date + " " + phd->collect.collection_start_time, QString(DATATIME_FORMAT));
    QDateTime collectStop = QDateTime::fromString(phd->collect.collection_stop_date + " " + phd->collect.collection_stop_time, QString(DATATIME_FORMAT));
    QDateTime acqStart = QDateTime::fromString(phd->acq.acquisition_start_date + " " + phd->acq.acquisition_start_time, QString(DATATIME_FORMAT));
    double Ts = collectStart.secsTo(collectStop);   // 采样时间
    double Td = collectStop.secsTo(acqStart);       // 衰变时间
    double Ta = phd->acq.acquisition_real_time;     // 能谱获取实时间
    double Tl = phd->acq.acquisition_live_time;     // 能谱获取活时间
    double Svol = phd->collect.air_volume;          // 样品采样体积
    double DCF1, DCF2, DCF3;

    if ( Ts == 0 ) DCF1 = 1;
    else DCF1 = lambda * Ts / (1-exp(-lambda*Ts));
    if ( Td == 0 ) DCF2 = 1;
    else DCF2 = exp(lambda*Td);
    if ( Ta == 0 ) DCF3 = 1;
    else DCF3 = lambda * Ta / (1-exp(-lambda*Ta));

    double DCF_act = exp(lambda * phd->usedSetting.refTime_act.secsTo(acqStart));
    double DCF_conc = exp(lambda * phd->usedSetting.refTime_conc.secsTo(collectStart));
    //double DCF = DCF1 * DCF2 * DCF3;

    // calculate line activities
    const PeakInfo &peak = phd->vPeak[mainPeakIdx];
    double netKeyPeakArea = peak.area;
    double linePeakAreaErr = peak.areaErr;

    double detectorEfficiency = peak.efficiency;
    double detectorEfficiencyError = DETECTOR_EFFICIENCY_ERROR_RATIO * peak.efficiency;
    double fwhmc = peak.fwhmc;
    int index = (int)(peak.peakCentroid + 0.5) - 1;
    double lcc = phd->vLc[index];
    double baseline = phd->vBase[index];

    double activity, mda;
    if ( keyLineYield == 0 )
    {
        activity = 0;
        mda = 0;
    }
    else
    {
        activity = ( netKeyPeakArea * CCF * DCF3*DCF_act ) / ( keyLineYield * detectorEfficiency * Tl ) ;//Bq,参考时间为referenceTimeActivity
        mda = (2 * fwhmc * (lcc - baseline) / 0.8591 ) * CCF * DCF3*DCF_act / ( keyLineYield * detectorEfficiency * Tl ) ;//Bq,参考时间为referenceTimeActivity
    }
    double concentration = ( netKeyPeakArea * CCF*DCF1*DCF2*DCF3*DCF_conc ) / ( keyLineYield * detectorEfficiency * Tl ) *1e6/ Svol;//mBq/m3,参考时间为referenceTimeConc
    double mdc = (2 * fwhmc * (lcc - baseline) / 0.8591 ) * CCF*DCF1*DCF2*DCF3*DCF_conc / ( keyLineYield * detectorEfficiency * Tl ) *1e6 / Svol;
    double activityErr = sqrt( pow(linePeakAreaErr/netKeyPeakArea, 2) + pow(detectorEfficiencyError/detectorEfficiency, 2) ) * activity;
    //double concentrationErr = activityErr / Svol;

    nucActMda.activity = activity;
    nucActMda.act_err = activityErr;
    nucActMda.mda = mda;
    nucActMda.mdc = mdc;
    nucActMda.efficiency = peak.efficiency;
    nucActMda.effi_err = detectorEfficiencyError;
    nucActMda.concentration = concentration;

    return (activity > mda);
}

double CalculateMDC(PHDFile *phd, stdvec vMdcInfo, double CCF)
{
    if(vMdcInfo.size() < 3 || vMdcInfo[2] == 0) return 0;

    QDateTime collectStart = QDateTime::fromString(phd->collect.collection_start_date + " " + phd->collect.collection_start_time, QString(DATATIME_FORMAT));
    QDateTime collectStop = QDateTime::fromString(phd->collect.collection_stop_date + " " + phd->collect.collection_stop_time, QString(DATATIME_FORMAT));
    QDateTime acqStart = QDateTime::fromString(phd->acq.acquisition_start_date + " " + phd->acq.acquisition_start_time, QString(DATATIME_FORMAT));
    double Ts = collectStart.secsTo(collectStop);   // 采样时间
    double Td = collectStop.secsTo(acqStart);       // 衰变时间
    double Ta = phd->acq.acquisition_real_time;     // 能谱获取实时间
    double Tl = phd->acq.acquisition_live_time;     // 能谱获取活时间
    double Svol = phd->collect.air_volume;          // 样品采样体积
    double DCF1, DCF2, DCF3;

    double lambda = log(2.0) / (vMdcInfo[2] * 86400);
    if ( Ts == 0 ) DCF1 = 1;
    else DCF1 = lambda * Ts / (1-exp(-lambda*Ts));
    if ( Td == 0 ) DCF2 = 1;
    else DCF2 = exp(lambda*Td);
    if ( Ta == 0 ) DCF3 = 1;
    else DCF3 = lambda * Ta / (1-exp(-lambda*Ta));

    //double DCF = DCF1 * DCF2 * DCF3;
    //double DCF_act = exp(lambda * phd->setting.refTime_act.secsTo(acqStart));
    double DCF_conc = exp(lambda * phd->usedSetting.refTime_conc.secsTo(collectStart));
    // calculate line activities
    //PeakInfo &peak = phd->vPeak[mainPeakIdx];

    vec energy; energy << vMdcInfo[0];
    stdvec channel = energyToChannel(energy, phd->usedEnerPara.p);
    vec dE = Independ::calDerivEval(channel, phd->usedEnerPara.p);
    double fwhmc = as_scalar(Independ::calFcnEval(energy, phd->usedResoPara.p) / dE);
    double effi = as_scalar( Independ::calFcnEval(energy, phd->usedEffiPara.p) );

    size_t index = 0;
    for(size_t i=1; i<phd->vEnergy.size(); ++i)
    {
        if(phd->vEnergy[i] >= vMdcInfo[0])
        {
            index = i;
            if(phd->vEnergy[i] - vMdcInfo[0] > vMdcInfo[0] - phd->vEnergy[i-1]) index = i-1;
            break;
        }
    }
    double lcc = phd->vLc[index];
    double baseline = phd->vBase[index];

    /*double mda = 0;
    if ( vMdcInfo[1] != 0 )
    {
        mda = (2 * fwhmc * (lcc - baseline) / 0.8591 ) * CCF * DCF / ( vMdcInfo[1] * effi * Tl ) * 1e06;
        mda = (2 * fwhmc * (lcc - baseline) / 0.8591 ) * CCF * DCF3*DCF_act / ( vMdcInfo[1] * effi * Tl ) ;
    }
    double concentration = activity / Svol;
    double mdc = mda / Svol;*/
    double mdc = (2 * fwhmc * (lcc - baseline) / 0.8591 ) * CCF*DCF1*DCF2*DCF3*DCF_conc / ( vMdcInfo[1] * effi * Tl ) *1e6 / Svol;

    return mdc;
}

bool ReadQCLimit(QMap<QString, QcCheckItem>& qcItems, double &ener_Be7, stdvec &vMdcInfo, QString system_type, QString strFPath)
{
    qcItems.clear();
    vMdcInfo.clear();

    //QFile file(qApp->applicationDirPath() + "/SystemManager.xml");
	QFile file(strFPath + "/SystemManager.xml");
    if(!file.open(QIODevice::ReadOnly))
    {
        return false;
    }

    QDomDocument doc;
    if(!doc.setContent(&file))
    {
        file.close();
        return false;
    }
    file.close();

    QDomNode node = doc.documentElement().firstChild();
    while(!node.isNull())
    {
        if(node.nodeName() == "QCFlags-P" && system_type.toUpper() == "P")
        {
            QDomNode child = node.firstChild();
            while(!child.isNull())
            {
                QString name = child.nodeName();
                if(name == QC_COL_TIME || name == QC_ACQ_TIME || name == QC_DECAY_TIME || name == QC_SAMP_VOL
                        || name == QC_AIR_FLOW || name == QC_Be7_FWHM || name == QC_Ba140_MDC)
                {
                    qcItems[name].standard = child.toElement().attribute("green");
                }
                else if(name == "Be7")
                {
                    ener_Be7 = child.toElement().attribute("energy").toDouble();
                }
                else if(name == "Ba140")
                {
                    vMdcInfo.push_back(child.toElement().attribute("energy").toDouble());
                    vMdcInfo.push_back(child.toElement().attribute("yield").toDouble());
                    vMdcInfo.push_back(child.toElement().attribute("halflife").toDouble());
                }
                child = child.nextSibling();
            }
            break;
        }
        else if(node.nodeName() == "QCFlags-G" && system_type.toUpper() == "G")
        {
            QDomNode child = node.firstChild();
            while(!child.isNull())
            {
                QString name = child.nodeName();
                if(name == QC_COL_TIME || name == QC_ACQ_TIME || name == QC_DECAY_TIME || name == QC_SAMP_VOL
                        || name == QC_AIR_FLOW || name == QC_Xe133_MDC)
                {
                    qcItems[name].standard = child.toElement().attribute("green");
                }
                else if(name == "Xe133")
                {
                    vMdcInfo.push_back(child.toElement().attribute("energy").toDouble());
                    vMdcInfo.push_back(child.toElement().attribute("yield").toDouble());
                    vMdcInfo.push_back(child.toElement().attribute("halflife").toDouble());
                }
                child = child.nextSibling();
            }
            break;
        }
        node = node.nextSibling();
    }

    return true;
}

void RunQC(PHDFile *phd, QString strFPath)
{
    QDateTime start = QDateTime::fromString(phd->collect.collection_start_date + " " + phd->collect.collection_start_time, QString(DATATIME_FORMAT));
    QDateTime end = QDateTime::fromString(phd->collect.collection_stop_date + " " + phd->collect.collection_stop_time, QString(DATATIME_FORMAT));
    QDateTime acq = QDateTime::fromString(phd->acq.acquisition_start_date + " " + phd->acq.acquisition_start_time, QString(DATATIME_FORMAT));

    double collect_hour = start.secsTo(end) / 3600.0;
    double acq_hour = phd->acq.acquisition_real_time / 3600.0;
    double Decay_hour = end.secsTo(acq) / 3600.0;

    double ener_Be7;
    stdvec vMdcInfo;
    QMap<QString, QcCheckItem> qcItems;
    if(!ReadQCLimit(qcItems, ener_Be7, vMdcInfo, phd->header.system_type, strFPath))
    {
        qDebug() << "Read QC Flags from SystemManager.xml Failed!";
    }

    qcItems[QC_COL_TIME].value = collect_hour;
    qcItems[QC_ACQ_TIME].value = acq_hour;
    qcItems[QC_DECAY_TIME].value = Decay_hour;
    qcItems[QC_SAMP_VOL].value = phd->collect.air_volume;
    qcItems[QC_AIR_FLOW].value = phd->collect.air_volume / collect_hour;

    if(phd->isValid && phd->vBase.size() == phd->Spec.num_g_channel)
    {
        if(phd->header.system_type == "P")
        {
            vec energy; energy << ener_Be7;
            vec fwhm = Independ::calFcnEval(energy, phd->usedResoPara.p);
            qcItems[QC_Be7_FWHM].value = fwhm(0);
            qcItems[QC_Ba140_MDC].value = CalculateMDC(phd, vMdcInfo, 1.0);
        }
        else qcItems[QC_Xe133_MDC].value = CalculateMDC(phd, vMdcInfo, 1.0);
    }

    for(QMap<QString, QcCheckItem>::iterator iter = qcItems.begin(); iter != qcItems.end(); ++iter)
    {
        if(iter.value().standard.isEmpty()) continue;
        QStringList lists = iter.value().standard.split(',', QString::SkipEmptyParts);
        bool bSatisfy = true;
        foreach(QString str, lists)
        {
            if(str.contains('-')) continue;
            else if(str.contains('('))
            {
                if(iter.value().value <= str.remove('(').toDouble()) { bSatisfy = false; break; }
            }
            else if(str.contains(')'))
            {
                if(iter.value().value >= str.remove(')').toDouble()) { bSatisfy = false; break; }
            }
            else if(str.contains('['))
            {
                if(iter.value().value < str.remove('[').toDouble()) { bSatisfy = false; break; }
            }
            else if(str.contains(']'))
            {
                if(iter.value().value > str.remove(']').toDouble()) { bSatisfy = false; break; }
            }
        }
        iter.value().bPass = bSatisfy;
    }
    phd->QcItems = qcItems;
}

void FilterNuclideLine(NuclideLines &lines, double lowE)
{
    stdvec vE = lines.vEnergy;
    for(size_t i=0, j=0; i<vE.size(); ++i)
    {
        if(vE[i] < lowE)
        {
            lines.fullNames.removeAt(j);
            lines.vEnergy.erase(lines.vEnergy.begin()+j);
            lines.vYield.erase(lines.vYield.begin()+j);
            lines.vUncertE.erase(lines.vUncertE.begin()+j);
            lines.vUncertY.erase(lines.vUncertY.begin()+j);
            if(i == lines.key_flag) lines.key_flag = -1;
            else if(i < lines.key_flag) --lines.key_flag;
        }
        else ++j;
    }
    if(lines.key_flag < 0 && lines.vEnergy.size() > 0)
    {
        stdvec &vY = lines.vYield;
        lines.maxYeildIdx = 0;
        double maxYield = vY[0];
        for(size_t ii=1; ii<vY.size(); ++ii)
        {
            if(vY[ii] > maxYield)
            {
                maxYield = vY[ii];
                lines.maxYeildIdx = ii;
            }
        }
    }
    else lines.maxYeildIdx = lines.key_flag;
}

void ReadSpecialNuclides(QMap<QString, double> &mapHalflife, QStringList &vNuclides,QString xmlPath = "")
{
    QString fileName = xmlPath + "/nuclide_ActMdc.txt";
    QFile t_file(fileName);
    if(!t_file.open(QIODevice::ReadOnly))
    {
        return;
    }

    QTextStream in(&t_file);
    while(!in.atEnd())
    {
        QString line = in.readLine();
        if(line.contains("#MDA"))
        {
            line = in.readLine();
            while(!line.isEmpty() && !line.contains('#'))
            {
                QStringList strList = line.split(QRegExp("\\s"), QString::SkipEmptyParts);
                if(strList.size() == 3)
                {
                    mapHalflife[ strList[0] ] = strList[2].toDouble() * 86400;
                }
                line = in.readLine();
            }
        }
        if(line.contains("#Identify"))
        {
            line = in.readLine();
            while(!line.isEmpty() && !line.contains('#'))
            {
                QStringList strList = line.split(QRegExp("\\s"), QString::SkipEmptyParts);
                vNuclides.append(strList);
                line = in.readLine();
            }
        }
    }
}

void NuclidesIdent(PHDFile* phd, QMap<QString, NuclideLines> &mapLines, QString xmlPath)
{
    // 当重新分析时先清除上一次的分析结果
    phd->mapNucActMda.clear();
    for(size_t i=0; i<phd->vPeak.size(); ++i) phd->vPeak[i].nuclides.clear();

    // 获取需特殊处理的核素
    QMap<QString, double> mapHalflife; // 用其他核素半衰期计算活度/浓度的核素
    QStringList vNuclides;             // 只识别不计算活度/浓度的核素
    ReadSpecialNuclides(mapHalflife, vNuclides, xmlPath);

    double energyWidth = phd->usedSetting.EnergyTolerance;
    int peakNum = phd->vPeak.size();
    for (QMap<QString, NuclideLines>::iterator iter = mapLines.begin(); iter != mapLines.end(); ++iter)
    {
        if(iter.value().halflife <= 0) continue;
        FilterNuclideLine(iter.value(), phd->usedSetting.ECutAnalysis_Low); // 过滤核素能量小于ECutAnalysis_Low的射线

        const stdvec &vEnergy = iter.value().vEnergy; // 该核素的所有γ射线能量
        const stdvec &vYield = iter.value().vYield;
        stdvec vEffi = calValuesOut(vEnergy, phd->usedEffiPara.p); // 该核素所有γ射线能量处的探测效率
        std::vector<int> vFit(vEnergy.size(), -1); // γ射线能量与峰中心道能量匹配标识

        double sum_total = 0; // 该核素所有γ射线能量处效率乘以分支比的和
        double sum_found = 0; // 所有匹配的γ射线能量处效率乘以分支比的和
        int mainPeakIdx = -1; // 记录核素主γ峰的索引下标
        for (size_t i=0, j=0; i<vEnergy.size(); ++i)
        {
            for(; j<peakNum; ++j)
            {
                if(phd->vPeak[j].energy >= 510 && phd->vPeak[j].energy <= 512)
                {
                    continue; // 峰中心道能量为511的峰不进行核素识别
                }
                if(vEnergy[i] < phd->vPeak[j].energy - energyWidth) break;
                else if(vEnergy[i] <= phd->vPeak[j].energy + energyWidth)
                {
                    sum_found += vEffi[i] * vYield[i];
                    vFit[i] = j;
                    if(iter.value().maxYeildIdx == i) mainPeakIdx = j;
                    break;
                }
            }
            sum_total += vEffi[i] * vYield[i];
        }

        // 核素匹配到峰
        if(sum_total > 0)
        {
            // 如果该核素属特殊核素，则用“特殊核素配置文件”中指明的其他核素的半衰期
            double halflife = (mapHalflife.find(iter.key()) == mapHalflife.end() ? iter.value().halflife : mapHalflife[iter.key()]);
            double lambda = log(2.0) / halflife;
            double rate = sum_found / sum_total;
            if(rate > NUCL_IDENT_VALUE2)
            {
                // 获取用于计算Activity、MDC的主γ峰和最大分支比
                double maxFoundYield = vYield[iter.value().maxYeildIdx];
                if(mainPeakIdx < 0)
                {
                    maxFoundYield = 0;
                    for(size_t ii=0; ii<vFit.size(); ++ii)
                    {
                        if(vFit[ii] >= 0 && vYield[ii] > maxFoundYield)
                        {
                            mainPeakIdx = vFit[ii];
                            maxFoundYield = vYield[ii];
                        }
                    }
                    if(mainPeakIdx < 0) continue;
                }

                NuclideActMda ActMda;
                bool bActBigger = CalculateMDCs(phd, ActMda, mainPeakIdx, lambda, maxFoundYield, 1);

                if(rate > NUCL_IDENT_VALUE1 || bActBigger)
                {
                    ActMda.halflife = halflife;
                    ActMda.key_flag = -1;
                    if(!vNuclides.contains(iter.key())) // 需要计算活度浓度的核素
                    {
                        ActMda.bCalculateMDA = true;
                    }
                    else {
                        ActMda.activity = 0;
                        ActMda.act_err = 0;
                        ActMda.efficiency = 0;
                        ActMda.effi_err = 0;
                        ActMda.mda = 0;
                        ActMda.mdc = 0;
                        ActMda.concentration = 0;
                    }

                    int fitLineNum = 0, peakIdx = -1;
                    for(size_t ii=0; ii<vFit.size(); ++ii)
                    {
                        if(vFit[ii] >= 0)
                        {
                            // 向峰信息表中添加识别到的核素
                            if(vFit[ii] != peakIdx)
                            {
                                phd->vPeak[ vFit[ii] ].nuclides.push_back( iter.key() );
                                peakIdx = vFit[ii];
                            }
                            // 添加匹配的γ射线的信息
                            if(vFit[ii] == mainPeakIdx) ActMda.calculateIdx = fitLineNum;
                            if(iter.value().maxYeildIdx == ii && iter.value().key_flag >= 0) ActMda.key_flag = fitLineNum;
                            ActMda.vPeakIdx.push_back(peakIdx+1);
                            ActMda.fullNames.push_back(iter.value().fullNames[ii]);
                            ActMda.vEnergy.push_back(vEnergy[ii]);
                            ActMda.vUncertE.push_back(iter.value().vUncertE[ii]);
                            ActMda.vYield.push_back(vYield[ii]);
                            ActMda.vUncertY.push_back(iter.value().vUncertY[ii]);

                            ++fitLineNum;
                        }
                    }
                    phd->mapNucActMda[iter.key()] = ActMda;
                }
            }
        }
    }
}

void CalcNuclideMDA(PHDFile* phd, NuclideLines &lines, const QString &nucName, std::vector<int> &vPeakIdx)
{
    if(lines.halflife <= 0) return;

    // 过滤核素能量小于ECutAnalysis_Low的射线
    FilterNuclideLine(lines, phd->usedSetting.ECutAnalysis_Low);

    // 获取需特殊处理的核素
    QMap<QString, double> mapHalflife; // 用其他核素半衰期计算活度/浓度的核素
    QStringList vNuclides;             // 只识别不计算活度/浓度的核素
    ReadSpecialNuclides(mapHalflife, vNuclides);

    double energyWidth = phd->usedSetting.EnergyTolerance;
    const stdvec &vEnergy = lines.vEnergy; // 该核素的所有γ射线能量
    double maxYield = 0;
    int mainPeakIdx = -1; // 记录核素主γ峰的索引下标

    NuclideActMda ActMda;
    ActMda.halflife = (mapHalflife.find(nucName) == mapHalflife.end() ? lines.halflife : mapHalflife[nucName]);
    for (size_t i=0, j=0; i<vEnergy.size(); ++i)
    {
        for(; j<vPeakIdx.size(); ++j)
        {
            double energy = phd->vPeak.at(vPeakIdx[j]).energy;
            if(vEnergy[i] < energy - energyWidth) break;
            else if(vEnergy[i] <= energy + energyWidth)
            {
                ActMda.vEnergy.push_back(vEnergy[i]);
                ActMda.vUncertE.push_back(lines.vUncertE[i]);
                ActMda.vYield.push_back(lines.vYield[i]);
                ActMda.vUncertY.push_back(lines.vUncertY[i]);
                ActMda.fullNames.push_back(lines.fullNames[i]);
                ActMda.vPeakIdx.push_back(vPeakIdx[j]+1);
                if(lines.key_flag == i) ActMda.key_flag = ActMda.vEnergy.size()-1;
                break;
            }
        }
    }
    for(size_t i=0; i<ActMda.vYield.size(); ++i)
    {
        if(ActMda.vYield[i] > maxYield)
        {
            maxYield = ActMda.vYield[i];
            mainPeakIdx = ActMda.vPeakIdx[i]-1;
            ActMda.calculateIdx = i;
        }
    }
    if(mainPeakIdx < 0) return;

    // 如果该核素属特殊核素，则用“特殊核素配置文件”中指明的其他核素的半衰期
    double halflife = (mapHalflife.find(nucName) == mapHalflife.end() ? lines.halflife : mapHalflife[nucName]);
    double lambda = log(2.0) / halflife;

    CalculateMDCs(phd, ActMda, mainPeakIdx, lambda, maxYield, 1);

    ActMda.bCalculateMDA = true;
    phd->mapNucActMda[nucName] = ActMda;
}

bool LoadSpectrum(PHDFile *phd, QString fileContents)
{
	std::string s = arma_version::as_string();
    bool bRet = true;

    RadionuclideMessage sample_data;
    if(fileContents.isNull())
    {
        bRet = sample_data.AnalysePHD_File(phd->filepath);
    }
    else
    {
        bRet = sample_data.AnalysePHD_Msg(fileContents);
    }
    if(!bRet) return bRet;

    phd->msgInfo = sample_data.GetMessageInfo();

    QString block_name = QLatin1String("#Header");
    QVariant header = sample_data.GetBlockData(block_name);
    if (!header.isNull())
    {
        phd->header = header.value<HeaderBlock>();
    }

    block_name = QLatin1String("#Comment");
    QVariant comment = sample_data.GetBlockData(block_name);
    if(!comment.isNull())
    {
        phd->oriTotalCmt = comment.value<CommentBlock>();
        //phd->totalCmt = phd->oriTotalCmt;
    }

    block_name = QLatin1String("#Collection");
    QVariant collection = sample_data.GetBlockData(block_name);
    if (!collection.isNull())
    {
        phd->collect = collection.value<CollectionBlock>();
        if(!phd->collect.collection_start_time.contains('.'))
        {
            phd->collect.collection_start_time.append(".0");
        }
        if(!phd->collect.collection_stop_time.contains('.'))
        {
            phd->collect.collection_stop_time.append(".0");
        }
    }
    else {
        phd->collect.air_volume = 0.0;
    }

    block_name = QLatin1String("#Acquisition");
    QVariant acquisition = sample_data.GetBlockData(block_name);
    if (!acquisition.isNull())
    {
        phd->acq = acquisition.value<AcquisitionBlock>();
        if(!phd->acq.acquisition_start_time.contains('.'))
        {
            phd->acq.acquisition_start_time.append(".0");
        }
    }
    else {
        phd->acq.acquisition_live_time = 0.0;
        phd->acq.acquisition_real_time = 0.0;
    }

    block_name = QLatin1String("#Processing");
    QVariant process = sample_data.GetBlockData(block_name);
    if (!process.isNull())
    {
        phd->process = process.value<ProcessingBlock>();
    }
    else {
        phd->process.sample_volume_of_Xe = 0.0;
        phd->process.Xe_collection_yield = 0.0;
        phd->process.uncertainty_1 = 0.0;
        phd->process.uncertainty_2 = 0.0;
    }

    block_name = QLatin1String("#Sample");
    QVariant sample = sample_data.GetBlockData(block_name);
    if(!sample.isNull())
    {
        phd->sampleBlock = sample.value<SampleBlock>();
    }
    else {
        phd->sampleBlock.dimension_1 = 0.0;
        phd->sampleBlock.dimension_2 = 0.0;
    }

    block_name = QLatin1String("#Calibration");
    QVariant calibra = sample_data.GetBlockData(block_name);
    if (!calibra.isNull())
    {
        phd->calibration = calibra.value<CalibrationBlock>();
    }

    block_name = QLatin1String("#Certificate");
    QVariant certify = sample_data.GetBlockData(block_name);
    if(!certify.isNull())
    {
        phd->certificate = certify.value<CertificateBlock>();
    }

    block_name = QLatin1String("#g_Spectrum");
    QVariant var_spec = sample_data.GetBlockData(block_name);
    if(!var_spec.isNull())
    {
        phd->Spec = var_spec.value<G_SpectrumBlock>();
        int i=0;
        for(; i<phd->Spec.num_g_channel; ++i)
        {
            if(phd->Spec.counts[i] > 0) break;
        }
        if(i == phd->Spec.num_g_channel) phd->isValid = false;
    }

    block_name = QLatin1String("#g_Energy");
    QVariant g_ener = sample_data.GetBlockData(block_name);
    if(!g_ener.isNull())
    {
        phd->mapEnerKD[CalPHD] = g_ener.value<G_EnergyBlock>();
    }

    block_name = QLatin1String("#g_Resolution");
    QVariant g_reso = sample_data.GetBlockData(block_name);
    if(!g_reso.isNull())
    {
        phd->mapResoKD[CalPHD] = g_reso.value<G_ResolutionBlock>();
    }

    block_name = QLatin1String("#g_Efficiency");
    QVariant g_effi = sample_data.GetBlockData(block_name);
    if(!g_effi.isNull())
    {
        if(phd->header.system_type == 'C')
        {
            phd->nMapEffiKD[CalPHD] = g_effi.value<n_G_EfficiencyBlock>();
        }
        else {
            phd->mapEffiKD[CalPHD] = g_effi.value<G_EfficiencyBlock>();
        }
    }

    block_name = QLatin1String("#TotalEff");
    QVariant g_tote = sample_data.GetBlockData(block_name);
    if(!g_tote.isNull())
    {
        phd->mapTotEKD[CalPHD] = g_tote.value<TotaleffBlock>();
    }


    block_name = QLatin1String("#b_self_Attenuation");
    QVariant bsel = sample_data.GetBlockData(block_name);
    if(!bsel.isNull())
    {
        phd->bSelfAttenuation = bsel.value<BSelfAttenuationBlock>();
    }

    block_name = QLatin1String("#b-gEfficiency");
    QVariant bgeff = sample_data.GetBlockData(block_name);
    if(!bgeff.isNull() && phd->header.system_type != 'C')
    {
        phd->mapbgEffiKD[CalPHD] = bgeff.value<BG_EfficiencyBlock>();
    }

    block_name = QLatin1String("#b_Efficiency");
    QVariant beff = sample_data.GetBlockData(block_name);
    if(!beff.isNull() && phd->header.system_type == 'C')
    {
        phd->nMapbgEffiKD[CalPHD] = beff.value<NBG_EfficiencyBlock>();
    }

    // 初始化默认分析设置
    if(phd->header.system_type == "P")
    {
        phd->setting.ECutAnalysis_Low = 35.0;
        phd->setting.bUpdateCal = true;
    }
    phd->setting.refTime_conc = QDateTime::fromString(phd->collect.collection_start_date + " " + phd->collect.collection_start_time, QString(DATATIME_FORMAT));
    phd->setting.refTime_act = QDateTime::fromString(phd->acq.acquisition_start_date + " " + phd->acq.acquisition_start_time, QString(DATATIME_FORMAT));
    phd->usedSetting = phd->setting;

    phd->bAnalyed = false;
    phd->analy_start_time = QDateTime::currentDateTimeUtc().toString(DATATIME_FORMAT_SPACE_SECONDS);

    qDebug() << QString("Load %1:").arg(phd->filepath) << LINE_END
             << QString("\tChannel number: %1").arg(phd->Spec.num_g_channel) << LINE_END
             << QString("\tStartChannel: %1").arg(phd->Spec.begin_channel) << LINE_END
             << QString("\tEnergy Calibration number: %1").arg(phd->mapEnerKD[CalPHD].record_count) << LINE_END
             << QString("\tResolution Calibration number: %1").arg(phd->mapResoKD[CalPHD].record_count) << LINE_END
             << QString("\tEfficiency Calibration number: %1").arg(phd->mapEffiKD[CalPHD].record_count) << LINE_END
             << QString("\tTotal Efficiency Calibration number: %1").arg(phd->mapTotEKD[CalPHD].record_count) << LINE_END;

    return bRet;
}

QStringList UserNuclide(QString sample_type, QString filename)
{
    if(filename.isEmpty())
    {
        if(sample_type.toUpper() == "P")
        {
            filename = qApp->applicationDirPath() + "/setup/P_default.nuclide";
        }
        else if(sample_type.toUpper() == "G")
        {
            filename = qApp->applicationDirPath() + "/setup/G_default.nuclide";
        }
    }

    QStringList nuclideList;
    QFile file(filename);
    if(file.open(QIODevice::ReadOnly))
    {
        QTextStream in(&file);
        QString nuclide = in.readLine();
        while(!nuclide.isNull())
        {
            nuclideList.push_back(nuclide);
            nuclide = in.readLine();
        }
        file.close();
    }
    else
    {
        //QString logName = GammaLogInfo::GetSysLogName(Module_G_Interactive);
        //MyLog4qt::MakeNewFileLogger(logName);
        //MyLog4qt::SetLoggerText(logName, GammaLogInfo::AddTime("Open default user nuclide library failed When Analyse ") + m_phd->filepath);
    }

    return nuclideList;
}

void ComputePeakRange(std::vector<PeakInfo> &vPeaks, int n, vec c, vec fwhmch, vec tlo, vec thi)
{
    vec fL = 3 + tlo;
    vec fH = 3 + thi;
    vec l = max(1, floor(c - fL % fwhmch));
    vec h = min(n, ceil(c + fH % fwhmch));

    ivec Bool = zeros<ivec>(n+1);
    for(uword i=0; i<l.size(); ++i)
    {
        Bool.rows(l(i)-1, h(i)-1) = ones<ivec>(h(i)-l(i)+1);
    }
    Bool(0) = 0;
    Bool(n) = 0;
    ivec d = diff(Bool);

    uvec rg_low = find(d == 1);
    uvec rg_high = find(d == -1);

    // 求出每个峰的左右边界及重峰信息
    uvec idx;
    int multiPeaks = 0, L, R, size_idx;
    for(uword i=0; i<rg_low.n_rows; ++i)
    {
        L = rg_low(i) + 1; // Ranges返回的是数组下标，但道从1开始
        R = rg_high(i) + 1;
        idx = find(c >= L && c <= R);
        size_idx = idx.size();
        if(size_idx > 1)
        {
            ++multiPeaks;
            for(uword j=0; j<size_idx; ++j)
            {
                vPeaks[idx(j)].left = L;
                vPeaks[idx(j)].right = R;
                vPeaks[idx(j)].multiIndex = multiPeaks;
                //vPeaks[idx(j)].no_peaks = j+1;
                //vPeaks[idx(j)].num_peaks = size_idx;
            }
        } else {
            vPeaks[idx(0)].left = L;
            vPeaks[idx(0)].right = R;
            //vPeaks[idx(0)].no_peaks = 0;
            //vPeaks[idx(0)].num_peaks = 0;
            vPeaks[idx(0)].multiIndex = 0;
        }
    }
}

stdvec energyToChannel(rowvec e, vec p)
{
    // maximal iterations in newton search
    int maxiter = 100;
    // relative error tolerance
    double tol = 1e-6;

    // get linear approximation
    rowvec x_test("0, 8000"); // x_test = [0, 8000];
    rowvec y_test = Independ::calFcnEval(x_test, p);
    double d = y_test(0);
    rowvec k = diff(y_test) / diff(x_test);

//    y_test.print("y_test = ");
//    k.print("k = ");

    // find bad energies
    uvec k2 = find_nonfinite(e);
    k2.insert_rows(k2.n_elem, find(e<1));
    // calculate with 1 and assing NaN afterwards
    e(k2).fill(1);

    // first estimate for channel
    rowvec c = (e - d) / k;
    rowvec eApp = Independ::calFcnEval(c, p);
    uvec idx = find(abs(eApp - e) > tol*e);
    //idx.print("idx = ");

    int iter = 1;
    while (iter < maxiter && !idx.is_empty())
    {
       iter = iter + 1;

       // newton step
       rowvec dE = Independ::calDerivEval(c(idx), p);

//       c.print("c = ");
//       e.print("e = ");
//       eApp.print("eApp = ");
//       dE.print("dE = ");

       c(idx) = c(idx) + (e(idx) - eApp(idx)) / dE;

       // check converged
       eApp(idx) = Independ::calFcnEval(c(idx), p);
       uvec subidx = find( abs(eApp(idx) - e(idx)) > tol*e(idx) );
       idx = idx(subidx);
    }

    if (!idx.is_empty())
    {
       qDebug() << QString("Could not find Channel for %1 energies").arg(idx.n_elem);
       c(idx).fill(gNaN);
    }

    // bad Energies
    c(k2).fill(gNaN);

    return conv_to<stdvec>::from(c);
}

void UpdateEfficiency(PHDFile *phd)
{
    size_t peakNum = phd->vPeak.size();
    stdvec vEner(peakNum);
    for (size_t i=0; i<peakNum; ++i)
    {
        vEner[i] = phd->vPeak[i].energy;
    }

    stdvec vEffi = calValuesOut(vEner, phd->usedEffiPara.p);
    for (size_t i=0; i<peakNum; ++i)
    {
        phd->vPeak[i].efficiency = vEffi[i];
    }
}

bool WriteBaseInfo(BaseControls &baseCtrl, const QString &filename)
{
    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly | QIODevice::Text)) return false;

    int numPerLine = 5;

    QTextStream out(&file);
    out.setFieldWidth(15);
    out.setFieldAlignment(QTextStream::AlignLeft);
    out.setPadChar(' ');

    out << "#AnalyseRange" << LINE_END;
    out << baseCtrl.rg_low << baseCtrl.rg_high << LINE_END;

    size_t i, nCP = baseCtrl.XCtrl.size(), nGroupCP = nCP / numPerLine * numPerLine;

    out << "#XCtrl" << LINE_END;
    stdvec &cx = baseCtrl.XCtrl;
    for(i=0; i<nGroupCP; i+=numPerLine)
    {
        out << cx[i] << cx[i+1] << cx[i+2] << cx[i+3] << cx[i+4] << LINE_END;
    }
    if(i < nCP)
    {
        for(; i<nCP; ++i) out << cx[i];
        out << LINE_END;
    }

    out << "#YCtrl" << LINE_END;
    stdvec &cy = baseCtrl.YCtrl;
    for(i=0; i<nGroupCP; i+=numPerLine)
    {
        out << cy[i] << cy[i+1] << cy[i+2] << cy[i+3] << cy[i+4] << LINE_END;
    }
    if(i < nCP)
    {
        for(; i<nCP; ++i) out << cy[i];
        out << LINE_END;
    }

    out << "#YSlope" << LINE_END;
    stdvec &cdy = baseCtrl.YSlope;
    for(i=0; i<nGroupCP; i+=numPerLine)
    {
        out << cdy[i] << cdy[i+1] << cdy[i+2] << cdy[i+3] << cdy[i+4] << LINE_END;
    }
    if(i < nCP)
    {
        for(; i<nCP; ++i) out << cdy[i];
        out << LINE_END;
    }

    size_t nBL = baseCtrl.Baseline.size(), nGroupBL = nBL / numPerLine * numPerLine;

    out << "#Baseline" << LINE_END;
    out << nBL << LINE_END;
    stdvec &bl = baseCtrl.Baseline;
    for(i=0; i<nGroupBL; i+=numPerLine)
    {
        out << bl[i] << bl[i+1] << bl[i+2] << bl[i+3] << bl[i+4] << LINE_END;
    }
    if(i < nBL)
    {
        for(; i<nBL; ++i) out << bl[i];
        out << LINE_END;
    }

    out << "#StepCounts" << LINE_END;
    out << nBL << LINE_END;
    stdvec &sc = baseCtrl.StepCounts;
    for(i=0; i<nGroupBL; i+=numPerLine)
    {
        out << sc[i] << sc[i+1] << sc[i+2] << sc[i+3] << sc[i+4] << LINE_END;
    }
    if(i < nBL)
    {
        for(; i<nBL; ++i) out << sc[i];
        out << LINE_END;
    }

    //Add By XYL 20170927
    file.close();

    return true;
}

bool WriteLcScac(const stdvec &vData, QString type, const QString &filename)
{
    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly | QIODevice::Text)) return false;

    int numPerLine = 5;

    QTextStream out(&file);
    out.setFieldWidth(15);
    out.setFieldAlignment(QTextStream::AlignLeft);
    out.setPadChar(' ');

    out << QString("#%1").arg(type) << LINE_END;
    out << vData.size() << LINE_END;

    size_t i, n = vData.size(), nGroup = n / numPerLine * numPerLine;
    for(i=0; i<nGroup; i+=numPerLine)
    {
        out << vData[i] << vData[i+1] << vData[i+2] << vData[i+3] << vData[i+4] << LINE_END;
    }
    if(i < n)
    {
        for(; i<n; ++i) out << vData[i];
        out << LINE_END;
    }

    //Add By XYL 20170927
    file.close();

    return true;
}

QString EquationDescription(int funcId)
{
    QString desc;
    switch (funcId)
    {
        case 1: desc = "y=yi+(y(i+1)-yi)*(x-xi)/(x(i+1)-xi) for xi<=x<x(i+1)"; break;
        case 2: desc = "y=a0+a1*x+a2*x^2+a3*x^3"; break;
        case 3: desc = "y=a0+a1*x^(1/2)+a2*x+a3*x^(3/2)"; break;
        case 4: desc = "y=sqrt(a0+a1*x+a2*x^2+a3*x^3)"; break;
        case 5: desc = "y=A*exp(-(E1/x)^k)*(1-exp(-(E2/x)^n))"; break;
        case 6: desc = "log(y)=a0+a1*log(x)+a2*log(x)^2+a3*log(x)^3"; break;
        case 7: desc = "log(y)=a0*x+a1+a2/x+a3/(x^2)"; break;
        case 8: desc = "log(y)=a0+a1*log(c/x)+a2*log(c/x)^2+a3*log(c/x)^3+a4*log(c/x)^4"; break;
        case 9: desc = "y=1/(a*x^(-k)+b*x^(-n))"; break;
        case 93: desc = "y=S*exp(-(E1/x)^k)*(1-exp(-(2*E3/(x-E3))^n"; break;
        case 94: desc = "y=S*exp(-(E1/x)^k)*(1-exp(-b*(1/(x-E2))^m))"; break;
        case 95: desc = "y=S*exp(-(E1/x)^k)*(1-exp(-b*(1/(x-E2))^m))*(1-exp(-(2*E3/(E-E3))^n))"; break;
        case 96: desc = "y=A+E1*x^(-n)"; break;
        case 97: desc = "y=a0+a1*exp(-lam1*x)"; break;
        case 98: desc = "y=c"; break;
        case 99: desc = "y=t_0+k_t*x if x>ECut"; break;
        default: desc = ""; break;
    }
    return desc;
}

QString EquationName(int funcId)
{
    QString name;
    switch (funcId)
    {
        case 1: name = "Interpolation"; break;
        case 2: name = "Polynomial"; break;
        case 3: name = "Square root polynomial"; break;
        case 4: name = "Square root of polynomial"; break;
        case 5: name = "HT Efficiency"; break;
        case 6: name = "Polynomial in log(y) against log(x)"; break;
        case 7: name = "Polynomial in log(y)"; break;
        case 8: name = "Polynomial in log(y) against log(1/x)"; break;
        case 9: name = "Inverse exponential"; break;
        case 93: name = "HAE Efficiency (1-3)"; break;
        case 94: name = "HAE Efficiency (1-2)"; break;
        case 95: name = "HAE Efficiency (1-2-3)"; break;
        case 96: name = "Inverse power"; break;
        case 97: name = "Exponential sum"; break;
        case 98: name = "Constant"; break;
        case 99: name = "Linear with cut-off"; break;
        default: name = ""; break;
    }
    return name;
}

}


void PrintMat2fvec(QString name, field<vec>& vecc)
{
    QFile file("C:/Result_Debug3/" + name);
    if (file.open(QIODevice::WriteOnly))
    {
        QTextStream out(&file);

        for (uword row = 0; row < vecc.size(); ++row)
        {
            for (uword col = 0; col < vecc(row).size(); ++col)
            {
                out << QString::number(vecc(row)(col), 'f', 20) << "\t";
            }
            out << "\r\n";
        }
        file.close();
    }
    else
    {
    }
}
void PrintMat22(QString name, mat& dataMat)
{
    QFile file("C:/Result_Debug3/" + name);
	if (file.open(QIODevice::WriteOnly))
	{
		QTextStream out(&file);
		if (dataMat.is_vec())
		{
			for (int i = 0; i < dataMat.size(); ++i)
			{
				out << QString::number(dataMat(i), 'g', 10) << "\r\n";
			}
		}
		else {
			for (uword row = 0; row < dataMat.n_rows; ++row)
			{
				for (uword col = 0; col < dataMat.n_cols; ++col)
				{
					out << QString::number(dataMat(row, col), 'g', 10) << "\t";
				}
				out << "\r\n";
			}
		}
		file.close();
	}
}

void PrintColvec(QString name, colvec& dataMat)
{
    QFile file("C:/Result_Debug3/" + name);
	if (file.open(QIODevice::WriteOnly))
	{
		QTextStream out(&file);
		for (int i = 0; i < dataMat.size(); ++i)
		{
			out << QString::number(dataMat(i), 'g', 10) << "\r\n";
		}
		
		file.close();
	}
}

void PrintFieldVec(QString name, field<vec>& vecc)
{
	QFile file(QCoreApplication::applicationDirPath() + "/Result_Debug3/" + name);
	if (file.open(QIODevice::WriteOnly))
	{
		QTextStream out(&file);

		for (uword row = 0; row < vecc.size(); ++row)
		{
			for (uword col = 0; col < vecc(row).size(); ++col)
			{
				out << QString::number(vecc(row)(col), 'g', 20) << "\t";
			}
			out << "\r\n";
		}
		file.close();
	}
}


/* ***********************************************************导出的类成员函数实现************************************************************ */
GammaAnalyALG::GammaAnalyALG()
{
    Init();
}

GammaAnalyALG::GammaAnalyALG(SpecSetup setup)
{
    Init();
    specSetup = setup;
}

GammaAnalyALG::GammaAnalyALG(PHDFile *phd)
{
    specSetup = phd->setting;
    SetBaseInfo(phd->Spec.counts, phd->Spec.begin_channel);

    QStringList names;
    names << phd->newEner << phd->newReso;
    SetCalPara(names, phd->mapEnerPara[names[0]], phd->mapResoPara[names[1]], phd->mapEffiPara[phd->newEffi], phd->mapTotEPara[phd->newTotE]);

    SetPAT(phd->vPeak);

    baseInfo(s_AnalysisRange) << phd->baseCtrls.rg_low << phd->baseCtrls.rg_high;
    baseInfo(s_XControl) = phd->baseCtrls.XCtrl;
    baseInfo(s_YControl) = phd->baseCtrls.YCtrl;
    baseInfo(s_YSlope) = phd->baseCtrls.YSlope;
    baseInfo(s_BaseLine) = phd->baseCtrls.Baseline;
    baseInfo(s_Steps) = phd->baseCtrls.StepCounts;
}

GammaAnalyALG::~GammaAnalyALG()
{
}

void GammaAnalyALG::Init()
{
	setJniBl(NULL, NULL, "");
	m_strFilePath = "";
    m_curEner = CalPHD;
    m_curReso = CalPHD;

    baseInfo.set_size(MAX_INDEX_SPEC);

    field<vec> f(2);
    f(0) << 2 << 0 << 0.3 << 0;
    f(1) = gNaN * ones<vec>(f(0).size()-1u);
    Acal_energy_para[CalPHD]  = f;

    f(0) << 4 << 1 << 0.003 << 0;
    f(1) = gNaN * ones<vec>(f(0).size()-1u);
    Acal_resolution_para[CalPHD] = f;

    f(0) << 94 << 0.1 << 30 << 3 << 50 << 70 << 0.8;
    f(1) = gNaN * ones<vec>(f(0).size()-1u);
    Acal_efficiency_para = f;
    Acal_tot_efficiency_para = f;

    Acal_tail_para << 99 << 0 << 0 << 0;
    Acal_tail_alpha_para << 98 << 1;
    Acal_tail_right_para << 98 << 0.35;
    Acal_tail_right_alpha_para << 98 << 1;
    Acal_step_ratio_para << 98 << 0;
}

bool GammaAnalyALG::Process(bool bUpdate, QChar dataType, vec certEne)
{
    if( !any(baseInfo(s_Spectrum) > 0) )
    {
        qDebug() << "The counts of spectrum are all zero!";
        return false;
    }
	
	callGammaProcess(0);

	stdvec vy;
	vy = getBaseInfo(s_YControl);

/* **********************************刻度更新*********************************************** */
    if(bUpdate)
    {
        calUpdate(dataType, certEne, true, true, true);
		callGammaProcess(1);
    }
	vy = getBaseInfo(s_YControl);
/* **********************************寻峰*********************************************** */
    PeakSearch();
	callGammaProcess(2);
    ////("____ps_peaks.txt", pat.Peaks);
	vy = getBaseInfo(s_YControl);
/* **********************************基线优化*********************************************** */
    baseImprove(specSetup.BaseImprovePSS);
	callGammaProcess(3);
    ////("____bi_peaks.txt", pat.Peaks);
	vy = getBaseInfo(s_YControl);
/* **********************************净面积拟合*********************************************** */
    vec t_idx = dmps(QStringList("NetAreaFree"))(0);
    //t_idx.print("t_idx");
    uvec idx = find(t_idx), ep;
	vy = getBaseInfo(s_YControl);
    //idx.print("idx = ");
    fitPeakFull(idx, ep, ep, ep);
	vy = getBaseInfo(s_YControl);
	callGammaProcess(4);

    baseInfo(s_Energy) = calValues(Cal_Energy, RangeVec2(1, m_nChans));
    baseInfo(s_FwhmcAll) = GetFwhmcAll();
    baseInfo(s_Lc) = calculateLC(baseInfo(s_BaseLine), baseInfo(s_FwhmcAll), specSetup.RiskLevelK);
    baseInfo(s_Scac) = calculateSCAC(baseInfo(s_Spectrum), baseInfo(s_BaseLine), baseInfo(s_FwhmcAll));
	vy = getBaseInfo(s_YControl);
    /*//("____na_peaks.txt", pat.Peaks);
    ////("____na_bl.txt", baseInfo(s_BaseLine));
    ////("____na_cx.txt", baseInfo(s_XControl));
    ////("____na_cy.txt", baseInfo(s_YControl));
    ////("____na_cdy.txt", baseInfo(s_YSlope));
    ////("____na_steps.txt", baseInfo(s_Steps));*/

    return true;
}

void GammaAnalyALG::callGammaProcess(int npro)
{
	if (m_pJniEnv == NULL || m_pSendObj == NULL)
		return;

	jclass jcls = m_pJniEnv->GetObjectClass(m_pSendObj);
 	if (jcls == NULL)
 	{
 		return;
 	}
 		
 	jmethodID jmth = m_pJniEnv->GetMethodID(jcls, "gammaProcess", "(Ljava/lang/String;Ljava/lang/String;)V");
 	if (jmth == NULL)
 	{
 		return;
 	}
 	jobject obj = m_pJniEnv->AllocObject(jcls);
 	if (obj == NULL)
 	{
 		return;
 	}

 	jstring juser = m_pJniEnv->NewStringUTF(m_strUserId.toStdString().c_str());

	char szbuf[32] = { 0 };
	itoa(npro, szbuf, 10);
	jstring jpro = m_pJniEnv->NewStringUTF(szbuf);
	m_pJniEnv->CallVoidMethod(obj, jmth, juser, jpro);
}

bool GammaAnalyALG::AnalyseSpectrum(PHDFile *phd, QMap<QString, NuclideLines> mapLines)
{
    // 分析处理调用
    //GammaAnalyALG alg(phd->setting);
    specSetup = phd->setting;
    SetBaseInfo(phd->Spec.counts, phd->Spec.begin_channel);
    QStringList names;
    names << phd->newEner << phd->newReso;
    bool bUpdate = false;
    if(phd->bAnalyed)
    {
        SetCalPara(names, phd->mapEnerPara[ names[0] ], phd->mapResoPara[ names[1] ], phd->mapEffiPara[phd->newEffi], phd->mapTotEPara[phd->newTotE]);
    }
    else {
        SetCalData(names, phd->mapEnerKD[ names[0] ], phd->mapResoKD[ names[1] ], phd->mapEffiKD[phd->newEffi], phd->mapTotEKD[phd->newTotE]);
        if(phd->setting.bUpdateCal) bUpdate = true;
    }
    if(!Process(bUpdate, phd->msgInfo.data_type.at(0), phd->certificate.g_energy.toStdVector())) return false;
    // 获取分析结果
    if(!phd->bAnalyed)
    {
        phd->bAnalyed = true;
//?
        if(!phd->mapEnerKD.isEmpty()) phd->mapEnerPara[phd->newEner] = GetCalPara(Cal_Energy, phd->newEner);
        if(!phd->mapResoKD.isEmpty()) phd->mapResoPara[phd->newReso] = GetCalPara(Cal_Resolution, phd->newReso);
        if(!phd->mapEffiKD.isEmpty()) phd->mapEffiPara[phd->newEffi] = GetCalPara(Cal_Efficiency);
        if(!phd->mapTotEKD.isEmpty()) phd->mapTotEPara[phd->newTotE] = GetCalPara(Cal_Tot_efficiency);
        phd->para_stepRatio = GetCalPara(Cal_Step_ratio);
        phd->para_tail = GetCalPara(Cal_Tail);
        phd->para_tailAlpha = GetCalPara(Cal_Tail_alpha);
        phd->para_tailRight = GetCalPara(Cal_Tail_right);
        phd->para_tailRightAlpha = GetCalPara(Cal_Tail_right_alpha);

        stdvec2 datas = GetCalData(Cal_Energy, CalUpdate1);
        if(datas.size() == 3)
        {
            G_EnergyBlock g_ener;
            g_ener.centroid_channel = QVector<double>::fromStdVector(datas[0]);
            g_ener.g_energy = QVector<double>::fromStdVector(datas[1]);
            g_ener.uncertainty = QVector<double>::fromStdVector(datas[2]);
            g_ener.record_count = g_ener.centroid_channel.size();

            phd->newEner = CalUpdate1;
            phd->mapEnerKD[CalUpdate1] = g_ener;
            phd->mapEnerPara[CalUpdate1] = GetCalPara(Cal_Energy, CalUpdate1);
        }

        datas = GetCalData(Cal_Energy, CalUpdate2);
        if(datas.size() == 3)
        {
            G_EnergyBlock g_ener;
            g_ener.centroid_channel = QVector<double>::fromStdVector(datas[0]);
            g_ener.g_energy = QVector<double>::fromStdVector(datas[1]);
            g_ener.uncertainty = QVector<double>::fromStdVector(datas[2]);
            g_ener.record_count = g_ener.centroid_channel.size();

            phd->newEner = CalUpdate2;
            phd->mapEnerKD[CalUpdate2] = g_ener;
            phd->mapEnerPara[CalUpdate2] = GetCalPara(Cal_Energy, CalUpdate2);
        }

        datas = GetCalData(Cal_Resolution, CalResUpdate);
        if(datas.size() == 3)
        {
            G_ResolutionBlock g_reso;
            g_reso.g_energy = QVector<double>::fromStdVector(datas[0]);
            g_reso.FWHM = QVector<double>::fromStdVector(datas[1]);
            g_reso.uncertainty = QVector<double>::fromStdVector(datas[2]);
            g_reso.record_count = g_reso.FWHM.size();

            phd->newReso = CalResUpdate;
            phd->mapResoKD[CalResUpdate] = g_reso;
            phd->mapResoPara[CalResUpdate] = GetCalPara(Cal_Resolution, CalResUpdate);
        }
    }
    phd->vEnergy = getBaseInfo(s_Energy);
    phd->vBase = getBaseInfo(s_BaseLine);
    phd->vLc = getBaseInfo(s_Lc);
    phd->vScac = getBaseInfo(s_Scac);

    phd->vPeak = getPAT();
    stdvec AnalyRg = getBaseInfo(s_AnalysisRange);
    if(AnalyRg.size() > 1)
    {
        phd->baseCtrls.rg_low = AnalyRg[0];
        phd->baseCtrls.rg_high = AnalyRg[1];
    }
    else {
        phd->baseCtrls.rg_low = 1;
        phd->baseCtrls.rg_high = phd->Spec.num_g_channel;
    }
    phd->baseCtrls.Baseline = phd->vBase;
    phd->baseCtrls.XCtrl = getBaseInfo(s_XControl);
    phd->baseCtrls.YCtrl = getBaseInfo(s_YControl);
    phd->baseCtrls.YSlope = getBaseInfo(s_YSlope);
    phd->baseCtrls.StepCounts = getBaseInfo(s_Steps);

    BaseCtrlStack bcStack;
    bcStack.cx = phd->baseCtrls.XCtrl;
    bcStack.cy = phd->baseCtrls.YCtrl;
    bcStack.cdy = phd->baseCtrls.YSlope;
    phd->baseCtrls.BaseStack.clear();
    phd->baseCtrls.BaseStack.push_back(bcStack);
    // 保存当前分析所用的参数信息
    phd->usedSetting = phd->setting;
    if(!phd->mapEnerKD.isEmpty())
    {
        phd->usedEner = phd->newEner;
        phd->usedEnerKD = phd->mapEnerKD[phd->newEner];
        phd->usedEnerPara = phd->mapEnerPara[phd->newEner];
    }
    if(!phd->mapResoKD.isEmpty())
    {
        phd->usedReso = phd->newReso;
        phd->usedResoKD = phd->mapResoKD[phd->newReso];
        phd->usedResoPara = phd->mapResoPara[phd->newReso];
    }
    if(!phd->mapEffiKD.isEmpty())
    {
        phd->usedEffi = phd->newEffi;
        phd->usedEffiKD = phd->mapEffiKD[phd->newEffi];
        phd->usedEffiPara = phd->mapEffiPara[phd->newEffi];
    }
    if(!phd->mapTotEKD.isEmpty())
    {
        phd->usedTotE = phd->newTotE;
        phd->usedTotEKD = phd->mapTotEKD[phd->newTotE];
        phd->usedTotEPara = phd->mapTotEPara[phd->newTotE];
    }
	callGammaProcess(5);
    AlgFunc::NuclidesIdent(phd, mapLines, this->m_strFilePath);
    AlgFunc::RunQC(phd, m_strFilePath);
	callGammaProcess(6);
    return true;
}

void GammaAnalyALG::SetSpecSetup(SpecSetup setup)
{
    specSetup = setup;
}

void GammaAnalyALG::SetBaseInfo(stdvec spec)
{
    m_nChans = spec.size();
    baseInfo(s_Spectrum) = spec;
    baseInfo(s_Approximation) = zeros<rowvec>(m_nChans);
    baseInfo(s_BaseLine) = zeros<rowvec>(m_nChans);
}

void GammaAnalyALG::SetBaseInfo(const QVector<long> &spec, int begin_channel)
{
    m_nChans = spec.size();
    baseInfo(s_Spectrum) = zeros<rowvec>(m_nChans);
    baseInfo(s_Approximation) = zeros<rowvec>(m_nChans);
    baseInfo(s_BaseLine) = zeros<rowvec>(m_nChans);
    for(int i=0; i<m_nChans; ++i)
    {
        baseInfo(s_Spectrum)[i] = spec[i];
    }
    if(begin_channel == 0)
    {
        baseInfo(s_Spectrum).insert_cols(m_nChans, 1);
        baseInfo(s_Spectrum)(m_nChans) = gNaN;
        baseInfo(s_Spectrum).shed_col(0);
    }
}

void GammaAnalyALG::SetBaseInfo(stdvec2 info)
{
    if(info.size() != MAX_INDEX_SPEC)
    {
        qDebug() << QString("The size of BaseInfo must be %1.").arg(MAX_INDEX_SPEC);
        return;
    }
    if(baseInfo.size() < MAX_INDEX_SPEC) baseInfo.set_size(MAX_INDEX_SPEC);
    for(int i=0; i<MAX_INDEX_SPEC; ++i)
    {
        baseInfo(i) = rowvec( info[i] );
    }
    m_nChans = baseInfo(s_Spectrum).size();
    /*baseInfo(s_BaseLine) = info.BaseLine;
    baseInfo(s_AnalysisRange) = info.AnalysisRange;
    baseInfo(s_XControl) = info.XControl;
    baseInfo(s_YControl) = info.YControl;
    baseInfo(s_YSlope) = info.YSlope;
    baseInfo(s_Residual) = info.Residual;
    baseInfo(s_Stripped) = info.Stripped;
    baseInfo(s_Steps) = info.Steps;
    baseInfo(s_ROISmooth) = info.ROISmooth;
    baseInfo(s_SmoothStripped) = info.SmoothStripped;
    baseInfo(s_Approximation) = info.Approx;
    baseInfo(s_ChiSquare) = info.ChiSquare;
    baseInfo(s_ArtXControl) = info.ArtXControl;
    baseInfo(s_ArtYControl) = info.ArtYControl;
    baseInfo(s_ArtYSlope) = info.ArtYSlope;
    baseInfo(s_ArtificialBase) = info.ArtificialBase;
    baseInfo(s_BaseFittingWeights) = info.BaseFittingWeights;
    baseInfo(s_PSS) = info.PSS;*/
}

stdvec2 GammaAnalyALG::getBaseInfo()
{
    stdvec2 info;
    for(uword i=0; i<baseInfo.size(); ++i)
    {
        info.push_back( conv_to< stdvec >::from( baseInfo(i) ) );
    }
    return info;
    /*AlgBaseInfo info;
    info.BaseLine = conv_to< vector<double> >::from( baseInfo(s_BaseLine) );
    info.AnalysisRange = conv_to< vector<double> >::from( baseInfo(s_AnalysisRange) );
    info.XControl = conv_to< vector<double> >::from( baseInfo(s_XControl) );
    info.YControl = conv_to< vector<double> >::from( baseInfo(s_YControl) );
    info.YSlope = conv_to< vector<double> >::from( baseInfo(s_YSlope) );
    info.Residual = conv_to< vector<double> >::from( baseInfo(s_Residual) );
    info.Stripped = conv_to< vector<double> >::from( baseInfo(s_Stripped) );
    info.Steps = conv_to< vector<double> >::from( baseInfo(s_Steps) );
    info.ROISmooth = conv_to< vector<double> >::from( baseInfo(s_ROISmooth) );
    info.SmoothStripped = conv_to< vector<double> >::from( baseInfo(s_SmoothStripped) );
    info.Approx = conv_to< vector<double> >::from( baseInfo(s_Approximation) );
    info.ChiSquare = conv_to< vector<double> >::from( baseInfo(s_ChiSquare) );
    info.ArtXControl = conv_to< vector<double> >::from( baseInfo(s_ArtXControl) );
    info.ArtYControl = conv_to< vector<double> >::from( baseInfo(s_ArtYControl) );
    info.ArtYSlope = conv_to< vector<double> >::from( baseInfo(s_ArtYSlope) );
    info.ArtificialBase = conv_to< vector<double> >::from( baseInfo(s_ArtificialBase) );
    info.BaseFittingWeights = conv_to< vector<double> >::from( baseInfo(s_BaseFittingWeights) );
    info.PSS = conv_to< vector<double> >::from( baseInfo(s_PSS) );
    return info;*/
}

stdvec GammaAnalyALG::getBaseInfo(SpecBaseIdx idx)
{
    return conv_to<stdvec>::from(baseInfo(idx));
}

void GammaAnalyALG::SetPAT(const std::vector<PeakInfo> &vPeaks)
{
    pat.Peaks.set_size(vPeaks.size(), MAX_INDEX_PEAK);
    pat.Peaks.fill(gNaN);

    uword i_row = 0;
    foreach(PeakInfo peak, vPeaks)
    {
        pat.Peaks(i_row, i_Centroid) = peak.peakCentroid;
        pat.Peaks(i_row, i_Energy) = peak.energy;
        pat.Peaks(i_row, i_FWHM) = peak.fwhm;
        pat.Peaks(i_row, i_FWHM_Ch) = peak.fwhmc;
        pat.Peaks(i_row, i_NetArea) = peak.area;
        pat.Peaks(i_row, i_AreaError) = peak.areaErr;
        pat.Peaks(i_row, i_Efficiency) = peak.efficiency;
        pat.Peaks(i_row, i_Significance) = peak.significance;
        pat.Peaks(i_row, i_Sensitivity) = peak.sensitivity;
        pat.Peaks(i_row, i_StepRatio) = peak.stepRatio;
        pat.Peaks(i_row, i_Tail) = peak.tail;
        pat.Peaks(i_row, i_TailAlpha) = peak.tailAlpha;
        pat.Peaks(i_row, i_UpperTail) = peak.upperTail;
        pat.Peaks(i_row, i_UpperTailAlpha) = peak.upperTailAlpha;
        pat.Peaks(i_row, i_LC) = peak.lc;
        pat.Peaks(i_row, i_LD) = peak.ld;
        pat.Peaks(i_row, i_MeanBackCount) = peak.meanBackCount;
        pat.Peaks(i_row, i_BackgroundArea) = peak.backgroundArea;
        ++i_row;
    }
}

std::vector<PeakInfo> GammaAnalyALG::getPAT()
{
    // 向 PAT 表中存在但尚未写入值的项写入值(这里只写入需要的）
    dmps();

    // pat表中没有的项只能通过字符串获取
    QStringList items;
    items << "BWWidthChan" << "RecoilBetaChan" << "RecoilDeltaChan";
    field<vec> f = dmps(items);

    std::vector<PeakInfo> vPI;
    for(uword row=0; row<pat.Peaks.n_rows; ++row)
    {
        PeakInfo pi;
        pi.index = row + 1;
        pi.peakCentroid = pat.Peaks(row, i_Centroid);
        pi.energy = pat.Peaks(row, i_Energy);
        pi.fwhm = pat.Peaks(row, i_FWHM);
        pi.fwhmc = pat.Peaks(row, i_FWHM_Ch);
        pi.area = pat.Peaks(row, i_NetArea);
        pi.areaErr = pat.Peaks(row, i_AreaError);
        pi.significance = pat.Peaks(row, i_Significance);
        pi.sensitivity = pat.Peaks(row, i_Sensitivity);
        pi.multiIndex = pat.Peaks(row, i_Multiplet);
        pi.efficiency = pat.Peaks(row, i_Efficiency);
        pi.stepRatio = pat.Peaks(row, i_StepRatio);
        pi.tail = pat.Peaks(row, i_Tail);
        pi.tailAlpha = pat.Peaks(row, i_TailAlpha);
        pi.upperTail = pat.Peaks(row, i_UpperTail);
        pi.upperTailAlpha = pat.Peaks(row, i_UpperTailAlpha);
        pi.lc = pat.Peaks(row, i_LC);
        pi.ld = pat.Peaks(row, i_LD);
        pi.meanBackCount = pat.Peaks(row, i_MeanBackCount);
        pi.backgroundArea = pat.Peaks(row, i_BackgroundArea);
        pi.BWWidthChan = f(0)(row);
        pi.recoilBetaChan = f(1)(row);
        pi.recoilDeltaChan = f(2)(row);
        vPI.push_back(pi);
    }

    // 求出每个峰的左右边界及重峰信息
    uvec idx;
    umat Ranges = lGetPlotRanges(m_nChans, pat.Peaks.col(i_Centroid), pat.Peaks.col(i_FWHM_Ch), pat.Peaks.col(i_Tail), pat.Peaks.col(i_UpperTail));
    int multiPeaks = 0, L, R, size_idx;
    for(uword i=0; i<Ranges.n_rows; ++i)
    {
        L = Ranges(i, 0) + 1; // Ranges返回的是数组下标，但道从1开始
        R = Ranges(i, 1) + 1;
        idx = find(pat.Peaks.col(i_Centroid) >= L && pat.Peaks.col(i_Centroid) <= R);
        size_idx = idx.size();
        if(size_idx > 1)
        {
            ++multiPeaks;
            for(uword j=0; j<size_idx; ++j)
            {
                vPI[idx(j)].left = L;
                vPI[idx(j)].right = R;
                vPI[idx(j)].multiIndex = multiPeaks;
                //vPI[idx(j)].no_peaks = j+1;
                //vPI[idx(j)].num_peaks = size_idx;
            }
        } else {
            vPI[idx(0)].left = L;
            vPI[idx(0)].right = R;
            vPI[idx(0)].multiIndex = 0;
            //vPI[idx(0)].num_peaks = 0;
        }
    }

    return vPI;
}

stdvec GammaAnalyALG::getPatInfo(PeakIndex idx)
{
    if( !pat.Peaks.col(idx).is_finite() )
    {
        return conv_to<stdvec>::from( dmps(idx) );
    }
    else return conv_to<stdvec>::from( pat.Peaks.col(idx) );
}

void GammaAnalyALG::setFilePath(QString strPath)
{
    m_strFilePath = strPath;
}

void GammaAnalyALG::setJniBl(JNIEnv* pEnv, jobject jobj, QString strUid)
{
	m_pJniEnv = pEnv;
	m_pSendObj = jobj;
	m_strUserId = strUid;
}

void GammaAnalyALG::SetCalData(CalibrationType cal, stdvec x, stdvec y, stdvec err, QString calName)
{
    int s = x.size();
    if(s != y.size())
    {
        qDebug() << "The size of x must be same as the size of y in function SetCalData.";
        return;
    }
    if(err.size() != s) { err.assign(s, gNaN); }

    mat calDatas(s, 3);
    calDatas.col(0) = vec(x);
    calDatas.col(1) = vec(y);
    calDatas.col(2) = vec(err);

    switch (cal)
    {
    case Cal_Energy: Cal_Ener_Data[calName] = calDatas; break;
    case Cal_Resolution: Cal_Reso_Data[calName] = calDatas; break;
    case Cal_Efficiency: Cal_Effi_Data = calDatas; break;
    case Cal_Tot_efficiency: Cal_Tot_Effi_Data = calDatas; break;
    default: break;
    }
}

void GammaAnalyALG::SetCalData(QStringList names, G_EnergyBlock EnerKD, G_ResolutionBlock ResoKD, G_EfficiencyBlock EffiKD, TotaleffBlock TotEKD)
{
    if(names.size() < 2)
    {
        qDebug() << "The parameter names must has at least 2 element in GammaAnalyALG::SetCalData.";
        return;
    }

    m_curEner = names[0];
    m_curReso = names[1];

    if(EnerKD.record_count > 0)
    {
        mat calDatas(EnerKD.record_count, 3);
        calDatas.col(0) = vec( EnerKD.centroid_channel.toStdVector() );
        calDatas.col(1) = vec( EnerKD.g_energy.toStdVector() );
        calDatas.col(2) = vec( EnerKD.uncertainty.toStdVector() );
        Cal_Ener_Data[m_curEner] = calDatas;

        Acal_energy_para[m_curEner] = FitCalPara(Cal_Energy, calDatas);
        //Acal_energy_para[m_curEner].print("Acal_energy_para = ");
    }
    if(ResoKD.record_count > 0)
    {
        mat calDatas(ResoKD.record_count, 3);
        calDatas.col(0) = vec( ResoKD.g_energy.toStdVector() );
        calDatas.col(1) = vec( ResoKD.FWHM.toStdVector() );
        calDatas.col(2) = vec( ResoKD.uncertainty.toStdVector() );
        Cal_Reso_Data[m_curReso] = calDatas;

        Acal_resolution_para[m_curReso] = FitCalPara(Cal_Resolution, calDatas);

		PrintFieldVec("Acal_resolution_para[m_curReso].txt", Acal_resolution_para[m_curReso]);

        //Acal_resolution_para[m_curReso].print("Acal_resolution_para = ");
    }
    if(EffiKD.record_count > 0)
    {
        mat calDatas(EffiKD.record_count, 3);
        calDatas.col(0) = vec( EffiKD.g_energy.toStdVector() );
        calDatas.col(1) = vec( EffiKD.efficiency.toStdVector() );
        calDatas.col(2) = vec( EffiKD.uncertainty.toStdVector() );
        Cal_Effi_Data = calDatas;

        Acal_efficiency_para = FitCalPara(Cal_Efficiency, Cal_Effi_Data);
        //Acal_efficiency_para.print("Acal_efficiency_para = ");
    }
    if(TotEKD.record_count > 0)
    {
        mat calDatas(TotEKD.record_count, 3);
        calDatas.col(0) = vec( TotEKD.g_energy.toStdVector() );
        calDatas.col(1) = vec( TotEKD.total_efficiency.toStdVector() );
        calDatas.col(2) = vec( TotEKD.uncertainty.toStdVector() );
        Cal_Tot_Effi_Data = calDatas;

        Acal_tot_efficiency_para = FitCalPara(Cal_Tot_efficiency, Cal_Tot_Effi_Data);
        //Acal_tot_efficiency_para.print("Acal_tot_efficiency_para = ");
    }

    /*QString name;
    setCalibrationData(name, Cal_Energy, "PHD", tmpEner);
    setCalibrationData(name, Cal_Resolution, "PHD", QVec2ToMat(ResoKD));
    setCalibrationData(name, Cal_Efficiency, "PHD", QVec2ToMat(EffiKD));
    setCalibrationData(name, Cal_Tot_efficiency, "PHD", QVec2ToMat(TotEKD)); */
}

stdvec2 GammaAnalyALG::GetCalData(CalibrationType cal, QString calName)
{
    stdvec2 calData;
    mat matData;
    switch (cal)
    {
        case Cal_Energy: matData = Cal_Ener_Data[calName]; break;
        case Cal_Resolution: matData = Cal_Reso_Data[calName]; break;
        case Cal_Efficiency: matData = Cal_Effi_Data; break;
        case Cal_Tot_efficiency: matData = Cal_Tot_Effi_Data; break;
    }
    if(matData.n_cols == 3)
    {
        calData.push_back( conv_to<stdvec>::from(matData.col(0u)) );
        calData.push_back( conv_to<stdvec>::from(matData.col(1u)) );
        calData.push_back( conv_to<stdvec>::from(matData.col(2u)) );
    }
    return calData;
}

void GammaAnalyALG::SetCalPara(CalibrationType cal, ParameterInfo para, QString calName)
{
    field<vec> fPara(2);
    fPara(0) = para.p;
    if(para.perr.empty()) { fPara(1) = gNaN * ones<vec>(fPara(0).size()-1); }
    else { fPara(1) = para.perr; }

    switch (cal) {
        case Cal_Energy: Acal_energy_para[calName] = fPara; break;
        case Cal_Resolution: Acal_resolution_para[calName] = fPara; break;
        case Cal_Efficiency: Acal_efficiency_para = fPara; break;
        case Cal_Tot_efficiency: Acal_tot_efficiency_para = fPara; break;
        case Cal_Tail: Acal_tail_para = vec(para.p); break;
        case Cal_Tail_alpha: Acal_tail_alpha_para = vec(para.p); break;
        case Cal_Tail_right: Acal_tail_right_para = vec(para.p); break;
        case Cal_Tail_right_alpha: Acal_tail_right_alpha_para = vec(para.p); break;
        case Cal_Step_ratio: Acal_step_ratio_para = vec(para.p);
        default: break;
    }
}

void GammaAnalyALG::SetCalPara(QStringList paraNames, ParameterInfo EnerPara, ParameterInfo ResoPara, ParameterInfo EffiPara, ParameterInfo TotEPara)
{
    if(paraNames.size() < 2)
    {
        qDebug() << "The parameter paraNames must has at least 2 element in GammaAnalyALG::SetCalPara.";
        return;
    }

    m_curEner = paraNames[0];
    m_curReso = paraNames[1];

    field<vec> f(2);
    f(0) = EnerPara.p;
    f(1) = EnerPara.perr;
    Acal_energy_para[m_curEner] = f;

    f(0) = ResoPara.p;
    f(1) = ResoPara.perr;
    Acal_resolution_para[m_curReso] = f;

    f(0) = EffiPara.p;
    f(1) = EffiPara.perr;
    Acal_efficiency_para = f;

    f(0) = TotEPara.p;
    f(1) = TotEPara.perr;
    Acal_tot_efficiency_para = f;
}

ParameterInfo GammaAnalyALG::GetCalPara(CalibrationType cal, QString calName)
{
    ParameterInfo param;
    field<vec> matPara;
    switch(cal)
    {
        case Cal_Energy: matPara = Acal_energy_para[calName]; break;
        case Cal_Resolution: matPara = Acal_resolution_para[calName]; break;
        case Cal_Efficiency: matPara = Acal_efficiency_para; break;
        case Cal_Tot_efficiency: matPara = Acal_tot_efficiency_para; break;
        case Cal_Step_ratio: matPara << Acal_step_ratio_para; break;
        case Cal_Tail: matPara << Acal_tail_para; break;
        case Cal_Tail_alpha: matPara << Acal_tail_alpha_para; break;
        case Cal_Tail_right: matPara << Acal_tail_right_para; break;
        case Cal_Tail_right_alpha: matPara << Acal_tail_right_alpha_para; break;
    }
    if(matPara.size() > 0)
    {
        param.p = conv_to<stdvec>::from(matPara(0));
        if(matPara.size() == 2) param.perr = conv_to<stdvec>::from(matPara(1));
    }
    return param;
}

field<vec> GammaAnalyALG::FitCalPara(CalibrationType cal, mat calData)
{
    field<vec> resField(2);
    if(calData.n_cols < 2) {
        qDebug() << "Calibration data need at least two columns!";
        return resField;
    }
    else if(calData.n_cols < 3)
    {
        calData.insert_cols(1, 2);
        calData.col(2).fill(gNaN);
    }

    vec p, perr;
    if(!Independ::calFitPara(p, perr, cal, calData.col(0u), calData.col(1u), calData.col(2u)))
    {
        bool b = Independ::lMakeInterp(p, perr, calData.col(0u), calData.col(1u), calData.col(2u));
        if(!b || p.size()<2 || !Independ::calParaCheck(p))
        {
            p.clear();
            perr.clear();
            switch (cal)
            {
                case Cal_Energy:
                {
                    p << 2 << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0;
                    if (m_nChans < 3000) p(2) = 1;
                    else if (m_nChans < 5000) p(2) = 0.5;
                    else if (m_nChans < 10000) p(2) = 0.33;
                    else p(2) = 0.2;
                    break;
                }
                case Cal_Resolution: p << 4 << 1 << 0.003 << 0; break;
                case Cal_Efficiency:
                case Cal_Tot_efficiency: p << 5 << 0.3 << 80 << 100 << 2.95 << 0.8; break;
                default: qDebug() << "invalid calibration string in calFitPara";
            }
        }
    }

    resField(0) = p;
    resField(1) = perr;
    return resField;
}

/* =====================================================寻峰模块======================================================== */
uvec GammaAnalyALG::PeakSearch(double PSSs, PeakSearchType PStype, rowvec ELim)
{
    vec C, NA, CNL, CNR, PSSfound;
    return PeakSearch(C, NA, CNL, CNR, PSSfound, PSSs, PStype, ELim);
}

uvec GammaAnalyALG::PeakSearch(vec &C, vec &NA, vec &CNL, vec &CNR,
          vec &PSSfound, double PSSs, PeakSearchType PStype, rowvec ELim)
{
    uvec idx;
    C.clear();
    NA.clear();
    CNL.clear();
    CNR.clear();
    PSSfound.clear();

    // lower cutoff energy
    double ECutLow = specSetup.ECutAnalysis_Low;
    double ECutHigh = qInf();//specSetup.ECutAnalysis_High;
    double deltaE = specSetup.EnergyTolerance;

    // least net area estimate, lower NA estimates will be increased
    double NAMin = 1; // min net area
    //double fROI = 2;
    double NWidth = 3;

    // get default variables for parameters
    if( qIsNaN(PSSs) ) PSSs = specSetup.PSS_low;
    //if(PSSb < 0.000001) PSSb = PSSs;

    // set parameter strings and flags for full or residual peak search
    bool fullSearch;
    char Src;
    colvec COld;    // 原来的Centroid
    QString type, approxType;
    if(PStype == search_residual)
    {
        fullSearch = false;
        type = "Stripped";
        approxType = "BaseLine";
        Src = 'R';
    }
    else {
        fullSearch = true;
        // full peak search: get old centroids, don't insert old peaks again
        if(pat.Peaks.size() > 0) COld = pat.Peaks.col(0); // COld = dmps("Centroid");
        type = "Spectrum";
        approxType = "Approximation";
        Src = 'M';
    }

    //get spectrum counts, analysis range and baseline control points
    //[spec, analysisRg, cx, cy, approx] = dmspec(sn, type, 'AnalysisRange', 'XControl', 'YControl', approxType);
    QStringList Opts;
    Opts << type << "AnalysisRange" << "XControl" << approxType;
    field<rowvec> t_f = dmspec(Opts);
    rowvec& spec = t_f(0);
    rowvec analysisRg = t_f(1);
    rowvec& cx = t_f(2);
    //cx.print("cx");
    //rowvec cy = t_f(3);
    rowvec& approx = t_f(3);
    //("peakSearch_spec.txt", spec);
    //("peakSearch_approx.txt", approx);

    // flag if we need to initialize
    bool doBaseInit;
    rowvec missedCounts;
    if(cx.is_empty()) doBaseInit = true;
    else {
        doBaseInit = false;
        missedCounts = spec - approx;
        NWidth = 0;
    }

    uword NChan = spec.size();
    rowvec Chan = RangeVec2(1, NChan);

    //get energy, resolution and minimum distance in channel units for every channel
    rowvec E = calValues(Cal_Energy, Chan);
    //("peakSearch_E.txt", E);

    rowvec dE = calDerivative(Cal_Energy, Chan);
    //("peakSearch_dE.txt", dE);

    rowvec res = calValues(Cal_Resolution, E);
    //("peakSearch_res.txt", res);

    if(res.size() != dE.size())
    {
        qDebug() << "The size of Energy Calibration must be Same as The size of Resolution Calibration!";
        return idx;
    }

    rowvec res_ch = res / dE;
    rowvec deltaC = deltaE / dE;
    ////("peakSearch_res_ch.txt", res_ch);
    ////("peakSearch_deltaC.txt", deltaC);

    // transform ELim to allowed range (channels)
    rowvec searchRg;
    if(!ELim.is_empty())
    {
        searchRg = lGetAnalysisRange(ELim(0), ELim(1), spec, E, res_ch);
        // make ELim the analysis range if not defined
        if(analysisRg.is_empty()) analysisRg = searchRg;
    }

    // no analysis range defined; get it
    if(analysisRg.is_empty())
    {
        analysisRg = lGetAnalysisRange(ECutLow, ECutHigh, spec, E, res_ch);
    }

    if(searchRg.is_empty()) searchRg = analysisRg;
    else {
        if(searchRg(0) < analysisRg(0)) searchRg(0) = analysisRg(0);
        if(searchRg(1) > analysisRg(1)) searchRg(1) = analysisRg(1);
    }

    int ChanLow = searchRg(0);
    int ChanHigh = searchRg(1);

//    qDebug() << "ChanLow: " << ChanLow;
//    qDebug() << "ChanHigh: " << ChanHigh;

    // channels where peak search shall be performed
    rowvec SearchChan = RangeVec2(ChanLow, ChanHigh);
//SearchChan.print("SearchChan");
    // get peak search sensitivity in every search channel (sign so pss > 0 at peaks)
    mat pss = - pksens(spec, res_ch.cols(ChanLow-1, ChanHigh-1), SearchChan, SearchChan);
    pss.reshape(1, pss.size());
    //("peakSearch_pss.txt", pss);

    // look for regions of interest (ROIs)
    int NPeaks = 0;
    int i = 1; // i = 1;
    int NSearch = SearchChan.size();
    int ChanOffset = ChanLow - 1;
    int l, L, r, R;
    //pss.print("pss");
    while(i < NSearch)
    {
        if(pss(i-1) > PSSs)
        {
            // we found a left ROI channel
            l = i;
            L = l + ChanOffset;
            // find the right ROI channel
            while(i < NSearch && pss(i-1) > PSSs)
            {
                i = i + 1;
            }
            r = i-1;
            R = r + ChanOffset;
            rowvec ROIidx = RangeVec2(l, r);

            rowvec ROI = ROIidx + ChanOffset;
            // Centroid ~ Center of Gravity
            // COG = sum( ROI.*pss(ROIidx) ) / sum(pss(ROIidx));
            rowvec temp_pss = pss.cols(l-1, r-1);
            double COG = sum(ROI % temp_pss) / sum(temp_pss);
            double maxPSS = temp_pss.max();
            // accept, if residual peak search, or not too near existing peak
            //... Centroid deviation    max delta at channel floor(Centroid)
            if( !fullSearch || all( abs(COld - COG) > deltaC( qFloor(COG) ) ) )
            {
                // add peak information to list
                CNL.resize(NPeaks+1);
                CNR.resize(NPeaks+1);
                C.resize(NPeaks+1);
                PSSfound.resize(NPeaks+1);

                CNL(NPeaks) = L;
                CNR(NPeaks) = R;
                C(NPeaks) = COG;
                PSSfound(NPeaks) = maxPSS;
                NPeaks = NPeaks + 1;
            }
        }
        i = i + 1;
    }
/* ---------------------------------打印输出-------------------------------- */
    //("peakSearch_CNL.txt", CNL);
    //("peakSearch_CNR.txt", CNR);
    //("peakSearch_C.txt", C);
    //("peakSearch_PSSfound.txt", PSSfound);

    NA = zeros<vec>(C.size());
    // Get Net Area:
    if(doBaseInit)
    {
        // get peak fwhm in channel units
        vec nE = calValues(Cal_Energy, C);
        vec ndE = calDerivative(Cal_Energy, C);
        vec nFE = calValues(Cal_Resolution, nE);
        vec FC = nFE / ndE; // mat FC = nFE ./ ndE;
/* ---------------------------------打印输出-------------------------------- */
        //("peakSearch_FC.txt", FC);

        // strip peaks from spectrum
        rowvec baseCounts, roiLo, roiHi;
        specCutROIs(baseCounts, roiLo, roiHi, spec, C, FC);
        uword j = 0;
        for(uword i = 0, j = 0; i<roiLo.size(); ++i)
        {
            while(j < C.size() && C(j) < roiHi(i))
            {
                // WORK: sumIdx may be wider, but should not exceed halfpoint to next peak
                double yplus = sum( spec.cols(CNL(j)-1, CNR(j)-1) ); // double yplus = sum(spec(sumIdx));
                double ybase = sum( baseCounts.cols(CNL(j)-1, CNR(j)-1) ); // double ybase = sum(baseCounts(sumIdx));
                NA(j) = max(NAMin, yplus - ybase); // NA(j) = max(NAMin, yplus - ybase);
                ++j; //j = j+1;
            }
        }
    } else {
        for(uword i=0; i<C.size(); ++i)
        {
            NA(i) = max(NAMin, sum(missedCounts.cols(CNL(i)-1, CNR(i)-1))); // NA(i) = max(NAMin, temp_sum);
        }
    }
/* ---------------------------------打印输出-------------------------------- */
    //("peakSearch_NA.txt", NA);

    if(!C.is_empty())
    {
        //idx = patAddPeaks(sn, 'Source', Src, 'Centroid', C, 'NetArea', NA, 'Sensitivity', PSSfound, 'CentroidMethod', 'M', 'NetAreaMethod', 's');
        uvec PropIdx;
        PropIdx << i_Centroid << i_NetArea << i_Sensitivity;
        mat Datas(C.size(), 3);
        Datas.col(0) = C;
        Datas.col(1) = NA;
        Datas.col(2) = PSSfound;
        uvec oldIdx;
        patAddPeaks(PropIdx, Datas, pat.Peaks, idx, oldIdx);
    }

    if(doBaseInit)
    {
        rowvec tmp; tmp << 0 << 0;
        tmpbase(analysisRg, tmp, tmp, analysisRg);
        baseInit();
    }
    return idx;
}

mat GammaAnalyALG::calValues(CalibrationType cal, mat Chan)
{
    mat y;
    colvec p;
    if(calGetPATPara(p, cal)) // get calibration parameter for named PAT
    {
        y = Independ::calFcnEval(Chan, p); // evaluate calibration function
    }

    return y;
}

mat GammaAnalyALG::calDerivative(CalibrationType cal, mat Chan)
{
    mat dy = zeros(Chan.n_rows, Chan.n_cols);
    colvec p;
    if(calGetPATPara(p, cal)) // get calibration parameter for named PAT
    {
        dy = Independ::calDerivEval(Chan, p); // evaluate calibration function
    }
    return dy;
}

bool GammaAnalyALG::calGetPATPara(colvec& p, CalibrationType cal)
{
    bool bRet;
    if(cal == Cal_Default) // return default parameter, if name is invalid
    {
        bRet = calDefaultPara(p, cal);
    }
    else { bRet = getCalibrationPara(p, cal); }

    return bRet;
}

bool GammaAnalyALG::getCalParaGlob(colvec &out, CalibrationType cal)
{
    bool bRet = true;
    // get global variable depending on calibration type
    QString curCal = CalPHD;
    switch (cal)
    {
        case Cal_Energy:
            if(Acal_energy_para[m_curEner].size() > 0) out = Acal_energy_para[m_curEner](0);
            break;
        case Cal_Resolution:
            if(Acal_resolution_para[m_curReso].size() > 0) out = Acal_resolution_para[m_curReso](0);
            break;
        case Cal_Efficiency:
            if(Acal_efficiency_para.size() > 0) out = Acal_efficiency_para(0);
            break;
        case Cal_Tot_efficiency:
            if(Acal_tot_efficiency_para.size() > 0) out = Acal_tot_efficiency_para(0);
            break;
        case Cal_Tail: out = Acal_tail_para; break;
        case Cal_Tail_alpha: out = Acal_tail_alpha_para; break;
        case Cal_Tail_right: out = Acal_tail_right_para; break;
        case Cal_Tail_right_alpha: out = Acal_tail_right_alpha_para; break;
        case Cal_Step_ratio: out = Acal_step_ratio_para; break;
        default:
            qDebug() << "invalid calibration string " << cal;
            bRet = false;
    }

    // get list for given spectrum number
    if (out.n_elem < 2)
    {
        bRet = false;
    }

    return bRet;
}

bool GammaAnalyALG::calDefaultPara(colvec& p, CalibrationType cal)
{
    bool bRet = true;

    // assign parameter p depending on calibration string
    switch(cal)
    {
    case Cal_Energy:
        p << 2 << 0 << 0.3 << 0; // linear
        break;
    case Cal_Resolution:
        p << 4 << 1 << 0.003 << 0; // square root of polynomial
        break;
    case Cal_Efficiency:
    case Cal_Tot_efficiency:
        p << 94 << 0.1 << 30 << 3 << 50 << 70 << 0.8; // HAE Efficiency (1-2)
        break;
    case Cal_Tail:
        p << 99 << 0 << 0 << 0; // cutoff linear, 0
        break;
    case Cal_Tail_alpha:
        p << 98 << 1; // constant 1
        break;
    case Cal_Tail_right:
        p << 98 << 0.35; // constant 0.35
        break;
    case Cal_Tail_right_alpha:
        p << 98 << 1; // constant 1
        break;
    case Cal_Step_ratio:
        p << 98 << 0; // constant 0
        break;
    case Cal_HalflifeAnalysis:
        p << 97 << 0 << 1 << log(2); // exponential sum, no constant part, halflife 1 unit
        break;
    default:
        bRet = false;
        qDebug() << "invalid calibration string";
    }

    /* error estimates not reasonable for default calibration
       perr = NaN * ones(1, length(p) - 1); */

    return bRet;
}

bool GammaAnalyALG::getCalibrationPara(colvec& p, CalibrationType cal)
{
    bool bRet = true;

    // get all calibration parameters plus names for this spectrum
    // [s, msg, calentries] = getCalParaGlob(sn, cal);

    bRet = getCalParaGlob(p, cal);
    if(!bRet) bRet = calDefaultPara(p, cal);
    return bRet;
}

rowvec GammaAnalyALG::lGetAnalysisRange(double ECutLow, double ECutHigh, mat spec, mat E, mat res_ch)
{
    // upper cutoff channels
    int CutHighChans = 3;
    // additional cutoff in resolution units
    int CutC = 2;

    int NChan = spec.n_elem;

    // lower cutoff: at least ECutLow, but skip zero's
    int ChanLow = 1;
    while(ChanLow < NChan && spec(ChanLow-1) == 0)
    {
        ChanLow = ChanLow + 1;
    }
    // get additional distance from first non-zero channel
    int t1 = ceil( CutC * res_ch(ChanLow-1) );
    if(t1 > 0) ChanLow = ChanLow + t1;
    // skip more channels if ChanLow is below ECutLow
    if(ChanLow < NChan)
    {
        while( E(ChanLow-1) < ECutLow)
        {
            ChanLow = ChanLow + 1;
            if(E.size()<=ChanLow)
            {
                throw std::string(" Mat::operator(): index out of bounds");
            }
        }
    }

    // ignore certain amount of channels in the end (may be checksum)
    int ChanHigh = NChan - CutHighChans;
    // skip zeros
    while( ChanHigh > ChanLow && spec(ChanHigh-1) == 0)
    {
        ChanHigh = ChanHigh - 1;
    }
    // additional distance for peak search filter
    int t2 = floor( CutC * res_ch(ChanHigh-1) );
    if(t2 > 0) ChanHigh = ChanHigh - t2;

    while(ChanHigh > ChanLow && E(ChanHigh-1) > ECutHigh)
    {
        ChanHigh = ChanHigh - 1;
    }

    if(ChanHigh <= ChanLow)
    {
        qDebug() << "Left Channel bigger than right channel, no peak search!";
        ChanHigh = min(NChan, ChanLow + 1);
        ChanLow = max(1, ChanHigh - 1);
    }
    rowvec rg;
    rg << ChanLow << ChanHigh;
    return rg;
}

colvec GammaAnalyALG::pksens(rowvec spec, mat res_c, mat LROI, mat RROI)
{
    // Cut off filter at 1percent of total area
    // int FilterCutoff = 1; // FilterCutoff=1;

    // check input parameters
    uword nROI = res_c.size(); // nROI = length(res_c);
    if(LROI.size() != nROI || RROI.size() != nROI)
    {
        qDebug() << "2nd to 4th argument must be double vectors of same length";
    }
    res_c.reshape(nROI, 1);
    LROI.reshape(nROI, 1);
    RROI.reshape(nROI, 1);

    uword NSpec = spec.size(); // NSpec=length(spec);
    // myNaN=NaN;
    colvec pss(nROI); // pss=myNaN(ones(size(res_c)));
    pss.fill(gNaN);
    // go through all ROIs
    for(uword i=0; i<nROI; ++i) // for i=1:nROI
    {
        // get filter for that ROI
        rowvec filt;
        uword FilterMean = sdfilt(filt, res_c(i)); // [filt,FilterMean]=sdfilt(res_c(i));
        // find the least sensitivity value in the ROI
        double sens=qInf(); // sens=Inf;
        for(uword ch=LROI(i); ch<=RROI(i); ++ch)
        {
            // find the channels that meet the filter, if centered in channel ch
            urowvec idx = RangeVec(1, filt.size()); // idx=1:length(filt);
            uvec iidx = find((ch+idx-FilterMean>0) && (ch+idx-FilterMean<=NSpec));
            idx = IdxMat(idx, iidx); // idx=idx(iidx);
            // apply the filter
            urowvec idx_temp = idx - 1u;
            urowvec idx_temp_2 = ch+idx_temp-FilterMean;
            rowvec cc = IdxMat(filt, idx_temp) % IdxMat(spec,idx_temp_2); // cc=filt(idx).*spec(ch+idx-FilterMean);
            double dd = sum(cc);
            double sd = sqrt(sum(cc % IdxMat(filt,idx_temp)));
            if(sd>0) // idx not empty
            {
                double ss = dd / sd;
            // lower sensitivity has been found
                if(ss < sens)
                {
                    sens = ss;
                }
            }
        }
        // record pss for this ROI
        pss(i) = sens; // pss(i)=sens;
    }
    return pss;
}

int GammaAnalyALG::sdfilt(rowvec &f, double RC, double CO)
{
    // if resolution value bad, return zero filter
    // if isempty(RC) | ~isfinite(RC) | RC < eps
    uword k = 1; // index of filter mid point
    if(qIsInf(RC) || RC < eps(1.0))
    {
        f << 0; // f = 0;
        // k = 1;
        return k;
    }

    double minPrec = 0.001;
    // catch bad CO values
    if(qIsInf(CO) || CO <= 0) // if isempty(CO) | ~isfinite(CO) | CO <= 0
    {
        CO = 1;
    }
    else if(CO < minPrec)
    {
       // limit precision (and thereby filter width)
       CO = minPrec;
    }

    // transform resolution to sigma
    double f2s = 4 * log(2);        // f2s=4*log(2);
    double c2 = pow(RC, 2) / f2s;   // c2=RC^2/f2s;
    // build one side of filter: tf
    rowvec tf; tf << -100; // tf=-100;
    // k=1;
    double k2 = 1; // k2=1;
    double v = 100*(k2-c2)*exp(-k2/(2*c2))/c2; // v=100*(k2-c2)*exp(-k2/(2*c2))/c2;
    // build untill filter decreases and value is less than cutoff
    while(v > min(tf(k-1u), CO)) // while v>min(tf(k),CO)
    {
        tf.resize(k+1u);
        tf(k) = v;        // tf(k+1)=v;
        ++k;        // k=k+1;
        k2 = pow(k, 2); // k2=k^2;
        v = 100*(k2-c2)*exp(-k2/(2*c2))/c2;
    }
    // make it a 0-sum filter by applying deficiency to entry beside mid-point

    double sum_tf;
    if(k < 3) sum_tf = 0;
    else sum_tf = 50-sum(tf.cols(2, k-1));

    if(k < 2) tf.resize(2);
    tf(1) = sum_tf;    // tf(2)=50-sum(tf(3:k));
    // create the symmetric filter arround mid-point

    f = zeros<rowvec>(2*k-1);             // f=zeros(1,2*k-1);
    for(uword i=0; i<k; ++i)
    {
        f(i+k-1u) = tf(i); // f(k:2*k-1)=tf;
        f(i) = tf(k-1u-i); // f(1:k-1)=tf(k:-1:2);
    }
    return k;
}

void GammaAnalyALG::specCutROIs(rowvec &base, rowvec &roiLo, rowvec &roiHi, rowvec spec, vec C, vec FC)
{
    // number of channels for interpolation
    int NWidth = 3; // NWidth = 3;
    // width of roi in fwhm units
    double fROI = 2.5; // fROI = 2.5;

    uword len = spec.size();
    base = spec; // base = spec;
    // define ROIs
    vec pkLo = max(1, floor(C - fROI*FC) - NWidth); // pkLo = max(1, floor(C - fROI*FC) - NWidth);
    vec pkHi = min(len, ceil(C + fROI*FC) + NWidth); // pkHi = min(length(spec), ceil(C + fROI*FC) + NWidth);
    rowvec roiBool = zeros<rowvec>(len); // roiBool = zeros(size(spec));
    for(uword i=0; i<C.size(); ++i) // for i = 1 : length(C)
    {
        // roiBool(pkLo(i) : pkHi(i)) = 1;
        for(int j=pkLo(i)-1; j<pkHi(i); ++j)
        {
            if(j < 0) continue;
            if(j >= len) roiBool.resize(j+1);
            roiBool(j) = 1;
        }
    }
    // roiBool([1,end+1]) = 0;
    roiBool.resize(roiBool.size()+1);
    roiBool(0) = 0;
    roiBool(roiBool.size()-1) = 0;

    rowvec diff_roiBool = diff(roiBool);
    uvec temp_1 = find(diff_roiBool == 1);
    uvec temp_2 = find(diff_roiBool == -1);
    roiLo.set_size(temp_1.size());
    roiHi.set_size(temp_2.size());
    for(int i=0; i<temp_1.size(); ++i)
    {
        roiLo(i) = temp_1(i);
    }
    for(int i=0; i<temp_2.size(); ++i)
    {
        roiHi(i) = temp_2(i);
    }

    for(uword i=0; i<roiLo.size() && i<roiHi.size(); ++i) // for i = 1 : length(roiLo)
    {
        double yLo = sum(spec.cols(roiLo(i), roiLo(i)+NWidth-1))/NWidth; // yLo = sum(spec(roiLo(i) : roiLo(i)+NWidth-1))/NWidth;
        double yHi = sum(spec.cols(roiHi(i)-NWidth+1, roiHi(i)))/NWidth; // yHi = sum(spec(roiHi(i)-NWidth+1 : roiHi(i)))/NWidth;
        rowvec localSum = cumsum(spec.cols(roiLo(i), roiHi(i))); // localSum = cumsum(spec(roiLo(i) : roiHi(i)));
        // base(roiLo(i) : roiHi(i)) = yLo + (yHi - yLo) * localSum / localSum(end);
        rowvec temp = yLo + (yHi - yLo) * localSum / localSum(localSum.size()-1);
        for(int j=roiLo(i); j<=roiHi(i); ++j)
        {
            base(j) = temp(j-roiLo(i));
        }
    }
}

bool GammaAnalyALG::patAddPeaks(uvec Idx, mat Datas, mat &peaks, uvec &nidx, uvec &oidx)
{
    // 如果索引的个数与数据的个数不对应则返回
    if(Datas.n_cols != Idx.size())
    {
        qDebug() << "The size of Datas and Idx must be same!";
        return false;
    }

    /* 创建一个新的存储峰信息的矩阵，然后把传入的数据 Datas 插入到矩阵对应列（由 Idx 指定）中
     * nold 峰列表中峰的个数， nnew 是新插入峰的个数 */
    int nold = peaks.n_rows;
    int nnew = Datas.n_rows;
    mat newPeaks(nnew, MAX_INDEX_PEAK);
    newPeaks.fill(gNaN);
    for(uword i=0; i<Idx.size(); ++i)
    {
        if(Idx(i) > MAX_INDEX_PEAK) continue;
        newPeaks.col(Idx(i)) = Datas.col(i);
    }

    if(!newPeaks.col(i_Centroid).is_finite())
    {
        qDebug() << "Centroid required";
        return false;
    }
    if(!newPeaks.col(i_FWHM_Ch).is_finite())
    {
		mat mt = newPeaks.col(i_Centroid);
//		PrintMat22("newPeaks.col(i_Centroid).txt", mt);

        colvec E = calValues(Cal_Energy, newPeaks.col(i_Centroid));
        colvec dE = calDerivative(Cal_Energy, newPeaks.col(i_Centroid));
        colvec F = calValues(Cal_Resolution, E);

        newPeaks.col(i_FWHM_Ch) = F / dE; // FC = F ./ dE;

//		PrintColvec("E.txt", E);
//		PrintColvec("dE.txt", dE);
//		PrintColvec("F.txt", F);

		vec pFC = newPeaks.col(i_FWHM_Ch);
//		PrintMat22("newPeaks.col(i_FWHM_Ch).txt", pFC);
    }
    if(!newPeaks.col(i_NetArea).is_finite())
    {
        qDebug() << "NetArea required";
        return false;
    }
    if(!newPeaks.col(i_Sensitivity).is_finite())
    {
        qDebug() << "Sensitivity unknown, setting to NaN";
    }

    peaks = join_vert(peaks, newPeaks);
    uvec index = sortrows(peaks, peaks); // [patjoint, index] = sortrows(patjoint);

    //[sindex, order] = sort(index);
    mat sindex;
    uvec order = sortrows(sindex, UMatToMat(index));
    nidx = order.rows(nold, order.size()-1);// nidx = order(nold + 1 : end);
    if(nold > 0) oidx = order.rows(0, nold-1); // oidx = order(1 : nold);

    colvec sqi, uqi;
    patBaseVar(sqi, uqi, nidx);

    return true;
}

/* =====================================================基线初始化======================================================= */
void GammaAnalyALG::patBaseVar(colvec &sqi, colvec &uqi, uvec idx, int updr)
{
    // patin and patout always assigned, even when nargin == 4
    /*if nargin < 3
       if nargin < 2
          idx = 'all';
       end
       patin = 'Current';
       patout = 'Temporary';
    else
       patin = pat;
       patout = pat;
    end*/

    // setup parameters: back-counts width in fwhm-units, sensitivity factor
    //[k, k_alpha, k_beta] = getSpecSetup(sn, 'k_back', 'k_alpha', 'k_beta');
    double k = specSetup.k_back;
    double k_alpha = specSetup.k_alpha;
    double k_beta = specSetup.k_beta;

    // get baseline and peak parameters
    field<rowvec> t_f1;
    if(updr)
    {
        //[BL, spec, resid] = dmspec(sn, 'BaseLine', 'Spectrum', 'Residual');
        QStringList specItems;
        specItems << "BaseLine" << "Residual";
        t_f1 = dmspec(specItems);
/* -------------------------------打印输出 ---------------------------------- */
        //("patBaseVar_BaseLine.txt", t_f1(0));
        //("patBaseVar_Residual.txt", t_f1(1));
    }
    else {
        t_f1 = dmspec(QStringList("BaseLine")); //BL = dmspec(sn, 'BaseLine');
/* -------------------------------打印输出 ---------------------------------- */
        //("patBaseVar_BaseLine.txt", t_f1(0));
    }

    //[Ct, Fc, NA, Sc] = dmps(sn, patin, 'Centroid', 'FWHM_Ch', 'NetArea', 'SigmaCh');
    QStringList dmpsItems;
    dmpsItems << "Centroid" << "FWHM_Ch" << "NetArea" << "SigmaCh";
    field<vec> t_f2 = dmps(dmpsItems);

    /* // modify all peaks, if idx is string
    if isstr(idx)
       idx = [1 : length(Ct)];
    end*/
    if(t_f2(0).is_empty()) return; // no peaks in PAT or called needlessly
    if(idx.is_empty())
    {
        idx = vectorise( RangeVec(0, t_f2(0).size()-1) );
    }

    vec MBC, LC, LD, BC;
    if(t_f1(0).is_empty()) //if isempty(BL)
    {
        // no baseline, baseline parameters are undefined
        MBC << gNaN;
        //S = gNaN;
        //Err = gNaN;
        LC << gNaN;
        LD << gNaN;
        sqi << gNaN;
        uqi << gNaN;
    }
    else {
        // restrict to indexed peaks
        t_f2(0) = t_f2(0)(idx); //Ct = Ct(idx);
        t_f2(1) = t_f2(1)(idx); //Fc = Fc(idx);
        t_f2(2) = t_f2(2)(idx); //NA = NA(idx);
        t_f2(3) = t_f2(3)(idx); //Sc = Sc(idx);

        vec W = 2 * k * t_f2(1); //W = 2 * k * Fc;
/* -------------------------------打印输出 ---------------------------------- */

        // total background counts from Ct +- k*Fc

        rowvec tmp1 = max(0, t_f1(0));
        BC = Independ::abkcnt( tmp1, t_f2(0), t_f2(1), k); //BC = abkcnt(max(BL, 0), Ct, Fc, k);
/* -------------------------------打印输出 ---------------------------------- */


        // Mean Back Counts
        MBC = BC / W;
/* -------------------------------打印输出 ---------------------------------- */

        /* standard deviation of A is:
         * D(A) = c * f * sqrt(w*b) with
         * w - sigma of gaussian peak (channels)
         * b - mean baseline/channel under peak
         * f - algorithm-dependent factor if baseline known
         * c - baseline estimation correction factor;
         * formerly: D(A)=sqrt(m*2*kb*fc*b), thus
         * f = sqrt(m*2*kb*sqrt(8*log(2))),
         * yields: f=1.9562 for m=0.65, kb=1.25 */
        double c = 1.1;
        double f = 1.9562;
        vec DA = c*f*sqrt(t_f2(3) % MBC); //DA = c*f*sqrt(Sc .* MBC);
/* -------------------------------打印输出 ---------------------------------- */
        //("patBaseVar_DA.txt", DA);

        // Currie's LC
        LC = k_alpha * DA;

        // Currie's LD
        // old formula: LD = k_alpha^2 + 2* LC;
        double kb2 = pow(k_beta, 2);
        // assymmetric formula, ref TBD
        LD = LC + (kb2/2) * ( 1 + sqrt(1 + (4/kb2)*( LC+pow(DA,2) ) ) );
/* -------------------------------打印输出 ---------------------------------- */
        //("patBaseVar_LD.txt", LD);

        if(updr)
        {
            // total counts
            vec STot = Independ::abkcnt(baseInfo(s_Spectrum), t_f2(0), t_f2(1), k); //STot = abkcnt(spec, Ct, Fc, k);
            // signed residual
            vec signedResidual = Independ::abkcnt(t_f1(1), t_f2(0), t_f2(1), k);
            // unsigned residual
            vec unsignedResidual = Independ::abkcnt(abs( t_f1(1) ), t_f2(0), t_f2(1), k);

/* -------------------------------打印输出 ---------------------------------- */
            //("patBaseVar_STot.txt", STot);
            //("patBaseVar_signedResidual.txt", signedResidual);
            //("patBaseVar_unsignedResidual.txt", unsignedResidual);

            vec sqrtCT = sqrt(STot);
            sqi = signedResidual / sqrtCT;
            uqi = unsignedResidual / sqrtCT;
        }
    }

    // write results to temporary PAT
    if(updr) lSetPat(idx, MBC, BC, LC, LD, sqi, uqi);
    else lSetPat(idx, MBC, BC, LC, LD);
}

void GammaAnalyALG::lSetPat(uvec idx, vec mbc, vec cba, vec lc, vec ld, vec sqi, vec uqi)
{
    /*% get PAT table
    [s, msg, patEntry] = getPATByName(sn, patin);
    if s
       pd = patEntry{2};
       % write PAT entries, compare helpApat_list */
       //MatAssign(pat.Peaks.col(i_MeanBackCount), idx, mbc); //pd(idx, 14) = mbc;
       //MatAssign(pat.Peaks.col(i_LC), idx, lc); //pd(idx, 34) = lc;
       //MatAssign(pat.Peaks.col(i_LD), idx, ld); //pd(idx, 35) = ld;
       //MatAssign(pat.Peaks.col(i_BackgroundArea), idx, cba); //pd(idx, 41) = cba;

       for(uword i=0; i<idx.size(); ++i)
       {
           pat.Peaks.col(i_MeanBackCount)(idx(i)) = mbc(i);
           pat.Peaks.col(i_LC)(idx(i)) = lc(i);
           pat.Peaks.col(i_LD)(idx(i)) = ld(i);
           pat.Peaks.col(i_BackgroundArea)(idx(i)) = cba(i);
           pat.Peaks.col(i_Significance)(idx(i)) = gNaN;
       }

       // either update residuals, or raise flag that they are not up to date
       if(!sqi.is_empty() && !uqi.is_empty()) //if nargin == 10
       {
          //MatAssign(pat.Peaks.col(i_SignedResidual), idx, sqi); //pd(idx, 36) = sqi;
          //MatAssign(pat.Peaks.col(i_UnsignedResidual), idx, uqi); //pd(idx, 37) = uqi;
          //MatAssign(pat.Peaks.col(i_ResidualUpToDate), idx, 1); //pd(idx, 38) = 1;
          for(uword i=0; i<idx.size(); ++i)
          {
              pat.Peaks.col(i_SignedResidual)(idx(i)) = sqi(i);
              pat.Peaks.col(i_UnsignedResidual)(idx(i)) = uqi(i);
              pat.Peaks.col(i_ResidualUpToDate)(idx(i)) = 1;
          }
       }
       else {
          //MatAssign(pat.Peaks.col(i_ResidualUpToDate), idx, 0); //pd(idx, 38) = 0;
          for(uword i=0; i<idx.size(); ++i)
          {
              pat.Peaks.col(i_ResidualUpToDate)(idx(i)) = 0;
          }
       }

       // significance overwritten
       //MatAssign(pat.Peaks.col(i_Significance), idx, gNaN); //pd(idx, 16) = NaN;

       /*patEntry{1} = patout;
       patEntry{2} = pd;
       [s, msg] = setPATEntry(sn, patEntry);
       if ~s
          logWrite('Failed writing Baseline PAT variables: %s', msg);
       end
    else
       logWrite('Problem while updating Baseline PAT variables: %s', msg);
    end*/
}

void GammaAnalyALG::baseInit(double PSSBreak, double PSSfwhm)
{
    rowvec cx, cy, cdy;
    //% least number of data points to be fitted for 1 CP
    double NMin = 5;
    //% absolute minimum distance of CPs
    double NCritical = 10;
    //% minimum number of data points between ROIs to use 3 CPs
    double NTwo = 50; //% should come from data condition

    rowvec data;
    rowvec step;
    rowvec rg;
    rowvec roiMask;
    uvec L;
    uvec H;
	stdvec vy;
	vy = getBaseInfo(s_YControl);
    lGetData(data, step, rg, roiMask, L, H, PSSBreak, PSSfwhm);
//	PrintMat22("baseInit_data.txt", data);
//	PrintMat22("baseInit_step.txt", step);
//	PrintMat22("baseInit_rg.txt", rg);
//	PrintMat22("baseInit_roiMask.txt", roiMask);
	vy = getBaseInfo(s_YControl);
    lInitBase(data, step, rg, roiMask, L, H, NMin, NCritical, NTwo, cx, cy, cdy);
//	PrintMat22("baseInit_cx.txt", cx);
//	PrintMat22("baseInit_cy.txt", cy);
//	PrintMat22("baseInit_cdy.txt", cdy);
	vy = getBaseInfo(s_YControl);
    tmpbase(cx,cy,cdy);
	vy = getBaseInfo(s_YControl);
    return;
}

void GammaAnalyALG::tmpbase(mat cx, mat cy, mat cdy, mat ar)
{
    //error(nargchk(4,5,nargin));
    /*if(ar.is_empty()) //if nargin < 5
    {
        if(!cx.is_empty()) // isempty(cx)
        {
            ar = baseInfo(s_AnalysisRange);
        }
    }*/
    if(!ar.is_empty())
    {
        if(ar.is_colvec()) baseInfo(s_AnalysisRange) = ar.t();
        baseInfo(s_AnalysisRange) = ar;
    }

    // create baseline structure
    /*global Ac_baseline_tmp
    Ac_baseline_tmp=struct('AnalysisRange', {ar}, ...
        'XControl', {cx(:)'}, ...
        'YControl', {cy(:)'}, ...
        'YSlope', {cdy(:)'}, ...
        'SpectrumNumber', {sn});*/
    if(cx.is_colvec()) baseInfo(s_XControl) = cx.t();
    else baseInfo(s_XControl) = cx;
    if(cy.is_colvec()) baseInfo(s_YControl) = cy.t();
    else baseInfo(s_YControl) = cy;
    if(cdy.is_colvec()) baseInfo(s_YSlope) = cdy.t();
    else baseInfo(s_YSlope) = cdy;

    // update baseline dependent variables
    vec sqi, uqi;
    patBaseVar(sqi, uqi);
}

void GammaAnalyALG::findROI(urowvec& l, urowvec& h, rowvec& roiMask,
             vec fLo, uvec idx, vec fHi)
{
    if(fHi.is_empty()) fHi = fLo;

    QStringList Opts; Opts << "Centroid" << "FWHM_Ch";
    field<vec> t_dmps = dmps(Opts);
    vec& C = t_dmps(0);
    vec& FWHM_CH = t_dmps(1);

    rowvec rg = baseInfo(s_AnalysisRange);

    if(idx.is_empty()) idx = RangeVec(0, C.size()-1).t();

    vec temp_FWHM = FWHM_CH(idx);
    vec temp_C = C(idx);
    vec tempL = fLo % temp_FWHM;
    vec tempH = fHi % temp_FWHM;

    vec pkLo = max(rg(0), floor(temp_C - tempL));
    vec pkHi = min(rg(1), ceil(temp_C + tempH));
    //("findROI_pkLo.txt", pkLo);
    //("findROI_pkHi.txt", pkHi);

    roiMask = zeros<rowvec>(m_nChans);

    for(size_t i=0; i!=pkLo.size(); ++i)
    {
         roiMask(span(pkLo(i)-1, pkHi(i)-1)) = ones<rowvec>(pkHi(i)-pkLo(i)+1);
    }

    roiMask[0] = 0;
    roiMask[roiMask.size()-1] = 0;

    rowvec db = diff(roiMask);
    l = find(db == 1).t() + 1;
    h = find(db == -1).t();
}

void GammaAnalyALG::lGetData(rowvec& data, rowvec& step, rowvec& rg, rowvec& roiMask, uvec& L, uvec& H,
              double PSSBreak, double PSSfwhm)
{
    double maxFWHMfact = 5;

    // same as matlab [stripped, step, rg] = dmspec(sn, 'Stripped', 'Steps', 'AnalysisRange');
    QStringList Opts;
    Opts << "Stripped" << "Steps" << "AnalysisRange";
    field<rowvec> t_f = dmspec(Opts);
    rowvec& stripped = t_f(0);
    step = t_f(1);
//	PrintMat22("stripped.txt", stripped);
//	PrintMat22("step.txt", step);
    rg = t_f(2);
    if(rg.is_empty())
    {
        rg = RangeVec2(1, m_nChans);
    }

    // same as matlab [pC, pPSS, pFC] = dmps(sn, 'Centroid', 'Sensitivity', 'FWHM_Ch');
    Opts.clear();
    Opts << "Centroid" << "Sensitivity";
    field<vec> f_dmps = dmps(Opts);
    vec& pC = f_dmps(0);
    vec& pPSS = f_dmps(1);
    //vec pFC = pat.Peaks.col(i_FWHM_Ch);

    data = stripped - step;

    // get peak ranges
    colvec factFWHM = 1 + (maxFWHMfact -1)*pPSS/PSSfwhm;
    factFWHM = min(maxFWHMfact,factFWHM); // if factFWHM's members bigger then maxFWHMfact set them equal maxFWHMfact
    factFWHM = NaNto1(factFWHM); //factFWHM(isnan(factFWHM)) = 1;
    //("lGetData_factFWHM.txt", factFWHM);

    urowvec l;
    urowvec h;

    findROI(l, h, roiMask, factFWHM); // l,h,roiMask are out params
    l = l + 1u;
    h = h + 1u;
    //("lGetData_roiMask.txt", roiMask);

    // get borders of ranges containing big peaks
    urowvec isbig = zeros<urowvec>(l.size()); // isbig = false(size(l));
    uvec idx;
    for(size_t i=0; i<l.size(); ++i) // for i = 1 : length(l)
    {
        idx = find(pC > l(i) && pC < h(i)); // idx = find(pC > l(i) & pC < h(i));
        if(any(pPSS(idx) > PSSBreak)) isbig(i) = 1; //if any(pPSS(idx) > PSSBreak), isbig(i) = true; end
    }

    idx = find(isbig > 0);
    L = l(idx); // L = l(isbig);
    H = h(idx); // H = h(isbig);

    // remove ROIs that end to near lower cutoff
    idx = find(H <= (rg(0) + 2)); // idx = find(H <= (rg(1)+2));
    if(!idx.is_empty()) // if any(idx), L(idx)=[]; H(idx)=[]; end
    {
        for(int i=idx.size()-1; i>=0; --i)
        {
            L.shed_col(idx(i));
            H.shed_col(idx(i));
        }
    }

    // shift ROI limit if first ROI goes beyond lower cutoff
    if((!L.is_empty()) && (L[0] <= rg[0]))
    {
        L[0] = rg[0] + 1;
    }

    // add lower cutoff as a ROI
    uvec tmp;
    tmp << 1u;
    L = join_vert(tmp, L); // L = [1; L(:)];
    tmp << rg(0);
    H = join_vert(tmp, H); // H = [rg(1); H(:)];
}

mat GammaAnalyALG::restpolyfit(vec& r, vec x, vec y, vec p0, vec w)
{
    if( x.size() != y.size() )
    {
        qDebug() << "X and Y vectors must be the same size."; // messagebox
    }
    if(w.is_empty())
    {
        w = ones<vec>(x.size());
    }
    else {
        y = y % w;
    }

    int n = p0.size()-1;

    // Construct 'weighted' Vandermonde matrix.
    mat V = zeros<mat>(w.size(),n+1);
    V.col(n) = w;
    for(int i=n-1; i>=0; --i)
    {
        V.col(i) = x % V.col(i+1);
    }

    // get free and fixed component
    uvec freeIdx = find_nonfinite(p0);
    uvec fixIdx = find_finite(p0);

    mat p;
    vec yfix = MatCols(V, fixIdx) * p0(fixIdx);//yfix = V(:, fixIdx)*p0(fixIdx);
    vec yfree = y - yfix;

    mat Q, R;
    qr_econ(Q, R, MatCols(V, freeIdx));
    p = solve(R, trans(Q) * yfree); // 待测试，P应为行向量

    colvec yhat = yfix + MatCols(V, freeIdx) * p;
    r = y - yhat;
    p0(freeIdx) = p;
    p = p0.t();

    return p;
}

void GammaAnalyALG::lNorm(rowvec yy, rowvec mask, double xc, double xh, double yh, double np,
      double& nrm, double& lbda, double& yc, double& dyc, double& xl, double& yho, double& dyh)
{
    // check argument, get polynomial mask
    mat p;
    if(np == 4)
    {
        p << gNaN << gNaN << gNaN << yh;
    }
    else if(np == 2)
    {
        p << gNaN << yh;
    }
    else qDebug() << "invalid np";

    // get whole fitting interval
    double k = xh - xc;
    xl = max(1, xh-2*k); //xl = max(1, xh-2*k);
    rowvec xx = RangeVec2(xl-1,xh-1); //xx = [xl : xh];

    // remove flagged points
    uvec idx = find( IdxMat(mask,xx) == true );
    for(int i=idx.size()-1; i>=0; --i)
    {
        xx.shed_col(idx(i));
    }
    if(xx.size()<5)
    {
        nrm = 0;
        lbda = 1;
        yc = yy(xc);
        dyc = 0;
        yh = yy(xh);
        dyh = 0;
        return;
    }

    // limit to actual points
    rowvec y = IdxMat(yy,xx);
    rowvec x = xx - xh + 1;

    // perform fit
    colvec r;
    p = restpolyfit(r, vectorise(x), vectorise(y), vectorise(p));

    // polyder.m
    rowvec dp, t_null;
    Matlab::polyder(dp, t_null, 1, p);

    nrm = accu(r % r);
    lbda = max(1, mean(sqrt(abs(y))));

    rowvec rg; rg << -k << 0;
    rowvec yval = Matlab::polyval(p, rg);
    rowvec dyval = Matlab::polyval(dp, rg);

    yc = yval(0);
    yho = yval(1);
    dyc = dyval(0);
    dyh = dyval(1);
}

void GammaAnalyALG::lSelectNext(double x, double y, rowvec data, double xb, rowvec imask, double np,
                  double& xn, double& yn, double& dyn, bool& hit, double& y1, double& dy1)
{
    // selection algorithm parameters:
    // lambda multiplier for norm threshold
    double NMax = 4000;

    // minimum interval length
    double MinDist = 10;

    // get next step
    double hi = max(1, x - MinDist);
    xn = max(1, hi - 1);

    // binary search
    // get norms for initial new x
    double norm, lambda, x0, y0, dy0;
    //[norm, lambda, yn, dyn, x0, y0, dy0] = lNorm(data, imask, xn, x, y,  np);
    lNorm(data, imask, xn, x, y, np, norm, lambda, yn, dyn, x0, y0, dy0);

    while(hi != xn)
    {
        if(norm < NMax * lambda)
        {
            double tmp = xn;
            xn = max(1, 3*xn - 2*hi);
            hi = tmp;
        }
        else
        {
            xn = hi - floor((hi-xn)/2);
        }

        lNorm(data, imask, xn, x, y, np, norm, lambda, yn, dyn, x0, y1, dy1);
    }

    if(xn < xb)
    {
        xn = xb;
        lNorm(data, imask, xn, x, y, np, norm, lambda, yn, dyn, x0, y1, dy1);
    }

    hit = (x0 <= xb);
}

void GammaAnalyALG::lFitRange(double x0, double y0, double x1, rowvec xv, rowvec data, rowvec mask, int npara,
               rowvec& yv, rowvec& dyv)
{
    rowvec xx = RangeVec2(min(x0, x1), max(x0,x1));
    //xx(mask(xx)==true) = [];
    rowvec t_xx = xx - 1;
    uvec idx = find( IdxMat(mask,t_xx) == true );
    for(int i=idx.size()-1; i>=0; --i)
    {
        if(idx(i) < xx.size())
        {
            xx.shed_col(idx(i));
            t_xx.shed_col(idx(i));
        }
    }

    if(xx.size()<5)
    {
        qDebug() << QString("Warning: only %1 points between %2 and %3")
                    .arg(xx.size()).arg(min(x0,x1)).arg(max(x0,x1));
        rowvec t_xv = xv - 1;
        yv = IdxMat(data, t_xv);
        dyv = zeros( size(xv) );
        return;
    }

    // fit linear or cubic polynomial
    rowvec p0;
    switch(npara)
    {
        case 2: p0 << gNaN << y0; break;
        case 4: p0 << gNaN << gNaN << gNaN << y0; break;
        default: qDebug() << "Bug: invalid npara value"; break;
    }

    // fit polynomial
    rowvec yy = IdxMat(data, t_xx);
    xx = xx - x0;
    colvec r;
    mat p = restpolyfit(r, vectorise(xx), vectorise(yy), vectorise(p0));

    // evaluate
    yv = Matlab::polyval(p, xv - x0);

    mat a, b;
    Matlab::polyder(a, b, 1, p);
    dyv = Matlab::polyval(a, xv - x0);
}

void GammaAnalyALG::lInitBase(const rowvec& data, const rowvec& step, const rowvec& rg,
               const rowvec& roiMask, const uvec& L, const uvec& H,
               double N1, double N2, double N3,
               rowvec& cx, rowvec& cy, rowvec& cdy)
{
    size_t roiIdx = L.size();
    double xb = H[roiIdx-1];

    cx.set_size(0);
    cy.set_size(0);
    cdy.set_size(0);

    rowvec yn, dyn, yo, dyo, tmp;
    int nP = 0;
    bool pingold = true, pingNew = false;
    double x = rg[1], xn, y = gNaN, dy = gNaN;
    double t_yn, t_dyn, t_yo, t_dyo;
    while(x > rg[0])
    {
        // get next controlpoints
        if(pingold) nP = 2;
        else nP = 4;

        //[xn, yn, dyn, pingNew, yo, dyo] = lSelectNext(x,y,data,xb,roiMask,nP);
        lSelectNext(x, y, data, xb, roiMask, nP, xn, t_yn, t_dyn, pingNew, t_yo, t_dyo);
        yn << t_yn; dyn << t_dyn; yo << t_yo; dyo = t_dyo;

        if(pingold)
        {
            // starting left of a ROI
            if(pingNew)
            {
                // delta = x - xb + 1;
                double delta = accu( roiMask.cols(xb,x) == false );

                if(delta < N2)
                {
                    if(cx.is_empty())
                    {
                        // highest controlpoint, must be introduced
                        //[cy, cdy] = lFitRange(x,NaN,xb, x, data, roiMask, 2);
                        tmp << x;
                        lFitRange(x, gNaN, xb, tmp, data, roiMask, 2, cy, cdy);
                        cx << x;
                    }
                    else if(roiIdx == 1)
                    {
                        // lowest controlpoint, must be introduced
                        //[yn, dyn] = lFitRange(xb,NaN,x, xb, data, roiMask, 2);
                        tmp << xb;
                        lFitRange(xb, gNaN, x, tmp, data, roiMask, 2, yn, dyn);
                        cx = join_horiz(tmp, cx); //cx = [xb, cx];
                        cy = join_horiz(yn, cy); //cy = [yn, cy];
                        cdy = join_horiz(dyn, cdy); //cdy = [dyn, cdy];
                    }
                    else
                    {
                        if(delta < N1)
                        {
                            // no control point introduced
                        }
                        else
                        {
                            // introducing only one control point
                            double xc = round((x + xb)/2);
                            tmp << xc;
                            lFitRange(xb,gNaN,x, tmp, data, roiMask, 2, yn, dyn);
                            cx = join_horiz(tmp, cx); //cx = [xc, cx];
                            cy = join_horiz(yn, cy); //cy = [yn, cy];
                            cdy = join_horiz(dyn, cdy); //cdy = [dyn, cdy];
                        }
                    }
                }
                else if(delta < N3)
                {
                    //[yn, dyn] = lFitRange(x, y, xb, [xb, x], data, roiMask, 2);
                    tmp << xb << x;
                    lFitRange(x, y, xb, tmp, data, roiMask, 2, yn, dyn);
                    cx = join_horiz(tmp, cx); //cx = [xb, x, cx];
                    cy = join_horiz(yn, cy); //cy = [yn, cy];
                    tmp << gNaN << gNaN;
                    cdy = join_horiz(tmp, cdy); //cdy = [NaN, NaN, cdy];
                }
                else
                {
                    // introducing three control points
                    double xc = round((x + xb)/2);
                    //[yn, dyn] = lFitRange(xb,NaN,x,[xb,xc,x],data,roiMask,4);
                    tmp << xb << xc << x;
                    lFitRange(xb, gNaN, x, tmp, data, roiMask, 4, yn, dyn);
                    cx = join_horiz(tmp, cx); //cx = [xb,xc,x,cx];
                    cy = join_horiz(yn, cy); //cy = [yn, cy];
                    tmp << dyn(0) << gNaN << gNaN;
                    cdy = join_horiz(tmp, cdy); //cdy = [dyn(1), NaN, NaN, cdy];
                }

                if(roiIdx == 1)
                {
                    break;
                }

                x = L(roiIdx -1);
                y = gNaN;
                dy = gNaN;
                roiIdx = roiIdx -1;
                xb = H(roiIdx -1);
            }
            else {
                // normal algorithm left of ROI
                tmp << x;
                cx = join_horiz(tmp, cx); //cx = [x, cx];
                cy = join_horiz(yo, cy); //cy = [yo, cy];
                cdy = join_horiz(dyo, cdy); //cdy = [dyo, cdy];
                // continue from here
                x = xn;
                y = t_yn;
                dy = t_dyn;
            }
        }
        else
        {
            // starting from normal CP
            if(pingNew)
            {
                // hitting a ROI
                if(roiIdx == 1)
                {
                    // hitting lower analysis range
                    if(xn == xb)
                    {
                        // next control point at end
                        //[yn, dyn] = lFitRange(x, y, xb, xb, data,roiMask, 2);
                        tmp << xb;
                        lFitRange(x, y, xb, tmp, data, roiMask, 2, yn, dyn);
                        cx = join_horiz(tmp, cx); //cx = [xb, cx];
                        cy = join_horiz(yn, cy); //cy = [yn, cy];
                        cdy = join_horiz(dyn, cdy); //cdy = [dyn, cdy];
                    }
                    else
                    {
                        double xm = round((x+xb)/2);
                        tmp << xb << xm;
                        lFitRange(x, y, xb, tmp, data, roiMask, 4, yn, dyn);

                        cx = join_horiz(tmp, cx); //cx = [xb, xm, cx];
                        cy = join_horiz(yn, cy); //cy = [yn, cy];
                        tmp << dyn(0) << gNaN;
                        cdy = join_horiz(tmp, cdy); //cdy = [dyn(1), NaN, cdy];
                    }
                    break;
                }
                else if(xn == xb)
                {
                    // hitting ROI completely
                    tmp << xn; cx = join_horiz(tmp, cx); //cx = [xn, cx];
                    cy = join_horiz(yn, cy); //cy = [yn, cy];
                    cdy = join_horiz(dyn, cdy); //cdy = [dyn, cdy];
                }
                else {
                    double xc = ceil((x+xb)/2);
                    tmp << xb << xc;
                    lFitRange(x, y, xb, tmp, data, roiMask, 4, yn, dyn);
                    //add righthand CP, with fixed dy;
                    cx = join_horiz(tmp, cx); //cx = [xb, xc, cx];
                    cy = join_horiz(yn, cy); //cy = [yn, cy];
                    tmp << dyn(0) << gNaN;
                    cdy = join_horiz(tmp, cdy); //cdy = [dyn(1), NaN, cdy];
                }
                // continue on lefthand, starting with free slope
                x = L(roiIdx-1);
                y = gNaN;
                dy = gNaN;
                roiIdx = roiIdx - 1;
                xb = H(roiIdx-1);
            }
            else
            {
                // normal procedure
                tmp << xn; cx = join_horiz(tmp, cx); //cx = [xn, cx];
                cy = join_horiz(yn, cy); //cy = [yn, cy];
                tmp << gNaN; cdy = join_horiz(tmp, cdy); //cdy = [NaN, cdy];
                // continue from here
                x = xn;
                y = t_yn;
                dy = t_dyn;
            }
        }
        pingold = pingNew;
    }

    if( qIsNaN( cdy(0) ) )
    {
        cdy(0) = 0;
    }
    if( qIsNaN( cdy( cdy.size()-1 ) ) )
    {
        cdy( cdy.size()-1 ) = 0;
    }
    if(cx.size() == 1)
    {
        cx = rg;
        double tmp = cy(0);
        cy << tmp << tmp;
        cdy << 0 << 0;
        qDebug() << "Whole baseline was fitted at once";
    }
    rowvec t_cx = cx - 1;
    cy = cy + IdxMat(step, t_cx);
}

/* =====================================================基线优化模块======================================================== */
void GammaAnalyALG::baseImprove(double pss_r)
{
    // accept old temporary results to make sure they are not lost
    //patCommit(sn);
    vec sqi, uqi;
    patBaseVar(sqi, uqi); //tbaseac(sn);

    // search peaks in residual, written to temporary PAT
    uvec idx = PeakSearch(pss_r, search_residual);
    //idx.print("idx = ");

    // fit net areas of new peaks
    if(!idx.is_empty())
    {
        uvec ep;
        fitPeakFull(idx, ep, ep, ep);
    }

    // fit the baseline to the residual
    FitBase();

    // undo insertion of residual peaks
    //patRollback;
    int size_idx = idx.size();
    for(int i=size_idx-1; i>=0; --i)
    {
        pat.Peaks.shed_row(idx(i));
    }
}

void GammaAnalyALG::fitPeakFull(uvec Af, uvec Cf, uvec Ff, uvec SRf, uvec CPf, bool bVerbose)
{
    QString addArg; //addArg = {};
    if(bVerbose) addArg = "Peaks above spline baseline";

    // get index to all peaks where any parameter is free
    // get common min and max
    uvec allidx = join_vert( join_vert(Af, Cf), join_vert(Ff, SRf) );
    if(allidx.is_empty()) return;
    //uword LIdx = min(allidx);
    uword RIdx = max(allidx);
    // get unique allidx
    vec Bool = zeros<vec>(RIdx+1u);
    MatAssign(Bool, allidx, 1);
    allidx = find(Bool);

    /*% get peak and baseline parameters, and spectrum data
    [C, FC, NA, Sr, T, Ta, UT, UTa, BW, RB, RD] = dmps(sn, ...
     'Centroid', 'FWHM_Ch', 'NetArea', 'StepRatio', 'Tail', 'TailAlpha', ...
     'UpperTail', 'UpperTailAlpha', 'BWWidthChan', ...
     'RecoilBetaChan', 'RecoilDeltaChan');*/
    QStringList Opt;
    Opt << "Centroid" << "FWHM_Ch" << "NetArea" << "StepRatio" << "Tail" << "TailAlpha"
        << "UpperTail" << "UpperTailAlpha" << "BWWidthChan" << "RecoilBetaChan" << "RecoilDeltaChan";
    field<vec> t_p = dmps(Opt);
    vec& C = t_p(0);
    vec& FC = t_p(1);
    vec& NA = t_p(2);
    vec& Sr = t_p(3);
    vec& T = t_p(4);
    vec& Ta = t_p(5);
    vec& UT = t_p(6);
    vec& UTa = t_p(7);
    vec& BW = t_p(8);
    vec& RB = t_p(9);
    vec& RD = t_p(10);

    //[S, rg, cx, cy, cdy] = dmspec(sn, 'ShortSpectrum', 'AnalysisRange', 'XControl', 'YControl', 'YSlope');
    Opt.clear();
    Opt << "ShortSpectrum" << "AnalysisRange" << "XControl" << "YControl" << "YSlope";
    field<rowvec> t_s = dmspec(Opt);
    rowvec& S = t_s(0);
    rowvec& rg = t_s(1);
    rowvec& cx = t_s(2);
    rowvec cy = t_s(3);
    rowvec cdy = t_s(4);

    // get fitting region
    rowvec x;
    if( rg.size() != 2 ) //if length(rg) ~= 2
    {
        qDebug() << "AnalysisRange not an interval, fitting peaks might result in error";
        x = RangeVec2(1, S.size()); // x = [1 : length(S)];
    }
    else
    {
        x = RangeVec2(rg(0), rg(1)); // x = [rg(1) : rg(2)];
    }
    rowvec t_x = x - 1;
    rowvec y = IdxMat(S, t_x);

    /* % fit parameters
    [s, CXo, cy, cdy, Ao, Co, Fo, Sro, To, Tao, UTo, UTao, BWo, RBo, RDo, X2] = ...
       fitFunction(x, y, 'wrapStepRatio', ...
       cx, [], cy, CPf, cdy, [], NA, Af, C, Cf, FC, Ff, Sr, SRf, ...
       T, [], Ta, [], UT, [], UTa, [], BW, [], RB, [], RD, [], ...
       'Restrict', addArg{:});*/
    double X2;
    uvec ep;
    field<vec> Ps;
    field<uvec> free_idx;
    Ps << vectorise(cx) << vectorise(cy) << vectorise(cdy) << NA << C << FC << Sr << T << Ta << UT << UTa << BW << RB << RD;
    free_idx << ep << CPf << ep << Af << Cf << Ff << SRf << ep << ep << ep << ep << ep << ep << ep;
    bool s = Independ::fitFunction(Ps, X2, free_idx, vectorise(x), vectorise(y), "wrapStepRatio", QStringList("Restrict"), addArg);
    //mat& CXo = Ps(0);
    cy = Ps(1).t();
    cdy = Ps(2).t();
    mat& Ao = Ps(3);
    mat& Co = Ps(4);
    mat& Fo = Ps(5);
    mat& Sro = Ps(6);
    //mat& To = Ps(7);
    //mat& Tao = Ps(8);
    //mat& UTo = Ps(9);
    //mat& UTao = Ps(10);
    //mat& BWo = Ps(11);
    //mat& RBo = Ps(12);
    //mat& RDo = Ps(13);

    // only write to temp pat if fitting converged
    if(s)
    {
        tmpbase(cx, cy, cdy); //tmpbase(sn, cx, cy, cdy);
        //("cy.txt", cy);
        //("cdy.txt", cdy);
        /*patSetPeaks(allidx, 'FittingMethod', 'F', 'NetArea', Ao(allidx), ...
          'Centroid', Co(allidx), 'FWHM_Ch', Fo(allidx), ...
          'NetAreaFree', Af, 'FWHMFree', Ff, 'CentroidFree', Cf);*/
        QStringList Descs;
        Descs << "NetArea" << "Centroid" << "FWHM_Ch";// << "NetAreaFree" << "FWHMFree" << "CentroidFree";
        field<vec> Vals;
        Vals << vectorise(Ao(allidx)) << vectorise(Co(allidx)) << vectorise(Fo(allidx));
             //<< UMatToMat(Af) << UMatToMat(Ff) << UMatToMat(Cf);
        patSetPeaks(allidx, Descs, Vals);

        Vals << vectorise(Sro(SRf));
        patSetPeaks(SRf, QStringList("StepRatio"), Vals);

        vec sqi, uqi;
        patBaseVar(sqi, uqi, allidx);
    }

    // return only the values of free parameters, if any are requested
    //if nargin > 1
    //Ao = Ao(Af);
    //Co = Co(Cf);
    //Fo = Fo(Ff);
    //Sro = Sro(SRf);
    //end
}

void GammaAnalyALG::patSetPeaks(uvec idx, QStringList Descs, field<vec> Values, bool bAll)
{
    // get spectrum number and peak indizes, expand idx 'all' to indizes for all peaks
    if(idx.is_empty() && bAll) //if isstr(idx) & strcmp(idx, 'all')
    {
        idx = RangeVec(0, pat.Peaks.n_rows); //idx = dmps(sn, 'Index');
    }

    // list of possible descriptors and their number in the list
    // list must be unique, variable numbers must match
    /*QString Descts[40] = {"FittingMethod", "Centroid", "FWHM_Ch", "NetArea", "AreaError", "StepRatio", "Sensitivity",
                         "CentroidMethod", "NetAreaMethod", "FWHMMethod", "NuclideSoft", "CentroidFree", "NetAreaFree",
                         "FWHMFree", "Spurious", "Reviewed", "MeanBackCount", "Significance", "ManualTail", "TailChanged",
                         "ResChanged", "BWWidth", "LineEnergy", "LineEnergyErr", "Yield", "YieldErr", "ReferenceID",
                         "SumPeakArea", "SumPeakAreaErrP", "NuclideHard", "Nuclide", "BWWidthErr", "CCF", "CCFErr",
                         "RefEfficiencyErr", "TailAlpha", "UpperTail", "UpperTailAlpha", "RecoilBeta", "RecoilDeltaE"};*/
    //int nold = pat.Peaks.n_rows;
    // parse descriptor-value pairs, write values to numbered positions in Value cell
    uword i=0;
    foreach(QString item, Descs) //for i= 3 : 2 : nargin-1
    {
        if(item == "Centroid")
        {
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_Centroid)(idx(ii)) = Values(i)(ii);
            }
        }
        else if(item == "FWHM_Ch")
        {
            vec FC = Values(i);
            if(any(FC < 0))
            {
                FC = abs(FC);
                qDebug() << "Negative FWHM forbidden, using absolute value";
            }
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_FWHM_Ch)(idx(ii)) = FC(ii);
            }
        }
        else if(item == "StepRatio")
        {
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_StepRatio)(idx(ii)) = Values(i)(ii);
            }
        }
        else if(item == "NetArea")
        {
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_NetArea)(idx(ii)) = Values(i)(ii);
            }
            /* reset area error, if not specified
            if isempty(Value{AreaError_n})
                pat(idx, 13) = NaN;
            end*/
        }
        else if(item == "AreaError")
        {
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_AreaError)(idx(ii)) = Values(i)(ii);
            }
        }
        else if(item == "MeanBackCount")
        {
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_MeanBackCount)(idx(ii)) = Values(i)(ii);
            }
        }
        else if(item == "Sensitivity")
        {
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_Sensitivity)(idx(ii)) = Values(i)(ii);
            }
        }
        else if(item == "Significance")
        {
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_Significance)(idx(ii)) = Values(i)(ii);
            }
        }
        else if(item == "BWWidth")
        {
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_BWWidth)(idx(ii)) = Values(i)(ii);
            }
        }
        else if(item == "BWWidthErr")
        {
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_BWWidthErr)(idx(ii)) = Values(i)(ii);
            }
        }
        else if(item == "CCF")
        {
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_CCF)(idx(ii)) = Values(i)(ii);
            }
        }
        else if(item == "CCFErr")
        {
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_CCFErr)(idx(ii)) = Values(i)(ii);
            }
        }
        else if(item == "LineEnergy")
        {
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_LineEnergy)(idx(ii)) = Values(i)(ii);
            }
        }
        else if(item == "LineEnergyErr")
        {
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_LineEnergyErr)(idx(ii)) = Values(i)(ii);
            }
        }
        else if(item == "Yield")
        {
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_Yield)(idx(ii)) = Values(i)(ii);
            }
        }
        else if(item == "YieldErr")
        {
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_YieldErr)(idx(ii)) = Values(i)(ii);
            }
        }
        else if(item == "ReferenceID")
        {
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_ReferenceID)(idx(ii)) = Values(i)(ii);
            }
        }
        else if(item == "SumPeakArea")
        {
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_SumPeakArea)(idx(ii)) = Values(i)(ii);
            }
        }
        else if(item == "SumPeakAreaErrP")
        {
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_SumPeakAreaErr)(idx(ii)) = Values(i)(ii);
            }
        }
        else if(item == "RefEfficiencyErr")
        {
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_RefEfficiencyErr)(idx(ii)) = Values(i)(ii);
            }
        }
        else if(item == "ManualTail")
        {
            /*% if tail parameter changes, raise a flag
            Tl = Value{ManualTail_n};
            if ~isempty(Tl)
               oldTl = pat(idx, 10);
               pat(idx, 10) = Tl;
               subIdx = find(Tl ~= oldTl);
               if any(subIdx)
                  flags(idx(subIdx), 28) = 'C';
               end
            end*/
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_Tail)(idx(ii)) = Values(i)(ii);
            }
        }
        else if(item == "TailAlpha")
        {
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_TailAlpha)(idx(ii)) = Values(i)(ii);
            }
        }
        else if(item == "UpperTail")
        {
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_UpperTail)(idx(ii)) = Values(i)(ii);
            }
        }
        else if(item == "UpperTailAlpha")
        {
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_UpperTailAlpha)(idx(ii)) = Values(i)(ii);
            }
        }
        else if(item == "RecoilBeta")
        {
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_RecoilBeta)(idx(ii)) = Values(i)(ii);
            }
        }
        else if(item == "RecoilDeltaE")
        {
            for(uword ii=0; ii<idx.size(); ++ii)
            {
                pat.Peaks.col(i_RecoilDeltaE)(idx(ii)) = Values(i)(ii);
            }
        }
        ++i;
    }

    // reset area error, if not specified
    if( Descs.contains("NetArea") && !Descs.contains("AreaError") )//if isempty(Value{AreaError_n})
    {
        for(uword ii=0; ii<idx.size(); ++ii)
        {
            pat.Peaks.col(i_AreaError)(idx(ii)) = gNaN;
        }
    }
}

void GammaAnalyALG::FitBase(uvec cidx, uvec sidx)
{
    // get spectrum
    /* [R, baseWeights, Step, rg, cx, cy, cdy] = dmspec(sn, ...
        'ROISmooth', 'BaseFittingWeights', 'Steps', 'AnalysisRange', 'XControl', 'YControl', 'YSlope');*/
    QStringList descs;
    descs << "ROISmooth" << "BaseFittingWeights" << "Steps" << "AnalysisRange" << "XControl" << "YControl" << "YSlope";
    field<rowvec> t_f = dmspec(descs);
    rowvec &R = t_f(0);
    rowvec &baseWeights = t_f(1);
    rowvec &Step = t_f(2);
    rowvec &rg = t_f(3);
    rowvec cx = t_f(4);
    rowvec cy = t_f(5);
    rowvec cdy = t_f(6);

    // fit all control points, if no index given
    if(cidx.is_empty())
    {
        cidx = RangeVec(0, cy.size()-1).t();
    }
    /*if(sidx.is_empty())
    {
        sidx = [];
    }*/

    // get interesting channels
    rowvec x;
    if(rg.size() != 2) //if length(rg) ~= 2
    {
       qDebug() << "AnalysisRange not an interval";
       x = RangeVec2(1, R.size()); //x = [1 : length(R)];
    } else {
       x = RangeVec2(rg(0), rg(1)); //x = [rg(1) : rg(2)];
    }
    rowvec t_x = x - 1;
    rowvec y = IdxMat(R,t_x);
    rowvec w = IdxMat(baseWeights,t_x);

    uvec idx = find_nonfinite( cdy(sidx) );
    if( any(idx) )
    {
        MatAssign(cdy, sidx(idx), 0);
    }
    // get stepfree control points
    rowvec t_cx = cx - 1;
    cy = cy - IdxMat(Step, t_cx);
    // fit control points to residual
    /*[s, cx, cy, cdy] = fitFunction(x, y, 'wrapBaseSpline', ...
       cx, [], cy, cidx, cdy, sidx, ...
       'Weights', w, 'Waitbar', 'Baseline');*/
    mat t_cidx = UMatToMat(cidx);
    mat t_sidx = UMatToMat(sidx);

    double X2;
    field<vec> Ps;
    field<uvec> free_idx;
    Ps << vectorise(cx) << vectorise(cy) << vectorise(cdy);
    free_idx << uvec() << cidx << sidx;
    QStringList Opts; Opts << "Weights" << "Waitbar";
    Independ::fitFunction(Ps, X2, free_idx, vectorise(x), vectorise(y), "wrapBaseSpline", Opts, "Baseline", w);
    cx = Ps(0).t();
    cy = Ps(1).t();
    cdy = Ps(2).t();
    // add step back on
    t_cx = cx-1;
    cy = vectorise(cy).t() + IdxMat(Step,t_cx);

    //("FitBase_cx.txt", cx);
    //("FitBase_cy.txt", cy);
    //("FitBase_cdy.txt", cdy);

    // create temporary baseline
    tmpbase(cx, cy, cdy);
}

/* =====================================================谱基线结构======================================================== */
void GammaAnalyALG::msrd(mat &d, int &nzero, mat Y, mat F)
{
    // convert to column vectors
    vec yt = vectorise(Y);
    vec ft = vectorise(F);

    // check parameter size
    int ntotal = yt.size();
    if (ft.size() != ntotal)
    {
        qDebug() << "Sizes of parameters must match";
    }

    // find "non-zero" data points
    uvec idx = find(abs(yt) > eps(1.0));
    if(idx.is_empty()) // isempty(idx)
    {
        d << gNaN;
        nzero = ntotal;
    }
    else {
        int n = idx.size();
        vec y = yt(idx);
        vec r = ft(idx) - y;

        // calculate the msrd
        d = ( (r.t() / abs(y).t() ) * r) / n;
        nzero = ntotal - n;
    }
}

mat GammaAnalyALG::lSmoothen(rowvec y)
{
    // filter half width in units of FWHM
    double f = 1.5;

    rowvec s = y;

    // get resolution in channel units
    int nchan = y.size();
    rowvec c = RangeVec2(1, nchan); //c = cumsum(ones(size(y)));
    rowvec e = calValues(Cal_Energy, c);
    rowvec de = calDerivative(Cal_Energy, c);
    rowvec r = calValues(Cal_Resolution, e);

    // filter half width: factor times FWHM/channels
    rowvec n = max(2, ceil(f * r / de));

    // find positions where half width changes
    uvec crossover = find(diff(n) != 0) + 1u;   //crossover = find(diff(n) ~= 0);
    crossover.resize(crossover.size() + 1u);
    crossover( crossover.size() ) = nchan;      //crossover(end+1) = nchan;

    // filter piecewise
    rowvec smooth;
    int chlow = 1, chhigh, hw, dlow, dhigh, width;
    for(int i=0; i<crossover.size(); ++i) //for i = 1 : length(crossover)
    {
        // next upper border
        chhigh = crossover(i);
        hw = n(chhigh-1);

        // use true data instead of transient at borders
        dlow = min(hw, chlow-1);
        dhigh = min(hw, nchan - chhigh);

        width = 2*hw+1;
        if( (chhigh+dhigh - (chlow-dlow)) < width )
        {
            if(nchan < width)
            {
                // don't smoothen, but warn, should hardly ever happen
                // smooth = y(chlow-dlow : chhigh+dhigh);
                smooth = y.cols(chlow-dlow-1, chhigh+dhigh-1);
                qDebug() << QString("Can't smoothen: spectrum length: %1, frame width: %2").arg(nchan).arg(width);
            }
            else // increase data width
            {
                while( (chhigh+dhigh - (chlow-dlow)) < width )
                {
                    if(dhigh < nchan - chhigh) dhigh = dhigh+1;
                    if(dlow < chlow - 1) dlow = dlow+1;
                }
                // smoothen with wider data
                smooth = sgderive(y.cols(chlow-dlow-1, chhigh+dhigh-1), width, 0);
            }
        }
        else // smoothen
        {
            smooth = sgderive(y.cols(chlow-dlow-1, chhigh+dhigh-1), width, 0);
        }

        // take component where filter width constant
        // s(chlow:chhigh) = smooth(1+dlow : end-dhigh);
        s.cols(chlow-1,chhigh-1) = smooth.cols(dlow, smooth.size()-dhigh-1);

        // next region
        chlow = chhigh;
    }
    return s;
}

mat GammaAnalyALG::sgderive(mat y, int np, int order)
{
    mat s = zeros(size(y));
    y = vectorise(y); //y = y(:);

    // calculate filter
    rowvec f;
    int k;
    sgfilter(f, k, np, order); //[f, k] = sgfilter(np, order);
    //("sgderive_f.txt", f);

    // continue data symetrically
    // dd = [y(k+1:-1:2); y; y(end-1:-1:end-k-1)];
    int size_y = y.size();
    vec dd = join_vert( join_vert( y(RangeVec(k, 1, -1)), y ),
                        y( RangeVec(size_y-2, size_y-k-2, -1) ) );

    // fold
    urowvec idx = RangeVec(0, 2*k); //idx = [0:2*k];
    for(uword i=0; i<y.size(); ++i) //for i = 1:length(y)
    {
        s(i) = as_scalar(f * dd(i + idx)); //s(i) = f*(dd(i+idx));
    }
    return s;
}

void GammaAnalyALG::sgfilter(rowvec &f, int &k, int frameWidth, int deriveOrder)
{
    k = (frameWidth-1) / 2;
    rowvec F = RangeVec2(-k, k); //F = [-k : k];

    double S0; //S0 = lQS(k, 0);
    double S2 = accu( pow( RangeVec2(-k, k), 2) ); //S2 = lQS(k, 2);
    double S4 = accu( pow( RangeVec2(-k, k), 4) ); //S4 = lQS(k, 4);
    double S6;
    double N;
    switch(deriveOrder)
    {
    case 0: // f = ( S2*F2 - S4*F0 ) / ( S2^2 - S4*S0 )
            S0 = accu( pow( RangeVec2(-k, k), 0) );
            f = S2 * pow(F,2) - S4 * pow(F,0);
            N = -( S4*S0 - pow(S2,2) );
            f = f/N;
            break;

    case 2: // 2*f = ( S0*F2 - S2*F0 ) / ( S4*S0 - S2^2 )
            S0 = accu( pow( RangeVec2(-k, k), 0) );
            f = S0 * pow(F,2) - S2 * pow(F,0);
            N = ( S4*S0 - pow(S2,2) )/2;
            f = f/N;
            break;

    case 3: // 6*f = ( S2*F3 - S4*F1 ) / ( S6*S2 - S4^2 )
            S6 = accu( pow( RangeVec2(-k, k), 6) ); //S6 = lQS(k, 6);
            f = S2 * pow(F,3) - S4 * F;
            N = ( S6*S2 - pow(S4,2) )/6;
            f = f/N;
            break;

    default: qDebug() << "derivation Order not implemented yet";
    }
}

mat GammaAnalyALG::lSmoothROIs(rowvec y, mat c, mat w)
{
    // roi width in FWHM units
    int f_roi = 2;
    // smoothing filter width (9-point smoothing)
    int NSmooth = 19;
    //NMean = 11;

    int n = y.size();
    irowvec bool_smooth = zeros<irowvec>(n); //bool_smooth = false(size(y));
    mat low_roi = max(1, floor( c - f_roi * w ) );
    mat high_roi = min(n, ceil( c + f_roi * w ) );
    //("lSmoothROIs_low_roi.txt", low_roi);
    //("lSmoothROIs_high_roi.txt", high_roi);

    for(int i=0; i<low_roi.size(); ++i)
    {
        //bool_smooth(low_roi(i):high_roi(i)) = true;
        bool_smooth.cols(low_roi(i)-1, high_roi(i)-1)
                = ones<irowvec>(high_roi(i)-low_roi(i)+1);
    }
    bool_smooth(0) = false; // bool_smooth([1,end]) = false;
    bool_smooth(bool_smooth.size()-1) = false;

    irowvec d_smooth = diff(bool_smooth);
    uvec low_smooth = find(d_smooth == 1);
    uvec high_smooth = find(d_smooth == -1);

    //ss = sgderive(y(~bool_smooth), NMean, 0);
    rowvec smoothSpec = sgderive(y, NSmooth, 0);
    //("lSmoothROIs_smoothSpec.txt", smoothSpec);

    rowvec s = y;
    int skip = 0;
    for(int i=0; i<low_smooth.size(); ++i) // i = 1 : length(low_smooth)
    {
        int dx = high_smooth(i) - low_smooth(i) + 1;
        int l = max(1, low_smooth(i)-skip);
        skip = skip + dx;

        rowvec localSmoothData = smoothSpec.cols(low_smooth(i), high_smooth(i));
        s.cols(low_smooth(i), high_smooth(i)) = localSmoothData;// - mean(localSmoothData) + ss(l);
    }
    return s;
}

field<rowvec> GammaAnalyALG::dmspec(QStringList items)
{
    field<rowvec> Results(items.size());
    int i=0;
    rowvec BaseChan = cumsum(ones<rowvec>(m_nChans));

    /*[pC, pNA, pFC, pSt, pT, pTa, pUT, pUTa, pBWW, pRBC, pRDC] = dmps(sn, ...
        'Centroid', 'NetArea', 'FWHM_Ch', 'Step', ...
        'Tail', 'TailAlpha', 'UpperTail', 'UpperTailAlpha', 'BWWidthChan', ...
        'RecoilBetaChan', 'RecoilDeltaChan');*/
    vec pC, pNA, pFC, pSt, pT, pTa, pUT, pUTa, pBWW, pRBC, pRDC, z, o;
    QStringList Opts;
    if(items.contains("BaseLine") || items.contains("Approximation") || items.contains("Stripped")
            || items.contains("Steps") || items.contains("ROISmooth") || items.contains("BaseFittingWeights"))
    {
        Opts << "Centroid" << "NetArea" << "FWHM_Ch" << "Step" << "Tail" << "TailAlpha" << "UpperTail"
             << "UpperTailAlpha" << "BWWidthChan" << "RecoilBetaChan" << "RecoilDeltaChan";
    }
    if(!Opts.isEmpty())
    {
        field<vec> f = dmps(Opts);
        pC = f(0);  pNA = f(1);  pFC = f(2);  pSt = f(3);  pT = f(4);    pTa = f(5);
        pUT = f(6); pUTa = f(7); pBWW = f(8); pRBC = f(9); pRDC = f(10);
        z = zeros(size(pC));
        o = ones(size(pC));
    }

    if(items.contains("Approximation") || items.contains("Residual") || items.contains("ChiSquare"))
    {
        PiecePoly bpp = Independ::makeBasePP(vectorise(baseInfo(s_XControl)), vectorise(baseInfo(s_YControl)),
                                   vectorise(baseInfo(s_YSlope)), pC, pFC, pSt);
        baseInfo(s_Approximation) = Independ::peakBaseSpline(BaseChan, bpp, pNA, pC, pFC, pSt, pT, pTa, pUT, pUTa, pBWW, pRBC, pRDC);
    }
    if(items.contains("Stripped") || items.contains("SmoothStripped") || items.contains("ROISmooth"))
    {
        if(!pC.is_empty())
        {
            PiecePoly bpp;
            rowvec pk = Independ::peakBaseSpline(BaseChan, bpp, pNA, pC, pFC, z, pT, pTa, pUT, pUTa, pBWW, pRBC, pRDC);

//			PrintMat22("BaseChan.txt", BaseChan);
//			PrintMat22("pNA.txt", pNA);
//			PrintMat22("pC.txt", pC);
//			PrintMat22("pFC.txt", pFC);
//			PrintMat22("z.txt", z);
//			PrintMat22("pT.txt", pT);
//			PrintMat22("pTa.txt", pTa);
//			PrintMat22("pUT.txt", pUT);
//			PrintMat22("pUTa.txt", pUTa);
//			PrintMat22("pBWW.txt", pBWW);
//			PrintMat22("pRBC.txt", pRBC);
//			PrintMat22("pRDC.txt", pRDC);

//			PrintMat22("baseInfo(s_Spectrum).txt", baseInfo(s_Spectrum));
//			PrintMat22("pk.txt", pk);

            baseInfo(s_Stripped) = baseInfo(s_Spectrum) - pk;
        }
        else baseInfo(s_Stripped) = baseInfo(s_Spectrum);
    }
    if(items.contains("Steps") || items.contains("ROISmooth"))
    {
        baseInfo(s_Steps) = Independ::peakSteps(BaseChan, pC, pFC, pSt);
    }

    foreach (QString item, items)
    {
        if(item == "Spectrum" || item == "LongSpectrum" || item == "ShortSpectrum")
        {
            Results(i) = baseInfo(s_Spectrum);
        }
        else if(item == "XControl")
        {
            Results(i) = baseInfo(s_XControl);
        }
        else if(item == "YControl")
        {
            Results(i) = baseInfo(s_YControl);
        }
        else if(item == "YSlope")
        {
            Results(i) = baseInfo(s_YSlope);
        }
        else if(item == "AnalysisRange")
        {
            Results(i) = baseInfo(s_AnalysisRange);
        }
        else if(item == "AnalysisRangeSave")
        {
            if(baseInfo(s_AnalysisRange).is_empty())
            {
                rowvec r; r << 1 << m_nChans;
                Results(i) = r;
            }
            else {
                Results(i) = baseInfo(s_AnalysisRange);
            }
        }
        else if(item == "BaseLine")
        {
            if(pat.Peaks.n_cols == MAX_INDEX_PEAK)
            {
                PiecePoly bpp = Independ::makeBasePP(vectorise(baseInfo(s_XControl)), vectorise(baseInfo(s_YControl)),
                                           vectorise(baseInfo(s_YSlope)), pC, pFC, pSt);
                Results(i) = Matlab::ppval(BaseChan, bpp)+Independ::peakSteps(BaseChan, pC, pFC, pSt);
                baseInfo(s_BaseLine) = Results(i);
            }
        }
        else if(item == "BaseFittingWeights")
        {
            //BaseFittingWeights = lGetWeights(BaseChan);
            // get peak data
            //[C, NA, FC, St] = dmps(sn, 'Centroid', 'NetArea', 'FWHM_Ch', 'Step');

            // get fitting weights: punish regions near peaks
            mat w = ones(size(BaseChan));
            for(int i=0; i<pC.size(); ++i) // for i = 1 : length(C)
            {
                if(pNA(i) > 1) // if NA(i) > 1
                {
                    // w = w + sqrt(NA(i)) ./ ((max(1, abs((x-C(i))/FC(i)))).^3 );
                    w = w + sqrt(pNA(i)) / pow( max(1, abs((BaseChan-pC(i))/pFC(i))), 3);
                }
            }
            Results(i) = 1 / (abs(w) + 1);
            baseInfo(s_BaseFittingWeights) = Results(i);
        }
        else if(item == "Approximation")
        {
            Results(i) = baseInfo(s_Approximation);
        }
        else if(item == "Residual")
        {
            Results(i) = baseInfo(s_Spectrum) - baseInfo(s_Approximation);
            baseInfo(s_Residual) = Results(i);
        }
        else if(item == "Stripped")
        {
            Results(i) = baseInfo(s_Stripped);
        }
        else if(item == "PSS")
        {
            rowvec e = calValues(Cal_Energy, BaseChan);
            rowvec dE = calDerivative(Cal_Energy, BaseChan);
            rowvec r = calValues(Cal_Resolution, e);
            Results(i) = -pksens(baseInfo(s_Spectrum), r/dE, BaseChan, BaseChan);
            baseInfo(s_PSS) = Results(i);
        }
        else if(item == "SmoothStripped")
        {
            Results(i) = lSmoothen(baseInfo(s_Stripped));
            baseInfo(s_SmoothStripped) = Results(i);
        }
        else if(item == "ChiSquare")
        {
            if(baseInfo(s_AnalysisRange).is_empty()) // isempty(AnalysisRange)
            {
                Results(i) << gNaN << gNaN; // ChiSquare = [NaN, NaN];
            }
            else {
                //ch = AnalysisRange(1) : AnalysisRange(2);
                rowvec F = baseInfo(s_Approximation).cols(baseInfo(s_AnalysisRange)(0), baseInfo(s_AnalysisRange)(1));
                rowvec S = baseInfo(s_Spectrum).cols(baseInfo(s_AnalysisRange)(0), baseInfo(s_AnalysisRange)(1));

                mat MSRD;
                int nzero;
                msrd(MSRD, nzero, S, F); // [MSRD, nzero] = msrd(S, F);

                // return ChiSquare and number of channels with zero counts
                rowvec tmp; tmp << nzero;
                Results(i) = join_horiz(MSRD, tmp); //ChiSquare=[ MSRD, nzero];
            }
            baseInfo(s_ChiSquare) = Results(i);
        }
        /*else if(item == "ArtXControl")
        {
            Results(i) = baseInfo(s_ArtXControl);
        }
        else if(item == "ArtYControl")
        {
            Results(i) = baseInfo(s_ArtYControl);
        }
        else if(item == "ArtYSlope")
        {
            Results(i) = baseInfo(s_ArtYSlope);
        }
        else if(item == "ArtificialBase")
        {
            if(baseInfo(s_ArtXControl).is_empty())
            {
                Results(i) = zeros(size(baseInfo(s_Spectrum)));
            }
            else //if ~isempty(artbasestruc)
            {
                mat empty = mat(0, 0);
                PiecePoly ArtBaselinePP = makeBasePP(baseInfo(s_ArtXControl), baseInfo(s_ArtYControl),
                                                     baseInfo(s_ArtYSlope), empty, empty, empty);
                Results(i) = peakBaseSpline(BaseChan, ArtBaselinePP, empty, empty, empty, empty,
                                            empty, empty, empty, empty, empty, empty, empty);
            }
        }*/
        else if(item == "Steps")
        {
            Results(i) = baseInfo(s_Steps);
        }
        /*else if(item == "SteplessBase")
        {
            Results(i) = baseInfo(s_Stripped) - baseInfo(s_Steps);
        }*/
        else if(item == "ROISmooth")
        {
            Results(i) = lSmoothROIs(baseInfo(s_Stripped) - baseInfo(s_Steps), pC, pFC);
            baseInfo(s_ROISmooth) = Results(i);
        }
        ++i;
    }

    /*for(int i=0; i<Results.size(); ++i)
    {
        QString name = QString("dmspec_%1.txt").arg(items.at(i));
        //(name, Results(i));
    }*/

    return Results;
}

/* ====================================================PAT结构======================================================== */
mat GammaAnalyALG::calValuesByName(CalibrationType cal, mat x)
{
    mat y = zeros(size(x));
    // get named calibration parameter
    vec p;
    if(getCalibrationPara(p, cal)) //[s, msg, p] = getCalibrationPara(p, cal, name);
    {
        // evaluate calibration function
        y = Independ::calFcnEval(x, p); //[y, s] = calFcnEval(x, p);
    }
    return y;
}

mat GammaAnalyALG::calDerivativeByName(CalibrationType cal, mat x)
{
    mat dy = zeros(size(x));
    // get named calibration parameter
    colvec p;
    if(getCalibrationPara(p, cal)) //[s, msg, p] = getCalibrationPara(sn, cal, name);
    {
        // evaluate slope
        dy = Independ::calDerivEval(x, p); //[dy, s] = calDerivEval(x, p);
    }
    return dy;
}

void GammaAnalyALG::dmps()
{
    QStringList items;
    items << "Energy" << "FWHM" << "Multiplet" << "AreaError" << "Significance" << "Tail"
          << "TailAlpha" << "UpperTail" << "UpperTailAlpha" << "StepRatio" << "Efficiency";
    field<vec> results = dmps(items);
    //results.print("results");
    pat.Peaks.col(i_Energy) = results(0);
    pat.Peaks.col(i_FWHM) = results(1);
    pat.Peaks.col(i_Multiplet) = results(2);
    pat.Peaks.col(i_AreaError) = results(3);
    pat.Peaks.col(i_Significance) = results(4);
    pat.Peaks.col(i_Tail) = results(5);
    pat.Peaks.col(i_TailAlpha) = results(6);
    pat.Peaks.col(i_UpperTail) = results(7);
    pat.Peaks.col(i_UpperTailAlpha) = results(8);
    pat.Peaks.col(i_StepRatio) = results(9);
    pat.Peaks.col(i_Efficiency) = results(10);
}

vec GammaAnalyALG::dmps(PeakIndex idx)
{
    vec vRet;
    vec Centroid = pat.Peaks.col(i_Centroid);
    if(idx == i_Energy)
    {
        vRet = pat.Peaks.col(i_Energy); //Energy = EnergyPAT;
        uvec Bool = find_nonfinite(vRet);
        if( !Bool.is_empty() )
        {
            //Energy(bool) = calValuesByName(sn, 'Energy', calNames{1}, Centroid(bool));
            vRet(Bool) = calValuesByName(Cal_Energy, Centroid(Bool));
        }
    }
    else if(idx == i_Multiplet)
    {
        vec Multiplet = pat.Peaks.col(i_Multiplet);
        uvec Bool = find_nonfinite(Multiplet);
        if ( !Bool.is_empty() )
        {
            vec FWHM_Ch = pat.Peaks.col(i_FWHM_Ch); //FWHM_Ch = FWHM_ChPAT;
            uvec Bool2 = find_nonfinite(FWHM_Ch); //bool = isnan(FWHM_Ch);
            if( !Bool2.is_empty() ) //if any(bool)
            {
                vec Energy = pat.Peaks.col(i_Energy); //Energy = EnergyPAT;
                uvec Bool3 = find_nonfinite(Energy);
                if( !Bool3.is_empty() )
                {
                    Energy(Bool3) = calValuesByName(Cal_Energy, Centroid(Bool3));
                }
                vec Res = calValuesByName(Cal_Resolution, Energy);
                vec dE = calDerivativeByName(Cal_Energy, Centroid);
                vec Res_Ch = Res / dE; //Res_Ch=Res./dE;
                FWHM_Ch(Bool2) = Res_Ch(Bool2); //FWHM_Ch(bool) = Res_Ch(bool);
            }
            // where the PAT entry is NaN, set calculated flag; for calculation,
            // all peaks are needed !!
            vec tmpMultiplet = peakFindMultiplets(Centroid, pat.Peaks.col(i_NetArea), FWHM_Ch, 2);
            Multiplet(Bool) = tmpMultiplet(Bool);
        }
        vRet = Multiplet;
        //vRet.print("vRet = ");
    }
    else if(idx == i_FWHM)
    {
        vec FWHM = pat.Peaks.col(i_FWHM);
        uvec Bool = find_nonfinite(FWHM);
        if(!Bool.is_empty())
        {
            vec FWHM_Ch = pat.Peaks.col(i_FWHM_Ch); //FWHM_Ch = FWHM_ChPAT;
            vec dE = calDerivativeByName(Cal_Energy, Centroid);
            uvec Bool2 = find_nonfinite(FWHM_Ch); //bool = isnan(FWHM_Ch);
            if( !Bool2.is_empty() ) //if any(bool)
            {
                vec Energy = pat.Peaks.col(i_Energy); //Energy = EnergyPAT;
                uvec Bool3 = find_nonfinite(Energy);
                if( !Bool3.is_empty() )
                {
                    Energy(Bool3) = calValuesByName(Cal_Energy, Centroid(Bool3));
                }
                vec Res = calValuesByName(Cal_Resolution, Energy);
                vec Res_Ch = Res / dE; //Res_Ch=Res./dE;
                FWHM_Ch(Bool2) = Res_Ch(Bool2); //FWHM_Ch(bool) = Res_Ch(bool);
            }
            FWHM(Bool) = FWHM_Ch(Bool) % dE(Bool);
        }
        vRet = FWHM;
    }
    else if(idx == i_AreaError)
    {
        //AreaError = lCalcAreaError(AreaErrorPAT, NetArea, BackgroundArea, LC);
        vec dNA = pat.Peaks.col(i_AreaError);
        uvec Bool = find_nonfinite(dNA);
        if(!Bool.is_empty())
        {
            vec NA = pat.Peaks.col(i_NetArea);
            vec LC = pat.Peaks.col(i_LC);
            vec BA = pat.Peaks.col(i_BackgroundArea);
            dNA(Bool) = sqrt( max( NA(Bool), LC(Bool) ) + BA(Bool) );
        }
        vRet = dNA;
    }
    else if(idx == i_Significance)
    {
        //Significance = lCalcSignificance(SignificancePAT, NetArea, LC);
        vec signif = pat.Peaks.col(i_Significance);
        uvec Bool = find_nonfinite(signif);
        if(!Bool.is_empty())
        {
            vec NA = pat.Peaks.col(i_NetArea);
            vec LC = pat.Peaks.col(i_LC);
            signif(Bool) = NA(Bool) / LC(Bool);
        }
        // Do not report significance for neutron lumps
        //Significance(Source == 'G') = NaN;
        vRet = signif;
    }
    else if(idx == i_Efficiency)
    {
        vec vRet = pat.Peaks.col(i_Efficiency);
        uvec Bool = find_nonfinite(vRet);
        if(!Bool.is_empty())
        {
            vec Energy = pat.Peaks.col(i_Energy); //Energy = EnergyPAT;
            uvec Bool2 = find_nonfinite(Energy);
            if( !Bool2.is_empty() )
            {
                //Energy(bool) = calValuesByName(sn, 'Energy', calNames{1}, Centroid(bool));
                Energy(Bool2) = calValuesByName(Cal_Energy, Centroid(Bool2));
            }
            vRet(Bool) = calValuesByName(Cal_Efficiency, Energy(Bool));
        }
    }

    return vRet;
}

field<vec> GammaAnalyALG::dmps(QStringList items)
{
    field<vec> Results(items.size());
    vec FWHM_Ch, StepRatio, Res, dE, Res_Ch, Energy, Centroid = pat.Peaks.col(i_Centroid);

    bool bFwhmCh = false, bdE = false;
    /* 需要 FWHM_Ch 的项：FWHM、SigmaCh、UpperTailBeta、Multiplet */
    if(items.contains("FWHM_Ch") || items.contains("FWHM") || items.contains("SigmaCh") || items.contains("Multiplet"))
    {
        bFwhmCh = true;
    }

    /* 需要 dE 的项：Res_Ch、FWHM、UpperTailBeta、BWWidthChan、RecoilBetaChan、RecoilDeltaChan
     * 1、Res_Ch：FWHM_Ch
     * 2、FWHM_Ch：FWHM、SigmaCh、UpperTailBeta、Multiplet*/
    if(bFwhmCh || items.contains("BWWidthChan") || items.contains("RecoilBetaChan") || items.contains("RecoilDeltaChan"))
    {
        bdE = true;
    }

    /* 需要 Energy 的项：Res、Efficiency、EfficiencyErr、Tail、TailAlpha、UpperTail、UpperTailAlpha、StepRatio、LineDeviation
     * 1、Res：Res_Ch
     * 2、Res_Ch：FWHM_Ch */
    if(bdE || items.contains("Energy") || items.contains("Res") || items.contains("Tail")
           || items.contains("TailAlpha") || items.contains("UpperTail") || items.contains("UpperTailAlpha")
           || items.contains("StepRatio") || items.contains("Efficiency"))
    {
        Energy = pat.Peaks.col(i_Energy); //Energy = EnergyPAT;
        uvec Bool = find_nonfinite(Energy);
        if( !Bool.is_empty() )
        {
            //Energy(bool) = calValuesByName(sn, 'Energy', calNames{1}, Centroid(bool));
            Energy(Bool) = calValuesByName(Cal_Energy, Centroid(Bool));
        }
    }

    if(bdE)
    {
        //dE = calDerivativeByName(sn, 'Energy', calNames{1}, Centroid);
        dE = calDerivativeByName(Cal_Energy, Centroid);
    }

    if(bFwhmCh)
    {
        FWHM_Ch = pat.Peaks.col(i_FWHM_Ch); //FWHM_Ch = FWHM_ChPAT;
//		PrintMat22("FWHM_CH.txt", FWHM_Ch);
        uvec Bool = find_nonfinite(FWHM_Ch); //bool = isnan(FWHM_Ch);
        if( !Bool.is_empty() ) //if any(bool)
        {
            //Res = calValuesByName(sn, 'Resolution', calNames{2}, Energy);
            Res = calValuesByName(Cal_Resolution, Energy);
            Res_Ch = Res / dE; //Res_Ch=Res./dE;
            FWHM_Ch(Bool) = Res_Ch(Bool); //FWHM_Ch(bool) = Res_Ch(bool);
        }
    }

    if(items.contains("StepRatio") || items.contains("Step"))
    {
        StepRatio = pat.Peaks.col(i_StepRatio); //StepRatio = StepRatioPAT;
        uvec Bool = find_nonfinite(StepRatio); //bool = isnan(StepRatio);
        if( !Bool.is_empty() ) //if any(bool)
        {
            StepRatio(Bool) = calValuesByName(Cal_Step_ratio, Energy(Bool));
        }
    }

    int i=0;
    foreach (QString item, items)
    {
        if(item == "Centroid")
        {
            Results(i) = Centroid;
        }
        else if(item == "Energy")
        {
            Results(i) = Energy;
        }
        else if(item == "Res")
        {
            if(Res.is_empty()) Results(i) = calValuesByName(Cal_Resolution, Energy);
            else Results(i) = Res;
        }
        else if(item == "FWHM")
        {
            vec FWHM = pat.Peaks.col(i_FWHM);
            uvec Bool = find_nonfinite(FWHM);
            if(!Bool.is_empty())
            {
                FWHM(Bool) = FWHM_Ch(Bool) % dE(Bool);
            }
            Results(i) = FWHM;
        }
        else if(item == "FWHM_Ch")
        {
            Results(i) = FWHM_Ch;
        }
        else if(item == "NetArea")
        {
            Results(i) = pat.Peaks.col(i_NetArea);
        }
        else if(item == "AreaError")
        {
            //AreaError = lCalcAreaError(AreaErrorPAT, NetArea, BackgroundArea, LC);
            vec dNA = pat.Peaks.col(i_AreaError);
            uvec Bool = find_nonfinite(dNA);
            if(!Bool.is_empty())
            {
                vec NA = pat.Peaks.col(i_NetArea);
                vec LC = pat.Peaks.col(i_LC);
                vec BA = pat.Peaks.col(i_BackgroundArea);
                dNA(Bool) = sqrt( max( NA(Bool), LC(Bool) ) + BA(Bool) );
            }
            Results(i) = dNA;
        }
        else if(item == "Significance")
        {
            //Significance = lCalcSignificance(SignificancePAT, NetArea, LC);
            vec signif = pat.Peaks.col(i_Significance);
            uvec Bool = find_nonfinite(signif);
            if(!Bool.is_empty())
            {
                vec NA = pat.Peaks.col(i_NetArea);
                vec LC = pat.Peaks.col(i_LC);
                signif(Bool) = NA(Bool) / LC(Bool);
            }
            // Do not report significance for neutron lumps
            //Significance(Source == 'G') = NaN;
            Results(i) = signif;
        }
        else if(item == "Sensitivity")
        {
            Results(i) = pat.Peaks.col(i_Sensitivity);
        }
        else if(item == "Multiplet")
        {
            vec Multiplet = pat.Peaks.col(i_Multiplet);
            uvec Bool = find_nonfinite(Multiplet);
            if ( !Bool.is_empty() )
            {
                // where the PAT entry is NaN, set calculated flag; for calculation,
                // all peaks are needed !!
                vec tmpMultiplet = peakFindMultiplets(Centroid, pat.Peaks.col(i_NetArea), FWHM_Ch, 2);
                Multiplet(Bool) = tmpMultiplet(Bool);
            }
            Results(i) = Multiplet;
        }
        else if(item == "Index")
        {
            Results(i) = RangeVec2(1, pat.Peaks.n_rows).t(); //Index = [1 : NPeaks]';
        }
        else if(item == "NetAreaFree")
        {
            Results(i) = ones<vec>(Centroid.size());
        }
        else if(item == "SigmaCh")
        {
            Results(i) = FWHM_Ch/sqrt(8*log(2)); //SigmaCh = FWHM_Ch/sqrt(8*log(2));
        }
        else if(item == "Step")
        {
            Results(i) = pat.Peaks.col(i_NetArea) % StepRatio; //Step = NetArea .* StepRatio;
        }
        else if(item == "StepRatio")
        {
            Results(i) = StepRatio;
        }
        else if(item == "Tail")
        {
            Results(i) = pat.Peaks.col(i_Tail); //Tail = ManualTail;
            uvec Bool = find_nonfinite(Results(i)); //bool = isnan(Tail);
            if( !Bool.is_empty() )
            {
                Results(i)(Bool) = calValuesByName(Cal_Tail, Energy(Bool));
            }
        }
        else if(item == "TailAlpha")
        {
            Results(i) = pat.Peaks.col(i_TailAlpha); //TailAlpha = TailAlphaPAT;
            uvec Bool = find_nonfinite(Results(i));
            if( !Bool.is_empty() )
            {
                Results(i)(Bool) = calValuesByName(Cal_Tail_alpha, Energy(Bool));
            }
        }
        else if(item == "UpperTail")
        {
            Results(i) = pat.Peaks.col(i_UpperTail); //TailAlpha = TailAlphaPAT;
            uvec Bool = find_nonfinite(Results(i));
            if( !Bool.is_empty() )
            {
                Results(i)(Bool) = calValuesByName(Cal_Tail_right, Energy(Bool));
            }
        }
        else if(item == "UpperTailAlpha")
        {
            Results(i) = pat.Peaks.col(i_UpperTailAlpha); //TailAlpha = TailAlphaPAT;
            uvec Bool = find_nonfinite(Results(i));
            if( !Bool.is_empty() )
            {
                Results(i)(Bool) = calValuesByName(Cal_Tail_right_alpha, Energy(Bool));
            }
        }
        else if(item == "BWWidthChan")
        {
            vec BWWidth = pat.Peaks.col(i_BWWidth);
            uvec Bool = find_nonfinite(BWWidth);
            if( !Bool.is_empty() )
            {
                MatAssign(BWWidth, Bool, 0);
            }
            Results(i) = BWWidth / dE;
        }
        else if(item == "RecoilBetaChan")
        {
            Results(i) = pat.Peaks.col(i_RecoilBeta) % dE; //RecoilBetaChan = RecoilBeta .* dE;
        }
        else if(item == "RecoilDeltaChan")
        {
            Results(i) = pat.Peaks.col(i_RecoilDeltaE) / dE; //RecoilDeltaChan = RecoilDeltaE ./ dE;
        }
        else if(item == "Efficiency")
        {
            vec effi = pat.Peaks.col(i_Efficiency);
            uvec Bool = find_nonfinite(effi);
            if(!Bool.is_empty())
            {
                effi(Bool) = calValuesByName(Cal_Efficiency, Energy(Bool));
            }
            Results(i) = effi;
        }
        else if(item == "EfficiencyErr")
        {
            vec eff_err = pat.Peaks.col(i_EfficiencyErr);
            uvec Bool = find_nonfinite(eff_err);
            if(!Bool.is_empty())
            {
                //[s, msg, p, perr] = calGetPATPara(sn, 'Efficiency', patName);
                vec p;
                if( calGetPATPara(p, Cal_Efficiency) )
                {
                    // relative error in percent
                    //eff_err(Bool) = 100 * calErrorEval(Energy(Bool), p, perr);
                }
            }
        }
        ++i;
    }

    return Results;
}

field<vec> GammaAnalyALG::dmps(QList<PeakIndex> list_idx)
{
    //pat.Peaks.print("132465");
    field<vec> Results(list_idx.size());
    vec FWHM_Ch, Res, dE, Res_Ch, Energy, Centroid = pat.Peaks.col(i_Centroid);

    bool bFwhmCh = false, bdE = false;
    /* 需要 FWHM_Ch 的项：FWHM、SigmaCh、UpperTailBeta、Multiplet */
    if(list_idx.contains(i_FWHM_Ch) || list_idx.contains(i_FWHM) || list_idx.contains(i_Multiplet))
    {
        bFwhmCh = true;
        bdE = true;
    }

    /* 需要 Energy 的项：Res、Efficiency、EfficiencyErr、Tail、TailAlpha、UpperTail、UpperTailAlpha、StepRatio、LineDeviation
     * 1、Res：Res_Ch    2、Res_Ch：FWHM_Ch */
    if(bdE || list_idx.contains(i_Energy))
    {
        Energy = pat.Peaks.col(i_Energy); //Energy = EnergyPAT;
        uvec Bool = find_nonfinite(Energy);
        if( !Bool.is_empty() )
        {
            //Energy(bool) = calValuesByName(sn, 'Energy', calNames{1}, Centroid(bool));
            Energy(Bool) = calValuesByName(Cal_Energy, Centroid(Bool));
        }
        //("dmps_Energy.txt", Energy);
    }

    if(bdE)
    {
        //dE = calDerivativeByName(sn, 'Energy', calNames{1}, Centroid);
        dE = calDerivativeByName(Cal_Energy, Centroid);
        //("dmps_dE.txt", dE);
    }

    if(bFwhmCh)
    {
        FWHM_Ch = pat.Peaks.col(i_FWHM_Ch); //FWHM_Ch = FWHM_ChPAT;
        uvec Bool = find_nonfinite(FWHM_Ch); //bool = isnan(FWHM_Ch);
        if( !Bool.is_empty() ) //if any(bool)
        {
            //Res = calValuesByName(sn, 'Resolution', calNames{2}, Energy);
            Res = calValuesByName(Cal_Resolution, Energy);
            Res_Ch = Res / dE; //Res_Ch=Res./dE;
            FWHM_Ch(Bool) = Res_Ch(Bool); //FWHM_Ch(bool) = Res_Ch(bool);
        }
    }

    int i=0;
    foreach (PeakIndex idx, list_idx)
    {
        if(idx == i_Centroid)
        {
            Results(i) = Centroid;
        }
        else if(idx == i_Energy)
        {
            Results(i) = Energy;
        }
        else if(idx == i_Multiplet)
        {
            vec Multiplet = pat.Peaks.col(i_Multiplet);
            uvec Bool = find_nonfinite(Multiplet);
            if ( !Bool.is_empty() )
            {
                // where the PAT entry is NaN, set calculated flag; for calculation,
                // all peaks are needed !!
                vec tmpMultiplet = peakFindMultiplets(Centroid, pat.Peaks.col(i_NetArea), FWHM_Ch, 2);
                Multiplet(Bool) = tmpMultiplet(Bool);
            }
            Results(i) = Multiplet;
        }
        else if(idx == i_FWHM)
        {
            vec FWHM = pat.Peaks.col(i_FWHM);
            uvec Bool = find_nonfinite(FWHM);
            if(!Bool.is_empty())
            {
                FWHM(Bool) = FWHM_Ch(Bool) % dE(Bool);
            }
            Results(i) = FWHM;
        }
        else if(idx == i_NetArea)
        {
            Results(i) = pat.Peaks.col(i_NetArea);
        }
        else if(idx == i_AreaError)
        {
            //AreaError = lCalcAreaError(AreaErrorPAT, NetArea, BackgroundArea, LC);
            vec dNA = pat.Peaks.col(i_AreaError);
            uvec Bool = find_nonfinite(dNA);
            if(!Bool.is_empty())
            {
                vec NA = pat.Peaks.col(i_NetArea);
                vec LC = pat.Peaks.col(i_LC);
                vec BA = pat.Peaks.col(i_BackgroundArea);
                dNA(Bool) = sqrt( max( NA(Bool), LC(Bool) ) + BA(Bool) );
            }
            Results(i) = dNA;
        }
        else if(idx == i_Significance)
        {
            //Significance = lCalcSignificance(SignificancePAT, NetArea, LC);
            vec signif = pat.Peaks.col(i_Significance);
            uvec Bool = find_nonfinite(signif);
            if(!Bool.is_empty())
            {
                vec NA = pat.Peaks.col(i_NetArea);
                vec LC = pat.Peaks.col(i_LC);
                signif(Bool) = NA(Bool) / LC(Bool);
            }
            // Do not report significance for neutron lumps
            //Significance(Source == 'G') = NaN;
            Results(i) = signif;
        }
        else if(idx == i_Sensitivity)
        {
            Results(i) = pat.Peaks.col(i_Sensitivity);
        }
        ++i;
    }

    return Results;
}

vec GammaAnalyALG::peakFindMultiplets(vec C, vec NA, vec FC, double RW)
{
    uword NPeaks = C.size();
    vec Mult = zeros<vec>(NPeaks);

    vec CL = floor(C - RW * FC);
    vec CR = ceil(C + RW * FC);

    for (uword i=0; i<NPeaks; ++i)
    {
        uvec idx = find(CL <= CR(i) && CR >= CL(i));
        if (idx.size() > 1)
        {
            int L = min(CL(idx));
            int R = max(CR(idx));
            CL(idx) = ones<vec>(size(idx)) * L; //CL(idx) = L;
            CR(idx) = ones<vec>(size(idx)) * R; //CR(idx) = R;
            Mult(idx) = ones<vec>(size(idx));    //Mult(idx) = 1;
        }
    }
    return Mult;
}

/* ====================================================刻度更新======================================================== */
void GammaAnalyALG::calUpdate(QChar dataType, vec certEne, bool E1, bool R, bool E2)
{
    // if true, peaks found in peak search will be kept
    //write_peaks_to_pat = getSpecSetup(sn, 'KeepCalPeakSearchPeaks');
    bool write_peaks_to_pat = specSetup.KeepCalPeakSearchPeaks;

    // for calibration spectra, use energies from certificate instead of setup file
    //[dataType, certEne] = dminfo(sn, 'DataType', 'CertificateGEnergies');
    vec ELibEne, ELibRes;
    if(dataType == 'C' && !certEne.is_empty())
    {
        ELibEne = certEne;
        ELibRes = certEne;
    }
    else // getting library energy lines
    {
        //QString filename = qApp->applicationDirPath() + "/setup/ene_cal_update.txt";
        QString filename = m_strFilePath + "/ene_cal_update.txt";

        QFile file(filename);
        if(file.open(QIODevice::ReadOnly))
        {
            QTextStream in(&file);

            mat tab;
            rowvec t_r;
            QString line = in.readLine().trimmed();
            uword n_row = 0;
            while(!in.atEnd())
            {
                if(!line.isEmpty() && line[0] != '%')
                {
                    QStringList list = line.split(" ", QString::SkipEmptyParts);
                    if(list.size() >= 5)
                    {
                        t_r << list[0].toDouble() << list[1].toDouble() << list[2].toDouble() << list[3].toDouble() << list[4].toDouble();
                        tab.resize(++n_row, 5u);
                        tab.row(n_row-1u) = t_r;
                    }
                }
                line = in.readLine().trimmed();
            }
            file.close();

            // lines for energy calibration: flag ~= 0
            uvec eneIdx = find(tab.col(0));
            ELibEne = tab(eneIdx, uvec("1"));
            // lines for resolution calibration: flag == 2
            uvec resIdx = find(tab.col(0) == 2);
            ELibRes = tab(resIdx, uvec("1"));
        }

        if(ELibEne.is_empty() || ELibRes.is_empty())
        {
            return;
        }
    }

    ////("calUpdate_ELibEne.txt", ELibEne);
    ////("calUpdate_ELibRes.txt", ELibRes);

    if(E1)
    {
        // first energy calibration update
        vec ELibEne1;
        uvec ene1PIdx = calPeakSearch(ELibEne1, ELibEne);
        calEnergyUpdate1(ene1PIdx, ELibEne1);
    }

    if(R)
    {
        // update resolution calibration
        vec ELibRes1;
        uvec resPIdx = calPeakSearch(ELibRes1, ELibRes);
        calResUpdate(resPIdx);
    }

    if(E2)
    {
        // second energy calibration update
        vec ELibEne2;
        uvec ene2PIdx = calPeakSearch(ELibEne2, ELibEne);
        calEnergyUpdate2(ene2PIdx, ELibEne2);
    }

    // keep changes (quickfit) in PAT, if requested
    if(write_peaks_to_pat)
    {
        //patCommit(sn);
        //tbaseac(sn);
        colvec sqi, uqi;
        patBaseVar(sqi, uqi, uvec(), 1);
    }
    else
    {
        //patRollback;
        pat.Peaks.clear();
        pat.Peaks.set_size(0, MAX_INDEX_PEAK);

        //tbaseac('reject');
        baseInfo(s_XControl).clear();
        baseInfo(s_YControl).clear();
        baseInfo(s_YSlope).clear();
        baseInfo(s_AnalysisRange).clear();
    }

}

uvec GammaAnalyALG::calPeakSearch(vec &EF, vec ELib)
{
    // max deviation (oldCalibEnergy(PATCentroid) - LibraryLine) < EDevMaxGain*Res
    double EDevMaxGain = 0.7;
    // first try to find NPeaksDesired peaks, possibly at lower PSS
    double NPeaksDesired = 6;
    // if not enough peaks, restart at high PSS, and search NPeaksMin peaks
    double NPeaksMin = 3;
    // Peak search sensitivity
    double PSSs = 4;

    // get PSS tresholds: high use level, minimum use level
    double PSSh = specSetup.CalibrationPSS_high;
    double PSSl = specSetup.CalibrationPSS_low;

    // get peaks from PAT or do peak search
    //[allIdx, C, E, NA, pss, Mult, res] = dmps(sn, 'Index', 'Centroid', 'Energy', 'NetArea', 'Sensitivity', 'Multiplet', 'Res');
    QStringList Opts;
    Opts << "Index" << "Centroid" << "Energy" << "NetArea" << "Sensitivity" << "Multiplet" << "Res";
    field<vec> t_f = dmps(Opts);

    uvec allIdx;
    vec E = t_f(2), res = t_f(6), pss = t_f(4), Mult = t_f(5);
    if (t_f(0).is_empty())
    {
        vec C, NA, CNL, CNR;
        allIdx = PeakSearch(C, NA, CNL, CNR, pss, PSSs, search_full);
        E = calValues(Cal_Energy, C);
        vec dE = calDerivative(Cal_Energy, C);
        res = calValues(Cal_Resolution, E);
        vec FC = res / dE;
        Mult = peakFindMultiplets(C, NA, FC, 2.25);
    }
    //("calPeakSearch_E.txt", E);
    //("calPeakSearch_res.txt", res);
    //("calPeakSearch_pss.txt", pss);
    //("calPeakSearch_Mult.txt", Mult);

    uvec Bool = zeros<uvec>(E.size());
    EF = zeros<vec>(E.size());

    uvec idx;
    // find out for every singlet peak, if it is close to a library peak
    for(int i=0; i<E.size(); ++i) //for i = 1 : length(E)
    {
        if (!Mult(i))
        {
            vec EDev = abs(ELib - E(i));
            double EDevMin = EDev.min();
            if (EDevMin <= EDevMaxGain * res(i))
            {
                // peak is close, therefor mark it
                idx = find(EDev == EDevMin);
                Bool(i) = 1u;
                EF(i) = ELib(idx(0));
            }
        }
    }

    // define steps to reduce PSS level
    double PSS = PSSh;
    double step = (PSSh - PSSl) / 10;
    if (step <= eps(1.0))
    {
        step = 1.0;
    }

    // first look, if we can find at least NPeaksDesired peaks at a PSS level
    // between PSSh and PSSl
    while (PSS >= PSSl)
    {
        idx = find(Bool && pss > PSS);
        if (idx.size() >= NPeaksDesired) break;
        PSS = PSS - step;
    }

    // if failed, loof if we can find at least NPeaksMin peaks at a PSS level
    // between PSSh and PSSl
    if (idx.size() < NPeaksDesired)
    {
        PSS = PSSh;
        while (PSS >= PSSl)
        {
            idx = find(Bool && pss > PSS);
            if (idx.size() >= NPeaksMin) break;
            PSS = PSS - step;
        }
    }

    // reduce vector from all peaks to usable peaks (other entries are 0)
    EF = EF(idx);
    //EF.print("EF = ");
    //idx.print("idx = ");
    return idx;
}

void GammaAnalyALG::calEnergyUpdate1(uvec idx, vec ELib)
{
    // get pat centroids
    vec C = pat.Peaks.col(i_Centroid);
    C = C(idx);

    // fit parameter
    vec p, perr;
    bool s = Independ::calFitPara(p, perr, Cal_Energy, C, ELib);

    //p.print("p = ");
    //perr.print("perr = ");
    if (s) // success, set parameter, data points and flag
    {
        //s = setCalibrationData("CalUpdate1", C, ELib);
        //s = setCalibrationPara("CalUpdate1", p, ep);
        mat data;
        if(lCalDataCheck(data, C, ELib) && Independ::calParaCheck(p))
        {
            vec F0 = Independ::calFcnEval(C, Acal_energy_para[m_curEner](0));
            double sum0 = sum(pow(F0 - ELib, 2)) / C.size();
            vec F1 = Independ::calFcnEval(C, p);
            double sum1 = sum(pow(F1 - ELib, 2)) / C.size();
            if(sum1 < sum0)
            {
                m_curEner = CalUpdate1;
                Cal_Ener_Data[m_curEner] = data;

                field<vec> f_p(2);
                f_p(0) = p; f_p(1) = perr;
                Acal_energy_para[m_curEner] = f_p;
            }
        }
        else { qDebug() << QString("First energy calibration update failed!"); }
    }
    else { qDebug() << QString("First energy calibration update failed!"); }

}

void GammaAnalyALG::calResUpdate(uvec idx)
{
    // Peaks required for fitting
    int NPeaksMin = 4;
    uword NPeaks = idx.size();

    //[NA, C, FC] = dmps(sn, 'NetArea', 'Centroid', 'FWHM_Ch');
    QStringList Opts; Opts << "NetArea" << "Centroid" << "FWHM_Ch";
    field<vec> t_f = dmps(Opts);
    vec NA = t_f(0)(idx);
    vec C = t_f(1)(idx);
    vec FC = t_f(2)(idx);
    //NA.print("NA = ");
    //C.print("C = ");
//    FC.print("FC = ");

    // fit peaks
    uvec Bool = zeros<uvec>(NPeaks);
    uvec tmp;
    vec ret1, ret2, ret3;
    double X2;
    for(uword i=0; i<NPeaks; ++i)
    {
        uword j = idx(i);
        // quick fixed centroid fit, no multiplets!!
        //[bool(i), NA(i), garbage, FC(i), X2] = fitPeakQuick(sn, j, [], j);
        tmp << j;
        Bool(i) = fitPeakQuick(ret1, ret2, ret3, X2, tmp, uvec(), tmp);
        NA(i) = as_scalar(ret1);
        FC(i) = as_scalar(ret3);
    }
    //NA.print("NA = ");
    //FC.print("FC = ");

    // get peak energies
    vec E = calValues(Cal_Energy, C);
    vec dE = calDerivative(Cal_Energy, C);
    //E.print("E = ");
    //dE.print("dE = ");

    // fwhm in energy units
    vec FE = FC % dE;
    //FE.print("FE = ");

    idx = find(Bool);
    //idx.print("idx = ");
    // save data points in any case
    //[stat_sv, msg, name] = setCalibrationData(sn, 'Resolution', 'CalUpdate', E(idx), FE(idx), []);
    /*bool stat_sv = setCalibrationData("CalUpdate", E(idx), FE(idx));
    if (!stat_sv)
    {
        qDebug() << "Can't set resolution calibration update data in calResUpdate";
    }*/

    vec p, perr;
    bool s;
    if (idx.size() > NPeaksMin)
    {
        s = Independ::calFitPara(p, perr, Cal_Resolution, E(idx), FE(idx)); // fit new resolution
        //p.print("p = ");
        //perr.print(" perr = ");
    } else {
        qDebug() << "Not enough peaks for fitting";
        s = false;
    }

    // write results to global
    mat data;
    if (s && lCalDataCheck(data, E(idx), FE(idx)) && Independ::calParaCheck(p) && lCalCheckPerr(p, perr))
    {
        //s = setCalibrationPara("CalUpdate", p, perr);
        //setCurrentCalibration(Cal_Resolution, pname);
        m_curReso = CalResUpdate;
        Cal_Reso_Data[m_curReso] = data;

        field<vec> f_p(2);
        f_p(0) = p; f_p(1) = perr;
        Acal_resolution_para[m_curReso] = f_p;
    }
    else { qDebug() << "Resolution Calibration Update Failed"; }
}

void GammaAnalyALG::calEnergyUpdate2(uvec idx, vec ELib)
{
    uword NPeaks = idx.size();
    vec C = zeros<vec>(NPeaks);

    // fit peaks
    vec ret1, ret2, ret3;
    uvec tmp;
    double X2;
    bool s;
    for(uword i=0; i<NPeaks; ++i) //for i = 1 : NPeaks
    {
        uword j = idx(i);
        // quick fix width fit
        //[s, NA, C(i), FC, X2] = fitPeakQuick(sn, j, j, []);
        tmp << j;
        s = fitPeakQuick(ret1, ret2, ret3, X2, tmp, tmp, uvec());
        C(i) = as_scalar(ret2);
    }
    //C.print("C = ");

    // fit new centroids to library energies
    vec p, perr;
    s = Independ::calFitPara(p, perr, Cal_Energy, C, ELib);
    //p.print("p = ");
    //perr.print("perr = ");
    // write results to global
    if (s)
    {
        //s = setCalibrationData("CalUpdate2", C, ELib);
        //s = setCalibrationPara("CalUpdate2", p, perr);
        mat data;
        if(lCalDataCheck(data, C, ELib) && Independ::calParaCheck(p) && lCalCheckPerr(p, perr))
        {
            vec F02 = Independ::calFcnEval(C, Acal_energy_para[m_curEner](0));
            double sum02 = sum(pow(F02 - ELib, 2)) / C.size();

            vec F2 = Independ::calFcnEval(C, p);
            double sum2 = sum(pow(F2 - ELib, 2)) / C.size();

            if(sum2 < sum02)
            {
                m_curEner = CalUpdate2;
                Cal_Ener_Data[m_curEner] = data;

                field<vec> f_p(2);
                f_p(0) = p; f_p(1) = perr;
                Acal_energy_para[m_curEner] = f_p;
            }
        }
        else { qDebug() << QString("Calibration update failed in calEnergyUpdate2"); }
    }
    else // failed
    {
        qDebug() << "Calibration update failed in calEnergyUpdate2";
    }
}

bool GammaAnalyALG::lCalDataCheck(mat &data, mat x, mat y, mat err)
{
    bool bRet = false;
    uword l = x.size();
    if(y.size() != l)
    {
        return bRet;
    }
    else {
        x = vectorise(x);
        y = vectorise(y);
        if(err.is_empty())
        {
            err = gNaN * ones(l, 1);
        }
        else { err = vectorise(err); }
    }
    data = join_horiz(x, y);
    data.insert_cols(2, err);

    if(!data.col(0).is_finite())
    {
        qDebug() << "Independent variable must be finite";
    }
    else if (!data.col(1).is_finite())
    {
        qDebug() << "Dependent variable must be finite";
    }
    else { bRet = true; }

    return bRet;
}

bool GammaAnalyALG::lCalCheckPerr(mat p, mat &perr)
{
    bool s = true;
    if (perr.is_empty()) perr = gNaN * ones(p.size()-1u, 1u);
    else {
        if (perr.size() != p.size() - 1u)
        {
            qDebug() << "error estimate of wrong size";
            s = false;
        }
        else if (any(vectorise(perr) < 0))
        {
            qDebug() << "negative error estimate";
            s = false;
        }
    }
    return s;
}

/*bool GammaAnalyALG::setCalibrationData(QString src, mat x, mat y, mat err)
{
    bool bRet = true;

    // Check if new data is valid, reshape if necessary
    mat data;
    bRet = lCalDataCheck(data, x, y, err);
    if (!bRet)
    {
        qDebug() << QString("Invalid calibration data from source %1").arg(src);
        return bRet;
    }

    if(Cal_update_data.size() != 3u) Cal_update_data.set_size(3u);
    if(src == "CalUpdate1")
    {
        Cal_update_data(0u) = data;
    }
    else if(src == "CalUpdate")
    {
        Cal_update_data(1u) = data;
    }
    else if(src == "CalUpdate2")
    {
        Cal_update_data(2u) = data;
    }
    else bRet = false;

    return bRet;
}

bool GammaAnalyALG::setCalibrationPara(QString src, mat p, mat perr)
{
    bool bRet = true;
    if(p.is_empty()) return false;

    //[s, msg] = lCalParaCheck(cal, p);
    bRet = Independ::calParaCheck(p);
    if (!bRet) return bRet;

    //[s, msg, perr] = lCalCheckPerr(p, perr);
    bRet = lCalCheckPerr(p, perr);
    if(!bRet) return bRet;

    if(Cal_update_para.size() != 3u) Cal_update_para.set_size(3u);
    if(src == "CalUpdate1")
    {
        Cal_update_para(0u) = vectorise(p);
    }
    else if(src == "CalUpdate")
    {
        Cal_update_para(1u) = vectorise(p);
    }
    else if(src == "CalUpdate2")
    {
        Cal_update_para(2u) = vectorise(p);
    }
    else bRet = false;

    return bRet;
}*/

/*bool GammaAnalyALG::setCurrentCalibration(CalibrationType cal, QString name)
{
    // check only one spectrum is affected
    bool s = true;
    // get index in list of calibration names
    vec sqi, uqi;
    field<vec> f; f << vec("1");
    int n;
    switch (cal)
    {
    case Cal_Energy: n = 1; break;
    case Cal_Resolution:
        patSetPeaks(uvec(0u), QStringList("ResChanged"), f, true);
        n = 2;
        break;
    case Cal_Efficiency: n = 3; break;
    case Cal_Tot_efficiency: n = 4; break;
    case Cal_Tail:
        patSetPeaks(uvec(0u), QStringList("TailChanged"), f, true);
        n = 5;
        break;
    case Cal_Tail_alpha: n = 6; break;
    case Cal_Tail_right: n = 7; break;
    case Cal_Tail_right_alpha: n = 8; break;
    case Cal_Step_ratio:
        patBaseVar(sqi, uqi);
        n = 9;
        break;
    default: s = false; qDebug() << "invalid calibration string in setCurrentCalibration";
    }
    pat.CalibName[n-1] = name;
}*/

bool GammaAnalyALG::fitPeakQuick(vec &Ao, vec &Co, vec &Fo, double &X2, uvec Af, uvec Cf, uvec Ff, QString vb, double FitWidth)
{
    /*% if called with parameter vb, define additional parameters for fitFunction
    if strcmpi(vb, 'verbose')
       addArg = {'Waitbar', 'peaks above linear baseline'};
    else
       addArg = {};
    end */
    bool s = true;
    QStringList addArg("");

    // get index of leftmost and rightmost peak
    //allidx = [Af; Cf; Ff];
    uvec allidx = join_vert(Af, Cf);
    allidx = join_vert(allidx, Ff);
    if (allidx.is_empty())
    {
        s = false;
        Ao.clear();
        Co.clear();
        Fo.clear();
        X2 = gNaN;
        return s;
    }

    uword LIdx = min(allidx);
    uword RIdx = max(allidx);

    // get unique allidx
    uvec Bool = zeros<uvec>(RIdx + 1u);
    Bool(allidx) = ones<uvec>(allidx.n_elem);
    allidx = find(Bool);

    // get peak parameters
    //[C, FC, NA] = dmps(sn, 'Centroid', 'FWHM_Ch', 'NetArea');
    QStringList Opts; Opts << "Centroid" << "FWHM_Ch" << "NetArea";
    field<vec> t_f = dmps(Opts);
    vec C = t_f(0);
    vec FC = t_f(1);
    vec NA = t_f(2);

    // interval for fitting: left centroid - width left FWHM to right centroid + ...
    int CL = floor(C(LIdx) - FitWidth * FC(LIdx));
    int CR = ceil(C(RIdx) + FitWidth * FC(RIdx));

    // spectrum and fitting data points (x,y)
    rowvec S = baseInfo(s_Spectrum); //S = dmspec(sn, 'Spectrum');
    rowvec x = RangeVec2(CL, CR);
    rowvec y = S.cols(CL-1, CR-1);

    // initial baseline y = k*x + d through extremal points (spectrum data)
    double k = (S(CR) - S(CL)) / (CR -CL);
    double d = S(CL) - k * CL;

    // put seed to too small net areas
    double minNA = 0.1;
    if (any(NA(Af) < minNA))
    {
        uvec zeroIdx = find(NA(Af) < minNA);
        MatAssign(NA, Af(zeroIdx), minNA);
    }

    // fit parameters
    //[s, BL, Ao, Co, Fo, X2] = fitFunction(x, y, 'peakBaseLin', [k d], [1 2], NA, Af, C, Cf, FC, Ff, addArg{:});
    field<vec> para;
    field<uvec> free_idx;
    vec tmp; tmp << k << d;
    para << tmp << NA << C << FC;
    free_idx << uvec("0 1") << Af << Cf << Ff;
    s = Independ::fitFunction(para, X2, free_idx, vectorise(x), vectorise(y), "peakBaseLin", addArg);
    if(para.size() > 3)
    {
        Ao = para(1);
        Co = para(2);
        Fo = para(3);
    }

    // write results to temporary PAT
    /*patSetPeaks(allidx, 'FittingMethod', 'Q', 'NetArea', Ao(allidx), ...
       'Centroid', Co(allidx), 'FWHM_Ch', Fo(allidx), 'NetAreaFree', Af, ...
       'FWHMFree', Ff, 'CentroidFree', Cf);*/
    QStringList items;
    items << "NetArea" << "Centroid" << "FWHM_Ch";// << "NetAreaFree" << "FWHMFree" << "CentroidFree";
    field<vec> vals;
    //allidx.print("allidx = ");
    //Ao.print("Ao = ");
    //Co.print("Co = ");
    //Fo.print("Fo = ");
    vals << Ao(allidx) << Co(allidx) << Fo(allidx);
    patSetPeaks(allidx, items, vals);
    vec sqi, uqi;
    patBaseVar(sqi, uqi, allidx);

    // return only fitted parameters
    Ao = Ao(Af);
    Co = Co(Cf);
    Fo = Fo(Ff);
    return s;
}

rowvec GammaAnalyALG::GetFwhmcAll()
{
    rowvec Chan = RangeVec2(1, m_nChans);
    rowvec E = calValues(Cal_Energy, Chan);
    rowvec dE = calDerivative(Cal_Energy, Chan);
    rowvec res = calValues(Cal_Resolution, E);
    rowvec FwhmC = res / dE;
    return FwhmC;
}

stdvec GammaAnalyALG::calculateSCAC(rowvec &spec, rowvec &baseline, rowvec &FwhmC)
{
    rowvec scac;
    int num = spec.size();
    if(num > 0 && baseline.size() == num && FwhmC.size() == num)
    {
        int iSum1,is;
        double s, ds, Sum1, Var1, Var2;
        scac = baseline;
        /* changed by AWST - the ported routine was insp38, while the correct routine is insp70 */
        for (int i = 1; i < num-1; ++i)
        {
            s = (1.25 * FwhmC(i) - 1) / 2 + 0.0000001;
            is = qFloor(s);
            ds = s - is;
            if ( i > is && ( i + is < num -1 ) && s > 0 )
            {
                Sum1 = 0;
                for (iSum1=i-is; iSum1<=i+is; ++iSum1)
                {
                    Sum1 += spec(iSum1);
                }
                Var1 = ds * spec(i - is - 1);
                Var2 = ds * spec(i + is + 1);
                scac(i) = (Sum1 + Var1 + Var2) / (1.25 * FwhmC(i));
            }
        }
    }
    return conv_to<stdvec>::from(scac);
}

stdvec GammaAnalyALG::calculateLC(rowvec &baseline, rowvec &FwhmC, double RiskLevelK)
{
    rowvec lc;
    int num = baseline.size();
    if(num > 0 && FwhmC.size() == num)
    {
        int is;
        double s;
        /*  Initilize the the SCAC array with the Spectrum Array */
        lc = baseline;
        for (int i=1; i<num-1; ++i)
        {
            s = (1.25 * FwhmC(i) - 1) / 2 + 0.00000001;
            is = qFloor(s);
            if (i-is>0 && i+is<num-1)
            {
                if ( baseline(i) > 0 && FwhmC(i) > 0 ) lc(i) = baseline(i) + RiskLevelK * sqrt(baseline(i) / (1.25 * FwhmC(i)));
                else lc(i) = 0;
            }
        }
    }
    return conv_to<stdvec>::from(lc);
}

umat GammaAnalyALG::lGetPlotRanges(int n, const vec &c, const vec &fwhmch, const vec &tlo, const vec &thi/*, vec bhi, vec rdelta*/)
{
    vec fL = 3 + tlo;
    vec fH = 3 + thi;
    /*uvec idx = find_finite(bhi);
    if (!idx.is_empty())
    {
        fH(idx) = fH(idx) + 1.5 / (fwhm(idx) % bhi(idx));
        fL(idx) = fL(idx) + rdelta(idx) / fwhm(idx);
    }*/

    vec l = max(1, floor(c - fL % fwhmch));
    vec h = min(n, ceil(c + fH % fwhmch));

    ivec Bool = zeros<ivec>(n+1);
    //l.print("l = ");
    //h.print("h = ");
    for(uword i=0; i<l.size(); ++i) //for i = 1 : length(l)
    {
        Bool.rows(l(i)-1, h(i)-1) = ones<ivec>(h(i)-l(i)+1);
    }
    //Bool([1, end+1]) = false;
    Bool(0) = 0;
    Bool(n) = 0;

    ivec d = diff(Bool);
    umat r = join_horiz(find(d == 1), find(d == -1));
    return r;
}

uvec GammaAnalyALG::insertpeak(vec iC, PeakInsertMethod src, double minNA)
{
    //uword peakNum = phd->vPeak.size();
    /*vec pC(peakNum), pFC(peakNum), pT(peakNum), pUT(peakNum), pNA(peakNum), pSt(peakNum), pTa(peakNum), pUTa(peakNum);
    vec pBWW(peakNum), pRBC(peakNum), pRDC(peakNum);
    pBWW.fill(gNaN);
    pRBC.fill(gNaN);
    pRDC.fill(gNaN);
    uword t_i = 0;
    foreach (PeakInfo peak, phd->vPeak)
    {
        pC(t_i) = peak.peakCentroid;
        pFC(t_i) = peak.fwhmc;
        pT(t_i) = peak.tail;
        pTa(t_i) = peak.tailAlpha;
        pUT(t_i) = peak.upperTail;
        pUTa(t_i) = peak.upperTailAlpha;
        pNA(t_i) = peak.area;
        pSt(t_i) = peak.area * peak.stepRatio;
        ++t_i;
    }*/

    uvec idx;
    rowvec rg = dmspec(QStringList("AnalysisRange"))(0);

    //accept only peaks in range
    //uvec gidx = find(iC > phd->baseCtrls.rg_low + 1 && iC < phd->baseCtrls.rg_high - 1);
    uvec gidx = find(iC > rg(0)+1 && iC < rg(1)-1);
    iC = iC(gidx);

    // are any peaks to be inserted?
    if(iC.is_empty())
    {
        return idx;
    }

    /*char src_fl, cent_fl;
    switch(src) // peak source flag and centroid source flag
    {
        case Ins_Manual: src_fl = 'I'; cent_fl = 'U'; break;
        case Ins_H0Test:
        case Ins_MultFit: src_fl = 'L'; cent_fl = 'L'; break;
        case Ins_SumPeakModel: src_fl = 'S'; cent_fl = 'L'; break;
        default: qDebug() << "Unkonwn source type"; break;
    }*/

    // get resolution in channels
    /*vec iE = Independ::calFcnEval(Cal_Energy, iC, phd->mapEnerPara[phd->usedEner].p);
    vec idE;
    Independ::calDerivEval(idE, iC, phd->mapEnerPara[phd->usedEner].p);
    vec iR = Independ::calFcnEval(Cal_Resolution, iE, phd->mapResoPara[phd->usedReso].p);*/
    vec iE = calValues(Cal_Energy, iC);
    vec idE = calDerivative(Cal_Energy, iC);
    vec iR = calValues(Cal_Resolution, iE);
    vec iRC = iR / idE;

    // residual above 0
    /*rowvec BaseChan = cumsum(ones<rowvec>(phd->Spec.num_g_channel));
    PiecePoly bpp = Independ::makeBasePP(phd->baseCtrls.XCtrl, phd->baseCtrls.YCtrl, phd->baseCtrls.YSlope, pC, pFC, pSt);
    rowvec approx = Independ::peakBaseSpline(BaseChan, bpp, pNA, pC, pFC, pSt, pT, pTa, pUT, pUTa, pBWW, pRBC, pRDC);
    rowvec resid = spec - approx;*/
    rowvec resid = dmspec(QStringList("Residual"))(0);
    uvec idx_resid = find(resid<0);
    resid(idx_resid) = zeros<rowvec>(idx_resid.size());

    // get net area estimate
    vec iNA(iC.size());
    for(int i=0; i<iC.size(); ++i) //for i=1:length(iC);
    {
        // ch=[floor(iC(i)-iRC(i)): ceil(iC(i)+iRC(i))];
        urowvec ch = RangeVec(floor(iC(i)-iRC(i))-1, ceil(iC(i)+iRC(i))-1);
        //iNA(i) = eval('max(minNA,sum(resid(ch)))','minNA');
        iNA(i) = max(minNA, sum(resid(ch)));
    }

    /*idx.set_size(iC.n_elem);
    for(int i=0; i<iC.size(); ++i)
    {
        int j = 0;
        while (j < peakNum && iC(i) > phd->vPeak[j].peakCentroid) ++j;

        PeakInfo peak;
        peak.peakCentroid = iC(i);
        peak.area = iNA(i);
        peak.sensitivity = gNaN;
        phd->vPeak.insert(phd->vPeak.begin() + j, peak);

        idx(i) = j;
    }*/
    // write peaks to temporary pat
    /*[idx, oidx] = patAddPeaks(sn, 'Centroid', iC, 'NetArea', iNA, ...
      'Source', src_fl, 'CentroidMethod', cent_fl, 'NetAreaMethod', 'S', ...
      'Sensitivity', NaN);*/
    uvec PropIdx;
    PropIdx << i_Centroid << i_NetArea << i_Sensitivity;
    mat Datas(iC.size(), 3);
    Datas.col(0) = iC;
    Datas.col(1) = iNA;
    Datas.col(2) = gNaN * ones<vec>(iC.size());
    uvec oldIdx;
    patAddPeaks(PropIdx, Datas, pat.Peaks, idx, oldIdx);

    return idx;

    // write peaks to temporary pat
    /*[idx, oidx] = patAddPeaks(sn, 'Centroid', iC, 'NetArea', iNA, ...
      'Source', src_fl, 'CentroidMethod', cent_fl, 'NetAreaMethod', 'S', ...
      'Sensitivity', NaN);*/
}

mat GammaAnalyALG::findPeakRange(uvec pnr, vec C, vec F_Ch, double rg_low, double rg_high)
{
    // get all peak centroids and FWHM in channel units
    //[C, F_Ch] = dmps(sn, pat, 'Centroid', 'FWHM_Ch');
    rowvec ar; //ar = dmspec(sn, 'AnalysisRange');
    ar << rg_low << rg_high;

    uword n = pnr.size();
    mat rg = zeros(n, 2);
    vec iL = zeros<vec>(n);
    vec iH = zeros<vec>(n);
    // left and right borders of multiplet
    for(int i=0; i<n; ++i) //for i = 1 : n
    {
        //[rg(i, 1), rg(i, 2), iL(i), iH(i)] = findVPeakRange(pnr(i), C, F_Ch, ar, varargin{:});
        rowvec temp = findVPeakRange(pnr(i), C, F_Ch, ar);
        rg(i, 0) = temp(0);
        rg(i, 1) = temp(1);
        iL(i) = temp(2);
        iH(i) = temp(3);
    }
    return rg;
}

rowvec GammaAnalyALG::findVPeakRange(int pNr, vec C, vec FC, rowvec rg, double LD, double RD)
{
    uword NP = C.size();

    // find leftmost peak in multiplet, and minimum left border
    int i = pNr;
    double l = C(i) - LD * FC(i);
    while (i > 0 && l < C(i-1) + RD * FC(i-1))
    {
        i = i - 1;
        l = min(l, C(i) - LD * FC(i));
    }

    // find rightmost peak in multiplet, and maximum right border
    int j = pNr;
    double h = C(j) + RD * FC(j);
    while (j < NP-1 && h > C(j+1) - LD * FC(j+1))
    {
        j = j + 1;
        h = max(h, C(j) + RD * FC(j));
    }

    // left and right borders of multiplet
    l = max(rg(0), floor(l));
    h = min(rg(1), ceil(h));

    rowvec results;
    results << l << h << i << j;
	return results;
}

/*void GammaAnalyALG::peakfitdlg(QString action, mat ELim)
{
    // write manual changes to PAT
    lUseChanges;

    // spectrum number and index to regarded peaks
    idx = ud.Index;

    // unpack fix/free selections
    if iscell(NAbool), NAbool=[NAbool{:}]; end
    if iscell(Cbool), Cbool=[Cbool{:}]; end
    if iscell(Fbool), Fbool=[Fbool{:}]; end

    iNA = idx(find(~NAbool)); //net areas not fixed
    iCT =idx(find(~Cbool)); //centroid not fixed
    iFC =idx(find(~Fbool)); //fwhm not fixed

    if(action == "Fit")
    {
        fitPeakFull(sn, iNA, iCT, iFC, [], 'Verbose');
    }
    else if(action == "QuickFit")
    {
        fitPeakQuick(sn, iNA, iCT, iFC, 'Verbose');
    }
    lPromptAccept( action);
}*/