//========================================================================
// simulation of hole and electron trajectories of  ATLAS SCT barrel strips
// 2017.7.17 updated
// 2016.3.30 combine with induced_charge.h
// 2016.2.15 updated for ROOT 6.02/12
// 2010.1.20 T. Kondo(KEK)
// y = 0 for HV plane, y = d for strip plane
// with library of sub_functions.cxx
//========================================================================
#include "Ramo.h"

void plot_eh_trajectory() {
//------------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
for (m_EfieldModel = 0; m_EfieldModel<3; m_EfieldModel++) {
m_VD = 70.;
m_VB = 150.;
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=2000;
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, ymax;
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;
//
TH1F *hdummy = new TH1F("hdummy","LorentzAngle",400, -160., 240.); 
hdummy->SetTitle("test");
hdummy->GetXaxis()->SetTitle("X [#mum]");
hdummy->GetXaxis()->SetTitleOffset(1.1);
hdummy->GetYaxis()->SetTitle("Y [#mum]");
hdummy->GetYaxis()->SetTitleOffset(1.1);
hdummy->SetMaximum(285.0);
hdummy->SetMinimum(0.0);
hdummy->Draw();
//---------------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------------------------------
TLine *z1;
for (int k=0; k<2; k++){
   x_origin = -0.0080;
   theta_0 = 0.;
   ntry = 15;
   if (k==1){
     x_origin = 0.000;
     theta_0 = 30.;
     ntry = 70;
   }
   theta_0 = k*30.;
   tantheta_0 = tan (theta_0/radian);
   z1 = new TLine(x_origin*10000.,0., x_origin*10000.+285. * tantheta_0, 285.); 
   z1->SetLineColor(1);   z1->SetLineWidth(1);   z1->Draw();

//---------loop over depth distribution -----------------------
efactor = e/ntry ;
for (int i=0; i<ntry ; i++) {
   if(fmod(i,100)==0) cout<<"i="<<i<<endl;
   y_origin = 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<50. ; 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++ ;
        X[j] = x * 10000.;
        Y[j] = y * 10000.;
        //cout << "time="<<time<<" j="<<j<<" (x,y)=("<<x<<" "<<y<<")"<<endl;
        if ( y > m_bulk_depth ) ipass = 1;
   }
gStyle->SetOptStat(0);
TGraph *h2 = new TGraph(j+1, X, Y);
h2->SetLineColor(fmod(i,9.)+1);
h2 -> Draw("SAME L");
//-------------------trace electrons------------------------------ 
   j = 0;
   x = x_origin + y_origin * tantheta_0;
   y = y_origin;
   X[j] = x * 10000.;
   Y[j] = y * 10000.;
   ipass = 0;
   for (time=0. ; time<50. ; time+=m_transportTimeStep ) {
      if ( ipass ==1 ) continue;
      if (!electron(x, y, vx, vy, D)) continue ;
      delta_y = vy * m_transportTimeStep *1.E-9; 
      y += delta_y;
      double dt = m_transportTimeStep;
      if ( y < y_origin_min ) {
          ipass = 1 ;
          dt = ( y_origin_min - (y-delta_y) )/delta_y * m_transportTimeStep; 
          y = y_origin_min;
      }
      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 < y_origin_min ) ipass = 1;
      j++ ;
      X[j] = x * 10000.;
      Y[j] = y * 10000.;
      //cout << "time="<<time<<" j="<<j<<" (x,y)=("<<x<<" "<<y<<")"<<endl;
   }
//
gStyle->SetOptStat(0);
h2 = new TGraph(j+1, X, Y);
h2->SetLineColor(fmod(i,9.)+1);
h2->SetLineStyle(2);
h2 -> Draw("SAME L");
   z1 = new TLine(x_origin*10000.,0., x_origin*10000.+285. * tantheta_0, 285.); 
   z1->SetLineColor(1);   z1->SetLineWidth(1);   z1->Draw();
}
}
//
if(m_EfieldModel==2) {
z1=new TLine(-130.,286, -110., 286.); z1->SetLineColor(1); 
z1->SetLineWidth(4);z1->Draw();
z1=new TLine(-50.,286, -30., 286.); z1->SetLineColor(1); 
z1->SetLineWidth(4);z1->Draw();
z1=new TLine(30.,286, 50, 286.); z1->SetLineColor(1); 
z1->SetLineWidth(4);z1->Draw();
z1=new TLine(110.,286, 130, 286.); z1->SetLineColor(1); 
z1->SetLineWidth(4);z1->Draw();
z1=new TLine(190.,286, 210, 286.); z1->SetLineColor(1); 
z1->SetLineWidth(4);z1->Draw();
} else {
z1=new TLine(-160.,286, -121, 286.);
z1->SetLineColor(2);z1->SetLineWidth(4);z1->Draw();
z1=new TLine(-119.,286, -41., 286.); 
z1->SetLineColor(1); z1->SetLineWidth(4);z1->Draw();
z1=new TLine(-39.,286, 39, 286.); 
z1->SetLineColor(2);z1->SetLineWidth(4);z1->Draw();
z1=new TLine( 41.,286, 119, 286.); 
z1->SetLineColor(1); z1->SetLineWidth(4);z1->Draw();
z1=new TLine(121.,286, 199, 286.);
z1->SetLineColor(2);z1->SetLineWidth(4);z1->Draw();
z1=new TLine(201.,286, 240, 286.);
z1->SetLineColor(1);z1->SetLineWidth(4);z1->Draw();
}
//
//
TLatex *t;
double xx = 90.;
TPave *p = new TPave( xx-4.,20., 233., 90.); p->SetLineColor(1);
p->SetFillColor(10);p->SetBorderSize(1);p->Draw(); 
char model[3][30]={"Uniform field model","Flat diode model","FEM field model"};
t=new TLatex( xx, 70., model[m_EfieldModel]);
t->SetTextSize(0.04); t->Draw();
sprintf(m_cid,"V_{B}=%4.0f V, V_{D}=%4.0f V",m_VB,m_VD);
t=new TLatex( xx, 50.,m_cid);
t->SetTextSize(0.04); t->Draw();
sprintf(m_cid,"B=%4.1f Tesla, T=%4.0f K",m_B,m_T);
t=new TLatex( xx, 30.,m_cid);
t->SetTextSize(0.04); t->Draw();
t=new TLatex( -60, 150.,"0^{o} track");
t->SetTextAngle(90);
t->SetTextSize(0.04); t->Draw();
t=new TLatex( 170, 200.,"30^{o} track");
t->SetTextAngle(60); 
t->SetTextSize(0.04); t->Draw();
//
sprintf(m_cid,"eh_trajectory_%d.png",m_EfieldModel);
c1->Print(m_cid);
//--------------------------------
}   // end of m_EfieldModel loop
return;
}