//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

#include "globals.hh"
#include "G4ios.hh"
#include "sr1PhysicsList.hh"
#include "G4VPhysicsConstructor.hh"
#include "G4ProcessManager.hh"
#include "G4ParticleTypes.hh"

#include "G4MuonDecayChannel.hh"
#include "G4DecayTable.hh"
//#include "sr1MuonDecayChannel.hh"
//#include "sr1Decay.hh"
//cks Added to have Geant default muon decay with spin
#include "G4MuonDecayChannelWithSpin.hh"
//#include "G4AtRestSpinPrecession.hh"

//TS Muonium "particle"
#include "sr1Muonium.hh"
#include "MuDecayChannel.hh"
#include "MuDecayChannelWithSpin.hh"

#include "sr1Parameters.hh"

// Constructor and Destructor  
sr1PhysicsList::sr1PhysicsList():G4VUserPhysicsList()
{
  defaultCutValue = 0.1*mm; // 1.0*cm; 
  SetVerboseLevel(0);
}

sr1PhysicsList::~sr1PhysicsList()
{}


// Construction of the different particles (calls different construction methods)

void sr1PhysicsList::ConstructParticle()
{
  ConstructBosons();
  ConstructLeptons();
  ConstructBaryons();
}

// Bosons

void sr1PhysicsList::ConstructBosons()
{
  // pseudo-particles
  G4Geantino::GeantinoDefinition();
  G4ChargedGeantino::ChargedGeantinoDefinition();

  // gamma
  G4Gamma::GammaDefinition();
}

// Leptons

void sr1PhysicsList::ConstructLeptons()
{
  //  e+/-
  G4Electron::ElectronDefinition();
  G4Positron::PositronDefinition();
  // mu+/-
  G4MuonPlus::MuonPlusDefinition();
  G4MuonMinus::MuonMinusDefinition();
  // Muonium - TS
  sr1Muonium::MuoniumDefinition();
  
  G4NeutrinoE::NeutrinoEDefinition();
  G4AntiNeutrinoE::AntiNeutrinoEDefinition();
  // nu_mu
  G4NeutrinoMu::NeutrinoMuDefinition();
  G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();

  // Decays of muon and muonium
  //cks:  Trial to use Geant4 muon decay with spin
  G4DecayTable* MuonPlusDecayTable = new G4DecayTable();
  MuonPlusDecayTable -> Insert(new G4MuonDecayChannelWithSpin("mu+",1.00));
  G4MuonPlus::MuonPlusDefinition() -> SetDecayTable(MuonPlusDecayTable);
  
  //TS: Using the muonium decay with and without spin
  G4DecayTable* MuoniumDecayTable = new G4DecayTable();
  MuoniumDecayTable -> Insert(new MuDecayChannel("Mu",0.50));
  MuoniumDecayTable -> Insert(new MuDecayChannelWithSpin("Mu",0.5));
  sr1Muonium::MuoniumDefinition() -> SetDecayTable(MuoniumDecayTable);
  //MuoniumDecayTable ->DumpInfo(); // Info on muonium decay channels

}

// Mesons (only some light mesons)

void sr1PhysicsList::ConstructMesons()
{
  G4PionPlus::PionPlusDefinition();
  G4PionMinus::PionMinusDefinition();
  G4PionZero::PionZeroDefinition();
}

// Baryons

void sr1PhysicsList::ConstructBaryons()
{
  G4Proton::ProtonDefinition();
  G4AntiProton::AntiProtonDefinition();
  G4Neutron::NeutronDefinition();
  G4AntiNeutron::AntiNeutronDefinition();
}




// Physical processes - general

void sr1PhysicsList::ConstructProcess()
{
  AddTransportation();
  ConstructEM();
  ConstructGeneral();
}

// Declaration of the different processes
#include "G4ComptonScattering.hh"
#include "G4GammaConversion.hh"
#include "G4PhotoElectricEffect.hh"

#include "G4MultipleScattering.hh"

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

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

#include "G4hIonisation.hh"

#include "G4UserSpecialCuts.hh"

//#include "sr1AtRestSpinRotation.hh"

// For low energy physics processes:
#include "G4LowEnergyCompton.hh"
//#include "G4LowEnergyPolarizedCompton.hh"
#include "G4LowEnergyGammaConversion.hh"
#include "G4LowEnergyPhotoElectric.hh"
#include "G4LowEnergyRayleigh.hh"
#include "G4LowEnergyBremsstrahlung.hh"
#include "G4LowEnergyIonisation.hh"
#include "G4hLowEnergyIonisation.hh"


// For Penelope processes:
#include "G4PenelopeCompton.hh"
#include "G4PenelopeGammaConversion.hh"
#include "G4PenelopePhotoElectric.hh"
#include "G4PenelopeRayleigh.hh"
#include "G4PenelopeIonisation.hh"
#include "G4PenelopeBremsstrahlung.hh"
#include "G4PenelopeAnnihilation.hh"

// For Coulomb scattering instead of multiple scattering
#include "G4CoulombScattering.hh"

// For Muonium formation in the Carbon foil
#include "sr1MuFormation.hh" // includes the yield function Y = Y(E).

// For a simple Muonium "scattering" when Mu hits solid materials
#include "sr1MuScatter.hh"


// Construction methods for the different processes

void sr1PhysicsList::ConstructEM()
{
  G4String myTypeOfProcesses  = sr1Parameters::GetInstance()->GetMyTypeOfProcesses();
  G4cout<<"   uuuuuuuuu myTypeOfProcesses="<<myTypeOfProcesses<<G4endl;
  
  G4bool inclMuProcesses = sr1Parameters::includeMuoniumProcesses;
  
  theParticleIterator->reset();
  while( (*theParticleIterator)() ){
    G4ParticleDefinition* particle = theParticleIterator->value();
    G4ProcessManager* pmanager = particle->GetProcessManager();
    G4String particleName = particle->GetParticleName();
    
    if (particleName == "gamma") {               // gamma
      if (myTypeOfProcesses=="lowenergy") {
	pmanager->AddDiscreteProcess(new G4LowEnergyPhotoElectric);
	pmanager->AddDiscreteProcess(new G4LowEnergyCompton);
	pmanager->AddDiscreteProcess(new G4LowEnergyGammaConversion);
	pmanager->AddDiscreteProcess(new G4LowEnergyRayleigh);
      }
      else if (myTypeOfProcesses=="penelope") {
	pmanager->AddDiscreteProcess(new G4PenelopePhotoElectric);
	pmanager->AddDiscreteProcess(new G4PenelopeCompton);
	pmanager->AddDiscreteProcess(new G4PenelopeGammaConversion);
	pmanager->AddDiscreteProcess(new G4PenelopeRayleigh);
      }
      else { // Standard G4 (default)
	pmanager->AddDiscreteProcess(new G4PhotoElectricEffect);
	pmanager->AddDiscreteProcess(new G4ComptonScattering);
	pmanager->AddDiscreteProcess(new G4GammaConversion);
      }

    }
    else if (particleName == "e-") {  //electron
      if (myTypeOfProcesses=="lowenergy") {
	pmanager->AddProcess(new G4MultipleScattering,     -1, 1,1);
	pmanager->AddProcess(new G4LowEnergyIonisation,    -1, 2,2);
	pmanager->AddProcess(new G4LowEnergyBremsstrahlung,-1,-1,3);    
      }
      else if (myTypeOfProcesses=="penelope") {
	pmanager->AddProcess(new G4MultipleScattering,     -1, 1, 1);
	pmanager->AddProcess(new G4PenelopeIonisation,     -1, 2, 2);
	pmanager->AddProcess(new G4PenelopeBremsstrahlung, -1,-1, 3);
      }
      else if (myTypeOfProcesses=="coulomb") {
	pmanager->AddDiscreteProcess(new G4CoulombScattering);
	pmanager->AddProcess(new G4eIonisation,       -1, 2,2);
	pmanager->AddProcess(new G4eBremsstrahlung,   -1, 3,3);      
      }
      else { // Standard G4 (default)
	pmanager->AddProcess(new G4MultipleScattering,-1, 1,1);
	pmanager->AddProcess(new G4eIonisation,       -1, 2,2);
	pmanager->AddProcess(new G4eBremsstrahlung,   -1, 3,3);      
      }
    }
    
    else if (particleName == "e+") {
      
      if (myTypeOfProcesses=="penelope") {
	pmanager->AddProcess(new G4MultipleScattering,     -1, 1, 1);
	pmanager->AddProcess(new G4PenelopeIonisation,     -1, 2, 2);
	pmanager->AddProcess(new G4PenelopeBremsstrahlung, -1,-1, 3);
	pmanager->AddProcess(new G4PenelopeAnnihilation,    0,-1, 4);
      }
      else if (myTypeOfProcesses=="coulomb") {
	pmanager->AddDiscreteProcess(new G4CoulombScattering);
	pmanager->AddProcess(new G4eIonisation,       -1, 2,2);
	pmanager->AddProcess(new G4eBremsstrahlung,   -1, 3,3);
	pmanager->AddProcess(new G4eplusAnnihilation,  0,-1,4);
      }
      else {  // Standard G4 (default)
	pmanager->AddProcess(new G4MultipleScattering,-1, 1,1);
	pmanager->AddProcess(new G4eIonisation,       -1, 2,2);
	pmanager->AddProcess(new G4eBremsstrahlung,   -1, 3,3);
	pmanager->AddProcess(new G4eplusAnnihilation,  0,-1,4);
      }
    }
    
    else if ((particleName=="mu+")||(particleName=="mu-")) {  //muon  
      
      if (myTypeOfProcesses=="coulomb") {
          pmanager->AddDiscreteProcess(new G4CoulombScattering);  // TS - Attention use Meyer!!
	  }
      else if (myTypeOfProcesses=="coulombAndMultiple") {
	  std::cout <<"musrPhysicsList:  Switching between Coulomb and Multiple scattering for mu+"<<std::endl;
	  pmanager->AddProcess(new G4MultipleScattering,-1, 1,1);
	  pmanager->AddDiscreteProcess(new G4CoulombScattering);
	  }
      else {
	  pmanager->AddProcess(new G4MultipleScattering,-1, 1,1); // TS - Attention use Meyer!!
	  }
	  
      if (inclMuProcesses){ // This accounts ONLY for Mu formation!
          G4VProcess* aMuoniumFormation = new sr1MuFormation();
	  pmanager->AddProcess(aMuoniumFormation);
	  pmanager->SetProcessOrdering(aMuoniumFormation,idxPostStep,2);
	  //pmanager->AddProcess(new sr1MuFormation, -1,-1,2); // TS Muonium formation process - shorthand
	  }
      
      pmanager->AddProcess(new G4MuIonisation,      -1, 2,2);
      pmanager->AddProcess(new G4MuBremsstrahlung,  -1, 3,3);
      pmanager->AddProcess(new G4MuPairProduction,  -1, 4,4);       
      //cks Change to Geant 4 default muon decay with spin (instead of using the code of Taofiq).
      //      pmanager->AddProcess(new sr1AtRestSpinRotation, 1, -1, -1);
      //      sr1Decay* theDecayProcess = new sr1Decay();
      //      pmanager->AddProcess(theDecayProcess);
      //      pmanager ->SetProcessOrderingToLast(theDecayProcess, idxAtRest);
      //      pmanager ->SetProcessOrdering(theDecayProcess, idxPostStep);
      //
      //      // The spin precession (rotation) at rest is already included in 
      //      // the G4DecayWithSpin();
      //      G4AtRestSpinPrecession* theSpinPrecession = new G4AtRestSpinPrecession();
      //      pmanager->AddRestProcess(theSpinPrecession);
      
      
      
      if (particle->GetParticleName()=="mu+"){
      G4DecayWithSpin* theDecayProcess = new G4DecayWithSpin();
      //sr1DecayWithSpin* theDecayProcess = new sr1DecayWithSpin();
      pmanager->AddProcess(theDecayProcess);
      pmanager->SetProcessOrderingToLast(theDecayProcess, idxAtRest);
      pmanager->SetProcessOrdering(theDecayProcess, idxPostStep);
      }
      }
 
    
    else if ((particle->GetParticleName()=="Mu")&&(inclMuProcesses)) { // other Mu processes
      G4cout<<"\nAccounting for Muonium formation and other Mu processes.\n"<<G4endl;
      //pmanager->AddProcess(new G4MultipleScattering, -1, 1, 1); 
      //pmanager->AddProcess(new sr1MuScatter, -1, 1, 1); <--- ERROR ?!?
      //pmanager->AddProcess(new sr1MuScatter, -1, -1, 1); //<--- ERROR ?!?
      
      // Muonium "scattering"
      G4VProcess* aMuScatt = new sr1MuScatter();
      pmanager->AddProcess(aMuScatt);
      pmanager->SetProcessOrdering(aMuScatt, idxPostStep, 1);
            
      // Muonium decay process
      G4Decay* theDecayProcess = new G4Decay();
      //sr1DecayWithSpin* theDecayProcess = new sr1DecayWithSpin();
      pmanager->AddProcess(theDecayProcess);
      pmanager->SetProcessOrderingToLast(theDecayProcess, idxAtRest);
      pmanager->SetProcessOrdering(theDecayProcess, idxPostStep);
      }
      
      
            
      /*
      else if (particle->GetParticleName()=="Mu")
	{
	  //!
	    Half of the muonium has an isotropic decay.
	  //
	  G4DecayTable* MuoniumDecayTable = new G4DecayTable();
	  MuoniumDecayTable -> Insert(new G4MuonDecayChannelWithSpin("Mu",0.5)); // ?? with spin? LEMuSRMuonDecayChannel("Mu",0.5));
	  MuoniumDecayTable -> Insert(new G4MuonDecayChannel("Mu",0.5)); 
	  MuoniumDecayTable ->DumpInfo();
	  MuoniumParticle::MuoniumDefinition() -> SetDecayTable(MuoniumDecayTable);
	  
	  G4ProcessManager* pmanager = particle->GetProcessManager();
	  if (theLEMuSRDecayProcess->IsApplicable(*particle)) 
	    { 
	      pmanager ->AddProcess(theLEMuSRDecayProcess);
	      pmanager ->AddProcess(aStepLimiter);
	      
	      // set ordering 
	      //	      pmanager ->SetProcessOrdering(aStepLimiter, idxPostStep, 1);  // not needed for the muonium
	      pmanager ->SetProcessOrderingToLast(theLEMuSRDecayProcess, idxPostStep);
	      pmanager ->SetProcessOrderingToLast(theLEMuSRDecayProcess, idxAtRest);
	    }

	      pmanager->DumpInfo();


	} */
      
      //csk 
    //}
    
    else if ((!particle->IsShortLived()) &&
	       (particle->GetPDGCharge() != 0.0) && 
	       (particle->GetParticleName() != "chargedgeantino")) {
      // All stable, charged particles except geantino
      pmanager->AddProcess(new G4MultipleScattering,     -1,1,1);
      if (myTypeOfProcesses=="lowenergy") {
	pmanager->AddProcess(new G4hLowEnergyIonisation, -1,2,2);
      }
      else { 
	pmanager->AddProcess(new G4hIonisation,          -1,2,2);
      }
      ///pmanager->AddProcess(new G4UserSpecialCuts,  -1,-1,3);      
    }
  }
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......



#include "G4Decay.hh"
void sr1PhysicsList::ConstructGeneral()
{
}


//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#include "G4Region.hh"
#include "G4RegionStore.hh"
#include "G4ProductionCuts.hh"

void sr1PhysicsList::SetCuts()
{
  //G4VUserPhysicsList::SetCutsWithDefault method sets 
  //the default cut value for all particle types 
  //
  SetCutsWithDefault();

  //cks  Set some cuts specific to a given region:
  //  G4Region* region;
  //  G4String regName;
  //  G4ProductionCuts* cuts;
  //  
  //  regName="Target_region";
  //  region = G4RegionStore::GetInstance()->GetRegion(regName);
  //  if (region==NULL) {
  //    G4cout << "sr1PhysicsList.cc: ERROR!  Region "<<regName<<" not found."<<G4endl;
  //    G4cout << "    Program is forced to stop!"<<G4endl;
  //    exit(EXIT_FAILURE);
  //  }
  //  cuts = new G4ProductionCuts;
  //
  //  cuts->SetProductionCut(0.005*mm,G4ProductionCuts::GetIndex("mu+"));
  //  cuts->SetProductionCut(0.005*mm,G4ProductionCuts::GetIndex("e-"));
  //  region->SetProductionCuts(cuts);
  //csk
    
   if (verboseLevel>0) DumpCutValuesTable();
   DumpCutValuesTable();
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......