//Renata Kopecna


#include <pulls.hh>

#include <string>
#include <sstream>
#include <iostream>

#include "folder.hh"
#include <fitter.hh>
#include <bu2kstarmumu_plotter.hh>
#include <bu2kstarmumu_generator.hh>
#include <bu2kstarmumu_pdf.hh>
#include <bu2kstarmumu_parameters.hh>
#include <paths.hh>
#include <helpers.hh>
#include <event.hh>

#include <TTree.h>

#include <spdlog.h>

//////////////////////////
// INIT VALUES FOR FITS //
//////////////////////////

  //Do they change anywhere or can I set them in a header?
//TODO replace by the parameters->init
double f_m_B[4]      = {5282.32, 5282.08, 5278.57, 5279.48};
double f_lambda1[4]  = {-0.00386, -0.00610, -0.00362, -0.00578};
double f_lambda2[4]  = {-0.03320, -0.00120, -0.04600, -0.00100};
double f_alpha_1[4]  = { 1.7910,  1.8000,  1.6510,  1.5400};
double f_alpha_2[4]  = { 2.3200,  2.4300,  2.0100,  2.2500};
double f_n_1[4]      = { 1.4400,  1.4900,  1.6800,  1.6100};
double f_n_2[4]      = { 2.8900,  3.7000,  3.8700,  3.7000};
double f_sigma_1[4]  = {14.8300, 15.4800, 14.8400, 14.4};
double f_sigma_2[4]  = {14.0700, 13.1000, 15.6200, 13.3};
double f_f_CB[4]     = { 0.6360,  0.3600,  0.3040,  0.6200};
double f_m_scale[4]  = { 1.1814,  1.0648,  1.1547,  1.1168};


//TODO posibly to be deleted so ignoring for now!
void do_pulls(fcnc::options opts, bool doToyPulls, bool doMCPulls, int job_id, double SwaveFraction, UInt_t nPulls, Double_t nToyEvents,  UInt_t nMCEvents ,std::string angularsuffix, bool Fit1bin, bool Fit2bins, bool FitAllbins ){ //TODO: change to int and return some errors possibly
    //TODO: figure out jobs_to_do and job_id, but there is no manipulation with it, so so far so good
    //Why is it not opts.job_id?
    if(doToyPulls)std::cout << "[INFO]\t\tGetting pull distributions and test fits from Toy Events!" << std::endl;

    //disable eps output if job_id is larger than -1. This is needed for HD cluster jobs:
    if(job_id > -1)opts.write_eps = false;

    //have extra option for the usage of binning:
    const bool UseBinnedFit = false;

    opts.fit_fl = false;
    opts.fit_afb = false;
    opts.only_angles = false;
    opts.only_Bmass = false;
    if(opts.only_Bmass)  opts.only_angles = false;

    opts.use_angular_acc = true;
    opts.use_weighted_bkg = true;
    opts.fit_full_angular_bkg = true;
    opts.always_generate_full_angular = true;
    opts.weighted_fit = true;

    opts.swave = true;
    opts.shift_lh = false;

    opts.squared_hesse = true;
    opts.minos_errors = false;

    opts.generate_mkpi = true;
    opts.fit_mkpi = false;
    opts.use_mkpi = true;
    opts.flat_bkg = false;

    //change the label of the produced plots:
    opts.plot_label = "LHCb toys";

    //SM
    std::vector<double> s1s, s1c, s2s, s2c, s3, s4, s5, s6s, s6c, s7, s8, s9;
    std::vector<double> f_signal, f_bckgnd, f_bckcoeff;
    double f_sig_norm = 0.0;
    std::vector<UInt_t> events_per_bin;

    //determine number of bins:
    const UInt_t nBins = UseBinnedFit ? opts.TheQ2binsmin.size() : 1;
    if(!UseBinnedFit){
        opts.TheQ2binsmin = {8.68};
        opts.TheQ2binsmax = {10.09};
    }

    //create TTree to save pull values
    double treeValue, treeStart, treeError, treeErrorUp, treeErrorDown;
    UInt_t treeBin, treeFitNumber, treeFitResult, treeIndex;

    TTree * T = nullptr;
    if(doToyPulls || doMCPulls){
        T = new TTree("PullTree", "PullTree");
        T->Branch("v", &treeValue, "v/D");
        T->Branch("i", &treeIndex, "i/i");
        T->Branch("s", &treeStart, "s/D");
        T->Branch("e", &treeError, "e/D");
        T->Branch("eu", &treeErrorUp, "eu/D");
        T->Branch("ed", &treeErrorDown, "ed/D");
        T->Branch("b", &treeBin  , "b/i");
        T->Branch("f", &treeFitNumber, "f/i");
        T->Branch("r", &treeFitResult, "r/i");
        T->Branch("fs", &SwaveFraction, "fs/D");
    }

    //TODO: fix the name of the genLvl root file
    basic_params params_genLvl = basic_params();
    params_genLvl.Run = 2;
    params_genLvl.year = 2017;
    params_genLvl.nBins = nBins;
    std::string genPath = final_result_name_genLvlMC(params_genLvl, nBins, false, true, false);

    s1s = load_param_values_into_vector("S1s",genPath);
    s3  = load_param_values_into_vector("S3", genPath);
    s4  = load_param_values_into_vector("S4", genPath);
    s5  = load_param_values_into_vector("S5", genPath);
    s6s = load_param_values_into_vector("S6s",genPath);
    s7  = load_param_values_into_vector("S7", genPath);
    s8  = load_param_values_into_vector("S8", genPath);
    s9  = load_param_values_into_vector("S9", genPath);

    if(UseBinnedFit){
        if(nBins == 8){
            f_signal   = {101.,  61.,  62.,  96., 125., 124., 129., 69.};
            f_bckgnd   = {131., 206., 279., 358., 358., 247., 158., 93.};
            f_bckcoeff = {-0.0052, -0.0084, -0.0088, -0.0031, -0.0028, -0.0158, -0.0079, -0.0058}; //for notation: exp(lambda * x)

            //results from Signal channel generator level MC fits:
            s1c = { 0.2450, 0.6996, 0.7993, 0.7384, 0.6414, 0.4201, 0.3513, 0.3344};
            s2s = { 0.1448, 0.0749, 0.0514, 0.0662, 0.0900, 0.1449, 0.1620, 0.1662};
            s2c = {-0.1953,-0.6639,-0.7774,-0.7254,-0.6334,-0.4171,-0.3497,-0.3332};
            s3  = { 0.0000, 0.0020,-0.0080,-0.0120,-0.0240,-0.0620,-0.1520,-0.2330};
            s4  = { 0.0870, 0.0040,-0.1080,-0.1950,-0.2510,-0.2770,-0.2900,-0.3040};
            s5  = { 0.2440, 0.0950,-0.1590,-0.3050,-0.4140,-0.4210,-0.3400,-0.2490};
            s6s = {-0.1270,-0.2070,-0.0860, 0.0950, 0.3020, 0.5140, 0.5570, 0.4700};
            s6c = { 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000};
            s7  = {-0.0060,-0.0100,-0.0030, 0.0020, 0.0030,-0.0020,-0.0020, 0.0050};
            s8  = { 0.0010,-0.0030,-0.0040,-0.0030,-0.0010, 0.0030,-0.0010, 0.0030};
            s9  = {-0.0020, 0.0000, 0.0020, 0.0040,-0.0020,-0.0010,-0.0010, 0.0010};

        }
        else if(nBins == 4){
            spdlog::error("YOU NEED TO UPDATE THESE FIRTS!");
            assert(0);
            f_signal   = {101.,  61.,  62.,  96.};
            f_bckgnd   = {131., 206., 279., 358.};
            f_bckcoeff = {-0.0052, -0.0084, -0.0088, -0.0031}; //for notation: exp(lambda * x)

            //results from Signal channel generator level MC fits:
            s1s = { 0.5120, 0.1800, 0.1370, 0.1930};
            s1c = { 0.2450, 0.6996, 0.7993, 0.7384};
            s2s = { 0.1448, 0.0749, 0.0514, 0.0662};
            s2c = {-0.1953,-0.6639,-0.7774,-0.7254};
            s3  = { 0.0000, 0.0020,-0.0080,-0.0120};
            s4  = { 0.0870, 0.0040,-0.1080,-0.1950};
            s5  = { 0.2440, 0.0950,-0.1590,-0.3050};
            s6s = {-0.1270,-0.2070,-0.0860, 0.0950};
            s6c = { 0.0000, 0.0000, 0.0000, 0.0000};
            s7  = {-0.0060,-0.0100,-0.0030, 0.0020};
            s8  = { 0.0010,-0.0030,-0.0040,-0.0030};
            s9  = {-0.0020, 0.0000, 0.0020, 0.0040};
        }
        else if(nBins == 2){
            f_signal   = { 61.+ 62.+ 96.+125., 129.+69.};
            f_bckgnd   = {206.+279.+358.+358., 158.+93.};
            f_bckcoeff = {-0.0064, -0.0065};

            s1s = { 0.1959,  0.4917};
            s1c = { 0.7441,  0.3447};
            s2s = { 0.0646,  0.1636};
            s2c = {-0.7217, -0.3432};
            s3  = {-0.0042, -0.1747};
            s4  = { 0.1033,  0.2929};
            s5  = {-0.1458, -0.3157};
            s6s = { 0.0652, -0.5338};
            s6c = { 0.0000,  0.0000};
            s7  = { 0.0340, -0.0000};
            s8  = {-0.0115, -0.0002};
            s9  = { 0.0003,  0.0003};
        }
        else {
            spdlog::error("No SM values given for binning scheme with {0:d} q2 bins. Exit", nBins);
            assert(0);
        }
    }
    else{
        f_signal   = {0.5};
        f_bckgnd   = {0.5};
        f_bckcoeff = {-2.835e-3}; //for notation: exp(lambda * x)
        s1s = { 0.2372};
        s1c = { 0.6863};
        s2s = { 0.0};
        s2c = {-0.6760};
        s3  = {-0.0107};
        s4  = { 0.2163};
        s5  = {-0.3806};
        s6s = {-0.1957};
        s6c = { 0.0};
        s7  = { 0.0290};
        s8  = {-0.0113};
        s9  = { 0.0003};
    }

    assert(f_signal.size() == f_bckgnd.size());
    f_sig_norm = std::accumulate(f_signal.begin(), f_signal.end(), 0.0)
            + std::accumulate(f_bckgnd.begin(), f_bckgnd.end(), 0.0);

    for(UInt_t i = 0; i < f_signal.size(); i++){
        f_signal.at(i) /= f_sig_norm;
        f_bckgnd.at(i) /= f_sig_norm;
    }

    fcnc::folder fldr(&opts);
    fcnc::fitter f(&opts);
    fcnc::bu2kstarmumu_plotter thePlotter(&opts);

    std::vector< std::vector<fcnc::bu2kstarmumu_parameters*> >theParameters;
    std::vector< std::vector<fcnc::bu2kstarmumu_pdf*> >       thePDFs;
    for(UInt_t c = 0; c < nBins; c++){
        //add vector of params to matrix:
        std::vector<fcnc::bu2kstarmumu_parameters*> params_per_bin;
        params_per_bin.clear();
        theParameters.push_back(params_per_bin);

        //add vector of pdfs to matrix:
        std::vector<fcnc::bu2kstarmumu_pdf*> pdf_per_bin;
        pdf_per_bin.clear();
        thePDFs.push_back(pdf_per_bin);
    }
    assert(theParameters.size() == nBins);
    assert(thePDFs.size()       == nBins);

    //pointers to current parameter, pdf and generator:
    fcnc::bu2kstarmumu_parameters * params[nBins];
    fcnc::bu2kstarmumu_pdf        * prob[nBins];
    fcnc::bu2kstarmumu_generator  * gen    = nullptr;

    std::vector<int> fitresults;
    std::vector<fcnc::event>selection;

    //number of fits is either determined by the larger number of pulls or fits requested
    UInt_t nFits = nPulls;

    //vectors to save values for pull plots:
    std::vector<int>var_indexs;
    std::vector< std::vector< std::vector<double> > >pull_values;
    std::vector< std::vector< std::vector<double> > >pull_errors;
    std::vector< std::vector<double> >pull_starts;


    //initialize parameters, events and fitter for every FitTest:
    for(UInt_t n = 0; n < nFits; n++){

        //if binning is used, loop over the number of q2 bins, else only run loop once!
        for(UInt_t b = 0; b < nBins; b++){

            std::cout << std::endl;
            std::cout << "****************************************" << std::endl;
            std::cout << "[FIT]\t\tStarting fit >> " << n << " <<" << std::endl;
            std::cout << "[FIT]\t\tIn bin >> " << b << " <<" << std::endl;
            std::cout << "****************************************" << std::endl;
            std::cout << std::endl;

            opts.swave = true;

            opts.q2_min = UseBinnedFit ? opts.TheQ2binsmin.at(b) : opts.TheQ2binsmin.front();
            opts.q2_max = UseBinnedFit ? opts.TheQ2binsmax.at(b) : opts.TheQ2binsmax.back();

            if(n == 0){
                //initiate parameters for mass fit
                params[b] =  new fcnc::bu2kstarmumu_parameters(&opts);
                params[b]->f_sig.init(f_signal.at(b)/(f_signal.at(b)+f_bckgnd.at(b)), 0.0, 1.0, 0.01);
                if(params[b]->f_sig.get_value() == 0.0)
                    params[b]->m_b.init(PDGMASS_B, B_MASS_LOW, B_MASS_HIGH, 0.0);
                else
                    params[b]->m_b.init(PDGMASS_B, B_MASS_LOW, B_MASS_HIGH, 0.01);
                if(opts.twotailedcrystalball)
                    params[b]->m_res_1.init(1.0, 0.0, 1.0, 0.0);
                else
                    params[b]->m_res_1.init(0.391, 0.0, 1.0, 0.0);// 0.01);
                params[b]->fm_tau.init(1.0, 0.0, 1.0, 0.0);
                if(opts.fit_lambda){
                    if(params[b]->f_sig.get_value() == 1.0)
                        params[b]->m_lambda.init(f_bckcoeff.at(b), -1.0e-1, -1.0e-6, 0);
                    else
                        params[b]->m_lambda.init(f_bckcoeff.at(b), -1.0e-1, -1.0e-6, 1.0e-5);
                    params[b]->m_lambda_2.init(-1.664e-3, -1.0e-1, -1.0e-6, 0.0);
                    params[b]->m_tau.init( 1./2.835e-3, 100.0, 1.0e+4, 0.0);
                    params[b]->m_tau_2.init(1./1.664e-3, 0.0, 1.0e+4, 0.0);
                }
                else{
                    params[b]->m_lambda.init(-1.0e-3, -1.0e-1, -1.0e-6, 0.0);
                    params[b]->m_lambda_2.init(-1.5e-3, -1.0e-1, -1.0e-6, 0.0);
                    params[b]->m_tau.init(-1./f_bckcoeff.at(b), 100.0, 1.0e+4, 1.0);
                    params[b]->m_tau_2.init(6.01e+2, 0.0, 1.0e+4, 0.0);
                }
                if(opts.twotailedcrystalball)
                    params[b]->m_sigma_1.init(f_sigma_1[n], 5.0, 200.0, 0.0);//, 0.1);
                else{
                    params[b]->m_sigma_1.init(f_sigma_1[n], 5.0, 200.0, 0.0);//, 0.1);
                    params[b]->m_sigma_2.init(f_sigma_2[n], 5.0, 200.0, 0.0);//, 0.1);
                }
                params[b]->alpha_1.init(f_alpha_1[n], 0.1, 10.0, 0.0);//, 0.15);
                params[b]->alpha_2.init(f_alpha_2[n], 0.1, 10.0, 0.0);//, 0.20);
                params[b]->n_1.init(f_n_1[n], 0.1, 15.0, 0.0);//, 2.0);
                params[b]->n_2.init(f_n_2[n], 0.1, 10.0, 0.0);//, 0.16);
                params[b]->m_scale.init(f_m_scale[n], 0.0, 2.0, 0.0);

                params[b]->load_only_Bmass_param_values("fitresult_SignalFit_MC_SimultaneousFit_2Dfit_bin"+std::to_string(b)+"_pdf"+std::to_string(n)+".txt");

                //S-wave
                double sWaveStepSize = 0.01;
                params[b]->FS.init(SwaveFraction, 0.0, 1.0, sWaveStepSize);

                //B0 fit results
                params[b]->SS1.init(-0.231, -1.0, 1.0, (opts.full_angular || opts.folding != 4 ? sWaveStepSize : 0.0));
                params[b]->SS2.init( 0.023, -1.0, 1.0, (opts.full_angular || opts.folding == 1 ? sWaveStepSize : 0.0));
                params[b]->SS3.init( 0.003, -1.0, 1.0, (opts.full_angular || opts.folding == 2 ? sWaveStepSize : 0.0));
                params[b]->SS4.init( 0.001, -1.0, 1.0, (opts.full_angular || opts.folding > 2 ? sWaveStepSize : 0.0));
                params[b]->SS5.init(-0.068, -1.0, 1.0, (opts.full_angular ? sWaveStepSize : 0.0));

                if(opts.fit_mkpi || opts.use_mkpi){
                    params[b]->gammakstar.init(0.0503, 0.01, 0.8, 0.001); //PDG value
                    params[b]->mkstar.init(PDGMASS_K_STAR_PLUS / 1000., 0.7, 1.2, 0.001); //PDG value
                    params[b]->gammakstarplus.init(0.236, 0.1, 1.0, 0.0); //PDG value
                    params[b]->mkstarplus.init(1.41, 1.2, 1.6, 0.0); //PDG value
                    params[b]->asphase.init(TMath::Pi(), 0.0, 2.0*TMath::Pi(), 0.0);
                    params[b]->a.init(1.95, 0.001, 20.0, 0.0);
                    params[b]->r.init(1.78, 0.001, 10.0, 0.0);
                    //background parameter
                    if(params[b]->f_sig() != 1.0){
                        params[b]->cbkgmkpi0.init(1.0, -1.0, 1.0, 0.0);
                        params[b]->cbkgmkpi1.init(0.0, -1.0, 1.0, 0.01);
                        params[b]->cbkgmkpi2.init(0.0, -1.0, 1.0, 0.0);
                        params[b]->cbkgmkpi3.init(0.0, -1.0, 1.0, 0.0);
                        params[b]->cbkgmkpi4.init(0.0, -1.0, 1.0, 0.0);
                        opts.mkpi_threshold = false;
                        params[b]->nthreshold.init(1.5, 0.0, 15.0, 0.0);
                    }
                }

                //define center of q2bin as effective q2:
                params[b]->eff_q2.init((UseBinnedFit ? 0.5*(opts.TheQ2binsmin.at(b)+opts.TheQ2binsmax.at(b)) : 6.0), opts.TheQ2binsmin.front(), opts.TheQ2binsmax.back(), 0.0);

                //if only mass fit: do not active angles in fit, i.e. set stepsize to 0.0
                double angleStepSize = opts.only_Bmass ? 0.0 : 0.001;
                if(params[b]->f_sig.get_value() == 0.0)
                    angleStepSize = 0.0;
                if(opts.fit_fl)
                    //params[b]->Fl.init(3.0/4.0*(s1c[b]-s2c[b]/3.0), 0.0, 1.0, angleStepSize);
                    params[b]->Fl.init(1.0-4.0/3.0*s1s[b], 0.0, 1.0, angleStepSize);
                else
                    params[b]->S1s.init(s1s[b], -1.0, 1.0, angleStepSize);
                params[b]->S3.init(s3[b], -1.0, 1.0, angleStepSize/10.);
                params[b]->S4.init(s4[b], -1.0, 1.0, (opts.full_angular || opts.folding == 1 ? angleStepSize : 0.0));
                params[b]->S5.init(s5[b], -1.0, 1.0, (opts.full_angular || opts.folding == 2 ? angleStepSize : 0.0));
                if(opts.fit_afb)
                    params[b]->Afb.init(3.0/4.0*s6s[b], -0.75, 0.75, (opts.full_angular || opts.folding == 0 ? angleStepSize: 0.0));
                else
                    params[b]->S6s.init(s6s[b], -1.0, 1.0, (opts.full_angular || opts.folding == 0 ? angleStepSize : 0.0));
                params[b]->S7.init(s7[b], -1.0, 1.0, (opts.full_angular || opts.folding == 3 ? angleStepSize/10. : 0.0));
                params[b]->S8.init(s8[b], -1.0, 1.0, (opts.full_angular || opts.folding == 4 ? angleStepSize/10. : 0.0));
                params[b]->S9.init(s9[b], -1.0, 1.0, (opts.full_angular || opts.folding == 0 ? angleStepSize/100. : 0.0));

                if(!opts.flat_bkg){
                    //ctl
                    params[b]->cbkgctl0.init(1.0, -1.0, 1.0, 0.0);
                    params[b]->cbkgctl2.init(0.0, -1.0, 1.0, 0.01);
                    params[b]->cbkgctl4.init(0.0, -1.0, 1.0, 0.0);
                    if(opts.full_angular || opts.folding == 0){
                        params[b]->cbkgctl1.init(0.0, -1.0, 1.0, 0.01);
                        params[b]->cbkgctl3.init(0.0, -1.0, 1.0, 0.0);
                    }
                    else{
                        params[b]->cbkgctl1.init(0.0, -1.0, 1.0, 0.0);
                        params[b]->cbkgctl3.init(0.0, -1.0, 1.0, 0.0);
                    }

                    //ctk
                    params[b]->cbkgctk0.init(1.0, -1.0, 1.0, 0.0);
                    params[b]->cbkgctk2.init(0.05, -1.0, 1.0, 0.01);
                    params[b]->cbkgctk4.init(0.0, -1.0, 1.0, 0.0);
                    if(opts.full_angular || opts.folding != 4){
                        params[b]->cbkgctk1.init(0.1, -1.0, 1.0, 0.01);
                        params[b]->cbkgctk3.init(0.0, -1.0, 1.0, 0.0);
                    }
                    else{
                        params[b]->cbkgctk1.init(0.0, -1.0, 1.0, 0.0);
                        params[b]->cbkgctk3.init(0.0, -1.0, 1.0, 0.0);
                    }

                    //phi
                    params[b]->cbkgphi0.init(1.0, -1.0, 1.0, 0.0);
                    params[b]->cbkgphi2.init(0.0, -1.0, 1.0, 0.0);
                    params[b]->cbkgphi4.init(0.0, -1.0, 1.0, 0.0);
                    if(opts.full_angular || opts.folding == 0 || opts.folding == 3 || opts.folding == 4){
                        params[b]->cbkgphi1.init(0.0, -1.0, 1.0, 0.0);
                        params[b]->cbkgphi3.init(0.0, -1.0, 1.0, 0.0);
                    }
                    else{
                        params[b]->cbkgphi1.init(0.0, -1.0, 1.0, 0.0);
                        params[b]->cbkgphi3.init(0.0, -1.0, 1.0, 0.0);
                    }

                    params[b]->load_only_bckgnd_param_values(("fitresult_OnlyBckgnd_bin"+std::to_string(b)+".txt").c_str());

                }
            }
            else{
                //S-wave
                double sWaveStepSize = 0.01;
                params[b]->FS.init(SwaveFraction, 0.0, 1.0, sWaveStepSize);
                /*
        //JPSI Fit results
        params[b]->SS1.init( 0.549, -1.0, 1.0, (opts.full_angular || opts.folding != 4 ? sWaveStepSize : 0.0));
        params[b]->SS2.init(-0.128, -1.0, 1.0, (opts.full_angular || opts.folding == 1 ? sWaveStepSize : 0.0));
        params[b]->SS3.init( 0.003, -1.0, 1.0, (opts.full_angular || opts.folding == 2 ? sWaveStepSize : 0.0));
        params[b]->SS4.init(-0.006, -1.0, 1.0, (opts.full_angular || opts.folding > 2 ? sWaveStepSize : 0.0));
        params[b]->SS5.init(-0.149, -1.0, 1.0, (opts.full_angular ? sWaveStepSize : 0.0));
*/

                //B0 fit results
                params[b]->SS1.init(-0.231, -1.0, 1.0, (opts.full_angular || opts.folding != 4 ? sWaveStepSize : 0.0));
                params[b]->SS2.init( 0.023, -1.0, 1.0, (opts.full_angular || opts.folding == 1 ? sWaveStepSize : 0.0));
                params[b]->SS3.init( 0.003, -1.0, 1.0, (opts.full_angular || opts.folding == 2 ? sWaveStepSize : 0.0));
                params[b]->SS4.init( 0.001, -1.0, 1.0, (opts.full_angular || opts.folding > 2 ? sWaveStepSize : 0.0));
                params[b]->SS5.init(-0.068, -1.0, 1.0, (opts.full_angular ? sWaveStepSize : 0.0));
                params[b]->reset_parameters();

            }

            //create vectors to save parameter values for pull plots:
            if(n == 0){
                if(doToyPulls || doMCPulls){
                    UInt_t pp = 0;
                    for(UInt_t p = 0; p < params[b]->nparameters(); p++){
                        if(params[b]->get_parameter(p)->get_step_size() != 0.0){
                            if(b == 0){
                                std::vector< std::vector<double> >pulls_per_parameter;
                                pull_values.push_back(pulls_per_parameter);

                                std::vector< std::vector<double> >errors_per_parameter;
                                pull_errors.push_back(errors_per_parameter);

                                std::vector<double>starts_per_parameter;
                                pull_starts.push_back(starts_per_parameter);

                                var_indexs.push_back(p);
                            }
                            //check that vectors have same size and that no alloc error occurs:
                            assert(pp < pull_values.size());
                            assert(var_indexs.size() == pull_values.size());
                            assert(pp < pull_errors.size());
                            assert(var_indexs.size() == pull_errors.size());
                            assert(pp < pull_starts.size());
                            assert(var_indexs.size() == pull_starts.size());

                            //save vectors for values and start values to matrices
                            std::vector<double>pulls_per_param_per_bin;
                            pull_values.at(pp).push_back(pulls_per_param_per_bin);
                            std::vector<double>errors_per_param_per_bin;
                            pull_errors.at(pp).push_back(errors_per_param_per_bin);
                            pull_starts.at(pp).push_back(params[b]->get_parameter(p)->get_start_value());

                            //increase parameter index:
                            pp++;
                        }
                    }
                }
            }

            //pdf
            prob[b] = new fcnc::bu2kstarmumu_pdf(&opts, params[b]);
            if(opts.use_angular_acc || opts.weighted_fit){
                prob[b]->load_coeffs_eff_phsp_4d();
            }
            prob[b]->update_cached_normalization(params[b]);

            //create vector with events according to the requested fits/pulls
            selection.clear();

            Double_t nEvents = 0.;

            if(UseBinnedFit){
                nEvents *= f_signal.at(b)+f_bckgnd.at(b);
            }
            else
                nEvents = nToyEvents;

            std::cout << "Generate " << nEvents << " events!" << std::endl;
            gen = new fcnc::bu2kstarmumu_generator(&opts);
            selection = gen->generate(nEvents, params[b], prob[b]);
            //delete gen;
            if(!opts.full_angular && opts.always_generate_full_angular){
                for(UInt_t ee = 0; ee < selection.size(); ee++){
                    fldr.fold(&selection.at(ee));
                }
            }

            //check if all angles are almost 0.0
            if(false){
                for(UInt_t i = 0; i < selection.size(); i++)
                    if((fabs(selection.at(i).costhetak) < 0.001) && (fabs(selection.at(i).costhetal) < 0.001) && (fabs(selection.at(i).phi) < 0.001))
                        std::cout << "[WARNING]\tEvent" << i << "\tctk: " << selection.at(i).costhetak << "\tctl: " << selection.at(i).costhetal << "\tphi: " << selection.at(i).phi << std::endl;
            }

            if(n == 0)events_per_bin.push_back(selection.size());

            //deactive the S-wave for the fit!
            /*
                  opts.swave = false;
                  params[b]->FS.init(0.0, 0.0, 1.0, 0.0);
                  params[b]->SS1.init( 0.0, -1.0, 1.0, 0.0);
                  params[b]->SS2.init( 0.0, -1.0, 1.0, 0.0);
                  params[b]->SS3.init( 0.0, -1.0, 1.0, 0.0);
                  params[b]->SS4.init( 0.0, -1.0, 1.0, 0.0);
                  params[b]->SS5.init( 0.0, -1.0, 1.0, 0.0);
            */
            //set values to jpsi results for the fit!
            params[b]->FS.init(0.0, 0.0, 1.0, 0.0);
            params[b]->SS1.init( 0.549, -1.0, 1.0, 0.0);
            params[b]->SS2.init(-0.128, -1.0, 1.0, 0.0);
            params[b]->SS3.init( 0.003, -1.0, 1.0, 0.0);
            params[b]->SS4.init(-0.006, -1.0, 1.0, 0.0);
            params[b]->SS5.init(-0.149, -1.0, 1.0, 0.0);

            //make sure that the correct coefficients are chosen for the fit:
            if(!opts.full_angular && (opts.use_angular_acc || opts.weighted_fit) && opts.always_generate_full_angular){
                prob[b]->load_coeffs_eff_phsp_4d();
                prob[b]->update_cached_normalization(params[b]);
            }

            //fit the events:
            bool do_fit = true;
            int fitresult = 0;
            if(do_fit){
                fitresult = f.fit(prob[b], params[b], &selection);
            }
            else{ //don't fit, just update the efficiencies
                fitresult = 300;
                if(opts.use_angular_acc || opts.weighted_fit){
                    prob[b]->update_cached_efficiencies(params[b], &selection);
                }
            }

            fitresults.push_back(fitresult);
            //nametag per bin!
            std::string plotname;
            plotname = "ToyGen_"+std::to_string(nToyEvents)+"ToyEvents_FS_"+std::to_string(SwaveFraction)+"_";

            plotname.append(angularsuffix);
            if(UseBinnedFit){
                plotname.append("_bin"+std::to_string(b));
                if(Fit1bin)
                    plotname.append("_1BIN");
                if(Fit2bins)
                    plotname.append("_2BINS");
                if(FitAllbins)
                    plotname.append("_9BINS");
            }

            //plot the current pdf with event data:
            if(n == 0 && opts.write_eps){
                std::cout << "PLOT: " << selection.size() << " events" << std::endl;
                thePlotter.plot_data(prob[b], params[b], &selection, get_PullPlot_path(), plotname+"_Fit"+std::to_string(n), true);

                //add pdfs for all binning:
                thePDFs.at(b).push_back(prob[b]);
                theParameters.at(b).push_back(params[b]);

            }

            //save param values to pull_values vector:
            if(doToyPulls){
                for(UInt_t p = 0; p < var_indexs.size(); p++){
                    int index = var_indexs.at(p);

                    double val  = params[b]->get_parameter(index)->get_value();
                    double err  = params[b]->get_parameter(index)->get_error();

                    pull_values.at(p).at(b).push_back(val);
                    pull_errors.at(p).at(b).push_back(err);

                    //save values to tree:
                    treeValue = params[b]->get_parameter(index)->get_value();
                    treeStart = params[b]->get_parameter(index)->get_start_value();
                    treeError = params[b]->get_parameter(index)->get_error();
                    treeErrorUp = params[b]->get_parameter(index)->get_error_up();
                    treeErrorDown = params[b]->get_parameter(index)->get_error_down();

                    treeIndex = index;
                    treeBin = b;
                    treeFitNumber = n;
                    treeFitResult = fitresult;

                    //treeVarName = params[b]->get_parameter(index)->get_name();
                    //treeVarDesc = params[b]->get_parameter(index)->get_description();

                    T->Fill();

                }
            }
        }
    }

    //save TTree to .root file:
    if(doToyPulls){
        TFile * F = new TFile(("PullResults_"
                               +std::to_string(job_id == -10 ? nMCEvents : nToyEvents)
                               +(job_id == -10 ? "_MC_" : "_Toys_")
                               +std::to_string(nPulls)+"_Fits_"
                               +angularsuffix
                               +"_F_S_"+std::to_string(SwaveFraction)
                               +(job_id >= 0 ? "_job"+std::to_string(job_id) : "")
                               +".root").c_str(), "RECREATE");
        F->cd();
        T->Write();
        F->Close();
        delete T;
    }
}