#include "MyExN01PhysicsList.hh"
#include "G4ParticleTypes.hh"
#include "G4ParticleTable.hh"
#include "globals.hh"
#include "G4ParticleDefinition.hh"
#include "G4ProcessManager.hh"
#include "G4ProcessVector.hh"
#include "G4ios.hh"
#include <iomanip>  
#include "G4UserSpecialCuts.hh" 
#include "G4StepLimiter.hh"

MyExN01PhysicsList::MyExN01PhysicsList()
{;}

MyExN01PhysicsList::~MyExN01PhysicsList()
{;}

void MyExN01PhysicsList::ConstructParticle()
{
  // In this method, static member functions should be called
  // for all particles which you want to use.
  // This ensures that objects of these particle types will be
  // created in the program. 

  ConstructIons();
  ConstructBaryons();
  ConstructLeptons();
  ConstructBosons();
  ConstructMesons();
}

#include "G4IonConstructor.hh"

void MyExN01PhysicsList::ConstructIons()
{ //Vou construir TODOS os Ions!!!
    G4IonConstructor iConstructor;
    iConstructor.ConstructParticle();
}

void MyExN01PhysicsList::ConstructBosons()
{
  // gamma
  G4Gamma::GammaDefinition();

  // optical photon
  G4OpticalPhoton::OpticalPhotonDefinition();
}

#include "G4LeptonConstructor.hh"
void MyExN01PhysicsList::ConstructLeptons()
{
  // Construct all leptons
  G4LeptonConstructor pConstructor;
  pConstructor.ConstructParticle();
}

#include "G4MesonConstructor.hh"
void MyExN01PhysicsList::ConstructMesons()
{
  //  Construct all mesons
  G4MesonConstructor pConstructor;
  pConstructor.ConstructParticle();
}

#include "G4BaryonConstructor.hh"
void MyExN01PhysicsList::ConstructBaryons()
{
  //  Construct all barions
  G4BaryonConstructor  pConstructor;
  pConstructor.ConstructParticle();  
}

void MyExN01PhysicsList::ConstructProcess()
{
  // Define transportation process
  AddTransportation();
  ConstructEM();
  ConstructHadronic();
  //  G4ProcessManager* pManager= G4GenericIon::GenericIon()->GetProcessManager();
}

void MyExN01PhysicsList::AddTransportation()
{
  G4VUserPhysicsList::AddTransportation();
}

#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 "G4MultipleScattering.hh"
#include "G4ionIonisation.hh"
void MyExN01PhysicsList::ConstructEM()
{
  theParticleIterator->reset();
  while( (*theParticleIterator)() ){
    G4ParticleDefinition* particle = theParticleIterator->value();
    G4ProcessManager* pmanager = particle->GetProcessManager();
    G4String particleName = particle->GetParticleName();
     
    if (particleName == "gamma") {
    // gamma
      // Construct processes for gamma
      pmanager->AddDiscreteProcess(new G4GammaConversion());
      pmanager->AddDiscreteProcess(new G4ComptonScattering());      
      pmanager->AddDiscreteProcess(new G4PhotoElectricEffect());

    } else if (particleName == "e-") {
    //electron
      // Construct processes for electron
      G4VProcess* theeminusMultipleScattering = new G4MultipleScattering();
      G4VProcess* theeminusIonisation = new G4eIonisation();
      G4VProcess* theeminusBremsstrahlung = new G4eBremsstrahlung();
      // add processes
      pmanager->AddProcess(theeminusMultipleScattering);
      pmanager->AddProcess(theeminusIonisation);
      pmanager->AddProcess(theeminusBremsstrahlung);      
      //step limit
      pmanager->AddProcess(new G4StepLimiter,       -1,-1,3);         
      //cuts
      pmanager->AddProcess(new G4UserSpecialCuts,   -1,-1,4);  

      // set ordering for AlongStepDoIt
      pmanager->SetProcessOrdering(theeminusMultipleScattering, idxAlongStep,  1);
      pmanager->SetProcessOrdering(theeminusIonisation, idxAlongStep,  2);
      // set ordering for PostStepDoIt
      pmanager->SetProcessOrdering(theeminusMultipleScattering, idxPostStep, 1);
      pmanager->SetProcessOrdering(theeminusIonisation, idxPostStep, 2);
      pmanager->SetProcessOrdering(theeminusBremsstrahlung, idxPostStep, 3);

    } else if (particleName == "e+") {
    //positron
      // Construct processes for positron
      G4VProcess* theeplusMultipleScattering = new G4MultipleScattering();
      G4VProcess* theeplusIonisation = new G4eIonisation();
      G4VProcess* theeplusBremsstrahlung = new G4eBremsstrahlung();
      G4VProcess* theeplusAnnihilation = new G4eplusAnnihilation();
      // add processes
      pmanager->AddProcess(theeplusMultipleScattering);
      pmanager->AddProcess(theeplusIonisation);
      pmanager->AddProcess(theeplusBremsstrahlung);
      pmanager->AddProcess(theeplusAnnihilation);
      // set ordering for AtRestDoIt
      pmanager->SetProcessOrderingToFirst(theeplusAnnihilation, idxAtRest);
      // set ordering for AlongStepDoIt
      pmanager->SetProcessOrdering(theeplusMultipleScattering, idxAlongStep,  1);
      pmanager->SetProcessOrdering(theeplusIonisation, idxAlongStep,  2);
      // set ordering for PostStepDoIt
      pmanager->SetProcessOrdering(theeplusMultipleScattering, idxPostStep, 1);
      pmanager->SetProcessOrdering(theeplusIonisation, idxPostStep, 2);
      pmanager->SetProcessOrdering(theeplusBremsstrahlung, idxPostStep, 3);
      pmanager->SetProcessOrdering(theeplusAnnihilation, idxPostStep, 4);
  
    } else if( particleName == "mu+" || 
               particleName == "mu-"    ) {
    //muon  
     // Construct processes for muon+
     G4VProcess* aMultipleScattering = new G4MultipleScattering();
     G4VProcess* aBremsstrahlung = new G4MuBremsstrahlung();
     G4VProcess* aPairProduction = new G4MuPairProduction();
     G4VProcess* anIonisation = new G4MuIonisation();
      // add processes
     pmanager->AddProcess(anIonisation);
     pmanager->AddProcess(aMultipleScattering);
     pmanager->AddProcess(aBremsstrahlung);
     pmanager->AddProcess(aPairProduction);
     // set ordering for AlongStepDoIt
     pmanager->SetProcessOrdering(aMultipleScattering, idxAlongStep,  1);
     pmanager->SetProcessOrdering(anIonisation, idxAlongStep,  2);
     // set ordering for PostStepDoIt
     pmanager->SetProcessOrdering(aMultipleScattering, idxPostStep, 1);
     pmanager->SetProcessOrdering(anIonisation, idxPostStep, 2);
     pmanager->SetProcessOrdering(aBremsstrahlung, idxPostStep, 3);
     pmanager->SetProcessOrdering(aPairProduction, idxPostStep, 4);
     
   } else if ((!particle->IsShortLived()) &&
	      (particle->GetPDGCharge() != 0.0) && 
	      (particle->GetParticleName() != "chargedgeantino")) {
     // all others charged particles except geantino
     G4VProcess* aMultipleScattering = new G4MultipleScattering();
     G4VProcess* anIonisation = new G4hIonisation();
     // add processes
     pmanager->AddProcess(anIonisation);
     pmanager->AddProcess(aMultipleScattering);
     // set ordering for AlongStepDoIt
     pmanager->SetProcessOrdering(aMultipleScattering, idxAlongStep,  1);
     pmanager->SetProcessOrdering(anIonisation, idxAlongStep,  2);
     // set ordering for PostStepDoIt
     pmanager->SetProcessOrdering(aMultipleScattering, idxPostStep, 1);
     pmanager->SetProcessOrdering(anIonisation, idxPostStep, 2);
   } else if( particleName == "deuteron" ||
              particleName == "triton" ||
              particleName == "He3" ||
              particleName == "alpha" ||
	      particleName == "GenericIon"){
     pmanager->AddProcess(new G4MultipleScattering, -1, 1, 1);
     pmanager->AddProcess(new G4ionIonisation, -1, 2, 2);
   }
  }
}

 // Neutron high-precision models: <20 MeV
#include "G4NeutronHPElastic.hh"
#include "G4NeutronHPElasticData.hh"
#include "G4NeutronHPorLElasticModel.hh"
#include "G4NeutronHPCapture.hh"
#include "G4NeutronHPCaptureData.hh"
#include "G4NeutronHPorLCapture.hh"
#include "G4NeutronHPFission.hh"
#include "G4NeutronHPFissionData.hh"
#include "G4NeutronHPorLFission.hh"
#include "G4NeutronInelasticProcess.hh"
#include "G4NeutronHPInelastic.hh"
#include "G4NeutronHPInelasticData.hh"
#include "G4NeutronHPorLEInelastic.hh"
#include "G4NeutronHPorLEInelasticData.hh"
#include "G4NeutronHPThermalScatteringData.hh"
#include "G4NeutronHPThermalScattering.hh"
#include "G4NeutronHPJENDLHEData.hh"
#include "G4LCapture.hh"

#include "G4Decay.hh"
#include "G4WilsonAbrasionModel.hh"
#include "G4WilsonAblationModel.hh"
#include "G4HadronElasticProcess.hh"
#include "G4LElastic.hh"
#include "G4IonInelasticProcess.hh"
#include "G4TripathiCrossSection.hh"
#include "G4IonsShenCrossSection.hh"
void MyExN01PhysicsList::ConstructHadronic()
{
  //decay for short lived particles
  G4Decay* theDecayProcess = new G4Decay();
  theParticleIterator->reset();

  // Construct processes
  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
  G4LElastic* theElasticModel = new G4LElastic;
  theElasticProcess->RegisterMe(theElasticModel);

  // Model for coherent elastic proton-proton or neutron-neutron scattering. Valid for particles: proton, neutron, energy range: 0 - 1200 GeV
  /*
      G4LEpp* elasticModel = new G4LEpp();
      G4HadronElasticProcess* hadElastProc = new G4HadronElasticProcess();
      hadElastProc->RegisterMe(elasticModel);
      //      pManager->AddDiscreteProcess(hadElastProc);
  
      //Fission process: This process handles neutron-induced fission.
      G4HadronFissionProcess* fissionProc = new G4HadronFissionProcess();
      G4HadronicProcess::AddDataSet(new G4HadronFissionDataSet);
      //      neutronProcessManager->AddDiscreteProcess(fissionProc);
      */
      //Capture process: This process handles neutron capture.
  /*G4LCapture* captureModel = new G4LCapture();
  G4HadronCaptureProcess* capProcess = new G4HadronCaptureProcess();
  capProcess->RegisterMe(captureModel);*/
      //      neutronProcessManager->AddDiscreteProcess(capProcess);

  while( (*theParticleIterator)() ){
    G4ParticleDefinition* particle = theParticleIterator->value();
    G4ProcessManager* pmanager = particle->GetProcessManager();
    G4String particleName = particle->GetParticleName();

    if( particle->GetParticleName() == "neutron" ) {
      // elastic scattering                                                   
      G4HadronElasticProcess* NeutronElasticProcess = new G4HadronElasticProcess;                                             

      G4NeutronHPElasticData* HPElasticData = new 
      G4NeutronHPElasticData();

      G4NeutronHPJENDLHEData* JENDLHEElasticData = new 
      G4NeutronHPJENDLHEData();

      G4NeutronHPorLElasticModel* NeutronElasticModel = new G4NeutronHPorLElasticModel();
      NeutronElasticProcess->AddDataSet(NeutronElasticModel->GiveHPXSectionDataSet());
      NeutronElasticProcess->AddDataSet(HPElasticData);
      NeutronElasticProcess->AddDataSet(JENDLHEElasticData);
      NeutronElasticModel->SetMinEnergy(4.0*eV);
      NeutronElasticProcess->RegisterMe(NeutronElasticModel);
      G4NeutronHPThermalScatteringData* HPThermalScatteringData = new
      G4NeutronHPThermalScatteringData();
      NeutronElasticProcess->AddDataSet(HPThermalScatteringData);
      G4NeutronHPThermalScattering* NeutronThermalElasticModel = new
      G4NeutronHPThermalScattering();
      NeutronThermalElasticModel->SetMaxEnergy(4.0*eV);
      NeutronElasticProcess->RegisterMe(NeutronThermalElasticModel); 
      pmanager->AddDiscreteProcess(NeutronElasticProcess);

      //Inelastic scattering

      G4NeutronInelasticProcess* NeutronInelasticProcess = new G4NeutronInelasticProcess();
      G4NeutronHPInelasticData* HPInelasticData = new
      G4NeutronHPInelasticData();
      G4NeutronHPorLEInelastic* NeutronInelasticModel = new G4NeutronHPorLEInelastic();
      NeutronInelasticProcess->AddDataSet(NeutronInelasticModel->GiveXSectionDataSet());
      NeutronInelasticProcess->AddDataSet(HPInelasticData);
      NeutronInelasticProcess->RegisterMe(NeutronInelasticModel);
      pmanager->AddDiscreteProcess(NeutronInelasticProcess);
      

      /*G4LElastic* theElasticModel1 = new G4LElastic;                           
      G4NeutronHPElastic * theElasticNeutron = new G4NeutronHPElastic;         
      theNeutronElasticProcess->RegisterMe(theElasticModel1);                
      theElasticModel1->SetMinEnergy(19*MeV);  
      theNeutronElasticProcess->RegisterMe(theElasticNeutron);                 
      G4NeutronHPElasticData * theNeutronData = new G4NeutronHPElasticData;    
      theNeutronElasticProcess->AddDataSet(theNeutronData);                    
      pmanager->AddDiscreteProcess(theNeutronElasticProcess);     */

      G4HadronCaptureProcess* NeutronCapture = new G4HadronCaptureProcess();
      G4NeutronHPCaptureData* HPCaptureData = new G4NeutronHPCaptureData();
      G4NeutronHPorLCapture*  NeutronCapModel = new G4NeutronHPorLCapture();
      NeutronCapture->AddDataSet(NeutronCapModel->GiveXSectionDataSet());
      NeutronCapture->AddDataSet(HPCaptureData);
      NeutronCapture->RegisterMe(NeutronCapModel); 
      pmanager->AddDiscreteProcess(NeutronCapture);

      G4HadronFissionProcess* NeutronFission = new G4HadronFissionProcess();
      G4NeutronHPorLFission* FissionModel = new G4NeutronHPorLFission();
      G4NeutronHPFissionData* FissionData = new G4NeutronHPFissionData();
      NeutronFission->AddDataSet(FissionModel->GiveXSectionDataSet());
      NeutronFission->AddDataSet(FissionData);
      NeutronFission->RegisterMe(FissionModel);
      pmanager->AddDiscreteProcess(NeutronFission);


      /*G4HadronCaptureProcess* neutroncapture = new G4HadronCaptureProcess();
      G4LCapture* captureModel = new G4LCapture();
      G4NeutronHPCapture* neutroncapProcess = new G4NeutronHPCapture();
      neutroncapture->RegisterMe(neutroncapProcess);
      //      capProcess->RegisterMe(captureModel);
      pmanager->AddDiscreteProcess(neutroncapture);*/
      

      // set ordering for AlongStepDoIt
      //      pmanager->SetProcessOrdering(capProcess, idxAlongStep,  1);
      //      pmanager->SetProcessOrdering(theNeutronElasticProcess, idxAlongStep,  2);
      // set ordering for PostStepDoIt
      //      pmanager->SetProcessOrdering(capProcess, idxPostStep, 1);
      //      pmanager->SetProcessOrdering(theNeutronElasticProcess, idxPostStep, 2);
    }   

    //Elastic processes for hadrons and generic ions: Model for elastic hadron scattering, valid for all energy range 
    else if ( particle->GetParticleName() == "anti_neutron" ||
	      particle->GetParticleName() == "proton"       ||
	      particle->GetParticleName() == "anti_proton"  ||
	      particle->GetParticleName() == "deuteron"     ||
	      particle->GetParticleName() == "triton"       ||
	      particle->GetParticleName() == "alpha"        || 
	      particle->GetParticleName() == "He3"          ||
	      particle->GetParticleName() == "GenericIon"   || 
	      (particle->GetParticleType() == "nucleus" && particle->GetPDGCharge() != 0.)) {
      pmanager->AddDiscreteProcess(theElasticProcess);
    }
    
    if( particleName == "GenericIon"){

	//Inelastic
	// Construct processes
	//Wilson abrasion/ablation model for generic ions
      G4WilsonAblationModel* AblModel = new G4WilsonAblationModel();
      G4WilsonAbrasionModel* AbrModel = new G4WilsonAbrasionModel(AblModel);
      G4IonInelasticProcess* ionInelProcess = new G4IonInelasticProcess();
      AbrModel->SetMinEnergy(0.0);
      ionInelProcess->RegisterMe(AbrModel);
         // add processes
      pmanager->AddDiscreteProcess(ionInelProcess);
	// Load Cross Sections .. Last loaded is used first for energy range
	// Shen 0-10 GeV/nuc .. Tripathi 0-1 GeV/nuc
      G4IonsShenCrossSection* shenCS = new G4IonsShenCrossSection();
      ionInelProcess->AddDataSet(shenCS);
      G4TripathiCrossSection* tripathiData = new G4TripathiCrossSection();
      ionInelProcess->AddDataSet(tripathiData);
    }else  if (theDecayProcess->IsApplicable(*particle)){
      // add processes
      pmanager->AddProcess(theDecayProcess);
      // set ordering for PostStepDoIt
      pmanager->SetProcessOrdering(theDecayProcess, idxPostStep);
      pmanager->SetProcessOrdering(theDecayProcess, idxAtRest);
    } 
  }
}


void MyExN01PhysicsList::SetCuts()
{
  // uppress error messages even in case e/gamma/proton do not exist            
  G4int temp = GetVerboseLevel();                                                SetVerboseLevel(0);                                                           
  //  " G4VUserPhysicsList::SetCutsWithDefault" method sets 
  //   the default cut value for all particle types 
  SetCutsWithDefault();   
  //  G4ProcessManager* pmanager = G4Electron::Electron->GetProcessManager();
  /*  theParticleIterator->reset();
  while( (*theParticleIterator)() ){
    G4ParticleDefinition* particle = theParticleIterator->value();
    G4ProcessManager* pmanager = particle->GetProcessManager();
    G4String particleName = particle->GetParticleName();

    if( particle->GetParticleName() == "electron" ){
      //step limit
      pmanager->AddProcess(new G4StepLimiter,       -1,-1,3);         
      //cuts
      G4UserSpecialCuts* meuCut = new G4UserSpecialCuts();
      pmanager->AddProcess(meuCut,-1,-1,1);
      //      pmanager->SetProcessOrdering(meuCut, idxAlongStep,  1);
      ///pmanager->AddProcess(new G4UserSpecialCuts,   -1,-1,4);  
    }else
        SetCutsWithDefault();   
	}*/
  // Retrieve verbose level
  SetVerboseLevel(temp);  
}

/*#include "G4ComptonScattering.hh"
#include "G4GammaConversion.hh"
#include "G4PhotoElectricEffect.hh"

#include "G4MultipleScattering.hh"

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

#include "G4hIonisation.hh"
void ExN05PhysicsList::ConstructEM()
{
  theParticleIterator->reset();
  while( (*theParticleIterator)() ){
    G4ParticleDefinition* particle = theParticleIterator->value();
    G4ProcessManager* pmanager = particle->GetProcessManager();
    G4String particleName = particle->GetParticleName();
     
    if (particleName == "gamma") {
    // gamma
      // Construct processes for gamma
      pmanager->AddDiscreteProcess(new G4GammaConversion());
      pmanager->AddDiscreteProcess(new G4ComptonScattering());      
      pmanager->AddDiscreteProcess(new G4PhotoElectricEffect());

    } else if (particleName == "e-") {
    //electron
      // Construct processes for electron
      G4VProcess* theeminusMultipleScattering = new G4MultipleScattering();
      G4VProcess* theeminusIonisation = new G4eIonisation();
      G4VProcess* theeminusBremsstrahlung = new G4eBremsstrahlung();
      // add processes
      pmanager->AddProcess(theeminusMultipleScattering);
      pmanager->AddProcess(theeminusIonisation);
      pmanager->AddProcess(theeminusBremsstrahlung);      
      // set ordering for AlongStepDoIt
      pmanager->SetProcessOrdering(theeminusMultipleScattering, idxAlongStep,  1);
      pmanager->SetProcessOrdering(theeminusIonisation, idxAlongStep,  2);
      // set ordering for PostStepDoIt
      pmanager->SetProcessOrdering(theeminusMultipleScattering, idxPostStep, 1);
      pmanager->SetProcessOrdering(theeminusIonisation, idxPostStep, 2);
      pmanager->SetProcessOrdering(theeminusBremsstrahlung, idxPostStep, 3);

    } else if (particleName == "e+") {
    //positron
      // Construct processes for positron
      G4VProcess* theeplusMultipleScattering = new G4MultipleScattering();
      G4VProcess* theeplusIonisation = new G4eIonisation();
      G4VProcess* theeplusBremsstrahlung = new G4eBremsstrahlung();
      G4VProcess* theeplusAnnihilation = new G4eplusAnnihilation();
      // add processes
      pmanager->AddProcess(theeplusMultipleScattering);
      pmanager->AddProcess(theeplusIonisation);
      pmanager->AddProcess(theeplusBremsstrahlung);
      pmanager->AddProcess(theeplusAnnihilation);
      // set ordering for AtRestDoIt
      pmanager->SetProcessOrderingToFirst(theeplusAnnihilation, idxAtRest);
      // set ordering for AlongStepDoIt
      pmanager->SetProcessOrdering(theeplusMultipleScattering, idxAlongStep,  1);
      pmanager->SetProcessOrdering(theeplusIonisation, idxAlongStep,  2);
      // set ordering for PostStepDoIt
      pmanager->SetProcessOrdering(theeplusMultipleScattering, idxPostStep, 1);
      pmanager->SetProcessOrdering(theeplusIonisation, idxPostStep, 2);
      pmanager->SetProcessOrdering(theeplusBremsstrahlung, idxPostStep, 3);
      pmanager->SetProcessOrdering(theeplusAnnihilation, idxPostStep, 4);
  
    }
  }
}*/