//LMFH

#include <iostream>
#include <fstream>

#include "TMinuit.h"
#include "TMath.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "Math/Integrator.h"
#include "TGraphErrors.h"
#include "Math/GSLIntegrator.h"
#include "Math/IFunction.h"

#include "Math/Polynomial.h"
#include "Math/Interpolator.h"

#include <gsl/gsl_integration.h>


const Int_t theSize = 27;
Float_t R_half_arcmin[theSize];
Float_t R_half_err[theSize];
Float_t R_bin_arcmin[theSize];
Float_t R_bin_arcmin_err[theSize];
Float_t Mean_V[theSize];
Float_t DeltaV3_err[theSize];
Float_t Sigmalos[theSize];
Float_t Sigmalos_err[theSize];
Float_t DeltaV5_err[theSize];
Float_t Totalstars_bin[theSize];
Float_t Members_bin[theSize];
Float_t R_pc[theSize];
Float_t R_pc_err[theSize];
Float_t Standardeviation_pc_bin[theSize];



Double_t RhoDM(float x, Double_t *par)
{
    Double_t value1 = par[0]/((1.-TMath::Pi())*(1.+TMath::Power((x/par[1]),2)));
        return value1;
}

Double_t MDM(float x, Double_t *par)
{
    Double_t value2 = (par[0]/(1.-(TMath::Pi())/4))*((x/par[1]) - TMath::ATan(x/par[1]));
    return value2;
}


Double_t f1(float x,  Double_t *par){

    Double_t value3 =(x/par[1] - TMath::ATan(x/par[1]))*(TMath::Sqrt(TMath::Pi())*TMath::Gamma(-0.5+par[2])/TMath::Gamma(par[2]) - TMath::BetaIncomplete(TMath::Power(1./x,2), -0.5 + par[2], 0.5) + par[2]*TMath::BetaIncomplete(TMath::Power(1/x,2), 0.5 + par[2], 0.5)  -TMath::Sqrt(TMath::Pi())*TMath::Gamma(0.5+par[2])/TMath::Gamma(par[2]))/TMath::Power(1.0 + (TMath::Power(x/par[3],2)), 2.5);
    return value3;
}


double Integrand(double x)
{
  
   return ROOT::Math::GSLIntegrator::IntegralUp(&f1,1) ;
    }




Double_t sigma(float R, Double_t *par){

    Double_t par1[3]={par[1]/R, par[2], par[3]/R};
    Double_t value7 = (1.5/par[3])*(TMath::Power(1.0 + (TMath::Power(R/par[3],2)), 2))*(par[0]/(1.- TMath::Pi()/4))*(0.5*TMath::Power(R,2.-2.*par[2]))*Integrand(par1);
    return value7;
}



Double_t func(float R, Double_t *par)
{
 
    Double_t value = TMath::Sqrt((sigma(R,par)));
    return value;
}

void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
    const Int_t nbins = theSize;
    Int_t i;
    
    // calculate chisquare
    Double_t chisq = 0;
    Double_t delta;
    for (i = 0; i < nbins; i++) {
        delta  = (Sigmalos[i] - func(R_pc[i], par))/Sigmalos_err[i];
        chisq += delta*delta;
    }
    f = chisq;
}

Double_t profile(Double_t *R, Double_t *par)
{
    return func(R[0], par);
}

void fit() {
    double valueR_half_arcmin,valueR_bin_arcmin,valueMean_V,valueDeltaV3_err,valueSigmalos,valueDeltaV5_err,valueTotalstars_bin,valueMembers_bin,valueR_pc,valueStandardeviation_pc_bin;
    std::ifstream theInput("car.txt");
    unsigned int theIndex = 0;
    while (theInput >> valueR_half_arcmin >> valueR_bin_arcmin >> valueMean_V >> valueDeltaV3_err >> valueSigmalos >> valueDeltaV5_err >> valueTotalstars_bin >> valueMembers_bin >> valueR_pc >> valueStandardeviation_pc_bin) {
        R_half_arcmin[theIndex]     = valueR_half_arcmin;
        R_half_err[theIndex]     = 0.;
        R_bin_arcmin[theIndex] = valueR_bin_arcmin;
        R_bin_arcmin_err[theIndex] = 0;
        Mean_V[theIndex] = valueMean_V;
        DeltaV3_err[theIndex] = valueDeltaV3_err;
        R_pc[theIndex] = valueR_pc;
        R_pc_err[theIndex] = 0;
        Sigmalos[theIndex] = valueSigmalos;
        Sigmalos_err[theIndex] = valueDeltaV5_err;
        Totalstars_bin[theIndex] = valueTotalstars_bin;
        Members_bin[theIndex] = valueMembers_bin;
        Standardeviation_pc_bin[theIndex] = valueStandardeviation_pc_bin;
        theIndex++;
    }

    
    TMinuit *gMinuit = new TMinuit(4);
    gMinuit->SetFCN(fcn);
    
    Double_t arglist[10];
    Int_t ierflg = 0;
    
    arglist[0] = 1;
    gMinuit->mnexcm("SET ERR", arglist , 1, ierflg);
    
    // Set starting values and step sizes for parameters
    
    static Double_t step[4] = {0.01 , 0.01, 0.01, 0.01};
    static Double_t vstart[4] = {1000.0, 1000., 1.0, 1000.};
    std::vector<std::string> parName = {"Ms", "rs", "Beta","rh"};
    vstart[3]=241.0;
    static Double_t vMin[4] = {1.0, 1., -1.0, 0.0 };
    static Double_t vMax[4] = {100000, 10000., 1.0, 250.};
    
    for (theIndex = 0; theIndex < 4; theIndex++) {
        gMinuit->mnparm(theIndex, parName.at(theIndex).c_str(), vstart[theIndex], step[theIndex], vMin[theIndex], vMax[theIndex], ierflg);
    }
    
    gMinuit->FixParameter(3);
    
//    static Double_t vstart[3] = {0.0001, 100., 0.1};
//    static Double_t step[3] = {0.1 , 0.1, 0.1};
//    gMinuit->mnparm(0, "rsch", vstart[0], step[0], 0,0, ierflg);
//    gMinuit->mnparm(1, "rho_sch", vstart[1], step[1], 0,0, ierflg);
//    gMinuit->mnparm(2, "Rc", vstart[2], step[2], 0,0, ierflg);
    
    // Now ready for minimization step
    arglist[0] = 10000;
    arglist[1] = 1.;
    gMinuit->mnexcm("MIGRAD", arglist, 1, ierflg);
    
    // Print results
    Double_t amin,edm,errdef;
    Int_t nvpar,nparx,icstat;
    
    gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
    gMinuit->mnprin(3, amin);
    
    double Ms, Ms_err;
    double rs, rs_err;
    double Beta, Beta_err;
    double rh, rh_err;
    gMinuit->GetParameter(0, Ms, Ms_err);
    gMinuit->GetParameter(1, rs, rs_err);
    gMinuit->GetParameter(2, Beta, Beta_err);
    gMinuit->GetParameter(3, rh, rh_err);
    
    if (abs(Ms) < 3.*Ms_err )
        std::cout << "\033[1;31mrsch = " << Ms << "+/-" << Ms_err << "\033[0m" << std::endl;
    if (abs(rs) < 3.*rs_err )
        std::cout << "\033[1;31mrho_sch = " << rs << "+/-" << rs_err << "\033[0m" << std::endl;
    if (abs(Beta) < 3.*Beta_err )
        std::cout << "\033[1;31mRc = " << Beta << '\t' << Beta_err << "\033[0m" << std::endl;

    
    
    std::cout << Ms << '\t' << Ms_err << std::endl;
    std::cout << rs << '\t' << rs_err << std::endl;
    std::cout << Beta << '\t' << Beta_err << std::endl;
    std::cout << rh << '\t' << rh_err << std::endl;
    
    TF1 *f = new TF1("f1", profile, 29.4, 865.2, 4);
    f->SetLineColor(kRed);
    f->SetNpx(500);
    f->SetParameters(Ms, rs, Beta, rh);
    
    
    double theChiSquare = 0.0;
    double temp;
    for (unsigned int index = 0; index < theSize; index++) {
        temp = (Sigmalos[index] - f->Eval(R_pc[index]))/Sigmalos_err[index];
        theChiSquare += temp*temp;
    }
    std::cout  << "Ms=" << Ms << "+-" << std::setprecision(4) << Ms_err << "," << std::setprecision(4) << "Beta=" << Beta << "+-" << std::setprecision(4) << Beta_err << ","
    << std::setprecision(4)<< "rs="<< (0.000232407/(4.0*(TMath::Pi())))*rs  << "+-" << std::setprecision(2)<< (0.000232407/(4.0*(TMath::Pi())))*rs_err  << std::setprecision(4) << ","
    << "chi^2=" << theChiSquare << std::setprecision(4) << ","
    << "chi^2_red=" << theChiSquare/(theSize - 3.0)
    << std::endl;
    
    TCanvas *c1 = new TCanvas("c1", "Data points", 1000, 600);
    
    const Int_t n = theSize;
    
    TGraphErrors *gr = new TGraphErrors(n, R_pc, Sigmalos, R_pc_err, Sigmalos_err);
    gr->SetTitle("CamB.dat");
    gr->SetMarkerColor(4);
    gr->SetMarkerStyle(21);
    gr->Draw("AP");
    f->Draw("same");
    c1->SaveAs("vcCamB.pdf");


}