#include "G4ios.hh"
#include <iomanip>

#include "LGPhysicsList.hh"
#include "LGPhysicsListMessenger.hh"

#include "G4ParticleTypes.hh"
#include "G4ParticleTable.hh"
#include "G4ProcessManager.hh"

#include "G4Cerenkov.hh"
#include "G4OpAbsorption.hh"
#include "G4OpRayleigh.hh"
#include "G4OpBoundaryProcess.hh"

LGPhysicsList::LGPhysicsList() : G4VUserPhysicsList(), 
Physics_List_Messenger(new LGPhysicsListMessenger(this))
{
  Cerenkov_Process	= 0;
  Absorption_Process	= 0;
  Rayleigh_Process	= 0;
  Boundary_Process	= 0;
  
  SetVerboseLevel(0);
}

LGPhysicsList::~LGPhysicsList()
{
  delete Physics_List_Messenger;
}

//----------------------------------------------------------------------
//-------- 1. Construction of Particles
//----------------------------------------------------------------------
void LGPhysicsList::ConstructParticle()
{
  ConstructBosons();
  ConstructLeptons();
  ConstructMesons();
  ConstructBaryons();
}

//			-------- 1.1. Bosons
void LGPhysicsList::ConstructBosons()
{
  G4Geantino::GeantinoDefinition();	// pseudo-particle: geantino
  G4ChargedGeantino::ChargedGeantinoDefinition();		// pseudo-particle: charged geantino
  G4Gamma::GammaDefinition();		// gamma
  G4OpticalPhoton::OpticalPhotonDefinition();		// optical photon
}

//			-------- 1.2. Leptons
void LGPhysicsList::ConstructLeptons()
{
  G4Electron::ElectronDefinition();		
  G4Positron::PositronDefinition();		
  G4NeutrinoE::NeutrinoEDefinition();		
  G4AntiNeutrinoE::AntiNeutrinoEDefinition();		
  G4MuonPlus::MuonPlusDefinition();
  G4MuonMinus::MuonMinusDefinition();
  G4NeutrinoMu::NeutrinoMuDefinition();
  G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();
}

//			-------- 1.3. Mesons
void LGPhysicsList::ConstructMesons()
{
  G4PionPlus::PionPlusDefinition();
  G4PionMinus::PionMinusDefinition();
  G4PionZero::PionZeroDefinition();
}

//			-------- 1.4. Baryons
void LGPhysicsList::ConstructBaryons()
{
  G4Proton::ProtonDefinition();
  G4AntiProton::AntiProtonDefinition();
  G4Neutron::NeutronDefinition();
  G4AntiNeutron::AntiNeutronDefinition();
}


//			----------------------------------------------------------------------
//			-------- 2. Construction and Registration of Physics Processes
//			----------------------------------------------------------------------
void LGPhysicsList::ConstructProcess()
{
  AddTransportation();
  ConstructGeneral();
  ConstructEM();
  ConstructOp();
}

//			----------------------------------------------------------------------
//			-------- 2.1. Decay process
//			----------------------------------------------------------------------
#include "G4Decay.hh"

void LGPhysicsList::ConstructGeneral()
{
  G4Decay* Decay_Process = new G4Decay();
  theParticleIterator->reset();
  while( (*theParticleIterator)() )
    {
      G4ParticleDefinition* particle = theParticleIterator->value();
      G4ProcessManager* processManager = particle->GetProcessManager();
      if (Decay_Process->IsApplicable(*particle)) 
	{
	  processManager->AddProcess(Decay_Process);
	  processManager->SetProcessOrdering(Decay_Process, idxPostStep);		// set ordering for PostStepDoIt and AtRestDoIt
	  processManager->SetProcessOrdering(Decay_Process, idxAtRest);
	}
    }
}

//			----------------------------------------------------------------------
//			-------- 2.2. E&M process
//			----------------------------------------------------------------------
#include "G4ComptonScattering.hh"
#include "G4GammaConversion.hh"
#include "G4PhotoElectricEffect.hh"

#include "G4eIonisation.hh"
#include "G4eBremsstrahlung.hh"
#include "G4eplusAnnihilation.hh"

#include "G4MuIonisation.hh"
#include "G4MuBremsstrahlung.hh"
#include "G4MuPairProduction.hh"

#include "G4hIonisation.hh"

void LGPhysicsList::ConstructEM()
{
  theParticleIterator->reset();
  while( (*theParticleIterator)() )
    {
      G4ParticleDefinition* particle = theParticleIterator->value();
      G4ProcessManager* processManager = particle->GetProcessManager();
      G4String particleName = particle->GetParticleName();
      
      if (particleName == "gamma")		// construction of processes for gammas
	{
	  processManager->AddDiscreteProcess(new G4GammaConversion());
	  processManager->AddDiscreteProcess(new G4ComptonScattering());
	  processManager->AddDiscreteProcess(new G4PhotoElectricEffect());
	} 
      
      else if (particleName == "e-")			// construction of processes for electrons
	{
	  //processManager->AddProcess(new G4MultipleScattering(), -1, 1, 1);
	  processManager->AddProcess(new G4eIonisation(), -1, 2, 2);
	  processManager->AddProcess(new G4eBremsstrahlung(), -1, 3, 3);
	} 
      
      else if (particleName == "e+")		// construction of processes for positrons
	{
	  //processManager->AddProcess(new G4MultipleScattering(), -1, 1, 1);
	  processManager->AddProcess(new G4eIonisation(), -1, 2, 2);
	  processManager->AddProcess(new G4eBremsstrahlung(), -1, 3, 3);
	  processManager->AddProcess(new G4eplusAnnihilation(),  0,-1, 4);
	} 
      
      else if (particleName == "mu+" || particleName == "mu-")		// construction of processes for muons
	{
	  //processManager->AddProcess(new G4MultipleScattering(), -1, 1, 1);
	  processManager->AddProcess(new G4MuIonisation(), -1, 2, 2);
	  processManager->AddProcess(new G4MuBremsstrahlung(), -1, 3, 3);
	  processManager->AddProcess(new G4MuPairProduction(), -1, 4, 4);
	} 
      
      else 
	{
	  if ((!particle->IsShortLived()) && (particle->GetPDGCharge() != 0.0) && (particle->GetParticleName() != "chargedgeantino"))
	    // all others charged particles except geantino
	    {
	      //processManager->AddProcess(new G4MultipleScattering(), -1, 1, 1);
	      processManager->AddProcess(new G4hIonisation(), -1, 2, 2);
	    }
	}
    }
}

//			----------------------------------------------------------------------
//			-------- 2.3. Optical process
//			----------------------------------------------------------------------
void LGPhysicsList::ConstructOp()
{
  Cerenkov_Process	= new G4Cerenkov("Cerenkov");
  Absorption_Process	= new G4OpAbsorption();
  Rayleigh_Process	= new G4OpRayleigh();
  Boundary_Process	= new G4OpBoundaryProcess();
  
  SetVerbose(0);
  
  Cerenkov_Process->SetMaxNumPhotonsPerStep(300);	
// set the maximum number of Cherenkov photons allowed to be generated during a tracking step
  Cerenkov_Process->SetTrackSecondariesFirst(true);		
// If true, the primary particle tracking is interrupted and any produced Cherekov 
// photons are tracked next. When all have been tracked, the tracking of the primary resumes.
  
  G4OpticalSurfaceModel themodel = unified;
  Boundary_Process->SetModel(themodel);		// set the optical surface model
  
  theParticleIterator->reset();
  while( (*theParticleIterator)() )
    {
      G4ParticleDefinition* particle = theParticleIterator->value();
      G4ProcessManager* processManager = particle->GetProcessManager();
      G4String particleName = particle->GetParticleName();

      if (Cerenkov_Process->IsApplicable(*particle)) 
	{
	  // processManager->AddContinuousProcess(Cerenkov_Process);
	  processManager->AddProcess(Cerenkov_Process);          
	  processManager->SetProcessOrdering(Cerenkov_Process,idxPostStep);
	  
	}
      
      if (particleName == "opticalphoton") 
	{
	  G4cout << " AddDiscreteProcess to OpticalPhoton " << G4endl;
	  processManager->AddDiscreteProcess(Absorption_Process);
	  processManager->AddDiscreteProcess(Rayleigh_Process);
	  processManager->AddDiscreteProcess(Boundary_Process);
	}
    }
}


//			----------------------------------------------------------------------
//			-------- 3. Set the verbose level for the processes
//			----------------------------------------------------------------------
void LGPhysicsList::SetVerbose(G4int verbose)
{
  Cerenkov_Process		->SetVerboseLevel(verbose);	
  // set control flag for output message, Ex) If verbose=0, then Silent, if 1, Warning message, if 2, More
  Absorption_Process	->SetVerboseLevel(verbose);
  Rayleigh_Process		->SetVerboseLevel(verbose);
  Boundary_Process	->SetVerboseLevel(verbose);
}


//			----------------------------------------------------------------------
//			-------- 4. Set the number of photons(Cerenkov)
//			----------------------------------------------------------------------
void LGPhysicsList::SetNbOfPhotonsCerenkov(G4int MaxNumber)
{
  Cerenkov_Process		->SetMaxNumPhotonsPerStep(MaxNumber);
}


//			----------------------------------------------------------------------
//			-------- 4. Set the cuts
//			----------------------------------------------------------------------
void LGPhysicsList::SetCuts()
{
  SetCutsWithDefault();
  if (verboseLevel>0) DumpCutValuesTable();
}