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: G4PreCompoundIon 00034 // 00035 // Author: V.Lara 00036 // 00037 // Modified: 00038 // 10.02.2009 J. M. Quesada fixed bug in density level of light fragments 00039 // 20.08.2010 V.Ivanchenko added G4Pow and G4PreCompoundParameters pointers 00040 // use int Z and A and cleanup 00041 // 00042 00043 #include "G4PreCompoundIon.hh" 00044 #include "G4PhysicalConstants.hh" 00045 00046 G4PreCompoundIon:: 00047 G4PreCompoundIon(const G4ParticleDefinition* part, 00048 G4VCoulombBarrier* aCoulombBarrier) 00049 : G4PreCompoundFragment(part,aCoulombBarrier) 00050 { 00051 G4double r0 = theParameters->Getr0(); 00052 fact = 0.75*CLHEP::millibarn/(CLHEP::pi*r0*r0*r0); 00053 } 00054 00055 G4PreCompoundIon::~G4PreCompoundIon() 00056 {} 00057 00058 G4double G4PreCompoundIon:: 00059 ProbabilityDistributionFunction(G4double eKin, 00060 const G4Fragment& aFragment) 00061 { 00062 if ( !IsItPossible(aFragment) ) { return 0.0; } 00063 G4double efinal = eKin + GetBindingEnergy(); 00064 //G4cout << "Efinal= " << efinal << " Ekin= " << eKin << G4endl; 00065 if(efinal <= 0.0 ) { return 0.0; } 00066 00067 G4double U = aFragment.GetExcitationEnergy(); 00068 G4int P = aFragment.GetNumberOfParticles(); 00069 G4int H = aFragment.GetNumberOfHoles(); 00070 G4int A = GetA(); 00071 G4int N = P + H; 00072 00073 G4double g0 = (6.0/pi2)*aFragment.GetA_asInt()*theParameters->GetLevelDensity(); 00074 G4double g1 = (6.0/pi2)*GetRestA()*theParameters->GetLevelDensity(); 00075 00076 //JMQ 06/02/209 This is THE BUG that was killing cluster emission 00077 // G4double gj = (6.0/pi2)*GetA() * 00078 // G4PreCompoundParameters::GetAddress()->GetLevelDensity(); 00079 00080 G4double gj = g1; 00081 00082 G4double A0 = G4double(P*P+H*H+P-3*H)/(4.0*g0); 00083 G4double A1 = std::max(0.0,(A0*g0 + A*(A-2*P-1)*0.25)/g1); 00084 00085 G4double E0 = U - A0; 00086 //G4cout << "E0= " << E0 << G4endl; 00087 if (E0 <= 0.0) { return 0.0; } 00088 00089 G4double E1 = (std::max(0.0,GetMaximalKineticEnergy() - eKin - A1)); 00090 00091 G4double Aj = A*(A+1)/(4.0*gj); 00092 G4double Ej = std::max(0.0,efinal - Aj); 00093 00094 G4double rj = GetRj(P, aFragment.GetNumberOfCharged()); 00095 G4double xs = CrossSection(eKin); 00096 00097 //G4cout << "rj= " << rj << " xs= " << xs << G4endl; 00098 00099 // JMQ 10/02/09 reshaping of the formula (unnecessary std::pow elimitated) 00100 /* 00101 G4double r0 = theParameters->Getr0(); 00102 G4double pA = (3.0/4.0) * std::sqrt(std::max(0.0, 2.0/(GetReducedMass()* 00103 (eKin+GetBindingEnergy()))))/(pi * r0 * r0 *r0* GetRestA())* 00104 eKin*CrossSection(eKin)*millibarn* 00105 CoalescenceFactor(aFragment.GetA_asInt()) * FactorialFactor(N,P)* 00106 GetRj(aFragment.GetNumberOfParticles(), aFragment.GetNumberOfCharged()); 00107 00108 G4double pB = std::pow((g1*E1)/(g0*E0),N-GetA()-1.0)*(g1/g0); 00109 G4double pC = std::pow((gj*Ej)/(g0*E0),GetA()-1.0)*(gj/g0)/E0; 00110 pA *= pB * pC; 00111 */ 00112 00113 G4double pA = fact*eKin*xs*rj 00114 * CoalescenceFactor(aFragment.GetA_asInt()) * FactorialFactor(N,P) 00115 * std::sqrt(2.0/(GetReducedMass()*efinal)) 00116 * g4pow->powN(g1*E1/(g0*E0), N-A-1) 00117 * g4pow->powN(gj*Ej/(g0*E0), A-1)*gj*g1/(g0*g0*E0*GetRestA()); 00118 00119 return pA; 00120 }