/*

FakePhysicist.com
admin@fakephysicist.com
19 April 2017

Compile with:
g++ GravityFallingBody.cpp -o GravityFallingBody `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 "TGraph.h"
#include "TGraphErrors.h"
#include "TPostScript.h"
#include "TLatex.h"
#include "TGaxis.h"
#include "TAxis.h"

using namespace std;

const int pPoints = 26;		// No of Strip Edges
const double p_Gap = 5.;	// Gap in mm

int main() {

  ifstream infile("/media/surya/Surya_1/Online/FakePhysicist/6.GravityFallingBody/GravityFallingBodyData.txt");

  double p_Pos[pPoints];	// Positions
  for(int ij=0;ij<pPoints;ij++) {
    p_Pos[ij] = ij*p_Gap;
  }
  
  vector<double> p_Time[pPoints]; // Time vector
  for(int ij=0;ij<pPoints;ij++) { p_Time[ij].clear();}

  double startTime;
  while(!infile.eof()) {
    int tm_cnt;
    double tm_time;
    infile >> tm_cnt >> tm_time;
    if(infile.eof()) {break;}
    if(tm_cnt == 0) { startTime = tm_time;}
    p_Time[tm_cnt].push_back((tm_time - startTime)/1.e6); // Pushing the time value to corresponding vector array in Second
  }
  
  
  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);

  TCanvas *can1 = new TCanvas("gravityplot1","gravityplot1",900,1200);
  can1->Divide(4,7);
  TCanvas *can2 = new TCanvas("gravityplot2","gravityplot2",900,1200);
  can2->Divide(3,5);
  TCanvas *can3 = new TCanvas("gravityplot3","gravityplot3",600,600);
  TPostScript *pss = new TPostScript("gravityplot.ps",111);
  

  TH1F *timeHisto[pPoints];	// For time
  for(int ij=0;ij<pPoints;ij++) {
    char tempchar[100];
    sprintf(tempchar,"TimeHisto_%.0fmm",p_Pos[ij]);
    timeHisto[ij] = new TH1F(tempchar,tempchar,50, p_Time[ij][0]-0.05, p_Time[ij][0]+0.05);

    // Filling the histograms
    for(int jk=0;jk<p_Time[ij].size();jk++) {
      timeHisto[ij]->Fill(p_Time[ij][jk]);
    }
  }
  double errorTime[pPoints], meanTime[pPoints];	// To store the error on measurement and the mean for each TIME value.
  pss->NewPage();
  for(int ij=0;ij<pPoints;ij++) {
    can1->cd(ij+1);
    timeHisto[ij]->Draw();
    timeHisto[ij]->Fit("gaus","SQ");
    TF1 *fitfn1 = (TF1*)timeHisto[ij]->GetFunction("gaus");
    meanTime[ij] = fitfn1->GetParameter(1);
    errorTime[ij] = fitfn1->GetParameter(2);    
  }
  can1->Update();


  /******** Getting Error on Distnace is little bit tricky ****/
  TH1F *posHisto[pPoints];	// For position
  for(int ij=0;ij<pPoints;ij++) {
    char tempchar[100];
    sprintf(tempchar,"PosHisto_%.0fmm",p_Pos[ij]);
    posHisto[ij] = new TH1F(tempchar,tempchar, 50, p_Pos[ij]-1.5, p_Pos[ij]+1.5);    
  }

  TF1 *sfunc = new TF1("sfunc","[0] + [1]*x + 0.5*[2]*x*x",0.,150.); // Function: s = ut + (1/2)*g*t^2 : [1] = u, [2] = g
  double tempTime[pPoints], tempPar[3];
  for(int jk=0;jk<p_Time[0].size();jk++) {
    for(int ij=0;ij<pPoints;ij++) {
      tempTime[ij] = p_Time[ij][jk];
    }
    TGraph *g1 = new TGraph(pPoints,&tempTime[0],&p_Pos[0]); // Distance VS Time Plot
    can2->cd(1);
    g1->Fit("sfunc","SQ");
    TF1 *fitfn1 = (TF1*)g1->GetFunction("sfunc");
    for(int ki=0;ki<3;ki++) {
      tempPar[ki] = fitfn1->GetParameter(ki);
    }
    // Filling the expected position. This is because of the manual cutiing of the PVC tape
    for(int ij=0;ij<pPoints;ij++) {
      posHisto[ij]->Fill(tempPar[0] + tempPar[1]*tempTime[ij] + 0.5*tempPar[2]*tempTime[ij]*tempTime[ij]);
    }
  }
  double errorPos[pPoints], meanPos[pPoints];	// To store the error and the mean for each POSITION value.
  pss->NewPage();
  for(int ij=0;ij<pPoints;ij++) {
    can1->cd(ij+1);
    posHisto[ij]->Draw();
    posHisto[ij]->Fit("gaus","SQ");
    TF1 *fitfn1 = (TF1*)posHisto[ij]->GetFunction("gaus");
    meanPos[ij] = fitfn1->GetParameter(1);
    errorPos[ij] = fitfn1->GetParameter(2);    
  }
  can1->Update();

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

  
  // Getting the Statistical error
  TGraph *g1 = new TGraph(pPoints,&meanTime[0],&meanPos[0]); // Distance VS Time Plot
  can2->cd(1);
  g1->Fit("sfunc","SQ");
  TF1 *fitfn1 = (TF1*)g1->GetFunction("sfunc");
  double gErrorStat = fitfn1->GetParError(2);

  // Getting g-value and systemic error
  pss->NewPage();
  TGraphErrors *gravityplot = new TGraphErrors(pPoints,&meanTime[0],&meanPos[0],&errorTime[0],&errorPos[0]);
  can3->cd();
  gravityplot->Draw("AP*");
  gravityplot->SetTitle("Distance VS Time Plot");
  gravityplot->GetXaxis()->SetTitle("Time in Sec");
  gravityplot->GetXaxis()->SetTitleSize(0.025);
  gravityplot->GetXaxis()->SetTitleOffset(1.0);
  gravityplot->GetXaxis()->CenterTitle();
  gravityplot->GetXaxis()->SetLabelSize(0.025);
  gravityplot->GetXaxis()->SetLabelOffset(0.001);
  gravityplot->GetYaxis()->SetTitle("Distance in mm");
  gravityplot->GetYaxis()->SetTitleSize(0.025);
  gravityplot->GetYaxis()->SetTitleOffset(1.1);
  gravityplot->GetYaxis()->CenterTitle();
  gravityplot->GetYaxis()->SetLabelSize(0.025);
  gravityplot->GetYaxis()->SetLabelOffset(0.001);
  gravityplot->Fit("sfunc","SQ");
  TF1 *fitfn2 = (TF1*)gravityplot->GetFunction("sfunc");
  double gValue = fitfn2->GetParameter(2);
  double gErrorSys = fitfn2->GetParError(2);
  char tempchar[100];
  sprintf(tempchar,"g = (%.0f #pm %.1f (Stat) #pm %.1f (Sys)) mm-s^{-2}",gValue,gErrorStat,gErrorSys);
  TLatex *Tl = new TLatex(0.2,0.2,tempchar);
  Tl->SetTextSize(0.025);
  Tl->SetNDC(kTRUE);
  Tl->Draw();
  can3->Update();

  pss->Close();

  return 0;

}