/**
 * @file bu2kstarmumu_pdf.hh
 * @author Christoph Langenbruch, Renata Kopecna
 * @date 2009-03-18
 *
 */

#ifndef BU2KSTARMUMU_PDF_H
#define BU2KSTARMUMU_PDF_H

#include <bu2kstarmumu_parameters.hh>
#include <pdf.hh>
#include <parameters.hh>
#include <folder.hh>

//Create vectors with the values of the polynomials
std::vector<double> init_ch_ctl(const fcnc::bu2kstarmumu_parameters* params);
std::vector<double> init_ch_ctk(const fcnc::bu2kstarmumu_parameters* params);
std::vector<double> init_ch_phi(const fcnc::bu2kstarmumu_parameters* params);

namespace fcnc {

  struct q2bin {
    double min;
    double max;
    q2bin(double a, double b): min(a), max(b) {};
  };

  class options;
  class event;
  class values;

  static const double mumass = 105.6583668; //TODO: OMFG

  class bu2kstarmumu_pdf;

  ///class that implements the pdf ot the signal decay Bs -> JPsi Phi
  class bu2kstarmumu_pdf: public pdf {
      ///These encapsulation functions used by ROOT/Minuit need access to class internals.
      friend class bu2kstarmumu_plotter;
      friend class bu2kstarmumu_generator;
  private:
      folder* fldr;
  public:
      static bu2kstarmumu_pdf* current_pdf;
      //double shift_lh;
      bu2kstarmumu_pdf(options* o, const bu2kstarmumu_parameters* params);
      virtual ~bu2kstarmumu_pdf();
      void init(bu2kstarmumu_parameters* params)
      {
      };
      void init(parameters* params);
      void update_cached_normalization(const bu2kstarmumu_parameters* params);
      void update_cached_normalization(const parameters* params);
      void update_cached_integrated_fis(const bu2kstarmumu_parameters* params);
      void update_cached_integrated_fis(const parameters* params);
      //void update_cached_normalization(const bu2kstarmumu_parameters* params, std::vector<event>* events);
      //void update_cached_normalization(const parameters* params, std::vector<event>* events);
      void update_cached_efficiencies(const bu2kstarmumu_parameters* params, std::vector<event>* events);
      void update_cached_efficiencies(const parameters* params, std::vector<event>* events);
      void update_cached_xis(const bu2kstarmumu_parameters* params, std::vector<event>* events);
      void update_cached_xis(const parameters* params, std::vector<event>* events);
  public:
      ///normalization
      double fnorm_sig;
      //double fnorm_sig_swave;
      double fnorm_bkg;
      double fbuffer_f1, fbuffer_f2, fbuffer_f3, fbuffer_f4, fbuffer_f5, fbuffer_f6, fbuffer_f7, fbuffer_f8, fbuffer_f9, fbuffer_f10, fbuffer_f11, fbuffer_f12;
      double fbuffer_fs1, fbuffer_fs2, fbuffer_fs3, fbuffer_fs4, fbuffer_fs5, fbuffer_fs6;
      virtual void get_swave_integrated_fj(double& f1, double& f2, double& f3, double& f4, double& f5, double& f6) const {
          f1 = fbuffer_fs1;
          f2 = fbuffer_fs2;
          f3 = fbuffer_fs3;
          f4 = fbuffer_fs4;
          f5 = fbuffer_fs5;
          f6 = fbuffer_fs6;
          return;
      };
      virtual void get_integrated_fj(double& f1, double& f2, double& f3, double& f4, double& f5, double& f6,
                                     double& f7, double& f8, double& f9, double& f10, double& f11, double& f12) const {
          f1 = fbuffer_f1;
          f2 = fbuffer_f2;
          f3 = fbuffer_f3;
          f4 = fbuffer_f4;
          f5 = fbuffer_f5;
          f6 = fbuffer_f6;
          f7 = fbuffer_f7;
          f8 = fbuffer_f8;
          f9 = fbuffer_f9;
          f10 = fbuffer_f10;
          f11 = fbuffer_f11;
          f12 = fbuffer_f12;
          return;
      };
      /*
      unsigned int neff_costhetal;
      unsigned int neff_costhetak;
      unsigned int neff_phi;

      unsigned int nbkg_costhetal;
      unsigned int nbkg_costhetak;
      unsigned int nbkg_phi;
      */
      ///coefficients for 4D angular acceptance (Sliding acceptance which is controlled by a parameter eff_q2range in the parameter set)
      std::vector<double> coeffs_eff_4d;

      ///chebyshev coefficients for factorizing bkg modeling
      //std::vector<double> coeffs_bkg_costhetal;
      //std::vector<double> coeffs_bkg_costhetak;
      //std::vector<double> coeffs_bkg_phi;

      //std::vector<event> bkg_events;
      std::vector<event> eff_events;
      const bu2kstarmumu_parameters* eff_params;

      //static void bkg_fcn(Int_t &npar, Double_t *grad, Double_t &lh, Double_t *params, Int_t iflag);
      static void eff_fcn(Int_t &npar, Double_t *grad, Double_t &lh, Double_t *params, Int_t iflag);
      static void eff_fcn_phsp(Int_t &npar, Double_t *grad, Double_t &lh, Double_t *params, Int_t iflag);
      static void eff_fcn_phsp_4d(Int_t &npar, Double_t *grad, Double_t &lh, Double_t *params, Int_t iflag);

      //////////////////
      ///Full angular///
      //////////////////

      //Angular correction polynomials are converted into Sum (c_ijkl * ctk^i * ctl^j * phi^k * q2^l)
      //S-wave
      void swave_fj(double ctl, double ctk, double phi,
                    double& f1, double& f2, double& f3, double& f4, double& f5, double& f6) const;
      void swave_integrated_fj_noacc(double& f1, double& f2, double& f3, double& f4, double& f5, double& f6) const;
      void swave_integrated_fj_noacc(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b,
                                     double& f1, double& f2, double& f3, double& f4, double& f5, double& f6) const;
      void swave_integrated_fj_chebyshev(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b,
                                         double& f1, double& f2, double& f3, double& f4, double& f5, double& f6, double q2) const;
      void swave_integrated_fj_chebyshev(double& f1, double& f2, double& f3, double& f4, double& f5, double& f6, double q2) const;

      //P-wave
      void fj(double ctl, double ctk, double phi,
              double& f1, double& f2, double& f3, double& f4, double& f5, double& f6,
              double& f7, double& f8, double& f9, double& f10, double& f11, double& f12) const;
      void integrated_fj_noacc(double& f1, double& f2, double& f3, double& f4, double& f5, double& f6,
                               double& f7, double& f8, double& f9, double& f10, double& f11, double& f12) const;
      void integrated_fj_noacc(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b,
                               double& f1, double& f2, double& f3, double& f4, double& f5, double& f6,
                               double& f7, double& f8, double& f9, double& f10, double& f11, double& f12) const;
      void integrated_fj_chebyshev(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b,
                                                     double& f1, double& f2, double& f3, double& f4, double& f5, double& f6,
                                                     double& f7, double& f8, double& f9, double& f10, double& f11, double& f12, double q2) const;
      void integrated_fj_chebyshev(double& f1, double& f2, double& f3, double& f4, double& f5, double& f6,
                                   double& f7, double& f8, double& f9, double& f10, double& f11, double& f12, double q2) const;

      /////////////////////
      ///Angular folding///
      /////////////////////

      //S-wave
      void folded_swave_fj(double ctl, double ctk, double phi, double& f1, double& f2, double& f3, double& f4, double& f5, double& f6) const;
      void folded_swave_integrated_fj_noacc(double& f1, double& f2, double& f3, double& f4, double& f5, double& f6) const;
      void folded_swave_integrated_fj_noacc(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b,
                                            double& f1, double& f2, double& f3, double& f4, double& f5, double& f6) const;
      void folded_swave_integrated_fj_chebyshev(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b,
                                                double& f1, double& f2, double& f3, double& f4, double& f5, double& f6, double q2) const;
      void folded_swave_integrated_fj_chebyshev(double& f1, double& f2, double& f3, double& f4, double& f5, double& f6, double q2) const;

      //P-wave
      void folded_fj(double ctl, double ctk, double phi, double& f1, double& f2, double& f3, double& f4, double& f5, double& f6, double& f7, double& f8, double& f9, double& f10, double& f11, double& f12) const;
      void folded_integrated_fj_noacc(double& f1, double& f2, double& f3, double& f4, double& f5, double& f6, double& f7, double& f8, double& f9, double& f10, double& f11, double& f12) const;
      void folded_integrated_fj_noacc(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b,
                                      double& f1, double& f2, double& f3, double& f4, double& f5, double& f6, double& f7, double& f8, double& f9, double& f10, double& f11, double& f12) const;
      void folded_integrated_fj_chebyshev(double& f1, double& f2, double& f3, double& f4, double& f5, double& f6, double& f7, double& f8, double& f9, double& f10, double& f11, double& f12, double q2) const;
      void folded_integrated_fj_chebyshev(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b,
                                          double& f1, double& f2, double& f3, double& f4, double& f5, double& f6, double& f7, double& f8, double& f9, double& f10, double& f11, double& f12, double q2) const;

      //background
      double integral_bkg(const bu2kstarmumu_parameters* params, double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b) const;
      double integral_bkg_chebyshev(const bu2kstarmumu_parameters* params, double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, double q2) const;

      //moments
      void get_moments(std::vector<fcnc::event> events, std::vector<double>& moments, std::vector<double>& covariance) const;
      void save_moments_to_obs(bu2kstarmumu_parameters* params, std::vector<double> obs, std::vector<double> obscov) const;

      //mkpi dependencies
      double mkpi_pwave_2(const bu2kstarmumu_parameters* params, const event& meas) const;
      double mkpi_swave_2(const bu2kstarmumu_parameters* params, const event& meas) const;
      double mkpi_re_swavepwave(const bu2kstarmumu_parameters* params, const event& meas) const;
      double mkpi_im_swavepwave(const bu2kstarmumu_parameters* params, const event& meas) const;
      double mkpi_bkg(const bu2kstarmumu_parameters* params, const event& meas) const;
      double mkpi_swave_2_norm;
      double mkpi_pwave_2_norm;
      double mkpi_re_swavepwave_norm;
      double mkpi_im_swavepwave_norm;
      double mkpi_bkg_norm;

      ///determine parametrization for efficiencies in angles and q2, needs reconstructed phsp MC events, needs to change weights so non-const
      void parametrize_eff_phsp_4d(const std::vector<event>& events, values* globalvalues, int tagNumber, bool assumePhiEven, bool checkSignifcance, bool runMinuit, bool checkFactorization, bool do3Dmoments);

      //Get the efficiency from the angular corrections
      double get_ang_eff(fcnc::options *opts, double q2, double ctl, double ctk, double phi) const;
      double get_ang_eff(fcnc::options *opts, fcnc::event evt, bool fa) const; //TBH, no clue why there has to be the const, but it compiles, so I'll take it

      //save and load the 4d coefficients to a file
      void save_coeffs_eff_phsp_4d();
      void load_coeffs_eff_phsp_4d();
      std::vector<double>read_coeffs_eff_phsp_4d();
      ///determine parametrization for efficiencies, needs the reconstructed (MC) events and the set of theory paramters
      void parametrize_eff(const std::vector<event>& events, const bu2kstarmumu_parameters* theory_params);
      ///parametrize the angular background distribution from the bd mass sidebands
      //void parametrize_bkg(const std::vector<event>& events);
      //can do this unbinned in principle ?
      //yes, do not forget that normalization is not free -> -1 parameter

      ///the options used for this pdf e.g. per event resolution or not etc.
      options* opts;
      ///The signal mass probability (Currently double gaussian). //TODO: check
      double m_sig_prob(const bu2kstarmumu_parameters* params, const event& meas) const;
      double integral_m_sig_prob(const bu2kstarmumu_parameters* params, double ma, double mb) const;
      ///The prompt background mass probability (Currently exponential).
      double m_bkg_prob(const bu2kstarmumu_parameters* params, const event& meas) const;
      double integral_m_bkg_prob(const bu2kstarmumu_parameters* params, double ma, double mb) const;
      ///The signal m(Kpi) mass probability
      double mkpi_sig_prob(const bu2kstarmumu_parameters* params, const event& meas) const;
      ///The prompt background in m(Kpi) spectrum probability //TODO: normalization?
      double mkpi_bkg_prob(const bu2kstarmumu_parameters* params, const event& meas) const;
      ///The signal angular probability
      double angular_sig_prob(const bu2kstarmumu_parameters* params, const event& meas) const;
      ///The prompt background angular probability (flat or Chebyshev)
      double angular_bkg_prob(const bu2kstarmumu_parameters* params, const event& meas) const;
      ///Simple background probabilily:
      double simple_angular_bkg_prob(const bu2kstarmumu_parameters* params, const event& meas) const;

      ///The total probability f_sig*p_sig + (1-f_sig)*p_bkg for event meas.
      virtual double prob(const parameters* params, const event& meas) const{
          return prob(static_cast<const bu2kstarmumu_parameters*>(params), meas);//still have the initialization as safety, this is faster
      };
      ///The total probability f_sig*p_sig + (1-f_sig)*p_bkg for event meas.
      double prob(const bu2kstarmumu_parameters* params, const event& meas) const;
      void calculate_sweights(const bu2kstarmumu_parameters* parms, std::vector<event>* ev) const;

      double angular_bkg_norm(const bu2kstarmumu_parameters* params, bool weighted);
      double angular_sig_norm(const bu2kstarmumu_parameters* params, bool weighted);
      double angular_swave_norm(const bu2kstarmumu_parameters* params, bool weighted);
      double angular_pwave_norm(const bu2kstarmumu_parameters* params, bool weighted);
      double angular_pswave_norm(const bu2kstarmumu_parameters* params, bool weighted);
  };


}

struct npolynom{
    unsigned int ctl;
    unsigned int ctk;
    unsigned int phi;
    unsigned int q2;

    unsigned int getSize(){
        return q2*ctl*ctk*phi;
    }

    unsigned int getSize2(){
        return getSize()*getSize();
    }

    npolynom(fcnc::options *opts){
        //Set orders of chebyshev polynomials
        ctl = opts->eff_order_costhetal;
        ctk = opts->eff_order_costhetak;
        phi       = opts->eff_order_phi;
        q2        = opts->eff_order_q2;
        //spdlog::trace("phi: {0:d}", phi);
        //spdlog::trace("opts->eff_order_phi: {0:d}", opts->eff_order_phi);
    }
    unsigned int get_bin_in4D(int h, int i, int j, int k);
};


#endif
