//-------------------------------------
//  2016.2.17 recovere from old ROOT program by Taka Kondo
//            using ROOT 6.02/12
//
//------------------------------------
#include <iostream>
#include <math.h>
#include "sub_functions_v7.C"
#include "induced_charge.h"
#include "GetExEy_150V.cxx"
#include "GetPotential.cxx"

void plot_Ey() {
initialize();
// Electric field calculation,  Taka Kondo (KEK) 2009.11.27
// unit: nsec, cm, sec, K, 
TCanvas *c1 = new TCanvas("c1","A Simple Graph with error bars",0,0,700,500);
c1->SetFillStyle(4100);
c1->SetGrid();
c1->GetFrame()->SetBorderSize(12);
char cid[50];
double d = 0.0285;           // [cm]
double Vd = 70.0;            // [V] 
double Vb;                      // [V]
double EY, EYU, Emax, y, wd, EX, EyFEM;
int nhist=200;
double deltay=d/nhist;
double XX[500], YY[3][500], YFEM[9][500];
TGraph *hE, *hEU, *g1;
//-----------------------------------------
Vb=150.;
int npoint = 0;
for (y=deltay*0.5; y<d ; y+=deltay){
  if (Vb>Vd) {
     Emax = (Vb+Vd)/d ;
     EY = Emax - 2.*Vd*y/(d*d); //flat diode model
     EYU = Vb/d; //uniform field model
  } else {
     wd = sqrt(Vb/Vd) * d;
     Emax = 2.* wd * Vd /(d*d);  
     EY = Emax*(1-y/wd);   //flat diode model
     EYU = Vb/wd; //uniform field model
  }
  YY[1][npoint] = EY;
  double depth = (d-y)*10000.;
  XX[npoint] = depth;
  YY[0][npoint] = EYU;
  for (int k=0; k<9; k++){
    double x = k*5.+40.;
    EField(x/10000., depth/10000., EX, EyFEM);    // [V/cm]
    YFEM[k][npoint] = EyFEM;
  }
  cout<<(d-y)*10000.<<" "<<EY<<" "<<EYU<<" "<<EyFEM<<endl;
  npoint++;
}
hE = new TGraph(npoint, XX, YY[1]) ;  // flat diode model
hEU = new TGraph(npoint, XX, YY[0]) ; // uniform field model
gStyle->SetOptStat(0); 
hE->SetMaximum(10000.);
hE->SetMinimum(0.);
hE->SetTitle(" ");
hE->GetXaxis()->SetTitle("depth [#mum]");
hE->GetYaxis()->SetTitle("Vertical Electric Field [V/cm]");
hE->GetYaxis()->SetTitleOffset(1.5);
hE->GetXaxis()->SetTitleOffset(1.0);
hE->GetXaxis()->SetLimits(0.,285.);
  //hE ->Draw("L");  
hE ->Draw("AL");  
hEU->SetLineColor(2);
hEU ->Draw("SAME L"); 
//
TPave *b = new TPave(15.,6600, 145., 9600.); b->SetFillColor(10);
b->SetBorderSize(2); b->SetLineColor(1); b->Draw();
sprintf(cid,"V_{D} = %4.0f V ,", Vd);
TLatex *t=new TLatex(30., 9000.,cid);
//TLatex *t=new TLatex(20., 9200.,"V_{D} = 65 V , V_{B}= 150 V");
t->SetTextSize(0.04); t->Draw();
sprintf(cid,"V_{B} = %4.0f V", Vb);
t=new TLatex(85., 9000.,cid);
t->SetTextSize(0.04); t->Draw();
t=new TLatex(20., 8300.,"N_{eff}   [cm^{-3}]");
t->SetTextSize(0.04); t->Draw();
t=new TLatex(25., 7600.,"1.03x10^{12} (Flat diode model)");
t->SetTextSize(0.04); t->Draw();
t=new TLatex(25., 7000.,"8.86x10^{11} (FEM field model)");
t->SetTextSize(0.04); t->Draw();
//
t=new TLatex(30., 5500.,"Uniform field model");
t->SetTextAngle(0); t->SetTextColor(2); t->SetTextSize(0.04); t->Draw();
t=new TLatex(30., 3700.,"Flat diode model");
t->SetTextAngle(18); t->SetTextSize(0.04); t->Draw();
t=new TLatex(110.,3700.,"FEM solutions");
t->SetTextColor(4); t->SetTextAngle(14); t->SetTextSize(0.04); t->Draw();
//
TLegend *leg, *leg1, *leg2;
leg = new TLegend(0.58 ,0.88,0.81,0.93);
leg->SetHeader("distance from strip center");leg->SetBorderSize(0);
leg->SetFillColor(10);leg->SetTextSize(0.035);leg->Draw();
leg1 = new TLegend(0.58 ,0.67 ,0.71,0.88);
leg1->SetFillColor(10);leg1->SetBorderSize(0);leg1->SetTextFont(132);
leg2 = new TLegend(0.71 ,0.71 ,0.84,0.88);
leg2->SetFillColor(10);leg2->SetBorderSize(0);leg2->SetTextFont(132);
//
//c1->Print("E_test.pdf");
//
for (int k=0; k<9; k++){
  g1= new TGraph(npoint,XX,YFEM[k]);
  g1->SetMarkerStyle(7);
  g1->SetMarkerColor(k+1);
  g1->Draw("SAME P");
  sprintf(cid," %d #mum",k*5);
  if (k<5)leg1->AddEntry(g1,cid,"P"); 
  if (k>4)leg2->AddEntry(g1,cid,"P"); 
}
leg1->Draw(); 
leg2->Draw();
c1->Print("Ey_Vd70_Vb150.png");
}