#include <fstream>
#include <string>
#include <sstream>
#include <iostream>
#include <iomanip>
#include <stdio.h>
#include <stdlib.h> // for exit function

#include "TMinuit.h"
#include "TROOT.h"

using std::cerr;
using std::cout;
using std::endl;
using std::ifstream;
using std::ofstream;

using namespace std;

//______________________________________________________________________________
//Cross section energy foils - energy groups

//double cs[nfoils][nEnGroups];
string fileList="/Users/nelli/sw/root/list_files.txt";
string fileList2="/Users/nelli/sw/root/list_param.txt";
//string fileList[10000] ;

int num_files;
int numline=0;
int numline2=0;
int numline3=0;
int numline4=0;
string fname, line;
//ifstream indata, indata2,indata3;
string fileNames[100];
string fileNames_par[100];
Float_t cs[15][73], Ameas_s[15], Error_s[15];//, Acalc_Guess[15][1];

int numparam=0;
int numparam2=0;





    //______________________________________________________________________________
void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
    //calculate chisquare
    Double_t chisq = 0;
    Double_t  delta =0;
    Double_t dd = 0.;
    
    for (int i = 0; i < 15 ; i++ )       // i :  number of reactions (15)
    {
        dd = 0.;
        for (int j = 0; j <= 72 ; j++ )   // j : number of energy groups (73)
        {
            dd +=cs[i][j]*par[j];
        }
        delta = 0;
        delta =(Ameas_s[i]-dd)/Error_s[i];
        chisq +=delta*delta;
    }
        f = chisq;  
}
  

//______________________________________________________________________________
void seghour()
{
    ifstream indata, indata2,indata3,indata4,indata5;
    // read file list
    
    indata.open("list_files.txt");
    while (! indata.eof() )
    {
        getline (indata,line);
        cout << line << endl;
        numline++;
    }
    
    indata.close();
    indata.clear();
    
    cout << "number of files: "<<numline << endl;
    cout << "=====================================" << endl;
    
    indata.open("list_files.txt");
    for (int i = 0; i < numline ; i++ ) // i number of foils
    {
        indata >> fileNames[i];
        cout << " read cross sections from "<<fileNames[i] << endl;
        // start reading here
        indata2.open(fileNames[i].c_str());
        while (! indata2.eof() )
        {
            getline (indata2,line);
            numline2++;
        }
        indata2.close();
        indata2.clear();
        cout << "number of lines" <<" "<<numline2 << endl;
        indata2.open(fileNames[i].c_str());
        for (int j = 0; j < numline2 ; j++ ) // j number of energy - groups
        {
            indata2 >> cs[i][j];
           // cout << cs[i][j]<< endl;
        }
    }



       
//***********************************a_meas********************************************
    
    indata3.open("Ameas_s.txt");
    
    if (indata3.is_open () ){
        while (! indata3.eof() ){
            getline (indata3,line);
            numline3++;}
        indata3.close();
        indata3.clear();
        cout << "number of lines in measured activity file " << numline3 << endl;
        indata3.open("Ameas_s.txt");
        
        for (int ir = 0; ir < numline3 ; ir++ ){
            indata3 >> Ameas_s[ir];
            cout << Ameas_s[ir] << "\n";}
    }
    
    indata4.open("Error_s.txt");
    
    if (indata4.is_open () ){
        while (! indata4.eof() ){
            getline (indata4,line);
            numline4++;}
        indata4.close();
        indata4.clear();
        cout << "number of lines in error file " << numline4 << endl;
        indata4.open("Error_s.txt");
        
        for (int iir = 0; iir < numline4 ; iir++ ){
            indata4 >> Error_s[iir];
            cout << Error_s[iir] << "\n";}
    }

    //*****************************************error********************************************
  


//*************************************************************************************

    
    TMinuit *gMinuit = new TMinuit(73);  //initialize TMinuit with a maximum of 15 params
   
    gMinuit->SetFCN(fcn);
    
    //select verbose level:
    //-1 : minimum
    // 0 : low
    // 1 : medium
    // 2 : high
    // 3 : maximum
    
   // gMinuit->SetPrintLevel(1);
    Double_t arglist[10];
    Int_t ierflg = 0;
    
    //Minimization strategy
    // 1 standard
    // 2 try to improve minimum (slower)
    
    arglist[0] = 1;
    gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
    //gMinuit->mnexcm("SET NOW", arglist ,1,ierflg);
    //gMinuit->mnexcm("SET STR", arglist, 1, ierflg) ;
    //arglist [0] = 1;
    //gMinuit->mnexcm("FIX ", arglist ,1,ierflg);
    //arglist[0] = 0;
    //gMinuit->mnexcm("SCAN", arglist,1,ierflg);
    //gMinuit->mnexcm("CALL FCN", arglist , 1,ierflg);
    
    //_____________________________________________________
    
   /*
    
    indata4.open("list_param.txt");
    while (! indata4.eof() )
    {
        getline (indata4,line);
        cout << line << endl;
        numparam++;
    }
    
    indata4.close();
    indata4.clear();
    
    cout << "number of files in list_param folder: "<<numparam << endl;
    cout << "=====================================" << endl;
    
    indata4.open("list_param.txt");
    for (int ii = 0; ii < numparam ; ii++ ) // ii number of files in folder
    {
        indata4 >> fileNames_par[ii];
        cout << "read parameters from "<<fileNames_par[ii] << endl;
        // start reading here
        indata5.open(fileNames_par[ii].c_str());
        while (! indata5.eof() )
        {
            getline (indata5,line);
            numparam2++;
            cout << numparam2 << endl;
        }
        indata5.close();
        indata5.clear();
        cout << "number of lines" <<" "<<numparam2 << endl;
        indata5.open(fileNames_par[ii].c_str());
        for (int jj = 0; jj < 15 ; jj++ ) //
        {
            indata5 >> Acalc_Guess[jj][ii];
            cout << "parameters: " << Acalc_Guess[jj][ii] << endl;
        }
    }

*/
        //_____________________________________________________
        
   
    // Set starting values and step sizes for parameters
   
    static Double_t Acalc_Guess[73] = {1.07E+08,	1.22E+08,	2.22e+08,	3.0e+08,	4.18E+08,	8.80E+08,	1.33E+09,	6.29E+08,	2.74E+07,	166.05,	61.05,	26.3	,15.575,	10.61,	7.81,	6.065,	4.89,	4.05,	3.425,	1.93,	0.83225,	0.4925,	0.335,	0.247,	0.192,	0.15475,	0.12825,	0.10825,	0.06105,	0.0263,	0.015575,	0.01061,	0.00781,	0.006065	,0.00489	,0.00405	,0.003425	,0.001871,	0.00083225	,0.0004925,	0.000335	,0.000247,	0.000192	,0.00015475,	0.00012825,	0.00010825,	5.92E-05	,2.63E-05	,1.56E-05,	1.06E-05	,7.81E-06,	6.07E-06,	4.89E-06	,4.05E-06,	3.43E-06	,4.89E-06	,5.96E-06	,6.60E-06	,6.99E-06	,7.22E-06,	7.33E-06,	7.34E-06	,7.29E-06	,7.19E-06,	6.08E-06	,3.84E-06	,2.16E-06	,1.14E-06	,5.79E-07,	2.85E-07,	1.37E-07,	6.50E-08,	1.81E-08};
   
    static Double_t step[73] = {100,	100,	100,	100,	100	, 100,	100	,100	,100,	100	,100,	100	,100,100	,100,  100,	100	,100,	100,	100,	100,	100,	100,	100,	100,	100,	100,100,100	,100,	100,	100,	100,       100,	100,    100,	100	,100	,100,	100	,100,	100,100	,100,100,	100,	100,   100,  	100,	100,	100	,100,	100	,100,	100,	100, 100,	100,100,100,100,   100,	100,100,	100,	    100,	100,	100,	100,	100, 100,	100	,100};
    
    gMinuit->mnparm(0, "a1", Acalc_Guess[0], step[0], 0,0,ierflg);
    gMinuit->mnparm(1, "a2", Acalc_Guess[1], step[1], 0,0,ierflg);
    gMinuit->mnparm(2, "a3", Acalc_Guess[2], step[2], 0,0,ierflg);
    gMinuit->mnparm(3, "a4", Acalc_Guess[3], step[3], 0,0,ierflg);
    gMinuit->mnparm(4, "a5", Acalc_Guess[4], step[4], 0,0,ierflg);
    gMinuit->mnparm(5, "a6", Acalc_Guess[5], step[5], 0,0,ierflg);
    gMinuit->mnparm(6, "a7", Acalc_Guess[6], step[6], 0,0,ierflg);
    gMinuit->mnparm(7, "a8", Acalc_Guess[7], step[7], 0,0,ierflg);
    gMinuit->mnparm(8, "a9", Acalc_Guess[8], step[8], 0,0,ierflg);
    gMinuit->mnparm(9, "a10", Acalc_Guess[9], step[9], 0,0,ierflg);
    gMinuit->mnparm(10, "a11", Acalc_Guess[10], step[10], 0,0,ierflg);
    gMinuit->mnparm(11, "a12", Acalc_Guess[11], step[11], 0,0,ierflg);
    gMinuit->mnparm(12, "a13", Acalc_Guess[12], step[12], 0,0,ierflg);
    gMinuit->mnparm(13, "a14", Acalc_Guess[13], step[13], 0,0,ierflg);
    gMinuit->mnparm(14, "a15", Acalc_Guess[14], step[14], 0,0,ierflg);
    
    gMinuit->mnparm(15, "a16", Acalc_Guess[15], step[15], 0,0,ierflg);
    gMinuit->mnparm(16, "a17", Acalc_Guess[16], step[16], 0,0,ierflg);
    gMinuit->mnparm(17, "a18", Acalc_Guess[17], step[17], 0,0,ierflg);
    gMinuit->mnparm(18, "a19", Acalc_Guess[18], step[18], 0,0,ierflg);
    gMinuit->mnparm(19, "a20", Acalc_Guess[19], step[19], 0,0,ierflg);
    gMinuit->mnparm(20, "a21", Acalc_Guess[20], step[20], 0,0,ierflg);
    gMinuit->mnparm(21, "a22", Acalc_Guess[21], step[21], 0,0,ierflg);
    gMinuit->mnparm(22, "a23", Acalc_Guess[22], step[22], 0,0,ierflg);
    gMinuit->mnparm(23, "a24", Acalc_Guess[23], step[23], 0,0,ierflg);
    gMinuit->mnparm(24, "a25", Acalc_Guess[24], step[24], 0,0,ierflg);
    gMinuit->mnparm(25, "a26", Acalc_Guess[25], step[25], 0,0,ierflg);
    gMinuit->mnparm(26, "a27", Acalc_Guess[26], step[26], 0,0,ierflg);
    gMinuit->mnparm(27, "a28", Acalc_Guess[27], step[27], 0,0,ierflg);
    gMinuit->mnparm(28, "a29", Acalc_Guess[28], step[28], 0,0,ierflg);
    gMinuit->mnparm(29, "a30", Acalc_Guess[29], step[29], 0,0,ierflg);
    
    gMinuit->mnparm(30, "a31", Acalc_Guess[30], step[30], 0,0,ierflg);
    gMinuit->mnparm(31, "a32", Acalc_Guess[31], step[31], 0,0,ierflg);
    gMinuit->mnparm(32, "a33", Acalc_Guess[32], step[32], 0,0,ierflg);
    gMinuit->mnparm(33, "a34", Acalc_Guess[33], step[33], 0,0,ierflg);
    gMinuit->mnparm(34, "a35", Acalc_Guess[34], step[34], 0,0,ierflg);
    gMinuit->mnparm(35, "a36", Acalc_Guess[35], step[35], 0,0,ierflg);
    gMinuit->mnparm(36, "a37", Acalc_Guess[36], step[36], 0,0,ierflg);
    gMinuit->mnparm(37, "a38", Acalc_Guess[37], step[37], 0,0,ierflg);
    gMinuit->mnparm(38, "a39", Acalc_Guess[38], step[38], 0,0,ierflg);
    gMinuit->mnparm(39, "a40", Acalc_Guess[39], step[39], 0,0,ierflg);
    gMinuit->mnparm(40, "a41", Acalc_Guess[40], step[40], 0,0,ierflg);
    gMinuit->mnparm(41, "a42", Acalc_Guess[41], step[41], 0,0,ierflg);
    gMinuit->mnparm(42, "a43", Acalc_Guess[42], step[42], 0,0,ierflg);
    gMinuit->mnparm(43, "a44", Acalc_Guess[43], step[43], 0,0,ierflg);
    gMinuit->mnparm(44, "a45", Acalc_Guess[44], step[44], 0,0,ierflg);
    gMinuit->mnparm(45, "a46", Acalc_Guess[45], step[45], 0,0,ierflg);
    
    gMinuit->mnparm(46, "a47", Acalc_Guess[46], step[46], 0,0,ierflg);
    gMinuit->mnparm(47, "a48", Acalc_Guess[47], step[47], 0,0,ierflg);
    gMinuit->mnparm(48, "a49", Acalc_Guess[48], step[48], 0,0,ierflg);
    gMinuit->mnparm(49, "a50", Acalc_Guess[49], step[49], 0,0,ierflg);
    gMinuit->mnparm(50, "a51", Acalc_Guess[50], step[50], 0,0,ierflg);
    gMinuit->mnparm(51, "a52", Acalc_Guess[51], step[51], 0,0,ierflg);
    gMinuit->mnparm(52, "a53", Acalc_Guess[52], step[52], 0,0,ierflg);
    gMinuit->mnparm(53, "a54", Acalc_Guess[53], step[53], 0,0,ierflg);
    gMinuit->mnparm(54, "a55", Acalc_Guess[54], step[54], 0,0,ierflg);
    gMinuit->mnparm(55, "a56", Acalc_Guess[55], step[55], 0,0,ierflg);
    gMinuit->mnparm(56, "a57", Acalc_Guess[56], step[56], 0,0,ierflg);
    gMinuit->mnparm(57, "a58", Acalc_Guess[57], step[57], 0,0,ierflg);
    gMinuit->mnparm(58, "a59", Acalc_Guess[58], step[58], 0,0,ierflg);
    gMinuit->mnparm(59, "a60", Acalc_Guess[59], step[59], 0,0,ierflg);
    gMinuit->mnparm(60, "a61", Acalc_Guess[60], step[60], 0,0,ierflg);
    
    gMinuit->mnparm(61, "a62", Acalc_Guess[61], step[61], 0,0,ierflg);
    gMinuit->mnparm(62, "a63", Acalc_Guess[62], step[62], 0,0,ierflg);
    gMinuit->mnparm(63, "a64", Acalc_Guess[63], step[63], 0,0,ierflg);
    gMinuit->mnparm(64, "a65", Acalc_Guess[64], step[64], 0,0,ierflg);
    gMinuit->mnparm(65, "a66", Acalc_Guess[65], step[65], 0,0,ierflg);
    gMinuit->mnparm(66, "a67", Acalc_Guess[66], step[66], 0,0,ierflg);
    gMinuit->mnparm(67, "a68", Acalc_Guess[67], step[67], 0,0,ierflg);
    gMinuit->mnparm(68, "a69", Acalc_Guess[68], step[68], 0,0,ierflg);
    gMinuit->mnparm(69, "a70", Acalc_Guess[69], step[69], 0,0,ierflg);
    gMinuit->mnparm(70, "a71", Acalc_Guess[70], step[70], 0,0,ierflg);
    gMinuit->mnparm(71, "a72", Acalc_Guess[71], step[71], 0,0,ierflg);
    gMinuit->mnparm(72, "a73", Acalc_Guess[72], step[72], 0,0,ierflg);
       
    // Now ready for minimization step
    arglist[0] = 500000; // max number of calls to fcn
    arglist[1] = 1e-20; // tolerance
    gMinuit->mnexcm("SIMPLEX", arglist ,2,ierflg);
    
    
    //gMinuit->SetMaxIteration (10000);
    //gMinuit->Migrad ();
    //gMinuit->GetStatus();
    

    // Print results
    cout << "\nPrint results from minuit\n";
    double ParamVal;
    double ParamErr;
    
    
    // Access to these parameters, use:
    Double_t amin,edm,errdef;
    Int_t nvpar,nparx,icstat;
    gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
    //gMinuit->mnprin(3,amin);
    //show print 3
    
    
    
    
    cout << "\n";
    cout << " Minimum chi square = " << amin << "\n";
    cout << " Estimated vert. distance to min. = " << edm << "\n";
    cout << " Number of variable parameters = " << nvpar << "\n";
    cout << " Highest number of parameters defined by user = " << nparx << "\n";
    cout << " Status of covariance matrix = " << icstat << "\n";
    
    cout << "\n";
    gMinuit->mnprin(3,amin);
}