#include <iostream>
#include <fstream>
#include <cstdlib>

#include <TCanvas.h>
#include <TROOT.h>
#include <TApplication.h>

#include "Garfield/ViewDrift.hh"

#include "Garfield/ComponentAnalyticField.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/AvalancheMC.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/TrackHeed.hh"
#include "Garfield/TrackSrim.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/ViewDrift.hh"

using namespace Garfield;

int main(int argc, char * argv[]) {
  Time_t start=time(0);
  TApplication app("app", &argc, argv);
 
  // Make a gas medium.
  MediumMagboltz gas;
  gas.LoadGasFile("Ar90_CH4_10_1bar.gas");
  auto installdir = std::getenv("GARFIELD_INSTALL");
  if (!installdir) {
    std::cerr << "GARFIELD_INSTALL variable not set.\n";
    return 1;
  }
  const std::string path = installdir;
  gas.LoadIonMobility(path + "/share/Garfield/Data/IonMobility_Ar+_Ar.txt");

  // Make a component with analytic electric field.
  ComponentAnalyticField cmp;
  cmp.SetMedium(&gas);
  // Wire radius [cm]
  const double rWire = 25.e-4;
  // Outer radius of the tube [cm]
  const double rTube = 0.71;
  // Voltages
  const double vWire = 1500.;
  const double vTube = 0.;
  // Add the wire in the centre.
  cmp.AddWire(0, 0, 2 * rWire, vWire, "s");
  // Add the tube.
  cmp.AddTube(rTube, vTube, 0);

  // Make a sensor.
  Sensor sensor;
  sensor.AddComponent(&cmp);
  sensor.AddElectrode(&cmp, "s");
  const double tstep = 2.;
  const double tmin = -0.5 * tstep;
  const unsigned int nbins = 1000;
  sensor.SetTimeWindow(tmin, tstep, nbins);
  
  TCanvas *cDrift=new TCanvas("","",1200,1200);

  ViewDrift driftview;
  driftview.SetCanvas(cDrift);

  // Set up Heed.
/*
  TrackSrim track;
  track.ReadFile("Alpha_in_Ar_90_CH4_10_1bar.txt");
  track.SetSensor(&sensor);
  track.SetKineticEnergy(2.1e6);  
  track.SetWorkFunction(30.0);
  track.SetFanoFactor(0.3);
  track.SetAtomicMassNumbers(37.6, 17.2);
//  cellview.SetComponent(&fm);
//  cellview.Plot2d();
*/
  TrackHeed track;
  track.SetParticle("electron");
  track.SetKineticEnergy(1.e11);
  track.SetSensor(&sensor);
  track.EnablePlotting(&driftview);

  
  

  // RKF integration.
  DriftLineRKF drift(&sensor);
  // 有雪崩时，平均放大倍数为8倍
  // 放大倍数从形状参数为0的Polya函数中抽样
  // drift.SetGainFluctuationsPolya(0., 1.);
  drift.EnablePlotting(&driftview);
  drift.EnableSignalCalculation();
  drift.EnableIonTail();
 

  const double rTrack = 0.3;
  const double x0 = rTrack;
  const double y0 = -sqrt(rTube * rTube - rTrack * rTrack);
  const unsigned int nTracks = 1;
  for (unsigned int j = 0; j < nTracks; ++j) {
    sensor.ClearSignal();
    track.NewTrack(x0, y0, 0, 0, 0, 1, 0);
    for (const auto& cluster : track.GetClusters()) {
      for (const auto&electrn:cluster.electrons) {
          drift.DriftElectron(electrn.x, electrn.y, electrn.z, electrn.t);
          double ne,ni;
          drift.GetAvalancheSize(ne,ni);
          std::cout<<"ne==   "<<ne<<"ni==   "<<ni<<std::endl;
      }
      for(const auto&ion:cluster.ions)
      {
          drift.DriftIon(ion.x,ion.y,ion.z,ion.t);
      }
    } 
      cmp.PlotCell(cDrift);
      constexpr bool twod = true;
      constexpr bool drawaxis = false;
      driftview.Plot(twod, drawaxis);
      cDrift->Print("drift.png");
    }  
  TCanvas* cS = new TCanvas("cS", "", 1200, 1200);  
  sensor.PlotSignal("s", cS);
  cS->Print("signal.png");

  app.Run(kTRUE);
  time_t end=time(0);
  for(int i=0;i<10;i++)std::cout<<"********************************************"<<std::endl;
  std::cout<<"the application cost:   "<<end-start<<" s"<<std::endl;
  for(int i=0;i<10;i++)std::cout<<"********************************************"<<std::endl;
  return 100000000;  

}