ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NonEquilibriumEvaporator.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4NonEquilibriumEvaporator.cc
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 //
27 // 20100114 M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
28 // 20100309 M. Kelsey -- Use new generateWithRandomAngles for theta,phi stuff;
29 // eliminate some unnecessary std::pow()
30 // 20100412 M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
31 // 20100413 M. Kelsey -- Pass buffers to paraMaker[Truncated]
32 // 20100517 M. Kelsey -- Inherit from common base class
33 // 20100617 M. Kelsey -- Remove "RUN" preprocessor flag and all "#else" code
34 // 20100622 M. Kelsey -- Use local "bindingEnergy()" function to call through.
35 // 20100701 M. Kelsey -- Don't need to add excitation to nuclear mass; compute
36 // new excitation energies properly (mass differences)
37 // 20100713 M. Kelsey -- Add conservation checking, diagnostic messages.
38 // 20100714 M. Kelsey -- Move conservation checking to base class
39 // 20100719 M. Kelsey -- Simplify EEXS calculations with particle evaporation.
40 // 20100724 M. Kelsey -- Replace std::vector<> D with trivial D[3] array.
41 // 20100914 M. Kelsey -- Migrate to integer A and Z: this involves replacing
42 // a number of G4double terms with G4int, with consequent casts.
43 // 20110214 M. Kelsey -- Follow G4InuclParticle::Model enumerator migration
44 // 20110922 M. Kelsey -- Follow G4InuclParticle::print(ostream&) migration
45 // 20120608 M. Kelsey -- Fix variable-name "shadowing" compiler warnings.
46 // 20121009 M. Kelsey -- Add some high-verbosity debugging output
47 // 20130622 Inherit from G4CascadeDeexciteBase, move to deExcite() interface
48 // with G4Fragment
49 // 20130808 M. Kelsey -- Use new object-version of paraMaker, for thread safety
50 // 20130924 M. Kelsey -- Replace std::pow with G4Pow::powN() for CPU speed
51 // 20150608 M. Kelsey -- Label all while loops as terminating.
52 
53 #include <cmath>
54 
55 #include "G4NonEquilibriumEvaporator.hh"
56 #include "G4SystemOfUnits.hh"
57 #include "G4CollisionOutput.hh"
58 #include "G4Fragment.hh"
59 #include "G4InuclElementaryParticle.hh"
60 #include "G4InuclNuclei.hh"
61 #include "G4InuclSpecialFunctions.hh"
62 #include "G4LorentzConvertor.hh"
63 #include "G4Pow.hh"
64 
65 using namespace G4InuclSpecialFunctions;
66 
67 
68 G4NonEquilibriumEvaporator::G4NonEquilibriumEvaporator()
69  : G4CascadeDeexciteBase("G4NonEquilibriumEvaporator"),
70  theParaMaker(verboseLevel), theG4Pow(G4Pow::GetInstance()) {}
71 
72 
73 void G4NonEquilibriumEvaporator::deExcite(const G4Fragment& target,
74  G4CollisionOutput& output) {
75  if (verboseLevel) {
76  G4cout << " >>> G4NonEquilibriumEvaporator::deExcite" << G4endl;
77  }
78 
79  if (verboseLevel>1) G4cout << " evaporating target:\n" << target << G4endl;
80 
81  const G4int a_cut = 5;
82  const G4int z_cut = 3;
83 
84  const G4double eexs_cut = 0.1;
85 
86  const G4double coul_coeff = 1.4;
87  const G4int itry_max = 1000;
88  const G4double small_ekin = 1.0e-6;
89  const G4double width_cut = 0.005;
90 
91  getTargetData(target);
92  G4LorentzVector pin = PEX; // Save original four-vector for later
93 
94  G4ExitonConfiguration config(target);
95  G4int QPP = config.protonQuasiParticles;
96  G4int QNP = config.neutronQuasiParticles;
97  G4int QPH = config.protonHoles;
98  G4int QNH = config.neutronHoles;
99 
100  G4int QP = QPP + QNP;
101  G4int QH = QPH + QNH;
102  G4int QEX = QP + QH;
103 
104  G4InuclElementaryParticle dummy(small_ekin, 1);
105  G4LorentzConvertor toTheExitonSystemRestFrame;
106  //*** toTheExitonSystemRestFrame.setVerbose(verboseLevel);
107  toTheExitonSystemRestFrame.setBullet(dummy);
108 
109  G4double EFN = FermiEnergy(A, Z, 0);
110  G4double EFP = FermiEnergy(A, Z, 1);
111 
112  G4int AR = A - QP;
113  G4int ZR = Z - QPP;
114  G4int NEX = QEX;
115  G4LorentzVector ppout;
116  G4bool try_again = (NEX > 0);
117 
118  // Buffer for parameter sets
119  std::pair<G4double, G4double> parms;
120 
121  while (try_again) { /* Loop checking 08.06.2015 MHK */
122  if (A >= a_cut && Z >= z_cut && EEXS > eexs_cut) { // ok
123  // update exiton system (include excitation energy!)
124  G4double nuc_mass = G4InuclNuclei::getNucleiMass(A, Z, EEXS);
125  PEX.setVectM(PEX.vect(), nuc_mass);
126  toTheExitonSystemRestFrame.setTarget(PEX);
127  toTheExitonSystemRestFrame.toTheTargetRestFrame();
128 
129  if (verboseLevel > 2) {
130  G4cout << " A " << A << " Z " << Z << " mass " << nuc_mass
131  << " EEXS " << EEXS << G4endl;
132  }
133 
134  G4double MEL = getMatrixElement(A);
135  G4double E0 = getE0(A);
136  G4double PL = getParLev(A, Z);
137  G4double parlev = PL / A;
138  G4double EG = PL * EEXS;
139 
140  if (QEX < std::sqrt(2.0 * EG)) { // ok
141  if (verboseLevel > 3)
142  G4cout << " QEX " << QEX << " < sqrt(2*EG) " << std::sqrt(2.*EG)
143  << " NEX " << NEX << G4endl;
144 
145  theParaMaker.getTruncated(Z, parms);
146  const G4double& AK1 = parms.first;
147  const G4double& CPA1 = parms.second;
148 
149  G4double VP = coul_coeff * Z * AK1 / (G4cbrt(A-1) + 1.0) /
150  (1.0 + EEXS / E0);
151  G4double DM1 = bindingEnergy(A,Z);
152  G4double BN = DM1 - bindingEnergy(A-1,Z);
153  G4double BP = DM1 - bindingEnergy(A-1,Z-1);
154  G4double EMN = EEXS - BN;
155  G4double EMP = EEXS - BP - VP * A / (A-1);
156  G4double ESP = 0.0;
157 
158  if (verboseLevel > 3) {
159  G4cout << " AK1 " << AK1 << " CPA1 " << " VP " << VP
160  << "\n bind(A,Z) " << DM1 << " dBind(N) " << BN
161  << " dBind(P) " << BP
162  << "\n EMN " << EMN << " EMP " << EMP << G4endl;
163  }
164 
165  if (EMN > eexs_cut) { // ok
166  G4int icase = 0;
167 
168  if (NEX > 1) {
169  G4double APH = 0.25 * (QP * QP + QH * QH + QP - 3 * QH);
170  G4double APH1 = APH + 0.5 * (QP + QH);
171  ESP = EEXS / QEX;
172  G4double MELE = MEL / ESP / (A*A*A);
173 
174  if (verboseLevel > 3)
175  G4cout << " APH " << APH << " APH1 " << APH1 << " ESP " << ESP
176  << G4endl;
177 
178  if (ESP > 15.0) {
179  MELE *= std::sqrt(15.0 / ESP);
180  } else if(ESP < 7.0) {
181  MELE *= std::sqrt(ESP / 7.0);
182  if (ESP < 2.0) MELE *= std::sqrt(ESP / 2.0);
183  };
184 
185  G4double F1 = EG - APH;
186  G4double F2 = EG - APH1;
187 
188  if (verboseLevel > 3)
189  G4cout << " MELE " << MELE << " F1 " << F1 << " F2 " << F2
190  << G4endl;
191 
192  if (F1 > 0.0 && F2 > 0.0) {
193  G4double F = F2 / F1;
194  G4double M1 = 2.77 * MELE * PL;
195  G4double D[3] = { 0., 0., 0. };
196  D[0] = M1 * F2 * F2 * theG4Pow->powN(F, NEX-1) / (QEX+1);
197  if (verboseLevel > 3) {
198  G4cout << " D[0] " << D[0] << " with F " << F
199  << " powN(F,NEX-1) " << theG4Pow->powN(F, NEX-1)
200  << G4endl;
201  }
202 
203  if (D[0] > 0.0) {
204 
205  if (NEX >= 2) {
206  D[1] = 0.0462 / parlev / G4cbrt(A) * QP * EEXS / QEX;
207 
208  if (EMP > eexs_cut)
209  D[2] = D[1] * theG4Pow->powN(EMP/EEXS, NEX) * (1.0 + CPA1);
210  D[1] *= theG4Pow->powN(EMN/EEXS, NEX) * getAL(A);
211 
212  if (verboseLevel > 3) {
213  G4cout << " D[1] " << D[1] << " with powN(EMN/EEXS, NEX) "
214  << theG4Pow->powN(EMN/EEXS, NEX) << G4endl
215  << " D[2] " << D[2] << " with powN(EMP/EEXS, NEX) "
216  << theG4Pow->powN(EMP/EEXS, NEX) << G4endl;
217  }
218 
219  if (QNP < 1) D[1] = 0.0;
220  if (QPP < 1) D[2] = 0.0;
221 
222  try_again = NEX > 1 && (D[1] > width_cut * D[0] ||
223  D[2] > width_cut * D[0]);
224 
225  if (try_again) {
226  G4double D5 = D[0] + D[1] + D[2];
227  G4double SL = D5 * inuclRndm();
228  G4double S1 = 0.;
229 
230  if (verboseLevel > 3)
231  G4cout << " D5 " << D5 << " SL " << SL << G4endl;
232 
233  for (G4int i = 0; i < 3; i++) {
234  S1 += D[i];
235  if (SL <= S1) {
236  icase = i;
237  break;
238  }
239  }
240 
241  if (verboseLevel > 3)
242  G4cout << " got icase " << icase << G4endl;
243  } // if (try_again)
244  } // if (NEX >= 2)
245  } else try_again = false; // if (D[0] > 0)
246  } else try_again = false; // if (F1>0 && F2>0)
247  } // if (NEX > 1)
248 
249  if (try_again) {
250  if (icase > 0) { // N -> N-1 with particle escape
251  if (verboseLevel > 3)
252  G4cout << " try_again icase " << icase << G4endl;
253 
254  G4double V = 0.0;
255  G4int ptype = 0;
256  G4double B = 0.0;
257 
258  if (A < 3.0) try_again = false;
259 
260  if (try_again) {
261 
262  if (icase == 1) { // neutron escape
263  if (verboseLevel > 3)
264  G4cout << " trying neutron escape" << G4endl;
265 
266  if (QNP < 1) icase = 0;
267  else {
268  B = BN;
269  V = 0.0;
270  ptype = 2;
271  };
272  } else { // proton esape
273  if (verboseLevel > 3)
274  G4cout << " trying proton escape" << G4endl;
275 
276  if (QPP < 1) icase = 0;
277  else {
278  B = BP;
279  V = VP;
280  ptype = 1;
281 
282  if (Z-1 < 1) try_again = false;
283  };
284  };
285 
286  if (try_again && icase != 0) {
287  if (verboseLevel > 3)
288  G4cout << " ptype " << ptype << " B " << B << " V " << V
289  << G4endl;
290 
291  G4double EB = EEXS - B;
292  G4double E = EB - V * A / (A-1);
293 
294  if (E < 0.0) icase = 0;
295  else {
296  G4double E1 = EB - V;
297  G4double EEXS_new = -1.;
298  G4double EPART = 0.0;
299  G4int itry1 = 0;
300  G4bool bad = true;
301 
302  /* Loop checking 08.06.2015 MHK */
303  while (itry1 < itry_max && icase > 0 && bad) {
304  itry1++;
305  G4int itry = 0;
306 
307  /* Loop checking 08.06.2015 MHK */
308  while (EEXS_new < 0.0 && itry < itry_max) {
309  itry++;
310  G4double R = inuclRndm();
311  G4double X;
312 
313  if (NEX == 2) {
314  X = 1.0 - std::sqrt(R);
315 
316  } else {
317  G4double QEX2 = 1.0 / QEX;
318  G4double QEX1 = 1.0 / (QEX-1);
319  X = theG4Pow->powA(0.5*R, QEX2);
320  if (verboseLevel > 3) {
321  G4cout << " R " << R << " QEX2 " << QEX2
322  << " powA(R, QEX2) " << X << G4endl;
323  }
324 
325  for (G4int i = 0; i < 1000; i++) {
326  G4double DX = X * QEX1 *
327  (1.0 + QEX2 * X * (1.0 - R / theG4Pow->powN(X, NEX)) / (1.0 - X));
328  if (verboseLevel > 3) {
329  G4cout << " NEX " << NEX << " powN(X, NEX) "
330  << theG4Pow->powN(X, NEX) << G4endl;
331  }
332 
333  X -= DX;
334 
335  if (std::fabs(DX / X) < 0.01) break;
336 
337  };
338  };
339  EPART = EB - X * E1;
340  EEXS_new = EB - EPART * A / (A-1);
341  } // while (EEXS_new < 0.0...
342 
343  if (itry == itry_max || EEXS_new < 0.0) {
344  icase = 0;
345  continue;
346  }
347 
348  if (verboseLevel > 2)
349  G4cout << " particle " << ptype << " escape " << G4endl;
350 
351  EPART /= GeV; // From MeV to GeV
352 
353  G4InuclElementaryParticle particle(ptype);
354  particle.setModel(G4InuclParticle::NonEquilib);
355 
356  // generate particle momentum
357  G4double mass = particle.getMass();
358  G4double pmod = std::sqrt(EPART * (2.0 * mass + EPART));
359  G4LorentzVector mom = generateWithRandomAngles(pmod,mass);
360 
361  // Push evaporated paricle into current rest frame
362  mom = toTheExitonSystemRestFrame.backToTheLab(mom);
363 
364  // Adjust quasiparticle and nucleon counts
365  G4int QPP_new = QPP;
366  G4int QNP_new = QNP;
367 
368  G4int A_new = A-1;
369  G4int Z_new = Z;
370 
371  if (ptype == 1) {
372  QPP_new--;
373  Z_new--;
374  };
375 
376  if (ptype == 2) QNP_new--;
377 
378  if (verboseLevel > 3) {
379  G4cout << " nucleus px " << PEX.px() << " py " << PEX.py()
380  << " pz " << PEX.pz() << " E " << PEX.e() << G4endl
381  << " evaporate px " << mom.px() << " py " << mom.py()
382  << " pz " << mom.pz() << " E " << mom.e() << G4endl;
383  }
384 
385  // New excitation energy depends on residual nuclear state
386  G4double mass_new = G4InuclNuclei::getNucleiMass(A_new, Z_new);
387 
388  EEXS_new = ((PEX-mom).m() - mass_new)*GeV;
389  if (EEXS_new < 0.) continue; // Sanity check for new nucleus
390 
391  if (verboseLevel > 3)
392  G4cout << " EEXS_new " << EEXS_new << G4endl;
393 
394  PEX -= mom;
395  EEXS = EEXS_new;
396 
397  A = A_new;
398  Z = Z_new;
399 
400  NEX--;
401  QEX--;
402  QP--;
403  QPP = QPP_new;
404  QNP = QNP_new;
405 
406  particle.setMomentum(mom);
407  output.addOutgoingParticle(particle);
408  ppout += mom;
409  if (verboseLevel > 3) {
410  G4cout << particle << G4endl
411  << " ppout px " << ppout.px() << " py " << ppout.py()
412  << " pz " << ppout.pz() << " E " << ppout.e() << G4endl;
413  }
414 
415  bad = false;
416  } // while (itry1<itry_max && icase>0
417 
418  if (itry1 == itry_max) icase = 0;
419  } // if (E < 0.) [else]
420  } // if (try_again && icase != 0)
421  } // if (try_again)
422  } // if (icase > 0)
423 
424  if (icase == 0 && try_again) { // N -> N + 2
425  if (verboseLevel > 3) G4cout << " adding excitons" << G4endl;
426 
427  G4double TNN = 1.6 * EFN + ESP;
428  G4double TNP = 1.6 * EFP + ESP;
429  G4double XNUN = 1.0 / (1.6 + ESP / EFN);
430  G4double XNUP = 1.0 / (1.6 + ESP / EFP);
431  G4double SNN1 = csNN(TNP) * XNUP;
432  G4double SNN2 = csNN(TNN) * XNUN;
433  G4double SPN1 = csPN(TNP) * XNUP;
434  G4double SPN2 = csPN(TNN) * XNUN;
435  G4double PP = (QPP * SNN1 + QNP * SPN1) * ZR;
436  G4double PN = (QPP * SPN2 + QNP * SNN2) * (AR - ZR);
437  G4double PW = PP + PN;
438  NEX += 2;
439  QEX += 2;
440  QP++;
441  QH++;
442  AR--;
443 
444  if (AR > 1) {
445  G4double SL = PW * inuclRndm();
446 
447  if (SL > PP) {
448  QNP++;
449  QNH++;
450  } else {
451  QPP++;
452  QPH++;
453  ZR--;
454  if (ZR < 2) try_again = false;
455  }
456  } else try_again = false;
457  } // if (icase==0 && try_again)
458  } // if (try_again)
459  } else try_again = false; // if (EMN > eexs_cut)
460  } else try_again = false; // if (QEX < sqrg(2*EG)
461  } else try_again = false; // if (A > a_cut ...
462  } // while (try_again)
463 
464  // everything finished, set output fragment
465 
466  if (output.numberOfOutgoingParticles() == 0) {
467  output.addRecoilFragment(target);
468  } else {
469  G4LorentzVector pnuc = pin - ppout;
470  output.addRecoilFragment(makeFragment(pnuc, A, Z, EEXS));
471 
472  if (verboseLevel>3)
473  G4cout << " remaining nucleus\n" << output.getRecoilFragment() << G4endl;
474  }
475 
476  validateOutput(target, output); // Check energy conservation, etc.
477  return;
478 }
479 
480 G4double G4NonEquilibriumEvaporator::getMatrixElement(G4int a) const {
481  if (verboseLevel > 3) {
482  G4cout << " >>> G4NonEquilibriumEvaporator::getMatrixElement" << G4endl;
483  }
484 
485  G4double me;
486 
487  if (a > 150) me = 100.0;
488  else if (a > 20) me = 140.0;
489  else me = 70.0;
490 
491  return me;
492 }
493 
494 G4double G4NonEquilibriumEvaporator::getE0(G4int ) const {
495  if (verboseLevel > 3) {
496  G4cout << " >>> G4NonEquilibriumEvaporator::getEO" << G4endl;
497  }
498 
499  const G4double e0 = 200.0;
500 
501  return e0;
502 }
503 
504 G4double G4NonEquilibriumEvaporator::getParLev(G4int a, G4int ) const {
505  if (verboseLevel > 3) {
506  G4cout << " >>> G4NonEquilibriumEvaporator::getParLev" << G4endl;
507  }
508 
509  // const G4double par = 0.125;
510  G4double pl = 0.125 * a;
511 
512  return pl;
513 }