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

Last change on this file since 962 was 819, checked in by garnier, 17 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.