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

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

#include "Garfield/ViewCell.hh"
#include "Garfield/ViewDrift.hh"
#include "Garfield/ViewSignal.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/SolidBox.hh"
#include "Garfield/GeometrySimple.hh"
#include "Garfield/ComponentConstant.hh"
#include "Garfield/TrackSrim.hh"
#include "Garfield/Random.hh"
#include "Garfield/Plotting.hh"


using namespace Garfield;

//电流前端电子学卷积函数
double Transferfunction(double t){     
  constexpr double risistance= 100 ;   //电阻设置
  constexpr double capacitance=  100;    //电容设置
  constexpr double tau=risistance*capacitance;
  return exp(t/tau);
}

//电流视为delat冲击电流的响应函数
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;
}

//入射离子产生的cluster的位置，时刻沉积能量输出到txt文件
void printinformationofclusters(std::vector<std::array<double,7>> informationofclusters,\
     float name[]){     
  char title[60];
  sprintf(title,"cluster information Ar_%.1f_CO2_%.1f_300K_%.2fAtm_%.0fV.txt",name[0],name[1],name[2],name[3]);
  std::ofstream outfile;
  outfile.open(title,std::ios::out);
  outfile <<  "the informations of every cluster in a track of projectile:\n";
  outfile <<  "order    X(cm)   Y(cm)    Z(cm)     moment(ns)     number of electrons     desposited energy(eV)   "
      "    enenrgy before collision(Mev):\n";
  const std::size_t numberofcluster=informationofclusters.size();
  for(size_t i=0;i<numberofcluster;++i){
  outfile <<i+1 <<"   "<<informationofclusters[i][0] << "   " <<informationofclusters[i][1] <<"   "\
                         <<informationofclusters[i][2] << "   " <<informationofclusters[i][3] <<"   "\
                         <<informationofclusters[i][4] << "   " <<informationofclusters[i][5] <<"   "\
                         <<informationofclusters[i][6] << "\n";
  }
  outfile.close();
}

//cluster的电子和离子倍增数,位置，时刻沉积能量及产生的输出到txt文件表头打印
void printtitleofavalancheinclusters(float name[]){ 
  char title[80];
  sprintf(title,"infoofavalancheinclusters Ar_%.1f_CO2_%.1f_300K_%.2fAtm_%.0fV.txt",name[0],name[1],name[2],name[3]);
  std::ofstream outfile;
  outfile.open(title,std::ios::out);
  outfile <<  "the avarange electron and ion in each cluster:\n";
  outfile <<  "order    X(cm)     Y(cm)    Z(cm)    moment(ns)     enenrgy before collision(Mev)"
              "    avalancheelectrons    avalancheions:\n";
  outfile.close();
}
//cluster的电子和离子倍增数,位置，时刻沉积能量及产生的输出到txt文件
void printavalancheinclusters(int j,double positionx,double positiony,double positionz,\
                              double moment,double enenrgybeforecollision,\
                              double avalancheelectrons,double avalancheions,float name[]){     
  char title[80];
  sprintf(title,"infoofavalancheinclusters Ar_%.1f_CO2_%.1f_300K_%.2fAtm_%.0fV.txt",name[0],name[1],name[2],name[3]);
  std::ofstream outfile;
  outfile.open(title,std::ios::app);   //追加打印
  outfile <<j+1 <<"    "<<positionx<<"    "<<positiony<<"    "<<positionz<<"    "<<moment<<"    "\
                <<enenrgybeforecollision<<"    "<<avalancheelectrons<<"    "<<avalancheions<<"\n";
  outfile.close();
}

//原始电流信号打印到txt文件函数
void printsignaltotxt(Sensor& sensor,const std::string label,const int nTimebins,const double tstep,const double tmin,float name[]){     
  char title[60];
  sprintf(title,"currentsignal Ar_%.1f_CO2_%.1f_300K_%.2fAtm_%.0fV.txt",name[0],name[1],name[2],name[3]); 
  std::ofstream outfile;
  outfile.open(title,std::ios::out);
  outfile << "time     "<< "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";
  }
}

//卷积计算信号信号打印到txt文件函数
void printconvolutesignaltotxt(Sensor& sensor,const std::string label,const int nTimebins,const double tstep){     
  std::ofstream outfile;
  outfile.open("convolutionsignal.txt",std::ios::out);
  for(unsigned int i=0;i<nTimebins;++i){
    const double t=(i-1)*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);
  
  float name[4];
  // ---创建气体介质
  double ratioofAr=0.85;              //**********
  double ratioofCO2=1-ratioofAr;
  name[0]=ratioofAr*100;
  name[1]=ratioofCO2*100;
  
  MediumMagboltz gas;
  const double pressure = 0.26;   //**********
  name[2]=pressure;
  const double temperature = 300;                     //**********
  gas.SetComposition("Ar", ratioofAr*100, "CO2", ratioofCO2*100);          //     自行设置气体组成压强和温度用于Initialise（）时
  gas.SetTemperature(temperature);                   //     可不用加载气体中电子输运参数文件
  gas.SetPressure(pressure*AtmosphericPressure);                        
  
  gas.EnableAutoEnergyLimit(true);
   gas.SetMaxElectronEnergy(600);
  gas.EnableCrossSectionOutput();    //打印一个Cs文件包含碰撞截面等信息
  gas.Initialise(true);  
  
  const std::string path = std::getenv("GARFIELD_HOME");
  gas.LoadIonMobility(path + "/Data/IonMobility_Ar+_Ar.txt");  //设置气体中离子输运参数

  // ---使用解析场构建Component
  ComponentAnalyticField cmp;
  cmp.SetMedium(&gas);
  const double rWire = 15.e-4;    //**********  //中心阳极丝半径 （cm）
  const double rTube = 3.5;      //**********  //阴极外壳半径   （cm）
  // 电压设置
  const double vWire = 1100;      //********** //阳极丝电压
  name[3]=vWire;
  const double vTube = 0.;        //**********  //阴极外壳电压
  // Add the wire in the centre.
  cmp.AddWire(0, 0, 2 * rWire, vWire, "s");//cmp.AddWire(0[X位置], 0[Y位置],半径, vWire[电压], "s"[标记]，长度‘可选’，张力‘可选，单位为g’，密度‘可选’);
  // Add the tube.
  cmp.AddTube(rTube, vTube, 0, "t");//第三个参数控制管的多边形，0为圆柱，可选3-8
  // Request calculation of the weighting field. 
  cmp.AddReadout("s");         //信号从标记‘s’处读出

  // ---make a sensor.
  Sensor sensor;
  sensor.AddComponent(&cmp);      //与组成Component“cmp”链接
  sensor.AddElectrode(&cmp, "s");     //以“cmp”为权电场计算标记为“s”上的信号
  // Set the signal time window.//设置信号计算的时间下时间限，步长及时间步长个数，单位ns
  const double tstep = 0.2;      //**********  //信号采样步长 （ns）
  const double tmin = 0;     //********** //信号采样时间下限
  const double tmax=10000;         //**********      
  const unsigned int nbins = (tmax-tmin)/tstep;   //信号采样时间间隔数 
  sensor.SetTimeWindow(tmin, tstep, nbins);
  sensor.SetTransferFunction(Transferfunction);//设置电流的响应函数
  sensor.ClearSignal();

  //---设置入射粒子径迹计算模块
  TrackSrim tr;
  tr.SetSensor(&sensor); // 与sensor链接.
  const std::string file = "Li_in_Ar_CO2_gas.txt";   //**********  //读取入射离子的srim计算文件
  if (!tr.ReadFile(file)) {
    std::cerr << "Reading SRIM file failed.\n";
  }

  tr.SetKineticEnergy(0.84e6);     //**********   //设置动能（eV）  alpha 1.47e6 eV   Li+  0.84e6 eV
  tr.SetWorkFunction(27.7);      //**********   // 设置材料功函数  
  tr.SetFanoFactor(0.172);      //********** // 设置法诺因子
  const double za = ratioofAr * (18. / 39.948) +ratioofCO2 * (22. / 44.0095);  // Set A and Z of the gas (not sure what's the correct mixing law).
  tr.SetAtomicMassNumbers(22. / za, 22);
  //tr.SetTargetClusterSize(500);     //**********  //设置多少个电子并为一簇   默认一个径迹产生约100个cluster，
  //int numbersofcluster=100;            //**********
  //tr.SetClustersMaximum(2*numbersofcluster);    //**********  //指定平均产生n/2个cluster    与tr.SetTargetClusterSize(500)二选一式使用
  tr.EnableTransverseStraggling(false);
  tr.EnableLongitudinalStraggling(false);

  // 绘制加载的srim文件信息
  //tr.PlotEnergyLoss();
  //tr.PlotRange();
  //tr.PlotStraggling(); 
  tr.Print();  //打印SRim计算文件的相关信息

  // 设置统计直方图
  TH1F* hX = new TH1F("hX", "x-end", 500, -10, 10);  //0-0.01为横坐标区间
  TH1F* hY = new TH1F("hY", "y-end", 500, -10, 10);
  TH1F* hZ = new TH1F("hZ", "z-end", 500, -10, 10);
  TH1F* hNe = new TH1F("hNe", "Electrons", 100, 1000, 50000);
 
  TCanvas* cD = nullptr;     //定义一个指向TCanvas的空指针以备用
  ViewCell cellView;   //  定义探测器结构绘制类
  ViewDrift driftView;  // 定义漂移线绘制类
  constexpr bool plotDrift = false;   //定义漂移线绘图控制开关
  if (plotDrift) {
    cD = new TCanvas("cD", "", 600, 600);  //定义新绘图界面
    cellView.SetCanvas(cD);             //设置探测器绘制在cD图层中   
    cellView.SetComponent(&cmp);        //设置指定所要绘制的探测器     
    driftView.SetCanvas(cD);            //设置探测器绘制在cD图层中
    tr.EnablePlotting(&driftView);   //  入射粒子轨迹绘制数据传递（TrackSRIM中控制项）     
  }
 
  TCanvas* cS = nullptr;
  ViewSignal signalView;
  constexpr bool plotSignal = true;    //电流信号绘制控制开关
  if (plotSignal) {
    cS = new TCanvas("cS", "", 600, 600);   //定义新绘图界面
    signalView.SetCanvas(cS);                //设置绘制信号在cS图层
    signalView.SetSensor(&sensor);            //链接sensor
    signalView.SetLabelY("signal uA");       //设置y轴标签为"signal [fC]"
  } 
  
  // ---Generate tracks.
  //  定义初始位置
  const double rTrack =0.3;     //**********
  const double projectilex= rTrack;
  const double projectiley= -sqrt(rTube * rTube - rTrack * rTrack);
  const double projectilez= 0.;  //**********
  //  定义初始方向
  const double projectiledx=0.;    //**********
  const double projectiledy=1.; //**********
  const double projectiledz=0.; //**********
  //  定义入射时刻及总的粒子数
  const double projectileT=0.;   //**********
  const unsigned int nTracks =1;    //**********    //nTracks设为入射粒子个数
  //std::vector<std::array<double,4>> primaryelectrons;
  //std::vector<std::array<double,4>> avalancheions;
  //std::vector<std::array<double,4>> primaryions;
  for (unsigned int n = 0; n < nTracks; ++n) {
    //对一个入射离子的径迹进行计算，得到多个clusters //设置入射粒子的位置方向，前三个参量为位置，后三个参量为方向
    const bool notnewtrack= !tr.NewTrack(projectilex, projectiley, projectilez, projectileT, projectiledx, projectiledy, projectiledz);    
    if (notnewtrack) {
      std::cerr << "Generating clusters failed; skipping this particle's track.\n";
      continue;
    }
    //对单个入射粒子计算产生的clusters信息
    unsigned int netot = 0;
    std::vector<std::array<double,7>> informationofclusters;
    while (true) {                          
      int ne = 0;
      double xc, yc, zc, tc, ec, ekin;
      const bool done = tr.GetCluster(xc, yc, zc, tc, ne, ec, ekin);  //GetCluster()为得到下一个cluster的信息：xc, yc, zc为cluster坐标，
      if(done) {                                                                  //tc为时刻，ne为单次碰撞产生的电子数，ec为单次碰撞中能量损失，ekin为入射离子此时此位置的动能
      informationofclusters.push_back({xc, yc, zc, tc, ne, ec,ekin});          //将每一个cluster的产生位置，时刻及沉积能量传到informationofclusters向量变量
      netot += ne;        //计算单个入射粒子产生的总的电子数
      }   
      if (!done) {
        hX->Fill(xc);   //每一个入射离子最后一个cluster中的电子X位置
        hY->Fill(yc);   //每一个入射离子最后一个cluster中的电子Y位置
        hZ->Fill(zc);   //每一个入射离子最后一个cluster中的电子Z位置
        hNe->Fill(netot);    //每一个入射离子产生的总的电子数
        std::cout  << informationofclusters.size() <<  " clusters infornation is passed to file.\n\n";
        break;                    //跳出while循环的控制
      }
    }
 
    //对cluster产生的时刻进行修正
    double MLi=6.9405*1.66053886e-27;     //**********  //Li离子的原子量换算为kg
    double Mev=1.602176565e-13;              //1MeV等量焦耳能
    for(int i=0;i<informationofclusters.size()-1;++i){
      double speedofprojectile=sqrt(2*informationofclusters[i+1][6]*Mev/MLi)/1.e7;    //计算那cluster间的输运速度并转换单位为cm/ns
      double distanceofclusters=sqrt(pow(informationofclusters[i+1][0]-informationofclusters[i][0],2) \
                                    +pow(informationofclusters[i+1][1]-informationofclusters[i][1],2) \
                                    +pow(informationofclusters[i+1][2]-informationofclusters[i][2],2));  //计算两相邻cluster间距离（cm）
      double duation=distanceofclusters/speedofprojectile;                  //计算入射离子在两cluster间的输运时间
      informationofclusters[i+1][3]=informationofclusters[i][3]+duation;     //在第i个cluster时刻的基础上加上输运时间对第i+1个cluster的产生时刻进行修正
    }
    printinformationofclusters(informationofclusters,name);  //将修正完成的入射离子的位置，时刻沉积能量，碰撞前能量等信息打印到txt文件
    printtitleofavalancheinclusters(name);//打印倍增信息表头
   //拉出informationofclusters向量数组的cluster位置，时刻信息进行并行电子/离子对雪崩输运计算
    omp_set_num_threads(55);    //并行计算线程数
    #pragma omp parallel for
      //for(int j=60; j<61; ++j){  // 测试命令
      for(int j=0; j<informationofclusters.size(  ); ++j){  // informationofclusters.size()
        double avalancheelectrons=0;
        double avalancheions=0;   
        AvalancheMicroscopic driftelectron;
        driftelectron.SetSensor(&sensor);
        driftelectron.EnableSignalCalculation();      //开启AvalancheMicroscopic漂移模拟的电流信号计算
        driftelectron.EnablePlotting(&driftView);   //电子漂移线绘制数据传递 (DriftLineRkF中控制项，会将指向一个ViewDrift类的指针‘&driftView’值传递给一个名为m_iew的变量供后续操作)
        
        AvalancheMC driftion;           //设置一个AvalancheMC用以计算初电离离子漂移
        driftion.SetSensor(&sensor);     //链接Sensor
        driftion.SetDistanceSteps(2.e-4);   //设置蒙卡计算离子漂移时的步长    
        driftion.EnableSignalCalculation();    //开启AvalancheMC漂移模拟的电流信号计算

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

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

        //for (int k = 0; k <5; ++k) {            //命令测试
        for (int k = 0; k <informationofclusters[j][4]; ++k) {            //  此循环为对一个cluster中的每一个电子离子对输运  
          std::cout << k+1<< "/" << informationofclusters[j][4] <<  " Transport of electron in " \
                    << j+1<< "/" << informationofclusters.size() << " cluster begin:\n"; 
          std::cout <<  "AvalancheElectron begin;                     "; 
          //#pragma omp critical
          {
          driftelectron.AvalancheElectron(informationofclusters[j][0], informationofclusters[j][1], \
                                          informationofclusters[j][2], informationofclusters[j][3], 0.00001, 0., 0., 0.);    //对(xe, ye, ze, te)位置处电子进行的漂移雪崩进行计算
          }
          int avalanchee,avalanchei;      
          driftelectron.GetAvalancheSize(avalanchee,avalanchei);   
          std::cout << "numberofavalancheelectron is " <<avalanchee <<"        " << "numberofavalancheion is " <<avalanchei <<"\n";     
          avalancheelectrons=avalancheelectrons+avalanchee;
          avalancheions=avalancheions+avalanchei;                         //统计每一个cluster中最终倍增的电子
          std::cout <<  "AvalancheElectron end\n";
          std::cout << k+1<< "/" << informationofclusters[j][4] <<  " Transport of electron in "\
                    << j+1<< "/" << informationofclusters.size() << " cluster is done.\n";  
        }
        //输出电子和离子的雪崩倍增系数
        avalancheelectrons=avalancheelectrons/informationofclusters[j][4];//;
        avalancheions=avalancheions/informationofclusters[j][4];                  //每个cluster中进行倍增数对cluster中电子离子对平均化
        printavalancheinclusters(j,informationofclusters[j][0],informationofclusters[j][1],informationofclusters[j][2],\
                                 informationofclusters[j][3],informationofclusters[j][6],avalancheelectrons,avalancheions,name);
        
        std::cout <<  "Ion tranport in Avalanche begin;             ";
        int np=0;
        np=driftelectron.GetNumberOfElectronEndpoints();  //获取一个cluster中最后输运的电子的漂移线碰撞点数目
        std::cout <<  "number of Ions in last electron drift line    "<< np <<"\n";
        for(int m=0; m<np; ++m){                         //此循环对一个cluster中每一个初电离电子雪崩中产生的离子进行漂移计算
          double xestart, yestart, zestart, testart, eestart;
          double xeend, yeend, zeend, teend, eeend;
          int status;
          //获取一个cluster中最后输运的初电离电子引起的雪崩中各倍增电子的初始位置
          driftelectron.GetElectronEndpoint(m, xestart, yestart, zestart, testart, eestart, xeend, yeend, zeend, teend, eeend, status);
          //#pragma omp critical
          {
          driftion.SetIonSignalScalingFactor(informationofclusters[j][4]);//对雪崩中产生的电子相对应的离子进行输运计算
          driftion.DriftIon(xestart, yestart, zestart, testart);
          }
        } 
        std::cout <<  "Ion tranport in Avalanche end\n";
        //cluster中的初电离离子输运
        //#pragma omp critical
        {
        driftprimaryion.SetIonSignalScalingFactor(informationofclusters[j][4]);
        driftprimaryion.DriftIon(informationofclusters[j][0], informationofclusters[j][1], \
                                 informationofclusters[j][2], informationofclusters[j][3]);       //对一个cluster中初电离中的离子进行输运计算
        }

      }
    
  }

  // 在同一图层中按2*2绘图,绘制每一个入射粒子出射最后一个cluster的位置分布及产生的电子总数的分布
  TCanvas* c = new TCanvas("c", "SRIM", 100, 100, 800, 800);
  c->Divide(2, 2);      //将绘图平面变为2*2的布局绘图
  c->cd(1); hX->Draw();
  c->cd(2); hY->Draw();
  c->cd(3); hZ->Draw();
  c->cd(4); hNe->Draw();
  c->Update();

  //  漂移线绘制
  if (plotDrift) {
    cD->Clear();
    cellView.Plot2d();          //绘制二维探测器结构图
    constexpr bool twod = true;
    constexpr bool drawaxis = false;
    driftView.Plot(twod, drawaxis);           //绘制二维漂移线图
  }
  
  //  信号计算
  printsignaltotxt(sensor,"s",nbins,tstep,tmin,name);      //打印原始总电流，电子电流，离子电流信号到txt文件
  //sensor.ConvoluteSignals();        //开启电流传递函数，进行卷积积分计算
  //printconvolutesignaltotxt(senor,"s",nbins,tstep);      //打印卷积计算后的信号到txt文件
  int nt = 0;                      //对记录越过阈值的信号数目计数器初始化
  //if (!sensor.ComputeThresholdCrossings(-2., "s", nt)) continue;      //设置“S”电极上的阈值，返回nt为越过阈值的信号数目
  if (plotSignal) signalView.PlotSignal("s");     //绘制电极“s”上的信号
  std::cout << "Calculation is finished.\n";    //计算完成标志
  app.Run(kTRUE);
}