#include "TVirtualFitter.h"
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <cmath>

#include <TROOT.h>
#include <TStyle.h>
#include <TRint.h>
#include <TApplication.h>
#include <TMinuit.h>
#include <TLorentzVector.h>
#include <TChain.h>
#include "TTree.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include <iostream>
#include <TMath.h>
#include <fstream>
#include <iomanip>
//#include <"Math/Minimizer.h">
//#include <"Math/Factory.h">
//#include <"Math/Functor">
//#include <"TRandom2.h">
//#include <"TError.h">

using namespace std;


double myFunction(double x1, double x2, double x3,double x4,double x5, double x6, double x7, double x8, double x9, double x10, double x11, double x12,
    double x13, double x14, double x15, double x16, double x17, double x18, double x19, double x20, double x21)  // x1 - x12 (normalizations)  x13 - x21 (parameters)
{

    float norms[12];  //normalization variables (used when reading in data points)
    norms[0]=x1;
    norms[1]=x2;
    norms[2]=x3;
    norms[3]=x4;
    norms[4]=x5;
    norms[5]=x6;
    norms[6]=x7;
    norms[7]=x8;
    norms[8]=x9;
    norms[9]=x10;
    norms[10]=x11;
    norms[11]=x12;

    //////////////////////////////////////////////////////////////// start read in data points
    ifstream inFile;
    inFile.open("allDataForMinuit.txt");

    float allData[1329][7];
    int tempN = 0;
    float final;

    if(!inFile)
    {
       cout << "Unable to open file" << endl;
       exit(0);
    }
    else
    {
        for(int i1 = 0 ; i1 <= 1328; i1++)   // W2, Q2, Sigma_N_e, sigmaSquare, A1data, A1fit=0, expNumber
        {                                    // 0 ,  1,         2,            3,      4,      5,         6
            double bpoints[7];
            for(int b = 0; b <= 6 ; b++)
            {
                inFile >> bpoints[b];
                allData[i1][b] = bpoints[b];
            }
            tempN = bpoints[6];
            allData[i1][6]=norms[tempN]; // assign normalization parameter value based on exp. number
        }
    }
    inFile.close();
    //////////////////////////////////////////////////////////////// end read in data points



    //////////////////////////////////////////////////////////////// start calculate A1fit
    //        x13 x14 x15 x16 x17 x18 x19 x20 x21
    //        p1  p2  p3  p4  p5  p6  p7  p8  p9

    float phoQ2 = 0.27150;
    float pi = 3.14159265;
    float W2;
    float Q2;
    float nu;
    float XI;
    float A;
    float B;
    float C;
    float D;

    for(int i = 0 ; i <= 1328; i++)  // W2, Q2, Sigma_N_e, sigmaSquare, A1data, A1fit=0, expNumber
    {
        W2 = allData[i][0];
        Q2 = allData[i][1];
        nu = (W2 - 0.93827*0.93827 + Q2) / 1.87654;
        XI = (Q2 + phoQ2) / 0.93827 / (nu + sqrt(nu*nu+Q2));
        A = x13 + x14*atan(x15*x15*Q2);
        B = x16 + x17*atan(x18*x18*Q2);
        C = sin( pi*pow(XI,x19) );
        D = x20*pow((1.0-XI),x21);
        allData[i][5] =  pow(XI,A) * (1.0 + B*C + D);
    }
    //////////////////////////////////////////////////////////////// end calculate A1fit



    //////////////////////////////////////////////////////////////// start calculate chisquare

    // W2, Q2, Sigma_N_e, sigmaSquare, A1data, A1fit=0, expNumber/norm
    // 0 ,  1,         2,            3,      4,      5,         6
    float numenator;
    float extraTerm;
    float chiSquare2 = 0.0;

    for(int i = 0 ; i <= 1328; i++)   // W2, Q2, Sigma_N_e, sigmaSquare, A1data, A1fit=0, expNumber
    {
        numenator  = (allData[i][4]-allData[i][6]*allData[i][5])*(allData[i][4]-allData[i][6]*allData[i][5]);
        extraTerm  = (allData[i][6]-1)*(allData[i][6]-1)/(allData[i][2])*(allData[i][2]);
        chiSquare2 = chiSquare2+numenator/allData[i][3]+extraTerm;
        //chiSquare2.Print();
    }
    
    //////////////////////////////////////////////////////////////// end calculate chisquare

   return chiSquare2;


}   // end myFunction







///////////////////////////////////////////////////////////////////////////////////////////////////////////////
      /////////////////////////////          START MINUIT STUFF        ///////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////







void minuitFunction(int& nDim, double* gout, double&result, double par[], int flg)
{
    result = myFunction(par[0],par[1],par[2],par[3],par[4],par[5],par[6], par[7], par[8] , par[9], par[10], par[11],
             par[12],par[13],par[14],par[15],par[16],par[17],par[18], par[19], par[20]);
}




void minuit1()
{
    TFitter* minimizer = new TFitter(21);
    TMinuit *minuit = minimizer -> GetMinuit();
   
    double bestx1 = 1.073, bestx2 = 0.959, bestx3 = 1.291, bestx4 = 1.026, bestx5 = 0.990, bestx6 = 1.218, bestx7 = 0.963;
    double bestx8 = 1.086, bestx9 = 0.881, bestx10 = 0.834, bestx11 = 1.009, bestx12 = 0.963, bestx13 = 2.227, bestx14 = -1.037;
    double bestx15 = -10.212, bestx16 = 1.371, bestx17 = -0.899, bestx18 = -2.089, bestx19 = -0.241, bestx20 = -0.469, bestx21 = 3.656;

    for(int i = 0; i < 21; i++)
    {
       
/*
  TMatrixD matrix0(9,9);
       gMinuit->mnemat(matrix0.GetMatrixArray(),9);
        matrix0.Print();
    
  */
        
        minimizer->SetFCN(minuitFunction);
        
                ////set print level ==========================
                ////arglist[0] = -1 for quiet mode
            ////arglist[0] =  1 for verbose mode
        double p1[10];
        {
            p1[0] = 1;
            minimizer->ExecuteCommand("SET PRINT",p1,1);
        }
        ROOT::Math::Minimizer* min = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
        ROOT::Math::MinimizerOptions::SetDefaultMaxFunctionCalls(100000000);
        ROOT::Math::MinimizerOptions::SetDefaultMaxIterations(100000);
        ROOT::Math::MinimizerOptions::SetDefaultTolerance(0.001);
        ROOT::Math::MinimizerOptions::SetDefaultEPS(0.001);
        ROOT::Math::MinimizerOptions::SetDefaultPrintLevel(3);
        
        
       
//parameter number, parameter name, initial parameter value, step size, lower limit for the parameter, upper limit for the parameter
               minimizer->SetParameter(0,"x1",bestx1,0.00001,0,2.146);
               minimizer->SetParameter(1,"x2",bestx2,0.00001,0,1.918);
               minimizer->SetParameter(2,"x3",bestx3,0.00001,0,2.582);
               minimizer->SetParameter(3,"x4",bestx4,0.00001,0,2.052);
               minimizer->SetParameter(4,"x5",bestx5,0.00001,0,1.98);
               minimizer->SetParameter(5,"x6",bestx6,0.00001,0,2.436);
               minimizer->SetParameter(6,"x7",bestx7,0.00001,0,1.926);
               minimizer->SetParameter(7,"x8",bestx8,0.00001,0,2.172);
               minimizer->SetParameter(8,"x9",bestx9,0.00001,0,1.762);
               minimizer->SetParameter(9,"x10",bestx10,0.00001,0,1.668);
               minimizer->SetParameter(10,"x11",bestx11,0.00001,0,2.018);
               minimizer->SetParameter(11,"x12",bestx12,0.00001,0,1.926);
               minimizer->SetParameter(12,"x13",bestx13,0.00001,0,4.454);
               minimizer->SetParameter(13,"x14",bestx14,0.00001,-2.074,0);
               minimizer->SetParameter(14,"x15",bestx15,0.00001,-20.424,0);
               minimizer->SetParameter(15,"x16",bestx16,0.00001,0,2.742);
               minimizer->SetParameter(16,"x17",bestx17,0.00001,-1.798,0);
               minimizer->SetParameter(17,"x18",bestx18,0.00001,-4.178,0);
               minimizer->SetParameter(18,"x19",bestx19,0.00001,-0.482,0);
               minimizer->SetParameter(19,"x20",bestx20,0.00001,-0.938,0);
               minimizer->SetParameter(20,"x21",bestx21,0.00001,0,7.312);


      minimizer->ExecuteCommand("SIMPLEX",p1,0);
        minimizer->ExecuteCommand("MINIMIZE",p1,0);
     // minimizer->ExecuteCommand("MINOS",0,0);
        bestx1 = minimizer->GetParameter(0);
        bestx2 = minimizer->GetParameter(1);
        bestx3 = minimizer->GetParameter(2);
        bestx4 = minimizer->GetParameter(3);
        bestx5 = minimizer->GetParameter(4);
        bestx6 = minimizer->GetParameter(5);
        bestx7 = minimizer->GetParameter(6);
        bestx8 = minimizer->GetParameter(7);
        bestx9 = minimizer->GetParameter(8);
        bestx10 = minimizer->GetParameter(9);
        bestx11 = minimizer->GetParameter(10);
        bestx12 = minimizer->GetParameter(11);
        bestx13 = minimizer->GetParameter(12);
        bestx14 = minimizer->GetParameter(13);
        bestx15 = minimizer->GetParameter(14);
        bestx16 = minimizer->GetParameter(15);
        bestx17 = minimizer->GetParameter(16);
        bestx18 = minimizer->GetParameter(17);
        bestx19 = minimizer->GetParameter(18);
        bestx20 = minimizer->GetParameter(19);
        bestx21 = minimizer->GetParameter(20);
        
     double amin, chi2, edm, errdef; // edm = Expected distance from minimum
        int nvpar, nparx, icstat;
    minimizer->GetStats(chi2,edm,errdef,nvpar,nparx);
       // minimizer->PrintResults(21,chi2);
       /* double call, init_chisq, chisq, result;
        if(call==1)
        {
            init_chisq = chisq;
            cout<<"\n Initial chisq: "<<chisq<<"\n\n";
        }
       */
       // TFitterMinuit::GetCovarianceMatrix();
   TMatrixD matrix0(9,9);
             gMinuit->mnemat(matrix0.GetMatrixArray(),9);
              matrix0.Print();
        
    //  double err_mat[9][9];
    //   gMinuit->mnemat(&err_mat[0][0],9);*/
        
       /*
        double call, init_chisq, chisq, result;
        
           if(call==1)
           {
               init_chisq = chisq;
               cout<<"\n Initial chisq: "<<chisq<<"\n\n";
           }
          // if(call%10==0) cout<<" chisq="<<result;
        */
       

        
        cout << "Results: " << endl;
        cout << bestx1 << endl;
        cout << bestx2 << endl;
        cout << bestx3 << endl;
        cout << bestx4 << endl;
        cout << bestx5 << endl;
        cout << bestx6 << endl;
        cout << bestx7 << endl;
        cout << bestx8 << endl;
        cout << bestx9 << endl;
        cout << bestx10 << endl;
        cout << bestx11 << endl;
        cout << bestx12 << endl;
        cout << bestx13 << endl;
        cout << bestx14 << endl;
        cout << bestx15 << endl;
        cout << bestx16 << endl;
        cout << bestx17 << endl;
        cout << bestx18 << endl;
        cout << bestx19 << endl;
        cout << bestx20 << endl;
        cout << bestx21 << endl;
        cout<< chi2 << '\n';
//cout << minimizer->GetCovarianceMatrix(9,9)  << endl;
 
    }
}