| 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 #include "G4ios.hh" 26 #include "G4ios.hh"
27 #include "G4Clebsch.hh" 27 #include "G4Clebsch.hh"
28 #include "G4Pow.hh" 28 #include "G4Pow.hh"
29 #include "G4Exp.hh" 29 #include "G4Exp.hh"
30 #include "G4Log.hh" 30 #include "G4Log.hh"
31 #include "Randomize.hh" 31 #include "Randomize.hh"
32 32
33 const G4int G4POWLOGFACTMAX = 512; 33 const G4int G4POWLOGFACTMAX = 512;
34 34
35 using namespace std; 35 using namespace std;
36 36
37 G4double G4Clebsch::ClebschGordanCoeff(G4int t 37 G4double G4Clebsch::ClebschGordanCoeff(G4int twoJ1, G4int twoM1,
38 G4int t 38 G4int twoJ2, G4int twoM2,
39 G4int t 39 G4int twoJ)
40 { 40 {
41 if(twoJ1 < 0 || twoJ2 < 0 || twoJ < 0 || 41 if(twoJ1 < 0 || twoJ2 < 0 || twoJ < 0 ||
42 ((twoJ1-twoM1) % 2) || ((twoJ2-twoM2) % 2 42 ((twoJ1-twoM1) % 2) || ((twoJ2-twoM2) % 2)) { return 0; }
43 43
44 G4int twoM = twoM1 + twoM2; 44 G4int twoM = twoM1 + twoM2;
45 if(twoM1 > twoJ1 || twoM1 < -twoJ1 || 45 if(twoM1 > twoJ1 || twoM1 < -twoJ1 ||
46 twoM2 > twoJ2 || twoM2 < -twoJ2 || 46 twoM2 > twoJ2 || twoM2 < -twoJ2 ||
47 twoM > twoJ || twoM < -twoJ) { return 0; 47 twoM > twoJ || twoM < -twoJ) { return 0; }
48 48
49 // Checks limits on J1, J2, J3 49 // Checks limits on J1, J2, J3
50 G4double triangle = TriangleCoeff(twoJ1, two 50 G4double triangle = TriangleCoeff(twoJ1, twoJ2, twoJ);
51 if(triangle == 0) { return 0; } 51 if(triangle == 0) { return 0; }
52 52
53 G4Pow* g4pow = G4Pow::GetInstance(); 53 G4Pow* g4pow = G4Pow::GetInstance();
54 G4double factor = g4pow->logfactorial((twoJ1 54 G4double factor = g4pow->logfactorial((twoJ1 + twoM1)/2) +
55 g4pow->logfactorial((twoJ1 55 g4pow->logfactorial((twoJ1 - twoM1)/2);
56 factor += g4pow->logfactorial((twoJ2 + twoM2 56 factor += g4pow->logfactorial((twoJ2 + twoM2)/2) +
57 g4pow->logfactorial((twoJ2 - twoM2 57 g4pow->logfactorial((twoJ2 - twoM2)/2);
58 factor += g4pow->logfactorial((twoJ + twoM)/ 58 factor += g4pow->logfactorial((twoJ + twoM)/2) +
59 g4pow->logfactorial((twoJ - twoM)/ 59 g4pow->logfactorial((twoJ - twoM)/2);
60 factor *= 0.5; 60 factor *= 0.5;
61 61
62 G4int kMin = 0; 62 G4int kMin = 0;
63 G4int sum1 = (twoJ1 - twoM1)/2; 63 G4int sum1 = (twoJ1 - twoM1)/2;
64 G4int kMax = sum1; 64 G4int kMax = sum1;
65 G4int sum2 = (twoJ - twoJ2 + twoM1)/2; 65 G4int sum2 = (twoJ - twoJ2 + twoM1)/2;
66 if(-sum2 > kMin) kMin = -sum2; 66 if(-sum2 > kMin) kMin = -sum2;
67 G4int sum3 = (twoJ2 + twoM2)/2; 67 G4int sum3 = (twoJ2 + twoM2)/2;
68 if(sum3 < kMax) kMax = sum3; 68 if(sum3 < kMax) kMax = sum3;
69 G4int sum4 = (twoJ - twoJ1 - twoM2)/2; 69 G4int sum4 = (twoJ - twoJ1 - twoM2)/2;
70 if(-sum4 > kMin) kMin = -sum4; 70 if(-sum4 > kMin) kMin = -sum4;
71 G4int sum5 = (twoJ1 + twoJ2 - twoJ)/2; 71 G4int sum5 = (twoJ1 + twoJ2 - twoJ)/2;
72 if(sum5 < kMax) kMax = sum5; 72 if(sum5 < kMax) kMax = sum5;
73 73
74 // sanity / boundary checks 74 // sanity / boundary checks
75 if(kMin < 0) { 75 if(kMin < 0) {
76 G4Exception("G4Clebsch::ClebschGordanCoeff 76 G4Exception("G4Clebsch::ClebschGordanCoeff()", "Clebsch001",
77 JustWarning, "kMin < 0"); 77 JustWarning, "kMin < 0");
78 return 0; 78 return 0;
79 } 79 }
80 if(kMax < kMin) { 80 if(kMax < kMin) {
81 G4Exception("G4Clebsch::ClebschGordanCoeff 81 G4Exception("G4Clebsch::ClebschGordanCoeff()", "Clebsch002",
82 JustWarning, "kMax < kMin"); 82 JustWarning, "kMax < kMin");
83 return 0; 83 return 0;
84 } 84 }
85 if(kMax >= G4POWLOGFACTMAX) { 85 if(kMax >= G4POWLOGFACTMAX) {
86 G4Exception("G4Clebsch::ClebschGordanCoeff 86 G4Exception("G4Clebsch::ClebschGordanCoeff()", "Clebsch003",
87 JustWarning, "kMax too big for G4Pow"); 87 JustWarning, "kMax too big for G4Pow");
88 return 0; 88 return 0;
89 } 89 }
90 90
91 // Now do the sum over k 91 // Now do the sum over k
92 G4double kSum = 0.; 92 G4double kSum = 0.;
93 for(G4int k = kMin; k <= kMax; k++) { 93 for(G4int k = kMin; k <= kMax; k++) {
94 G4double sign = (k % 2) ? -1 : 1; 94 G4double sign = (k % 2) ? -1 : 1;
95 kSum += sign * G4Exp(factor - g4pow->logfa 95 kSum += sign * G4Exp(factor - g4pow->logfactorial(sum1-k) -
96 g4pow->logfactorial(sum2+k) - 96 g4pow->logfactorial(sum2+k) -
97 g4pow->logfactorial(sum3-k) - 97 g4pow->logfactorial(sum3-k) -
98 g4pow->logfactorial(sum4+k) - 98 g4pow->logfactorial(sum4+k) -
99 g4pow->logfactorial(k) - 99 g4pow->logfactorial(k) -
100 g4pow->logfactorial(sum5-k)); 100 g4pow->logfactorial(sum5-k));
101 } 101 }
102 102
103 return triangle*sqrt(twoJ+1)*kSum; 103 return triangle*sqrt(twoJ+1)*kSum;
104 } 104 }
105 105
106 G4double G4Clebsch::ClebschGordan(G4int twoJ1, 106 G4double G4Clebsch::ClebschGordan(G4int twoJ1, G4int twoM1,
107 G4int twoJ2, G4int twoM2, 107 G4int twoJ2, G4int twoM2,
108 G4int twoJ) 108 G4int twoJ)
109 { 109 {
110 // ClebschGordanCoeff() will do all input ch 110 // ClebschGordanCoeff() will do all input checking
111 G4double clebsch = ClebschGordanCoeff(twoJ1, 111 G4double clebsch = ClebschGordanCoeff(twoJ1, twoM1, twoJ2, twoM2, twoJ);
112 return clebsch*clebsch; 112 return clebsch*clebsch;
113 } 113 }
114 114
115 std::vector<G4double> 115 std::vector<G4double>
116 G4Clebsch::GenerateIso3(G4int twoJ1, G4int two 116 G4Clebsch::GenerateIso3(G4int twoJ1, G4int twoM1,
117 G4int twoJ2, G4int twoM2, 117 G4int twoJ2, G4int twoM2,
118 G4int twoJOut1, G4int twoJOut2) 118 G4int twoJOut1, G4int twoJOut2)
119 { 119 {
120 std::vector<G4double> temp; 120 std::vector<G4double> temp;
121 121
122 // ---- Special cases first ---- 122 // ---- Special cases first ----
123 123
124 // Special case, both Jin are zero 124 // Special case, both Jin are zero
125 if (twoJ1 == 0 && twoJ2 == 0) { 125 if (twoJ1 == 0 && twoJ2 == 0) {
126 G4Exception("G4Clebsch::GenerateIso3()", " 126 G4Exception("G4Clebsch::GenerateIso3()", "Clebsch010",
127 JustWarning, "both twoJ are zero"); 127 JustWarning, "both twoJ are zero");
128 temp.push_back(0.); 128 temp.push_back(0.);
129 temp.push_back(0.); 129 temp.push_back(0.);
130 return temp; 130 return temp;
131 } 131 }
132 132
133 G4int twoM3 = twoM1 + twoM2; 133 G4int twoM3 = twoM1 + twoM2;
134 134
135 // Special case, either Jout is zero 135 // Special case, either Jout is zero
136 if (twoJOut1 == 0) { 136 if (twoJOut1 == 0) {
137 temp.push_back(0.); 137 temp.push_back(0.);
138 temp.push_back(twoM3); 138 temp.push_back(twoM3);
139 return temp; 139 return temp;
140 } 140 }
141 if (twoJOut2 == 0) { 141 if (twoJOut2 == 0) {
142 temp.push_back(twoM3); 142 temp.push_back(twoM3);
143 temp.push_back(0.); 143 temp.push_back(0.);
144 return temp; 144 return temp;
145 } 145 }
146 146
147 // Number of possible states, in 147 // Number of possible states, in
148 G4int twoJMinIn = std::max(std::abs(twoJ1 - 148 G4int twoJMinIn = std::max(std::abs(twoJ1 - twoJ2), std::abs(twoM3));
149 G4int twoJMaxIn = twoJ1 + twoJ2; 149 G4int twoJMaxIn = twoJ1 + twoJ2;
150 150
151 // Number of possible states, out 151 // Number of possible states, out
152 G4int twoJMinOut = 9999; 152 G4int twoJMinOut = 9999;
153 for(G4int i=-1; i<=1; i+=2) { 153 for(G4int i=-1; i<=1; i+=2) {
154 for(G4int j=-1; j<=1; j+=2) { 154 for(G4int j=-1; j<=1; j+=2) {
155 G4int twoJTmp= std::abs(i*twoJOut1 + j*t 155 G4int twoJTmp= std::abs(i*twoJOut1 + j*twoJOut2);
156 if(twoJTmp < twoJMinOut) twoJMinOut = tw 156 if(twoJTmp < twoJMinOut) twoJMinOut = twoJTmp;
157 } 157 }
158 } 158 }
159 twoJMinOut = std::max(twoJMinOut, std::abs(t 159 twoJMinOut = std::max(twoJMinOut, std::abs(twoM3));
160 G4int twoJMaxOut = twoJOut1 + twoJOut2; 160 G4int twoJMaxOut = twoJOut1 + twoJOut2;
161 161
162 // Possible in and out common states 162 // Possible in and out common states
163 G4int twoJMin = std::max(twoJMinIn, twoJMi 163 G4int twoJMin = std::max(twoJMinIn, twoJMinOut);
164 G4int twoJMax = std::min(twoJMaxIn, twoJMa 164 G4int twoJMax = std::min(twoJMaxIn, twoJMaxOut);
165 if (twoJMin > twoJMax) { 165 if (twoJMin > twoJMax) {
166 G4Exception("G4Clebsch::GenerateIso3()", " 166 G4Exception("G4Clebsch::GenerateIso3()", "Clebsch020",
167 JustWarning, "twoJMin > twoJMax"); 167 JustWarning, "twoJMin > twoJMax");
168 return temp; 168 return temp;
169 } 169 }
170 170
171 // Number of possible isospins 171 // Number of possible isospins
172 G4int nJ = (twoJMax - twoJMin) / 2 + 1; 172 G4int nJ = (twoJMax - twoJMin) / 2 + 1;
173 173
174 // A few consistency checks 174 // A few consistency checks
175 175
176 if ( (twoJ1 == 0 || twoJ2 == 0) && twoJMin ! 176 if ( (twoJ1 == 0 || twoJ2 == 0) && twoJMin != twoJMax ) {
177 G4Exception("G4Clebsch::GenerateIso3()", " 177 G4Exception("G4Clebsch::GenerateIso3()", "Clebsch021",
178 JustWarning, "twoJ1 or twoJ2 = 0, but twoJ 178 JustWarning, "twoJ1 or twoJ2 = 0, but twoJMin != JMax");
179 return temp; 179 return temp;
180 } 180 }
181 181
182 // MGP ---- Shall it be a warning or an exce 182 // MGP ---- Shall it be a warning or an exception?
183 if (nJ == 0) { 183 if (nJ == 0) {
184 G4Exception("G4Clebsch::GenerateIso3()", " 184 G4Exception("G4Clebsch::GenerateIso3()", "Clebsch022",
185 JustWarning, "nJ is zero, no overlap betwe 185 JustWarning, "nJ is zero, no overlap between in and out");
186 return temp; 186 return temp;
187 } 187 }
188 188
189 // Loop over all possible combinations of tw 189 // Loop over all possible combinations of twoJ1, twoJ2, twoM11, twoM2, twoJTot
190 // to get the probability of each of the in- 190 // to get the probability of each of the in-channel couplings
191 191
192 std::vector<G4double> clebsch; 192 std::vector<G4double> clebsch;
193 G4double sum = 0.0; 193 G4double sum = 0.0;
194 for(G4int twoJ=twoJMin; twoJ<=twoJMax; twoJ+ 194 for(G4int twoJ=twoJMin; twoJ<=twoJMax; twoJ+=2) {
195 sum += ClebschGordan(twoJ1, twoM1, twoJ2, 195 sum += ClebschGordan(twoJ1, twoM1, twoJ2, twoM2, twoJ);
196 clebsch.push_back(sum); 196 clebsch.push_back(sum);
197 } 197 }
198 198
199 // Consistency check 199 // Consistency check
200 if (static_cast<G4int>(clebsch.size()) != nJ 200 if (static_cast<G4int>(clebsch.size()) != nJ) {
201 G4Exception("G4Clebsch::GenerateIso3()", " 201 G4Exception("G4Clebsch::GenerateIso3()", "Clebsch023",
202 JustWarning, "nJ inconsistency"); 202 JustWarning, "nJ inconsistency");
203 return temp; 203 return temp;
204 } 204 }
205 205
206 // Consistency check 206 // Consistency check
207 if (sum <= 0.) { 207 if (sum <= 0.) {
208 G4Exception("G4Clebsch::GenerateIso3()", " 208 G4Exception("G4Clebsch::GenerateIso3()", "Clebsch024",
209 JustWarning, "Sum of Clebsch-Gordan probab 209 JustWarning, "Sum of Clebsch-Gordan probabilities <=0");
210 return temp; 210 return temp;
211 } 211 }
212 212
213 // Generate a random twoJTot according to th 213 // Generate a random twoJTot according to the Clebsch-Gordan pdf
214 sum *= G4UniformRand(); 214 sum *= G4UniformRand();
215 G4int twoJTot = twoJMin; 215 G4int twoJTot = twoJMin;
216 for (G4int i=0; i<nJ; ++i) { 216 for (G4int i=0; i<nJ; ++i) {
217 if (sum < clebsch[i]) { 217 if (sum < clebsch[i]) {
218 twoJTot += 2*i; 218 twoJTot += 2*i;
219 break; 219 break;
220 } 220 }
221 } 221 }
222 222
223 // Generate twoM3Out 223 // Generate twoM3Out
224 224
225 std::vector<G4double> mMin; 225 std::vector<G4double> mMin;
226 mMin.push_back(-twoJOut1); 226 mMin.push_back(-twoJOut1);
227 mMin.push_back(-twoJOut2); 227 mMin.push_back(-twoJOut2);
228 228
229 std::vector<G4double> mMax; 229 std::vector<G4double> mMax;
230 mMax.push_back(twoJOut1); 230 mMax.push_back(twoJOut1);
231 mMax.push_back(twoJOut2); 231 mMax.push_back(twoJOut2);
232 232
233 // Calculate the possible |J_i M_i> combinat 233 // Calculate the possible |J_i M_i> combinations and their probability
234 234
235 std::vector<G4double> m1Out; 235 std::vector<G4double> m1Out;
236 std::vector<G4double> m2Out; 236 std::vector<G4double> m2Out;
237 237
238 const G4int size = 20; 238 const G4int size = 20;
239 G4double prbout[size][size]; 239 G4double prbout[size][size];
240 240
241 G4int m1pos(0), m2pos(0); 241 G4int m1pos(0), m2pos(0);
242 G4int j12; 242 G4int j12;
243 G4int m1pr(0), m2pr(0); 243 G4int m1pr(0), m2pr(0);
244 244
245 sum = 0.; 245 sum = 0.;
246 for(j12 = std::abs(twoJOut1-twoJOut2); j12<= 246 for(j12 = std::abs(twoJOut1-twoJOut2); j12<=(twoJOut1+twoJOut2); j12+=2)
247 { 247 {
248 m1pos = -1; 248 m1pos = -1;
249 for (m1pr = static_cast<G4int>(mMin[0]+.00 249 for (m1pr = static_cast<G4int>(mMin[0]+.00001); m1pr <= mMax[0]; m1pr+=2)
250 { 250 {
251 m1pos++; 251 m1pos++;
252 if (m1pos >= size) { 252 if (m1pos >= size) {
253 G4Exception("G4Clebsch::GenerateIso3() 253 G4Exception("G4Clebsch::GenerateIso3()", "Clebsch025",
254 JustWarning, "m1pos > size"); 254 JustWarning, "m1pos > size");
255 return temp; 255 return temp;
256 } 256 }
257 m1Out.push_back(m1pr); 257 m1Out.push_back(m1pr);
258 m2pos = -1; 258 m2pos = -1;
259 for (m2pr = static_cast<G4int>(mMin[1]+. 259 for (m2pr = static_cast<G4int>(mMin[1]+.00001); m2pr <= mMax[1]; m2pr+=2)
260 { 260 {
261 m2pos++; 261 m2pos++;
262 if (m2pos >= size) 262 if (m2pos >= size)
263 { 263 {
264 G4Exception("G4Clebsch::GenerateIso3 264 G4Exception("G4Clebsch::GenerateIso3()", "Clebsch026",
265 JustWarning, "m2pos > size"); 265 JustWarning, "m2pos > size");
266 return temp; 266 return temp;
267 } 267 }
268 m2Out.push_back(m2pr); 268 m2Out.push_back(m2pr);
269 269
270 if(m1pr + m2pr == twoM3) 270 if(m1pr + m2pr == twoM3)
271 { 271 {
272 G4int m12 = m1pr + m2pr; 272 G4int m12 = m1pr + m2pr;
273 G4double c12 = ClebschGordan(twoJOut 273 G4double c12 = ClebschGordan(twoJOut1, m1pr, twoJOut2,m2pr, j12);
274 G4double c34 = ClebschGordan(0,0,0,0 274 G4double c34 = ClebschGordan(0,0,0,0,0);
275 G4double ctot = ClebschGordan(j12, m 275 G4double ctot = ClebschGordan(j12, m12, 0, 0, twoJTot);
276 G4double cleb = c12*c34*ctot; 276 G4double cleb = c12*c34*ctot;
277 prbout[m1pos][m2pos] = cleb; 277 prbout[m1pos][m2pos] = cleb;
278 sum += cleb; 278 sum += cleb;
279 } 279 }
280 else 280 else
281 { 281 {
282 prbout[m1pos][m2pos] = 0.; 282 prbout[m1pos][m2pos] = 0.;
283 } 283 }
284 } 284 }
285 } 285 }
286 } 286 }
287 287
288 if (sum <= 0.) { 288 if (sum <= 0.) {
289 G4Exception("G4Clebsch::GenerateIso3()", " 289 G4Exception("G4Clebsch::GenerateIso3()", "Clebsch027",
290 JustWarning, "sum (out) <=0"); 290 JustWarning, "sum (out) <=0");
291 return temp; 291 return temp;
292 } 292 }
293 293
294 for (G4int i=0; i<size; i++) { 294 for (G4int i=0; i<size; i++) {
295 for (G4int j=0; j<size; j++) { 295 for (G4int j=0; j<size; j++) {
296 prbout[i][j] /= sum; 296 prbout[i][j] /= sum;
297 } 297 }
298 } 298 }
299 299
300 G4double rand = G4UniformRand(); 300 G4double rand = G4UniformRand();
301 301
302 G4int m1p, m2p; 302 G4int m1p, m2p;
303 303
304 for (m1p=0; m1p<m1pos; m1p++) { 304 for (m1p=0; m1p<m1pos; m1p++) {
305 for (m2p=0; m2p<m2pos; m2p++) { 305 for (m2p=0; m2p<m2pos; m2p++) {
306 if (rand < prbout[m1p][m2p]) { 306 if (rand < prbout[m1p][m2p]) {
307 temp.push_back(m1Out[m1p]); 307 temp.push_back(m1Out[m1p]);
308 temp.push_back(m2Out[m2p]); 308 temp.push_back(m2Out[m2p]);
309 return temp; 309 return temp;
310 } 310 }
311 else rand -= prbout[m1p][m2p]; 311 else rand -= prbout[m1p][m2p];
312 } 312 }
313 } 313 }
314 314
315 G4Exception("G4Clebsch::GenerateIso3()", "Cl 315 G4Exception("G4Clebsch::GenerateIso3()", "Clebsch028",
316 JustWarning, "Should never get here"); 316 JustWarning, "Should never get here");
317 return temp; 317 return temp;
318 } 318 }
319 319
320 G4double G4Clebsch::Weight(G4int twoJ1, G4int 320 G4double G4Clebsch::Weight(G4int twoJ1, G4int twoM1,
321 G4int twoJ2, G4int twoM2, 321 G4int twoJ2, G4int twoM2,
322 G4int twoJOut1, G4int twoJOut2) 322 G4int twoJOut1, G4int twoJOut2)
323 { 323 {
324 G4double value = 0.; 324 G4double value = 0.;
325 325
326 G4int twoM = twoM1 + twoM2; 326 G4int twoM = twoM1 + twoM2;
327 327
328 G4int twoJMinIn = std::max(std::abs(twoJ1 - 328 G4int twoJMinIn = std::max(std::abs(twoJ1 - twoJ2), std::abs(twoM));
329 G4int twoJMaxIn = twoJ1 + twoJ2; 329 G4int twoJMaxIn = twoJ1 + twoJ2;
330 330
331 G4int twoJMinOut = std::max(std::abs(twoJOut 331 G4int twoJMinOut = std::max(std::abs(twoJOut1 - twoJOut2), std::abs(twoM));
332 G4int twoJMaxOut = twoJOut1 + twoJOut2; 332 G4int twoJMaxOut = twoJOut1 + twoJOut2;
333 333
334 G4int twoJMin = std::max(twoJMinIn,twoJMinOu 334 G4int twoJMin = std::max(twoJMinIn,twoJMinOut);
335 G4int twoJMax = std::min(twoJMaxIn,twoJMaxOu 335 G4int twoJMax = std::min(twoJMaxIn,twoJMaxOut);
336 336
337 for (G4int twoJ=twoJMin; twoJ<=twoJMax; twoJ 337 for (G4int twoJ=twoJMin; twoJ<=twoJMax; twoJ+=2) {
338 // ClebschGordan() will do all input check 338 // ClebschGordan() will do all input checking
339 value += ClebschGordan(twoJ1, twoM1, twoJ2 339 value += ClebschGordan(twoJ1, twoM1, twoJ2, twoM2, twoJ);
340 } 340 }
341 341
342 return value; 342 return value;
343 } 343 }
344 344
345 G4double G4Clebsch::Wigner3J(G4double j1, G4do 345 G4double G4Clebsch::Wigner3J(G4double j1, G4double j2, G4double j3,
346 G4double m1, G4double m2, G4double 346 G4double m1, G4double m2, G4double m3)
347 { 347 {
348 // G4Exception("G4Clebsch::Wigner3J()", "Cl 348 // G4Exception("G4Clebsch::Wigner3J()", "Clebsch030", JustWarning,
349 // "G4Clebsch::Wigner3J with double argumen 349 // "G4Clebsch::Wigner3J with double arguments is deprecated. Please use G4int version.");
350 G4int twoJ1 = (G4int) (2.*j1); 350 G4int twoJ1 = (G4int) (2.*j1);
351 G4int twoJ2 = (G4int) (2.*j2); 351 G4int twoJ2 = (G4int) (2.*j2);
352 G4int twoJ3 = (G4int) (2.*j3); 352 G4int twoJ3 = (G4int) (2.*j3);
353 G4int twoM1 = (G4int) (2.*m1); 353 G4int twoM1 = (G4int) (2.*m1);
354 G4int twoM2 = (G4int) (2.*m2); 354 G4int twoM2 = (G4int) (2.*m2);
355 G4int twoM3 = (G4int) (2.*m3); 355 G4int twoM3 = (G4int) (2.*m3);
356 return Wigner3J(twoJ1, twoM1, twoJ2, twoM2, 356 return Wigner3J(twoJ1, twoM1, twoJ2, twoM2, twoJ3, twoM3);
357 } 357 }
358 358
359 G4double G4Clebsch::Wigner3J(G4int twoJ1, G4in 359 G4double G4Clebsch::Wigner3J(G4int twoJ1, G4int twoM1,
360 G4int twoJ2, G4in 360 G4int twoJ2, G4int twoM2,
361 G4int twoJ3) 361 G4int twoJ3)
362 { 362 {
363 G4double clebsch = ClebschGordanCoeff(twoJ1, 363 G4double clebsch = ClebschGordanCoeff(twoJ1, twoM1, twoJ2, twoM2, twoJ3);
364 if(clebsch == 0) return clebsch; 364 if(clebsch == 0) return clebsch;
365 if( (twoJ1-twoJ2+twoM1+twoM2)/2 % 2) clebsch 365 if( (twoJ1-twoJ2+twoM1+twoM2)/2 % 2) clebsch = -clebsch;
366 return clebsch / sqrt(twoJ3+1); 366 return clebsch / sqrt(twoJ3+1);
367 } 367 }
368 368
369 G4double G4Clebsch::Wigner3J(G4int twoJ1, G4in 369 G4double G4Clebsch::Wigner3J(G4int twoJ1, G4int twoM1,
370 G4int twoJ2, G4in 370 G4int twoJ2, G4int twoM2,
371 G4int twoJ3, G4in 371 G4int twoJ3, G4int twoM3)
372 { 372 {
373 if(twoM1 + twoM2 != -twoM3) return 0; 373 if(twoM1 + twoM2 != -twoM3) return 0;
374 G4double clebsch = ClebschGordanCoeff(twoJ1, 374 G4double clebsch = ClebschGordanCoeff(twoJ1, twoM1, twoJ2, twoM2, twoJ3);
375 if(clebsch == 0) return clebsch; 375 if(clebsch == 0) return clebsch;
376 if( (twoJ1-twoJ2-twoM3)/2 % 2) clebsch = -cl 376 if( (twoJ1-twoJ2-twoM3)/2 % 2) clebsch = -clebsch;
377 return clebsch / sqrt(twoJ3+1); 377 return clebsch / sqrt(twoJ3+1);
378 } 378 }
379 379
380 G4double G4Clebsch::NormalizedClebschGordan(G4 380 G4double G4Clebsch::NormalizedClebschGordan(G4int twoJ, G4int twoM,
381 G4int twoJ1, G4int twoJ2, 381 G4int twoJ1, G4int twoJ2,
382 G4int twoM1, G4int twoM2) 382 G4int twoM1, G4int twoM2)
383 { 383 {
384 // Calculate the normalized Clebsch-Gordan c 384 // Calculate the normalized Clebsch-Gordan coefficient, that is the prob
385 // of isospin decomposition of (J,m) into J1 385 // of isospin decomposition of (J,m) into J1, J2, m1, m2
386 386
387 G4double cleb = 0.; 387 G4double cleb = 0.;
388 if(twoJ1 == 0 || twoJ2 == 0) return cleb; 388 if(twoJ1 == 0 || twoJ2 == 0) return cleb;
389 389
390 // Loop over all J1,J2,Jtot,m1,m2 combinatio 390 // Loop over all J1,J2,Jtot,m1,m2 combinations
391 G4double sum = 0.0; 391 G4double sum = 0.0;
392 for(G4int twoM1Current=-twoJ1; twoM1Current< 392 for(G4int twoM1Current=-twoJ1; twoM1Current<=twoJ1; twoM1Current+=2) {
393 G4int twoM2Current = twoM - twoM1Current; 393 G4int twoM2Current = twoM - twoM1Current;
394 // ClebschGordan() will do all further inp 394 // ClebschGordan() will do all further input checking
395 G4double prob = ClebschGordan(twoJ1, twoM1 395 G4double prob = ClebschGordan(twoJ1, twoM1Current, twoJ2,
396 twoM2Current, twoJ); 396 twoM2Current, twoJ);
397 sum += prob; 397 sum += prob;
398 if (twoM2Current == twoM2 && twoM1Current 398 if (twoM2Current == twoM2 && twoM1Current == twoM1) cleb += prob;
399 } 399 }
400 400
401 // Normalize probs to 1 401 // Normalize probs to 1
402 if (sum > 0.) cleb /= sum; 402 if (sum > 0.) cleb /= sum;
403 403
404 return cleb; 404 return cleb;
405 } 405 }
406 406
407 G4double G4Clebsch::TriangleCoeff(G4int twoA, 407 G4double G4Clebsch::TriangleCoeff(G4int twoA, G4int twoB, G4int twoC)
408 { 408 {
409 // TC(ABC) = sqrt[ (A+B-C)! (A-B+C)! (-A+B+C 409 // TC(ABC) = sqrt[ (A+B-C)! (A-B+C)! (-A+B+C)! / (A+B+C+1)! ]
410 // return 0 if the triad does not satisfy th 410 // return 0 if the triad does not satisfy the triangle inequalities
411 G4Pow* g4pow = G4Pow::GetInstance(); 411 G4Pow* g4pow = G4Pow::GetInstance();
412 412
413 double val = 0; 413 double val = 0;
414 G4int i = twoA+twoB-twoC; 414 G4int i = twoA+twoB-twoC;
415 // only have to check that i is even the fir 415 // only have to check that i is even the first time
416 if(i<0 || (i%2)) return 0; 416 if(i<0 || (i%2)) return 0;
417 else val += g4pow->logfactorial(i/2); 417 else val += g4pow->logfactorial(i/2);
418 418
419 i = twoA-twoB+twoC; 419 i = twoA-twoB+twoC;
420 if(i<0) return 0; 420 if(i<0) return 0;
421 else val += g4pow->logfactorial(i/2); 421 else val += g4pow->logfactorial(i/2);
422 422
423 i = -twoA+twoB+twoC; 423 i = -twoA+twoB+twoC;
424 if(i<0) return 0; 424 if(i<0) return 0;
425 else val += g4pow->logfactorial(i/2); 425 else val += g4pow->logfactorial(i/2);
426 426
427 i = twoA+twoB+twoC+2; 427 i = twoA+twoB+twoC+2;
428 if(i<0) return 0; 428 if(i<0) return 0;
429 return G4Exp(0.5*(val - g4pow->logfactorial( 429 return G4Exp(0.5*(val - g4pow->logfactorial(i/2)));
430 } 430 }
431 431
432 G4double G4Clebsch::Wigner6J(G4int twoJ1, G4in 432 G4double G4Clebsch::Wigner6J(G4int twoJ1, G4int twoJ2, G4int twoJ3,
433 G4int twoJ4, G4int twoJ5, G4i 433 G4int twoJ4, G4int twoJ5, G4int twoJ6)
434 { 434 {
435 if(twoJ1 < 0 || twoJ2 < 0 || twoJ3 < 0 || 435 if(twoJ1 < 0 || twoJ2 < 0 || twoJ3 < 0 ||
436 twoJ4 < 0 || twoJ5 < 0 || twoJ6 < 0) retu 436 twoJ4 < 0 || twoJ5 < 0 || twoJ6 < 0) return 0;
437 437
438 // There is a fast calculation (no sums or e 438 // There is a fast calculation (no sums or exps) when twoJ6 = 0,
439 // so permute to use it when possible 439 // so permute to use it when possible
440 if(twoJ6 == 0) { 440 if(twoJ6 == 0) {
441 if(twoJ1 != twoJ5) return 0; 441 if(twoJ1 != twoJ5) return 0;
442 if(twoJ2 != twoJ4) return 0; 442 if(twoJ2 != twoJ4) return 0;
443 if(twoJ1+twoJ2 < twoJ3) return 0; 443 if(twoJ1+twoJ2 < twoJ3) return 0;
444 if((twoJ1 > twoJ2) && (twoJ3 < (twoJ1-twoJ 444 if((twoJ1 > twoJ2) && (twoJ3 < (twoJ1-twoJ2))) return 0;
445 if((twoJ2 > twoJ1) && (twoJ3 < (twoJ2-twoJ 445 if((twoJ2 > twoJ1) && (twoJ3 < (twoJ2-twoJ1))) return 0;
446 if((twoJ1+twoJ2+twoJ3) % 2) return 0; 446 if((twoJ1+twoJ2+twoJ3) % 2) return 0;
447 return (((twoJ1+twoJ2+twoJ3)/2) % 2 ? -1. 447 return (((twoJ1+twoJ2+twoJ3)/2) % 2 ? -1. : 1.) /sqrt((twoJ1+1)*(twoJ2+1));
448 } 448 }
449 if(twoJ1 == 0) return Wigner6J(twoJ6, twoJ2, 449 if(twoJ1 == 0) return Wigner6J(twoJ6, twoJ2, twoJ4, twoJ3, twoJ5, 0);
450 if(twoJ2 == 0) return Wigner6J(twoJ1, twoJ6, 450 if(twoJ2 == 0) return Wigner6J(twoJ1, twoJ6, twoJ5, twoJ4, twoJ3, 0);
451 if(twoJ3 == 0) return Wigner6J(twoJ4, twoJ2, 451 if(twoJ3 == 0) return Wigner6J(twoJ4, twoJ2, twoJ6, twoJ1, twoJ5, 0);
452 if(twoJ4 == 0) return Wigner6J(twoJ3, twoJ2, 452 if(twoJ4 == 0) return Wigner6J(twoJ3, twoJ2, twoJ1, twoJ6, twoJ5, 0);
453 if(twoJ5 == 0) return Wigner6J(twoJ1, twoJ3, 453 if(twoJ5 == 0) return Wigner6J(twoJ1, twoJ3, twoJ2, twoJ4, twoJ6, 0);
454 454
455 // Check triangle inequalities and calculate 455 // Check triangle inequalities and calculate triangle coefficients.
456 // Also check evenness of sums 456 // Also check evenness of sums
457 G4Pow* g4pow = G4Pow::GetInstance(); 457 G4Pow* g4pow = G4Pow::GetInstance();
458 double triangles = 0; 458 double triangles = 0;
459 G4int i; 459 G4int i;
460 i = twoJ1+twoJ2-twoJ3; if(i<0 || i%2) ret 460 i = twoJ1+twoJ2-twoJ3; if(i<0 || i%2) return 0; else triangles += g4pow->logfactorial(i/2);
461 i = twoJ1-twoJ2+twoJ3; if(i<0 || i%2) ret 461 i = twoJ1-twoJ2+twoJ3; if(i<0 || i%2) return 0; else triangles += g4pow->logfactorial(i/2);
462 i = -twoJ1+twoJ2+twoJ3; if(i<0 || i%2) ret 462 i = -twoJ1+twoJ2+twoJ3; if(i<0 || i%2) return 0; else triangles += g4pow->logfactorial(i/2);
463 i = twoJ1+twoJ2+twoJ3+2; if(i<0 || i%2) ret 463 i = twoJ1+twoJ2+twoJ3+2; if(i<0 || i%2) return 0; else triangles -= g4pow->logfactorial(i/2);
464 i = twoJ1+twoJ5-twoJ6; if(i<0 || i%2) ret 464 i = twoJ1+twoJ5-twoJ6; if(i<0 || i%2) return 0; else triangles += g4pow->logfactorial(i/2);
465 i = twoJ1-twoJ5+twoJ6; if(i<0 || i%2) ret 465 i = twoJ1-twoJ5+twoJ6; if(i<0 || i%2) return 0; else triangles += g4pow->logfactorial(i/2);
466 i = -twoJ1+twoJ5+twoJ6; if(i<0 || i%2) ret 466 i = -twoJ1+twoJ5+twoJ6; if(i<0 || i%2) return 0; else triangles += g4pow->logfactorial(i/2);
467 i = twoJ1+twoJ5+twoJ6+2; if(i<0 || i%2) ret 467 i = twoJ1+twoJ5+twoJ6+2; if(i<0 || i%2) return 0; else triangles -= g4pow->logfactorial(i/2);
468 i = twoJ4+twoJ2-twoJ6; if(i<0 || i%2) ret 468 i = twoJ4+twoJ2-twoJ6; if(i<0 || i%2) return 0; else triangles += g4pow->logfactorial(i/2);
469 i = twoJ4-twoJ2+twoJ6; if(i<0 || i%2) ret 469 i = twoJ4-twoJ2+twoJ6; if(i<0 || i%2) return 0; else triangles += g4pow->logfactorial(i/2);
470 i = -twoJ4+twoJ2+twoJ6; if(i<0 || i%2) ret 470 i = -twoJ4+twoJ2+twoJ6; if(i<0 || i%2) return 0; else triangles += g4pow->logfactorial(i/2);
471 i = twoJ4+twoJ2+twoJ6+2; if(i<0 || i%2) ret 471 i = twoJ4+twoJ2+twoJ6+2; if(i<0 || i%2) return 0; else triangles -= g4pow->logfactorial(i/2);
472 i = twoJ4+twoJ5-twoJ3; if(i<0 || i%2) ret 472 i = twoJ4+twoJ5-twoJ3; if(i<0 || i%2) return 0; else triangles += g4pow->logfactorial(i/2);
473 i = twoJ4-twoJ5+twoJ3; if(i<0 || i%2) ret 473 i = twoJ4-twoJ5+twoJ3; if(i<0 || i%2) return 0; else triangles += g4pow->logfactorial(i/2);
474 i = -twoJ4+twoJ5+twoJ3; if(i<0 || i%2) ret 474 i = -twoJ4+twoJ5+twoJ3; if(i<0 || i%2) return 0; else triangles += g4pow->logfactorial(i/2);
475 i = twoJ4+twoJ5+twoJ3+2; if(i<0 || i%2) ret 475 i = twoJ4+twoJ5+twoJ3+2; if(i<0 || i%2) return 0; else triangles -= g4pow->logfactorial(i/2);
476 triangles = G4Exp(0.5*triangles); 476 triangles = G4Exp(0.5*triangles);
477 477
478 // Prepare to sum over k. If we have made it 478 // Prepare to sum over k. If we have made it this far, all of the following
479 // sums must be non-negative and divisible b 479 // sums must be non-negative and divisible by two
480 480
481 // k must be >= all of the following sums: 481 // k must be >= all of the following sums:
482 G4int sum1 = (twoJ1 + twoJ2 + twoJ3)/2; 482 G4int sum1 = (twoJ1 + twoJ2 + twoJ3)/2;
483 G4int kMin = sum1; 483 G4int kMin = sum1;
484 G4int sum2 = (twoJ1 + twoJ5 + twoJ6)/2; 484 G4int sum2 = (twoJ1 + twoJ5 + twoJ6)/2;
485 if(sum2 > kMin) kMin = sum2; 485 if(sum2 > kMin) kMin = sum2;
486 G4int sum3 = (twoJ4 + twoJ2 + twoJ6)/2; 486 G4int sum3 = (twoJ4 + twoJ2 + twoJ6)/2;
487 if(sum3 > kMin) kMin = sum3; 487 if(sum3 > kMin) kMin = sum3;
488 G4int sum4 = (twoJ4 + twoJ5 + twoJ3)/2; 488 G4int sum4 = (twoJ4 + twoJ5 + twoJ3)/2;
489 if(sum4 > kMin) kMin = sum4; 489 if(sum4 > kMin) kMin = sum4;
490 490
491 // and k must be <= all of the following sum 491 // and k must be <= all of the following sums:
492 G4int sum5 = (twoJ1 + twoJ2 + twoJ4 + twoJ5) 492 G4int sum5 = (twoJ1 + twoJ2 + twoJ4 + twoJ5)/2;
493 G4int kMax = sum5; 493 G4int kMax = sum5;
494 G4int sum6 = (twoJ2 + twoJ3 + twoJ5 + twoJ6) 494 G4int sum6 = (twoJ2 + twoJ3 + twoJ5 + twoJ6)/2;
495 if(sum6 < kMax) kMax = sum6; 495 if(sum6 < kMax) kMax = sum6;
496 G4int sum7 = (twoJ1 + twoJ3 + twoJ4 + twoJ6) 496 G4int sum7 = (twoJ1 + twoJ3 + twoJ4 + twoJ6)/2;
497 if(sum7 < kMax) kMax = sum7; 497 if(sum7 < kMax) kMax = sum7;
498 498
499 // sanity / boundary checks 499 // sanity / boundary checks
500 if(kMin < 0) { 500 if(kMin < 0) {
501 G4Exception("G4Clebsch::Wigner6J()", "Cleb 501 G4Exception("G4Clebsch::Wigner6J()", "Clebsch040",
502 JustWarning, "kMin < 0"); 502 JustWarning, "kMin < 0");
503 return 0; 503 return 0;
504 } 504 }
505 if(kMax < kMin) { 505 if(kMax < kMin) {
506 G4Exception("G4Clebsch::Wigner6J()", "Cleb 506 G4Exception("G4Clebsch::Wigner6J()", "Clebsch041",
507 JustWarning, "kMax < kMin"); 507 JustWarning, "kMax < kMin");
508 return 0; 508 return 0;
509 } 509 }
510 if(kMax >= G4POWLOGFACTMAX) { 510 if(kMax >= G4POWLOGFACTMAX) {
511 G4Exception("G4Clebsch::Wigner6J()", "Cleb 511 G4Exception("G4Clebsch::Wigner6J()", "Clebsch041",
512 JustWarning, "kMax too big for G4Pow"); 512 JustWarning, "kMax too big for G4Pow");
513 return 0; 513 return 0;
514 } 514 }
515 515
516 // Now do the sum over k 516 // Now do the sum over k
517 G4double kSum = 0.; 517 G4double kSum = 0.;
518 G4double sign = (kMin % 2) ? -1 : 1; 518 G4double sign = (kMin % 2) ? -1 : 1;
519 for(G4int k = kMin; k <= kMax; k++) { 519 for(G4int k = kMin; k <= kMax; k++) {
520 kSum += sign * G4Exp(g4pow->logfactorial(k 520 kSum += sign * G4Exp(g4pow->logfactorial(k+1) -
521 g4pow->logfactoria 521 g4pow->logfactorial(k-sum1) -
522 g4pow->logfactoria 522 g4pow->logfactorial(k-sum2) -
523 g4pow->logfactoria 523 g4pow->logfactorial(k-sum3) -
524 g4pow->logfactoria 524 g4pow->logfactorial(k-sum4) -
525 g4pow->logfactoria 525 g4pow->logfactorial(sum5-k) -
526 g4pow->logfactoria 526 g4pow->logfactorial(sum6-k) -
527 g4pow->logfactoria 527 g4pow->logfactorial(sum7-k));
528 sign *= -1; 528 sign *= -1;
529 } 529 }
530 return triangles*kSum; 530 return triangles*kSum;
531 } 531 }
532 532
533 G4double G4Clebsch::Wigner9J(G4int twoJ1, G4in 533 G4double G4Clebsch::Wigner9J(G4int twoJ1, G4int twoJ2, G4int twoJ3,
534 G4int twoJ4, G4int twoJ5, G4i 534 G4int twoJ4, G4int twoJ5, G4int twoJ6,
535 G4int twoJ7, G4int twoJ8, G4i 535 G4int twoJ7, G4int twoJ8, G4int twoJ9)
536 { 536 {
537 if(twoJ1 < 0 || twoJ2 < 0 || twoJ3 < 0 || 537 if(twoJ1 < 0 || twoJ2 < 0 || twoJ3 < 0 ||
538 twoJ4 < 0 || twoJ5 < 0 || twoJ6 < 0 || 538 twoJ4 < 0 || twoJ5 < 0 || twoJ6 < 0 ||
539 twoJ7 < 0 || twoJ8 < 0 || twoJ9 < 0) retu 539 twoJ7 < 0 || twoJ8 < 0 || twoJ9 < 0) return 0;
540 540
541 if(twoJ9 == 0) { 541 if(twoJ9 == 0) {
542 if(twoJ3 != twoJ6) return 0; 542 if(twoJ3 != twoJ6) return 0;
543 if(twoJ7 != twoJ8) return 0; 543 if(twoJ7 != twoJ8) return 0;
544 G4double sixJ = Wigner6J(twoJ1, twoJ2, two 544 G4double sixJ = Wigner6J(twoJ1, twoJ2, twoJ3, twoJ5, twoJ4, twoJ7);
545 if(sixJ == 0) return 0; 545 if(sixJ == 0) return 0;
546 if((twoJ2+twoJ3+twoJ4+twoJ7)/2 % 2) sixJ = 546 if((twoJ2+twoJ3+twoJ4+twoJ7)/2 % 2) sixJ = -sixJ;
547 return sixJ/sqrt((twoJ3+1)*(twoJ7+1)); 547 return sixJ/sqrt((twoJ3+1)*(twoJ7+1));
548 } 548 }
549 if(twoJ1 == 0) return Wigner9J(twoJ9, twoJ6, 549 if(twoJ1 == 0) return Wigner9J(twoJ9, twoJ6, twoJ3, twoJ8, twoJ5, twoJ2, twoJ7, twoJ4, twoJ1);
550 if(twoJ2 == 0) return Wigner9J(twoJ7, twoJ9, 550 if(twoJ2 == 0) return Wigner9J(twoJ7, twoJ9, twoJ8, twoJ4, twoJ6, twoJ5, twoJ1, twoJ3, twoJ2);
551 if(twoJ4 == 0) return Wigner9J(twoJ3, twoJ2, 551 if(twoJ4 == 0) return Wigner9J(twoJ3, twoJ2, twoJ1, twoJ9, twoJ8, twoJ7, twoJ6, twoJ5, twoJ4);
552 if(twoJ5 == 0) return Wigner9J(twoJ1, twoJ3, 552 if(twoJ5 == 0) return Wigner9J(twoJ1, twoJ3, twoJ2, twoJ7, twoJ9, twoJ8, twoJ4, twoJ6, twoJ5);
553 G4int twoS = twoJ1+twoJ2+twoJ3+twoJ4+twoJ5+t 553 G4int twoS = twoJ1+twoJ2+twoJ3+twoJ4+twoJ5+twoJ6+twoJ7+twoJ8+twoJ9;
554 if(twoS % 2) return 0; 554 if(twoS % 2) return 0;
555 G4double sign = (twoS/2 % 2) ? -1 : 1; 555 G4double sign = (twoS/2 % 2) ? -1 : 1;
556 if(twoJ3 == 0) return sign*Wigner9J(twoJ7, t 556 if(twoJ3 == 0) return sign*Wigner9J(twoJ7, twoJ8, twoJ9, twoJ4, twoJ5, twoJ6, twoJ1, twoJ2, twoJ3);
557 if(twoJ6 == 0) return sign*Wigner9J(twoJ1, t 557 if(twoJ6 == 0) return sign*Wigner9J(twoJ1, twoJ2, twoJ3, twoJ7, twoJ8, twoJ9, twoJ4, twoJ5, twoJ6);
558 if(twoJ7 == 0) return sign*Wigner9J(twoJ3, t 558 if(twoJ7 == 0) return sign*Wigner9J(twoJ3, twoJ2, twoJ1, twoJ6, twoJ5, twoJ4, twoJ9, twoJ8, twoJ7);
559 if(twoJ8 == 0) return sign*Wigner9J(twoJ1, t 559 if(twoJ8 == 0) return sign*Wigner9J(twoJ1, twoJ3, twoJ2, twoJ4, twoJ6, twoJ5, twoJ7, twoJ9, twoJ8);
560 560
561 // No element is zero: check triads now for 561 // No element is zero: check triads now for speed
562 G4int i; 562 G4int i;
563 i = twoJ1+twoJ2-twoJ3; if(i<0 || i%2) retur 563 i = twoJ1+twoJ2-twoJ3; if(i<0 || i%2) return 0;
564 i = twoJ1-twoJ2+twoJ3; if(i<0 || i%2) retur 564 i = twoJ1-twoJ2+twoJ3; if(i<0 || i%2) return 0;
565 i = -twoJ1+twoJ2+twoJ3; if(i<0 || i%2) retur 565 i = -twoJ1+twoJ2+twoJ3; if(i<0 || i%2) return 0;
566 i = twoJ4+twoJ5-twoJ6; if(i<0 || i%2) retur 566 i = twoJ4+twoJ5-twoJ6; if(i<0 || i%2) return 0;
567 i = twoJ4-twoJ5+twoJ6; if(i<0 || i%2) retur 567 i = twoJ4-twoJ5+twoJ6; if(i<0 || i%2) return 0;
568 i = -twoJ4+twoJ5+twoJ6; if(i<0 || i%2) retur 568 i = -twoJ4+twoJ5+twoJ6; if(i<0 || i%2) return 0;
569 i = twoJ7+twoJ8-twoJ9; if(i<0 || i%2) retur 569 i = twoJ7+twoJ8-twoJ9; if(i<0 || i%2) return 0;
570 i = twoJ7-twoJ8+twoJ9; if(i<0 || i%2) retur 570 i = twoJ7-twoJ8+twoJ9; if(i<0 || i%2) return 0;
571 i = -twoJ7+twoJ8+twoJ9; if(i<0 || i%2) retur 571 i = -twoJ7+twoJ8+twoJ9; if(i<0 || i%2) return 0;
572 i = twoJ1+twoJ4-twoJ7; if(i<0 || i%2) retur 572 i = twoJ1+twoJ4-twoJ7; if(i<0 || i%2) return 0;
573 i = twoJ1-twoJ4+twoJ7; if(i<0 || i%2) retur 573 i = twoJ1-twoJ4+twoJ7; if(i<0 || i%2) return 0;
574 i = -twoJ1+twoJ4+twoJ7; if(i<0 || i%2) retur 574 i = -twoJ1+twoJ4+twoJ7; if(i<0 || i%2) return 0;
575 i = twoJ2+twoJ5-twoJ8; if(i<0 || i%2) retur 575 i = twoJ2+twoJ5-twoJ8; if(i<0 || i%2) return 0;
576 i = twoJ2-twoJ5+twoJ8; if(i<0 || i%2) retur 576 i = twoJ2-twoJ5+twoJ8; if(i<0 || i%2) return 0;
577 i = -twoJ2+twoJ5+twoJ8; if(i<0 || i%2) retur 577 i = -twoJ2+twoJ5+twoJ8; if(i<0 || i%2) return 0;
578 i = twoJ3+twoJ6-twoJ9; if(i<0 || i%2) retur 578 i = twoJ3+twoJ6-twoJ9; if(i<0 || i%2) return 0;
579 i = twoJ3-twoJ6+twoJ9; if(i<0 || i%2) retur 579 i = twoJ3-twoJ6+twoJ9; if(i<0 || i%2) return 0;
580 i = -twoJ3+twoJ6+twoJ9; if(i<0 || i%2) retur 580 i = -twoJ3+twoJ6+twoJ9; if(i<0 || i%2) return 0;
581 581
582 // Okay, have to do the full sum over 6J's 582 // Okay, have to do the full sum over 6J's
583 // Find limits for K sum 583 // Find limits for K sum
584 G4int twoKMax = twoJ1+twoJ9; 584 G4int twoKMax = twoJ1+twoJ9;
585 if(twoJ4+twoJ8 < twoKMax) twoKMax = twoJ4+tw 585 if(twoJ4+twoJ8 < twoKMax) twoKMax = twoJ4+twoJ8;
586 if(twoJ2+twoJ6 < twoKMax) twoKMax = twoJ2+tw 586 if(twoJ2+twoJ6 < twoKMax) twoKMax = twoJ2+twoJ6;
587 G4int twoKMin = twoJ1-twoJ9; 587 G4int twoKMin = twoJ1-twoJ9;
588 if(twoJ9-twoJ1 > twoKMin) twoKMin = twoJ9-tw 588 if(twoJ9-twoJ1 > twoKMin) twoKMin = twoJ9-twoJ1;
589 if(twoJ4-twoJ8 > twoKMin) twoKMin = twoJ4-tw 589 if(twoJ4-twoJ8 > twoKMin) twoKMin = twoJ4-twoJ8;
590 if(twoJ8-twoJ4 > twoKMin) twoKMin = twoJ8-tw 590 if(twoJ8-twoJ4 > twoKMin) twoKMin = twoJ8-twoJ4;
591 if(twoJ2-twoJ6 > twoKMin) twoKMin = twoJ2-tw 591 if(twoJ2-twoJ6 > twoKMin) twoKMin = twoJ2-twoJ6;
592 if(twoJ6-twoJ2 > twoKMin) twoKMin = twoJ6-tw 592 if(twoJ6-twoJ2 > twoKMin) twoKMin = twoJ6-twoJ2;
593 if(twoKMin > twoKMax) return 0; 593 if(twoKMin > twoKMax) return 0;
594 594
595 G4double sum = 0; 595 G4double sum = 0;
596 for(G4int twoK = twoKMin; twoK <= twoKMax; t 596 for(G4int twoK = twoKMin; twoK <= twoKMax; twoK += 2) {
597 G4double value = Wigner6J(twoJ1, twoJ4, tw 597 G4double value = Wigner6J(twoJ1, twoJ4, twoJ7, twoJ8, twoJ9, twoK);
598 if(value == 0) continue; 598 if(value == 0) continue;
599 value *= Wigner6J(twoJ2, twoJ5, twoJ8, two 599 value *= Wigner6J(twoJ2, twoJ5, twoJ8, twoJ4, twoK, twoJ6);
600 if(value == 0) continue; 600 if(value == 0) continue;
601 value *= Wigner6J(twoJ3, twoJ6, twoJ9, two 601 value *= Wigner6J(twoJ3, twoJ6, twoJ9, twoK, twoJ1, twoJ2);
602 if(value == 0) continue; 602 if(value == 0) continue;
603 if(twoK % 2) value = -value; 603 if(twoK % 2) value = -value;
604 sum += value*G4double(twoK+1); 604 sum += value*G4double(twoK+1);
605 } 605 }
606 return sum; 606 return sum;
607 } 607 }
608 608
609 G4double G4Clebsch::WignerLittleD(G4int twoJ, 609 G4double G4Clebsch::WignerLittleD(G4int twoJ, G4int twoM, G4int twoN,
610 G4double cosTheta) 610 G4double cosTheta)
611 { 611 {
612 if(twoM < -twoJ || twoM > twoJ || twoN < -tw 612 if(twoM < -twoJ || twoM > twoJ || twoN < -twoJ || twoN > twoJ
613 || ((twoM % 2) != (twoJ % 2)) || ((twoN % 613 || ((twoM % 2) != (twoJ % 2)) || ((twoN % 2) != (twoJ % 2)))
614 { return 0; } 614 { return 0; }
615 615
616 if(cosTheta == 1.0) { return G4double(two 616 if(cosTheta == 1.0) { return G4double(twoM == twoN); }
617 617
618 G4int kMin = 0; 618 G4int kMin = 0;
619 if(twoM > twoN) kMin = (twoM-twoN)/2; 619 if(twoM > twoN) kMin = (twoM-twoN)/2;
620 G4int kMax = (twoJ + twoM)/2; 620 G4int kMax = (twoJ + twoM)/2;
621 if((twoJ-twoN)/2 < kMax) kMax = (twoJ-twoN)/ 621 if((twoJ-twoN)/2 < kMax) kMax = (twoJ-twoN)/2;
622 622
623 G4double lnCosHalfTheta = G4Log((cosTheta+1. 623 G4double lnCosHalfTheta = G4Log((cosTheta+1.)*0.5) * 0.5;
624 G4double lnSinHalfTheta = G4Log((1.-cosTheta 624 G4double lnSinHalfTheta = G4Log((1.-cosTheta)*0.5) * 0.5;
625 625
626 G4Pow* g4pow = G4Pow::GetInstance(); 626 G4Pow* g4pow = G4Pow::GetInstance();
627 G4double d = 0; 627 G4double d = 0;
628 for(G4int k = kMin; k <= kMax; k++) { 628 for(G4int k = kMin; k <= kMax; k++) {
629 G4double logSum = 0.5*(g4pow->logfactorial 629 G4double logSum = 0.5*(g4pow->logfactorial((twoJ+twoM)/2) +
630 g4pow->logfactorial 630 g4pow->logfactorial((twoJ-twoM)/2) +
631 g4pow->logfactorial 631 g4pow->logfactorial((twoJ+twoN)/2) +
632 g4pow->logfactorial 632 g4pow->logfactorial((twoJ-twoN)/2));
633 logSum += -g4pow->logfactorial((twoJ+twoM) 633 logSum += -g4pow->logfactorial((twoJ+twoM)/2 - k) -
634 g4pow->logfactorial((twoJ-twoN) 634 g4pow->logfactorial((twoJ-twoN)/2 - k) -
635 g4pow->logfactorial(k) - 635 g4pow->logfactorial(k) -
636 g4pow->logfactorial(k+(twoN-two 636 g4pow->logfactorial(k+(twoN-twoM)/2);
637 logSum += (twoJ+(twoM-twoN)/2 - 2*k)*lnCos 637 logSum += (twoJ+(twoM-twoN)/2 - 2*k)*lnCosHalfTheta +
638 (2*k + (twoN-twoM)/2)*lnSinHalfTheta; 638 (2*k + (twoN-twoM)/2)*lnSinHalfTheta;
639 G4double sign = (k % 2) ? -1 : 1; 639 G4double sign = (k % 2) ? -1 : 1;
640 d += sign * G4Exp(logSum); 640 d += sign * G4Exp(logSum);
641 } 641 }
642 return d; 642 return d;
643 } 643 }
644 644