//-----------------------------------------------
//   2017.7.18  for ICM (induced charge model) by Taka Kondo (KEK)
//-----------------------------------------------
#include "Ramo.h"
void execute(double& max_clsize) ;
TProfile *h, *m_h_phi ;
TLegend *leg, *leg1;
TLatex *t;

void plot_LA_simulation(){
  TCanvas *c1 = new TCanvas("c1","A Simple Graph ",300,0,750,560); 
  c1->SetGrid(1);
  int HV0[3]={75,100,150};
//------------set main parameters and initialization----------------
  m_isSCTDigiModel = true;  // (old) SCT Digitization model
  m_isSCTDigiModel = false; // induced charge model
  m_isNewCrossTalk = false; 
  m_is25nsGate  = true ;
  m_is25nsGate  = false ;
  m_coutLevel = 3;  // all details
  m_coutLevel = 1;  // all details
  m_coutLevel = 0;  // 0(print every 2000 event)
  m_coutLevel = -1;  // 0(print every 2000 event)
//------------local_phi and eta samplings -------------
  int TotalEventNumber =200000;
  TotalEventNumber =2000;
  double newClusterSize[100]={}, Vbias[100]={};
  int nVbias=0;
//-------after initialization changes -----
  m_phi_min = -20.;
  m_phi_max = 10. ;
  m_eta_min = 1.4;
  m_eta_max = 1.6 ;
  m_eh_pairs = 108.;
  double absEtaL[4]={0.0, 0.4, 0.9, 1.4};
  double absEtaH[4]={0.1, 0.6, 1.1, 1.6};
//m_isLandau=true;
//-----------------
//  leg = new TLegend(0.55,0.73,0.95,0.88);
  double ymax=3.0, ymin=0.9; 
  double xmin = -20., xmax = 10.;
  int nVD = 0;
  int VDlist[8] = {+30,-30,-30};
  //for (int kVD=0; kVD<8; kVD++) {
  for (int kVD = 0; kVD<2; kVD++) {
  //for (int imodel=0; imodel<4; imodel++) {
  for (int imodel=1; imodel<=1; imodel++) {
  c1->Divide(3,1);
    nVD++;
    m_VD = VDlist[kVD];
    switch(imodel){
      case 0: 
        m_isNewCrossTalk = false; m_is25nsGate  = false ;  break;
      case 1: 
        m_isNewCrossTalk = false; m_is25nsGate  = true ;  break;
      case 2: 
        m_isNewCrossTalk = true; m_is25nsGate  = false ;  break;
      case 3: 
        m_isNewCrossTalk = true; m_is25nsGate  = true ;  break;
    }
//----------------------------
    for (int iVB=0; iVB<3; iVB++) {
      c1->cd(iVB+1);
      m_VB = HV0[iVB];
      initialize();
      cout<<"imodel, VB = "<<imodel<<" "<<m_VB<<endl;
   leg = new TLegend(0.54 ,0.76,0.92,0.94," |track eta|");
   for (int ieta=3; ieta>=0; ieta--) {
     m_eta_min = absEtaL[ieta];
     m_eta_max = absEtaH[ieta];
      sprintf(m_cid,"VB= %d, %d",iVB, ieta);
      h = new TProfile(m_cid,"HVscan",60, -20., 10., 0., 50.); 
//------------------------------------
      for (int n=0; n<TotalEventNumber ; n++)  execute(); 
//------get results from h----------------------
    double xx[200]={}, yy[200]={};
    int n=0;
    for (int i=0; i<60; i++){
      //if(h->GetBinContent(i+1) <=0.) {
      xx[i] = h->GetXaxis()->GetBinCenter(i+1);
      yy[i] = h->GetBinContent(i+1);
      cout<<i<<" = "<<xx[i]<<" "<<yy[i]<<endl;
      n++;
    }
//----------plot charges -----------------
    TGraph *h1 = new TGraph(n,xx,yy);
    h1 -> SetMaximum(ymax);
    h1 -> SetMinimum(ymin);
    h1 -> SetMarkerColor(nVD);
    //h -> SetLineColor(1);
    //h1 -> SetLineColor(nVD);
    h1 -> GetXaxis()->SetTitle("inc. angle #phi [degree]");
    h1 -> GetYaxis()->SetTitle("new cluster size");
    h1 -> GetXaxis()->SetLimits(xmin,xmax);
    h1->GetXaxis()->SetLabelSize(0.07);
    h1->GetYaxis()->SetLabelSize(0.07);
    h1->GetXaxis()->SetTitleSize(0.07);
    h1->GetYaxis()->SetTitleSize(0.07);
    h1->GetXaxis()->SetTitleOffset(0.60);
    h1->GetYaxis()->SetTitleOffset(1.25);
    h1->GetXaxis()->SetLabelOffset(-0.02);
    h1->GetYaxis()->SetLabelOffset(0.01);
    h1->SetMarkerColor(ieta+1);
    h1->SetMarkerSize(0.5);
    //sprintf(m_cid,"check only at 30ns");
    //if(imodel==1)sprintf(m_cid,"check for full 25ns");
    //if(imodel==2)sprintf(m_cid,"new crosstalk shape");
    //if(imodel==3)sprintf(m_cid,"25 ns + new crosstalk");
    if( ieta==3) h1 -> Draw("AP");
    else h1 -> Draw("SAME P");
    sprintf(m_cid," %3.1f - %3.1f", m_eta_min, m_eta_max) ;
    leg->AddEntry(h1,m_cid,"p");
}
// leg->SetFillColor(10);leg->SetBorderSize(0);
  //leg->SetTextSize(0.03); leg->Draw(); 
  //gPad->RedrawAxis();
//-------------------------------------
  leg->SetFillColor(10);leg->SetBorderSize(1);
  leg->SetTextSize(0.07); leg->Draw(); 
  TLatex *t;
  t=new TLatex(xmin, (ymax-ymin)*1.01+ymin,"SCT Digitization Model");
  if(!m_isSCTDigiModel)
  t=new TLatex(xmin, (ymax-ymin)*1.01+ymin,"Induced Charge Model");
  t->SetTextSize(0.08); t->Draw(); 
  double xx1 = (xmax-xmin)*0.08+xmin;
  sprintf(m_cid,"V_{FD}=%dV", int(m_VD));
  t=new TLatex(xx1, (ymax-ymin)*0.92+ymin,m_cid);
  t->SetTextSize(0.070); t->Draw(); 
  xx1=(xmax-xmin)*0.50+xmin;
  sprintf(m_cid,"HV=%dV", int(m_VB));
  t=new TLatex(xx1, (ymax-ymin)*0.70+ymin,m_cid);
  t->SetTextSize(0.090); t->Draw(); 
}
//  int Layer=3;
//  sprintf(m_cid,"B%d Run%d HV=%dV",Layer, run,HV);
//  TLatex *t= new TLatex(xmin,(ymax-ymin)*1.01+ymin,m_cid);
//  t->SetTextSize(0.08); t->Draw();

//-------store results-----------
  sprintf(m_cid,"LA_ICM_%dV_%d.png",int(m_VD),imodel);
  if(m_isSCTDigiModel) sprintf(m_cid,"LA_sctDigi_%dV_%d.png",int(m_VD),imodel);
  c1->Print(m_cid);
  c1->Clear();

}  // loop for imodel
}   // loop for VFD
  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 < 5 ) 
    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][200] = {{}};

//---------loop over depth distribution -----------------------
for (int ievent=0; ievent<m_ydivision ; ievent++) {    
   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;
//----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;
   if(m_coutLevel>1) 
     std::cout<<ievent<<" (x_eh,y_eh)="
     <<" ("<<x_eh*1.e4<<","<<y_eh*1.e4<<") strip_offset="<<strip_offset
     <<" (x0,y0)=("<<x0*1e4<<","<<y0*1e4<<")-----"<<std::endl;
//------------- one event loop ---------------------------------
//  temporal storage arrays for every m_transportTimeStep for 5 strips nearby
    double Q_m2[200]={}, Q_m1[200]={}, Q_00[200]={}, Q_p1[200]={}, Q_p2[200]={} ;
//  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<200; 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];
    
//-----------------print one track result----------------
    } // --- end of time bin loop

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

//-----------------print one track result----------------
if( m_eventNumber < 0 || m_coutLevel > 2 ) {
    double totalCharge = 0.;
    for (int it=0; it<200; 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;
}

//-----------put through amplifiers associated with crosstalks-------------
//  storage arrays of pulses at every m_transportTimeStep
int clusterSize = 0;   // 
if( m_is25nsGate ) {
  double Pulse[21][400]={{}}, nsec[400]={};
  for (int it=0; it<400; it++) nsec[it] = (it + 0.5) * m_transportTimeStep;
//  scan over all strips 
  for (int istrip = 0 ; istrip < 21; istrip++ ){
 // scan over all times  
    for (int it = 0 ; it < 200 ; it++ ) {
      double delta_q = qstrip[istrip][it];
      if (delta_q == 0.) continue;  // skip zero cahrge cases
      delta_q *= qfactor;  // change from electron to fC
      double time0 = (it + 0.5) * m_transportTimeStep; // timing of input
      if (time0 > m_GateStartTiming + m_GateWidth) break; //skip later charges
 // add amplified pulses
      for (int jt = 1 ; jt < 400 ; jt++ ) { // trace every m_transportTimeStap
        double timeAmp = jt * m_transportTimeStep  ;
        double time = time0 + timeAmp;
        if (time > m_GateStartTiming + m_GateWidth) break;
        int itime = time / m_transportTimeStep ;
        if( itime >= 400) break;
        double response =  Amp_response(timeAmp)  * delta_q;
        double crosstalk = Amp_crosstalk(timeAmp) * delta_q;
  //cout<<istrip<<" "<<it<<" q="<<delta_q<<" "<<jt<<" -> "<<itime<<endl;
        Pulse[istrip][itime] += response;
        if( istrip > 0 ) Pulse[istrip-1][itime] += crosstalk;
        if( istrip < 19) Pulse[istrip+1][itime] += crosstalk;
      } 
    }  
  }   
  //-------find maximum cluster size within the gate --------
  int hit[21]={} ;
  clusterSize = 0;   // 
  for (int it=0; it<400; it++) {
    double timeOfPulse = (it+0.5) * m_transportTimeStep;
    if ( timeOfPulse < m_GateStartTiming ) continue;
    if ( timeOfPulse > m_GateStartTiming + m_GateWidth ) break;
//-------print center strip pulses -------------
    //if( m_coutLevel > 0 && m_eventNumber < 3) {
    if( m_coutLevel > 0) {
      sprintf(m_cid,"%d(%6.3fns) : %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f",it,
      (it+0.5)*m_transportTimeStep,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;
    }
//----------------------------------------
    for (int istrip=0; istrip<21; istrip++) {
      if( Pulse[istrip][it] < m_threshold ) continue;
      hit[istrip] = 1;
    }
  }
  for (int istrip=0; istrip<21; istrip++) clusterSize += hit[istrip];
  if( m_coutLevel >0 || m_eventNumber%2000==0) {
    for (int i=0; i<21 ; i++ )  cout<<hit[i];
    cout<< "  clusterSize = "<<clusterSize<< " event="<<m_eventNumber<<endl;
  }

//------------  case for fixed timing threshold --------------------
} else {  //  case for fixed timing threshold
  double pulse[21][3]={{}};   // pulse height at m_timeOfThreshold -25,0,+25 ns
  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 < 200 ; it++) {
      double delta_q = qstrip[istrip][it] * qfactor;
      if (delta_q == 0.) continue;
      double time0 = (it + 0.5) * m_transportTimeStep; 
      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;
       }
    } 
  }
//----------- find hit strips in the current event------------
  int tbin[21][3] = {{}};
  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;
  }
  clusterSize = 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;
}

//----------------------------------
  h->Fill(phi, clusterSize);
  return ;
}