//script for root to visualize output files of CalcEField
//NEEDS TO BE COMPILED WITH ".L scriptname+"

#include <fstream>
#include <vector>
#include <string>
#include <cstdlib>
#include <iostream>
#include "TpcEFieldCyl.h"
#include "TVector3.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TH2D.h"
#include "TCanvas.h"
#include "TColor.h"


void plotEField() {


  //gROOT->Reset();
gROOT->SetStyle("Plain");
gStyle->SetPalette(1);

//define variables
double nomEField;
double minR;
double segR;
int nSegR;
double minZ;
double segZ;
int nSegZ;

TH2D* histoR;
TH2D* histoZ;

//read file
std::cout<<"Enter filename: ";
std::string filename;
getline(std::cin,filename);
std::cout<<std::endl;
std::string answer;
std::cout<<"Do you want to use the original binning? (y/n): ";
getline(std::cin, answer);
std::cout.flush();

std::ifstream infile(filename.c_str(), std::fstream::in);

while(answer!="n" && answer!= "y") {
  std::cout<<"\ninvalid input! Please repeat: ";
  getline(std::cin, answer);
}
 
if(answer=="y") {
if (!infile.good())
{
  std::cerr<<"Not Existing"<<std::endl;
}
//read geometry
infile>>nomEField>>minR>>segR>>nSegR>>minZ>>segZ>>nSegZ;
//read el field and
//fill histogramms

// WATCH OUT FOR PROPER VALUES IN YOUR FIELD FILES!!!

// segR = ((nSegR-1)*segR)/nSegR;
// segZ = ((nSegZ-1)*segZ)/nSegZ;

histoR = new TH2D(filename.c_str(),"Electric field in r direction [V/cm]",
			nSegZ,minZ,minZ+nSegZ*segZ,nSegR,minR,minR+nSegR*segR);
histoZ = new TH2D(filename.c_str(),"Electric field in z direction [V/cm]",
			nSegZ,minZ,minZ+nSegZ*segZ,nSegR,minR,minR+nSegR*segR);
for (int nr=0; nr<nSegR; nr++)
{
  for (int nz=0; nz<nSegZ; nz++)
  {
    double eFieldR; 
    double eFieldZ;
    infile>>eFieldR>>eFieldZ;
    histoR->SetBinContent(nz+1,nr+1,eFieldR);
    histoZ->SetBinContent(nz+1,nr+1,eFieldZ);
  }
}

}

if(answer=="n") {

std::ifstream infile(filename.c_str(), std::fstream::in);
if (!infile.good())
{
  std::cerr<<"Not Existing"<<std::endl;
}

//read geometry
infile>>nomEField>>minR>>segR>>nSegR>>minZ>>segZ>>nSegZ;

// eval max coordinates:
double maxR = minR + segR*(nSegR-1);
double maxZ = minZ + segZ*(nSegZ-1);

//get new binning:
std::cout<<"\nEnter bins in r direction: ";
getline(std::cin, answer);
nSegR=atoi(answer.c_str());
std::cout<<"Enter bins in z direction: ";
getline(std::cin, answer);
nSegZ=atoi(answer.c_str());

//eval new binning
segR = (maxR - minR)/nSegR;
segZ = (maxZ - minZ)/nSegZ;


//create and fill histograms
TString titleR("Electric field in r direction (V/cm), ");
TString titleZ("Electric field in z direction (V/cm), ");
titleR.Append(filename.c_str());
titleZ.Append(filename.c_str());


histoR = new TH2D(filename.c_str(),titleR,
			nSegZ,minZ,minZ+(nSegZ-1)*segZ,nSegR,minR,minR+(nSegR-1)*segR);
histoZ = new TH2D(filename.c_str(),titleZ,
			nSegZ,minZ,minZ+(nSegZ-1)*segZ,nSegR,minR,minR+(nSegR-1)*segR);

TpcEFieldCyl* field = new TpcEFieldCyl(filename.c_str());

for (int nr=0; nr<nSegR; nr++)
{
  for (int nz=0; nz<nSegZ; nz++)
  {
    double eFieldR; 
    double eFieldZ;
    TVector3 pos = TVector3(minR + (nr+0.5)*segR, 0., minZ + (nz+0.5)*segZ);
    TVector3 value=field->value(pos);
    eFieldR=value.X();
    eFieldZ=value.Z();
    histoR->SetBinContent(nz+1,nr+1,eFieldR);
    histoZ->SetBinContent(nz+1,nr+1,eFieldZ);
  }
}
}

infile.close();

Int_t nb=200;
//TColor::CreateGradientColorTable(Number,Stops,Red,Green,Blue,nb);
TColor::SetPalette(1,0);

gStyle->SetOptStat(0);

//draw
 TString cfile(filename.c_str());
 TString rfile = cfile;
 rfile.Append("_r");
TCanvas* c1 = new TCanvas(rfile,rfile,1000,350);
histoR->GetXaxis()->SetTitle("z (cm)");
histoR->GetYaxis()->SetTitle("r (cm)");
histoR->SetContour(nb);
histoR->Draw("COLZ");
 TString zfile = cfile;
 rfile.Append("_z");
 TCanvas* c2= new TCanvas(zfile,zfile,1000,350);
histoZ->GetXaxis()->SetTitle("z (cm)");
histoZ->GetYaxis()->SetTitle("r (cm)");
histoZ->SetContour(nb);
histoZ->Draw("COLZ");


}