#include <iostream>
#include <sstream>
#include <string>
#include <sys/time.h>
#include <stdlib.h>
#include <stdio.h>
#include <TFile.h>
#include <TH1.h>
#include <TMath.h>
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/SolidBox.hh"
#include "Garfield/ComponentAnalyticField.hh"
#include "Garfield/ComponentGrid.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/GeometrySimple.hh"
#include "Garfield/DriftLineRKF.hh"
#include <TRandom3.h>
#include <TROOT.h>
#include <TRint.h>
#include <fstream>
#include <sys/types.h>
#include <unistd.h>
#include <fstream>
#include <vector>
#include <chrono>
#include "Garfield/ViewIsochrons.hh"

using namespace Garfield;

int main(int argc, char * argv[]) {
  auto start = std::chrono::high_resolution_clock::now();
  int pid = getpid();
  timeval t;
  gettimeofday(&t, NULL);
  int seed = pid*t.tv_usec;
  std::cout << "Random Seed: " << seed << std::endl;
  gRandom->TRandom::GetSeed();
  MediumMagboltz * gas = new MediumMagboltz();
  // Setup the gas
  const double pressure = 14.616*51.7149; // in torr- we were at 14.6 psi (1 psi = 51.7149 torr) 
  const double temperature = 299.15;
  gas->SetTemperature(temperature); // from CLI
  gas->SetPressure(pressure);
  gas->SetComposition("CO2", 90.,"AR", 10.);
  const int nFields = 5;
  const double E_not = 984.25;
  const double emin = E_not-E_not;
  const double emax = E_not+E_not;
  // Flag to request logarithmic spacing.
  const bool useLog = false;
  const double bmin=0;
  const double bmax=0; // do we need magnetic field on?
  const int nBFields=1;
  gas->SetFieldGrid(emin, emax, nFields, useLog, bmin,bmax,nBFields,TMath::Pi()/2.0,TMath::Pi()/2.0,1); 
  // Turn on penning transfer?
  gas->EnablePenningTransfer();
  gas->SetMaxElectronEnergy(200);
  std::cout << "number of levels: " << gas->GetNumberOfLevels();
  const int ncoll = 5;
  gas->GenerateGasTable(ncoll);
  char * IonData = getenv("GARFIELD_IONDATA") ;
  gas->LoadIonMobility(IonData);
  gas->PrintGas();
  ComponentAnalyticField * cmpE = new ComponentAnalyticField();
  GeometrySimple * geo = new GeometrySimple();
  SolidBox * enclosure = new SolidBox(0,0,0,5,5,5);
  geo->AddSolid(enclosure, gas);
  cmpE->SetGeometry(geo);
  const double vAnode= 0;
  const double rAnode= 20e-4;
  cmpE->AddWire(0,0,2 * rAnode, vAnode, "a");
  cmpE->AddPlaneX(-7,-7500,"cP1"); 
  cmpE->AddPlaneX(7,-7500,"cP2");

  // Load the field map.
  ComponentGrid * cmpB = new ComponentGrid();
  cmpB->SetGeometry(geo);
  cmpB->LoadMagneticField("garfield_randomB.txt", "XYZ"); // come up with a file that has x,y,z in cm and bx,by,bz in Tesla
  // add both cmpB and cmpE to sensor object which is given to DriftLineRKF
  Sensor * sensor = new Sensor;
  sensor->AddComponent(cmpE);
  sensor->AddComponent(cmpB);
  sensor->SetTimeWindow(0,2,20000); // might need to change this, its start, step size, number of steps
  cmpE->AddReadout("a");
  sensor->AddElectrode(cmpE,"a");
  sensor->EnableComponent(0, true);
  sensor->EnableComponent(1, true);

  std::vector<std::array<double, 3> > points;
  unsigned int nPoints = 20;
  const double ymin=-2.0;
  const double ymax=2.0;
  double start_drift =-2;
  for (unsigned int i = 0; i < nPoints; ++i) {
    const double x0 = start_drift;
    const double y0 = ymin + i * (ymax - ymin) / nPoints;
    std::array<double, 3> p0 = {x0, y0, 0.};
    points.push_back(std::move(p0));
  } 
  std::vector<int> statusCodes;
  // steps when running through ComputeDriftLines
  DriftLineRKF drift; 
  drift.SetSensor(sensor);
  drift.SetMaximumStepSize(); // need to think about this potentially
  drift.EnableSignalCalculation(false);
// everything else is one big for loop around the starting points
  for (const auto& point : points) {
    // always use DriftElectron 
    drift.DriftElectron(point[0], point[1], point[2], 0.);
    const unsigned int nu = drift.GetNumberOfDriftLinePoints();
    // Check that the drift line has enough points.
    if (nu < 3) continue;
    int status = 0;
    double xf = 0., yf = 0., zf = 0., tf = 0.;
    drift.GetEndPoint(xf, yf, zf, tf, status);
    statusCodes.push_back(status);
  }
  std::cout << "Size Statuses: " << statusCodes.size() << std::endl;
  std::cout << "And now min and max status codes are: ";
  auto minmax_it = std::minmax_element(statusCodes.begin(), statusCodes.end());
  int min_value = *minmax_it.first;
  int max_value = *minmax_it.second;
  double step = 0.5; // for ints..
  int nbins=std::round(((double)(max_value)-(double)(min_value))/(step/1.0))+1;
  double low_edge=min_value-(step/2.0);
  double high_edge=max_value+(step/2.0);
  TH1F * status_H = new TH1F("statusCodes","statusCodes", nbins,low_edge,high_edge);  
  for (int thisstatus : statusCodes) {
      status_H->Fill(thisstatus);
  }
  auto end = std::chrono::high_resolution_clock::now();
  std::chrono::duration<double> elapsed_seconds = end - start;
  std::cout << "Elapsed time: " << elapsed_seconds.count() << " seconds\n";
  std::stringstream outrootfilename;
  //if(BFieldValue==0.0)  outrootfilename << "DriftLineAllwires4800_" << std::fixed << std::setprecision(1) << BFieldValue << "T_" << temperature << "K_" << input_tstep << "ns_wakely.root";
  outrootfilename << "StatusCodesMRE_output.root";
  TFile * Outfile = new TFile(outrootfilename.str().c_str(),"recreate");
  Outfile->cd();
  status_H->Write();
  Outfile->Close();
  return 0;
}