G4PreCompoundTransitions.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 // $Id$
00027 //
00028 // -------------------------------------------------------------------
00029 //
00030 // GEANT4 Class file
00031 //
00032 //
00033 // File name:     G4PreCompoundTransitions
00034 //
00035 // Author:         V.Lara
00036 //
00037 // Modified:  
00038 // 16.02.2008 J.M.Quesada fixed bugs 
00039 // 06.09.2008 J.M.Quesada added external choices for:
00040 //                      - "never go back"  hipothesis (useNGB=true) 
00041 //                      -  CEM transition probabilities (useCEMtr=true)
00042 // 30.10.2009 J.M.Quesada: CEM transition probabilities have been renormalized 
00043 //                       (IAEA benchmark)
00044 // 20.08.2010 V.Ivanchenko move constructor and destructor to the source and 
00045 //                         optimise the code
00046 // 30.08.2011 M.Kelsey - Skip CalculateProbability if no excitons
00047 
00048 #include "G4PreCompoundTransitions.hh"
00049 #include "G4PhysicalConstants.hh"
00050 #include "G4SystemOfUnits.hh"
00051 #include "Randomize.hh"
00052 #include "G4Pow.hh"
00053 #include "G4HadronicException.hh"
00054 #include "G4PreCompoundParameters.hh"
00055 #include "G4Proton.hh"
00056 
00057 G4PreCompoundTransitions::G4PreCompoundTransitions() 
00058 {
00059   proton = G4Proton::Proton();
00060   FermiEnergy = G4PreCompoundParameters::GetAddress()->GetFermiEnergy();
00061   r0 = G4PreCompoundParameters::GetAddress()->GetTransitionsr0();
00062   aLDP = G4PreCompoundParameters::GetAddress()->GetLevelDensity();
00063   g4pow = G4Pow::GetInstance();
00064 }
00065 
00066 G4PreCompoundTransitions::~G4PreCompoundTransitions() 
00067 {}
00068 
00069 // Calculates transition probabilities with 
00070 // DeltaN = +2 (Trans1) -2 (Trans2) and 0 (Trans3)
00071 G4double G4PreCompoundTransitions::
00072 CalculateProbability(const G4Fragment & aFragment)
00073 {
00074   // Number of holes
00075   G4int H = aFragment.GetNumberOfHoles();
00076   // Number of Particles 
00077   G4int P = aFragment.GetNumberOfParticles();
00078   // Number of Excitons 
00079   G4int N = P+H;
00080   // Nucleus 
00081   G4int A = aFragment.GetA_asInt();
00082   G4int Z = aFragment.GetZ_asInt();
00083   G4double U = aFragment.GetExcitationEnergy();
00084 
00085   //G4cout << aFragment << G4endl;
00086   
00087   if(U < 10*eV || 0==N) { return 0.0; }
00088   
00089   //J. M. Quesada (Feb. 08) new physics
00090   // OPT=1 Transitions are calculated according to Gudima's paper 
00091   //       (original in G4PreCompound from VL) 
00092   // OPT=2 Transitions are calculated according to Gupta's formulae
00093   //
00094   if (useCEMtr){
00095 
00096     // Relative Energy (T_{rel})
00097     G4double RelativeEnergy = 1.6*FermiEnergy + U/G4double(N);
00098     
00099     // Sample kind of nucleon-projectile 
00100     G4bool ChargedNucleon(false);
00101     G4double chtest = 0.5;
00102     if (P > 0) { 
00103       chtest = G4double(aFragment.GetNumberOfCharged())/G4double(P); 
00104     }
00105     if (G4UniformRand() < chtest) { ChargedNucleon = true; }
00106     
00107     // Relative Velocity: 
00108     // <V_{rel}>^2
00109     G4double RelativeVelocitySqr(0.0);
00110     if (ChargedNucleon) { 
00111       RelativeVelocitySqr = 2.0*RelativeEnergy/CLHEP::proton_mass_c2; 
00112     } else { 
00113       RelativeVelocitySqr = 2.0*RelativeEnergy/CLHEP::neutron_mass_c2; 
00114     }
00115     
00116     // <V_{rel}>
00117     G4double RelativeVelocity = std::sqrt(RelativeVelocitySqr);
00118     
00119     // Proton-Proton Cross Section
00120     G4double ppXSection = 
00121       (10.63/RelativeVelocitySqr - 29.92/RelativeVelocity + 42.9)
00122       * CLHEP::millibarn;
00123     // Proton-Neutron Cross Section
00124     G4double npXSection = 
00125       (34.10/RelativeVelocitySqr - 82.20/RelativeVelocity + 82.2)
00126       * CLHEP::millibarn;
00127     
00128     // Averaged Cross Section: \sigma(V_{rel})
00129     //  G4double AveragedXSection = (ppXSection+npXSection)/2.0;
00130     G4double AveragedXSection(0.0);
00131     if (ChargedNucleon)
00132       {
00133         //JMQ: small bug fixed
00134         //AveragedXSection=((Z-1.0) * ppXSection + (A-Z-1.0) * npXSection)/(A-1.0);
00135         AveragedXSection = ((Z-1)*ppXSection + (A-Z)*npXSection)/G4double(A-1);
00136       }
00137     else 
00138       {
00139         AveragedXSection = ((A-Z-1)*ppXSection + Z*npXSection)/G4double(A-1);
00140         //AveragedXSection = ((A-Z-1)*npXSection + Z*ppXSection)/G4double(A-1);
00141       }
00142     
00143     // Fermi relative energy ratio
00144     G4double FermiRelRatio = FermiEnergy/RelativeEnergy;
00145     
00146     // This factor is introduced to take into account the Pauli principle
00147     G4double PauliFactor = 1.0 - 1.4*FermiRelRatio;
00148     if (FermiRelRatio > 0.5) {
00149       G4double x = 2.0 - 1.0/FermiRelRatio;
00150       PauliFactor += 0.4*FermiRelRatio*x*x*std::sqrt(x);
00151       //PauliFactor += 
00152       //(2.0/5.0)*FermiRelRatio*std::pow(2.0 - (1.0/FermiRelRatio), 5.0/2.0);
00153     }
00154     // Interaction volume 
00155     //  G4double Vint = (4.0/3.0)
00156     //*pi*std::pow(2.0*r0 + hbarc/(proton_mass_c2*RelativeVelocity) , 3.0);
00157     G4double xx = 2.0*r0 + hbarc/(CLHEP::proton_mass_c2*RelativeVelocity);
00158     //    G4double Vint = (4.0/3.0)*CLHEP::pi*xx*xx*xx;
00159     G4double Vint = CLHEP::pi*xx*xx*xx/0.75;
00160     
00161     // Transition probability for \Delta n = +2
00162     
00163     TransitionProb1 = AveragedXSection*PauliFactor
00164       *std::sqrt(2.0*RelativeEnergy/CLHEP::proton_mass_c2)/Vint;
00165 
00166     //JMQ 281009  phenomenological factor in order to increase 
00167     //   equilibrium contribution
00168     //   G4double factor=5.0;
00169     //   TransitionProb1 *= factor;
00170     //
00171     if (TransitionProb1 < 0.0) { TransitionProb1 = 0.0; } 
00172     
00173     // GE = g*E where E is Excitation Energy
00174     G4double GE = (6.0/pi2)*aLDP*A*U;
00175     
00176     //G4double Fph = ((P*P+H*H+P-H)/4.0 - H/2.0);
00177     G4double Fph = G4double(P*P+H*H+P-3*H)/4.0;
00178     
00179     G4bool NeverGoBack(false);
00180     if(useNGB) { NeverGoBack=true; }
00181     
00182     //JMQ/AH  bug fixed: if (U-Fph < 0.0) NeverGoBack = true;
00183     if (GE-Fph < 0.0) { NeverGoBack = true; }
00184     
00185     // F(p+1,h+1)
00186     G4double Fph1 = Fph + N/2.0;
00187     
00188     G4double ProbFactor = g4pow->powN((GE-Fph)/(GE-Fph1),N+1);
00189     
00190     if (NeverGoBack)
00191       {
00192         TransitionProb2 = 0.0;
00193         TransitionProb3 = 0.0;
00194       }
00195     else 
00196       {
00197         // Transition probability for \Delta n = -2 (at F(p,h) = 0)
00198         TransitionProb2 = 
00199           TransitionProb1 * ProbFactor * (P*H*(N+1)*(N-2))/((GE-Fph)*(GE-Fph));
00200         if (TransitionProb2 < 0.0) { TransitionProb2 = 0.0; } 
00201         
00202         // Transition probability for \Delta n = 0 (at F(p,h) = 0)
00203         TransitionProb3 = TransitionProb1*(N+1)* ProbFactor  
00204           * (P*(P-1) + 4.0*P*H + H*(H-1))/(N*(GE-Fph));
00205         if (TransitionProb3 < 0.0) { TransitionProb3 = 0.0; }
00206       }
00207   } else {
00208     //JMQ: Transition probabilities from Gupta's work    
00209     // GE = g*E where E is Excitation Energy
00210     G4double GE = (6.0/pi2)*aLDP*A*U;
00211  
00212     G4double Kmfp=2.;
00213         
00214     //TransitionProb1=1./Kmfp*3./8.*1./c_light*1.0e-9*(1.4e+21*U-2./(N+1)*6.0e+18*U*U);
00215     TransitionProb1 = 3.0e-9*(1.4e+21*U - 1.2e+19*U*U/G4double(N+1))
00216       /(8*Kmfp*CLHEP::c_light);
00217     if (TransitionProb1 < 0.0) { TransitionProb1 = 0.0; }
00218 
00219     TransitionProb2=0.;
00220     TransitionProb3=0.;
00221     
00222     if (!useNGB && N > 1) {
00223       // TransitionProb2=1./Kmfp*3./8.*1./c_light*1.0e-9*(N-1.)*(N-2.)*P*H/(GE*GE)*(1.4e+21*U - 2./(N-1)*6.0e+18*U*U);      
00224       TransitionProb2 = 
00225         3.0e-9*(N-2)*P*H*(1.4e+21*U*(N-1) - 1.2e+19*U*U)/(8*Kmfp*c_light*GE*GE);      
00226       if (TransitionProb2 < 0.0) TransitionProb2 = 0.0; 
00227     }
00228   }
00229   //  G4cout<<"U = "<<U<<G4endl;
00230   //  G4cout<<"N="<<N<<"  P="<<P<<"  H="<<H<<G4endl;
00231   //  G4cout<<"l+ ="<<TransitionProb1<<"  l- ="<< TransitionProb2
00232   //   <<"  l0 ="<< TransitionProb3<<G4endl; 
00233   return TransitionProb1 + TransitionProb2 + TransitionProb3;
00234 }
00235 
00236 void G4PreCompoundTransitions::PerformTransition(G4Fragment & result)
00237 {
00238   G4double ChosenTransition = 
00239     G4UniformRand()*(TransitionProb1 + TransitionProb2 + TransitionProb3);
00240   G4int deltaN = 0;
00241   //  G4int Nexcitons = result.GetNumberOfExcitons();
00242   G4int Npart     = result.GetNumberOfParticles();
00243   G4int Ncharged  = result.GetNumberOfCharged();
00244   G4int Nholes    = result.GetNumberOfHoles();
00245   if (ChosenTransition <= TransitionProb1) 
00246     {
00247       // Number of excitons is increased on \Delta n = +2
00248       deltaN = 2;
00249     } 
00250   else if (ChosenTransition <= TransitionProb1+TransitionProb2) 
00251     {
00252       // Number of excitons is increased on \Delta n = -2
00253       deltaN = -2;
00254     }
00255 
00256   // AH/JMQ: Randomly decrease the number of charges if deltaN is -2 and  
00257   // in proportion to the number charges w.r.t. number of particles,  
00258   // PROVIDED that there are charged particles
00259   deltaN /= 2;
00260 
00261   //G4cout << "deltaN= " << deltaN << G4endl;
00262 
00263   // JMQ the following lines have to be before SetNumberOfCharged, otherwise the check on 
00264   // number of charged vs. number of particles fails
00265   result.SetNumberOfParticles(Npart+deltaN);
00266   result.SetNumberOfHoles(Nholes+deltaN); 
00267 
00268   if(deltaN < 0) {
00269     if( Ncharged >= 1 && G4int(Npart*G4UniformRand()) <= Ncharged) 
00270       { 
00271         result.SetNumberOfCharged(Ncharged+deltaN); // deltaN is negative!
00272       }
00273 
00274   } else if ( deltaN > 0 ) {
00275     // With weight Z/A, number of charged particles is increased with +1
00276     G4int A = result.GetA_asInt();
00277     G4int Z = result.GetZ_asInt();
00278     if( G4int(std::max(1, A - Npart)*G4UniformRand()) <= Z) 
00279       {
00280         result.SetNumberOfCharged(Ncharged+deltaN);
00281       }
00282   }
00283   
00284   // Number of charged can not be greater that number of particles
00285   if ( Npart < Ncharged ) 
00286     {
00287       result.SetNumberOfCharged(Npart);
00288     }
00289   //G4cout << "### After transition" << G4endl;
00290   //G4cout << result << G4endl;
00291 }
00292 

Generated on Mon May 27 17:49:25 2013 for Geant4 by  doxygen 1.4.7