I've been trying to turn the thermal neutron treatment on with the FTFP_BERT_HP physics list. I chose as a test case a monoenergetic, isotropic 2.5 MeV neutron source at the center of a 10 cm diameter polyethylene sphere. I then histogram the energy spectrum of the neutrons that escape the sphere.

To do so, I've stared with a modified version of Basic example B4a. I modified the DetectorConstruction and PrimaryGeneratorAction, and created my own version of the FTFP_BERT_HP physics list (called FTFP_BERT_HP_Thermal) where the only modifications are a modified version of G4HadronElasticPhysicsHP (called G4HadronElasticPhysicsHP_Thermal). The only modification to G4HadronElasticPhysicsHP is in ConstructProcess() where I've added a few lines. Here it is:

void G4HadronElasticPhysicsHP_Thermal::ConstructProcess()
{
	if (wasActivated)
	{
		return;
	}
	wasActivated = true;

	mainElasticBuilder->ConstructProcess();
	mainElasticBuilder->GetNeutronModel()->SetMinEnergy(19.5*MeV);

	G4HadronicProcess* hel = mainElasticBuilder->GetNeutronProcess();

	G4NeutronHPElastic* hp = new G4NeutronHPElastic();
	hp->SetMinEnergy(4.0*eV);
	hel->RegisterMe(hp);
	hel->AddDataSet(new G4NeutronHPElasticData());

	G4NeutronHPThermalScattering* theThermalModel = new G4NeutronHPThermalScattering();
	theThermalModel->SetMaxEnergy(4.0*eV);
	hel->RegisterMe(theThermalModel);
	hel->AddDataSet(new G4NeutronHPThermalScatteringData());

	if (verbose > 1)
	{
		G4cout << "### G4HadronElasticPhysicsHP_Thermal is constructed " << G4endl;
	}
}

When I histogram the energies of the escaping neutrons, I get the attached histogram (in blue). Also plotted is the result I get using the standard FTFP_BERT_HP physics list (magenta), as well as the results of the same geometry with MCNP 1.6 with the thermal treatment 'MT poly.01t' (black) and with no thermal treatment (orange.) Note that the upper vertical scale limit is cutting off much of the un-scattered component. I've attached the same plot with a logarithmic y-axis to show everything.

I'm obviously doing something wrong here, as the FTFP_BERT_HP_Thermal plot has a crazy jump at 4.0 eV.

What is the best way to add G4NeutronHPThermalScattering to an existing physics list? I've looked back in the forum here, and it seems that the way physics lists are constructed has changed significantly over the past couple of years so I'm not sure the older suggestions (or the examples in the PowerPoint slides from tutorials) still apply?

In short, what am I doing wrong?

Any help the community can give would be greatly appreciated.