#include <G4EmPenelopePhysics.hh>
#include <G4ProcessManager.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/EmPhysicsListV3.hh>
#include <RAT/PhysicsList.hh>
#include <CLHEP/Units/SystemOfUnits.h>
using namespace CLHEP;
#include <RAT/Log.hh>

void RAT::EmPhysicsListV3::ConstructProcess()
{
  G4ParticleDefinition* particle = 0;
  G4String particleName = "";

  theParticleTable->GetIterator()->reset();
  while( ( *(theParticleTable->GetIterator()) )() )
    {
      particle = theParticleTable->GetIterator()->value();
      G4ProcessManager* pmanager = particle->GetProcessManager();
      particleName = particle->GetParticleName();

	  const double kPenelopeHighEnergyLimit = 1.0*GeV;

      if (particleName == "gamma")
		{
		  //Photo-electric effect
		  G4PhotoElectricEffect* thePhotoElectricEffect = new G4PhotoElectricEffect();
		  G4PenelopePhotoElectricModel* thePEPenelopeModel = new G4PenelopePhotoElectricModel();
		  thePEPenelopeModel->SetHighEnergyLimit(kPenelopeHighEnergyLimit);
		  thePhotoElectricEffect->AddEmModel(0,thePEPenelopeModel);
		  pmanager->AddDiscreteProcess(thePhotoElectricEffect);

		  //Compton scattering
		  G4ComptonScattering* theComptonScattering = new G4ComptonScattering();
		  G4PenelopeComptonModel* theComptonPenelopeModel = new G4PenelopeComptonModel();
		  theComptonPenelopeModel->SetHighEnergyLimit(kPenelopeHighEnergyLimit);
		  theComptonScattering->AddEmModel(0,theComptonPenelopeModel);
		  pmanager->AddDiscreteProcess(theComptonScattering);

		  //Gamma conversion
		  G4GammaConversion* theGammaConversion = new G4GammaConversion();
		  G4PenelopeGammaConversionModel* theGCPenelopeModel = new G4PenelopeGammaConversionModel();
		  theGammaConversion->AddEmModel(0,theGCPenelopeModel);
		  pmanager->AddDiscreteProcess(theGammaConversion);

		  //Rayleigh scattering
		  G4RayleighScattering* theRayleigh = new G4RayleighScattering();
		  G4PenelopeRayleighModel* theRayleighPenelopeModel = new G4PenelopeRayleighModel();
		  theRayleighPenelopeModel->SetHighEnergyLimit(kPenelopeHighEnergyLimit);
		  theRayleigh->AddEmModel(0,theRayleighPenelopeModel);
		  pmanager->AddDiscreteProcess(theRayleigh);
		}
      else if (particleName == "e-")
		{
		  G4eMultipleScattering* msc = new G4eMultipleScattering();
		  msc->AddEmModel(0, new G4GoudsmitSaundersonMscModel());
		  msc->SetStepLimitType(fUseDistanceToBoundary);
		  pmanager->AddProcess(msc, -1, 1, 1);

		  //Ionisation
		  G4eIonisation* eIoni = new G4eIonisation();
		  G4PenelopeIonisationModel* theIoniPenelope = new G4PenelopeIonisationModel();
		  theIoniPenelope->SetHighEnergyLimit(kPenelopeHighEnergyLimit);
		  eIoni->AddEmModel(0,theIoniPenelope,new G4UniversalFluctuation());
		  eIoni->SetStepFunction(0.2, 100*um);
		  pmanager->AddProcess(eIoni, -1, 2, 2);

		  //Bremsstrahlung
		  G4eBremsstrahlung* eBrem = new G4eBremsstrahlung();
		  G4PenelopeBremsstrahlungModel* theBremPenelope = new G4PenelopeBremsstrahlungModel();
		  theBremPenelope->SetHighEnergyLimit(kPenelopeHighEnergyLimit);
		  eBrem->AddEmModel(0,theBremPenelope);
		  pmanager->AddProcess(eBrem, -1,-3, 3);
		}
      else if (particleName == "e+")
		{
		  G4eMultipleScattering* msc = new G4eMultipleScattering();
		  msc->AddEmModel(0, new G4GoudsmitSaundersonMscModel());
		  msc->SetStepLimitType(fUseDistanceToBoundary);
		  pmanager->AddProcess(msc, -1, 1, 1);

		  //Ionisation
		  G4eIonisation* eIoni = new G4eIonisation();
		  G4PenelopeIonisationModel* theIoniPenelope = new G4PenelopeIonisationModel();
		  theIoniPenelope->SetHighEnergyLimit(kPenelopeHighEnergyLimit);
		  eIoni->AddEmModel(0,theIoniPenelope,new G4UniversalFluctuation());
		  eIoni->SetStepFunction(0.2, 100*um);
		  pmanager->AddProcess(eIoni, -1, 2, 2);

		  //Bremsstrahlung
		  G4eBremsstrahlung* eBrem = new G4eBremsstrahlung();
		  G4PenelopeBremsstrahlungModel* theBremPenelope = new G4PenelopeBremsstrahlungModel();
		  theBremPenelope->SetHighEnergyLimit(kPenelopeHighEnergyLimit);
		  eBrem->AddEmModel(0,theBremPenelope);
		  pmanager->AddProcess(eBrem, -1,-3, 3);

		  //Annihilation
		  G4eplusAnnihilation* eAnni = new G4eplusAnnihilation();
		  G4PenelopeAnnihilationModel* theAnnPenelope = new G4PenelopeAnnihilationModel();
		  theAnnPenelope->SetHighEnergyLimit(kPenelopeHighEnergyLimit);
		  eAnni->AddEmModel(0,theAnnPenelope);
		  pmanager->AddProcess(eAnni,0,-1, 4);
		}
      else if (particleName == "mu-" || particleName == "mu+")
		{
		  if ( !PhysicsList::OmitMuonic() )
            info << "EmPhysicsListV3: muon processes omitted" << newline;
		  else
			{
			  G4MuMultipleScattering* msc = new G4MuMultipleScattering();
			  msc->AddEmModel(0, new G4WentzelVIModel());
			  pmanager->AddProcess(msc, -1, 1, 1);

			  G4MuIonisation* muIoni = new G4MuIonisation();
			  muIoni->SetStepFunction(0.2, 50*um);

			  pmanager->AddProcess(muIoni, -1, 2, 2);
			  pmanager->AddProcess(new G4MuBremsstrahlung, -1,-3, 3);
			  pmanager->AddProcess(new G4MuPairProduction, -1,-4, 4);
			  pmanager->AddDiscreteProcess(new G4CoulombScattering());
			}

		}
	  else if (particleName == "proton" || particleName == "antiproton" || particleName == "pi+" || particleName == "pi-")
		{
		  pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);

		  G4hIonisation* hIoni = new G4hIonisation();
		  hIoni->SetStepFunction(0.2, 50*um);

		  pmanager->AddProcess(hIoni,                     -1, 2, 2);
		  pmanager->AddProcess(new G4hBremsstrahlung,     -1,-3, 3);
		  pmanager->AddProcess(new G4hPairProduction,     -1,-4, 4);
		}
	  else if (particleName == "alpha" || particleName == "GenericIon" || particleName == "He3" || particleName == "deuteron" || particleName == "triton" )
		{
		  pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);

		  G4ionIonisation* ionIoni = new G4ionIonisation();
		  ionIoni->SetStepFunction(0.1, 20*um);
		  pmanager->AddProcess(ionIoni,                   -1, 2, 2);
		  pmanager->AddProcess(new G4NuclearStopping(),   -1, 3,-1);
		  pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
		}
    }

  // 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);
}