//-------------------------------------
//  2016.4.1 recovere from old ROOT program by Taka Kondo
//           using ROOT 6.04/12
//  Electric field calculation,  Taka Kondo (KEK) 2009.11.27
//------------------------------------
#include "induced_charge.h"

void plot_Ey() {
//------------set main parameters and initialization----------------
int set_m_EfieldModel = 2;
double set_m_VD = 70.;
double set_m_VB = 150.;
initialize(set_m_EfieldModel, set_m_VD, set_m_VB);
//------------------------------
// 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);
//--------------------------------
double EX, EY, EYU, EyFEM;
double x, y;
int nhist=200;
double deltay = m_bulk_depth/nhist;
double XX[500], YY[3][500], YFEM[9][500];
//-----------------------------------------
int npoint = 0;
for (y=deltay*0.5; y<m_bulk_depth ; y+=deltay){
  x = 40./10000.;
  XX[npoint] = y * 10000.;
//---------------------Uniform field model-----------
  m_EfieldModel = 0;
  EField(x, y, EX, EYU);    // [V/cm]
  YY[0][npoint] = EYU;
//---------------------Flat Diode Model -------------
  m_EfieldModel = 1;
  EField(x, y, EX, EY);    // [V/cm]
  YY[1][npoint] = EY;
//---------------------FEM field --------------------
  m_EfieldModel = 2;
  for (int k=0; k<9; k++){
    x = (k*5.+40.)/10000.;
    EField(x, y, EX, EyFEM);    // [V/cm]
    YFEM[k][npoint] = EyFEM;
    if(k==0)cout<<"x,y[cm]=("<<x<<","<<y<<") : Ey(0)="<<EYU<<" Ey(1)="<<EY<<" Ey(2)="<<EyFEM<<endl;
  }
  npoint++;
}
//--------------plot Ey ----------------------
TGraph *hE, *hEU, *g1;
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("height [#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->GetXaxis()->SetTitleSize(0.04);
hE->GetXaxis()->SetLabelSize(0.04);
hE->GetYaxis()->SetTitleSize(0.04);
hE->GetYaxis()->SetLabelSize(0.04);
  //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(m_cid,"V_{D} = %4.0f V ,", m_VD);
TLatex *t=new TLatex(30., 9000.,m_cid);
//TLatex *t=new TLatex(20., 9200.,"V_{D} = 65 V , V_{B}= 150 V");
t->SetTextSize(0.04); t->Draw();
sprintf(m_cid,"V_{B} = %4.0f V", m_VB);
t=new TLatex(85., 9000.,m_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(m_cid," %d #mum",k*5);
  if (k<5)leg1->AddEntry(g1,m_cid,"P"); 
  if (k>4)leg2->AddEntry(g1,m_cid,"P"); 
}
leg1->Draw(); 
leg2->Draw();
sprintf(m_cid,"Ey_VD%d_VB%d.png",int(m_VD),int(m_VB));
c1->Print(m_cid);
}