Geant4  10.01.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 
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 
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 
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
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 
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 
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 
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 
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 
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 
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)
static const G4double d1
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
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:196
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
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:110
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
static const G4double d2
G4double bindingEnergy(G4int A, G4int Z)
virtual void deExcite(const G4Fragment &target, G4CollisionOutput &output)
Definition: G4BigBanger.cc:68
void setTarget(const G4InuclParticle *target)
CLHEP::HepLorentzVector G4LorentzVector