/*

FakePhysicist.com
admin@fakephysicist.com
13 April 2017

Compile with:
g++ v-i_advance.cpp -o v-i_advance `root-config --cflags --libs`

*/

#include <iostream>
#include <iomanip>
#include <cmath>
#include <fstream>
#include <vector>
#include <cstring>

/****** Include the required ROOT Classes *****/
#include "TF1.h"
#include "TH1.h"
#include "TStyle.h"
#include "TPad.h"
#include "TCanvas.h"
#include "TGraphErrors.h"
#include "TPostScript.h"
#include "TLatex.h"
#include "TGaxis.h"
#include "TAxis.h"

using namespace std;

const int vPoints = 11;		// No of voltage reading
const double v_voltage[vPoints] = {
  0.,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0
};

int main() {

  ifstream infile("/media/surya/Surya_1/Online/FakePhysicist/2.v-i_advance/vidata_advance.txt");

  vector<double> v_current[vPoints];
  for(int ij=0;ij<vPoints;ij++) {v_current[ij].clear();}

  while(!infile.eof()) {
    double tm_voltage, tm_current;
    infile >> tm_voltage >> tm_current;
    if(infile.eof()) {break;}
    
    v_current[int(2.*tm_voltage)].push_back(tm_current); // Pushing the current value to corresponding vector array
  }

  /************** Now starts the ROOT part ************/

  gStyle->SetPadBottomMargin(0.14);
  gStyle->SetPadTopMargin(0.09);
  gStyle->SetPadLeftMargin(0.14);
  gStyle->SetPadRightMargin(0.14);

  gStyle->SetOptFit(1111);
  gStyle->SetOptStat(0);  
  TGaxis::SetMaxDigits(3);

  // We will use histograms here.
  // 1-D histograms - TH1F which can be filled with float data.
  // For vPoints of dataset we will create vPoints no of TH1F
  TH1F *currentHisto[vPoints];

  // Now, it is time to set name, title and bin for each histogram
  for(int ij=0;ij<vPoints;ij++) {
    char tempchar[100];
    sprintf(tempchar,"Histogram_%.1fv",v_voltage[ij]);
    currentHisto[ij] = new TH1F(tempchar,tempchar,50,v_current[ij][0] - 0.0005,v_current[ij][0] + 0.0005);
    // Arguments: name, title, bin no, xmin, ymin
    // Since each histogram is for different current values, the ranges for each histogram was set around the value of the first set.

    // Now, we fill the histogram
    for(int jk=0;jk<v_current[ij].size();jk++) {
      currentHisto[ij]->Fill(v_current[ij][jk]);
      // For histogram, each element needs to be filled one at a time. Not like TGraph where only the array is needed.
    }
  }

  double errorCurrent[vPoints], meanCurrent[vPoints];	// To store the error on measurement and the mean position for each voltage value.

  TCanvas *can1 = new TCanvas("viplot1","viplot1",900,1200);
  can1->Divide(3,4);		// Divide the canvas into 12 pads to draw those 12 histograms
  
  TCanvas *can2 = new TCanvas("viplot2","viplot2",600,600);

  TPostScript *pss = new TPostScript("viplot.ps",111);

  pss->NewPage();
  
  for(int ij=0;ij<vPoints;ij++) {
    can1->cd(ij+1);		// Call the pad on which the histogram is printed, pad no starts from 1 (not 0)
    currentHisto[ij]->Draw();
    currentHisto[ij]->Fit("gaus","SQ"); // Fit the histogram with gauss : [0]*exp(-0.5*((x-[1])/[2])**2).
    TF1 *fitfn1 = (TF1*)currentHisto[ij]->GetFunction("gaus");
    meanCurrent[ij] = fitfn1->GetParameter(1);
    errorCurrent[ij] = fitfn1->GetParameter(2);

    cout << "\t" << v_voltage[ij]
	 << "\t" << meanCurrent[ij]
	 << "\t" << errorCurrent[ij]
	 << endl;
    
  }
  can1->Update();

  gStyle->SetOptFit(0);
  gStyle->SetOptStat(0);

  pss->NewPage();

  // Now, we create a TGraph with errors
  TGraphErrors *viplot = new TGraphErrors(vPoints,&meanCurrent[0],&v_voltage[0],&errorCurrent[0],0);
  // Last two Arguments: x-error, y-error (here zero)

  can2->cd();
  viplot->Draw("AP*");
  viplot->SetTitle("Voltage VS Current Plot");
  viplot->GetXaxis()->SetTitle("Current in A");
  viplot->GetXaxis()->SetTitleSize(0.06);
  viplot->GetXaxis()->SetTitleOffset(1.0);
  viplot->GetXaxis()->CenterTitle();
  viplot->GetXaxis()->SetLabelSize(0.06);
  viplot->GetXaxis()->SetLabelOffset(0.001);
  viplot->GetYaxis()->SetTitle("Voltage in V");
  viplot->GetYaxis()->SetTitleSize(0.06);
  viplot->GetYaxis()->SetTitleOffset(0.9);
  viplot->GetYaxis()->CenterTitle();
  viplot->GetYaxis()->SetLabelSize(0.06);
  viplot->GetYaxis()->SetLabelOffset(0.001);
  
  viplot->Fit("pol1","SQ");
  TF1 *fitfn1 = (TF1*)viplot->GetFunction("pol1");
  double resistance = fitfn1->GetParameter(1);
  double reserror = fitfn1->GetParError(1);

  char tempchar[100];
  sprintf(tempchar,"R = (%.0f #pm %.1f) #Omega",resistance,reserror);
  TLatex *Tl = new TLatex(0.4,0.25,tempchar);
  Tl->SetNDC(kTRUE);
  Tl->Draw();

  can2->Update();

  pss->Close();

  return 0;

}