I am trying to simulate neutron capture process using Geant 4.3.1.
I chose exampleN03 as a base, with exampleN04 PhysicsList and other *Physcis
classes.  After changing the class name prefixes from ExN04 to ExN03, I have
assigned the following interaction to neutron in the HadronPhysics class so
that the neutron capture process would be implemented.

   G4HadronElasticProcess theElasticProcess1;
   G4LElastic* theElasticModel1;
   G4NeutronHPElastic*   theElasticNeutron;
   G4NeutronHPInelastic* theLENeutronInelasticModel;
   G4NeutronHPFission*   theLENeutronFissionModel;
   G4NeutronHPCapture*   theLENeutronCaptureModel;

  // Neutron
  pManager = G4Neutron::Neutron()->GetProcessManager();

    // elastic scattering

  theElasticModel1 = new G4LElastic;
  theElasticProcess1.RegisterMe(theElasticModel1);
  theElasticModel1->SetMinEnergy(19*MeV);

  theElasticNeutron = new G4NeutronHPElastic;
  theElasticProcess1.RegisterMe(theElasticNeutron);
  G4CrossSectionDataStore * theStore =
    ((G4HadronElasticProcess)theElasticProcess1).GetCrossSectionDataStore();
  G4NeutronHPElasticData * theNeutronData = new G4NeutronHPElasticData;
  theStore->AddDataSet(theNeutronData);

  pManager->AddDiscreteProcess(&theElasticProcess1);


    // inelastic scattering

  theLENeutronModel = new G4LENeutronInelastic();
  theHENeutronModel = new G4HENeutronInelastic();
  theLENeutronModel->SetMinEnergy(19*MeV);
  theNeutronInelastic.RegisterMe(theLENeutronModel);
  theNeutronInelastic.RegisterMe(theHENeutronModel);

  theLENeutronInelasticModel = new G4NeutronHPInelastic;
  theNeutronInelastic.RegisterMe(theLENeutronInelasticModel);
  G4CrossSectionDataStore * theStore1 =
    ((G4HadronInelasticProcess)theNeutronInelastic).GetCrossSectionDataStore();
  G4NeutronHPInelasticData * theNeutronData1 = new G4NeutronHPInelasticData;
  theStore1->AddDataSet(theNeutronData1);

  pManager->AddDiscreteProcess(&theNeutronInelastic);


     // fission

  theNeutronFissionModel = new G4LFission;
  theNeutronFissionModel->SetMinEnergy(19*MeV);
  theNeutronFission.RegisterMe(theNeutronFissionModel);

  theLENeutronFissionModel = new G4NeutronHPFission;
  theNeutronFission.RegisterMe(theLENeutronFissionModel);
  G4CrossSectionDataStore * theStore2 =
    ((G4HadronFissionProcess)theNeutronFission).GetCrossSectionDataStore();
  G4NeutronHPFissionData * theNeutronData2 = new G4NeutronHPFissionData;
  theStore2->AddDataSet(theNeutronData2);

  pManager->AddDiscreteProcess(&theNeutronFission);

     // capture

  theNeutronCaptureModel = new G4LCapture;
  theNeutronCaptureModel->SetMinEnergy(19*MeV);
  theNeutronCapture.RegisterMe(theNeutronCaptureModel);

  theLENeutronCaptureModel = new G4NeutronHPCapture;
  theNeutronCapture.RegisterMe(theLENeutronCaptureModel);
  G4CrossSectionDataStore * theStore3 =
    ((G4HadronCaptureProcess)theNeutronCapture).GetCrossSectionDataStore();
  G4NeutronHPCaptureData * theNeutronData3 = new G4NeutronHPCaptureData;
  theStore3->AddDataSet(theNeutronData3);

  pManager->AddDiscreteProcess(&theNeutronCapture);


When 1 eV - 100 keV neutrons were shot against quartz (SiO2), which was
prepared in exampleN03 DetectorConstruction class by default, my program ends
with segmentation fault before it emits any gamma-ray. Here is the gdb output:


Program received signal SIGSEGV, Segmentation fault.
#0  G4RWTPtrOrderedVector<G4ReactionProduct>::length (this=0x0)
    at /home/naoyuki/G4UG/geant4.3.1/source/global/STLInterface/g4rw/
tpordvec.icc:312
   312       return rwsize;

#1  0x80db1e4 in G4NeutronHPCapture::ApplyYourself (this=0x115877f8,
    aTrack=@0x8da2e78, aTargetNucleus=@0x8e0ea0c)
    at src/G4NeutronHPCapture.cc:134
   134       G4int nSecondaries = thePhotons->length;

#2  0x824dff5 in G4HadronicProcess::GeneralPostStepDoIt (this=0x8e0e78c,
    aTrack=@0x8da2e78, aStep=@0x8db1c38) at src/G4HadronicProcess.cc:76
#3  0x8205c0d in G4HadronCaptureProcess::PostStepDoIt (this=0x8e0e78c,
    aTrack=@0x8da2e78, aStep=@0x8db1c38)
    at include/G4HadronCaptureProcess.hh:76


When shot against Aluminium, the resultant gamma-ray energy spectrum is hard to
understand.

I haven't modified anything other than the physics lists, and I have set the
environment variables NeutronHPCrossSections properly.  What could cause my
problem?  Any comment or advise will be greatly appreciated.

Best Regards,

Naoyuki Yamashita
Waseda University