#include <iostream>
#include <vector>

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

#include "Garfield/MediumMagboltz.hh"
#include "Garfield/ComponentAnalyticField.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/TrackHeed.hh"
#include "Garfield/TrackSrim.hh"
#include "Garfield/FundamentalConstants.hh"
#include "Garfield/Random.hh"
#include "Garfield/ViewCell.hh"

using namespace Garfield;

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

  TApplication app("app", &argc, argv);
  
  randomEngine.Seed(245);
  
  // Setup the gas.
  MediumMagboltz gas;
  // Set the temperature [K] and pressure [Torr].
  gas.SetTemperature(293.15);
  gas.SetPressure(150.);
  // Set the gas mixture.
  gas.SetComposition("isobutane", 100.);

  constexpr double lx = 3.81;
  constexpr double ly = 5.;
  // Cell layout.
  ComponentAnalyticField cmp;
  cmp.SetMedium(&gas);
  cmp.AddPlaneY(0., -700., "cathode_plate");
  cmp.AddPlaneY(ly,    0., "top_plane");

  // Make a sensor.
  Sensor sensor;
  sensor.AddComponent(&cmp);
  sensor.SetArea(-lx, 0, -40., lx, ly, 40.); 
	// Create a track class.
  TrackSrim track;
  // Connect the track to a sensor.
  track.SetSensor(&sensor);
  // Read SRIM output from file.
  const std::string file = "Hydrogen_in_IsoButane_gas.txt";
  track.ReadFile(file);
  if (!track.ReadFile(file)) {
    std::cerr << "Reading SRIM file failed.\n";
    return 0;
  }
  track.SetTargetClusterSize(10);
  //track.SetClustersMaximum(3000);
  // Set the initial kinetic energy of the particle (in eV).
  track.SetKineticEnergy(10.e6);
  // Set the W value and Fano factor of the gas.
  track.SetWorkFunction(26.52);
  track.SetFanoFactor(0.27);
  // Set A and Z of the gas.
  track.SetAtomicMassNumbers(58, 34);
  track.SetDensity(gas.GetMassDensity());

  int NofProton = 10;
  for (int i = 0; i < NofProton; i++) {
    track.NewTrack(lx, 1.0, 0., 0., -1., 0., 0.);
    int nsum = 0;
    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)) {
      nsum += nc;
		}
		std::cout << "Proton number: " << i << std::endl;
		std::cout << "Number of electrons: " << nsum << std::endl;
  }
	
  //app.Run(true);
  
}