Message: Neutron elastic problem Not Logged In (login)
 Next-in-Thread Next-in-Thread
 Next-in-Forum Next-in-Forum

Question Neutron elastic problem 

Keywords: neutron
Forum: Hadronic Processes
Date: 02 Mar, 2007
From: Justin Dingley <dinglj@rpi.edu>

Hello all,

I am currently simulating a beam of low energy (2.5 MeV) neutrons incident
upon a NE-213 detector. There should be a distribution of energy given 
to protons via elastic scattering, but what I'm getting is almost all 
the protons having max neutron energy.

I am using the LHEP_PRECO_HP library for hadron interactions and the 
low energy EM processes for everything else. Is there some peculiarity 
to the (n,p) interaction that might be responsible for this?

For reference, the following is the physics list I'm using:

exrdmPhysicsList::exrdmPhysicsList() : G4VModularPhysicsList()
{
  G4LossTableManager::Instance();
  defaultCutValue = 1.*mm;
  cutForGamma     = defaultCutValue;
  cutForElectron  = defaultCutValue;
  cutForPositron  = defaultCutValue;
  DetectorCuts = 0;
  TargetCuts   = 0;

  pMessenger = new exrdmPhysicsListMessenger(this);

  // Particles
  particleList = new G4DecayPhysics("decays");

  // EM physics
  emPhysicsList = new G4LowEmPhysics("lowEnergy");
}

exrdmPhysicsList::~exrdmPhysicsList()
{
  delete pMessenger;
  delete particleList;
  delete emPhysicsList;
  for(size_t i=0; i<hadronPhys.size(); i++) {
    delete hadronPhys[i];
  }
}

void exrdmPhysicsList::ConstructParticle()
{
  particleList->ConstructParticle();
}

void exrdmPhysicsList::ConstructProcess()
{
  AddTransportation();
  emPhysicsList->ConstructProcess();
  particleList->ConstructProcess();
  for(size_t i=0; i<hadronPhys.size(); i++) {
    hadronPhys[i]->ConstructProcess();
  }

  G4HadronProcessStore::Instance()->Dump(1);
}

void exrdmPhysicsList::AddPhysicsList(const G4String& name)
{
  if (verboseLevel>0)
    G4cout << "PhysicsList::AddPhysicsList: <" << name << ">" << G4endl;

  if (name == "LHEP") {

    hadronPhys.push_back( new G4EmExtraPhysics("gamma_nuc"));
    hadronPhys.push_back( new G4HadronElasticPhysics("LElastic",
						     verboseLevel,false));
    hadronPhys.push_back( new HadronPhysicsLHEP());
    hadronPhys.push_back( new G4IonPhysics("ion"));
    
  } else if (name == "LHEP_PRECO_HP") {

    hadronPhys.push_back( new G4EmExtraPhysics("gamma_nuc"));
    hadronPhys.push_back( new G4HadronElasticPhysics("LElastic",
						     verboseLevel,false));
    hadronPhys.push_back( new HadronPhysicsLHEP_PRECO_HP());
    hadronPhys.push_back( new G4QStoppingPhysics("stopping"));
    hadronPhys.push_back( new G4IonPhysics("ion"));
} else {

    G4cout << "PhysicsList::AddPhysicsList: <" << name << ">"
           << " is not defined"
           << G4endl;
  }
.........
}


void exrdmPhysicsList::SetStandardList(G4bool flagHP)
{
  hadronPhys.push_back( new G4EmExtraPhysics("gamma_nuc"));
  hadronPhys.push_back( new G4HadronElasticPhysics("elastic",verboseLevel,
						   flagHP));
  hadronPhys.push_back( new G4QStoppingPhysics("stopping"));
  hadronPhys.push_back( new G4IonBinaryCascadePhysics("binary_ion"));
}

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

void exrdmPhysicsList::SetCuts()
{

  if (verboseLevel >0){
    G4cout << "PhysicsList::SetCuts:";
    G4cout << "CutLength : " << G4BestUnit(defaultCutValue,"Length") << G4endl;
  }

  // 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(cutForGamma, "gamma");
  SetCutValue(cutForElectron, "e-");
  SetCutValue(cutForPositron, "e+");
  
  if( !TargetCuts ) SetTargetCut(cutForElectron);
  G4Region* region = (G4RegionStore::GetInstance())->GetRegion("Target");
  region->SetProductionCuts(TargetCuts);
  G4cout << "Target cuts are set" << G4endl;

  if( !DetectorCuts ) SetDetectorCut(cutForElectron);
  region = (G4RegionStore::GetInstance())->GetRegion("Detector");
  region->SetProductionCuts(DetectorCuts);
  G4cout << "Detector cuts are set" << G4endl;

  if (verboseLevel>0) DumpCutValuesTable();
}

void exrdmPhysicsList::SetCutForGamma(G4double cut)
{
  cutForGamma = cut;
  SetParticleCuts(cutForGamma, G4Gamma::Gamma());
}

void exrdmPhysicsList::SetCutForElectron(G4double cut)
{
  cutForElectron = cut;
  SetParticleCuts(cutForElectron, G4Electron::Electron());
}

void exrdmPhysicsList::SetCutForPositron(G4double cut)
{
  cutForPositron = cut;
  SetParticleCuts(cutForPositron, G4Positron::Positron());
}

void exrdmPhysicsList::SetTargetCut(G4double cut)
{
  if( !TargetCuts ) TargetCuts = new G4ProductionCuts();

  TargetCuts->SetProductionCut(cut, idxG4GammaCut);
  TargetCuts->SetProductionCut(cut, idxG4ElectronCut);
  TargetCuts->SetProductionCut(cut, idxG4PositronCut);

}

void exrdmPhysicsList::SetDetectorCut(G4double cut)
{
  if( !DetectorCuts ) DetectorCuts = new G4ProductionCuts();

  DetectorCuts->SetProductionCut(cut, idxG4GammaCut);
  DetectorCuts->SetProductionCut(cut, idxG4ElectronCut);
  DetectorCuts->SetProductionCut(cut, idxG4PositronCut);
}

Here is a histogram of energy deposited in the detector:

http://img525.imageshack.us/my.php?image=testyo8.jpg

Thank you in advance,

Justin Dingley
Rensselaer Polytechnic
Nuclear Engineering

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

1 None: Re: Neutron elastic problem   (Vladimir IVANTCHENKO - 04 Mar, 2007)
 Add Message Add Message
to: "Neutron elastic problem"

 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 ]