Geant4  10.02.p01
G4Clebsch.cc
Go to the documentation of this file.
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 
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 
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>
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 }
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
static G4double Weight(G4int twoJ1, G4int twoM1, G4int twoJ2, G4int twoM2, G4int twoJOut1, G4int twoJOut2)
Definition: G4Clebsch.cc:320
Definition: G4Pow.hh:56
static const double m3
Definition: G4SIunits.hh:130
const G4int G4POWLOGFACTMAX
Definition: G4Clebsch.cc:33
static G4double ClebschGordanCoeff(G4int twoJ1, G4int twoM1, G4int twoJ2, G4int twoM2, G4int twoJ)
Definition: G4Clebsch.cc:37
int G4int
Definition: G4Types.hh:78
static G4double TriangleCoeff(G4int twoA, G4int twoB, G4int twoC)
Definition: G4Clebsch.cc:407
static G4double NormalizedClebschGordan(G4int twoJ, G4int twom, G4int twoJ1, G4int twoJ2, G4int twom1, G4int twom2)
Definition: G4Clebsch.cc:380
#define G4UniformRand()
Definition: Randomize.hh:97
static const double m2
Definition: G4SIunits.hh:129
static G4double Wigner9J(G4int twoJ1, G4int twoJ2, G4int twoJ3, G4int twoJ4, G4int twoJ5, G4int twoJ6, G4int twoJ7, G4int twoJ8, G4int twoJ9)
Definition: G4Clebsch.cc:533
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static std::vector< G4double > GenerateIso3(G4int twoJ1, G4int twoM1, G4int twoJ2, G4int twoM2, G4int twoJOut1, G4int twoJOut2)
Definition: G4Clebsch.cc:116
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static G4double Wigner3J(G4double j1, G4double j2, G4double j3, G4double m1, G4double m2, G4double m3)
Definition: G4Clebsch.cc:345
static const G4double factor
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double logfactorial(G4int Z) const
Definition: G4Pow.hh:269
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static G4double ClebschGordan(G4int twoJ1, G4int twoM1, G4int twoJ2, G4int twoM2, G4int twoJ)
Definition: G4Clebsch.cc:106
static G4double Wigner6J(G4int twoJ1, G4int twoJ2, G4int twoJ3, G4int twoJ4, G4int twoJ5, G4int twoJ6)
Definition: G4Clebsch.cc:432
double G4double
Definition: G4Types.hh:76
G4int sign(const T t)
A simple sign function that allows us to port fortran code to c++ more easily.
static G4double WignerLittleD(G4int twoJ, G4int twoM, G4int twoN, G4double cosTheta)
Definition: G4Clebsch.cc:609