Geant4: G4GEMProbability.cc Source 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 // $Id$
00027 //
00028 //---------------------------------------------------------------------
00029 //
00030 // Geant4 class G4GEMProbability
00031 //
00032 //
00033 // Hadronic Process: Nuclear De-excitations
00034 // by V. Lara (Sept 2001) 
00035 //
00036 //
00037 // Hadronic Process: Nuclear De-excitations
00038 // by V. Lara (Sept 2001)
00039 //
00040 // J. M. Quesada : several fixes in total GEM width 
00041 // J. M. Quesada 14/07/2009 bug fixed in total emission width (hbarc)
00042 // J. M. Quesada (September 2009) several fixes:
00043 //       -level density parameter of residual calculated at its right excitation energy.
00044 //       -InitialLeveldensity calculated according to the right conditions of the 
00045 //        initial excited nucleus.
00046 // J. M. Quesada 19/04/2010 fix in  emission probability calculation.
00047 // V.Ivanchenko  20/04/2010 added usage of G4Pow and use more safe computation
00048 // V.Ivanchenko 18/05/2010 trying to speedup the most slow method
00049 //            by usage of G4Pow, integer Z and A; moved constructor, 
00050 //            destructor and virtual functions to source
00051 
00052 #include "G4GEMProbability.hh"
00053 #include "G4PairingCorrection.hh"
00054 #include "G4PhysicalConstants.hh"
00055 #include "G4SystemOfUnits.hh"
00056 
00057 G4GEMProbability:: G4GEMProbability(G4int anA, G4int aZ, G4double aSpin) : 
00058   theA(anA), theZ(aZ), Spin(aSpin), theCoulombBarrierPtr(0), 
00059   Normalization(1.0)
00060 {
00061   theEvapLDPptr = new G4EvaporationLevelDensityParameter;
00062   fG4pow = G4Pow::GetInstance(); 
00063   fPlanck= CLHEP::hbar_Planck*fG4pow->logZ(2);
00064   fPairCorr = G4PairingCorrection::GetInstance();
00065 }
00066     
00067 G4GEMProbability::~G4GEMProbability()
00068 {
00069   delete theEvapLDPptr;
00070 }
00071 
00072 G4double G4GEMProbability::EmissionProbability(const G4Fragment & fragment,
00073                                                G4double MaximalKineticEnergy)
00074 {
00075   G4double probability = 0.0;
00076     
00077   if (MaximalKineticEnergy > 0.0 && fragment.GetExcitationEnergy() > 0.0) {
00078     G4double CoulombBarrier = GetCoulombBarrier(fragment);
00079       
00080     probability = 
00081       CalcProbability(fragment,MaximalKineticEnergy,CoulombBarrier);
00082 
00083     // Next there is a loop over excited states for this channel 
00084     // summing probabilities
00085     size_t nn = ExcitEnergies.size();
00086     if (0 < nn) {
00087       G4double SavedSpin = Spin;
00088       for (size_t i = 0; i <nn; ++i) {
00089         Spin = ExcitSpins[i];
00090         // substract excitation energies
00091         G4double Tmax = MaximalKineticEnergy - ExcitEnergies[i];
00092         if (Tmax > 0.0) {
00093           G4double width = CalcProbability(fragment,Tmax,CoulombBarrier);
00094           //JMQ April 2010 added condition to prevent reported crash
00095           // update probability
00096           if (width > 0. && fPlanck < width*ExcitLifetimes[i]) {
00097             probability += width;
00098           }
00099         }
00100       }
00101       // Restore Spin
00102       Spin = SavedSpin;
00103     }
00104   }
00105   Normalization = probability;
00106   return probability;
00107 }
00108 
00109 
00110 G4double G4GEMProbability::CalcProbability(const G4Fragment & fragment, 
00111                                            G4double MaximalKineticEnergy,
00112                                            G4double V)
00113 
00114 // Calculate integrated probability (width) for evaporation channel
00115 {
00116   G4int A = fragment.GetA_asInt();
00117   G4int Z = fragment.GetZ_asInt();
00118 
00119   G4int ResidualA = A - theA;
00120   G4int ResidualZ = Z - theZ;
00121   G4double U = fragment.GetExcitationEnergy();
00122   
00123   G4double NuclearMass = fragment.ComputeGroundStateMass(theZ, theA);
00124   
00125   G4double Alpha = CalcAlphaParam(fragment);
00126   G4double Beta = CalcBetaParam(fragment);
00127     
00128   //                       ***RESIDUAL***
00129   //JMQ (September 2009) the following quantities refer to the RESIDUAL:
00130   
00131   G4double delta0 = fPairCorr->GetPairingCorrection(ResidualA, ResidualZ);  
00132   
00133   G4double a = theEvapLDPptr->
00134     LevelDensityParameter(ResidualA,ResidualZ,MaximalKineticEnergy+V-delta0);
00135   G4double Ux = (2.5 + 150.0/G4double(ResidualA))*MeV;
00136   G4double Ex = Ux + delta0;
00137   G4double T  = 1.0/(std::sqrt(a/Ux) - 1.5/Ux);
00138   //JMQ fixed bug in units
00139   G4double E0 = Ex - T*(std::log(T/MeV) - std::log(a*MeV)/4.0 
00140         - 1.25*std::log(Ux/MeV) + 2.0*std::sqrt(a*Ux));
00141   //                      ***end RESIDUAL ***
00142   
00143   //                       ***PARENT***
00144   //JMQ (September 2009) the following quantities refer to the PARENT:
00145      
00146   G4double deltaCN = fPairCorr->GetPairingCorrection(A, Z);                                    
00147   G4double aCN     = theEvapLDPptr->LevelDensityParameter(A, Z, U-deltaCN);
00148   G4double UxCN    = (2.5 + 150.0/G4double(A))*MeV;
00149   G4double ExCN    = UxCN + deltaCN;
00150   G4double TCN     = 1.0/(std::sqrt(aCN/UxCN) - 1.5/UxCN);
00151   //                     ***end PARENT***
00152 
00153   G4double Width;
00154   G4double InitialLevelDensity;
00155   G4double t = MaximalKineticEnergy/T;
00156   if ( MaximalKineticEnergy < Ex ) {
00157     //JMQ 190709 bug in I1 fixed (T was  missing)
00158     Width = (I1(t,t)*T + (Beta+V)*I0(t))/std::exp(E0/T);
00159     //JMQ 160909 fix:  InitialLevelDensity has been taken away 
00160     //(different conditions for initial CN..) 
00161   } else {
00162 
00163     //VI minor speedup
00164     G4double expE0T = std::exp(E0/T);
00165     const G4double sqrt2 = std::sqrt(2.0);
00166 
00167     G4double tx = Ex/T;
00168     G4double s0 = 2.0*std::sqrt(a*(MaximalKineticEnergy-delta0));
00169     G4double sx = 2.0*std::sqrt(a*(Ex-delta0));
00170     Width = I1(t,tx)*T/expE0T + I3(s0,sx)*std::exp(s0)/(sqrt2*a);
00171     // For charged particles (Beta+V) = 0 beacuse Beta = -V
00172     if (theZ == 0) {
00173       Width += (Beta+V)*(I0(tx)/expE0T + 2.0*sqrt2*I2(s0,sx)*std::exp(s0));
00174     }
00175   }
00176   
00177   //JMQ 14/07/2009 BIG BUG : NuclearMass is in MeV => hbarc instead of hbar_planck must be used
00178   //    G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbar_Planck*hbar_Planck);
00179   G4double gg = (2.0*Spin+1.0)*NuclearMass/(pi2* hbarc*hbarc);
00180   
00181   //JMQ 190709 fix on Rb and  geometrical cross sections according to Furihata's paper 
00182   //                      (JAERI-Data/Code 2001-105, p6)
00183   //    G4double RN = 0.0;
00184   G4double Rb = 0.0;
00185   if (theA > 4) 
00186     {
00187       G4double Ad = fG4pow->Z13(ResidualA);
00188       G4double Aj = fG4pow->Z13(theA);
00189       Rb = 1.12*(Aj + Ad) - 0.86*((Aj+Ad)/(Aj*Ad))+2.85;
00190       Rb *= fermi;
00191     }
00192   else if (theA>1)
00193     {
00194       G4double Ad = fG4pow->Z13(ResidualA);
00195       G4double Aj = fG4pow->Z13(theA);
00196       Rb=1.5*(Aj+Ad)*fermi;
00197     }
00198   else
00199     {
00200       G4double Ad = fG4pow->Z13(ResidualA);
00201       Rb = 1.5*Ad*fermi;
00202     }
00203   //   G4double GeometricalXS = pi*RN*RN*std::pow(ResidualA,2./3.); 
00204   G4double GeometricalXS = pi*Rb*Rb; 
00205   //end of JMQ fix on Rb by 190709
00206   
00207   //JMQ 160909 fix: initial level density must be calculated according to the 
00208   // conditions at the initial compound nucleus 
00209   // (it has been removed from previous "if" for the residual) 
00210 
00211   if ( U < ExCN ) 
00212     {
00213       //JMQ fixed bug in units
00214       //VI moved the computation here
00215       G4double E0CN = ExCN - TCN*(std::log(TCN/MeV) - std::log(aCN*MeV)/4.0 
00216                                   - 1.25*std::log(UxCN/MeV) 
00217                                   + 2.0*std::sqrt(aCN*UxCN));
00218       InitialLevelDensity = (pi/12.0)*std::exp((U-E0CN)/TCN)/TCN;
00219     } 
00220   else 
00221     {
00222       //VI speedup
00223       G4double x  = U-deltaCN;
00224       G4double x1 = std::sqrt(aCN*x);
00225       InitialLevelDensity = (pi/12.0)*std::exp(2*x1)/(x*std::sqrt(x1));
00226     }
00227 
00228   //JMQ 190709 BUG : pi instead of sqrt(pi) must be here according 
00229   // to Furihata's report:
00230   Width *= pi*gg*GeometricalXS*Alpha/(12.0*InitialLevelDensity); 
00231    
00232   return Width;
00233 }
00234 
00235 G4double G4GEMProbability::I3(G4double s0, G4double sx)
00236 {
00237   G4double s2 = s0*s0;
00238   G4double sx2 = sx*sx;
00239   G4double S = 1.0/std::sqrt(s0);
00240   G4double S2 = S*S;
00241   G4double Sx = 1.0/std::sqrt(sx);
00242   G4double Sx2 = Sx*Sx;
00243   
00244   G4double p1 = S *(2.0 + S2 *( 4.0 + S2 *( 13.5 + S2 *( 60.0 + S2 * 325.125 ))));
00245   G4double p2 = Sx*Sx2 *(
00246                          (s2-sx2) + Sx2 *(
00247                                           (1.5*s2+0.5*sx2) + Sx2 *(
00248                                                                    (3.75*s2+0.25*sx2) + Sx2 *(
00249                                                                                               (12.875*s2+0.625*sx2) + Sx2 *(
00250                                                                                                                             (59.0625*s2+0.9375*sx2) + Sx2 *(324.8*s2+3.28*sx2))))));
00251   
00252   p2 *= std::exp(sx-s0);
00253   
00254   return p1-p2; 
00255 }