Context Navigation

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

Last change on this file since 1315 was 1315, checked in by garnier, 14 years ago
update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 12.7 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.29 2010/05/21 17:56:34 mkelsey Exp $
26	// Geant4 tag: $Name: geant4-09-04-beta-cand-01 $
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
35	#define RUN
36
37	#include <cmath>
38	#include "G4NonEquilibriumEvaporator.hh"
39	#include "G4CollisionOutput.hh"
40	#include "G4InuclElementaryParticle.hh"
41	#include "G4InuclNuclei.hh"
42	#include "G4InuclSpecialFunctions.hh"
43	#include "G4LorentzConvertor.hh"
44	#include "G4NucleiProperties.hh"
45	#include "G4HadTmpUtil.hh"
46
47	using namespace G4InuclSpecialFunctions;
48
49
50	G4NonEquilibriumEvaporator::G4NonEquilibriumEvaporator()
51	: G4VCascadeCollider("G4NonEquilibriumEvaporator") {}
52
53
54	void G4NonEquilibriumEvaporator::collide(G4InuclParticle* /bullet/,
55	G4InuclParticle* target,
56	G4CollisionOutput& output) {
57
58	if (verboseLevel > 3) {
59	G4cout << " >>> G4NonEquilibriumEvaporator::collide" << G4endl;
60	}
61
62	// Sanity check
63	G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target);
64	if (!nuclei_target) {
65	G4cerr << " NonEquilibriumEvaporator -> target is not nuclei " << G4endl;
66	return;
67	}
68
69	const G4double a_cut = 5.0;
70	const G4double z_cut = 3.0;
71
72	#ifdef RUN
73	const G4double eexs_cut = 0.1;
74	#else
75	const G4double eexs_cut = 100000.0;
76	#endif
77
78	const G4double coul_coeff = 1.4;
79	const G4int itry_max = 1000;
80	const G4double small_ekin = 1.0e-6;
81	const G4double width_cut = 0.005;
82
83	G4double A = nuclei_target->getA();
84	G4double Z = nuclei_target->getZ();
85	G4LorentzVector PEX = nuclei_target->getMomentum();
86	G4LorentzVector pin = PEX;
87	G4double EEXS = nuclei_target->getExitationEnergy();
88	pin.setE(pin.e() + 0.001 * EEXS);
89	G4InuclNuclei dummy_nuc;
90	G4ExitonConfiguration config = nuclei_target->getExitonConfiguration();
91
92	G4double QPP = config.protonQuasiParticles;
93
94	G4double QNP = config.neutronQuasiParticles;
95
96	G4double QPH = config.protonHoles;
97
98	G4double QNH = config.neutronHoles;
99
100	G4double QP = QPP + QNP;
101	G4double QH = QPH + QNH;
102	G4double QEX = QP + QH;
103	G4InuclElementaryParticle dummy(small_ekin, 1);
104	G4LorentzConvertor toTheExitonSystemRestFrame;
105
106	toTheExitonSystemRestFrame.setBullet(dummy);
107
108	G4double EFN = FermiEnergy(A, Z, 0);
109	G4double EFP = FermiEnergy(A, Z, 1);
110
111	G4double AR = A - QP;
112	G4double ZR = Z - QPP;
113	G4int NEX = G4int(QEX + 0.5);
114	G4LorentzVector ppout;
115	G4bool try_again = (NEX > 0);
116
117	// Buffer for parameter sets
118	std::pair<G4double, G4double> parms;
119
120	while (try_again) {
121
122	if (A >= a_cut && Z >= z_cut && EEXS > eexs_cut) { // ok
123
124	if (verboseLevel > 2) {
125	G4cout << " A " << A << " Z " << Z << " EEXS " << EEXS << G4endl;
126	}
127
128	// update exiton system
129	G4double nuc_mass = dummy_nuc.getNucleiMass(A, Z);
130	PEX.setVectM(PEX.vect(), nuc_mass);
131	toTheExitonSystemRestFrame.setTarget(PEX, nuc_mass);
132	toTheExitonSystemRestFrame.toTheTargetRestFrame();
133	G4double MEL = getMatrixElement(A);
134	G4double E0 = getE0(A);
135	G4double PL = getParLev(A, Z);
136	G4double parlev = PL / A;
137	G4double EG = PL * EEXS;
138
139	if (QEX < std::sqrt(2.0 * EG)) { // ok
140
141	paraMakerTruncated(Z, parms);
142	const G4double& AK1 = parms.first;
143	const G4double& CPA1 = parms.second;
144
145	G4double VP = coul_coeff * Z * AK1 / (G4cbrt(A - 1.0) + 1.0) /
146	(1.0 + EEXS / E0);
147	// G4double DM1 = bindingEnergy(A, Z);
148	// G4double BN = DM1 - bindingEnergy(A - 1.0, Z);
149	// G4double BP = DM1 - bindingEnergy(A - 1.0, Z - 1.0);
150	G4double DM1 = G4NucleiProperties::GetBindingEnergy(G4lrint(A), G4lrint(Z));
151	G4double BN = DM1 - G4NucleiProperties::GetBindingEnergy(G4lrint(A-1.0), G4lrint(Z));
152	G4double BP = DM1 - G4NucleiProperties::GetBindingEnergy(G4lrint(A-1.0), G4lrint(Z-1.0));
153	G4double EMN = EEXS - BN;
154	G4double EMP = EEXS - BP - VP * A / (A - 1.0);
155	G4double ESP = 0.0;
156
157	if (EMN > eexs_cut) { // ok
158	G4int icase = 0;
159
160	if (NEX > 1) {
161	G4double APH = 0.25 * (QP * QP + QH * QH + QP - 3.0 * QH);
162	G4double APH1 = APH + 0.5 * (QP + QH);
163	ESP = EEXS / QEX;
164	G4double MELE = MEL / ESP / (AAA);
165
166	if (ESP > 15.0) {
167	MELE *= std::sqrt(15.0 / ESP);
168
169	} else if(ESP < 7.0) {
170	MELE *= std::sqrt(ESP / 7.0);
171
172	if (ESP < 2.0) MELE *= std::sqrt(ESP / 2.0);
173	};
174
175	G4double F1 = EG - APH;
176	G4double F2 = EG - APH1;
177
178	if (F1 > 0.0 && F2 > 0.0) {
179	G4double F = F2 / F1;
180	G4double M1 = 2.77 * MELE * PL;
181	std::vector<G4double> D(3, 0.0);
182	D[0] = M1 * F2 * F2 * std::pow(F, NEX - 1) / (QEX + 1.0);
183
184	if (D[0] > 0.0) {
185
186	if (NEX >= 2) {
187	D[1] = 0.0462 / parlev / G4cbrt(A) * QP * EEXS / QEX;
188
189	if (EMP > eexs_cut)
190	D[2] = D[1] * std::pow(EMP / EEXS, NEX) * (1.0 + CPA1);
191	D[1] = std::pow(EMN / EEXS, NEX) getAL(A);
192
193	if (QNP < 1.0) D[1] = 0.0;
194
195	if (QPP < 1.0) D[2] = 0.0;
196
197	try_again = NEX > 1 && (D[1] > width_cut * D[0] \|\|
198	D[2] > width_cut * D[0]);
199
200	if (try_again) {
201	G4double D5 = D[0] + D[1] + D[2];
202	G4double SL = D5 * inuclRndm();
203	G4double S1 = 0.;
204
205	for (G4int i = 0; i < 3; i++) {
206	S1 += D[i];
207
208	if (SL <= S1) {
209	icase = i;
210
211	break;
212	};
213	};
214	};
215	};
216
217	} else {
218	try_again = false;
219	};
220
221	} else {
222	try_again = false;
223	};
224	};
225
226	if (try_again) {
227
228	if (icase > 0) { // N -> N - 1 with particle escape
229	G4double V = 0.0;
230	G4int ptype = 0;
231	G4double B = 0.0;
232
233	if (A < 3.0) try_again = false;
234
235	if (try_again) {
236
237	if (icase == 1) { // neutron escape
238
239	if (QNP < 1.0) {
240	icase = 0;
241
242	} else {
243	B = BN;
244	V = 0.0;
245	ptype = 2;
246	};
247
248	} else { // proton esape
249	if (QPP < 1.0) {
250	icase = 0;
251
252	} else {
253	B = BP;
254	V = VP;
255	ptype = 1;
256
257	if (Z - 1.0 < 1.0) try_again = false;
258	};
259	};
260
261	if (try_again && icase != 0) {
262	G4double EB = EEXS - B;
263	G4double E = EB - V * A / (A - 1.0);
264
265	if (E < 0.0) {
266	icase = 0;
267
268	} else {
269	G4double E1 = EB - V;
270	G4double EEXS_new = -1.;
271	G4double EPART = 0.0;
272	G4int itry1 = 0;
273	G4bool bad = true;
274
275	while (itry1 < itry_max && icase > 0 && bad) {
276	itry1++;
277	G4int itry = 0;
278
279	while (EEXS_new < 0.0 && itry < itry_max) {
280	itry++;
281	G4double R = inuclRndm();
282	G4double X;
283
284	if (NEX == 2) {
285	X = 1.0 - std::sqrt(R);
286
287	} else {
288	G4double QEX2 = 1.0 / QEX;
289	G4double QEX1 = 1.0 / (QEX - 1.0);
290	X = std::pow(0.5 * R, QEX2);
291
292	for (G4int i = 0; i < 1000; i++) {
293	G4double DX = X * QEX1 *
294	(1.0 + QEX2 * X * (1.0 - R / std::pow(X, NEX)) / (1.0 - X));
295	X -= DX;
296
297	if (std::fabs(DX / X) < 0.01) break;
298
299	};
300	};
301	EPART = EB - X * E1;
302	EEXS_new = EB - EPART * A / (A - 1.0);
303	};
304
305	if (itry == itry_max \|\| EEXS_new < 0.0) {
306	icase = 0;
307
308	} else { // real escape
309	G4InuclElementaryParticle particle(ptype);
310
311	particle.setModel(5);
312	G4double mass = particle.getMass();
313	EPART *= 0.001; // to the GeV
314	// generate particle momentum
315	G4double pmod = std::sqrt(EPART * (2.0 * mass + EPART));
316	G4LorentzVector mom = generateWithRandomAngles(pmod,mass);
317	G4LorentzVector mom_at_rest;
318
319	G4double QPP_new = QPP;
320	G4double Z_new = Z;
321
322	if (ptype == 1) {
323	QPP_new -= 1.;
324	Z_new -= 1.0;
325	};
326
327	G4double QNP_new = QNP;
328
329	if (ptype == 2) QNP_new -= 1.0;
330
331	G4double A_new = A - 1.0;
332	G4double new_exiton_mass =
333	dummy_nuc.getNucleiMass(A_new, Z_new);
334	mom_at_rest.setVectM(-mom.vect(), new_exiton_mass);
335
336	G4LorentzVector part_mom =
337	toTheExitonSystemRestFrame.backToTheLab(mom);
338	part_mom.setVectM(part_mom.vect(), mass);
339
340	G4LorentzVector ex_mom =
341	toTheExitonSystemRestFrame.backToTheLab(mom_at_rest);
342	ex_mom.setVectM(ex_mom.vect(), new_exiton_mass);
343
344	// check energy conservation and set new exitation energy
345	EEXS_new = 1000.0 * (PEX.e() + 0.001 * EEXS -
346	part_mom.e() - ex_mom.e());
347
348	if (EEXS_new > 0.0) { // everything ok
349	particle.setMomentum(part_mom);
350	output.addOutgoingParticle(particle);
351	ppout += part_mom;
352
353	A = A_new;
354	Z = Z_new;
355	PEX = ex_mom;
356	EEXS = EEXS_new;
357	NEX -= 1;
358	QEX -= 1;
359	QP -= 1.0;
360	QPP = QPP_new;
361	QNP = QNP_new;
362	bad = false;
363	};
364	};
365	};
366
367	if (itry1 == itry_max) icase = 0;
368	};
369	};
370	};
371	};
372
373	if (icase == 0 && try_again) { // N -> N + 2
374	G4double TNN = 1.6 * EFN + ESP;
375	G4double TNP = 1.6 * EFP + ESP;
376	G4double XNUN = 1.0 / (1.6 + ESP / EFN);
377	G4double XNUP = 1.0 / (1.6 + ESP / EFP);
378	G4double SNN1 = csNN(TNP) * XNUP;
379	G4double SNN2 = csNN(TNN) * XNUN;
380	G4double SPN1 = csPN(TNP) * XNUP;
381	G4double SPN2 = csPN(TNN) * XNUN;
382	G4double PP = (QPP * SNN1 + QNP * SPN1) * ZR;
383	G4double PN = (QPP * SPN2 + QNP * SNN2) * (AR - ZR);
384	G4double PW = PP + PN;
385	NEX += 2;
386	QEX += 2.0;
387	QP += 1.0;
388	QH += 1.0;
389	AR -= 1.0;
390
391	if (AR > 1.0) {
392	G4double SL = PW * inuclRndm();
393
394	if (SL > PP) {
395	QNP += 1.0;
396	QNH += 1.0;
397
398	} else {
399	QPP += 1.0;
400	QPH += 1.0;
401	ZR -= 1.0;
402
403	if (ZR < 2.0) try_again = false;
404
405	};
406
407	} else {
408	try_again = false;
409	};
410	};
411	};
412
413	} else {
414	try_again = false;
415	};
416
417	} else {
418	try_again = false;
419	};
420
421	} else {
422	try_again = false;
423	};
424	};
425	// everything finished, set output nuclei
426	// the exitation energy has to be re-set properly for the energy
427	// conservation
428
429	G4LorentzVector pnuc = pin - ppout;
430	G4InuclNuclei nuclei(pnuc, A, Z);
431
432	nuclei.setModel(5);
433	nuclei.setEnergy();
434
435	pnuc = nuclei.getMomentum();
436	G4double eout = pnuc.e() + ppout.e();
437	G4double eex_real = 1000.0 * (pin.e() - eout);
438
439	nuclei.setExitationEnergy(eex_real);
440	output.addTargetFragment(nuclei);
441
442	return;
443	}
444
445	G4double G4NonEquilibriumEvaporator::getMatrixElement(G4double A) const {
446
447	if (verboseLevel > 3) {
448	G4cout << " >>> G4NonEquilibriumEvaporator::getMatrixElement" << G4endl;
449	}
450
451	G4double me;
452
453	if (A > 150.0) {
454	me = 100.0;
455
456	} else if(A > 20.0) {
457	me = 140.0;
458
459	} else {
460	me = 70.0;
461	};
462
463	return me;
464	}
465
466	G4double G4NonEquilibriumEvaporator::getE0(G4double ) const {
467
468	if (verboseLevel > 3) {
469	G4cout << " >>> G4NonEquilibriumEvaporator::getEO" << G4endl;
470	}
471
472	const G4double e0 = 200.0;
473
474	return e0;
475	}
476
477	G4double G4NonEquilibriumEvaporator::getParLev(G4double A,
478	G4double ) const {
479
480	if (verboseLevel > 3) {
481	G4cout << " >>> G4NonEquilibriumEvaporator::getParLev" << G4endl;
482	}
483
484	// const G4double par = 0.125;
485	G4double pl = 0.125 * A;
486
487	return pl;
488	}

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

Download in other formats: