Geant4  10.02.p02
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 "G4BigBanger.hh"
80 #include "G4CascadeInterpolator.hh"
81 #include "G4CollisionOutput.hh"
82 #include "G4Fissioner.hh"
83 #include "G4Fragment.hh"
84 #include "G4InuclNuclei.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 
153  : G4CascadeDeexciteBase("G4EquilibriumEvaporator"),
154  theParaMaker(verboseLevel), QFinterp(XREP) {
155  parms.first.resize(6,0.);
156  parms.second.resize(6,0.);
157 }
158 
160 
163  theFissioner.setVerboseLevel(verbose);
164  theBigBanger.setVerboseLevel(verbose);
165 }
166 
167 
168 // Main operation
169 
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  // simple implementation of the equilibium evaporation a la Dostrowski
179  const G4double huge_num = 50.0;
180  const G4double small = -50.0;
181  const G4double prob_cut_off = 1.0e-15;
182  const G4double Q1[6] = { 0.0, 0.0, 2.23, 8.49, 7.72, 28.3 };
183  const G4int AN[6] = { 1, 1, 2, 3, 3, 4 };
184  const G4int Q[6] = { 0, 1, 1, 1, 2, 2 };
185  const G4double G[6] = { 2.0, 2.0, 6.0, 6.0, 6.0, 4.0 };
186  const G4double BE = 0.0063;
187  const G4double fisssion_cut = 1000.0;
188  const G4double cut_off_energy = 0.1;
189 
190  const G4double BF = 0.0242;
191  const G4int itry_max = 1000;
192  const G4int itry_global_max = 1000;
193  const G4double small_ekin = 1.0e-6;
194  const G4int itry_gam_max = 100;
195 
196  G4double W[8], u[6], V[6], TM[6];
197  G4int A1[6], Z1[6];
198 
199  getTargetData(target);
200  if (verboseLevel > 3) G4cout << " after noeq: eexs " << EEXS << G4endl;
201 
202  G4InuclElementaryParticle dummy(small_ekin, proton);
203  G4LorentzConvertor toTheNucleiSystemRestFrame;
204  //*** toTheNucleiSystemRestFrame.setVerbose(verboseLevel);
205  toTheNucleiSystemRestFrame.setBullet(dummy);
206 
207  G4LorentzVector ppout;
208 
209  // See if fragment should just be dispersed
210  if (explosion(A, Z, EEXS)) {
211  if (verboseLevel > 1) G4cout << " big bang in eql start " << G4endl;
212  theBigBanger.deExcite(target, output);
213 
214  validateOutput(target, output); // Check energy conservation
215  return;
216  }
217 
218  // If nucleus is in ground state, no evaporation
219  if (EEXS < cut_off_energy) {
220  if (verboseLevel > 1) G4cout << " no energy for evaporation" << G4endl;
221  output.addRecoilFragment(target);
222 
223  validateOutput(target, output); // Check energy conservation
224  return;
225  }
226 
227  // Initialize evaporation attempts
228  G4double coul_coeff = (A >= 100.0) ? 1.4 : 1.2;
229 
230  G4LorentzVector pin = PEX; // Save original target for testing
231 
232  G4bool try_again = true;
233  G4bool fission_open = true;
234  G4int itry_global = 0;
235 
236  /* Loop checking 08.06.2015 MHK */
237  while (try_again && itry_global < itry_global_max) {
238  itry_global++;
239 
240  // Set rest frame of current (recoiling) nucleus
241  toTheNucleiSystemRestFrame.setTarget(PEX);
242  toTheNucleiSystemRestFrame.toTheTargetRestFrame();
243 
244  if (verboseLevel > 2) {
246  G4cout << " A " << A << " Z " << Z << " mass " << nuc_mass
247  << " EEXS " << EEXS << G4endl;
248  }
249 
250  if (explosion(A, Z, EEXS)) { // big bang
251  if (verboseLevel > 2)
252  G4cout << " big bang in eql step " << itry_global << G4endl;
253 
255 
256  validateOutput(target, output); // Check energy conservation
257  return;
258  }
259 
260  if (EEXS < cut_off_energy) { // Evaporation not possible
261  if (verboseLevel > 2)
262  G4cout << " no energy for evaporation in eql step " << itry_global
263  << G4endl;
264 
265  try_again = false;
266  break;
267  }
268 
269  // Normal evaporation chain
270  G4double E0 = getE0(A);
271  G4double parlev = getPARLEVDEN(A, Z);
272  G4double u1 = parlev * A;
273 
275  const std::vector<G4double>& AK = parms.first;
276  const std::vector<G4double>& CPA = parms.second;
277 
278  G4double DM0 = bindingEnergy(A,Z);
279  G4int i(0);
280 
281  for (i = 0; i < 6; i++) {
282  A1[i] = A - AN[i];
283  Z1[i] = Z - Q[i];
284  u[i] = parlev * A1[i];
285  V[i] = 0.0;
286  TM[i] = -0.1;
287 
288  if (goodRemnant(A1[i], Z1[i])) {
289  G4double QB = DM0 - bindingEnergy(A1[i],Z1[i]) - Q1[i];
290  V[i] = coul_coeff * Z * Q[i] * AK[i] / (1.0 + EEXS / E0) /
291  (G4cbrt(A1[i]) + G4cbrt(AN[i]));
292  TM[i] = EEXS - QB - V[i] * A / A1[i];
293  if (verboseLevel > 3) {
294  G4cout << " i=" << i << " : QB " << QB << " u " << u[i]
295  << " V " << V[i] << " TM " << TM[i] << G4endl;
296  }
297  }
298  }
299 
300  G4double ue = 2.0 * std::sqrt(u1 * EEXS);
301  G4double prob_sum = 0.0;
302 
303  W[0] = 0.0;
304  if (TM[0] > cut_off_energy) {
305  G4double AL = getAL(A);
306  G4double A13 = G4cbrt(A1[0]);
307  W[0] = BE * A13*A13 * G[0] * AL;
308  G4double TM1 = 2.0 * std::sqrt(u[0] * TM[0]) - ue;
309 
310  if (TM1 > huge_num) TM1 = huge_num;
311  else if (TM1 < small) TM1 = small;
312 
313  W[0] *= G4Exp(TM1);
314  prob_sum += W[0];
315  }
316 
317  for (i = 1; i < 6; i++) {
318  W[i] = 0.0;
319  if (TM[i] > cut_off_energy) {
320  G4double A13 = G4cbrt(A1[i]);
321  W[i] = BE * A13*A13 * G[i] * (1.0 + CPA[i]);
322  G4double TM1 = 2.0 * std::sqrt(u[i] * TM[i]) - ue;
323 
324  if (TM1 > huge_num) TM1 = huge_num;
325  else if (TM1 < small) TM1 = small;
326 
327  W[i] *= G4Exp(TM1);
328  prob_sum += W[i];
329  }
330  }
331 
332  // fisson part
333  W[6] = 0.0;
334  if (A >= 100.0 && fission_open) {
335  G4double X2 = Z * Z / A;
336  G4double X1 = 1.0 - 2.0 * Z / A;
337  G4double X = 0.019316 * X2 / (1.0 - 1.79 * X1 * X1);
338  G4double EF = EEXS - getQF(X, X2, A, Z, EEXS);
339 
340  if (EF > 0.0) {
341  G4double AF = u1 * getAF(X, A, Z, EEXS);
342  G4double TM1 = 2.0 * std::sqrt(AF * EF) - ue;
343 
344  if (TM1 > huge_num) TM1 = huge_num;
345  else if (TM1 < small) TM1 = small;
346 
347  W[6] = BF * G4Exp(TM1);
348  if (W[6] > fisssion_cut*W[0]) W[6] = fisssion_cut*W[0];
349 
350  prob_sum += W[6];
351  }
352  }
353 
354  // again time to decide what next
355  if (verboseLevel > 2) {
356  G4cout << " Evaporation probabilities: sum = " << prob_sum
357  << "\n\t n " << W[0] << " p " << W[1] << " d " << W[2]
358  << " He3 " << W[3] << "\n\t t " << W[4] << " alpha " << W[5]
359  << " fission " << W[6] << G4endl;
360  }
361 
362  G4int icase = -1;
363 
364  if (prob_sum < prob_cut_off) { // photon emission chain
365  G4double UCR0 = 2.5 + 150.0 / A;
366  G4double T00 = 1.0 / (std::sqrt(u1 / UCR0) - 1.25 / UCR0);
367 
368  if (verboseLevel > 3)
369  G4cout << " UCR0 " << UCR0 << " T00 " << T00 << G4endl;
370 
371  G4int itry_gam = 0;
372  while (EEXS > cut_off_energy && try_again) {
373  itry_gam++;
374  G4int itry = 0;
375  G4double T04 = 4.0 * T00;
376  G4double FMAX;
377 
378  if (T04 < EEXS) {
379  FMAX = (T04*T04*T04*T04) * G4Exp((EEXS - T04) / T00);
380  } else {
381  FMAX = EEXS*EEXS*EEXS*EEXS;
382  };
383 
384  if (verboseLevel > 3)
385  G4cout << " T04 " << T04 << " FMAX (EEXS^4) " << FMAX << G4endl;
386 
387  G4double S(0), X1(0);
388  while (itry < itry_max) {
389  itry++;
390  S = EEXS * inuclRndm();
391  X1 = (S*S*S*S) * G4Exp((EEXS - S) / T00);
392 
393  if (X1 > FMAX * inuclRndm()) break;
394  };
395 
396  if (itry == itry_max) { // Maximum attempts exceeded
397  try_again = false;
398  break;
399  }
400 
401  if (verboseLevel > 2) G4cout << " photon escape ?" << G4endl;
402 
403  if (S < EEXS) { // Valid evaporate
404  S /= GeV; // Convert to Bertini units
405 
406  G4double pmod = S;
408 
409  // Push evaporated particle into current rest frame
410  mom = toTheNucleiSystemRestFrame.backToTheLab(mom);
411 
412  if (verboseLevel > 3) {
413  G4cout << " nucleus px " << PEX.px() << " py " << PEX.py()
414  << " pz " << PEX.pz() << " E " << PEX.e() << G4endl
415  << " evaporate px " << mom.px() << " py " << mom.py()
416  << " pz " << mom.pz() << " E " << mom.e() << G4endl;
417  }
418 
419  PEX -= mom; // Remaining four-momentum
420  EEXS -= S*GeV; // New excitation energy (in MeV)
421 
422  // NOTE: In-situ construction will be optimized away (no copying)
425 
426  if (verboseLevel > 3)
427  G4cout << output.getOutgoingParticles().back() << G4endl;
428 
429  ppout += mom;
430  } else {
431  if (itry_gam == itry_gam_max) try_again = false;
432  }
433  } // while (EEXS > cut_off
434  try_again = false;
435  } else { // if (prob_sum < prob_cut_off)
436  G4double SL = prob_sum * inuclRndm();
437  if (verboseLevel > 3) G4cout << " random SL " << SL << G4endl;
438 
439  G4double S1 = 0.0;
440  for (i = 0; i < 7; i++) { // Select evaporation scenario
441  S1 += W[i];
442  if (SL <= S1) {
443  icase = i;
444  break;
445  };
446  };
447 
448  if (icase < 0) continue; // Failed to choose scenario, try again
449 
450  if (icase < 6) { // particle or light nuclei escape
451  if (verboseLevel > 2) {
452  G4cout << " particle/light-ion escape ("
453  << (icase==0 ? "neutron" : icase==1 ? "proton" :
454  icase==2 ? "deuteron" : icase==3 ? "He3" :
455  icase==4 ? "triton" : icase==5 ? "alpha" : "ERROR!")
456  << ")?" << G4endl;
457  }
458 
459  G4double uc = 2.0 * std::sqrt(u[icase] * TM[icase]);
460  G4double ur = (uc > huge_num ? G4Exp(huge_num) : G4Exp(uc));
461  G4double d1 = 1.0 / ur;
462  G4double d2 = 1.0 / (ur - 1.0);
463  G4int itry1 = 0;
464  G4bool bad = true;
465 
466  while (itry1 < itry_max && bad) {
467  itry1++;
468  G4int itry = 0;
469  G4double EPR = -1.0;
470  G4double S = 0.0;
471 
472  while (itry < itry_max && EPR < 0.0) {
473  itry++;
474  G4double uu = uc + G4Log((1.0 - d1) * inuclRndm() + d2);
475  S = 0.5 * (uc * uc - uu * uu) / u[icase];
476  EPR = TM[icase] - S * A / (A - 1.0) + V[icase];
477  };
478 
479  if (verboseLevel > 3) {
480  G4cout << " EPR " << EPR << " V " << V[icase]
481  << " S " << S << " EEXS " << EEXS << G4endl;
482  }
483 
484  if (EPR > 0.0 && S > V[icase] && S < EEXS) { // real escape
485  if (verboseLevel > 2)
486  G4cout << " escape itry1 " << itry1 << " icase "
487  << icase << " S (MeV) " << S << G4endl;
488 
489  S /= GeV; // Convert to Bertini units
490 
491  if (icase < 2) { // particle escape
492  G4int ptype = 2 - icase;
493  if (verboseLevel > 2)
494  G4cout << " particle " << ptype << " escape" << G4endl;
495 
496  // generate particle momentum
497  G4double mass =
499 
500  G4double pmod = std::sqrt((2.0 * mass + S) * S);
501  G4LorentzVector mom = generateWithRandomAngles(pmod, mass);
502 
503  // Push evaporated particle into current rest frame
504  mom = toTheNucleiSystemRestFrame.backToTheLab(mom);
505 
506  if (verboseLevel > 2) {
507  G4cout << " nucleus px " << PEX.px() << " py " << PEX.py()
508  << " pz " << PEX.pz() << " E " << PEX.e() << G4endl
509  << " evaporate px " << mom.px() << " py " << mom.py()
510  << " pz " << mom.pz() << " E " << mom.e() << G4endl;
511  }
512 
513  // New excitation energy depends on residual nuclear state
514  G4double mass_new =
515  G4InuclNuclei::getNucleiMass(A1[icase],Z1[icase]);
516 
517  G4double EEXS_new = ((PEX-mom).m() - mass_new)*GeV;
518  if (EEXS_new < 0.0) continue; // Sanity check for new nucleus
519 
520  PEX -= mom; // Remaining four-momentum
521  EEXS = EEXS_new;
522 
523  A = A1[icase];
524  Z = Z1[icase];
525 
526  // NOTE: In-situ construction optimizes away (no copying)
528  ptype, G4InuclParticle::Equilib));
529 
530  if (verboseLevel > 3)
531  G4cout << output.getOutgoingParticles().back() << G4endl;
532 
533  ppout += mom;
534  bad = false;
535  } else { // if (icase < 2)
536  if (verboseLevel > 2) {
537  G4cout << " nucleus A " << AN[icase] << " Z " << Q[icase]
538  << " escape icase " << icase << G4endl;
539  }
540 
541  G4double mass =
542  G4InuclNuclei::getNucleiMass(AN[icase],Q[icase]);
543 
544  // generate particle momentum
545  G4double pmod = std::sqrt((2.0 * mass + S) * S);
546  G4LorentzVector mom = generateWithRandomAngles(pmod,mass);
547 
548  // Push evaporated particle into current rest frame
549  mom = toTheNucleiSystemRestFrame.backToTheLab(mom);
550 
551  if (verboseLevel > 2) {
552  G4cout << " nucleus px " << PEX.px() << " py " << PEX.py()
553  << " pz " << PEX.pz() << " E " << PEX.e() << G4endl
554  << " evaporate px " << mom.px() << " py " << mom.py()
555  << " pz " << mom.pz() << " E " << mom.e() << G4endl;
556  }
557 
558  // New excitation energy depends on residual nuclear state
559  G4double mass_new =
560  G4InuclNuclei::getNucleiMass(A1[icase],Z1[icase]);
561 
562  G4double EEXS_new = ((PEX-mom).m() - mass_new)*GeV;
563  if (EEXS_new < 0.0) continue; // Sanity check for new nucleus
564 
565  PEX -= mom; // Remaining four-momentum
566  EEXS = EEXS_new;
567 
568  A = A1[icase];
569  Z = Z1[icase];
570 
571  // NOTE: In-situ constructor optimizes away (no copying)
573  AN[icase], Q[icase], 0.*GeV,
575 
576  if (verboseLevel > 3)
577  G4cout << output.getOutgoingNuclei().back() << G4endl;
578 
579  ppout += mom;
580  bad = false;
581  } // if (icase < 2)
582  } // if (EPR > 0.0 ...
583  } // while (itry1 ...
584 
585  if (itry1 == itry_max || bad) try_again = false;
586  } else { // if (icase < 6)
587  if (verboseLevel > 2) {
588  G4cout << " fission: A " << A << " Z " << Z << " eexs " << EEXS
589  << " Wn " << W[0] << " Wf " << W[6] << G4endl;
590  }
591 
592  // Catch fission output separately for verification
595 
596  if (fission_output.numberOfFragments() == 2) { // fission ok
597  if (verboseLevel > 2) G4cout << " fission done in eql" << G4endl;
598 
599  // Move fission fragments to lab frame for processing
600  fission_output.boostToLabFrame(toTheNucleiSystemRestFrame);
601 
602  // Now evaporate the fission fragments individually
603  this->deExcite(fission_output.getRecoilFragment(0), output);
604  this->deExcite(fission_output.getRecoilFragment(1), output);
605 
606  validateOutput(target, output); // Check energy conservation
607  return;
608  } else { // fission forbidden now
609  fission_open = false;
610  }
611  } // End of fission case
612  } // if (prob_sum < prob_cut_off)
613  } // while (try_again
614 
615  // this time it's final nuclei
616 
617  if (itry_global == itry_global_max) {
618  if (verboseLevel > 3) {
619  G4cout << " ! itry_global " << itry_global_max << G4endl;
620  }
621  }
622 
623  G4LorentzVector pnuc = pin - ppout;
624 
625  // NOTE: In-situ constructor will be optimized away (no copying)
626  output.addOutgoingNucleus(G4InuclNuclei(pnuc, A, Z, 0.,
628 
629  if (verboseLevel > 3) {
630  G4cout << " remaining nucleus \n" << output.getOutgoingNuclei().back()
631  << G4endl;
632  }
633 
634 
635  validateOutput(target, output); // Check energy conservation
636  return;
637 }
638 
640  G4int z,
641  G4double e) const {
642  if (verboseLevel > 3) {
643  G4cout << " >>> G4EquilibriumEvaporator::explosion? ";
644  }
645 
646  const G4double be_cut = 3.0;
647 
648  // Different criteria from base class, since nucleus more "agitated"
649  G4bool bigb = (!(a >= 12 && z >= 0 && z < 3*(a-z)) &&
650  (e >= be_cut * bindingEnergy(a,z))
651  );
652 
653  if (verboseLevel > 3) G4cout << bigb << G4endl;
654 
655  return bigb;
656 }
657 
659  G4int z) const {
660  if (verboseLevel > 3) {
661  G4cout << " >>> G4EquilibriumEvaporator::goodRemnant(" << a << "," << z
662  << ")? " << (a>1 && z>0 && a>z) << G4endl;
663  }
664 
665  return a > 1 && z > 0 && a > z;
666 }
667 
669  G4double x2,
670  G4int a,
671  G4int /*z*/,
672  G4double ) const {
673  if (verboseLevel > 3) {
674  G4cout << " >>> G4EquilibriumEvaporator::getQF ";
675  }
676 
677  const G4double G0 = 20.4;
678  const G4double XMIN = 0.6761;
679  const G4double XMAX = 0.8274;
680 
681  G4double QFF = 0.0;
682 
683  if (x < XMIN || x > XMAX) {
684  G4double X1 = 1.0 - 0.02 * x2;
685  G4double FX = (0.73 + (3.33 * X1 - 0.66) * X1) * (X1*X1*X1);
686  G4double A13 = G4cbrt(a);
687  QFF = G0 * FX * A13*A13;
688  } else {
689  QFF = QFinterp.interpolate(x, QFREP);
690  }
691 
692  if (QFF < 0.0) QFF = 0.0;
693 
694  if (verboseLevel>3) G4cout << " returns " << QFF << G4endl;
695 
696  return QFF;
697 }
698 
700  G4int /*a*/,
701  G4int /*z*/,
702  G4double e) const {
703 
704  if (verboseLevel > 3) {
705  G4cout << " >>> G4EquilibriumEvaporator::getAF" << G4endl;
706  }
707 
708  // ugly parameterisation to fit the experimental fission cs for Hg - Bi nuclei
709  G4double AF = 1.285 * (1.0 - e / 1100.0);
710 
711  if(AF < 1.06) AF = 1.06;
712  // if(AF < 1.08) AF = 1.08;
713 
714  return AF;
715 }
716 
718  G4int /*z*/) const {
719 
720  if (verboseLevel > 3) {
721  G4cout << " >>> G4EquilibriumEvaporator::getPARLEVDEN" << G4endl;
722  }
723 
724  const G4double par = 0.125;
725 
726  return par;
727 }
728 
730 
731  if (verboseLevel > 3) {
732  G4cout << " >>> G4EquilibriumEvaporator::getE0" << G4endl;
733  }
734 
735  const G4double e0 = 200.0;
736 
737  return e0;
738 }
virtual void setVerboseLevel(G4int verbose=0)
static const G4double d1
#define A13
double S(double temp)
G4double z
Definition: TRTMaterials.hh:39
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)
G4double a
Definition: TRTMaterials.hh:39
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)
G4GLOB_DLL std::ostream G4cout
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)
static const double GeV
Definition: G4SIunits.hh:214
G4double G4Log(G4double x)
Definition: G4Log.hh:230
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
const G4double x[NPOINTSGL]
static const G4double * SL[nLA]
const std::vector< G4InuclElementaryParticle > & getOutgoingParticles() const
#define G4endl
Definition: G4ios.hh:61
G4InuclSpecialFunctions::paraMaker theParaMaker
static const double m
Definition: G4SIunits.hh:128
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:66
static const G4double d2
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