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

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

#ifndef BU2KSTARMUMU_PARAMETERS_H
#define BU2KSTARMUMU_PARAMETERS_H

#include <parameters.hh>
#include <options.hh>
#include <assert.h>

//Simple struct to see whether fix or constrain a parameter
struct fixConstr{
    bool fix = false;
    bool constrain = false;
    //Constructor
    fixConstr(bool b_fix, bool b_constr);
};

namespace fcnc {

///class that implements the different parameter sets for Bd -> Kstar mu mu
class bu2kstarmumu_parameters: public parameters {
private:
    ///adds all parameters to the params vector
    void add_parameters();
    options* opts;

public: 
    ///constructor uses toy challenge version 1 data as default
    bu2kstarmumu_parameters(options* o);

    ///use default B+ values
    /// //THIS IS SOMETHING EVERY ANALYSIS HAS TO CHANGE ON THEIR OWN!
    void use_default_bkg(); //Just reset bkg, useful in the code down the line
    void use_default_observables(); //Reset S-wave observables, useful in the code down the line
    void use_default();

    //simple and quick loading and saving of parameters to/from txt file
    void load_param_values(std::string filename);
    void load_param(std::string filename, std::string parname);

    int  get_param_from_rootfile(std::string fileName, std::vector<std::string> names, int PDF, int bin, fixConstr FC);

    //Fix parameter to a value from a root file
    int fix_param_from_rootfile(std::string fileName, std::vector<std::string> names, int PDF, int bin);
    //Constrain parameter to a value from a root file
    int constrain_param_from_rootfile(std::string fileName, std::vector<std::string>names, int PDF, int bin);
    //Initialize parameter from a root file and let if float as it wishes
    int load_param_from_rootfile(std::string fileName,std::vector<std::string> names, int PDF, int bin);

    //TODO: remove as it is obsolete when done with pulls.cc
    void load_only_bckgnd_param_values(std::string filename);
    void load_only_Bmass_param_values(std::string filename);

    //Initialize parameters in the fit
    void init_mass_parameters(int PDF, int nBins, int bin, double defaultStepSize);    
    void init_angular_background_parameters(bool fitReference, double stepsize);
    void init_kpi_background_parameters(bool fitReference, double stepsize);
    void init_mkpi_pWave_parameters(bool fitReference, double stepsize);
    void init_mkpi_sWave_parameters(bool fitReference, double stepsize);
    void init_mass_background_parameters(int nBins, int bin, bool useLambda);
    void init_angular_parameters(int nBins, int bin, double stepsize, double scale, bool blind);
    void init_Bmass(std::string fileName, int PDF, double stepSize, fixConstr FC);
    void init_ang_parameters_fromRefDavid(double stepsize,  double scale, bool blind);
    void init_sWave_parameters(double stepsize);

    void save_param_values(std::string filename);

    //Effective q2
    parameter eff_q2   = parameter("eff_q2",  "{q^{2}}_eff");

    //Number of signal/background events
    parameter n_sig = parameter("n_sig", "N_\\mathrm{sig}");
    parameter n_bkg = parameter("n_bkg", "N_\\mathrm{bkg}");

    //Fraction of signal events, effectivelly (nsig/(nsig+nbkg)))
    parameter f_sig = parameter("f_sig",   "f_\\mathrm{sig}");

    //The B massneeds the full range of the B-mass fit,
    //as it  decides the range of the data used in the fit
    parameter m_b       = parameter("m_b",     "m_{B^{+}}");

    //All mass parameters
    parameter m_res_1   = parameter("m_res_1", "f_{CB}");
    parameter m_scale   = parameter("m_scale", "s_\\mathrm{m}");
    parameter m_sigma_1 = parameter("m_sigma_1", "#sigma_{1}^{CB}");
    parameter m_sigma_2 = parameter("m_sigma_2", "#sigma_{2}^{CB}");
    parameter alpha_1   = parameter("alpha_1", "#alpha_{1}^{CB}");
    parameter alpha_2   = parameter("alpha_2", "#alpha_{2}^{CB}");
    parameter n_1 = parameter("n_1", "n_{1}^{CB}");
    parameter n_2 = parameter("n_2", "n_{2}^{CB}");

    //These are for exponential background, either exp(-lambda*x) or exp(x/tau) 
    parameter fm_tau   = parameter("fm_tau",       "f_{m,#tau}");
    parameter m_lambda = parameter("m_lambda",     "\\lambda_\\mathrm{m}");
    parameter m_lambda_2 = parameter("m_lambda_2", "\\lambda_\\mathrm{m,2}");
    parameter m_tau    = parameter("m_tau",        "#tau_{m,1}");
    parameter m_tau_2  = parameter("m_tau_2",      "#tau_{m,2}");

    //Kstar related parameters /

    parameter mkstar      = parameter("mkst","m(K^{*+})");
    parameter mkstarplus = parameter("mkstz","m(K^{*+}_{x})"); //K1+ resonances
    parameter mkstp      = parameter("mkstp", "m(K^+\\pi^0)"); //TODO

    parameter gammakstar = parameter("gammakstar", "\\Gamma(K^{*+})"); //Width of K*
    parameter gammakstarplus = parameter("gamkstp", "\\Gamma(K^+\\pi^0)"); //Width of K1
    //parameter gamkstp    = parameter("gamkstp", "\\Gamma(K^+\\pi^0)");


    //Afb and FL, S- and P- parameters are in S_pars(int) and  P_pars(int)
    parameter fl  = parameter("Fl",  "F_{L}");
    parameter afb = parameter("Afb", "A_{FB}");

    //double exponential bkg
    parameter Fl  = parameter("Fl",  "F_{L}");
    parameter S1s = parameter("S1s", "S_{1s}");
    parameter S3  = parameter("S3", "S_{3}");
    parameter S4  = parameter("S4", "S_{4}");
    parameter S5  = parameter("S5", "S_{5}");
    parameter Afb = parameter("Afb", "A_{FB}");
    parameter S6s = parameter("S6s", "S_{6s}");
    parameter S7  = parameter("S7", "S_{7}");
    parameter S8  = parameter("S8", "S_{8}");
    parameter S9  = parameter("S9", "S_{9}");

    //less FF dependent variables P_i
    parameter P1 = parameter("P1", "P_{1}");//=2*S3/(1-Fl)
    parameter P2 = parameter("P2", "P_{2}");//=0.5*S6/(1-Fl)
    parameter P3 = parameter("P3", "P_{3}");//=-S9/(1-Fl)
    parameter P4 = parameter("P4", "P_{4}^{\\prime}");//=S4/sqrt(Fl Ft)
    parameter P5 = parameter("P5", "P_{5}^{\\prime}");//=S5/sqrt(Fl Ft)
    parameter P6 = parameter("P6", "P_{6}^{\\prime}");//=S7/sqrt(Fl Ft)
    parameter P8 = parameter("P8", "P_{8}^{\\prime}");//=S8/sqrt(Fl Ft)

    //S-wave parameters
    parameter FS = parameter("FS" ,"F_{S}");
    parameter SS1 = parameter("SS1","S_{S1}");
    parameter SS2 = parameter("SS2","S_{S2}");
    parameter SS3 = parameter("SS3","S_{S3}");
    parameter SS4 = parameter("SS4","S_{S4}");
    parameter SS5 = parameter("SS5","S_{S5}");

    //More Kstar related parameters
    parameter asphase = parameter("asphase", "arg(A_{S})");
    //Control the K* shape; one is ~p(K*) and one is ~1/p(K*),
    parameter a       = parameter("a", "a");
    parameter r       = parameter("r", "r");

    //Related to kappa(800), or f(800), it changed name apparently
    //Essentially a wide resonance
    parameter mf800     = parameter("mf800","m(f_{800})");
    parameter gammaf800 = parameter("gamf800","\\Gamma(f_{800})");
    parameter f800mag   = parameter("f800mag","|f_{800}|");
    parameter f800phase = parameter("f800phase","arg(f_{800})");

    //R (also refered to as d in the literature)
    //is a parameter in the mkpi parameterisation.
    //It refers to the effective hadronic size
    //See L861f (page 66) of https://cds.cern.ch/record/2318554/files/LHCb-ANA-2018-022.pdf
    parameter R = parameter("R", "R");
    //As a standard this is the radius of a B messon

    //Background parameters
    //Polynomial of maximal order of four
    //Each parameter correspons to x^n
    //Eg cbkgctl1 and cbkgctl3 are non-zero, rest is zero means polynomial
    // x^1 + x^3

    //ctl
    parameter cbkgctl0 = parameter("cbkgctl0", "c_{bkg}^{0}(cos#Theta_{L})");
    parameter cbkgctl1 = parameter("cbkgctl1", "c_{bkg}^{1}(cos#Theta_{L})");
    parameter cbkgctl2 = parameter("cbkgctl2", "c_{bkg}^{2}(cos#Theta_{L})");
    parameter cbkgctl3 = parameter("cbkgctl3", "c_{bkg}^{3}(cos#Theta_{L})");
    parameter cbkgctl4 = parameter("cbkgctl4", "c_{bkg}^{4}(cos#Theta_{L})");
    //ctk
    parameter cbkgctk0 = parameter("cbkgctk0", "c_{bkg}^{0}(cos#Theta_{K})");
    parameter cbkgctk1 = parameter("cbkgctk1", "c_{bkg}^{1}(cos#Theta_{K})");
    parameter cbkgctk2 = parameter("cbkgctk2", "c_{bkg}^{2}(cos#Theta_{K})");
    parameter cbkgctk3 = parameter("cbkgctk3", "c_{bkg}^{3}(cos#Theta_{K})");
    parameter cbkgctk4 = parameter("cbkgctk4", "c_{bkg}^{4}(cos#Theta_{K})");
    parameter cbkgctk5 = parameter("cbkgctk5", "c_{bkg}^{5}(cos#Theta_{K})");
    parameter cbkgctk6 = parameter("cbkgctk6", "c_{bkg}^{6}(cos#Theta_{K})");
    //phi
    parameter cbkgphi0 = parameter("cbkgphi0", "c_{bkg}^{0}(#phi)");
    parameter cbkgphi1 = parameter("cbkgphi1", "c_{bkg}^{0}(#phi)");
    parameter cbkgphi2 = parameter("cbkgphi2", "c_{bkg}^{0}(#phi)");
    parameter cbkgphi3 = parameter("cbkgphi3", "c_{bkg}^{0}(#phi)");
    parameter cbkgphi4 = parameter("cbkgphi4", "c_{bkg}^{0}(#phi)");

    //mkpi
    parameter cbkgmkpi0 = parameter("cbkgmkpi0", "c_{bkg}^{0}(m_{K#pi})");
    parameter cbkgmkpi1 = parameter("cbkgmkpi1", "c_{bkg}^{1}(m_{K#pi})");
    parameter cbkgmkpi2 = parameter("cbkgmkpi2", "c_{bkg}^{2}(m_{K#pi})");
    parameter cbkgmkpi3 = parameter("cbkgmkpi3", "c_{bkg}^{3}(m_{K#pi})");
    parameter cbkgmkpi4 = parameter("cbkgmkpi4", "c_{bkg}^{4}(m_{K#pi})");

    //Fancy background to mkpi in the p^2 dimension
    parameter cswavemkpi0 = parameter("cswavemkpi0", "c_{S}^{0}(m_{K#pi})");
    parameter cswavemkpi1 = parameter("cswavemkpi1", "c_{S}^{1}(m_{K#pi})");
    parameter cswavemkpi2 = parameter("cswavemkpi2", "c_{S}^{2}(m_{K#pi})");
    parameter cswavemkpi3 = parameter("cswavemkpi3", "c_{S}^{3}(m_{K#pi})");
    parameter cswavemkpi4 = parameter("cswavemkpi4", "c_{S}^{4}(m_{K#pi})");

    parameter cbkgp20 = parameter("cbkgp20", "c_{bkg}^{0}(p^{2})");
    parameter cbkgp21 = parameter("cbkgp21", "c_{bkg}^{1}(p^{2})");
    parameter cbkgp22 = parameter("cbkgp22", "c_{bkg}^{2}(p^{2})");
    parameter cbkgp23 = parameter("cbkgp23", "c_{bkg}^{3}(p^{2})");
    parameter cbkgp24 = parameter("cbkgp24", "c_{bkg}^{4}(p^{2})");

    //Nobody used this ever, probably not even in the PDF
    parameter nthreshold = parameter("nthreshold", "n_{thr.}");

    double J1s() const {
        return 3.0/4.0*(1.0-J1c()); // 3/4*(1-F_L)
    };
    double J1c() const {
        if (opts->fit_pprimes) assert(opts->fit_fl);
        if (opts->fit_fl)      return Fl();
        else return 1.0-4.0/3.0*S1s();
    };
    double J2s() const {
        return J1s()/3.0; //1/4(1-F_L)
    };
    double J2c() const {
        return -J1c(); //-F_L
    };
    double J3() const {//P1;//=S3/(1-Fl) correction P1=2S3/(1-FL)
        if (opts->fit_pprimes) return 0.5*P1()*(1.0-J1c());
        else                   return S3();
    };
    double J4() const {//P4;//=S4/sqrt(Fl Ft)
        if (opts->fit_pprimes) return P4()*sqrt(J1c()*(1.0-J1c()));
        else                   return S4();
    };
    double J5() const {//P5;//=S5/sqrt(Fl Ft)
        if (opts->fit_pprimes)  return P5()*sqrt(J1c()*(1.0-J1c()));
        else                    return S5();
    };
    double J6s() const {//P2;//=S6/(1-Fl) correction P2=0.5*S6s/(1-FL)
        if (opts->fit_pprimes)  return 2.0*P2()*(1.0-J1c());
        else if (opts->fit_afb) return 4.0/3.0*Afb();
        else                    return S6s();
    };
    double J6c() const {
        return 0.0;
    };
    double J7() const {//P6;//=S7/sqrt(Fl Ft)
        if (opts->fit_pprimes)  return P6()*sqrt(J1c()*(1.0-J1c()));
        else                    return S7();
    };
    double J8() const {//P8;//=S8/sqrt(Fl Ft)
        if (opts->fit_pprimes)  return P8()*sqrt(J1c()*(1.0-J1c()));
        else                    return S8();
    };
    double J9() const {//P3;//=S9/(1-Fl) correction P3=-S9/(1-FL)
        if (opts->fit_pprimes)  return -P3()*(1.0-J1c());
        else                    return S9();
    };
};


}

//Returns value or error of a parameter from a rootFile
double get_param_value_from_rootfile(std::string fileName, std::string names, int PDF, int bin);
double get_param_error_from_rootfile(std::string fileName, std::string names, int PDF, int bin);

//Returns a scaled number of events in given bin from total number of events in data
double eventsInBin_fraction(int bin, int run, int nBins, bool fromRef);

//Returns a scaled number of events in given bin from total number of events
//Eg I want to have 1000 events in 4 bins and 2 pdfs, it generates 155 signal + 100 bkg events in bin1, 125+120 in bin2, 180 + 50 in bin3 and 200+70 in bin4 as random example
void EventNumbers(unsigned int bin, unsigned int pdf,
                  double & fraction, unsigned int & eventnumbers,unsigned int TotalEvents,
                  unsigned int bins, unsigned int pdfs);

#endif