#include <iostream>
#include <fstream>
#include <cstdlib>
#include<string>

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

#include "Garfield/ViewCell.hh"
#include "Garfield/ViewDrift.hh"
#include "Garfield/ViewSignal.hh"

#include "Garfield/ComponentAnalyticField.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/TrackHeed.hh"

using namespace Garfield;
using namespace std;

//double transfer(double t) {
//constexpr double tau = 25.;
//return (t / tau) * exp(1 - t / tau);
}



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

  TApplication app("app", &argc, argv);
 
  // Make a gas medium.
  MediumMagboltz gas;
  gas.LoadGasFile("./gascard/ar_96_co2_4_2200v_0224.gas");
  const std::string path = std::getenv("GARFIELD_INSTALL");
  gas.LoadIonMobility(path + "/share/Garfield/Data/IonMobility_Ar+_Ar.txt");

  // Make a component with analytic electric field.
  ComponentAnalyticField cmp;
  cmp.SetMedium(&gas);
  // Wire radius [cm]
  const double rWire = 25.e-4;
 //const double rWire = 20.e-4;
  // Outer radius of the tube [cm]
  const double rTube = 1.46;
  //const double rTube = 2.46;
  // Voltages
  const double vWire = 2950.;
  const double vTube = 0.;
  // Add the wire in the centre.
  cmp.AddWire(0, 0, 2 * rWire, vWire, "s");
  // Add the tube.
  cmp.AddTube(rTube, vTube, 0, "t");
  // Request calculation of the weighting field. 
  cmp.AddReadout("s");

  // Make a sensor.
  Sensor sensor;
  sensor.AddComponent(&cmp);
  sensor.AddElectrode(&cmp, "s");
  // Set the signal time window.
  const double tstep = 0.5;
  const double tmin = -0.5 * tstep;
  const unsigned int nbins = 5000;
  sensor.SetTimeWindow(tmin, tstep, nbins);
  // Set the delta reponse function.
    
    sensor.ClearSignal();
    //sensor.SetTransferFunction(transfer);
    

  // Set up Heed.
  TrackHeed track;
  track.SetParticle("muon");
  track.SetEnergy(3.e9);
  track.SetSensor(&sensor);


  // RKF integration.
  DriftLineRKF drift;
  drift.SetSensor(&sensor);
  drift.SetGainFluctuationsPolya(0., 20000.);
    //drift.EnableIonTail();
 

  int loop;
  string filename;
  string filename1="distance_";
  string filename2=".txt";

  for(loop=0;loop<14;loop++)
{        cout<<loop<<endl;                
	    filename=filename1+to_string(loop)+filename2;                      
            ofstream out(filename);
  //////////////////////////////////////////////////////////////////////////
  const double rTrack = loop*0.1;   
  const double x0 = rTrack;
  const double y0 = -sqrt(rTube * rTube - rTrack * rTrack);
  const unsigned int nTracks = 1000;
  ////////////////////////////////////////////////////////////////////////////

  for (unsigned int j = 0; j < nTracks; ++j) {
    sensor.ClearSignal();
    track.NewTrack(x0, y0, 0, 0, 0, 1, 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)) {
  
      for (int k = 0; k < nc; ++k) {
        double xe = 0., ye = 0., ze = 0., te = 0., ee = 0.;
        double dx = 0., dy = 0., dz = 0.;
        
        track.GetElectron(k, xe, ye, ze, te, ee, dx, dy, dz);
        
        drift.DriftElectron(xe, ye, ze, te);  
        }
       }

    sensor.ConvoluteSignals();
    int nt = 0;
    if (!sensor.ComputeThresholdCrossings(-1, "s", nt)) continue;
 	int i;
 	double time,level; bool rise;
 	sensor.GetThresholdCrossing(i,time,level,rise);
 	 std::cout<<"r="<<rTrack<<"     "<<"No."<<j+1<<"     "<<time<<"ns"<<std::endl;

   	out<<time<<std::endl;
  	}  
  	out.close();
 }
app.Run(kTRUE);
}