41 for (i=1; i<nLogs; i++)
45 logs.push_back(value);
79 G4int an_m = iso3In1 + iso3In2;
81 G4int jMinIn = std::max(std::abs(isoIn1 - isoIn2), std::abs(an_m));
82 G4int jMaxIn = isoIn1 + isoIn2;
84 G4int jMinOut = std::max(std::abs(isoOut1 - isoOut2), std::abs(an_m));
85 G4int jMaxOut = isoOut1 + isoOut2;
87 G4int jMin = std::max(jMinIn,jMinOut);
88 G4int jMax = std::min(jMaxIn,jMaxOut);
91 for (j=jMin; j<=jMax; j+=2)
117 throw G4HadronicException(__FILE__, __LINE__,
"G4Clebsch::ClebschGordan - sqrt of negative argument");
118 G4double coeff = std::sqrt(argument) / (std::pow(-1.,n));
139 std::vector<G4double> n;
140 n.push_back(-j1 + j2 + j3);
141 n.push_back(j1 - m_1);
142 n.push_back(j1 + m_1);
143 n.push_back(j1 - j2 + j3);
144 n.push_back(j2 - m_2);
145 n.push_back(j2 + m_2);
146 n.push_back(j1 + j2 - j3);
147 n.push_back(j3 - m_3);
148 n.push_back(j3 + m_3);
156 G4double sum1 = n[i-1] + n[i+2] + n[i+5];
157 G4double sum2 = n[3*i-1] + n[3*i-2] + n[3*i-3];
158 if (sum1 != sigma || sum2 != sigma) ok =
false;
162 if (n[i+3*j-4] < 0.) ok =
false;
178 if (n[i+3*j-4] < smallest)
180 smallest = n[i+3*j-4];
191 for(
G4int j=1; j<=3; ++j)
194 n[j*3-3] = n[iMin+j*3-4];
197 sign = (
G4int) std::pow(-1.,sigma);
205 n[i-1] = n[i+jMin*3-4];
208 sign *= (
G4int) std::pow(-1.,sigma);
211 const std::vector<G4double> logVector =
GetLogs();
215 G4int logEntries = logVector.size() - 1;
216 for (i=0; i<n.size(); i++)
218 if (n[i] < 0. || n[i] > logEntries)
219 throw G4HadronicException(__FILE__, __LINE__,
"G4Clebsch::Wigner3J - Outside logVector boundaries, n");
243 if (sigma1 < 0. || sigma1 > logEntries)
244 throw G4HadronicException(__FILE__, __LINE__,
"G4Clebsch::Wigner3J - Outside logVector boundaries, sigma");
247 G4double hlp1 = (l2 + l3 + l4 +l7 -ls -l1 -l5 -l9 -l6 -l8) / 2.;
248 G4int expon =
static_cast<G4int>(r6 + r8+.00001);
249 G4double sgn = std::pow(-1., expon);
250 G4double coeff = std::exp(hlp1) * sgn;
252 G4int n61 =
static_cast<G4int>(r6 - r1+.00001);
253 if (n61 < 0. || n61 > logEntries)
254 throw G4HadronicException(__FILE__, __LINE__,
"G4Clebsch::Wigner3J - Outside logVector boundaries, n61");
255 G4int n81 =
static_cast<G4int>(r8 - r1+.00001);
256 if (n81 < 0. || n81 > logEntries)
257 throw G4HadronicException(__FILE__, __LINE__,
"G4Clebsch::Wigner3J - Outside logVector boundaries, n81");
259 G4double hlp2 = l6 - logVector[n61] + l8 - logVector[n81];
261 std::vector<G4double> S;
264 for (i=1; i<=n1; i++)
267 G4double den = i * (r6 - r1 + i) * (r8 - r1 + i);
270 G4double data = -last * (r1 + 1.0 - i) * (r5 + 1.0 - i) * (r9 + 1. - i) / den;
274 value = coeff * sum * sign;
296 std::vector<G4double> temp;
301 if (isoIn1 == 0 && isoIn2 == 0)
303 G4cout <<
"WARNING: G4Clebsch::GenerateIso3 - both isoIn are zero" <<
G4endl;
309 G4int iso3 = iso3In1 + iso3In2;
315 temp.push_back(iso3);
320 temp.push_back(iso3);
326 G4int jMinIn = std::max(std::abs(isoIn1 - isoIn2), std::abs(iso3));
327 G4int jMaxIn = isoIn1 + isoIn2;
331 G4int jMinOut = 9999;
334 for(i=-1; i<=1; i+=2)
336 for(j=-1; j<=1; j+=2)
338 jTmp= std::abs(i*isoA + j*isoB);
339 if(jTmp < jMinOut) jMinOut = jTmp;
342 jMinOut = std::max(jMinOut, std::abs(iso3));
343 G4int jMaxOut = isoA + isoB;
346 G4int jMin = std::max(jMinIn, jMinOut);
347 G4int jMax = std::min(jMaxIn, jMaxOut);
354 G4int nJ = (jMax - jMin) / 2 + 1;
358 if ( (isoIn1 == 0 || isoIn2 == 0) && jMin != jMax )
359 throw G4HadronicException(__FILE__, __LINE__,
"G4Clebsch::GenerateIso3 - J1 or J2 = 0, but jMin != JMax");
363 throw G4HadronicException(__FILE__, __LINE__,
"G4Clebsch::GenerateIso3 - nJ is zero, no overlap between in and out");
368 std::vector<G4double> clebsch;
370 for(j=jMin; j<=jMax; j+=2)
373 clebsch.push_back(cg);
377 if (
static_cast<G4int>(clebsch.size()) != nJ)
378 throw G4HadronicException(__FILE__, __LINE__,
"G4Clebsch::GenerateIso3 - nJ inconsistency");
388 throw G4HadronicException(__FILE__, __LINE__,
"G4Clebsch::GenerateIso3 - Sum of Clebsch-Gordan probabilities <=0");
392 std::vector<G4double> clebschPdf;
394 clebschPdf.push_back(previous/sum);
397 previous += clebsch[j];
399 clebschPdf.push_back(prob);
408 if (rand < clebschPdf[j])
418 std::vector<G4double> mMin;
419 mMin.push_back(-isoA);
420 mMin.push_back(-isoB);
422 std::vector<G4double> mMax;
423 mMax.push_back(isoA);
424 mMax.push_back(isoB);
428 std::vector<G4double> m1Out;
429 std::vector<G4double> m2Out;
431 const G4int size = 20;
434 G4int m1pos(0), m2pos(0);
436 G4int m1pr(0), m2pr(0);
439 for(j12 = std::abs(isoA-isoB); j12<=(isoA+isoB); j12+=2)
442 for (m1pr =
static_cast<G4int>(mMin[0]+.00001); m1pr <= mMax[0]; m1pr+=2)
447 m1Out.push_back(m1pr);
449 for (m2pr =
static_cast<G4int>(mMin[1]+.00001); m2pr <= mMax[1]; m2pr+=2)
456 m2Out.push_back(m2pr);
458 if(m1pr + m2pr == iso3)
460 G4int m12 = m1pr + m2pr;
465 prbout[m1pos][m2pos] = cleb;
470 prbout[m1pos][m2pos] = 0.;
479 for (i=0; i<size; i++)
481 for (j=0; j<size; j++)
491 for (m1p=0; m1p<m1pos; m1p++)
493 for (m2p=0; m2p<m2pos; m2p++)
495 if (rand < prbout[m1p][m2p])
497 temp.push_back(m1Out[m1p]);
498 temp.push_back(m2Out[m2p]);
503 rand -= prbout[m1p][m2p];
522 if(J1 == 0 || J2 == 0)
return cleb;
528 for(
G4int m1Current=-J1; m1Current<=J1; m1Current+=2)
530 G4int m2Current = M - m1Current;
534 if (m2Current == m_2 && m1Current == m_1) cleb += prob;
538 if (sum > 0.) cleb /= sum;
G4DLLIMPORT std::ostream G4cout
const std::vector< G4double > & GetLogs() const
G4double ClebschGordan(G4int isoIn1, G4int iso3In1, G4int isoIn2, G4int iso3In2, G4int jOut) const
G4bool operator==(const G4Clebsch &right) const
std::vector< G4double > GenerateIso3(G4int isoIn1, G4int iso3In1, G4int isoIn2, G4int iso3In2, G4int isoOut1, G4int isoOut2) const
G4double Weight(G4int isoIn1, G4int iso3In1, G4int isoIn2, G4int iso3In2, G4int isoOut1, G4int isoOut2) const
G4bool operator!=(const G4Clebsch &right) const
G4double NormalizedClebschGordan(G4int J, G4int m, G4int J1, G4int J2, G4int m1, G4int m2) const
G4double Wigner3J(G4double j1, G4double j2, G4double j3, G4double m1, G4double m2, G4double m3) const