| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4PreCompoundTransitions
33 //
34 // Author: V.Lara
35 //
36 // Modified:
37 // 16.02.2008 J.M.Quesada fixed bugs
38 // 06.09.2008 J.M.Quesada added external choices for:
39 // - "never go back" hipothesis (useNGB=true)
40 // - CEM transition probabilities (useCEMtr=true)
41 // 30.10.2009 J.M.Quesada: CEM transition probabilities have been renormalized
42 // (IAEA benchmark)
43 // 20.08.2010 V.Ivanchenko move constructor and destructor to the source and
44 // optimise the code
45 // 30.08.2011 M.Kelsey - Skip CalculateProbability if no excitons
46
47 #include "G4PreCompoundTransitions.hh"
48 #include "G4PhysicalConstants.hh"
49 #include "G4SystemOfUnits.hh"
50 #include "Randomize.hh"
51 #include "G4NuclearLevelData.hh"
52 #include "G4DeexPrecoParameters.hh"
53 #include "G4Fragment.hh"
54 #include "G4Proton.hh"
55 #include "G4Exp.hh"
56 #include "G4Log.hh"
57
58 G4PreCompoundTransitions::G4PreCompoundTransitions()
59 {
60 proton = G4Proton::Proton();
61 fNuclData = G4NuclearLevelData::GetInstance();
62 G4DeexPrecoParameters* param = fNuclData->GetParameters();
63 FermiEnergy = param->GetFermiEnergy();
64 r0 = param->GetTransitionsR0();
65 }
66
67 // Calculates transition probabilities with
68 // DeltaN = +2 (Trans1) -2 (Trans2) and 0 (Trans3)
69 G4double G4PreCompoundTransitions::
70 CalculateProbability(const G4Fragment & aFragment)
71 {
72 // Number of holes
73 G4int H = aFragment.GetNumberOfHoles();
74 // Number of Particles
75 G4int P = aFragment.GetNumberOfParticles();
76 // Number of Excitons
77 G4int N = P+H;
78 // Nucleus
79 G4int A = aFragment.GetA_asInt();
80 G4int Z = aFragment.GetZ_asInt();
81 G4double U = aFragment.GetExcitationEnergy();
82 TransitionProb2 = 0.0;
83 TransitionProb3 = 0.0;
84 /*
85 G4cout << "G4PreCompoundTransitions::CalculateProbability H/P/N/Z/A= "
86 << H << " " << P << " " << N << " " << Z << " " << A <<G4endl;
87 G4cout << aFragment << G4endl;
88 */
89 if(U < 10*eV || 0==N) { return 0.0; }
90
91 //J. M. Quesada (Feb. 08) new physics
92 // OPT=1 Transitions are calculated according to Gudima's paper
93 // (original in G4PreCompound from VL)
94 // OPT=2 Transitions are calculated according to Gupta's formulae
95 //
96 static const G4double sixdpi2 = 6.0/CLHEP::pi2;
97 G4double GE = sixdpi2*U*fNuclData->GetLevelDensity(Z,A,U);
98 if (useCEMtr) {
99 // Relative Energy (T_{rel})
100 G4double RelativeEnergy = 1.6*FermiEnergy + U/G4double(N);
101
102 // Sample kind of nucleon-projectile
103 G4bool ChargedNucleon(false);
104 if(G4lrint(P*G4UniformRand()) <= aFragment.GetNumberOfCharged()) {
105 ChargedNucleon = true;
106 }
107
108 // Relative Velocity:
109 // <V_{rel}>^2
110 G4double RelativeVelocitySqr;
111 if (ChargedNucleon) {
112 RelativeVelocitySqr = 2*RelativeEnergy/CLHEP::proton_mass_c2;
113 } else {
114 RelativeVelocitySqr = 2*RelativeEnergy/CLHEP::neutron_mass_c2;
115 }
116 // <V_{rel}>
117 G4double RelativeVelocity = std::sqrt(RelativeVelocitySqr);
118
119 // Proton-Proton Cross Section
120 G4double ppXSection =
121 (10.63/RelativeVelocitySqr - 29.92/RelativeVelocity + 42.9)
122 * CLHEP::millibarn;
123 // Proton-Neutron Cross Section
124 G4double npXSection =
125 (34.10/RelativeVelocitySqr - 82.20/RelativeVelocity + 82.2)
126 * CLHEP::millibarn;
127
128 // Averaged Cross Section: \sigma(V_{rel})
129 G4double AveragedXSection;
130 if (ChargedNucleon)
131 {
132 //JMQ: small bug fixed
133 AveragedXSection = ((Z-1)*ppXSection + (A-Z)*npXSection)/G4double(A-1);
134 }
135 else
136 {
137 AveragedXSection = ((A-Z-1)*ppXSection + Z*npXSection)/G4double(A-1);
138 }
139
140 // Fermi relative energy ratio
141 G4double FermiRelRatio = FermiEnergy/RelativeEnergy;
142
143 // This factor is introduced to take into account the Pauli principle
144 G4double PauliFactor = 1.0 - 1.4*FermiRelRatio;
145 if (FermiRelRatio > 0.5) {
146 G4double x = 2.0 - 1.0/FermiRelRatio;
147 PauliFactor += 0.4*FermiRelRatio*x*x*std::sqrt(x);
148 }
149 // Interaction volume
150 G4double xx = 2*r0 + CLHEP::hbarc/(CLHEP::proton_mass_c2*RelativeVelocity);
151 G4double Vint = CLHEP::pi*xx*xx*xx/0.75;
152
153 // Transition probability for \Delta n = +2
154
155 TransitionProb1 = std::max(0.0, AveragedXSection*PauliFactor
156 *std::sqrt(2.0*RelativeEnergy/CLHEP::proton_mass_c2)/Vint);
157
158 //JMQ 281009 phenomenological factor in order to increase
159 // equilibrium contribution
160 // G4double factor=5.0;
161 // TransitionProb1 *= factor;
162
163 // GE = g*E where E is Excitation Energy
164 G4double Fph = G4double(P*P+H*H+P-3*H)*0.25;
165
166 if(!useNGB) {
167
168 // F(p+1,h+1)
169 G4double Fph1 = Fph + N*0.5;
170
171 static const G4double plimit = 100;
172
173 //JMQ/AH bug fixed: if (U-Fph < 0.0)
174 if (GE > Fph1) {
175 G4double x0 = GE-Fph;
176 G4double x1 = (N+1)*G4Log(x0/(GE-Fph1));
177 if(x1 < plimit) {
178 x1 = G4Exp(x1)*TransitionProb1/x0;
179
180 // Transition probability for \Delta n = -2 (at F(p,h) = 0)
181 TransitionProb2 = std::max(0.0, (P*H*(N+1)*(N-2))*x1/x0);
182
183 // Transition probability for \Delta n = 0 (at F(p,h) = 0)
184 TransitionProb3 = std::max(0.0,((N+1)*(P*(P-1) + 4*P*H + H*(H-1)))*x1
185 /static_cast<G4double>(N));
186 }
187 }
188 }
189
190 } else {
191 //JMQ: Transition probabilities from Gupta's work
192 // GE = g*E where E is Excitation Energy
193 TransitionProb1 = std::max(0.0, U*(4.2e+12 - 3.6e+10*U/G4double(N+1)))
194 /(16*CLHEP::c_light);
195
196 if (!useNGB && N > 1) {
197 TransitionProb2 = ((N-1)*(N-2)*P*H)*TransitionProb1/(GE*GE);
198 }
199 }
200 // G4cout<<"U = "<<U<<G4endl;
201 // G4cout<<"N="<<N<<" P="<<P<<" H="<<H<<G4endl;
202 // G4cout<<"l+ ="<<TransitionProb1<<" l- ="<< TransitionProb2
203 // <<" l0 ="<< TransitionProb3<<G4endl;
204 return TransitionProb1 + TransitionProb2 + TransitionProb3;
205 }
206
207 void G4PreCompoundTransitions::PerformTransition(G4Fragment & result)
208 {
209 G4double ChosenTransition =
210 G4UniformRand()*(TransitionProb1 + TransitionProb2 + TransitionProb3);
211 G4int deltaN = 0;
212 G4int Npart = result.GetNumberOfParticles();
213 G4int Ncharged = result.GetNumberOfCharged();
214 G4int Nholes = result.GetNumberOfHoles();
215 if (ChosenTransition <= TransitionProb1)
216 {
217 // Number of excitons is increased on \Delta n = +2
218 deltaN = 2;
219 }
220 else if (ChosenTransition <= TransitionProb1+TransitionProb2)
221 {
222 // Number of excitons is increased on \Delta n = -2
223 deltaN = -2;
224 }
225
226 // AH/JMQ: Randomly decrease the number of charges if deltaN is -2 and
227 // in proportion to the number charges w.r.t. number of particles,
228 // PROVIDED that there are charged particles
229 deltaN /= 2;
230
231 //G4cout << "deltaN= " << deltaN << G4endl;
232
233 // JMQ the following lines have to be before SetNumberOfCharged,
234 // otherwise the check on number of charged vs. number of particles fails
235 result.SetNumberOfParticles(Npart+deltaN);
236 result.SetNumberOfHoles(Nholes+deltaN);
237
238 if(deltaN < 0) {
239 if( (Ncharged == Npart) ||
240 (Ncharged >= 1 && G4int(Npart*G4UniformRand()) <= Ncharged))
241 {
242 result.SetNumberOfCharged(Ncharged+deltaN); // deltaN is negative!
243 }
244
245 } else if ( deltaN > 0 ) {
246 // With weight Z/A, number of charged particles is increased with +1
247 G4int A = result.GetA_asInt() - Npart;
248 G4int Z = result.GetZ_asInt() - Ncharged;
249 if((Z == A) || (Z > 0 && G4lrint(A*G4UniformRand()) <= Z))
250 {
251 result.SetNumberOfCharged(Ncharged+deltaN);
252 }
253 }
254
255 // Number of charged can not be greater that number of particles
256 if ( Npart < Ncharged )
257 {
258 result.SetNumberOfCharged(Npart);
259 }
260 //G4cout << "### After transition" << G4endl;
261 //G4cout << result << G4endl;
262 }
263
264