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

#include <G4EmPenelopePhysics.hh>
#include <G4ProcessManager.hh>
#include <G4PhysicsListHelper.hh>
// gamma
#include <G4PhotoElectricEffect.hh>
#include <G4PenelopePhotoElectricModel.hh>
#include <G4ComptonScattering.hh>
#include <G4PenelopeComptonModel.hh>
#include <G4GammaConversion.hh>
#include <G4PenelopeGammaConversionModel.hh>
#include <G4RayleighScattering.hh>
#include <G4PenelopeRayleighModel.hh>
// e- and e+
#include <G4eMultipleScattering.hh>
#include <G4UniversalFluctuation.hh>
#include <G4eIonisation.hh>
#include <G4PenelopeIonisationModel.hh>
#include <G4eBremsstrahlung.hh>
#include <G4PenelopeBremsstrahlungModel.hh>
// e+ only
#include <G4eplusAnnihilation.hh>
#include <G4PenelopeAnnihilationModel.hh>
// mu
#include <G4MuMultipleScattering.hh>
#include <G4MuIonisation.hh>
#include <G4MuBremsstrahlung.hh>
#include <G4MuPairProduction.hh>
// hadrons
#include <G4hMultipleScattering.hh>
#include <G4MscStepLimitType.hh>
#include <G4hBremsstrahlung.hh>
#include <G4hPairProduction.hh>
#include <G4hIonisation.hh>
#include <G4ionIonisation.hh>
#include <G4alphaIonisation.hh>
#include <G4IonParametrisedLossModel.hh>
#include <G4NuclearStopping.hh>
// msc models
#include <G4WentzelVIModel.hh>
#include <G4UrbanMscModel.hh>
#include <G4CoulombScattering.hh>
// Em Options
#include <G4EmProcessOptions.hh>
// Deexcitation
#include <G4UAtomicDeexcitation.hh>

#include <RAT/EmPhysicsListV5.hh>
#include <RAT/PhysicsList.hh>
#include <RAT/Log.hh>

void RAT::EmPhysicsListV5::ConstructProcess()
{
  // SUMMARY:
  // GAMMA:
  //     - Photoelectric effect  | Penelope below 1 GeV
  //                             | Standard above
  //     - Compton scattering    | Penelope below 1 GeV
  //                             | Standard above
  //     - Gamma conversion      | Penelope below 1 GeV
  //                             | Standard above
  //     - Rayleigh scattering   | Penelope below 1 GeV
  //                             | Standard above
  // ELECTRON:
  //     - Multiple scattering   | Urban
  //     - Ionisation            | Penelope with universal fluctuation below 1 GeV
  //                             | Standard above
  //     - Bremsstrahlung        | Penelope below 1 GeV
  //                             | Standard above
  // POSITRON:
  //     - Multiple scattering   | Urban
  //     - Ionisation            | Penelope with universal fluctuation below 1 GeV
  //                             | Standard above
  //     - Bremsstrahlung        | Penelope below 1 GeV
  //                             | Standard above
  //     - Annihilation          | Penelope below 1 GeV
  //                             | Standard above
  // MUONS (+/-)
  //     - Mulitple scattering   | Wentzel VI
  //     - Single scattering     | Coulomb Scattering
  //     - Pair Production       | Standard
  //     - Bremsstrahlung        | Standard
  //     - Ionisation            | Standard with bespoke step function
  //
  // ALPHA
  //    - Multiple scattering    | Standard
  //    - Ionisation             | Standard with bespoke step function
  //    - Nuclear Stopping       | Standard
  //
  // PROTON
  //    - Multiple scattering    | Standard
  //    - Ionisation             | Standard with bespoke step function
  //    - Bremsstrahlung         | Standard
  //    - Pair Production        | Standard
  //    - Nuclear Stopping       | Standard
  //
  // PIONS, KAONS
  //    - Multiple scattering    | Standard
  //    - Ionisation             | Standard with Penelope step function
  //    - Bremsstrahlung         | Standard
  //    - Pair Production        | Standard
  //
  // DEUTERON, TRITON, GenericIon, He3
  //    - Multiple scattering    | Standard
  //    - Ionisation             | Standard
  //    - Nuclear Stopping       | Standard
  //
  const G4double penenlopeHighEnergyLimit = 1.0 * GeV; // The penelope models should cut out at 1GeV

  G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
  theParticleTable->GetIterator()->reset();
  while( ( *(theParticleTable->GetIterator()) )() )
    {
      G4ParticleDefinition* particle = theParticleTable->GetIterator()->value();
      G4String particleName = particle->GetParticleName();

      if( particleName == "gamma" )
        {
          G4PhotoElectricEffect* thePhotoElectricEffect = new G4PhotoElectricEffect();
          G4PenelopePhotoElectricModel* thePEPenelopeModel = new G4PenelopePhotoElectricModel();
          thePEPenelopeModel->SetHighEnergyLimit( penenlopeHighEnergyLimit );
          thePhotoElectricEffect->SetEmModel( thePEPenelopeModel, 1 );
          ph->RegisterProcess( thePhotoElectricEffect, particle );

          G4ComptonScattering* theComptonScattering = new G4ComptonScattering();
          G4PenelopeComptonModel* theComptonPenelopeModel = new G4PenelopeComptonModel();
          theComptonPenelopeModel->SetHighEnergyLimit( penenlopeHighEnergyLimit );
          theComptonScattering->SetEmModel( theComptonPenelopeModel, 1 );
          ph->RegisterProcess( theComptonScattering, particle );

          G4GammaConversion* theGammaConversion = new G4GammaConversion();
          G4PenelopeGammaConversionModel* theGammaPenelopeModel = new G4PenelopeGammaConversionModel();
          theGammaPenelopeModel->SetHighEnergyLimit( penenlopeHighEnergyLimit );
          theGammaConversion->SetEmModel( theGammaPenelopeModel, 1 );
          ph->RegisterProcess( theGammaConversion, particle );

          G4RayleighScattering* theRayleigh = new G4RayleighScattering();
          G4PenelopeRayleighModel* theRayleighPenelopeModel = new G4PenelopeRayleighModel();
          theRayleighPenelopeModel->SetHighEnergyLimit( penenlopeHighEnergyLimit );
          theRayleigh->SetEmModel( theRayleighPenelopeModel, 1 );
          ph->RegisterProcess( theRayleigh, particle );
        }
      else if( particleName == "e-" )
        {
          G4eMultipleScattering* msc = new G4eMultipleScattering();
          msc->AddEmModel( 0, new G4UrbanMscModel() );
          ph->RegisterProcess(msc, particle);

          G4eIonisation* eIoni = new G4eIonisation();
          G4PenelopeIonisationModel* theIoniPenelope = new G4PenelopeIonisationModel();
          theIoniPenelope->SetHighEnergyLimit( penenlopeHighEnergyLimit );
          eIoni->AddEmModel( 0, theIoniPenelope, new G4UniversalFluctuation() );
          ph->RegisterProcess( eIoni, particle );

          G4eBremsstrahlung* eBrem = new G4eBremsstrahlung();
          G4PenelopeBremsstrahlungModel* theBremPenelope = new G4PenelopeBremsstrahlungModel();
          theBremPenelope->SetHighEnergyLimit( penenlopeHighEnergyLimit );
          eBrem->AddEmModel( 0, theBremPenelope );
          ph->RegisterProcess( eBrem, particle );
        }
      else if( particleName == "e+" )
        {
          G4eMultipleScattering* msc = new G4eMultipleScattering();
          msc->AddEmModel( 0, new G4UrbanMscModel() );
          ph->RegisterProcess(msc, particle);

          G4eIonisation* eIoni = new G4eIonisation();
          G4PenelopeIonisationModel* theIoniPenelope = new G4PenelopeIonisationModel();
          theIoniPenelope->SetHighEnergyLimit( penenlopeHighEnergyLimit );
          eIoni->AddEmModel( 0, theIoniPenelope, new G4UniversalFluctuation() );
          ph->RegisterProcess( eIoni, particle );

          G4eBremsstrahlung* eBrem = new G4eBremsstrahlung();
          G4PenelopeBremsstrahlungModel* theBremPenelope = new G4PenelopeBremsstrahlungModel();
          theBremPenelope->SetHighEnergyLimit( penenlopeHighEnergyLimit );
          eBrem->AddEmModel( 0, theBremPenelope );
          ph->RegisterProcess( eBrem, particle );

          G4eplusAnnihilation* eAnni = new G4eplusAnnihilation();
          G4PenelopeAnnihilationModel* theAnnPenelope = new G4PenelopeAnnihilationModel();
          theAnnPenelope->SetHighEnergyLimit( penenlopeHighEnergyLimit );
          eAnni->AddEmModel( 0, theAnnPenelope );
          ph->RegisterProcess( eAnni, particle );
        }
      else if( !PhysicsList::OmitMuonic() && ( particleName == "mu-" || particleName == "mu+" ) )
        {
          G4MuMultipleScattering* msc = new G4MuMultipleScattering();
          msc->AddEmModel(0, new G4WentzelVIModel());
          ph->RegisterProcess( msc, particle );

          G4MuIonisation* muIoni = new G4MuIonisation();
          muIoni->SetStepFunction( 0.01, 50*um );
          ph->RegisterProcess( muIoni, particle );

          ph->RegisterProcess( new G4MuBremsstrahlung(), particle );
          ph->RegisterProcess( new G4MuPairProduction(), particle );
          ph->RegisterProcess( new G4CoulombScattering(), particle );
        }
      else if( !PhysicsList::OmitHadronic() && particleName == "alpha" )
        {
          ph->RegisterProcess( new G4hMultipleScattering(), particle );
          G4ionIonisation* ionIoni = new G4ionIonisation();
          ionIoni->SetStepFunction( 0.1, 0.01*um );
          ph->RegisterProcess(ionIoni, particle);
          ph->RegisterProcess( new G4NuclearStopping(), particle );
        }
      else if( !PhysicsList::OmitHadronic() && ( particleName == "proton" || particleName == "anti_proton" ) )
        {
          ph->RegisterProcess( new G4hMultipleScattering, particle );
          G4hIonisation* ioni = new G4hIonisation();
          ioni->SetStepFunction( 0.1, 0.01*um );
          ph->RegisterProcess( ioni, particle );
          ph->RegisterProcess( new G4hBremsstrahlung, particle );
          ph->RegisterProcess( new G4hPairProduction, particle );
          ph->RegisterProcess( new G4NuclearStopping, particle );
        }
      else if( !PhysicsList::OmitHadronic() && ( particleName == "pi+" || particleName == "pi-" || particleName == "kaon+" || particleName == "kaon-" ) )
        {
          ph->RegisterProcess( new G4hMultipleScattering, particle );
          G4ionIonisation* ionIoni = new G4ionIonisation();
          ionIoni->SetStepFunction( 0.2, 50*um );
          ph->RegisterProcess(ionIoni, particle);
          ph->RegisterProcess( new G4hBremsstrahlung, particle );
          ph->RegisterProcess( new G4hPairProduction, particle );
        }
      else if( !PhysicsList::OmitHadronic() && ( particleName == "deuteron" || particleName == "triton" ||  particleName == "GenericIon" || particleName == "He3" ) )
        {
          ph->RegisterProcess( new G4hMultipleScattering(), particle );
          ph->RegisterProcess( new G4ionIonisation(), particle );
          ph->RegisterProcess( new G4NuclearStopping(), particle );
        }
    }
  // Em options
  //
  G4EmProcessOptions opt;
  opt.SetVerbose( fVerbose );
  opt.SetMinEnergy(100*eV);
  opt.SetMaxEnergy(10*TeV);
  opt.SetDEDXBinning(220);
  opt.SetLambdaBinning(220);

  //opt.SetSplineFlag(true);
  opt.SetPolarAngleLimit(CLHEP::pi);

  // Deexcitation
  G4VAtomDeexcitation* deexcitation = new G4UAtomicDeexcitation();
  G4LossTableManager::Instance()->SetAtomDeexcitation(deexcitation);
  deexcitation->SetFluo(true);
}