{
#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",300,0,700,500);
c1->SetFillStyle(4100);
//c1->SetGrid();
c1->GetFrame()->SetBorderSize(12);
double d = 0.0285;           // [cm]
double d0 = 0.0300;           // [cm]
double Vd = 70.0;            // [V] 
double Vb;                      // [V]
double Emax, wd;
int nhist=100;
double deltay=d0/nhist;
double h[nhist], EY[nhist], EYU[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 (int i=0; i<nhist; i++) {
      h[i] = deltay * i * 10000.;
      double dd = d - h[i]/10000.;
      if (Vb>Vd) {
           Emax = (Vb+Vd)/d ;
           EY[i] = Emax - 2.*Vd*dd/(d*d); //flat diode model
           EYU[i] = Vb/d; //uniform field model
      } else {
           wd = sqrt(Vb/Vd) * d;
           Emax = 2.* wd * Vd /(d*d);  
           EY[i] = Emax*(1-dd/wd);   //flat diode model
           EYU[i] = Vb/wd; //uniform field model
      }
   }
  gStyle->SetOptStat(0); 
  hE = new TGraph(nhist,h,EY);
  hEU = new TGraph(nhist,h,EYU);
  hE->SetMaximum(12000.);
  hE->SetMinimum(0.);
  hE->GetXaxis()->SetLimits(0., 285.);
  hE->SetTitle();
  hE->GetXaxis()->SetTitle("height [#mum]");
  hE->GetYaxis()->SetTitle("E_{y} [V/cm]");
  hE->GetYaxis()->SetTitleOffset(1.5);
  hE->GetXaxis()->SetTitleOffset(1.0);
  hE ->SetLineStyle(2);
  hEU ->SetLineStyle(1);
  hE ->Draw("AL");  
  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->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 field model");
t->SetTextColor(1); 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.040);leg->Draw();
leg1 = new TLegend(0.58 ,0.70 ,0.71,0.88);
leg1->SetFillColor(10);leg1->SetBorderSize(0);leg1->SetTextFont(132);
leg2 = new TLegend(0.71 ,0.75 ,0.84,0.88);
leg2->SetFillColor(10);leg2->SetBorderSize(0);leg2->SetTextFont(132);
//
//
int n=115;
//int m=17;
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,4.0)>0) continue;
int kk=k/4.0;
g1= new TGraph(n,Y,Ey);
g1->SetMarkerStyle(kk+20);
g1->SetMarkerSize(0.4);
g1->Draw("SAME P");
cid = new char[9];sprintf(cid," %d #mum", 40-kk*10);cid[8]=0; 
if (kk<3)leg1->AddEntry(g1,cid,"P"); 
if (kk>2)leg2->AddEntry(g1,cid,"P");}
leg1->Draw(); 
leg2->Draw(); 
c1->Print("E_150V_1D.eps");
}