G4PolarizationTransition Class Reference

#include <G4PolarizationTransition.hh>

Collaboration diagram for G4PolarizationTransition:

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

Private Types
typedef std::vector< std::vector< G4complex > >	POLAR

Private Member Functions
G4double	LnFactorial (int k) const

Private Attributes
G4int	fTwoJ1
G4int	fTwoJ2
G4int	fLbar
G4int	fL
G4double	fDelta
G4PolynomialPDF	kPolyPDF
G4LegendrePolynomial	fgLegendrePolys

Detailed Description

Definition at line 58 of file G4PolarizationTransition.hh.

Member Typedef Documentation

POLAR

typedef std::vector< std::vector<G4complex> > G4PolarizationTransition::POLAR

private

Definition at line 60 of file G4PolarizationTransition.hh.

Constructor & Destructor Documentation

G4PolarizationTransition()

G4PolarizationTransition::G4PolarizationTransition

(

)

Definition at line 50 of file G4PolarizationTransition.cc.

  : fTwoJ1(0), fTwoJ2(0), fLbar(1), fL(0), fDelta(0), 
    kPolyPDF(0, nullptr, -1, 1)
{}

~G4PolarizationTransition()

G4PolarizationTransition::~G4PolarizationTransition

(

)

Definition at line 55 of file G4PolarizationTransition.cc.

{}

Member Function Documentation

DumpTransitionData()

void G4PolarizationTransition::DumpTransitionData

(

const POLAR &

pol

)

const

Definition at line 312 of file G4PolarizationTransition.cc.

{
  G4cout << "G4PolarizationTransition transition: ";
  (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:

F3Coefficient()

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

Definition at line 69 of file G4PolarizationTransition.cc.

{
  G4double fCoeff = G4Clebsch::Wigner3J(LL, Lprime, K, 1, -1, 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:

FCoefficient()

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

Definition at line 58 of file G4PolarizationTransition.cc.

{
  G4double fCoeff = G4Clebsch::Wigner3J(LL, Lprime, K, 1, -1, 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:

GammaTransF3Coefficient()

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

Definition at line 102 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:

GammaTransFCoefficient()

G4double G4PolarizationTransition::GammaTransFCoefficient

(

G4int

K

)

const

Definition at line 93 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:

GenerateGammaCosTheta()

G4double G4PolarizationTransition::GenerateGammaCosTheta

(

const POLAR &

pol

)

Definition at line 112 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) {
      G4cout << "Warning: 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:

GenerateGammaPhi()

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

Definition at line 136 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
  vector<G4double> amp(length, 0.0);
  vector<G4double> phase(length, 0.0);
  for(size_t kappa = 0; kappa < length; ++kappa) {
    G4complex cAmpSum = 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.*exp(0.5*(LnFactorial(k-kappa) - LnFactorial(k+kappa)));
      cAmpSum += ((pol)[k])[kappa]*tmpAmp;
    }
    if(kappa == 0 && std::abs(cAmpSum.imag()) > kEps) {
      G4cout << "G4PolarizationTransition::GenerateGammaPhi: WARNING "
         << " got complex amp for kappa = 0! A = " << cAmpSum.real() 
         << " + " << cAmpSum.imag() << "*i" << G4endl;
    }
    amp[kappa] = std::abs(cAmpSum);
    if(amp[kappa] < 0) {
      G4cout << "G4PolarizationTransition::GenerateGammaPhi: WARNING "
         << "got negative abs for kappa = " << kappa << G4endl;
    }
    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) {
    G4cout << "G4PolarizationTransition::GenerateGammaPhi: WARNING "
       << "got pdfMax = 0 for ";
    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) {
      G4cout << "G4PolarizationTransition::GenerateGammaPhi: WARNING "
             << "got pdfSum (" << pdfSum << ") > pdfMax (" 
         << pdfMax << ") at phi = " << phi << G4endl;
    }
  }

  G4cout << "G4PolarizationTransition::GenerateGammaPhi: WARNING "
     << "no phi generated in 1000 throws! Returning isotropic phi..." << G4endl;
  return G4UniformRand()*CLHEP::twopi;
}

Here is the call graph for this function:

LnFactorial()

G4double G4PolarizationTransition::LnFactorial ( int  k ) const

inlineprivate

Definition at line 89 of file G4PolarizationTransition.hh.

{ return G4Pow::GetInstance()->logfactorial(k); }

Here is the call graph for this function:

Here is the caller graph for this function:

SetGammaTransitionData()

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

Definition at line 82 of file G4PolarizationTransition.cc.

{
  fTwoJ1 = twoJ1;
  fTwoJ2 = twoJ2;
  fLbar = Lbar;
  fDelta = delta;
  fL = Lprime;
}

UpdatePolarizationToFinalState()

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

Definition at line 211 of file G4PolarizationTransition.cc.

{
  if(fTwoJ2 == 0) {
    frag->SetNuclearPolarization(nullptr);
    return;
  }
  POLAR pol;
  G4NuclearPolarization* nucpol = frag->GetNuclearPolarization();
  if(nucpol) { pol = nucpol->GetPolarization(); }

  size_t length = pol.size();

  size_t newlength = fTwoJ2+1;
  POLAR newPol;
  newPol.resize(newlength);

  for(size_t k2=0; k2<newlength; ++k2) {
    (newPol[k2]).assign(k2+1, 0);
    for(size_t k1=0; k1<length; ++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(k1, k, k2, -kappa1, kappa, kappa2);
            if(std::abs(tmpAmp) < kEps) continue;
            if(recalcTF3) {
              tF3 = GammaTransF3Coefficient(k, k2, k1);
              recalcTF3 = false;
            }
            if(std::abs(tF3) < kEps) break;
            tmpAmp *= GammaTransF3Coefficient(k, k2, k1);
            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]) {
    frag->SetNuclearPolarization(nullptr);
    return;
  }
  if(std::abs((newPol[0])[0].imag()) > kEps) {
    G4cout << "G4PolarizationTransition::UpdatePolarizationToFinalState WWARNING:"
       << " P[0][0] has a non-zero imaginary part! Unpolarizing..." << G4endl;
    frag->SetNuclearPolarization(nullptr);
    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;
  }
  while(newPol.size() != lastNonEmptyK2+1) { newPol.pop_back(); }
  (newPol[0])[0] = 1.0;

  if(!nucpol) { 
    nucpol = new G4NuclearPolarization(); 
    nucpol->SetPolarization(newPol);
    frag->SetNuclearPolarization(nucpol);
  } else {
    nucpol->SetPolarization(newPol);
  }
}

Here is the call graph for this function:

Member Data Documentation

fDelta

G4double G4PolarizationTransition::fDelta

private

Definition at line 93 of file G4PolarizationTransition.hh.

fgLegendrePolys

G4LegendrePolynomial G4PolarizationTransition::fgLegendrePolys

private

Definition at line 95 of file G4PolarizationTransition.hh.

fL

G4int G4PolarizationTransition::fL

private

Definition at line 92 of file G4PolarizationTransition.hh.

fLbar

G4int G4PolarizationTransition::fLbar

private

Definition at line 92 of file G4PolarizationTransition.hh.

fTwoJ1

G4int G4PolarizationTransition::fTwoJ1

private

Definition at line 91 of file G4PolarizationTransition.hh.

fTwoJ2

G4int G4PolarizationTransition::fTwoJ2

private

Definition at line 91 of file G4PolarizationTransition.hh.

kPolyPDF

G4PolynomialPDF G4PolarizationTransition::kPolyPDF

private

Definition at line 94 of file G4PolarizationTransition.hh.

---

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

• G4PolarizationTransition.hh
• G4PolarizationTransition.cc