Geant4  10.02.p02
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 
46 using namespace std;
47 
48 static const G4double kEps = 1.e-15;
49 
51  : fTwoJ1(0), fTwoJ2(0), fLbar(1), fL(0), fDelta(0),
52  kPolyPDF(0, nullptr, -1, 1)
53 {}
54 
56 {}
57 
59  G4int twoJ2, G4int twoJ1) const
60 {
61  G4double fCoeff = G4Clebsch::Wigner3J(LL, Lprime, K, 1, -1, 0);
62  if(fCoeff == 0) return 0;
63  fCoeff *= G4Clebsch::Wigner6J(2*LL, 2*Lprime, 2*K, twoJ1, twoJ1, twoJ2);
64  if(fCoeff == 0) return 0;
65  if(((twoJ1+twoJ2)/2 - 1) %2) fCoeff = -fCoeff;
66  return fCoeff*sqrt((2*K+1)*(twoJ1+1)*(2*LL+1)*(2*Lprime+1));
67 }
68 
70  G4int LL, G4int Lprime,
71  G4int twoJ2, G4int twoJ1) const
72 {
73  G4double fCoeff = G4Clebsch::Wigner3J(LL, Lprime, K, 1, -1, 0);
74  if(fCoeff == 0) return 0;
75  fCoeff *= G4Clebsch::Wigner9J(twoJ2, 2*LL, twoJ1, twoJ2, 2*Lprime, twoJ1,
76  2*K2, 2*K, 2*K1);
77  if(fCoeff == 0) return 0;
78  if((Lprime+K2+K1+1) % 2) fCoeff = -fCoeff;
79  return fCoeff*sqrt((twoJ1+1)*(twoJ2+1)*(2*LL+1)*(2*Lprime+1)*(2*K+1)*(2*K1+1)*(2*K2+1));
80 }
81 
83  G4int Lbar, G4double delta,
84  G4int Lprime)
85 {
86  fTwoJ1 = twoJ1;
87  fTwoJ2 = twoJ2;
88  fLbar = Lbar;
89  fDelta = delta;
90  fL = Lprime;
91 }
92 
94 {
95  double transFCoeff = FCoefficient(K, fLbar, fLbar, fTwoJ2, fTwoJ1);
96  if(fDelta == 0) return transFCoeff;
97  transFCoeff += 2.*fDelta*FCoefficient(K, fLbar, fL, fTwoJ2, fTwoJ1);
98  transFCoeff += fDelta*fDelta*FCoefficient(K, fL, fL, fTwoJ2, fTwoJ1);
99  return transFCoeff;
100 }
101 
103  G4int K1) const
104 {
105  double transF3Coeff = F3Coefficient(K, K2, K1, fLbar, fLbar, fTwoJ2, fTwoJ1);
106  if(fDelta == 0) return transF3Coeff;
107  transF3Coeff += 2.*fDelta*F3Coefficient(K, K2, K1, fLbar, fL, fTwoJ2, fTwoJ1);
108  transF3Coeff += fDelta*fDelta*F3Coefficient(K, K2, K1, fL, fL, fTwoJ2, fTwoJ1);
109  return transF3Coeff;
110 }
111 
113 {
114  size_t length = pol.size();
115  // Isotropic case
116  if(length == 1) return G4UniformRand()*2.-1.;
117 
118  // kappa > 0 terms integrate out to zero over phi: 0->2pi, so just need (k,0)
119  // terms to generate cos theta distribution
120  vector<G4double> polyPDFCoeffs(length, 0.0);
121  for(size_t k = 0; k < length; k += 2) {
122  if(std::abs(((pol)[k])[0].imag()) > kEps) {
123  G4cout << "Warning: fPolarization[" << k << "][0] has imag component: = "
124  << ((pol)[k])[0].real() << " + "
125  << ((pol)[k])[0].imag() << "*i" << G4endl;
126  }
127  G4double a_k = sqrt(2*k+1)*GammaTransFCoefficient(k)*((pol)[k])[0].real();
128  for(size_t iCoeff=0; iCoeff < fgLegendrePolys.GetNCoefficients(k); ++iCoeff) {
129  polyPDFCoeffs[iCoeff] += a_k*fgLegendrePolys.GetCoefficient(iCoeff, k);
130  }
131  }
132  kPolyPDF.SetCoefficients(polyPDFCoeffs);
133  return kPolyPDF.GetRandomX();
134 }
135 
137  const POLAR& pol)
138 {
139  size_t length = pol.size();
140  // Isotropic case
141  G4bool phiIsIsotropic = true;
142  for(size_t i=0; i<length; ++i) {
143  if(((pol)[i]).size() > 1) {
144  phiIsIsotropic = false;
145  break;
146  }
147  }
148  if(phiIsIsotropic) { return G4UniformRand()*CLHEP::twopi; }
149 
150  // Otherwise, P(phi) can be written as a sum of cos(kappa phi + phi_kappa).
151  // Calculate the amplitude and phase for each term
152  vector<G4double> amp(length, 0.0);
153  vector<G4double> phase(length, 0.0);
154  for(size_t kappa = 0; kappa < length; ++kappa) {
155  G4complex cAmpSum = 0.;
156  for(size_t k = kappa + (kappa % 2); k < length; k += 2) {
157  if(kappa >= length || std::abs(((pol)[k])[kappa]) < kEps) { continue; }
158  G4double tmpAmp = GammaTransFCoefficient(k);
159  if(tmpAmp == 0) { continue; }
160  tmpAmp *= sqrt(2*k+1) * fgLegendrePolys.EvalAssocLegendrePoly(k, kappa, cosTheta);
161  if(kappa > 0) tmpAmp *= 2.*exp(0.5*(LnFactorial(k-kappa) - LnFactorial(k+kappa)));
162  cAmpSum += ((pol)[k])[kappa]*tmpAmp;
163  }
164  if(kappa == 0 && std::abs(cAmpSum.imag()) > kEps) {
165  G4cout << "G4PolarizationTransition::GenerateGammaPhi: WARNING "
166  << " got complex amp for kappa = 0! A = " << cAmpSum.real()
167  << " + " << cAmpSum.imag() << "*i" << G4endl;
168  }
169  amp[kappa] = std::abs(cAmpSum);
170  if(amp[kappa] < 0) {
171  G4cout << "G4PolarizationTransition::GenerateGammaPhi: WARNING "
172  << "got negative abs for kappa = " << kappa << G4endl;
173  }
174  phase[kappa] = arg(cAmpSum);
175  }
176 
177  // Normalize PDF and calc max (note: it's not the true max, but the max
178  // assuming that all of the phases line up at a max)
179  G4double pdfMax = 0;
180  for(size_t kappa = 0; kappa < amp.size(); ++kappa) { pdfMax += amp[kappa]; }
181  if(pdfMax < kEps) {
182  G4cout << "G4PolarizationTransition::GenerateGammaPhi: WARNING "
183  << "got pdfMax = 0 for ";
184  DumpTransitionData(pol);
185  G4cout << "I suspect a non-allowed transition! Returning isotropic phi..."
186  << G4endl;
187  return G4UniformRand()*CLHEP::twopi;
188  }
189 
190  // Finally, throw phi until it falls in the pdf
191  for(size_t i=0; i<1000; ++i) {
193  G4double prob = G4UniformRand()*pdfMax;
194  G4double pdfSum = amp[0];
195  for(size_t kappa = 1; kappa < amp.size(); ++kappa) {
196  pdfSum += amp[kappa]*cos(phi*kappa + phase[kappa]);
197  }
198  if(prob < pdfSum) return phi;
199  if(pdfSum > pdfMax) {
200  G4cout << "G4PolarizationTransition::GenerateGammaPhi: WARNING "
201  << "got pdfSum (" << pdfSum << ") > pdfMax ("
202  << pdfMax << ") at phi = " << phi << G4endl;
203  }
204  }
205 
206  G4cout << "G4PolarizationTransition::GenerateGammaPhi: WARNING "
207  << "no phi generated in 1000 throws! Returning isotropic phi..." << G4endl;
208  return G4UniformRand()*CLHEP::twopi;
209 }
210 
212  G4double phi,
213  G4Fragment* frag)
214 {
215  if(fTwoJ2 == 0) {
216  frag->SetNuclearPolarization(nullptr);
217  return;
218  }
219  POLAR pol;
221  if(nucpol) { pol = nucpol->GetPolarization(); }
222 
223  size_t length = pol.size();
224 
225  size_t newlength = fTwoJ2+1;
226  POLAR newPol;
227  newPol.resize(newlength);
228 
229  for(size_t k2=0; k2<newlength; ++k2) {
230  (newPol[k2]).assign(k2+1, 0);
231  for(size_t k1=0; k1<length; ++k1) {
232  for(size_t k=0; k<=k1+k2; k+=2) {
233  // TransF3Coefficient takes the most time. Only calculate it once per
234  // (k, k1, k2) triplet, and wait until the last possible moment to do
235  // so. Break out of the inner loops as soon as it is found to be zero.
236  G4double tF3 = 0.;
237  G4bool recalcTF3 = true;
238  for(size_t kappa2=0; kappa2<=k2; ++kappa2) {
239  G4int ll = (pol[k1]).size();
240  for(G4int kappa1 = 1 - ll; kappa1<ll; ++kappa1) {
241  if(k+k2<k1 || k+k1<k2) continue;
242  G4complex tmpAmp = (kappa1 < 0) ?
243  conj((pol[k1])[-kappa1])*(kappa1 % 2 ? -1.: 1.) : (pol[k1])[kappa1];
244  if(std::abs(tmpAmp) < kEps) continue;
245  G4int kappa = kappa1-(G4int)kappa2;
246  tmpAmp *= G4Clebsch::Wigner3J(k1, k, k2, -kappa1, kappa, kappa2);
247  if(std::abs(tmpAmp) < kEps) continue;
248  if(recalcTF3) {
249  tF3 = GammaTransF3Coefficient(k, k2, k1);
250  recalcTF3 = false;
251  }
252  if(std::abs(tF3) < kEps) break;
253  tmpAmp *= GammaTransF3Coefficient(k, k2, k1);
254  if(std::abs(tmpAmp) < kEps) continue;
255  tmpAmp *= ((kappa1+(G4int)k1)%2 ? -1. : 1.)
256  * sqrt((2.*k+1.)*(2.*k1+1.)/(2.*k2+1.));
257  tmpAmp *= fgLegendrePolys.EvalAssocLegendrePoly(k, kappa, cosTheta);
258  if(kappa != 0) {
259  tmpAmp *= G4Exp(0.5*(LnFactorial(((G4int)k)-kappa)
260  - LnFactorial(((G4int)k)+kappa)));
261  tmpAmp *= polar(1., phi*kappa);
262  }
263  (newPol[k2])[kappa2] += tmpAmp;
264  }
265  if(!recalcTF3 && std::abs(tF3) < kEps) break;
266  }
267  }
268  }
269  }
270 
271  // sanity checks
272  if(0.0 == newPol[0][0]) {
273  frag->SetNuclearPolarization(nullptr);
274  return;
275  }
276  if(std::abs((newPol[0])[0].imag()) > kEps) {
277  G4cout << "G4PolarizationTransition::UpdatePolarizationToFinalState WWARNING:"
278  << " P[0][0] has a non-zero imaginary part! Unpolarizing..." << G4endl;
279  frag->SetNuclearPolarization(nullptr);
280  return;
281  }
282 
283  // Normalize and trim
284  size_t lastNonEmptyK2 = 0;
285  for(size_t k2=0; k2<newlength; ++k2) {
286  G4int lastNonZero = -1;
287  for(size_t kappa2=0; kappa2<(newPol[k2]).size(); ++kappa2) {
288  if(k2 == 0 && kappa2 == 0) {
289  lastNonZero = 0;
290  continue;
291  }
292  if(std::abs((newPol[k2])[kappa2]) > 0.0) {
293  lastNonZero = kappa2;
294  (newPol[k2])[kappa2] /= (newPol[0])[0];
295  }
296  }
297  while((newPol[k2]).size() != size_t (lastNonZero+1)) (newPol[k2]).pop_back();
298  if((newPol[k2]).size() > 0) lastNonEmptyK2 = k2;
299  }
300  while(newPol.size() != lastNonEmptyK2+1) { newPol.pop_back(); }
301  (newPol[0])[0] = 1.0;
302 
303  if(!nucpol) {
304  nucpol = new G4NuclearPolarization();
305  nucpol->SetPolarization(newPol);
306  frag->SetNuclearPolarization(nucpol);
307  } else {
308  nucpol->SetPolarization(newPol);
309  }
310 }
311 
313 {
314  G4cout << "G4PolarizationTransition transition: ";
315  (fTwoJ1 % 2) ? G4cout << fTwoJ1 << "/2" : G4cout << fTwoJ1/2;
316  G4cout << " --(" << fLbar;
317  if(fDelta > 0) G4cout << " + " << fDelta << "*" << fL;
318  G4cout << ")--> ";
319  (fTwoJ2 % 2) ? G4cout << fTwoJ2 << "/2" : G4cout << fTwoJ2/2;
320  G4cout << ", P = [ { ";
321  for(size_t k=0; k<pol.size(); ++k) {
322  if(k>0) G4cout << " }, { ";
323  for(size_t kappa=0; kappa<(pol[k]).size(); ++kappa) {
324  if(kappa > 0) G4cout << ", ";
325  G4cout << (pol[k])[kappa].real() << " + " << (pol[k])[kappa].imag() << "*i";
326  }
327  }
328  G4cout << " } ]" << G4endl;
329 }
330 
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
G4double LnFactorial(int k) const
static size_t GetNCoefficients(size_t order)
G4double GenerateGammaCosTheta(const POLAR &)
int G4int
Definition: G4Types.hh:78
std::vector< std::vector< G4complex > > POLAR
G4LegendrePolynomial fgLegendrePolys
std::complex< G4double > G4complex
Definition: G4Types.hh:81
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
void SetNuclearPolarization(G4NuclearPolarization *)
Definition: G4Fragment.hh:213
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
static const double twopi
Definition: G4SIunits.hh:75
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]
static const G4double kEps
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
G4NuclearPolarization * GetNuclearPolarization() const
Definition: G4Fragment.hh:430
G4double GetCoefficient(size_t i, size_t order)
void SetPolarization(std::vector< std::vector< G4complex > > &p)