//   eTransport4.C
//   02 Aug 2019
//   Use "matrix" output from ESTAT rather than .EOU output
//   matrix files are regular x-y grids
//   Geometry of XP2020 phototube near photocathode
//   Track photoelectrons from photocathode to wherever
//   Read in E-fields for each electrode at 1 Volt
//   Add them together, weighted with the proper boundary
//   potentials and compare results to those from eTransport.C

//#include "Riostream.h"
#include <Riostream.h>
#include <string>
#include <fstream> 
#include <iomanip> 
#include <cmath> 
#include <TCanvas.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TLatex.h>
#include <TFrame.h>
#include <TLine.h>
#include <TTimeStamp.h>

struct focusData {
  int KMax,LMax,NReg,ICylin;
  double XMin,XMax,YMin,YMax;
  double DUnit;
  double dxg,dyg;
  int RgNo[651][251];
  double x[651][251],y[651][251],phi[651][251],Ex[651][251],Ey[651][251];
  double divE[650][250];
} F,D1,G1,G2,K1,K2,K3,K4,K5,K6;

struct Rrings {   // Resistance between photocathode rings (MOhm)
  double Resist[6]={1591.6, 
                     795.8,
                     530.5,
                     397.9,
                     256.2,
                      98.5};
} R;

double KResist(int k) {
  // Return resistance between photocathode rings
  // Scale these resistors to correcpond to different photocathode resistivities
  // other than the 1 GOhm/square nominal value.
  return (R.Resist[k]);
}

void ReadFile(char* filename, focusData& C) {

  // Read in files created in Field Precision EStat 8.0.  In EStat the potentials
  // are independent of the scale of the physical dimensions.  In EStat I specified
  // the physical dimensions to be meters when, in fact, the physical dimensions are
  // mm.  So, when EStat calculates E-fields it assumes the physical dimensions
  // are in meters.  It then takes the difference of two potentials separated by a
  // step of 0.1 and divides that difference by 0.1, thinking the units are meters
  // when, in fact, the step is in mm.  In eTransport we want the units of E-field
  // to be in V/m so we must multiply the E-fields given by EStat by 1000.

  int kk,ll,ru,rd;
  char asc[200];
  ifstream inp;
  inp.open(filename);
  if(!inp.is_open()) printf("Could not open %s\n",filename);
  inp.getline(asc,200); // skip title lines
  inp >> asc >> C.XMin >> asc >> C.YMin; inp.getline(asc,100); // skip title lines
  inp >> asc >> C.XMax >> asc >> C.YMax; inp.getline(asc,100); // skip title lines
  inp >> asc >> C.KMax >> asc >> C.LMax; inp.getline(asc,100); // skip title lines
  for(int i=0;i<2;i++) {inp.getline(asc,200);} // Skip headers
  for(int l=0;l<=C.LMax;l++) {
    for(int k=0;k<=C.KMax;k++) {
      inp >> C.x[k][l] >> C.y[k][l] >> C.RgNo[k][l] 
          >> C.phi[k][l] >> C.Ex[k][l] >> C.Ey[k][l];
      inp.getline(asc,100);
      C.Ex[k][l]*=1000.; C.Ey[k][l]*=1000.; // Correct units so distances are in mm
    }
  }
  inp.close();
  C.dxg=(C.XMax-C.XMin)/double(C.KMax); C.dyg=(C.YMax-C.YMin)/double(C.LMax);
  //printf("KMax, LMax, dxg,dyg: %3d %3d %8.3f %8.3f\n",C.KMax,C.LMax,C.dxg,C.dyg);
}

void ReadVfields(int HVbase) {
  // HVbase is the HV potential (V) supplied to the phototube base
  char cHVbase[]="    ", ecHVbase[6];
  sprintf(ecHVbase,"%d",HVbase+10000);
  for(int i=0;i<4;i++) {cHVbase[i]=ecHVbase[i+1];}
  cHVbase[4]='\0';
  char filename[100];
  sprintf(filename,"DYNODE-1.MTX");
  ReadFile(filename,D1);
  sprintf(filename,"GRID-1.MTX");
  ReadFile(filename,G1);
  sprintf(filename,"GRID-2.MTX");
  ReadFile(filename,G2);
  sprintf(filename,"KATHODE-1.MTX");
  ReadFile(filename,K1);
  sprintf(filename,"KATHODE-2.MTX");
  ReadFile(filename,K2);
  sprintf(filename,"KATHODE-3.MTX");
  ReadFile(filename,K3);
  sprintf(filename,"KATHODE-4.MTX");
  ReadFile(filename,K4);
  sprintf(filename,"KATHODE-5.MTX");
  ReadFile(filename,K5);
  sprintf(filename,"KATHODE-6.MTX");
  ReadFile(filename,K6);
}

void CombineVfields(int HVint, double Vring[]) {
  double Vdyn1=40/10/21.5*double(HVint); // Potential at first dynode relative to photocathode
  double Vgrd1=14/10/21.5*double(HVint); // Potential at first grid relative to photocathode
  double Vgrd2=90/10/21.5*double(HVint); // Potential at second grid relative to photocathode
  // Vdyn1=158.; Vgrd1=55.4; Vgrd2=355.; // **************** Temp test
  F.XMin=K1.XMin;
  F.XMax=K1.XMax;
  F.KMax=K1.KMax;
  F.YMin=K1.YMin;
  F.YMax=K1.YMax;
  F.LMax=K1.LMax;
  F.DUnit=K1.DUnit;
  F.NReg=K1.NReg;
  F.ICylin=K1.ICylin;
  for(int l=0;l<=F.LMax;l++) {
    for(int k=0;k<=F.KMax;k++) {
      F.RgNo[k][l]=K1.RgNo[k][l];
      F.x[k][l]=K1.x[k][l]; F.y[k][l]=K1.y[k][l];
      F.phi[k][l]=Vdyn1*D1.phi[k][l]+Vgrd1*G1.phi[k][l]+Vgrd2*G2.phi[k][l]+
         Vring[0]*K1.phi[k][l]+Vring[1]*K2.phi[k][l]+Vring[2]*K3.phi[k][l]+
         Vring[3]*K4.phi[k][l]+Vring[4]*K5.phi[k][l]+Vring[5]*K6.phi[k][l];
      F.Ex[k][l]=Vdyn1*D1.Ex[k][l]+Vgrd1*G1.Ex[k][l]+Vgrd2*G2.Ex[k][l]+
         Vring[0]*K1.Ex[k][l]+Vring[1]*K2.Ex[k][l]+Vring[2]*K3.Ex[k][l]+
         Vring[3]*K4.Ex[k][l]+Vring[4]*K5.Ex[k][l]+Vring[5]*K6.Ex[k][l];
      F.Ey[k][l]=Vdyn1*D1.Ey[k][l]+Vgrd1*G1.Ey[k][l]+Vgrd2*G2.Ey[k][l]+
         Vring[0]*K1.Ey[k][l]+Vring[1]*K2.Ey[k][l]+Vring[2]*K3.Ey[k][l]+
         Vring[3]*K4.Ey[k][l]+Vring[4]*K5.Ey[k][l]+Vring[5]*K6.Ey[k][l];
    }
  }
  F.dxg=K1.dxg; F.dyg=K1.dyg;
}

void findE(double x, double y, double& Ex, double& Ey) {
  double ay=fabs(y);
  double rx,ry,a,b,c,z1,z2;
  double xk,yl,dx,dy;
  int k,l;
  xk=(x-F.XMin)/F.dxg; yl=(ay-F.YMin)/F.dyg; k=xk; l=yl;
  if(k<=0) k=0; if(k>=F.KMax) k=F.KMax-1;
  if(l<=0) l=0; if(l>=F.LMax) l=F.LMax-1;
  /*
  dx=x-F.x[k][l]; dy=ay-F.y[k][l];
  rx=dx/F.dxg;    ry=dy/F.dyg;
  //Ex=F.Ex[k][l]+(F.Ex[k+1][l]-F.Ex[k][l])*dx/F.dxg;
  //Ey=F.Ey[k][l]+(F.Ey[k][l+1]-F.Ey[k][l])*dy/F.dyg;
  //Ex=F.Ex[k][l]; Ey=F.Ey[k][l]; // Temporary ****************
  z1=F.Ex[k+1][l]; z2=F.Ex[k][l+1];
  if((rx+ry)>1.0) {c=F.Ex[k+1][l]+F.Ex[k][l+1]-F.Ex[k+1][l+1];}
  else            {c=F.Ex[k][l];}
  Ex=c-(c-z1)*rx-(c-z2)*ry;
  z1=F.Ey[k+1][l]; z2=F.Ey[k][l+1];
  if((rx+ry)>1.0) {c=F.Ey[k+1][l]+F.Ey[k][l+1]-F.Ey[k+1][l+1];}
  else            {c=F.Ey[k][l];}
  Ey=c-(c-z1)*rx-(c-z2)*ry;
  */
  Ex=-(F.phi[k+1][l]-F.phi[k][l])/F.dxg*1000.; // Make E-field
  Ey=-(F.phi[k][l+1]-F.phi[k][l])/F.dyg*1000.; // units in V/m
  if(y<0.0) Ey=-Ey;
}

int findRgNo(double x, double y) {
  double ay=fabs(y);
  int k,l;
  k=(x-F.XMin)/F.dxg; l=(ay-F.YMin)/F.dyg;
  if(k<=0) k=0; if(k>=F.KMax) k=F.KMax-1;
  if(l<=0) l=0; if(l>=F.LMax) l=F.LMax-1;
  return (F.RgNo[k][l]);
}

void divergE(focusData& C) {
  // Calculate divergence of E-field in file C
  for(int l=0;l<250;l++) {
    for(int k=0;k<650;k++) {
      C.divE[k][l]=
        (C.y[k][l+1]*C.Ey[k][l+1]-C.y[k][l]*C.Ey[k][l])/C.dyg/(C.y[k][l+1]+C.y[k][l])*2.
        +(C.Ex[k+1][l]-C.Ex[k][l])/C.dxg;
    }
  }
  FILE * divFile;
  divFile = fopen("divergE.txt","w");
  for(int l=0;l<250;l++) {
    for(int k=0;k<650;k++) {
      fprintf(divFile,"%3d %3d %14.6f\n",k,l,C.divE[k][l]);
    }
  }
  fclose(divFile);
}

// ****************************************************************************** //
// ****************************************************************************** //

void eTransport4() {

  int HVint;
  double QKphoto; // Charge (in fC) in photoelectron pulses (at rate of 1 MHz)
  // Using 0.75 photoelectrons/proton (See "Calibration.docx" of 5 Mar 2019)
  // we get the following useful numbers:
  // 1 fC = 6241 photoelectrons = 8321 protons
  // Max charge I've used here is 30 fC which corresponds to 2.5e5 protons
  // Run 1137 has bunches in this range with some bunches up to 7 times higher
  // Of course the phototube is well saturated so the actual number of protons
  // is much higher.
  double Vring[]={0.0,0.0,0.0,0.0,0.0,0.0};
  double Ikat[6];
  printf("\nEnter Phototube HV (in Volts) as an integer: ");
  scanf("%d",&HVint);
  printf("Enter charge (in fC) in photoelectron 1 MHz pulses: ");
  scanf("%lf",&QKphoto);
  printf("\n");
  ReadVfields(HVint); // Read in potential fields

  char PlotTitle1[]="photoelectron transport";
  TCanvas* c1=new TCanvas("c1",PlotTitle1,900,40,1000,770);
  c1->SetGrid();          //put grid on graph
  c1->SetFillColor(10);   //white
  c1->SetGrid();

  double zlo,zhi,rlo,rhi;
  zlo=0.0; zhi=65.0; rlo=-25.00; rhi=25.00;
  TH2F *hpx1=new TH2F("hpx1","",10,zlo,zhi,10,rlo,rhi);
  hpx1->SetStats(kFALSE);
  hpx1->GetXaxis()->SetTitle("  z (mm)");
  hpx1->GetXaxis()->SetTitleSize(0.07);
  hpx1->GetXaxis()->SetTitleOffset(0.6);
  hpx1->GetXaxis()->CenterTitle();
  hpx1->GetXaxis()->SetLabelSize(0.045);
  hpx1->GetYaxis()->SetLabelSize(0.045);
  hpx1->GetYaxis()->SetNdivisions(1005);
  hpx1->GetYaxis()->SetTitle("   r (mm)");
  hpx1->GetYaxis()->SetTitleSize(0.07);
  hpx1->GetYaxis()->SetTitleOffset(0.6);
  hpx1->GetYaxis()->CenterTitle();
  hpx1->SetFillColor(10);  //white

  TGraph *gr1=new TGraph();
  gr1->SetLineColor(4);
  gr1->SetLineWidth(1);  // Linewidth of error bars
  gr1->SetLineStyle(1);  // Linestyle solid

  int niter=1; // Requested number of iterations
  int iter=0;  // Total number of iterations

  while(1) { // Loop over iterations

    if(niter<=0) {
      printf("How many more iterations? ");
      scanf("%d",&niter);
      printf("\n");
    }
    if(niter<=0) break;
    iter++; // Keep track of how many iterations

    CombineVfields(HVint,Vring);
    //    printf("Ex[15][22],Ey[15][22]: %14.6e %14.6e\n",F.Ex[15][22],F.Ey[15][22]);
    divergE(F); // ************** Test E-field divergence

    int RgNoDyn1=2; // Region number for first Dynode
    int RgNoPk1=4, RgNoPk2=9;   // Region number limits for PhotoKathode

    if(niter==1) {c1->Clear(); hpx1->DrawCopy();}

    double vz,vr,az,ar,dvz,dvr,dz,dr;
    double z[2000],r[2000];
    double zint,rint;
    int npts;
    double Ez,Er;
    double Eline; // Potential drop along photoelectron track
    double me=9.11E-31;  // Electron mass in kg
    double mc2=0.511E6;  // Electron rest mass energy in eV
    double c=3E11;       // Speed of light in mm/s
    double eq=1.602E-19; // Electron charge in Coulombs
    double t,dt=0.05E-9;  // Time step in seconds
    double T=0.0, T2=0.0, sigmaT=0.0;
    int RgNoC;
    bool vac;
    int nGen=5000;
    //    int nGen=5; // This is for a test
    int nGnt[6];  // How many generated at kathode ring i
    int nDyn1[6]; // How many from kathode ring i reached first dynode
    int nPk[6];   // How many from kathode ring i returned to any kathode ring
    int kPk[6][6];// How many from kathode ring i returned to kathode ring j
    int mPk[6];   // How many from any kathode ring returned to kathode ring j
    int nDyn1T;   // Total number of photoelectrons reaching first dynode
    double RingWidth=4.4; // Width of each kathode ring
    double Dr=RingWidth*5/double(nGen); // Generate photoelectrons only out of first 5 rings
    //    double z0=1.5; // Starting position
    double z0=1.1; // Starting position
    double KE=1.0; // Starting kinetic energy in eV

    for(int j=0;j<6;j++)
      {nGnt[j]=0; nDyn1[j]=0; nPk[j]=0; mPk[j]=0; for(int k=0;k<6;k++) kPk[j][k]=0;}

    for(int i=0;i<nGen;i++) { // Loop over generated photoelectrons
      int j;
      t=0.0;
      z[0]=z0; r[0]=Dr*(0.5+double(i));
      j=r[0]/RingWidth;
      nGnt[j]++;
      Eline=0.0;
      vz =c*sqrt(2.*KE/mc2); // starting speed in mm/s
      vr =0.0;
      npts=0;
      vac=true;
      while (vac&&(npts<1999)) {
        zint=z[npts]+vz*dt/2.; rint=r[npts]+vr*dt/2.;
        findE(zint,rint,Ez,Er);
        npts++;
        az=-eq*Ez*1E3/me;     ar=-eq*Er*1E3/me; // Convert acceleration to mm/s^2
        dvz=az*dt;            dvr=ar*dt;
        vz+=dvz;              vr+=dvr;
        dz=vz*dt;             dr=vr*dt;
        Eline+=(Ez*dz+Er*dr)/1000.; // Accumulate line integral of E-field along electron path
        z[npts]=z[npts-1]+dz; r[npts]=r[npts-1]+dr;
        t+=dt;
        RgNoC=findRgNo(z[npts],r[npts]);
        if((z[npts]<F.XMin)||(z[npts]>F.XMax)||(r[npts]<-F.YMax)||(r[npts]>F.YMax)) RgNoC=12;
	// printf("npts,t,z,r,RgNo,Ez,Er: %4d %14.6e %9.4f %9.4f %3d %14.6e %14.6e\n",npts,t,z[npts],r[npts],RgNoC,Ez,Er);
        vac=RgNoC==1; // True if electron is still in the vacuum
      } // end while (vac) loop over photoelectron time steps
      //        printf("Eline,RgNo,npts: %10.2f V %5d %5d\n",Eline,RgNoC,npts); // Print out E-field line integral (potential difference)
      if(j<5&&niter==1) gr1->DrawGraph(npts,z,r,"L");
      //    printf("i,z,r: %3d %9.2f %9.2f\n",i,z[npts-1],r[npts-1]);
      //    printf("i,npts,RgNoC,t: %4d %5d %4d %8.1f\n",i,npts,RgNoC,1e9*t);
      if(RgNoC==RgNoDyn1) {  // Region number for first dynode
        nDyn1[j]++;          // Number of photoelectrons reaching the first dynode
        T+=1E9*t;          // Transit time in ns to reach the first dynode
        T2+=1E18*t*t;      // Square of transit time
      }
      if((RgNoC>=RgNoPk1)&&(RgNoC<=RgNoPk2)) {
        nPk[j]++;
        kPk[j][RgNoC-RgNoPk1]++;
        mPk[RgNoC-RgNoPk1]++;
      }
    } // end loop over generated events
    nDyn1T=0;
    for(int j=0;j<5;j++) {nDyn1T+=nDyn1[j];}
    if(T>1.01*1E9*dt&&nDyn1T>0) { 
      T/=double(nDyn1T);    // average transit time in ns
      T2/=double(nDyn1T);   // 
      sigmaT=sqrt(T2-T*T); // rms spread of transit times
    }
    double IKath[6],sIKath[6],sIkat=0.0,Idyn1=0.0;
    double rin,rout;
    double area[5]; // Area of ring as a fraction of the whole
    for(int i=0;i<5;i++) {
      rin=0.2*double(i); rout=0.2*double(i+1);
      area[i]=rout*rout-rin*rin; // Sum over ring areas is unity
      Ikat[i]=area[i]*QKphoto/1000.; // Time averaged current in uA
    }
    Ikat[5]=0.0;
    for(int j=0;j<6;j++) {
      IKath[j]=Ikat[j];
      for(int i=0;i<5;i++) {if(nGnt[i]!=0) {
	  IKath[j]-=Ikat[i]*double(kPk[i][j])/double(nGnt[i]);}}
      if(j==0) sIKath[j]=IKath[j];
      else     sIKath[j]=sIKath[j-1]+IKath[j];
      sIkat+=Ikat[j];
      if(nGnt[j]!=0) Idyn1+=Ikat[j]*double(nDyn1[j])/double(nGnt[j]);
    }
    double Vnew[6];
    for(int i=5;i>=0;i--) {
      if(i==5) Vnew[i]=KResist(i)*sIKath[i];
      else     Vnew[i]=Vnew[i+1]+KResist(i)*sIKath[i];
    }
    // Now estimate new ring potentials
    for(int i=0;i<6;i++) Vring[i]+=(Vnew[i]-Vring[i])/10.; // Soften approach ???

    printf("HV: %4d V   photoK: %6.2f fC   sIkat: %8.3f nA   Idyn1: %8.3f nA\n",
	   HVint,QKphoto,sIkat*1000.,Idyn1*1000.);
    printf("aveT,sigT: %8.2f %8.2f       Iteration number %4d\n",T,sigmaT,iter);
    printf("\n  jth   Potntl   photoI   nGnt  nDyn1");
    printf("     No. from the ith ring to the jth ring");
    printf("  Ikathod   NewPot\n");
    printf(" ring      (V)     (nA)              ");
    printf("   i= 1      2      3      4      5      6");
    printf("     (nA)      (V)\n");
    for(int j=0;j<5;j++) {
      printf("%5d %8.3f %8.3f %6d %6d ",
         j+1,Vring[j],Ikat[j]*1000.,nGnt[j],nDyn1[j]);
      for(int i=0;i<6;i++) printf("%6d ",kPk[i][j]);
      printf("%8.3f %8.3f\n",IKath[j]*1000.,Vnew[j]);
    }
    int jj=5;
    char space[]="             ";
    printf("%5d %8.3f %8.3f %s ",
      jj+1,Vring[jj],Ikat[jj],space);
    for(int i=0;i<6;i++) printf("%6d ",kPk[i][jj]);
    printf("%8.3f %8.3f\n\n",IKath[jj]*1000.,Vnew[jj]);
    niter--;
  } // End while(1) loop over iterations

  // Complete the plot

  char sQKphoto[40];
  sprintf(sQKphoto,"HV %4d V, PhotoK Charge Pulse %6.2f fC",HVint,QKphoto);

  TLatex *sctex = new TLatex();
  sctex->SetTextAlign(13);
  sctex->SetTextFont(42);
  sctex->SetTextSize(0.03);
  sctex->SetLineWidth(1);
  sctex->DrawLatex(30.,24.,sQKphoto);

  TLine *OutLin = new TLine();
  OutLin->SetLineWidth(12);
  OutLin->SetLineColor(29);
  OutLin->DrawLine(60.50, -8.00,60.50,  8.00);  // First dynode  RgNo  2
  OutLin->SetLineColor(1);
  OutLin->DrawLine(50.50,  4.00,50.50, 21.00);  // Grid 2        RgNo 11
  OutLin->DrawLine(50.50, -4.00,50.50,-21.00);
  OutLin->DrawLine(33.00, 17.50,50.10, 17.50);  
  OutLin->DrawLine(33.00,-17.50,50.10,-17.50);
  OutLin->DrawLine( 2.00, 22.50,24.00, 22.50);  // Shield        RgNo  3
  OutLin->DrawLine( 2.00,-22.50,24.00,-22.50);
  OutLin->DrawLine(20.00, 13.50,37.00, 11.50);  // Grid 1        RgNo 10
  OutLin->DrawLine(20.00,-13.50,37.00,-11.50);
  OutLin->DrawLine(28.50, 12.10,28.50, 10.00);
  OutLin->DrawLine(28.50,-12.10,28.50,-10.00);
  OutLin->SetLineColor(2);
  OutLin->DrawLine( 0.50, -4.32, 0.50,  4.32);  // Kathode 1     RgNo  4  Inner most ring
  OutLin->SetLineColor(3);
  OutLin->DrawLine( 0.50,  4.48, 0.50,  8.72);  // Kathode 2     RgNo  5
  OutLin->DrawLine( 0.50, -4.48, 0.50, -8.72);
  OutLin->SetLineColor(4);
  OutLin->DrawLine( 0.50,  8.88, 0.50, 13.12);  // Kathode 3     RgNo  6
  OutLin->DrawLine( 0.50, -8.88, 0.50,-13.12);
  OutLin->SetLineColor(5);
  OutLin->DrawLine( 0.50, 13.28, 0.50, 17.52);  // Kathode 4     RgNo  7
  OutLin->DrawLine( 0.50,-13.28, 0.50,-17.52);
  OutLin->SetLineColor(6);
  OutLin->DrawLine( 0.50, 17.68, 0.50, 21.92);  // Kathode 5     RgNo  8
  OutLin->DrawLine( 0.50,-17.68, 0.50,-21.92);
  OutLin->SetLineColor(7);
  OutLin->DrawLine( 0.50, 22.08, 0.50, 25.02);  // Kathode 6     RgNo  9  Outer most ring
  OutLin->DrawLine( 0.50,-22.08, 0.50,-25.02);

}  // End of eTransport2