JGlobalfit.hh

Go to the documentation of this file.

#ifndef __JACOUSTICS__JGLOBALFIT__
#define __JACOUSTICS__JGLOBALFIT__
 
#include "JFit/JEstimator.hh"
 
#include "JAcoustics/JKatoomba_t.hh"
#include "JAcoustics/JSoundVelocity.hh"
#include "JAcoustics/JGeometry.hh"
#include "JAcoustics/JModel.hh"
#include "JAcoustics/JHit.hh"
#include "JAcoustics/JFitParameters.hh"
#include "JAcoustics/JTimeRange.hh"
 
 
/**
 * \file
 *
 * Fit function of acoustic model.
 * \author mdejong
 */
namespace JACOUSTICS {}
namespace JPP { using namespace JACOUSTICS; }
 
namespace JACOUSTICS {
 
  /**
   * Global fit of prameterised detector geometry to acoustics data.
   */
  struct JGlobalfit {
    /**
     * Constructor.
     *
     * \param  geometry           detector geometry
     * \param  V                  sound velocity
     * \param  parameters         parameters
     */
    JGlobalfit(const JGeometry& geometry, const JSoundVelocity& V, const JFitParameters& parameters) :
      parameters(parameters),
      estimator(geometry, V, parameters.option),
      evaluator(geometry, V, parameters.option),
      gandalf  (geometry, V, parameters.option)
    {
      using namespace JPP;
 
      evaluator.estimator.reset(getMEstimator(parameters.mestimator));
      gandalf  .estimator.reset(getMEstimator(parameters.mestimator));
    }
 
 
    /**
     * Auxiliary data structure for return value of fit.
     */
    template<class T>
    struct result_type {
      JModel& value;
      double  chi2;
      double  ndf;
      T       begin;
      T       end;
 
      /**
       * Get number of hits used in fit (after outlier removal).
       *
       * \return                size
       */
      int size() const
      {
        return std::distance(begin, end);
      }
      
      /**
       * Get time range of data.
       *
       * \return                time range
       */
      JTimeRange getTimeRange() const
      {
        JTimeRange range(JTimeRange::DEFAULT_RANGE());
 
        for (T i = begin; i != end; ++i) {
          range.include(i->getValue());
        }
 
        return range;
      }
    };
 
 
    /**
     * Fit.
     *
     * \param  __begin        begin of data
     * \param  __end          end   of data
     * \return                fit result
     */
    template<class T>
    result_type<T> operator()(T __begin, T __end)
    {
      using namespace std;
 
      JModel  model;
      JModel& result = gandalf.value;                                 // start value
 
      result = estimator(__begin, __end);                             // pre-fit
 
      if (parameters.stdev > 0.0) {                                   // remove outliers
 
        double chi2[] = { numeric_limits<double>::max(),
                          numeric_limits<double>::max() };
 
        for ( ; ; ) {
 
          T out = __end;
 
          double xmax = 0.0;
 
          for (T hit = __begin; hit != __end; ++hit) {
 
            const double x = fabs(hit->getValue() - estimator.getToA(result, *hit)) / hit->getSigma();
 
            if (x > xmax) {
              xmax = x;
              out  = hit;
            }
          }
 
          if (xmax > parameters.stdev) {
 
            if (chi2[0] == numeric_limits<double>::max()) {
              chi2[0] = evaluator(result, __begin, __end);
            }
 
            model = result;
 
            iter_swap(out, --__end);
 
            result  = estimator(__begin, __end);
            chi2[1] = evaluator(result, __begin, __end);
 
            if (chi2[0] - chi2[1]  >  parameters.stdev * parameters.stdev) {
 
              chi2[0] = chi2[1];
 
              continue;
 
            } else {
 
              iter_swap(out, __end++);
 
              result = model;
            }
          }
 
          break;
        }
      }
 
      const double chi2  =  gandalf  (__begin, __end)  /  gandalf.estimator->getRho(1.0);
      const double ndf   =  getWeight(__begin, __end)  -  gandalf.value.getN();
 
      return { gandalf.value, chi2, ndf, __begin, __end };            // return values
    }
 
 
    JFitParameters                parameters;
    JKatoomba<JEstimator>         estimator;
    JKatoomba<JAbstractMinimiser> evaluator;
    JKatoomba<JGandalf>           gandalf;
  };
}
 
#endif