/***************************************************************************
		Copyright (C) 2002-2015 Kentaro Kitagawa
		                   kitag@issp.u-tokyo.ac.jp
		
		This program is free software; you can redistribute it and/or
		modify it under the terms of the GNU Library General Public
		License as published by the Free Software Foundation; either
		version 2 of the License, or (at your option) any later version.
		
		You should have received a copy of the GNU Library General 
		Public License and a list of authors along with this program; 
		see the files COPYING and AUTHORS.
***************************************************************************/
//---------------------------------------------------------------------------
#include "nmrspectrumbase.h"
#include "nmrpulse.h"

#include <graph.h>
#include <graphwidget.h>
#include <xwavengraph.h>

#include <QPushButton>
#include <QComboBox>
#include <QCheckBox>

//---------------------------------------------------------------------------
template <class FRM>
XNMRSpectrumBase<FRM>::XNMRSpectrumBase(const char *name, bool runtime,
	Transaction &tr_meas, const shared_ptr<XMeasure> &meas)
	: XSecondaryDriver(name, runtime, ref(tr_meas), meas),
	m_pulse(create<XItemNode<XDriverList, XNMRPulseAnalyzer> >(
	  "PulseAnalyzer", false, ref(tr_meas), meas->drivers(), true)),
	m_bandWidth(create<XDoubleNode>("BandWidth", false)),
	m_bwList(create<XComboNode>("BandWidthList", false, true)),
	m_autoPhase(create<XBoolNode>("AutoPhase", false)),
	m_phase(create<XDoubleNode>("Phase", false, "%.2f")),
	m_clear(create<XTouchableNode>("Clear", true)),
	m_solverList(create<XComboNode>("SpectrumSolver", false, true)),
	m_windowFunc(create<XComboNode>("WindowFunc", false, true)),
	m_windowWidth(create<XDoubleNode>("WindowWidth", false)),
	m_solver(create<SpectrumSolverWrapper>("SpectrumSolver", true, m_solverList, m_windowFunc, m_windowWidth)),
    m_form(new FRM),
	m_statusPrinter(XStatusPrinter::create(m_form.get())),
    m_spectrum(create<XWaveNGraph>("Spectrum", true, m_form->m_graph, m_form->m_edDump, m_form->m_tbDump, m_form->m_btnDump,
            m_form->m_tlbMath, meas, static_pointer_cast<XDriver>(shared_from_this()))) {
    m_form->m_btnClear->setIcon(QApplication::style()->standardIcon(QStyle::SP_DialogResetButton));
    
	connect(pulse());

	iterate_commit([=](Transaction &tr){
		const char *labels[] = {"X", "Re [V]", "Im [V]", "Weights", "Abs [V]", "Dark [V]"};
		tr[ *m_spectrum].setColCount(6, labels);
        if( !tr[ *m_spectrum].insertPlot(tr, labels[4], 0, 4, -1, 3)) return;
        if( !tr[ *m_spectrum].insertPlot(tr, labels[1], 0, 1, -1, 3)) return;
        if( !tr[ *m_spectrum].insertPlot(tr, labels[2], 0, 2, -1, 3)) return;
        if( !tr[ *m_spectrum].insertPlot(tr, labels[5], 0, 5, -1, 3)) return;
		tr[ *tr[ *m_spectrum].axisy()->label()] = i18n("Intens. [V]");
		tr[ *tr[ *m_spectrum].plot(1)->label()] = i18n("real part");
        tr[ *tr[ *m_spectrum].plot(3)->lineColor()] = clLime; //QColor(0xa0, 0xa0, 0x00).rgb();
        tr[ *tr[ *m_spectrum].plot(2)->lineColor()] = (unsigned int)tr[ *tr[ *m_spectrum].plot(1)->lineColor()];
        tr[ *tr[ *m_spectrum].plot(1)->lineColor()] = (unsigned int)tr[ *tr[ *m_spectrum].plot(0)->lineColor()];
        tr[ *tr[ *m_spectrum].plot(0)->lineColor()] = clWhite;
        tr[ *tr[ *m_spectrum].plot(0)->barColor()] = QColor(0xf0, 0xf0, 0xc0).rgb();
//		tr[ *tr[ *m_spectrum].plot(0)->barColor()] = QColor(0x60, 0x60, 0xc0).rgb();
//		tr[ *tr[ *m_spectrum].plot(0)->lineColor()] = QColor(0x60, 0x60, 0xc0).rgb();
        tr[ *tr[ *m_spectrum].plot(1)->drawPoints()] = false;
		tr[ *tr[ *m_spectrum].plot(2)->label()] = i18n("imag. part");
		tr[ *tr[ *m_spectrum].plot(2)->drawPoints()] = false;
		tr[ *tr[ *m_spectrum].plot(0)->label()] = i18n("abs.");
		tr[ *tr[ *m_spectrum].plot(0)->drawPoints()] = false;
		tr[ *tr[ *m_spectrum].plot(0)->drawLines()] = true;
		tr[ *tr[ *m_spectrum].plot(0)->drawBars()] = true;
		tr[ *tr[ *m_spectrum].plot(0)->intensity()] = 0.5;
		tr[ *tr[ *m_spectrum].plot(3)->label()] = i18n("dark");
		tr[ *tr[ *m_spectrum].plot(3)->drawBars()] = false;
		tr[ *tr[ *m_spectrum].plot(3)->drawLines()] = true;
		tr[ *tr[ *m_spectrum].plot(3)->drawPoints()] = false;
		tr[ *tr[ *m_spectrum].plot(3)->intensity()] = 0.5;
		{
			shared_ptr<XXYPlot> plot = m_spectrum->graph()->plots()->template create<XXYPlot>(
				tr, "Peaks", true, tr, m_spectrum->graph());
            if( !plot) return;
			m_peakPlot = plot;
			tr[ *plot->label()] = i18n("Peaks");
			tr[ *plot->axisX()] = tr[ *m_spectrum].axisx();
			tr[ *plot->axisY()] = tr[ *m_spectrum].axisy();
			tr[ *plot->drawPoints()] = false;
			tr[ *plot->drawLines()] = false;
			tr[ *plot->drawBars()] = true;
            tr[ *plot->intensity()] = 0.3;
			tr[ *plot->displayMajorGrid()] = false;
            tr[ *plot->pointColor()] = clWhite; //QColor(0xa0, 0x00, 0xa0).rgb();
            tr[ *plot->barColor()] = clWhite; //QColor(0xa0, 0x00, 0xa0).rgb();
			tr[ *plot->clearPoints()].setUIEnabled(false);
			tr[ *plot->maxCount()].setUIEnabled(false);
		}
		tr[ *m_spectrum].clearPoints();

		tr[ *bandWidth()] = 50;
		tr[ *bwList()].add("50%");
		tr[ *bwList()].add("100%");
		tr[ *bwList()].add("200%");
		tr[ *bwList()] = 1;
		tr[ *autoPhase()] = true;

		tr[ *windowFunc()].str(XString(SpectrumSolverWrapper::WINDOW_FUNC_DEFAULT));
		tr[ *windowWidth()] = 100.0;
    });
  
    //Ranges should be preset in prior to connectors.
    m_form->m_dblPhase->setRange(-360.0, 360.0);
    m_form->m_dblPhase->setSingleStep(10.0);
    m_form->m_dblWindowWidth->setRange(0.1, 200.0);
    m_form->m_dblWindowWidth->setSingleStep(1.0);

    m_conBaseUIs = {
        xqcon_create<XQLineEditConnector>(m_bandWidth, m_form->m_edBW),
        xqcon_create<XQComboBoxConnector>(m_bwList, m_form->m_cmbBWList, Snapshot( *m_bwList)),
        xqcon_create<XQDoubleSpinBoxConnector>(m_phase, m_form->m_dblPhase, m_form->m_slPhase),
        xqcon_create<XQToggleButtonConnector>(m_autoPhase, m_form->m_ckbAutoPhase),
        xqcon_create<XQComboBoxConnector>(m_pulse, m_form->m_cmbPulse, ref(tr_meas)),
        xqcon_create<XQButtonConnector>(m_clear, m_form->m_btnClear),
        xqcon_create<XQComboBoxConnector>(m_solverList, m_form->m_cmbSolver, Snapshot( *m_solverList)),
        xqcon_create<XQDoubleSpinBoxConnector>(m_windowWidth, m_form->m_dblWindowWidth, m_form->m_slWindowWidth),
        xqcon_create<XQComboBoxConnector>(m_windowFunc, m_form->m_cmbWindowFunc, Snapshot( *m_windowFunc)),
    };

	iterate_commit([=](Transaction &tr){
		m_lsnOnClear = tr[ *m_clear].onTouch().connectWeakly(
			shared_from_this(), &XNMRSpectrumBase<FRM>::onClear);
		m_lsnOnCondChanged = tr[ *bandWidth()].onValueChanged().connectWeakly(
			shared_from_this(), &XNMRSpectrumBase<FRM>::onCondChanged);
		tr[ *autoPhase()].onValueChanged().connect(m_lsnOnCondChanged);
		tr[ *phase()].onValueChanged().connect(m_lsnOnCondChanged);
		tr[ *solverList()].onValueChanged().connect(m_lsnOnCondChanged);
		tr[ *windowWidth()].onValueChanged().connect(m_lsnOnCondChanged);
		tr[ *windowFunc()].onValueChanged().connect(m_lsnOnCondChanged);
		tr[ *bwList()].onValueChanged().connect(m_lsnOnCondChanged);
    });
}
template <class FRM>
XNMRSpectrumBase<FRM>::~XNMRSpectrumBase() {
}
template <class FRM>
void
XNMRSpectrumBase<FRM>::showForms() {
    m_form->showNormal();
	m_form->raise();
}
template <class FRM>
void
XNMRSpectrumBase<FRM>::onCondChanged(const Snapshot &shot, XValueNodeBase *node) {
//    if((node == phase()) && *autoPhase()) return;
	if((node == bandWidth().get()) || onCondChangedImpl(shot, node))
        trans( *this).m_timeClearRequested = XTime::now();
    requestAnalysis();
}
template <class FRM>
void
XNMRSpectrumBase<FRM>::onClear(const Snapshot &shot, XTouchableNode *) {
    trans( *this).m_timeClearRequested = XTime::now();
    requestAnalysis();
}
template <class FRM>
bool
XNMRSpectrumBase<FRM>::checkDependency(const Snapshot &shot_this,
	const Snapshot &shot_emitter, const Snapshot &shot_others,
	XDriver *emitter) const {
    shared_ptr<XNMRPulseAnalyzer> pulse__ = shot_this[ *pulse()];
    if( !pulse__) return false;
    if(emitter == this) return true;
    return (emitter == pulse__.get()) && checkDependencyImpl(shot_this, shot_emitter, shot_others, emitter);
}
template <class FRM>
void
XNMRSpectrumBase<FRM>::analyze(Transaction &tr, const Snapshot &shot_emitter, const Snapshot &shot_others,
	XDriver *emitter) {
	const Snapshot &shot_this(tr);

	shared_ptr<XNMRPulseAnalyzer> pulse__ = shot_this[ *pulse()];
	assert( pulse__ );
	const Snapshot &shot_pulse((emitter == pulse__.get()) ? shot_emitter : shot_others);
 
	if(shot_pulse[ *pulse__->exAvgIncr()]) {
		m_statusPrinter->printWarning(i18n("Do NOT use incremental avg. Skipping."));
		throw XSkippedRecordError(__FILE__, __LINE__);
	}

    bool clear = (shot_this[ *this].m_timeClearRequested.isSet());
    tr[ *this].m_timeClearRequested = {};
  
//	double interval = shot_pulse[ *pulse__].interval();
	double df = shot_pulse[ *pulse__].dFreq();
	
	double res = getFreqResHint(shot_this);
	res = df * std::max(1L, lrint(res / df - 0.5));

	double max__ = getMaxFreq(shot_this);
	double min__ = getMinFreq(shot_this);
	
	if(max__ <= min__) {
		throw XSkippedRecordError(i18n("Invalid min. and max."), __FILE__, __LINE__);
	}
    constexpr ssize_t MAX_ACCUM_LEN = 65536 * 2;
    if(res * MAX_ACCUM_LEN < max__ - min__) {
        //restricts a size of accumulation buffers, due to memory consumption and friendliness.
        res = df * lrint((max__ - min__) / MAX_ACCUM_LEN / df - 0.5);
//		throw XSkippedRecordError(i18n("Too small resolution."), __FILE__, __LINE__);
    }
    //small change is discarded.
	if(fabs(log(shot_this[ *this].res() / res)) < log(2.0))
		res = shot_this[ *this].res();
	if((shot_this[ *this].res() != res) || clear) {
		tr[ *this].m_res = res;
		for(int bank = 0; bank < Payload::ACCUM_BANKS; bank++) {
			tr[ *this].m_accum[bank].clear();
			tr[ *this].m_accum_weights[bank].clear();
			tr[ *this].m_accum_dark[bank].clear();
		}
	}
	else {
        //expands/shrinks the begining of buffers.
        int diff = lrint(shot_this[ *this].min() / res) - lrint(min__ / res);
		for(int bank = 0; bank < Payload::ACCUM_BANKS; bank++) {
            auto &accum = tr[ *this].m_accum[bank];
            auto &accum_weights = tr[ *this].m_accum_weights[bank];
            auto &accum_dark = tr[ *this].m_accum_dark[bank];
            for(int i = 0; i < diff; i++) {
                accum.push_front(0.0);
                accum_weights.push_front(0);
                accum_dark.push_front(0.0);
			}
			for(int i = 0; i < -diff; i++) {
                if( !accum.empty()) {
                    accum.pop_front();
                    accum_weights.pop_front();
                    accum_dark.pop_front();
				}
			}
		}
	}
	tr[ *this].m_min = min__;
    //expands/shrinks the end of buffers.
    int length = lrint((max__ - min__) / res);
	for(int bank = 0; bank < Payload::ACCUM_BANKS; bank++) {
		tr[ *this].m_accum[bank].resize(length, 0.0);
		tr[ *this].m_accum_weights[bank].resize(length, 0);
		tr[ *this].m_accum_dark[bank].resize(length, 0.0);
	}
	tr[ *this].m_wave.resize(length);
	std::fill(tr[ *this].m_wave.begin(), tr[ *this].m_wave.end(), std::complex<double>(0.0));
	tr[ *this].m_weights.resize(length);
	std::fill(tr[ *this].m_weights.begin(), tr[ *this].m_weights.end(), 0.0);
	tr[ *this].m_darkPSD.resize(length);
	std::fill(tr[ *this].m_darkPSD.begin(), tr[ *this].m_darkPSD.end(), 0.0);

	if(clear) {
		tr[ *m_spectrum].clearPoints();
		tr[ *this].m_peaks.clear();
		trans( *pulse__->avgClear()).touch();
		throw XSkippedRecordError(__FILE__, __LINE__);
	}

	if(emitter == pulse__.get()) {
		fssum(tr, shot_pulse, shot_others);
		m_isInstrumControlRequested = true;
	}
	else
		m_isInstrumControlRequested = false;
	
	analyzeIFT(tr, shot_pulse);
	std::vector<std::complex<double> > &wave(tr[ *this].m_wave);
	const std::vector<double> &weights(shot_this[ *this].weights());
	int wave_size = shot_this[ *this].wave().size();
	if(shot_this[ *autoPhase()]) {
		std::complex<double> csum(0.0, 0.0);
		for(unsigned int i = 0; i < wave_size; i++)
			csum += wave[i] * weights[i];
		double ph = 180.0 / M_PI * atan2(std::imag(csum), std::real(csum));
		if(fabs(ph) < 180.0)
			tr[ *phase()] = ph;
		tr.unmark(m_lsnOnCondChanged); //avoiding recursive signaling.
	}
	double ph = shot_this[ *phase()] / 180.0 * M_PI;
	std::complex<double> cph = std::polar(1.0, -ph);
	for(unsigned int i = 0; i < wave_size; i++)
		wave[i] *= cph;
}
template <class FRM>
void
XNMRSpectrumBase<FRM>::visualize(const Snapshot &shot) {
	if( !shot[ *this].time()) {
		iterate_commit([=](Transaction &tr){
			tr[ *m_spectrum].clearPoints();
			tr[ *m_peakPlot->maxCount()] = 0;
        });
		return;
	}

    if(m_isInstrumControlRequested.compare_set_strong((int)true, (int)false))
		rearrangeInstrum(shot);

	int length = shot[ *this].wave().size();
	std::vector<double> values;
	getValues(shot, values);
	assert(values.size() == length);
    m_spectrum->iterate_commit([=](Transaction &tr){
		double th = FFT::windowFuncHamming(0.1);
		const std::complex<double> *wave( &shot[ *this].wave()[0]);
		const double *weights( &shot[ *this].weights()[0]);
		const double *darkpsd( &shot[ *this].darkPSD()[0]);
		tr[ *m_spectrum].setRowCount(length);
        std::vector<double> colx(length);
        std::vector<float> colr(length), coli(length), colw(length),
            colabs(length), coldark(length);
		for(int i = 0; i < length; i++) {
			colx[i] = values[i];
			colr[i] = std::real(wave[i]);
			coli[i] = std::imag(wave[i]);
			colw[i] = (weights[i] > th) ? weights[i] : 0.0;
			colabs[i] = std::abs(wave[i]);
			coldark[i] = sqrt(darkpsd[i]);
		}
        tr[ *m_spectrum].setColumn(0, std::move(colx), 9);
        tr[ *m_spectrum].setColumn(1, std::move(colr), 5);
        tr[ *m_spectrum].setColumn(2, std::move(coli), 5);
        tr[ *m_spectrum].setColumn(3, std::move(colw), 4);
        tr[ *m_spectrum].setColumn(4, std::move(colabs), 5);
        tr[ *m_spectrum].setColumn(5, std::move(coldark), 4);
        const auto &peaks(shot[ *this].m_peaks);
		int peaks_size = peaks.size();
		tr[ *m_peakPlot->maxCount()] = peaks_size;
        auto &points(tr[ *m_peakPlot].points());
		points.resize(peaks_size);
		for(int i = 0; i < peaks_size; i++) {
			double x = peaks[i].second;
			int j = lrint(x - 0.5);
			j = std::min(std::max(0, j), length - 2);
			double a = values[j] + (values[j + 1] - values[j]) * (x - j);
			points[i] = XGraph::ValPoint(a, peaks[i].first);
		}
		m_spectrum->drawGraph(tr);
    });
}

template <class FRM>
void
XNMRSpectrumBase<FRM>::fssum(Transaction &tr, const Snapshot &shot_pulse, const Snapshot &shot_others) {
	const Snapshot &shot_this(tr);
	shared_ptr<XNMRPulseAnalyzer> pulse__ = shot_this[ *pulse()];

	int len = shot_pulse[ *pulse__].ftWidth();
	double df = shot_pulse[ *pulse__].dFreq();
	if((len == 0) || (df == 0)) {
		throw XRecordError(i18n("Invalid waveform."), __FILE__, __LINE__);
	}
	//bw *= 1.8; // for Hamming.
	//	bw *= 3.6; // for FlatTop.
	//bw *= 2.2; // for Kaiser3.
	int bw = abs(lrint(shot_this[ *bandWidth()] * 1000.0 / df * 2.2));
	if(bw >= len) {
		throw XRecordError(i18n("BW beyond Nyquist freq."), __FILE__, __LINE__);
	}
	double cfreq = getCurrentCenterFreq(shot_this, shot_others);
	std::vector<std::complex<double> > ftwavein(len, 0.0), ftwaveout(len);
	if( !shot_this[ *this].m_preFFT || (shot_this[ *this].m_preFFT->length() != len)) {
		tr[ *this].m_preFFT.reset(new FFT(-1, len));
	}
	int wlen = std::min(len, (int)shot_pulse[ *pulse__].wave().size());
	int woff = -shot_pulse[ *pulse__].waveFTPos()
		+ len * ((shot_pulse[ *pulse__].waveFTPos() > 0) ? ((int)shot_pulse[ *pulse__].waveFTPos() / len + 1) : 0);
	const std::complex<double> *pulse_wave( &shot_pulse[ *pulse__].wave()[0]);
	for(int i = 0; i < wlen; i++) {
		int j = (i + woff) % len;
		ftwavein[j] = pulse_wave[i];
	}
	tr[ *this].m_preFFT->exec(ftwavein, ftwaveout);
	bw /= 2.0;
	double normalize = 1.0 / (double)shot_pulse[ *pulse__].wave().size();
	double darknormalize = shot_this[ *this].res() / df;
	for(int bank = 0; bank < Payload::ACCUM_BANKS; bank++) {
		double min = shot_this[ *this].min();
		double res = shot_this[ *this].res();
		int size = (int)shot_this[ *this].m_accum[bank].size();
		std::deque<std::complex<double> > &accum_wave(tr[ *this].m_accum[bank]);
		std::deque<double> &accum_weights(tr[ *this].m_accum_weights[bank]);
		std::deque<double> &accum_dark(tr[ *this].m_accum_dark[bank]);
		const double *pulse_dark( &shot_pulse[ *pulse__].darkPSD()[0]);
		for(int i = -bw / 2; i <= bw / 2; i++) {
			double freq = i * df;
			int idx = lrint((cfreq + freq - min) / res);
			if((idx >= size) || (idx < 0))
				continue;
			double w = FFT::windowFuncKaiser1((double)i / bw);
			int j = (i + len) % len;
			accum_wave[idx] += ftwaveout[j] * w * normalize;
			accum_weights[idx] += w;
			accum_dark[idx] += pulse_dark[j] * w * w * darknormalize;
		}
		bw *= 2.0;
	}
}
template <class FRM>
void
XNMRSpectrumBase<FRM>::analyzeIFT(Transaction &tr, const Snapshot &shot_pulse) {
	const Snapshot &shot_this(tr);
	int bank = shot_this[ *bwList()];
	if((bank < 0) || (bank >= Payload::ACCUM_BANKS))
		throw XSkippedRecordError(__FILE__, __LINE__);
	double bw_coeff = 0.5 * pow(2.0, (double)bank);
	
	double th = FFT::windowFuncHamming(0.49);
	int max_idx = 0;
	const std::deque<std::complex<double> > &accum_wave(shot_this[ *this].m_accum[bank]);
	const std::deque<double> &accum_weights(shot_this[ *this].m_accum_weights[bank]);
	const std::deque<double> &accum_dark(shot_this[ *this].m_accum_dark[bank]);
	int accum_size = accum_wave.size();
	int min_idx = accum_size - 1;
	int taps_max = 0; 
	for(int i = 0; i < accum_size; i++) {
		if(accum_weights[i] > th) {
			min_idx = std::min(min_idx, i);
			max_idx = std::max(max_idx, i);
			taps_max++;
		}
	}
	if(max_idx <= min_idx)
		throw XSkippedRecordError(__FILE__, __LINE__);
	shared_ptr<XNMRPulseAnalyzer> pulse__ = shot_this[ *pulse()];
	double res = shot_this[ *this].res();
	int iftlen = max_idx - min_idx + 1;
	double wave_period = shot_pulse[ *pulse__].waveWidth() * shot_pulse[ *pulse__].interval();
	int npad = lrint(
		6.0 / (res * wave_period) + 0.5); //# of pads in frequency domain.
	//Truncation factor for IFFT.
	int trunc2 = lrint(pow(2.0, ceil(log(iftlen * 0.03) / log(2.0))));
	if(trunc2 < 1)
		throw XSkippedRecordError(__FILE__, __LINE__);
	iftlen = ((iftlen * 3 / 2 + npad) / trunc2 + 1) * trunc2;
	int tdsize = lrint(wave_period * res * iftlen);
	int iftorigin = lrint(shot_pulse[ *pulse__].waveFTPos() * shot_pulse[ *pulse__].interval() * res * iftlen);
	int bwinv = abs(lrint(1.0 / (shot_this[ *bandWidth()] * bw_coeff * 1000.0 * shot_pulse[ *pulse__].interval() * res * iftlen)));
	
	if( !shot_this[ *this].m_ift || (shot_this[ *this].m_ift->length() != iftlen)) {
		tr[ *this].m_ift.reset(new FFT(1, iftlen));
	}
	
	std::vector<std::complex<double> > fftwave(iftlen), iftwave(iftlen);
	std::fill(fftwave.begin(), fftwave.end(), std::complex<double>(0.0));
	for(int i = min_idx; i <= max_idx; i++) {
        int k = (i - (max_idx + min_idx) / 2 + iftlen) % iftlen;
        assert(k >= 0);
		if(accum_weights[i] > th)
			fftwave[k] = accum_wave[i] / accum_weights[i];
	}
	tr[ *this].m_ift->exec(fftwave, iftwave);
	
	SpectrumSolver &solver(tr[ *m_solver].solver());
	std::vector<std::complex<double> > solverin;
	FFT::twindowfunc wndfunc = m_solver->windowFunc(shot_this);
	double wndwidth = shot_this[ *windowWidth()] / 100.0;
	double psdcoeff = 1.0;
	if(solver.isFT()) {
		std::vector<double> weight;
		SpectrumSolver::window(tdsize, -iftorigin, wndfunc, wndwidth, weight);
		double w = 0;
		for(int i = 0; i < tdsize; i++)
			w += weight[i] * weight[i];
		psdcoeff = w / (double)tdsize;
		//Compensate broadening due to convolution.
		solverin.resize(iftlen);
		double wlen = SpectrumSolver::windowLength(tdsize, -iftorigin, wndwidth);
		wlen += bwinv * 2; //effect of convolution.
		wndwidth = wlen / solverin.size();
		iftorigin = solverin.size() / 2;
	}
	else {
		solverin.resize(tdsize);
	}
	for(int i = 0; i < (int)solverin.size(); i++) {
		int k = (-iftorigin + i + iftlen) % iftlen;
		assert(k >= 0);
		solverin[i] = iftwave[k];
	}
	try {
		solver.exec(solverin, fftwave, -iftorigin, 0.1e-2, wndfunc, wndwidth);
	}
	catch (XKameError &e) {
		throw XSkippedRecordError(e.msg(), __FILE__, __LINE__);
	}

	std::vector<std::complex<double> > &wave(tr[ *this].m_wave);
	std::vector<double> &weights(tr[ *this].m_weights);
	std::vector<double> &darkpsd(tr[ *this].m_darkPSD);
	psdcoeff /= wave_period;
	for(int i = min_idx; i <= max_idx; i++) {
		int k = (i - (max_idx + min_idx) / 2 + iftlen) % iftlen;
        assert(k >= 0);
		wave[i] = fftwave[k] / (double)iftlen;
		double w = accum_weights[i];
		weights[i] = w;
		darkpsd[i] = accum_dark[i] / (w * w) * psdcoeff;
	}
	th = FFT::windowFuncHamming(0.1);
	tr[ *this].m_peaks.clear();
	int weights_size = shot_this[ *this].weights().size();
	std::deque<std::pair<double, double> > &peaks(tr[ *this].m_peaks);
	for(int i = 0; i < solver.peaks().size(); i++) {
		double k = solver.peaks()[i].second;
		double j = (k > iftlen / 2) ? (k - iftlen) : k;
		j += (max_idx + min_idx) / 2;
		int l = lrint(j);
		if((l >= 0) && (l < weights_size) && (weights[l] > th))
            peaks.emplace_back(solver.peaks()[i].first / (double)iftlen, j);
	}
}
