| Geant4 Cross Reference |
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
39 #include "G4PolarizationTransition.hh"
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
49 G4PolarizationTransition::G4PolarizationTransition()
50 : fVerbose(1), fTwoJ1(0), fTwoJ2(0), fLbar(1), fL(0), fDelta(0),
51 kEps(1.e-15), kPolyPDF(0, nullptr, -1, 1)
52 {}
53
54 G4double G4PolarizationTransition::FCoefficient(G4int K, G4int LL, G4int Lprime,
55 G4int twoJ2, G4int twoJ1) const
56 {
57 G4double fCoeff = G4Clebsch::Wigner3J(2*LL, 2, 2*Lprime, -2, 2*K, 0);
58 if(fCoeff == 0) return 0;
59 fCoeff *= G4Clebsch::Wigner6J(2*LL, 2*Lprime, 2*K, twoJ1, twoJ1, twoJ2);
60 if(fCoeff == 0) return 0;
61 if(((twoJ1+twoJ2)/2 - 1) %2) fCoeff = -fCoeff;
62 return fCoeff*std::sqrt(G4double((2*K+1)*(twoJ1+1)*(2*LL+1)*(2*Lprime+1)));
63 }
64
65 G4double G4PolarizationTransition::F3Coefficient(G4int K, G4int K2, G4int K1,
66 G4int LL, G4int Lprime,
67 G4int twoJ2, G4int twoJ1) const
68 {
69 G4double fCoeff = G4Clebsch::Wigner3J(2*LL, 2, 2*Lprime, -2, 2*K, 0);
70 if(fCoeff == 0) return 0;
71 fCoeff *= G4Clebsch::Wigner9J(twoJ2, 2*LL, twoJ1, twoJ2, 2*Lprime, twoJ1,
72 2*K2, 2*K, 2*K1);
73 if(fCoeff == 0) return 0;
74 if((Lprime+K2+K1+1) % 2) fCoeff = -fCoeff;
75
76 //AR-13Jun2017 : apply Jason Detwiler's conversion to double
77 // in the argument of sqrt() to avoid integer overflows.
78 return fCoeff*std::sqrt(G4double((twoJ1+1)*(twoJ2+1)*(2*LL+1))
79 *G4double((2*Lprime+1)*(2*K+1)*(2*K1+1)*(2*K2+1)));
80 }
81
82 G4double G4PolarizationTransition::GammaTransFCoefficient(G4int K) const
83 {
84 double transFCoeff = FCoefficient(K, fLbar, fLbar, fTwoJ2, fTwoJ1);
85 if(fDelta == 0) return transFCoeff;
86 transFCoeff += 2.*fDelta*FCoefficient(K, fLbar, fL, fTwoJ2, fTwoJ1);
87 transFCoeff += fDelta*fDelta*FCoefficient(K, fL, fL, fTwoJ2, fTwoJ1);
88 return transFCoeff;
89 }
90
91 G4double G4PolarizationTransition::GammaTransF3Coefficient(G4int K, G4int K2,
92 G4int K1) const
93 {
94 double transF3Coeff = F3Coefficient(K, K2, K1, fLbar, fLbar, fTwoJ2, fTwoJ1);
95 if(fDelta == 0) return transF3Coeff;
96 transF3Coeff += 2.*fDelta*F3Coefficient(K, K2, K1, fLbar, fL, fTwoJ2, fTwoJ1);
97 transF3Coeff += fDelta*fDelta*F3Coefficient(K, K2, K1, fL, fL, fTwoJ2, fTwoJ1);
98 return transF3Coeff;
99 }
100
101 G4double G4PolarizationTransition::GenerateGammaCosTheta(const POLAR& pol)
102 {
103 std::size_t length = pol.size();
104 // Isotropic case
105 if(length <= 1) return G4UniformRand()*2.-1.;
106
107 // kappa > 0 terms integrate out to zero over phi: 0->2pi, so just need (k,0)
108 // terms to generate cos theta distribution
109 std::vector<G4double> polyPDFCoeffs(length, 0.0);
110 for(G4int k = 0; k < (G4int)length; k += 2) {
111 if ((pol[k]).size() > 0 ) {
112 if(fVerbose > 1 && std::abs(((pol)[k])[0].imag()) > kEps) {
113 G4cout << "G4PolarizationTransition::GenerateGammaCosTheta WARNING: \n"
114 << " fPolarization["
115 << k << "][0] has imag component: = "
116 << ((pol)[k])[0].real() << " + "
117 << ((pol)[k])[0].imag() << "*i" << G4endl;
118 }
119 G4double a_k = std::sqrt((G4double)(2*k+1))*GammaTransFCoefficient(k)*((pol)[k])[0].real();
120 std::size_t nCoeff = fgLegendrePolys.GetNCoefficients(k);
121 for(std::size_t iCoeff=0; iCoeff < nCoeff; ++iCoeff) {
122 polyPDFCoeffs[iCoeff] += a_k*fgLegendrePolys.GetCoefficient(iCoeff, k);
123 }
124 } else {
125 G4cout << "G4PolarizationTransition::GenerateGammaCosTheta: WARNING: \n"
126 << " size of pol[" << k << "] = " << (pol[k]).size()
127 << " returning isotropic " << G4endl;
128 return G4UniformRand()*2.-1.;
129 }
130 }
131 if(fVerbose > 1 && polyPDFCoeffs[polyPDFCoeffs.size()-1] == 0) {
132 G4cout << "G4PolarizationTransition::GenerateGammaCosTheta: WARNING: "
133 << "got zero highest-order coefficient." << G4endl;
134 DumpTransitionData(pol);
135 }
136 kPolyPDF.SetCoefficients(polyPDFCoeffs);
137 return kPolyPDF.GetRandomX();
138 }
139
140 G4double G4PolarizationTransition::GenerateGammaPhi(G4double& cosTheta,
141 const POLAR& pol)
142 {
143 G4int length = (G4int)pol.size();
144 // Isotropic case
145 G4bool phiIsIsotropic = true;
146 for(G4int i=0; i<length; ++i) {
147 if(((pol)[i]).size() > 1) {
148 phiIsIsotropic = false;
149 break;
150 }
151 }
152 if(phiIsIsotropic) { return G4UniformRand()*CLHEP::twopi; }
153 // Otherwise, P(phi) can be written as a sum of cos(kappa phi + phi_kappa).
154 // Calculate the amplitude and phase for each term
155 std::vector<G4double> amp(length, 0.0);
156 std::vector<G4double> phase(length, 0.0);
157 for(G4int kappa = 0; kappa < length; ++kappa) {
158 G4complex cAmpSum(0.,0.);
159 for(G4int k = kappa + (kappa % 2); k < length; k += 2) {
160 G4int kmax = (G4int)pol[k].size();
161 if(kmax > 0) {
162 if(kappa >= kmax || std::abs(((pol)[k])[kappa]) < kEps) { continue; }
163 G4double tmpAmp = GammaTransFCoefficient(k);
164 if(tmpAmp == 0) { continue; }
165 tmpAmp *= std::sqrt((G4double)(2*k+1))
166 * fgLegendrePolys.EvalAssocLegendrePoly(k, kappa, cosTheta);
167 if(kappa > 0) tmpAmp *= 2.*G4Exp(0.5*(LnFactorial(k-kappa) - LnFactorial(k+kappa)));
168 cAmpSum += ((pol)[k])[kappa]*tmpAmp;
169 } else {
170 if(fVerbose > 1) {
171 G4cout << "G4PolarizationTransition::GenerateGammaPhi: WARNING: \n"
172 << " size of pol[" << k << "] = " << (pol[k]).size()
173 << " returning isotropic " << G4endl;
174 }
175 return G4UniformRand()*CLHEP::twopi;
176 }
177 }
178 if(fVerbose > 1 && kappa == 0 && std::abs(cAmpSum.imag()) > kEps) {
179 G4cout << "G4PolarizationTransition::GenerateGammaPhi: WARNING: \n"
180 << " Got complex amp for kappa = 0! A = " << cAmpSum.real()
181 << " + " << cAmpSum.imag() << "*i" << G4endl;
182 }
183 amp[kappa] = std::abs(cAmpSum);
184 phase[kappa] = arg(cAmpSum);
185 }
186
187 // Normalize PDF and calc max (note: it's not the true max, but the max
188 // assuming that all of the phases line up at a max)
189 G4double pdfMax = 0.;
190 for(G4int kappa = 0; kappa < length; ++kappa) { pdfMax += amp[kappa]; }
191 if(fVerbose > 1 && pdfMax < kEps) {
192 G4cout << "G4PolarizationTransition::GenerateGammaPhi: WARNING "
193 << "got pdfMax = 0 for \n";
194 DumpTransitionData(pol);
195 G4cout << "I suspect a non-allowed transition! Returning isotropic phi..."
196 << G4endl;
197 return G4UniformRand()*CLHEP::twopi;
198 }
199
200 // Finally, throw phi until it falls in the pdf
201 for(std::size_t i=0; i<100; ++i) {
202 G4double phi = G4UniformRand()*CLHEP::twopi;
203 G4double prob = G4UniformRand()*pdfMax;
204 G4double pdfSum = amp[0];
205 for(G4int kappa = 1; kappa < length; ++kappa) {
206 pdfSum += amp[kappa]*std::cos(phi*kappa + phase[kappa]);
207 }
208 if(fVerbose > 1 && pdfSum > pdfMax) {
209 G4cout << "G4PolarizationTransition::GenerateGammaPhi: WARNING: \n"
210 << "got pdfSum (" << pdfSum << ") > pdfMax ("
211 << pdfMax << ") at phi = " << phi << G4endl;
212 }
213 if(prob <= pdfSum) { return phi; }
214 }
215 if(fVerbose > 1) {
216 G4cout << "G4PolarizationTransition::GenerateGammaPhi: WARNING: \n"
217 << "no phi generated in 1000 throws! Returning isotropic phi..."
218 << G4endl;
219 }
220 return G4UniformRand()*CLHEP::twopi;
221 }
222
223 void G4PolarizationTransition::SampleGammaTransition(
224 G4NuclearPolarization* nucpol,
225 G4int twoJ1, G4int twoJ2,
226 G4int L0, G4int Lp, G4double mpRatio,
227 G4double& cosTheta, G4double& phi)
228 {
229 if(nucpol == nullptr) {
230 if(fVerbose > 1) {
231 G4cout << "G4PolarizationTransition::SampleGammaTransition ERROR: "
232 << "cannot update NULL nuclear polarization" << G4endl;
233 }
234 return;
235 }
236 fTwoJ1 = twoJ1; // add abs to remove negative J
237 fTwoJ2 = twoJ2;
238 fLbar = L0;
239 fL = Lp;
240 fDelta = mpRatio;
241 if(fVerbose > 2) {
242 G4cout << "G4PolarizationTransition: 2J1= " << fTwoJ1 << " 2J2= " << fTwoJ2
243 << " Lbar= " << fLbar << " delta= " << fDelta << " Lp= " << fL
244 << G4endl;
245 G4cout << *nucpol << G4endl;
246 }
247
248 const POLAR& pol = nucpol->GetPolarization();
249
250 if(fTwoJ1 == 0) {
251 nucpol->Unpolarize();
252 cosTheta = 2*G4UniformRand() - 1.0;
253 phi = CLHEP::twopi*G4UniformRand();
254 }
255 else {
256 cosTheta = GenerateGammaCosTheta(pol);
257 phi = GenerateGammaPhi(cosTheta, pol);
258 if (fVerbose > 2) {
259 G4cout << "G4PolarizationTransition::SampleGammaTransition: cosTheta= "
260 << cosTheta << " phi= " << phi << G4endl;
261 }
262 }
263
264 if(fTwoJ2 == 0) {
265 nucpol->Unpolarize();
266 return;
267 }
268
269 std::size_t newlength = fTwoJ2+1;
270 //POLAR newPol(newlength);
271 POLAR newPol;
272
273 for(G4int k2=0; k2<(G4int)newlength; ++k2) {
274 std::vector<G4complex> npolar;
275 npolar.resize(k2+1, 0);
276 //(newPol[k2]).assign(k2+1, 0);
277 for(G4int k1=0; k1<(G4int)pol.size(); ++k1) {
278 for(G4int k=0; k<=k1+k2; k+=2) {
279 // TransF3Coefficient takes the most time. Only calculate it once per
280 // (k, k1, k2) triplet, and wait until the last possible moment to do
281 // so. Break out of the inner loops as soon as it is found to be zero.
282 G4double tF3 = 0.;
283 G4bool recalcTF3 = true;
284 for(G4int kappa2=0; kappa2<=k2; ++kappa2) {
285 G4int ll = (G4int)pol[k1].size();
286 for(G4int kappa1 = 1 - ll; kappa1<ll; ++kappa1) {
287 if(k+k2<k1 || k+k1<k2) continue;
288 G4complex tmpAmp = (kappa1 < 0) ?
289 conj((pol[k1])[-kappa1])*(kappa1 % 2 ? -1.: 1.) : (pol[k1])[kappa1];
290 if(std::abs(tmpAmp) < kEps) continue;
291 G4int kappa = kappa1-(G4int)kappa2;
292 tmpAmp *= G4Clebsch::Wigner3J(2*k1, -2*kappa1, 2*k, 2*kappa, 2*k2, 2*kappa2);
293 if(std::abs(tmpAmp) < kEps) continue;
294 if(recalcTF3) {
295 tF3 = GammaTransF3Coefficient(k, k2, k1);
296 recalcTF3 = false;
297 }
298 if(std::abs(tF3) < kEps) break;
299 tmpAmp *= tF3;
300 if(std::abs(tmpAmp) < kEps) continue;
301 tmpAmp *= ((kappa1+(G4int)k1)%2 ? -1. : 1.)
302 * sqrt((2.*k+1.)*(2.*k1+1.)/(2.*k2+1.));
303 //AR-13Jun2017 Useful for debugging very long computations
304 //G4cout << "G4PolarizationTransition::UpdatePolarizationToFinalState : k1=" << k1
305 // << " ; k2=" << k2 << " ; kappa1=" << kappa1 << " ; kappa2=" << kappa2 << G4endl;
306 tmpAmp *= fgLegendrePolys.EvalAssocLegendrePoly(k, kappa, cosTheta);
307 if(kappa != 0) {
308 tmpAmp *= G4Exp(0.5*(LnFactorial(((G4int)k)-kappa)
309 - LnFactorial(((G4int)k)+kappa)));
310 tmpAmp *= polar(1., phi*kappa);
311 }
312 npolar[kappa2] += tmpAmp;
313 }
314 if(!recalcTF3 && std::abs(tF3) < kEps) break;
315 }
316 }
317 }
318 newPol.push_back(std::move(npolar));
319 }
320
321 // sanity checks
322 if(fVerbose > 2 && 0.0 == newPol[0][0]) {
323 G4cout << "G4PolarizationTransition::SampleGammaTransition WARNING:"
324 << " P[0][0] is zero!" << G4endl;
325 G4cout << "Old pol is: " << *nucpol << G4endl;
326 DumpTransitionData(newPol);
327 G4cout << "Unpolarizing..." << G4endl;
328 nucpol->Unpolarize();
329 return;
330 }
331 if(fVerbose > 2 && std::abs((newPol[0])[0].imag()) > kEps) {
332 G4cout << "G4PolarizationTransition::SampleGammaTransition WARNING: \n"
333 << " P[0][0] has a non-zero imaginary part! Unpolarizing..."
334 << G4endl;
335 nucpol->Unpolarize();
336 return;
337 }
338 if(fVerbose > 2) {
339 G4cout << "Before normalization: " << G4endl;
340 DumpTransitionData(newPol);
341 }
342 // Normalize and trim
343 std::size_t lastNonEmptyK2 = 0;
344 for(std::size_t k2=0; k2<newlength; ++k2) {
345 G4int lastNonZero = -1;
346 for(std::size_t kappa2=0; kappa2<(newPol[k2]).size(); ++kappa2) {
347 if(k2 == 0 && kappa2 == 0) {
348 lastNonZero = 0;
349 continue;
350 }
351 if(std::abs((newPol[k2])[kappa2]) > 0.0) {
352 lastNonZero = (G4int)kappa2;
353 (newPol[k2])[kappa2] /= (newPol[0])[0];
354 }
355 }
356 while((newPol[k2]).size() != std::size_t (lastNonZero+1)) (newPol[k2]).pop_back();
357 if((newPol[k2]).size() > 0) lastNonEmptyK2 = k2;
358 }
359
360 // Remove zero-value entries
361 while(newPol.size() != lastNonEmptyK2+1) { newPol.pop_back(); }
362 (newPol[0])[0] = 1.0;
363
364 nucpol->SetPolarization(newPol);
365 if(fVerbose > 2) {
366 G4cout << "Updated polarization: " << *nucpol << G4endl;
367 }
368 }
369
370 void G4PolarizationTransition::DumpTransitionData(const POLAR& pol) const
371 {
372 G4cout << "G4PolarizationTransition: ";
373 (fTwoJ1 % 2) ? G4cout << fTwoJ1 << "/2" : G4cout << fTwoJ1/2;
374 G4cout << " --(" << fLbar;
375 if(fDelta != 0) G4cout << " + " << fDelta << "*" << fL;
376 G4cout << ")--> ";
377 (fTwoJ2 % 2) ? G4cout << fTwoJ2 << "/2" : G4cout << fTwoJ2/2;
378 G4cout << ", P = [ { ";
379 for(std::size_t k=0; k<pol.size(); ++k) {
380 if(k>0) G4cout << " }, { ";
381 for(std::size_t kappa=0; kappa<(pol[k]).size(); ++kappa) {
382 if(kappa > 0) G4cout << ", ";
383 G4cout << (pol[k])[kappa].real() << " + " << (pol[k])[kappa].imag() << "*i";
384 }
385 }
386 G4cout << " } ]" << G4endl;
387 }
388
389