#ifndef NeutronMulti_H
#define NeutronMulti_H
#include<vector>
#include<iostream>
#include<cmath>
#include<cassert>
#if !defined(__CINT__) || defined(G__DICTIONARY)
#include "Minuit2/FCNBase.h"
#endif
#include "NeutronMulti.h"
#include <fstream>
#include <iostream>
using namespace ROOT::Minuit2;
class NeutronMultiFcn : public FCNBase{


public:
NeutronMultiFcn(const std::vector<double>& en, 
	    const std::vector<double>& ye,
		const std::vector<double>& er_ye): energy(en), 
    yield(ye), erroryield(er_ye){
    //cout<<en[0]<<"   "<<en[1]<<endl;
  }
        

~NeutronMultiFcn() {}

virtual double Up() const{return theerror;}
virtual double operator()(const std::vector<double>&) const;

std::vector<double> measurements() const{return energy;}
std::vector<double> theyield() const{return yield;}
std::vector<double> error() const{return erroryield;}

private:
std::vector<double> energy;
std::vector<double> yield;
std::vector<double> erroryield;
double theerror;

};
#endif


double NeutronMultiFcn::operator()(const std::vector<double>& par) const{
  
  assert(par.size()==5);
  cout<<par[0]<<"   "<<par[1]<<"   "<<par[2]<<"   "<<par[3]<<"   "<<par[4]<<endl;
NeutronMulti nm(par[0],par[1],par[2],par[3],par[4]);
 //NeutronMulti nm(2,2.1,2.2,2,2);
  double chi2=0;
  for(unsigned int i=0; i!=yield.size(); i++){
    double num = (nm(energy[i])-yield[i]);
    if(i ==0){
       std::cout<<"Theory: "<<num<<std::endl;
       std::cout<<"Esti: "<<yield[i]<<std::endl;
       std::cout<<"EstiE: "<<erroryield[i]<<std::endl;
       }
    double den = erroryield[i];
    double del = num/den;
     if(del==del)
    chi2 = chi2+del*del;
     if(i ==0) 
       {    
       cout<<"chi2: "<<chi2<<endl; }
  }
    cout<<"chi2: "<<chi2<<endl;
  return chi2;
  
}
/*  int main(){ */
  
/*      std::ifstream inpfile("F7.dat"); */
/*      // TH1D* h1 = new TH1D("h1","h2",100,0,20); */
/*   std::vector<double> the_E; */
/*   std::vector<double> the_Y; */
/*   std::vector<double> the_Yer; */
/*   int i=0; */

/*   while(inpfile.good()) */
/*     { i++; */
/*       double en =0; */
/*       double yi =0; */
/*       double yerr =0; */
/*       inpfile >> en >> yi >> yerr ; */
/*       the_E.push_back(en); */
/*       the_Y.push_back(yi); */
/*       the_Yer.push_back(yerr); */
/*       cout<<i<<endl; */
/*     } */
/*   //Double yxxxx =  nm(2,2.1,2.1,2,2); */
/*   } */