Context Navigation

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

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

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

Download in other formats: