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

#include <TCanvas.h>
#include <TROOT.h>
#include <TApplication.h>
#include <TFile.h>
#include <TTree.h>
#include<TH1D.h>

#include "Garfield/ViewCell.hh"
#include "Garfield/ViewDrift.hh"
#include "Garfield/ViewSignal.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/ViewMedium.hh"

#include "Garfield/ComponentAnalyticField.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/Medium.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/AvalancheMC.hh"
#include "Garfield/TrackSrim.hh"

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

// preamplifier transfer function
double transfer(double t) 
{
	const double tau = 50.;
  const double fC_to_mV = 0.5; // mV per fC
  // Impulse response of the amplifier.
  return -fC_to_mV * exp(4) * pow((t / (4 * tau)), 4) * exp(- 0.2 * t / tau);     
//	return (t / tau) * exp(1 - t / tau);
}

//
bool readTransferFunction(Sensor& sensor) {

  std::ifstream infile;
  infile.open("mdt_elx_delta.txt", std::ios::in);
  if (!infile) {
    std::cerr << "Could not read delta response function.\n";
    return false;
  }
  std::vector<double> times;
  std::vector<double> values;
  while (!infile.eof()) {
    double t = 0., f = 0.;
    infile >> t >> f;
    if (infile.eof() || infile.fail()) break;
    times.push_back(1.e3 * t);
    values.push_back(f);
  }
  infile.close();
  sensor.SetTransferFunction(times, values);  
  return true;
}


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

  TApplication app("app", &argc, argv);
 
/*
void checkhwf() {
  const double work = 30.;
  const double fano = 0.3;
  TH1F* hwf = new TH1F("hwf", "Heed work Fano", 200, 0, 200);
  for (unsigned int i = 0; i < 10000000; ++i) {
    double rnd = 1.e6 * (RndmHeedWF(1.e-6 * work, fano));
    hwf->Fill(rnd);
  }
  TCanvas* chwf = new TCanvas("chwf", "Heed Work Fano", 100, 100, 800, 800);
  chwf->cd();
  hwf->Draw();
  chwf->Update();
  double mean = hwf->GetMean();
  double rms = hwf->GetRMS();
  const double r = rms / mean;
  std::cout << "Histogram mean: " << mean << ", Fano: " << r * r << "\n";
}
*/

   MediumMagboltz gas;
  // Read the electron transport coefficients from a .gas file.
  gas.LoadGasFile("ar_c4h10.gas");
  // Read the ion mobility table from file.
  const std::string path = std::getenv("GARFIELD_HOME");
  gas.LoadIonMobility(path + "/Data/IonMobility_Ar+_Ar.txt");  //Ar离子的迁移率   文件的路径和名字


  //----------------------------
  double amp=0;
  int e_counts=0;
  TFile *file2 = new TFile("Max_info.root", "RECREATE");
  TTree *Max_info = new TTree("Max_info", "Maxium information");
  Max_info->Branch("amp", &amp, "amp/D");
  Max_info->Branch("e_counts", &e_counts, "e_counts/I");
  TH1D *hmax=new TH1D("hmax","The peak information",600,4000,10000); 
  //TH1D *hcns=new TH1D("hmax","The electron counts",10000,0,50000);  

  //---------------------------
  //double amp=0;
  //TFile *file2 = new TFile("Max_info.root", "RECREATE");
  //TTree *Max_info = new TTree("Max_info", "Maxium information");
  //Max_info -> Branch("amp" ,"&amp" ,"amp/D");


  /*
  
  ViewMedium mediumView;
  mediumView.SetMedium(&gas);
  mediumView.PlotElectronVelocity('e');
  */


  // Make a component with analytic electric field.
  ComponentAnalyticField cmp;
  cmp.SetMedium(&gas);   
  // Distance between grid and anode wires [cm].
  constexpr double gap = 0.4;
  // Periodicity (grid spacing) [cm].
  constexpr double periodg = 0.05;  
  // grid diameters.
  constexpr double dg = 0.0050;    
   //anode diameters.
  constexpr double da = 0.0020;     
//-----------------------------------------------------------------------------
  constexpr double vg = 0.;        
  constexpr double va = 1000.;      
//----------------------------------------------------------------------------   

// Periodicity (anode wire spacing) [cm].
  constexpr double perioda = 0.2;   
  cmp.SetPeriodicityX(perioda);    
  // Add the grid wires.
  constexpr double yg = 1 * gap;
  constexpr double xg0 = 0 * periodg;
  constexpr double xg1 = 1 * periodg;
  constexpr double xg2 = 2 * periodg;
  constexpr double xg3 = 3 * periodg;
  cmp.AddWire(xg0, yg, dg, vg, "g");
  cmp.AddWire(xg1, yg, dg, vg, "g");
  cmp.AddWire(xg2, yg, dg, vg, "g");
  cmp.AddWire(xg3, yg, dg, vg, "g");


// Add the anode wires. 
  constexpr double xa = 0.;
  constexpr double ya = 0.;
  cmp.AddWire(xa, ya, da, va, "a");


 // Add the cathode and pad planes.
  constexpr double ycathode = 3.4;
  constexpr double vcathode = -390.;
  constexpr double ypad = -2.6;
  constexpr double vpad = 0;        
  cmp.AddPlaneY(ycathode, vcathode, "cathode");
  cmp.AddPlaneY(ypad, vpad, "pad");


// Request calculation of the weighting field. 
  //cmp.AddReadout("g");
  cmp.AddReadout("a");      
  //cmp.AddReadout("pad");

  const double xmin = -10 * periodg; 
  const double xmax =  10 * periodg;
//----------------------------
  // Make a sensor.
  Sensor sensor;
  sensor.AddComponent(&cmp);
  sensor.SetArea(xmin, ypad, -1., xmax, ycathode, 1.);
  //sensor.AddElectrode(&cmp, "g");
  sensor.AddElectrode(&cmp, "a");
  //sensor.AddElectrode(&cmp,"pad");
   
/*
//-------------------------
  
  ViewField fieldView;
  fieldView.SetComponent(&cmp);
  fieldView.SetArea(xmin, -0.2, xmax, 0.2);
  fieldView.SetVoltageRange(-600., 2000.);
  fieldView.SetNumberOfContours(100);
  fieldView.PlotContour();


//-----------------------
  TCanvas* cf = nullptr;
  constexpr bool plotField = true;
  if (plotField) {
    cf = new TCanvas("cf", "", 600, 600);
    ViewField* fieldView = new ViewField();
    fieldView->SetComponent(&cmp);
    fieldView->SetSensor(&sensor);
    fieldView->SetVoltageRange(-600., 2000.);
    //fieldView->SetElectricFieldRange(0., 20.);
//  fieldView->SetArea(xmin, -1 * gap -h, xmax, 3 * gap + h);
    fieldView->SetArea(xmin, -0.5, xmax, 0.5);
    fieldView->SetCanvas(cf);
    fieldView->Plot("e","CONT4");
}
*/

//------------------------
  //Set the time window and binning for the signal calculation.
  const double tstep = 0.5;
  const double tmin = 0. * tstep;
  const unsigned int nbins = 12000;
  sensor.SetTimeWindow(0, tstep, nbins);
  // Set the delta response function.
  sensor.SetTransferFunction(transfer);
  //if (!readTransferFunction(sensor)) return 0;
  sensor.ClearSignal();


/*
//---------------------------
  // Set up srim.
  TrackSrim track;
  // Connect the track to a sensor.
  track.SetSensor(&sensor);
  //Read Srim output from file
  const std::string file1 = "srim.txt";
  if (!track.ReadFile(file1)) {
    std::cerr << "Reading SRIM file failed.\n";
    return 0;
  }
  // Set the initial kinetic energy of the particle (in eV).
  track.SetKineticEnergy(5.486e5);
  // Set the W value and Fano factor of the gas.
  track.SetWorkFunction(30.0);
  track.SetFanoFactor(0.3);
  // Set A and Z of the gas (not sure what's the correct mixing law).
  const double za = 0.9 * (18. / 40.) + 0.1 * (10. / 16.);
  track.SetAtomicMassNumbers(16. / za, 16);
  // Specify how many electrons we want to be grouped to a cluster.
  track.SetTargetClusterSize(500);
  // tr.SetClustersMaximum(1000);

  // Make some plots of the SRIM data.
  tr.PlotEnergyLoss();
  tr.PlotRange();
  tr.PlotStraggling();
  // Print a table of the SRIM data.
  //tr.Print();

*/



//-----------------------
/*
// RKF integration. 
  DriftLineRKF drift;
  drift.SetSensor(&sensor);
  drift.SetMaximumStepSize(5.e-4);  //Set the maximum step size that is allowed. By default, there is no limit
  //drift.SetGainFluctuationsPolya(0., 10.);   
  drift.EnableIonTail();            //Enable/disable simulation of the ion tail (default: off). 
  drift.EnableSignalCalculation();
*/


//only can calculate electron
  AvalancheMicroscopic aval;
  aval.SetSensor(&sensor);
  //aval.EnableAvalancheSizeLimit(3000.); 
  aval.EnableSignalCalculation();


//calculate electron and ion
  AvalancheMC mc;
  mc.SetSensor(&sensor);
  mc.SetDistanceSteps(2.e-4);
  mc.EnableSignalCalculation();



//----------------------------------
  //const double ymin = -1.5; 
  //const double ymax =  2.5;
  TCanvas* cD = nullptr;
  ViewCell cellView;//Visualize the "cell" defined in an analytic-field component(ComponentAnalyticField)
  ViewDrift driftView;
  cD = new TCanvas("cD", "", 600, 600);
  constexpr bool plotDrift = true;
  if (plotDrift) {
    cellView.SetCanvas(cD);
    cellView.SetComponent(&cmp);
    //cellView.SetArea(xmin, ypad, xmax, ycathode);   
    //cellView.Plot2d();
    driftView.SetArea(xmin, ypad, xmax, ycathode);
    driftView.SetCanvas(cD);        
    //drift.EnablePlotting(&driftView);
    //track.EnablePlotting(&driftView);
    //aval.EnablePlotting(&driftView);
    mc.EnablePlotting(&driftView);
  }




//------------------------- 
  TCanvas* cS = nullptr;
  ViewSignal signalView;//信号可视化
  constexpr bool plotSignal = true;
  if (plotSignal) {
    cS = new TCanvas("cS", "", 600, 600);
    signalView.SetCanvas(cS);
    signalView.SetSensor(&sensor);
    signalView.SetLabelY("signal [fC]");
  } 



//------------------------
  constexpr unsigned int nEvents = 1;
  for (unsigned int i = 0; i < nEvents; ++i) 
  {
    std::cout << i << "/" << nEvents << "\n";
    sensor.ClearSignal();
    const double x0 = 0.025;
    const double y0 = 0.3;
    const double z0 = 0;
    const double t0 = 0;
    const double e0 = 0.1;


     
    aval.AvalancheElectron(x0, y0, z0, t0, e0, 0., 0., 0.);
    //mc.AvalancheElectron(x0, y0, z0, t0);
    double xe1 = 0., ye1 = 0., ze1 = 0., te1 = 0., e1 = 0.;
    int status = 0;
    double xe2 = 0., ye2 = 0., ze2 = 0., te2 = 0., e2 = 0.;          
    const unsigned int np = aval.GetNumberOfElectronEndpoints();
    //const unsigned int np = mc.GetNumberOfElectronEndpoints();
    cout<<"the number of electron trajectories: "<<np<<endl; 
    for(unsigned int j = 0; j < np; ++j)
     {
        aval.GetElectronEndpoint(j, xe1, ye1, ze1, te1, e1, xe2, ye2, ze2, te2, e2, status);
        //mc.GetElectronEndpoint(j, xe1, ye1, ze1, te1, xe2, ye2, ze2, te2, status);
        //aval.DriftElectron(xe1, ye1, ze1, te1, 0., 0., 0.);
        mc.DriftElectron(xe1, ye1, ze1, te1);
        //drift.DriftElectron(xe1, ye1, ze1, te1);
        mc.DriftIon(xe1, ye1, ze1, te1);
        //drift.DriftElectron(xe1, ye1, ze1, te1);
     }


//-----------------------------
    if (plotDrift) {
        cD->Clear();
        cellView.SetArea(xmin, ypad, xmax, ycathode);
        driftView.SetArea(xmin, ypad, xmax, ycathode);
        constexpr bool twod = true;
        constexpr bool drawaxis = false;     
        //driftView.Plot(twod, drawaxis);
        driftView.Plot(twod, twod);    
        //driftView.Plot2D(twod);
        cellView.Plot2d();             
    }


//-------------------------------
   //sensor.AddWhiteNoise("pad", 400000, true, 100);
//   
    //sensor.ConvoluteSignals();   //
    //sensor.ConvoluteSignal("g");
    sensor.ConvoluteSignal("a");
    //sensor.ConvoluteSignal("pad");
    int nt = 1;
    if (!sensor.ComputeThresholdCrossings(-0.1, "a", nt)) continue;
    //if (!sensor.ComputeThresholdCrossings(-0.1, "pad", nt)) continue;
    //if (!sensor.ComputeThresholdCrossings(-0.1, "g", nt)) continue;

//------------------------------
   double f_signal=0;
   double t_signal=0;
   //double tsignal_min = std::numeric_limits<double>::max();
   for (unsigned int i = 0; i < nbins ; ++i) {
   t_signal=i*tstep;
   const double fsignal = sensor.GetSignal("a", i);
   //const double fsignal = sensor.GetSignal("g", i);
   //const double fsignal = sensor.GetSignal("pad", i);
   f_signal = std::max(f_signal, fsignal);
    }
   amp = f_signal;
   Max_info -> Fill();
   hmax -> Fill(amp);
   //hcns ->Fill(e_counts);
   cout << "Max_value: " << amp << endl;
   


//-----------------------------
    if (plotSignal) 
    {
      //signalView.PlotSignal("g");
      signalView.PlotSignal("a");    
      //signalView.PlotSignal("pad");     
      //TFile *hist = new TFile("Simu_info.root", "RECREATE");
      //cS->Write("signal");            
      //cD->Write("drift");             
      //hist->Close();
    }


  }

   Max_info->Write();
   hmax->Write();
   //hcns->Write();
   //Max_info->StartViewer();
   file2->Close();

  cout << " process completed! " << endl;
  app.Run(true);

}