| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 t
38 G4int t
39 G4int t
40 {
41 if(twoJ1 < 0 || twoJ2 < 0 || twoJ < 0 ||
42 ((twoJ1-twoM1) % 2) || ((twoJ2-twoM2) % 2
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, two
51 if(triangle == 0) { return 0; }
52
53 G4Pow* g4pow = G4Pow::GetInstance();
54 G4double factor = g4pow->logfactorial((twoJ1
55 g4pow->logfactorial((twoJ1
56 factor += g4pow->logfactorial((twoJ2 + twoM2
57 g4pow->logfactorial((twoJ2 - twoM2
58 factor += g4pow->logfactorial((twoJ + twoM)/
59 g4pow->logfactorial((twoJ - twoM)/
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
77 JustWarning, "kMin < 0");
78 return 0;
79 }
80 if(kMax < kMin) {
81 G4Exception("G4Clebsch::ClebschGordanCoeff
82 JustWarning, "kMax < kMin");
83 return 0;
84 }
85 if(kMax >= G4POWLOGFACTMAX) {
86 G4Exception("G4Clebsch::ClebschGordanCoeff
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->logfa
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,
107 G4int twoJ2, G4int twoM2,
108 G4int twoJ)
109 {
110 // ClebschGordanCoeff() will do all input ch
111 G4double clebsch = ClebschGordanCoeff(twoJ1,
112 return clebsch*clebsch;
113 }
114
115 std::vector<G4double>
116 G4Clebsch::GenerateIso3(G4int twoJ1, G4int two
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()", "
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 -
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*t
156 if(twoJTmp < twoJMinOut) twoJMinOut = tw
157 }
158 }
159 twoJMinOut = std::max(twoJMinOut, std::abs(t
160 G4int twoJMaxOut = twoJOut1 + twoJOut2;
161
162 // Possible in and out common states
163 G4int twoJMin = std::max(twoJMinIn, twoJMi
164 G4int twoJMax = std::min(twoJMaxIn, twoJMa
165 if (twoJMin > twoJMax) {
166 G4Exception("G4Clebsch::GenerateIso3()", "
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 !
177 G4Exception("G4Clebsch::GenerateIso3()", "
178 JustWarning, "twoJ1 or twoJ2 = 0, but twoJ
179 return temp;
180 }
181
182 // MGP ---- Shall it be a warning or an exce
183 if (nJ == 0) {
184 G4Exception("G4Clebsch::GenerateIso3()", "
185 JustWarning, "nJ is zero, no overlap betwe
186 return temp;
187 }
188
189 // Loop over all possible combinations of tw
190 // to get the probability of each of the in-
191
192 std::vector<G4double> clebsch;
193 G4double sum = 0.0;
194 for(G4int twoJ=twoJMin; twoJ<=twoJMax; twoJ+
195 sum += ClebschGordan(twoJ1, twoM1, twoJ2,
196 clebsch.push_back(sum);
197 }
198
199 // Consistency check
200 if (static_cast<G4int>(clebsch.size()) != nJ
201 G4Exception("G4Clebsch::GenerateIso3()", "
202 JustWarning, "nJ inconsistency");
203 return temp;
204 }
205
206 // Consistency check
207 if (sum <= 0.) {
208 G4Exception("G4Clebsch::GenerateIso3()", "
209 JustWarning, "Sum of Clebsch-Gordan probab
210 return temp;
211 }
212
213 // Generate a random twoJTot according to th
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> combinat
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<=
247 {
248 m1pos = -1;
249 for (m1pr = static_cast<G4int>(mMin[0]+.00
250 {
251 m1pos++;
252 if (m1pos >= size) {
253 G4Exception("G4Clebsch::GenerateIso3()
254 JustWarning, "m1pos > size");
255 return temp;
256 }
257 m1Out.push_back(m1pr);
258 m2pos = -1;
259 for (m2pr = static_cast<G4int>(mMin[1]+.
260 {
261 m2pos++;
262 if (m2pos >= size)
263 {
264 G4Exception("G4Clebsch::GenerateIso3
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(twoJOut
274 G4double c34 = ClebschGordan(0,0,0,0
275 G4double ctot = ClebschGordan(j12, m
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()", "
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()", "Cl
316 JustWarning, "Should never get here");
317 return temp;
318 }
319
320 G4double G4Clebsch::Weight(G4int twoJ1, G4int
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 -
329 G4int twoJMaxIn = twoJ1 + twoJ2;
330
331 G4int twoJMinOut = std::max(std::abs(twoJOut
332 G4int twoJMaxOut = twoJOut1 + twoJOut2;
333
334 G4int twoJMin = std::max(twoJMinIn,twoJMinOu
335 G4int twoJMax = std::min(twoJMaxIn,twoJMaxOu
336
337 for (G4int twoJ=twoJMin; twoJ<=twoJMax; twoJ
338 // ClebschGordan() will do all input check
339 value += ClebschGordan(twoJ1, twoM1, twoJ2
340 }
341
342 return value;
343 }
344
345 G4double G4Clebsch::Wigner3J(G4double j1, G4do
346 G4double m1, G4double m2, G4double
347 {
348 // G4Exception("G4Clebsch::Wigner3J()", "Cl
349 // "G4Clebsch::Wigner3J with double argumen
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,
357 }
358
359 G4double G4Clebsch::Wigner3J(G4int twoJ1, G4in
360 G4int twoJ2, G4in
361 G4int twoJ3)
362 {
363 G4double clebsch = ClebschGordanCoeff(twoJ1,
364 if(clebsch == 0) return clebsch;
365 if( (twoJ1-twoJ2+twoM1+twoM2)/2 % 2) clebsch
366 return clebsch / sqrt(twoJ3+1);
367 }
368
369 G4double G4Clebsch::Wigner3J(G4int twoJ1, G4in
370 G4int twoJ2, G4in
371 G4int twoJ3, G4in
372 {
373 if(twoM1 + twoM2 != -twoM3) return 0;
374 G4double clebsch = ClebschGordanCoeff(twoJ1,
375 if(clebsch == 0) return clebsch;
376 if( (twoJ1-twoJ2-twoM3)/2 % 2) clebsch = -cl
377 return clebsch / sqrt(twoJ3+1);
378 }
379
380 G4double G4Clebsch::NormalizedClebschGordan(G4
381 G4int twoJ1, G4int twoJ2,
382 G4int twoM1, G4int twoM2)
383 {
384 // Calculate the normalized Clebsch-Gordan c
385 // of isospin decomposition of (J,m) into J1
386
387 G4double cleb = 0.;
388 if(twoJ1 == 0 || twoJ2 == 0) return cleb;
389
390 // Loop over all J1,J2,Jtot,m1,m2 combinatio
391 G4double sum = 0.0;
392 for(G4int twoM1Current=-twoJ1; twoM1Current<
393 G4int twoM2Current = twoM - twoM1Current;
394 // ClebschGordan() will do all further inp
395 G4double prob = ClebschGordan(twoJ1, twoM1
396 twoM2Current, twoJ);
397 sum += prob;
398 if (twoM2Current == twoM2 && twoM1Current
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,
408 {
409 // TC(ABC) = sqrt[ (A+B-C)! (A-B+C)! (-A+B+C
410 // return 0 if the triad does not satisfy th
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 fir
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(
430 }
431
432 G4double G4Clebsch::Wigner6J(G4int twoJ1, G4in
433 G4int twoJ4, G4int twoJ5, G4i
434 {
435 if(twoJ1 < 0 || twoJ2 < 0 || twoJ3 < 0 ||
436 twoJ4 < 0 || twoJ5 < 0 || twoJ6 < 0) retu
437
438 // There is a fast calculation (no sums or e
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-twoJ
445 if((twoJ2 > twoJ1) && (twoJ3 < (twoJ2-twoJ
446 if((twoJ1+twoJ2+twoJ3) % 2) return 0;
447 return (((twoJ1+twoJ2+twoJ3)/2) % 2 ? -1.
448 }
449 if(twoJ1 == 0) return Wigner6J(twoJ6, twoJ2,
450 if(twoJ2 == 0) return Wigner6J(twoJ1, twoJ6,
451 if(twoJ3 == 0) return Wigner6J(twoJ4, twoJ2,
452 if(twoJ4 == 0) return Wigner6J(twoJ3, twoJ2,
453 if(twoJ5 == 0) return Wigner6J(twoJ1, twoJ3,
454
455 // Check triangle inequalities and calculate
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) ret
461 i = twoJ1-twoJ2+twoJ3; if(i<0 || i%2) ret
462 i = -twoJ1+twoJ2+twoJ3; if(i<0 || i%2) ret
463 i = twoJ1+twoJ2+twoJ3+2; if(i<0 || i%2) ret
464 i = twoJ1+twoJ5-twoJ6; if(i<0 || i%2) ret
465 i = twoJ1-twoJ5+twoJ6; if(i<0 || i%2) ret
466 i = -twoJ1+twoJ5+twoJ6; if(i<0 || i%2) ret
467 i = twoJ1+twoJ5+twoJ6+2; if(i<0 || i%2) ret
468 i = twoJ4+twoJ2-twoJ6; if(i<0 || i%2) ret
469 i = twoJ4-twoJ2+twoJ6; if(i<0 || i%2) ret
470 i = -twoJ4+twoJ2+twoJ6; if(i<0 || i%2) ret
471 i = twoJ4+twoJ2+twoJ6+2; if(i<0 || i%2) ret
472 i = twoJ4+twoJ5-twoJ3; if(i<0 || i%2) ret
473 i = twoJ4-twoJ5+twoJ3; if(i<0 || i%2) ret
474 i = -twoJ4+twoJ5+twoJ3; if(i<0 || i%2) ret
475 i = twoJ4+twoJ5+twoJ3+2; if(i<0 || i%2) ret
476 triangles = G4Exp(0.5*triangles);
477
478 // Prepare to sum over k. If we have made it
479 // sums must be non-negative and divisible b
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 sum
492 G4int sum5 = (twoJ1 + twoJ2 + twoJ4 + twoJ5)
493 G4int kMax = sum5;
494 G4int sum6 = (twoJ2 + twoJ3 + twoJ5 + twoJ6)
495 if(sum6 < kMax) kMax = sum6;
496 G4int sum7 = (twoJ1 + twoJ3 + twoJ4 + twoJ6)
497 if(sum7 < kMax) kMax = sum7;
498
499 // sanity / boundary checks
500 if(kMin < 0) {
501 G4Exception("G4Clebsch::Wigner6J()", "Cleb
502 JustWarning, "kMin < 0");
503 return 0;
504 }
505 if(kMax < kMin) {
506 G4Exception("G4Clebsch::Wigner6J()", "Cleb
507 JustWarning, "kMax < kMin");
508 return 0;
509 }
510 if(kMax >= G4POWLOGFACTMAX) {
511 G4Exception("G4Clebsch::Wigner6J()", "Cleb
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
521 g4pow->logfactoria
522 g4pow->logfactoria
523 g4pow->logfactoria
524 g4pow->logfactoria
525 g4pow->logfactoria
526 g4pow->logfactoria
527 g4pow->logfactoria
528 sign *= -1;
529 }
530 return triangles*kSum;
531 }
532
533 G4double G4Clebsch::Wigner9J(G4int twoJ1, G4in
534 G4int twoJ4, G4int twoJ5, G4i
535 G4int twoJ7, G4int twoJ8, G4i
536 {
537 if(twoJ1 < 0 || twoJ2 < 0 || twoJ3 < 0 ||
538 twoJ4 < 0 || twoJ5 < 0 || twoJ6 < 0 ||
539 twoJ7 < 0 || twoJ8 < 0 || twoJ9 < 0) retu
540
541 if(twoJ9 == 0) {
542 if(twoJ3 != twoJ6) return 0;
543 if(twoJ7 != twoJ8) return 0;
544 G4double sixJ = Wigner6J(twoJ1, twoJ2, two
545 if(sixJ == 0) return 0;
546 if((twoJ2+twoJ3+twoJ4+twoJ7)/2 % 2) sixJ =
547 return sixJ/sqrt((twoJ3+1)*(twoJ7+1));
548 }
549 if(twoJ1 == 0) return Wigner9J(twoJ9, twoJ6,
550 if(twoJ2 == 0) return Wigner9J(twoJ7, twoJ9,
551 if(twoJ4 == 0) return Wigner9J(twoJ3, twoJ2,
552 if(twoJ5 == 0) return Wigner9J(twoJ1, twoJ3,
553 G4int twoS = twoJ1+twoJ2+twoJ3+twoJ4+twoJ5+t
554 if(twoS % 2) return 0;
555 G4double sign = (twoS/2 % 2) ? -1 : 1;
556 if(twoJ3 == 0) return sign*Wigner9J(twoJ7, t
557 if(twoJ6 == 0) return sign*Wigner9J(twoJ1, t
558 if(twoJ7 == 0) return sign*Wigner9J(twoJ3, t
559 if(twoJ8 == 0) return sign*Wigner9J(twoJ1, t
560
561 // No element is zero: check triads now for
562 G4int i;
563 i = twoJ1+twoJ2-twoJ3; if(i<0 || i%2) retur
564 i = twoJ1-twoJ2+twoJ3; if(i<0 || i%2) retur
565 i = -twoJ1+twoJ2+twoJ3; if(i<0 || i%2) retur
566 i = twoJ4+twoJ5-twoJ6; if(i<0 || i%2) retur
567 i = twoJ4-twoJ5+twoJ6; if(i<0 || i%2) retur
568 i = -twoJ4+twoJ5+twoJ6; if(i<0 || i%2) retur
569 i = twoJ7+twoJ8-twoJ9; if(i<0 || i%2) retur
570 i = twoJ7-twoJ8+twoJ9; if(i<0 || i%2) retur
571 i = -twoJ7+twoJ8+twoJ9; if(i<0 || i%2) retur
572 i = twoJ1+twoJ4-twoJ7; if(i<0 || i%2) retur
573 i = twoJ1-twoJ4+twoJ7; if(i<0 || i%2) retur
574 i = -twoJ1+twoJ4+twoJ7; if(i<0 || i%2) retur
575 i = twoJ2+twoJ5-twoJ8; if(i<0 || i%2) retur
576 i = twoJ2-twoJ5+twoJ8; if(i<0 || i%2) retur
577 i = -twoJ2+twoJ5+twoJ8; if(i<0 || i%2) retur
578 i = twoJ3+twoJ6-twoJ9; if(i<0 || i%2) retur
579 i = twoJ3-twoJ6+twoJ9; if(i<0 || i%2) retur
580 i = -twoJ3+twoJ6+twoJ9; if(i<0 || i%2) retur
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+tw
586 if(twoJ2+twoJ6 < twoKMax) twoKMax = twoJ2+tw
587 G4int twoKMin = twoJ1-twoJ9;
588 if(twoJ9-twoJ1 > twoKMin) twoKMin = twoJ9-tw
589 if(twoJ4-twoJ8 > twoKMin) twoKMin = twoJ4-tw
590 if(twoJ8-twoJ4 > twoKMin) twoKMin = twoJ8-tw
591 if(twoJ2-twoJ6 > twoKMin) twoKMin = twoJ2-tw
592 if(twoJ6-twoJ2 > twoKMin) twoKMin = twoJ6-tw
593 if(twoKMin > twoKMax) return 0;
594
595 G4double sum = 0;
596 for(G4int twoK = twoKMin; twoK <= twoKMax; t
597 G4double value = Wigner6J(twoJ1, twoJ4, tw
598 if(value == 0) continue;
599 value *= Wigner6J(twoJ2, twoJ5, twoJ8, two
600 if(value == 0) continue;
601 value *= Wigner6J(twoJ3, twoJ6, twoJ9, two
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,
610 G4double cosTheta)
611 {
612 if(twoM < -twoJ || twoM > twoJ || twoN < -tw
613 || ((twoM % 2) != (twoJ % 2)) || ((twoN %
614 { return 0; }
615
616 if(cosTheta == 1.0) { return G4double(two
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)/
622
623 G4double lnCosHalfTheta = G4Log((cosTheta+1.
624 G4double lnSinHalfTheta = G4Log((1.-cosTheta
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
630 g4pow->logfactorial
631 g4pow->logfactorial
632 g4pow->logfactorial
633 logSum += -g4pow->logfactorial((twoJ+twoM)
634 g4pow->logfactorial((twoJ-twoN)
635 g4pow->logfactorial(k) -
636 g4pow->logfactorial(k+(twoN-two
637 logSum += (twoJ+(twoM-twoN)/2 - 2*k)*lnCos
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