Context Navigation

G4NonEquilibriumEvaporator.cc @ 1340

Last change on this file since 1340 was 1340, checked in by garnier, 14 years ago
update ti head
File size: 14.0 KB

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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: trunk/source/processes/hadronic/models/cascade/cascade/src/G4NonEquilibriumEvaporator.cc @ 1340

Download in other formats: