#include <iostream>
#include <fstream>
#include <cmath>

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

#include "Garfield/TrackHeed.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/SolidTube.hh"
#include "Garfield/GeometrySimple.hh"
#include "Garfield/ComponentConstant.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/FundamentalConstants.hh"
#include "Garfield/Random.hh"
#include "Garfield/Plotting.hh"
#include "Garfield/ComponentAnalyticField.hh"

#include "Garfield/ViewDrift.hh"
#include "Garfield/DriftLineRKF.hh"

#include "Garfield/AvalancheMicroscopic.hh"

using namespace Garfield;

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

    TApplication app("app", &argc, argv);
    plottingEngine.SetDefaultStyle();

    MediumMagboltz gas;
//    gas.LoadGasFile("Ar97p7_iso2p3_1atm.gas");
//    gas.LoadGasFile("ArIso5Oxygen0p1.gas");
//    gas.LoadGasFile("ArIso5new.gas");
    gas.LoadGasFile("ArIso5_1pcO2.gas");

const double rP = 0.20;
const double lambdaP = 0.;
gas.EnablePenningTransfer(rP, lambdaP, "ar");
    // Build the Micromegas geometry.
    // We use separate components for drift and amplification regions.
    // 50 um amplification gap.
    ComponentAnalyticField cmpA;
    cmpA.SetMedium(&gas);
    constexpr double yMesh = 0.005;
    constexpr double vMesh = 0.;
    constexpr double vPad = 320.;
    cmpA.AddPlaneY(yMesh, vMesh, "mesh1");
    cmpA.AddPlaneY(0., vPad, "pad");
    cmpA.AddReadout("mesh1");
    cmpA.AddReadout("pad");
    Sensor sensorA;
    sensorA.AddComponent(&cmpA);
    sensorA.SetArea(-5.0, 0., -5.0, 5.0, yMesh, 5.0);
    // 10 mm drift gap.
    ComponentAnalyticField cmpD;
    cmpD.SetMedium(&gas);
    constexpr double yDrift = yMesh + 1.0;
    constexpr double vDrift = -300.;
    cmpD.AddPlaneY(yDrift, vDrift, "drift");
    cmpD.AddPlaneY(yMesh, vMesh, "mesh2");
    cmpD.AddReadout("mesh2");
    cmpD.AddReadout("drift");
    Sensor sensorD;
    sensorD.AddComponent(&cmpD);
    sensorD.SetArea(-5.0, yMesh, -5.0, 5.0, yDrift, 5.0);
    AvalancheMicroscopic driftA, driftD;
    driftA.SetSensor(&sensorA);
    driftD.SetSensor(&sensorD);
    // Use Heed for simulating the photon absorption.
    TrackHeed track(&sensorD);   
    track.EnableElectricField();

gas.SetMaxElectronEnergy(200.); // not print messages!

    // Histogram
/*    const int nBins = 500;
    TH1::StatOverflows(true);
    TH1F hElectrons("hElectrons", ";Number of electrons;",
                    nBins, -0.5, nBins - 0.5);

    TH1F hElectrons2("hElectrons2", ";Number of electrons2;",
                    200, -0.5, 500000 - 0.5);
//                    200, -0.5, 80000 - 0.5);*/

    // Simulate 100 events
    const unsigned int nEvents = 1000;
    for (unsigned int i = 0; i < nEvents; ++i) {
//        if (i % 10 == 0) std::cout << i << "/" << nEvents << "\n";
        int nsum = 0;
        int ne = 0;

        // Initial coordinates of the photon.
        const double x0 = 0.;
        const double y0 = 1.0;
        const double z0 = 0.;
        const double t0 = 0.;
        // Sample the photon energy, using the relative intensities according to XDB.
        const double r = 167. * RndmUniform();
        const double egamma = r < 100. ? 5898.8 : r < 150. ? 5887.6 : 6490.4;

        track.TransportPhoton(x0, y0, z0, t0, egamma, 0., -1.0, 0., ne);

        // Fill number of electrons in ionization
        // if (ne > 0) hElectrons.Fill(ne);

        // Loop over electrons in ionization
        for (int k1 = 0; k1 < ne; ++k1) {
            double xe, ye, ze, te, ee, dx, dy, dz;
            track.GetElectron(k1, xe, ye, ze, te, ee, dx, dy, dz);

            // Avalanche electrons in drifting area
            // (should not create much of an avalanche)
            driftD.AvalancheElectron(xe, ye, ze, te, 0, 0, 0, 0);

            // Loop over electrons that arrived at the mesh
            int n_mesh = driftD.GetNumberOfElectronEndpoints();
            for (int k2 = 0; k2 < n_mesh; ++k2) {
                double x1, y1, z1, t1, e1, x2, y2, z2, e2, t2;
                int status;
                driftD.GetElectronEndpoint(k2, x1, y1, z1, t1, e1, x2, y2, z2, t2, e2, status);

                // Resimulate - actual avalanche
                driftA.AvalancheElectron(x2, yMesh - 1e-10, z2, t2, 0, 0, 0, 0);

                int n_pad = driftA.GetNumberOfElectronEndpoints();
                nsum += n_pad;
            }
        }
        // Fill number of electrons at the pad
        //if (nsum > 0) hElectrons2.Fill(nsum);
 
        // if (nsum > 0 && ne > 0) std::cout << "nsum = " << nsum << " ne = " << ne << "\n";
        if (nsum > 0 && ne > 0) std::cout << nsum << "  " << ne << "\n"; 
    }

/*    TCanvas c("c", "", 600, 600);
    c.cd();
    hElectrons.SetFillColor(kBlue + 2);
    hElectrons.SetLineColor(kBlue + 2);
    hElectrons.Draw();
    c.Update();

    TCanvas c2("c2", "", 600, 600);
    c2.cd();
    hElectrons2.SetFillColor(kBlue + 2);
    hElectrons2.SetLineColor(kBlue + 2);
    hElectrons2.Draw();

    c2.Update();*/
    
    
 //   app.Run(true);
    }