diff --git a/lib/Physics/CorrelatorFitter.cpp b/lib/Physics/CorrelatorFitter.cpp index c997272..23d9ba9 100644 --- a/lib/Physics/CorrelatorFitter.cpp +++ b/lib/Physics/CorrelatorFitter.cpp @@ -227,6 +227,77 @@ DVec CorrelatorModels::parameterGuess(const DMatSample &corr, return init; } +/****************************************************************************** + * Correlator utilities * + ******************************************************************************/ +DMatSample CorrelatorUtils::shift(const DMatSample &c, const Index ts) +{ + if (ts != 0) + { + const Index nt = c[central].rows(); + DMatSample buf = c; + + FOR_STAT_ARRAY(buf, s) + { + for (Index t = 0; t < nt; ++t) + { + buf[s]((t - ts + nt)%nt) = c[s](t); + } + } + + return buf; + } + else + { + return c; + } +} + +DMatSample CorrelatorUtils::fold(const DMatSample &c) +{ + const Index nt = c[central].rows(); + DMatSample buf = c; + + FOR_STAT_ARRAY(buf, s) + { + for (Index t = 0; t < nt; ++t) + { + buf[s](t) = 0.5*(c[s](t) + c[s]((nt - t) % nt)); + } + } + + return buf; +} + +DMatSample CorrelatorUtils::fourierTransform(const DMatSample &c, FFT &fft, + const unsigned int dir) +{ + const Index nSample = c.size(); + const Index nt = c[central].rows(); + bool isComplex = (c[central].cols() > 1); + CMatSample buf(nSample, nt, 1); + DMatSample out(nSample, nt, 2); + + fft.resize(nt); + FOR_STAT_ARRAY(buf, s) + { + buf[s].real() = c[s].col(0); + if (isComplex) + { + buf[s].imag() = c[s].col(1); + } + else + { + buf[s].imag() = DVec::Constant(nt, 0.); + } + fft(buf[s], dir); + out[s].col(0) = buf[s].real(); + out[s].col(1) = buf[s].imag(); + } + + return out; +} + /****************************************************************************** * CorrelatorFitter implementation * ******************************************************************************/ diff --git a/lib/Physics/CorrelatorFitter.hpp b/lib/Physics/CorrelatorFitter.hpp index 9acadb7..694fd9b 100644 --- a/lib/Physics/CorrelatorFitter.hpp +++ b/lib/Physics/CorrelatorFitter.hpp @@ -22,6 +22,7 @@ #include #include +#include #include BEGIN_LATAN_NAMESPACE @@ -49,6 +50,17 @@ namespace CorrelatorModels DVec parameterGuess(const DMatSample &corr, const ModelPar par); }; +/****************************************************************************** + * Correlator utilities * + ******************************************************************************/ +namespace CorrelatorUtils +{ + DMatSample shift(const DMatSample &c, const Index ts); + DMatSample fold(const DMatSample &c); + DMatSample fourierTransform(const DMatSample &c, FFT &fft, + const unsigned int dir = FFT::Forward); +}; + /****************************************************************************** * Correlator fit utility class * ******************************************************************************/ diff --git a/physics/2pt-fit.cpp b/physics/2pt-fit.cpp index 4e86ee1..c322ada 100644 --- a/physics/2pt-fit.cpp +++ b/physics/2pt-fit.cpp @@ -105,28 +105,14 @@ int main(int argc, char *argv[]) DMatSample tmp, corr; Index nSample, nt; - tmp = Io::load(corrFileName); - nSample = tmp.size(); - nt = tmp[central].rows(); - tmp = tmp.block(0, 0, nt, 1); - corr = tmp; - FOR_STAT_ARRAY(corr, s) - { - for (Index t = 0; t < nt; ++t) - { - corr[s]((t - shift + nt)%nt) = tmp[s](t); - } - } + corr = Io::load(corrFileName); + nSample = corr.size(); + nt = corr[central].rows(); + corr = corr.block(0, 0, nt, 1); + corr = CorrelatorUtils::shift(corr, shift); if (fold) { - tmp = corr; - FOR_STAT_ARRAY(corr, s) - { - for (Index t = 0; t < nt; ++t) - { - corr[s](t) = 0.5*(tmp[s](t) + tmp[s]((nt - t) % nt)); - } - } + corr = CorrelatorUtils::fold(corr); } // make model //////////////////////////////////////////////////////////////