Geant4  10.00.p01
G4Clebsch.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 
28 #include "globals.hh"
29 #include "G4ios.hh"
30 #include "G4HadronicException.hh"
31 #include "G4Clebsch.hh"
32 #include "Randomize.hh"
33 #include "G4Proton.hh"
34 #include "G4HadTmpUtil.hh"
35 
37 {
38  G4int nLogs = 101;
39  logs.push_back(0.);
40  G4int i;
41  for (i=1; i<nLogs; i++)
42  {
43  G4double previousLog = logs.back();
44  G4double value = previousLog + std::log((G4double)i);
45  logs.push_back(value);
46  }
47 }
48 
49 
51 { }
52 
53 
55 {
56  return (this == (G4Clebsch *) &right);
57 }
58 
59 
61 {
62  return (this != (G4Clebsch *) &right);
63 }
64 
65 
66 const std::vector<G4double>& G4Clebsch::GetLogs() const
67 {
68  return logs;
69 }
70 
71 
72 
74  G4int isoIn2, G4int iso3In2,
75  G4int isoOut1, G4int isoOut2) const
76 {
77  G4double value = 0.;
78 
79  G4int an_m = iso3In1 + iso3In2;
80 
81  G4int jMinIn = std::max(std::abs(isoIn1 - isoIn2), std::abs(an_m));
82  G4int jMaxIn = isoIn1 + isoIn2;
83 
84  G4int jMinOut = std::max(std::abs(isoOut1 - isoOut2), std::abs(an_m));
85  G4int jMaxOut = isoOut1 + isoOut2;
86 
87  G4int jMin = std::max(jMinIn,jMinOut);
88  G4int jMax = std::min(jMaxIn,jMaxOut);
89 
90  G4int j;
91  for (j=jMin; j<=jMax; j+=2)
92  {
93  value += ClebschGordan(isoIn1,iso3In1, isoIn2,iso3In2, j);
94  }
95 
96  return value;
97 }
98 
99 
101  G4int isoIn2, G4int iso3In2,
102  G4int jOut) const
103 {
104  // Calculates Clebsch-Gordan coefficient
105 
106  G4double j1 = isoIn1 / 2.0;
107  G4double j2 = isoIn2 / 2.0;
108  G4double j3 = jOut / 2.0;
109 
110  G4double m_1 = iso3In1 / 2.0;
111  G4double m_2 = iso3In2 / 2.0;
112  G4double m_3 = - (m_1 + m_2);
113 
114  G4int n = G4lrint(m_3+j1+j2+.1);
115  G4double argument = 2. * j3 + 1.;
116  if (argument < 0.)
117  throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::ClebschGordan - sqrt of negative argument");
118  G4double coeff = std::sqrt(argument) / (std::pow(-1.,n));
119  G4double clebsch = coeff * Wigner3J(j1,j2,j3, m_1,m_2,m_3);
120  G4double value = clebsch * clebsch;
121 
122 // G4cout << "ClebschGordan("
123 // << isoIn1 << "," << iso3In1 << ","
124 // << isoIn2 << "," << iso3In2 << "," << jOut
125 // << ") = " << value << G4endl;
126 
127  return value;
128 }
129 
130 
132  G4double m_1, G4double m_2, G4double m_3) const
133 {
134  // Calculates Wigner 3-j symbols
135 
136  G4double value = 0.;
137 
138  G4double sigma = j1 + j2 + j3;
139  std::vector<G4double> n;
140  n.push_back(-j1 + j2 + j3); // n0
141  n.push_back(j1 - m_1); // n1
142  n.push_back(j1 + m_1); // n2
143  n.push_back(j1 - j2 + j3); // n3
144  n.push_back(j2 - m_2); // n4
145  n.push_back(j2 + m_2); // n5
146  n.push_back(j1 + j2 - j3); // n6
147  n.push_back(j3 - m_3); // n7
148  n.push_back(j3 + m_3); // n8
149 
150  // Some preliminary checks
151 
152  G4bool ok = true;
153  size_t i;
154  for(i=1; i<=3; i++)
155  {
156  G4double sum1 = n[i-1] + n[i+2] + n[i+5];
157  G4double sum2 = n[3*i-1] + n[3*i-2] + n[3*i-3];
158  if (sum1 != sigma || sum2 != sigma) ok = false;
159  G4int j;
160  for(j=1; j<=3; j++)
161  {
162  if (n[i+3*j-4] < 0.) ok = false;
163  }
164  }
165 
166  if (ok)
167  {
168  G4int iMin = 1;
169  G4int jMin = 1;
170  G4double smallest = n[0];
171 
172  // Find the smallest n
173  for (i=1; i<=3; i++)
174  {
175  G4int j;
176  for (j=1; j<=3; j++)
177  {
178  if (n[i+3*j-4] < smallest)
179  {
180  smallest = n[i+3*j-4];
181  iMin = i;
182  jMin = j;
183  }
184  }
185  }
186 
187  G4int sign = 1;
188 
189  if(iMin > 1)
190  {
191  for(G4int j=1; j<=3; ++j)
192  {
193  G4double tmp = n[j*3-3];
194  n[j*3-3] = n[iMin+j*3-4];
195  n[iMin+j*3-4] = tmp;
196  }
197  sign = (G4int) std::pow(-1.,sigma);
198  }
199 
200  if (jMin > 1)
201  {
202  for(i=1; i<=3; i++)
203  {
204  G4double tmp = n[i-1];
205  n[i-1] = n[i+jMin*3-4];
206  n[i+jMin*3-4] = tmp;
207  }
208  sign *= (G4int) std::pow(-1.,sigma);
209  }
210 
211  const std::vector<G4double> logVector = GetLogs();
212  size_t n1 = G4lrint(n[0]);
213 
214  // Some boundary checks
215  G4int logEntries = logVector.size() - 1;
216  for (i=0; i<n.size(); i++)
217  {
218  if (n[i] < 0. || n[i] > logEntries)
219  throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::Wigner3J - Outside logVector boundaries, n");
220  }
221 
222  G4double r1 = n[0];
223  G4double r2 = n[3];
224  G4double r3 = n[6];
225  G4double r4 = n[1];
226  G4double r5 = n[4];
227  G4double r6 = n[7];
228  G4double r7 = n[2];
229  G4double r8 = n[5];
230  G4double r9 = n[8];
231 
232  G4double l1 = logVector[(G4int)r1];
233  G4double l2 = logVector[(G4int)r2];
234  G4double l3 = logVector[(G4int)r3];
235  G4double l4 = logVector[(G4int)r4];
236  G4double l5 = logVector[(G4int)r5];
237  G4double l6 = logVector[(G4int)r6];
238  G4double l7 = logVector[(G4int)r7];
239  G4double l8 = logVector[(G4int)r8];
240  G4double l9 = logVector[(G4int)r9];
241 
242  G4double sigma1 = sigma + 1.;
243  if (sigma1 < 0. || sigma1 > logEntries)
244  throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::Wigner3J - Outside logVector boundaries, sigma");
245 
246  G4double ls = logVector[static_cast<G4int>(sigma1+.00001)];
247  G4double hlp1 = (l2 + l3 + l4 +l7 -ls -l1 -l5 -l9 -l6 -l8) / 2.;
248  G4int expon = static_cast<G4int>(r6 + r8+.00001);
249  G4double sgn = std::pow(-1., expon);
250  G4double coeff = std::exp(hlp1) * sgn;
251 
252  G4int n61 = static_cast<G4int>(r6 - r1+.00001);
253  if (n61 < 0. || n61 > logEntries)
254  throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::Wigner3J - Outside logVector boundaries, n61");
255  G4int n81 = static_cast<G4int>(r8 - r1+.00001);
256  if (n81 < 0. || n81 > logEntries)
257  throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::Wigner3J - Outside logVector boundaries, n81");
258 
259  G4double hlp2 = l6 - logVector[n61] + l8 - logVector[n81];
260  G4double sum = std::exp(hlp2);
261  std::vector<G4double> S;
262  S.push_back(sum);
263  n1 = (size_t)r1;
264  for (i=1; i<=n1; i++)
265  {
266  G4double last = S.back();
267  G4double den = i * (r6 - r1 + i) * (r8 - r1 + i);
268  if (den == 0)
269  throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::Wigner3J - divide by zero");
270  G4double data = -last * (r1 + 1.0 - i) * (r5 + 1.0 - i) * (r9 + 1. - i) / den;
271  S.push_back(data);
272  sum += data;
273  }
274  value = coeff * sum * sign;
275  } // endif ok
276  else
277  {
278  }
279 
280 
281 // G4cout << "Wigner3j("
282 // << j1 << "," << j2 << "," << j3 << ","
283 // << m1 << "," << m2 << "," << m3 << ") = "
284 // << value
285 // << G4endl;
286 
287  return value;
288 }
289 
290 
291 
292 std::vector<G4double> G4Clebsch::GenerateIso3(G4int isoIn1, G4int iso3In1,
293  G4int isoIn2, G4int iso3In2,
294  G4int isoA, G4int isoB) const
295 {
296  std::vector<G4double> temp;
297 
298  // ---- Special cases first ----
299 
300  // Special case, both Jin are zero
301  if (isoIn1 == 0 && isoIn2 == 0)
302  {
303  G4cout << "WARNING: G4Clebsch::GenerateIso3 - both isoIn are zero" << G4endl;
304  temp.push_back(0.);
305  temp.push_back(0.);
306  return temp;
307  }
308 
309  G4int iso3 = iso3In1 + iso3In2;
310 
311  // Special case, either Jout is zero
312  if (isoA == 0)
313  {
314  temp.push_back(0.);
315  temp.push_back(iso3);
316  return temp;
317  }
318  if (isoB == 0)
319  {
320  temp.push_back(iso3);
321  temp.push_back(0.);
322  return temp;
323  }
324 
325  // Number of possible states, in
326  G4int jMinIn = std::max(std::abs(isoIn1 - isoIn2), std::abs(iso3));
327  G4int jMaxIn = isoIn1 + isoIn2;
328 
329  // Number of possible states, out
330 
331  G4int jMinOut = 9999;
332  G4int jTmp, i, j;
333 
334  for(i=-1; i<=1; i+=2)
335  {
336  for(j=-1; j<=1; j+=2)
337  {
338  jTmp= std::abs(i*isoA + j*isoB);
339  if(jTmp < jMinOut) jMinOut = jTmp;
340  }
341  }
342  jMinOut = std::max(jMinOut, std::abs(iso3));
343  G4int jMaxOut = isoA + isoB;
344 
345  // Possible in and out common states
346  G4int jMin = std::max(jMinIn, jMinOut);
347  G4int jMax = std::min(jMaxIn, jMaxOut);
348  if (jMin > jMax)
349  {
350  throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::GenerateIso3 - jMin > JMax");
351  }
352 
353  // Number of possible isospins
354  G4int nJ = (jMax - jMin) / 2 + 1;
355 
356  // A few consistency checks
357 
358  if ( (isoIn1 == 0 || isoIn2 == 0) && jMin != jMax )
359  throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::GenerateIso3 - J1 or J2 = 0, but jMin != JMax");
360 
361  // MGP ---- Shall it be a warning or an exception?
362  if (nJ == 0)
363  throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::GenerateIso3 - nJ is zero, no overlap between in and out");
364 
365  // Loop over all possible combinations of isoIn1, isoIn2, iso3In11, iso3In2, jTot
366  // to get the probability of each of the in-channel couplings
367 
368  std::vector<G4double> clebsch;
369 
370  for(j=jMin; j<=jMax; j+=2)
371  {
372  G4double cg = ClebschGordan(isoIn1, iso3In1, isoIn2, iso3In2, j);
373  clebsch.push_back(cg);
374  }
375 
376  // Consistency check
377  if (static_cast<G4int>(clebsch.size()) != nJ)
378  throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::GenerateIso3 - nJ inconsistency");
379 
380  G4double sum = clebsch[0];
381 
382  for (j=1; j<nJ; j++)
383  {
384  sum += clebsch[j];
385  }
386  // Consistency check
387  if (sum <= 0.)
388  throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::GenerateIso3 - Sum of Clebsch-Gordan probabilities <=0");
389 
390  // Generate a normalized pdf
391 
392  std::vector<G4double> clebschPdf;
393  G4double previous = clebsch[0];
394  clebschPdf.push_back(previous/sum);
395  for (j=1; j<nJ; j++)
396  {
397  previous += clebsch[j];
398  G4double prob = previous / sum;
399  clebschPdf.push_back(prob);
400  }
401 
402  // Generate a random jTot according to the Clebsch-Gordan pdf
403  G4double rand = G4UniformRand();
404  G4int jTot = 0;
405  for (j=0; j<nJ; j++)
406  {
407  G4bool found = false;
408  if (rand < clebschPdf[j])
409  {
410  found = true;
411  jTot = jMin + 2*j;
412  }
413  if (found) break;
414  }
415 
416  // Generate iso3Out
417 
418  std::vector<G4double> mMin;
419  mMin.push_back(-isoA);
420  mMin.push_back(-isoB);
421 
422  std::vector<G4double> mMax;
423  mMax.push_back(isoA);
424  mMax.push_back(isoB);
425 
426  // Calculate the possible |J_i M_i> combinations and their probability
427 
428  std::vector<G4double> m1Out;
429  std::vector<G4double> m2Out;
430 
431  const G4int size = 20;
432  G4double prbout[size][size];
433 
434  G4int m1pos(0), m2pos(0);
435  G4int j12;
436  G4int m1pr(0), m2pr(0);
437 
438  sum = 0.;
439  for(j12 = std::abs(isoA-isoB); j12<=(isoA+isoB); j12+=2)
440  {
441  m1pos = -1;
442  for (m1pr = static_cast<G4int>(mMin[0]+.00001); m1pr <= mMax[0]; m1pr+=2)
443  {
444  m1pos++;
445  if (m1pos >= size)
446  throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::GenerateIso3 - m1pos > size");
447  m1Out.push_back(m1pr);
448  m2pos = -1;
449  for (m2pr = static_cast<G4int>(mMin[1]+.00001); m2pr <= mMax[1]; m2pr+=2)
450  {
451  m2pos++;
452  if (m2pos >= size)
453  {
454  throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::GenerateIso3 - m2pos > size");
455  }
456  m2Out.push_back(m2pr);
457 
458  if(m1pr + m2pr == iso3)
459  {
460  G4int m12 = m1pr + m2pr;
461  G4double c12 = ClebschGordan(isoA, m1pr, isoB,m2pr, j12);
462  G4double c34 = ClebschGordan(0,0,0,0,0);
463  G4double ctot = ClebschGordan(j12, m12, 0, 0, jTot);
464  G4double cleb = c12*c34*ctot;
465  prbout[m1pos][m2pos] = cleb;
466  sum += cleb;
467  }
468  else
469  {
470  prbout[m1pos][m2pos] = 0.;
471  }
472  }
473  }
474  }
475 
476  if (sum <= 0.)
477  throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::GenerateIso3 - sum (out) <=0");
478 
479  for (i=0; i<size; i++)
480  {
481  for (j=0; j<size; j++)
482  {
483  prbout[i][j] /= sum;
484  }
485  }
486 
487  rand = G4UniformRand();
488 
489  G4int m1p, m2p;
490 
491  for (m1p=0; m1p<m1pos; m1p++)
492  {
493  for (m2p=0; m2p<m2pos; m2p++)
494  {
495  if (rand < prbout[m1p][m2p])
496  {
497  temp.push_back(m1Out[m1p]);
498  temp.push_back(m2Out[m2p]);
499  return temp;
500  }
501  else
502  {
503  rand -= prbout[m1p][m2p];
504  }
505  }
506  }
507 
508  throw G4HadronicException(__FILE__, __LINE__, "Should never get here ");
509  return temp;
510 }
511 
512 
514  G4int J1, G4int J2,
515  G4int m_1, G4int m_2) const
516 {
517  // Calculate the normalized Clebsch-Gordan coefficient, that is the prob
518  // of isospin decomposition of (J,m) into J1, J2, m1, m2
519 
520  G4double cleb = 0.;
521 
522  if(J1 == 0 || J2 == 0) return cleb;
523 
524  G4double sum = 0.0;
525 
526  // Loop over all J1,J2,Jtot,m1,m2 combinations
527 
528  for(G4int m1Current=-J1; m1Current<=J1; m1Current+=2)
529  {
530  G4int m2Current = M - m1Current;
531 
532  G4double prob = ClebschGordan(J1, m1Current, J2, m2Current, J);
533  sum += prob;
534  if (m2Current == m_2 && m1Current == m_1) cleb += prob;
535  }
536 
537  // Normalize probs to 1
538  if (sum > 0.) cleb /= sum;
539 
540  return cleb;
541 }
const std::vector< G4double > & GetLogs() const
Definition: G4Clebsch.cc:66
G4double NormalizedClebschGordan(G4int J, G4int m, G4int J1, G4int J2, G4int m1, G4int m2) const
Definition: G4Clebsch.cc:513
int G4int
Definition: G4Types.hh:78
std::vector< G4double > GenerateIso3(G4int isoIn1, G4int iso3In1, G4int isoIn2, G4int iso3In2, G4int isoOut1, G4int isoOut2) const
Definition: G4Clebsch.cc:292
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
const G4int n
virtual ~G4Clebsch()
Definition: G4Clebsch.cc:50
G4double Wigner3J(G4double j1, G4double j2, G4double j3, G4double m1, G4double m2, G4double m3) const
Definition: G4Clebsch.cc:131
int G4lrint(double ad)
Definition: templates.hh:163
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
std::vector< G4double > logs
Definition: G4Clebsch.hh:91
G4double Weight(G4int isoIn1, G4int iso3In1, G4int isoIn2, G4int iso3In2, G4int isoOut1, G4int isoOut2) const
Definition: G4Clebsch.cc:73
#define G4endl
Definition: G4ios.hh:61
G4bool operator!=(const G4Clebsch &right) const
Definition: G4Clebsch.cc:60
G4bool operator==(const G4Clebsch &right) const
Definition: G4Clebsch.cc:54
double G4double
Definition: G4Types.hh:76
G4int sign(const T t)
A simple sign function that allows us to port fortran code to c++ more easily.
G4double ClebschGordan(G4int isoIn1, G4int iso3In1, G4int isoIn2, G4int iso3In2, G4int jOut) const
Definition: G4Clebsch.cc:100