/** * @file bu2kstarmumu_pdf.hh * @author Christoph Langenbruch, Renata Kopecna * @date 2009-03-18 * */ #ifndef BU2KSTARMUMU_PDF_H #define BU2KSTARMUMU_PDF_H #include #include #include #include //Create vectors with the values of the polynomials std::vector init_ch_ctl(const fcnc::bu2kstarmumu_parameters* params); std::vector init_ch_ctk(const fcnc::bu2kstarmumu_parameters* params); std::vector 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* events); //void update_cached_normalization(const parameters* params, std::vector* events); void update_cached_efficiencies(const bu2kstarmumu_parameters* params, std::vector* events); void update_cached_efficiencies(const parameters* params, std::vector* events); void update_cached_xis(const bu2kstarmumu_parameters* params, std::vector* events); void update_cached_xis(const parameters* params, std::vector* 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 coeffs_eff_4d; ///chebyshev coefficients for factorizing bkg modeling //std::vector coeffs_bkg_costhetal; //std::vector coeffs_bkg_costhetak; //std::vector coeffs_bkg_phi; //std::vector bkg_events; std::vector 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 events, std::vector& moments, std::vector& covariance) const; void save_moments_to_obs(bu2kstarmumu_parameters* params, std::vector obs, std::vector 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& 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::vectorread_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& events, const bu2kstarmumu_parameters* theory_params); ///parametrize the angular background distribution from the bd mass sidebands //void parametrize_bkg(const std::vector& 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(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* 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