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