{
#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 = 70.0;            // [V] 
double Vb;                      // [V]
double EY, EYU, 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=150.;
   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(12000.);
  hE->SetMinimum(0.);
  hE->SetTitle();
  hE->SetXTitle("height [#mum]");
  hE->SetYTitle("E_{y} [V/cm]");
  hE->GetYaxis()->SetTitleOffset(1.5);
  hE->GetXaxis()->SetTitleOffset(1.0);
  hE ->Draw("L");  
  hEU->SetLineColor(2);
  hEU ->Draw("SAME L"); 
//
b = new TPave(15.,8600, 145., 11600.); 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(20., 11000.,cid);
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., 11000.,cid);
t->SetTextSize(0.04); t->Draw();
TLatex *t=new TLatex(20., 10300.,"N_{eff}   [ cm^{-3}] ");
t->SetTextSize(0.04); t->Draw();
TLatex *t=new TLatex(25., 9500.,"1.1x10^{12} (Flat diode model)");
t->SetTextSize(0.04); t->Draw();
TLatex *t=new TLatex(25., 8900.,"9.6x10^{11} (FEM field model)");
t->SetTextSize(0.04); t->Draw();
//
TLatex *t=new TLatex(30., 5500.,"Uniform field model");
t->SetTextAngle(0); t->SetTextColor(2); t->SetTextSize(0.04); t->Draw();
TLatex *t=new TLatex(20., 3500.,"Flat diode model");
t->SetTextAngle(15); t->SetTextSize(0.04); t->Draw();
TLatex *t=new TLatex(110.,3500.,"FEM solutions");
t->SetTextColor(4); t->SetTextAngle(14); t->SetTextSize(0.04); t->Draw();
//
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);
//
//

int n=115;
int m=17;
double Y[n], Ey[n], Ex[n]; 

ifstream fin("E_150V.DAT"); //150V=Vb
for (int k=0; k<m; k++){
 double x = k * 2.5;
 for (int i=0; i<n; i++){
   Y[i] = i * 2.5;
   EField( x/10000., Y[i]/10000.,Ex[i], Ey[i] );
 }
if(fmod(k,2.0)>0) continue;
int kk=k/2.0;
g1= new TGraph(n,Y,Ey);
g1->SetMarkerStyle(7);
g1->SetMarkerColor(kk+1);
g1->Draw("SAME P");
cid = new char[9];sprintf(cid," %d #mum", 40-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();
c1->Print("E_150V_1D.png");
}