Geant4 11.1.1
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 ()
 
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 ( )
explicit

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 ( )

Definition at line 54 of file G4PolarizationTransition.cc.

55{}

Member Function Documentation

◆ DumpTransitionData()

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

Definition at line 385 of file G4PolarizationTransition.cc.

386{
387 G4cout << "G4PolarizationTransition: ";
388 (fTwoJ1 % 2) ? G4cout << fTwoJ1 << "/2" : G4cout << fTwoJ1/2;
389 G4cout << " --(" << fLbar;
390 if(fDelta != 0) G4cout << " + " << fDelta << "*" << fL;
391 G4cout << ")--> ";
392 (fTwoJ2 % 2) ? G4cout << fTwoJ2 << "/2" : G4cout << fTwoJ2/2;
393 G4cout << ", P = [ { ";
394 for(std::size_t k=0; k<pol.size(); ++k) {
395 if(k>0) G4cout << " }, { ";
396 for(std::size_t kappa=0; kappa<(pol[k]).size(); ++kappa) {
397 if(kappa > 0) G4cout << ", ";
398 G4cout << (pol[k])[kappa].real() << " + " << (pol[k])[kappa].imag() << "*i";
399 }
400 }
401 G4cout << " } ]" << G4endl;
402}
G4endl
#define G4endl
Definition: G4ios.hh:57
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 68 of file G4PolarizationTransition.cc.

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}
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 57 of file G4PolarizationTransition.cc.

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}
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 94 of file G4PolarizationTransition.cc.

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

Referenced by SampleGammaTransition().

◆ GammaTransFCoefficient()

G4double G4PolarizationTransition::GammaTransFCoefficient ( G4int  K) const

Definition at line 85 of file G4PolarizationTransition.cc.

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

◆ SampleGammaTransition()

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

Definition at line 229 of file G4PolarizationTransition.cc.

234{
235 if(nucpol == nullptr) {
236 if(fVerbose > 1) {
237 G4cout << "G4PolarizationTransition::SampleGammaTransition ERROR: "
238 << "cannot update NULL nuclear polarization" << G4endl;
239 }
240 return;
241 }
242 fTwoJ1 = twoJ1; // add abs to remove negative J
243 fTwoJ2 = twoJ2;
244 fLbar = L0;
245 fL = Lp;
246 fDelta = mpRatio;
247 if(fVerbose > 2) {
248 G4cout << "G4PolarizationTransition: 2J1= " << fTwoJ1 << " 2J2= " << fTwoJ2
249 << " Lbar= " << fLbar << " delta= " << fDelta << " Lp= " << fL
250 << G4endl;
251 G4cout << *nucpol << G4endl;
252 }
253
254 const POLAR& pol = nucpol->GetPolarization();
255
256 if(fTwoJ1 == 0) {
257 nucpol->Unpolarize();
258 cosTheta = 2*G4UniformRand() - 1.0;
259 phi = CLHEP::twopi*G4UniformRand();
260 }
261 else {
262 cosTheta = GenerateGammaCosTheta(pol);
263 if(fVerbose > 2) {
264 G4cout << "G4PolarizationTransition::SampleGammaTransition: cosTheta= "
265 << cosTheta << G4endl;
266 }
267 phi = GenerateGammaPhi(cosTheta, pol);
268 if(fVerbose > 2) {
269 G4cout << "G4PolarizationTransition::SampleGammaTransition: phi= "
270 << phi << G4endl;
271 }
272 }
273
274 if(fTwoJ2 == 0) {
275 nucpol->Unpolarize();
276 return;
277 }
278
279 std::size_t newlength = fTwoJ2+1;
280 //POLAR newPol(newlength);
281 POLAR newPol;
282
283 //map<G4int, map<G4int, G4double> > cache;
284 map<G4int, map<G4int, G4double> >* cachePtr = nullptr;
285 //if(newlength > 10 || pol.size() > 10) cachePtr = &cache;
286
287 for(G4int k2=0; k2<(G4int)newlength; ++k2) {
288 std::vector<G4complex> npolar;
289 npolar.resize(k2+1, 0);
290 //(newPol[k2]).assign(k2+1, 0);
291 for(G4int k1=0; k1<(G4int)pol.size(); ++k1) {
292 for(G4int k=0; k<=k1+k2; k+=2) {
293 // TransF3Coefficient takes the most time. Only calculate it once per
294 // (k, k1, k2) triplet, and wait until the last possible moment to do
295 // so. Break out of the inner loops as soon as it is found to be zero.
296 G4double tF3 = 0.;
297 G4bool recalcTF3 = true;
298 for(G4int kappa2=0; kappa2<=k2; ++kappa2) {
299 G4int ll = (G4int)pol[k1].size();
300 for(G4int kappa1 = 1 - ll; kappa1<ll; ++kappa1) {
301 if(k+k2<k1 || k+k1<k2) continue;
302 G4complex tmpAmp = (kappa1 < 0) ?
303 conj((pol[k1])[-kappa1])*(kappa1 % 2 ? -1.: 1.) : (pol[k1])[kappa1];
304 if(std::abs(tmpAmp) < kEps) continue;
305 G4int kappa = kappa1-(G4int)kappa2;
306 tmpAmp *= G4Clebsch::Wigner3J(2*k1, -2*kappa1, 2*k, 2*kappa, 2*k2, 2*kappa2);
307 if(std::abs(tmpAmp) < kEps) continue;
308 if(recalcTF3) {
309 tF3 = GammaTransF3Coefficient(k, k2, k1);
310 recalcTF3 = false;
311 }
312 if(std::abs(tF3) < kEps) break;
313 tmpAmp *= tF3;
314 if(std::abs(tmpAmp) < kEps) continue;
315 tmpAmp *= ((kappa1+(G4int)k1)%2 ? -1. : 1.)
316 * sqrt((2.*k+1.)*(2.*k1+1.)/(2.*k2+1.));
317 //AR-13Jun2017 Useful for debugging very long computations
318 //G4cout << "G4PolarizationTransition::UpdatePolarizationToFinalState : k1=" << k1
319 // << " ; k2=" << k2 << " ; kappa1=" << kappa1 << " ; kappa2=" << kappa2 << G4endl;
320 tmpAmp *= fgLegendrePolys.EvalAssocLegendrePoly(k, kappa, cosTheta, cachePtr);
321 if(kappa != 0) {
322 tmpAmp *= G4Exp(0.5*(LnFactorial(((G4int)k)-kappa)
323 - LnFactorial(((G4int)k)+kappa)));
324 tmpAmp *= polar(1., phi*kappa);
325 }
326 //(newPol[k2])[kappa2] += tmpAmp;
327 npolar[kappa2] += tmpAmp;
328 }
329 if(!recalcTF3 && std::abs(tF3) < kEps) break;
330 }
331 }
332 }
333 newPol.push_back(npolar);
334 }
335
336 // sanity checks
337 if(fVerbose > 2 && 0.0 == newPol[0][0]) {
338 G4cout << "G4PolarizationTransition::SampleGammaTransition WARNING:"
339 << " P[0][0] is zero!" << G4endl;
340 G4cout << "Old pol is: " << *nucpol << G4endl;
341 DumpTransitionData(newPol);
342 G4cout << "Unpolarizing..." << G4endl;
343 nucpol->Unpolarize();
344 return;
345 }
346 if(fVerbose > 2 && std::abs((newPol[0])[0].imag()) > kEps) {
347 G4cout << "G4PolarizationTransition::SampleGammaTransition WARNING: \n"
348 << " P[0][0] has a non-zero imaginary part! Unpolarizing..."
349 << G4endl;
350 nucpol->Unpolarize();
351 return;
352 }
353 if(fVerbose > 2) {
354 G4cout << "Before normalization: " << G4endl;
355 DumpTransitionData(newPol);
356 }
357 // Normalize and trim
358 std::size_t lastNonEmptyK2 = 0;
359 for(std::size_t k2=0; k2<newlength; ++k2) {
360 G4int lastNonZero = -1;
361 for(std::size_t kappa2=0; kappa2<(newPol[k2]).size(); ++kappa2) {
362 if(k2 == 0 && kappa2 == 0) {
363 lastNonZero = 0;
364 continue;
365 }
366 if(std::abs((newPol[k2])[kappa2]) > 0.0) {
367 lastNonZero = (G4int)kappa2;
368 (newPol[k2])[kappa2] /= (newPol[0])[0];
369 }
370 }
371 while((newPol[k2]).size() != std::size_t (lastNonZero+1)) (newPol[k2]).pop_back();
372 if((newPol[k2]).size() > 0) lastNonEmptyK2 = k2;
373 }
374
375 // Remove zero-value entries
376 while(newPol.size() != lastNonEmptyK2+1) { newPol.pop_back(); }
377 (newPol[0])[0] = 1.0;
378
379 nucpol->SetPolarization(newPol);
380 if(fVerbose > 2) {
381 G4cout << "Updated polarization: " << *nucpol << G4endl;
382 }
383}
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
G4LegendrePolynomial::EvalAssocLegendrePoly
G4double EvalAssocLegendrePoly(G4int l, G4int m, G4double x, std::map< G4int, std::map< G4int, G4double > > *cache=NULL)
Definition: G4LegendrePolynomial.cc:49
G4PolarizationTransition::DumpTransitionData
void DumpTransitionData(const POLAR &pol) const
Definition: G4PolarizationTransition.cc:385
G4PolarizationTransition::GammaTransF3Coefficient
G4double GammaTransF3Coefficient(G4int K, G4int K2, G4int K1) const
Definition: G4PolarizationTransition.cc:94

Referenced by G4GammaTransition::SampleDirection().

◆ SetVerbose()

void G4PolarizationTransition::SetVerbose ( G4int  val)
inline

Definition at line 83 of file G4PolarizationTransition.hh.

83{ fVerbose = val; };

Referenced by G4GammaTransition::SetVerbose().

The documentation for this class was generated from the following files: