//-----------------------------------------------
//   2017.7.18 updated
//   2016.4.1  for ICM (induced charge model) by Taka Kondo (KEK)
//-----------------------------------------------
#include "Ramo.h"
TProfile *m_h_phi ;
TLegend *leg; 
TLatex *t;

void plot_track(){
  setupCanvas(1,1);
//------------set main parameters and initialization----------------
  m_EfieldModel = 2;
  m_VD = 70.;
  m_VB = 150.;
//-----------------
  m_coutLevel = 2;  // 2(print origin of eh pair and hole/elec results)
  m_coutLevel = 1;  // 0(print ievery 2000 event)
  m_coutLevel = 3;  // 3(print t,x,y of all time steps)
  m_coutLevel = 0;  // 0(print ievery 2000 event)
//-------------------------------
  m_phi_min = -16.;
  m_phi_max = 8. ;
  m_eta_min = -1.5 ;
  m_eta_max = 1.5 ;
  int TotalEventNumber = 500000;
  double ymin= 0.9, ymax=2.0;
  leg = new TLegend(0.7,0.8,0.9,0.9);
  m_isNewCrossTalk = false;
  m_isNewCrossTalk = true;
//---------digitization models---------------
for (int model=0; model<2; model++) {
  m_isSCTDigiModel = false; // induced charge model
  if( model==0) m_isSCTDigiModel = true;  // (old) SCT Digitization model
  initialize();
//-----------------
  m_h_phi = new TProfile("phi_dist","phi_distrib",48, -16., 8., 0., 20.); 
//------------------------------------
//------------local_phi and eta samplings -------------
  for (int n=0; n<TotalEventNumber ; n++)  execute();
//--------------------------
//----------plots -----------------
  m_h_phi -> GetXaxis()->SetLabelSize(0.04);
  m_h_phi -> GetYaxis()->SetLabelSize(0.04);
  m_h_phi -> GetXaxis()->SetTitleSize(0.04);
  m_h_phi -> GetYaxis()->SetTitleSize(0.04);
  m_h_phi -> GetXaxis()->SetTitle("Incident angle [#circ]"); 
  m_h_phi -> GetXaxis()->SetTitleOffset(1.0); 
  m_h_phi -> GetYaxis()->SetTitle("New Cluster Width"); 
  m_h_phi -> GetYaxis()->SetTitleOffset(1.1); 
  m_h_phi -> GetYaxis()->SetLabelOffset(0.01); 
  m_h_phi -> SetMaximum(ymax);
  m_h_phi -> SetMinimum(ymin);
  m_h_phi -> SetMarkerColor(model+1);
  if(model==0) m_h_phi -> Draw();
  else m_h_phi -> Draw("SAME P");
//----------------------------------
  sprintf(m_cid,"Ramo model");
  if(m_isSCTDigiModel) sprintf(m_cid,"SCTDigi model");
  leg->AddEntry(m_h_phi,m_cid,"p"); 
}
  leg->SetFillColor(10);leg->SetBorderSize(0);
  leg->SetTextSize(0.04); leg->Draw(); 
  sprintf(m_cid,"V_{D}=%5.1fV, V_{B}=%5.1fV",m_VD, m_VB);
  t=new TLatex(-15.0,(ymax-ymin)*0.9+ymin,m_cid); 
  t->SetTextSize(0.04);   t->Draw(); 
  sprintf(m_cid,"Eta ranges : %4.1f - %4.1f ",m_eta_min, m_eta_max);
  t=new TLatex(-15.0,(ymax-ymin)*0.84+ymin,m_cid); 
  t->SetTextSize(0.04);   t->Draw(); 
  sprintf(m_cid,"with original cross-talk shape");
  if(m_isNewCrossTalk) sprintf(m_cid,"with new cross-talk shape");
  t=new TLatex(-15.0,(ymax-ymin)*0.78+ymin,m_cid); 
  t->SetTextSize(0.04);   t->Draw(); 
//-------------output plot-----
  sprintf(m_cid,"Lorentz_VD%d_VB%d.png",int(m_VD), int(m_VB));
  canvas->Print(m_cid);
//----------------------------------------
return;
}

void execute() {

    m_eventNumber++ ;

//--------------- dE/dX fluctuation ------------------------------
double Landau = 108.;
if(m_isLandau) {
  for (int j=0; j<100; j++) {
      Landau = gRandom->Landau(72., 11.8) ;
     if (Landau < 500.) break;
  }
}
       // eh_pairs = number of eh_pairs per micros
double eh_pairs = m_eh_pairs * Landau/108. ;
//------------------particle track generation------------------------------
double radian = 180./3.141593;
double width = m_phi_max - m_phi_min;
double phi = m_phi_min + width * gRandom->Uniform(0.0,1.0);
       width = m_eta_max - m_eta_min;
double eta = m_eta_min + width * gRandom->Uniform(0.0,1.0);
double z_theta = 2.* atan(exp(-eta));
       width = m_x_origin_max - m_x_origin_min;
double x_origin = m_x_origin_min + width * gRandom->Uniform(0.0,1.0);
       //  charge = total eh_pair charges in  fC
double charge = m_bulk_depth * 1.E4 * eh_pairs *1.6E-19/1.E-15;  
//cout<<"eh_pairs = "<<eh_pairs<<" /micron, charge = "<<charge<<" fC"<<endl;
double qfactor = ( charge / cos(phi/radian) / sin(z_theta) ) /m_ydivision ;
//-------- first 50 events ----------------------------
if(m_eventNumber < 10 ) 
    std::cout << "ev#="<<m_eventNumber<<" eh="<< eh_pairs
    << " phi="<<phi<<" eta="<<eta<<" x0="<<x_origin*1.e4
    << "um,  qfactor="<<qfactor<<std::endl;
//----- accumulate charges for strips of -10,-9,....,+10 -------
double qstrip[21][50] = {{}};
//---------loop over depth distribution -----------------------
for (int ievent=0; ievent<m_ydivision ; ievent++) {    
   //double y_origin = m_bulk_depth * gRandom->Uniform(0.0,1.0);
   double y_origin = m_bulk_depth * (ievent+0.5)/m_ydivision;
   if (y_origin <= m_bulk_depth - m_depletion_depth ) continue;
//-------------------get origin of e-h pairs------------------------ 
   double x_eh = x_origin + y_origin * tan (phi/radian);
   double y_eh = y_origin;
   if(m_coutLevel>1) 
     std::cout<<"ievent,(x_eh,y_eh)="<<ievent
     <<" ("<<x_eh*1.e4<<","<<y_eh*1.e4<<")"<<std::endl;
//----find strip number containing x_ehpair---------------------
   int strip_offset = int((x_eh + m_strip_pitch*9999)/m_strip_pitch) - 9999;
   double x0 = x_eh - strip_offset * m_strip_pitch;
   double y0 = y_eh;
//------------- one event loop ---------------------------------
//  temporal storage arrays for every 0.5 nsec for 5 strips nearby
    double Q_m2[50]={}, Q_m1[50]={}, Q_00[50]={}, Q_p1[50]={}, Q_p2[50]={} ;
//  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 oldSCTDigitization(x0, y0, Q_m2, Q_m1, Q_00, Q_p1, Q_p2 );
//
    for (int it=0; it<50; it++) {
//  m2
       int jstrip = 10 + strip_offset - 2;
       if (jstrip >=0 && jstrip < 21) qstrip[jstrip][it] += Q_m2[it];
//  m1
       jstrip++;
       if (jstrip >=0 && jstrip < 21) qstrip[jstrip][it] += Q_m1[it];
//  00
       jstrip++;
       if (jstrip >=0 && jstrip < 21) qstrip[jstrip][it] += Q_00[it];
//  p1
       jstrip++;
       if (jstrip >=0 && jstrip < 21) qstrip[jstrip][it] += Q_p1[it];
//  p2
       jstrip++;
       if (jstrip >=0 && jstrip < 21) qstrip[jstrip][it] += Q_p2[it];
    } // --- end of time bin loop

}  // -----------------end of one track loop ---------------

//-----------------print one track result----------------
if( m_eventNumber < 3 ) {
  double totalCharge = 0.;
  for (int it=0; it<50; it++) {
    sprintf(m_cid,"%d : %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f",it,
    qstrip[7][it],qstrip[8][it],qstrip[9][it],qstrip[10][it],qstrip[11][it],
    qstrip[12][it],qstrip[13][it],qstrip[14][it]);
    cout<<m_cid<<endl;
  for (int istrip=0; istrip<21; istrip++) totalCharge += qstrip[istrip][it];
  } 
  cout<<"Total induced charge = "<<totalCharge*qfactor<<" fC"<<endl;
}
//-------------------amplify the charges ----------------------
double pulse[21][3];   // pulse height at m_timeOfThreshold -25,0,+25 ns
for (int i=0 ; i<21; i++) { 
    for(int j=0; j<3 ;j++) pulse[i][j] = 0.;
}
double timeOfThreshold[3];
for (int i=0 ; i<3 ;i++) timeOfThreshold[i] = m_timeOfThreshold + (i-1)*25.0;

for (int istrip = 0; istrip < 21 ; istrip++) {
   for (int it = 0 ; it < 50 ; it++) {
      double delta_q = qstrip[istrip][it] * qfactor;
      if (delta_q == 0.) continue;
      double time0 = (it + 0.5) * 0.5; 
      for (int j=0; j<3 ;j++) {
         double timeAmp = timeOfThreshold[j] - time0;
         if ( timeAmp <= 0.) continue;
         double response =  Amp_response(timeAmp)  * delta_q;
         double crosstalk = Amp_crosstalk(timeAmp) * delta_q;
         pulse[istrip][j] += response;
         if( istrip > 0) pulse[istrip-1][j] += crosstalk;
         if( istrip < 19) pulse[istrip+1][j] += crosstalk;
       }
    } 
}
//-------------------------------
if(m_coutLevel>1) {
   std::cout<<"tbin\ttime\tpulse[8]\tpulse[9]\tpulse[10]\tpulse[11]"
   <<std::endl;
   for (int i=0;i<3;i++){
      std::cout<<i<<"\t"<<timeOfThreshold[i] <<"\t"
      <<pulse[8][i]<<"\t" <<pulse[9][i]<<"\t"
      <<pulse[10][i]<<"\t" <<pulse[11][i]<<"\t" <<std::endl;
   }
}
//-----------------print one track result----------------
if( m_eventNumber < 3 ) {
  for (int it=0; it<3; it++) {
    sprintf(m_cid,"%d : %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f",it,
    pulse[7][it],pulse[8][it],pulse[9][it],pulse[10][it],pulse[11][it],
    pulse[12][it],pulse[13][it],pulse[14][it]);
    cout<<m_cid<<endl;
  }
} 

//----------- find hit strips in the current event------------
 int tbin[21][3];
 for (int i=0;i<21;i++)  { for (int j=0;j<3;j++) tbin[i][j] = 0;} 
 for (int istrip=0 ; istrip<21;istrip++) {
    for(int j=0 ; j<3; j++) {
      if(pulse[istrip][j] < m_threshold) continue;
      tbin[istrip][j]=1;
  }
 }
// ----- find maximum cluster size ---------------
 double max_clsize = 0.;
 bool previous = false;
 double cluster_size = 0.;
 for(int istrip=0 ;istrip<21;istrip++ ) {
     if( tbin[istrip][1] == 1 ) {
        if (!previous) {
          cluster_size = 1.;
          previous = true;
        } else cluster_size += 1.;
     } else { 
        previous = false;
        cluster_size = 0.;
     }
     max_clsize = (cluster_size>max_clsize) ? cluster_size : max_clsize;
  }

  if( m_coutLevel>0 || (m_eventNumber<10 ||fmod(m_eventNumber,2000) == 0)) 
  std::cout<<m_eventNumber<<" p0,e0,x0=" << phi 
  <<" "<<eta<<" "<<x_origin*1.E4 <<"  tbin[8-12]= "
           <<tbin[8][0]<<tbin[8][1]<<tbin[8][2]<<" "
           <<tbin[9][0]<<tbin[9][1]<<tbin[9][2]<<" "
           <<tbin[10][0]<<tbin[10][1]<<tbin[10][2]<<" "
           <<tbin[11][0]<<tbin[11][1]<<tbin[11][2]<<" "
           <<tbin[12][0]<<tbin[12][1]<<tbin[12][2]<<" size= "
           <<max_clsize<<std::endl;
  //if( max_clsize > 0) {
  m_h_phi -> Fill(phi,max_clsize,1);
  //}
//

  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++;
  }}
}