diff --git a/docs/rmg-commands.rst b/docs/rmg-commands.rst index 904943b..5d5dca4 100644 --- a/docs/rmg-commands.rst +++ b/docs/rmg-commands.rst @@ -14,6 +14,7 @@ Command directory path : /RMG/ * */RMG/Processes/*: ``Commands for controlling physics processes`` * */RMG/Geometry/*: ``Commands for controlling geometry definitions`` * */RMG/Generator/*: ``Commands for controlling generators`` + * */RMG/GrabmayrGammaCascades/*: ``Control Peters gamma cascade model`` * */RMG/Confinement/*: ``...Title not available...`` Command directory path : /RMG/Manager/ @@ -340,6 +341,7 @@ Commands for controlling physics processes * *EnableGammaAngularCorrelation*: ``Set correlated gamma emission flag`` * *GammaTwoJMAX*: ``Set max 2J for sampling of angular correlations`` * *StoreICLevelData*: ``Store e- internal conversion data`` + * *UseGrabmayrsGammaCascades*: ``Use custom RMGNeutronCapture to apply Grabmayrs gamma cascades.`` Command /RMG/Processes/Realm ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -422,6 +424,15 @@ Store e- internal conversion data * **Omittable**: ``False`` * **Candidates**: ``0 1`` +Command /RMG/Processes/UseGrabmayrsGammaCascades +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Use custom RMGNeutronCapture to apply Grabmayrs gamma cascades. + +* **Parameter**: ``value`` + * **Parameter type**: ``b`` + * **Omittable**: ``False`` + Command directory path : /RMG/Processes/Stepping/ ------------------------------------------------- @@ -1112,6 +1123,48 @@ Set Z length * **Default value**: ``cm`` * **Candidates**: ``pc km m cm mm um nm Ang fm parsec kilometer meter centimeter millimeter micrometer nanometer angstrom fermi`` +Command directory path : /RMG/GrabmayrGammaCascades/ +---------------------------------------------------- + +Control Peters gamma cascade model + +* **Commands**: + * *SetGammaCascadeRandomStartLocation*: ``Set the whether the start location in the gamma cascade file is random or not`` + * *SetGammaCascadeFile*: ``Set a gamma cascade file for neutron capture on a specified isotope`` + +Command /RMG/GrabmayrGammaCascades/SetGammaCascadeRandomStartLocation +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Set the whether the start location in the gamma cascade file is random or not + +0 = don't + +1 = do + +* **Parameter**: ``arg0`` + * **Parameter type**: ``i`` + * **Omittable**: ``False`` + * **Default value**: ``0`` + * **Candidates**: ``0 1`` + +Command /RMG/GrabmayrGammaCascades/SetGammaCascadeFile +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Set a gamma cascade file for neutron capture on a specified isotope + +* **Parameter**: ``Z`` +* Z of isotope + * **Parameter type**: ``i`` + * **Omittable**: ``False`` +* **Parameter**: ``A`` +* A of isotope + * **Parameter type**: ``i`` + * **Omittable**: ``False`` +* **Parameter**: ``file`` +* /path/to/file of gamma cascade + * **Parameter type**: ``s`` + * **Omittable**: ``False`` + Command directory path : /RMG/Confinement/ ------------------------------------------ diff --git a/examples/06-NeutronCapture/CMakeLists.txt b/examples/06-NeutronCapture/CMakeLists.txt new file mode 100644 index 0000000..f0e87ba --- /dev/null +++ b/examples/06-NeutronCapture/CMakeLists.txt @@ -0,0 +1,10 @@ +cmake_minimum_required(VERSION 3.8) +project(06-NeutronCapture) + +set(CMAKE_BUILD_TYPE RelWithDebInfo) + +find_package(remage) + +add_executable(06-NeutronCapture main.cc HPGeTestStand.hh HPGeTestStand.cc IsotopeOutputScheme.cc + IsotopeOutputScheme.hh) +target_link_libraries(06-NeutronCapture PUBLIC RMG::remage) diff --git a/examples/06-NeutronCapture/GammaCascades/156Gd-100Entries.txt b/examples/06-NeutronCapture/GammaCascades/156Gd-100Entries.txt new file mode 100644 index 0000000..3983682 --- /dev/null +++ b/examples/06-NeutronCapture/GammaCascades/156Gd-100Entries.txt @@ -0,0 +1,124 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% listing of gamma cascades by MAURINA/GAMMOC +%% +%% GD155+N KRK-D, 3LO GAMMA-RAY STRENGTH FUNCTIONS +%% MATNUM: 15511 +%% E_n= 5 keV +%% +%% date: 20221101 time: 212412.007/pgrabmayr +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% En: neutron energy [keV] +%% Ex: excitation energy [keV] +%% M: multiplicity of gamma cascade; M_max=20 +%% Em: missing energy [keV] +%% :=0 ground state after capture +%% :=E_iso isomeric state after capture +%% & M=M_max :=E_left left out due too many gammas +%% Eg1: energy of first photon [keV] +%% Eg2: energy of 2nd photon [keV] ..... EgM +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% variables: +%% En/keV Ex/kev M Em/keV Eg1/keV Eg2/kev Eg3/kev etc.. EgM +%% format : +%% i6 i6 i3 i6 i6 i6 i6 i6 i6 i6 i6 i6 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + 5 8540 6 0 1895 1065 2233 1980 1278 89 + 5 8540 5 0 992 5156 1025 1278 89 + 5 8540 4 0 2765 2985 2701 89 + 5 8540 3 0 3360 3938 1242 + 5 8540 4 0 1800 3485 3166 89 + 5 8540 1 0 8540 + 5 8540 4 0 5310 1988 1153 89 + 5 8540 3 0 3262 2361 2917 + 5 8540 4 0 1700 5711 1040 89 + 5 8540 4 0 1618 3282 3551 89 + 5 8540 6 0 1516 3285 1805 780 1065 89 + 5 8540 4 0 4097 2252 1037 1154 + 5 8540 4 0 3600 3664 1187 89 + 5 8540 4 0 6359 939 1153 89 + 5 8540 5 0 2704 3576 940 1231 89 + 5 8540 5 0 1589 4711 1111 1040 89 + 5 8540 4 0 2732 3591 2128 89 + 5 8540 4 0 3620 2087 2744 89 + 5 8540 4 0 2979 3217 977 1367 + 5 8540 6 0 1928 1191 3910 1223 199 89 + 5 8540 4 0 3135 2227 2024 1154 + 5 8540 4 0 3611 2976 711 1242 + 5 8540 7 0 2320 1918 901 1389 736 1187 89 + 5 8540 3 0 4459 3992 89 + 5 8540 7 0 1337 3023 2352 580 960 199 89 + 5 8540 3 0 2500 2641 3399 + 5 8540 4 0 5705 565 2181 89 + 5 8540 5 0 2306 1640 3318 1187 89 + 5 8540 3 0 6239 2212 89 + 5 8540 4 0 5317 2935 199 89 + 5 8540 5 0 1792 3511 1939 1209 89 + 5 8540 7 0 1840 2958 765 1125 604 1159 89 + 5 8540 4 0 5441 1857 1153 89 + 5 8540 6 0 1356 2053 842 2922 1278 89 + 5 8540 7 0 2423 2431 1734 704 960 199 89 + 5 8540 3 0 3822 4629 89 + 5 8540 5 0 4878 1345 1075 1153 89 + 5 8540 5 0 3986 3278 988 199 89 + 5 8540 6 0 3901 2609 675 1067 199 89 + 5 8540 4 0 3287 3053 2111 89 + 5 8540 3 0 4254 4197 89 + 5 8540 5 0 3808 2629 1815 199 89 + 5 8540 5 0 2217 1396 3679 1159 89 + 5 8540 5 0 1073 4355 2824 199 89 + 5 8540 5 0 2100 2251 1998 2102 89 + 5 8540 3 0 2912 942 4686 + 5 8540 4 0 2298 4875 1278 89 + 5 8540 3 0 6034 2417 89 + 5 8540 4 0 2833 3855 698 1154 + 5 8540 5 0 3252 2173 1760 1266 89 + 5 8540 6 0 2682 1768 2156 567 1278 89 + 5 8540 4 0 1295 4192 2964 89 + 5 8540 6 0 5452 1043 534 263 1159 89 + 5 8540 4 0 1984 4435 2032 89 + 5 8540 5 0 2259 4053 861 1278 89 + 5 8540 4 0 2560 3722 2169 89 + 5 8540 5 0 2910 1971 2383 1187 89 + 5 8540 5 0 4593 1844 949 1065 89 + 5 8540 4 0 5404 1894 1153 89 + 5 8540 6 0 4208 1509 1565 970 199 89 + 5 8540 5 0 2864 2649 1898 1040 89 + 5 8540 5 0 4617 1953 694 1187 89 + 5 8540 4 0 3722 2793 486 1539 + 5 8540 7 0 2771 1472 2122 877 1010 199 89 + 5 8540 4 0 4131 3133 1187 89 + 5 8540 7 0 2752 1428 1495 1242 1335 199 89 + 5 8540 6 0 3119 1518 1643 893 1278 89 + 5 8540 7 0 2827 1275 1391 1578 1181 199 89 + 5 8540 6 0 1430 781 2105 1832 1150 1242 + 5 8540 5 0 2368 4069 949 1065 89 + 5 8540 6 0 3002 1668 2331 1251 199 89 + 5 8540 4 0 2263 1231 2076 2970 + 5 8540 4 0 3676 3544 1231 89 + 5 8540 5 0 3082 2850 1288 1231 89 + 5 8540 3 0 5109 1284 2147 + 5 8540 3 0 3614 3684 1242 + 5 8540 5 0 4121 2075 977 1278 89 + 5 8540 6 0 2928 3760 497 1067 199 89 + 5 8540 4 0 4300 2920 1231 89 + 5 8540 5 0 1874 3372 2036 1169 89 + 5 8540 5 0 3941 2587 692 1231 89 + 5 8540 5 0 2498 2964 1924 1065 89 + 5 8540 2 0 6393 2147 + 5 8540 4 0 1970 3725 1691 1154 + 5 8540 3 0 6218 2233 89 + 5 8540 5 0 4378 879 2041 1153 89 + 5 8540 3 0 4760 1334 2446 + 5 8540 7 0 4015 2495 407 1038 297 199 89 + 5 8540 4 0 3966 3254 1231 89 + 5 8540 5 0 5017 1743 532 1159 89 + 5 8540 4 0 1502 5718 1231 89 + 5 8540 2 0 4525 4015 + 5 8540 3 0 2571 3646 2323 + 5 8540 6 0 1595 1430 2344 1050 2032 89 + 5 8540 4 0 3475 3698 1278 89 + 5 8540 5 0 5540 1107 1605 199 89 + 5 8540 5 0 3397 2873 1141 1040 89 + 5 8540 8 0 1504 1512 3349 767 823 297 199 89 + 5 8540 4 0 4060 1143 3248 89 + 5 8540 4 0 1739 5481 1231 89 diff --git a/examples/06-NeutronCapture/GammaCascades/158Gd-100Entries.txt b/examples/06-NeutronCapture/GammaCascades/158Gd-100Entries.txt new file mode 100644 index 0000000..605cae2 --- /dev/null +++ b/examples/06-NeutronCapture/GammaCascades/158Gd-100Entries.txt @@ -0,0 +1,124 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% listing of gamma cascades by MAURINA/GAMMOC +%% +%% GD157+N KRK-D 3LO GAMMA-RAY STRENGTH (PAPER) +%% MATNUM: 15711 +%% E_n= 5 keV +%% +%% date: 20221218 time: 180456.105/pgrabmayr +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%$$%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% En: neutron energy [keV] +%% Ex: excitation energy [keV] +%% M: multiplicity of gamma cascade; M_max=20 +%% Em: missing energy [keV] +%% :=0 ground state after capture +%% :=E_iso isomeric state after capture +%% & M=M_max :=E_left left out due too many gammas +%% Eg1: energy of first photon [keV] +%% Eg2: energy of 2nd photon [keV] ..... EgM +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% variables: +%% En/keV Ex/kev M Em/keV Eg1/keV Eg2/kev Eg3/kev etc.. EgM +%% format : +%% i6 i6 i3 i6 i6 i6 i6 i6 i6 i6 i6 i6 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + 5 7941 4 0 1830 2292 3739 80 + 5 7941 3 0 2443 5418 80 + 5 7941 5 0 3856 2107 936 962 80 + 5 7941 4 0 4174 3506 181 80 + 5 7941 5 0 2329 3345 1290 897 80 + 5 7941 4 0 4112 2787 962 80 + 5 7941 2 0 5657 2284 + 5 7941 4 0 1868 3075 2918 80 + 5 7941 5 0 4050 1771 1078 962 80 + 5 7941 3 0 6917 944 80 + 5 7941 5 0 2026 3967 924 944 80 + 5 7941 4 0 3559 3358 944 80 + 5 7941 3 0 3700 3264 977 + 5 7941 5 0 1380 2895 1730 1856 80 + 5 7941 4 0 2047 4631 1183 80 + 5 7941 3 0 5299 2562 80 + 5 7941 4 0 2829 2787 1348 977 + 5 7941 4 0 2154 1296 1992 2499 + 5 7941 4 0 3389 2759 816 977 + 5 7941 4 0 2298 4203 1360 80 + 5 7941 3 0 3574 3104 1263 + 5 7941 5 0 1108 5014 777 962 80 + 5 7941 5 0 2931 3062 924 944 80 + 5 7941 5 0 3905 2243 751 962 80 + 5 7941 4 0 2423 4494 944 80 + 5 7941 5 0 2874 2742 1301 944 80 + 5 7941 4 0 1830 3901 2130 80 + 5 7941 5 0 1264 5100 1316 181 80 + 5 7941 5 0 3627 2911 361 962 80 + 5 7941 5 0 2307 3309 1301 944 80 + 5 7941 5 0 2063 3116 1738 944 80 + 5 7941 5 0 2831 1939 2147 944 80 + 5 7941 5 0 1774 4809 1097 181 80 + 5 7941 3 0 2827 4137 977 + 5 7941 3 0 3166 2880 1895 + 5 7941 6 0 3618 2504 660 898 181 80 + 5 7941 4 0 2502 3786 466 1187 + 5 7941 4 0 3139 3166 1556 80 + 5 7941 2 0 5601 2340 + 5 7941 4 0 1493 2501 3867 80 + 5 7941 4 0 3311 3367 1183 80 + 5 7941 4 0 2970 3335 1556 80 + 5 7941 6 0 1471 1680 2772 1757 181 80 + 5 7941 5 0 2770 1285 2446 1360 80 + 5 7941 4 0 5062 2618 181 80 + 5 7941 5 0 1017 2320 2555 1072 977 + 5 7941 6 0 1884 4103 540 1153 181 80 + 5 7941 4 0 3438 3526 897 80 + 5 7941 4 0 3897 3020 944 80 + 5 7941 4 0 1094 2447 4320 80 + 5 7941 5 0 2051 2014 2852 944 80 + 5 7941 6 0 1895 1904 2983 898 181 80 + 5 7941 4 0 1794 2136 2751 1260 + 5 7941 5 0 1767 3728 1469 897 80 + 5 7941 5 0 2124 2025 1851 1861 80 + 5 7941 5 0 3069 3011 884 897 80 + 5 7941 5 0 2181 4601 898 181 80 + 5 7941 3 0 5977 1884 80 + 5 7941 5 0 3598 2121 819 1323 80 + 5 7941 4 0 3409 3555 897 80 + 5 7941 5 0 1517 4304 706 1334 80 + 5 7941 3 0 5681 2180 80 + 5 7941 4 0 5726 1954 181 80 + 5 7941 4 0 4120 2779 962 80 + 5 7941 4 0 5657 2023 181 80 + 5 7941 4 0 3934 2983 944 80 + 5 7941 5 0 3108 2855 1001 897 80 + 5 7941 5 0 3706 2832 1142 181 80 + 5 7941 4 0 3894 2784 1183 80 + 5 7941 5 0 1749 648 3608 1856 80 + 5 7941 4 0 2950 3130 884 977 + 5 7941 4 0 4062 1901 1898 80 + 5 7941 5 0 3072 1441 2451 897 80 + 5 7941 6 0 2695 3385 819 781 181 80 + 5 7941 4 0 1794 5170 897 80 + 5 7941 3 0 2423 5438 80 + 5 7941 5 0 3227 3555 898 181 80 + 5 7941 4 0 2218 4320 1323 80 + 5 7941 4 0 2948 4016 897 80 + 5 7941 5 0 2708 1356 2463 1334 80 + 5 7941 5 0 2318 3804 795 944 80 + 5 7941 5 0 2295 4232 372 962 80 + 5 7941 5 0 2580 3708 611 962 80 + 5 7941 4 0 3801 3098 962 80 + 5 7941 5 0 3808 3091 781 181 80 + 5 7941 6 0 2316 3806 643 915 181 80 + 5 7941 6 0 3384 2603 912 781 181 80 + 5 7941 4 0 2389 4289 1183 80 + 5 7941 2 0 7861 80 + 5 7941 5 0 2015 4474 1191 181 80 + 5 7941 3 0 3954 3010 977 + 5 7941 4 0 3229 2760 1872 80 + 5 7941 5 0 2317 1880 3483 181 80 + 5 7941 6 0 1842 2133 2924 781 181 80 + 5 7941 5 0 2697 2977 1243 944 80 + 5 7941 4 0 2326 4212 1323 80 + 5 7941 3 0 2810 3868 1263 + 5 7941 4 0 2869 3436 1556 80 + 5 7941 5 0 3266 2727 906 962 80 + 5 7941 4 0 3883 2655 1323 80 diff --git a/examples/06-NeutronCapture/HPGeTestStand.cc b/examples/06-NeutronCapture/HPGeTestStand.cc new file mode 100644 index 0000000..fb61020 --- /dev/null +++ b/examples/06-NeutronCapture/HPGeTestStand.cc @@ -0,0 +1,106 @@ +#include "HPGeTestStand.hh" + +#include "G4Box.hh" +#include "G4LogicalSkinSurface.hh" +#include "G4LogicalVolume.hh" +#include "G4Material.hh" +#include "G4NistManager.hh" +#include "G4OpticalSurface.hh" +#include "G4PVPlacement.hh" +#include "G4UnitsTable.hh" + +namespace u = CLHEP; + +G4VPhysicalVolume* HPGeTestStand::DefineGeometry() { + + G4State state; + std::string name, symbol; + std::vector elements; + std::vector mass_fraction; + double density; + double temperature; + double pressure; + double abundance; + int n_isotopes; + int n_components; + int n_atoms; + + auto man = G4NistManager::Instance(); + + // define enriched germanium + auto Ge70 = new G4Isotope(name = "Ge70", 32, 70, 69.92 * u::g / u::mole); + auto Ge72 = new G4Isotope(name = "Ge72", 32, 72, 71.92 * u::g / u::mole); + auto Ge73 = new G4Isotope(name = "Ge73", 32, 73, 73.00 * u::g / u::mole); + auto Ge74 = new G4Isotope(name = "Ge74", 32, 74, 74.00 * u::g / u::mole); + auto Ge76 = new G4Isotope(name = "Ge76", 32, 76, 76.00 * u::g / u::mole); + + auto el_enr_ge = new G4Element(name = "EnrichedGermanium", symbol = "EnrGe", n_isotopes = 5); + el_enr_ge->AddIsotope(Ge70, abundance = 0.0 * u::perCent); + el_enr_ge->AddIsotope(Ge72, abundance = 0.1 * u::perCent); + el_enr_ge->AddIsotope(Ge73, abundance = 0.2 * u::perCent); + el_enr_ge->AddIsotope(Ge74, abundance = 13.1 * u::perCent); + el_enr_ge->AddIsotope(Ge76, abundance = 86.6 * u::perCent); + + auto LAr = man->FindOrBuildMaterial("G4_lAr"); + auto water = man->FindOrBuildMaterial("G4_WATER"); + + density = 3.01 * u::g / u::cm3; + + G4Element* elGd = new G4Element("Gadolinium", "Gd", 64, 157.25 * u::g / u::mole); + G4Element* elS = new G4Element("Sulfur", "S", 16., 32.066 * u::g / u::mole); + G4Element* O = new G4Element("Oxygen", "O", 8., 16.00 * u::g / u::mole); + + G4Material* gadoliniumSulfate = new G4Material("GadoliniumSulfate", density, 3); // Gd2(SO4)3 + gadoliniumSulfate->AddElement(elGd, 2); + gadoliniumSulfate->AddElement(elS, 3); + gadoliniumSulfate->AddElement(O, 12); + + G4Material* Gdwater = new G4Material("GdLoadedWater", 1.000000 * u::g / u::cm3, 2); + Gdwater->AddMaterial(water, 1. - 0.002); + Gdwater->AddMaterial(gadoliniumSulfate, 0.002); + + + // Place Volumes + // ooooooooooooooooooooooOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOoooooooooooooooooooo + + auto mat_enr_ge = new G4Material("CryogenicEnrichedGermanium", + density = 5.56 * u::g / (u::cm * u::cm * u::cm), n_components = 1, state = G4State::kStateSolid, + temperature = LAr->GetTemperature(), pressure = LAr->GetPressure()); + + mat_enr_ge->AddElement(el_enr_ge, n_atoms = 1); + + auto world_s = new G4Box("WoldWater", 0.5 * u::m, 0.5 * u::m, 0.5 * u::m); + + auto world_l = new G4LogicalVolume(world_s, Gdwater, "WoldWater"); + + auto world_p = new G4PVPlacement(nullptr, G4ThreeVector(), world_l, "World", nullptr, false, 0); + + auto hpge_s = new G4Box("HPGe", 10 * u::cm, 10 * u::cm, 10 * u::cm); + + + auto hpge1_l = + new G4LogicalVolume(hpge_s, G4Material::GetMaterial("CryogenicEnrichedGermanium"), "HPGe1"); + auto hpge2_l = + new G4LogicalVolume(hpge_s, G4Material::GetMaterial("CryogenicEnrichedGermanium"), "HPGe2"); + auto hpge3_l = + new G4LogicalVolume(hpge_s, G4Material::GetMaterial("CryogenicEnrichedGermanium"), "HPGe3"); + auto hpge4_l = + new G4LogicalVolume(hpge_s, G4Material::GetMaterial("CryogenicEnrichedGermanium"), "HPGe4"); + + + // Place HPGe Detectors + auto spacing = 10.5 * u::cm; + new G4PVPlacement(nullptr, G4ThreeVector(+spacing, +spacing, 0), hpge1_l, "HPGe1", world_l, false, + 0); + new G4PVPlacement(nullptr, G4ThreeVector(+spacing, -spacing, 0), hpge2_l, "HPGe2", world_l, false, + 0); + new G4PVPlacement(nullptr, G4ThreeVector(-spacing, +spacing, 0), hpge3_l, "HPGe3", world_l, false, + 0); + new G4PVPlacement(nullptr, G4ThreeVector(-spacing, -spacing, 0), hpge4_l, "HPGe4", world_l, false, + 0); + + + return world_p; +} + +// vim: tabstop=2 textwidth=2 expandtab diff --git a/examples/06-NeutronCapture/HPGeTestStand.hh b/examples/06-NeutronCapture/HPGeTestStand.hh new file mode 100644 index 0000000..aa21780 --- /dev/null +++ b/examples/06-NeutronCapture/HPGeTestStand.hh @@ -0,0 +1,13 @@ +#ifndef _HPGE_TEST_STAND_HH_ +#define _HPGE_TEST_STAND_HH_ + +#include "RMGHardware.hh" + +class HPGeTestStand : public RMGHardware { + + public: + + G4VPhysicalVolume* DefineGeometry() override; +}; + +#endif diff --git a/examples/06-NeutronCapture/IsotopeOutputScheme.cc b/examples/06-NeutronCapture/IsotopeOutputScheme.cc new file mode 100644 index 0000000..cef87ed --- /dev/null +++ b/examples/06-NeutronCapture/IsotopeOutputScheme.cc @@ -0,0 +1,74 @@ + +#include "IsotopeOutputScheme.hh" + +#include "G4AnalysisManager.hh" +#include "G4Event.hh" +#include "G4EventManager.hh" + +#include "RMGLog.hh" +#include "RMGManager.hh" + +namespace u = CLHEP; + +IsotopeOutputScheme::IsotopeOutputScheme() {} + +void IsotopeOutputScheme::ClearBeforeEvent() { + zOfEvent.clear(); + aOfEvent.clear(); +} + +// invoked in RMGRunAction::SetupAnalysisManager() +void IsotopeOutputScheme::AssignOutputNames(G4AnalysisManager* ana_man) { + + auto rmg_man = RMGManager::Instance(); + + auto id = rmg_man->RegisterNtuple(OutputRegisterID, + ana_man->CreateNtuple("IsotopeOutput", "Event data")); + + ana_man->CreateNtupleIColumn(id, "evtid"); + // Could be done with particle PDG but this way it is human readable + ana_man->CreateNtupleIColumn(id, "Z"); + ana_man->CreateNtupleIColumn(id, "A"); + + ana_man->FinishNtuple(id); +} + +void IsotopeOutputScheme::TrackingActionPre(const G4Track* aTrack) { + const auto particle = aTrack->GetParticleDefinition(); + // if (!particle->IsGeneralIon()) return; + const int z = particle->GetAtomicNumber(); + const int a = particle->GetAtomicMass(); + if (z == 0 || a == 0) return; // not an isotope, but any other particle. + + G4String CreatorName = aTrack->GetCreatorProcess()->GetProcessName(); + G4int id = G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID(); + + // In Theory there should always be only one capture per event. But the more general approach is to be able to store multiple + if (CreatorName == "RMGnCapture" || CreatorName == "nCapture") { + zOfEvent.push_back(z); + aOfEvent.push_back(a); + } +} + +// invoked in RMGEventAction::EndOfEventAction() +void IsotopeOutputScheme::StoreEvent(const G4Event* event) { + auto rmg_man = RMGManager::Instance(); + if (rmg_man->IsPersistencyEnabled()) { + RMGLog::OutDev(RMGLog::debug, "Filling persistent data vectors"); + const auto ana_man = G4AnalysisManager::Instance(); + + for (size_t i = 0; i < zOfEvent.size(); ++i) { + auto ntupleid = rmg_man->GetNtupleID(OutputRegisterID); + + int col_id = 0; + ana_man->FillNtupleIColumn(ntupleid, col_id++, event->GetEventID()); + ana_man->FillNtupleIColumn(ntupleid, col_id++, zOfEvent[i]); + ana_man->FillNtupleIColumn(ntupleid, col_id++, aOfEvent[i]); + + // NOTE: must be called here for hit-oriented output + ana_man->AddNtupleRow(ntupleid); + } + } +} + +// vim: tabstop=2 shiftwidth=2 expandtab diff --git a/examples/06-NeutronCapture/IsotopeOutputScheme.hh b/examples/06-NeutronCapture/IsotopeOutputScheme.hh new file mode 100644 index 0000000..8772183 --- /dev/null +++ b/examples/06-NeutronCapture/IsotopeOutputScheme.hh @@ -0,0 +1,39 @@ +#ifndef _ISOTOPE_OUTPUT_SCHEME_HH_ +#define _ISOTOPE_OUTPUT_SCHEME_HH_ + +#include +#include + +#include "G4AnalysisManager.hh" +#include "G4GenericMessenger.hh" + +#include "RMGVOutputScheme.hh" + +class G4Event; +class IsotopeOutputScheme : public RMGVOutputScheme { + + public: + + IsotopeOutputScheme(); + + void ClearBeforeEvent() override; + void AssignOutputNames(G4AnalysisManager* ana_man) override; + void StoreEvent(const G4Event*) override; + // bool ShouldDiscardEvent(const G4Event*) override; + // std::optional StackingActionNewStage(int) override; + // std::optional StackingActionClassify(const G4Track*, int) override; + + void TrackingActionPre(const G4Track* aTrack) override; + + private: + + // Make sure the ID will not end up being used by other detectors (should be unique) + G4int OutputRegisterID = 1000; + + std::vector zOfEvent; + std::vector aOfEvent; +}; + +#endif + +// vim: tabstop=2 shiftwidth=2 expandtab diff --git a/examples/06-NeutronCapture/main.cc b/examples/06-NeutronCapture/main.cc new file mode 100644 index 0000000..e904405 --- /dev/null +++ b/examples/06-NeutronCapture/main.cc @@ -0,0 +1,32 @@ +#include "RMGLog.hh" +#include "RMGManager.hh" + +#include "HPGeTestStand.hh" +#include "IsotopeOutputScheme.hh" + +int main(int argc, char** argv) { + + // RMGLog::SetLogLevel(RMGLog::debug); + + RMGManager manager("06-NeutronCapture", argc, argv); + manager.SetUserInit(new HPGeTestStand()); + + auto user_init = manager.GetUserInit(); + user_init->AddOptionalOutputScheme("IsotopeOutputScheme"); + + // Need at least 1 active detector, otherwise persistency not enabled (even if manually setting it) + manager.GetDetectorConstruction()->RegisterDetector(RMGHardware::kGermanium, "HPGe1", 0); + // manager.GetDetectorConstruction()->RegisterDetector(RMGHardware::kGermanium, "HPGe2", 1); + // manager.GetDetectorConstruction()->RegisterDetector(RMGHardware::kGermanium, "HPGe3", 2); + // manager.GetDetectorConstruction()->RegisterDetector(RMGHardware::kGermanium, "HPGe4", 3); + + std::string macro = argc > 1 ? argv[1] : ""; + if (!macro.empty()) manager.IncludeMacroFile(macro); + else manager.SetInteractive(true); + + manager.SetNumberOfThreads(1); + manager.Initialize(); + manager.Run(); + + return 0; +} diff --git a/examples/06-NeutronCapture/runNormalnCapture.mac b/examples/06-NeutronCapture/runNormalnCapture.mac new file mode 100644 index 0000000..d6a482e --- /dev/null +++ b/examples/06-NeutronCapture/runNormalnCapture.mac @@ -0,0 +1,22 @@ +/RMG/Manager/Logging/LogLevel summary + +#/tracking/verbose 1 + +/RMG/Processes/HadronicPhysics Shielding +/RMG/Processes/ThermalScattering 1 + +/RMG/Output/ActivateOutputScheme IsotopeOutputScheme + +/run/initialize + +/RMG/Output/NtuplePerDetector false +/RMG/Output/FileName Normal.root + +/RMG/Generator/Select GPS +/gps/particle neutron +/gps/ang/type iso +/gps/energy 0.5 MeV + +#/RMG/Processes/Realm CosmicRays + +/run/beamOn 10000 diff --git a/examples/06-NeutronCapture/runRMGnCapture.mac b/examples/06-NeutronCapture/runRMGnCapture.mac new file mode 100644 index 0000000..bcfc3f3 --- /dev/null +++ b/examples/06-NeutronCapture/runRMGnCapture.mac @@ -0,0 +1,26 @@ +/RMG/Manager/Logging/LogLevel summary + +#/tracking/verbose 1 + +/RMG/Processes/HadronicPhysics Shielding +/RMG/Processes/ThermalScattering 1 +/RMG/Processes/UseGrabmayrsGammaCascades true + +/RMG/Output/ActivateOutputScheme IsotopeOutputScheme + +/RMG/GrabmayrGammaCascades/SetGammaCascadeFile 64 155 GammaCascades/156Gd-100Entries.txt +/RMG/GrabmayrGammaCascades/SetGammaCascadeFile 64 157 GammaCascades/158Gd-100Entries.txt +/RMG/GrabmayrGammaCascades/SetGammaCascadeRandomStartLocation 1 +/run/initialize + +/RMG/Output/NtuplePerDetector false +/RMG/Output/FileName Custom.root + +/RMG/Generator/Select GPS +/gps/particle neutron +/gps/ang/type iso +/gps/energy 0.5 MeV + +#/RMG/Processes/Realm CosmicRays + +/run/beamOn 10000 diff --git a/examples/06-NeutronCapture/vis-traj.mac b/examples/06-NeutronCapture/vis-traj.mac new file mode 100644 index 0000000..4eca000 --- /dev/null +++ b/examples/06-NeutronCapture/vis-traj.mac @@ -0,0 +1,27 @@ +#/RMG/Processes/OpticalPhysics true +/RMG/Processes/HadronicPhysics Shielding +/RMG/Processes/ThermalScattering 1 +/RMG/Output/ActivateOutputScheme IsotopeFilter +/run/initialize +/vis/open OGL + +/vis/scene/create +/vis/sceneHandler/attach + +/vis/drawVolume + +/vis/viewer/set/defaultColour black +/vis/viewer/set/background white + +/vis/scene/add/trajectories smooth +/vis/scene/add/hits +/vis/scene/endOfEventAction accumulate + +/RMG/Output/IsotopeFilter/AddIsotope 77 32 +/RMG/Output/IsotopeFilter/DiscardPhotonsIfIsotopeNotProduced true +/RMG/Output/NtuplePerDetector false + +/RMG/Generator/Select GPS +/gps/particle neutron +/gps/ang/type iso +/gps/energy 0.5 MeV diff --git a/include/RMGGrabmayrGCReader.hh b/include/RMGGrabmayrGCReader.hh new file mode 100644 index 0000000..77678a9 --- /dev/null +++ b/include/RMGGrabmayrGCReader.hh @@ -0,0 +1,88 @@ +// Copyright (C) 2022 Luigi Pertoldi +// +// This program is free software: you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// This program is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with this program. If not, see . + +#ifndef _RMG_GRABMAYR_GC_READER_HH_ +#define _RMG_GRABMAYR_GC_READER_HH_ + +#include +#include +#include +#include +#include + +#include "G4GenericMessenger.hh" +#include "globals.hh" + +// Modified from WLGDPetersGammaCascadeReader originally contributed by Moritz Neuberger +struct GammaCascadeLine { + G4int en; // neutron energy [keV] + G4int ex; // excitation energy [keV] + G4int m; // multiplicity of gamma cascade + G4int em; // missing energy [keV] + std::vector eg; // energies of photons [keV] +}; + + +class RMGGrabmayrGCReader { + public: + + static G4ThreadLocal RMGGrabmayrGCReader* GetInstance(); + ~RMGGrabmayrGCReader(); + // RMGGrabmayrGCReader& operator=(const RMGGrabmayrGCReader&) = delete; + + G4bool IsApplicable(G4int a, G4int z); + void CloseFiles(); + + GammaCascadeLine GetNextEntry(G4int z, G4int a); + + private: + + static G4ThreadLocal RMGGrabmayrGCReader* instance; + RMGGrabmayrGCReader(); + // std::vector> files; + // map holding the corresponding file for each isotope + std::map, std::unique_ptr> fCascadeFiles; + std::unique_ptr fGenericMessenger; + G4int fGammaCascadeRandomStartLocation = 0; + + void SetGammaCascadeFile(const G4int z, const G4int a, const G4String file_name); + void SetGammaCascadeRandomStartLocation(const int answer); + void SetStartLocation(std::ifstream& file); + + void RandomizeFiles(); + void DefineCommands(); + + // Nested messenger class + class GCMessenger : public G4UImessenger { + public: + + GCMessenger(RMGGrabmayrGCReader* reader); + ~GCMessenger(); + + void SetNewValue(G4UIcommand* command, G4String newValues) override; + + private: + + RMGGrabmayrGCReader* fReader; + G4UIcommand* fGammaFileCmd; + + void GammaFileCmd(const std::string& parameters); + }; + + std::unique_ptr fUIMessenger; +}; + + +#endif diff --git a/include/RMGNeutronCaptureProcess.hh b/include/RMGNeutronCaptureProcess.hh new file mode 100644 index 0000000..26ee062 --- /dev/null +++ b/include/RMGNeutronCaptureProcess.hh @@ -0,0 +1,35 @@ +// Copyright (C) 2022 Luigi Pertoldi +// +// This program is free software: you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// This program is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with this program. If not, see . + + +#ifndef _RMG_NEUTRON_CAPTURE_PROCESS_HH_ +#define _RMG_NEUTRON_CAPTURE_PROCESS_HH_ + +#include "G4HadronicProcess.hh" +#include "globals.hh" + +class RMGNeutronCaptureProcess : public G4HadronicProcess { + public: + + explicit RMGNeutronCaptureProcess(const G4String& processName = "RMGnCapture"); + + virtual ~RMGNeutronCaptureProcess(); + + G4bool IsApplicable(const G4ParticleDefinition& aParticleType) final; + + // Override PostStepDoIt from G4HadronicProcess + G4VParticleChange* PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) final; +}; +#endif diff --git a/include/RMGPhysics.hh b/include/RMGPhysics.hh index aa7f38f..308fd25 100644 --- a/include/RMGPhysics.hh +++ b/include/RMGPhysics.hh @@ -103,6 +103,7 @@ class RMGPhysics : public G4VModularPhysicsList { StepCutStore fStepCutsSensitive; bool fConstructOptical = false; bool fUseThermalScattering = false; + bool fUseGrabmayrGammaCascades = false; LowEnergyEMOption fLowEnergyEMOption = LowEnergyEMOption::kLivermore; HadronicPhysicsListOption fHadronicPhysicsListOption = HadronicPhysicsListOption::kNone; diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 981ce31..6ced21f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -13,10 +13,12 @@ set(PROJECT_PUBLIC_HEADERS ${_root}/include/RMGGeneratorUtil.hh ${_root}/include/RMGGermaniumDetector.hh ${_root}/include/RMGGermaniumOutputScheme.hh + ${_root}/include/RMGGrabmayrGCReader.hh ${_root}/include/RMGIsotopeFilterOutputScheme.hh ${_root}/include/RMGLog.hh ${_root}/include/RMGManager.hh ${_root}/include/RMGMasterGenerator.hh + ${_root}/include/RMGNeutronCaptureProcess.hh ${_root}/include/RMGNavigationTools.hh ${_root}/include/RMGOpticalDetector.hh ${_root}/include/RMGOpticalOutputScheme.hh @@ -49,11 +51,13 @@ set(PROJECT_SOURCES ${_root}/src/RMGGeneratorUtil.cc ${_root}/src/RMGGermaniumDetector.cc ${_root}/src/RMGGermaniumOutputScheme.cc + ${_root}/src/RMGGrabmayrGCReader.cc ${_root}/src/RMGIsotopeFilterOutputScheme.cc ${_root}/src/RMGLog.cc ${_root}/src/RMGManager.cc ${_root}/src/RMGMasterGenerator.cc ${_root}/src/RMGNavigationTools.cc + ${_root}/src/RMGNeutronCaptureProcess.cc ${_root}/src/RMGOpticalDetector.cc ${_root}/src/RMGOpticalOutputScheme.cc ${_root}/src/RMGPhysics.cc diff --git a/src/RMGGrabmayrGCReader.cc b/src/RMGGrabmayrGCReader.cc new file mode 100644 index 0000000..98cfe54 --- /dev/null +++ b/src/RMGGrabmayrGCReader.cc @@ -0,0 +1,180 @@ +#include "RMGGrabmayrGCReader.hh" + +#include "G4Tokenizer.hh" +#include "Randomize.hh" + +#include "RMGLog.hh" + + +G4ThreadLocal RMGGrabmayrGCReader* RMGGrabmayrGCReader::instance = nullptr; + +RMGGrabmayrGCReader* RMGGrabmayrGCReader::GetInstance() { + if (instance == nullptr) { instance = new RMGGrabmayrGCReader(); } + return instance; +} + +RMGGrabmayrGCReader::RMGGrabmayrGCReader() { DefineCommands(); } + +RMGGrabmayrGCReader::~RMGGrabmayrGCReader() { CloseFiles(); } + +void RMGGrabmayrGCReader::CloseFiles() { + RMGLog::Out(RMGLog::detail, "Closing gamma cascade files"); + for (const auto& el : fCascadeFiles) { + if (el.second->is_open()) { el.second->close(); } + } +} + +// Returns true if there exists a cascade file for the Isotope Z, A +G4bool RMGGrabmayrGCReader::IsApplicable(G4int z, G4int a) { + std::pair key = std::make_pair(z, a); + auto it = fCascadeFiles.find(key); + return it != fCascadeFiles.end(); +} + +// Returns the next cascade for the Isotope Z, A. +GammaCascadeLine RMGGrabmayrGCReader::GetNextEntry(G4int z, G4int a) { + // Find the corresponding file + std::pair key = std::make_pair(z, a); + auto it = fCascadeFiles.find(key); + if (it == fCascadeFiles.end()) + RMGLog::OutFormat(RMGLog::fatal, "Isotope Z: {} A: {} does not exist.", z, a); + // read next line from file + std::string line; + do { + if (!std::getline(*(it->second), line)) { + // if end-of-file is reached, reset file and read first line + RMGLog::Out(RMGLog::detail, "Gamma cascade file EOF reached, re-opening the file"); + it->second->clear(); // clear EOF flag + it->second->seekg(0, std::ios::beg); // move to beginning of file + if (!std::getline(*(it->second), line)) { + RMGLog::Out(RMGLog::fatal, "Failed to read next line after re-opening the file. Exit!"); + } + } + } while (line[0] == '%' || (line.find("version") != + std::string::npos)); // This could be outsourced to SetStartLocation + + // parse line and return as struct + GammaCascadeLine gamma_cascade; + std::istringstream iss(line); + iss >> gamma_cascade.en >> gamma_cascade.ex >> gamma_cascade.m >> gamma_cascade.em; + gamma_cascade.eg.reserve(gamma_cascade.m); + int eg_value; + for (int i = 0; i < gamma_cascade.m; i++) { + if (!(iss >> eg_value)) { + RMGLog::Out(RMGLog::fatal, "Failed to read gamma energy from file. Exit!"); + } + gamma_cascade.eg.push_back(eg_value); + } + return gamma_cascade; +} + + +void RMGGrabmayrGCReader::SetStartLocation(std::ifstream& file) { + if (!file.is_open()) + RMGLog::Out(RMGLog::fatal, "The file is not open to set start location! Exit."); + file.clear(); // clear EOF flag + file.seekg(0, std::ios::beg); // move to beginning of file + // Skip Header + std::string line; + int header_length = 0; + do { + std::getline(file, line); + header_length++; + } while (line[0] == '%' || (line.find("version") != std::string::npos)); + + // In case the Random start location macro is set + if (fGammaCascadeRandomStartLocation) { + int n_entries_in_file = 0; + // Seems excessiv to read through the entire file, there has to be a quicker way? + while (std::getline(file, line)) n_entries_in_file++; + + file.clear(); // clear EOF flag + file.seekg(0, std::ios::beg); // move to beginning of file + + int start_location = (int)(n_entries_in_file * G4UniformRand() + header_length); + + RMGLog::Out(RMGLog::detail, "Random start location: ", start_location); + for (int j = 0; j < start_location; j++) std::getline(file, line); + } +} + +void RMGGrabmayrGCReader::SetGammaCascadeFile(const G4int z, const G4int a, const G4String file_name) { + + RMGLog::Out(RMGLog::detail, "Opening file ", file_name); + std::unique_ptr file = std::make_unique(file_name); + + if (z == 0 || a == 0) + RMGLog::OutFormat(RMGLog::fatal, "Isotope Z: {} A: {} does not exist.", z, a); + if (!file || !file->is_open()) + RMGLog::Out(RMGLog::fatal, "Gamma cascade file: " + file_name + " not found! Exit."); + + SetStartLocation(*file); + + fCascadeFiles.insert({{z, a}, std::move(file)}); +} + +void RMGGrabmayrGCReader::RandomizeFiles() { + RMGLog::Out(RMGLog::detail, "(Un)-Randomizing start locations"); + for (auto& el : fCascadeFiles) { SetStartLocation(*el.second); } +} + +void RMGGrabmayrGCReader::SetGammaCascadeRandomStartLocation(const int answer) { + fGammaCascadeRandomStartLocation = answer; + RMGLog::Out(RMGLog::detail, + "setting fGammaCascadeRandomStartLocation to: ", fGammaCascadeRandomStartLocation); + RandomizeFiles(); +} + +void RMGGrabmayrGCReader::DefineCommands() { + fGenericMessenger = std::make_unique(this, "/RMG/GrabmayrGammaCascades/", + "Control Peters gamma cascade model"); + + fGenericMessenger + ->DeclareMethod("SetGammaCascadeRandomStartLocation", + &RMGGrabmayrGCReader::SetGammaCascadeRandomStartLocation) + .SetGuidance("Set the whether the start location in the gamma cascade file is random or not") + .SetGuidance("0 = don't") + .SetGuidance("1 = do") + .SetCandidates("0 1") + .SetDefaultValue("0") + .SetStates(G4State_PreInit, G4State_Idle); + + // SetGammaCascadeFile cannot be defined with the G4GenericMessenger (it has too many parameters). + fUIMessenger = std::make_unique(this); +} + + +RMGGrabmayrGCReader::GCMessenger::GCMessenger(RMGGrabmayrGCReader* reader) : fReader(reader) { + fGammaFileCmd = new G4UIcommand("/RMG/GrabmayrGammaCascades/SetGammaCascadeFile", this); + fGammaFileCmd->SetGuidance("Set a gamma cascade file for neutron capture on a specified isotope"); + + auto p_Z = new G4UIparameter("Z", 'i', false); + p_Z->SetGuidance("Z of isotope"); + fGammaFileCmd->SetParameter(p_Z); + + auto p_A = new G4UIparameter("A", 'i', false); + p_A->SetGuidance("A of isotope"); + fGammaFileCmd->SetParameter(p_A); + + auto p_file = new G4UIparameter("file", 's', false); + p_file->SetGuidance("/path/to/file of gamma cascade"); + fGammaFileCmd->SetParameter(p_file); + + fGammaFileCmd->AvailableForStates(G4State_PreInit, G4State_Idle); +} + +RMGGrabmayrGCReader::GCMessenger::~GCMessenger() { delete fGammaFileCmd; } + +void RMGGrabmayrGCReader::GCMessenger::SetNewValue(G4UIcommand* command, G4String newValues) { + if (command == fGammaFileCmd) GammaFileCmd(newValues); +} + +void RMGGrabmayrGCReader::GCMessenger::GammaFileCmd(const std::string& parameters) { + G4Tokenizer next(parameters); + + auto Z = std::stoi(next()); + auto A = std::stoi(next()); + auto file = next(); + + fReader->SetGammaCascadeFile(Z, A, file); +} diff --git a/src/RMGNeutronCaptureProcess.cc b/src/RMGNeutronCaptureProcess.cc new file mode 100644 index 0000000..ee9a34a --- /dev/null +++ b/src/RMGNeutronCaptureProcess.cc @@ -0,0 +1,124 @@ +// Copyright (C) 2022 Luigi Pertoldi +// +// This program is free software: you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// This program is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with this program. If not, see . + + +#include "RMGNeutronCaptureProcess.hh" + +#include "G4Element.hh" +#include "G4Gamma.hh" +#include "G4IonTable.hh" +#include "G4Isotope.hh" +#include "G4Neutron.hh" +#include "G4NeutronCaptureXS.hh" +#include "G4NistManager.hh" +#include "G4Nucleus.hh" +#include "G4ParticleChange.hh" +#include "G4ParticleDefinition.hh" +#include "G4RandomDirection.hh" +#include "G4Step.hh" +#include "G4Track.hh" + +#include "RMGGrabmayrGCReader.hh" +#include "RMGLog.hh" + +namespace u = CLHEP; + +RMGNeutronCaptureProcess::RMGNeutronCaptureProcess(const G4String& processName) + : G4HadronicProcess(processName, fCapture) { + AddDataSet(new G4NeutronCaptureXS()); +} + +RMGNeutronCaptureProcess::~RMGNeutronCaptureProcess() {} + +G4bool RMGNeutronCaptureProcess::IsApplicable(const G4ParticleDefinition& aParticleType) { + return (&aParticleType == G4Neutron::Neutron()); +} + +// Do exactly as in G4HadronicProcess, unless the capture happens on certain isotopes +G4VParticleChange* RMGNeutronCaptureProcess::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) { + // Get the proposed result from the "normal" process + G4VParticleChange* ProposedResult = G4HadronicProcess::PostStepDoIt(aTrack, aStep); + + RMGGrabmayrGCReader* CascadeReader = RMGGrabmayrGCReader::GetInstance(); + + G4int nSec = ProposedResult->GetNumberOfSecondaries(); + G4int z; + G4int a; + G4bool IsApplicable = false; + // Search through the proposed secondaries for our Isotopes of interest + if (nSec > 0) { + for (G4int i = 0; i < nSec; ++i) { + const auto particle = ProposedResult->GetSecondary(i)->GetParticleDefinition(); + z = particle->GetAtomicNumber(); + a = particle->GetAtomicMass() - + 1; // -1 because this is after the neutron capture, but the functions expect the target nucleus + if (z == 0 || a == 0) continue; // not an isotope, but any other particle. Quick skip + if (CascadeReader->IsApplicable(z, a)) { + IsApplicable = true; + // As only one Isotope is possible per Capture, break here so z and a remains correct. + break; + } + } + } + // If not Applicable we leave everything untouched + if (!IsApplicable) return ProposedResult; + + // If Applicable we have to propose our own result + RMGLog::OutDev(RMGLog::debug, "Proposing own neutron capture result"); + // If we do it ourselves, we are responsible that everything is done correctly. + // Most of this is copied from G4HadronicProcess PostStepDoIt() or FillResult(), as we are replacing that method. + theTotalResult->Clear(); + theTotalResult->Initialize(aTrack); + fWeight = aTrack.GetWeight(); + theTotalResult->ProposeWeight(fWeight); + + // Can skip all of the Checks for illegal track status, as in that case no Applicable isotope would have been produced! + + // Do our Physics. + GammaCascadeLine input = CascadeReader->GetNextEntry(z, a); + G4ThreeVector location = aTrack.GetPosition(); + G4double time = aTrack.GetGlobalTime(); + + // Expect no deposition as the total energy will be distributed + theTotalResult->ProposeLocalEnergyDeposit(0.0); + theTotalResult->ProposeTrackStatus(fStopAndKill); + theTotalResult->ProposeEnergy(0.0); + + // Produce the created Isotope + the Gammas + theTotalResult->SetNumberOfSecondaries(input.m + 1); + G4IonTable* theTable = G4IonTable::GetIonTable(); + + G4ParticleDefinition* particleDef_nuc = theTable->GetIon(z, a + 1, (double)(input.em * u::keV)); + G4DynamicParticle* particle_nuc = new G4DynamicParticle(particleDef_nuc, G4RandomDirection()); + G4Track* secondary_nuc = new G4Track(particle_nuc, time, location); + + // secondary_nuc->SetCreatorModelID(idModel); No idea what this should be. Not relevant? + secondary_nuc->SetWeight(fWeight); + secondary_nuc->SetTouchableHandle(aTrack.GetTouchableHandle()); + theTotalResult->AddSecondary(secondary_nuc); + + for (const G4int energy : input.eg) { + G4ParticleDefinition* particleDef_gamma = G4Gamma::Gamma(); + G4DynamicParticle* particle_gamma = + new G4DynamicParticle(particleDef_gamma, energy * u::keV, G4RandomDirection()); + G4Track* secondary_gamma = new G4Track(particle_gamma, time, location); + secondary_gamma->SetWeight(fWeight); + secondary_gamma->SetTouchableHandle(aTrack.GetTouchableHandle()); + secondary_gamma->SetKineticEnergy(energy * u::keV); + theTotalResult->AddSecondary(secondary_gamma); + } + + return theTotalResult; +} diff --git a/src/RMGPhysics.cc b/src/RMGPhysics.cc index 0818b54..8fcc79c 100644 --- a/src/RMGPhysics.cc +++ b/src/RMGPhysics.cc @@ -42,6 +42,7 @@ #include "G4IonTable.hh" #include "G4LeptonConstructor.hh" #include "G4MesonConstructor.hh" +#include "G4NeutronCaptureProcess.hh" #include "G4NuclearLevelData.hh" #include "G4OpAbsorption.hh" #include "G4OpBoundaryProcess.hh" @@ -49,6 +50,7 @@ #include "G4OpWLS.hh" #include "G4OpticalParameters.hh" #include "G4ParticleDefinition.hh" +#include "G4ParticleHPCaptureData.hh" #include "G4ParticleHPElastic.hh" #include "G4ParticleHPElasticData.hh" #include "G4ParticleHPThermalScattering.hh" @@ -64,6 +66,7 @@ #include "G4ThermalNeutrons.hh" #include "RMGLog.hh" +#include "RMGNeutronCaptureProcess.hh" #include "RMGTools.hh" namespace u = CLHEP; @@ -222,6 +225,33 @@ void RMGPhysics::ConstructProcess() { RMGLog::Out(RMGLog::detail, "Adding hadronic inelastic physics"); hPhysics->ConstructProcess(); + if (fUseGrabmayrGammaCascades) { + // Apply RMG custom neutron capture + // Mostly similar to examples/extended/Hadr04 + auto pManager = G4Neutron::Neutron()->GetProcessManager(); + auto processVector = pManager->GetProcessList(); + // Find the existing neutron capture process + auto neutronCaptureProcess = + dynamic_cast(pManager->GetProcess("nCapture")); + + // Overwrite the old Process, keeping all of the interactions + if (neutronCaptureProcess) { + RMGLog::Out(RMGLog::detail, "Overwriting NeutronCaptureProcess"); + pManager->RemoveProcess(neutronCaptureProcess); + auto RMGProcess = new RMGNeutronCaptureProcess(); + // HP cross section data set not naturally in G4NeutronCaptureProcess + auto dataSet = new G4ParticleHPCaptureData(); + RMGProcess->AddDataSet(dataSet); + // Move all interactions to the new process + for (auto& el : neutronCaptureProcess->GetHadronicInteractionList()) { + RMGProcess->RegisterMe(el); + } + pManager->AddDiscreteProcess(RMGProcess); + } else { + RMGLog::Out(RMGLog::error, "Could not apply custom neutron capture model"); + } + } + RMGLog::Out(RMGLog::detail, "Adding stopping physics"); G4VPhysicsConstructor* stoppingPhysics = new G4StoppingPhysics(G4VModularPhysicsList::verboseLevel); @@ -466,6 +496,10 @@ void RMGPhysics::DefineCommands() { .SetGuidance("Store e- internal conversion data") .SetCandidates("0 1") .SetStates(G4State_PreInit); + + fMessenger->DeclareProperty("UseGrabmayrsGammaCascades", fUseGrabmayrGammaCascades) + .SetGuidance("Use custom RMGNeutronCapture to apply Grabmayrs gamma cascades.") + .SetStates(G4State_PreInit); } // vim: shiftwidth=2 tabstop=2 expandtab diff --git a/src/RMGUserAction.cc b/src/RMGUserAction.cc index bb9ede8..4a5013d 100644 --- a/src/RMGUserAction.cc +++ b/src/RMGUserAction.cc @@ -21,6 +21,7 @@ #include "G4MultiTrackingAction.hh" #include "RMGEventAction.hh" +#include "RMGGrabmayrGCReader.hh" #include "RMGManager.hh" #include "RMGMasterGenerator.hh" #include "RMGRunAction.hh" @@ -35,6 +36,7 @@ void RMGUserAction::BuildForMaster() const { auto generator_primary = new RMGMasterGenerator(); this->SetUserAction( new RMGRunAction(generator_primary, RMGManager::Instance()->IsPersistencyEnabled())); + RMGGrabmayrGCReader::GetInstance(); } void RMGUserAction::Build() const { @@ -70,6 +72,7 @@ void RMGUserAction::Build() const { this->SetUserAction(new RMGStackingAction(run_action)); this->SetUserAction(stepping_action); this->SetUserAction(tracking_action); + RMGGrabmayrGCReader::GetInstance(); } // vim: tabstop=2 shiftwidth=2 expandtab