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