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

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

#include "Garfield/ViewDrift.hh"
#include "Garfield/ViewSignal.hh"
#include "Garfield/ViewFEMesh.hh"
#include "Garfield/ViewField.hh"

#include "Garfield/DriftLineRKF.hh"
#include "Garfield/ComponentComsol.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/AvalancheMC.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/SolidBox.hh"
#include "Garfield/GeometrySimple.hh"
#include "Garfield/ComponentConstant.hh"
#include "Garfield/TrackSrim.hh"
#include "Garfield/Random.hh"
#include "Garfield/Plotting.hh"
#include "Garfield/TrackHeed.hh"

using namespace Garfield;

//电流前端电子学卷积函数
double Transferfunction(double t){     
  /* constexpr double risistance= 100 ;   //电阻设置
  constexpr double capacitance=  100;    //电容设置 */
  //constexpr double tau=25;//ns
  //return exp(-t/tau);
  constexpr double tau = 160;
  constexpr double fC_to_mV = 12.7;
  return fC_to_mV * exp(4) * pow((t / tau), 4) * exp(-4 * t / tau); 
}

//原始电流信号打印到txt文件函数
void printsignaltotxt(Sensor& sensor,const std::string label,const int nTimebins,const double tstep,const double tmin){     
  char title[60];
  sprintf(title,"currentsignal Xe_%.1f_CO2_%.1f_300K_%.2fAtm_%.0fV.txt",85.0,15.0,1.0,1.0); 
  std::ofstream outfile;
  outfile.open(title,std::ios::out);
  outfile << "time（ns）     "<< "total current     "<< "electron current     " << "ion current\n";
  for(unsigned int i=0;i<nTimebins;++i){
    const double t=tstep/2+tmin+i*tstep;
    const double f=sensor.GetSignal(label,i);
    const double fe=sensor.GetElectronSignal(label,i);
    const double fion=sensor.GetIonSignal(label,i);
    outfile << t << "      "  << f << "      "  <<fe << "      "   << fion << "\n";
  }
}

//原始voltage信号打印到txt文件函数
void printvol(Sensor& sensor,const std::string label,const int nTimebins,const double tstep,const double tmin){     
  char title[60];
  sprintf(title, "signal_vol.txt"); 
  std::ofstream outfile;
  outfile.open(title,std::ios::out);
  outfile << "time(ns)     "<< "total vol     "<< "electron vol   " << "ion vol\n";
  for(unsigned int i=0;i<nTimebins;++i){
    const double t=tstep/2+tmin+i*tstep;
    const double f=sensor.GetSignal(label,i);
    const double fe=sensor.GetElectronSignal(label,i);
    const double fion=sensor.GetIonSignal(label,i);
    outfile << t << "      "  << f << "      "  <<fe << "      "   << fion << "\n";
  }
}



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

  TApplication app("app", &argc, argv);

  //导入电场文件COSMOL、材料关联
  ComponentComsol fm;
  fm.Initialise("mesh.mphtxt", "dielectrics.dat", "field.txt", "mm");
  fm.SetWeightingField("weight.txt", "wf");
  fm.PrintRange();

  //电离室长宽高(cm)10mm
  const double xhalf = 1;  

  // Setup the gas.
  MediumMagboltz gas;
  gas.SetComposition("xe", 100.); //100%Xe
  gas.SetTemperature(293.15);
  gas.SetPressure(3*760.);//yaqiang
  gas.EnableAutoEnergyLimit(false);
  gas.SetMaxElectronEnergy(1200);
  gas.Initialise(true);  
  // Load the ion mobilities.
  const std::string path = std::getenv("GARFIELD_HOME");
  //gas.LoadIonMobility(path + "/Data/IonMobility_Ar+_Ar.txt");
  gas.LoadIonMobility(path + "/Data/IonMobility_Xe+_P12_Xe.txt");//Xe离子种类需要改进，看论文
  // Associate the gas with the corresponding field map material. 
  const unsigned int nMaterials = fm.GetNumberOfMaterials();
  for (unsigned int i = 0; i < nMaterials; ++i) {
    const double eps = fm.GetPermittivity(i);
    if (eps == 1.) fm.SetMedium(i, &gas);//金属的介电常数需要查询，不能为1 ，可以设置为1e10
  }
  fm.PrintMaterials();

  // Create the sensor.
  Sensor sensor;
  sensor.AddComponent(&fm);
  sensor.AddElectrode(&fm,"wf");     
  sensor.SetArea(-4.5 * xhalf, -0.275 * xhalf, -0.8* xhalf, 
                  4.5 * xhalf, 0.275 * xhalf, 0.8 * xhalf);

  //设置信号计算的时间下时间限，步长及时间步长个数，单位ns
  const double tstep = 100;      //**********  //信号采样步长，0.1微秒 （ns,s=1e9ns）
  const double tmin = 0;     //********** //信号采样时间下限
  const double tmax=1000000;         //*****上限，根据结果调整  1ms   
  const unsigned int nbins = (tmax-tmin)/tstep;   //信号采样时间间隔数 ,随机取样
  sensor.SetTimeWindow(tmin, tstep, nbins);
  sensor.ClearSignal();

  AvalancheMC drift;           //设置一个AvalancheMC用以计算初电离离子漂移
  drift.SetSensor(&sensor);     //链接Sensor
  drift.SetDistanceSteps(2.e-4);   //设置蒙卡计算离子漂移时的步长    
  drift.EnableSignalCalculation();    //开启AvalancheMC漂移模拟的电流信号计算

  // Set up the charged particle track.
  TrackHeed track;
  track.SetSensor(&sensor);
  track.SetParticle("muon");
  track.SetEnergy(170.e9);

  ViewDrift driftView;
  constexpr bool plotDrift = true;
  if (plotDrift) {
    track.EnablePlotting(&driftView);//电子漂移线绘制数据传递 (DriftLineRkF中控制项，会将指向一个ViewDrift类的指针‘&driftView’值传递给一个名为m_iew的变量供后续操作)
    drift.EnablePlotting(&driftView);//离子漂移线
  }

/*   //设置离子与电子    1mm=0.1cm
  constexpr unsigned int nEvents = 100;//离子数与电子数
  for (unsigned int i = 0; i <= nEvents; ++i) { 
    std::cout << i << "/" << nEvents << "\n";
    // Randomize the initial position. 距离单位cm
    const double x0 = -4.5* RndmUniform();//45mm
    const double y0 = 0.198* RndmUniform()+0.026;//0.2249
    const double z0 = 0.9* RndmUniform();  //9.19mm
    const double t0 = 0.;
    driftion.DriftIon(x0, y0, z0, t0);
  } */
  //设置离子与电子    1mm=0.1cm
  const double x0 = 0;//45mm
  const double y0 = 0.225;//0.275;//0.2249;//0.199* RndmUniform()+0.026;
  const double z0 = 0.75;//0.75;//7.5mm
  constexpr unsigned int nEvents = 1;//lizishu
  for (unsigned int i = 0; i < nEvents; ++i) { 
    sensor.ClearSignal();
    //track.NewTrack(x0, y0, z0, 0, 0, 0,-1);
    track.NewTrack(x0, y0, z0, 0, 0, -0.570,-0.821);//方向矢量
    // Randomize the initial position. 距离单位cm
    double xc = 0., yc = 0., zc = 0., tc = 0., ec = 0., extra = 0.;
    int nc = 0;
    while(track.GetCluster(xc, yc, zc, tc, nc, ec, extra) ){
      for (int j = 0; j < nc; ++j){
        double xe0, ye0, ze0, te0, ee0, dxe0, dye0, dze0;
        //double x1, y1, z1, t1, e1, dx1, dy1, dz1;
        track.GetElectron(j, xe0, ye0, ze0, te0, ee0, dxe0, dye0, dze0);
        drift.DriftElectron(xe0, ye0, ze0, 0.);
        //track.GetIon(j, x1, y1, z1);
        //drift.DriftIon(x1, y1, z1, 0.);
      }
    }
  }
   

  //漂移线绘制
  if (plotDrift) {
    TCanvas* cd = new TCanvas();
    ViewFEMesh* meshView = new ViewFEMesh();
    meshView->SetArea(-4.5 * xhalf, -0.42 * xhalf, -0.92* xhalf, 
                        4.5 * xhalf, 0.42 * xhalf, 0.92 * xhalf);
    meshView->SetCanvas(cd);
    meshView->SetComponent(&fm);
    // y-z projection.
    meshView->SetPlane(1, 0, 0, 0, 0, 0);
    meshView->SetFillMesh(true);
    // Set the color of the 304.map ID,不锈钢上色
    meshView->SetColor(1, kRed);
    meshView->EnableAxes();
    meshView->SetViewDrift(&driftView);
    meshView->Plot();
    cd->SaveAs("dirftline.jpg");//输出漂移线图片
  }   

  // Plot the induced current.
  ViewSignal signalView;
  signalView.SetSensor(&sensor);
  TCanvas* cg = new TCanvas("cg", "", 600, 600);
  signalView.SetCanvas(cg);
  signalView.SetLabelY("signal [fc/ns]]");
  signalView.PlotSignal("wf");
  //电流信号输出
  //打印原始总电流，电子电流，离子电流信号到txt文件
  printsignaltotxt(sensor,"wf",nbins,tstep,tmin);   
  
  // Convolute with the transfer function and plot again.感应电压
  sensor.SetTransferFunction(Transferfunction);
  constexpr bool fft = true;
  sensor.ConvoluteSignals(fft);
  TCanvas* cv = new TCanvas("cv", "", 600, 600);
  signalView.SetCanvas(cv);
  signalView.SetLabelY("signal [mV]");
  signalView.PlotSignal("wf"); 
  //导出电压信号
  printvol(sensor,"wf",nbins,tstep,tmin);  

  std::cout << "Calculation is finished.\n";    //计算完成标志

  app.Run(true);

}