//Renata Kopecna

#include <angularcorr.hh>
#include <sstream>      // std::istringstream

#include <options.hh>
#include <event.hh>
#include <values.hh>
#include <helpers.hh>
#include <paths.hh>
#include <design.hh>

#include <bu2kstarmumu_loader.hh>
#include <bu2kstarmumu_parameters.hh>
#include <bu2kstarmumu_pdf.hh>

#include <spdlog.h>

#include <TFitResult.h>
#include <TH2D.h>
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
#include <TRandom3.h>
#include <TMatrixD.h>
#include <TStyle.h>


//-------------------------------
//     scan angular acceptance max. order of legendre
//-------------------------------


void logYears(std::vector<int>years, std::vector< std::vector<fcnc::event> > eves, int loadedSize){

    for(UInt_t y = 0; y < years.size(); y++){
        spdlog::info("Events in run {0:d} : {1:d}", years.at(y) , eves.at(y).size());
    }
    spdlog::info( "------------------------");
    spdlog::info("Events in total: {0:d}", loadedSize);
    return;
}

void logYears(int run, int evtsSize, int loadedSize){
    spdlog::info("Events in run {0:d} : {1:d}", run, evtsSize);
    spdlog::info("------------------------");
    spdlog::info("Events in total : {0:d}", loadedSize);
}

//-------------------------------
//     scan angular acceptance
//-------------------------------

std::vector<int> get_scan_low(bool quickTest, bool only_1D_chi2){
    //q2, thetak, thetal, phi
    if (quickTest) return {6,5,5,2};//use small range for quick tests on format or algorithms:
    if (only_1D_chi2) return {2,2,2,2};//use huge range for quick tests of 1D chi2:
    return {6,5,2,1};
}

std::vector<int> get_scan_high(bool quickTest, bool only_1D_chi2){
    //q2, thetak, thetal, phi
    if (quickTest) return {6,10,5,2};//use small range for quick tests on format or algorithms:
    if (only_1D_chi2) return {9,9,9,9};//use huge range for quick tests of 1D chi2:
    return {8,7,4,3};
}

std::vector<int> get_scan_range(std::vector<int> scan_low, std::vector<int> scan_high){
    std::vector<int> scan_range;
    for_indexed(auto item: scan_low){
        //If high >= low, raise an exception
        assert(scan_high.at(i) >= item);
        scan_range.push_back(scan_high.at(i) -item + 1);
    }
    return scan_range;
}

int sanityCheck(std::vector<fcnc::event>events, fcnc::options opts, bool PHSP, int bin, bool norm){
    //events: input events
    //opts: general options
    //PHSP: using PHSP or sigMC?
    //bin: # of bin in q2, when all, set to -1
    //norm: Normalize the plots?

    //Create parameter set
    fcnc::bu2kstarmumu_parameters * leParameters = new fcnc::bu2kstarmumu_parameters(&opts);
    //create PDF
    fcnc::bu2kstarmumu_pdf * lePDF = new fcnc::bu2kstarmumu_pdf(&opts, leParameters);
    std::vector<fcnc::event> *leEvents= new std::vector<fcnc::event>;

    opts.weighted_fit = true; //Jesus efin christ
    int nAngleBins = 20;

    for(auto meas: events){
        leEvents->push_back(meas);
    }
    leParameters->use_default();
    leParameters->take_current_as_start();

    lePDF->load_coeffs_eff_phsp_4d();
    opts.multiply_eff = true;
    opts.weighted_fit = true;
    lePDF->update_cached_efficiencies(leParameters, leEvents);
    lePDF->update_cached_normalization(leParameters);

    //check the distributions of weighted events
    TH1D * ctl = new TH1D ("h_ctl", "h_ctl", nAngleBins, CTL_MIN, CTL_MAX);
    TH1D * ctk = new TH1D ("h_ctk", "h_ctk", nAngleBins, CTK_MIN, CTK_MAX);
    TH1D * phi = new TH1D ("h_phi", "h_phi", nAngleBins, PHI_MIN, PHI_MAX);
    TH1D * q2  = new TH1D ("h_q2" , "h_q2" , nAngleBins, Q2_MIN_RANGE, Q2_MAX_RANGE);
    TH2D * ctkq2  = new TH2D ("h_ctkq2" , "h_ctkq2",  nAngleBins, CTK_MIN, CTK_MAX, nAngleBins, Q2_MIN_RANGE, Q2_MAX_RANGE);
    TH2D * ctlq2  = new TH2D ("h_ctlq2" , "h_ctlq2",  nAngleBins, CTL_MIN, CTL_MAX, nAngleBins, Q2_MIN_RANGE, Q2_MAX_RANGE);
    TH2D * phiq2  = new TH2D ("h_phiq2" , "h_phiq2",  nAngleBins, PHI_MIN, PHI_MAX, nAngleBins, Q2_MIN_RANGE, Q2_MAX_RANGE);
    TH1D * ctl_noW = new TH1D ("ctl_noW", "ctl_noW", nAngleBins, CTL_MIN, CTL_MAX);
    TH1D * ctk_noW = new TH1D ("ctk_noW", "ctk_noW", nAngleBins, CTK_MIN, CTK_MAX);
    TH1D * phi_noW = new TH1D ("phi_noW", "phi_noW", nAngleBins, PHI_MIN, PHI_MAX);
    TH1D * q2_noW  = new TH1D ("q2_noW" , "q2_noW" , nAngleBins, Q2_MIN_RANGE, Q2_MAX_RANGE);
    TH1D* hweight  = new TH1D("hweight",  ";w;",     100,        0.0, 1.0/EFF_CUTOFF);


    for(auto &meas: *leEvents){
        if(meas.weight < 0.0){
            spdlog::warn("negative weight event:\tm={0:f} \t\tctl={1:f} \tctk={2:f}\tphi={3:f} \tweight={4:f}", meas.m, meas.costhetal, meas.costhetak, meas.phi, meas.weight);
        }        
        ctl->Fill(meas.costhetal, meas.weight);
        ctk->Fill(meas.costhetak, meas.weight);
        phi->Fill(meas.phi,       meas.weight);
        q2 ->Fill(meas.q2,        meas.weight);
        ctkq2->Fill(meas.costhetak, meas.q2, meas.weight);
        ctlq2->Fill(meas.costhetal, meas.q2, meas.weight);
        phiq2->Fill(meas.phi, meas.q2, meas.weight);

        ctl_noW->Fill(meas.costhetal);
        ctk_noW->Fill(meas.costhetak);
        phi_noW->Fill(meas.phi);
        q2_noW ->Fill(meas.q2);
        if (meas.weight > 30){
            spdlog::warn("Event with very high weight {0:f} found!", meas.weight);
            print_event(meas);
        }

        hweight->Fill(meas.weight);

    }
    //phi->GetYaxis()->SetRangeUser(5000,7000);


    int tagNumber =1*opts.eff_order_q2+
                   100*opts.eff_order_costhetak+
                   10000* opts.eff_order_costhetal+
                   1000000*opts.eff_order_phi ;


    //Normalize to one

    if (norm){
        ctl->Scale(nAngleBins/ctl->Integral());
        ctk->Scale(nAngleBins/ctk->Integral());
        phi->Scale(nAngleBins/phi->Integral());
        q2->Scale(nAngleBins/q2->Integral());
        ctkq2->Scale(nAngleBins*nAngleBins/ctkq2->Integral());
        ctlq2->Scale(nAngleBins*nAngleBins/ctlq2->Integral());
        phiq2->Scale(nAngleBins*nAngleBins/phiq2->Integral());
    }

    std::string subfolder = get_angCorrPath(false);
    std::string tag = (PHSP ? "" : "MC_") + (bin==-1 ? "" : "_bin"+std::to_string(bin)) + tag_Run(opts.run) +"_" + std::to_string(tagNumber);
    if (PHSP){
        //Get the fitted lines
        TF1 *fitLine_ctl = fitStraightLine(ctl, "fitLine_ctl", CTL_MIN, CTL_MAX);
        TF1 *fitLine_ctk = fitStraightLine(ctk, "fitLine_ctk", CTK_MIN, CTK_MAX);
        TF1 *fitLine_phi = fitStraightLine(phi, "fitLine_phi", PHI_MIN, PHI_MAX);
        TF1 *fitLine_q2  = fitStraightLine(q2,  "fitLine_q2",  Q2_MIN_RANGE, Q2_MAX_RANGE);
        //and plot them
        plotAndSaveWithLine(ctl,fitLine_ctl,"c_ctl",subfolder+"debug_ctl_"+tag,"eps");
        plotAndSaveWithLine(ctk,fitLine_ctk,"c_ctk",subfolder+"debug_ctk_"+tag,"eps");
        plotAndSaveWithLine(phi,fitLine_phi,"c_phi",subfolder+"debug_phi_"+tag,"eps");
        plotAndSaveWithLine(q2, fitLine_q2, "c_q2" ,subfolder+"debug_q2_"+tag, "eps");
    }
    else{
        plotAndSave(ctl,"c_ctl",subfolder+"debug_ctl_"+tag,"eps");
        plotAndSave(ctk,"c_ctk",subfolder+"debug_ctk_"+tag,"eps");
        plotAndSave(phi,"c_phi",subfolder+"debug_phi_"+tag,"eps");
        plotAndSave(q2, "c_q2" ,subfolder+"debug_q2_"+tag, "eps");
    }
    plotAndSave(ctkq2, "c_ctkq2" ,subfolder+"debug_ctkq2_"+tag, "eps");
    plotAndSave(ctlq2, "c_ctlq2" ,subfolder+"debug_ctlq2_"+tag, "eps");
    plotAndSave(phiq2, "c_phiq2" ,subfolder+"debug_phiq2_"+tag, "eps");
    plotAndSave(hweight, "c_weight" ,subfolder+"debug_weights_"+tag, "eps");

    std::string outputPath = subfolder + "debug_angles"+ tag + ".root";
    TFile *output = new TFile(outputPath.c_str(),"RECREATE");
    output->cd();
    ctl->Write("",TObject::kWriteDelete);
    ctk->Write("",TObject::kWriteDelete);
    phi->Write("",TObject::kWriteDelete);
    q2->Write("",TObject::kWriteDelete);
    ctkq2->Write("",TObject::kWriteDelete);
    ctlq2->Write("",TObject::kWriteDelete);
    phiq2->Write("",TObject::kWriteDelete);
    hweight->Write("",TObject::kWriteDelete);
    output->Close();


    if (norm){ //save only when the option of normalization is on
        plotAndSave(ctl_noW,"c_ctl_noW",subfolder+"debug_ctl_noW_"+tag,"eps");
        plotAndSave(ctk_noW,"c_ctk_noW",subfolder+"debug_ctk_noW_"+tag,"eps");
        plotAndSave(phi_noW,"c_phi_noW",subfolder+"debug_phi_noW_"+tag,"eps");
        plotAndSave(q2_noW, "c_q2_noW" ,subfolder+"debug_q2_noW_"+tag, "eps");
    }

    delete ctl;
    delete ctk;
    delete phi;
    delete q2;
    delete ctlq2;
    delete ctkq2;
    delete phiq2;
    delete ctl_noW;
    delete ctk_noW;
    delete phi_noW;
    delete q2_noW;
    return 1;
}

int sanityCheck_MC(fcnc::options opts){

    spdlog::info("Checking MC distribution after correction");

    //load MC
    std::vector<fcnc::event> events;

    if (opts.run == 1){
        events = fcnc::load_events(get_theFCNCpath(1,1), "Events", -1);
    }
    else if (opts.run == 2){
        events = fcnc::load_events(get_theFCNCpath(1,2), "Events", -1);
    }
    else {
        std::vector<fcnc::event> events1 = fcnc::load_events(get_theFCNCpath(1,1), "Events", -1);
        std::vector<fcnc::event> events2 = fcnc::load_events(get_theFCNCpath(1,2), "Events", -1);
        events = {merge_two_evtVecs(events1,events2)};
    }
    for (unsigned int b = 0; b < opts.get_nQ2bins(); b++){
        std::vector<fcnc::event> filtered_events;
        for (auto &evt: events){

            if (evt.costhetal > opts.ctl_max) continue;
            if (evt.costhetal < opts.ctl_min) continue;
            if (evt.costhetak > opts.ctk_max) continue;
            if (evt.costhetak < opts.ctk_min) continue;
            if (evt.phi > opts.phi_max) continue;
            if (evt.phi < opts.phi_min) continue;
            if (evt.q2 > opts.TheQ2binsmax[b]) continue;
            if (evt.q2 < opts.TheQ2binsmin[b]) continue;
            filtered_events.push_back(evt);
        }
        sanityCheck(filtered_events, opts, false, b, false);
    }

    return 0;
}
std::vector< std::vector<fcnc::event> > select_event(fcnc::options opts,std::vector<int> years){
    //Load and check size of events in both runs

    std::vector<fcnc::event>events[2] = {
        fcnc::load_events(get_theFCNCpath(3,1), "Events", -1),
        fcnc::load_events(get_theFCNCpath(3,2), "Events", -1)
    };
    assert(events[0].size());
    assert(events[1].size());
    spdlog::debug("Load angular corrections per two years?\t{0:b}",opts.angacccorrpertwoyears);
    spdlog::debug("Load angular corrections for both runs?\t{0:b}",opts.angacccorrbothruns);

    std::vector<fcnc::event> loaded_events;
    if (opts.angacccorrbothruns)         loaded_events = merge_two_evtVecs(events[0], events[1]);
    else if (opts.angacccorrpertwoyears) loaded_events = events[1];
    else                                 loaded_events = events[opts.run-1];

    std::vector< std::vector<fcnc::event> > eves;

    if (opts.angacccorrperyear || opts.angacccorrpertwoyears){
        spdlog::debug("Seperate sample into {0:d} sub-samples. One per year!",years.size());
        //create indidividual vectors for every year
        std::vector<fcnc::event>eve_year(years.size());
        //fill events into the correct vector
        for(auto ev: loaded_events){
            for(uint y = 0; y < years.size(); y++){
                if(ev.year == years.at(y)) eves.at(y).push_back(ev);
            }
        }

        //print the corresponding numbers in every vector:
        logYears(years, eves, loaded_events.size());
    }
    else eves.push_back(loaded_events);

    //If the correction is done by two years only, merge the two vectors into one
    if (opts.angacccorrpertwoyears){
        eves = {merge_two_evtVecs(eves.at(0),eves.at(1))};
        spdlog::debug("Size of the final vector is {0:d}",eves.size());
        spdlog::debug("Size of the final vector at(0) is {0:d}",eves.at(0).size());
    }

    return eves;

}

int get_flat_Q2_weights(){

    spdlog::info("Start preparing weights for flat-flat PHSP.");
    std::string PHSPpath = getSelectedTuplePath(4,2017,2);

    TFile *output = new TFile(CONFIG_PHSP_WEIGHT.c_str(),"RECREATE");
    output->cd();
    spdlog::debug("Created " + std::string(output->GetPath()));
    TH1D *flatq2weights = new TH1D("flatq2weights","flatq2weights",100,Q2_MIN_RANGE,Q2_MAX_RANGE);
    TH1D *originalQ2 = new TH1D("originalQ2","originalQ2",100,Q2_MIN_RANGE,Q2_MAX_RANGE);


    spdlog::debug("Reading tree DecayTree from " + PHSPpath);
    TChain *tree = new TChain("DecayTree");
    tree->Add(PHSPpath.c_str());
    spdlog::debug("Events in tree DecayTree:\t{0:d}", tree->GetEntries());

    double Q2;
    tree->SetBranchStatus("Q2",1);
    tree->SetBranchAddress("Q2",&Q2);

    for (int e = 0; e<tree->GetEntries(); e++){
        //spdlog::trace("Q2 for event {0:d}:\t{1:f}",e,Q2);
        tree->GetEntry(e);
        originalQ2->Fill(Q2/1e6); //Fill q2 in GeV instead of MeV
    }


    for (int b =1; b <= originalQ2->GetNbinsX(); b++){
        flatq2weights->SetBinContent(b,1.0/originalQ2->GetBinContent(b)); //Create the weights
    }
    flatq2weights->Scale(flatq2weights->GetNbinsX()/ flatq2weights->Integral()); //Normalize the weights
    output->cd();
    originalQ2->Write("",TObject::kWriteDelete);
    flatq2weights->Write("",TObject::kWriteDelete);
    output->Close();
    spdlog::info("Done preparing weights for flat-flat PHSP.");

    return 0;
}

int scan_max_order_angular_ccorrection(fcnc::options opts,
                                       bool assumePhiEven,
                                       bool quickTest, bool test_4times1D,
                                       bool checkSignificance, bool runMinuit,
                                       bool checkFactorization, bool do3Dmoments){
    get_flat_Q2_weights();
    //only use chi2 from 1D projections
    opts.only_4_1D_chi2    = test_4times1D;
    //prevent any plot printing for the scan
    opts.write_eps         = quickTest || test_4times1D;
    //Set year/run options
    set_ang_year_options(opts);

    //define the range to scan through the parameters:
    //format: {q2, theta_K, theta_L, phi}
    std::vector<int> scan_low  = get_scan_low(quickTest, opts.only_4_1D_chi2);
    std::vector<int> scan_high = get_scan_high(quickTest, opts.only_4_1D_chi2);
    std::vector<int> scan_range = get_scan_range(scan_low, scan_high);

    spdlog::debug("Scan range:" + convert_vector_to_string(scan_range));

    //construct results in 2D histogram:
    const UInt_t nBinsX = scan_range.at(1)*scan_range.at(2); // cos(theta_K) and cos(theta_L)
    const UInt_t nBinsY = scan_range.at(0)*scan_range.at(3); // q^2 and phi

    //load phase-space MC events
    std::vector<int> years = get_years(opts.run,false, false);
    std::vector< std::vector<fcnc::event> > eves = select_event(opts, years);

    //object for exporting values from the parametrization process:
    fcnc::values valu;
    spdlog::debug("Eves vector size: {0:d}", eves.size());       

    for_indexed(auto evts: eves){

        if(opts.angacccorrperyear ){
            opts.year = years.at(i);
            spdlog::debug("[YEAR]\t\t{0:d}", opts.year);
        }
        else opts.year = -1;

        std::string sample_suffix;
        if (opts.angacccorrbothruns)   sample_suffix = "Run1and2";
        else if(opts.angacccorrperyear)sample_suffix = std::to_string(opts.year);
        else                           sample_suffix = "Run"+std::to_string(opts.run);

        TH2D * chi2histo = new TH2D(("ScanLegendreOrders_Chi2result_"+sample_suffix).c_str(),
                                    ("#chi^{2} scan of max. order of polynoms ("+sample_suffix+");#bf{#theta_{K}} - #theta_{L};#bf{q^{2}} - #phi").c_str(),
                                    nBinsX, -0.5, nBinsX - 0.5,
                                    nBinsY, -0.5, nBinsY - 0.5);

        TH2D * pvaluehisto = new TH2D(("ScanLegendreOrders_Pvalueresult_"+sample_suffix).c_str(),
                                      ("p-value scan of max. order of polynoms ("+sample_suffix+");#bf{#theta_{K}} - #theta_{L};#bf{q^{2}} - #phi").c_str(),
                                      nBinsX, -0.5, nBinsX - 0.5,
                                      nBinsY, -0.5, nBinsY - 0.5);


        //label the axese
        labelAngularScan(scan_low, scan_range, chi2histo);
        labelAngularScan(scan_low, scan_range, pvaluehisto);

        /*
        //individual histograms for each 1D projection:
        TH2D * dist_4_1D[4];
        std::string svar[4] = {"ctk", "ctl", "phi", "q2"};

        if(opts.only_4_1D_chi2){
            //was only needed to proof that all four 1D projections are indipendent.
            //keep the code for possible re-checks:
            for(int h = 0; h < 4; h++){
                dist_4_1D[h] = new TH2D(("ScanLegendreOrders_Chi2result_"+svar[h]+"_"+sample_suffix).c_str(),
                                        (svar[h]+": #chi^{2} scan of max. order of polynoms ("+sample_suffix+");#bf{#theta_{K}} - #theta_{L};#bf{q^{2}} - #phi").c_str(),
                                        nBinsX, -0.5, nBinsX - 0.5,
                                        nBinsY, -0.5, nBinsY - 0.5);
                //label the x-axis
                labelAngularScan(scan_low, scan_range, dist_4_1D[h]);
            }
        }
        */



        //determine max scan range:
        int range_low  = *min_element(std::begin(scan_low),std::end(scan_low));
        //for_indexed(auto item: scan_low) if(scan_low[c] < range_low)range_low = scan_low[c];
        int range_high = *max_element(std::begin(scan_high),std::end(scan_high));

        //get the histogram for all four varialbes:
        TH2D * chi2_4_1D = new TH2D(("Scan_Chi2result_"+sample_suffix).c_str(),
                                    ("#chi^{2} scan of max. order of polynoms ("+sample_suffix+");max. order;variable").c_str(),
                                    range_high - range_low + 1, range_low - 0.5, range_high + 0.5,
                                    4, -0.5, 3.5);
        //label the x-axis
        for(int b = range_low; b <= range_high; b++){
            chi2_4_1D->GetXaxis()->SetBinLabel(b - range_low + 1, (std::to_string(b)).c_str());
        }
        //label the y-axis
        chi2_4_1D->GetYaxis()->SetBinLabel(1, LATEX_Q2_CH());
        chi2_4_1D->GetYaxis()->SetBinLabel(2, LATEX_TK_CH());
        chi2_4_1D->GetYaxis()->SetBinLabel(3, LATEX_TL_CH());
        chi2_4_1D->GetYaxis()->SetBinLabel(4, LATEX_PHI_CH());

        if(opts.only_4_1D_chi2){  //Screams in python
            for(int i = 0; i < 4; i++){
                for(int j = scan_low.at(i); j <= scan_high.at(i); j++){

                    //set to max order to current scan value and all others to min value
                    if(i == 0)opts.eff_order_q2        = j;
                    else      opts.eff_order_q2        = scan_low.at(0);
                    if(i == 1)opts.eff_order_costhetak = j;
                    else      opts.eff_order_costhetak = scan_low.at(1);
                    if(i == 2)opts.eff_order_costhetal = j;
                    else      opts.eff_order_costhetal = scan_low.at(2);
                    if(i == 3)opts.eff_order_phi       = j;
                    else      opts.eff_order_phi       = scan_low.at(3);

                    spdlog::debug("[SCAN]\t\tq2 {0:d} \tcos(t_K) {1:d} \tcos(t_L)  {2:d} \tphi  {3:d}", opts.eff_order_q2, opts.eff_order_costhetak, opts.eff_order_costhetal , opts.eff_order_phi);

                    //All Unbiased parameters!
                    fcnc::bu2kstarmumu_parameters params(&opts);
                    //pdf
                    fcnc::bu2kstarmumu_pdf prob(&opts, &params);

                    //prob.parametrize_eff_phsp(selected);
                    int tagNumber = i*10+j;
                    spdlog::debug("Tag number:\t{0:d}", tagNumber);
                    prob.parametrize_eff_phsp_4d(evts, &valu,tagNumber, assumePhiEven, checkSignificance, runMinuit, checkFactorization, do3Dmoments);
                    //Tag for all possible combinations when testing

                    //save the one chi2 value:
                    double chi2 = 0.;
                    if(i == 0)     chi2 = valu.chi2_4_1D_q2;
                    else if(i == 1)chi2 = valu.chi2_4_1D_ctk;
                    else if(i == 2)chi2 = valu.chi2_4_1D_ctl;
                    else           chi2 = valu.chi2_4_1D_phi;
                    int ndof = 1;
                    if(i == 0)     ndof = valu.ndof_4_1D_q2;
                    else if(i == 1)ndof = valu.ndof_4_1D_ctl;
                    else if(i == 2)ndof = valu.ndof_4_1D_ctl;
                    else           ndof = valu.ndof_4_1D_phi;
                    if(chi2 / ndof < 2.) chi2_4_1D->SetBinContent(j - range_low + 1, i + 1, chi2 / ndof);
                }
            }
        }
        else{
            for(int qq = 0; qq < scan_range.at(0); qq++){
                for(int tk = 0; tk < scan_range.at(1); tk++){
                    for(int tl = 0; tl < scan_range.at(2); tl++){
                        for(int fi = 0; fi < scan_range.at(3); fi++){

                            //change max order of polynoms according to the scan:
                            opts.eff_order_q2        = qq+scan_low.at(0);
                            opts.eff_order_costhetak = tk+scan_low.at(1);
                            opts.eff_order_costhetal = tl+scan_low.at(2);
                            opts.eff_order_phi       = fi+scan_low.at(3);

                             int tagNumber = 1*opts.eff_order_q2+
                                            10*opts.eff_order_costhetak+
                                            100* opts.eff_order_costhetal+
                                            1000*opts.eff_order_phi ;
                                            //This works fine assuming all the orders are <10!

                            spdlog::info("[SCAN]\t\tq2 {0:d} \tcos(t_K) {1:d} \tcos(t_L)  {2:d} \tphi  {3:d}", opts.eff_order_q2, opts.eff_order_costhetak, opts.eff_order_costhetal , opts.eff_order_phi);
                            spdlog::debug("Tag number:\t{0:d}", tagNumber);

                            //All Unbiased parameters!
                            fcnc::bu2kstarmumu_parameters params(&opts);
                            //pdf
                            fcnc::bu2kstarmumu_pdf prob(&opts, &params);

                            //prob.parametrize_eff_phsp(selected);
                            prob.parametrize_eff_phsp_4d(evts, &valu, tagNumber, assumePhiEven, checkSignificance, runMinuit, checkFactorization, do3Dmoments);

                            double chi2 = valu.angacccorr_chi2;
                            int ndof    = valu.angacccorr_ndof;

                            chi2histo->SetBinContent(tk * scan_range.at(2) + tl + 1,
                                                     qq * scan_range.at(3) + fi + 1,  chi2 / ndof);

                            pvaluehisto->SetBinContent(tk * scan_range.at(2) + tl + 1,
                                                       qq * scan_range.at(3) + fi + 1,
                                                       TMath::Prob(chi2, ndof));
                            //P(chi2,ndof) represents the probability that the observed Chi-squared
                            //for a correct model should be less than the value chi2.
                            //The returned probability corresponds to 1-P(chi2,ndof),
                            //which denotes the probability that an observed Chi-squared exceeds
                            //the value chi2 by chance, even for a correct model.
                        }
                    }
                }
            }//loop over max. orders
        }

        //save the histogram with chi2/ndof output
        TCanvas * cChi2Ndof = new TCanvas(("cChi2Ndof_"+sample_suffix).c_str(), ("cChi2Ndof_"+sample_suffix).c_str());
        TCanvas * cPvalue   = new TCanvas(("cPvalue_"+sample_suffix).c_str(), ("cPvalue_"+sample_suffix).c_str());

        gStyle->SetPaintTextFormat("1.3f");
        std::string plotName = "plots/angular/test/ScanChi2MaxOrderPolynom" + sample_suffix+std::string(quickTest ? "_quickTest" : "");
        if(opts.only_4_1D_chi2){
            replace(plotName,"Polynom","Polynom_4_1D_");
            cChi2Ndof->cd();
            cChi2Ndof->SetLogz(false);
            chi2_4_1D->GetZaxis()->SetRangeUser(-1.0,2.5);
            chi2_4_1D->Draw("TEXTCOLZ");
            cChi2Ndof->SaveAs((plotName+".eps").c_str(), "eps");
        }
        else{
            cChi2Ndof->cd();
            chi2histo->GetZaxis()->SetRangeUser(0.98, 1.05);
            chi2histo->Draw("TEXTCOLZ");
            cChi2Ndof->SaveAs((plotName + ".eps").c_str(), "eps");
            cChi2Ndof->SetLogz();
            cChi2Ndof->SaveAs((plotName+"_log.eps").c_str(), "eps");

            replace(plotName,"Chi2Max","PvalueMax");
            cPvalue->cd();
            //pvaluehisto->GetZaxis()->SetRangeUser(1e-8, 2.5e-2);
            pvaluehisto->Draw("TEXTCOLZ");

            cPvalue->SaveAs((plotName+".eps").c_str(), "eps");
            cPvalue->SetLogz();
            cPvalue->SaveAs((plotName+"_log.eps").c_str(), "eps");


            replace(plotName,"PvalueMax","Max");
            TFile * f = new TFile((plotName+".root").c_str(), "UPDATE");
            f->cd();
            chi2histo  ->Write("", TObject::kWriteDelete);
            pvaluehisto->Write("", TObject::kWriteDelete);
            f->Close();
        }

        delete cChi2Ndof;
        delete cPvalue;

        delete chi2histo;
        delete pvaluehisto;

    }

    //end program after scan
    return 0;

}

//-------------------------------
//     fit angular acceptance
//-------------------------------
int get_angular_acceptance(fcnc::options opts, bool testSaving){

    //Mostly useless here, but in case it is needed quickly get the correction to PHSP from genLvl MC
    //Needed to correct for flat q2
    get_flat_Q2_weights();

    spdlog::info("[START] Determination of angular acceptance corrections");
    //Set year/run options
    set_ang_year_options(opts);

    //load phase-space MC
    std::vector<int> years = get_years(opts.run,false, false);
    std::vector< std::vector<fcnc::event> > eves = select_event(opts, years);

    //object for exporting values from the parametrization process:
    fcnc::values valu;

    for(UInt_t y = 0; y < (opts.angacccorrperyear ? years.size() : 1); y++){

        opts.year = opts.angacccorrperyear ? years.at(y) : -1;
        if(opts.angacccorrperyear) spdlog::debug("[YEAR]\t\t{0:d}",years.at(y));

        //All Unbiased parameters!
        fcnc::bu2kstarmumu_parameters params(&opts);
        //pdf
        fcnc::bu2kstarmumu_pdf prob(&opts, &params);

        //prob.parametrize_eff_phsp(selected);
        prob.parametrize_eff_phsp_4d(eves.at(y), &valu, 0, IS_PHI_EVEN, false, false, false, false);

        //save the coefficients to .dat file
        prob.save_coeffs_eff_phsp_4d();


        //Plot the corrected distributions
        sanityCheck(eves.at(y), opts, true, -1, true);

        //test succesfull saving by loading and checking the coefficients:
        if(testSaving){
            //load the coefficients from file:
            std::vector<double>testCoeffs;
            testCoeffs = prob.read_coeffs_eff_phsp_4d();
            //check if the correct number of coeffs is saved and reloaded
            assert(testCoeffs.size() == prob.coeffs_eff_4d.size());

            //test if all coeffs are identical
            for(UInt_t c = 0; c < prob.coeffs_eff_4d.size(); c++){
                if(TMath::Abs(testCoeffs.at(c) - prob.coeffs_eff_4d.at(c)) > 1.0e-14){
                    spdlog::warn(" Coeff #{0:d} unequal.",c+1);
                    spdlog::warn(" Saved: {0:f}", prob.coeffs_eff_4d.at(c));
                    spdlog::warn(" Loaded: {0:f}", testCoeffs.at(c));
                    spdlog::warn(" Difference: {0:f}", testCoeffs.at(c) - prob.coeffs_eff_4d.at(c));
                }
            }
        }


        npolynom npol(&opts);
        assert(prob.coeffs_eff_4d.size() == npol.getSize());

        //this can be varying parameter, or can be fixed for example to the middle of the bin
        //double q2 = params->eff_q2();
        TRandom3* rnd = new TRandom3(432759);
        std::vector<fcnc::event> negativeeff;
        const bool testdata = false;
        const bool testtoy = false;
        unsigned int ntest = 100000;
        if (testdata){
            std::vector<fcnc::event> dataevents = fcnc::filterResonances(fcnc::load_events( get_theFCNCpath(0, opts.run), "Events", ntest));
            ntest = 0;//dataevents.size();
            for (unsigned int l=0; l<dataevents.size(); l++){
                fcnc::event meas = dataevents.at(l);
                if ( (meas.q2 < Q2_MIN_RANGE)
                     || (meas.q2 > Q2_MAX_RANGE)
                     || (meas.m < opts.m_low || meas.m > opts.m_high) )
                    continue;

                continue;
                ntest++;

                double eff = 0.0;
                unsigned int bin = 0;
                double costhetal = meas.costhetal;
                double costhetak = meas.costhetak;
                double phi = meas.phi;
                double q2 = meas.q2;
                for (unsigned int h = 0; h < npol.q2; h++){
                    for (unsigned int i = 0; i < npol.ctl; i++){
                        for (unsigned int j = 0; j < npol.ctk; j++){
                            for (unsigned int k = 0; k < npol.phi; k++){
                                bin = npol.get_bin_in4D(h,i,j,k);
                                eff += prob.coeffs_eff_4d.at(bin)*pow(q2, h)*pow(costhetal, i)*pow(costhetak, j)*pow(phi, k);
                            }
                        }
                    }
                }
                if (eff < 0.0){
                    meas.weight = eff;
                    negativeeff.push_back(meas);
                }
            }
        }
        else if (testtoy){
            for (unsigned int l=0; l<ntest; l++){
                if (l % (ntest/100) == 0) spdlog::debug("{0:d}%", l/(ntest/100));
                double eff = 0.0;
                unsigned int bin = 0;
                double costhetal = CTL_MIN + rnd->Rndm()*CTL_RANGE();
                double costhetak = CTK_MIN + rnd->Rndm()*CTK_RANGE();
                double phi       = PHI_MIN + rnd->Rndm()*PHI_RANGE();
                double q2        = Q2_MIN_RANGE + rnd->Rndm()*Q2_RANGE_FULL();
                for (unsigned int h = 0; h < npol.q2; h++){
                    for (unsigned int i = 0; i < npol.ctl; i++){
                        for (unsigned int j = 0; j < npol.ctk; j++){
                            for (unsigned int k = 0; k < npol.phi; k++){
                                bin = npol.get_bin_in4D(h,i,j,k);
                                eff += prob.coeffs_eff_4d.at(bin)*pow(q2, h)*pow(costhetal, i)*pow(costhetak, j)*pow(phi, k);
                            }
                        }
                    }
                }
                if (eff < 0.0){
                    fcnc::event meas;
                    meas.costhetal = costhetal;
                    meas.costhetak = costhetak;
                    meas.phi = phi;
                    meas.q2 = q2;
                    meas.weight = eff;
                    negativeeff.push_back(meas);
                }
                else if (eff > 3.0) spdlog::debug("LARGE {0:f} {1:f} {2:f} {3:f} {4:f}", costhetal, costhetak, phi, q2, eff);
            }
        }
        if(testtoy || testdata){
            spdlog::info( "looping over {0:d}  events"+ std::string(testtoy ? " (TOYS)" : " (DATA)"),ntest);
            for (unsigned int i=0; i<negativeeff.size();i++){
                const fcnc::event& meas = negativeeff.at(i);
                spdlog::debug(  "{0:f} {1:f} {2:f} {3:f} {4:f} {5:f}", meas.costhetal, meas.costhetak, meas.phi, meas.q2, meas.m, meas.weight);
            }
            spdlog::debug( "negative fraction {0:f}", double(negativeeff.size())/ntest);
        }

        if(opts.only_4_1D_chi2){
            spdlog::debug( "Results from >>>> " + std::string(opts.year == -1 ? (opts.angacccorrbothruns ? "Run1+2" : "Run"+std::to_string(opts.run)) : std::to_string(opts.year)) + " <<<< data set:");
            spdlog::debug( "Reduced chi2 of 4*1D projections:\t {0:f} ({1:f}/{2:f})",
                           valu.chi2_4_1D/ valu.ndof_4_1D, valu.chi2_4_1D,valu.ndof_4_1D);
            spdlog::debug( "Reduced chi2 of 4*1D projections: ctk \t {0:f} ({1:f}/{2:f})",
                           valu.chi2_4_1D_ctk/valu.ndof_4_1D_ctk, valu.chi2_4_1D_ctk, valu.ndof_4_1D_ctk);
            spdlog::debug( "Reduced chi2 of 4*1D projections: ctl\t {0:f} ({1:f}/{2:f})",
                           valu.chi2_4_1D_ctl/valu.ndof_4_1D_ctl, valu.chi2_4_1D_ctl, valu.ndof_4_1D_ctl);
            spdlog::debug( "Reduced chi2 of 4*1D projections: phi\t {0:f} ({1:f}/{2:f})",
                           valu.chi2_4_1D_phi/valu.ndof_4_1D_phi, valu.chi2_4_1D_phi, valu.ndof_4_1D_phi);
            spdlog::debug( "Reduced chi2 of 4*1D projections: q2\t {0:f} ({1:f}/{2:f})",
                           valu.chi2_4_1D_q2/valu.ndof_4_1D_q2, valu.chi2_4_1D_q2, valu.ndof_4_1D_q2 );
        }
    }

    return 0;
}


//-------------------------------
//     fit angular resolution
//-------------------------------
int get_angular_resolution(fcnc::options opts, std::vector<int> years){
    //Don't forget to set the full ctk range when doing this!

    //create resolution histograms for all three angles (x - x_true)
    const unsigned int leBins = 50;
    const double leRange = 0.12;


    spdlog::info("Determine angular resolutions from signal MC true information!");
    gROOT->SetBatch(kTRUE);
    gROOT->SetStyle("Plain");
    set_gStyle();

    assert(years.size());
    opts.reject_identical_muon_TRUEID = true;
    fcnc::bu2kstarmumu_loader loader(&opts);

    double resolution[3][years.size()];

    for_indexed(auto yr: years){
        spdlog::info("Evaluate year {0:d}", yr);
        //load tuples:
        std::vector<fcnc::event> MCevents = loader.read_full_tuple(yr,
                                                                   getSelectedTuplePath(1, get_yearFromRun(yr),
                                                                                        yr),
                                                                   get_inputTree_name(1),
                                                                   true, false, false, -1);

        std::vector<fcnc::event> MCtrue   = loader.read_full_tuple(yr,
                                                                   getSelectedTuplePath(1, get_yearFromRun(yr),
                                                                                        yr),
                                                                   get_inputTree_name(1),
                                                                   true, true, false, -1);

        assert(MCevents.size() == MCtrue.size());

        //create resolution histograms for all three angles (x - x_true)

        TCanvas * cReso = nullptr;
        std::string sangle[3] = {"ctk", "ctl", "phi"};
        std::string latexangle[3] = {LATEX_CTK,LATEX_CTL,LATEX_PHI};

        TH1D * h[3];
        for(unsigned int i = 0; i < 3; i++){
            h[i] = new TH1D(("h"+sangle[i]+"_"+std::to_string(yr)).c_str(),
                            (";"+latexangle[i]+" - "+latexangle[i]+"_{true};Events").c_str(),
                            leBins, -leRange, +leRange);
        }

        //fill histograms:
        for(unsigned int e = 0; e < MCevents.size(); e++){
            fcnc::event measMC   = MCevents.at(e);
            fcnc::event measTRUE = MCtrue.at(e);

            h[0]->Fill(measMC.costhetak - measTRUE.costhetak);
            h[1]->Fill(measMC.costhetal - measTRUE.costhetal);
            h[2]->Fill(measMC.phi - measTRUE.phi);
        }

        spdlog::info("Entries in histograms: {0:f}", h[0]->Integral());

        TF1 * fGauss = new TF1("theGaussian", "[0]*([2]*exp(-0.5*((x-[1])/[3])**2)+(1.-[2])*exp(-0.5*((x-[1])/[4])**2))", -leRange, +leRange);

        fGauss->SetLineStyle(kDashed);
        fGauss->SetLineColor(kMagenta - 3);

        //Set init parameters differently for each variable
        double initSigma[3]      = {0.0125, 0.005, 0.025};
        double sigmaLowRange[3]  = {0.008,  0.001, 0.01};
        double sigmaHighRange[3] = {0.020,  0.01,  0.04};
        for(int j = 0; j < 3; j++){
            fGauss->SetParameter(0, MCevents.size());
            fGauss->SetParameter(1, 0.);
            fGauss->SetParameter(2, 0.5);
            fGauss->SetParLimits(2, 0.0, 1.0);
            //Set the init resolution extra for ctl, as it is very narrow compared to phi and ctk

            fGauss->SetParameter(3, initSigma[j]);
            fGauss->SetParLimits(3, sigmaLowRange[j], sigmaHighRange[j]);
            fGauss->SetParameter(4, initSigma[j]);
            fGauss->SetParLimits(4, sigmaLowRange[j], sigmaHighRange[j]);
            TFitResultPtr fr;
            do{
                h[j]->GetListOfFunctions()->Clear();
                fr = h[j]->Fit(fGauss, "MS+Q");
                spdlog::info("Fit status = {0:d}", fr->Status());
            } while (!fr->IsValid());
            TMatrixD covmat = fr->GetCovarianceMatrix();
            if (spdlog::default_logger_raw()->level()==spdlog::level::trace) covmat.Print();

            double rel_frac  = fGauss->GetParameter(2);
            double drel_frac = fGauss->GetParError(2);
            double sigma1  = fGauss->GetParameter(3);
            double dsigma1 = fGauss->GetParError(3);
            double sigma2  = fGauss->GetParameter(4);
            double dsigma2 = fGauss->GetParError(4);

            double eff_sigma  = sqrt(rel_frac*sigma1*sigma1 + (1.-rel_frac)*sigma2*sigma2);
            double deff_sigma = 1./(2.*eff_sigma) * sqrt(
                        TMath::Power(2*rel_frac*sigma1*dsigma1, 2)
                        + TMath::Power(2*(1.-rel_frac)*sigma2*dsigma2, 2)
                        + TMath::Power((sigma1*sigma1-sigma2*sigma2)*drel_frac, 2)
                        + 2 * 2*rel_frac*sigma1 * 2*(1.-rel_frac)*sigma2 * covmat[3][4]
                    + 2 * 2*rel_frac*sigma1 * (sigma1*sigma1-sigma2*sigma2) * covmat[3][2]
                    + 2 * 2*(1.-rel_frac)*sigma2 * (sigma1*sigma1-sigma2*sigma2) * covmat[4][2]);

            cReso = new TCanvas((sangle[j]+"_reso_"+std::to_string(yr)).c_str(),
                                (sangle[j]+"_reso_"+std::to_string(yr)).c_str(), 1600, 1200);
            TLatex * leg = getPrettyTex(0.06, 13);

            std::ostringstream sRMS;
            sRMS << std::fixed << std::setprecision(4) << "#sigma = " << eff_sigma << " #pm " << deff_sigma;
            resolution[j][i] = eff_sigma;

            cReso->SetMargin(0.125,0.05,0.115,0.05);
            cReso->cd();
            h[j]->GetYaxis()->SetTitleOffset(1.15);
            h[j]->GetYaxis()->SetTitleSize(0.05);
            h[j]->GetYaxis()->SetLabelSize(0.05);
            h[j]->GetXaxis()->SetTitleSize(0.05);
            h[j]->GetXaxis()->SetLabelSize(0.05);
            h[j]->Draw("E1");
            leg->DrawLatex(0.15,0.84, sRMS.str().c_str());
            leg->DrawLatex(0.7,0.88, "LHCb MC");
            leg->DrawLatex(0.7,0.82, std::to_string(yr).c_str());

            cReso->Print(get_angResoPlot_path(sangle[j],yr,"eps").c_str(), "eps");
            cReso->SaveAs(get_angResoPlot_path(sangle[j],yr,"root").c_str());

            delete h[j];
        }
    }

    //print out resolutions:
    spdlog::info("resolution: {");
    for(int i = 0; i < 3; i++){
        std::cout << "{";
        for(UInt_t y = 0; y < years.size(); y++){
            std::cout << std::fixed << std::setprecision(6) << resolution[i][y];
            if(y < years.size() - 1)
                std::cout << ", ";
        }
        std::cout << "}" << std::endl;
    }
    return 0;
}