G4Clebsch.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
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   // Calculates Clebsch-Gordan coefficient
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 //   G4cout << "ClebschGordan(" 
00123 //       << isoIn1 << "," << iso3In1 << ","
00124 //       << isoIn2 << "," << iso3In2 << "," << jOut
00125 //       << ") = " << value << G4endl;
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   // Calculates Wigner 3-j symbols
00135 
00136   G4double value = 0.;
00137 
00138   G4double sigma = j1 + j2 + j3;
00139   std::vector<G4double> n;
00140   n.push_back(-j1 + j2 + j3);      // n0
00141   n.push_back(j1 - m_1);            // n1
00142   n.push_back(j1 + m_1);            // n2
00143   n.push_back(j1 - j2 + j3);       // n3
00144   n.push_back(j2 - m_2);            // n4
00145   n.push_back(j2 + m_2);            // n5
00146   n.push_back(j1 + j2 - j3);       // n6
00147   n.push_back(j3 - m_3);            // n7
00148   n.push_back(j3 + m_3);            // n8
00149 
00150   // Some preliminary checks
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     // Find the smallest n
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     // Some boundary checks
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   } // endif ok
00276   else
00277   {
00278   }
00279 
00280 
00281 //  G4cout << "Wigner3j(" 
00282 //       << j1 << "," << j2 << "," << j3 << "," 
00283 //       << m1 << "," << m2 << "," << m3 << ") = " 
00284 //       << value
00285 //       << G4endl;
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   // ---- Special cases first ----
00299 
00300   // Special case, both Jin are zero
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   // Special case, either Jout is zero
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   // Number of possible states, in 
00326   G4int jMinIn = std::max(std::abs(isoIn1 - isoIn2), std::abs(iso3));
00327   G4int jMaxIn = isoIn1 + isoIn2;
00328 
00329   // Number of possible states, out
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    // Possible in and out common states 
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    // Number of possible isospins
00354    G4int nJ = (jMax - jMin) / 2 + 1;
00355 
00356    // A few consistency checks
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    // MGP ---- Shall it be a warning or an exception?
00362    if (nJ == 0)
00363      throw G4HadronicException(__FILE__, __LINE__,  "G4Clebsch::GenerateIso3 - nJ is zero, no overlap between in and out");
00364 
00365    // Loop over all possible combinations of isoIn1, isoIn2, iso3In11, iso3In2, jTot
00366    // to get the probability of each of the in-channel couplings
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    // Consistency check
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    // Consistency check
00387    if (sum <= 0.)
00388      throw G4HadronicException(__FILE__, __LINE__,  "G4Clebsch::GenerateIso3 - Sum of Clebsch-Gordan probabilities <=0");
00389 
00390    // Generate a normalized pdf 
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    // Generate a random jTot according to the Clebsch-Gordan pdf
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    // Generate iso3Out
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    // Calculate the possible |J_i M_i> combinations and their probability
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   // Calculate the normalized Clebsch-Gordan coefficient, that is the prob 
00518   // of isospin decomposition of (J,m) into J1, J2, m1, m2
00519 
00520   G4double cleb = 0.;
00521 
00522   if(J1 == 0 || J2 == 0) return cleb; 
00523   
00524   G4double sum = 0.0;
00525 
00526   // Loop over all J1,J2,Jtot,m1,m2 combinations
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   // Normalize probs to 1 
00538   if (sum > 0.) cleb /= sum; 
00539 
00540   return cleb;
00541 }

Generated on Mon May 27 17:47:54 2013 for Geant4 by  doxygen 1.4.7