/* -*- C++ -*- *************************************************************
 *
 *   CbmRichRingFinderImp.h, Thu Jul 15 11:43:23 CEST 2004
 *   lsf4rich.h, Tue Jun 22 17:11:21 MSD 2004
 *   Copyright (C) 2004 Soloviev Alexei <solovjev@cv.jinr.ru>
 *
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#ifndef LSF4RICH_H
#define LSF4RICH_H

/***************************************************************************/
class LSF4RICH;
typedef LSF4RICH CbmRichRingFinderImp;
/***************************************************************************/

#include "TMinuit.h"

#include <complex>
#include <vector>
#include <iostream>
using namespace std;

/***************************************************************************
 *                                                                         *
 *  Least Square Function (LSF) for RICH                                   *
 *                                                                         *
 *  Usage:                                                                 *
 *                                                                         *
 *  LSF4RICH sigma;                     // create object                   *
 *  sigma.clear();                      // clean it if necessary           *
 *  while ( ... ) {                                                        *
 *    sigma.add_data_point ( x, y );    // fill if with data in a loop     *
 *  }                                                                      *
 *  if ( ... ) {                                                           *
 *    sigma.set_guidance ( x, y, r );   // set ring guidance manually...   *
 *  } else {                                                               *
 *    sigma.find_guidance();            // or find ring guidance over data *
 *  }                                                                      *
 *  sigma.minimize ( robust, verbose ); // perform fitting                 *
 *  x = sigma.ring_center_x();          // output ring center x coordinate *
 *  y = sigma.ring_center_y();          // output ring center y coordinate *
 *  R = sigma.ring_radius();            // output ring radius              *
 *                                                                         *
 ***************************************************************************/

class LSF4RICH : private TMinuit {

public: /** C-tor & d-tor */

  LSF4RICH();
  ~LSF4RICH();

public: /** Interface */

  /** Clear data storage and ring guidance */
  void clear();
  /** Add point to data storage */
  void add_data_point ( const double & x, const double & y );
  /** Ring center x coordinate access */
  const double & ring_center_x() const { return __ring_center_x; }
  /** Ring center y coordinate access */
  const double & ring_center_y() const { return __ring_center_y; }
  /** Ring radius access */
  const double & ring_radius() const { return __ring_radius; }
  /** Set sought ring guidance */
  void set_guidance ( const double & x, const double & y, const double & r );
  /** Set sought ring guidance over data */
  void find_guidance();
  /** Set threshold value for norm (sum of weights) */
  void set_threshold ( const double & t ) { __threshold = t; }
  /** LSF value (to be minimized) */
  const double value() const;
  /** MINUIT interface */
  bool minimize ( bool robust = false, const int verbose = -1 );

private: /** Functions to be minimized */

  /** Standard LSF value */
  const double standard() const;
  /** Weighted LSF value */
  const double weighted() const;

private: /** Storages */

  /** Data storage */
  typedef complex<double> point;
  vector <point> __data;
  /** Ring storages */
  mutable double __ring_center_x;
  mutable double __ring_center_y;
  mutable double __ring_radius;
  /** Pointer to function to be minimized */
  typedef const double (LSF4RICH::*fun4min) () const;
  fun4min __fun_ptr;

private: /** Robustness stuff */

  /** Weights for weighted LSF value */
  vector<double> __weight;
  /** Norm (sum of weights) for weighted LSF value */
  double __norm;
  /** Calculate Tukey's bi-square weights for weighted LSF value */
  void set_weights();
  /** Set weights = 1 for weighted LSF value 
      (should be called before any weights uasge!) */
  void init_weights();
  /** Constant for Tukey's bi-weights calculation */
  double C_TUKEY;

private: /** Guidance stuff */

  /** Has guidance? */
  bool __guided;
  /** Set weights = 0 for points out of guidance */
  void guide_weights();
  /** Threshold for norm (sum of weights) */
  double __threshold;

private: /** MINUIT stuff (defined in lsf4rich_minuit.cpp) */

  /** Static pointer to this for FCN */
  static const LSF4RICH * __lsf_ptr;
  /** User FCN function for MINUIT */
  static void fcn ( int &, double *, double &, double *, int );
  /** Set printout level for MINUIT */
  void set_mn_printout ( const int level = -1 );
  /** Set parameters for MINUIT */
  void set_mn_parms();
  /** Get parameters from MINUIT */
  void get_mn_parms() const;

};

/** LSF object output */
ostream & operator<< ( ostream & out, const LSF4RICH & lsf );

#endif // LSF4RICH_H