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