AnalysisSystemForRadionuclide_qt_alg/FitFunc.txt

Matlab:
    polyfit.m


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

    // fit parameter
    //[s, msg, p, perr] = calFitPara('Energy', C, ELib, []);
    mat p, perr, ep;
    bool s = calFitPara(p, perr, Cal_Energy, C, ELib);
    if (s) // success, set parameter, data points and flag
    {
        //[s, msg, name] = setCalibrationData(sn, 'Energy', 'CalUpdate1', C, ELib, []);
        QString name, pname;
        s = setCalibrationData(name, Cal_Energy, "CalUpdate1", C, ELib);
        if(s)
        {
            //[s, pname] = setCalibrationPara(sn, 'Energy', 'CalUpdate1', p, [], name);
            s = setCalibrationPara(pname, Cal_Energy, "CalUpdate1", p, ep, name);
            //lPrintSetFail(sn, s, msg);
            if(!s) { qDebug() << "Failed saving energy calibration parameter in calEnergyUpdate1"; }
        }
        else {
            qDebug() << "Failed saving energy calibration data points";
            //[s, pname] = setCalibrationPara(sn, 'Energy', 'CalUpdate1', p, [], '');
            s = setCalibrationPara(pname, Cal_Energy, "CalUpdate1", p, ep, "");
            //lPrintSetFail(sn, s, msg);
            if(!s) { qDebug() << "Failed saving energy calibration parameter in calEnergyUpdate1"; }
        }
        //if (s) { s = setCurrentCalibration(Cal_Energy, pname); }
    }
    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);

    // 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(0u), tmp);
        NA(i) = as_scalar(ret1);
        FC(i) = as_scalar(ret3);
        /*qDebug() << QString("calResUpdate: C: %1 FC: %2 Status: %3 X2: %4")
                    .arg(C(i), 0, 'f', 8)
                    .arg(FC(i), 0, 'f', 8)
                    .arg(Bool(i))
                    .arg(X2, 0, 'f', 8); */
    }

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

    // fwhm in energy units
    vec FE = FC % dE;

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

    mat p, perr;
    bool s;
    if (idx.size() > NPeaksMin)
    {
        // fit new resolution
        //[s, msg, p, perr] = calFitPara('Resolution', E(idx), FE(idx), []);
        s = calFitPara(p, perr, Cal_Resolution, E(idx), FE(idx));
    } else {
        qDebug() << "Not enough peaks for fitting";
        s = false;
    }

    // write results to global
    if (s)
    {
        // success, set parameters
        QString pname;
        if (stat_sv)
        {
            //[s, pname] = setCalibrationPara(sn, 'Resolution', 'CalUpdate', p, perr, name);
            s = setCalibrationPara(pname, Cal_Resolution, "CalUpdate", p, perr, name);
            //lPrintSetFail(sn, s, msg)
        } else {
            //[s, pname] = setCalibrationPara(sn, 'Resolution', 'CalUpdate', p, perr, '');
            s = setCalibrationPara(pname, Cal_Resolution, "CalUpdate", p, perr, "");
            //lPrintSetFail(sn, s, msg)
        }
        if (s)
        {
            //setCurrentCalibration(sn, 'Resolution', pname);
            setCurrentCalibration(Cal_Resolution, pname);
        }
    }
    else { qDebug() << "Resolution Calibration Update Failed"; }
}

void GammaAnalyALG::calEnergyUpdate2(uvec idx, mat 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(0u));
        C(i) = as_scalar(ret2);
        /*qDebug() << QString("calEnergyUpdate2: C: %1 ELib: %2 X2: %3")
                    .arg(C(i), 0, 'f', 8)
                    .arg(ELib(i), 0, 'f', 8)
                    .arg(X2, 0, 'f', 8); */
    }

    // fit new centroids to library energies
    mat p, perr;
    s = calFitPara(p, perr, Cal_Energy, C, ELib);

    // write results to global
    if (s)
    {
        // success, set parameter, data points and flag
        //[s, msg, name] = setCalibrationData(sn, 'Energy', 'CalUpdate2', C, ELib, []);
        QString name, pname;
        s = setCalibrationData(name, Cal_Energy, "CalUpdate2", C, ELib);
        if (s)
        {
            //[s, pname] = setCalibrationPara(sn, 'Energy', 'CalUpdate2', p, [], name);
            s = setCalibrationPara(pname, Cal_Energy, "CalUpdate2", p, perr, name);
            //lPrintSetFail(sn, s, msg);
            if(!s) qDebug() << "Failed saving energy calibration parameter!";
        }
        else {
            qDebug() << "Failed saving energy calibration data points";
            //[s, pname] = setCalibrationPara(sn, 'Energy', 'CalUpdate2', p, [], '');
            s = setCalibrationPara(pname, Cal_Energy, "CalUpdate2", p, perr, "");
            //lPrintSetFail(sn, s, msg);
            if(!s) qDebug() << "Failed saving energy calibration parameter!";
        }
        if (s)
        {
            //setCurrentCalibration(sn, 'Energy', pname);
            setCurrentCalibration(Cal_Energy, pname);
        }
    }
    else // failed
    {
        qDebug() << "Calibration update failed in calEnergyUpdate2";
    }
}

bool GammaAnalyALG::setCalibrationData(QString &name, CalibrationType cal, QString src, mat x, mat y, mat err)
{
    bool bRet = true;

    QMap2Data::iterator it = m_CalData.find(cal);
    QStringList names;
    int num = 1;
    if(it != m_CalData.end())
    {
        names = it.value().keys();
        num = names.filter(src).size() + 1;
    }
    if(src.contains(' ') || num > 99) return false;
    name = QString("%1 %2").arg(src).arg(num);

    // 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;
    }

    // pack and write list back to global
    /*caldatalist(pos, 1) = {name};
    caldatalist(pos, 2) = {data};
    s = setCalDataGlob(sn, cal, caldatalist);*/
    QMap<QString, mat>& subMap = m_CalData[cal];
    subMap[name] = data;

    return bRet;
}

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

    /*[s, msg, calentries] = getCalParaGlob(sn, cal);
    if ~s
        logWrite('%s', msg);
    end

    names = calentries(:, 1);

    [s, msg, name, pos] = lGetNewName(src, names);
    if ~s
       logWrite('invalid %s calibration source for spectrum nr %d: %s', ...
          cal, sn, msg);
       return
    end*/
    QMap2Para::iterator it = m_CalPara.find(cal);
    QStringList names;
    int num = 1;
    if(it != m_CalPara.end())
    {
        names = it.value().keys();
        num = names.filter(src).size() + 1;
    }
    if(src.contains(' ') || num > 99) return false;
    pname = QString("%1 %2").arg(src).arg(num);

    // handle interpolation: internal p also contains x,y data
    if(p.size() == 1 && p(0u) == 1) //if isequal(p, 1)
    {
        bRet = lMakeInterp(p, perr, cal, m_CalData[cal][dnm]);
        if (!bRet)
        {
            return bRet;
        }
    }

    //[s, msg] = lCalParaCheck(cal, p);
    bRet = Independ::calParaCheck(p(0u), p.cols(1u, p.size()-1u));
    if (!bRet) return bRet;

    //[s, msg, perr] = lCalCheckPerr(p, perr);
    if (perr.is_empty()) perr = gNaN * ones(p.size()-1u, 1u);
    else {
        if (perr.size() != p.size() - 1u)
        {
            qDebug() << "error estimate of wrong size";
            return false;
        }
        else if (any(vectorise(perr) < 0))
        {
            qDebug() << "negative error estimate";
            return false;
        }
    }

    /*if isempty(data)
       dnm = '';
    elseif isstr(data)
       [s, msg, x, y, err, dnm] = getCalibrationData(sn, cal, data);
    else
       [s, msg, dnm] = setCalibrationData(sn, cal, src, data);
       if ~s
          logWrite(...
        'failed storing data set, calib [%s], source [%s], sn [%d]: %s', ...
        cal, src, sn, msg);
          return
       end
    end

    [s, msg] = lCalCheckData(sn, cal, dnm);
    if ~s
       logWrite( ...
          'invalid data name, calib [%s], source [%s], sn [%d]: %s', ...
          cal, src, sn, msg);
       return
    end*/

    /*calentries(pos, 1) = {name};
    calentries(pos, 2) = {p};
    calentries(pos, 3) = {perr};
    calentries(pos, 4) = {dnm};
    s = setCalParaGlob(sn, cal, calentries);*/
    CalibPara cp;
    cp.name = pname;
    cp.dnm = dnm;
    cp.p = p;
    cp.perr = perr;
    QMap<QString, CalibPara>& subMap = m_CalPara[cal];
    subMap[pname] = cp;

    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";
    }

    // get old PAT
    /*[s, msg, patEntry] = getPATByName(sn, 'TrueCurrent');
    patName = patEntry{1};
    calData = patEntry{5};
    nat = patEntry{7};
    if ~s
      return
    end

    // missing: check that name is valid calibration string
    calData{n} = name;
    // raise update-flag in NAT
    nat{6} = 1;

    // write change to global
    patEntry([5,7]) = {calData, nat};
    setPATEntry(sn, patEntry);

    // change in temporary as well, if there is a temporary pat
    [s, msg, tpatEntry] = getPATByName(sn, 'Temporary');
    tpatName = tpatEntry{1};
    tcalData = tpatEntry{5};
    tnat = tpatEntry{7};
    if ~s | isempty(tpatName)
      return
    end
    tcalData{n} = name;
    tnat{6} = 1;
    tpatEntry([5,7]) = {tcalData, tnat};
    setPATEntry(sn, tpatEntry);*/
    pat.CalibName[n-1] = name;
}