//-------------------------------------
//  2022.1.16  reading from  /cvmfs/atlas.cern.ch/repo/sw/database/GroupData/
//             SCT_Digitization/SCT_InducedChargedModel.root
//  2017.7.17 updated
//  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
//------------------------------------
//----------for plotting------
TCanvas *c1;
TPad *pad[10];
void setupCanvasNXNY(const int Nx, const int Ny);
const int TextFont=43, TextSize=16;
double GetNDCX(double x) ;
double GetNDCY(double y) ;
//------------------------------



#include "../Ramo_class/toyMC_InducedChargeModel.h"
#include "../Ramo_class/toyMC_InducedChargeModel.cxx"

//const std::vector<float> InducedChargeModel::s_VFD0 =
//  { -180.,  -150.,  -120.,  -90.,  -60.,  -30.,  0.,  30.,  70.};
//const std::vector<std::string> InducedChargeModel::s_VFD0str =
//  {"M180", "M150", "M120", "M90", "M60", "M30", "0", "30", "70"};
const std::vector<float> InducedChargeModel::s_VFD0 =
  {-420,-390, -360, -330, -300, -270, -240., -210., -180.,  -150.,  -120.,  -90.,  -60.,  -30.,  0.,  30.,  70.};
const std::vector<std::string> InducedChargeModel::s_VFD0str =
  {"M420","M390","M360","M330","M300","M270","M240", "M210", "M180", "M150", "M120", "M90", "M60", "M30", "0", "30", "70"};


void plot_ExEy(){
  double VFD=-400; 
  double HV= 400;

//sprintf(ROOT_FILE, "/cvmfs/atlas.cern.ch/repo/sw/database/GroupData/SCT_Digitization/SCT_InducedChargedModel.root");
sprintf(ROOT_FILE, "~/public/SCT_InducedChargedModel.root");
//------------------------------
// unit: nsec, cm, sec, K, 
setupCanvasNXNY(2,1);
c1->GetFrame()->SetBorderSize(12);
char cid[200];

//--------------------------------
  float vdepl = VFD;
  float vbias = HV;
  
  InducedChargeModel obj(0);
  InducedChargeModel::SCT_InducedChargeModelData* data{nullptr};
  
double ymin[2]={-2000., -2000.},  ymax[2]={20000., 20000.};
double EX, EY ;
double x, y;
int nhist=200;
double bulk_depth = 0.0285;
double deltay = bulk_depth/nhist;
double XX[500], YY[3][500], YFEM[9][500];
double XY[3][500], XFEM[9][500];
//-----------------------------------------
//---------------------Uniform field model-----------
  auto model = obj.UniformE;
  auto m_InducedChargeModel = std::make_unique<InducedChargeModel>(4088,model);
  m_InducedChargeModel->showEFieldModel(model);
  data = m_InducedChargeModel->setWaferData(vdepl,vbias,bulk_depth,model);
int npoint = 0;
for (double y=deltay*0.5; y<bulk_depth ; y+=deltay){
  x = 40./10000.;
  XX[npoint] = y * 10000.;
  m_InducedChargeModel->getEField(*data, x, y, EX, EY); // [V/cm]
  if(npoint%50==0) cout<<"("<<x<<","<<y<<")\tEx="<<EX<< "\tEy="<<EY<<endl;
  YY[0][npoint] = EY;
  XY[0][npoint] = EX;
  npoint++;
}

//---------------------Flat Diode Model -------------
  npoint = 0;
  model = obj.FlatDiodeModel;
  m_InducedChargeModel = std::make_unique<InducedChargeModel>(4088,model);
  m_InducedChargeModel->showEFieldModel(model);
  data = m_InducedChargeModel->setWaferData(vdepl,vbias,bulk_depth,model);
for (double y=deltay*0.5; y<bulk_depth ; y+=deltay){
  x = 40./10000.;
  XX[npoint] = y * 10000.;
  m_InducedChargeModel->getEField(*data, x, y, EX, EY); // [V/cm]
  if(npoint%50==0) cout<<"("<<x<<","<<y<<")\tEx="<<EX<< "\tEy="<<EY<<endl;
  YY[1][npoint] = EY;
  XY[1][npoint] = EX;
  npoint++;
}

//---------------------FEM field --------------------
  npoint = 0;
  model = obj.FEMsolutions;
  m_InducedChargeModel = std::make_unique<InducedChargeModel>(4088,model);
  m_InducedChargeModel->showEFieldModel(model);
  data = m_InducedChargeModel->setWaferData(vdepl,vbias,bulk_depth,model);
for (double y=deltay*0.5; y<bulk_depth ; y+=deltay){
  for (int k=0; k<9; k++){
    x = (k*5.)/10000.;
    m_InducedChargeModel->getEField(*data, x, y, EX, EY); // [V/cm]
    if(k==8 && npoint%50==0) cout<<"("<<x<<","<<y<<")\tEx="<<EX<< "\tEy="<<EY<<endl;
    YFEM[k][npoint] = EY;
    XFEM[k][npoint] = EX;
    //if(k==0)cout<<"x,y[cm]=("<<x<<","<<y<<") : Ey(0)="<<EYU<<" Ey(1)="<<EY<<" Ey(2)="<<EyFEM<<endl;
  }
  npoint++;
}

//--------------plot EX(iplot==0) Ey(iplot==1) ----------------------
double xmin=0, xmax=285.;
TGraph *hE, *hEU, *g1;
TLatex *t;

for (int iplot=0; iplot<2; iplot++) {
pad[iplot+1]->cd();
pad[iplot+1]->SetGrid();

hEU = new TGraph(npoint, XX, XY[0]) ; // uniform field model
hE = new TGraph(npoint, XX, XY[1]) ;  // flat diode model
if(iplot==1) {
  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(ymax[iplot]);
hE->SetMinimum(ymin[iplot]);
hE->SetTitle(" ");
hE->GetXaxis()->SetTitle("height y [#mum]");
hE->GetYaxis()->SetTitle("Ex [V/cm]");
if(iplot==1) hE->GetYaxis()->SetTitle("Ey [V/cm]");
hE->GetYaxis()->SetTitleOffset(2.4);
hE->GetXaxis()->SetTitleOffset(1.0);
hE->GetXaxis()->SetLimits(xmin,xmax);
hE->GetXaxis()->SetTitleFont(TextFont);
hE->GetXaxis()->SetLabelFont(TextFont);
hE->GetYaxis()->SetTitleFont(TextFont);
hE->GetYaxis()->SetLabelFont(TextFont);
hE->GetXaxis()->SetTitleSize(TextSize);
hE->GetXaxis()->SetLabelSize(TextSize);
hE->GetYaxis()->SetTitleSize(TextSize);
hE->GetYaxis()->SetLabelSize(TextSize);
  //hE ->Draw("L");  
hE->SetLineStyle(3);
hE ->Draw("AL");  
hEU->SetLineColor(2);
//hEU ->Draw("SAME L"); 
//
TLegend *leg, *leg1, *leg2;
leg1 = new TLegend(GetNDCX(30.) ,0.70,GetNDCX(100.),0.86);
leg1->SetFillColor(10);leg1->SetBorderSize(0);
leg1->SetTextFont(TextFont); leg1->SetTextSize(TextSize);
leg2 = new TLegend(GetNDCX(110.) ,0.70,GetNDCX(180.),0.86);
leg2->SetTextFont(TextFont); leg2->SetTextSize(TextSize);
leg2->SetFillColor(10);leg2->SetBorderSize(0);
//
//c1->Print("E_test.pdf");
//
for (int k=0; k<9; k++){
  g1= new TGraph(npoint,XX,XFEM[k]);
  if(iplot==1) g1= new TGraph(npoint,XX,YFEM[k]);
  g1->SetMarkerStyle(7);
  g1->SetMarkerColor(k+1);
  g1->Draw("SAME P");
  if (k<5)leg1->AddEntry(g1,Form(" %d #mum",k*5),"P");
  if (k>4 && k!=8)leg2->AddEntry(g1,Form(" %d #mum",k*5),"P");
  if (k==8)leg2->AddEntry(g1,Form(" %d (strip center)",k*5),"P");
}
leg1->Draw(); 
leg2->Draw();

sprintf(cid,"Distance from the cell edge");
t=new TLatex(30.,(ymax[iplot]-ymin[iplot])*0.94+ymin[iplot],cid); 
t->SetTextFont(TextFont);
t->SetTextSize(TextSize); t->Draw();

TLine *z = new TLine(xmin, 0., xmax, 0.);
z->Draw();

}  // iplot loop

//--------top line title------------------------
pad[0] -> cd();
sprintf(cid,"Electric field of SCT strips : FEM(colored), FDM(dotted line) for V_{FD}=%5.1fV, HV=%4.1fV", vdepl, vbias);
t=new TLatex(0.10, 0.95,cid); t->SetTextFont(TextFont);
t->SetTextSize(TextSize); t->Draw();
//

//---------------------------plot printout------------
sprintf(cid,"ExEy_%d_%d.png",int(vdepl), int(vbias));
c1->Print(cid); 
}



void setupCanvasNXNY(const int Nx, const int Ny) {
  //--------Taka Kondo (KEK) 2018.3.12-------------------
  //---- Adust following five parameters--------
  double gapx=0.100;           // gap width between plots
  double gapy=0.015;           // gap width between plots
  double left_margin=0.10;     // left_margin >= right_margin      
  double right_margin=0.03;    // left_margin >= right_margin     
  double top_margin=0.070;     // top margin      
  double bottom_margin=0.090;  // bottom margin
  //-----------------------------------------------
  c1 = new TCanvas("c1", "canvas", 700, 500);
  double apron=0.05; // pad aprons at top and bottom for label overlaps
  apron = 0.;
  double sum=0.;
  double H[9] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00};
  for (int i=0; i<Ny; i++)  { sum += H[i]; }
  for (int i=0; i<Ny; i++)  { H[i] /= sum; }
  double plotspace=1.0-(Ny-1)*gapy-top_margin-bottom_margin;
  double height, pad_topM, pad_bottomM;
   // pad[0] for full screen, pad[1] to pad[Nx] for individual plots
  double width=(1.0-(Nx-1)*gapx-left_margin-right_margin)/Nx;
  double padWidth=1.0-right_margin-(Nx-1)*width-(Nx-2)*gapx;
  if(Nx==1) padWidth=1.0;
  double x1=0.0, x2=1.0, y1=0.0, y2=1.0;
  //for (int ipad=Nx; ipad>=0; ipad--) {
  int ipad=0;
  for (int iy=0; iy<=Ny; iy++) {
  for (int ix=1; ix<=Nx; ix++) {
    if(iy==0 && ix!=1) continue;
    if(x1<0) x1=0;
    x2=1.0-right_margin-width-(Nx-ix-1)*(gapx+width);
    if(ix==Nx)x2=1.;
    if(ipad!=0 && ix==1) {
       pad_topM= gapy+apron; pad_bottomM=gapy+apron;
       if(iy==1) pad_topM=top_margin;
       if(iy==Ny) pad_bottomM=bottom_margin;
       height=pad_topM+H[iy-1]*plotspace+pad_bottomM;
       y1 = y2 - height;
       if(y1<0.01) y1=0.;
    }
    x1=x2-padWidth;
    if(x1<0.) x1=0.;
    if(ipad==0) pad[ipad]=new TPad("pad","pad",0.0,0.0,1.0,1.0);
    else pad[ipad]=new TPad("pad","pad",x1,y1,x2,y2);
    //cout<<ipad<<" "<<x1<<" "<<x2<<"    "<<y1<<" "<<y2<<endl;
    if(ipad!=0) {
      pad[ipad]->SetTopMargin(pad_topM/height);
      pad[ipad]->SetBottomMargin(pad_bottomM/height);
      pad[ipad]->SetRightMargin(gapx/padWidth);
      pad[ipad]->SetLeftMargin((padWidth-width-gapx)/padWidth);
      if(ix==Nx) y2 = y1 + gapy + 2*apron;
    }   
    if(ix==Nx){
       pad[ipad]->SetRightMargin((right_margin)/padWidth);
       pad[ipad]->SetLeftMargin((padWidth-width-right_margin)/padWidth);
    }
    //pad[ipad]->SetFillColorAlpha(ipad+1,0.10); //make pad transparent;
    pad[ipad]->SetFillColorAlpha(ipad,0.00); //make pad transparent;
    //pad[ipad]->SetFillColor(0);
    //pad[ipad]->SetFillColorAlpha(0,0.10); //make pad transparent;
    pad[ipad]->Draw();
    ipad++;
  }
  }
  return;   
}

double GetNDCX(double x) {
  //  get NDC from user-coordinate y 
  //  see https://root-forum.cern.ch/t/from-user-to-ndc/9943/4
  gPad->Update();   //this is necessary!
  return (x - gPad->GetX1())/(gPad->GetX2()-gPad->GetX1());
}

double GetNDCY(double y) {
  //  get NDC from user-coordinate y 
  //  see https://root-forum.cern.ch/t/from-user-to-ndc/9943/4
  gPad->Update();   //this is necessary!
  return (y - gPad->GetY1())/(gPad->GetY2()-gPad->GetY1());
}
  //