diff --git a/config/AR23_20i/CommonParam.xml b/config/AR23_20i/CommonParam.xml index e0a74f912f..6ad7145a75 100644 --- a/config/AR23_20i/CommonParam.xml +++ b/config/AR23_20i/CommonParam.xml @@ -34,6 +34,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/EX00_00a/CommonParam.xml b/config/EX00_00a/CommonParam.xml index 2cff53c44f..5714e92a45 100644 --- a/config/EX00_00a/CommonParam.xml +++ b/config/EX00_00a/CommonParam.xml @@ -34,6 +34,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_01a/G18_01a_02_11a/CommonParam.xml b/config/G18_01a/G18_01a_02_11a/CommonParam.xml index c2ecd246b3..a257fb085e 100644 --- a/config/G18_01a/G18_01a_02_11a/CommonParam.xml +++ b/config/G18_01a/G18_01a_02_11a/CommonParam.xml @@ -42,7 +42,9 @@ University of Liverpool - 0.501716712132 + 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_01a/G18_01a_02_11b/CommonParam.xml b/config/G18_01a/G18_01a_02_11b/CommonParam.xml index 210900c27a..5c6a3c940e 100644 --- a/config/G18_01a/G18_01a_02_11b/CommonParam.xml +++ b/config/G18_01a/G18_01a_02_11b/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_01b/G18_01b_02_11a/CommonParam.xml b/config/G18_01b/G18_01b_02_11a/CommonParam.xml index 669767b91f..204e5f8673 100644 --- a/config/G18_01b/G18_01b_02_11a/CommonParam.xml +++ b/config/G18_01b/G18_01b_02_11a/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_01b/G18_01b_02_11b/CommonParam.xml b/config/G18_01b/G18_01b_02_11b/CommonParam.xml index ed226bb787..e68d7dec5a 100644 --- a/config/G18_01b/G18_01b_02_11b/CommonParam.xml +++ b/config/G18_01b/G18_01b_02_11b/CommonParam.xml @@ -43,6 +43,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_02a/G18_02a_02_11a/CommonParam.xml b/config/G18_02a/G18_02a_02_11a/CommonParam.xml index 7ebf0768e4..375b1dd2f7 100644 --- a/config/G18_02a/G18_02a_02_11a/CommonParam.xml +++ b/config/G18_02a/G18_02a_02_11a/CommonParam.xml @@ -41,7 +41,9 @@ University of Liverpool - 0.501716712132 + 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_02a/G18_02a_02_11b/CommonParam.xml b/config/G18_02a/G18_02a_02_11b/CommonParam.xml index 3348d5ac36..486fc8c9e1 100644 --- a/config/G18_02a/G18_02a_02_11b/CommonParam.xml +++ b/config/G18_02a/G18_02a_02_11b/CommonParam.xml @@ -48,6 +48,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_02a/G18_02a_03_320/CommonParam.xml b/config/G18_02a/G18_02a_03_320/CommonParam.xml index 1681a084ac..b34ad7eb9f 100644 --- a/config/G18_02a/G18_02a_03_320/CommonParam.xml +++ b/config/G18_02a/G18_02a_03_320/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_02a/G18_02a_03_330/CommonParam.xml b/config/G18_02a/G18_02a_03_330/CommonParam.xml index 1681a084ac..b34ad7eb9f 100644 --- a/config/G18_02a/G18_02a_03_330/CommonParam.xml +++ b/config/G18_02a/G18_02a_03_330/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_02b/G18_02b_02_11a/CommonParam.xml b/config/G18_02b/G18_02b_02_11a/CommonParam.xml index 2aa6870caf..6f2450f5c2 100644 --- a/config/G18_02b/G18_02b_02_11a/CommonParam.xml +++ b/config/G18_02b/G18_02b_02_11a/CommonParam.xml @@ -42,7 +42,9 @@ University of Liverpool - 0.501716712132 + 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_02b/G18_02b_02_11b/CommonParam.xml b/config/G18_02b/G18_02b_02_11b/CommonParam.xml index 1681a084ac..b34ad7eb9f 100644 --- a/config/G18_02b/G18_02b_02_11b/CommonParam.xml +++ b/config/G18_02b/G18_02b_02_11b/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_02c/G18_02c_02_11a/CommonParam.xml b/config/G18_02c/G18_02c_02_11a/CommonParam.xml index 7ebf0768e4..375b1dd2f7 100644 --- a/config/G18_02c/G18_02c_02_11a/CommonParam.xml +++ b/config/G18_02c/G18_02c_02_11a/CommonParam.xml @@ -41,7 +41,9 @@ University of Liverpool - 0.501716712132 + 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_02c/G18_02c_02_11b/CommonParam.xml b/config/G18_02c/G18_02c_02_11b/CommonParam.xml index 1681a084ac..b34ad7eb9f 100644 --- a/config/G18_02c/G18_02c_02_11b/CommonParam.xml +++ b/config/G18_02c/G18_02c_02_11b/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_02c/G18_02c_03_320/CommonParam.xml b/config/G18_02c/G18_02c_03_320/CommonParam.xml index 1681a084ac..b34ad7eb9f 100644 --- a/config/G18_02c/G18_02c_03_320/CommonParam.xml +++ b/config/G18_02c/G18_02c_03_320/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_02c/G18_02c_03_330/CommonParam.xml b/config/G18_02c/G18_02c_03_330/CommonParam.xml index 1681a084ac..b34ad7eb9f 100644 --- a/config/G18_02c/G18_02c_03_330/CommonParam.xml +++ b/config/G18_02c/G18_02c_03_330/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_02d/G18_02d_02_11a/CommonParam.xml b/config/G18_02d/G18_02d_02_11a/CommonParam.xml index 7ebf0768e4..375b1dd2f7 100644 --- a/config/G18_02d/G18_02d_02_11a/CommonParam.xml +++ b/config/G18_02d/G18_02d_02_11a/CommonParam.xml @@ -41,7 +41,9 @@ University of Liverpool - 0.501716712132 + 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_02d/G18_02d_02_11b/CommonParam.xml b/config/G18_02d/G18_02d_02_11b/CommonParam.xml index 1681a084ac..b34ad7eb9f 100644 --- a/config/G18_02d/G18_02d_02_11b/CommonParam.xml +++ b/config/G18_02d/G18_02d_02_11b/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_02d/G18_02d_03_320/CommonParam.xml b/config/G18_02d/G18_02d_03_320/CommonParam.xml index 1681a084ac..b34ad7eb9f 100644 --- a/config/G18_02d/G18_02d_03_320/CommonParam.xml +++ b/config/G18_02d/G18_02d_03_320/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_02d/G18_02d_03_330/CommonParam.xml b/config/G18_02d/G18_02d_03_330/CommonParam.xml index 1681a084ac..b34ad7eb9f 100644 --- a/config/G18_02d/G18_02d_03_330/CommonParam.xml +++ b/config/G18_02d/G18_02d_03_330/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_10a/G18_10a_02_11a/CommonParam.xml b/config/G18_10a/G18_10a_02_11a/CommonParam.xml index bdd16dbd44..1f48bd7cb1 100644 --- a/config/G18_10a/G18_10a_02_11a/CommonParam.xml +++ b/config/G18_10a/G18_10a_02_11a/CommonParam.xml @@ -42,7 +42,9 @@ University of Liverpool - 0.501716712132 + 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_10a/G18_10a_02_11b/CommonParam.xml b/config/G18_10a/G18_10a_02_11b/CommonParam.xml index 1681a084ac..b34ad7eb9f 100644 --- a/config/G18_10a/G18_10a_02_11b/CommonParam.xml +++ b/config/G18_10a/G18_10a_02_11b/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_10a/G18_10a_03_320/CommonParam.xml b/config/G18_10a/G18_10a_03_320/CommonParam.xml index 1681a084ac..b34ad7eb9f 100644 --- a/config/G18_10a/G18_10a_03_320/CommonParam.xml +++ b/config/G18_10a/G18_10a_03_320/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_10a/G18_10a_03_330/CommonParam.xml b/config/G18_10a/G18_10a_03_330/CommonParam.xml index 1681a084ac..b34ad7eb9f 100644 --- a/config/G18_10a/G18_10a_03_330/CommonParam.xml +++ b/config/G18_10a/G18_10a_03_330/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_10b/G18_10b_02_11a/CommonParam.xml b/config/G18_10b/G18_10b_02_11a/CommonParam.xml index ccd795ec71..74885f348c 100644 --- a/config/G18_10b/G18_10b_02_11a/CommonParam.xml +++ b/config/G18_10b/G18_10b_02_11a/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_10b/G18_10b_02_11b/CommonParam.xml b/config/G18_10b/G18_10b_02_11b/CommonParam.xml index 1681a084ac..b34ad7eb9f 100644 --- a/config/G18_10b/G18_10b_02_11b/CommonParam.xml +++ b/config/G18_10b/G18_10b_02_11b/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_10c/G18_10c_02_11a/CommonParam.xml b/config/G18_10c/G18_10c_02_11a/CommonParam.xml index bdd16dbd44..1f48bd7cb1 100644 --- a/config/G18_10c/G18_10c_02_11a/CommonParam.xml +++ b/config/G18_10c/G18_10c_02_11a/CommonParam.xml @@ -42,7 +42,9 @@ University of Liverpool - 0.501716712132 + 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_10c/G18_10c_02_11b/CommonParam.xml b/config/G18_10c/G18_10c_02_11b/CommonParam.xml index 1681a084ac..b34ad7eb9f 100644 --- a/config/G18_10c/G18_10c_02_11b/CommonParam.xml +++ b/config/G18_10c/G18_10c_02_11b/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_10c/G18_10c_03_320/CommonParam.xml b/config/G18_10c/G18_10c_03_320/CommonParam.xml index 1681a084ac..b34ad7eb9f 100644 --- a/config/G18_10c/G18_10c_03_320/CommonParam.xml +++ b/config/G18_10c/G18_10c_03_320/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_10c/G18_10c_03_330/CommonParam.xml b/config/G18_10c/G18_10c_03_330/CommonParam.xml index 1681a084ac..b34ad7eb9f 100644 --- a/config/G18_10c/G18_10c_03_330/CommonParam.xml +++ b/config/G18_10c/G18_10c_03_330/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_10d/G18_10d_02_11a/CommonParam.xml b/config/G18_10d/G18_10d_02_11a/CommonParam.xml index bdd16dbd44..1f48bd7cb1 100644 --- a/config/G18_10d/G18_10d_02_11a/CommonParam.xml +++ b/config/G18_10d/G18_10d_02_11a/CommonParam.xml @@ -42,7 +42,9 @@ University of Liverpool - 0.501716712132 + 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_10d/G18_10d_02_11b/CommonParam.xml b/config/G18_10d/G18_10d_02_11b/CommonParam.xml index 1681a084ac..b34ad7eb9f 100644 --- a/config/G18_10d/G18_10d_02_11b/CommonParam.xml +++ b/config/G18_10d/G18_10d_02_11b/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_10d/G18_10d_03_320/CommonParam.xml b/config/G18_10d/G18_10d_03_320/CommonParam.xml index 1681a084ac..b34ad7eb9f 100644 --- a/config/G18_10d/G18_10d_03_320/CommonParam.xml +++ b/config/G18_10d/G18_10d_03_320/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_10d/G18_10d_03_330/CommonParam.xml b/config/G18_10d/G18_10d_03_330/CommonParam.xml index 3a145e73b8..35991ff996 100644 --- a/config/G18_10d/G18_10d_03_330/CommonParam.xml +++ b/config/G18_10d/G18_10d_03_330/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_10i/G18_10i_02_11a/CommonParam.xml b/config/G18_10i/G18_10i_02_11a/CommonParam.xml index bdd16dbd44..1f48bd7cb1 100644 --- a/config/G18_10i/G18_10i_02_11a/CommonParam.xml +++ b/config/G18_10i/G18_10i_02_11a/CommonParam.xml @@ -42,7 +42,9 @@ University of Liverpool - 0.501716712132 + 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_10i/G18_10i_02_11b/CommonParam.xml b/config/G18_10i/G18_10i_02_11b/CommonParam.xml index 1681a084ac..b34ad7eb9f 100644 --- a/config/G18_10i/G18_10i_02_11b/CommonParam.xml +++ b/config/G18_10i/G18_10i_02_11b/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_10i/G18_10i_03_320/CommonParam.xml b/config/G18_10i/G18_10i_03_320/CommonParam.xml index 1681a084ac..b34ad7eb9f 100644 --- a/config/G18_10i/G18_10i_03_320/CommonParam.xml +++ b/config/G18_10i/G18_10i_03_320/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_10i/G18_10i_03_330/CommonParam.xml b/config/G18_10i/G18_10i_03_330/CommonParam.xml index 1681a084ac..b34ad7eb9f 100644 --- a/config/G18_10i/G18_10i_03_330/CommonParam.xml +++ b/config/G18_10i/G18_10i_03_330/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_10j/G18_10j_02_11a/CommonParam.xml b/config/G18_10j/G18_10j_02_11a/CommonParam.xml index bdd16dbd44..1f48bd7cb1 100644 --- a/config/G18_10j/G18_10j_02_11a/CommonParam.xml +++ b/config/G18_10j/G18_10j_02_11a/CommonParam.xml @@ -42,7 +42,9 @@ University of Liverpool - 0.501716712132 + 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_10j/G18_10j_02_11b/CommonParam.xml b/config/G18_10j/G18_10j_02_11b/CommonParam.xml index 1681a084ac..b34ad7eb9f 100644 --- a/config/G18_10j/G18_10j_02_11b/CommonParam.xml +++ b/config/G18_10j/G18_10j_02_11b/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_10j/G18_10j_03_320/CommonParam.xml b/config/G18_10j/G18_10j_03_320/CommonParam.xml index 1681a084ac..b34ad7eb9f 100644 --- a/config/G18_10j/G18_10j_03_320/CommonParam.xml +++ b/config/G18_10j/G18_10j_03_320/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G18_10j/G18_10j_03_330/CommonParam.xml b/config/G18_10j/G18_10j_03_330/CommonParam.xml index 1681a084ac..b34ad7eb9f 100644 --- a/config/G18_10j/G18_10j_03_330/CommonParam.xml +++ b/config/G18_10j/G18_10j_03_330/CommonParam.xml @@ -42,6 +42,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G21_11a/CommonParam.xml b/config/G21_11a/CommonParam.xml index 27b0717f3b..e1d93ea807 100644 --- a/config/G21_11a/CommonParam.xml +++ b/config/G21_11a/CommonParam.xml @@ -35,7 +35,9 @@ University of Liverpool - 0.501716712132 + 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G21_11b/CommonParam.xml b/config/G21_11b/CommonParam.xml index bfbee303c2..3ba12eaa55 100644 --- a/config/G21_11b/CommonParam.xml +++ b/config/G21_11b/CommonParam.xml @@ -35,7 +35,9 @@ University of Liverpool - 0.501716712132 + 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G21_11c/CommonParam.xml b/config/G21_11c/CommonParam.xml index 27b0717f3b..e1d93ea807 100644 --- a/config/G21_11c/CommonParam.xml +++ b/config/G21_11c/CommonParam.xml @@ -35,7 +35,9 @@ University of Liverpool - 0.501716712132 + 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/G21_11d/CommonParam.xml b/config/G21_11d/CommonParam.xml index 27b0717f3b..e1d93ea807 100644 --- a/config/G21_11d/CommonParam.xml +++ b/config/G21_11d/CommonParam.xml @@ -35,7 +35,9 @@ University of Liverpool - 0.501716712132 + 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/GEM21_11a/CommonParam.xml b/config/GEM21_11a/CommonParam.xml index fa2629d926..62e3c7d3c7 100644 --- a/config/GEM21_11a/CommonParam.xml +++ b/config/GEM21_11a/CommonParam.xml @@ -35,7 +35,9 @@ University of Liverpool - 0.501716712132 + 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/GEM21_11b/CommonParam.xml b/config/GEM21_11b/CommonParam.xml index 423ba90290..0d3a999a17 100644 --- a/config/GEM21_11b/CommonParam.xml +++ b/config/GEM21_11b/CommonParam.xml @@ -35,7 +35,9 @@ University of Liverpool - 0.501716712132 + 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/GEM21_11c/CommonParam.xml b/config/GEM21_11c/CommonParam.xml index f551fcf1f6..15f7c64b73 100644 --- a/config/GEM21_11c/CommonParam.xml +++ b/config/GEM21_11c/CommonParam.xml @@ -35,7 +35,9 @@ University of Liverpool - 0.501716712132 + 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/GEM21_11d/CommonParam.xml b/config/GEM21_11d/CommonParam.xml index f551fcf1f6..15f7c64b73 100644 --- a/config/GEM21_11d/CommonParam.xml +++ b/config/GEM21_11d/CommonParam.xml @@ -35,7 +35,9 @@ University of Liverpool - 0.501716712132 + 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/GPRD18_10a/GPRD18_10a_02_11b/CommonParam.xml b/config/GPRD18_10a/GPRD18_10a_02_11b/CommonParam.xml index fb629fbb74..85cf5fc3c7 100644 --- a/config/GPRD18_10a/GPRD18_10a_02_11b/CommonParam.xml +++ b/config/GPRD18_10a/GPRD18_10a_02_11b/CommonParam.xml @@ -37,6 +37,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/GVLE18_01a/CommonParam.xml b/config/GVLE18_01a/CommonParam.xml index 3536a9147b..56ff0a910f 100644 --- a/config/GVLE18_01a/CommonParam.xml +++ b/config/GVLE18_01a/CommonParam.xml @@ -34,6 +34,8 @@ University of Liverpool 0.501716712132 + 1.16639E-5 + 7.2973525332858855E-3 diff --git a/config/HG4BertCascIntranuke.xml b/config/HG4BertCascIntranuke.xml index 237d78e3c1..1e59072a07 100644 --- a/config/HG4BertCascIntranuke.xml +++ b/config/HG4BertCascIntranuke.xml @@ -25,8 +25,24 @@ DelRNucleon double Yes mult. factor for nucleon de-Broglie wavelength - NUCL,INUKE,HNINUKE + NUCL genie::NuclearModelMap/Default + + 0.00 + 0.05 + 1.0 + 1.0 + 1.0 + 0.041 + 1.0 + 0.250 + true + true + true + + 0.0 + 0.0 + true diff --git a/config/NucBindEnergyAggregator.xml b/config/NucBindEnergyAggregator.xml index 40c9d6f8ba..90feee4743 100644 --- a/config/NucBindEnergyAggregator.xml +++ b/config/NucBindEnergyAggregator.xml @@ -18,7 +18,6 @@ AllowNuclRecombination bool Yes See [1] true --> - false diff --git a/src/Physics/HadronTransport/HG4BertCascIntranuke.cxx b/src/Physics/HadronTransport/HG4BertCascIntranuke.cxx index cd66d95893..5469ffcb73 100644 --- a/src/Physics/HadronTransport/HG4BertCascIntranuke.cxx +++ b/src/Physics/HadronTransport/HG4BertCascIntranuke.cxx @@ -51,6 +51,9 @@ #include "Geant4/G4InuclElementaryParticle.hh" #include "Geant4/G4InuclNuclei.hh" #include "Geant4/G4KineticTrackVector.hh" +#include "Geant4/G4Diproton.hh" +#include "Geant4/G4Dineutron.hh" +#include "Geant4/G4UnboundPN.hh" using std::ostringstream; @@ -63,14 +66,14 @@ using namespace genie::controls; //___________________________________________________________________________ HG4BertCascIntranuke::HG4BertCascIntranuke() - : EventRecordVisitorI("genie::HG4BertCascIntranuke") +: EventRecordVisitorI("genie::HG4BertCascIntranuke") { InitG4Particles(); } //___________________________________________________________________________ HG4BertCascIntranuke::HG4BertCascIntranuke(string config) - : EventRecordVisitorI("genie::HG4BertCascIntranuke", config) +: EventRecordVisitorI("genie::HG4BertCascIntranuke", config) { InitG4Particles(); } @@ -157,8 +160,7 @@ int HG4BertCascIntranuke::G4BertCascade(GHepRecord * evrec) const{ npdg = outgoingHadrons[l].getDefinition()->GetPDGEncoding(); G4LorentzVector hadP = outgoingHadrons[l].getMomentum(); TLorentzVector tempP(hadP.px(), hadP.py(), hadP.pz(), hadP.e() ); - GHepParticle new_particle(npdg, kIStStableFinalState, - 0, 1,-1,-1,tempP, remX); + GHepParticle new_particle(npdg, kIStStableFinalState, 0, 1,-1,-1,tempP, remX); evrec->AddParticle(new_particle); } @@ -179,19 +181,31 @@ int HG4BertCascIntranuke::G4BertCascade(GHepRecord * evrec) const{ // Get remnant momentum from cascade G4LorentzVector nucP = outgoingFragments[rem_index].getMomentum(); TLorentzVector remP(nucP.px(), nucP.py(), nucP.pz(), nucP.e() ); - npdg = outgoingFragments[rem_index].getDefinition()->GetPDGEncoding(); - GHepParticle largest_Fragment(npdg, kIStFinalStateNuclearRemnant, - 1,0,-1,-1, remP, remX); - evrec->AddParticle(largest_Fragment); - + npdg = outgoingFragments[rem_index].getDefinition()->GetPDGEncoding(); + //Checks if the remnant is present i PDGLibrary + TParticlePDG * pdgRemn=PDGLibrary::Instance()->Find(npdg,false); + + if(!pdgRemn) + { + LOG("HG4BertCascIntranuke", pINFO) + << "NO Particle with pdg = " << npdg << " in PDGLibrary!"; + // Add the particle with status id=15 and change it to HadroBlob + GHepParticle largest_Fragment(kPdgHadronicBlob, kIStFinalStateNuclearRemnant, + 1,0,-1,-1, remP, remX); + evrec->AddParticle(largest_Fragment); + } + else + { + GHepParticle largest_Fragment(npdg, kIStStableFinalState,1,-1,-1,-1, remP, remX); + evrec->AddParticle(largest_Fragment); + } // If any nuclear fragments left, add them to the event for (G4int k = 0; k < Nfrag; k++) { if (k != rem_index) { npdg = outgoingFragments[k].getDefinition()->GetPDGEncoding(); nucP = outgoingFragments[k].getMomentum(); // need to boost by fRemnP4 TLorentzVector tempP(nucP.px(), nucP.py(), nucP.pz(), nucP.e() ); - GHepParticle nuclear_Fragment(npdg, kIStStableFinalState, 0, 1,-1,-1, - tempP, remX); + GHepParticle nuclear_Fragment(npdg, kIStStableFinalState, 0, 1,-1,-1, tempP, remX); evrec->AddParticle(nuclear_Fragment); } } @@ -202,17 +216,14 @@ int HG4BertCascIntranuke::G4BertCascade(GHepRecord * evrec) const{ TLorentzVector p4h (0.,0.,probe->Pz(),probe->E()); TLorentzVector x4null(0.,0.,0.,0.); TLorentzVector p4tgt (0.,0.,0.,tgtNucl->Mass()); - evrec->AddParticle(probe->Pdg(), kIStStableFinalState, - 0,-1,-1,-1, p4h,x4null); - evrec->AddParticle(tgtNucl->Pdg(),kIStFinalStateNuclearRemnant, - 1,-1,-1,-1,p4tgt,x4null); + evrec->AddParticle(probe->Pdg(), kIStStableFinalState, 0,-1,-1,-1, p4h,x4null); + evrec->AddParticle(tgtNucl->Pdg(), kIStStableFinalState, 1,-1,-1,-1,p4tgt,x4null); } delete sp; return 0; } //___________________________________________________________________________ -double HG4BertCascIntranuke::GenerateStep(GHepRecord* /*evrec*/, - GHepParticle* p) const { +double HG4BertCascIntranuke::GenerateStep(GHepRecord* /*evrec*/, GHepParticle* p) const { // Generate a step (in fermis) for particle p in the input event. // Computes the mean free path L and generate an 'interaction' distance d // from an exp(-d/L) distribution @@ -235,8 +246,7 @@ void HG4BertCascIntranuke::ProcessEventRecord(GHepRecord* evrec) const // Check the event generation mode: should be lepton-nucleus fGMode = evrec->EventGenerationMode(); if ( fGMode== kGMdLeptonNucleus) { - LOG("HG4BertCascIntranuke", pINFO) - << " Lepton-nucleus event generation mode "; + LOG("HG4BertCascIntranuke", pINFO) << " Lepton-nucleus event generation mode "; GHepParticle* nucltgt = evrec->TargetNucleus(); if (nucltgt) { // Decide tracking radius for the current nucleus (few * R0 * A^1/3) @@ -249,12 +259,10 @@ void HG4BertCascIntranuke::ProcessEventRecord(GHepRecord* evrec) const // Collect hadrons from initial interaction and track them through the // nucleus using Bertini cascade TransportHadrons(evrec); - } else if(fGMode == kGMdHadronNucleus || - fGMode == kGMdPhotonNucleus){ - G4BertCascade(evrec); + } else if(fGMode == kGMdHadronNucleus || fGMode == kGMdPhotonNucleus){ + G4BertCascade(evrec); } else{ - LOG("HG4BertCascIntranuke", pINFO) - << " Inappropriate event generation mode - exit "; + LOG("HG4BertCascIntranuke", pINFO) << " Inappropriate event generation mode - exit "; return; } @@ -368,10 +376,8 @@ void HG4BertCascIntranuke::TransportHadrons(GHepRecord * evrec) const // In frozen nucleus approximation it is the remnant nucleus corrected for // final state lepton charge - //rwh-unused//GHepParticle* probe = evrec->Probe(); // incoming neutrino GHepParticle* tgtNucl = evrec->TargetNucleus(); GHepParticle* remNucl = evrec->RemnantNucleus(); - GHepParticle* outLept = evrec->FinalStatePrimaryLepton(); // outgoing lepton GHepParticle* struckNucleon = evrec->HitNucleon(); int inucl = evrec->RemnantNucleusPosition(); @@ -384,9 +390,9 @@ void HG4BertCascIntranuke::TransportHadrons(GHepRecord * evrec) const GHepParticle* incidentBaryon = 0; TObjArrayIter piter(evrec); TObjArrayIter pitter(evrec); + TObjArrayIter pittter(evrec); int icurr =-1; bool has_incidentBaryon(false), has_secondaries(false); - //rwh-unused// bool transparents(false), bool has_remnant(false), has_incidentparticle(false); fRemnA=remNucl->A(); @@ -425,16 +431,8 @@ void HG4BertCascIntranuke::TransportHadrons(GHepRecord * evrec) const //rwh-unused//transparents=true; } if ( ! has_incidentBaryon && sp->Status() == kIStHadronInTheNucleus ) { - /* - if (sp->Pdg() == kPdgProton || - sp->Pdg() == kPdgNeutron|| - sp->Pdg() == kPdgLambda || - sp->Pdg() == kPdgSigmaP || - sp->Pdg() == kPdgSigma0 || - sp->Pdg() == kPdgSigmaM ) { - */ if ( IsBaryon(sp) ) { - incidentBaryon = sp; + incidentBaryon = new GHepParticle(*sp); incidentDef = PDGtoG4Particle(sp->Pdg() ); has_incidentBaryon=true; } else { @@ -452,21 +450,13 @@ void HG4BertCascIntranuke::TransportHadrons(GHepRecord * evrec) const delete sp; } + int Nsec(0); + std::vector Postion_evrec; if ( has_secondaries ) { if ( ! incidentBaryon ) LOG("G4BertCascInterface::TransportHadrons", pINFO) << "Unrecognized baryon in nucleus"; - int Zinit = remNucl->Z() - outLept->Charge()/3; - if (incidentBaryon) { - Zinit += (struckNucleon->Charge() - incidentBaryon->Charge() )/3; - } - - //std::cout << " Zinit = " << Zinit <<" remNucl->Z()= " - // <Z()<< std::endl; - - //rwh-noused//int Ainit = remNucl->A(); - //std::cout << " Ainit = " << Ainit << std::endl; - + delete incidentBaryon; G4Fancy3DNucleus* g4Nucleus = new G4Fancy3DNucleus(); TLorentzVector pIncident; @@ -474,7 +464,7 @@ void HG4BertCascIntranuke::TransportHadrons(GHepRecord * evrec) const g4Nucleus->Init(remNucl->A(),remNucl->Z()); double EE = struckNucleon->E() - tgtNucl->Mass() + g4Nucleus->GetMass()*units::MeV; TLorentzVector struckMomentum(struckNucleon->Px(), struckNucleon->Py(), struckNucleon->Pz(), EE); - Double_t PxI(0),PyI(0),PzI(0),EEI(0); + Double_t PxI(0),PyI(0),PzI(0),EEI(0), Q(0), P(0), N(0); int icccur=-1; int pos_in_evrec(0); while( (p = (GHepParticle*) pitter.Next()) ) { @@ -484,16 +474,13 @@ void HG4BertCascIntranuke::TransportHadrons(GHepRecord * evrec) const PyI+=p->P4()->Py(); PzI+=p->P4()->Pz(); EEI+=p->P4()->E(); + Q+=p->Charge()/3; + if ( pos_in_evrec==0 ) pos_in_evrec = icccur; + Postion_evrec.push_back(icccur); + if (genie::pdg::IsProton(p->Pdg())) P++; + if (genie::pdg::IsNeutron(p->Pdg())) N++; if ( pos_in_evrec==0 ) pos_in_evrec = icccur; if ( ! has_incidentparticle ) { // take the baryon as incident particle - /* - if (p->Pdg() == kPdgProton || - p->Pdg() == kPdgNeutron || - p->Pdg() == kPdgLambda || - p->Pdg() == kPdgSigmaP || - p->Pdg() == kPdgSigma0 || - p->Pdg() == kPdgSigmaM ) { - */ if ( IsBaryon(p) ) { incidentDef = PDGtoG4Particle(p->Pdg() ); has_incidentparticle=true; @@ -501,10 +488,15 @@ void HG4BertCascIntranuke::TransportHadrons(GHepRecord * evrec) const } } } - GHepParticle * pinN = evrec->Particle(pos_in_evrec); if ( ! has_incidentparticle) { + GHepParticle * pinN = evrec->Particle(pos_in_evrec); incidentDef=PDGtoG4Particle(pinN->Pdg()); // if no baryon among the secondaries } + if (P == 2) incidentDef = PDGtoG4Particle(kPdgClusterPP); + else if (N == 2) incidentDef = PDGtoG4Particle(kPdgClusterNN); + else if (P == 1 && N == 1) incidentDef = PDGtoG4Particle(kPdgClusterNP); + + pIncident.SetPxPyPzE(PxI,PyI,PzI,EEI); @@ -517,6 +509,7 @@ void HG4BertCascIntranuke::TransportHadrons(GHepRecord * evrec) const // Create pseudo-particle to supply Bertini collider with bullet G4DynamicParticle dp(incidentDef, incidentDir, incidentKE/units::MeV, dynamicMass/units::MeV); + dp.SetCharge(Q); G4InuclElementaryParticle* incident = new G4InuclElementaryParticle(dp,G4InuclParticle::bullet); @@ -524,7 +517,8 @@ void HG4BertCascIntranuke::TransportHadrons(GHepRecord * evrec) const G4KineticTrackVector* g4secondaries = ConvertGenieSecondariesToG4(evrec); - int Nsec = g4secondaries->size(); + Nsec = g4secondaries->size(); + // Set up output and start the cascade G4CollisionOutput cascadeOutput; G4InuclCollider bertCollider; @@ -552,13 +546,18 @@ void HG4BertCascIntranuke::TransportHadrons(GHepRecord * evrec) const // Now the single hadrons int Nhad = cascadeOutput.numberOfOutgoingParticles(); const std::vector& outgoingHadrons = cascadeOutput.getOutgoingParticles(); + + int mother1=Postion_evrec.at(0); + int mother2(0); + if(Nsec==1)mother2=-1; + else if(Nsec>1)mother2=Postion_evrec.at(Nsec-1); for (int l = 0; l < Nhad; l++) { npdg = outgoingHadrons[l].getDefinition()->GetPDGEncoding(); G4LorentzVector hadP = outgoingHadrons[l].getMomentum(); TLorentzVector tempP(hadP.px(), hadP.py(), hadP.pz(), hadP.e() ); - GHepParticle new_particle(npdg, kIStStableFinalState, -1, -1,-1,-1,tempP, remX); + GHepParticle new_particle(npdg, kIStStableFinalState,mother1, mother2,-1,-1,tempP, remX); evrec->AddParticle(new_particle); } @@ -586,7 +585,7 @@ void HG4BertCascIntranuke::TransportHadrons(GHepRecord * evrec) const nucP = outgoingFragments[k].getMomentum(); // need to boost by fRemnP4 TLorentzVector tempP(nucP.px(), nucP.py(), nucP.pz(), nucP.e() ); - GHepParticle nuclear_Fragment(npdg, kIStStableFinalState, rem_nucl, 0,-1,-1,tempP, remX); + GHepParticle nuclear_Fragment(npdg, kIStStableFinalState, rem_nucl,-1,-1,-1,tempP, remX); evrec->AddParticle(nuclear_Fragment); } } @@ -597,8 +596,22 @@ void HG4BertCascIntranuke::TransportHadrons(GHepRecord * evrec) const remP.SetPy(remP.Py()+remNucl->P4()->Py()); remP.SetPz(remP.Pz()+remNucl->P4()->Pz()); - GHepParticle largest_Fragment(npdg, kIStFinalStateNuclearRemnant,rem_nucl,-1,-1,-1, remP, remX); - evrec->AddParticle(largest_Fragment); + //Checks if the remnant is present i PDGLibrary + TParticlePDG * pdgRemn=PDGLibrary::Instance()->Find(npdg,false); + if(!pdgRemn) + { + LOG("HG4BertCascIntranuke", pINFO) + << "NO Particle with pdg = " << npdg << " in PDGLibrary!"; + // Add the particle with status id=15 and change it to HadroBlob + GHepParticle largest_Fragment(kPdgHadronicBlob, kIStFinalStateNuclearRemnant, + rem_nucl,-1,-1,-1, remP, remX); + evrec->AddParticle(largest_Fragment); + } + else + { + GHepParticle largest_Fragment(npdg, kIStStableFinalState,rem_nucl,-1,-1,-1, remP, remX); + evrec->AddParticle(largest_Fragment); + } } // Nfrag > 0 has_remnant=true; } @@ -606,11 +619,23 @@ void HG4BertCascIntranuke::TransportHadrons(GHepRecord * evrec) const if(!has_remnant){ GHepParticle * sp = new GHepParticle(*evrec->Particle(inucl)); sp->SetFirstMother(inucl); - sp->SetStatus(kIStFinalStateNuclearRemnant); + sp->SetStatus(kIStStableFinalState); evrec->AddParticle(*sp); delete sp; } evrec->Particle(inucl)->SetStatus(kIStIntermediateState); + // Tests + int dau1(0), dau2(0); + if(Nsec>1){ + GHepParticle * pinN = evrec->Particle(Postion_evrec.at(0)); + dau1=pinN->FirstDaughter(); + dau2=pinN->LastDaughter(); + for(int ii=1;iiParticle(Postion_evrec.at(ii))->SetFirstDaughter(dau1); + evrec->Particle(Postion_evrec.at(ii))->SetLastDaughter(dau2); + } + } + // Geant4 conservation test // this probably isn't configured right ... skip it for now @@ -676,12 +701,16 @@ const G4ParticleDefinition* HG4BertCascIntranuke::PDGtoG4Particle(int pdg) const { const G4ParticleDefinition* pDef = 0; + if (pdg == kPdgClusterPP) return G4Diproton::Diproton(); + if (pdg == kPdgClusterNN) return G4Dineutron::Dineutron(); + if (pdg == kPdgClusterNP) return G4UnboundPN::UnboundPN(); + if ( abs(pdg) < 1000000000 ) { pDef = G4ParticleTable::GetParticleTable()->FindParticle(pdg); } else if ( pdg < 2000000000 ) { pDef = G4IonTable::GetIonTable()->GetIon(pdg); } - + if ( ! pDef ) { LOG("HG4BertCascIntranuke", pWARN) << "Unrecognized Bertini particle type: " << pdg; @@ -772,12 +801,10 @@ bool HG4BertCascIntranuke::Conserve4Momentum(GHepRecord* evrec) const p = (GHepParticle*) (*evrec)[j]; if (p->FirstMother() == evrec->RemnantNucleusPosition() ) { sum4mom += *(p->P4() ); - /* LOG("HG4BertCascIntranuke", pINFO) << " Final state 4-momentum = (" << p->P4()->Px() << ", " << p->P4()->Py() << ", " << p->P4()->Pz() << ", " << p->P4()->E() << ") "; - */ } } sum4mom += *(finalLepton->P4() ); diff --git a/src/Physics/HadronTransport/HINCLCascadeIntranuke.cxx b/src/Physics/HadronTransport/HINCLCascadeIntranuke.cxx index 63e254fc0e..5f21a313fc 100644 --- a/src/Physics/HadronTransport/HINCLCascadeIntranuke.cxx +++ b/src/Physics/HadronTransport/HINCLCascadeIntranuke.cxx @@ -74,20 +74,20 @@ using std::ostringstream; using namespace std; HINCLCascadeIntranuke::HINCLCascadeIntranuke() : - EventRecordVisitorI("genie::HINCLCascadeIntranuke"), - theINCLConfig(0), theINCLModel(0), theDeExcitation(0) +EventRecordVisitorI("genie::HINCLCascadeIntranuke"), +theINCLConfig(0), theINCLModel(0), theDeExcitation(0) { LOG("HINCLCascadeIntranuke", pDEBUG) - << "default ctor"; + << "default ctor"; } //______________________________________________________________________________ HINCLCascadeIntranuke::HINCLCascadeIntranuke(string config) : - EventRecordVisitorI("genie::HINCLCascadeIntranuke", config), - theINCLConfig(0), theINCLModel(0), theDeExcitation(0) +EventRecordVisitorI("genie::HINCLCascadeIntranuke", config), +theINCLConfig(0), theINCLModel(0), theDeExcitation(0) { LOG("HINCLCascadeIntranuke", pDEBUG) - << "ctor from config " << config; + << "ctor from config " << config; } //______________________________________________________________________________ @@ -158,7 +158,7 @@ void HINCLCascadeIntranuke::LoadConfig(void) } LOG("HINCLCascadeIntranuke", pDEBUG) - << "LoadConfig() create theINCLConfig"; + << "LoadConfig() create theINCLConfig"; theINCLConfig = theParser.parse(nflags,flags); // there's currently no way to update *all* of the data paths in the Config @@ -171,7 +171,7 @@ void HINCLCascadeIntranuke::LoadConfig(void) } LOG("HINCLCascadeIntranuke", pDEBUG) - << "doCascade new G4INCL::INCL"; + << "doCascade new G4INCL::INCL"; theINCLModel = new G4INCL::INCL(theINCLConfig); // code copied fomr INCL's main/src/INCLCascade.cc @@ -181,43 +181,43 @@ void HINCLCascadeIntranuke::LoadConfig(void) case G4INCL::DeExcitationABLAv3p: theDeExcitation = new G4INCLAblaInterface(theINCLConfig); LOG("HINCLCascadeIntranuke", pINFO) - << "using ABLAv3p for DeExcitation"; + << "using ABLAv3p for DeExcitation"; break; #endif #ifdef INCL_DEEXCITATION_ABLA07 case G4INCL::DeExcitationABLA07: theDeExcitation = new ABLA07CXX::Abla07Interface(theINCLConfig); LOG("HINCLCascadeIntranuke", pINFO) - << "using ABLA07CXX for DeExcitation"; + << "using ABLA07CXX for DeExcitation"; break; #endif #ifdef INCL_DEEXCITATION_SMM case G4INCL::DeExcitationSMM: theDeExcitation = new SMMCXX::SMMInterface(theINCLConfig); LOG("HINCLCascadeIntranuke", pINFO) - << "using SMMCXX for DeExcitation"; + << "using SMMCXX for DeExcitation"; break; #endif #ifdef INCL_DEEXCITATION_GEMINIXX case G4INCL::DeExcitationGEMINIXX: theDeExcitation = new G4INCLGEMINIXXInterface(theINCLConfig); LOG("HINCLCascadeIntranuke", pINFO) - << "using GEMINIXX for DeExcitation"; + << "using GEMINIXX for DeExcitation"; break; #endif default: std::stringstream ss; ss << "########################################################\n" - << "### WARNING WARNING WARNING ###\n" - << "### ###\n" - << "### You are running the code without any coupling to ###\n" - << "### a de-excitation model! ###\n" - << "### Results will be INCOMPLETE and UNPHYSICAL! ###\n" - << "### Are you sure this is what you want to do? ###\n" - << "########################################################\n"; + << "### WARNING WARNING WARNING ###\n" + << "### ###\n" + << "### You are running the code without any coupling to ###\n" + << "### a de-excitation model! ###\n" + << "### Results will be INCOMPLETE and UNPHYSICAL! ###\n" + << "### Are you sure this is what you want to do? ###\n" + << "########################################################\n"; INCL_INFO(ss.str()); LOG("HINCLCascadeIntranuke", pWARN) - << '\n' << ss.str(); + << '\n' << ss.str(); //std::cout << ss.str(); break; } @@ -242,61 +242,61 @@ bool HINCLCascadeIntranuke::AddDataPathFlags(size_t& nflags, char** flags) { std::vector datapaths; LOG("HINCLCascadeIntranuke", pINFO) - << "check for data location of INCL"; + << "check for data location of INCL"; datapaths.clear(); datapaths.push_back(theINCLConfig->getINCLXXDataFilePath()); datapaths.push_back("!INCL-incl-data-alt1"); datapaths.push_back("!INCL-incl-data-alt2"); needed_update |= - LookForAndAddValidPath(datapaths,0,"--inclxx-datafile-path",nflags,flags); + LookForAndAddValidPath(datapaths,0,"--inclxx-datafile-path",nflags,flags); switch(theINCLConfig->getDeExcitationType()) { #ifdef INCL_DEEXCITATION_ABLAXX case G4INCL::DeExcitationABLAv3p: LOG("HINCLCascadeIntranuke", pINFO) - << "using ABLAv3p for DeExcitation -- check for data location"; + << "using ABLAv3p for DeExcitation -- check for data location"; datapaths.clear(); datapaths.push_back(theINCLConfig->getABLAv3pCxxDataFilePath()); datapaths.push_back("!INCL-ablav3p-data-alt1"); datapaths.push_back("!INCL-ablav3p-data-alt2"); needed_update |= - LookForAndAddValidPath(datapaths,0,"--ablav3p-cxx-datafile-path",nflags,flags); + LookForAndAddValidPath(datapaths,0,"--ablav3p-cxx-datafile-path",nflags,flags); break; #endif #ifdef INCL_DEEXCITATION_ABLA07 case G4INCL::DeExcitationABLA07: LOG("HINCLCascadeIntranuke", pINFO) - << "using ABLA07 for DeExcitation -- check for data location"; + << "using ABLA07 for DeExcitation -- check for data location"; datapaths.clear(); datapaths.push_back(theINCLConfig->getABLA07DataFilePath()); datapaths.push_back("!INCL-abla07-data-alt1"); datapaths.push_back("!INCL-abla07-data-alt2"); needed_update |= - LookForAndAddValidPath(datapaths,0,"--abla07-datafile-path",nflags,flags); + LookForAndAddValidPath(datapaths,0,"--abla07-datafile-path",nflags,flags); break; #endif #ifdef INCL_DEEXCITATION_SMM case G4INCL::DeExcitationSMM: LOG("HINCLCascadeIntranuke", pINFO) - << "using SMMCXX for DeExcitation -- no data files to check for"; + << "using SMMCXX for DeExcitation -- no data files to check for"; break; #endif #ifdef INCL_DEEXCITATION_GEMINIXX case G4INCL::DeExcitationGEMINIXX: LOG("HINCLCascadeIntranuke", pINFO) - << "using GEMINIXX for DeExcitation -- check for data location"; + << "using GEMINIXX for DeExcitation -- check for data location"; datapaths.clear(); datapaths.push_back(theINCLConfig->getGEMINIXXDataFilePath()); datapaths.push_back("!INCL-gemini-data-alt1"); datapaths.push_back("!INCL-gemini-data-alt2"); needed_update |= - LookForAndAddValidPath(datapaths,0,"--geminixx-datafile-path",nflags,flags); + LookForAndAddValidPath(datapaths,0,"--geminixx-datafile-path",nflags,flags); break; #endif default: LOG("HINCLCascadeIntranuke", pINFO) - << "using no DeExcitation -- no data files to check for"; + << "using no DeExcitation -- no data files to check for"; break; } @@ -307,17 +307,17 @@ bool HINCLCascadeIntranuke::AddDataPathFlags(size_t& nflags, char** flags) { ss << "[" << setw(3) << i << "] " << flags[i] << '\n'; } LOG("HINCLCascadeIntranuke", pNOTICE) - << "Flags passed to INCLConfigParser" - << '\n' << ss.str(); + << "Flags passed to INCLConfigParser" + << '\n' << ss.str(); return needed_update; } //______________________________________________________________________________ bool HINCLCascadeIntranuke::LookForAndAddValidPath(std::vector& datapaths, - size_t defaultIndx, - const char* optString, - size_t& nflags, char** flags) { + size_t defaultIndx, + const char* optString, + size_t& nflags, char** flags) { // okay, we have a series of paths _OR_ parameter names ("!" as first char) // loop over possibilities @@ -332,15 +332,15 @@ bool HINCLCascadeIntranuke::LookForAndAddValidPath(std::vector& dat bool needed_update = false; LOG("HINCLCascadeIntranuke", pINFO) - << "looking for a valid path for " << optString - << " (default [" << defaultIndx << "]"; + << "looking for a valid path for " << optString + << " (default [" << defaultIndx << "]"; size_t foundIndx = SIZE_MAX; // flag as unfound size_t npaths = datapaths.size(); for (size_t ipath = 0; ipath < npaths; ++ipath) { std::string& apath = datapaths[ipath]; LOG("HINCLCascadeIntranuke", pINFO) - << "looking at " << apath; + << "looking at " << apath; if ( apath.at(0) == '!' ) { apath.erase(0,1); // remove the "!" // parameter name, returned value, default value @@ -349,19 +349,19 @@ bool HINCLCascadeIntranuke::LookForAndAddValidPath(std::vector& dat std::string notfoundvalue = std::string("no-such-param-") + apath; GetParamDef(apath,newpath,notfoundvalue); LOG("HINCLCascadeIntranuke", pINFO) - << "fetch param "<< "[" << ipath << "] "<< apath << " got " << newpath; + << "fetch param "<< "[" << ipath << "] "<< apath << " got " << newpath; apath = newpath; } const char* expandedPath = gSystem->ExpandPathName(apath.c_str()); if ( ! expandedPath ) { LOG("HINCLCascadeIntranuke", pINFO) - << "expandedPath was NULL"; + << "expandedPath was NULL"; continue; } bool valid = utils::system::DirectoryExists(expandedPath); LOG("HINCLCascadeIntranuke", pINFO) - << "expandedPath " << expandedPath << " " - << ((valid)?"valid":"doesn't exist"); + << "expandedPath " << expandedPath << " " + << ((valid)?"valid":"doesn't exist"); if ( valid ) { apath = std::string(expandedPath); foundIndx = ipath; @@ -378,8 +378,8 @@ bool HINCLCascadeIntranuke::LookForAndAddValidPath(std::vector& dat ss << "[" << ipath << "] " << apath << "\n"; } LOG("HINCLCascadeIntranuke", pWARN) - << "no valid path found for " << optString - << ", tried: \n" << ss.str(); + << "no valid path found for " << optString + << ", tried: \n" << ss.str(); } else { // add the flag with the valid path flags[nflags] = strdup(optString); ++nflags; @@ -428,23 +428,45 @@ int HINCLCascadeIntranuke::doCascade(GHepRecord * evrec) const { if ( result.transparent ) { evrec->AddParticle(pdg_codeProbe, ist1, 0,-1,-1,-1, p4h,x4null); - evrec->AddParticle(pdg_codeTarget,kIStFinalStateNuclearRemnant,1,-1,-1,-1,p4tgt,x4null); + evrec->AddParticle(pdg_codeTarget,kIStStableFinalState,1,-1,-1,-1,p4tgt,x4null); INCL_DEBUG("Transparent event" << std::endl); } else { INCL_DEBUG("Number of produced particles: " << result.nParticles << "\n"); if ( theDeExcitation != 0 ) { theDeExcitation->deExcite(&result); } - int mom = 1; + for (int nP = 0; nP < result.nParticles; nP++){ if ( nP == result.nParticles-1 ) { - GHepParticle *p_outR = - INCLtoGenieParticle(result,nP,kIStFinalStateNuclearRemnant,mom,-1); - evrec->AddParticle(*p_outR); - delete p_outR; + int pdgRem=INCLtopdgcode(result.A[nP],result.Z[nP]); + TParticlePDG * pdgRemn=PDGLibrary::Instance()->Find(pdgRem,false); + if(!pdgRemn) + { + LOG("HINCLCascadeIntranuke", pINFO) + << "NO Particle with pdg = " << pdgRem << " in PDGLibrary!"; + // Add the particle with status id=15 and change it to HadroBlob + TVector3 p3M(result.px[nP]/1000,result.py[nP]/1000,result.pz[nP]/1000); + + double MassRem=0.5*((result.px[nP])*(result.px[nP]) + + (result.py[nP])*(result.py[nP]) + + (result.pz[nP])*(result.pz[nP]) - + result.EKin[nP]*result.EKin[nP]) / (result.EKin[nP]); + float ERem=result.EKin[nP]+MassRem; + TLorentzVector p4tgtf(p3M,ERem/1000); + GHepParticle p_outR(kPdgHadronicBlob, kIStFinalStateNuclearRemnant, + 1,-1,-1,-1, p4tgtf, x4null);; + evrec->AddParticle(p_outR); + } + else + { + GHepParticle *p_outR = + INCLtoGenieParticle(result,nP,kIStStableFinalState,1,-1); + evrec->AddParticle(*p_outR); + delete p_outR; + } } else { GHepParticle *p_outR = - INCLtoGenieParticle(result,nP,kIStStableFinalState,0,-1); + INCLtoGenieParticle(result,nP,kIStStableFinalState,0,-1); evrec->AddParticle(*p_outR); delete p_outR; } @@ -457,17 +479,17 @@ int HINCLCascadeIntranuke::doCascade(GHepRecord * evrec) const { void HINCLCascadeIntranuke::ProcessEventRecord(GHepRecord * evrec) const { LOG("HINCLCascadeIntranuke", pNOTICE) - << "************ Running HINCLCascadeIntranuke MODE INTRANUKE ************"; + << "************ Running HINCLCascadeIntranuke MODE INTRANUKE ************"; fGMode = evrec->EventGenerationMode(); if ( fGMode == kGMdHadronNucleus || - fGMode == kGMdPhotonNucleus ) { + fGMode == kGMdPhotonNucleus ) { HINCLCascadeIntranuke::doCascade(evrec); - } else if ( fGMode == kGMdLeptonNucleus ) { - HINCLCascadeIntranuke::TransportHadrons(evrec); - } +} else if ( fGMode == kGMdLeptonNucleus ) { + HINCLCascadeIntranuke::TransportHadrons(evrec); +} - LOG("HINCLCascadeIntranuke", pINFO) << "Done with this event"; +LOG("HINCLCascadeIntranuke", pINFO) << "Done with this event"; } bool HINCLCascadeIntranuke::CanRescatter(const GHepParticle * p) const { @@ -476,76 +498,82 @@ bool HINCLCascadeIntranuke::CanRescatter(const GHepParticle * p) const { // rescattered by this cascade MC assert(p); return ( p->Pdg() == kPdgPiP || - p->Pdg() == kPdgPiM || - p->Pdg() == kPdgPi0 || - p->Pdg() == kPdgProton || - p->Pdg() == kPdgNeutron - ); + p->Pdg() == kPdgPiM || + p->Pdg() == kPdgPi0 || + p->Pdg() == kPdgProton || + p->Pdg() == kPdgNeutron + ); } void HINCLCascadeIntranuke::TransportHadrons(GHepRecord * evrec) const { int inucl = -1; if ( fGMode == kGMdHadronNucleus || - fGMode == kGMdPhotonNucleus ) { + fGMode == kGMdPhotonNucleus ) { inucl = evrec->TargetNucleusPosition(); - } else { - if ( fGMode == kGMdLeptonNucleus || - fGMode == kGMdNucleonDecay || - fGMode == kGMdNeutronOsc ) { - inucl = evrec->RemnantNucleusPosition(); - } - } +} else { + if ( fGMode == kGMdLeptonNucleus || + fGMode == kGMdNucleonDecay || + fGMode == kGMdNeutronOsc ) { + inucl = evrec->RemnantNucleusPosition(); +} +} - LOG("HINCLCascadeIntranuke", pNOTICE) - << "Propagating hadrons within nucleus found in position = " << inucl; - int tpos = evrec->TargetNucleusPosition(); - GHepParticle * target = evrec->Particle(tpos); - GHepParticle * nucl = evrec->Particle(inucl); - if ( ! nucl ) { - LOG("HINCLCascadeIntranuke", pERROR) - << "No nucleus found in position = " << inucl; - LOG("HINCLCascadeIntranuke", pERROR) - << *evrec; - return; - } - fRemnA = nucl->A(); - fRemnZ = nucl->Z(); - int A_f(0), Z_f(0), Aft(0), A_i(target->A()),Z_i(target->Z()); +LOG("HINCLCascadeIntranuke", pNOTICE) +<< "Propagating hadrons within nucleus found in position = " << inucl; +int tpos = evrec->TargetNucleusPosition(); +GHepParticle * target = evrec->Particle(tpos); +GHepParticle * nucl = evrec->Particle(inucl); +if ( ! nucl ) { + LOG("HINCLCascadeIntranuke", pERROR) + << "No nucleus found in position = " << inucl; + LOG("HINCLCascadeIntranuke", pERROR) + << *evrec; + return; +} +fRemnA = nucl->A(); +fRemnZ = nucl->Z(); +GHepParticle * Incident_particle = evrec->Particle(0); - LOG("HINCLCascadeIntranuke", pNOTICE) - << "Nucleus (A,Z) = (" << fRemnA << ", " << fRemnZ << ")"; +int A_f(0), Z_f(0), Aft(0), A_i(target->A()),Z_i(0), Charge_probe(Incident_particle->Charge()); +if(Charge_probe==0) Z_i=target->Z(); +else if(Charge_probe<0) Z_i=target->Z()-1; +else if(Charge_probe>0)Z_i=target->Z()+1; - const TLorentzVector & p4nucl = *(nucl->P4()); - TLorentzVector x4null(0.,0.,0.,0.); - fRemnP4 = p4nucl; - TObjArrayIter piter(evrec); - GHepParticle * p = 0; +LOG("HINCLCascadeIntranuke", pNOTICE) +<< "Nucleus (A,Z) = (" << fRemnA << ", " << fRemnZ << ")"; + +const TLorentzVector & p4nucl = *(nucl->P4()); +TLorentzVector x4null(0.,0.,0.,0.); +fRemnP4 = p4nucl; + +TObjArrayIter piter(evrec); +GHepParticle * p = 0; - int icurr = -1; +int icurr = -1; - bool is_QE = evrec->Summary()->ProcInfo().IsQuasiElastic(); +bool is_QE = evrec->Summary()->ProcInfo().IsQuasiElastic(); - TLorentzVector * p_4 = nucl->P4(); +TLorentzVector * p_4 = nucl->P4(); // momentum of the remnant nucleus. - double pxRemn = p_4->Px(); - double pyRemn = p_4->Py(); - double pzRemn = p_4->Pz(); - int pdg_codeTargetRemn= genie::pdg::IonPdgCode(nucl->A(),nucl->Z()); - TLorentzVector p4tgf(p_4->Px(),p_4->Py(),p_4->Pz(),0.0); +double pxRemn = p_4->Px(); +double pyRemn = p_4->Py(); +double pzRemn = p_4->Pz(); +int pdg_codeTargetRemn= genie::pdg::IonPdgCode(nucl->A(),nucl->Z()); +TLorentzVector p4tgf(p_4->Px(),p_4->Py(),p_4->Pz(),0.0); // Loop over GHEP and run intranuclear rescattering on handled particles - std::vectorListeOfINCLresult; - std::vector Postion_evrec; +std::vectorListeOfINCLresult; +std::vector Postion_evrec,num_of_AZexception; GHepParticle * fsl = evrec->FinalStatePrimaryLepton(); // primary Lepton double ExcitaionE(0), the_pxRemn(0), the_pyRemn(0), the_pzRemn(0); - int Zl(0), Aresult(0), Zresult(0), Aexception(0), Zexception(0), - Pos(0), theA_Remn(0), theZ_Remn(0); + int Zl(0), Aresult(0), Zresult(0), Aexception(0), Zexception(0),Pos(0), + mother1(0),mother2(0),theA_Remn(0), theZ_Remn(0); if ( fsl->Charge() == 0. ) Zl = 0; - if ( fsl->Charge() < 0. ) Zl = -1; - if ( fsl->Charge() > 0. ) Zl = 1; + else if ( fsl->Charge() < 0. ) Zl = -1; + else if ( fsl->Charge() > 0. ) Zl = 1; bool has_remnant = false; while ( (p = (GHepParticle *) piter.Next() ) ) { @@ -563,11 +591,25 @@ void HINCLCascadeIntranuke::TransportHadrons(GHepRecord * evrec) const { // if I can't rescatter it, I will just take it out of the nucleus LOG("HINCLCascadeIntranuke", pNOTICE) - << "... Current version can't rescatter a " << sp->Name(); + << "... Current version can't rescatter a " << sp->Name(); sp->SetFirstMother(icurr); sp->SetStatus(kIStStableFinalState); - evrec->AddParticle(*sp); - delete sp; + if ( sp->Charge() == 0. ) { + Zl+=0; + Aft+=1; + } + else if ( sp->Charge() < 0. ) { + Zl-=1; + Aft+=1; + } + else if ( sp->Charge() > 0. ) { + Zl+=1; + Aft+=1; + } + + evrec->AddParticle(*sp); + evrec->Particle(sp->FirstMother())->SetRescatterCode(1); + delete sp; continue; // <-- skip to next GHEP entry } @@ -593,6 +635,7 @@ void HINCLCascadeIntranuke::TransportHadrons(GHepRecord * evrec) const { sp->SetFirstMother(icurr); sp->SetStatus(kIStStableFinalState); evrec->AddParticle(*sp); + evrec->Particle(sp->FirstMother())->SetRescatterCode(1); delete sp; continue; } @@ -608,14 +651,15 @@ void HINCLCascadeIntranuke::TransportHadrons(GHepRecord * evrec) const { theSpecies.theZ=pdgcpiontoZ(sp->Pdg()); G4INCL::Particle *pincl = - new G4INCL::Particle( theType , E , momentum, thePosition); + new G4INCL::Particle( theType , E , momentum, thePosition); G4INCL::Random::SeedVector const theInitialSeeds = - G4INCL::Random::getSeeds(); + G4INCL::Random::getSeeds(); G4INCL::Random::saveSeeds(); G4INCL::EventInfo result; + result=theINCLModel->processEvent(theSpecies,pincl,EKin,fRemnA,fRemnZ,"FSI"); // Exception get remnant with Z and A <0 @@ -631,35 +675,51 @@ void HINCLCascadeIntranuke::TransportHadrons(GHepRecord * evrec) const { sp->SetFirstMother(icurr); sp->SetStatus(kIStStableFinalState); evrec->AddParticle(*sp); + evrec->Particle(sp->FirstMother())->SetRescatterCode(1); delete sp; for (size_t it=0; itParticle(Pos); - theA_Remn+= (fRemnA + pdgcpiontoA(pinthenucleus->Pdg())- ListeOfINCLresult.at(it).ARem[0]); - theZ_Remn+= (fRemnZ + pdgcpiontoZ(pinthenucleus->Pdg())- ListeOfINCLresult.at(it).ZRem[0]); - if ( (A_i-theA_Remn-Aft) < (Z_i-theZ_Remn-Zl) ) { - theA_Remn-= (fRemnA + pdgcpiontoA(pinthenucleus->Pdg())- ListeOfINCLresult.at(it).ARem[0]); - theZ_Remn-= (fRemnZ + pdgcpiontoZ(pinthenucleus->Pdg())- ListeOfINCLresult.at(it).ZRem[0]); - Zl+=pdgcpiontoZ(pinthenucleus->Pdg()); - Aft+=pdgcpiontoA(pinthenucleus->Pdg()); - pinthenucleus->SetFirstMother(Pos); - pinthenucleus->SetStatus(kIStStableFinalState); - evrec->AddParticle(*pinthenucleus); - AZexception=true; - } else { - the_pxRemn+=ListeOfINCLresult.at(it).pxRem[0]; - the_pyRemn+=ListeOfINCLresult.at(it).pyRem[0]; - the_pzRemn+=ListeOfINCLresult.at(it).pzRem[0]; - ExcitaionE+=ListeOfINCLresult.at(it).EStarRem[0]; - } - } - if ( AZexception ) ListeOfINCLresult.pop_back(); + Pos=Postion_evrec.at(it);// Position of the mother in evrec + + GHepParticle * pinthenucleus = evrec->Particle(Pos); + theA_Remn+= (fRemnA + pdgcpiontoA(pinthenucleus->Pdg())- ListeOfINCLresult.at(it).ARem[0]); + theZ_Remn+= (fRemnZ + pdgcpiontoZ(pinthenucleus->Pdg())- ListeOfINCLresult.at(it).ZRem[0]); + if ( (A_i-theA_Remn-Aft) < (Z_i-theZ_Remn-Zl) ) { + theA_Remn-= (fRemnA + pdgcpiontoA(pinthenucleus->Pdg())- ListeOfINCLresult.at(it).ARem[0]); + theZ_Remn-= (fRemnZ + pdgcpiontoZ(pinthenucleus->Pdg())- ListeOfINCLresult.at(it).ZRem[0]); + Zl+=pdgcpiontoZ(pinthenucleus->Pdg()); + Aft+=pdgcpiontoA(pinthenucleus->Pdg()); + + pinthenucleus->SetFirstMother(Pos); + pinthenucleus->SetStatus(kIStStableFinalState); + evrec->AddParticle(*pinthenucleus); + evrec->Particle(pinthenucleus->FirstMother())->SetRescatterCode(1); + AZexception=true; + num_of_AZexception.push_back(it); + } else { + the_pxRemn+=ListeOfINCLresult.at(it).pxRem[0]; + the_pyRemn+=ListeOfINCLresult.at(it).pyRem[0]; + the_pzRemn+=ListeOfINCLresult.at(it).pzRem[0]; + ExcitaionE+=ListeOfINCLresult.at(it).EStarRem[0]; + } + } + if (AZexception) { + for (size_t it=0;it1) mom2=Postion_evrec.at(Number_of_Sec-1); + for (size_t it=0; itAddParticle(*p_outD); delete p_outD; } //Add result without the remnant nucleus @@ -673,39 +733,54 @@ void HINCLCascadeIntranuke::TransportHadrons(GHepRecord * evrec) const { ListeOfINCLresult.at(it).EStarRem[0]=ExcitaionE; theDeExcitation->deExcite(&ListeOfINCLresult.at(it)); for (int nP=0;nPAddParticle(*p_outFinal); - delete p_outFinal; - has_remnant=true; + int rem_index=FindlargestFragment(ListeOfINCLresult.at(it)); + if(nP==rem_index||ListeOfINCLresult.at(it).A[nP]>1) { + mom1=inucl; + mom2=-1; + } + else{ + mom1=Postion_evrec.at(0); + if(Number_of_Sec==1) mom2=-1; + if(Number_of_Sec>1) mom2=Postion_evrec.at(Number_of_Sec-1); + } + GHepParticle *p_outFinal = INCLtoGenieParticle(ListeOfINCLresult.at(it), + nP,kIStStableFinalState,mom1,mom2); + evrec->AddParticle(*p_outFinal); + delete p_outFinal; + has_remnant=true; }//Add all the result with the correct remnant nucleus } } ListeOfINCLresult.clear(); - } else { + num_of_AZexception.clear(); + }// + } else { //if result give a transparent event FSI=1 // Store *sp - if ( result.transparent ) { - Zl+=pdgcpiontoZ(sp->Pdg()); - Aft+=pdgcpiontoA(sp->Pdg()); - sp->SetFirstMother(icurr); - sp->SetStatus(kIStStableFinalState); - evrec->AddParticle(*sp); - delete sp; - } else { - Postion_evrec.push_back(icurr); - ListeOfINCLresult.push_back(result); - delete sp; - } - + if ( result.transparent ) { + Zl+=pdgcpiontoZ(sp->Pdg()); + Aft+=pdgcpiontoA(sp->Pdg()); + //sp->SetFirstMother(icurr); + sp->SetStatus(kIStStableFinalState); + evrec->AddParticle(*sp); + evrec->Particle(sp->FirstMother())->SetRescatterCode(1); + delete sp; + } else { + Postion_evrec.push_back(icurr); + ListeOfINCLresult.push_back(result); + delete sp; } + } + } //Ghep-entry if ( ListeOfINCLresult.size() != 0 ) { - bool AZexception(false); + bool AZexception=false; for (size_t it=0; itParticle(Pos); theA_Remn+= (fRemnA + pdgcpiontoA(pinthenucleus->Pdg())- ListeOfINCLresult.at(it).ARem[0]); theZ_Remn+= (fRemnZ + pdgcpiontoZ(pinthenucleus->Pdg())- ListeOfINCLresult.at(it).ZRem[0]); @@ -714,10 +789,12 @@ void HINCLCascadeIntranuke::TransportHadrons(GHepRecord * evrec) const { theZ_Remn-= (fRemnZ + pdgcpiontoZ(pinthenucleus->Pdg())- ListeOfINCLresult.at(it).ZRem[0]); Zl+=pdgcpiontoZ(pinthenucleus->Pdg()); Aft+=pdgcpiontoA(pinthenucleus->Pdg()); + pinthenucleus->SetFirstMother(Pos); pinthenucleus->SetStatus(kIStStableFinalState); evrec->AddParticle(*pinthenucleus); AZexception=true; + num_of_AZexception.push_back(it); } else { the_pxRemn+=ListeOfINCLresult.at(it).pxRem[0]; the_pyRemn+=ListeOfINCLresult.at(it).pyRem[0]; @@ -725,45 +802,77 @@ void HINCLCascadeIntranuke::TransportHadrons(GHepRecord * evrec) const { ExcitaionE+=ListeOfINCLresult.at(it).EStarRem[0]; } } - if (AZexception) ListeOfINCLresult.pop_back(); - for (size_t it=0; it < ListeOfINCLresult.size(); it++) { - if ( is_QE) { + if (AZexception) { + for (size_t it=0;itdeExcite(&ListeOfINCLresult.at(it)); - for (int nP=0; nP < ListeOfINCLresult.at(it).nParticles; nP++ ) { - GHepParticle *p_out = INCLtoGenieParticle(ListeOfINCLresult.at(it), - nP,kIStStableFinalState,Pos,inucl); - evrec->AddParticle(*p_out); - delete p_out; + ListeOfINCLresult.at(it).pxRem[0] += pxRemn*1000; + ListeOfINCLresult.at(it).pyRem[0] += pyRemn*1000; + ListeOfINCLresult.at(it).pzRem[0] += 1000*pzRemn; + if ( theDeExcitation != 0 ) theDeExcitation->deExcite(&ListeOfINCLresult.at(it)); + int rem_index=FindlargestFragment(ListeOfINCLresult.at(it)); + + for (int nP=0; nP < ListeOfINCLresult.at(it).nParticles; nP++ ) { + // if(nP==ListeOfINCLresult.at(it).nParticles-1) Pos=inucl; + if(nP==rem_index) Pos=inucl; + else if(nP!=rem_index) Pos=Postion_evrec.at(it); + GHepParticle *p_out = INCLtoGenieParticle(ListeOfINCLresult.at(it), + nP,kIStStableFinalState,Pos,-1); + evrec->AddParticle(*p_out); + delete p_out; }// Add to evrec the result } else { if ( it < ListeOfINCLresult.size()-1 ) { for (int nP=0; nP < ListeOfINCLresult.at(it).nParticles; nP++ ) { A_f+=ListeOfINCLresult.at(it).A[nP]; Z_f+=ListeOfINCLresult.at(it).Z[nP]; - GHepParticle *p_outD = - INCLtoGenieParticle(ListeOfINCLresult.at(it),nP,kIStStableFinalState,Pos,inucl); - evrec->AddParticle(*p_outD); - delete p_outD; - } //Add result without the remnant nucleus - } else { - ListeOfINCLresult.at(it).ARem[0]=A_i-theA_Remn-Aft; - ListeOfINCLresult.at(it).ZRem[0]=Z_i-theZ_Remn-Zl; - ListeOfINCLresult.at(it).pxRem[0]= the_pxRemn + (pxRemn*1000); - ListeOfINCLresult.at(it).pyRem[0]= the_pyRemn + (pyRemn*1000); - ListeOfINCLresult.at(it).pzRem[0]= the_pzRemn + (1000*pzRemn); - ListeOfINCLresult.at(it).EStarRem[0]=ExcitaionE; - if ( theDeExcitation != 0 ) theDeExcitation->deExcite(&ListeOfINCLresult.at(it)); - for (int nP=0; nP < ListeOfINCLresult.at(it).nParticles; nP++ ) { - A_f+=ListeOfINCLresult.at(it).A[nP]; - Z_f+=ListeOfINCLresult.at(it).Z[nP]; - GHepParticle *p_outR = INCLtoGenieParticle(ListeOfINCLresult.at(it), - nP,kIStStableFinalState,Pos,inucl); - evrec->AddParticle(*p_outR); - delete p_outR; + if(ListeOfINCLresult.at(it).A[nP]>1) { // if the event emits cluster during cascade + mother1=inucl; + mother2=-1; + } + else{ + mother1=Postion_evrec.at(0); + if(ListeOfINCLresult.size()==1) mother2=-1; + else if(ListeOfINCLresult.size()>1){ + mother2=Postion_evrec.at(ListeOfINCLresult.size()-1); + } + } + GHepParticle *p_outD = + INCLtoGenieParticle(ListeOfINCLresult.at(it),nP,kIStStableFinalState,mother1,mother2); + evrec->AddParticle(*p_outD); + delete p_outD; + } //Add result without the remnant nucleus + } else { + ListeOfINCLresult.at(it).ARem[0]=A_i-theA_Remn-Aft; + ListeOfINCLresult.at(it).ZRem[0]=Z_i-theZ_Remn-Zl; + ListeOfINCLresult.at(it).pxRem[0]= the_pxRemn + (pxRemn*1000); + ListeOfINCLresult.at(it).pyRem[0]= the_pyRemn + (pyRemn*1000); + ListeOfINCLresult.at(it).pzRem[0]= the_pzRemn + (1000*pzRemn); + ListeOfINCLresult.at(it).EStarRem[0]=ExcitaionE; + if ( theDeExcitation != 0 ) theDeExcitation->deExcite(&ListeOfINCLresult.at(it)); + int rem_index=FindlargestFragment(ListeOfINCLresult.at(it)); + for (int nP=0; nP < ListeOfINCLresult.at(it).nParticles; nP++ ) { + A_f+=ListeOfINCLresult.at(it).A[nP]; + Z_f+=ListeOfINCLresult.at(it).Z[nP]; + if(nP==rem_index||ListeOfINCLresult.at(it).A[nP]>1) { + mother1=inucl; + mother2=-1; + } + else{ + mother1=Postion_evrec.at(0); + if(ListeOfINCLresult.size()==1) mother2=-1; + else if(ListeOfINCLresult.size()>1){ + mother2=Postion_evrec.at(ListeOfINCLresult.size()-1); + } + } + GHepParticle *p_outR = INCLtoGenieParticle(ListeOfINCLresult.at(it), + nP,kIStStableFinalState,mother1,mother2); + evrec->AddParticle(*p_outR); + has_remnant=true; + delete p_outR; } //Add all the result with the correct remnant nucleus } } @@ -771,10 +880,31 @@ void HINCLCascadeIntranuke::TransportHadrons(GHepRecord * evrec) const { } if ( ListeOfINCLresult.size() == 0 && !has_remnant) { - GHepParticle remnant(pdg_codeTargetRemn, kIStFinalStateNuclearRemnant, inucl,-1,-1,-1, fRemnP4, x4null); + TParticlePDG * pdgRemn=PDGLibrary::Instance()->Find(pdg_codeTargetRemn,false); + if(!pdgRemn) + { + LOG("HINCLCascadeIntranuke", pINFO) + << "NO Particle with pdg = " << pdg_codeTargetRemn << " in PDGLibrary!"; + GHepParticle remnant(kPdgHadronicBlob, kIStFinalStateNuclearRemnant, inucl,-1,-1,-1, fRemnP4, x4null); + evrec->AddParticle(remnant); + }else{ + GHepParticle remnant(pdg_codeTargetRemn, kIStStableFinalState, inucl,-1,-1,-1, fRemnP4, x4null); evrec->AddParticle(remnant); } - evrec->Particle(inucl)->SetStatus(kIStIntermediateState); +} +evrec->Particle(inucl)->SetStatus(kIStIntermediateState); + + +int dau1(0), dau2(0),Nsec=ListeOfINCLresult.size(); +if(Nsec>1){ + GHepParticle * pinN = evrec->Particle(Postion_evrec.at(0)); + dau1=pinN->FirstDaughter(); + dau2=pinN->LastDaughter(); + for(int ii=1;iiParticle(Postion_evrec.at(ii))->SetFirstDaughter(dau1); + evrec->Particle(Postion_evrec.at(ii))->SetLastDaughter(dau2); + } +} } //______________________________________________________________________________ @@ -810,8 +940,8 @@ bool HINCLCascadeIntranuke::NeedsRescattering(const GHepParticle * p) const { void HINCLCascadeIntranuke::Configure(const Registry & config) { LOG("HINCLCascadeIntranuke", pDEBUG) - << "Configure from Registry: '" << config.Name() << "'\n" - << config; + << "Configure from Registry: '" << config.Name() << "'\n" + << config; Algorithm::Configure(config); this->LoadConfig(); @@ -822,7 +952,7 @@ void HINCLCascadeIntranuke::Configure(const Registry & config) { void HINCLCascadeIntranuke::Configure(string param_set) { LOG("HINCLCascadeIntranuke", pDEBUG) - << "Configure from param_set name: " << param_set; + << "Configure from param_set name: " << param_set; Algorithm::Configure(param_set); this->LoadConfig(); diff --git a/src/Physics/HadronTransport/HNIntranuke2018.cxx b/src/Physics/HadronTransport/HNIntranuke2018.cxx index 269f8fc489..09f5eb144b 100644 --- a/src/Physics/HadronTransport/HNIntranuke2018.cxx +++ b/src/Physics/HadronTransport/HNIntranuke2018.cxx @@ -720,7 +720,7 @@ void HNIntranuke2018::ElasHN( t->SetPdgCode(tcode); double Mt = t->Mass(); //t->SetMomentum(TLorentzVector(0,0,0,Mt)); - + t->SetRemovalEnergy(0); // handle fermi momentum if(fDoFermi) { @@ -773,6 +773,9 @@ void HNIntranuke2018::InelasticHN(GHepRecord* ev, GHepParticle* p) const GHepParticle s1(*p); GHepParticle s2(*p); GHepParticle s3(*p); + s2.SetRemovalEnergy(0); + s3.SetRemovalEnergy(0); + if (utils::intranuke2018::PionProduction(ev,p,&s1,&s2,&s3,fRemnA,fRemnZ,fRemnP4,fDoFermi,fFermiFac,fFermiMomentum,fNuclmodel)) diff --git a/src/Physics/HadronTransport/INCLConvertParticle.hh b/src/Physics/HadronTransport/INCLConvertParticle.hh index 206d2f4569..2c49bcbe5e 100644 --- a/src/Physics/HadronTransport/INCLConvertParticle.hh +++ b/src/Physics/HadronTransport/INCLConvertParticle.hh @@ -40,40 +40,58 @@ namespace G4INCL { } } - G4INCL::ParticleType toINCLparticletype(int pdgc) { + int FindlargestFragment(G4INCL::EventInfo result){ + int maxA(0), rem_index(0); + for (size_t ll = 0; ll < result.nParticles; ll++) { + if (result.A[ll] > maxA) { + maxA = result.A[ll]; + rem_index = ll; + } + } + return rem_index; + } - if (pdgc == 2212) return G4INCL::Proton; - else if (pdgc == 2112) return G4INCL::Neutron; - else if (pdgc == 211) return G4INCL::PiPlus; - else if (pdgc == -211) return G4INCL::PiMinus; - else if (pdgc == 111) return G4INCL::PiZero; - else if (pdgc == 1000010020) return G4INCL::Composite; - else if (pdgc == 1000010030) return G4INCL::Composite; - else if (pdgc == 1000020030) return G4INCL::Composite; - else if (pdgc == 1000020040) return G4INCL::Composite; - else return G4INCL::UnknownParticle; + G4INCL::ParticleType toINCLparticletype(int pdgc) { - } + if (pdgc == 2212) return G4INCL::Proton; + else if (pdgc == 2112) return G4INCL::Neutron; + else if (pdgc == 211) return G4INCL::PiPlus; + else if (pdgc == -211) return G4INCL::PiMinus; + else if (pdgc == 111) return G4INCL::PiZero; + else if (pdgc == 1000010020) return G4INCL::Composite; + else if (pdgc == 1000010030) return G4INCL::Composite; + else if (pdgc == 1000020030) return G4INCL::Composite; + else if (pdgc == 1000020040) return G4INCL::Composite; + else return G4INCL::UnknownParticle; + + } - GHepParticle *INCLtoGenieParticle(G4INCL::EventInfo result, - int nP, GHepStatus_t ist, int mom1, int mom2) { - int pdg_codeP(0); - double E_pnP(0), EKinP = result.EKin[nP]; - TVector3 p3M(result.px[nP]/1000,result.py[nP]/1000,result.pz[nP]/1000); - TLorentzVector x4null(0.,0.,0.,0.); - double m_pnP = 0.5*((result.px[nP])*(result.px[nP]) + - (result.py[nP])*(result.py[nP]) + - (result.pz[nP])*(result.pz[nP]) - - EKinP*EKinP) / (EKinP); + GHepParticle *INCLtoGenieParticle(G4INCL::EventInfo result, + int nP, GHepStatus_t ist, int mom1, int mom2) { + int pdg_codeP(0); + double E_pnP(0), EKinP = result.EKin[nP]; + TVector3 p3M(result.px[nP]/1000,result.py[nP]/1000,result.pz[nP]/1000); + TLorentzVector x4null(0.,0.,0.,0.); + double m_pnP = 0.5*((result.px[nP])*(result.px[nP]) + + (result.py[nP])*(result.py[nP]) + + (result.pz[nP])*(result.pz[nP]) - + EKinP*EKinP) / (EKinP); if (m_pnP < 10 && result.A[nP] == 0 && result.Z[nP] == 0) { // photon pdg_codeP = 22; E_pnP = TMath::Sqrt((result.px[nP])*(result.px[nP]) + - (result.py[nP])*(result.py[nP]) + - (result.pz[nP])*(result.pz[nP]) ); + (result.py[nP])*(result.py[nP]) + + (result.pz[nP])*(result.pz[nP]) ); } else { pdg_codeP = INCLtopdgcode(result.A[nP],result.Z[nP]); - double Mass_prodPar = PDGLibrary::Instance()->Find(pdg_codeP)->Mass(); - E_pnP = EKinP + Mass_prodPar*1000; + TParticlePDG * pdgRemn=PDGLibrary::Instance()->Find(pdg_codeP,false); + if(!pdgRemn) + { + LOG("HINCLCascadeIntranuke", pINFO) + << "NO Particle with pdg = " << pdg_codeP << " in PDGLibrary!"; + pdg_codeP=kPdgHadronicBlob; + ist=kIStFinalStateNuclearRemnant; + } + E_pnP = EKinP + m_pnP; } TLorentzVector p4tgtf(p3M,E_pnP/1000); return new GHepParticle(pdg_codeP,ist,mom1,mom2,-1,-1,p4tgtf,x4null); diff --git a/src/Physics/HadronTransport/NucBindEnergyAggregator.cxx b/src/Physics/HadronTransport/NucBindEnergyAggregator.cxx index 8d67821d0d..1987a9ff70 100644 --- a/src/Physics/HadronTransport/NucBindEnergyAggregator.cxx +++ b/src/Physics/HadronTransport/NucBindEnergyAggregator.cxx @@ -88,7 +88,6 @@ void NucBindEnergyAggregator::ProcessEventRecord(GHepRecord * evrec) const LOG("Nuclear", pINFO) << "Kinetic energy before subtraction = " << KE; KE -= bindE; - KE = TMath::Max(0.,KE); LOG("Nuclear", pINFO) << "Kinetic energy after subtraction = " << KE; @@ -114,6 +113,7 @@ void NucBindEnergyAggregator::ProcessEventRecord(GHepRecord * evrec) const p->SetPx ( pxn ); p->SetPy ( pyn ); p->SetPz ( pzn ); + p->SetRemovalEnergy(0); //-- and add a GHEP entry to record this in the event record and // conserve energy/momentum @@ -128,38 +128,50 @@ void NucBindEnergyAggregator::ProcessEventRecord(GHepRecord * evrec) const LOG("Nuclear", pNOTICE) << "Recombining remnant nucleus + f/s nucleon"; - LOG("Nuclear", pERROR) - << "*** This functionality is temporarily disabled"; -/* - LOG("Nuclear", pNOTICE) << *evrec; - // find the remnant nucleus + // find the remnant nucleus if corresponding cascade modules provide it int rnucpos = evrec->RemnantNucleusPosition(); - assert(rnucpos); - - GHepParticle * rnucl = evrec->Particle(rnucpos); - - // mark both the remnant nucleus and the final state nucleon as - // intermediate states - rnucl -> SetStatus(kIStIntermediateState); - p -> SetStatus(kIStIntermediateState); - - // figure out the recombined nucleus PDG code - int Z = rnucl->Z(); - int A = rnucl->A(); - if(pdg::IsProton(p->Pdg())) Z++; - A++; - int ipdgc = pdg::IonPdgCode(A,Z); - - // add-up their 4-momenta - double pxnuc = rnucl->Px() + p->Px(); - double pynuc = rnucl->Py() + p->Py(); - double pznuc = rnucl->Pz() + p->Pz(); - double Enuc = rnucl->E() + p->E(); - - evrec->AddParticle(ipdgc, kIStStableFinalState, - rnucpos,-1,-1,-1, pxnuc,pynuc,pznuc,Enuc, 0,0,0,0); -*/ + + if (rnucpos != -1) + { + GHepParticle * rnucl = evrec->Particle(rnucpos); + + // mark both the remnant nucleus and the final state nucleon as + // intermediate states + rnucl -> SetStatus(kIStIntermediateState); + p -> SetStatus(kIStIntermediateState); + + // figure out the recombined nucleus PDG code + int Z = rnucl->Z(); + int A = rnucl->A(); + if(pdg::IsProton(p->Pdg())) Z++; + A++; + int ipdgc = pdg::IonPdgCode(A,Z); + + // add-up their 4-momenta + double pxnuc = rnucl->Px() + p->Px(); + double pynuc = rnucl->Py() + p->Py(); + double pznuc = rnucl->Pz() + p->Pz(); + double Enuc = rnucl->E() + p->E(); + + evrec->AddParticle(ipdgc, kIStStableFinalState, + rnucpos,-1,-1,-1, pxnuc,pynuc,pznuc,Enuc, 0,0,0,0); + } + else + { + // Find hadronic blob which is provided by Intranuke models + GHepParticle *hblob = evrec->FindParticle(kPdgHadronicBlob , kIStFinalStateNuclearRemnant, 0); + // If there is no hadronic blob do nothing + if (hblob == nullptr) return; + p -> SetStatus(kIStIntermediateState); + + // add-up their 4-momenta + double pxhblob = hblob->Px() + p->Px(); + double pyhblob = hblob->Py() + p->Py(); + double pzhblob = hblob->Pz() + p->Pz(); + double Ehblob = hblob->E() + p->E(); + hblob->SetMomentum(pxhblob, pyhblob, pzhblob, Ehblob); + } } } } diff --git a/src/Physics/QuasiElastic/XSection/LinkDef.h b/src/Physics/QuasiElastic/XSection/LinkDef.h index 10eb46b728..b2ad5c3c49 100644 --- a/src/Physics/QuasiElastic/XSection/LinkDef.h +++ b/src/Physics/QuasiElastic/XSection/LinkDef.h @@ -46,6 +46,6 @@ // Wrappers for GSL/MathMore lib -#pragma link C++ class genie::utils::gsl::d2Xsec_dQ2dv::d2Xsec_dQ2dv; +#pragma link C++ class genie::utils::gsl::d2Xsec_dQ2dv; #endif diff --git a/src/Physics/Resonance/EventGen/LinkDef.h b/src/Physics/Resonance/EventGen/LinkDef.h index 6d2e3bbcab..541dda71ac 100644 --- a/src/Physics/Resonance/EventGen/LinkDef.h +++ b/src/Physics/Resonance/EventGen/LinkDef.h @@ -17,6 +17,6 @@ #pragma link C++ class genie::RSPPInteractionListGenerator; // Wrappers for GSL/MathMore lib -#pragma link C++ class genie::utils::gsl::d4XSecMK_dWQ2CosThetaPhi_E::d4XSecMK_dWQ2CosThetaPhi_E; +#pragma link C++ class genie::utils::gsl::d4XSecMK_dWQ2CosThetaPhi_E; #endif