ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4EquilibriumEvaporator.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4EquilibriumEvaporator.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 // 20100308 M. Kelsey -- Bug fix for setting masses of evaporating nuclei
29 // 20100319 M. Kelsey -- Use new generateWithRandomAngles for theta,phi stuff;
30 // eliminate some unnecessary std::pow()
31 // 20100319 M. Kelsey -- Bug fix in new GetBindingEnergy() use right after
32 // goodRemnant() -- Q1 should be outside call.
33 // 20100412 M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
34 // 20100413 M. Kelsey -- Pass buffers to paraMaker[Truncated]
35 // 20100419 M. Kelsey -- Handle fission output list via const-ref
36 // 20100517 M. Kelsey -- Use G4CascadeInterpolator for QFF
37 // 20100520 M. Kelsey -- Inherit from common base class, make other colliders
38 // simple data members. Rename timeToBigBang() to override
39 // base explosion().
40 // 20100617 M. Kelsey -- Remove "RUN" preprocessor flag and all "#else" code,
41 // pass verbosity to colliders.
42 // 20100620 M. Kelsey -- Use local "bindingEnergy()" function to call through.
43 // 20100701 M. Kelsey -- Don't need to add excitation to nuclear mass; compute
44 // new excitation energies properly (mass differences)
45 // 20100702 M. Kelsey -- Simplify if-cascades, indentation
46 // 20100712 M. Kelsey -- Add conservation checking
47 // 20100714 M. Kelsey -- Move conservation checking to base class. Use
48 // _generated_ evaporate energy (S) to adjust EEXS directly,
49 // and test for S < EEXS before any particle generation; shift
50 // nucleus momentum (PEX) by evaporate momentum directly
51 // 20100719 M. Kelsey -- Remove duplicative EESX_new calculation.
52 // 20100923 M. Kelsey -- Migrate to integer A and Z
53 // 20110214 M. Kelsey -- Follow G4InuclParticle::Model enumerator migration
54 // 20110728 M. Kelsey -- Fix Coverity #22951: check for "icase = -1" after
55 // loop which is supposed to set it. Otherwise indexing is wrong.
56 // 20110801 M. Kelsey -- Move "parms" buffer to data member, allocate in
57 // constructor.
58 // 20110809 M. Kelsey -- Move "foutput" to data member, get list by reference;
59 // create final-state particles within "push_back" to avoid
60 // creation of temporaries.
61 // 20110922 M. Kelsey -- Follow G4InuclParticle::print(ostream&) migration
62 // 20111007 M. Kelsey -- Add G4InuclParticleNames, replace hardcoded numbers
63 // 20120608 M. Kelsey -- Fix variable-name "shadowing" compiler warnings.
64 // 20121009 M. Kelsey -- Improve debugging output
65 // 20130129 M. Kelsey -- Move QF interpolation to global statics for
66 // multi-threaded shared memory.
67 // 20130621 M. Kelsey -- Remove turning off conservation checks after fission
68 // 20130622 Inherit from G4CascadeDeexciteBase, move to deExcite() interface
69 // with G4Fragment
70 // 20130628 Fissioner produces G4Fragment output directly in G4CollisionOutput
71 // 20130808 M. Kelsey -- Use new object-version of paraMaker, for thread safety
72 // 20130924 M. Kelsey -- Use G4Log, G4Exp for CPU speedup
73 // 20150608 M. Kelsey -- Label all while loops as terminating.
74 // 20150622 M. Kelsey -- For new G4cbrt(int), move powers of A outside.
75 
76 #include "G4EquilibriumEvaporator.hh"
77 #include "G4SystemOfUnits.hh"
78 #include "Randomize.hh"
79 #include "G4BigBanger.hh"
80 #include "G4CascadeInterpolator.hh"
81 #include "G4CollisionOutput.hh"
82 #include "G4Fissioner.hh"
83 #include "G4Fragment.hh"
84 #include "G4InuclNuclei.hh"
85 #include "G4InuclSpecialFunctions.hh"
86 #include "G4InuclParticleNames.hh"
87 #include "G4LorentzConvertor.hh"
88 #include "G4LorentzVector.hh"
89 #include "G4ThreeVector.hh"
90 #include "G4Log.hh"
91 #include "G4Exp.hh"
92 
93 using namespace G4InuclParticleNames;
94 using namespace G4InuclSpecialFunctions;
95 
96 namespace { // Interpolation arrays for QF
97  const G4double QFREP[72] = {
98  // TL201 * * * *
99  // 1 2 3 4 5
100  22.5, 22.0, 21.0, 21.0, 20.0,
101  // BI209 BI207 PO210 AT213 * TH234
102  // 6 7 8 9 10 11
103  20.6, 20.6, 18.6, 15.8, 13.5, 6.5,
104  // TH233 TH232 TH231 TH230 TX229 PA233 PA232 PA231 PA230 U240
105  // 12 13 14 15 16 17 18 19 20 21
106  6.65, 6.22, 6.27, 6.5, 6.7, 6.2, 6.25, 5.9, 6.1, 5.75,
107  // U239 U238 U237 U236 U235 U234 U233 U232 U231
108  // 22 23 24 25 26 27 28 29 30
109  6.46, 5.7, 6.28, 5.8, 6.15, 5.6, 5.8, 5.2, 5.8,
110  // NP238 NP237 NP236 NP235 PU245 NP234 PU244 NP233
111  // 31 32 33 34 35 36 37 38
112  6.2 , 5.9 , 5.9, 6.0, 5.8, 5.7, 5.4, 5.4,
113  // PU242 PU241 PU240 PU239 PU238 AM247 PU236 AM245 AM244 AM243
114  // 39 40 41 42 43 44 45 46 47 48
115  5.6, 6.1, 5.57, 6.3, 5.5, 5.8, 4.7, 6.2, 6.4, 6.2,
116  // AM242 AM241 AM240 CM250 AM239 CM249 CM248 CM247 CM246
117  // 49 50 51 52 53 54 55 56 57
118  6.5, 6.2, 6.5, 5.3, 6.4, 5.7, 5.7, 6.2, 5.7,
119  // CM245 CM244 CM243 CM242 CM241 BK250 CM240
120  // 58 59 60 61 62 63 64
121  6.3, 5.8, 6.7, 5.8, 6.6, 6.1, 4.3,
122  // BK249 CF252 CF250 CF248 CF246 ES254 ES253 FM254
123  // 65 66 67 68 69 70 71 72
124  6.2, 3.8, 5.6, 4.0, 4.0, 4.2, 4.2, 3.5 };
125 
126  static const G4double XREP[72] = {
127  // 1 2 3 4 5
128  0.6761, 0.677, 0.6788, 0.6803, 0.685,
129  // 6 7 8 9 10 11
130  0.6889, 0.6914, 0.6991, 0.7068, 0.725, 0.7391,
131  // 12 13 14 15 16 17 18 19 20 21
132  0.74, 0.741, 0.742, 0.743, 0.744, 0.7509, 0.752, 0.7531, 0.7543, 0.7548,
133  // 22 23 24
134  0.7557, 0.7566, 0.7576,
135  // 25 26 27 28 29 30 31 32 33 34
136  0.7587, 0.7597, 0.7608, 0.762, 0.7632, 0.7644, 0.7675, 0.7686, 0.7697, 0.7709,
137  // 35 36 37 38 39 40 41
138  0.7714, 0.7721, 0.7723, 0.7733, 0.7743, 0.7753, 0.7764,
139  // 42 43 44 45 46 47 48 49
140  0.7775, 0.7786, 0.7801, 0.781, 0.7821, 0.7831, 0.7842, 0.7852,
141  // 50 51 52 53 54 55 56 57 58
142  0.7864, 0.7875, 0.7880, 0.7887, 0.7889, 0.7899, 0.7909, 0.7919, 0.7930,
143  // 59 60 61 62 63 64
144  0.7941, 0.7953, 0.7965, 0.7977, 0.7987, 0.7989,
145  // 65 66 67 68 69 70 71 72
146  0.7997, 0.8075, 0.8097, 0.8119, 0.8143, 0.8164, 0.8174, 0.8274 };
147 }
148 
149 
150 // Constructor and destructor
151 
152 G4EquilibriumEvaporator::G4EquilibriumEvaporator()
153  : G4CascadeDeexciteBase("G4EquilibriumEvaporator"),
154  theParaMaker(verboseLevel), QFinterp(XREP) {
155  parms.first.resize(6,0.);
156  parms.second.resize(6,0.);
157 }
158 
159 G4EquilibriumEvaporator::~G4EquilibriumEvaporator() {}
160 
161 void G4EquilibriumEvaporator::setVerboseLevel(G4int verbose) {
162  G4CascadeDeexciteBase::setVerboseLevel(verbose);
163  theFissioner.setVerboseLevel(verbose);
164  theBigBanger.setVerboseLevel(verbose);
165 }
166 
167 
168 // Main operation
169 
170 void G4EquilibriumEvaporator::deExcite(const G4Fragment& target,
171  G4CollisionOutput& output) {
172  if (verboseLevel) {
173  G4cout << " >>> G4EquilibriumEvaporator::deExcite" << G4endl;
174  }
175 
176  if (verboseLevel>1) G4cout << " evaporating target: \n" << target << G4endl;
177 
178  // Implementation of the equilibium evaporation according to Dostrovsky
179  // (Phys. Rev. 116 (1959) 683.
180  // Note that pairing energy which is of order 1 MeV for odd-even and even-odd
181  // nuclei is not included
182 
183  // Element in array: 0 : neutron
184  // 1 : proton
185  // 2 : deuteron
186  // 3 : triton
187  // 4 : 3He
188  // 5 : alpha
189 
190  const G4double huge_num = 50.0;
191  const G4double small = -50.0;
192  const G4double prob_cut_off = 1.0e-15;
193  const G4double Q1[6] = { 0.0, 0.0, 2.23, 8.49, 7.72, 28.3 }; // binding energy
194  const G4int AN[6] = { 1, 1, 2, 3, 3, 4 }; // mass number
195  const G4int Q[6] = { 0, 1, 1, 1, 2, 2 }; // charge
196  const G4double G[6] = { 2.0, 2.0, 6.0, 6.0, 6.0, 4.0 };
197  const G4double BE = 0.0063;
198  const G4double fisssion_cut = 1000.0;
199  const G4double cut_off_energy = 0.1;
200 
201  const G4double BF = 0.0242;
202  const G4int itry_max = 1000;
203  const G4int itry_global_max = 1000;
204  const G4double small_ekin = 1.0e-6;
205  const G4int itry_gam_max = 100;
206 
207  G4double W[8];
208  G4double u[6]; // Level density for each particle type: (Atarg - Aemitted)*parlev
209  G4double V[6]; // Coulomb energy
210  G4double TM[6]; // Maximum possible kinetic energy of emitted particle
211  G4int A1[6]; // A of remnant nucleus
212  G4int Z1[6]; // Z of remnant nucleus
213 
214  getTargetData(target);
215  if (verboseLevel > 3) G4cout << " after noeq: eexs " << EEXS << G4endl;
216 
217  G4InuclElementaryParticle dummy(small_ekin, proton);
218  G4LorentzConvertor toTheNucleiSystemRestFrame;
219  //*** toTheNucleiSystemRestFrame.setVerbose(verboseLevel);
220  toTheNucleiSystemRestFrame.setBullet(dummy);
221 
222  G4LorentzVector ppout;
223 
224  // See if fragment should just be dispersed
225  if (explosion(A, Z, EEXS)) {
226  if (verboseLevel > 1) G4cout << " big bang in eql start " << G4endl;
227  theBigBanger.deExcite(target, output);
228 
229  validateOutput(target, output); // Check energy conservation
230  return;
231  }
232 
233  // If nucleus is in ground state, no evaporation
234  if (EEXS < cut_off_energy) {
235  if (verboseLevel > 1) G4cout << " no energy for evaporation" << G4endl;
236  output.addRecoilFragment(target);
237 
238  validateOutput(target, output); // Check energy conservation
239  return;
240  }
241 
242  // Initialize evaporation attempts
243  G4double coul_coeff = (A >= 100.0) ? 1.4 : 1.2;
244 
245  G4LorentzVector pin = PEX; // Save original target for testing
246 
247  G4bool try_again = true;
248  G4bool fission_open = true;
249  G4int itry_global = 0;
250 
251  /* Loop checking 08.06.2015 MHK */
252  while (try_again && itry_global < itry_global_max) {
253  itry_global++;
254 
255  // Set rest frame of current (recoiling) nucleus
256  toTheNucleiSystemRestFrame.setTarget(PEX);
257  toTheNucleiSystemRestFrame.toTheTargetRestFrame();
258 
259  if (verboseLevel > 2) {
260  G4double nuc_mass = G4InuclNuclei::getNucleiMass(A, Z, EEXS);
261  G4cout << " A " << A << " Z " << Z << " mass " << nuc_mass
262  << " EEXS " << EEXS << G4endl;
263  }
264 
265  if (explosion(A, Z, EEXS)) { // big bang
266  if (verboseLevel > 2)
267  G4cout << " big bang in eql step " << itry_global << G4endl;
268 
269  theBigBanger.deExcite(makeFragment(PEX,A,Z,EEXS), output);
270 
271  validateOutput(target, output); // Check energy conservation
272  return;
273  }
274 
275  if (EEXS < cut_off_energy) { // Evaporation not possible
276  if (verboseLevel > 2)
277  G4cout << " no energy for evaporation in eql step " << itry_global
278  << G4endl;
279 
280  try_again = false;
281  break;
282  }
283 
284  // Normal evaporation chain
285  G4double E0 = getE0(A);
286  G4double parlev = getPARLEVDEN(A, Z);
287  G4double u1 = parlev * A;
288 
289  theParaMaker.getParams(Z, parms);
290  const std::vector<G4double>& AK = parms.first;
291  const std::vector<G4double>& CPA = parms.second;
292 
293  G4double DM0 = bindingEnergy(A,Z);
294  G4int i(0);
295 
296  for (i = 0; i < 6; i++) {
297  A1[i] = A - AN[i];
298  Z1[i] = Z - Q[i];
299  u[i] = parlev * A1[i];
300  V[i] = 0.0;
301  TM[i] = -0.1;
302 
303  if (goodRemnant(A1[i], Z1[i])) {
304  G4double QB = DM0 - bindingEnergy(A1[i],Z1[i]) - Q1[i];
305  V[i] = coul_coeff * Z * Q[i] * AK[i] / (1.0 + EEXS / E0) /
306  (G4cbrt(A1[i]) + G4cbrt(AN[i]));
307  TM[i] = EEXS - QB - V[i] * A / A1[i];
308  if (verboseLevel > 3) {
309  G4cout << " i=" << i << " : QB " << QB << " u " << u[i]
310  << " V " << V[i] << " TM " << TM[i] << G4endl;
311  }
312  }
313  }
314 
315  G4double ue = 2.0 * std::sqrt(u1 * EEXS);
316  G4double prob_sum = 0.0;
317 
318  W[0] = 0.0;
319  if (TM[0] > cut_off_energy) {
320  G4double AL = getAL(A);
321  G4double A13 = G4cbrt(A1[0]);
322  W[0] = BE * A13*A13 * G[0] * AL;
323  G4double TM1 = 2.0 * std::sqrt(u[0] * TM[0]) - ue;
324 
325  if (TM1 > huge_num) TM1 = huge_num;
326  else if (TM1 < small) TM1 = small;
327 
328  W[0] *= G4Exp(TM1);
329  prob_sum += W[0];
330  }
331 
332  for (i = 1; i < 6; i++) {
333  W[i] = 0.0;
334  if (TM[i] > cut_off_energy) {
335  G4double A13 = G4cbrt(A1[i]);
336  W[i] = BE * A13*A13 * G[i] * (1.0 + CPA[i]);
337  G4double TM1 = 2.0 * std::sqrt(u[i] * TM[i]) - ue;
338 
339  if (TM1 > huge_num) TM1 = huge_num;
340  else if (TM1 < small) TM1 = small;
341 
342  W[i] *= G4Exp(TM1);
343  prob_sum += W[i];
344  }
345  }
346 
347  // fisson part
348  W[6] = 0.0;
349  if (A >= 100.0 && fission_open) {
350  G4double X2 = Z * Z / A;
351  G4double X1 = 1.0 - 2.0 * Z / A;
352  G4double X = 0.019316 * X2 / (1.0 - 1.79 * X1 * X1);
353  G4double EF = EEXS - getQF(X, X2, A, Z, EEXS);
354 
355  if (EF > 0.0) {
356  G4double AF = u1 * getAF(X, A, Z, EEXS);
357  G4double TM1 = 2.0 * std::sqrt(AF * EF) - ue;
358 
359  if (TM1 > huge_num) TM1 = huge_num;
360  else if (TM1 < small) TM1 = small;
361 
362  W[6] = BF * G4Exp(TM1);
363  if (W[6] > fisssion_cut*W[0]) W[6] = fisssion_cut*W[0];
364 
365  prob_sum += W[6];
366  }
367  }
368 
369  // again time to decide what next
370  if (verboseLevel > 2) {
371  G4cout << " Evaporation probabilities: sum = " << prob_sum
372  << "\n\t n " << W[0] << " p " << W[1] << " d " << W[2]
373  << " He3 " << W[3] << "\n\t t " << W[4] << " alpha " << W[5]
374  << " fission " << W[6] << G4endl;
375  }
376 
377  G4int icase = -1;
378 
379  if (prob_sum < prob_cut_off) { // photon emission chain
380  G4double UCR0 = 2.5 + 150.0 / A;
381  G4double T00 = 1.0 / (std::sqrt(u1 / UCR0) - 1.25 / UCR0);
382 
383  if (verboseLevel > 3)
384  G4cout << " UCR0 " << UCR0 << " T00 " << T00 << G4endl;
385 
386  G4int itry_gam = 0;
387  while (EEXS > cut_off_energy && try_again) {
388  itry_gam++;
389  G4int itry = 0;
390  G4double T04 = 4.0 * T00;
391  G4double FMAX;
392 
393  if (T04 < EEXS) {
394  FMAX = (T04*T04*T04*T04) * G4Exp((EEXS - T04) / T00);
395  } else {
396  FMAX = EEXS*EEXS*EEXS*EEXS;
397  };
398 
399  if (verboseLevel > 3)
400  G4cout << " T04 " << T04 << " FMAX (EEXS^4) " << FMAX << G4endl;
401 
402  G4double S(0), X1(0);
403  while (itry < itry_max) {
404  itry++;
405  S = EEXS * inuclRndm();
406  X1 = (S*S*S*S) * G4Exp((EEXS - S) / T00);
407 
408  if (X1 > FMAX * inuclRndm()) break;
409  };
410 
411  if (itry == itry_max) { // Maximum attempts exceeded
412  try_again = false;
413  break;
414  }
415 
416  if (verboseLevel > 2) G4cout << " photon escape ?" << G4endl;
417 
418  if (S < EEXS) { // Valid evaporate
419  S /= GeV; // Convert to Bertini units
420 
421  G4double pmod = S;
422  G4LorentzVector mom = generateWithRandomAngles(pmod, 0.);
423 
424  // Push evaporated particle into current rest frame
425  mom = toTheNucleiSystemRestFrame.backToTheLab(mom);
426 
427  if (verboseLevel > 3) {
428  G4cout << " nucleus px " << PEX.px() << " py " << PEX.py()
429  << " pz " << PEX.pz() << " E " << PEX.e() << G4endl
430  << " evaporate px " << mom.px() << " py " << mom.py()
431  << " pz " << mom.pz() << " E " << mom.e() << G4endl;
432  }
433 
434  PEX -= mom; // Remaining four-momentum
435  EEXS -= S*GeV; // New excitation energy (in MeV)
436 
437  // NOTE: In-situ construction will be optimized away (no copying)
438  output.addOutgoingParticle(G4InuclElementaryParticle(mom, photon,
439  G4InuclParticle::Equilib));
440 
441  if (verboseLevel > 3)
442  G4cout << output.getOutgoingParticles().back() << G4endl;
443 
444  ppout += mom;
445  } else {
446  if (itry_gam == itry_gam_max) try_again = false;
447  }
448  } // while (EEXS > cut_off
449  try_again = false;
450  } else { // if (prob_sum < prob_cut_off)
451  G4double SL = prob_sum * inuclRndm();
452  if (verboseLevel > 3) G4cout << " random SL " << SL << G4endl;
453 
454  G4double S1 = 0.0;
455  for (i = 0; i < 7; i++) { // Select evaporation scenario
456  S1 += W[i];
457  if (SL <= S1) {
458  icase = i;
459  break;
460  };
461  };
462 
463  if (icase < 0) continue; // Failed to choose scenario, try again
464 
465  if (icase < 6) { // particle or light nuclei escape
466  if (verboseLevel > 2) {
467  G4cout << " particle/light-ion escape ("
468  << (icase==0 ? "neutron" : icase==1 ? "proton" :
469  icase==2 ? "deuteron" : icase==3 ? "He3" :
470  icase==4 ? "triton" : icase==5 ? "alpha" : "ERROR!")
471  << ")?" << G4endl;
472  }
473 
474  G4double Xmax = (std::sqrt(u[icase]*TM[icase] + 0.25) - 0.5)/u[icase];
475  G4int itry1 = 0;
476  G4bool bad = true;
477 
478  while (itry1 < itry_max && bad) {
479  itry1++;
480  G4int itry = 0;
481  G4double S = 0.0;
482  G4double X = 0.0;
483  G4double Ptest = 0.0;
484 
485  while (itry < itry_max) {
486  itry++;
487 
488  // Sampling from eq. 17 of Dostrovsky
489  X = G4UniformRand()*TM[icase];
490  Ptest = (X/Xmax)*G4Exp(-2.*u[icase]*Xmax + 2.*std::sqrt(u[icase]*(TM[icase] - X)));
491  if (G4UniformRand() < Ptest) {
492  S = X + V[icase];
493  break;
494  }
495  }
496 
497  if (S > V[icase] && S < EEXS) { // real escape
498 
499  if (verboseLevel > 2)
500  G4cout << " escape itry1 " << itry1 << " icase "
501  << icase << " S (MeV) " << S << G4endl;
502 
503  S /= GeV; // Convert to Bertini units
504 
505  if (icase < 2) { // particle escape
506  G4int ptype = 2 - icase;
507  if (verboseLevel > 2)
508  G4cout << " particle " << ptype << " escape" << G4endl;
509 
510  // generate particle momentum
511  G4double mass =
512  G4InuclElementaryParticle::getParticleMass(ptype);
513 
514  G4double pmod = std::sqrt((2.0 * mass + S) * S);
515  G4LorentzVector mom = generateWithRandomAngles(pmod, mass);
516 
517  // Push evaporated particle into current rest frame
518  mom = toTheNucleiSystemRestFrame.backToTheLab(mom);
519 
520  if (verboseLevel > 2) {
521  G4cout << " nucleus px " << PEX.px() << " py " << PEX.py()
522  << " pz " << PEX.pz() << " E " << PEX.e() << G4endl
523  << " evaporate px " << mom.px() << " py " << mom.py()
524  << " pz " << mom.pz() << " E " << mom.e() << G4endl;
525  }
526 
527  // New excitation energy depends on residual nuclear state
528  G4double mass_new =
529  G4InuclNuclei::getNucleiMass(A1[icase],Z1[icase]);
530 
531  G4double EEXS_new = ((PEX-mom).m() - mass_new)*GeV;
532  if (EEXS_new < 0.0) continue; // Sanity check for new nucleus
533 
534  PEX -= mom; // Remaining four-momentum
535  EEXS = EEXS_new;
536 
537  A = A1[icase];
538  Z = Z1[icase];
539 
540  // NOTE: In-situ construction optimizes away (no copying)
541  output.addOutgoingParticle(G4InuclElementaryParticle(mom,
542  ptype, G4InuclParticle::Equilib));
543 
544  if (verboseLevel > 3)
545  G4cout << output.getOutgoingParticles().back() << G4endl;
546 
547  ppout += mom;
548  bad = false;
549  } else { // if (icase < 2)
550  if (verboseLevel > 2) {
551  G4cout << " nucleus A " << AN[icase] << " Z " << Q[icase]
552  << " escape icase " << icase << G4endl;
553  }
554 
555  G4double mass =
556  G4InuclNuclei::getNucleiMass(AN[icase],Q[icase]);
557 
558  // generate particle momentum
559  G4double pmod = std::sqrt((2.0 * mass + S) * S);
560  G4LorentzVector mom = generateWithRandomAngles(pmod,mass);
561 
562  // Push evaporated particle into current rest frame
563  mom = toTheNucleiSystemRestFrame.backToTheLab(mom);
564 
565  if (verboseLevel > 2) {
566  G4cout << " nucleus px " << PEX.px() << " py " << PEX.py()
567  << " pz " << PEX.pz() << " E " << PEX.e() << G4endl
568  << " evaporate px " << mom.px() << " py " << mom.py()
569  << " pz " << mom.pz() << " E " << mom.e() << G4endl;
570  }
571 
572  // New excitation energy depends on residual nuclear state
573  G4double mass_new =
574  G4InuclNuclei::getNucleiMass(A1[icase],Z1[icase]);
575 
576  G4double EEXS_new = ((PEX-mom).m() - mass_new)*GeV;
577  if (EEXS_new < 0.0) continue; // Sanity check for new nucleus
578 
579  PEX -= mom; // Remaining four-momentum
580  EEXS = EEXS_new;
581 
582  A = A1[icase];
583  Z = Z1[icase];
584 
585  // NOTE: In-situ constructor optimizes away (no copying)
586  output.addOutgoingNucleus(G4InuclNuclei(mom,
587  AN[icase], Q[icase], 0.*GeV,
588  G4InuclParticle::Equilib));
589 
590  if (verboseLevel > 3)
591  G4cout << output.getOutgoingNuclei().back() << G4endl;
592 
593  ppout += mom;
594  bad = false;
595  } // if (icase < 2)
596  } // if (EPR > 0.0 ...
597  } // while (itry1 ...
598 
599  if (itry1 == itry_max || bad) try_again = false;
600  } else { // if (icase < 6)
601  if (verboseLevel > 2) {
602  G4cout << " fission: A " << A << " Z " << Z << " eexs " << EEXS
603  << " Wn " << W[0] << " Wf " << W[6] << G4endl;
604  }
605 
606  // Catch fission output separately for verification
607  fission_output.reset();
608  theFissioner.deExcite(makeFragment(A,Z,EEXS), fission_output);
609 
610  if (fission_output.numberOfFragments() == 2) { // fission ok
611  if (verboseLevel > 2) G4cout << " fission done in eql" << G4endl;
612 
613  // Move fission fragments to lab frame for processing
614  fission_output.boostToLabFrame(toTheNucleiSystemRestFrame);
615 
616  // Now evaporate the fission fragments individually
617  this->deExcite(fission_output.getRecoilFragment(0), output);
618  this->deExcite(fission_output.getRecoilFragment(1), output);
619 
620  validateOutput(target, output); // Check energy conservation
621  return;
622  } else { // fission forbidden now
623  fission_open = false;
624  }
625  } // End of fission case
626  } // if (prob_sum < prob_cut_off)
627  } // while (try_again
628 
629  // this time it's final nuclei
630 
631  if (itry_global == itry_global_max) {
632  if (verboseLevel > 3) {
633  G4cout << " ! itry_global " << itry_global_max << G4endl;
634  }
635  }
636 
637  G4LorentzVector pnuc = pin - ppout;
638 
639  // NOTE: In-situ constructor will be optimized away (no copying)
640  output.addOutgoingNucleus(G4InuclNuclei(pnuc, A, Z, 0.,
641  G4InuclParticle::Equilib));
642 
643  if (verboseLevel > 3) {
644  G4cout << " remaining nucleus \n" << output.getOutgoingNuclei().back()
645  << G4endl;
646  }
647 
648 
649  validateOutput(target, output); // Check energy conservation
650  return;
651 }
652 
653 G4bool G4EquilibriumEvaporator::explosion(G4int a,
654  G4int z,
655  G4double e) const {
656  if (verboseLevel > 3) {
657  G4cout << " >>> G4EquilibriumEvaporator::explosion? ";
658  }
659 
660  const G4double be_cut = 3.0;
661 
662  // Different criteria from base class, since nucleus more "agitated"
663  G4bool bigb = (!(a >= 12 && z >= 0 && z < 3*(a-z)) &&
664  (e >= be_cut * bindingEnergy(a,z))
665  );
666 
667  if (verboseLevel > 3) G4cout << bigb << G4endl;
668 
669  return bigb;
670 }
671 
672 G4bool G4EquilibriumEvaporator::goodRemnant(G4int a,
673  G4int z) const {
674  if (verboseLevel > 3) {
675  G4cout << " >>> G4EquilibriumEvaporator::goodRemnant(" << a << "," << z
676  << ")? " << (a>1 && z>0 && a>z) << G4endl;
677  }
678 
679  return a > 1 && z > 0 && a > z;
680 }
681 
682 G4double G4EquilibriumEvaporator::getQF(G4double x,
683  G4double x2,
684  G4int a,
685  G4int /*z*/,
686  G4double ) const {
687  if (verboseLevel > 3) {
688  G4cout << " >>> G4EquilibriumEvaporator::getQF ";
689  }
690 
691  const G4double G0 = 20.4;
692  const G4double XMIN = 0.6761;
693  const G4double XMAX = 0.8274;
694 
695  G4double QFF = 0.0;
696 
697  if (x < XMIN || x > XMAX) {
698  G4double X1 = 1.0 - 0.02 * x2;
699  G4double FX = (0.73 + (3.33 * X1 - 0.66) * X1) * (X1*X1*X1);
700  G4double A13 = G4cbrt(a);
701  QFF = G0 * FX * A13*A13;
702  } else {
703  QFF = QFinterp.interpolate(x, QFREP);
704  }
705 
706  if (QFF < 0.0) QFF = 0.0;
707 
708  if (verboseLevel>3) G4cout << " returns " << QFF << G4endl;
709 
710  return QFF;
711 }
712 
713 G4double G4EquilibriumEvaporator::getAF(G4double ,
714  G4int /*a*/,
715  G4int /*z*/,
716  G4double e) const {
717 
718  if (verboseLevel > 3) {
719  G4cout << " >>> G4EquilibriumEvaporator::getAF" << G4endl;
720  }
721 
722  // ugly parameterisation to fit the experimental fission cs for Hg - Bi nuclei
723  G4double AF = 1.285 * (1.0 - e / 1100.0);
724 
725  if(AF < 1.06) AF = 1.06;
726  // if(AF < 1.08) AF = 1.08;
727 
728  return AF;
729 }
730 
731 G4double G4EquilibriumEvaporator::getPARLEVDEN(G4int /*a*/,
732  G4int /*z*/) const {
733 
734  if (verboseLevel > 3) {
735  G4cout << " >>> G4EquilibriumEvaporator::getPARLEVDEN" << G4endl;
736  }
737 
738  const G4double par = 0.125;
739 
740  return par;
741 }
742 
743 G4double G4EquilibriumEvaporator::getE0(G4int /*a*/) const {
744 
745  if (verboseLevel > 3) {
746  G4cout << " >>> G4EquilibriumEvaporator::getE0" << G4endl;
747  }
748 
749  const G4double e0 = 200.0;
750 
751  return e0;
752 }