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