#ifndef ULMAPS_TEST_H_
#define ULMAPS_TEST_H_
// C++ libraries
#include <cstdio>

// Other libraries
#include <TROOT.h>
#include <TStyle.h>

// C++ libraries
#include <stdio.h>
#include <iostream>
#include <sys/stat.h>

// Other libraries
#include <TH1.h>
#include <TH2.h>
#include <TLegend.h>
#include <TCanvas.h>
#include <TChain.h>
#include <TFile.h>
#include <TLine.h>
#include <TGraphErrors.h>
#include <TGraphAsymmErrors.h>
#include <TH1.h>
#include <TH2.h>
#include <TLorentzVector.h>
#include <TMultiGraph.h>
#include <TROOT.h>
#include <TSelector.h>
#include <TGaxis.h>
#include <time.h>
#include <TMath.h>
#include <TPaletteAxis.h>
#include <TF1Convolution.h>
#include <TStyle.h>
#include <TFitResultPtr.h>
#include <TFitResult.h>
#include "TVirtualFitter.h"
#include "../../myPDGmasses.h"
#include <TArrow.h>
#include <TMatrixTSym.h>
#include <TMatrixD.h>
#include <stdio.h>
#include <Math/MinimizerOptions.h>
#include <TF1Convolution.h>
#include <TComplex.h>
#include <RooHistError.h>
#include <TArc.h>

#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
#include <TSelector.h>
#include <TTree.h>
#include <TTreeReader.h>
#include <TTreeReaderValue.h>
#include <TTreeReaderArray.h>
#include <TLegend.h>
#include <TLatex.h>
#include <TPad.h>
#include <TEllipse.h>
#include <string.h>
#include <assert.h>
#include "Fit/Fitter.h"
#include "Fit/BinData.h"
#include "Fit/Chi2FCN.h"
#include "Fit/PoissonLikelihoodFCN.h"
#include "Fit/LogLikelihoodFCN.h"
#include "TList.h"
#include "Math/WrappedMultiTF1.h"
#include "HFitInterface.h"
#include <stdarg.h>
#include <iostream>
#include <fstream>
#include <sstream>


TFile * output_file;

void ulmaps_test(Int_t, Int_t);
void ULcoh(Int_t, Int_t);
void MinuitFCN_ContinuumBWcoh(int&, double *, double&, double *, int);
const int nPostana = 15;

double ecms[nPostana] = {3.9,3.95,4.0,4.05,4.1,4.15,4.2,4.25,4.3,4.35,4.4,4.45,4.5,4.55,4.6};
double acc[nPostana] = {0.3,0.29,0.28,0.27,0.26,0.25,0.24,0.23,0.22,0.21,0.20,0.19,0.18,0.17,0.16};
double mean[nPostana] = {100,95,85,80,75,70,72,76,73,65,60,55,50,45,40};
double sig_left[nPostana] = {100*0.1,95*0.1,85*0.1,80*0.1,75*0.1,70*0.1,72*0.1,76*0.1,73*0.1,65*0.1,60*0.1,55*0.1,50*0.1,45*0.1,40*0.1};
double sig_right[nPostana] = {100*0.1,95*0.1,85*0.1,80*0.1,75*0.1,70*0.1,72*0.1,76*0.1,73*0.1,65*0.1,60*0.1,55*0.1,50*0.1,45*0.1,40*0.1};


Double_t F_Continuum_BW_coh(Double_t * x, Double_t * par)
{
  Double_t E = x[0];
  Double_t mass  = par[2];
  Double_t width = par[3];
  
  TComplex bw = TComplex(0, 0);

  bw = mass * width / TComplex(E * E - mass * mass, mass * width);

  TComplex phase(cos(par[5]), sin(par[5]));
  TComplex se(0, 0);
  se += TMath::Sqrt(pow(par[0] / x[0], par[1]));
  
  if (par[4] != 0.) { se += TMath::Sqrt(par[4]) * bw * phase; }
  Double_t f = se.Rho2();
  return f;
}

void MinuitFCN_ContinuumBWcoh(int& nDim, double * gout, double& result, double * par, int flg)
{  
  Double_t L = 1.;
  Double_t x[1] = {0.};
  Double_t lumi = 500;
  for (int i = 0; i < nPostana; i++)
  {
    x[0] = ecms[i];
    Double_t f_sigma = F_Continuum_BW_coh(x, par);
    Double_t br = 0.492;
    Double_t n = f_sigma * lumi * acc[i] * br;
    
    Double_t G = 1.; 
   
    if (n < mean[i]){
      G *= 1. / (TMath::Sqrt(2. * TMath::Pi()) * TMath::Sqrt(pow(sig_left[i],2)));
      G *= TMath::Exp(-0.5 * pow(n - mean[i] , 2.) / (pow(sig_left[i],2.)));
    } 
    else{
      G *= 1. / (TMath::Sqrt(2. * TMath::Pi()) * TMath::Sqrt(pow(sig_right[i],2)));
      G *= TMath::Exp(-0.5 * pow(n - mean[i] , 2.) / (pow(sig_right[i],2.)));
    }
    L *= G;
  }
  result = -2. * TMath::Log(L);
}


TH1D * NewTH1D(TString name,
                Int_t bins_x, Double_t x_min, Double_t x_max,
                TString x_entry)
 {
   Double_t x_binning = (x_max - x_min) / bins_x;
   TH1D * h1d = new TH1D(name, "", bins_x, x_min, x_max);
   h1d->GetXaxis()->SetTitle(Form("%s", x_entry.Data()));
   h1d->GetYaxis()->SetTitle(Form("dN/d%s / (%g)^{-1}", x_entry.Data(), x_binning));

   /*
   if (x_binning == 1.)
   {
     h1d->GetXaxis()->SetTitle(Form("%s", x_entry.Data()));
     h1d->GetYaxis()->SetTitle(Form("dN/d%s", x_entry.Data()));
   }
   else
   {
     h1d->GetXaxis()->SetTitle(Form("%s", x_entry.Data()));
     h1d->GetYaxis()->SetTitle(Form("dN/d%s / (%g)", x_entry.Data(), x_binning));
   }
   */
   return h1d;
 }
TH1D * NewTH1D(TString name,
               Int_t bins_x, Double_t x_min, Double_t x_max,
               TString x_entry, TString x_unit)
{
  Double_t x_binning = (x_max - x_min) / bins_x;
  TH1D * h1d = new TH1D(name, "", bins_x, x_min, x_max);
  h1d->GetXaxis()->SetTitle(Form("%s / (%s)", x_entry.Data(), x_unit.Data()));
  h1d->GetYaxis()->SetTitle(Form("dN/d%s / (%g %s)^{-1}", x_entry.Data(), x_binning, x_unit.Data()));
  return h1d;
}
TH2D * NewTH2D(TString name,
               Int_t bins_x, Double_t x_min, Double_t x_max,
               Int_t bins_y, Double_t y_min, Double_t y_max,
               TString x_label, TString y_label, TString z_label = "",
               TString draw_options = "COLZ")
{
  TH2D *h2d = new TH2D(name, "", bins_x, x_min, x_max, bins_y, y_min, y_max);
  h2d->GetXaxis()->SetTitle(x_label);
  h2d->GetYaxis()->SetTitle(y_label);
  h2d->GetZaxis()->SetTitle(z_label);
  h2d->SetOption(draw_options);
  return h2d;
}

Double_t CalculateUpperLimit(TH1 * h, Double_t conf)
{
  Int_t min = h->GetXaxis()->FindBin(0.);
  Int_t max = h->GetNbinsX();
  Double_t integral = h->Integral(min, max);

  Int_t n = 0;
  Double_t value = 0.;
  while(value < conf)
  {
    n++;
    value = h->Integral(min, n) / integral;
  }
  return h->GetBinCenter(n);
}

void PrintProcessStatus(double start_time, double current_time,
                        long long entry, long long number_of_entries,
                        int bar_width = 50, int update_frequency = 0,
                        bool print_eta = true)
{
  if (entry % (number_of_entries / bar_width) == 0 ||
      entry == number_of_entries - 1 ||
      (update_frequency && entry % update_frequency == 0))
  {
    Double_t progress = (entry + 1.) / number_of_entries;
    Int_t position = (Int_t)(bar_width * progress);

    printf("|");
    for (int i = 0; i < bar_width; i++)
    {
      if (i < position)
      {
        printf("█");
      }
      else
      {
        printf(" ");
      }
    }
    printf("| %3i%%", (Int_t)(progress * 100));

    Int_t eta = 0;
    Double_t dt = current_time - start_time;
    if (dt)
    {
      Double_t v = entry / (current_time - start_time);
      eta = (Int_t)((number_of_entries - (entry + 1)) / v);
    }

    if (print_eta)
    {
      if (eta)
      {
        Int_t eta_hours = eta / 3600;
        Int_t eta_minutes = (eta % 3600) / 60;
        Int_t eta_seconds = ((eta % 3600) % 60);

        printf("  ETA:");
        if (eta_hours)
        {
          printf(" %ih", eta_hours);
        }
        if (eta_minutes)
        {
          printf(" %imin", eta_minutes);
        }
        if (eta_seconds)
        {
          printf(" %is", eta_seconds);
        }
      }
      printf("                     ");
    }

    if (entry == number_of_entries - 1)
    {
      printf("\n");
    }
    else
    {
      printf("\r");
      fflush(stdout);
    }
  }
}

#endif