#include <iostream>
#include <cmath>

#include <TROOT.h>
#include <TApplication.h>
#include <TH1F.h>
#include <TFile.h>

#include "Garfield/MediumMagboltz.hh"
#include "Garfield/ComponentConstant.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/AvalancheMicroscopic.hh"

using namespace Garfield;

int main() {

  TFile outfile("arco2.root", "RECREATE");

  // Electric field [V / cm].
  constexpr double field = 570. * 40; 
  constexpr double gap = 0.025;
  // Initial electron energy [eV].
  constexpr double e0 = 1.;

  // Make a gas medium.
  MediumMagboltz gas("ar", 90., "co2", 10.);
  gas.SetTemperature(293.15);
  gas.SetPressure(760.);
  gas.SetMaxElectronEnergy(1000.);
  gas.Initialise();

  // Make a component with constant drift field.
  ComponentConstant cmp;
  cmp.SetMedium(&gas);
  cmp.SetArea(-2., -2., -gap, 2., 2., gap);
  cmp.SetElectricField(0, 0, field);

  // Make a sensor.
  Sensor sensor;
  sensor.AddComponent(&cmp);
  
  // Microscopic tracking.
  AvalancheMicroscopic aval;
  aval.SetSensor(&sensor);
  aval.SetCollisionSteps(10);
    
  // Histograms
  TH1::StatOverflows(true);
  TH1F hElectrons("hElectrons", "Avalanche Size", 200, 0, 200);
  TH1F hIons("hIons", "Avalanche Size", 200, 0, 200);

  constexpr unsigned int nEvents = 10000;
  for (unsigned int j = 0; j < nEvents; ++j) {
    aval.AvalancheElectron(0., 0., 0., 0., e0);
    int ne, ni;
    aval.GetAvalancheSize(ne, ni);
    if (j % 100 == 0) { 
      std::cout << j << "/" << nEvents << "\n"
                << "    " << ne << " electrons\n";
	  }
    hElectrons.Fill(ne);
    hIons.Fill(ni);
  }
  outfile.cd();
  hElectrons.Write();
  hIons.Write();

  const double mean = hElectrons.GetMean();
  const double rms = hElectrons.GetRMS();
  std::cout << "mean = " << mean 
            << ", f = " << rms * rms / (mean * mean) << "\n";
  outfile.Close();
}