#include<TMath.h>
#include<fstream>
#include<string>
#include<sstream>
#include<iostream>
#include<Math/Polynomial.h>
#include<Math/GSLRootFinder.h>
#include<complex>
#include<vector>

//max should be 670
void graphit(Char_t* filename1, Char_t* filename2, const double thick1, const double thick2){
  gROOT->Reset();
  gInterpreter->ExecuteMacro("../macros/daya_style.C");

  TCanvas *c1 = new TCanvas("c1", "Drift Plot", 0, 0, 600, 400);
  c1->Clear();
  
  const int max = 670;
  Double_t IOR[max];
  Double_t wavelength[max];

  Double_t wavelength1[max];
  Double_t transmittance1[max];

  Double_t wavelength2[max];
  Double_t transmittance2[max];

  Double_t difference[max];
  Double_t Gamma[max];
  Int_t j = 0;

  ifstream datafile1(filename1, ios::in);
  ifstream datafile2(filename2, ios::in);

  for (int i=0; i<max ; i++) {
    
    datafile1 >> wavelength1[i] >> transmittance1[i];
    datafile2 >> wavelength2[i] >> transmittance2[i];
        
    if(transmittance2[i]<.075)
      Gamma[i]=0;
    else if(transmittance1[i]/transmittance2[i]<.00001)
	Gamma[i]=0;
    else 
      Gamma[i]=(thick2-thick1)/TMath::Log(transmittance1[i]/transmittance2[i]);
    
    if(Gamma[i]==0)
      IOR[i] = 0;
    else
      {
	ROOT::Math::Polynomial poly(4); 
	double coeffs[5] = {1,4,6,4,1}; 
	poly.SetParameters(coeffs); 
	vector<complex<double> > roots; 
	roots = poly.FindNumRoots(); 
	std::cout << "number of roots = " << roots.size() << std::endl;
	for (unsigned int j = 0; j < roots.size(); ++j ) 
	  std::cout << roots[j] << std::endl;
      
	
	IOR[i] = roots[3];
      }

      cout<<"2in"<<"    "<<wavelength1[i]<<"   "<<transmittance2[i]<<"    "<<Gamma[i]<<"   "<<IOR[i]<<endl;

  }
  
  TGraph *g1 = new TGraph(max,wavelength1,Gamma);
  
  g1->SetTitle("Attenuation Length: Acrylic UVT:  Comparing 3/8\" and 2\" Sample");
  g1->SetMarkerStyle(7);
  g1->SetMarkerColor(kBlack);
  g1->SetLineColor(kBlack);

  g1->GetXaxis()->SetTitle("Wavelength (nm)");
  g1->GetYaxis()->SetTitle("Absorption Lenth (In)");
  g1->Draw("APL");
  c1->SaveAs("AbsorptionLength.ps");
}