//Renata Kopecna

#include <event.hh>
#include <funcs.hh>
#include <bu2kstarmumu_generator.hh>
#include <helpers.hh>

#include <thread>
#include <mutex>

#include <spdlog.h>

//using namespace std;
using namespace fcnc;

//boost::mutex generator_mutex;
std::mutex generator_mutex;

//Function that returns a parameter with given range
double genValueFromRange(double min, double max, TRandom3* rnd){
    return min+rnd->Rndm()*(max-min);
}
double genValueFromPar(fcnc::parameter *par,  TRandom3* rnd){
    return genValueFromRange(par->get_min(),par->get_max(),rnd);
}


std::vector<event> bu2kstarmumu_generator::generate(unsigned int nevents, bu2kstarmumu_parameters* parameters, bu2kstarmumu_pdf* probability){
    prob = probability;
    params = parameters;

    //even if folding (non-full-angular) is chosen, use full_angular mode to create events!
    assert(opts->full_angular == prob->opts->full_angular);
    assert(opts->always_generate_full_angular == prob->opts->always_generate_full_angular);    
    //always_generate_full_angular is true everywhere anyway

    bool only_generate_full_angular = prob->opts->always_generate_full_angular && !prob->opts->full_angular;
    spdlog::debug("only_generate_full_angular: " + boolToString(only_generate_full_angular));
    if(only_generate_full_angular){
        prob->opts->full_angular = true;
        opts->full_angular = true;
        prob->opts->update_angle_ranges();
        opts->update_angle_ranges();
        if(opts->use_angular_acc){
            prob->load_coeffs_eff_phsp_4d();
            spdlog::debug("Loaded full angular coefficients.");
        }
    }

    prob->init(params);

    //account for the bias in f_sig by rescaling the f_sig for production, using efficiency corrected normalization and non-corrected normalization of background and signal PDF
    //TODO: ask David wth this comment means
    if((opts->weighted_fit || (opts->use_weighted_bkg && opts->use_angular_acc))){
        //Get normalizations
        double bkg_norm     = prob->angular_bkg_norm(params, false);
        double bkg_eff_norm = prob->angular_bkg_norm(params, true);
        double sig_norm     = prob->angular_sig_norm(params, false);
        double sig_eff_norm = prob->angular_sig_norm(params, true);
        double f_sig_scale = (bkg_norm / bkg_eff_norm) / (sig_norm / sig_eff_norm);

        spdlog::info("BKG_NORM={0:f}", bkg_norm);
        spdlog::info("BKG_NORM_EFF={0:f}", bkg_eff_norm);

        //make sure the scale is not too large, because this might lead to unwanted effects
        if(f_sig_scale > 2.0){
            spdlog::error("The rescaling of the f_sig in toy generation is larger than 2.0 (value is: {0:f}), which seems very odd. Check your parameters!",f_sig_scale);
            assert(f_sig_scale < 2.0);
        }
        if(f_sig_scale < 0.5){
            spdlog::error("The rescaling of the f_sig in toy generation is smaller 0.5 (value is: {0:f}), which seems very odd. Check your parameters!",f_sig_scale);
            assert(f_sig_scale > 0.5);
        }

        //set effective f_sig:
        if(params->f_sig() != 0.0 && params->f_sig() != 1.0){
            effective_f_sig = params->f_sig() * f_sig_scale / (1 + params->f_sig() * (f_sig_scale - 1));
        }
        else{
            effective_f_sig = params->f_sig();
        }

        //test if resulting 'effective' f_sig is out of the range [0.0, 1.0]
        if(effective_f_sig > 1.0 || effective_f_sig < 0.0){
            spdlog::error("Rescaled f_sig in generation is out of the allowed range (0-1). Use either pure background or pure signal toys or a more equal mixture of both!");
            assert(0);
        }
        spdlog::debug("[INFO]\t\tSignal fraction: f_sig={0:0.4f} and eff_f_sig={1:f}", params->f_sig(), effective_f_sig);
    }
    else{
        effective_f_sig = params->f_sig();
    }

    prob->update_cached_normalization(params);
    spdlog::info("Generating Toy MC");
    //Create empty space for the events
    std::vector< std::vector<event> > generated_events;
    for (unsigned int i = 0; i < opts->ncores; i++) {
        generated_events.push_back(std::vector<event>());
        generated_events.at(i).reserve(nevents/opts->ncores);
    }

    spdlog::info("Generation started");
    //Divide in threads if possible
    std::vector<std::thread*> threads;
    if (opts->ncores > 1 && nevents > 1000) {
        for (unsigned int i = 0; i < opts->ncores; i++) {
            //boost::thread* t = new boost::thread(boost::bind(&bu2kstarmumu_generator::generation_thread, this, i, nevents/opts->ncores, &generated_events.at(i)));
            std::thread* t = new std::thread(std::bind(&bu2kstarmumu_generator::generation_thread, this, i, nevents/opts->ncores, &generated_events.at(i)));
            threads.push_back(t);
        }
        for (unsigned int i = 0; i < opts->ncores; i++) threads[i]->join();
    }
    else generation_thread(0, nevents, &generated_events.at(0));

    //Copy the generated events into the final vector
    std::vector< event > evtVec;
    for (unsigned int i = 0; i < opts->ncores; i++){
        copy(generated_events.at(i).begin(), generated_events.at(i).end(), back_inserter(evtVec));
    }
    if (opts->ncores > 1 && nevents > 1000) {
        for (unsigned int i = 0; i < opts->ncores; i++) delete threads[i];
    }

    //generate last events on single thread if number of requested event is not multiple of nCPUcores
    //and add it to the final vector
    if(evtVec.size() < nevents){
        std::vector<fcnc::event> last_events;
        generation_thread(0, nevents - evtVec.size(), &last_events);
        evtVec.insert(evtVec.end(), last_events.begin(), last_events.end());
    }

    if(evtVec.size() != nevents){
        spdlog::critical("Number of requested events not reached: {0:d}/{1:d} ({2:f}%) ",
                         evtVec.size(), nevents, 100. * evtVec.size() / nevents);
        assert(0);
    }
    
    //TODO: ask David what is this supposed to be
    //opts->use_angular_res = angular_res;
    spdlog::info("Generation of {0:d} pseudo-events completed", nevents);

    //reset full_angular to false, if it is only needed for the generation of events
    //Aka when fititng and folding and you need a generator before that, restore the
    //folding back to being used
    if(only_generate_full_angular){
        prob->opts->full_angular = false;
        opts->full_angular = false;
        prob->opts->update_angle_ranges();
        opts->update_angle_ranges();
    }

    return evtVec;
};

void bu2kstarmumu_generator::generation_thread(unsigned int thread_id, unsigned long nevents, std::vector<event>* eventVec){
    //nevents is a number of actually generated events

    //TRandom3 has an internal state!
    //we therefore need one TRandom3 per thread
    TRandom3* rnd = new TRandom3();
    if (opts->always_static_seed){//always use the same static seed
        rnd->SetSeed(71598+thread_id*12); //TODO
    }
    else if (opts->static_seed){//this needs to be rerun dependent otherwise every rerun you get the same data!
        spdlog::info("Thread={0:d}\tJobID={1:d}\tSeed={2:d}", thread_id, opts->get_job_id(), opts->get_job_id()*opts->ncores+thread_id+12345);
        rnd->SetSeed(opts->get_job_id()*opts->ncores+thread_id+12345);
    }
    else rnd->SetSeed(0);//this is correct, see documentation //WHY IS THERE NO LINK AND ONE HAS TO GOOGLE?!

    spdlog::debug("Generator OK");

    double mkpi_min = opts->mkpi_min/1.0e+3;
    double mkpi_max = opts->mkpi_max/1.0e+3;
    double asphase = params->asphase();
    double a = params->a();
    double r = params->r();
    double gammakstar     = params->gammakstar();
    double mkstar         = params->mkstar();
    double gammakstarplus = params->gammakstarplus();
    double mkstarplus     = params->mkstarplus();
    double gammaf800      = params->gammaf800();
    double mf800          = params->mf800();
    double f800mag        = params->f800mag();
    double f800phase      = params->f800phase();
    double q2             = params->eff_q2();

    double R              = params->R();

    bool swave = opts->swave;

    //integrate over mkpi
    const double mk   = PDGMASS_KST_KAON/1.0e+3;
    const double mpi  = PDGMASS_KST_PION/1.0e+3;
    const double mb   = PDGMASS_B/1.0e+3;
    const double from = opts->mkpi_full_range_norm ? (mk+mpi) : mkpi_min;
    const double to   = opts->mkpi_full_range_norm ?  mb      : mkpi_max;
    double mkpi_swave_2_norm=1.0, mkpi_pwave_2_norm=1.0;
    double int_mkpi_re_swavepwave=1.0, int_mkpi_im_swavepwave=1.0;
    //TODO: check whether I can turn off the mkpi generation all together
    if(opts->generate_mkpi){
        spdlog::debug("Generating m(Kpi)!");
        if(opts->simple_mkpi){
            mkpi_swave_2_norm = fcnc::int_mkpi_simple_kstarzero_amp_squared(from, to, q2, f800mag, f800phase, gammaf800, mf800, gammakstarplus, mkstarplus);
            mkpi_pwave_2_norm = fcnc::int_mkpi_simple_bw_kstar_amp_squared(from, to, q2, gammakstar, mkstar);
            int_mkpi_re_swavepwave = fcnc::int_mkpi_simple_bw_kstar_amp_bar_kstarzero_amp(from, to, q2, gammakstar, mkstar, asphase, f800mag, f800phase, gammaf800, mf800, gammakstarplus, mkstarplus).real();
            int_mkpi_im_swavepwave = fcnc::int_mkpi_simple_bw_kstar_amp_bar_kstarzero_amp(from, to, q2, gammakstar, mkstar, asphase, f800mag, f800phase, gammaf800, mf800, gammakstarplus, mkstarplus).imag();
        }
        else if(opts->isobar){
            mkpi_swave_2_norm = fcnc::int_mkpi_kstarzero_isobar_amp_squared(from, to, q2, f800mag, f800phase, gammaf800, mf800, gammakstarplus, mkstarplus, R);
            mkpi_pwave_2_norm = fcnc::int_mkpi_bw_kstar_amp_squared(from, to, q2, gammakstar, mkstar, R);
            int_mkpi_re_swavepwave = fcnc::int_mkpi_bw_kstar_amp_bar_kstarzero_isobar_amp(from, to, q2, gammakstar, mkstar, asphase, f800mag, f800phase, gammaf800, mf800, gammakstarplus, mkstarplus, R).real();
            int_mkpi_im_swavepwave = fcnc::int_mkpi_bw_kstar_amp_bar_kstarzero_isobar_amp(from, to, q2, gammakstar, mkstar, asphase, f800mag, f800phase, gammaf800, mf800, gammakstarplus, mkstarplus, R).imag();
        }
        else{
            mkpi_swave_2_norm = fcnc::int_mkpi_kstarzero_lass_amp_squared(from, to, q2, asphase, a, r, gammakstarplus, mkstarplus, R);
            mkpi_pwave_2_norm = fcnc::int_mkpi_bw_kstar_amp_squared(from, to, q2, gammakstar, mkstar, R);
            int_mkpi_re_swavepwave = fcnc::int_mkpi_bw_kstar_amp_bar_kstarzero_lass_amp(from, to, q2, gammakstar, mkstar, asphase, a, r, gammakstarplus, mkstarplus, R).real();
            int_mkpi_im_swavepwave = fcnc::int_mkpi_bw_kstar_amp_bar_kstarzero_lass_amp(from, to, q2, gammakstar, mkstar, asphase, a, r, gammakstarplus, mkstarplus, R).imag();
        }
    }
    else spdlog::debug("Skipping generation of m(Kpi)!");

    //get all P- and S-wave angular observables
    const double j1s = params->J1s();
    const double j1c = params->J1c();
    const double j2s = params->J2s();
    const double j2c = params->J2c();
    const double j3 = params->J3();
    const double j4 = params->J4();
    const double j5 = params->J5();
    const double j6s = params->J6s();
    const double j6c = params->J6c();
    const double j7 = params->J7();
    const double j8 = params->J8();
    const double j9 = params->J9();
    const double fs = params->FS();
    const double js1 = params->SS1();
    const double js2 = params->SS2();
    const double js3 = params->SS3();
    const double js4 = params->SS4();
    const double js5 = params->SS5();

    unsigned int n_gen = 0;  //Number of generated events
    unsigned int n_neg = 0;  //Number of events with negative angular probability
    unsigned int n_rej = 0;  //Number of rejected events in the generation
    unsigned int n_exc = 0;  //Number of accepted events

    //main event loop
    for (unsigned int ev = 0; ev < nevents; ev++){
        event meas;
        meas.cp_conjugate = rnd->Rndm() < 0.5;//no overall cp asymmetry, so 50/50 chance of a 1 and 0

        //give a progress report
        if (ev % 1000 == 0 || ev == nevents - 1){
            std::lock_guard<std::mutex> lock(generator_mutex);
            spdlog::debug("\rThread {0:d} generating event no. {1:d} \r", thread_id , ev + 1 );
            spdlog::default_logger_raw()->flush();
        }

        //determine if event is signal or background
        bool is_signal = (rnd->Rndm() < effective_f_sig);
        if(opts->generate_only_bkg) is_signal = false;
        if(is_signal) meas.event_type = 0;
        else          meas.event_type = 1;

        //m_res_1 should be the ratio between the two CBs
        //In case of 2-tailed CB, m_res_1 is always 1
        if(rnd->Rndm() < params->m_res_1()) meas.sigma_m = params->m_sigma_1();
        else                                meas.sigma_m = params->m_sigma_2();

        //generate mass (hit and miss), no real need to speed this up for now
        bool finished_mass = false; //Means done with generating the mass
        while(!finished_mass){
            //get a random mass value
            meas.m = genValueFromPar(&params->m_b, rnd); //get random number between the range of B mass parameter
            //determine probability from signal PDF
            if(is_signal){
                event mmax;
                mmax.m = params->m_b();
                //The probability is calculated in bu2kstarmumu_pdf, and it takes into account twotailedCB if opts.twotailedcrystalball is true
                if(rnd->Rndm() < prob->m_sig_prob(params, meas)/prob->m_sig_prob(params, mmax)){
                    finished_mass = true;
                }
            }
            //or determine background probability
            else{
                event mmax;
                if(opts->fit_lambda){
                    if (params->m_lambda() < 0.0) mmax.m = params->m_b.get_min();
                    else                          assert(params->m_lambda() < 0.0);
                }
                else{
                    if (params->m_tau() > 0.0) mmax.m = params->m_b.get_min();
                    else                       mmax.m = params->m_b.get_max();
                }
                double m_probmax = prob->m_bkg_prob(params, mmax);
                //systematic 11: remove a gaussian shaped hole in the upper mass sideband:
                if(opts->systematic == 11){ //TODO
                    double gaussian_width = 30.;
                    double gaussian_mean = 5520.;
                    double gaussian_frac = 0.02;
                    double gaussian_hole = TMath::Exp(-(meas.m - gaussian_mean)*(meas.m - gaussian_mean)/(2.*gaussian_width*gaussian_width))/TMath::Sqrt(2.*TMath::Pi()*gaussian_width*gaussian_width);
                    double exp_norm = 1./params->m_lambda()*(exp(opts->m_max*params->m_lambda()) - exp(opts->m_min*params->m_lambda()));
                    gaussian_hole *= gaussian_frac * exp_norm;

                    if (rnd->Rndm()*m_probmax < prob->m_bkg_prob(params, meas) - gaussian_hole) finished_mass = true;
                }
                else{
                    if (rnd->Rndm()*m_probmax < prob->m_bkg_prob(params, meas)) finished_mass = true;
                    //If your random shot on the y-axis is smaller than the wanted distribution/comparison function, keep the event!
                }
            }
        }//mass loop end

        //generate angles
        if(is_signal){//generate signal
            double max_gen = opts->generate_mkpi ? 10.0 : 2.0;
            //This is the 'upper limit' or 'comparison funcition', pretty much the top of your generation. It always has to be bigger than dGamma/d(ctl)d(ctk)d(phi)!!!

            bool finished = false;
            while (!finished){

                n_gen++;
                double result = 1.0;

                //get q2 value
                meas.q2 = params->eff_q2();
                if (opts->weighted_fit){ //TODO eh?
                    assert(opts->q2_min > 0.0);
                    assert(opts->q2_max > opts->q2_min);
                    meas.q2 = genValueFromRange(opts->q2_min, opts->q2_max, rnd);
                }

                //get angles
                meas.costhetak = genValueFromRange(opts->ctk_min, opts->ctk_max, rnd);
                meas.costhetal = genValueFromRange(opts->ctl_min, opts->ctl_max, rnd);
                meas.phi       = genValueFromRange(opts->phi_min, opts->phi_max, rnd);

                double f1=0.0, f2=0.0, f3=0.0, f4 =0.0, f5 =0.0, f6 =0.0,
                        f7=0.0, f8=0.0, f9=0.0, f10=0.0, f11=0.0, f12=0.0;
                if(opts->full_angular){
                    prob->fj(meas.costhetal, meas.costhetak, meas.phi,
                             f1, f2, f3, f4,  f5,  f6,
                             f7, f8, f9, f10, f11, f12);
                }
                else{
                    prob->folded_fj(meas.costhetal, meas.costhetak, meas.phi,
                                    f1, f2, f3, f4,  f5,  f6,
                                    f7, f8, f9, f10, f11, f12);
                }

                //calculate P-wave probability
                result = f1*j1s + f2*j1c + f3*j2s + f4*j2c + f5*j3  + f6*j4
                        + f7*j5  + f8*j6s + f9*j6c + f10*j7 + f11*j8 + f12*j9;

                if(opts->fit_asymmetries && meas.cp_conjugate){
                    result = f1*j1s + f2*j1c + f3*j2s + f4*j2c - f5*j3  - f6*j4
                            - f7*j5  - f8*j6s - f9*j6c - f10*j7 - f11*j8 - f12*j9;
                }

                //generate m(Kpi)
                double prob_mkpi_pwave=0.0, prob_mkpi_swave=0.0;
                std::complex<double> prob_mkpi_swavepwave={0.0, 0.0};
                if (opts->generate_mkpi){ //TODO: here it is!
                    assert(mkpi_min > 0.0);
                    assert(mkpi_max > mkpi_min);
                    meas.mkpi =  genValueFromRange(mkpi_min, mkpi_max, rnd)*1.0e+3;
                    if(opts->simple_mkpi) prob_mkpi_pwave = fcnc::mkpi_simple_bw_kstar_amp_squared(meas.mkpi/1.0e+3, meas.q2, gammakstar, mkstar);
                    else if(opts->isobar) prob_mkpi_pwave = fcnc::mkpi_bw_kstar_amp_squared(meas.mkpi/1.0e+3, meas.q2, gammakstar, mkstar, R);
                    else                  prob_mkpi_pwave = fcnc::mkpi_bw_kstar_amp_squared(meas.mkpi/1.0e+3, meas.q2, gammakstar, mkstar, R);
                    result *= prob_mkpi_pwave/mkpi_pwave_2_norm;
                }
                else{
                    meas.mkpi = PDGMASS_K_STAR_PLUS;
                }

                //calculate S-wave probability
                double result_swave = 0.0;
                if(swave){
                    //P-wave gets factor (1-FS)
                    result *= (1.0-fs);
                    double fs1=0.0, fs2=0.0, fs3=0.0, fs4=0.0, fs5=0.0, fs6=0.0;
                    if(opts->full_angular){
                        prob->swave_fj(meas.costhetal, meas.costhetak, meas.phi,
                                       fs1, fs2, fs3, fs4, fs5, fs6);
                    }
                    else{
                        prob->folded_swave_fj(meas.costhetal, meas.costhetak, meas.phi,
                                              fs1, fs2, fs3, fs4, fs5, fs6);
                    }

                    result_swave = fs*fs1 + js1*fs2 + js2*fs3 + js3*fs4 + js4*fs5 + js5*fs6;
                    if(opts->generate_mkpi){
                        if(opts->simple_mkpi){
                            prob_mkpi_swave = fcnc::mkpi_simple_kstarzero_amp_squared(meas.mkpi/1.0e+3, meas.q2,
                                                                                      f800mag, f800phase, gammaf800, mf800, gammakstarplus, mkstarplus);
                            prob_mkpi_swavepwave = fcnc::mkpi_simple_bw_kstar_amp_bar_kstarzero_amp(meas.mkpi/1.0e+3, meas.q2,
                                                                                                    gammakstar, mkstar, asphase,
                                                                                                    f800mag, f800phase, gammaf800, mf800, gammakstarplus, mkstarplus);
                        }
                        else if(opts->isobar){
                            prob_mkpi_swave = fcnc::mkpi_kstarzero_isobar_amp_squared(meas.mkpi/1.0e+3, meas.q2,
                                                                                      f800mag, f800phase, gammaf800, mf800, gammakstarplus, mkstarplus, R);
                            prob_mkpi_swavepwave = fcnc::mkpi_bw_kstar_amp_bar_kstarzero_isobar_amp(meas.mkpi/1.0e+3, meas.q2,
                                                                                                    gammakstar, mkstar, asphase,
                                                                                                    f800mag, f800phase, gammaf800, mf800, gammakstarplus, mkstarplus, R);
                        }
                        else{
                            prob_mkpi_swave = fcnc::mkpi_kstarzero_lass_amp_squared(meas.mkpi/1.0e+3, meas.q2, asphase, a, r, gammakstarplus, mkstarplus, R);
                            prob_mkpi_swavepwave = fcnc::mkpi_bw_kstar_amp_bar_kstarzero_lass_amp(meas.mkpi/1.0e+3, meas.q2, gammakstar, mkstar, asphase, a, r, gammakstarplus, mkstarplus, R);
                        }
                        result_swave =  fs*fs1*prob_mkpi_swave/mkpi_swave_2_norm
                                +  (js1*fs2 + js2*fs3 + js3*fs4)*prob_mkpi_swavepwave.real()/int_mkpi_re_swavepwave
                                +  (js4*fs5 + js5*fs6)*prob_mkpi_swavepwave.imag()/int_mkpi_im_swavepwave;
                    }
                    //add S-wave probability to P-wave probability
                    result += result_swave;
                }

                //safety check on nan or inf
                if(std::isnan(result) || std::isinf(result)){
                    spdlog::critical("The calculated probability is = {0:f}", result);
                    assert(!std::isnan(result));
                    assert(!std::isinf(result));
                }

                //safety check to avoid negative probabilities at toy generation
                if(result < 0.0){
                    n_neg++;
                    continue;
                }

                //determine angular acceptance
                if (opts->use_angular_acc || opts->weighted_fit){
                    result *=prob->get_ang_eff(opts, meas, false);
                }

                if (result > max_gen){
                    double eff = prob->get_ang_eff(opts, meas, false);
                    //Yes, the eff is calculated again, but it doesn't cost any time as it crashes here anyway
                    spdlog::error("Determined probability is too large: {0:e}. With efficiency: {1:f}", result, eff);
                    spdlog::error("swave={0:f}\tpwave={1:f}", result_swave, result/eff - result_swave);
                    spdlog::critical("If you are seing this, increase the max_gen!!!");
                    assert(0);
                }
                else if (rnd->Rndm() < result/max_gen){ //If your random shot on the y-axis is smaller than the wanted distribution/comparison function, keep the event!
                    finished = true;
                    n_exc++;
                }
                else{
                    n_rej++;
                }
            };
        }
        else{//generate bkg
            bool finished = false;
            double max_gen = 25.;//This is the 'upper limit' or 'comparison funcition', pretty much the top of your generation. It always has to be bigger than dGamma/d(ctl)d(ctk)d(phi)!!!

            while(!finished){

                //generate angles
                meas.costhetak = genValueFromRange(opts->ctk_min, opts->ctk_max, rnd);
                meas.costhetal = genValueFromRange(opts->ctl_min, opts->ctl_max, rnd);
                meas.phi       = genValueFromRange(opts->phi_min, opts->phi_max, rnd);

                //generate q2
                meas.q2 = params->eff_q2();
                if(opts->weighted_fit){
                    assert(opts->q2_min > 0.0);
                    assert(opts->q2_max > opts->q2_min);
                    meas.q2 = genValueFromRange(opts->q2_min, opts->q2_max, rnd);
                }

                double result = 1.0;
                if(!opts->flat_bkg){
                    //use this (3 * 1D)*eff generation instead of 4D(Omega, eff) generation from prob->angular_bkg_prob
                    std::vector<double> ch_ctl = init_ch_ctl(params);
                    std::vector<double> ch_ctk = init_ch_ctk(params);
                    std::vector<double> ch_phi = init_ch_phi(params);

                    std::vector<double> poly_ctl(ch_ctl.size(), 0.0);
                    std::vector<double> poly_ctk(ch_ctk.size(), 0.0);
                    std::vector<double> poly_phi(ch_phi.size(), 0.0);

                    std::vector<double> temp_ctl(ch_ctl.size(), 0.0);
                    std::vector<double> temp_ctk(ch_ctk.size(), 0.0);
                    std::vector<double> temp_phi(ch_phi.size(), 0.0);

                    chebyshev_to_poly(ch_ctl, temp_ctl);
                    chebyshev_to_poly(ch_ctk, temp_ctk);
                    chebyshev_to_poly(ch_phi, temp_phi);

                    for (unsigned int i=0; i<temp_ctl.size(); i++) poly_ctl.at(i) = temp_ctl.at(i);
                    for (unsigned int i=0; i<temp_ctk.size(); i++) poly_ctk.at(i) = temp_ctk.at(i);
                    correct_poly(temp_phi, poly_phi, -TMath::Pi(), +TMath::Pi());

                    double p_costhetal = 0.0;
                    for(unsigned int i=0; i<ch_ctl.size(); i++) p_costhetal += poly_ctl.at(i)*pow(meas.costhetal,i);
                    double p_costhetak = 0.0;
                    for(unsigned int i=0; i<ch_ctk.size(); i++) p_costhetak += poly_ctk.at(i)*pow(meas.costhetak,i);
                    double p_phi = 0.0;
                    for(unsigned int i=0; i<ch_phi.size(); i++) p_phi += poly_phi.at(i)*pow(meas.phi,i);

                    result = p_costhetal*p_costhetak*p_phi;///prob->fnorm_bkg;
                }

                //efficiency correction in generation:

                if((opts->use_angular_acc && opts->use_weighted_bkg) || opts->weighted_fit){//applies efficiency to bkg for weighted fit
                    result *= prob->get_ang_eff(opts, meas, false);
                }

                if(result > max_gen){
                    spdlog::error("The probability is too large: {0:f}", result);
                    spdlog::error("It has to be smaller than: {0:f}", max_gen);
                    spdlog::error("The normalization of the PDF is: {0:f}", prob->fnorm_bkg);
                    spdlog::critical("If you are seing this, increase the max_gen!!!");
                    assert(0);
                }
                else if(max_gen*rnd->Rndm() < result) finished = true;
            }//end angular background generation

            finished = false;
            if(opts->generate_mkpi){
                double mkpi_min = opts->mkpi_min/1.0e+3;
                double mkpi_max = opts->mkpi_max/1.0e+3;
                std::vector<double> ch_mkpi(5, 0.0);
                ch_mkpi.at(0) = params->cbkgmkpi0();
                ch_mkpi.at(1) = params->cbkgmkpi1();
                ch_mkpi.at(2) = params->cbkgmkpi2();
                ch_mkpi.at(3) = params->cbkgmkpi3();
                ch_mkpi.at(4) = params->cbkgmkpi4();
                std::vector<double> poly_mkpi(ch_mkpi.size(), 0.0);
                chebyshev_to_poly(ch_mkpi, poly_mkpi);
                std::vector<double> poly_corr_mkpi(ch_mkpi.size(), 0.0);
                correct_poly(poly_mkpi, poly_corr_mkpi, mkpi_min, mkpi_max);

                while(!finished){
                    double mkpi = mkpi_min+rnd->Rndm()*(mkpi_max-mkpi_min);
                    double probbkg = 0.0;
                    if (opts->mkpi_threshold){
                        probbkg = fcnc::threshold(mkpi, 0.633, params->nthreshold());
                    }
                    else{
                        for(unsigned int i=0; i<ch_mkpi.size(); i++) probbkg += poly_corr_mkpi.at(i)*pow(mkpi,i);
                    }
                    double probmax = 10.0;
                    assert(probbkg < probmax);
                    if(rnd->Rndm()*probmax < probbkg){
                        meas.mkpi = mkpi*1.0e+3;
                        finished = true;
                    }
                }
            }//end m(Kpi) background generation
            else{
                meas.mkpi = PDGMASS_K_STAR_PLUS;
            }
        }

        //save angles also in backup/full angular variables:
        meas.costhetal_fa = meas.costhetal;
        meas.costhetak_fa = meas.costhetak;
        meas.phi_fa = meas.phi;
        meas.year = opts->year;

        eventVec->push_back(meas);

    }//end for all n events that are to be constructed
    spdlog::debug("Thread #{0:d}: Angular signal events: n_gen={1:d}\tn_neg={2:d}\tn_rej={3:d}\tn_exc={4:d}", thread_id, n_gen, n_neg, n_rej, n_exc);
    delete rnd;
}
