#define ULMAPS_TEST_CXX_


#include "ulmaps_test.h"
// C++ libraries
#include <iostream>
#include <fstream>

// Other libraries
#include <TCanvas.h>
#include <TColor.h>
#include <TComplex.h>
#include <TFitResult.h>
#include <TF1.h>
#include <TGraph.h>
#include <TH1.h>
#include <TH2.h>
#include <TLatex.h>
#include <TLegend.h>
#include <TLine.h>
#include <TMath.h>
#include <TMinuit.h>
#include <TROOT.h>
#include <TStyle.h>

#include <Math/DistFuncMathCore.h>

void ulmaps_test(Int_t i_mass, Int_t i_width)
{
  TH1::AddDirectory(kFALSE);
  output_file = new TFile("test.root", "RECREATE");

  ULcoh(i_mass, i_width);

  output_file->Write();
  delete output_file;

  printf("root ulmaps.C+(%d, %d) done\n", i_mass, i_width);
}

void ULcoh(Int_t i_mass, Int_t i_width)
{
  const int mass_n = 380;
  double mass_min = 4.010;
  double mass_max = 4.390;

  const int width_n = 270;
  double width_min = 0.030;
  double width_max = 0.300;

  long long fits_n = mass_n * width_n / 12.;
  long long fits_i = 0;
  double start_time = time(0);

  TH2D * h2d_scan_coh = NewTH2D("h2d_scan_coh",
                                 mass_n,  mass_min,  mass_max,
                                width_n, width_min, width_max,
                                "mass m / GeV", "width #Gamma / GeV", "#sigma_{coh}^{UL} / pb");

  const int n_params = 6;
  double best_minfcn = 0.;

  double best_c_par = 0.;
  double best_c_err = 0.;
  double best_lambda_par = 0.;
  double best_lambda_err = 0.;

  double lastPhi=0;
  double lastC=0;
  double lastLambda=0;

  double par[n_params];
  double err[n_params];
  double min[n_params];
  double max[n_params];
// loop over mass and width range for possible resonance
  for (int i = 0; i < mass_n; i++)
  {
    //for time reasons this is divided into block. The map has 12 blocks. I want to start 12 processes simultaneously. For test only 0,0 as input for ulmaps.C
    double mass = h2d_scan_coh->GetXaxis()->GetBinCenter(i+1);
    if ((mass < (4.000 + 0.100 * i_mass)) ||
        (mass > (4.100 + 0.100 * i_mass))) continue;

    for (int j = 0; j < width_n; j++)
    {
      double width = h2d_scan_coh->GetYaxis()->GetBinCenter(j+1);
      if ((width < (0.000 + 0.100 * i_width)) ||
          (width > (0.100 + 0.100 * i_width))) continue;

      //do a quick scan for "phi" parameter so have better start parameters for the real scan!
      {
        TMinuit * minuit = new TMinuit(n_params);
        minuit->SetFCN(MinuitFCN_ContinuumBWcoh);

        minuit->SetPrintLevel(-1);
        minuit->SetMaxIterations(5000);

        par[0] = 6.08484e+00;
        par[1] = 5.92729e+00;
        par[2] = mass;
        par[3] = width;
        par[4] = 3.11338e-02;
        par[5] = 3.41615e+00;

        err[0] = 0.01;
        err[1] = 0.01;
        err[2] = 0.01;
        err[3] = 0.01;
        err[4] = 0.01;
        err[5] = 0.01;

        min[0] = 1.;
        min[1] = 1.;
        min[2] = 0.;
        min[3] = 0.;
        min[4] = 0.;
        min[5] = 0. * TMath::Pi();

        max[0] = 10.;
        max[1] = 10.;
        max[2] = 5.;
        max[3] = 0.5;
        max[4] = 100.;
        max[5] = 2. * TMath::Pi();

        int ierflg = 0;
        double arglist[10] = {0.};

        // Define start parameters
        minuit->mnparm(0, "C",      par[0], err[0], min[0], max[0], ierflg);
        minuit->mnparm(1, "lambda", par[1], err[1], min[1], max[1], ierflg);
        minuit->mnparm(2, "mass",   par[2], err[2], min[2], max[2], ierflg);
        minuit->mnparm(3, "width",  par[3], err[3], min[3], max[3], ierflg);
        minuit->mnparm(4, "A",      par[4], err[4], min[4], max[4], ierflg);
        minuit->mnparm(5, "phase",  par[5], err[5], min[5], max[5], ierflg);

        minuit->mnexcm("SET NOW", arglist, 1, ierflg);
        arglist[0] = 1;
        minuit->mnexcm("SET ERR", arglist, 1, ierflg);
        arglist[0] = 2;
        minuit->mnexcm("SET STR", arglist, 1, ierflg);
        minuit->mnexcm("CALL FCN", arglist, 1, ierflg);

        // Fix some parameters
        minuit->FixParameter(2);
        minuit->FixParameter(3);

        // Scan all free parameters to set a new phi window
        minuit->mncomd("scan 0 1000", ierflg);
        minuit->GetParameter(5, par[5], err[5]);
        min[5] = par[5] - TMath::Pi();
        max[5] = par[5] + TMath::Pi();
        minuit->mnparm(5, "phase", par[5], err[5], min[5], max[5], ierflg);
        
        // Perform maximum loglikelihood fit
        minuit->Migrad();
        
        int npari, nparx, istat;
        double minfcn, edm, errdef;
        minuit->mnstat(minfcn, edm, errdef, npari, nparx, istat);

        best_minfcn = minfcn;

        //get bkg parameters (c and lambda) for maximum LL in Amplitude scan
        minuit->GetParameter(0, best_c_par, best_c_err);
        minuit->GetParameter(1, best_lambda_par, best_lambda_err);
        minuit->GetParameter(5, lastPhi, err[5]);
        lastC=best_c_par;
        lastLambda=best_lambda_par;
    
        delete minuit;
      }
     //assumption: we scan whole Likelihood when we set scan_max to 10 (this is true, but i had to test that first)
      const int scan_n = 1; //scan_n is how many scan points I have for the likelihood scan; set it to 100 because we want to have good resolution
      double scan_min = 0.;
      double scan_max = 10.;
      PrintProcessStatus(start_time, time(0), fits_i++, fits_n, 50, 10, true);
      TH1D * h_scan = NewTH1D("h_scan_coh", scan_n, scan_min, scan_max, "");
      double A_i = 0;
      int npari, nparx, istat;
      double minfcn, edm, errdef;
      int ierflg = 0;
      double arglist[10] = {0.};
      double L_i = 0;
      //scan starts here
      for (int k = 0; k < scan_n; k++)
      {
        A_i = h_scan->GetXaxis()->GetBinCenter(k+1);

        TMinuit * minuit = new TMinuit(n_params);
        minuit->SetFCN(MinuitFCN_ContinuumBWcoh);

        // set print level (-1: no output, 0: standard output)
        minuit->SetPrintLevel(-1);
        minuit->SetMaxIterations(5000);

        par[0] = lastC;
        par[1] = lastLambda;
        par[2] = mass;
        par[3] = width;
        par[4] = A_i;
        par[5] = lastPhi;

        err[0] = 0.01;
        err[1] = 0.01;
        err[2] = 0.01;
        err[3] = 0.01;
        err[4] = 0.01;
        err[5] = 0.01;

        min[0] = 1.;
        min[1] = 1.;
        min[2] = 0.;
        min[3] = 0.;
        min[4] = 0.;
        min[5] = 0. * TMath::Pi();

        max[0] = 10.;
        max[1] = 10.;
        max[2] = 5.;
        max[3] = 0.5;
        max[4] = 100.;
        max[5] = 2. * TMath::Pi();

        // Define start parameters
        minuit->mnparm(0, "C",      par[0], err[0], min[0], max[0], ierflg);
        minuit->mnparm(1, "lambda", par[1], err[1], min[1], max[1], ierflg);
        minuit->mnparm(2, "mass",   par[2], err[2], min[2], max[2], ierflg);
        minuit->mnparm(3, "width",  par[3], err[3], min[3], max[3], ierflg);
        minuit->mnparm(4, "A",      par[4], err[4], min[4], max[4], ierflg);
        minuit->mnparm(5, "phase",  par[5], err[5], min[5], max[5], ierflg);

        minuit->mnexcm("SET NOW", arglist, 1, ierflg);
        arglist[0] = 1;
        minuit->mnexcm("SET ERR", arglist, 1, ierflg);
        arglist[0] = 2;
        minuit->mnexcm("SET STR", arglist, 1, ierflg);
        minuit->mnexcm("CALL FCN", arglist, 1, ierflg);

        // Fix some parameters
        minuit->FixParameter(2);
        minuit->FixParameter(3);
        minuit->FixParameter(4);

        // Scan all free parameters to set a new phi window .. maybe more scans are better for fit convergence, thats why i implemented a loop
        for(int l = 0; l<1; l++){
          minuit->mncomd("scan 0 1000", ierflg);
          minuit->GetParameter(5, par[5], err[5]);
          min[5] = par[5] - TMath::Pi();
          max[5] = par[5] + TMath::Pi();
          minuit->mnparm(5, "phase", par[5], err[5], min[5], max[5], ierflg);
        }
        
        minuit->Migrad();
        //save parameters to use them as input start parameters for next fit in next loop iteration
        minuit->mnstat(minfcn, edm, errdef, npari, nparx, istat);
        minuit->GetParameter(0, lastC, best_c_err);
        minuit->GetParameter(1, lastLambda, best_lambda_err);
        minuit->GetParameter(5, lastPhi, err[5]);

        delete minuit;
        minuit = 0;

        L_i = TMath::Exp(-0.5 * (minfcn - best_minfcn));
        //fill my h_scan histogram with the likelihood, normalized to best_minfcn (one histogram is for one mass,width pair)
        h_scan->SetBinContent(k+1, L_i);
      }

      // Calculate upper limit for this hypothesis, so for certain mass,width pair
      Double_t upp_likeli = CalculateUpperLimit(h_scan, 0.9);
      delete h_scan;

        //fill upper limit into the mass,width map - repeat everything for next mass,width pair
      h2d_scan_coh->SetBinContent(i+1, j+1, upp_likeli);
    }
  }

  output_file->cd();
  h2d_scan_coh->Write();
}