//-----------------------------------------------
//   2017.7.18 updated
//   2016.3.31, 4.1  for ICM (induced charge model) by Taka Kondo (KEK)
//-----------------------------------------------
#include "Ramo.h"

void plot_pulse(){
  setupCanvas(5,1);
//------------set main parameters and initialization----------------
 m_EfieldModel = 1;
 m_VD = 70.;
 m_VB = 150.;
 m_isSCTDigiModel = false; // induced charge model
 m_isSCTDigiModel = true;  // (old) SCT Digitization model
 initialize();
//------------------------------------
  m_coutLevel = 3;
  for (int n=0; n<9; n++) {
    m_x0 = int(n/3.) * 40.;  // ionization point in x [micron]
    m_y0 = n%3 * 130.;       // ionization point in y[micron]
    m_eventNumber++ ;
    execute();
  }
//----------plot charges -----------------
  if(m_isSCTDigiModel) canvas->Print("SCTDigi_pulse.png");
  else canvas->Print("Ramo_pulse.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 );
}
//-----------put through amplifiers associated with crosstalks-------------
//  storage arrays for every 0.5 nsec for amplified pulses
double Pulse[5][200], nsec[200];
double delta_q;
for (int it=0; it<200; it++) {
  nsec[it] = (it + 0.5) * 0.5;
  for (int istrip=0; istrip<5; istrip++) {
    Pulse[istrip][it]=0.;
} }
for (int istrip = 0 ; istrip < 5 ; istrip++ ){
  for (int it = 0 ; it < 50 ; it++ ) {
    switch(istrip){
      case 0: delta_q = Q_m2[it];  break;
      case 1: delta_q = Q_m1[it];  break;
      case 2: delta_q = Q_00[it];  break;
      case 3: delta_q = Q_p1[it];  break;
      case 4: delta_q = Q_p2[it];  break;
    }
    if (delta_q == 0.) continue;
    double time0 = (it + 0.5) * 0.5;     // timing of current input
    for (int jt = 1 ; jt < 200 ; jt++ ) {   // trace every 0.5ns step
      double timeAmp = jt * 0.5  ;
      double time = time0 + timeAmp;
      int itime = time / 0.5;
      if( itime >= 200) break;
      double response =  Amp_response(timeAmp)  * delta_q;
      double crosstalk = Amp_crosstalk(timeAmp) * delta_q;
      Pulse[istrip][itime] += response;
      if( istrip > 0) Pulse[istrip-1][itime] += crosstalk;
      if( istrip < 3) Pulse[istrip+1][itime] += crosstalk;
    } 
  }  
}
//-------------------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 = 1.2, 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(200, nsec, Pulse[0]); 
            ylabel = 0.15; ytitle = 0.13; break;
    case 1: h0 = new TGraph(200, nsec, Pulse[1]); break;
    case 2: h0 = new TGraph(200, nsec, Pulse[2]); break;
    case 3: h0 = new TGraph(200, nsec, Pulse[3]); break;
    case 4: h0 = new TGraph(200, nsec, Pulse[4]); 
            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.,95.);
  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", 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", 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++;
  }}
}