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