Message: Testing thermal neutron creation Not Logged In (login)
 Next-in-Thread Next-in-Thread
 Next-in-Forum Next-in-Forum

Feedback Testing thermal neutron creation 

Forum: Physics List
Date: 16 Oct, 2008
From: Sean Turnbull <Sean Turnbull>

Hello, I have tested a 3"x3" NaI detector and it's working. Now I am trying to simulate neutron capture in Na and I producing respective isotopes. I am using the HP models for capture. I am also firing thermal neutron energies at the NaI detector, around 0.025 ev from 10cm away from NaI target.

I am posting my PhysicsList in case I am missing something as I can't seem to produce any istopes or capture processes.

// // // PhysicsList.cc, NaI Project // #include "PhysicsList.hh"

#include "G4ParticleDefinition.hh"
#include "G4ParticleWithCuts.hh"
#include "G4ProcessManager.hh"
#include "G4ProcessVector.hh"
#include "G4ParticleTypes.hh"
#include "G4ParticleTable.hh"
#include "G4ios.hh"

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

PhysicsList::PhysicsList():  G4VUserPhysicsList()
{
  defaultCutValue = 1.0*mm;		// cut values
  cutForElectron  = 6.0*mm;
  cutForGamma     = 5.0*mm;
  cutForNeutron   = 1.0*mm;
  cutForMuon	  = defaultCutValue;
  SetVerboseLevel(1);
}

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

PhysicsList::~PhysicsList() { }

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

void PhysicsList::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. 

  ConstructBosons();
  ConstructLeptons();
  ConstructHadrons();

//  ConstructMesons();
//  ConstructBaryons();
//  ConstructIons();
  ConstructShortLiveds();
}

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

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

  // gamma
  G4Gamma::GammaDefinition();
}

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

void PhysicsList::ConstructLeptons()
{
  // leptons
  //  e+/-
  G4Electron::ElectronDefinition();
  G4Positron::PositronDefinition();
  // mu+/-
  G4MuonPlus::MuonPlusDefinition();
  G4MuonMinus::MuonMinusDefinition();
  // nu_e
  G4NeutrinoE::NeutrinoEDefinition();
  G4AntiNeutrinoE::AntiNeutrinoEDefinition();
  // nu_mu
  G4NeutrinoMu::NeutrinoMuDefinition();
  G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();
}

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

void PhysicsList::ConstructMesons()
{
  // mesons
  G4PionPlus::PionPlusDefinition();
  G4PionMinus::PionMinusDefinition();
  G4PionZero::PionZeroDefinition();
  G4KaonPlus::KaonPlusDefinition();
  G4KaonMinus::KaonMinusDefinition();
  G4Eta::EtaDefinition();
  G4EtaPrime::EtaPrimeDefinition();
  G4KaonZero::KaonZeroDefinition();
  G4AntiKaonZero::AntiKaonZeroDefinition();
  G4KaonZeroLong::KaonZeroLongDefinition();
  G4KaonZeroShort::KaonZeroShortDefinition();
}

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

void PhysicsList::ConstructBaryons()
{
  //  baryons
  G4Proton::ProtonDefinition();
  G4AntiProton::AntiProtonDefinition();

  G4Neutron::NeutronDefinition();
  G4AntiNeutron::AntiNeutronDefinition();
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PhysicsList::ConstructIons()
{
  G4Deuteron::DeuteronDefinition();
  G4Alpha::AlphaDefinition();
  G4GenericIon::GenericIonDefinition();
}
#include "G4MesonConstructor.hh"
#include "G4BaryonConstructor.hh"
#include "G4IonConstructor.hh"

void PhysicsList::ConstructHadrons()
{
 //  mesons
  G4MesonConstructor mConstructor;
  mConstructor.ConstructParticle();

  //  baryons
  G4BaryonConstructor bConstructor;
  bConstructor.ConstructParticle();

  //  ions
  G4IonConstructor iConstructor;
  iConstructor.ConstructParticle();

}
void PhysicsList::ConstructShortLiveds()
{
  // ShortLiveds
  ;
}
// Construct Processes ///////////////////////////////////////////////////    	       
void PhysicsList::ConstructProcess()
{
  AddTransportation();
  ConstructEM();
  ConstructHad();
  ConstructGeneral();
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... #include "G4MultipleScattering.hh"

//muon

#include "G4MuIonisation.hh"
#include "G4MuBremsstrahlung.hh"
#include "G4MuPairProduction.hh"
#include "G4MuonMinusCaptureAtRest.hh"
//gamma
#include "G4PhotoElectricEffect.hh" 
#include "G4ComptonScattering.hh" 
#include "G4GammaConversion.hh" 
#include "G4GammaConversionToMuons.hh"
// low energy
#include "G4LowEnergyPhotoElectric.hh"
#include "G4LowEnergyCompton.hh"
#include "G4LowEnergyRayleigh.hh" 
#include "G4LowEnergyGammaConversion.hh"
// e-
#include "G4LowEnergyIonisation.hh"
#include "G4LowEnergyBremsstrahlung.hh"
// e+
#include "G4eIonisation.hh"
#include "G4eBremsstrahlung.hh"
#include "G4eplusAnnihilation.hh"
//ion #include "G4ionIonisation.hh" // alpha and GenericIon and deuterons, triton, He3: #include "G4EnergyLossTables.hh" // this uses Ziegler 1988 and is the default - requires detailed testing // particularly with split between NuclearStopping and electron ionisation

#include "G4hIonisation.hh" // standard hadron ionisation

//Penelope Processes //#include "G4PenelopePhotoElectric.hh" //#include "G4PenelopeCompton.hh" //#include "G4PenelopeRayleigh.hh" //#include "G4PenelopeGammaConversion.hh"

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

void PhysicsList::ConstructEM() // low energy EM processes
{
  theParticleIterator->reset();
  while( (*theParticleIterator)() )
    {
      G4ParticleDefinition* particle = theParticleIterator->value();
      G4ProcessManager* pmanager = particle->GetProcessManager();
      G4String particleName = particle->GetParticleName();
      G4String particleType = particle->GetParticleType();
//      G4double charge = particle->GetPDGCharge();
      //processes

      if (particleName == "gamma") 
        { 
	  //gamma
//	  lowePhot = new G4LowEnergyPhotoElectric("LowEnPhotElec");    
//	  lowePhot->ActivateAuger(true);
//	  pmanager->AddDiscreteProcess(new G4LowEnergyRayleigh);
//	  pmanager->AddDiscreteProcess(lowePhot);
//	  pmanager->AddDiscreteProcess(new G4LowEnergyCompton);
//	  pmanager->AddDiscreteProcess(new G4LowEnergyGammaConversion);

	  pmanager->AddDiscreteProcess(new G4PhotoElectricEffect);
	  pmanager->AddDiscreteProcess(new G4ComptonScattering);
	  pmanager->AddDiscreteProcess(new G4GammaConversion);
	  pmanager->AddDiscreteProcess(new G4GammaConversionToMuons);

// pmanager->AddDiscreteProcess(new G4PenelopePhotoElectric()); // pmanager->AddDiscreteProcess(new G4PenelopeCompton()); // pmanager->AddDiscreteProcess(new G4PenelopeRayleigh()); // pmanager->AddDiscreteProcess(new G4PenelopeGammaConversion());

        }
        else if (particleName == "e-") {
          //electron
//	  loweIon  = new G4LowEnergyIonisation("LowEnergyIoni");
//    	  loweBrem = new G4LowEnergyBremsstrahlung("LowEnBrem");
//    	  loweBrem->SetAngularGenerator("tsai");

// pmanager->AddProcess(new G4MultipleScattering, -1, 1,1); // pmanager->AddProcess(loweIon, -1, 2,2); // pmanager->AddProcess(loweBrem, -1,-1,3);

	  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+") {
          // positron
//          pmanager->AddProcess(new G4MultipleScattering,-1, 1,1);
//          pmanager->AddProcess(new G4eIonisation,       -1, 2,2);
//          pmanager->AddProcess(new G4eBremsstrahlung,   -1,-1,3);// -1,3,3
//          pmanager->AddProcess(new G4eplusAnnihilation,  0,-1,4);

	 G4MultipleScattering* aMultipleScattering = new G4MultipleScattering();
         pmanager->AddProcess(aMultipleScattering,      -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  
          G4MultipleScattering* aMultipleScattering = new G4MultipleScattering();
          pmanager->AddProcess(aMultipleScattering,             -1, 1, 1);
          pmanager->AddProcess(new G4MuIonisation,      	-1, 2, 2);
          pmanager->AddProcess(new G4MuBremsstrahlung,  	-1,-1, 3);
          pmanager->AddProcess(new G4MuPairProduction,  	-1,-1, 4);
	  if(particleName == "mu-")
            pmanager->AddProcess(new G4MuonMinusCaptureAtRest(), 0,-1,-1);
	}
	else if( particleName == "GenericIon" ) 
	{
          pmanager->AddProcess(new G4MultipleScattering,-1, 1,1);
          pmanager->AddProcess(new G4ionIonisation,     -1, 2,2);
        } 
        else if ((!particle->IsShortLived()) &&
               (particle->GetPDGCharge() != 0.0));
        {
          //all others charged particles except geantino
          pmanager->AddProcess(new G4MultipleScattering,-1,1,1);
//	  pmanager->AddProcess(new G4hIonisation,       -1,2,2);
        } 
    }
}

// Hadronic processes ///////////////////////////////////////////////

// Elastic processes: #include "G4HadronElasticProcess.hh"

// Inelastic processes:

#include "G4PionPlusInelasticProcess.hh"
#include "G4PionMinusInelasticProcess.hh"
#include "G4KaonPlusInelasticProcess.hh"
#include "G4KaonZeroSInelasticProcess.hh"
#include "G4KaonZeroLInelasticProcess.hh"
#include "G4KaonMinusInelasticProcess.hh"
#include "G4ProtonInelasticProcess.hh"
#include "G4AntiProtonInelasticProcess.hh"
#include "G4NeutronInelasticProcess.hh"
#include "G4AntiNeutronInelasticProcess.hh"
#include "G4DeuteronInelasticProcess.hh"
#include "G4TritonInelasticProcess.hh"
#include "G4AlphaInelasticProcess.hh"

// Low-energy Models: < 20GeV

#include "G4LElastic.hh"
#include "G4LEPionPlusInelastic.hh"
#include "G4LEPionMinusInelastic.hh"
#include "G4LEKaonPlusInelastic.hh"
#include "G4LEKaonZeroSInelastic.hh"
#include "G4LEKaonZeroLInelastic.hh"
#include "G4LEKaonMinusInelastic.hh"
#include "G4LEProtonInelastic.hh"
#include "G4LEAntiProtonInelastic.hh"
#include "G4LENeutronInelastic.hh"
#include "G4LEAntiNeutronInelastic.hh"
#include "G4LEDeuteronInelastic.hh"
#include "G4LETritonInelastic.hh"
#include "G4LEAlphaInelastic.hh"

// High-energy Models: >20 GeV

#include "G4HEPionPlusInelastic.hh"
#include "G4HEPionMinusInelastic.hh"
#include "G4HEKaonPlusInelastic.hh"
#include "G4HEKaonZeroInelastic.hh"
#include "G4HEKaonZeroInelastic.hh"
#include "G4HEKaonMinusInelastic.hh"
#include "G4HEProtonInelastic.hh"
#include "G4HEAntiProtonInelastic.hh"
#include "G4HENeutronInelastic.hh"
#include "G4HEAntiNeutronInelastic.hh"

// Neutron high-precision models: <20 MeV

#include "G4NeutronHPElastic.hh"
#include "G4NeutronHPElasticData.hh"
#include "G4NeutronHPCapture.hh"
#include "G4NeutronHPCaptureData.hh"
#include "G4NeutronHPInelastic.hh"
#include "G4NeutronHPInelasticData.hh"
#include "G4LCapture.hh"

// Stopping processes

#include "G4PiMinusAbsorptionAtRest.hh"
#include "G4KaonMinusAbsorptionAtRest.hh"
#include "G4AntiProtonAnnihilationAtRest.hh"
#include "G4AntiNeutronAnnihilationAtRest.hh"

// Isotope production #include "G4NeutronIsotopeProduction.hh"

// ConstructHad() // Makes discrete physics processes for the hadrons, at present limited // to those particles with GHEISHA interactions (INTRC > 0). // The processes are: Elastic scattering and Inelastic scattering. // F.W.Jones 09-JUL-1998

void PhysicsList::ConstructHad()
{
  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
  G4LElastic* theElasticModel = new G4LElastic;
  theElasticProcess->RegisterMe(theElasticModel);

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

      if (particleName == "pi+")
        {
          pmanager->AddDiscreteProcess(theElasticProcess);
          G4PionPlusInelasticProcess* theInelasticProcess =
            new G4PionPlusInelasticProcess("inelastic");
          G4LEPionPlusInelastic* theLEInelasticModel =
            new G4LEPionPlusInelastic;
          theInelasticProcess->RegisterMe(theLEInelasticModel);
          G4HEPionPlusInelastic* theHEInelasticModel =
            new G4HEPionPlusInelastic;
          theInelasticProcess->RegisterMe(theHEInelasticModel);
          pmanager->AddDiscreteProcess(theInelasticProcess);
        }

      else if (particleName == "pi-")
        {
          pmanager->AddDiscreteProcess(theElasticProcess);
          G4PionMinusInelasticProcess* theInelasticProcess =
            new G4PionMinusInelasticProcess("inelastic");
          G4LEPionMinusInelastic* theLEInelasticModel =
            new G4LEPionMinusInelastic;
          theInelasticProcess->RegisterMe(theLEInelasticModel);
          G4HEPionMinusInelastic* theHEInelasticModel =
            new G4HEPionMinusInelastic;
          theInelasticProcess->RegisterMe(theHEInelasticModel);
          pmanager->AddDiscreteProcess(theInelasticProcess);
          G4String prcNam;
          pmanager->AddRestProcess(new G4PiMinusAbsorptionAtRest, ordDefault);
        }

      else if (particleName == "kaon+")
        {
          pmanager->AddDiscreteProcess(theElasticProcess);
          G4KaonPlusInelasticProcess* theInelasticProcess =
            new G4KaonPlusInelasticProcess("inelastic");
          G4LEKaonPlusInelastic* theLEInelasticModel =
            new G4LEKaonPlusInelastic;
          theInelasticProcess->RegisterMe(theLEInelasticModel);
          G4HEKaonPlusInelastic* theHEInelasticModel =
            new G4HEKaonPlusInelastic;
          theInelasticProcess->RegisterMe(theHEInelasticModel);
          pmanager->AddDiscreteProcess(theInelasticProcess);
        }

      else if (particleName == "kaon0S")
        {
          pmanager->AddDiscreteProcess(theElasticProcess);
          G4KaonZeroSInelasticProcess* theInelasticProcess =
            new G4KaonZeroSInelasticProcess("inelastic");
          G4LEKaonZeroSInelastic* theLEInelasticModel =
            new G4LEKaonZeroSInelastic;
          theInelasticProcess->RegisterMe(theLEInelasticModel);
          G4HEKaonZeroInelastic* theHEInelasticModel =
            new G4HEKaonZeroInelastic;
           theInelasticProcess->RegisterMe(theHEInelasticModel);
          pmanager->AddDiscreteProcess(theInelasticProcess);
        }

      else if (particleName == "kaon0L")
        {
          pmanager->AddDiscreteProcess(theElasticProcess);
          G4KaonZeroLInelasticProcess* theInelasticProcess =
            new G4KaonZeroLInelasticProcess("inelastic");
          G4LEKaonZeroLInelastic* theLEInelasticModel =
            new G4LEKaonZeroLInelastic;
          theInelasticProcess->RegisterMe(theLEInelasticModel);
          G4HEKaonZeroInelastic* theHEInelasticModel =
            new G4HEKaonZeroInelastic;
          theInelasticProcess->RegisterMe(theHEInelasticModel);
          pmanager->AddDiscreteProcess(theInelasticProcess);
        }

      else if (particleName == "kaon-")
        {
          pmanager->AddDiscreteProcess(theElasticProcess);
          G4KaonMinusInelasticProcess* theInelasticProcess =
            new G4KaonMinusInelasticProcess("inelastic");
          G4LEKaonMinusInelastic* theLEInelasticModel =
            new G4LEKaonMinusInelastic;
          theInelasticProcess->RegisterMe(theLEInelasticModel);
          G4HEKaonMinusInelastic* theHEInelasticModel =
            new G4HEKaonMinusInelastic;
          theInelasticProcess->RegisterMe(theHEInelasticModel);
          pmanager->AddDiscreteProcess(theInelasticProcess);
          pmanager->AddRestProcess(new G4KaonMinusAbsorptionAtRest, ordDefault);
        }

      else if (particleName == "proton")
        {
          pmanager->AddDiscreteProcess(theElasticProcess);
          G4ProtonInelasticProcess* theInelasticProcess =
            new G4ProtonInelasticProcess("inelastic");
          G4LEProtonInelastic* theLEInelasticModel = new G4LEProtonInelastic;
          theInelasticProcess->RegisterMe(theLEInelasticModel);
          G4HEProtonInelastic* theHEInelasticModel = new G4HEProtonInelastic;
          theInelasticProcess->RegisterMe(theHEInelasticModel);
          pmanager->AddDiscreteProcess(theInelasticProcess);
        }

      else if (particleName == "anti_proton")
        {
          pmanager->AddDiscreteProcess(theElasticProcess);
          G4AntiProtonInelasticProcess* theInelasticProcess =
            new G4AntiProtonInelasticProcess("inelastic");
          G4LEAntiProtonInelastic* theLEInelasticModel =
            new G4LEAntiProtonInelastic;
          theInelasticProcess->RegisterMe(theLEInelasticModel);
          G4HEAntiProtonInelastic* theHEInelasticModel =
            new G4HEAntiProtonInelastic;
          theInelasticProcess->RegisterMe(theHEInelasticModel);
          pmanager->AddDiscreteProcess(theInelasticProcess);
        }

      else if (particleName == "neutron") {
        // elastic scattering
        G4HadronElasticProcess* theNeutronElasticProcess =
          new G4HadronElasticProcess;
        G4LElastic* theElasticModel1 = new G4LElastic;
        G4NeutronHPElastic * theElasticNeutron = new G4NeutronHPElastic;
        theNeutronElasticProcess->RegisterMe(theElasticModel1);
        theElasticModel1->SetMinEnergy(0.00001*eV);
        theNeutronElasticProcess->RegisterMe(theElasticNeutron);
        G4NeutronHPElasticData * theNeutronData = new G4NeutronHPElasticData;
        theNeutronElasticProcess->AddDataSet(theNeutronData);
        pmanager->AddDiscreteProcess(theNeutronElasticProcess);
        // inelastic scattering
        G4NeutronInelasticProcess* theInelasticProcess =
          new G4NeutronInelasticProcess("inelastic");
        G4LENeutronInelastic* theInelasticModel = new G4LENeutronInelastic;
        theInelasticModel->SetMinEnergy(0.00001*eV);
        theInelasticProcess->RegisterMe(theInelasticModel);
        G4NeutronHPInelastic * theLENeutronInelasticModel =
          new G4NeutronHPInelastic;
        theInelasticProcess->RegisterMe(theLENeutronInelasticModel);
        G4NeutronHPInelasticData * theNeutronData1 =
          new G4NeutronHPInelasticData;
        theInelasticProcess->AddDataSet(theNeutronData1);
        pmanager->AddDiscreteProcess(theInelasticProcess);
        // capture
        G4HadronCaptureProcess* theCaptureProcess =
          new G4HadronCaptureProcess;
        G4LCapture* theCaptureModel = new G4LCapture;
        theCaptureModel->SetMinEnergy(0.00001*eV);
        theCaptureProcess->RegisterMe(theCaptureModel);
        G4NeutronHPCapture * theLENeutronCaptureModel = new G4NeutronHPCapture;
        theCaptureProcess->RegisterMe(theLENeutronCaptureModel);
        G4NeutronHPCaptureData * theNeutronData3 = new G4NeutronHPCaptureData;
        theCaptureProcess->AddDataSet(theNeutronData3);
        pmanager->AddDiscreteProcess(theCaptureProcess);
        //  G4ProcessManager* pmanager = G4Neutron::Neutron->GetProcessManager(
      //  pmanager->AddProcess(new G4UserSpecialCuts(),-1,-1,1);

      // isotope production
        G4NeutronIsotopeProduction * aProductionModel = new
	G4NeutronIsotopeProduction;
	theCaptureProcess->RegisterIsotopeProductionModel(aProductionModel);
	//theInelasticProcess->RegisterIsotopeProductionModel(aProductionModel);
	G4HadronicProcess::EnableIsotopeProductionGlobally();
      }
      else if (particleName == "anti_neutron")
        {
          pmanager->AddDiscreteProcess(theElasticProcess);
          G4AntiNeutronInelasticProcess* theInelasticProcess =
            new G4AntiNeutronInelasticProcess("inelastic");
          G4LEAntiNeutronInelastic* theLEInelasticModel =
            new G4LEAntiNeutronInelastic;
          theInelasticProcess->RegisterMe(theLEInelasticModel);
          G4HEAntiNeutronInelastic* theHEInelasticModel =
            new G4HEAntiNeutronInelastic;
          theInelasticProcess->RegisterMe(theHEInelasticModel);
          pmanager->AddDiscreteProcess(theInelasticProcess);
        }

      else if (particleName == "deuteron")
        {
          pmanager->AddDiscreteProcess(theElasticProcess);
          G4DeuteronInelasticProcess* theInelasticProcess =
            new G4DeuteronInelasticProcess("inelastic");
          G4LEDeuteronInelastic* theLEInelasticModel =
            new G4LEDeuteronInelastic;
          theInelasticProcess->RegisterMe(theLEInelasticModel);
          pmanager->AddDiscreteProcess(theInelasticProcess);
        }

      else if (particleName == "triton")
        {
          pmanager->AddDiscreteProcess(theElasticProcess);
          G4TritonInelasticProcess* theInelasticProcess =
            new G4TritonInelasticProcess("inelastic");
          G4LETritonInelastic* theLEInelasticModel =
            new G4LETritonInelastic;
          theInelasticProcess->RegisterMe(theLEInelasticModel);
          pmanager->AddDiscreteProcess(theInelasticProcess);
        }

      else if (particleName == "alpha")
        {
          pmanager->AddDiscreteProcess(theElasticProcess);
          G4AlphaInelasticProcess* theInelasticProcess =
            new G4AlphaInelasticProcess("inelastic");
          G4LEAlphaInelastic* theLEInelasticModel =
            new G4LEAlphaInelastic;
          theInelasticProcess->RegisterMe(theLEInelasticModel);
          pmanager->AddDiscreteProcess(theInelasticProcess);
        }
    }
}

// Decays ///////////////////////////////////////////////////////////////////

#include "G4Decay.hh"
#include "G4RadioactiveDecay.hh"
#include "G4IonTable.hh"
#include "G4Ions.hh"

void PhysicsList::ConstructGeneral() {

  // Add Decay Process
  G4Decay* theDecayProcess = new G4Decay();
  theParticleIterator->reset();
  while( (*theParticleIterator)() )
    {
      G4ParticleDefinition* particle = theParticleIterator->value();
      G4ProcessManager* pmanager = particle->GetProcessManager();

      if (theDecayProcess->IsApplicable(*particle) && !particle->IsShortLived())
        {
          pmanager ->AddProcess(theDecayProcess);
          // set ordering for PostStepDoIt and AtRestDoIt
          pmanager ->SetProcessOrdering(theDecayProcess, idxPostStep);
          pmanager ->SetProcessOrdering(theDecayProcess, idxAtRest);
        }
    }

  // Declare radioactive decay to the GenericIon in the IonTable.
  const G4IonTable *theIonTable =
    G4ParticleTable::GetParticleTable()->GetIonTable();
  G4RadioactiveDecay *theRadioactiveDecay = new G4RadioactiveDecay();

  for (G4int i=0; i<theIonTable->Entries(); i++)
    {
      G4String particleName = theIonTable->GetParticle(i)->GetParticleName();
      G4String particleType = theIonTable->GetParticle(i)->GetParticleType();

      if (particleName == "GenericIon")
        {
          G4ProcessManager* pmanager =
            theIonTable->GetParticle(i)->GetProcessManager();
          pmanager->SetVerboseLevel(VerboseLevel);
          pmanager ->AddProcess(theRadioactiveDecay);
          pmanager ->SetProcessOrdering(theRadioactiveDecay, idxPostStep);
          pmanager ->SetProcessOrdering(theRadioactiveDecay, idxAtRest);
        }
    }
}

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

void PhysicsList::SetCuts()

  //G4VUserPhysicsList::SetCutsWithDefault method sets 
  //the default cut value for all particle types 
  //
{
  if (verboseLevel >0){
    G4cout << "PhysicsList::SetCuts:";
    G4cout << "CutLength : " << defaultCutValue/mm << " (mm)" << G4endl;
    G4cout << "CutLength : " << G4BestUnit(defaultCutValue,"Length") << G4endl;
  }

//special for low energy physics
 // G4double lowlimit=250*eV;
  G4double lowlimit=0.00150*eV;
  G4ProductionCutsTable::GetProductionCutsTable()->SetEnergyRange(lowlimit,100.*GeV);

  // set cut values for gamma at first and for e- second and next for e+,
  // because some processes for e+/e- need cut values for gamma 
  SetCutValue(cutForMuon, "mu-");
  SetCutValue(cutForGamma, "gamma");
  SetCutValue(cutForElectron, "e-");
  SetCutValue(cutForNeutron, "neutron");
  SetCutValue(cutForElectron, "e+");

  if (verboseLevel>0) DumpCutValuesTable();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

Thanks!!

Inline Depth:
 1 1
 All All
Outline Depth:
 1 1
 2 2
 All All
Add message: (add)

1 More: Re: Testing thermal neutron creation   (Sean Turnbull - 21 Oct, 2008)
 Add Message Add Message
to: "Testing thermal neutron creation"

 Subscribe Subscribe

This site runs SLAC HyperNews version 1.11-slac-98, derived from the original HyperNews


[ Geant 4 Home | Geant 4 HyperNews | Search | Request New Forum | Feedback ]