//
// Created by Dario Stocco (stoccod@ethz.ch) on 02.08.2023.
//
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <numeric>

#include <TApplication.h>
#include <TCanvas.h>
#include <TH1F.h>
#include <TSystem.h>

#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/AvalancheGridSpaceCharge.hh"
#include "Garfield/ComponentParallelPlate.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/Plotting.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/ViewSignal.hh"
#include "Garfield/TrackHeed.hh"

using namespace Garfield;

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

   TApplication app("app", &argc, argv);
   plottingEngine.SetDefaultStyle();

   double voltage = 2000;

   MediumMagboltz gas;
   gas.EnableAutoEnergyLimit(false);
   gas.SetMaxElectronEnergy(100.);
   gas.LoadGasFile("c2h2f4_94_sf6_1_ic4h10_5.gas");
   gas.Initialise(true);

   // Dimensions of RPC (cm)
   double d_dlc = 0.00001;    // 100 nm
   double d_pol = 0.005;      // 50 um
   double d_gas = 0.0375;     // 375 um
   std::vector<double> layers = {d_pol,
      d_dlc,
      d_gas,
      d_dlc,
      d_pol};
   double y_mid = d_pol + d_dlc + d_gas / 2;

   // Relative permittivity
   double e_dlc = 4.;
   double e_pol = 3.3;
   double e_gas = 1.;
   std::vector<double> eps = {e_pol,
      e_dlc,
      e_gas,
      e_pol,
      e_dlc};

   // ComponentParallelPlate
   ComponentParallelPlate cmp;
   cmp.Setup(int(layers.size()), eps, layers, voltage, {});
   std::string label = "readout";
   cmp.AddPlane(label);
   cmp.SetMedium(&gas);

   // Sensor
   Sensor sens;
   sens.AddComponent(&cmp);
   sens.AddElectrode(&cmp, label);
   sens.SetTimeWindow(0, (5. - 0) / 200., 200);

   // AvalancheMicroscopic
   AvalancheMicroscopic aval;
   aval.SetSensor(&sens);

   // AvalancheGridSpace Charge
   AvalancheGridSpaceCharge avalsc;
   avalsc.EnableDiffusion(true);
   avalsc.EnableStickyAnode(false);
   avalsc.EnableAdaptiveTimeStepping(true);
   avalsc.SetStopAtK(true);
   avalsc.EnableSpaceChargeEffect(false);
   avalsc.SetSensor(&sens);
   avalsc.Set2dGrid(y_mid - 0.5 * d_gas + 1.e-8,
                    y_mid + 0.5 * d_gas - 1.e-8, 400, 0.05, 100);
   //avalsc.AvalancheElectron(xi, yi, zi, ti, nEvents);
   //avalsc.StartGridAvalanche();
   avalsc.ImportElectronsFromAvalancheMicroscopic(&aval);

   // Set the electron start parameters.
   const double zi = 0., xi = 0., yi = y_mid - 0.99 * d_gas / 2;
   const double ti = 0., ei = 0.;

   // Avalanche electron
   TH1D *hist = new TH1D("energy", Form("Avalanche energy at %.0f V", voltage), 200, -1, 50);
   constexpr unsigned int nEvents = 1000;
   for ( unsigned int i = 0; i < nEvents; i++) {
      if (i%100 == 0)   std::cout << i << "/" << nEvents << std::endl;
      aval.AvalancheElectron(xi, yi, zi, ti, ei, 0., 0., 0.);
      for ( const auto &electron : aval.GetElectrons() ) {
         const auto &p = electron.path[1];
         hist->Fill(p.energy);
      }
   }

   TCanvas *canv = new TCanvas("canv", "canv", 700, 600);
   canv->SetGrid();
   hist->Draw();

   // view recorded signals from plane electrode
   const bool signal = false;
   if (signal) {
      ViewSignal *signal_view = new ViewSignal(&sens);
      TCanvas *c_signal = new TCanvas(label.c_str(), label.c_str(), 600, 600);
      signal_view->SetCanvas(c_signal);
      signal_view->PlotSignal(label);
      c_signal->SetTitle(label.c_str());
      gSystem->ProcessEvents();
   }

   std::cout << "Done" << std::endl;

   app.Run(true);
   return 0;
}