#include <iostream>
#include <fstream>
#include <vector>
#include <algorithm>
#include <string>

#include "TFile.h"
#include "TCollection.h"
#include "TKey.h"
#include "TGraphErrors.h"
#include "TCanvas.h"
#include "TString.h"
#include "TF1.h"
#include "TVectorT.h"

#define A_RE43_B 1.966
#define A_RE43_TN 1.313
#define A_RE43_TW 0.555

#define A_RE42_B 1.169
#define A_RE42_TN 0.644
#define A_RE42_TW 0.458

#define Line0Low_TW 0.
#define Line0High_TW 1950.
#define Line1Low_TW 2100.
#define Line1High_TW 2300.

#define Line0Low_TN 0.
#define Line0High_TN 1820.
#define Line1Low_TN 2050.
#define Line1High_TN 2400.

#define Line0Low_B 0.
#define Line0High_B 1900.
#define Line1Low_B 2060.
#define Line1High_B 2180.

using namespace std;

void Argon_904( string inputfilename, char* type, int size )
{
  // strip off .root from filename
  // USAGE  .x Argon_904.cc+("inputfilename.dat", "B", 2);
  if ( inputfilename.substr(inputfilename.find_last_of(".")) == ".dat" )
    inputfilename = inputfilename.erase(inputfilename.find_last_of("."));

  vector<Float_t> voltages;
  vector<Float_t> currents;
  char cdummy[100];
  float voltage, current, fdummy;
  ifstream inputfile((inputfilename+".dat").c_str());
  while ( inputfile.good() ) {
    //    cdummy = inputfile.get();
    inputfile >> cdummy >> cdummy >> cdummy >> cdummy >> cdummy >> voltage >> current
	      >> fdummy >> fdummy >> fdummy >> fdummy >> fdummy;
    voltages.push_back(voltage);
    currents.push_back(current);
  }
  
  TGraph *g1 = new TGraph(voltages.size(), &voltages[0], &currents[0]);

  TCanvas *c1 = new TCanvas();
  c1->Draw();

  TF1 *Line0, *Line1;
  if ( strcmp(type, "B") == 0 ) {
    Line0 = new TF1("Line0","[0] + [1]*x", Line0Low_B, Line0High_B);
    Line1 = new TF1("Line1","[0] + [1]*x", Line1Low_B, Line1High_B);   
  } else if ( strcmp(type, "TW") == 0 ) {
    Line0 = new TF1("Line0","[0] + [1]*x", Line0Low_TW, Line0High_TW);
    Line1 = new TF1("Line1","[0] + [1]*x", Line1Low_TW, Line1High_TW);
  } else if ( strcmp(type, "TN") == 0 ) {
    Line0 = new TF1("Line0","[0] + [1]*x", Line0Low_TN, Line0High_TN);
    Line1 = new TF1("Line1","[0] + [1]*x", Line1Low_TN, Line1High_TN);
  }
  
  g1->Fit("Line0","NR");
  g1->Fit("Line1","NR");

  // Fit parameters for y=m*x+q
  Double_t q_left = Line0->GetParameter(0);
  Double_t m_left = Line0->GetParameter(1);
  Double_t q_right = Line1->GetParameter(0);
  Double_t m_right = Line1->GetParameter(1);

  cout << "q_left = " << q_left << endl;
  cout << "m_left = " << m_left << endl;
  cout << "q_right = " << q_right << endl;
  cout << "m_right = " << m_right << endl;

  // Useful values
  Double_t Rtot = 1./m_right;
  Double_t HV_Onset = (q_right-q_left)/(m_left-m_right);
  Double_t Resist = -10.;
  if ( strcmp(type, "B") == 0 ) {
    if ( size == 3 )
      Resist = 10.*(Rtot*A_RE43_B)/4.;
    else if ( size == 2 )
      Resist = 10.*(Rtot*A_RE42_B)/4.;
  }
  else if ( strcmp(type, "TW") == 0 ) {
    if ( size == 3 )
      Resist = 10.*(Rtot*A_RE43_TW)/4.;
    else if ( size == 2 )
      Resist = 10.*(Rtot*A_RE42_TW)/4.;
  }
  else if ( strcmp(type, "TN") == 0 ) {
    if ( size == 3 )
      Resist = 10.*(Rtot*A_RE43_TN)/4.;
    else if ( size == 2 )
      Resist = 10.*(Rtot*A_RE42_TN)/4.;
  }

  cout << "Rtot = " << Rtot << " [10^6 Ohm]" << endl;
  cout << "HV_Onset = " << HV_Onset << " [V]" << endl;
  cout << "Resist = " << Resist << " [10^10 Ohm cm]" << endl;

  g1->Draw("AP*");
  Line0->SetRange(0.,2600.);
  Line1->SetRange(1800.,2600.);
  Line0->SetLineColor(kRed);
  Line1->SetLineColor(kRed);
  Line0->Draw("same");
  Line1->Draw("same");

  c1->Print((inputfilename+".jpg").c_str());

}