#include "Minuit2/Minuit2Minimizer.h"
#include "Math/Functor.h"
#include "TFile.h"
#include "TTree.h"
#include "TROOT.h"
#include <iostream>
using namespace std;

TFile *tf = NULL;
TTree *gentof = NULL;

   // define structure

   struct Events_t {

      // generated values
      Float_t dE_V;
      Float_t dE_H;
      Float_t t0_H;
      Float_t t0_V;
      Float_t x;
      Float_t y;
      Float_t tw_R;
      Float_t tw_L;
      Float_t tw_T;
      Float_t tw_B;
      Float_t t_offset;
      Float_t sigma_t;
      Float_t veff;
      Float_t lambda;
      Float_t Vmax;
      Float_t Vthr;

      // reconstructed values
      Float_t P_R;
      Float_t P_L;
      Float_t P_T;
      Float_t P_B;
      Float_t t_R;
      Float_t t_L;
      Float_t t_T;
      Float_t t_B;
      Float_t tcor_R;
      Float_t tcor_L;
      Float_t tcor_T;
      Float_t tcor_B;
      };
   Events_t Events;

double tof_chi2(const double *xx ){

   // constain constants for all pmts to be the same

   Double_t walk0L = xx[0];
   Double_t walk0R = xx[0];
   Double_t walk0B = xx[0];
   Double_t walk0T = xx[0];
   Double_t walk1L = xx[1];
   Double_t walk1R = xx[1];
   Double_t walk1B = xx[1];
   Double_t walk1T = xx[1];
   Double_t chi2=0;

   /*char string[256];
   Int_t jj;
   Double_t pi=3.14159;*/
   Int_t nevents = 500;
 
   /*char filename[132];
   sprintf(filename,"tof_tree.root");
   printf ("chi2: Input filename = %s\n",filename);
   if (!tf) TFile *tf = new TFile(filename);
   gentof = (TTree*) tf->Get("gentof");
   gentof->SetBranchAddress("Events",&Events.dE_V);
   gentof->Print();*/
  
   Int_t nentries = gentof->GetEntries();
   for (Int_t j=0; j<nentries && j<nevents; j++){
      gentof->GetEntry(j);

      Double_t tLc = Events.t_L - walk0L/(1+walk1L*sqrt(Events.P_L));
      Double_t tRc = Events.t_R - walk0R/(1+walk1R*sqrt(Events.P_R));
      Double_t tBc = Events.t_B - walk0B/(1+walk1B*sqrt(Events.P_B));
      Double_t tTc = Events.t_T - walk0T/(1+walk1T*sqrt(Events.P_T));

      Double_t t0H = 0.5*(tLc+tRc);
      Double_t t0V = 0.5*(tBc+tTc);
      // printf ("t0H =%f, t0V=%f\n",t0H,t0V);

     chi2 = chi2  + (t0H-t0V)*(t0H-t0V)/(Events. sigma_t*Events.sigma_t);
   }

  // printf ("chi2 = %f walk0L=%f, walk0R=%f, walk0B=%f, walk0T=%f \n\n",chi2,walk0L,walk0R,walk0B,walk0T);
  return chi2;
}

int tof_calibration_fit2()
{

   // read root tree 
   char filename[132];

   sprintf(filename,"tof_tree.root");
   printf ("tof_calibration_fit: Input filename = %s\n",filename);
   tf = new TFile(filename);
   gentof = (TTree*) tf->Get("gentof");
   gentof->SetBranchAddress("Events",&Events.dE_V);
   gentof->Print();

   // Choose method upon creation between:
   // kMigrad, kSimplex, kCombined, 
   // kScan, kFumili
   ROOT::Minuit2::Minuit2Minimizer min ( ROOT::Minuit2::kMigrad );

   min.SetMaxFunctionCalls(1000000);
   min.SetMaxIterations(100000);
   min.SetTolerance(0.001);

   // use same constants for all PMTs

   ROOT::Math::Functor f(&tof_chi2,2); 
   double step[2] = {0.01,0.01};
   double variable[2] = { 10,0.8};

   min.SetFunction(f);

   // Set the free variables to be minimized!
   min.SetVariable(0,"walk0",variable[0], step[0]);
   min.SetVariable(1,"walk1",variable[1], step[1]);

   min.Minimize(); 

   const double *xs = min.X();
   const double *xerr = min.Errors();
   cout << "Minimum: Chi2(" << xs[0] << "," << xs[0] << "," << xs[0] << "," << xs[0] << "," << xs[1] << "," << xs[1] << "," << xs[1] << "," << xs[1] << "): " 
       << tof_chi2(xs) << endl;

   cout << "Minimum: Errors(" << xerr[0] << "," << xerr[0] << "," << xerr[0] << "," << xerr[0] << "," << xerr[1] << "," << xerr[1] << "," << xerr[1] << "," << xerr[1] << "): "  << endl;


  delete tf;

   return 0;
}