Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Static Public Member Functions | List of all members
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< G4double > GenerateIso3 (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

◆ ClebschGordan()

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

Referenced by GenerateIso3(), NormalizedClebschGordan(), and Weight().

◆ ClebschGordanCoeff()

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}
G4POWLOGFACTMAX
const G4int G4POWLOGFACTMAX
Definition: G4Clebsch.cc:33
JustWarning
@ JustWarning
Definition: G4ExceptionSeverity.hh:65
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:59
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
G4int
int G4int
Definition: G4Types.hh:85
G4Clebsch::TriangleCoeff
static G4double TriangleCoeff(G4int twoA, G4int twoB, G4int twoC)
Definition: G4Clebsch.cc:407
G4Pow
Definition: G4Pow.hh:49
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition: G4Pow.cc:41
G4Pow::logfactorial
G4double logfactorial(G4int Z) const
Definition: G4Pow.hh:237
G4INCL::Math::sign
G4int sign(const T t)
Definition: G4INCLGlobals.hh:107

Referenced by ClebschGordan(), and Wigner3J().

◆ GenerateIso3()

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}
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:52
G4Clebsch::ClebschGordan
static G4double ClebschGordan(G4int twoJ1, G4int twoM1, G4int twoJ2, G4int twoM2, G4int twoJ)
Definition: G4Clebsch.cc:106

Referenced by G4VXResonance::IsospinCorrection().

◆ NormalizedClebschGordan()

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}

◆ RacahWCoeff()

static G4double G4Clebsch::RacahWCoeff ( G4int  twoJ1,
G4int  twoJ2,
G4int  twoJ,
G4int  twoJ3,
G4int  twoJ12,
G4int  twoJ23 
)
static

◆ TriangleCoeff()

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}

Referenced by ClebschGordanCoeff().

◆ Weight()

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}

Referenced by G4VXResonance::DetailedBalance(), and G4VXResonance::IsospinCorrection().

◆ Wigner3J() [1/3]

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}
G4Clebsch::Wigner3J
static G4double Wigner3J(G4double j1, G4double j2, G4double j3, G4double m1, G4double m2, G4double m3)
Definition: G4Clebsch.cc:345

Referenced by G4PolarizationTransition::F3Coefficient(), G4PolarizationTransition::FCoefficient(), G4PolarizationTransition::SampleGammaTransition(), and Wigner3J().

◆ Wigner3J() [2/3]

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}

◆ Wigner3J() [3/3]

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}

◆ Wigner6J()

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}
G4Clebsch::Wigner6J
static G4double Wigner6J(G4int twoJ1, G4int twoJ2, G4int twoJ3, G4int twoJ4, G4int twoJ5, G4int twoJ6)
Definition: G4Clebsch.cc:432

Referenced by G4PolarizationTransition::FCoefficient(), Wigner6J(), and Wigner9J().

◆ Wigner9J()

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}
G4Clebsch::Wigner9J
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

Referenced by G4PolarizationTransition::F3Coefficient(), and Wigner9J().

◆ WignerLittleD()

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}
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:227
The documentation for this class was generated from the following files: