rest-for-physics / restG4

A generic Geant4 code project that uses the REST Geant4 library to define the simulation conditions and collect the simulated results as a ROOT file, in the form of REST c++ objects. restG4 just requires two inputs, and RML file defining the simulation conditions and a GDML geometry file.
3 stars 5 forks source link

How to deactivate the neutron killer (nKiller) process? #118

Closed tsiledak closed 10 months ago

tsiledak commented 1 year ago

thermaln.rml.txt geometry.gdml.txt setup.gdml.txt

Dear experts,

I am sorry for bothering again... I am still puzzled to find why I have so few energy deposition entries in an Argon gas from a thermal neutron capture in LiF target where I expect alphas and tritons. The example based on the 08Alpha from the directory is executed as:

export REST_FOIL=0.35
export REST_ENERGY=0.000000025 restG4 thermaln.rml

for 100000 events and for a LiF foil of 0.35 microns thick.

For neutrons I have:

                       Hadronic Processes for neutron

Process: neutronInelastic Model: QGSP: 12 GeV ---> 100 TeV Model: FTFP: 3 GeV ---> 25 GeV Model: Binary Cascade: 19.9 MeV ---> 6 GeV Model: NeutronHPInelastic: 0 eV ---> 20 MeV Cr_sctns: NeutronHPInelasticXS: 0 eV ---> 20 MeV Cr_sctns: G4NeutronInelasticXS: 0 eV ---> 100 TeV

Process: nCapture Model: NeutronHPCapture: 0 eV ---> 20 MeV Model: nRadCapture: 19.9 MeV ---> 100 TeV Cr_sctns: NeutronHPCaptureXS: 0 eV ---> 20 MeV Cr_sctns: G4NeutronCaptureXS: 0 eV ---> 100 TeV

Process: nFission Model: NeutronHPFission: 0 eV ---> 20 MeV Model: G4LFission: 19.9 MeV ---> 100 TeV Cr_sctns: NeutronHPFissionXS: 0 eV ---> 20 MeV Cr_sctns: ZeroXS: 0 eV ---> 100 TeV

Process: hadElastic Model: hElasticCHIPS: 19.5 MeV ---> 100 TeV Model: NeutronHPElastic: 0 eV ---> 20 MeV Cr_sctns: NeutronHPElasticXS: 0 eV ---> 20 MeV Cr_sctns: G4NeutronElasticXS: 0 eV ---> 100 TeV

Process: nKiller


It seems that this process kill the very big majority of the neutrons before they are captured and give alphas and tritons... As a result, out of 100000 events I have only 21 entries!

I implement the same geometry using the standalone Geant4 (same version) and the physics list: QGSP_BIC_HP. Running the same number of events (the nKiller is absent) results in the following pict where I have 10 times more entries (in MeV). g4-run

I think that this nKiller process might related to the physics list used? Any way to deactivate it in the main rml file? Something like: /process/inactivate nKiller after the initialization?

Thanks a lot in advance for your time

jgalan commented 1 year ago

Hi Giorgos I am sorry you are experiencing problems with this issue. And I think it is also a relevant question for ourselves.

I have found the following Geant4 post: https://geant4-forum.web.cern.ch/t/deactivation-of-the-neutron-killer-nkiller-process/6125

It seems there is a way to disable the process using macros, perhaps this might give a clue on how to disable it on restG4? @lobis?

Perhaps one could do a quick test by adding the following line inside src/Application.cxx ?

UI->ApplyCommand("/process/inactivate nKiller");

Just before the line:

        UI->ApplyCommand("/run/beamOn " + to_string(nEvents));
jgalan commented 1 year ago

Then, if that solves the issue it would be nice to understand how the physics list decides that this process is enabled/disabled.

tsiledak commented 1 year ago

Hi Javier, thanks so much for your comments. I have implemented as u suggested as follows:

In the ~/rest-framework/source/packages/restG4/src$ I edited the Application.cxx by adding: if (nEvents > 0) // batch mode { UI->ApplyCommand("/tracking/verbose 0"); UI->ApplyCommand("/run/initialize"); UI->ApplyCommand("/physics_engine/neutron/energyLimit 0. GeV"); UI->ApplyCommand("/physics_engine/neutron/timeLimit 10. s"); UI->ApplyCommand("/process/inactivate nKiller"); UI->ApplyCommand("/run/beamOn " + to_string(nEvents)); }

either the inactivate line or the two lines above or all the 3 together...

Then I compiled as: gt226504@dphnpct97:~/rest-framework/build$ make [ 23%] Built target RestFramework [ 23%] Built target restManager [ 24%] Built target restRoot [ 37%] Built target RestAxion [ 44%] Built target RestTrack [ 64%] Built target RestDetector [ 72%] Built target RestGeant4 [ 86%] Built target RestRaw [ 90%] Built target RestConnectors [ 92%] Built target RestLegacy [ 94%] Built target RestWimp Consolidate compiler generated dependencies of target RestRestG4 [ 94%] Building CXX object source/packages/restG4/CMakeFiles/RestRestG4.dir/src/Application.cxx.o [ 94%] Linking CXX shared library libRestRestG4.so [ 97%] Built target RestRestG4 Consolidate compiler generated dependencies of target restG4 [ 97%] Building CXX object source/packages/restG4/CMakeFiles/restG4.dir/src/Application.cxx.o [ 97%] Linking CXX executable restG4 [100%] Built target restG4 gt226504@dphnpct97:~/rest-framework/build$

and then re-run the simulation BUT always get the same very low number of events and the nKiller process is always present! In standalone Geant4 using G4:ftfp-bert-hp or even QGSP_BIC_HP I do not observe such behavior.... Or might be anything wrong in the input rml file I am using?

lobis commented 1 year ago

@tsiledak I think the nKiller process is due to the fact you are using <physicsList name="G4NeutronTrackingCut"/>.

I think this physic list introduces a process that kills neutron track after some time in order to speed up simulations. Try to remove this line from your rml and let us know if the problem still exists.

tsiledak commented 1 year ago

Hi, I have just removed the line you proposed from the thermaln.rml file but unfortunately I am getting still the same number of events :-( Only 22 entries instead of 10 times more as expected from fluka or geant used as a standalone code... But at least the line with the process nKiller is not showed up on the screen:


                       Hadronic Processes for neutron

Process: neutronInelastic Model: QGSP: 12 GeV ---> 100 TeV Model: FTFP: 3 GeV ---> 25 GeV Model: Binary Cascade: 19.9 MeV ---> 6 GeV Model: NeutronHPInelastic: 0 eV ---> 20 MeV Cr_sctns: NeutronHPInelasticXS: 0 eV ---> 20 MeV Cr_sctns: G4NeutronInelasticXS: 0 eV ---> 100 TeV

Process: nCapture Model: NeutronHPCapture: 0 eV ---> 20 MeV Model: nRadCapture: 19.9 MeV ---> 100 TeV Cr_sctns: NeutronHPCaptureXS: 0 eV ---> 20 MeV Cr_sctns: G4NeutronCaptureXS: 0 eV ---> 100 TeV

Process: nFission Model: NeutronHPFission: 0 eV ---> 20 MeV Model: G4LFission: 19.9 MeV ---> 100 TeV Cr_sctns: NeutronHPFissionXS: 0 eV ---> 20 MeV Cr_sctns: ZeroXS: 0 eV ---> 100 TeV

Process: hadElastic Model: hElasticCHIPS: 19.5 MeV ---> 100 TeV Model: NeutronHPElastic: 0 eV ---> 20 MeV Cr_sctns: NeutronHPElasticXS: 0 eV ---> 20 MeV Cr_sctns: G4NeutronElasticXS: 0 eV ---> 100 TeV


Screenshot from 2023-09-20 20-12-17

You could see in the picture attached the factor of 10 between restG4 (10M events) and geant4 (1M events)... This was from a run before deleting the line you proposed that anyway doesnt change the results...

jgalan commented 1 year ago

Hi Giorgos, just for making sure, the GDML geometry you are using is exactly the same one? same material definition?

Sensitive volume properties:
    - Material: G4_Ar
    - Temperature: 293.15 K
    - Density: 0.00166201 g/cm3

Could you share the Geant4 standalone code to try to reproduce it myself? Thanks!

tsiledak commented 1 year ago

Hi Javier, I have just shared a folder in google drive at your email where the geant standalone piece of code is. Let me know if you got it since I did not receive any confirmation mail in my sent messages... Anyway, the geometry I used is the same with the exception that the volumes of LiF and Cu are cylinders of the same thickness as in GDML in rest... But since I use a point thermal neutron source in the centre this does not play role. I confirmed the results using FLUKA with exactly the same geometry as defined in rest test and the results are always a factor of 10 more statistics... This is what I am getting with standalone geant4 for 100000 events (with rest I have only ~21 entries). I am grateful for your help. In case you did not receive the shared code as I pointed out before, let me know to re-send it...

Best regards G4run

Material: G4_Ar density: 1.662 mg/cm3 RadL: 117.621 m Nucl.Int.Length: 719.888 m
Imean: 188.000 eV temperature: 293.15 K pressure: 1.00 atm

jgalan commented 10 months ago

Hi Giorgos, I am really not sure about the problem here. But I thought it might be an issue with units on GDML? I suspect because we had problems in previous occasions, but in principle it was only affecting visualization.

If I increase the size of the foil to 3.5um, then do you think that matches the official Geant4/FLUKA result for 0.35um?

What do you think? @lobis?

tsiledak commented 10 months ago

Thanks so much Javier, So you think it might be a wrong interpretation of the units used... that might explain the very low efficiency... I will upgrade and re-run it according to your suggestions and let u know... I am really grateful for the invitations to be a member of rest-for-physics.

Best regards

Georgios

tsiledak commented 10 months ago

Thanks so much Javier, So you think it might be a wrong interpretation of the units used... that might explain the very low efficiency... I will upgrade and re-run it according to your suggestions and let u know... I am really grateful for the invitations to be a member of rest-for-physics.

Best regards

Georgios

tsiledak commented 10 months ago

Dear Javier,

I tried to re-run the simulation with 3.5um instead of 0.35um. Screenshot from 2023-11-17 16-02-06 Unfortunately it does not match the distribution as expected by FLUKA or GEANT4 neither the shape nor the number of entries... In the attached pict is the result with FLUKA (top plot) and REST (bottom plot) that were run for 1M thermal neutrons and 3.5um thick LiF... I have no idea why... Is it the physics list? but I am using the same version of GEANT4 and can reproduce with the standalone version the same with FLUKA result...

Best regards

Georgios