G4PolarizationTransition Class Reference

#include <G4PolarizationTransition.hh>

Public Member Functions
	G4PolarizationTransition ()
	~G4PolarizationTransition ()
G4double	FCoefficient (G4int K, G4int L, G4int Lprime, G4int twoJ2, G4int twoJ1) const
G4double	F3Coefficient (G4int K, G4int K2, G4int K1, G4int L, G4int Lprime, G4int twoJ2, G4int twoJ1) const
G4double	GammaTransFCoefficient (G4int K) const
G4double	GammaTransF3Coefficient (G4int K, G4int K2, G4int K1) const
void	SetGammaTransitionData (G4int twoJ1, G4int twoJ2, G4int Lbar, G4double delta=0, G4int Lprime=1)
G4double	GenerateGammaCosTheta (const POLAR &)
G4double	GenerateGammaPhi (G4double cosTheta, const POLAR &)
void	UpdatePolarizationToFinalState (G4double cosTheta, G4double phi, G4Fragment *)
void	DumpTransitionData (const POLAR &pol) const
void	SetVerbose (G4int val)

Detailed Description

Definition at line 58 of file G4PolarizationTransition.hh.

Constructor & Destructor Documentation

G4PolarizationTransition::G4PolarizationTransition ( )

explicit

Definition at line 49 of file G4PolarizationTransition.cc.

  : fVerbose(0), fTwoJ1(0), fTwoJ2(0), fLbar(1), fL(0), fDelta(0), kEps(1.e-15), 
    kPolyPDF(0, nullptr, -1, 1)
{}

G4PolarizationTransition::~G4PolarizationTransition

(

)

Definition at line 54 of file G4PolarizationTransition.cc.

{}

Member Function Documentation

void G4PolarizationTransition::DumpTransitionData

(

const POLAR &

pol

)

const

Definition at line 319 of file G4PolarizationTransition.cc.

{
  G4cout << "G4PolarizationTransition: ";
  (fTwoJ1 % 2) ? G4cout << fTwoJ1 << "/2" : G4cout << fTwoJ1/2;
  G4cout << " --(" << fLbar;
  if(fDelta != 0) G4cout << " + " << fDelta << "*" << fL;
  G4cout << ")--> ";
  (fTwoJ2 % 2) ? G4cout << fTwoJ2 << "/2" : G4cout << fTwoJ2/2;
  G4cout << ", P = [ { ";
  for(size_t k=0; k<pol.size(); ++k) {
    if(k>0) G4cout << " }, { ";
    for(size_t kappa=0; kappa<(pol[k]).size(); ++kappa) {
      if(kappa > 0) G4cout << ", ";
      G4cout << (pol[k])[kappa].real() << " + " << (pol[k])[kappa].imag() << "*i";
    }
  }
  G4cout << " } ]" << G4endl;
}

Here is the caller graph for this function:

G4double G4PolarizationTransition::F3Coefficient	(	G4int	K,
		G4int	K2,
		G4int	K1,
		G4int	L,
		G4int	Lprime,
		G4int	twoJ2,
		G4int	twoJ1
	)		const

Definition at line 68 of file G4PolarizationTransition.cc.

{
  G4double fCoeff = G4Clebsch::Wigner3J(2*LL, 2, 2*Lprime, -2, 2*K, 0);
  if(fCoeff == 0) return 0;
  fCoeff *= G4Clebsch::Wigner9J(twoJ2, 2*LL, twoJ1, twoJ2, 2*Lprime, twoJ1, 
                2*K2, 2*K, 2*K1);
  if(fCoeff == 0) return 0;
  if((Lprime+K2+K1+1) % 2) fCoeff = -fCoeff;
  return fCoeff*sqrt((twoJ1+1)*(twoJ2+1)*(2*LL+1)*(2*Lprime+1)*(2*K+1)*(2*K1+1)*(2*K2+1));
}

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4PolarizationTransition::FCoefficient	(	G4int	K,
		G4int	L,
		G4int	Lprime,
		G4int	twoJ2,
		G4int	twoJ1
	)		const

Definition at line 57 of file G4PolarizationTransition.cc.

{
  G4double fCoeff = G4Clebsch::Wigner3J(2*LL, 2, 2*Lprime, -2, 2*K, 0);
  if(fCoeff == 0) return 0;
  fCoeff *= G4Clebsch::Wigner6J(2*LL, 2*Lprime, 2*K, twoJ1, twoJ1, twoJ2);
  if(fCoeff == 0) return 0;
  if(((twoJ1+twoJ2)/2 - 1) %2) fCoeff = -fCoeff;
  return fCoeff*sqrt((2*K+1)*(twoJ1+1)*(2*LL+1)*(2*Lprime+1));
}

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4PolarizationTransition::GammaTransF3Coefficient	(	G4int	K,
		G4int	K2,
		G4int	K1
	)		const

Definition at line 105 of file G4PolarizationTransition.cc.

{
  double transF3Coeff = F3Coefficient(K, K2, K1, fLbar, fLbar, fTwoJ2, fTwoJ1);
  if(fDelta == 0) return transF3Coeff;
  transF3Coeff += 2.*fDelta*F3Coefficient(K, K2, K1, fLbar, fL, fTwoJ2, fTwoJ1);
  transF3Coeff += fDelta*fDelta*F3Coefficient(K, K2, K1, fL, fL, fTwoJ2, fTwoJ1);
  return transF3Coeff;
}

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4PolarizationTransition::GammaTransFCoefficient

(

G4int

K

)

const

Definition at line 96 of file G4PolarizationTransition.cc.

{
  double transFCoeff = FCoefficient(K, fLbar, fLbar, fTwoJ2, fTwoJ1);
  if(fDelta == 0) return transFCoeff;
  transFCoeff += 2.*fDelta*FCoefficient(K, fLbar, fL, fTwoJ2, fTwoJ1);
  transFCoeff += fDelta*fDelta*FCoefficient(K, fL, fL, fTwoJ2, fTwoJ1);
  return transFCoeff;
}

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4PolarizationTransition::GenerateGammaCosTheta

(

const POLAR &

pol

)

Definition at line 115 of file G4PolarizationTransition.cc.

{
  size_t length = pol.size();
  // Isotropic case
  if(length == 1) return G4UniformRand()*2.-1.;

  // kappa > 0 terms integrate out to zero over phi: 0->2pi, so just need (k,0)
  // terms to generate cos theta distribution
  vector<G4double> polyPDFCoeffs(length, 0.0);
  for(size_t k = 0; k < length; k += 2) {
    if(std::abs(((pol)[k])[0].imag()) > kEps && fVerbose > 0) {
      G4cout << "G4PolarizationTransition::GenerateGammaCosTheta WARNING: \n"
         << "          fPolarization[" 
         << k << "][0] has imag component: = " 
         << ((pol)[k])[0].real() << " + " 
         << ((pol)[k])[0].imag() << "*i" << G4endl;
    }
    G4double a_k = sqrt(2*k+1)*GammaTransFCoefficient(k)*((pol)[k])[0].real();
    for(size_t iCoeff=0; iCoeff < fgLegendrePolys.GetNCoefficients(k); ++iCoeff) {
      polyPDFCoeffs[iCoeff] += a_k*fgLegendrePolys.GetCoefficient(iCoeff, k);
    }
  }
  kPolyPDF.SetCoefficients(polyPDFCoeffs);
  return kPolyPDF.GetRandomX();
}

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4PolarizationTransition::GenerateGammaPhi	(	G4double	cosTheta,
		const POLAR &	pol
	)

Definition at line 141 of file G4PolarizationTransition.cc.

{
  size_t length = pol.size();
  // Isotropic case
  G4bool phiIsIsotropic = true;
  for(size_t i=0; i<length; ++i) {
    if(((pol)[i]).size() > 1) {
      phiIsIsotropic = false;
      break;
    }
  }
  if(phiIsIsotropic) { return G4UniformRand()*CLHEP::twopi; }

  // Otherwise, P(phi) can be written as a sum of cos(kappa phi + phi_kappa).
  // Calculate the amplitude and phase for each term
  std::vector<G4double> amp(length, 0.0);
  std::vector<G4double> phase(length, 0.0);
  for(size_t kappa = 0; kappa < length; ++kappa) {
    G4complex cAmpSum(0.,0.);
    for(size_t k = kappa + (kappa % 2); k < length; k += 2) {
      if(kappa >= length || std::abs(((pol)[k])[kappa]) < kEps) { continue; }
      G4double tmpAmp = GammaTransFCoefficient(k);
      if(tmpAmp == 0) { continue; }
      tmpAmp *= sqrt(2*k+1) * fgLegendrePolys.EvalAssocLegendrePoly(k, kappa, cosTheta);
      if(kappa > 0) tmpAmp *= 2.*G4Exp(0.5*(LnFactorial(k-kappa) - LnFactorial(k+kappa)));
      cAmpSum += ((pol)[k])[kappa]*tmpAmp;
    }
    if(kappa == 0 && std::abs(cAmpSum.imag()) > kEps && fVerbose > 0) {
      G4cout << "G4PolarizationTransition::GenerateGammaPhi: WARNING: \n"
         << "    Got complex amp for kappa = 0! A = " << cAmpSum.real() 
         << " + " << cAmpSum.imag() << "*i" << G4endl;
    }
    amp[kappa] = std::abs(cAmpSum);
    phase[kappa] = arg(cAmpSum);
  }

  // Normalize PDF and calc max (note: it's not the true max, but the max
  // assuming that all of the phases line up at a max)
  G4double pdfMax = 0.;
  for(size_t kappa = 0; kappa < amp.size(); ++kappa) { pdfMax += amp[kappa]; }
  if(pdfMax < kEps && fVerbose > 0) {
    G4cout << "G4PolarizationTransition::GenerateGammaPhi: WARNING "
       << "got pdfMax = 0 for \n";
    DumpTransitionData(pol);
    G4cout << "I suspect a non-allowed transition! Returning isotropic phi..." 
       << G4endl;
    return G4UniformRand()*CLHEP::twopi;
  }

  // Finally, throw phi until it falls in the pdf
  for(size_t i=0; i<1000; ++i) {
    G4double phi = G4UniformRand()*CLHEP::twopi;
    G4double prob = G4UniformRand()*pdfMax;
    G4double pdfSum = amp[0];
    for(size_t kappa = 1; kappa < amp.size(); ++kappa) {
      pdfSum += amp[kappa]*cos(phi*kappa + phase[kappa]);
    }
    if(prob < pdfSum) return phi;
    if(pdfSum > pdfMax && fVerbose > 0) {
      G4cout << "G4PolarizationTransition::GenerateGammaPhi: WARNING: \n"
             << "got pdfSum (" << pdfSum << ") > pdfMax (" 
         << pdfMax << ") at phi = " << phi << G4endl;
    }
  }
  if(fVerbose > 0) {
    G4cout << "G4PolarizationTransition::GenerateGammaPhi: WARNING: \n"
       << "no phi generated in 1000 throws! Returning isotropic phi..." << G4endl;
  }
  return G4UniformRand()*CLHEP::twopi;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void G4PolarizationTransition::SetGammaTransitionData	(	G4int	twoJ1,
		G4int	twoJ2,
		G4int	Lbar,
		G4double	delta = 0,
		G4int	Lprime = 1
	)

Definition at line 81 of file G4PolarizationTransition.cc.

{
  fTwoJ1 = std::abs(twoJ1);  // add abs to remove negative J
  fTwoJ2 = std::abs(twoJ2);
  fLbar = Lbar;
  fDelta = delta;
  fL = Lprime;
  if(fVerbose > 1) {
    G4cout << "SET G4PolarizationTransition: J1= " << fTwoJ1 << " J2= " << fTwoJ2
       << " Lbar= " << fLbar << " delta= " << fDelta << " Lp= " << fL << G4endl;
  }
}

Here is the caller graph for this function:

void G4PolarizationTransition::SetVerbose ( G4int  val )

inline

Definition at line 87 of file G4PolarizationTransition.hh.

{ fVerbose = val; };

Here is the caller graph for this function:

void G4PolarizationTransition::UpdatePolarizationToFinalState	(	G4double	cosTheta,
		G4double	phi,
		G4Fragment *	frag
	)

Definition at line 213 of file G4PolarizationTransition.cc.

{
  G4NuclearPolarization* nucpol = frag->GetNuclearPolarization();
  if(nucpol == nullptr) {
    if(fVerbose > 0) {
      G4cout << "G4PolarizationTransition::UpdatePolarizationToFinalState ERROR: "
             << "cannot update NULL nuclear polarization" << G4endl;
    }
    return;
  }

  if(fTwoJ2 == 0) {
    nucpol->Unpolarize();
    return;
  }

  const POLAR& pol = nucpol->GetPolarization();
  size_t newlength = fTwoJ2+1;
  POLAR newPol(newlength);

  for(size_t k2=0; k2<newlength; ++k2) {
    (newPol[k2]).assign(k2+1, 0);
    for(size_t k1=0; k1<pol.size(); ++k1) {
      for(size_t k=0; k<=k1+k2; k+=2) {
        // TransF3Coefficient takes the most time. Only calculate it once per
        // (k, k1, k2) triplet, and wait until the last possible moment to do
        // so. Break out of the inner loops as soon as it is found to be zero.
        G4double tF3 = 0.;
        G4bool recalcTF3 = true;
        for(size_t kappa2=0; kappa2<=k2; ++kappa2) {
          G4int ll = (pol[k1]).size();
          for(G4int kappa1 = 1 - ll; kappa1<ll; ++kappa1) {
            if(k+k2<k1 || k+k1<k2) continue;
            G4complex tmpAmp = (kappa1 < 0) ?  
          conj((pol[k1])[-kappa1])*(kappa1 % 2 ? -1.: 1.) : (pol[k1])[kappa1];
            if(std::abs(tmpAmp) < kEps) continue;
            G4int kappa = kappa1-(G4int)kappa2;
            tmpAmp *= G4Clebsch::Wigner3J(2*k1, -2*kappa1, 2*k, 2*kappa, 2*k2, 2*kappa2);
            if(std::abs(tmpAmp) < kEps) continue;
            if(recalcTF3) {
              tF3 = GammaTransF3Coefficient(k, k2, k1);
              recalcTF3 = false;
            }
            if(std::abs(tF3) < kEps) break;
            tmpAmp *= tF3;
            if(std::abs(tmpAmp) < kEps) continue;
            tmpAmp *= ((kappa1+(G4int)k1)%2 ? -1. : 1.) 
          * sqrt((2.*k+1.)*(2.*k1+1.)/(2.*k2+1.));
            tmpAmp *= fgLegendrePolys.EvalAssocLegendrePoly(k, kappa, cosTheta);
            if(kappa != 0) {
              tmpAmp *= G4Exp(0.5*(LnFactorial(((G4int)k)-kappa) 
                   - LnFactorial(((G4int)k)+kappa)));
              tmpAmp *= polar(1., phi*kappa);
            }
            (newPol[k2])[kappa2] += tmpAmp;
          }
          if(!recalcTF3 && std::abs(tF3) < kEps) break;
        }
      }
    }
  }

  // sanity checks
  if(0.0 == newPol[0][0] && fVerbose > 1) {
    G4cout << "G4PolarizationTransition::UpdatePolarizationToFinalState WARNING:"
       << " P[0][0] is zero!" << G4endl;
    G4cout << "Old pol is: " << *nucpol << G4endl;
    DumpTransitionData(newPol);
    G4cout << "Unpolarizing..." << G4endl;
    nucpol->Unpolarize();
    return;
  }
  if(std::abs((newPol[0])[0].imag()) > kEps && fVerbose > 1) {
    G4cout << "G4PolarizationTransition::UpdatePolarizationToFinalState WARNING: \n"
       << " P[0][0] has a non-zero imaginary part! Unpolarizing..." << G4endl;
    nucpol->Unpolarize();
    return;
  }

  // Normalize and trim
  size_t lastNonEmptyK2 = 0;
  for(size_t k2=0; k2<newlength; ++k2) {
    G4int lastNonZero = -1;
    for(size_t kappa2=0; kappa2<(newPol[k2]).size(); ++kappa2) {
      if(k2 == 0 && kappa2 == 0) {
        lastNonZero = 0;
        continue;
      }
      if(std::abs((newPol[k2])[kappa2]) > 0.0) {
        lastNonZero = kappa2;
        (newPol[k2])[kappa2] /= (newPol[0])[0];
      }
    }
    while((newPol[k2]).size() != size_t (lastNonZero+1)) (newPol[k2]).pop_back();
    if((newPol[k2]).size() > 0) lastNonEmptyK2 = k2;
  }

  // Remove zero-value entries 
  while(newPol.size() != lastNonEmptyK2+1) { newPol.pop_back(); }
  (newPol[0])[0] = 1.0;

  nucpol->SetPolarization(newPol);
}

Here is the call graph for this function:

Here is the caller graph for this function:

---

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

• source/geant4.10.03.p02/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4PolarizationTransition.hh
• source/geant4.10.03.p02/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4PolarizationTransition.cc