Changeset 1340 for trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4HETCNeutron.cc
- Timestamp:
- Nov 5, 2010, 3:45:55 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4HETCNeutron.cc
r819 r1340 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4HETCNeutron.cc,v 1.3 2010/08/28 15:16:55 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-ref-09 $ 28 // 29 // by V. Lara 30 // 31 // Modified: 32 // 23.08.2010 V.Ivanchenko general cleanup, move constructor and destructor 33 // the source, use G4Pow 34 26 35 #include "G4HETCNeutron.hh" 36 #include "G4Neutron.hh" 27 37 38 G4HETCNeutron::G4HETCNeutron() 39 : G4HETCFragment(G4Neutron::Neutron(), &theNeutronCoulombBarrier) 40 {} 41 42 G4HETCNeutron::~G4HETCNeutron() 43 {} 44 45 G4double G4HETCNeutron::GetAlpha() 46 { 47 return 0.76+2.2/g4pow->Z13(GetRestA()); 48 } 49 50 G4double G4HETCNeutron::GetBeta() 51 { 52 return (2.12/g4pow->Z23(GetRestA())-0.05)*MeV/GetAlpha(); 53 } 54 55 G4double G4HETCNeutron::GetSpinFactor() 56 { 57 // (2s+1) 58 return 2.0; 59 } 60 28 61 G4double G4HETCNeutron::K(const G4Fragment & aFragment) 29 62 { 30 if (GetStage() != 1) return 1.0; 31 // Number of protons in projectile 32 G4double Pa = static_cast<G4int>(aFragment.GetParticleDefinition()->GetPDGCharge()); 33 // Number of neutrons in projectile 34 G4double Na = aFragment.GetParticleDefinition()->GetBaryonNumber(); 35 G4double TargetA = aFragment.GetA() - Na; 36 G4double TargetZ = aFragment.GetZ() - Pa; 37 Na -= Pa; 38 G4double r = TargetZ/TargetA; 63 // Number of protons in emitted fragment 64 G4int Pa = GetZ(); 65 // Number of neutrons in emitted fragment 66 G4int Na = GetA() - Pa; 67 68 G4int TargetZ = GetRestZ(); 69 G4int TargetA = GetRestA(); 70 G4double r = G4double(TargetZ)/G4double(TargetA); 39 71 40 G4 doubleP = aFragment.GetNumberOfParticles();41 G4 doubleH = aFragment.GetNumberOfHoles();72 G4int P = aFragment.GetNumberOfParticles(); 73 G4int H = aFragment.GetNumberOfHoles(); 42 74 43 75 G4double result = 0.0; 44 76 if (P > 0) 45 77 { 46 result = (H*(1.0-r)+Na)/P; 47 48 result /= (TargetA-TargetZ)/TargetA; 78 result = (H + Na/(1.0-r))/P; 49 79 } 50 80 … … 52 82 } 53 83 54 55 84 G4double G4HETCNeutron::GetKineticEnergy(const G4Fragment & aFragment) 56 85 { 57 G4double H = aFragment.GetNumberOfHoles(); 58 G4double Pb = aFragment.GetNumberOfParticles() - GetA(); 59 G4double Nb = Pb + H; 86 G4int H = aFragment.GetNumberOfHoles(); 87 G4int Pb = aFragment.GetNumberOfParticles(); 88 G4int Nb = Pb + H; 89 G4double g0 = (6.0/pi2)*aFragment.GetA_asInt()*theParameters->GetLevelDensity(); 60 90 61 G4double Ab = std::max(0.0, (Pb*Pb+H*H+Pb-3*H)/4.0);91 G4double Ab = std::max(0.0,G4double(Pb*Pb+H*H+Pb-3*H)/(4.0*g0)); 62 92 G4double Emax = GetMaximalKineticEnergy() - Ab; 63 93 64 G4double cut = GetBeta() / (GetBeta()+Emax/ (Nb+1));94 G4double cut = GetBeta() / (GetBeta()+Emax/G4double(Nb+1)); 65 95 G4double x(0.0); 66 96 if (G4UniformRand() <= cut) 67 97 { 68 x = BetaRand( static_cast<G4int>(Nb),1);98 x = BetaRand(Nb,1); 69 99 } 70 100 else 71 101 { 72 x = BetaRand( static_cast<G4int>(Nb),2);102 x = BetaRand(Nb,2); 73 103 } 74 104
Note: See TracChangeset
for help on using the changeset viewer.