Geant4 10.7.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 ()
 
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 387 of file G4PolarizationTransition.cc.

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

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}
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:179
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:387
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: