41 if(twoJ1 < 0 || twoJ2 < 0 || twoJ < 0 ||
42 ((twoJ1-twoM1) % 2) || ((twoJ2-twoM2) % 2)) {
return 0; }
44 G4int twoM = twoM1 + twoM2;
45 if(twoM1 > twoJ1 || twoM1 < -twoJ1 ||
46 twoM2 > twoJ2 || twoM2 < -twoJ2 ||
47 twoM > twoJ || twoM < -twoJ) {
return 0; }
50 G4double triangle = TriangleCoeff(twoJ1, twoJ2, twoJ);
51 if(triangle == 0) {
return 0; }
63 G4int sum1 = (twoJ1 - twoM1)/2;
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;
76 G4Exception(
"G4Clebsch::ClebschGordanCoeff()",
"Clebsch001",
81 G4Exception(
"G4Clebsch::ClebschGordanCoeff()",
"Clebsch002",
86 G4Exception(
"G4Clebsch::ClebschGordanCoeff()",
"Clebsch003",
103 return triangle*sqrt(twoJ+1)*kSum;
111 G4double clebsch = ClebschGordanCoeff(twoJ1, twoM1, twoJ2, twoM2, twoJ);
112 return clebsch*clebsch;
115 std::vector<G4double>
120 std::vector<G4double> temp;
125 if (twoJ1 == 0 && twoJ2 == 0) {
126 G4Exception(
"G4Clebsch::GenerateIso3()",
"Clebsch010",
133 G4int twoM3 = twoM1 + twoM2;
138 temp.push_back(twoM3);
142 temp.push_back(twoM3);
148 G4int twoJMinIn =
std::max(std::abs(twoJ1 - twoJ2), std::abs(twoM3));
149 G4int twoJMaxIn = twoJ1 + twoJ2;
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;
159 twoJMinOut =
std::max(twoJMinOut, std::abs(twoM3));
160 G4int twoJMaxOut = twoJOut1 + twoJOut2;
165 if (twoJMin > twoJMax) {
166 G4Exception(
"G4Clebsch::GenerateIso3()",
"Clebsch020",
172 G4int nJ = (twoJMax - twoJMin) / 2 + 1;
176 if ( (twoJ1 == 0 || twoJ2 == 0) && twoJMin != twoJMax ) {
177 G4Exception(
"G4Clebsch::GenerateIso3()",
"Clebsch021",
178 JustWarning,
"twoJ1 or twoJ2 = 0, but twoJMin != JMax");
184 G4Exception(
"G4Clebsch::GenerateIso3()",
"Clebsch022",
185 JustWarning,
"nJ is zero, no overlap between in and out");
192 std::vector<G4double> clebsch;
194 for(
G4int twoJ=twoJMin; twoJ<=twoJMax; twoJ+=2) {
195 sum += ClebschGordan(twoJ1, twoM1, twoJ2, twoM2, twoJ);
196 clebsch.push_back(sum);
200 if (static_cast<G4int>(clebsch.size()) != nJ) {
201 G4Exception(
"G4Clebsch::GenerateIso3()",
"Clebsch023",
208 G4Exception(
"G4Clebsch::GenerateIso3()",
"Clebsch024",
209 JustWarning,
"Sum of Clebsch-Gordan probabilities <=0");
215 G4int twoJTot = twoJMin;
216 for (
G4int i=0; i<nJ; ++i) {
217 if (sum < clebsch[i]) {
225 std::vector<G4double> mMin;
226 mMin.push_back(-twoJOut1);
227 mMin.push_back(-twoJOut2);
229 std::vector<G4double> mMax;
230 mMax.push_back(twoJOut1);
231 mMax.push_back(twoJOut2);
235 std::vector<G4double> m1Out;
236 std::vector<G4double> m2Out;
238 const G4int size = 20;
241 G4int m1pos(0), m2pos(0);
243 G4int m1pr(0), m2pr(0);
246 for(j12 = std::abs(twoJOut1-twoJOut2); j12<=(twoJOut1+twoJOut2); j12+=2)
249 for (m1pr = static_cast<G4int>(mMin[0]+.00001); m1pr <= mMax[0]; m1pr+=2)
253 G4Exception(
"G4Clebsch::GenerateIso3()",
"Clebsch025",
257 m1Out.push_back(m1pr);
259 for (m2pr = static_cast<G4int>(mMin[1]+.00001); m2pr <= mMax[1]; m2pr+=2)
264 G4Exception(
"G4Clebsch::GenerateIso3()",
"Clebsch026",
268 m2Out.push_back(m2pr);
270 if(m1pr + m2pr == twoM3)
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);
277 prbout[m1pos][m2pos] = cleb;
282 prbout[m1pos][m2pos] = 0.;
289 G4Exception(
"G4Clebsch::GenerateIso3()",
"Clebsch027",
294 for (
G4int i=0; i<size; i++) {
295 for (
G4int j=0; j<size; j++) {
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]);
311 else rand -= prbout[m1p][m2p];
315 G4Exception(
"G4Clebsch::GenerateIso3()",
"Clebsch028",
326 G4int twoM = twoM1 + twoM2;
328 G4int twoJMinIn =
std::max(std::abs(twoJ1 - twoJ2), std::abs(twoM));
329 G4int twoJMaxIn = twoJ1 + twoJ2;
331 G4int twoJMinOut =
std::max(std::abs(twoJOut1 - twoJOut2), std::abs(twoM));
332 G4int twoJMaxOut = twoJOut1 + twoJOut2;
337 for (
G4int twoJ=twoJMin; twoJ<=twoJMax; twoJ+=2) {
339 value += ClebschGordan(twoJ1, twoM1, twoJ2, twoM2, twoJ);
356 return Wigner3J(twoJ1, twoM1, twoJ2, twoM2, twoJ3, twoM3);
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);
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);
388 if(twoJ1 == 0 || twoJ2 == 0)
return cleb;
392 for(
G4int twoM1Current=-twoJ1; twoM1Current<=twoJ1; twoM1Current+=2) {
393 G4int twoM2Current = twoM - twoM1Current;
395 G4double prob = ClebschGordan(twoJ1, twoM1Current, twoJ2,
398 if (twoM2Current == twoM2 && twoM1Current == twoM1) cleb += prob;
402 if (sum > 0.) cleb /= sum;
414 G4int i = twoA+twoB-twoC;
416 if(i<0 || (i%2))
return 0;
427 i = twoA+twoB+twoC+2;
435 if(twoJ1 < 0 || twoJ2 < 0 || twoJ3 < 0 ||
436 twoJ4 < 0 || twoJ5 < 0 || twoJ6 < 0)
return 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));
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);
458 double triangles = 0;
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);
482 G4int sum1 = (twoJ1 + twoJ2 + twoJ3)/2;
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;
492 G4int sum5 = (twoJ1 + twoJ2 + twoJ4 + twoJ5)/2;
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;
530 return triangles*kSum;
537 if(twoJ1 < 0 || twoJ2 < 0 || twoJ3 < 0 ||
538 twoJ4 < 0 || twoJ5 < 0 || twoJ6 < 0 ||
539 twoJ7 < 0 || twoJ8 < 0 || twoJ9 < 0)
return 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));
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;
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);
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;
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;
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;
612 if(twoM < -twoJ || twoM > twoJ || twoN < -twoJ || twoN > twoJ
613 || ((twoM % 2) != (twoJ % 2)) || ((twoN % 2) != (twoJ % 2)))
616 if(cosTheta == 1.0) {
return G4double(twoM == twoN); }
619 if(twoM > twoN) kMin = (twoM-twoN)/2;
621 if((twoJ-twoN)/2 <
kMax) kMax = (twoJ-twoN)/2;
637 logSum += (twoJ+(twoM-twoN)/2 - 2*k)*lnCosHalfTheta +
638 (2*k + (twoN-twoM)/2)*lnSinHalfTheta;
640 d += sign *
G4Exp(logSum);
static G4Pow * GetInstance()
static G4double Weight(G4int twoJ1, G4int twoM1, G4int twoJ2, G4int twoM2, G4int twoJOut1, G4int twoJOut2)
const G4int G4POWLOGFACTMAX
static G4double ClebschGordanCoeff(G4int twoJ1, G4int twoM1, G4int twoJ2, G4int twoM2, G4int twoJ)
static G4double TriangleCoeff(G4int twoA, G4int twoB, G4int twoC)
static G4double NormalizedClebschGordan(G4int twoJ, G4int twom, G4int twoJ1, G4int twoJ2, G4int twom1, G4int twom2)
static G4double Wigner9J(G4int twoJ1, G4int twoJ2, G4int twoJ3, G4int twoJ4, G4int twoJ5, G4int twoJ6, G4int twoJ7, G4int twoJ8, G4int twoJ9)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
static std::vector< G4double > GenerateIso3(G4int twoJ1, G4int twoM1, G4int twoJ2, G4int twoM2, G4int twoJOut1, G4int twoJOut2)
G4double G4Log(G4double x)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
static G4double Wigner3J(G4double j1, G4double j2, G4double j3, G4double m1, G4double m2, G4double m3)
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
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)
static G4double Wigner6J(G4int twoJ1, G4int twoJ2, G4int twoJ3, G4int twoJ4, G4int twoJ5, G4int twoJ6)
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)