//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
/// \file GarfieldPhysics.cc
/// \brief Implementation of the GarfieldPhysics class
#include "GarfieldPhysics.hh"

#include "Garfield/AvalancheMC.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "GarfieldAnalysis.hh"
#include "Garfield/ComponentConstant.hh"
#include "Garfield/DriftLineRKF.hh"

#include "G4PhysicalConstants.hh"
#include "G4RandomDirection.hh"
#include "G4Threading.hh"
#include "G4ThreeVector.hh"


GarfieldPhysics* GarfieldPhysics::fGarfieldPhysics = nullptr;

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

GarfieldPhysics* GarfieldPhysics::GetInstance() {
  if (!fGarfieldPhysics) fGarfieldPhysics = new GarfieldPhysics();
  return fGarfieldPhysics;
}

void GarfieldPhysics::Dispose() {
  delete fGarfieldPhysics;
  fGarfieldPhysics = nullptr;
}

GarfieldPhysics::~GarfieldPhysics() {
  delete fMediumMagboltzDrift;
  delete fMediumMagboltzItf;
  delete fMediumMagboltzEL;
  // delete fSensor;
  delete fSensorDrift;
  delete fSensorItf;
  delete fSensorEL;
  // delete fTrackHeed;
  delete fTrackHeedDrift;
  delete fTrackHeedItf;
  delete fTrackHeedEL;

  std::cout << "Deconstructor GarfieldPhysics" << std::endl;
}

std::string GarfieldPhysics::GetIonizationModel() { return fIonizationModel; }

double GarfieldPhysics::GetELPhotonYieldPerMm(double E_kVcm, double pressure_bar) {
	const double Eth = 1.0;		// cut for kV/(cm*bar)
	const double S   = 140.0;	// constant for photons/(e*cm*(kV/cm)*bar)
	if (pressure_bar <= 0.0) return 0.0;
	const double Ep = E_kVcm / pressure_bar;
	const double over = std::max(0.0, Ep - Eth);
	const double per_cm = S * over * pressure_bar;
	return per_cm / 10.0;		// per cm to per mm
}

void GarfieldPhysics::SetIonizationModel(std::string model, bool useDefaults) {
  if (model != "PAIPhot" && model != "PAI" && model != "Heed") {
    std::cout << "Unknown ionization model " << model << std::endl;
    std::cout << "Using PAIPhot as default model!" << std::endl;
    model = "PAIPhot";
  }
  fIonizationModel = model;

  if (fIonizationModel == "PAIPhot" || fIonizationModel == "PAI") {
    if (useDefaults) {
      // Particle types and energies for which the G4FastSimulationModel with
      // Garfield++ is valid
      this->AddParticleName("e-", 0, 1e-3, "garfield");
      this->AddParticleName("gamma", 1e-6, 1e+8, "garfield");

      // Particle types and energies for which the PAI or PAIPhot model is valid
      this->AddParticleName("e-", 0, 1e+8, "geant4");
      this->AddParticleName("e+", 0, 1e+8, "geant4");
      this->AddParticleName("mu-", 0, 1e+8, "geant4");
      this->AddParticleName("mu+", 0, 1e+8, "geant4");
      this->AddParticleName("proton", 0, 1e+8, "geant4");
      this->AddParticleName("pi+", 0, 1e+8, "geant4");
      this->AddParticleName("pi-", 0, 1e+8, "geant4");
      this->AddParticleName("alpha", 0, 1e+8, "geant4");
      this->AddParticleName("He3", 0, 1e+8, "geant4");
      this->AddParticleName("GenericIon", 0, 1e+8, "geant4");
    }

  } else if (fIonizationModel == "Heed") {
    if (useDefaults) {
      // Particle types and energies for which the G4FastSimulationModel with
      // Garfield++ is valid
      this->AddParticleName("gamma", 1e-6, 1e+8, "garfield");
      this->AddParticleName("e-", 6e-2, 1e+7, "garfield");
      this->AddParticleName("e+", 6e-2, 1e+7, "garfield");
      this->AddParticleName("mu-", 1e+1, 1e+8, "garfield");
      this->AddParticleName("mu+", 1e+1, 1e+8, "garfield");
      this->AddParticleName("pi-", 2e+1, 1e+8, "garfield");
      this->AddParticleName("pi+", 2e+1, 1e+8, "garfield");
      this->AddParticleName("kaon-", 1e+1, 1e+8, "garfield");
      this->AddParticleName("kaon+", 1e+1, 1e+8, "garfield");
      this->AddParticleName("proton", 9.e+1, 1e+8, "garfield");
      this->AddParticleName("anti_proton", 9.e+1, 1e+8, "garfield");
      this->AddParticleName("deuteron", 2.e+2, 1e+8, "garfield");
      // this->AddParticleName("alpha", 4.e+2, 1e+8, "garfield");
      this->AddParticleName("alpha", 0.001, 1.e+8, "garfield"); // to ionise using alpha
    }
  }
}

void GarfieldPhysics::AddParticleName(const std::string particleName,
                                      double ekin_min_MeV, double ekin_max_MeV,
                                      std::string program) {
  if (ekin_min_MeV >= ekin_max_MeV) {
    std::cout << "Ekin_min=" << ekin_min_MeV
              << " keV is larger than Ekin_max=" << ekin_max_MeV << " keV"
              << std::endl;
    return;
  }

  if (program == "garfield") {
    std::cout << "Garfield model (Heed) is applicable for G4Particle "
              << particleName << " between " << ekin_min_MeV << " MeV and "
              << ekin_max_MeV << " MeV" << std::endl;

    fMapParticlesEnergyGarfield.insert(std::make_pair(
        particleName, std::make_pair(ekin_min_MeV, ekin_max_MeV)));
  } else {
    std::cout << fIonizationModel << " is applicable for G4Particle "
              << particleName << " between " << ekin_min_MeV << " MeV and "
              << ekin_max_MeV << " MeV" << std::endl;
    fMapParticlesEnergyGeant4.insert(std::make_pair(
        particleName, std::make_pair(ekin_min_MeV, ekin_max_MeV)));
  }
}

bool GarfieldPhysics::FindParticleName(std::string name, std::string program) {
  if (program == "garfield") {
    auto it = fMapParticlesEnergyGarfield.find(name);
    if (it != fMapParticlesEnergyGarfield.end()) return true;
  } else {
    auto it = fMapParticlesEnergyGeant4.find(name);
    if (it != fMapParticlesEnergyGeant4.end()) return true;
  }
  return false;
}

bool GarfieldPhysics::FindParticleNameEnergy(std::string name, double ekin_MeV,
                                             std::string program) {
  if (program == "garfield") {
    auto it = fMapParticlesEnergyGarfield.find(name);
    if (it != fMapParticlesEnergyGarfield.end()) {
      EnergyRange_MeV range = it->second;
      if (range.first <= ekin_MeV && range.second >= ekin_MeV) {
        return true;
      }
    }
  } else {
    auto it = fMapParticlesEnergyGeant4.find(name);
    if (it != fMapParticlesEnergyGeant4.end()) {
      EnergyRange_MeV range = it->second;
      if (range.first <= ekin_MeV && range.second >= ekin_MeV) {
        return true;
      }
    }
  }
  return false;
}

double GarfieldPhysics::GetMinEnergyMeVParticle(std::string name,
                                                std::string program) {
  if (program == "garfield") {
    auto it = fMapParticlesEnergyGarfield.find(name);
    if (it != fMapParticlesEnergyGarfield.end()) {
      EnergyRange_MeV range = it->second;
      return range.first;
    }
  } else {
    auto it = fMapParticlesEnergyGeant4.find(name);
    if (it != fMapParticlesEnergyGeant4.end()) {
      EnergyRange_MeV range = it->second;
      return range.first;
    }
  }
  return -1;
}

double GarfieldPhysics::GetMaxEnergyMeVParticle(std::string name,
                                                std::string program) {
  if (program == "garfield") {
    auto it = fMapParticlesEnergyGarfield.find(name);
    if (it != fMapParticlesEnergyGarfield.end()) {
      EnergyRange_MeV range = it->second;
      return range.second;
    }
  } else {
    auto it = fMapParticlesEnergyGeant4.find(name);
    if (it != fMapParticlesEnergyGeant4.end()) {
      EnergyRange_MeV range = it->second;
      return range.second;
    }
  }
  return -1;
}

void GarfieldPhysics::InitializePhysics() {
  // Define the gas mixture.

  // Set the Penning transfer efficiency.
  const double rPenning = 0.00;
  const double lambdaPenning = 0.00;

  G4cout << "System starts to read GarfieldPhysics::InitializePhysics()" << G4endl;

  if (!fInitDone) {
  	  G4cout << "[INIT] this = " << this
		 << " tid = " << G4Threading::G4GetThreadId()
		 << " isMaster = " << G4Threading::IsMasterThread() << G4endl;

	  fMediumMagboltzDrift = new Garfield::MediumMagboltz();
	  fMediumMagboltzDrift->SetComposition("xe", 100.);
  	  fMediumMagboltzDrift->SetTemperature(273);		// Kelvin
  	  fMediumMagboltzDrift->SetPressure(760 * 1);	// Torr
  	  fMediumMagboltzDrift->Initialise(true);
  	  fMediumMagboltzDrift->EnablePenningTransfer(rPenning, lambdaPenning, "xe");
  	  bool ok = fMediumMagboltzDrift->LoadGasFile("xe_100_1bar_440.gas");
	  G4cout << "[CHK] LoadGasFile = " << ok << G4endl;
  	  fComponentConstantDrift = new Garfield::ComponentConstant();
  	  fComponentConstantDrift->SetMedium(fMediumMagboltzDrift);
  	  fComponentConstantDrift->SetElectricField(-440, 0, 0);
  	  fSensorDrift = new Garfield::Sensor();
  	  fSensorDrift->Clear();
	  fSensorDrift->AddComponent(fComponentConstantDrift);
  	  fSensorDrift->SetArea(0, 16, -6, 6, -6, 6);
	  fSensorDrift->SetTimeWindow(0., 2e5, 40000);
	  int n = fSensorDrift->GetNumberOfComponents();
	  std::cout << "Number of Components in Drift Region= " << n << G4endl;
  	  fTrackHeedDrift = new Garfield::TrackHeed();
  	  fTrackHeedDrift->EnableDeltaElectronTransport();
  	  fTrackHeedDrift->SetSensor(fSensorDrift);
  	  G4cout << "[Init]GarfieldPhysics::InitializePhysics()\n"
  	         << "fSensorDrift = " << fSensorDrift << G4endl;
          G4cout << "[Init]GarfieldPhyscis::InitializePhysics()\n"
		  <<"fComponentConstantDrift = " << fComponentConstantDrift << G4endl;

  	  fMediumMagboltzItf = new Garfield::MediumMagboltz("xe", 100.);
  	  fMediumMagboltzItf->SetTemperature(273);		// Kelvin
  	  fMediumMagboltzItf->SetPressure(760 * 1);	// Torr
  	  fMediumMagboltzItf->Initialise(true);
	  fMediumMagboltzDrift->SetFieldGrid(-471e-7, -471e-7, 20,
			                     0, 0, 1,
			                     0, 0, 1);
  	  fMediumMagboltzItf->EnablePenningTransfer(rPenning, lambdaPenning, "xe");
  	  // fMediumMagboltzItf->LoadGasFile("xe_100_1bar_471.gas");
  	  fComponentConstantItf = new Garfield::ComponentConstant();
  	  fComponentConstantItf->SetMedium(fMediumMagboltzItf);
  	  fComponentConstantItf->SetElectricField(-471, 0, 0);
  	  fSensorItf = new Garfield::Sensor(fComponentConstantItf);
  	  fSensorItf->SetArea(16, 18.1,
  	        	      -6, 6   ,
  	        	      -6, 6   );
  	  fSensorItf->AddComponent(fComponentConstantItf);
  	  fTrackHeedItf = new Garfield::TrackHeed(fSensorItf);
  	  fTrackHeedItf->EnableDeltaElectronTransport();

  	  fMediumMagboltzEL = new Garfield::MediumMagboltz("xe", 100.);
  	  fMediumMagboltzEL->SetTemperature(273);		// Kelvin
  	  fMediumMagboltzEL->SetPressure(760 * 1);	// Torr
  	  fMediumMagboltzEL->Initialise(true);
	  fMediumMagboltzEL->SetFieldGrid(-140e-5, -140e-5, 20,
			                     0, 0, 1,
			                     0, 0, 1);
  	  fMediumMagboltzEL->EnablePenningTransfer(rPenning, lambdaPenning, "xe");
  	  // fMediumMagboltzEL->LoadGasFile("xe_100_1bar_14000.gas");
  	  fComponentConstantEL = new Garfield::ComponentConstant();
  	  fComponentConstantEL->SetMedium(fMediumMagboltzEL);
  	  fComponentConstantEL->SetElectricField(-14000, 0, 0);
  	  fSensorEL = new Garfield::Sensor(fComponentConstantEL);
  	  fSensorEL->SetArea(18.1, 18.81,
  	        	     -6  , 6    ,
  	        	     -6  , 6    );
  	  fSensorEL->AddComponent(fComponentConstantEL);
  	  fTrackHeedEL = new Garfield::TrackHeed(fSensorEL);
  	  fTrackHeedEL->EnableDeltaElectronTransport();

	  fInitDone = true;
  }
  constexpr double radius = 3.60000;		// [cm]
  constexpr double pPlate1 = 0.0000;		// [cm]
  constexpr double pPlate2 = 16.000;		// [cm]
  constexpr double pPlate3 = 18.100;		// [cm]
  constexpr double pPlate4 = 18.810;		// [cm]
  // constexpr double vPlate1 = 0; 
  // constexpr double vPlate2 = vPlate1 + 7.04e+3;
  // constexpr double vPlate3 = vPlate2 + 9.90e+3;
  // constexpr double vPlate4 = vPlate3 + 9.98e+3;

  Garfield::Medium* medium = nullptr;
  int status;
  double ex, ey, ez, v;
  std::cout << "====== E field Output (Component) ======" << std::endl;
  fComponentConstantDrift->ElectricField(8, 0, 0, ex, ey, ez, v, medium, status);
  std::cout << "Drift Zone Centre (8, 0, 0)" << std::endl;
  std::cout << "Ex = " << ex << " V/cm" << std::endl;
  std::cout << "Ey = " << ey << " V/cm" << std::endl;
  std::cout << "Ez = " << ez << " V/cm" << std::endl;
  fComponentConstantItf->ElectricField(17.05, 0, 0, ex, ey, ez, v, medium, status);
  std::cout << "Interface Zone Centre (17.05, 0, 0)" << std::endl;
  std::cout << "Ex = " << ex << " V/cm" << std::endl;
  std::cout << "Ey = " << ey << " V/cm" << std::endl;
  std::cout << "Ez = " << ez << " V/cm" << std::endl;
  fComponentConstantEL->ElectricField(18.455, 0, 0, ex, ey, ez, v, medium, status);
  std::cout << "EL Zone Centre (18.455, 0, 0)" << std::endl;
  std::cout << "Ex = " << ex << " V/cm" << std::endl;
  std::cout << "Ey = " << ey << " V/cm" << std::endl;
  std::cout << "Ez = " << ez << " V/cm" << std::endl;

  std::cout << "===== Memory Address (Component) =====" << std::endl;
  std::cout << "[MT]fSensorDrift = " << fComponentConstantDrift << G4endl;
  std::cout << "[MT]fSensorItf = " << fComponentConstantItf << G4endl;
  std::cout << "[MT]fSensorEL = " << fComponentConstantEL << G4endl;
  
  std::cout << "====== E field Output (Sensor) ======" << std::endl;
  fSensorDrift->ElectricField(8, 0, 0, ex, ey, ez, v, medium, status);
  std::cout << "Drift Zone Centre (8, 0, 0)" << std::endl;
  std::cout << "Ex = " << ex << " V/cm" << std::endl;
  std::cout << "Ey = " << ey << " V/cm" << std::endl;
  std::cout << "Ez = " << ez << " V/cm" << std::endl;
  fSensorItf->ElectricField(17.05, 0, 0, ex, ey, ez, v, medium, status);
  std::cout << "Interface Zone Centre (17.05, 0, 0)" << std::endl;
  std::cout << "Ex = " << ex << " V/cm" << std::endl;
  std::cout << "Ey = " << ey << " V/cm" << std::endl;
  std::cout << "Ez = " << ez << " V/cm" << std::endl;
  fSensorEL->ElectricField(18.455, 0, 0, ex, ey, ez, v, medium, status);
  std::cout << "EL Zone Centre (18.455, 0, 0)" << std::endl;
  std::cout << "Ex = " << ex << " V/cm" << std::endl;
  std::cout << "Ey = " << ey << " V/cm" << std::endl;
  std::cout << "Ez = " << ez << " V/cm" << std::endl;

  std::cout << "===== Memory Address (Sensor) =====" << std::endl;
  std::cout << "[MT]fSensorDrift = " << fSensorDrift << G4endl;
  std::cout << "[MT]fSensorItf = " << fSensorItf << G4endl;
  std::cout << "[MT]fSensorEL = " << fSensorEL << G4endl;
}

void GarfieldPhysics::DoIt(std::string particleName, double ekin_MeV,
                           double time, double x_cm, double y_cm, double z_cm,
                           double dx, double dy, double dz) {

  G4cout << "[GarfieldPhysics::DoIt] is activated\n"
	 << "Particle = " << particleName
	 << ", E = " << ekin_MeV / CLHEP::MeV << " MeV\n"
	 << "at (" << x_cm << ", " << y_cm << ", " << z_cm << ") cm\n" << G4endl;

  constexpr double radius  = 3.600;		// [cm]
  constexpr double pPlate1 = 0.000;		// [cm]
  constexpr double pPlate2 = 16.00;		// [cm]
  constexpr double pPlate3 = 18.10;		// [cm]
  constexpr double pPlate4 = 18.81;		// [cm]

  if (pPlate1 <= x_cm && x_cm < pPlate2) {
	  G4cout << "Particle is at drift zone" << G4endl;
          G4cout << "[DoIt] GarfieldPhysics::InitializePhysics()\n"
        	 << "fSensorDrift = " << fSensorDrift << G4endl;
	  // electron transfer ----------
          double ex, ey, ez, v;
  Garfield::Medium* med = nullptr;
  int status = 0;
  fComponentConstantDrift->ElectricField(x_cm, y_cm, z_cm, ex, ey, ez, v, med, status);
  G4cout << "[CHK] COMP  E(V/cm)=(" << ex << "," << ey << "," << ez
         << ") status=" << status << " med=" << med << G4endl;

  // 3) 센서 레벨에서도 필드를 가져와봄(컴포넌트 매핑 확인)
  fSensorDrift->ElectricField(x_cm, y_cm, z_cm, ex, ey, ez, v, med, status);
  G4cout << "[CHK] SENSOR E(V/cm)=(" << ex << "," << ey << "," << ez
         << ") status=" << status << " med=" << med << G4endl;
  G4cout << "fComponentConstantDrift = " << fComponentConstantDrift << G4endl;
	  Garfield::AvalancheMC drift;
	  drift.SetSensor(fSensorDrift);
	  drift.SetDistanceSteps(1e-1);		// mm
	  drift.DriftElectron(x_cm, y_cm, z_cm, time);
	  // electron transfer end -----
	  // check coordinate of electron in drift zone -----
	  double x0, y0, z0, t0, x1, y1, z1, t1;
	  // int status;
	  const int nEnd = drift.GetNumberOfElectronEndpoints();
	  G4cout << "number of segments: " << nEnd << G4endl;
	  drift.GetElectronEndpoint(0, x0, y0, z0, t0, x1, y1, z1, t1, status);
	  G4cout << "status: " << status << G4endl;
	  G4cout << "electron drift starts in drift zone at "
		 << "(" << x0 << ", " << y0 << ", " << z0 << ") cm" << G4endl;
	  drift.GetElectronEndpoint(nEnd - 1, x0, y0, z0, t0, x1, y1, z1, t1, status);
	  G4cout << "electron drift ends in drift zone at "
		 << "(" << x1 << ", " << y1 << ", " << z1 << ") cm" << G4endl;
	  // check coordinate of electron in drift zone end ----
  } else if (pPlate2 <= x_cm && x_cm < pPlate3) {
	  G4cout << "Particle is at interface zone" << G4endl;
  } else if (pPlate3 <= x_cm && x_cm < pPlate4) {
	  G4cout << "Particle is at El zone" << G4endl;
  }
}