00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #include "globals.hh"
00029 #include "G4ios.hh"
00030 #include "G4HadronicException.hh"
00031 #include "G4Clebsch.hh"
00032 #include "Randomize.hh"
00033 #include "G4Proton.hh"
00034 #include "G4HadTmpUtil.hh"
00035
00036 G4Clebsch::G4Clebsch()
00037 {
00038 G4int nLogs = 101;
00039 logs.push_back(0.);
00040 G4int i;
00041 for (i=1; i<nLogs; i++)
00042 {
00043 G4double previousLog = logs.back();
00044 G4double value = previousLog + std::log((G4double)i);
00045 logs.push_back(value);
00046 }
00047 }
00048
00049
00050 G4Clebsch::~G4Clebsch()
00051 { }
00052
00053
00054 G4bool G4Clebsch::operator==(const G4Clebsch &right) const
00055 {
00056 return (this == (G4Clebsch *) &right);
00057 }
00058
00059
00060 G4bool G4Clebsch::operator!=(const G4Clebsch &right) const
00061 {
00062 return (this != (G4Clebsch *) &right);
00063 }
00064
00065
00066 const std::vector<G4double>& G4Clebsch::GetLogs() const
00067 {
00068 return logs;
00069 }
00070
00071
00072
00073 G4double G4Clebsch::Weight(G4int isoIn1, G4int iso3In1,
00074 G4int isoIn2, G4int iso3In2,
00075 G4int isoOut1, G4int isoOut2) const
00076 {
00077 G4double value = 0.;
00078
00079 G4int an_m = iso3In1 + iso3In2;
00080
00081 G4int jMinIn = std::max(std::abs(isoIn1 - isoIn2), std::abs(an_m));
00082 G4int jMaxIn = isoIn1 + isoIn2;
00083
00084 G4int jMinOut = std::max(std::abs(isoOut1 - isoOut2), std::abs(an_m));
00085 G4int jMaxOut = isoOut1 + isoOut2;
00086
00087 G4int jMin = std::max(jMinIn,jMinOut);
00088 G4int jMax = std::min(jMaxIn,jMaxOut);
00089
00090 G4int j;
00091 for (j=jMin; j<=jMax; j+=2)
00092 {
00093 value += ClebschGordan(isoIn1,iso3In1, isoIn2,iso3In2, j);
00094 }
00095
00096 return value;
00097 }
00098
00099
00100 G4double G4Clebsch::ClebschGordan(G4int isoIn1, G4int iso3In1,
00101 G4int isoIn2, G4int iso3In2,
00102 G4int jOut) const
00103 {
00104
00105
00106 G4double j1 = isoIn1 / 2.0;
00107 G4double j2 = isoIn2 / 2.0;
00108 G4double j3 = jOut / 2.0;
00109
00110 G4double m_1 = iso3In1 / 2.0;
00111 G4double m_2 = iso3In2 / 2.0;
00112 G4double m_3 = - (m_1 + m_2);
00113
00114 G4int n = G4lrint(m_3+j1+j2+.1);
00115 G4double argument = 2. * j3 + 1.;
00116 if (argument < 0.)
00117 throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::ClebschGordan - sqrt of negative argument");
00118 G4double coeff = std::sqrt(argument) / (std::pow(-1.,n));
00119 G4double clebsch = coeff * Wigner3J(j1,j2,j3, m_1,m_2,m_3);
00120 G4double value = clebsch * clebsch;
00121
00122
00123
00124
00125
00126
00127 return value;
00128 }
00129
00130
00131 G4double G4Clebsch::Wigner3J(G4double j1, G4double j2, G4double j3,
00132 G4double m_1, G4double m_2, G4double m_3) const
00133 {
00134
00135
00136 G4double value = 0.;
00137
00138 G4double sigma = j1 + j2 + j3;
00139 std::vector<G4double> n;
00140 n.push_back(-j1 + j2 + j3);
00141 n.push_back(j1 - m_1);
00142 n.push_back(j1 + m_1);
00143 n.push_back(j1 - j2 + j3);
00144 n.push_back(j2 - m_2);
00145 n.push_back(j2 + m_2);
00146 n.push_back(j1 + j2 - j3);
00147 n.push_back(j3 - m_3);
00148 n.push_back(j3 + m_3);
00149
00150
00151
00152 G4bool ok = true;
00153 size_t i;
00154 for(i=1; i<=3; i++)
00155 {
00156 G4double sum1 = n[i-1] + n[i+2] + n[i+5];
00157 G4double sum2 = n[3*i-1] + n[3*i-2] + n[3*i-3];
00158 if (sum1 != sigma || sum2 != sigma) ok = false;
00159 G4int j;
00160 for(j=1; j<=3; j++)
00161 {
00162 if (n[i+3*j-4] < 0.) ok = false;
00163 }
00164 }
00165
00166 if (ok)
00167 {
00168 G4int iMin = 1;
00169 G4int jMin = 1;
00170 G4double smallest = n[0];
00171
00172
00173 for (i=1; i<=3; i++)
00174 {
00175 G4int j;
00176 for (j=1; j<=3; j++)
00177 {
00178 if (n[i+3*j-4] < smallest)
00179 {
00180 smallest = n[i+3*j-4];
00181 iMin = i;
00182 jMin = j;
00183 }
00184 }
00185 }
00186
00187 G4int sign = 1;
00188
00189 if(iMin > 1)
00190 {
00191 for(G4int j=1; j<=3; ++j)
00192 {
00193 G4double tmp = n[j*3-3];
00194 n[j*3-3] = n[iMin+j*3-4];
00195 n[iMin+j*3-4] = tmp;
00196 }
00197 sign = (G4int) std::pow(-1.,sigma);
00198 }
00199
00200 if (jMin > 1)
00201 {
00202 for(i=1; i<=3; i++)
00203 {
00204 G4double tmp = n[i-1];
00205 n[i-1] = n[i+jMin*3-4];
00206 n[i+jMin*3-4] = tmp;
00207 }
00208 sign *= (G4int) std::pow(-1.,sigma);
00209 }
00210
00211 const std::vector<G4double> logVector = GetLogs();
00212 size_t n1 = G4lrint(n[0]);
00213
00214
00215 G4int logEntries = logVector.size() - 1;
00216 for (i=0; i<n.size(); i++)
00217 {
00218 if (n[i] < 0. || n[i] > logEntries)
00219 throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::Wigner3J - Outside logVector boundaries, n");
00220 }
00221
00222 G4double r1 = n[0];
00223 G4double r2 = n[3];
00224 G4double r3 = n[6];
00225 G4double r4 = n[1];
00226 G4double r5 = n[4];
00227 G4double r6 = n[7];
00228 G4double r7 = n[2];
00229 G4double r8 = n[5];
00230 G4double r9 = n[8];
00231
00232 G4double l1 = logVector[(G4int)r1];
00233 G4double l2 = logVector[(G4int)r2];
00234 G4double l3 = logVector[(G4int)r3];
00235 G4double l4 = logVector[(G4int)r4];
00236 G4double l5 = logVector[(G4int)r5];
00237 G4double l6 = logVector[(G4int)r6];
00238 G4double l7 = logVector[(G4int)r7];
00239 G4double l8 = logVector[(G4int)r8];
00240 G4double l9 = logVector[(G4int)r9];
00241
00242 G4double sigma1 = sigma + 1.;
00243 if (sigma1 < 0. || sigma1 > logEntries)
00244 throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::Wigner3J - Outside logVector boundaries, sigma");
00245
00246 G4double ls = logVector[static_cast<G4int>(sigma1+.00001)];
00247 G4double hlp1 = (l2 + l3 + l4 +l7 -ls -l1 -l5 -l9 -l6 -l8) / 2.;
00248 G4int expon = static_cast<G4int>(r6 + r8+.00001);
00249 G4double sgn = std::pow(-1., expon);
00250 G4double coeff = std::exp(hlp1) * sgn;
00251
00252 G4int n61 = static_cast<G4int>(r6 - r1+.00001);
00253 if (n61 < 0. || n61 > logEntries)
00254 throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::Wigner3J - Outside logVector boundaries, n61");
00255 G4int n81 = static_cast<G4int>(r8 - r1+.00001);
00256 if (n81 < 0. || n81 > logEntries)
00257 throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::Wigner3J - Outside logVector boundaries, n81");
00258
00259 G4double hlp2 = l6 - logVector[n61] + l8 - logVector[n81];
00260 G4double sum = std::exp(hlp2);
00261 std::vector<G4double> S;
00262 S.push_back(sum);
00263 n1 = (size_t)r1;
00264 for (i=1; i<=n1; i++)
00265 {
00266 G4double last = S.back();
00267 G4double den = i * (r6 - r1 + i) * (r8 - r1 + i);
00268 if (den == 0)
00269 throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::Wigner3J - divide by zero");
00270 G4double data = -last * (r1 + 1.0 - i) * (r5 + 1.0 - i) * (r9 + 1. - i) / den;
00271 S.push_back(data);
00272 sum += data;
00273 }
00274 value = coeff * sum * sign;
00275 }
00276 else
00277 {
00278 }
00279
00280
00281
00282
00283
00284
00285
00286
00287 return value;
00288 }
00289
00290
00291
00292 std::vector<G4double> G4Clebsch::GenerateIso3(G4int isoIn1, G4int iso3In1,
00293 G4int isoIn2, G4int iso3In2,
00294 G4int isoA, G4int isoB) const
00295 {
00296 std::vector<G4double> temp;
00297
00298
00299
00300
00301 if (isoIn1 == 0 && isoIn2 == 0)
00302 {
00303 G4cout << "WARNING: G4Clebsch::GenerateIso3 - both isoIn are zero" << G4endl;
00304 temp.push_back(0.);
00305 temp.push_back(0.);
00306 return temp;
00307 }
00308
00309 G4int iso3 = iso3In1 + iso3In2;
00310
00311
00312 if (isoA == 0)
00313 {
00314 temp.push_back(0.);
00315 temp.push_back(iso3);
00316 return temp;
00317 }
00318 if (isoB == 0)
00319 {
00320 temp.push_back(iso3);
00321 temp.push_back(0.);
00322 return temp;
00323 }
00324
00325
00326 G4int jMinIn = std::max(std::abs(isoIn1 - isoIn2), std::abs(iso3));
00327 G4int jMaxIn = isoIn1 + isoIn2;
00328
00329
00330
00331 G4int jMinOut = 9999;
00332 G4int jTmp, i, j;
00333
00334 for(i=-1; i<=1; i+=2)
00335 {
00336 for(j=-1; j<=1; j+=2)
00337 {
00338 jTmp= std::abs(i*isoA + j*isoB);
00339 if(jTmp < jMinOut) jMinOut = jTmp;
00340 }
00341 }
00342 jMinOut = std::max(jMinOut, std::abs(iso3));
00343 G4int jMaxOut = isoA + isoB;
00344
00345
00346 G4int jMin = std::max(jMinIn, jMinOut);
00347 G4int jMax = std::min(jMaxIn, jMaxOut);
00348 if (jMin > jMax)
00349 {
00350 throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::GenerateIso3 - jMin > JMax");
00351 }
00352
00353
00354 G4int nJ = (jMax - jMin) / 2 + 1;
00355
00356
00357
00358 if ( (isoIn1 == 0 || isoIn2 == 0) && jMin != jMax )
00359 throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::GenerateIso3 - J1 or J2 = 0, but jMin != JMax");
00360
00361
00362 if (nJ == 0)
00363 throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::GenerateIso3 - nJ is zero, no overlap between in and out");
00364
00365
00366
00367
00368 std::vector<G4double> clebsch;
00369
00370 for(j=jMin; j<=jMax; j+=2)
00371 {
00372 G4double cg = ClebschGordan(isoIn1, iso3In1, isoIn2, iso3In2, j);
00373 clebsch.push_back(cg);
00374 }
00375
00376
00377 if (static_cast<G4int>(clebsch.size()) != nJ)
00378 throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::GenerateIso3 - nJ inconsistency");
00379
00380 G4double sum = clebsch[0];
00381
00382 for (j=1; j<nJ; j++)
00383 {
00384 sum += clebsch[j];
00385 }
00386
00387 if (sum <= 0.)
00388 throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::GenerateIso3 - Sum of Clebsch-Gordan probabilities <=0");
00389
00390
00391
00392 std::vector<G4double> clebschPdf;
00393 G4double previous = clebsch[0];
00394 clebschPdf.push_back(previous/sum);
00395 for (j=1; j<nJ; j++)
00396 {
00397 previous += clebsch[j];
00398 G4double prob = previous / sum;
00399 clebschPdf.push_back(prob);
00400 }
00401
00402
00403 G4double rand = G4UniformRand();
00404 G4int jTot = 0;
00405 for (j=0; j<nJ; j++)
00406 {
00407 G4bool found = false;
00408 if (rand < clebschPdf[j])
00409 {
00410 found = true;
00411 jTot = jMin + 2*j;
00412 }
00413 if (found) break;
00414 }
00415
00416
00417
00418 std::vector<G4double> mMin;
00419 mMin.push_back(-isoA);
00420 mMin.push_back(-isoB);
00421
00422 std::vector<G4double> mMax;
00423 mMax.push_back(isoA);
00424 mMax.push_back(isoB);
00425
00426
00427
00428 std::vector<G4double> m1Out;
00429 std::vector<G4double> m2Out;
00430
00431 const G4int size = 20;
00432 G4double prbout[size][size];
00433
00434 G4int m1pos(0), m2pos(0);
00435 G4int j12;
00436 G4int m1pr(0), m2pr(0);
00437
00438 sum = 0.;
00439 for(j12 = std::abs(isoA-isoB); j12<=(isoA+isoB); j12+=2)
00440 {
00441 m1pos = -1;
00442 for (m1pr = static_cast<G4int>(mMin[0]+.00001); m1pr <= mMax[0]; m1pr+=2)
00443 {
00444 m1pos++;
00445 if (m1pos >= size)
00446 throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::GenerateIso3 - m1pos > size");
00447 m1Out.push_back(m1pr);
00448 m2pos = -1;
00449 for (m2pr = static_cast<G4int>(mMin[1]+.00001); m2pr <= mMax[1]; m2pr+=2)
00450 {
00451 m2pos++;
00452 if (m2pos >= size)
00453 {
00454 throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::GenerateIso3 - m2pos > size");
00455 }
00456 m2Out.push_back(m2pr);
00457
00458 if(m1pr + m2pr == iso3)
00459 {
00460 G4int m12 = m1pr + m2pr;
00461 G4double c12 = ClebschGordan(isoA, m1pr, isoB,m2pr, j12);
00462 G4double c34 = ClebschGordan(0,0,0,0,0);
00463 G4double ctot = ClebschGordan(j12, m12, 0, 0, jTot);
00464 G4double cleb = c12*c34*ctot;
00465 prbout[m1pos][m2pos] = cleb;
00466 sum += cleb;
00467 }
00468 else
00469 {
00470 prbout[m1pos][m2pos] = 0.;
00471 }
00472 }
00473 }
00474 }
00475
00476 if (sum <= 0.)
00477 throw G4HadronicException(__FILE__, __LINE__, "G4Clebsch::GenerateIso3 - sum (out) <=0");
00478
00479 for (i=0; i<size; i++)
00480 {
00481 for (j=0; j<size; j++)
00482 {
00483 prbout[i][j] /= sum;
00484 }
00485 }
00486
00487 rand = G4UniformRand();
00488
00489 G4int m1p, m2p;
00490
00491 for (m1p=0; m1p<m1pos; m1p++)
00492 {
00493 for (m2p=0; m2p<m2pos; m2p++)
00494 {
00495 if (rand < prbout[m1p][m2p])
00496 {
00497 temp.push_back(m1Out[m1p]);
00498 temp.push_back(m2Out[m2p]);
00499 return temp;
00500 }
00501 else
00502 {
00503 rand -= prbout[m1p][m2p];
00504 }
00505 }
00506 }
00507
00508 throw G4HadronicException(__FILE__, __LINE__, "Should never get here ");
00509 return temp;
00510 }
00511
00512
00513 G4double G4Clebsch::NormalizedClebschGordan(G4int J, G4int M,
00514 G4int J1, G4int J2,
00515 G4int m_1, G4int m_2) const
00516 {
00517
00518
00519
00520 G4double cleb = 0.;
00521
00522 if(J1 == 0 || J2 == 0) return cleb;
00523
00524 G4double sum = 0.0;
00525
00526
00527
00528 for(G4int m1Current=-J1; m1Current<=J1; m1Current+=2)
00529 {
00530 G4int m2Current = M - m1Current;
00531
00532 G4double prob = ClebschGordan(J1, m1Current, J2, m2Current, J);
00533 sum += prob;
00534 if (m2Current == m_2 && m1Current == m_1) cleb += prob;
00535 }
00536
00537
00538 if (sum > 0.) cleb /= sum;
00539
00540 return cleb;
00541 }