Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PolarizationTransition.cc
Go to the documentation of this file.
1 // ********************************************************************
2 // * License and Disclaimer *
3 // * *
4 // * The Geant4 software is copyright of the Copyright Holders of *
5 // * the Geant4 Collaboration. It is provided under the terms and *
6 // * conditions of the Geant4 Software License, included in the file *
7 // * LICENSE and available at http://cern.ch/geant4/license . These *
8 // * include a list of copyright holders. *
9 // * *
10 // * Neither the authors of this software system, nor their employing *
11 // * institutes,nor the agencies providing financial support for this *
12 // * work make any representation or warranty, express or implied, *
13 // * regarding this software system or assume any liability for its *
14 // * use. Please see the license in the file LICENSE and URL above *
15 // * for the full disclaimer and the limitation of liability. *
16 // * *
17 // * This code implementation is the result of the scientific and *
18 // * technical work of the GEANT4 collaboration. *
19 // * By using, copying, modifying or distributing the software (or *
20 // * any work based on the software) you agree to acknowledge its *
21 // * use in resulting scientific publications, and indicate your *
22 // * acceptance of all terms of the Geant4 Software license. *
23 // ********************************************************************
24 //
25 //
26 // -------------------------------------------------------------------
27 // GEANT4 Class file
28 //
29 // File name: G4PolarizationTransition
30 //
31 // Author: Jason Detwiler (jasondet@gmail.com)
32 //
33 // Creation date: Aug 2012
34 // -------------------------------------------------------------------
35 #include <iostream>
36 #include <vector>
37 #include "Randomize.hh"
38 
40 #include "G4Clebsch.hh"
41 #include "G4PolynomialPDF.hh"
42 #include "G4Fragment.hh"
43 #include "G4NuclearPolarization.hh"
44 #include "G4SystemOfUnits.hh"
45 #include "G4Exp.hh"
46 
47 using namespace std;
48 
50  : fVerbose(0), fTwoJ1(0), fTwoJ2(0), fLbar(1), fL(0), fDelta(0), kEps(1.e-15),
51  kPolyPDF(0, nullptr, -1, 1)
52 {}
53 
55 {}
56 
58  G4int twoJ2, G4int twoJ1) const
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 }
67 
69  G4int LL, G4int Lprime,
70  G4int twoJ2, G4int twoJ1) const
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 }
80 
82  G4int Lbar, G4double delta,
83  G4int Lprime)
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 }
95 
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 }
104 
106  G4int K1) const
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 }
114 
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 }
140 
142  const POLAR& pol)
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 }
212 
214  G4double phi,
215  G4Fragment* frag)
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 }
318 
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 }
337 
G4double FCoefficient(G4int K, G4int L, G4int Lprime, G4int twoJ2, G4int twoJ1) const
void SetCoefficients(const std::vector< G4double > &v)
G4double GetRandomX()
G4double GenerateGammaPhi(G4double cosTheta, const POLAR &)
G4double GammaTransF3Coefficient(G4int K, G4int K2, G4int K1) const
static size_t GetNCoefficients(size_t order)
G4double GenerateGammaCosTheta(const POLAR &)
int G4int
Definition: G4Types.hh:78
std::complex< G4double > G4complex
Definition: G4Types.hh:81
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
std::vector< std::vector< G4complex > > & GetPolarization()
bool G4bool
Definition: G4Types.hh: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
G4double GammaTransFCoefficient(G4int K) const
void SetGammaTransitionData(G4int twoJ1, G4int twoJ2, G4int Lbar, G4double delta=0, G4int Lprime=1)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double F3Coefficient(G4int K, G4int K2, G4int K1, G4int L, G4int Lprime, G4int twoJ2, G4int twoJ1) const
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 UpdatePolarizationToFinalState(G4double cosTheta, G4double phi, G4Fragment *)
static const G4int LL[nN]
void DumpTransitionData(const POLAR &pol) const
#define G4endl
Definition: G4ios.hh:61
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
static constexpr double twopi
Definition: SystemOfUnits.h:55
G4NuclearPolarization * GetNuclearPolarization() const
Definition: G4Fragment.hh:458
G4double GetCoefficient(size_t i, size_t order)
void SetPolarization(std::vector< std::vector< G4complex > > &p)