G4Clebsch Class Reference

#include <G4Clebsch.hh>

Static Public Member Functions
static G4double	ClebschGordan (G4int twoJ1, G4int twoM1, G4int twoJ2, G4int twoM2, G4int twoJ)
static G4double	ClebschGordanCoeff (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	NormalizedClebschGordan (G4int twoJ, G4int twom, G4int twoJ1, G4int twoJ2, G4int twom1, G4int twom2)
static G4double	RacahWCoeff (G4int twoJ1, G4int twoJ2, G4int twoJ, G4int twoJ3, G4int twoJ12, G4int twoJ23)
static G4double	TriangleCoeff (G4int twoA, G4int twoB, G4int twoC)
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	Wigner6J (G4int twoJ1, G4int twoJ2, G4int twoJ3, G4int twoJ4, G4int twoJ5, G4int twoJ6)
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.

{
  // ClebschGordanCoeff() will do all input checking
  G4double clebsch = ClebschGordanCoeff(twoJ1, twoM1, twoJ2, twoM2, twoJ);
  return clebsch*clebsch;
}

References ClebschGordanCoeff().

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.

{
  if(twoJ1 < 0 || twoJ2 < 0 || twoJ < 0 ||
     ((twoJ1-twoM1) % 2) || ((twoJ2-twoM2) % 2)) { return 0; }
 
  G4int twoM = twoM1 + twoM2;
  if(twoM1 > twoJ1 || twoM1 < -twoJ1 ||
     twoM2 > twoJ2 || twoM2 < -twoJ2 ||
     twoM > twoJ || twoM < -twoJ) { return 0; }
 
  // Checks limits on J1, J2, J3
  G4double triangle = TriangleCoeff(twoJ1, twoJ2, twoJ);
  if(triangle == 0) { return 0; }
 
  G4Pow* g4pow = G4Pow::GetInstance();
  G4double factor = g4pow->logfactorial((twoJ1 + twoM1)/2) +
                    g4pow->logfactorial((twoJ1 - twoM1)/2);
  factor += g4pow->logfactorial((twoJ2 + twoM2)/2) +
            g4pow->logfactorial((twoJ2 - twoM2)/2);
  factor += g4pow->logfactorial((twoJ + twoM)/2) +
            g4pow->logfactorial((twoJ - twoM)/2);
  factor *= 0.5;
 
  G4int kMin = 0;
  G4int sum1 = (twoJ1 - twoM1)/2;
  G4int kMax = sum1;
  G4int sum2 = (twoJ - twoJ2 + twoM1)/2;
  if(-sum2 > kMin) kMin = -sum2;
  G4int sum3 = (twoJ2 + twoM2)/2;
  if(sum3 < kMax) kMax = sum3;
  G4int sum4 = (twoJ - twoJ1 - twoM2)/2;
  if(-sum4 > kMin) kMin = -sum4;
  G4int sum5 = (twoJ1 + twoJ2 - twoJ)/2;
  if(sum5 < kMax) kMax = sum5;
 
  // sanity / boundary checks
  if(kMin < 0) {
    G4Exception("G4Clebsch::ClebschGordanCoeff()", "Clebsch001", 
        JustWarning, "kMin < 0");
    return 0;
  }
  if(kMax < kMin) {
    G4Exception("G4Clebsch::ClebschGordanCoeff()", "Clebsch002", 
        JustWarning, "kMax < kMin");
    return 0;
  }
  if(kMax >= G4POWLOGFACTMAX) {
    G4Exception("G4Clebsch::ClebschGordanCoeff()", "Clebsch003", 
        JustWarning, "kMax too big for G4Pow");
    return 0;
  }
 
  // Now do the sum over k
  G4double kSum = 0.;
  for(G4int k = kMin; k <= kMax; k++) {
    G4double sign = (k % 2) ? -1 : 1;
    kSum += sign * G4Exp(factor - g4pow->logfactorial(sum1-k) -
             g4pow->logfactorial(sum2+k) -
             g4pow->logfactorial(sum3-k) -
             g4pow->logfactorial(sum4+k) -
             g4pow->logfactorial(k) -
             g4pow->logfactorial(sum5-k));
  }
 
  return triangle*sqrt(twoJ+1)*kSum;
}

References G4Exception(), G4Exp(), G4POWLOGFACTMAX, G4Pow::GetInstance(), JustWarning, BoundingBox::kMax, BoundingBox::kMin, G4Pow::logfactorial(), G4INCL::Math::sign(), and TriangleCoeff().

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.

{
  std::vector<G4double> temp;
 
  // ---- Special cases first ----
 
  // Special case, both Jin are zero
  if (twoJ1 == 0 && twoJ2 == 0) {
    G4Exception("G4Clebsch::GenerateIso3()", "Clebsch010", 
        JustWarning, "both twoJ are zero");
    temp.push_back(0.);
    temp.push_back(0.);
    return temp;
  }
 
  G4int twoM3 = twoM1 + twoM2;
 
  // Special case, either Jout is zero
  if (twoJOut1 == 0) {  
    temp.push_back(0.);
    temp.push_back(twoM3);
    return temp;
  }
  if (twoJOut2 == 0) {
    temp.push_back(twoM3);
    temp.push_back(0.);
    return temp;
  }
  
  // Number of possible states, in 
  G4int twoJMinIn = std::max(std::abs(twoJ1 - twoJ2), std::abs(twoM3));
  G4int twoJMaxIn = twoJ1 + twoJ2;
 
  // Number of possible states, out
  G4int twoJMinOut = 9999;
  for(G4int i=-1; i<=1; i+=2) {
    for(G4int j=-1; j<=1; j+=2) {
      G4int twoJTmp= std::abs(i*twoJOut1 + j*twoJOut2);
      if(twoJTmp < twoJMinOut) twoJMinOut = twoJTmp;
    }
  }
  twoJMinOut = std::max(twoJMinOut, std::abs(twoM3));
  G4int twoJMaxOut = twoJOut1 + twoJOut2;
 
  // Possible in and out common states 
  G4int twoJMin  =  std::max(twoJMinIn, twoJMinOut);
  G4int twoJMax  =  std::min(twoJMaxIn, twoJMaxOut);
  if (twoJMin > twoJMax) {
    G4Exception("G4Clebsch::GenerateIso3()", "Clebsch020", 
        JustWarning, "twoJMin > twoJMax");
    return temp;
  }
  
  // Number of possible isospins
  G4int nJ = (twoJMax - twoJMin) / 2 + 1;
 
  // A few consistency checks
  
  if ( (twoJ1 == 0 || twoJ2 == 0) && twoJMin != twoJMax ) {
    G4Exception("G4Clebsch::GenerateIso3()", "Clebsch021", 
        JustWarning, "twoJ1 or twoJ2 = 0, but twoJMin != JMax");
    return temp;
  }
 
  // MGP ---- Shall it be a warning or an exception?
  if (nJ == 0) {
    G4Exception("G4Clebsch::GenerateIso3()", "Clebsch022", 
        JustWarning, "nJ is zero, no overlap between in and out");
    return temp;
  }
 
  // Loop over all possible combinations of twoJ1, twoJ2, twoM11, twoM2, twoJTot
  // to get the probability of each of the in-channel couplings
 
  std::vector<G4double> clebsch;
  G4double sum = 0.0;
  for(G4int twoJ=twoJMin; twoJ<=twoJMax; twoJ+=2) {
    sum += ClebschGordan(twoJ1, twoM1, twoJ2, twoM2, twoJ);
    clebsch.push_back(sum);
  }     
 
  // Consistency check
  if (static_cast<G4int>(clebsch.size()) != nJ) {
    G4Exception("G4Clebsch::GenerateIso3()", "Clebsch023", 
        JustWarning, "nJ inconsistency");
    return temp;
  }
 
  // Consistency check
  if (sum <= 0.) {
    G4Exception("G4Clebsch::GenerateIso3()", "Clebsch024", 
        JustWarning, "Sum of Clebsch-Gordan probabilities <=0");
    return temp;
  }
 
  // Generate a random twoJTot according to the Clebsch-Gordan pdf
  sum *= G4UniformRand();
  G4int twoJTot = twoJMin;
  for (G4int i=0; i<nJ; ++i) {
    if (sum < clebsch[i]) {
      twoJTot += 2*i;
      break;
    }
  }
 
  // Generate twoM3Out
 
  std::vector<G4double> mMin;
  mMin.push_back(-twoJOut1);
  mMin.push_back(-twoJOut2);
 
  std::vector<G4double> mMax;
  mMax.push_back(twoJOut1);
  mMax.push_back(twoJOut2);
 
  // Calculate the possible |J_i M_i> combinations and their probability
 
  std::vector<G4double> m1Out;
  std::vector<G4double> m2Out;
 
  const G4int size = 20;
  G4double prbout[size][size];
 
  G4int m1pos(0), m2pos(0);
  G4int j12;
  G4int m1pr(0), m2pr(0);
 
  sum = 0.;
  for(j12 = std::abs(twoJOut1-twoJOut2); j12<=(twoJOut1+twoJOut2); j12+=2)
  {
    m1pos = -1;
    for (m1pr = static_cast<G4int>(mMin[0]+.00001); m1pr <= mMax[0]; m1pr+=2)
    {
      m1pos++;
      if (m1pos >= size) {
        G4Exception("G4Clebsch::GenerateIso3()", "Clebsch025", 
            JustWarning, "m1pos > size");
        return temp;
      }
      m1Out.push_back(m1pr);
      m2pos = -1;
      for (m2pr = static_cast<G4int>(mMin[1]+.00001); m2pr <= mMax[1]; m2pr+=2)
      {
        m2pos++;
        if (m2pos >= size)
        {
          G4Exception("G4Clebsch::GenerateIso3()", "Clebsch026", 
              JustWarning, "m2pos > size");
          return temp;
        }
        m2Out.push_back(m2pr);
 
        if(m1pr + m2pr == twoM3) 
        {
          G4int m12 = m1pr + m2pr;
          G4double c12 = ClebschGordan(twoJOut1, m1pr, twoJOut2,m2pr, j12);
          G4double c34 = ClebschGordan(0,0,0,0,0);
          G4double ctot = ClebschGordan(j12, m12, 0, 0, twoJTot);
          G4double cleb = c12*c34*ctot;
          prbout[m1pos][m2pos] = cleb;
          sum += cleb;
        }
        else
        {
          prbout[m1pos][m2pos] = 0.;
        }
      }
    }
  }
  
  if (sum <= 0.) {
    G4Exception("G4Clebsch::GenerateIso3()", "Clebsch027", 
        JustWarning, "sum (out) <=0");
    return temp;
  }
 
  for (G4int i=0; i<size; i++) {
    for (G4int j=0; j<size; j++) {
      prbout[i][j] /= sum;
    }
  }
 
  G4double rand = G4UniformRand();
 
  G4int m1p, m2p;
 
  for (m1p=0; m1p<m1pos; m1p++) {
    for (m2p=0; m2p<m2pos; m2p++) {
      if (rand < prbout[m1p][m2p]) {
        temp.push_back(m1Out[m1p]);
        temp.push_back(m2Out[m2p]);
        return temp;
      }   
      else rand -= prbout[m1p][m2p];
    }     
  }   
 
  G4Exception("G4Clebsch::GenerateIso3()", "Clebsch028", 
          JustWarning, "Should never get here");
  return temp;
}

References ClebschGordan(), G4Exception(), G4UniformRand, JustWarning, G4INCL::Math::max(), and G4INCL::Math::min().

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.

{
  // Calculate the normalized Clebsch-Gordan coefficient, that is the prob 
  // of isospin decomposition of (J,m) into J1, J2, m1, m2
 
  G4double cleb = 0.;
  if(twoJ1 == 0 || twoJ2 == 0) return cleb; 
  
  // Loop over all J1,J2,Jtot,m1,m2 combinations
  G4double sum = 0.0;
  for(G4int twoM1Current=-twoJ1; twoM1Current<=twoJ1;  twoM1Current+=2) {
    G4int twoM2Current = twoM - twoM1Current;
    // ClebschGordan() will do all further input checking
    G4double prob = ClebschGordan(twoJ1, twoM1Current, twoJ2, 
                  twoM2Current, twoJ);
    sum += prob;
    if (twoM2Current == twoM2 && twoM1Current == twoM1) cleb += prob;
  }
 
  // Normalize probs to 1 
  if (sum > 0.) cleb /= sum; 
 
  return cleb;
}

References ClebschGordan().

Referenced by G4XAnnihilationChannel::NormalizedClebsch().

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.

{
  // TC(ABC) = sqrt[ (A+B-C)! (A-B+C)! (-A+B+C)! / (A+B+C+1)! ]
  // return 0 if the triad does not satisfy the triangle inequalities
  G4Pow* g4pow =  G4Pow::GetInstance();
 
  double val = 0;
  G4int i = twoA+twoB-twoC;
  // only have to check that i is even the first time
  if(i<0 || (i%2)) return 0;
  else val += g4pow->logfactorial(i/2);
 
  i = twoA-twoB+twoC;
  if(i<0) return 0;
  else val += g4pow->logfactorial(i/2);
 
  i = -twoA+twoB+twoC;
  if(i<0) return 0;
  else val += g4pow->logfactorial(i/2);
 
  i = twoA+twoB+twoC+2;
  if(i<0) return 0;
  return G4Exp(0.5*(val - g4pow->logfactorial(i/2)));
}

References G4Exp(), G4Pow::GetInstance(), and G4Pow::logfactorial().

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.

{
  G4double value = 0.;
  
  G4int twoM = twoM1 + twoM2;
 
  G4int twoJMinIn = std::max(std::abs(twoJ1 - twoJ2), std::abs(twoM));
  G4int twoJMaxIn = twoJ1 + twoJ2;
 
  G4int twoJMinOut = std::max(std::abs(twoJOut1 - twoJOut2), std::abs(twoM));
  G4int twoJMaxOut = twoJOut1 + twoJOut2;
 
  G4int twoJMin = std::max(twoJMinIn,twoJMinOut);
  G4int twoJMax = std::min(twoJMaxIn,twoJMaxOut);
 
  for (G4int twoJ=twoJMin; twoJ<=twoJMax; twoJ+=2) {
    // ClebschGordan() will do all input checking
    value += ClebschGordan(twoJ1, twoM1, twoJ2, twoM2, twoJ);
  }
 
  return value;
}

References ClebschGordan(), G4INCL::Math::max(), and G4INCL::Math::min().

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.

{
  //  G4Exception("G4Clebsch::Wigner3J()", "Clebsch030", JustWarning, 
  //  "G4Clebsch::Wigner3J with double arguments is deprecated. Please use G4int version.");
  G4int twoJ1 = (G4int) (2.*j1);
  G4int twoJ2 = (G4int) (2.*j2);
  G4int twoJ3 = (G4int) (2.*j3);
  G4int twoM1 = (G4int) (2.*m1);
  G4int twoM2 = (G4int) (2.*m2);
  G4int twoM3 = (G4int) (2.*m3);
  return Wigner3J(twoJ1, twoM1, twoJ2, twoM2, twoJ3, twoM3);
}

References m2, m3, and Wigner3J().

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.

{
  G4double clebsch = ClebschGordanCoeff(twoJ1, twoM1, twoJ2, twoM2, twoJ3);
  if(clebsch == 0) return clebsch;
  if( (twoJ1-twoJ2+twoM1+twoM2)/2 % 2) clebsch = -clebsch;
  return clebsch / sqrt(twoJ3+1);
}

References ClebschGordanCoeff().

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.

{
  if(twoM1 + twoM2 != -twoM3) return 0;
  G4double clebsch = ClebschGordanCoeff(twoJ1, twoM1, twoJ2, twoM2, twoJ3);
  if(clebsch == 0) return clebsch;
  if( (twoJ1-twoJ2-twoM3)/2 % 2) clebsch = -clebsch;
  return clebsch / sqrt(twoJ3+1);
}

References ClebschGordanCoeff().

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.

{
  if(twoJ1 < 0 || twoJ2 < 0 || twoJ3 < 0 || 
     twoJ4 < 0 || twoJ5 < 0 || twoJ6 < 0) return 0;
 
  // There is a fast calculation (no sums or exps) when twoJ6 = 0, 
  // so permute to use it when possible
  if(twoJ6 == 0) {
    if(twoJ1 != twoJ5) return 0;
    if(twoJ2 != twoJ4) return 0;
    if(twoJ1+twoJ2 < twoJ3) return 0;
    if((twoJ1 > twoJ2) && (twoJ3 < (twoJ1-twoJ2))) return 0;
    if((twoJ2 > twoJ1) && (twoJ3 < (twoJ2-twoJ1))) return 0;
    if((twoJ1+twoJ2+twoJ3) % 2) return 0;
    return (((twoJ1+twoJ2+twoJ3)/2) % 2 ? -1. : 1.) /sqrt((twoJ1+1)*(twoJ2+1));
  }
  if(twoJ1 == 0) return Wigner6J(twoJ6, twoJ2, twoJ4, twoJ3, twoJ5, 0);
  if(twoJ2 == 0) return Wigner6J(twoJ1, twoJ6, twoJ5, twoJ4, twoJ3, 0);
  if(twoJ3 == 0) return Wigner6J(twoJ4, twoJ2, twoJ6, twoJ1, twoJ5, 0);
  if(twoJ4 == 0) return Wigner6J(twoJ3, twoJ2, twoJ1, twoJ6, twoJ5, 0);
  if(twoJ5 == 0) return Wigner6J(twoJ1, twoJ3, twoJ2, twoJ4, twoJ6, 0);
 
  // Check triangle inequalities and calculate triangle coefficients.
  // Also check evenness of sums
  G4Pow* g4pow =  G4Pow::GetInstance();
  double triangles = 0;
  G4int i;
  i =  twoJ1+twoJ2-twoJ3;   if(i<0 || i%2) return 0; else triangles += g4pow->logfactorial(i/2);
  i =  twoJ1-twoJ2+twoJ3;   if(i<0 || i%2) return 0; else triangles += g4pow->logfactorial(i/2);
  i = -twoJ1+twoJ2+twoJ3;   if(i<0 || i%2) return 0; else triangles += g4pow->logfactorial(i/2);
  i =  twoJ1+twoJ2+twoJ3+2; if(i<0 || i%2) return 0; else triangles -= g4pow->logfactorial(i/2);
  i =  twoJ1+twoJ5-twoJ6;   if(i<0 || i%2) return 0; else triangles += g4pow->logfactorial(i/2);
  i =  twoJ1-twoJ5+twoJ6;   if(i<0 || i%2) return 0; else triangles += g4pow->logfactorial(i/2);
  i = -twoJ1+twoJ5+twoJ6;   if(i<0 || i%2) return 0; else triangles += g4pow->logfactorial(i/2);
  i =  twoJ1+twoJ5+twoJ6+2; if(i<0 || i%2) return 0; else triangles -= g4pow->logfactorial(i/2);
  i =  twoJ4+twoJ2-twoJ6;   if(i<0 || i%2) return 0; else triangles += g4pow->logfactorial(i/2);
  i =  twoJ4-twoJ2+twoJ6;   if(i<0 || i%2) return 0; else triangles += g4pow->logfactorial(i/2);
  i = -twoJ4+twoJ2+twoJ6;   if(i<0 || i%2) return 0; else triangles += g4pow->logfactorial(i/2);
  i =  twoJ4+twoJ2+twoJ6+2; if(i<0 || i%2) return 0; else triangles -= g4pow->logfactorial(i/2);
  i =  twoJ4+twoJ5-twoJ3;   if(i<0 || i%2) return 0; else triangles += g4pow->logfactorial(i/2);
  i =  twoJ4-twoJ5+twoJ3;   if(i<0 || i%2) return 0; else triangles += g4pow->logfactorial(i/2);
  i = -twoJ4+twoJ5+twoJ3;   if(i<0 || i%2) return 0; else triangles += g4pow->logfactorial(i/2);
  i =  twoJ4+twoJ5+twoJ3+2; if(i<0 || i%2) return 0; else triangles -= g4pow->logfactorial(i/2);
  triangles = G4Exp(0.5*triangles);
  
  // Prepare to sum over k. If we have made it this far, all of the following
  // sums must be non-negative and divisible by two
 
  // k must be >= all of the following sums:
  G4int sum1 = (twoJ1 + twoJ2 + twoJ3)/2;
  G4int kMin = sum1;
  G4int sum2 = (twoJ1 + twoJ5 + twoJ6)/2;
  if(sum2 > kMin) kMin = sum2;
  G4int sum3 = (twoJ4 + twoJ2 + twoJ6)/2;
  if(sum3 > kMin) kMin = sum3;
  G4int sum4 = (twoJ4 + twoJ5 + twoJ3)/2;
  if(sum4 > kMin) kMin = sum4;
 
  // and k must be <= all of the following sums:
  G4int sum5 = (twoJ1 + twoJ2 + twoJ4 + twoJ5)/2;
  G4int kMax = sum5;
  G4int sum6 = (twoJ2 + twoJ3 + twoJ5 + twoJ6)/2;
  if(sum6 < kMax) kMax = sum6;
  G4int sum7 = (twoJ1 + twoJ3 + twoJ4 + twoJ6)/2;
  if(sum7 < kMax) kMax = sum7;
 
  // sanity / boundary checks
  if(kMin < 0) {
    G4Exception("G4Clebsch::Wigner6J()", "Clebsch040", 
        JustWarning, "kMin < 0");
    return 0;
  }
  if(kMax < kMin) {
    G4Exception("G4Clebsch::Wigner6J()", "Clebsch041", 
        JustWarning, "kMax < kMin");
    return 0;
  }
  if(kMax >= G4POWLOGFACTMAX) {
    G4Exception("G4Clebsch::Wigner6J()", "Clebsch041", 
        JustWarning, "kMax too big for G4Pow");
    return 0;
  }
 
  // Now do the sum over k
  G4double kSum = 0.;
  G4double sign = (kMin % 2) ? -1 : 1;
  for(G4int k = kMin; k <= kMax; k++) {
    kSum += sign * G4Exp(g4pow->logfactorial(k+1) -
                            g4pow->logfactorial(k-sum1) -
                            g4pow->logfactorial(k-sum2) -
                            g4pow->logfactorial(k-sum3) -
                            g4pow->logfactorial(k-sum4) -
                            g4pow->logfactorial(sum5-k) -
                            g4pow->logfactorial(sum6-k) -
                            g4pow->logfactorial(sum7-k));
    sign *= -1;
  }
  return triangles*kSum;
}

References G4Exception(), G4Exp(), G4POWLOGFACTMAX, G4Pow::GetInstance(), JustWarning, BoundingBox::kMax, BoundingBox::kMin, G4Pow::logfactorial(), G4INCL::Math::sign(), and Wigner6J().

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.

{
  if(twoJ1 < 0 || twoJ2 < 0 || twoJ3 < 0 || 
     twoJ4 < 0 || twoJ5 < 0 || twoJ6 < 0 ||
     twoJ7 < 0 || twoJ8 < 0 || twoJ9 < 0) return 0;
 
  if(twoJ9 == 0) {
    if(twoJ3 != twoJ6) return 0;
    if(twoJ7 != twoJ8) return 0;
    G4double sixJ = Wigner6J(twoJ1, twoJ2, twoJ3, twoJ5, twoJ4, twoJ7);
    if(sixJ == 0) return 0;
    if((twoJ2+twoJ3+twoJ4+twoJ7)/2 % 2) sixJ = -sixJ;
    return sixJ/sqrt((twoJ3+1)*(twoJ7+1));
  }
  if(twoJ1 == 0) return Wigner9J(twoJ9, twoJ6, twoJ3, twoJ8, twoJ5, twoJ2, twoJ7, twoJ4, twoJ1);
  if(twoJ2 == 0) return Wigner9J(twoJ7, twoJ9, twoJ8, twoJ4, twoJ6, twoJ5, twoJ1, twoJ3, twoJ2);
  if(twoJ4 == 0) return Wigner9J(twoJ3, twoJ2, twoJ1, twoJ9, twoJ8, twoJ7, twoJ6, twoJ5, twoJ4);
  if(twoJ5 == 0) return Wigner9J(twoJ1, twoJ3, twoJ2, twoJ7, twoJ9, twoJ8, twoJ4, twoJ6, twoJ5);
  G4int twoS = twoJ1+twoJ2+twoJ3+twoJ4+twoJ5+twoJ6+twoJ7+twoJ8+twoJ9;
  if(twoS % 2) return 0;
  G4double sign = (twoS/2 % 2) ? -1 : 1;
  if(twoJ3 == 0) return sign*Wigner9J(twoJ7, twoJ8, twoJ9, twoJ4, twoJ5, twoJ6, twoJ1, twoJ2, twoJ3);
  if(twoJ6 == 0) return sign*Wigner9J(twoJ1, twoJ2, twoJ3, twoJ7, twoJ8, twoJ9, twoJ4, twoJ5, twoJ6);
  if(twoJ7 == 0) return sign*Wigner9J(twoJ3, twoJ2, twoJ1, twoJ6, twoJ5, twoJ4, twoJ9, twoJ8, twoJ7);
  if(twoJ8 == 0) return sign*Wigner9J(twoJ1, twoJ3, twoJ2, twoJ4, twoJ6, twoJ5, twoJ7, twoJ9, twoJ8);
 
  // No element is zero: check triads now for speed
  G4int i;
  i =  twoJ1+twoJ2-twoJ3; if(i<0 || i%2) return 0;
  i =  twoJ1-twoJ2+twoJ3; if(i<0 || i%2) return 0;
  i = -twoJ1+twoJ2+twoJ3; if(i<0 || i%2) return 0;
  i =  twoJ4+twoJ5-twoJ6; if(i<0 || i%2) return 0;
  i =  twoJ4-twoJ5+twoJ6; if(i<0 || i%2) return 0;
  i = -twoJ4+twoJ5+twoJ6; if(i<0 || i%2) return 0;
  i =  twoJ7+twoJ8-twoJ9; if(i<0 || i%2) return 0;
  i =  twoJ7-twoJ8+twoJ9; if(i<0 || i%2) return 0;
  i = -twoJ7+twoJ8+twoJ9; if(i<0 || i%2) return 0;
  i =  twoJ1+twoJ4-twoJ7; if(i<0 || i%2) return 0;
  i =  twoJ1-twoJ4+twoJ7; if(i<0 || i%2) return 0;
  i = -twoJ1+twoJ4+twoJ7; if(i<0 || i%2) return 0;
  i =  twoJ2+twoJ5-twoJ8; if(i<0 || i%2) return 0;
  i =  twoJ2-twoJ5+twoJ8; if(i<0 || i%2) return 0;
  i = -twoJ2+twoJ5+twoJ8; if(i<0 || i%2) return 0;
  i =  twoJ3+twoJ6-twoJ9; if(i<0 || i%2) return 0;
  i =  twoJ3-twoJ6+twoJ9; if(i<0 || i%2) return 0;
  i = -twoJ3+twoJ6+twoJ9; if(i<0 || i%2) return 0;
  
  // Okay, have to do the full sum over 6J's
  // Find limits for K sum
  G4int twoKMax = twoJ1+twoJ9;
  if(twoJ4+twoJ8 < twoKMax) twoKMax = twoJ4+twoJ8;
  if(twoJ2+twoJ6 < twoKMax) twoKMax = twoJ2+twoJ6;
  G4int twoKMin = twoJ1-twoJ9;
  if(twoJ9-twoJ1 > twoKMin) twoKMin = twoJ9-twoJ1;
  if(twoJ4-twoJ8 > twoKMin) twoKMin = twoJ4-twoJ8;
  if(twoJ8-twoJ4 > twoKMin) twoKMin = twoJ8-twoJ4;
  if(twoJ2-twoJ6 > twoKMin) twoKMin = twoJ2-twoJ6;
  if(twoJ6-twoJ2 > twoKMin) twoKMin = twoJ6-twoJ2;
  if(twoKMin > twoKMax) return 0;
 
  G4double sum = 0;
  for(G4int twoK = twoKMin; twoK <= twoKMax; twoK += 2) {
    G4double value = Wigner6J(twoJ1, twoJ4, twoJ7, twoJ8, twoJ9, twoK);
    if(value == 0) continue;
    value *= Wigner6J(twoJ2, twoJ5, twoJ8, twoJ4, twoK, twoJ6);
    if(value == 0) continue;
    value *= Wigner6J(twoJ3, twoJ6, twoJ9, twoK, twoJ1, twoJ2);
    if(value == 0) continue;
    if(twoK % 2) value = -value;
    sum += value*G4double(twoK+1);
  }
  return sum;
}

References G4INCL::Math::sign(), Wigner6J(), and Wigner9J().

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.

{
  if(twoM < -twoJ || twoM > twoJ || twoN < -twoJ || twoN > twoJ 
     || ((twoM % 2) != (twoJ % 2)) || ((twoN % 2) != (twoJ % 2))) 
    { return 0; }
 
     if(cosTheta == 1.0) { return G4double(twoM == twoN); }
 
  G4int kMin = 0;
  if(twoM > twoN) kMin = (twoM-twoN)/2;
  G4int kMax = (twoJ + twoM)/2;
  if((twoJ-twoN)/2 < kMax) kMax = (twoJ-twoN)/2;
 
  G4double lnCosHalfTheta = G4Log((cosTheta+1.)*0.5) * 0.5;
  G4double lnSinHalfTheta = G4Log((1.-cosTheta)*0.5) * 0.5;
 
  G4Pow* g4pow =  G4Pow::GetInstance();
  G4double d = 0;
  for(G4int k = kMin; k <= kMax; k++) {
    G4double logSum = 0.5*(g4pow->logfactorial((twoJ+twoM)/2) +
                           g4pow->logfactorial((twoJ-twoM)/2) +
                           g4pow->logfactorial((twoJ+twoN)/2) +
                           g4pow->logfactorial((twoJ-twoN)/2));
    logSum += -g4pow->logfactorial((twoJ+twoM)/2 - k) -
               g4pow->logfactorial((twoJ-twoN)/2 - k) -
               g4pow->logfactorial(k) - 
               g4pow->logfactorial(k+(twoN-twoM)/2);
    logSum += (twoJ+(twoM-twoN)/2 - 2*k)*lnCosHalfTheta + 
      (2*k + (twoN-twoM)/2)*lnSinHalfTheta;
    G4double sign = (k % 2) ? -1 : 1;
    d += sign * G4Exp(logSum);
  }
  return d;
}

References G4Exp(), G4Log(), G4Pow::GetInstance(), BoundingBox::kMax, BoundingBox::kMin, G4Pow::logfactorial(), and G4INCL::Math::sign().

---

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

• source/processes/hadronic/util/include/G4Clebsch.hh
• source/processes/hadronic/util/src/G4Clebsch.cc