/*

FakePhysicist.com
admin@fakephysicist.com
18 April 2017

Compile with:
g++ rootread.cpp testTree.C -o rootread `root-config --cflags --libs`

*/

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

/****** Include the required ROOT Classes *****/
#include "TFile.h"
#include "TTree.h"
#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"

#include "testTree.h"		// Structure of current root file

using namespace std;

int main() {

  TFile *fileIn = new TFile("/media/surya/Surya_1/Online/FakePhysicist/5.ROOTFiles/test.root","READ");
  TTree *treeIn = (TTree*)fileIn->Get("testTree"); // creating a pointer of the current TTree
  testTree *event = new testTree(treeIn);	   // creating a pointer of testTree to read data
  event->Loop();				   // Looping over all data entries

  int nentries = int(treeIn->GetEntries()); // Getting number of entries in the tree

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

  for(int ij=0;ij<nentries;ij++) {

    fileIn->cd();		// Calling the current file
    event->GetEntry(ij);	// Prepare the entry for reading
    
    for(int jk=0;jk<vPoints;jk++) {
      if(ij == 0) {
	v_voltage[jk] = event->reading[jk]; // Reading voltage
      } else {
	v_current[int(2.*v_voltage[jk])].push_back(event->reading[jk]); // Reading current
      }
    }
  }
  

  /************** Same as V-I Advanced Plotting ************/

  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);

  TH1F *currentHisto[vPoints];

  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);

    for(int jk=0;jk<v_current[ij].size();jk++) {
      currentHisto[ij]->Fill(v_current[ij][jk]);
    }
  }

  double errorCurrent[vPoints], meanCurrent[vPoints];

  TCanvas *can1 = new TCanvas("viplot1","viplot1",900,1200);
  can1->Divide(3,4);
  
  TCanvas *can2 = new TCanvas("viplot2","viplot2",600,600);

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

  pss->NewPage();
  
  for(int ij=0;ij<vPoints;ij++) {
    can1->cd(ij+1);
    currentHisto[ij]->Draw();
    currentHisto[ij]->Fit("gaus","SQ");
    TF1 *fitfn1 = (TF1*)currentHisto[ij]->GetFunction("gaus");
    meanCurrent[ij] = fitfn1->GetParameter(1);
    errorCurrent[ij] = fitfn1->GetParameter(2);    
  }
  can1->Update();

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

  pss->NewPage();

  TGraphErrors *viplot = new TGraphErrors(vPoints,&meanCurrent[0],&v_voltage[0],&errorCurrent[0],0);

  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;

}