/**
 * @file feldman_cousins.hh
 * @author Christoph Langenbruch, David Gerick, Renata Kopecna
 * @date 2021-01-12
 *
 */

#include <feldman_cousins.hh>

//TODO: removed what all is unused
#include <iostream>
#include <complex>
#include <assert.h>

#include <event.hh>
#include <parameters.hh>
#include <funcs.hh>
#include <pdf.hh>
#include <plotter.hh>
#include <options.hh>
#include <bu2kstarmumu_generator.hh>
#include <fitter.hh>
#include "folder.hh"

#include <thread>
#include <mutex>

#include <TMinuit.h>
#include <TFile.h>
#include <TTree.h>
#include <TRandom3.h>
#include <TMatrix.h>
#include <TH1D.h>
#include <TCanvas.h>
#include <TROOT.h>
#include <TStyle.h>
#include <TVector.h>
#include <TMatrixD.h>
#include <TMatrixDSym.h>
#include <TMatrixTSym.h>
#include <TMatrixT.h>
#include <TDecompChol.h>


#include <spdlog.h>

bool fcnc::feldman_cousins::fc_1d(std::string parname, unsigned int nsteps, double parmin, double parmax, unsigned int ntoys, unsigned int pdfidx,
                            std::vector<pdf*> probs, std::vector<parameters*> params, generator* gen, fitter* f,
                            std::vector<std::vector<event> *> events, unsigned int bin, int from, int to){
    std::ostringstream sout;
    bool TwoBins = false;
    if(opts->TheQ2binsmin.size() == 8)      TwoBins = false;
    else if(opts->TheQ2binsmin.size() == 2) TwoBins = true;
    else assert(opts->TheQ2binsmin.size() == 2 || opts->TheQ2binsmin.size() == 8);

    sout << "FCresult" << (TwoBins ? "_2BINS_" : "") << bin << "_" << parname << "_" << from << "_" << to << ".root";
    struct stat buffer;
    if(stat (sout.str().c_str(), &buffer) == 0){
        spdlog::warn("[SKIP]\t\tFC job already run as root-file '" + sout.str() + "' exists, continue to next job!");
        return false;
    }
    else{
        spdlog::info("Feldman-Cousins study for parameter="+ parname + " in q2bin={0:d}", bin);
    }
    //    if(!opts->full_angular)
    //      opts->always_generate_full_angular = false;
    bool cache_minos = opts->minos_errors;
    bool cache_hesse = opts->hesse_postrun;
    bool cache_sq_hesse = opts->squared_hesse;
    bool cache_shift = opts->shift_lh;
    //int cache_print = spdlog::default_logger_raw()->level();

    const unsigned int nPDF = probs.size();
    assert(params.size() == nPDF);
    assert(events.size() == nPDF);
    assert(pdfidx < nPDF);

    //temp save parameter set
    unsigned int nparameters = params.at(pdfidx)->nparameters();
    std::vector<std::vector<double> > smvalues (nPDF, std::vector<double> (nparameters, 0.0));
    for(unsigned int n = 0; n < nPDF; n++){
        for (unsigned int k=0; k<nparameters; k++){
            smvalues.at(n).at(k) = params.at(n)->get_parameter(k)->get_value();
        }
    }

    parameter* par = params.at(pdfidx)->get_parameter(parname);
    if (par != 0){
        opts->minos_errors = false;
        opts->squared_hesse = false;
        opts->hesse_postrun = false;
        opts->simplex_prerun = false;
        opts->shift_lh = false;

        opts->minuit_strategy = 2.0;
        //fitter f(opts);
        int statusdatafloated = f->fit(probs, params, events);
        double lhdatafloated = f->likelihood();
        double valuedatafloated = par->get_value();
        spdlog::info("Data fit finished with -2logLh = {0:f}",lhdatafloated );
        std::vector<std::vector<double> > floatingvalues(nPDF, std::vector<double> (nparameters, 0.0));
        std::vector<std::vector<double> > prev_error_low(nPDF, std::vector<double> (nparameters, 0.0));
        std::vector<std::vector<double> > prev_error_high(nPDF, std::vector<double> (nparameters, 0.0));
        std::vector<std::vector<double> > prev_meas(nPDF, std::vector<double> (nparameters, 0.0));
        std::vector<std::vector<bool> > isConstraint(nPDF, std::vector<bool> (nparameters, false));
        for(unsigned int n = 0; n < nPDF; n++){
            for (unsigned int k=0; k<nparameters; k++){
                floatingvalues.at(n).at(k)  = params.at(n)->get_parameter(k)->get_value();
                prev_error_low.at(n).at(k)  = params.at(n)->get_parameter(k)->get_previous_error_low();
                prev_error_high.at(n).at(k) = params.at(n)->get_parameter(k)->get_previous_error_high();
                prev_meas.at(n).at(k)       = params.at(n)->get_parameter(k)->get_previous_measurement();
                isConstraint.at(n).at(k)    = params.at(n)->get_parameter(k)->get_gaussian_constraint();
            }
        }

        for(unsigned int n = 0; n < nPDF; n++){
            for (unsigned int k=0; k<nparameters; k++){
                if(isConstraint.at(n).at(k)){
                    spdlog::info("[PDF{0:d}][PAR{1:d}]\tpar={2:s}: prev_measurement={3:f} prev_low={4:f}prev_high={5:f}",
                                 n, k, params.at(n)->get_parameter(k)->get_name(),
                                 prev_meas.at(n).at(k),prev_error_low.at(n).at(k) ,prev_error_high.at(n).at(k));
                }
            }
        }

        //opts->minuit_strategy = 1.0;//this should be 2 for the fixed data fit I think
        if (from < 0) from = 0;
        if (to > int(nsteps) || to < 0) to = nsteps;
        std::vector<double> lhdatafixed(to-from, 0.0);
        std::vector<int> statusdatafixed(to-from, 0);
        std::vector<std::vector<double> > lhtoysfixed;
        std::vector<std::vector<double> > lhtoysfloated;

        std::vector<std::vector<int> > statustoysfixed;
        std::vector<std::vector<int> > statustoysfloated;

        std::vector<std::vector< double> > params_dfloated (nPDF, std::vector< double> (nparameters, 0.0));
        std::vector<std::vector<std::vector<double>>>params_dfixed(nPDF,std::vector<std::vector<double>>
                                                                   (to-from, std::vector<double>(nparameters, 0.0)));
        std::vector<std::vector<std::vector< std::vector< double> > > > params_tfixed (nPDF,
                                                                                       std::vector<std::vector< std::vector< double> > > (to - from,
                                                                                                                                          std::vector< std::vector< double> > (ntoys,
                                                                                                                                                                               std::vector< double> (nparameters, 0.0))));
        std::vector<std::vector<std::vector< std::vector< double> > > > params_tfloated (nPDF,
                                                                                         std::vector<std::vector< std::vector< double> > > (to - from,
                                                                                                                                            std::vector< std::vector< double> > (ntoys,
                                                                                                                                                                                 std::vector< double> (nparameters, 0.0))));

        //buffer floated data values
        for(unsigned int n = 0; n < nPDF; n++){
            for (unsigned int k=0; k<nparameters; k++){
                params_dfloated.at(n).at(k) = params.at(n)->get_parameter(k)->get_value();
            }
        }

        std::string pardesc = par->get_description();
        for (int j=from; j<to; j++){//loop over steps
            double parcur = parmin+(j+0.5)*(parmax-parmin)/nsteps;
            spdlog::info("[STEP]\tEvaluating step {0:d}out of {1:d} steps!",j-from+1, to-from);
            spdlog::info("[INFO]\tStudying point "+parname+"={0:f}",parcur);

            //init parameters from floating fit values
            for(unsigned int n = 0; n < nPDF; n++){
                for(unsigned int k=0; k<nparameters; k++){
                    if(isConstraint.at(n).at(k)){
                        params.at(n)->get_parameter(k)->init(floatingvalues.at(n).at(k), params.at(n)->get_parameter(k)->get_min(), params.at(n)->get_parameter(k)->get_max(),
                                                             params.at(n)->get_parameter(k)->get_step_size(), prev_error_low.at(n).at(k), prev_error_high.at(n).at(k));
                        params.at(n)->get_parameter(k)->set_previous_measurement(prev_meas.at(n).at(k));
                    }
                    else{
                        params.at(n)->get_parameter(k)->init(floatingvalues.at(n).at(k), params.at(n)->get_parameter(k)->get_min(), params.at(n)->get_parameter(k)->get_max(),
                                                             params.at(n)->get_parameter(k)->get_step_size());
                    }
                }
                params.at(n)->reset_parameters();
            }
            for(unsigned int n = 0; n < nPDF; n++){
                params.at(n)->get_parameter(parname)->init(parcur, parmin, parmax, 0.0);
            }
            opts->minuit_strategy = 2.0;//more precise for fixed data fit (used for lh)
            opts->hesse_postrun = false;
            opts->simplex_prerun = false;
            spdlog::info("Start fit with fixed parameter");
            statusdatafixed.at(j-from) = f->fit(probs, params, events);
            lhdatafixed.at(j-from) = f->likelihood();
            spdlog::info("[DONE]\tFinished fit with fixed parameter");
            opts->minuit_strategy = 1.0;//less precise for toys -> speed more important

            opts->hesse_postrun = false;
            opts->simplex_prerun = false;

            std::vector<std::vector<double> > fitvalues (nPDF, std::vector<double> (nparameters, 0.0));
            for(unsigned int n = 0; n < nPDF; n++){
                //parameters from fixed point
                for (unsigned int k=0; k<nparameters; k++){
                    fitvalues.at(n).at(k) = params.at(n)->get_parameter(k)->get_value();
                }
                //buffer fixed data values
                for (unsigned int k=0; k<nparameters; k++){
                    params_dfixed.at(n).at(j-from).at(k) = params.at(n)->get_parameter(k)->get_value();
                }
            }

            //make all toys with these params!
            std::vector<std::vector<std::vector<event> > >toys(ntoys, std::vector<std::vector<event> > ());
            spdlog::info("Start generating toys");
            for (unsigned int k = 0; k < ntoys; k++){
                spdlog::debug("[TOYS]\tGenerating toy sample {0:d}/{1:d}",k+1,ntoys );
                for(unsigned int n = 0; n < nPDF; n++){
                    std::vector<event> toy = gen->generate(events.at(n)->size(), params.at(n), probs.at(n));
                    if(!opts->full_angular){
                        fcnc::folder fldr(opts);
                        for(UInt_t e = 0; e < toy.size(); e++){
                            fldr.fold(&toy.at(e));
                            if(toy.at(e).m < 1000.){
                                spdlog::warn("[CORRUPT]\tEvent={0:d} pdf={0:f} toy={0:f}",e, n, k );
                                spdlog::warn("[CORRUPT]\tctl={0:f} ctk={0:f} phi={0:f} q2={0:f}",
                                             toy.at(e).costhetal, toy.at(e).costhetak, toy.at(e).phi, toy.at(e).q2 );
                            }
                        }
                    }
                    toys.at(k).push_back(toy);
                }
            }
            spdlog::info("[DONE]\tGenerating toys");

            //init parameters from fixed point
            for(unsigned int n = 0; n < nPDF; n++){
                for(unsigned int k=0; k<nparameters; k++){
                    if(isConstraint.at(n).at(k)){
                        params.at(n)->get_parameter(k)->init(fitvalues.at(n).at(k), params.at(n)->get_parameter(k)->get_min(), params.at(n)->get_parameter(k)->get_max(),
                                                             params.at(n)->get_parameter(k)->get_step_size(), prev_error_low.at(n).at(k), prev_error_high.at(n).at(k));
                        params.at(n)->get_parameter(k)->set_previous_measurement(prev_meas.at(n).at(k));
                    }
                    else{
                        params.at(n)->get_parameter(k)->init(fitvalues.at(n).at(k), params.at(n)->get_parameter(k)->get_min(), params.at(n)->get_parameter(k)->get_max(),
                                                             params.at(n)->get_parameter(k)->get_step_size());
                    }
                }
            }
            //fit toys repeatedly
            std::vector<double> toysfixed;
            std::vector<double> toysfloated;
            std::vector<int> statusfixed;
            std::vector<int> statusfloated;
            spdlog::info("Start fitting all toys");
            for (unsigned int k=0; k<ntoys; k++){
                for(unsigned int n = 0; n < nPDF; n++){
                    params.at(n)->reset_parameters();//these are not the fitted/generated parameters->fix, they should float however, so should be ok?
                }
                for(unsigned int n = 0; n < nPDF; n++){
                    params.at(n)->get_parameter(parname)->init(parcur, parmin, parmax, 0.0);
                }

                //save pointers to toy events into vector:
                std::vector<std::vector<event> * > leToys;
                for(unsigned int n = 0; n < nPDF; n++){
                    leToys.push_back(&toys.at(k).at(n));
                }

                //fit with fixed parameter
                spdlog::info("[TOYS]\tFitting fixed toy sample {0:d}/{1:d} at step {2:d}/{3:d}",
                             k+1, ntoys, j-from+1, to-from );
                int resultfixed = f->fit(probs, params, leToys);
                toysfixed.push_back(f->likelihood());
                statusfixed.push_back(resultfixed);

                //buffer fixed toy values
                for(unsigned int n = 0; n < nPDF; n++){
                    for (unsigned int l=0; l<nparameters; l++){
                        params_tfixed.at(n).at(j-from).at(k).at(l) = params.at(n)->get_parameter(l)->get_value();
                    }
                    params.at(n)->reset_parameters();

                }
                for(unsigned int n = 0; n < nPDF; n++){
                    params.at(n)->get_parameter(parname)->init(parcur, parmin, parmax, 0.01);
                }

                //fit with floated parameter
                spdlog::info("[TOYS]\tFitting freed toy sample0:d}/{1:d} at step {2:d}/{3:d}",
                             k+1, ntoys, j-from+1, to-from );
                int resultfloated = f->fit(probs, params, leToys);
                toysfloated.push_back(f->likelihood());
                statusfloated.push_back(resultfloated);

                //buffer fixed toy values
                for(unsigned int n = 0; n < nPDF; n++){
                    for (unsigned int l=0; l<nparameters; l++){
                        params_tfloated.at(n).at(j-from).at(k).at(l) = params.at(n)->get_parameter(l)->get_value();
                    }
                }
            }
            lhtoysfixed.push_back(toysfixed);
            lhtoysfloated.push_back(toysfloated);
            statustoysfixed.push_back(statusfixed);
            statustoysfloated.push_back(statusfloated);
        }// end loop over steps
        spdlog::info("[DONE]\tFitting all toys");
        spdlog::debug("tLoop over all {0:d} steps completed. Writing results to file!",nsteps );

        //save to root file
        TFile* fout = new TFile(sout.str().c_str(), "RECREATE");
        fout->cd();
        TTree* t = new TTree(parname.c_str(), pardesc.c_str());
        int pdf=0, step=0, statustfixed=0, statustfloated=0, toy=0, statusdfixed=0, statusdfloated=0;
        double value=0.0, lhdfloated=lhdatafloated, lhdfixed=0.0, lhtfixed=0.0, lhtfloated=0.0;
        t->Branch("bin",&bin,"bin/I");
        t->Branch("pdf",&pdf,"pdf/I");
        t->Branch("step",&step,"step/I");
        t->Branch("value",&value,"value/D");
        t->Branch("toy",&toy,"toy/I");
        t->Branch("statusdatafixed",&statusdfixed,"statusdatafixed/I");
        t->Branch("statusdatafloated",&statusdfloated,"statusdatafloated/I");
        t->Branch("statustoyfixed",&statustfixed,"statustoyfixed/I");
        t->Branch("statustoyfloated",&statustfloated,"statustoyfloated/I");
        t->Branch("lhdatafloated",&lhdfloated,"lhdatafloated/D");
        t->Branch("valuedatafloated",&valuedatafloated,"valuedatafloated/D");
        t->Branch("lhdatafixed",&lhdfixed,"lhdatafixed/D");
        t->Branch("lhtoyfixed",&lhtfixed,"lhtoyfixed/D");
        t->Branch("lhtoyfloated",&lhtfloated,"lhtoyfloated/D");
        //also save the fitted values
        std::vector<double> tfloated(nparameters, 0.0), tfixed(nparameters, 0.0), dfloated(nparameters, 0.0), dfixed(nparameters, 0.0);
        for (unsigned int k = 0; k < nparameters; k++){
            fcnc::parameter* param = params.at(pdfidx)->get_parameter(k);
            if (param->get_step_size() != 0.0){
                std::string parname(param->get_name());
                std::string pardesc(param->get_description());
                //TTree* t = new TTree(parname.c_str(), pardesc.c_str());
                t->Branch((std::string("datafixed_")+parname).c_str(),&dfixed.at(k),(std::string("datafixed_")+parname+std::string("/D")).c_str());
                t->Branch((std::string("datafloated_")+parname).c_str(),&dfloated.at(k),(std::string("datafloated_")+parname+std::string("/D")).c_str());
                t->Branch((std::string("toyfixed_")+parname).c_str(),&tfixed.at(k),(std::string("toyfixed_")+parname+std::string("/D")).c_str());
                t->Branch((std::string("toyfloated_")+parname).c_str(),&tfloated.at(k),(std::string("toyfloated_")+parname+std::string("/D")).c_str());
            }
        }
        for(unsigned int n = 0; n < nPDF; n++){
            pdf = n; //TODO: what is this supposed to do?
            for (int j=from; j<to; j++){
                step = j;
                for (unsigned int k=0; k<ntoys; k++){
                    toy = k;
                    value          = parmin+(j+0.5)*(parmax-parmin)/nsteps;
                    statustfixed   = statustoysfixed.at(j-from).at(k);
                    statustfloated = statustoysfloated.at(j-from).at(k);
                    statusdfloated = statusdatafloated;
                    statusdfixed   = statusdatafixed.at(j-from);
                    lhtfixed       = lhtoysfixed.at(j-from).at(k);
                    lhtfloated     = lhtoysfloated.at(j-from).at(k);
                    lhdfixed       = lhdatafixed.at(j-from);
                    //lhdfloated = lhdatafloated; //always the same
                    for (unsigned int l = 0; l < nparameters; l++){
                        tfloated.at(l) = params_tfloated.at(n).at(j-from).at(k).at(l);
                        tfixed.at(l) = params_tfixed.at(n).at(j-from).at(k).at(l);
                        dfloated.at(l) = params_dfloated.at(n).at(l);
                        dfixed.at(l) = params_dfixed.at(n).at(j-from).at(l);
                    }
                    t->Fill();
                }
            }
        }
        t->Write();
        fout->Write();
        fout->Close();
        delete fout;
        /*
      //analyse results
      std::string suffix = "";
      TH1D* cl = new TH1D((std::string("cl")+parname+suffix).c_str(), (std::string(";")+pardesc+std::string(";Confidence level")).c_str(), nsteps, parmin, parmax);
      TH1D* lh = new TH1D((std::string("lh")+parname+suffix).c_str(), (std::string(";")+pardesc+std::string(";Confidence level")).c_str(), nsteps, parmin, parmax);
      for (unsigned int j=0; j<nsteps; j++)
      {
      unsigned int nlarger = 0;
      for (unsigned int k=0; k<ntoys; k++)
      {
      double dtoy = lhtoysfixed.at(j).at(k)-lhtoysfloated.at(j).at(k);
      double ddata = lhdatafixed.at(j)-lhdatafloated;
      if (dtoy > ddata)
      nlarger++;
      }
      lh->SetBinContent(j+1, TMath::Erf(sqrt(lhdatafixed.at(j)-lhdatafloated)/sqrt(2.0)));
      cl->SetBinContent(j+1, 1.0-double(nlarger)/ntoys);
      }
      double low = parmin;
      const double thecl = 0.683;
      if (cl->GetBinContent(1) > thecl)//???
      {
      for (int i=0; i<nsteps; i++)
      if (cl->GetBinContent(i+1) > thecl && cl->GetBinContent(i+2) < thecl)
      {
      low = cl->GetBinCenter(i+1) + (thecl-cl->GetBinContent(i+1))*(cl->GetBinCenter(i+2)-cl->GetBinCenter(i+1))/(cl->GetBinContent(i+2)-cl->GetBinContent(i+1));
      break;
      }
      }
      double high = parmax;
      if (cl->GetBinContent(nsteps) > thecl)
      {
      for (int i=nsteps-1; i>=0; i--)
      if (cl->GetBinContent(i+1) < thecl && cl->GetBinContent(i+2) > thecl)
      {
      high = cl->GetBinCenter(i+1) + (thecl-cl->GetBinContent(i+1))*(cl->GetBinCenter(i+2)-cl->GetBinCenter(i+1))/(cl->GetBinContent(i+2)-cl->GetBinContent(i+1));
      break;
      }
      }
      TCanvas* c1 = new TCanvas("c1", "c1", 1600, 1200);
      TLine* l = new TLine();
      c1->cd();
      cl->SetMinimum(0.0);
      cl->SetMaximum(1.1);
      cl->Draw();
      lh->SetLineColor(4);
      lh->Draw("lsame");
      l->SetLineColor(2);
      l->DrawLine(low, 0.0, low, 1.1);
      l->DrawLine(high, 0.0, high, 1.1);
      c1->Print((std::string("plots/cl")+parname+suffix+std::string(".eps")).c_str(),"eps");
      //delete cl;
      delete l;
      delete c1;
      delete cl;
      */
    }
    else{
        spdlog::critical("Could not find parameter "+parname);
        assert(0);
    }
    spdlog::info("Finished Feldman-Cousins study for parameter "+ parname);
    opts->hesse_postrun = cache_hesse;
    opts->squared_hesse = cache_sq_hesse;
    opts->minos_errors = cache_minos;
    opts->shift_lh = cache_shift;
    return true;
}

bool fcnc::feldman_cousins::fc_1d(std::string parname, unsigned int nsteps, double parmin, double parmax, unsigned int ntoys,
                            pdf* prob, parameters* params, generator* gen, fitter* f,
                            std::vector<event> * events, unsigned int bin, int from, int to){
    std::vector<pdf *> the_probs;
    the_probs.push_back(prob);

    std::vector<parameters *> the_params;
    the_params.push_back(params);

    std::vector<std::vector<event>* > the_events;
    the_events.push_back(events);

    return fc_1d(parname, nsteps, parmin, parmax, ntoys, 0,
                 the_probs, the_params, gen, f,
                 the_events, bin, from, to);
}

bool fcnc::feldman_cousins::fc_2d(std::string parxname, unsigned int nstepsx, double parxmin, double parxmax,
                            std::string paryname, unsigned int nstepsy, double parymin, double parymax,
                            unsigned int ntoys, unsigned int pdfidx, //this parameter just changes output file
                            std::vector<pdf*> probs, std::vector<parameters*> params, generator* gen, fitter* f,
                            std::vector<std::vector<event> *> events, unsigned int bin, int from, int to){
    spdlog::info("Feldman-Cousins study for parameters" + parxname +", " +paryname );
    /*
    bool cache_minos = opts->minos_errors;
    bool cache_hesse = opts->hesse_postrun;
    bool cache_shift = opts->shift_lh;
    int cache_print = opts->print_level;

    const unsigned int nPDF = probs.size();
    assert(params.size() == nPDF);
    assert(events.size() == nPDF);
    assert(pdfidx < nPDF);

    parameter* parx = params.at(pdfidx)->get_parameter(parxname);
    parameter* pary = params.at(pdfidx)->get_parameter(paryname);
    if (parx != 0 && pary != 0)
      {
      opts->minos_errors = false;
      opts->hesse_postrun = false;
      opts->shift_lh = false;

      //fitter f(opts);
      f->fit(probs, params, events);
      double lhdatafloated = f.likelihood();
      double valuexdatafloated = parx->get_value();
      double valueydatafloated = pary->get_value();
      spdlog::info("Data fit finished with -2logLh = {0:f}",lhdatafloated );

      opts->print_level = -1;
      //std::cout << from << " {0:f}",to << " {0:f}",from - to << std::endl;
      if (from < 0)
        from = 0;
      if (to > int(nstepsx*nstepsy) || to < 0)
        to = nstepsx*nstepsy;
      std::vector<double> lhdatafixed(to-from, 0.0);
      std::vector<std::vector<double> > lhtoysfixed;
      std::vector<std::vector<double> > lhtoysfloated;

      std::vector<std::vector<int> > statustoysfixed;
      std::vector<std::vector<int> > statustoysfloated;

      std::string parxdesc = parx->get_description();
      std::string parydesc = pary->get_description();
      for (int j=from; j<to; j++)
        {
          //xidx = j%nstepsx
          //yidx = j/nstepsx
          double parxcur = parxmin+((j%nstepsx)+0.5)*(parxmax-parxmin)/nstepsx;
          double parycur = parymin+((j/nstepsx)+0.5)*(parymax-parymin)/nstepsy;

          spdlog::info("Studying point" +parxname +"="+ parxcur + paryname + "=" +parycur );
          params.reset_parameters();

          parx->init(parxname, parxdesc, parxcur, parxmin, parxmax, 0.0);
          pary->init(paryname, parydesc, parycur, parymin, parymax, 0.0);
          f->fit(probs, params, events);
          lhdatafixed.at(j-from) = f->likelihood();

          //make all toys with these params!
          std::vector<std::vector<event> > toys;
          for (unsigned int k=0; k<ntoys; k++)
            {
        std::vector<event> toy = gen->generate(events.size(), params, prob);
        toys.push_back(toy);
            }
          //fit toys repeatedly
          std::vector<double> toysfixed;
          std::vector<double> toysfloated;
          std::vector<int> statusfixed;
          std::vector<int> statusfloated;
          for (unsigned int k=0; k<ntoys; k++)
            {
        params.reset_parameters();//these are not the fitted/generated parameters->fix, they should float however, so should be ok?
        parx->init(parxname, parxdesc, parxcur, parxmin, parxmax, 0.0);
        pary->init(paryname, parydesc, parycur, parymin, parymax, 0.0);
        int resultfixed = f->fit(probs, params, &toys.at(k));
        toysfixed.push_back(f->likelihood());
        statusfixed.push_back(resultfixed);

        params.reset_parameters();
        parx->init(parxname, parxdesc, parxcur, parxmin, parxmax, 0.01);
        pary->init(paryname, parydesc, parycur, parymin, parymax, 0.01);
        int resultfloated = f->fit(probs, params, &toys.at(k));
        toysfloated.push_back(f->likelihood());
        statusfloated.push_back(resultfloated);
            }
          lhtoysfixed.push_back(toysfixed);
          lhtoysfloated.push_back(toysfloated);
          statustoysfixed.push_back(statusfixed);
          statustoysfloated.push_back(statusfloated);
        }
      //std::cout <<"after loop" );

      //save to root file
      std::ostringstream sout;
      sout << "FCresult_" << bin << "_" << parxname << "_" << paryname << "_" << from << "_" << to << ".root";
      TFile* fout = new TFile(sout.str().c_str(), "RECREATE");
      fout->cd();
      TTree* t = new TTree((parxname+std::string("_")+paryname).c_str(), (std::string(";")+parxdesc+std::string(";")+parydesc).c_str());
      int step=0, stepx=0, stepy=0, statustfixed=0, statustfloated=0, toy=0;
      double valuex=0.0, valuey, lhdfloated=lhdatafloated, lhdfixed=0.0, lhtfixed=0.0, lhtfloated=0.0;
      t->Branch("bin",&bin,"bin/I");
      t->Branch("step",&step,"step/I");
      t->Branch("stepx",&stepx,"stepx/I");
      t->Branch("stepy",&stepy,"stepy/I");
      t->Branch("valuex",&valuex,"valuex/D");
      t->Branch("valuey",&valuey,"valuey/D");
      t->Branch("valuexdatafloated",&valuexdatafloated,"valuexdatafloated/D");
      t->Branch("valueydatafloated",&valueydatafloated,"valueydatafloated/D");
      t->Branch("toy",&toy,"toy/I");
      t->Branch("statustoyfixed",&statustfixed,"statustoyfixed/I");
      t->Branch("statustoyfloated",&statustfloated,"statustoyfloated/I");
      t->Branch("lhdatafloated",&lhdfloated,"lhdatafloated/D");
      t->Branch("lhdatafixed",&lhdfixed,"lhdatafixed/D");
      t->Branch("lhtoyfixed",&lhtfixed,"lhtoyfixed/D");
      t->Branch("lhtoyfloated",&lhtfloated,"lhtoyfloated/D");
      for (int j=from; j<to; j++)
        {
          for (unsigned int k=0; k<ntoys; k++)
            {
        step = j;
        stepx = j%nstepsx;
        stepy = j/nstepsx;
        toy = k;
        valuex = parxmin+((j%nstepsx)+0.5)*(parxmax-parxmin)/nstepsx;
        valuey = parymin+((j/nstepsx)+0.5)*(parymax-parymin)/nstepsy;
        statustfixed = statustoysfixed.at(j-from).at(k);
        statustfloated = statustoysfloated.at(j-from).at(k);
        lhtfixed = lhtoysfixed.at(j-from).at(k);
        lhtfloated = lhtoysfloated.at(j-from).at(k);
        lhdfixed = lhdatafixed.at(j-from);
        lhdfloated = lhdatafloated;
        t->Fill();
            }
        }
      t->Write();
      fout->Write();
      fout->Close();
      delete fout;
      }
    else{
      if (parx == 0) spdlog::error("Could not find parameter" +parxname );
      if (pary == 0) spdlog::error("Could not find parameter" +paryname );
        assert(0);
      }
    spdlog::info("Finished Feldman-Cousins study for parameters "+parxname + "," +paryname );
    opts->hesse_postrun = cache_hesse;
    opts->minos_errors = cache_minos;
    opts->shift_lh = cache_shift;
    opts->print_level = cache_print;
*/
    return true;
}

bool fcnc::feldman_cousins::fc_2d(std::string parxname, unsigned int nstepsx, double parxmin, double parxmax,
                            std::string paryname, unsigned int nstepsy, double parymin, double parymax,
                            unsigned int ntoys, //this parameter just changes output file
                            pdf* prob, parameters* params, generator* gen, fitter* f,
                            std::vector<event> * events, unsigned int bin, int from, int to){
    std::vector<pdf *> the_probs;
    the_probs.push_back(prob);

    std::vector<parameters *> the_params;
    the_params.push_back(params);

    std::vector<std::vector<event>* > the_events;
    the_events.push_back(events);

    return fc_2d(parxname, nstepsx, parxmin, parxmax,
                 paryname, nstepsy, parymin, parymax,
                 ntoys, 0,
                 the_probs, the_params, gen, f,
                 the_events, bin, from, to);
}