EWP-BplusToKstMuMu-AngAna/Code/FCNCFitter/sources/Core/fitter.hh

/**
 * @file fitter.hh
 * @author Christoph Langenbruch, Renata Kopecna
 * @date 2010-10-13
 *
 */

#ifndef FITTER_H
#define FITTER_H

#include <string>
#include <vector>
#include <mutex>

#include <TMinuit.h>

namespace fcnc {
    class event;
    class parameter;
    class parameters;
    class pdf;
    class plotter;
    class options;
class fitter;

///gives wether a i (absolutely) larger than b
template<class T> struct GreaterThan{
  bool operator()(const T a, const T b) const
  { return std::abs(a) > std::abs(b); }
};

///adds up all probability values in an efficient and numerically optimal manner
template<typename iterator> typename iterator::value_type add_results(iterator begin, iterator end){
  typedef typename iterator::value_type T;
  if (end == begin) return 0.;
  std::make_heap(begin, end, GreaterThan<T>());

  for (; begin + 1 < end; ) {
    T x = *begin;
    std::pop_heap(begin, end--, GreaterThan<T>());
    do {
      x += *begin;
      std::pop_heap(begin, end--, GreaterThan<T>());
    } while (std::abs(*begin) >= std::abs(x) && begin < end);
    *end++ = x;
    std::push_heap(begin, end, GreaterThan<T>());
  }
  return *begin;
}

///The fitter class implements an unbinned multithreaded maximum likelihood fit
class fitter {
public:
  ///global mutex for use of threads
  //static boost::mutex fitter_mutex;
  static std::mutex fitter_mutex;
  ///global fitter object for minuit
  static fitter* minuit_one_fitter;
  ///empirical constant which is subtracted to improve the numerical stability of the fit
  double empirical_constant;
  ///the minuit 1 minimizer object
  TMinuit* minuit_one;
  ///pointer to events to fit
  std::vector< std::vector<event>* > events;
  ///pointer to parameter sets to fit
  std::vector< parameters* > params;
  ///common parameters in different parameter sets
  std::vector< std::string > common_params;
  //plotters
  std::vector< plotter* > plots;
  ///pointer to the pdf(s)
  std::vector<pdf*> pdfs;
  ///pointer to the options object
  options* opts;
  ///if true, this fitter has already been initialized
  bool initialized;
  ///every thread adds the result of their likelihood calculation to this
  long double lh_tot;
  ///use square of weights for error correction in weighted fits
  bool square_weights;
public:
  ///minuit one FCN
  static void minuit_one_fcn(Int_t &npar, Double_t *grad, Double_t &lh, Double_t *params, Int_t iflag);
  ///constructor
  fitter(options* o);
  ///destructor
  ~fitter();
  /** fit the parameter set to the events using the pdf
     * @param probability pointer to pdf to use in the likelihood fit
     * @param parameters pointer to parameterset to use in the likelihood fit
     * @param events pointer to a vector of events for which to maximise the likelihood
     * @param index optional parameter to give the current run of a toystudy **/
  int fit(std::vector<pdf*> probs, std::vector<parameters*> parms, std::vector< std::vector<event> * > ev, std::string index="");
  /** fit the parameter set to the events using the pdf
     * @param vector of probability pointers to pdfs to use in the likelihood fit
     * @param vector of pointers to parametersets to use in the likelihood fit
     * @param vector of pointers to vector of events for which to maximise the likelihood
     * @param index optional parameter to give the current run of a toystudy **/
  int fit(pdf* probs, parameters* parms, std::vector<event>* ev, std::string index="");
  //set pdfs, parameters and events to extract likelihoods manually:
  void init(std::vector<pdf*> probs, std::vector<parameters*> parms, std::vector< std::vector<event> * > ev, std::string index="");
  void init(pdf* probs, parameters* parms, std::vector<event>* ev, std::string index="");
  ///add a parameter to minuit
  void define_parameter(int i, const parameter & p);
  ///reset all parameters to the start values
  void reset_param_values();
  ///calculates log (prod(p_i)) = sum lh(p_i), calls likelihood_thread for every thread
  double likelihood();
  ///calculates the log pdf for a single thread
  void likelihood_thread(int thread_id);
  //set plotter
  void set_plotters(plotter* plot);
  ///set plotters
  void set_plotters(std::vector<plotter*> plotters);
  ///define common parameters for all pdfs
  void set_common_parameters(std::vector<std::string> common_pars);
  ///is the parameter with name one of the common parameters
  bool is_common(std::string name);
  ///is this the first occurence of the parameter?
  bool is_first(std::string name, unsigned int idx_i, unsigned int idx_j);
  ///returns the nth minuit index for parameter with name
  int nth_minuit_index(std::string name, unsigned int n);
  ///returns the first minuit index for parameter with name
  int first_minuit_index(std::string name);

  ///returns wether this is an analysis with blinded parameters
  bool is_blind();
};


}//end namespace

//Print migrad status with explanation
void printMigradStatus(int status);
#endif