Geant4 11.3.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4PolarizationTransition Class Reference

#include <G4PolarizationTransition.hh>

Public Member Functions

 G4PolarizationTransition ()
 
 ~G4PolarizationTransition ()=default
 
void SampleGammaTransition (G4NuclearPolarization *np, G4int twoJ1, G4int twoJ2, G4int L0, G4int Lp, G4double mpRatio, G4double &cosTheta, G4double &phi)
 
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 DumpTransitionData (const POLAR &pol) const
 
void SetVerbose (G4int val)
 

Detailed Description

Definition at line 58 of file G4PolarizationTransition.hh.

Constructor & Destructor Documentation

◆ G4PolarizationTransition()

G4PolarizationTransition::G4PolarizationTransition ( )

Definition at line 49 of file G4PolarizationTransition.cc.

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

◆ ~G4PolarizationTransition()

G4PolarizationTransition::~G4PolarizationTransition ( )
default

Member Function Documentation

◆ DumpTransitionData()

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

Definition at line 370 of file G4PolarizationTransition.cc.

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}
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout

Referenced by SampleGammaTransition().

◆ F3Coefficient()

G4double G4PolarizationTransition::F3Coefficient ( G4int K,
G4int K2,
G4int K1,
G4int L,
G4int Lprime,
G4int twoJ2,
G4int twoJ1 ) const

Definition at line 65 of file G4PolarizationTransition.cc.

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}
G4double
double G4double
Definition G4Types.hh:83
G4Clebsch::Wigner3J
static G4double Wigner3J(G4double j1, G4double j2, G4double j3, G4double m1, G4double m2, G4double m3)
Definition G4Clebsch.cc:345
G4Clebsch::Wigner9J
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

Referenced by GammaTransF3Coefficient().

◆ FCoefficient()

G4double G4PolarizationTransition::FCoefficient ( G4int K,
G4int L,
G4int Lprime,
G4int twoJ2,
G4int twoJ1 ) const

Definition at line 54 of file G4PolarizationTransition.cc.

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}
G4Clebsch::Wigner6J
static G4double Wigner6J(G4int twoJ1, G4int twoJ2, G4int twoJ3, G4int twoJ4, G4int twoJ5, G4int twoJ6)
Definition G4Clebsch.cc:432

Referenced by GammaTransFCoefficient().

◆ GammaTransF3Coefficient()

G4double G4PolarizationTransition::GammaTransF3Coefficient ( G4int K,
G4int K2,
G4int K1 ) const

Definition at line 91 of file G4PolarizationTransition.cc.

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}
G4PolarizationTransition::F3Coefficient
G4double F3Coefficient(G4int K, G4int K2, G4int K1, G4int L, G4int Lprime, G4int twoJ2, G4int twoJ1) const
Definition G4PolarizationTransition.cc:65

Referenced by SampleGammaTransition().

◆ GammaTransFCoefficient()

G4double G4PolarizationTransition::GammaTransFCoefficient ( G4int K) const

Definition at line 82 of file G4PolarizationTransition.cc.

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}
G4PolarizationTransition::FCoefficient
G4double FCoefficient(G4int K, G4int L, G4int Lprime, G4int twoJ2, G4int twoJ1) const
Definition G4PolarizationTransition.cc:54

◆ SampleGammaTransition()

void G4PolarizationTransition::SampleGammaTransition ( G4NuclearPolarization * np,
G4int twoJ1,
G4int twoJ2,
G4int L0,
G4int Lp,
G4double mpRatio,
G4double & cosTheta,
G4double & phi )

Definition at line 223 of file G4PolarizationTransition.cc.

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}
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition G4Exp.hh:180
G4bool
bool G4bool
Definition G4Types.hh:86
G4complex
std::complex< G4double > G4complex
Definition G4Types.hh:88
G4int
int G4int
Definition G4Types.hh:85
G4UniformRand
#define G4UniformRand()
Definition Randomize.hh:52
G4PolarizationTransition::DumpTransitionData
void DumpTransitionData(const POLAR &pol) const
Definition G4PolarizationTransition.cc:370
G4PolarizationTransition::GammaTransF3Coefficient
G4double GammaTransF3Coefficient(G4int K, G4int K2, G4int K1) const
Definition G4PolarizationTransition.cc:91

◆ SetVerbose()

void G4PolarizationTransition::SetVerbose ( G4int val)
inline

Definition at line 83 of file G4PolarizationTransition.hh.

83{ fVerbose = val; };
The documentation for this class was generated from the following files: