#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 <G4GoudsmitSaundersonMscModel.hh>
#include <G4CoulombScattering.hh>
// Em Options
#include <G4EmProcessOptions.hh>

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

void RAT::EmPhysicsListV4::ConstructProcess()
{
  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();
          thePhotoElectricEffect->SetEmModel( new G4PenelopePhotoElectricModel() );
          ph->RegisterProcess( thePhotoElectricEffect, particle );

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

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

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

          G4eIonisation* eIoni = new G4eIonisation();
          eIoni->SetEmModel( new G4PenelopeIonisationModel() );
          eIoni->SetFluctModel( new G4UniversalFluctuation() );
          ph->RegisterProcess( eIoni, particle );

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

          G4eIonisation* eIoni = new G4eIonisation();
          eIoni->SetEmModel( new G4PenelopeIonisationModel() );
          eIoni->SetFluctModel( new G4UniversalFluctuation() );
          ph->RegisterProcess( eIoni, particle );

          G4eBremsstrahlung* eBrem = new G4eBremsstrahlung();
          eBrem->SetEmModel( new G4PenelopeBremsstrahlungModel() );
          ph->RegisterProcess( eBrem, particle );

          G4eplusAnnihilation* eAnni = new G4eplusAnnihilation();
          eAnni->SetEmModel( new G4PenelopeAnnihilationModel() ); // GEANT4 Bug?
          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 == "deuteron" || particleName == "triton" ||  particleName == "GenericIon" || particleName == "He3" ) )
        {
          ph->RegisterProcess( new G4hMultipleScattering(), particle );
          ph->RegisterProcess( new G4ionIonisation(), particle );
          ph->RegisterProcess( new G4NuclearStopping(), 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" )
        {
          ph->RegisterProcess( new G4hMultipleScattering, particle );
          ph->RegisterProcess( new G4hIonisation(), particle );
          ph->RegisterProcess( new G4hBremsstrahlung, particle );
          ph->RegisterProcess( new G4hPairProduction, 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(0.2);
}