//========================================================================
// simulation of hole and electron trajectories of  ATLAS SCT barrel strips
// 2017.7.17 updated
// 2016.10.14 combine with induced_charge.h updated for ROOT 6.02/12
// 2010.1.20 T. Kondo(KEK)
//========================================================================

//#include "induced_charge.h"
#include "Ramo.h"

void plot_charge(){
  setupCanvas(5,1);
//------------set main parameters and initialization----------------
  m_EfieldModel = 2; // use 2 for Ramo model
  m_VD = 70.;
  m_VB = 150.;
  m_isSCTDigiModel = false ;
  m_isSCTDigiModel = true ;
//-------------------------
  initialize();
//-------------------------
  m_x0 = 40.0 ;  // ionization point [micron]
  m_y0 = 000.0 ; // ionization point [micron]
  m_coutLevel = 3;
  for (int n=0; n<9; n++) {
    m_x0 = int(n/3.) * 40.;
    m_y0 = n%3 * 130.;
    m_eventNumber++ ;
    execute();
  }
  if(m_isSCTDigiModel) canvas->Print("SCTDigi_charge.png");
  else canvas->Print("Ramo_charge.png");
return;
}


void execute() {

//  storage arrays for every 0.5 nsec
double Q_m2[50], Q_m1[50], Q_00[50], Q_p1[50], Q_p2[50] ;
double x0 = m_x0 / 10000. ; // conversion from [micron] to [cm] 
double y0 = m_y0 / 10000. ; // conversion from [micron] to [cm] 

//------------- one event loop ---------------------------------
//  clear charge arrays for induced current
for (int it=0; it<50; it++) {
    Q_m2[it]=0.; Q_m1[it]=0.; Q_00[it]=0.; Q_p1[it]=0.; Q_p2[it]=0.;
}
//  perform electron and hole transportation starting x0 and y0
if( !m_isSCTDigiModel) {
   holeTransport    (x0, y0, Q_m2, Q_m1, Q_00, Q_p1, Q_p2) ;
   electronTransport(x0, y0, Q_m2, Q_m1, Q_00, Q_p1, Q_p2) ;
} else {
   cout<<"old model"<<endl;
   oldSCTDigitization(x0, y0, Q_m2, Q_m1, Q_00, Q_p1, Q_p2 );
}
//  printout one event 
double nsec[50];
std::cout<<"-------------------------------------------"<<std::endl;
std::cout<<"EventNumber = "<<m_eventNumber<<", (x0,y0)=("
  <<x0*1e4<<","<<y0*1.e4<<") [microns]"<<std::endl;
std::cout<< std::setprecision(5)<<std::setiosflags(std::ios::fixed)
  <<"   time :         Q_m2      Q_m1      Q_00      Q_p1      Q_p2"<<std::endl;
double qm2=0., qm1=0., q00=0., qp1=0., qp2=0.;
for (int it = 0; it< 50; it++) {
   double time = (it +1) * 0.5 ;
   nsec[it] = time;
   std::cout <<"  "<<time<<" ns : "
   <<std::setw(10)<<Q_m2[it]<<std::setw(10)<<Q_m1[it] 
   <<std::setw(10)<<Q_00[it]<<std::setw(10)<<Q_p1[it] 
   <<std::setw(10)<<Q_p2[it]<<std::endl;
   qm2 += Q_m2[it];
   qm1 += Q_m1[it];
   q00 += Q_00[it];
   qp1 += Q_p1[it];
   qp2 += Q_p2[it];
}

double qt = qm2 + qm1 + q00 + qp1 + qp2;
std::cout <<"  TOTAL SUM   : "
<<std::setw(10)<<qm2<<std::setw(10)<<qm1 
<<std::setw(10)<<q00<<std::setw(10)<<qp1 
<<std::setw(10)<<qp2<<std::setw(10)<<qt<<std::endl;
//-------------------plot----------------------
TGraph *h0;
TLegend *leg;
TLatex *txt;
if(m_eventNumber==1) leg = new TLegend(0.03,0.5,0.93,0.90,"eh pair origin [#mum]");
double ymax = 0.6, ymin=-0.2;
if(m_isSCTDigiModel) ymax=1.2;
for (int k=0; k<5; k++) {
  canv[k]->cd();
  double ylabel = 0.0, xlabel = 0.15;
  double ytitle = 0.0, xtitle = 0.00;
  switch(k) {
    case 0: h0 = new TGraph(50, nsec, Q_m2); 
            ylabel = 0.15; ytitle = 0.13; break;
    case 1: h0 = new TGraph(50, nsec, Q_m1); break;
    case 2: h0 = new TGraph(50, nsec, Q_00); break;
    case 3: h0 = new TGraph(50, nsec, Q_p1); break;
    case 4: h0 = new TGraph(50, nsec, Q_p2); 
            xtitle = 0.15; break;
  }
  h0->SetMarkerSize(0.7);
  h0->SetMarkerColor(m_eventNumber);
  h0->SetLineColor(m_eventNumber);
  h0->GetXaxis()->SetLabelSize(xlabel);
  h0->GetXaxis()->SetLabelOffset(-0.07),
  h0->GetXaxis()->SetTitle("time [nsec]");
  h0->GetXaxis()->SetTitleSize(xtitle);
  h0->GetXaxis()->SetTitleOffset(0.35),
  h0->GetYaxis()->SetLabelSize(ylabel);
  h0->GetYaxis()->SetLabelOffset(0.1);
  h0->GetYaxis()->SetTitle("Induced charge [electrons]");
  h0->GetYaxis()->SetTitleSize(ytitle);
  h0->GetYaxis()->SetTitleOffset(1.95),
  h0->SetMaximum(ymax); h0->SetMinimum(ymin);
  //h0->GetXaxis()->SetLimits(0., 32.);
  h0->GetXaxis()->SetLimits(0., 18.);
  h0->GetXaxis()->SetNdivisions(505);
  if(m_eventNumber==1) h0->Draw("ALP");
  else h0->Draw("SAME LP");
//------------write conditions at the top 
  if(m_eventNumber==1 && k ==0) {
     sprintf(m_cid,"Ramo model : B=%5.1fTesla, V_{D}=%4.1fV, V_{B}=%5.1fV, E-field and Potential from FEM", m_B, m_VD, m_VB);
     if( m_isSCTDigiModel ) 
     sprintf(m_cid,"SCT Digitization Model : B=%5.1fTesla, V_{D}=%4.1fV, V_{B}=%5.1fV, Efield:Flat Diode Model", m_B, m_VD, m_VB);
     txt = new TLatex(0., (ymax-ymin)*1.025+ymin, m_cid);
     txt->SetTextSize(0.13); txt-> Draw();
  }
//-----------write strip number----------
  if(m_eventNumber==1) {
     sprintf(m_cid,"strip %d", k-2);
     txt = new TLatex(5., (ymax-ymin)*0.94+ymin, m_cid);
     txt->SetTextSize(0.15); txt-> Draw();
  }
//---------write legends----------------------
  if (k !=4 ) continue;  
  sprintf(m_cid,"x,y=%4.0f,%4.0f",m_x0,m_y0);
  leg -> AddEntry(h0,m_cid,"lp"); 
  if(m_eventNumber != 9) continue;
  leg->SetTextSize(0.10); leg->Draw();
}
//------------- end of one event loop -------------------------
 return;
}


void setupCanvas(const int NX, const int NY) {
//-------------set up for NX x NY plots--------------------
//TCanvas canvas("canvas");
  canvas = new TCanvas("canvas","12 plots",0,0,1400,1000);
  canvas->Divide(NX, NY);
  double rightM = 0.02, leftM = 0.10, topM = 0.07, bottomM=0.08;
  double gapY=0.03, gapX=0.02;
  double height[NY], width[NX];
  for(int i=0;i<NX;i++){width[i]=1./NX;}
  for(int i=0;i<NY;i++){height[i]=1./NY;}
  double deltaW = 1.-leftM-rightM;
  double deltaH = 1.-topM-bottomM;
  double X1[NX], X2[NX], Y1[NY], Y2[NY]; 
  for (int i=0; i<NX; i++) {
    if(i==0) X1[i] = leftM;
    else X1[i] = X2[i-1];
    X2[i] = X1[i]+deltaW*width[i];
    cout<<"SetupCanvas_x: i="<<i<<" width="<<width[i]<<" x1="<<X1[i]<<" x2="<<X2[i]<<endl;
    }
  for (int i=0; i<NY; i++) {
    if(i==0) Y2[i] = 1. - topM;
    else Y2[i] = Y1[i-1];
    Y1[i] = Y2[i]-deltaH*height[i];
    cout<<"SetupCanvas_y: i="<<i<<" height="<<height[i]<<" y1="<<Y1[i]<<" y2="<<Y2[i]<<endl;
  }
//----------------------------------------
  int ixy = 0;  
  for (int iy=0;iy<NY;iy++) {  
  for (int ix=0;ix<NX;ix++) {
      sprintf(m_cid,"canvas_%d",ixy+1);
      canv[ixy] = (TPad*) canvas->GetListOfPrimitives()->FindObject(m_cid);
      canv[ixy]->SetPad(X1[ix],Y1[iy],X2[ix],Y2[iy]);
      canv[ixy]->SetFillColor(0);
      canv[ixy]->SetFrameFillColor(0);
      canv[ixy]->SetTopMargin(0.);
      canv[ixy]->SetBottomMargin(gapY);
      canv[ixy]->SetLeftMargin(0.);
      canv[ixy]->SetRightMargin(gapX);
      ixy++;
  }}
}