// A simple C++ program to illustrate the use of ROOT class TMinuit 
// for function minimization.  The example shows a Maximum Likelihood
// fit in which TMinuit minimizes - 2*log(L).   
// fcn passes back f = -2*ln L by reference; this is the function to minimize.
// The factor of -2 allows MINUIT to get the errors using the same
// recipe as for least squares, i.e., go up from the minimum by 1.

// TMinuit does not allow fcn to be a member function, and the function
// arguments are fixed, so the one of the only ways to bring the data  
// into fcn is to declare a pointer to the data (xVecPtr) as global.

// For more info on TMinuit see root.cern.ch/root/html/TMinuit.html .

// Glen Cowan
// RHUL Physics
// 9 August, 2011

#include <iostream>
#include <fstream>
#include <sstream>
#include <cstdlib>
#include <cmath>
#include <string>
#include <vector>

#include <TMinuit.h>
#include <TApplication.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TF1.h>
#include <TH1F.h>
#include <TAxis.h>
#include <TLine.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TROOT.h>

using namespace std;

// Declare pointers to data as global (not elegant but TMinuit needs this).

vector<double>* Rdata = new vector<double>();
vector<double>* VLongCables = new vector<double>();
vector<double>* VShortCables = new vector<double>();
vector<double>* dVLC = new vector<double>();
vector<double>* dVSC = new vector<double>();

void plotFit(TCanvas* canvas, string outfile, double xMin, double xMax,
     double yMin, double yMax, double par[], 
     vector<double>* xVec, vector<double>* yVec, vector<double>* dyVec);

double BandpassIntegral(double R, double C){

  const double pi = 3.14159265;

  static vector<double> fVec;
  static vector<double> g2Vec;
  static bool firstCall = true; 
  if (firstCall) {

    firstCall = false;
    string fileName = "./FilterData.txt";
    ifstream inFile;
    inFile.open(fileName.c_str());
    if (inFile.fail()) { 
      cout << "Couldn't open filter data file!" << endl;
      exit(1); 
    }

    istringstream instream;       // Declare an input string stream
    string line;
    int lineNum = 0;
    while ( getline(inFile, line) ) {
      bool useLine = !( line.substr(0,1) == "#" || line.substr(0,1) == "!");
      if ( useLine ) {
        instream.clear();          // Reset from possible previous errors
        instream.str(line);        // Use line as source of input
        double f, v0, vi, g2, b;
        lineNum++;
        instream >> f >> v0 >> vi >> g2 >> b;
        fVec.push_back(f);
        g2Vec.push_back(g2);
      }
    }

  }            // end of initialization

  double bpi = 0.;
  for (int i=0; i<fVec.size()-1; i++){
    double fAve = 0.5 * (fVec[i] + fVec[i+1]);
    double alpha = 2.*pi*fAve*R*C;
    bpi += 0.5*(g2Vec[i] + g2Vec[i+1])*(fVec[i+1] - fVec[i]) /
           (1. + alpha*alpha);
  }

  return bpi;

}

// The function to be fitted.
// First argument needs to be a pointer in order to plot with the TF1 class.

double Vrms(double* xArray, double par[]){        
  double R = xArray[0];
  double kB = par[0];     // Boltzmann's const.
  double C  = par[1];
  double T  = par[2];     // temperature
  double B = BandpassIntegral(R,C);
  double val = sqrt(4*R*kB*T*B);
  return val;
}

//-------------------------------------------------------------------------

//  function to read in the data from a file

void getData(vector<double>* Rdata, 
             vector<double>* VLongCables, 
             vector<double>* VShortCables){

  //  string fileName;
  //  cout << "Enter file name: ";
  //  cin >> fileName;
  string fileName = "VrmsData.txt";

  ifstream inFile;
  inFile.open(fileName.c_str());
  if (inFile.fail()) { 
    cout << "Couldn't open file!" << endl;
    exit(1); 
  }

  // KLUDGE -- for now set 2% relative errors by hand

  double relErrLC = 0.02;
  double relErrSC = 0.02;

  istringstream instream;       // Declare an input string stream
  string line;
  int lineNum = 0;
  while ( getline(inFile, line) ) {
    bool useLine = !( line.substr(0,1) == "#" || line.substr(0,1) == "!");
    if ( useLine ) {
      instream.clear();          // Reset from possible previous errors
      instream.str(line);        // Use line as source of input
      double r, vlc, vsc;
      lineNum++;
      instream >> r >> vlc >> vsc;
      Rdata->push_back(r);
      VLongCables->push_back(vlc*1.e-6);     // convert nV to V
      VShortCables->push_back(vsc*1.e-6);    // convert nV to V
      dVLC->push_back(relErrLC*vlc*1.e-6);
      dVSC->push_back(relErrSC*vsc*1.e-6);
    }
  }
  int numPoints = Rdata->size();

  cout << "Read info for " << numPoints << "    data points" << endl;
  for (int i=0; i<numPoints; i++){
    cout << (*Rdata)[i] << "  " << (*VLongCables)[i] << "  " 
         << (*VShortCables)[i] << endl;
  }

}

//-------------------------------------------------------------------------

// fcn passes back f = - 2*ln(L), the function to be minimized.

void fcn(int& npar, double* deriv, double& f, double par[], int flag){

  double parLC[3];
  parLC[0] = par[0];       // kB
  parLC[1] = par[1];       // CL
  parLC[2] = par[3];       // T

  double parSC[3];
  parSC[0] = par[0];       // kB
  parSC[1] = par[2];       // CS
  parSC[2] = par[3];       // T

  int n = (*Rdata).size();
  double chi2 = 0.0;
  for (int i=0; i<n; i++){
    double R = (*Rdata)[i];
    double VLC = (*VLongCables)[i];
    double VSC = (*VShortCables)[i];
    double VpredLC = Vrms(&R, parLC);
    double VpredSC = Vrms(&R, parSC);
    double deltaLC = (VLC - VpredLC)/(*dVLC)[i];
    double deltaSC = (VSC - VpredSC)/(*dVSC)[i];
    chi2 += deltaLC*deltaLC  +  deltaSC*deltaSC;   
  }
  f = chi2;
 
}                         // end of fcn

//-------------------------------------------------------------------------

int main(int argc, char **argv) {

  TApplication theApp("App", &argc, argv);
  TCanvas* canvas = new TCanvas("canvas", "canvas",  10, 10, 500, 500);

// Set a bunch of parameters to make the plot look nice

  canvas->SetFillColor(0);
  canvas->UseCurrentStyle();
  canvas->SetBorderMode(0);       // still leaves red frame bottom and right
  canvas->SetFrameBorderMode(0);   // need this to turn off red hist frame!
  gROOT->SetStyle("Plain");
  canvas->UseCurrentStyle();
  gROOT->ForceStyle();

  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(0);
  gStyle->SetTitleBorderSize(0);
  gStyle->SetTitleSize(0.04);
  gStyle->SetTitleFont(42, "hxy");      // for histogram and axis titles
  gStyle->SetLabelFont(42, "xyz");      // for axis labels (values)
  gROOT->ForceStyle();

// Read in the data (declared global)

  getData(Rdata, VLongCables, VShortCables);   

// Initialize minuit, set initial values etc. of parameters.

  const int npar = 4;              // the number of parameters
  TMinuit minuit(npar);
  minuit.SetFCN(fcn);

  double par[npar];               // the start values
  double stepSize[npar];          // step sizes 
  double minVal[npar];            // minimum bound on parameter 
  double maxVal[npar];            // maximum bound on parameter
  string parName[npar];

  par[0] = 1e-23;            // a guess
  stepSize[0] = 1e-24; 
  minVal[0] = 0.0;
  maxVal[0] = 0.0;
  parName[0] = "kB";

  par[1] = 72.8e-12;            // Farads
  stepSize[1] = 3.e-12;        // 
  minVal[1] = 0.0;
  maxVal[1] = 0.0;
  parName[1] = "CL";

  par[2] = 49.2e-12;            // Farads
  stepSize[2] = 3.e-12;        // 
  minVal[2] = 0.0;
  maxVal[2] = 0.0;
  parName[2] = "CS";

  par[3] = 300.;            // Kelvin
  stepSize[3] = 0.;        // 
  minVal[3] = 0.0;
  maxVal[3] = 0.0;
  parName[3] = "T";

  for (int i=0; i<npar; i++){
    minuit.DefineParameter(i, parName[i].c_str(), 
      par[i], stepSize[i], minVal[i], maxVal[i]);
  }

  // minuit.FixParameter(1);
  // minuit.FixParameter(2);

// Do the minimization!

  minuit.Migrad();       // Minuit's best minimization algorithm
  double outpar[npar], err[npar];
  for (int i=0; i<npar; i++){
    minuit.GetParameter(i, outpar[i], err[i]);
  }

  cout << "covariance and correlation coefficients..." << endl;
  double covmat[npar][npar];
  minuit.mnemat(&covmat[0][0], npar);
  for (int i=0; i<npar; i++){
    for (int j=0; j<npar; j++){
      double sigma_i = sqrt(covmat[i][i]);
      double sigma_j = sqrt(covmat[j][j]);
      double rho = covmat[i][j]/(sigma_i*sigma_j);
      cout << i << "  " << j << "  " << covmat[i][j] << "  " << rho << endl;
    }
  }

// Plot the result

  double xMin = 0.0;
  double xMax = 1e6;        // Ohms
  double yMin = 0.0;
  double yMax = 1.e-5;      // Volts

  double outparLC[3];
  outparLC[0] = outpar[0];
  outparLC[1] = outpar[1];
  outparLC[2] = outpar[3];

  double outparSC[3];
  outparSC[0] = outpar[0];
  outparSC[1] = outpar[2];
  outparSC[2] = outpar[3];

  string LCplot = "lcplot";
  string SCplot = "scplot";

  plotFit(canvas, LCplot, xMin, xMax, yMin, yMax, outparLC, 
          Rdata, VLongCables, dVLC);

  plotFit(canvas, SCplot, xMin, xMax, yMin, yMax, outparSC, 
          Rdata, VShortCables, dVSC);

  cout << "To exit, quit ROOT from the File menu of the plot" << endl;
  theApp.Run(true);

  canvas->Close();
  delete canvas;
  return 0;

}

void plotFit(TCanvas* canvas, string outfile, double xMin, double xMax,
     double yMin, double yMax, double par[], 
     vector<double>* xVec, vector<double>* yVec, vector<double>* dyVec){

  gPad->SetLeftMargin(0.15);
  gPad->SetRightMargin(0.15);
  gPad->SetTopMargin(0.1);
  gPad->SetBottomMargin(0.2);

  TF1* func = new TF1("func", Vrms, xMin, xMax, 3);
  func->SetLineStyle(1);             //  1 = solid, 2 = dashed, 3 = dotted
  func->SetLineColor(1);             //  black (default)
  func->SetLineWidth(1);

  func->SetParameters(par);
  func->Draw();

  TH1* axhist = func->GetHistogram(); 
  axhist->SetTitle("");
  axhist->SetXTitle("R (Ohms)");
  axhist->SetYTitle("V_{rms} (Volts)");
  axhist->SetMinimum(yMin); 
  axhist->SetMaximum(yMax);
  axhist->SetAxisRange(xMin, xMax);

  TAxis* xa = axhist->GetXaxis();
  TAxis* ya = axhist->GetYaxis();

  xa->SetTitleOffset(1.4);    //  factor multiplies default offset
  ya->SetTitleOffset(1.4);

  xa->SetLabelOffset(0.015);
  ya->SetLabelOffset(0.015);

  xa->SetTickLength(0.015);  // default  = 0.03
  ya->SetTickLength(0.015);  // default  = 0.03

  xa->SetTitleSize(0.04);
  ya->SetTitleSize(0.04);

  //  gPad->SetLogx(1);
  //  xa->SetLimits(90., 700.);

  xa->SetNdivisions(-5); // negative value should force number of divisions?
  ya->SetNdivisions(-5);

  xa->SetLabelSize(0.04);
  ya->SetLabelSize(0.04);

// add the data points

  int numPoints = (*xVec).size();
  double xArray[numPoints];
  double yArray[numPoints];
  double* dxArray = 0;
  double dyArray[numPoints];
  for (int i=0; i<numPoints; i++){
    xArray[i] = (*xVec)[i];
    yArray[i] = (*yVec)[i];
    dyArray[i] = (*dyVec)[i];
    
  }
  TGraphErrors* tg = new TGraphErrors(numPoints, xArray, yArray,
                                dxArray, dyArray);

  tg->SetMarkerColor(1);
  tg->SetMarkerSize(0.8);
  tg->SetMarkerStyle(20);
  tg->Draw("P,same");
  canvas->Show();
  string psfile = outfile + ".eps";
  string giffile = outfile + ".gif";
  canvas->Print(psfile.c_str(), "eps");
  canvas->Print(giffile.c_str(), "gif");

  string command;
  cout << "Enter any key to continue: ";
  getline(cin, command);

}