#include <stdio.h>
#include <iostream>
#include <fstream>
#include <string>
#include <TLegend.h>
#include <TMath.h>


using namespace std;

#include <TCanvas.h>
#include <TROOT.h>
#include <TH1.h>
#include <TF1.h>
#include <TFile.h>
#include <TTree.h>
#include <TChain.h>
#include <TNtuple.h>
#include <TLegend.h>
#include <TGraphErrors.h>



int fact(int n){

     return (n==0) || (n==1) ? 1 : n* fact(n-1);
}


double alpha_shab(double* Energy, double *p) {
	
  double E = Energy[0] ;


 double I_k_eq_1 = p[0];  // Intensity for k =1  
 double E_k_eq_1 = 37640.8*1E-4;  // zero-phonon line for k =1  
 double omega1_k_eq_1 = 13.9541*1E-4;  //mode frequency for k = 1
 double S_k_eq_1 = p[1]  ;     //Huang Rhees parameter for k=1  
 double e_sigma_k_eq_1 = p[2]*1E-4;  // electronic sigma for k = 1  
 double p_sigma_k_eq_1 = p[3]*1E-4;  // phonon sigma for k = 2
  double offset = 0.170281461750221;  // phonon sigma for k = 2
  
  


  return offset + 
   I_k_eq_1*((E_k_eq_1 + 0*omega1_k_eq_1)/(E_k_eq_1))*exp(-S_k_eq_1)*(pow(S_k_eq_1,0)/fact(0))*exp(-pow((E_k_eq_1+0*omega1_k_eq_1-E),2)/(2*(e_sigma_k_eq_1*e_sigma_k_eq_1+0*p_sigma_k_eq_1*p_sigma_k_eq_1)))/sqrt(2*(e_sigma_k_eq_1*e_sigma_k_eq_1+0*p_sigma_k_eq_1*p_sigma_k_eq_1))
  +I_k_eq_1*((E_k_eq_1 + 1*omega1_k_eq_1)/(E_k_eq_1))*exp(-S_k_eq_1)*(pow(S_k_eq_1,1)/fact(1))*exp(-pow((E_k_eq_1+1*omega1_k_eq_1-E),2)/(2*(e_sigma_k_eq_1*e_sigma_k_eq_1+1*p_sigma_k_eq_1*p_sigma_k_eq_1)))/sqrt(2*(e_sigma_k_eq_1*e_sigma_k_eq_1+1*p_sigma_k_eq_1*p_sigma_k_eq_1))
  +I_k_eq_1*((E_k_eq_1 + 2*omega1_k_eq_1)/(E_k_eq_1))*exp(-S_k_eq_1)*(pow(S_k_eq_1,2)/fact(2))*exp(-pow((E_k_eq_1+2*omega1_k_eq_1-E),2)/(2*(e_sigma_k_eq_1*e_sigma_k_eq_1+2*p_sigma_k_eq_1*p_sigma_k_eq_1)))/sqrt(2*(e_sigma_k_eq_1*e_sigma_k_eq_1+2*p_sigma_k_eq_1*p_sigma_k_eq_1))
  +I_k_eq_1*((E_k_eq_1 + 3*omega1_k_eq_1)/(E_k_eq_1))*exp(-S_k_eq_1)*(pow(S_k_eq_1,3)/fact(3))*exp(-pow((E_k_eq_1+3*omega1_k_eq_1-E),2)/(2*(e_sigma_k_eq_1*e_sigma_k_eq_1+3*p_sigma_k_eq_1*p_sigma_k_eq_1)))/sqrt(2*(e_sigma_k_eq_1*e_sigma_k_eq_1+3*p_sigma_k_eq_1*p_sigma_k_eq_1)) ;
}

 
int main(){
	
	 gROOT->cd(); 
  // gROOT->SetStyle("Plain");   //set plain TStyle
   gStyle->SetOptStat(0); //draw statistics on plots, 0 for no output
 // gStyle->SetOptStat(111111); //draw statistics on plots, 0 for no output
   gStyle->SetOptFit(1);    //draw fit results on plot, 0 for no output
   gStyle->SetPalette(57);     //set color map
   gStyle->SetOptTitle(0);    //suppress title box


  const int data_points = 375 ;               // no. of data points
  double x_data[data_points] = {0};	           // data array x_data
  double y_od_data[data_points] = {0};               // data array optical density
  double y_abs_data[data_points] = {0};               // data array absorption intensity
  double y_abs_data_off[data_points] = {0};               // data array absorption intensity subtract offset
  double err_x[data_points] = {.000};               // data array absorption intensity subtract offset
  double err_y[data_points] = {0.000};               // data array absorption intensity subtract offset
  	
  FILE *fp;
  char *fname = "C:/root_v5.34.36/perfil/fitting.txt"; 
 fp = fopen( fname, "r" );
  gROOT->cd(); 
  int ret;

  char buf[3][14];
  double data[4];
  Int_t nlines;
  char c;
 
  
  if( fp == NULL ){
    printf( "%sファイルが開けません\n", fname );
    return -1;
  }

  printf("\n");

for (Int_t nlines=0; nlines < 0; ++nlines){             // number of lines to skip

  do {
    c = fgetc(fp);
    } while (c !=  '\n');

}


 
  while( (ret=fscanf(fp, "%lf,%lf", &data[0], &data[1])) != EOF){
    printf("%le %le\n", data[0], data[1]);
    
   x_data[nlines] = data[0];	           // data array x_data
   y_od_data[nlines] = data[1];               // 
   y_abs_data[nlines] = 1-pow(10,-data[1]);                  // absorption data   
   y_abs_data_off[nlines] = 1-pow(10,-data[1])-y_abs_data[0];    // absorption data- offset
   
   
   
 ++nlines;
  }
  
  cout << y_abs_data[0] << endl;


  printf("\n");
  fclose( fp );
  
  
  
    TCanvas *c1 = new TCanvas("c1","Good Vibrations", 700, 500);
    graph1 = new TGraphErrors(data_points, x_data, y_abs_data,  err_y, err_x);
     graph1->GetXaxis()->SetRangeUser(3.76255, 3.78035);
	 graph1->GetYaxis()->SetRangeUser(y_abs_data[0], 1);
	 graph1->GetXaxis()->SetTitle("Wave Number (10^{4} cm^{-1})");
	 graph1->GetXaxis()->SetTitleOffset(0.8);
     graph1->GetYaxis()->SetTitle("Absorption Intensity");
     graph1->GetYaxis()->SetTitleOffset(0.55);
     graph1->GetYaxis()->SetDecimals(3); 
     graph1->GetXaxis()->SetDecimals(2); 
     graph1->GetXaxis()->CenterTitle();
     graph1->GetYaxis()->CenterTitle();
	 graph1->SetMarkerSize(0);
	 graph1->SetMarkerColor(kRed);	 
	 graph1->SetMarkerStyle(20);
     graph1->SetLineColor(kRed);
     graph1->SetLineWidth(2);
     graph1->Draw();
     graph1->SetName("Good Vibrations");
	
	
 double I_k_eq_1 = 0.00229662;  // Intensity for k =1  
 double S_k_eq_1 =    3.3496 ;     //Huang Rhees parameter for k=1  
 double e_sigma_k_eq_1 =  .000470457E4;  // electronic sigma for k = 1 
 double p_sigma_k_eq_1 = 0.00016E4;  // phonon sigma for k = 1  

   

    TF1* alpha= new TF1("alpha", alpha_shab, 3.76255, 3.7685, 4);   // numerator function  first 3 peaks
     alpha->SetParameter(0, I_k_eq_1);
	 alpha->SetParameter(1, S_k_eq_1); 
	 alpha->SetParameter(2, e_sigma_k_eq_1);
 	 alpha->SetParameter(3,p_sigma_k_eq_1);
	  
     alpha->SetParName(0, "I_{0}");
	 alpha->SetParName(1, "S"); 
	 alpha->SetParName(2, "#sigma/cm^{-1}");
	 alpha->SetParName(3,"#sigma_{0}/cm^{-1}");	  
	//  alpha->SetParName(6,"offset");	 
	//TCanvas *c1 = new TCanvas("c1","alpha", 700, 500);
	  
    graph1->Fit("alpha", "R");                        // fit but not yet
    alpha->SetLineColor(kBlack);
    alpha-> Draw("SAME");
    
 
	  
    return 0;
    
   

	
  
}