#include <unistd.h>
#include <TH1D.h>
#include <TH2F.h>
#include <TF1.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TGraphErrors.h>
#include <TMultiGraph.h>
#include <TText.h>

#include <fstream>
#include <string>
#include <sstream>

float HV[2][44][2][10];// [plane][paddle][side][run]
int RUNS[10] = {10646,  10648,  10838,  10839,  10840,  
		10841,  10842,  10843,  10844,  10845};
float AttenPar[2][44][2][10][6]; //[plane][paddle][side][run][attenpar];
float NewHV[2][44][2];

int DEBUG  = 0;

int AmpGoalmV = 35;    // 45mV == 90 ADC : this is 0-2V dynamic range ADC
float AmpGoalADC = 70; // 35mV == 70 ADC  
float XDistIn_ns =7.;  // Distance from PMT in ns

void readHV(){
  //read HV from data files (generated from EPCIS archiver data)

  // run 10646 nomial HV
  // run 10648 -200 V
  // run 10838 -150 V
  // run 10839  -50 V
  // run 10840 -200 V
  // run 10841  -75 V
  // run 10842 -125 V
  // run 10843 -250 V
  // run 10844  -25 V
  // run 10845 -300 V
  
  char fnam[128];
  for (int n=0;n<10;n++){
    
    sprintf(fnam,"tof_hv_run%d.dat",RUNS[n]);
    cout<<fnam<<endl;
    ifstream INF(fnam);
    string dummy;
    getline(INF,dummy);
    //cout<<dummy;
    getline(INF,dummy);
    //cout<<dummy;
    istringstream strstr(dummy);
    string date,time;
    strstr>>date>>time;
    for (int k=0;k<44;k++){
      strstr>>HV[1][k][0][n];
      //    cout<<HV[1][k][0][n]<<endl;
    }
    for (int k=0;k<44;k++){
      strstr>>HV[1][k][1][n];
      //    cout<<HV[1][k][1][n]<<endl;
    }
    for (int k=0;k<44;k++){
      strstr>>HV[0][k][0][n];
      //    cout<<HV[0][k][0][n]<<endl;
    }
    for (int k=0;k<44;k++){
      strstr>>HV[0][k][1][n];
      //    cout<<HV[0][k][1][n]<<endl;
    }
    INF.close();
  }
}

void readAttenPar(){

  // read fit parameters from Attenuation fit to each PMT

  char fnam[128];

  for (int n=0;n<10;n++){
    
    sprintf(fnam,"genroot/hists/%06d/attenuation_run%d.dat",RUNS[n],RUNS[n]);
    cout<<fnam<<endl;
    ifstream INF(fnam);
    int plane,paddle,side;
    while (!INF.eof()){
      INF>>plane>>side>>paddle;
      INF>>AttenPar[plane][paddle][side][n][0];
      INF>>AttenPar[plane][paddle][side][n][1];
      INF>>AttenPar[plane][paddle][side][n][2];
      INF>>AttenPar[plane][paddle][side][n][3];
      INF>>AttenPar[plane][paddle][side][n][4];
      INF>>AttenPar[plane][paddle][side][n][5];
    }
  
    INF.close();
  }


}

float FunctionF(int plane, int paddle, int side, int rid){

  // evaluate atteniation FunctionF for PMT # = rid

  TF1 *f1 = new TF1("f1","[0]/2.*(exp(-(x*15.7+126.)/[1])+exp(-(x*15.7+126)/[2]))",-5.,5.);
  f1->SetParameter(0,AttenPar[plane][paddle][side][rid][0]);
  f1->SetParameter(1,AttenPar[plane][paddle][side][rid][2]);
  f1->SetParameter(2,AttenPar[plane][paddle][side][rid][4]);

  // negative value is close the the PMT poisitive value is away to the pmt

  float val = XDistIn_ns;
  float res = f1->Eval(val);

  return res;
}

void gains(){

  readHV();
  readAttenPar();


  int PLANE,PADDLE,SIDE;

  PLANE = 0;
  PADDLE = 15;
  SIDE = 0;

  float fitResults[176][4];
  float dfitResults[176][4];

  for (int PLA=0;PLA<2;PLA++){
    for (int SI=0;SI<2;SI++){
      for (int PAD=0;PAD<44;PAD++){
	PADDLE = PAD;
	PLANE = PLA;
	SIDE = SI;
	if ((PADDLE>20) && (PADDLE<23)){
	  continue;// skip over half paddles
	}
	
	int PMTID = PLA*88 + SI*44 + PAD;

	float HighVolate[10];
	float Amp[10];
	for (int k=0;k<10;k++){
	  HighVolate[k] = HV[PLANE][PADDLE][SIDE][k];
	  Amp[k] = FunctionF(PLANE,PADDLE,SIDE,k);
	  cout<<HighVolate[k]<<"   "<<Amp[k]<<endl;
	}
	
	TGraph *gr = new TGraph(10,HighVolate,Amp);
	gr->SetMarkerStyle(20);
	gr->SetMarkerColor(2);
	
	char tit[128];
	sprintf(tit,"Plane %d, Paddle %d, Side %d   (At %dns*15.7cm/ns+126cm)",PLANE+1, PADDLE+1, SIDE,(int)XDistIn_ns);
	gr->SetTitle(tit);
	
	gr->Draw("AP");
	gr->GetXaxis()->SetTitle("High Volate [V]");
	gr->GetYaxis()->SetTitle("MPV 1mV=2ADC [ADC]");
	gPad->SetGrid();
	
	float min = 9999.;
	float max = 0.;
	for (int j=0;j<10;j++){
	  if (HighVolate[j]<min){
	    min = HighVolate[j];
	  }
	  if (HighVolate[j]>max){
	    max = HighVolate[j];
	  }
	}
	min -=10.;
	max +=10.;
	
	TF1 *f1 = new TF1("f1","[0]+[1]*exp([2]*x+[3])",min, max);
	f1->SetParameter(0,1.);
	f1->SetParameter(1,0.00001);
	f1->SetParameter(2,0.0001);
	f1->SetParameter(3,1.);
	f1->SetLineColor(4);
	
	gr->Fit(f1,"","R",min, max);
	TF1 *pf = gr->GetFunction("f1");
	pf->SetLineColor(4);
	gPad->Update();
	
	TF1 *finv = new TF1("finv","1./[2]*(log(x-[0])-log([1])-[3])",0.,2000.);
	finv->SetParameter(0,pf->GetParameter(0));
	finv->SetParameter(1,pf->GetParameter(1));
	finv->SetParameter(2,pf->GetParameter(2));
	finv->SetParameter(3,pf->GetParameter(3));
	float SigAmpHV = finv->Eval(AmpGoalADC);
	NewHV[PLANE][PADDLE][SIDE] = SigAmpHV;

	for (int j=0;j<4;j++){
	  fitResults[PMTID][j] = pf->GetParameter(j);
	  dfitResults[PMTID][j] = pf->GetParError(j);
	}

	
	char tx1[128];
	sprintf(tx1,"at %dmV HV is %6.1f",AmpGoalmV,SigAmpHV);
	TText *t1 = new TText(0.2,0.6,tx1);
	t1->SetTextColor(2);
	t1->SetNDC();
	t1->SetTextSize(0.03);
	t1->Draw();
	gPad->Update();
	
	char anam[256];
	sprintf(anam,"plots/PMTgain_plane%d_paddle%d_side%d_amp%dmV_at%dns.pdf",PLANE+1, PADDLE+1, SIDE,AmpGoalmV,(int)XDistIn_ns);
	gPad->SaveAs(anam);
	sprintf(anam,"plots/PMTgain_plane%d_paddle%d_side%d_amp%dmV_at%dns.gif",PLANE+1, PADDLE+1, SIDE,AmpGoalmV,(int)XDistIn_ns);
	gPad->SaveAs(anam);
    
	if (DEBUG){
	  getchar();
	}

      }
      
    }
  }

  char ofnam[128];
  sprintf(ofnam, "newhv_for%dmV_at7ns.dat",AmpGoalmV);
  ofstream OUTF(ofnam);
  for (int k=0; k<2; k++){
    for (int s=0; s<2; s++){
      for (int p=0; p<44; p++){
	OUTF<<k<<"  "<<s<<"  "<<p<<"    "<<NewHV[k][p][s]<<endl;
      }
    }
  }
  OUTF.close();

  for (int k=0;k<176;k++){
    cout<<k<<" ";
    for (int n=0;n<4;n++){
      cout<< fitResults[k][n]<<"  "<<dfitResults[k][n]<<"  ";
    }
    cout<<endl;
  }
  
  
}