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