Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
nf_angularMomentumCoupling.cc File Reference
#include <stdlib.h>
#include <cmath>
#include "nf_specialFunctions.h"
Include dependency graph for nf_angularMomentumCoupling.cc:

Go to the source code of this file.

Functions

static int parity (int x)
 
static int max3 (int a, int b, int c)
 
static int max4 (int a, int b, int c, int d)
 
static int min3 (int a, int b, int c)
 
static double w6j0 (int, int *)
 
static double w6j1 (int *)
 
static double cg1 (int, int, int)
 
static double cg2 (int, int, int, int, int, int, int, int)
 
static double cg3 (int, int, int, int, int, int)
 
double nf_amc_log_factorial (int n)
 
double nf_amc_factorial (int n)
 
double nf_amc_wigner_3j (int j1, int j2, int j3, int j4, int j5, int j6)
 
double nf_amc_wigner_6j (int j1, int j2, int j3, int j4, int j5, int j6)
 
double nf_amc_wigner_9j (int j1, int j2, int j3, int j4, int j5, int j6, int j7, int j8, int j9)
 
double nf_amc_racah (int j1, int j2, int l2, int l1, int j3, int l3)
 
double nf_amc_clebsh_gordan (int j1, int j2, int m1, int m2, int j3)
 
double nf_amc_z_coefficient (int l1, int j1, int l2, int j2, int s, int ll)
 
double nf_amc_zbar_coefficient (int l1, int j1, int l2, int j2, int s, int ll)
 
double nf_amc_reduced_matrix_element (int lt, int st, int jt, int l0, int j0, int l1, int j1)
 

Variables

static const int MAX_FACTORIAL = 200
 
static const double nf_amc_log_fact [] = {0.0, 0.0, 0.69314718056, 1.79175946923, 3.17805383035, 4.78749174278, 6.57925121201, 8.52516136107, 10.6046029027, 12.8018274801, 15.1044125731, 17.5023078459, 19.9872144957, 22.5521638531, 25.1912211827, 27.8992713838, 30.6718601061, 33.5050734501, 36.395445208, 39.3398841872, 42.3356164608, 45.3801388985, 48.4711813518, 51.6066755678, 54.7847293981, 58.003605223, 61.261701761, 64.557538627, 67.8897431372, 71.2570389672, 74.6582363488, 78.0922235533, 81.5579594561, 85.0544670176, 88.5808275422, 92.1361756037, 95.7196945421, 99.3306124548, 102.968198615, 106.631760261, 110.320639715, 114.034211781, 117.7718814, 121.533081515, 125.317271149, 129.123933639, 132.952575036, 136.802722637, 140.673923648, 144.565743946, 148.477766952, 152.409592584, 156.360836303, 160.331128217, 164.320112263, 168.327445448, 172.352797139, 176.395848407, 180.456291418, 184.533828861, 188.628173424, 192.739047288, 196.866181673, 201.009316399, 205.168199483, 209.342586753, 213.532241495, 217.736934114, 221.956441819, 226.190548324, 230.439043566, 234.701723443, 238.978389562, 243.268849003, 247.572914096, 251.89040221, 256.22113555, 260.564940972, 264.921649799, 269.291097651, 273.673124286, 278.06757344, 282.474292688, 286.893133295, 291.323950094, 295.766601351, 300.220948647, 304.686856766, 309.16419358, 313.65282995, 318.15263962, 322.663499127, 327.185287704, 331.717887197, 336.261181979, 340.815058871, 345.379407062, 349.954118041, 354.539085519, 359.13420537, 363.739375556, 368.354496072, 372.979468886, 377.614197874, 382.258588773, 386.912549123, 391.575988217, 396.248817052, 400.930948279, 405.622296161, 410.322776527, 415.032306728, 419.7508056, 424.478193418, 429.214391867, 433.959323995, 438.712914186, 443.475088121, 448.245772745, 453.024896238, 457.812387981, 462.608178527, 467.412199572, 472.224383927, 477.044665493, 481.87297923, 486.709261137, 491.553448223, 496.405478487, 501.265290892, 506.132825342, 511.008022665, 515.890824588, 520.781173716, 525.679013516, 530.584288294, 535.49694318, 540.416924106, 545.344177791, 550.278651724, 555.220294147, 560.169054037, 565.124881095, 570.087725725, 575.057539025, 580.034272767, 585.017879389, 590.008311976, 595.005524249, 600.009470555, 605.020105849, 610.037385686, 615.061266207, 620.091704128, 625.128656731, 630.172081848, 635.221937855, 640.27818366, 645.340778693, 650.409682896, 655.484856711, 660.566261076, 665.653857411, 670.747607612, 675.84747404, 680.953419514, 686.065407302, 691.183401114, 696.307365094, 701.437263809, 706.573062246, 711.714725802, 716.862220279, 722.015511874, 727.174567173, 732.339353147, 737.509837142, 742.685986874, 747.867770425, 753.05515623, 758.248113081, 763.446610113, 768.6506168, 773.860102953, 779.07503871, 784.295394535, 789.521141209, 794.752249826, 799.988691789, 805.230438804, 810.477462876, 815.729736304, 820.987231676, 826.249921865, 831.517780024, 836.790779582, 842.068894242, 847.35209797, 852.640365001, 857.933669826, 863.231987192}
 

Function Documentation

static double cg1 ( int  x1,
int  x2,
int  x3 
)
static

Definition at line 339 of file nf_angularMomentumCoupling.cc.

339  {
340 
341  int p1, p2, p3, p4, q1, q2, q3, q4;
342  double a;
343 
344  p1 = x1 + x2 + x3 - 1; if ( ( p1 % 2 ) != 0 ) return( 0.0 );
345  p2 = x1 + x2 - x3;
346  p3 =-x1 + x2 + x3;
347  p4 = x1 - x2 + x3;
348  if ( p2 <= 0 || p3 <= 0 || p4 <= 0 ) return( 0.0 );
349  if ( p1 >= MAX_FACTORIAL ) return( INFINITY );
350 
351  q1 = ( p1 + 1 ) / 2 - 1; p1--;
352  q2 = ( p2 + 1 ) / 2 - 1; p2--;
353  q3 = ( p3 + 1 ) / 2 - 1; p3--;
354  q4 = ( p4 + 1 ) / 2 - 1; p4--;
355 
357  + 0.5 * ( nf_amc_log_fact[ 2 * x3 - 1 ] - nf_amc_log_fact[ 2 * x3 - 2 ]
359 
360  return( ( ( ( q1 + x1 - x2 ) % 2 == 0 ) ? 1.0 : -1.0 ) * G4Exp( a ) );
361 }
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
static const int MAX_FACTORIAL
static const double nf_amc_log_fact[]
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183

Here is the call graph for this function:

Here is the caller graph for this function:

static double cg2 ( int  k,
int  q0,
int  z1,
int  z2,
int  w1,
int  w2,
int  w3,
int  mm 
)
static

Definition at line 365 of file nf_angularMomentumCoupling.cc.

365  {
366 
367  int q1, q2, q3, q4, p1, p2, p3, p4;
368  double a;
369 
370  p1 = z1 + q0 + 2;
371  p2 = z1 - q0 + 1;
372  p3 = z2 + q0 + 1;
373  p4 = -z2 + q0 + 1;
374  if ( p2 <= 0 || p3 <= 0 || p4 <= 0) return( 0.0 );
375  if ( p1 >= MAX_FACTORIAL ) return( INFINITY );
376 
377  q1 = ( p1 + 1 ) / 2 - 1; p1--;
378  q2 = ( p2 + 1 ) / 2 - 1; p2--;
379  q3 = ( p3 + 1 ) / 2 - 1; p3--;
380  q4 = ( p4 + 1 ) / 2 - 1; p4--;
381 
382  a = nf_amc_log_fact[q1] - ( nf_amc_log_fact[ q2 ] + nf_amc_log_fact[ q3 ] + nf_amc_log_fact[ q4 ] )
383  + 0.5 * ( nf_amc_log_fact[ w3 + 1 ] - nf_amc_log_fact[ w3 ]
384  + nf_amc_log_fact[ w1 ] - nf_amc_log_fact[ w1 + 1 ]
385  + nf_amc_log_fact[ w2 ] - nf_amc_log_fact[ w2 + 1 ]
386  + nf_amc_log_fact[ p2 ] + nf_amc_log_fact[ p3 ] + nf_amc_log_fact[ p4 ] - nf_amc_log_fact[ p1 ] );
387 
388  return( ( ( ( q4 + k + ( mm > 0 ) * ( p1 + 2 ) ) % 2 == 0 ) ? -1.0 : 1.0 ) * 2.0 * G4Exp( a ) );
389 }
static constexpr double mm
Definition: G4SIunits.hh:115
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
static const int MAX_FACTORIAL
static const double nf_amc_log_fact[]
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183

Here is the call graph for this function:

Here is the caller graph for this function:

static double cg3 ( int  x1,
int  x2,
int  x3,
int  y1,
int  y2,
int  y3 
)
static

Definition at line 393 of file nf_angularMomentumCoupling.cc.

393  {
394 
395  int nx, i, k1, k2, q1, q2, q3, q4, p1, p2, p3, z1, z2, z3;
396  double a, cg;
397 
398  nx = x1 + x2 + x3 - 1;
399  if ( ( z1 = nx - x1 - y1 ) < 0 ) return( 0.0 );
400  if ( ( z2 = nx - x2 - y2 ) < 0 ) return( 0.0 );
401  if ( ( z3 = nx - x3 - y3 ) < 0 ) return( 0.0 );
402 
403  k1 = x2 - y3;
404  k2 = y1 - x3;
405 
406  q1 = max3( k1, k2, 0 );
407  q2 = min3( y1, x2, z3 + 1 ) - 1;
408  q3 = q1 - k1;
409  q4 = q1 - k2;
410 
411  p1 = y1 - q1 - 1;
412  p2 = x2 - q1 - 1;
413  p3 = z3 - q1;
414 
415  a = cg = G4Exp( 0.5 * ( nf_amc_log_fact[ x3 + y3 - 1 ] - nf_amc_log_fact[ x3 + y3 - 2 ] - nf_amc_log_fact[ nx - 1 ]
416  + nf_amc_log_fact[ z1 ] + nf_amc_log_fact[ z2 ] + nf_amc_log_fact[ z3 ]
417  + nf_amc_log_fact[ x1 - 1 ] + nf_amc_log_fact[ x2 - 1 ] + nf_amc_log_fact[ x3 - 1 ]
418  + nf_amc_log_fact[ y1 - 1 ] + nf_amc_log_fact[ y2 - 1 ] + nf_amc_log_fact[ y3 - 1 ] )
419  - nf_amc_log_fact[ p1 ] - nf_amc_log_fact[ p2 ] - nf_amc_log_fact[ p3 ]
420  - nf_amc_log_fact[ q1 ] - nf_amc_log_fact[ q3 ] - nf_amc_log_fact[ q4 ] ) * ( ( ( q1 % 2 ) == 0 ) ? 1 : -1 );
421  if( cg == INFINITY ) return( INFINITY );
422 
423  if ( q1 != q2 ){
424  q3 = q2 - k1;
425  q4 = q2 - k2;
426  p1 = y1 - q2;
427  p2 = x2 - q2;
428  p3 = z3 - q2 + 1;
429  for( i = 0; i < ( q2 - q1 ); i++ )
430  cg = a - cg * ( ( p1 + i ) * ( p2 + i ) * ( p3 + i ) ) / ( ( q2 - i ) * ( q3 - i ) * ( q4 - i ) );
431  }
432  return( cg );
433 }
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
static int min3(int a, int b, int c)
static const double nf_amc_log_fact[]
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static int max3(int a, int b, int c)

Here is the call graph for this function:

Here is the caller graph for this function:

static int max3 ( int  a,
int  b,
int  c 
)
static

Definition at line 530 of file nf_angularMomentumCoupling.cc.

530  {
531 
532  if( a < b ) a = b;
533  if( a < c ) a = c;
534  return( a );
535 }
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
tuple b
Definition: test.py:12
tuple c
Definition: test.py:13

Here is the caller graph for this function:

static int max4 ( int  a,
int  b,
int  c,
int  d 
)
static

Definition at line 539 of file nf_angularMomentumCoupling.cc.

539  {
540 
541  if( a < b ) a = b;
542  if( a < c ) a = c;
543  if( a < d ) a = d;
544  return( a );
545 }
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
tuple b
Definition: test.py:12
tuple c
Definition: test.py:13

Here is the caller graph for this function:

static int min3 ( int  a,
int  b,
int  c 
)
static

Definition at line 549 of file nf_angularMomentumCoupling.cc.

549  {
550 
551  if( a > b ) a = b;
552  if( a > c ) a = c;
553  return( a );
554 }
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
tuple b
Definition: test.py:12
tuple c
Definition: test.py:13

Here is the caller graph for this function:

double nf_amc_clebsh_gordan ( int  j1,
int  j2,
int  m1,
int  m2,
int  j3 
)

Definition at line 288 of file nf_angularMomentumCoupling.cc.

288  {
289 /*
290 * Clebsh-Gordan coefficient
291 * = <j1,j2,m1,m2|j3,m1+m2>
292 * = (-)^(j1-j2+m1+m2) * std::sqrt(2*j3+1) * / j1 j2 j3 \
293 * \ m1 m2 -m1-m2 /
294 *
295 * Note: Last value m3 is preset to m1+m2. Any other value will evaluate to 0.0.
296 */
297 
298  int m3, x1, x2, x3, y1, y2, y3;
299  double cg = 0.0;
300 
301  if ( j1 < 0 || j2 < 0 || j3 < 0) return( 0.0 );
302  if ( j1 + j2 + j3 > 2 * MAX_FACTORIAL ) return( INFINITY );
303 
304  m3 = m1 + m2;
305 
306  if ( ( x1 = ( j1 + m1 ) / 2 + 1 ) <= 0 ) return( 0.0 );
307  if ( ( x2 = ( j2 + m2 ) / 2 + 1 ) <= 0 ) return( 0.0 );
308  if ( ( x3 = ( j3 - m3 ) / 2 + 1 ) <= 0 ) return( 0.0 );
309 
310  if ( ( y1 = x1 - m1 ) <= 0 ) return( 0.0 );
311  if ( ( y2 = x2 - m2 ) <= 0 ) return( 0.0 );
312  if ( ( y3 = x3 + m3 ) <= 0 ) return( 0.0 );
313 
314  if ( j3 == 0 ){
315  if ( j1 == j2 ) cg = ( 1.0 / std::sqrt( (double)j1 + 1.0 ) * ( ( y1 % 2 == 0 ) ? -1:1 ) );
316  }
317  else if ( (j1 == 0 || j2 == 0 ) ){
318  if ( ( j1 + j2 ) == j3 ) cg = 1.0;
319  }
320  else {
321  if( m3 == 0 && std::abs( m1 ) <= 1 ){
322  if( m1 == 0 ) cg = cg1( x1, x2, x3 );
323  else cg = cg2( x1 + y1 - y2, x3 - 1, x1 + x2 - 2, x1 - y2, j1, j2, j3, m2 );
324  }
325  else if ( m2 == 0 && std::abs( m1 ) <=1 ){
326  cg = cg2( x1 - y2 + y3, x2 - 1, x1 + x3 - 2, x3 - y1, j1, j3, j3, m1 );
327  }
328  else if ( m1 == 0 && std::abs( m3 ) <= 1 ){
329  cg = cg2( x1, x1 - 1, x2 + x3 - 2, x2 - y3, j2, j3, j3, -m3 );
330  }
331  else cg = cg3( x1, x2, x3, y1, y2, y3 );
332  }
333 
334  return( cg );
335 }
static double cg3(int, int, int, int, int, int)
static constexpr double m3
Definition: G4SIunits.hh:131
static double cg1(int, int, int)
static const int MAX_FACTORIAL
static double cg2(int, int, int, int, int, int, int, int)
static constexpr double m2
Definition: G4SIunits.hh:130

Here is the call graph for this function:

Here is the caller graph for this function:

double nf_amc_factorial ( int  n)

Definition at line 96 of file nf_angularMomentumCoupling.cc.

96  {
97 /*
98 * returns n! for pre-computed table. INFINITY is return if n is negative or too large.
99 */
100  return G4Exp( nf_amc_log_factorial( n ) );
101 }
double nf_amc_log_factorial(int)
const G4int n
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183

Here is the call graph for this function:

double nf_amc_log_factorial ( int  n)

Definition at line 85 of file nf_angularMomentumCoupling.cc.

85  {
86 /*
87 * returns ln( n! ).
88 */
89  if( n > MAX_FACTORIAL ) return( INFINITY );
90  if( n < 0 ) return( INFINITY );
91  return nf_amc_log_fact[n];
92 }
static const int MAX_FACTORIAL
static const double nf_amc_log_fact[]
const G4int n

Here is the caller graph for this function:

double nf_amc_racah ( int  j1,
int  j2,
int  l2,
int  l1,
int  j3,
int  l3 
)

Definition at line 253 of file nf_angularMomentumCoupling.cc.

253  {
254 /*
255 * Racah coefficient definition in Edmonds (AR Edmonds, "Angular Momentum in Quantum Mechanics", Princeton (1980) is
256 * W(j1, j2, l2, l1 ; j3, l3) = (-1)^(j1+j2+l1+l2) * { j1 j2 j3 }
257 * { l1 l2 l3 }
258 * The call signature of W(...) appears jumbled, but hey, that's the convention.
259 *
260 * This convention is exactly that used by Blatt-Biedenharn (Rev. Mod. Phys. 24, 258 (1952)) too
261 */
262 
263  double sig;
264 
265  sig = ( ( ( j1 + j2 + l1 + l2 ) % 4 == 0 ) ? 1.0 : -1.0 );
266  return sig * nf_amc_wigner_6j( j1, j2, j3, l1, l2, l3 );
267 }
double nf_amc_wigner_6j(int, int, int, int, int, int)

Here is the call graph for this function:

Here is the caller graph for this function:

double nf_amc_reduced_matrix_element ( int  lt,
int  st,
int  jt,
int  l0,
int  j0,
int  l1,
int  j1 
)

Definition at line 473 of file nf_angularMomentumCoupling.cc.

473  {
474 /*
475 * Reduced Matrix Element for Tensor Operator
476 * = < l1j1 || T(YL,sigma_S)J || l0j0 >
477 *
478 * M.B.Johnson, L.W.Owen, G.R.Satchler
479 * Phys. Rev. 142, 748 (1966)
480 * Note: definition differs from JOS by the factor sqrt(2j1+1)
481 */
482  int llt;
483  double x1, x2, x3, reduced_mat, clebsh_gordan;
484 
485  if ( parity( lt ) != parity( l0 ) * parity( l1 ) ) return( 0.0 );
486  if ( std::abs( l0 - l1 ) > lt || ( l0 + l1 ) < lt ) return( 0.0 );
487  if ( std::abs( ( j0 - j1 ) / 2 ) > jt || ( ( j0 + j1 ) / 2 ) < jt ) return( 0.0 );
488 
489  llt = 2 * lt;
490  jt *= 2;
491  st *= 2;
492 
493  if( ( clebsh_gordan = nf_amc_clebsh_gordan( j1, j0, 1, -1, jt ) ) == INFINITY ) return( INFINITY );
494 
495  reduced_mat = 1.0 / std::sqrt( 4 * M_PI ) * clebsh_gordan / std::sqrt( jt + 1.0 ) /* BRB jt + 1 <= 0? */
496  * std::sqrt( ( j0 + 1.0 ) * ( j1 + 1.0 ) * ( llt + 1.0 ) )
497  * parity( ( j1 - j0 ) / 2 ) * parity( ( -l0 + l1 + lt ) / 2 ) * parity( ( j0 - 1 ) / 2 );
498 
499  if( st == 2 ){
500  x1 = ( l0 - j0 / 2.0 ) * ( j0 + 1.0 );
501  x2 = ( l1 - j1 / 2.0 ) * ( j1 + 1.0 );
502  if ( jt == llt ){
503  x3 = ( lt == 0 ) ? 0 : ( x1 - x2 ) / std::sqrt( lt * ( lt + 1.0 ) );
504  }
505  else if ( jt == ( llt - st ) ){
506  x3 = ( lt == 0 ) ? 0 : -( lt + x1 + x2 ) / std::sqrt( lt * ( 2.0 * lt + 1.0 ) );
507  }
508  else if ( jt == ( llt + st ) ){
509  x3 = ( lt + 1 - x1 - x2 ) / std::sqrt( ( 2.0 * lt + 1.0 ) * ( lt + 1.0 ) );
510  }
511  else{
512  x3 = 1.0;
513  }
514  }
515  else x3 = 1.0;
516  reduced_mat *= x3;
517 
518  return( reduced_mat );
519 }
#define M_PI
Definition: SbMath.h:34
static int parity(int x)
double nf_amc_clebsh_gordan(int, int, int, int, int)

Here is the call graph for this function:

double nf_amc_wigner_3j ( int  j1,
int  j2,
int  j3,
int  j4,
int  j5,
int  j6 
)

Definition at line 105 of file nf_angularMomentumCoupling.cc.

105  {
106 /*
107 * Wigner's 3J symbol (similar to Clebsh-Gordan)
108 * = / j1 j2 j3 \
109 * \ j4 j5 j6 /
110 */
111  double cg;
112 
113  if( ( j4 + j5 + j6 ) != 0 ) return( 0.0 );
114  if( ( cg = nf_amc_clebsh_gordan( j1, j2, j4, j5, j3 ) ) == 0.0 ) return ( 0.0 );
115  if( cg == INFINITY ) return( cg );
116  return( ( ( ( j1 - j2 - j6 ) % 4 == 0 ) ? 1.0 : -1.0 ) * cg / std::sqrt( j3 + 1.0 ) ); /* BRB j3 + 1 <= 0? */
117 }
double nf_amc_clebsh_gordan(int, int, int, int, int)

Here is the call graph for this function:

double nf_amc_wigner_6j ( int  j1,
int  j2,
int  j3,
int  j4,
int  j5,
int  j6 
)

Definition at line 121 of file nf_angularMomentumCoupling.cc.

121  {
122 /*
123 * Wigner's 6J symbol (similar to Racah)
124 * = { j1 j2 j3 }
125 * { j4 j5 j6 }
126 */
127  int i, x[6];
128 
129  x[0] = j1; x[1] = j2; x[2] = j3; x[3] = j4; x[4] = j5; x[5] = j6;
130  for( i = 0; i < 6; i++ ) if ( x[i] == 0 ) return( w6j0( i, x ) );
131 
132  return( w6j1( x ) );
133 }
tuple x
Definition: test.py:50
static double w6j0(int, int *)
static double w6j1(int *)

Here is the call graph for this function:

Here is the caller graph for this function:

double nf_amc_wigner_9j ( int  j1,
int  j2,
int  j3,
int  j4,
int  j5,
int  j6,
int  j7,
int  j8,
int  j9 
)

Definition at line 226 of file nf_angularMomentumCoupling.cc.

226  {
227 /*
228 * Wigner's 9J symbol
229 * / j1 j2 j3 \
230 * = | j4 j5 j6 |
231 * \ j7 j8 j9 /
232 *
233 */
234  int i, i0, i1;
235  double rac;
236 
237  i0 = max3( std::abs( j1 - j9 ), std::abs( j2 - j6 ), std::abs( j4 - j8 ) );
238  i1 = min3( ( j1 + j9 ), ( j2 + j6 ), ( j4 + j8 ) );
239 
240  rac = 0.0;
241  for ( i = i0; i <= i1; i += 2 ){
242  rac += nf_amc_racah( j1, j4, j9, j8, j7, i )
243  * nf_amc_racah( j2, j5, i, j4, j8, j6 )
244  * nf_amc_racah( j9, i, j3, j2, j1, j6 ) * ( i + 1 );
245  if( rac == INFINITY ) return( INFINITY );
246  }
247 
248  return( ( ( (int)( ( j1 + j3 + j5 + j8 ) / 2 + j2 + j4 + j9 ) % 4 == 0 ) ? 1.0 : -1.0 ) * rac );
249 }
static int min3(int a, int b, int c)
double nf_amc_racah(int, int, int, int, int, int)
static int max3(int a, int b, int c)

Here is the call graph for this function:

double nf_amc_z_coefficient ( int  l1,
int  j1,
int  l2,
int  j2,
int  s,
int  ll 
)

Definition at line 437 of file nf_angularMomentumCoupling.cc.

437  {
438 /*
439 * Biedenharn's Z-coefficient coefficient
440 * = Z(l1 j1 l2 j2 | S L )
441 */
442  double z, clebsh_gordan = nf_amc_clebsh_gordan( l1, l2, 0, 0, ll ), racah = nf_amc_racah( l1, j1, l2, j2, s, ll );
443 
444  if( ( clebsh_gordan == INFINITY ) || ( racah == INFINITY ) ) return( INFINITY );
445  z = ( ( ( -l1 + l2 + ll ) % 8 == 0 ) ? 1.0 : -1.0 )
446  * std::sqrt( l1 + 1.0 ) * std::sqrt( l2 + 1.0 ) * std::sqrt( j1 + 1.0 ) * std::sqrt( j2 + 1.0 ) * clebsh_gordan * racah;
447 
448  return( z );
449 }
const XML_Char * s
Definition: expat.h:262
double nf_amc_racah(int, int, int, int, int, int)
tuple z
Definition: test.py:28
double nf_amc_clebsh_gordan(int, int, int, int, int)

Here is the call graph for this function:

double nf_amc_zbar_coefficient ( int  l1,
int  j1,
int  l2,
int  j2,
int  s,
int  ll 
)

Definition at line 453 of file nf_angularMomentumCoupling.cc.

453  {
454 /*
455 * Lane & Thomas's Zbar-coefficient coefficient
456 * = Zbar(l1 j1 l2 j2 | S L )
457 * = (-i)^( -l1 + l2 + ll ) * Z(l1 j1 l2 j2 | S L )
458 *
459 * Lane & Thomas Rev. Mod. Phys. 30, 257-353 (1958).
460 * Note, Lane & Thomas define this because they did not like the different phase convention in Blatt & Biedenharn's Z coefficient. They changed it to get better time-reversal behavior.
461 * Froehner uses Lane & Thomas convention as does T. Kawano.
462 */
463  double zbar, clebsh_gordan = nf_amc_clebsh_gordan( l1, l2, 0, 0, ll ), racah = nf_amc_racah( l1, j1, l2, j2, s, ll );
464 
465  if( ( clebsh_gordan == INFINITY ) || ( racah == INFINITY ) ) return( INFINITY );
466  zbar = std::sqrt( l1 + 1.0 ) * std::sqrt( l2 + 1.0 ) * std::sqrt( j1 + 1.0 ) * std::sqrt( j2 + 1.0 ) * clebsh_gordan * racah;
467 
468  return( zbar );
469 }
const XML_Char * s
Definition: expat.h:262
double nf_amc_racah(int, int, int, int, int, int)
double nf_amc_clebsh_gordan(int, int, int, int, int)

Here is the call graph for this function:

static int parity ( int  x)
static

Definition at line 523 of file nf_angularMomentumCoupling.cc.

523  {
524 
525  return( ( ( x / 2 ) % 2 == 0 ) ? 1 : -1 );
526 }
tuple x
Definition: test.py:50

Here is the caller graph for this function:

static double w6j0 ( int  i,
int x 
)
static

Definition at line 137 of file nf_angularMomentumCoupling.cc.

137  {
138 
139  switch( i ){
140  case 0: if ( ( x[1] != x[2] ) || ( x[4] != x[5] ) ) return( 0.0 );
141  x[5] = x[3]; x[0] = x[1]; x[3] = x[4]; break;
142  case 1: if ( ( x[0] != x[2] ) || ( x[3] != x[5] ) ) return( 0.0 );
143  x[5] = x[4]; break;
144  case 2: if ( ( x[0] != x[1] ) || ( x[3] != x[4] ) ) return( 0.0 );
145  break;
146  //TK fix bug and add comment on 17-05-23
147  //This is the case of 6.3.2 of A. R. Edmonds, Angular Momentum in Quantum Mechanics, Princeton University Press 1974.
148  case 3: if ( ( x[1] != x[5] ) || ( x[2] != x[4] ) ) return( 0.0 );
149  x[5] = x[0]; x[0] = x[4]; x[3] = x[1]; break;
150  case 4: if ( ( x[0] != x[5] ) || ( x[2] != x[3] ) ) return( 0.0 );
151  x[5] = x[1]; break;
152  case 5: if ( ( x[0] != x[4] ) || ( x[1] != x[3] ) ) return( 0.0 );
153  x[5] = x[2]; break;
154  }
155 
156  if( ( x[5] > ( x[0] + x[3] ) ) || ( x[5] < std::abs( x[0] - x[3] ) ) ) return( 0.0 );
157  if( x[0] > MAX_FACTORIAL || x[3] > MAX_FACTORIAL ) { /* BRB Why this test? Why not x[5]? */
158  return( INFINITY );
159  }
160 
161  return( 1.0 / std::sqrt( (double) ( ( x[0] + 1 ) * ( x[3] + 1 ) ) ) * ( ( ( x[0] + x[3] + x[5] ) / 2 ) % 2 != 0 ? -1 : 1 ) );
162 }
tuple x
Definition: test.py:50
static const int MAX_FACTORIAL

Here is the caller graph for this function:

static double w6j1 ( int x)
static

Definition at line 166 of file nf_angularMomentumCoupling.cc.

166  {
167 
168  double w6j, w;
169  int i, k, k1, k2, n, l1, l2, l3, l4, n1, n2, n3, m1, m2, m3, x1, x2, x3, y[4];
170  static int a[3][4] = { { 0, 0, 3, 3},
171  { 1, 4, 1, 4},
172  { 2, 5, 5, 2} };
173 
174  w6j = 0.0;
175 
176  for ( k = 0; k < 4; k++ ){
177  x1 = x[ ( a[0][k] ) ];
178  x2 = x[ ( a[1][k] ) ];
179  x3 = x[ ( a[2][k] ) ];
180 
181  n = ( x1 + x2 + x3 ) / 2;
182  if( n > MAX_FACTORIAL ) {
183  return( INFINITY ); }
184  else if( n < 0 ) {
185  return( 0.0 );
186  }
187 
188  if ( ( n1 = n - x3 ) < 0 ) return( 0.0 );
189  if ( ( n2 = n - x2 ) < 0 ) return( 0.0 );
190  if ( ( n3 = n - x1 ) < 0 ) return( 0.0 );
191 
192  y[k] = n + 2;
193  w6j += nf_amc_log_fact[n1] + nf_amc_log_fact[n2] + nf_amc_log_fact[n3] - nf_amc_log_fact[n+1];
194  }
195 
196  n1 = ( x[0] + x[1] + x[3] + x[4] ) / 2;
197  n2 = ( x[0] + x[2] + x[3] + x[5] ) / 2;
198  n3 = ( x[1] + x[2] + x[4] + x[5] ) / 2;
199 
200  k1 = max4( y[0], y[1], y[2], y[3] ) - 1;
201  k2 = min3( n1, n2, n3 ) + 1;
202 
203  l1 = k1 - y[0] + 1; m1 = n1 - k1 + 1;
204  l2 = k1 - y[1] + 1; m2 = n2 - k1 + 1;
205  l3 = k1 - y[2] + 1; m3 = n3 - k1 + 1;
206  l4 = k1 - y[3] + 1;
207 
208  w6j = w = G4Exp( 0.5 * w6j + nf_amc_log_fact[k1] - nf_amc_log_fact[l1] - nf_amc_log_fact[l2] - nf_amc_log_fact[l3] - nf_amc_log_fact[l4]
209  - nf_amc_log_fact[m1] - nf_amc_log_fact[m2] - nf_amc_log_fact[m3] ) * ( ( k1 % 2 ) == 0 ? -1: 1 );
210  if( w6j == INFINITY ) return( INFINITY );
211 
212  if( k1 != k2 ){
213  k = k2 - k1;
214  m1 -= k-1; m2 -= k-1; m3 -= k-1;
215  l1 += k ; l2 += k ; l3 += k ; l4 += k;
216 
217  for ( i = 0; i < k; i++ )
218  w6j = w - w6j * ( ( k2 - i ) * ( m1 + i ) * ( m2 + i ) * ( m3 + i ) )
219  / ( ( l1 - i ) * ( l2 - i ) * ( l3 - i ) * ( l4 - i ) );
220  }
221  return( w6j );
222 }
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
static int max4(int a, int b, int c, int d)
tuple x
Definition: test.py:50
static constexpr double m3
Definition: G4SIunits.hh:131
static const int MAX_FACTORIAL
static int min3(int a, int b, int c)
static const double nf_amc_log_fact[]
const G4int n
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static constexpr double m2
Definition: G4SIunits.hh:130

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

const int MAX_FACTORIAL = 200
static

Definition at line 68 of file nf_angularMomentumCoupling.cc.

const double nf_amc_log_fact[] = {0.0, 0.0, 0.69314718056, 1.79175946923, 3.17805383035, 4.78749174278, 6.57925121201, 8.52516136107, 10.6046029027, 12.8018274801, 15.1044125731, 17.5023078459, 19.9872144957, 22.5521638531, 25.1912211827, 27.8992713838, 30.6718601061, 33.5050734501, 36.395445208, 39.3398841872, 42.3356164608, 45.3801388985, 48.4711813518, 51.6066755678, 54.7847293981, 58.003605223, 61.261701761, 64.557538627, 67.8897431372, 71.2570389672, 74.6582363488, 78.0922235533, 81.5579594561, 85.0544670176, 88.5808275422, 92.1361756037, 95.7196945421, 99.3306124548, 102.968198615, 106.631760261, 110.320639715, 114.034211781, 117.7718814, 121.533081515, 125.317271149, 129.123933639, 132.952575036, 136.802722637, 140.673923648, 144.565743946, 148.477766952, 152.409592584, 156.360836303, 160.331128217, 164.320112263, 168.327445448, 172.352797139, 176.395848407, 180.456291418, 184.533828861, 188.628173424, 192.739047288, 196.866181673, 201.009316399, 205.168199483, 209.342586753, 213.532241495, 217.736934114, 221.956441819, 226.190548324, 230.439043566, 234.701723443, 238.978389562, 243.268849003, 247.572914096, 251.89040221, 256.22113555, 260.564940972, 264.921649799, 269.291097651, 273.673124286, 278.06757344, 282.474292688, 286.893133295, 291.323950094, 295.766601351, 300.220948647, 304.686856766, 309.16419358, 313.65282995, 318.15263962, 322.663499127, 327.185287704, 331.717887197, 336.261181979, 340.815058871, 345.379407062, 349.954118041, 354.539085519, 359.13420537, 363.739375556, 368.354496072, 372.979468886, 377.614197874, 382.258588773, 386.912549123, 391.575988217, 396.248817052, 400.930948279, 405.622296161, 410.322776527, 415.032306728, 419.7508056, 424.478193418, 429.214391867, 433.959323995, 438.712914186, 443.475088121, 448.245772745, 453.024896238, 457.812387981, 462.608178527, 467.412199572, 472.224383927, 477.044665493, 481.87297923, 486.709261137, 491.553448223, 496.405478487, 501.265290892, 506.132825342, 511.008022665, 515.890824588, 520.781173716, 525.679013516, 530.584288294, 535.49694318, 540.416924106, 545.344177791, 550.278651724, 555.220294147, 560.169054037, 565.124881095, 570.087725725, 575.057539025, 580.034272767, 585.017879389, 590.008311976, 595.005524249, 600.009470555, 605.020105849, 610.037385686, 615.061266207, 620.091704128, 625.128656731, 630.172081848, 635.221937855, 640.27818366, 645.340778693, 650.409682896, 655.484856711, 660.566261076, 665.653857411, 670.747607612, 675.84747404, 680.953419514, 686.065407302, 691.183401114, 696.307365094, 701.437263809, 706.573062246, 711.714725802, 716.862220279, 722.015511874, 727.174567173, 732.339353147, 737.509837142, 742.685986874, 747.867770425, 753.05515623, 758.248113081, 763.446610113, 768.6506168, 773.860102953, 779.07503871, 784.295394535, 789.521141209, 794.752249826, 799.988691789, 805.230438804, 810.477462876, 815.729736304, 820.987231676, 826.249921865, 831.517780024, 836.790779582, 842.068894242, 847.35209797, 852.640365001, 857.933669826, 863.231987192}
static

Definition at line 70 of file nf_angularMomentumCoupling.cc.