#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/AvalancheMC.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("ArIso5new.gas");
    gas.LoadGasFile("ArIso5Oxygen0p1.gas");
//    gas.LoadGasFile("ArIso5_1pcO2.gas");

    const std::string path = std::getenv("GARFIELD_INSTALL");
  	gas.LoadIonMobility(path + "/share/Garfield/Data/IonMobility_Ar+_Ar.txt");
  	gas.LoadIonMobility(path + "/share/Garfield/Data/IonMobility_C8Hn+_iC4H10.txt");

    gas.EnableDrift(true);
    gas.EnableThermalMotion(true);
    gas.EnablePrimaryIonisation(true);

  // Set the Penning transfer efficiency.
  constexpr double rPenning = 0.26; // "3D simulation of micromegas detector performance" GUO Jun-Jun
  constexpr double lambdaPenning = 0.;
    gas.EnablePenningTransfer(rPenning, lambdaPenning, "ar");
    gas.SetMaxElectronEnergy(200.);
    gas.Initialise(true);
    // 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);

    AvalancheMC driftDe(&sensorD);
    AvalancheMC driftAe(&sensorA);

    driftDe.DisableAvalancheSizeLimit(); 
    driftAe.DisableAvalancheSizeLimit(); 

    driftDe.SetDistanceSteps(2.e-4);
    driftAe.SetDistanceSteps(2.e-4);
  
    // Use Heed for simulating the photon absorption.
    TrackHeed track(&sensorD);   
    track.EnableElectricField();


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

// drift: field, eta
std::cout << "\nField intensity in the middle of drift region:\n";
double y1 = 0.5 * (yMesh + yDrift); 
auto edrift = cmpD.ElectricField(0., y1, 0.);
double magnitude1 = sqrt(edrift[0] * edrift[0] + edrift[1] * edrift[1]);
std::cout << "At y = " << y1 << " cm:\n";
std::cout << "  E = " << magnitude1/1000 << " kV/cm\n";
double eta1 = 0.;
gas.ElectronAttachment(edrift[0], edrift[1], edrift[2], 0., 0., 0., eta1);
std::cout <<" eta1 = " << eta1 << "\n";
// avalanche: field, eta
std::cout << "\nField intensity in the middle of avalanche region:\n";
double y2 = 0.5 * (yMesh); 
auto edrift2 = cmpA.ElectricField(0., y2, 0.);
double magnitude2 = sqrt(edrift2[0] * edrift2[0] + edrift2[1] * edrift2[1]);
std::cout << "At y2 = " << y2 << " cm:\n";
std::cout << "  E = " << magnitude2/1000 << " kV/cm\n";
double eta2 = 0.;
gas.ElectronAttachment(edrift2[0], edrift2[1], edrift2[2], 0., 0., 0., eta2);
std::cout <<" eta2 = " << eta2 << "\n";

    // Histogram
     const int nBins = 500;
    TH1::StatOverflows(true);
    TH1F hElectrons("hElectrons", ";electrons after interaction with g;",
                    200, -0.5, 300 - 0.5);

    TH1F hElectrons1("hElectrons1", ";#electrons at the end of drift ;",
                    200, -0.5, 10 - 0.5);

    TH1F hElectrons11("hElectrons11", ";electrons at the end of drift ;",
                    200, -0.5, 300 - 0.5);

    TH1F hElectrons2("hElectrons2", ";electrons at the end of pad;",
                    200, -0.5, 2000000 - 0.5); 



    // Simulate 1000 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

int printCount = 0;
        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);


//    Drift all produced electrons
        driftDe.DriftElectron(xe, ye, ze, te);
        double xe1, ye1, ze1, te1;
        double xe2, ye2, ze2, te2;
        int status;
        driftDe.GetElectronEndpoint(0, xe1, ye1, ze1, te1, xe2, ye2, ze2, te2,status);

// status = -7 -> att, -5 ->  particle not inside a drift medium

if (std::abs(ye2 - yMesh) < 1e-6) {
//    std::cout << printCount + 1 << "  ye1 = " << ye1 << ", ye2 = " << ye2 << ", status = " << status << "\n";
    ++printCount;
}


  if (status != -5) continue;  // Skip if not status -5




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

            // Loop over electrons that arrived at the mesh
            int n_mesh = driftDe.GetNumberOfElectronEndpoints();
            if(n_mesh > 0) hElectrons1.Fill(n_mesh);

            for (int k2 = 0; k2 < n_mesh; ++k2) {
                double x1, y1, z1, t1, e1, x2, y2, z2, e2, t2;
                int status;

                 driftDe.GetElectronEndpoint(0, x1, y1, z1, t1, x2, y2, z2, t2, status);               
                 driftAe.AvalancheElectron(x2, yMesh - 1e-14, z2, t2, false);

            //     int n_pad = driftA.GetNumberOfElectronEndpoints();
                int n_pad = driftAe.GetNumberOfElectronEndpoints();


                nsum += n_pad;

            }

        }



//  Fill number of electrons at the pad
  if (nsum > 0) hElectrons2.Fill(nsum);

// Fill a single number of total electrons at the pad
 if (printCount > 0) hElectrons11.Fill(printCount);

 
 
        if (nsum > 0 && ne > 0 && printCount > 0) std::cout << ne << "  " << printCount << "  " << nsum << "\n"; 


    }

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

    TCanvas c1("c1", "", 600, 600);
    c1.cd();
    hElectrons1.SetFillColor(kBlue + 2);
    hElectrons1.SetLineColor(kBlue + 2);
    hElectrons1.Draw();
    c1.Update();

    TCanvas c11("c11", "", 600, 600);
    c11.cd();
    hElectrons11.SetFillColor(kBlue + 2);
    hElectrons11.SetLineColor(kBlue + 2);
    hElectrons11.Draw();
    c11.Update();


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

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