//========================================================================
// 2019.3.6  Taka Kondo (KEK)
//========================================================================
#include "Ramo.h"

void plot_2D() {
//------------set main parameters and initialization----------------
m_EfieldModel = 0;  // Uniform Electric Field model;
m_EfieldModel = 1;  // Flat Diode Model;
m_EfieldModel = 2;  // Electric field by FEM solusion
m_VD = 70.;
m_VB = 150.;
m_transportTimeStep = 0.01;    // one step side in time [nsec]
initialize();

//--------------------------------------------
TCanvas *c1 = new TCanvas("c1","A Simple Graph",0,0,700,500);
c1->SetFillStyle(4100);
//c1->SetLeftMargin(0.12); 
//c1->SetBottomMargin(0.12); 
int istrip, ipass ; 
int    ntry = 2, ipoints, npoints=900;
double X[npoints], Y[npoints];
double strip_width=0.008, pi = 3.141592 ;
double x, delta_x, y, delta_y, q, delta_q, vx, vy, D, sigma;
double x_origin, y_origin, y_origin_min,  theta_0;
double tantheta_0, radian = 180./pi;
double e = m_bulk_depth*1.E4*72.*1.6E-19/1.E-15, efactor ;  // [FC] 
double time ; // [ns]
double delta_y_new, LorentzAngle; 
double diffusion;
//
//---------------setting basic parameters---------------------------
bool isDiffusion = true;
x_origin = 0.000;
y_origin_min = 0.0000;

// ----------- find for the case of model = 2----------------------
double Ex, Ey;
if (m_EfieldModel == 2) {
  for ( double y = 0.; y < m_depletion_depth ; y += 0.00025 ) {
     EField(0., y, Ex, Ey);
     if ( Ey > 0.1) continue;
     y_origin_min = y ;
  }
} else y_origin_min = m_bulk_depth - m_depletion_depth ; 
std::cout << "model, y_origin_min = " <<m_EfieldModel<<"    "
<< y_origin_min * 1.E4 <<" um"<< std::endl;

//--------particle trajectory-  k=0 0-degree, k=1 30 degree--------
TLine *z1;
x_origin = 0.000;
theta_0 = 30.;
tantheta_0 = tan (theta_0/radian);
ntry = 500;
ntry = 500000;

//----------setup histograms
double tmax = 3.1, checktime[5] = {0.1, 0.5, 1.0, 2.0, 3.0};
double xmin=40., xmax=200., ymin=200., ymax=285.;
TH2F *h2D[5];
for (int i=0;i<5;i++) 
  h2D[i] = new TH2F(Form("%d",i),"2D;X [#mum];Y [#mum]",200,xmin,xmax,200,ymin,ymax); 
//---------loop over depth distribution -----------------------
efactor = e/ntry ;
for (int i=0; i<ntry ; i++) {
   if(fmod(i,1000)==0) cout<<"i="<<i<<endl;
   y_origin = gRandom->Rndm()*m_bulk_depth;
   //y_origin = (0.6+0.4*gRandom->Rndm())*m_bulk_depth;
   if(y_origin <= y_origin_min) continue;

//-------------------trace a hole------------------------------ 
   int j = 0;
   x = x_origin + y_origin * tantheta_0;
   y = y_origin;
   X[j] = x * 10000.;
   Y[j] = y * 10000.;
   TMarker *m=new TMarker(x*10000.,y*10000.,8); 
   m->SetMarkerSize(0.75); m->SetMarkerColor(fmod(i,9.)+1);m->Draw(); 
   ipass = 0;
   //for (time=0. ; time<50. ; time+=m_transportTimeStep ) {
   for (time=0. ; time<=tmax ; time+=m_transportTimeStep ) {
        if ( ipass ==1 ) continue;
        if (!hole(x, y, vx, vy, D)) continue ;
        delta_y = vy * m_transportTimeStep *1.E-9; 
        y += delta_y;
        double dt = m_transportTimeStep;
        if ( y > m_bulk_depth ) {
            ipass = 1 ;
            dt = (m_bulk_depth - (y-delta_y))/delta_y * m_transportTimeStep; 
            y = m_bulk_depth;
        }
        x += vx * dt *1.E-9;
        if (isDiffusion) {
           diffusion = sqrt (D * dt*1.E-9);
           y += diffusion * gRandom->Gaus(0, 1.);
           x += diffusion * gRandom->Gaus(0, 1.);
        }
        if( y > m_bulk_depth) ipass = 1;
        j++ ;
        for (int ktime=0; ktime<5; ktime++) {
           if( abs(time - checktime[ktime]) < m_transportTimeStep*0.2) {
           h2D[ktime] ->Fill(x*10000., y*10000.);
           //cout<<i<<" "<<time<<" "<<x*10000.<<" "<<y*10000.<<endl;
           }
        }
        //cout << "time="<<time<<" j="<<j<<" (x,y)=("<<x<<" "<<y<<")"<<endl;
        if ( y > m_bulk_depth ) ipass = 1;
   }
}
//
//------------------------------------------
h2D[0]->Draw("COL");
for (int i=1; i<5; i++) h2D[i]->Draw("SAME COL");
// draw track trajectory
z1 = new TLine(x_origin*10000.+ymin*tantheta_0,ymin, 
               x_origin*10000.+ymax * tantheta_0, ymax); 
z1->SetLineColor(1);   z1->SetLineWidth(1);   z1->Draw();
for (int i=0; i<5; i++) {
  TLatex *t=new TLatex((xmax-xmin)*0.8+xmin,(ymax-ymin)*(0.5-0.07*i)+ymin,
    Form("%5.2f ns",checktime[i]));
  t->Draw();
}

TBox *b = new TBox(40.,285.,50.,286.); b->SetFillColor(1); b->Draw(); 
b = new TBox(110.,285.,130.,286.); b->SetFillColor(1); b->Draw(); 
b = new TBox(190.,285.,200.,286.); b->SetFillColor(1); b->Draw(); 

sprintf(m_cid,"30degree_track.png");
c1->Print(m_cid);
//--------------------------------
return;
}