Geant4_10
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 
49 #include "G4Fissioner.hh"
50 #include "G4CollisionOutput.hh"
51 #include "G4Fragment.hh"
52 #include "G4HadTmpUtil.hh"
53 #include "G4InuclNuclei.hh"
54 #include "G4InuclParticle.hh"
55 #include "G4FissionStore.hh"
58 
59 using namespace G4InuclSpecialFunctions;
60 
61 
63  G4CollisionOutput& output) {
64  if (verboseLevel) {
65  G4cout << " >>> G4Fissioner::deExcite" << G4endl;
66  }
67 
68  if (verboseLevel > 1)
69  G4cout << " Fissioner input\n" << target << G4endl;
70 
71  // Initialize buffer for fission possibilities
72  fissionStore.setVerboseLevel(verboseLevel);
73  fissionStore.clear();
74 
75  getTargetData(target);
76  G4double mass_in = PEX.m();
77  G4double e_in = mass_in; // Mass includes excitation
78  G4double PARA = 0.055 * G4cbrt(A*A) * (G4cbrt(A-Z) + G4cbrt(Z));
79  G4double TEM = std::sqrt(EEXS / PARA);
80  G4double TETA = 0.494 * G4cbrt(A) * TEM;
81 
82  TETA = TETA / std::sinh(TETA);
83 
84  if (A < 246) PARA += (nucleiLevelDensity(A) - PARA) * TETA;
85 
86  G4int A1 = A/2 + 1;
87  G4int Z1;
88  G4int A2 = A - A1;
89 
90  G4double ALMA = -1000.0;
91  G4double DM1 = bindingEnergy(A,Z);
92  G4double EVV = EEXS - DM1;
94  G4double DTEM = (A < 220 ? 0.5 : 1.15);
95 
96  TEM += DTEM;
97 
98  G4double AL1[2] = { -0.15, -0.15 };
99  G4double BET1[2] = { 0.05, 0.05 };
100 
101  G4double R12 = G4cbrt(A1) + G4cbrt(A2);
102 
103  for (G4int i = 0; i < 50 && A1 > 30; i++) {
104  A1--;
105  A2 = A - A1;
106  G4double X3 = 1.0 / G4cbrt(A1);
107  G4double X4 = 1.0 / G4cbrt(A2);
108  Z1 = G4lrint(getZopt(A1, A2, Z, X3, X4, R12) - 1.);
109  G4double EDEF1[2];
110  G4int Z2 = Z - Z1;
111  G4double VPOT, VCOUL;
112 
113  potentialMinimization(VPOT, EDEF1, VCOUL, A1, A2, Z1, Z2, AL1, BET1, R12);
114 
115  G4double DM3 = bindingEnergy(A1,Z1);
116  G4double DM4 = bindingEnergyAsymptotic(A1, Z1);
117  G4double DM5 = bindingEnergy(A2,Z2);
118  G4double DM6 = bindingEnergyAsymptotic(A2, Z2);
119  G4double DMT1 = DM4 + DM6 - DM2;
120  G4double DMT = DM3 + DM5 - DM1;
121  G4double EZL = EEXS + DMT - VPOT;
122 
123  if(EZL > 0.0) { // generate fluctuations
124  // faster, using randomGauss
125  G4double C1 = std::sqrt(getC2(A1, A2, X3, X4, R12) / TEM);
126  G4double DZ = randomGauss(C1);
127 
128  DZ = DZ > 0.0 ? DZ + 0.5 : -std::fabs(DZ - 0.5);
129  Z1 += G4int(DZ);
130  Z2 -= G4int(DZ);
131 
132  G4double DEfin = randomGauss(TEM);
133  G4double EZ = (DMT1 + (DMT - DMT1) * TETA - VPOT + DEfin) / TEM;
134 
135  if (EZ >= ALMA) ALMA = EZ;
136  G4double EK = VCOUL + DEfin + 0.5 * TEM;
137  G4double EV = EVV + bindingEnergy(A1,Z1) + bindingEnergy(A2,Z2) - EK;
138 
139  if (EV > 0.0) fissionStore.addConfig(A1, Z1, EZ, EK, EV);
140  };
141  };
142 
143  G4int store_size = fissionStore.size();
144  if (store_size == 0) return; // No fission products
145 
147  fissionStore.generateConfiguration(ALMA, inuclRndm());
148 
149  A1 = G4int(config.afirst);
150  A2 = A - A1;
151  Z1 = G4int(config.zfirst);
152 
153  G4int Z2 = Z - Z1;
154 
155  G4double mass1 = G4InuclNuclei::getNucleiMass(A1,Z1);
156  G4double mass2 = G4InuclNuclei::getNucleiMass(A2,Z2);
157  G4double EK = config.ekin;
158  G4double pmod = std::sqrt(0.001 * EK * mass1 * mass2 / mass_in);
159 
160  G4LorentzVector mom1 = generateWithRandomAngles(pmod, mass1);
161  G4LorentzVector mom2; mom2.setVectM(-mom1.vect(), mass2);
162 
163  G4double e_out = mom1.e() + mom2.e();
164  G4double EV = 1000.0 * (e_in - e_out) / A;
165  if (EV <= 0.0) return; // No fission energy
166 
167  G4double EEXS1 = EV*A1;
168  G4double EEXS2 = EV*A2;
169 
170  // Pass only last two nuclear fragments
171  output.addRecoilFragment(makeFragment(mom1, A1, Z1, EEXS1));
172  output.addRecoilFragment(makeFragment(mom2, A2, Z2, EEXS2));
173 }
174 
175 G4double G4Fissioner::getC2(G4int A1,
176  G4int A2,
177  G4double X3,
178  G4double X4,
179  G4double R12) const {
180 
181  if (verboseLevel > 3) {
182  G4cout << " >>> G4Fissioner::getC2" << G4endl;
183  }
184 
185  G4double C2 = 124.57 * (1.0 / A1 + 1.0 / A2) + 0.78 * (X3 + X4) - 176.9 *
186  ((X3*X3*X3*X3) + (X4*X4*X4*X4)) + 219.36 * (1.0 / (A1 * A1) + 1.0 / (A2 * A2)) - 1.108 / R12;
187 
188  return C2;
189 }
190 
191 G4double G4Fissioner::getZopt(G4int A1,
192  G4int A2,
193  G4int ZT,
194  G4double X3,
195  G4double X4,
196  G4double R12) const {
197 
198  if (verboseLevel > 3) {
199  G4cout << " >>> G4Fissioner::getZopt" << G4endl;
200  }
201 
202  G4double Zopt = (87.7 * (X4 - X3) * (1.0 - 1.25 * (X4 + X3)) +
203  ZT * ((124.57 / A2 + 0.78 * X4 - 176.9 * (X4*X4*X4*X4) + 219.36 / (A2 * A2)) - 0.554 / R12)) /
204  getC2(A1, A2, X3, X4, R12);
205 
206  return Zopt;
207 }
208 
209 void G4Fissioner::potentialMinimization(G4double& VP,
210  G4double( &ED)[2],
211  G4double& VC,
212  G4int AF,
213  G4int AS,
214  G4int ZF,
215  G4int ZS,
216  G4double AL1[2],
217  G4double BET1[2],
218  G4double& R12) const {
219 
220  if (verboseLevel > 3) {
221  G4cout << " >>> G4Fissioner::potentialMinimization" << G4endl;
222  }
223 
224  const G4double huge_num = 2.0e35;
225  const G4int itry_max = 2000;
226  const G4double DSOL1 = 1.0e-6;
227  const G4double DS1 = 0.3;
228  const G4double DS2 = 1.0 / DS1 / DS1;
229  G4int A1[2] = { AF, AS };
230  G4int Z1[2] = { ZF, ZS };
231  G4double D = 1.01844 * ZF * ZS;
232  G4double D0 = 1.0e-3 * D;
233  G4double R[2];
234  R12 = 0.0;
235  G4double C[2];
236  G4double F[2];
237  G4double Y1;
238  G4double Y2;
239  G4int i;
240 
241  for (i = 0; i < 2; i++) {
242  R[i] = G4cbrt(A1[i]);
243  Y1 = R[i] * R[i];
244  Y2 = Z1[i] * Z1[i] / R[i];
245  C[i] = 6.8 * Y1 - 0.142 * Y2;
246  F[i] = 12.138 * Y1 - 0.145 * Y2;
247  };
248 
249  G4double SAL[2];
250  G4double SBE[2];
251  G4double X[2];
252  G4double X1[2];
253  G4double X2[2];
254  G4double RAL[2];
255  G4double RBE[2];
256  G4double AA[4][4];
257  G4double B[4];
258  G4int itry = 0;
259 
260  while (itry < itry_max) {
261  itry++;
262  G4double S = 0.0;
263 
264  for (i = 0; i < 2; i++) {
265  S += R[i] * (1.0 + AL1[i] + BET1[i] - 0.257 * AL1[i] * BET1[i]);
266  };
267  R12 = 0.0;
268  Y1 = 0.0;
269  Y2 = 0.0;
270 
271  for (i = 0; i < 2; i++) {
272  SAL[i] = R[i] * (1.0-0.257 * BET1[i]);
273  SBE[i] = R[i] * (1.0-0.257 * AL1[i]);
274  X[i] = R[i] / S;
275  X1[i] = X[i] * X[i];
276  X2[i] = X[i] * X1[i];
277  Y1 += AL1[i] * X1[i];
278  Y2 += BET1[i] * X2[i];
279  R12 += R[i] * (1.0 - AL1[i] * (1.0 - 0.6 * X[i]) + BET1[i] * (1.0 - 0.429 * X1[i]));
280  };
281 
282  G4double Y3 = -0.6 * Y1 + 0.857 * Y2;
283  G4double Y4 = (1.2 * Y1 - 2.571 * Y2) / S;
284  G4double R2 = D0 / (R12 * R12);
285  G4double R3 = 2.0 * R2 / R12;
286 
287  for (i = 0; i < 2; i++) {
288  RAL[i] = -R[i] * (1.0 - 0.6 * X[i]) + SAL[i] * Y3;
289  RBE[i] = R[i] * (1.0 - 0.429 * X1[i]) + SBE[i] * Y3;
290  };
291 
292  G4double DX1;
293  G4double DX2;
294 
295  for (i = 0; i < 2; i++) {
296 
297  for (G4int j = 0; j < 2; j++) {
298  G4double DEL1 = i == j ? 1.0 : 0.0;
299  DX1 = 0.0;
300  DX2 = 0.0;
301 
302  if (std::fabs(AL1[i]) >= DS1) {
303  G4double XXX = AL1[i] * AL1[i] * DS2;
304  G4double DEX = XXX > 100.0 ? huge_num : std::exp(XXX);
305  DX1 = 2.0 * (1.0 + 2.0 * AL1[i] * AL1[i] * DS2) * DEX * DS2;
306  };
307 
308  if (std::fabs(BET1[i]) >= DS1) {
309  G4double XXX = BET1[i] * BET1[i] * DS2;
310  G4double DEX = XXX > 100.0 ? huge_num : std::exp(XXX);
311  DX2 = 2.0 * (1.+2.0 * BET1[i] * BET1[i] * DS2) * DEX * DS2;
312  };
313 
314  G4double DEL = 2.0e-3 * DEL1;
315  AA[i][j] = R3 * RBE[i] * RBE[j] -
316  R2 * (-0.6 *
317  (X1[i] * SAL[j] +
318  X1[j] * SAL[i]) + SAL[i] * SAL[j] * Y4) +
319  DEL * C[i] + DEL1 * DX1;
320  G4int i1 = i + 2;
321  G4int j1 = j + 2;
322  AA[i1][j1] = R3 * RBE[i] * RBE[j]
323  - R2 * (0.857 *
324  (X2[i] * SBE[j] +
325  X2[j] * SBE[i]) + SBE[i] * SBE[j] * Y4) +
326  DEL * F[i] + DEL1 * DX2;
327  AA[i][j1] = R3 * RAL[i] * RBE[j] -
328  R2 * (0.857 *
329  (X2[j] * SAL[i] -
330  0.6 * X1[i] * SBE[j]) + SBE[j] * SAL[i] * Y4 -
331  0.257 * R[i] * Y3 * DEL1);
332  AA[j1][i] = AA[i][j1];
333  };
334  };
335 
336  for (i = 0; i < 2; i++) {
337  DX1 = 0.0;
338  DX2 = 0.0;
339 
340  if (std::fabs(AL1[i]) >= DS1) DX1 = 2.0 * AL1[i] * DS2 * std::exp(AL1[i] * AL1[i] * DS2);
341 
342  if (std::fabs(BET1[i]) >= DS1) DX2 = 2.0 * BET1[i] * DS2 * std::exp(BET1[i] * BET1[i] * DS2);
343  B[i] = R2 * RAL[i] - 2.0e-3 * C[i] * AL1[i] + DX1;
344  B[i + 2] = R2 * RBE[i] - 2.0e-3 * F[i] * BET1[i] + DX2;
345  };
346 
347  G4double ST = 0.0;
348  G4double ST1 = 0.0;
349 
350  for (i = 0; i < 4; i++) {
351  ST += B[i] * B[i];
352 
353  for (G4int j = 0; j < 4; j++) ST1 += AA[i][j] * B[i] * B[j];
354  };
355 
356  G4double STEP = ST / ST1;
357  G4double DSOL = 0.0;
358 
359  for (i = 0; i < 2; i++) {
360  AL1[i] += B[i] * STEP;
361  BET1[i] += B[i + 2] * STEP;
362  DSOL += B[i] * B[i] + B[i + 2] * B[i + 2];
363  };
364  DSOL = std::sqrt(DSOL);
365 
366  if (DSOL < DSOL1) break;
367  };
368 
369  if (verboseLevel > 3) {
370  if (itry == itry_max)
371  G4cout << " maximal number of iterations in potentialMinimization " << G4endl
372  << " A1 " << AF << " Z1 " << ZF << G4endl;
373 
374  };
375 
376  for (i = 0; i < 2; i++) ED[i] = F[i] * BET1[i] * BET1[i] + C[i] * AL1[i] * AL1[i];
377 
378  VC = D / R12;
379  VP = VC + ED[0] + ED[1];
380 }
Double_t Z2
Definition: plot.C:266
G4double randomGauss(G4double sigma)
Double_t X1
Definition: plot.C:266
G4double bindingEnergyAsymptotic(G4int A, G4int Z)
Double_t Y2
Definition: plot.C:266
Double_t X2
Definition: plot.C:266
const XML_Char * target
Definition: expat.h:268
int G4int
Definition: G4Types.hh:78
Float_t X
Definition: plot.C:39
void setVectM(const Hep3Vector &spatial, double mass)
Hep3Vector vect() const
G4GLOB_DLL std::ostream G4cout
Float_t Z
Definition: plot.C:39
G4double getNucleiMass() const
G4LorentzVector generateWithRandomAngles(G4double p, G4double mass=0.)
Double_t Y1
Definition: plot.C:266
#define C1
int G4lrint(double ad)
Definition: templates.hh:163
#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:62
G4double bindingEnergy(G4int A, G4int Z)
G4double nucleiLevelDensity(G4int A)
Double_t Z1
Definition: plot.C:266