{
#include <iostream>
#include <math.h>
// Electric field calculation,  Taka Kondo (KEK) 2009.11.27
// unit: nsec, cm, sec, K, 
c1 = new TCanvas("c1","A Simple Graph with error bars",0,0,700,500);
c1->SetFillStyle(4100);
c1->SetGrid();
c1->GetFrame()->SetBorderSize(12);
double d = 0.0285;           // [cm]
double Vd = 65.0;            // [V] 
double Vb;                      // [V]
double EY, Emax, y, wd;
int nhist=200;
double deltay=d/nhist;

TH1F *hE = new TH1F("hE","Electric field",nhist,0,d*10000.); 
TH1F *hEU = new TH1F("hEU","Uniform Electric field",nhist,0,d*10000.); 

   Vb=65.;
   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
      }
      hE->Fill((d-y)*10000.,EY);
      hEU->Fill((d-y)*10000.,EYU);
   }
  gStyle->SetOptStat(0); 
  hE->SetMaximum(10000.);
  hE->SetMinimum(0.);
  hE->SetTitle();
  hE->SetXTitle("depth [#mum]");
  hE->SetYTitle("Vertical Electric Field [V/cm]");
  hE->GetYaxis()->SetTitleOffset(1.5);
  hE->GetXaxis()->SetTitleOffset(1.0);
  hE ->Draw("L");  
  hEU->SetLineStyle(2);
  hEU ->Draw("SAME L"); 
//
b = new TPave(15.,6600, 145., 9600.); b->SetFillColor(10);
b->SetBorderSize(2); b->SetLineColor(1); b->Draw();
cid = new char[15];sprintf(cid,"V_{D} = %d V ,", Vd);cid[14]=0;
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();
cid = new char[15];sprintf(cid,"V_{B} = %d V", Vb);cid[14]=0;
TLatex *t=new TLatex(85., 9000.,cid);
t->SetTextSize(0.04); t->Draw();
TLatex *t=new TLatex(20., 8300.,"N_{eff}   [cm^{-3}]");
t->SetTextSize(0.04); t->Draw();
TLatex *t=new TLatex(25., 7600.,"1.03x10^{12} (Flat diode model)");
t->SetTextSize(0.04); t->Draw();
TLatex *t=new TLatex(25., 7000.,"8.86x10^{11} (FEM field model)");
t->SetTextSize(0.04); t->Draw();
//
TLatex *t=new TLatex(30., 2500.,"Uniform field model");
t->SetTextAngle(0); t->SetTextSize(0.04); t->Draw();
TLatex *t=new TLatex(160., 2900.,"Flat diode model");
t->SetTextAngle(18); t->SetTextSize(0.04); t->Draw();
TLatex *t=new TLatex(100., 700.,"FEM solutions");
t->SetTextColor(4); t->SetTextAngle(14); t->SetTextSize(0.04); t->Draw();
//
leg = new TLegend(0.58 ,0.87,0.81,0.92);
leg->SetHeader("distance from strip center");leg->SetBorderSize(0);
leg->SetFillColor(10);leg->SetTextSize(0.035);leg->Draw();
leg1 = new TLegend(0.58 ,0.66 ,0.71,0.87);
leg1->SetFillColor(10);leg1->SetBorderSize(0);leg1->SetTextFont(132);
leg2 = new TLegend(0.71 ,0.70 ,0.84,0.87);
leg2->SetFillColor(10);leg2->SetBorderSize(0);leg2->SetTextFont(132);
//
//
int n=115;
int m=17;
double X[m][n], Y[m][n], Potential[m][n], E[m][n],j[m][n], PDens[m][n], Ey[m][n], Ex[m][n], Jx[m][n], Jy[m][n], Epsi[m][n], Rho[m][n]; 
ifstream fin("E_65V.DAT"); //150V=Vb
for (int k=0; k<m; k++){
for (int i=0; i<n; i++){
fin >> X[k][i]>>Y[k][i]>>Potential[k][i]>>E[k][i]>>j[k][i]>>PDens[k][i]>>Ex[k][i]>>Ey[k][i]>>Jx[k][i]>>Jy[k][i]>>Epsi[k][i]>>Rho[k][i]; 
  E[k][i] = E[k][i]/100.;
  Ex[k][i] = Ex[k][i]/100.;
  Ey[k][i] = Ey[k][i]/100.;
  cout <<"k,i="<<k<<" "<<i<<" x,y="<<X[k][i] << " "<<Y[k][i]<<" "<<E[k][i]<<" Ex="<<Ex[k][i]<< " Ey=" <<Ey[k][i]<<endl; 
}
if(fmod(k,2.0)>0) continue;
int kk=k/2.0;
g1= new TGraph(n,&Y[k][0],&E[k][0]);
g1->SetMarkerStyle(7);
g1->SetMarkerColor(kk+1);
g1->Draw("SAME P");
cid = new char[9];sprintf(cid," %d #mum",kk*5);cid[8]=0; 
if (kk<5)leg1->AddEntry(g1,cid,"P"); 
if (kk>4)leg2->AddEntry(g1,cid,"P"); }
leg1->Draw(); 
leg2->Draw();
}
fin.close();
c1->Print("E_65V.pdf");
}