source: trunk/source/processes/hadronic/models/im_r_matrix/src/G4Clebsch.cc

Last change on this file was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 13.9 KB
Line 
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
36G4Clebsch::G4Clebsch()
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
50G4Clebsch::~G4Clebsch() 
51{  }
52
53
54G4bool G4Clebsch::operator==(const G4Clebsch &right) const
55{
56  return (this == (G4Clebsch *) &right);
57}
58
59
60G4bool G4Clebsch::operator!=(const G4Clebsch &right) const
61{
62  return (this != (G4Clebsch *) &right);
63}
64
65
66const std::vector<G4double>& G4Clebsch::GetLogs() const
67{
68  return logs;
69}
70
71
72
73G4double G4Clebsch::Weight(G4int isoIn1,  G4int iso3In1, 
74                           G4int isoIn2,  G4int iso3In2, 
75                           G4int isoOut1, G4int isoOut2) const
76{
77  G4double value = 0.;
78 
79  G4int m = iso3In1 + iso3In2;
80
81  G4int jMinIn = std::max(std::abs(isoIn1 - isoIn2), std::abs(m));
82  G4int jMaxIn = isoIn1 + isoIn2;
83
84  G4int jMinOut = std::max(std::abs(isoOut1 - isoOut2), std::abs(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
100G4double G4Clebsch::ClebschGordan(G4int isoIn1, G4int iso3In1, 
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 m1 = iso3In1 / 2.0;
111  G4double m2 = iso3In2 / 2.0;
112  G4double m3 = - (m1 + m2);
113
114  G4int n = G4lrint(m3+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, m1,m2,m3);
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
131G4double G4Clebsch::Wigner3J(G4double j1, G4double j2, G4double j3, 
132                             G4double m1, G4double m2, G4double m3) 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 - m1);            // n1
142  n.push_back(j1 + m1);            // n2
143  n.push_back(j1 - j2 + j3);       // n3
144  n.push_back(j2 - m2);            // n4
145  n.push_back(j2 + m2);            // n5
146  n.push_back(j1 + j2 - j3);       // n6
147  n.push_back(j3 - m3);            // n7
148  n.push_back(j3 + m3);            // 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
292std::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
513G4double G4Clebsch::NormalizedClebschGordan(G4int J, G4int m, 
514                                            G4int J1, G4int J2,
515                                            G4int m1, G4int m2) 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 == m2 && m1Current == m1) cleb += prob;
535    }
536
537  // Normalize probs to 1
538  if (sum > 0.) cleb /= sum; 
539
540  return cleb;
541}
Note: See TracBrowser for help on using the repository browser.