#include <iostream>

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

#include <Garfield/ComponentConstant.hh>
#include <Garfield/FundamentalConstants.hh>
#include <Garfield/MediumMagboltz.hh>
#include <Garfield/Sensor.hh>
#include <Garfield/TrackHeed.hh>

using namespace Garfield;

int main(int argc, char **argv) {
    TApplication app("app", &argc, argv);

    TH1::StatOverflows(true);
    TH1F pElectrons("pElectrons", "Number of primary electrons", 100, 0, 100);
    TH1F dElectrons("dElectrons", "Number of delta electrons", 100, 0, 500);

    double temperature = 293.15; // Gas temperature [K]
    double pressure = 1.; // Gas pressure [atm]

    MediumMagboltz gas("Ar", 75., "iC4H10", 25.);
    gas.SetTemperature(temperature);
    gas.SetPressure(pressure * AtmosphericPressure);

    double lx = 5.; // Half-width along the x-axis [cm]
    double ly = 5.; // Half-width along the y-axis [cm]
    double lz = 1.; // Half-width along the z-axis [cm]
    double enorm = 100.; // Electric field norm [V/cm]

    ComponentConstant cmp;
    cmp.SetArea(-lx, -ly, 0., 
                 lx,  ly, lz);
    cmp.SetMedium(&gas);
    cmp.SetElectricField(enorm, 0., 0.);

    Sensor sensor(&cmp);

    // Transport of charged relativistic particles
    TrackHeed track(&sensor);
    track.SetParticle("muon");
    track.SetKineticEnergy(1e9);
    track.DisableDeltaElectronTransport();

    // Transport of delta electrons
    TrackHeed delta(&sensor);

    const unsigned int nEvents = 10000;
    for (unsigned int i = 0; i < nEvents; ++i) {
        if (i % 1000 == 0) { std::cout << i << "/" << nEvents << "\n"; }

        // Initial position and direction of the particle
        double x0 = 0., y0 = 0., z0 = 0., t0 = 0.;
        double dx0 = 0., dy0 = 0., dz0 = 1.;
        // Simulate the particle track
        track.NewTrack(x0, y0, z0, t0, dx0, dy0, dz0);

        int np = 0; // # of primary electrons
        int nd = 0; // # of delta electrons

        // Loop over the clusters
        for (const auto &cluster : track.GetClusters()) {
            np += cluster.electrons.size();
            // Loop over the cluster electrons
            for (const auto &electron : cluster.electrons) {
                int tmp = 0;
                delta.TransportDeltaElectron(electron.x, electron.y, electron.z, electron.t, electron.e,
                                             electron.dx, electron.dy, electron.dz, tmp);
                nd += tmp;
                /* Some I-O here */
                for (const auto &dcluster : delta.GetClusters()) {
                    for (const auto &delectron : dcluster.electrons) {
                        /* Some more I-O here */
                    }
                }
            }
        }
        pElectrons.Fill(np);
        dElectrons.Fill(nd);
    }
    std::cout << nEvents << "/" << nEvents << "\n";

    TCanvas c1("c1", "", 600, 600);
    c1.SetLeftMargin(0.16);
    pElectrons.GetXaxis()->SetTitle("number of primary electrons");
    pElectrons.Draw();

    TCanvas c2("c2", "", 600, 600);
    c2.SetLeftMargin(0.16);
    dElectrons.GetXaxis()->SetTitle("number of delta electrons");
    dElectrons.Draw();

    app.Run();
}