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