| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // neutron_hp -- source file 26 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996 27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron trans 28 // A prototype of the low energy neutron transport model.
29 // 29 //
30 // 080801 Protect div0 error, when theCompundF 30 // 080801 Protect div0 error, when theCompundFraction is 1 by T. Koi
31 // 31 //
32 // P. Arce, June-2014 Conversion neutron_hp to 32 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
33 // 33 //
34 // June-2019 - E. Mendoza --> perform some cor 34 // June-2019 - E. Mendoza --> perform some corrections
35 35
36 #include "G4ParticleHPKallbachMannSyst.hh" 36 #include "G4ParticleHPKallbachMannSyst.hh"
37 37
38 #include "G4Exp.hh" 38 #include "G4Exp.hh"
39 #include "G4HadronicException.hh" 39 #include "G4HadronicException.hh"
40 #include "G4Log.hh" 40 #include "G4Log.hh"
41 #include "G4Pow.hh" 41 #include "G4Pow.hh"
42 #include "G4SystemOfUnits.hh" 42 #include "G4SystemOfUnits.hh"
43 #include "Randomize.hh" 43 #include "Randomize.hh"
44 44
45 G4double G4ParticleHPKallbachMannSyst::Sample( 45 G4double G4ParticleHPKallbachMannSyst::Sample(G4double anEnergy)
46 { 46 {
47 G4double result = 0.; 47 G4double result = 0.;
48 48
49 G4double zero = GetKallbachZero(anEnergy); 49 G4double zero = GetKallbachZero(anEnergy);
50 if (zero > 1) zero = 1.; 50 if (zero > 1) zero = 1.;
51 if (zero < -1) zero = -1.; 51 if (zero < -1) zero = -1.;
52 G4double max = Kallbach(zero, anEnergy); 52 G4double max = Kallbach(zero, anEnergy);
53 G4double upper = Kallbach(1., anEnergy); 53 G4double upper = Kallbach(1., anEnergy);
54 G4double lower = Kallbach(-1., anEnergy); 54 G4double lower = Kallbach(-1., anEnergy);
55 if (upper > max) max = upper; 55 if (upper > max) max = upper;
56 if (lower > max) max = lower; 56 if (lower > max) max = lower;
57 G4double value, random; 57 G4double value, random;
58 58
59 G4int icounter = 0; 59 G4int icounter = 0;
60 G4int icounter_max = 1024; 60 G4int icounter_max = 1024;
61 do { 61 do {
62 icounter++; 62 icounter++;
63 if (icounter > icounter_max) { 63 if (icounter > icounter_max) {
64 G4cout << "Loop-counter exceeded the thr 64 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of "
65 << __FILE__ << "." << G4endl; 65 << __FILE__ << "." << G4endl;
66 break; 66 break;
67 } 67 }
68 result = 2. * G4UniformRand() - 1; 68 result = 2. * G4UniformRand() - 1;
69 value = Kallbach(result, anEnergy) / max; 69 value = Kallbach(result, anEnergy) / max;
70 random = G4UniformRand(); 70 random = G4UniformRand();
71 } while (random > value); // Loop checking, 71 } while (random > value); // Loop checking, 11.05.2015, T. Koi
72 72
73 return result; 73 return result;
74 } 74 }
75 75
76 G4double G4ParticleHPKallbachMannSyst::Kallbac 76 G4double G4ParticleHPKallbachMannSyst::Kallbach(G4double cosTh, G4double anEnergy)
77 { 77 {
78 // Kallbach-Mann systematics without normali 78 // Kallbach-Mann systematics without normalization.
79 G4double result; 79 G4double result;
80 G4double theX = A(anEnergy) * cosTh; 80 G4double theX = A(anEnergy) * cosTh;
81 result = 81 result =
82 0.5 * (G4Exp(theX) * (1 + theCompoundFract 82 0.5 * (G4Exp(theX) * (1 + theCompoundFraction) + G4Exp(-theX) * (1 - theCompoundFraction));
83 return result; 83 return result;
84 } 84 }
85 85
86 G4double G4ParticleHPKallbachMannSyst::GetKall 86 G4double G4ParticleHPKallbachMannSyst::GetKallbachZero(G4double anEnergy)
87 { 87 {
88 G4double result; 88 G4double result;
89 // delta 2.0e-16 in not good. 89 // delta 2.0e-16 in not good.
90 // delta 4.0e-16 is OK 90 // delta 4.0e-16 is OK
91 // safety factor of 2 91 // safety factor of 2
92 G4double delta = 8.0e-16; 92 G4double delta = 8.0e-16;
93 if (std::abs(theCompoundFraction - 1) < delt 93 if (std::abs(theCompoundFraction - 1) < delta) {
94 theCompoundFraction = 1.0 - delta; 94 theCompoundFraction = 1.0 - delta;
95 } 95 }
96 result = 0.5 * (1. / A(anEnergy)) * G4Log((1 96 result = 0.5 * (1. / A(anEnergy)) * G4Log((1 - theCompoundFraction) / (1 + theCompoundFraction));
97 return result; 97 return result;
98 } 98 }
99 99
100 G4double G4ParticleHPKallbachMannSyst::A(G4dou 100 G4double G4ParticleHPKallbachMannSyst::A(G4double anEnergy)
101 { 101 {
102 G4double result; 102 G4double result;
103 G4double C1 = 0.04 / MeV; 103 G4double C1 = 0.04 / MeV;
104 G4double C2 = 1.8E-6 / (MeV * MeV * MeV); 104 G4double C2 = 1.8E-6 / (MeV * MeV * MeV);
105 G4double C3 = 6.7E-7 / (MeV * MeV * MeV * Me 105 G4double C3 = 6.7E-7 / (MeV * MeV * MeV * MeV);
106 106
107 G4double epsa = anEnergy * theTargetMass / ( 107 G4double epsa = anEnergy * theTargetMass / (theTargetMass + theIncidentMass);
108 G4int Ac = theTargetA + theProjectileA; 108 G4int Ac = theTargetA + theProjectileA;
109 G4int Nc = Ac - theTargetZ - theProjectileZ; 109 G4int Nc = Ac - theTargetZ - theProjectileZ;
110 G4int AA = theTargetA; 110 G4int AA = theTargetA;
111 G4int ZA = theTargetZ; 111 G4int ZA = theTargetZ;
112 G4double ea = epsa + SeparationEnergy(Ac, Nc 112 G4double ea = epsa + SeparationEnergy(Ac, Nc, AA, ZA, theProjectileA, theProjectileZ);
113 G4double Et1 = 130 * MeV; 113 G4double Et1 = 130 * MeV;
114 G4double R1 = std::min(ea, Et1); 114 G4double R1 = std::min(ea, Et1);
115 // theProductEnergy is still in CMS!!! 115 // theProductEnergy is still in CMS!!!
116 G4double epsb = theProductEnergy * (theProdu 116 G4double epsb = theProductEnergy * (theProductMass + theResidualMass) / theResidualMass;
117 G4int AB = theResidualA; 117 G4int AB = theResidualA;
118 G4int ZB = theResidualZ; 118 G4int ZB = theResidualZ;
119 G4double eb = epsb + SeparationEnergy(Ac, Nc 119 G4double eb = epsb + SeparationEnergy(Ac, Nc, AB, ZB, theProductA, theProductZ);
120 G4double X1 = R1 * eb / ea; 120 G4double X1 = R1 * eb / ea;
121 G4double Et3 = 41 * MeV; 121 G4double Et3 = 41 * MeV;
122 G4double R3 = std::min(ea, Et3); 122 G4double R3 = std::min(ea, Et3);
123 G4double X3 = R3 * eb / ea; 123 G4double X3 = R3 * eb / ea;
124 124
125 G4double Ma = 1; 125 G4double Ma = 1;
126 G4double mb = 1; 126 G4double mb = 1;
127 if (theProjectileA == 1 || (theProjectileZ = 127 if (theProjectileA == 1 || (theProjectileZ == 1 && theProjectileA == 2)) {
128 Ma = 1; 128 Ma = 1;
129 } // neutron,proton,deuteron 129 } // neutron,proton,deuteron
130 else if (theProjectileA == 4 && theProjectil 130 else if (theProjectileA == 4 && theProjectileZ == 2) {
131 Ma = 0; 131 Ma = 0;
132 } // alpha 132 } // alpha
133 else if (theProjectileA == 3 && (theProjecti 133 else if (theProjectileA == 3 && (theProjectileZ == 1 || theProjectileZ == 2)) {
134 Ma = 0.5; 134 Ma = 0.5;
135 } // tritum,He3 : set intermediate value 135 } // tritum,He3 : set intermediate value
136 else { 136 else {
137 throw G4HadronicException(__FILE__, __LINE 137 throw G4HadronicException(__FILE__, __LINE__,
138 "Severe error in 138 "Severe error in the sampling of Kallbach-Mann Systematics");
139 } 139 }
140 if (theProductA == 1 && theProductZ == 0) { 140 if (theProductA == 1 && theProductZ == 0) {
141 mb = 1. / 2.; 141 mb = 1. / 2.;
142 } // neutron 142 } // neutron
143 else if (theProductA == 4 && theProductZ == 143 else if (theProductA == 4 && theProductZ == 2) {
144 mb = 2; 144 mb = 2;
145 } // alpha 145 } // alpha
146 else { 146 else {
147 mb = 1; 147 mb = 1;
148 } 148 }
149 149
150 result = C1 * X1 + C2 * G4Pow::GetInstance() 150 result = C1 * X1 + C2 * G4Pow::GetInstance()->powN(X1, 3)
151 + C3 * Ma * mb * G4Pow::GetInstance 151 + C3 * Ma * mb * G4Pow::GetInstance()->powN(X3, 4);
152 return result; 152 return result;
153 } 153 }
154 154
155 G4double G4ParticleHPKallbachMannSyst::Separat 155 G4double G4ParticleHPKallbachMannSyst::SeparationEnergy(G4int Ac, G4int Nc, G4int AA, G4int ZA,
156 156 G4int Abinding, G4int Zbinding)
157 { 157 {
158 G4double result; 158 G4double result;
159 G4int NA = AA - ZA; 159 G4int NA = AA - ZA;
160 G4int Zc = Ac - Nc; 160 G4int Zc = Ac - Nc;
161 result = 15.68 * (Ac - AA); 161 result = 15.68 * (Ac - AA);
162 result += -28.07 * ((Nc - Zc) * (Nc - Zc) / 162 result += -28.07 * ((Nc - Zc) * (Nc - Zc) / (G4double)Ac - (NA - ZA) * (NA - ZA) / (G4double)AA);
163 result += 163 result +=
164 -18.56 * (G4Pow::GetInstance()->A23(G4doub 164 -18.56 * (G4Pow::GetInstance()->A23(G4double(Ac)) - G4Pow::GetInstance()->A23(G4double(AA)));
165 result += 33.22 165 result += 33.22
166 * ((Nc - Zc) * (Nc - Zc) / G4Pow:: 166 * ((Nc - Zc) * (Nc - Zc) / G4Pow::GetInstance()->powA(G4double(Ac), 4. / 3.)
167 - (NA - ZA) * (NA - ZA) / G4Pow 167 - (NA - ZA) * (NA - ZA) / G4Pow::GetInstance()->powA(G4double(AA), 4. / 3.));
168 result += -0.717 168 result += -0.717
169 * (Zc * Zc / G4Pow::GetInstance()- 169 * (Zc * Zc / G4Pow::GetInstance()->A13(G4double(Ac))
170 - ZA * ZA / G4Pow::GetInstance( 170 - ZA * ZA / G4Pow::GetInstance()->A13(G4double(AA)));
171 result += 1.211 * (Zc * Zc / (G4double)Ac - 171 result += 1.211 * (Zc * Zc / (G4double)Ac - ZA * ZA / (G4double)AA);
172 G4double totalBinding(0); 172 G4double totalBinding(0);
173 if (Zbinding == 0 && Abinding == 1) totalBin 173 if (Zbinding == 0 && Abinding == 1) totalBinding = 0;
174 if (Zbinding == 1 && Abinding == 1) totalBin 174 if (Zbinding == 1 && Abinding == 1) totalBinding = 0;
175 if (Zbinding == 1 && Abinding == 2) totalBin 175 if (Zbinding == 1 && Abinding == 2) totalBinding = 2.224596;
176 if (Zbinding == 1 && Abinding == 3) totalBin 176 if (Zbinding == 1 && Abinding == 3) totalBinding = 8.481798;
177 if (Zbinding == 2 && Abinding == 3) totalBin 177 if (Zbinding == 2 && Abinding == 3) totalBinding = 7.718043;
178 if (Zbinding == 2 && Abinding == 4) totalBin 178 if (Zbinding == 2 && Abinding == 4) totalBinding = 28.29566;
179 result += -totalBinding; 179 result += -totalBinding;
180 result *= MeV; 180 result *= MeV;
181 return result; 181 return result;
182 } 182 }
183 183