#include <CLHEP/Units/SystemOfUnits.h>
using namespace CLHEP;

#include <G4ProcessManager.hh>
#include <G4HadronElasticProcess.hh>
#include <G4ProtonInelasticProcess.hh>
#include <G4NeutronHPElastic.hh>
#include <G4NeutronHPElasticData.hh>
#include <G4NeutronHPThermalScattering.hh>
#include <G4NeutronHPThermalScatteringData.hh>
#include <G4NeutronHPInelastic.hh>
#include <G4NeutronHPInelasticData.hh>
#include <G4NeutronInelasticProcess.hh>
#include <G4HadronElastic.hh>
#include <G4NeutronRadCapture.hh>
#include <G4NeutronHPCaptureData.hh>
#include <G4HadronCaptureProcess.hh>
#include <G4HadronFissionProcess.hh>
#include <G4NeutronHPFission.hh>
#include <G4NeutronHPFissionData.hh>
#include <G4LFission.hh>
#include <G4DeuteronInelasticProcess.hh>
#include <G4TritonInelasticProcess.hh>
#include <G4AlphaInelasticProcess.hh>

#include <G4CascadeInterface.hh>
#include <G4BinaryLightIonReaction.hh>
#include <G4NeutronHPCapture.hh>
#include <G4HadronicAbsorptionBertini.hh>

#include <RAT/HadronicPhysicsListV4.hh>
#include <RAT/Log.hh>

void RAT::HadronicPhysicsListV4::ConstructProcess()
{
  // SUMMARY:
  // More details at http://geant4.cern.ch/support/proc_mod_catalog/processes/
  // PROTON:
  //         - Low Energy Elastic model            | Generic elastic low energy
  //         - Low Energy Inelastic model          | Proton energies 0 - 25 GeV, final states
  //
  // NEUTRON:
  //         - High precision Elastic model        | neutron energies 4eV - 20 MeV
  //         - Low Energy Elastic model            | above 19MeV only
  //         - Thermal model                       | below 4eV
  //         - High precision Inelastic model      | neutron energies 0 - 20 MeV
  //         - Low energy inelastic model          | above 19MeV only
  //         - High precision Capture model        | neutron capture for neutron energies 0 - 20 MeV
  //         - High precision Fission model        | neutron fission for neutron energies 0 - 20 MeV
  //         - Low energy fission model            | above 19MeV only
  //         - Neutron capture at Rest             | At rest neutrons specialty
  //
  // DEUTERON:
  //         - Low Energy Elastic model            | Generic elastic low energy
  //         - Low Energy Deuteron Inelastic model | deuteron energies 0 - 100 MeV
  //
  // TRITON:
  //         - Low Energy Elastic model            | Generic elastic low energy
  //         - Low Energy Triton Inelastic model   | triton energies 0 - 100 MeV
  //
  // ALPHA:
  //         - Low Energy Elastic model            | Generic elastic low energy
  //         - Low Energy Alpha Inelastic model    | alpha energies 0 - 100 MeV

  G4ProcessManager* pManager = NULL;
  ///////////////////////////////////////////////////////////////////////////////////
  // Proton physics
  ///////////////////////////////////////////////////////////////////////////////////
  {
    pManager = G4Proton::Proton()->GetProcessManager();

    G4HadronElasticProcess* protonElasticProcess = new G4HadronElasticProcess();
    G4HadronElastic* protonElasticModel = new G4HadronElastic();
    protonElasticProcess->RegisterMe( protonElasticModel );
    pManager->AddDiscreteProcess( protonElasticProcess );

    G4ProtonInelasticProcess* protonInelasticProcess = new G4ProtonInelasticProcess();
    G4CascadeInterface* protonInelasticModel = new G4CascadeInterface();
    protonInelasticProcess->RegisterMe( protonInelasticModel );
    pManager->AddDiscreteProcess( protonInelasticProcess );
  }
  ///////////////////////////////////////////////////////////////////////////////////
  // Neutron Physics
  ///////////////////////////////////////////////////////////////////////////////////
  {
    pManager = G4Neutron::Neutron()->GetProcessManager();

    G4HadronElasticProcess* neutronElasticProcess = new G4HadronElasticProcess();
    G4HadronElastic* neutronElasticModel = new G4HadronElastic();
    neutronElasticModel->SetMinEnergy( 19.0 * MeV ); // HP below 19MeV, le above
    G4NeutronHPElastic* neutronElasticModelHP = new G4NeutronHPElastic();
    neutronElasticModelHP->SetMinEnergy( 4.0 * eV ); // Thermal below 4eV
    G4NeutronHPElasticData * neutronElasticData = new G4NeutronHPElasticData();
    G4NeutronHPThermalScattering* neutronThermalElasticModel = new G4NeutronHPThermalScattering();
    G4NeutronHPThermalScatteringData* neutronThermalScatteringData = new G4NeutronHPThermalScatteringData();
    neutronElasticProcess->RegisterMe( neutronThermalElasticModel );
    neutronElasticProcess->RegisterMe( neutronElasticModelHP );
    neutronElasticProcess->RegisterMe( neutronElasticModel );
    neutronElasticProcess->AddDataSet( neutronThermalScatteringData );
    neutronElasticProcess->AddDataSet( neutronElasticData );
    pManager->AddDiscreteProcess( neutronElasticProcess );

    G4NeutronInelasticProcess* neutronInelasticProcess = new G4NeutronInelasticProcess();
    G4CascadeInterface* neutronInelasticModel = new G4CascadeInterface();
    neutronInelasticModel->SetMinEnergy( 19.0 * MeV); // HP below 19MeV, le above
    G4NeutronHPInelastic* neutronInelasticModelHP = new G4NeutronHPInelastic();
    G4NeutronHPInelasticData * neutronInelasticData = new G4NeutronHPInelasticData();
    neutronInelasticProcess->RegisterMe( neutronInelasticModelHP );
    neutronInelasticProcess->RegisterMe( neutronInelasticModel );
    neutronInelasticProcess->AddDataSet( neutronInelasticData );
    pManager->AddDiscreteProcess( neutronInelasticProcess );

    G4HadronCaptureProcess* neutronCaptureProcess = new G4HadronCaptureProcess();
    G4NeutronHPCapture* neutronCaptureModelHP = new G4NeutronHPCapture();
    G4NeutronHPCaptureData* neutronCaptureData = new G4NeutronHPCaptureData();
    G4NeutronRadCapture* neutronCaptureModel = new G4NeutronRadCapture();
    neutronCaptureModel->SetMinEnergy( 19.0 * MeV ); // HP below 19MeV, le above
    neutronCaptureProcess->RegisterMe( neutronCaptureModelHP );
    neutronCaptureProcess->AddDataSet( neutronCaptureData );
    neutronCaptureProcess->RegisterMe( neutronCaptureModel );
    pManager->AddDiscreteProcess( neutronCaptureProcess );

    G4HadronFissionProcess* neutronFissionProcess = new G4HadronFissionProcess();
    G4LFission* neutronFissionModel = new G4LFission();
    neutronFissionModel->SetMinEnergy( 19.0 * MeV); // HP below 19MeV, le above
    G4NeutronHPFission* neutronFissionModelHP = new G4NeutronHPFission();
    G4NeutronHPFissionData* neutronFissionData = new G4NeutronHPFissionData();
    neutronFissionProcess->RegisterMe( neutronFissionModelHP );
    neutronFissionProcess->AddDataSet( neutronFissionData );
    neutronFissionProcess->RegisterMe( neutronFissionModel );
    pManager->AddDiscreteProcess( neutronFissionProcess );

    G4HadronicAbsorptionBertini* neutronAtRestProcess = new G4HadronicAbsorptionBertini();
    pManager->AddRestProcess( neutronAtRestProcess );
  }
  ///////////////////////////////////////////////////////////////////////////////////
  // Deuteron Physics
  ///////////////////////////////////////////////////////////////////////////////////
  {
    pManager = G4Deuteron::Deuteron()->GetProcessManager();

    G4HadronElasticProcess* thedueElasticProcess = new G4HadronElasticProcess();
    G4HadronElastic* thedueElasticModel = new G4HadronElastic();
    thedueElasticProcess->RegisterMe(thedueElasticModel);
    pManager->AddDiscreteProcess(thedueElasticProcess);

    G4DeuteronInelasticProcess* theDeuteronInelasticProcess = new G4DeuteronInelasticProcess();
    G4BinaryLightIonReaction* theDeuteronLEPModel = new G4BinaryLightIonReaction();
    theDeuteronInelasticProcess->RegisterMe(theDeuteronLEPModel);
    pManager->AddDiscreteProcess(theDeuteronInelasticProcess);
  }
  ///////////////////////////////////////////////////////////////////////////////////
  // Triton Physics
  ///////////////////////////////////////////////////////////////////////////////////
  {
    pManager = G4Triton::Triton()->GetProcessManager();

    G4HadronElasticProcess* thetriElasticProcess = new G4HadronElasticProcess();
    G4HadronElastic* thetriElasticModel = new G4HadronElastic();
    thetriElasticProcess->RegisterMe(thetriElasticModel);
    pManager->AddDiscreteProcess(thetriElasticProcess);

    G4TritonInelasticProcess* theTritonInelasticProcess = new G4TritonInelasticProcess();
    G4BinaryLightIonReaction* theTritonLEPModel = new G4BinaryLightIonReaction();
    theTritonInelasticProcess->RegisterMe(theTritonLEPModel);
    pManager->AddDiscreteProcess(theTritonInelasticProcess);
  }
  ///////////////////////////////////////////////////////////////////////////////////
  // Alpha Physics
  ///////////////////////////////////////////////////////////////////////////////////
  {
    pManager = G4Alpha::Alpha()->GetProcessManager();

    G4HadronElasticProcess* thealElasticProcess = new G4HadronElasticProcess();
    G4HadronElastic* thealElasticModel = new G4HadronElastic();
    thealElasticProcess->RegisterMe(thealElasticModel);
    pManager->AddDiscreteProcess(thealElasticProcess);

    G4AlphaInelasticProcess* theAlphaInelasticProcess = new G4AlphaInelasticProcess();
    G4BinaryLightIonReaction* theAlphaLEPModel = new G4BinaryLightIonReaction();
    theAlphaInelasticProcess->RegisterMe(theAlphaLEPModel);
    pManager->AddDiscreteProcess(theAlphaInelasticProcess);
  }
}