Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4Clebsch Class Reference

#include <G4Clebsch.hh>

Static Public Member Functions

static G4double ClebschGordanCoeff (G4int twoJ1, G4int twoM1, G4int twoJ2, G4int twoM2, G4int twoJ)
 
static G4double ClebschGordan (G4int twoJ1, G4int twoM1, G4int twoJ2, G4int twoM2, G4int twoJ)
 
static std::vector< G4doubleGenerateIso3 (G4int twoJ1, G4int twoM1, G4int twoJ2, G4int twoM2, G4int twoJOut1, G4int twoJOut2)
 
static G4double Weight (G4int twoJ1, G4int twoM1, G4int twoJ2, G4int twoM2, G4int twoJOut1, G4int twoJOut2)
 
static G4double Wigner3J (G4double j1, G4double j2, G4double j3, G4double m1, G4double m2, G4double m3)
 
static G4double Wigner3J (G4int twoJ1, G4int twoM1, G4int twoJ2, G4int twoM2, G4int twoJ3)
 
static G4double Wigner3J (G4int twoJ1, G4int twoM1, G4int twoJ2, G4int twoM2, G4int twoJ3, G4int twoM3)
 
static G4double NormalizedClebschGordan (G4int twoJ, G4int twom, G4int twoJ1, G4int twoJ2, G4int twom1, G4int twom2)
 
static G4double TriangleCoeff (G4int twoA, G4int twoB, G4int twoC)
 
static G4double Wigner6J (G4int twoJ1, G4int twoJ2, G4int twoJ3, G4int twoJ4, G4int twoJ5, G4int twoJ6)
 
static G4double RacahWCoeff (G4int twoJ1, G4int twoJ2, G4int twoJ, G4int twoJ3, G4int twoJ12, G4int twoJ23)
 
static G4double Wigner9J (G4int twoJ1, G4int twoJ2, G4int twoJ3, G4int twoJ4, G4int twoJ5, G4int twoJ6, G4int twoJ7, G4int twoJ8, G4int twoJ9)
 
static G4double WignerLittleD (G4int twoJ, G4int twoM, G4int twoN, G4double cosTheta)
 

Detailed Description

Definition at line 62 of file G4Clebsch.hh.

Member Function Documentation

G4double G4Clebsch::ClebschGordan ( G4int  twoJ1,
G4int  twoM1,
G4int  twoJ2,
G4int  twoM2,
G4int  twoJ 
)
static

Definition at line 106 of file G4Clebsch.cc.

109 {
110  // ClebschGordanCoeff() will do all input checking
111  G4double clebsch = ClebschGordanCoeff(twoJ1, twoM1, twoJ2, twoM2, twoJ);
112  return clebsch*clebsch;
113 }
static G4double ClebschGordanCoeff(G4int twoJ1, G4int twoM1, G4int twoJ2, G4int twoM2, G4int twoJ)
Definition: G4Clebsch.cc:37
double G4double
Definition: G4Types.hh:76
G4double G4Clebsch::ClebschGordanCoeff ( G4int  twoJ1,
G4int  twoM1,
G4int  twoJ2,
G4int  twoM2,
G4int  twoJ 
)
static

Definition at line 37 of file G4Clebsch.cc.

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 }
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
Definition: G4Pow.hh:56
const G4int G4POWLOGFACTMAX
Definition: G4Clebsch.cc:33
int G4int
Definition: G4Types.hh:78
static G4double TriangleCoeff(G4int twoA, G4int twoB, G4int twoC)
Definition: G4Clebsch.cc:407
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double logfactorial(G4int Z) const
Definition: G4Pow.hh:269
double G4double
Definition: G4Types.hh:76
G4int sign(const T t)

Here is the call graph for this function:

std::vector< G4double > G4Clebsch::GenerateIso3 ( G4int  twoJ1,
G4int  twoM1,
G4int  twoJ2,
G4int  twoM2,
G4int  twoJOut1,
G4int  twoJOut2 
)
static

Definition at line 116 of file G4Clebsch.cc.

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 }
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:97
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
T max(const T t1, const T t2)
brief Return the largest of the two arguments
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
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4Clebsch::NormalizedClebschGordan ( G4int  twoJ,
G4int  twom,
G4int  twoJ1,
G4int  twoJ2,
G4int  twom1,
G4int  twom2 
)
static

Definition at line 380 of file G4Clebsch.cc.

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 }
int G4int
Definition: G4Types.hh:78
static G4double ClebschGordan(G4int twoJ1, G4int twoM1, G4int twoJ2, G4int twoM2, G4int twoJ)
Definition: G4Clebsch.cc:106
double G4double
Definition: G4Types.hh:76
static G4double G4Clebsch::RacahWCoeff ( G4int  twoJ1,
G4int  twoJ2,
G4int  twoJ,
G4int  twoJ3,
G4int  twoJ12,
G4int  twoJ23 
)
static
G4double G4Clebsch::TriangleCoeff ( G4int  twoA,
G4int  twoB,
G4int  twoC 
)
static

Definition at line 407 of file G4Clebsch.cc.

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 }
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
Definition: G4Pow.hh:56
int G4int
Definition: G4Types.hh:78
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double logfactorial(G4int Z) const
Definition: G4Pow.hh:269

Here is the call graph for this function:

G4double G4Clebsch::Weight ( G4int  twoJ1,
G4int  twoM1,
G4int  twoJ2,
G4int  twoM2,
G4int  twoJOut1,
G4int  twoJOut2 
)
static

Definition at line 320 of file G4Clebsch.cc.

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 }
int G4int
Definition: G4Types.hh:78
const XML_Char int const XML_Char * value
Definition: expat.h:331
T max(const T t1, const T t2)
brief Return the largest of the two arguments
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
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4Clebsch::Wigner3J ( G4double  j1,
G4double  j2,
G4double  j3,
G4double  m1,
G4double  m2,
G4double  m3 
)
static

Definition at line 345 of file G4Clebsch.cc.

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 }
static constexpr double m3
Definition: G4SIunits.hh:131
int G4int
Definition: G4Types.hh:78
static G4double Wigner3J(G4double j1, G4double j2, G4double j3, G4double m1, G4double m2, G4double m3)
Definition: G4Clebsch.cc:345
static constexpr double m2
Definition: G4SIunits.hh:130

Here is the caller graph for this function:

G4double G4Clebsch::Wigner3J ( G4int  twoJ1,
G4int  twoM1,
G4int  twoJ2,
G4int  twoM2,
G4int  twoJ3 
)
static

Definition at line 359 of file G4Clebsch.cc.

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 }
static G4double ClebschGordanCoeff(G4int twoJ1, G4int twoM1, G4int twoJ2, G4int twoM2, G4int twoJ)
Definition: G4Clebsch.cc:37
double G4double
Definition: G4Types.hh:76
G4double G4Clebsch::Wigner3J ( G4int  twoJ1,
G4int  twoM1,
G4int  twoJ2,
G4int  twoM2,
G4int  twoJ3,
G4int  twoM3 
)
static

Definition at line 369 of file G4Clebsch.cc.

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 }
static G4double ClebschGordanCoeff(G4int twoJ1, G4int twoM1, G4int twoJ2, G4int twoM2, G4int twoJ)
Definition: G4Clebsch.cc:37
double G4double
Definition: G4Types.hh:76
G4double G4Clebsch::Wigner6J ( G4int  twoJ1,
G4int  twoJ2,
G4int  twoJ3,
G4int  twoJ4,
G4int  twoJ5,
G4int  twoJ6 
)
static

Definition at line 432 of file G4Clebsch.cc.

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 }
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
Definition: G4Pow.hh:56
const G4int G4POWLOGFACTMAX
Definition: G4Clebsch.cc:33
int G4int
Definition: G4Types.hh:78
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double logfactorial(G4int Z) const
Definition: G4Pow.hh:269
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)

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4Clebsch::Wigner9J ( G4int  twoJ1,
G4int  twoJ2,
G4int  twoJ3,
G4int  twoJ4,
G4int  twoJ5,
G4int  twoJ6,
G4int  twoJ7,
G4int  twoJ8,
G4int  twoJ9 
)
static

Definition at line 533 of file G4Clebsch.cc.

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 }
int G4int
Definition: G4Types.hh:78
const XML_Char int const XML_Char * value
Definition: expat.h:331
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
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)

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4Clebsch::WignerLittleD ( G4int  twoJ,
G4int  twoM,
G4int  twoN,
G4double  cosTheta 
)
static

Definition at line 609 of file G4Clebsch.cc.

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
Definition: G4Pow.hh:56
int G4int
Definition: G4Types.hh:78
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double logfactorial(G4int Z) const
Definition: G4Pow.hh:269
double G4double
Definition: G4Types.hh:76
G4int sign(const T t)

Here is the call graph for this function:


The documentation for this class was generated from the following files: