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