Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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.

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

Definition at line 54 of file G4PolarizationTransition.cc.

55 {}

Member Function Documentation

void G4PolarizationTransition::DumpTransitionData ( const POLAR &  pol) const

Definition at line 319 of file G4PolarizationTransition.cc.

320 {
321  G4cout << "G4PolarizationTransition: ";
322  (fTwoJ1 % 2) ? G4cout << fTwoJ1 << "/2" : G4cout << fTwoJ1/2;
323  G4cout << " --(" << fLbar;
324  if(fDelta != 0) G4cout << " + " << fDelta << "*" << fL;
325  G4cout << ")--> ";
326  (fTwoJ2 % 2) ? G4cout << fTwoJ2 << "/2" : G4cout << fTwoJ2/2;
327  G4cout << ", P = [ { ";
328  for(size_t k=0; k<pol.size(); ++k) {
329  if(k>0) G4cout << " }, { ";
330  for(size_t kappa=0; kappa<(pol[k]).size(); ++kappa) {
331  if(kappa > 0) G4cout << ", ";
332  G4cout << (pol[k])[kappa].real() << " + " << (pol[k])[kappa].imag() << "*i";
333  }
334  }
335  G4cout << " } ]" << G4endl;
336 }
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61

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.

71 {
72  G4double fCoeff = G4Clebsch::Wigner3J(2*LL, 2, 2*Lprime, -2, 2*K, 0);
73  if(fCoeff == 0) return 0;
74  fCoeff *= G4Clebsch::Wigner9J(twoJ2, 2*LL, twoJ1, twoJ2, 2*Lprime, twoJ1,
75  2*K2, 2*K, 2*K1);
76  if(fCoeff == 0) return 0;
77  if((Lprime+K2+K1+1) % 2) fCoeff = -fCoeff;
78  return fCoeff*sqrt((twoJ1+1)*(twoJ2+1)*(2*LL+1)*(2*Lprime+1)*(2*K+1)*(2*K1+1)*(2*K2+1));
79 }
static G4double Wigner9J(G4int twoJ1, G4int twoJ2, G4int twoJ3, G4int twoJ4, G4int twoJ5, G4int twoJ6, G4int twoJ7, G4int twoJ8, G4int twoJ9)
Definition: G4Clebsch.cc:533
static G4double Wigner3J(G4double j1, G4double j2, G4double j3, G4double m1, G4double m2, G4double m3)
Definition: G4Clebsch.cc:345
static const G4int LL[nN]
double G4double
Definition: G4Types.hh:76

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.

59 {
60  G4double fCoeff = G4Clebsch::Wigner3J(2*LL, 2, 2*Lprime, -2, 2*K, 0);
61  if(fCoeff == 0) return 0;
62  fCoeff *= G4Clebsch::Wigner6J(2*LL, 2*Lprime, 2*K, twoJ1, twoJ1, twoJ2);
63  if(fCoeff == 0) return 0;
64  if(((twoJ1+twoJ2)/2 - 1) %2) fCoeff = -fCoeff;
65  return fCoeff*sqrt((2*K+1)*(twoJ1+1)*(2*LL+1)*(2*Lprime+1));
66 }
static G4double Wigner3J(G4double j1, G4double j2, G4double j3, G4double m1, G4double m2, G4double m3)
Definition: G4Clebsch.cc:345
static const G4int LL[nN]
static G4double Wigner6J(G4int twoJ1, G4int twoJ2, G4int twoJ3, G4int twoJ4, G4int twoJ5, G4int twoJ6)
Definition: G4Clebsch.cc:432
double G4double
Definition: G4Types.hh:76

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.

107 {
108  double transF3Coeff = F3Coefficient(K, K2, K1, fLbar, fLbar, fTwoJ2, fTwoJ1);
109  if(fDelta == 0) return transF3Coeff;
110  transF3Coeff += 2.*fDelta*F3Coefficient(K, K2, K1, fLbar, fL, fTwoJ2, fTwoJ1);
111  transF3Coeff += fDelta*fDelta*F3Coefficient(K, K2, K1, fL, fL, fTwoJ2, fTwoJ1);
112  return transF3Coeff;
113 }
G4double F3Coefficient(G4int K, G4int K2, G4int K1, G4int L, G4int Lprime, G4int twoJ2, G4int twoJ1) const

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.

97 {
98  double transFCoeff = FCoefficient(K, fLbar, fLbar, fTwoJ2, fTwoJ1);
99  if(fDelta == 0) return transFCoeff;
100  transFCoeff += 2.*fDelta*FCoefficient(K, fLbar, fL, fTwoJ2, fTwoJ1);
101  transFCoeff += fDelta*fDelta*FCoefficient(K, fL, fL, fTwoJ2, fTwoJ1);
102  return transFCoeff;
103 }
G4double FCoefficient(G4int K, G4int L, G4int Lprime, G4int twoJ2, G4int twoJ1) const

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.

116 {
117  size_t length = pol.size();
118  // Isotropic case
119  if(length == 1) return G4UniformRand()*2.-1.;
120 
121  // kappa > 0 terms integrate out to zero over phi: 0->2pi, so just need (k,0)
122  // terms to generate cos theta distribution
123  vector<G4double> polyPDFCoeffs(length, 0.0);
124  for(size_t k = 0; k < length; k += 2) {
125  if(std::abs(((pol)[k])[0].imag()) > kEps && fVerbose > 0) {
126  G4cout << "G4PolarizationTransition::GenerateGammaCosTheta WARNING: \n"
127  << " fPolarization["
128  << k << "][0] has imag component: = "
129  << ((pol)[k])[0].real() << " + "
130  << ((pol)[k])[0].imag() << "*i" << G4endl;
131  }
132  G4double a_k = sqrt(2*k+1)*GammaTransFCoefficient(k)*((pol)[k])[0].real();
133  for(size_t iCoeff=0; iCoeff < fgLegendrePolys.GetNCoefficients(k); ++iCoeff) {
134  polyPDFCoeffs[iCoeff] += a_k*fgLegendrePolys.GetCoefficient(iCoeff, k);
135  }
136  }
137  kPolyPDF.SetCoefficients(polyPDFCoeffs);
138  return kPolyPDF.GetRandomX();
139 }
void SetCoefficients(const std::vector< G4double > &v)
G4double GetRandomX()
static size_t GetNCoefficients(size_t order)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
G4double GammaTransFCoefficient(G4int K) const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4double GetCoefficient(size_t i, size_t order)

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.

143 {
144  size_t length = pol.size();
145  // Isotropic case
146  G4bool phiIsIsotropic = true;
147  for(size_t i=0; i<length; ++i) {
148  if(((pol)[i]).size() > 1) {
149  phiIsIsotropic = false;
150  break;
151  }
152  }
153  if(phiIsIsotropic) { return G4UniformRand()*CLHEP::twopi; }
154 
155  // Otherwise, P(phi) can be written as a sum of cos(kappa phi + phi_kappa).
156  // Calculate the amplitude and phase for each term
157  std::vector<G4double> amp(length, 0.0);
158  std::vector<G4double> phase(length, 0.0);
159  for(size_t kappa = 0; kappa < length; ++kappa) {
160  G4complex cAmpSum(0.,0.);
161  for(size_t k = kappa + (kappa % 2); k < length; k += 2) {
162  if(kappa >= length || std::abs(((pol)[k])[kappa]) < kEps) { continue; }
163  G4double tmpAmp = GammaTransFCoefficient(k);
164  if(tmpAmp == 0) { continue; }
165  tmpAmp *= sqrt(2*k+1) * fgLegendrePolys.EvalAssocLegendrePoly(k, kappa, cosTheta);
166  if(kappa > 0) tmpAmp *= 2.*G4Exp(0.5*(LnFactorial(k-kappa) - LnFactorial(k+kappa)));
167  cAmpSum += ((pol)[k])[kappa]*tmpAmp;
168  }
169  if(kappa == 0 && std::abs(cAmpSum.imag()) > kEps && fVerbose > 0) {
170  G4cout << "G4PolarizationTransition::GenerateGammaPhi: WARNING: \n"
171  << " Got complex amp for kappa = 0! A = " << cAmpSum.real()
172  << " + " << cAmpSum.imag() << "*i" << G4endl;
173  }
174  amp[kappa] = std::abs(cAmpSum);
175  phase[kappa] = arg(cAmpSum);
176  }
177 
178  // Normalize PDF and calc max (note: it's not the true max, but the max
179  // assuming that all of the phases line up at a max)
180  G4double pdfMax = 0.;
181  for(size_t kappa = 0; kappa < amp.size(); ++kappa) { pdfMax += amp[kappa]; }
182  if(pdfMax < kEps && fVerbose > 0) {
183  G4cout << "G4PolarizationTransition::GenerateGammaPhi: WARNING "
184  << "got pdfMax = 0 for \n";
185  DumpTransitionData(pol);
186  G4cout << "I suspect a non-allowed transition! Returning isotropic phi..."
187  << G4endl;
188  return G4UniformRand()*CLHEP::twopi;
189  }
190 
191  // Finally, throw phi until it falls in the pdf
192  for(size_t i=0; i<1000; ++i) {
194  G4double prob = G4UniformRand()*pdfMax;
195  G4double pdfSum = amp[0];
196  for(size_t kappa = 1; kappa < amp.size(); ++kappa) {
197  pdfSum += amp[kappa]*cos(phi*kappa + phase[kappa]);
198  }
199  if(prob < pdfSum) return phi;
200  if(pdfSum > pdfMax && fVerbose > 0) {
201  G4cout << "G4PolarizationTransition::GenerateGammaPhi: WARNING: \n"
202  << "got pdfSum (" << pdfSum << ") > pdfMax ("
203  << pdfMax << ") at phi = " << phi << G4endl;
204  }
205  }
206  if(fVerbose > 0) {
207  G4cout << "G4PolarizationTransition::GenerateGammaPhi: WARNING: \n"
208  << "no phi generated in 1000 throws! Returning isotropic phi..." << G4endl;
209  }
210  return G4UniformRand()*CLHEP::twopi;
211 }
std::complex< G4double > G4complex
Definition: G4Types.hh:81
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
G4double GammaTransFCoefficient(G4int K) const
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double EvalAssocLegendrePoly(G4int l, G4int m, G4double x)
void DumpTransitionData(const POLAR &pol) const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
static constexpr double twopi
Definition: SystemOfUnits.h:55

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.

84 {
85  fTwoJ1 = std::abs(twoJ1); // add abs to remove negative J
86  fTwoJ2 = std::abs(twoJ2);
87  fLbar = Lbar;
88  fDelta = delta;
89  fL = Lprime;
90  if(fVerbose > 1) {
91  G4cout << "SET G4PolarizationTransition: J1= " << fTwoJ1 << " J2= " << fTwoJ2
92  << " Lbar= " << fLbar << " delta= " << fDelta << " Lp= " << fL << G4endl;
93  }
94 }
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61

Here is the caller graph for this function:

void G4PolarizationTransition::SetVerbose ( G4int  val)
inline

Definition at line 87 of file G4PolarizationTransition.hh.

87 { 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.

216 {
218  if(nucpol == nullptr) {
219  if(fVerbose > 0) {
220  G4cout << "G4PolarizationTransition::UpdatePolarizationToFinalState ERROR: "
221  << "cannot update NULL nuclear polarization" << G4endl;
222  }
223  return;
224  }
225 
226  if(fTwoJ2 == 0) {
227  nucpol->Unpolarize();
228  return;
229  }
230 
231  const POLAR& pol = nucpol->GetPolarization();
232  size_t newlength = fTwoJ2+1;
233  POLAR newPol(newlength);
234 
235  for(size_t k2=0; k2<newlength; ++k2) {
236  (newPol[k2]).assign(k2+1, 0);
237  for(size_t k1=0; k1<pol.size(); ++k1) {
238  for(size_t k=0; k<=k1+k2; k+=2) {
239  // TransF3Coefficient takes the most time. Only calculate it once per
240  // (k, k1, k2) triplet, and wait until the last possible moment to do
241  // so. Break out of the inner loops as soon as it is found to be zero.
242  G4double tF3 = 0.;
243  G4bool recalcTF3 = true;
244  for(size_t kappa2=0; kappa2<=k2; ++kappa2) {
245  G4int ll = (pol[k1]).size();
246  for(G4int kappa1 = 1 - ll; kappa1<ll; ++kappa1) {
247  if(k+k2<k1 || k+k1<k2) continue;
248  G4complex tmpAmp = (kappa1 < 0) ?
249  conj((pol[k1])[-kappa1])*(kappa1 % 2 ? -1.: 1.) : (pol[k1])[kappa1];
250  if(std::abs(tmpAmp) < kEps) continue;
251  G4int kappa = kappa1-(G4int)kappa2;
252  tmpAmp *= G4Clebsch::Wigner3J(2*k1, -2*kappa1, 2*k, 2*kappa, 2*k2, 2*kappa2);
253  if(std::abs(tmpAmp) < kEps) continue;
254  if(recalcTF3) {
255  tF3 = GammaTransF3Coefficient(k, k2, k1);
256  recalcTF3 = false;
257  }
258  if(std::abs(tF3) < kEps) break;
259  tmpAmp *= tF3;
260  if(std::abs(tmpAmp) < kEps) continue;
261  tmpAmp *= ((kappa1+(G4int)k1)%2 ? -1. : 1.)
262  * sqrt((2.*k+1.)*(2.*k1+1.)/(2.*k2+1.));
263  tmpAmp *= fgLegendrePolys.EvalAssocLegendrePoly(k, kappa, cosTheta);
264  if(kappa != 0) {
265  tmpAmp *= G4Exp(0.5*(LnFactorial(((G4int)k)-kappa)
266  - LnFactorial(((G4int)k)+kappa)));
267  tmpAmp *= polar(1., phi*kappa);
268  }
269  (newPol[k2])[kappa2] += tmpAmp;
270  }
271  if(!recalcTF3 && std::abs(tF3) < kEps) break;
272  }
273  }
274  }
275  }
276 
277  // sanity checks
278  if(0.0 == newPol[0][0] && fVerbose > 1) {
279  G4cout << "G4PolarizationTransition::UpdatePolarizationToFinalState WARNING:"
280  << " P[0][0] is zero!" << G4endl;
281  G4cout << "Old pol is: " << *nucpol << G4endl;
282  DumpTransitionData(newPol);
283  G4cout << "Unpolarizing..." << G4endl;
284  nucpol->Unpolarize();
285  return;
286  }
287  if(std::abs((newPol[0])[0].imag()) > kEps && fVerbose > 1) {
288  G4cout << "G4PolarizationTransition::UpdatePolarizationToFinalState WARNING: \n"
289  << " P[0][0] has a non-zero imaginary part! Unpolarizing..." << G4endl;
290  nucpol->Unpolarize();
291  return;
292  }
293 
294  // Normalize and trim
295  size_t lastNonEmptyK2 = 0;
296  for(size_t k2=0; k2<newlength; ++k2) {
297  G4int lastNonZero = -1;
298  for(size_t kappa2=0; kappa2<(newPol[k2]).size(); ++kappa2) {
299  if(k2 == 0 && kappa2 == 0) {
300  lastNonZero = 0;
301  continue;
302  }
303  if(std::abs((newPol[k2])[kappa2]) > 0.0) {
304  lastNonZero = kappa2;
305  (newPol[k2])[kappa2] /= (newPol[0])[0];
306  }
307  }
308  while((newPol[k2]).size() != size_t (lastNonZero+1)) (newPol[k2]).pop_back();
309  if((newPol[k2]).size() > 0) lastNonEmptyK2 = k2;
310  }
311 
312  // Remove zero-value entries
313  while(newPol.size() != lastNonEmptyK2+1) { newPol.pop_back(); }
314  (newPol[0])[0] = 1.0;
315 
316  nucpol->SetPolarization(newPol);
317 }
G4double GammaTransF3Coefficient(G4int K, G4int K2, G4int K1) const
int G4int
Definition: G4Types.hh:78
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4GLOB_DLL std::ostream G4cout
std::vector< std::vector< G4complex > > & GetPolarization()
bool G4bool
Definition: G4Types.hh:79
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static G4double Wigner3J(G4double j1, G4double j2, G4double j3, G4double m1, G4double m2, G4double m3)
Definition: G4Clebsch.cc:345
G4double EvalAssocLegendrePoly(G4int l, G4int m, G4double x)
void DumpTransitionData(const POLAR &pol) const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4NuclearPolarization * GetNuclearPolarization() const
Definition: G4Fragment.hh:458
void SetPolarization(std::vector< std::vector< G4complex > > &p)

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: