#include <iostream>
#include <TROOT.h>
#include <TApplication.h>
#include <math.h>
#include <bits/stdc++.h>
#include <omp.h>
#include <vector>
#include <fstream>
#include <stdio.h> // for log files

#include "Garfield/MediumMagboltz.hh"
#include "Garfield/ComponentAnalyticField.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/TrackHeed.hh"
#include "Garfield/TrackSrim.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/FundamentalConstants.hh"
#include "Garfield/Random.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/ViewCell.hh"
using namespace Garfield;

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

  //freopen("init_track_log.txt","w",stdout);

  TApplication app("app", &argc, argv);
  
  randomEngine.Seed(245);
  
  time_t start,end;
  double time_taken;
  time(&start);
  
  // Wire radii [cm].
  // Active wire
  const double rActive = .00254; //1 mil
  // Guard wire
  const double rGuard = .0127; //5 mil 

  // 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.);
  // Read the ion mobility table from file.
  gas.LoadIonMobility("/home/khang/FPD_Electric_Field/field_shaping_cage/version8_SRIM/isobutane_ion.txt");
  gas.SetMaxElectronEnergy(3000.);
  gas.Initialise(true);

  // Setup the electric field.
  ComponentAnalyticField cmp;
  cmp.SetMedium(&gas);

  // Calculate the current of the field shaping cage.
  // Note that the cathode plate is also part of this
  // circuit. There are a total of 12 10MOhm resistors.
  double Vcathode = -700.;
  double R = 10.e6;
  double I = Vcathode/(12*R);
  double vdrop = I*R;

  // Adding the cathode plate.
  double vloop = Vcathode - vdrop;
  double ypos = 0.;
  cmp.AddPlaneY(ypos, vloop, "cathode_plate");

  // Adding ten rectangular wire 'loops'.
  constexpr double x = 3.81; // the depth is 3 in. (~7.62 cm) centered at 0
  const double inch = 2.54; // cm
  const double loop_gap = 0.136*inch;
  double diameter_active = 2.*rActive;
  for (int i = 0; i < 10; i++)
  {
	vloop -= vdrop;
	ypos += loop_gap;
	cmp.AddWire(-x, ypos, diameter_active, vloop, "rl"); // left side
	cmp.AddWire( x, ypos, diameter_active, vloop, "rl"); // right side
  }

  //Adding Frisch grid
  int NofWires = 28;
  double left_edge = -x + 0.2;
  double right_edge = x - 0.2;
  double spacing = (right_edge - left_edge) / (NofWires - 1);
  vloop -= vdrop;
  ypos += loop_gap;
  for (int i = 0; i < NofWires; i++) 
  {
    cmp.AddWire(left_edge + i*spacing, ypos, diameter_active, vloop, "f");
  }
  
  //Adding the anode wires and the guard wires
  //Left side
  constexpr double vG = 1000.;
  constexpr double vA = 1400.;
  const double center_left = -0.5*x;
  const double center_right = 0.5*x;
  const double wire_gap = .25*inch;
  ypos += .375*inch;
  cmp.AddWire(center_left - wire_gap, ypos, diameter_active, vA, "al-1");
  cmp.AddWire(center_left     			, ypos, diameter_active, vA, "al");
  cmp.AddWire(center_left + wire_gap, ypos, diameter_active, vA, "al+1");
  //Right side
  cmp.AddWire(center_right - wire_gap, ypos, diameter_active, vA, "ar-1");
  cmp.AddWire(center_right     			 , ypos, diameter_active, vA, "ar");
  cmp.AddWire(center_right + wire_gap, ypos, diameter_active, vA, "ar+1");
  
  //"Square" blocks on the sides of anode wires
  const double dBlock = .25*2.54;
  cmp.AddWire( -x + 0.5*dBlock, ypos, dBlock, 0., "Block");
  cmp.AddWire(-0.5*dBlock - .1, ypos, dBlock, 0., "Block");
  cmp.AddWire( 0.5*dBlock + .1, ypos, dBlock, 0., "Block");
  cmp.AddWire(  x - 0.5*dBlock, ypos, dBlock, 0., "Block");
  
  //Adding neutral plane above the wires
  ypos += .125*inch + .15; // the .15 accounts for the plastic washers b/w supports and pickup pads
  cmp.AddPlaneY(ypos, 0., "top_plane");

  // Make a sensor.
  Sensor sensor;
  sensor.SetArea(-x, 0, -40., x, ypos + 0.2, 40.); //setting drift region
  sensor.AddComponent(&cmp);

  // Set the time bins.
  const unsigned int nTimeBins = 10000;
  const double tmin =  0.;
  const double tmax = 10000.;
  const double tstep = (tmax - tmin) / nTimeBins;
  sensor.SetTimeWindow(tmin, tstep, nTimeBins);

	// 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;
  }
  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 (not sure what's the correct mixing law).
  const double za = 34. / 58.;
  track.SetAtomicMassNumbers(34. / za, 34);
  track.SetDensity(4.769e-4); // derived using ideal gas law with 150 torr and 293.15K

  double xc = 0., yc = 0., zc = 0., tc = 0., ec = 0., extra = 0.;
  int nc = 0;
  std::vector<double> electrons;
  double xe = 0., ye = 0., ze = 0., te = 0., ee = 0.;
  double dx = 0., dy = 0., dz = 0.;
  int index = 0;
  int left=0, right = 0;
  int NofProton = 10;
  int noe = 0;
  for (int i = 0; i < NofProton; i++) {
    // Simulate a track
    // x-dir goes along the width of the detector (i.e. perpendicular to the PSAs)
    // y-dir goes along the height of the detector
    // z-dir goes along the length of the detector (i.e. parallel to the PSAs)
    //track.NewTrack(x, 1.0, 0.6, 100., -1., 0., 0.); //90 degrees
    track.NewTrack(-x, 1.0, -x, 100., 1., 0., 1.); //45 degrees
    //track.NewTrack(-x, 1.0, -(x+center_right)/.900404040443, 100., .6691290616, 0., .74314311); //42 degrees
    noe = 0;
		while (track.GetCluster(xc, yc, zc, tc, nc, ec, extra)) {
			for (int k = 0; k < nc; ++k) {	
			  electrons.push_back(xc);
			  electrons.push_back(yc);
			  electrons.push_back(zc);
			  electrons.push_back(tc);
			  ++index;
			  ++noe;
			}
		}
		std::cout << "Proton number: " << i << std::endl;
		std::cout << "Number of electrons: " << noe << std::endl;
  }
	
  // Sanity check to make sure that the electrons array has the correct size
  // std::cout << "index " << index << std::endl;
  // std::cout << "size " << electrons.size() << std::endl;
  int size = electrons.size()/4;
  if (size != index) {
  	std::cout << "Something is wrong with electrons!!" << std::endl;
  	return 0;
  }

  //total output for testing
  std::ofstream total("totalfile.txt", std::ios::out);
  for (int j = 0; j < index; j++) {
    total << electrons.at(j*4 + 0) << "\t" << electrons.at(j*4 + 1) << "\t" <<
               electrons.at(j*4 + 2) << "\t" << electrons.at(j*4 + 3) << "\n";
  }
  total.close();  

  time(&end);
  time_taken = double(end-start);
  std::cout << "Time taken by program is: " << std::fixed << time_taken << std::setprecision(5) << " sec" << std::endl;

	//fclose(stdout);  
  //app.Run(true);
  
}