//Content-Disposition: inline; filename="Vitt.C"

//
// Example of a fitting program
// ============================
//
// The fitting function fcn is a simple chisquare function
// The data consists of 5 data points (arrays x,y,z) + the errors in errorsz
// More details on the various functions or parameters for these functions
// can be obtained in an interactive ROOT session with:
// Root > TMinuit *minuit = new TMinuit(11);
// Root > minuit->mnhelp("*",0) to see the list of possible keywords
// Root > minuit->mnhelp("SET",0) explains most parameters
//
#include "TROOT.h"
#include "TMinuit.h"
#include "TRandom.h"
#include "TMath.h"
#include <iostream>
#include <fstream>
#include <cstdlib>
#include <cmath>
#include <string>
#include <vector>
#include <TMinuit.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TROOT.h>
#include <TF1.h>






Double_t bincent[300],content[300],err[300];
extern Double_t func(Double_t *Ep,Double_t *par);
extern void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
Int_t ncount = 0;


//TRandom r;

//______________________________________________________________________________
minimise()
{  
      TFile *f1 = new TFile("spectrum.root");
   f1->ls();
   TH1F *hmass = (TH1F*)f1->Get("hx2");
   Int_t size = hx2->GetXaxis()->GetNbins();
   cout <<size<<endl;
//    Double_t * bincent = new Double_t [size];
//     Double_t * content = new Double_t [size];
//      Double_t * err = new Double_t [size];
  
   for (int i = 1; i <=  size; i++){
      bincent[i] = hmass->GetBinCenter(i);
      content[i] = hmass->GetBinContent(i);
      err[i]= hmass->GetBinError(i);
//    cout <<bincent[i]<<"\t"<<content[i]<<"\t"<<err[i]<<endl;
     }



//TROOT minexam("TestMinuit","Test of Minuit");
const Int_t npars = 3;
TMinuit *gMinuit = new TMinuit(npars); //initialize TMinuit with a maximum of 5 params
gMinuit->SetFCN(fcn);

Double_t arglist[300];
//  Double_t * arglist = new Double_t [size];
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 vstart[3] = {220.,122.,28.6};
//static
 Double_t step[3] = {0.1,0.1,0.1};
gMinuit->mnparm(0, "Constant", vstart[0], step[0], 0,0,ierflg);
gMinuit->mnparm(1, "Mean", vstart[1], step[1], 0,0,ierflg);
gMinuit->mnparm(2, "Sigma", vstart[2], step[2], 0,0,ierflg);



// Now ready for minimization step
arglist[0] = 500;
arglist[1] = 1;//1
gMinuit->mnexcm("MIGRAD", arglist ,2,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);

//_____________________________________________________________________
   Int_t iuext;
   TString chnam;   // The name of the parameter
   Double_t val;    // The current (external) value of the parameter 
   Double_t errl;   // The current estimate of the parameter uncertainty  
   Double_t xlolim; // The lower bound (or zero if no limits)
   Double_t xuplim; // The upper bound (or zero if no limits)
   Int_t iuint;     // The internal parameter number 
   
   
   Double_t currentPar[npars] = {0};
   for (Int_t i=0; i< npars;i++) {
      gMinuit->mnpout(i, chnam, currentPar[i], errl, xlolim, xuplim, iuint);
   }


     Int_t gcount = 300;
   TGraphErrors *gr = new TGraphErrors(gcount,bincent,content,0,err);
   gr->Draw("AP");
//   gr->SetMarkerStyle(4);
//   gr->SetMinimum(0.0001);
//   gr->SetMaximum(0.5);


//     hmass->Draw(); 
     TF1 *fun_1=new TF1("fun_1",func,55,195,3); 	
     fun_1->SetParameters(currentPar);
     fun_1->SetLineColor(kBlue);
     fun_1->SetLineStyle(1);
     fun_1->SetLineWidth(2);
     fun_1->Draw("same");
}
//______________________________________________________________________________
void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
 Int_t nbins;
 nbins = 300;

//calculate chisquare
Double_t chisq = 0;
Double_t delta;
for (Int_t i=0;i<nbins; i++) {
delta = (content[i]-func(&bincent[i],par))/err[i];

chisq += delta*delta;
}
f = chisq;
ncount++;
}


Double_t func(Double_t *Ep,Double_t *par)
{
//    MOVING SOURCE FITTING FUNCTION !!NOW O + Au
  Double_t arg=0;
  Double_t E=*Ep;

//____________________________________________________________________

   if(par[1] !=0 && par[2] !=0)
     {
     
     arg = (E - par[1])/par[2];


Double_t fitval = 
par[0]*TMath::Exp(-0.5*((E-par[1])/par[2])*((E-par[1])/par[2]));}

//(par[0]*TMath::Exp(-0.5*arg*arg))/(sqrt(2.*TMath::Pi()*par[1]*par[1])); }



     return fitval;
}