#include <iostream>
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/FundamentalConstants.hh"
#include <cstdlib>
#include <fstream>
#include <TROOT.h>
#include <TApplication.h>
#include <TCanvas.h>
#include <TH1F.h>
#include <TTree.h>
#include <TFile.h>
#include "Garfield/ViewDrift.hh"
#include "Garfield/ComponentAnsys123.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/ViewFEMesh.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/AvalancheMC.hh"
#include "Garfield/Random.hh"
#include "Garfield/TrackHeed.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/ViewMedium.hh"

using namespace Garfield;
using namespace ROOT;
using namespace std;

int main(int argc, char * argv[]) {

  TApplication app("app", &argc, argv);
  //SetDefaultStyle();
// Load the field map.
  ComponentAnsys123 fm;
  fm.Initialise("ELIST.lis", "NLIST.lis", "MPLIST.lis", "PRNSOL.lis", "cm");
  fm.PrintRange();
  fm.PrintMaterials();
// Get the extent of the field map.
  double x0, y0, z0, x1, y1, z1;
  fm.GetBoundingBox(x0, y0, z0, x1, y1, z1);

  ViewField fieldView;
  constexpr bool plotField = true;
  if (plotField) {
    fieldView.SetComponent(&fm);
    // Set the normal vector of the viewing plane (xz plane).
    fieldView.SetPlane(0, 0, -1, 0, 0, 0.5 * (z0 + z1));
    fieldView.PlotContour();
    std::vector<double> xf;
    std::vector<double> yf;
    std::vector<double> zf;
    fieldView.PlotFieldLines(xf, yf, zf, true, true, kOrange-3); 
  }

  Sensor sensor;
  sensor.AddComponent(&fm);//calculate the field
  //sensor.AddElectrode(&fm, "s")//require calculate the electrode signal

  // Setup the gas.
  MediumMagboltz gas;
  gas.LoadGasFile("4He.gas");
  gas.PrintGas();
// Associate the gas with the corresponding field map material. 
  /*const size_t nMaterials = fm.GetNumberOfMaterials();
  for (size_t i = 0; i < nMaterials; ++i) {
    if (fm.GetPermittivity(i) == 1.) fm.SetMedium(i, &gas);
  }*/
  fm.SetMedium(0, &gas);
  fm.PrintMaterials();


  ViewMedium view;
  view.SetMedium(&gas);
  view.PlotElectronVelocity('e');
  //view.PlotElectronDiffusion('e');
  //view.PlotElectronTownsend('e');
  //view.PlotElectronAttachment('e');
  
  // Track class
  TrackHeed track;
  track.SetParticle("proton");
  track.SetEnergy(100.e9);
  track.SetSensor(&sensor);

   //RKF integration
  DriftLineRKF drift;
  drift.SetSensor(&sensor); //chuanganqidaoru
  drift.EnableSignalCalculation();
  
  ViewDrift driftview;
  //the visualize of the electron drift
  TCanvas* CD= nullptr;
  constexpr bool plotDrift = true;
  if (plotDrift)
  {
    CD = new TCanvas("CD", "", 600, 600);
    const bool twod = false;
  	const bool axis = true;
    driftview.SetCanvas(CD);
    driftview.Plot(twod,axis);
    drift.EnablePlotting(&driftview);//plot the track of electron
    track.EnablePlotting(&driftview);//plot the track of proton
  }

  //the calculation of the transportation
  const unsigned int nTracks = 1;//one proton
  double counts,counts_e=0;
  for (unsigned int j=0; j < nTracks; ++j)
  {
    sensor.ClearSignal();
    //const double rTrack = 0.033;
    //const double rTube = 0.053;
    //const double y0 = rTrack;
    //const double z0 = -0.07;
    //const double x0 = -sqrt(rTube * rTube - rTrack * rTrack-0.0049);
    double x0 = -0.1, y0 = 0.204, z0 = -0.07, t0 = 0.;  //make sure the start point of the proton and the direction
    //double t0 = 0.;
    double dx0 = 0., dy0 = 1., dz0 = 0.; 
    track.NewTrack(x0, y0, z0, t0, dx0, dy0, dz0);
    double xc = 0., yc = 0., zc = 0., tc = 0., ec = 0., extra = 0.;
    int nc=0;

    //loop clusters
    while (track.GetCluster(xc, yc, zc, tc, nc, ec, extra))//loop clusters,gain the information of each cluster
    {
      counts++;
      cout<<"cluster position: "<< xc <<"\t"<< yc <<"\t"<< zc<<endl;

      //loop electrons in a cluster
      for (int k = 0; k < nc; k++)//loop electrons in a cluster,one cluster has nc number of electrons
      {
        double xe = 0., ye = 0., ze = 0., te = 0., ee = 0.;
        double dx = 0., dy = 0., dz = 0.;
        track.GetElectron(k, xe, ye, ze, te, ee, dx, dy, dz);//gain the electron information in one cluster
        drift.DriftElectron(xe, ye, ze, te);//calculate the drift of the electron
        
      }
    }

  }

  // Set the signal time window.
  const double tstep = 0.5;
  const double tmin = -0.5 * tstep;
  const unsigned int nbins = 1000;
  sensor.SetTimeWindow(tmin, tstep, nbins);

  
  app.Run(true);
}