#include <iostream>
#include <iomanip>
#include <fstream>
#include <iostream>
#include <sstream>

// #include "Riostream.h" // what is this?

#include "TString.h"
#include "TFile.h"
#include "TNtuple.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TPad.h"
#include "TAxis.h"
#include "TLatex.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TSystem.h"
#include "TROOT.h" 
#include "TStyle.h"
#include "TProfile.h"
#include "TF1.h"
#include "TPaveStats.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TVectorD.h"
#include "TPad.h"

void savePlot(TCanvas*, TString);
void TownsendCoefficient_doublefit_PScan_Automatic() {

   gROOT->SetStyle("Plain");
   // TStyle *gStyle = gROOT->GetStyle("Plain");
   gStyle->SetPalette(1,0);
   gStyle->SetOptStat(000000);
   gStyle->SetOptFit(000000);
   // gStyle->SetOptFit(111111);
   // gStyle->SetOptStat(10);
   // gStyle->SetOptFit(1011);
   gStyle->SetOptTitle(0);
   gStyle->SetGridStyle(1);
   gROOT->UseCurrentStyle();
   // gPad->SetTicks();

   string filenames[] = {"ArCO2_70_30_94000Pa_Townsend.txt"};  // test with just one file
   // string filenames[] = {"ArCO2_70_30_94000Pa_Townsend.txt", "ArCO2_70_30_95000Pa_Townsend.txt", "ArCO2_70_30_96000Pa_Townsend.txt", "ArCO2_70_30_97000Pa_Townsend.txt", "ArCO2_70_30_98000Pa_Townsend.txt","ArCO2_70_30_99000Pa_Townsend.txt", "ArCO2_70_30_100000Pa_Townsend.txt", "ArCO2_70_30_101000Pa_Townsend.txt", "ArCO2_70_30_101325Pa_Townsend.txt", "ArCO2_70_30_102000Pa_Townsend.txt" };
   double pressure[] = {94000}; // test with just one file
   // double pressure[] = {94000, 9500, 96000, 97000, 98000, 99000,  100000, 101000, 101325, 102000};


   // declare some vectors here to save data
   // ======================================
   const int size   = 1;
   const int n      = 20;
   const int digits = 2;
   
   std::vector< double > dummy; 
   for(int i=0; i<n; ++i) { dummy.push_back(0); }
   std::vector< std::vector<double> > x, y, p, e;
   for(int i=0; i<size; ++i) { x.push_back(dummy); y.push_back(dummy); p.push_back(dummy); e.push_back(dummy); }

   std::vector< std::string > plotname;
   std::vector< std::vector< TGraphErrors* > > tgraphs;
   std::vector< TF1* > theor_korff;
   std::vector< std::vector< TF1* > > fit_tf1s;
   std::vector< std::vector< std::string > > fit_names;
   std::vector< std::vector< TGraphErrors* > > tgraphs_res;
   std::vector< std::vector< TVectorD > > fit_params;

   // declare canvas here for combined data
   // =====================================
   // TCanvas *d = new TCanvas("d","d",0,0,800,600);
   // TCanvas *e = new TCanvas("e","e",0,0,800,600);


   // declare common settings (colors,markers,ranges, etc)
   // ====================================================

   // Fit to Low E-field values: technical trick: reduce the size n1 
   // such that first n1 elements of array n are considered
   int n1 = 10;
   // Fit to High E-field values: technical trick: select subarrray starting at element n1
   // int *arr2 = arr1 + n1;
   int overlap = 0;
   int n2 = n-n1+1+overlap;
   std::vector< int > nk; nk.push_back(n); nk.push_back(n1); nk.push_back(n2);

   // plot & fit ranges
   int RangeTot1  = 1.5E6, RangeTot2  = 1.65E7;
   int RangeLow1  = 1.5E6, RangeLow2  = 7E6;
   int RangeHigh1 = 7E6,   RangeHigh2 = 1.65E7;

   // colors and markers for single pressure
   // 1) tgrapherrors with data
   int dfit_tgraph_linecolor[] = {1,1,1};   int dfit_tgraph_linewidth[] = {2,2,2};
   int dfit_tgraph_markcolor[] = {1,1,1}; double dfit_tgraph_marksize[] = {1.25, 1.75, 2.75}; int dfit_tgraph_markstyle[] = {20,25,27};
   // 2) theoretical korff curve
   int dfit_tkorff_linecolor = 1; int dfit_tkorff_linewidth = 1; int dfit_tkorff_linestyle = 2;
   // 3) fitted korff curves
   int dfit_fit_linecolor[] = {2,4,6}; int dfit_fit_linewidth[] = {2,2,2}; int dfit_fit_linestyle[] = {1,1,1};

   // colors and markers for overview plot



   // Loop over all input files
   // and save data in vectors
   // =========================
   for(int i=0; i<size; ++i) {
     // read inputfile and save info
     ifstream in;
     in.open((filenames[i]).c_str());
     Int_t nlines = 0;
     std::cout<<"here still ok"<<std::endl;
     while (1) {
       in >> x[i][nlines] >> y[i][nlines] >> p[i][nlines];
       y[i][nlines]=y[i][nlines]*100;
       e[i][nlines]=y[i][nlines]*p[i][nlines]*1.0/100;
       if (!in.good()) break;
       if (nlines < 5) printf("x=%8f, y=%8f, e=%8f\n",x[i][nlines],y[i][nlines],e[i][nlines]);
       nlines++;
     }
     in.close();
   }
    

   // Loop over all vectors
   // and make tgraphs and fits
   // =========================
   for(int i=0; i<size; ++i) {

     // declare canvas
     TCanvas *c = new TCanvas("c","c",0,0,800,600);

     // declare name
     std::stringstream name_ss; name_ss<<"TownsendCoefficient_20k_160k_"<<pressure[i]<<"Pa_doublefit"; std::string name_s = name_ss.str();
     plotname.push_back(name_s);

     // tgraph errors with data points
     TGraphErrors *garfield_tot  = new TGraphErrors(n, &x[i][0],&y[i][0],0,&e[i][0]);
     TGraphErrors *garfield_low  = new TGraphErrors(n1,&x[i][0],&y[i][0],0,&e[i][0]);
     double *x_high = &x[i][0] + n1-1-overlap, *y_high = &y[i][0] + n1-1-overlap, *e_high = &e[i][0] + n1-1-overlap;
     TGraphErrors *garfield_high = new TGraphErrors(n2,x_high,y_high,0,e_high);
     std::vector< TGraphErrors* > tgraph_loop; 
     tgraph_loop.push_back(garfield_tot); tgraph_loop.push_back(garfield_low); tgraph_loop.push_back(garfield_high);
     tgraphs.push_back(tgraph_loop);
     for(int j=0; j<3; ++j) {
       tgraphs[i][j]->SetLineColor(dfit_tgraph_linecolor[j]);   tgraphs[i][j]->SetLineWidth(dfit_tgraph_linewidth[j]);
       tgraphs[i][j]->SetMarkerColor(dfit_tgraph_markcolor[j]); tgraphs[i][j]->SetMarkerSize(dfit_tgraph_marksize[j]); tgraphs[i][j]->SetMarkerStyle(dfit_tgraph_markstyle[j]);
     }

     // theoretical korff curve
     std::stringstream korff_ArCO2_ss; korff_ArCO2_ss << "korff_ArCO2_" << pressure[i]; std::string korff_ArCO2_s = korff_ArCO2_ss.str();
     TF1 *korff_ArCO2 = new TF1(korff_ArCO2_s.c_str(),"[0]*[2]*exp(-[1]*[2]/x)",RangeTot1,RangeTot2);
     korff_ArCO2->SetParameter(0,11.85);   // Ar/CO2 70/30 :: A = 11.85  Pa^-1 m^-1
     korff_ArCO2->SetParameter(1,199.37);  // Ar/CO2 70/30 :: B = 199.37 Pa^-1 m^-1 V
     korff_ArCO2->FixParameter(2,pressure[i]);
     korff_ArCO2->SetLineColor(dfit_tkorff_linecolor); korff_ArCO2->SetLineWidth(dfit_tkorff_linewidth); korff_ArCO2->SetLineStyle(dfit_tkorff_linestyle); 
     theor_korff.push_back(korff_ArCO2);

     // fitted korff curves
     std::stringstream korff_garf_tot_ss;  korff_garf_tot_ss << "korff_garfield_tot_" << pressure[i];   std::string korff_garf_tot_s = korff_garf_tot_ss.str();
     std::stringstream korff_garf_low_ss;  korff_garf_low_ss << "korff_garfield_low_" << pressure[i];   std::string korff_garf_low_s = korff_garf_low_ss.str();
     std::stringstream korff_garf_high_ss; korff_garf_high_ss << "korff_garfield_high_" << pressure[i]; std::string korff_garf_high_s = korff_garf_high_ss.str();
     TF1 *korff_garfield_tot  = new TF1(korff_garf_tot_s.c_str(), "[0]*[2]*exp(-[1]*[2]/x)",RangeTot1,RangeTot2);
     TF1 *korff_garfield_low  = new TF1(korff_garf_low_s.c_str(), "[0]*[2]*exp(-[1]*[2]/x)",RangeLow1,RangeLow2);
     TF1 *korff_garfield_high = new TF1(korff_garf_high_s.c_str(),"[0]*[2]*exp(-[1]*[2]/x)",RangeHigh1,RangeHigh2);
     std::vector< std::string > str_loop;
     str_loop.push_back(korff_garf_tot_s); str_loop.push_back(korff_garf_low_s); str_loop.push_back(korff_garf_high_s); 
     fit_names.push_back(str_loop);
     std::vector< TF1* > tf1_loop; 
     tf1_loop.push_back(korff_garfield_tot); tf1_loop.push_back(korff_garfield_low); tf1_loop.push_back(korff_garfield_high);
     fit_tf1s.push_back(tf1_loop);
     for(int j=0; j<3; ++j) {
       fit_tf1s[i][j]->SetParameter(0,10); fit_tf1s[i][j]->SetParameter(1,130); // initial values
       fit_tf1s[i][j]->FixParameter(2,pressure[i]); // P = 95000 Pa
       fit_tf1s[i][j]->SetLineColor(dfit_fit_linecolor[j]); fit_tf1s[i][j]->SetLineWidth(dfit_fit_linewidth[j]); fit_tf1s[i][j]->SetLineStyle(dfit_fit_linestyle[j]);
       std::cout<<" Fit tgraphs["<<i<<"]["<<j<<"] with TF1 named "<<fit_names[i][j]<<std::endl;
       // tgraphs[i][j]->Fit((fit_names[i][j]).c_str(),"R"); // puzzling "Warning in <Fit>: Fit data is empty"
       tgraphs[i][j]->Fit(fit_tf1s[i][j],"R");               // puzzling "Warning in <Fit>: Fit data is empty"
       // tgraphs[i][j]->Fit(fit_tf1s[i][j],"WW");           // experimenting 
       std::cout<<" Check tgraphs["<<i<<"]["<<j<<"] :: bins = "<<nk[j]<<" bincontent: ";
       for(int k=0; k<nk[j]; ++k) {
	 if(j<2) std::cout<<" bin ["<<k<<"][ x = "<<x[i][k]<<" y = "<<tgraphs[i][j]->Eval(x[i][k])<<"]";
	 else std::cout<<" bin ["<<k<<"][ x = "<<x[i][k+n1-1-overlap]<<" y = "<<tgraphs[i][j]->Eval(x[i][k+n1-1-overlap])<<"]";
       } std::cout<<std::endl;
     }

     // fit residuals
     double garfield_tot_res_val[n], garfield_tot_res_err[n], garfield_low_res_val[n1], garfield_low_res_err[n1], garfield_high_res_val[n2], garfield_high_res_err[n2];
     for(int k=0; k<n; ++k)  {garfield_tot_res_val[k]  = (y[i][k]   - fit_tf1s[i][0]->Eval(x[i][k]))/y[i][k];}
     for(int k=0; k<n1; ++k) {garfield_low_res_val[k]  = (y[i][k]   - fit_tf1s[i][1]->Eval(x[i][k]))/y[i][k];}
     for(int k=0; k<n2; ++k) {garfield_high_res_val[k] = (y_high[k] - fit_tf1s[i][2]->Eval(x_high[k]))/y_high[k];}
     TGraphErrors *garf_res_tot  = new TGraphErrors(n, &x[i][0],     garfield_tot_res_val,0,0);
     TGraphErrors *garf_res_low  = new TGraphErrors(n1,&x[i][0],     garfield_low_res_val,0,0);
     TGraphErrors *garf_res_high = new TGraphErrors(n2,x_high,garfield_high_res_val,0,0);
     std::vector< TGraphErrors* > tgraph_res_loop; tgraph_res_loop.push_back(garf_res_tot); tgraph_res_loop.push_back(garf_res_low); tgraph_res_loop.push_back(garf_res_high);
     tgraphs_res.push_back(tgraph_res_loop);
     for(int j=0; j<3; ++j) {
       tgraphs_res[i][j]->SetLineColor(dfit_tgraph_linecolor[j]);   tgraphs_res[i][j]->SetLineWidth(dfit_tgraph_linewidth[j]);
       tgraphs_res[i][j]->SetMarkerColor(dfit_tgraph_markcolor[j]); tgraphs_res[i][j]->SetMarkerSize(dfit_tgraph_marksize[j]); tgraphs_res[i][j]->SetMarkerStyle(dfit_tgraph_markstyle[j]);
     }

     // fit params
     TVectorD c_graph1(7), c_graph2(7), c_graph3(7);
     std::vector< TVectorD > param_loop; param_loop.push_back(c_graph1); param_loop.push_back(c_graph2); param_loop.push_back(c_graph3);
     fit_params.push_back(param_loop);
     for(int j=0; j<3; ++j) {
       fit_params[i][j](0) = fit_tf1s[i][j]->GetParameter(0);
       fit_params[i][j](1) = fit_tf1s[i][j]->GetParameter(1);
       fit_params[i][j](2) = fit_tf1s[i][j]->GetParameter(2);
       fit_params[i][j](3) = fit_tf1s[i][j]->GetChisquare();
       fit_params[i][j](4) = fit_tf1s[i][j]->GetNDF();
       fit_params[i][j](5) = fit_tf1s[i][j]->GetParError(0);
       fit_params[i][j](6) = fit_tf1s[i][j]->GetParError(1);
     }

     // legends
     std::vector< std::string > legend_strings;
     for(int j=0; j<3; ++j) {
       std::stringstream legss1, legss2, legss3; legss1.precision(digits); legss2.precision(digits); legss3.precision(digits);
       legss1<<fixed<<"A = "<<fit_params[i][j](0)<<" #pm "<<fit_params[i][j](5)<<" [ m^{-1} Pa^{-1} ]";
       legss2<<fixed<<"B = "<<fit_params[i][j](1)<<" #pm "<<fit_params[i][j](6)<<" [ V m^{-1} Pa^{-1} ]";
       legss3<<"P = "<<fit_params[i][j](2)<<" [ Pa ] #chi^{2}/N_{dof} = "<<fixed<<fit_params[i][j](3)/fit_params[i][j](4);
       std::string legstr1 = legss1.str(); legend_strings.push_back(legstr1);
       std::string legstr2 = legss2.str(); legend_strings.push_back(legstr2);
       std::string legstr3 = legss3.str(); legend_strings.push_back(legstr3);
     }

     // make plot
     c->cd();
     TPad *topPad = new TPad("TopPad","",0,0,1,1);
     topPad->SetFillStyle(4000);
     topPad->SetFillColor(4000);
     topPad->SetFrameFillColor(4000);
     topPad->SetFrameFillStyle(4000);
     topPad->SetFrameBorderMode(0);
     topPad->SetTicks(1,1);
     topPad->SetGrid(0,1);
     topPad->Draw();
     topPad->cd();
     topPad->SetTitle("Townsend Coefficient (m^{-1})");
     topPad->SetLogy(1);
     // c->SetLogx(1);
     tgraphs[i][0]->Draw("AP");
     tgraphs[i][0]->GetYaxis()->SetTitle("#eta (m^{-1})");
     tgraphs[i][0]->GetYaxis()->SetTitleOffset(1.25);
     tgraphs[i][0]->GetYaxis()->SetRangeUser(50,1E6);
     // tgraphs[i][0]->GetXaxis()->SetTitle("E (V m^{-1})");
     tgraphs[i][0]->GetXaxis()->SetLabelSize(0.001);
     tgraphs[i][0]->GetXaxis()->SetRangeUser(RangeTot1,RangeTot2);
     // tgraphs[i][0]->GetXaxis()->SetMoreLogLabels(1);
     // tgraphs[i][0]->GetYaxis()->SetMoreLogLabels(1);
     theor_korff[i]->Draw("same");
     topPad->Update();
     tgraphs[i][0]->Draw("sameEP");
     tgraphs[i][1]->Draw("sameEP");
     tgraphs[i][2]->Draw("sameEP");
     topPad->Update();
     // leg->Draw();
     topPad->Update();

     double canvasratio = 0.2; // so ratio canvas is last 20% of your canvas
     topPad->SetBottomMargin(canvasratio + (1-canvasratio)*topPad->GetBottomMargin()-canvasratio*topPad->GetTopMargin());
     TPad *bottomPad = new TPad("BottomPad","",0,0,1,1);
     bottomPad->SetTopMargin((1-canvasratio) - (1-canvasratio)*bottomPad->GetBottomMargin()+canvasratio*bottomPad->GetTopMargin());
     bottomPad->SetFillStyle(4000);
     bottomPad->SetFillColor(4000);
     bottomPad->SetFrameFillColor(4000);
     bottomPad->SetFrameFillStyle(4000);
     bottomPad->SetFrameBorderMode(0);
     bottomPad->SetTicks(1,1);
     bottomPad->SetGrid(0,1);
     bottomPad->Draw();
     bottomPad->cd();
     tgraphs_res[i][0]->GetYaxis()->SetLabelSize(0.03);
     tgraphs_res[i][0]->GetYaxis()->SetNdivisions(504);
     tgraphs_res[i][0]->GetYaxis()->SetTitle("#frac{data-fit}{data}");
     // tgraphs_res[i][0]->GetYaxis()->SetTitleSize(1.00);
     tgraphs_res[i][0]->GetXaxis()->SetTicks("+");
     tgraphs_res[i][0]->GetXaxis()->SetTitle("E (V m^{-1})");
     tgraphs_res[i][0]->GetXaxis()->SetTitleOffset(1.15);
     tgraphs_res[i][0]->GetXaxis()->SetRangeUser(RangeTot1,RangeTot2);
     tgraphs_res[i][0]->GetYaxis()->SetRangeUser(-0.075,+0.075);
     tgraphs_res[i][0]->Draw("ACP");
     tgraphs_res[i][1]->Draw("CPsame");
     tgraphs_res[i][2]->Draw("CPsame");
     c->cd();
     c->Update();
     savePlot(c,plotname[i].c_str());
   }


   /*
  TLegend * leg = new TLegend(0.45,0.15,0.85,0.65);
  // leg->SetHeader("Townsend Coefficient");
  leg->AddEntry("korff_ArCO2_0950",   "Korff Ar/CO2 (70/30) 950 mbar","l");
  leg->AddEntry(garfield_0950,        "GARFIELD 950 mbar", "p");
  leg->AddEntry(korff_garfield_0950,  " 950 mbar Fit All", "l");
  leg->AddEntry(korff_garfield_0950A, " 950 mbar Fit Low", "l");
  leg->AddEntry(korff_garfield_0950B, " 950 mbar Fit High", "l");
  leg->AddEntry("korff_ArCO2_0980",   "Korff Ar/CO2 (70/30) 980 mbar","l");
  leg->AddEntry(garfield_0980,        "GARFIELD 980 mbar", "p");
  leg->AddEntry(korff_garfield_0980,  " 980 mbar Fit All", "l");
  leg->AddEntry(korff_garfield_0980A, " 980 mbar Fit Low", "l");
  leg->AddEntry(korff_garfield_0980B, " 980 mbar Fit High", "l");
  leg->AddEntry("korff_ArCO2_1010",   "Korff Ar/CO2 (70/30) 1010 mbar","l");
  leg->AddEntry(garfield_1010,        "GARFIELD 1010 mbar", "p");
  leg->AddEntry(korff_garfield_1010,  "1010 mbar Fit All", "l");
  leg->AddEntry(korff_garfield_1010A, "1010 mbar Fit Low", "l");
  leg->AddEntry(korff_garfield_1010B, "1010 mbar Fit High", "l");
  leg->Draw();

  savePlot(c,"TownsendCoefficient_20k_160k_PScan_doublefit");
   */
  /*
  d->cd();
  d->SetTitle("Korff Param A (m^{-1} Pa^{-1})");
  // d->SetLogy(1);
  // d->SetLogx(1);
  A_All->Draw("ACP");
  A_All->GetXaxis()->SetTitle("P (Pa)");
  A_All->GetYaxis()->SetTitle("Korff Param A (m^{-1} Pa^{-1})");
  // A_All->GetXaxis()->SetMoreLogLabels(1);
  // A_All->GetYaxis()->SetMoreLogLabels(1);
  A_Low->Draw("CPsame");
  A_High->Draw("CPsame");
  A_All->GetYaxis()->SetRangeUser(0,10);
  d->Update();
  savePlot(d,"TownsendCoefficient_20k_160k_PScan_doublefit_Korff_A");

  e->cd();
  e->SetTitle("Korff Param B (V m^{-1} Pa^{-1})");
  // e->SetLogy(1);
  // e->SetLogx(1);
  B_All->Draw("ACP");
  B_All->GetXaxis()->SetTitle("P (Pa)");
  B_All->GetYaxis()->SetTitle("Korff Param B (V m^{-1} Pa^{-1})");
  // B_All->GetXaxis()->SetMoreLogLabels(1);
  // B_All->GetYaxis()->SetMoreLogLabels(1);
  B_Low->Draw("CPsame");
  B_High->Draw("CPsame");
  B_All->GetYaxis()->SetRangeUser(80,160);
  e->Update();
  savePlot(e,"TownsendCoefficient_20k_160k_PScan_doublefit_Korff_B");
  */
}

void savePlot(TCanvas * c, TString name) {
  c->SaveAs(name+".png");
  c->SaveAs(name+".eps");
  gSystem->Exec("epstopdf "+name+".eps");
  std::cout<<"Plot "<<name<<" saved"<<std::endl;
}