Geant4  10.02.p02
G4Fissioner.cc
Go to the documentation of this file.
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 // $Id: G4Fissioner.cc 71954 2013-06-29 04:40:40Z mkelsey $
27 //
28 // 20100114 M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
29 // 20100318 M. Kelsey -- Bug fix setting mass with G4LV
30 // 20100319 M. Kelsey -- Use new generateWithRandomAngles for theta,phi stuff;
31 // eliminate some unnecessary std::pow()
32 // 20100413 M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
33 // 20100517 M. Kelsey -- Inherit from common base class
34 // 20100622 M. Kelsey -- Use local "bindingEnergy()" to call through
35 // 20100711 M. Kelsey -- Add energy-conservation checking, reduce if-cascades
36 // 20100713 M. Kelsey -- Don't add excitation energy to mass (already there)
37 // 20100714 M. Kelsey -- Move conservation checking to base class
38 // 20100728 M. Kelsey -- Make fixed arrays static, move G4FissionStore to data
39 // member and reuse
40 // 20100914 M. Kelsey -- Migrate to integer A and Z
41 // 20110214 M. Kelsey -- Follow G4InuclParticle::Model enumerator migration
42 // 20110801 M. Kelsey -- Replace constant-size std::vector's w/C arrays
43 // 20110922 M. Kelsey -- Follow G4InuclParticle::print(ostream&) migration
44 // 20120517 A. Ribon -- Removed static in local vectors for reproducibility
45 // 20130622 Inherit from G4CascadeDeexciteBase, move to deExcite() interface
46 // with G4Fragment
47 // 20130628 Replace local list of fragments with use of output G4Fragments
48 // 20150608 M. Kelsey -- Label all while loops as terminating.
49 // 20150619 M. Kelsey -- Replace std::exp with G4Exp
50 // 20150622 M. Kelsey -- For new G4cbrt(int), move powers of A outside.
51 
52 #include "G4Fissioner.hh"
53 #include "G4CollisionOutput.hh"
54 #include "G4Exp.hh"
55 #include "G4Fragment.hh"
56 #include "G4HadTmpUtil.hh"
57 #include "G4InuclNuclei.hh"
58 #include "G4InuclParticle.hh"
59 #include "G4FissionStore.hh"
62 
63 using namespace G4InuclSpecialFunctions;
64 
65 
66 void G4Fissioner::deExcite(const G4Fragment& target,
67  G4CollisionOutput& output) {
68  if (verboseLevel) {
69  G4cout << " >>> G4Fissioner::deExcite" << G4endl;
70  }
71 
72  if (verboseLevel > 1)
73  G4cout << " Fissioner input\n" << target << G4endl;
74 
75  // Initialize buffer for fission possibilities
76  fissionStore.setVerboseLevel(verboseLevel);
77  fissionStore.clear();
78 
79  getTargetData(target);
80  G4double A13 = G4cbrt(A);
81  G4double mass_in = PEX.m();
82  G4double e_in = mass_in; // Mass includes excitation
83  G4double PARA = 0.055 * A13*A13 * (G4cbrt(A-Z) + G4cbrt(Z));
84  G4double TEM = std::sqrt(EEXS / PARA);
85  G4double TETA = 0.494 * A13 * TEM;
86 
87  TETA = TETA / std::sinh(TETA);
88 
89  if (A < 246) PARA += (nucleiLevelDensity(A) - PARA) * TETA;
90 
91  G4int A1 = A/2 + 1;
92  G4int Z1;
93  G4int A2 = A - A1;
94 
95  G4double ALMA = -1000.0;
96  G4double DM1 = bindingEnergy(A,Z);
97  G4double EVV = EEXS - DM1;
99  G4double DTEM = (A < 220 ? 0.5 : 1.15);
100 
101  TEM += DTEM;
102 
103  G4double AL1[2] = { -0.15, -0.15 };
104  G4double BET1[2] = { 0.05, 0.05 };
105 
106  G4double R12 = G4cbrt(A1) + G4cbrt(A2);
107 
108  for (G4int i = 0; i < 50 && A1 > 30; i++) {
109  A1--;
110  A2 = A - A1;
111  G4double X3 = 1.0 / G4cbrt(A1);
112  G4double X4 = 1.0 / G4cbrt(A2);
113  Z1 = G4lrint(getZopt(A1, A2, Z, X3, X4, R12) - 1.);
114  G4double EDEF1[2];
115  G4int Z2 = Z - Z1;
116  G4double VPOT, VCOUL;
117 
118  potentialMinimization(VPOT, EDEF1, VCOUL, A1, A2, Z1, Z2, AL1, BET1, R12);
119 
120  G4double DM3 = bindingEnergy(A1,Z1);
121  G4double DM4 = bindingEnergyAsymptotic(A1, Z1);
122  G4double DM5 = bindingEnergy(A2,Z2);
123  G4double DM6 = bindingEnergyAsymptotic(A2, Z2);
124  G4double DMT1 = DM4 + DM6 - DM2;
125  G4double DMT = DM3 + DM5 - DM1;
126  G4double EZL = EEXS + DMT - VPOT;
127 
128  if(EZL > 0.0) { // generate fluctuations
129  // faster, using randomGauss
130  G4double C1 = std::sqrt(getC2(A1, A2, X3, X4, R12) / TEM);
131  G4double DZ = randomGauss(C1);
132 
133  DZ = DZ > 0.0 ? DZ + 0.5 : -std::fabs(DZ - 0.5);
134  Z1 += G4int(DZ);
135  Z2 -= G4int(DZ);
136 
137  G4double DEfin = randomGauss(TEM);
138  G4double EZ = (DMT1 + (DMT - DMT1) * TETA - VPOT + DEfin) / TEM;
139 
140  if (EZ >= ALMA) ALMA = EZ;
141  G4double EK = VCOUL + DEfin + 0.5 * TEM;
142  G4double EV = EVV + bindingEnergy(A1,Z1) + bindingEnergy(A2,Z2) - EK;
143 
144  if (EV > 0.0) fissionStore.addConfig(A1, Z1, EZ, EK, EV);
145  };
146  };
147 
148  G4int store_size = fissionStore.size();
149  if (store_size == 0) return; // No fission products
150 
152  fissionStore.generateConfiguration(ALMA, inuclRndm());
153 
154  A1 = G4int(config.afirst);
155  A2 = A - A1;
156  Z1 = G4int(config.zfirst);
157 
158  G4int Z2 = Z - Z1;
159 
160  G4double mass1 = G4InuclNuclei::getNucleiMass(A1,Z1);
161  G4double mass2 = G4InuclNuclei::getNucleiMass(A2,Z2);
162  G4double EK = config.ekin;
163  G4double pmod = std::sqrt(0.001 * EK * mass1 * mass2 / mass_in);
164 
165  G4LorentzVector mom1 = generateWithRandomAngles(pmod, mass1);
166  G4LorentzVector mom2; mom2.setVectM(-mom1.vect(), mass2);
167 
168  G4double e_out = mom1.e() + mom2.e();
169  G4double EV = 1000.0 * (e_in - e_out) / A;
170  if (EV <= 0.0) return; // No fission energy
171 
172  G4double EEXS1 = EV*A1;
173  G4double EEXS2 = EV*A2;
174 
175  // Pass only last two nuclear fragments
176  output.addRecoilFragment(makeFragment(mom1, A1, Z1, EEXS1));
177  output.addRecoilFragment(makeFragment(mom2, A2, Z2, EEXS2));
178 }
179 
181  G4int A2,
182  G4double X3,
183  G4double X4,
184  G4double R12) const {
185 
186  if (verboseLevel > 3) {
187  G4cout << " >>> G4Fissioner::getC2" << G4endl;
188  }
189 
190  G4double C2 = 124.57 * (1.0 / A1 + 1.0 / A2) + 0.78 * (X3 + X4) - 176.9 *
191  ((X3*X3*X3*X3) + (X4*X4*X4*X4)) + 219.36 * (1.0 / (A1 * A1) + 1.0 / (A2 * A2)) - 1.108 / R12;
192 
193  return C2;
194 }
195 
197  G4int A2,
198  G4int ZT,
199  G4double X3,
200  G4double X4,
201  G4double R12) const {
202 
203  if (verboseLevel > 3) {
204  G4cout << " >>> G4Fissioner::getZopt" << G4endl;
205  }
206 
207  G4double Zopt = (87.7 * (X4 - X3) * (1.0 - 1.25 * (X4 + X3)) +
208  ZT * ((124.57 / A2 + 0.78 * X4 - 176.9 * (X4*X4*X4*X4) + 219.36 / (A2 * A2)) - 0.554 / R12)) /
209  getC2(A1, A2, X3, X4, R12);
210 
211  return Zopt;
212 }
213 
215  G4double( &ED)[2],
216  G4double& VC,
217  G4int AF,
218  G4int AS,
219  G4int ZF,
220  G4int ZS,
221  G4double AL1[2],
222  G4double BET1[2],
223  G4double& R12) const {
224 
225  if (verboseLevel > 3) {
226  G4cout << " >>> G4Fissioner::potentialMinimization" << G4endl;
227  }
228 
229  const G4double huge_num = 2.0e35;
230  const G4int itry_max = 2000;
231  const G4double DSOL1 = 1.0e-6;
232  const G4double DS1 = 0.3;
233  const G4double DS2 = 1.0 / DS1 / DS1;
234  G4int A1[2] = { AF, AS };
235  G4int Z1[2] = { ZF, ZS };
236  G4double D = 1.01844 * ZF * ZS;
237  G4double D0 = 1.0e-3 * D;
238  G4double R[2];
239  R12 = 0.0;
240  G4double C[2];
241  G4double F[2];
242  G4double Y1;
243  G4double Y2;
244  G4int i;
245 
246  for (i = 0; i < 2; i++) {
247  R[i] = G4cbrt(A1[i]);
248  Y1 = R[i] * R[i];
249  Y2 = Z1[i] * Z1[i] / R[i];
250  C[i] = 6.8 * Y1 - 0.142 * Y2;
251  F[i] = 12.138 * Y1 - 0.145 * Y2;
252  };
253 
254  G4double SAL[2];
255  G4double SBE[2];
256  G4double X[2];
257  G4double X1[2];
258  G4double X2[2];
259  G4double RAL[2];
260  G4double RBE[2];
261  G4double AA[4][4];
262  G4double B[4];
263  G4int itry = 0;
264 
265  while (itry < itry_max) { /* Loop checking 08.06.2015 MHK */
266  itry++;
267  G4double S = 0.0;
268 
269  for (i = 0; i < 2; i++) {
270  S += R[i] * (1.0 + AL1[i] + BET1[i] - 0.257 * AL1[i] * BET1[i]);
271  };
272  R12 = 0.0;
273  Y1 = 0.0;
274  Y2 = 0.0;
275 
276  for (i = 0; i < 2; i++) {
277  SAL[i] = R[i] * (1.0-0.257 * BET1[i]);
278  SBE[i] = R[i] * (1.0-0.257 * AL1[i]);
279  X[i] = R[i] / S;
280  X1[i] = X[i] * X[i];
281  X2[i] = X[i] * X1[i];
282  Y1 += AL1[i] * X1[i];
283  Y2 += BET1[i] * X2[i];
284  R12 += R[i] * (1.0 - AL1[i] * (1.0 - 0.6 * X[i]) + BET1[i] * (1.0 - 0.429 * X1[i]));
285  };
286 
287  G4double Y3 = -0.6 * Y1 + 0.857 * Y2;
288  G4double Y4 = (1.2 * Y1 - 2.571 * Y2) / S;
289  G4double R2 = D0 / (R12 * R12);
290  G4double R3 = 2.0 * R2 / R12;
291 
292  for (i = 0; i < 2; i++) {
293  RAL[i] = -R[i] * (1.0 - 0.6 * X[i]) + SAL[i] * Y3;
294  RBE[i] = R[i] * (1.0 - 0.429 * X1[i]) + SBE[i] * Y3;
295  };
296 
297  G4double DX1;
298  G4double DX2;
299 
300  for (i = 0; i < 2; i++) {
301 
302  for (G4int j = 0; j < 2; j++) {
303  G4double DEL1 = i == j ? 1.0 : 0.0;
304  DX1 = 0.0;
305  DX2 = 0.0;
306 
307  if (std::fabs(AL1[i]) >= DS1) {
308  G4double XXX = AL1[i] * AL1[i] * DS2;
309  G4double DEX = XXX > 100.0 ? huge_num : G4Exp(XXX);
310  DX1 = 2.0 * (1.0 + 2.0 * AL1[i] * AL1[i] * DS2) * DEX * DS2;
311  };
312 
313  if (std::fabs(BET1[i]) >= DS1) {
314  G4double XXX = BET1[i] * BET1[i] * DS2;
315  G4double DEX = XXX > 100.0 ? huge_num : G4Exp(XXX);
316  DX2 = 2.0 * (1.+2.0 * BET1[i] * BET1[i] * DS2) * DEX * DS2;
317  };
318 
319  G4double DEL = 2.0e-3 * DEL1;
320  AA[i][j] = R3 * RBE[i] * RBE[j] -
321  R2 * (-0.6 *
322  (X1[i] * SAL[j] +
323  X1[j] * SAL[i]) + SAL[i] * SAL[j] * Y4) +
324  DEL * C[i] + DEL1 * DX1;
325  G4int i1 = i + 2;
326  G4int j1 = j + 2;
327  AA[i1][j1] = R3 * RBE[i] * RBE[j]
328  - R2 * (0.857 *
329  (X2[i] * SBE[j] +
330  X2[j] * SBE[i]) + SBE[i] * SBE[j] * Y4) +
331  DEL * F[i] + DEL1 * DX2;
332  AA[i][j1] = R3 * RAL[i] * RBE[j] -
333  R2 * (0.857 *
334  (X2[j] * SAL[i] -
335  0.6 * X1[i] * SBE[j]) + SBE[j] * SAL[i] * Y4 -
336  0.257 * R[i] * Y3 * DEL1);
337  AA[j1][i] = AA[i][j1];
338  };
339  };
340 
341  for (i = 0; i < 2; i++) {
342  DX1 = 0.0;
343  DX2 = 0.0;
344 
345  if (std::fabs(AL1[i]) >= DS1) DX1 = 2.0 * AL1[i] * DS2 * G4Exp(AL1[i] * AL1[i] * DS2);
346 
347  if (std::fabs(BET1[i]) >= DS1) DX2 = 2.0 * BET1[i] * DS2 * G4Exp(BET1[i] * BET1[i] * DS2);
348  B[i] = R2 * RAL[i] - 2.0e-3 * C[i] * AL1[i] + DX1;
349  B[i + 2] = R2 * RBE[i] - 2.0e-3 * F[i] * BET1[i] + DX2;
350  };
351 
352  G4double ST = 0.0;
353  G4double ST1 = 0.0;
354 
355  for (i = 0; i < 4; i++) {
356  ST += B[i] * B[i];
357 
358  for (G4int j = 0; j < 4; j++) ST1 += AA[i][j] * B[i] * B[j];
359  };
360 
361  G4double STEP = ST / ST1;
362  G4double DSOL = 0.0;
363 
364  for (i = 0; i < 2; i++) {
365  AL1[i] += B[i] * STEP;
366  BET1[i] += B[i + 2] * STEP;
367  DSOL += B[i] * B[i] + B[i + 2] * B[i + 2];
368  };
369  DSOL = std::sqrt(DSOL);
370 
371  if (DSOL < DSOL1) break;
372  };
373 
374  if (verboseLevel > 3) {
375  if (itry == itry_max)
376  G4cout << " maximal number of iterations in potentialMinimization " << G4endl
377  << " A1 " << AF << " Z1 " << ZF << G4endl;
378 
379  };
380 
381  for (i = 0; i < 2; i++) ED[i] = F[i] * BET1[i] * BET1[i] + C[i] * AL1[i] * AL1[i];
382 
383  VC = D / R12;
384  VP = VC + ED[0] + ED[1];
385 }
const double C2
G4double randomGauss(G4double sigma)
const double C1
#define A13
double S(double temp)
G4double getZopt(G4int A1, G4int A2, G4int ZT, G4double X3, G4double X4, G4double R12) const
Definition: G4Fissioner.cc:196
G4double bindingEnergyAsymptotic(G4int A, G4int Z)
G4double getC2(G4int A1, G4int A2, G4double X3, G4double X4, G4double R12) const
Definition: G4Fissioner.cc:180
double C(double temp)
double B(double temperature)
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
G4double getNucleiMass() const
G4LorentzVector generateWithRandomAngles(G4double p, G4double mass=0.)
void potentialMinimization(G4double &VP, G4double(&ED)[2], G4double &VC, G4int AF, G4int AS, G4int ZF, G4int ZS, G4double AL1[2], G4double BET1[2], G4double &R12) const
Definition: G4Fissioner.cc:214
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
int G4lrint(double ad)
Definition: templates.hh:163
double D(double temp)
#define G4endl
Definition: G4ios.hh:61
void addRecoilFragment(const G4Fragment *aFragment)
double G4double
Definition: G4Types.hh:76
virtual void deExcite(const G4Fragment &target, G4CollisionOutput &output)
Definition: G4Fissioner.cc:66
G4double bindingEnergy(G4int A, G4int Z)
G4double nucleiLevelDensity(G4int A)
CLHEP::HepLorentzVector G4LorentzVector