qavrgfitter.cpp

Go to the documentation of this file.

#include "qavrgfitter.h"
#include "qavrgacquisition.h"
#include "qavrgsettings.h"
#include "qavrgmatrix.h"

#include <QTimer>
#include <QThread>
#include <QMetaProperty>
#include <stdio.h>

QavrgFitter::QavrgFitter(QavrgAcquisition *acq, int chan)
  : QObject(),
    m_Mutex(QMutex::Recursive),
    m_Acquisition(acq),
    m_Channel(chan),
    m_NSamples(0),
    m_InputGain(m_Acquisition->saver(), this, "inputGain", 0, "Channel Input Gain"),
    m_InputOffset(m_Acquisition->saver(), this, "inputOffset", 0, "Channel Input Offset"),
    m_InputBandwidth(m_Acquisition->saver(), this, "inputBandwidth", 0, "Channel Input Bandwidth"),
    m_InputCoupling(m_Acquisition->saver(), this, "inputCoupling", 3, "Channel Input Coupling"),
    m_FittingOffset(m_Acquisition->saver(), this, "fittingOffset", 0, "Channel Fitting Offset"),
    m_FittingStart(m_Acquisition->saver(), this, "fittingStart", 0, "Channel Fitting Start"),
    m_FittingEnd(m_Acquisition->saver(), this, "fittingEnd", 32, "Channel Fitting End"),
    m_DisplayData(m_Acquisition->saver(), this, "displayData", false, "Display Channel Data"),
    m_DisplayRaw(m_Acquisition->saver(), this, "displayRaw", false, "Display Channel Raw Data"),
    m_DisplayReference(m_Acquisition->saver(), this, "displayReference", false, "Display Channel Reference Data"),
    m_DisplayDark(m_Acquisition->saver(), this, "displayDark", false, "Display Channel Dark Data"),
    m_DisplayFit(m_Acquisition->saver(), this, "displayFit", false, "Display Channel Fit Data"),
    m_DarkAvailable(QcepSettingsSaverWPtr(), this, "darkAvailable", false, "Channel Dark Data Available"),
    m_ReferenceAvailable(QcepSettingsSaverWPtr(), this, "referenceAvailable", false, "Channel Reference Data Available")
//    m_DarkChanged(this, "darkChanged", false),
//    m_ReferenceChanged(this, "referenceChanged", false),
{
//  printf("Fitter channel %d started\n", chan);

  prepareCalculation();
}

QavrgFitter::~QavrgFitter()
{
//  printf("Fitter channel %d stopped\n", m_Channel);
}

int QavrgFitter::resultSize()
{
  QMutexLocker lock(&m_Mutex);

  prepareCalculation();

  return m_ResultSize;
}

QVector<double> QavrgFitter::readResult(int parm, int start, int nbins)
{
  QMutexLocker lock(&m_Mutex);

  return m_Results.value(parm).mid(start,nbins);
}

QVector<double> QavrgFitter::readResult(int parm)
{
  QMutexLocker lock(&m_Mutex);

  return m_Results.value(parm);
}

double QavrgFitter::readResult(int parm, int bin)
{
  QMutexLocker lock(&m_Mutex);

  return m_Results.value(parm).value(bin);
}

double QavrgFitter::readResultAverage(int parm, int start, int nbins)
{
  QMutexLocker lock(&m_Mutex);

  double res = 0;
  QVector<double> vec = m_Results.value(parm);

  for (int i=0; i<nbins; i++) {
    res += vec.value(i+start);
  }

  return res / nbins;
}

double QavrgFitter::readResultBunchAverage(int parm, int bunch, int norbits)
{
  QMutexLocker lock(&m_Mutex);

  double res = 0;
  QVector<double> vec = m_Results.value(parm);

  for (int i=0; i<norbits; i++) {
    int nb = i*m_Acquisition -> get_FilledBucketsPerOrbit() + bunch;

    res += vec.value(nb);
  }

  return res/norbits;
}

void QavrgFitter::prepareCalculation()
{
  QMutexLocker lock(&m_Mutex);

  double bucketsperorbit = m_Acquisition -> get_BucketsPerOrbit();
  double orbitperiod = m_Acquisition -> get_SamplesPerOrbit();
  double bunchperiod = orbitperiod/bucketsperorbit;

  QVector<int> filledBuckets = m_Acquisition -> get_FilledBuckets();

  if (filledBuckets.size() == 0) {
//    printf("Constructing fake fill pattern");

    for (int i=0; i<24; i++) {
      filledBuckets.append(i*bucketsperorbit/24.0);
    }
  }

  int sz = m_Acquisition -> get_NSamples();

  m_ResultSize = 0;

  double t0 = 0;

//  printf("Ch%d: Samples per orbit %g, samples per bunch %g, buckets per orbit %g samples %d\n",
//         m_Channel, orbitperiod, bunchperiod, bucketsperorbit, sz);
//  printf("%d filled buckets\n", filledBuckets.size());

//   if (filledBuckets.size() == 24) {
//     for (int i=0; i<24; i++) {
//       printf("%d\t%d\t%g\n", i, filledBuckets.value(i), (i*bucketsperorbit/24.0));
//     }
//   }

  while (t0 < sz) {
    int b;
    foreach(b,filledBuckets) {
      double t = t0 + b*bunchperiod;

//      int i0 = (int) (t + get_FittingStart());
      int i1 = (int) (t + get_FittingEnd());

      if (i1 < sz) {
        m_ResultSize += 1;
      }
    }

    t0 += orbitperiod;
  }
}

void QavrgFitter::performCalculation()
{
  QMutexLocker lock(&m_Mutex);

//  printf("Fitter channel %d calculation started\n", m_Channel);
//  printf("Current thread %p, this->thread() %p\n", QThread::currentThread(), thread());

  QTime tm;
  tm.start();

  int sz = m_Acquisition -> get_NSamples();

  if (sz > m_RawData.size()) {
    printf("data record for channel %d shorter than expected (expected %d, actual %d, adjusting\n",
           m_Channel, sz, m_RawData.size());

    sz = m_RawData.size();
  }

  if (!get_ReferenceAvailable()) {
    printf("can't perform fitting if no reference is available\n");
    return;
  }

  if (sz > m_ReferenceData.size()) {
    printf("reference data record for channel %d shorter than expected (expected %d, actual %d, adjusting\n",
           m_Channel, sz, m_ReferenceData.size());

    sz = m_ReferenceData.size();
  }

  if (get_DarkAvailable() && sz > m_DarkData.size()) {
    printf("dark data record for channel %d shorter than expected (expected %d, actual %d, adjusting\n",
           m_Channel, sz, m_DarkData.size());

    sz = m_DarkData.size();
  }

  m_FitData.resize(sz);
  m_FitData.fill(0);

  double bucketsperorbit = m_Acquisition -> get_BucketsPerOrbit();
  double orbitperiod = m_Acquisition -> get_SamplesPerOrbit();
  double bunchperiod = orbitperiod/bucketsperorbit;

  QVector<int> filledBuckets = m_Acquisition -> get_FilledBuckets();

  if (filledBuckets.size() == 0) {
//    printf("Constructing fake fill pattern");

    for (int i=0; i<24; i++) {
      filledBuckets.append(i*bucketsperorbit/24.0);
    }
  }

  double* d = m_RawData.data();
  double* ref = m_ReferenceData.data();
  double* drk = m_DarkData.data();
  double* fit = m_FitData.data();

  m_Results.resize(3);
  m_Results[0].resize(0);
  m_Results[1].resize(0);
  m_Results[2].resize(0);

  double t0 = 0;

//  printf("Ch%d: Samples per orbit %g, samples per bunch %g\n", m_Channel, orbitperiod, bunchperiod);
//  printf("%d filled buckets\n", filledBuckets.size());

  int fitStart = get_FittingStart();
  int fitEnd   = get_FittingEnd();
  double fitOffset = get_FittingOffset();

  while (t0 < sz) {
    int b;
    foreach(b,filledBuckets) {
      double t = t0 + b*bunchperiod;

      int i0 = (int) (t + fitStart);
      int i1 = (int) (t + fitEnd);

      if (i1 < sz) {
        QavrgMatrix m(3,3);
        QVector<double> v(3);

        for (int i=0; i<3; i++) {
          for (int j=0; j<3; j++) {
            m(i,j) = 0;
          }
          v[i]=0;
        }

        if (get_DarkAvailable()) {
          for (int i=i0; i<i1; i++) {
            double x=i-i0;
            double r=ref[i] - drk[i] - fitOffset;
            double y=(d[i] - drk[i] - fitOffset);

            m(0,0) += 1;
            m(0,1) += x;
            m(0,2) += r;
            v[0]   += y;

            m(1,0) += x;
            m(1,1) += x*x;
            m(1,2) += x*r;
            v[1]   += y*x;

            m(2,0) += r;
            m(2,1) += r*x;
            m(2,2) += r*r;
            v[2]   += y*r;
          }
        } else {
          for (int i=i0; i<i1; i++) {
            double x=i-i0;
            double r=ref[i] - fitOffset;
            double y=(d[i] - fitOffset);

            m(0,0) += 1;
            m(0,1) += x;
            m(0,2) += r;
            v[0]   += y;
            
            m(1,0) += x;
            m(1,1) += x*x;
            m(1,2) += x*r;
            v[1]   += y*x;
            
            m(2,0) += r;
            m(2,1) += r*x;
            m(2,2) += r*r;
            v[2]   += y*r;
          }
        }

        QavrgMatrix::gaussj(m,v);

//        printf("Ch%d:Bnch%d:v0=%g:v1=%g:v2=%g\n", m_Channel, b, v[0], v[1], v[2]);

        if (get_DarkAvailable()) {
          for (int i=i0; i<i1; i++) {
            double x=i-i0;
            double r=ref[i] - drk[i] - fitOffset;

            fit[i] = v[0]+x*v[1]+r*v[2] + fitOffset;
          }
        } else {
          for (int i=i0; i<i1; i++) {
            double x=i-i0;
            double r=ref[i] - fitOffset;

            fit[i] = v[0]+x*v[1]+r*v[2] + fitOffset;
          }
        }

        m_Results[0].append(v[0]);
        m_Results[1].append(v[1]);
        m_Results[2].append(v[2]);

//        printf("Ch%d: Fit %g + %g*x +%g*ref\n", m_Channel, v[0], v[1], v[2]);
      }
    }

    t0 += orbitperiod;
  }
}

QVector<double> QavrgFitter::get_RawData(int start, int nbins)
{
  QMutexLocker lock(&m_Mutex);

  return m_RawData.mid(start,nbins);
}

QVector<double> QavrgFitter::get_RawData()
{
  QMutexLocker lock(&m_Mutex);

  return m_RawData;
}

double QavrgFitter::get_RawData(int bin)
{
  QMutexLocker lock(&m_Mutex);

  return m_RawData.value(bin);
}

QVector<double> QavrgFitter::get_ReferenceData(int start, int nbins)
{
  QMutexLocker lock(&m_Mutex);

  return m_ReferenceData.mid(start,nbins);
}

QVector<double> QavrgFitter::get_ReferenceData()
{
  QMutexLocker lock(&m_Mutex);

  return m_ReferenceData;
}

double QavrgFitter::get_ReferenceData(int bin)
{
  QMutexLocker lock(&m_Mutex);

  return m_ReferenceData.value(bin);
}

QVector<double> QavrgFitter::get_DarkData(int start, int nbins)
{
  QMutexLocker lock(&m_Mutex);

  return m_DarkData.mid(start,nbins);
}

QVector<double> QavrgFitter::get_DarkData()
{
  QMutexLocker lock(&m_Mutex);

  return m_DarkData;
}

double QavrgFitter::get_DarkData(int bin)
{
  QMutexLocker lock(&m_Mutex);

  return m_DarkData.value(bin);
}

QVector<double> QavrgFitter::get_FitData(int start, int nbins)
{
  QMutexLocker lock(&m_Mutex);

  return m_FitData.mid(start,nbins);
}

QVector<double> QavrgFitter::get_FitData()
{
  QMutexLocker lock(&m_Mutex);

  return m_FitData;
}

double QavrgFitter::get_FitData(int bin)
{
  QMutexLocker lock(&m_Mutex);

  return m_FitData.value(bin);
}

void QavrgFitter::readSettings(QSettings *settings)
{
  QMutexLocker lock(&m_Mutex);

  QcepProperty::readSettings(this, &staticMetaObject, tr("channels%1").arg(m_Channel), settings);
}

void QavrgFitter::writeSettings(QSettings *settings)
{
  QMutexLocker lock(&m_Mutex);

  QcepProperty::writeSettings(this, &staticMetaObject, tr("channels%1").arg(m_Channel), settings);
}

double QavrgFitter::inputFullScale(int index) const
{
  const double vals[] = {0.05, 0.1, 0.2, 0.5, 1.0, 2.0, 5.0};

  if (index < 0 || index > 6) {
    return 0;
  } else {
    return vals[index];
  }
}

double QavrgFitter::get_InputFullScale() const
{
//  QMutexLocker lock(&m_Mutex);

  return inputFullScale(get_InputGain());
}

void QavrgFitter::set_InputFullScale(double val)
{
  QMutexLocker lock(&m_Mutex);

  for (int i=0; i<7; i++) {
    if (inputFullScale(i) >= val) {
      set_InputGain(i);
      return;
    }
  }
}

void QavrgFitter::resize(int nsamples)
{
  QMutexLocker lock(&m_Mutex);

  if (nsamples != m_NSamples) {
    m_NSamples = nsamples;

    m_RawData.resize(m_NSamples);
    m_DarkData.resize(m_NSamples);
    m_ReferenceData.resize(m_NSamples);
    m_FitData.resize(m_NSamples);
  }
}

void QavrgFitter::setReferenceData()
{
  QMutexLocker lock(&m_Mutex);

  m_ReferenceData = m_RawData;

  set_ReferenceAvailable(true);
}

QString QavrgFitter::referenceDataPath()
{
  return QavrgSettings::settingsDirectory().filePath(tr("ref%1").arg(m_Channel));
}

void QavrgFitter::saveReferenceData()
{
  QMutexLocker lock(&m_Mutex);

  QFile file(referenceDataPath());

  if (!file.open(QIODevice::WriteOnly)) {
    printf("Couldn't save reference data for channel %d\n", m_Channel);
    return;
  }

  //   printf("File size = %ld\n", file.size());

  qint64 sz = m_ReferenceData.size()*sizeof(double);

  //   printf("Reference data size = %ld elems\n", m_ReferenceData.size());

  qint64 szwrt = file.write((char*) m_ReferenceData.data(), sz);

  //   printf("New file size = %ld\n", file.size());

  if (szwrt != sz) {
    printf("Incomplete write, wrote %ld of %ld.\n", (long) szwrt, (long) sz);
  }
}

void QavrgFitter::loadReferenceData()
{
  QMutexLocker lock(&m_Mutex);

  QFile file(referenceDataPath());

  if (!file.open(QIODevice::ReadOnly)) {
    printf("Couldn't load reference data for channel %d\n", m_Channel);
    return;
  }

  int sz = file.size();

  int nsamp = sz / sizeof(double);
  int nread = nsamp*sizeof(double);

  //   printf("File size %d, nsamp %d, nread %d\n", sz, nsamp, nread);

  if (nsamp > m_NSamples) {
    resize(nsamp);
  }

  m_ReferenceData.fill(0);

  //   printf("File size %d, nsamp %d, nread %d\n", sz, nsamp, nread);

  int szrd = file.read((char*) m_ReferenceData.data(), nread);

  //   printf("File size %d, nsamp %d, nread %d\n", sz, nsamp, nread);

  if (szrd != nread) {
    printf("Incomplete read, read %d of %d.\n", szrd, nread);
  }

  set_ReferenceAvailable(true);
}

void QavrgFitter::setDarkData()
{
  QMutexLocker lock(&m_Mutex);

  m_DarkData = m_RawData;

  set_DarkAvailable(true);
}

QString QavrgFitter::darkDataPath()
{
  return QavrgSettings::settingsDirectory().filePath(tr("dark%1").arg(m_Channel));
}

void QavrgFitter::saveDarkData()
{
  QMutexLocker lock(&m_Mutex);

  QFile file(darkDataPath());

  if (!file.open(QIODevice::WriteOnly)) {
    printf("Couldn't save dark data for channel %d\n", m_Channel);
    return;
  }

  //   printf("File size = %ld\n", file.size());

  qint64 sz = m_DarkData.size()*sizeof(double);

  //   printf("Dark data size = %ld elems\n", m_ReferenceData.size());

  qint64 szwrt = file.write((char*) m_DarkData.data(), sz);

  //   printf("New file size = %ld\n", file.size());

  if (szwrt != sz) {
    printf("Incomplete write, wrote %ld of %ld.\n", (long) szwrt, (long) sz);
  }
}

void QavrgFitter::loadDarkData()
{
  QMutexLocker lock(&m_Mutex);

  QFile file(darkDataPath());

  if (!file.open(QIODevice::ReadOnly)) {
    printf("Couldn't load dark data for channel %d\n", m_Channel);
    return;
  }

  int sz = file.size();

  int nsamp = sz / sizeof(double);
  int nread = nsamp*sizeof(double);

  //   printf("File size %d, nsamp %d, nread %d\n", sz, nsamp, nread);

  if (nsamp > m_NSamples) {
    resize(nsamp);
  }

  m_DarkData.fill(0);

  //   printf("File size %d, nsamp %d, nread %d\n", sz, nsamp, nread);

  int szrd = file.read((char*) m_DarkData.data(), nread);

  //   printf("File size %d, nsamp %d, nread %d\n", sz, nsamp, nread);

  if (szrd != nread) {
    printf("Incomplete read, read %ld of %ld.\n", (long) szrd, (long) nread);
  }

  set_DarkAvailable(true);
}

QMutex *QavrgFitter::mutex()
{
  return &m_Mutex;
}

double *QavrgFitter::rawDataPtr()
{
//  if (!m_Mutex.locked()) {
//    printf("Warning: QavrgFitter::rawDataPtr called on unlocked object\n");
//  }

  return m_RawData.data();
}

double *QavrgFitter::referenceDataPtr()
{
//  if (!m_Mutex.locked()) {
//    printf("Warning: QavrgFitter::referenceDataPtr called on unlocked object\n");
//  }

  return m_ReferenceData.data();
}

double *QavrgFitter::darkDataPtr()
{
//  if (!m_Mutex.locked()) {
//    printf("Warning: QavrgFitter::darkDataPtr called on unlocked object\n");
//  }

  return m_DarkData.data();
}